Open Access
Add to BibSonomy
Abstract
This article, Part II of an article series on GPU-accelerated Monte Carlo simulation of photon transport through turbid media, focuses on the validation of the online software Multi-Scattering. While Part I detailed the implementation of the computational model, simulated and experimental results are now compared for the distribution of the scattered light intensity. The scattering phantoms prepared here are aqueous dispersions of polystyrene microspheres of diameter D = 0.5, 2 and 5 μm and at various concentrations, resulting in optical depth ranging from OD = 1 to 17.5. The Lorenz-Mie scattering phase functions used in the simulations have been verified experimentally at low particle concentrations by analyzing the angular light intensity distribution at the Fourier plane of a collecting lens. The validation approach herein accounts for the specific light collection and image formation by the camera. The front and side surfaces of the medium are imaged and the corresponding light intensity distributions are compared qualitatively and quantitatively. It is concluded that the model enables reliable simulations over the tested parameters, offering predictive simulations of transmitted intensities with a mean relative error ≤~19% over the full range. The online software is available at: https://multi-scattering.com/
© 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
Examples of turbid media are countless and includes human tissues (e.g. blood, skin layers, muscles, brain samples, etc.), turbid liquids (e.g. waste water, wine, milk, unfiltered beer, etc.), as well as various natural environments (e.g. clouds, fog, smoke, undersea, etc.). Despite being widely different in structure and dimension, those media all share a common trait: they are challenging to image through and/or analyze, due to the occurrence of multiple light scattering. In other words, photons are interacting several times with randomly distributed scattering/absorbing centers, resulting to blur effects that reduce or fully suppresses visibility.
Turbid media are characterized by their own optical properties. Those includes the scattering coefficient $\mu _s$ (mm$^{-1}$), the absorption coefficient $\mu _a$ (mm$^{-1}$) and the optical depth $OD=(\mu _s+\mu _a)\cdot l$ ; also commonly denoted as $\tau$ in the literature, where $l$ is the length traveled by the light through the medium. The optical depth quantifies the turbidity of a medium and corresponds to the average number of scattering/absorbing events occurring along the distance $l$. As explained in Part I of this article series [1] and also described in [2,3], based on the optical depth three scattering regimes can be identified. The single scattering regime corresponds to $OD \leqslant 1$ where single scattering events are dominant. The intermediate scattering regime applies when $2 \leqslant OD \leqslant 9$, where visibility is reduced and where the dominant number of scattering event, also called scattering order, is usually close to the $OD$ value. The multiple scattering regime for $OD \geqslant 10$, where visibility is lost with standard imaging approaches.
Attempts in retrieving visibility in the range $10 \leqslant OD \leqslant 15$ have been done by means of advanced filtering strategies such as time gating [4,5] or Structured Illumination Fourier filtering [6]. However, at $OD \geqslant 15$ this task is becoming extremely challenging. Under the assumption that the number of scattering events experienced by photon packets are equally distributed and that no absorption occurs, the diffusion approximation applies within the multiple scattering regime. This is the case, for example, when a near-infrared light beam is crossing a few millimeters of skin tissues. In 1983, Ishimaru et al. [7] investigated the transmission of an optical beam through randomly distributed particles and compared their experimental data with the diffusion approximation. While a good agreement was obtained for $D=0.1\;\mu$m, disagreements were found for the $D=2\;\mu$m particles (especially at optical depth below $OD=20$), demonstrating the need for a more predictive modeling approach, such as those based on the Monte Carlo method.
The Monte Carlo method for simulation of photon transport is used in a wide range of applications, such as atmospheric and marine optics [8,9], fruit inspection [10] and spray diagnostics [11,12]. However, it is particularly relevant within the field of Near Infrared Spectroscopy (NIRS), in medicine, where the diffuse light is recorded in order to deduce the optical properties of human tissues. Those techniques include: Continuous-Wave (CW-NIRS) [13], Frequency Domain (FD-NIRS) [14], Time Domain (TD-NIRS) [15,16], Diffuse Correlation Spectroscopy (DSC) [17,18] and Gas in Scattering Media Absorption Spectroscopy (GASMAS) [19]. The common point and limitation of those techniques is that they rely on the average distance photons are propagating in the probed tissue between the source and the detector. In such case, knowledge of photon transport and absorption is very important; a quantity that can be predicted by means of a Monte Carlo simulator that is validated experimentally.
One of the first comparison of results between Monte Carlo and the diffusion equation was provided in 1989, in the Part I article from Flock et al. [20]. In Part II [21] of this study, the authors provided a comparison of experimental and Monte Carlo estimates of the relative fluence versus depth within various solutions of intralipids. The predictions of the Monte Carlo model were found to be in good agreement with those measurements. In 1993, Keijzer made a comparison with the numerical results from her Monte Carlo Model [22] using the experimental results from Ishimaru et al. [7] on light transmission through latex particles. It was also shown that her Monte Carlo calculations were providing much closer agreements than the diffusion approximation, especially for the case of highly forward scattering phase functions. From this validation work, Keijzer further used her Monte Carlo model for a variety of studies, including calculating the light distributions in arterial tissue [23] as well as optimizing the beam diameter in port wine stain treatment [24].
In those studies light transport in medical tissues was investigated. Within this context the authors relied on the use of Henyey-Greenstein phase functions [25], which remains extensively used today [26,27]. This approach is reasonable for intralipids solutions where the scatterers are of size smaller than the wavelength. However, in her thesis [22], Keijzer mentioned that using Lorenz-Mie scattering phase functions [28,29] would have lead to more predictive numerical results. For the case of spherical scatterers of size ranging from $D=0.1\lambda$ to 1000$\lambda$ and assuming independent scattering, the correct scattering phase function used in a Monte Carlo simulator can be derived by means of the Lorenz-Mie theory.
One of the earliest validation work of a Monte Carlo model using Lorenz-Mie scattering phase functions has been reported in 1991 by Donelli et al. [30]. In this article the authors have compared the measured and calculated irradiance on the image plane of a point source after crossing various suspensions of polystyrene microspheres in water. Results for different values of sphere diameter ($D=0.3$, 1 and 15.7 $\mu$m) and particle concentration were investigated, where the transmitted light intensity was measured using a single detector. By varying the active detection area, the point spread function and the modulation transfer function could be deduced at optical depth up to $OD=6$. Another validation work has been performed in 2007 by Berrocal et al. [3] where the spatial light intensity distribution was considered for both transmitted and side scattering detection at optical depth up to $OD=10$ and for particles of size $D=1$ and 5 $\mu$m. While an acceptance angle was, in this study, used to consider the portion of light being collected by the camera objective, the modeling of image formation was not included in the simulation.
Thanks to its detailed description and the accessibility to its source code, MCML (Monte Carlo Multi-Layered code developed in 1994 [31]) has over the years become the most popular and widely used model for the simulation of photon transport through skin tissues. MCML was upgraded in 2014 into the multi-voxels method [32]. The code, named mcxyz, is written in C and supports programs written in MATLAB. Nowadays, many modern GPU-based Monte Carlo models [33–35], are being validated by comparing their simulated results with those from the MCML or mcxyz. However, it should be noted that while those comparisons provide some level of confidence for accurately predicting photon transport and energy deposition, a complete validation requires the use of experimental data covering a variety of scattering media with different optical depths and scattering phase functions.
This article aims at validating the open-access software Multi-Scattering by comparing experimental and simulated images of the light intensity distribution for the case of a laser beam illuminating a 30 mm side cubic scattering volume of well-controlled optical properties. These phantoms are aqueous dispersions containing polystyrene microspheres of diameter, D = 0.50 $\mu$m, 2.07 $\mu$m or 4.94 $\mu$m with standard deviations of less than 10%. Thus, each liquid mixture is considered monodispersed, and hereafter referred to as D = 0.5 $\mu$m, 2 $\mu$m and 5 $\mu$m respectively. In the first part of this validation work, the Lorenz-Mie scattering phase function of each microsphere size are experimentally verified. This is achieved using diluted dispersions respecting the single scattering regime. Images of the Fourier plane of the collecting lens are generated both experimentally and via simulation and compared.
The second part of this article concerns the validation of the model in the intermediate- and multiple scattering regimes, where images of the light intensity distribution are recorded both for the forward and side scattering detection. In this case, a total of 21 different phantoms have been used, involving the three particle sizes mentioned above and seven optical depths ranging from $OD=2$ to $OD=17.5$. Note that modeling of photon collection and image formation by the camera system is implemented in this study. This approach provides a faithful experiment/simulation comparison and differs from what has been presented in the past by the authors [2,3], where ray tracing of image formation was omitted. Additionally, the range of optical depth has been extended from the intermediate scattering regime to the multiple scattering regime where $OD > 10$. Thus, the comparison results presented here serves to affirm the validity of the updated and extended GPU-accelerated Monte Carlo model.
2. Experimental verification of the scattering phase functions
2.1 Lorenz-Mie scattering phase functions
For any Monte Carlo simulation of light propagation through a turbid medium, the scattering phase function is a key parameter as it defines the angular distribution of scattered light after each scattering event. When using experimental phantom media, scattering phase functions are usually approximated by inserting the desired anisotropy factor $g$ in the Henyey-Greenstein equation [25]. However, a more accurate way of deducing the scattering phase function for model validation purposes is to use phantoms consisting of spherical scatterers of well-known size and refractive index where the Lorenz-Mie theory can be applied. This theory provides the solution to Maxwell’s equations when considering a plane wave interacting with a homogeneous sphere [28,29]. For spheres of size comparable to or larger than the wavelength the solution predicts an angular distribution which, due to diffraction and interference, exhibits lobes of higher intensity for some preferential directions. For micrometer scale particles illuminated with visible light (380-750 nm) the forward scattering lobe is strongly dominating the scattering phase function, which then is typically plotted in logarithmic scale. As shown in Fig. 1, the Lorenz-Mie scattering phase functions of the polystyrene spheres used in this study have been calculated for $D=0.5$, 2 and 5 $\mu$m, at unpolarized $\lambda$ = 473 nm illumination. The corresponding complex refractive index of the microspheres was set to $n = 1.594 - 0.00033i$ following the measurement data reported in [36] and the refractive index of the surrounding water was set to $n = 1.333 - 0.0i$.
Fig. 1. Lorenz-Mie scattering phase functions (in logarithmic scale) showing the distribution of the light intensity [a.u.] for the polystyrene spheres of diameter $D=0.5\;\mu$m, 2 $\mu$m and 5 $\mu$m illuminated at $\lambda$ = 473 nm with unpolarized light. The indicated maximum intensities are normalized to the maximum value of the 5 $\mu$m scattering phase function.
2.2 Description of the experimental and simulated setups
The light source used in the experimental setup is a continuous wave diode laser, emitting at $\lambda$ = 473 nm. It illuminates a 30 mm cubic glass cell containing a solution of monodisperse polystyrene spheres manufactured by Bangs Laboratories, Inc. The number density of spheres of each size has been calculated to fix the optical depth through the cuvette to $OD=1$, as detailed in Table 1. A plano-convex spherical lens of diameter 45 mm and 200 mm focal distance collects the light transmitted through the cuvette. At the focal distance of the lens $f$ a flat white screen displays the scattered light, while a 2 mm pinhole located at the focal point allows the ballistic i.e., non-scattered, component to pass through it. The diffuse reflection from the screen is imaged using a 16 bit 5.5 Mpixel Andor Zyla sCMOS camera. The experimental setup is illustrated in Fig. 2.
Fig. 2. Illustration of the experimental setup: The laser beam illuminates an aqueous dispersion of $OD=1$ and the scattered light is collected by a positive spherical lens. A white screen is positioned at the Fourier plane of the collecting lens where a pinhole lets the ballistic light pass through. The light intensity profile reflected from the white screen is imaged by a sCMOS camera.
Table 1. Extinction cross-sections [$\sigma _e$] and number densities [N] of the aqueous polystyrene sphere dilutions corresponding to $OD=1$ through $L=30$ mm path length.
From the recorded images the relative light intensity distribution as a function of scattering angle can be deduced. The pinhole location corresponds to the 0$^{\circ }$ scattering angle and each given radial distance $d$ from the pinhole corresponds to a specific scattering angle $\theta _s$, given by:
Due to the experimental configuration, the angle density on the screen is itself a function of $d$ as both average distance and angle of incident varies. Thus, the intensity $I_{\theta _s}$ is related to the camera pixel values through a scaling factor, such that:
where $P$ is the pixel value recorded by the camera and $A_p$ is the pixel area. Most characteristics of the experimental setup described above can accurately be simulated by the Multi-Scattering software as listed below:- • The incident beam is modeled using a matrix corresponding to the image of the laser beam crossing the cuvette containing distilled water without scattering particles (nearly non-absorbing scattering medium).
- • The dimension and geometry of the scattering medium corresponds to the 30 mm cubic cuvette.
- • The spherical lens collecting the scattered light is located 10 mm away from the cuvette and has a diameter of 45 mm.
- • The scattering phase functions are those presented in Fig. 1 assuming a light beam of 473 nm wavelength illuminating polystyrene spheres in water.
- • For each particle size the extinction cross-section $\sigma _e$ is deduced from the Lorenz-Mie theory and given in Table 1. The number density $N$ of particles are also indicated such that an optical depth $OD=1$ is obtained over the 30 mm the cuvette, resulting in an extinction coefficient of $\mu _e = 0.033$ mm$^{-1}$.
The angular distribution of the scattered light is saved directly for photons reaching the collecting lens. Thus $I_{\theta _s}$ is obtained directly which differs from the experimental data where Eq. (2) must be applied. To consider the pinhole in the screen, the contribution of the ballistic photons is removed from the simulated data. The imaging system, involving the performance of the objective lens and the response of the camera sensor have not been considered in the simulation and may result in differences between the simulated and experimental data. Details of the simulated setup are given in Fig. 3. For each simulation 10 billion photons were launched corresponding to $\sim$30 seconds duration.
Fig. 3. Description of the simulated optical configuration. A matrix of the experimental laser profile is used as input data to launch 10 billion photons into the 30 mm side cubic scattering medium. Light is collected by a 45 mm lens positioned at 10 mm away from the scattering volume. Photons collected by the lens are then binned as a function of their scattering angle in order to deduce a photon count per angle proportional to $I_{\theta _s}$. The dimensions/geometry of the setup as well as the optical properties of the scattering medium considered in the simulation are mimicking the experiment described Fig. 2.
2.3 Comparison of results
The experimental light intensity distribution at the Fourier plane of the collecting lens is shown at the top Fig. 4(a) for each size of the polystyrene spheres. Due to symmetry, only half of the scattering pattern is displayed while the other half has been generated via simulation using Multi-Scattering. The data have been normalized by dividing each pixel value (in photon counts) with the respective mean value of all pixels for each image. The images are plotted in logarithmic scale with a common color coding, ranging from the joint minimum to the joint maximum. The corresponding radial profiles of the light intensity are shown below the images, and it is seen that a good agreement between the experimental and simulated profiles is obtained, especially for the case of the 2 and 5 $\mu$m polystyrene spheres where the ring patterns have similar position, frequency, and relative intensity. The minor disparities may be related to variations of particle diameter around the nominal size. However for the 0.5 $\mu$m spheres (and to a lesser degree for the 2 $\mu$m spheres), the main disparity is the larger intensity level observed at scattering angles close to 0$^{\circ }$ for the experimental data. The exact reasons for those differences are not fully known by the authors, but may be explained as follows: Despite the attempt of removing the ballistic light using a pinhole, this signal may originate from the angular spread of the incident laser beam. In the simulation, perfectly collimated rays are instead assumed without considering any small divergence of the incident laser beam. The portion of the measured light originating from such divergence is expected to decrease as particle size increase, due to the increase of the highly forward scattering lobe.
Fig. 4. (a): Experimental and simulated images at the Fourier plane of a collecting lens, for a laser beam ($\lambda$ = 473 nm) crossing an aqueous dispersion of polystyrene spheres (0.5, 2 and 5 $\mu$m diameter and $OD=1$ optical depth). Intensity is in units of mean pixel value in each image $\overline {I_{px}}$. Below the images are the corresponding angular intensity profiles derived from the images. (b): Contribution of each scattering orders from the simulated image.
The light intensity contribution of each scattering order has been obtained from the simulated data and is given in Fig. 4(b). It is observed here that the relative distribution is quite similar for each sphere diameter. As the optical thickness is set to $OD=1$, the ballistic and single scattering photons are the most dominant and here contribute equally to $\sim$37.5% of the total detected light intensity. Most of the multiple light scattering corresponds to scattering orders ranging from 2 to 5. This contribution deforms and blurs the ring pattern from single light scattering. As ballistic photons are rejected the formed imaged consist of $\sim$60% and $\sim$40% of light intensity from single and multiple scattering, respectively.
From the comparison of these results, it is concluded that the scattering phase function generated by the Lorenz-Mie theory depict the spreading of light intensity well, especially for the 2 and 5 $\mu$m spheres.
3. Validation of the model in the intermediate- and multiple scattering regimes
3.1 Description of the experimental setup
As illustrated in in Fig. 5 the experimental setup used for the validation of the Monte Carlo model generates images of the scattered light intensity distribution from the phantom medium, for both the forward (180$^{\circ }$) and side (90$^{\circ }$) scattering directions. While side scattering is low within the single scattering regime, this arrangement allows for a more complete comparison between experiment and simulation in the higher scattering regimes. The 473 nm laser beam is coupled into an optical fiber that is fixed on a vibrating device to reduce speckle effects, producing a more homogeneous beam profile. At the exit of the fiber, a collimating lens is forming a 15 mm wide beam and the central part of this beam is selected using a 3.8 mm circular aperture. The resulting illumination beam has 135 mW power and its Full Width at Half Maximum equals to FWHM = 2.5 mm.
Fig. 5. Experimental setup used for generating images of the light intensity distribution at the exiting front and side faces (indicated by the red lines) of the phantom medium. The laser power is continuously monitored by directing 10% of the beam into a power meter while the light intensity is accurately adjusted using the desired combinations of neutral density filters (21 attenuation levels ranging between 1 and $10^{-6}$). The illuminated phantom is a 30 mm cubic glass cuvette containing an aqueous dispersion of scattering microspheres of known size and concentration. The front and side face of the cuvette are imaged by the two identical sCMOS cameras and objective lenses.
To continuously monitor laser stability a small portion of the light is directed onto a power meter using a 90/10 cube splitter. This is used to correct for laser irradiance fluctuations. To avoid saturation and optimize the dynamic range of the camera sensor the incident laser intensity is set to a desired value by means of two neutral density filter wheels. This provides a combined attenuation factor that can be adjusted between 1 and $10^{-6}$ within 21 different configurations of filter pairs. The largest attenuation factor was used to image the beam crossing purified water only (without any scatterers), serving as the incident light distribution $I_{i}$ and irradiance reference case, shown in Fig. 6. The measured scattered light intensity distributions from the phantoms $I_{raw}$ can then be scaled to be comparable to $I_{i}$ through:
where $A_p$, $P_p$, $A_i$ and $P_i$ are the attenuation factors and power meter values for the phantom measurement and the incident light measurement respectively.Fig. 6. Experimental incident light distributions $I_i$ for the three particles sizes $D=0.5$, 2 and 5 $\mu$m. The maximum intensity indicated has been used to normalize the image results given in Fig. 8, 9 and 10. The data is available in Dataset 1 [37], Dataset 2 [38], and Dataset 3 [39], respectively.
A large variety of aqueous monodisperse solutions of polystyrene spheres have been imaged in this study. For each particle size seven phantoms have been prepared accordingly, where the optical depth $OD$ ranges from 2 to 14 for both the 0.5 and 2 $\mu$m spheres. However, as the highly forward scattering lobe of the 5 $\mu$m spheres, shown in Fig. 1, reduces light spreading, the optical depth was for this case increased and ranging from 2.5 to 17.5. The scattering solution at the highest particle concentration was first prepared and then diluted with purified water to obtain each respective solutions at lower optical depth. Note that the number density $N$ of polystyrene spheres, which is given at $OD=1$ in Table 1, can be deduced for any other optical depth as: $N_{OD}=N_{OD=1} \cdot OD$.
As shown in Fig. 5, two identical imaging systems (consisting each of a 5.5 Mpixel sCMOS Andor Zyla camera mounted with a 50 mm focal Nikon objective) collect simultaneously the forward and side scattered light respectively; and image the cuvette walls. This configuration matches the virtual configuration used in the Multi-Scattering software, where the specific light collection and image formation are considered in the model. By adding extension rings the optical arrangement was close to a 4$f$ imaging system, resulting in an image magnification of 1.17 and a pixel resolution of 5.5 $\mu$m/pixel. Note that for an imaging system of magnification close to 1 the light-gathering capacity is not well represented by f-number (here, F$\#$ 1.4). Instead, the working f-number F$\#_w \approx (1+M) \cdot$F$\#$ can be used, which is equivalent to a maximum collection angle of:
where $D_l$ is the diameter of the collecting lens, $f$ is the focal distance and $M$ the magnification [40]. From the characteristics of the imaging system used here, it is deduced that $\theta _w=9.5^{\circ }$. For comparison purpose, the exposure time of both cameras are identical and set to $t=30$ ms. The full dynamic range (16-bit) of the front view camera sensor was optimized by choosing the appropriate pair of neutral density filter. However, this configuration resulted in a low signal-to-noise ratio for the side view at low values of $OD$.3.2 Description of the simulated setup
Most aspects of the experimental setup described in section 3.1 can be replicated virtually using the Multi-Scattering software as visualized in Fig. 7. Image formation is simulated here by defining the characteristics of a virtual imaging lens (a disk of given size) and a collection array (a 2D matrix of given dimension). Images are formed from x- and y-coordinates of the collected photons in the array. By matching the input parameters to those of the experimental setup, the simulated images can directly be compared to their experimental counterparts. The characteristics and optical properties of the simulations are listed below:
- • The experimental incident light distribution image $I_i$ of the beam crossing purified water is used in the simulation to accurately replicate the experimental incident photon distribution.
- • The scattering medium is a 30 mm cubic volume where the optical depth is fixed to match those from the experimental phantoms.
- • The imaging lenses have a diameter of 37 mm and a working distance of 92 mm exactly matching the characteristics of the objective lenses.
- • An acceptance angle of 9$^{\circ }$ is applied in the Fourier plane to closely mimic the working f-number and the maximum collection angle $\theta _w$ defined in Eq. (4).
- • A magnification factor of 1.17 of the photon coordinates ensures that the simulation field of view matches the experimental one.
- • Each collection array is 2160 $\times$ 2160 pixels, which equals the pixel area of interest from the camera sensors.
- • The number of photons launched per simulation ranges from 10 billion for the most dilute cases ($\sim$45 seconds running time) to 310 billion for the most concentrated cases ($\sim$3000 seconds running time). A total number of 21 simulations have been set up and run.
Fig. 7. Description of the setup used to generate the simulated images using Multi-Scattering. This is a virtual replica of the experimental setup in Fig. 5, with the indicated dimensions. The experimental light incident distribution $I_i$ is used as incident photon distributions. Image formation of the front and side of the cuvette, indicated with red lines, is simulated by means of virtual lenses and detector arrays.
All bullet points given above, aim at obtaining simulated data as close as possible from the experimental data. This includes selecting and binning the simulated photons into pixels in the same way as for the experiment where identical collection characteristics, number of pixels and magnification factor have been employed. However, three notable differences remain between the experiment and the Monte Carlo simulation used here: First, the glass cuvette surrounding the scattering medium is not considered in the model and any reflections from the glass walls are not simulated. Second, the top surface of the experimental liquid phantom is curved due to surface tension, while the simulated phantom is a perfect cube. Finally, the simulation method does not account for any wave properties of light, such as beam divergence or interference. These effects may cause disparities in comparison with the experimental images.
3.3 Comparison of results
To obtain quantitative light intensities that are comparable between experiment and simulation, the scattered light distribution $I$ is normalized by the maximum incident photon count per pixel such as: $I / (\max I_i)$. The values used as $(\max I_i)$ for each particle size have been averaged to suppress noise variations and are given in Fig. 6. This normalization is applied to both the front- and the side view data ($I_t$ and $I_s$ respectively). For the simulated results, the same normalization is performed, but by considering the number of photons initially launched per simulation.
The experimental and simulated images for the polystyrene spheres of 0.5, 2 and 5 $\mu$m diameter are presented in Fig. 8, 9 and 10 respectively. In each figure, six cases of optical depth are presented transiting from the intermediate to the multiple scattering regimes. Images of both the front and the side face of the scattering cuvette are given, where photon spreading, and light transmission can be analyzed. This extensive set of data aims at comparing results: 1) between simulated and experimental data, 2) between the side and front views, 3) between optical depths, and 4) between scatterers of different sizes having their own scattering phase functions. The colormap for all images given in Fig. 8, 9 and 10 have been scaled to their own respective maximum in order to emphasize the spatial distribution of the intensity. As expected, and well predicted by the simulation, the incident laser beam is gradually broadening when increasing the optical depth. This effect is more important for smaller particles which have a scattering phase function with a less pronounced forward scattering lobe. A detailed analysis of this broadening can be deduced from Fig. 11 where the 1D profile of the transmitted light is given for all cases. The profiles are shown in logarithmic scale and are derived by vertically integrating the relative intensity over a 3 mm wide horizontal band covering the central part of the front view images. For the 0.5 $\mu$m particles and at low $OD$, it is seen that the intensity profile is dominated by a broad central peak with a clear cut-off on its edges. This signal, that depicts the exact shape of the incident beam, reveals the contribution of the ballistic light. As optical depth increases, the profile of the laser beam becomes less and less perceivable, ultimately disappearing at $OD=12$. This transition is well predicted by the simulated profiles. For larger scatterers, such as the 5 $\mu$m spheres, the profiles are smoother with no clear cut-off on the edge of the beam. Thus, the differentiation between the ballistic and scattered components is not as distinct as for the 0.5 $\mu$m spheres. This is explained by the smooth deviation of photon propagation after each scattering event, due to the more pronounced forward scattering lobe for the larger spheres (see Fig. 4). Those experimental observations are well mirrored by the simulations. Some noticeable differences between the experimental and simulated profiles can, however, been identified. For instance, the peak intensity is larger in the experimental data for the case of the 0.5 $\mu$m particles. The authors draw parallels between these characteristics and the disparities observed during the verification of the scattering phase functions (see sub-section 2.3) where an unexpected light peak was experimentally measured near the 0$^{\circ}$ scattering angle. As this is a systematic discrepancy it seems to most likely originate from an unaccounted divergence of the incident beam or internal reflections within the cuvette. Additionally, towards the edges of the cuvette the experimental profiles flatten out as a result of the limited dynamic range of the camera detector.
Fig. 8. Experimental and simulated images of the front (180$^{\circ }$) and side (90$^{\circ }$) view of the scattering cuvette illuminated with a 473 nm laser beam. Here the cuvette is containing an aqueous dispersion of 0.5 $\mu$m polystyrene spheres at various optical depths, as indicated on the left-hand side. The data are the scattered light intensities $I_t$ and $I_s$ divided by the maximum of the incident distribution $\max I_i$. Thus, the intensity values given here are quantitative and directly comparable between each other. The underlying data $I_t$ and $I_s$, both experimental and simulated, and $I_i$ are available in Dataset 1 [37].
Fig. 9. Experimental and simulated images of the front (180$^{\circ }$) and side (90$^{\circ }$) view of the scattering cuvette illuminated with a 473 nm laser beam. Here the cuvette is containing an aqueous dispersion of 2 $\mu$m polystyrene spheres at various optical depths, as indicated on the left-hand side. The data correspond to the scattered light intensities $I_t$ and $I_s$ divided by the maximum of the incident distribution $\max I_i$. Thus, the intensity values given here are quantitative and directly comparable between each other. The underlying data $I_t$ and $I_s$, both experimental and simulated, and $I_i$ are available in Dataset 2 [38].
Fig. 10. Experimental and simulated images of the front (180$^{\circ }$) and side (90$^{\circ }$) view of the scattering cuvette illuminated with a 473 nm laser beam. Here the cuvette is containing an aqueous dispersion of 5 $\mu$m polystyrene spheres at various optical depths, as indicated on the left-hand side. The data correspond to the scattered light intensities $I_t$ and $I_s$ divided by the maximum of the incident distribution $\max I_i$. Thus, the intensity values given here are quantitative and directly comparable between each other. The underlying data $I_t$ and $I_s$, both experimental and simulated, and $I_i$ are available in Dataset 3 [39].
Fig. 11. Experimental and simulated intensity profiles for the front(180$^{\circ }$) view at various $OD$ for the three particles sizes 0.5, 2 and 5 $\mu$m. The profiles are extracted from a 3 mm band at the center of the front view images shown in Fig. 8, 9 and 10.
To analyze and quantify the broadening of the transmitted light, the relative FWHM of the profiles given in Fig. 11 have been calculated in Fig. 12. Here, the ratio between the FWHM of the profiles and of the profile of $I_i$ is plotted against the optical depth. On the overall, it is seen that the broadening is larger in the simulation than in the experiment. For $D=0.5\;\mu$m the ratio initially increases very slowly with $OD$, as the ballistic component dominates the profiles. However, once the ballistic component becomes negligible, at $OD > 6$, the broadening rapidly increases with $OD$. For $D=2\;\mu$m a similar transition is observed but at $OD > 8$. However, for $D=5\;\mu$m the transition is more gradual. A relative FWHM of 6, 5 and 2 are experimentally measured for the 0.5, 2 and 5 $\mu$m particles respectively, confirming the significant influence of the shape of the scattering phase function on light broadening.
Fig. 12. Experimental and simulated FWHM in the front (180$^{\circ }$) view, relative to the FWHM of the incident beam at various $OD$ for the three particles sizes $D=0.5$, 2 and 5 $\mu$m. The values are calculated from the front view profiles shown in Fig. 11.
Figure 13 shows the one-dimensional normalized intensity profiles of all measured cases for the side (90$^{\circ }$) view detection. The graphs are derived by vertically integrating each image and normalizing the resulting curves. Note that the top part of all images was ignored to suppress the unwanted intensity contribution from the liquid surface observed experimentally. It is seen that when increasing the $OD$ the curves are shifted towards the entrance of the cuvette. This trend is accurately predicted by Multi-Scattering as indicated by the dashed lines marking the top of each profile. Note that those shifts depend also on particle size. For instance, at $OD = 10$, the maximum intensity occurs at 12.22 mm, 10.76 mm and 12.02 mm experimentally; and 12.22 mm, 10.96 mm and 12.01 mm in the simulation, for the 0.5, 2 and 5 $\mu$m particles respectively.
Fig. 13. Experimental and simulated intensity profiles for the side (90$^{\circ }$) view detection at various $OD$ for the three particles sizes 0.5, 2 and 5 $\mu$m. The profiles are extracted from the vertical mean of the side view images shown in Fig. 8, 9 and 10.
Despite these agreements on peak position, the main disparity between the experimental and simulated data is the flattening of the experimental profiles due to the presence of an intensity "base-line". A possible explanation for this additional intensity contribution is light being reflected by the glass walls back into the cuvette. Reflections explain why the intensity contribution is strongly enhanced at the exit of the cuvette for low $OD$ where more power is transmitted through the sample. Conversely at the entrance of the cuvette, the intensity contribution is enhanced when $OD$ is increased, due to a larger amount of back-scattered light.
Figure 14 shows how the experimental and simulated transmissions vary from the Beer-Lambert law [41]. The values are derived from the image data through $T = \sum _{AOI} I_t / \sum _{AOI} I_i$, where the transmission $T$ is calculated over two different areas of interest (AOI): The blue curves $T^{all}$ consider the AOI of the full image; while the red curves $T^{beam}$ consider the AOI of only the central beam. It is seen from those results that the deviation of the transmitted light intensity from the prediction of the Beer-Lambert law is well depicted by the simulations. The transmission gradually deviates from the Beer-Lambert law with:
- • Increasing $OD$, as the more multiple scattering occurs and is collected.
- • Increasing particle size. This is the result of the more forward scattering phase functions, as seen in Fig. 1, allowing photons to be transmitted despite multiple scattering events.
- • Increasing the collection area. $T^{all}$ deviates more strongly than $T^{beam}$ due to the collection of additional scattered photons.
Fig. 14. Experimental ($T_{exp}$) and simulated ($T_{sim}$) transmission comparison and relative error ($RE$), at all $OD$ for the three particles sizes $D = 0.5$, 2 and 5 $\mu$m. The red curves show the mean transmission intensity at the location of the incident beam only while the blue curves show the total transmission. The values are extracted from the front view images given in Fig. 8, 9 and 10.
The diagrams in the lower part of Fig. 14 show the relative error $RE$ between the experimental and simulated transmission, as well as the corresponding mean relative error $MRE$ defined as:
4. Conclusions
In this study the Monte Carlo software, Multi-Scattering, has been rigorously validated by means of a variety of controlled scattering media containing polystyrene spheres of well-known size and concentration. The Lorenz-Mie scattering phase functions used in the simulation have been verified experimentally where differences have been observed for the 0.5 $\mu$m where a high intensity peak, near 0$^{\circ }$ scattering angle, was visible in the experiment. For the 2 $\mu$m and 5 $\mu$m this effect was significantly reduced, and good agreements were found confirming the suitability of the scattering phase functions used here.
To cover as many scattering situations as possible the optical depth was spanning from $OD=2$ to 14 for the 0.5 and 2 $\mu$m spheres and from $OD=2.5$ to 17.5 for the 5 $\mu$m spheres, for both font and side scattering detection. Light transmission has been quantified through comparison with the same non-scattered incident light used for both simulation and experiment in a faithful comparison. It is found that the Monte Carlo model is predicting the transmitted light intensity well, with a mean absolute relative error over the complete range of optical depth of 16% for light confined in the central beam (13%, 24% and 11% for the 0.5, 2 and 5 $\mu$m spheres respectively) and 19% for all transmitted light. Additionally, the broadening of the transmitted beam has been analyzed. It was observed that the simulated results tend to slightly overestimate light broadening. Such effects are in-line with the differences observed in the scattering phase functions. For the side scattering detection, the location of the maximum light intensity has been well predicted. However, differences originating from reflections in the glass cuvette were observed. Such reflection effects are affecting all experimental data, directly contributing to the data disparities between experiment and simulation.
To conclude, the validation work performed here has proven the high-degree of reliability of Multi-Scattering for qualitatively and quantitatively predicting photon transport through scattering media. The software can therefore confidently be used for the prediction of photon transport for a variety of research and education studies.
Funding
Stiftelsen för Strategisk Forskning (ITM17-0152); H2020 European Research Council (638546).
Disclosures
The authors declare no conflicts of interest
Data availability
Data underlying the results presented in this paper are available in Dataset 1 [37] Dataset 2 [38], and Dataset 3 [39].
References
1. J. Jönsson and E. Berrocal, “Multi-scattering software: part i: online accelerated monte carlo simulation of light transport through scattering media,” Opt. Express 28(25), 37612–37638 (2020). [CrossRef]
2. E. Berrocal, “Multiple scattering of light in optical diagnostics of dense sprays and other complex turbid media,” Ph.D. thesis, Cranfield University (2006).
3. E. Berrocal, D. L. Sedarsky, M. E. Paciaroni, I. V. Meglinski, and M. A. Linne, “Laser light scattering in turbid media Part I: Experimental and simulated results for the spatial intensity distribution,” Opt. Express 15(17), 10649–10665 (2007). [CrossRef]
4. M. Rahm, M. Paciaroni, Z. Wang, D. Sedarsky, and M. Linne, “Evaluation of optical arrangements for ballistic imaging in sprays,” Opt. Express 23(17), 22444–22462 (2015). [CrossRef]
5. D. Sedarsky, E. Berrocal, and M. Linne, “Quantitative image contrast enhancement in time-gated transillumination of scattering media,” Opt. Express 19(3), 1866–1883 (2011). [CrossRef]
6. E. Berrocal, S.-G. Pettersson, and E. Kristensson, “High-contrast imaging through scattering media using structured illumination and Fourier filtering,” Opt. Lett. 41(23), 5612–5615 (2016). [CrossRef]
7. A. Ishimaru, Y. Kuga, R. L.-T. Cheung, and K. Shimizu, “Scattering and diffusion of a beam wave in randomly distributed scatterers,” J. Opt. Soc. Am. 73(2), 131–136 (1983). [CrossRef]
8. D. Doxaran, E. Leymarie, B. Nechad, A. Dogliotti, K. Ruddick, P. Gernez, and E. Knaeps, “Improved correction methods for field measurements of particulate light backscattering in turbid waters,” Opt. Express 24(4), 3615–3637 (2016). [CrossRef]
9. H. S. Ugulen, H. Sandven, B. Hamre, A. S. Kristoffersen, and C. Sætre, “Analysis of multiple scattering errors in lisst-vsf volume scattering function measurements using monte carlo simulations and experimental data,” Opt. Express 29(8), 12413–12428 (2021). [CrossRef]
10. M. L. Askoura, F. Vaudelle, and J.-P. L’Huillier, “Multispectral measurement of scattering-angular light distribution in apple skin and flesh samples,” Appl. Opt. 55(32), 9217–9225 (2016). [CrossRef]
11. E. Berrocal, D. Y. Churmakov, V. P. Romanov, M. C. Jermy, and I. V. Meglinski, “Crossed source–detector geometry for a novel spray diagnostic: Monte carlo simulation and analytical results,” Appl. Opt. 44(13), 2519–2529 (2005). [CrossRef]
12. E. Berrocal, I. Meglinski, and M. Jermy, “New model for light propagation in highly inhomogeneous polydisperse turbid media with applications in spray diagnostics,” Opt. Express 13(23), 9181–9195 (2005). [CrossRef]
13. F. Scholkmann, S. Kleiser, A. J. Metz, R. Zimmermann, J. M. Pavia, U. Wolf, and M. Wolf, “A review on continuous wave functional near-infrared spectroscopy and imaging instrumentation and methodology,” NeuroImage 85, 6–27 (2014). Celebrating 20 Years of Functional Near Infrared Spectroscopy (fNIRS). [CrossRef]
14. S. Fantini and M. A. Franceschini, Chapter 7: Frequency-Domain Techniques for Tissue Spectroscopy and Imaging (2001).
15. D. T. Delpy, M. Cope, P. van der Zee, S. Arridge, S. Wray, and J. Wyatt, “Estimation of optical pathlength through tissue from direct time of flight measurement,” Phys. Med. Biol. 33(12), 1433–1442 (1988). [CrossRef]
16. A. Torricelli, D. Contini, A. Pifferi, M. Caffini, R. Re, L. Zucchelli, and L. Spinelli, “Time domain functional NIRS imaging for human brain mapping,” NeuroImage 85, 28–50 (2014). Celebrating 20 Years of Functional Near Infrared Spectroscopy (fNIRS). [CrossRef]
17. E. M. Buckley, A. B. Parthasarathy, P. E. Grant, A. G. Yodh, and M. A. Franceschini, “Diffuse correlation spectroscopy for measurement of cerebral blood flow: future prospects,” Neurophotonics 1(1), 011009 (2014). [CrossRef]
18. W. Liu, R. Qian, S. Xu, P. Chandra Konda, J. Jönsson, M. Harfouche, D. Borycki, C. Cooke, E. Berrocal, Q. Dai, H. Wang, and R. Horstmeyer, “Fast and sensitive diffuse correlation spectroscopy with highly parallelized single photon detection,” APL Photonics 6(2), 026106 (2021). [CrossRef]
19. S. Svanberg, “Gas in scattering media absorption spectroscopy - from basic studies to biomedical applications,” Laser Photonics Rev. 7(5), 779–796 (2013). [CrossRef]
20. S. Flock, M. Patterson, B. Wilson, and D. Wyman, “Monte carlo modeling of light propagation in highly scattering tissues. i. model predictions and comparison with diffusion theory,” IEEE Trans. Biomed. Eng. 36(12), 1162–1168 (1989). [CrossRef]
21. S. Flock, B. Wilson, and M. Patterson, “Monte carlo modeling of light propagation in highly scattering tissues. ii. comparison with measurements in phantoms,” IEEE Trans. Biomed. Eng. 36(12), 1169–1173 (1989). [CrossRef]
22. M. Keijzer, “Light transport for medical laser treatments,” Ph.D. thesis, Technische Univ., Delft (Netherlands). (1993).
23. M. Keijzer, S. L. Jacques, S. A. Prahl, and A. J. Welch, “Light distributions in artery tissue: Monte Carlo simulations for finite-diameter laser beams,” Lasers Surg. Med. 9(2), 148–154 (1989). [CrossRef]
24. S. A. Prahl, “A Monte Carlo model of light propagation in tissue,” in Dosimetry of Laser Radiation in Medicine and Biology, vol. 10305G. J. Mueller, D. H. Sliney, and R. F. Potter, eds., International Society for Optics and Photonics (SPIE, 1989), pp. 105–114.
25. L. C. Henyey and J. Greenstein, “Diffuse radiation in the galaxy,” Astrophys. J. 93, 70–83 (1941). [CrossRef]
26. F. Bevilacqua, D. Piguet, P. Marquet, J. D. Gross, B. J. Tromberg, and C. Depeursinge, “In vivo local determination of tissue optical properties: applications to human brain,” Appl. Opt. 38(22), 4939–4950 (1999). [CrossRef]
27. M. Ewerlöf, E. G. Salerud, T. Strömberg, and M. Larsson, “Estimation of skin microcirculatory hemoglobin oxygen saturation and red blood cell tissue fraction using a multispectral snapshot imaging system: a validation study,” J. Biomed. Opt. 26(2), 026002 (2021). [CrossRef]
28. L. Lorenz, “Lysbevægelsen i og uden for en af plane Lysbølger belyst Kugle,” Vidensk. Selk. Skr. 6, 1–62 (1890).
29. G. Mie, “Beiträge zur Optik trüber Medien, speziell kolloidaler Metallösungen,” Ann. Phys. 330(3), 377–445 (1908). [CrossRef]
30. P. Donelli, P. Bruscaglioni, A. Ismaelli, and G. Zaccanti, “Experimental validation of a Monte Carlo procedure for the evaluation of the effect of a turbid medium on the point spread function of an optical system,” J. Mod. Opt. 38(11), 2189–2201 (1991). [CrossRef]
31. L. Wang, S. L. Jacques, and L. Zheng, “Mcml—monte carlo modeling of light transport in multi-layered tissues,” Comput. Methods Programs Biomed. 47(2), 131–146 (1995). [CrossRef]
32. S. L. Jacques, “Coupling 3d monte carlo light transport in optically heterogeneous tissues to photoacoustic signal generation,” Photoacoustics 2(4), 137–142 (2014). [CrossRef]
33. E. Alerstam, W. C. Y. Lo, T. D. Han, J. Rose, S. Andersson-Engels, and L. Lilge, “Next-generation acceleration and code optimization for light transport in turbid media using gpus,” Biomed. Opt. Express 1(2), 658–675 (2010). [CrossRef]
34. A. Doronin and I. Meglinski, “Online object oriented monte carlo computational tool for the needs of biomedical optics,” Biomed. Opt. Express 2(9), 2461–2469 (2011). [CrossRef]
35. T. Young-Schultz, S. Brown, L. Lilge, and V. Betz, “Fullmontecuda: a fast, flexible, and accurate gpu-accelerated monte carlo simulator for light propagation in turbid media,” Biomed. Opt. Express 10(9), 4711–4726 (2019). [CrossRef]
36. X. Ma, J. Lu, R. Brock, K. Jacobs, P. Yang, and X. Hu, “Determination of complex refractive index of polystyrene microspheres from 370 to 1610 nm,” Phys. Med. Biol. 48(24), 4165–4172 (2003). [CrossRef]
37. D. Frantz, “Experimental Monte Carlo validation data (0.5 μm polystyrene spheres phantoms),” figshare: Version 1, 2021, https://doi.org/10.6084/m9.figshare.17152868.
38. D. Frantz, “Experimental Monte Carlo validation data (2 μm polystyrene spheres phantoms),” figshare: Version 1, 2021, https://doi.org/10.6084/m9.figshare.17152865.
39. D. Frantz, “Experimental Monte Carlo validation data (5 μm polystyrene spheres phantoms),” figshare: Version 1, 2021, https://doi.org/10.6084/m9.figshare.17152871.
40. J. E. Greivenkamp, Field guide to geometrical optics., SPIE field guides: FG01 (SPIE Optical Engineering Press, 2004).
41. D. F. Swinehart, “The beer-lambert law,” J. Chem. Educ. 39(7), 333 (1962). [CrossRef]
