基于曲率的GPU光線投射

袁 斌

（北京應用物理與計算數學研究所，北京 100088）

基于曲率的繪制十分有用，它可以幫助科學家分析曲率分布，設計新的計算模型。曲率可以用來表現線狀特征和粒狀特征。K. Gordon等在CPU中實現了基于B樣條濾波的曲率計算[1]，但計算速度較慢。采用 GPU可以大大加速曲率的計算。

現代 GPU已經發展成為眾核處理器，具有強大的計算能力和高內存帶寬。在一般的圖形卡上，采用多條繪制流水線，可以對頂點(vertex)計算和片元(fragment)計算進行并行處理。采用GPU可以大大加速體繪制，如基于切片的體繪制方法[2-4]和GPU的光線投射方法[5-7]。預先計算一階導數、二階導數，需要占用大量的存儲空間（存儲梯度和Hessian矩陣），這顯然是不能接受的。因此需要實時計算他們。

本文在GPU上實現了基于線性濾波的曲率計算方法和并將C.Sigg方法結合到GPU光線投射方法中，按需計算導數，曲率和光照效果。在線性濾波情況下，一階導數和二階導數無法精確計算，只能采用差分方法。在三三次B樣條濾波條件下，可以精確計算標量，一節導數、二階導數。本文給出有效的計算曲率的方法。設計多種關于曲率的函數，實現了交互選擇曲率函數，突出感興趣特征的方法。

1 相關工作

S. Stegmaier給出單趟GPU光線投射方法[5]。該方法預先計算梯度并保存起來，這樣必然占用大量的內存。在透視投影下，光柵化時如果直接采用線性插值，就不能正確計算片元的顏色、紋理坐標等屬性，為此OpenGL對插值算法作了修正,這樣它可以正確計算片元的屬性[8]，這是基于代理面的GPU光線投射的基礎。T. Klein給出結合預積分的GPU光線投射方法[6]，通過繪制體的前表面得到光線的起點，起點減去相機的位置（在紋理空間）得到光線的方向。Scharsach的論文[7]把體包圍盒頂點坐標轉換成顏色，并設置為頂點顏色。繪制前面和后面到不同的緩沖區，后面的圖像減前面的圖像，得到光線方向，并根據這個方向進行光線積分。該文還設計基于塊的空區域跳躍方法，通過深度測試，把第一前面和最后一個后面繪制到不同的緩沖區中，這樣在光線積分時可以跳過前面和的后面的空區域，提高了光線投射的速度。鄭杰等在紋理切片體繪制方法中，采用實時計算梯度的方法[9]，大大節省紋理內存，適合于大規模數據場，是很好的方法。K.Gordon等[1]在CPU中實現了基于B樣條一階導數、二階導數以及曲率計算，但計算速度較慢。C.Sigg等[10]在GPU中實現基于B樣條的一階導數、二階導數快速計算，并應用到最大主曲率的計算；并把該方法與等值面繪制結合，繪制速度很快；但未把該方法與光線投射體繪制結合。

2 基于曲率GPU光線投射

本文實現基于曲率的可視化，要計算隱式曲面的曲率，需要計算Hessian矩陣和梯度。

一階導數，二階導數計算采用實時計算方法。如果預先計算各種導數，每個網格節點需要保存3個一階導數，3個二階導數，和3個二階偏導數，需要多存儲9個數據，若采用單字節表示，另外需要 1個字節保存梯度量，則需要 10個字節，存儲量將是本文方法的 11倍。采用這種方法，精度很低，會降低繪制質量。若采用浮點數表示標量和導數，存儲量為本文方法的 37倍。對較大的數據，如 Xmas-Tree，不可能裝入到紋理內存。因此，不能采用預計算一階導數，二階導數的方法。

2.1 基于B樣條濾波的按需計算導數的方法

在線性濾波條件，采用近似的導數計算公式，而在B樣條濾波條件下，存在精確的導數計算公式。當用B樣條濾波重構數據場時，本文采用 C.Sigg方法快速計算標量和導數[10]。這里分析一下C.Sigg的方法的鄰點數并進行優化處理。

如圖1所示，計算標量需要8個鄰點，計算梯度需要24個鄰點，計算二階偏導數需要24個鄰點，計算二階導數需要計算36個鄰點。用Tex指令取出鄰點上的標量值。

圖1 C.Sigg 方法

計算導數和曲率需要大量的計算，本文按需計算導數，曲率和光照。這樣可以大大加速計算。在GPU光線投射算法中,用C.Sigg方法計算標量值，需要用8個Tex指令取出當前采樣的（8個）鄰點上標量值，結合前一采樣點的標量值，得到當前積分步的不透明度，如果不透明度為非0，用 C.Sigg方法計算計算導數，需要用 84個Tex指令取出（84個）鄰點上標量值，然后計算曲率以及光照效果。這樣有效地加速繪制過程。

2.2 基于線性濾波的按需計算導數的方法

目前的 GPU已經實現基于線性濾波的紋理映射，沒有實現基于高次濾波的紋理映射。在線性濾波條件下，近似計算各種導數。

圖2 領域和鄰點示意圖

考慮以采樣點為中心的2×2×2軸向長方形網格D，其網格間隔與被處理的網格一致，D被稱為采樣點的鄰域，其節點數為3×3×3，節點上的物理量記為其中f111等于光線段上當前采樣點的值。因為光線段上的采樣點一般不在網格節點上，所以，D的節點一般情況下，位于原網格單元的內部，而不是在節點上。采用中心差分計算一階導數和二階導數需要18個鄰點，以光線采樣點為中心，取18個鄰點,如圖2(b)。鄰點的物理值通過三線性插值計算，通過紋理映射實現。采用向量指令計算梯度和Hessian矩陣，用一個向量保存對角元素( fxx,fyy,fzz)，另一個向量保存( fxy,fyz,fzx)

在光線投射算法中，先用1個指令取出當前的標量值，結合前一采樣點的標量值，得到當前積分步的不透明度，如果不透明度非0，采用18鄰點方法，計算梯度，Hessian矩陣和曲率，是一個實用的方法。

2.3 曲率計算方法

稱為總曲率. 將 trace(P HP(P HP )T)直接展開比較麻煩。采用矩陣計算會簡捷一些，但仍然需要采用較多步數，因此，不同于文獻[1]，本文不直接計算，而是先計算高斯曲率，下面給出高斯曲率公式的簡單推導過程

這樣，可以用向量乘法和DP3指令實現以上公式[11]。導數、二階導數和曲率均在物理空間計算，而不是在圖像空間計算。在計算出b和c后，就可以計算主曲率。

主曲率可以描述曲面的局部特征，根據主曲率，曲面上的點可以分為橢圓點，雙曲點，和拋物點。曲面的幾何特征如“溝谷”，“山脊”，“山峰”，“洼坑”，線狀特征和粒狀特征等與曲率密切相關。根據需要設計關于主曲率的函數F(κ1,κ2)，來突出感興趣的特征。

2.4 曲率GPU光線投射的實現

實現基于代理面的 GPU光線投射體繪制的基本步驟是：

1）把數據場裝入GPU內存；

2）裝入寫緩沖區和光線積分片元程序；3）連接寫緩沖區的片元程序；

4）繪制數據體的后表面到FBO；

5）連接GPU光線投射的片元程序；

6）設置GPU光線投射程序的局部參數；7）繪制數據體的前表面。

開始時，實現基于B樣條濾波的曲率GPU光線投射，未考慮按需計算導數，繪制速度較慢。為了加速繪制，采用按需計算導數的方法。

F為關于主曲率的函數，向量v的x,y分量分別保存當前積分步兩端的標量值。基于曲率的GPU光線投射算法如下：

算法1 基于曲率的光線投射算法

把當前紋理坐標(即當前光線與前面的交點)保存到tex0;

置光線累積透明度為0;

從FBO中取出背面上的紋理坐標，即當前光線與數據體后表面的交點;

計算光線方向;

修正采樣步長;

計算當前光線總的積分步數;

計算單步增量;

LOOP (255,0,1){

LOOP (255,0,1){

將當前采樣點的坐標變換到基本紋理空間;

計算當前采樣點的標量值scal;

if(當前采樣點是第一個采樣點)v.x = scal;

否則{

當前積分步終點標量值v.y = scal;

用v查預積分表得到不透明度;

v.x = scal;

如果當前積分步不透明度非0{

用SIMD方法計算標量、法向量、

偏導數和二階導數;

計算法向量的模;

采用SIMD方法計算高斯曲率和主曲率之和；

計算F;

規格化梯度量；

根據梯度量查調制因子;

將F轉換為顏色;

用當前積分步的不透明度預乘顏色；

計算光照；

用當前積分步的不透明度修正鏡面光；

用梯度量調制環境光、漫反射光和鏡面光；

累積顏色；

調制當前積分步的不透明度；

累積光線段的透明度；

如果累積透明度小于0.02,退出循環；

}

}

前進一步;

如果累積步數等于總步數，退出循環；

}

}

計算光線不透明度；

輸出顏色和不透明度；

注：LOOP (255,0,1)表示循環255次。

在GPU光線投射的片元程序中實現算法1，算法采用按需計算導數的方法。根據當前采樣點的標量和前一采樣點的標量，查預積分表取得當前積分步的不透明度，只有不透明度非0時，才計算梯度、Hessian矩陣、曲率和光照。計算標量時，如果采用線性濾波，用一個Tex指令取出當前采樣點的標量值[11]；若采用B樣條濾波，用8個Tex指令取出8個鄰點。在計算導數時，若采用線性濾波需要18個鄰點；若采用B樣條濾波，需要 84個鄰點。顏色轉換函數把關于主曲率的函數F(κ1,κ2)轉換成顏色。顏色轉換函數和預積分表用紋理實現。用 GPU快速計算預積分[12]，這樣可以交互修改轉換函數。函數F的值緊縮到[0.02,0.98]，如果它大于 0.98,置為 0.98;如果小于0.02,置為0.02。

本算法充分利用SIMD技術進行計算。采用光線早中止加速繪制過程。在本文中，采用SIMD方法計算高斯曲率和主曲率之和，只需要 17條指令。計算高斯曲率的一些中間結果可以用來計算主曲率之和。temp3保存二階導數，temp5保存二階偏導，normal保存法向量。temp2.x為法向量模的平方。下面是相關的程序段。

MUL temp8,normal,normal;

ADD temp9,temp8.yzxw,temp8.zxyw;

DP3 temp9.x,temp9,temp3;

MUL temp7,temp3,temp3.yzxw;

DP3 temp4.x,temp7,temp8.zxyw;

MUL temp7,temp5,temp5;

DP3 temp4.y,temp7,temp8.zxyw;

MUL temp8,normal.yzxw,normal.zxyw;

DP3 temp9.y,temp5,temp8.zxyw;

MUL temp7,temp5,temp5.zxyw;

DP3 temp4.z,temp7,temp8;

MUL temp7,temp3,temp5.yzxw;

DP3 temp4.w,temp7,temp8;

DP4 temp6,temp4,{1,-1,2,-2};

DIV temp6,temp6,temp2.x;// 計算高斯曲率;

DP2 temp10,temp9,{1,-2,0,0};

MUL temp10,temp10,temp1.x;;// 計算主曲率之和;

我們已經基于 VTK[13]實現了一個可視化原型系統，將本文算法集成到原型系統中進行測試。在原型系統中，實現了修改轉換函數以及保存、恢復交互參數的功能。

在GPU光線投射中，采用高次B-樣條濾波可以起到光順作用，消除走樣現象，繪制質量很好。但是由于采用B樣條濾波，計算量較大，繪制速度較慢。本文同時實現了基于線性濾波和B樣條濾波的按需計算導數的GPU光線投射方法。在原型系統中，增加了選擇曲率函數的交互功能；修改保存和恢復功能，使之包括曲率函數。

在交互過程中，先采用基于線性濾波的繪制方法交互選擇必要的相機參數、轉換函數和曲率函數，再用基于B樣條濾波的方法繪制高質量圖像。

3 實驗與分析

本文在GPU上實現了基于線性濾波和B樣條濾波的按需計算導數、曲率和光照效果的光線投射方法，分別記為LRT和BRT,為了便于比較，還實現了基于B樣條濾波的非按需計算導數（曲率和光照效果按需計算）的光線投射方法，記為BRTN。本文在裝有NVIDIA NVS 4200M 顯卡的intel i5機器（便攜式電腦）上測試算法。以下測試中，除非有特殊聲明，采樣步長為最小網格間隔的1/4光線積分計算中，采用光線早中止技術，不透明度決定了實際的采樣步數，顏色轉換函數不影響采樣步數。在測試繪制速度時，對一種數據，選擇一個曲率函數即可。在測試過程中，本文利用基于線性濾波的光線投射方法 LRT交互選擇相機設置、各種轉換函數和曲率函數，突出感興趣的特征，然后用基于 B樣條濾波的方法BRT繪制高質量的圖像。

圖3(a)與圖3(b)繪制質量相同，圖3(b)速度更快。B樣條濾波可以更好地突出特征，如圖3(b)所示，而線性濾波淡化了特征，如圖3(c)所示。從圖3可以看出，+0.4可以較好地表現線狀特征（綠色）。圖4為head數據，采樣步長為最小網格間隔的1/16，采用B樣條濾波，管狀特征很光滑；而采用線性濾波，圖像質量可以接受。圖5(a)與圖5(b)繪制質量相同，圖5(b)速度更快。B樣條濾波可以更好地突出特征，如圖5(b)所示，而線性濾波淡化了特征，如圖5(c)所示，適當調整顏色轉換函數可以增強特征。從圖5可以看出，高斯曲率可以較好地表現粒狀特征（紅色或綠色）。這些粒狀特征是由噪聲引起的。高斯曲率可以分析噪聲的分布。

圖5 foot 數據256×256×256，高斯曲率Fg

如圖6(a)，圖6(b)所示，繪制質量相同。基于B樣條濾波F0函數可以很好表現松針。如圖6(c)所示，基于線性濾波，F0也可以表現松針，但質量較差。F0可以用于提取線狀特征。

圖6 Xmas-Tree 數據 512×512×999，F0

如圖 7(a)所示，平均曲率 Fm既可以表現線狀特征，又可以表現粒狀特征；圖7(b)采用總曲率Ft突出了線狀特征。

由表1可以看出，按需計算可以大大加快繪制速度，線性濾波快于B樣條濾波。采用線性濾波可以進行交互繪制。

LRT先用1個Tex指令取出當前的標量值，結合前一采樣點的標量值，得到當前積分步的不透明度，如果不透明度非0，采用18鄰點方法，計算梯度，Hessian矩陣和曲率，是一個實用的方法。BRT用C.Sigg方法計算標量值，需要用8個Tex指令取出當前采樣的(8個)鄰點上標量值，結合前一采樣點的標量值，得到當前積分步的不透明度，如果不透明度為非0，用C.Sigg方法計算計算導數，需要用84個Tex指令取出（84個）鄰點上標量值，然后計算曲率以及光照效果，從而有效地加速繪制過程。BRT繪制質量比基于LRT好，LRT速度比BRT快。本文的方法可以交互選擇曲率函數，突出感興趣的特征，從而幫助科學家有效地分析數據。

圖7 平均曲率和總曲率的繪制效果(BRT繪制)

表1 繪制時間比較

[1] Kindlmann,G,Whitaker R,Tasdizen T,et al. 2003 Curvature-based transfer functions for direct volume rendering: Methods and applications[C]//Proceedings of IEEE Visualization,2003: 513-520.

[2] Cabral B,Cam N,Foran J. Accelerated volume rendering and tomographic reconstruction using texture mapping hardware [C]//Proceedings of the 1994 Symposium on Volume Visualization,New York,NY,USA,ACM Press,1994: 91-98.

[3] Engel K,Kraus M,Ertl T. High-quality pre-integrated volume rendering using hard-ware-accelerated pixel shading [C]//Proceedings of the ACM SIGGRAPH/EUROGRAPHICS Workshop on Graphics hardware,2001: 9-16.

[4] 童 欣,唐澤圣. 基于空間跳躍的三維紋理硬件體繪制算法[J]. 計算機學報,1998,21(9): 807-812.

[5] Stegmaier S,Strengert M,Klein T,et al. A simple and flexible volume rendering framework for graphics-hardware-based raycasting [C]//Proceedings of the International Workshop on Volume Graphics'05,Stony Brook,New York,USA,2005: 187-195

[6] Klein T,Strengert M,Stegmaier S,et al. Exploiting frame-to-frame coherence for accelerating highquality volume raycasting on graphics hardware [C]//Proceedings of IEEE Visualization,2005: 223-230.

[7] Scharsach H. Advanced GPU raycasting [C]//Proceedings of Central European Seminar on Computer Graphics 2005,Budmerice Castle,2005:69-76.

[8] Segal M,Akeley K. The OpenGL graphics system: a specification [EB/OL]. 2004.http://www.opengl.org/documentation/specs/

[9] 鄭 杰,姬紅兵. 基于動態紋理載入的大規模數據場體繪制[J]. 中國圖象圖形學報,2008,13(2):316-321.

[10] Sigg C,Hadwiger M. Fast third-order texture filtering [C]//Book Chapter in GPU Gems II,Matt Pharr(ed.),Addison Wesley,2005: 313-329

[11] NVIDIA. NVIDIA OpenGL extension specifications. 2011. https://developer.nvidia.com/nvidia-opengl-specs

[12] 袁 斌. 改進的均勻數據場 GPU 光線投射[J]. 中國圖象圖形學報,2011,16(7): 1269-1275.

[13] Martin K,Schroeder W,Lorensen B. The Visualization ToolKit(VTK)[EB/OL]. 2008. http://www.vtk.org/