function snpm_cp(CWD)
% Nonparametric Perm/Rando statistical analysis with General linear model
% FORMAT snpm_cp(CWD)
%
% CWD -	Directory containing SnPMcfg.mat configuration file
%_______________________________________________________________________
% 
% snpm_cp is the engine of the SnPM toolbox and implements the general
% linear model for a set of design matrices, each design matrix 
% constituting one permutation.  First the "correct" permutation
% is calculated in its entirety, then all subsequent permutations are 
% calculated, possibly on a plane-by-plane basis.  
%
% The output of snpm_cp parallels spm_spm: for the correct permutation
% image files containing parameter estimates, statistic values, and F
% values are saved (this is in distinction from SnPM2 and previous
% versions, where this information was saved in .mat files); the permutation 
% distribution of the statistic interest and (optionally) suprathreshold
% stats are also saved.  All results are written to the directory
% that CfgFile resides in.  IMPORTANT: Existing results are overwritten
% without prompting.
%
% Unlike spm_spm, voxels are not discarded on the basis of the F statistic.
% All gray matter voxels (as defined by the gray matter threshold) are
% retained for analysis; note that this will increase the size of all .mat
% files.
%
%
%-----------------------------------------------------------------------
%
% Output File Descriptions:
%
%   XYZ.mat contains a 3 x N matrix of the x,y and z location of the
% voxels in SPMF in mm (usually referring the the standard anatomical
% space (Talairach and Tournoux 1988)} (0,0,0) corresponds to the
% centre of the voxel specified by ORIGIN in the *.hdr of the original
% and related data.
%
%   BETA.mat contains a p x S matrix of the p parameter estimates at
% each of the S voxels for the correct permutation.  These parameters 
% include all effects specified by the design matrix.
%
%   SnPMt.mat contains a 1 x S matrix of the statistic of interest (either
% t or pseudo-t if variance smoothing is used) supplied for all S voxels at
% locations XYZ.
%
%   SnPMucp.mat contains a 1 x S matrix of the nonparametric P values of
% the statistic of interest supplied for all S voxels at locations XYZ.
%
%   SnPM.mat contains a collection of strings and matrices that pertain 
% to the analysis.  In contrast to spm_spm's SPM.mat, most of the essential
% matrices are in the any of the matrices stored here in the CfgFile
% and hence are not duplicated here.   Included are the number of voxels
% analyzed (S) and the image and voxel dimensions [V].  See below
% for complete listing.
%
% snpm_cp writes out the following image files (for each image, there are
% two files: .img and .hdr files)  
%  
% beta_**** (from 0001 to p): p images of p parameter estimates at each
% voxel for the correct permutation. These p parameters include all
% effects specified by the design matrix. 
%
% ResMS: One image of residual mean square errors at each voxel. 
% 
% (SnPM, like SPM, only implements single tailed tests. In the following
% files, '+' or '-' correspond to 'positive' or 'negative' effects (as in
% snpm_pp.m). Here, '+' images are the images for large values,
% indicating evidence against the null hypothesis in favour of a positive
% alternative (activation, or positive slope in a covariate analysis))
% 
% snpmT+ & snpmT-: Images of the statistic of interest (either t or
% pseduo-t if variance smoothing is used), positive or negative. 
% The numbers (i.e. not NaN) saved in snpmT+ images are also saved in the
% SnPMt.mat file. 
% 
% lP+ & lP-: Images of -log10(uncorrected non-parametric P-values,
% positive or negative).
% 
% lP_FWE+ & lP_FWE-: Images of -log10(FWE-corrected non-parametric
% P-values, positive or negative). Here, FWE-corrected non-parametric
% P-values are the proportion of the permutation distribution for the
% maximal statistic which exceeds the statistic image at the voxel. 
%
% lP_FDR+ & lP_FDR-: Images of -log10(FDR-corrected non-parametric
% P-values, positive or negative). 
%
% The following is an example of matlab codes for reading in an image file. 
% P='.../.../beta_0001.img';
% V=spm_vol(P);
% Y=spm_read_vols(V);
% Y(~isfinite(Y))=[]; %delete NaN values from vector Y.
%
%
%-----------------------------------------------------------------------
%
% As an "engine", snpm_cp does not produce any graphics; if the SPM windows
% are open, a progress thermometer bar will be displayed.  
%
% If out-of-memory problems are encountered, the first line of defense is to
% run snpm_cp in a virgin matlab session with out first starting SPM.
%
%
% Variables saved in SnPM.mat
%=======================================================================
%
% S              Volume analyzed (in voxels)
% V              Volume handles (see spm_vol)
% df             Residual degrees of freedom of raw t-statistic
% MaxT           nPerm x 2 matrix of [max;min] t-statistics per perm
% ST_Ut          Threshold above which suprathreshold info was collected.
%                Voxel locations, t and perm are saved in SnPM_ST.mat for
%                t's greater than ST_Ut. ST_Ut=Inf if not saving STCdata
%
% s_SnPM_save    List of variables saved in SnPM.mat file
% CfgFile        SnPM config sile used (full pathname)
% s_SnPMcfg_save List of variables saved in SnPMcfg.mat file
% 
% Data structure of SnPM_ST.mat: suprathreshold stats (if collected)
%-----------------------------------------------------------------------
% 5xn matrix, each column containing:
%       [x, y, z, abs(T), perm]'
%       perm is negative if T was negative
%
%_______________________________________________________________________
% % @(#)snpm_cp_0811.m	3.38 	3.5 Thomas Nichols, Andrew Holmes 04/06/15 05/09/26
%	$Id: snpm_cp.m,v 1.2 2007/06/23 16:44:00 huizhang Exp $	



%-----------------------------functions-called------------------------
% spm_append_96
% spm_conv
% spm_figure
% spm_select
% spm_hread
% spm_invTcdf
% spm_matrix
% spm_progress_bar
% spm_smooth
% spm_str_manip
%-----------------------------functions-called------------------------


% Programmers / Hackers help / Code notes...
%=======================================================================
% snpm_cp is modeled after SPM95's spm_spm, which is different from the
% current spm_spm (which is similar to SPM99).  While the current
% (SPM99-SPM5) spm_spm reads data in 'planks' (a partial or whole slice),
% snpm_cp reads either a plane at a time or reads the entire dataset into
% memory ('Volumetric mode', bVolm=1).
%
% If bVolm is true, this function will load the entire dataset (all
% planes, all subjects) into memory.  If bVolm is false the dataset
% will be loaded a plane at a time (all subjects), but smoothing the
% variance in the z-direction is not permitted.  A version that 
% does not load the whole volume but allows volumetric smoothing
% is under development.
% 
% For designs with large number of permutations AND bVolm true, a possible
% approach would be have a stopping feature where the user had decided
% "enough" permutations had run.  Not sure if this is useful.

%-Variable "decoder"
%-----------------------------------------------------------------------
% bWin    - Do we have windows?
% bVarSm  - Variance Smoothing?
% bVolm   - Work on whole volume at once?
% q       - number of observations
% p       - number of predictors
% r       - Model degrees of freedom
% df      - Residual degrees of freedom
% nPerm   - number of permutations
% WorkDim - Number of voxels read in (either a plane's worth or the whole lot)
% MaxT    - Permutation Distribution of intensity maximum
% nP      - WorkDim vector of nonparametric P-values
%
%-Supratreshold Threshold
% 
% STalpha - if parametric T, critical val for STalpha ised
% STprop  - if pseudo T, 100*(1-STprop)%ile of correct perm's values used


%-Setup
%=======================================================================
global defaults
if isempty(defaults), spm_defaults; end


fprintf('\nSnPM: snpm_cp\n'),fprintf('%c','='*ones(1,72)),fprintf('\n')
disp('Initialising...');
bWin = ~isempty(spm_figure('FindWin','Interactive'));
s_SnPM_save = 's_SnPM_save CfgFile s_SnPMcfg_save S V df1 df MaxT ST_Ut STAT';

%-Check arguments & parameters from CfgFile
%-----------------------------------------------------------------------
if nargin == 0
	tmp = spm_select(1,'SnPMcfg.mat','Select SnPMcfg.mat CfgFile...');
	drawnow
	CWD = spm_str_manip(tmp,'hd');
end
if strcmp(CWD, '.')
  CWD=pwd;
end  
if ~strcmp(pwd,CWD)
	cd(CWD), CWD=pwd; fprintf('Changing directory to %s\n',CWD)
end
CfgFile = fullfile(CWD,'SnPMcfg.mat');

%-Load config file & catch all problem cases now
%-----------------------------------------------------------------------
load(CfgFile);
if isempty([H C])
    error('No model specified; [H C] empty'); end
if ~isempty(H) & ~isempty(C)
    error('Cannot have both heirachical and covariate effects'); end
if size(CONT,2) ~= size([H C B G],2)
    error('Contrast problem; wrong number of columns'); end
if size(CONT,1) > 1
    warning('F contrast!  F statistic images are being created.'); 
    STAT = 'F';
    if (CONT(1,:) == -CONT(2,:))
      CONT = CONT(1,:);
    end
else
    STAT = 'T';
end
if rank(CONT)<size(CONT,1)
  [u s] = spm_svd(CONT'); % Kill zero-rank components
  CONT = full(u*sqrt(s))';
end
if ~bVolm & bVarSm & vFWHM(3)
    error('Cannot z-smooth variance in non-volumetric mode'); end
if exist('bVarAlph')~=1, bVarAlph=0; end
if bVarAlph & ~(~bVarSm & bVolm)
    error('No pseudo t or nonvolumetric w/ variable alpha');
end
if ~bVolm & pU_ST_Ut>=0
    error('Must work volumetrically to computer STCS on-the-fly');
end
% Re-map files to avoid Endian headaches
if defaults.analyze.multivol
  for i = 1:length(V)
    V(i) = spm_vol([V(i).fname ',' num2str(V(i).n)]);
  end
else
  V = spm_vol(str2mat(V(:).fname));
end

%-Delete files from previous analyses, if they exist
%-----------------------------------------------------------------------
files = {	'^ResMS\..{3}$','^beta_.{4}\..{3}', '^lP_.{4}\..{3}',...
		'^lP.{1}\..{3}','^snpm.{2}\..{3}','^snpm.{1}\..{3}'};

for i=1:length(files)
    j = spm_select('List',pwd,files{i});
    for k=1:size(j,1)
        spm_unlink(deblank(j(k,:)));
    end
end

spm_unlink SnPM.mat SnPM_ST.mat SnPMt.mat SnPMucp.mat XYZ.mat SnPM_pp.mat ...
    SnPM_pp_Neg.mat STCS.mat


%-Parameters & Initialisation
%=======================================================================

%-Suprathreshold parameters
%-----------------------------------------------------------------------
global SnPMdefs
if isempty(SnPMdefs), snpm_defaults; end
STalpha = SnPMdefs.STalpha; 
STprop  = SnPMdefs.STprop;

s_SnPM_save = [s_SnPM_save ' STalpha STprop'];  % Save for PP

%-Work out degrees of freedom
%-----------------------------------------------------------------------
q       = size([H C B G],1);		%-# observations
p       = size([H C B G],2);		%-# predictors
r       = rank([H C B G]);		%-Model degrees of freedom
df      = q - r;			%-Residual degrees of freedom
nPerm   = size(PiCond,1);		%-# permutations

%-Get ORIGIN, etc
%-----------------------------------------------------------------------
DIM    = [V(1).dim(1)   V(1).dim(2)   V(1).dim(3)];
M=V(1).mat(1:3, 1:3);
VOX=sqrt(diag(M'*M))';
MAT    = V(1).mat;
IMAT   = inv(MAT);
ORIGIN = IMAT(1:3,4);


%-Var-alpha stuff
%-----------------------------------------------------------------------
bMask = 0;
if exist('bVarAlph')==1 & bVarAlph
  Vwt    = spm_vol(Pwt);
  MinwP  = repmat(Inf,nPerm,2);
  s_SnPM_save = [s_SnPM_save ' MinwP Pwt Vwt'];
  bMask = 1;
elseif exist('Pwt')
  Vwt    = spm_vol(Pwt);
  s_SnPM_save = [s_SnPM_save ' Pwt Vwt'];
  bMask = 1;
end


%-Useful quantities - handy for later
%-----------------------------------------------------------------------
xdim     = DIM(1);			%-X dimension
ydim     = DIM(2);			%-Y dimension
zdim     = DIM(3);			%-Z dimension
PlDim    = xdim*ydim;			%-Plane size in voxels
VolDim   = xdim*ydim*zdim;		%-Volume size in voxels
if bVolm, WorkDim = VolDim;		%-Working dimension (if volumetric)
	else, WorkDim = PlDim; end	%-Working dimension (if plane by plane)

%-Location vectors --> In units of mm <--
%-----------------------------------------------------------------------
[y x] = meshgrid([1:ydim],[1:xdim]');
x     = x(:);
y     = y(:);
z     = (1:zdim)'; %-Note diff from spm_spm.m

%-Initialize variables
%-----------------------------------------------------------------------
TH    = TH*ones(1,WorkDim);			%-Global activities
S     = 0;					%-Volume analyzed
MaxT  = repmat(-Inf,nPerm,2);			%-Max t
nP    = zeros(1,WorkDim);			%-Nonparam P's
XYZ_total=[];                                   %-the variable for keeping all XYZ
%-If working plane by plane, preallocate Q & XYZ for speed/mem. efficiency
if ~bVolm, 
  Q    = zeros(1,PlDim);
  XYZ  = zeros(3,PlDim);
end

SmTime = 0;					%-Smoothing time
perm   = 0;

% Initialize structure template
%----------------------------------------
Vt=V(1);

%
%-Initialize image structures.
%
for ii=1:p
  fname= sprintf('beta_%04d.img',ii);
  descrip=sprintf('beta_%04d hats',ii);
  Vbeta(ii)=snpm_clone_vol(Vt,fname,descrip);
end  
Vbeta = spm_create_vol(Vbeta);

VResMS=snpm_clone_vol(Vt,'ResMS.img','Residual sum-of-squares');
VResMS=spm_create_vol(VResMS);
if bVarSm==0
  str = sprintf('%c_{%d} statistic',STAT,df);
else
  if STAT=='T'
    str = sprintf('SmVar T_{%d} statistic, %fx%fx%f VarSm',df,vFWHM);
  elseif STAT=='F'
    str = sprintf('SmVar F_{%d,%d} statistic, %fx%fx%f VarSm',...
		  [rank(CONT) df],vFWHM);
  end
end
%%% Vst = sf_InitStatImg(Vt,STAT,str);  Vst.T_pos Vst.T_pos = ...;Vst.T_neg = [];
if STAT=='T'
  VT_pos=snpm_clone_vol(Vt,'snpmT+.img',[str,' (+ve)']);
  VT_pos=spm_create_vol(VT_pos);
  VT_neg=snpm_clone_vol(Vt,'snpmT-.img',[str,' (-ve)']);
  VT_neg=spm_create_vol(VT_neg);
elseif STAT=='F'
  VF=snpm_clone_vol(Vt,'snpmF.img',str);
  VF=spm_create_vol(VF);
end

VlP_pos=snpm_clone_vol(Vt, 'lP+.img', '-log10(uncor. non-para. P, +ve)');
VlP_pos=spm_create_vol(VlP_pos);
VlP_FWE_pos=snpm_clone_vol(Vt, 'lP_FWE+.img','-log10(FWE-corr. P, +ve)');
VlP_FWE_pos=spm_create_vol(VlP_FWE_pos);
VlP_FDR_pos=snpm_clone_vol(Vt, 'lP_FDR+.img','-log10(FDR-corr. P, +ve)');
VlP_FDR_pos=spm_create_vol(VlP_FDR_pos);

if STAT=='T'
  VlP_neg=snpm_clone_vol(Vt, 'lP-.img', '-log10(uncor. non-para. P, -ve)');
  VlP_neg=spm_create_vol(VlP_neg);
  VlP_FWE_neg=snpm_clone_vol(Vt, 'lP_FWE-.img','-log10(FWE-corr. P, -ve)');
  VlP_FWE_neg=spm_create_vol(VlP_FWE_neg);
  VlP_FDR_neg=snpm_clone_vol(Vt, 'lP_FDR-.img','-log10(FDR-corr. P, -ve)');
  VlP_FDR_neg=spm_create_vol(VlP_FDR_neg);
end

if bVarAlph
  VlwP=snpm_clone_vol(Vt, 'lwP.img','-log10(weighted p-value)');
  VlwP=spm_create_vol(VlwP);
end  

%	
%-Initialize image data.
%
%BETA_image=repmat(NaN,p,WorkDim);
%ResSS_image=repmat(NaN,1,WorkDim);
%T_pos_image=repmat(NaN,1,WorkDim);
%T_neg_image=repmat(NaN,1,WorkDim);
lP_pos_image=repmat(NaN,1,VolDim);
lP_FWE_pos_image=repmat(NaN,1, VolDim);
lP_FDR_pos_image=repmat(NaN,1, VolDim);
if STAT=='T'
  lP_neg_image=repmat(NaN,1,VolDim);
  lP_FWE_neg_image=repmat(NaN,1, VolDim);	
  lP_FDR_neg_image=repmat(NaN,1, VolDim);	
end

%if bVarAlph
%  lwP_image=repmat(NaN, 1, WorkDim); 
%end
	
%=======================================================================
% - C O R R E C T   P E R M U T A T I O N
%=======================================================================
% Work out correct permuation completely. Separating the first
% permutation simplifies the permutation loop (fewer conditionals) and
% allows determination of pseudo-t threshold when saving supratheshold
% statistics.
disp('Working on correct permutation...');

SnPMt=[]; %Initialzie SnPMt,which will store the t's from correct permutation.

for i = 1:zdim
  
    %-Initialize the image data for this slice/volume
    %---------------------------------------------------------------------
    BETA_image=repmat(NaN,p,WorkDim);
    ResSS_image=repmat(NaN,1,WorkDim);
    if STAT=='T'
      T_pos_image=repmat(NaN,1,WorkDim);
      T_neg_image=repmat(NaN,1,WorkDim);
    elseif STAT=='F'
      F_image=repmat(NaN,1,WorkDim);
    end
    
    if bVarAlph
	lwP_image=repmat(NaN, 1, WorkDim); 
    end
    
    %-Form data matrix for this slice/volume
    %---------------------------------------------------------------------
    X     = zeros(q,WorkDim);
    if bMask, 
      Wt = zeros(1,WorkDim); 
    else
      Wt = 1;
    end
    if bVolm
	for j = 1:q
	    for k = 1:zdim
		tmp    = spm_slice_vol(V(j),spm_matrix([0 0 k]), ...
		    [xdim ydim],0);
		X(j,(k-1)*PlDim+1:k*PlDim) = tmp(:)';
	    end
	end
	if bMask
	  for k = 1:zdim
	    tmp    = spm_slice_vol(Vwt,spm_matrix([0 0 k]), ...
				   [xdim ydim],0);
	    tmp(~isfinite(tmp) | tmp<0) = 0;
	    Wt(1,(k-1)*PlDim+1:k*PlDim) = tmp(:)';
	  end
	end
    else
	for j = 1:q
	    tmp    = spm_slice_vol(V(j),spm_matrix([0 0 i]), ...
		[xdim ydim],0);
	    X(j,:) = tmp(:)';
	end
	if bMask
	  tmp = spm_slice_vol(Vwt,spm_matrix([0 0 i]),[xdim ydim],0);
	  tmp(~isfinite(tmp) | tmp<0) = 0;
	  Wt  = tmp(:)';
	end
    end

    %-Eliminate background voxels (based on threshold TH), and
    % eliminate voxels where there are no differences across scans.
    %---------------------------------------------------------------------
    Q = find(all(X>TH) & any(diff(X)) & Wt);
    
    if length(Q)
   
        X     = X(:,Q);
	S     = S + length(Q); 			%-Volume 
	if bVolm
	    XYZ   = [ x(rem(Q-1,PlDim)+1)          ...
		      y(rem(Q-1,PlDim)+1)            ...
		      z(ceil(Q/PlDim))      ]'; %-Locations
	else
	    XYZ   = [ x(rem(Q-1,PlDim)+1)          ...
		      y(rem(Q-1,PlDim)+1)            ...
		      z(i)*ones(length(Q),1)]';	%-Locations
	end 

	% Convert Voxels to mm's
	XYZ = MAT*[XYZ;ones(1,length(Q))]; XYZ(4,:) = [];

	if (bMask)
	  Wt    = Wt(1,Q);
	end

	perm = 1;

	%-Estimate parameters and sum of squares due to error.
	% Use pseudo inverse rather than BETA=inv(D'*D)*D'*X for 
	% D = DesMtx, to allow for non-unique designs. See matlab help.
	%-----------------------------------------------------------------
	BETA  = pinv([H C B G])*X;
	ResSS = sum((X - [H C B G]*BETA).^2);
	
	%-Variance smoothing.
	% Blurred mask is used to truncate kernal to brain; if not
	% used variance at edges would be underestimated due to
	% convolution with zero activity out side the brain.
	%-----------------------------------------------------------------
	if bVarSm
	    if bVolm
		SmResSS   = zeros(xdim, ydim, zdim);
		SmMask    = zeros(xdim, ydim, zdim);
		TmpVol    = zeros(xdim, ydim, zdim);
		TmpVol(Q) = ones(size(Q));
		spm_smooth(TmpVol,SmMask,vFWHM./VOX);
		TmpVol(Q) = ResSS;
		spm_smooth(TmpVol,SmResSS,vFWHM./VOX);

		ResSS     = SmResSS(Q)./SmMask(Q);
	    else
		TmpPl     = zeros(xdim,ydim);
		TmpPl(Q)  = ones(size(Q));
		SmMask    = spm_conv(TmpPl, vFWHM(1)/VOX(1),vFWHM(2)/VOX(2));
		TmpPl(Q)  = ResSS;
		SmResSS   = spm_conv(TmpPl, vFWHM(1)/VOX(1),vFWHM(2)/VOX(2));
		
		ResSS     = SmResSS(Q)./SmMask(Q);
	    end
	end
	
	%-Compute t-statistics for specified compounds of parameters
	%-----------------------------------------------------------
	T      = zeros(1,size(BETA,2));
	Co     = CONT;
	if STAT=='T'
	  % t, as usual
	  T(1,:) = Co*BETA./sqrt((ResSS*(Co*pinv([H C B G]'*[H C B G])*Co'))/df);
	else
	  % F!
	  pX   = pinv([H C B G]);
	  T(1,:) = (sum(((Co*BETA)'*inv(Co*pinv([H C B G]'*[H C B G])*Co'))' .* ...
			(Co*BETA),1)/rank(Co)) ./ (ResSS/df);
	end	

	%-Save Max T statistic
	%-----------------------------------------------------------
	MaxT(perm,:) = max([ max(T(1,:)), -min(T(1,:));   ...
	                    MaxT(perm,1),  MaxT(perm,2) ]);
	
	%-Save min weighted p-value
	%-----------------------------------------------------------
	if bVarAlph,
	  MinwP(perm,:) = min([ min(Wt.*(1-spm_Tcdf(T(1,:),df))),     ...
		                min(Wt.*(1-spm_Tcdf(-T(1,:),df)));    ...
		                MinwP(perm,1), MinwP(perm,2) ]);
	end

	%-Save weighted p-value (later converted into corr'd wt'd p-val)
	%-----------------------------------------------------------
	if bVarAlph,
	  wP = [Wt.*(1-spm_Tcdf( T(1,:),df));  ...
		Wt.*(1-spm_Tcdf(-T(1,:),df))];
	end

	%-Adjustment (remove effects of no interest) & save
	%-----------------------------------------------------------
	XA = X - [zeros(size([H C])) B G]*BETA;
		
	%
	%- New! Write out data images.
	%- Input image data.
	BETA_image(:,Q)=BETA;
	ResSS_image(:,Q)=ResSS;
	if STAT=='T'
	  T_pos_image(:,Q)=T;
	  T_neg_image(:,Q)=-T;
	elseif STAT=='F'
	  F_image(:,Q)=T;
	end
	
	if bVarAlph
	  lwP_image(:,Q)=-log10(wP);
	end  
	
        if bVolm
	  SnPMt=T; % save T's
	else
	  SnPMt=[SnPMt,T]; % save T's.
	end

	XYZ_total=[XYZ_total, XYZ];
	
	%if bVarAlph,
	%  spm_append_96('SnPMwP',wP);      % wt'd p-val of Stat du jour
	%end
	    
     
    end %(if length(Q))
    
    %-The image of the volume or the slice should be written out no matter length(Q)=1
    %or 0. 
    if bVolm
          for ii=1:p
	    BETA_vol=reshape(BETA_image(ii,:),DIM(1),DIM(2),DIM(3));
	    spm_write_vol(Vbeta(ii),BETA_vol);
	  end
	  
	  ResSS_vol=reshape(ResSS_image,DIM(1),DIM(2),DIM(3));
	  spm_write_vol(VResMS, ResSS_vol);
	  
	  if STAT=='T'
	    T_pos_vol=reshape(T_pos_image,DIM(1),DIM(2),DIM(3));
	    spm_write_vol(VT_pos,T_pos_vol);
	  
	    T_neg_vol=reshape(T_neg_image,DIM(1),DIM(2),DIM(3));
	    spm_write_vol(VT_neg,T_neg_vol);
	  
	  elseif STAT=='F'
	    F_vol=reshape(F_image,DIM(1),DIM(2),DIM(3));
	    spm_write_vol(VF,F_vol);
	  end
	  
	  if bVarAlph
	    lwP_vol=reshape(lwP_image, DIM(1), DIM(2), DIM(3));
	    spm_write_vol(VlwP, lwP_vol);
	  end
	  
    else  
     	  for ii=1:p
	    BETA_plate=reshape(BETA_image(ii,:), DIM(1), DIM(2));
	    spm_write_plane(Vbeta(ii),BETA_plate,i);
	  end
	  
	  ResSS_plate=reshape(ResSS_image, DIM(1), DIM(2));
	  spm_write_plane(VResMS,ResSS_plate,i);
	    
	  if STAT=='T'
	    T_pos_plate=reshape(T_pos_image, DIM(1), DIM(2));
	    spm_write_plane(VT_pos,T_pos_plate,i);
	  
	    T_neg_plate=reshape(T_neg_image, DIM(1), DIM(2));
	    spm_write_plane(VT_neg,T_neg_plate,i);
	  
	  elseif STAT=='F'
	    F_plate=reshape(F_image, DIM(1), DIM(2));
	    spm_write_plane(VF,F_plate,i);
	  end  
	    
	  if bVarAlph
	    lwP_plate=reshape(lwP_image, DIM(1), DIM(2));
	    spm_write_plane(VlwP, lwP_plate,i);
	  end
    end	  
    
    %-Whole volume complete in one pass if volumetric
    if bVolm, break; end
end

%close the images.
% for ii=1:p
%   spm_close_vol(Vbeta(ii));
% end  
% spm_close_vol(VResMS);

% if STAT=='T'
%   spm_close_vol(VT_pos);
%   spm_close_vol(VT_neg);
% elseif STAT=='F'
%   spm_close_vol(VF);
% end

% if bVarAlph
% spm_close_vol(VlwP);
% end

% Make an error if actually 'no voxels in brain'. 
if perm==0, error('No voxels in brain'); end

save SnPMt SnPMt

%save XYZ in a XYZ.mat file.
%===============
XYZ=XYZ_total;
save XYZ XYZ


%-Set SupraThreshold t-threshold
%=======================================================================
if bST 
  if pU_ST_Ut==-1 % No threshold has been set yet.
    if bVarSm
	load SnPMt
	SnPMt = sort(SnPMt')';
	ST_Ut = SnPMt(round((1-STprop)*length(SnPMt)));
	% clear SnPMt
    else
	ST_Ut = spm_invTcdf(1 - STalpha, df);
    end
  else            % A threshold has been set.    
    if bVarSm
        ST_Ut=pU_ST_Ut;
    else
       if (pU_ST_Ut>1)
         ST_Ut=pU_ST_Ut;
       else
         ST_Ut=spm_invTcdf(1-pU_ST_Ut, df);
       end   
    end 
  end  
  StartPerm = 1;   % redo 1st perm for ST stats
else
  ST_Ut = Inf;
  StartPerm = 2;
end

%-Save correctly labeled T's
if bVolm & (StartPerm==2)
  T0 = T;
  nPtmp = ones(size(T));
else
  StartPerm = 1;
  nPtmp=[];
end


%=======================================================================
% - C O M P U T E   F O R   P E R M U T A T I O N S
%=======================================================================
%-Cycle over planes (or just once for volumetric mode)

%-If working plane by plane, preallocate Q & XYZ for speed/mem. efficiency
if ~bVolm, 
  Q    = zeros(1,PlDim);
  XYZ  = zeros(3,PlDim);
end

%-Setup progress bar
if bWin & ~bVolm
    spm_progress_bar('Init',zdim,'Looping over (perms within) planes...','Plane')
elseif bWin
    spm_progress_bar('Init',nPerm,'Volumetric mode...','Permutation')
    spm_progress_bar('Set',StartPerm-1)
end
tic %-Start the clock: Timing code is commented with "clock" symbol: (>)

%-Loop over planes (breaks out after first loop if bVolm)
%-----------------------------------------------------------------------
nP = [];
for i = 1:zdim
    
    PlStart=toc;SmTime=0; %-Timestamp (>) 
	
    if bVolm, disp('Working on the whole volume');
	else, fprintf('\tPlane %3d: ',i); end

    %-Form data matrix for this slice (done in correctPerm code above if bVolm)
    %---------------------------------------------------------------------
    if ~bVolm
	X     = zeros(q,PlDim);
	for j = 1:q
	    tmp = spm_slice_vol(V(j),spm_matrix([0 0 i]),[xdim ydim],0);
	    X(j,:) = tmp(:)';
	end

	%-Eliminate background voxels (based on global threshold TH),
	% and eliminate voxels where there are no differences across scans.
	%-----------------------------------------------------------------
	Q = find(all(X > TH) & any(diff(X)));
    end % (if ~bVolm)
    
    if length(Q)
	if ~bVolm, X = X(:,Q); end	%-Already done if bVolm

	if bST & ~bVolm			%-XYZ already done if bVolm
	    XYZ   = [ x(rem(Q-1,PlDim)+1)          ...
		      y(rem(Q-1,PlDim)+1)          ...
		      z(i)*ones(length(Q),1)]';	%-Locations
	    XYZ = MAT*[XYZ;ones(1,length(Q))]; 
            XYZ(4,:) = [];
	end 
		
	if bVarSm & ~bVolm			%-Smoothing & plane-by-plane
	    SmStart = toc;			%-Timestamp (>)
	    TmpPl     = zeros(xdim,ydim);
	    TmpPl(Q)  = ones(size(Q));
	    SmMask     = spm_conv(TmpPl, vFWHM(1)/VOX(1),vFWHM(2)/VOX(2));
	    SmTime = SmTime + toc-SmStart;	%-Timestamp (>)
	end
	
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Initialize structure STCS 
	if bST & pU_ST_Ut>=0
	  STCS = snpm_STcalc('init',nPerm); 
	end
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	
	
	%-Loop over permutations
	%-----------------------------------------------------------------
	for perm = StartPerm:nPerm
	    PmStart = toc;			%-Timestamp (>)
	    
	    if bVolm, SmTime=0;			%-Timestamp (>)
		else, clear T BETA ResSS; end	%-Clean up
	    
	    %-Rebuild H C for current permuation
	    %-----------------------------------------------------------
	    HC = eval(sHCform);
	    
	    %-Estimate parameters and sum of squares due to error
	    %-Use pseudo inverse rather than BETA=inv(D'*D)*D'*X
	    % for D = DesMtx, to allow for non-unique designs. 
	    % See matlab help.
	    %-----------------------------------------------------------
	    BETA  = pinv([HC B G])*X;
	    ResSS = sum((X - [HC B G]*BETA).^2);
	    
	    if bVarSm
		SmStart=toc;			%-Timestamp (>)
		if bVolm
		    TmpVol(Q) = ResSS;
		    spm_smooth(TmpVol,SmResSS,vFWHM./VOX);
		    ResSS     = SmResSS(Q)./SmMask(Q);
		else
		    TmpPl(Q)  = ResSS;
		    SmResSS   = spm_conv(TmpPl,vFWHM(1)/VOX(1),vFWHM(2)/VOX(2));
		    ResSS     = SmResSS(Q)./SmMask(Q);
		end
		SmTime = SmTime + toc-SmStart;	%-Timestamp (>)
	    end
	    
	    %-Compute t-statistics for specified contrast of parameters
	    %-----------------------------------------------------------
	    T      = zeros(1,size(BETA,2));
	    Co     = CONT;
	    if STAT=='T'
	      % t, as usual
	      T(1,:) = Co*BETA./sqrt((ResSS*(Co*pinv([HC B G]'*[HC B G])*Co'))/df);
	    else
	      % F!
	      pX   = pinv([HC B G]);
	      T(1,:) = (sum(((Co*BETA)'*inv(Co*pinv([HC B G]'*[HC B G])*Co'))' .* ...
			    (Co*BETA),1)/size(Co,1)) ./ (ResSS/df);
	    end	
	    	   
	    
	    %-Save Max T statistic
	    %-----------------------------------------------------------
	    MaxT(perm,:) = max([ max(T(1,:)), -min(T(1,:));      ...
		                MaxT(perm,1), MaxT(perm,2) ]);
	    
	    %-Update nonparametric P-value
	    %-----------------------------------------------------------
	    if (perm==1)
	      T0 = T;
	      nPtmp = ones(size(T));
	      if bhPerms
		nPtmp = nPtmp + (T0<0);
	      end
	    else
	      if bhPerms
		nPtmp = nPtmp + (T>=T0) + (-T>=T0);   % NB: Worry if T0=T=0
	        % if STAT=='T', then T, T0 >=0, so (-T>=T0) will be empty.
	      else
		nPtmp = nPtmp + (T>=T0);
	      end
	    end
	    
	    %-Save min weighted p-value
	    %-----------------------------------------------------------
	    if bVarAlph,
	      MinwP(perm,:) = min([ min(Wt.*(1-spm_Tcdf( T(1,:),df))),     ...
		                    min(Wt.*(1-spm_Tcdf(-T(1,:),df)));    ...
		                    MinwP(perm,1), MinwP(perm,2) ]);
	    end

	    %-Save T,XYZ,perm for suprathreshold analysis
	    %-----------------------------------------------------------
	    %if bST & pU_ST_Ut==-1
            if bST & (pU_ST_Ut~= -1),
	        clear d1 d2
		d1 = find(T(1,:) >  ST_Ut);
		d2 = find(T(1,:) < -ST_Ut);
		spm_append_96('SnPM_ST',[                            ...
		    XYZ(:,d1),               XYZ(:,d2);              ...
		    T(1,d1),                 -T(1,d2);               ...
		    perm*ones(1,length(d1)), -perm*ones(1,length(d2)) ...
			]);
        end
	          
	    if bST & pU_ST_Ut>=0
		clear d1 d2 SnPM_ST_Pos SnPM_ST_Neg
		d1 = find(T(1,:) >  ST_Ut);
		d2 = find(T(1,:) < -ST_Ut);
		SnPM_ST_Pos=[               ...
		    XYZ(:,d1);              ...
		    T(1,d1)];  
		
		SnPM_ST_Neg=[               ...
		    XYZ(:,d2);              ...
		    -T(1,d2)];
		
		% Negate perm numbers if looking at negative contrast
		%SnPM_ST_Neg=SnPM_ST;
		%SnPM_ST_Neg(5,:)=-SnPM_ST_Neg(5,:);
		
		%if bhPerms
		%%-Renumber negative perms according to -flipud PiCond
		%tQ = SnPM_ST(5,:)<0;
		%SnPM_ST(5,tQ) = 2*nPerm +1 +SnPM_ST(5,tQ);
		%tQ_Neg = SnPM_ST_Neg(5,:)<0;
		%SnPM_ST_Neg(5,tQ_Neg) = 2*nPerm +1 +SnPM_ST_Neg(5,tQ_Neg);
		
		%else
		%%-Not bhPerms: Lose entries for negative excursions
		%SnPM_ST = SnPM_ST(:,SnPM_ST(5,:)>0);
		%SnPM_ST_Neg = SnPM_ST_Neg(:,SnPM_ST_Neg(5,:)>0);
		%end
		
		if STAT== 'F'
		  loop = 1;
		else
		  loop = 1:2;
		end
				
		for isPos = loop %1 for positive; 2 for negative
		  if isPos==1
		    SnPM_ST = SnPM_ST_Pos;
		  else
		    SnPM_ST = SnPM_ST_Neg;
		  end
		
		  % consider Permuation NO. perm
		  %  tQ1=(SnPM_ST(5,:)==perm);
		  if any(SnPM_ST)
		    Locs_mm=SnPM_ST(1:3,:);
		    Locs_mm (4,:) = 1;
		    Locs_vox = IMAT * Locs_mm;
		    
		    %STCS = snpm_STcalc('update',STCS, SnPM_ST(4,:),...
			%	       Locs_vox(1:3,:),isPos,perm);
            STCS = snpm_STcalc('update',STCS, SnPM_ST(4,:),...
                       Locs_vox(1:3,:),isPos,perm,pU_ST_Ut,df);
		    		    
		    %save perm 1 stats for use later -[X;Y;Z;T;perm;STCno]
		    if (perm==1)
		      tmp = spm_clusters(Locs_vox(1:3,:));
		      STCstats=[SnPM_ST;perm*ones(1,size(SnPM_ST,2));tmp];
		      if isPos==1
			save SnPM_pp STCstats
		      else
			STCstats_Neg = STCstats;
			save SnPM_pp_Neg STCstats_Neg
		      end
		    end			
		  end  %if any(SnPM_ST) 

		end  %for isPos=loop    
		  		
	    
	    end
	    
	    %-Print status at each perm if bVolm (& stop maybe)
	    %-----------------------------------------------------------
	    if bVolm
		if bWin, spm_progress_bar('Set',perm), end
		%-Printout timing information (>)
		fprintf('\tPerm %4d: %3d''%2d" (%3d%% Sm)\n', perm, ...
		    floor((toc-PmStart)/60),round(rem(toc-PmStart,60)), ...
		    round(100*SmTime/(toc-PmStart)))

	    end % (if bVolm)
	end 	% (for perm = StartPerm:nPerm) - Perm loop
	
	  
	  %- save STCS
	  if bST & pU_ST_Ut>=0
	    if bhPerms %Double the STCS variables.
	      STCS = snpm_STcalc('double',STCS);
	    end
	     	    
	    save STCS STCS
	  end
	  
    end 	% (length(Q)) - Conditional on non-zero voxels
    nP = [nP, nPtmp];
    nPtmp = [];

    if bVolm 
	break; end
    
    %-Print status at each plane if ~bVolm
    %-----------------------------------------------------------
    if bWin, spm_progress_bar('Set',i), end
    %-Printout timing information (>)
    fprintf('%3d''%2d" (%3d%% Sm)\n', ...
	floor((toc-PlStart)/60),round(rem(toc-PlStart,60)), ...
	round(100*SmTime/(toc-PlStart)))
    
end		% (for i = 1:zdim) - loop over planes

fprintf('\n\nPermutations are done. Writing out images.\n')

if bhPerms
  nP = nP/(2*nPerm);
else
  nP = nP/nPerm;
end
SnPMucp=nP;
save SnPMucp SnPMucp

%
% - write out lP+ and lP- images;
%
nP_pos=nP;
lP_pos=-log10(nP_pos);
lP_pos_image(spm_xyz2e(XYZ_total, Vt))=lP_pos;
lP_pos_vol=reshape(lP_pos_image,DIM(1),DIM(2),DIM(3));
spm_write_vol(VlP_pos, lP_pos_vol);

if STAT == 'T'
  if bhPerms
    nP_neg=1+1/(2*nPerm)-nP;
  else 
    nP_neg=1+1/nPerm-nP;
  end

  lP_neg=-log10(nP_neg);
  lP_neg_image(spm_xyz2e(XYZ_total, Vt))=lP_neg;
  lP_neg_vol=reshape(lP_neg_image,DIM(1),DIM(2),DIM(3));
  spm_write_vol(VlP_neg, lP_neg_vol);
end

%
% - write out lP_FWE+ and lP_FWE- images;
%
tol = 1e-4;	% Tolerance for comparing real numbers

cP_pos=zeros(size(nP));

if bhPerms
	MaxT   = [ MaxT; flipud(fliplr(MaxT)) ];
end
MaxT_pos=MaxT(:,1);

for t = MaxT_pos'
	%-FEW-corrected p is proportion of randomisation greater or
	% equal to statistic.
	%-Use a > b -tol rather than a >= b to avoid comparing
	% two reals for equality.
	cP_pos = cP_pos + (t > SnPMt -tol);
end

if STAT =='T'
  cP_neg=zeros(size(nP));
  MaxT_neg=MaxT(:,2);

  for t = MaxT_neg'
        cP_neg = cP_neg + (t > -SnPMt -tol);
  end
end

if bhPerms
  cP_pos = cP_pos / (2* nPerm);
  if STAT=='T'	
    cP_neg = cP_neg / (2* nPerm);
  end
else 
  cP_pos = cP_pos / nPerm;  
  if STAT=='T'
    cP_neg = cP_neg / nPerm;
  end
end

lP_FWE_pos=-log10(cP_pos);
lP_FWE_pos_image(spm_xyz2e(XYZ_total, Vt))=lP_FWE_pos;
lP_FWE_pos_vol=reshape(lP_FWE_pos_image,DIM(1),DIM(2),DIM(3));
spm_write_vol(VlP_FWE_pos, lP_FWE_pos_vol);

if STAT=='T'
  lP_FWE_neg=-log10(cP_neg);
  lP_FWE_neg_image(spm_xyz2e(XYZ_total, Vt))=lP_FWE_neg;
  lP_FWE_neg_vol=reshape(lP_FWE_neg_image,DIM(1),DIM(2),DIM(3));
  spm_write_vol(VlP_FWE_neg, lP_FWE_neg_vol);
end

%
% - write out lP_FDR+ and lP_FDR- images;
%
[snP_pos,I_pos]=sort(nP_pos);

Pfdr_pos=snpm_P_FDR([],[],'P',[],snP_pos');

Pfdr_pos(I_pos) = Pfdr_pos;

lP_FDR_pos=-log10(Pfdr_pos);

lP_FDR_pos_image(spm_xyz2e(XYZ_total, Vt))=lP_FDR_pos;

lP_FDR_pos_vol=reshape(lP_FDR_pos_image,DIM(1),DIM(2),DIM(3));
spm_write_vol(VlP_FDR_pos, lP_FDR_pos_vol);

if STAT =='T'
  [snP_neg,I_neg]=sort(nP_neg);
  Pfdr_neg=snpm_P_FDR([],[],'P',[],snP_neg');
  Pfdr_neg(I_neg) = Pfdr_neg;
  lP_FDR_neg=-log10(Pfdr_neg);
  lP_FDR_neg_image(spm_xyz2e(XYZ_total, Vt))=lP_FDR_neg;

  lP_FDR_neg_vol=reshape(lP_FDR_neg_image,DIM(1),DIM(2),DIM(3));
  spm_write_vol(VlP_FDR_neg, lP_FDR_neg_vol);
end


if bWin, spm_progress_bar('Clear'), end
%-Printout final timing information (>)
fprintf('\n\nThe run took %0.02f minutes\n', toc/60);

%-Cleanup
clear X

%-Save key variables
%=======================================================================
eval(['save SnPM ',s_SnPM_save])

%-Print quick summary info (allowing for STOPing)
%=======================================================================
if bhPerms
    Rank = sum([MaxT(1:perm,1);MaxT(1:perm,2)] >= MaxT(1,1));
else
    Rank = sum([MaxT(1:perm,1)] >= MaxT(1,1));
end
fprintf(['\nCorrect Perm has max t %g & rank %d out of %d ', ...
	'completed permutations\n'],MaxT(1,1),Rank,perm*(bhPerms+1));
fprintf('\n\tRun snpm_pp for full results\n\n');



function Vs = sf_close_vol(Vs)
% Don't need to close images in SPM5
return