Monte Carlo eXtreme for OpenCL (MCX-CL)
MCX-CL is the OpenCL version of Monte Carlo eXtreme (MCX),
a GPU-accelerated Monte Carlo simulator of light transport in 3D turbid media such as
biological tissues. Where MCX needs an NVIDIA (CUDA) or AMD (HIP) GPU, MCX-CL runs on
almost any processor with an OpenCL driver: NVIDIA, AMD and Intel GPUs, integrated
graphics, Apple silicon, and multi-core CPUs. It reads the same input files, accepts
the same MATLAB/Python inputs and writes the same outputs as MCX, and for a given seed
it traces the same photon paths.
MCX-CL can be used in four ways:
| Mode | What it is | Best for |
|---|
mcxcl command line | a standalone executable reading a JSON input file | batch jobs, clusters, scripting, reproducible inputs |
| MCXLAB-CL (MATLAB/Octave) | mcxlabcl(cfg), a mex function with in-memory inputs/outputs | interactive analysis, image reconstruction, plotting |
| PMCXCL (Python) | pmcxcl.run(cfg), a Python module installed with pip | Python/NumPy workflows, notebooks, macOS |
| MCX Studio (GUI) | a graphical editor and viewer that runs mcx, mcxcl and mmc | learning MC, designing and inspecting simulations |
Table of Contents
What's New
MCX-CL v2026 (3.0) is a milestone release that brings MCX-CL to near feature
parity with CUDA-based MCX, and makes its ray-tracing results verifiable
step-by-step against MCX. Notable features ported to MCX-CL in this release
include
- directly building adjoint Jacobians via
cfg.outputtype='adjoint', including
the adjoint_dcoeff/mus/musp and dual Jacobian outputs, saving both forward
and Jacobian in the adjoint mode
- frequency-domain (FD/RF) forward simulations and adjoint-based FD-Jacobians
(complex-valued), and RF replay
- split-voxel MC (SVMC) for hybrid voxel/surface domains (MCX v2021.2)
- polarized light simulations (MCX v2023)
- mua/mus/g/n 4-float media format (MCX v2024.2)
- wedge/fan line sources and slit sources with Gaussian broadening
- international language support (
--lang)
- 1x1x1 volumes, and
.jdt/.jdb detected photon and trajectory files
with separate srcid records
- per-thread random number seeds derived exactly as in MCX, so that MCX and
MCX-CL trace the same photon paths for the same seed, on every OS
- exact output sums:
make fixed accumulates the output as 64-bit fixed-point
integers instead of a float with a second float buffer (MCX v2026); the sums
do not depend on the order of the additions, so a run with -A 0 is
bit-reproducible. Each sum is a pair of 32-bit words, low and high words in
separate planes, so only 32-bit atomics are needed and a deposit touches
4 bytes per element, as the float buffers do. The speed is about that of the
default build: 1.06x on an RTX 4090, 1.01x on an A100, 0.99x on a TITAN V and
0.97x on a CPU (POCL), from 0.89x to 1.44x over 13 built-in examples
In addition, this release contains the following fixes and improvements
- the long-standing OpenCL GPU memory leak on NVIDIA devices, first reported in 2020,
is finally fixed (#66, see also mmc#127): on NVIDIA's OpenCL runtime, every
simulation lost GPU memory, so repeated mcxlabcl/pmcxcl calls in one
MATLAB/Octave/Python session eventually ran out of GPU memory; leaked kernel-launch
events kept a reference to the OpenCL context, so the GPU buffers were never freed.
AMD, Intel and CPU (POCL) devices were not affected
- 1.2-1.5x faster simulations on OpenCL devices, from dynamic per-block photon
dispatch, 1-byte label volumes and
-cl-unsafe-math-optimizations; -A 0
keeps the static photon split and strict math for reproducible trajectories
- fixed an off-by-one photon index in wp/dcs/wptof/rfmus replay that used the
time-of-flight and weight of the next photon; wp/dcs/jacobian/wltof replay
outputs now agree with MCX to ~1e-3 (relative L2)
- ported the MCX replay normalization fixes (
-U 2, time-gate clamping,
SVMC mus division, ns time-of-flight accumulation, double-buffering)
- fixed all-zero RF replay (
-O r, -O s) output; RF and rfmus replay now
agree with MCX (relative L2 4e-5 and 2e-3)
- fixed adjoint Jacobian outputs (
-O A/D/U/V/W/Q) on the command line and a
segfault reading JSON detectors with a Dir field; adjoint detector sources
are now launched as disks so that detdir(4) shapes the detector beam
- pattern/pattern3d sources now work in wp/dcs/wptof/rfmus replay and carry
per-pattern weights in the RF imaginary parts; fixed a segfault with
multi-pattern pattern3d sources
- fixed forcing g=1 in continuous media (#65), a macOS OpenCL runtime error (#61)
and a double-free in
mcx_clearfluence
- zmat is upgraded to 1.0.0 (miniz 3.1.1, lz4 1.10.0)
- new step-by-step trajectory tests (
test/testtraj.sh) compare MCX-CL
against the MCX golden trajectories, and run, together with testmcx.sh,
on a CPU OpenCL device (POCL) in the CI with -A 0 -t <photons> -T 1
- polarized simulations now match MCX: the Mie scattering matrix is tabulated at
1000 angles as in MCX instead of 181 (the coarse table over-weighted forward
scattering), the source Stokes vector (
Optode.Source.IQUV) is read from JSON
input, and Stokes trajectories (--trajstokes 1) are sized and saved correctly
--maxvoidstep and --trajstokes are accepted on the command line as in MCX- the host-generated RNG seeds had ~75% of their bits set on Linux and macOS, which
biased the first photon of each thread; the per-thread seeds are now derived by the
xorshift128+ generator of the kernels instead of the C library
rand(), which differs
between Linux, macOS and Windows, so a fixed seed simulates the same photons on every
OS (the random streams for a given seed change, replay of saved seeds is unaffected)
We would like to thank @hanno-gmg (#60) and @ptrhbt (#67) for reporting issues.
The most important changes since v2025.10 are summarized below; ChangeLog.txt
lists all of them, with their commit hashes.
Features ported from MCX
- [7c7c6d5] split-voxel MC (SVMC); [00fcf09] split-voxel sides and exits tracked as in MCX
- [68f3211] polarized light simulations; [b3af30b] Mie matrix tabulated at 1000 angles as in MCX
- [8b05259] the mua/mus/g/n 4-float (
asgn_float) media format
- [5756a54, 7e92747] RF replay; [00bdfae, bcd8067] RF forward and RF/CW adjoint Jacobians, including
the dcoeff/mus/musp and dual outputs; [0f692be] forward fluence and Jacobians both returned
- [4d53ce3] wedge/fan line sources and slit sources with Gaussian broadening
- [5c2816e] international language support (
--lang)
- [e0b526a] 1x1x1 volumes
- [0b17a15] detected photons and trajectories saved as
.jdt/.jdb, with separate srcid records
- [95b9dd1] per-thread RNG seeds derived exactly as in MCX: the same seed gives the same photon paths
- [1d2af50]
make fixed: exact, order-independent 64-bit fixed-point output sums, as in MCX
- [2a14c42] pattern sources in wp/dcs/wptof/rfmus replay; [a2bc0a4] per-pattern RF imaginary parts
- [4fa20b7]
--maxvoidstep and --trajstokes; [c06cdf0] issave2pt, flog and lang in
mcxlabcl/pmcxcl, as in mcxlab/pmcx
Performance
- [8680d3f] dynamic per-block photon dispatch, 1-byte label volumes and the unsafe-math JIT
option in autopilot mode: 1.2-1.5x faster
- [2964165] parallel zlib/gzip compression of JData outputs, with a block index for parallel decoding
Bug fixes
- [7bb12a9] off-by-one photon index in wp/dcs/wptof/rfmus replay; replay now agrees with MCX
- [7e92747] all-zero RF replay (
-O r, -O s) output
- [57545b0] adjoint outputs (
-O A/D/U/V/W/Q) on the command line; detector JSON with Dir
- [d7ea87b]
-U 2 skips only the detected-weight normalization in replay; [b5d3b0a] double
preprocessing with -j
- [8c9c3f0] GPU memory leak on NVIDIA devices (#66)
- [0ad93bb, 48e38ac] host RNG seeds had ~75% of their bits set; now xorshift128+, uniform and
the same on every OS
- [5acc84f] photons entering label 0 end as in MCX, fixing a hang with interior label-0 regions
- [5e139f0]
cl_ext_float_atomics kernels built as OpenCL C 3.0, fixing Intel Arc GPUs
- [7b8f5f2] strict math is used when saving or replaying seeds
- [00d55fc]
-r repeats or splits the photons as documented; [25ea487] a time step of 0 is one gate
- [614e647] RF forward output keeps its imaginary part in
.jnii/.bnii files
- [04a7734]
Session.OutputType names read in full, not just their first letter
- [e14dbef, ef4ec6b, aca5a92]
--dumpjson, mcx2json and json2mcx round-trip all settings
- [68b84f8] per-voxel media scaled mua/mus by
unitinmm twice when cfg.unitinmm was set before
cfg.vol, and labelplus scaled mus/g instead of mua/mus
- [3ff79ff] mcxlabcl: an empty field of a cfg struct array keeps its default; it was read as 0
- [552d641] mcxlabcl accepts any numeric class (
gpuid=int32(1) failed), logical volumes and column
vectors; pmcxcl accepts debuglevel='T'
- [bd6078a]
cfg.flog also controls the echo of the cfg fields: 'none' silences everything
- [cbb5436]
cfg.flog='none' and -l none silence the log on every OS ('NUL' and '/dev/null' are accepted too; 'NUL' wrote a file on Linux, #280)
- [ee5dbda] the mex files are no longer UPX-packed: a packed mex aborted MATLAB/Octave whenever it raised an error
- [1346c53] the device can be picked by name:
-G 4070 or cfg.gpuid='4070' uses the first device whose name contains 4070 (ignoring case), -G titan:2 the 2nd one
- [6f2f2bf] on Apple Silicon GPUs, every simulation after the first in a MATLAB/Octave/Python session could start from the previous output; the output buffer is now cleared by a write
Testing and robustness
- [253b22f] step-by-step trajectory tests against MCX's golden data, run on a CPU (POCL) in CI
- [7fe8bea] host unit tests and short GPU feature tests; [c9de165] host coverage about 95%
- [55da8c8] four rounds of fuzzing of inputs and options: all crashes, hangs and memory errors
fixed [70ef956, ef650f2]; input values now validated
- [73a7aea]
test/fuzz/dfuzz.py and the built-in inputs of MCX in test/builtin; six more rounds
of fuzzing found five bugs shared with MCX, fixed in [7c1fcba]
- [27effdc, 162ebd4] unit tests of the mcxlabcl and pmcxcl interfaces: input types, every output and
source type, and short absorbed-fraction checks against references shared with MCX
Packaging
- [a3aefe2] Linux releases built in manylinux2014 (glibc 2.17+); Octave 10 mex files for macOS/Windows
- pmcxcl wheels for Linux aarch64 (#67); [ee9125f] macOS wheels tagged with their real
architecture (#60); [b429887]
pmcxcl.mcxlab runs with pmcxcl
- [65dfa31] bundled zmat updated to 1.2.0; [87dd73e] cJSON updated to 1.7.19
- [3caf4bc] macOS builds use the Khronos OpenCL headers (src/khronos), free of the
deprecation warnings of Apple's OpenCL framework headers
Introduction
Monte Carlo eXtreme simulates how light propagates in 3D heterogeneous, scattering
media by tracing a large number of photon packets through a voxelated domain. Each
packet is attenuated by absorption, randomly redirected by scattering, and
reflected or refracted at interfaces between media; accumulating their weights gives
a statistically exact solution of the radiative transfer equation (RTE), the model
behind diffuse optics, fluorescence imaging, photobiomodulation dosimetry and many
other biophotonics applications. Tracing the photons in parallel on a GPU makes this
hundreds to thousands of times faster than a single-threaded CPU code.
MCX-CL implements the MCX algorithm in OpenCL [Yu2018], so the same simulation runs on
GPUs and CPUs from any vendor. Since v2026, MCX-CL supports nearly all MCX features and
traces the same photon paths as MCX for the same seed; the differences are listed in
MCX-CL and MCX: differences.
What MCX-CL can do
- arbitrary 3D heterogeneous media given as label volumes, continuously varying
optical properties, or built from simple shapes; 2D simulations
- time-resolved, continuous-wave (CW) and frequency-domain (RF) simulations
- 18 source types, from pencil beams to wide-field patterns and spatial-frequency
sources; many sources in one run
- disk detectors and bounding-box faces as detectors; per-photon partial path
lengths, scattering counts, momentum transfer, exit positions/directions
- Fresnel reflection, total absorption, mirror and cyclic boundary conditions
- Jacobians by photon replay or by the adjoint method, in CW, time and frequency
domain
- photon sharing, polarized light, split-voxel MC (SVMC), custom phase functions
- several devices - even from different vendors - in one simulation
- JSON-based inputs and JData/NeuroJSON-based outputs
The MCX family
| Software | Description |
|---|
| MCX | voxel-based MC in CUDA (NVIDIA) or HIP (AMD); MCXLAB, PMCX and MCX Studio |
| MCX-CL | this package: the OpenCL version of MCX; MCXLAB-CL (mcxlabcl) and PMCXCL (pmcxcl) |
| MMC | mesh-based MC on tetrahedral meshes, on CPUs, OpenCL, CUDA and OptiX; MMCLAB and PMMC |
| MCX Cloud | run MCX and MMC in a web browser on cloud GPUs |
| Redbird / Redbirdpy | diffuse optical tomography solvers that can use MCX/MMC as the forward model |
Scripts move between MCX and MCX-CL by changing mcxlab to mcxlabcl, or pmcx to
pmcxcl; in MATLAB, mcxlab(cfg, 'opencl') also runs MCX-CL.
Installation and configuration
What you need
MCX-CL needs an OpenCL runtime, the part of a graphics or CPU driver that runs
OpenCL programs (OpenCL.dll on Windows, libOpenCL.so on Linux, built into macOS).
No SDK or compiler is needed: MCX-CL compiles its kernel for your device when it
starts.
| Processor | OpenCL runtime |
|---|
| NVIDIA GPU | included in the NVIDIA driver; CUDA is not needed |
| AMD GPU | included in the AMD Adrenalin driver (Windows); the ROCm OpenCL runtime (Linux) |
| Intel GPU (integrated or Arc) | included in the Intel graphics driver (Windows); intel-opencl-icd (Linux) |
| Apple silicon / Intel Mac | built into macOS |
| CPU | POCL (pocl-opencl-icd on Debian/Ubuntu, works on any x86_64 or ARM CPU), or Intel's CPU runtime |
On Linux, an ICD loader (ocl-icd-libopencl1 on Debian/Ubuntu) dispatches OpenCL
calls to the installed runtimes, and clinfo lists the devices they provide. Several
runtimes can be installed together, for example a GPU driver and POCL.
As a rough guide, high-end dedicated GPUs are 20-60x faster than a multi-core CPU, and
integrated GPUs 3-4x; see https://mcx.space/computebench/ for measured speeds.
Getting MCX-CL
PMCXCL wheels are built for Linux (x86_64 and aarch64), macOS (arm64 and x86_64) and
Windows; pmcx provides the utility functions PMCXCL uses. The released Linux binaries
and mex files run on glibc 2.17 or newer; the Octave mex files for Windows and macOS
need Octave 10 or newer.
The MCX Studio package
The all-in-one package unzips into one folder:
MCXStudio/
├── mcxstudio the GUI (MCXStudio.app on macOS)
├── mcxsuite_addpath.m adds all the MATLAB/Octave toolboxes to the path
├── MATLAB/ mcxlab, mcxlabcl, mmclab, iso2mesh and redbird
└── MCXSuite/ mcx, mcxcl and mmc, each with its bin/ folder and helper scripts
Add MCXSuite/mcxcl/bin (and mcx/bin, mmc/bin) to your PATH. In MATLAB or
Octave, cd into the MCXStudio folder, run mcxsuite_addpath, and savepath to keep
the paths for later sessions. On macOS, unzip the package under your home folder.
Checking your devices
Unzip the package, add its bin/ folder to your PATH, and run
mcxcl -L # list the OpenCL devices
mcxcl -Q cube60 # run the "cube60" benchmark on the first device
mcxcl -L prints each platform (runtime) and its devices, numbered from 1 across all
platforms:
Platform [0] Name NVIDIA CUDA
============ GPU device ID 0 [1 of 2]: NVIDIA TITAN V ============
Device 1 of 2: NVIDIA TITAN V
Compute units : 80 core(s)
Global memory : 12625903616 B
...
Auto-thread : 327680
Auto-block : 64
Use the device number with -G: -G 2 runs on the 2nd device, and -G 0110 on the
2nd and 3rd together. If no device is listed, the OpenCL runtime is missing or broken:
install or update the driver (or POCL for a CPU) and check with clinfo. Without
clinfo, you can check that the runtime library exists:
- Windows:
OpenCL.dll in C:\Windows\System32;
- macOS:
/System/Library/Frameworks/OpenCL.framework exists;
- Linux:
ldconfig -p | grep libOpenCL lists libOpenCL.so.1; if the library is in
a non-standard folder, add that folder to LD_LIBRARY_PATH, and check that the
vendor's .icd file is in /etc/OpenCL/vendors.
Intel CPUs older than Ivy Bridge (3rd generation Core, 2012) are not well supported by
Intel's OpenCL runtimes; POCL is an alternative for them.
Setting up MCXLAB-CL and PMCXCL
addpath('/path/to/mcxlabcl'); % mcxcl.mex*, mcxlabcl.m
addpath('/path/to/mcxcl/utils'); % helper functions
info = mcxlabcl('gpuinfo') % list the devices
python3 -m pip install pmcxcl pmcx jdata bjdata matplotlib
python3 -c "import pmcxcl; print(pmcxcl.gpuinfo())"
The MATLAB helpers in mcxcl/utils are a subset of those of MCX; for the full set
(mcxpreview, mcxplotshapes, mcxsvmc, ...), also add MCX's utils folder, which
works with MCXLAB-CL's outputs.
Operating system notes
Windows: Windows resets a GPU that is busy for more than a few seconds. The MCX
Studio installer raises this limit; with the zip packages, right-click
setup/win64/apply_timeout_registry_fix.bat, choose Run as administrator, and
reboot.
Linux laptops with hybrid NVIDIA graphics: if the screen freezes when MCX-CL runs,
use the NVIDIA GPU as the primary GPU (sudo prime-select nvidia). If simulations fail
after the laptop wakes up from suspension, close MATLAB/Python and run
sudo rmmod nvidia-uvm; sudo modprobe nvidia-uvm.
macOS: Apple's OpenCL runs MCX-CL on Apple silicon and Intel Macs; PMCXCL is the
way to run simulations from Python on a Mac. MCX Studio is provided as a universal
MCXStudio.app.
Building from source
git clone https://github.com/fangq/mcxcl.git
cd mcxcl/src
make # the mcxcl binary in ../bin
make mex # MCXLAB-CL for MATLAB: ../mcxlabcl/mcxcl.mex*
make oct # MCXLAB-CL for Octave: ../mcxlabcl/mcxcl.mex
make unittest # host unit tests
make fixed # exact 64-bit fixed-point output sums (-DUSE_FIXEDPOINT)
cd ../pmcxcl && python3 -m pip install . # PMCXCL
Building needs a C/C++ compiler and the OpenCL headers and loader (opencl-headers
and ocl-icd-opencl-dev on Debian/Ubuntu); make oct needs liboctave-dev, and the
Python module needs cmake. The kernel source, src/mcx_core.cl, is embedded in the
binary.
Quick start: one simulation, four ways
The same simulation - a pencil beam entering a 60x60x60 mm homogeneous cube with four
detectors - in each mode.
Command line. Save this as qtest.json:
{
"Domain": {
"Dim": [60, 60, 60],
"OriginType": 1,
"LengthUnit": 1,
"Media": [
{"mua": 0.00, "mus": 0.0, "g": 1.00, "n": 1.0},
{"mua": 0.005, "mus": 1.0, "g": 0.01, "n": 1.37}
]
},
"Shapes": [
{"Grid": {"Tag": 1, "Size": [60, 60, 60]}}
],
"Session": {
"Photons": 1000000,
"RNGSeed": 29012392,
"ID": "qtest"
},
"Forward": {
"T0": 0.0e+00,
"T1": 5.0e-09,
"Dt": 5.0e-10
},
"Optode": {
"Source": {
"Type": "pencil",
"Pos": [29.0, 29.0, 0.0],
"Dir": [0.0, 0.0, 1.0]
},
"Detector": [
{"Pos": [29.0, 19.0, 0.0], "R": 1.0},
{"Pos": [29.0, 39.0, 0.0], "R": 1.0},
{"Pos": [19.0, 29.0, 0.0], "R": 1.0},
{"Pos": [39.0, 29.0, 0.0], "R": 1.0}
]
}
}
and run
mcxcl -f qtest.json -G 1
MCX-CL writes qtest.jnii (the fluence rate, 60x60x60 voxels x 10 time gates) and
qtest_detp.jdt (the detected photons).
MATLAB/Octave:
cfg.nphoton = 1e6;
cfg.vol = uint8(ones(60, 60, 60));
cfg.prop = [0 0 1 1; 0.005 1 0.01 1.37]; % [mua(1/mm), mus(1/mm), g, n] per label
cfg.issrcfrom0 = 1;
cfg.srcpos = [29 29 0];
cfg.srcdir = [0 0 1];
cfg.detpos = [29 19 0 1; 29 39 0 1; 19 29 0 1; 39 29 0 1];
cfg.tstart = 0; cfg.tend = 5e-9; cfg.tstep = 5e-10;
cfg.gpuid = 1; % device number from mcxlabcl('gpuinfo')
[flux, detp] = mcxlabcl(cfg);
imagesc(log10(squeeze(sum(flux.data(:, 30, :, :), 4)))); % CW fluence rate at y=30
Python:
import numpy as np, pmcxcl, matplotlib.pyplot as plt
cfg = {
'nphoton': 1000000,
'vol': np.ones([60, 60, 60], dtype='uint8'),
'prop': [[0, 0, 1, 1], [0.005, 1, 0.01, 1.37]],
'issrcfrom0': 1,
'srcpos': [29, 29, 0], 'srcdir': [0, 0, 1],
'detpos': [[29, 19, 0, 1], [29, 39, 0, 1], [19, 29, 0, 1], [39, 29, 0, 1]],
'tstart': 0, 'tend': 5e-9, 'tstep': 5e-10,
'gpuid': 1,
}
res = pmcxcl.mcxlab(cfg)
plt.imshow(np.log10(res['flux'][:, 29, :, :].sum(axis=2)))
plt.show()
GUI: start mcxstudio, open qtest.json, choose the intel, amd or cpu
processor (which run mcxcl), and press Run.
Key concepts
Grid and units. The domain is a grid of Nx*Ny*Nz cubic voxels of edge length
LengthUnit (JSON), -u (command line) or cfg.unitinmm (MATLAB/Python), in mm.
Positions of sources and detectors are given in voxel units. With
OriginType: 1 (-z 1, cfg.issrcfrom0=1), the corner of the first voxel is at
[0,0,0]; otherwise (the default in the command line and MCXLAB-CL) it is at [1,1,1],
matching MATLAB's 1-based indices. Absorption (mua) and scattering (mus)
coefficients are in 1/mm regardless of the voxel size.
Media. Each voxel holds a medium label. The optical properties of label i are
the (i+1)-th row of Media (cfg.prop): [mua, mus, g, n], with mus the scattering
(not reduced scattering) coefficient, mus=mus'/(1-g). Label 0 is the background
(outside of the domain): a photon moving from a non-zero voxel into a 0 voxel
leaves the domain. Instead of labels, a voxel can also store its own optical
properties, see Volume and media formats.
Time gates. A simulation covers the time window [T0, T1], divided into gates of
width Dt; the output has one volume per gate. Photons are tracked until T1, so for
a CW solution choose T1 long enough for the light to decay (a few ns in tissue),
and sum the gates.
Photons and statistics. Results are averages over the simulated photon
packets; their noise falls as 1/sqrt(number of photons). 10^6-10^8 photons are
typical. Russian roulette (-e) terminates packets of low weight without biasing
the result.
Normalization. By default the output is normalized to a unit source, so MCX
produces the Green's function of the fluence rate (1/(mm^2 s), or W/mm^2 per watt
of source power), directly comparable with analytical solutions. -U 0
(cfg.isnormalized=0) returns the raw accumulated weights.
Boundaries. At the domain boundary, photons are reflected according to the
Fresnel equations when the refractive indices differ (-b 1, the default), or
leave the domain (-b 0). -B/cfg.bc sets each face separately, including
mirror and cyclic boundaries.
Running the mcxcl command line
A simulation needs a JSON input file (-f) or a built-in benchmark (-Q).
Command-line options override the corresponding settings in the file:
mcxcl -f input.json [options]
mcxcl -Q cube60 -n 1e7 -G 1
mcxcl (or mcxcl -h) prints all options. They are the same as those of mcx,
except for the device options; the most useful ones are below.
Input and simulation control
| Option | Meaning |
|---|
-f file, --input | read a JSON input file; a string starting with { is read as inline JSON, -f - reads stdin |
-Q name, --bench | run a built-in benchmark (-Q alone lists them): cube60, cube60b, cube60planar, cubesph60b, onecube, twocube, skinvessel, sphshells, spherebox, colin27 |
-N name, --net | download and run a simulation shared on NeuroJSON.io |
-j '{...}', --json | modify any setting of the input with a JSON fragment |
-P '{...}', --shapes | add shapes to the domain, see Shape constructs |
-n num, --photon | number of photons |
-r num, --repeat | repeat the simulation num times; negative: split the photons into |num| runs |
-E seed, --seed | RNG seed; a _detp.jdt file name replays its detected photons |
-u mm, --unitinmm | voxel edge length in mm |
-z 0/1, --srcfrom0 | 1: the first voxel's corner is [0,0,0]; 0: [1,1,1] |
-K fmt, --mediabyte | volume data format, see Volume and media formats |
--gscatter n | after n scattering events, scatter isotropically (faster) |
Physics, boundaries and output
| Option | Meaning |
|---|
-b 1/0, --reflect | Fresnel reflection at the domain boundary (1) or none (0) |
-B '______', --bc | per-face boundary conditions: r Fresnel, a absorbing, m mirror, c cyclic, _ use -b; 6 more digits mark faces as detectors |
-V 0/1, --specular | include the specular reflection when the source is outside the domain |
--voidtime 1/0 | for a source outside the domain, start the clock at launch (1) or at entry (0) |
-e w, --minenergy | Russian roulette threshold |
-U 1/0/2, --normalize | normalize the output (1), raw (0), or in replay skip the detected-weight normalization (2) |
-s name, --session | name of the output files |
-O X, --outputtype | the output quantity, see Output quantities |
-F fmt, --outputformat | jnii (default), bnii, nii, hdr, mc2, tx3 |
-Z zlib, --zip | compression for jnii/bnii: zlib, gzip, base64, lzip, lzma, lz4, lz4hc |
-d 1/0/3, --savedet | save detected photons (1), not (0), or stop when -H photons are detected (3) |
-w DP, --savedetflag | detected photon data to save, see Detectors |
Devices and the OpenCL kernel
| Option | Meaning |
|---|
-L, --listgpu | list the OpenCL devices and exit |
-l [1|0|none], --log | 1 (or -l alone): write the messages to session.log; 0: print them without the banner; none: print nothing, not even the progress bar |
-G n, -G '0110' or -G 'name:N', --gpu | use device n (from 1, as numbered by -L), or several devices given as a 0/1 mask, or by name: -G 'name:N' is the N-th device whose name contains name (ignoring case), -G name the first, e.g. -G 4070, -G radeon:2 |
-W '50,30,20', --workload | relative share of photons for each selected device |
-A 1/0, --autopilot | choose the thread count, dispatch photons dynamically and use fast math (1, default), or use -t/-T with strict math (0) |
-t num, --thread / -T num, --blocksize | total threads and threads per work group when -A 0 |
-o 3, --optlevel | kernel optimization level, 0 (none) to 3 |
-J '-D...', --compileropt | extra OpenCL compiler options, e.g. -DMCX_GPU_DEBUG (print each step) or -cl-nv-verbose (register counts on NVIDIA) |
-k file.cl, --kernel | run your own copy of the OpenCL kernel instead of the built-in one |
--showkernel 1 | print the built-in (or -k) kernel source |
--internalsrc 1 | the source is inside the domain: skip the search for the entry point |
-I, --printgpu | print the device information, then run |
Recipes
mcxcl -L # which devices do I have?
mcxcl -Q cube60b -n 1e8 -G 0111 -W 40,40,20 # three devices, unequal shares
mcxcl -f input.json -s run1 -F bnii -Z lz4 # compressed binary output
mcxcl -Q cube60 -A 0 -t 4096 -T 64 -E 1234 # bit-wise reproducible run
mcxcl -Q cube60 --showkernel 1 > mykernel.cl # get the kernel, edit it, then
mcxcl -Q cube60 -k mykernel.cl # ... run the edited kernel
mcxcl -f input.json -E run1_detp.jdt -O J # replay run1's detected photons for a Jacobian
JSON input file format
An MCX/MCX-CL input file is a JSON object with four required sections, Domain,
Session, Forward and Optode, and an optional Shapes section. Field names are
case sensitive; the order of fields does not matter, and unknown sections (such as
Help) are ignored. A field can also be written in the root object with a dotted
name, e.g. "Session.Photons": 1e7. Command-line options take precedence over the
file, and -j edits any field without touching the file.
The fields MCX-CL reads:
Domain - the grid and its media
| Field | Meaning |
|---|
VolumeFile | a binary volume file (format set by MediaFormat or -K), a JSON/JNIfTI (.json, .jnii, .bnii) volume, or a shape file; omit it when Shapes builds the volume |
Dim | [Nx, Ny, Nz] |
MediaFormat | voxel format of the volume, "byte" by default (see below) |
Media | array of {"mua":..., "mus":..., "g":..., "n":...}, the first for label 0 |
LengthUnit | voxel size in mm (-u) |
OriginType | 1: first voxel's corner at [0,0,0]; 0: at [1,1,1] (-z) |
MieScatter | polarized MC: array of {"mua", "radius" (micron), "rho" (1/micron^3), "nsph", "nmed"}, one per label from 1 |
InverseCDF | a user-defined scattering phase function, the inverse CDF of cos(theta) |
AngleInverseCDF | a user-defined launch-angle distribution (also accepted under Optode.Source) |
Session - what to compute and save
| Field | Meaning |
|---|
ID | name of the output files (-s) |
Photons | number of photons (-n) |
RNGSeed | RNG seed (-E) |
OutputType | output quantity: "flux", "fluence", "energy", "jacobian", ... or the -O letter |
OutputFormat | "jnii", "bnii", "nii", "hdr", "mc2", "tx3" (-F) |
DoMismatch | Fresnel reflection at the boundary (-b) |
BCFlags | per-face boundary conditions (-B) |
DoNormalize | normalize the output (-U) |
DoSaveVolume | save the volume output (-S) |
DoPartialPath | save detected photons (-d) |
SaveDataMask | detected photon data to save, e.g. "DPX" (-w) |
DoSaveExit, DoSaveRef, DoSaveSeed, DoDCS | as -x, -X, -q and -m |
DoSpecular | specular reflection at entry (-V) |
DoAutoThread | autopilot (-A) |
MaxDetPhoton | the most detected photons kept (-H) |
MinEnergy | Russian roulette threshold (-e) |
DebugFlag | debug flags such as "P" (-D) |
RootPath | folder of the input files () |
Forward - time window: T0 (start), T1 (end) and Dt (gate width), in seconds.
Optode - sources and detectors
| Field | Meaning |
|---|
Source.Type | source type, see Sources |
Source.Pos | [x, y, z] in voxel units |
Source.Dir | [vx, vy, vz] unit direction, optionally a 4th value: focal length (nan isotropic, -inf Lambertian) |
Source.Param1, Source.Param2 | 4 numbers each, depending on Type |
Source.Pattern | the pattern of pattern/pattern3d sources: {"Nx", "Ny", "Nz", "Data"}, or a JData array |
Source.SrcNum | number of patterns simulated together (photon sharing) |
Source.Weight | initial photon weight |
Source.Frequency | modulation frequency (Hz) for frequency-domain (RF) outputs |
Source.IQUV, Source.WaveLength | Stokes vector and wavelength (nm) for polarized MC |
Detector | array of {"Pos": [x,y,z], "R": radius}, and "Dir" for the adjoint outputs |
SDMap | channels of the adjoint Jacobians: an array of [srcid, detid] arrays (1-based; extra values ignored); default: all source-detector pairs |
Optode.Source can also be an array of sources, each with its own Type, Pos,
Dir and parameters. Several examples are in the example/ folder, for example
example/quicktest, example/skinvessel and example/multisrc. mcxcl --dumpjson
prints the complete input of any simulation, including the defaults, which is a
good template.
Legacy files: mcx2json.m and json2mcx.m (MATLAB) or pmcxcl.mcx2json /
pmcxcl.json2mcx convert between an MCXLAB-CL/PMCXCL cfg and a JSON input file, so a
simulation set up in MATLAB or Python can be run with the command line, and vice
versa.
Shape constructs
Instead of supplying a volume file, the domain can be built from simple shapes. The
Shapes array is processed in order, each shape overwriting the voxels it covers,
and can be placed in the input file, in a separate .json file named by
Domain.VolumeFile, or on the command line with -P:
"Shapes": [
{"Name": "head model"},
{"Origin": [0, 0, 0]},
{"Grid": {"Tag": 1, "Size": [60, 60, 60]}},
{"ZLayers": [[1, 5, 1], [6, 8, 2], [9, 60, 3]]},
{"Sphere": {"Tag": 4, "O": [30, 30, 30], "R": 10}},
{"Box": {"Tag": 0, "O": [10, 10, 10], "Size": [10, 10, 10]}},
{"Cylinder": {"Tag": 2, "C0": [0, 30, 20], "C1": [60, 30, 20], "R": 3}}
]
| Shape | Parameters | Result |
|---|
Grid | Tag, Size: [Nx,Ny,Nz] | (re)creates the grid, filled with Tag |
Origin | [x,y,z] | origin of the coordinates of the shapes that follow |
Sphere | Tag, O (center), R | a sphere |
Box | Tag, O (lower corner), Size | a box, floating-point coordinates |
Subgrid | Tag, O, Size | a box of whole voxels, integer indices, both ends included |
Cylinder | Tag, C0, C1 (axis ends), R | a finite cylinder |
XLayers, YLayers, ZLayers | [[start, end, tag], ...] | layers normal to an axis, 1-based voxel indices, ends included |
XSlabs, YSlabs, ZSlabs | Tag, Bound: [[start, end], ...] | slabs normal to an axis, floating-point bounds |
UpperSpace | Tag, Coef: [A,B,C,D] | the half space A*x+B*y+C*z>D |
Name | a string | a comment |
Tag is the medium label, and must be smaller than the number of media. Coordinates
start from the corner of the first voxel, at [0,0,0] or [1,1,1] depending on
OriginType, unless an Origin is given. In MATLAB/Python, the same JSON string can
be passed as cfg.shapes; mcxplotshapes.m (MCX's utils) and pmcxcl.plotshapes draw it.
MCX's example/shapetest has more examples.
Volume and media formats
The domain is usually a label volume, but MCX-CL accepts several voxel formats, set by
Domain.MediaFormat in JSON or -K on the command line. In MCXLAB-CL and PMCXCL, the
format follows from the type and shape of cfg.vol:
-K | MediaFormat | MCXLAB-CL/PMCXCL cfg.vol | Each voxel stores |
|---|
| 1 | byte | 3D uint8 | a label, 0-255 |
| 2 | short | 3D uint16 | a label, 0-65535 |
| 4 | integer | 3D uint32 (or single/double) | a label |
| 101 | mua_float | 1xNxxNyxNz single | mua; mus, g, n from label 1 |
| 100 | muamus_float | - | mua and mus; g, n from label 1 |
| 102 | muamus_half | 2xNxxNyxNz single | mua and mus, stored as half floats |
| 104 | muamus_short | 2xNxxNyxNz uint16 | mua and mus as 0-65535 levels between labels 1 and 2 |
| 103 | asgn_byte | 4xNxxNyxNz uint8 | mua, mus, g, n as 0-255 levels between labels 1 and 2 |
| 96 | asgn_float | 4xNxxNyxNz single | mua, mus, g, n |
| 99 | labelplus | 3xNxxNyxNz single | a label, plus one property of it replaced by a value |
| 98 | mixlabel | 3xNxxNyxNz uint16 | two labels and the fraction of the first |
| 97 | svmc | 8xNxxNyxNz uint8 | split-voxel MC: two labels and the interface plane in the voxel |
Volumes with per-voxel properties (all but the label and svmc formats) carry no
per-medium partial path data. Binary volume files are read in column-major (MATLAB)
order unless -a 1 is given. Examples: demo_continuous_mua_mus.m and the
demo_svmc_*.m scripts.
Sources
The source is set by its type, position (Pos/srcpos), direction (Dir/srcdir)
and two 4-element parameter vectors (Param1/srcparam1, Param2/srcparam2).
The 4th element of the direction is a focal length for the sources marked [*]:
0 for a collimated beam, positive to converge to a focus that far along the
direction, negative to diverge from a virtual point behind the source, nan for an
isotropic launch and -inf for a Lambertian (cosine) launch (3D only).
| Type | Description and parameters |
|---|
pencil | a pencil beam (default) |
isotropic | an isotropic point source |
arcsine | like isotropic, with a uniformly distributed zenith angle |
cone | a uniform cone; Param1[0]: half angle (rad) |
zgaussian | a Gaussian distribution of the zenith angle; Param1[0]: variance |
gaussian [*] | a collimated Gaussian beam; Param1[0]: waist radius |
hyperboloid [*] | a hyperboloid Gaussian beam; Param1: [waist radius, distance to focus, Rayleigh range] |
disk [*] | a uniform disk normal to Dir; Param1[0]: radius, Param1[1]: inner radius of an annulus |
ring [*] | a ring or ring sector; Param1: [outer radius, inner radius, start angle, end angle] |
planar [*] | a uniform quadrilateral with corners Pos, Pos+Param1[0:3], Pos+Param2[0:3] |
pattern [*] | the planar area with intensities from Pattern (Param1[3] x Param2[3] pixels, values 0-1) |
pattern3d [*] | a 3D intensity volume of size Param1[0:3] starting at Pos |
fourier [*] | a spatial-frequency (SFDI) pattern 0.5*[1+M*cos(2*pi*(fx*x+fy*y)+phi0)] over the planar area; fx, fy and phi0, M from Param1[3], Param2[3] |
fourierx [*] | a general Fourier source; Param1: [v1, |v2|], Param2: [kx, ky, phi0, M] |
fourierx2d [*] | a general 2D Fourier basis; Param1: [v1, |v2|], Param2: [kx, ky, phix, phiy] |
pencilarray | an array of pencil beams over the planar area; Param1[3] x Param2[3] beams |
line | a line source from Pos to Pos+Param1[0:3], emitting normal to the line; Param2[0]: half angle of a wedge (fan) beam |
slit [*] | a collimated slit beam along Dir from the same segment; Param2[0:1]: Gaussian broadening |
Several sources are given as an array of Optode.Source objects, or as one row
per source in srcpos/srcdir/srcparam1/srcparam2. By default they are
simulated together into one output; --srcid -1 (cfg.srcid=-1) gives each source
its own output, and --srcid N simulates only the N-th source. Example:
demo_multisrc.m.
Initial weight. A 4th element of srcpos (or Source.Weight) sets the initial
photon weight, which may be negative; the output scales linearly with it.
Photon sharing. A pattern source with SrcNum/cfg.srcnum > 1 simulates
many illumination patterns in a single run, each photon carrying a weight for
every pattern; the output gains a pattern dimension. Example:
demo_photon_sharing.m [Yan2020-OL].
Custom angular distributions. Source.AngleInverseCDF (cfg.angleinvcdf)
gives the launch angle distribution as the inverse CDF of theta/pi.
Example: demo_mcxlab_launchangle.m.
mcxpreview(cfg) (MATLAB, from MCX's utils) and pmcxcl.preview(cfg) (Python) draw the domain,
sources and detectors, which is the quickest way to check a source setup;
demo_mcxlab_srctype.m shows every source type.
Detectors and detected photon data
Detectors are disks given by a center and a radius (Optode.Detector,
cfg.detpos = [x y z r]); a photon leaving the domain within a detector is
recorded. Faces of the bounding box can also act as detectors, through the 7th-12th
characters of -B/cfg.bc (see demo_bc_det.m).
For each detected photon, MCX-CL can save the fields selected by -w
(Session.SaveDataMask, cfg.savedetflag), default DP:
| Letter | Field | Columns |
|---|
D | detector ID | 1 |
S | scattering events in each medium | number of media |
P | partial path length in each medium (voxel units) | number of media |
M | momentum transfer in each medium (for DCS) | number of media |
X | exit position | 3 |
V | exit direction | 3 |
W | initial weight | 1 |
I | Stokes vector (polarized MC) | 4 |
These data make it possible to compute detector readings for different optical
properties without rerunning the simulation (perturbation Monte Carlo): the
weight of a detected photon is w0*exp(-sum(mua_i*L_i)). The helper functions do
this for you: mcxdetweight/pmcxcl.detweight (weights), mcxmeanpath/
pmcxcl.meanpath (mean partial path lengths), mcxdettpsf/pmcxcl.dettpsf (temporal
point spread function), mcxcwdref/pmcxcl.cwdref (CW reflectance) and
mcxdcsg1/pmcxcl.dcsg1 (DCS autocorrelation g1).
Saving the RNG seed of each detected photon (-q 1, or requesting the 4th output in
MCXLAB-CL) enables photon replay, see Advanced features.
Output quantities, files and formats
Output quantities
-O (Session.OutputType, cfg.outputtype) selects what the volume output holds:
-O | Name | Quantity |
|---|
X | flux | fluence rate, 1/(mm^2 s) per unit source (default) |
F | fluence | fluence, 1/mm^2, integrated over each time gate |
E | energy | energy deposited per voxel |
L | length | total path length per voxel |
J | jacobian (wl) | replay: mua Jacobian |
P | nscat (wp) | replay: weighted scattering counts, for the mus Jacobian |
M | wm | replay: weighted momentum transfer |
R | rf | replay: frequency-domain mua Jacobian; without replay: complex fluence rate |
S | rfmus | replay: frequency-domain mus Jacobian |
T | wltof | replay: time of flight x path length |
B | wptof | replay: time of flight x scattering count |
A | adjoint | adjoint mua Jacobian |
D | adjoint_dcoeff | adjoint diffusion coefficient Jacobian |
U | adjoint_mus | adjoint mus Jacobian |
V | adjoint_musp | adjoint mus' Jacobian |
W | adjoint_mua_d | adjoint mua and D Jacobians together |
Q | adjoint_mua_musp | adjoint mua and mus' Jacobians together |
The adjoint Jacobians have one channel per source-detector pair, [Nx,Ny,Nz,Ns*Nd]
with the detector index varying the fastest; Optode.SDMap (cfg.sdmap), a list of
1-based [srcid, detid] rows, computes only those channels, in that order.
Interpreting the fluence
The volume output is an Nx*Ny*Nz*Ng array of the chosen quantity (with an extra
dimension for photon-sharing patterns, separate sources or replayed detectors). The
default fluence rate is the Green's function for a unit-energy impulse at t=0: the
fraction of the source energy (or photons) passing through a unit area per unit
time. Each gate stores the average over its time window, i.e. approximately the
value at the center of the window, t0+(k-1/2)*Dt. To obtain
- the fluence of each gate: multiply the fluence rate by Dt (or use
-O F);
- the CW fluence: sum the fluence rate over all gates and multiply by Dt, with
T1 long enough for the light to decay;
- the energy deposition: use
-O E, or multiply the fluence by mua
(mcxfluence2energy.m).
The diffuse reflectance can be obtained from the detected photons (the detected
weight divided by the detector area gives the current density), or with -X 1
(cfg.issaveref=1), which records the escaping weight in the background voxels
next to the domain (stored as negative values in the raw volume; MCXLAB-CL/PMCXCL
return it separately as dref). Pad the volume with a layer of 0 voxels where you
want the reflectance. The volume output keeps the voxel order of the input volume
(column-major, or row-major with -a 1). example/validation compares MCX-CL with
analytical solutions (cwdiffusion.m, tddiffusion.m, pmcxcl.cwdiffusion, ...).
Units of the output
With the default normalization (-U 1), the output is the response to a source that
emits one unit of light (one joule, or one photon) in an impulse at t=0, so each
quantity can be read in energy or in photon units:
| Output | Energy units | Photon units |
|---|
fluence rate (flux, default) | W/mm^2 per J emitted | 1/(mm^2 s) per photon emitted |
fluence (fluence) | J/mm^2 per J emitted | 1/mm^2 per photon emitted |
energy (energy) | J absorbed in the voxel per J emitted | the fraction of the photons absorbed there |
To obtain absolute values, multiply by the energy (or number of photons) your source
delivers; for a continuous source of power P, the CW fluence rate is P times the
time-integrated (summed over gates and multiplied by Dt) fluence rate. Divide the
energy output by the voxel volume (LengthUnit^3) for the absorbed energy density;
mcxfluence2energy.m converts a fluence output into energy deposition. See [Fang2009] and Lihong Wang's Biomedical Optics: Principles
and Imaging, chapter 5, for the definitions.
Fluence rate and current density. The fluence rate is a scalar: the light passing
through a small sphere around a point from all directions. The current density J is a
vector: the net light flowing through a surface in a given direction; it has the same
units. What a detector or the diffuse reflectance measures is the current density
flowing out of the domain: the total weight of the detected photons divided by the
detector area gives J at the detector, as a function of time when the photons are
binned by their time of flight.
Output files of the command line
| File | Content |
|---|
session.jnii (-F jnii, default) | the volume output, JSON-based JNIfTI with compressed binary data |
session.bnii (-F bnii) | the same in binary JSON (BJData), smaller and faster to read |
session.nii / .hdr+.img (-F nii/hdr) | NIfTI-1 / Analyze 7.5 |
session.mc2 (-F mc2) | raw float32 array, column-major, no header (legacy) |
session_detp.jdt / .jdb | detected photons (and seeds with -q 1), with jnii/bnii |
session.mch | detected photons in the legacy binary format, with mc2/nii/hdr |
session_traj.jdt / .jdb (or .mct) | photon trajectories, with -D M |
session.log | the messages, with -l |
JNIfTI (.jnii, .bnii) and JData (.jdt, .jdb) files follow the
NeuroJSON specifications: they are plain JSON or binary
JSON, readable in any language, with the array size, type and compression stored
alongside the data. Detected photon files have the structure
{"MCXData": {
"Info": {"Version", "MediaNum", "DetNum", "Media", ...},
"PhotonData": {"detid", "srcid", "nscat", "ppath", "mom", "p", "v", "w0"},
"Seed": [...]
}}
with only the fields selected by -w (srcid only with several sources, Seed
only with -q 1). Trajectory files hold a Trajectory object with photonid,
p, w0, srcid and, with --trajstokes 1, iquv. Files with the older .jdat
suffix are still accepted for replay.
The legacy .mch/.mct files have a 256-byte header (MCXHistoryHeader: magic
MCXH, version, number of media, detectors, columns, photons simulated, detected
and saved, voxel size, seed bytes, normalizer, respin, srcnum, savedetflag, total
sources) followed by a float32 array of saved photons x columns, and with -q 1
the seeds; loadmch.m reads them.
Messages
MCX-CL prints its settings, the device(s) used, the time spent in each phase (ms) and,
at the end, the simulation speed (photons/ms) and the absorbed fraction. -D P
(cfg.debuglevel='P') adds a progress bar.
Reading MCX-CL data in MATLAB and Python
MATLAB/Octave - install JSONLab (or
octave-jsonlab on Debian/Ubuntu/Fedora) for JNIfTI/JData files; the other loaders
are in MCX's utils folder:
nii = loadjd('qtest.jnii'); % JSONLab: .jnii, .bnii, .jdt, .jdb, ...
flux = nii.NIFTIData; % Nx x Ny x Nz x Ng
detp = loadjd('qtest_detp.jdt');
ppath = detp.MCXData.PhotonData.ppath; % one row per detected photon
traj = loadjd('qtest_traj.jdt');
flux = loadmc2('qtest.mc2', [60 60 60 10]); % utils: .mc2 needs the dimensions
[data, header] = loadmch('qtest.mch'); % utils: legacy detected photons
nii = mcxloadnii('qtest.nii'); % utils (MATLAB), or niftiread
Python - install jdata (pip install jdata bjdata, or python3-jdata on
Debian/Ubuntu):
import jdata as jd
nii = jd.load('qtest.jnii') # also .bnii, .jdt, .jdb (needs bjdata)
flux = nii['NIFTIData'] # numpy array, Nx x Ny x Nz x Ng x 1
detp = jd.load('qtest_detp.jdt')['MCXData']['PhotonData'] # dict of arrays
import pmcxcl
flux = pmcxcl.loadmc2('qtest.mc2', [60, 60, 60, 10])
data, header = pmcxcl.loadmch('qtest.mch')[:2] # a 3rd value holds the seeds, if saved
.nii files can be read with nibabel.
Using MCXLAB-CL in MATLAB and Octave
MCXLAB-CL runs MCX-CL inside MATLAB or GNU Octave. It takes the same cfg struct as
MCXLAB and returns the same outputs:
[flux, detp, vol, seeds, traj] = mcxlabcl(cfg);
Request only the outputs you need. help mcxlabcl documents every cfg field, its
default and the outputs; the OpenCL-specific fields are cfg.gpuid (a device number,
or a 0/1 string selecting several), cfg.workload, cfg.optlevel,
cfg.compileropt, cfg.kernelfile and cfg.internalsrc (cfg.isatomic is still
accepted but ignored: the output is always accumulated with atomic operations).
| Output | Content |
|---|
flux.data | the volume output, Nx x Ny x Nz x Ng (plus extra dimensions for patterns, sources, detectors or RF parts) |
flux.dref | diffuse reflectance, with cfg.issaveref=1 |
flux.stat | runtime, photons, launched and absorbed energy, normalizer |
flux.jmua, .jd, .jmus, .jmusp | adjoint Jacobians |
detp | detected photons: detid, ppath, nscat, mom, p, v, w0, ... |
vol | the preprocessed volume |
seeds | RNG seeds of the detected photons, for replay |
traj | photon trajectories: id, pos, data |
Other calls: mcxlabcl('gpuinfo') lists the devices, mcxlabcl('version') returns the
version, mcxlabcl(cfg, 'preview') draws the domain, and mcxlabcl(cfg, 'cuda') runs
the simulation with MCX (mcx.mex) instead. A variable USE_MCXCL in the base
workspace switches both mcxlab and mcxlabcl: 1 uses MCX-CL, 0 uses MCX, and an
integer N>1 (or a device string such as '0101') also selects that device.
Helper functions in utils: mcxdetphoton, mcxdetweight, mcxmeanpath,
mcxmeanscat, mcxdettpsf, mcxdettime, mcxdcsg1 (detected photons);
cwdiffusion, tddiffusion (analytical solutions); mcxrfreplay; mcxplotphotons,
mcxplotvol, islicer, slice3i (plotting); mcx2json, json2mcx, loadmc2,
loadmch, mcxloadnii, mcxfluence2energy (files). MCX's utils folder adds
mcxpreview, mcxplotshapes, mcxsvmc, mcxcwdref, mcxcreate and more.
Examples: about 40 scripts in mcxlabcl/examples, from demo_mcxlabcl_basic.m to
atlas head models, SFDI, replay and adjoint Jacobians, polarized light and SVMC; they
are the MCXLAB examples running with MCX-CL.
Using PMCXCL in Python
PMCXCL brings MCX-CL to Python 3.6+, with the same interface as PMCX:
import pmcxcl
res = pmcxcl.run(cfg) # cfg: a dict with the MCXLAB field names
res = pmcxcl.run(nphoton=1e6, vol=vol, ...) # or keyword arguments
res = pmcxcl.mcxlab(cfg) # run(), then decode detp/traj into dicts
print(pmcxcl.gpuinfo()) # the OpenCL devices
res is a dict with flux, stat and, when produced, detp, seeds, traj,
dref, vol and the adjoint Jacobians. The utility functions - detweight,
meanpath, dettpsf, cwdref, dcsg1, cwdiffusion, tddiffusion, rfreplay,
loadmc2, loadmch, mcx2json, json2mcx, preview, plotshapes, plotphotons,
plotvol, svmc and more - come from pmcx and are available as pmcxcl.<name>.
All MCXLAB examples have Python versions in MCX's
pmcx/examples; they run with
MCX-CL after replacing import pmcx by import pmcxcl as pmcx. The
PMCX tutorial notebook
applies to PMCXCL as well. See pmcxcl/README.md for building from source.
Using MCX Studio
MCX Studio 2, included in the MCX Studio packages, is the graphical front end for
MCX, MCX-CL and MMC (source code in MCX's mcxstudio2/ folder). A simulation is the
same .json file mcxcl reads; the settings are grouped in the order a simulation is
described, with F1 help for each; a 3-D view shows the domain, sources, detectors,
fluence and photon paths. Choosing the intel, amd or cpu processor runs mcxcl:
mcxstudio qtest.json # edit and run
mcxstudio --run --processor intel qtest.json # run with mcxcl
See the MCX README and mcxstudio2/README.md in MCX for details.
Advanced features
MCX-CL provides the same advanced features as MCX; the MCX README describes them in
detail, and the MCXLAB-CL examples show them in use.
- Photon replay and Jacobians [Yao2018]: save the seeds of the detected photons
(
-q 1, or the 4th output of mcxlabcl), then replay them (-E file_detp.jdt, or
cfg.seed=seeds.data with cfg.detphotons) with -O J (mua), P (mus), R
(frequency domain), S, T or B, for one detector (-Y n) or all of them.
Examples: demo_mcxlab_replay.m, demo_replay_timedomain.m,
demo_replay_frequencydomain.m; example/rfreplay.
- Adjoint Jacobians: give each detector a direction (
Optode.Detector.Dir,
cfg.detdir) and use output type adjoint, adjoint_dcoeff, adjoint_mus,
adjoint_musp, adjoint_mua_d or adjoint_mua_musp, in CW or the frequency
domain. Examples: demo_validate_mcx_adjoint_*.m.
- Frequency-domain (RF) simulations: set the modulation frequency
(
Optode.Source.Frequency in Hz, cfg.omega in rad/s) with output type rf to get
the complex fluence rate. Example: demo_validate_mcx_rf_redbird.m.
- Photon sharing [Yan2020-OL]: a
pattern source with SrcNum/cfg.srcnum > 1
simulates many illumination patterns at once. Example: demo_photon_sharing.m.
- Polarized light [Yan2022-pol]:
Domain.MieScatter/cfg.polprop. Example:
demo_polarized_photon.m.
- Split-voxel MC [Yan2020] for curved boundaries, and per-voxel optical
properties. Examples:
demo_svmc_*.m, demo_continuous_mua_mus.m.
- Custom phase functions and launch angles:
cfg.invcdf, cfg.angleinvcdf.
Examples: demo_mcxlab_phasefun.m, demo_mcxlab_launchangle.m.
- Several devices:
-G 0111 -W 40,40,20 splits one simulation across devices,
which may come from different vendors; give a faster device a larger share.
- Custom kernels:
--showkernel 1 prints the kernel source; edit it and run it
with -k mykernel.cl to experiment with the algorithm without recompiling mcxcl.
MCX-CL and MCX: differences
MCX-CL and MCX share their input format, cfg fields, outputs and helper functions,
and are tested against each other step by step. The differences are:
| MCX | MCX-CL |
|---|
| hardware | NVIDIA GPUs (CUDA), AMD GPUs (HIP) | any OpenCL device: NVIDIA, AMD, Intel GPUs, Apple silicon, CPUs |
| kernel | compiled ahead of time | compiled for your device at startup (-o, -J, -k) |
-k option | --voidtime | --kernel (use --voidtime for the void timer) |
| 2D disk sources, Lambertian launch in 2D | supported | not yet (Lambertian launch is 3D only) |
| MATLAB | mcxlab | mcxlabcl, or mcxlab(cfg, 'opencl') |
| Python | pmcx | pmcxcl (utility functions from pmcx) |
Known device issues in this release:
- on the POCL CPU device, a pattern3d source run with the default dynamic photon
dispatch underestimates the absorbed energy at 1e6 photons or more; add
-A 0
(static split) for such runs. GPUs are not affected.
- on Apple silicon GPUs, about one in 20 to 60 runs may return an empty result
(a total simulated energy of 0); repeat the simulation if that happens.
Performance and best practices
- Pick the right device. Run
mcxcl -L and choose a dedicated GPU with -G; the
first device listed is not always the fastest. A discrete GPU is typically tens of
times faster than a CPU, and several devices can work together (-G, -W).
- Keep the autopilot on (
-A 1, the default): it sizes the thread count for the
device, hands out photons dynamically, and compiles the kernel with fast math
(-cl-unsafe-math-optimizations).
- Use
-A 0 to reproduce a run exactly. With dynamic dispatch and fast math, two
runs with the same seed agree statistically but not bit by bit. For identical
results, use -A 0 and keep the seed, -t and -T unchanged; saving or replaying
seeds always uses strict math and static dispatch.
- Expect a short startup delay. The kernel is compiled for the device each time
mcxcl starts; this takes from a fraction of a second to a few seconds.
- Avoid GPU timeouts. A GPU that also drives a display may be reset by the OS
after a few seconds of computing (see Operating system notes); use a dedicated GPU,
raise the timeout, or split the photons into shorter runs with
-r -N.
- CPU runs (POCL or Intel's runtime) use all CPU cores and are useful for testing
without a GPU.
- Use 1-byte labels, save only the outputs you need, and remember that the noise
falls as 1/sqrt(photons).
Testing
| Test | What it checks | How to run |
|---|
| host unit tests | parsing, JData, shapes, Mie and output writers, without a device | make unittest in src |
| regression tests | the command-line options and outputs | test/testmcx.sh |
| built-in examples | MCX's 141 example cases, against golden values | test/testbuiltin.sh |
| trajectory tests | step-by-step photon trajectories against MCX's golden data | test/testtraj.sh |
| PMCXCL tests | the Python module | python3 -m unittest discover -s test in pmcxcl/ |
| fuzzing | mutated built-in inputs and options: crashes, hangs, sanitizer and OpenCL kernel errors | test/fuzz/dfuzz.py |
The continuous integration runs all but the fuzzer on a CPU OpenCL device (POCL) with -A 0.
Getting help
Contributing
MCX-CL is open-source software under the GPLv3 license: you may use, modify and share
it, as long as derived work is also released under the GPLv3. Bug fixes and new
features are welcome:
- fork https://github.com/fangq/mcxcl, make your change in your fork, keep it small
and local, and open a pull request; describe what it fixes and how you tested it;
- or send the change as a
git diff to the mailing list or a GitHub issue.
Please format C/C++ code with make pretty (astyle) in src, MATLAB code with
mh_style --fix (miss_hit) and Python code with black before submitting, and add a
test when you can. Contributions are acknowledged in the What's New section and the
ChangeLog.
Publications and how to cite
If you use MCX-CL in your research, please cite [Yu2018], which describes MCX-CL, and
[Fang2009], the original MCX paper, as well as the papers of the specific features
you use. The complete
list is maintained at https://mcx.space/#publications.
Core papers
- [Fang2009] Qianqian Fang and David A. Boas, "Monte Carlo Simulation of Photon
Migration in 3D Turbid Media Accelerated by Graphics Processing Units," Opt.
Express, 17(22), 20178-20190 (2009).
https://opg.optica.org/oe/abstract.cfm?uri=oe-17-22-20178
- [Yu2018] Leiming Yu, Fanny Nina-Paravecino, David Kaeli, Qianqian Fang, "Scalable
and massively parallel Monte Carlo photon transport simulations for
heterogeneous computing platforms," J. Biomed. Opt. 23(1), 010504 (2018).
https://doi.org/10.1117/1.JBO.23.1.010504
- [Yan2020] Shijie Yan and Qianqian Fang, "Hybrid mesh and voxel based Monte Carlo
algorithm for accurate and efficient photon transport modeling in complex
bio-tissues," Biomed. Opt. Express, 11(11), 6262-6270 (2020).
https://opg.optica.org/boe/abstract.cfm?uri=boe-11-11-6262
Features
- [Hirvi2026] Pauliina Hirvi, Jaakko Olkkonen, Qianqian Fang, Ilkka Nissilä, "Computation
of frequency- and time-domain Jacobians in optical tomography with Monte Carlo
simulations," arXiv:2604.27558 (2026), preprint. https://arxiv.org/abs/2604.27558
- [Yao2018] Ruoyang Yao, Xavier Intes, Qianqian Fang, "A direct approach to compute
Jacobians for diffuse optical tomography using perturbation Monte Carlo-based
photon 'replay'," Biomed. Opt. Express, 9(10), 4588-4603 (2018).
https://opg.optica.org/boe/abstract.cfm?uri=boe-9-10-4588
- [Yan2020-OL] Shijie Yan, Ruoyang Yao, Xavier Intes, Qianqian Fang, "Accelerating
Monte Carlo modeling of structured-light-based diffuse optical imaging via
'photon sharing'," Opt. Lett., 45, 2842-2845 (2020).
- [Yan2022-pol] Shijie Yan, Steven L. Jacques, Jessica C. Ramella-Roman, Qianqian
Fang, "Graphics processing unit-accelerated Monte Carlo simulation of polarized
light in complex three-dimensional media," J. Biomed. Opt. 27(8), 083015 (2022).
https://doi.org/10.1117/1.JBO.27.8.083015
- [Fang2022-cloud] Qianqian Fang and Shijie Yan, "MCX Cloud - a modern, scalable,
high-performance and in-browser Monte Carlo simulation platform with cloud
computing," J. Biomed. Opt. 27(8), 083008 (2022).
https://doi.org/10.1117/1.JBO.27.8.083008
- [Yuan2018] Yaoshen Yuan, Leiming Yu, Zafer Doğan, Qianqian Fang, "Graphics
processing units-accelerated adaptive nonlocal means filter for denoising
three-dimensional Monte Carlo photon transport simulations," J. Biomed. Opt.
23(12), 121618 (2018).
- Matin Raayai Ardakani, Leiming Yu, David R. Kaeli, Qianqian Fang, "Framework for
denoising Monte Carlo photon transport simulations using deep learning," J.
Biomed. Opt. 27(8), 083019 (2022). https://doi.org/10.1117/1.JBO.27.8.083019
Mesh-based Monte Carlo (MMC), for the related mesh-based simulator
- Qianqian Fang, "Mesh-based Monte Carlo method using fast ray-tracing in Plücker
coordinates," Biomed. Opt. Express, 1(1), 165-175 (2010).
https://opg.optica.org/boe/abstract.cfm?uri=boe-1-1-165
- Qianqian Fang and Shijie Yan, "Graphics processing unit-accelerated mesh-based
Monte Carlo photon transport simulations," J. Biomed. Opt. 24(11), 115002 (2019).
https://doi.org/10.1117/1.JBO.24.11.115002
- Shijie Yan, Anh Phong Tran, Qianqian Fang, "A dual-grid mesh-based Monte Carlo
algorithm for efficient photon transport simulations in complex 3D media," J.
Biomed. Opt. 24(2), 020503 (2019). https://doi.org/10.1117/1.JBO.24.2.020503
- Yaoshen Yuan, Shijie Yan, Qianqian Fang, "Light transport modeling in highly
complex tissues using the implicit mesh-based Monte Carlo algorithm," Biomed.
Opt. Express, 12(1), 147-161 (2021).
https://opg.optica.org/boe/fulltext.cfm?uri=boe-12-1-147
- Ruoyang Yao, Xavier Intes, Qianqian Fang, "Generalized mesh-based Monte Carlo for
wide-field illumination and detection via mesh retessellation," Biomed. Opt.
Express, 7(1), 171-184 (2016).
https://opg.optica.org/boe/abstract.cfm?uri=boe-7-1-171
- Yuxuan Zhang and Qianqian Fang, "BlenderPhotonics - an integrated open-source
software environment for three-dimensional meshing and photon simulations in
complex tissues," J. Biomed. Opt. 27(8), 083014 (2022).
https://doi.org/10.1117/1.JBO.27.8.083014
Applications: MCX has been used to model light in the brain (fNIRS and
photobiomodulation dosimetry), neonatal heads, bone, and many other tissues; see
https://mcx.space/#publications for examples, and Google Scholar for the many
works citing MCX.
Acknowledgement
The MCX project, including MCX-CL, is funded by the US National Institutes of Health (NIH) / National
Institute of General Medical Sciences (NIGMS) under grant R01-GM114365. We thank
the users and contributors who reported bugs, suggested features and contributed
code; they are credited in the What's New section and the ChangeLog.
MCX-CL contains modified versions of the below source codes from other open-source
projects (with a compatible license).
cJSON library by Dave Gamble
ZMat data compression unit
miniz data compression library
- Files: src/zmat/miniz/*
- Copyright: 2010-2014 Rich Geldreich and Tenacious Software LLC, 2013-2014 RAD Game Tools and Valve Software
- URL: https://github.com/richgel999/miniz
- License: MIT License, src/zmat/miniz/LICENSE
LZ4 data compression library
LZMA SDK and Easylzma data compression libraries
- Files: src/zmat/easylzma/*
- Copyright: 2009, Lloyd Hilaiel, 2008, Igor Pavlov
- License: public-domain
- Comment:
All the cruft you find here is public domain. You don't have to
credit anyone to use this code, but my personal request is that you mention
Igor Pavlov for his hard, high quality work.
UBJ (UBJSON) library by Steven Braeger
Pybind11
myslicer toolbox by Anders Brun
