ft_prepare_neighbours.m

function [neighbours, cfg] = ft_prepare_neighbours(cfg, data)

% FT_PREPARE_NEIGHBOURS finds the channel neighbours for spatial clustering
% or interpolation of bad channels. Using the 'distance' method, neighbours
% are based on a minimum neighbourhood distance (in cfg.neighbourdist).
% Using the 'triangulation' method calculates a triangulation based on a 2D
% projection of the sensor positions. The 'template' method loads a default
% template for the given data type. Alternatively, using the 'parcellation'
% method, in combination with an atlas as input data, spatial neighbours
% of parcels are determined, based on the spatial relationship between the
% labeled mesh vertices. Currently, only atlases defined on a triangular
% mesh are supported.
%
% Use as
%   neighbours = ft_prepare_neighbours(cfg)
% or
%   neighbours = ft_prepare_neighbours(cfg, data)
% with an input data structure with the channels of interest and that
% contains a sensor description, or represents an atlas, see FT_READ_ATLAS
%
% The configuration can contain
%   cfg.channel       = channels in the data for which neighbours should be determined
%   cfg.method        = 'distance', 'triangulation' or 'template'
%   cfg.template      = name of the template file, e.g. CTF275_neighb.mat
%   cfg.neighbourdist = number, maximum distance between neighbouring sensors
%                       (only for 'distance', default is 40 mm)
%   cfg.compress      = 'yes' or 'no', add extra edges by compressing in the
%                       x- and y-direction (only for 'triangulation', default is yes)
%   cfg.feedback      = 'yes' or 'no' (default = 'no')
%
% The 3D sensor positions can be present in the data or can be specified as
%   cfg.elec          = structure with electrode positions or filename, see FT_READ_SENS
%   cfg.grad          = structure with gradiometer definition or filename, see FT_READ_SENS
%
% The 2D channel positions can be specified as
%   cfg.layout        = filename of the layout, see FT_PREPARE_LAYOUT
%
% With an atlas in the input, the method 'parcellation' has the additional
% options
%   cfg.parcellation  = string that denotes the field in the atlas that is to be used
%
% The output is an array of structures with the "neighbours" which is
% structured like this:
%        neighbours(1).label = 'Fz';
%        neighbours(1).neighblabel = {'Cz', 'F3', 'F3A', 'FzA', 'F4A', 'F4'};
%        neighbours(2).label = 'Cz';
%        neighbours(2).neighblabel = {'Fz', 'F4', 'RT', 'RTP', 'P4', 'Pz', 'P3', 'LTP', 'LT', 'F3'};
%        neighbours(3).label = 'Pz';
%        neighbours(3).neighblabel = {'Cz', 'P4', 'P4P', 'Oz', 'P3P', 'P3'};
%        etc.
%
% Note that a channel is not considered to be a neighbour of itself.
%
% See also FT_NEIGHBOURPLOT, FT_PREPARE_LAYOUT, FT_DATATYPE_SENS,
% FT_READ_SENS, FT_READ_ATLAS

% Copyright (C) 2006-2021, Eric Maris, Jorn M. Horschig, Robert Oostenveld,
% 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    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 arguments or can be read from disk
hasdata  = exist('data', 'var');
hasatlas = hasdata && (ft_datatype(data, 'mesh+label') || ft_datatype(data, 'volume+label') || ft_datatype(data, 'source+label'));

% these undocumented methods are needed to support some of the high-level FT functions that call this function
if ~isfield(cfg, 'method')
  if ~isfield(cfg, 'neighbours')
    ft_error('cannot figure out how to construct neighbours, please specify cfg.neighbours or cfg.method and call this function directly')
  else
    if ischar(cfg.neighbours)
      ft_notice('reading neighbours from file %s', cfg.neighbours);
      cfg.neighbours = loadvar(cfg.neighbours);
    elseif isstruct(cfg.neighbours) && ~isempty(cfg.neighbours)
      ft_notice('using specified neighbours for the channels');
    elseif isempty(cfg.neighbours) && ~isfield(cfg, 'neighbourdist')
      ft_notice('not using neighbours for the channels');
      % make an empty neighbour structure
      tmp = struct('label', [], 'neighblabel', []);
      cfg.neighbours = tmp([]);
    end
    % this is a hack to get past the case-statement
    cfg.method = 'specified';
  end
end

if hasatlas && ~strcmp(cfg.method, 'parcellation')
  ft_error('an atlas as input data argument requires cfg.method = ''parcellation''');
elseif strcmp(cfg.method, 'parcellation') && ~hasatlas
  ft_error('cfg.method = ''parcellation'' requires an atlas as input data');
end

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'forbidden',  {'channels'}); % prevent accidental typos, see issue 1729
cfg = ft_checkconfig(cfg, 'required',   {'method'});
cfg = ft_checkconfig(cfg, 'renamed',    {'elecfile', 'elec'});
cfg = ft_checkconfig(cfg, 'renamed',    {'gradfile', 'grad'});
cfg = ft_checkconfig(cfg, 'renamed',    {'optofile', 'opto'});
cfg = ft_checkconfig(cfg, 'renamedval', {'method', 'tri', 'triangulation'});

% set the defaults
cfg.feedback     = ft_getopt(cfg, 'feedback',     'no');
cfg.channel      = ft_getopt(cfg, 'channel',      'all');
cfg.compress     = ft_getopt(cfg, 'compress',     'yes');
cfg.parcellation = ft_getopt(cfg, 'parcellation', 'parcellation');

if hasdata && ~hasatlas
  % check if the input data is valid for this function
  data = ft_checkdata(data);
  % set the default for senstype depending on the data
  if isfield(data, 'grad')
    cfg.senstype = ft_getopt(cfg, 'senstype', 'meg');
  elseif isfield(data, 'elec')
    cfg.senstype = ft_getopt(cfg, 'senstype', 'eeg');
  elseif isfield(data, 'opto')
    cfg.senstype = ft_getopt(cfg, 'senstype', 'opto');
  else
    cfg.senstype = ft_getopt(cfg, 'senstype', []);
  end
end

if strcmp(cfg.method, 'distance') || strcmp(cfg.method, 'triangulation')
  % these methods require channel positions in either 3D or in 2D
  
  if isfield(cfg, 'layout') && ~isempty(cfg.layout)
    % get 2D positions from the layout
    tmpcfg  = keepfields(cfg, {'layout', 'channel',   'rows', 'columns', 'commentpos', ...
               'skipcomnt', 'scalepos',  'skipscale', 'projection', 'viewpoint', ...
               'rotate',    'width',     'height',    'elec', 'grad', 'opto', 'showcallinfo', 'trackcallinfo', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo', 'checksize'});
    tmpcfg.skipscale = 'yes';
    tmpcfg.skipcomnt = 'yes';
    layout  = ft_prepare_layout(tmpcfg);
    chanpos = layout.pos;
    label   = layout.label;
  else
    % get 3D positions from the sensor description
    if hasdata
      sens = ft_fetch_sens(cfg, data);
    else
      sens = ft_fetch_sens(cfg);
    end
    chanpos = sens.chanpos;
    label   = sens.label;
  end
  
  if hasdata
    % remove channels that are not in data
    [dum, sensidx] = match_str(data.label, label);
    chanpos = chanpos(sensidx, :);
    label   = label(sensidx);
  end
  
  % select the desired channels
  desired = ft_channelselection(cfg.channel, label);
  [sensidx] = match_str(label, desired);
  chanpos = chanpos(sensidx, :);
  label   = label(sensidx);
  
  if strcmp(ft_senstype(label), 'neuromag306')
    ft_warning('Please be aware of the different sensor types in neuromag306 system, see http://www.fieldtriptoolbox.org/faq/why_are_there_multiple_neighbour_templates_for_the_neuromag306_system');
  end
  
end % if distance or triangulation

switch cfg.method
  case 'specified'
    % use the neighbours as specified by the user
    neighbours = cfg.neighbours;
    
  case 'template'
    fprintf('Trying to load sensor neighbours from a template\n');
    % determine from where to load the neighbour template
    if ~isfield(cfg, 'template')
      % if data has been put in, try to estimate the sensor type
      if hasdata
        fprintf('Estimating sensor type of data to determine the layout filename\n');
        senstype = ft_senstype(data.label);
        fprintf('Data is of sensor type ''%s''\n', senstype);
        if ~exist([senstype '_neighb.mat'], 'file')
          if exist([senstype '.lay'], 'file')
            cfg.layout = [senstype '.lay'];
          else
            fprintf('Name of sensor type does not match name of layout- and template-file\n');
          end
        else
          cfg.template = [senstype '_neighb.mat'];
        end
      end
    end
    % if that failed
    if ~isfield(cfg, 'template')
      % check whether a layout can be used
      if ~isfield(cfg, 'layout')
        % error if that fails as well
        ft_error('You need to define a template or layout or give data as an input argument when ft_prepare_neighbours is called with cfg.method=''template''');
      end
      fprintf('Using the 2-D layout filename to determine the template filename\n');
      cfg.template = [strtok(cfg.layout, '.') '_neighb.mat'];
    end
    % adjust filename
    if ~exist(cfg.template, 'file')
      cfg.template = lower(cfg.template);
    end
    % add necessary extensions
    if numel(cfg.template) < 4 || ~isequal(cfg.template(end-3:end), '.mat')
      if numel(cfg.template) < 7 || ~isequal(cfg.template(end-6:end), '_neighb')
        cfg.template = [cfg.template, '_neighb'];
      end
      cfg.template = [cfg.template, '.mat'];
    end
    % check for existence
    if ~exist(cfg.template, 'file')
      ft_error('Template file could not be found - please check spelling or see http://www.fieldtriptoolbox.org/faq/how_can_i_define_my_own_neighbourhood_template (please consider sharing it with others via the FT mailing list)');
    end
    load(cfg.template);
    fprintf('Successfully loaded neighbour structure from %s\n', cfg.template);
    
  case 'distance'
    % use a smart default for the distance
    if ~isfield(cfg, 'neighbourdist')
      if exist('sens', 'var')
        sens = ft_determine_units(sens);
        cfg.neighbourdist = 40 * ft_scalingfactor('mm', sens.unit);
      elseif exist('layout', 'var')
        layout = ft_determine_units(layout);
        cfg.neighbourdist = 40 * ft_scalingfactor('mm', layout.unit);
      end
      fprintf('using a distance threshold of %g\n', cfg.neighbourdist);
    end
    neighbours = compneighbstructfrompos(chanpos, label, cfg.neighbourdist);
    
  case 'triangulation'
    if size(chanpos, 2)==2 || all(chanpos(:,3)==0)
      % the sensor positions are already on a 2D plane
      prj = chanpos(:,1:2);
    else
      % project sensor on a 2D plane
      prj = elproj(chanpos);
    end
    % make a 2d delaunay triangulation of the projected points
    tri = delaunay(prj(:,1), prj(:,2));
    if strcmp(cfg.compress, 'yes')
      tri_x = delaunay(prj(:,1)./2, prj(:,2)); % compress in the x-direction
      tri_y = delaunay(prj(:,1), prj(:,2)./2); % compress in the y-direction
      tri = [tri; tri_x; tri_y];
    end
    neighbours = compneighbstructfromtri(chanpos, label, tri);
    
  case 'parcellation'
    
    if ~isfield(data, cfg.parcellation)
      ft_error('required field %s not present in the input atlas', cfg.parcellation);
    end
    
    switch ft_datatype(data)
      case 'mesh+label'
        
        tri = data.tri;
        
        % ensure that the vertices are indexed starting from 1
        if min(tri(:))==0
          tri = tri + 1;
        end
        
        % ensure that the vertices are indexed according to 1:number of unique vertices
        tri = tri_reindex(tri);
        
        % create the unique edges from the triangulation
        edges  = [tri(:,[1 2]); tri(:,[1 3]); tri(:,[2 3])];
        edges  = double(unique(sort([edges; edges(:,[2 1])],2), 'rows'));
        
        % get the parcel values for the edges that 'go across parcels'
        boundary = data.(cfg.parcellation)(edges);
        boundary = boundary(boundary(:,1)~=boundary(:,2),:);
        
        % unique parcel-crossing edges
        uboundary = unique(boundary, 'rows');
        %uboundary = unique([boundary; boundary(:,[2 1])], 'rows');
        
        % count
        nboundary = zeros(size(uboundary, 1), 1);
        for k = 1:numel(nboundary)
          nboundary(k) = sum(boundary(:,1)==uboundary(k,1) & boundary(:,2)==uboundary(k,2));
        end
        
        nthresh   = 2; % at least nthresh edges need to be there, in order to count as a neighbour, e.g.: touching at the corner does not count
        uboundary = uboundary(nboundary>nthresh, :);
        uboundary = [uboundary; uboundary(:,[2 1])]; % to make the adjacency matrix symmetric
        
        % fill the adjacency matrix
        n   = size(uboundary,1);
        adj = full(sparse(uboundary(:,1),uboundary(:,2),ones(n,1)));
        
        label = data.([cfg.parcellation 'label']);
        
        neighbours = struct([]);
        for k = 1:size(adj,1)
          neighbours(k).label       = label{k};
          neighbours(k).neighblabel = label(adj(k,:)>0);
        end
        
      case {'volume+label' 'source+label'}
        ft_error('not yet implemented');
      otherwise
        % this shouldn't happen, but throw an error just to be sure
        ft_error('the input data should be an atlas with cfg.method = ''parcellation''');
    end
    
  otherwise
    ft_error('unsupported method ''%s''', cfg.method);
end

% only select those channels that are in the data
if isfield(cfg, 'channel') && ~isempty(cfg.channel)
  if hasdata && ~hasatlas
    desired = ft_channelselection(cfg.channel, data.label);
  else
    desired = ft_channelselection(cfg.channel, {neighbours(:).label});
  end
elseif hasdata
  desired = data.label;
else
  desired = {neighbours(:).label};
end

if ~isempty(desired)
  complete = struct;
  for i = 1:numel(desired)
    complete(i).label = desired{i};
    sel = find(strcmp({neighbours(:).label}, desired{i}));
    if numel(sel)==1
      % take the set of neighbours from the definition
      complete(i).neighblabel = neighbours(sel).neighblabel;
    else
      % there are no neighbours defined for this channel
      complete(i).neighblabel = {};
    end
  end
  neighbours = complete;
end

for i = 1:length(neighbours)
  % convert them into row-arrays for a nicer code representation with PRINTRSTRUCT
  neighbours(i).neighblabel = neighbours(i).neighblabel(:)';
end
  
k = 0;
for i = 1:length(neighbours)
  if isempty(neighbours(i).neighblabel)
    ft_warning('no neighbours found for %s', neighbours(i).label);
  end
  k = k + length(neighbours(i).neighblabel);
end

if k==0
  ft_warning('No neighbouring channels were specified or found');
else
  fprintf('there are on average %.1f neighbours per channel\n', k/length(neighbours));
end

if strcmp(cfg.feedback, 'yes') && ~hasatlas
  % give some graphical feedback
  tmpcfg = keepfields(cfg, {'layout', 'rows', 'columns', 'commentpos', 'skipcomnt', 'scalepos', 'skipscale', 'projection', 'viewpoint', 'rotate', 'width', 'height', 'elec', 'grad', 'opto', 'showcallinfo', 'trackcallinfo', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo', 'checksize'});
  tmpcfg.neighbours = neighbours;
  if hasdata
    tmpcfg.senstype = cfg.senstype;
    ft_neighbourplot(tmpcfg, data);
  else
    ft_neighbourplot(tmpcfg);
  end
elseif strcmp(cfg.feedback, 'yes')
  ft_notice('no visual feedback about the neighbours is possible with the current input');
end

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble previous   data
ft_postamble provenance neighbours
ft_postamble history    neighbours
ft_postamble savevar    neighbours


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION that compute the neighbourhood geometry from the
% gradiometer/electrode positions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [neighbours] = compneighbstructfrompos(chanpos, label, neighbourdist)

nchan = length(label);

% compute the distance between all sensors
dist = zeros(nchan,nchan);
for i = 1:nchan
  dist(i,:) = sqrt(sum((chanpos(1:nchan,:) - repmat(chanpos(i,:), nchan, 1)).^2,2))';
end

% find the neighbouring electrodes based on distance
% later we have to restrict the neighbouring electrodes to those actually selected in the dataset
channeighbstructmat = (dist<neighbourdist);

% electrode istelf is not a neighbour
channeighbstructmat = (channeighbstructmat .* ~eye(nchan));

% convert back to logical
channeighbstructmat = logical(channeighbstructmat);

% construct a struct-array with all neighbours
neighbours=struct;
for i = 1:nchan
  neighbours(i).label       = label{i};
  neighbours(i).neighblabel = label(channeighbstructmat(i,:));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION that computes the neighbourhood geometry from the
% triangulation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [neighbours] = compneighbstructfromtri(chanpos, label, tri)

nchan = length(label);

channeighbstructmat = zeros(nchan,nchan);
% mark neighbours according to triangulation
for i=1:size(tri, 1)
  channeighbstructmat(tri(i, 1), tri(i, 2)) = 1;
  channeighbstructmat(tri(i, 1), tri(i, 3)) = 1;
  channeighbstructmat(tri(i, 2), tri(i, 1)) = 1;
  channeighbstructmat(tri(i, 3), tri(i, 1)) = 1;
  channeighbstructmat(tri(i, 2), tri(i, 3)) = 1;
  channeighbstructmat(tri(i, 3), tri(i, 2)) = 1;
end

% construct a structured cell-array with all neighbours
neighbours = struct;
alldist = [];
for i = 1:nchan
  neighbours(i).label       = label{i};
  neighbidx                 = find(channeighbstructmat(i,:));
  neighbours(i).dist        = sqrt(sum((repmat(chanpos(i, :), numel(neighbidx), 1) - chanpos(neighbidx, :)).^2, 2));
  alldist                   = [alldist; neighbours(i).dist];
  neighbours(i).neighblabel = label(neighbidx);
end

% remove neighbouring channels that are too far away (imporntant e.g. in
% case of missing sensors)
neighbdist = mean(alldist)+3*std(alldist);
for i=1:nchan
  idx = neighbours(i).dist > neighbdist;
  neighbours(i).dist(idx)         = [];
  neighbours(i).neighblabel(idx)  = [];
end
neighbours = rmfield(neighbours, 'dist');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION that ensures a consistent triangulation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [newtri] = tri_reindex(tri)

% this subfunction reindexes tri such that they run from 1:number of unique vertices
newtri       = tri;
[srt, indx]  = sort(tri(:));
tmp          = cumsum(double(diff([0;srt])>0));
newtri(indx) = tmp;