袁 斌

（北京應(yīng)用物理與計算數(shù)學(xué)研究所，北京 100088）

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

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

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

1 相關(guān)工作

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

2 基于曲率GPU光線投射

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

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

2.1 基于B樣條濾波的按需計算導(dǎo)數(shù)的方法

在線性濾波條件，采用近似的導(dǎo)數(shù)計算公式，而在B樣條濾波條件下，存在精確的導(dǎo)數(shù)計算公式。當(dāng)用B樣條濾波重構(gòu)數(shù)據(jù)場時，本文采用 C.Sigg方法快速計算標(biāo)量和導(dǎo)數(shù)[10]。這里分析一下C.Sigg的方法的鄰點(diǎn)數(shù)并進(jìn)行優(yōu)化處理。

如圖1所示，計算標(biāo)量需要8個鄰點(diǎn)，計算梯度需要24個鄰點(diǎn)，計算二階偏導(dǎo)數(shù)需要24個鄰點(diǎn)，計算二階導(dǎo)數(shù)需要計算36個鄰點(diǎn)。用Tex指令取出鄰點(diǎn)上的標(biāo)量值。

圖1 C.Sigg 方法

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

2.2 基于線性濾波的按需計算導(dǎo)數(shù)的方法

目前的 GPU已經(jīng)實現(xiàn)基于線性濾波的紋理映射，沒有實現(xiàn)基于高次濾波的紋理映射。在線性濾波條件下，近似計算各種導(dǎo)數(shù)。

圖2 領(lǐng)域和鄰點(diǎn)示意圖

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

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

2.3 曲率計算方法

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

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

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

2.4 曲率GPU光線投射的實現(xiàn)

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

1）把數(shù)據(jù)場裝入GPU內(nèi)存；

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

4）繪制數(shù)據(jù)體的后表面到FBO；

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

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

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

F為關(guān)于主曲率的函數(shù)，向量v的x,y分量分別保存當(dāng)前積分步兩端的標(biāo)量值?；谇实腉PU光線投射算法如下：

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

把當(dāng)前紋理坐標(biāo)(即當(dāng)前光線與前面的交點(diǎn))保存到tex0;

置光線累積透明度為0;

從FBO中取出背面上的紋理坐標(biāo)，即當(dāng)前光線與數(shù)據(jù)體后表面的交點(diǎn);

計算光線方向;

修正采樣步長;

計算當(dāng)前光線總的積分步數(shù);

計算單步增量;

LOOP (255,0,1){

LOOP (255,0,1){

將當(dāng)前采樣點(diǎn)的坐標(biāo)變換到基本紋理空間;

計算當(dāng)前采樣點(diǎn)的標(biāo)量值scal;

if(當(dāng)前采樣點(diǎn)是第一個采樣點(diǎn))v.x = scal;

否則{

當(dāng)前積分步終點(diǎn)標(biāo)量值v.y = scal;

用v查預(yù)積分表得到不透明度;

v.x = scal;

如果當(dāng)前積分步不透明度非0{

用SIMD方法計算標(biāo)量、法向量、

偏導(dǎo)數(shù)和二階導(dǎo)數(shù);

計算法向量的模;

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

計算F;

規(guī)格化梯度量；

根據(jù)梯度量查調(diào)制因子;

將F轉(zhuǎn)換為顏色;

用當(dāng)前積分步的不透明度預(yù)乘顏色；

計算光照；

用當(dāng)前積分步的不透明度修正鏡面光；

用梯度量調(diào)制環(huán)境光、漫反射光和鏡面光；

累積顏色；

調(diào)制當(dāng)前積分步的不透明度；

累積光線段的透明度；

如果累積透明度小于0.02,退出循環(huán)；

}

}

前進(jìn)一步;

如果累積步數(shù)等于總步數(shù)，退出循環(huán)；

}

}

計算光線不透明度；

輸出顏色和不透明度；

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

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

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

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;;// 計算主曲率之和;

我們已經(jīng)基于 VTK[13]實現(xiàn)了一個可視化原型系統(tǒng)，將本文算法集成到原型系統(tǒng)中進(jìn)行測試。在原型系統(tǒng)中，實現(xiàn)了修改轉(zhuǎn)換函數(shù)以及保存、恢復(fù)交互參數(shù)的功能。

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

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

3 實驗與分析

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

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

圖5 foot 數(shù)據(jù)256×256×256，高斯曲率Fg

如圖6(a)，圖6(b)所示，繪制質(zhì)量相同?；贐樣條濾波F0函數(shù)可以很好表現(xiàn)松針。如圖6(c)所示，基于線性濾波，F(xiàn)0也可以表現(xiàn)松針，但質(zhì)量較差。F0可以用于提取線狀特征。

圖6 Xmas-Tree 數(shù)據(jù) 512×512×999，F(xiàn)0

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

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

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

圖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]. 計算機(jī)學(xué)報,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] 鄭 杰,姬紅兵. 基于動態(tài)紋理載入的大規(guī)模數(shù)據(jù)場體繪制[J]. 中國圖象圖形學(xué)報,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] 袁 斌. 改進(jìn)的均勻數(shù)據(jù)場 GPU 光線投射[J]. 中國圖象圖形學(xué)報,2011,16(7): 1269-1275.

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