function nbrhood = cosmo_meeg_chan_neighborhood(ds, varargin)
% determine neighborhood of channels in MEEG dataset
%
% nbrhood=cosmo_meeg_chan_neighborhood(ds, ...)
%
% Inputs:
% ds MEEG dataset struct.
% ft_nbrs Optional Fieldtrip-like neighborhood struct with
% fields .label and .neighblabel such as produced by
% FieldTrip's ft_prepare_neighbours. In this case none
% of the other options mentioned below can be
% provided; and the neighborhood structure is based
% on the neighbors from ft_nbrs. The struct ft_nbrs
% must be structured as in this simple example:
% ft_nbrs(1).name='ch1'
% ft_nbrs(1).neighblabel={'ch2','ch6','ch9'}
% ft_nbrs(2).name='ch2'
% ft_nbrs(2).neighblabel={'ch1','ch6'}
% ft_nbrs(3).name='ch6'
% ft_nbrs(3).neighblabel={'ch1','ch2'}
% ft_nbrs(4).name='ch9'
% ft_nbrs(4).neighblabel={'ch1'}
% 'label', lab Optional labels to return in output, one of:
% 'layout' : determine neighbors based on layout
% associated with ds (default). All
% labels in the layout are used as
% center labels.
% 'dataset' : determine neighbors based on labels
% present in ds. Only labels present in
% ds are used as center labels
% {x1,...,xn} : use labels x1 ... xn
% 'chantype', tp (optional) channel type of neighbors, can be one of
% 'eeg', 'meg_planar', 'meg_axial', or
% 'meg_combined_from_planar'.
% Use 'all' to use all channel types associated with
% lab, and 'all_combined' to use
% 'meg_combined_from_planar' with all other channel
% types in ds except for 'meg_planar'.
% If there is only one channel type associated with
% lab, then this argument is not required.
% 'radius', r } select neighbors either within radius r, grow
% 'count', c } the radius to get neighbors at c locations,
% 'delaunay', true } or use Delaunay triangulation to find direct
% } neighbors for each channel.
% } These three options are mutually exclusive
%
% Output:
% nbrhood struct with fields:
% .neighbors Kx1 cell with feature indices of neighbors of
% k-th channel
% .fa.chan channel indices, equal to 1:K
% .a.fdim.values set to a cell with as only element a cell with
% center channel labels
% .a.fdim.labels set to {'chan'}
% .origin.fa } set to the corresponding values from
% .origin.a } the input dataset, ds.fa and ds.a
%
% Examples:
% % (This example requires FieldTrip)
% cosmo_skip_test_if_no_external('fieldtrip');
% %
% % get neighbors at 4 neighboring sensor location for
% % planar neuromag306 channels
% ds=cosmo_synthetic_dataset('type','meeg','size','big');
% nbrhood=cosmo_meeg_chan_neighborhood(ds,...
% 'chantype','meg_planar','count',4);
% cosmo_disp(nbrhood,'edgeitems',1);
% %|| .neighbors
% %|| { [ 2 ... 2e+03 ]@1x28
% %|| :
% %|| [ 149 ... 2.14e+03 ]@1x28 }@204x1
% %|| .fa
% %|| .chan
% %|| [ 1 ... 204 ]@1x204
% %|| .a
% %|| .fdim
% %|| .labels
% %|| { 'chan' }
% %|| .values
% %|| { { 'MEG0113' ... 'MEG2643' }@1x204 }
% %|| .meeg
% %|| .samples_type
% %|| 'timelock'
% %|| .samples_field
% %|| 'trial'
% %|| .samples_label
% %|| 'rpt'
% %|| .origin
% %|| .fa
% %|| .chan
% %|| [ 1 ... 306 ]@1x2142
% %|| .time
% %|| [ 1 ... 7 ]@1x2142
% %|| .a
% %|| .fdim
% %|| .labels
% %|| { 'chan'
% %|| 'time' }
% %|| .values
% %|| { { <char>@1x7 ... <char>@1x7 }@1x306
% %|| [ -0.2 ... 0.1 ]@1x7 }
% %|| .meeg
% %|| .samples_type
% %|| 'timelock'
% %|| .samples_field
% %|| 'trial'
% %|| .samples_label
% %|| 'rpt'
%
% % (This example requires FieldTrip)
% cosmo_skip_test_if_no_external('fieldtrip');
% %
% % get neighbors with radius of .1 for
% % planar neuromag306 channels, but with the center labels
% % the set of combined planar channels
% % (there are 8 channels in the .neighblabel fields, because
% % there are two planar channels per combined channel)
% ds=cosmo_synthetic_dataset('type','meeg','size','big');
% nbrhood=cosmo_meeg_chan_neighborhood(ds,...
% 'chantype','meg_combined_from_planar','radius',.1);
% cosmo_disp(nbrhood,'edgeitems',1);
% %|| .neighbors
% %|| { [ 2 ... 2e+03 ]@1x42
% %|| :
% %|| [ 149 ... 2.14e+03 ]@1x56 }@102x1
% %|| .fa
% %|| .chan
% %|| [ 1 ... 102 ]@1x102
% %|| .a
% %|| .fdim
% %|| .labels
% %|| { 'chan' }
% %|| .values
% %|| { { 'MEG0112+0113' ... 'MEG2642+2643' }@1x102 }
% %|| .meeg
% %|| .samples_type
% %|| 'timelock'
% %|| .samples_field
% %|| 'trial'
% %|| .samples_label
% %|| 'rpt'
% %|| .origin
% %|| .fa
% %|| .chan
% %|| [ 1 ... 306 ]@1x2142
% %|| .time
% %|| [ 1 ... 7 ]@1x2142
% %|| .a
% %|| .fdim
% %|| .labels
% %|| { 'chan'
% %|| 'time' }
% %|| .values
% %|| { { <char>@1x7 ... <char>@1x7 }@1x306
% %|| [ -0.2 ... 0.1 ]@1x7 }
% %|| .meeg
% %|| .samples_type
% %|| 'timelock'
% %|| .samples_field
% %|| 'trial'
% %|| .samples_label
% %|| 'rpt'
%
% % (This example requires FieldTrip)
% cosmo_skip_test_if_no_external('fieldtrip');
% %
% % As above, but combine the two types of channels
% % Here the axial channels only have axial neighbors, and the planar
% % channels only have planar neighbors. With 7 timepoints and 10
% % neighboring channels, the meg_axial channels all have 70 axial
% % neighbors while the meg_planar_combined channels all have
% % 140 neighbors based on the planar channel pairs in the original
% % dataset
% ds=cosmo_synthetic_dataset('type','meeg','size','big');
% nbrhood=cosmo_meeg_chan_neighborhood(ds,...
% 'chantype','all_combined','count',10);
% cosmo_disp(nbrhood,'edgeitems',1);
% %|| .neighbors
% %|| { [ 1 ... 2.07e+03 ]@1x70
% %|| :
% %|| [ 71 ... 2.14e+03 ]@1x140 }@204x1
% %|| .fa
% %|| .chan
% %|| [ 1 ... 204 ]@1x204
% %|| .a
% %|| .fdim
% %|| .labels
% %|| { 'chan' }
% %|| .values
% %|| { { 'MEG0111' ... 'MEG2642+2643' }@1x204 }
% %|| .meeg
% %|| .samples_type
% %|| 'timelock'
% %|| .samples_field
% %|| 'trial'
% %|| .samples_label
% %|| 'rpt'
% %|| .origin
% %|| .fa
% %|| .chan
% %|| [ 1 ... 306 ]@1x2142
% %|| .time
% %|| [ 1 ... 7 ]@1x2142
% %|| .a
% %|| .fdim
% %|| .labels
% %|| { 'chan'
% %|| 'time' }
% %|| .values
% %|| { { <char>@1x7 ... <char>@1x7 }@1x306
% %|| [ -0.2 ... 0.1 ]@1x7 }
% %|| .meeg
% %|| .samples_type
% %|| 'timelock'
% %|| .samples_field
% %|| 'trial'
% %|| .samples_label
% %|| 'rpt'
%
%
% Notes:
%
% See also: cosmo_meeg_chan_neighbors
%
% # For CoSMoMVPA's copyright information and license terms, #
% # see the COPYING file distributed with CoSMoMVPA. #
[dim, unused, attr_name, dim_name, ds_label] = cosmo_dim_find(ds, ...
'chan', true);
nbrs = get_neighbors(ds, varargin{:});
ds_chan = ds.(attr_name).chan;
[nbrs_label, chan_idxs] = get_neighbor_indices(nbrs, ds_label, ds_chan);
nbrhood = struct();
nbrhood.neighbors = chan_idxs;
nbrhood.(attr_name).chan = 1:numel(nbrs);
nbrhood.a = ds.a;
nbrhood.a.(dim_name) = struct();
nbrhood.a.(dim_name).labels = {'chan'};
nbrhood.a.(dim_name).values = {nbrs_label(:)'};
if dim == 2
other_dim_name = 'sdim';
else
other_dim_name = 'fdim';
end
if isfield(nbrhood.a, other_dim_name)
nbrhood.a = rmfield(nbrhood.a, other_dim_name);
end
origin = struct();
origin.(attr_name) = ds.(attr_name);
origin.a = ds.a;
nbrhood.origin = origin;
function [nbrs_label, chan_idxs] = get_neighbor_indices(nbrs, ds_label, ds_chan)
nbrs_label = {nbrs.label};
ds_label_cell = cellfun(@(x){x}, ds_label, 'UniformOutput', false);
ncenter = numel(nbrs_label);
ow = cosmo_overlap({nbrs.neighblabel}, ds_label_cell);
[chan_pos, chan_cell] = cosmo_index_unique({ds_chan});
chan = chan_cell{1};
%
chan_idxs = cell(ncenter, 1);
for k = 1:ncenter
% find the indices in ds where the neighbors match
idx = find(ow(k, :));
chan_pos_msk = cosmo_match(chan, idx);
chan_idx = cat(1, chan_pos{chan_pos_msk});
if isempty(chan_idx)
chan_idx = zeros(0, 1);
end
chan_idxs{k} = chan_idx';
end
function nbrs = get_neighbors(ds, varargin)
if numel(varargin) == 1 && is_neighborhood_struct(varargin{1})
nbrs = varargin{1};
else
nbrs = cosmo_meeg_chan_neighbors(ds, varargin{:});
end
verify_neighbors(nbrs);
function tf = is_neighborhood_struct(nbrs)
tf = isstruct(nbrs) && ...
isfield(nbrs, 'label') && ...
isfield(nbrs, 'neighblabel');
function verify_neighbors(nbrs)
if ~is_neighborhood_struct(nbrs)
error('neighbor struct must be neighborhood struct');
end
