【PBRT】《基于物理的渲染:从理论到实践》梳理III

四百万个三角形面片的模型

上面这张图主要是用来吸引人的,跟本文内容没啥关系(偷笑~)

还是跟之前一样,本文是阅读《基于物理的渲染:从理论到实践》的总结,文中不会涉及到方方面面,只会整理笔者认为重要的一两个点。本文关注的重点是误差,尤其是舍入误差(Rounding Error),主要来源书中的3.9节。

由于作者水平所限,文中如果有错误或者没有解释到位的地方,还请读者不吝指正。

误差的来源

首先我们必须认识到,没有任何方法可以避免误差,我们能做的就只能是减少误差。为啥?从误差的来源我们就可以看出来了。

误差会来自两个方面:

  • 第一,用于计算的数据(即输入数据)本身就不是精确的,即输入数据本身就带有误差。

  • 第二,保存数据的工具会产生误差。这个怎么理解?我们用计算机来保存数据,计算机是二进制的,我们用来保存数据的大小也有限制,不可能表示数学上类似于“实数”这种无限精度的数字。

从误差的来源我们就可以看出,误差无法避免,无论是哪个方面的误差我们都无法避免,只能尽可能的减少误差。这就意味着,当我们进行计算的时候(尤其是要进行光线和表面交点这种精度要求非常高的计算),我们必须把误差考虑进去,否则渲染出来的场景会和理想中的场景相差很大。而输入数据本身的精确性我们没法控制,就当它是精确的吧!我们要让第二个误差,也就是保存数据的工具产生的误差减少!

舍入误差

舍入误差,英文名叫rounding error,是用计算机保存实数的时候,由于精度有限,导致的保存的数据和原来的数据的差距。这个差距,就是舍入误差。

就拿32位浮点数来说。IEEE的标准(下面讲的内容都是IEEE标准)是32位浮点数,首位用来保存符号+/-(用s表示),接下来的8位用来保存指数值e,最后的23位用来保存有效数字m(二进制0或者1),即:

32位浮点数

于是,我们的数字可以表示成。如果你的观察够敏锐的话,会发现,这种表示方式无法表示0,这是不能容忍的,所以IEEE又规定,如果指数e是0,那么表示的数字将是,这样,如果m等于0的话,我们就能表示0了,虽然有+0和-0之分。当然,这时候IEEE又跳出来了,规定说+0和-0的表现必须是一样的,这样才有了唯一的0值。

绝大多数情况下,我们的指数e不会等于0,所以我们可以放心大胆地只考虑第一种情况。考虑这样一个问题,浮点数1和离1最近的比1大的浮点数之间的间隔是多少呢?

我们来看一下,表示1的时候,m值是00...00,然后e是127,这样求出来1.0*2^0=1。而比这个数稍微大点的数就是m=0...01在第23位的地方是1,这就比1大了,也是比1大的最小数,它比1大了2^{-23}。这也就说明了,我们没法表示11+2^{-23}之间的数字。而这个2^{-23}就是通常所说的机器精度(machine epsilon),我们用\epsilon表示。

很少有数字正好落在浮点数能精确表示的位置上,不在这个位置上的数怎么办呢?就比如说1+2^{-25}怎么表示?有两种方法可以采用,一是多出来的精度我就不管了,直接舍弃。二是我先看看这个精确的数字更接近哪个数,然后表示成那个数。11+2^{-23}的中间量是1+2^{-24}1+2^{-25}显然要比这个中间量小,所以它会被表示成1

那么,正好落在中间的数怎么办?比如说1+2^{-24}?这又是一个麻烦事,还是得寄出IEEE大法:

IEEE舍入最近法则:
如果第24位,为1,第24位后的所有已知位是0,那么如果第23位是1,就往第23位加1,如果第23位为0,就不加。
如果第24位为1,第24位之后还有至少一个1,那么就往第23位上加1.
如果第24位为0,那么就舍弃之后的所有位数,保持第23位的数不变。

很绕的一个方法,解释一下。

情况一、表示1+2^{-25}
因为这个数第24位是0,第25位是1,根据IEEE舍入最近法则,第24位为0的话,保持第23位数不变的原则,第25位会被舍弃,而这个数就会被舍入成1

情况二、表示1+2^{-24}
这个数第24位是0,之后的所有数位都是0,而1的第23位是0,所以不会进位,导致它也会被表示成1

情况三、表示1+2^{-23}+2^{-24}
这个数的第24位是1,并且之后所有的数位都是0,而1+2^{-23}的第23位是1,所以它会进位,导致最后表示成1+2^{-22}

神奇吧?

误差的表示

知道舍入误差产生的原因之后,接下来就要表示误差了。在分析误差的时候,我们需要两个概念:绝对误差(absolute error)和相对误差(relative error)。

a是精确值,而\tilde{a}是浮点值
绝对误差\delta_a表示为
\delta_a=|\tilde{a}-a|
即,浮点值与精确值之差的绝对值。
相对误差\delta_r表示为
\delta_r=|\frac{\tilde{a}-a}{a}|
即,绝对误差与精确值绝对值之比。

由于舍入的规则,精确值和浮点值之间的相对误差不会超过半个机器精度,即\frac{\epsilon}{2},大小是2^{-24}

fl(x)表示将精确值x舍入成浮点值,舍入后的浮点值通常表示成\tilde{x},就是在上面加个波浪线。然后,用符号来表示浮点算术运算:
a\oplus{b}=fl(a+b)
a\ominus{b}=fl(a-b)
a\otimes{b}=fl(a\times{b})
a\oslash{b}=fl(a/b)

误差计算

确定了符号之后,紧接着进行误差的计算。标准的误差计算模型是
fl(x\quad op\quad y)=(x\quad op\quad y)(1+\mu), \quad \quad |\mu|\leq \frac{\epsilon}{2}, \quad\quad op=+,-,*,/
每当发生一次浮点运算的时候,引入一个\mu的相对误差,这个误差的大小不超过机器精度的一半。这是一个非常合理的模型,用来分析误差大小非常合适。

误差的累积

考虑这样一个计算:a\oplus b\oplus c\oplus d,假设a,b,c,d都是精确值。在计算机中,我们可以认为它的计算顺序是这样的(((a\oplus b)\oplus c)\oplus d),于是它的误差可以表示成:
(((a+b)(1+\mu)+c)(1+\mu)+d)(1+\mu)
展开式子可得
(a+b)(1+\mu)^3+c(1+\mu)^2+d(1+\mu)
这个式子的最大误差是
|a+b|(1+\mu)^3+|c|(1+\mu)^2+|d|(1+\mu)
但是这个表示过于繁琐,我们考虑误差不想这么复杂大计算这么一个长式子,所以把它简化为:
(|a+b|+|c|+|d|)(1+\mu)^3
虽然扩大了误差,但是简化了计算。还有就是要注意,这里必须要有绝对值,因为如果a+b, c, d中有负数的话,那么相加之后的误差会减少,只有当这些数的符号都相同的时候,它才是最大误差,所以我们在这里加上绝对值符号。比如说,考虑如下的情况:

假设有两个数相加,a=3,b=-2,\mu=0.001,计算表达式为:a(1+\mu)^3+b(1+\mu)^2,我们可以得到最终结果是1.005007003,误差是0.005007003,用绝对值限制的误差结果是0.015015005,真实误差在绝对值限制误差之内,而如果没有绝对值,得到的误差会是0.003003001,显然真实误差要超过限制了,这是不允许的。

这里要指出书上的一处错误,书上说如果a(1+\mu)^i\oplus b(1+\mu)^j中,a和b的符号相同的话,那么最大误差为|a+b|(1+\mu)^{i+j+1},这很明显太大了。从它的实现代码来看,它用的也不是这个范围,而是我上面说的那个。

对乘法的计算来说,考虑误差的计算非常直观:
a(1+\mu)^i\otimes b(1+\mu)^j\subset ab(1+\mu)^{i+j+1}
因为本身a\otimes b就要引入1+\mu的误差,所以最后的误差指数需要i+j再加1

除法的计算就不那么直观了:
\frac{a(1+\mu)^i}{b(1+\mu)^j}\subset \frac{a}{b}(1+\mu)^{i+j+1}
你没看错,是i+j而不是i-j,因为这里的\mu有正有负,它是一个范围,所以不能把它消掉,它的误差和相乘的误差是一样的!

一个简洁的记号

根据从其他书上借来的引理,如果|\mu|<\frac{\epsilon}{2},并且对正整数n有n\mu<1,可以得到以下式子:
(1+\mu)^n=1+\theta_n,\quad 且|\theta_n|\leq\frac{n\mu}{1-n\mu}=\gamma_n
这样,a(1+\mu)^n就可以简单地记作a(1+\gamma_n)了。

误差计算的实例

根据上面的原理,我们可以看到pbrt中有正确的误差,有错误的误差,还有不知道是不是正确的误差(-_-)!

先看正确的误差:
用计算公式计算出球表面和光线的交点后,还需要执行一步,把交点投射到球表面上去,因为由于误差的缘故,计算出来的交点可能不在表面上。这个投射的过程是
x'=x\frac{r}{\sqrt{x^2+y^2+z^2}}
x'是投射后的x坐标,投射y和z坐标的方式一样。我们来分析一下这个操作的误差
\begin{equation} \begin{split} x'&=x\otimes r\oslash sqrt((x\otimes x)\oplus(y\otimes y)\oplus(z\otimes z))\\ &\subset \frac{xr(1+\gamma_1)}{\sqrt{x^2(1+\gamma_3)+y^2(1+\gamma_3)+z^2(1+\gamma_2)}(1+\gamma_1)}\\ &\subset\frac{xr(1+\gamma_1)}{\sqrt{(x^2+y^2+z^2)(1+\gamma_4)}(1+\gamma_1)}\\ &=\frac{xr(1+\gamma_1)}{\sqrt{x^2+y^2+z^2}(1+\gamma_3)}\\ &\subset \frac{xr}{\sqrt{x^2+y^2+z^2}}(1+\gamma_5) \end{split} \end{equation}
其中,\sqrt{}操作引入的误差为\gamma_1,然后为了开根号,我们把分母内的误差扩大了\gamma_1,使得最后的式子非常简洁。可以查看pbrt-v3的代码,我们可以看到:

// Compute error bounds for sphere intersection
Vector3f pError = gamma(5) * Abs((Vector3f)pHit);

没错,就是\gamma_5

再来看个错误的误差:
在用矩阵转换顶点的时候,代码是这样的

T xp = m.m[0][0] * x + m.m[0][1] * y + m.m[0][2] * z + m.m[0][3];

而书中的式子是这样的



注意这里在后面加了个括号,代码中是没有这个括号的,有了括号之后,计算出来的误差和没有括号计算的误差是不一样的。

我们先来看有括号之后的计算过程是什么
\begin{equation} \begin{split} x'&=((m_{0,0}\otimes x)\oplus(m_{0,1}\otimes y))\oplus((m_{0,2}\otimes z)\oplus m_{0,3})\\ &\subset m_{0,0}x(1+\gamma_3)+m_{0,1}y(1+\gamma_3)+m_{0,2}z(1+\gamma_3)+m_{0,3}(1+\gamma_2)\\ &\subset (|m_{0,0}x|+|m_{0,1}y|+|m_{0,2}z|+|m_{0,3}|)(1+\gamma_3) \end{split} \end{equation}
最终引入的误差是\gamma_3

不加括号是什么样子?
\begin{equation} \begin{split} x'&=(m_{0,0}\otimes x)\oplus (m_{0,1}\otimes y)\oplus (m_{0,2}\otimes z)\oplus m_{0,3}\\ &\subset m_{0,0}x(1+\gamma_1)\oplus m_{0,1}y(1+\gamma_1)\oplus m_{0,2}z(1+\gamma_1)\oplus m_{0,3}\\ &\subset (((m_{0,0}x(1+\gamma_1)+m_{0,1}y(1+\gamma_1))(1+\gamma_1)+m_{0,2}z(1+\gamma_1))(1+\gamma_1)+m_{0,3})(1+\gamma_1)\\ &=m_{0,0}x(1+\gamma_4)+m_{0,1}y(1+\gamma_4)+m_{0,2}z(1+\gamma_2)+m_{0,3}(1+\gamma_1)\\ &\subset (|m_{0,0}x|+|m_{0,1}y|+|m_{0,2}z|+|m_{0,3}|)(1+\gamma_4) \end{split} \end{equation}
没错,最终引入的误差是\gamma_4。而它的代码中,引入的误差却仍然是\gamma_3

T xAbsSum = (std::abs(m.m[0][0] * x) + std::abs(m.m[0][1] * y) +
            std::abs(m.m[0][2] * z) + std::abs(m.m[0][3]));
*pError = gamma(3) * Vector3<T>(xAbsSum, yAbsSum, zAbsSum);

最后是不知道正确不正确的误差:
计算光线和三角形交点的时候,涉及到的计算非常复杂,最终的计算公式是


但是里面的所有参数,都是从别的地方计算出来的,也就是说这些参数本身就包含误差在里面。误差的计算过程也弄不清楚,为啥要分别计算x,y,z分量的绝对误差值,然后再计算t值的误差,完全搞不懂啊!

总结

说了这么多,你一定觉得误差非常恐怖。其实不然,因为不管误差怎么累积,它的量级也不会太大,除非你放大误差。跟分析误差相比,更应该关注的是算法的稳定性和精确度,因为如果算法不好,它会放大某些误差,导致结果完全不一样。所以,更多地关注算法,在算法出现问题的时候,从误差的角度去思考一下,是一个不错的方法。

参考资料

Physically Based Rendering
pbrt源码,版本3
Accuracy and Stability of Numerical Algorithm
数值分析(第二版)
Rounding Errors in Algebraic Processes

推荐阅读更多精彩内容