1.

Introduction

Polarized light has been found to be highly sensitive to medium structures and hence has been widely adopted in optical imaging to probe microstructural features inside biological tissues.^1–4 For example, the polarization status of the backscattered light can be measured to characterize the superficial layer of skin for cancer diagnostic purposes.^5,6 The measurements of tissue birefringence permit quantification of abnormalities of the retinal nerve fiber layer⁷ and cornea,⁸ as well as three-dimensional (3D) reconstruction of nerve fiber orientations inside human brains⁹ and orientations of collagen within the uterine cervix.¹⁰ By measuring the unequal absorption of left-handed and right-handed circularly polarized light, circular dichroism can rapidly determine the folding properties of proteins.¹¹ Quantification of collagen and birefringent media alignment can improve evaluation of a therapeutic strategy and its outcome in scar management.¹² Polarized light imaging (PLI) uses linearly co-polarized images subtracted by those of cross-polarized light to create a differential image based on the small population of superficially scattered photons that still retain much of the incident polarization state.⁵ PLI subtracts the large randomized population of multiple-scattered photons that produce a blinding background of diffuse light. The resulting difference image enhances the contrast of superficial tissue layers and rejects deeper tissue structures, enabling wide-field screening of epidermal or epithelial layers. Mueller matrix polarimetry and the use of various decomposition methods can also be used to pinpoint different regions and structures within biological tissue.^2,3 Accurately simulating polarized light transport inside complex tissues allow quantitative investigations of the depth response of polarized light and the perturbations produced by local tissue abnormalities.

The propagation of polarized light inside scattering media can be described by the vector radiative transfer equation (VRTE).¹³ Analytically solving the VRTE is not viable in complex media such as human tissues. Owing to its high flexibility and simplicity in programming, the Monte Carlo (MC) method, among other numerical techniques,^14–19 has been one of the limited approaches available to quantitatively model interactions of polarized light with scattering media. Depending on the vectorial representations of polarization states, polarized light MC algorithms can be largely categorized into two formalisms—Jones calculus and Mueller calculus.^2,4,20 The Jones calculus used in the electric-field MC (EMC) algorithm traces the amplitudes and phases of two orthogonal electric-field components (Jones vector) and is therefore well suited for simulating light coherence effects.²¹ On the other hand, the Mueller calculus describes the state of polarization using the Stokes vector.^22–26 The Stokes vector does not contain the absolute phase of the electrical field but allows to model unpolarized, partially polarized and fully polarized light. It can be obtained by measuring four intensity values. Similarly, Mueller matrices can be obtained through 16 intensity measurements and using Mueller matrix decomposition, quantities, such as tissue retardation, depolarization and de-attenuation can be obtained.⁴

A well-known limitation of MC methods is the long computation time. Due to the rapid emergence of massively parallel computing techniques, benefit largely from the fast advances in many core processors such as graphics processing units (GPUs), MC simulations of polarized light have seen significant speed improvement over the past decade. Several groups have reported parallel EMC implementations.^27–29 Wang et al.²⁷ presented a compute unified device architecture (CUDA)³⁰-based EMC to model coherent light in a single homogeneous slab and achieved over $370\times$ speedup compared to the central processing unit (CPU) counterpart.²¹ Ding et al.²⁹ extended the GPU-based EMC algorithm to consider multi-layered media at the expense of reduced speedup ($\sim 45\times$). In addition, Li et al.³¹ presented a CUDA-based polarized light MC algorithm to model interstitial media embedded with spherical and cylindrical scatterers. It employed a single-kernel scheme and was hundreds of times faster than its CPU version.²⁶ In 2019, Oulhaj et al.³² reported a GPU-accelerated MC algorithm to efficiently compute the sensitivity profile for polarized light inside homogeneous media. The reported GPU implementation was verified against a widely used CPU-based code developed by Ramella-Roman et al.²⁴ and reported over $150\times$ speedup.

Although these studies have demonstrated significantly improved simulation speed, most of these simulators only support layered domains and can not address the needs in modeling increasingly complex media.⁴ In the simulation of biological tissues with irregular-shaped structures, employing simplifications in domain geometries could introduce significant errors. For example, in the MC modeling of human brains, noticeable differences have been observed between layered-slab models and more anatomically realistic models such as voxel-based and mesh-based brain models.³³

In this work, we present an open-source and GPU-accelerated MC simulator to model polarized light inside 3D heterogeneous media. This MC algorithm utilizes a 3D voxelated grid to represent spatially varying distributions of spherical scatterers, characterized by their radii and densities. We use Muller calculus to update the Stokes vectors of simulated photon packets, from which we can compute various polarimetry related measurements along the surface of the domain. In the remainder of this paper, we first briefly review the steps of the meridian-plane MC algorithm.²⁴ Then we detail our GPU-implementation of this algorithm, as part of our enhanced open-source MC software—Monte Carlo eXtreme (MCX),³⁴ including the preprocessing steps to encode the distribution of particles into a 3D array data structure and optimization strategies to better use GPU resources. In Sec. 3, we validate the proposed GPU-based polarization-enabled MCX (pMCX) against the widely used CPU MC simulator “meridianMC” written by Ramella-Roman et al.³⁵ and quantify the speed improvement using several benchmarks of homogeneous and heterogeneous domains. Finally, we summarize the key findings and discuss future directions.

2.

Methods

2.2.

Implementing Meridian Plane MC in MCX

The original CPU-based meridian-plane MC program,³⁵ referred to as “mcMeridian” (stok1.c) hereinafter, is dedicated to modeling homogeneous infinite slab geometries. In contrast, the CUDA-based MCX is capable of modeling arbitrarily heterogeneous media represented by a 3D voxelated domain.³⁴ In a non-polarization MCX simulation, the domain is represented by a 3D integer array with each number representing the index or label of the tissue at each voxel. The actual optical properties of the tissue label are stored in a look-up table, with four element per tissue type: absorption coefficient ${\mu }_a$ (1/mm), scattering coefficient ${\mu }_s$ (1/mm), anisotropy $g$, and refractive index $n$. To simulate polarized light, the scattering properties of each type of scatterer must be included in addition. To simply the computation, here we only consider spherical scatterers. The radius $r$ ($\hspace{0.17em}\mu \mathrm{m}$), refractive index $n_{\mathrm{sph}}$, and volumetric number density $\rho$ ($1/\hspace{0.17em}\mu {\mathrm{m}}^3$) of the spherical particle scatters can be specified for each tissue type. When spherical particle properties are specified, the corresponding ${\mu }_s$, $g$, and elements of the scattering matrix $M(\theta )$ are pre-computed using the Mie theory³⁹ on the host (i.e., CPU). As shown in Eq. (2), the scattering matrix of homogeneous spherical scatterers consists of four independent floating-point numbers $s_{11}(\theta ),s_{12}(\theta ),s_{33}(\theta )$, and $s_{34}(\theta )$. In our implementation, the scattering parameters are sampled at 1000 evenly spaced points between 0 and $\pi$, as done in mcMeridian.³⁵ The pre-computed optical properties and scattering matrix data are then transferred to the device (i.e., GPU), as shown in Fig. 1.

Fig. 1

Media representation and media data preprocessing in a polarization-enabled Monte Carlo simulation.

3.

Results and Discussions

In this section, we first validate the aforementioned pMCX in a homogeneous slab using the single-threaded CPU-based implementation (mcMeridian²⁴) as a reference, which has been extensively used by the community and validated by experimental studies.²⁵ The speed improvement is also quantified. It is noteworthy that mcMeridian³⁵ simulates an infinite slab media geometry in the $x/y$ directions, whereas in an MCX simulation, a photon is confined inside a bounding box with user-specified dimensions.³⁴ To ensure that our speed comparison is valid, we modified the source code of mcMeridian and added an implicit bounding box to match the dimension settings in pMCX. In the first benchmark, the simulation domain is a $20\times 20\times 10{\mathrm{mm}}^3$ homogeneous slab, the Mie scattering parameters of the embedded spherical scatterer are ${\mu }_a=0{\mathrm{mm}}^{-1}$, $r=1.015\mu \mathrm{m}$, $\rho =1.152\times {10}^{-4}\mu {\mathrm{m}}^{-3}$, and $n_{\mathrm{sph}}=1.59$. The refractive index of the background medium is $n=1.33$. A monochromatic pencil beam source is positioned at the bottom center $(\mathrm{10,10},0)\mathrm{mm}$ of the domain, pointing toward the $+z$ axis and emitting horizontally polarized light at wavelength $\lambda =632.8$ nm. The initial Stokes vector of the incident beam is $\overrightarrow{\mathbf{S}}=[\mathrm{1,1},\mathrm{0,0}]$. The backscattered photons are collected by a square-shaped area detector ($20\times 20{\mathrm{mm}}^2$) placed on the boundary at $z=0\mathrm{mm}$. In this benchmark, ${10}^8$ photon packets are simulated on a desktop running Ubuntu 18.04 with an Intel i7-6700K CPU and an NVIDIA RTX 2080 GPU.

In Fig. 2, we compare the distribution of backscattered [$I,Q,U,V$] components using contour plots in MATLAB (MathWorks, Inc., Natick, Massachusetts, United States) and observed excellent agreement between mcMeridian and pMCX solutions. For further quantitative analysis, the root-mean-square errors of $I,Q,U,V$ (in ${\log }_{10}$ scale) between mcMeridian and pMCX are computed on the matching detector area ($20\times 20\hspace{0.17em}{\mathrm{mm}}^2$), reporting 0.0076, 0.0881, 0.1015, and 0.0938, respectively. We measure the total runtimes, including input data preprocessing, photon transport simulation, and output image generation, with mcMeridian reporting 18111.44 s using the Intel CPU and pMCX reporting 19.45 s on the NVIDIA RTX 2080 GPU, suggesting a $931\times$ speedup. In addition, we also benchmark simulation speeds when storing the scattering matrix data over different GPU memory locations, including global, shared, and constant memories.³⁰ The global memory implementation reports the fastest speed at 8401 photons/ms, followed by the shared memory (4965 photons/ms) and constant memory (2182 photons/ms) implementations. Although the shared memory is known to be the fastest among the three memory types, its has a very small size, up to 48 KB per block.³⁰ For storing the scattering matrix of a single species of scatterer at 1000 angular steps, a total of 16 KB memory is needed. Allocating a large amount of shared memory can lead to drastically reduced active block number, which explains the lower speed compared to the global memory case. On the other hand, constant memory also has a small size (64 KB).³⁰ It is most efficient when a memory value is being reused many times after a single read. However, the use of the rejection method requires random access to the buffer which fails to be accelerated by the constant memory due to high “cache-miss.”³⁰

Fig. 2

Contour plots of the absolute backscattered (a) $I$, (b) $Q$, (c) $U$ and (d) $V$ (in ${\log }_{10}$ scale) generated by mcMeridian (black solid lines) and pMCX (white dashed lines).

In the next benchmark, we further validate our pMCX simulator by comparing with an extended mcMeridian (with added support of layered media) in a two-layer domain. The slab-shaped simulation domain has a size of $100\times 100\times 50{\mathrm{mm}}^3$ with the thickness of the superficial layer, $d_e$, ranging between 0 and 10 mm. The Mie scattering parameters are $\lambda =632.8\mathrm{nm}$, ${\mu }_a=0.001{\mathrm{mm}}^{-1}$, sphere radius $r=0.05\mu \mathrm{m}$, number density $\rho =19.11\mu {\mathrm{m}}^{-3}$, $n_{\mathrm{sph}}=1.59$, and $n=1.33$ for the superficial layer. The bottom layer has $r=0.3\mu \mathrm{m}$, $\rho =2.198\times {10}^{-2}\mu {\mathrm{m}}^{-3}$, and the other parameters are the same as the superficial layer. These choices of $r$ and $\rho$ yield a reduced scattering coefficient ${\mu }_s^{\prime }={\mu }_s·(1-g)=1{\mathrm{mm}}^{-1}$ for both layers, along with the same absorption ${\mu }_a$ and hence approximately the same reflected intensity $I$ for all $d_e$. A pencil beam is located on the surface of the slab at $(\mathrm{50,50},0)\mathrm{mm}$, pointing toward the $+z$ axis and emitting horizontally polarized light, with the initial Stokes vector $\overrightarrow{\mathbf{S}}=[\mathrm{1,1},\mathrm{0,0}]$. A total of ${10}^7$ photon packets are simulated for both mcMeridian and pMCX.

In Fig. 3, we plot the total reflected $I$ and $Q$ components as a function of the superficial layer thickness $d_e$. We do not include total $U$ and $V$ plots because they are nearly zeros across all $d_e$ values, this is expected as $U$ and $V$ distributions sum to zeros due to symmetric positive and negative components. The outputs from mcMeridian (black solid lines) and those from pMCX (red circles) once again show excellent agreements. The plot of $Q$ increases with the thickness $d_e$ of the superficial layer, which contains smaller spherical scatters and hence stronger back-scattering than the deeper layer. With thickness $d_e$ growing from 0 to 1 cm, the total $I$ value shows a minuscule increase by 0.17% (from 0.9080 to 0.9095), as shown in the inset in Fig. 3(a), as a result of sub-diffusive scattering. The two-phase transition of $Q$ matches our expectations: when the superficial layer is very thin, the reflectance values are close to the value as if the domain is entirely filled with the bottom medium (green dashed line); as we increase $d_e$, the reflectance values asymptotically approach those determined by the media in the superficial layer (blue dashed line).

Fig. 3

Validation of pMCX in a two-layer domain. We plot the backscattered (a) $I$ and (b) $Q$ components computed by pMCX and mcMeridian as the superficial layer thickness ($d_e$) increases from 0 to 10 mm. Two dashed lines in (b) indicate back-scattered $Q$ values computed from a homogeneous slab filled only with the medium of the bottom layer (green) and that of the superficial layer (blue). The inset in (a) shows a zoom-in view of the $y$ axis of the $I$ component obtained by pMCX to demonstrate subtle variations due to sub-diffusive scattering effect.

Finally, we show simulation of a slab-shaped medium embedded with a spherical inclusion, showcasing pMCX’s capability of modeling heterogeneous domains. In this benchmark (Fig. 4), the simulation domain is a $10\times 10\times 1.2{\mathrm{mm}}^3$ slab with a spherical inclusion of radius 0.5 mm centered at $(\mathrm{6,6},0.6)\mathrm{mm}$. The inclusion and the slab share identical absorption coefficient ${\mu }_a=0.005{\mathrm{mm}}^{-1}$, reduced scattering coefficient ${\mu }_s^{\prime }=1{\mathrm{mm}}^{-1}$, and refractive index $n=1.33$. However, the Mie scatters inside both domains are different. The background medium is filled with scatterers of radius $r=0.05\mu \mathrm{m}$ and volume density $\rho =19.11\mu {\mathrm{m}}^{-3}$; the inclusion is filled with scatterers of radius $r=1\mu \mathrm{m}$ and volume density $\rho =1.11\times {10}^{-3}\mu {\mathrm{m}}^{-3}$. The choices of $r$ and $\rho$ values in either domain was computed based on the Mie theory to ensure their reduced scattering coefficients are the same. All spherical scatterings have a refractive index of, $n_{\mathrm{sph}}=1.59$. A $10\times 10{\mathrm{mm}}^2$ uniform planar light source is placed at the bottom ($z=0\mathrm{mm}$) surface, pointing toward the $+z$ axis and emitting horizontally polarized light with the wavelength $\lambda =632.8$ nm. A cyclic boundary condition is applied to the four bounding box facets at $\pm x/\pm y$ directions to approximate an infinite slab and infinite-plane source. The incident Stokes vector is $\overrightarrow{\mathbf{S}}=[\mathrm{1,1},\mathrm{0,0}]$. A total of $2\times {10}^8$ photons are simulated on an NVIDIA RTX 2080 GPU. We compare the distributions of backscattered $[I,Q,U,V]$ at $z=0\hspace{0.17em}\mathrm{mm}$, as shown in Fig. 4. Because the inclusion and background slab share the same ${\mu }_a$, ${\mu }_s^{\text{'}}$, and $n$, a regular diffuse optics forward model without polarization capability would generate no contrast to the inclusion. However, our pMCX simulation has revealed distinct image contrasts in $I$, $Q$, and $V$ images at the correct inclusion locations, suggesting the potential to detect tissue microstructure differences using polarized light. It is noteworthy that the noise level in each of the images is partially related to the anisotropy $g$ determined based on the background scatterer parameters—a larger spherical radius results in a higher $g$ value and less back-scattered photons. The significantly higher amplitude of inclusion contrast in $Q$ image compared to that in $I$ further demonstrates the advantages of using polarized imaging in detecting scattering differences compared to traditional diffuse optics where only $I$ is typically measured.

Fig. 4

Distributions of (a) $I$, (b) $Q$, (c) $U$ and (d) $V$ backscattered from a $10\times 10\times 1.2{\mathrm{mm}}^3$ slab. A spherical inclusion of radius 0.5 mm is centered at (6,6,0.6) mm.

4.

Conclusion

In summary, we report a massively parallel implementation of polarized MC algorithm in our MCX simulator for modeling the propagation of polarized light inside complex media filled with spherical scatterers. Enabled by its built-in voxel-based geometric representation, the pMCX can handle arbitrarily heterogeneous media. We have described the preprocessing steps to encode the scattering properties of spherical particles with various radii and volume densities into a 3D voxel-based data structure. We provide validation and speed benchmarks ranging from simple homogeneous to complex heterogeneous domains. In all benchmarks, our pMCX solver reports excellent match with the widely used reference solver mcMeridian, while providing a speedup nearly three orders of magnitude. In addition, we observed different GPU memory utilization efficiency among global, constant, and shared memories, with the global memory implementation yielding the highest speed and least restriction. It is noteworthy that a limitation in both mcMeridian and pMCX simulations is that all media boundaries are assumed to have matched refractive indices. We plan to further extend this work to update the Stokes vector across mismatched boundaries.