2.1.1.

General structure of the open-source package PyXOpto

The general outline of the package PyXOpto is shown in Fig. 1. While the light propagation core is based on the PyOpenCL library and written in OpenCL C extension, the supporting modules are entirely Python-based and form a highly modular structure of the package PyXOpto. The PyOpenCL library enables massively parallel computations on OpenCL-enabled devices such as AMD, Intel, Nvidia, Qualcomm, ARM, and more. The light propagation core offers customizable kernel types that can perform single- or double-precision computations and use 32- or 64-bit photon packets counters. This allows more than the usual $2^{32}$ ($4.29\times {10}^9$) photon packets to be processed in a single run.

Fig. 1

Schematic of the modular structure of the open-source package PyXOpto. Boxes represent available modules within the subpackages xopto.mcml and xopto.mcvox or helper utilities.

The package PyXOpto is divided into two subpackages for modeling the light propagation in multilayered (xopto.mcml) and voxel-based (xopto.mcvox) configurations of tissues. Both xopto.mcml and xopto.mcvox share the same principle by using a central module mc from which the constructor Mc can be passed several modules that define the geometry, phase functions, surface layout of the top or bottom tissue boundary, sources, detectors, fluence accumulators, and trace accumulators. This configuration makes the PyXOpto package highly flexible since various combinations of modules can be used interchangeably and seamlessly without modifying the central light propagation core. In addition, the package PyXOpto offers various utilities for easier calculations of SVs, a transformation of reflectance to a spatial-frequency domain, scattering phase functions, and calculations of their subdiffusive quantifiers^32–34 (e.g., Legendre moments, $\gamma$, $\delta$, and $\sigma$) and various materials (e.g., wavelength-dependent refractive indices). The utilities also include a module for computing and accessing the presented MCDataset or generating new datasets for custom source, detector, and tissue configurations.

Both xopto.mcml and xopto.mcvox subpackages follow the standard propagation logic of the photon packets as used in the MCML.¹⁶ The details regarding each of the previously mentioned modules to the Mc constructor are described in the sections below.

2.1.2.

Geometry

The geometry description depends on the selection of the PyXOpto subpackages.

For the xopto.mcml, layers are defined in a similar fashion to the MCML except for the phase function that is further detailed in Sec. 2.1.3. For each layer, the user defines optical properties including the absorption and scattering coefficients, refractive index, the thickness of the layer, and phase function. The topmost and bottommost layers correspond to the media surrounding the tissue and supply only the refractive index mismatch, while other layers are used for photon packet propagation.

For the xopto.mcvox, the three-dimensional (3D) volume is labeled using a 3D matrix and for each label, the optical properties are defined. Generally, each voxel can be assigned a different absorption and scattering coefficient, refractive index, and phase function.

2.1.3.

Phase function

One of the stronger suits of the proposed light propagation model is that we do not strictly enforce the generic Henyey–Greenstein (HG) phase function.³⁵ While most of the currently available free and open-source implementations of the MC method indeed offer only sampling of the scattering angles according to the HG phase function, such a restriction can significantly limit the applicability of the light propagation model for tissues exhibiting microscopic refractive index variations that cannot be accounted for by a single-parameter HG phase function. We have noted in several of our previous publications that the phase functions significantly influence the acquired reflectance at small source–detector separations when using optical fiber probes.^33,36 Therefore, we have implemented a general lookup-table-based sampling scheme for phase function with cumulative probability distributions that have no analytical inverse.³⁷ Such an approach supports the utilization of various scattering phase functions, such as the modified versions of the HG,³² Gegenbauer kernel,³⁸ and Power of Cosines,³² measured phase functions, and more complex scattering phase functions based on the Mie theory.³⁹ The latter can be used to simulate light propagation through water-based suspension of spherical microparticles, such as polystyrene microspheres and is therefore ideally suited for the fabrication of optical phantoms.

The phase function is used as an input to the geometry module that describes the biological tissue in which the light propagation is modeled. In general, for each layer or voxel, a different phase function can be applied thus enabling various complex configurations of tissues with spatially varying microstructure.

2.1.4.

Surface layout

The surface layout is an optional module and offers a realistic description of the interface between the tissue and the surrounding medium such as a Lambertian surface or a realistic layout of the optical fibers and materials comprising the tip of common commercially available optical fiber probes. We have found that especially the optical probe tip can significantly affect the acquired reflectance and thus the probe-tissue interface must be precisely described.¹⁴ To the best of our knowledge, the existing MC method implementations lack support for this important feature. If the surface layout is not set, the simulations are run with laterally uniform boundaries, which is the “standard” way of performing MC simulations in multilayered tissues.

2.1.5.

Source

The source is a mandatory input to the Mc constructor. We have implemented an isotropic source, line source, collimated Gaussian and uniform beams, optical fiber source, and rectangular sources that can be patched together to form more complex sources. All the sources require their specific parameters. For example, the optical fiber source depends on the diameter, position, direction, refractive index, and numerical aperture. Furthermore, the optical fiber and rectangular source can both accept a custom angular emittance profile measured with a goniometric setup. More realistic modeling conditions, especially custom emittance profiles, are useful for, e.g., health monitoring devices.

2.1.6.

Detector and fluence

Detection of photon packets can be performed using a surface detector or a 3D accumulator which we refer to as fluence.

Implemented surface detector can be divided into groups of specular, top, and bottom. Specular detectors merely record the fractions of the photon packets that have reflected due to the surrounding medium-tissue or source-tissue refractive index mismatch. Top and bottom sample surface allow different configurations of the detectors which can be selected from total detectors for quantities such as reflectance and transmittance, radial detectors for spatially resolved reflectance and transmittance, Cartesian detectors, detector configurations related to commercially available optical fiber probes such as the six-around-one and linear-arrays and finally a detector based on custom arrays of optical fibers for special use cases. As in the case of the source definitions, detectors accept various parameters including minimum acceptance cosines in the case of radial and Cartesian detectors and numerical aperture for fiber-based detectors. Finally, detectors can be specified with a custom angular acceptance measured with a goniometric setup.

The absorbed weight of the photon packets can also be recorded volumetrically using the so-called fluence accumulator. The latter can record energy deposition or true fluence depending on the additional flag that is specified to the constructor.

2.1.7.

Trace and sampling volume

The photon packets can optionally be traced by recording the photon packet position, direction, and weight at each absorption, scattering, or refraction event. The traces can be recorded independently of the selected detector and easily filtered according to the desired final direction and position of the photon packets. The filtered traces can be transformed into a voxelated SV, where each voxel contains path lengths of the photon packets multiplied by their respective terminal weights. The SV can be used to represent a 3D heat map indicating the probability of a photon packet to travel through a particular voxel from the source to the detector.⁴⁰ The processing of the photon packet traces in our implementation utilizes the PyOpenCL implementation that significantly accelerates the computations of the SVs.

2.2.1.

MCML comparison dataset

The MCML comparison dataset is used for comparison to the MCML by Wang et al.¹⁶ It comprises three subsets, namely one-layer-1-mm, one-layer-100-mm, and two-layer-100-$\mu m$-1-mm corresponding to different layer arrangements and thicknesses, and optical properties (see Fig. 2, geometry). The source is an infinitely narrow perpendicular line (pencil beam). The optical properties of the single or top layer are varied across all combinations given by the outer product of the anisotropy factor $g$, absorption coefficient ${\mu }_a$, and reduced scattering coefficient ${\mu }_s^{\prime }$ (other optical properties and parameters remain constant). The specular detector stores total specular reflectance relative to the initial number of photon packets $N$. The radially accumulated photon packets at the top and bottom surface of the sample hold the spatially resolved reflectance and transmittance provided relative to $N$ and per surface area of the annular ring (${1/\mathrm{m}}^2$). Note that the transmittance is different from zero only for the one-layer-1-mm and two-layer-100-$\mu m$-1-mm subsets. The photon packets are collected for all exit angles (numerical aperture equal to 1.0). Finally, the energy deposition is stored in the fluence accumulator and is given relative to $N$ and the volume of the voxel (${1/\mathrm{m}}^3$). The MCML comparison dataset comprises 81 different configurations. Each configuration was simulated with ${10}^8$ photon packets.

Fig. 2

Configuration of the MCML comparison dataset. The colored dots correspond to different layer arrangements, thicknesses, and optical properties (see geometry).

2.2.2.

MCML dataset

The MCML dataset is considerably more extensive than the MCML comparison dataset and includes various phase functions, surface layouts, sources, and detectors. The dataset comprises seven groups of subsets that are color-coded and presented in Fig. 3. Each color represents a different source type. The number of subsets in each group depends on the number of different parameters defined for each source type. For example, the group corresponding to the collimated Gaussian beam source includes a single subset with the FWHM of the beam equal to $100\mu \mathrm{m}$, while the group corresponding to the single-fiber source includes four subsets where the diameter of the fiber core and cladding varies from 100 and $120\mu \mathrm{m}$ to 800 and $820\mu \mathrm{m}$. In total, the MCML dataset comprises 11 subsets. The subsets were selected according to various source and detector configurations that we have found in the literature.^41–43 Each subset is assigned its corresponding geometry, surface layout, source, and detector according to the colored dots given in Fig. 3.

Fig. 3

Configuration of the MCML dataset. The colored dots represent subsets of the MCML dataset with specifically assigned geometry, surface layout, source, and detector types.

For example, the fiber dataset is based on a single semi-infinite layer geometry with the refractive index $n_{\text{above}}$ above the tissue equal to 1.462 representing a laterally uniform fiber-tissue boundary (no surface layout). Furthermore, the fiber dataset includes all the phase functions, a source with a fiber core and cladding diameter of 200 and $220\mu \mathrm{m}$, the numerical aperture of 0.22 and refractive index of the core equal to 1.462, a detector of total specular reflectance, and top and bottom radial detectors with numerical apertures of 0.22 and 1.0, respectively, spanning from 0 to 5 mm with 500 bins. It should be noted that the refractive index $n_{\text{above}}$ is for the fiber subset deliberately set to 1.462 (fused silica). This way, the spatially resolved reflectance given by the radial detector can then be integrated to obtain an estimate of reflectance for a detector fiber at a given source–detector separation. In other cases, such as six-around-one, six-linear-array, or single-fiber, the proper probe-tissue boundary is achieved by the surface layout selection and thus the refractive index $n_{\text{above}}$ is set to 1.0 as it represents the refractive index of the air surrounding the metal housing of the probe. Photon packets incident within the area of the probe tip reflects from the metal housing with the reflectivity set to 60%. Finally, for the six-linear-array or single-fiber, the parameter $N_{\text{fib}}$ represents the number of fibers in the layout. For example, $N_{\text{fib}}=6$ corresponds to one fiber being used simultaneously as a source and detector and the remaining five fibers as detectors in a tightly packed manner (i.e., the smallest source–detector separation is equal to the diameter of the fiber cladding). As expected, the value of $N_{\text{fib}}$ for the single-fiber setup is set to 1.

The explanation for the parameters of various phase functions used in the dataset can be found in one of our recent publications.³⁷ The Mie phase function was derived for polystyrene and fused silica microspheres ranging in diameter from 0.25 to $4\mu \mathrm{m}$ as given in Fig. 3. The refractive index of polystyrene and fused silica was obtained from Nikolov and Ivanov⁴⁴ and Malitson⁴⁵ for the wavelength of light equal to 500 nm. The surrounding medium was set to water with a refractive index of 1.337 obtained from Daimon and Masumura.⁴⁶ All the refractive indices can be retrieved from the utility materials in the package PyXOpto.

The MCML dataset comprises 11 subsets corresponding to different sources, each simulated for 30,870 combinations of optical properties and phase functions, resulting in a total number of 339,570 configurations. Each configuration was simulated with ${10}^8$ photon packets except for the linear-array subset configuration, where ${10}^9$ photon packets were used to achieve an acceptable signal-to-noise ratio (SNR) at the most distant optical fiber.

2.2.3.

MCVOX dataset

The MCVOX dataset is based on a two-layered skin model with an embedded blood vessel. Light propagation modeling for this dataset is performed using the xopto.mcvox subpackage from PyXOpto. The optical properties of the epidermis, dermis, and blood vessel were obtained from Jacques⁴⁷ and are given in Fig. 4 along with dimensional parameters of the skin model. In this dataset, a line (pencil beam) was used as a source. The energy deposition was accumulated using the fluence module and the reflectance and transmittance using a Cartesian type detector. The details of the voxelization of the energy deposition are given in Fig. 4. The MCVOX dataset comprises 26 geometrical configurations in which the depth of the embedded vessel is varied uniformly from 0.2 to 0.8 mm as is illustrated in Fig. 4. Each configuration was simulated using ${10}^9$ photon packets.

Fig. 4

Configuration of the MCVOX dataset. Energy deposition is normalized by the product of the voxel volume and the total weight of the launched packets (${\mathrm{m}}^{-3}$).

2.2.4.

SV dataset

The SV dataset is based on a single geometrical configuration given in Fig. 5 and comprises an SV for an optical fiber probe with the source and detector fibers separated by $500\mu \mathrm{m}$. The two fibers are embedded in a metal housing with reflectivity of 60% that is implemented through the surface layout module. The dimensions of the probe and optical properties of the tissue are further given in Fig. 5. The photon packet traces were recorded by storing up to 1000 events per photon packet. It should be noted that the capacity of the trace can be changed and was in this case determined beforehand since the photon packets that reach the detector rarely undergo $>1000$ events. The traces of the photon packets were filtered using the specifications of the detector fiber, which were the same as for the source fiber, i.e., numerical aperture equal to 0.22, core refractive index equal to 1.462, and core and cladding diameters equal to 200 and $220\mu \mathrm{m}$. The SV was then prepared using the embedded utility of the PyXOpto package. The spans and number of bins for the voxelization of the SV are given in Fig. 5.

Fig. 5

Configuration of the SV dataset.

2.2.5.

SFDI dataset

The SFDI dataset comprises two subsets of simulated spatial frequency domain reflectance spanning spatial frequencies from 0 to $0.8{\mathrm{mm}}^{-1}$ using 81 uniformly distributed points. The two subsets are color-coded in Fig. 6 and differ in terms of the detectors. Detectors in the first subset include total specular reflectance and top and bottom spatially resolved reflectance acquired perpendicularly to the tissue with a 10 deg acceptance cone using 4000 logarithmically spaced radially symmetric bins spanning from 0 to 150 mm. We have recently shown that such a configuration is necessary for the accurate discrete transformation of the spatially resolved reflectance into the spatial frequency domain.⁴⁸ While the detectors in the second subset also include the total specular reflectance, the spatially resolved reflectance is acquired through a detector tilted by 20 deg with an acceptance cone equal to 10 deg using 8000 logarithmically spaced $x$-symmetric bins (infinite in the $y$-direction) spanning from $-150$ to 150 mm. The latter configuration more accurately represents the realistic case where the tissue is illuminated perpendicularly, and the backscattered light is acquired with a tilted camera. Finally, it should be noted that for the first subset the spatial frequency domain reflectance can be obtained using the Hankel transformation since the spatially resolved reflectance is radially symmetric. For the second subset, one-dimensional Fourier transformation must be performed along the $x$-direction, since the spatially resolved reflectance is symmetric across the $x$-axis. The spatial frequency domain reflectance is unitless. The SFDI dataset comprises 2205 different configurations of optical properties for each subset resulting in a total number of 4410 configurations. Each configuration was simulated with ${10}^9$ photon packets to yield an acceptable SNR at long distances from the source.

Fig. 6

Configuration of the SFDI dataset.

2.2.6.

Organization of the MCDataset

The MCDataset is publicly available on GitHub ( https://github.com/xopto/mcdataset) and can be generated, customized, and accessed using the dataset utility within the PyXOpto package.

3.1.1.

Total reflectance and transmittance

Figures 7 and 8 show relative differences between the total reflectance and transmittance for the subset two-layer-100-$\mu m$-1-mm with various combinations of ${\mu }_a$, ${\mu }_s^{\prime }$, and $g$ varied in the top layer. In principle, we can observe excellent agreement between the estimations is given by the xopto.mcml and the MCML for the total reflectance and transmittance. The absolute relative errors never exceed 0.32%. However, we do observe a slight positive and negative bias in the total reflectance and transmittance, which we attribute to the choice of the random number generator (RNG) in the MCML code, which is further discussed in Sec. 3.1.2. Figures 17 and 18 in the Appendix show that replacing the RNG in the MCML code with a cryptographic-grade RNG reduces the relative differences to negligible values. Note that subsets one-layer-1-mm and one-layer-100-mm exhibited similar errors.

Fig. 7

Relative errors between the xopto.mcml and MCML simulations of the total reflectance for the dataset MCML comparison and subset two-layer-100-$\mu \mathrm{m}$-1-mm.

Fig. 8

Relative errors between the xopto.mcml and MCML simulations of the total transmittance for the dataset MCML comparison and subset two-layer-100-$\mu m$-1-mm.

3.1.2.

Spatially resolved reflectance and transmittance

Figures 9 and 10 show spatially resolved reflectance and transmittance for the subset one-layer-1-mm for various combinations of optical properties ${\mu }_a$, ${\mu }_s^{\prime }$, and $g$. The spatially resolved reflectance and transmittance completely overlap with the calculated relative errors indicating no bias. The relative errors are higher at longer distances due to the simulation noise.

Fig. 9

Spatially resolved reflectance as simulated by the xopto.mcml (black, blue, and green solid lines) and MCML (red dashed lines) for various combinations of optical properties ${\mu }_a$, ${\mu }_s^{\prime }$, and $g$ in the dataset MCML comparison and subset one-layer-1-mm. Relative errors are denoted with the same colors. Scale is provided at the right-hand side of the graphs.

Fig. 10

Spatially resolved transmittance as simulated by the xopto.mcml (black, blue, and green solid lines) and MCML (red dashed lines) for various combinations of optical properties ${\mu }_a$, ${\mu }_s^{\prime }$, and $g$ in the dataset MCML comparison and subset one-layer-1-mm. Relative errors are denoted with the same colors. Scale is provided at the right-hand side of the graphs.

3.1.3.

Energy deposition and effect of the random number generator

The energy deposition maps as simulated by the xopto.mcml and MCML are given in Fig. 11. While the maps seem similar and show no apparent differences, Fig. 12(a) clearly shows a pronounced positive and negative bias. This bias is even more pronounced in Fig. 13 (black solid line), which shows the relative errors of the energy deposition maps projected onto the $z$ or $r$ axis.

Fig. 11

Deposited energy map simulated by (a) xopto.mcml and (b) MCML for ${\mu }_a=2.5{\mathrm{cm}}^{-1}$, ${\mu }_s^{\prime }=30{\mathrm{cm}}^{-1}$, and $g=0.9$.

Fig. 12

The relative errors between the energy deposition maps simulated by (a) xopto.mcml and MCML and (b) xopto.mcml and MCML with a modified RNG based on the cryptographic-grade RNG of the OpenSSL library. Deposited energy was simulated for a medium with ${\mu }_a=2.5{\mathrm{cm}}^{-1}$, ${\mu }_s^{\prime }=30{\mathrm{cm}}^{-1}$, and $g=0.9$.

Further investigation of the MCML code showed that the bias is apparently due to the RNG ran3 from the Numerical Recipes in C.⁴⁹ To confirm our hypothesis, we have modified the MCML code to accept large sets of random bytes that can be generated using for example cryptographic-grade RNG of the OpenSSL library.⁵⁰ Using the random numbers generated by the latter RNG and repeating the simulations of the dataset MCML comparison removed all the apparent bias present in the total reflectance and transmittance (Figs. 7 and 8) and energy deposition in Fig. 11. This is confirmed by Figs. 12(b) and 13 (blue line) which show that the remaining errors are mostly governed by the simulation noise.

Fig. 13

Relative errors between the energy deposition projected onto the (a) $x$ and (b) $r$ axis. The black solid line represents relative differences between xopto.mcml and MCML, while the blue solid line represents relative differences between the xopto.mcml and MCML with a modified RNG based on the cryptographic-grade RNG of the OpenSLL library.