create_spm99_contrast_struct

 [contrasts] = create_contrast_struct()

 Utility function for populating the contrasts structure that is used by SPM99 

 Previously, this script was part of the *_contrasts.m file associated with an experiment

0001 function [contrasts] = create_contrast_struct(sinfo,cont_name,cont_type,incl_code)
0002 % [contrasts] = create_contrast_struct()
0003 %
0004 % Utility function for populating the contrasts structure that is used by SPM99
0005 %
0006 % Previously, this script was part of the *_contrasts.m file associated with an experiment
0007 
0008 % 12/31/03 PJ
0009 % 09/09/04 PJ  Added analysis path variable
0010 
0011 global rootpath model_proto_name proc_subs analysis_path
0012 
0013 subject_ids = cellstr(char(sinfo(:).id));
0014 nsub_proc = length(proc_subs);
0015 nsub = length(sinfo);
0016 nc = length(cont_name);
0017 
0018 for isub = 1:nsub_proc
0019   co = {};
0020   sub_idx = proc_subs(isub);
0021   
0022   subpath = fullfile(rootpath, subject_ids{sub_idx}, analysis_path, model_proto_name);  
0023   
0024   %
0025   % Get the column names for this subject
0026   %
0027   % Selection of column names can happen in one of two different ways.
0028   % If conditions are modeled with several basis functions, and the contrast is
0029   % specified by the names of the conditions, then the column
0030   % indices must be grabbed from the Sess structure array.
0031   %
0032   % If contrasts are specified with individual column names, the column indices
0033   % are retrieved by matching the names specified in the contrast with the
0034   % names in xX.Xnames
0035   %
0036   % Those names ending with a number in parentheses refer to individual
0037   % columns/basis functions, as listed in xX.Xnames
0038   %
0039   % Column indices associated with names that don't end with a number in
0040   % parentheses are *safely* retrieved from the Sess structure. Note that there
0041   % is a Sess structure for each run
0042   
0043   x = load(fullfile(subpath,'SPM.mat'));
0044   cond_names_detail = x.xX.Xnames;
0045   Sess = x.Sess;
0046   clear x
0047   
0048   % Deal with the Sess mess
0049   cond_names_grouped = [];
0050   tot_ngroup = 0;
0051   for isess = 1:length(Sess)
0052     cond_names_grouped = [cond_names_grouped Sess{isess}.name];
0053     ngroup = length(Sess{isess}.name);
0054     for igroup = 1:ngroup
0055       tot_ngroup = tot_ngroup+1;
0056       grouping_idxs{tot_ngroup} = Sess{isess}.col(Sess{isess}.ind{igroup});
0057     end
0058   end
0059   
0060   cond_names_detail = strvcat(cond_names_detail);
0061   cond_names_detail = cond_names_detail(:,7:end);    % strip leading session identifiers
0062   cond_names_detail = cellstr(cond_names_detail);
0063 
0064   good_cont_vect = [];
0065   ngood = 0;
0066   for icont = 1:nc
0067     switch cont_type{icont}
0068       case 'T'
0069     cond_a_idx = [];
0070     cond_b_idx = [];
0071     % First, try to grab the column index from the detailed list
0072     % This will pick up any regressors and condition components
0073     cond_a_idx = find(ismember(cond_names_detail, incl_code{icont,1}));
0074     cond_b_idx = find(ismember(cond_names_detail, incl_code{icont,2}));
0075     
0076     % Now, grab column indices based on group specifiers
0077     tmp_idx = cat(2,grouping_idxs{find(ismember(cond_names_grouped, ...
0078           incl_code{icont,1}))});
0079     cond_a_idx = [cond_a_idx tmp_idx];
0080     
0081     tmp_idx = cat(2,grouping_idxs{find(ismember(cond_names_grouped, ...
0082           incl_code{icont,2}))});
0083     cond_b_idx = [cond_b_idx tmp_idx];
0084     
0085     cond_a_vect = zeros(1, length(cond_names_detail));
0086     cond_b_vect = zeros(1, length(cond_names_detail));
0087     
0088     cond_a_vect(cond_a_idx) = 1;
0089     cond_b_vect(cond_b_idx) = 1;
0090     
0091     % Weight the contrasts appropriately, i.e. have sum=0
0092     suma = sum(cond_a_vect);
0093     sumb = sum(cond_b_vect);
0094     
0095     if suma & sumb
0096       ngood = ngood+1;
0097       good_cont_vect(icont) = 1;
0098       
0099       weight_b = suma/sumb*-1;
0100       
0101       cond_b_vect = cond_b_vect * weight_b;
0102       
0103       co{ngood} = cond_a_vect + cond_b_vect;
0104     elseif suma
0105       %disp('Only numerator was specified in T contrast')
0106       ngood = ngood+1;
0107       good_cont_vect(icont) = 1;
0108       co{ngood} = cond_a_vect;
0109     else
0110       good_cont_vect(icont) = 0;
0111       disp(sprintf('Skipping T contrast (%s): num_num: %d; num_denom: %d', cont_name{icont},suma, sumb))
0112     end
0113     
0114       case 'F'
0115     col_idx = [];
0116     
0117     col_idx = [col_idx find(ismember(cond_names_detail, incl_code{icont,1}))];
0118     
0119     tmp_idx = cat(2,grouping_idxs{find(ismember(cond_names_grouped, ...
0120           incl_code{icont,1}))});
0121     
0122     col_idx = [col_idx tmp_idx];
0123     
0124     cols = zeros(1,length(cond_names_detail));
0125     cols(col_idx) = 1;
0126     
0127     if any(cols)
0128       ngood = ngood+1;
0129       good_cont_vect(icont) = 1;
0130       co{ngood} = full(sparse(1:sum(cols),find(cols),1));
0131       co{ngood}(end,length(cond_names_detail)) = 0;
0132     else
0133       good_cont_vect(icont) = 0;
0134       disp(sprintf('Skipping F contrast (%s)', cont_name{icont}))
0135     end
0136     
0137       otherwise
0138     error(sprintf('Unknown contrast type: %s', cont_type{icont}))
0139     
0140     end % switch cont_type{icont}
0141     
0142   end % for icont
0143 
0144   good_idxs = find(good_cont_vect);
0145   contrasts(sub_idx).names = cont_name(good_idxs);
0146   contrasts(sub_idx).types = cont_type(good_idxs);
0147   contrasts(sub_idx).values = co;
0148 end