文章

Introduction to Lighting Part 2

发表于
作者 褴褛飞旋
26 分钟阅读

1 Area Lights: Math Foundations

1.1 Why Are Area Lights Ideal?

Punctual Light因其简单性和较低的计算成本，至今仍广泛应用。然而，这类光源存在固有的概念缺陷，使其无法精确还原真实世界的光照效果。相较于面光源，其主要缺陷可以概括为两点：

  • Punctual Light的反射不具备物理真实性

    • 我们之所以能”看见”物体，是因为光线从被其他光源照亮的物体表面（间接光照）或直接从光源本身（直接光照）反射而来。场景中的物体具有体积，因此模拟其表面反射通常较为直观。然而，若光源本身没有体积尺寸，如何能在表面产生反射？

    • 如下图所示，我们可以观察到三个区域光在平面上的反射，其反射尺寸与场景中作为三维物体（本例为简单圆盘）的光源大小成正比。图中的最左侧为使用点光源的效果，虽然平面上可见该光源的反射，但反射尺寸与物理现实不符。因此，尽管存在反射现象，其物理真实感缺失且尺寸不可控，这在追求拟真图像时会造成显著挑战。

  • Punctual Light无法直接生成软阴影：区域光尺寸可以通过通过遮挡球体影响地面投影的半影区大小。光源越小，阴影边缘越锐利；反之光源越大，阴影过渡越柔和。由于点光源不具备物理尺寸，其设计原理决定了只能生成完全锐利的阴影，这在追求真实感的渲染中通常是不可取的。

1.2 A Quick, Friendly Intro to Concept of Monte Carlo Integration

让我们代入被区域光照明的点的视角。假设你就是下图中这个点，法线方向垂直向上，此时观察到的是一个球体部分遮挡了区域光表面。从这个视角不难推断：由于球体遮挡，仅有部分区域光发出的光线能到达你所在位置。

你的任务是量化”区域光表面发出的光线中有多少比例被球体遮挡？”虽然球体具有规整的几何特性（可能诱使你利用其对称性计算遮挡率），但如果换一个几何体呢？所以任务的核心是开发一种普适方法——无论遮挡物是平面、球体、兔子甚至导演临时起意要求的独角兽，都能精确计算光能损耗，也就是说：如何实现这种几何无关的遮挡计算？

设想一张512x512像素的纯白底图，其中心位置有一个半径为128像素的黑色圆盘。这足以模拟我们前面所讨论的情景：从着色点视角观察到的区域光，其中黑色圆盘代表遮挡光源的球体。为解决遮挡计算问题，我们将采用辐射度启发式方法：沿两个维度递归细分图像，首轮划分生成四个象限。持续细分直至达到最大层级，或当前象限完全位于球体内/外。若象限同时包含黑白像素则判定为与球体重叠。最终统计纯白象限内的像素总数（排除与球体重叠或完全被遮挡的象限），将其除以图像总像素数即可获得光源未被遮挡的比例。

所谓的辐射度启发式方法，是指将着色点上方的空间”细化”为越来越小的面片，通过识别独立的光照区块，将光能传输计算分解为面片间的能量传递。其设计原理导致难以有效处理镜面互反射，这成为其逐渐被光线追踪取代的主因——后者在足够算力支持下能卓越地模拟此类效果。

但，在着手统计之前，我们不难想到，通过数学方法，我们可精确计算区域光未被球体遮挡的比例。即，基于圆面积公式可以直接计算遮挡比例。这种解析计算法所得数值即为理论期望值，expected value。

如上图所示，随着细分次数增加，计算出的遮挡比例逐渐逼近80.37%，也就是理论期望值。这说明该技术有效，且通过增加四叉树细分层级可提升估算精度。值得注意的是，虽然我们以球体作为遮挡物进行演示，但该方法适用于任意形状物体。很好。

现在，我们回头思考一下当前的核心目标：确定着色点可见区域光的比例。从数学角度，这可以被形式化为积分问题——即通过对积分域（如区域光表面）进行求和/积分运算，获取可见性、光照面积或光强等总体度量。此前我们通过辐射度启发的空间细分法处理该问题，现在则将引入另一种经典方法：蒙特卡洛积分（Monte Carlo Integration）。

我在其他博客中有详细介绍蒙特卡洛积分的数学原理与实现过程，所以在这篇博客中，我们将从简单、非正式但实用的案例出发，逐步形式化核心概念。

蒙特卡洛积分法基于概率论，其核心思想是对区域光表面进行随机采样。若要计算某国成年人平均身高，测量全体人口显然不现实，但随机选取数千人测量后求均值即可获得统计意义上的合理近似。同理，我们可将此方法应用于区域光：想象向光源表面随机抛撒豆子（类比随机抽取人口样本），统计落在遮挡圆盘外（即从着色点可见区域）的豆子比例。数学实现非常简单：在512x512图像上生成随机点，判定各点是否位于球体外，统计可见点数量后除以总采样数。如下图所示：

以下为对应的Python代码实现：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49

import numpy as np
import matplotlib.pyplot as plt

def draw_disk(image, center, radius):
    for x in range(image.shape[0]):
        for y in range(image.shape[1]):
            if (x - center[0])**2 + (y - center[1])**2 <= radius**2:
                image[x, y] = 0  # Black for inside the disk

def is_inside_disk(center, radius, point):
    return (point[0] - center[0]) ** 2 + (point[1] - center[1]) ** 2 <= radius ** 2

# Initialize parameters
width, height = 512, 512
radius = 128
center = (width // 2, height // 2)
num_points = 128

# Create a white image
image = np.ones((width, height)) * 255

# Draw a black disk
draw_disk(image, center, radius)

# Monte Carlo simulation to estimate the area outside the disk
outside_count = 0
points = np.random.randint(0, high=width, size=(num_points, 2))

# Set up the plot
fig, ax = plt.subplots()
ax.imshow(image, cmap='gray', origin='upper')

# Check each point and adjust markersize
for point in points:
    inside = is_inside_disk(center, radius, point)
    if not inside:
        outside_count += 1
    color = 'ro' if inside else 'go'
    ax.plot(point[1], point[0], color, markersize=3)  # Smaller point size

# Calculate percentage of points outside the disk
percentage_outside = (outside_count / num_points) * 100

# Display the result
ax.set_title(f"Monte Carlo Estimation of Area Outside Disk: {percentage_outside:.2f}%")
ax.set_xticks([])
ax.set_yticks([])

plt.show()

如下图所示，我们可视化了通过蒙特卡洛积分计算这档比例的过程。注意在最后一张图像中（使用2048个样本），计算结果(81%)超过了理论期望值(80.37%)。蒙特卡洛方法的结果可能低于或高于期望值，这种估计值与理论值的偏差称为方差（Variance）。由于估算结果依赖于随机样本选择，即使保持相同样本数，更改随机数生成器的种子也可能导致完全不同的结果（可能高于、低于或等于期望值）。具体结果不可预测，因为这取决于随机过程——样本的选择。估算值是否偏离期望值并不关键，重要的是随着样本数量增加，结果接近理论值的概率会提升（若持续增加样本量可验证此现象）。样本量越大，我们对结果逼近理论值的置信度越高。

让我们总结一下本小节的主要内容：通过在积分域（光源表面）随机采样并评估目标函数（区域光可见性），最终将样本结果均值作为积分估计值。

1.3 Rendering Equation and Hemisphere vs Area Sampling

为深入探究区域光，我们需回溯计算机图形学中光照模拟的基本原理，由此引出渲染领域最重要的方程之一，渲染方程。本篇博客将仅讨论适用于仅漫反射表面的简化版本。该方程的原理并不复杂：假设你是一个着色点，下方为不透明表面，因此光照仅能来自上方空间。此时，以法线为轴的半球空间（如下图所示）囊括了所有可能入射方向的光线。在本示例场景中，着色点的光照来源于两个区域光。

我们需要对来自半球方向的光线进行”积分“。这构成积分问题的原因在于：半球是一个连续表面，而非离散元素集合。在计算机图形学中，当涉及方向集合（此处为半球）的积分时，需引入立体角Solid Angle概念——它类似于我们之前所讨论的微分面积，立体角可视为单位球面（此处为半球）上的面积表征，其单位为球面度，Steradian。单位球表面积为$4\pi$，而半球对应$2\pi$的球面度。

换种方式理解，可参考下图展示的从着色点视角出发的180度鱼眼视图。这本质上等同于该点”所见”场景，所有到达该点的光线均来自其”视野”范围。我们的任务是对该视图进行”扫描”并收集光线，但由于涉及连续空间，必须借助积分进行数学建模。

渲染方程正是为此而生。它对从半球空间内所有可能方向到达着色点P的光量进行积分，并在每个方向上乘以表面BRDF与Lambertian余弦定律。其数学表述如下：

\[L(x)=\frac{R(x)}{\pi}\int_{H^2}Li(x,\omega)cos(\theta)d\omega\]

其中：

  • $L(x)$表示着色点的radiance，也就是颜色值
  • $Li(x,\omega)$表示来自于方向$\omega$射向着色点$x$的入射光线
  • $\theta$为表面法线与光照方向之间的夹角
  • $H^2$表示积分域，在我们的例子中为整个半球，也可以用$\Omega$表示

渲染方程通过积分半球方向上入射到某点的光线量来求解。被积函数（即需要积分的部分）表示从半球内所有可能方向到达该点的光能。我们使用蒙特卡洛积分法进行求解，这与前文解决“可见光源表面比例”问题的思路一致。具体实现时，通过在半球表面随机选取若干方向，沿这些方向发射射线，并检测其是否与光源（如面光源）相交。将每条射线返回的颜色值（若未与光源相交则为0，相交则为光源的辐射亮度）累加，再根据采样总数取平均值。下图展示了地面某一点的采样过程：

观察上图可见，在投射的64条光线中，仅有2条与面光源表面相交，这种有效采样率极低的情况揭示了该方法的核心缺陷——噪点。为了准确计算面光源对当前着色点的光照贡献，我们不仅需要确保光线能检测到光源的存在，还需通过足够多的有效采样来充分覆盖光源表面。

更明确地说，当采样光线数量不足时，不同像素对应的有效击中面光源的光线数量会产生显著波动（例如某像素2次命中，相邻像素可能5次甚至0次）。这种采样结果的剧烈空间变化，会导致各像素对同一点光照强度的估算值差异极大，从而在图像中形成噪点。问题的根源不仅是估算方差较大，更在于相邻像素间的方差剧烈跳变。要缓解此问题，唯一途径是大幅增加采样数量。

虽然增加光线数量可缓解此问题，但会大幅延长渲染时间。因此，如何在保持采样数不变的前提下降低噪点？

解决方案在于逆向思考问题：不再对半球方向进行采样，而是直接对面光源表面进行采样。如下图所示，我们不再沿法线方向在半球空间生成随机方向，而是在面光源表面均匀生成采样点（该方法通过在光源表面均匀分布随机位置实现）。随后，将这些采样点与当前着色点连线，得到对应的光照方向向量。

只要稍作数学调整，这种方法是完全可行的。这里所说的数学调整是指在光源表面与其对应的立体角表示之间建立转换关系。可能的思路是将面光源投影到半球空间并在该区域采样，但这会引发一个复杂问题：多边形在球面上的投影会形成球面多边形（如下图所示），而此类几何体缺乏简洁的数学表征方法。

那我们应该如何解决？可将面光源上的每个采样点视为其表面的微分面积$dA$。我们需要将此微分面积$dA$转换为对应的微分立体角$d\omega$。幸运的是，两者间的转换关系存在且极为简洁：

\[d\omega=\frac{cos(\theta')}{|x-x'|^2}dA\]

其中：

  • $\omega$表示采样光线的方向，即$\vec{x’ - x}$

  • $\theta’$表示面光源法线与采样光线方向$\omega$之间的夹角
  • $x$与$x’$分别表示着色点与面光源表面上的采样点

如下图所示：

注意分母中的距离平方项体现了平方反比定律。该定律源于光线从光源辐射后呈球面扩散的特性。假设面光源向各方向均匀辐射（各向同性），那么光线覆盖的表面积会随距离光源的平方增长，导致单位面积能量按距离平方衰减。

那么，我们就给出新的形式的渲染方程：

\[\] \[L(x)=\frac{R(x)}{\pi}\int_{A}L_evis(x, x')\frac{cos(\theta')}{|x-x'|^2}dA\]

其中，$vis(x, x’)$这个项表示两点之间的可见性。注意此处积分域已从半球方向空间转换到了光源表面区域。技术实现上，若要将该数学表达式转化为代码，需再次讨论概率密度函数PDF的构建，但这部分将在稍后展开。我们可以给出伪代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13

uint32_t num_samples = 64;
Vec3f L = 0;
for (uint32_t n = 0; n < num_samples; ++n) 
{
    // Pick a random direction on the surface of the light
    Vec3f wi;
	float pdf;
	float t_max;
	Vec3f Le = light->Sample(x, r1, r2, &wi, &pdf, &t_max); 
	if (!Occluded(x, wi, t_max, scene)) // don't look for intersections beyond t_max
		L += 1 / M_PI * Le * max(0, wi.Dot(isect_info.N) * max(0, -wi.Dot(light->normal)) / (t_max * t_max * pdf)); 
}
L /= num_samples;

注意上述伪代码中，我们不再将光源视为具有自发光的几何体，而是从Light类派生的光源对象。因此，在大多数生产级渲染器中，面光源默认不会在图像中可见——它们不被当作几何体处理。代码中直接从光源获取其辐射能量（Le），并利用光源法线。t_max如前所述，存储了从着色点到光源采样点p_light的距离。光源表面采样点通过两个范围在[0:1]的随机数（伪代码中的r1和r2）计算生成，具体方法后续讨论。下图展示了此方法生成的图像效果。关于概率密度函数（pdf）的进一步解析将在后文展开。

1.4 Intro to PDF, Probability Distributions and Monte Carlo Estimators

我们将通过几个案例来直观理解PDF的定义及其在蒙特卡洛积分中的应用场景。以区域光采样为例：假设存在一个1×1单位的区域光，使用16个采样点时，我们会在光源表面均匀分布这些采样点，并计算每个采样点的光照方向。但如果现在需要采样一个2×2单位的区域光（面积扩大4倍），您是否意识到光源尺寸扩大后其贡献会显著增加？请注意，我们此时不应单纯将其视为大尺寸光源，而应理解为将4个原始光源紧密拼接组成超大光源——这意味着总能量变成了原始光源的4倍，对吗？

然而，若直接应用我们提供的区域光采样代码，会发现其中并未直接体现光源尺寸差异的影响（假设保持原始计算逻辑）。这意味着除非调整每次计算的概率密度函数（pdf）值，否则两种尺寸光源的采样结果将趋同。通过调整可实现：采样2×2大光源时，通过更低的pdf值（如原始1×1光源pdf=1，放大后pdf=0.25），使得最终贡献值通过”辐射值/0.25”操作获得真实四倍增量。这种数学修正完美对应物理规律，验证了路径的正确性。

[TODO]

Graphics, Scratchpixel
本文由作者按照 CC BY 4.0 进行授权