2.1.

Similarities to mcxyz

MCmatlab, in the same way as mcxyz, operates in Cartesian coordinates and makes no assumptions as to symmetry in the simulated volume. Its boundary is a rectangular cuboid and its internal volume is uniformly divided into voxels, which are themselves also rectangular cuboids. The user assigns to each voxel a medium or tissue type, described by its absorption coefficient ${\mu }_{\mathrm{a}}$, its scattering coefficient ${\mu }_{\mathrm{s}}$, and its Henyey–Greenstein scattering anisotropy factor $g$ at the applied wavelength.

An input light beam is simulated by launching photon packets and calculating their paths in the simulated volume, using pseudorandomly generated numbers to determine initial photon packet position and trajectory as well as path lengths between scattering events and scattering angles. As photon packets propagate from one voxel to the next, they deposit some of their energy (“weight”) into the voxel depending on the voxel’s absorption coefficient ${\mu }_{\mathrm{a}}$. This deposited energy is numerically accumulated in a 3-D matrix, which, upon conclusion of the simulation is normalized to the input power, providing the local absorption rate per watt of incident light, in units of $\mathrm{W}/{\mathrm{cm}}^3/\mathrm{W}$.incident. Dividing this by the voxels’ absorption coefficients yields $F(x,y,z)$, the local fluence rate (irradiance or intensity) per watt of incident light, in units of $\mathrm{W}/{\mathrm{cm}}^2/\mathrm{W}$.incident.

2.2.

Differences to mcxyz

There are three main differences between MCmatlab and mcxyz (as of mcxyz version June 1, 2017), mostly centered around improved usability for nonexperts in programming. The most apparent of these is the rewriting of the C code to be compilable and callable in Windows or on Mac from MATLAB as a MEX function, combining the speed of C with the flexibility and ease-of-use of MATLAB. This means that one can perform the entire simulation initialization, execution, postprocessing, visualization, and further interfacing in MATLAB. Therefore, users on Windows or Mac no longer need, e.g., a separate UNIX-emulating console shell to run a compiler and the MC executable, nor is there a need for any intermediate input/output files, since all data are passed in and out of the MEX function through memory. The simulation inputs and the fluence rate matrix output (the solution to the radiative transport equation) are saved as MATLAB “.mat” files, a compressed binary file format allowing for flexible and easy saving and loading and ensuring smaller file sizes compared to uncompressed binary or plain-text file formats.

The second main difference to the mcxyz code is the introduction of a wider variety of beam type definitions. The full selection of beam types is currently pencil beams; isotropically emitting points; infinite plane waves; Gaussian focus, Gaussian far-field beams; Gaussian focus, top-hat far-field beams; top-hat focus, Gaussian far-field beams; top-hat focus, top-hat far-field beams; and Laguerre–Gaussian LG01 beams.

The focal plane of the beams is defined as the plane containing the user-specified focus point with a normal vector given by the user-specified beam center trajectory. For the beams with Gaussian and/or top-hat beam profiles, the initial position and trajectory of each photon are calculated by first randomly sampling, according to the chosen focus beam profile, Gaussian, or top-hat, the point in the focal plane that the photon would hit in the absence of scattering and then randomly sampling the initial trajectory according to the chosen far-field beam profile. The initial position of the photon is then defined as the intersection between the $z=0$ plane and the line going through the focal plane target point along the photon’s initial trajectory. For the Laguerre–Gaussian LG01 beam, the sampling is slightly different; for a focal plane target point at a position $r$ relative to the focus point, the trajectory is chosen only among those directions that are orthogonal to $r$. This ensures that the intensity is zero everywhere along the center axis of the beam, as expected of an LG01 beam. The focus and far-field widths can be freely specified and could, for instance, be chosen to match that of a diffraction-limited beam.

The third major difference is the visualizations. It can be challenging to properly visualize data in a 3-D volume, especially if the volume contains fine-structured heterogeneities in multiple layers. Therefore, MCmatlab uses interactive 3-D slice plotting, in which the color scale can be switched between linear and logarithmic by checking a box.

Additional changes from mcxyz include: (a) changing pseudorandom number generator from Donald Knuth’s subtractive algorithm from 1969²⁰ to the Fast Mersenne Twister algorithm from 2008,²¹ (b) a change of coordinate system from $yxz$-coordinates to $xyz$-coordinates, (c) CPU multithread parallelization of the simulation on Windows using OpenMP, (d) minor bug fixes, and (e) significant general speed optimizations.

MCmatlab can also optionally simulate refraction and Fresnel reflection at interfaces between media with different refractive index, provided that the geometry is oriented so that the interface is parallel to the $xy$-plane. Since MCmatlab does not track polarization, the reflection probabilities are calculated assuming that the light is unpolarized.

2.3.

Example and Comparison of Results

To test whether the results of our MC RTE solver agrees with prior work, we compared MCmatlab with mcxyz using the example shown on the mcxyz website,⁴ a $200\times 200\times 200$ voxel model of a cylindrical blood vessel of radius $100\mu \mathrm{m}$ embedded $400\mu \mathrm{m}$ deep in dermis, with a $60\text{-}\mu \mathrm{m}$ layer of epidermis on the surface. Figure 1 shows the geometry, visualized using the interactive 3-D volumetric slice plotting used for many different plots in MCmatlab. Table 1 shows the parameters used for the demonstration. For comparison purposes, all ${\mu }_{\mathrm{a}}$, ${\mu }_{\mathrm{s}}$, and $g$ values were set equal to those used in the mcxyz demonstration example. The thermal properties are used only for the thermal simulations, as shown in Sec. 3. Note that the values are not necessarily accurate for real tissue but are for demonstration purposes only. The illumination beam was a collimated top-hat beam of radius $300\mu \mathrm{m}$.

Fig. 1

Screenshot of MCmatlab’s illustration of the simulation geometry definition with slices shown at $x=0.05\mathrm{cm}$, $x=0\mathrm{cm}$, $y=0.05\mathrm{cm}$, $y=0.01\mathrm{cm}$, and $z=0.1\mathrm{cm}$.

Table 1

An overview of the optical properties, the physical thermal properties, and the chemical thermal properties of the four tissues/media used in the example. The values are for demonstration purposes only.

For this comparison, $6.36\times {10}^7$ photon packets were launched in both programs. In Fig. 2, the results of the MCmatlab simulations are shown in terms of the fluence rate and the absorbed power per unit volume. Figure 3 shows a contour-plot comparison to the results of mcxyz in the $y=0$ and $x=0$ planes. It can be seen from Fig. 3 that there is good agreement between the results of mcxyz and MCmatlab. There is a bug in mcxyz that affects fluence rates in voxels bordering some of the cuboid boundaries, but this does not affect all inner voxels.

Fig. 2

Screenshots of MCmatlab’s illustration of (a) the fluence rate and (b) the logarithm of the absorbed power per unit volume. Slices are shown at the same positions as in Fig. 1.

Fig. 3

A comparison of the values of the logarithm of the fluence rate calculated by MCmatlab and mcxyz after launching $6.36\times {10}^7$ photon packets. The geometry definition is shown in Fig 1. Black contour curves: MCmatlab, orange: mcxyz. (a) $x=0$ and (b) $y=0$.

Running in a UNIX-emulating Cygwin terminal on a mid-range laptop PC from 2012 (Dell Latitude E6530) with 64-bit Windows, mcxyz launched $6.29\times {10}^4$ photons per minute. When run in MATLAB on the same PC, MCmatlab launched $1.09\times {10}^6$ photon packets per minute (a factor of 17 more). MCmatlab’s multithreading to the laptop’s 8 CPU processors accounts for roughly a factor of 4.5 of the speed increase compared to mcxyz, whereas the remaining factor of 4 is due partly to switching to the Fast Mersenne Twister pseudorandom number generator and partly to many small optimizations throughout the code base.

3.1.

Theory

The heat deposition and diffusion solver calculates the time-resolved temperature distribution in the simulation volume based on the fluence rate $F$ and the geometry definition of the MC simulation, this time with additional material parameters; the media’s physical thermal properties $c_V$ and $k$ and, optionally, chemical thermal properties $E_a$ and $A$ of the considered chemical reaction. Although the physical thermal parameters must be specified for all involved media, the chemical properties may be specified for all, some, or none of them.

The temperature evolution can be simulated for continuous light exposure or for single- or multipulse light exposure with user-specified duty cycle and period. At the end of the simulation, a plot is shown that illustrates the spatial distribution of chemically changed media in the simulation volume. A video can also be automatically generated showing the temperature evolution.

Denoting the local volumetric heat capacity by $c_V$ and thermal conductivity by $k$, the heat equation in continuous form is given as

The thermal chemical change is quantified for each voxel in terms of the dimensionless Arrhenius $\mathrm{\Omega }$ value, which starts at zero and increases over time as the voxel is exposed to elevated temperatures $T$, depending on the Arrhenius activation energy $E_{\mathrm{a}}$ and Arrhenius pre-exponential factor $A$:

3.2.

Example

In this demonstration, the blood vessel defined in Sec. 2 was used as an input to the thermal solver, setting the input power to 4 W, running for 5 light pulses of 1 ms on-time and 4 ms off-time each. Video 1 is provided as supplementary material, and a still image is shown here as Fig. 4. The resulting illustration of damaged tissue is shown in Fig. 5.

Fig. 4

Still image from a video showing the temperature evolution during five pulses of light incident on the tissue defined in Sec. 2 (Video 1, MP4, 360 KB [URL: https://doi.org/10.1117/1.JBO.23.12.121622.1]).

Fig. 5

Screenshot of MCmatlab’s thermal damage illustration. The volume of coagulated blood is shown in green.