mri = ft_read_mri('002-t1_mprage_sag_p2_iso_1.0_20ch_head/PAUGAA_20160825_MEGANATOMY.MR.DCCN_PRISMAFIT.0002.0184.2016.08.25.16.54.09.114696.482176288.IMA');

%%

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

%%

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

cfg = [];
cfg.location = ini_vox;
cfg.locationcoordinates = 'voxel';
cfg.flip = 'no';
ft_sourceplot(cfg, mri_flipped);

%%

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

%%

nas_orig = ft_warp_apply(mri_realigned.transform, nas_vox);
lpa_orig = ft_warp_apply(mri_realigned.transform, lpa_vox);
rpa_orig = ft_warp_apply(mri_realigned.transform, rpa_vox);
ini_orig = ft_warp_apply(mri_realigned.transform, ini_vox);

%%

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.location = nas_orig;
cfg.locationcoordinates = 'head';
cfg.flip = 'no';
ft_sourceplot(cfg, mri_resliced);

%%

ft_write_mri('individual.nii', mri_realigned.anatomy, 'transform', mri_realigned.transform, 'dataformat', 'nifti');

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

mri_realigned = ft_read_mri('individual.nii');

nsubj = 10clos;

mkdir population

for i=1:nsubj
  fprintf('------------------------------- %d -------------------------------\n', i);

  % construct transformations with an X amount of scaling/rotation/translation
  % implemented as a uniform distribution of -X/2 up to +X/2

  s = ones(1,3)  + (rand(1,3)-1.0) * 0.1; % not bigger, but up to 10% smaller
  r = zeros(1,3) + (rand(1,3)-0.5) * 10; % 10 degree
  t = zeros(1,3) + (rand(1,3)-0.5) * 10; % 10 mm

  Ts = ft_warp_apply(s, eye(3), 'scale'); Ts(4,4) = 1;
  Tr = ft_warp_apply(r, eye(3), 'rotate'); Tr(4,4) = 1;
  Tt = eye(4); Tt(1:3,4) = ft_warp_apply(t, zeros(1,3), 'translate');

  % apply them in a random order
  T = {Ts, Tr, Tt};
  T = T(randperm(3));
  T = T{1} * T{2} * T{3};

  mri_population = ft_transform_geometry(T, mri_realigned);

  nas = ft_warp_apply(T, nas_orig);
  ini = ft_warp_apply(T, ini_orig);
  lpa = ft_warp_apply(T, lpa_orig);
  rpa = ft_warp_apply(T, rpa_orig);

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

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

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

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

  cfg = [];
  cfg.method = 'ortho';
  cfg.location = nas;
  cfg.locationcoordinates = 'head';
  % ft_sourceplot(cfg, mri_population);

  mri_population.anatomy = int16(16384 * mri_population.anatomy/max(mri_population.anatomy(:)));

  filename = sprintf('population/subject%03d.nii', i);
  ft_write_mri(filename, mri_population.anatomy, 'transform', mri_population.transform, 'dataformat', 'nifti');

  filename = sprintf('population/fiducial%03d.mat', i);
  save(filename, 'nas', 'ini', 'lpa', 'rpa');

end

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nsubj = 10;

mri = cell(1, nsubj);
fiducial = cell(1, nsubj);

for i=1:nsubj
  fprintf('------------------------------- %d -------------------------------\n', i);
  mrifile = sprintf('population/subject%03d.nii', i);
  fidfile = sprintf('population/fiducial%03d.mat', i);
  mri{i} = ft_read_mri(mrifile);
  mri{i}.coordsys = 'ctf';
  fiducial{i} = load(fidfile);
end

%%

ft_determine_coordsys(mri{1}, 'interactive', 'no')

%%

cfg = [];
cfg.method = 'ortho';
cfg.location = fiducial{1}.nas;
cfg.locationcoordinates = 'head';
ft_sourceplot(cfg, mri{1});

%%

mri_resliced = cell(1, nsubj);
mri_segmented = cell(1, nsubj);

for i=1:nsubj
  fprintf('------------------------------- %d -------------------------------\n', i);
  cfg = [];
  cfg = [];
  cfg.xrange = [ -97.5000 157.5000] - 30;
  cfg.yrange = [-127.5000 127.5000];
  cfg.zrange = [ -87.5000 167.5000] - 5;
  mri_resliced{i} = ft_volumereslice(cfg, mri{i});

  cfg = [];
  cfg.output = 'scalp';
  mri_segmented{i} = ft_volumesegment(cfg, mri_resliced{i});
end

%%

cfg = [];
cfg.method = 'ortho';
cfg.location = fiducial{1}.nas;
cfg.funparameter = 'scalp';
ft_sourceplot(cfg, mri_segmented{1});

%%

mri_averaged = rmfield(mri_segmented{1}, 'cfg');
for i=2:nsubj
  mri_averaged.scalp = mri_averaged.scalp + mri_segmented{i}.scalp;
end
mri_averaged.scalp = mri_averaged.scalp/nsubj;

cfg = [];
cfg.method = 'ortho';
cfg.location = [0 0 0];
cfg.funparameter = 'scalp';
ft_sourceplot(cfg, mri_averaged);

%%

mri_90percentile       = mri_averaged;
mri_90percentile.scalp = mri_averaged.scalp>0.1;

% make it perfectly symmetric
mri_90percentile.scalp =  mri_90percentile.scalp | flip(mri_90percentile.scalp, 2);

cfg = [];
cfg.method = 'ortho';
cfg.location = [0 0 0];
cfg.funparameter = 'scalp';
ft_sourceplot(cfg, mri_90percentile);

%%

mri_segmented = mri_90percentile;

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.method = 'isosurface';
cfg.numvertices = inf;

cfg.tissue = 'scalp';
headshape = ft_prepare_mesh(cfg, mri_segmented);

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);

%%

nas_avg = fiducial{1}.nas;
ini_avg = fiducial{1}.ini;
lpa_avg = fiducial{1}.lpa;
rpa_avg = fiducial{1}.rpa;

for i=2:nsubj
  nas_avg = nas_avg + fiducial{i}.nas;
  ini_avg = ini_avg + fiducial{i}.ini;
  lpa_avg = lpa_avg + fiducial{i}.lpa;
  rpa_avg = rpa_avg + fiducial{i}.rpa;
end

nas_avg = nas_avg/nsubj;
ini_avg = ini_avg/nsubj;
lpa_avg = lpa_avg/nsubj;
rpa_avg = rpa_avg/nsubj;

% don't immediately overwrite the variables
tmp1 = (lpa_avg + [1 -1 1] .* rpa_avg)/2;
tmp2 = (rpa_avg + [1 -1 1] .* lpa_avg)/2;
% now we can overwrite them
lpa_avg = tmp1;
rpa_avg = tmp2;

%%

headshape.fid.pos = [
  nas_avg
  ini_avg
  lpa_avg
  rpa_avg
  ];

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

ft_plot_headshape(headshape, 'facecolor', 'skin', 'facealpha', 0.5, 'fidmarker', '.', 'fidcolor', 'k', 'fidlabel', true, 'fidsize', 24)
ft_headlight

%%

cfg = [];
cfg.fiducial.nas = nas_avg;
cfg.fiducial.ini = ini_avg;
cfg.fiducial.lpa = lpa_avg;
cfg.fiducial.rpa = rpa_avg;
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);

%%

mkdir population

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

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