Couple suggested fixes to replay modes

[HirviP]

Dear Dr. Fang and MCX team,

We observed a couple of issues with the replay modes. We are not sure how these should be fixed, but we are happy to help where we can.

All test scripts can be found at scripts. We used the latest release (v2025.10 Kilo-Kelvin) of MCXLAB. We consider continuous-wave (CW) and frequency-domain (FD) modes, thus all scripts use only one long time window.

1) cfg.replaydet = -1 for replaying all detectors and saving outputs in separate volumes (output has 5 dimensions) does not work for the replay types with scattering counts, i.e., 'nscat', ''wp', 'wptof' and 'rfmus'. This can be observed with the test script replaydet_test.m with two detectors. As shown below for 'wptof', it seems that the results for both detectors are summed together along the first output in the fifth dimension, and the second output is just zeros. This resembles cfg.replaydet = 0.

2) For a pattern/pattern3d source, the mua derivative for the sine-weighted part (dYdmua = rfjac.data(:, :, :, :, :, 2)) in the FD ('rf') mode seems to be wrong on line 2272 in mcx/src/mcx_core.cu Line2272

The weight (first term) should be (-1) times the final detected weight of the photon packet ('-replayweight'); not the current weight ('p.w'), similarly as on Line2212 and Line2255.

We are unsure how to fix this for the pattern source, since there is not a similar line for the cosine-weighted term (dXdmua = rfjac.data(:, :, :, :, :, 1)) .

The effect of this bug can be observed with the script pattern_test.m. The Jacobians for CW intensity and FD amplitude should look very similar. If you change the source type to gaussian or disk, the Jacobians match, but not for the pattern source, as shown below. The Jacobian for X looks as expected.

Furthermore, USE DOUBLE seems to compute something different, but we guess it is not in use?

3) If the voxel grid unit is not 1 mm (cfg.unitinmm), are the path lengths correct when computing the detected weights with mcx_replayprep() in the standalone mode in mcx/src/mcx_utils.c Line3718? In the function mcx_replayinit() above for MATLAB/Python, 'plen' is scaled after computing the weight on Line3672, since 'mua' is already scaled to account for cfg.unitinmm. This was fixed for the embedded mode in orig_fix

4) Loosing one photon in replay has been discussed now and then in mcx-users. We observed that with cfg.voidtime = 0, the lost photon has total time-of-flight (TOF) that slightly exceeds the maximum simulation time (cfg.tend). As a result, the photon is excluded from the replay simulation in the original time window [cfg.start, cfg.tend]. In most test cases, we observed that all photons can be re-detected by extending cfg.tend for replay, but in some cases the detection still fails. Since this only affects one(/couple) photons with very long TOF, the numerical effect is presumably negligible, but we decided to still share this observation. Edit: This has also been discussed previously in mcx-users

The above is based on the script voidtime_TOFlimit_replay_test.m. Since photon loss occurs randomly, one might need to iterate this code multiple times. Below is the relevant output from one of our simulations with cfg.tend = 1e-8.

Many thanks in advance for any help!
Best regards,
Pauliina Hirvi and Jaakko Olkkonen