Detected photon discrepancy between pmcx and mcxlab

[rajsinghmn]

Hi,

I'm finding a pretty substantial discrepancy for detected photon counts and partial paths between MCXlab and pMCX and was wondering if someone might be able to help me pin down the source. Below are the scripts I'm using for comparison:

Python Validation Script:

import pmcx
import numpy as np

num_photons = 1e7
edge_length_in_voxels = 6
vol = np.ones([edge_length_in_voxels, edge_length_in_voxels, edge_length_in_voxels],dtype="uint8")

cfg = {
    "seed": 123456,           
    "nphoton": num_photons, 
    "vol": vol,          
    "prop":[[0,0,1,1],[0.000, 5.0, 0.9, 1.4]], 
    "srcpos": [edge_length_in_voxels/10,edge_length_in_voxels/10,0], 
    "srcdir":[0,0,1], 
    "detpos":[[edge_length_in_voxels/2,edge_length_in_voxels/2,0,edge_length_in_voxels]],
    "tstart":0, 
    "tend":5e-5, 
    "tstep":5e-5,
    "autopilot":1,
    "gpuid":1,
    "issavedet":1,                        
    "issrcfrom0":1,                       
    "isreflect": 1,
    "isrefint": 1,
    "maxdetphoton": 1e8,    
    "savedetflag" : "dpx",
    "unitinmm": 100                        
}

res=pmcx.run(cfg)
print(res.keys())

remittedDetections = res["detp"][:,res["detp"][-1]<0]
remittedDetections = remittedDetections.T
print(f"{remittedDetections.shape[0]} remitted photons detected")
mua=0.001
dref = np.sum(np.exp(-remittedDetections[:,1]*mua*cfg["unitinmm"]))/num_photons
print(f"Diffuse Reflectance: {dref:.7f}")

Python output:

7342212 remitted photons detected
Diffuse Reflectance: 0.6375916

Matlab Validation Script:

clear;

cfg.seed=123456;
Np = 1e7;
edge = 6;

cfg.nphoton=Np;
cfg.vol=ones(edge,edge,edge);
cfg.vol=uint8(cfg.vol);
cfg.prop=[0 0 1 1           
   0.001 5.0   0.9 1.4];

cfg.srcpos=[edge/2, edge/2,0];
cfg.srcdir=[0 0 1];
cfg.detpos=[edge/2,edge/2,0,edge];

cfg.tstart=0;
cfg.tend=5e-5;
cfg.tstep=5e-5;
cfg.autopilot=1;
cfg.gpuid=1;
cfg.issavedet = 1;
cfg.issrcfrom0=1;
cfg.isreflect=1;
cfg.isrefint=1;
cfg.maxdetphoton = 1e8;
cfg.savedetflag ='dpx';
cfg.unitinmm = 100;

[f1ux,detection]=mcxlab(cfg);

%%
remittedIdxs = detection.data(5,:)<0;
mua=0.001;
dref = sum(exp(-detection.ppath(remittedIdxs)*mua*cfg.unitinmm))/Np;
display( string(sum(remittedIdxs)) + ' remitted photons detected')
display('Diffuse Reflectance: ' + string(dref))

Matlab output:

9768689 remitted photons detected
Diffuse Reflectance: 0.80187

[fangq]

@rajsinghmn, have you tried the latest nightly build mcxlab? what was the version of your pmcx? (use pmcx.__version__)

from my test, the two solvers have similar results. see below screenshot

also, I noticed that the two simulations have some differences, for example,

• cfg.prop(2,1) is 0.001 while that in pmcx is 0.0,
• cfg.srcpos in mcxlab is [3 3 0] while that in pmcx is [0.6, 0.6, 0].

if I set pmcx's srcpos to match that in mcxlab, I am able to get a closer match.

I am closing this issue because it seems no longer a problem in the latest code. Feel free to reopen if you still have difficulty to match.

[rajsinghmn]

Thanks for checking @fangq . I believe this narrows down the issue to a difference between the latest builds(both pMCX and MCXLab) and v2020 of MCXLab. The previous output was using v2020 MCXLab; after downloading the latest v2023 MCXLab build from here I now get something closer to 0.6236 for calculated reflectance.

The MCXLab V2020 result is closer to the diffusion equation ( 0.7932 ) and MCML for the same optical property set :

Any thoughts on what might explain the difference between v2020 and v2023? Would you recommend using an earlier version of pMCX (v0.0.13 used above) to obtain a result consistent with v2020?

[fangq]

@rajsinghmn, both v2020 and the latest github version (to be released as v2023) have a built-in validation script against diffusion based diffuse reflectance, see

https://github.com/fangq/mcx/blob/master/mcxlab/examples/demo_diffuse_reflectance_validation.m

I ran both v2020 and v2023 with this script and both got matched results against analytical models. I suggest you check your diffuse reflectance calculations with the formula used in our test script and see if there is any difference.

regarding the differences in detected photon number you noticed between v2020 and v2023, I believe this is because the way how you selected the "reflected/remitted" photons. In your code, I see you used the z-coordinate of the exit position to find those reflected from z=0. this is not robust and is influenced by floating-point precision of the ouput.

the better way to do this is to set cfg.bc='______001000'. Please see my previous post and mcxlab help on this flag.

by removing cfg.detpos and adding cfg.bc, I was able to get matched detected photons at z=0 plane, see

>> [f1ux,detection]=mcxlab(cfg);
Launching MCXLAB - Monte Carlo eXtreme for MATLAB & GNU Octave ...
Running simulations for configuration #1 ...
mcx.seed=123456;
mcx.nphoton=1e+07;
mcx.dim=[6 6 6];
mcx.mediabyte=1;
mcx.medianum=2;
mcx.srcpos=[3 3 0];
mcx.srcdir=[0 0 1 0];
mcx.tstart=0;
mcx.tend=5e-05;
mcx.tstep=5e-05;
mcx.autopilot=1;
mcx.gpuid=1;
WARNING: redundant field 'issavedet'
mcx.issrcfrom0=1;
mcx.isreflect=1;
mcx.isrefint=1;
mcx.maxdetphoton=1e+08;
mcx.savedetflag=21;
mcx.unitinmm=100;
mcx.bc='______001000';
###############################################################################
#                      Monte Carlo eXtreme (MCX) -- CUDA                      #
#          Copyright (c) 2009-2020 Qianqian Fang <q.fang at neu.edu>          #
#                             http://mcx.space/                               #
...
###############################################################################
$Rev::af3d38$ v2020 $Date::2023-08-12 18:01:11 -04$ by $Author::Qianqian Fang $
###############################################################################
- variant name: [Fermi] compiled by nvcc [11.5] with CUDA [11050]
- compiled with: RNG [xorshift128+] with Seed Length [4]

GPU=1 (NVIDIA GeForce RTX 2060) threadph=162 extra=46720 np=10000000 nthread=61440 maxgate=1 repetition=1
initializing streams ...    init complete : 0 ms
requesting 1280 bytes of shared memory
launching MCX simulation for time window [0.00e+00ns 5.00e+04ns] ...
simulation run# 1 ... 
kernel complete:    12865 ms
retrieving fields ...   detected 9768179 photons, total: 9768179    transfer complete:  13513 ms
normalizing raw data ...    source 1, normalization factor alpha=0.000000
data normalization complete : 13514 ms
simulated 10000000 photons (10000000) with 61440 threads (repeat x1)
MCX simulation speed: 778.33 photon/ms
total simulated energy: 10000000.00 absorbed: 19.80699%
(loss due to initial specular reflection is excluded in the total)
>> detection

detection = 

  struct with fields:

       detid: [9768179×1 int32]
       ppath: [9768179×1 single]
           p: [9768179×3 single]
        prop: [2×4 double]
    unitinmm: 100
        data: [5×9768179 single]

comparing that with the output from v2023

>> cfg=rmfield(cfg,'detpos');
>> cfg.bc='______001000';
>> [f1ux,detection]=mcxlab(cfg);
Launching MCXLAB - Monte Carlo eXtreme for MATLAB & GNU Octave ...
Running simulations for configuration #1 ...
mcx.seed=123456;
mcx.nphoton=1e+07;
mcx.dim=[6 6 6];
mcx.mediabyte=1;
mcx.medianum=2;
mcx.srcpos=[3 3 0 1];
mcx.srcdir=[0 0 1 0];
mcx.tstart=0;
mcx.tend=5e-05;
mcx.tstep=5e-05;
mcx.autopilot=1;
mcx.gpuid=1;
mcx.issavedet=1;
mcx.issrcfrom0=1;
mcx.isreflect=1;
mcx.isrefint=1;
mcx.maxdetphoton=1e+08;
mcx.savedetflag=21;
mcx.unitinmm=100;
mcx.bc='______001000';
###############################################################################
#                      Monte Carlo eXtreme (MCX) -- CUDA                      #
#          Copyright (c) 2009-2023 Qianqian Fang <q.fang at neu.edu>          #
#                             http://mcx.space/                               #
...
###############################################################################
$Rev::af3d38$ v2023  $Date::2023-08-12 18:01:11 -04$ by $Author::Qianqian Fang$
###############################################################################
- variant name: [Fermi] compiled by nvcc [11.5] with CUDA [11050]
- compiled with: RNG [xorshift128+] with Seed Length [4]

GPU=1 (NVIDIA GeForce RTX 2060) threadph=162 extra=46720 np=10000000 nthread=61440 maxgate=1 repetition=1
initializing streams ...    init complete : 0 ms
requesting 1280 bytes of shared memory
launching MCX simulation for time window [0.00e+00ns 5.00e+04ns] ...
simulation run# 1 ... 
kernel complete:    12694 ms
retrieving fields ...   detected 9766735 photons, total: 9766735    transfer complete:  13335 ms
normalizing raw data ...    source 1, normalization factor alpha=0.000000
data normalization complete : 13336 ms
simulated 10000000 photons (10000000) with 61440 threads (repeat x1)
MCX simulation speed: 789.08 photon/ms
total simulated energy: 10000000.00 absorbed: 19.79982%
(loss due to initial specular reflection is excluded in the total)
>> detection

detection = 

  struct with fields:

       detid: [9766735×1 int32]
       ppath: [9766735×1 single]
           p: [9766735×3 single]
        prop: [2×4 double]
    unitinmm: 100
        data: [5×9766735 single]

again, the results from the two versions are largely matching. as to why it does not match your analytical model, I refer you to our validation script. please verify your diffusion calculation.

[rajsinghmn]

cfg.bc='__001000'

Thanks so much for your help @fangq ! The discrepancies I noticed was indeed due to my approach for collecting photons at z<=0. I'm now seeing matched results using the cfg.bc approach.