企业官网建设流程全解析

做图形学渲染相关实验时，大概率都会撞上这样一个需求：要在单位圆盘上生成均匀分布的采样点。我第一次正经面对它，是在复现PBRT（Physically Based Rendering Toolkit）里的光照积分器时——不管是模拟景深镜头的薄透镜采样，还是给圆形面光源做重要性采样，圆盘均匀采样都是绕不开的底层工具。当时我图省事，直接用Python写了个"随机半径+随机角度"的版本，点画出来全挤在圆心附近，渲染结果整张图都在闪，噪点密集到我怀疑是光照公式写错了。后来把PBRT里那套同心映射方案搬过来用Python重写了一遍，配合可视化验证，才彻底搞明白这个"看似简单"的采样问题背后藏着多少细节。

这篇文章就来完整拆解圆盘均匀采样从数学原理到Python落地的全过程，把四种主流方案挨个对比，给出可直接运行的代码和验证方法，再把我在实际项目里踩过的坑一并交代清楚。适合正在学路径追踪、刚开始啃PBRT源码、或者做图形学作业需要在Python里补一个可靠采样器的同学参考。如果你只想抄代码，直接去第3节；如果你想弄明白为什么代码要这么写，建议从头看。

1. 圆盘均匀采样在PBRT里解决什么问题

1.1 采样器在渲染系统中的位置

现代渲染基本都是在算积分，而且是高维、无解析解的积分。路径追踪采用蒙特卡洛方法，把积分拆成无数个采样点，每个采样点代表一条光路，最终通过对这些点的亮度求平均来逼近真实结果。既然是蒙特卡洛，采样点是怎么生成的、分布是什么样的，就直接决定了收敛速度——采样分布越匹配被积函数，方差越小，画面越干净。

在PBRT这类物理渲染框架里，采样这件事被拆成两层：底层是采样器（Sampler），负责产生多维的低差异序列或伪随机数；上层是一系列映射函数，把[0,1)^2这个标准单位正方形上的均匀样本，转换到我们真正需要的几何域上，比如半球方向、球面坐标、三角形内部，以及今天的主题——圆盘上的坐标点。圆盘采样就是这些映射函数中很基础的一种：输入两个[0,1)之间的数u1, u2，输出一个落在圆盘内的(x, y)。

PBRT中所有采样器产出的多维样本，最终由积分器决定每个维度怎么用。有的维度用于挑选光源，有的用于BRDF方向采样，有的用于透镜位置采样。圆盘采样通常就是透镜位置或者圆形光源位置选择的最后一跳。所以它看起来不起眼，但因为被调用的频率极高，一个小缺陷会被放大成整张画面的噪声。

1.2 "均匀"的准确含义与常见误解

很多人以为均匀采样就是"每个点等概率被选到"，对连续分布来说，这句话不够精确。真正的均匀，是指概率密度函数p(x, y)在圆盘内任意位置的微小面积元dA上都是常数。对单位圆盘来说，p(x, y) = 1/π，因为圆盘面积是π，密度的面积分应该等于1。

问题在于"每单位面积密度恒定"不等于"每单位半径密度恒定"，也不等于"每单位角度密度恒定"。最典型的错误写法是：

theta = 2 * np.pi * u1 r = u2 x = r * np.cos(theta) y = r * np.sin(theta)

这段代码生成的点，半径在[0,1]上均匀，角度在[0,2π)上均匀。写成直方图，半径均匀确实没毛病，但画在平面上，点会明显聚集在圆心附近。原因是几何直觉欺骗了你：半径为r、宽度为dr的圆环，面积大约是2πr·dr，半径越小，同样dr对应的面积越小。如果半径本身均匀采样，落在小半径圆环内的概率密度就更高，最终得到的点密度就不是每单位面积恒定，而是中心密、边缘疏。

这也是为什么说"均匀"必须先明确是在哪个度量空间里的均匀。我们在圆盘上做积分，自然的度量是面积，所以必须保证面积均匀。

1.3 雅可比行列式：绕不开的数学账

从映射的角度看，圆盘采样是一个从参数域到几何域的映射D: (u1, u2) → (x, y)。要保证面积均匀，映射必须满足一个条件：雅可比行列式的绝对值在定义域内是常数。

简单回顾一下：雅可比行列式衡量的是"参数空间的一小块面积，映射到几何空间后变成了多大"。如果参数域上处处面积均匀，但映射后几何域上的面积元在某些地方被拉伸了、某些地方被压缩了，那几何域上的点密度自然就不均匀。要让输出均匀，就必须让雅可比行列式等于几何区域面积的倒数，也就是常数1/π。

用极坐标映射来算这笔账：x = r·cosθ, y = r·sinθ，如果直接令r = u2, θ = 2π·u1，雅可比行列式|det(J)| = r，这是变量r的函数，不是常数。所以输出不均匀，问题就出在这里。要让行列式变成常数，必须令r与u2之间满足平方根关系。具体推导如下：圆盘均匀分布下，半径R的边缘概率密度为f_R(r) = 2r/R²（因为半径为r的圆环面积占比与其周长成正比），其累积分布函数为F_R(r) = r²/R²。逆变换采样告诉我们，令F_R(r) = u，解出r = R·√u，这样再配一个独立的均匀角度θ = 2πv，输出就是均匀的。此时雅可比行列式恒为1/(2π)的常数，乘以角度范围后正好满足单位面积均匀。

这一节如果没看明白，后面就不用往下看了，因为所有正确实现都离不开这条平方根关系。

2. 四种采样方案对比：从毛坯到精装

2.1 拒绝采样法：简单粗暴，天生浪费

最直观的方法是拒绝采样：在[-1,1]²的外接正方形里生成点(x,y)，如果x²+y² > 1就丢掉重来。代码三行：

def rejection_sample_disk(rng): while True: x = 2.0 * rng.random() - 1.0 y = 2.0 * rng.random() - 1.0 if x * x + y * y <= 1.0: return x, y

从数学上讲，这个方法是完全正确的，均匀性天然成立，没有任何偏倚。接受率是圆面积除以正方形面积，即π/4 ≈ 78.5%，也就是说平均每生成100个点，有21个多会被扔掉。

但它在渲染管线里有三个致命缺点。

第一，浪费。路径追踪里每次采样都很贵，白白扔掉五分之一以上的随机数，等于渲染时间白白多出一截。

第二，破坏序列结构。PBRT这类框架如果用低差异序列（比如Sobol序列）来降低方差，拒绝采样会打乱序列的分布规律，低差异特性被严重破坏。低差异序列讲究的是"尽量均匀铺满整个空间"，拒绝掉一部分点，剩下的点在圆盘上的覆盖就不那么均匀了。

第三，不优雅。虽然这不是工程问题，但当你手头有一堆高质量的[0,1)^2序列样本时，你希望的是把它们整体映射成圆盘上的高质量样本，而不是丢三落四。

拒绝采样适合做快速验证，或者在非性能敏感、非低差异场景下临时顶一下。我自己的项目里，它只出现在早期写原型验证算法思路的时候，从来没有进过正式渲染管线。

2.2 极坐标直接采样：经典翻车现场

开头提到的直接采样法，也就是令θ = 2π·u1, r = u2的写法，属于典型的翻车操作。很多人刚开始都会这么写，因为它看起来很有道理：角度均匀、半径均匀，圆盘不就是由这两个量决定的吗？

翻了车才明白问题所在：(r, θ)的均匀和(x, y)的均匀不是一回事。从面积度量的角度，圆盘上以极坐标表示的面积元是dxdy = r·dr·dθ，这里多了一个r的权重。如果(r, θ)均匀采样，等于在面积元上额外乘了权重r，直观表现就是圆心附近被过度采样。

用蒙特卡洛估计来算会更清楚。假设我们要估计圆盘上某个函数f(x,y)的积分，蒙卡估计器写作(1/N)·Σ f(x_i, y_i)/p(x_i, y_i)。如果实际采样密度p在圆心附近比理论值大，而在边缘比理论值小，那么估计量和真实积分之间就会出现系统性偏差——不是方差大，而是直接偏了。渲染里表现成整体亮度偏高或偏低，而且怎么增加采样数都修不回来。

从雅可比行列式看：x = r·cosθ, y = r·sinθ，当r = u2时行列式是u2，在u2=0附近面积元趋近于0，映射极度收缩，样本挤成一团；在u2=1附近面积元膨胀，样本稀疏。这正是"中心聚集、边缘稀疏"的数学根源。

2.3 平方根补偿后的教科书方案

把上一节的教训用上，就得到了最经典的教科书实现：

def uniform_sample_disk_polar(u1, u2): theta = 2.0 * np.pi * u1 r = np.sqrt(u2) x = r * np.cos(theta) y = r * np.sin(theta) return x, y

注意这里只有一处变化：r = np.sqrt(u2)。这个平方根完美补偿了极坐标面积元中的r因子。此时(x,y)在圆盘上是面积均匀的，p(x,y) = 1/π。

这个方案在数学上完全正确，实现也足够简单，很多入门教程都用它。但它有一个不那么容易察觉的短板：映射在圆心附近存在奇异性。参数域[0,1)^2中接近u2=0的一小片矩形区域，会被映射到圆心附近一个极小的扇形区域内，这里的局部拉伸非常剧烈。如果用散点图观察低差异序列产生的样本，会看到圆心附近出现放射状的条纹，某些方向上点明显更密。

对于纯伪随机数来说，这个短板不明显；对于追求极致均匀性的低差异序列来说，这种局部扭曲会在渲染结果里表现为奇怪的低频花斑。所以PBRT没有止步于此，而是选用了下一节的方案。

2.4 PBRT采用的同心映射方案

PBRT实际使用的是同心映射（Concentric Disk Mapping）。思路很巧妙：与其在极坐标的"半径—角度"体系里打转，不如把单位正方形按照对角线划分成四个扇区，然后将每个扇区平滑地映射到圆盘的对应区域。这样避免了在圆心处出现一个"参数域小方块被拉伸成扇形"的奇点，整个映射在圆盘内部处处连续、扭曲更均匀。

实现逻辑如下：

import numpy as np def concentric_sample_disk(u1, u2): sx = 2.0 * u1 - 1.0 sy = 2.0 * u2 - 1.0 if sx == 0.0 and sy == 0.0: return 0.0, 0.0 if sx >= -sy: if sx > sy: r = sx theta = np.pi * 0.25 * (sy / sx) else: r = sy theta = np.pi * 0.25 * (2.0 - sx / sy) else: if sx < sy: r = -sx theta = np.pi * 0.25 * (4.0 + sy / sx) else: r = -sy theta = np.pi * 0.25 * (6.0 - sx / sy) return r * np.cos(theta), r * np.sin(theta)

这个代码读起来稍绕，但拆开看就清楚了：(sx, sy)是把[0,1)^2平移到[-1,1)^2后的坐标。四个分支对应四条对角线划分出的四个扇区，每个扇区都分别计算半径r和角度theta。角度被限制在每一份π/2范围里，半径r取值在[0,1]之间，但不同扇区的半径来源不同——有的来自|sx|，有的来自|sy|，目的是让映射在扇区边界处光滑衔接。

为什么这个映射比极坐标方案好？核心在于雅可比行列式的变化模式。极坐标方案在圆心处行列式趋近于零，局部扭曲极大；同心映射的行列式在整个圆盘内只是平滑变化，没有奇点，所以从参数域到圆盘的映射更"温和"。PBRT的代码设计哲学是：宁可把扭曲分散到整个域上，也不要集中在一处爆炸。

四种方案对比如下表：

3. Python实现与验证代码逐段解析

3.1 完整代码与接口设计

我一般在Python项目里封装一个采样器模块，对外暴露统一接口，输入是两个[0,1)的样本值，输出是圆盘上的坐标。这样做的好处是，上游无论是用numpy.random还是低差异序列，接口都不用变。

import numpy as np class DiskSampler: def __init__(self, radius=1.0): self.radius = radius def sample(self, u1, u2): x, y = concentric_sample_disk(u1, u2) return x * self.radius, y * self.radius def pdf(self): return 1.0 / (np.pi * self.radius * self.radius)

这里我额外加了一个pdf()方法，返回圆盘采样的概率密度。路径追踪里做重要性采样或者计算蒙特卡洛权重时，这个密度值是必不可少的。很多初学者写采样器只写sample不写pdf，结果算权重时到处找公式，实际上这里就是最简单的常数。

sample方法输入的u1, u2可以是标量，也可以是任意形状的numpy数组。为了支持这一点，sample内部用的是上一节的concentric_sample_disk，但我们会把它改造成支持向量运算的版本，见下一节。

3.2 向量化版本与边界条件

逐点调用Python函数生成十万个样本，性能不可接受。真正写进渲染循环里的采样器，必须支持numpy数组批量计算。把同心映射改造成向量化版本，关键是去掉if-else分支，改用掩码筛选和np.where逐层选择：

def concentric_sample_disk_vec(u1, u2): u1 = np.asarray(u1, dtype=np.float64) u2 = np.asarray(u2, dtype=np.float64) sx = 2.0 * u1 - 1.0 sy = 2.0 * u2 - 1.0 # 四个扇区掩码 c1 = (sx > 0.0) & (sx > sy) & (-sy < sx) # 右下扇区 c2 = (sy > 0.0) & (sx <= sy) & (-sy <= sx) # 右上扇区 c3 = (sx < 0.0) & (sx < sy) & (-sy > sx) # 左上扇区 c4 = (sy < 0.0) & (sx >= sy) & (-sy <= sx) # 左下扇区 # 这里还需要处理主对角线边界。为清晰起见，直接用原始条件。

说句实话，上面这种手写掩码的方式很容易把扇区边界搞错，我实际项目里踩过这个坑。更稳妥的做法是保持原始的四分支结构，用np.where按优先级逐层选择，再配合np.errstate抑制除零警告：

def concentric_sample_disk_vec(u1, u2): u1 = np.asarray(u1, dtype=np.float64) u2 = np.asarray(u2, dtype=np.float64) sx = 2.0 * u1 - 1.0 sy = 2.0 * u2 - 1.0 r = np.zeros_like(sx) theta = np.zeros_like(sx) with np.errstate(divide='ignore', invalid='ignore'): q1 = (sx >= -sy) & (sx > sy) r = np.where(q1, sx, r) theta = np.where(q1, np.pi * 0.25 * (sy / sx), theta) q2 = (sx >= -sy) & (sx <= sy) r = np.where(q2, sy, r) theta = np.where(q2, np.pi * 0.25 * (2.0 - sx / sy), theta) q3 = (sx < -sy) & (sx < sy) r = np.where(q3, -sx, r) theta = np.where(q3, np.pi * 0.25 * (4.0 + sy / sx), theta) q4 = (sx < -sy) & (sx >= sy) r = np.where(q4, -sy, r) theta = np.where(q4, np.pi * 0.25 * (6.0 - sx / sy), theta) center = (sx == 0.0) & (sy == 0.0) r = np.where(center, 0.0, r) theta = np.where(center, 0.0, theta) return r * np.cos(theta), r * np.sin(theta)

这里np.where会先完整计算所有分支的表达式，包括那些会得到nan或inf的项，但np.errstate把警告压掉了，最后np.where只挑选对应扇区的有效值。中心点单独处理，因为(0,0)处无论哪个扇区公式都会除零。

这个向量化版本经过测试，和逐点循环版结果完全一致，性能有几十倍的提升。如果要在GPU上跑，也可以直接把相同的逻辑翻译成TensorFlow或PyTorch的张量运算。

3.3 可视化检查：散点、径向与角向分布

写完采样器，第一件事是画图。用一万个样本，分别做散点图、半径平方直方图和角度直方图：

import matplotlib.pyplot as plt rng = np.random.default_rng(42) u1 = rng.random(10000) u2 = rng.random(10000) x, y = concentric_sample_disk_vec(u1, u2) fig, axes = plt.subplots(1, 3, figsize=(15, 4)) axes[0].scatter(x, y, s=1, alpha=0.5) axes[0].set_aspect('equal') axes[0].set_title('Scatter') r2 = x * x + y * y axes[1].hist(r2, bins=50) axes[1].set_title('r^2 histogram') theta = np.arctan2(y, x) axes[2].hist(theta, bins=48) axes[2].set_title('angle histogram')

散点图应该均匀铺满单位圆，没有明显聚集或空洞；半径平方的直方图应该近似水平；角度的直方图应该近似水平。这三个检查跑一遍，基本就能判断采样器有没有写对。

4. 采样质量的三层检验手段

4.1 统计直方图与KS检验

直方图只能给个直观印象，要量化判断，得用统计检验。前面提到，单位圆盘均匀采样的关键判据是r²在[0,1]上均匀分布，角度在[0,2π)上均匀分布。所以我们可以对r²做KS检验（Kolmogorov-Smirnov test）：

from scipy import stats r2 = x * x + y * y ks = stats.kstest(r2, 'uniform', args=(0.0, 1.0)) print(f'KS statistic = {ks.statistic:.4f}, p-value = {ks.pvalue:.4f}')

p-value大于0.05时，不能拒绝"来自均匀分布"的原假设，说明采样器的径向分布没有明显问题。注意KS检验不是万能的，它对"局部聚集但不影响整体分布"的模式不敏感，所以还需要配合区域计数检验。

角度同样可以检验，把角度值归一化到[0,1]后做KS检验。不过我实际经验是，角度分布一般不会出问题，最容易出问题的是半径关系，所以先紧着r²检验。

4.2 卡方检验与区域计数

卡方检验是另一个更细的检查方式。把单位圆盘划分成若干个同心环和扇区组成的格子，统计每个格子里的样本数量，和理论期望数量做卡方比对：

def chi_square_disk_test(x, y, n_rings=5, n_sectors=12, seed=42): r = np.sqrt(x * x + y * y) theta = np.arctan2(y, x) % (2 * np.pi) ring_idx = np.minimum((r * n_rings).astype(int), n_rings - 1) sector_idx = np.minimum((theta / (2 * np.pi) * n_sectors).astype(int), n_sectors - 1) observed = np.zeros((n_rings, n_sectors), dtype=int) for i, j in zip(ring_idx, sector_idx): observed[i, j] += 1 # 每个格子的面积占比决定期望计数 expected = np.zeros_like(observed, dtype=float) for i in range(n_rings): area_frac = ((i + 1) ** 2 - i ** 2) / (n_rings ** 2) expected[i, :] = area_frac * len(x) / n_sectors chi2 = np.sum((observed - expected) ** 2 / expected) dof = n_rings * n_sectors - 1 p_value = 1.0 - stats.chi2.cdf(chi2, dof) return chi2, p_value

卡方检验比KS检验更细致，因为它要求每个空间格子的计数都和理论值接近，能捕捉到"整体均匀但局部聚集"的问题。实测下来，正确的采样器两种检验都能通过，而极坐标直接采样法会在r²的KS检验上直接挂掉，p-value小到浮点数精度都不够展示。

4.3 蒙特卡洛积分误差对比实验

统计检验通过，只能说明采样点分布符合均匀分布，但采样方法好不好，最终还要看它在实际蒙卡积分里的表现。我做一个简单实验：估算单位圆盘上函数f(x,y) = exp(-(x²+y²))的积分。这个积分的真值是π·(1 - e⁻¹) ≈ 1.9859。

true_value = np.pi * (1.0 - np.exp(-1.0)) def mc_estimate(sample_fn, n_samples=10000, n_trials=200): errors = [] for _ in range(n_trials): u1 = rng.random(n_samples) u2 = rng.random(n_samples) x, y = sample_fn(u1, u2) estimate = np.mean(np.exp(-(x * x + y * y)) / (1.0 / np.pi)) errors.append((estimate - true_value) ** 2) return np.mean(errors)

分别传入极坐标直接采样、极坐标+开根号、同心映射三种采样函数，计算各自的均方误差。以我实测的结果看：极坐标直接采样法因为样本分布不均匀又强行按均匀密度加权，均方误差比正确方法大两个数量级，而且均值有系统性偏差；极坐标+开根号和同心映射在纯随机数上表现接近，均方误差基本就在理论方差附近；但换成Sobol低差异序列后，同心映射的误差会明显低于极坐标+开根号，这正体现了同心映射对低差异序列更友好这一特点。

这个实验很值得自己跑一遍，因为只有亲眼看到错误版本产生的偏差，才会对那条平方根关系有刻骨铭心的认识。

5. 实践场景、踩坑记录与优化建议

5.1 渲染管线里的典型用途

圆盘均匀采样在渲染管线里最常见的用途有三个。

第一个是薄透镜相机模型。真实相机的光圈是一个圆，聚焦时从同一物点出发的光线经过镜头后，会在成像面上汇聚；没有完美聚焦的情况下（即景深效果），光线会穿过光圈的不同位置，在像面上形成弥散圆。PBRT的薄透镜模型在计算景深时，需要在镜头光圈这个圆盘上采样一个位置，这个位置就是通过圆盘采样得到的。

第二个是圆形面光源采样。场景里如果放了一个圆形的面光源，光源上任一点对某个着色点的贡献是不同的（因为距离和角度不同），做光源重要性采样时需要在这个圆形区域上均匀采样一个点，然后把它作为阴影射线的起点。

第三个是余弦加权半球方向的采样，这也是我最常遇到的一处应用。兰伯特漫反射BRDF的采样方向分布正比于余弦值cosθ，而Malley方法提供了一个漂亮的构造：先在单位圆盘上均匀采样一个点(x, y)，然后把它垂直投影到半球面上，得到的方向自动就服从余弦加权分布。这个技巧直接复用了圆盘采样器，不需要额外写一个半球采样器。

5.2 我踩过的四个坑

第一个坑就是开头说的忘记开根号。这个坑的特点是症状非常隐蔽——图像不是全黑，而是比正确结果暗一些，噪点集中在中心亮部区域。我当时排查了很久，甚至怀疑是BRDF公式写错了。后来把采样点的散点图画出来才恍然大悟。从此我的规则是：任何采样器写完第一件事，画散点图。

第二个坑是低差异序列配合拒绝采样。我最初为了省事，在PBRT的复现里用了拒绝采样来处理圆盘，结果Sobol序列被破坏得非常厉害，渲染画面出现了一堆规律性的条纹噪声。后来换成同心映射，同样的样本预算，画面立刻干净了一个档次。低差异序列和拒绝采样天生八字不合，这个教训很深刻。

第三个坑是pdf()写错。圆盘半径为R时，pdf = 1/(πR²)，不是1/π。我在某个实验里把半径设成了2.0，但密度忘了除以4，结果蒙特卡洛估计整体偏大了4倍。这类错误无法通过散点图发现，因为散点图看起来还是均匀的，只有在算积分结果时才会露馅。所以接口设计里必须强制返回pdf，并写单测验证∫pdf = 1。

第四个坑是角度兼容问题。np.arctan2(y, x)返回的范围是[-π, π)，而采样代码里的角度范围是[0, 2π)。如果做角度直方图时不把theta归一化到一致的范围，直方图画出来会看到在±π边界处出现一个诡异的双峰，让人误以为角度分布不均匀。不是采样器的问题，是测试代码自身的问题。

5.3 性能优化与工程化建议

如果采样器在渲染循环的最内层被调用，性能就需要认真对待。三个优化方向：

第一是向量化，这在上文已经实现过了。批量生成样本时，向量化版本比标量循环快几十倍，这是性价比最高的优化。

第二是预计算。如果你的采样器参数固定（半径、圆心不变），能提前算的东西就提前算。不过圆盘采样本身计算量很小，预计算收益通常不大，真正收益大的场景是配合着色器时避免重复计算三角函数。

第三是缓存。如果你多次需要同一组样本，比如相机采样里每个像素都用同一套镜头样本做低差异抖动复用，那就应该把采样结果缓存下来，而不是每次重新采样。PBRT的Sampler类里也有类似的样本复用机制。

另外提一点接口设计建议：尽量让采样器只依赖u1, u2两个输入，不要偷偷调用内部随机数生成器。这样好处是便于复现和调试——同样的输入永远得到同样的输出，配合低差异序列时也更容易控制全局样本的维度分配。我之前犯过把u1, u2的生成逻辑写进采样器内部、导致外部无法注入Sobol序列的错误，重构之后才真正体会到"采样器只做映射、不做随机数生成"这句工程格言的含义。

最后分享一个小技巧：如果你需要在极坐标方案和同心映射方案之间快速切换调试，可以抽一个公共接口，内部用参数控制走哪条路径。我一般在开发阶段会把两种实现都留着，测试代码同时跑两个版本，确保它们在数学上等价。等验证完毕，再在正式管线里只保留同心映射版本。

我个人在实际操作中的体会是：圆盘采样看起来简单，但它是渲染管线里被调用次数最多的基础函数之一，值得花一晚上把原理彻底搞透。以后无论遇到圆面光源、景深还是余弦采样，都能一眼看穿它背后用的是哪种方案、有没有写错密度。这种底层的确定性知识，比盲目堆采样数更能提升渲染质量。