1 function [SPM,xSPM] = csl_getSPM2(xSPM)
   2 % Adapted version of spm_getSPM.m with few changes to handle the xSPM structure.
   3 % allowing automatic processing.      
   4 % computes a specified and thresholded SPM/PPM following parameter estimation
   5 % FORMAT [SPM,xSPM] = csl_getSPM2;
   6 %
   7 % xSPM      - structure containing SPM, distribution & filtering details
   8 % .swd      - SPM working directory - directory containing current SPM.mat
   9 % .title    - title for comparison (string)
  10 % .Z        - minimum of n Statistics {filtered on u and k}
  11 % .n        - number of conjoint tests        
  12 % .STAT     - distribution {Z, T, X, F or P}     
  13 % .df       - degrees of freedom [df{interest}, df{residual}]
  14 % .STATstr  - description string     
  15 % .Ic       - indices of contrasts (in SPM.xCon)
  16 % .Im       - indices of masking contrasts (in xCon)
  17 % .pm       - p-value for masking (uncorrected)
  18 % .Ex       - flag for exclusive or inclusive masking
  19 % .u        - height threshold
  20 % .k        - extent threshold {voxels}
  21 % .XYZ      - location of voxels {voxel coords}
  22 % .XYZmm    - location of voxels {mm}
  23 % .S        - search Volume {voxels}
  24 % .R        - search Volume {resels}
  25 % .FWHM     - smoothness {voxels}     
  26 % .M        - voxels -> mm matrix
  27 % .iM       - mm -> voxels matrix
  28 % .VOX      - voxel dimensions {mm} - column vector
  29 % .DIM      - image dimensions {voxels} - column vector
  30 % .Vspm     - Mapped statistic image(s)
  31 % .Ps       - list of P values for voxels at SPM.xVol.XYZ (used by FDR)
  32 %
  33 % Required feilds of SPM
  34 %
  35 % xVol   - structure containing details of volume analysed
  36 %
  37 % xX     - Design Matrix structure
  38 %        - (see spm_spm.m for structure)
  39 %
  40 % xCon   - Contrast definitions structure array
  41 %        - (see also spm_FcUtil.m for structure, rules & handling)
  42 % .name  - Contrast name
  43 % .STAT  - Statistic indicator character ('T', 'F' or 'P')
  44 % .c     - Contrast weights (column vector contrasts)
  45 % .X0    - Reduced design matrix data (spans design space under Ho)
  46 %          Stored as coordinates in the orthogonal basis of xX.X from spm_sp
  47 %          (Matrix in SPM99b)  Extract using X0 = spm_FcUtil('X0',...
  48 % .iX0   - Indicates how contrast was specified:
  49 %          If by columns for reduced design matrix then iX0 contains the
  50 %          column indices. Otherwise, it's a string containing the
  51 %          spm_FcUtil 'Set' action: Usuall one of {'c','c+','X0'}
  52 % .X1o   - Remaining design space data (X1o is orthogonal to X0)
  53 %          Stored as coordinates in the orthogonal basis of xX.X from spm_sp
  54 %          (Matrix in SPM99b)  Extract using X1o = spm_FcUtil('X1o',...
  55 % .eidf  - Effective interest degrees of freedom (numerator df)
  56 %        - Or effect-size threshold for Posterior probability
  57 % .Vcon  - Name of contrast (for 'T's) or ESS (for 'F's) image
  58 % .Vspm  - Name of SPM image
  59 %
  60 % In addition, the xCon.mat file is updated. For newly evaluated
  61 % contrasts, SPM images (spmT_????.{img,hdr}) are written, along with
  62 % contrast (con_????.{img,hdr}) images for SPM{T}'s, or Extra
  63 % Sum-of-Squares images (ess_????.{img,hdr}) for SPM{F}'s.
  64 % 
  65 % The contrast images are the weighted sum of the parameter images,
  66 % where the weights are the contrast weights, and are uniquely
  67 % estimable since contrasts are checked for estimability by the
  68 % contrast manager. These contrast images (for appropriate contrasts)
  69 % are suitable summary images of an effect at this level, and can be
  70 % used as input at a higher level when effecting a random effects
  71 % analysis. (Note that the ess_????.{img,hdr} and
  72 % SPM{T,F}_????.{img,hdr} images are not suitable input for a higher
  73 % level analysis.) See spm_RandFX.man for further details.
  74 %
  75 %_______________________________________________________________________
  76 %
  77 % spm_getSPM prompts for an SPM and applies thresholds {u & k}
  78 % to a point list of voxel values (specified with their locations {XYZ})
  79 % This allows the SPM be displayed and characterized in terms of regionally 
  80 % significant effects by subsequent routines.
  81 % 
  82 % For general linear model Y = XB + E with data Y, desgin matrix X,
  83 % parameter vector B, and (independent) errors E, a contrast c'B of the
  84 % parameters (with contrast weights c) is estimated by c'b, where b are
  85 % the parameter estimates given by b=pinv(X)*Y.
  86 % 
  87 % Either single contrasts can be examined or conjunctions of different
  88 % contrasts. Contrasts are estimable linear combinations of the
  89 % parameters, and are specified using the SPM contrast manager
  90 % interface [spm_conman.m]. SPMs are generated for the null hypotheses
  91 % that the contrast is zero (or zero vector in the case of
  92 % F-contrasts). See the help for the contrast manager [spm_conman.m]
  93 % for a further details on contrasts and contrast specification.
  94 % 
  95 % A conjunction assesses the conjoint expression of two or more
  96 % effects. The conjunction SPM is the minimum of the component SPMs
  97 % defined by the multiple contrasts. The distributional results used
  98 % for minimum fileds require the SPMs to be identically distributed and
  99 % independent. Thus, all component SPMs must be either SPM{t}'s, or
 100 % SPM{F}'s with the same degrees of freedom. Independence is roughly
 101 % guaranteed for large degrees of freedom (and independent data) by
 102 % ensuring that the contrasts are "orthogonal". Note that it is *not*
 103 % the contrast weight vectors per se that are required to be
 104 % orthogonal, but the subspaces of the data space implied by the null
 105 % hypotheses defined by the contrasts (c'pinv(X)). Furthermore,
 106 % this assumes that the errors are i.i.d. (i.e. the estimates are
 107 % maximum likelihood or Gauss-Markov. This is the default in spm_spm). 
 108 % 
 109 % To ensure approximate independence of the component SPMs in a
 110 % conjunction, non-orthogonal contrasts are serially orthogonalised
 111 % in the order specified, possibly generating new contrasts, such that the
 112 % second is orthogonal to the first, the third to the first two, and so on.
 113 %
 114 % Masking simply eliminates voxels from the current contrast if they
 115 % do not survive an uncorrected p value (based on height) in one or
 116 % more further contrasts.  No account is taken of this masking in the
 117 % statistical inference pertaining to the masked contrast.
 118 % 
 119 % The SPM is subject to thresholding on the basis of height (u) and the
 120 % number of voxels comprising its clusters {k}. The height threshold is
 121 % specified as above in terms of an [un]corrected p value or
 122 % statistic.  Clusters can also be thresholded on the basis of their
 123 % spatial extent. If you want to see all voxels simply enter 0.  In this
 124 % instance the 'set-level' inference can be considered an 'omnibus test'
 125 % based on the number of clusters that obtain.
 126 %
 127 % BAYESIAN INFERENCE AND PPMS - POSTERIOR PROBABILITY MAPS
 128 %
 129 % If conditional estimates are available (and your contrast is a T
 130 % contrast) then you are asked whether the inference should be 'Bayesian'
 131 % or 'classical' (using GRF).  If you choose Bayesian the contrasts are of
 132 % conditional (i.e. MAP) estimators and the inference image is a
 133 % posterior probability map (PPM).  PPMs encode the probability that the
 134 % contrast exceeds a specified threshold.  This threshold is stored in
 135 % the xCon.eidf.  Subsequent plotting and tables will use the conditional
 136 % estimates and associated posterior or conditional probabilities.
 137 % 
 138 % see spm_results_ui.m for further details of the SPM results section.
 139 % see also spm_contrasts.m
 140 %_______________________________________________________________________
 141 % @(#)spm_getSPM.m	2.51 Andrew Holmes, Karl Friston & Jean-Baptiste Poline 03/05/22
 142 %- KHERIF Ferath csl, Cambridge, UK.
 143 
 144 SCCSid = '2.51';
 145 
 146 %-GUI setup
 147 %-----------------------------------------------------------------------
 148 SPMid  = spm('SFnBanner',mfilename,SCCSid);
 149 spm_help('!ContextHelp',mfilename)
 150 
 151 %-Select SPM.mat & note SPM results directory
 152 %-----------------------------------------------------------------------
 153 %swd    = spm_str_manip(spm_get(1,'SPM.mat','Select SPM.mat'),'H');
 154 
 155 try
 156 	swd = spm_str_manip(xSPM.swd,'H');
 157 catch	
 158 	swd    = spm_str_manip(spm_get(1,'SPM.mat','Select SPM.mat'),'H');
 159 
 160 end
 161 %-Preliminaries...
 162 %=======================================================================
 163 
 164 %-Load SPM.mat
 165 %-----------------------------------------------------------------------
 166 load(fullfile(swd,'SPM.mat'));
 167 SPM.swd = swd;
 168 
 169 %-Get volumetric data from SPM.mat
 170 %-----------------------------------------------------------------------
 171 try
 172 	xX   = SPM.xX;				%-Design definition structure
 173 	XYZ  = SPM.xVol.XYZ;			%-XYZ coordinates
 174 	S    = SPM.xVol.S;			%-search Volume {voxels}
 175 	R    = SPM.xVol.R;			%-search Volume {resels}
 176 	M    = SPM.xVol.M(1:3,1:3);		%-voxels to mm matrix
 177 	VOX  = sqrt(diag(M'*M))';		%-voxel dimensions
 178 catch
 179 
 180 	% check the model has been estimated
 181 	%---------------------------------------------------------------
 182 	str = {	'This model has not been estimated.';...
 183 		'Would you like to estimate it now?'};
 184 	if spm_input(str,1,'bd','yes|no',[1,0],1)
 185 		[SPM] = spm_spm(SPM);
 186 	else
 187 		return
 188 	end
 189 end
 190 
 191 
 192 %-Contrast definitions
 193 %=======================================================================
 194 
 195 %-Load contrast definitions (if available)
 196 %-----------------------------------------------------------------------
 197 try
 198 	xCon = SPM.xCon;
 199 catch
 200 	xCon = {};
 201 end
 202 
 203 
 204 %=======================================================================
 205 % - C O N T R A S T S ,  S P M   C O M P U T A T I O N ,   M A S K I N G
 206 %=======================================================================
 207 
 208 %-Get contrasts
 209 %-----------------------------------------------------------------------
 210 try 
 211 	Ic=xSPM.Ic
 212 	
 213 catch 
 214 [Ic,xCon] = spm_conman(xX,xCon,'T&F',Inf,...
 215 		'	Select contrasts...',' for conjunction',1);
 216 end
 217 
 218 
 219 
 220 %-Enforce orthogonality of multiple contrasts for conjunction
 221 % (Orthogonality within subspace spanned by contrasts)
 222 %-----------------------------------------------------------------------
 223 if length(Ic) > 1 & ~spm_FcUtil('|_?',xCon(Ic), xX.xKXs)
 224 
 225     
 226     %-Successively orthogonalise
 227     %-NB: This loop is peculiarly controlled to account for the
 228     %     possibility that Ic may shrink if some contrasts diasppear
 229     %     on orthogonalisation (i.e. if there are colinearities)
 230     %-------------------------------------------------------------------
 231     i = 1; while(i < length(Ic)), i = i + 1;
 232     
 233 	%-Orthogonalise (subspace spanned by) contrast i wirit previous
 234 	%---------------------------------------------------------------
 235 	oxCon = spm_FcUtil('|_',xCon(Ic(i)), xX.xKXs, xCon(Ic(1:i-1)));
 236     
 237 	%-See if this orthogonalised contrast has already been entered
 238 	% or is colinear with a previous one. Define a new contrast if
 239 	% neither is the case.
 240 	%---------------------------------------------------------------
 241 	d     = spm_FcUtil('In',oxCon,xX.xKXs,xCon);
 242 
 243 	if spm_FcUtil('0|[]',oxCon,xX.xKXs)
 244 
 245 	    %-Contrast was colinear with a previous one - drop it
 246 	    %-----------------------------------------------------------
 247 	    Ic(i) = [];
 248 	    i     = i - 1;
 249 
 250 	elseif any(d)
 251 
 252 	    %-Contrast unchanged or already defined - note index
 253 	    %-----------------------------------------------------------
 254 	    Ic(i) = min(d);
 255 
 256 	else
 257 
 258 	    %-Define orthogonalised contrast as new contrast
 259 	    %-----------------------------------------------------------
 260 	    oxCon.name = [xCon(Ic(i)).name,' (orth. w.r.t {',...
 261     	    	sprintf('%d,',Ic(1:i-2)), sprintf('%d})',Ic(i-1))];
 262     	    xCon  = [xCon, oxCon];
 263 	    Ic(i) = length(xCon); 
 264 	end
 265 
 266     end % while...
 267 end % if length(Ic)...
 268 
 269 
 270 %-Get contrasts for masking
 271 %-----------------------------------------------------------------------
 272 maskCon= 0; %by defualt don't use masking. 
 273 try
 274 	Im = xSPM.Im;
 275 	if ~isempty(Im) maskCon=1; end 
 276 catch
 277 	maskCon=spm_input('mask with other contrast(s)','+1','y/n',[1,0],2);
 278 
 279 end
 280 
 281 if maskCon
 282 	
 283 	try
 284 		Im = xSPM.Im;
 285 	
 286 	catch
 287 	
 288 	[Im,xCon] = spm_conman(xX,xCon,'T&F',-Inf,...
 289 		'Select contrasts for masking...',' for masking',1);
 290 	end
 291 
 292 	%-Threshold for mask (uncorrected p-value)
 293 	%---------------------------------------------------------------
 294 	try
 295 		pm = xSPM.pm;
 296 	catch
 297 		pm = spm_input('uncorrected mask p-value','+1','r',0.05,1,[0,1]);
 298 	end
 299 
 300 	%-Inclusive or exclusive masking
 301 	%---------------------------------------------------------------
 302 	try
 303 		Ex = xSPM.Ex;
 304 	catch
 305 		Ex = spm_input('nature of mask','+1','b','inclusive|exclusive',[0,1]);
 306 	end
 307 else
 308 	Im = [];
 309 	pm = [];
 310 	Ex = [];
 311 end
 312 
 313 
 314 %-Create/Get title string for comparison
 315 %-----------------------------------------------------------------------
 316 if length(Ic) == 1
 317     try
 318         str  = xCon(Ic).name;
 319     catch
 320        str='RFX' % djm added 
 321     end
 322 else
 323 	str  = [sprintf('contrasts {%d',Ic(1)),sprintf(',%d',Ic(2:end)),'}'];
 324 end
 325 if Ex
 326 	mstr = 'masked [excl.] by';
 327 else
 328 	mstr = 'masked [incl.] by';
 329 end
 330 if length(Im) == 1
 331 	str  = sprintf('%s (%s %s at p=%g)',str,mstr,xCon(Im).name,pm);
 332 
 333 elseif ~isempty(Im)
 334 	str  = [sprintf('%s (%s {%d',str,mstr,Im(1)),...
 335 		sprintf(',%d',Im(2:end)),...
 336 		sprintf('} at p=%g)',pm)];
 337 end
 338 
 339 titlestr='';
 340 if isfield(xSPM,'title')
 341 titlestr     = xSPM.title;
 342 end
 343 if isempty(titlestr) 
 344 titlestr  = str;
 345 end   
 346 
 347 
 348 
 349 %-Bayesian or classical Inference?
 350 %-----------------------------------------------------------------------
 351 if isfield(SPM,'PPM') & xCon(Ic(1)).STAT == 'T'
 352 
 353     if length(Ic) == 1 & isempty(xCon(Ic).Vcon)
 354 
 355 	if spm_input('Inference',1,'b',{'Bayesian','classical'},[1 0]);
 356 
 357 	% set STAT to 'P'
 358         %---------------------------------------------------------------
 359 	xCon(Ic).STAT = 'P';
 360 
 361 	%-Get Bayesian threshold (Gamma) stored in xCon(Ic).eidf
 362 	% The default is one conditional s.d. of the contrast
 363         %---------------------------------------------------------------
 364 	str           = 'threshold {default: prior s.d.}';
 365 	Gamma         = sqrt(xCon(Ic).c'*SPM.PPM.Cb*xCon(Ic).c);
 366 	xCon(Ic).eidf = spm_input(str,'+1','e',sprintf('%0.2f',Gamma));
 367 
 368 	end
 369     end
 370 end
 371 
 372 
 373 %-Compute & store contrast parameters, contrast/ESS images, & SPM images
 374 %=======================================================================
 375 SPM.xCon = xCon;
 376 SPM      = spm_contrasts(SPM,unique([Ic,Im]));
 377 xCon     = SPM.xCon;
 378 VspmSv   = cat(1,xCon(Ic).Vspm);
 379 STAT     = xCon(Ic(1)).STAT;
 380 n        = length(Ic);
 381 
 382 %-Check conjunctions - Must be same STAT w/ same df
 383 %-----------------------------------------------------------------------
 384 if (n > 1) & (any(diff(double(cat(1,xCon(Ic).STAT)))) | ...
 385 	      any(abs(diff(cat(1,xCon(Ic).eidf))) > 1))
 386 	error('illegal conjunction: can only conjoin SPMs of same STAT & df')
 387 end
 388 
 389 
 390 %-Degrees of Freedom and STAT string describing marginal distribution
 391 %-----------------------------------------------------------------------
 392 df          = [xCon(Ic(1)).eidf xX.erdf];
 393 if n > 1
 394 	str = sprintf('^{%d}',n);
 395 else
 396 	str = '';
 397 end
 398 
 399 switch STAT
 400 case 'T'
 401 	STATstr = sprintf('%c%s_{%.0f}','T',str,df(2));
 402 case 'F'
 403 	STATstr = sprintf('%c%s_{%.0f,%.0f}','F',str,df(1),df(2));
 404 case 'P'
 405 	STATstr = sprintf('%s^{%0.2f}','PPM',df(1));
 406 end
 407 
 408 
 409 %-Compute (unfiltered) SPM pointlist for masked conjunction requested
 410 %=======================================================================
 411 fprintf('\t%-32s: %30s\n','SPM computation','...initialising')         %-#
 412 
 413 
 414 %-Compute conjunction as minimum of SPMs
 415 %-----------------------------------------------------------------------
 416 Z         = Inf;
 417 for i     = Ic
 418 	Z = min(Z,spm_get_data(xCon(i).Vspm,XYZ));
 419 end
 420 
 421 % P values for False Discovery FDR rate computation (all search voxels)
 422 %=======================================================================
 423 switch STAT
 424 case 'T'
 425 	Ps = (1 - spm_Tcdf(Z,df(2))).^n;
 426 case 'P'
 427 	Ps = (1 - Z).^n;
 428 case 'F'
 429 	Ps = (1 - spm_Fcdf(Z,df)).^n;
 430 end
 431 
 432 
 433 %-Compute mask and eliminate masked voxels
 434 %-----------------------------------------------------------------------
 435 for i = Im
 436 	fprintf('%s%30s',sprintf('\b')*ones(1,30),'...masking')
 437 
 438 	Mask = spm_get_data(xCon(i).Vspm,XYZ);
 439 	um   = spm_u(pm,[xCon(i).eidf,xX.erdf],xCon(i).STAT);
 440 	if Ex
 441 		Q = Mask <= um;
 442 	else
 443 		Q = Mask >  um;
 444 	end
 445 	XYZ       = XYZ(:,Q);
 446 	Z         = Z(Q);
 447 	if isempty(Q)
 448 		fprintf('\n')                                           %-#
 449     		warning(sprintf('No voxels survive masking at p=%4.2f',pm))
 450 		break
 451 	end
 452 end
 453 
 454 %-clean up interface
 455 %-----------------------------------------------------------------------
 456 fprintf('\t%-32s: %30s\n','SPM computation','...done')         %-#
 457 spm('Pointer','Arrow')
 458 
 459 
 460 
 461 %=======================================================================
 462 % - H E I G H T   &   E X T E N T   T H R E S H O L D S
 463 %=======================================================================
 464 
 465 %-Height threshold - classical inference
 466 %-----------------------------------------------------------------------
 467 u      = -Inf;
 468 k      = 0;
 469 if STAT ~= 'P'
 470 
 471 
 472     %-Get height threshold
 473     %-------------------------------------------------------------------
 474     str = 'FWE|FDR|none';
 475     % str = 'FWE|none';      % Use this line to disable FDR threshold
 476    try
 477 	 Mcp = xSPM.Mcp;
 478    catch  
 479 	 
 480 	Mcp = spm_input('p value adjustment to control','+1','b',str,[],1);
 481     end
 482     switch Mcp
 483 
 484 
 485 	case 'FWE' % family-wise false positive rate
 486         %---------------------------------------------------------------
 487 	try 
 488 		u = xSPM.u;
 489 	catch
 490 		u  = spm_input('p value (family-wise error)','+0','r',0.05,1,[0,1]);
 491 	end 
 492 	u  = spm_uc(u,df,STAT,R,n,S);
 493 
 494 	case 'FDR' % False discovery rate
 495 	%---------------------------------------------------------------	
 496 	try 
 497 		u = xSPM.u;
 498 	catch
 499 		u  = spm_input('p value (false discovery rate)','+0','r',0.05,1,[0,1]);
 500 	end
 501 	u  = spm_uc_FDR(u,df,STAT,n,VspmSv,0);
 502 
 503 	otherwise  %-NB: no adjustment
 504 	% p for conjunctions is p of the conjunction SPM
 505    	%---------------------------------------------------------------
 506 	try 
 507 		u = xSPM.u;
 508 	catch
 509 		u  = spm_input(['threshold {',STAT,' or p value}'],'+0','r',0.001,1);
 510 	end	
 511 	if u <= 1; u = spm_u(u^(1/n),df,STAT); end
 512 
 513     end
 514 
 515 %-Height threshold - Bayesian inference
 516 %-----------------------------------------------------------------------
 517 elseif STAT == 'P'
 518 	try 
 519 	  u = xSPM.u;
 520 	catch
 521 	  u  = spm_input(['p value threshold for PPM'],'+0','r',.95,1);
 522 	end
 523 
 524 end % (if STAT)
 525 
 526 %-Calculate height threshold filtering
 527 %-------------------------------------------------------------------
 528 Q      = find(Z > u);
 529 
 530 %-Apply height threshold
 531 %-------------------------------------------------------------------
 532 Z      = Z(:,Q);
 533 XYZ    = XYZ(:,Q);
 534 if isempty(Q)
 535 	warning(sprintf('No voxels survive height threshold u=%0.2g',u))
 536 end
 537 
 538 
 539 %-Extent threshold (disallowed for conjunctions)
 540 %-----------------------------------------------------------------------
 541 if ~isempty(XYZ) & length(Ic) == 1
 542 
 543     %-Get extent threshold [default = 0]
 544     %-------------------------------------------------------------------
 545     try 
 546       k  = xSPM.k;
 547     catch
 548       k     = spm_input('& extent threshold {voxels}','+1','r',0,1,[0,Inf]);
 549     end
 550 
 551     %-Calculate extent threshold filtering
 552     %-------------------------------------------------------------------
 553     A     = spm_clusters(XYZ);
 554     Q     = [];
 555     for i = 1:max(A)
 556         j = find(A == i);
 557         if length(j) >= k; Q = [Q j]; end
 558     end
 559 
 560     % ...eliminate voxels
 561     %-------------------------------------------------------------------
 562     Z     = Z(:,Q);
 563     XYZ   = XYZ(:,Q);
 564     if isempty(Q)
 565 	warning(sprintf('No voxels survive extent threshold k=%0.2g',k))
 566     end
 567 
 568 else
 569 
 570     k = 0;
 571 
 572 end % (if ~isempty(XYZ))
 573 
 574 
 575 %=======================================================================
 576 % - E N D
 577 %=======================================================================
 578 fprintf('\t%-32s: %30s\n','SPM computation','...done')         %-#
 579 
 580 %-Assemble output structures of unfiltered data
 581 %=======================================================================
 582 xSPM   = struct('swd',		swd,...
 583 		'title',	titlestr,...
 584 		'Z',		Z,...
 585 		'n',		n,...
 586 		'STAT',		STAT,...
 587 		'df',		df,...
 588 		'STATstr',	STATstr,...
 589 		'Ic',		Ic,...
 590 		'Im',		Im,...
 591 		'pm',		pm,...
 592 		'Ex',		Ex,...
 593 		'u',		u,...
 594 		'k',		k,...
 595 		'XYZ',		XYZ,...
 596 		'XYZmm',	SPM.xVol.M(1:3,:)*[XYZ; ones(1,size(XYZ,2))],...
 597 		'S',		SPM.xVol.S,...
 598 		'R',		SPM.xVol.R,...
 599 		'FWHM',		SPM.xVol.FWHM,...
 600 		'M',		SPM.xVol.M,...
 601 		'iM',		SPM.xVol.iM,...
 602 		'DIM',		SPM.xVol.DIM,...
 603 		'VOX',		VOX,...
 604 		'Vspm',		VspmSv,...
 605 		'Ps',		Ps, ...
 606 		'Mcp',		Mcp);
 607 
 608 % RESELS per voxel (density) if it exists
 609 %-----------------------------------------------------------------------
 610 if isfield(SPM,'VRpv'), xSPM.VRpv = SPM.VRpv; end