Unexpected photon trajectories for very low scattering coefficients

[dkalsan]

Issue

Large number of scattering events take place in a medium with a very low scattering coefficient. To me this behaviour seems counterintuitive as I would expect little to no scattering events.

Setup

A 32x32 pixel 2D volume containing one circular artery in the middle surrounded by background tissue is being simulated. A pencil illumination source located at the top-middle of the volume is used and directed towards the artery. The background tissue is assigned a very low scattering coefficient $\mu_s=10^{-11}$ which leads to unexpected photon trajectories.

Observed behaviour

When plotting the photon trajectories for $\mu_s \in [10^{-11}, 10^0]$ the following can be observed:

• There is a large number of scattering events happening in the BG tissue for $\mu_s \in [10^{-11}, 10^{-8}]$
• Photon trajectories are traced outside of the volume (black = background tissue)
• Photons are terminated in the middle of the volume with weights much higher than the termination threshold (e.g. 0.8 weight, 0.0 threshold)
• The behaviour stabilizes at approximately $\mu_s \geq 10^{-6}$

Expected behaviour

There should be little to no scattering events in the background tissue for very low $\mu_s$, i.e. photon traces should be straight lines from the source origin (top-middle) to the artery, and then they should be scattered by the latter.

Code to reproduce

import pmcx
import numpy as np


def _circle(h, w, center, radius):
    Y, X = np.ogrid[:h, :w]
    dist_from_center = np.sqrt((X - center[0])**2 + (Y-center[1])**2)
    mask = dist_from_center <= radius
    return mask


optical_properties= [
    [0, 0, 1, 1],
    [0.02, 1e-11, 0.9, 1],
    [0.30, 80.25, 0.9, 1]
]

volume_height=32
volume_width=32
artery_center=(16, 16)
artery_radius=7
source_position=[0, 0, 16]
source_direction=[0, 1, 0]

volume = np.ones([1, volume_height, volume_width], dtype='uint8')
volume[:, _circle(volume_height,
                  volume_width,
                  artery_center,
                  artery_radius)] = 2

out = pmcx.run(debuglevel="M",
               seed=42,
               nphoton=10000,
               tstart=0,
               tend=5e-9,
               tstep=5e-9,
               outputtype="energy",
               issrcfrom0=1,
               unitinmm=0.1,
               vol=volume,
               srcpos=source_position,
               srcdir=source_direction,
               prop=optical_properties)

[fangq]

thanks, I am able to reproduce this issue in both python and mcxlab. I agree with you both the scattering behavior and the out-of-bbx trajectories are incorrect. I will look into this and will post my updates here.

mcxlab test code is attached

%%octave

cfg.nphoton=1e4;
cfg.seed=42;

cfg.vol=permute(uint8(ones(32,32)), [3,1,2]); % from 2d to 3d
cfg.vol(1,10:22, 10:22)=2;
cfg.issrcfrom0=1;
cfg.srctype='pencil';

cfg.srcpos=[0,0,16];   % src position must be located in the 2D plane
cfg.srcdir=[0 1 0];    % src dir must align in the plan (y-z in this case)

cfg.gpuid=1;
cfg.autopilot=1;

cfg.prop=[0, 0, 1, 1,
    0.02, 1e-11, 0.9, 1,
    0.30, 80.25, 0.9, 1];
cfg.unitinmm=0.1;

cfg.tstart=0;
cfg.tend=5e-9;
cfg.tstep=5e-9;
cfg.outputtype='energy';

[flux,detp,vol,seeds,traj]=mcxlab(cfg);
sortedlines=mcxplotphotons(traj);

[fangq]

@dkalsan, again, thanks for reporting this issue.

I looked into it further, and was able to identify the cause - the handling of near-zero mus causes issues when combined with the new DDA ray-tracing approach implemented last year. The photon voxel indices (flipdir(1:3)) no longer follow the floating-point positions (p.{x,y,z}). this was not a problem before using the DDA based photon advancing.

after a one-line change, I believe this behavior has been fixed. see my above commit.

the above matlab demo code now returns correct trajectories, see below.

feel free to reopen the ticket if you see anything else is not fixed.

[fangq]

Interestingly, mcxcl/mcxlabcl does not have this behavior - it uses native_divide for the affected line, which appears to be able to handle divided by 0.

[dkalsan]

Perfect, I will give it a test. Thanks a lot for taking the time for this!