key: cord-0162542-f7cy24jo
authors: McMillan, Lewis; Bruce, Graham D.; Dholakia, Kishan
title: Meshless Monte Carlo Radiation Transfer Method for Curved Geometries using Signed Distance Functions
date: 2021-12-15
journal: nan
DOI: nan
sha: 149d02ed34cf70091c18cc7f490c6246c2a84cde
doc_id: 162542
cord_uid: f7cy24jo

Significance: Monte Carlo radiation transfer (MCRT) is the gold standard of modeling light transport in turbid media. Typical MCRT models use voxels or meshes to approximate experimental geometry. A voxel based geometry does not allow for the accurate modeling of smooth curved surfaces, such as may be found in biological systems or food and drink packaging. Aim: We present our algorithm which we term signedMCRT (sMCRT), a new geometry-based method which uses signed distance functions (SDF) to represent the geometry of the model. SDFs are capable of modeling smooth curved surfaces accurately whilst also modeling complex geometries. Approach: We show that using SDFs to represent the problem's geometry is more accurate and can be faster than voxel based methods. sMCRT, can easily be incorporated into existing voxel based models. Results: sMCRT is validated against theoretical expressions, and other voxel based MCRT codes. We show that sMCRT can accurately model arbitrary complex geometries such as microvascular vessel network using SDFs. In comparison to the current state-of-the-art in MCRT methods specifically for curved surfaces, sMCRT is up-to forty-five times more accurate. Conclusions: sMCRT is a highly accurate, fast MCRT method that outperforms comparable voxel based models due to its ability to model smooth curved surfaces. sMCRT is up-to three times faster than a voxel model for equivalent scenarios. sMCRT is publicly available at https://github.com/lewisfish/signedMCRT

The modeling of light transport is important to our understanding of how light interacts with turbid media. It allows us to make predictions of the viability of treatment modalities [1, 2] , simulate the behavior of complex shaped light in highly scattering media [3] , retrieve images of objects in highly scattering media [4] and optimize light sensors in the food and drink industry [5] among other applications. The current "gold standard" of modeling light transport in turbid media is the Monte Carlo radiative transfer (MCRT) method. MCRT can model light transport in arbitrary 3D geometries and model several micro-physics phenomena such as Raman scattering [6, 7] , fluorescence [8, 9] , polarization [10, 11, 12] , and has been applied to problems ranging from light propagation in dynamic fluid systems [13, 14] to simulating thermal gradients in illuminated tissue [15, 16] .

In order to simulate the transport of light through a medium, the geometry of the problem must modeled. Most Monte Carlo codes rely on voxels [17, 18, 19] , or meshes [20, 21] to approximate the geometry of the problem. Voxel models are only suitable for the simplest problems that do not require accurate treatment of curved surfaces, due to their cubic nature [22] . Curved surfaces arise in many problems where MCRT maybe used, for example the propagation of light through: an optical system, the anatomy of animals or humans such as the brain or vascular networks, among many other possible examples. In contrast, mesh based models can accurately treat curved surfaces but usually require specialist software to create high quality meshes which can also be computationally expensive to produce, and MCRT codes require extensive software engineering to incorporate meshes in a computationally efficient manner.

A number of previous methods have been introduced to tackle the problem of smooth surfaces in voxel based models. Tran and Jacques preprocessed the voxels to determine where the material interfaces are, and computed surface normals for each voxel which can then be smoothed via interpolation to create curved surfaces [23] . Whilst this method, on the whole, accurately models curved surfaces, it can have an increased memory footprint, and is more computationally intensive. Alternatively, implicit surfaces can be defined using mathematical formulae [24, 25] . However, while this method of defining mathematical surfaces allows the modeling of smooth surfaces, it has the drawback that each surface needs an accompanying intersection and surface normal routine which can be computationally costly to evaluate and increases the workload on the programmer.

In this work, we present a novel Monte Carlo radiative transfer model where we eschew the common voxel or mesh based approaches for an approach based upon signed distance functions (SDFs), which we call signedMCRT (sMCRT). SDFs have been commonly used to define implicit surfaces in computer graphics [26] , video games [27] , and computer vision [28] . Recently, there has been considerable interest in using neural networks to define SDFs from point clouds and meshes. This interest has been led by computer graphics and deep learning researchers, looking for memory efficient representations of meshes and/or point clouds at high spatial resolutions [28, 29, 30] .

SDFs allow the efficient transport of photon packets through the modeled geometry using sphere tracing, which is faster, in most cases, compared to traditional ray tracing methods used in MCRT [26] . SDFs, while being similar to the approach of mathematical surfaces, do not need individual intersection routines as they are naturally included in the SDF definition. Moreover, we can use numerical differentiation to provide the surface normals. SDFs also require little effort to incorporate into existing voxel based codes, requiring tweaks to the optical depth integration and geometry initialization routines. These features of SDFs make them an attractive alternative to voxel models.

We show in this paper that sMCRT is more malleable and accurate the traditional voxel models for curved surfaces and more accurate than current solutions (by up-to 45 times more accurate) to the voxel discretization problem. Moreover, we demonstrate that sMCRT is faster (for certain problems). MCRT can be faster, or slower than traditional voxel based methods depending on the output required. When MCRT is used to retrieve the intensity at each point of a volume of interest, there is fractional increase in runtime of sMCRT (up-to 0.6 times slower) when compared with a traditional voxel based code. However, for problems where the desired output is only the light distribution at a single plane of interest, then sMCRT is fractionally faster than the traditional voxel-based methods (up-to 3 times faster).

The radiation transfer equation (RTE) describes the transfer of energy in a medium. However, analytical solutions for the RTE only exist for simple geometries. Therefore, numerical methods such as the diffusion method [31] or the Monte Carlo radiation transfer method (MCRT) must be used to compute a solution.

MCRT uses interaction probabilities and probability distribution functions that describe the physics of light transport, to model light transport though turbid and non-turbid media. Each photon is propagated a distance τ /µ t , where τ is the optical distance [-] and µ t [cm −1 ] is the extinction coefficient, before it interacts with the medium. The value signedMCRT A PREPRINT of τ is sampled from the probability distribution function for the mean free path of a photon using the Monte Carlo method [17] , as shown in Equation 1, where ξ is a random number drawn in the range [0, 1]. For voxel models the extinction coefficient can differ from voxel to voxel. MCRT is highly accurate (as the number of photons simulated tends to infinity) and can be used to model the light transport in arbitrary 3D geometries provided these geometries can be modeled.

The MCRT code presented in this work is broadly based upon previous MCRT codes used in various astronomical, medical, and bio-photonics applications [3, 32, 33, 34] . We use the same routines for releasing photons, input/output, scattering, random number generation, and helper routines. What differs in this work is the optical depth integration routine, and the geometry modeling method, which are accomplished by the use of signed distance functions.

Signed distance functions determine the distance from a point p, to the boundary of a specified shape. The function returns a positive value if p is outside the boundary, and a negative value if inside the boundary. Formally, this can be described using level set representation. In level set representation contours are modeled at the zero-level set (φ = 0) of a function defined in a higher dimension. Let Φ : Ω → R 3 be a Lipchitz function that refers to a level set representation for a given shape S [35] then:

An example of an SDF is shown in Equation 4 for a sphere, where r is the radius of the sphere, and p is the position of a photon:

SDFs can easily be translated, rotated, twisted and scaled among many other operations. Constructive solid geometry (CSG) operations such as union, intersection and difference can also be used on the SDFs. Figure 1 shows a subset of shapes, and possible operations on SDFs [36, 37] .

To incorporate SDFs into a pre-existing voxel based MCRT code, requires only relativity small adjustments: modifications to the geometry initialization routine, and to the optical depth integration routine. An overview of the complete MCRT algorithm is shown in the left panel of Figure 2 .

To create the geometry in voxel based models, each voxel is independently assigned a set of optical properties (scattering and absorption coefficients, refractive index, and anisotropy g value). In sMCRT the geometry is initialized by selecting the functional form, size, and location of SDF(s) required to model the problem, applying any CSG operations required to generate more complex shapes, and finally setting the optical properties for each SDF. Each SDF has its own set of optical properties, which include scattering and absorption coefficient, refractive index and the anisotropy g value. We then create a bounding box around all the SDFs, which gives us a simulation volume of interest.

In voxel based MCRT codes, each photon packet is randomly ascribed a specific optical path length that it travels before an interaction, such as scattering or absorption, according to Equation 1 and is scaled by µ t (µ t can be different for each voxel). The photon packet is then propagated through the voxel grid using ray tracing until it reaches that interaction point or leaves the voxel grid.

In our SDF based MCRT algorithm, the first step in the SDF optical depth integration routine is the same as in the voxel case, i.e. randomly assign an optical depth. As before, this is calculated using Equation 1 and µ t can be different for each SDF. The next step is to acquire the distance from the current position of the photon packet to the nearest boundary. This is computed by using the SDFs to calculate the distance to each boundary in the modeled geometry, and taking the minimum value (d sdf ). This process is called sphere tracing [36] , and is illustrated in Figure 3 .

signedMCRT A PREPRINT Figure 1 : Several examples surfaces that can be created by SDFs, rendered in Blender. For illustrative purposes, SDFs are voxelized in sMCRT then transformed to a mesh using Skimage's [38] marching cubes algorithm and then rendered using Blender [39] . The left two panels show a subset of basic shapes calculated using SDFs (a to d and j to m). The right two panels show a subset of possible operations on SDFs: smooth (e) and non-smooth union (f), intersection (g), subtraction (h), repetition (i), displacement (n), elongation (o), bend (p), and twist (q).

If the remaining optical depth for the photon packet is less than d sdf , the photon packet undergoes some interaction, and the optical depth integration routine restarts. If the optical depth is not reached, then we move the full distance d sdf , and then recalculate the distances to all boundaries. If the SDF of the bounding box returns a positive value we are outside the volume of interest, so we terminate the packet and start a new packet.

If the SDF for the bounding box returns a negative value, we then check if the smallest distance, d sdf , is less than some threshold, δ. In this case, the photon packet is on a boundary so we need to check if there is a change in refractive index.

If there is a change in refractive index we calculate the Fresnel coefficients and the surface normals, then reflect or refract the photon packet. If d sdf is larger than δ, and all distances to the SDFs are not positive then we start this whole process again until one of the exit conditions has been met. The above algorithm is shown in the right panel of Figure 2 .

To ensure that our novel SDF-based geometry method works accurately, we validate our algorithm against a theoretical expression, and two other MCRT codes. All simulations are fully parallelized with OpenMP and were run on a workstation with an AMD Ryzen 9 3950X 16-Core Processor with 64 GB RAM utilizing the full 32 threads available.

We first compare sMCRT's accuracy by computing the average number of scattering events occurring to a photon in an isotropic sphere [40] .

For a photon's random walk from the center to the edge of a uniformly scattering sphere of radius r, the average number of scatterings that take place can be written as (see Supporting Information):

To compare Equation 5 to sMCRT, we model a sphere of radius 0.5 cm, vary the radial optical depth between 0.1 and 100, and release 10 million photons isotropically from its center. The agreement of the code and analytical expression can be seen in Figure 4 .

We also validated sMCRT against our previous voxel based MCRT code [3, 15, 41] , and S. Jacques et al. MCRT code [42] . We validate against S. Jacques code as it incorporates all the relevant physics we need in an MCRT code; scattering, absorption, and refractive index mismatches. For this validation, the medium is set up as a semi-infinite slab and light is uniformly incident on the surface of the slab (negative z direction) and propagates until it is absorbed or escapes via the top surface (positive z direction). We then fit against Equation 6 to compare between codes.

Here Ψ(z) is the penetration of the incident light or equivalently the fluence rate [W cm −2 ], Ψ 0 is a normalization constant [W cm −2 ], C n and k n are fitted coefficients [-] , and δ is the optical penetration depth, defined as δ = 1/ 3µ a (µ a + µ s (1 − g)) [cm]. The optical properties for the slab are shown in Table 1 , where we use the Henyey-Greenstein phase function [43] with a g of 0.9, and we model two wavelengths in separate simulations. The refractive index for the medium was set to 1.38 to mimic the rat skin used in S. Jacques code [42] , and for the surrounding medium a refractive index of 1.0 was set, to mimic air. As evidenced in Figure 5 , sMCRT gives a better match to the results in S. Jacques et al. work than our previous voxel model. An exact match is not possible, due to the difference in the code underlying workings such as cylindrical fluence bin shape used by S. Jacques et al. versus our rectangular bin shape. sMCRT is more accurate than the voxel model at 420 nm as it incorporates the left hand-side of the peak near the surface whereas the voxel model levels off after the peak fluence value.

In this section we compare the accuracy and speed of a voxel based MCRT code to that of our new SDF based MCRT code, sMCRT. We devise two tests based upon speed and accuracy.

The first test is a simple test of modeling a smooth surface: for this we use the test case in Tran and Jacques surface normal approach [23] . They model a glass sphere (refractive index 1.46) and let a beam of light, with a waist of 0.3 cm, be incident on the sphere (in the negative z direction). Figure 6 shows a slice of the fluence through the sphere for sMCRT, a voxel based MCRT code, and Tran and Jacques code. The difference in accuracy between the surface normal approach and sMCRT can be seen in the lower row of Figure 6 , which shows a line profile across each of the reflected and refracted beams for the two approaches. The beam profiles produced by sMCRT exhibit a reduction in the root mean square deviation from the expected beam profile by a factor of 10 for the reflected beam and by a factor of 45 for the refracted beam, when compared to the surface normal approach.

sMCRT also outperforms the surface normal approach on speed for this test case. sMCRT took on average 61 seconds to run ∼ 27 million photons compared to 250 seconds for the surface normal approach. The speedup factor here is 4. Both codes were run on the full 32 threads available on our workstation.

The second test is a comparison of speed between the SDF and voxel models. Here we use the same setup from the validation section, an isotropic scattering sphere which releases photons isotropically from its center. We vary the scattering coefficient and time the simulation for 10 million photons. As shown in Figure 7 , the SDF model is up to 3 times faster when run on the same workstation.

This however may not be representative of the true speed difference between the SDF and voxel models in real use scenarios. In another test, we calculated the average number of scatterings per photon, which requires no fluence computation and thus no grid to record fluence. The most common use of MCRT codes is to compute the fluence throughout the simulated medium. The voxel model provides this with little extra computation. However, sMCRT needs no such grid, so the introduction of a fluence counting grid will slow down the simulation. We therefore, test the two models whilst computing the fluence for two different scenarios: the S. Jacques validations case, and the aforementioned isotropic sphere test. Table 2 summarizes the results of these tests. For the S. Jacques et al. validation case, we use the same geometry as before and a fluence grid of size 200 × 200 × 400 voxels for both the SDF and voxel models. We run 10 million photons and use 32 threads on our workstation for both models. The isotropic sphere test was run for 1 million photons and a grid of 200 3 voxels for both the SDF and voxel models. The radial optical depth of the sphere was set to 100.

The results show that the sMCRT is slower than the voxel model for both test cases by a factor of ≈ 1.5. It should be noted that sMCRT is a new code, whereas the voxel model has had several years of development and optimizations applied to it, meaning that sMCRT could reach the same level of performance as the voxel model with a similar effort applied. sMCRT does however, give more accurate results for both of these problems at the expense of speed. 

SDFs can be used in place of current voxel based geometries as they offer increased accuracy in the modeling of complex geometries, as well as the accurate treatment of smooth surfaces. We have recently used sMCRT in our simulations of Raman spectroscopy of alcoholic beverages through glass bottles [44] . Accurate modeling of a bottle in a voxel geometry was not possible, as this did not allow accurate reflection/refractions at the glass/air and glass/alcohol interfaces. The lack of these accurate reflections meant a massive decrease in the accuracy of the simulation of the light distribution in the bottle, as the bottle acts as a cylindrical lens on the light, see SI 8. Thus far we have shown that sMCRT can model simple geometries, so to demonstrate that sMCRT can model complex shapes, we model a blood vessel network embedded in tissue, and show that sMCRT can model arbitrary SDF generated neural networks, such as DeepSDF or SIREN, see SI Figure 1 .

The vessels are a 3D synthetic microvascular network from data published in [45] and preprocessed by Yuan et al. [46] . We model the slab of tissue using a box SDF (second shape in bottom left panel of Figure 1 ), and the vessels as a collection of capsule SDFs (third shape in top left panel of Figure 1 ) which are then joined together using the CSG operator union (bottom left operation in top right panel of Figure 1 ). The simulation volume is 326 × 305 × 611 µm 3 and we use a voxel grid of 400 3 to record the fluence. The optical properties of the slab of tissue and the vessels are taken from [18] , and are shown in Table 3 . The slab is uniformly illuminated on its top surface by 10 million photons, which are allowed to propagate until they are absorbed or leave the simulated medium. Figure 9 shows the fluence on the vessel network, and slices of fluence through the tissue slab.

We have shown a novel meshless, geometrical method for Monte Carlo radiation transport, using signed distance functions. SDF based models achieve higher accuracy than voxel based models, particularly for modeling smooth Therefore, SDFs offer an improvement on voxel based models, whilst retaining high accuracy and ease of implementation. However, there are a number of potential downsides to using SDFs.

In certain configurations the number of steps needed to be taken by a photon packet can be extremely large, see SI Figure 3 . This occurs when the photon is approximately parallel to a surface whilst the distance between the photon and the surface is small. Recent work has been undertaken to alleviate this problem. This includes segment tracing which accelerates the sphere tracing method by improving the marching step computation and enhanced sphere tracing which uses an over-relaxation-based method for accelerating sphere tracing [47, 48] .

Large collections of SDFs can also cause massive slow downs due having to evaluate every SDF each time the photon needs to be moved. This is analogous to the same issue in Monte Carlo models which use triangular meshes. As in the triangular mesh case, this can be diminished by using a space-partitioning data structure.

The final potential downside is the combination of multiple CSG operations can lead to non-bounded SDFs. Nonbounded SDFs can pose a problem in terms of accuracy and speed [49] . In terms of speed, non-bounded SDFs only give a conservative distance to the surface, resulting in more SDF evaluations which can cause computational slow down. The accuracy problem only affects the computation of surface normals, and is therefore only applicable at refractive index interfaces.

Despite these potential downsides, we envision MCRT codes using SDFs to model the geometry to probe problems such as: effect of skin color on pulse oximtry accuracy, fluence calculation of droplets with viral loads such as Covid-19, and accurate simulations of light propagation in fruit. In these problems, using SDF over voxels would allow accurate modeling of curved surfaces allowing better accuracy in the simulations. 1 Derivation of Average number of scattering in an isotropic sphere For a photon's random walk from the center to the edge of a uniformly scattering sphere of radius r. Consider the net displacement for a single photon from a starting point p after N mean free paths is:

The mean square displacement traveled by a photon (l * ) is thus: l 2 * = R 2 = r 2 1 + r 2 2 + ... + r 2 N + 2 r 1 · r 2 + ...

where 2 r 2 · r 2 = 2 |r 1 ||r 2 | cos δ

For isotropic scattering cos δ = 0 therefore all cross terms vanish. As each term involving the square of the displacement averages to the mean square of the free path of a photon, therefore

Then multiplying both sides by µ 2 , the scattering coefficient, gives an expression for the average number of scatterings for optically thick media. By definition, the scattering coefficient multiplied by the radius of the sphere is τ max :

Then using the following identity

we arrive at an expression for the average number of scatterings as a function of maximum optical depth for an isotropic sphere in the optically thick limit:

For the optically thin limit, the number of scatterings is small, on the order of 1 − e −τ ≈ τ .

2 Examples of sMCRT on complex geometries SI Figure 1 : Example of an arbitrary SDF generated via a neural SDF method (SIREN [1] ). Image shows 3D fluence inside the Stanford bunny [2] for a point source inside the Bunny. Fluence outside the bunny has been removed for clarity.

SI Figure 2 : Fluence on the University of St Andrews crest. Crest is converted from a simplified Scalable Vector Graphic (SVG) of the crest, and translated into 2D line segments. Line segments are then extruded in the z-axis and assigned optical properties, where µ a =0.1 cm −1 , µ s =10.0 cm −1 , g=0.9, and n=1.5. Left shows a slice of fluence in the x-y plane in the middle of the simulated medium. Photons are incident on the crest's surface (into the page). Right shows the simplified SVG where each line has a different color. SVG was simplified by removing some elements, and converting Bézier curves into line segments using Inkscape. 

Monte Carlo modelling of daylight activated photodynamic therapy

Could psoralen plus ultraviolet A1 ("PUVA1") work? Depth penetration achieved by phototherapy lamps

Imaging in thick samples, a phased Monte Carlo radiation transfer algorithm

Computational time-of-flight diffuse optical tomography

Monte Carlo simulation for quantification of light transport features in apples

Monte Carlo model of spatially offset Raman spectroscopy for breast tumor margin analysis

Monte Carlo simulation of in vivo Raman spectral measurements of human skin with a multi-layered tissue optical model

Experimental validation of Monte Carlo modeling of fluorescence in tissues in the UV-visible spectrum

Propagation of fluorescent light

Polarised light sheet tomography

Three Monte Carlo programs of polarized light transport into scattering media: part I

Three Monte Carlo programs of polarized light transport into scattering media: part II

The Monte Carlo photoionization and moving-mesh radiation hydrodynamics code CMacIonize

The TORUS radiation transfer code

Development of a predictive Monte Carlo radiative transfer model for ablative fractional skin lasers

Monte Carlo simulations of heat deposition during photothermal skin cancer therapy using nanoparticles

Mcml-Monte Carlo modeling of light transport in multi-layered tissues

MCmatlab: an open-source, user-friendly, MATLABintegrated three-dimensional Monte Carlo light transport solver with heat diffusion and tissue damage

Three dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult human head

Mesh-based Monte Carlo method using fast ray-tracing in Plücker coordinates

Mesh-based Monte Carlo code for fluorescence modeling in complex tissues with irregular boundaries

Light transport in tissue by 3d Monte Carlo: influence of boundary voxelization

Modeling voxel-based Monte Carlo light transport with curved and oblique boundary surfaces

Monte Carlo simulation of light transport in turbid medium with embedded object-spherical, cylindrical, ellipsoidal, or cuboidal objects embedded within multilayered tissues

Monte Carlo simulation of radiation transport in human skin with rigorous treatment of curved tissue boundaries

Ray tracing deterministic 3-D fractals

Learning from failure: a survey of promising, unconventional and mostly abandoned renderers for 'Dreams PS4', a geometrically dense, painterly UGC game

Deepsdf: Learning continuous signed distance functions for shape representation

Implicit neural representations with periodic activation functions

Neural geometric level of detail: Real-time rendering with implicit 3d shapes

Diffusion theory of light transport

A model for the scattered light contribution and polarization of the diffuse Hα galactic background

Quantifying direct DNA damage in the basal layer of skin exposed to uv radiation from sunbeds

Depth penetration of light into skin as a function of wavelength from 200 nm-1000 nm

Matching distance functions: A shape-to-area variational approach for global-to-local registration

Sphere tracing: A geometric method for the antialiased ray tracing of implicit surfaces

Distance functions

scikit-image: image processing in python

Blender -a 3D modelling and rendering package. Blender Foundation

Radiative processes in astrophysics

Advanced 3D Monte Carlo algorithms for biophotonic and medical applications

How photobleaching affects dosimetry and fluorescence monitoring of PDT in turbid media

Diffuse radiation in the galaxy

To focus-match or not to focus-match inverse spatially offset Raman spectroscopy: a question of light penetration

Deepvesselnet: Vessel segmentation, centerline prediction, and bifurcation detection in 3-D angiographic volumes

Light transport modeling in highly complex tissues using the implicit mesh-based Monte Carlo algorithm

Segment tracing using local Lipschitz bounds

Enhanced sphere tracing

Zippered polygon meshes from range images

Implicit neural representations with periodic activation functions

Zippered polygon meshes from range images