CUDA Volume Rendering [Part.4 Render(GPU part)]

前面已經把 CPU 要做的事都講完了，接下來，就是最重要的 GPU 的 Kernel 程式了！這份 kernel 的定義如下：

__global__ void 
d_render(uint *d_output, uint imageW, uint imageH,
          float density, float brightness,
          float transferOffset, float transferScale) 

他要接受的參數有：

• d_output
  儲存輸出結果的 1D array，對應到 pixel buffer object
• imageW, imageH
  輸出影像的大小，也就是 d_output 的大小
• density
  調整每一點顏色的透明度的比例
• brightness
  最後結果的亮度調整比例
• transferOffset, transferScale
  把 volume 資料對應到 transfer function 的調整參數

而 Heresy 把這個 kernel 程式的動作，分成四大部分：

1. 設定前置參數、計算 index 和視線在空間中的位置
2. 計算視線和 volume 的外框相交的點和攝影機的距離
3. 由後往前，計算視線上的顏色
4. 將結果寫入到 pixel buffer object

第一部分就是設定一些需要的變數了∼這裡包含了 maxSteps、tstep、boxMin、boxMax 四個變數，Heresy 會在用到的時候再去解釋。而接下來，就是根據 thread index 和 block index，計算出這個 thread 的索引值了∼這邊還是很單純的，用 index = (block index) * (block size) (thread index) 來做計算。

uint x = __umul24(blockIdx.x, blockDim.x)   threadIdx.x; 
uint y = __umul24(blockIdx.y, blockDim.y)   threadIdx.y;

不過，這邊他改用 24bit 的整數乘法函式 __umul24() 來取代 32bit 的 operator*；在目前的硬體上，這樣是會比較快一點的。

算出來的 x 和 y，就是代表現在這個點的結果，是要儲存在輸出畫面的哪一點上；而要開始計算這一點的顏色，還要再把他轉換成空間座標的一條線。而這邊就要先講一下這個範例程式的空間定義了∼首先，他把 volume 視為一個 (-1, -1, -1) 到 (1, 1, 1) 的一個正方體，並把這兩個座標儲存在 boxMin、boxMax。而攝影機的位置在預設會有一個 z 方向的位移 4，所以位置會是 (0, 0, 4)；而視平面則是建立在相對攝影機的 (0,0,-2) 的位置，大小是 2*2。而在旋轉時，實際上都是將攝影機的位置做移動，不會動到代表 volume 的 box。

整個空間關係，可以畫成下方的示意圖。

實際上在程式的部分，是寫成下面的樣子：

float u = (x / (float) imageW)*2.0f-1.0f; 
float v = (y / (float) imageH)*2.0f-1.0f;
  
// calculate eye ray in world space 
Ray eyeRay; 
eyeRay.o = make_float3(mul(c_invViewMatrix, make_float4(0.0f, 0.0f, 0.0f, 1.0f))); 
eyeRay.d = normalize(make_float3(u, v, -2.0f)); 
eyeRay.d = mul(c_invViewMatrix, eyeRay.d);

第一段在計算 u, v 的值，就是在把他們由 [0, imageW] 和 [0, imageH] 轉換到 [-1, 1] 的座標定義；而 eyeRay 就是用來儲存攝影機的位置（eyeRay.o）和這條視線所看的方向（eyeRay.d）。eyeRay.o 就是直接帶入矩陣 c_invViewMatrix；而 eyeRay.d 則是根據計算出的 u, v，再考慮視平面 z 軸的位移 -2，建立出長度為 1 的向量，最後再考慮矩陣 c_invViewMatrix。如此，就可以建立出符合上述空間定義的資料了。

而接下來第二部分， 就是要計算 eyeRay 這條視線，和 volume 的 box 的交點了。在這邊，他是呼叫這個 intersectBox() 函式；傳入 eyeRay 和 boxMin、boxMax 後，他會傳回這條視線是否有和 volume 相交，同時把和 volume 相交點到攝影機位置距離的最小值和最大值記錄在 tnear、tfar 中。而 tnear、tfar 也就代表了這個 volume 在這條視線上，對攝影機來說對短距離與最遠距離；這是之後用來計算 ray casting 要走過哪些區域的依據。這部分的程式如下：

// find intersection with box 
float tnear, tfar;
     int hit = intersectBox(eyeRay, boxMin, boxMax, &tnear, &tfar); 
if (!hit) return;
     if (tnear < 0.0f) tnear = 0.0f;     // clamp to near plane

而這邊 intersectBox() 用的計算方法，算是 Heresy 覺得相當有趣的一個方法，他是參考 Siggraph 的一篇 education 來做的；基本上，就是透過將 box 限制在和座標軸同方向，來簡化計算。在這個方法裡，他可以簡單的向量運算，就快速的求出一條射線和六個平面焦點對於某個點的距離。不過，在這邊細節就先不提了。

再來的第三部分，就是實際要計算 volume 在 (x, y) 這點的顏色的部分了！這段的程式碼如下：

// march along ray from back to front, accumulating color 
float4 sum = make_float4(0.0f);; 
float t = tfar; 
for(int i=0; i<maxSteps; i  ) {
     float3 pos = eyeRay.o   eyeRay.d*t;
     pos = pos*0.5f 0.5f;    // map position to [0, 1] coordinates
 
     // read from 3D texture
     float sample = tex3D(tex, pos.x, pos.y, pos.z);
 
     // lookup in transfer function texture
     float4 col = tex1D(transferTex, (sample-transferOffset)*transferScale);
 
     // accumulate result
     sum = lerp(sum, col, col.w*density);
 
     t -= tstep;
     if (t < tnear) break; 
} 
sum *= brightness;

基本上，他的計算方法就是由一條視線和 volume 相交最遠的點（距離為 tfar），以 tstep 為間距向攝影機前進，直到和攝影機的距離小於 tnear，或累積次數大於 maxSteps。每前進一步，都可以透過「pos = eyeRay.o eyeRay.d*t」計算出在空間座標中的位置；不過由於這是計算一個在 volume box 中的點，所以他的座標會是介於 (-1,-1,-1) 和 (1,1,1) 之間；而由於 texture 的 normailized 座標是介於 [0, 1] 之間，所以必須再把他做個轉換。

而接下來，就是透過 tex3D() 到 volume data 的 3D texture（tex）去取得 (pos.x, pos.y, pos.z) 這一點的值 sample。而這裡取出來的值會是灰階的，要再用 tex1D() 到代表 transfer function 的 1D texture transferTex 中，去找到對應的色彩；不過，這邊會先把取得的 sample 透過 transferOffset 和 transferScale 做一個調整。

取得這一點的色彩後，接下來就是要和目前的結果 sum 做累加的動作；在這邊，他是用一個定義在 cutil_math.h 裡的函式 lerp() 來進行這項計算。而「sum = lerp(sum, col, col.w*density);」實際的動作，就是：「sum = (col.w * density) * (col – sum);」；其中，density 這個參數算是在調整每一點的 alpha 值了∼

而最後，就是計算計數器 t 的值，並確認他還沒有超過 tnear；如果已經超過的話，就提前結束整個迴圈。而最後，他還會再把最後的顏色，根據 brightness 做一的調整，然後才是最後的 (x, y) 這個點的顏色 sum。

最後的第四部分，就是將計算出來的顏色值 sum，寫入到代表 pixel buffer object 的 d_output 裡了。而他的第一個動作，就是確認 (x, y) 的質是否有超過 (imageW, imageH)（其實這個動作應該可以更早做？）；之後，就是計算出 (x, y) 在一維陣列中 d_output 中的索引值 i，並把值寫進去。
不過，由於 d_output 的型別是 uint*，而 sum 是 float4，所以這篇他是先透過 rgbaFloatToInt() 來進行轉換，把 float4 的色彩值用一個 unsigned int 來儲存，再存到 d_output。

這段程式的內容如下：

if ((x < imageW) && (y < imageH)) {
     // write output color
     uint i = __umul24(y, imageW)   x;
     d_output[i] = rgbaFloatToInt(sum); 
}

而到此，也就是 CUDA Volume Rendering 的 kernel 的全部內容了∼而 Heresy 關於 CUDA 提供的這個 Volume Render 範例，大概也就先講到這了∼ 🙂

原始發表：CUDA Volume Rendering [Part.4 Render(GPU part)]