% SIMBIO SOURCE RECONSTRUCTION EEG FEM

%%%

% The image preprocessing steps are skipped in this tutorial, the image can
% be directly loaded from the file mri.mat

%%%

% HEADMODEL
% read mri
% mri = ft_read_mri();
load mri.mat

disp(mri)
%               dim: [256 256 192]
%           anatomy: [256x256x192 double]
%               hdr: [1x192 struct]
%         transform: [4x4 double]
%              unit: 'mm'
%     transformorig: [4x4 double]
%          coordsys: 'ctf'

% realign mri

% Using the n/l/r keys the fiducials can be specified, the
%   landmarks can be specified with a/p/z. When pressing q the interactive
%   mode will stop and the transformation matrix is computed.

cfg = [];
cfg.method = 'interactive';

% alternatively pre-defined fiducials (from the paper) can be used:

% cfg = [];
% cfg.method = 'fiducial';
% cfg.fiducial.nas = [36,93,92];
% cfg.fiducial.lpa = [143,142,184];
% cfg.fiducial.rpa = [150,148,10];
% cfg.fiducial.zpoint = [140,1,96];

mri = ft_volumerealign(cfg, mri);

disp(mri)
%               dim: [256 256 192]
%           anatomy: [256x256x192 double]
%               hdr: [1x192 struct]
%         transform: [4x4 double]
%              unit: 'mm'
%               cfg: [1x1 struct]
%     transformorig: [4x4 double]
%          coordsys: 'ctf'

% reslice and plot
%
% We want to keep as much of the MRI here as possible (especially also in
% the neck area, therefore we only apply cfg.method = 'flip').
% For visualization purposes, one can later apply ft_volumereslice
% regularly

cfg = [];
cfg.method = 'flip';
mri = ft_volumereslice(cfg,mri);

%save mri_resliced.mat mri

%%%

% segmentation
cfg = [];
cfg.output = {'gray','white','csf','skull','scalp'};
cfg.spmversion = 'spm12';
cfg.scalpthreshold = 0.06;
segmentedmri = ft_volumesegment(cfg, mri);

%save segmentedmri segmentedmri

disp(segmentedmri)
%           dim: [256 256 192]
%     transform: [4x4 double]
%      coordsys: 'ctf'
%          unit: 'mm'
%          gray: [256x256x192 logical]
%         white: [256x256x192 logical]
%           csf: [256x256x192 logical]
%         skull: [256x256x192 logical]
%         scalp: [256x256x192 logical]
%           cfg: [1x1 struct]

% plot segmentation image
seg_i            = ft_datatype_segmentation(segmentedmri,'segmentationstyle',...
                                            'indexed');
cfg              = [];
cfg.funparameter = 'seg';
cfg.location     = 'center';

%ft_sourceplot(cfg,seg_i);

t(3) = toc

% create mesh
cfg=[];
cfg.shift = 0.3;
cfg.method = 'hexahedral';
mesh = ft_prepare_mesh(cfg,segmentedmri);

%save mesh.mat mesh

disp(mesh)
%             hex: [5196230x8 double]
%             pos: [5317016x3 double]
%          tissue: [5196230x1 double]
%     tissuelabel: {'csf'  'gray'  'scalp'  'skull'  'white'}
%            unit: 'mm'
%        coordsys: 'ctf'
%             cfg: [1x1 struct]

% create headmodel
cfg = [];
cfg.method = 'simbio';
cfg.conductivity = [0.33 0.14 1.79 0.01 0.43]; % same as tissuelabel in vol
vol = ft_prepare_headmodel(cfg, mesh);

%save vol.mat vol

disp(vol)
%             pos: [5317016x3 double]
%             hex: [5196230x8 double]
%          tissue: [5196230x1 double]
%            cond: [0.3300 0.1400 1.7900 0.0100 0.4300]
%     tissuelabel: {'csf'  'gray'  'scalp'  'skull'  'white'}
%           stiff: [5317016x5317016 double]
%            type: 'simbio'
%            unit: 'mm'
%             cfg: [1x1 struct]

%%%

% The EEG data preprocessing steps are skipped in this tutorial, the image can
% be directly loaded from the file mri.mat

%%%

% SENSORS
% read eeg-information from the header of the combined meg/eeg-file
%hdr = ft_read_header('', ...
%                     'headerformat', 'ctf_ds');
%elec.elecpos = [hdr.orig.eeg.pos]';
%elec.chanpos = [hdr.orig.eeg.pos]';
%elec.label = {hdr.orig.eeg.name}';
%elec.type = 'eeg';
%elec.unit = 'cm';

% convert units
%elec = ft_convert_units(elec,'mm');

% plot sensors and headmodel
%figure;
%ft_plot_sens(elec,'style','r');
%hold on;
%ft_plot_mesh(mesh,'surfaceonly','yes','vertexcolor','none','facecolor',...
%             'skin','facealpha',0.5,'edgealpha',0.1)
%ft_plot_sens(elec,'style', 'k');    
%hold off;

% align sensors
%cfg          = [];
%cfg.method   = 'interactive';
%cfg.elec     = elec;
%cfg.headshape = vol;
%elec_aligned = ft_electroderealign(cfg);
%cfg          = [];
%cfg.method   = 'project';
%cfg.elec     = elec_aligned;
%cfg.headshape = vol;
%elec_projected = ft_electroderealign(cfg);
%save elec_projected.mat elec_projected

load('elec_projected.mat');

% FUNCTIONAL DATA
% trial definition
%cfg                         = [];
%cfg.dataset                 = ''; % enter your dataset name here
%cfg.channel                 = ''; % enter your channel layout here
%cfg.trialfun                = 'ft_trialfun_general'; % this is the default
%cfg.trialdef.eventtype      = 'left';
%cfg.trialdef.prestim        = 0.15; % in seconds
%cfg.trialdef.poststim       = 0.15; % in seconds
%cfg = ft_definetrial(cfg);
% found 5835 events
% created 967 trials

% preprocessing
%cfg.channel    = {'EEG'};
%cfg.continuous = 'no';
%dataLeft = ft_preprocessing(cfg);
%disp(dataLeft)
%            hdr: [1x1 struct]
%          label: {81x1 cell}
%           time: {1x967 cell}
%          trial: {1x967 cell}
%        fsample: 1200
%     sampleinfo: [967x2 double]
%           grad: [1x1 struct]
%            cfg: [1x1 struct]

% filtering
%cfg = [];
%cfg.hpfilter = 'yes';
%cfg.lpfilter = 'yes';
%cfg.dftfilter = 'yes';
%cfg.hpfreq = 20;
%cfg.lpfreq = 250;
%cfg.dftfreq = [50 100 150];
%cfg.demean     = 'yes';
%cfg.baselinewindow  = [-0.2 0];
%dataLeftFiltered = ft_preprocessing(cfg,dataLeft);

%load('dataLeftFiltered.mat');

% example signal at the 30th electrode in the 17th trial
%plot(dataLeftFiltered.time{17}, dataLeftFiltered.trial{17}(30,:))
%save dataLeftFiltered dataLeftFiltered

% artifact rejection

% delete some channels with 'summary' and 'trial' option
%cfg          = [];
%cfg.method   = 'summary';
%cfg.alim     = 1e-12; 
%dummy        = ft_rejectvisual(cfg,dataLeftFiltered); 

%cfg          = [];
%cfg.method   = 'trial';
%dummy        = ft_rejectvisual(cfg,dummy);

% delete some trials with databrowser option
%cfg=[];
%artf=ft_databrowser(cfg,dummy);

% do timelock analysis
%cfg = [];
%cfg.covariance = 'yes';
%cfg.covariancewindow = [-inf 0];
%tlaLeft = ft_timelockanalysis(cfg, dummy);
%save tlaLeft tlaLeft;

load('tlaLeft');

% plot timelocked analysis averaged solution
%plot(tlaLeft.time,tlaLeft.avg);

% FORWARD SOLUTION

% precompute transfer matrix (computaion for 81 electrodes takes about 6
% hours and requiers a memory of 2.2 GB)

% compute transfer matrix
vol_transfer = ft_prepare_vol_sens(vol,elec_projected);

%save vol_transfer.mat vol_transfer -v7.3

% compute leadfield / forward solution
cfg = [];

% original grid resolution
%cfg.grid.resolution = 2.0;

% reduced grid resolution for example
cfg.grid.resolution = 7.5;

cfg.mri = segmentedmri;
cfg.threshold = 0.1;
cfg.smooth = 5;
cfg.elec = elec_projected;
cfg.vol = vol_transfer;
[grid, cfg] = ft_prepare_sourcemodel(cfg, vol_transfer, elec_projected);
% 3514 dipoles inside, 5913 dipoles outside brain

%leadfield:
cfg = [];
cfg.grid = grid;
cfg.headmodel = vol_transfer;
cfg.elec = elec_projected;
cfg.reducerank = 3;
[grid, cfg] = ft_prepare_leadfield(cfg, tlaLeft);

% save grid with leadfield:
%save grid.mat grid;

% INVERSE SOLUTION
% dipole fit
cfg = [];
cfg.numdipoles = 1;
cfg.model = 'regional';
cfg.latency = [0.0235,0.0245];
cfg.grid = grid;
cfg.headmodel = vol_transfer;
cfg.elec = elec_projected;
cfg.nonlinear = 'no';
source = ft_dipolefitting(cfg, tlaLeft);

% VISUALIZATION

% figure;
% hold on
% 
% ft_plot_dipole(source.dip.pos(1,:),mean(source.dip.mom(1:3,:),2),'color','b','unit','mm');
% ft_plot_slice(mri.anatomy,'transform',mri.transform,'location',pos,'orientation',[1 0 0],'resolution',0.1)
% ft_plot_slice(mri.anatomy,'transform',mri.transform,'location',pos,'orientation',[0 1 0],'resolution',0.1)
% ft_plot_slice(mri.anatomy,'transform',mri.transform,'location',pos,'orientation',[0 0 1],'resolution',0.1)
% 
% ft_plot_crosshair(pos, 'color', [1 1 1]/2);
% 
% axis tight
% axis off
% 
% view(120,20)

%% coronal slice
%view(90,0)

%% sagittal slice
%view(180,0)

%% axial slice
%view(0,90)

% ________________________________________________________________
