function [normalise] = volumenormalise(cfg, interp)

% VOLUMENORMALISE normalises anatomical and functional volume data
% to a template anatomical MRI.
%
% This function uses the SPM2-toolbox, see http://www.fil.ion.ucl.ac.uk/spm/
%
% Use as
%   [volume] = volumenormalise(cfg, volume)
%
% The input volume should be the result from SOURCEINTERPOLATE.
% Alternatively, the input can contain a single anatomical MRI that
% was read with READ_FCDC_MRI, or you can specify a filename of an
% anatomical MRI.
%
% Configuration options are:
%   cfg.template    = '/opt/spm2/templates/T1.mnc' or string
%   cfg.parameter   = cell-array with the functional data which has to
%                     be normalised, can be 'all'
%   cfg.downsample  = number (default = 1, meaning no downsampling)
%   cfg.name        = string for output filename
%   cfg.write       = 'no' (default) or 'yes', writes the segmented volumes to SPM2
%                     compatible analyze-file, with the suffix
%                     _anatomy for the anatomical MRI volume
%                     _param   for each of the functional volumes
%   cfg.nonlinear   = 'yes' (default) or 'no', estimates a nonlinear transformation
%                     in addition to the linear affine registration. If a reasonably
%                     accurate normalisation is sufficient, a purely linearly transformed
%                     image allows for 'reverse-normalisation', which might come in handy
%                     when for example a region of interest is defined on the normalised
%                     group-average.                     

% Undocumented options
%   cfg.keepintermediate = 'yes' or 'no'
%   cfg.intermediatename = prefix of the the coregistered images and of the
%                          original images in the original headcoordinate system
%   cfg.coordinates      = 'ctf' or 'spm' (prevent interactive questions)
%   cfg.spmparams        = one can feed in parameters from a prior normalisation

% Copyright (C) 2004, Jan-Mathijs Schoffelen
%
% $Log: volumenormalise.m,v $
% Revision 1.8  2006/02/24 16:46:23  roboos
% switched to fixvolume function, with inside as 3d volume
% not neccessary any more to treat inside/outside seperately
% not neccessary any more to reshape volumes
% not neccessary any more to change from xgrid/etc into transform
%
% Revision 1.7  2006/01/31 10:06:23  jansch
% changed the order of volumedownsample and the conversion of inside into a
% volume
%
% Revision 1.6  2006/01/31 09:48:39  jansch
% fixed bug regarding last change
%
% Revision 1.5  2006/01/31 09:45:32  jansch
% changed else cfg.keepinside='no' into an elseif to prevent crashes
%
% Revision 1.4  2006/01/31 09:27:58  roboos
% in case of keepinside: instead of setting outside to [], remove it from structure
%
% Revision 1.3  2006/01/30 13:55:19  roboos
% changed the order in which the initial preparatory steps are taken
% related to inside-conversion, parameterselection and grid2transform
%
% Revision 1.2  2006/01/24 14:20:35  roboos
% removed the obsolete option cfg.voxelcoord, new behaviour is that it is always 'yes'
%
% Revision 1.1  2006/01/05 12:58:20  roboos
% This function (VOLUMExxx) replaces a function with the name xxxVOLUME.
% The fields xgrid/ygrid/zgrid are removed (this is from now on handled by
% grid2transform and the VOLUMExxx function should only work on volumes that
% are described using a transformation matrix).
% Writing of spm/analyze volumes is handled by private/volumewrite_spm.
%
% Revision 1.12  2005/11/10 17:20:45  roboos
% switched to using the general parameterselection subfunction
% changed the order of handling the downsampling and the inside handling
%
% Revision 1.11  2005/11/10 10:02:14  roboos
% minor change in whitespace and help
%
% Revision 1.10  2005/10/20 12:44:45  jansch
% included the option to do only linear normalisation
%
% Revision 1.9  2005/09/12 09:44:49  jansch
% also normalise 'ref'
%
% Revision 1.8  2005/08/19 16:55:17  roboos
% added the option to re-use the warping parameters that were computed in a previous run (cfg.spmparams)
%
% Revision 1.7  2005/06/29 12:46:29  roboos
% the following changes were done on a large number of files at the same time, not all of them may apply to this specific file
% - added try-catch around the inclusion of input data configuration
% - moved cfg.version, cfg.previous and the assignment of the output cfg to the end
% - changed help comments around the configuration handling
% - some changes in whitespace
%
% Revision 1.6  2005/05/17 17:50:37  roboos
% changed all "if" occurences of & and | into && and ||
% this makes the code more compatible with Octave and also seems to be in closer correspondence with Matlab documentation on shortcircuited evaluation of sequential boolean constructs
%
% Revision 1.5  2005/05/10 09:16:46  roboos
% fixed inside/outside brain assignment for voxels that could not be interpolated by SPM (nan values)
%
% Revision 1.4  2005/05/06 10:34:52  roboos
% fixed bug in assignment of voxels that were outside of the scanned volume
%
% Revision 1.3  2005/05/04 16:33:29  roboos
% added cfg.coordinates=ctf|spm to the "undocumented options" documentation
% changed the order of adding the parameters to the list in order to fix problem with inside volume
%
% Revision 1.2  2005/03/14 10:10:18  jansch
% fixed small bug in the handling of file-names for functional images
%
% Revision 1.1  2005/03/03 10:44:38  roboos
% Cleaned up the sourcenormalise function and moved everything to a
% new function with the name normalisevolume.  The new function has a
% different way of handling the subfields that can be normalised, and uses
% a different (and more robust) initial pre-alignment of the input volume
% with the SPM template. Also the type-detection (i.e. either SPM or CTF)
% of the input volume has been changed.
%
% Revision 1.6  2005/02/09 13:08:55  roboos
% changed the order in which the parameters were assigned, there was a probloem if cfg.parameter was not 'all'
%
% Revision 1.5  2005/02/08 14:04:45  roboos
% implemented keepinside, changed handling of parameters to use getsubfield/issubfield
% separated default transformation as a separate variable (no functional change yet)
%
% Revision 1.4  2004/09/10 11:47:01  jansch
% fixed bug in checking for the presence of specified parameters in
% configuration
%
% Revision 1.3  2004/09/09 05:48:13  jansch
% removed bug using an ambiguous function for the rotation of the transform
% matrix.
%
% Revision 1.2  2004/09/08 11:15:04  jansch
% added version-information to output.
%

% check if spm2 is in your path, see http://www.fil.ion.ucl.ac.uk/spm
hasspm = (exist('spm_vol') & exist('spm_write_vol') & exist('spm_normalise'));
if ~hasspm
  error('the SPM2 toolbox is not installed, see: http://www.fil.ion.ucl.ac.uk/spm/');
end

% set the defaults
if ~isfield(cfg,'template'),         cfg.template = '/opt/spm2/templates/T1.mnc'; end;
if ~isfield(cfg,'parameter'),        cfg.parameter = 'all';                       end;
if ~isfield(cfg,'downsample'),       cfg.downsample = 1;                          end;
if ~isfield(cfg,'write'),            cfg.write = 'no';                            end;
if ~isfield(cfg,'keepinside'),       cfg.keepinside = 'yes';                      end;
if ~isfield(cfg,'keepintermediate'), cfg.keepintermediate = 'no';                 end;
if ~isfield(cfg,'coordinates'),      cfg.coordinates = [];                        end;
if ~isfield(cfg,'nonlinear'),        cfg.nonlinear = 'yes';                       end;

if strcmp(cfg.keepinside, 'yes')
  % add inside to the list of parameters
  if ~iscell(cfg.parameter),
    cfg.parameter = {cfg.parameter 'inside'};
  else
    cfg.parameter(end+1) = {'inside'};
  end
end

if ~isfield(cfg,'intermediatename')
  cfg.intermediatename = tempname;
end

if ~isfield(cfg,'name') && strcmp(cfg.write,'yes')
  error('you must specify the output filename in cfg.name');
end

if isempty(cfg.template),
  error('you must specify a template anatomical MRI');
end

if ~isfield(interp,'anatomy'),
  error('no anatomical information available, this is required for normalisation');
end

if ischar(interp),
  % read the anatomical MRI data from file
  filename = interp;
  fprintf('reading MRI from file\n');
  interp = read_fcdc_mri(filename);
  if filetype(filename, 'ctf_mri') && isempty(cfg.coordinates)
    % based on the filetype assume that the coordinates correspond with CTF convention
    cfg.coordinates = 'ctf';
  end
end

% ensure that the structure correctly describes a volume
interp = fixvolume(interp);

% select the parameters that should be normalised
cfg.parameter = parameterselection(cfg.parameter, interp);

% the anatomy should always be normalised as the first volume
cfg.parameter(strcmp(cfg.parameter, 'anatomy')) = [];
cfg.parameter = {'anatomy' cfg.parameter{:}};

% downsample the volume
tmpcfg = [];
tmpcfg.downsample = cfg.downsample;
tmpcfg.parameter  = cfg.parameter;
interp = volumedownsample(tmpcfg, interp);

% remember the original transformation matrix
original.transform = interp.transform;
original.nosedir = [];
original.origin = [];

if isempty(cfg.coordinates)
  % determine in which direction the nose is pointing
  while isempty(original.nosedir)
    response = input('Is the nose pointing in the positive X- or Y-direction? [x/y] ','s');
    if strcmp(response, 'x'),
      original.nosedir = 'positivex';
    elseif strcmp(response, 'y'),
      original.nosedir = 'positivey';
    end
  end

  % determine where the origin is
  while isempty(original.origin)
    response = input('Is the origin on the Anterior commissure or between the Ears? [a/e] ','s');
    if strcmp(response, 'e'),
      original.origin = 'interauricular';
    elseif strcmp(response, 'a'),
      original.origin = 'ant_comm';
    end
  end

  % determine the coordinatesystem of the input MRI volume
  if strcmp(original.origin, 'interauricular') && strcmp(original.nosedir, 'positivex')
    cfg.coordinates = 'ctf';
  elseif strcmp(original.origin, 'ant_comm') && strcmp(original.nosedir, 'positivey')
    cfg.coordinates = 'spm';
  end
end

% ensure that the input MRI volume is approximately aligned with the SPM template
if strcmp(cfg.coordinates, 'ctf')
  fprintf('assuming CTF coordinates for input, i.e. positive X-axis towards nasion and Y-axis through ears\n');
  % flip, rotate and translate the CTF mri so that it approximately corresponds with SPM coordinates
  interp = align_ctf2spm(interp);
elseif strcmp(cfg.coordinates, 'spm')
  fprintf('assuming that the input MRI is already approximately alligned with SPM coordinates\n');
  % nothing needs to be done
else
  error('cannot determine the (approximate) alignmenmt of the input MRI with the SPM template');
end

warning off;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% here the normalisation starts
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% create an spm-compatible header for the anatomical volume data
Va = volumewrite_spm([cfg.intermediatename,'_anatomy.img'], interp.anatomy, interp.transform);

% create an spm-compatible file for each of the functional volumes
for parlop=2:length(cfg.parameter)  % skip the anatomy
  tmp  = cfg.parameter{parlop};
  data = reshape(getsubfield(interp, tmp), interp.dim);
  tmp(find(tmp=='.')) = '_';
  volumewrite_spm([cfg.intermediatename,'_' tmp '.img'], data, interp.transform);
end

% read the template anatomical volume data
if strcmp(char(cfg.template(end-2:end)),'mnc'),
  Vtempl    = spm_vol_minc(cfg.template);
elseif strcmp(char(cfg.template(end-2:end)),'img'),
  Vtempl    = spm_vol(cfg.template);
else
  error('Unknown template');
end

fprintf('performing the normalisation\n');
% do spatial normalisation according to these steps
% step 1: read header information for template and source image
% step 2: compute transformation parameters
% step 3: write the results to a file with prefix 'w'

if ~isfield(cfg, 'spmparams') && strcmp(cfg.nonlinear, 'yes'),
  fprintf('warping the invdividual anatomy to the template anatomy\n');
  % compute the parameters by warping the individual anatomy
  Va        = spm_vol([cfg.intermediatename,'_anatomy.img']);
  params    = spm_normalise(Vtempl,Va);
elseif ~isfield(cfg, 'spmparams') && strcmp(cfg.nonlinear, 'no'),
  fprintf('warping the invdividual anatomy to the template anatomy, using only linear transformations\n');
  % compute the parameters by warping the individual anatomy
  Va         = spm_vol([cfg.intermediatename,'_anatomy.img']);
  flags.nits = 0; %put number of non-linear iterations to zero
  params     = spm_normalise(Vtempl,Va,[],[],[],flags);
else
  fprintf('using the parameters specified in the configuration\n');
  % use the externally specified parameters
  params = cfg.spmparams;
end
flags.vox = [cfg.downsample,cfg.downsample,cfg.downsample];
files     = {};

% remember the normalisation parameters in the configuration
cfg.spmparams = params;

% apply the normalisation parameters to each of the volumes
for parlop=1:length(cfg.parameter)
  fprintf('creating normalised analyze-file for %s\n', cfg.parameter{parlop});
  tmp = cfg.parameter{parlop};
  tmp(find(tmp=='.')) = '_';
  files{parlop} = sprintf('%s_%s.img', cfg.intermediatename, tmp);
  [p, f, x] = fileparts(files{parlop});
  wfiles{parlop} = fullfile(p, ['w' f x]);
end
spm_write_sn(char(files),params,flags);  % this creates the 'w' prefixed files

normalise = [];

% read the normalised results from the 'w' prefixed files
V = spm_vol(char(wfiles));
for vlop=1:length(V)
  normalise = setsubfield(normalise, cfg.parameter{vlop}, spm_read_vols(V(vlop)));
end

normalise.transform = V(1).mat;
normalise.dim       = size(normalise.anatomy);

if isfield(normalise, 'inside')
  % convert back to a logical volume
  normalise.inside  = abs(normalise.inside-1)<=10*eps;
end

% flip and permute the dimensions to align the volume with the headcoordinate axes
normalise = align_ijk2xyz(normalise);

if strcmp(cfg.write,'yes')
  % create an spm-compatible file for each of the normalised volumes
  for parlop=1:length(cfg.parameter)  % include the anatomy
    tmp  = cfg.parameter{parlop};
    data = reshape(getsubfield(normalise, tmp), normalise.dim);
    tmp(find(tmp=='.')) = '_';
    volumewrite_spm([cfg.name,'_' tmp '.img'], data, normalise.transform);
  end
end

if strcmp(cfg.keepintermediate,'no')
  % remove the intermediate files
  for flop=1:length(files)
    [p, f, x] = fileparts(files{flop});
    delete(fullfile(p, [f, '.*']));
    [p, f, x] = fileparts(wfiles{flop});
    delete(fullfile(p, [f, '.*']));
  end
end

% add version information to the configuration
try
  % get the full name of the function
  cfg.version.name = mfilename('fullpath');
catch
  % required for compatibility with Matlab versions prior to release 13 (6.5)
  [st, i] = dbstack;
  cfg.version.name = st(i);
end
cfg.version.id = '$Id: volumenormalise.m,v 1.8 2006/02/24 16:46:23 roboos Exp $';
% remember the configuration details of the input data
try, cfg.previous = interp.cfg; end
% remember the exact configuration details in the output 
normalise.cfg = cfg;

