%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% eegchan2src.m
% version 1.0
% 24 August 2017
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This script for MATLAB relies on Brainstorm and OpenMEEG:
% Tadel F, Baillet S, Mosher JC, Pantazis D, Leahy RM (2011)
% Brainstorm: A User-Friendly Application for MEG/EEG Analysis
% Computational Intelligence and Neuroscience, vol. 2011, ID 879716
% http://neuroimage.usc.edu/brainstorm/
% Gramfort A, Papadopoulo T, Olivi E, Clerc M (2010)
% OpenMEEG: opensource software for quasistatic bioelectromagnetics
% Biomedical engineering online, 9(1), 45
% http://openmeeg.github.io
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1    EXPORT FROM BRAINSTORM AND INITIAL PREPARATION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load('cortex.mat') % default anatomy (ICBM152) cortex previously exported
                   % from Brainstorm
load('bem.mat') % OpenMEEG BEM previously exported from BrainStorm

% convert the gain matrix to a constrained source model
% (ref. Brainstorm support thread 918)
gain_constrained = bst_gain_orient(bem.Gain, bem.GridOrient);
clear bem

% remove non-EEG channels (if present)
gain_constrained = gain_constrained([1:18, 21:65], :); % for 10-10

% compute absolute values of each element in gain_constrained and then find the
% maximum value (amplitude) and its index (electrode) per each column (vertex)
[~, index] = max(abs(gain_constrained));

% set the number of electrodes
nChannels = 63; % 10-10 = 63; M10 = 61;

% duplicate the Scouts structure (rows e.g. from cortex.Atlas(2).Scouts)
elec_scouts = cortex.Atlas(2).Scouts(1:nChannels);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2    OBTAIN THE CHANNEL FILE AND FIND MAXIMA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load channel names exported from Brainstorm:
load('elec_1010.mat');

% remove non-EEG channels (if present)
chan = elec_1010.Channel([1:18, 21:65]);

% use the most active vertices for each electrode
% to repopulate the elec_scouts structure and
% set the first vertex as the seed (see Brainstorm documentation)
for i = 1:nChannels
  elec_scouts(i).Vertices = find(index==i);
  elec_scouts(i).Label = chan(i).Name;
  if isempty(elec_scouts(i).Vertices)
    elec_scouts(i).Seed = 0;
  else
    elec_scouts(i).Seed = elec_scouts(i).Vertices(1);
  end
end;

% remove unnecessary variables
clear i index value gain_constrained cortex chan

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3    FIND INTERSECTIONS BETWEEN THE MAXIMA FOR EACH CHANNEL AND CORTICAL AREAS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load a scout file (atlas) from Brainstorm (e.g. Destrieux or Mindboggle) and
% save it to a file (e.g. scout_Mindboggle_62.mat)

% each valid scout file (.mat) must contain the following variables:
% Name (string)
% Scouts (struct [1 x channels])
% TessNbVertices (int number of vertices)

load('scout_Mindboggle_62.mat'); % or any other atlas

% create a structure containing electrode labels and their respective
% intersections between the forward model maxima and areas defined by the atlas

[~, nAreas] = size(Scouts);

chan_areas = struct('channel', [], 'areas', []);
for i = 1:nChannels
  chan_areas(i).channel = elec_scouts(i).Label;
  areas = cell(nAreas, 2);
  for j = 1:nAreas
    areas{j, 1} = Scouts(j).Label;
    areas{j, 2} = length(intersect(elec_scouts(i).Vertices, ...
      Scouts(j).Vertices));
  end
  chan_areas(i).areas = sortrows(areas, -2);
  chan_areas(i).num = sum(cellfun(@double,chan_areas(i).areas(:,2)));
  chan_areas(i).areas(:, 2) = num2cell(cellfun(@(x) ...
    round(x ./ chan_areas(i).num * 100, 1), chan_areas(i).areas(:,2)));
end

% remove unnecessary variables
clear i j areas Name TessNbVertices

% Same as above, only by areas instead of channels
area_chans = struct('area', [], 'channels', []);
for i = 1:nAreas
    area_chans(i).area = Scouts(i).Label;
    channels = cell(nChannels, 2);
    for j = 1:nChannels
        channels{j, 1} = elec_scouts(j).Label;
        channels{j, 2} = length(intersect(Scouts(i).Vertices, ...
          elec_scouts(j).Vertices));
    end
    area_chans(i).channels = sortrows(channels, -2);
    area_chans(i).num = sum(cellfun(@double,area_chans(i).channels(:,2)));
    area_chans(i).channels(:, 2) = num2cell(cellfun(@(x) ...
      round(x ./ area_chans(i).num * 100, 1), area_chans(i).channels(:,2)));
end

% remove unnecessary variables
clear i j areas Name TessNbVertices

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 4    EXPORT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% export the correspondence data to the file output.txt sorted by channel names
fileID = fopen('output-byChannels.csv','w');
threshold = 10; % arbitrary threshold to omit highly insignificant results
for i = 1:nChannels
  fprintf(fileID, '%s', chan_areas(i).channel);
  for j = 1:length(find(cell2mat(chan_areas(i).areas(:,2)) > threshold))
    fprintf(fileID, ',%30s, %3.1f%%\n', chan_areas(i).areas{j, 1}, ...
      chan_areas(i).areas{j, 2});
  end
end
fclose(fileID);

% remove unnecessary variables
clear i j fileID threshold

%%% Same as above, only by areas instead of channels
fileID = fopen('output-byAreas.csv','w');
threshold = 10; % arbitrary threshold to omit highly insignificant results
for i = 1:nAreas
  fprintf(fileID, '%s', area_chans(i).area);
  for j = 1:length(find(cell2mat(area_chans(i).channels(:,2)) > threshold))
    fprintf(fileID, ',%30s, %3.1f%%\n', area_chans(i).channels{j, 1}, ...
      area_chans(i).channels{j, 2});
  end
end
fclose(fileID);

% remove unnecessary variables
clear i j fileID threshold
