Fixed Pattern Noise due to xorshift128 RNG

[albertoQUEL]

We are observing fixed spatial pattern noise in simulations when using the disk source type. Although some starting seed numbers can minimize the observed pattern noise, it is generally pervasive across different seed numbers.

The current RNG implementation appears to use an xorshift128-style RNG to initialize each photon’s RNG state by combining the global seed words with the photon index

ran.reseed(seeds.x ^ i, seeds.y | i, seeds.z ^ i, seeds.w | i);

This implementation of the xorshift128 appears to be causing the fixed pattern noise. Note that pattern/noise generated by xorshift128 has been previously reported.[1]

I was able to eliminate the fixed pattern noise by preconditioning the seed and photon index with SplitMix64 before initializing the xorshift-style RNG state.[2] This approach preserves deterministic behavior for a given seed and photon index while causing adjacent photon indices to produce substantially different initial RNG states.

The attached figure shows top-down reflectance together with two-dimensional sagittal and coronal deposited-energy slices for a 20 × 20 × 20 mm block containing an embedded vessel. It compares the fixed pattern noise observed with the current implementation against the results obtained using SplitMix64-based reseeding.

Example .json that should reproduce the fix-pattern noise in the figure:
channel_block.json

Preliminary timing tests did not show a significant increase in simulation time with this implementation. Below is the code I used to implement the proposed change. It is worth noting that I used OpenAI Codex to assist with writing the code.
Replace the current per-photon reseeding line with:

MCX_reseed_photon(ran, seeds, i);

where

inline uint64_t MCX_splitmix64(uint64_t x) {
x += 0x9E3779B97F4A7C15ULL;
x = (x ^ (x >> 30)) * 0xBF58476D1CE4E5B9ULL;
x = (x ^ (x >> 27)) * 0x94D049BB133111EBULL;
return x ^ (x >> 31);
}

inline void MCX_reseed_photon(MCX_rand& ran, const dim4& seeds, uint64_t photonid) {
uint64_t base0 = ((uint64_t)seeds.x << 32) | seeds.y;
uint64_t base1 = ((uint64_t)seeds.z << 32) | seeds.w;
uint64_t state0 = MCX_splitmix64(base0 ^ photonid);
uint64_t state1 = MCX_splitmix64(base1 ^ (photonid + 0xD1B54A32D192ED03ULL));

if ((state0 | state1) == 0ULL) {
state1 = 0x9E3779B97F4A7C15ULL;
}

ran.reseed((uint32_t)(state0 >> 32), (uint32_t)state0,
(uint32_t)(state1 >> 32), (uint32_t)state1);
}

Citations:

1. Haramoto, H., Matsumoto, M., and Saito, M. “Unveiling patterns in xorshift128+ pseudorandom number generators.” Journal of Computational and Applied Mathematics, 2022.
2. Steele, G. L. Jr., Lea, D., and Flood, C. H. “Fast Splittable Pseudorandom Number Generators.” OOPSLA, 2014.

[fangq]

can you share your script for producing the patterned output?

this RNG was used by default in mcx/mcxcl/mmc, and I haven't seen such pattern. I would like reproduce first.

[albertoQUEL]

I had attached the .json file that was used to run the simulation. Would this be enough to try to reproduce the output?

I have written a Python “wrapper” for internal use that produces the .json files and runs them on the unmodified uMCX code. However, the Python script itself would not technically run without our internal wrapper functions.

Below is the re-attached .json file and its contents. Note that I use a z-layer of high absorption to serve as a “sensor” for simulating a reflectance image. This seems to allow for better CPU/memory usage on my laptop compared to running a true detector-based reflectance simulation.

channel_block.json

{
  "Session": {
    "ID": "channel_block",
    "Photons": 50000000,
    "RNGSeed": 1648335518,
    "DoMismatch": true,
    "DoSaveVolume": true,
    "DoNormalize": true,
    "DoPartialPath": false,
    "OutputType": "e",
    "MaxDetPhoton": 50000000,
    "SaveDetFlag": 1,
    "BlockSize": 64,
    "DeviceID": 1
  },
  "Forward": {
    "T0": 0.0,
    "T1": 5e-09,
    "Dt": 5e-09
  },
  "Domain": {
    "Dim": [
      200,
      200,
      201
    ],
    "LengthUnit": 0.1,
    "MediaFormat": "byte",
    "Media": [
      {
        "mua": 0.0,
        "mus": 0.0,
        "g": 1.0,
        "n": 1.0
      },
      {
        "mua": 0.01,
        "mus": 10.000000000000002,
        "g": 0.9,
        "n": 1.37
      },
      {
        "mua": 0.05,
        "mus": 10.000000000000002,
        "g": 0.9,
        "n": 1.37
      },
      {
        "mua": 1000000000.0,
        "mus": 0.0,
        "g": 1.0,
        "n": 1.0
      }
    ]
  },
  "Optode": {
    "Source": {
      "Type": "disk",
      "Pos": [
        100.0,
        100.0,
        1.0
      ],
      "Dir": [
        0,
        0,
        1
      ],
      "Param1": [
        200.0,
        0.0,
        0,
        0
      ]
    }
  },
  "Shapes": [
    {
      "Box": {
        "Tag": 1,
        "O": [
          0.0,
          0.0,
          1.0
        ],
        "Size": [
          200.0,
          200.0,
          200.0
        ]
      }
    },
    {
      "Cylinder": {
        "Tag": 2,
        "C0": [
          0.0,
          100.0,
          21.0
        ],
        "C1": [
          200.0,
          100.0,
          21.0
        ],
        "R": 10.0
      }
    },
    {
      "ZLayers": [
        [
          0.0,
          1.0,
          3
        ]
      ]
    }
  ]
}

[fangq]

@albertoQUEL, my apologies, I did not read your initial post carefully, where you have already included the .json file.

I ran your input file using my locally compiled umcx, unfortunately I did not see the pattern shown in your screenshot. see my screenshot

my umcx binary was compiled using "make nvc" using nvc version 25.3-0; the test machine is a laptop with NVIDIA 4070 mobile GPU (driver version 580.95.05) and running on Ubuntu 22.04.

I am curious how you compiled umcx? what hardware, driver and OS version did you run the above test?

[fangq]

from reading the Haramoto paper, it appears that the pattern only appears if the sequence of random numbers are consumed in a particular way (triplets) and the scaling of the domain to an extremely small level - neither of such happens in this particular MC code.

but I am still curious if the patterned simulation you showed above was a result of umcx, or LLM mock output. If it is latter, I would close this ticket until real reproducible artifacts manifest.

[albertoQUEL]

I wonder if the fixed pattern noise is coming from my use of the CPU (Windows OS) rather than GPU (was running on my laptop with no GPU at the time). Let me gather some more information on the code compilation and if I can also replicate your results (no fixed pattern) with a GPU implementation.

I should be able to get back to you this week!