function [interp] = ft_sourceinterpolate(cfg, functional, anatomical)

% FT_SOURCEINTERPOLATE interpolates source activity or statistical maps onto the
% voxels or vertices of an anatomical description of the brain.  Both the functional
% and the anatomical data can either describe a volumetric 3D regular grid, a
% triangulated description of the cortical sheet or an arbitrary cloud of points.
%
% The functional data in the output data will be interpolated at the locations at
% which the anatomical data are defined. For example, if the anatomical data was
% volumetric, the output data is a volume-structure, containing the resliced source
% and the anatomical volume that can be visualized using FT_SOURCEPLOT or written to
% file using FT_SOURCEWRITE.
%
% The following scenarios can be considered:
%
% - Both functional data and anatomical data are defined on 3D regular grids, for
%   example with a low-res grid for the functional data and a high-res grid for the
%   anatomy.
%
% - The functional data is defined on a 3D regular grid and the anatomical data is
%   defined on an irregular point cloud, which can be a 2D triangulated surface mesh.
%
% - The functional data is defined on an irregular point cloud, which can be a 2D
%   triangulated surface mesh, and the anatomical data is defined on a 3D regular grid.
%
% - Both the functional and the anatomical data are defined on an irregular point
%   cloud, which can be a 2D triangulated mesh.
%
% - The functional data is defined on a low-resolution 2D triangulated surface mesh and the
%   anatomical data is defined on a high-resolution 2D triangulated surface mesh, where the
%   low-res vertices form a subset of the high-res vertices. This allows for mesh-based
%   interpolation. The algorithm currently implemented is so-called 'smudging' as it is
%   also applied by the MNE-suite software.
%
% Use as
%   [interp] = ft_sourceinterpolate(cfg, source, anatomy)
%   [interp] = ft_sourceinterpolate(cfg, stat,   anatomy)
% where
%   source  is the output of FT_SOURCEANALYSIS
%   stat    is the output of FT_SOURCESTATISTICS
%   anatomy is the output of FT_READ_MRI, or one of the FT_VOLUMExxx functions,
%           or a cortical sheet that was read with FT_READ_HEADSHAPE,
%           or a regular 3D grid created with FT_PREPARE_SOURCEMODEL.
%
% The configuration should contain:
%   cfg.parameter     = string or cell-array with the functional parameter(s) to be interpolated
%   cfg.downsample    = integer number (default = 1, i.e. no downsampling)
%   cfg.interpmethod  = string, can be 'nearest', 'linear', 'cubic',  'spline', 'sphere_avg', 'sphere_weighteddistance', or 'smudge' (default = 'linear for interpolating two 3D volumes, 'nearest' for all other cases)
%
% For interpolating two 3D regular grids or volumes onto each other the supported
% interpolation methods are 'nearest', 'linear', 'cubic' or 'spline'. For all other
% cases the supported interpolation methods are 'nearest', 'sphere_avg',
% 'sphere_weighteddistance' or 'smudge'.
%
% The functional and anatomical data should be expressed in the same
% coordinate sytem, i.e. either both in MEG headcoordinates (NAS/LPA/RPA)
% or both in SPM coordinates (AC/PC).
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
%   cfg.outputfile  =  ...
% If you specify one of these (or both) the input data will be read from a *.mat
% file on disk and/or the output data will be written to a *.mat file. These mat
% files should contain only a single variable, corresponding with the
% input/output structure.
%
% See also FT_READ_MRI, FT_READ_HEADSHAPE, FT_SOURCEPLOT, FT_SOURCEANALYSIS,
% FT_SOURCEWRITE

% Copyright (C) 2003-2007, Robert Oostenveld
% Copyright (C) 2011-2014, Jan-Mathijs Schoffelen
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar functional anatomical
ft_preamble provenance functional anatomical

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% this is not supported any more as of 26/10/2011
if ischar(anatomical)
  ft_error('please use cfg.inputfile instead of specifying the input variable as a sting');
end

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'unused',     {'keepinside' 'voxelcoord'});
cfg = ft_checkconfig(cfg, 'deprecated', {'sourceunits', 'mriunits'});
cfg = ft_checkconfig(cfg, 'required',   'parameter');
cfg = ft_checkconfig(cfg, 'renamedval', {'parameter', 'avg.pow', 'pow'});
cfg = ft_checkconfig(cfg, 'renamedval', {'parameter', 'avg.coh', 'coh'});
cfg = ft_checkconfig(cfg, 'renamedval', {'parameter', 'avg.mom', 'mom'});

% set the defaults
cfg.downsample   = ft_getopt(cfg, 'downsample', 1);
cfg.feedback     = ft_getopt(cfg, 'feedback', 'text');
cfg.interpmethod = ft_getopt(cfg, 'interpmethod', []);   % cfg.interpmethod depends on how the interpolation should be done and actual defaults will be specified below

% replace pnt by pos
anatomical = fixpos(anatomical);
functional = fixpos(functional);

% ensure the functional data to be in double precision, the maxdepth parameter ensure double precision up to the content of functional.avg.mom{:}, avoiding too much recursion
functional = ft_struct2double(functional, 3);

if (strcmp(cfg.interpmethod, 'nearest') || strcmp(cfg.interpmethod, 'mode')) && (ft_datatype(functional, 'volume+label') || ft_datatype(functional, 'source+label') || ft_datatype(functional, 'mesh+label'))
  % the first input argument describes a parcellation or segmentation with tissue labels
  isAtlasFun = true;
else
  isAtlasFun = false;
end

if isfield(anatomical, 'transform') && isfield(anatomical, 'dim')
  % anatomical volume
  isUnstructuredAna  = false;
elseif isfield(anatomical, 'pos') && isfield(anatomical, 'dim')
  % positions that can be mapped onto a 3D regular grid
  isUnstructuredAna  = false;
elseif isfield(anatomical, 'pos')
  % anatomical data that consists of a mesh, but no smudging possible
  isUnstructuredAna  = true;
end

if isfield(functional, 'transform') && isfield(functional, 'dim')
  % functional volume
  isUnstructuredFun  = false;
elseif isfield(functional, 'pos') && isfield(functional, 'dim')
  % positions that can be mapped onto a 3D regular grid
  isUnstructuredFun  = false;
else
  isUnstructuredFun  = true;
end

if isUnstructuredAna
  anatomical = ft_checkdata(anatomical, 'datatype', {'source', 'source+label', 'mesh'}, 'insidestyle', 'logical', 'feedback', 'yes', 'hasunit', 'yes');
else
  anatomical = ft_checkdata(anatomical, 'datatype', {'volume', 'volume+label'}, 'insidestyle', 'logical', 'feedback', 'yes', 'hasunit', 'yes');
end

if isUnstructuredFun
  functional = ft_checkdata(functional, 'datatype', 'source', 'insidestyle', 'logical', 'feedback', 'yes', 'hasunit', 'yes');
else
  functional = ft_checkdata(functional, 'datatype', 'volume', 'insidestyle', 'logical', 'feedback', 'yes', 'hasunit', 'yes');
end

% select the parameters from the data, this needs to be done here, because after running checkdata, the parameterselection fails if the numeric data has nfreq/ntime/etc>1
cfg.parameter = parameterselection(cfg.parameter, functional);

% ensure that the functional data has the same unit as the anatomical data
functional = ft_convert_units(functional, anatomical.unit);

if isfield(functional, 'coordsys') && isfield(anatomical, 'coordsys') && ~isequal(functional.coordsys, anatomical.coordsys)
  % FIXME is this different when smudged or not?
  % ft_warning('the coordinate systems are not aligned');
  % ft_error('the coordinate systems are not aligned');
end

if ~isUnstructuredAna && cfg.downsample~=1
  % downsample the anatomical volume
  tmpcfg = keepfields(cfg, {'downsample', 'showcallinfo', 'trackcallinfo', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo', 'checksize'});
  tmpcfg.parameter = 'anatomy';
  anatomical = ft_volumedownsample(tmpcfg, anatomical);
  % restore the provenance information and put back cfg.parameter
  tmpparameter = cfg.parameter;
  [cfg, anatomical] = rollback_provenance(cfg, anatomical);
  cfg.parameter = tmpparameter;
end

% collect the functional volumes that should be converted
dat_name = {};
dat_array = {};
for i=1:length(cfg.parameter)
  if ~iscell(getsubfield(functional, cfg.parameter{i}))
    dat_name{end+1} = cfg.parameter{i};
    dat_array{end+1} = getsubfield(functional, cfg.parameter{i});
  else
    fprintf('not interpolating %s, since it is represented in a cell-array\n', cfg.parameter{i});
  end
end

% hmmmm, if the input data contains a time and/or freq dimension, then the output
% may be terribly blown up; most convenient would be to output only the
% smudging matrix, and project the data when plotting

if isUnstructuredFun && isUnstructuredAna && isfield(anatomical, 'orig') && isfield(anatomical.orig, 'pos') && isfield(anatomical.orig, 'tri')
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % functional data defined on subset of vertices in an anatomical mesh
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % FIXME this should not be decided on the basis of the data structures but on the basis of the cfg.interpmethod option
  % FIXME the distribution of 3 geometries over the 2 structures is weird
  % FIXME a (perhaps extreme) application of this would be to interpolate data from parcels on the sheet, i.e. an inverse parcellation

  % anatomical data consists of a decimated triangulated mesh, containing
  % the original description, allowing for smudging.

  % smudge the low resolution functional data according to the strategy in
  % MNE-suite (chapter 8.3 of the manual)

  interpmat = interp_ungridded(anatomical.pos, anatomical.orig.pos, 'projmethod', 'smudge', 'triout', anatomical.orig.tri);
  interpmat(~anatomical.inside(:), :) = 0;

  % start with an empty structure, keep only some fields
  interp = keepfields(functional, {'time', 'freq'});
  interp = copyfields(anatomical, interp, {'unit', 'coordsys'});
  interp = copyfields(anatomical.orig, interp, {'pos', 'tri', 'dim'});

  % identify the inside voxels after interpolation
  nzeros     = sum(interpmat~=0,2);
  newinside  = (nzeros~=0);
  newoutside = (nzeros==0);

  interp.inside = false(size(anatomical.pos,1),1);
  interp.inside(newinside) = true;

  % interpolate all functional data
  for i=1:length(dat_name)
    fprintf('interpolating %s\n', dat_name{i});

    dimord = getdimord(functional, dat_name{i});
    dimtok = tokenize(dimord, '_');
    dimf   = getdimsiz(functional, dat_name{i});
    dimf(end+1:length(dimtok)) = 1; % there can be additional trailing singleton dimensions

    % should be 3-D array, can have trailing singleton dimensions
    if numel(dimf)<2
      dimf(2) = 1;
    end
    if numel(dimf)<3
      dimf(3) = 1;
    end

    allav = zeros([size(anatomical.orig.pos,1), dimf(2:end)]);
    for k=1:dimf(2)
      for m=1:dimf(3)
        fv     = double_ifnot(dat_array{i}(:,k,m));
        av     = interpmat*fv;
        av(newoutside) = nan;
        allav(:,k,m)   = av;
      end
    end
    interp = setsubfield(interp, dat_name{i}, allav);
  end


elseif isUnstructuredFun && isUnstructuredAna
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % functional data defined on a point cloud/mesh, anatomy on a point cloud/mesh
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % set default interpmethod for this situation
  cfg.interpmethod = ft_getopt(cfg, 'interpmethod', 'nearest');
  cfg.sphereradius = ft_getopt(cfg, 'sphereradius', 0.5);
  cfg.power        = ft_getopt(cfg, 'power',        1);

  interpmat = interp_ungridded(functional.pos, anatomical.pos, 'projmethod', cfg.interpmethod, 'sphereradius', cfg.sphereradius, 'power', cfg.power); % FIXME include other key-value pairs as well
  interpmat(~anatomical.inside(:), :) = 0;

  % start with an empty structure, keep only some fields
  interp = keepfields(functional, {'time', 'freq'});
  interp = copyfields(anatomical, interp, {'pos', 'tri', 'dim', 'transform', 'unit', 'coordsys'});

  % identify the inside voxels after interpolation
  nzeros     = sum(interpmat~=0,2);
  newinside  = (nzeros~=0);
  newoutside = (nzeros==0);

  interp.inside = false(size(anatomical.pos,1),1);
  interp.inside(newinside) = true;

  % interpolate all functional data
  for i=1:length(dat_name)
    fprintf('interpolating %s\n', dat_name{i});

    dimord = getdimord(functional, dat_name{i});
    dimtok = tokenize(dimord, '_');
    dimf   = getdimsiz(functional, dat_name{i});
    dimf(end+1:length(dimtok)) = 1; % there can be additional trailing singleton dimensions

    % should be 3-D array, can have trailing singleton dimensions
    if numel(dimf)<2
      dimf(2) = 1;
    end
    if numel(dimf)<3
      dimf(3) = 1;
    end

    allav = zeros([size(anatomical.pos,1), dimf(2:end)]);
    for k=1:dimf(2)
      for m=1:dimf(3)
        fv     = double_ifnot(dat_array{i}(:,k,m));
        av     = interpmat*fv;
        av(newoutside) = nan;
        allav(:,k,m)   = av;
      end
    end
    interp = setsubfield(interp, dat_name{i}, allav);
  end


elseif isUnstructuredFun && ~isUnstructuredAna
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % functional data defined on a point cloud/mesh, anatomy on a volume
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
  % set default interpmethod for this situation
  cfg.interpmethod = ft_getopt(cfg, 'interpmethod', 'nearest');
  cfg.sphereradius = ft_getopt(cfg, 'sphereradius', 0.5);
  cfg.power        = ft_getopt(cfg, 'power',        1);

  [ax, ay, az] = voxelcoords(anatomical.dim, anatomical.transform);
  anatomical.pos = [ax(:) ay(:) az(:)];
  clear ax ay az

  tmp = interp_ungridded(functional.pos, anatomical.pos(anatomical.inside,:), 'projmethod', cfg.interpmethod, 'sphereradius', cfg.sphereradius, 'power', cfg.power); % FIXME include other key-value pairs as well
  interpmat( anatomical.inside(:), :) = tmp;
  interpmat(~anatomical.inside(:), :) = 0;
  
  % start with an empty structure, keep only some fields
  interp = keepfields(functional, {'time', 'freq'});
  interp = copyfields(anatomical, interp, {'pos', 'tri', 'dim', 'transform', 'unit', 'coordsys', 'anatomy'});

  % identify the inside voxels after interpolation
  nzeros     = sum(interpmat~=0,2);
  newinside  = (nzeros~=0);
  newoutside = (nzeros==0);

  interp.inside = false(anatomical.dim);
  interp.inside(newinside) = true;

  % interpolate all functional data
  for i=1:length(dat_name)
    fprintf('interpolating %s\n', dat_name{i});

    dimord = getdimord(functional, dat_name{i});
    dimtok = tokenize(dimord, '_');
    dimf   = getdimsiz(functional, dat_name{i});
    dimf(end+1:length(dimtok)) = 1; % there can be additional trailing singleton dimensions

    % should be 3-D array, can have trailing singleton dimensions
    if numel(dimf)<2
      dimf(2) = 1;
    end
    if numel(dimf)<3
      dimf(3) = 1;
    end

    av    = zeros([anatomical.dim            ]);
    allav = zeros([anatomical.dim dimf(2:end)]);

    for k=1:dimf(2)
      for m=1:dimf(3)
        fv     = double_ifnot(dat_array{i}(:,k,m));
        av(:)  = interpmat*fv;
        av(newoutside)   = nan;
        allav(:,:,:,k,m) = av;
      end
    end
    if isfield(interp, 'freq') || isfield(interp, 'time')
      % the output should be a source representation, not a volume
      allav = reshape(allav, prod(anatomical.dim), dimf(2), dimf(3));
    end
    interp = setsubfield(interp, dat_name{i}, allav);
  end

  if ~isfield(interp, 'freq') && ~isfield(interp, 'time')
      % the output should be a volume representation, not a source
    interp = rmfield(interp, 'pos');
  end

elseif ~isUnstructuredFun && isUnstructuredAna
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % functional data defined on a volume, anatomy on a point cloud/mesh
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % set default interpmethod for this situation
  cfg.interpmethod = ft_getopt(cfg, 'interpmethod', 'nearest');
  cfg.sphereradius = ft_getopt(cfg, 'sphereradius', []);
  cfg.power        = ft_getopt(cfg, 'power',        1);

  % interpolate the 3D volume onto the anatomy
  if ~strcmp(cfg.interpmethod, 'project')
    % use interp_gridded
    [interpmat, dummy] = interp_gridded(functional.transform, zeros(functional.dim), anatomical.pos, 'projmethod', cfg.interpmethod, 'sphereradius', cfg.sphereradius, 'inside', functional.inside, 'power', cfg.power);

    % use interp_ungridded
    % interpmat = interp_ungridded(functional.pos, anatomical.pos, 'projmethod', cfg.interpmethod, 'sphereradius', cfg.sphereradius, 'inside', functional.inside, 'power', cfg.power);
  else
    % do the interpolation below, the current implementation of the
    % 'project' method does not output an interpmat (and is therefore quite
    % inefficient

    % set the defaults
    cfg.projvec        = ft_getopt(cfg, 'projvec',       1);
    cfg.projweight     = ft_getopt(cfg, 'projweight',    ones(size(cfg.projvec)));
    cfg.projcomb       = ft_getopt(cfg, 'projcomb',      'mean'); % or max
    cfg.projthresh     = ft_getopt(cfg, 'projthresh',    []);
  end

  % start with an empty structure, keep some fields
  interp = keepfields(functional, {'time', 'freq'});
  interp = copyfields(anatomical, interp, {'pos', 'tri', 'dim', 'transform', 'unit', 'coordsys'});

  % identify the inside voxels after interpolation
  interp.inside    = true(size(anatomical.pos,1),1);

  % interpolate all functional data
  for i=1:length(dat_name)
    fprintf('interpolating %s\n', dat_name{i});

    dimord = getdimord(functional, dat_name{i});
    dimtok = tokenize(dimord, '_');
    dimf   = getdimsiz(functional, dat_name{i});
    dimf(end+1:length(dimtok)) = 1; % there can be additional trailing singleton dimensions

    if prod(functional.dim)==dimf(1)
      % convert into 3-D, 4-D or 5-D array
      dimf        = [functional.dim dimf(2:end)];
      dat_array{i} = reshape(dat_array{i}, dimf);
    end

    % should be 5-D array, can have trailing singleton dimensions
    if numel(dimf)<4
      dimf(4) = 1;
    end
    if numel(dimf)<5
      dimf(5) = 1;
    end

    allav = zeros([size(anatomical.pos,1), dimf(4:end)]);
    if ~strcmp(cfg.interpmethod, 'project')
      for k=1:dimf(4)
        for m=1:dimf(5)
          fv    = double_ifnot(dat_array{i}(:,:,:,k,m)); % ensure double precision to allow sparse multiplication
          fv    = fv(functional.inside(:));
          av    = interpmat*fv;
          allav(:,k,m) = av;
        end
      end
    else
      for k=1:dimf(4)
        for m=1:dimf(5)
          fv   = double_ifnot(dat_array{i}(:,:,:,k,m));
          av   = interp_gridded(functional.transform, fv, anatomical.pos, 'dim', functional.dim, 'projmethod', 'project', 'projvec', cfg.projvec, 'projweight', cfg.projweight, 'projcomb', cfg.projcomb, 'projthresh', cfg.projthresh);
          allav(:,k,m) = av;
        end
      end
    end
    interp = setsubfield(interp, dat_name{i}, allav);
  end

elseif ~isUnstructuredFun && ~isUnstructuredAna
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % functional data defined on a volume, anatomy on a differently sampled volume
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % set default interpmethod for this situation
  cfg.interpmethod = ft_getopt(cfg, 'interpmethod', 'linear');
  if isequal(cfg.interpmethod, 'mode') && isAtlasFun
    % use a mode-based interpolation, i.e. a majority vote of the nearby
    % voxel locations.
    
    % first to a nearest interpolation of the voxel coordinates the other
    % way around
    tmp = anatomical;
    [tx, ty, tz] = voxelcoords(tmp.dim, tmp.transform);
    tmp.tx = tx;
    tmp.ty = ty;
    tmp.tz = tz;
    
    tmpcfg = [];
    tmpcfg.interpmethod = 'nearest';
    tmpcfg.parameter = {'tx' 'ty' 'tz'};
    tmpint = ft_sourceinterpolate(tmpcfg, tmp, functional);
    
    [ix,i1,i2] = intersect([tx(:) ty(:) tz(:)],[tmpint.tx(:) tmpint.ty(:) tmpint.tz(:)],'rows');
    
    interp = keepfields(anatomical, {'pos', 'tri', 'dim', 'transform', 'unit', 'coordsys', 'anatomy'});
    interp.inside = false(interp.dim);
    interp.inside(i1) = true;
    
    for k = 1:numel(cfg.parameter)
      interp.(cfg.parameter{k}) = nan(interp.dim);
      fun = functional.(cfg.parameter{k});
      for m = 1:numel(i1)
        values = reshape(fun(tmpint.tx==tx(i1(m)) & tmpint.ty==ty(i1(m)) & tmpint.tz==tz(i1(m))),[],1);
        [M,f,c] = mode(values);
        if numel(c)==1
          interp.(cfg.parameter{k})(i1(m)) = M;
        else
          ft_warning('multiple modes per voxel, returning NaN');
          interp.(cfg.parameter{k})(i1(m)) = nan;
        end
      end
    end
    
  elseif isequal(cfg.interpmethod, 'mode') && ~isAtlasFun
    ft_error('the interpolation method ''mode'' is only supported for parcellations');

  else
    % start with an empty structure, keep some fields
    interp = keepfields(functional, {'time', 'freq'});
    interp = copyfields(anatomical, interp, {'pos', 'tri', 'dim', 'transform', 'unit', 'coordsys', 'anatomy'});
    
    % convert the anatomical voxel positions into voxel indices into the functional volume
    anatomical.transform = functional.transform \ anatomical.transform;
    functional.transform = eye(4);
    
    [fx, fy, fz] = voxelcoords(functional.dim, functional.transform);
    [ax, ay, az] = voxelcoords(anatomical.dim, anatomical.transform);
    
    % estimate the subvolume of the anatomy that is spanned by the functional volume
    minfx = 1;
    minfy = 1;
    minfz = 1;
    maxfx = functional.dim(1);
    maxfy = functional.dim(2);
    maxfz = functional.dim(3);
    sel = ax(:)>=minfx & ...
      ax(:)<=maxfx & ...
      ay(:)>=minfy & ...
      ay(:)<=maxfy & ...
      az(:)>=minfz & ...
      az(:)<=maxfz;
    fprintf('selecting subvolume of %.1f%%\n', 100*sum(sel)./prod(anatomical.dim));
    
    if all(functional.inside(:))
      % keep all voxels marked as inside
      interp.inside = true(anatomical.dim);
    else
      % reslice and interpolate inside
      interp.inside = zeros(anatomical.dim);
      % interpolate with method nearest
      interp.inside( sel) = my_interpn(double(functional.inside), ax(sel), ay(sel), az(sel), 'nearest', cfg.feedback);
      interp.inside(~sel) = 0;
      interp.inside = logical(interp.inside);
    end
    
    % prepare the grid that is used in the interpolation
    fg = [fx(:) fy(:) fz(:)];
    clear fx fy fz
    
    % reslice and interpolate all functional volumes
    for i=1:length(dat_name)
      fprintf('reslicing and interpolating %s\n', dat_name{i});
      
      dimord = getdimord(functional, dat_name{i});
      dimtok = tokenize(dimord, '_');
      dimf   = getdimsiz(functional, dat_name{i});
      dimf(end+1:length(dimtok)) = 1; % there can be additional trailing singleton dimensions
      
      if prod(functional.dim)==dimf(1)
        % convert into 3-D, 4-D or 5-D array
        dimf = [functional.dim dimf(2:end)];
        dat_array{i} = reshape(dat_array{i}, dimf);
      end
      
      % should be 5-D array, can have trailing singleton dimensions
      if numel(dimf)<4
        dimf(4) = 1;
      end
      if numel(dimf)<5
        dimf(5) = 1;
      end
      
      av    = zeros([anatomical.dim            ]);
      allav = zeros([anatomical.dim dimf(4:end)]);
      functional.inside = functional.inside(:,:,:,1,1);
      
      if any(dimf(4:end)>1) && ~strcmp(cfg.feedback, 'none')
        % this is needed to prevent feedback to be displayed for every time-frequency point
        ft_warning('disabling feedback');
        cfg.feedback = 'none';
      end
      
      for k=1:dimf(4)
        for m=1:dimf(5)
          fv = double_ifnot(dat_array{i}(:,:,:,k,m));
          % av( sel) = my_interpn(fx, fy, fz, fv, ax(sel), ay(sel), az(sel), cfg.interpmethod, cfg.feedback);
          if islogical(dat_array{i})
            % interpolate always with method nearest
            av( sel) = my_interpn(fv, ax(sel), ay(sel), az(sel), 'nearest', cfg.feedback);
            av = logical(av);
          else
            if ~all(functional.inside(:))
              % extrapolate the outside of the functional volumes for better interpolation at the edges
              fv(~functional.inside) = griddatan(fg(functional.inside(:), :), fv(functional.inside(:)), fg(~functional.inside(:), :), 'nearest');
            end
            % interpolate functional onto anatomical grid
            av( sel) = my_interpn(fv, ax(sel), ay(sel), az(sel), cfg.interpmethod, cfg.feedback);
            av(~sel) = nan;
            av(~interp.inside) = nan;
          end
          allav(:,:,:,k,m) = av;
        end
      end
      if (isfield(interp, 'freq') && numel(interp.freq)>1) || (isfield(interp, 'time') && numel(interp.time)>1)
        % the output should be a source representation, not a volume
        allav = reshape(allav, prod(anatomical.dim), dimf(4), dimf(5));
      end
      interp = setsubfield(interp, dat_name{i}, allav);
      % keep the description of the labels in the segmentation/parcellation
      if strcmp(cfg.interpmethod, 'nearest') && isfield(functional, [dat_name{i} 'label'])
        interp.([dat_name{i} 'label']) = functional.([dat_name{i} 'label']);
      end
    end
    
  end

end % computing the interpolation according to the input data

if isfield(interp, 'freq') || isfield(interp, 'time')
  % the output should be a source representation, not a volumetric representation
  if ~isfield(interp, 'pos')
    [x, y, z] = voxelcoords(interp.dim, interp.transform);
    interp.pos = [x(:) y(:) z(:)];
  end
end

if isAtlasFun
  for i=1:numel(dat_name)
    % keep the labels that describe the different tissue types
    interp = copyfields(functional, interp, [dat_name{i} 'label']);
    % replace NaNs that fall outside the labeled area with zero
    tmp = interp.(dat_name{i});
    tmp(isnan(tmp)) = 0;
    interp.(dat_name{i}) = tmp;
  end
  % remove the inside field if present
  interp = removefields(interp, 'inside');
end

if exist('interpmat', 'var')
  cfg.interpmat = interpmat;
  cfg.interpmat; % access it once to fool the cfg-tracking
end

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble previous   functional anatomical
ft_postamble provenance interp
ft_postamble history    interp
ft_postamble savevar    interp


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION this function computes the location of all voxels in head
% coordinates in a memory efficient manner
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [x, y, z] = voxelcoords(dim, transform)

xgrid     = 1:dim(1);
ygrid     = 1:dim(2);
zgrid     = 1:dim(3);
npix      = prod(dim(1:2));  % number of voxels in a single slice

x = zeros(dim);
y = zeros(dim);
z = zeros(dim);
X = zeros(1,npix);
Y = zeros(1,npix);
Z = zeros(1,npix);
E = ones(1,npix);
% determine the voxel locations per slice
for i=1:dim(3)
  [X(:), Y(:), Z(:)] = ndgrid(xgrid, ygrid, zgrid(i));
  tmp = transform*[X; Y; Z; E];
  x((1:npix)+(i-1)*npix) = tmp(1,:);
  y((1:npix)+(i-1)*npix) = tmp(2,:);
  z((1:npix)+(i-1)*npix) = tmp(3,:);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION for memory efficient interpolation
% the only reason for this function is that it does the interpolation in smaller chuncks
% this prevents memory problems that I often encountered here
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function [av] = my_interpn(fx, fy, fz, fv, ax, ay, az, interpmethod, feedback);
function [av] = my_interpn(fv, ax, ay, az, interpmethod, feedback)

num = numel(ax);            % total number of voxels
blocksize = floor(num/20);  % number of voxels to interpolate at once, split it into 20 chuncks
lastblock = 0;              % boolean flag for while loop
sel = 1:blocksize;          % selection of voxels that are interpolated, this is the first chunck
av  = zeros(size(ax));
ft_progress('init', feedback, 'interpolating');
while (1)
  ft_progress(sel(1)/num, 'interpolating %.1f%%\n', 100*sel(1)/num);
  if sel(end)>=num
    sel = sel(1):num;
    lastblock = 1;
  end
  av(sel) = interpn(fv, ax(sel), ay(sel), az(sel), interpmethod);
  if lastblock
    break
  end
  sel = sel + blocksize;
end
ft_progress('close');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION to cast array to double precision, only if needed
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function out = double_ifnot(in)

if ~isa(in, 'double')
  out = double(in);
else
  out = in;
end
