function m_slice(SPM,VOL,hReg,varargin)

% Displays multiple (up to 24) transverse slices.
% based on spm_transverse.m from the SPM99 distribution
% ---------------------------
% Usage: Type "m_slice(SPM,VOL,hReg)" at the matlab prompt after
%        displaying the results through the "RESULTS" button.
%
%        Slices are peeled downward from the slice on which the
%        cursor is positioned (only the z dimension is relevant). 
%
%        The default max T value for the color bar is 7.
%        If a different T value is desired, it can be entered as
%        the fourth argument, e.g. "m_slice(SPM,VOL,hReg,9)"
%        for a max T value of 9.
% ---------------------------
% Kalina Christoff 2000-02-16
%
%
% for more details, please see the original spm_transverse.m file
%_______________________________________________________________________
% @(#)spm_transverse.m	2.15 Karl Friston - modified by John Ashburner 99/10/14

if nargin < 3,

        fprintf('\nPlease type "m_slice(SPM,VOL,hReg)" or\n');
        fprintf('"m_slice(SPM,VOL,hReg,<maxT>)", e.g. for a maxT value of 9, type\n');
        fprintf('"m_slice(SPM,VOL,hReg,9)"\n\n');
        error ('Insufficient number of input arguments');

end

%set orientation to default = 'axial'
%---------------------------------


% if max T is not specified, set default d=maxT=7
%------------------------------------------------
if nargin == 3 & isstruct(SPM),
        d=7; 
        axial=1;
        fprintf('\nThe maxT value on the color bar is set to a default of T=7.');
        fprintf('\nTo use a different value, type ');
        fprintf('"m_slice(SPM,VOL,hReg,<Tmax-value>)"\n');
        fprintf('e.g. "m_slice(SPM,VOL,hReg,9)"\n\n');
        init(SPM,VOL,hReg,d,axial);
end;

% if max T or orientation are specified, assign them to d=maxT and axial=[0 or 1]
%--------------------------------------------------------------------------------
if      nargin == 4 & ~ischar(varargin{1}) & isstruct(SPM),
        d=varargin{1};
        axial=1;
        init(SPM,VOL,hReg,d,axial);
elseif...
	nargin == 4 & ischar(varargin{1}) & isstruct(SPM),
        d=7;
        if strcmp(varargin{1},'coronal'); axial=0; else axial=1; end
        init(SPM,VOL,hReg,d,axial);
end

% if both max T and rientation are specified, assign them approriately
%---------------------------------------------------------------------
if      nargin == 5 & ~ischar(varargin{1}) & isstruct(SPM), 
        d=varargin{1};
        if strcmp(varargin{2},'coronal'); axial=0; else axial=1; end
        init(SPM,VOL,hReg,d,axial);
elseif... 
        nargin == 5 & ischar(varargin{1}) & isstruct(SPM),   
        d=varargin{2}; 
        if strcmp(varargin{1},'coronal'); axial=0; else axial=1; end
        init(SPM,VOL,hReg,d,axial);
end;


return;

%_______________________________________________________________________
%_______________________________________________________________________
function init(SPM,VOL,hReg,d,axial)

%-Get figure handles
%-----------------------------------------------------------------------
Fgraph = spm_figure('GetWin','Graphics');

%-Get the image on which to render
%-----------------------------------------------------------------------
%spms   = '/.automount/sperry/vol/disk56/kalina/analysis/sem_epi/GROUP_anat/group_anat.img';
spms   = spm_get(1,'.img','select an image for rendering');

spm('Pointer','Watch');

% Warn that it might take a while
%--------------------------------
fprintf('\n\nPlease wait, preparing images for rendering ...\n\n');

%-Delete previous axis and their pagination controls (if any)
%-----------------------------------------------------------------------
spm_results_ui('Clear',Fgraph);

% kch - delete all the glass brains and design matrix picture to make
% space for the slices
%---------------------------------------------------------------------
spm_figure('Clear',Fgraph)


% transv.blob is the functional image and transv.V is the anatomy
%-----------------------------------------------------------------
global transv
transv = struct('blob',[],'V',spm_vol(spms),'h',[],'hReg',hReg,'fig',Fgraph);
transv.blob = struct('xyz', round(SPM.XYZ), 't',SPM.Z, 'dim',VOL.DIM(1:3),...
	             'iM',VOL.iM,...
		     'vox', sqrt(sum(VOL.M(1:3,1:3).^2)), 'u', SPM.u);


%-Get current location and convert to pixel co-ordinates
%-----------------------------------------------------------------------
xyzmm  = spm_XYZreg('GetCoords',transv.hReg);
xyz    = round(transv.blob.iM(1:3,:)*[xyzmm; 1]);


% extract data from SPM [at one plane separation]
% and get background slices
%-----------------------------------------------------------------------------
dim    = ceil(transv.blob.dim(1:3)'.*transv.blob.vox); % the bounding box in mm
A      = transv.blob.iM*transv.V.mat;
hld    = 0;

% orientation case
% ---------------
switch axial

  case 1; % for axial orientation
  inpl=[1 2]; thrpl=3; inpl_crd='z';

  case 0; % for coronal orientation
  inpl=[1 3]; thrpl=2; inpl_crd='y';

end

if axial==1,
  zoomMfunc = inv(spm_matrix([0 0 -1  0 0 0  transv.blob.vox(1) transv.blob.vox(2) 1]));
  zoomManat =     spm_matrix([0 0  0  0 0 0  transv.blob.vox(1) transv.blob.vox(2) 1]);

elseif axial==0,
  zoomMfunc  = inv(spm_matrix([0 -1 0  0 0 0  transv.blob.vox(1) 1 transv.blob.vox(3)]));
  zoomManat  =     spm_matrix([0 0  0  pi/2 0 0  transv.blob.vox(1) 1 transv.blob.vox(3)]);
end

%-Determine paramters
%----------------------
slice_gap  = transv.blob.vox(thrpl);

%-Get bounding box borders
%------------------------
global sptl_BB

bb_min=sptl_BB(1,:); % gets the lower limit in the throuplane directions

% number of slices down or left of cursor
range      = ceil( (xyzmm(thrpl) - bb_min(thrpl)) / slice_gap )  +  1; 

% number of slices to be peeled (max of 24)
num_slices = min (24, range); 


maxD=0;

for i = 1 : num_slices

   % functional image
   % ----------------
   Q       = find(abs(transv.blob.xyz(thrpl,:) - (xyz(thrpl) -(i-1) ) ) < 0.5);
   T_tmp= full(sparse(transv.blob.xyz(inpl(1),Q),transv.blob.xyz(inpl(2),Q), ...
                      transv.blob.t(Q),transv.blob.dim(inpl(1)),transv.blob.dim(inpl(2))));
   T_tmp   = spm_slice_vol(T_tmp,zoomMfunc,dim(inpl),hld);
   T(:,:,i)= T_tmp;

   % anatomy image
   % -------------
   if axial==1,     mtrx =[1 0 0 0;0 1 0 0; 0 0 1 -(xyz(thrpl)-(i-1)); 0 0 0 1];
   elseif axial==0, mtrx =[1 0 0 0;0 1 0 -(xyz(thrpl)-(i-1)); 0 0 1 0; 0 0 0 1];
   end
   D            = zoomManat*mtrx*A; 
   D_tmp        = spm_slice_vol(transv.V,inv(D),dim(inpl),1);
   Dslice(:,:,i)= D_tmp;

   maxD         = max([maxD ; D_tmp(:)]);

end


for i= 1:num_slices
  Dslice(:,:,i)=Dslice(:,:,i)/maxD;
end


%-Configure {128 level} colormap
%-----------------------------------------------------------------------
cmap   = get(Fgraph,'Colormap');
if size(cmap,1) ~= 128
	figure(Fgraph)
	spm_figure('Colormap','gray-hot')
	cmap   = get(Fgraph,'Colormap');
end

D      = length(cmap)/2;


for i = 1:num_slices
   T_tmp          = T(:,:,i);
   D_tmp          = Dslice(:,:,i);

              Q        = T_tmp(:) > transv.blob.u; 
              T_tmp    = T_tmp(Q)/d;  
              D_tmp(Q) = 1 + T_tmp;
              T_tmp    = D*D_tmp; 

   Tslice(:,:,i) = T_tmp;
   Dslice(:,:,i) = D_tmp;
end

set(Fgraph,'Units','normalized')
siz    = get(Fgraph,'Position');
siz    = siz(3:4);


%-Render activation foci on background images
%-----------------------------------------------------------------------
	P = xyz.*transv.blob.vox';

% Matrix of positions, ([left, bottom, width, height] x 24)
%----------------------------------------------------------
RECT(1,:) = [repmat([.81 .61 .41 .21 .01],1,2) .61 .41 .21 .01 repmat([.81 .61 .41 .21 .01],1,2)];
RECT(2,:) = [repmat(-.02,1,5) repmat(.16,1,5) repmat(.34,1,4) repmat(.52,1,5) repmat(.7,1,5)];
RECT(3,:) = .1667; % 100
RECT(4,:) = .2; % 120

for i = 1:num_slices        % index for slices   
   
        j=2*i-1;             % index for transv.h

	transv.h(j) = axes('Units','normalized','Parent',Fgraph, ... 
                           'Position',RECT(:,i)');
	transv.h(j+1) = image(rot90(Tslice(:,:,i)));
	axis image; axis off;
        title( [sprintf(inpl_crd) , sprintf(' = %0.0fmm',(xyzmm(thrpl) - (i-1)*slice_gap))]  );
      
end

clear j

	% colorbar
	%-----------------------------------------------------------------------
	q      = [ .87 .37 .04 .15];
	transv.h(num_slices*2) = axes('Units','normalized','Parent',Fgraph,'Position',q,'Visible','off');
	transv.h(num_slices*2+1) = image([0 d/32],[0 d],[1:D]' + D);
	str    = [SPM.STAT ' value'];
	axis xy; title(str,'FontSize',9);
	set(gca,'XTickLabel',[]);

        % annotation
        %-----------------------------------------------------------------------
        q      = [ 0 .9 1 .18];
        transv.h(num_slices*2+2) = axes('Units','normalized','Parent',Fgraph,'Position',q,'Visible','off');

        AnBox  = get(transv.h(num_slices*2+2),'Position'); 
        set(transv.h(num_slices*2+2),'YLim',[0,AnBox(4)])

        transv.h(num_slices*2+3) = text(0.03,0.045,SPM.title,'Parent', ...
                                     transv.h(num_slices*2+2),'FontWeight','Bold','FontSize',16);

        transv.h(num_slices*2+4) = text(0.03,0.02, ['...' SPM.swd(length(SPM.swd)-40 : length(SPM.swd))], ...
                            'Parent',transv.h(num_slices*2+2),'FontWeight','Bold','FontSize',12);

        transv.h(num_slices*2+5) = text(0.6, 0.04, ...
                sprintf('Height threshold %c = %0.2f',SPM.STAT,SPM.u),'Parent', transv.h(num_slices*2+2), ...
                         'FontSize',12);

        transv.h(num_slices*2+6) = text(0.6 ,0.02, ...
                sprintf('Extent threshold k = %0.0f voxels',SPM.k),'Parent', transv.h(num_slices*2+2), ...
                         'FontSize',12);



%spm_print;




