function [data_vc] = ft_virtualchannel(cfg, data, source, parcellation)

% FT_VIRTUALCHANNEL creates virtual channel data, combining numeric data from a data

% structure defined at the channel level with spatial filter information from a

% source data structure, and optional parcellation information.

%

% Use as

% output = ft_virtualchannel(cfg, data, source)

% or

% output = ft_virtualchannel(cfg, data, source, parcellation)

%

% where the input "data" is a channel-level data structure that can be linearly

% mapped onto the virtual channel level, e.g. a raw data structure obtained with

% FT_PREPROCESSING, a timelock structure, obtained with FT_TIMELOCKANALYSIS, or a

% freq structure with fourierspectra, obtained with FT_FREQANALYSIS.

%

% The input "source" is a source structure that has been obtained with

% FT_SOURCEANALYSIS, and which contains spatial filter information for at least one

% dipole location, in the source.filter, or source.avg.filter field.

%

% The optional input "parcellation" is described in detail in

% FT_DATATYPE_PARCELLATION (2-D) or FT_DATATYPE_SEGMENTATION (3-D) and can be

% obtained from FT_READ_ATLAS or from a custom parcellation/segmentation for your

% individual subject. Alternatively, the input "source" can already contain a

% parcellation.

%

% The configuration "cfg" is a structure that should either contain

% cfg.pos = Nx3 matrix containing the dipole positions for the virtual

% channel(s). These positions should match the entries in

% the source.pos field. (default = [])

% or

% cfg.parcellation = string, name of the field that is used for the

% parcel labels. (default = [])

% cfg.parcel = string, or cell-array of strings, specifying for which

% parcels to return the output. (default = 'all')

%

% Moreover, the cfg structure can contain

% cfg.method = string, determines how the components of the specified virtual

% channel(s) are to to be combined. 'svd' (default), 'none', 'pca',

% 'runica', 'fastica', 'dss'.

% cfg.numcomponent = scalar (or 'all'), determines the number of components per virtual

% channel in the output. (default = 1)

%

% See also FT_SOURCEANALYSIS, FT_DATATYPE_PARCELLATION, FT_DATATYPE_SEGMENTATION,

% FT_SOURCEPARCELLATE, FT_COMPONENTANALYSIS

% Copyright (C) 2020, Jan-Mathijs Schoffelen and Robert Oostenveld

%

% 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 data source parcellation

ft_preamble provenance data source parcellation

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

if ft_abort

return

end

% FIXME for now support only raw-type data structures, this should be extended to

% also allow freq (with Fourier) and timelock structures, although the latter is

% probably already covered by the call to checkdata, it should be checked how this is

% done in ft_megplanar.

% store the original input representation of the data, this is used later on to convert it back

isfreq = ft_datatype(data, 'freq');

istlck = ft_datatype(data, 'timelock'); % this will be temporary converted into raw

data = ft_checkdata(data, 'datatype', {'raw+comp', 'raw', 'mvar' 'freq'}, 'feedback', 'yes');

% ensure that the source input is a source structure , not a volume structure

% this will also return source.filter, rather than source.avg.filter

source = ft_checkdata(source, 'datatype', 'source', 'insidestyle', 'logical', 'hasunit', 'yes');

% get the defaults

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

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

cfg.parcel = ft_getopt(cfg, 'parcel', 'all');

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

cfg.feedback = ft_getopt(cfg, 'feedback', 'text');

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

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

hasparcellation = exist('parcellation', 'var');

if hasparcellation

% the parcellation is specified as separate structure

elseif ~isempty(cfg.parcellation)

% the parcellation is represented in the source structure itself

parcellation = source;

end

useparcellation = ~isempty(cfg.parcellation) && isfield(parcellation, cfg.parcellation);

usepos = ~isempty(cfg.pos);

if usepos && useparcellation

ft_error('you should either specify cfg.pos, or cfg.parcellation, not both');

elseif usepos

% extract the virtual channels from the specified positions

elseif useparcellation

% extract the virtual channels from the specified parcels

else

ft_error('you should either specify cfg.pos, or a valid cfg.parcellation');

end

if ~isfield(source, 'filter')

ft_error('the input source needs a ''filter'' field');

% FIXME how about the output of ft_dipolefitting?

end

if isfreq

% some restrictions apply to frequency data

if ~isfield(data, 'fourierspctrm'), ft_error('freq data should contain Fourier spectra'); end

%if numel(data.freq)>1, ft_error('with spectral input data only a single frequency bin is supported'); end

if ~any(strcmp({'svd' 'none'}, cfg.method)), ft_error('with spectral input data only ''svd'' or ''none'' are supported methods'); end

end

% select channels and trials of interest, by default this will select all channels and trials

[i1, i2] = match_str(data.label, source.label);

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

tmpcfg.channel = data.label(i1);

data = ft_selectdata(tmpcfg, data);

% restore the provenance information

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

% check and fix the channel order in the spatial filters

if ~isequal(i2, (1:numel(source.label))')

ft_warning('the order of the channels in the spatialf filters is different from the order of channels in the data, reordering\n');

for i = find(source.inside(:)')

source.filter{i} = source.filter{i}(:, i2);

end

source.label = source.label(i2);

end

% do the actual work

if usepos

% identify the indices of the positions in the source structure that are

% to be used, and check whether a spatial filter exists

npos = size(source.pos,1);

nvc = size(cfg.pos,1);

indx = zeros(nvc,1);

mindist = zeros(nvc,1);

for i = 1:nvc

dpos = sqrt(sum( (source.pos - cfg.pos(i.*ones(npos,1),:)).^2, 2));

[mindist(i), indx(i)] = min(dpos);

% check that the requested positions are at most 1 mm away from the actual

% positions

if mindist(i)>1*ft_scalingfactor('mm', source.unit)

ft_error('the requested dipole position is > 1 mm away from the closest source position');

end

% check that the spatial filter exists

if isempty(source.filter{indx(i)})

ft_error(sprintf('the spatial filter information is missing for dipole with position [%d %d %d]', cfg.pos(i,:)));

end

end

bname = 'virtualchannel'; % will be used as name for the balancing scheme for the sensor array

elseif useparcellation

% do some checks and balances on the input data and cfg

% keep the transformation matrix

if isfield(parcellation, 'transform')

transform = parcellation.transform;

else

transform = [];

end

% ensure it is a parcellation, not a segmentation

parcellation = ft_checkdata(parcellation, 'datatype', 'parcellation', 'parcellationstyle', 'indexed');

% keep the transformation matrix

if ~isempty(transform)

parcellation.transform = transform;

end

% ensure that the source and the parcellation are anatomically consistent

tolerance = 0.1 * ft_scalingfactor('mm', source.unit);

if ~isalmostequal(source.pos, parcellation.pos, 'abstol', tolerance)

ft_error('the source positions are not consistent with the parcellation, please use FT_SOURCEINTERPOLATE');

end

if isequal(cfg.parcel, 'all')

cfg.parcel = parcellation.(sprintf('%slabel',cfg.parcellation));

end

indx = match_str(parcellation.(sprintf('%slabel',cfg.parcellation)), cfg.parcel);

cfg.parcel = parcellation.(sprintf('%slabel',cfg.parcellation))(indx);

nvc = numel(cfg.parcel);

bname = cfg.parcellation;

end

% Create a montage for each of the dipoles/parcels, and use ft_apply_montage,

% on the data, followed by a dimensionality reduction step, using

% ft_componentanalysis. Also, keep in mind that this is essentially a

% two-step montage application, so keep the individual steps in mind, to be

% able to combine it also to the sensor description (which will broadcast

% the unmixing to the output data structure, so that it can be used later

% on as well

unmixing1 = cell(1, nvc);

unmixing2 = cell(1, nvc);

label_in = cell(1, nvc);

label_out = cell(1, nvc);

compmethods = {'pca' 'runica' 'fastica' 'dss'}; % FIXME in principle all methods supported in ft_componentanalysis should work

for i = 1:nvc

montage = [];

if usepos

montage.tra = source.filter{indx(i)};

montage.labelold = source.label;

montage.labelnew = cell(size(montage.tra,1), 1);

for k = 1:size(montage.tra,1)

montage.labelnew{k} = sprintf('virtualchannel%03d_orientation%03d', i, k);

end

elseif useparcellation

sel = find(parcellation.(cfg.parcellation)==indx(i));

montage.tra = cat(1, source.filter{sel});

montage.labelold = source.label;

montage.labelnew = cell(size(montage.tra,1), 1);

cnt = 0;

for k = 1:numel(sel)

for kk = 1:size(source.filter{sel(k)},1)

cnt = cnt+1;

montage.labelnew{cnt} = sprintf('%s_dipole%03d_orientation%03d', cfg.parcel{i}, k, kk);

end

end

end

unmixing1{1,i} = montage.tra;

label_in{1,i} = data.label;

% Next, apply the montage to the numeric data, but make an exception for the

% svd method: this one is much more efficient, in that its implementation

% operates on the covariance matrix directly, bypasses ft_componentanalysis,

% and does not require the intermediate (possibly memory greedy, and slow)

% data representation.

switch cfg.method

case 'svd'

if ~exist('C', 'var')

% create a covariance, or csd matrix that is to be used for the

% svd. this needs to be computed only once.

if isfreq

tmp = ft_checkdata(data, 'cmbstyle', 'fullfast');

C = tmp.crsspctrm;

if ndims(C)>2

siz = size(C);

C = mean(reshape(C, [siz(1) siz(2) prod(siz(3:end))]),3);

end

else

tmpcfg = [];

tmpcfg.covariance = 'yes';

tlck = ft_timelockanalysis(tmpcfg, data);

C = tlck.cov;

end

% create a matricial square root of the full covariance/csd matrix

[u, s, v] = svd(real(C));

Csqrtm = u*sqrt(s);

end

% don't do the sandwiching for efficiency, we will use only the u

% matrix, this yields the same u matrix as the svd on the tlck.cov

[u, s, v] = svd(unmixing1{1, i} * Csqrtm, 'econ');

if isequal(cfg.numcomponent, 'all')

ncomp = size(u,2);

elseif ischar(cfg.numcomponent) && isequal(cfg.numcomponent(end), '%')

% return number of components that explain x% of variance

perc = sscanf(cfg.numcomponent, '%f%%')./100;

diags = diag(s);

ncomp = find(cumsum(diags)./sum(diags) > perc, 1, 'first');

else

if size(u,2)>=cfg.numcomponent

ncomp = cfg.numcomponent;

else

ncomp = size(u,2);

ft_warning(sprintf('using %d components,rather than the requested %d\n', ncomp, cfg.numcomponent));

end

end

unmixing2{1, i} = u(:,1:ncomp)';

label_out{1, i} = cell(ncomp,1);

if usepos

str = sprintf('virtualchannel%03d', i);

else

str = cfg.parcel{i};

end

if ncomp == 1

label_out{1, i}{1} = str;

else

for k = 1:ncomp

label_out{1, i}{k} = sprintf('%s_%03d', str, k);

end

end

case 'none'

% do nothing, i.e. stick to the original 'montage', which is the

% 3D, or concatenated spatial filter

unmixing2{1, i} = eye(size(unmixing1{1, i}, 1));

label_out{1, i} = montage.labelnew;

case compmethods

% apply the montage to the data

tmpdata = ft_apply_montage(data, montage, 'feedback', 'none');

tmpcfg = keepfields(cfg, {'method' cfg.method 'numcomponent' 'cellmode' 'demean' 'doscale'});

if isequal(tmpcfg.numcomponent, 'all')

tmpcfg = rmfield(tmpcfg, 'numcomponent');

end

tmpdata = ft_componentanalysis(tmpcfg, tmpdata);

unmixing2{i} = tmpdata.unmixing;

for k = 1:numel(tmpdata.label)

if usepos

str = sprintf('virtualchannel%03d', i);

else

str = cfg.parcel{i};

end

if numel(tmpdata.label)==1

tmpdata.label{k} = sprintf('%s', str);

else

tmpdata.label{k} = sprintf('%s_%03d', str, k);

end

end

label_out{1, i} = tmpdata.label;

otherwise

ft_error('currently not yet supported');

% the idea could be to support a custom function here, with a

% function(cfg.(cfg.method), tmpdata{i}) API

end % reduction of components

end % for i = # of virtual channels

% create and apply the compound montage

unmixing = zeros(0, size(unmixing1{1},2));

for i = 1:numel(unmixing1)

unmixing = cat(1, unmixing, unmixing2{i} * unmixing1{i});

end

montage = [];

montage.tra = unmixing;

montage.labelnew = cat(1, label_out{:});

montage.labelold = label_in{1};

data_vc = ft_apply_montage(data, montage, 'feedback', 'none');

% apply the montage to the sensor description

sensfields = {'grad' 'elec' 'opto'};

for k = 1:numel(sensfields)

if isfield(data_vc, sensfields{k})

ft_info(sprintf('applying the montage to the %s structure\n', sensfields{k}));

data_vc.(sensfields{k}) = ft_apply_montage(data.(sensfields{k}), montage, 'feedback', 'none', 'keepunused', 'yes', 'balancename', bname);

end

end

if istlck

% convert the raw structure back into a timelock structure

data_vc = ft_checkdata(data_vc, 'datatype', 'timelock');

end

if usepos

brainordinate = keepfields(source, {'pos', 'tri', 'dim', 'transform', 'unit', 'coordsys', 'inside'});

brainordinate.index = zeros(npos, 1);

brainordinate.index(indx) = indx;

brainordinate.indexlabel = data_vc.label; % FIXME this only works with one component per vc

elseif useparcellation

brainordinate = parcellation;

end

data_vc.brainordinate = brainordinate;

% do the general cleanup and bookkeeping at the end of the function

ft_postamble debug

ft_postamble previous data source parcellation

ft_postamble provenance data_vc

ft_postamble history data_vc

ft_postamble savevar data_vc