學校課程/計算機圖學

探討平行與優化技術於熱輻射法 (上)

contents

  1. 1. Parallel Computing and Optimization Skills for Radiosity
  2. 2. 介紹 Introduction
  3. 3. 優化技術 Code Review & Optimization
    1. 3.1. 減少光線投射計算量 Strength Reduction for Ray Casting
    2. 3.2. 減少 FF 計算量 Strength Reduction for Form-Factor
    3. 3.3. 改善資料局部性 Improve Data Locality
    4. 3.4. 其他優化 Other Optimization
      1. 3.4.1. Improve I-cache Miss
      2. 3.4.2. Short Circuit Design
      3. 3.4.3. Clipping Algorithm
      4. 3.4.4. Strength Reduction for Float-Point
      5. 3.4.5. Shrink the Scope of Variables
      6. 3.4.6. Reduce the Number of Arguments
      7. 3.4.7. Remove Implications of Pointer Aliasing
  4. 4. 小結
    1. 4.1. Experiment

此為計算機圖學課程中的一環，自由挑選主題。而這個 Radiosity 有好幾年前的原始碼，程式碼部分由老師提供。看到幾處寫得不是很好的地方，於是就拿來加速。在 tracing code 耗費了相當大的力氣，雖然才幾千行的 C 程序，完全不熟的領域，慢慢做起來也別有一番風趣。morris821028/hw-radiosity

Parallel Computing and Optimization Skills for Radiosity

  • R04922067 楊翔雲、R04922133 古君葳

介紹 Introduction

熱輻射法 (Radiosity) 是一種渲染的技術。相較於光線追蹤法 (Ray Tracing)，熱輻射法可以產生更接近於現實場景中光亮的變化。當場景使用光線追蹤法時，物體的陰影的邊緣相對銳利，但在現實情況下，物體陰影漸層呈現，因此使用熱輻射法可以更貼近我們想要的。

$$\begin{align*} B_i dA_i = E_i dA_i + R_i \int_j B_j F_{ji} dA_j \end{align*}$$
  • $A_i$ : Area of element i (computable)
  • $B_i$ : Radiosity of element i (unknown)
  • $E_i$ : Radient emitted flux density of element i (given)
  • $R_i$ : Refletance of element i (given)
  • $F_{ji}$ : Form Factor from j to i (computable)

假設整個場景中有 $N$ 個三角形，每一次迭代選擇一個最亮的三角形當作光源，由這個光源計算與場景中其他三角形的 Radiosity 之值。其中，判斷光源是否可以輻射到某個三角形之複雜度介於 $O(\log N)$ 到 $O(N)$ (視Data structure而定)，而計算 Form-Factor的花費可以視為常數 $O(1)$ ，因此每次迭代的複雜度介於 $O(N \log N)$ 到 $O(N^2)$。

其中佔據效能的因素是 Form-Factor 估算，因此有像 Hemicube之類的近似逼近，大幅度減少計算量，但投影回到原本物件上會失真。

$$\begin{align*} F_{ij} = \frac{1}{A_i} \int_{A_i}\int_{A_j} \frac{\cos \phi_2 \cos \phi_1}{\pi r^2} dA_i dA_j \end{align*}$$

進入優化主題吧

優化技術 Code Review & Optimization

首先，我們先對助教提供的程式碼加速，分成以下幾個部分討論

  • 減少光線投射計算量 Strength Reduction for Ray Casting
  • 減少 Form-Factor計算量 Strength Reduction for Form-Factor
  • 改善資料局部性 Improve Data Locality
  • 其他優化 Other Optimization:
    • Improve I-cache Miss
    • Short Circuit Design
    • Clipping Algorithm
    • Strength Reduction for Float-Point
    • Shrink the Scope of Variables
    • Reduce the Number of Arguments
    • Remove Implications of Pointer Aliasing
    • Copy Optimization

減少光線投射計算量 Strength Reduction for Ray Casting

判斷射線 (Ray) 是否能打到三角形 $A$ 上，先用 bounding box 包住 $A$ ，計算 $p$ 到 bounding box 的時間 $t$ ，若 $t$ 大於目前的最小 $t_{\min}$ ，則退出。相反地，再計算更精準的 $t$ 。加入利用已知結果 $v = p + t_{\min} \cdot d, t_{\min} > 0$

1
2
3
4
5
6
7
8
9
10
int TriangleHitted(Vector p, Vector d, TrianglePtr tp, float *t_min) {
    float t = /* time t from p to bounding box of Triangle tp */;
    if (t < eps)
        return false;
    if (t >= *t_min)    /* important !! */
        return false;
    /* ... */
    *t_min = t;
    return true;
}

減少 FF 計算量 Strength Reduction for Form-Factor

根據公式 $F_{ij} = \frac{1}{A_i} \int_{A_i}\int_{A_j} \frac{\cos \phi_2 \cos \phi_1}{\pi r^2} dA_i dA_j$ ，一般我們的判斷順序會得如下：

1
2
3
4
5
6
7
8
9
float computeFormFactor(TrianglePtr srcTri, int logSrc, TrianglePtr desTri, int logDes, Vector p) {
    Vector dir = srcTri->c - p;
    float ff = 0;
    if (RayHitted(p, dir, logDes) == logDes) {
        float theta1 = CosTheta(dir, srcTri->normal), theta2 = CosTheta(dir, desTri->normal);
        ff = theta1 * theta2 * srcTri->area / (norm2(dir) * PI);
    }
    return max(ff, 0.f);
}

效能考量的因素：

  • RayHitted() 需要大量的計算
  • Form-Factor 在 float 儲存格式下可能無法得到貢獻，改採優先計算 Form-Factor 的值，再運行 RayHitted 判斷。調整加速了 2 倍多。
1
2
3
4
5
6
7
8
9
10
11
float computeFormFactor(TrianglePtr srcTri, int logSrc, TrianglePtr desTri, int logDes, Vector p)
    Vector dir = srcTri->c - p;
    float theta1 = CosTheta(dir, srcTri->normal), 
theta2 = CosTheta(dir, desTri->normal);
    float ff = theta1 * theta2;
    if (ff <= 0)    return 0.f;
    ff *= srcTri->area / (norm2(dir) * PI);
    if (ff <= 0)    return 0.f;
    if (RayHitted(p, dir, logDes) == logDes)
        return ff;
    return 0.f;

改善資料局部性 Improve Data Locality

程式碼中使用 3D-DDA Line Algorithm/Bresenham’s Line Algorithm 搭配 Octree，在找尋某個射線與哪個最近三角形相交。

  • 只需要儲存葉節點代表的立方體中，所有可能相交的三角形編號
  • 移除掉中間產生的編號，讓每一次 access 的 cache-miss 下降

在 3D-DDA 中，我們需要反查找空間中一點 p 在哪一個葉節點中，藉由固定的長寬高切割長度，可以在 $O(1)$ 時間內得知 [i][j][k] 各別的值。若限制大小，則建立陣列 [i][j][k] 查找。若自適應大小，則建立 hash 表查找，但根據實驗結果，效能並沒有改善，因為三角形個數過多導致命中機率過低。

對於 Static Tree 的 Memory Layout，大致上分成四種DFS Layout、Inorder Layout、BFS Layout、和 van Emde Boas Layout，目前程式使用的是 DFS Layout，這方面會影響到存取的效能。若有更多的時間，我們也可以測試平行處理的細粒度與這些 Memory Layout 的影響。

其他優化 Other Optimization

Improve I-cache Miss

壓縮程式碼長度以改善 I-cache miss，因為大部分的初始化只會運行一次，不應該交錯在時常被呼叫的函數之間，指令載入效能才會提高，同時也要做好 Code Layout，就能改善執行效能。

原始版本如下：

1
2
3
4
5
6
7
8
9
x, y = compute(0)
buildTree(x, y)
x, y = compute(1)
buildTree()
x, y = compute(2)
buildTree()
x, y = compute(3)
buildTree()
...

壓縮和整理後

1
2
3
for i from 0 to n
    x, y = compute(i)
    buildTree()

Short Circuit Design

判斷三角形與一個正交立方體是否相交，使用投影到二維空間中與一線段的交點是否皆存在。投影方法有 3 種 x-y, y-z, z-x，線段投影共有 2 種，共計 6 種情況。原先程式沒有做好短路設計，只要其中一種不符就應退出。

1
2
3
4
5
6
7
8
9
10
11
12
int CrossOver(TrianglePtr tri, 
            Vector g0, Vector g1) {
    for (xyz = 0; xyz < 3; xyz++) {
        // front face project
        if (!test())
            return false;
        // back face project
        if (!test())
            return false;
    }
    return true;
}

Clipping Algorithm

我們實作課程中所描述的 Cohen–Sutherland Algorithm 降低 branch 次數，使用 bitwise 操作引出 SSE (Streaming SIMD Extensions)。儘管 compiler -O2 替我們優化，為減少 stack push/pop 的次數，實作時請不要使用的 procedure call，否則會慢上許多。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
void computeRadiosity(TrianglePtr srcTri, TrianglePtr desTri, 
float ff[3]) {
    char mask1 = (p0[x] < g0[x])<<0 |
        (p0[x] > g1[x])<<1 |
        (p0[y] < g0[y])<<2 |
        (p0[y] > g1[y])<<3 ;
    char mask2 = (p1[x] < g0[x])<<0 |
        (p1[x] > g1[x])<<1 |
        (p1[y] < g0[y])<<2 |
        (p1[y] > g1[y])<<3 ;
    if (mask1&mask2)    
        return false;
    if (!(mask1|mask2))
    return true;
    // ... test

Strength Reduction for Float-Point

兩個外積結果相乘小於零，減少 instruction cycle 量，盡量用整數作為運算型態。

在現在的 Intel CPU 中，32-bit 浮點數運算基本上跟整數一樣快

1
2
3
4
5
6
7
8
float a = cross(/* */);
float b = cross(/* */);
if (a * b < 0)
    return false;
b = cross(/* */);
if (a * b < 0)
    return false;
...

事先判斷正負號，同時也防止溢位。

1
2
3
4
5
6
7
8
int a = cross(/* */) < 0;
int b = cross(/* */) < 0;
if (a != b)
    return false;
b = cross(/* */) < 0;
if (a != b)
    return false;
...

Shrink the Scope of Variables

減少變數生命週期的長度以增加放入暫存器的機會，而非 stack 上。

1
2
3
4
5
6
7
8
9
10
float rgb[3];
for (int i = 0; i < 3; i++)
    rgb[f(i)] = g(i);
/* ... */
if (maybe) {
    for (int i = 0; i < 3; i++) {
        rgb[h(i)] = g(i);
    }
    /* ... */
}

當邏輯很複雜時，編譯器不太能幫忙做分析，所以自己手動優化的效果會比較好，在 C/C++ 語言中，可以利用大括弧進行區域變數的設定。

1
2
3
4
5
6
7
8
9
10
11
12
13
{
    float rgb[3];
    for (int i = 0; i < 3; i++)
        rgb[f(i)] = g(i);
    /* ... */
}
if (maybe) {
    float rgb[3];
    for (int i = 0; i < 3; i++) {
        rgb[h(i)] = g(i);
    }
    /* ... */
}

Reduce the Number of Arguments

減少 stack push/pop 次數

1
2
3
4
5
6
struct Arg {
    int a0, a1;
};    // p1.a1 = p2.a0
int f(Arg p1, Arg p2) {
    /* ... */
}

如何修改合適的參數個數，必須看使用的機率和次數，才能達到最大效益。

1
2
3
4
5
6
struct Arg {
    int a0, a1, a2;
};
int f(Arg p1p2) {
    /* ... */
}

Remove Implications of Pointer Aliasing

移除指標 Aliasing，意指可能會指向相同記憶體位址，導致每次計算都要重新載入，不能放進暫存器中。如下述的寫法，編譯器無法判定 srcTri 是否與 desTri 相同，在累加時則重新載入 srcTri->deltaB[] 的數值，計算上可能會產生數次的 cache miss，隨著迴圈次數不斷突顯效能差異。

1
2
3
4
5
6
7
8
9
10
11
void computeRadiosity(TrianglePtr srcTri, TrianglePtr desTri,
 float ff[3]) {
    for (int v = 0; v < 3; v++) {    // vertex
        for (int c = 0; c < 3; c++) {    // color RGB
            flaot deltaB = desTri->Frgb[c]/255.0*RefRatio*srcTri->deltaB[c]*ff[v]/3;
            desTri->deltaB[c] += deltaB; 
            desTri->accB[v][c] += deltaB; 
            desTri->deltaAccB[v][c] += deltaB; 
        }
    }
}
  • 方法 1: 加入 if (srcTri != desTri) 判斷，讓編譯器在 Function Pass 階段著手的 Dependency Analysis 更好
  • 方法 2: 使用 Copy Optimization，同時把重複計算搬到 stack 上，或者使用 Polyhedal 表示法進行 Reordering Accesses for Loop Nest。這裡我們選擇前者，更容易引出 SSE
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
void computeRadiosity(TrianglePtr srcTri, TrianglePtr desTri, 
float ff[3]) {
    const float k = RefRatio / 255.0;
    float lo[3] = { desTri->Frgb[0]*k*(srctri->deltaB[0]), 
                    desTri->Frgb[1]*k*(srctri->deltaB[1]), 
                    desTri->Frgb[2]*k*(srctri->deltaB[2])};
    for (int v = 0; v < 3; v++)    {    // vertex
        for (int c = 0; c < 3; c++) {    // color RGB
            /* calculate the reflectiveness */
            float deltaB = lo[c] * ff[v] / 3;
            desTri->deltaB[c] += deltaB;
            desTri->accB[v][c] += deltaB;
            desTri->deltaaccB[v][c] = deltaB;
        }
    }
}

小結

  • 使用洽當的編譯器參數可加速 2 倍
  • 減少 Form-Factor 計算加速 2 倍
  • 剩餘優化部份改善 10% ~ 20% 效能。

至今，我們加速了 4 倍，在下一篇文章中，我們將繼續探討平行處理。

Experiment

Model Origin (sec.) Our v0.1 (sec.) Speedup
room.tri 10.27 4.70 2.18
hall.tri 176.92 38.50 4.59
church.tri 72.32 42.64 1.69