clear cfg;
close all
clear all


% Comment the below line if you have an NVIDIA GPU
USE_MCXCL=1;

% Load data (0.1 mm^3/voxel)
load digimouse;
figure,imagesc(squeeze(digimouse(:,:,50)))

% Simulated photon number: 1e6 photon
cfg.nphoton=3e6;

% Downsample by 4 (0.4 mm^3/voxel)
volume = digimouse(1:4:end,1:4:end,1:4:end);
volume = volume(:,90:190,:);
cfg.vol=volume;

% Tissue properties
cfg.prop=[0 0 0 1 1
1 0.0191 6.6 0.9 1.37
2 0.0136 8.6 0.9 1.37
3 0.0026 0.01 0.9 1.37
4 0.0186 11.1 0.9 1.37
5 0.0186 11.1 0.9 1.37
6 0.0186 11.1 0.9 1.37
7 0.0186 11.1 0.9 1.37
8 0.0186 11.1 0.9 1.37
9 0.0240 8.9  0.9 1.37
10 0.0026 0.01 0.9 1.37
11 0.0240 8.9  0.9 1.37
12 0.0240 8.9  0.9 1.37
13 0.0240 8.9  0.9 1.37
14 0.0240 8.9  0.9 1.37
15 0.0240 8.9  0.9 1.37
16 0.072  5.6 0.9 1.37
17 0.072  5.6 0.9 1.37
18 0.072  5.6 0.9 1.37
19 0.050  5.4 0.9 1.37
20 0.024  8.9 0.9 1.37
21 0.076 10.9 0.9 1.37];
cfg.prop=cfg.prop(:,2:5);

% Define source position: (x,y,z)=(26 38 52) in voxel unit
cfg.srcpos=[26 38 52];
% Define source direction: unitary vector
cfg.srcdir=[0 0 -1];

% info=mcxlab('gpuinfo')
% GPU id
cfg.gpuid=1;
% Automatically set threads and blocks
cfg.autopilot=1;

% Set starting time
cfg.tstart=0;
% TODO: Set simulation end time (4 ns)
cfg.tend=4e-9;
% TODO: Time step (0.4 ns)
cfg.tstep=4e-10;

%cfg.printnum=10;
cfg.seed=99999;

% Tell the unit in mm
% Original data (0.1mm/voxel) is downsampled by 4, so currently it is
% 0.4mm/voxel
cfg.unitinmm=0.1*4;
% Define source type (pencil)
cfg.srctype='pencil';

% Define detector
cfg.issaveexit = 1;
cfg.detpos = [28 40 50 2];

% Run MCXLAB
[flux1,detphoton]=mcxlab(cfg);

% TODO: compute the CW solution for flux1
CW = sum(flux1.data,4);

% Planar source
cfg.srctype='planar';
cfg.srcparam1=[10 0 0 0];
cfg.srcparam2=[0 10 0 0];
[flux2,detphoton]=mcxlab(cfg);
CW2 = sum(flux2.data,4);

%% What you can do
% Change photon number
% Change volume
% Change source type
% Plot TPSF

%% Visualization
vol = single(cfg.vol);
vol = log(vol)/log(4)+1;

V = vol;
G1 = squeeze(flux2.data(:,:,:,1));
mcx = log10(G1);
mcx(mcx<0) = 0;
idx = find(mcx>0);
V(idx) = mcx(idx);
slice = squeeze(V(51,:,:));
figure,imagesc(rot90(slice)),axis equal,colormap jet

V = vol;
G2 = squeeze(flux2.data(:,:,:,2));
mcx = log10(G2);
mcx(mcx<0) = 0;
idx = find(mcx>0);
V(idx) = mcx(idx);
slice = squeeze(V(51,:,:));
figure,imagesc(rot90(slice)),axis equal,colormap jet

V = vol;
G3 = squeeze(flux2.data(:,:,:,3));
mcx = log10(G3);
mcx(mcx<0) = 0;
idx = find(mcx>0);
V(idx) = mcx(idx);
slice = squeeze(V(51,:,:));
figure,imagesc(rot90(slice)),axis equal,colormap jet