1、第14章 蒙特卡羅積分I：基本概念在介紹計算沿光線到達相機的輻射亮度的SurfaceIntergrator類和VolumeIntegrator類之前，我們先做一些關于求解散射積分 方程的基本工作。這些積分方程通常沒有解析解，所以我們必須求助于數值方法。雖然諸如梯形積分法和高斯積分法這類標準的數值積分技術在求低維光滑積分時很 有效，但它們對于在渲染過程中所遇到的高維不連續積分而言，其收斂速度還是太慢了。蒙特卡羅積分技術為這個問題提供了一個解決方案。這類方法使用了隨機性做積分求值，且收斂速度跟被積函數的維的個數無關。在本章中，我們將回顧一些概率論 中的重要概念，并為使用蒙特卡羅解決渲染過程中的積分

2、問題打下一個基礎。下一章將描述提高收斂速度的計算，并描述pbrt所用到的技術。對隨機性的明智運用使算法設計領域發生了革命性的變化。隨機性算法大致分兩大類：拉斯維加斯算法和蒙特卡羅算法。拉斯維加斯算法使用了隨機性但最終會產生 相同的結果（如在快速排序算法中隨機地選擇一個數組項作為主元素）。而蒙特卡羅算法依賴于在算法過程中所使用的不同的隨機數列而得到不同的結果，但結果大 致上（在平均意義上）是正確的。所以，對蒙特卡羅算法的多次運行結果做平均，就有可能找到跟正確答案很接近的統計意義上的結果。蒙特卡羅積分是一項用隨機采樣來估算積分值的技術。蒙特卡羅積分的一個非常有用的特性是：為了估算積分值，我們只需保

3、證能夠在定義域中的任意點對被積函數求值。這個特性不僅使得蒙特卡羅技術易于實現，而且可以應用于非常廣泛的被積函數類型，包括那些不連續函數。在渲染過程中所遇到的積分是很難或不可能直接求積分的。例如，為了計算表面上某一個點上的反射光的量，我們要用方程(5.6)來對入射輻射亮度和BSDF 的乘積做單位球面上的積分。因為物體的可見性使得在復雜場景中的入射輻射亮度函數會有難以預測的變化，對于這個乘積中的所有項就很難找到封閉形式的表達 式，即使表達式存在，對之進行解析方式的積分也是通常不可能的。因為有了蒙特卡羅積分，計算反射輻射亮度就成為可能，因為我們只需選擇球面上的一組方向， 計算相應的入射輻射亮度，分別

4、乘以相應方向上的BSDF值，再用上一個權值項。這樣一來就可以處理任意的BSDF、光源和幾何體；所要做的只是在任意點上 對這些函數求值。蒙特卡羅積分的主要缺點是，如果用了n個采樣點來估算積分，算法收斂到正確解的速度是O(n-1/2)。換句話來說，如果將誤差減 半，就要使用4倍的采樣點。在渲染過程中，為了計算被積函數的值，每個采樣點需要追蹤一條或多條光線，當用蒙特卡羅積分做圖像合成時，這真是很痛苦的代 價。在圖像中，蒙特卡羅采樣所產生的人為缺陷會以噪聲的形式出現，即像素隨機性地出現太黑或太亮的情況。當前許多關于蒙特卡羅積分的研究著重于如何盡可能 地減少誤差又要減少采樣點的數量。14.4 背景和概率

5、論知識的回顧我們先定義一些基本的術語并回顧一下概率論的基本知識。我們先假定讀者已經熟悉了概率論的基本概念；需要對這一領域做近一步的全面了解的讀者可以參考相應 的教科書，例如Sheldon Ross的Introduction to Probability Models(2002)。一個隨機變量X是在某個隨機過程中所選定的一個值。我們通用用大寫字母來代表隨機變量，特殊情況下我們使用希臘字母代碼一些特殊的隨機變量。隨機變量總是 在某些定義域中抽取出來的，它既可以是離散的（例如一個固定的概率集合），也可以是連續的（如實數集R）。對隨機變量X使用函數f就得到另一個新的隨機變 量Y = f(X)。例如，隨

6、機擲一個骰子的結果就是一個離散的隨機變量，它被采樣于一個事件集合Xi = 1, 2, 3, 4, 5, 6 。每個事件有一個概率pi = 1/6， 這樣概率總和pi一定是1。我們可以取一個連續的均勻分布的隨機變量0,1，并將之映射為一個隨機變量， 如果下式滿足，則取值Xi：        j=1,i-1 Pj  < j=1,i Pj注意所有的Pi之和為1。一個隨機變量的累積分布函數（CDF）P(x)是變量分布中的一個值小于或等于某個值x的概率。       &

7、#160;P(x) = Pr X x以擲骰子為例，P(2) = 1/3， 因為小于或等于2的情況占總數6個中的2個。14.1.1 連續隨機變量在渲染過程中，跟離散隨機變量相比，會更多地用到連續隨機變量，這些變量在一個連續的定義域內取值（如實數域，或單位球上的方向）。其中一個非常重要的隨機變量是典型均勻隨機變量，記為。該變量之所以重要有兩個原因。首先，在軟件中生成這個分布的變量很容易，因為大多數運行庫中的偽 隨機數發生器就可以做到這一點。第二，我們可以先生成典型均勻隨機變量再施加上某種變換來得到任意的隨機分布。前面所提到的將映射到骰子中6個面就是這 個技術在離散情況中的應用。我們來看另一個定義在

8、0,2上的連續隨機變量的例子：它在值x的概率值跟2-x成正比：即在0附近取值的可能性是在1處附近取值時的兩倍。概率密度函數(PDF)對此有正規的表述：它描述了隨機變量在某個值上取值的相對概率。PDF p(x)是CDF的導數：        p(x) = dP(x) / dx均勻隨機變量的p(x)是個常數。對于，我們有：        p(x) = 1 （當x 0,1)  或 0 (其它情況）PDF必須是非負的，并在其定義域上的積分是1。給定一個定義域中的一

9、個區間a,b，PDF給出了隨機變量處于該區間的概率：        P (x a,b)  = a,b p(x) dx上式可以根據微積分第一基本定理和PDF的定義得出。14.1.2 期望值和方差函數f的期望值Epf(x)被定義成該函數在其定義域的上某種p(x)分布的平均值。在下一節中我們將會看到如何用蒙特卡羅積分計算任意積分的期望值的。在域D中的期望值定義如下：        Epf(x) = D f(x)p(x) dx舉個例子，考慮一下如何求余弦函數在0,

10、上的期望值，其中p是均勻的。因為PDF p(x)在定義域上的積分必須為1，所以 p (x) = 1 / 。故有：        Ecosx = 0,   p(x) dx = (-sin + sin 0) / = 0這正是所期望得到的值，考慮一下cosx在0,的曲線就明白了。一個函數的方差是函數值跟其期望值的誤差函數。方差是用來量化蒙特卡羅算法估算值的誤差的基本概念。它提供了誤差量化的精確方法，并可以測量蒙特卡羅算法 為減少誤差所做的改進程度。第15章的大部分內容將用來介紹減少方差從而改進pbrt的計算結果的技術。一

11、個函數f的方差定義為：         Vf(x) = E (f(x) - E f(x)2 從期望值和方差定義中，我們馬上可以看到它們有三個重要的性質。        E a f(x) = a E f(x)        E i f(Xi) = i  E f(Xi)        V a f(x) = a2 V f(x)利用這些性質做些簡

12、單的代數推導，可以得到方差的一個更簡單的表達式：        Vf(x) = E (f(x) 2  - E f(x) 2這樣，方差只是函數平方的期望值減去期望值的平方。如果隨機變量是相互獨立的，那么方差還有這樣的性質，即方差的和等于和的方差：        i Vf(Xi) = V i  f(Xi)14.2 Monte Carlo估計量現在我們可以定義基本的Monte Carlo估計量了，它可以近似地估算任意積分的值。這是第16，17

13、章中光傳輸算法的基礎。假定我們想要對一個一維積分a, b f(x) dx求值。對于一個均勻分布的隨機變量Xia,b，Monte Carlo估計量用下式給出期望值：        FN = (b - a) / N  1,n f(Xi)E(FN)就是這個積分值。這可以用幾步推導來得出。注意隨機變量Xi所對應的PDF p(x)必須等于1/(b-a)，因為p必須是常量，并且在a,b上的積分為1。推導如下：我們做一些一般化推廣，就可以去掉隨機變量必須是均勻的限制。這是極其重要的一步，因為對生成采樣的PDF的選取是減少蒙特卡羅

14、方差（第15.4節）的重要技術。如果隨機變量是某個任意PDF p(x)來抽取的，那么下面的估計量可用來估算積分：        FN = 1/ N 1,Nf(Xi) / p(Xi)p(x)的唯一的限制是它必須在所有 滿足F(x) > 0的x為非零。類似地，我們可以看出這個估計量的期望值計算所要求的積分：我們可以直接將這個估計量推廣到高維或復雜的積分域中。從多維PDF生成的N個采樣也可以照樣使用這個估計量。例如，考慮一個三維積分：        xo,x1 yo,y1&#

15、160; yo,y1  f(x,y,z) dxdydz如果采樣Xi = (xi, yi, zi)從三維盒(x0,y0,z0)-(x1,y1,z1)中均勻地取出，則PDF p(X)就是常量：        1/(x1-x0)(y1-y0)(z1-z0)估計量為：        (x1-x0)(y1-y0)(z1-z0) / N f(Xi)注意采樣個數N是任意選取的，跟被積函數的維數無關。這是另一個Monte Carlo對傳統積分的優勢。蒙特卡羅所使

16、用的采樣個數完全獨立于積分維數，而標準的數值積分計算所需要的采樣個數跟維數呈指數級增長。只指出蒙特卡羅估計量可以收斂于正確結果還不夠，因為良好的收斂速度也是很重要的。我們這里不做推導，只指出其誤差按照O(N1/2) 的速度收斂(N為采樣個數)。雖然標準積分技術在一維的情況下的收斂速度要高于這個速度，但隨著被積函數的維數增加，其性能會有指數級的下降，而蒙特卡羅 的收斂速度跟維數無關，這樣蒙特卡羅就成為唯一可行的高維積分的數值積分技術。我們已經遇到一些高維積分的例子，在第16章的路徑跟蹤公式實際上是個無限 維的積分！14.3 隨機變量的采樣為了對蒙特卡羅估計量求值，必須能夠從所選定的概率分布中抽取

17、隨機采樣。本節介紹這個過程的基本知識，并舉一些簡單的例子。下一節將介紹一般性的多維情況，下一章的第15.5節和15.6節將利用這些技術在BSDF和光源的分布中生成采樣。14.3.1 逆轉法逆轉法使用一個或多個均勻隨機變量，并將之映射到所需分布的隨機變量中。為了解釋這個過程，我們先看一個簡單的離散的例子我們有一個4種結果的過程，每個結果的概率分別為p1,p2,p3,p4，并且p1+p2+p3+p4=1。相應的PDF如下：    為了在這個分布中取樣，我們先求CDF P(x)。在連續的情況下，P是p的不定積分。在離散的情況下，我們可以直接從左邊開始，將相應的直方條疊在一起，如圖

18、：（注意最右的一疊的高度為1，因為所有的概率之和為1）。為了在分布中抽樣，我們使用一個均勻隨機數，用它通過CDF來選擇其中一個可能的結果，并且要求該結果的概率等于它自己的概率。如圖，事件的概率被投影 到垂直軸上，隨機變量從中選取概率值。應該清楚這個抽樣的分布是正確的，即均勻采樣的概率恰好等于它所碰到的特定直方條的高度。為了將這個技術推廣到一 般的連續分布上，我們可以考慮當離散的概率個數趨于無窮大的情況。則上面的PDF變成一個光滑的曲線，而CDF就是它的積分。上面所表示的投影過程仍然一 樣，這個投影有一個很方便的數學解釋，即它表示求CDF的反函數，并在處求反函數的值。這個技術被稱為逆轉法。更精確

19、地講，我們可以從一個任意的PDF p(x)中按照下列步驟抽取采樣Xi:1. 計算CDF P(x) = 0, x f(x') dx'。2. 求反函數P-1(x)。3. 取一個均勻分布的隨機數。4. 計算Xi = P-1()。14.3.2舉例：冪分布作為一個該過程的一個例子，我們考慮從一個冪分布抽樣的過程。我們在對Blinn微平面模型采樣時就是這種情況。冪分布的PDF為：        p(x) = c x n第一個任務是求PDF。在大多數情況下，這值涉及到求比例常數c。我們可以利用p(x)dx = 1就可以得到：所以，p

20、(x) = (n+1)xn。我們對之積分，就可以得到CDF:    P(x) = 0, x  f(x') dx' = x n+1求反函數很簡單： P-1(x) = x (1/n+1)。這樣，給定一個均勻隨機變量，從冪分布抽樣所得到的結果是：        X = (1/n+1)14.3.3 舉例：指數分布當我們渲染有參與介質的圖像時，經常要從一個指數分布中抽樣。第一步是對該分布進行歸一化，使其積分值為1。在這種情況下，我們希望所生成的采樣能夠覆蓋0, infnite)，而不是0,

21、1，故有：所以c = a，而PDF為p(x) = ae-ax。對之積分得到P(x):         該函數的反函數為：        P-1(x) = - ln(1-x) / a故我們可以得到下列的抽樣：        X = - ln(1- ) / a我們可以對該式做近一步地簡化，因為是均勻分布的隨機數，故1- 也是均勻分布的，故而我們可以將1- 替換成而不影響分布，故有：    

22、   X = - ln() / a14.3.4 舉例：一維分段函數另一個很有趣味的練習是在一個一維分段函數上取樣。不失一般性，我們只考慮定義在0,1上的分段函數假定該函數的定義域被分為N等分，每等分為=1/N。這些區間的端點為xi = i ，其中i的范圍是0到N。在每個區間中，f(x)為一個常量。如圖：            f(x) = vi ,     其中 xi <=x < xi+1 (i = 0, N - 1)積分f(x)dx

23、為：    c = 0,1f(x)dx = i=1,N-1 vi / N        我們很容易構造出f(x)的PDF p(x)， 即f(x)/c。CDF P(x)是分段的線性函數（上面圖b）。    在點xi和xi+1之間，CDF按照斜率vi/c線性增加回憶一下為了對f(x)采樣，我們需要求CDF的翻函數，并找到值x，使得：         = 0,xf(x')dx' = P(x)由于CDF是單調遞增的，故x必在xi

24、到xi+1之間，并滿足P(xi) <= <= P(xi+1)為了有效地進行采樣，我們提供一個函數，它從f(x)的vi中取值，計算在xi的CDF值，并返回f(x)的積分值。<MC Function Definitions> =    void ComputeStep1dCDF(float *f, int nSteps, float *c, float *cdf)         <Compute integral of step function at xi>  

25、0;     <Transform step function integral into cdf>    該函數先計算f(x)的積分。它將結果存放到cdf數組，故無需再申請額外的內存。調用者必須為cdf數組申請nStep+1個浮點數的內存空間，因為如 果f(x)有N個階梯值，則需要存放CDF在每個N+1個xi處的CDF值。在數組后面存放CDF值1是冗余的，但有利于簡化代碼。<Compute integral of step function at Xi>    int i;  &

26、#160; cdf0 = 0.;    for (i = 1; i < nSteps+1; +i)        cdf i = cdf i - 1 + f i - 1 / nSteps;由于在0,1上的積分值存放到在cdfnSteps中，CDF可以用這個值進行歸一化：<Transform step function integral into cdf> =    *c = cdfnSteps;        for

27、(i = 1; i < nSteps+1; +i)            cdf i /= *c;下列函數用來對函數采樣：    <MC Function Definitions> +=    float SampleStepld(float *f, float *cdf, float c,            int nSteps, float u, float *pdf)

28、        <Find surrounding cdf segments>        <Return offset along current cdf segment>    首先我們需要找到覆蓋的一對CDF值。因為cdf數組是單調遞增的，我們可以用二分查找來得到它。其中用到c+標準函數lower_bound()。<Find surrounding cdf segments> =    flo

29、at *ptr = std:lower_bound(cdf, cdf+nSteps+1, u);    int offset = (int) (ptr-cdf-1);有了這一對CDF值，我們就可以計算x值了。我們確定在cdfoffset和cdfoffset+1中的位置。因為CDF是線性的，x位于xi 和xi+1之間。這個采樣的PDF p(x)很容易計算，因為我們有歸一化值c，它滿足p(x) = f(x) / c。<Return offset along current cdf segment> =    u = (u - cdfoffset)

30、 / (cdfoffset+1 - cdfoffset);    *pdf = foffset / c;    return (offset + u) / nSteps;14.3.5 拒絕法對于某些函數，我們無法使用積分的方法求它們的PDF，或者無法解析地得到CDF的反函數。拒絕法(rejection method)可以不做這兩步也可以按照函數的分布來生成采樣，這實際上是擲飛鏢方法。假定我們要在函數f(x)上采樣，但我們有一個PDF p(x)滿足f(x) < cp(x)，c為某個常數，并且假定我們知道如何從p上采樣。這時拒絕法就很簡單： 

31、   loop forever:        sample X from p's direction        if < f(X) / (cp(X) then            return X這個例程重復不斷地選擇一對隨機變量(X, )。如果點(X, cp(X)位于f(X)之下，則接受采樣X，否則就拒絕它，并選擇下一對采樣。無需深入討論，我們就可以明白該方法的效

32、率依賴于cp(x)界定f(x)的緊密程度。該技術適用于任何維數。在實際應用中， pbrt中的任何蒙特卡羅算法都沒有使用拒絕法。我們通常需要找跟f(x)類似的可以直接進行采樣的分布，下一章我們再介紹其中的原因。然而，拒絕法是一 項值得關注的重要技術，特別是在調試蒙特卡羅的程序時。例如，如果我們懷疑使用逆轉法抽樣程序有錯誤，我們可以用更直接的基于拒絕法的實現代替之，再看蒙 特卡羅估計量是否計算出相同的結果。當然，這要進行大量的采樣。14.3.6 舉例：用拒絕法采樣單位圓假定我們想要在一個單位圓內選一個均勻分布的點。如果使用拒絕法，只需隨機地在外切正方形內選擇一個隨機位置(x,y)，再看它是否落在圓

33、內。如圖：    下面的RejectionSampleDisk()函數實現了這個算法。我們可以用類似的方法來生成任意復雜形狀之內的均勻分布采樣，只要該形狀有一個測試點在形狀之內或之外的方法。<Monte Carlo Function Definitions> +=    void RejectionSampleDisk(float *x, float *y)             float sx, sy;      &#

34、160; do             sx = 1.f - 2.f * RandomFloat();            sy = 1.f - 2.f * RandomFloat();         while (sx*sx + sy*sy > 1.f);        *x = sx; 

35、;       *y = sy;    一般地說，拒絕法的效率依賴于被拒絕的采樣所占的比例。對于2D的均勻采樣的情況，還是很容易計算的，即是圓面積除以外切正方形面積：/4，約等于78.5%。對于一般的n維情況，在超球內采樣的效率會隨著n值的增大而降低（因為n維的超球體積趨于0）。14.4 分布之間的變換為了描述逆轉法，我們介紹了一種按某種分布生成采樣的技術，它將典型均勻分布的隨機變量以特定的方式進行變換。這里我們將研究將采樣從任意一種分布變換到另一種分布的更一般化的問題。假定我們有一個從某個PDF Px(x)抽樣的隨機變量X

36、i。現在，如果我們要計算Yi = y(Xi)，我們想要找到新隨機變量Yi的分布。這看上去是個很奇怪的問題，但我們將體會到理解這種變換對從多維分布函數中抽樣至關重要。函數y(x)必須是一一映射。如果多個x值被映射到同一個y值，那么就無法清楚地描述一個y值的概率密度。由于y必須是一一映射，那么它的導數必須嚴格大于0或嚴格小于0，這就意味著下式成立：        Pr y <= y(x) = Pr X <= x所以：        Py(y) = Py(y(x) =

37、Px(x)利用這個CDF之間的關系可以直接導出它們PDF之間的關系。如果我們假定y的導數大于0，通過求微分，可知：        py(y) dy/dx = px(x)故有：        py(y) = (dy/dx) -1 px(x)一般而言，y的導數要么嚴格為正，要么嚴格為負，它們的密度關系式如下：        py(y) = |dy/dx| -1 px(x)怎么用這個公式呢？假定px(x) = 2x，定義域

38、為0,1，令Y = sinX。隨機變量Y的PDF是什么呢？因為我們知道dy/dx = cosx，        py(x) = px(x) / |cosx| = 2x / cosx = 2sin-1y/(1-y2)1/2)這個例程似乎有些倒退-通常我們已知某個采樣的PDF，而不是某個變換。例如，我們想從某個px(x)抽樣的X，想要從某個分布py(y)計算出Y。那么我們用什么變換呢？我們只需要求它們的CDF相等，即Py(y) = Px(x)， 那么就可以得到變換：        

39、;y(x) = Py-1(Px(x)這是逆轉法的一般化情況，因為如果X是在0,1均勻分布的，就有Px(x)=x，就會得到上一節的相同的例程。14.4.1 多維情況下的變換對于一般的n維情況，我們可以有類似的推導來得到不同密度之間的類似關系。這里不再列出推導過程了，因為它跟一維情況有相同的形式。假定我們有一個n維隨機變量X，其密度函數為px(x)。另Y = T(X)， 其中T為一個雙射函數。這時，密度關系式為：        py(y) = py(T(x) = px(x) / |JT(x)|其中|JT|是T的雅克比矩陣的行列式值的絕對值

40、。14.4.2 舉例： 極坐標極坐標變換公式為：        x = r cos        y = r sin假定我們要從每個密度函數p(r, )抽樣，相應的密度函數p(x,y)應該是怎樣的？這個變換的雅克比矩陣為：        它的行列式值為r(cos2 + sin2) = r。所以p(x,y) = p(r, ) / r。當然并非我們想要得到的結果-因為我們先從笛卡爾坐標系開始采樣再將之變換到極坐標，這時我們有p(r, ) =

41、 r p(x,y)。14.4.3 舉例：球面坐標球面坐標變換公式為：        x = r cos cos        y = r sinsin        z = r cos這個變換的雅克比行列式值為| JT| = r2sin，所以相應的密度函數為：        p(r, , ) = r2sin p(x,y,z)這個變換很重要，因為它可以將方向用

42、單位球上的點(x,y,z)來表示。我們記得立體角是用單位球面上的點集面積來定義的。在球面坐標系中，我們在前面推導過：        d = sin d d所以，如果我們有一個在立體角上的密度函數，這意味著下式成立：        Pr    = p()dw我們就可以推導出關于 和的密度函數:        p(, ) d d = p() d       &#

43、160; p(, ) = sin p()14.5 使用多維變換的2D采樣假定我們有一個2D聯合密度函數p(x,y)做(X,Y)的采樣。有時多維函數是可分離的，可以被表示為多個一維密度函數的乘積，例如：        p(x,y) = px(x) . py(y)在這種情況下，我們可以獨立地從px采樣X，從py中采樣Y來得到(X,Y)。但是，許多有用的密度函數并不是可分離的，所以我們要介紹如何在一般情況下進行多維分布的采樣理論。給定一個2D密度函數，我們可以“積分掉”兩維中的一個而得到邊緣密度函數:   &#

44、160;    p(x) = p(x,y)dy這個函數可以視為X單獨的密度函數。更精確的說法是，它是特定的x對所有可能的y值的平均密度。條件密度函數p(y|x)是選定某個x的關于y的密度函數：        p(y|x) = p(x,y) / p(x)從聯合分布中進行2D采樣的基本思想是，先計算出邊緣密度函數，將一個變量分離出來，再用標準的一維技術利用該密度函數采樣。一旦得到這個采樣，再計算這個采樣值所對應的條件密度函數，然后再利用標準的一維技術繼續采樣。14.5.1 舉例：在半球上的均勻采樣做為一個例子

45、，我們在半球上根據立體角均勻地對方向采樣。我們還記得均勻分布的意義是密度函數是一個常量，所以我們有p() = c。由于密度函數在其定義域上積分必須為1，我們就得到：    故而p() = 1/(2)，或者p(,) = sin/（2）。注意這個密度函數是可分離的。然而，我們仍使用邊緣密度和條件密度來演示多維采樣技術。我們先來采樣。為此，我們需要知道的密度函數p():    p() = 0,2 p(,)d = 0,2sin d = sin現在我們計算的條件密度：        p(|) = p(,)

46、/p() = 1/(2)注意到的密度函數本身是均勻的，這可以從半球的對稱性上直觀地看出來。現在我們用一維逆轉技術分別對這些PDF進行采樣：    P() = 0, p(')d' = 1 - cos    P(|) =0, 1/2 d ' = /2 我們可以很方便地得到它們的反函數，注意到和1-在0,1都是均勻分布的隨機變量，我們可以用代替1-，就得到：         = cos -1 1         = 2

47、 2將這些值變換會笛卡爾坐標系中，就得到最后的采樣公式：         我們用下面的代碼來實現這個采樣策略。我們傳入兩個均勻分布的隨機數u1，u2，得到一個半球上的向量：<MC Function Definitions>    COREDLL Vector UniformSampleHemisphere(float u, float u2)         float z = u;       

48、; float r = sqrtf(max(0.f, .f - z*z);        float phi = 2 * M_PI * u2;        float x = r * cosf(phi);        float y = r * sinf(phi);        return Vector(x, y, z);  &#

49、160; 在pbrt中，對于每個類似這樣的采樣例程，都有一個對應的返回特定采樣的PDF值的函數。對于這樣的函數，我們需要搞清楚到底是用哪一個PDF求值 -例如，在對半球上方向上采樣時，我們已經看到這些密度函數即可以用立體角表示，也可以用（，）表示。對半球（包括所有其它關于方向的采樣）而 言，這些函數返回關于立體角的值，立體角的PDF是一個常數 p() = 1/(2)。<MC Function Definitions> +=    COREDLL float UniformHemispherePdf(float theta, float phi)   &

50、#160;     return INV TWOPI;    在整個球面上的均勻采樣公式可以用相同的方式推導出來，這里略去了。公式如下：<MC Function Definitions> +=    COREDLL Vector UniformSampleSphere(float u1 , float u2)        float z = 1.f - 2.f * u1;       

51、0;float r = sqrtf(max(0.f, 1.f - z*z);        float phi = 2.f * M_PI * u2;        float x = r * cosf(phi);        float y = r * sinf(phi);        return Vector(x, y, z);  

52、0; <MC Function Definitions>+=    COREDLL float UniformSpherePdf()         return 1.f / (4.f * M_PI);    14.5.2 舉例：在單位圓盤上的采樣雖然圓盤看上去比半球更簡單，但要進行均勻采樣卻更麻煩一些。一個直觀但不正確的解決方案是： r = 1,   = 2 2。雖然所得到的結果是隨機的且落在圓盤之內，但它并不是均勻分布的，而是采樣點擠在靠近圓心的位置。因為我們

53、是根據面積來進行均勻采樣，PDF p(x,y)必須是個常數。通過歸一化約束我們得到 p(x,y) = 1 / 。如果我們變換到極坐標系中，則有p(r, ) = r / 。現在我們使用前面介紹的方法計算邊緣密度和條件密度：    跟半球上采樣一樣，由于圓的對稱性，我們知道p( | r)是一個常數。通過積分和求反函數，就得到了P(r), P-1(r), P() 和p-1()。我們可以求得在圓盤上均勻采樣的正確解：        r = 11/2        

54、 = 2 2<Monte Carlo Function Definitions> +=    void UniformSampleDisk(float u1, float u2, float *x, float *y)         float r = sqrtf(u1);        float theta = 2.0f * M_PI * u2;        *x = r *

55、cosf(theta);        *y = r * sinf(theta);    雖然解決了目前的問題，但圓盤上的面積卻產生了變形（如圖），第15.2.3節再介紹這種變形所帶來的壞處。Peter Shirley(1997)開發出一種從單位正方形到單位圓的“同心”映射來避免這個問題（如圖）。    同心映射將正方形-1,1x-1,1上的點均勻地按照“同心正方形-同心圓”的方式進行映射，如圖：正方形上的三角楔形被映射為餅條形狀，映射公式為：      

56、  r = x         = y / x<Monte Carlo Function Definitions> +=    void ConcentricSampleDisk(float u1, float u2, float *dx, float *dy)         float r, theta;        <Map uniform rando

57、m numbers to 1, 1x-1,1>        <Map square to (r , )>        *dx = r * cosf(theta);        *dy = r * sinf(theta);    <Map uniform random numbers to 1, 1x-1,1>    float sx = 2

58、* u1 - 1;    float sy = 2 * u2 - 1;<Map square to (r , )> =    <Handle degeneracy at the origin>    if (sx >= -sy)         if (sx > sy)             <Handle first region of disk>

59、;                else             <Handle second region of disk>                else         if (sx <= sy) 

60、0;           <Handle third region of disk>                else             <Handle fourth region of disk>                theta *= M_PI / 4.f;<Handl