CUDA Volume Renderering

Volume Render 是 CUDA 2.0 Beta 版新加入的一個範例，主要是直接透過 3D Texture 來做 Volume 的 Ray marching。

關於 Volume Rendering 這項技術，可以先去參考維基百科的說明：「立體渲染」。不過還是大概提一下，Volume Rendering 實際上在醫學方面其實相當普遍；主要是一般的斷層掃描都是一張一張的平面影像，把這些平面影像按照順序堆疊起來，就會變成一個 3D 的影像資料（也就是所謂的 Volume data）。Volume Rendering 就是一種用來繪製這種 3D 影像資料的技術。（應用參考）

而由於傳統的電腦圖學都是以多邊形的方式來建構、繪製 3D 場景；這點和呈現 Volume 所需要的技術是不同的！所以 Volume Rendering 沒辦法直接使用傳統的電腦圖學方法來呈現。

而一般來說，常見的 Volume Rendering 有兩種概念：第一種是將 Volume data 建立出多邊形的資料，然後再用這些資料還繪圖；第二種則是直接拿整個 Volume 資料去畫，這種方法一般叫做 Direct Volume Render。

在 Direct Volume Rendering 又有數種不同的技術，下面列舉常見的兩種：

Ray castingSlice base

由觀測點（eye）往 volume 看，每一條「視線」，都根據某個固定的間距在 Volume 中取樣，然後依此累算出這條視線最後會呈現的顏色。

實際上，每一條線會對應到最後呈現的圖上的一個點；所以最後要呈現的解析度要多高，就得做幾次這樣的運算。

根據和視線的垂直方向，把 Vloume data 重新取樣，產生出多張和視線垂直的 slice；接著再由後而前，依序把這些 Slice 用傳統的多邊形方法來繪製。

上圖資料來源：Volume Rendering For Games

而在 CUDA 2.0 的 SDK 中所提供的 Volume Rendering 範例（專案名稱就是 volumeRender，檔案在 C:\Program Files\NVIDIA Corporation\NVIDIA CUDA SDK\projects\volumeRender），就是 ray casting 的方法；不過，他藉由 CUDA 這套 SDK，可以發揮 nVidia GPU 的大量平行化的好處，來做 GPU 的 ray casting。

程式的基本概念，是把 Volume 的資料讀到電腦的主記憶體後，當成 3D Texture 來 bind 到 device memory 中，然後再透過 CUDA kernel 來透過存取 3D texture，進行 ray casting 的計算。而運算的結果，會直接當成 OpenGL 的 Buffer Object，直接畫出來；如此也避免了必須要先把結果由 device memory 複製回 host memory，再送到 OpenGL render 的傳輸時間。

他的程式檔有兩個：

• volumeRender.cu
  主程式、資料讀取和 OpenGL 等相關的部分。
• volumeRender_kernel.cu
  實際做 ray casting 的 CUDA kernel。

在 volumeRender.cu 中的 function 列表如下：

mainmain function。loadRawFile用來讀取 Volume data 的函式。
在這個範例中，Volume 的資料是用 3D 的 RAW 檔來儲存的。initPixelBuffer建立 OpenGL 的 pixel buffer object，用來對應到 CUDA ray casting 儲存結果的記憶體空間。也直接拿來畫出結果。cleanup清理資料用的。iDivUp做除法的無條件進位，用來算 grid 大小用的。initCuda起始化 CUDA 的資料。
主要是把讀進來的資料建立成 3D Texture，以及建立顏色對應用的 transfer function（包含陣列資料以及 1D texture）。render把旋轉矩陣由 host memory 複製成 device 上的 constant 變數。
然後呼叫 kernel 函式，進行計算。displayglut 的 callback function，用來顯示用的，每次要更新畫面，都是執行這個函式。
他會計算物體的旋轉矩陣，然後再呼叫 render() 來做 ray casting 計算，最後再把儲存了結果的 pixel buffer object 畫出來。reshape這四個 function 是 glut 的 callback function。
reshape 是當視窗大小位置改變時會被執行到、motion 是滑鼠移動時會被呼叫的、mouse 是滑鼠按鈕會執行的函式、keyboard 則是鍵盤的。
這邊不多加解釋了。motionmousekeyboard

而在 volumeRender_kernel.cu 中，則是：

intersectBox計算一條視線和 Volume 的 box 的交點，並傳回相交的最近點和最遠點。mul計算矩陣乘上一個向量。rgbaFloatToInt把 rgba 四項的 color，轉換成一個 int 來儲存。d_renderCUDA 的 device kernel function。
為每一個像素，用 Ray casting 的方法來計算他的顏色。

 

 

在前面的部分，已經先做了一些簡單的介紹；接下來，開始看程式碼吧～這一部分，主要是看 volumeRender.cu 中 main function，並大概講一下 volume data 的前置處理。

在 main() 函式裡，主要分成幾個部分：

1. 透過 cutil 來處理執行時的參數。
2. 讀取 volume data，並初始化 volume 資料
3. 設定 glut 環境
4. 起始設定 pixel buffer object
5. 執行 glut 的 main loop

第一部分主要就是透過 cutil 的 cutGetCmdLineArgumenti() 這個函式，來處理main() 的 argc 以及 argv。而這個程式是設計成可以指定檔案的檔名以及大小，例如：「-file=brain.raw -xsize=512 -ysize=512 -zsize=125」就是指定檔名為 brain.raw，x, y, z 三軸的大小依序為 512, 512, 125。

而第二部分則是先透過 cutFindFilePath() 去找出檔案完整的路徑，然後再透過loadRawFile() 來讀取 RAW 資料成為一個 unsigned char 的陣列；這邊的 RAW 資料，原則上就是一張一張 2D 灰階圖組合成的單一檔案。接著，就是透過initCuda()，來把讀取進來 uchar 陣列（本程式一開始就把 uchar 定義為 unsigned char，以下將以 uchar 沿用），轉換成 ray casting 需要的 3D Texture 了～
（由於他是使用 cudaArray 來 bind 到 3D texture，所以建議可以先參考看看之前寫的《CUDA Texture Part.1 簡介》）

接下來，就是直接看 initCuda() 這個函式了～這個函式裡做的事，主要包含了兩部分：

1. 建立 volume 資料本身的 3D texture
2. 建立將 volume 資料的灰階，對應到彩色的 transfer function，以及他的 1D texture

將 volume data 建立成 3D texture 的程式如下：

// create 3D array cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc<uchar>(); CUDA_SAFE_CALL( cudaMalloc3DArray(&d_volumeArray, &channelDesc, volumeSize) ); // The sample does not work with pageable memory // This is a known issue for beta that will be fixed for the public release #define USE_PAGE_LOCKED_MEMORY 1 #if USE_PAGE_LOCKED_MEMORY // copy to page-locked mem cudaPitchedPtr pagelockedPtr; pagelockedPtr.pitch = volumeSize.width*sizeof(uchar); pagelockedPtr.xsize = volumeSize.width; pagelockedPtr.ysize = volumeSize.height; size_t size = volumeSize.width*volumeSize.height*volumeSize.depth*sizeof(uchar); CUDA_SAFE_CALL( cudaMallocHost(&(pagelockedPtr.ptr), size) ); memcpy(pagelockedPtr.ptr, h_volume, size); #endif // copy data to 3D array cudaMemcpy3DParms copyParams = {0}; #if USE_PAGE_LOCKED_MEMORY copyParams.srcPtr = pagelockedPtr; #else copyParams.srcPtr = make_cudaPitchedPtr((void*)h_volume, volumeSize.width*sizeof(uchar), volumeSize.width, volumeSize.height); #endif copyParams.dstArray = d_volumeArray; copyParams.extent = volumeSize; copyParams.kind = cudaMemcpyHostToDevice; CUDA_SAFE_CALL( cudaMemcpy3D(&copyParams) ); // set texture parameters tex.normalized = true; // access with normalized texture coordinates tex.filterMode = cudaFilterModeLinear; // linear interpolation tex.addressMode[0] = cudaAddressModeClamp; // wrap texture coordinates tex.addressMode[1] = cudaAddressModeClamp; // bind array to 3D texture CUDA_SAFE_CALL(cudaBindTextureToArray(tex, d_volumeArray, channelDesc));

在這邊，讀進來的資料是 h_volume，他是一個 uchar 的陣列。而和 2D CUDA Array 時一樣，先透過一個 cudaChannelFormatDesc 來描述cudaArray 的資料型別，然後再透過 cudaMalloc3DArray() 來建立一個大小為 volumeSize 的 3D cudaArray：d_volumeArray。

不過，在註解裡有寫，目前的 2.0 beta 版，還不能直接把一般的記憶體空間的資料，直接複製到 3D 的 cudaArray 裡；必須要先透過cudaMallocHost() 宣告一個 page lock 的記憶體空間，並把h_volume 的資料複製到這個 page lock 的變數 pagelockedPtr。然後再設定 cudaMemcpy3DParms，把資料用 cudaMemcpy3D() 複製到d_volumeArray。比起來，是多了一個先轉成 page lock 變數的過程，不過這個問題在正式版發佈時，應該是會解決的掉的。

而把資料複製到 cudaArray 後，就是設定一下 3D texture 的參數，然後再透過cudaBindTextureToArray()，把 cudaArrayd_volumeArray bind 到 texture<uchar, 3, cudaReadModeNormalizedFloat> tex 了。

而在把 Volume 資料本身處裡完了之後，接下來還有一份額外的陣列。由於在醫學的 volume rendering 中，CT、MRI 這些資料都是灰階的；如果要用彩色來呈現、凸顯某些部位的話，大部分都是用一個transfer function 來做色彩的對應。

在 initCuda() 這個函式中的後段，就是在處理這份資料。

// create transfer function texture float4 transferFunc[] = { { 0.0, 0.0, 0.0, 0.0, }, { 1.0, 0.0, 0.0, 1.0, }, { 1.0, 0.5, 0.0, 1.0, }, { 1.0, 1.0, 0.0, 1.0, }, { 0.0, 1.0, 0.0, 1.0, }, { 0.0, 1.0, 1.0, 1.0, }, { 0.0, 0.0, 1.0, 1.0, }, { 1.0, 0.0, 1.0, 1.0, }, { 0.0, 0.0, 0.0, 0.0, }, }; cudaChannelFormatDesc channelDesc2 = cudaCreateChannelDesc<float4>(); cudaArray* d_transferFuncArray; CUDA_SAFE_CALL(cudaMallocArray( &d_transferFuncArray, &channelDesc2, sizeof(transferFunc)/sizeof(float4), 1)); CUDA_SAFE_CALL(cudaMemcpyToArray( d_transferFuncArray, 0, 0, transferFunc, sizeof(transferFunc), cudaMemcpyHostToDevice)); transferTex.filterMode = cudaFilterModeLinear; transferTex.normalized = true; // access with normalized texture coordinates transferTex.addressMode[0] = cudaAddressModeClamp; // wrap texture coordinates // Bind the array to the texture CUDA_SAFE_CALL( cudaBindTextureToArray( transferTex, d_transferFuncArray, channelDesc2));

在這個範例中，他是先宣告一個 float4 的陣列 transferFunc，裡面算是有九個控制點，分別代表灰階值不同所要呈現的階段顏色（介於中間的值，會透過texture 用線性內插來算）。而這邊他就是用一般的 1D cudaArray 來做，也就不多加介紹了；最後會拿來用的，就是已經 bind 好資料的texture<float4, 1, cudaReadModeElementType> transferTex。

再來的第三部分，則是透過 glut 來建立 OpenGL 的環境。

// initialize GLUT callback functions glutInit(&argc, argv); glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE); glutInitWindowSize(width, height); glutCreateWindow("CUDA volume rendering"); glutDisplayFunc(display); glutKeyboardFunc(keyboard); glutMouseFunc(mouse); glutMotionFunc(motion); glutReshapeFunc(reshape); glutIdleFunc(idle); glewInit(); if (!glewIsSupported("GL_VERSION_2_0 GL_ARB_pixel_buffer_object")) { fprintf(stderr, "Required OpenGL extensions missing."); exit(-1); }

前面兩行，是在起始化整個 OpenGL 的環境，而第三、四行則是在建立一個 width * height，標題為「CUDA volume rendering」的視窗。而之後的glut*Func() 則是在設定不同的 callback function。而之後，則是做 glew 的起始化，並確認目前的 OpenGL 環境，是否有支援必要的 pixel buffer object。

而第四部分，就是在 initPixelBuffer() 這個函式中，建立用來當輸出結果的 pixel buffer object。

if (pbo) { // delete old buffer CUDA_SAFE_CALL(cudaGLUnregisterBufferObject(pbo)); glDeleteBuffersARB(1, &pbo); } // create pixel buffer object for display glGenBuffersARB(1, &pbo); glBindBufferARB(GL_PIXEL_UNPACK_BUFFER_ARB, pbo); glBufferDataARB(GL_PIXEL_UNPACK_BUFFER_ARB, width*height*sizeof(GLubyte)*4, 0, GL_STREAM_DRAW_ARB); glBindBufferARB(GL_PIXEL_UNPACK_BUFFER_ARB, 0); CUDA_SAFE_CALL(cudaGLRegisterBufferObject(pbo)); // calculate new grid size gridSize = dim3(iDivUp(width, blockSize.x), iDivUp(height, blockSize.y));

第一段的 if 是判斷 pbo 這個 buffer object（型別是 GLuint）是否已經被建立，如果有的話，就先把現有的刪除掉。而第二段就是透過 OpenGL 的函式，來建立大小為width * height，每一點的資料是 4 個 GLubyte 的一塊GL_PIXEL_UNPACK_BUFFER_ARB 了～

而在 CUDA 的部分，則是需要透過 cudaGLRegisterBufferObject() 這個函式來註冊pbo 這個 buffer object；如此，之後才能在 kernel 程式中存取透過 cudaGLMapBufferObject() 所取得的記憶體位址。

當使用完後，如果在 if 的區段中所做，除了必須要透過 OpenGL 的 glDeleteBuffersARB() 來把 buffer object 刪除外，在之前也需要使用cudaGLUnegisterBufferObject() 來取消這份 buffer object 的註冊。

而在 main() 中最後的 glutMainLoop()，就是開始執行 OpenGL 的 main loop；之後，要顯示內容、控制程式，就是要靠之前設定的 callback function 囉～

 

 

而接下來這一部分，主要就來看顯示的部分了～

在這個程式中，要去把 volume 畫出來的部分，主要就是在 OpenGL 的 display callback function，也就是 display() 這個函式；當每次需要 render 的時候，就會呼叫到這個函式。而整個繪製的流程，大致上如下圖所示：

上圖中，左邊的區域，就是 display() 的大致流程。中間「Render to pixel buffer object」的大框，則是render() 這個函式在做的事；CUDA 的 kernel 程式 d_Render()，也是在 render() 中被呼叫的。右邊紅色的區域，則是代表 CUDA 的 kerneld_Render() 的流程。

接下來，就先來看 display() 的程式內容：

// use OpenGL to build view matrix GLfloat modelView[16]; glMatrixMode(GL_MODELVIEW); glPushMatrix(); glLoadIdentity(); glRotatef(-viewRotation.x, 1.0, 0.0, 0.0); glRotatef(-viewRotation.y, 0.0, 1.0, 0.0); glTranslatef(-viewTranslation.x, -viewTranslation.y, -viewTranslation.z); glGetFloatv(GL_MODELVIEW_MATRIX, modelView); glPopMatrix(); invViewMatrix[0] = modelView[0]; invViewMatrix[1] = modelView[4]; invViewMatrix[2] = modelView[8]; invViewMatrix[3] = modelView[12]; invViewMatrix[4] = modelView[1]; invViewMatrix[5] = modelView[5]; invViewMatrix[6] = modelView[9]; invViewMatrix[7] = modelView[13]; invViewMatrix[8] = modelView[2]; invViewMatrix[9] = modelView[6]; invViewMatrix[10] = modelView[10]; invViewMatrix[11] = modelView[14]; render(); // display results glClear(GL_COLOR_BUFFER_BIT); // draw image from PBO glDisable(GL_DEPTH_TEST); glRasterPos2i(0, 0); glBindBufferARB(GL_PIXEL_UNPACK_BUFFER_ARB, pbo); glDrawPixels(width, height, GL_RGBA, GL_UNSIGNED_BYTE, 0); glBindBufferARB(GL_PIXEL_UNPACK_BUFFER_ARB, 0); glutSwapBuffers(); glutReportErrors();

如同前面所述，display() 的內容分成三部分：

1. 處理攝影機的位置矩陣
2. 呼叫 render() 進行繪製
3. 繪製 Pixel Buffer Object

第一部分矩陣計算的部分，基本上這個範例程式在 motion() 中，可以透過滑鼠來控制攝影機（眼睛的位置）的旋轉和平移，並把資料存在viewRotation 和 viewTranslation 裡；所以在進行 ray casting 前，要先算出來這些影響所形成的矩陣（viewTranslation 在 z 軸預設位移是 -4）。而他的方法是透過 OpenGL 的 model view matrix 來計算，並把結果讀到modelView 中，然後進行一個轉置，轉換成 kernel 中要用的形式（三乘四的矩陣），也就是 invViewMatrix。

而接下來的第二部分，是把最主要、也就是呼叫 ray casting 的 CUDA kernel 程式，都包到 render() 這個函式；裡面的內容有四步：

1. 將矩陣複製到 device constant memory
2. 將 pixel buffer object pbo 對應到 device memory d_output
3. 呼叫 kernel d_Render()
4. 解除 pbo 和 d_output 的對應

其程式碼如下：

CUDA_SAFE_CALL( cudaMemcpyToSymbol(c_invViewMatrix, invViewMatrix, sizeof(float4)*3) ); // map PBO to get CUDA device pointer uint *d_output; CUDA_SAFE_CALL(cudaGLMapBufferObject((void**)&d_output, pbo)); CUDA_SAFE_CALL(cudaMemset(d_output, 0, width*height*4)); // call CUDA kernel, writing results to PBO d_render<<<gridsize , blockSize>>>(d_output, width, height,density, brightness, transferOffset, transferScale); CUT_CHECK_ERROR("kernel failed"); CUDA_SAFE_CALL(cudaGLUnmapBufferObject(pbo));

因為在 ray casting 時要用到 invViewMatrix 這個矩陣的內容，而且是要在 kernel 中讀取，所以要先透過 cudaMemcpyToSymbol() 把資料複製到 c_invViewMatrix 這個 constant 變數。而c_invViewMatrix 是一個自訂的 structure，他的定義和宣告如下：

typedef struct { float4 m[3]; } float3x4; __constant__ float3x4 c_invViewMatrix; // inverse view matrix  

接著，就是透過 cudaGLMapBufferObject() 這個函式，來取得之前在part.2 時所建立的 pixel buffer object pbo 的記憶體位址：d_output；並再透過cudaMemset() 把他的資料都填 0。

再來就是呼叫 kernel 函式 d_render() 來進行最重要的 ray casting 的動作了！這個 kernel 的參數包含了要寫入結果的 PBO：d_output，以及 pixel buffer object 的寬和高（width,height）；而剩下的 density、brightness、transferOffset、transferScale 這四個變數，主要是用來控制呈現的效果的。關於 kernel 的部分，會在下一部分的 part.4 再來講。

在 kernel 結束後，ray casting 的結果也都會存到 d_output 裡了～而由於此時已經不會再對 pixel buffer object 做寫入的動作，所以也需要透過cudaGLUnmapBufferObject()，解除 pbo 和d_output 之間的對應。

而當 render() 的動作都結束後，就要回到 display() 執行最後一的動作：把 pixel buffer object 畫出來了～基本上，這邊就全都是 OpenGL 的東西了。主要，就是先透過glBindBufferARB() 來指定要使用的 buffer object（這邊當然是pbo），然後再透過 glDrawPixels() 把他畫出來了～

 

前面已經把 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 texturetransferTex 中，去找到對應的色彩；不過，這邊會先把取得的 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 範例，大概也就先講到這了～ :)