%% import MRI data

% setup directories and files:
D = load([fn(1:end-11) '_ICAfffd.mat']);
D = D.D;
D_ft = load(fn_ft);
D_ft = D_ft.D;
mrifile = D.other.inv{1}.mesh.sMRI(1:end-2);

% load basic data:
mri_orig = ft_read_mri(mrifile);
grad     = ft_read_sens(fn,'senstype','meg');
elec     = D_ft.elec; % ft_read_sens(dataset,'senstype','eeg');
shape    = ft_read_headshape(fn,'unit','mm');

% check it fits okay:
figure, subplot(121)
ft_plot_headshape(shape);
ft_plot_sens(grad, 'style', '*b');
ft_plot_sens(elec, 'style', '*g');
view([1 0 0])


%% neuromag coord-sys

% Do the ICP
mri_orig.coordsys = 'acpc';
cfg           = [];
cfg.method    = 'headshape';
cfg.headshape.interactive = 'no';
cfg.headshape.icp = 'yes';
cfg.headshape.headshape = shape; 
cfg.coordsys  = 'neuromag';
cfg.spmversion = 'spm12';
cfg.viewresult    = 'no';
mri = ft_volumerealign(cfg, mri_orig); % note - to avoid a popup/interactive window, i had to set line 52 of ft_determine_coordsys to interactive = no
mri.coordsys = 'neuromag';

% reslice at 1mm for niceness:
cfg = [];
cfg.resolution = 1;
cfg.xrange = [-100 100];
cfg.yrange = [-110 140];
cfg.zrange = [-80 180];
mri = ft_volumereslice(cfg, mri);
mri = ft_convert_units(mri, 'cm');

% save:
save([fileparts(fn) '/FT_MRI_realign_reslice.mat'],'mrifile','mri','shape')


%% segmenting and meshes

rmpath('/imaging/local/software/spm_cbu_svn/releases/spm8_latest')
addpath('/imaging/local/software/spm_cbu_svn/releases/spm12_latest')
cfg           = [];
cfg.output    = {'brain', 'skull', 'scalp'};
mri_segmented = ft_volumesegment(cfg, mri);

% copy the anatomy into the segmented mri
mri_segmented.anatomy = mri.anatomy;

% check it:
cfg = [];
cfg.funparameter = 'brain';
ft_sourceplot(cfg, mri_segmented);
cfg.funparameter = 'skull';
ft_sourceplot(cfg, mri_segmented);
cfg.funparameter = 'scalp';
ft_sourceplot(cfg, mri_segmented);

% save it:
save([fileparts(fn) '/FT_MRI_segm_mesh.mat'],'mri_segmented')


%% meshes:

rmpath('/imaging/local/software/spm_cbu_svn/releases/spm12_latest')
addpath('/imaging/local/software/spm_cbu_svn/releases/spm8_latest')

cfg = [];
cfg.method = 'projectmesh';
cfg.tissue = 'brain';
cfg.numvertices = 3000;
mesh_brain = ft_prepare_mesh(cfg, mri_segmented);

cfg = [];
cfg.method = 'projectmesh';
cfg.tissue = 'skull';
cfg.numvertices = 2000;
mesh_skull = ft_prepare_mesh(cfg, mri_segmented);

cfg = [];
cfg.method = 'projectmesh';
cfg.tissue = 'scalp';
cfg.numvertices = 1000;
mesh_scalp = ft_prepare_mesh(cfg, mri_segmented);

cfg = [];
cfg.method = 'isosurface';
cfg.tissue = 'scalp';
cfg.numvertices = 16000;
highres_scalp = ft_prepare_mesh(cfg, mri_segmented);

% check it:
figure
ft_plot_mesh(highres_scalp, 'edgecolor', 'none', 'facecolor', 'skin')
material dull, camlight, lighting phong

% save it:
save([fileparts(fn) '/FT_MRI_segm_mesh.mat'],'mesh_brain','mesh_skull','mesh_scalp','-append')


%% headmodel

cfg = [];
cfg.method = 'singleshell';
headmodel_meg = ft_prepare_headmodel(cfg, mesh_brain);
headmodel_meg = ft_convert_units(headmodel_meg,'cm');

% check it:
figure
hold on
ft_plot_headshape(shape);
ft_plot_sens(grad, 'style', 'ob');
ft_plot_sens(elec, 'style', 'og');
ft_plot_headmodel(headmodel_meg, 'facealpha', 0.5, 'edgecolor', 'none'); % "lighting phong" does not work with opacity
material dull
camlight
view([1 0 0])

% save it:
save([fileparts(fn) '/FT_MRI_headmod.mat'],'headmodel_meg','shape')


%% getting atlas coords inside you headmodel space:

% load your atlas
atlas = ft_read_atlas('/home/na01/MATLAB TA/mfiles/downloaded on path w subs/ROI_MNI_V4.nii'); % HCPMMP1_on_MNI152_ICBM2009a_nlin_hd_TA.nii'); % ROI_MNI_V4.nii');
atlas_shortname = 'AAL'; % 'HCPMMP1';
atlas = ft_convert_units(atlas,'cm');

% % make a template grid:
% vol = load('/home/na01/MATLAB TA/mfiles/downloaded on path w subs/standard_singleshell.mat');
% vol = vol.vol;
% cfg = [];
% cfg.xgrid  = -20:.5:20;
% cfg.ygrid  = -20:.5:20;
% cfg.zgrid  = -20:.5:20;
% cfg.unit   = 'cm';
% cfg.tight  = 'yes';
% cfg.inwardshift = -.5;
% cfg.headmodel   = vol;
% template_grid   = ft_prepare_sourcemodel(cfg);
% template_grid = ft_convert_units(template_grid,'cm');

% !! or use downloaded template: !!
template_filename = '/home/na01/Downloads/standard_sourcemodel3d4mm.mat';
load(template_filename)
tmp = whos('-file',template_filename);
template_grid = eval(tmp.name);
template_grid = ft_convert_units(template_grid, 'cm');

% finds which template coords are inside ALL atlas rois:
cfg = [];
cfg.atlas      = atlas;
cfg.inputcoord = 'mni';
cfg.roi        = atlas.tissuelabel;  % here you can also specify a single label, i.e. single ROI
mask_all       = ft_volumelookup(cfg, template_grid);

% remove vertices not in the template:
template_grid_all = template_grid; % template_grid; % 
template_grid_all.inside = false(template_grid.dim);
template_grid_all.inside(mask_all==1) = true;

% inverse warp the subject-specific grid to the template grid:
cfg             = [];
cfg.spmversion  = 'spm12';
cfg.warpmni     = 'yes';
cfg.template    = template_grid;
cfg.nonlinear   = 'yes'; % use non-linear normalization
cfg.mri         =  mri;
sourcemodel_all = ft_prepare_sourcemodel(cfg);
sourcemodel_all = ft_convert_units(sourcemodel_all, 'cm');
sourcemodel_all.pos(~sourcemodel_all.inside,:) = [];
sourcemodel_all = rmfield(sourcemodel_all,'inside');

% check it:
figl; subplot(131), hold on
ft_plot_mesh(template_grid.pos(template_grid.inside,:));
ft_plot_headmodel(vol,  'facecolor', 'cortex', 'edgecolor', 'none');
alpha 0.5, camlight
subplot(132), hold on
ft_plot_mesh(template_grid_all.pos(template_grid_all.inside,:));
ft_plot_headmodel(vol,  'facecolor', 'cortex', 'edgecolor', 'none');
alpha 0.5, camlight
subplot(133), hold on
ft_plot_headmodel(headmodel_meg,  'facecolor', 'cortex', 'edgecolor', 'none');
ft_plot_axes(headmodel_meg); alpha 0.4
ft_plot_mesh(sourcemodel_all.pos);


%%