Incorrect exit position of detected photons

[devhello145]

Some detected photons (many) have incorrect exit position. In this example I created a sphere of R=1, and detector with R=2 (centers of them are [0,0,0]). At the end of simulation, position of detected photons are inside sphere of R=1 and even equal to [0, 0, 0].

Code to reproduce (Octave, run on Linux):

%% add path to iso2mesh and mmclab
%
%

%% Create a sphere with R=1
[node, facet, elem] = meshasphere([0, 0, 0], 1, 0.1);
% plotmesh(node, elem);
elem(:, 5) = 1;

cfg = struct();
cfg.seed = 1;

cfg.srctype = 'isotropic';
cfg.srcpos = [0.0001, 0., 0.];  % bias x coord due to error
cfg.srcdir = [0 0 1];

%% set time of flight  (set to maximum) 
cfg.tstart = 0;
cfg.tend = 1e-6;
cfg.tstep = 1e-6;

cfg.isreflect = 1;
cfg.issaveexit = 1;
cfg.issaveseed = 1;  % overrides in mmclab by number of outputs

cfg.outputtype = 'fluence';

%% pernament
cfg.unitinmm = 1.;

%% tuned
cfg.minenergy = 1e-9;
cfg.maxdetphoton = 1e+7; % 1e+8 - broken
cfg.method = 'elem';
cfg.gpuid = -1;  % -1 - CPU mode, 1 - GPU mode

%% additional
cfg.debuglevel = 'TP';  % M - movement info;

cfg.prop = [0, 0, 1, 1; 1e-4, 1e-4, 1, 1];
cfg.detpos = [0 0 0 2]; % x y z radius

%% pernament (mesh)
cfg.node = node;
cfg.elem = elem(:, 1:4);
cfg.elemprop = elem(:, 5);

%% tuned
cfg.nphoton = 1e5;
cfg.gpuid = -1;

[fluence, detp, ncfg, seeds] = mmclab(cfg);
x = detp.p(:, 1);
y = detp.p(:, 2);
z = detp.p(:, 3);

r = sqrt(x.*x + y.*y + z.*z);

fprintf('Min detp.p (radius): %f\n', min(r));
fprintf('Max detp.p (radius): %f\n', max(r));

fprintf('Min detp.ppath: %f\n', min(detp.ppath));
fprintf('Max detp.ppath: %f\n', max(detp.ppath));

Output:

```txt Surface Mesh Extraction Utility (Based on CGAL 3.8) (modified for iso2mesh by Qianqian Fang) http://iso2mesh.sf.net RNG seed 1648335518 set initial mesh size to 50 Refining... maximum node number is set to 40000 done. Final number of points: 793 generating tetrahedral mesh from closed surfaces ... creating volumetric mesh from a surface mesh ... volume mesh generation is complete Launching MMCLAB - Mesh-based Monte Carlo for MATLAB & GNU Octave ... Running simulations for configuration #1 ... mmc.seed=1; mmc.srctype='isotropic'; mmc.srcpos=[0.0001 0 0]; mmc.srcdir=[0 0 1 0]; mmc.tstart=0; mmc.tend=1e-06; mmc.tstep=1e-06; mmc.isreflect=1; mmc.issaveexit=1; mmc.issaveseed=1; mmc.outputtype='fluence'; mmc.unitinmm=1; mmc.minenergy=1e-09; mmc.maxdetphoton=1e+07; mmc.method='elem'; mmc.gpuid=-1; mmc.debuglevel='TP'; mmc.prop=1; mmc.detnum=1; mmc.nn=1762; mmc.elem=[8545,4]; mmc.ne=8545; mmc.nphoton=100000; mmc.facenb=[8545,4]; mmc.evol=8545; mmc.e0=4803; done 2 simulating ... ############################################################################### # Mesh-based Monte Carlo (MMC) - OpenCL # # Copyright (c) 2010-2020 Qianqian Fang # # http://mcx.space/#mmc # # # #Computational Optics & Translational Imaging (COTI) Lab [http://fanglab.org]# # Department of Bioengineering, Northeastern University, Boston, MA, USA # # # # Research funded by NIH/NIGMS grant R01-GM114365 # ############################################################################### $Rev:: $v2021.2$Date:: $ by $Author:: $ ############################################################################### seed=1 simulating ... Progress: [==============================================================] 100% done 60 speed ... 1666.67 photon/ms, 2619822 ray-tetrahedron tests (0 overhead, 43663.70 test/ms) detected 100000 photons source 1 total simulated energy: 100000.000000 absorbed: 0.00997% normalizor=1e-05 saving fluence ...saving detected photons ... done 64 done 66 Min detp.p (radius): 0.000000 <-- See HERE!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Max detp.p (radius): 0.999994 Min detp.ppath: 0.994178 Max detp.ppath: 1.000085 ```

[fangq]

@devhello145, I am not sure I understood what was the problem.

The detector sphere covers the entire surface of the domain, so, all photons are detected and are located on the mesh surface. This to me is expected.

[devhello145]

As I understand, photons are detected after leaving a mesh. In this example, all detected photons should be located on the surface of sphere (or r=sqrt(x^2+y^2+z^2) =1). In this example, many detected photons have r < 1 and even r=0. Why??? Last 4 lines of output:

Min detp.p (radius): 0.000000    <-- See HERE!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Max detp.p (radius): 0.999994
Min detp.ppath: 0.994178
Max detp.ppath: 1.000085

Gpu (opencl) version of code (gpuid=1) doesnot have this problem (all photons have r=1).

[devhello145]

Can you say when you plan to fix this bug?

[fangq]

I have already been working on it. all invalid positions are [0,0,0] (while ppath is nearly 1). inside the ray-tracer thread, no such position appears, but somehow they are set to 0 when copied to the host.

you can use their position to remove these photons or use OpenCL output before I find a fix.

[devhello145]

Thank you for your fast feedback. Yes you are right, incorrect position is only [0,0,0], other is good.

I more often use cpu version due to gpu version is less stable. // Note: When i run gpu version with large number photons, it's quickly run till 98% and stuck for 4 hours (or more). It's hard to find bug

[fangq]

@devhello145, sorry for taking a long time for me to finally fix this issue.

if you still have needs for using the CPU based simulations, please download the automatically compiled binaries from http://mcx.space/nightly/github/ once this github action script is complete (show green checks)

https://github.com/fangq/mmc/actions/runs/6044719776/job/16403710100