Real-time rendering of layered materials with anisotropic normal distributions

This paper proposes a lightweight bidirectional scattering distribution function (BSDF) model for layered materials with anisotropic reflection and refraction properties. In our method, each layer of the materials can be described by a microfacet BSDF using an anisotropic normal distribution function (NDF). Furthermore, the NDFs of layers can be defined on tangent vector fields, which differ from layer to layer. Our method is based on a previous study in which isotropic BSDFs are approximated by projecting them onto base planes. However, the adequateness of this previous work has not been well investigated for anisotropic BSDFs. In this paper, we demonstrate that the projection is also applicable to anisotropic BSDFs and that the BSDFs are approximated by elliptical distributions using covariance matrices.


Introduction
In the last several decades, the visual quality of computer graphics has been improved significantly due to the long-standing efforts of both the research and industrial communities. In particular, success in reflectance modeling has enabled representation of a surprisingly wide variety of real-world materials in computer graphics.
Among such materials, those comprising of thin layers of different material components have recently attracted much attention due to the demand for surface-painted man-made objects. For example, a characteristic appearance of cars is generated by the surface coating process in which the car body is coated many times with different types of paints.
While accurate representation [1,2] and sampling of light transport paths [3] for layered materials have been proposed in the context of offline rendering, light transport in layered materials is usually approximated using analytic models particularly in real-time rendering. Weidlich and Wilkie [4] and the extension of their work by Elek [5] linearly combined the bidirectional scattering distribution functions (BSDFs) of layers using a transmission factor. Guo et al. [6] extended normal distribution functions (NDFs) using von Mises-Fisher (vMF) distributions to consider multiple reflection lobes and internal scattering. However, the vMF distributions cannot capture the heavy tails of directional distributions, which is often required for modeling metallic materials. Recently, Belcour [7] considered directional statistics of light rays by projecting the directional distribution on a base plane. The reflection or refraction property of each layer is then defined as an operator of changing the projected distribution. He referred to the operator as an atomic operator. Although the atomic operator was practically simple and powerful, its applicability to anisotropic reflection and refraction has not been well investigated.
In this paper, we extend the atomic operator for anisotropic reflection and refraction properties of layers. The previous method limited its application to isotropic reflection due to the use of a scalar value to define the variance of an isotropic distribution. We replace this scalar variance with a 2 × 2 covariance matrix to define the anisotropy of the distribution.
However, this extension is non-trivial because changes in the shapes of directional distributions have not been well investigated for rough boundary surfaces between layers with anisotropic scattering properties. Even in the most related study [8], the changes for reflection and refraction are calculated only on the center of directional distributions. In contrast, we derive the distribution shapes for entire directional distribution by considering the coordinate transform between those for NDFs and projected directional distributions. We implement this extended atomic operator for anisotropic reflections/refractions on a real-time rendering system [9] by following a publicly available implementation of the previous method [10]. The experimental results demonstrate that our extension synthesizes almost identical appearances to those obtained by offline Monte Carlo path tracing while its computational overhead from the previous method is as small as only 2.5%.

Background
In the original method [7], a behavior of light interaction with layered materials was represented by energy of light e and two statistical parameters, that is, the mean μ ∈ [−1, 1] 2 and variance σ ∈ [0, ∞] of the distribution of light directions. The symbol σ represents variance rather than standard deviation following the original paper [7]. The property of a surface between two neighboring layers is defined by three functions each of which modifies one of the three parameters above. For rough reflection and refraction, the parameters are transformed as follows: Reflection : Refraction : where h(α) = α 1.
where ω i ∈ S 2 denotes an incident direction; ω t ∈ S 2 refers to a refracted direction; n ∈ S 2 denotes a surface normal; (e i , μ i , σ i ) refer to the parameters of incident light; (e {R,T } , μ {R,T } , σ {R,T } ) denote the parameters of reflected or transmitted light; α ∈ [0, 1] and η refer to the roughness parameter and relative refractive index on a boundary surface between layers, respectively. F GD ∞ represents an integral of the product of Fresnel term F , shadowingmasking function G, and NDF D. In the original method, Belcour [7] precomputed the F GD ∞ values while considering multiple scattering effects [11] and stored them in a lookup table. On the other hand, another implementation in Unity [10] ignored multiple scattering effects and approximated the integral using a simple product of F and G for a direction of perfect reflection or refraction. For detailed definitions of function h(α) and roughness scaling factor s, refer to the original paper [7]. By successively applying the above transformations by the layers, we can obtain e q , μ q , and σ q of outgoing light for a configuration q of successive light interactions. For instance, q = T RT represents a transmission-reflection-transmission path. Let Q be a set of valid sequences of light interactions. Then, a bidirectional reflectance distribution function (BRDF) ρ is defined as follows: where shadowing-masking function, and reflect(μ q ) represents the direction of perfect reflection for μ q .
The above formulas are only applicable to isotropic reflection and refraction because the variance is modeled with a single variance parameter σ to define a radially symmetric distribution.

Layered materials with anisotropic normal distributions
The proposed method is an extension of Belcour's method [7] which approximated BSDFs by projecting them onto the base plane. The previous study restricted their applicability only to isotropic NDFs. In contrast, the proposed method extends the approach to anisotropic NDFs, as shown in Fig. 1.

Fig. 1
Our method works for layers with anisotropic NDFs defined on varying tangent vector fields, whereas the previous method [7] can be applied only to isotropic NDFs.

Covariance of projected distribution
To represent anisotropic BSDFs projected on the base plane, we employ a 2 × 2 covariance matrix Σ rather than a scalar variance σ. However, the relationship between the tangent vector field and the covariance matrix is non-trivial. Let t x ∈ S 1 and t y ∈ S 1 be the tangent and binormal vectors, respectively, that are the orthogonal vectors on the 2D local coordinate system P of the tangent vector field. When an NDF is projected on the plane, the principal directions of the projected distribution coincide with t x and t y following the definitions of GGX and Beckmann distributions (see Fig. 2).
Because an anisotropic BSDF can be approximated by an anisotropic spherical Gaussian [12], its projection to the base plane is also approximated by an anisotropic Gaussian function on the region near the distribution center, as illustrated in Figs. 2 and 3. This preservation of projected elliptic shapes is satisfied for both BRDF and BTDF with varying incident zenith angles as can be seen in Fig. 3.  Next, let us consider the relationship between a halfvector h = (ω i + ω o )/ ω i + ω o and outgoing direction ω o for reflection. Let (x h , y h ) be the projection of h in P, and (x o , y o ) be the projection of ω o in P. As discussed in Appendix B of a previous study [8], we assume that ω i = (sin θ i , 0, cos θ i ) in the tangent space using an incident zenith angle θ i ∈ [0, π/2] and an azimuthal angle of zero without loss of generality. Then, the relationship between (x h , y h ) and (x o , y o ) can be written as follows: h Assuming x h , y h , and θ i are small enough and we can ignore the second-and higher-order terms of x h , y h , and sin θ i , we approximate the above Jacobian matrix as The same representation was derived by Stam [8] as an exact solution at the perfect reflection vector (which corresponds to μ R ). Unlike his solution, Eq. (4) is the approximation over the region near μ R . For refraction, we also approximate the Jacobian matrix J t using the same assumption, as follows: where θ t is the zenith angle for the direction of refraction ω t . For both cases, the Jacobian matrix is a simple scaling matrix. For the derivation, please refer to Appendix A.1.
Although the assumption of the small zenith angle θ i causes a large error in the grazing angle, we prioritize the simplicity of implementation over physical strictness. While this problem of the grazing angle was also observed in the previous study [7], a compromise is allowable in practice, as we demonstrate later. We also need to consider the effects of the Fresnel term, shadowing-masking function, and cosine term to define a BSDF [13]. Nevertheless, these effects are low-frequency and can be negligible when the roughness parameters are relatively small. In the following, we discuss the property of the coordinate transform using the above Jacobian matrices.
As we can see in Eqs. (4) and (5), the Jacobian matrices for both reflection and refraction are diagonal, and their two diagonal entries are equal. The diagonality implies that the directions of the orthogonal basis vectors of P are preserved, as shown in Fig. 2. Therefore, we only need to transform anisotropic roughness parameters (α x , α y ) ∈ [0, 1] 2 along the tangent vector t x and binormal vector t y to define a covariance matrix for the projected BSDF. The uniformity of the diagonal entries implies that the stretch of the variances along t x and t y depends on neither the definition of the tangent vector field nor the difference in the roughness parameters. Therefore, we transform the roughness parameters (α x , α y ) to corresponding scalar variances (σ x , σ y ) ∈ [0, ∞] 2 using h(α). Accordingly, the covariance matrix for a BSDF is given as for refraction For energy e and mean μ, we use the same representations as those in the previous study because the anisotropy of BSDFs does not affect these terms significantly. Therefore, we obtain an extended BSDF with anisotropic NDFs by substituting the above covariance matrix Σ into Eqs. (1) and (2). To build a global BRDF using the adding-doubling method as in the original paper [7], we take exactly the same procedure introduced in it.

Results and discussion
The following experiments were conducted on a computer with an Intel R Core TM i7-8700 3.2 GHz CPU and NVIDIA R GeForce R RTX 2080 Ti GPU. We use a two-layer material in which the bottom conductor layer is coated with a clear dielectric layer unless otherwise specified. The formulas for two-layer materials are obtained by the adding-doubling method, and their derivations are described in Appendix A.2. We implement the proposed method using Marmoset Toolbag 3 [9].
In the rendering pipeline, we follow an approximation for F GD ∞ in the implementation of Unity [10] to avoid the lookup table being memory consuming for anisotropic materials. In addition, we calculate an average covariance matrix, which can differ from channel to channel, for three color channels following the public implementation of the previous studies [7,10].
While we mainly show the results of image-based lighting, the computation time of our method for a trivial scene with a directional light is 1.01 ms, which is sufficiently short for real-time applications such as interactive material editing. Figure 4 shows the rendering results obtained using our method for various layered materials with isotropic/anisotropic NDFs defined on the same/different tangent vector fields. These results include only direct illumination from environment maps. For this image-based lighting, we compute the Monte Carlo integration using visible NDF importance sampling [14] for each term of Eq. (3). We compare the computation time between Belcour's method [7] for isotropic BSDFs and our extension to anisotropic BSDFs. In this comparison, we evaluate both methods by sampling the BRDFs in Eq. (3) separately to focus on the overhead incurred by changing scalar variance σ to our covariance matrix Σ. Figure 5 shows the rendering time using varying numbers of samples. Although our method increases the ALU overhead and register pressure, this experimental result demonstrates that the performance degradation when using our method for anisotropic BSDFs is negligible.
The visual comparisons of the results for different layer configurations and different tangent vector directions are shown in Figs. 6 and 7, respectively. In these figures, pixel-wise root-mean-square errors (RMSEs) are visualized in the column to the right. We find that the RMSEs are rather large on the rim regions of the sphere, where the viewing angles are comparatively small. However, the overall computation time to render a single frame with 2048 samples per pixel was 587.8 ms for naive simulation to obtain reference images, whereas that for our method is only 45.1 ms. Additional results using different Fig. 5 The chart to the left compares the computation time of Belcour [7] for isotropic BSDFs and our extension to anisotropic BSDFs. The chart to the right shows the additional computation time in percentage terms of our extension over Belcour's method. roughness parameters and rotation angles for the local coordinate system appear in the Electronic Supplementary Material.
Our extension to anisotropic BSDFs does not depend on the number of material layers. Therefore, it is also available with materials with more than three layers as shown in Fig. 8(b). In this material, the layered material shown in Fig. 8(a) (which is equivalent to Fig. 4(d)) is further coated by a smooth dielectric layer. In addition, the atomic operators for anisotropic reflection/refraction can be combined also with those for isotropic ones. For instance, Fig. 8(c) shows a rendering result for a layered material including an interface of a participating medium. In Belcour's method, he  assumed the medium exhibits only forward scattering, and approximated the change of the directional distribution using an atomic operator as well. In our result in Fig. 8(c), we inserted this forward scattering layer in the middle of dielectric and conductor layers of the material shown in Fig. 8(a). For the forward scattering layer, an atomic operator is calculated by the parameter of the medium, i.e., scattering coefficient k s , absorption coefficient k a , and average cosine of the scattering angle g (see the literature [15]). For the details, refer to Belcour's original paper [7].
While our method renders anisotropic layered materials in interactive frame rates by extending Belcour's method [7], it is inevitable that our method inherits the limitations of this previous method. For example, for layers with very high roughness and indexes of refraction, the directional distributions projected on the base plane may not be elliptical. Despite this fact, the applicability of our method to a broad range of layered materials is beneficial to practical graphics production.

Conclusions
In this paper, we introduced a real-time approach for rendering layered materials wherein the layers are modeled by anisotropic NDFs defined on varying tangent vector fields. The proposed method is easily implemented on the top of the original approach proposed by Belcour [7] for isotropic NDFs, and works with minor additional computation cost.

A.1 Derivation of Jacobian matrices
To derive Jacobian matrices, we partly followed the derivation by Stam [8]. Different from his derivation, we derived an approximate solution for the Jacobian matrices over the region near to the center of the NDF, while Stam derived the exact solution only at the center. Without loss of generality, we can assume incident direction ω i as (θ i , 0). Let ω r and ω t be directions for reflection and refraction, respectively. We denote the directions ω i , ω r , ω t , and h as follows: ω i = (sin θ i , 0, cos θ i ) ω r = (x r , y r , z r ) Let η be a relative refractive index between two interfaces, we can write ω r and ω t as follows: Using these equations, we can obtain the projected 2D coordinates (x r , y r ) and (x t , y t ) of ω r and ω t : h . Therefore, for reflection, the Jacobian matrix is obtained as in the main body of the paper. For refraction, the Jacobian matrix is calculated as follows: As we wrote in the main body of the paper, we assume x h , y h , and θ i are small enough and we can ignore the second-and higher-order terms of x h , y h , and sin θ i . Then, we can approximate J t as follows: Thus, the Jacobian matrix for refraction is also diagonal and its diagonal entries are the same.