Skip to content

Commit

Permalink
add benchmark scripts for mmcl
Browse files Browse the repository at this point in the history
  • Loading branch information
@fangq
fangq committed Jul 1, 2019
1 parent c1ec5f1 commit f5aeac9
Show file tree
Hide file tree
Showing 4 changed files with 411 additions and 0 deletions.
107 changes: 107 additions & 0 deletions mmclab/example/demo_mmcl_b1_b1d.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,107 @@
%%-----------------------------------------------------------------
%% validate MMC with a homogeneous cubic domain
% (see MMC paper Fig. 2)
%%-----------------------------------------------------------------
%
% In this example, we validate the MMCM algorithm using a homogeneous
% cubic domain. This simulation generates the results
% shown in Fig. 2 in the paper.
%
% The cubic domain has a dimension of 60x60x60 mm with optical properties
% mua=0.001, mus=1, n=1.0 and g=0.01. The analytical solution
% can be computed by the cwdiffusion (for CW solution) and
% tddiffusion (for time-domain solutions) functions in the
% package of Monte Carlo eXtreme (MCX) under the mcx/utils directory.
%
%%-----------------------------------------------------------------

%addpath('/Please/add/path/to/mcx/utils/')
addpath('../../matlab/');

clear cfg
cfg.nphoton=1e7;
cfg.seed=1648335518;
[cfg.node,cfg.elem]=genT5mesh(0:2:60,0:2:60,0:2:60);
cfg.elemprop=ones(size(cfg.elem,1),1);
cfg.srcpos=[30.1,30.2,0];
cfg.srcdir=[0 0 1];
cfg.tstart=0;
cfg.tend=5e-9;
cfg.tstep=1e-10;
cfg.prop=[0 0 1 1;0.005 1.0 0.01 1.0];
cfg.debuglevel='TP';
cfg.isreflect=0;
cfg.method='elem';

cfg.gpuid=1;

b1gpu=mmclab(cfg);
b1gpu=sum(b1gpu.data,2);

%%
cfg.gpuid=-1;

b1cpu=mmclab(cfg);
b1cpu=sum(b1cpu.data,2);


%%-----------------------------------------------------------------
%% add paths to the necessary toolboxes
%%-----------------------------------------------------------------

c0=299792458000;

twin=[cfg.tstart+cfg.tstep/2:cfg.tstep:cfg.tend];
gates=length(twin);

[xi,yi]=meshgrid(0:2:60,0:2:60);

[cutpos1,cutvalue1,facedata1]=qmeshcut(cfg.elem(:,1:4),cfg.node,b1cpu,[0 30.2 0; 0 30.2 1; 1 30.2 0]);
[cutpos2,cutvalue2,facedata2]=qmeshcut(cfg.elem(:,1:4),cfg.node,b1gpu,[0 30.2 0; 0 30.2 1; 1 30.2 0]);
Phicpu=griddata(cutpos1(:,1),cutpos1(:,3),cutvalue1,xi,yi);
Phigpu=griddata(cutpos2(:,1),cutpos2(:,3),cutvalue2,xi,yi);


%% DMMC

[cfg.node,cfg.elem]=meshgrid6(0:60:60,0:60:60,0:60:60);
cfg.elem(:,1:4)=meshreorient(cfg.node,cfg.elem(:,1:4));

cfg.elemprop=ones(size(cfg.elem,1),1);
cfg.method='grid';

cfg.gpuid=1;

flux=mmclab(cfg);

Phidgpu=sum(flux.data,4);

%% DMMC
cfg.gpuid=-1;

flux=mmclab(cfg);
Phidcpu=sum(flux.data,4);

%%-----------------------------------------------------------------
%% generate a contour plot along y=30.2
%%-----------------------------------------------------------------
figure
clines = 0:-0.5:-8;
[xi,yi]=meshgrid(0:2:60,0:2:60);
srcs=[30.1,30.2,0];
dets=[xi(:) 30.2*ones(size(xi(:))) yi(:)];

hold on
[c h2]=contourf(xi,yi, log10(max(squeeze(Phicpu*cfg.tstep),1e-8)), clines, 'k-' );
contour(xi,yi,log10(abs(Phigpu*cfg.tstep)),clines,'r:')
contour(xi,yi,log10(abs(squeeze(Phidcpu(1:2:end,30,1:2:end))'*cfg.tstep)),clines,'b--')
contour(xi,yi,log10(abs(squeeze(Phidgpu(1:2:end,30,1:2:end))'*cfg.tstep)),clines,'y--')

axis equal
set(gca,'xlim',[1 60])
set(gca,'fontsize',20)
xlabel('x (mm)')
ylabel('z (mm)')
legend('MMC','MMCL','D-MMC','D-MMCL')
legend boxoff;
box on;
132 changes: 132 additions & 0 deletions mmclab/example/demo_mmcl_b2_b2d.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,132 @@
%%-----------------------------------------------------------------
%% validate DMMC with MMC in a heterogeneous cubic domain
% (see DMMC paper Fig. 1)
%%-----------------------------------------------------------------
%
% In this example, we validate the DMMC algorithm with MMC using a
% heterogeneous cubic domain. We also compare the DMMC/MMC solution with
% that of voxel-based MC (MCX). This simulation generates the results
% shown in Fig. 1(c-e) in the paper.
%
%%-----------------------------------------------------------------

% preparation

% 1. you need to add the path to iso2mesh toolbox
% addpath('/path/to/iso2mesh/toolbox/');

% 2. you need to add the path to MMC matlab folder
% addpath('../../matlab');

% create a surface mesh for a 10 mm radius sphere
[nsph_10,fsph_10]=meshasphere([30 30 30],10,1.0);
[nsph_10,fsph_10]=removeisolatednode(nsph_10,fsph_10);

% create a surface mesh for a 23 mm radius sphere
[nsph_23,fsph_23]=meshasphere([30 30 30],23,2.0);
[nsph_23,fsph_23]=removeisolatednode(nsph_23,fsph_23);

% create a surface mesh for a 25 mm radius sphere
[nsph_25,fsph_25]=meshasphere([30 30 30],25,2.0);
[nsph_25,fsph_25]=removeisolatednode(nsph_25,fsph_25);

%%-----------------------------------------------------------------
%% create simulation parameters
%%-----------------------------------------------------------------
clear cfg

cfg.nphoton=3e6;
cfg.seed=1648335518;
cfg.srcpos=[30,30.1,0];
cfg.srcdir=[0 0 1];
cfg.tstart=0;
cfg.tend=5e-9;
cfg.tstep=5e-10;
cfg.prop=[0,0,1,1;0.02 7.0 0.89 1.37;0.004 0.009 0.89 1.37;0.02 9.0 0.89 1.37;0.05 0.0 1.0 1.37];
cfg.debuglevel='TP';
cfg.isreflect=1;

cfgs=[cfg cfg];

%%-----------------------------------------------------------------
%% tetrahedral mesh generation
%%-----------------------------------------------------------------

% generate a coarse CDT volumetric mesh from triangular surfaces of
% three spheres with an additional bounding box for DMMC

ISO2MESH_SESSION='dmmc1_';

[nbox,ebox]=meshgrid6(0:60:60,0:60:60,0:60:60);
fbox=volface(ebox);
[no,fc]=mergemesh(nsph_10,fsph_10,nsph_23,fsph_23,nsph_25,fsph_25,nbox,fbox);
[no,fc]=removeisolatednode(no,fc);

ISO2MESH_TETGENOPT='-A -q -Y'
[cfgs(1).node,cfgs(1).elem]=surf2mesh(no,fc,[0 0 0],[60.1 60.1 60.1],1,100,[1 1 1;30 30 6;30 30 15;30 30 30]);%thin layer
% [cfgs(1).node,cfgs(1).elem]=surf2mesh(no,fc,[0 0 0],[60.1 60.1 60.1],1,100,[1 1 1;30 30 6;30 30 17;30 30 30]);%thick layer
cfgs(1).method='grid';

% generate a refined volumetric mesh from triangular surfaces of
% three spheres with an additional bounding box for MMC
clear ISO2MESH_TETGENOPT
ISO2MESH_SESSION='dmmc2_';

[no,fc]=mergemesh(nsph_10,fsph_10,nsph_23,fsph_23,nsph_25,fsph_25);
[no,fc]=removeisolatednode(no,fc);
srcpos=cfgs(1).srcpos;
[cfgs(2).node,cfgs(2).elem,face2]=surf2mesh(no,fc,[0 0 0],[60 60 60],1,0.25,[1 1 1 0.1;30 30 6 0.1;30 30 17 0.1;30 30 30 0.1],[],[1 1 1 1 1 1 1 1]);
[cfgs(2).node,cfgs(2).elem]=sortmesh(srcpos,cfgs(2).node,cfgs(2).elem,1:4);
cfgs(2).method='elem';

%%

cfgs(1).gpuid=1;
b2dgpu=mmclab(cfgs(1));
b2dgpu=sum(b2dgpu.data,4);

%%
cfgs(1).gpuid=-1;
b2dcpu=mmclab(cfgs(1));
b2dcpu=sum(b2dcpu.data,4);

%%
cfgs(2).gpuid=1;
b2gpu=mmclab(cfgs(2));
b2gpu=sum(b2gpu.data,2);

%%
cfgs(2).gpuid=-1;
b2cpu=mmclab(cfgs(2));
b2cpu=sum(b2cpu.data,2);

%%
[xi,yi]=meshgrid(0:2:60,0:2:60);
[cutpos1,cutvalue1,facedata1]=qmeshcut(cfgs(2).elem(:,1:4),cfgs(2).node,b2cpu,[0 30.2 0; 0 30.2 1; 1 30.2 0]);
[cutpos2,cutvalue2,facedata2]=qmeshcut(cfgs(2).elem(:,1:4),cfgs(2).node,b2gpu,[0 30.2 0; 0 30.2 1; 1 30.2 0]);
b2cpu=griddata(cutpos1(:,1),cutpos1(:,3),cutvalue1,xi,yi);
b2gpu=griddata(cutpos2(:,1),cutpos2(:,3),cutvalue2,xi,yi);

%%-----------------------------------------------------------------
%% generate a contour plot along y=30.2
%%-----------------------------------------------------------------
figure
clines = 0:-0.5:-8;
[xi,yi]=meshgrid(0:2:60,0:2:60);
srcs=[30.1,30.2,0];
dets=[xi(:) 30.2*ones(size(xi(:))) yi(:)];

hold on
[c h2]=contourf(xi,yi, log10(max(squeeze(b2cpu*cfgs(1).tstep),1e-8)), clines, 'k-' );
contour(xi,yi,log10(abs(b2gpu*cfgs(1).tstep)),clines,'r:')
contour(xi,yi,log10(abs(squeeze(b2dcpu(1:2:end,30,1:2:end))'*cfgs(1).tstep)),clines,'b--')
contour(xi,yi,log10(abs(squeeze(b2dgpu(1:2:end,30,1:2:end))'*cfgs(1).tstep)),clines,'y--')

axis equal
set(gca,'xlim',[1 60])
set(gca,'fontsize',20)
xlabel('x (mm)')
ylabel('z (mm)')
legend('MMC','MMCL','D-MMC','D-MMCL')
legend boxoff;
box on;
75 changes: 75 additions & 0 deletions mmclab/example/demo_mmcl_b3.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,75 @@
addpath('../../../matlab/');

if(~exist('MMC_Collins_Atlas_Mesh_Version_2L.mat','file'))
if(~exist('MMC_Collins_Atlas_Mesh_Version_2L.tar.gz','file'))
urlwrite('http://downloads.sourceforge.net/project/mcx/mmc/AdultBrain%20Atlas%20FEM%20Mesh/Version%202/MMC_Collins_Atlas_Mesh_Version_2L.tar.gz?r=http%3A%2F%2Fsourceforge.net%2Fprojects%2Fmcx%2Ffiles%2Fmmc%2FAdultBrain%2520Atlas%2520FEM%2520Mesh%2FVersion%25202%2F&ts=1451344236&use_mirror=tcpdiag','MMC_Collins_Atlas_Mesh_Version_2L.tar.gz');
end
if(~exist('MMC_Collins_Atlas_Mesh_Version_2L.tar','file'))
gunzip('MMC_Collins_Atlas_Mesh_Version_2L.tar.gz');
end
if(~exist('MMC_Collins_Atlas_Mesh_Version_2L.mat','file'))
untar('MMC_Collins_Atlas_Mesh_Version_2L.tar');
end
end

load MMC_Collins_Atlas_Mesh_Version_2L.mat
cfg.node=node;
cfg.elem=elem(:,1:4);
cfg.elemprop=elem(:,5);

cfg.tstart=0;
cfg.tend=5e-09;
cfg.tstep=2e-10;

cfg.srcpos=[75 67.38 167.5];
cfg.srcdir=[0.1636 0.4569 -0.8743];
cfg.srcdir=cfg.srcdir/norm(cfg.srcdir);

cfg.detpos=[ 75.0000 77.1900 170.3000 3.0000
75.0000 89.0000 171.6000 3.0000
75.0000 97.6700 172.4000 3.0000
75.0000 102.4000 172.0000 3.0000];

cfg.prop=[ 0 0 1.0000 1.0000 % background/air
0.0190 7.8182 0.8900 1.3700 % scalp
0.0190 7.8182 0.8900 1.3700 % skull
0.0040 0.0090 0.8900 1.3700 % csf
0.0200 9.0000 0.8900 1.3700 % gray matters
0.0800 40.9000 0.8400 1.3700 % white matters
0 0 1.0000 1.0000]; % air pockets

cfg.seed=29012392;
cfg.nphoton=1e7;
cfg.method='grid';

%%
cfg.gpuid=1;

b1gpu=mmclab(cfg);
b1gpu=sum(b1gpu.data,4);
%%
cfg.gpuid=-1;

b1cpu=mmclab(cfg);
b1cpu=sum(b1cpu.data,4);


%%-----------------------------------------------------------------
%% generate a contour plot along y=30.2
%%-----------------------------------------------------------------
figure
clines = 0:-0.5:-8;

hold on
[c h2]=contourf(log10(squeeze(b1cpu*cfg.tstep)), clines, 'k-' );
contour(log10(abs(b1gpu*cfg.tstep)),clines,'r:')

axis equal
set(gca,'xlim',[1 60])
set(gca,'fontsize',20)
xlabel('x (mm)')
ylabel('z (mm)')
legend('MMC','MMCL')
legend boxoff;
box on;

Loading

0 comments on commit f5aeac9

Please sign in to comment.