Skip to content

Commit

Permalink
[format] reformat all MATLAB codes with miss_hit
Browse files Browse the repository at this point in the history
  • Loading branch information
@fangq
fangq committed Feb 29, 2024
1 parent ce1d374 commit 580cb27
Show file tree
Hide file tree
Showing 53 changed files with 2,380 additions and 2,366 deletions.
163 changes: 81 additions & 82 deletions example/validation/plotsimudata.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,128 +2,127 @@
% matlab script to produce the validation plots (Fig. 4)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if(exist('loadjson')~=2)
error('you need to first download jsonlab (https:/fangq/jsonlab) and add it to the path in order to run this script');
if (exist('loadjson') ~= 2)
error('you need to first download jsonlab (https:/fangq/jsonlab) and add it to the path in order to run this script');
end

c0=299792458000;
srcpos=[30 30 1];
detpos=[30 15 10];
twin=[5e-11:1e-10:5e-9];
c0 = 299792458000;
srcpos = [30 30 1];
detpos = [30 15 10];
twin = [5e-11:1e-10:5e-9];

%----------------------------------------------------------------
% ----------------------------------------------------------------
% Plot time-domain comparison at a given position
%----------------------------------------------------------------
% ----------------------------------------------------------------

figure;
hold on;

mcx=loadjson('semi_infinite.jnii',[60 60 60 50]);
semilogy((1:50)/10,tddiffusion(0.005, 1, c0, 0, srcpos, detpos,twin),'r');
semilogy((1:50)/10,squeeze(mcx(30,14,9,:)),'o');
mcx = loadjson('semi_infinite.jnii', [60 60 60 50]);
semilogy((1:50) / 10, tddiffusion(0.005, 1, c0, 0, srcpos, detpos, twin), 'r');
semilogy((1:50) / 10, squeeze(mcx(30, 14, 9, :)), 'o');

mcxb=loadjson('semi_infinite_b.jnii',[60 60 60 50]);
semilogy((1:50)/10,tddiffusion(0.005, 1, c0/1.37, 0.493, srcpos, detpos,twin),'r--');
semilogy((1:50)/10,squeeze(mcxb(30,14,9,:)),'+');
mcxb = loadjson('semi_infinite_b.jnii', [60 60 60 50]);
semilogy((1:50) / 10, tddiffusion(0.005, 1, c0 / 1.37, 0.493, srcpos, detpos, twin), 'r--');
semilogy((1:50) / 10, squeeze(mcxb(30, 14, 9, :)), '+');

set(gca,'ylim',[1 10^8.5]);
set(gca, 'ylim', [1 10^8.5]);

legend('Diffusion (no reflection)','MCX (no reflection)','Diffusion (with reflection)','MCX (with reflection)');
legend('Diffusion (no reflection)', 'MCX (no reflection)', 'Diffusion (with reflection)', 'MCX (with reflection)');
legend boxoff;
set(gca,'yscale','log');
set(gca,'fontsize',18);
xlabel('time (ns)')
ylabel('Fluence in 1/(m^2s)')
set(gca, 'yscale', 'log');
set(gca, 'fontsize', 18);
xlabel('time (ns)');
ylabel('Fluence in 1/(m^2s)');
box on;
saveas(gcf,'td_validate.fig')
print -depsc2 td_validate.eps
saveas(gcf, 'td_validate.fig');
print -depsc2 td_validate.eps;

%----------------------------------------------------------------
% Plot CW comparison along a line
%----------------------------------------------------------------
% ----------------------------------------------------------------
% Plot CW comparison along a line
% ----------------------------------------------------------------

figure;
hold on;

cwmcx=sum(mcx,4);
srcpos=[0 0 0];
detpos=[0.5+(0:29)', 0.5*ones(30,1),0.5*ones(30,1)];
phicw=cwdiffusion(0.005, 1, 0, srcpos, detpos);
phimcx=cwmcx(30,:,1);
cwmcx = sum(mcx, 4);
srcpos = [0 0 0];
detpos = [0.5 + (0:29)', 0.5 * ones(30, 1), 0.5 * ones(30, 1)];
phicw = cwdiffusion(0.005, 1, 0, srcpos, detpos);
phimcx = cwmcx(30, :, 1);

cwmcxb=sum(mcxb,4);
phicwb=cwdiffusion(0.005, 1, 0.493, srcpos, detpos);
phimcxb=(cwmcxb(30,:,1));
h=semilogy(1:30, phicw, 'r-', 1:30, phimcx(30:59)*1e-10, 'o', 1:30, phicwb, 'r--', 1:30, phimcxb(30:59)*1e-10, 'b+');
cwmcxb = sum(mcxb, 4);
phicwb = cwdiffusion(0.005, 1, 0.493, srcpos, detpos);
phimcxb = (cwmcxb(30, :, 1));
h = semilogy(1:30, phicw, 'r-', 1:30, phimcx(30:59) * 1e-10, 'o', 1:30, phicwb, 'r--', 1:30, phimcxb(30:59) * 1e-10, 'b+');
box on;
set(gca,'yscale','log')
set(gca, 'yscale', 'log');

legend('Diffusion (no reflection)','MCX (no reflection)','Diffusion (with reflection)','MCX(with reflection)');
legend('Diffusion (no reflection)', 'MCX (no reflection)', 'Diffusion (with reflection)', 'MCX(with reflection)');
legend boxoff;
set(gca,'fontsize',18);
xlabel('x (mm)')
ylabel('Fluence in 1/(m^2s)')
saveas(gcf,'td_validate_cw.fig')
print -depsc2 td_validate_cw.eps
set(gca, 'fontsize', 18);
xlabel('x (mm)');
ylabel('Fluence in 1/(m^2s)');
saveas(gcf, 'td_validate_cw.fig');
print -depsc2 td_validate_cw.eps;

%----------------------------------------------------------------
% ----------------------------------------------------------------
% Plot CW comparison at a cross-section
%----------------------------------------------------------------
% ----------------------------------------------------------------

figure
figure;
hold on;

clines = -1.5:-0.5:-4;
srcpos=[29.5 29.5 0.5];
[xi,yi]=meshgrid(0.5:59.5,0.5:59.5);
detpos=[xi(:) 30*ones(size(xi(:))) yi(:)];
phicw=reshape(cwdiffusion(0.005, 1.0, 0, srcpos, detpos),size(xi));

[c h2]=contour(xi,yi, log10(max(squeeze(phicw),1e-8)), clines, 'k-' );
[c h1]=contour(log10(max(squeeze(cwmcx(:,30,:)*1e-10)',1e-8)), clines, 'k:' );
legend('Diffusion','MCX')
srcpos = [29.5 29.5 0.5];
[xi, yi] = meshgrid(0.5:59.5, 0.5:59.5);
detpos = [xi(:) 30 * ones(size(xi(:))) yi(:)];
phicw = reshape(cwdiffusion(0.005, 1.0, 0, srcpos, detpos), size(xi));

[c h2] = contour(xi, yi, log10(max(squeeze(phicw), 1e-8)), clines, 'k-');
[c h1] = contour(log10(max(squeeze(cwmcx(:, 30, :) * 1e-10)', 1e-8)), clines, 'k:');
legend('Diffusion', 'MCX');
legend boxoff;
box on;
set(gca,'fontsize',18);
xlabel('x (mm)')
ylabel('z (mm)')
saveas(gcf,'td_validate_cw_2d.fig')
print -depsc2 td_validate_cw_2d.eps
set(gca, 'fontsize', 18);
xlabel('x (mm)');
ylabel('z (mm)');
saveas(gcf, 'td_validate_cw_2d.fig');
print -depsc2 td_validate_cw_2d.eps;

%----------------------------------------------------------------
% ----------------------------------------------------------------
% Plot time-domain comparison at a cross-section
%----------------------------------------------------------------
% ----------------------------------------------------------------

c0=299792458000;
twin=5e-11:1e-10:5e-9;
mua=0.005;
musp=1;
Reff=0;
c0 = 299792458000;
twin = 5e-11:1e-10:5e-9;
mua = 0.005;
musp = 1;
Reff = 0;

srcpos=[29.5 29.5 0];
[xi,yi]=meshgrid(0.5:59.5,0.5:59.5);
srcpos = [29.5 29.5 0];
[xi, yi] = meshgrid(0.5:59.5, 0.5:59.5);

detpos=[xi(:) 29.5*ones(size(xi(:))) yi(:)];
for i=1:length(twin)
phitd(:,:,i)=reshape(tddiffusion(mua,musp,c0,Reff,srcpos,detpos,twin(i)),size(xi));
detpos = [xi(:) 29.5 * ones(size(xi(:))) yi(:)];
for i = 1:length(twin)
phitd(:, :, i) = reshape(tddiffusion(mua, musp, c0, Reff, srcpos, detpos, twin(i)), size(xi));
end

figure;
hold on
hold on;

for i=1:5:21
halfmax=max(max(phitd(:,:,i)))/2;
[c h2]=contour(xi, yi,max(squeeze(phitd(:,:,i)),1e-12),[halfmax halfmax], 'k-');
halfmax=max(max(mcx(:,30,:,i)))/2;
[c h2]=contour(max(squeeze(mcx(:,30,:,i))',1e-12),[halfmax halfmax], 'k:');
for i = 1:5:21
halfmax = max(max(phitd(:, :, i))) / 2;
[c h2] = contour(xi, yi, max(squeeze(phitd(:, :, i)), 1e-12), [halfmax halfmax], 'k-');
halfmax = max(max(mcx(:, 30, :, i))) / 2;
[c h2] = contour(max(squeeze(mcx(:, 30, :, i))', 1e-12), [halfmax halfmax], 'k:');
end
legend('Diffusion','MCX')
legend('Diffusion', 'MCX');
legend boxoff;
box on;
set(gca,'fontsize',18);
xlabel('x (mm)')
ylabel('z (mm)')
saveas(gcf,'td_validate_2d.fig')

print -depsc2 td_validate_2d.eps
set(gca, 'fontsize', 18);
xlabel('x (mm)');
ylabel('z (mm)');
saveas(gcf, 'td_validate_2d.fig');

print -depsc2 td_validate_2d.eps;
104 changes: 52 additions & 52 deletions mcxlabcl/examples/demo_4layer_head.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
% MCXLAB - Monte Carlo eXtreme for MATLAB/Octave by Qianqina Fang
%
% In this example, we simulate a 4-layer brain model using MCXLAB.
% We will investigate the differences between the solutions with and
% We will investigate the differences between the solutions with and
% witout boundary reflections (both external and internal) and show
% you how to display and analyze the resulting data.
%
Expand All @@ -14,94 +14,94 @@

%% preparing the input data
% set seed to make the simulation repeatible
cfg.seed=hex2dec('623F9A9E');
cfg.seed = hex2dec('623F9A9E');

cfg.nphoton=5e7;
cfg.nphoton = 5e7;

% define a 4 layer structure
cfg.vol=ones(100,100,50);
cfg.vol(:,:,21:25)=2;
cfg.vol(:,:,26:35)=3;
cfg.vol(:,:,36:end)=4;
cfg.vol=uint8(cfg.vol);
cfg.vol = ones(100, 100, 50);
cfg.vol(:, :, 21:25) = 2;
cfg.vol(:, :, 26:35) = 3;
cfg.vol(:, :, 36:end) = 4;
cfg.vol = uint8(cfg.vol);

% define the source position
cfg.srcpos=[50,50,0]+1;
cfg.srcdir=[0 0 1];
cfg.srcpos = [50, 50, 0] + 1;
cfg.srcdir = [0 0 1];

% use the brain optical properties defined at
% http://mcx.sourceforge.net/cgi-bin/index.cgi?MMC/CollinsAtlasMesh
% format: [mua(1/mm) mus(1/mm) g n]

cfg.prop=[0 0 1 1 % medium 0: the environment
0.019 7.8 0.89 1.37 % medium 1: skin & skull
0.004 0.009 0.89 1.37 % medium 2: CSF
0.02 9.0 0.89 1.37 % medium 3: gray matter
0.08 40.9 0.84 1.37]; % medium 4: white matter
cfg.prop = [0 0 1 1 % medium 0: the environment
0.019 7.8 0.89 1.37 % medium 1: skin & skull
0.004 0.009 0.89 1.37 % medium 2: CSF
0.02 9.0 0.89 1.37 % medium 3: gray matter
0.08 40.9 0.84 1.37]; % medium 4: white matter

% time-domain simulation parameters
cfg.tstart=0;
cfg.tend=5e-9;
cfg.tstep=5e-10;
cfg.tstart = 0;
cfg.tend = 5e-9;
cfg.tstep = 5e-10;

% GPU thread configuration
cfg.autopilot=1;
cfg.gpuid=1;
cfg.autopilot = 1;
cfg.gpuid = 1;

%% running simulation without boundary reflection
fprintf('running simulation ... this takes about 35 seconds on a GTX 470\n');
cfg.isreflect=0; % disable reflection at exterior boundary
cfg.isreflect = 0; % disable reflection at exterior boundary
tic;
f1=mcxlabcl(cfg);
f1 = mcxlabcl(cfg);
toc;

%% running simulation with boundary reflection enabled
fprintf('running simulation ... this takes about 50 seconds on a GTX 470\n');
cfg.isreflect=1; % enable reflection at exterior boundary
cfg.isrefint=1; % enable reflection at interior boundary too
cfg.issavedet=1; % enable recording partial pathlength of detected photons
cfg.detpos=[31 51 1 2];
cfg.isreflect = 1; % enable reflection at exterior boundary
cfg.isrefint = 1; % enable reflection at interior boundary too
cfg.issavedet = 1; % enable recording partial pathlength of detected photons
cfg.detpos = [31 51 1 2];
tic;
[f2,det2]=mcxlabcl(cfg);
[f2, det2] = mcxlabcl(cfg);
toc;

%% plot the results
figure
figure;
subplot(221);
contourf(log10(squeeze(sum(f1.data(:,51,:,:),4))'),1:0.5:8);
hold on
plot([0 100],[21 21],'--',[0 100],[26 26],'--',[0 100],[36 36],'--');
contourf(log10(squeeze(sum(f1.data(:, 51, :, :), 4))'), 1:0.5:8);
hold on;
plot([0 100], [21 21], '--', [0 100], [26 26], '--', [0 100], [36 36], '--');
title('Fluence rate with no reflection');
set(gca,'clim',[1 8]);
set(gca, 'clim', [1 8]);

subplot(222);
contourf(log10(squeeze(sum(f2.data(:,51,:,:),4))'),1:0.5:8);
hold on
plot([0 100],[21 21],'--',[0 100],[26 26],'--',[0 100],[36 36],'--');
contourf(log10(squeeze(sum(f2.data(:, 51, :, :), 4))'), 1:0.5:8);
hold on;
plot([0 100], [21 21], '--', [0 100], [26 26], '--', [0 100], [36 36], '--');
title('Fluence rate with reflection at boundaries');
set(gca,'clim',[1 8]);
set(gca, 'clim', [1 8]);

subplot(223);
hold on;
xi=1e9*((cfg.tstart:cfg.tstep:cfg.tend-cfg.tstep)+0.5*cfg.tstep);
for i=31:49
semilogy(xi,squeeze(f2.data(51,i,1,:)),'color',[1-(i-25)/25 1-(i-25)/25 1]);
xi = 1e9 * ((cfg.tstart:cfg.tstep:cfg.tend - cfg.tstep) + 0.5 * cfg.tstep);
for i = 31:49
semilogy(xi, squeeze(f2.data(51, i, 1, :)), 'color', [1 - (i - 25) / 25 1 - (i - 25) / 25 1]);
end
xlabel('time (ns)')
ylabel('Fluence rate (1/mm^2/s)')
set(gca,'yscale','log');
xlabel('time (ns)');
ylabel('Fluence rate (1/mm^2/s)');
set(gca, 'yscale', 'log');
title('time point spread functions (TPSF)');

subplot(224);
hold on;
[hs1,c1]=hist(det2.ppath(find(det2.ppath(:,1)),1),200);
[hs2,c2]=hist(det2.ppath(find(det2.ppath(:,2)),2),200);
[hs3,c3]=hist(det2.ppath(find(det2.ppath(:,3)),3),200);
[hs4,c4]=hist(det2.ppath(find(det2.ppath(:,4)),4),200);
bar(c1,hs1,'edgecolor','none','facecolor','r');
bar(c2,hs2,'edgecolor','none','facecolor','b');
bar(c3,hs3,'edgecolor','none','facecolor','g');
bar(c4,hs4,'edgecolor','none','facecolor','k');
legend('tissue 1','tissue 2','tissue 3','tissue 4');
[hs1, c1] = hist(det2.ppath(find(det2.ppath(:, 1)), 1), 200);
[hs2, c2] = hist(det2.ppath(find(det2.ppath(:, 2)), 2), 200);
[hs3, c3] = hist(det2.ppath(find(det2.ppath(:, 3)), 3), 200);
[hs4, c4] = hist(det2.ppath(find(det2.ppath(:, 4)), 4), 200);
bar(c1, hs1, 'edgecolor', 'none', 'facecolor', 'r');
bar(c2, hs2, 'edgecolor', 'none', 'facecolor', 'b');
bar(c3, hs3, 'edgecolor', 'none', 'facecolor', 'g');
bar(c4, hs4, 'edgecolor', 'none', 'facecolor', 'k');
legend('tissue 1', 'tissue 2', 'tissue 3', 'tissue 4');
xlabel('partial pathlength (mm)');
title(sprintf('detected %d photons',size(det2.ppath,1)));
title(sprintf('detected %d photons', size(det2.ppath, 1)));
Loading

0 comments on commit 580cb27

Please sign in to comment.