method to convert SPM2 design to SPM99 design Heavily based with thanks on work by Jeff Cooper: Written June 2003 by Jeff Cooper (jcooper@stanford.edu) $Id$
0001 function o = mardo_2(o) 0002 % method to convert SPM2 design to SPM99 design 0003 % 0004 % Heavily based with thanks on work by Jeff Cooper: 0005 % Written June 2003 by Jeff Cooper (jcooper@stanford.edu) 0006 % 0007 % $Id$ 0008 0009 % Process design 0010 params = paramfields(o); 0011 SPM99 = params.des_struct; 0012 xX = SPM99.xX; 0013 0014 SPM.xY = struct('P', char(image_names(o)), ... 0015 'VY', get_images(o)); 0016 RT = mars_struct('getifthere', SPM99, 'xX', 'RT'); 0017 if ~isempty(RT), SPM.xY.RT = RT; end 0018 0019 SPM = mars_struct('merge', SPM, ... 0020 mars_struct('split', SPM99, {'xM', 'xGX'})); 0021 SPM.xCon = mars_struct('getifthere', SPM99, 'xCon'); 0022 0023 % This is the first substructure that we have to do some 0024 % actual work to assemble. Much of the details of the 0025 % condition structures (U) can be gleaned from the SPM99 0026 % Sess variable - onset vectors, stick functions, etc. - but 0027 % not quite everything. Trial durations are difficult to figure out, but 0028 % we'll give it a go here with the event_onsets method 0029 0030 if isfield(SPM99, 'Sess'), Sess = SPM99.Sess; else Sess = []; end 0031 for i = 1:length(Sess) 0032 currsess = Sess{i}; 0033 0034 % indices for this session's row in design matrix 0035 SPM.Sess(i).row = currsess.row; 0036 0037 % indices for this session's cols in design matrix 0038 SPM.Sess(i).col = currsess.col; 0039 0040 % number of scans in this session 0041 SPM.nscan(i) = length(currsess.row); 0042 0043 % We can use the length of the onset cell array as an 0044 % effective number of conditions variable - only 'true' 0045 % conditions have onset vectors, not Volterra items or 0046 % user-specified covariates. 0047 for j = 1:length(currsess.ons) 0048 0049 % Try getting onsets and durations from design 0050 [ons dur] = event_onsets(D, [i j]); 0051 0052 % time bin length in seconds 0053 SPM.Sess(i).U(j).dt = xX.dt; 0054 0055 % name of trial type 0056 SPM.Sess(i).U(j).name = {currsess.name{j}}; 0057 0058 % vector of onsets 0059 SPM.Sess(i).U(j).ons = ons; 0060 0061 % duration of trials 0062 SPM.Sess(i).U(j).dur = dur; 0063 0064 % actual stick functions - SPM2 uses 32-bin offset. 0065 SPM.Sess(i).U(j).u = [zeros(32,1); currsess.sf{j}]; 0066 0067 % peristimulus time course (seconds) 0068 SPM.Sess(i).U(j).pst = currsess.pst{j}; 0069 0070 % parameter name 0071 SPM.Sess(i).U(j).P.name = currsess.Pname{j}; 0072 if isempty(currsess.Pname{j}) 0073 % SPM2 includes these values for 'none' 0074 % parameters - just mimicking what they do... 0075 SPM.Sess(i).U(j).P.name = 'none'; 0076 SPM.Sess(i).U(j).P.P = currsess.ons{j}; % values of parameter 0077 SPM.Sess(i).U(j).P.h = 0; % order of polynomial expansion 0078 else 0079 SPM.Sess(i).U(j).P.name = currsess.Pname{j}; % parameter name 0080 SPM.Sess(i).U(j).P.P = currsess.Pv{j}; % parameter values 0081 % leave h and i blank for now... 0082 end 0083 end 0084 % SPM99 doesn't save regressors in independent 0085 % structures as SPM2 does, but it does always append 0086 % them to the end of the session and append 'user 0087 % specified covariates' to the description string of the 0088 % session if they exist, which we can use. 0089 % 0090 % First check and see if there are any: 0091 if ~(isempty(strfind(lower(currsess.DSstr), 'user specified covariates'))) 0092 % i.e., some are in here, so find out what their indices 0093 % are - all the columns that don't correspond to 0094 % specified conditions with onsets and such. 0095 % Note that Volterra items aren't included in 0096 % onsets, either, but they are included in 'the 'name' 0097 % field - so we want to use that as our chopping 0098 % point for finding covariate columns - everything 0099 % after name is a covariate. 0100 reg_indices = currsess.col((length(currsess.name)+1):end); 0101 % extract values from design matrix 0102 SPM.Sess(i).C.C = xX.X(currsess.row, reg_indices); % user-specified covariate values 0103 SPM.Sess(i).name = xX.Xnames(reg_indices); % user-specified covariate names 0104 else 0105 % there are no user specified covariates 0106 SPM.Sess(i).C.C = []; 0107 SPM.Sess(i).C.name = {}; 0108 end 0109 % Everything which is in the 'name' field is either an 0110 % actual condition or a Volterra item, so that's our 0111 % effective number for the Fc array. 0112 for j = 1:length(currsess.name) 0113 SPM.Sess(i).Fc(j).i = currsess.ind{j}; % index for input j and interactions 0114 SPM.Sess(i).Fc(j).name = char(currsess.name(j)); % name for input j and interactions 0115 end 0116 end % Session loop 0117 0118 if ~isempty(Sess) 0119 %%%%%%%%%% 0120 % SPM.xBF - basis function structure 0121 %%%%%%%%%% 0122 0123 % NOTE on the BF structure - where SPM99 allows the 0124 % individual specification of basis functions for each trial 0125 % type, SPM2 just has one basis function set for the whole 0126 % experiment. Since it's difficult to figure out which of 0127 % the various SPM99 basis functions to use on any kind of 0128 % reasoned basis, this program takes the arbitrary step of 0129 % just taking the first one from the first session. This 0130 % should work fine for most experiments, but those with 0131 % mixed basis functions that want to select a different one 0132 % among the ones save are encouraged to modify the lines 0133 % below to get the one to their liking... 0134 sess_idx_for_bf = 1; 0135 cond_idx_for_bf = 1; 0136 0137 SPM.xBF.T = xX.RT / xX.dt; 0138 0139 % # of time bins per scan 0140 % T0 - the reference time bin - isn't saved by SPM99, and 0141 % it's not extractable from the design matrix in any way I 0142 % can figure out. So we'll take a stab at loading it from 0143 % the site's defaults, and if you re-set it, you'd better be 0144 % paying attention... 0145 % first run local defaults file - that'll pull something 0146 % up uniquely. 0147 spm_defaults; 0148 global defaults; %SPM2-style 0149 global fMRI_T0; %SPM99-style 0150 if ~isempty(defaults) 0151 SPM.xBF.T0 = defaults.stats.fmri.t0; % reference time bin 0152 elseif ~isempty(fMRI_T0) 0153 SPM.xBF.T0 = fMRI_T0; 0154 else 0155 % shouldn't get here - there should be some defaults on 0156 % this system - but not a killer error. One they should 0157 % know about, though... 0158 disp('Warning: Defaults don''t contain reference slice value!'); 0159 SPM.xBF.T0 = 1; 0160 end 0161 % SPM99 onsets were reconstructed in scan units 0162 % units that onsets are in ['scans', 'seconds'] 0163 SPM.xBF.UNITS = 'scans'; 0164 if length(SPM.Sess(1).U) < length(SPM.Sess(1).Fc) 0165 % Volterra flag [1=no Volterra, 2=Volterra modeled] 0166 SPM.xBF.Volterra = 2; 0167 else 0168 SPM.xBF.Volterra = 1; 0169 end 0170 SPM.xBF.dt = xX.dt; % time bin length in seconds 0171 SPM.xBF.name = ''; % string: type of basis function - not saved by SPM99 batching 0172 SPM.xBF.bf = Sess{sess_idx_for_bf}.bf(cond_idx_for_bf); 0173 % the actual basis function matrix 0174 % SPM99 doesn't explicitly save length or order of basis 0175 % functions, but this below should work for all non-Fourier, 0176 % non-Gamma basis fns, and maybe even some of them... 0177 SPM.xBF.length = xX.dt*size(SPM.xBF.bf, 1); % window length of basis fn in secs 0178 SPM.xBF.order = size(SPM.xBF.bf,2); % order of basis fn 0179 0180 end 0181 0182 %%%%%%%%%% 0183 % SPM.xVi - temporal non-sphericity struct 0184 %%%%%%%%%% 0185 0186 % We'll be assuming that SPM99 results, in general, made the 0187 % i.i.d. assumption. Estimation of results with the AR(1) 0188 % option specified proceeds very differently in SPM2 than in 0189 % SPM99, and so those that used AR(1) in SPM99 are warned 0190 % that their results will not translate perfectly to SPM2, 0191 % as the whitening matrix W will always be set to identity 0192 % in importing these results. 0193 0194 switch xX.xVi.Form 0195 case 'none' 0196 SPM.xVi.form = 'i.i.d'; % string description of xVi 0197 % covariance constraints - sparse identity matrix in this option for both SPM99 and SPM2. 0198 SPM.xVi.Vi = {xX.xVi.Vi}; 0199 SPM.xVi.V = speye(size(xX.xVi)); % estimated non-sphericity itself. 0200 case 'AR(1)' 0201 SPM.xVi.form = 'AR(0.2)'; % string description of xVi - this is parallel to what SPM2 does 0202 SPM.xVi.Vi = {xX.xVi.Vi} % covariance constraints 0203 0204 % I don't think this will actually make a 0205 % difference, as the whitening matrix W is always 0206 % going to be identity. But it may be useful to 0207 % have. True SPM2 results will also have xVi.V and 0208 % so forth, but not these imported ones... 0209 end 0210 0211 %%%%%%%%%% 0212 % SPM.xX - design matrix structure 0213 %%%%%%%%%% 0214 0215 % Here's the big 'un. Much of this stuff will be translated 0216 % directly from the SPM99 results, and it remains to be seen 0217 % how well that's going to work. We will assume that 0218 % ordinary least squares estimation has been used in SPM99, 0219 % and therefore set W to the identity matrix. 0220 0221 SPM.xX = mars_struct('split', xX, ... 0222 {'X','iH', 'iC', 'iB', 'iG', 'V', 'I' ... 0223 'xKXs', 'pKX', 'Bcov', 'trRV', 'trRVRV', 'erdf', 'nKX'}); 0224 SPM.xX.name = xX.Xnames; 0225 if isfield(xX, 'K') 0226 if iscell(xX.K) 0227 for i = 1:length(xX.K) 0228 SPM.xX.K(i).RT = SPM.xY.RT; % experimental RT 0229 SPM.xX.K(i).row = xX.K{i}.row; % row indices in design matrix for this filter 0230 SPM.xX.K(i).HParam = xX.K{i}.HParam; % high-pass cutoff period in secs 0231 SPM.xX.K(i).X0 = full(xX.K{i}.KH); % frequencies to be removed by filter. 0232 end 0233 else 0234 SPM.xX.K = xX.K; 0235 end 0236 end 0237 SPM.xX.W = speye(size(xX.X,1)); % whitening matrix - identity for ordinary least squares estimates 0238 0239 if isfield(SPM99, 'M') 0240 SPM.xVol = mars_struct('split', SPM99, ... 0241 {'M', 'DIM', 'XYZ', 'S', 'R', 'FWHM'}); 0242 SPM.xVol.iM = inv(SPM.xVol.M); % inverse of M 0243 % We will deal with the RPV image later on... 0244 end 0245 0246 %%%%%%%%%% 0247 % SPM.miscellaneous stuff 0248 %%%%%%%%%% 0249 0250 % various string descriptions of experimental design 0251 SPM.xsDes = SPM99.xsDes; 0252 if isfield(SPM.xsDes, 'Conditions_per_session') 0253 SPM.xsDes.Trials_per_session = SPM.xsDes.Conditions_per_session; 0254 rmfield(SPM.xsDes, 'Conditions_per_session'); 0255 end 0256 SPM.xsDes.Serial_correlations = SPM.xVi.form; 0257 0258 SPM.SPMid = ['SPM2: Results imported from SPM99 design: ' SPM99.SPMid]; 0259 0260 % Try and deal with estimated result volumes 0261 if isfield(SPM99, 'swd') 0262 SPM.swd = SPM99.swd; 0263 if ~exist(SPM99.swd, 'dir') 0264 warning(['Could not find directory: ' SPM.swd]); 0265 else 0266 o_pwd = pwd; 0267 cd(SPM99.swd); 0268 if isfield(SPM99, 'Vbeta') 0269 for i = 1:length(SPM99.Vbeta) 0270 SPM.Vbeta(i) = spm_vol(SPM99.Vbeta{i}); 0271 end 0272 SPM.VResMS = spm_vol(SPM99.VResMS); 0273 SPM.VM = spm_vol(SPM99.VM); 0274 end 0275 for i = 1:length(SPM.xCon) 0276 if ~isempty(SPM.xCon(i).Vcon) 0277 SPM.xCon(i).Vcon = spm_vol(SPM.xCon(i).Vcon); 0278 end 0279 if ~isempty(SPM.xCon(i).Vspm) 0280 SPM.xCon(i).Vspm = spm_vol(SPM.xCon(i).Vspm); 0281 end 0282 end 0283 if isfield(SPM99, 'M') 0284 % Now deal with RPV image 0285 if exist('RPV.img', 'file') 0286 % filehandle of resels per voxel image 0287 SPM.xVol.VRpv = spm_vol('RPV.img'); 0288 else 0289 % Don't really understand this one... 0290 SPM.xVol.VRpv.fname = ['../' SPM.xVol.VRpv.fname]; 0291 end 0292 end 0293 cd(o_pwd); 0294 end 0295 end 0296 0297 % put into parent object 0298 params.des_struct = SPM; 0299 o = mardo_2(params); 0300 0301 return