function [df, p]=fmrilm(input_file,output_file_base,X_cache, ...
    contrast,exclude,which_stats,fwhm_rho,n_poly, ...
    reference_delay, contrast_is_delay, confounds, numlags)

%FMRILM
%
% Reads in an fMRI minc file and writes out a T statistic minc file. 
% The method is based on linear models with correlated AR(p) errors:
% Y=hrf*Xb+e, e_t=a_1 e_(t-1) + ... + a_p e_(t-p) + white noise_t. 
% 
% [DF,P] = FMRILM( INPUT_FILE, OUTPUT_FILE_BASE, X_CACHE,  
%    CONTRAST [, EXCLUDE [, WHICH_STATS [, FWHM_RHO [, N_POLY 
%    [, REFERENCE_DELAY [, CONTRAST_IS_DELAY [, CONFOUNDS [,NUMLAGS]]]]]]]]] )
% 
% INPUT_FILE is the input fmri file, the dependent variable, with extension
% .img for ANALYZE format, or .mnc for MINC foramt.
% Scans must be equally spaced in time, unless FWHM_RHO=Inf (see below).
% 
% OUTPUT_FILE_BASE: matrix whose rows are the base for output statistics,
% one base for each row of CONTRAST, padded with extra blanks if 
% necessary; they will be removed before use.
%
% X_CACHE: The design matrices, usually from fmridesign.m, as a 4D array. 
% Dim 1: frames; Dim 2: response variables; 
% Dim 3: derivative + 1; Dim 4: slices. 
%
% CONTRAST is a matrix whose rows are contrasts for the
% response variables, and columns are the contrast. Extra columns can 
% be added to estimate contrasts in the polynomial terms and the confounds 
% (in that order - see N_POLY and CONFOUNDS below). 
% 
% EXCLUDE is a list of frames that should be excluded from the
% analysis. This must be used with Siemens EPI scans to remove the
% first few frames, which do not represent steady-state images.
% If NUMLAGS=1, the excluded frames can be arbitrary, otherwise they 
% must be from the beginning and/or end. Default is [1].
% 
% WHICH_STATS: logical matrix indicating output statistics by 1.
% Rows correspond to rows of CONTRAST, columns correspond to:
%   1: T statistic image, OUTPUT_FILE_BASE_mag_t.mnc or .img (for contrasts 
%      in the delays, mag is replaced by del in the extension).  
%      The degrees of freedom is DF. Note that tstat=effect/sdeffect. 
%   2: effect image, OUTPUT_FILE_BASE_mag_ef.mnc or .img.
%   3: standard deviation of the effect, OUTPUT_FILE_BASE_mag_sd.mnc or .img. 
%   4: F-statistic image, OUTPUT_FILE_BASE_mag_F.mnc or .img, for testing
%      all rows of CONTRAST that deal with magnitudes (see CONTRAST_IS_DELAY
%      below). The degrees of freedom are [DF, P].
%   5: the temporal autocorrelation(s), OUTPUT_FILE_BASE_rho.mnc or .img.
%   6: the residuals from the model, OUTPUT_FILE_BASE_resid.mnc or .img,
%      only for the non-excluded frames (warning: uses lots of memory). 
%   7: the whitened residuals from the model, OUTPUT_FILE_BASE_wresid.mnc or
%      .img, only for the non-excluded frames (warning: uses lots of memory).
%   8: the AR parameter(s) a_1 ... a_p, in OUTPUT_FILE_BASE_A.mnc or .img.
% If WHICH_STATS to a row vector of length 8, then it is used for all
% contrasts; if the number of columns is less than 8 it is padded
% with zeros. Default is 1. 
%   
% FWHM_RHO is the fwhm in mm of a 3D Gaussian kernel used to smooth the
% autocorrelation. This is NOT the fwhm of the image data as used
% by tstat_threshold, but should be fairly large. Setting it to Inf
% smooths the autocorrelation to 0, i.e. it assumes the scans
% are uncorrelated (useful for TR>10 seconds). Default is 15. 
% If FWHM_RHO is a character string, then this file is used for the
% autocorrelations, e.g. the _rho.mnc or .img file created by a previous run.
% 
% N_POLY is the order of polynomial trend terms to remove
% (N_POLY>=0); =0 will remove just the constant term and no trends.
% -1 will not remove anything, in which case the design matrix is
% completely determined by X. Default is 3.
%
% REFERENCE_DELAY: A row vector of reference delays for the reponses.
% If it is >0 then the model allows for variable delays; if it is =0
% then the delay is fixed by the reference hrf and delays are not
% estimable. Theory: We start with a reference hrf h(t) close
% to the true hrf, e.g. the Glover (1999) hrf used by FMRIDESIGN.
% The hrf for each variable is then modeled as  h( t / c ) / c where 
% c = DELAY / REFERENCE_DELAY, i.e. the time is re-scaled by c.
% The program then estimates contrasts in the c's (specified by 
% CONTRAST; contrasts for the polynomial drift terms are ignored),
% their standard deviations, and T statistics for testing c=1 
% (no change). Output is DELAY = c * REFERENCE_DELAY in 
% OUTPUT_FILE_BASE_del_ef.mnc, sd is in OUTPUT_FILE_BASE_del_sd.mnc,
% and t statistic in OUTPUT_FILE_BASE_del__t.mnc (or .img). Note that
% REFERENCE_DELAY is only used to multiply c and its sd for
% output. If REFERENCE_DELAY is set equal to the first 
% element of HRF_PARAMETERS (5.4 seconds by default) then DELAY is 
% the estimated time to the first peak of the hrf (in seconds).
% To do the estimation, the stimuli must be convolved with first 
% and second derivatives of h with respect to log(c) at c=1
% and added to X_CACHE. FMRIDESIGN does this for the Glover (1999) hrf.
% If REFERENCE_DELAY is a scalar, it is used for all variables.
% Note that F statistics are not yet available with this option.
% Default is 0, i.e. all delays are fixed and not estimable.
%
% CONTRAST_IS_DELAY is a logical vector indicating if the row of CONTRAST
% is a contrast in the delays (1) or magnitudes (0). Note that you cannot
% estimate delays of the polynomial terms or confounds, and that if
% you estimate delays (1) then each non-zero element of CONTRAST must
% be estimable i.e. reference delay >0 for that response. If the length
% of CONTRAST_IS_DELAY is less then the number of contrasts, it is padded
% with zeros. Note that the t test for contrasts with a single number tests 
% that the delay is equal to the reference delay (not that it equals zero). 
% Default is 0, i.e. all contrasts refer to magnitudes.
% 
% CONFOUNDS: Extra columns for the design matrix that are the same for
% every slice, e.g. movement artifacts. Default is [], i.e. no confounds.
% 
% NUMLAGS is the order (p) of the autoregressive model. Default is 1.
%
% DF is the (residual) degrees of freedom of the statistics.
%
% P is the degrees of freedom of the set of CONTRAST_MAGNITUDES.

%############################################################################
% COPYRIGHT:   Copyright 2001 K.J. Worsley and C. Liao, 
%              Department of Mathematics and Statistics,
%              McConnell Brain Imaging Center, 
%              Montreal Neurological Institute,
%              McGill University, Montreal, Quebec, Canada. 
%              worsley@math.mcgill.ca, liao@math.mcgill.ca
%
%              Permission to use, copy, modify, and distribute this
%              software and its documentation for any purpose and without
%              fee is hereby granted, provided that the above copyright
%              notice appear in all copies.  The author and McGill University
%              make no representations about the suitability of this
%              software for any purpose.  It is provided "as is" without
%              express or implied warranty.
%############################################################################

% Defaults:

if nargin < 5
    exclude=[1]
end
if nargin < 6
    which_stats=1
end
if nargin < 7
    fwhm_rho= 15
end
if nargin < 8
    n_poly=3
end
if nargin < 9
    reference_delay=0.0
end
if nargin < 10
    contrast_is_delay=0
end
if nargin < 11
    confounds=[]
end
if nargin < 12
    numlags=1
end

% Open image:

d=fmris_read_image(input_file,0,0);
d.dim
numframes=d.dim(4);
numslices=d.dim(3);
numys=d.dim(2);
numxs=d.dim(1);
numpix=numxs*numys;
Steps=d.vox
file_path=d.file_path;
ext=input_file((length(input_file)-2):length(input_file));

if strcmp(computer,'PCWIN')
   if strcmp(output_file_base(2:3),':\');file_path='';end
else
   if strcmp(output_file_base(1),'/');file_path='';end
end


% Keep time points that are not excluded:

allpts = 1:numframes;
allpts(exclude) = zeros(1,length(exclude));
keep = allpts( find( allpts ) );
n=length(keep)

% Create polynomial trends:

if n_poly>=0 
    trend=((2*keep-(max(keep)+min(keep)))./(max(keep)-min(keep)))';
    Trend=(trend*ones(1,n_poly+1)).^(ones(n,1)*(0:n_poly));
else
    Trend=[];
end 

% Add confounds:

if ~isempty(confounds)
    Trend=[Trend confounds(keep,:)];
end
numtrends=size(Trend,2)
numresponses=size(X_cache,2)

% Make full contrasts:

numcontrasts=size(contrast,1)
contrast_is_delay=[contrast_is_delay ...
        zeros(1,numcontrasts-length(contrast_is_delay))]
contrasts=[contrast zeros(numcontrasts,numresponses+numtrends-size(contrast,2))]
p=rank(contrast(~contrast_is_delay,:))

% Check for estimability:

tolerance=0.0000001;
for slice=1:numslices
    if size(X_cache,4)==1;
        X=[squeeze(X_cache(keep,:,1)) Trend];
    else
        X=[squeeze(X_cache(keep,:,1,slice)) Trend];
    end
    NullSpaceX=null(X);
    Cmhalf=diag(1./sqrt(diag(contrasts*contrasts')));
    ContrastNullSpaceX=Cmhalf*contrasts*NullSpaceX;
    nonest=sum(abs(ContrastNullSpaceX)>tolerance,2);
    if sum(nonest)>0
        fprintf(['Error: the following contrasts are nonestimable in slice ' ...
                num2str(slice) ':']);
        RowsNonEstContrast=find(nonest>0)
        NonEstContrast=contrasts(RowsNonEstContrast,:)
        NullSpaceX
        ContrastNullSpaceX
        return
    end
end

% Open files for output:

if size(which_stats,1)==1
    which_stats=repmat(which_stats,numcontrasts,1);
end
which_stats=[which_stats zeros(numcontrasts,8-size(which_stats,2))]

typestat=['_t. '; '_ef.'; '_sd.'];
for i=1:numcontrasts
    if contrast_is_delay(i)
        param='_del';
    else
        param='_mag';
    end
    for stat=1:3
        if which_stats(i,stat) 
            out_handle(i,stat).file_name=[deblank(output_file_base(i,:)) param deblank(typestat(stat,:)) ext];
            out_handle(i,stat).file_path=file_path;
            out_handle(i,stat).dim=[numxs numys numslices 0];
            out_handle(i,stat).parent_file=input_file;
            out_handle(i,stat).precision='float';  
            out_handle(i,stat).descrip='';
        end
    end
end

if which_stats(1,4)
    out_handle_Fstat.file_name=[deblank(output_file_base(1,:)) '_mag_F.' ext];
    out_handle_Fstat.file_path=file_path;
    out_handle_Fstat.dim=[numxs numys numslices 0];
    out_handle_Fstat.parent_file=input_file;
    out_handle_Fstat.precision='float';  
    out_handle_Fstat.descrip='';
end
if which_stats(1,6)
    out_handle_resid.file_name=[deblank(output_file_base(1,:)) '_resid.' ext];
    out_handle_resid.file_path=file_path;
    out_handle_resid.dim=[numxs numys numslices n];
    out_handle_resid.parent_file=input_file;
    out_handle_resid.precision='float';  
    out_handle_resid.descrip='';
end
if which_stats(1,7)
    out_handle_wresid.file_name=[deblank(output_file_base(1,:)) '_wresid.' ext];
    out_handle_wresid.file_path=file_path;
    out_handle_wresid.dim=[numxs numys numslices n];
    out_handle_wresid.parent_file=input_file;
    out_handle_wresid.precision='float';  
    out_handle_wresid.descrip='';
end
if which_stats(1,8)
    out_handle_A.file_name=[deblank(output_file_base(1,:)) '_A.' ext];
    out_handle_A.file_path=file_path;
    out_handle_A.dim=[numxs numys numslices numlags];
    out_handle_A.parent_file=input_file;
    out_handle_A.precision='float';  
    out_handle_A.descrip='';
end

% Setup for finding rho:

rho_slice=zeros(numxs,numys);
rho_vol=squeeze(zeros(numpix, numslices, numlags));

indk1=((keep(2:n)-keep(1:n-1))==1);
k1=find(indk1)+1;

if length(reference_delay)==1
    reference_delay=ones(1,numresponses)*reference_delay
end
response_is_delay=(reference_delay>0)

% If fwhm_rho is a file name, then read in rho, otherwise estimte it:

if isnumeric(fwhm_rho) & fwhm_rho < Inf 
    
    if which_stats(1,5)
        out_handle_rho.file_name=[deblank(output_file_base(1,:)) '_rho.' ext];
        out_handle_rho.file_path=file_path;
        out_handle_rho.dim=[numxs numys numslices numlags-(numlags==1)];
        out_handle_rho.parent_file=input_file;
        out_handle_rho.precision ='float';  
        out_handle_rho.descrip='';
    end
    
    % Smoothing rho in slice is done using conv2 with a kernel ker_xy.
    
    if fwhm_rho>0
        fwhm_x=fwhm_rho/abs(Steps(1));
        ker_x=exp(-(-ceil(fwhm_x):ceil(fwhm_x)).^2*4*log(2)/fwhm_x^2);
        ker_x=ker_x/sum(ker_x);
        fwhm_y=fwhm_rho/abs(Steps(2));
        ker_y=exp(-(-ceil(fwhm_y):ceil(fwhm_y)).^2*4*log(2)/fwhm_y^2);
        ker_y=ker_y/sum(ker_y);
        ker_xy=ker_x'*ker_y;
    else
        ker_xy=1;
    end
    
    % First loop over slices, then pixels, to get the AR parameter:
    
    Diag1=diag(indk1,1)+diag(indk1,-1);
    
    for slice=1:numslices
        slice
        if size(X_cache,4)==1
            X=[squeeze(X_cache(keep,:,1)) Trend ... 
                    squeeze(X_cache(keep,response_is_delay,2))];
        else
            X=[squeeze(X_cache(keep,:,1,slice)) Trend ... 
                    squeeze(X_cache(keep,response_is_delay,2,slice))];
        end    
        pinvX=pinv(X);
        
        % Preliminary calculations for unbiased estimates of autocorrelation:
        
        R=eye(n)-X*pinvX;
        if numlags==1
            M11=trace(R);
            M12=trace(R*Diag1);
            M21=M12/2;
            M22=trace(R*Diag1*R*Diag1)/2;
            M=[M11 M12; M21 M22];
        else
            M=zeros(numlags+1);
            for i=1:(numlags+1)
                for j=1:(numlags+1)
                    Di=(diag(ones(1,n-i+1),i-1)+diag(ones(1,n-i+1),-i+1))/(1+(i==1));
                    Dj=(diag(ones(1,n-j+1),j-1)+diag(ones(1,n-j+1),-j+1))/(1+(j==1));
                    M(i,j)=trace(R*Di*R*Dj)/(1+(i>1));
                end
            end
        end
        invM=inv(M);
        %invM=eye(numlags+1); this will undo the correction.
        
        % Loop over pixels in slice to find rho: 
        
        d=fmris_read_image(input_file,slice,1:numframes);
        img=reshape(d.data,numpix,numframes);
        Y=img(:,keep)';
        betahat_ls=pinvX*Y;
        resid=Y-X*betahat_ls;
        if numlags==1
            Cov0=sum(resid.*resid,1);
            Cov1=sum(resid(k1,:).*resid(k1-1,:),1);
            Covadj=invM*[Cov0; Cov1];
            rho_slice(:)=(Covadj(2,:)./ ...
                (Covadj(1,:)+(Covadj(1,:)<=0)).*(Covadj(1,:)>0))';
            rho_vol(:,slice)=reshape(conv2(rho_slice,ker_xy,'same'),numpix,1);   
        else
            for lag=0:numlags
                Cov(lag+1,:)=sum(resid(1:(n-lag),:).*resid((lag+1):n,:));
            end
            Covadj=invM*Cov;
            Coradj_slice= ( Covadj(2:(numlags+1),:) ...
                .*( ones(numlags,1)*((Covadj(1,:)>0)./ ...
                (Covadj(1,:)+(Covadj(1,:)<=0)))) )';
            for lag=1:numlags
                Cor_xy=reshape(Coradj_slice(:,lag),numxs,numys);
                rho_vol(:,slice,lag)=reshape(conv2(Cor_xy,ker_xy,'same'),numpix,1);   
            end
        end 
    end
    
    % Smoothing rho betwen slices is done by straight matrix multiplication
    % by a toeplitz matrix K normalized so that the column sums are 1.
    
    if fwhm_rho>0
        fwhm_z=fwhm_rho/abs(Steps(3));
        ker_z=exp(-(0:(numslices-1)).^2*4*log(2)/fwhm_z^2);
        K=toeplitz(ker_z);
        K=K./(ones(numslices)*K);
        for lag=1:numlags
            rho_vol(:,:,lag)=squeeze(rho_vol(:,:,lag))*K;
        end
    end
    
    if which_stats(1,5)         
        out_handle_rho.data=squeeze(reshape(rho_vol,numxs,numys,numslices,numlags));
        fmris_write_image(out_handle_rho);
    end
    
end

if ~isnumeric(fwhm_rho)
    out_handle_rho=fmris_read_image(fwhm_rho);
    rho_vol=squeeze(reshape(out_handle_rho.data,numpix,numslices,numlags));
end

% Make space for results. Last dim of stat_slice: 1=Tstat, 2=Effect, 3=Sd:

stat_slice=zeros(numpix,numcontrasts,3);
if which_stats(1,4)
    Fstat_slice=zeros(numpix,1);
end
if which_stats(1,6)
    resid_slice=zeros(numpix,n);
end
if which_stats(1,7)
    wresid_slice=zeros(numpix,n);
end
if which_stats(1,8)
    A_slice=zeros(numpix,numlags);
end

% Set up for second loop:

colsXnodel=[find(~response_is_delay) numresponses+(1:numtrends)]
colsX0=find(response_is_delay)
numdelays=length(colsX0)
colsX1=numresponses+numtrends+(1:numdelays)
colsX01=[colsX0 colsX1]
colsX0T=1:(numresponses+numtrends)
numcolX=numresponses+numtrends+numdelays

find_contrast_is_delay=find(contrast_is_delay)
contr_delay=contrast(find_contrast_is_delay,colsX0)
contr_delay2=repmat(contr_delay,1,2)
numcontr_delay=length(find_contrast_is_delay)
contr_delay_is_1_col=zeros(numcontr_delay,1);
find_delay_is_1_col=ones(numcontr_delay,1);
for i=1:numcontr_delay
    pos=find(contr_delay(i,:)~=0);
    if length(pos)==1
        contr_delay_is_1_col(i)=1;
        find_delay_is_1_col(i)=pos;
    end
end
contr_delay_is_1_col 
find_delay_is_1_col

find_contrast_is_mag=find(~contrast_is_delay)
contr_mag=contrasts(find_contrast_is_mag,:)
numcontr_mag=length(find_contrast_is_mag)
contr_mag_is_1_col=zeros(numcontr_mag,1);
find_mag_is_1_col=ones(numcontr_mag,1);
for i=1:numcontr_mag
    pos=find(contr_mag(i,:)~=0);
    if length(pos)==1
        contr_mag_is_1_col(i)=1;
        find_mag_is_1_col(i)=pos;
    end
end
contr_mag_is_1_col
find_mag_is_1_col

% Second loop over voxels to get statistics:

drho=0.01;
for slice=1:numslices
    slice
    if size(X_cache,4)==1       
        X=[squeeze(X_cache(keep,:,1)) Trend ... 
                squeeze(X_cache(keep,response_is_delay,2))];
        X2=squeeze(X_cache(keep,response_is_delay,3));
    else
        X=[squeeze(X_cache(keep,:,1,slice)) Trend ... 
                squeeze(X_cache(keep,response_is_delay,2,slice))];
        X2=squeeze(X_cache(keep,response_is_delay,3,slice));
    end    
    Xstar=X;
    X2star=X2;
    df=n-rank(X);
    
    d=fmris_read_image(input_file,slice,1:numframes);
    img=reshape(d.data,numpix,numframes);
    Y=img(:,keep)';
    
    if numlags==1
        % bin rho to intervals of length drho, avoiding -1 and 1:
        irho=round(rho_vol(:,slice)/drho)*drho;
        irho=min(irho,1-drho);
        irho=max(irho,-1+drho);
    else
        % use dummy unique values so every pixel is analysed seperately:
        irho=(1:numpix)';
    end
    
    for rho=unique(irho)'
        pix=find(irho==rho);
        numrhos=length(pix);
        Ystar=Y(:,pix);
        if numlags==1
            factor=1./sqrt(1-rho^2);
            Ystar(k1,:)=(Y(k1,pix)-rho*Y(k1-1,pix))*factor;
            Xstar(k1,:)=(X(k1,:)-rho*X(k1-1,:))*factor;
            X2star(k1,:)=(X2(k1,:)-rho*X2(k1-1,:))*factor;
        else
            Coradj_pix=squeeze(rho_vol(pix,slice,:));
            [Ainvt posdef]=chol(toeplitz([1 Coradj_pix']));
            nl=size(Ainvt,1);
            A=inv(Ainvt');
            if which_stats(1,8)
                A_slice(pix,1:(nl-1))=-A(nl,(nl-1):-1:1)/A(nl,nl);
            end
            B=ones(n-nl,1)*A(nl,:);
            Vmhalf=spdiags(B,1:nl,n-nl,n);
            Ystar=zeros(n,1);
            Ystar(1:nl)=A*Y(1:nl,pix);
            Ystar((nl+1):n)=Vmhalf*Y(:,pix);
            Xstar(1:nl,:)=A*X(1:nl,:);
            Xstar((nl+1):n,:)=Vmhalf*X;
            X2star(1:nl,:)=A*X2(1:nl,:);
            X2star((nl+1):n,:)=Vmhalf*X2;
        end
        pinvXstar=pinv(Xstar);
        betahat=pinvXstar*Ystar;
        resid=Ystar-Xstar*betahat;
        if which_stats(1,6)
            resid_slice(pix,:)=(Y(:,pix)-X*betahat)';
        end
        if which_stats(1,7)
            wresid_slice(pix,:)=resid';
        end
        SSE=sum(resid.^2,1);
        sd=sqrt(SSE/df);
        V=pinvXstar*pinvXstar';
        sdbetahat=sqrt(diag(V))*sd;
        T0=betahat./(sdbetahat+(sdbetahat<=0)).*(sdbetahat>0);
        
        if numdelays>0
            
            % delays, columns colsX0 (beta) and colsX1 (delta*beta):
            
            beta=betahat(colsX0,:);
            delta_beta=betahat(colsX1,:);
            inv_beta=(beta~=0)./(beta+(beta==0));
            delta_R=delta_beta.*inv_beta;
            D1=repmat(sum(pinvXstar(colsX1,:).*X2star',2),1,numrhos);
            delta_U=delta_R-0.5*delta_R.^2.*D1;
            T0_2=T0(colsX0,:).^2;
            c0=T0_2./(T0_2+1);
            delta=delta_U.*c0;
            delay=exp(delta).*repmat(reference_delay',1,numrhos);
            del_ef=contr_delay*delay;
            
            % estimates of sd of delay:
            
            gdot1=(c0./(T0_2+1)-(1-D1.*delta_R).*c0.^2).*delta_R.*inv_beta.*delay;
            gdot2=(1-D1.*delta_R).*c0.*inv_beta.*delay;
            gdot=[gdot1; gdot2]; 
            gdotc=kron(gdot,ones(1,numcontr_delay)).*repmat(contr_delay2',1,numrhos);
            Vcontr=sum((V(colsX01,colsX01)*gdotc).*gdotc,1);
            del_sd=reshape(sqrt(Vcontr),numcontr_delay,numrhos) ...
                .*repmat(sd,numcontr_delay,1);
            
            % write out delays:
            
            stat_slice(pix,find_contrast_is_delay,2)=del_ef';
            stat_slice(pix,find_contrast_is_delay,3)=del_sd';
            contr_delay_is_1_cols=repmat(contr_delay_is_1_col,1,numrhos);
            sddelta_beta=sdbetahat(colsX1,:);
            T1=delta_beta./(sddelta_beta+(sddelta_beta<=0)).*(sddelta_beta>0);
            stat_slice(pix,find_contrast_is_delay,1)= ...
                (contr_delay_is_1_cols.*T1(find_delay_is_1_col,:)+ ...
                ~contr_delay_is_1_cols.*del_ef./(del_sd+(del_sd<=0)).*(del_sd>0))';
        end
        
        if ~isempty(find_contrast_is_mag)
            
            % magnitudes:
            
            mag_ef=contr_mag*betahat(colsX0T,:);
            
            % estimates of sd of magnitudes:
            
            mag_sd=sqrt(diag(contr_mag*V(colsX0T,colsX0T)*contr_mag'))*sd;
            
            % write out magnitudes:
            
            stat_slice(pix,find_contrast_is_mag,2)=mag_ef';
            stat_slice(pix,find_contrast_is_mag,3)=mag_sd';
            contr_mag_is_1_cols=repmat(contr_mag_is_1_col,1,numrhos);
            stat_slice(pix,find_contrast_is_mag,1)= ...
                (contr_mag_is_1_cols.*T0(find_mag_is_1_col,:)+ ...
                ~contr_mag_is_1_cols.*mag_ef./(mag_sd+(mag_sd<=0)).*(mag_sd>0))';
            
            if which_stats(1,4)
                cVcinv=pinv(contr_mag*V(colsX0T,colsX0T)*contr_mag');
                SST=sum((cVcinv*mag_ef).*mag_ef,1);
                Fstat_slice(pix,1)=(SST./(SSE+(SSE<=0)).*(SSE>0)/p*df)';
            end      
        end
    end
    
    % Write out results:
    
    stat_slice(:,:,1)=min(stat_slice(:,:,1),100);
    
    for i=1:numcontrasts
        for stat=1:3
            if which_stats(i,stat) 
                out_handle(i,stat).data=reshape(stat_slice(:,i,stat),numxs,numys);
                fmris_write_image(out_handle(i,stat),slice,1);
            end
        end
    end
    if which_stats(1,4)
        Fstat_slice=min(Fstat_slice,10000);
        out_handle_Fstat.data=reshape(Fstat_slice,numxs,numys);
        fmris_write_image(out_handle_Fstat,slice,1);
    end  
    if which_stats(1,6)
        out_handle_resid.data=reshape(resid_slice,numxs,numys,n);
        fmris_write_image(out_handle_resid,slice,1:n);
    end  
    if which_stats(1,7)
        out_handle_wresid.data=reshape(wresid_slice,numxs,numys,n);
        fmris_write_image(out_handle_wresid,slice,1:n);
    end  
    if which_stats(1,8)
        out_handle_A.data=reshape(A_slice,numxs,numys,numlags);
        fmris_write_image(out_handle_A,slice,1:numlags);
    end  
    
end

% End.