function spm_adjmean_fmri_ui
% Adjusted means of fMRI via box-car General Linear Model with confounds
% FORMAT spm_adjmean_fmri_ui
%_______________________________________________________________________
%
% spm_adjmean_fmri_ui uses the General Linear Model to produce adjusted
% condition images, adjusted for global effects and confounds.
%
% This program is designed for collapsing data from a single session
% fMRI epoch-related design into a set of representative condition
% images, one for each condition, adjusted for global effects and with
% low frequency drifts removed via a discrete cosine basis set high
% pass filter. The resulting data sets are suitable for a (2nd level)
% random effects analysis of multiple sessions or subjects, or group
% comparisons.
%
% See spm_RandFX.man for further details on implementing random effects
% analyses in SPM96 using a multi-level approach.
%
%
% Overview
% ----------------------------------------------------------------------
% The program works with a single fMRI session, fitting a General
% Linear Model consisting of simple box-cars (optionally convolved with
% an estimated haemodynamic response function) and an (optional) high
% pass filter of discrete cosine basis functions. The effects of
% interest (the box-cars) are orthogonalised (residualised) with
% respect to the confounds (the high pass filter) to ensure that *all*
% confounds are removed from the data, even if they're correlated with
% the effects of interest. (For well designed experiments this makes
% little difference.) Proportional scaling and AnCova global
% normalisation are supported, the latter including the option to scale
% all the images such that their the grand mean (GM) is a specified
% value.
%
% The interface is similar to a cut-down SPM-fMRI, and the adjusted
% means are the parameter estimates from the model. The user is first
% prompted to select the scans for a single fMRI session. Then the
% epoch condition order is specified. This should be a r-vector, where
% r is the number of epochs, of integers 1:n or 0:n-1 where n is the
% number of conditions (0 can be used to indicate the baseline or
% "rest" condition. Then the number of scans per epoch is specified:
% This can be a single integer (all epochs have the same number of
% scans), or an r-vector of integers specifying number of scans for the
% corresponding epoch.
%
% Once the experimental design has been specified, the user is given
% various options: The box-cars can be convolved with an approximate
% haemodynamic reponse function; High-pass filter components can be
% added to the confounds (at a user-specified cut-off which defaults to
% twice the maximum period of the experiment); Global normalisation can
% be omitted, or implemented by proportional scaling or AnCova; and the
% grand mean scaling options specified.
%
% With the design and adjustments specified, the model is constructed,
% and the user prompted to enter/confirm filenames for the adjusted
% condition images.
%
% The model, filenames, global values and options are saved to a MatLab
% *.mat file named SPMadj.mat in the current working directory.
%
% Implicit masking is carried out: Zero voxels are implicitly assummed
% to be masked out. Thus, the adjusted mean is calculated at voxels
% which are non-zero in *all* the input images pertaining to the
% adjusted mean (usually those from the appropriate subject). (This is
% *not* a softmean.) Data realigned in a single session with SPM'96 (or
% later) are automatically implicitly zero masked with a consistent
% mask in this way.
%
% GM, the value for grand mean scaling, is user specified.
% The default value is 100.
%
% If computing adjusted means for subsequent (2nd level) modelling, as
% with a random effects analysis, then it is important to use a
% seperable model, such that the adjustment for one subject is
% independent of other subjects entered into the model. Thus,
% proportional scaling or subject-specific AnCova adjustment must be
% used. Further, multiple runs *must* use the same GM value, and should
% scale Grand mean *by subject*.
%
% ( A separate program (spm_adjmean_ui) is available for computing       )
% ( adjusted condition means of PET data. The functionality is similar   )
% ( to this code, but the two routines have been separated for           )
% ( algorithmic clarity.                                                 )
%
% Diagnostic output
% ----------------------------------------------------------------------
% Diagnostic output consists of two sections:
%
% The first page lists the filenames, various parameters (Grand mean
% scaling etc.), and gives a plot of the image global means against
% scan number, overlaid on an "image" of the condition effects. Watch
% out for condition dependent global changes!
%
% The second part is a single page depicting the design matrix, effect
% names, parameter contrasts used, and the corresponding image files
% written.
%
% As always, look at the resulting mean images to make sure they look OK!
%
%
% Algorithm
% ----------------------------------------------------------------------
% The model at each voxel is Y = X*B + e, with a set of least squares
% estimates for the vector of parameters B as b = pinv(X)*Y. For c a
% vector of contrast weights extracting the appropriate parameter, the
% contrast of the parameter estimates is c'*b = c'*pinv(X)*Y, a
% weighted sum (or weighted mean) of the data at that voxel. These
% weights are identical for all voxels, so the image of the parameter
% estimate can be computed as a weighted mean of the images.
%
% The design matrix is split into effects of interest [C], a constant
% term [B], and confounds [G]. The columns of G are centered so that
% the confound cannot model any of the mean. The effects of interest
% are orthogonalised wirit. the confounds, using C=C-G*pinv(G)*C; This
% ensures that *all* confound effects are removed from the data, even
% if they are correlated with the effects of interest.
%
% Once the weights have been worked out for each adjusted mean image,
% computation proceeds by passing appropriate weights and image
% filenames to spm_mean, which writes out the appropriate parameter
% image as an Analyze format image of the same type (see spm_type) as
% the input images.
%
%
% Variables saved in SPMadj.mat data file
% ----------------------------------------------------------------------
% Des           Structure containing design parameters & specification
%   .DesName    Design name
%   .HForm      Form of DesMtx H partition
%   .iSubj      Subject indicator vector
%   .iCond      Condition indicator vector
%   .iRepl      Replication indicator vector
%   .iGloNorm   Global normalisation option
%   .sGloNorm   Global normalisation description
%   .iGMsca     Grand mean scaling option
%   .sGMsca     Grand mean scaling description
%   .GM         Grand Mean used for scaling
%   .iAdjTo     Adjustment (for AnCova) option
%   .sAdjTo     Adjustment (for AnCova) description
%   .aGM        AnCova adjustment value (subtracted from GX before AnCova)
%   .gSF        Image scale factors for global scaling
%   .X          Design matrix
%   .nX         Normalised (for imaging) design matrix
%   .Xnames     Effects corresponding to cols of X (cellstr)
%   .aPMap      Additional parameter to effect name mappings (see spm_desMtx)
%   .EXnames    English effect names corresponding to TeX parameters of Xnames 
%   .iX         Structure defining design matrix subpartitions
%       .H      Columns of X corresponding to H partition
%       .C      Columns of X corresponding to C partition
%       .B      Columns of X corresponding to B partition
%       .G      Columns of X corresponding to G partition
% c             Matrix of contrasts, contrasts in rows
% cNames        Names associated with contrasts
% W             Weights for images corresponding to contrasts
% Fnames        Filenames of adjusted mean images written (cellstr)
% rGX           raw global means (before any scaling)
% GX            Global means after scaling
%
% P     String matrix of filenames
% iCond     Condition indicator vector
% iGloNorm  Global normalisation option
% sGloNorm  Global normalisation description
% iGMsca    Grand mean scaling option
% sGMsca    Grand mean scaling description
% HPFc      High pass filter cut-off period (s)
% HPF       High pass filter
% sHPF      Description of high-pass filter
% rC        raw C partition of design matrix, prior to orthogonalisation
% C     C (covariates of interest) partition of design matrix
% Cnames    Names of parameters corresponding to columns of C
% B     B (block) partition of the design matrix
% Bnames    Names of parameters corresponding to columns of B
% G     G (confounding covariates) partition of design matrix
% Gnames    Names of parameters corresponding to columns of G
% rX        raw design matrix, prior to orthogonalisation of C partition
% X     design matrix (=[C,B,G])
% nrX       raw design matrix, normalised for display
% nX        design matrix, normalised for display
% c     Matrix of contrasts, contrasts in rows
% cNames    Names associated with contrasts
% W     Weights for images corresponding to contrasts
% CWD       Current Working Directory (when run)
% Fnames    Filenames of adjusted mean images written
% rGX       raw global means (before any scaling)
% gSF       Image scale factors for global scaling (inc. grand mean scaling)
% GX        Global means after scaling
% GM        Grans Mean used for scaling
%
%
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Andrew Holmes
% $Id: spm_adjmean_fmri_ui.m 4418 2011-08-03 12:00:13Z guillaume $



%=======================================================================
% - S E T U P
%=======================================================================
SCCSid = '2.9';
SPMid = spm('FnBanner',mfilename,SCCSid);
[Finter,Fgraph,CmdLine] = spm('FnUIsetup','AdjMean/fMRI',1);
spm_help('!ContextHelp',[mfilename,'.m'])


%=======================================================================
% - D E S I G N   P A R A M E T E R S
%=======================================================================

%-Get filenames
%-----------------------------------------------------------------------
P     = spm_select(Inf,'image','select scans (single session)');
nScan = size(P,1);
if nScan==1, error('Only one image - gimme more!'), end

%-Get condition information & construct vector of conditions indices
%-----------------------------------------------------------------------
%-Epoch condition order
eCond = spm_input('epoch order Eg 121... {0=null}',1,'c')';
Conds = unique(eCond);

%-Epoch length(s)
eLen   = [];
guiPos = spm_input('!NextPos');
while any(~eLen) | sum(eLen) ~= nScan
    eLen = spm_input('#scans/epoch Eg 8 or 8 6 8... ',guiPos);
    if length(eLen)==1, eLen=eLen*ones(1,length(eCond)); end
end
eLen    = eLen(:)';

%-Epoch onsets, in scans, starting at 1
eOns = cumsum(eLen) - eLen + 1;

%-iCond condition indicator vector
iCond = zeros(nScan,1);
for i=1:length(eLen), iCond(eOns(i)+[1:eLen(i)]-1)=eCond(i); end


%-Get Repeat time, construct approximate haemodynamic response function
%-----------------------------------------------------------------------
RT  = spm_input('interscan interval (in seconds)','+1');
hrf = spm_hrf(RT);


%-Construct design matrix C partition
% (box-cars convolved with approximate HRF if requested)
% null condition with iCond==0 is ignored by spm_DesMtx with '~' constraint
%-----------------------------------------------------------------------
[C,Cnames,Ci] = spm_DesMtx(iCond,'~','\alpha');
% C = spm_detrend(C);

%-Convolve with hemodynamic response function, if requested
bHRFconv = spm_input('Conv. box-cars w/ approx HRF?','+1','y/n',[1,0],1);
if bHRFconv
    d = length(hrf);
    C = [ones(d,1)*C(1,:); C];
    C = spm_sptop(hrf,nScan + d,1)*C;
    C = C([1:nScan]+d,:);
end


%-Construct design matrix block (B) partition
%-----------------------------------------------------------------------
B = ones(nScan,1); Bnames = '\mu';


%-Construct design matrix confound (G) partition?
%-----------------------------------------------------------------------
%*** Allow user specified confounds?
%*** ...and user specified convolution of confounds with approximate HRF?
G = []; Gnames = '';


%-Add high-pass filter using discrete cosine set
%-----------------------------------------------------------------------
bHPF = spm_input('Use high pass filter?','+1','y/n',[1 0],1);
if bHPF
    HPFc=nScan;
    for i=Conds', HPFc = min([HPFc,max(diff(eOns(eCond==i)))]); end
    HPFc=spm_input('HPF cut-off period {in seconds}','0','e',2*HPFc*RT);
    %-Find max order for discrete cosine set, HPFk
    % (from period>HPFc, period=nScan*RT/(k/2) for order k; k<nScan/2)
    HPFk = fix(min(2*(nScan*RT)/HPFc,nScan/2));
    HPF = []; HPFnames = '';
    for k = 1:HPFk
        HPF = [HPF, cos(k*pi*([1:nScan]-1)'/(nScan-1))];
        HPFnames = ...
            strvcat(HPFnames,sprintf('LowHz (%4.0fs)',2*RT*nScan/k));
    end
    HPF = HPF - ones(nScan,1)*mean(HPF);    %-Mean correct (by column)
    sHPF = sprintf('High-pass filter of %d components, cut off = %ds',...
        HPFk,HPFc);

    G = [G, HPF]; Gnames = strvcat(Gnames,HPFnames);
else
    HPFc=0; HPF=[]; sHPF='No high-pass filter';
end

%-Global normalization options
%-----------------------------------------------------------------------
sGloNorm = strvcat('None','Proportional scaling','AnCova');
iGloNorm = spm_input('Select global normalisation',...
    '+1','m',sGloNorm,[],2);
sGloNorm = deblank(sGloNorm(iGloNorm,:));


%-Grand mean scaling
%-----------------------------------------------------------------------
sGMsca = strvcat('None','Scaling of overall Grand Mean',...
    '(Implicitly via PropSca global normalisation)');
if iGloNorm==2, iGMsca=3; else, iGMsca=2; end
sGMsca = deblank(sGMsca(iGMsca,:));
GM = 100;


%-Temporal smoothing
%-----------------------------------------------------------------------
%**** What to do with temporal smoothing?
%SIGMA = spm_input('Temporal smoothing FWHM {secs}','+1','e',6);
%SIGMA = SIGMA/sqrt(8*log(2))/RT;



%=======================================================================
% - C O M P U T A T I O N
%=======================================================================
spm('FigName','AdjMean/fMRI: configuring',Finter,CmdLine);
fprintf('\tconfiguring: ')
spm('Pointer','Watch');

%-Memory map files
%-----------------------------------------------------------------------
V = spm_vol(char(P));
spm_check_orientations(V);

%-Work out required Analyze header info from handles
%-----------------------------------------------------------------------
DIM    = V(1).dim(1:3);
VOX    = sqrt(sum(V(1).mat(1:3,1:3).^2));
ORIGIN = (V(1).mat\[0 0 0 1]')';
ORIGIN = round(ORIGIN(1:3));

%-Compute global values
%-----------------------------------------------------------------------
fprintf('(globals)')
GX     = zeros(nScan,1);
for i  = 1:nScan, GX(i) = spm_global(V(i)); end
fprintf('\b - done)\n')

%-Scaling: compute global scaling factors required to implement proportional
% scaling global normalisation or Grand mean scaling, as requested
%-----------------------------------------------------------------------
rGX = GX;
if iGloNorm==2
    %-Proportional scaling global normalisation
    gSF = GM./GX;
    GX  = GM*ones(nScan,1);
    sGloNorm = sprintf('%s, to %g',sGloNorm,GM);
elseif iGMsca==2
    %-Grand mean scaling (overall)
    gSF = GM/mean(GX);
    GX  = GX*gSF;
    sGMsca = sprintf('%s, to %g',sGMsca,GM);
else    %-No scaling
    gSF = ones(nScan,1);
end


%-AnCova options: Construct Global covariates of no interest partition
%-----------------------------------------------------------------------
if iGloNorm==3
    G      = [G, GX-mean(GX)];
    Gnames = strvcat(Gnames,'Global');
    sGloNorm = [sGloNorm,', to Grand Mean'];
end

%-Design matrix - raw & with effects of interest orthogonalised wirit G
%-----------------------------------------------------------------------
rC = C;
rX = [rC B G];
[nrX,Xnames] = spm_DesMtx('sca',C,Cnames,B,Bnames,G,Gnames);
if size(G,2)
    C  = C - G*pinv(G)*C;
end
X  = [C B G];
nX = spm_DesMtx('sca',X,Xnames);
EXnames = spm_DesMtx('ETeXNames',Xnames);
iX = struct(    'H',    [],     'C',    [1:size(C,2)],...
        'B',    size(C,2)+1,    'G',    size(C,2) +1 +[1:size(G,2)]);

%-Temporal smoothing?
%-----------------------------------------------------------------------
%**** Include temporal smoothing? Theoretically, shouldn't make any
% difference to parameter estimates, but would improve efficiency and
% match SPM-fMRI.
%**** Q:Orthogonalisation of C wirit G pre or post temporal smoothing?

%-Parameter estimation matrix & data weights matrix for contrasts
%-----------------------------------------------------------------------
pX  = pinv(X);

%-Contrasts (c) & associated weight matrix (W)
%-----------------------------------------------------------------------
nc     = size(C,2);
c      = [eye(nc), ones(nc,size(B,2)), zeros(nc,size(G,2))];
cNames = spm_DesMtx('Fnames',EXnames(iX.C));
%-if we have a null condition, then constant *is* null condition effect
if any(Conds==0)
    c      = [zeros(1,nc), 1, zeros(1,size(G,2)); c];
    cNames = [{'baseline'};cNames];
    nc     = nc+1;
end

W   = c * pX;
Fnames = cNames;


%=======================================================================
% - D I A G N O S T I C   O U T P U T
%=======================================================================
fprintf('\tdisplaying diagnostic output...\n')
figure(Fgraph);


%-Display files and variables
%=======================================================================

%-Display parameters
%-----------------------------------------------------------------------
figure(Fgraph); spm_clf; axis off
text(0.30,1.00,'Adjusted means (fMRI)','Fontsize',16,'Fontweight','Bold')
text(-.10,0.96,'Files:','Fontweight','Bold')
text(-.10,0.94,sprintf('1 to %4d: %s...',nScan,P{1}),'FontSize',10)
text(-.10,0.92,sprintf('Repeat time = %.4f seconds',RT),'FontSize',10)
text(-.10,0.88,'Image Global means (& conditions):','Fontweight','Bold')

text(-.1,0.20,['Grand mean scaling: ',sGMsca])
text(-.1,0.17,['Global normalisation: ',sGloNorm])
text(-.1,0.14,sHPF)
text(-.1,0.11,'Effects of interest orthogonalised (residualised) to confounds')
text(-.1,0.08,'No temporal smoothing (only required for inference)')
text(-.1,0.05,sprintf('Parameters saved to: %s/SPMadj.mat',pwd),'FontSize',10)

%-Global mean values plot, with underlayed conditions
%-----------------------------------------------------------------------
%-Conditions plot (for comparison with global means plot)
axes('Position',[0.1,0.5,0.8,0.3])
image(32+(rC'-min(rC(:)))*32*(max(rC(:))-min(rC(:))))
set(gca,'Visible','off')

axes('Position',[0.1,0.5,0.8,0.3])
if iGloNorm==2
    plot(rGX)
    ylabel('global mean (pre-scaling)')
else
    plot(GX)
    ylabel('global mean (post-scaling)')
end
set(gca,'XLim',[0.5,nScan+0.5],'Color','none')
xlabel('scan index')

spm_print


%-Depict and label design matrix, depict & label contrasts
%=======================================================================
spm_clf(Fgraph); axes('Position',[0 0 1 .95],'Visible','off')
text(.2,1,'Design Matrix, contrasts & contrast files',...
    'Fontsize',16,'Fontweight','Bold');

%-Image scaled design matrix, label the effects, add filenames (if <=32)
%-----------------------------------------------------------------------
hDesMtx = axes('Position',[.07 .5 .6 .3]);
image((nX + 1)*32);
set(hDesMtx,'TickDir','out')
ylabel('scans')
xlabel('parameters')
axes('Position',[.07 .8 .6 .1],'Visible','off')
dx    = 1/size(nX,2);
tXnames = spm_DesMtx('TeXnames',Xnames);
for i = 1:size(nX,2)
    text((i - .5)*dx+.01,.3,EXnames{i},'Fontsize',9,'Rotation',90)
    if ~strcmp(Xnames{i},EXnames{i})
        text((i - .5)*dx,.1,tXnames{i},'Fontsize',10,'Interpreter','TeX')
    end
end
if nScan<=32
    set(hDesMtx,'YTick',1:nScan)
    axes('Position',[.68 .5 .3 .3],'Visible','off')
    dy = 1/nScan;
    tP = spm_str_manip(P,'k40');
    for i=1:nScan
        text(0,(nScan-i+0.5)*dy,tP{i},'FontSize',8)
    end
end

%-Depict contrasts and associated (preliminary) filenames
%-----------------------------------------------------------------------
dy = .4/nc;
axes('Position',[.025 .05 .05 .4],'Visible','off')
text(0,.5,'contrasts','HorizontalAlignment','Center','Rotation',90,...
    'FontSize',14,'FontWeight','Bold')
axes('Position',[.6 .44 .40 .02],'Visible','off')
text(0,1,'Contrast files...','FontSize',10,'FontWeight','Bold')
text(0,0,sprintf('...in %s',pwd),'FontSize',8)
hFnames = zeros(1,nc);
for i = 1:nc
    axes('Position',[.1 (.45 -dy*i) .6 .9*dy])
    h = bar(c(i,:),1);
    set(h,'FaceColor',[1 1 1]*.8)
    set(gca,'XLim',[.5,size(c,2)+.5],'Visible','off')
    text(0,0,num2str(i),'HorizontalAlignment','Right','FontSize',10)
    hFnames(i) = text(size(c,2)+.55,.1,Fnames{i},'FontSize',10,...
        'Color',[1 1 1]*.5,'FontAngle','Italic');
end



%=======================================================================
% - W R I T E   C O N T R A S T   I M A G E S
%=======================================================================
fprintf('\tspecify filenames for contrast images...\n')
spm('Pointer');
guiPos = spm_input('!NextPos');
for i = 1:nc
    Fnames{i} = spm_input(sprintf('Contrast %d: filename?',i),...
            guiPos,'s',Fnames{i});
    set(hFnames(i),'string',Fnames{i},'Color','k','FontAngle','Normal')
end
spm_print

%-Parameter images (of interest) - Adjusted mean images
%-----------------------------------------------------------------------
spm('FigName','AdjMean/fMRI - writing',Finter,CmdLine);
spm('Pointer','Watch');

%-Computation - calculations handled by spm_add.c
%-Using implicit zero masking feature of spm_add: The resultant image
% will be zero at voxels where *any* of the input images are zero.

%-Create handle template for output images as 16bit Int16's
Vo = struct(    'fname',    '',...
        'dim',      V(1).dim(1:3),...
        'dt',           [4 spm_platform('bigend')],...
        'mat',      V(1).mat,...
        'pinfo' ,   [1,0,0]',...
        'descrip',  '');

%-Loop over contrasts
for i = 1:nc
    fprintf('\t...writing image %d/%d: %-20s',i,nc,Fnames{i})
    %-Implement weighted sum by weighting scalefactors in image handles
    w  = W(i,:).*gSF';
    Q  = find(abs(w)>0);
    w  = w(Q); wV = V(Q);
    for j=1:length(Q), wV(j).pinfo(1:2,:)=wV(j).pinfo(1:2,:)*w(j); end
    %-Write header
    Vo.fname   = [Fnames{i},'.img'];
    Vo.descrip = sprintf('Adjusted mean (spm_adjmean_fmri) - %s',Fnames{i});
    Vo.pinfo   = [1,0,0]';
    Vo         = spm_create_vol(Vo);
    %-Compute & rewrite header scalefactor
    Vo.pinfo(1)  = spm_add(wV,Vo,'m');
    Vo           = spm_create_vol(Vo);

    fprintf(' (done)\n')
end

%-Prepend PWD to Fnames
Fnames = cellstr([repmat([pwd,filesep],nc,1),char(Fnames)]);


%-Save parameters to SPMadj.mat in current directory
%-----------------------------------------------------------------------
%-Pack design parameters up in a structure
Des = struct(...
        'DesName',  'Single sess. fMRI box-cars',...
        'HForm',    [],...
        'iSubj',    [],...
        'iCond',    iCond,...
        'iGloNorm', iGloNorm,...
        'sGloNorm', sGloNorm,...
        'iGMsca',   iGMsca,...
        'sGMsca',   sGMsca,...
        'GM',       GM,...
        'gSF',      gSF,...
        'iAdjTo',   [],...
        'sAdjTo',   [],...
        'aGM',      [],...
        'X',        X,...
        'nX',       nX,...
        'Xnames',   Xnames,...
        'aPMap',    [],...
        'EXnames',  EXnames,...
        'iX',       iX      );
if spm_check_version('matlab','7') >= 0
    save('SPMadj.mat','-V6',...
        'SPMid','Des','V','c','cNames','W','Fnames','HPFc','HPF','sHPF','rX','nrX','rGX','GX');
else
    save('SPMadj.mat',...
        'SPMid','Des','V','c','cNames','W','Fnames','HPFc','HPF','sHPF','rX','nrX','rGX','GX');
end;


%=======================================================================
% - E N D
%=======================================================================
spm('FigName','AdjMean/fMRI - done',Finter,CmdLine);
spm('Pointer','Arrow')
fprintf('\n\n')