function [SPM,VOL] = batch_spm_getSPM(spm_file, xCon, con_idx, p_corrected_flag, p_thresh, ext_thresh)
% computs a specified and thresholded SPM following parameter estimation
% batched for use with batch_spm_svc.m
%
% FORMAT [SPM,VOL] = spm_getSPM(spm_file, con_idx);
%
% spm_file - full path to SPM.mat file from which the SPM
%            variable will be extracted
% xCon     - actual xCon variable from xCon.mat file
% con_idx -  index in xCon for contrast to extract from -
%            xCon.mat is presumed to be in same directory as
%            spm_file
% p_corrected_flag - logical; 1 if thresh is corrected for
%            multiple comparison, 0 otherwise
% p_thresh - p-threshold to use
% ext_thresh - extent threshold to use (in voxels)
%
% SPM    - structure containing SPM, distribution & filtering details
% .swd   - SPM working directory - directory containing current SPM.mat
% .title - title for comparison (string)
% .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}]
% .Ic    - indices of contrasts (in xCon)
% .Im    - indices of masking contrasts (in xCon)
% .pm    - p-value for masking (uncorrected)
% .Ex    - flag for exclusive or inclusive masking
% .u     - height threshold
% .k     - extent threshold {voxels}
% .XYZ   - location of voxels {voxel coords}
% .XYZmm - location of voxels {mm}
%
%
% VOL    - structure containing details of volume analysed
% .S     - search Volume {voxels}
% .R     - search Volume {resels}
% .FWHM  - smoothness {voxels}     
% .M     - voxels - > mm matrix
% .iM    - mm -> voxels matrix
% .VOX   - voxel dimensions {mm} - column vector
% .DIM   - image dimensions {voxels} - column vector
% 

%__________________________________________________________________
%
% spm_getSPM prompts for an SPM and applies thresholds {u & k}
% to a point list of voxel values (specified with their locations {XYZ})
% This allows the SPM be displayed and characterized in terms of regionally 
% significant effects by subsequent routines.
% 
% This version is batched and written specifically for use
% with batch_spm_svc.m.  See the original spm_getSPM.m for
% more details.
%_______________________________________________________________________
% batch_spm_getSPM.m 1.0 Jeff Cooper July 2003
% based on
% @(#)spm_getSPM.m	2.37 Andrew Holmes, Karl Friston & Jean-Baptiste Poline 01/06/23



swd = fileparts(spm_file);

%-Load SPM.mat
%-----------------------------------------------------------------------
xSDM   = load(spm_file);	%-Contents of SPM.mat (in a struct)



%-Map beta and ResMS images, canonicalising relative pathnames
% (SPM result images stored with relative pathnames, for robustness.)
%-----------------------------------------------------------------------
xSDM.Vbeta  = ...
	spm_vol([repmat([swd,filesep],length(xSDM.Vbeta),1),char(xSDM.Vbeta)]);
xSDM.VResMS = spm_vol(fullfile(swd,xSDM.VResMS));


%-Get Stats data from SPM.mat
%-----------------------------------------------------------------------
xX     = xSDM.xX;			%-Design definition structure
XYZ    = xSDM.XYZ;			%-XYZ coordinates
S      = xSDM.S;			%-search Volume {voxels}
R      = xSDM.R;			%-search Volume {resels}
xSDM   = rmfield(xSDM,{'xX','XYZ'});	%-Remove (large) duplicate fields

%-Get/Compute mm<->voxel matrices & image dimensions from SPM.mat
%-----------------------------------------------------------------------
 M     = xSDM.M;
iM     = inv(M);
DIM    = xSDM.DIM;

%=======================================================================
% - C O N T R A S T S ,  S P M   C O M P U T A T I O N 
%=======================================================================


%-Compute & store contrast parameters, contrast/ESS images, & SPM images
%=======================================================================

% i is index of contrast


	%-Already got contrast/ESS image - remap it w/ full pathname
        %---------------------------------------------------------------
	xCon(con_idx).Vcon = spm_vol(fullfile(swd,xCon(con_idx).Vcon));


	%-Already got statistic image - remap it w/ full pathname
        %---------------------------------------------------------------
	xCon(con_idx).Vspm = spm_vol(fullfile(swd,xCon(con_idx).Vspm));



%-Compute (unfiltered) SPM pointlist for requested contrast
%=======================================================================



%-Compute conjunction SPM (as minimum SPM) within mask...
%-----------------------------------------------------------------------
    Z     = Inf;
    Z = min(Z,spm_sample_vol(xCon(con_idx).Vspm,XYZ(1,:),XYZ(2,:),XYZ(3,:),0));
    


%-Degrees of Fredom and STAT string describing marginal distribution
%-----------------------------------------------------------------------
edf   = [xCon(con_idx).eidf xX.erdf];

	str = '';


switch xCon(con_idx).STAT
case 'T'
	STATstr = sprintf('%c%s_{%.4g}','T',str,edf(2));
case 'F'
	STATstr = sprintf('%c%s_{[%.4g,%.4g]}','F',str,edf(1),edf(2));
end


%-Various parameters...
%-----------------------------------------------------------------------
STAT = xCon(con_idx).STAT;
n    = 1;
u    = -Inf;
k    = 0;


%=======================================================================
% - H E I G H T   &   E X T E N T   T H R E S H O L D S
%=======================================================================

%-Height threshold
%-----------------------------------------------------------------------
if ~isempty(XYZ)

    %-Get height threshold
    %-------------------------------------------------------------------
    if p_corrected_flag
	    u  = spm_uc(p_thresh,edf,STAT,xSDM.R,n,xSDM.S);
    else
        u = spm_u(p_thresh,edf,STAT);
    end

    %-Calculate height threshold filtering
    %-------------------------------------------------------------------
    Q     = find(Z > u);

    %-Apply height threshold
    %-------------------------------------------------------------------
    Z     = Z(:,Q);
    XYZ   = XYZ(:,Q);

    if isempty(Q)
        warning(sprintf('No voxels survive height threshold u=%0.2g',u));
    end

end % (if ~isempty(XYZ))


%-Extent threshold (only for allowed cases)
%-----------------------------------------------------------------------
if ~isempty(XYZ) & STAT == 'T'


    %-Calculate extent threshold filtering
    %-------------------------------------------------------------------
    A     = spm_clusters(XYZ);
    Q     = [];
    for i = 1:max(A)
        j = find(A == i);
        if length(j) >= ext_thresh; Q = [Q j]; end
    end

    % ...eliminate voxels
    %-------------------------------------------------------------------
    Z     = Z(:,Q);
    XYZ   = XYZ(:,Q);

    if isempty(Q)
	    warning(sprintf('No voxels survive extent threshold k=%0.2g',k));
    end

else

    ext_thresh = 0;

end % (if ~isempty(XYZ))


%=======================================================================
% - E N D
%=======================================================================

%-Assemble output structures of unfiltered data
%=======================================================================
SPM    = struct('swd',		swd,...
		'Z',		Z,...
		'n',		n,...
		'STAT',		STAT,...
		'df',		edf,...
		'STATstr',	STATstr,...
		'Ic',		con_idx,...
		'Im',		[],...
		'pm',		[],...
		'Ex',		[],...
		'u',		u,...
		'k',		ext_thresh,...
		'XYZ',		XYZ,...
		'XYZmm',	M(1:3,:)*[XYZ; ones(1,size(XYZ,2))]);

VOL    = struct('S',		S,...
		'R',		R,...
		'FWHM',		xSDM.FWHM,...
		'M',		M,...
		'iM',		iM,...
		'VOX',		sqrt(sum(M(1:3,1:3).^2))',...
		'DIM',		DIM);

% RESELS per voxel (density) if it exists
%-----------------------------------------------------------------------
if isfield(xSDM,'VRVP'), VOL.VRVP = spm_vol(xSDM.VRVP); end

