function [cfg] = ft_sourceplot(cfg, functional, anatomical)

% FT_SOURCEPLOT plots functional source reconstruction data on slices or on a surface,

% optionally as an overlay on anatomical MRI data, where statistical data can be used to

% determine the opacity of the mask. Input data comes from FT_SOURCEANALYSIS,

% FT_SOURCEGRANDAVERAGE or statistical values from FT_SOURCESTATISTICS.

%

% Use as

% ft_sourceplot(cfg, anatomical)

% ft_sourceplot(cfg, functional)

% ft_sourceplot(cfg, functional, anatomical)

% where the input data can contain either anatomical, functional or statistical data,

% or a combination of them.

%

% The input data can be in a 3-D volumetric representation or in a 2-D cortical sheet

% representation. If both anatomical and functional/statistical data is provided as input,

% they should be represented in the same coordinate system or interpolated on the same

% geometrical representation, e.g. using FT_SOURCEINTERPOLATE.

%

% The slice and ortho visualization plot the data in the input data voxel arrangement, i.e.

% the three ortho views are the 1st, 2nd and 3rd dimension of the 3-D data matrix, not of

% the head coordinate system. The specification of the coordinate for slice intersection

% is specified in head coordinates, i.e. relative to anatomical landmarks or fiducials and

% expressed in mm or cm. If you want the visualisation to be consistent with the head

% coordinate system, you can reslice the data using FT_VOLUMERESLICE. See http://bit.ly/1OkDlVF

%

% The configuration should contain:

% cfg.method = 'ortho', plots the data on three orthogonal slices

% 'slice', plots the data on a number of slices in the same plane

% 'surface', plots the data on a 3D brain surface

% 'glassbrain', plots a max-projection through the brain

% 'vertex', plots the grid points or vertices scaled according to the functional value

% 'cloud', plot the data as clouds, spheres, or points scaled according to the functional value

% and

% cfg.anaparameter = string, field in data with the anatomical data (default = 'anatomy' if present in data)

% cfg.funparameter = string, field in data with the functional parameter of interest (default = [])

% cfg.maskparameter = string, field in the data to be used for opacity masking of fun data (default = [])

% If values are between 0 and 1, zero is fully transparant and one is fully opaque.

% If values in the field are not between 0 and 1 they will be scaled depending on the values

% of cfg.opacitymap and cfg.opacitylim (see below)

% You can use masking in several ways, f.i.

% - use outcome of statistics to show only the significant values and mask the insignificant

% NB see also cfg.opacitymap and cfg.opacitylim below

% - use the functional data itself as mask, the highest value (and/or lowest when negative)

% will be opaque and the value closest to zero transparent

% - Make your own field in the data with values between 0 and 1 to control opacity directly

%

% The following parameters can be used in all methods:

% cfg.downsample = downsampling for resolution reduction, integer value (default = 1) (orig: from surface)

% cfg.atlas = string, filename of atlas to use (default = []) see FT_READ_ATLAS

% for ROI masking (see 'masking' below) or for orthogonal plots (see method='ortho' below)

% cfg.visible = string, 'on' or 'off' whether figure will be visible (default = 'on')

% cfg.figure = 'yes' or 'no', whether to open a new figure. You can also specify a figure handle from FIGURE, GCF or SUBPLOT. (default = 'yes')

% cfg.position = location and size of the figure, specified as [left bottom width height] (default is automatic)

% cfg.renderer = string, 'opengl', 'zbuffer', 'painters', see RENDERERINFO. The OpenGL renderer is required when using opacity (default = 'opengl')

%

% The following parameters can be used for the functional data:

% cfg.funcolormap = colormap for functional data, see COLORMAP (default = 'auto')

% 'auto', depends structure funparameter, or on funcolorlim

% - funparameter: only positive values, or funcolorlim:'zeromax' -> 'hot'

% - funparameter: only negative values, or funcolorlim:'minzero' -> 'cool'

% - funparameter: both pos and neg values, or funcolorlim:'maxabs' -> 'default'

% - funcolorlim: [min max] if min & max pos-> 'hot', neg-> 'cool', both-> 'default'

% cfg.funcolorlim = color range of the functional data (default = 'auto')

% [min max]

% 'maxabs', from -max(abs(funparameter)) to +max(abs(funparameter))

% 'zeromax', from 0 to max(funparameter)

% 'minzero', from min(funparameter) to 0

% 'auto', if funparameter values are all positive: 'zeromax',

% all negative: 'minzero', both possitive and negative: 'maxabs'

% cfg.colorbar = 'yes' or 'no' (default = 'yes')

% cfg.colorbartext = string indicating the text next to colorbar

%

% The 'ortho' method can also plot time and/or frequency, the other methods can not.

% If your functional data has a time and/or frequency dimension, you can use

% cfg.latency = scalar or string, can be 'all', 'prestim', 'poststim', or [beg end], specify time range in seconds

% cfg.avgovertime = string, can be 'yes' or 'no' (default = 'no')

% cfg.frequency = scalar or string, can be 'all', or [beg end], specify frequency range in Hz

% cfg.avgoverfreq = string, can be 'yes' or 'no' (default = 'no')

%

% The following parameters can be used for the masking data:

% cfg.maskstyle = 'opacity', or 'colormix'. If 'opacity', low-level

% graphics opacity masking is applied, if

% 'colormix', the color data is explicitly

% expressed as a single RGB value, incorporating

% the opacitymask. Yields faster and more robust

% rendering in general.

% cfg.opacitymap = opacitymap for mask data, see ALPHAMAP (default = 'auto')

% 'auto', depends structure maskparameter, or on opacitylim

% - maskparameter: only positive values, or opacitylim:'zeromax' -> 'rampup'

% - maskparameter: only negative values, or opacitylim:'minzero' -> 'rampdown'

% - maskparameter: both pos and neg values, or opacitylim:'maxabs' -> 'vdown'

% - opacitylim: [min max] if min & max pos-> 'rampup', neg-> 'rampdown', both-> 'vdown'

% - NB. to use p-values use 'rampdown' to get lowest p-values opaque and highest transparent

% cfg.opacitylim = range of mask values to which opacitymap is scaled (default = 'auto')

% [min max]

% 'maxabs', from -max(abs(maskparameter)) to +max(abs(maskparameter))

% 'zeromax', from 0 to max(abs(maskparameter))

% 'minzero', from min(abs(maskparameter)) to 0

% 'auto', if maskparameter values are all positive: 'zeromax',

% all negative: 'minzero', both positive and negative: 'maxabs'

% cfg.roi = string or cell of strings, region(s) of interest from anatomical atlas (see cfg.atlas above)

% everything is masked except for ROI

%

% When cfg.method='ortho', three orthogonal slices will be rendered. You can click in any

% of the slices to update the display in the other two. You can also use the arrow keys on

% your keyboard to navigate in one-voxel steps. Note that the slices are along the first,

% second and third voxel dimension, which do not neccessarily correspond to the axes of the

% head coordinate system. See http://bit.ly/1OkDlVF

%

% The following parameters apply when cfg.method='ortho'

% cfg.location = location of cut, (default = 'auto')

% 'auto', 'center' if only anatomy, 'max' if functional data

% 'min' and 'max' position of min/max funparameter

% 'center' of the brain

% [x y z], coordinates in voxels or head, see cfg.locationcoordinates

% cfg.locationcoordinates = coordinate system used in cfg.location, 'head' or 'voxel' (default = 'head')

% 'head', headcoordinates as mm or cm

% 'voxel', voxelcoordinates as indices

% cfg.crosshair = 'yes' or 'no' (default = 'yes')

% cfg.axis = 'on' or 'off' (default = 'on')

% cfg.queryrange = number, in atlas voxels (default = 1)

% cfg.clim = lower and upper anatomical MRI limits (default = [0 1])

%

% When cfg.method='slice', a NxM montage with a large number of slices will be rendered.

% All slices are evenly spaced and along the same dimension.

%

% The following parameters apply for cfg.method='slice'

% cfg.nslices = number of slices, (default = 20)

% cfg.slicerange = range of slices in data, (default = 'auto')

% 'auto', full range of data

% [min max], coordinates of first and last slice in voxels

% cfg.slicedim = dimension to slice 1 (x-axis) 2(y-axis) 3(z-axis) (default = 3)

% cfg.title = string, title of the plot

% cfg.figurename = string, title of the figure window

%

% When cfg.method='surface', the functional data will be rendered onto a cortical mesh

% (can be an inflated mesh). If the input source data contains a tri-field (i.e. a

% description of a mesh), no interpolation is needed. If the input source data does not

% contain a tri-field, an interpolation is performed onto a specified surface. Note that

% the coordinate system in which the surface is defined should be the same as the coordinate

% system that is represented in the pos-field.

%

% The following parameters apply to cfg.method='surface' when an interpolation is required

% cfg.surffile = string, file that contains the surface (default = 'surface_white_both.mat')

% 'surface_white_both.mat' contains a triangulation that corresponds with the

% SPM anatomical template in MNI coordinates

% cfg.surfinflated = string, file that contains the inflated surface (default = [])

% may require specifying a point-matching (uninflated) surffile

% cfg.surfdownsample = number (default = 1, i.e. no downsampling)

% cfg.projmethod = projection method, how functional volume data is projected onto surface

% 'nearest', 'project', 'sphere_avg', 'sphere_weighteddistance'

% cfg.projvec = vector (in mm) to allow different projections that

% are combined with the method specified in cfg.projcomb

% cfg.projcomb = 'mean', 'max', method to combine the different projections

% cfg.projweight = vector of weights for the different projections (default = 1)

% cfg.projthresh = implements thresholding on the surface level

% for example, 0.7 means 70% of maximum

% cfg.sphereradius = maximum distance from each voxel to the surface to be

% included in the sphere projection methods, expressed in mm

% cfg.distmat = precomputed distance matrix (default = [])

%

% The following parameters apply to cfg.method='surface' irrespective of whether an interpolation is required

% cfg.camlight = 'yes' or 'no' (default = 'yes')

% cfg.facecolor = [r g b] values or string, for example 'brain', 'cortex', 'skin', 'black', 'red', 'r',

% or an Nx3 or Nx1 array where N is the number of faces

% cfg.vertexcolor = [r g b] values or string, for example 'brain', 'cortex', 'skin', 'black', 'red', 'r',

% or an Nx3 or Nx1 array where N is the number of vertices

% cfg.edgecolor = [r g b] values or string, for example 'brain', 'cortex', 'skin', 'black', 'red', 'r'

%

% When cfg.method = 'cloud', the functional data will be rendered as as clouds (groups of points),

% spheres, or single points at each sensor position. These spheres or point clouds can either be

% viewed in 3D or as 2D slices. The 'anatomical' input may also consist of a single or multiple

% triangulated surface meshes in an Nx1 cell-array to be plotted with the interpolated functional

% data (see FT_PLOT_CLOUD).

%

% The following parameters apply to cfg.method='cloud'

% cfg.cloudtype = 'point' plots a single point at each sensor position

% 'cloud' (default) plots each a group of spherically arranged points at each sensor position

% 'surf' plots a single spherical surface mesh at each sensor position

% cfg.radius = scalar, maximum radius of cloud (default = 4)

% cfg.colorgrad = 'white' or a scalar (e.g. 1), degree to which color of points in cloud

% changes from its center

% cfg.slice = requires 'anatomical' as input (default = 'none')

% '2d', plots 2D slices through the cloud with an outline of the mesh

% '3d', draws an outline around the mesh at a particular slice

% cfg.ori = 'x', 'y', or 'z', specifies the orthogonal plane which will be plotted (default = 'y')

% cfg.slicepos = 'auto' or Nx1 vector specifying the position of the

% slice plane along the orientation axis (default = 'auto': chooses slice(s) with

% the most data)

% cfg.nslices = scalar, number of slices to plot if 'slicepos' = 'auto (default = 1)

% cfg.minspace = scalar, minimum spacing between slices if nslices>1 (default = 1)

%

% To facilitate data-handling and distributed computing you can use

% cfg.inputfile = ...

% If you specify this option the input data will be read from a *.mat file on

% disk. This mat files should contain only a single variable corresponding to the

% input structure.

%

% See also FT_SOURCEMOVIE, FT_SOURCEANALYSIS, FT_SOURCEGRANDAVERAGE, FT_SOURCESTATISTICS,

% FT_VOLUMELOOKUP, FT_READ_ATLAS, FT_READ_MRI

% TODO have to be built in:

% cfg.marker = [Nx3] array defining N marker positions to display (orig: from sliceinterp)

% cfg.markersize = radius of markers (default = 5)

% cfg.markercolor = [1x3] marker color in RGB (default = [1 1 1], i.e. white) (orig: from sliceinterp)

% white background option

% undocumented TODO

% slice in all directions

% surface also optimal when inside present

% come up with a good glass brain projection

%

% undocumented option

% cfg.intersectmesh = cell-array of mesh(es) to be plotted along with the

% anatomy, useful for evaluating coregistration. Does

% at present not check for coordinate system

% Copyright (C) 2007-2022, Robert Oostenveld, Ingrid Nieuwenhuis, J.M.

% 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(functional)

ft_error('please use cfg.inputfile instead of specifying the input variable as a sting');

end

% ensure that old and unsupported options are not being relied on by the end-user's script

cfg = ft_checkconfig(cfg, 'renamedval', {'funparameter', 'avg.pow', 'pow'});

cfg = ft_checkconfig(cfg, 'renamedval', {'funparameter', 'avg.coh', 'coh'});

cfg = ft_checkconfig(cfg, 'renamedval', {'funparameter', 'avg.mom', 'mom'});

cfg = ft_checkconfig(cfg, 'renamedval', {'maskparameter', 'avg.pow', 'pow'});

cfg = ft_checkconfig(cfg, 'renamedval', {'maskparameter', 'avg.coh', 'coh'});

cfg = ft_checkconfig(cfg, 'renamedval', {'maskparameter', 'avg.mom', 'mom'});

cfg = ft_checkconfig(cfg, 'renamedval', {'location', 'interactive', 'auto'});

% instead of specifying cfg.coordsys, the user should specify the coordsys in the data

cfg = ft_checkconfig(cfg, 'forbidden', {'units', 'coordsys', 'inputcoord', 'inputcoordsys', 'coordinates'});

% see http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=2837

cfg = ft_checkconfig(cfg, 'renamed', {'viewdim', 'axisratio'});

cfg = ft_checkconfig(cfg, 'renamed', {'newfigure', 'figure'});

if isfield(cfg, 'atlas') && ~isempty(cfg.atlas)

% the atlas lookup requires the specification of the coordsys

functional = ft_checkdata(functional, 'datatype', {'source', 'volume'}, 'feedback', 'yes', 'hasunit', 'yes', 'hascoordsys', 'yes');

else

% check if the input functional is valid for this function, a coordsys is not directly needed

functional = ft_checkdata(functional, 'datatype', {'source', 'volume'}, 'feedback', 'yes', 'hasunit', 'yes');

end

% set the defaults for all methods

cfg.method = ft_getopt(cfg, 'method', 'ortho');

cfg.funparameter = ft_getopt(cfg, 'funparameter', []);

cfg.maskparameter = ft_getopt(cfg, 'maskparameter', []);

cfg.maskstyle = ft_getopt(cfg, 'maskstyle', 'opacity');

cfg.downsample = ft_getopt(cfg, 'downsample', 1);

cfg.flip = ft_getopt(cfg, 'flip', []); % the default is set below

cfg.title = ft_getopt(cfg, 'title', []);

cfg.figurename = ft_getopt(cfg, 'figurename', []);

cfg.atlas = ft_getopt(cfg, 'atlas', []);

cfg.marker = ft_getopt(cfg, 'marker', []);

cfg.markersize = ft_getopt(cfg, 'markersize', 5);

cfg.markercolor = ft_getopt(cfg, 'markercolor', [1 1 1]);

cfg.colorbar = ft_getopt(cfg, 'colorbar', 'yes');

cfg.colorbartext = ft_getopt(cfg, 'colorbartext', '');

cfg.voxelratio = ft_getopt(cfg, 'voxelratio', 'data'); % display size of the voxel, 'data' or 'square'

cfg.axisratio = ft_getopt(cfg, 'axisratio', 'data'); % size of the axes of the three orthoplots, 'square', 'voxel', or 'data'

cfg.visible = ft_getopt(cfg, 'visible', 'on');

cfg.clim = ft_getopt(cfg, 'clim', [0 1]); % this is used to scale the orthoplot

cfg.intersectmesh = ft_getopt(cfg, 'intersectmesh');

cfg.renderer = ft_getopt(cfg, 'renderer', 'opengl');

cfg.interactive = ft_getopt(cfg, 'interactive', 'yes'); % used to disable interaction for method=glassbrain

if isempty(cfg.flip)

if strcmp(cfg.method, 'ortho')

cfg.flip = 'yes';

else

cfg.flip = 'no';

end

end

if ~isfield(cfg, 'anaparameter')

if isfield(functional, 'anatomy')

cfg.anaparameter = 'anatomy';

else

cfg.anaparameter = [];

end

end

% set the common defaults for the functional data

cfg.funcolormap = ft_getopt(cfg, 'funcolormap', 'auto');

cfg.funcolorlim = ft_getopt(cfg, 'funcolorlim', 'auto');

% set the common defaults for the statistical data

cfg.opacitymap = ft_getopt(cfg, 'opacitymap', 'auto');

cfg.opacitylim = ft_getopt(cfg, 'opacitylim', 'auto');

cfg.roi = ft_getopt(cfg, 'roi', []);

cfg.maskstyle = ft_getopt(cfg, 'maskstyle', 'opacity');

% select the functional and the mask parameter

cfg.funparameter = parameterselection(cfg.funparameter, functional);

cfg.maskparameter = parameterselection(cfg.maskparameter, functional);

% only a single parameter should be selected

try, cfg.funparameter = cfg.funparameter{1}; end

try, cfg.maskparameter = cfg.maskparameter{1}; end

if isfield(functional, 'time') || isfield(functional, 'freq')

% make a selection of the time and/or frequency dimension

tmpcfg = keepfields(cfg, {'frequency', 'avgoverfreq', 'keepfreqdim', 'latency', 'avgovertime', 'keeptimedim', 'showcallinfo', 'trackcallinfo', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo', 'checksize'});

functional = ft_selectdata(tmpcfg, functional);

% restore the provenance information

[cfg, functional] = rollback_provenance(cfg, functional);

end

% the data can be passed as input argument or can be read from disk

hasanatomical = exist('anatomical', 'var');

if hasanatomical && ~strcmp(cfg.method, 'cloud') % cloud method should be able to take multiple surfaces and does not require interpolation

% interpolate on the fly, this also does the downsampling if requested

tmpcfg = keepfields(cfg, {'downsample', 'interpmethod', 'sphereradius', 'showcallinfo', 'trackcallinfo', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo', 'checksize'});

tmpcfg.parameter = cfg.funparameter;

functional = ft_sourceinterpolate(tmpcfg, functional, anatomical);

[cfg, functional] = rollback_provenance(cfg, functional);

cfg.anaparameter = 'anatomy';

elseif ~hasanatomical && cfg.downsample~=1

% optionally downsample the functional volume

tmpcfg = keepfields(cfg, {'downsample', 'showcallinfo', 'trackcallinfo', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo', 'checksize'});

tmpcfg.parameter = {cfg.funparameter, cfg.maskparameter, cfg.anaparameter};

functional = ft_volumedownsample(tmpcfg, functional);

[cfg, functional] = rollback_provenance(cfg, functional);

end

if isfield(functional, 'dim') && isfield(functional, 'transform')

% this is a regular 3D functional volume

isUnstructuredFun = false;

elseif isfield(functional, 'dim') && isfield(functional, 'pos')

% these are positions that can be mapped onto a 3D regular grid

isUnstructuredFun = false;

% construct the transformation matrix from the positions

functional.transform = pos2transform(functional.pos, functional.dim);

else

% this is functional data on irregular positions, such as a cortical sheet

isUnstructuredFun = true;

end

if ~isUnstructuredFun && strcmp(cfg.flip, 'yes')

% align the anatomical volume approximately to coordinate system, this puts it upright

origmethod = cfg.method;

tmpcfg = [];

tmpcfg.method = 'flip';

tmpcfg.trackcallinfo = 'no';

tmpcfg.showcallinfo = 'no';

functional = ft_volumereslice(tmpcfg, functional);

[cfg, functional] = rollback_provenance(cfg, functional);

cfg.method = origmethod;

end

% this only relates to the dimensions of the geometry, which is npos*1 or nx*ny*nz

if isUnstructuredFun

dim = [size(functional.pos,1) 1];

else

dim = functional.dim;

end

%% get the elements that will be plotted

hasatlas = ~isempty(cfg.atlas);

if hasatlas

[atlas, functional] = handle_atlas_input(cfg.atlas, functional);

end

hasroi = ~isempty(cfg.roi);

if hasroi

if ~hasatlas

ft_error('specify cfg.atlas which specifies cfg.roi')

else

% get the mask

tmpcfg = keepfields(cfg, {'roi', 'atlas'});

roi = ft_volumelookup(tmpcfg, functional);

end

end

hasana = isfield(functional, cfg.anaparameter);

if hasana

ana = getsubfield(functional, cfg.anaparameter);

if isa(ana, 'uint8') || isa(ana, 'uint16') || isa(ana, 'int8') || isa(ana, 'int16')

ana = double(ana);

end

fprintf('scaling anatomy to [0 1]\n');

dmin = min(ana(:));

dmax = max(ana(:));

ana = (ana-dmin)./(dmax-dmin);

ana = reshape(ana, dim);

end

%%% funparameter

hasfun = isfield(functional, cfg.funparameter);

if hasfun

fun = getsubfield(functional, cfg.funparameter);

dimord = getdimord(functional, cfg.funparameter);

dimtok = tokenize(dimord, '_');

% replace the cell-array functional with a normal array

if startsWith(dimord, '{pos}')

tmpdim = getdimsiz(functional, cfg.funparameter);

tmpfun = nan(tmpdim);

insideindx = find(functional.inside);

for i=insideindx(:)'

tmpfun(i,:) = fun{i};

end

fun = tmpfun;

clear tmpfun

dimord = ['pos' dimord(6:end)]; % update the description of the dimensions

dimtok = tokenize(dimord, '_');

elseif startsWith(dimord, '{dim1_dim2_dim3}')

tmpdim = getdimsiz(functional, cfg.funparameter);

tmpfun = nan(tmpdim);

for i1=1:functional.dim(1)

for i2=1:functional.dim(2)

for i3=1:functional.dim(3)

if functional.inside(i1, i2, i3)

tmpfun(i1,i2,i3,:) = fun{i1, i2, i3};

end

end

end

end

fun = tmpfun;

clear tmpfun

dimord = ['dim1_dim2_dim3' dimord(17:end)]; % update the description of the dimensions

dimtok = tokenize(dimord, '_');

end

% ensure that the functional data is real

if ~isreal(fun)

ft_warning('functional data is complex, taking absolute value');

fun = abs(fun);

end

if ~isa(fun, 'double')

ft_warning('converting functional data to double precision');

fun = double(fun);

end

if strcmp(dimord, 'pos_rgb') || (ndims(fun)>3 && size(fun,4)==3)

% treat functional data as rgb values

if any(fun(:)>1 | fun(:)<0)

% scale

tmpdim = size(fun);

nvox = prod(tmpdim(1:end-1));

tmpfun = reshape(fun,[nvox tmpdim(end)]);

m1 = max(tmpfun,[],1);

m2 = min(tmpfun,[],1);

tmpfun = (tmpfun-m2(ones(nvox,1),:))./(m1(ones(nvox,1),:)-m2(ones(nvox,1),:));

fun = reshape(tmpfun, tmpdim);

clear tmpfun

end

qi = 1;

hasfreq = 0;

hastime = 0;

doimage = 1;

fcolmin = 0;

fcolmax = 1;

cfg.funcolormap = 'rgb';

else

% determine scaling min and max (fcolmin fcolmax) and funcolormap

if ~isa(fun, 'logical')

funmin = min(fun(:));

funmax = max(fun(:));

else

funmin = 0;

funmax = 1;

end

% smart automatic limits

if isequal(cfg.funcolorlim, 'auto')

if sign(funmin)>-1 && sign(funmax)>-1

cfg.funcolorlim = 'zeromax';

elseif sign(funmin)<1 && sign(funmax)<1

cfg.funcolorlim = 'minzero';

else

cfg.funcolorlim = 'maxabs';

end

end

if ischar(cfg.funcolorlim)

% limits are given as string

if isequal(cfg.funcolorlim, 'maxabs')

fcolmin = -max(abs([funmin,funmax]));

fcolmax = max(abs([funmin,funmax]));

if isequal(cfg.funcolormap, 'auto'); cfg.funcolormap = 'default'; end

elseif isequal(cfg.funcolorlim, 'zeromax')

fcolmin = 0;

fcolmax = funmax;

if isequal(cfg.funcolormap, 'auto'); cfg.funcolormap = 'hot'; end

elseif isequal(cfg.funcolorlim, 'minzero')

fcolmin = funmin;

fcolmax = 0;

if isequal(cfg.funcolormap, 'auto'); cfg.funcolormap = 'cool'; end

else

ft_error('do not understand cfg.funcolorlim');

end

else

% limits are numeric

fcolmin = cfg.funcolorlim(1);

fcolmax = cfg.funcolorlim(2);

% smart colormap

if isequal(cfg.funcolormap, 'auto')

if sign(fcolmin) == -1 && sign(fcolmax) == 1

cfg.funcolormap = 'default';

else

if fcolmin < 0

cfg.funcolormap = 'cool';

else

cfg.funcolormap = 'hot';

end

end

end

end % if ischar

clear funmin funmax

% what if fun is 4D?

if ndims(fun)>3 || prod(dim)==size(fun,1)

if strcmp(dimord, 'pos_freq_time') || strcmp(dimord, 'dim1_dim2_dim3_freq_time')

% functional contains time-frequency representation

qi = [1 1];

hasfreq = numel(functional.freq)>1;

hastime = numel(functional.time)>1;

fun = reshape(fun, [dim numel(functional.freq) numel(functional.time)]);

elseif strcmp(dimord, 'pos_time') || strcmp(dimord, 'dim1_dim2_dim3_time')

% functional contains evoked field

qi = 1;

hasfreq = 0;

hastime = numel(functional.time)>1;

fun = reshape(fun, [dim numel(functional.time)]);

elseif strcmp(dimord, 'pos_ori_time') || strcmp(dimord, 'dim1_dim2_dim3_ori_time')

% functional contains evoked field

qi = 1;

hasfreq = 0;

hastime = numel(functional.time)>1;

% the following will fail if the number of orientations is larger than 1

fun = reshape(fun, [dim numel(functional.time)]);

elseif strcmp(dimord, 'pos_freq') || strcmp(dimord, 'dim1_dim2_dim3_freq')

% functional contains frequency spectra

qi = 1;

hasfreq = numel(functional.freq)>1;

hastime = 0;

fun = reshape(fun, [dim numel(functional.freq)]);

else

% functional contains scalar value for each position

qi = 1;

hasfreq = 0;

hastime = 0;

fun = reshape(fun, dim);

end

else

% do nothing

qi = 1;

hasfreq = 0;

hastime = 0;

end

doimage = 0;

end % if dimord has rgb or something else

else

% there is no functional data

qi = 1;

hasfreq = 0;

hastime = 0;

doimage = 0;

fcolmin = 0; % needs to be defined for callback

fcolmax = 1;

end

hasmsk = issubfield(functional, cfg.maskparameter);

if hasmsk

if ~hasfun

ft_error('you can not have a mask without functional data')

else

msk = getsubfield(functional, cfg.maskparameter);

if islogical(msk) % otherwise sign() not posible

msk = double(msk);

end

end

% reshape to match fun

if strcmp(dimord, 'pos_freq_time')

% functional contains timefrequency representation

msk = reshape(msk, [dim numel(functional.freq) numel(functional.time)]);

elseif strcmp(dimord, 'pos_time')

% functional contains evoked field

msk = reshape(msk, [dim numel(functional.time)]);

elseif strcmp(dimord, 'pos_freq')

% functional contains frequency spectra

msk = reshape(msk, [dim numel(functional.freq)]);

else

msk = reshape(msk, dim);

end

% determine scaling and opacitymap

mskmin = min(msk(:));

mskmax = max(msk(:));

% determine the opacity limits and the opacity map

% smart limits: make from auto other string, or equal to funcolorlim if funparameter == maskparameter

if isequal(cfg.opacitylim, 'auto')

if isequal(cfg.funparameter,cfg.maskparameter)

cfg.opacitylim = cfg.funcolorlim;

else

if sign(mskmin)>-1 && sign(mskmax)>-1

cfg.opacitylim = 'zeromax';

elseif sign(mskmin)<1 && sign(mskmax)<1

cfg.opacitylim = 'minzero';

else

cfg.opacitylim = 'maxabs';

end

end

end

if ischar(cfg.opacitylim)

% limits are given as string

switch cfg.opacitylim

case 'zeromax'

opacmin = 0;

opacmax = mskmax;

if isequal(cfg.opacitymap, 'auto'), cfg.opacitymap = 'rampup'; end

case 'minzero'

opacmin = mskmin;

opacmax = 0;

if isequal(cfg.opacitymap, 'auto'), cfg.opacitymap = 'rampdown'; end

case 'maxabs'

opacmin = -max(abs([mskmin, mskmax]));

opacmax = max(abs([mskmin, mskmax]));

if isequal(cfg.opacitymap, 'auto'), cfg.opacitymap = 'vdown'; end

otherwise

ft_error('incorrect specification of cfg.opacitylim');

end % switch opacitylim

else

% limits are numeric

opacmin = cfg.opacitylim(1);

opacmax = cfg.opacitylim(2);

if isequal(cfg.opacitymap, 'auto')

if sign(opacmin)>-1 && sign(opacmax)>-1

cfg.opacitymap = 'rampup';

elseif sign(opacmin)<1 && sign(opacmax)<1

cfg.opacitymap = 'rampdown';

else

cfg.opacitymap = 'vdown';

end

end

end % handling opacitylim and opacitymap

clear mskmin mskmax

else

opacmin = [];

opacmax = [];

end

% prevent outside fun from being plotted

if hasfun && ~hasmsk && isfield(functional, 'inside')

hasmsk = 1;

msk = zeros(dim);

cfg.opacitymap = 'rampup';

opacmin = 0;

opacmax = 1;

% make intelligent mask

if isequal(cfg.method, 'surface')

msk(functional.inside&isfinite(functional.(cfg.funparameter))) = 1;

if any(functional.inside&~isfinite(functional.(cfg.funparameter)))

ft_warning('functional data contains %d NaNs labeled as inside', sum(functional.inside&~isfinite(functional.(cfg.funparameter))));

end

else

if hasana

msk(functional.inside) = 0.5; % so anatomy is visible

else

msk(functional.inside) = 1;

end

end

end

% if region of interest is specified, mask everything besides roi

if hasfun && hasroi && ~hasmsk

hasmsk = 1;

msk = roi;

cfg.opacitymap = 'rampup';

opacmin = 0;

opacmax = 1;

elseif hasfun && hasroi && hasmsk

msk = roi .* msk;

opacmin = [];

opacmax = []; % has to be defined

elseif hasroi

ft_error('you can not have a roi without functional data')

end

%% give some feedback

if ~hasfun && ~hasana

% this seems to be a problem that people often have due to incorrect specification of the cfg

ft_error('no anatomy is present and no functional data is selected, please check your cfg.funparameter');

end

if ~hasana

fprintf('not plotting anatomy\n');

end

if ~hasfun

fprintf('not plotting functional data\n');

end

if ~hasmsk

fprintf('not applying a mask on the functional data\n');

end

if ~hasatlas

fprintf('not using an atlas\n');

end

if ~hasroi

fprintf('not using a region-of-interest\n');

end

%% start building the figure

% open a new figure with the specified settings

h = open_figure(keepfields(cfg, {'figure', 'position', 'visible', 'renderer', 'figurename', 'title'}));

set(h, 'color', [1 1 1]);

%% set color and opacity mapping for this figure

if hasfun

if ischar(cfg.funcolormap) && ~strcmp(cfg.funcolormap, 'rgb')

cfg.funcolormap = ft_colormap(cfg.funcolormap);

elseif iscell(cfg.funcolormap)

cfg.funcolormap = ft_colormap(cfg.funcolormap{:});

end

end

if hasmsk

cfg.opacitymap = alphamap(cfg.opacitymap);

alphamap(cfg.opacitymap);

if ndims(fun)>3 && ndims(msk)==3 && ~isequal(cfg.funcolormap, 'rgb')

siz = size(fun);

msk = repmat(msk, [1 1 1 siz(4:end)]);

end

end

switch cfg.method

case 'slice'

assert(~hastime, 'method "%s" does not support time', cfg.method);

assert(~hasfreq, 'method "%s" does not support freq', cfg.method);

% set the defaults for method=slice

cfg.nslices = ft_getopt(cfg, 'nslices', 20);

cfg.slicedim = ft_getopt(cfg, 'slicedim', 3);

cfg.slicerange = ft_getopt(cfg, 'slicerange', 'auto');

% ADDED BY JM: allow for slicedim different than 3

switch cfg.slicedim

case 1

if hasana, ana = permute(ana,[2 3 1]); end

if hasfun, fun = permute(fun,[2 3 1]); end

if hasmsk, msk = permute(msk,[2 3 1]); end

cfg.slicedim=3;

case 2

if hasana, ana = permute(ana,[3 1 2]); end

if hasfun, fun = permute(fun,[3 1 2]); end

if hasmsk, msk = permute(msk,[3 1 2]); end

cfg.slicedim=3;

otherwise

% nothing needed

end

%%%%% select slices

if ~ischar(cfg.slicerange)

ind_fslice = cfg.slicerange(1);

ind_lslice = cfg.slicerange(2);

elseif isequal(cfg.slicerange, 'auto')

if hasfun % default

if isfield(functional, 'inside')

insideMask = false(size(fun));

insideMask(functional.inside) = true;

ind_fslice = find(max(max(insideMask,[],1),[],2), 1, 'first');

ind_lslice = find(max(max(insideMask,[],1),[],2), 1, 'last');

else

ind_fslice = find(~isnan(max(max(fun,[],1),[],2)), 1, 'first');

ind_lslice = find(~isnan(max(max(fun,[],1),[],2)), 1, 'last');

end

elseif hasana % if only ana, no fun

ind_fslice = find(max(max(ana,[],1),[],2), 1, 'first');

ind_lslice = find(max(max(ana,[],1),[],2), 1, 'last');

else

ft_error('no functional parameter and no anatomical parameter, can not plot');

end

else

ft_error('do not understand cfg.slicerange');

end

ind_allslice = linspace(ind_fslice,ind_lslice,cfg.nslices);

ind_allslice = round(ind_allslice);

% make new ana, fun, msk, mskana with only the slices that will be plotted (slice dim is always third dimension)

if hasana; new_ana = ana(:,:,ind_allslice); clear ana; ana=new_ana; clear new_ana; end

if hasfun; new_fun = fun(:,:,ind_allslice); clear fun; fun=new_fun; clear new_fun; end

if hasmsk; new_msk = msk(:,:,ind_allslice); clear msk; msk=new_msk; clear new_msk; end

% update the dimensions of the volume

if hasana

dim=size(ana);

else

dim=size(fun);

end

%%%%% make a "quilt", that contain all slices on 2D patched sheet

% Number of patches along sides of Quilt (M and N)

% Size (in voxels) of side of patches of Quilt (m and n)

% take care of a potential singleton 3rd dimension

if numel(dim)<3

dim(end+1:3) = 1;

end

m = dim(1);

n = dim(2);

M = ceil(sqrt(dim(3)));

N = ceil(sqrt(dim(3)));

num_patch = N*M;

num_slice = (dim(cfg.slicedim));

% put empty slides on ana, fun, msk, mskana to fill Quilt up

if hasana; ana(:,:,end+1:num_patch)=0; end

if hasfun; fun(:,:,end+1:num_patch)=0; end

if hasmsk; msk(:,:,end+1:num_patch)=0; end

% if hasmskana; mskana(:,:,end:num_patch)=0; end

% put the slices in the quilt

for iSlice = 1:num_slice

xbeg = floor((iSlice-1)./M);

ybeg = mod(iSlice-1, M);

if hasana

quilt_ana(ybeg*m+1:(ybeg+1)*m, xbeg*n+1:(xbeg+1)*n)=ana(:,:,iSlice);

end

if hasfun

quilt_fun(ybeg*m+1:(ybeg+1)*m, xbeg*n+1:(xbeg+1)*n)=fun(:,:,iSlice);

end

if hasmsk

quilt_msk(ybeg*m+1:(ybeg+1)*m, xbeg*n+1:(xbeg+1)*n)=msk(:,:,iSlice);

end

end

% make vols and scales, containes volumes to be plotted (fun, ana, msk), added by ingnie

if hasana; vols2D{1} = quilt_ana; scales{1} = []; end % needed when only plotting ana

if hasfun; vols2D{2} = quilt_fun; scales{2} = [fcolmin fcolmax]; end

if hasmsk; vols2D{3} = quilt_msk; scales{3} = [opacmin opacmax]; end

% the transpose is needed for displaying the matrix using the MATLAB image() function

if hasana; ana = vols2D{1}'; end

if hasfun && ~doimage; fun = vols2D{2}'; end

if hasfun && doimage; fun = permute(vols2D{2},[2 1 3]); end

if hasmsk; msk = vols2D{3}'; end

if hasana

% scale anatomy between 0 and 1

fprintf('scaling anatomy\n');

amin = min(ana(:));

amax = max(ana(:));

ana = (ana-amin)./(amax-amin);

clear amin amax;

% convert anatomy into RGB values

ana = cat(3, ana, ana, ana);

ha = imagesc(ana);

end

hold on

if hasfun

if doimage

hf = image(fun);

else

hf = imagesc(fun);

try

caxis(scales{2});

end

% apply the opacity mask to the functional data

if hasmsk

% set the opacity

set(hf, 'AlphaData', msk)

set(hf, 'AlphaDataMapping', 'scaled')

try

alim(scales{3});

end

elseif hasana

set(hf, 'AlphaData', 0.5)

end

end

end

axis equal

axis tight

axis xy

axis off

if istrue(cfg.colorbar)

if hasfun

% use a normal MATLAB colorbar

hc = colorbar;

set(hc, 'YLim', [fcolmin fcolmax]);

ylabel(hc, cfg.colorbartext);

else

ft_warning('no colorbar possible without functional data');

cfg.colorbar = 'no'; % prevent the warning from apearing twice

end

end

case 'ortho'

% set the defaults for method=ortho

cfg.location = ft_getopt(cfg, 'location', 'auto');

cfg.locationcoordinates = ft_getopt(cfg, 'locationcoordinates', 'head');

cfg.crosshair = ft_getopt(cfg, 'crosshair', 'yes');

cfg.axis = ft_getopt(cfg, 'axis', 'on');

cfg.queryrange = ft_getopt(cfg, 'queryrange', 1);

if ~ischar(cfg.location)

if strcmp(cfg.locationcoordinates, 'head')

% convert the headcoordinates location into voxel coordinates

loc = inv(functional.transform) * [cfg.location(:); 1];

loc = round(loc(1:3));

elseif strcmp(cfg.locationcoordinates, 'voxel')

% the location is already in voxel coordinates

loc = round(cfg.location(1:3));

else

ft_error('you should specify cfg.locationcoordinates');

end

else

if isequal(cfg.location, 'auto')

if hasfun

if isequal(cfg.funcolorlim, 'maxabs')

loc = 'max';

elseif isequal(cfg.funcolorlim, 'zeromax')

loc = 'max';

elseif isequal(cfg.funcolorlim, 'minzero')

loc = 'min';

else % if numerical

loc = 'max';

end

else

loc = 'center';

end

else

loc = cfg.location;

end

end

% determine the initial intersection of the cursor (xi yi zi)

if ischar(loc) && strcmp(loc, 'min')

if isempty(cfg.funparameter)

ft_error('cfg.location is min, but no functional parameter specified');

end

[dummy, minindx] = min(fun(:));

[xi, yi, zi] = ind2sub(dim, minindx);

elseif ischar(loc) && strcmp(loc, 'max')

if isempty(cfg.funparameter)

ft_error('cfg.location is max, but no functional parameter specified');

end

[dummy, maxindx] = max(fun(:));

[xi, yi, zi] = ind2sub(dim, maxindx);

elseif ischar(loc) && strcmp(loc, 'center')

xi = round(dim(1)/2);

yi = round(dim(2)/2);

zi = round(dim(3)/2);

elseif ~ischar(loc)

% using nearest instead of round ensures that the position remains within the volume

xi = nearest(1:dim(1), loc(1));

yi = nearest(1:dim(2), loc(2));

zi = nearest(1:dim(3), loc(3));

end

if numel(dim)<3