低差異序列 (low-discrepancy sequences)之Hammerysley在半球中取樣點方法的介紹

半球上的Hammersley

源作者：Holger Dammertz

一組關於如何在2D中使用Hammersley點集以在著色器程式中快速實用地生成半球方向的筆記。如果你發現任何錯誤或有意見，不要猶豫，請聯絡我或在我的部落格上留言。

1. 概述

在編寫與光照有關的著色器時，人們經常需要一些以表面法線為方向的半球方向。通常的做法是預先計算這些方向，並將它們儲存在一個靜態/統一的陣列中，或者建立一個包含這些方向的可查詢的紋理。在這一頁，我研究瞭如何在著色器中直接使用2d中的Hammersley點集來快速計算合理的分佈方向。當然，這個點集並不侷限於取樣方向，它還可以用於陰影圖過濾、螢幕空間環境遮蔽（SSAO）或其他任何需要2D取樣模式的應用。

二維的Hammersley點集是最簡單的低差異序列之一，在計算機圖形中的偽蒙特卡洛（Quasi Monte Carlo）整合方面有許多可能的應用。對各種不同的偽蒙特卡羅序列的總體概述可以在Niederreiter [niederreiter92]的書中找到。例如在 [wong97]中可以找到明確使用Hammerysley點集進行球體取樣的方法。

下面我首先簡要介紹一下如何從數學上構造Hammersley點集，然後展示這些點在2d中的樣子，以及它們在對映到半球時的樣子。之後，我向你展示了一個可以直接插入GLSL著色器中的實現（或者很容易翻譯成任何其他語言）。

2. Hammersley點集

我們考慮一個點集

與點的數量\(N≥1\)在二維單位面積上\([0,1)\)（原本中關於此處均使用左閉右開）。Hammersley點集\(H_N\)，在2D中，現在被定義為

其中\(\varPhi _2(i)\)是Van der Corput序列

Van der Corput序列的想法是在小數點處映象\(i\)的二進位制表示，以得到區間[0,1]中的一個數字。

Van der Corput序列\(\varPhi _2\)的數學定義如下:

其中\(a_0a_1...a_r\)是\(i\)用二進位制表示的各個數字(i.e. \(i=a_0+a_12+a_{2}2^{2}+a_{3}2^{3}+...a_{r}2^r\)).

2.1 Radical Inverse Examples

這一點最好通過一個例子來理解。在下面的表格中，我給出了\(i\)的十進位制和二進位制表示，以及相關的Van der Corput 基數逆（radical inverse）（如公式2中的定義），在小數點處二進位制表示。為了便於理解，我還在\(i\)的二進位制表示中加入了小數點，儘管在考慮二進位制整數時通常會省略這一部分。

該表顯示了Van der Corput序列中的前幾個數字。它只是i在小數點處映象的二進位制表示。

這種簡單的映象二進位制表示法使得Hammersley點集稱為在著色器中快速和容易實現的點集的理想人選。在我討論程式碼之前，我想向你展示點集在2D中的實際情況。

2.2 Hammersley點集影象

在這個互動式影象中，你可以改變點的數量\(N\)，並得到由公式1定義的、在單位方格\([0,1)^2\)中視覺化的結果點集。你可以看到，點\(\left( \begin{array}{c}
0\\
0\\
\end{array} \right)\)總是這個序列的一部分。你還可以看到，對於任何數量的N個點來說，這些點的分佈都很好

3. 生成半球上的點

為了從\(H_N\)建立半球上的方向，可以使用從\([0,1)\)到半球的任何對映，但必須注意所選擇的對映要滿足應用的要求。在光傳輸背景下進行積分時，通常需要球面上的均勻或餘弦加權分佈。

讓\(x_i=
\left( \begin{array}{c}
u\\
v\\
\end{array} \right)
∈H_N\)是我們選擇的Hammersley點集的一個點。讓兩個座標\(u\)和\(v\)，我們現在可以計算處以下的對映。把球座標系轉成直接座標系即可。

3.1 統一對映

3.2 餘弦對映

4. 原始碼

為了實現Hammersley點集，我們只需要一個有效的方法來實現Van der Corput基數逆（radical inverse）\(Φ_2(i)\)。由於它是以2為基數的，我們可以使用一些基本的位操作來實現。Hacker's Delight [warren01]這本書為我們提供了一個簡單的方法來反轉一個給定的32位整數的位。使用這個方法，下面的程式碼就可以在GLSL中實現\(Φ_2(i)\)

 float radicalInverse_VdC(uint bits) {

     bits = (bits << 16u) | (bits >> 16u);

     bits = ((bits & 0x55555555u) << 1u) | ((bits & 0xAAAAAAAAu) >> 1u);

     bits = ((bits & 0x33333333u) << 2u) | ((bits & 0xCCCCCCCCu) >> 2u);

     bits = ((bits & 0x0F0F0F0Fu) << 4u) | ((bits & 0xF0F0F0F0u) >> 4u);

     bits = ((bits & 0x00FF00FFu) << 8u) | ((bits & 0xFF00FF00u) >> 8u);

     return float(bits) * 2.3283064365386963e-10; // / 0x100000000

 }

然後通過以下方式計算出第\(i\)個點\(x_i\)

 vec2 hammersley2d(uint i, uint N) {

     return vec2(float(i) / float(N), radicalInverse_VdC(i));

 }

一個小的優化是將上面的程式碼中除以N的部分改為乘以它的逆數，因為無論如何N都需要被固定。

為了完整起見，這裡還有從Hammersley點集\({x_i=
\left( \begin{array}{c}
u\\
v\\
\end{array} \right)
}\)建立均勻或餘弦分佈方向(z-up)的程式碼。請注意，程式碼中的\(1-u\)是必要的，以便將序列中的第一個點對映到向上的方向而不是切平面。

 const float PI = 3.14159265358979;

 vec3 hemisphereSample_uniform(float u, float v) {

     float phi = v * 2.0 * PI;

     float cosTheta = 1.0 - u;

     float sinTheta = sqrt(1.0 - cosTheta * cosTheta);

     return vec3(cos(phi) * sinTheta, sin(phi) * sinTheta, cosTheta);

 }

 vec3 hemisphereSample_cos(float u, float v) {

     float phi = v * 2.0 * PI;

     float cosTheta = sqrt(1.0 - u);

     float sinTheta = sqrt(1.0 - cosTheta * cosTheta);

     return vec3(cos(phi) * sinTheta, sin(phi) * sinTheta, cosTheta);

 }

請注意，返回的方向當然是在預設的座標框架中，並帶有z-up（如果你願意，你也可以稱之為切線空間）。因此，在應用程式中可能需要執行的最後一步是通過使用由本地法向量構建的正交法向量基矩陣將返回的方向轉換為世界空間。

5. 補充說明

在這裡，我收集了一些關於如何使用這些點的進一步說明，以及其他一些值得了解的事情。

5.1. 這些點有多好呢

質量當然取決於應用，但首先要注意的是，對於少量的點，2D中的Hammersley點集在半球上的分佈不是很好。這一點可以在一節中很容易得到驗證。當需要很多點的時候，這個問題就會減輕，例如在進行交錯取樣的時候。雖然有許多其他的準隨機序列可以使用，但到目前為止，我還沒有發現任何可以像Hammeysley點集那樣有效實施的序列（除了秩-1的格子，但他們不容易用於不同數量的點）。

這種取樣模式的另一個問題是，在開始計算任何一個點之前，需要預先知道\(N\)個點的數量。這在實時渲染應用中幾乎不是問題，但一旦進行某種形式的漸進式渲染，使用一個準蒙特卡洛序列（而不是一個點集）可能是一個更好的選擇。

5.2. 交錯式取樣

交錯取樣[keller01]是一個經常使用的優化方法 [segovia06] ，用於實時渲染，其中一個積分的計算被分散到幾個畫素上，然後在一個額外的過濾步驟中合併。

當然，任何點的集合都可以用於此，但是簡單的構造和使用使它特別容易將點集合分割成幾個片段和/或半球。

6. 補充原文

原文中雖然有視覺化生成Hammerslep隨機點的例子，但是後續解釋公式的時候，並未給出詳細的例子，所以此處參考低差異序列（一）- 常見序列的定義及性質補充更為詳細的例子

6.1. Radical Inversion與Van der Corput序列

接下來在介紹這些序列定義之前，先介紹一個基本的運算，Radical Inversion。這篇專欄將會介紹的所有序列都會用到這個運算過程。

這個操作非常直觀，\(b\)時一個正整數，則任何一個整數\(i\)，如果先將\(i\)表示成\(b\)進位制的數，然後把得到的數中的每一個位上的數字\(a_l(i)\)排成一個向量，和一個生成矩陣（genetator matrix）\(C\)相乘得到一個新的向量，最後再把這個新向量中映象到小數點右邊去就能得到這個數以\(b\)為底數，以\(C\)為生成矩陣的radical inversion（基數逆）\(\varPhi_{b,C}(i)\)

上面的描述可能略微有些複雜，但如果矩陣\(C\)是單位矩陣（Identity Matrix）的時候，這個過程會簡化很多，即是直接把這個\(b\)進位制數映象到小數點右邊去即可。同事也是Van der Corput序列的定義。

舉個例子，正整數8以2為底數的radical inverse的計算過程如下，首先算出8的2進製表示，1000。此處假設\(C\)為單位矩陣，所以直接將1000映象到小數點右邊，0.0001.這個二進位制數的值就是最終結果，把它轉換回10進位制就得到1/16，既\(\varPhi_2(8)=\frac{1}{16}=0.0625\)。下面的表給出了更多的以2為底數的Van der Corput序列的例子。

Van der Corput序列有幾個屬性：

• 每一個樣本點都會落在當前已經有的點裡"最沒有被覆蓋"的區域。例如\(\varPhi_2(3)=\frac{3}{4}\)時剛好落在了\([0,1)\)區間中被覆蓋最少的區域（\(\varPhi_2(0)=0,\varPhi_2(1)=\frac{1}{2},\varPhi_2(2)=\frac{1}{4}\)）
• 樣本個數到達\(b^m\)個點時對\([0,1)\)會形成uniform的劃分。例如\(\varPhi_2(0)=0,\varPhi_2(1)=\frac{1}{2},\varPhi_2(2)=\frac{1}{4},\varPhi_2(3)=\frac{3}{4}\)
• 很多時候並不能代替偽隨機數，因為點的位置和索引有很強的關係，例如在以2為底的Van der Corput序列中，索引為奇數時候序列的值大於等於0.5，偶數時則小於0.5

6.2. 使用Hammersley產生隨機數的虛擬碼

產生二維隨機數，效率較高

// 此處以2為底數產生二維隨機數所用，直接用位操作效率極高，如要拓展到多維，則看下一種做法

float RadicalInverse_VdC(uint bits)

{

     bits = (bits << 16u) | (bits >> 16u);

     bits = ((bits & 0x55555555u) << 1u) | ((bits & 0xAAAAAAAAu) >> 1u);

     bits = ((bits & 0x33333333u) << 2u) | ((bits & 0xCCCCCCCCu) >> 2u);

     bits = ((bits & 0x0F0F0F0Fu) << 4u) | ((bits & 0xF0F0F0F0u) >> 4u);

     bits = ((bits & 0x00FF00FFu) << 8u) | ((bits & 0xFF00FF00u) >> 8u);

     return float(bits) * 2.3283064365386963e-10; // / 0x100000000

}

// ----------------------------------------------------------------------------

vec2 Hammersley(uint i, uint N)

{

    return vec2(float(i)/float(N), RadicalInverse_VdC(i));

}

const uint sample_count = 1024u;

for(uint i = 0u; i < sample_count; ++i)

{

	Vector2f u_v = Hammersley(i, sample_count);

	.....

	// random_x = u_v.x;

	// random_y = u_v.y;

}

擴充套件到多維

// 0.35

double make_hammersley_sequence(int index, int base) {

    double f = 1, r = 0;

    while (index > 0) {

        f = f / base;

        r = r + f * (index % base);

        index = index / base;

    }

    return r;

}

void hammersley_random()

{

    for (uint i = 0u; i < max_num; ++i)

    {

        draw_points.emplace_back(float(i) / float(max_num), make_hammersley_sequence(i, 2));

        // 三維

        draw_points.emplace_back(float(i) / float(max_num),make_hammersley_sequence(i, 2)

         , make_hammersley_sequence(i, 3));

    }

}

6.3. 效果對比圖

可以看到Hammersley序列比c++自身標準庫的random實現更加的分佈均勻

翻譯自：http://holger.dammertz.org/stuff/notes_HammersleyOnHemisphere.html

引用部落格：

• 低差異序列（一）- 常見序列的定義及性質
• 低差異序列（二）- 高效實現以及應用
• 低差異序列的構造

可拓展閱讀：

• 常見的低差異序列 (low-discrepancy sequences ) 之間的區別？