【RAY TRACING THE REST OF YOUR LIFE 超詳解】 光線追蹤 3-2 蒙特卡羅(二) 重要性採樣

 

書本內容：見相冊

注：鑑於不少網站隨意爬取數據，可能致使內容殘缺以及引用失效等問題，影響閱讀，請認準原創網址：

https://www.cnblogs.com/lv-anchoret/category/1368696.html

 preface

還記的咱們上一篇說的Monte Carlo 維度詛咒嗎

上一篇算是二維的例子吧，你們看了以後是否想着寫一個一維的Monte Carlo模擬積分？（我想了，沒寫出來）

爲何要整這個嘞

光照渲染中多處涉及重積分，最終結果是要求取一個近似值，於是須要對其值進行數值估計，Monte Carlo方法就是一個較爲理想的方案。

其實咱們的光線追蹤器不只用了不少向量運算，還用了不少數值分析的知識，好比以前的背景用的是最簡單的插值，柏林噪聲紋理用的是三線性插值等

 

 ready

你須要對如下知識有了解：

二重定積分

機率密度函數

反函數

 

 Chapter 2：One Dimensional MC Integration

積分是用來求取面積或者體積的東東，常見形式爲：I=∫x² dx

假設積分區間爲0->2,咱們用計算機符號表示爲I = area(x²,0,2)

那麼咱們來模擬如下這個積分

咱們先來看一個普通的：

 

那麼，咱們就能夠寫代碼了（程序2-1）

void MC_integration(double(*f)(double x), double low, double high)
{
    size_t all = 1000000;
    double sum{ 0 };
    for (int i = 0; i < all; ++i)
    {
        double x = (high - low) * rand01() + low;
        sum += x*x;
    }
    stds cout << "I = " << (high - low) * sum / all << stds endl;
}

int main()
{
    MC_integration([](double x) {return x*x; }, 0, 2);
}

正如咱們所指望的，但上述方法也能夠計算那些咱們用解析法沒法解決的積分，如被積函數爲log（sin（x）。

在圖形學中，咱們常常有一些咱們能夠評估但不能明確寫下來的函數，或者咱們只能機率評估的函數。 事實上，前面兩本書的光線追蹤代碼中的lerp函數，其實咱們不知道在每一個方向上看到的是什麼顏色，但咱們能夠在任何給定的維度上統計估算它。

 

我 又想起來上本書剛剛寫過的區域光照渲染圖片，裏面的黑色噪點真的把好精美的Cornell box 盒子給糟蹋了，那時由於咱們對全部的東西作了統一的採樣，而沒有對光源進行足夠的光源採樣，因此咱們能夠對光源進行更多的隨機採樣，可是咱們又須要控制這個比重，以防止過採樣的發生，因此，如何控制這個比例，已達到更好的效果，這就須要引入一個很是重要的概念——重要性採樣

可是在講重要性採樣以前，咱們須要瞭解一些關於機率密度函數的東東

　　

 

先看書上這張圖 

若是咱們添加更多的tree，那麼這張圖的長條將會增高（變長），由於他們的數量更多了，若是Height軸細分爲更多的小塊，那麼，每一個長條的高度會下降（變短）。

離散密度函數與直方圖的不一樣之處在於它將頻率y軸歸一化爲分數或百分比。

而連續直方圖中，咱們將Height軸細分爲無數個小塊，那麼，縱軸將不會是一個分數，由於全部小格對應的長條基本沒有高度（每一個小格表明的樹高區間基本只有一棵樹或者沒有樹，出現頻率幾乎爲0）

因此對於密度函數，咱們調整小格對應的影響因素，使得咱們添加更多的小格時，長條不會過低

對於上述條形圖，咱們採起以下：

小格的高 = 樹高介於H~H'的比例 / (H - H')

它將是頗有用的！ 咱們能夠將其解釋爲樹高的統計預測器：

即：一個介於H~H'的隨機樹高的機率 = 小格的高 * (H - H')

若是咱們須要知道多個小格對應的機率，那麼加和便可

其實上述公式就構成了機率密度函數，即，面積表明着對應橫座標區間的機率，也就是上面的第二個公式

而機率密度函數就是一種連續的分數直方圖，簡稱爲pdf(probability density function)

 

讓咱們來整一個pdf，或許會有更深刻的理解。若是我想要一個0~2的隨機數r它的出現機率和r本身成比例，咱們將指望pdf爲以下圖像：

 

 

r位於（x0,x1）的機率爲C*area(p(r),x0,x1),其中C爲常量,爲了簡便，C通常取1

而area(p(r),0,2) = 1（機率密度函數的總面積就是總機率）

由於p(r)和r成比例，因此，設另外一個常量C',則p(r) = C'r

故咱們解積分：

1 = ∫_0->2 C'r dr = C'r²/2|_0->2 = 2C' =>C' = 1/2  公式2-1

因此 p(r) = r/2

咱們如何用pdf p(r)生成一個隨機數？爲此咱們須要更多的機器。不要擔憂這不會永遠持續下去！（不會像上一篇最後那個程序，永無止境。。）

 

給一個隨機數d = rand01(),咱們知道咱們寫的rand01()產生的隨機數對於0~1的每一個浮點數的機率都是均勻一致的，咱們須要找個函數f(d)獲得咱們想要的。假定e = f(d) = d²,此次將再也不是一個均勻一致的pdf了，之前是的，以後咱們講（把此處記爲pos1遺留問題）

咱們爲了觀察這個函數還須要累積機率分佈函數Cpdf(Cumulative probability distribution function)

記爲P(x),則P(x) = area(p,-inf,x),幾何意義就是從負無窮到x的累積機率和，也就是x左邊部分的積分造成的面積

並且咱們規定在沒有定義的定義域內p(x) = 0,若採用pdf函數p(r) = r/2，那麼，P(x)以下：

P(x) = 0,　　x < 0

P(x) = x²/4, 0 < x < 2

P(x) = 1,　　x > 2

當x∈（0,2）,P(x) = ∫_0->x p(r) dr = ∫_0->x r/2 dr = x²/4

那麼，x和r是什麼關係？其實就是兩個虛擬變量，相似於程序中的參數，若是咱們把x調整爲有效區間的中點，那麼

P(1.0) = 1/4

這就是說基於咱們的pdf機率密度函數p(r) = r/2的設定下產生一個小於1的隨機數的機率爲25%

這產生了一種巧妙的觀察結果，這種觀察結果是產生非均勻隨機數的許多方法的基礎。

咱們想要一個函數f（），它知足：

當咱們調用f（drand48（））時，咱們獲得一個基於設定的Cpdf（x²/4）的一個返回值。 咱們不知道那是什麼，但咱們知道當x爲25％時應該小於1，而參數爲75％應該高於1。 若是f（）爲增函數，那麼咱們指望f（0.25）= 1.0。 

上面說了一堆，意思就是x和r無非就是函數參數，即P(1.0) = 0.25的時候，咱們一樣想獲得一個函數f()，它知足f(0.25)=1.0

即：根據Cpdf得知隨機數小於1.0的機率爲0.25，而存在一個f()，使得以0.25爲參數的時候可以獲得那個1.0（即累積機率分佈函數右端的邊界線）

這能夠推廣到每一個可能的輸入

即：f(P(x)) = x

也就是說f(x) = P^-1(x)

也就是P(x)的反函數。對於咱們的目的而言，上述的意思就是：若是咱們有了pdf函數p()和對應的Cpdf函數P(),咱們爲f()輸入一個隨機數，那麼它就會給出咱們想要的：

e = P^-1(rand01())

對於咱們的pdf函數p(x) = x/2，以及對應的P(x),咱們須要計算P^-1

即：若是咱們有 y = x²/4

那麼，x = sqrt(4y)

所以基於pdf,咱們輸入隨機數，獲得的就是

e = sqrt(4*rand01())

那麼e∈（0,2），正如咱們所指望的，咱們輸入一個0.25，它就返回一個1.0

 

到了這裏，其實上面重複了不少廢話，上面一大堆，表面上幹了一個什麼事情呢，用一個0~1的隨機數映射0~2這個積分區間

咱們想一下咱們要作的是什麼，模擬x的二次方曲線的積分，咱們要用隨機數模擬隨機選取積分區間的值，而咱們的積分區間就是0~2，

可是不必這麼麻煩啊，0~1隨機數隨機模擬0~2取值，直接2 * rand01()不就完啦，也就是程序2-1

 

不不不，它們有着天差地別，咱們上述一頓操做是爲了引入一個概念——Important Sample

咱們經過設置pdf使得，對於積分區間的自變量取值有了權重比例，也就是說有些地方的採樣要多一點，有些地方要少一點。

也就是：I = 1/n * Σ(F(x_i)/pdf(x_i))

由此可知，每一個積分採樣 I_i = 採樣高度（函數值）/該採樣點佔的權重，全部的積分採樣取均值，實現了不一樣採樣點爲整個積分作的貢獻是不一樣的

constexpr double pdf(const double x) {return 0.5 * x; }

void pdf1()
{
    size_t all{ 100000 };
    double sum{ 0. };
    for (int i = 0; i < all; ++i)
    {
        double x = sqrt(4 * rand01());
        sum += x*x / pdf(x);
    }
    stds cout << "I = " << sum / all << stds endl;
}

 

咱們來解釋pos1處的遺留問題

以前的程序2-1所指的積分方案中，是均勻一致採樣，即對於任意採樣點同等對待

即,設p(x) = C

則公式2-1 寫爲 

1 = ∫_0->2 C dr = Cr|_0->2 = 2C =>C = 1/2

即，開篇的方法中的pdf函數爲p(x) = 1/2

因此程序2-1能夠改寫爲：

constexpr double pdf(const double x) {return 0.5; }

void pdf1()
{
    size_t all{ 100000 };
    double sum{ 0. };
    for (int i = 0; i < all; ++i)
    {
        double x = 2 * rand01();
        sum += x*x / pdf(x);
    }
    stds cout << "I = " << sum / all << stds endl;
}

 

在important sample中咱們以爲函數峯高的部分對於整個積分的貢獻更多，低峯對於積分貢獻不大，因此咱們對於x²積分函數採起的pdf函數爲x/2

其實，pdf應該採起的就是積分函數自己，在知道答案的狀況下，咱們採起以下pdf函數最佳：

p(x) = (3/8)x²

則P(x) = x³/8

P^-1(x) = pow(8x,1/3)

則函數爲：

constexpr double pdf(const double x) { return 3 * x * x / 8; }

void pdf1()
{
    size_t all{ 1 };
    double sum{ 0. };
    for (int i = 0; i < all; ++i)
    {
        double x = pow(8 * rand01(), 1. / 3);
        sum += x*x / pdf(x);
    }
    stds cout << "I = " << sum / all << stds endl;
}

 

至此，全部的內容都已經闡述完畢

咱們回顧一下，由於這個是光線追蹤蒙特卡羅中的大部分概念

1. 你有一個被積函數f(x)的積分域【a,b】

2. 在域【a,b】中選擇一個非零的pdf函數p()

3. 對全部的積分採樣I_i = f(r)/p(r)取平均，pdf函數p()，隨機數r

這老是會收斂到正確答案，p越接近f，收斂越快

 

感謝您的閱讀，生活愉快~