|
1 | | -clear cfg |
| 1 | +clear cfg; |
2 | 2 |
|
3 | | -disp('== mesh generation (iso2mesh) ...') |
4 | | -%# Create a box-like homogeneous domain |
5 | | -boxsize = [60, 50, 40]; %# domain size |
6 | | -trisize = 4; %# max triangule size on the surface |
7 | | -maxvol = 4; %# max tet element volume |
| 3 | +disp('== mesh generation (iso2mesh) ...'); |
| 4 | +%#Create a box-like homogeneous domain |
| 5 | +boxsize = [60, 50, 40]; %#domain size |
| 6 | +trisize = 4; %#max triangule size on the surface |
| 7 | +maxvol = 4; %#max tet element volume |
8 | 8 |
|
9 | | -tic |
10 | | -[nbox1, fbox1] = meshabox([0,0,0], boxsize, trisize); %# create a big box |
11 | | -[nbox2, fbox2] = meshabox([10,10,10], [30,30,30], trisize); %# create a box inclusion |
12 | | -[nbox1, fbox1] = removeisolatednode(nbox1, fbox1(:,1:3)); %# clean the surface |
13 | | -[nbox2, fbox2] = removeisolatednode(nbox2, fbox2(:,1:3)); |
| 9 | +tic; |
| 10 | +[nbox1, fbox1] = meshabox([0, 0, 0], boxsize, trisize); %#create a big box |
| 11 | +[nbox2, fbox2] = meshabox([10, 10, 10], [30, 30, 30], trisize); %#create a box inclusion |
| 12 | +[nbox1, fbox1] = removeisolatednode(nbox1, fbox1(:, 1:3)); %#clean the surface |
| 13 | +[nbox2, fbox2] = removeisolatednode(nbox2, fbox2(:, 1:3)); |
14 | 14 |
|
15 | | -[no1, fc1] = mergemesh(nbox1, fbox1, nbox2, fbox2); %# combine the two non-intersecting surfaces |
16 | | -regionseed = [1, 1, 1 %# a seed point inside region 1 (large box) |
17 | | - 11,11,11]; %# a seed point inside region 2 (small box) |
18 | | -[cfg.node, cfg.elem] = s2m(no1,fc1,1,maxvol, 'tetgen', regionseed); %# generate tetrahedral mesh - outer box: label 1, inner box: label 2 |
| 15 | +[no1, fc1] = mergemesh(nbox1, fbox1, nbox2, fbox2); %#combine the two non-intersecting surfaces |
| 16 | +regionseed = [1, 1, 1 %#a seed point inside region 1 (large box) |
| 17 | + 11, 11, 11]; %#a seed point inside region 2 (small box) |
| 18 | +[cfg.node, cfg.elem] = s2m(no1, fc1, 1, maxvol, 'tetgen', regionseed); %#generate tetrahedral mesh - outer box: label 1, inner box: label 2 |
19 | 19 |
|
20 | | -cfg.seg = cfg.elem(:,5); %# cfg.seg is similar to cfg.elemprop in mmclab, but also supports node-based labels |
21 | | -cfg.elem(:,5)=[]; |
22 | | -toc |
| 20 | +cfg.seg = cfg.elem(:, 5); %#cfg.seg is similar to cfg.elemprop in mmclab, but also supports node-based labels |
| 21 | +cfg.elem(:, 5) = []; |
| 22 | +toc; |
23 | 23 |
|
24 | | -plotmesh(cfg.node, [cfg.elem, cfg.seg], 'y>20') |
| 24 | +plotmesh(cfg.node, [cfg.elem, cfg.seg], 'y>20'); |
25 | 25 |
|
26 | | -disp('== define settings in cfg ...') |
| 26 | +disp('== define settings in cfg ...'); |
27 | 27 |
|
28 | | -%# Creating forward simulation data structure cfg |
29 | | -%# properties:[mua(1/mm),mus(1/mm),g, n] % if both mus/g are given, mus'=mus/(1-g) will be used for diffusion, usually set g to 0 |
30 | | -cfg.prop = [0, 0, 1, 1 %# cfg.prop is the same as mcx/mmc, row-1 is for label 0, row-2 for label 1 etc |
31 | | - 0.006, 0.8, 0, 1.37 %# label 1 (background domain) |
32 | | - 0.02, 1, 0, 1.37]; %# label 2 (inclusion) |
| 28 | +%#Creating forward simulation data structure cfg |
| 29 | +%#properties:[mua(1/mm),mus(1/mm),g, n] % if both mus/g are given, mus'=mus/(1-g) will be used for diffusion, usually set g to 0 |
| 30 | +cfg.prop = [0, 0, 1, 1 %#cfg.prop is the same as mcx/mmc, row-1 is for label 0, row-2 for label 1 etc |
| 31 | + 0.006, 0.8, 0, 1.37 %#label 1 (background domain) |
| 32 | + 0.02, 1, 0, 1.37]; %#label 2 (inclusion) |
33 | 33 |
|
34 | | -%# Default redbird source is a pencil beam, same as mcx/mmc |
35 | | -cfg.srcpos = [25, 25, 0]; %# redbird srcpos can have multiple rows |
36 | | -cfg.srcdir = [0, 0, 1]; %# srcdir determines how sources are sunken into the mesh |
| 34 | +%#Default redbird source is a pencil beam, same as mcx/mmc |
| 35 | +cfg.srcpos = [25, 25, 0]; %#redbird srcpos can have multiple rows |
| 36 | +cfg.srcdir = [0, 0, 1]; %#srcdir determines how sources are sunken into the mesh |
37 | 37 |
|
38 | | -%# Redbird detector positions are point-like, directly sampling the output fluence; it is different from the disk-like shape in mcx/mmc |
39 | | -cfg.detpos = [35, 25, max(cfg.node(:,3))]; %# redbird detpos can have multiple rows |
40 | | -cfg.detdir = [0, 0, -1]; %# redbird automatically computes the adjoint solutions, treating detectors as sources |
| 38 | +%#Redbird detector positions are point-like, directly sampling the output fluence; it is different from the disk-like shape in mcx/mmc |
| 39 | +cfg.detpos = [35, 25, max(cfg.node(:, 3))]; %#redbird detpos can have multiple rows |
| 40 | +cfg.detdir = [0, 0, -1]; %#redbird automatically computes the adjoint solutions, treating detectors as sources |
41 | 41 |
|
42 | | -%# redbird cfg does not need cfg.{nphoton,tstart,tend,tstep,elemprop}, which are required in mmclab |
| 42 | +%#redbird cfg does not need cfg.{nphoton,tstart,tend,tstep,elemprop}, which are required in mmclab |
43 | 43 |
|
44 | | -disp('== preprocessing domain ...') |
45 | | -%# calling rbmeshprep() populates other needed mesh data, such as cfg.{evol,nvol,face,area,reff,deldotdel,cols,idxsum} |
46 | | -%# this is similar to `cfg = mmclab(cfg, 'prep')` |
| 44 | +disp('== preprocessing domain ...'); |
| 45 | +%#calling rbmeshprep() populates other needed mesh data, such as cfg.{evol,nvol,face,area,reff,deldotdel,cols,idxsum} |
| 46 | +%#this is similar to `cfg = mmclab(cfg, 'prep')` |
47 | 47 |
|
48 | | -tic |
| 48 | +tic; |
49 | 49 | cfg = rbmeshprep(cfg); |
50 | | -toc |
| 50 | +toc; |
51 | 51 |
|
52 | | -disp('== run simulation ...') |
53 | | -%# Run forward simulation (you can also call rbrun()) |
54 | | -tic |
| 52 | +disp('== run simulation ...'); |
| 53 | +%#Run forward simulation (you can also call rbrun()) |
| 54 | +tic; |
55 | 55 | [detphi, phi] = rbrunforward(cfg); |
56 | | -toc |
| 56 | +toc; |
57 | 57 |
|
58 | | -%# rbrunforward returned two outputs: |
59 | | -%# detphi - the sampled measurements (fluence, not diffuse reflectance!) at all src/det pairs |
60 | | -%# phi - the fluence at all nodes for all sources and detectors (column-dimension) |
| 58 | +%#rbrunforward returned two outputs: |
| 59 | +%#detphi - the sampled measurements (fluence, not diffuse reflectance!) at all src/det pairs |
| 60 | +%#phi - the fluence at all nodes for all sources and detectors (column-dimension) |
61 | 61 |
|
62 | | -disp('== plot results ...') |
| 62 | +disp('== plot results ...'); |
63 | 63 | figure; |
64 | | -colormap('jet') |
65 | | -plotmesh([cfg.node log10(phi(:,1))], cfg.elem, 'y > 25') %# first column is from the src |
66 | | -%shading interp |
67 | | -title('forward solution from source') |
68 | | - |
69 | | -figure |
70 | | -colormap('jet') |
71 | | -plotmesh([cfg.node log10(phi(:,2))], cfg.elem, 'y > 25') %# 2nd column is from the detector |
72 | | -%shading interp |
73 | | -title('forward solution from detector') |
| 64 | +colormap('jet'); |
| 65 | +plotmesh([cfg.node log10(phi(:, 1))], cfg.elem, 'y > 25'); %#first column is from the src |
| 66 | +% shading interp |
| 67 | +title('forward solution from source'); |
| 68 | + |
| 69 | +figure; |
| 70 | +colormap('jet'); |
| 71 | +plotmesh([cfg.node log10(phi(:, 2))], cfg.elem, 'y > 25'); % # 2nd column is from the detector |
| 72 | +% shading interp |
| 73 | +title('forward solution from detector'); |
0 commit comments