Physically Based Rendering - Lighting

As soon as we modeled the surface’s physical properties that covered a certain range of material in real life, we would need to emit light onto them, in order to finally get the outcome radiance from the surface. If you take a look back at the rendering equation, the outcome radiance is just an integral of the income radiance over the semi-sphere around the normal. This gives us a fundamental assumption that with a brutal algorithm such as path tracing we could get a numerical solution for the problem. We could just model the geometry shape of the light source, then emit single light through all the possible directions and evaluate whether they hit a surface, until the overall energy diminished to a certain threshold. But in real-time rendering, the computational source currently we had in hand won’t be sufficient for such cost of the algorithm. We need to find the analytical replacement of the integral, or at least some cheaper numerical integrals.

Remapping of reality

Before we dive into the detail of the analytical representation of different light shapes, I’d introduce another type of unit system to measure energy and related quantities. If you remembered what I introduced in a previous post about radiometry, they are identically related. Thus the photometry it is, basically it’s a weighted version of radiometry, but with the respect to the nature of human eyes’ sensitivity. Because when we deal with the light in real life, typically the light sources we would get are bulbs, candles, flashlights, and our old friends, sun, we perceive the electromagnetic wave from them by the cell on our retina. And because of the non-linearly of human nature (of course!), our eyes have different sensitivities among different light wavelengths, so finally people invented the photometry quantities to better measure the light we actually “see”.

Here are the common photometry quantities we would occur in real-time rendering:

As you could see, it’s basically another version of quantity measurement for the electromagnetic radiation, but with more emphasis on the visible light range. Because the two types of retina cell - rods and cones have quite different responses for different intensity and wavelengths of lights, the actual mapping between radiometry to photometry requires to be done in two different versions. The one related to cones cell typically has a luminance range from ${10 nt}$ to ${10^8 nt}$ and we call it photopic vision in biology, and another one related to rods cell is called scotopic vision with the luminance range from ${10^{-3} nt}$ to ${10^{-6} nt}$. The cones cell could receive the chromatic of the world and the rods cell receives the silhouette of the object. But you may wonder, how much is one $nt$? The SI defines the fundamental unit of the light measurement by ${cd}$, and then the other units are its deduced units. And with the standard definition, one candela means “A $540.0154×10^{12}Hz$ monochromatic light emits $1/683$ Watt per steradian”, then we could say how much a nit it is by adding “per square meter” to it. The reason to choose such a reference wavelength is that one of the three types of cone cells is most sensitive to it, and if you took a look at the visible light spectrum you’d see it’s a greenish color. The magic number 1/683 came from a historical compatibility requirement to transit the old unit system to the modern one.

“…Radiometry deals purely with physical quantities, without taking account of human

perception. A related field, photometry, is like radiometry, except that it weights everything by the sensitivity of the human eye. The results of radiometric computations are converted to photometric units by multiplying by the CIE photometric curve, bell-shaped curve centered around 555 nm that represents the eye’s response to various wavelengths of light…” - pg. 271, ”Real-Time Rendering”, 4th Edition.

So with a certain wavelength $\lambda$, the mapping from radiometry to photometry is expressed like ${I_v(\lambda)}=683.002\overline{y}(\lambda)I_e(\lambda)$, or we could integrate over the entire visible light wavelength range to get ${I_v}=683.002\int_{380}^{780}\overline{y}(\lambda)I_e(\lambda)d\lambda$. The $\overline{y}(\lambda)$ here is the luminosity function, it’s a bell-shaped distribution function and has different values in photopic and scotopic. The photopic luminosity function would reach 1 in 555nm or 540 THz (the definition frequency of candela) and the scotopic luminosity function is 1 in around 510nm, they demonstrate the eyes’ perception ratio of different wavelengths. And we could furthermore deduce two quantities called luminous efficacy and luminous efficiency from the mapping function. Luminous efficacy $\eta$ is defined as ${\eta}=683.002\frac{\int_{380}^{780}\overline{y}(\lambda)I_e(\lambda)d\lambda}{\int_{380}^{780}I_e(\lambda)d\lambda}$ and the unit is ${\frac{lm}{W}}$, it represents the light source’s ability to emit visible light; And if we normalized the luminous efficacy with ${\frac{1W}{683lm}}$ then we would get luminous efficiency, which is simply defined as $\frac{\int_{380}^{780}\overline{y}(\lambda)I_e(\lambda)d\lambda}{\int_{380}^{780}I_e(\lambda)d\lambda}$. With the help of them we could easily evaluate whether a light source is wasting or saving energy, or to directly convert from the radiometric quantities to photometric version. And we could get that an ideal 540 THz black body would have a 683 ${\frac{lm}{W}}$ luminous efficacy or 100% luminous efficiency.

Since what we’ve been talking so far is all about the real-time rendering, we would use the discrete RGB space instead of a continuous spectrum, and then we could simplify the continuous luminosity function to a discrete version. And we could then give the user a few parameters like the radiometric quantities and the luminous efficiency to adjust the light source. But when we build a physically-based rendering pipeline we’d often go to find references in real life, and almost all the light sources like those bulbs sold in supermarket would print their photometric quantities on the package, it would be more convenient to use photometric quantities in the light system for the final users.

Integration matters

The more crucial part of a light system is how to model the geometry appearance of the light source. When dealing with a non-physically based rendering pipelines, there are often 3 types of analytical light we’d like to use, directional light, point light, and spot light. I’d not discuss any of them here because all of them are not modeling around the actual light source geometry but rather around the appearance of the lighting result, another reason is that there are tons of information online about how to implement them and I’d better not repeat again. I’d indicate them as punctual light sources later since the actual shape of them is infinitesimal. The physically correct light sources all have volume, or if we ignore one of the 3 dimensions with the respect of the nature of projection they still have the area. I’d indicate them as area light source later on. The typical analytical area light source has shapes like sphere, disk, tube and rectangular, or with a trivial shape which requires some advanced approximations.

If we took a look once more at the reflectance equation, let’s write a slightly different version $L_o=\int\limits_{\Omega+} f(v,l) V(l) L_i \langle n·l \rangle dl$ instead, the $V(l)$ here is a Heaviside function to indicate whether a light source is visible along a certain direction, and it describes the shape and the shadow of the light. As long as the shape of the light source is non-trivial, it’s impossible to find a general analytical solution for such hemisphere integral. The punctual light sources all assume the shape is infinitesimal so the integral would be simplified by an integral of cosine over hemisphere and finally we’d get $L_o=\pi f(v,l) \langle n·l \rangle$, that’s the reason why we could write Lo = rho * NDotL when using the simple Lambert diffuse BTDF $f(v,l) = \frac{\rho}{\pi}$, the $\pi$ is just perfectly canceled. For area light we cannot have such one-line-to-rule-them-all enjoyment, there are some practical solutions so far, that all of them would solve the problem to a certain level which is acceptable enough. I’d introduce them one by one below.

Form Factor

The first approach is that we could still try to solve the integral directly. We can’t numerically solve it in runtime, but if we convert it to the integral over the light source’s area instead, then it would bring some certainties because we would have known the shape or the actual area of the light source when shading the object. The rewritten version of the reflectance equation is $L_o=\int\limits_{A} f(v,l) L_i \langle n·l \rangle \frac{\langle n_a·-l \rangle}{r^2} dA$, here $n_a$ is the normal of $dA$, $r$ is the distance from the lighting point to $dA$. Because the effective area along the $dl$ is proportional to the normal orientation of $dA$ and the distance from the shading point to the $dA$, so here we introduced a $\frac{\langle n_a·-l \rangle}{r^2}$ factor.

And then we could move the BSDF out of the integral by assuming it’s not interleaved with light direction so much, for example in a Lambert diffuse BTDF case. Then what we’d solve is $L_o = f(v,l) \int\limits_{\Omega+} L_i \langle n·l \rangle dl$, or $L_o = f(v,l) E(n)$ that we could treat all the income luminance as the illuminance. The analytic method to integrate illuminance has been already developed in problem domain of energy transformation, such like radio transfer or heat transfer, since fundamentally they are all the same type of problems about how to calculate the energy transferred from one object to another, or specifically speaking in our case, from a surface to another one.

“…The view factor from a general surface A1 to another general surface A2 is given by: $\displaystyle F_{1\rightarrow 2}={\frac {1}{A_{1}}}\int _{A_{1}}\int _{A_{2}}{\frac {\cos \theta _{1}\cos \theta _{2}}{\pi s^{2}}}\,{d}A_{2}\,{d}A_{1}$…” [Wiki1]

The actual analytic formulae are more complicated but well-documented in [Mar14] and [HSM10], I’d only give each link of the commonly used geometry shapes below:

• Sphere: “Patch to a sphere - Titled” in pg. 10 in [Mar14], and B-43 in [HSM10]
• Disk: “Parallel configurations - Patch to disc” in pg. 22 in [Mar14], and B-13 in [HSM10]
• Rectangular: “Parallel configurations - Patch to rectangular plate” in pg. 22 and “Perpendicular -configurations - Patch to rectangular plate” in pg. 26 in [Mar14], and B-5 in [HSM10]

The solution in [HSM10] covers more general cases than [Mar14]. The overall runtime cost if we’d implement them naively is not so acceptable for real products, but they bring us the perfect analytic forms, we could use them wisely in products with some optimizations and simplifications, or just leave them as a ground truth faster than the classic path tracing.

A code example for sphere light that modified from [LR14]:

	for (int i = 0; i < NR_SPHERE_LIGHTS; ++i)
	{
		float lightRadius = sphereLightCBuffer.data[i].luminousFlux.w;
		if (lightRadius > 0)
		{
			vec3 unormalizedL = sphereLightCBuffer.data[i].position.xyz - posWS;
			vec3 L = normalize(unormalizedL);
			vec3 H = normalize(V + L);

			float LdotH = max(dot(L, H), 0.0);
			float NdotH = max(dot(N, H), 0.0);
			float NdotL = max(dot(N, L), 0.0);

			float sqrDist = dot(unormalizedL, unormalizedL);

			float Beta = acos(NdotL);
			float H2 = sqrt(sqrDist);
			float h = H2 / lightRadius;
			float x = sqrt(max(h * h - 1, eps));
			float y = -x * (1 / tan(Beta));
			float illuminance = 0;

			if (h * cos(Beta) > 1)
			{
				illuminance = cos(Beta) / (h * h);
			}
			else
			{
				illuminance = (1 / max(PI * h * h, eps))
					* (cos(Beta) * acos(y) - x * sin(Beta) * sqrt(max(1 - y * y, eps)))
					+ (1 / PI) * atan((sin(Beta) * sqrt(max(1 - y * y, eps)) / x));
			}
		}
	}

As you could see here, the form factor approach is suitable for our low-frequency signals - the diffuse part of the BSDF, since it is modeled basically around the phenomenon of energy bounce between diffuse surfaces. But still we need to find a solution for the specular part, and if you want to use some more advanced BTDF than the Lambert model which involves more about the complex subsurface scattering phenomenon, then the form factor solution is not accurate enough.

Representative Point

Another approach is that we won’t solve any of the integral at all, instead, we would try to substitute the problem with our old method, use punctual light to approximate area light. If we still focused on the lighting result, area light is similar to a bunch of point lights. And then we could find a point on the surface of the light source that contributes the most to the lighting result, and use it as the location of a point light to continue our shading process. This method is often called Representative Point or Most Representative Point (abbr. MRP), which has been successfully used in lots of PBR pipelines like [Kar13] and [LR14].

“…These approaches resemble the idea of importance sampling in Monte Carlo integral, where we numerically compute the value of a definite integral by averaging samples over the integral domain. In order

to do so more efficiently, we can try to prioritize samples that have a large contribution to the overall average…” - pg. 385, ”Real-Time Rendering”, 4th Edition.

When using the MRP method, we could both apply it to diffuse and specular part of BSDF, while it depends on the actual requirement of the products. One of the general algorithms has been developed by [Dro14], it uses the halfway vector method to find MRP for diffuse and cone projection to find MRP for specular. For the diffuse part, the halfway vector could be expressed as $\vec{h} = \frac{\vec{pp_0} + \vec{pp_1}}{||\vec{pp_0} + \vec{pp_1}||}$, $p$ is the shading point, $p_0$ is the intersection point of a ray from $p$ along the reflection vector $\vec{r}$ of the view direction, and $p_1$ is the intersection point of a ray from $p$ along negative direction of the light plane normal $\vec{n'}$. While the demonstration here may sound not complex, in a real scenario with an actual light shape, we often need to adjust the halfway vector method’s MRP result onto the surface of the light source, since it won’t always fall inside the light area. In another word, we need to find the closest point of the diffuse MRP in the light area if it’s outside.

For the specular part, the calculation became more complicated with the respect of the probability distribution functions (abbr. PDF) in the BRDF, we need to generate a cone along the reflection direction based on the PDF and use the geometry center of the intersection area between the cone and the light surface as the MRP. Furthermore, we could pre-compute lookup tables for a specific BSDF and light texture combination we chose, by eliminating some variables we could limit the LUT to 3D tables and free us from runtime calculation. And after we use the MRP to calculate the illuminance we should weight it by the intersection area, to preserve the energy conservation. All the details were well demonstrated in [Dro14], I’d recommend reading if it’s available for you.

Another way to calculate specular MRP is that we could find the point on the light surface with the smallest distance to the reflection ray of view direction, which is demonstrated in [Kar13]. For example, in a sphere light case we could calculate $p_{cr} = (l · r)r − l$ then $p_{cs} = l + p_{cr} · min(1,\frac{radius}{||p_{cr}||})$, while $p_{cr}$ is the point on the reflection ray closest to the sphere center, and $p_{cs}$ is the point on the sphere surface closest to $p_{cr}$.

Linearly Transformed Cosines

The LTC approach tries to solve the integral from another point of view, or more precisely speaking, another space. The reflectance equation involves several variables, and we’d evaluate each direction of the incoming light with a BSDF function, and the complexity is $O(n^2)$ due to the variance of the $dl$ and the BSDF. But if we could control one of them to be some sort of constant, then we would optimize the complexity to $O(n)$. Thanks to the popular usage of parametric BSDF nowadays, we could pre-compute lots of them into a multi-dimensional lookup table and just sample from it in runtime to save our precious frame time. But unfortunately, when we implement such a solution and try to fit it into the rendering pipeline, we’d still need to integrate the illuminance over a semispherical region, with an uncertain shape and orientations of the light surface. But if we consider the solving procedure of the reflectance equation in a linear algebra point of view, then the integral of the incoming light in a “world” spherical space could be solved as the integral of the transformed incoming light in a “local” spherical space. If we could choose a suitable transformation to ensure the linearity then the above hypothesis would become a provable truth. And if this transformation could eliminate the shape factor or the orientation factor of the light then we would be one step closer to our expectation of $O(n)$. Luckily the BSDF is naturally a spherical distribution function, the only problem left is how to find the transformation for a BSDF and find a close form of the integral in that BSDF space, after all, if we can’t solve the reflectance equation easier then all of these works would be nonsense.

The most commonly used distribution function of BSDF today is the GGX approximation for the Cook-Torrance microfacet model, the isotropic version of it only depends on the view direction and the surface roughness $\alpha$, even we add the anisotropicity later on, the number of the variables is still only 3. But even 3 variables the integral is still too expensive to calculate in runtime, we have to continue eliminating variables and transforming the light source until we finally find a space cheaper enough to evaluate the integral in real-time. That’s how the LTC came from, we would use a clamped cosine spherical distribution function as the base space since the integral over it is both cheap and analytic, then find a transformation matrix to approximately transform it to a specific BSDF, then use the inverse of that matrix to transform the vertices of the light source into the base space, and finally calculate the integral by utilising the analytic form of line integral.

The theoretical details are well demonstrated in [Hei16] and the reference C++ source code has been provided on Github, but without too much explanation about what it is doing. Actually, if you have ever implemented any BSDF LUT, the approach behind is similar, we just generate the value pair of all possible parameter combination ahead of runtime, and store them in an easy-to-fetch data structure. What we want to get here is a table of transformation matrices that, we could transform our specific BSDF distribution function to the clamped cosine one bidirectionally. Or mathematically speaking, we need to find the $M$ in $D_{BSDF}(\omega) = M * D_o(\omega)$, where $D_o$ would be the original clamped cosine distribution function we chose, with a form of $D_o(\omega_o = (x, y, z)) = \frac{1}{\pi}max(0, z)$ that $\omega_o$ is the normalized direction vector. As you could see we would employ a 3D vector as usual, so the matrix should be a 3x3 matrix.

The next step would be how to calculate the $M$ for our specific BSDF, for example, what the [Hei16] did for GGX distribution. They use Nelder–Mead method to numerically search the $M$, with a cost/error function implemented by simply comparing the actual BSDF value and LTC fitted value. And because in the paper they only fitted isotropic version of GGX, finally the $M$ only has 5 effective elements on each diagonal.

And with normalization by one of the elements there would be only 4 elements left, which could just be stored into a 4 channels texture and linear sampled in runtime. In the paper, they use $e$ as the normalizer, but later in [Hill16], they admit that it would bring some serious numerical error, and somehow I found that they use $c$ later again in the repo’s history, that should be a measured decision maybe. Also in [Hill16], they mentioned multiple problems when implementing LTC solutions, from how to deal with polygon clipping effectively to how to ensure the robust of inverse trigonometric functions, highly recommend to read.

After you get the LUT for $M$, then could just use it in runtime for the shading. You just need to write some shader codes to transform the light vertices to LTC space and do a spherical line integral by $\frac{1}{2\pi}\sum_{i=1}^{n}acos(v_i · v_j)\frac{v_i \times v_j}{||v_i \times v_j||} · n$. But this numerical integral is not practical for light shapes such as single line, capsule, sphere and disk, since you can’t find a small but smooth enough vertices number for ellipse shapes. Moreover, it’s the natural solution for convex polygons like rectangular light. Luckily later they provided the corresponding integral form for other shapes in [Hei17], where the final appearance and performance are quite plausible for real products in contrast with other solutions. I’d recommend again to find details in their publications.

To be continued.

Bibliography：

[Wiki1] https://en.wikipedia.org/wiki/View_factor

[Mar14] http://webserver.dmt.upm.es/~isidoro/tc3/Radiation%20View%20factors.pdf

[HSM10] http://www.thermalradiation.net/tablecon.html

[Kar13] B. Karis. “Real Shading in Unreal Engine 4”. In: Physically Based Shading in Theory and Practice, ACM SIGGRAPH 2013 Courses. SIGGRAPH ’13. Anaheim, California: ACM, 2013, 22:1–22:8. isbn: 978-1-4503-2339-0. doi: 10.1145/2504435.2504457. url: http://selfshadow.com/publications/s2013-shading-course/.

[LR14] S. Lagarde and C. de Rousiers. “Moving Frostbite to PBR”. In: Physically Based Shading in Theory and Practice, ACM SIGGRAPH 2014 Courses. SIGGRAPH ’14. Vancouver, Canada: ACM, 2014, 23:1–23:8. isbn: 978-1-4503-2962-0. doi: 10.1145/2614028.2615431. url: http://www.frostbite.com/2014/11/moving-frostbite-to-pbr/.

[Dro14] M. Drobot. “Physically Based Area Lights”. In: GPU Pro 5 Advanced Rendering Techniques. ISBN: 978-1-4822-0864-1. url: https://www.crcpress.com/GPU-Pro-5-Advanced-Rendering-Techniques/Engel/p/book/9781482208634

[Hei16] E.Heitz, J.Dupuy, S.Hill and D.Neubelt. “Real-time polygonal-light shading with linearly transformed cosines”. In: ACM Transactions on Graphics (TOG) Volume 35 Issue 4, July 2016 Article No. 41. url: https://eheitzresearch.wordpress.com/415-2/

[Hill16] S.Hill and E.Heitz. “Real-Time Area Lighting: a Journey from Research to Production”. In Advances in Real-Time Rendering in Games, ACM SIGGRAPH 2016 Courses. SIGGRAPH ’16. Anaheim, California: ACM, 2016, url: https://blog.selfshadow.com/publications/s2016-advances/

[Hei17] E.Heitz and S.Hill. “Real-Time Line- and Disk-Light Shading”. In Physically Based Shading in Theory and Practice, ACM SIGGRAPH 2017 Courses. SIGGRAPH ’17. Los Angeles, California: ACM, 2017, url: https://blog.selfshadow.com/publications/s2017-shading-course/heitz/s2017_pbs_ltc_lines_disks.pdf