【PBRT】《基于物理的渲染：從理論到實踐》梳理III

四百萬個三角形面片的模型

上面這張圖主要是用來吸引人的考廉，跟本文內(nèi)容沒啥關(guān)系（偷笑~）

還是跟之前一樣秘豹，本文是閱讀《基于物理的渲染：從理論到實踐》的總結(jié)，文中不會涉及到方方面面昌粤，只會整理筆者認為重要的一兩個點既绕。本文關(guān)注的重點是誤差，尤其是舍入誤差（Rounding Error）涮坐，主要來源書中的3.9節(jié)凄贩。

由于作者水平所限，文中如果有錯誤或者沒有解釋到位的地方袱讹，還請讀者不吝指正怎炊。

誤差的來源

首先我們必須認識到，沒有任何方法可以避免誤差廓译，我們能做的就只能是減少誤差评肆。為啥？從誤差的來源我們就可以看出來了非区。

誤差會來自兩個方面：

第一瓜挽，用于計算的數(shù)據(jù)（即輸入數(shù)據(jù)）本身就不是精確的，即輸入數(shù)據(jù)本身就帶有誤差征绸。
第二久橙，保存數(shù)據(jù)的工具會產(chǎn)生誤差俄占。這個怎么理解？我們用計算機來保存數(shù)據(jù)淆衷，計算機是二進制的缸榄，我們用來保存數(shù)據(jù)的大小也有限制，不可能表示數(shù)學(xué)上類似于“實數(shù)”這種無限精度的數(shù)字祝拯。

從誤差的來源我們就可以看出甚带，誤差無法避免，無論是哪個方面的誤差我們都無法避免佳头，只能盡可能的減少誤差鹰贵。這就意味著，當(dāng)我們進行計算的時候（尤其是要進行光線和表面交點這種精度要求非常高的計算）康嘉，我們必須把誤差考慮進去碉输，否則渲染出來的場景會和理想中的場景相差很大。而輸入數(shù)據(jù)本身的精確性我們沒法控制亭珍，就當(dāng)它是精確的吧敷钾！我們要讓第二個誤差，也就是保存數(shù)據(jù)的工具產(chǎn)生的誤差減少肄梨！

舍入誤差

舍入誤差闰非，英文名叫rounding error，是用計算機保存實數(shù)的時候峭范，由于精度有限财松，導(dǎo)致的保存的數(shù)據(jù)和原來的數(shù)據(jù)的差距。這個差距纱控，就是舍入誤差辆毡。

就拿32位浮點數(shù)來說。IEEE的標準（下面講的內(nèi)容都是IEEE標準）是32位浮點數(shù)甜害，首位用來保存符號+/-（用s表示）舶掖，接下來的8位用來保存指數(shù)值e，最后的23位用來保存有效數(shù)字m（二進制0或者1）尔店，即：

32位浮點數(shù)

于是眨攘，我們的數(shù)字可以表示成。如果你的觀察夠敏銳的話嚣州，會發(fā)現(xiàn)鲫售，這種表示方式無法表示0，這是不能容忍的该肴，所以IEEE又規(guī)定情竹，如果指數(shù)e是0，那么表示的數(shù)字將是匀哄，這樣秦效，如果m等于0的話雏蛮，我們就能表示0了，雖然有+0和-0之分阱州。當(dāng)然挑秉，這時候IEEE又跳出來了，規(guī)定說+0和-0的表現(xiàn)必須是一樣的苔货，這樣才有了唯一的0值犀概。

絕大多數(shù)情況下，我們的指數(shù)e不會等于0蒲赂，所以我們可以放心大膽地只考慮第一種情況阱冶〉蟊铮考慮這樣一個問題滥嘴，浮點數(shù)1和離1最近的比1大的浮點數(shù)之間的間隔是多少呢？

我們來看一下至耻，表示1的時候若皱，m值是00...00，然后e是127尘颓，這樣求出來 $1.0*2^0=1$ 走触。而比這個數(shù)稍微大點的數(shù)就是 $m=0...01$ 在第23位的地方是1，這就比1大了疤苹，也是比1大的最小數(shù)互广，它比1大了 $2^{-23}$ 喻频。這也就說明了努咐，我們沒法表示 $1$ 和 $1+2^{-23}$ 之間的數(shù)字棍矛。而這個 $2^{-23}$ 就是通常所說的機器精度（machine epsilon）娃圆，我們用 $\epsilon$ 表示捂龄。

很少有數(shù)字正好落在浮點數(shù)能精確表示的位置上曹铃，不在這個位置上的數(shù)怎么辦呢绞佩？就比如說 $1+2^{-25}$ 怎么表示冬阳？有兩種方法可以采用颤霎，一是多出來的精度我就不管了媳谁，直接舍棄。二是我先看看這個精確的數(shù)字更接近哪個數(shù)友酱，然后表示成那個數(shù)晴音。 $1$ 和 $1+2^{-23}$ 的中間量是 $1+2^{-24}$ ， $1+2^{-25}$ 顯然要比這個中間量小缔杉，所以它會被表示成 $1$ 段多。

那么，正好落在中間的數(shù)怎么辦壮吩？比如說 $1+2^{-24}$ ?這又是一個麻煩事进苍，還是得寄出IEEE大法：

IEEE舍入最近法則：
如果第24位,為1加缘，第24位后的所有已知位是0，那么如果第23位是1觉啊，就往第23位加1拣宏，如果第23位為0，就不加杠人。
如果第24位為1勋乾，第24位之后還有至少一個1，那么就往第23位上加1.
如果第24位為0嗡善，那么就舍棄之后的所有位數(shù)辑莫，保持第23位的數(shù)不變。

很繞的一個方法罩引，解釋一下各吨。

情況一、表示 $1+2^{-25}$
因為這個數(shù)第24位是0袁铐，第25位是1揭蜒，根據(jù)IEEE舍入最近法則，第24位為0的話剔桨，保持第23位數(shù)不變的原則屉更，第25位會被舍棄，而這個數(shù)就會被舍入成 $1$ 洒缀。

情況二瑰谜、表示 $1+2^{-24}$
這個數(shù)第24位是0，之后的所有數(shù)位都是0树绩，而 $1$ 的第23位是0萨脑，所以不會進位，導(dǎo)致它也會被表示成 $1$ 葱峡。

情況三砚哗、表示 $1+2^{-23}+2^{-24}$
這個數(shù)的第24位是1，并且之后所有的數(shù)位都是0砰奕，而 $1+2^{-23}$ 的第23位是1蛛芥，所以它會進位，導(dǎo)致最后表示成 $1+2^{-22}$ 军援。

神奇吧仅淑？

誤差的表示

知道舍入誤差產(chǎn)生的原因之后，接下來就要表示誤差了胸哥。在分析誤差的時候涯竟，我們需要兩個概念：絕對誤差（absolute error）和相對誤差（relative error）。

令 $a$ 是精確值，而 $\tilde{a}$ 是浮點值
絕對誤差 $\delta_a$ 表示為
$\delta_a=|\tilde{a}-a|$
即庐船，浮點值與精確值之差的絕對值银酬。
相對誤差 $\delta_r$ 表示為
$\delta_r=|\frac{\tilde{a}-a}{a}|$
即，絕對誤差與精確值絕對值之比筐钟。

由于舍入的規(guī)則揩瞪，精確值和浮點值之間的相對誤差不會超過半個機器精度，即 $\frac{\epsilon}{2}$ 篓冲，大小是 $2^{-24}$ 李破。

用 $fl(x)$ 表示將精確值 $x$ 舍入成浮點值，舍入后的浮點值通常表示成 $\tilde{x}$ 壹将，就是在上面加個波浪線嗤攻。然后，用符號來表示浮點算術(shù)運算：
$a\oplus诽俯=fl(a+b)$
$a\ominus妇菱=fl(a-b)$
$a\otimes=fl(a\times惊畏)$
$a\oslash恶耽=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=+,-,*,/$
每當(dāng)發(fā)生一次浮點運算的時候，引入一個 $\mu$ 的相對誤差浪讳，這個誤差的大小不超過機器精度的一半缰盏。這是一個非常合理的模型，用來分析誤差大小非常合適淹遵。

誤差的累積

考慮這樣一個計算： $a\oplus b\oplus c\oplus d$ 口猜，假設(shè) $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)$
但是這個表示過于繁瑣，我們考慮誤差不想這么復(fù)雜大計算這么一個長式子辐真，所以把它簡化為：
$(|a+b|+|c|+|d|)(1+\mu)^3$
雖然擴大了誤差须尚，但是簡化了計算。還有就是要注意侍咱，這里必須要有絕對值耐床，因為如果 $a+b, c, d$ 中有負數(shù)的話，那么相加之后的誤差會減少楔脯，只有當(dāng)這些數(shù)的符號都相同的時候撩轰，它才是最大誤差，所以我們在這里加上絕對值符號。比如說堪嫂，考慮如下的情況：

假設(shè)有兩個數(shù)相加偎箫， $a=3,b=-2,\mu=0.001$ ，計算表達式為： $a(1+\mu)^3+b(1+\mu)^2$ 皆串，我們可以得到最終結(jié)果是 $1.005007003$ 镜廉，誤差是 $0.005007003$ ，用絕對值限制的誤差結(jié)果是 $0.015015005$ 愚战，真實誤差在絕對值限制誤差之內(nèi)娇唯，而如果沒有絕對值，得到的誤差會是 $0.003003001$ 寂玲，顯然真實誤差要超過限制了塔插，這是不允許的。

這里要指出書上的一處錯誤拓哟，書上說如果 $a(1+\mu)^i\oplus b(1+\mu)^j$ 中想许，a和b的符號相同的話，那么最大誤差為 $|a+b|(1+\mu)^{i+j+1}$ 断序，這很明顯太大了流纹。從它的實現(xiàn)代碼來看，它用的也不是這個范圍违诗，而是我上面說的那個漱凝。

對乘法的計算來說，考慮誤差的計算非常直觀:
$a(1+\mu)^i\otimes b(1+\mu)^j\subset ab(1+\mu)^{i+j+1}$
因為本身 $a\otimes b$ 就要引入 $1+\mu$ 的誤差诸迟，所以最后的誤差指數(shù)需要 $i+j$ 再加 $1$ 茸炒。

除法的計算就不那么直觀了：
$\frac{a(1+\mu)^i}{b(1+\mu)^j}\subset \frac{a}(1+\mu)^{i+j+1}$
你沒看錯阵苇，是 $i+j$ 而不是 $i-j$ 壁公，因為這里的 $\mu$ 有正有負，它是一個范圍绅项，所以不能把它消掉紊册，它的誤差和相乘的誤差是一樣的！

一個簡潔的記號

根據(jù)從其他書上借來的引理快耿，如果 $|\mu|<\frac{\epsilon}{2}$ 囊陡，并且對正整數(shù)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)$ 了。

誤差計算的實例

根據(jù)上面的原理铺浇，我們可以看到pbrt中有正確的誤差痢畜，有錯誤的誤差，還有不知道是不是正確的誤差（-_-）！

先看正確的誤差：
用計算公式計算出球表面和光線的交點后丁稀，還需要執(zhí)行一步吼拥，把交點投射到球表面上去，因為由于誤差的緣故线衫，計算出來的交點可能不在表面上凿可。這個投射的過程是
$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$ 枯跑，然后為了開根號，我們把分母內(nèi)的誤差擴大了 $\gamma_1$ 白热，使得最后的式子非常簡潔敛助。可以查看pbrt-v3的代碼屋确，我們可以看到：

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

沒錯纳击，就是 $\gamma_5$ 。

再來看個錯誤的誤差：
在用矩陣轉(zhuǎn)換頂點的時候攻臀，代碼是這樣的

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);

最后是不知道正確不正確的誤差：
計算光線和三角形交點的時候，涉及到的計算非常復(fù)雜纵东，最終的計算公式是

但是里面的所有參數(shù)粘招，都是從別的地方計算出來的，也就是說這些參數(shù)本身就包含誤差在里面偎球。誤差的計算過程也弄不清楚洒扎，為啥要分別計算x，y衰絮，z分量的絕對誤差值袍冷，然后再計算t值的誤差，完全搞不懂懊怠胡诗！

總結(jié)

說了這么多，你一定覺得誤差非常恐怖煌恢。其實不然骇陈，因為不管誤差怎么累積，它的量級也不會太大瑰抵，除非你放大誤差你雌。跟分析誤差相比，更應(yīng)該關(guān)注的是算法的穩(wěn)定性和精確度二汛，因為如果算法不好婿崭，它會放大某些誤差，導(dǎo)致結(jié)果完全不一樣肴颊。所以逛球，更多地關(guān)注算法，在算法出現(xiàn)問題的時候苫昌，從誤差的角度去思考一下颤绕，是一個不錯的方法。

參考資料

Physically Based Rendering
pbrt源碼祟身，版本3
Accuracy and Stability of Numerical Algorithm
數(shù)值分析（第二版）
Rounding Errors in Algebraic Processes

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者

人面猴
序言：七十年代末奥务，一起剝皮案震驚了整個濱河市，隨后出現(xiàn)的幾起案子袜硫，更是在濱河造成了極大的恐慌氯葬，老刑警劉巖，帶你破解...
沈念sama閱讀 219,539評論 6贊 508
死咒
序言：濱河連續(xù)發(fā)生了三起死亡事件婉陷，死亡現(xiàn)場離奇詭異帚称，居然都是意外死亡，警方通過查閱死者的電腦和手機秽澳，發(fā)現(xiàn)死者居然都...
沈念sama閱讀 93,594評論 3贊 396
救了他兩次的神仙讓他今天三更去死
文/潘曉璐我一進店門闯睹，熙熙樓的掌柜王于貴愁眉苦臉地迎上來，“玉大人担神，你說我怎么就攤上這事楼吃。” “怎么了妄讯？”我有些...
開封第一講書人閱讀 165,871評論 0贊 356
道士緝兇錄：失蹤的賣姜人
文/不壞的土叔我叫張陵孩锡，是天一觀的道長。經(jīng)常有香客問我亥贸，道長躬窜，這世上最難降的妖魔是什么？我笑而不...
開封第一講書人閱讀 58,963評論 1贊 295
?港島之戀（遺憾婚禮）
正文為了忘掉前任炕置，我火速辦了婚禮荣挨，結(jié)果婚禮上溜族，老公的妹妹穿的比我還像新娘。我一直安慰自己垦沉，他們只是感情好煌抒，可當(dāng)我...
茶點故事閱讀 67,984評論 6贊 393
惡毒庶女頂嫁案：這布局不是一般人想出來的
文/花漫我一把揭開白布。她就那樣靜靜地躺著厕倍，像睡著了一般寡壮。火紅的嫁衣襯著肌膚如雪。梳的紋絲不亂的頭發(fā)上讹弯，一...
開封第一講書人閱讀 51,763評論 1贊 307
城市分裂傳說
那天况既，我揣著相機與錄音，去河邊找鬼组民。笑死棒仍，一個胖子當(dāng)著我的面吹牛，可吹牛的內(nèi)容都是我干的臭胜。我是一名探鬼主播莫其，決...
沈念sama閱讀 40,468評論 3贊 420
雙鴛鴦連環(huán)套：你想象不到人心有多黑
文/蒼蘭香墨我猛地睜開眼，長吁一口氣：“原來是場噩夢啊……” “哼耸三！你這毒婦竟也來了乱陡？” 一聲冷哼從身側(cè)響起，我...
開封第一講書人閱讀 39,357評論 0贊 276
萬榮殺人案實錄
序言：老撾萬榮一對情侶失蹤仪壮，失蹤者是張志新（化名）和其女友劉穎憨颠，沒想到半個月后，有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體积锅，經(jīng)...
沈念sama閱讀 45,850評論 1贊 317
?護林員之死
正文獨居荒郊野嶺守林人離奇死亡爽彤，尸身上長有42處帶血的膿包…… 初始之章·張勛以下內(nèi)容為張勛視角年9月15日...
茶點故事閱讀 38,002評論 3贊 338
?白月光啟示錄
正文我和宋清朗相戀三年，在試婚紗的時候發(fā)現(xiàn)自己被綠了缚陷。大學(xué)時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片适篙。...
茶點故事閱讀 40,144評論 1贊 351
活死人
序言：一個原本活蹦亂跳的男人離奇死亡，死狀恐怖蹬跃，靈堂內(nèi)的尸體忽然破棺而出匙瘪，到底是詐尸還是另有隱情，我是刑警寧澤蝶缀，帶...
沈念sama閱讀 35,823評論 5贊 346
?日本核電站爆炸內(nèi)幕
正文年R本政府宣布，位于F島的核電站薄货，受9級特大地震影響翁都，放射性物質(zhì)發(fā)生泄漏。R本人自食惡果不足惜谅猾，卻給世界環(huán)境...
茶點故事閱讀 41,483評論 3贊 331
男人毒藥：我在死后第九天來索命
文/蒙蒙一柄慰、第九天我趴在偏房一處隱蔽的房頂上張望鳍悠。院中可真熱鬧，春花似錦坐搔、人聲如沸藏研。這莊子的主人今日做“春日...
開封第一講書人閱讀 32,026評論 0贊 22
一樁弒父案概行，背后竟有這般陰謀
文/蒼蘭香墨我抬頭看了看天上的太陽蠢挡。三九已至，卻和暖如春凳忙，著一層夾襖步出監(jiān)牢的瞬間业踏，已是汗流浹背。一陣腳步聲響...
開封第一講書人閱讀 33,150評論 1贊 272
情欲美人皮
我被黑心中介騙來泰國打工涧卵，沒想到剛下飛機就差點兒被人妖公主榨干…… 1. 我叫王不留勤家，地道東北人。一個月前我還...
沈念sama閱讀 48,415評論 3贊 373
代替公主和親
正文我出身青樓柳恐，卻偏偏與公主長得像伐脖，于是被迫代替她去往敵國和親。傳聞我的和親對象是個殘疾皇子乐设，可洞房花燭夜當(dāng)晚...
茶點故事閱讀 45,092評論 2贊 355