[Marsbar-users] Code to extract data from a single voxel

[Swann P. <swa...@gm...>]

Dear Matthew & MB users, 
The following code has been written with the help of the Marsbar tutorial. 
Although it works fine with ROI files, I am not sure that the way I use 
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 ROI file.mat
        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