本文主要是介绍《学一辈子光线追踪》 二 一维MC积分(一),希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!
蒙特卡洛光线追踪技术系列 见 蒙特卡洛光线追踪技术
积分是关于计算面积和体积的,所以如果我们想让第1章最大限度地混淆,我们可以将它以一个积分的形式组织起来。但有时积分是最自然、最干净的表达方式。渲染常常是这样一个问题。让我们看一个经典的积分:
在计算机科学符号中,我们可以这样写:
I = area( x^2, 0, 2 )
我们也可以这样写:
I = 2*average(x^2, 0, 2)
我们也可以这样写:
这表明了一种MC方法:
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include "time.h"double myRandom() {return rand() / (RAND_MAX + 1.0);
}int main() {srand(time(NULL));int inside_circle = 0;int inside_circle_stratified = 0;int N = 1000000;float sum = 0.0;for (int i = 0;i < N;i++) {float x = 2 * myRandom();sum += x*x;}printf("Regular Estimate of PI = %f \n", 2 * sum / N);system("pause");
}
正如所料,这产生了近似于代数的精确答案,I=8/3。但我们也可以这样计算像 log(sin(x)) 那样我们不能解析积分的函数。在图形学中,我们经常有我们可以计算但不能显式写下的函数,或者我们只能概率地计算的函数。实际上,这就是前两本书中的光线跟踪color() 函数——我们不知道在每个方向上看到的是什么颜色,但我们可以在任何给定维度上对其进行统计估计。
我们在前两本书中编写的随机程序的一个问题是,小光源会产生太多的噪声,因为我们的均匀采样不够频繁。如果我们发送更多的随机样本到光照下,我们可以减少这个问题,但是我们需要降低这些样本的权重,以适应过采样。我们如何调整?为此,我们需要概率密度函数的概念。
解释一下上面这一段,是什么意思呢,就是说,光源很小,我们发射的采样光线大多数根本打不到光源上,比如挨着的两个像素,而且这两个像素显示的内容并不是某物体的边界处,两个像素各自发射100个采样光线,在随机散射下,其中一个像素有80个采样光线最终打中了光源,比较亮,但是另一个像素很不巧,只打中了20个,就很暗,但按理来说它们挨着,又不是边界,明暗应该很接近才对,所以产生很强烈的噪声(当然如果你每个像素发射10000个光线,随机噪声就会大大降低,但是也更浪费时间,是只发射100个光线的渲染速度的100倍!)。所以我们应该让像素发射的光线尽可能在能打到光的方向进行散射,换成这里的积分的例子,即我们应该要让更多的样本产生在x比较大的地方,因为在积分中,大的x对积分的准确性影响更大,但是我们又不能将所有x都取大的,不然得到的积分值肯定大于正常值,所以需要有一个权重来衡量:即我们要产生更多对最后结果有用的样本,同时还要让这些样本通过某种权重防止其过采样。
这样说肯定还是有些抽象,我用这个例子再详细说明一下,比如,我要计算x^2从0到2的积分,我们真实地计算一下发现,从0到1积分起来的值只有1/3,而从1到2积分的值是7/3,所以说,从1-2的积分的值不就占比重更大吗?当我们要通过估计的方法来获得值的时候,假如我们随机产生200个数,其中100个在(0-1)之间,100个在(1-2)之间,假设通过100个数的估计偏差在2%左右,则从0-1积分的估计值与真实值的偏差是0.0067(即(1/3)*0.02),而从1-2积分值的估计值与真实值的偏差为 0.047(即(7/3)*0.02),两个偏差加起来的值是0.053左右。而如果我们进行一下分配,让随机数在1-2之间产生的更多,这样就会降低1-2之间积分的偏差。根据刚才的计算,我们知道1-2之间占积分值比重很大,所以这样就能很好的降低总的偏差。
但是这样分两段(即0-1和1-2)不够完美,更好的方式是占积分比越大的值越容易产生,即在0-2积分中,2产生的概率大于1.99大于1.98这样,但是我们产生的这个值必须要除以一定的权重,不然最后的总积分值一定是会远大于原值的(就是肯定是不对的)。
注意,我们以下的部分都是在研究怎么样通过均匀的随机数发生器来产生我们需要的非均匀随机数。
首先,什么是密度函数?这只是直方图的连续形式。以下是来自直方图 Wikipedia 页面的示例:
如果我们为更多的树添加数据,直方图将变得更高。如果我们把数据分成更多的箱子,它会变得更短。与直方图的不同之处在于,它将频率y轴规格化为分数或百分比 (只是分数乘以100) 。我们把箱子的数量取到无穷大,一个连续的柱状图,因为所有箱子的高度都会降到零不能是分数。密度函数是一个我们把箱子拿出来调整它们的函数,这样当我们增加更多的箱子时它们不会变短。对于上面的树直方图,我们可以尝试:
Bin-height = (Fraction of trees between height H and H’) / (H-H’)
那就行了!我们可以将其解释为树高度的统计预测:
Probability a random tree is between H and H’ = Bin-height*(H-H’)
如果我们想知道在多个箱子里的概率,我们可以求和,概率密度函数,从pdf开始,是连续的分数直方图。
让我们制作一个pdf,并用它了解更多。假设我想要一个介于0和2之间的随机数r,其概率与自身成正比:r。我们希望pdf p(r)看起来类似于下图。但是应该有多高?
高度只有p(2)。那应该是什么?我们可以按惯例合理地做任何事,我们应该挑选一些方便的东西。就像直方图一样,我们可以对区域进行求和(积分)以计算r在某个区间(x0,x1)中的概率:
Probability x0 < r < x1 = C * area(p(r), x0, x1)
其中C是标度常数。我们不妨为了整洁,设C=1,这正是概率论所做的。我们知道r的概率在某个地方有值1,所以对于这个例子:
area(p(r), 0, 2) = 1.
因为 p(r) 与 r 成正比,也就是说,p=C'r对于另一个常数C' ,(注意该C'只是表示一个常数,与前面的C没有任何联系与关系)
area(C’r, 0, 2) = INTEGRAL_0^2 C’r dr = (½)C’r^2_0^2 C’r = C’/2
所以 p(r) =r/2 :生成数 r (r在0-2之间)的概率密度是r/2,注意这不是概率,因为这是有区间的,你要生成的数字是连续的,所以生成某一个数的几率本质是0。
我们如何用这个pdf p(r)生成一个随机数?为此我们需要更多的数学原理,别担心这不会永远持续下去。
给定 d=drand48() 中的一个均匀且介于0和1之间的随机数,我们应该能够找到一些函数 f(d) 来满足我们的需要。假设e = f(d) = d^2, 这不再是一个统一的pdf。e的pdf在0附近比在1附近大(将0和1之间的数字平方使其变小,如果f(d)=d的话,则在0和1的pdf是一样的)。要对函数进行一般观测,我们需要累积概率分布函数 P(x) :
P(x) = area(p, -infinity, x)
注意,对于没有定义P(x)的x,P(x)=0,即x的概率为零。对于我们的示例pdf p(r)=r/2,P(x)是:
P(x) = 0 : x < 0
P(x) = ¼ x^2 : 0 < x < 2
P(x) = 1 : x > 2
一个问题是,x和r是怎么回事?它们是伪变量-类似于程序中的函数参数。如果我们在x=1.0(书中写的0.5,这里我强烈认为是1.0)时计算P,我们得到:
P(1.0) = ¼
这表示pdf中随机变量小于1的概率为25%。这就产生了一个聪明的观察结果,它是产生非均匀随机数的许多方法的基础。我们需要一个函数f() ,当我们使用f(drand48()) 时,我们得到一个返回值,它的pdf为1/4 x^2。我们不知道这是什么,但我们知道,比如我们要通过该函数产生10000个结果,其中25%的结果应该小于1,75%的结果应该大于1。如果 f() 是递增的,那么我们期望 f(0.25)=1.0。这可以推广到计算每个可能输入的 f() :
f(P(x)) = x
这意味着f只是撤销P所做的一切。所以,
f(x) = P^-1 (x)
-1表示“反函数”,符号很丑,但比较标准。就我们的目的而言,这意味着,如果我们有pdf p() 和它的累积分布函数P(),那么如果我们对一个随机数这样做,我们将得到我们想要的:
e = P^-1(drand48())
对于我们的pdf p(x)=x/2和对应的P(x),我们需要计算p的逆。如果我们有:
y = ¼ x^2
我们通过用y来解x来求逆:
x = sqrt(4y)
因此,我们得到密度为p的随机数:
e = sqrt(4*drand48())
请注意,范围确实是0到2,如果我们的drand48()发送0.25,我们将得到1。注意在这里我们要产生的随机数x的几率随着x值变大而线性变大。
这样,我们得到的随机数x的几率就随着x变大而线性变大了。
这篇关于《学一辈子光线追踪》 二 一维MC积分(一)的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!