
function batch_spm_svc(varargin)
% Extract p-values and coords for SPM filtered by region,
% and optionally save out active voxel counts and filtered
% image or coordinate files.  Based off spm_VOI in SPM99.

% FORMAT batch_spm_svc(expdir, names, resdir, con_list, roi_list,...
%           p_corrected_flag, p_thresh, ext_thresh, save_rois_flag, saveopt, save_voxcount_flag, outdir)
% FORMAT batch_spm_svc
%
% All the arguments are optional - any that are not supplied
% will be asked for - but they must come strictly in the
% order they're specified here, i.e. the 4th arg, if
% supplied, must be the con_list.  You can't, for example,
% supply an ROI_list and a threshold without also specifying
% a con_list and a p_corrected_flag.
%
% expdir -  full-path location of experiment's home
%           directory, in which all subject directories are
%           located.
% names  -  string matrix (cell or character), each line
%           containing the name of a single subject
%           directory, relative to expdir.
% resdir -  string containing path to results directory for
%           each subject dir, relative to the subject
%           directory.  So if data structure is
%           ./joe_subject/results/results_connect, and the
%           desired directory is results_connect, resdir
%           should be 'results/results_connect'.
% con_list - string matrix (cell or char), each line
%           containing the name of a single contrast file
%           from which activated voxels are to be taken.
%           This is a filename relative to resdir, not a
%           fullpath name.
% roi_list - string matrix (cell or char), each line
%           containing the full-path filename for an ROI
%           to serve as an image mask for a small volume
%           correction.
%
% 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}
%
%_______________________________________________________________________
%
% batch_spm_svc is based off the SPM99 function, spm_VOI,
% which computed p-values corrected for a specified volume
% of interest.  In practice, it was often used to create
% functional ROIs consisting of all active voxels within a
% certain external mask image, and so batch_spm_svc is a
% streamlined function for that purpose.  It allows the
% batched creation of functional ROI files for many
% subjects, ROIs and contrasts at once, in which each
% functional ROI file is a subject's active-voxel list,
% filtered both by p-threshold and by mask image.
%
% The user supplies a subject list, a contrast list, and an
% ROI list, as well as supplementary arguments, and a
% separate S.V.C. will be performed (and a separate ROI
% written, if desired) for each ROI, for each contrast, for
% each subject.
% 
% Optionally, the script also writes out a text file for
% each input subject containing active-voxel counts for each mask x
% contrast combination for that subject - another common use
% of spm_VOI.
%
% Mask images should be in the same orientation as the
% analysis images (i.e., the beta images and such).
%
% see also: spm_VOI.m.
%_______________________________________________________________________
% batch_spm_svc.m 1.0 Jeff Cooper Jul 2003
% based on
% @(#)spm_VOI.m	2.10 Karl Friston 99/11/17



%-Get input arguments
%-----------------------------------------------------------------------
if nargin < 1 
    expdir = spm_get(-1, '*', 'Select experiment home directory');
else
    expdir = varargin{1};
end

if nargin < 2
    names = spm_get(-Inf, '*', 'Select subject home directories');
    names = cellstr(names);
    names = deblank(names);
    for i = 1:length(names)
        curr = names{i};
        [dummy_dir, names{i}] = fileparts(curr);
    end
else
    names = varargin{2};
    if ischar(names)
        names = cellstr(names);
    end
end

if nargin < 3
    full_resdir = spm_get(-1, '*', 'Select results directory from one subject');
    [dummy_dir, resdir] = fileparts(full_resdir);
    resdir = deblank(resdir);
    if (exist(fullfile(expdir, names{1}, resdir, 'SPM.mat')))==0
        % results directory is not at first level of subject dir
        [test_dummy_dir, test_dir] = fileparts(dummy_dir);
        while(~(strcmp(test_dir, names{1})))
            [dummy_dir, newpart] = fileparts(dummy_dir);
            resdir = fullfile(newpart, resdir);
            [test_dummy_dir, test_dir] = fileparts(dummy_dir);
        end
    end
else
    resdir = varargin{3};
end 
resdir = deblank(resdir);  %I'm being compulsive; I apologize.

if nargin < 4
    con_list = spm_get(Inf, '*.img', 'Select contrast con*.img files from one subject');
    con_list = cellstr(con_list);
    con_list = deblank(con_list);
    for i = 1:length(con_list)
        con_list{i} = spm_str_manip(con_list{i}, 'td');
    end
else
    con_list = varargin{4};
    if ischar(con_list)
        con_list = cellstr(con_list);
    end
end

if nargin < 5
    roi_list = spm_get(Inf, '*.img', 'Select ROI .img files');
    roi_list = cellstr(roi_list);
else
    roi_list = varargin{5};
    if ischar(roi_list)
        roi_list = cellstr(roi_list);
    end
end

if nargin < 6
    p_corrected_flag = spm_input_ui('Correct p-threshold for multiple comparison?',1,'y/n',[1 0],2);
else
    p_corrected_flag = varargin{6};
end

if nargin < 7
    p_thresh = spm_input_ui('Enter threshold (p)','+1','r',0.001,1,[0 1]);
else
    p_thresh = varargin{7};
end

if nargin < 8
    ext_thresh = spm_input_ui('Enter extent threshold (voxels)','+1','i','0', 1);
else
    ext_thresh = varargin{8};
end

if nargin < 9
    save_rois_flag = spm_input_ui('Save small-vol corrected ROI files?', '+1', 'y/n', [1 0], 1);
else
    save_rois_flag = varargin{9};
end

if (nargin < 10) & save_rois_flag
    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',...
            [1:4],1);
elseif nargin < 10
    saveopt = 0;
else
    saveopt = varargin{10};
    if ischar(saveopt)
        switch saveopt
            case 'both'
                saveopt = 1;
            case 'tal'
                saveopt = 2;
            case 'analyze'
                saveopt = 3;
            case 'tal+t'
                saveopt = 4;
        end
    end
end

if nargin < 11
    save_voxcount_flag = spm_input_ui('Save voxel counts?', '+1', 'y/n', [1 0], 1);
else
    save_voxcount_flag = varargin{11};
end

if ~(save_rois_flag | save_voxcount_flag)
    disp('No output required - exiting program...');
    return;  % not asked to return any output, so we're outta here...
end

if nargin < 12
    outdir = spm_get(-1, '*', 'Select output directory');
else
    outdir = varargin{12};
end

% Begin loops
% ----------------------------------------------------------

% begin subject loop
% ------------------
for i = 1:length(names)
    curr_name = names{i};
    fprintf('------------------------\n');
    fprintf('------------------------\n');
    fprintf('Working on subject %s\n', curr_name);
    fprintf('------------------------\n');
    fprintf('------------------------\n\n');
    curr_res_path = fullfile(expdir, curr_name, resdir);
    if (exist(fullfile(curr_res_path, 'SPM.mat'), 'file'))==0;
        full_spm_path = spm_get(1, 'SPM.mat', ['Please locate SPM.mat file for ' curr_name '.'], curr_res_path);
    else
        full_spm_path = fullfile(curr_res_path, 'SPM.mat');
    end
    % initialize vox_count array
    curr_voxcount = zeros(length(roi_list), length(con_list));
    
    % begin contrast loop
    % -------------------
    for j = 1:length(con_list)
        curr_con = deblank(con_list{j});
        load(fullfile(curr_res_path, 'xCon.mat'));
        % this brings xCon into the workspace
        
        % now check out this slick move:
        %    strcmp on an array returns a logical array with
        % zeros everywhere the first argument string is
        % matched, so find should give us the index of
        % whichever contrast has a vCon that matches our
        % desired contrast image.
        curr_con_idx = find(strcmp(curr_con, {xCon.Vcon}));
        % ah, that's hottt.
        
        fprintf('------------------------\n');
        fprintf('Working on contrast %s\n', xCon(curr_con_idx).name);
        fprintf('------------------------\n\n');
        % now go get the SPM and VOL structures we'll need
        % from batch_spm_getSPM
        
        clear SPM VOL
        
        [SPM, VOL] = batch_spm_getSPM(full_spm_path, xCon, curr_con_idx, p_corrected_flag, p_thresh, ext_thresh);
        
        % begin ROI loop
        % --------------
        for k = 1:length(roi_list)
            curr_roi = deblank(roi_list{k});
            
            %these two are re-set for each ROI because their
            %pointlists are edited
            curr_SPM = SPM;
            curr_VOL = VOL;
            
            fprintf('Working on ROI %s\n\n', spm_str_manip(curr_roi, 'rtd'));
            
            % Get voxels in ROI and correct for differences
            % between it and SPM
	        %-----------------------------------------------
	        % 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.
            
        	curr_roi_vol = spm_vol(curr_roi);
            tM    = curr_roi_vol.mat \ curr_VOL.M;
   
            if any(tM([2:5,7:10,12])>eps)
		        spm('alert!',{'ROI image rotated/sheared!','(relative to SPM image)','Can''t use for SVC.'},curr_roi,0);
		        return
	        end
	        XYZ_roivox  = curr_roi_vol.mat \ [curr_SPM.XYZmm; ones(1, size(curr_SPM.XYZmm, 2))];
            %i.e., this is the complete list of XYZ coords
            %in the SPM, but rendered into the space of the
            %ROI image (and in voxels)
	        roi_locs = spm_sample_vol(curr_roi_vol, XYZ_roivox(1,:),XYZ_roivox(2,:),XYZ_roivox(3,:),0) > 0;	
	        
	        %-Compute in-mask volume S:
	        % Correct for differences in mask and SPM voxel sizes
	        roi_Y     = spm_read_vols(curr_roi_vol);
	        vsc   = sqrt(sum(curr_roi_vol.mat(1:3,1:3).^2)) ./ curr_VOL.VOX';
	        S     = sum(roi_Y(:)>0) * prod(vsc);



            %-Select voxels within subspace
            %------------------------------
            curr_SPM.Z     = curr_SPM.Z(roi_locs);
            curr_SPM.XYZ   = curr_SPM.XYZ(:,roi_locs);
            curr_SPM.XYZmm = curr_SPM.XYZmm(:,roi_locs);
            curr_VOL.R     = spm_resels(curr_VOL.FWHM./vsc,curr_roi_vol,'I');
            curr_VOL.S     = S;


            %-If this one's empty, continue
            %------------------------------
            if isempty(curr_SPM.XYZ)
                curr_voxcount(k,j) = 0;
                % rois in rows, contrasts in columns
                fprintf('!!!!!!!!!!!!!!!!!\n');
                fprintf('No voxels survive S.V.C. for this ROI/con/subject combo.\n');
                fprintf('No s.v.c. ROIS will be saved; moving on...\n');
                fprintf('!!!!!!!!!!!!!!!!!\n');
                continue;
            end
           
            %-Assemble voxel data structure
            %------------------------------
            XYZmm = curr_SPM.XYZmm';
            J(:,1) = curr_SPM.Z';
            J(:,2:4) = curr_SPM.XYZmm';  
            %Note that indices for SPM.Z and SPM.XYZ reference same voxel - thus they do 
            %here as well
            J(:,5) = (spm_clusters(curr_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 x = 1:p
                idx = find(J(:,5)==x);
                E(x) = {J(idx,:)};
                %Each E(x) now is a double array consisting of all and only the voxels in 
                %cluster x.  This will do us nicely.
                Q = sortrows(E{x},1);
                Q = flipud(Q);
                E{x} = Q;
                r(x) = 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, x 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 = [];
            XYZmm_tosave_plusT = [];

            for x = 1:p
                curr = E{r(x,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);


            
            % Output files - tal, img
            %-----------------------------------------------
            if save_rois_flag
                
                outroi_name = spm_str_manip(curr_roi, 'urtd');
                %remove leading pathname, suffix;
                %deblank
                outcon_name = spm_str_manip(xCon(curr_con_idx).name, 'ud');
                % remove illegal chars (spaces) and
                % deblank
                outfilename = [curr_name '_' outcon_name '_' outroi_name];
                if ((saveopt == 1) | (saveopt == 2))
                    outfile = fullfile(outdir, [outfilename '.tal']);
                    save(outfile, 'XYZmm_tosave', '-ASCII');
                    fprintf('\n\n----------------------------------------------');
                    fprintf('\nA new file %s was created in',[outfilename '.tal']);
                    fprintf('\n%s\n',outdir);
                end
                if ((saveopt == 1) | (saveopt == 3))
                    outfile = fullfile(outdir, outfilename);
                    Vo = spm_write_filtered(curr_SPM.Z,curr_SPM.XYZ,curr_VOL.DIM,curr_VOL.M,'SPM-filtered',outfile);
                end
                if (saveopt==4)
                    outfile = fullfile(outdir, [outfilename '.tal']);
                    save(outfile, 'XYZmm_tosave_plusT', '-ASCII');
                    fprintf('\n\n----------------------------------------------');
                    fprintf('\nA new file %s was created in',[outfilename '.tal']);
                    fprintf('\n%s\n',outdir);
                end
                
            end % save_rois
            
            
            if save_voxcount_flag
                curr_voxcount(k,j) = size(XYZmm_tosave, 1);
                % rois in rows, contrasts in columns
            end 
            
            
        clear J Q E r p;
        end % ROI loop
    end %contrast loop
    
    % save subject-specific vox-count file
    % ------------------------------------
    if save_voxcount_flag
        voxcount_fname = fullfile(outdir, [curr_name '_voxcount.txt']);
        fid = fopen(voxcount_fname, 'w');
        fprintf(fid, 'Voxel Counts: %s\n\n', curr_name);
        fprintf(fid, 'ROI Name\tContrast Name\tActive Voxels\n');
        for k = 1:size(curr_voxcount,1)
            for j = 1:size(curr_voxcount,2)
                curr_con = deblank(con_list{j});
                curr_con_idx = find(strcmp(curr_con, {xCon.Vcon}));
                fprintf(fid, '%s\t%s\t%g\t\n', spm_str_manip(roi_list{k}, 'td'), ...
                    xCon(curr_con_idx).name, curr_voxcount(k,j));
            end
        end
        fclose(fid);
    end
    
    
end %subject loop