function view_slices(input_file, mask_file, mask_thresh, slices, frames, zlim, d3)

%VIEW_SLICES displays the slices in a 3D volume.
%
% VIEW_SLICES( INPUT_FILE [, MASK_FILE [, MASK_THRESH [, SLICES 
%           [, FRAMES [, ZLIM]]]]])
%
% INPUT_FILE is the name of a single 4D image file with one or more  
% volumes, or a matrix of image file names, each with a single 3D volume.
%
% MASK_FILE is a mask file. If empty, it is ignored. Default is [].
%
% MASK_THRESH defines the search volume as the first frame of MASK_FILE 
% > MASK_THRESH. If MASK_THRESH is a vector [a b], a<=b, then mask 
% is a < MASK_FILE <= b. If empty (default), calls fmri_mask_thresh. 
%
% SLICES is a vector of desired slices, starting at 0, as in 'register',
% displayed at top left. Default is 0:(numslices-1), i.e. all the slices. 
%
% FRAMES is a vector of frames, displayed in the top right. Default is 1.
% If INPUT_FILE is a matrix, then FRAMES refers to the rows.
%
% If SLICES and FRAMES is a scalar, then axis values are given in voxels,
% starting at 0, as in 'register'. Since matlab can't reverse axis labels, 
% they are negative if the step sizes are negative - ignore the minus signs.
%
% ZLIM is the limits [min max] for the image, ignored if empty (default).

%############################################################################
% COPYRIGHT:   Copyright 2002 K.J. Worsley, 
%              Department of Mathematics and Statistics,
%              McConnell Brain Imaging Center, 
%              Montreal Neurological Institute,
%              McGill University, Montreal, Quebec, Canada. 
%              worsley@math.mcgill.ca
%
%              Permission to use, copy, modify, and distribute this
%              software and its documentation for any purpose and without
%              fee is hereby granted, provided that this copyright
%              notice appears in all copies. The author and McGill University
%              make no representations about the suitability of this
%              software for any purpose.  It is provided "as is" without
%              express or implied warranty.
%############################################################################

% Defaults:

if nargin < 2
   mask_file=[]
end
if nargin < 3
   mask_thresh=[]
end
if nargin < 5
   frames=1
end
if nargin < 6
   zlim=[]
end
if nargin < 7
   d3=3
end

numfiles=size(input_file,1)
d=fmris_read_image(deblank(input_file(1,:)),0,0);
ext=input_file(min(findstr('.',input_file))+(1:3));
isminc= all(lower(ext)=='mnc');
isanal= all(lower(ext)=='img');

if nargin < 4
   slices=(1:d.dim(d3))-1
end

perm=[2 3; 1 3; 1 2];
d1=perm(d3,1);
d2=perm(d3,2);
if ~isempty(mask_file)
   if isempty(mask_thresh)
      mask_thresh=fmri_mask_thresh(mask_file);
   end
   mask_thresh1=mask_thresh(1);
   if length(mask_thresh)>=2
      mask_thresh2=mask_thresh(2);
   else
      mask_thresh2=Inf;
   end
   switch d3;
   case 1;
      m=fmris_read_image(mask_file,1:d.dim(3),1);
      m.data=squeeze(permute(m.data(slices+1,1:d.dim(2),1:d.dim(3)),[2 3 1]));
   case 2;
      m=fmris_read_image(mask_file,1:d.dim(3),1);
      m.data=squeeze(permute(m.data(1:d.dim(1),slices+1,1:d.dim(3)),[1 3 2]));
   case 3;
      m=fmris_read_image(mask_file,slices+1,1);
   end
   m.data= (m.data>mask_thresh1 & m.data<=mask_thresh2);
   xf=find(max(max(m.data,[],2),[],3));
   xr=min(xf):max(xf);
   yf=find(max(max(m.data,[],1),[],3));
   yr=min(yf):max(yf);
else
   xr=1:d.dim(d1);
   yr=1:d.dim(d2);
end
xp=xr-1; signx=' '; 
yp=yr-1; signy=' '; 
if isminc & d.vox(d1)<0; xp=-(d.dim(d1)-xp-1); signx='-'; end
if isminc & d.vox(d2)<0 
   if d2==3
      yp=fliplr(-yp);
   else
      yp=-(d.dim(d2)-yp-1); 
   end;
   signy='-'; 
end
if isanal & d1==1; xp=fliplr(-xp); signx='-'; end

lab='xyz';
xlab=[signx lab(d1) '=' num2str(min(xp)) ':' num2str(max(xp)) ' voxels'];
ylab=[signy lab(d2) '=' num2str(min(yp)) ':' num2str(max(yp)) ' voxels'];

mask_value=NaN;
nx=length(xr);
ny=length(yr);
nslices=length(slices)*length(frames);
Slices=repmat(slices,1,length(frames));
Frames=repmat(frames',1,length(slices))';
Frames=Frames(:)';
mx=max(round(sqrt(nx*abs(d.vox(d1))*ny*abs(d.vox(d2))*nslices)/(nx*abs(d.vox(d1)))),1);
my=ceil(nslices/mx);
M=zeros(mx*nx,my*ny).*NaN;
for i=1:mx
   for j=1:my
      islice=i+(my-j)*mx;
      if islice<=nslices
         slice=Slices(islice);
         frame=Frames(islice);
         if numfiles==1
            file=input_file;
            f=frame;
         else
            file=deblank(input_file(frame,:));
            f=1;
         end
         switch d3;
         case 1;
            d=fmris_read_image(file,1:d.dim(3),f);
            d.data=squeeze(permute(d.data(slices+1,1:d.dim(2),1:d.dim(3)),[2 3 1]));
         case 2;
            d=fmris_read_image(file,1:d.dim(3),f);
            d.data=squeeze(permute(d.data(1:d.dim(1),slices+1,1:d.dim(3)),[1 3 2]));
         case 3;
            d=fmris_read_image(file,slice+1,f);
         end
         v=d.data(xr,yr);
         if ~isempty(mask_file)
            v(m.data(xr,yr,find(slices==slice))==0)=mask_value;
         end
         if isanal & d1==1
            v=flipdim(v,1);
         end
         if isminc & d.vox(d2)<0 & d2==3
            v=flipdim(v,2);
         end
         M((i-1)*nx+(1:nx),(j-1)*ny+(1:ny))=v;
      end
   end
end
if length(frames)>1 | length(slices)>1 
   xp=(1:size(M,1))*abs(d.vox(d1));
   yp=(1:size(M,2))*abs(d.vox(d2));
end
sx=abs(d.vox(d1))*length(xp)*7;
sy=abs(d.vox(d2))*length(yp)*8;
pos=get(gca,'Position');
h=pos(4)/max(sx,sy);
set(gca,'Position',[pos(1)+(pos(3)-sx*h)/2 pos(2)+(pos(4)-sy*h)/2 sx*h sy*h]);
if isempty(zlim)
   M(isnan(M))=min(min(M));
   imagesc(xp,yp,M');
else
   M(isnan(M))=min(min(min(M)),min(zlim));
   imagesc(xp,yp,M',zlim);
end
axis xy;
if length(frames)>1 | length(slices)>1 
   set(gca,'XTick',[]); set(gca,'YTick',[]); 
end
tit=strrep(input_file,'_',' ');
if length(frames)==1
   tit=[tit ', frame ' num2str(frames)];  
end
if length(slices)==1
   tit=[tit ', slice ' num2str(slices)];  
end
title(tit);
for i=1:mx
   for j=1:my
      islice=i+(my-j)*mx;
      if islice<=nslices
         slice=Slices(islice);
         frame=Frames(islice);
         if length(slices)>1
            text((i-1+0.05)*nx*abs(d.vox(d1)),(j-1+0.9)*ny*abs(d.vox(d2)), ...
               num2str(slice),'Color','w');
         end
         if length(frames)>1
            text((i-1+0.9)*nx*abs(d.vox(d1)),(j-1+0.9)*ny*abs(d.vox(d2)), ...
               num2str(frame),'Color','w');
         end
      end
   end
end
xlabel(xlab); ylabel(ylab);
colorbar;
colormap(spectral);

return