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