function [search_volume, voxel_volume]= ...
   glass_brain(input_file,input_thresh,mask_file,mask_thresh);

%GLASS_BRAIN
%
% GLASS_BRAIN( INPUT_FILE, INPUT_THRESH, [MASK_FILE, MASK_THRESH])
%
% Makes an SPM stlye maximum projection of INPUT_FILE above INPUT_THRESH
% and MASK_FILE 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. 
% Coordinates are in voxels, starting at 0, as in 'register'.
% If x or y step sizes are negative, then their voxel indices 
% are reversed so that they agree with 'register'. Since matlab can't
% reverse axis labels, they are negative - ignore the minus signs.
%
% SEARCH_VOLUME is the volume of the masked region in mm^3. 
%
% VOXEL_VOLUME is the volume of one voxel in mm^3.

if nargin<3
   mask_file=[]
   mask_thresh=[]
end

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

% glass brain:

if ~isempty(mask_file)
   d=fmris_read_image(mask_file,0,0);
   numslices=d.dim(3);
   d=fmris_read_image(mask_file,1:numslices,1);
   mask=d.data;
   mask= (mask>mask_thresh1 & mask<=mask_thresh2);
   voxel_volume=abs(prod(d.vox));
   search_volume=sum(sum(sum(mask)))*voxel_volume
   voxel_volume
   
   xf=find(max(max(mask,[],2),[],3));
   xr=min(xf):max(xf);
   yf=find(max(max(mask,[],1),[],3));
   yr=min(yf):max(yf);
   zf=find(max(max(mask,[],1),[],2));
   zr=min(zf):max(zf);
   mask=mask(xr,yr,zr);
   
   d=fmris_read_image(input_file);
   X=d.data(xr,yr,zr);
   X(X<=input_thresh)=input_thresh-1;
   X(~mask)=input_thresh-2;
else
   d=fmris_read_image(input_file);
   X=d.data;
   X(X<=input_thresh)=input_thresh-1;
   xr=1:d.dim(1);
   yr=1:d.dim(2);
   zr=1:d.dim(3);
   voxel_volume=abs(prod(d.vox));
   search_volume=prod(d.dim(1:3))*voxel_volume
   voxel_volume
end

ax=abs(d.vox(1))*length(xr);
ay=abs(d.vox(2))*length(yr);
az=abs(d.vox(3))*length(zr);

if d.vox(1)>0
   xr=xr-1;
   xlab='x';
else
   xr=-(d.dim(1)-xr);
   xlab='-x';
end
if d.vox(2)>0
   yr=yr-1;
   ylab='y';
else
   yr=-(d.dim(2)-yr);
   ylab='-y';
end

clf;
h=0.75;
axes('position',[.1 .15+ax/(ax+az)*h ay/(ay+ax)*h az/(ax+az)*h])
X1=squeeze(max(X,[],1));
imagesc(yr,zr-1,X1'); 
if d.vox(3)>0 
   axis xy; 
end 
ylabel('z');  grid on;
title([input_file ' > ' num2str(input_thresh)]);

axes('position',[.15+ay/(ay+ax)*h .15+ax/(ax+az)*h ax/(ay+ax)*h az/(ax+az)*h])
X2=squeeze(max(X,[],2));
imagesc(xr,zr-1,X2'); 
if d.vox(3)>0 
   axis xy; 
end 
xlabel(xlab); grid on;

axes('position',[.1 .1 ay/(ay+ax)*h*7/6 ax/(ax+az)*h])
X3=squeeze(max(X,[],3));
imagesc(yr,xr,X3); colorbar;
xlabel(ylab); ylabel(xlab); grid on;

return