1. Introduction

Translucent materials, i.e. weakly scattering and absorbing materials, where light can enter deeply and be scattered progressively along its path, are common among manufactured products (polymer inks in 3D printing [1], dental repair biomaterials [2,3], ceramics [4],…) as well as in human tissues [5], raising issues in the control and reproduction of their appearance. Indeed, color and general appearance of translucent materials significantly depend on their thickness, in a way that a simplified light scattering model such as the Kubelka-Munk (or 2-flux) model [6], or even the 4-flux model [7–9] cannot predict accurately. In the Kubelka-Munk model, the light is assumed to propagate along two opposite directions, upwards and downwards, inside a slab. The incident light is assumed lambertian, and all the reflected or transmitted flux is collected. This model assumes that light is either absorbed or scattered with a lambertian angular distribution. The Fresnel reflections at the interfaces of the slab are accounted for by the Saunderson correction [10], which corrects the reflectance and transmittance calculated by the Kubelka-Munk model, taking into account the observation and illumination geometries, and refractive index mismatch. The validity of this correction was discussed by several authors for different applications [11,12]. One reason for the lack of accuracy of 2-flux models with translucent layers is the poor estimation of the internal reflectance of the interfaces bordering the material layer. Indeed, the Saunderson correction [10,13] relies on the assumption that the light reaching the interfaces has a lambertian type angular distribution, which may not be the case in thin and weakly scattering layers. More importantly, the exact angular distribution depends on the scattering and absorption properties of each layer, as well as its thickness. Following guidelines from the CIE [14], the reflectance or transmittance measurement geometries 1 the d:8$^{\circ }$ geometry (where a sample is illuminated with diffuse light in an integrating sphere and the emerging light is captured at an 8$^{\circ }$ angle) or likewise the 8$^{\circ }$:d geometry (where a sample is illuminated with directional light at an 8$^{\circ }$ angle and the emerging light is collected by an integrating sphere). Thus the incident light (or the captured light) is directional, whereas the 2-flux model assumes measurements based on a bi-hemispherical geometry. This has almost no impact in strongly scattering layers but may raise difficulties in translucent materials. As an example, it has been shown recently that in dental repair biomaterials, using an internal reflectance value of 4% (instead of the expected lambertian value of 60% for a refractive index 1.5) improves a lot the color prediction accuracy for slabs of various thicknesses [15]. The 4-flux model, proposed by Beasley et al. [7], popularized by Maheu et al. [8] and generalized by Vargas [9], considers the propagation of upwards and downwards fluxes inside a slab, either absorbed or scattered. Compared to the 2-flux model, the 4-flux also considers the propagation of directional light and the fraction of scattered directional light which contributes to the diffuse flux. Therefore, it is in principle more suited to reproduce experiments based on measuring devices with the 8°:d or d:8°geometry as recommended by the CIE [14] than the 2-flux model, which assumes a d:d geometry. Let’s stress that such a d:d geometry is very difficult to set up experimentally, although it might be equivalent to a d:8$^{\circ }$ or 8$^{\circ }$:d geometry when measuring highly scattering media. In the 4-flux model, the angular distribution of diffused light is also considered with the forward scattering ratio parameter (denoted $\zeta$, it indicates which proportion of the diffused flux is scattered in the forward direction) and the average path length parameter (denoted $\epsilon$, it reflects the average increase of distance traveled in the slab by scattered light rays compared to directional light rays). For isotropic scattering, the forward scattering ratio is 0.5 (same amount of light is scattered in the forward and the backward hemispheres) and the average path length parameter is 2. However, this simplifying assumption does not hold because single scattering is often anisotropic forwards, as in the case of Mie scattering by spherical particles. In the formulation by Maheu et al.[8], a symmetry condition was introduced based on phenomenological arguments. Indeed, it was assumed that the forward-scattering ratio $\zeta$ under the diffuse anisotropic radiation is equal to the forward-scattering ratio under collimated radiation, and that the average path-length parameter $\epsilon$ is the same for forward and backward diffuse radiations. The formulation proposed by Vargas [9], directly calculated from the radiative transfer equation, is more general as it does not rely on this approximation.

In principle, these issues could be solved if a more sophisticated solution of the Radiative Transfer Equation [16] (either by direct solving or by Monte Carlo methods) is used instead of 2- or 4-flux methods. However, these approaches are more time consuming, making challenging the inverse problem of extraction of optical parameters from a large set of experiments.

Following preliminary results presented by our group in [17], the aim of the study is to investigate this point in details using the formalism of the radiative transfer equation, and to compute more accurate values of the internal reflectance at the interfaces to be used in the simplified 2-flux model (such as the Saunderson correction for the Kubelka-Munk model [10,13]), in function of thickness, refractive index, absorption and scattering coefficients.

The method used to solve the radiative transfer equation in the case of homogeneous and flat scattering layers is presented first, then used to simulate the angular radiance distribution at their interfaces. The impact of absorption coefficient, layer thickness, scattering anisotropy in both translucent and strongly scattering materials is then carefully investigated. Finally, the accuracy of conventional and improved 2- and 4-flux models is discussed.

2. Solving the radiative transfer equation in a slab

The formalism and solution of the Radiative Transfer Equation in a slab of thickness $d$ has been extensively investigated by many authors, and details can be found in the following reference books [16], [18]. In this section, the approach used in this work is briefly presented for completeness. The radiative transfer equation (RTE) links the spatial evolution of the radiance $L$ (also called reduced intensity in [18]) at each point $\mathbf{{r}}$ and light direction $\mathbf{{s}}$ to the medium’s optical properties (absorption $\mu _a$ and scattering $\mu _s$ coefficients and phase function $p(\mathbf{{s}},\mathbf{{s'}})$):

Considering the scattering angle between $\mathbf{{s}}$ and $\mathbf{{s'}}$, the phase function $p(\mathbf{{s}},\mathbf{{s'}})$ describes the dependence of the scattered radiance on the scattering angle. In this work, we limit ourselves to the commonly used Henyey-Greenstein (HG) phase function [19], known to reproduce the forward peak of Mie scattering quite well, although in principle any type of phase function can be used. The HG phase function has the advantage of depending on one parameter only (the anisotropy parameter $g$). Several authors (see for instance [20,21]) have reported that it tends to underestimate backward diffuse radiation, explaining that modifications of the HG phase function have been proposed in the literature, at the cost however of an additional parameter. In common applications where it is intended to characterize a material of unknown composition by a limited set of optical measurements, it is necessary to find a trade-off between a good description of light scattering and a limited number of parameters to evaluate in order to prevent underdetermined inverse problems. The HG phase function seems to be a good solution in many practical cases. Of course, it is possible to extend the approach presented in our study to other phase function, in particular when the diffusing elements in the medium studied are well known, and/or when the optical measurements are better resolved angularly (e.g. BRDF / BTDF measurements). The following normalization procedure [18] was used where $\omega _0$ is the albedo:

In this work, the incident light is assumed to be a directional flux at normal incidence over the whole surface, giving rise to a symmetry of revolution around the normal direction, thus a constant radiance according to the azimuthal angle $\phi$. In consequence, in this 1D geometry, physical quantities only depend on the depth $z$ and the zenithal angle $\theta$.Introducing $\mu =\mathrm {cos}(\theta )$ and $\tau = (\mu _a+\mu _s)z$, the RTE can be re-written as:

The radiance $L$ can be split into a directional radiance contribution $L_N$ and a diffused radiance contribution $L_d$. The directional component is the contribution to the radiance of light that has not been scattered. Accounting for multiples reflections at the front side $\tau =0$ and back side $\tau =\tau _d=(\mu _a+\mu _s)d$ of the directional component, it is given by [22]:

Furthermore, we wish to emphasize that the assumptions made in this model regarding illumination and observation suggest that we can disregard lateral losses due to edge effects. Subsequently, we would like to discuss this specific point: are these assumptions valid from an experimental standpoint? In this study, we consider diffusing materials, transparent materials, and translucent materials. There is no edge loss effect for transparent materials as they are illuminated under normal incidence, with light propagating straightforwardly. For the others, namely diffusing and translucent materials, the light undergoes numerous grazing reflections. It can be guided far from its entry point, and the edge loss can consequently be significant, especially if the material is additionally weakly absorbing. However, this point has been discussed previously in the literature [26]. It is recommended to choose the illumination and observation zones in a way that minimizes measurement error caused by the edge loss effect. Indeed, as demonstrated in [26], it is appropriate to select the illuminated area (or respectively, the observed area) larger than the convolution between the material’s diffusion point spread function and the observed area (or respectively, the illuminated area). It then falls upon the experimenter to use a suitable measurement setup. Under these conditions, the assumptions of the model described here are applicable.

3. Semi-infinite layer: results and discussion

We first consider the case of a semi-infinite scattering medium (optical thickness $\tau _d=(\mu _a+\mu _s) d >> 1$). To simplify the discussion, the phase function is first assumed isotropic, i.e., $g = 0$. The impact of scattering anisotropy is discussed later on, in the case of a slab only. The main mechanisms occurring during light transport in the scattering material are illustrated on Fig. 2. First, the incident directional flux can be either reflected by the top interface or transmitted into the medium, where it is progressively scattered or absorbed. Some light propagating upwards may reach the interface, where it can be either transmitted or reflected. Note that all light rays with incident angle higher than the critical angle $i_c = \mathrm {arcsin}(1/n)$ are necessarily reflected (total internal reflection). It must be noted that in the semi-infinite material, the light that propagates upwards has necessarily undergone one or several scattering events. Let’s us now discuss the collective impact of these mechanisms on the overall angular distribution of the radiance at the top interface. In Fig. 3 are plotted, in polar coordinates, cut-views of the relative angular radiance $L_z(\theta )$ at $z = 0$ (named $L_0(\theta )$ hereinafter) as a function of the angle $\theta$ defined in Fig. 1, for a non-absorbing medium (blue curve) or absorbing medium (red curve).

 

Fig. 1. Studied configuration where a layer of scattering medium with a refractive index different from the surrounding air ("slab") is illuminated by directional light at normal incidence. Physical quantities depend on depth $z$ and angle $\theta$.

Download Full Size | PDF

 

Fig. 2. Photon trajectory inside a semi-infinite layer of scattering medium under directional incident light (irradiance $F_0$). The incident light can be either reflected at the interface (A) or transmitted into the medium. After a scattering event (B), a photon may propagate until reaching the interface. If the incident angle is below the critical angle ic, light can exit the medium or be reflected (C). If it is higher than $i_c$, light is totally reflected back to the medium (D). Percentages are given for $n=n_2/n_1=1.5$.

Download Full Size | PDF

 

Fig. 3. Cut view (a) and zoom of the cut view (b) of the relative angular radiance diagram near the top interface for semi-infinite scattering layers of two media of relative refractive index $n = 1.5$, with $g = 0$: a non-absorbing medium of scattering coefficient $\mu _s =5.0\,\mathrm {mm}^{-1}$ and albedo $\omega _0 = 1$ (blue curve), and an absorbing medium of extinction coefficient $(\mu _a+\mu _s) = 5.0\,\mathrm {mm}^{-1}$ and albedo $\omega _0 = 0.9$, resp $\omega _0 = 0.5$, $\omega _0 = 0.2$ and $\omega _0 = 0.1$ (red curve, resp. orange curve, yellow curve and brown curve). The light illuminating the interface from inside the layer propagates with zenithal angles $\lvert \theta \rvert > 90 ^{\circ }$ whereas the reflected light propagates with angles $\lvert \theta \rvert < 90 ^{\circ }$.

Download Full Size | PDF

In the case of a non-absorbing semi-infinite material (blue curve), the light that propagates upwards, after backscattering, is lambertian: the radiance is a constant $L_u$, the graph draws a hemi-circle. This is expected, and supports the common assumptions made in the literature. In the lower hemisphere, as set by the boundary conditions, the downward radiance $L_0(\theta )$ is the result of the internal reflection of the radiance propagating upwards, therefore the product of radiance $L_u$ with the Fresnel reflectance. For a refractive index $n$ of 1.5, this latter is equal to 0.04 at $\theta = 0^{\circ }$, and suddenly reaches 1 beyond the critical angle $i_c = \mathrm {arcsin}(1/n) \approx 42^{\circ }$. Let’s examine the corresponding diffused internal reflection coefficient $r_{id}$. The latter is defined as:

In the case of an absorbing semi-infinite medium (red curve), due to absorption, the overall radiance is of course lower than in the non-absorbing medium (both are plotted with the same scale in Fig. 3). Moreover, we observe that the upward radiance (in the upper hemisphere) is not constant anymore: it takes higher values at $\lvert \theta \rvert = 90^{\circ }$ than at $\theta = 180^{\circ }$. The light propagating upwards is therefore no longer lambertian. Let us comment this discrepancy. First, remind that at $z = 0$, the upward radiance results from the backscattering of directional light propagating downwards and diffused light propagating upwards. As previously seen, the downward radiance is not lambertian, because of the angular dependence of the Fresnel angular reflectance of the interface. However, in the non-absorbing medium, the upward radiance results from an infinite number of isotropic scattering events, which allows to get a lambertian distribution. On the contrary, in presence of absorption, the upward radiance results from a limited number of scattering events only. In consequence, although isotropic, scattering fails to distribute light uniformly over the sphere, which explains why the contribution to the upward flux of nearly normal radiances is lower than the one of grazing radiances. According to Fig. 3(b), the more the absorption increases, the less the upward luminance distribution is lambertian. Moreover, in Fig. S1 of Supplement 1, normalized radiation for different optical indexes has been computed at the same level of absorption, confirming that the non-lambertian shape of the radiance is indeed a consequence of combined effects of absorption and Fresnel reflections. The anisotropy of $L_0(\theta )$ modifies the calculation of the diffused internal reflectance $r_{id}$. In this case, the numerator of Eq. (8) changes less than the denominator which increases, resulting in an increase of $r_{id}$ to 0.62 (see Fig. 3). In the following, these conclusions are extended to the case of a finite thickness layer, with first isotropic, and then anisotropic scattering.

4. Finite thickness layer (slab): results and discussion

We now consider the case of a layer of finite thickness. Global internal reflection may occur at both interfaces in a similar way, although the situation of each interface is not symmetrical, as the directional light enters from the top interface only. In this case, not only the top (denoted $r_{id}$) but also the bottom (denoted $r_{id}'$) internal reflectance require attention. It is worth to be noted that $r_{id}$ and $r_{id}'$ refer to the internal reflectance of the diffused radiance component $L_d$ only. For this reason, we introduced two other internal reflectances (named global internal reflectances): $r_{ig}$ for the front interface, and $r_{ig}'$ for the back interface, that account for both directional and diffused components, such that

As the 4-flux model treats the propagation of directional and diffused light separately, $r_{id}$ and $r_{id}'$ should be used for the diffused component in the 4-flux model. On the contrary, as the 2-flux model ignores directional light, the reflectance $r_{ig}$ and $r_{ig'}$ should be used as correct parameters for this latter approach.

4.1 Radiance angular variations and internal reflectances in the case of an isotropic scattering phase function

Figure 4(A) (resp. Figure 4(B)) represents diffused radiances $L_0(\theta )$ and $L_d(\theta )$ at the top and bottom interfaces as functions of the zenithal angle in isotropic cases in absence of absorption (respectively in presence of absorption) for various optical thicknesses (the sum $\mu _a+\mu _s$ of the absorption and scattering coefficients remaining constant in our simulations).

 

Fig. 4. Polar plot of the diffused radiance at the top interface, $L_0(\theta )$ (Fig. 4Aa), resp. (Fig 4Ba), and the bottom interface, $L_d(\theta )$ (Fig. 4Ab), resp. (Fig 4Bb), for isotropic layers of different optical thicknesses in absence of absorption (refractive index $n = 1.5$, albedo $\omega _0 = \mu _s / (\mu _a+\mu _s)=1.0$), resp. in presence of absorption (refractive index $n = 1.5$, albedo $\omega _0 = \mu _s / (\mu _a+\mu _s)=0.7$).

Download Full Size | PDF

It’s natural to question the relevance of scattering within an extremely thin optical layer. If the layer has a negligible thickness, there’s no space left for diffusers. Moreover, the model becomes dubious when the size of the diffusers approaches the same order of magnitude as the thickness. Therefore, one might ask at what optical thicknesses the simple scattering approximation applies. While the threshold for the simple scattering approximation remains a topic of debate [27], an optical thickness of $\tau _d=0.35$ is definitively above this threshold. Consequently, we confine our analysis to optical thicknesses exceeding 0.35.

For large optical thicknesses, i.e., $(\mu _a + \mu _s) d >> 1$ (red-orange curves), with or without absorption, for an isotropic phase function, the diffused radiance distribution at the top interface tends to the one observed for a semi-infinite medium discussed in the previous section (Fig. 3): the slab simply behaves as a semi-infinite medium. In absence of absorption, the radiance distribution at the bottom interface also tends to be a constant, and thus, $r_{ig}$ (resp. $r_{id}$) tends to be equal to $r_{ig}'$ (resp. $r_{id}'$). In presence of absorption, the diffused radiance distribution at the bottom interface tends to be negligible for large optical thicknesses. Even if it is possible (although difficult) to compute a diffused and global internal reflectance at this bottom interface, as the bottom radiance is negligible, this quantity does not really have any physical meaning. Consequently, this case will not be discussed hereafter.

For small optical distance, i.e., $(\mu _a + \mu _s) d << 1$ (blue curves), however, there is less scattering, but also less absorption, than large optical thicknesses. We thus observe a similar radiance pattern at both interfaces, in presence or in absence of absorption: radiances $L_0(\theta )$ and $L_d(\theta )$ become weak for incidence angles lower than the critical angle $i_c$. Indeed, the radiance at higher incidence angles is mainly caused by the total internal reflection at each interface. This effect strongly depends on the optical index value, as shown in Supplement 1 (Fig. S11). In consequence, upward and downward radiances tend to become almost identical at each interface. For this reason, the numerator and denominator of Eq. (8) tend to be equal for the diffused internal reflectance at both interfaces, resulting in in an increase in the diffuse internal reflectance value as $d$ decreases.

The diffused internal reflectance of the top and bottom interfaces, respectively $r_{id}$ and $r_{id}'$, are plotted in Fig. 5(Aa), resp. Figure 4(Ab), as functions of the optical thickness for various albedo values. We retrieve the fact that when the optical thickness $\tau _d$, is very small, the diffused internal reflectance value increases at both interfaces. As the optical thickness increases, the diffused internal reflectance decreases. It asymptotically approaches 0.60, the value usually considered in the Saunderson correction, for the top interface when the medium is non-absorbing ($\omega _0 = 1$), or a higher value than 0.60 if the medium is absorbing (0.63 when $\omega _0 = 1$).

 

Fig. 5. Internal reflectances versus optical thickness $(\mu _a+\mu _s)d$ for isotropic scattering ($g=0$) and several albedo $\omega _0 = \mu _s /(\mu _a+\mu _s)$ with $\mu _a+\mu _s = 5\,\mathrm {cm}^{-1}$. The bold dashed line at $r_i$ = 0.6 corresponds to the value usually considered in the Saunderson correction.

Download Full Size | PDF

Let’s now examine the case of the global internal reflection $r_{ig}$ and $r_{ig}'$ (Fig. 5(Ba) and Fig. 5(Bb)). First of all, it should be noted that the two coefficients may differ, in particular at lower optical thickness: when $r_{id}$ and $r_{id}'$ s.t.tend to 1increase, their counterparts $r_{ig}$ and $r_{ig}'$ take much lower values. Indeed, at low optical thickness, the contribution of diffused light to the total radiance decreases, explaining the decrease of the global internal reflectance, which tends to 4%; the Fresnel coefficient at normal incidence. This effect is particularly strong for low albedo values, as absorption prevents scattering to randomize the light direction. At high optical thickness however, both $r_{ig}$ and $r_{ig}'$ tend to $r_{id}$ and $r_{id}'$ as expected, justifying the use of the lambertian value 0.6 in the 2-flux model.

4.2 Radiance angular variations and internal reflectances in the case of an anisotropic scattering phase function

It is known that light scattering mechanisms are usually anisotropic, except in the particular case of scattering centers smaller than the wavelength (Rayleigh scattering). Living tissues, such as skin, are a well-known examples of turbid media where light is mostly scattered in the forward direction [28]. For this reason, we restrict ourselves to an anisotropic scattering phase function with positives values of the anisotropy parameter $g$. However, note that results with negative value of the anisotropy parameter $g$ are given in the Supplement 1 (Figs. S6-S9).

Figure 6(A) (resp. Figure 6(B)), represents the diffused radiances $L_0(\theta )$ and $L_d(\theta )$ at the top and bottom interfaces as functions of the zenithal angle in the anisotropic case ($g = 0.8$) in absence of absorption (respectively in presence of absorption), for various optical thicknesses of the medium. With or without absorption, for large optical thickness, i.e., $(\mu _a + \mu _s) d >> 1$ (red-orange curves), the diffused radiance distribution at the top interface tends to the one observed for a semi-infinite medium as discussed in the previous section (Fig. 3): the high number of scattering events occurring within the slab tends to compensate the effect of anisotropy in this case.

 

Fig. 6. Polar plot of the diffused radiance at the top $L_0(\theta )$ (Fig. 6Aa), resp (Fig 6Ba), and bottom $L_d(\theta )$ (Fig. 6Ab), resp. (Fig 6Bb); for layers of different optical thicknesses in presence of absorption (refractive index $n = 1.5$, albedo $\omega _0 = \mu _s / (\mu _a+\mu _s)=1.0$, anisotropy factor $g=0.8$), resp in presence of absorption (refractive index $n = 1.5$, albedo $\omega _0 = \mu _s / (\mu _a+\mu _s)=0.7$, anisotropy factor $g=0.8$).

Download Full Size | PDF

When $(\mu _a + \mu _s) d << 1$ (blue curves), with or without absorption, the impact of anisotropy becomes significant, in particular for the radiance at the bottom interface. As the value of the anisotropy factor $g$ is positive, $\theta = 0^{\circ }$ appears to be the privileged direction for scattering light. For this reason, a peak in the $\theta = 0^{\circ }$ direction clearly appears in the bottom diffused radiance. Moreover, the radiance overall values are weaker than in the isotropic case, as the scattered light, preferentially directed in the $\theta = 0^{\circ }$ direction, can emerge more easily from the back without being affected by any internal reflection. At the front interface, a small peak appears at $\theta = 180^{\circ }$, corresponding to the light reflected back by Fresnel reflection at the bottom interface.

Let’s consider the impact of scattering anisotropy on the diffused- and global internal reflectances versus optical thickness (Fig. 7). Note that the isotropic case ($g=0$) is also plotted (solid line) for comparison, together with a medium ($g=0.5$) and strong ($g=0.8$) anisotropy. As previously, we will not discuss hereafter the particular case of the bottom interface in presence of absorption, the radiance being negligible at the bottom interface in this case. For all the other cases, anisotropy does not seem to significantly impact the trends previously discussed for both diffused and global internal reflectances.

 

Fig. 7. Internal reflectances versus optical thickness $(\mu _a+\mu _s)d$ for several albedo $\omega _0 = \mu _s /(\mu _a+\mu _s)$ with $\mu _a+\mu _s = 5\,\mathrm {cm}^{-1}$ and for anisotropic scattering: colored solid lines correspond to $g=0$, colored dashed lines correspond to $g=0.5$ and colored dotted lines correspond to $g=0.8$. The bold dashed line at $r_i = 0.6$ corresponds to the value usually considered in the Saunderson correction.

Download Full Size | PDF

The most significant differences when comparing isotropic with anisotropic situations are a modest increase of the top internal reflectance $r_{id}$ which can be observed at high optical thickness (0.7 instead of 0.6 for $\omega _0=0.1$), and, due to the light scattered in the normal direction, a faster decrease of the top global internal reflectance $r_{ig}$ as the optical thickness tends to 0 in anisotropic situations.

In conclusion, despite significant changes in term of radiance directivity, scattering anisotropy does not introduce significant changes in term of internal reflectances, especially for large optical thickness (refer to solid lines of Figs. S2-S5 of Supplement 1).

5. Validity of the 2- and 4-flux models with and without internal reflectances calculated by the radiative transfer equation

Previous simulations have shown that internal reflectances may differ significantly from the lambertian value in translucent materials. In this section, we would like to assess its impact on the prediction accuracy of conventional optical models, and see if it is possible to improve it by simply replacing the lambertian one by the ones calculated by the radiative transfer equation.

For the purpose of this experiment, the propagation of light inside translucent slabs of variable thickness was simulated with the RTE, and the reflectance and transmittance factors of each slab as well as the internal reflectances $r_{id}, r_{id}', r_{ig}, r_{ig}'$ at their interfaces were calculated. The slabs are characterized by a scattering coefficient $\mu _s = 4.5\,\mathrm {mm}^{-1}$ and an absorption coefficient $\mu _a = 0.5\,\mathrm {mm}^{-1}$. These values are consistent with measurements performed on translucent skin tissues at 500 nm or 550 nm [5,28–30]. Scattering was modelled using the Henyey-Greenstein phase function with $g = 0$ (isotropic scattering), $g = 0.5$ and $g = 0.8$. Slabs with optical thicknesses ranging from 5.10^-4 to 7 (thickness ranging from 10^-4 mm to 1.4 mm) were simulated for each value of $g$.

Then, the absorption $k$ and scattering $s$ coefficients of the best available model, i.e. the 4-flux model, were extracted using a fitting algorithm to reproduce at best the RTE data (taken as reference) for all thicknesses. Note that here the $r_{id}$ and $r_{id}'$ values calculated with the RTE were used to extract the optical parameters of the 4-flux model. Here, we limit ourselves for simplicity to the approximated 4-flux approaches proposed by Rozé [31]. Consequently, the forward scattering ratio $\zeta$ is calculated according to formula (11) as proposed in [31] :

The average path length parameter $\epsilon$ was set to 2, which corresponds to semi-isotropic scattering, as usually done in 2-flux and 4-flux models. We assume that it is constant with respect to the layer thickness, and that its value is the same at both interfaces of the slab. This is a strong assumption which will necessarily limit the accuracy of all models, but it simplifies the study and allows to draw general conclusions on the models’ accuracy. Note that the RTE enables to calculate $\epsilon$ at any depth within the slab thanks to its integral definition given in [31], and use different values for the top and bottom interfaces, as highlighted in the 4-flux model from Vargas et al. [32,33]. Explicit expressions for both forward and backward average path-length parameter have been derived in [34], which is theoretically valid in the limit of dilute suspension of spherical particles.

The fitting algorithm is based on the GlobalSearch global optimization solver from Matlab, and searches for the parameters that minimize the average root mean square difference between the reflectance and transmittance factors predicted and simulated with the RTE (ground truth).

Finally, the absorption and scattering coefficients of the 2-flux model are deduced from those of the 4-flux model using the well-known equations(12,13) [35,36]:

For the 2-flux model, several possible values of global internal reflectance $r_{ig}$ and $r_{ig}'$ were considered:

For the 4-flux model, only two cases where $r_{id} = r_{id}' = 0.60$ and the $r_{id}$ and $r_{id}'$ values calculated by the RTE were considered. The $r_{id}, r_{id}',r_{ig}$ and $r_{ig}'$ values are given in the Fig. S10 of the Supplement 1. The reflectance and transmittance factor predictions of the 2- and 4-flux models are reported on Fig. 8. The 4-flux model with the $r_{id}$ and $r_{id}'$ values calculated by the RTE is found the most accurate for predicting the reflectance and the transmittance, regardless of the g value. Of course, it is the consequence of the construction of the experiment ($k$ and $s$ are extracted from RTE data using this model), but such agreement would not have been possible if the 4-flux model with the $r_{id}$ and $r_{id}'$ values calculated by the RTE was not able to reproduce accurately RTE results.

 

Fig. 8. Reflectance (left graphs) and transmittance (right graphs) factors predicted by the 2-flux and 4-flux models for different values of the internal reflectance at the interface. The red dotted curves represent the 2-flux model with $r_{ig} = r_{ig}' = 0$ (ie. no Saunderson correction). The magenta dotted curves represent the 2-flux model with $r_{ig} = r_{ig}' = 0.04$. The yellow (resp. blue) dotted curves represent the 2-flux (resp. 4-flux) model with $r_{ig} = r_{ig}' = 0.60$ (ie. lambertian flux). The green (resp. black) dotted curves represent the 2-flux (resp. 4-flux) model with $r_{ig}$ (resp. $r_{id}$) and $r_{ig}'$ (resp. $r_{id}'$) calculated by the RTE for each slab. The black solid lines represent the ground truth values simulated with the RTE. a) Simulations of the RTE are performed for $g = 0$. b) Simulations of the RTE are performed for $g = 0.5$. c) Simulations of the RTE are performed with $g = 0.8$.

Download Full Size | PDF

In reflectance, the 4-flux model with $r_{id}$ and $r_{id}'= 0.60$ is the most accurate for $g=0$, and for slabs thicker than $\tau _d = 3.5$ for $g=0.5$. For $g = 0,\,\epsilon = 2$ and $\zeta = 0.5$, the simulation corresponds to the assumptions of the 4-flux model, i.e., lambertian flux within the layer and at the interfaces, which explains why this model fits well with the ground truth. The only error comes from the use of different values of $r_{id}$ and $r_{id}'$. The 2-flux models are less accurate than the 4-flux models, and do not even converge towards the 4-flux model for high optical thickness which is inherent to formulae (12) and (13). For $g = 0.8$, the 2-flux model with $r_{ig}$ and $r_{ig}'$ values calculated with the RTE is the most accurate model. This tends to show that adapting the value of $r_{ig}$ and $r_{ig}'$ is more important for highly asymmetric scattering. 2-flux and 4-flux models with $r_i = 0.60$ are notably inaccurate for thin slabs, but their accuracy increases for slabs of high optical thickness. The 2-flux model with $r_{ig} = 0.04$ is inaccurate with $g = 0$ and $g = 0.5$, but is fairly accurate for very thin slabs, i.e., $\tau _d \leq 1$ when $g = 0.8$. The 2-flux model with $r_{ig}$ and $r_{ig}' = 0$ is inaccurate in all the tested cases.

In transmittance, regardless of the value of $g$, the 2-flux model with $r_{ig}$ and $r_{ig}' = 0.60$ is inaccurate. For $g = 0$, the 4-flux model with $r_{id}$ and $r_{id}'= 0.60$ is the most accurate as in reflectance mode. The 2-flux model with $r_{ig}$ and $r_{ig}'$ values calculated with the RTE is more accurate than the other 2-flux models, with the 2-flux model with $r_{ig}$ and $r_{ig}' = 0.60$ being the least accurate. For $g = 0.5$, the 2-flux model with $r_{ig}$ and $r_{ig}' = 0.04$ is the most accurate, but the other models (except the 2-flux model with $r_{ig}$ and $r_{ig}' = 0.60$) remain fairly accurate. For $g = 0.8$, the 2-flux model with $r_{ig}$ and $r_{ig}'$ values calculated with the RTE is the most accurate, as in reflectance mode. The other models are notably less accurate, except for the 4-flux model with $r_{id}$ and $r_{id}' = 0.60$ for very thick slabs and the 2-flux model with $r_{ig}$ and $r_{ig}' = 0.04$ for very thin slabs.

It is worth noting that the reflectance and transmittance factors are perfectly predicted by the 4-flux models when the optical thickness tends towards 0 regardless of the value of $r_{id}$, and the 2-flux model with the $r_{ig}$ values calculated with the RTE or with $r_{ig} = 0.04$. Indeed, for an extremely thin slab (with optical thickness almost 0), scattering and absorption become negligible and light is only reflected by the top and bottom interfaces. Therefore, the slab is equivalent to a layer of clear glass (non-scattering and non-absorbing medium with refractive index $n = 1.5$), and its reflectance and transmittance factors can be calculated easily with Fresnel laws. This trend, which must be accounted for at the interfaces of the layer, is only satisfied by the aforementioned 2-flux and 4-flux models.

Within the limits of this study where the 4-flux model with $r_{id}$ and $r_{id}'$ calculated by RTE has been taken as a reference, we conclude that the 4-flux model with $r_{id} = r_{id}' = 0.60$ is the most accurate model for $g = 0$ and $0.5$ in reflectance and transmittance mode. However, for higher values of $g$, the model is not accurate enough and specific values of $r_i$ must be used. This is why the 2-flux model with $r_{ig}$ and $r_{ig}'$ values calculated with the RTE is the most accurate model for $g = 0.8$ in this particular case. The other 2-flux models are generally less accurate than the 4-flux model, although the 2-flux model with $r_{ig} = r_{ig}' = 0.04$ is more accurate in transmittance mode.

In other words, these preliminary results suggest that in all cases, for an accurate material parameter extraction using 2- or 4-flux models in translucent material, internal reflectance should be modified in agreement with radiative transport theory, especially for highly anisotropic scattering. If conventional values of internal reflectance with $r_i= r_i' = 0.60$ have to be used for simplicity, it is highly recommended to prefer the 4-flux model rather than the 2-flux model. Indeed, the 4-flux model naturally takes into account the decrease of the diffuse flux compared to the specular flux in translucent materials, whereas the 2-flux models need to account for it by setting appropriate value the global internal reflectances.

We also investigated the capacity of the four-flux model with or without $r_i$ and $r_i'$ values calculated with the RTE to extract optical parameters for one thickness and then to use them to predict the spectral reflectance and transmittance for other thicknesses. Spectral reflectances and transmittances were generated with the RTE for $n = 1.5$, $\mu _a$ varying from 0.5 to 0.01 mm^-1, $\mu _s$ varying from 3.4 to 2.2 mm^-1, $g$ either equal to 0 (isotropic scattering) or $0.75$, and for varying layer thickness, in order to simulate measurements performed on translucent samples with a spectrophotometer. The measurements of the samples with thickness 1.2 mm were used to extract the optical parameters according to the four-flux model with $r_i = r_i' = 0.60$ and the four-flux model with $r_i$ and $r_i'$ given by the RTE by means of an optimization algorithm. Then, these coefficients were used to predict the spectral reflectances and transmittances of the other samples. The prediction accuracy of each model is then assessed by calculating the Root Mean Square Difference between the predicted and the generated spectra, averaged on all wavelengths. The prediction accuracy for the reflectance and transmittance is displayed in Fig. 9(a) for the isotropic case ($g = 0$) and Fig. 9(b) for the anisotropic case ($g = 0.75$). The reflectance and transmittance predictions obtained with the four-flux model with $r_i = r_i' = 60{\% }$ show discrepancies, especially for thin layers. However, the predictions obtained with the four-flux model with $r_i$ and $r_i'$ given by the RTE are much more accurate, especially for the reflectance and for both $g=0$ and $g=0.75$. This further highlights the accuracy improvement enabled by the use of $r_i$ and $r_i'$ values rigorously calculated with the RTE instead of using the Lambertian approximation in the case of translucent material.

 

Fig. 9. Root Mean Square Difference (RMSD) assessing the deviation between the generated spectral reflectance (on the left) and transmittance factors (on the right) of samples with various thickness and the corresponding predictions given by the four-flux models with or without $r_i$ and $r_i'$ values calculated with the RTE. The sample with 1.2 mm is used for calibration of the models. a) The spectral reflectance and transmittance are generated with $g = 0$ (isotropic scattering). b) The spectral reflectance and transmittance are generated with $g = 0.75$.

Download Full Size | PDF

6. Conclusion

Detailed simulations based on the numerical solving of the radiative transfer equation (RTE) were used to calculate the internal diffused reflectance, a parameter needed in the 2-flux or 4-flux models used for optical parameter extractions and appearance prediction. In this work, both isotropic and anisotropic scattering were considered, even though we found that anisotropy does not significantly impact the overall conclusions.

For large optical thicknesses, the internal diffused reflectances at both interfaces were found in agreement with the expectation of the lambertian approximation when absorption is negligible. In presence of absorption however, these internal reflectances slightly differ, the top one being a bit higher, and the bottom one a bit lower than the lambertian value.

For intermediate and low optical thicknesses (translucent material), it turns out that both internal reflectances tends to be approximately equal, and exceed the value for a lambertian distribution. On the opposite, global internal reflectances (including both directional and diffused light), required for 2-flux models, tends to directional values at low optical distances.

These results are the consequence of the competition between the directional reflection (Fresnel formulae), which induces anisotropy in the angular distribution of radiance (low incidence angle rays may escape from the slab while high incidence angle rays remain trapped by global internal reflection), and scattering, which tends to distribute light uniformly.

Finally, the impact of the internal reflectances on extraction performed using the 2- and Maheu’s 4-flux models was investigated, using RTE simulations as a reference. It turns out that the commonly used approximation, i.e., keeping constant the internal reflection equal to its lambertian value, can lead to significant errors (up to 100%) for 2-flux models when dealing with translucent materials. Moreover, parameter extraction with the 2- and 4-flux models can be significantly improved using the internal coefficient calculated using the radiative transfer equation, especially in translucent material with highly anisotropic scattering medium.

These results open the way to a better prediction of the spectral reflectance and transmittance, thereby the color appearance, of slices of translucent materials with respect to their thickness.

Appendix: benchmark of the radiative transfer equation solver

In order to validate the implemented model, benchmarks with a widely used Monte Carlo code [23] with different sets of parameters (relative refractive indices, absorption and scattering coefficients, layer thickness, and anisotropy parameter g of the Henyey-Greenstein phase function) were performed with success, considering the same geometrical configurations and physical parameters (not shown here). Note that GPU time of this Monte Carlo code has been found 6 times longer than our RTE solver at the same level of accuracy.

In addition, we report here in Fig. 10 and Fig. 11 successful comparison with literature data for two specific cases with and without index-matching and considering isotropic or anisotropic phase function.

For an anisotropic phase function and a slab geometry of finite thickness with index matching, van de Hulst’s tables [24] served as a reference for total reflection and transmission. The Henyey-Greenstein phase function was considered with an anisotropy parameter $g$ of 0.75. The albedo was 0.9. According to Table 1, one can notice that our work is consistent with existing values of van de Hulst et al. [24]. Moreover, the angular dependencies of relative radiance in a slab with index matching and with or without absorption have also been benchmarked with van de Hulst et al. data, and are reported in the following figure.

 

Fig. 10. Comparison of mean luminance given by the Monte Carlo code [23] (blue dashed curve) and our RTE solver (orange solid curve) for (a) for semi-infinite scattering layer of relative refractive index $n = 1$, with $g = 0.6$, albedo $\omega _0 = 0.9$ and an extinction coefficient $(\mu _a+\mu _s) = 10.0\,\mathrm {mm}^{-1}$ (b) slab scattering layer of relative refractive index $n = 1.5$, with $g = 0$, albedo $\omega _0 = 0.9$, optical thickness $\tau _d=5$ and an $(\mu _a+\mu _s) = 5.0\,\mathrm {mm}^{-1}$.

Download Full Size | PDF

 

Fig. 11. Comparison of cut views of the relative angular radiance diagram near the top (a) and bottom (b) interface between our solver (solid lines) and values given by Van de Hulst (crosses) for a slab scattering layer of optical thickness $\tau _d=4.00$ and relative refractive index $n = 1.0$, with $g = 0$ at various albedo $\omega _0 = 0.95, \omega _0 = 0.9, \omega _0 = 0.6$ and $\omega _0 = 0.2$.

Download Full Size | PDF

Table 1. Total Reflection and Transmission comparison between van de Hulst [24] and our work for an anisotropic scattering slab medium.

View Table | View all tables in this article

Table 2. Total Reflection comparison between Giovanelli [25] and our work for a semi-infinite medium which is not index matched with surrounding medium ($n = 1.5$) and the albedo is 0.9.

View Table | View all tables in this article

Finding exact solutions for media which are not index matched is more difficult. Giovanelli’s work [25] provides data only for a semi-infinite medium with isotropic scattering, an optical index mismatch of 1.5 and an albedo of 0.9. According to Table 2, one can notice that our work is consistent with existing values given by Giovanelli. [25]

Funding

EUR Manutech Sleight (ANR-17-EURE-0026); France Life Imaging (ANR-11-INBS-0006); LabEx PRIMES (ANR-11-IDEX-0007, ANR-11-LABX-0063).

Acknowledgments

This work has been funded by LABEX PRIMES (ANR-11-LABX-0063) of Université de Lyon, within the program "Investissements d’Avenir" (ANR-11-IDEX-0007) operated by the French National Research Agency (ANR) and carried out within the framework of France Life Imaging (ANR-11-INBS-0006). It has also been funded by a public grant from the French National Research Agency (ANR) under the "France 2030" investment plan, which has the reference EUR MANUTECH SLEIGHT - ANR-17-EURE-0026.

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

Supplemental document

See Supplement 1 for supporting content.

References

1. T. P. Van Song, C. Andraud, and M. V. Ortiz-Segovia, “Towards spectral prediction of 2.5D prints for soft-proofing applications,” in 2016 Sixth International Conference on Image Processing Theory, Tools and Applications (IPTA), (IEEE, 2016), pp. 1–6.

2. L. M. Schabbach, B. C. dos Santos, L. S. De Bortoli, et al., “Application of Kubelka-Munk model on the optical characterization of translucent dental zirconia,” Mater. Chem. Phys. 258, 123994 (2021). [CrossRef]  

3. I.-S. Pop-Ciutrila, R. Ghinea, M. d. M. Perez Gomez, et al., “Dentine scattering, absorption, transmittance and light reflectivity in human incisors, canines and molars,” J. Dent. 43(9), 1116–1124 (2015). [CrossRef]  

4. S. Rio, C. Andraud, P. Deniard, et al., “Study of the influence of crystalline phases on optical characteristics of a glass-ceramic in the visible range via simulations by the four-flux method,” J. Non-Cryst. Solids 551, 120446 (2021). [CrossRef]  

5. M. Doi and S. Tominaga, “Spectral estimation of human skin color using the Kubelka-Munk theory,” in Color Imaging VIII: Processing, Hardcopy, and Applications, vol. 5008 (SPIE, 2003), pp. 221–228.

6. P. Kubelka, “Ein Beitrag zur Optik der Farbanstriche (Contribution to the optic of paint),” Zeitschrift fur technische Physik 12, 593–601 (1931).

7. J. Beasley, J. Atkins, and F. Billmeyer Jr, “Scattering and absorption of light in turbid media,” ICES II: Electromagnetic Scattering p. 765 (1967).

8. B. Maheu, J. Letoulouzan, and G. Gouesbet, “Four-flux models to solve the scattering transfer equation in terms of lorenz-mie parameters,” Appl. Opt. 23(19), 3353–3362 (1984). [CrossRef]  

9. W. E. Vargas, “Generalized four-flux radiative transfer model,” Appl. Opt. 37(13), 2615–2623 (1998). [CrossRef]  

10. J. L. Saunderson, “Calculation of the Color of Pigmented Plastics*,” J. Opt. Soc. Am. 32(12), 727–736 (1942). [CrossRef]  

11. S. E. Orchard, “The missing variable: internal surface reflection,” Color Res. Appl. 2(1), 26–31 (1977). [CrossRef]  

12. W. M. Johnston, W. J. O’Brien, and T.-Y. Tien, “The determination of optical absorption and scattering in translucent porcelain,” Color Res. Appl. 11(2), 125–130 (1986). [CrossRef]  

13. M. Hébert, R. D. Hersch, and P. Emmel, “Fundamentals of Optics and Radiometry for Color Reproduction,” in Handbook of Digital Imaging, (John Wiley & Sons, Ltd, 2014), pp. 1–57.

14. CIE, Absolute methods for reflection measurements (CIE Technical Report, 1979).

15. V. Duveiller, L. Gevaux, R. Clerc, et al., “Reflectance and transmittance of flowable dental resin composite predicted by the two-flux model: on the importance of analyzing the effective measurement geometry,” Color. Imaging Conf. 28(1), 313–320 (2020). [CrossRef]  

16. S. Chandrasekhar, Radiative Transfer, Dover Books on Intermediate and Advanced Mathematics (Dover Publications, 1960).

17. A. Gautheron, R. Clerc, V. Duveiller, et al., “Light scattering in translucent layers: angular distribution and internal reflections at flat interfaces,” Electronic Imaging 34(15), 221-1–221-6 (2022). [CrossRef]  

18. A. Ishimaru, Wave propagation and scattering in random media, vol. 2 (Academic press New York, 1978).

19. L. G. Henyey and J. L. Greenstein, “Diffuse radiation in the Galaxy,” The Astrophys. J. 93, 70–83 (1941). [CrossRef]  

20. S. L. Jacques, C. Alter, and S. A. Prahl, “Angular dependence of HeNe laser light scattering by human dermis,” Lasers Life Sci 1, 309–333 (1987).

21. S. Leyre, Y. Meuret, G. Durinck, et al., “Estimation of the effective phase function of bulk diffusing materials with the inverse adding-doubling method,” Appl. Opt. 53(10), 2117–2125 (2014). [CrossRef]  

22. G. L. Rogers, “Optical dot gain in a halftone print,” J. Imaging Sci. Technol. 41(6), 643–656 (1997). [CrossRef]  

23. Q. Fang and D. A. Boas, “Monte Carlo Simulation of Photon Migration in 3D Turbid Media Accelerated by Graphics Processing Units,” Opt. Express 17(22), 20178–20190 (2009). [CrossRef]  

24. H. C. van de Hulst, “Multiple light scattering,” (1980).

25. R. Giovanelli, “Reflection by semi-infinite diffusers,” Opt. Acta: Int. J. Opt. 2(4), 153–162 (1955). [CrossRef]  

26. L. Gevaux, L. Simonot, R. Clerc, et al., “Evaluating edge loss in the reflectance measurement of translucent materials,” Appl. Opt. 59(28), 8939–8950 (2020). [CrossRef]  

27. L. Ma, J. Zhai, and C. Wang, “Investigation of the single scattering approximation through direct electromagnetic scattering simulation,” OSA Continuum 4(9), 2496–2509 (2021). [CrossRef]  

28. S. L. Jacques, “Corrigendum: Optical properties of biological tissues: a review,” Phys. Med. Biol. 58(14), 5007–5008 (2013). [CrossRef]  

29. E. Salomatina, B. Jiang, J. Novak, et al., “Optical properties of normal and cancerous human skin in the visible and near-infrared spectral range,” J. Biomed. Opt. 11(6), 064026 (2006). [CrossRef]  

30. S.-H. Tseng, C.-K. Hsu, J. Yu-Yun Lee, et al., “Noninvasive evaluation of collagen and hemoglobin contents and scattering property of in vivo keloid scars and normal skin using diffuse reflectance spectroscopy: pilot study,” J. Biomed. Opt. 17(7), 0770051 (2012). [CrossRef]  

31. C. Roze, T. Girasole, G. Gréhan, et al., “Average crossing parameter and forward scattering ratio values in four-flux model for multiple scattering media,” Opt. Commun. 194(4-6), 251–263 (2001). [CrossRef]  

32. W. E. Vargas, J. Wang, and G. A. Niklasson, “Scattering and absorption cross sections of light diffusing materials retrieved from reflectance and transmittance spectra of collimated radiation,” J. Mod. Opt. 67(11), 974–991 (2020). [CrossRef]  

33. W. E. Vargas, J. Wang, and G. A. Niklasson, “Effective backscattering and absorption coefficients of light diffusing materials retrieved from reflectance and transmittance spectra of diffuse radiation,” J. Mod. Opt. 68(12), 605–623 (2021). [CrossRef]  

34. W. E. Vargas, D. M. Jiménez, and M. L. Montero, “Model to subtract contributions of scattered radiation from measured direct transmittance and specular reflectance by light diffusing materials,” Appl. Opt. 61(34), 10197–10206 (2022). [CrossRef]  

35. L. F. Gate, “Comparison of the Photon Diffusion Model and Kubelka-Munk Equation with the Exact Solution of the Radiative Transport Equation,” Appl. Opt. 13(2), 236–238 (1974). [CrossRef]  

36. S. N. Thennadil, “Relationship between the Kubelka-Munk scattering and radiative transfer coefficients,” J. Opt. Soc. Am. A 25(7), 1480–1485 (2008). [CrossRef]