Home > spm_vbglmar_slice > spm_contrasts.m

spm_contrasts

PURPOSE ^

Fills in SPM.xCon and writes con_????.img and SPM?_????.img

SYNOPSIS ^

function [SPM] = spm_contrasts(SPM,Ic)

DESCRIPTION ^

 Fills in SPM.xCon and writes con_????.img and SPM?_????.img
 FORMAT [SPM] = spm_contrasts(SPM,Ic);
 Ic  - indices of xCon to compute
_______________________________________________________________________
 @(#)spm_contrasts.m    2.3 Andrew Holmes, Karl Friston & Jean-Baptiste Poline 02/12/30

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [SPM] = spm_contrasts(SPM,Ic)
0002 % Fills in SPM.xCon and writes con_????.img and SPM?_????.img
0003 % FORMAT [SPM] = spm_contrasts(SPM,Ic);
0004 % Ic  - indices of xCon to compute
0005 %_______________________________________________________________________
0006 % @(#)spm_contrasts.m    2.3 Andrew Holmes, Karl Friston & Jean-Baptiste Poline 02/12/30
0007 
0008 
0009 %-Get and change to results directory
0010 %-----------------------------------------------------------------------
0011 try
0012     swd   = SPM.swd;
0013     cd(swd)
0014 catch
0015     swd   = pwd;
0016 end
0017 
0018 %-Get contrast definitions (if available)
0019 %-----------------------------------------------------------------------
0020 try
0021     xCon  = SPM.xCon;
0022 catch
0023     xCon  = [];
0024 end
0025 
0026 %-set all contrasts by default
0027 %-----------------------------------------------------------------------
0028 if nargin < 2
0029     Ic    = 1:length(xCon);
0030 end
0031 
0032 
0033 % map parameter and hyperarameter files
0034 %-----------------------------------------------------------------------
0035 if xCon(Ic(1)).STAT == 'P' 
0036 
0037     % Conditional estimators and error variance hyperparameters
0038     %----------------------------------------------------------------
0039     if isfield(SPM.PPM,'VB')
0040         Vbeta = SPM.VCbeta;
0041         VHp   = SPM.VHp;
0042         VPsd  = SPM.VPsd;
0043         VAR   = SPM.VAR;
0044     else
0045         Vbeta = SPM.VCbeta;
0046         VHp   = SPM.VHp;
0047     end
0048 
0049 else
0050     % OLS estimators and error variance estimate
0051     %----------------------------------------------------------------
0052     Vbeta = SPM.Vbeta;
0053     VHp   = SPM.VResMS;
0054 end
0055 
0056 
0057 
0058 %-Compute & store contrast parameters, contrast/ESS images, & SPM images
0059 %=======================================================================
0060 spm('Pointer','Watch')
0061 XYZ   = SPM.xVol.XYZ;
0062 for i = 1:length(Ic)
0063 
0064 
0065     %-Canonicalise contrast structure with required fields
0066     %-------------------------------------------------------------------
0067     ic  = Ic(i);
0068     if isempty(xCon(ic).eidf)
0069     X1o           = spm_FcUtil('X1o',xCon(ic),SPM.xX.xKXs);
0070     [trMV,trMVMV] = spm_SpUtil('trMV',X1o,SPM.xX.V);
0071         xCon(ic).eidf = trMV^2/trMVMV;
0072     end
0073 
0074 
0075     %-Write contrast/ESS images?
0076     %-------------------------------------------------------------------
0077     if isempty(xCon(ic).Vcon)
0078  
0079     switch(xCon(ic).STAT)
0080 
0081         case {'T','P'} %-Implement contrast as sum of scaled beta images
0082             %-----------------------------------------------------------
0083             fprintf('\t%-32s: %-10s%20s',sprintf('contrast image %2d',i),...
0084                                       '(spm_add)','...initialising') %-#
0085 
0086         Q     = find(abs(xCon(ic).c) > 0);
0087         V     = Vbeta(Q);
0088         for j = 1:length(Q)
0089             V(j).pinfo(1:2,:) = V(j).pinfo(1:2,:)*xCon(ic).c(Q(j));
0090         end
0091         
0092         %-Prepare handle for contrast image
0093         %-----------------------------------------------------------
0094         xCon(ic).Vcon = struct(...
0095             'fname',  sprintf('con_%04d.img',ic),...
0096                 'dim',    [SPM.xVol.DIM',16],...
0097                 'mat',    SPM.xVol.M,...
0098                 'pinfo',  [1,0,0]',...
0099                 'descrip',sprintf('SPM contrast - %d: %s',ic,xCon(ic).name));
0100 
0101         %-Write image
0102         %-----------------------------------------------------------
0103             fprintf('%s%20s',sprintf('\b')*ones(1,20),'...computing')%-#
0104             xCon(ic).Vcon            = spm_create_vol(xCon(ic).Vcon);
0105             xCon(ic).Vcon.pinfo(1,1) = spm_add(V,xCon(ic).Vcon);
0106         xCon(ic).Vcon            = spm_close_vol(xCon(ic).Vcon);
0107             xCon(ic).Vcon            = spm_create_vol(xCon(ic).Vcon,'noopen');
0108         xCon(ic).Vcon            = spm_close_vol(xCon(ic).Vcon);
0109             
0110             fprintf('%s%30s\n',sprintf('\b')*ones(1,30),sprintf(...
0111                 '...written %s',spm_str_manip(xCon(ic).Vcon.fname,'t')))%-#
0112 
0113 
0114 
0115         case 'F'  %-Implement ESS as sum of squared weighted beta images
0116             %-----------------------------------------------------------
0117             fprintf('\t%-32s: %30s',sprintf('ESS image %2d',i),...
0118                                                      '...computing') %-#
0119 
0120             %-Residual (in parameter space) forming mtx
0121         %-----------------------------------------------------------
0122             h       = spm_FcUtil('Hsqr',xCon(ic),SPM.xX.xKXs);
0123 
0124         %-Prepare handle for ESS image
0125         %-----------------------------------------------------------
0126         xCon(ic).Vcon = struct(...
0127             'fname',  sprintf('ess_%04d.img',ic),...
0128                 'dim',    [SPM.xVol.DIM',16],...
0129                 'mat',    SPM.xVol.M,...
0130                 'pinfo',  [1,0,0]',...
0131                 'descrip',sprintf('SPM ESS -contrast %d: %s',ic,xCon(ic).name));
0132 
0133             %-Write image
0134         %-----------------------------------------------------------
0135             fprintf('%s',sprintf('\b')*ones(1,30))                   %-#
0136             xCon(ic).Vcon = spm_create_vol(xCon(ic).Vcon);
0137             xCon(ic).Vcon = spm_resss(Vbeta,xCon(ic).Vcon,h);
0138         xCon(ic).Vcon = spm_close_vol(xCon(ic).Vcon);
0139             xCon(ic).Vcon = spm_create_vol(xCon(ic).Vcon,'noopen');
0140         xCon(ic).Vcon = spm_close_vol(xCon(ic).Vcon);
0141 
0142 
0143     otherwise
0144         %---------------------------------------------------------------
0145         error(['unknown STAT "',xCon(ic).STAT,'"'])
0146 
0147     end % (switch(xCon...)
0148 
0149     end % (if ~isfield...)
0150 
0151     
0152     
0153     
0154     %-Write inference SPM/PPM
0155     %===================================================================
0156     if isempty(xCon(ic).Vspm)
0157     
0158         fprintf('\t%-32s: %30s',sprintf('spm{%c} image %2d',xCon(ic).STAT,i),...
0159                                                     '...computing')  %-#
0160 
0161     switch(xCon(ic).STAT)
0162 
0163     case 'T'                                  %-Compute SPM{t} image
0164         %---------------------------------------------------------------
0165     cB    = spm_get_data(xCon(ic).Vcon,XYZ);
0166     l     = spm_get_data(VHp,XYZ);
0167     VcB   = xCon(ic).c'*SPM.xX.Bcov*xCon(ic).c; 
0168     Z     = cB./sqrt(l*VcB);
0169     str   = sprintf('[%.1f]',SPM.xX.erdf);
0170     
0171 
0172     case 'P'                                  %-Compute PPM{P} image
0173         %---------------------------------------------------------------
0174     c     = xCon(ic).c;
0175     cB    = spm_get_data(xCon(ic).Vcon,XYZ);
0176     if isfield(SPM.PPM,'VB');
0177         % Get approximate posterior covariance at ic
0178         % using Taylor-series approximation
0179         
0180         % Get posterior SD beta's
0181         Nk=size(SPM.xX.X,2);
0182         for k=1:Nk,
0183             sd_beta(k,:) = spm_get_data(VPsd(k),XYZ);
0184         end
0185         
0186         % Get AR coefficients
0187         for p=1:SPM.PPM.AR_P
0188             a(p,:) = spm_get_data(VAR(p),XYZ);
0189         end
0190         
0191         % Get noise SD
0192         lambda = spm_get_data(VHp,XYZ);
0193         
0194         % Loop over voxels
0195         Nvoxels=size(XYZ,2);
0196         for v=1:Nvoxels,
0197             % Which slice are we in ?
0198             slice_index=XYZ(3,v);
0199             
0200             % Reconstuct approximation to voxel wise correlation matrix
0201             R=SPM.PPM.slice(slice_index).mean.R;
0202             dh=a(:,v)'-SPM.PPM.slice(slice_index).mean.a;
0203             dh=[dh lambda(v)-SPM.PPM.slice(slice_index).mean.lambda];
0204             for i=1:length(dh),
0205                 R=R+SPM.PPM.slice(slice_index).mean.dR(:,:,i)*dh(i);
0206             end 
0207             % Reconstuct approximation to voxel wise covariance matrix
0208             Sigma_post = (sd_beta(:,v)*sd_beta(:,v)').*R;
0209             
0210             VcB (v) = c'*Sigma_post*c;
0211         end
0212             
0213   
0214     else
0215         VcB   = c'*SPM.PPM.Cby*c;
0216         for j = 1:length(SPM.PPM.l)
0217             
0218             % hyperparameter and Taylor approximation
0219             %-------------------------------------------------------
0220             l   = spm_get_data(VHp(j),XYZ);
0221             VcB = VcB + (c'*SPM.PPM.dC{j}*c)*(l - SPM.PPM.l(j));
0222         end
0223     end
0224     
0225     % posterior probability cB > g
0226         %---------------------------------------------------------------
0227      Gamma         = xCon(ic).eidf;
0228     Z             = 1 - spm_Ncdf(Gamma,cB,VcB);
0229     str           = sprintf('[%.2f]',Gamma);
0230     xCon(ic).name = [xCon(ic).name ' ' str];
0231 
0232 
0233     case 'F'                                  %-Compute SPM{F} image
0234         %---------------------------------------------------------------
0235     MVM   = spm_get_data(xCon(ic).Vcon,XYZ)/trMV;
0236     RVR   = spm_get_data(VHp,XYZ);
0237     Z     = MVM./RVR;
0238         str   = sprintf('[%.1f,%.1f]',xCon(ic).eidf,SPM.xX.erdf);
0239 
0240     otherwise
0241         %---------------------------------------------------------------
0242         error(['unknown STAT "',xCon(ic).STAT,'"'])
0243     end
0244 
0245 
0246         %-Write SPM - statistic image
0247         %---------------------------------------------------------------
0248         fprintf('%s%30s',sprintf('\b')*ones(1,30),'...writing')      %-#
0249 
0250         xCon(ic).Vspm = struct(...
0251         'fname',  sprintf('spm%c_%04d.img',xCon(ic).STAT,ic),...
0252         'dim',    [SPM.xVol.DIM',16],...
0253         'mat',    SPM.xVol.M,...
0254         'pinfo',  [1,0,0]',...
0255         'descrip',sprintf('SPM{%c_%s} - contrast %d: %s',...
0256                            xCon(ic).STAT,str,ic,xCon(ic).name));
0257 
0258         tmp           = zeros(SPM.xVol.DIM');
0259     Q             = cumprod([1,SPM.xVol.DIM(1:2)'])*XYZ - ...
0260                    sum(cumprod(SPM.xVol.DIM(1:2)'));
0261     tmp(Q)        = Z;
0262     xCon(ic).Vspm = spm_write_vol(xCon(ic).Vspm,tmp);
0263 
0264 
0265     clear tmp Z
0266         fprintf('%s%30s\n',sprintf('\b')*ones(1,30),sprintf(...
0267             '...written %s',spm_str_manip(xCon(ic).Vspm.fname,'t')))  %-#
0268 
0269     end % (if ~isfield...)
0270 
0271 end % (for i = 1:length(Ic))
0272 
0273 
0274 % place xCon back in SPM
0275 %-----------------------------------------------------------------------
0276 SPM.xCon = xCon;
0277 save SPM SPM
Generated on Mon 23-Aug-2004 14:59:38 by m2html © 2003