Description of fill

0001 function D = fill(D, actions)
0002 % fills missing entries from SPM FMRI design matrix
0003 % FORMAT D = fill(D, actions)
0004 %
0005 % D          - mardo object containing spm design
0006 % actions    - string or cell array of strings with actions:
0007 %            'defaults' - fills empty parts of design with defaults
0008 %            (in fact this is always done)
0009 %            'filter'  - asks for and fills filter
0010 %            'autocorr' - asks for and fills autocorrelation
0011 %            'for_estimation' fill filter &| autocorr if not present
0012 %            'images'  - asks for and fills with images, mask, scaling
0013 %
0014 % Returns
0015 % D         - returned mardo SPM design
0016 %
0017 % Copied/pasted then rearranged from SPM2 spm_fmri_spm_ui
0018 % Matthew Brett - 17/11/01 - MRS2TH
0019 %
0020 % $Id: fill.m 543 2004-12-14 08:58:05Z matthewbrett $
0021 
0022 if nargin < 2
0023   actions = '';
0024 end
0025 if ~is_fmri(D), return, end
0026 if isempty(actions), actions = {'defaults'}; end
0027 if ischar(actions), actions = {actions}; end
0028 fe = find(ismember(actions, 'for_estimation'));
0029 if ~isempty(fe)
0030   A = [];
0031   if is_fmri(D)
0032     if ~has_filter(D), A = {'filter'}; end
0033     if ~has_autocorr(D), A = [A {'autocorr'}]; end
0034   end
0035   actions(fe) = [];
0036   actions = [actions(1:fe(1)-1) A actions(fe(1):end)];
0037 end
0038 actions = [{'defaults'}, actions];
0039 
0040 % Get design, put into some useful variables
0041 v_f = verbose(D);
0042 SPM = des_struct(D);
0043 xX  = SPM.xX;
0044 have_sess = isfield(SPM, 'Sess');
0045 if have_sess, Sess = SPM.Sess; end
0046 
0047 % get file indices
0048 %---------------------------------------------------------------
0049 row = block_rows(D);
0050 nsess  = length(row);
0051 nscan  = zeros(1,nsess);
0052 for  i = 1:nsess
0053   nscan(i) = length(row{i});
0054 end
0055 
0056 done_list = {};
0057 for a = 1:length(actions)
0058   if ismember(actions{a}, done_list), continue, end
0059   done_list = [actions(a) done_list];
0060   switch lower(actions{a})
0061    case 'defaults'
0062     
0063     % prepare various default settings, offer to design
0064     xM = [];             % masking
0065     xGX = [];            % globals
0066     sGXcalc  = 'none';   % global calculation description
0067     sGMsca   = 'none';   % grand mean scaling description
0068     Global = 'none';     % proportional scaling or no
0069  
0070     BFstr = ''; DSstr = ''; ntr = [];
0071     if have_sess
0072       % Number of trial types per session
0073       for i     = 1:nsess, ntr(i) = length(SPM.Sess(i).U); end
0074       BFstr = SPM.xBF.name;
0075     end
0076     
0077     xsDes = struct(...
0078     'Basis_functions',    BFstr,...
0079     'Number_of_sessions',    sprintf('%d',nsess),...
0080     'Trials_per_session',    sprintf('%-3d',ntr),...
0081     'Global_calculation',    sGXcalc,...
0082     'Grand_mean_scaling',    sGMsca,...
0083     'Global_normalisation',    Global);
0084 
0085     if isfield(SPM, 'xsDes')
0086       xsDes = mars_struct('fillafromb', SPM.xsDes, xsDes);
0087     end
0088     
0089     SPM.xsDes = xsDes;
0090     SPM = mars_struct('merge', SPM, ...
0091                struct('xGX', xGX,...
0092                   'xM',  xM));
0093                   
0094    case 'images'
0095     [Finter,Fgraph,CmdLine] = spm('FnUIsetup','fMRI stats model setup',0);
0096     % get filenames
0097     %---------------------------------------------------------------
0098     P     = [];
0099     for i = 1:nsess
0100       str = sprintf('select scans for session %0.0f',i);
0101       q   = spm_get(nscan(i),mars_veropts('get_img_ext'),str);
0102       P   = strvcat(P,q);
0103     end
0104     
0105     % place in data field
0106     %---------------------------------------------------------------
0107     SPM.xY.P = P;
0108     
0109     % Assemble remaining design parameters
0110     %=======================================================================
0111     
0112     % Global normalization
0113     %-----------------------------------------------------------------------
0114     spm_input('Global intensity normalisation...',1,'d',mfilename)
0115     str             = 'remove Global effects';
0116     SPM.xGX.iGXcalc = spm_input(str,'+1','scale|none',{'Scaling' 'None'});
0117     SPM.xGX.sGXcalc = 'mean voxel value';
0118     SPM.xGX.sGMsca  = 'session specific';
0119     
0120     % Assemble other design parameters
0121     %=======================================================================
0122     spm_help('!ContextHelp',mfilename)
0123     spm_input('Global intensity normalisation...',1,'d',mfilename);
0124     
0125     % get file identifiers and Global values
0126     %=======================================================================
0127     fprintf('%-40s: ','Mapping files')                                   %-#
0128     VY     = spm_vol(SPM.xY.P);
0129     fprintf('%30s\n','...done')                                          %-#
0130     
0131     %-Check compatability of images
0132     %-----------------------------------------------------------------------
0133     [samef msg] = mars_vol_check(VY);
0134     if ~samef, disp(char(msg)),error('Cannot use images'),end;
0135     
0136     %-place in xY
0137     %-----------------------------------------------------------------------
0138     SPM.xY.VY = VY;
0139     
0140     %-Compute Global variate
0141     %=======================================================================
0142     GM    = 100;
0143     q     = length(VY);
0144     g     = zeros(q,1);
0145     fprintf('%-40s: %30s','Calculating globals',' ')                     %-#
0146     for i = 1:q
0147       fprintf('%s%30s',repmat(sprintf('\b'),1,30),sprintf('%4d/%-4d',i,q)) %-#
0148       g(i) = spm_global(VY(i));
0149     end
0150     fprintf('%s%30s\n',repmat(sprintf('\b'),1,30),'...done')               %-#
0151     
0152     % scale if specified (otherwise session specific grand mean scaling)
0153     %-----------------------------------------------------------------------
0154     gSF   = GM./g;
0155     if strcmp(SPM.xGX.iGXcalc,'None')
0156       for i = 1:nsess
0157     gSF(SPM.Sess(i).row) = GM./mean(g(SPM.Sess(i).row));
0158       end
0159     end
0160     
0161     %-Apply gSF to memory-mapped scalefactors to implement scaling
0162     %-----------------------------------------------------------------------
0163     for i = 1:q
0164       SPM.xY.VY(i).pinfo(1:2,:) = SPM.xY.VY(i).pinfo(1:2,:)*gSF(i);
0165     end
0166     
0167     %-place global variates in global structure
0168     %-----------------------------------------------------------------------
0169     SPM.xGX.rg    = g;
0170     SPM.xGX.GM    = GM;
0171     SPM.xGX.gSF   = gSF;
0172     
0173     
0174     %-Masking structure
0175     %---------------------------------------------------------------
0176     SPM.xM     = struct('T',    ones(q,1),...
0177             'TH',    g.*gSF,...
0178             'I',    0,...
0179             'VM',    {[]},...
0180             'xs',    struct('Masking','analysis threshold'));
0181         
0182     xsDes = struct(...
0183     'Global_calculation',    SPM.xGX.sGXcalc,...
0184     'Grand_mean_scaling',    SPM.xGX.sGMsca,...
0185     'Global_normalisation',    SPM.xGX.iGXcalc);
0186       
0187     SPM.xsDes = mars_struct('ffillmerge',...
0188                  SPM.xsDes,...
0189                  xsDes);
0190 
0191    case 'filter'
0192     % Get filter
0193     if ~have_sess, return, end
0194 
0195     [Finter,Fgraph,CmdLine] = spm('FnUIsetup', 'FMRI model filter', 0);
0196 
0197     % TR if not set (it should be)
0198     if ~mars_struct('isthere', SPM, 'xY', 'RT')
0199       SPM.xY.RT  = spm_input('Interscan interval {secs}','+1');
0200     end
0201     SPM.xsDes.Interscan_interval = sprintf('%0.2f {s}',SPM.xY.RT);
0202 
0203     spm_input('High pass filter','+1','d',mfilename)
0204     [SPM.xX.K SPM.xsDes.High_pass_Filter] = ...
0205     pr_get_filter(SPM.xY.RT, SPM.Sess);
0206     
0207    case 'autocorr'
0208     [Finter,Fgraph,CmdLine] = ...
0209     spm('FnUIsetup','FMRI autocorrelation options',0);
0210 
0211     % Contruct Vi structure for non-sphericity ReML estimation
0212     %===============================================================
0213     str   = 'Correct for serial correlations?';
0214     cVi   = {'none', 'SPM AR(0.2)','SPM AR (specify)', 'fmristat AR(n)'};
0215     cVi   = spm_input(str,'+1','m',cVi, cVi);
0216     
0217     % create Vi struct
0218     %-----------------------------------------------------------------------
0219     vox_cov_possible = 0;
0220     switch lower(cVi{1})
0221 
0222      case 'fmristat ar(n)'
0223       % Fit fmristat model AR(n)
0224       %---------------------------------------------------------------
0225       cVi = spm_input('fmristat AR model order', '+1', 'e', 2);
0226       SPM.xVi.Vi = struct('type', 'fmristat', 'order', cVi);
0227       cVi        = sprintf('fmristat AR(%d)',cVi);
0228       f2cl       = 'V'; % Field to CLear
0229       
0230      case 'spm ar (specify)'
0231       % SPM AR coefficient(s) to be specified
0232       %---------------------------------------------------------------
0233       cVi = spm_input('AR rho parameter(s)', '+1', 'e', 0.2);
0234       SPM.xVi.Vi = pr_spm_ce(nscan,cVi);
0235       cVi        = sprintf('AR(%0.1f)',cVi(1));
0236       f2cl       = 'V'; 
0237       vox_cov_possible   = 1;
0238       
0239      case 'none'        
0240       %  xVi.V is i.i.d
0241       %---------------------------------------------------------------
0242       SPM.xVi.V  = speye(sum(nscan));
0243       cVi        = 'i.i.d';
0244       f2cl       = 'Vi'; 
0245                   
0246      otherwise        
0247       % otherwise assume AR(0.2) in xVi.Vi
0248       %---------------------------------------------------------------
0249       SPM.xVi.Vi = pr_spm_ce(nscan,0.2);
0250       cVi        = 'AR(0.2)';
0251       f2cl       = 'V'; 
0252       vox_cov_possible   = 1;
0253       
0254     end
0255 
0256     % If we've set V, need to clear Vi, because the
0257     % estimate method takes the presence of Vi to mean that
0258     % V can be cleared, with 'redo_covar' flag
0259     % Conversely V needs to be cleared if Vi was estimated
0260     if isfield(SPM.xVi, f2cl)
0261       SPM.xVi = rmfield(SPM.xVi, f2cl);
0262       if v_f, fprintf('Clearing previous %s matrix\n', f2cl); end
0263     end
0264     
0265     % Also: remove previous W matrices
0266     % Either will need to be recalculated or won't be used
0267     if isfield(SPM.xX, 'W')
0268       SPM.xX = rmfield(SPM.xX, 'W');
0269       if v_f, fprintf('Clearing previous W matrix\n'); end
0270     end
0271     
0272     % Whether to average covariance estimates over voxels
0273     SPM.xVi.cov_calc = 'summary';
0274     if vox_cov_possible
0275       if spm_input('Use voxel data for covariance','+1','y/n', [1 0], 1);
0276     SPM.xVi.cov_calc = 'vox';
0277       end
0278     end
0279     
0280     % fill into design
0281     SPM.xVi.form = cVi;
0282     xsDes = struct('Serial_correlations', SPM.xVi.form);
0283     SPM.xsDes = mars_struct('ffillmerge', SPM.xsDes, xsDes);
0284   
0285    otherwise
0286     error(['Unpredictable: ' actions{a}]);
0287   end
0288 end
0289 
0290 % put stuff into object
0291 D = des_struct(D,SPM);
fill

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
