1. INTRODUCTION

Predicting the visual appearance of objects by carrying out an acceptable computation effort is often a challenge because of the variety of materials and surface finishes and the complexity of the multiple optical phenomena occurring within the material layers. Accurate prediction is however, crucial for the digital design of objects. It requires modeling the multiple reflections and scattering of light at the interfaces and within the turbid media composing the material. A radiometric approach consists in writing the balance for the absorbed and scattered fluxes by introducing the properties of single scattering. In the general tridimensional case, multiple scattering is mostly solved by Monte Carlo methods where each individual photon event is described with probabilities. For specific material structures, especially stacks of planar layers, one can solve the radiative transfer equation (RTE) [1]. This integro-differential equation takes into account both the spatial position and the orientations of the incident and scattered fluxes. However, as it is complex to solve, many resolution techniques and approximations have been suggested. One common approximation is not to take into account explicitly lateral scattering within the material. Therefore, the equation depends only on one spatial dimension, namely the depth $z$ within the stack. Regarding the angular distribution of light, it can be discretized into $N$-annular solid angles. This $N$-flux model was first proposed by Mudgett and Richards [2] in the case of azimuthally isotropic scattering and generalized by Stamnes et al. [3] under the so-called discrete ordinate method. $N$ can exceed 20 [2] but small $N$ values present the advantage of simple expressions for the reflectance and the transmittance factors. For $N=2$, only two hemispherical fluxes with constant radiance propagate toward positive and negative $z$. In that case, the RTE has analytical solutions, which are well known as the Kubelka–Munk formulas [4,5]. This two-flux model was also extended to determine the diffuse reflectances and transmittances of stacks of scattering layers [6]. Whereas the Kubelka–Munk model is the result of the continuous integration of the RTE, the Kubelka 1954 model can be interpreted as the corresponding discrete summation. The correspondences between the continuous and discrete two-flux approaches have already been discussed [7,8].

However, considering two diffuse fluxes is not possible when the incident light is collimated and a part of it becomes diffuse. The four-flux model with two additional collimated fluxes propagating perpendicularly to the planar layers toward positive and negative $z$ directions improves the reflectance and transmittance predictions in case of collimated illumination. Resolutions of the RTE according to the four-flux approach with various boundary conditions were proposed by Beasley et al. [9], Mudgett and Richards [2], and Ishimaru [10]. The formulation proposed by Maheu et al. [11,12] became the main reference when compared with exact calculations for specific cases [13] or after comparisons with Monte Carlo simulations [14]. In their formulation Maheu et al. introduced an average path-length parameter that can take values from 1 for a collimated flux to 2 for a perfect isotropic radiation. They also introduce a forward scattering ratio. The determination of these parameters was discussed by several authors [15–20]. The four-flux model can be expressed by using a matrix formalism as suggested by Rozé et al. in the case of multilayer [21]. Recent formulations [22,23] also enable predicting interface effects. The four-flux model can be used for various scattering systems as illustrated by recent publications [24–27]. However, even if the four-flux approach is much easier to use than more elaborate models, the simplicity of the two-flux approximation is still often preferred. For this reason intermediate models between two flux and four flux have also been proposed [28–30].

In the present study, we use the four-flux model without specifically focusing on the resolution of the RTE. In Section 2 we adopt a matrix formalism to calculate the reflectance and transmittance factors of a superposition of optical components (interfaces and propagating media). This approach can be seen as an extension of the Kubelka model [6] to four fluxes. The main contribution of this study, presented in Section 3, consists in adapting the four-flux matrix model in order to generate families of bidirectional scattering distribution factors (BSDF) by considering an incident collimated flux in any direction and directional diffuse fluxes exiting the material. The bidirectional transfer matrices are described for the particular cases of highly scattering (Lambertian) and of nonscattering (transparent) components in Section 4. We consider flat or rough dielectric interfaces in Section 5. The complete method is finally presented in Section 6 to determine the BSDF in the special cases where a flat or a rough interface is at the top of an opaque Lambertian background.

2. FOUR-FLUX MATRIX MODEL

The four-flux model considers a parallel planar structure of material and therefore reduces the radiative-transfer equation to a problem with one spatial dimension. It can be presented as a special case of the $N$-flux model [2] where the radiation field for each position in the stack of layers is composed of two collimated beams $I_c$ and $J_c$ and two isotropic diffuse beams $I_d$ and $J_d$. The fluxes propagate perpendicularly to the plane, forward ($I_c$ and $I_d$) and backward ($J_c$ and $J_d$).

A multilayer material can be described as a succession of interfaces and media. Each component of the stack, interface or medium, gives rise to flux transfers: front side reflectance $r$, back side reflectance $r^{\prime }$, forward transmittance $t$, and backward transmittance $t^{\prime }$. They can be collimated-to-collimated (label cc), diffuse-to-diffuse (label dd), or collimated-to-diffuse (label cd) transfers. Figure 1 represents the flux transfers for a stack of two components.

 

Fig. 1. Flux transfers between two components represented by thin arrows. Bold arrows correspond to fluxes.

Download Full Size | PDF

Let us consider the component labeled $k$. The fluxes labeled by superscripts $k$ and $k-1$ are related according to the following equations where, for the sake of writing simplicity, we omit the label $k$ in the transfer factors:

A. Matrix Formulation

The system of Eq. (1) can be presented under two possible matrix equations. The first matrix equation used to solve the RTE [31] in a multi-angle approach is shown in Appendix A. We consider here the second matrix equation, which focuses on the transfer nature (collimated-to-collimated, diffuse-to-diffuse, and collimated-to-diffuse) and is therefore specific to the four-flux approach,

B. Matrix Multiplication

In order to obtain the transfer matrix representing a stack of components, the components’ individual transfer matrices are multiplied by respecting the stacking order of the components. With two components characterized by matrices ${\mathbf{M}}_1$ and ${\mathbf{M}}_2$, from front to back, the transfer matrix of the two components together is

Given the complex expressions of the collimated-to-diffuse reflectance and transmittance [Eqs. (10) and (11)], even with two components only, the matrix formalism is much more convenient. However, the matrix computation is valid only when the following condition is satisfied for each component:

3. EXTENSION FOUR-FLUX MATRIX MODEL TO EXPRESS THE BSDF OF COMPONENT STACKS

In its original expression [11], the four-flux model assumes collimated and isotropic hemispherical diffuse fluxes propagating perpendicular to the stack of layers. In order to define BSDF models from the four-flux matrix method presented in Section 2, some adaptations are needed. The BSDF expresses the bidirectional reflectance and transmittance distribution functions (BRDF and BTDF), for which the incident illumination is assumed to be a unique collimated flux in any direction $\mathbf{i}$ of the upper hemisphere not only at the normal incidence. Moreover, the diffuse fluxes exiting the material from the border components of the stack are not necessarily assumed Lambertian and can therefore depend on the output direction $\mathbf{o}$ of the upper hemisphere for BRDF and of the lower hemisphere for BTDF.

The radiometric definitions and relations used in this section are detailed in the literature, for example in [33,34].

A. BSDF Configuration

According to the definition of the BSDF, the incident illumination is assumed to be a collimated flux in the incident direction $\mathbf{i}$. There is no incident diffuse flux or upward incident flux. In the case of the system represented in Fig. 1, this means that $I_d^0=J_c^2=J_d^2=0$. The scattered light is captured in every direction $\mathbf{o}$. Figure 2 explains the notations.

 

Fig. 2. Useful notations for defining the BSDF.

Download Full Size | PDF

B. Collimated Fluxes Related to the Incident Direction

The collimated incident beam is defined for a freely chosen orientation $\mathbf{i}$ within the upper hemisphere. The collimated fluxes after multiple reflections and transmissions remain within the incident plane. Knowing the different refractive indices, the directions of the downward fluxes are defined from the incident direction $\mathbf{i}$ and by Snell’s refraction law. The direction of each upward flux is deduced from the direction of the corresponding downward flux according to Snell’s reflection law.

C. Directional Diffuse Output Transfers

In the original four-flux model, the diffuse fluxes are assumed to be hemispherical with angle-independent radiance (i.e., Lambertian). In our approach we make an exception for the diffuse fluxes exiting the material from the first or from the last component of the stack. Therefore, the reflectance and transmittance factors to be chosen depend on the position of the component in the stack. Directional light that exits the stack of components is diffused according to the BRDF and the BTDF of the bordering components. Table 1 defines the collimated-to-diffuse and diffuse-to-diffuse transfer factors to be used according to the definitions by Nicodemus et al. [33]. For the sake of simplicity, we specify vector $\mathbf{i}$ only when the incident light is collimated and/or vector $\mathbf{o}$ only when light is captured in one direction (see Tables 1 and 2).

Table 1. Collimated-to-Diffuse and Diffuse-to-Diffuse Reflectance and Transmittance Factors of a Component According to its Position in the Component Stack^a

View Table | View all tables in this article

Table 2. Expressions of the Different Reflectance or Transmittance Factors in Function of the BRDF or BTDF $\mathsf{f}$, where ${\int }_{\mathsf{2}\mathsf{\pi }}\mathsf{f}\mathsf{cos}\mathsf{\theta }\mathsf{d}\mathsf{\omega }={\int }_{\mathsf{\theta }=\mathsf{0}}^{\mathsf{\pi }/\mathsf{2}}{\int }_{\mathsf{\phi }=\mathsf{0}}^{\mathsf{\pi }/\mathsf{2}}\mathsf{f}\mathsf{cos}\mathsf{\theta }\mathsf{sin}\mathsf{\theta }\mathsf{d}\mathsf{\theta }\mathsf{d}\mathsf{\phi }$

View Table | View all tables in this article

By knowing the BRDF and the BTDF of a component, we can define its reflectance or transmittance factors according to the formulas given by Table 2 [33]. For a stack of components, according to Table 1 and after the matrix multiplications, one can deduce the BRDF $f_r(\mathbf{i},\mathbf{o})$ and BTDF $f_t(\mathbf{i},\mathbf{o})$ of the stack

Moreover, in the case where part of the incident collimated flux remains collimated after exiting the border components, the specular reflectance $r_{\mathrm{cc}}(\mathbf{i})$ and the regular transmittance $t_{\mathrm{cc}}(\mathbf{i})$, given by Eq. (6) with $xx=\mathrm{cc}$, have to be added to the BSDF formally by using Dirac delta functions [35].

In contrast with more elaborated models based on the multi-angle scattering approach [36], we assume that the diffuse fluxes lose their directionality within the stack. In many systems this limitation has a weak influence when at least one component of the stack is sufficiently scattering to make the assumption acceptable. But even when it is not the case, the directionality of the collimated fluxes is preserved within the stack, and the directionality of the diffuse fluxes is conserved at the extreme components of the stack. A more restrictive limitation, but intrinsic to the four-flux approach, is the fact that the angular spreading of the collimated fluxes is not rendered.

4. TRANSFER MATRICES FOR LAMBERTIAN AND NONSCATTERING COMPONENTS

Among the optical components, the Lambertian scattering as well as the nonscattering components are interesting limit cases. Their presence in a stack enables important simplifications of the prediction method, especially when the border components of the stack are nonscattering.

A. Lambertian Component

For a Lambertian component, light is uniformly scattered over the hemisphere independent of the orientation of the incident light. It is worth noting that it is an ideal case [37]. Consequently, an incident collimated light is entirely transformed into diffuse light and $t_{\mathrm{cc}}=r_{\mathrm{cc}}=0$ [Fig. 3(a)]. However, to fulfill the condition, Eq. (12), we artificially define a transmittance $t_{\mathrm{cc}}$, which we will set to zero hereafter, and the following matrix ${\mathbf{M}}_{\mathrm{cc}}$:

 

Fig. 3. Flux transfers (a) for a Lambertian component, (b) for a nonscattering component.

Download Full Size | PDF

B. Nonscattering Component

Without scattering, the collimated fluxes cannot be transferred into diffuse fluxes [Fig. 3(b)]. Therefore,

For the special case of a nonscattering component, the diffuse-to-diffuse reflectances or transmittances are the integrals of the collimated beams in all directions of the upper or lower hemisphere, and can then be directly expressed from the corresponding collimated-to-collimated transfer factors

It is worth noting that the transfer matrix for the superposition of nonscattering components verifies Eq. (18). The resulting component is therefore also nonscattering. The reflectances and transmittances can be determined by operating independently with the $2\times 2$ transfer matrices ${\mathbf{M}}_{\mathrm{cc}}$ and ${\mathbf{M}}_{\mathrm{dd}}$. However, the diffuse reflectances and transmittances obtained by diffuse-to-diffuse transfer matrix multiplication [Eqs. (9) with $xx=\mathrm{dd}$] are crude approximations for nonscattering components. One needs to calculate first the collimated-to-collimated reflectances and transmittances [Eqs. (9) with $xx=\mathrm{cc}$] and integrate them over the hemisphere [Eq. (19)].

C. Scattering Components Surrounded by Nonscattering Components

When several components are superposed, sub-stacks of nonscattering components are first regrouped and their corresponding transfer matrices are determined. The expressions of their diffuse-to-diffuse transmittance factors depend on if they are boundary component or not according to Table 1. In the case that the nonscattering sub-stack is the first component, the hemispherical-directional $t_{\mathrm{dd}}^{\prime }(\mathbf{o})$ is

If it is the last component of the stack, the hemispherical-directional $t_{\mathrm{dd}}(\mathbf{o})$ is

A nonscattering component presents the advantage to preserve the bidirectional reflectance factor, respectively, transmittance factor, when it is the first or the last component. Let us consider a nonscattering component (index 1) on any scattering component (index 2) as represented in Fig. 4(a). By using Eq. (11) and the properties of the nonscattering component 1 [Eqs. (18) and (20)], the resulting bidirectional reflectance factor $r_{\mathrm{cd}}(\mathbf{i},\mathbf{o})$ can be expressed in terms of the one of the component 2, $r_{\mathrm{cd}2}({\mathbf{i}}_2,{\mathbf{o}}_2)$:

 

Fig. 4. (a) Bidirectional reflectance factor $r_{\mathrm{cd}}(\mathbf{i},\mathbf{o})$ of a nonscattering component (index 1) on a scattering component (index 2). (b) Bidirectional transmittance factor $t_{\mathrm{cd}}(\mathbf{i},\mathbf{o})$ of a scattering component (index 1) on a nonscattering component (index 2).

Download Full Size | PDF

With similar considerations, we can obtain the bidirectional transmittance factor for any scattering component (index 1) on a nonscattering component (index 2) as represented in Fig. 4(b):

5. TRANSFER MATRICES FOR FLAT AND ROUGH INTERFACES

As the first and the last components mainly influence the BRDF and BTDF of a layered material, the scattering responses of these two components must be analyzed carefully. These border components are most often interfaces. We present the corresponding four-flux matrices for a flat or a rough interface between two dielectric media labeled 0 and 1, with respective refractive indices $n_0$ and $n_1$. The relative refractive index of the interface is denoted as $n=n_1/n_0$.

A. Flat Interface

A flat interface is a nonscattering component, which does not enable any collimated-to-diffuse light transfer. Therefore, the $2\times 2$ matrix ${\mathbf{M}}_{\mathrm{cd}}$ is a zero matrix. The collimated-to-collimated transfers are given by the Fresnel formulas as functions of the incident direction $\mathbf{i}$. By calling $r_{\mathrm{cc}}=R_{01}(\mathbf{i})$ and $t_{\mathrm{cc}}=T_{01}(\mathbf{i})$, and by considering

According to Eqs. (20) and (21), if the flat interface is the first or the last component of the stack, the hemispherical-directional transmittance factors can be calculated as

B. Rough Interface

Any model describing the BSDF of a rough interface can be used in the four-flux model presented in this paper. We adopt the micro-facet model described by Walter et al. [35]. The roughness parameter σ becomes an additional index for all transfer factors. We assume that the incident collimated flux is completely converted into diffuse fluxes. Hence, the collimated-to-collimated transfers are assumed to be zero [Eq. (14)]. The bidirectional reflectance factor is defined as

All entries of the transfer matrix are reflectance and transmittance factors that can be calculated according to the position of the component within the stack (Table 1) and to the angular distribution of the incident flux on both faces (relations of Table 2). For example, the following equation gives the expression of the bihemispherical reflectance factor $r_{01\sigma }$ as a function of the directional-hemispherical reflectance factor $r_{01\sigma }(\mathbf{i})$:

6. FLAT OR ROUGH INTERFACE ON AN OPAQUE LAMBERTIAN BACKGROUND

As case studies, we consider a flat or a rough interface on the top of an opaque Lambertian background [Figs. 5(a) and 5(b)]. These systems involve only two components and the matrix calculations result in compact analytical relations.

 

Fig. 5. (a) Perfectly flat interface on a Lambertian background, (b) micro-facet rough interface on a Lambertian background, (c) distribution of interfaced Lambertian micro-facets.

Download Full Size | PDF

The calculation consists in multiplying the transfer matrix of the interface described in Section 5 by the transfer matrix of the opaque Lambertian background defined in Section 4.B. As a Lambertian component, the substrate matrix must verify Eqs. (14) and (16). The background diffuse reflectance is denoted as $\rho =r_{\mathrm{d}d}=r_{\mathrm{cd}}$. Although the substrate is opaque, we artificially assume that $t_{\mathrm{cc}}\ne 0$ and $t_{\mathrm{dd}}\ne 0$ [Eq. (12)] while $t_{\mathrm{cd}}=t_{\mathrm{cd}}^{\prime }=t_{\mathrm{cc}}^{\prime }=t_{\mathrm{dd}}^{\prime }=0$ and $r_{\mathrm{cd}}^{\prime }=r_{\mathrm{cc}}^{\prime }=r_{\mathrm{dd}}^{\prime }=0$. The transfer matrices of the opaque Lambertian background can then be written as

A. Flat Interface on a Lambertian Background

As the first components of the stack, the transfer matrices for the plane interface are [Eqs. (25)–(28)]

B. Rough Interface on a Lambertian Background

By using the notations introduced in Section 5, as well as Eqs. (4) and (5), the transfer matrices for the rough interface as the first component of the stack

After multiplying the rough-interface matrix, Eq. (37), by the opaque Lambertian background matrix, Eq. (33), we obtain

C. Interfaced Lambertian Facets

It is interesting to compare the configuration described by Eq. (39) and Fig. 5(b) with the one described by Fig. 5(c), which was developed in a previous work [44], with the same micro-facet slope distribution. The corresponding BRDF can be written as

It can be first noted that both models are equivalent for a flat interface $(\sigma =0)$ with the expression given in Eq. (36). Moreover, the first term $r_{01\sigma }(\mathbf{i},\mathbf{o})$, due to surface scattering, is similar for both Eqs. (39) and (40). It corresponds to a Cook–Torrance-like specular lobe [45] whose expression is given in Eq. (30). We therefore focus the comparison on the second term of Eqs. (39) and (40), due to volume scattering.

When a Lambertian background has a flat interface, its volume BRDF trends rapidly toward zero at grazing incident angles (see Fig. 6 for $\sigma =0$). This effect is reduced when considering a rough interface superposed over the Lambertian background [Fig. 6(a)]. The resulting volume BRDF tends to the one of a Lambertian reflector when the roughness increases. This BRDF is azimuthally isotropic, which is not the case for the distribution of interfaced Lambertian facets [Fig. 6(b)]. The variation of the volume BRDF depends on the observation angle. It decreases when the observation is from the part of the hemisphere containing the specular direction [${\theta }_o>0$ in Fig. 6(b)] but it increases for grazing angles toward the backscattering direction [${\theta }_o<0$ in Fig. 6(b)]. This difference between both models is striking even when there is no refractive index change $(n=1)$ between media. In that specific case, Fig. 5(b) [Eq. (39)] is equivalent to the flat Lambert background while Fig. 5(c) [Eq. (40)] corresponds to a distribution of Lambertian micro-facets and is therefore equivalent to an Oren–Nayar-like model [46].

 

Fig. 6. Volume BRDF in the incident plane, assuming a Beckmann distribution for $D_{\sigma }$ and the corresponding Smith shadowing masking function for $G$ [47,48], with ${\theta }_i=60°$ (backscattering direction at ${\theta }_o=-60°$ and specular direction at ${\theta }_o=60°$), $n=1.5$ and different roughness parameters σ. (a) Rough interface on a Lambertian background [second term of Eq. (39) with $r_{10\sigma }=r_{10}$], (b) distribution of interfaced Lambertian facets [second term of Eq. (40)].

Download Full Size | PDF

7. CONCLUSION

The methodology presented in this paper opens new perspectives to solve radiative transfer problems. We adapt the four-flux model by describing a material as a stack of discrete components. We introduce a $4\times 4$ transfer matrix that describes each component of the stack, collimated-to-collimated, diffuse-to-diffuse, and collimated-to-diffuse flux transfers both in transmittance and reflectance. Stacks of components are built by carrying out the corresponding matrix multiplications. Special matrices describe a Lambertian or a nonscattering component, a flat or a rough interface, and a border component given by its BRDF and BTDF. By construction, the $4\times 4$ transfer matrix can be reduced to $2\times 2$ matrices, and the two-flux models for either collimated-only or diffuse-only beams are special cases of this four-flux model. The use of four-flux is justified as soon as at least one component induces collimated-to-diffuse light transfers, and particularly for translucent materials where the two-flux models fail. We also extend the four-flux to obtain the BSDF of multilayer systems. A restriction of the model is the directionality loss of the diffuse fluxes within the layer stack. However, in most common cases, this limitation has a negligible impact because the model accounts for the directionality of the collimated fluxes at every position within the stack and for the directionality of the diffuse fluxes for the most external scattering components of the stack. Another restriction is that the four-flux approach does not describe the progressive angular broadening of the collimated beams. However, the method is easy to use and offers compact matrix expressions and quick computations. For computer graphics, it allows us to generate families of virtual BSDF for a superposition of components knowing the BSDF of each component.

We apply the method for a flat and a rough interface on an opaque Lambertian background. Although systems with only two components are relatively basic, they offer a wide variety of physical-based reflectance models from matte to glossy materials. We intend to apply this methodology to more complex systems where the compact matrix formalism will be an attractive feature.

APPENDIX A: ANOTHER TRANSFER MATRIX

The transfer matrix chosen in the paper focuses on the nature of the light transfers (labeled cc, dd, or cd). After rearrangement of the flux order for the vectors, another matrix relation alternative to Eq. (2) can be obtained:

(A1)$$\left(\begin{array}{cccc}{t}_{\mathrm{cc}}& 0& 0& 0\\ t_{\mathrm{cd}}& t_{\mathrm{d}d}& 0& 0\\ -r_{\mathrm{cc}}& 0& 1& 0\\ -r_{\mathrm{cd}}& -r_{\mathrm{dd}}& 0& 1\end{array}\right)\left(\begin{array}{l}{I}_c^{k-1}\\ I_d^{k-1}\\ J_c^{k-1}\\ J_d^{k-1}\end{array}\right)=\left(\begin{array}{cccc}1& 0& -r_{\mathrm{cc}}^{\prime }& 0\\ 0& 1& -r_{\mathrm{cd}}^{\prime }& -r_{\mathrm{dd}}^{\prime }\\ 0& 0& t_{\mathrm{cc}}^{\prime }& 0\\ 0& 0& t_{\mathrm{cd}}^{\prime }& t_{\mathrm{dd}}^{\prime }\end{array}\right)\left(\begin{array}{l}{I}_c^k\\ I_d^k\\ J_c^k\\ J_d^k\end{array}\right).$$

By introducing the $2\times 2$ matrices $\mathbf{T}$, $\mathbf{R}$, ${\mathbf{T}}^{\prime }$, and ${\mathbf{R}}^{\prime }$, the matrix relation can be rewritten by block

(A2)$$\left(\begin{array}{cc}\mathbf{T}& {\mathbf{0}}_{2,2}\\ -\mathbf{R}& {\mathbf{1}}_{2,2}\end{array}\right)\left(\begin{array}{l}{I}_c^{k-1}\\ I_d^{k-1}\\ J_c^{k-1}\\ J_d^{k-1}\end{array}\right)=\left(\begin{array}{cc}{\mathbf{1}}_{2,2}& -{\mathbf{R}}^{\prime }\\ {\mathbf{0}}_{2,2}& {\mathbf{T}}^{\prime }\end{array}\right)\left(\begin{array}{l}{I}_c^k\\ I_d^k\\ J_c^k\\ J_d^k\end{array}\right).$$

By assuming $t_{\mathrm{cc}}{t}_{\mathrm{dd}}\ne 0$, the matrix of the left side is inverted and left multiplied to the matrix of the right side:

(A3)$$\begin{gathered} \left(\begin{array}{l}{I}_c^{k-1}\\ I_d^{k-1}\\ J_c^{k-1}\\ J_d^{k-1}\end{array}\right)=\mathbf{M}\left(\begin{array}{l}{I}_c^k\\ I_d^k\\ J_c^k\\ J_d^k\end{array}\right),\text{with}\mathbf{M}=\left(\begin{array}{cc}{\mathbf{T}}^{-1}& {\mathbf{0}}_{2,2}\\ {\mathbf{RT}}^{-1}& {\mathbf{1}}_{2,2}\end{array}\right)\left(\begin{array}{cc}{\mathbf{1}}_{2,2}& -{\mathbf{R}}^{\prime }\\ {\mathbf{0}}_{2,2}& {\mathbf{T}}^{\prime }\end{array}\right) \\ =\left(\begin{array}{cc}{\mathbf{T}}^{-1}& -{\mathbf{T}}^{-1}{\mathbf{R}}^{\prime }\\ {\mathbf{RT}}^{-1}& {\mathbf{T}}^{\prime }-{\mathbf{RT}}^{-1}{\mathbf{R}}^{\prime }\end{array}\right)=\left(\begin{array}{cc}{\mathbf{M}}_{11}& {\mathbf{M}}_{12}\\ {\mathbf{M}}_{21}& {\mathbf{M}}_{22}\end{array}\right). \end{gathered}$$

Given the $2\times 2$ block matrices ${\mathbf{M}}_{ij}$ of the transfer matrix, the following matrix relations can be obtained:

(A4)$$\{\begin{array}{l}\mathbf{R}={\mathbf{M}}_{21}{\mathbf{M}}_{11}^{-1}\\ \mathbf{T}={\mathbf{M}}_{11}^{-1}\\ {\mathbf{R}}^{\prime }=-{\mathbf{M}}_{11}^{-1}{\mathbf{M}}_{12}\\ {\mathbf{T}}^{\prime }={\mathbf{M}}_{12}-{\mathbf{M}}_{21}{\mathbf{M}}_{11}^{-1}{\mathbf{M}}_{12}\end{array}.$$

This transfer-matrix formalism was classically used to solve the radiative transfer equation with $N\times N$ block matrices [31]. Let us note that both matrix expressions [Eqs. (2) and (A1)] are related by permutation matrices.

APPENDIX B: ENERGY PRESERVATION FOR A MICRO-FACET ROUGH INTERFACE

We consider a rough interface between two media with $n=1.5$ whose bidirectional transfer factors are given by the Eqs. (30) and (31). The function $D_{\sigma }$ is the Beckmann distribution and the corresponding shadowing masking function $G$ is the one described by Smith [47] and generalized by Bourlier et al. [48]. Figure 7 shows the different directional-hemispherical factors. As expected, the rougher the surface is, the more Lambertian these resulting factors are. The differences with the flat interface $(\sigma =0)$ are more pronounced at grazing incident angles from the less to the more refractive medium [Fig. 7(a)]. The differences are much more important in the opposite direction due to the effect of total reflection for a large part of the incident hemisphere [Fig. 7(b)].

 

Fig. 7. Directional-hemispherical reflectance (solid lines) or transmittance (dashed lines) factors of rough interfaces as a function of the incident angle ${\theta }_i$ for different roughness parameters $\sigma$ with $n=1.5$ (a) from medium 0 to medium 1, (b) from medium 1 to medium 0.

Download Full Size | PDF

From relations similar to Eq. (32), the different bihemispherical factors can be calculated and are represented in Fig. 8 in terms of the roughness parameter. The conservation of energy should give $r+t=1$. This is however, not the case because the micro-facet models do not account for interactions with multiple facets. The energy loss increases steadily with the roughness parameter. It is relatively weak from medium 0 to medium 1 (3% loss for $\sigma =0.6$). It has a larger impact on the diffuse transmittance $t_{01\sigma }$ than on the diffuse reflectance $r_{01\sigma }$ [Fig. 8(a)]. From medium 1 to medium 0, the loss is much more important (20% loss for $\sigma =0.6$). Due to the total reflections, it mainly impacts the diffuse reflectance $r_{10\sigma }$ [Fig. 8(b)]. Let us note that because of energy loss, $r_{01\sigma }$ continuously decreases with the roughness parameter while $t_{10\sigma }$ is almost constant.

 

Fig. 8. Bihemispherical factors of a rough interface in terms of the roughness parameter $\sigma$ with $n=1.5$ (a) from medium 0 to medium 1, (b) from medium 1 to medium 0.

Download Full Size | PDF

As a physical description of the interactions with multiple facets seems to be very complex, these effects are mostly ignored or empirically corrected. For example, Jakob et al. [36] suggest to reintroduce the energy loss as a diffuse radiation in reflection and transmission so that the energy is conserved. With the same idea, we suggest to first calculate $r_{01\sigma }$, respectively, $t_{10\sigma }$, and then deduce $t_{01\sigma }=1-r_{01\sigma }$, respectively, $r_{10\sigma }=1-t_{10\sigma }$, to preserve energy conservation. The internal diffuse reflectance $r_{10\sigma }$ can then be considered as independent of the roughness and can be approximated by the internal diffuse reflectance $r_{10}$ of a flat interface.

Acknowledgment

L. Simonot thanks the Poitou-Charentes region for founding a visiting fellowship at the EPFL, Laboratoire des Systèmes Périphériques, from June to August 2015.

REFERENCES

1. S. Chandrasekhar, Radiative Transfer (Dover, 1960).

2. P. S. Mudgett and L. W. Richards, “Multiple scattering calculations for technology,” Appl. Opt. 10, 1485–1502 (1971). [CrossRef]  

3. K. Stamnes, S. Chee Tsay, W. Wiscombe, and K. Jayaweera, “Numerically stable algorithm for discrete-ordinate-method radiative transfer in multiple scattering and emitting layer media,” Appl. Opt. 27, 2502–2510 (1988). [CrossRef]  

4. P. Kubelka and F. Munk, “Ein Beitrag zur Optik der Farbanstriche,” Zeitschrift für technische Physik 12, 593–601 (1931).

5. P. Kubelka, “New contributions to the optics of intensely light-scattering material, part I,” J. Opt. Soc. Am. 38, 448–457 (1948). [CrossRef]  

6. P. Kubelka, “New contributions to the optics of intensely light-scattering material, part II. Non-homogeneous layers,” J. Opt. Soc. Am. 44, 330–334 (1954). [CrossRef]  

7. M. Vöge and K. Simon, “The Kubelka-Munk and Dyck paths,” J. Stat. Mech. 2007, P02018 (2007). [CrossRef]  

8. M. Hébert and J.-M. Becker, “Correspondence between continuous and discrete two-flux models for reflectance and transmittance of diffusing layers,” J. Opt. A 10, 035006 (2008). [CrossRef]  

9. J. K. Beasley, J. T. Atkins, and F. W. Billmeyer, “Scattering and absorption in turbid media,” in Electromagnetic Scattering, R. L. Rowell and R. S. Stein, eds. (Gordon and Breach, 1967), pp. 765–785.

10. A. Ishimaru, Wave Propagation and Scattering in Random Media (Academic, 1978).

11. B. Maheu, J. N. Le Toulouzan, and G. Gouesbet, “Four-flux models to solve the scattering transfer equation in terms of Lorenz–Mie parameters,” Appl. Opt. 23, 3353–3362 (1984). [CrossRef]  

12. B. Maheu and G. Gouesbet, “Four-flux models to solve the scattering transfer equation. Special cases,” Appl. Opt. 25, 1122–1228 (1986). [CrossRef]  

13. G. A. Niklasson, “Comparison between four flux theory and multiple scattering theory,” Appl. Opt. 26, 4034–4036 (1987). [CrossRef]  

14. B. Maheu, J. P. Briton, and G. Gouesbet, “Four-flux model and a Monte Carlo code: comparisons between two simple and complementary tools for multiple scattering calculations,” Appl. Opt. 28, 22–24 (1989). [CrossRef]  

15. Y. P. Wang, S. W. Zheng, and K. F. Ren, “Four-flux model with adjusted average crossing parameter to solve the scattering transfer equation,” Appl. Opt. 28, 24–26 (1989). [CrossRef]  

16. W. E. Vargas and G. A. Niklasson, “Forward average path-length parameter in four-flux radiative transfer models,” Appl. Opt. 36, 3735–3738 (1997). [CrossRef]  

17. W. E. Vargas and G. A. Niklasson, “Generalized method for evaluating scattering parameters used in radiative transfer models,” J. Opt. Soc. Am. A 14, 2243–2252 (1997). [CrossRef]  

18. W. E. Vargas, “Generalized four-flux radiative transfer model average path-length parameter in four-flux radiative transfer models,” Appl. Opt. 37, 2615–2623 (1998). [CrossRef]  

19. C. A. Arancibia-Bulnes and J. C. Ruiz-Suarez, “Average path-length parameter of diffuse light in scattering media,” Appl. Opt. 38, 1877–1883 (1999). [CrossRef]  

20. C. Rozé, T. Girasole, G. Gréhan, G. Gouesbet, and B. Maheu, “Average crossing parameter and forward scattering ratio values in four-flux model for multiple scattering media,” Opt. Commun. 194, 251–263 (2001). [CrossRef]  

21. C. Rozé, T. Girasole, and A. G. Tafforin, “Multilayer four-flux model of scattering, emitting and absorbing media,” Atmos. Environ. 35, 5125–5130 (2001). [CrossRef]  

22. N. Dong, J. Ge, and Y. Zhang, “Four-flux Kubelka–Munk model of the light reflectance for printing of rough surface,” Proc. SPIE 7241, 72411I (2009). [CrossRef]  

23. M. Hébert, R. D. Hersch, and P. Emmel, “Fundamentals of optics and radiometry for color reproduction,” in Handbook of Digital Imaging, M. Kriss, ed. (Wiley, 2015), pp. 1021–1077.

24. F. El Haber, X. Rocquefelte, C. Andraud, B. Amrani, S. Jobic, O. Chauvet, and G. Froyer, “Prediction of the transparency in the visible range of x-ray absorbing nanocomposites built upon the assembly of LaF₃ or LaPO₄ nanoparticles with poly(methyl methacrylate),” J. Opt. Soc. Am. B 29, 305–311 (2012). [CrossRef]  

25. S. Bayou, M. Mouzali, F. Aloui, L. Lecamp, and P. Lebaudy, “Simulation of conversion profiles inside a thick dental material photopolymerized in the presence of nanofillers,” Polym. J. 45, 863–870 (2013). [CrossRef]  

26. K. Laaksonen, S.-Y. Li, S. R. Puisto, N. K. J. Rostedt, T. Ala-Nissila, C. G. Granqvist, R. M. Nieminen, and G. A. C. Niklasson, “Nanoparticles of TiO₂ and VO₂ in dielectric media: Conditions for low optical scattering, and comparison between effective medium and four-flux theories,” Sol. Energy Mater. Sol. Cells 130, 132–137 (2014). [CrossRef]  

27. L. Wang, J. I. Eldridge, and S. M. Guo, “Comparison of different models for the determination of the absorption and scattering coefficients of thermal barrier coatings,” Acta Mater. 64, 402–410 (2014). [CrossRef]  

28. J. W. Ryde, “The scattering of light by turbid media—Part 1,” Proc. Roy. Soc. (London) A131, 451–464 (1931).

29. H. Pauli and D. Eitel, “Comparison of different theoretical models of multiple scattering for pigmented media,” Colour 73, 423–426 (1973).

30. A. B. Murphy, “Modified Kubelka–Munk model for calculation of the reflectance of coatings with optically-rough surfaces,” J. Phys. D 39, 3571–3581 (2006). [CrossRef]  

31. R. Aronson and D. L. Yarmush, “Transfer-matrix method for gamma-ray and neutron penetration,” J. Math. Phys. 7, 221–237 (1966). [CrossRef]  

32. G. G. Stokes, “On the intensity of the light reflected from or transmitted through a pile of plates,” Proc. R. Soc. London 11, 545–556 (1860). [CrossRef]  

33. F. E. Nicodemus, J. C. Richmond, J. J. Hsia, I. W. Ginsber, and T. Limperis, “Geometrical consideration and nomenclature for reflectance,” J. Res. Natl. Bur. Stand. 160, 1–52 (1977).

34. M. Hébert and P. Emmel, “Two-flux and multiflux matrix model for color surface,” in Handbook of Digital Imaging, M. Kriss, ed. (Wiley, 2015), pp. 1233–1277.

35. B. Walter, S. R. Marschner, H. Li, and K. E. Torrance, “Microfacet models for refraction through rough surfaces,” in Proceeding of Eurographics Symposium on Rendering (2007), pp. 195–206.

36. W. Jakob, E. d’Eon, O. Jakob, and S. Marschner, ”A comprehensive framework for rendering layered materials, “ ACM Trans. Graph. 33, 1–12 (2014).

37. A. Kienle and F. Forschum, “250 years Lambert surface: does it really exist,” Opt. Express 19, 3881–3889 (2011). [CrossRef]  

38. D. B. Judd, “Fresnel reflection of diffusely incident light,” J. Natl. Bur. Stand. 29, 329–332 (1942).

39. R. Molenaar, J. J. ten Bosch, and J. R. Zijp, “Determination of Kubelka–Munk scattering and absorption coefficient,” Appl. Opt. 38, 2068–2077 (1999). [CrossRef]  

40. M. Elias, L. Simonot, and M. Menu, “Bidirectional reflectance of a diffuse background covered by a partly absorbing layer,” Opt. Commun. 191, 1–7 (2001). [CrossRef]  

41. F. C. Williams and F. R. Clapper, “Multiple internal reflections in photographic color prints,” J. Opt. Soc. Am. 43, 595–597 (1953). [CrossRef]  

42. J. L. Saunderson, “Calculation of the color pigmented plastics,” J. Opt. Soc. Am. 32, 727–736 (1942). [CrossRef]  

43. L. Simonot, M. Hébert, and R. D. Hersch, “Extension of the Williams–Clapper model to stacked nondiffusing colored coatings with different refractive indices,” J. Opt. Soc. Am. A 23, 1432–1441 (2006). [CrossRef]  

44. L. Simonot, “A photometric model of diffuse surfaces described as a distribution of interfaced Lambertian facets,” Appl. Opt. 48, 5793–5801 (2009). [CrossRef]  

45. R. L. Cook and K. E. Torrance, “A reflectance model for computer graphics,” in SIGGRAPH ‘81 Proceedings of the 8th Annual Conference on Computer Graphics and Interactive Techniques (1981), Vol. 15, pp. 307–316. [CrossRef]  

46. M. Oren and S. K. Nayar, “Generalization of the Lambertian model and implications for machine vision,” Int. J. Comput. Vis. 14, 227–251 (1995). [CrossRef]  

47. B. G. Smith, “Geometrical shadowing of a random rough surface,” IEEE Trans. Antennas Propag. 15, 668–671 (1967). [CrossRef]  

48. C. Bourlier, G. Berginc, and J. Saillard, “One and two-dimensional shadowing functions for any height and slope stationary uncorrelated surface in the monostatic and bistatic configurations,” IEEE Trans. Antennas Propag. 50, 312–324 (2002). [CrossRef]