mri = ft_read_mri('individual.nii');

ft_determine_coordsys(mri, 'interactive', 'no')
rotate3d

%%

cfg = [];
cfg.method = 'ortho';
cfg.flip = 'no'; % important for identifying voxel indices
ft_sourceplot(cfg, mri_realigned);

nas_vox = [ 102 217 123 ];
ini_vox = [ 90   18 104 ];
lpa_vox = [ 20  135 103 ];
rpa_vox = [ 173 121  95 ];

%%

cfg = [];
cfg.method = 'fiducial';
cfg.coordsys = 'ctf';
cfg.fiducial.nas = nas_vox;
cfg.fiducial.lpa = lpa_vox;
cfg.fiducial.rpa = rpa_vox;
mri_realigned = ft_volumerealign(cfg, mri);

nas = ft_warp_apply(mri.transform, nas_vox);
ini = ft_warp_apply(mri.transform, ini_vox);
lpa = ft_warp_apply(mri.transform, lpa_vox);
rpa = ft_warp_apply(mri.transform, rpa_vox);

%%

cfg = [];
cfg.method = 'ortho';
cfg.location = nas; % nas, ini, lpa, rpa
cfg.locationcoordinates = 'head';
ft_sourceplot(cfg, mri_realigned);

%%

cfg = [];
mri_resliced = ft_volumereslice(cfg, mri_realigned);

cfg = [];
cfg.method = 'ortho';
cfg.location = [0 0 0];
cfg.locationcoordinates = 'head';
ft_sourceplot(cfg, mri_resliced);

disp(mri_resliced.cfg)

cfg = [];
cfg.xrange = [ -97.5000 157.5000] - 30;
cfg.yrange = [-127.5000 127.5000];
cfg.zrange = [ -87.5000 167.5000] - 5;
mri_resliced = ft_volumereslice(cfg, mri_realigned);

cfg = [];
cfg.method = 'ortho';
cfg.location = [0 0 0];
cfg.locationcoordinates = 'head';
ft_sourceplot(cfg, mri_resliced);

%%

cfg = [];
cfg.method = 'ortho';
cfg.location = nas; % nas, ini, lpa, rpa
cfg.locationcoordinates = 'head';
ft_sourceplot(cfg, mri_resliced);

%%

ft_determine_coordsys(mri_resliced, 'interactive', 'no')
rotate3d

%%

cfg = [];
cfg.output = 'scalp';
mri_segmented = ft_volumesegment(cfg, mri_resliced);

%%

cfg = [];
cfg.method = 'projectmesh';
cfg.numvertices = 4000;
headshape = ft_prepare_mesh(cfg, mri_segmented);

%%

headshape.fid.pos = [
  nas
  lpa
  rpa
  ini
  ];

headshape.fid.label = {
  'nas'
  'lpa'
  'rpa'
  'ini'
  };

figure
ft_plot_headshape(headshape, 'facecolor', 'skin', 'facealpha', 0.7, 'fidmarker', 'o', 'fidcolor', 'k', 'fidlabel', true, 'fidsize', 24)
ft_headlight

%%

cfg = [];
cfg.funparameter = 'scalp';
ft_sourceplot(cfg, mri_segmented)

%%

mri_segmented.scalp(:,:,1:50) = 0;
mri_segmented.airgap = imdilate(mri_segmented.scalp,  strel('sphere', 3));
mri_segmented.helmet = imdilate(mri_segmented.airgap, strel('sphere', 7));

mri_indexed = ft_checkdata(mri_segmented, 'segmentationstyle', 'indexed');

cfg = [];
cfg.funparameter = 'tissue';
cfg.location = [0 0 0];
cfg.locationcoordinates = 'head';
cfg.atlas = mri_indexed;
ft_sourceplot(cfg, mri_indexed)

%%

cfg = [];
cfg.method = 'isosurface';
cfg.numvertices = inf;

cfg.tissue = 'airgap'; % this makes a surface from the outside of the "airgap" part
inside = ft_prepare_mesh(cfg, mri_segmented);

cfg.tissue = 'helmet';
outside = ft_prepare_mesh(cfg, mri_segmented);

%%

figure
ft_plot_headshape(headshape, 'facecolor', 'skin', 'facealpha', 0.5, 'edgecolor', 'none');
ft_plot_mesh(inside, 'facecolor', 'r', 'facealpha', 0.5, 'edgecolor', 'none');
ft_headlight
rotate3d

figure
ft_plot_headshape(headshape, 'facecolor', 'skin', 'facealpha', 0.5, 'edgecolor', 'none');
ft_plot_mesh(outside, 'facecolor', 'g', 'facealpha', 0.5, 'edgecolor', 'none');
ft_headlight
rotate3d

%%

cfg = [];
cfg.fiducial.nas = nas;
cfg.fiducial.ini = ini;
cfg.fiducial.lpa = lpa;
cfg.fiducial.rpa = rpa;
cfg.method = '1020';
elec = ft_electrodeplacement(cfg, headshape);

%%

chansel = ft_channelselection({'eeg1020', '-Fpz', '-Oz'}, elec.label);

cfg = [];
cfg.elec = elec;
cfg.channel = chansel;
cfg.outwardshift = 10;

cfg.template = 'fieldline_sensor.stl';
sensor = ft_sensorplacement(cfg, headshape);
cfg.template = 'fieldline_holder.stl';
holder = ft_sensorplacement(cfg, headshape);
cfg.template = 'fieldline_cuttingtool.stl';
cuttingtool = ft_sensorplacement(cfg, headshape);

%%

figure
ft_plot_headshape(headshape, 'facecolor', 'skin', 'facealpha', 1, 'axes', 1);
ft_plot_mesh(outside, 'facecolor', 'lightgray', 'facealpha', 0.5);

for i=1:numel(sensor)
  ft_plot_headshape(sensor(i), 'facecolor', 'r', 'facealpha', 1);
  ft_plot_headshape(holder(i), 'facecolor', 'g', 'facealpha', 1);
  ft_plot_headshape(cuttingtool(i), 'facecolor', 'b', 'facealpha', 0.5);
end
ft_headlight

%%

mkdir individual

ft_write_headshape('individual/helmet-inside.stl', inside, 'fileformat', 'stl');
ft_write_headshape('individual/helmet-outside.stl', outside, 'fileformat', 'stl');

for i=1:numel(sensor)
  disp(chansel{i});
  filename = sprintf('individual/sensor-%s.stl', chansel{i});
  ft_write_headshape(filename, sensor(i), 'fileformat', 'stl');
  filename = sprintf('individual/holder-%s.stl', chansel{i});
  ft_write_headshape(filename, sensor(i), 'fileformat', 'stl');
  filename = sprintf('individual/cuttingtool-%s.stl', chansel{i});
  ft_write_headshape(filename, sensor(i), 'fileformat', 'stl');
end
