NSDT工具推荐： Three.js AI纹理开发包 - YOLO合成数据生成器 - GLTF/GLB在线编辑 - 3D模型格式在线转换 - 可编程3D场景编辑器 - REVIT导出3D模型插件 - 3D模型语义搜索引擎 - AI模型在线查看 - Three.js虚拟轴心开发包 - 3D模型在线减面 - STL模型在线切割
本文介绍在三维渲染中的蒙特卡洛积分（Monte Carlo Integration）的两个应用。
1、使用蒙特卡洛积分绘制 McBeth 图
在颜色和数字图像课程中，我们解释了颜色可以表示为曲线，给出可见光谱中每个波长（用希腊字母 lambda 表示）的光量（这些曲线称为光谱功率分布或 SPD）。
图 1 显示了三种不同材料的此类曲线示例。如你所见，这些曲线可以解释为波长的函数：f(λ)。通常，可见光谱被认为在 380 到 730 纳米的区间内。要将此曲线转换回三刺激值 X、Y 和 Z，我们需要计算三个一维积分：
其中函数X_, Y_, Z_ 是 CIE 标准观察者颜色匹配函数。然后可以将得到的三刺激值 X、Y 和 Z 转换为 RGB。我们在颜色和数字图像课程中详细解释了这个过程。
希望你已经猜到我们将使用蒙特卡洛积分来解决这些积分问题。我们将使用的 SPD 是 McBeth 图上使用的 24 种不同颜色的光谱数据。这些颜色的光谱数据在网上很容易找到。如果你还记得我们在上一章中关于蒙特卡罗积分的内容，那么这个估计器的定义如下：
在我们的例子中，区间 [a,b] 是定义可见光谱的波长范围 [380,730]。为了近似此颜色检查器每种颜色的 XYX 值，我们需要在此区间内选择一个随机波长，评估并相乘此随机波长的 spd 和颜色匹配函数，重复此过程 N 次（其中 N 是样本数），将所有结果相加，将最终数字除以 N，然后乘以 (b-a)。在伪代码中我们得到：
int N = 32;
float lambdaMin = 380, lambdaMax = 730;
float X = 0, Y = 0, Z = 0;
for (int i = 0; i < N; ++i) { // for each sample
// select a random wavelength
float lambda = lambdaMin + drand48() * (lambdaMax - lambdaMin);
// evaluate the sod and the color-matching function at this wavelength,
// multiply the results, and add the contribution of these samples to the others
X += spd(lambda) * CIE_X(lambda);
Y += spd(lambda) * CIE_Y(lambda);
Z += spd(lambda) * CIE_Z(lambda);
}
// complete the integral, multiply by (b-a) and divide by N
X *= (lambdaMax - lambdaMin) / N;
Y *= (lambdaMax - lambdaMin) / N;
Z *= (lambdaMax - lambdaMin) / N;
简单吗？确实如此。
为了说明方差是什么，我们将使用 MC 积分将 MacBeth 图表的每种颜色渲染为 64x64 正方形，该正方形的每个像素（颜色桶的每个像素都是通过运行新的蒙特卡罗模拟获得的，因此每个像素可能与其他像素“略有”不同。差异程度主要取决于模拟中使用的样本数 N）。
我们的最终目标是重新创建看起来像真实 McBeth 图表的东西，但是由 CG 生成（当然使用蒙特卡罗积分）。在第 5 课中，我们提供了一个小程序的源代码，该程序使用黎曼和来计算积分来渲染 McBeth 图表。我们将从这个程序开始，因为它已经包含 McBeth 图表的光谱数据以及 CIE 颜色匹配函数。代码的核心是进行蒙特卡罗积分的函数。
unsigned int N = 32;
inline float linerp(const float *f, const short &i, const float &t, const int &max)
{ return f[i] * (1 - t) + f[std::min(max, i + 1)] * t; }
void monteCarloIntegration(const short &curveIndex, float &X, float &Y, float &Z)
{
float S = 0; // sum, used to normalize XYZ values
for (int i = 0; i < N; ++i) {
float lambda = drand48() * (lambdaMax - lambdaMin);
float b = lambda / 10;
short j = (short)b;
float t = b - j;
// interpolate
float fx = linerp(spd[curveIndex], j, t, nbins - 1);
b = lambda / 5;
j = (short)b;
t = b - j;
X += linerp(CIE_X, j, t, 2 * nbins - 1) * fx;
Y += linerp(CIE_Y, j, t, 2 * nbins - 1) * fx;
Z += linerp(CIE_Z, j, t, 2 * nbins - 1) * fx;
S += linerp(CIE_Y, j, t, 2 * nbins - 1);
}
// sum, normalization factor
S *= (lambdaMax - lambdaMin) / N;
// integral = (b-a) * 1/N * sum_{i=0}^{N-1} f(X_i) and normalize
X *= (lambdaMax - lambdaMin) / N / S;
Y *= (lambdaMax - lambdaMin) / N / S;
Z *= (lambdaMax - lambdaMin) / N / S;
}从此代码中可以看出，我们首先对区间 [380,730] 内的波长进行随机采样（第 10 行）。我们需要将此随机波长转换为索引，以便查找 CIE 配色函数（每 5 nm 采样一次）和 McBeth 光谱数据（每 10 nm 采样一次）一维表。由于波长可能介于表中的两个样本之间，因此最好执行线性插值（插入两个最近的样本）。最后，使用 CIE 配色函数值和此波长的光谱数据更新 X、Y 和 Z（第 19-21 行）。
由于这些值需要归一化（查看第 5 课），我们还需要计算与 Y 分量相对应的 CIE 配色函数的积分。在这里，我们也可以使用蒙特卡洛积分（第 22 行）。最后，将 X、Y 和 Z 值乘以 (b-a)（本例中为最小和最大波长）并除以 N（第 27-29 行）。我们还必须对它们进行归一化，但这与 MC 积分无关（它只是将光谱数据转换为 XYZ 的过程的一部分）。
const double XYZ_to_RGB[][3] = {
{ 2.3706743, -0.9000405, -0.4706338},
{-0.5138850, 1.4253036, 0.0885814},
{ 0.0052982, -0.0146949, 1.0093968}
};
void XYZtoRGB(const float &X, const float &Y, const float &Z, float &r, float &g, float &b)
{
r = std::max(0., X * XYZ_to_RGB[0][0] + Y * XYZ_to_RGB[0][1] + Z * XYZ_to_RGB[0][2]);
g = std::max(0., X * XYZ_to_RGB[1][0] + Y * XYZ_to_RGB[1][1] + Z * XYZ_to_RGB[1][2]);
b = std::max(0., X * XYZ_to_RGB[2][0] + Y * XYZ_to_RGB[2][1] + Z * XYZ_to_RGB[2][2]);
}
int main(int argc, char **argv)
{
...
for (int y = 0, pixIndex = 0; y < height; ++y) {
short row = (short)(y / squareSize);
for (int x = 0; x < width; ++x, pixIndex++) {
short column = (short)(x / squareSize);
// compute each pixel of the square
float X = 0, Y = 0, Z = 0;
monteCarloIntegration(row * 6 + column, X, Y, Z);
float r, g, b;
XYZtoRGB(X, Y, Z, r, g, b);
buffer[pixIndex][0] += r;
buffer[pixIndex][1] += g;
buffer[pixIndex][2] += b;
}
}
...
}
程序的主要功能是循环遍历帧中的每个像素。对于每个像素，它计算该像素在 McBeth 图中所对应的行和列（图表的方块被组织成一系列 4 行和 6 列。您可以在图 2 中看到 McBeth 图的示例）。该颜色的光谱数据将用于计算 X、Y 和 Z 的值（第 22 行和第 23 行）。然后将它们转换为 RGB 值，我们也在“光、颜色和颜色空间简介”课程中解释了这个过程。最后，这些 RGB 值存储在缓冲区中，然后在处理完所有像素后将其保存到磁盘（你的屏幕可能希望图像为 sRGB，因此我们还会在将像素值存储到磁盘之前对其应用伽马）。图 3 显示了 32 个样本的结果。现在，正如您所见，结果远非糟糕，但图像相当嘈杂！这就是著名的方差，我们在蒙特卡罗方法的数学基础课程中已经讨论过。或者当你查看渲染时，可能更常称之为噪音。
2、渐进式渲染示例
为了演示蒙特卡罗积分的另一个非常好的属性，我们将修改我们的程序，以便它通过计算我们想要的尽可能多的图像版本（我们将这些图像称为通道）并平均其结果来不断优化结果。
在某种程度上，这与增加样本数量 N 相同。但是，如果将 N 增加到 4096，则图像可能需要更长的时间来渲染（实际上，从技术上讲，比 N = 32 的图像长 128 倍）。因此，当 N 较小时，我们可以立即看到结果，如果我们喜欢这个快速的第一遍之后看到的结果，那么我们可以让程序运行并计算更多遍，直到我们最终对图像质量感到满意（即当您认为噪音不再分散注意力时）。这种方法称为渐进式渲染（Progressive Rendering）。
假设你需要在一定时间内提供一系列图像。如果 X 是将总渲染时间除以需要渲染的图像总数后得到的分钟数，那么你可以让程序渲染这些图像中的每一个不超过 X 分钟。这样，你就可以确定所有图像都会按时渲染。它们可能会有一些噪音，但你将拥有完整的图像序列可供查看，并且在某些情况下，这通常比拥有不完整的图像序列要好，即使它们看起来要好得多。渲染时间是可预测的，这一事实在生产环境中通常很有用。
此外，如果你喜欢经过 X 小时渲染后获得的图像序列，则可以重新启动该序列并让其再细化 X 小时。程序无需从头开始渲染每个图像，只需从现有渲染开始，并通过添加更多遍历来不断细化它们即可。因此，不会浪费任何渲染资源！这些是渐进式渲染的一些巨大好处。
更正式地说，请记住，蒙特卡罗积分就像计算随机变量的样本均值一样。在蒙特卡罗方法的数学基础课程中，我们详细解释了这些样本的均值（它们本身是随机数）也可以取平均值以产生更准确的结果。在此示例中，遍历只不过是样本均值。并且这些遍历的结果可以取平均值以产生更准确的结果（与样本均值一样）。这就是蒙特卡罗积分使渐进式渲染成为可能的原因。
我们在程序中实现了渐进式渲染的基本形式。这次我们不会在 OpenGL 窗口中显示结果。该程序完全在 shell 中运行。它将继续渲染过程并将结果累积在缓冲区中，直到用户终止程序（使用 ctrl-c）。程序将捕获此信号并调用一个函数，我们将在退出之前将缓冲区的结果保存到磁盘中。请注意，实际上，在每次渲染之后，我们实际上会将缓冲区的内容转换为适当的图像，其中缓冲区的内容除以总渲染次数（缓冲区就像一个累积缓冲区）。保存到磁盘的是图像的内容，而不是缓冲区。以下是源代码：
void handler(int s)
{
printf("Saving image after %d passes\n", npasses);
std::ofstream ofs;
ofs.open("./mcbeth.ppm");
ofs << "P6\n" << width << " " << height << "\n255\n";
for (uint32_t i = 0; i < width * height; ++i) {
unsigned char r, g, b;
r = (unsigned char)(std::min(1.f, pixels[i][0]) * 255);
g = (unsigned char)(std::min(1.f, pixels[i][1]) * 255);
b = (unsigned char)(std::min(1.f, pixels[i][2]) * 255);
ofs << r << g << b;
}
ofs.close();
exit(0);
}
int main(int argc, char **argv)
{
pixels = new color[width * height];
memset(pixels, 0x0, sizeof(color) * width * height);
buffer = new color [width * height];
memset(buffer, 0x0, sizeof(color) * width * height);
struct sigaction sigIntHandler;
sigIntHandler.sa_handler = handler;
sigemptyset(&sigIntHandler.sa_mask);
sigIntHandler.sa_flags = 0;
sigaction(SIGINT, &sigIntHandler, NULL);
static const float gamma = 1 / 2.2;
while (1) {
for (int y = 0, pixIndex = 0; y < height; ++y) {
short row = (short)(y / squareSize);
for (int x = 0; x < width; ++x, pixIndex++) {
short column = (short)(x / squareSize);
// compute each pixel of the square
float X = 0, Y = 0, Z = 0;
monteCarloIntegration(row * 6 + column, X, Y, Z);
float r, g, b;
XYZtoRGB(X, Y, Z, r, g, b);
buffer[pixIndex][0] += r;
buffer[pixIndex][1] += g;
buffer[pixIndex][2] += b;
}
}
npasses++;
for (uint32_t i = 0; i < width * height; ++i) {
pixels[i][0] = powf(buffer[i][0] / npasses, gamma);
pixels[i][1] = powf(buffer[i][1] / npasses, gamma);
pixels[i][2] = powf(buffer[i][2] / npasses, gamma);
}
printf("npasses %3d (num samples %5d)\n", npasses, npasses * N);
}
delete [] pixels;
delete [] buffer;
return 0;
}让我们看一些结果。以下图像序列显示了 1、8、16 和 32 次传递后的结果（分别对应于 32、256、512 和 1024 个样本）：
显然，正如预期的那样，随着遍历次数（即样本数量）的增加，噪声会减少。通过缩小这些图像的一部分，我们可以更好地看到噪声。请记住，在蒙特卡罗中，你需要 4 倍的样本才能将噪声（或方差）减少 2。
正如本课中已经多次提到的，蒙特卡罗渲染的艺术主要是寻找减少这种噪音的方法。我们将在本课后面讨论方差减少技术。
3、结束语
希望通过这个练习，你能更好地理解什么是蒙特卡罗积分，正如你所见，它不仅是一个简单的想法，而且非常容易实现。
本文还提供一个简单的程序来演示渐进式渲染的概念（这只有在蒙特卡罗算法的性质下才有可能）。希望你现在能更好地理解渲染中的噪音（或方差，如果使用正确的技术术语）来自哪里。如果在我们的 3D 渲染中看到这种噪音，可以肯定生成此图像的渲染器使用了某种蒙特卡罗积分。
原文链接：Monte Carlo in Rendering (A Practical Example)
BimAnt翻译整理，转载请标明出处
