Pierre Mézières
email: pierre.mezieres@irit.fr

Mathias Paulin
email: mathias.paulin@irit.fr

Efficient spherical harmonic shading for separable BRDF

Keywords: Computing methodologies → Rendering, Rasterization, Realtime simulation

1441ms 1441ms

: Computing per-pixel shading using Spherical Harmonics (SH) is a costly task that requires to compute a triple SH product when evaluating specular contribution. In this work we derive a new formulation of this product for separable BRDFs. Our approach reduces computation timings by two orders of magnitude and improves the quality of the final image (as measured by RMSE). Top: images obtained with our approach using 9 SH bands. Bottom: Reference images generated by sampling 10.000 rays for each light.

INTRODUCTION

Spherical Harmonics (SH) are a well-known basis of functions used to solve efficiently the rendering equation. For instance, Precomputed Radiance Transfer (PRT) [Sloan et al. 2002] pre-computes and stores the lighting using SH at the vertices of a mesh and interpolate the SH over the triangles to compute per-pixel shading. 3D scenes are nowadays made of very refined meshes, and thus require long pre-computation time and large memory space to store and interpolate the SH coefficients.

In addition, SH have limited ability to represent high frequencies with limited memory requirements. As such, they first have been used to store and compute diffuse lighting contribution when using separable BRDFs. Recent approaches have been proposed to use SH for the specular contribution, but at the cost of increasing the rendering time. The main objective of this work is to reformulate SH operations to speed up the shading computation.

We observed that diffuse and specular contributions are always treated separately in the literature, but in practice there is a clear computational similarity between them. In this work we demonstrate the gain obtained by taking advantage of this similarity, and expose an efficient algorithm to compute the SH product between any function and the clamped cosine one.

At first, we quickly recall the SH properties we need for this paper (Sec. 2). Then, we present how to compute the diffuse contribution from the specular one (Sec. 3) and the way we can compute efficiently the SH product between an arbitrary function and a clamped cosine (Sec. 4). We conclude by an evaluation of our proposal and explore the different use cases for both methods (Sec. 5).

BACKGROUND

Spherical Harmonics. Real spherical harmonics are defined by

𝑌 𝑚 𝑙 (𝜃, 𝜙) =          √ 2𝐾 |𝑚 | 𝑙 𝑃 |𝑚 | 𝑙 (cos 𝜃 ) sin(|𝑚|𝜙) 𝑚 < 0 𝐾 0 𝑙 𝑃 0 𝑙 (cos 𝜃 ) 𝑚 = 0 √ 2𝐾 𝑚 𝑙 𝑃 𝑚 𝑙 (cos 𝜃 ) cos(𝑚𝜙) 𝑚 > 0 (1)
where 𝐾 is a normalization factor given by

𝐾 𝑚 𝑙 = √ 2𝑙 + 1 4𝜋 (𝑙 -𝑚)! (𝑙 + 𝑚)! . ( 2
)
Projection to SH obey the following equation

F 𝑚 𝑙 = ∫ Ω 𝐹 (𝜔) 𝑌 𝑚 𝑙 (𝜔) d𝜔 . (3)
Double and triple product. In literature, we find two classical operations on SH referenced as the double and triple product. Double product computes a spherical convolution between two functions, and is defined as

𝐷 = ∫ Ω 𝐿(𝜔)𝐹 (𝜔) d𝜔 = ∫ Ω ∑ 𝑙1,𝑚1 L 𝑚1 𝑙1 𝑌 𝑚1 𝑙1 (𝜔) ∑ 𝑙2,𝑚2 F 𝑚2 𝑙2 𝑌 𝑚2 𝑙2 (𝜔) d𝜔 = ∑ 𝑙1,𝑚1 ∑ 𝑙2,𝑚2 L 𝑚1 𝑙1 F 𝑚2 𝑙2 ∫ Ω 𝑌 𝑚1 𝑙1 (𝜔)𝑌 𝑚2 𝑙2 (𝜔) d𝜔 = ∑ 𝑙,𝑚 L 𝑚 𝑙 F 𝑚 𝑙 .
(4)

The triple product computes the SH coefficients of the multiplication of two spherical functions:

T 𝑚 𝑙 = ∫ Ω ∑ 𝑙1,𝑚1 L 𝑚1 𝑙1 𝑌 𝑚1 𝑙1 (𝜔) ∑ 𝑙2,𝑚2 F 𝑚2 𝑙2 𝑌 𝑚2 𝑙2 (𝜔) 𝑌 𝑚 𝑙 (𝜔) d𝜔 = ∑ 𝑙1,𝑚1 ∑ 𝑙2,𝑚2 L 𝑚1 𝑙1 F 𝑚2 𝑙2 ∫ Ω 𝑌 𝑚1 𝑙1 (𝜔)𝑌 𝑚2 𝑙2 (𝜔)𝑌 𝑚 𝑙 (𝜔) d𝜔 (5)
Usually, the integral in Eq.5 is precomputed to generate an efficient code that apply the SH product.

Estimating diffuse and specular contributions. Let's consider a typical rendering case with a light function 𝐿, a visibility 𝑉 (cosine weighted according to surface normal) and a separable BRDF 𝑀. The diffuse contribution is computed by using the double product between L and V and then multiplying by the albedo. The specular contribution is given by computing the triple product between L and V and then apply the double product between this result and the BRDF specular lobe. This process requires to compute double and triple products on the same two functions, which makes numerous computations redundant.

ESTIMATING DIFFUSE FROM SPECULAR

In this section, we show how to deduce the diffuse from specular contribution by computing the double product from the triple product with a single multiplication. Let's call 𝑂, the unit function over the sphere:

𝑂 (𝜔) = 1 = ∑ 𝑙,𝑚 O 𝑚 𝑙 𝑌 𝑚 𝑙 (𝜔) . (6)
Due to SH properties, only coefficient O 0 0 is non null and equal to

𝑂 0 0 = ∫ Ω 𝑌 0 0 d𝜔 = √ 1 4𝜋 4𝜋 = 2 √ 𝜋 ≈ 3.5449077 . (7)
Intuitively, if we multiply two spherical functions, i.e. apply the triple product, we can obtain the convolution between those two by applying the double product between the result of the triple product and 𝑂. Looking at the first SH coefficient obtained by triple product T 0 0 , we place the SH 𝑌 0 0 out of the integral as it's constant.

T 0 0 = 𝑌 0 0 ∑ 𝑙1,𝑚1 ∑ 𝑙2,𝑚2 L 𝑚1 𝑙1 F 𝑚2 𝑙2 ∫ Ω 𝑌 𝑚1 𝑙1 (𝜔)𝑌 𝑚2 𝑙2 (𝜔) d𝜔 = 𝑌 0 0 𝐷 . (8)
By reversing, we obtain

T 0 0 2 √ 𝜋 = T 0 0 O 0 0 = 𝐷 . (9)
This signifies that when computing the triple product for two functions, the computation of the double product between the same functions is done with a single multiplication (𝐷 ≈ 3.5449077T 0 0 ) which corresponds to the double product between T and O.

EFFICIENT CLAMPED COSINE PRODUCT

We present a method that drastically decreases the number of arithmetic operations needed to compute the SH triple product of a function L and a clamped cosine.

As the cosine is a circularly symmetric function, only the zonal harmonics (ZH) (𝑚 = 0) coefficients are non null. By replacing the function F in the triple product equation (Eq. 5) by the equation of rotation of ZH functions [START_REF] Pike | Local, deformable precomputed radiance transfer[END_REF], we obtain: where 𝑛 is the normal surface and 𝐶 𝑚 1 𝑚 2 𝑚 𝑙 1 𝑙 2 𝑙 are the Clebsch-Gordan coefficients corresponding to the integral in Eq.5.

T 𝑚 𝑙 = ∑ 𝑙 1 ,𝑚 1 ∑ 𝑙 2 ,𝑚 2 L 𝑚 1 𝑙 1 𝑌 𝑚 2 𝑙 2 (𝑛) √ 4𝜋 2𝑙 2 + 1 F 0 𝑙 2 𝐶 𝑚 1 𝑚 2 𝑚 𝑙 1 𝑙 2 𝑙 (10)
This solution is more efficient than the classical SH triple product (Fig. 2). Moreover, as most of the cosine energy is caught on 3 SH bands, we can limit the cosine projection to this number of bands (𝑙 2 < 3) hence improving its evaluation. Nevertheless, there is a more efficient solution that does not imply a low frequency approximation of the cosine function. We can instead express the cosine in the world coordinate system and precompute it directly in the integral:

T 𝑚 𝑙 = ∑ 𝑙 1 ,𝑚 1 L 𝑚 1 𝑙 1 ∫ 𝑌 𝑚 1 𝑙 1 (𝜔) max(0, 𝜔 • 𝑧)𝑌 𝑚 𝑙 (𝜔)𝑑𝜔 (11)
When the light function is expressed in the same coordinate system, this solution provides fewer arithmetic operations. If light is not in the same coordinate system, adding a rotation step, as the zxzxz method [START_REF] Kautz | Fast, arbitrary BRDF shading for low-frequency lighting using spherical harmonics[END_REF], adds few operations and remains more efficient than other methods (Fig. 2). As noted for triple product, the parts of the equations into brackets (Eq.10 and Eq.11) are precomputed to generate the code that will apply the SH product.

EVALUATION AND USE CASES

As pointed out in the introduction, the goal of our methods is to compute the shading with SH using as few arithmetic operations and data storage as possible, and thus to reduce the computation time for rendering.

We compiled in Fig. 3 the different algorithms one can implement to compute shading using SH, considering four functions that appear in the rendering equation: light (𝐿), visibility (𝑉 ), brdf (𝑀) and the clamped cosine (cos). Paths that do not use visibility are particularly adapted to render fully dynamic scenes because shading do not depend on mesh-specific precomputations, unlike when visibility is pre-computed at each vertex. Using traditional PRT, clamped cosine is packed with visibility in a precomputation step for each vertex. Similar but less common, cosine can be packed with light instead of visibility. This path can be used with dynamic visibility that do not take cosine into account, as it is included in the light function. Finally, path at the bottom, where 𝐿, 𝑉 and cosine are projected separately on SH is the more general. This path can be taken if everything is dynamic, light and visibility, for instance when using many polygonal lights [START_REF] Wu | Analytic spherical harmonic gradients for real-time rendering with many polygonal area lights[END_REF]] and with dynamic visibility [START_REF] Zhong Ren | Real-time soft shadows in dynamic scenes using spherical harmonic exponentiation[END_REF]]. In such case cosine cannot be packed with either function. Our contributions should be used for all of these algorithms to reduce their computational cost.

Eq. for SH product Fig. 2 Table 1: RMSE (×10 2 ) for the SH product. Our method always generates the lowest rmse because the cosine is not directly approximated by SH (Eq. 11). The important difference in magnitude between figures is due to the nature of materials used, mainly diffuse (Fig. 1b) or glossy (Fig. 1c). Specular is computed up to 9 SH bands while diffuse is computed with our method (Eq. 9).

Eq. for SH product Fig. 2 Table 2: Timings (ms) with the different solution to apply the SH product with a clamped cosine. Specular is computed up to the corresponding band while diffuse is computed using our proposal (Eq. 9).

Figure 2 :

 2 Figure 2: Number of arithmetic operations needed to compute different SH product up to a certain band.

Figure 3 :

 3 Figure 3: Flowchart of typical SH rendering with separable BRDF. We annotated four paths: 1. paths we used for our evaluation, 2. new paths from our first contribution, 3. Irradiance map from [Ramamoorthi and Hanrahan 2001], 4. Traditionnal PRT.

Figure 4 :

 4 Figure 4: Clamped cosine projected on different number of SH bands. Even though the majority of the cosine energy is caught with few bands, it produces non-negligible ringing artefacts as it is clamped.

Siggraph Asia Technical Communications, 2021, Tokyo Pierre Mézières and Mathias Paulin

ACKNOWLEDGMENTS

Sponza and Gallery taken from the McGuire Computer Graphics Archive. Living room taken from BlendSwap. Elephant taken from Sumner and Popovic.

 4). Scene contains only one light behind the elephant. Specular contribution is always computed on 9 SH bands. Diffuse contribution is computed with various number of SH bands, or evaluated from specular computation (our method). Negative color of shading result are clamped to zero, it explains why inset for 7 bands is darker than 9 while rmse is still higher. Ringing problem occurs too when using visibility with SH because of the frequency cut introduced by the clamping of the southern hemisphere of the unit sphere.

For all our experiments, we use the magenta path with spherical lights sources where lighting is computed at each pixel. We evaluate our results on a RTX 2080 with all renderings done at resolution 1920 × 968. To compare, we compute the traditional SH product (Eq. 5) using the same path for specular and Irr. map path for diffuse. To render dynamic scenes, cosine is precomputed on ZH and rotated on the fly [START_REF] Pike | Local, deformable precomputed radiance transfer[END_REF]] before computing the product.

Estimating diffuse from specular Our method to compute double product from triple product allows to avoid the convolution needed to compute diffuse contribution, thus replacing 𝑁 2 multiplications, where 𝑁 is the number of SH bands, by a single multiplication. This method creates a new derivation on each path (Fig. 3). Moreover, diffuse is classically computed on 3 SH bands as the majority of the cosine energy is caught. This restriction produces non-negligible ringing artifacts, visible when light is only on one object's side (Fig. 5), and is clearly explained by looking at the behaviour of clamped cosine projection on SH (Fig. 4). Our method removes these artifacts by computing diffuse lighting on the same number of bands as specular at no extra cost (Fig. 5).

Efficient clamped cosine product The efficient SH product with cosine is well adapted to render dynamic scenes. It fits into paths where the SH product with cosine is a specific step. We have shown that with our method the number of arithmetic operations is greatly reduced (Fig. 2), as well as computation time (Table 2). All our results were computed using Eq.11 with all SH coefficients expressed in the same frame. We also demonstrated (Table 2, last row) that when coefficients are not expressed in the same frame, thus requiring an additional SH rotation, our proposal still greatly improves the computation time. This last case is useful, e.g., when computing SH gradients as in [START_REF] Wu | Analytic spherical harmonic gradients for real-time rendering with many polygonal area lights[END_REF]. Moreover, since the cosine is not approximated with SH (Eq. 11), the rmse is lower than the other methods (Table 1).

CONCLUSION AND FUTURE WORK

We have presented a general method to avoid redundant computations when computing both diffuse and specular contributions in SH shading. We shown that diffuse contribution can be easily computed with a single multiplication from the specular computation. We also shown that the SH product between light, or an arbitrary function, and a clamped cosine can be done with only a few arithmetic operations. Our proposals then improve the computation time up to two orders of magnitude when using 9 SH bands.

It could be envisioned to apply the same principle to other functions that we did for cosine, i.e. replace cosine by any other function in Eq.11. Nevertheless, while improving the computation time for the triple product, the final number of required operations will be hard to predict.

Algorithms that fit into our classification (Fig. 3) should also benefit from our approaches by using the blue paths appropriately.