Description of pr

0001 function SPM = pr_estimate(SPM, VY, params)
0002 % Wavelet analysis enabled version of spm_spm.m
0003 % Estimation of the General Linear Model
0004 % FORMAT SPM = pr_estimate(SPM, VY, params)
0005 %
0006 % Based on:
0007 % @(#)spm_spm.m    2.37 Andrew Holmes, Jean-Baptiste Poline, Karl Friston 99/11/26
0008 %
0009 % The primary differences from the SPM99 version are:
0010 % 1) No smoothness estimation (and no RPV output image)
0011 % 2) Output file name is [wtprefix 'phiw_spm.mat']
0012 % 3) Beta, ResMS image have wtprefix appended
0013 % 4) xCon is not saved separately, but in output file
0014 % 5) Design returned in SPM output structure
0015 % 6) Option to write out residual images during estimation
0016 %
0017 % VY    - nScan x nVar struct array of mapped image volumes
0018 %         Images must have the same orientation, voxel size and data type
0019 %       - Any scaling should have already been applied via the image handle
0020 %         scalefactors.
0021 %       - defaults to images stored in design
0022 %
0023 % params - structure containing options
0024 %       'write_res'  if 1, write residual images
0025 %
0026 % For detailed help on the mathematics, structures etc, see spm_spm.m in
0027 % the SPM99 distribution - version string above.
0028 %
0029 % $Id: pr_estimate.m,v 1.9 2005/07/01 18:02:06 matthewbrett Exp $
0030 
0031 %-Condition arguments
0032 %-----------------------------------------------------------------------
0033 def_params = struct('write_res', 1);
0034 
0035 if nargin < 1
0036   error('Need SPM design structure');
0037 end
0038 if nargin < 2
0039   VY = [];
0040 end
0041 if isempty(VY)
0042   if ~isfield(SPM, 'VY')
0043     error('Need data to estimate for');
0044   end
0045   VY = SPM.VY;
0046 end
0047 if nargin < 3
0048   params = [];
0049 end
0050 params = mars_struct('ffillsplit', def_params, params);
0051 
0052 if ~isfield(SPM, 'xX'), error('Need design in SPM structure'); end
0053 xX = SPM.xX;
0054 if ~isfield(SPM, 'xM')
0055   xM = zeros(size(xX.X,1),1);
0056 else
0057   xM = SPM.xM;
0058 end
0059 if ~isfield(SPM, 'F_iX0')
0060   F_iX0 = []; 
0061 else
0062   F_iX0 = SPM.F_iX0;
0063 end
0064 
0065 % Check images are wt'ed
0066 wvobj = phiw_wvimg(VY(1), struct('noproc', 1));
0067 if isempty(wvobj), error('Images do not appear to have been WTed'), end
0068 wti = wtinfo(wvobj);
0069 wtp = wti.wtprefix;
0070 phiw_mat_name = [wtp 'phiw_spm.mat'];
0071 
0072 % We must set SPMid to contain SPM99 string in order for the mardo_99 object
0073 % to recognize this as an SPM99 design
0074 SPMid  = sprintf('SPM99: Phiwave estimation; %s version %s', ...
0075          mfilename, ...
0076          mars_cvs_version(mfilename, 'phido_99'));
0077 
0078 % Backspace macro
0079 bs30 = repmat(sprintf('\b'),1,30);
0080 
0081 %-Say hello
0082 %-----------------------------------------------------------------------
0083 Finter   = spm('FigName','Stats: estimation...'); spm('Pointer','Watch')
0084 
0085 %-Parameters
0086 %-----------------------------------------------------------------------
0087 maxMem   = 2^30;    %-Max data chunking size, in bytes
0088 
0089 %-Check required fields of xX structure exist - default optional fields
0090 %-----------------------------------------------------------------------
0091 for tmp = {'X','Xnames'}
0092     if ~isfield(xX,tmp)
0093              error(sprintf('xX doesn''t contain ''%s'' field',tmp{:}))
0094     end
0095 end
0096 if ~isfield(xX,'K')
0097     xX.K  = speye(size(xX.X,1));
0098 end
0099 if ~isfield(xX,'xVi')
0100     xX.xVi = struct(    'Vi',    speye(size(xX.X,1)),...
0101                 'Form',    'none'); 
0102 end
0103 
0104 %-If xM is not a structure then assumme it's a vector of thresholds
0105 %-----------------------------------------------------------------------
0106 if ~isstruct(xM), xM = struct(    'T',    [],...
0107                 'TH',    xM,...
0108                 'I',    0,...
0109                 'VM',    {[]},...
0110                 'xs',    struct('Masking','analysis threshold'));
0111 end
0112 
0113 %-Delete files from previous analyses
0114 %-----------------------------------------------------------------------
0115 if exist(fullfile('.',phiw_mat_name),'file')==2
0116     spm('alert!',{...
0117         'Current directory contains some phiwave stats results files!',...
0118             ['        (pwd = ',pwd,')'],...
0119         'Existing results are being overwritten!'},...
0120         mfilename,1);
0121     warning(sprintf('Overwriting existing results\n\t (pwd = %s) ',pwd))
0122     drawnow
0123 end
0124 
0125 files = { 'mask.???','ResMS.???','ResI_?????.???'...
0126       'beta_????.???','con_????.???',...
0127       'ess_????.???', 'spm?_????.???'};
0128 for i=1:length(files), files{i} = [wtp files{i}]; end
0129 files{end+1} = phiw_mat_name;
0130 
0131 for i=1:length(files)
0132     if any(files{i} == '*'|files{i} == '?' )
0133         [tmp,null] = spm_list_files(pwd,files{i});
0134         for i=1:size(tmp,1)
0135             spm_unlink(deblank(tmp(i,:)))
0136         end
0137     else
0138         spm_unlink(files{i})
0139     end
0140 end
0141 
0142 %-Check & note Y images dimensions, orientation & voxel size
0143 %-----------------------------------------------------------------------
0144 if any(any(diff(cat(1,VY(:).dim),1,1),1) & [1,1,1,0]) 
0145     error('images do not all have the same dimensions')
0146 end
0147 if any(any(any(diff(cat(3,VY(:).mat),1,3),3)))
0148     error('images do not all have same orientation & voxel size')
0149 end
0150 
0151 M      = VY(1,1).mat;
0152 DIM    = VY(1,1).dim(1:3)';
0153 N      = 3 - sum(DIM == 1);
0154 
0155 
0156 %=======================================================================
0157 % - A N A L Y S I S   P R E L I M I N A R I E S
0158 %=======================================================================
0159 
0160 %-Initialise design space
0161 %=======================================================================
0162 fprintf('%-40s: %30s','Initialising design space','...computing')    %-#
0163 
0164 %-Construct design parmameters, and store in design structure xX
0165 % Take care to apply temporal convolution - KX stored as xX.xKX.X
0166 %-Note that Vi may not be known exactly at this point, if it is to be
0167 % estimated. Parameters dependent on Vi are committed to xX at the end.
0168 %-Note that the default F-contrast (used to identify "interesting" voxels
0169 % to save raw data for) computation requires Vi. Thus, if Vi is to be
0170 % estimated, any F-threshold will only be have upper tail probability UFp.
0171 %-----------------------------------------------------------------------
0172 % V            - Autocorrelation matrix K*Vi*K'
0173 % xX.xKXs      - Design space structure of KX
0174 % xX.pKX       - Pseudoinverse of KX
0175 % trRV,trRVRV  - Variance expectations
0176 % erdf         - Effective residual d.f.
0177 %-----------------------------------------------------------------------
0178 [nScan nBeta] = size(xX.X);
0179 [nScan nVar]  = size(VY);
0180 KVi           = mardo_99('spm_filter', 'apply',xX.K, xX.xVi.Vi);
0181 V             = mardo_99('spm_filter', 'apply',xX.K,KVi');
0182 xX.xKXs       = spm_sp('Set',mardo_99('spm_filter', 'apply',xX.K, xX.X));
0183 xX.pKX        = spm_sp('x-',xX.xKXs);
0184 [trRV trRVRV] = spm_SpUtil('trRV',xX.xKXs,V);
0185 erdf          = trRV^2/trRVRV;
0186 
0187 
0188 %-Check estimability
0189 %-----------------------------------------------------------------------
0190 if  erdf < 0
0191     error(sprintf('This design is completely unestimable! (df=%-.2g)',erdf))
0192 elseif erdf == 0
0193     error('This design has no residuals! (df=0)')
0194 elseif erdf < 4
0195     warning(sprintf('Very low degrees of freedom (df=%-.2g)',erdf))
0196 end
0197 
0198 
0199 %-Default F-contrasts (in contrast structure) & Y.mad pointlist filtering
0200 %=======================================================================
0201 fprintf('%s%30s',bs30,'...F-contrast')           %-#
0202 
0203 if isempty(F_iX0)
0204     F_iX0 = struct(    'iX0',        [],...
0205             'name',        'all effects');
0206 elseif ~isstruct(F_iX0)
0207     F_iX0 = struct(    'iX0',        F_iX0,...
0208             'name',        'effects of interest');
0209 end
0210 
0211 %-Create Contrast structure array
0212 %-----------------------------------------------------------------------
0213 xCon  = spm_FcUtil('Set',F_iX0(1).name,'F','iX0',F_iX0(1).iX0,xX.xKXs);
0214 for i = 2:length(F_iX0)
0215   xcon = spm_FcUtil('Set',F_iX0(i).name,'F','iX0',F_iX0(i).iX0,xX.xKXs);
0216   xCon = [xCon xcon];
0217 end
0218 
0219 %-Parameters for saving in Y.mad (based on first F-contrast)
0220 %-----------------------------------------------------------------------
0221 [trMV trMVMV] = spm_SpUtil('trMV',spm_FcUtil('X1o',xCon(1),xX.xKXs),V);
0222 eidf          = trMV^2/trMVMV;
0223 h             = spm_FcUtil('Hsqr',xCon(1),xX.xKXs);
0224 
0225 %-Modify structures for multivariate inference
0226 %-----------------------------------------------------------------------
0227 if nVar > 1
0228 
0229   fprintf('%s%30s',bs30,'...multivariate prep')%-#
0230   
0231   % pseudoinverse of null partition KX0
0232   %---------------------------------------------------------------
0233   KX0         = spm_FcUtil('X0',xCon(1),xX.xKXs);
0234   pKX0        = pinv(KX0);
0235   
0236   %-Modify Contrast structure for multivariate inference
0237   %---------------------------------------------------------------
0238   str       = 'Canonical variate';
0239   xCon      = spm_FcUtil('Set',str,'F','iX0',xCon.iX0,xX.xKXs);
0240   
0241   %-Degrees of freedom (Rao 1951)
0242   %---------------------------------------------------------------
0243   h         = rank(spm_FcUtil('X1o',xCon(1),xX.xKXs));
0244   p         = nVar;
0245   r         = erdf;
0246   a         = r - (p - h + 1)/2;
0247   if (p + h) == 3;
0248     b = 1;
0249   else
0250     b = sqrt((p^2 * h^2 - 4)/(p^2 + h^2 - 5));
0251   end
0252   c         = (p*h - 2)/2;
0253   erdf      = a*b - c;
0254   eidf      = p*h;
0255   xCon.eidf = eidf;
0256 end
0257 
0258 fprintf('%s%30s\n',bs30,'...done')               %-#
0259 
0260 %-Initialise output images
0261 %=======================================================================
0262 fprintf('%-40s: %30s','Output images','...initialising')             %-#
0263 
0264 %-Image dimensions
0265 %-----------------------------------------------------------------------
0266 xdim    = DIM(1); ydim = DIM(2); zdim = DIM(3);
0267 YNaNrep = spm_type(VY(1,1).dim(4),'nanrep');
0268 
0269 %-Intialise the name of the new mask : current mask & conditions on voxels
0270 %-----------------------------------------------------------------------
0271 VM    = struct(        'fname',    [wtp 'mask.img'],...
0272             'dim',        [DIM',spm_type('uint8')],...
0273             'mat',        M,...
0274             'pinfo',    [1 0 0]',...
0275             'descrip',    'phiwave:resultant analysis mask');
0276 VM    = spm_create_image(VM);
0277 
0278 
0279 %-Intialise beta image files
0280 %-----------------------------------------------------------------------
0281 Vbeta_tmp = deal(struct(...
0282             'fname',    [],...
0283             'dim',        [DIM',spm_type('float')],...
0284             'mat',        M,...
0285             'pinfo',    [1 0 0]',...
0286             'descrip',    ''));
0287 
0288 for i = 1:nBeta
0289     Vbeta_tmp.fname   = sprintf('%sbeta_%04d.img', wtp, i);
0290     Vbeta_tmp.descrip = sprintf('phiwave:beta (%04d) - %s',i,xX.Xnames{i});
0291     spm_unlink(Vbeta_tmp.fname)
0292     Vbeta(i)         = spm_create_image(Vbeta_tmp);
0293 end
0294 
0295 
0296 %-Intialise residual sum of squares image file
0297 %-----------------------------------------------------------------------
0298 VResMS = struct(    'fname',    [wtp 'ResMS.img'],...
0299             'dim',        [DIM',spm_type('double')],...
0300             'mat',        M,...
0301             'pinfo',    [1 0 0]',...
0302             'descrip',    'phiwave:Residual sum-of-squares');
0303 VResMS = spm_create_image(VResMS);
0304 
0305 %-Intialise residual images
0306 %-----------------------------------------------------------------------
0307 if params.write_res
0308   VResI_shell = struct(...
0309       'fname',  [],...
0310       'dim',    [DIM',spm_type('double')],...
0311       'mat',    M,...
0312       'pinfo',  [1 0 0]',...
0313       'descrip','spm_spm:Residual image');
0314   
0315   for i = 1:nScan
0316     vr_i = VResI_shell;
0317     vr_i.fname   = sprintf('%sResI_%05d.img', wtp, i);
0318     vr_i.descrip = sprintf('spm_spm:ResI (%04d)', i);
0319     spm_unlink(vr_i.fname);
0320     VResI(i) = spm_create_image(vr_i);
0321   end
0322 else
0323   VResI = [];
0324 end
0325 
0326 %-Intialise multivariate SPMF file
0327 %-----------------------------------------------------------------------
0328 if nVar > 1
0329     Vspm   = struct('fname',    [wtp 'mvSPMF.img'],...
0330             'dim',        [DIM',spm_type('double')],...
0331             'mat',        M,...
0332             'pinfo',    [1 0 0]',...
0333             'descrip',    'phiwave:multivariate F');
0334     Vspm   = spm_create_image(Vspm);
0335 end
0336 
0337 fprintf('%s%30s\n',bs30,'...initialised')        %-#
0338 
0339 
0340 %=======================================================================
0341 % - F I T   M O D E L   &   W R I T E   P A R A M E T E R    I M A G E S
0342 %=======================================================================
0343 
0344 %-Find a suitable block size for the main loop, which proceeds a bunch
0345 % of lines at a time (minimum = one line; maximum = one plane)
0346 % (maxMem is the maximum amount of data that will be processed at a time)
0347 %-----------------------------------------------------------------------
0348 blksz    = maxMem/8/nScan/nVar;            %-block size (in bytes)
0349 if ydim < 2, error('ydim < 2'), end        %-need at least 2 lines
0350 nl     = max(min(round(blksz/xdim),ydim),1);     %-max # lines / block
0351 clines    = 1:nl:ydim;                %-bunch start line #'s
0352 blines  = diff([clines ydim+1]);        %-#lines per bunch
0353 nbch    = length(clines);            %-#bunches
0354 
0355 
0356 %-Intialise other variables used through the loop
0357 %=======================================================================
0358 BePm        = zeros(1,xdim*ydim);            %-below plane (mask)
0359 CrVox       = struct('res',[],'ofs',[],'ind',[]);   %-current voxels
0360                             % res : residuals
0361                                % ofs : mask offset
0362                             % ind : indices
0363 
0364 xords  = [1:xdim]'*ones(1,ydim); xords = xords(:)'; %-plane X coordinates
0365 yords  = ones(xdim,1)*[1:ydim];  yords = yords(:)'; %-plane Y coordinates
0366 
0367 %-Initialise XYZ matrix of in-mask voxel co-ordinates (real space)
0368 %-----------------------------------------------------------------------
0369 XYZ   = [];
0370 
0371 %-Smoothness estimation variables
0372 %-----------------------------------------------------------------------
0373 S      = 0;                                     % Volume analyzed (in voxels)
0374 
0375 %-parameter for estimation of intrinsic correlations AR(1) model
0376 %-----------------------------------------------------------------------
0377 A      = 0;                    %-regression coeficient
0378 
0379 %-Cycle over bunches of lines within planes (planks) to avoid memory problems
0380 %=======================================================================
0381 spm_progress_bar('Init',100,'model estimation','');
0382 
0383 for z = 1:zdim                %-loop over planes (2D or 3D data)
0384 
0385     zords   = z*ones(xdim*ydim,1)';    %-plane Z coordinates
0386     CrBl    = [];            %-current plane betas
0387     CrResI  = [];            %-normalized residuals
0388     CrmvF   = [];            %-current plane mvF-squared
0389     CrResSS = [];            %-current plane ResSS
0390 
0391     for bch = 1:nbch            %-loop over bunches of lines (planks)
0392 
0393     %-# Print progress information in command window
0394     %---------------------------------------------------------------
0395     fprintf('\r%-40s: %30s',sprintf('Plane %3d/%-3d, plank %3d/%-3d',...
0396         z,zdim,bch,nbch),' ')                                %-#
0397 
0398     cl    = clines(bch);          %-line index of first line of bunch
0399     bl    = blines(bch);          %-number of lines for this bunch
0400 
0401     %-construct list of voxels in this bunch of lines
0402     %---------------------------------------------------------------
0403     I     = ((cl-1)*xdim+1):((cl+bl-1)*xdim);    %-lines cl:cl+bl-1
0404     xyz   = [xords(I); yords(I); zords(I)];        %-voxel coords in bch
0405 
0406 
0407     %-Get data & construct analysis mask for this bunch of lines
0408     %===============================================================
0409     fprintf('%s%30s',bs30,'...read & mask data')%-#
0410     CrLm    = logical(ones(1,xdim*bl));        %-current lines mask
0411     CrLmxyz = zeros(size(CrLm));            %-and for smoothnes
0412 
0413     %-Compute explicit mask
0414     % (note that these may not have same orientations)
0415     %---------------------------------------------------------------
0416     for i = 1:length(xM.VM)
0417         tM   = inv(xM.VM(i).mat)*M;        %-Reorientation matrix
0418         tmp  = tM * [xyz;ones(1,size(xyz,2))];    %-Coords in mask image
0419 
0420         %-Load mask image within current mask & update mask
0421         %-------------------------------------------------------
0422         CrLm(CrLm) = spm_sample_vol(xM.VM(i),...
0423                 tmp(1,CrLm),tmp(2,CrLm),tmp(3,CrLm),0) > 0;
0424     end
0425     
0426     %-Get the data in mask, compute threshold & implicit masks
0427     %---------------------------------------------------------------
0428     Y     = zeros(nScan,xdim*bl);
0429     for i = 1:nScan
0430         if ~any(CrLm), break, end        %-Break if empty mask
0431         Y(i,CrLm)  = spm_sample_vol(VY(i,1),... %-Load data in mask
0432                 xyz(1,CrLm),xyz(2,CrLm),xyz(3,CrLm),0);
0433         CrLm(CrLm) = Y(i,CrLm) > xM.TH(i,1);    %-Threshold (& NaN) mask
0434         if xM.I & ~YNaNrep & xM.TH(i,1)<0    %-Use implicit 0 mask
0435             CrLm(CrLm) = abs(Y(i,CrLm))>eps;
0436         end
0437     end
0438 
0439     %-Mask out voxels where data is constant
0440     %---------------------------------------------------------------
0441     CrLm(CrLm) = any(diff(Y(:,CrLm),1));
0442 
0443     %-Apply mask
0444     %---------------------------------------------------------------
0445     Y          = Y(:,CrLm);            %-Data matrix within mask
0446     CrS        = sum(CrLm);            %-#current voxels
0447     CrVox.ofs  = I(1) - 1;            %-Current voxels line offset
0448     CrVox.ind  = find(CrLm);        %-Voxel indicies (within bunch)
0449 
0450 
0451     %-if any voxels
0452     %---------------------------------------------------------------
0453     nVox  = sum(CrLm);
0454     if nVox
0455 
0456     %-Proceed with General Linear Model & smoothness estimation
0457     %===============================================================
0458     if nVar == 1                % univariate
0459 
0460         %-Estimate intrinsic correlation structure AR(1) model
0461         %-------------------------------------------------------
0462         switch xX.xVi.Form
0463 
0464             case 'AR(1)'
0465             %---------------------------------------------------
0466             fprintf('%s%30s',bs30,...
0467                     '...AR(1) estimation')         %-#
0468 
0469             for i = 1:length(xX.xVi.row)
0470             y = spm_detrend(Y(xX.xVi.row{i},:));
0471             q = 1:(size(y,1) - 1);
0472             a = sum(y(q,:).*y(q + 1,:))./sum(y(q,:).*y(q,:));
0473             A = A + [1; -mean(a)];
0474             end
0475         end
0476 
0477         %-Temporal smoothing
0478         %-------------------------------------------------------
0479         fprintf('%s%30s',bs30,...
0480                     '...temporal smoothing')     %-#
0481 
0482         KY        = mardo_99('spm_filter', 'apply',xX.K, Y);
0483 
0484         %-General linear model: least squares estimation
0485         % (Using pinv to allow for non-unique designs            )
0486         % (Including temporal convolution of design matrix with K)
0487         %-------------------------------------------------------
0488         fprintf('%s%30s',bs30,...
0489                     '...parameter estimation')   %-#
0490         beta      = xX.pKX * KY;        %-Parameter estimates
0491         res       = spm_sp('r',xX.xKXs,KY);    %-Residuals
0492         ResSS     = sum(res.^2);        %-Res sun-of-squares
0493         clear KY                %-Clear to save memory
0494 
0495     %-ManCova (assuming no filtering)
0496     %===============================================================
0497     else
0498         
0499         %-get nVar-variate response variable
0500         %-------------------------------------------------------
0501         fprintf('%s%30s',bs30,...
0502                   '...Canonical Variates Analysis')   %-#
0503 
0504         y     = zeros(nScan,nVar,nVox);
0505         Y     = zeros(nScan,nVox);
0506         res   = zeros(nScan,nVox);
0507         beta  = zeros(nBeta,nVox);
0508         V     = zeros(nVar,nVox);
0509         for i = 1:nScan
0510             for j = 1:nVar
0511                 y(i,j,:)  = spm_sample_vol(VY(i,j),...
0512                 xyz(1,CrLm),xyz(2,CrLm),xyz(3,CrLm),0);
0513             end
0514         end
0515         for i = 1:nVox
0516             
0517             %-parameter estimates
0518             %-----------------------------------------------
0519             y(:,:,i)  = y(:,:,i) - KX0*(pKX0 * y(:,:,i));
0520             BETA      = xX.pKX * y(:,:,i);
0521             h         = xX.xKXs.X * BETA;
0522             r         = y(:,:,i) - h;
0523 
0524             %-Canonical variate analysis
0525             %-----------------------------------------------
0526             [u v]     = eig(h'*h,r'*r);
0527             v         = diag(v);
0528             CU        = u(:,find(v == max(v)));
0529             
0530             %-project onto first CV
0531             %-----------------------------------------------
0532             V(:,i)    = v;
0533             res(:,i)  = r*CU;
0534             Y(:,i)    = y(:,:,i)*CU;
0535             beta(:,i) = BETA*CU;
0536 
0537         end
0538 
0539         % conpute multivariate F transform for Wilk's Lambda
0540         %-------------------------------------------------------
0541         ResSS = sum(res.^2);        %-along canonical vector
0542         W     = prod(1./(1 + V)).^(1/b);
0543         mvF   = (1 - W)./W*erdf/eidf;
0544         CrmvF = [CrmvF,mvF];
0545 
0546     end
0547 
0548     clear Y                    %-Clear to save memory
0549 
0550 
0551     %-Save betas etc. for current plane in memory as we go along
0552     % (if there is not enough memory, could save directly as float)
0553     % (analyze images, or *.mad and retreive when plane complete )
0554     %---------------------------------------------------------------
0555     CrBl     = [CrBl,beta];
0556     if params.write_res, CrResI  = [CrResI, res]; end
0557     CrResSS = [CrResSS,ResSS];
0558 
0559     end % (nVox)
0560 
0561     %-Append new inmask voxel locations and volumes
0562     %---------------------------------------------------------------
0563     XYZ            = [XYZ,xyz(:,CrLm)];    %-InMask XYZ voxel coordinates
0564     S              = S + CrS;        %-Volume analysed (voxels)
0565                             % (equals size(XYZ,2))
0566 
0567     %-Roll...
0568     %---------------------------------------------------------------
0569     BePm(I)        = CrLm;            %-"below plane" mask
0570 
0571     end % (for bch = 1:nbch)
0572 
0573 
0574     %-Plane complete, write out plane data to image files
0575     %===================================================================
0576     fprintf('%s%30s',bs30,'...saving plane')     %-#
0577 
0578     %-Mask image
0579     %-BePm now contains a complete plane mask
0580     %-------------------------------------------------------------------
0581     VM    = spm_write_plane(VM, reshape(BePm,xdim,ydim), z);
0582 
0583     %-Construct voxel indices for BePm
0584     %-------------------------------------------------------------------
0585     Q     = find(BePm);
0586     tmp   = NaN*ones(xdim,ydim);
0587 
0588     %-Write beta images
0589     %-------------------------------------------------------------------
0590     for i = 1:nBeta
0591         if length(Q), tmp(Q) = CrBl(i,:); end
0592     Vbeta(i) = spm_write_plane(Vbeta(i),tmp,z);
0593     end
0594 
0595     %-Write variance images
0596     %-------------------------------------------------------------------
0597     if params.write_res
0598       for i = 1:nScan
0599         if length(Q), tmp(Q) = CrResI(i,:); end
0600     VResI(i) = spm_write_plane(VResI(i), tmp, z);
0601       end
0602     end
0603     
0604     %-Write ResSS into ResMS (variance) image
0605     % (Scaling of ResSS to ResMS by trRV is accomplished by adjusting the
0606     % (scalefactor at the end, once the intrinsic temporal autocorrelation
0607     % (Vi has been estimated.
0608     %-------------------------------------------------------------------
0609     if length(Q), tmp(Q) = CrResSS; end
0610     VResMS = spm_write_plane(VResMS,tmp,z);        
0611 
0612     %-Write SPM (multivariate inference only)
0613     %-------------------------------------------------------------------
0614     if nVar > 1
0615     if length(Q), tmp(Q) = CrmvF; end
0616     Vspm = spm_write_plane(Vspm,tmp,z);    
0617     end
0618 
0619     %-Report progress
0620     %-------------------------------------------------------------------
0621     fprintf('%s%30s',bs30,'...done')             %-#
0622     spm_progress_bar('Set',100*(bch + nbch*(z-1))/(nbch*zdim));
0623 
0624 end % (for z = 1:zdim)
0625 fprintf('\n')                                                        %-#
0626 
0627 
0628 %=======================================================================
0629 % - P O S T   E S T I M A T I O N   C L E A N U P
0630 %=======================================================================
0631 if S == 0, warning('No inmask voxels - empty analysis!'), end
0632 
0633 
0634 %-Intrinsic autocorrelations: Vi
0635 %=======================================================================
0636 fprintf('%-40s: %30s','Design parameters','...intrinsic autocorrelation') %-#
0637 
0638 %-Compute (session specific) intrinsic autocorrelation Vi
0639 %-----------------------------------------------------------------------
0640 switch xX.xVi.Form
0641     case 'AR(1)'
0642     %---------------------------------------------------
0643     p     = length(A) - 1;            % order AR(p)
0644     A     = A/A(1);
0645     for i = 1:length(xX.xVi.row)
0646         q     = xX.xVi.row{i};
0647         n     = length(q);
0648         Ki    = inv(spdiags(ones(n,1)*A',-[0:p],n,n));
0649         Ki    = Ki.*(Ki > 1e-6);
0650         Vi    = Ki*Ki';
0651         D     = spdiags(sqrt(1./diag(Vi)),0,n,n);
0652         Vi    = D*Vi*D;
0653         xX.xVi.Vi(q,q) = Vi;
0654     end
0655 end
0656 xX.xVi.Param = A;
0657 
0658 
0659 %-[Re]-enter Vi & derived values into design structure xX
0660 %-----------------------------------------------------------------------
0661 fprintf('%s%30s',bs30,'...V, & traces')          %-#
0662 
0663 KVi      = mardo_99('spm_filter', 'apply',xX.K, xX.xVi.Vi);
0664 xX.V     = mardo_99('spm_filter', 'apply',xX.K,KVi');     %-V matrix
0665 xX.pKXV  = xX.pKX*xX.V;                %-for contrast variance weight
0666 xX.Bcov  = xX.pKXV*xX.pKX';            %-Variance of est. param.
0667 [xX.trRV,xX.trRVRV] ...                %-Variance expectations
0668          = spm_SpUtil('trRV',xX.xKXs,xX.V);
0669 xX.erdf  = xX.trRV^2/xX.trRVRV;            %-Effective residual d.f.
0670 
0671 %-Compute scaled design matrix for display purposes
0672 %-----------------------------------------------------------------------
0673 fprintf('%s%30s',bs30,'...scaling DesMtx')       %-#
0674 xX.nKX   = spm_DesMtx('sca',xX.xKXs.X,xX.Xnames);
0675 
0676 %-Set VResMS scalefactor as 1/trRV (raw voxel data is ResSS)
0677 % Note that, due to a bug in the SPM2 vol utils, this scalefactor does
0678 % not get properly written if we just use pinfo, so we have to save,
0679 % reload to set it here.
0680 %-----------------------------------------------------------------------
0681 VResMS = spm_close_vol(VResMS);
0682 img    = spm_read_vols(VResMS);
0683 img    = img / xX.trRV;
0684 VResMS = spm_write_vol(VResMS, img);
0685 
0686 %-"close" written image files, updating scalefactor information
0687 %=======================================================================
0688 fprintf('%s%30s',bs30,'...closing image files')  %-#
0689 VM                      = spm_close_vol(VM);
0690 Vbeta                   = spm_close_vol(Vbeta);
0691 if nVar > 1,   Vspm     = spm_close_vol(Vspm);     end
0692 if params.write_res, VResI = spm_close_vol(VResI); end
0693 
0694 %-Unmap files, retaining image names, and reset erdf if MV
0695 %-----------------------------------------------------------------------
0696 fprintf('%s%30s',bs30,'...tidying file handles') %-#
0697 VM     = VM.fname;
0698 Vbeta  = {Vbeta.fname}';
0699 VResMS = VResMS.fname;
0700 if nVar > 1
0701     xCon.Vspm = Vspm.fname;
0702     xX.erdf   = erdf;
0703 end
0704 if params.write_res, VResI = {{VResI.fname}'}; end
0705 
0706 fprintf('%s%30s\n',bs30,'...done')               %-#
0707 
0708 
0709 
0710 %-Save remaining results files and analysis parameters
0711 %=======================================================================
0712 fprintf('%-40s: %30s','Saving results','...writing')                 %-#
0713 
0714 %-Save analysis parameters in phiw_spm.mat file
0715 %-----------------------------------------------------------------------
0716 SPM = mars_struct('ffillmerge', SPM, ...
0717           struct(...
0718               'SPMid', SPMid, ...
0719               'VY',    {VY}, ...
0720               'xX',    xX, ...
0721               'xM',    xM, ...
0722               'M',     M, ...
0723               'DIM',   DIM, ...
0724               'VM',    VM, ...
0725               'Vbeta', {Vbeta}, ...
0726               'VResI', VResI, ...
0727               'VResMS',VResMS, ...
0728               'XYZ',   XYZ, ...
0729               'F_iX0', F_iX0, ...
0730               'S',     S, ...
0731               'xPhi',  struct('estimated', 1, 'wave', thin(wvobj)), ...
0732               'swd',   pwd, ...
0733               'fname', phiw_mat_name, ...
0734               'xCon',  xCon));
0735 
0736 % save design
0737 savestruct(phiw_mat_name, SPM);
0738 
0739 fprintf('%s%30s\n',bs30,'...done')               %-#
0740 
0741 %=======================================================================
0742 %- E N D: Cleanup GUI
0743 %=======================================================================
0744 spm_progress_bar('Clear')
0745 spm('FigName','Stats: done',Finter); spm('Pointer','Arrow')
0746 fprintf('%-40s: %30s\n','Completed',spm('time'))                     %-#
0747 fprintf('...use the results section for assessment\n\n')             %-#
pr_estimate

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
