1. Introduction

Volume scattering, or bulk scattering, is the angular and spatial redistribution of light due to scattering particles inside a material. It is a crucial effect that needs to be taken into account when analysing the propagation of light in many materials in a diverse set of fields such as biomedicine [1, 2], astronomy [3] and lighting [4]. For the specific case of lighting, the volume scattering properties of materials are tailored by manufacturers and researchers in order to develop diffusers with specific transmitted and reflected scattered light distributions. Such diffusers, are ubiquitous in lighting luminaires and display backlights [4–6].

Ray tracing is the fundamental method used in lighting to design and optimize optical configurations [7–9]. Volume scattering is included in most ray-tracing tools by using a statistical model. This model uses several parameters to describe the probability of scattering and absorption interactions inside a certain material and the angular deviation at each interaction. The correct parameters for this model can be directly obtained if enough information is available about the microscopic properties of the scattering particle added to the material, such as the refractive index, size, shape and number density [10]. Unfortunately, this information is rarely known, so, an inverse or fitting procedure is typically needed to estimate these parameters.

Such procedures fit simulation results to experimental data, which allows deriving the volume scattering parameters that reproduce the measured information. Inverse adding-doubling (IAD) [11] is an example of a popular fitting method to characterize the volume scattering properties of samples. It is based on adding-doubling (AD) [12, 13], so it shares most of its advantages and drawbacks. Adding-doubling uses the volume scattering properties of a sample to analytically calculate the resulting total integrated reflection and transmission under diffuse or collimated illumination. Because AD uses an analytical formulation, it is much faster than stochastic ray-tracing methods. This becomes a significant advantage when using it to solve the inverse problem of calculating the scattering parameters from experimental data, which typically involves brute-force iterative methods. However, AD is limited to very simple configurations. Only domains with planar geometries can be simulated in AD and the domain’s properties are assumed to be homogeneous. In addition, the angular scattering probabilities inside the domain are limited to rotationally symmetric distributions. In most situations, these limitations are not relevant, as characterization samples can be manufactured to have planar geometries and homogeneous scattering properties.

Recently, an extended IAD method was proposed, that fits to the reflected and transmitted intensity distributions in addition to fitting to the total integrated transmission and reflection [14]. Because of the limitations of the AD method, the calculated intensity distributions are all rotationally symmetric. Despite this, it was shown that the extended IAD method estimates volume scattering parameters which predict the scattered intensity distributions of experimental samples with better accuracy.

For lighting luminaires and display backlights, luminance is the fundamental quantity to investigate and optimize [15, 16]. In addition to the angular information of the scattered light distribution, the luminance also encompasses the spatial information. This means that the luminance characteristics of the diffuser used not only determine the radiation pattern of the total system, but also determine how the system itself is perceived. Unfortunately, the adding-doubling method only allows the calculation of the angular light distribution of volume scattering materials. Thus, the IAD method cannot be extended to fit to the luminance data and no sufficiently fast alternative inverse methods exist that can fit to this type of data. Because the extended IAD method is limited to intensity information, it is not clear if the obtained volume scattering parameters with this method are truly complete. To be considered complete for lighting and display purposes, the obtained volume scattering parameters must not only reproduce the intensity distribution accurately, but also the luminance distribution.

This vital question is tackled in this paper by performing a simulation study that only considers virtual volume scattering samples with pre-established properties. Each sample’s volume scattering parameters are used to simulate intensity and luminance distributions, which are considered as experimentally measured quantities throughout this study. Then, the intensity information is used to obtain a set of fitted volume scattering parameters with the extended IAD method. These parameters are tested by simulating the resulting luminance distributions with a ray-tracer, which are compared to the luminance distribution of the original, virtual sample. Using a comprehensive sampling of the volume scattering properties to create the virtual samples and a general phase function model allows covering a broad range of scattering scenarios (e.g. dependent and independent scattering).

To assess the impact of measurement errors and bridge the gap between simulations and measurement, an additional study is performed. The simulation study is repeated, but artificial noise is applied to the original samples’ intensity distributions before fitting with IAD. This emulates the case of having noisy experimental measurements and is used to investigate the impact of noise in the volume scattering parameters obtained.

Investigating the above questions with a simulation study, greatly reduces its complexity and avoids being held back by experimental problems, while still allowing to draw definite conclusions, as the paradigm is the same regardless of considering experimental measurements or virtual measurements.

2. Methods

2.1. Samples

To model volume scattering inside a sample with a ray-tracer, it is necessary to specify the scattering and absorption coefficient, µ_s and µ_a, and the single event phase function p(θ) parameters. The µ_s and µ_a parameters describe the probability per path length of a scattering or absorption event, respectively, while the single event phase function p(θ), parametrized by certain variables, describes the probability distribution of the angular deflection at each scattering event. For this reason, phase function models are usually employed. In this study, a generalization of a widely used phase function model, the Henyey-Greenstein (HG) phase function [5,11], termed multiple term Henyey-Greenstein (MTHG) was used. Starting from an individual HG phase function, as defined in Eq. (1), a MTHG is expressed as a weighted sum of individual Henyey Greenstein functions as shown in Eq. (2).

Each asymmetry parameter g_i, describes the directionality bias of each individual Henyey Greenstein phase function (g = 1 for forwardly peaked distributions, g = −1 for backscattering distributions). The w_i parameter is simply the weight or contribution that each HG phase function has in the MTHG phase function. The MTHG phase function model was chosen to increase the variability of the virtual samples created, which improves the robustness of the study. If the sample’s true parameters can be recovered using the intensity distribution, that implies that, even despite the variability in the phase function, there is enough information encoded in the intensity distribution.

In this study, the thickness d of the sample and µ_s are the parameters which show the greatest influence on the results. The absorption coefficient was kept fixed at µ_a = 0.1 as was the refractive index with n = 1.5. Only 6 different samples are considered, each having a µ_s linearly sampled from [0.5, 3.0] mm⁻¹ and an individual MTHG phase function, with N_g = 3, randomly generated. For each of these samples, 6 different thickness values were considered, linearly sampled from [0.5, 3.0] mm. This results in a set of 36 virtual samples. To have a more intuitive way to characterize the scattering regime of the samples created, an approximation to the average number of events per sample was calculated as $\overline{\mathit{\text{event s}}}={\mu }_s.d$. With the ranges established for µ_s and d, the average number of events for each sample covers both low scattering interaction regimes $\left(\overline{\mathit{\text{event s}}}=0.25\right)$ and high scattering interaction regimes $\left(\overline{\mathit{\text{event s}}}=9\right)$. This approximation, while useful, only holds for very forward scattering samples (g ≈ 1).

2.2. Inverse adding-doubling

The inverse adding-doubling process is essentially an iterative error minimization technique. At each iteration, the adding-doubling method is used to simulate the reflected and transmitted scattered intensity distribution for a test set of volume scattering parameters and this is compared to ’experimental’ measurements. If the fit between the two is not within a specified tolerance, a new test set is chosen. This is repeated until a good fit between the experimental and simulated intensities is obtained.

Due to its nature, IAD can be implemented using any error minimization algorithm. The Nelder-Mead simplex method is traditionally used in IAD for this purpose [11]. However, it can result in low accuracy fits for this study because of the high dimensionality that results from the order of the MTHG phase function used. Additionally, for non-smooth solution spaces, which may occur when solving simultaneously for all scattering parameters, the algorithm’s performance is strongly dependent on the starting point chosen. To mitigate these problems, a hybrid approach that performs a coarse search followed by a refinement step was taken. A genetic algorithm (GA) was used for the coarse search and a sequential quadratic programming (SQP) method was applied to refine the best result from the coarse search. Both methods are comprehensively described elsewhere [17,18], nevertheless the rationale behind coupling them and details on its implementation is given below.

The GA family of methods is well suited to global and sparse searches, which makes it an ideal candidate for this problem, where the possibility of local minima is high. With the parametrization used, the method will converge towards the vicinity of several possible global minima, while never converging on the actual minimum. This way, the domain is globally sampled. The best results are used as a starting seed to the SQP method to refine them by searching around their immediate vicinity, increasing the likelihood of identifying the global minimum. The software was implemented using Matlab (The MathWorks Inc., MA, USA), taking advantage of the built-in software for the genetic algorithm and the SQP method. At each iteration, the hybrid method calculates the error by taking the root mean squared difference (RMSE), defined in Eq. (3), between the intensity distributions, both reflected and transmitted, of the original and the test samples, as is shown in Eq. (4).

In both equations, the superscript o denotes the original samples, while f denotes fitted samples. Similarly, the subscript R identifies reflected intensity distributions, while T denotes transmitted intensity distributions. This differs from the original application of the extended IAD method, because in that case the total integrated quantities were explicitly considered in addition to the intensity distributions. Since the total integrated quantities are derived from the intensity distributions, a perfect match to the intensity information guarantees a perfect match to the total integrated quantities. The reverse, however, is not true. For this reason, the total integrated quantities were not considered in this study.

In an experimental situation, for thin samples or samples with low absorption and scattering (i.e. low concentration of scattering particles) the regular transmission T_c can be measured. The regular transmission describes the amount of light that is transmitted without being scattered or absorbed. Using this measured quantity can greatly improve the IAD fitting process, as one of the crucial fitting variables, namely the optical path length, can be obtained directly from the regular transmission. Internally, IAD converts the scattering and absorption coefficients into the dimensionless optical path length τ = d (µ_a + µ_s) and albedo a = µ_a/(µ_a + µ_s). If an accurate value for T_c is available, along with the Fresnel reflection r and transmission t at the interfaces between the sample and air, it is possible to calculate τ using Eq. (5). This effectively couples µ_s and µ_a, removing one as a fitting variable for IAD, reducing the dimensionality of the problem and greatly increasing the smoothness of the solution space.

Because this study is based on virtual samples, the accurate value of the regular transmission is always available. However, in an experimental situation with thick or highly scattering samples that is not the case because the regular transmission becomes very small. The ability to measure the regular transmission depends on not only the sample, but also on the measuring equipment used. Hence, to increase the inclusiveness of this study, all virtual samples were also fitted with IAD without using T_c. This results in a total of three sets of parameters: original samples, fitted samples obtained using T_c and samples fitted without considering regular transmission.

2.3. Luminance simulation

To compute the luminance distribution of the virtual samples, LightTools (Synopsys, Inc., CA, USA) was used. In LightTools, each sample was modelled as a rectangular geometry to which a material with certain volume scattering parameters was assigned. For each sample, rays were traced from a collimated light source with a single wavelength, incident at an angle of 60° to the normal. Because the scattering parameters are fixed, this study does not consider spectral information explicitly. Using different wavelengths for the illuminating source or changing scattering parameters to mimic different wavelengths would not change the results of this study, due to its comparative nature. Luminance maps were calculated at different azimuth and altitude angles. In the simulations performed, the altitude angles were sampled in [0°,70°]_7° and the azimuth angles in [90°,270°]_45°. Although this study uses luminance maps instead of the entire luminance distribution, they are sampled to include the most significant part of the distribution. Hence, any effect seen using the collection of luminance maps, i.e. the luminance characteristics of a sample, would also be present in the luminance distribution.

Because ray tracing is a stochastic process, two different simulations of the same sample may result in two significantly different luminance distributions. While these deviations decrease as more rays are sampled, it is impossible to eliminate them without going to impractically long simulation times. For each simulation in this study, 1 × 10⁷ rays were traced, which is a good compromise between simulation speed and simulation accuracy. The random deviations caused by the limited sampling greatly hinder drawing conclusions from directly comparing the simulated luminance distributions of different samples. To bypass this problem, the random number generator’s seed was fixed for all simulations. This forces each simulation to use the same sequence of random numbers, hence, any observed deviations between the simulated maps of fitted and original samples can only be caused by differences in their volume scattering properties. However, if the differences between the fitted and original scattering parameters are not significant, regardless of their magnitude, the deviations in simulated luminance characteristics will also be insignificant.

To compare the luminance characteristics, the normalized root mean squared error (NRMSE) was chosen [19]. The NRMSE calculates an absolute value for the difference and then normalizes it to original sample’s value range. This metric is defined in Eq. (6) and it is applied to each individual simulated luminance map of each sample. A NRMSE of 0% implies no difference between the luminance map of original and fitted samples.

To have a single representative value for the difference in luminance characteristics between each original and fitted sample, the highest deviations are chosen by taking the maximum value of the NRMSE in all luminance maps. This results in a comparison which considers the worst case scenario; for the luminance characteristics of two samples to be considered significantly different, it is only necessary for one luminance map out of the entire luminance distribution to be significantly different.

2.4. Experimental error

To investigate the applicability of the simulation study described in the previous section in an experimental situation, an additional study was performed. Instead of using the simulated intensity distribution directly as an input to IAD, a random source of error was added to the intensity signal. This emulates the case of measuring a sample’s scattered intensity distribution with an experimental setup, such as a CDD mounted on a goniophotometer [7]. The added noise mimics shot noise, which can be modelled as a zero mean Poisson distribution with a variance that depends on the collected signal’s intensity. In practice, the Poisson distribution can be approximated by a Gaussian distribution [20], which results in Eq. (7) for the added error. No other sources of error, such as alignment errors, were considered.

In this equation, I^e is the signal with error added, I^o is the original signal, S is the standard deviation scaling factor and r_N is a random number sampled from a zero-mean, unity standard deviation normal distribution. By varying the standard deviation’s dependence on the intensity, using S ∈ {1%, 5%, 10%}, sources of error with higher distorting influence were tested, as can be seen in Fig. 1. This artificial error formulation was chosen to create a very distorting source of error, although in practice the error sources are much more subtle.

 

Fig. 1 Transmitted intensities simulated for a virtual sample and corresponding distorted transmitted intensity simulated using a (a) 1%, a (b) 5% and a (c) 10% standard deviation scaling factor S.

Download Full Size | PDF

The IAD algorithm is applied to fit the distorted intensity distributions. The resulting volume scattering parameters are then used in LightTools to calculate the different luminance maps which are again compared to the original luminance data.

3. Results

3.1. Inverse adding-doubling

By using the extended IAD algorithm, an excellent fit to the reflected and transmitted intensity distributions is always achieved, both when using the regular transmission value T_c and without. The mean fitting RMSE and standard deviation for all samples is 1.5 × 10⁻⁰⁶ ± 3.2 × 10⁻⁰⁶ for the T_c case and 2.8 × 10⁻⁰⁶ ± 5.4 × 10⁻⁰⁶ for the other case. The fact that the fit in both cases is excellent, means that the optimisation algorithms succeeds for all considered samples in finding a set of volume scattering parameters that perfectly reproduce the scattered intensity information of the original sample. However, this does not necessarily mean that the obtained volume scattering parameters are the same as the parameters of the original sample.

The absolute difference between the original and fitted volume scattering parameters are shown in Figs. 2 and 3. The differences between original and fitted µ_a and µ_s for fits using T_c are very small throughout the sample domain. When not including T_c in the fit, the differences in scattering parameters values are still very low in most of the domain, except for very thick and scattering samples, where they rise sharply.

 

Fig. 2 Absolute difference between original and fitted µ_s for the samples fitted (a) using T_c and (b) without using T_c.

Download Full Size | PDF

 

Fig. 3 Absolute difference between original and fitted µ_a for the samples fitted (a) with T_c and (b) without T_c.

Download Full Size | PDF

The order g and weights w of the MTHG phase function can have different values while resulting in similar angular probability distribution. To investigate the phase function difference, the RMSE between original and fitted phase functions curves is calculated for all samples and shown in Fig. 4 for the T_c and no T_c cases. The landscape of differences between the phase function of fitted and original samples is very similar to the landscape of differences of the other volume scattering parameters.

 

Fig. 4 Value of the RMSE between original and fitted phase functions for the fitting case (a) that included T_c and (b) without including T_c.

Download Full Size | PDF

The fitting results show that, except for very thick and scattering samples, a very good approximation to the original volume scattering parameters is always obtained even without using T_c.

3.2. Luminance simulation

In Fig. 5 some calculated luminance maps are shown to illustrate their variation with the altitude and the azimuth angles. The results from the luminance simulations performed is shown in Fig. 6. It is clear that the landscape of differences of the luminance characteristics of fitted and original samples is very similar for both T_c and no T_c cases, although the latter shows slightly larger deviations. Despite this, the no T_c case only shows significant deviations for the samples with the highest thickness values and scattering coefficients tested (µ_s = 3.0 and d ≥ 2.5). Both landscapes closely mimic the landscape of differences between original and fitted volume scattering parameters for the samples that showed the highest average number of scattering events. However, there are noticeable deviations from this in the no T_c case for samples thicker than 2 mm. While not significant, these deviations show that in the absence of regular transmission measurements, it is preferable to use a thinner sample for volume scattering characterization.

 

Fig. 5 Luminance maps collected for the original {µ_s = 2.0 mm⁻¹, d = 2.0 mm} sample. Each row shows the luminance map collected at altitude angles 70°, 28°, 7° for (a–c) an azimuth angle of 270° and for (d–f) an azimuth angle of 180°.

Download Full Size | PDF

 

Fig. 6 Differences computed with the NRMSE between the luminance characteristics of the 36 original and fitted virtual samples for the (a) T_c case and (b) no T_c case.

Download Full Size | PDF

To have a qualitative idea of the difference between the luminance distributions of the original and fitted samples Fig. 7 shows the luminance map with highest deviation for the sample with highest NRMSE: {µ_s = 3.0 mm⁻¹, d = 3.0 mm}. In Fig. 8, the luminance map with highest deviation is shown for a slightly thinner and less scattering sample: {µ_s = 2.5 mm⁻¹, d = 2.5 mm}.

 

Fig. 7 Luminance maps for {µ_s = 3.0 mm⁻¹, d = 3.0 mm}. The first row (a–c) shows a luminance map of the (a) original sample, (b) the fitted sample and (c) a vertical cross section comparison for the T_c case. The second row (e–f) shows the same information for the fits without using T_c.

Download Full Size | PDF

 

Fig. 8 Luminance maps for {µ_s = 2.5 mm⁻¹, d = 2.5 mm}. Each row shows a luminance map of the (a) original sample, (b) the fitted sample and (c) a vertical cross section comparison, with (a–c) for the fitting with T_c and (e–f) for the fitting without T_c.

Download Full Size | PDF

There are two important conclusions that can be drawn from these results. First, the fit to the intensity, even without considering the regular transmission, returns volume scattering parameters that reproduce the luminance very well for all samples that are not very thick or highly scattering (i.e. $\overline{events}\le 6.25$). This indicates that in this sample range, the value of T_c is not necessary to obtain the complete volume scattering characterization of a sample with the IAD method for lighting or display applications. Secondly, fitting to the intensity without using T_c for highly scattering and thick samples can result in a set of parameters that, while perfectly reproducing the intensity distribution, fails to produce a good approximation to the original luminance characteristics. This implies that for these samples, fitting to the intensity information alone is not sufficient to have the correct volume scattering characterization that reproduces the luminance characteristics of the sample under investigation. While using the T_c value in this simulation study does improve the results, using it is not an option in an experimental situation, since it is typically very difficult to measure T_C with sufficient accuracy for these samples.

This essentially shows that the range of applicability of the extended IAD method to determine the volume scattering parameters that reproduce the luminance characteristics of a diffuser is limited to samples with $\overline{events}\le 6.25$. For the ranges of parameters used in this study, this corresponds to samples outside the range defined by µ_s ≥ 3.0 mm⁻¹ and d ≥ 2.5 mm. This limitation is not necessarily a significant caveat, as samples can be manufactured to conform to that range. If the concentration is kept fixed, the thickness of the sample does not influence its volume scattering properties. Thus, the scattering parameters extracted using a sufficiently thin sample can be directly used to simulate thicker samples. As shown in this study, even without using the regular transmission measurement, the characterization obtained with thin samples reproduces the luminance characteristics better than if using thick samples.

3.3. Experimental error

Adding noise to the simulated intensities results in fits with lower accuracy, as expected. The luminance differences calculated with the NRMSE are shown in Fig. 9 for all samples. The noise added to the intensity distribution has a marginal impact in the characterization if the accurate T_c value is still measurable. When it is not, large deviations in both µ_a and µ_s may lead to a poor characterization, which in turn leads to higher deviations in the simulated luminance characteristics. The noise added has a more pronounced effect in low scattering and thicker samples, changing the landscape of differences shown in the previous section for both cases. Only for the most distorted intensity distribution are there significant errors outside the viability range previously established. To accommodate strong noise sources, the viability range should be compressed to exclude samples with d > 2.0 mm. This way, even using the most distorted intensity distribution, a viability region for applying the extended IAD method to determine the volume scattering parameters that reproduce the luminance characteristics can still be established.

 

Fig. 9 Luminance deviations calculated with the NRMSE for (a–c) fits with T_c and (d–f) fits without T_c. Each column shows (a,d) the 1% error case, (b,e) the 5% error case and (c,f) the 10% error case.

Download Full Size | PDF

These results indicate that the extended inverse adding-doubling method can be applied, with confidence, to experimentally measured intensity values. Even with a significant error margin, the characterization obtained is robust and can be used to simulate the luminance characteristics of the sample.

4. Conclusion

By using a simulation study, it is shown that an extension of the IAD method can be applied to scattered intensity information to efficiently obtain a volume scattering characterization that reproduces the luminance characteristics for an extensive range of samples. It is also shown that the regular transmission, an important quantity in both IAD and extended IAD methods, is only necessary to obtain a good characterization for very thick and highly scattering samples, when it can no longer be accurately measured. Therefore, its influence does not increase the range of samples to which the extended IAD method can be applied.

Moreover, it is shown that adding noise to the scattered intensities prior to applying the extended IAD method has only a marginal deteriorating effect. It also reinforces the idea that thin samples should be preferred for characterization purposes, as they can be accurately fitted even in the presence of highly distorting noise. This increases the confidence of applying the investigated method in an experimental context.

These results will allow for a faster and easier characterization of diffuser samples, as they show that, for lighting and display purposes, it is only necessary to measure the intensity information from a relatively thin sample in order to obtain volume scattering parameters that accurately reproduce the luminance characteristics.