【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}$ 。这也就说明了，我们没法表示 $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

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 203,271评论 5赞 476
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 85,275评论 2赞 380
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 150,151评论 0赞 336
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 54,550评论 1赞 273
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 63,553评论 5赞 365
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 48,559评论 1赞 281
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 37,924评论 3赞 395
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 36,580评论 0赞 257
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 40,826评论 1赞 297
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 35,578评论 2赞 320
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 37,661评论 1赞 329
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 33,363评论 4赞 318
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 38,940评论 3赞 307
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 29,926评论 0赞 19
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 31,156评论 1赞 259
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 42,872评论 2赞 349
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 42,391评论 2赞 342