function m=plotvol(input_file_x,input_file_y,xlims,ylims,...
   mask_file,mask_thresh,nx,ny)

%PLOTVOL
% M = PLOTVOL( INPUT_FILE_X, INPUT_FILE_Y, XLIMS, YLIMS, ...
%                  [MASK_FILE, MASK_THRESH [, NX NY]])
%
% Plots values of INPUT_FILE_Y in the range [YLIMS(1) YLIMS(2)] versus 
% values of INPUT_FILE_X in the range [XLIMS(1) XLIMS(2)] 
% where MASK_FILE is above MASK_THRESH (ignored if empty; empty by default).
% If MASK_THRESH is a vector [a b], a<=b, then mask is a<MASK_FILE<=b. 
% Output is a NX by NY (default 64 by 64) matrix M of counts.

if nargin<5
   mask_file=[]
   mask_thresh=[]
end
if nargin<7
   nx=64;
   ny=64;
end

mask_thresh1=mask_thresh(1);
if length(mask_thresh)>=2
   mask_thresh2=mask_thresh(2);
else
   mask_thresh2=Inf;
end

x=fmris_read_image(input_file_x,0,0);
h=zeros(1,(nx+2)*(ny+2));
mask=ones(x.dim(1),x.dim(2));
for slice=1:x.dim(3)
   x=fmris_read_image(input_file_x,slice,1);
   y=fmris_read_image(input_file_y,slice,1);
   if ~isempty(mask_file)
      d=fmris_read_image(mask_file,slice,1);
      mask= (d.data>mask_thresh1 & d.data<=mask_thresh2);
   end
   xx=round((x.data(mask)-xlims(1))/(xlims(2)-xlims(1))*nx+0.5);
   xx=max(min(xx,nx+1),0)+1;
   yy=round((y.data(mask)-ylims(1))/(ylims(2)-ylims(1))*ny+0.5);
   yy=max(min(yy,ny+1),0)+1;
   h=h+hist(xx+(yy-1)*(nx+2),1:((nx+2)*(ny+2)));
end
mm=reshape(h,nx+2,ny+2);
m=mm(2:(nx+1),2:(ny+1));
xl=((1:nx)-0.5)/nx*(xlims(2)-xlims(1))+xlims(1);
yl=((1:ny)-0.5)/ny*(ylims(2)-ylims(1))+ylims(1);

subplot(2,2,3)
imagesc(xl,yl,m'); axis xy; title('counts');
xlabel(input_file_x); ylabel(input_file_y);

subplot(2,2,4)
mt=repmat(max(m,[],1),nx,1);
mdx=m./(mt+(mt<=0)).*(mt>0);
imagesc(xl,yl,mdx'); axis xy; title('divided by max in x');

subplot(2,2,1)
mt=repmat(max(m,[],2),1,ny);
mdy=m./(mt+(mt<=0)).*(mt>0);
imagesc(xl,yl,mdy'); axis xy; title('divided by max in y');

subplot(2,2,2)
plot(xl,m*yl'./sum(m,2),xl*m./sum(m,1),yl); 
xlim([min(xl) max(xl)]) ; ylim([min(yl) max(yl)]);
title('means')

return

input_file_x1='c:/keith/results/ha_093923_wresid_fwhm.mnc';
input_file_x2='c:/keith/results/ha_100326_wresid_fwhm.mnc';
input_file_x3='c:/keith/results/ha_101410_wresid_fwhm.mnc';
input_file_x4='c:/keith/results/ha_102703_wresid_fwhm.mnc';

d=fmris_read_image(input_file_x1);
tot=d.data;
d=fmris_read_image(input_file_x2);
tot=tot+d.data;
d=fmris_read_image(input_file_x3);
tot=tot+d.data;
d=fmris_read_image(input_file_x4);
tot=tot+d.data;
d.data=tot/4;

d.file_name='c:/keith/results/ha_average_fwhm.mnc';
d.parent_file=input_file_x1;
d.dim(4)=0;
fmris_write_image(d);

d=fmris_read_image(input_file_y);
fr=d.data;
fa=tot/4;
r=fr./(fa+(fa<=0)).*(fa>0);
subplot(2,2,2)
imagesc(r(:,:,5)',[0 2]); colorbar; axis xy;


input_file_x='c:/keith/results/ha_average_fwhm.mnc';
input_file_y='c:/keith/results/ha_multi_hmw_wresid_fwhm.mnc';
mask_file='c:/keith/data/brian_ha_19971124_1_100326_mri_MC.mnc';
mask_thresh=500;
xlims=[4 14];
ylims=[4 20];
nx=64;
ny=64;
plotvol(input_file_x,input_file_y,xlims,ylims,mask_file,mask_thresh,nx,ny);   
colormap(spectral)   
   

X=[1 1 1 1]';
contrast=[1];
which_stats=[0 0 1];
input_files_effect=['c:/keith/results/ha_093923_hmw_mag_ef.mnc';
                    'c:/keith/results/ha_100326_hmw_mag_ef.mnc';
                    'c:/keith/results/ha_101410_hmw_mag_ef.mnc';
                    'c:/keith/results/ha_102703_hmw_mag_ef.mnc'];
input_files_sdeffect=['c:/keith/results/ha_093923_hmw_mag_sd.mnc';
                      'c:/keith/results/ha_100326_hmw_mag_sd.mnc';
                      'c:/keith/results/ha_101410_hmw_mag_sd.mnc';
                      'c:/keith/results/ha_102703_hmw_mag_sd.mnc'];
output_file_base='c:/keith/results/ha_multi_hmw_infsmooth'

input_files_df=112
input_files_fwhm=6
df=multistat(input_files_effect,input_files_sdeffect,input_files_df, ...
   input_files_fwhm,X,contrast,output_file_base,which_stats,Inf)

input_file_x='c:/keith/results/ha_multi_hmw_infsmooth_sd.mnc';
input_file_y='c:/keith/results/ha_multi_hmw_nosmooth_sd.mnc';
mask_file='c:/keith/data/brian_ha_19971124_1_100326_mri_MC.mnc';
mask_thresh=500;
xlims=[0 1.5];
ylims=[0 1.5];
nx=64;
ny=64;
plotvol(input_file_x,input_file_y,xlims,ylims,mask_file,mask_thresh,nx,ny);   
colormap(spectral)   
   
input_file_y='c:/keith/results/ha_average_fwhm.mnc';
input_file_y='c:/keith/results/ha_multi_hmw_wresid_fwhm.mnc';
input_file_x='c:/keith/results/ha_multi_hmw_infsmooth_sd.mnc';
mask_thresh=500;
xlims=[0 1.5];
ylims=[4 14];
nx=64;
ny=64;
plotvol(input_file_x,input_file_y,xlims,ylims,mask_file,mask_thresh,nx,ny);   
colormap(spectral)