Swann Pichon - 2007-11-06

Dear Matthew & MB users,
The following code as been written with the help of the Marsbar tutorial.
Although it works fine with ROI files, I am not sure that the way I wrote it to get stat for a single coord is correct.
The T/F values that I get are quite different from SPM's ones.
I am just beginning to use MB functions, there might be some evident stuff that I am still missing..
Thanks for your indulgence,
Swann

ANALYSIS_dir =  'I:\Data XP\FMRI BEE4\Analyse HRF - scaling slicetiming';
ANALYSIS_SRC_C = {'S1','S2'};
% ROI = {'I:\Data XP\fMRI BEE4\Analyse HRF - scaling slicetiming\Roi_mask\points_coordinate_mm_-58_26_12_roi.mat'};
ROI = {[-58 26 12]};    % MNI coords
UNIT = 'mm'; % coordinates are in voxel space 'vox' or mm space 'mm'
SUMFUNC = 'mean';   % one of 'mean', 'median', 'eigen1', 'wtmean'  to summarize data in the ROIs

res = [];
for roi_file_no=1:length(ROI)
    resStat = [];
    roi_file = ROI{roi_file_no};
   
    for analysis_no=1:length(ANALYSIS_SRC_C)
        work_dirs{analysis_no} = [ANALYSIS_dir filesep ANALYSIS_SRC_C{analysis_no}];
       
        spm_name = [work_dirs{analysis_no} filesep 'SPM.mat'];
        resStat(analysis_no).subj_name = spm_name;
       
        % Make marsbar design object
        D  = mardo(spm_name);
       
        % Make marsbar ROI object
        if ischar(ROI)
            R  = maroi(roi_file);   % From a ROI.mat file
        else
            R = maroi_pointlist(struct('XYZ', ROI{analysis_no}'), UNIT);    % From single MNI coord
        end
       
        % Fetch data into marsbar data object
        Y  = get_marsy(R, D, SUMFUNC);
       
        % Get contrasts from original design
        resStat(analysis_no).xCon = get_contrasts(D);
       
        % Estimate design on ROI data
        resStat(analysis_no).E = estimate(D, Y);
       
        % Put contrasts from original design back into design object
        resStat(analysis_no).E = set_contrasts(resStat(analysis_no).E, resStat(analysis_no).xCon);
       
        % get design betas
        resStat(analysis_no).b = betas(resStat(analysis_no).E);
       
        % get stats and stuff for all contrasts into statistics structure
        resStat(analysis_no).stats = compute_contrasts(resStat(analysis_no).E, 1:length(resStat(analysis_no).xCon));        
    end
    res(roi_file_no).roi_name = roi_file;
    res(roi_file_no).stat = resStat;
end