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

作者: 闪电的蓝熊猫 | 来源:发表于2019-11-30 11:19 被阅读0次

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

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

还是跟之前一样，本文是阅读《基于物理的渲染：从理论到实践》的总结，文中不会涉及到方方面面，只会整理笔者认为重要的一两个点。本文关注的重点是误差，尤其是舍入误差（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}$ 。这也就说明了，我们没法表示 $1$ 和 $1+2^{-23}$ 之间的数字。而这个 $2^{-23}$ 就是通常所说的机器精度（machine epsilon），我们用 $\epsilon$ 表示。

很少有数字正好落在浮点数能精确表示的位置上，不在这个位置上的数怎么办呢？就比如说 $1+2^{-25}$ 怎么表示？有两种方法可以采用，一是多出来的精度我就不管了，直接舍弃。二是我先看看这个精确的数字更接近哪个数，然后表示成那个数。 $1$ 和 $1+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

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

误差的来源

舍入误差

误差的表示

误差计算

误差的累积

一个简洁的记号

误差计算的实例

总结

参考资料

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读

OpenGL