clear cfg;
close all
clear all

% Load data (0.1 mm^3/voxel)
load digimouse;

%%%%%% TODO: simulated photon number: 1e6 photon %%%%%%
cfg.nphoton =

% 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];

% GPU id
cfg.gpuid=1;
% Automatically set threads and blocks
cfg.autopilot=1;

% Set starting time
cfg.tstart=0;
%%%%%% TODO: Set simulation end time (1 ns) %%%%%%
cfg.tend =
%%%%%% TODO: Time step (0.1 ns) %%%%%%
cfg.tstep =

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

%%%%%% TODO: define source type (pencil) %%%%%%
cfg.srctype =

%%%%%% TODO: define detector (x,y,z)=(26 46 50) with radius=2 %%%%%%
cfg.issaveexit = 1; % save exit position and direction of detected photons
cfg.detpos = 

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

%%%%%% TODO: compute the CW solution for flux1 %%%%%%
CW = 

% 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);
%% Visualization
vol = single(cfg.vol);
vol = log(vol)/log(4)+1;

V = vol;
G1 = squeeze(flux1.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(flux1.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(flux1.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