
function TabDat = spm_VOI(SPM,VOL,Dis,Num,hReg)
% List of local maxima and adjusted p-values for a small Volume of Interest
% FORMAT TabDat = spm_VOI(SPM,VOL,Dis,Num,hReg)
%
% SPM    - structure containing SPM, distribution & filtering details
%        - required fields are:
% .swd   - SPM working directory - directory containing current SPM.mat
% .Z     - minimum of n Statistics {filtered on u and k}
% .n     - number of conjoint tests
% .STAT  - distribution {Z, T, X or F}
% .df    - degrees of freedom [df{interest}, df{residual}]
% .u     - height threshold
% .k     - extent threshold {resels}
% .XYZ   - location of voxels {voxel coords}
% .XYZmm - location of voxels {mm}
%
% VOL    - structure containing details of volume analysed
%        - required fields are:
% .S     - search Volume {voxels}
% .R     - search Volume {resels}
% .FWHM  - smoothness {voxels}
% .M     - voxels - > mm matrix
% .VOX   - voxel dimensions {mm}
%
% Dis    - Minimum distance between maxima
%          (Passed to spm_list.m: Defaults on missing or empty)
% Num    - Maxiumum number of local maxima tabulated per cluster
%          (Passed to spm_list.m: Defaults on missing or empty)
% hReg   - Handle of results section XYZ registry (see spm_results_ui.m)
%
% TabDat - Structure containing table data
%        - see spm_list for definition
%
%_______________________________________________________________________
%
% spm_VOI is  called by the SPM results section and takes variables in
% SPM and VOL to compute p-values corrected for a specified volume of
% interest.
%
% The volume of interest may be defined as a box or sphere centred on
% the current voxel, by the cluster in which it is embedded, or by an
% external mask image.
%
% If the VOI is defined by the cluster in which it is embedded, this
% cluster must have been defined independently of the SPM using a mask
% based on an orthogonal contrast and u = -Inf (i.e. p = 1)
%
% External mask images should be in the same orientation as the SPM
% (i.e. as the input used in stats estimation). The VOI is defined by
% voxels with values greater than 0.
%
% See also: spm_list
%_______________________________________________________________________
% @(#)spm_VOI.m	2.10 Karl Friston 99/11/17

%-Parse arguments
%-----------------------------------------------------------------------
if nargin<2,     error('insufficient arguments'), end
if nargin<5,	 hReg = []; end
if nargin<4,     Num  = []; end
if isempty(Num), Num  = 16; end
if nargin<3,     Dis  = []; end
if isempty(Dis), Dis  = 04; end

%-Title
%-----------------------------------------------------------------------
spm('FigName',['SPM{',SPM.STAT,'}: Small Volume Correction']);

%-Get current location
%-----------------------------------------------------------------------
xyzmm   = spm_results_ui('GetCoords');

  


%-Specify search volume
%-----------------------------------------------------------------------
str     = sprintf(' at [%.0f,%.0f,%.0f]',xyzmm(1),xyzmm(2),xyzmm(3));
SPACE   = spm_input('Search volume...',-1,'m',...
		{['Sphere',str],['Box',str],'Nearest cluster',...
		'Image'},['S','B','V','I']);
Q       = ones(1,size(SPM.XYZmm,2));
vsc     = [1 1 1];				%-Voxel size scaling for FWHM

switch SPACE, case 'S'                                          % Sphere
	%---------------------------------------------------------------
	D     = spm_input('radius of spherical VOI {mm}',-2);
	str   = sprintf('%0.1fmm sphere',D);
	j     = find(sum((SPM.XYZmm - xyzmm*Q).^2) <= D^2);
	D     = D./VOL.VOX;
	S     = (4/3)*pi*prod(D);

case 'B'                                                           % Box
	%---------------------------------------------------------------
	D     = spm_input('box dimensions [k l m] {mm}',-2);
	str   = sprintf('%0.1f x %0.1f x %0.1f mm box',D(1),D(2),D(3));
	j     = find(all(abs(SPM.XYZmm - xyzmm*Q) <= D(:)*Q/2));
	D     = D(:)./VOL.VOX;
	S     = prod(D);

case 'V'                                                        %-Voxels
	%---------------------------------------------------------------
	if ~length(SPM.XYZ)
		spm('alert!','No suprathreshold clusters!',mfilename,0);
		spm('FigName',['SPM{',SPM.STAT,'}: Results']);
		return
	end

	[xyzmm,i] = spm_XYZreg('NearestXYZ',xyzmm,SPM.XYZmm);
	spm_results_ui('SetCoords',xyzmm);
	A     = spm_clusters(SPM.XYZ);
	j     = find(A == A(i));
	str   = sprintf('%0.0f voxel cluster',length(j));
	D     = SPM.XYZ(:,j);
	S     = length(j);

case 'I'                                                         % Image
	%---------------------------------------------------------------
	%-The VOI image must be in the same (real) space as the SPM:
	%-Although the masking code below does take account of the
	% .mat file for the VOI image, the orientation of the VOI image
	% .mat file is *not* used in the resel calculation,
	% (see spm_resels_vol.c) so any mat file should contain
	% no rotations / shears relative to the SPM, for accuracy.
	im    = spm_get(1,'img','Image defining volume subset');
	D     = spm_vol(im);
   tM    = D.mat \ VOL.M;
   %if any(tM([2:5,7:10,12]))  %bug fix - sue
   if any(tM([2:5,7:10,12])>eps)

		spm('alert!',{	'Mask image rotated/sheared!',...
				'(relative to SPM image)',...
				'Can''t use for SVC.'},mfilename,0);
		spm('FigName',['SPM{',SPM.STAT,'}: Results']);
		return
	end
	mXYZ  = D.mat \ [SPM.XYZmm; ones(1, size(SPM.XYZmm, 2))];
	j     = spm_sample_vol(D, mXYZ(1,:),mXYZ(2,:),mXYZ(3,:),0) > 0;	
	str   = sprintf('image mask: %s',spm_str_manip(im,'a30'));
	%-Compute in-mask volume S:
	% Correct for differences in mask and SPM voxel sizes
	Y     = spm_read_vols(D);
	vsc   = sqrt(sum(D.mat(1:3,1:3).^2)) ./ VOL.VOX';
	S     = sum(Y(:)>0) * prod(vsc);

end

spm('Pointer','Watch')

%-Select voxels within subspace
%-----------------------------------------------------------------------
SPM.Z     = SPM.Z(j);
SPM.XYZ   = SPM.XYZ(:,j);
SPM.XYZmm = SPM.XYZmm(:,j);
VOL.R     = spm_resels(VOL.FWHM./vsc,D,SPACE);
VOL.S     = S;


%-Tabulate p values
%-----------------------------------------------------------------------
str = sprintf('search volume: %s',str);
if any(strcmp(SPACE,{'S','B','V'}))
	str = sprintf('%s at [%.0f,%.0f,%.0f]',str,xyzmm(1),xyzmm(2),xyzmm(3));
end
TabDat = spm_list('List',SPM,VOL,Dis,Num,str,hReg);

%-Reset title
%-----------------------------------------------------------------------
spm('FigName',['SPM{',SPM.STAT,'}: Results']);


%%%%%%%%%%%%%%%%%%modified to update MIP, save coord and write filtered image  _Sue...%%%%%%%%%%%%%%%%

figure;
colormap(gray);
spm_mip(SPM.Z, SPM.XYZmm, VOL.M, VOL.DIM);

 
XYZmm = SPM.XYZmm';

% - this next bit added by Jeff 10/22/02 for Gayle to sort out .tal files by cluster and
% significance level, a la glassbrain3.
% **************************************
% **************************************

%Assemble voxel data structure--------------------

J(:,1) = SPM.Z';
J(:,2:4) = SPM.XYZmm';  
%Note that indices for SPM.Z and SPM.XYZ reference same voxel - thus they do 
%here as well
J(:,5) = (spm_clusters(SPM.XYZ))';

%J is now a double array of n x 5 size, where n is the number of voxels in
%the SPM.  Its first column is the voxel T-score (or whatever statistic is 
%listed in SPM.STAT); its 5th is the voxel cluster
%number.  The columns in between are the coordinates in mm in MNI space for
%each voxel.

J = sortrows(J,5);
%J is now sorted into contiguous cluster rows (which it should have been anyways,
%but it's nice to guarantee it); we can extract each cluster
%and treat it as a sub-array of voxels.


p = max(J(:,5));
%p is the number of clusters.

for i = 1:p
    k = find(J(:,5)==i);
    E(i) = {J(k,:)};
    %Each E(i) now is a double array consisting of all and only the voxels in 
    %cluster i.  This will do us nicely.
    Q = sortrows(E{i},1);
    Q = flipud(Q);
    E{i} = Q;
    r(i) = Q(1,1);
end

%E is now a cell array in which each cell is a sorted list of voxels in a given
%cluster, in descending order of Z-score.  That's nice.  We now need to get
%E in order of cluster rank, which is easy enough.  SPM ranks clusters in order
%of lowest single uncorrected probability in that cluster.  Conveniently, r now 
%indexes those peak Z-scores.

r(2,:) = double(1:p);
%tack on the cluster number
r = r';
r = sortrows(r,1);
r = flipud(r);

XYZmm_tosave = [];

% This XYZ plus T business added by Jeff 2/11/03 for Ali
% and Kevin - the sprlio study.

XYZmm_tosave_plusT = [];

for i = 1:p
    curr = E{r(i,2)};
    XYZmm_tosave = cat(1,XYZmm_tosave,curr(:,2:4));
    XYZmm_tosave_plusT = cat(1, XYZmm_tosave_plusT, curr(:,1:4));
end

t_stats = XYZmm_tosave_plusT(:,1);
XYZmm_tosave_plusT = XYZmm_tosave_plusT(:,2:4);
XYZmm_tosave_plusT = horzcat(XYZmm_tosave_plusT, t_stats);


% **************************************
% **************************************
% - end added by Jeff 10/22/02

%tal     = spm_str_manip(spm_input('output name for .tal',1,'s'),'sdv');
%keyboard;
saveopt = spm_input_ui('How would you like to save S.V.C image?...',...
            1,'m','TAL and ANALYZE|TAL only|ANALYZE only| TAL plus T | No save',...
            [],1);

if ((saveopt == 1) | (saveopt == 2))
    resdir = spm_get(-1,'.','Select directory to save TAL file.');
    cd(resdir)
    Q = spm_str_manip(spm_input('Choose TAL filename',1,'s'),'sdv');
    save(cat(2,Q,'.tal'), 'XYZmm_tosave', '-ASCII');
    fprintf('\n\n----------------------------------------------');
    fprintf('\nA new file %s.tal was created in',Q);
    fprintf('\n%s\n',pwd);
    fprintf('\nIt contans a list of the coordinates (in mm) for the voxels surviving SVC\n');
end
if ((saveopt == 1) | (saveopt == 3))
    resdir = spm_get(-1,'.','Select directory to save ANALYZE image.');
    cd(resdir)
    Vo = spm_write_filtered(SPM.Z,SPM.XYZ,VOL.DIM,VOL.M,'SPM-filtered');
end
if (saveopt==4)
    resdir = spm_get(-1,'.','Select directory to save TAL + T file.');
    cd(resdir)
    Q = spm_str_manip(spm_input('Choose filename',1,'s'),'sdv');
    save(cat(2,Q,'.tal'), 'XYZmm_tosave_plusT', '-ASCII');
    fprintf('\n\n----------------------------------------------');
    fprintf('\nA new file %s.tal was created in',Q);
    fprintf('\n%s\n',pwd);
    fprintf('\nIt contans a list of the coordinates (in mm) for the voxels surviving SVC\n');
    fprintf('as well as a fourth column, containing the SVC-ed T-statistic for that row''s voxel.\n');
end
