2.1.

Running Your First Simulation on Docker

After creating a Gitlab user account (replace user.name below with your account name), you can login to the registry and run the image by typing: docker login -u user.name registry.gitlab.com/fullmonte/fullmontesw/fullmonte-run docker run --pull --rm -ti -v /host/path:/output registry.gitlab.com/fullmonte/fullmontesw/fullmonte-run

Note that only a single login command is required per session. The options provided cause Docker to pull an updated image from the registry (--pull) when available, delete the container when it is finished running (--rm), make a host path accessible to the container (/host/path, which should be replaced by the location you want the output data to appear; data written to /output in the container will appear on the host at that path), provide a terminal (-t) and run interactively (-i). Within that container, you can run a simulation of a point source located inside the bladder by invoking the simulator inside the container and giving it a script name to run: cd /output tclmonte.sh /usr/local/FullMonteSW/share/Examples/BladderDirectional/sim.tcl

Simulation time will depend on machine performance but may be a few tens of seconds. The wiki page “directed surface scoring” gives a detailed walkthrough of the simulation results and visualization.

2.2.

Problem Definition

Inputs which define the problem—geometry, optical properties, and source configuration—are the key arguments to the program that change the expected result. The FullMonte kernel itself outputs its results in terms of packet weight. If $N$ packets are launched, then the total energy absorbed plus exited from the system will equal $N$ arbitrary units. To model total energy emission from the source $E_0$ in physical units (e.g., Joules), the results must be normalized by a factor of $E_0/N$, which is a task left to the user.

Technically, packet count then alters the expected output, but only by a proportionality constant that should be removed in postprocessing. Its most important effects are on result variance and run time, so we class it as a runtime option instead of part of the problem definition.

The units of length must be consistent between the mesh and the optical properties, which are given in inverse-length units, and results will be reported in the same units. For instance, if a mesh is given in mm, then the optical properties will be interpreted as ${\mathrm{mm}}^{-1}$ and after applying the EnergyToFluence filter on the data, results will be in $\frac{\text{packet weight}}{{\mathrm{mm}}^2}$ (which can be rescaled to ${\mathrm{Jmm}}^{-2}$ or ${\mathrm{Wmm}}^{-2}$ as desired). A mismatch between the mesh unit and optical-property unit will cause the simulator to use incorrect optical properties.

2.2.1.

Geometry

Delineation of medical images and conversion of the results into a volumetric tetrahedral mesh ( www.cgal.org) are complex tasks but both commercial and open-source (segmentation: ITK-Snap www.itksnap.org, meshing: CGAL www.cgal.org) tools exist, so we provide only file import-export facility in FullMonte. While there are many file formats, tetrahedral meshes are fundamentally composed of points and tetrahedra, with each tetrahedron being assigned to a region. Points are specified in a list with $(x,y,z)$ coordinates for each. Tetrahedra are given by the index in the point list for each of its four vertices and an associated material code. FullMonte can read and write the mesh formats used by other simulators (TIM-OS, MMC), the COMSOL multiphysics package, and the open-source Visualization Toolkit (VTK). The Colin27 brain mesh produced by Fang⁶ and distributed with MMC is used below as an example (Fig. 1) with over 400,000 tetrahedra.

Fig. 1

Cutway view of Fang’s Colin27⁶ mesh used for testing, a good example of adaptive mesh sizing and representation of curved surfaces.

When asked to handle MCML layered geometries, each layer is represented as a degenerate tetrahedron with two parallel faces in the $xy$-plane (the top and bottom layer bounds) and two faces at infinity. Layered cases may therefore be run through the unmodified tetrahedral kernel.

2.3.

Runtime Options

By runtime options, we mean settings which impact the behavior of the simulator in ways other than the first moment of the output distribution. The expected results are invariant to the run parameters, but the result variance, the distribution of variance across different output quantities, and the computational run time may change significantly.

The number of packets to simulate ($N_0$) is the most important such quantity, since run time is proportional to $N_0$ while output quantity uncertainty (standard deviation) decreases in proportion to ${\sqrt{N}}_0$. Users will need to assess qualitatively or quantitatively what level of accuracy is acceptable given the runtime cost. We opt for ${10}^6$ as a default value that often yields a fluence map of subjectively reasonable quality within seconds to tens of seconds for the BLI and PDT cases tested so far. The packet count required to achieve a given level of quality depends on a complex interaction of mesh fineness, optical properties, the output quantity of interest, tolerable level of variance, and region of interest. Unfortunately, a general method for determining an acceptable packet count for a given application does not exist so the user must rely on judgment and experiment. When examining linear combinations of elements, or when the output quantity otherwise contains some variance-reducing tendency (e.g., dose-volume histograms for PDT),¹⁴ then far fewer ($10\times$ or more) packets may be required with correspondingly faster run time. For high spatial resolution, low desired output variance, and low-fluence/low-interaction regions of interest, orders of magnitude more packets may be called for to get acceptable output statistics. The user can assess result variance by running multiple times with different random seeds. FullMonte provides both command-line and scripting methods to generate per-seed results or a mean–variance summary.

For finer control of time-quality tradeoffs, the roulette threshold¹⁵ (weight threshold below which packets are subjected to roulette elimination) modifies the run time-variance tradeoff in elements with very low fluence. Since roulette bundles a number of low-weight packets into a single higher-weight packet, increasing the threshold will speed up the simulation at the cost of fewer but higher-weighted absorption events (i.e., higher result variance) in low-fluence regions. MMC provides a command-line option with a default value of ${10}^{-6}$, whereas TIM-OS has an unnamed constant fixed at compile time (TIMOS.cpp line 1782), FullMonte has multiple ways to set the value which defaults to ${10}^{-5}$, and MCML is fixed at compile-time to ${10}^{-4}$ (MCML.h line 74). To compare simulator performance and quality fairly, this value should be set equally for all parts of the comparison. We use FullMonte’s value of ${10}^{-5}$ when comparing with MMC and TIM-OS, and ${10}^{-4}$ to match MCML.

2.4.

Simulation Kernel

FullMonte’s kernel is designed so that the same photon-propagation logic can be used to capture a number of different lifecycle events in various and configurable ways detailed in Sec. 3. The most common case—scoring surface exit and volume absorption events—is the one provided with the standard command-line interface. It outputs the total, non-normalized (nonunitized) packet weight exiting through each surface face and deposited in each volume element. These choices reflect a design preference for maximal simplicity so that the output from the kernel is exactly what was recorded in the simulation, with all normalization, unit conversion, and changes of physical quantity (e.g., energy to fluence) being deferred to a postprocessing stage. Conservation of energy can be checked by summing the total surface and volume scores, and comparing against the total launched.

When cross-verifying FullMonte, we found and reported a bug in MMC that appears under certain combinations of options (Havel raytracing, energy output, piecewise-constant scoring basis -O E -C 0 -M H) that is a direct artifact of a more-complex control path in which the control flow through photon tracing depends on multiple command-line options: ray-tracing method, output type, and output coordinate basis. Our election to keep photon propagation, event scoring, and output postprocessing code fully separate with well-defined interfaces reduces the probability of such errors.

2.5.

Output Data

When the result is normalized by the number of packets (a factor of $1/N_0$), it gives a Green’s function solution in units of dimensionless weight. Packet weight is analogous to photon arrival probability (directly proportional to energy via $E_{\gamma }=h\nu$), so multiplication of the Green’s function by the appropriate number of physical photons emitted $N_{\gamma }$, or the total energy $E_0$ emitted gives the result as a physical quantity. The results may also be interpreted per unit time as photon or radiant flux.

FullMonte also has the ability to score the energy crossing internal boundaries, which can be useful for modeling thin layers, such as the urothelium for PDT of the bladder. Meshing such structures directly is computationally expensive, may produce meshes with many small tetrahedra that degrade quality (in terms of output variance) and numerical stability. Separate logging of photon transit through the layer is preferable because it avoids these issues. A detailed example is provided with the source distribution and Docker image (Examples/BladderDirectional).

2.6.

Postprocessing

From the raw exit and absorption scores ($s$) and a specification of total energy emitted ($E_0$), FullMonte can produce a number of different outputs. Per volume element, it can provide physical energy absorbed by rescaling weight ($E=s/E_0$), average energy density ($E/V$), or average fluence ($\mathrm{\varphi }=E/V{\mu }_a$). Over the surface, typically, it is the emittance (exit energy per unit area) that is of interest. Other options include dose–volume and dose–surface histograms for regions partitioned by tissue type.

2.7.

Visualization

Visualization is made possible using VTK. Basic Tcl script users can export .vtk files for viewing in Paraview, as shown in the examples. More advanced users can automate generation of figures using VTK’s own Tcl scripting interface from the files generated. VTK and Paraview provide their own extensive documentation and tutorials showing how to use them.

2.8.

Result Export

We provide simple text-format export facilities so that the mesh and associated per-element data can be saved and read into other postprocessing tools, such as MATLAB. When running multiple simulations of the same problem, we define a mean–variance file format for future comparison. Data from TIM-OS, MMC, and MCML can also be imported for comparison or conversion purposes.

3.1.

MCML Support

For verification, we modified FullMonte to support MCML-like geometry definitions and scoring since it is widely used and trusted. Handling the problem description required creation of a new reader class to handle the.mci file format, extracting the material definitions, layer information, and scoring grid information. We created a facility to convert the layered geometry into a degenerate tetrahedral mesh in which the top and bottom faces of each tetrahedron are parallel in the $xy$-plane, and the other two are at infinity. Since FullMonte typically scores absorption per tetrahedral element, we created a new Scorer class that accepts absorption events and accumulates them to cylindrical bins defined by $(r,z)$ as MCML does. Last, the conversion from output weight score to fluence used the existing EnergyToFluence class and required only that the newly created layered geometry description provide the volume and material identifier for each bin. The core photon-tracing logic is completely unmodified compared to the normal 3-D tetrahedral mode, meaning the change set is confined to a small amount of new code and no modifications to the propagation core.

3.2.

Path Tracing

To visualize the operation of the simulator for debugging and educational purposes, we produced a version of the simulator to capture the full position trace of each photon, producing Fig. 2. We added a PathScorer that concatenates the positions from event notifications of launch, reflection, refraction, and scattering. To invoke that scorer, we created a TetraTraceKernel that binds the standard core loop to the PathScorer, without altering the core photon-propagation logic. After postprocessing to reformat the data, the result is saved as a .vtk file for visualization.

Fig. 2

A BLI testcase using an embedded source with 1k packet traces simulated by TetraTraceKernel, rendered by Paraview (script file mouse_paths.tcl, state file MouseTrace.pvsm)

3.3.

New Source Types

Physical light sources are represented by two classes in FullMonte: a class which specifies the logical description of the optode and an Emitter-derived class which has the machinery to use a random-number generator to emit photons matching that description into a specific mesh. We provide an EmitterFactory class that maps logical source descriptions to concrete Emitter derived classes. When launching photons, we use C++ virtual function mechanism to dispatch a launch request from the kernel, which knows only that it has a pointer to some kind of Emitter, to a specific derived class which handles the details. Since photon launch is infrequent compared to ray-tetra intersection, scattering, and other operations, the overhead imposed is negligible.

4.1.

Internal Consistency Checks

In addition to the primary outputs of exit and absorption photon scores, the simulator also produces other quantities of interest to assess performance and check conservation of energy. All provided FullMonte kernels return a result specifying the disposition of weight between total absorption, roulette increase/decrease, and exit. The values can be compared against the totals launched to verify that all energy is accounted for (Table 1 summarizes lifecycle events), where a conservation failure would indicate either a problem in the scoring logic, an unaccounted loss of weight within the kernel, or a numerical problem in weight accumulation.

Table 1

Events impacting conservation of energy in the simulation kernel, where w0 is the packet weight before the event, P is probability of roulette win, and μa,μs are the optical absorption/scattering coefficients.

MMC invests some extra computational effort in performing double-precision Kahan summation for each mesh element score to guard against underflow. We opted instead for higher-performance simple summation at double precision and check ex post if there are problems arising from that choice in the form of lost weight, an event which has not yet arisen. The total over all elements for absorption and exit can be checked against the sum of the per-element scores as well to identify whether numerical underflow is a significant issue. Such checks would also catch catastrophic failures of simulation logic in which packets are lost, which is primarily useful for developers to quickly identify logical errors or numerical precision issues when making changes in the core photon propagation loop (which ought to be infrequent, per Sec. 3).

Checks take the form of both rigid conservation laws which must be true to within rounding error and statistical checks which are zero-mean with an associated variance. Physically, photons entering the tissue will either exit or be absorbed, so the physical energy launched should equal the energy absorbed plus exiting. In simulation, roulette is a random process introduced as a computational shortcut that preserves expected conservation of weight (energy), so weight added via roulette wins less that dropped by roulette termination should equal zero within statistical variance. The results should account for all energy so the total weight sourced (by launch and roulette) and disposed should be equal to within rounding error; any imbalance in excess of the roulette variance can be attributed to numerical precision effects.

Regardless of the scoring system (spatial binning) employed, the total of each absorption or exit score should also match that shown in the conservation scorer to within rounding error. Intuitively, the total weight of all photon absorption or exit events ought to equal the sum over all volume and surface elements.

4.2.

Unit Testing

The source code provided for FullMonte includes a number of unit tests verifying key functions, such as the distribution of photons emitted by the various Emitter classes, the distributions of the random-number generators, and geometric primitives, such as ray-tetrahedron intersection tests crucial to the success of the kernel. After compiling the code, a complete set of unit tests can be run by the command ctest -R unit_ to provide a summary with the number of tests run and identification of any failing tests. Provision of unit tests also provides a quick check that any new feature or performance enhancement has not caused the existing code to become inaccurate, so that developers may work to improve performance or features but remain confident that key specifications are still met.

4.3.

Cross-Verification with Other Simulators

In the physical system being simulated, the radiance $L(\mathbf{r},\widehat{\mathbf{n}},t)$ as a function of position $\mathbf{r}$, direction $\widehat{\mathbf{n}}$, and time $t$ must be a solution to the radiative transfer equation. The MC scoring system traces photon packets through the geometry, reporting photon lifecycle events (e.g., launch, arrival at an element boundary, refraction, reflection, absorption, scattering, domain exit) to Scorer classes, which accumulate scores to approximate physical quantities of interest. Typical scoring arrangements include the energy absorbed within a volume or passing through a surface, though they can be further refined by subdividing into slices of time, angular resolution, etc. We consider the statistical distribution of a single packet’s contribution $\mathbf{S}$ to the output score vector in terms of its mean $\mu$ and variance ${\sigma }^2$. In any correct MC formulation, the mean $\mu$ must be proportional to the physical solution dictated by the radiative transfer equation. On the other hand, the per-packet variance and its distribution between the output elements depend on the exact MC formulation and variance-reduction schemes (most notably here, roulette form and parameters).

To conduct a simulation, we generate $N$ independent photon paths whose contributions to the score vector $\mathbf{s}$ are drawn from $\mathbf{S}$ and sum their contributions to produce an output score $\mathbf{x}$ with elements $x_j={\sum }_{i=1}^N{s}_{ij}$. Its expected value and variance are given by

When using an MC simulator, variance is typically a secondary consideration as users are concerned only to set a packet count $N$ to make the variance negligibly small for their purposes. To robustly compare two datasets, though, the tolerance permitted needs to account for the expected variance to assess whether a difference is statistically significant, or can plausibly be dismissed as random variation.

We pose equality comparison of two independent output samples $\mathbf{x},\mathbf{y}$ as a test of the null hypothesis that they are drawn from independent MC simulations of the same problem, which necessarily means independent identically distributed packet score contributions described by $\mu ,\sigma$. Considering a measure of the difference $\mathrm{\Delta }$ with different numbers of packets simulated $N_x,N_y$ yields

While the difference simulating identical systems is expected to be zero, it will, in practice, have some variance, which must be distinguished from systematic difference. Previous work has relied on subjective comparison, imposition of arbitrary tolerances, excluding low-value (which generally means high-variance) elements, or running very long large packet counts to drive the variance to nearly zero. By considering the MC result variance, we arrive at a more rigorous standard in which the simulated data itself is used to derive an appropriate comparison tolerance. Under the null hypothesis, the per-packet score contribution parameters are the same, so by substitution of Eq. (2), we know the variance of the difference vector:

Since the actual score contribution variance ${\sigma }^2$ is not known, it must be estimated. Let us consider the comparison of a reference arm $\mathbf{x}$ generated by simulating $N_r$ independent runs of $N_p$ packets to provide empirical estimates of variance via Eq. (2). Though the score contribution vector is not normally distributed, a sufficiently large summation of such event contributions will converge toward normal. With that estimate in hand, we may proceed to test whether an observed difference $\mathrm{\Delta }$ falls within the distribution proposed by the null hypothesis.

To assess whether the simulator in the test arm produces significantly different results, we use it to generate a single simulation output $\mathbf{y}$, possibly using a different number of packets. Using the variance estimated from the reference arm, we assume that the variance of the test arm is the same, scaled to the number of photons run by using Eq. (2). Since the convergence relied upon above for estimating $\sigma$ depends on the frequency of accumulation events in the score vector, different vector elements will require different numbers of packets to converge to a given level. We can speed up the test by considering only a subset of elements $\mathcal{C}$ that converge at a relatively small number of packets for comparison. Taking elements with a coefficient of variation ($\sigma /\mu$) that is under 10% in both the reference and test arms has worked well as a guideline.

Under those assumptions, we produce a vector of $z$-scores for the difference $\mathbf{\Delta }$ between the reference and test results:

The sum of squares of $k=|\mathcal{C}|$ such variables is ${\chi }^2$ distributed with $k$ degrees of freedom, leading to a standard statistical test of the null hypothesis with test statistic:

If the observed vector is sufficiently improbable given the variance established in the reference arm, then we conclude that the test and reference arms model different phenomena.

5.1.

Subjective Comparison

In all test cases, the difference appears subjectively zero-mean, with the magnitude increasing away from the source as elements receive fewer scoring events, and hence, have higher variance due to a lower frequency of packet arrival and hence smaller sample. Spatially correlated structure would indicate an underlying difference in the phenomena simulated, but we note no such structure in any of the comparisons. The figures below illustrate a 3-D case (Colin27-Point0 MMC vs. FullMonte) in Fig. 3 and a 2-D case in Fig. 4 with MCML.

Fig. 3

Colin27-Point0 case showing a cut view through a 3-D model. Samples with a coefficient of variation exceeding 20% have been excluded. The red sphere shows the point source location. From left to right: (a) fluence, (b) MMC/FM result ratio, and (c) difference $z$-score.

Fig. 4

Comparison of FullMonte versus MCML showing from left to right: (a) fluence, (b) percentage difference, and (c) $z$-score difference.

5.2.

Delta Statistics

Since our comparison test relies on assumptions about the statistical distribution of the difference between simulation runs, we must first validate those assumptions. If the assumptions and the null hypothesis are true, then the scaled difference vector $z$ of Eq. (8) will have a unit normal distribution. We checked this assumption using a Q-Q (quantile–quantile) plot, which shows each measured data point’s measured value on the $y$ axis against the value of the test distribution at the same quantile on the $x$. A unit-slope straight line through the origin indicates a perfect match. Figure 5 shows the result from the MCML-Sample0 case was the difference vector, which is indeed normally distributed with zero mean. The variance very slightly exceeds expected, which suggests that it may be slightly underestimated, likely due to low-fluence regions that have not quite converged to a normal distribution. Results for the other test cases were similar, indicating that the model is correct for the simulators, test cases, and parameters (run counts and packets per run) posed, so the results are well described by a normal distribution.

Fig. 5

A Q–Q plot of the scaled MCML-FullMonte difference vector ($z$-score) against the standard normal distribution, showing strong agreement for the MCML-Sample0 regression test.

5.3.

${\chi }^2$ Test Validation

Since the test is statistical, we need to check that our test is suitably sensitive and specific.

5.3.2.

Sensitivity

To check that we can indeed detect correctly when the null hypothesis is false, we altered the optical parameters subtly and ran another 64 simulations at ${10}^6$ packets to compare against the reference arm. With Colin27-Point0 and an increase of 1% in scattering in the gray matter, the test flagged a difference in absorbed energy 27/64 times. The distribution of ${\chi }^2$ statistics from the test runs is presented in Fig. 6, with a vertical line showing the $p=0.05$ critical value. Values to the left are accepted as identical to the reference within variance at $p=0.05$ while those to the right are significantly different. We conclude that the test as configured is acceptably sensitive to changes in optical properties, despite its modest run time (seconds to tens of seconds). Choosing longer simulations (larger packet counts) will leave the mean unchanged while reducing variance, which would increase sensitivity even further if needed.

Fig. 6

The empirical cumulative distribution (CDF) of observed ${\chi }^2$ values for the test covering 95% of absorbed energy in Colin27-Point0. The critical value ${\chi }_c^2=8911.91$ is shown as a vertical line and was computed for $p=0.05$ using the normal approximation for 8694 degrees of freedom (the number of tetras included).