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 casting	Slice 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 列表如下：

main	main 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 函式，進行計算。
display	glut 的 callback function，用來顯示用的，每次要更新畫面，都是執行這個函式。 他會計算物體的旋轉矩陣，然後再呼叫 render() 來做 ray casting 計算，最後再把儲存了結果的 pixel buffer object 畫出來。
reshape	這四個 function 是 glut 的 callback function。 reshape 是當視窗大小位置改變時會被執行到、motion 是滑鼠移動時會被呼叫的、mouse 是滑鼠按鈕會執行的函式、keyboard 則是鍵盤的。 這邊不多加解釋了。
motion
mouse
keyboard

而在 volumeRender_kernel.cu 中，則是：

intersectBox	計算一條視線和 Volume 的 box 的交點，並傳回相交的最近點和最遠點。
mul	計算矩陣乘上一個向量。
rgbaFloatToInt	把 rgba 四項的 color，轉換成一個 int 來儲存。
d_render	CUDA 的 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()，把 cudaArray d_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 的 kernel d_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 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 範例，大概也就先講到這了～ :)