-
Notifications
You must be signed in to change notification settings - Fork 0
/
MC_test0806.m
177 lines (168 loc) · 9.52 KB
/
MC_test0806.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
function [cfg, det1, seeds] = MC_test0806(config, exitangle, para)
% 0.033mm grids
% DetInd1 mistake corrected
%% preparing the input data
% set seed to make the simulation repeatible
% cfg.seed=hex2dec('623F9A9E');
cfg.seed = randi(2^32);
cfg.nphoton=1e7;
cfg.unitinmm = 1/30; % define units in mm
cfg.size = [300 301 300]; % if unit = 0.05: 15mm cube
% cfg.size = [400 400 500]; % if unit = 0.005: 2x2x2.5mm cube
% time-domain simulation parameters
cfg.tstart=0;
cfg.tend=5e-9;
cfg.tstep=5e-9;
cfg.issrcfrom0 = 0;
% ========== define source and detector ==========
% cfg.srctype = 'pencil';
cfg.srctype = 'cone';
NA = 0.03;
cfg.srcparam1=[asin(NA) 0 0 0];
% cfg.srcparam1 = 100; %radius for the disk
cfg.srcpos=[cfg.size(1)/2 + 1, cfg.size(2)/2 + 1, 0.99]; % floor this x/y number to get the index of the grid (center = 125.5, ind = 125, in this case)
cfg.srcdir=[0 0 1];
cfg.detpos=[cfg.size(1)/2 + 1, cfg.size(2)/2 + 1, 0, min(cfg.size(1)/sqrt(2), cfg.size(2)/sqrt(2))]; % [x y z radius]
cfg.maxdetphoton = cfg.nphoton;
% ========== define a 3 layer structure ==========
cfg.config = config;
switch cfg.config
% ========== no skull ==========
case 'without skull' % 1: gray; 2: white
cfg.vol = 2.*ones(cfg.size); % 15 mm volume
cfg.bounds = round(para.gm_thick/cfg.unitinmm);
cfg.vol(:,:,cfg.bounds+1:end)= 3; % the rest: white matter
cfg.vol=uint8(cfg.vol);
case 'without skull off focus' %
cfg.srcpos=[cfg.size(1)/2-0.5, cfg.size(2)/2-0.5, -round(0.5/cfg.unitinmm)];
cfg.vol = 2.*ones(cfg.size); % 15 mm volume
cfg.bounds = round(para.gm_thick/cfg.unitinmm);
cfg.vol(:,:,cfg.bounds+1:end)= 3; % the rest: white matter
cfg.vol=uint8(cfg.vol);
% ========== with skull ==========
case 'with skull'
cfg.vol = ones(cfg.size); % 15 mm volume
cfg.bounds = [round(para.skull_thick/cfg.unitinmm) round(para.skull_thick/cfg.unitinmm)+round(para.gm_thick/cfg.unitinmm)]; % skull, gray matter, white matter
cfg.vol(:,:,cfg.bounds(1)+1:cfg.bounds(2))=2; % 0.5mm: skull
cfg.vol(:,:,cfg.bounds(2)+1:end)=3; % the rest: white matter
cfg.vol=uint8(cfg.vol);
case 'with skull off focus'
cfg.srcpos=[cfg.size(1)/2-0.5, cfg.size(2)/2-0.5, -round(0.5/cfg.unitinmm)];
cfg.vol = ones(cfg.size); % 15 mm volume
cfg.bounds = [round(para.skull_thick/cfg.unitinmm) round(para.skull_thick/cfg.unitinmm)+round(para.gm_thick/cfg.unitinmm)]; % skull, gray matter, white matter
cfg.vol(:,:,cfg.bounds(1)+1:cfg.bounds(2))=2; % 0.5mm: skull, 1.3 gray matter
cfg.vol(:,:,cfg.bounds(2)+1:end)=3; % the rest: white matter
cfg.vol=uint8(cfg.vol);
otherwise
end
cfg.isreflect = 1; % disable reflection at exterior boundary
% GPU thread configuration
cfg.autopilot=1;
cfg.gpuid=1;
% ========== output control ==========
cfg.savedetflag = 'dpxv';
cfg.debuglevel = 'M';
cfg.maxjumpdebug = cfg.nphoton; % number of trajectory positions saved
%% forward simulation
f1 = cell(1,length(para.ind));
seeds = f1;
det1 = f1;
traj = f1;
%%
cfg.perturb = nan; % a number or nan (perturb by 10%)
% cfg.exitangle = [0 0.03];
cfg.exitangle = exitangle;
% [0.8548 0.8844] angled with 0.87NA, deviation = asin(0.03) (±1.7191
% degrees), 0.87 determined by objective radius = 16.32, working distance =
% 9.25
% [0.5757 0.6237] angled with 0.6NA, deviation = asin(0.03) (±1.7191 degrees)
figure('DefaultAxesFontSize',18, 'DefaultLineLineWidth', 2,'color','w')
for i = 1:length(para.ind)
% for i = 2
% switch cfg.config
% case {'without skull', 'without skull off focus'}
% if isnan(cfg.perturb)
% cfg.prop = [0 0 1 1 % medium 0: the environment
% para.brain_mua(para.ind(i)) para.gray_matter_mus(para.ind(i)) para.gray_matter_g(para.ind(i)), 1.37 % medium 1: gray matter
% para.brain_mua(para.ind(i)) para.white_matter_mus(para.ind(i)) para.white_matter_g(para.ind(i)), 1.37 % medium 2: white matter
% para.brain_mua(para.ind(i))*1.1 para.gray_matter_mus(para.ind(i)) para.gray_matter_g(para.ind(i)), 1.37 % medium 3: gray matter, perturbed
% para.brain_mua(para.ind(i))*1.1 para.white_matter_mus(para.ind(i)) para.white_matter_g(para.ind(i)), 1.37]; % medium 4: white matter, pertubed
% else
% cfg.prop = [0 0 1 1 % medium 0: the environment
% para.brain_mua(para.ind(i)) para.gray_matter_mus(para.ind(i)) para.gray_matter_g(para.ind(i)), 1.37 % medium 1: gray matter
% para.brain_mua(para.ind(i)) para.white_matter_mus(para.ind(i)) para.white_matter_g(para.ind(i)), 1.37 % medium 2: white matter
% para.brain_mua(para.ind(i))+cfg.perturb para.gray_matter_mus(para.ind(i)) para.gray_matter_g(para.ind(i)), 1.37 % medium 3: gray matter, perturbed
% para.brain_mua(para.ind(i))+cfg.perturb para.white_matter_mus(para.ind(i)) para.white_matter_g(para.ind(i)), 1.37]; % medium 4: white matter, perturbed
% end
% case {'with skull', 'with skull off focus'}
if isnan(cfg.perturb)
cfg.prop=[0 0 1 1 % medium 0: the environment
para.skull_mua(para.ind(i)) para.skull_mus(para.ind(i)) para.skull_g, 1.56 % medium 1: skull
para.brain_mua(para.ind(i)) para.gray_matter_mus(para.ind(i)) para.gray_matter_g(para.ind(i)), 1.37 % medium 2: gray matter
para.brain_mua(para.ind(i)) para.white_matter_mus(para.ind(i)) para.white_matter_g(para.ind(i)), 1.37 % medium 3: white matter
para.skull_mua(para.ind(i))*1.1 para.skull_mus(para.ind(i)) para.skull_g, 1.56 % medium 4: skull, perturbed
para.brain_mua(para.ind(i))*1.1 para.gray_matter_mus(para.ind(i)) para.gray_matter_g(para.ind(i)), 1.37 % medium 5: gray matter, perturbed
para.brain_mua(para.ind(i))*1.1 para.white_matter_mus(para.ind(i)) para.white_matter_g(para.ind(i)), 1.37]; % medium 6: white matter, perturbed
else
cfg.prop=[0 0 1 1 % medium 0: the environment
para.skull_mua(para.ind(i)) para.skull_mus(para.ind(i)) para.skull_g, 1.56 % medium 1: skull
para.brain_mua(para.ind(i)) para.gray_matter_mus(para.ind(i)) para.gray_matter_g(para.ind(i)), 1.37 % medium 2: gray matter
para.brain_mua(para.ind(i)) para.white_matter_mus(para.ind(i)) para.white_matter_g(para.ind(i)), 1.37 % medium 3: white matter
para.skull_mua(para.ind(i))+cfg.perturb para.skull_mus(para.ind(i)) para.skull_g, 1.56 % medium 4: skull, perturbed
para.brain_mua(para.ind(i))+cfg.perturb para.gray_matter_mus(para.ind(i)) para.gray_matter_g(para.ind(i)), 1.37 % medium 5: gray matter, perturbed
para.brain_mua(para.ind(i))+cfg.perturb para.white_matter_mus(para.ind(i)) para.white_matter_g(para.ind(i)), 1.37]; % medium 6: white matter, perturbed
end
% otherwise
% disp('Check property setting (cfg.prop)')
% end
fprintf('running simulation ... this takes about 35 seconds on a GTX 470\n');
tic;
% f1=mcxlab(cfg);
[f1{i},det1{i},vol,seeds{i},traj{i}]=mcxlab(cfg);
toc;
% ===== show light distribution at the surface ====
% surface = squeeze(sum(f1{i}.data(:,:,1,:),4))';
% surface_profile = surface(:,round(cfg.size(2)/2));
% ========== detected photons ==========
ExitAngle = vecnorm(det1{i}.v(:,1:2)'); % sin_theta (i.e. NA), v is a normal vector
% det1{i}.DetInd1 = find(det1{i}.p(:,3)<round(0.5/cfg.unitinmm) & det1{i}.v(:,3)<0); % collect photons only exit from the 0-plane
det1{i}.DetInd1 = find(det1{i}.p(:,3)<0 & det1{i}.v(:,3)<0); % collect photons only exit from the 0-plane
det1{i}.DetInd2 = find(ExitAngle < cfg.exitangle(2) & ExitAngle >= cfg.exitangle(1)); % exit angle selection
det1{i}.DetInd = intersect(det1{i}.DetInd1, det1{i}.DetInd2);
det1{i}.DetNumber = length(det1{i}.DetInd);
det1{i}.DetExitPos = det1{i}.p(det1{i}.DetInd,1:2);
% ========== plot figure ==========
subplot(2,3,i);
I1 = log10(squeeze(sum(f1{i}.data(:,round(cfg.size(2)./2),:,:),4))');
I1(I1<0) = 0;
% if i == 1
cmax1 = max(I1(:)); cmin1 = min(I1(:));
% end
contourf(I1, cmin1:0.5:cmax1);
axis square
hold on
for k = 1:length(cfg.bounds)
plot([0 size(cfg.vol,1)],[cfg.bounds(k) cfg.bounds(k)],'--');
end
title(['flux ', cfg.config, ', wavelength ',num2str(para.WLs(i)),'nm'],'fontsize',14);
set(gca,'clim',[cmin1 cmax1]);
colorbar
drawnow
% subplot(3,6,6+i);
% h1 = histogram2(DetExitPosition(:,1),DetExitPosition(:,2),100,...
% 'XBinLimits',[cfg.detpos(1)-10, cfg.detpos(1)+10],'YBinLimits',[cfg.detpos(2)-10, cfg.detpos(2)+10],...
% 'DisplayStyle','tile','ShowEmptyBins','on');
% DetExitProfile = sum(h1.Values,2);
% [mm, ii] = max(DetExitProfile); cutoff = mm./2;
% [~,xx1] = min(abs(DetExitProfile(1:ii)-cutoff));
% [~,xx2] = min(abs(DetExitProfile(ii:end)-cutoff));
% DetExitSize = (xx2 + ii - xx1).*cfg.unitinmm;
% title(['Exit position, FWHM: ',num2str(DetExitSize,3),'mm'])
%
% subplot(3,6,12+i);
% DetWeight = mcxdetweight(det1,cfg.prop);
% DetPathlength = det1.ppath.*repmat(DetWeight(:),1,size(det1.ppath,2));
% DetMeanPathlength = cfg.unitinmm*sum(DetPathlength(:,2)) / sum(DetWeight(:));
% h2 = histogram(DetPathlength);
% title(['Pathlength distribution, mean: ',num2str(DetMeanPathlength,3),'mm'])
end