如果你调过景深散景光斑大概率见过那种“中心亮成一团、外圈稀稀拉拉”的奇怪画面。一开始我总怀疑是对焦逻辑算错了直到某次把采样点直接画出来才意识到问题根本不在光照而在一个看起来毫不起眼的环节——圆盘均匀采样。这也是我从PBRTPhysical Based Rendering Techniques里翻到对应章节时最有感触的地方渲染器里大量“看起来高级”的效果最后都依赖“在一个圆盘上均匀取点”这个基本功。这篇文章就围绕PBRT里的圆盘均匀采样展开把极坐标映射和同心映射两种方式的数学原理讲透给出单点和numpy批量版本的Python实现再用散点图、环带统计、网格变形实验验证它们的差异最后扩展到半球采样和景深应用。适合正在学蒙特卡洛渲染、想实现软阴影或景深、或者只是想把采样算法搞明白的人。1. 圆盘均匀采样到底在渲染里做了什么从散景到软阴影都在靠它1.1 三个最典型的应用场景圆盘采样不是一个孤立的数学玩具。蒙地卡洛渲染器对“随机”这件事有极强的依赖而圆盘往往是各种积分域里最基础的一个。第一个场景是景深。真实相机镜头是有限大小的圆光线从镜头光圈上不同的点进入。渲染器模拟薄透镜模型时会在光圈圆盘上采样一个点再从这点出发连接焦平面上的焦点。采样点分布直接决定散景的形状和亮度分布如果样本往圆心堆积散景光斑中心会比边缘亮如果样本全是外圈点散景会变成空心环。PBRT在处理相机采样时用的正是圆盘采样。第二个场景是软阴影。面光源要模拟出柔和过渡的阴影需要在光源表面采样多个点逐个去做可见性判断。很多实现会把面光源近似成圆盘或需要从圆盘映射过去的形状。光源采样点分布不均匀阴影区域就会出现条带和明显噪点。第三个场景是模糊后处理和环境光遮蔽。屏幕空间的景深模糊、泛光效果常常需要二维模糊核AO采样虽然通常在半球上做但半球余弦采样的基础也是圆盘均匀采样。正因为到处都要用PBRT才专门花篇幅推导圆盘采样这个“小函数”。1.2 “均匀”的真正含义以面积为尺度的概率密度很多人以为均匀采样就是“在圆盘上随机撒点看起来哪里都有”。但“均匀”在蒙地卡洛里的定义要严格得多圆盘上每个单位面积被采到的概率必须相同。换句话说概率密度函数 (p(x, y)) 在圆盘内是常数 (1 / \pi)因为圆盘面积为 (\pi)。这里有个关键逻辑我们手里的随机数发生器只能产生 ([0, 1)) 均匀分布的数值要把这些数值映射到圆盘上映射后的概率密度必须保持平坦。任何“看起来均匀但实际密度不平坦”的采样都会让积分估计产生偏差或额外方差。后面要说的极坐标映射和同心映射本质都是在解决“怎么让密度平坦”这件事。2. 从极坐标映射说起公式看着没毛病但缺陷同样明显2.1 为什么要用sqrt雅可比行列式的一次完整推导最直觉的做法是采样半径 (r) 和角度 (\theta)。令 (r u)(\theta 2\pi v)其中 (u, v) 是 ([0,1)) 的均匀随机数。这个做法对吗不对。问题出在面积度量上。圆盘上的面积元素是[ dA r \cdot dr \cdot d\theta ]半径越小同样一段 (dr) 对应的实际面积就越小。如果让 (r) 在数值上均匀分布那么半径0附近的一小段距离只对应很小的圆环面积但样本落在这一段的概率却不低结果就是圆心附近样本密度很高。正确的做法是令 (r \sqrt{u})。我们来验证一下。设[ x \sqrt{u} \cos(2\pi v), \quad y \sqrt{u} \sin(2\pi v) ]从 ((u, v)) 到 ((x, y)) 的雅可比行列式为[ J \begin{vmatrix} \frac{\partial x}{\partial u} \frac{\partial x}{\partial v} \ \frac{\partial y}{\partial u} \frac{\partial y}{\partial v} \end{vmatrix} ]代入后计算[ \frac{\partial x}{\partial u} \frac{\cos(2\pi v)}{2\sqrt{u}},\quad \frac{\partial y}{\partial u} \frac{\sin(2\pi v)}{2\sqrt{u}} ][ \frac{\partial x}{\partial v} -2\pi\sqrt{u}\sin(2\pi v),\quad \frac{\partial y}{\partial v} 2\pi\sqrt{u}\cos(2\pi v) ]行列式展开[ |J| \frac{\cos(2\pi v)}{2\sqrt{u}} \cdot 2\pi\sqrt{u}\cos(2\pi v)\left(-2\pi\sqrt{u}\sin(2\pi v)\right) \cdot \frac{\sin(2\pi v)}{2\sqrt{u}} ]化简后[ |J| \pi\cos^2(2\pi v) \pi\sin^2(2\pi v) \pi ]雅可比是常数 (\pi)。这意味着从 ((u,v)) 空间到圆盘概率密度函数保持不变((u,v)) 是均匀的映射后 ((x,y)) 在面积上也均匀。很多人直接把这段推导跳过然后就会对“为什么半径要用开方”百思不解。2.2 极坐标映射的两个软肋中心聚集感与网格扭曲既然雅可比是常数为什么实际效果还是有问题第一正方形左边 (u0) 这一整条边在极坐标映射里全部收缩到圆心一点 ((0,0))。从连续概率角度看这条边面积为零不影响均匀性但离散采样时如果样本在 (v) 方向上分布得比较密就会有很多点“挤”在圆心附近视觉上产生中心聚集感。第二局部形变特别大。(u) 方向上很小的变化在圆心附近对应 (dr du / (2\sqrt{u}))当 (u) 接近0时这个增量会非常大而 (v) 方向始终对应完整的 (2\pi)。正方形里一个个规整的小格子映射到圆盘上变成圆心附近细长的楔形、外圈又被压扁的条状。这种“各向异性扭曲”在纯随机采样时不那么显眼一旦换成低差异序列或分层采样样本之间的相邻关系就会被严重破坏低差异序列该有的均匀性也打了折扣。2.3 PBRT为什么偏要另起炉灶PBRT在处理这类问题时往往不只是问“概率密度是否均匀”还会问“映射是否保留了几何结构”。极坐标映射在概率密度上是达标的但它在把正方形变成圆盘的过程中将正方形内部大量的网格线压扁、拉长、甚至把一整条边收成一个点。PBRT认为这种映射的“扭曲”太大因此推荐同心映射concentric mapping同样是把正方形映射到圆盘但让正方形里的“同心正方形”轮廓一一对应到圆盘上的“同心圆”。这样既没有“一条边收缩到圆心”的奇怪行为网格线的变形也小得多和分层采样、低差异序列配合时表现更好。3. 同心映射原理拆解把正方形“卷”成圆盘3.1 从切比雪夫半径到极角几何直觉同心映射的思路是先看正方形里的点离中心有多“远”。这里的“远”指的是切比雪夫距离[ r \max(|u|, |v|) ]注意单位正方形先被映射到 ([-1, 1] \times [-1, 1])。切比雪夫距离相等的点组成一个个正方形轮廓比如 (|u| 0.5, |v| \le 0.5) 就组成一个正方形边框。同心映射把这些正方形轮廓按比例放到半径为 (r) 的圆上正方形轮廓上的点则按它们所在的“边”线性映射到对应的圆弧。正方形有四条边正好对应圆盘的四段圆弧。于是整个映射被自然地分成四个三角区域左三角、右三角、上三角、下三角。每个三角区域里的点都映射到圆盘对应的90度扇区内。3.2 分区域映射的雅可比验证每一片都是等面积同心映射的代码看起来像是“如果 (|u| |v|) 就怎么怎么样”很多人以为这只是个分支技巧。实际上这个分支背后是严格的面积保持。以右三角区域为例此时 (u |v|)映射写为[ r u, \quad \theta \frac{\pi}{4} \cdot \frac{v}{u} ]于是[ x u \cos\left(\frac{\pi v}{4u}\right),\quad y u \sin\left(\frac{\pi v}{4u}\right) ]对这个变换求雅可比行列式结果正好是常数 (\pi / 4)。也就是说右三角区域到右扇区的映射是等面积的。其他三个三角区域完全对称也都是常数雅可比。四个区域合起来覆盖整个圆盘边界之间严丝合缝因此整个圆盘上的概率密度处处平坦。这就是同心映射比极坐标映射更强的地方它不仅在概率密度上均匀还在局部面积比例上均匀网格结构的变形也小得多。PBRT选它作为默认方案是有充分理由的。3.3 代码行为细节带符号半径与局部角偏移同心映射在代码实现里比较绕半径 (r) 是带符号的角度 (\theta) 也不直接是最终极角。以 (|u| |v|) 为例代码令 (r u)。如果 (u) 为负半径就是负的(\theta) 只在 (-\pi/4) 到 (\pi/4) 之间变化。最终笛卡尔坐标 (x r\cos(\theta), y r\sin(\theta))。当 (r 0) 时颜色坐标其实落在左半圆实际极角是 (\theta \pi)。这是PBRT源码里的精简写法避免了显式判断左右半圆但理解时最好还是回到“四个三角区域分别映射到四个扇区”的几何直觉上。import math import random def concentric_sample_disk(rngrandom): # 将 [0,1)^2 映射到 [-1,1)^2 u rng.random() * 2.0 - 1.0 v rng.random() * 2.0 - 1.0 # 原点是一个特殊点直接返回 if u 0.0 and v 0.0: return [0.0, 0.0] if abs(u) abs(v): r u theta (math.pi / 4.0) * (v / u) else: r v theta (math.pi / 2.0) - (math.pi / 4.0) * (u / v) return [r * math.cos(theta), r * math.sin(theta)]注意这里u 0.0 and v 0.0的判断。如果不处理(|u| |v|) 为假会走 else 分支(r v 0)(\theta \pi/2)最终结果也是原点看起来没问题。但浮点数除法里出现 (0/0) 这种不确定情况时不同平台可能有不同表现所以PBRT显式把原点单独挑出来。4. Python实现单点版、numpy批量版与随机源注入4.1 单点采样函数逐行解释上面的代码就是单点版。它可以接受任意提供random()方法的对象默认用 Python 的random模块方便在测试时固定种子。import math import random def concentric_sample_disk(rngrandom): u rng.random() * 2.0 - 1.0 v rng.random() * 2.0 - 1.0 if u 0.0 and v 0.0: return [0.0, 0.0] if abs(u) abs(v): r u theta (math.pi / 4.0) * (v / u) else: r v theta (math.pi / 2.0) - (math.pi / 4.0) * (u / v) return [r * math.cos(theta), r * math.sin(theta)]这个版本的优点是简单、直观适合学习和小规模测试。缺点也很明显每次调用都要执行 Python 函数调用、多个数学函数生成百万级样本时速度不够快。4.2 numpy批量版向量化与除零的工程技巧实际用 Python 做渲染实验或数值验证时建议直接写向量化版本。核心思路是把u和v变成数组用np.where做条件分支。import numpy as np def concentric_sample_disk_np(u, v): # u, v 为 [0,1) 的数组 uu 2.0 * u - 1.0 vv 2.0 * v - 1.0 mask np.abs(uu) np.abs(vv) # 半径带符号 r np.where(mask, uu, vv) # 防止除零mask 为 True 时 uu 一定非零mask 为 False 时 vv 一定非零 # 用 1 当占位符不会影响被选中分支的计算 theta np.where( mask, np.pi / 4.0 * (vv / np.where(mask, uu, 1.0)), np.pi / 2.0 - np.pi / 4.0 * (uu / np.where(mask, 1.0, vv)) ) x r * np.cos(theta) y r * np.sin(theta) return x, y这里的核心技巧在分母处理np.where(mask, uu, 1.0)表示当 mask 为真时取uu做分母否则用 1 占位。因为 mask 为真时必然有 (|uu| |vv|)所以uu不可能为零同理 mask 为假时vv不可能为零。这样既不产生除零警告也不影响结果。用 numpy 生成一百万个样本几毫秒就能完成比纯 Python 单点版快两个数量级以上。做环带统计、散点图验证时这个版本是主力。4.3 随机源设计可复现与线程安全函数签名里的rngrandom不是可有可无的。渲染器里随机数来源通常不是 Python 的random而是自己的生成器固定种子的伪随机数、Halton/Sobol 低差异序列、蓝噪声序列等。把随机源作为参数传入测试时方便复现接入其他序列时也不用改函数内部逻辑。import numpy as np from scipy.stats import qmc # 生成低差异序列代替 np.random.rand sampler qmc.Sobol(d2, scrambleTrue) samples sampler.random(10000) u_arr samples[:, 0] v_arr samples[:, 1] x_arr, y_arr concentric_sample_disk_np(u_arr, v_arr)这种设计也符合 PBRT 里“采样器sampler与映射函数分离”的思路采样器负责产出不同分布的原始二维坐标映射函数负责把坐标变换到目标域。5. 用可视化与统计检验说话均匀性不是“看起来均匀”5.1 散点图与环带计数先直观再定量光看散点图容易误判。比如极坐标映射在样本数不够多时视觉上也会有“貌似均匀”的效果。我的经验是先做散点图再做环带计数用数据说话。环带计数的办法很简单把圆盘按面积均分成 (n) 个环带第 (k) 个环带的外半径是[ r_k \sqrt{\frac{k}{n}} ]因为圆盘面积占比与半径平方成正比。统计每个环带里的样本数理论上应该接近 (N / n)。def ring_counts(x_arr, y_arr, n_rings10): r2 x_arr**2 y_arr**2 # 半径平方 ring_idx np.clip((r2 * n_rings).astype(int), 0, n_rings - 1) counts np.bincount(ring_idx, minlengthn_rings) return counts # 生成 200000 个样本做测试 rng np.random.default_rng(42) u rng.random(200000) v rng.random(200000) xc, yc concentric_sample_disk_np(u, v) # 极坐标映射做同样样本数对比 xr np.sqrt(u) * np.cos(2 * np.pi * v) yr np.sqrt(u) * np.sin(2 * np.pi * v) print(同心映射环带计数:, ring_counts(xc, yc)) print(极坐标映射环带计数:, ring_counts(xr, yr))实测中两种映射的环带计数都接近均匀差距在统计噪声范围内。这说明极坐标映射从概率密度角度看也是达标的环带计数不足以区分它们必须看下一节的网格变形实验。5.2 网格变形实验分层采样的底层差异要区分极坐标映射和同心映射最直观的方法是取一个规整网格把网格点当作二维样本送进映射画出映射后的连线图。极坐标映射下正方形左侧边缘整条线被吸到圆心网格单元在圆心附近被拉成极长的楔形越靠近边缘又被横向压缩。同心映射则不同网格轮廓保持成一个个近似圆的正方形同心轮廓每个局部网格单元虽然也有拉伸但远没有极坐标映射那么剧烈。因此当样本来自分层采样或低差异序列时同心映射能更好地保留样本之间的“间隔感”不会出现“同一行样本在圆心附近堆成一团”的问题。这也是 PBRT 优先选择同心映射最实在的理由。5.3 性能与代码可读性的取舍纯 Python 单点版适合理解原理numpy 批量版适合做实验。如果你在写一个真正的渲染器性能关键路径上其实很少用 Python 按点调用而是直接写 C/Rust 和 SIMD 版本或者把 numpy 批量生成的样本一次性传给光线追踪内核。我做这套验证时的体感是算法逻辑本身不到 20 行但在 Python 里做批量采样时分支处理才是最容易出 bug 的地方。尤其是np.where里分子分母同时为零的场景建议写一个小型断言提前把原点坐标过滤掉或单独处理避免出现 NaN 悄悄混进后续流程。6. 从二维圆盘到三维半球Malley方法与向下扩展6.1 Malley方法均匀圆盘天然生成余弦密度半球采样圆盘采样最有价值的扩展是把它变成半球上的余弦采样。这在漫反射表面的环境光照积分、AO 采样里非常常见。思路来自一个投影关系把圆盘上的均匀采样点 ((x, y))“垂直抬升”到半球面上取[ z \sqrt{1 - x^2 - y^2} ]这样就得到一个半球上的点 ((x, y, z))。这个看似简单的提升得到的样本方向分布恰好是余弦密度[ p(\omega) \frac{\cos\theta}{\pi} ]也就是越靠近法线方向(z) 轴的样本越多符合朗伯散射的特性。PBRT 里的CosineSampleHemisphere就是基于圆盘均匀采样实现的方法名叫 Malleys method。def cosine_sample_hemisphere(u, v): x, y concentric_sample_disk_np(u, v) z np.sqrt(np.maximum(0.0, 1.0 - x**2 - y**2)) return x, y, z注意这里对1 - x**2 - y**2做了maximum(0, ...)保护。浮点误差可能让圆盘边缘的点算出一个微小的负数开根号后变成 NaN。6.2 均匀半球采样与余弦采样的差异均匀半球采样是另外一条路直接令 (z u)角度 (\phi 2\pi v)再算 (r \sqrt{1 - z^2})。这样得到的半球表面单位面积概率相等。但渲染漫反射时一般用余弦采样因为朗伯模型的反射亮度正比于 (\cos\theta)余弦采样能有效降低方差。对比代码的话均匀半球和余弦半球差的只是“先在平面圆盘上取样再提升”这一步。从圆盘采样出发的 Malley 方法在工程上还有一个好处它可以复用圆盘采样的所有改进包括分层、低差异序列和同心映射。6.3 景深应用的简化代码示例把圆盘采样接进景深逻辑也很直接。薄透镜模型里从胶片点 (p) 出发的射线先指向焦平面上的点再在镜头光圈上采样一个点以这两个点为端点确定光线方向。def sample_ray_with_dof(camera, u_disk, v_disk): px, py concentric_sample_disk_np(u_disk, v_disk) lens_point ( camera.lens_center px * camera.lens_u py * camera.lens_v ) # 焦平面焦点镜头中心沿方向向量走焦距 focus_point camera.lens_center camera.focal_distance * camera.forward ray_dir focus_point - lens_point ray_dir ray_dir / np.linalg.norm(ray_dir) return lens_point, ray_dir这里光圈采样点分布越均匀散景光斑的亮度过渡就越自然。如果你只在圆心附近采样光斑中心会明显过曝。7. 落地时的坑与心得边界、低差异序列和扩展方向7.1 浮点边界与零坐标处理圆盘采样代码里最容易翻车的是边界条件。u和v从随机数发生器出来理论上不会等于 0但低差异序列比如 Sobol 序列的第一个样本经常就是 ([0, 0])。送到映射函数里半径计算和角度计算的零除问题就会暴露。处理方式有两种一种是在映射函数里像 PBRT 一样显式判断原点并返回 ([0, 0])另一种是在生成样本时把 0 替换成极小正值比如u np.clip(u, 1e-10, 1.0 - 1e-10) v np.clip(v, 1e-10, 1.0 - 1e-10)用裁剪而不是直接过滤能保持样本数量不变适合批量生成场景。7.2 低差异序列替换随机数采样点会“均匀得更有纪律”如果只是做后处理模糊用np.random.rand生成圆盘样本也够用。但在渲染器里特别是需要低方差积分的时候低差异序列几乎是必需品。Sobol 序列生成的二维点在单位正方形里分布均匀经过同心映射到圆盘后分布依然均匀且不会出现普通随机数那种“这里一团、那里一片空地”的聚簇现象。用scipy.stats.qmc.Sobol时要注意序列第一个点通常是全 0直接送去计算可能出现上面说的原点问题另外低差异序列的样本数一般取 2 的幂如 1024、4096、65536才能发挥最佳排列效果。7.3 实际渲染器里的使用方式与进一步扩展PBRT 里的圆盘采样不是单独存在的它和整个采样基础设施连在一起像素网格分层采样、低差异序列、自适应采样计数器最终才把二维样本送到ConcentricSampleDisk。如果你后面要做球体采样、圆环采样、圆柱采样思路也是一样的先写出“单位正方形到目标域”的映射算清楚雅可比是否为常数再写代码。圆盘均匀采样是这套方法论里最简单也最经典的一个例子吃透它再看别的采样域会轻松很多。这段时间在极坐标映射和同心映射之间来回切换最大的教训是看分布图不能只看散点密度还要看样本序列的内部结构。圆盘均匀采样看似只是几十行代码的小函数但在景深、软阴影、半球余弦采样这些地方它决定了光斑亮度是否自然、阴影过渡是否干净。如果你只是临时生成一个模糊核极坐标映射凑合能用但要做低差异采样、分层采样或者完整的渲染器采样管线还是直接上同心映射更稳。代码也就多几行收益却相当实在。