clear
%Reading in the segmentation information
data = loadnifti('datasets/2a/c2ANTS45-49Years_head.nii.gz');
seg.wm = data.NIFTIData;
data = loadnifti('datasets/2a/c1ANTS45-49Years_head.nii.gz');
seg.gm = data.NIFTIData;
data = loadnifti('datasets/2a/c3ANTS45-49Years_head.nii.gz');
seg.csf = data.NIFTIData;
data = loadnifti('datasets/2a/c4ANTS45-49Years_head.nii.gz');
seg.skull = data.NIFTIData;
data = loadnifti('datasets/2a/c5ANTS45-49Years_head.nii.gz');
seg.scalp = data.NIFTIData;

%Using coarser parameters than the default values
cfg.samplingwm = 2;
cfg.samplinggm = 2;
cfg.samplingcsf = 2.5;
cfg.samplingskull = 3.5;
cfg.samplingscalp = 3.5;
cfg.meshfield = 5;
cfg.maxvol = 50;

%Generating the mesh from the seg and cfg
[node,elem,face] = brain2mesh(seg,cfg);

% plotmesh(node,face)

%Inputs to the 10-10 tool: Reference Points in order nasion Nz, inion Iz and 
%two preauricular points, AR on the right side and AL on the left side.
sampleinputs = [83.77 191 98.05;
90.18 9.13 87.14;
156 104.1 73.86;
25.03 99.26 70.42];

face_scalp = volface(elem(:,1:4));
[node_scalp,face_scalp] = meshcheckrepair(node(:,1:3),face_scalp(:,1:3),'isolated');
cutoff = min(node_scalp(:,3));
node_scalp = node_scalp(node_scalp(:,3)>cutoff+1,:);

[P, N] = nfri_auto1010(sampleinputs, node_scalp(:,1:3), 'hemispheric', 'balanced', [], 1, 'PrintPlot');


%Let's say we use the F1 position - point 31 in a 10-10
POI = 31;

cfg2.nphoton=1e5;
cfg2.node = node(:,1:3);
cfg2.elem = elem(:,1:4);
cfg2.elemprop=elem(:,5);

% [a,b]=mmcraytrace(cfg2.node,cfg2.elem,cfg2.srcpos,cfg2.srcdir,cfg2.e0);
% cfg.srcpos=a;
cfg2.srcpos=P(POI,1:3)+rand*0.01;
cfg2.srcdir=[0 0 -1];
cfg2.prop=[0 0 1 1; 0.07 40.1 0.85 1.37; 0.02 8.4 0.9 1.37; 0.0004 1 0.99 1.37; 0.0101 100 0.99 1.37; 0.0101 80 0.99 1.37; 0 0 1 1];
cfg2.tstart=0;
cfg2.tend=5e-9;
cfg2.tstep=5e-9;
cfg2.debuglevel='TP';
cfg2.e0 = '-';

% T = tsearchn(node(:,1:3),elem(:,1:4),a)

flux=mmclab(cfg2);

qmeshcut(cfg2.elem(:,1:4),cfg2.node(:,1:3),log10(abs(flux.data(:))),sprintf('y=%i',P(POI,2)),'linestyle','none');
hold on;
[cutpos,cutvalue,facedata]=qmeshcut(elem(:,1:4),node(:,1:3),log10(abs(flux.data(:))),102); 

[bcutpos,bcutvalue,bcutedges]=qmeshcut(face(:,1:3),node(:,1:3),node(:,1),sprintf('y=%i',P(POI,2)));
[bcutpos,bcutedges]=removedupnodes(bcutpos,bcutedges);
bcutloop=extractloops(bcutedges);

% plot the plane-surface cuts
delimiter = isnan(bcutloop);
k=find(delimiter);
k = [0 k];
for i = 1:length(k)-1
    plot3(bcutpos(bcutloop(k(i)+1:k(i+1)-1),1),bcutpos(bcutloop(k(i)+1:k(i+1)-1),2),bcutpos(bcutloop(k(i)+1:k(i+1)-1),3),'r','LineWidth',1);
end
axis off
view([0 0])