function [montage, cfg] = ft_prepare_montage(cfg, data)

% FT_PREPARE_MONTAGE creates a referencing scheme based on the input configuration

% options and the channels in the data structure. The resulting montage can be

% given as input to FT_APPLY_MONTAGE, or as cfg.montage to FT_PREPROCESSING.

%

% Use as

% montage = ft_prepare_montage(cfg, data)

%

% The configuration can contain the following fields:

% cfg.refmethod = 'avg', 'comp', 'bipolar', 'laplace', 'doublebanana', 'longitudinal', 'circumferential', 'transverse' (default = 'avg')

% cfg.implicitref = string with the label of the implicit reference, or empty (default = [])

% cfg.refchannel = cell-array with new EEG reference channel(s), this can be 'all' for a common average reference

% cfg.groupchans = 'yes' or 'no', should channels be rereferenced in separate groups

% for bipolar and laplace methods, this requires channnels to be

% named using an alphanumeric code, where letters represent the

% group and numbers represent the order of the channel whithin

% its group (default = 'no')

%

% The implicitref option allows adding the implicit reference channel to the data as

% a channel with zeros.

%

% The resulting montage is a structure with the fields

% montage.tra = MxN matrix

% montage.labelold = Nx1 cell-array

% montage.labelnew = Mx1 cell-array

%

% As an example, an output bipolar montage could look like this

% bipolar.labelold = {'1', '2', '3', '4'}

% bipolar.labelnew = {'1-2', '2-3', '3-4'}

% bipolar.tra = [

% +1 -1 0 0

% 0 +1 -1 0

% 0 0 +1 -1

% ];

%

% See also FT_PREPROCESSING, FT_APPLY_MONTAGE

% Copyright (C) 2017-2022, 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 loadvar data

ft_preamble provenance data

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

if ft_abort

return

end

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

hasdata = exist('data', 'var') && ~isempty(data);

% basic check/initialization of input arguments

if ~hasdata

data = struct([]);

else

data = ft_checkdata(data, 'datatype', {'raw+comp', 'timelock+comp', 'raw', 'timelock'});

end

% do a sanity check for incompatible options which are used in ft_preprocessing and elsewhere

cfg = ft_checkconfig(cfg, 'forbidden', {'montage', 'method'});

% set default configuration options

cfg.refmethod = ft_getopt(cfg, 'refmethod', 'avg');

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

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

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

cfg.groupchans = ft_getopt(cfg, 'groupchans', 'no');

if isfield(cfg, 'reref')

assert(istrue(cfg.reref), 'cannot create a montage without cfg.reref=''yes''');

end

% here the actual work starts

if hasdata

cfg.channel = ft_channelselection(cfg.channel, data.label);

cfg.refchannel = ft_channelselection(cfg.refchannel, cat(1, data.label(:), cfg.implicitref));

else

if ischar(cfg.channel)

cfg.channel = {cfg.channel};

end

if ischar(cfg.refchannel)

cfg.refchannel = {cfg.refchannel};

end

end

% the first montage serves to add the implicit reference channel

% the last row for the implicitref will be all zero

montage1 = [];

montage1.labelold = cfg.channel(:);

montage1.labelnew = cfg.channel(:);

if ~isempty(cfg.implicitref)

montage1.labelnew{end+1} = cfg.implicitref;

end

montage1.tra = eye(numel(montage1.labelnew), numel(montage1.labelold));

% the second montage serves to subtract the selected reference channels

switch cfg.refmethod

case 'avg'

% make a montage that subtracts the mean of all channels specified in cfg.refchannel

montage2 = [];

montage2.labelold = montage1.labelnew;

montage2.labelnew = montage1.labelnew;

montage2.tra = eye(numel(montage2.labelnew));

refsel = match_str(montage2.labelold, cfg.refchannel);

montage2.tra(:,refsel) = montage2.tra(:,refsel) - 1/numel(refsel);

% apply montage2 to montage1, the result is the combination of both

montage = ft_apply_montage(montage1, montage2);

case 'comp'

assert(isempty(cfg.implicitref), 'implicitref not supported with refmethod=''%s''', cfg.refmethod);

% construct an linear projection from channels to components

montage = [];

montage.labelold = data.topolabel;

montage.labelnew = data.label;

montage.tra = data.unmixing;

case 'invcomp'

assert(isempty(cfg.implicitref), 'implicitref not supported with refmethod=''%s''', cfg.refmethod);

% construct an linear projection from components to channels

montage = [];

montage.labelold = data.label;

montage.labelnew = data.topolabel;

montage.tra = data.topo;

case 'bipolar'

if istrue(cfg.groupchans)

% this uses recursion to make the montage for each group

label = montage1.labelnew;

group = cell(size(label));

for i=1:numel(label)

group(i) = regexp(label{i}, '^[a-zA-Z]+', 'match');

end

ugroup = unique(group, 'stable');

tmpcfg = keepfields(cfg, {'refmethod', 'implicitref', 'refchannel', 'channel'});

tmpcfg.groupchans = 'no';

tmpcfg.showcallinfo = 'no';

groupmontage = cell(size(ugroup));

for g = 1:length(ugroup) % for each group of channels

% select the channels within the group

tmpcfg.channel = label(strcmp(group, ugroup{g}));

groupmontage{g} = ft_prepare_montage(tmpcfg);

end

% start with a montage that keeps all channels the same

montage2 = [];

montage2.labelold = label(:)';

montage2.labelnew = label(:);

montage2.tra = eye(length(label));

for g = length(ugroup):-1:1

% apply the montage for each groep (i.e., combine them), keep all other channels as they are

% unused channels end up at the end, hence doing it in reverse order keeps the channel order consistent

montage2 = ft_apply_montage(montage2, groupmontage{g}, 'keepunused', true);

end

else

% create a montage for the bipolar derivation of sequential channels

montage2 = [];

montage2.labelold = montage1.labelnew;

montage2.labelnew = strcat(montage1.labelnew(1:end-1),'-',montage1.labelnew(2:end));

tra_neg = diag(-ones(numel(montage1.labelnew)-1,1), 1);

tra_plus = diag( ones(numel(montage1.labelnew)-1,1),-1);

montage2.tra = tra_neg(1:end-1,:)+tra_plus(2:end,:);

end % if groupchans

% apply montage2 to montage1, the result is the combination of both

montage = ft_apply_montage(montage1, montage2);

case 'laplace'

if istrue(cfg.groupchans)

% this uses recursion to make the montage for each group

label = montage1.labelnew;

group = cell(size(label));

for i=1:numel(label)

group(i) = regexp(label{i}, '^[a-zA-Z]+', 'match');

end

group = unique(group, 'stable');

tmpcfg = keepfields(cfg, {'refmethod', 'implicitref', 'refchannel', 'channel'});

tmpcfg.groupchans = 'no';

tmpcfg.showcallinfo = 'no';

groupmontage = cell(size(group));

for g = 1:length(group) % for each group of channels

% select the channels within the group

tmpcfg.channel = label(startsWith(label, group{g}));

groupmontage{g} = ft_prepare_montage(tmpcfg);

end

% start with a montage that keeps all channels the same

montage2 = [];

montage2.labelold = label(:)';

montage2.labelnew = label(:);

montage2.tra = eye(length(label));

for g = length(group):-1:1

% apply the montage for each groep (i.e., combine them), keep all other channels as they are

% unused channels end up at the end, hence doing it in reverse order keeps the channel order consistent

montage2 = ft_apply_montage(montage2, groupmontage{g}, 'keepunused', true);

end

else

% create a montage for the laplacian derivation of neighboring channels where

% each channel is re-referenced against the mean of the two nearest channels,

% while channels at the extremities will be re-referenced against their closest

% neighbor; this is particularly useful for sEEG shafts

montage2 = [];

montage2.labelold = montage1.labelnew;

tra_neg1 = diag(-0.5*ones(numel(montage1.labelnew)-1,1), 1);

tra_neg2 = diag(-0.5*ones(numel(montage1.labelnew)-1,1), -1);

if numel(montage2.labelold) > 2

montage2.labelnew = montage1.labelnew;

tra_plus = diag(ones(numel(montage1.labelnew),1),0);

montage2.tra = tra_neg1+tra_neg2+tra_plus;

montage2.tra(1,2) = -1;

montage2.tra(end,end-1) = -1;

else

% The following commands are just meant to avoid errors, in case

% only one channel is given as input: the algorithm will behave

% exactly as it would in the 'bipolar' case

montage2.labelnew = cell(numel(montage1.labelnew),0);

tra_plus = diag(ones(numel(montage1.labelnew)-1,1),-1);

montage2.tra = tra_neg1+tra_neg2+tra_plus(2:end,:);

end

end % if groupchans

% apply montage2 to montage1, the result is the combination of both

montage = ft_apply_montage(montage1, montage2);

case {'longitudinal', 'doublebanana'}

% see https://www.learningeeg.com/montages-and-technical-components

montage = [];

montage.labelold = {'Fp1', 'Fpz', 'Fp2', 'F7', 'F3', 'Fz', 'F4', 'F8', 'T3', 'C3', 'Cz', 'C4', 'T4', 'T5', 'P3', 'Pz', 'P4', 'T6', 'O1', 'Oz', 'O2'};

montage.labelnew = {

% left temporal chain

'Fp1-F7'

'F7-T3'

'T3-T5'

'T5-O1'

% left parasagittal chain

'Fp1-F3'

'F3-C3'

'C3-P3'

'P3-O1'

% central chain

'Fz-Cz'

'Cz-Pz'

% right parasagittal chain

'Fp2-F4'

'F4-C4'

'C4-P4'

'P4-O2'

% right temporal chain

'Fp2-F8'

'F8-T4'

'T4-T6'

'T6-O2'

};

% construct the montage from new channel labels

montage.tra = zeros(length(montage.labelnew), length(montage.labelold));

for i=1:length(montage.labelnew)

lab = split(montage.labelnew{i}, '-');

montage.tra(i, strcmp(montage.labelold, lab{1})) = +1;

montage.tra(i, strcmp(montage.labelold, lab{2})) = -1;

end

% do a sanity check on the montage

assert(all(sum(montage.tra, 2)==0));

assert(all(max(montage.tra, [], 2)==+1));

assert(all(min(montage.tra, [], 2)==-1));

case 'circumferential'

% see https://www.learningeeg.com/montages-and-technical-components

montage = [];

montage.labelold = {'Fp1', 'Fpz', 'Fp2', 'F7', 'F3', 'Fz', 'F4', 'F8', 'T3', 'C3', 'Cz', 'C4', 'T4', 'T5', 'P3', 'Pz', 'P4', 'T6', 'O1', 'Oz', 'O2'};

montage.labelnew = {

'Fp1-F7'

'F7-T3'

'T3-T5'

'T5-O1'

'O1-O2'

'O2-T6'

'T6-T4'

'T4-F8'

'F8-Fp2'

'Fp2-Fp1'

};

% construct the montage from new channel labels

montage.tra = zeros(length(montage.labelnew), length(montage.labelold));

for i=1:length(montage.labelnew)

lab = split(montage.labelnew{i}, '-');

montage.tra(i, strcmp(montage.labelold, lab{1})) = +1;

montage.tra(i, strcmp(montage.labelold, lab{2})) = -1;

end

% do a sanity check on the montage

assert(all(sum(montage.tra, 2)==0));

assert(all(max(montage.tra, [], 2)==+1));

assert(all(min(montage.tra, [], 2)==-1));

case 'transverse'

% this is inspired by https://doi.org/10.1016/j.earlhumdev.2011.08.008 but probably better documented elsewhere

% see https://www.learningeeg.com/montages-and-technical-components

% this particular implementation only considers the 10% distances

montage = [];

montage.labelold = {'Fp1', 'Fpz', 'Fp2', 'F7', 'F3', 'Fz', 'F4', 'F8', 'T3', 'C3', 'Cz', 'C4', 'T4', 'T5', 'P3', 'Pz', 'P4', 'T6', 'O1', 'Oz', 'O2'};

montage.labelnew = {

'Fp1-Fp2'

'F7-F3'

'F3-Fz'

'Fz-F4'

'F4-F8'

'T3-C3'

'C3-Cz'

'Cz-C4'

'C4-T4'

'T5-P3'

'P3-Pz'

'Pz-P4'

'P4-T6'

'O1-O2'

};

% construct the montage from new channel labels

montage.tra = zeros(length(montage.labelnew), length(montage.labelold));

for i=1:length(montage.labelnew)

lab = split(montage.labelnew{i}, '-');

montage.tra(i, strcmp(montage.labelold, lab{1})) = +1;

montage.tra(i, strcmp(montage.labelold, lab{2})) = -1;

end

% do a sanity check on the montage

assert(all(sum(montage.tra, 2)==0));

assert(all(max(montage.tra, [], 2)==+1));

assert(all(min(montage.tra, [], 2)==-1));

otherwise

error('unsupported refmethod=''%s''', cfg.refmethod);

end

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

ft_postamble provenance

ft_postamble previous data

ft_postamble history montage