Implicit mesh based MC

[asassa01]

Hi Qianqian,

I am using mmclab and I modified the demo "demo_immc_basic" (only the node-based section). I have two questions:

1) If I use a different source:

% cfg0.srctype='disk';
% cfg0.srcpos = [0.5 0.5 1];
% cfg0.srcdir=[0 0 -1];
% cfg0.srcparam1=[1 0 0 0];

MATLAB will crush. 

2) if I decide to output the field of the photon paths 

[flux_nimmc, MC_det] = mmclab(cfg);

it seems that no path is travelled in "region 2", which I assumed it was the interior of the sphere. Two more details about the code. The optical properties:

cfg0.prop = [0 0 1 1; 0.01 0 1 1; 0.01 0 1 1.5];

and the detector features:

cfg.detpos=[0.5 0.5 1 1];   

Thank you.

Angelo

[asassa01]

I realized that I have to add a few more details:

For question 2) the source features are:

cfg0.srcpos = [0.5 0.21 1]; 
cfg0.srcdir = [0 0 -1];

I also considered:

nbox = [nbox; [0.5 0.5 0.5]];  % insert new nodes (node 9)
rr=0.3; %(mm)
cfg.noderoi(9) = rr;

Thanks!
Angelo

--
You received this message because you are subscribed to the Google Groups "mmc-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to mmc-users+...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/mmc-users/232a2a48-a835-4ea4-9261-921f9d1b4c8fn%40googlegroups.com.

[Qianqian Fang]

hi Angelo,

using spatially distributed sources (like disk, planar, pattern), either internal or outside of the meshed object requires an extra preprocessing step - mesh retessellation. the details on this approach is described in

https://opg.optica.org/boe/viewmedia.cfm?uri=boe-7-1-171&html=true

briefly, mmc requires the knowledge of the index of the tetrahedral element (cfg.e0) enclosing the photon's position at any given time. Imagine you have a internal disk source spanning across a few tetrahedral elements, the launch position of the photon is dynamically computed, if you set cfg.e0 using the srcpos position, but this element is incorrect if the launch position moves to a different element. This gets harder when the disk is outside of the mesh (i.e. widefield), because there is no element enclosing the launch position.

in the case of a widefield source, you must call mmcaddsrc() and expand the mesh to include the entire launch aperture of the source, and properly label the element (label=-1) that enclose the source. see examples

https://github.com/fangq/mmc/blob/master/mmclab/example/demo_sfdi_2layer.m

is your disk source outside of the original mesh? or inside?

Qianqian

To view this discussion visit https://groups.google.com/d/msgid/mmc-users/CAJNg7oyCtdEu_s6OZsatE6SLEkF4f51nn0uNy%3DonHOa2kAxj-w%40mail.gmail.com.

[asassa01]

Thank you so much for the quick reply.  I will look into the link you sent me  and call mmcaddsrc() following the example.
 To answer your question, it is easier if I send you the code. As a start I would be happy also to use the pencil beam (which works) but I would like to use a detector at the boundary as in the code below.

It is not clear to me why the field of partial path lengths in region 2 is "0".

Thanks in advance,

Angelo

addpath('..\..\iso2mesh-master\')
addpath('../')

cfg0.nphoton = 1e6;

%%% PENCIL BEAM SOURCE%%%%%%%%%%%%%%%%
cfg0.srcpos = [0.5 0.21 1]; %[0.5 0.5 1]

cfg0.srcdir = [0 0 -1];

%%%%% DISK SOURCE %%%%%%%%%%%%%%%%%%%%%%%%%
% cfg0.srctype='disk';
% cfg0.srcpos = [0.5 0.5 1];% MATLAB CRASHES

% cfg0.srcdir=[0 0 -1];
% cfg0.srcparam1=[1 0 0 0];

cfg0.prop = [0 0 1 1; 0.01 0 1 1; 0.01 0 1 1.5];

cfg0.tstart = 0;
cfg0.tend = 5e-9;
cfg0.tstep = 5e-9;
cfg0.debuglevel = 'TP';
cfg0.method = 'grid';
cfg0.steps = [0.001 0.001 0.001];
dx=cfg0.steps(1,1);
cfg0.isreflect = 1;
cfg0.gpuid = -1;

% (a) generate bounding box and insert edge
Lx=1; %(mm)
Ly=1;
Lz=1;
nx=Lx/dx;
ny=Ly/dx;
nz=cfg0.srcpos(1,3)/dx;
[nbox, ebox] = meshgrid6(0:Lx, 0:Ly, 0:Lz);
fbox = volface(ebox);

nbox = [nbox; [0.5 0.5 0.5]];  % insert new nodes (node 9)

cfg = cfg0;

% (b) generate mesh
[cfg.node, cfg.elem] = s2m(nbox, num2cell(fbox, 2), 1, 100, 'tetgen1.5', [], [], '-YY');
cfg.elemprop = ones(size(cfg.elem, 1), 1);

% (c) label the edge that has node 9  and add radii
cfg.noderoi = zeros(size(cfg.node, 1), 1);

rr=0.3; %(mm)
cfg.noderoi(9) = rr;

 cfg.detpos=[0.5 0.5 1 1];   %  an N by 4 array, each row specifying a detector: [x,y,z,radius]

% run node-based iMMC
[flux_nimmc, MC_det] = mmclab(cfg);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%% PLOTTING %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CIRCLE
CC=[Lx/2 Ly/2];
th = linspace(0, 2*pi, 100);
rad=rr;
xc = rad * cos(th) + CC(1,1);
yc = rad * sin(th) + CC(1,2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% yz SECTION OF FLUENCE
pl_mat=squeeze(flux_nimmc.data(nx/2+1, 1:ny, 1:nz));
figure(1)
xax=0:dx:(nx-1)*dx;
zax=xax;
cols=turbo(100);
cols(1, :)=0;
cols(end, :)=1;
pcolor(xax,zax,rot90(pl_mat));
set(gca, 'ColorScale', 'log');
set(gca, 'Colormap', cols);
colorbar;
clim([1, 1e4]);
axis equal;
shading flat;
xlabel('y (mm)')
ylabel('z (mm)')

colorbar
hold on
plot(xc, yc, 'r-', 'LineWidth', 2);

To view this discussion visit https://groups.google.com/d/msgid/mmc-users/eb2c0a9c-9a39-43f0-8b72-b1ef55529fd3%40neu.edu.

[Qianqian Fang]

hi Angelo,

you were correct. immc currently failed to track the partial-path lengths when photons enter the ROI (node/edge/face). instead, the label of the element that carries the ROI is used.

I just created a new issue on github

https://github.com/fangq/mmc/issues/116

and was able to fix it using the below patch

https://github.com/fangq/mmc/commit/7be94c0f48ba382b2fbbc3a7219ac4ac0db810ef

please redownload the automatically built github-action packages from

https://mcx.space/nightly/github/

I also found that you set cfg.steps to 0.001 mm spacing - making the output grid 1001^3 voxels. On my laptop with 16GB memory, this caused a crash because the memory needs for this output exceeds my memory (8GB per output, 2x to copy to matlab). I suggest you use a bigger steps size, say 0.01. having exceedingly refined output volume while photon number is small won't actually gain better accuracy.

Qianqian

To view this discussion visit https://groups.google.com/d/msgid/mmc-users/CAJNg7ozx2xcb7rin_tej6Rhk%2BKetL4zRAi1LucjXuEAj2NDQkQ%40mail.gmail.com.