1. Introduction

The volume, or bulk, scattering of light by materials is a crucial effect that needs to be considered in applications ranging from biomedical research to commercial lighting systems [1–4]. Typically, this effect is described by a model that relies on a small number of parameters. The scattering and absorption coefficient describe the likelihood of light scattering or absorption events occurring along the length of the optical path, while a phase function describes the angular probability distribution in which the light is scattered at each scattering event. The phase function can have any arbitrary shape, but for most materials a good balance between flexibility and accuracy can be obtained by using a phase function model, e.g. Henyey-Greenstein, Gegenbauer Kernel or Von-Mises [5, 6]. Choosing, for instance, the Henyey-Greenstein (HG) phase function model, which relies on a single parameter, it is possible to have a reliable volume scattering model for a specific material using only three parameters.

The accuracy of volume scattering simulations is heavily dependent on having accurate estimates for the volume scattering parameters. Typically, these are obtained using inverse, or fitting, methods. These methods iteratively test different candidate parameters by comparing simulated data to experimental measurements until a good match is found [1]. For instance, the inverse adding-doubling (IAD) method starts by using the adding-doubling method to simulate the total transmission and total reflection of a volume scattering sample. Then, it compares this simulated data to experimental data to identify the set of scattering parameters that provide the best match [7]. However, the scattering parameters obtained using only reflection and transmission information are usually not accurate enough to simulate more complex light scattering characteristics such as the intensity distribution. This was the motivation to extend the IAD method to fit to intensity information instead of using only total integrated quantities. This extension allowed retrieving volume scattering parameters that more accurately reproduced the angular patterns of light scattered from samples [8].

For lighting applications, it is important to be able to predict the luminance or radiance from light sources [9–11]. The radiance distribution includes more descriptive information about the bulk scattered light than the intensity distribution, as it describes both the spatial and angular distribution unlike the intensity, which only describes the angular distribution. This raised doubts regarding the intensity based IAD method’s ability of retrieving volume scattering parameters that could be used to reproduce both the intensity distribution and the radiance distribution. An in-depth simulation study tackled this question and showed that, using all the information available in the intensity distribution, the IAD method could accurately characterize samples in a way that could also reproduce the scattered radiance distribution. These results were obtained for samples with short to moderate optical path lengths, that is, samples with low to moderate absorption and scattering coefficients. However, as samples with longer optical path lengths were investigated, the method’s accuracy dropped; it was able to identify volume scattering parameters that reproduced the intensity information well, but the same parameters failed to describe the radiance distribution [12]. A possible explanation may lie in how differently the intensity and radiance distributions describe the scattered light. At longer optical path lengths, the descriptive power remaining in the intensity distribution may not be enough to accurately identify the true set of volume scattering parameters, while there may be more descriptive information available in the radiance distribution that could be used to obtain a more accurate characterization.

To investigate this, a Monte Carlo ray-tracer was implemented. Unlike the adding-doubling method, the ray-tracer can easily be adapted to simulate radiance distributions, instead of intensity distributions. This increased flexibility comes at the cost of much longer simulation times. This significant disadvantage was addressed by implementing the Monte Carlo method as a routine that exploits the parallel capabilities of modern graphics processing units (GPUs). With this implementation, it is possible to simulate the scattered radiance distribution of samples in a reasonable amount of time. Due to its efficiency, this implementation allowed to design and implement an inverse Monte Carlo (IMC) method to estimate the volume scattering properties using the radiance distribution of the light scattered by samples.

In this report, the implemented GPU based ray-tracer and the inverse algorithm are discussed and their performance is analysed. More importantly, the descriptive power differences between using radiance or intensity to identify accurate volume scattering parameters is studied for samples with moderate to long optical path lengths.

2. Methods

In this work, a simulation study is performed by considering only virtual samples, i.e. samples with pre-established volume scattering parameters. As the volume scattering parameters of these virtual samples are known, it is possible to investigate the performance of the inverse method in detail, which would be infeasible in a typical experimental scenario. The volume scattering model used throughout this paper consists of three parameters. The absorption coefficient μ_a, the scattering coefficient μ_s and the anisotropy parameter g of the Henyey-Greenstein phase function model. A sample can also be described using other variables derived from these parameters, such as the single-scattering albedo a = μ_s/(μ_s + μ_a) and the optical path length τ = l (μ_s + μ_a), for a material of thickness l. The optical path length is relevant to this work, as it describes how many interactions occur, on average, between a photon packet and the volume scattering material. These scattering parameters along with the sample’s geometry and its refractive index are the only necessary inputs for a ray-tracer to obtain the scattered light’s radiance distribution for a specific incoming light distribution. For all samples tested, the refractive index was kept fixed at n = 1.4, the sample’s thickness at l = 1 mm and, to facilitate comparison with the adding-doubling method, very large lateral dimensions were used with width w ≥ 500 mm. Only pencil beams that are incident perpendicularly to the surface of the sample are considered.

In this section, the implementation of the Monte Carlo ray-tracer on a GPU is briefly discussed. This method and the information provided with it is then analysed to understand its structure and properties. This analysis is used to develop an optimization method to solve the inverse problem of finding accurate estimates for the scattering parameters using the radiance distribution. Finally, the differences between using radiance and intensity distributions to obtain accurate volume scattering parameters is studied, along with the impact of experimental noise on the conclusions obtained.

2.1. GPU accelerated Monte Carlo ray-tracer

The original Monte Carlo algorithm for photon transport in scattering media was used as the basis for this GPU implementation. The method’s details can be found elsewhere [13], thus only the main differences in relation to this implementation are discussed. Unlike in the original Monte Carlo algorithm and the adding-doubling method, the sample’s geometry constraints are relaxed to better fit to typical experimental scenarios. The samples are not required to be fully planar nor infinitely wide [14]. Moreover, in lighting, characterization samples are generally restricted to a single homogeneous layer, which helps to simplify the method. Also, because the goal of this ray-tracing method is to produce the distributions of the light scattered by the sample, the internal distributions, such as the fluence rate, are not calculated. The implemented GPU ray-tracer can also simulate multiple wavelengths on each run, use arbitrary phase function shapes and handle luminescent materials. For this study however, only a single wavelength is considered and only volume scattering materials are analysed.

The output of the ray-tracing method is a ray-set that contains the necessary information to calculate, for instance, the intensity distribution or the total transmission of the light scattered by the sample. As this study focuses on radiance, that is the only quantity which is calculated from the dataset. The ray-set information is binned into a 4 dimensional structure; it bins the rays’ exit point at the sample’s surface (2D spatial data) and their azimuth and inclination (2D angular data). Each specific bin accumulates the flux of the rays with specific spatial and angular properties. This results in a 4D structure that encodes the spatial and angular distribution of the optical flux. To obtain a 4D radiance structure, it is necessary to scale each bin with its projected area and solid angle. That is, each bin corresponds to a certain area element on the emitting surface. At the same time, each bin also corresponds to light emitted within a specific solid angle. The corresponding projected area and solid angle are used to scale the accumulated flux at each bin, thus obtaining the required 4D discretized radiance structure.

The accuracy of this radiance distribution depends on the amount of rays traced and the sample’s volume scattering characteristics. A sample with a short optical path length, that is a sample with low absorption and scattering coefficients, requires tracing less rays to obtain an accurate simulation of the transmitted radiance distribution in comparison to a sample with a long optical path length. To avoid using heuristics to identify the adequate number of rays to be traced for each sample type all simulations use the same number of rays. The number of rays was determined by considering the simulated radiance of the sample with the highest absorption and scattering, i.e. the most challenging sample. The ray number was chosen so that the simulated distribution’s shape would be clearly discernible, while keeping the simulation time low. This guarantees that the simulated distributions for all other samples share, at least, the same quality. Increasing the number of bins, i.e. the resolution of the radiance distribution, also increases the noise within each bin. This trade-off between accuracy (low variance) and resolution forces a compromise, which in this work was settled by using 32 by 32 spatial bins, 16 inclination bins and 2 azimuth bins and by tracing 2³⁰ rays. Instead of using the full area of the sample, only the area corresponding to the emitting area is used, which corresponds to x, y ∈ [−2, 2] mm. The inclination angles bins are bound by θ ∈ [0, 60]° and the azimuth angles used are ϕ ∈ [0, 360]°

A routine that performs the ray-tracing and ray binning was implemented using CUDA C++, resulting in a very short simulation time when compared with typical ray-tracers used for illumination applications. To evaluate the method’s accuracy and speed, the GPU method was validated against LightTools (Synopsys, Inc., CA, USA), one of the most popular commercial ray-tracers available. A grid of virtual samples spanning short and long optical path lengths was generated with μ_s ∈ [5.0, 30.0]₇ mm⁻¹, μ_a ∈ [0.1, 5.0]₇ mm⁻¹ and g ∈ [0.5, 0.9]₃. For each set of parameters, the radiance distribution was simulated with both LightTools and the GPU ray-tracer. In these simulations only 2²⁵ rays were traced because LightTools’ simulation time and computational space requirements were too demanding to allow using higher values. The cosine distance was used to calculate the error between x, the radiance distribution simulated with LightTools, and y, the radiance distribution simulated with the GPU ray-tracer:

The cosine distance estimates the difference in shapes between the distributions, instead of their absolute difference. It is obtained from the cosine similarity [15] (also sometimes referred to normalized cross-correlation) by applying the arccos to reverse its domain. When the distributions are closely correlated, that is, have very similar shapes, the cosine distance will be close to 0. For completely uncorrelated distributions, the cosine distance will have a value of 0.5.

The volumetric error landscape is averaged along the phase function g axis, resulting in Fig. 1. The image shows that the error is low throughout the domain, demonstrating that both methods produce essentially the same radiance distributions. As a reference, running the same LightTools simulation using different random number streams results in cosine distance values ranging from 1 × 10⁻³ to 4 × 10⁻³. The computational speed performance is shown in Fig. 2. As expected, the GPU ray-tracer shows a speed-up compared to LightTools that increases with the number of rays traced and ranges from 25× to 190×, for the domain tested. This significant computational speed improvement is crucial for the GPU ray-tracer to be a good candidate for an inverse routine, which typically requires performing multiple iterations until converging on a solution.

 

Fig. 1 Cosine distance between the radiance distribution of samples simulated with LightTools and with the GPU ray-tracer. Left and right inserts show cross-sections (images) of the radiance distributions simulated with LightTools (top) and the GPU ray-tracer (bottom) for two specific samples.

Download Full Size | PDF

 

Fig. 2 Time taken over amount of rays traced for both the LightTools ray-tracer and the GPU ray-tracer for a long optical path length sample shown in (a) a logarithmic scale and (b) a linear scale. Results shown are the average of 5 runs, including all overheads, using a workstation with an Intel Xeon Processor E5-2630 v3 CPU and a NVIDIA Titan X (Maxwell) GPU.

Download Full Size | PDF

2.2. Inverse Monte Carlo - radiance fitting

The inverse procedure, similarly to IAD, is an iterative brute-force method. It receives the sample’s scattered radiance distribution as input and, in each iteration, chooses a test set of volume scattering parameters and simulates the corresponding scattered radiance distribution. It then compares the original and test distributions and, if there is no agreement, it chooses a new set of volume scattering parameters. This process is repeated until a good agreement is found. The search of the best set of parameters can essentially be framed as a numerical optimization problem whose objective or merit function is simply the error between the original radiance distribution and the test radiance distribution. As such, it can be solved with one of the many tools available.

Understanding the structure of the objective function domain is helpful in choosing an optimization tool that can best exploit it. This is the first step in developing the inverse procedure. This was done by building a grid of linearly spaced volume scattering parameters with μ_s ∈ [20.00, 35.00]₄₁ mm⁻¹, μ_a ∈ [1.00, 8.00]₄₁ mm⁻¹ and g ∈ [0.80, 0.90]₄₁. For each individual set of parameters (each sample) the radiance distribution was simulated using the GPU ray-tracer. An objective function domain was calculated by considering one sample as being the objective sample and comparing its scattered radiance distribution to that of all other samples in the grid using the cosine distance. Instead of depicting the full volumetric objective function domain, only three of its cross-sections intersecting at the set of volume scattering parameters corresponding to the objective sample are shown. An illustration of these cross-sections is shown in Fig. 3. The objective function domain was calculated for two samples at opposite ends of the volumetric domain: a short optical path length sample with μ_s = 21.69 mm⁻¹, μ_a = 1.79 mm⁻¹ and g = 0.81 and a long optical path length sample with μ_s = 33.69 mm⁻¹, μ_a = 7.30 mm⁻¹ and g = 0.89.

 

Fig. 3 Cross-sections of the objective function domain calculated for a short optical path length sample. The three cross-sections intersect at the set of volume scattering parameters of the objective sample.

Download Full Size | PDF

The different objective function domain cross-sections are shown in Fig. 4. For easier discrimination of the domain’s structure, the fourth root is used to scale all images shown. For the cross-sections perpendicular to either the μ_s or the g axis, there is a clear observable global minimum and a gradient towards it, for both the long and the short optical path length cases. Yet, the cross-sections that are perpendicular to μ_a show a non-convex domain, making it slightly less amenable to optimization. When the objective function domain is analysed at higher resolutions, it becomes clear that the domain is very noisy. This is a consequence of using a stochastic ray-tracer to simulate the radiance distributions. While this can be alleviated by tracing more rays, the simulation time also increases, which is undesirable. This implies that gradient based methods are ill advised, as discrete approximations to the gradient will be very noisy. A more attractive strategy would be to use a method based on a surrogate model, commonly used for expensive black-box optimization problems [16]. However, initial tests with the ray-tracer provided an important result that was key in identifying a good optimization method. By tracing an increasingly lower amount of rays and calculating the objective function cross-sections it became clear that globally these cross-sections retain their shape, as shown in Fig. 5. That is, the objective function domains calculated using the cosine distance between radiance distributions are robust to the statistical noise present in ray-tracer simulations. This results from the fact that the cosine distance is more sensitive to the differences in shape than absolute value.

 

Fig. 4 Radiance based objective function cross-sections perpendicular to the (a,d) μ_a axis, (b,e) the μ_s axis and (c,f) the g axis. The red cross indicates the true solution and correspond to a sample with (a–c) a short optical path length and (d–f) a long optical path length.

Download Full Size | PDF

 

Fig. 5 Radiance based objective function cross-sections for the long optical path length sample. The cross-sections were calculated by tracing (a) 2²⁶ rays, (b) 2²⁸ rays and (c) 2³⁰ rays.

Download Full Size | PDF

This property, along with the short simulation time, allows using the GPU ray-tracer to simulate noisy, but still descriptive, volumetric domains in a short amount of time. These initial domains can be polled for a useful, albeit noisy, first guess of the real solution. Hence, the IMC method implemented consists of first simulating, by tracing a low amount of rays, the radiance distribution of all samples in a specific domain. Then, the cosine distance between these distributions and the original radiance distribution is calculated to identify a good first guess. Finally, this guess is given to a direct-search method (derivative free) to refine it. For this, the Mesh Adaptive Direct Search (MADS) method [17] was chosen, parametrized to perform a search in the initial iteration with the Nelder-Mead simplex method [18]. The hybrid NM-MADS method was implemented in Matlab (The MathWorks Inc., MA, USA).

2.3. Radiance versus intensity for fitting

The method described above allows estimating the volume scattering model parameters of a sample using the radiance distribution of the light scattered from it. While it can also calculate and use intensity information, for this purpose the adding-doubling method is more efficient. The AD method was used to investigate if using the scattered intensity information of long optical path length samples also allows extracting accurate volume scattering parameters. To investigate this, it is necessary to study the shape of the objective function domain, as was done for the radiance case, but with intensity information. The same sample domain is used and for each sample the intensity distribution is simulated using the adding-doubling method and compared to the intensity distribution of the sample established as the objective sample.

The adding-doubling method uses the same inputs as the GPU ray-tracer, with the exception of the angular resolution that is fixed at 256 angular points for stability and accuracy reasons. Unlike the ray-tracer, the adding-doubling method simulates intensity distributions without stochastic noise, so any objective function domain calculated with it should be accurate to machine precision. It is important to note that the comparison done is independent of using two different tools to simulate the intensity and radiance distributions. Any conclusion drawn is the result of how each distribution encodes the information about the light scattered by the sample.

2.4. Radiance versus intensity for fitting with noise

An analysis performed with an accurate simulation study usually corresponds to the perfect experimental scenario, that is, one that is free of experimental noise. To generalize any conclusions obtained in the simulation study to a realistic experimental context, the effect of experimental noise was included. A simple model for both thermal noise (Gaussian noise) and shot noise (Poisson noise) was used [19]. The noisy signal S′, for both radiance and intensity distributions, can be calculated from the original signal S by adding zero-mean, normally distributed random numbers r_N, scaled by specific thermal and shot noise factors f_G and f_P using the following equation:

The values for the factors were chosen to imitate typical noise levels in experimental scenarios. An example of the original and noise added signals is shown in Fig. 6 for a long optical path length sample. When the signal is very low the effect of the noise is significant, which can be clearly seen in the transmitted intensity calculated for the sample. Similarly, the cross-section of the transmitted radiance, shown in the same figure, shows a very distorted signal. When the signal amplitude is higher, like the reflected intensity signal shown in Fig. 6, the effect of the noise is more subtle. This is also the case when looking at low optical path length samples, where the signal amplitude is usually high. Because both the transmitted and reflected parts of the distributions are required for an accurate fit, the effect of the experimental noise cannot be discounted by disregarding the noisy part of the distribution and keeping only the noiseless part.

 

Fig. 6 Signal before and after adding artificial noise for a sample with μ_s = 31.40 mm⁻¹, μ_a = 6.60 mm⁻¹ and g = 0.89. In (a) the reflected and transmitted intensity are shown, with the noisy versions in black, while a cross-section of the transmitted radiance distribution is shown (b) before adding noise and (c) after adding noise.

Download Full Size | PDF

Using a noisy version of the original signal, for both radiance and intensity distributions, the objective function cross-sections for the sample domain that was previously used were calculated.

3. Results

In this section three main results are discussed. First, the accuracy of the proposed IMC method is analysed and discussed. Then, the impact of using either intensity or radiance distributions to estimate volume scattering parameters of samples is elaborated. The previously performed investigation of the information available in the radiance distribution, which was necessary to develop the inverse method, is used in this comparison. Finally, the influence of experimental noise on the comparison done between the radiance and the intensity distributions is discussed.

3.1. Inverse Monte Carlo

To investigate the performance of the proposed inverse Monte Carlo method, an extensive sample domain was used with μ_s ∈ [10.00, 40.00]₃₁ mm⁻¹, μ_a ∈ [1.00, 10.00]₃₁ mm⁻¹ and g ∈ [0.65, 0.95]₂₁. In this test, each point in the volumetric sample domain was used to simulate a radiance distribution which was then supplied to the IMC method to estimate a set of volume scattering parameters that could reproduce it. The error between the estimated and original radiance distribution is calculated using the cosine distance, resulting in a volumetric error domain. Similar to the procedure employed for the results shown in Fig. 4, the error is averaged along the phase function anisotropy factor g axis. Because the true volume scattering parameters are known, the error between the estimates and the original parameters is also calculated and shown in Fig. 7. To calculate this error, the mean percent error between the original and estimated parameters are calculated for each sample in the grid and also averaged along the g axis.

 

Fig. 7 Error analysis of the results from the inverse Monte Carlo method. In (a) the average cosine distance between original and fitted radiance distribution is shown, while in (b) the average percent error between original and estimated volume scattering parameters is presented.

Download Full Size | PDF

Overall the fitting error, that is, the error between the original and the best radiance distribution found, is below 1×10⁻³ – a low error value. This indicates that the inverse method is generally successful in finding a radiance distribution that matches well with the supplied radiance distribution. More importantly, the average error between the estimated and the original volume scattering parameters is also low, generally under 2%, and does not change significantly throughout the sample domain. This result demonstrates that the proposed IMC method is capable of extracting the sample’s volume scattering parameters accurately even when long optical path length samples are considered. Moreover, the fitting error has a strong correlation to the parameter error, which shows that the cosine distance is indeed a good metric for a radiance based IMC method.

3.2. Radiance versus intensity for fitting

The radiance based objective function cross-sections calculated in Fig. 4 showed that a global solution existed and that the domain was generally well behaved, in a numerical optimization perspective. The analysis of the information in these images was key in developing the optimization method to solve the inverse problem detailed in this work. With this inverse method, it was possible to obtain accurate volume scattering parameters for the samples under test. Analysing the objective function calculated with intensity distributions, instead of radiance distributions, is sufficient to ascertain the performance an intensity based inverse method will have. The intensity based objective function cross-sections for the volumetric domain and two samples, as used for the radiance analysis, are shown in Fig. 8.

 

Fig. 8 Intensity based objective function cross-sections perpendicular to the (a,d) μ_a axis, (b,e) the μ_s axis and (c,f) the g axis. The red cross indicates the true solution and correspond to a sample with a (a–c) short optical path length and (d–f) a long optical path length.

Download Full Size | PDF

In the cross-sections perpendicular to the μ_s and g axis, the global minimum is easily identifiable, for both the long and short optical path length samples. These objective function landscapes all share similar non-convex structures, with very narrow low error regions around the global minimum and a sharp slope towards the valley. For the cross-section perpendicular to the μ_a axis, the non-convexity of the objective function is more clearly visible. In this case the solution also resides in a very narrow valley but the slope towards the solution inside this valley is very faint. Moreover, for the low-optical path length sample, the μ_a cross-section does not contain an easily identifiable global minimum, despite the cross-section’s relatively high resolution. Taken together, these issues result in a challenging numerical optimization problem, especially when compared to the objective function cross-sections calculated using the radiance distribution for the same sample space.

According to these results, the objective function domains obtained using radiance distributions are significantly more descriptive than the ones obtained using intensity distributions. In the radiance based maps, the global minimum is always clearly observable and there is always a gradient towards the true solution. This does not, however, necessarily imply that intensity information can not be used to obtain an accurate solution for these sample domains. It only implies that much greater care is required in choosing an optimal numerical optimization method to successfully solve the problem.

3.3. Radiance versus intensity for fitting with noise

The impact of experimental noise on the results of the simulation study performed is tested using the noise model previously described. This noise model has a much larger impact on signals with low amplitude. Hence, only the objective function maps calculated for long optical path length samples are significantly affected and shown here. Applying noise to the simulated intensity and radiance distributions for the long optical path length sample results in the cross-sections shown in Fig. 9.

 

Fig. 9 Objective function cross-sections perpendicular to the (a,d) μ_a axis, (b,e) the μ_s axis and (c,f) the g axis. The cross-sections were calculated using (a–c) radiance distributions and (d–f) intensity distributions. The red cross indicates the real solution and correspond to a sample with μ_s = 33.69 mm⁻¹, μ_a = 7.30 mm⁻¹ and g = 0.89.

Download Full Size | PDF

In the radiance based objective function cross-sections (shown at the top of Fig. 9), the global minimum is still clearly observable along with a gradient towards it. There is barely a noticeable difference between the noiseless and the noisy case, only a very small broadening of the low error region can be observed. This demonstrates that the results of an inverse method that uses radiance are quite robust to experimental noise. For the intensity based objective function cross-sections (shown at the bottom of Fig. 9), the noise has a much larger impact, which is clearly shown by the significant broadening of the low error region around the solution. Additionally, this low error region for the noisy case seems to slightly shift away from the true solution. This indicates that volume scattering parameters estimated using intensity data with experimental noise are less likely to correspond to the real solution. This shows that for similar levels of noise in the intensity and radiance distributions, at least for long optical path length samples, a radiance based fitting method is more robust to noise than an intensity based fitting approach.

4. Conclusion

In this work, a method that uses the radiance distribution of light scattered from samples to estimate their volume scattering parameters was presented and its accuracy analysed. For the domain tested the method accurately estimated the volume scattering properties.

This work also tackled the more fundamental question regarding the difference in descriptive power between the radiance and the intensity distributions to obtain the volume scattering parameters of different samples. For the domains analysed, the intensity information contains less descriptive information. This results in a harder, but still solvable, optimization problem when designing an intensity based inverse method. However, when experimental error is considered, this no longer holds and using the radiance distribution becomes necessary to accurately estimate the volume scattering properties of the samples considered.

In essence, this work shows that intensity based inverse methods can be successfully applied to estimate the volume scattering properties of samples with short optical path lengths. For longer optical path lengths, a radiance based inverse method, such as the inverse Monte Carlo method presented, is much better suited.