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

接續上一篇

平行設計 Parallel Design

回顧原始算法,為了加快收斂速度,每一次會挑選單位面積能量最多的三角形 $f$ ,之後便拿所有三角形 $t$ 進行傳導,直到傳導能量小於閥值,迭代將停止。算法如下所述:

1
2
3
4
5
6
while (not converge) {
f = PickMaxRadiosityTriangle()
foreach triangle t in model
shader(f, t);

clearRadiosity(f)
}

在嘗試傳遞的過程中,若三角形 $f$ 的三頂點的能量差異大,則選擇自適應切割三角形,直到三頂點能量差異小,這麼計算 Form-Factor 才會正確。在自適應部份,切割方法如下圖所示:

當偵測到綠色三角形 $A$ 頂點之間的 Form-Factor 差異過大時,使用最長邊的中點切割,這麼做是盡可能產生銳角三角形,為了圖形完整,必然也要對鄰居切割。為減少計算量,只算新增中心點 $P$ 的 Form-Factor。對於下方的三角形 $B$ 而言,分成兩種情況,已在這一輪完成計算,則重新計算 Form-Factor;相反地,不做任何事。

Single-Source Parallel Algorithm

我們發現計算 Form-Factor 相當獨立的,但自適應處理需要遞迴切割,因此選用多核心平台,而非通用圖形處理器平台,因為目前的 GPU 實作遞迴所需的 stack 使用 global memory 作為存取位址,所以一般多核心平台效果會更好。在這一次報告中,我們選用 OpenMP 這套跨平台多執行緒 API 進行實驗。

著手將多個三角形平行處理,意即每一個執行緒負責多個三角形的 Form-Factor 計算。

1
2
3
4
5
6
// Single-Source Parallel Algorithm
while (not converage) {
f = PickMaxRadiosityTriangle()
parallel foreach triangle t in model
shader(s, t);

clearRadiosity(f)

Multi-Source Parallel Algorithm

從原始算法中,我們發現到每一次迭代將只有一個熱源輻射到場景中,當場景有多個高能量熱源時,場景必須經過好幾次迭代才能近似最終結果。藉由平行處理的效能,我們可以一次迭代多熱源,便可將低執行緒之間分配工作不均的情況,不僅僅前幾次迭代就能近似最終結果,同時也能加速運算。

1
2
3
4
5
6
7
8
9
10
// Multi-Source Parallel Algorithm
while (not converge) {
set<Triangle> f = RadiosityTriangleCandiateCandidate();
parallel foreach triangle t in model
if (f.find(t))
continue;
foreach s in f
shader(s, t);
clearRadiosity(f);
}

我們所用的平行方無法搭配上述自適應的切割方案,其原因在於分裂過程中,同時也要對鄰居三角形分裂,整個圖形產生的節點與邊的關係才會正確,無法保證鄰居在同一執行緒內處理,若沒有做好空間切割,我們便無法處理這部份。若模型格式會是數個獨立的物體,而非單純的三角形資訊,可分配每一個執行緒處理多個獨立物體,我們預期可以達到更好的效果。

平行過程中每一個執行緒共享和衝突的區段越少越好,這意味著我們必須在運行輻射前就必須將模型切得相當細緻。特別注意到,切得細緻與否對於光線投射 (Ray Casting) 複雜度不變,因為邏輯上他們處理同一平面。

一旦切得細緻,傳遞的效果就不是這麼好,在邊界的陰影更加顯著。根據理論和實作層面推測,其一原因是能傳遞的總能量隨著迭代減少,那麼從分裂過程中傳遞能量採用較多的加法完成,相較於多個 32-bit floating point 誤差就少了許多。

我們也試著使用獨立的切割方案─重心切割,切割的結果不依賴鄰居,只需要在加入三角形清單部份使用 critical section 即可,效能影響並不大。

根據重心切割,下述實驗中,從 156 個三角形,自動分裂到 30000 個三角形後進行輻射的結果如下圖所示:

明顯地,根據重心的切割方法容易產生鈍角三角形,看起來就會像很多紡錘體。在眾多數學性質中,只使用重心也許不是好的解決方案,這是値得探討的一部份,由於製作上的時間限制,我們並沒有去探討各個不同切割方案,所對應的自適應的效果如何。

Longest Edge Center

展示結果 Demo

只使用優化技術渲染結果

blocks room
hall church

平行效能比較

在 Intel Xeon E5-2620 v3 上,我們測試不同平行度帶來的影響,由於只有兩個實體 CPU,每一個 CPU 有 6 個核心,每個核心皆有 Hyper-threading 技術,故可產生 24 個執行緒。

Single-Source Parallel Algorithm - Scalability

我們對模型 room.tri 以預先切割 14977 個三角形後,根據先前提到的平行算法 Single-Source Parallel Algorithm,即是迭代一次只取一個熱源,平行計算所有三角形到此熱源的 Form-Factor 値,針對不同的執行緒個數和運行時間記錄,結果如上圖。在由於過多的執行緒可能會帶來更多的 false sharing,造成資料在不同的 CPU 之間運行 data transmission,所以效果就逐漸不明顯。

Multi-Source Parallel Algorithm - Scalability

接著,我們測試 Multi-Source Parallel Algorithm,以 room.tri 預先切割 50017 個三角形後,每次迭代皆取數個熱源,平行計算所有三角形到所有被選取熱源的 Form-Factor 的總和。同樣地,因為查找的資料重複存取的模式造成不好的影響,類似上述所提到的 false sharing 影響,故會呈現一種陡坡。

參考論文 Reference

後記 Note

當我們進行平行效能比較時,要特別小心編譯器行為的差異,意即平行處理 $P$ 控制時,當 $P=1$ 時,使用平行版本比較合適。因為有時候,平行部分的函數被編譯器偵測到,由於分享記憶體的關係,部分代碼無法優化,導致效能瞬間慢個四到五倍都是有可能的,再加上 false sharing 的關係,更有可能在發生密集計算時,效能更加低落。

在不同平台上的情況也有所不同,例如在一般 server 上運行時,CPU 頻率通常都會比一般 PC 慢上許多,又因為很多個 CPU 導致總共的快取大小遠比一般 PC 多,所以平行效能將會受限於運行的應用行為,是否需要時常存取大量的資料。

致謝 Thanks

一開始在挑選主題相當困惑,畢竟使用 Unity 可以做出更生動的作品,但作品的元素和創意相當重要,相當迷惘於要選哪一種類型才好。如果要兩人以上一起做,那麼主題又不能太過狹隘。百般思慮下,還是由我拉選了這個主題下來做。

「我可能會扯你後腿,還是我們各別做?」組員擔心地說道

不用擔心,我自己也沒信心將所有程序都看完且修改更好更快,每天都焦頭爛額地煩惱整份程序的運作,深怕來不及在時限內完成足夠的報告份量。

「快來幫幫我啊,在身邊一起 trace 程序也好,咱的記憶體不足啊。」內心如此吶喊

一起修課的博班學長給了我們一些意見與鼓勵,而學弟們只會在一旁扯後腿問「學姊今天會來嗎?」最後,我們完成了整份程序的理解與討論。

Read More +

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

此為計算機圖學課程中的一環,自由挑選主題。而這個 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
Read More +