袁 斌
(北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所,北京 100088)
基于曲率的繪制十分有用,它可以幫助科學(xué)家分析曲率分布,設(shè)計(jì)新的計(jì)算模型。曲率可以用來(lái)表現(xiàn)線狀特征和粒狀特征。K. Gordon等在CPU中實(shí)現(xiàn)了基于B樣條濾波的曲率計(jì)算[1],但計(jì)算速度較慢。采用 GPU可以大大加速曲率的計(jì)算。
現(xiàn)代 GPU已經(jīng)發(fā)展成為眾核處理器,具有強(qiáng)大的計(jì)算能力和高內(nèi)存帶寬。在一般的圖形卡上,采用多條繪制流水線,可以對(duì)頂點(diǎn)(vertex)計(jì)算和片元(fragment)計(jì)算進(jìn)行并行處理。采用GPU可以大大加速體繪制,如基于切片的體繪制方法[2-4]和GPU的光線投射方法[5-7]。預(yù)先計(jì)算一階導(dǎo)數(shù)、二階導(dǎo)數(shù),需要占用大量的存儲(chǔ)空間(存儲(chǔ)梯度和Hessian矩陣),這顯然是不能接受的。因此需要實(shí)時(shí)計(jì)算他們。
本文在GPU上實(shí)現(xiàn)了基于線性濾波的曲率計(jì)算方法和并將C.Sigg方法結(jié)合到GPU光線投射方法中,按需計(jì)算導(dǎo)數(shù),曲率和光照效果。在線性濾波情況下,一階導(dǎo)數(shù)和二階導(dǎo)數(shù)無(wú)法精確計(jì)算,只能采用差分方法。在三三次B樣條濾波條件下,可以精確計(jì)算標(biāo)量,一節(jié)導(dǎo)數(shù)、二階導(dǎo)數(shù)。本文給出有效的計(jì)算曲率的方法。設(shè)計(jì)多種關(guān)于曲率的函數(shù),實(shí)現(xiàn)了交互選擇曲率函數(shù),突出感興趣特征的方法。
S. Stegmaier給出單趟GPU光線投射方法[5]。該方法預(yù)先計(jì)算梯度并保存起來(lái),這樣必然占用大量的內(nèi)存。在透視投影下,光柵化時(shí)如果直接采用線性插值,就不能正確計(jì)算片元的顏色、紋理坐標(biāo)等屬性,為此OpenGL對(duì)插值算法作了修正,這樣它可以正確計(jì)算片元的屬性[8],這是基于代理面的GPU光線投射的基礎(chǔ)。T. Klein給出結(jié)合預(yù)積分的GPU光線投射方法[6],通過(guò)繪制體的前表面得到光線的起點(diǎn),起點(diǎn)減去相機(jī)的位置(在紋理空間)得到光線的方向。Scharsach的論文[7]把體包圍盒頂點(diǎn)坐標(biāo)轉(zhuǎn)換成顏色,并設(shè)置為頂點(diǎn)顏色。繪制前面和后面到不同的緩沖區(qū),后面的圖像減前面的圖像,得到光線方向,并根據(jù)這個(gè)方向進(jìn)行光線積分。該文還設(shè)計(jì)基于塊的空區(qū)域跳躍方法,通過(guò)深度測(cè)試,把第一前面和最后一個(gè)后面繪制到不同的緩沖區(qū)中,這樣在光線積分時(shí)可以跳過(guò)前面和的后面的空區(qū)域,提高了光線投射的速度。鄭杰等在紋理切片體繪制方法中,采用實(shí)時(shí)計(jì)算梯度的方法[9],大大節(jié)省紋理內(nèi)存,適合于大規(guī)模數(shù)據(jù)場(chǎng),是很好的方法。K.Gordon等[1]在CPU中實(shí)現(xiàn)了基于B樣條一階導(dǎo)數(shù)、二階導(dǎo)數(shù)以及曲率計(jì)算,但計(jì)算速度較慢。C.Sigg等[10]在GPU中實(shí)現(xiàn)基于B樣條的一階導(dǎo)數(shù)、二階導(dǎo)數(shù)快速計(jì)算,并應(yīng)用到最大主曲率的計(jì)算;并把該方法與等值面繪制結(jié)合,繪制速度很快;但未把該方法與光線投射體繪制結(jié)合。
本文實(shí)現(xiàn)基于曲率的可視化,要計(jì)算隱式曲面的曲率,需要計(jì)算Hessian矩陣和梯度。
一階導(dǎo)數(shù),二階導(dǎo)數(shù)計(jì)算采用實(shí)時(shí)計(jì)算方法。如果預(yù)先計(jì)算各種導(dǎo)數(shù),每個(gè)網(wǎng)格節(jié)點(diǎn)需要保存3個(gè)一階導(dǎo)數(shù),3個(gè)二階導(dǎo)數(shù),和3個(gè)二階偏導(dǎo)數(shù),需要多存儲(chǔ)9個(gè)數(shù)據(jù),若采用單字節(jié)表示,另外需要 1個(gè)字節(jié)保存梯度量,則需要 10個(gè)字節(jié),存儲(chǔ)量將是本文方法的 11倍。采用這種方法,精度很低,會(huì)降低繪制質(zhì)量。若采用浮點(diǎn)數(shù)表示標(biāo)量和導(dǎo)數(shù),存儲(chǔ)量為本文方法的 37倍。對(duì)較大的數(shù)據(jù),如 Xmas-Tree,不可能裝入到紋理內(nèi)存。因此,不能采用預(yù)計(jì)算一階導(dǎo)數(shù),二階導(dǎo)數(shù)的方法。
在線性濾波條件,采用近似的導(dǎo)數(shù)計(jì)算公式,而在B樣條濾波條件下,存在精確的導(dǎo)數(shù)計(jì)算公式。當(dāng)用B樣條濾波重構(gòu)數(shù)據(jù)場(chǎng)時(shí),本文采用 C.Sigg方法快速計(jì)算標(biāo)量和導(dǎo)數(shù)[10]。這里分析一下C.Sigg的方法的鄰點(diǎn)數(shù)并進(jìn)行優(yōu)化處理。
如圖1所示,計(jì)算標(biāo)量需要8個(gè)鄰點(diǎn),計(jì)算梯度需要24個(gè)鄰點(diǎn),計(jì)算二階偏導(dǎo)數(shù)需要24個(gè)鄰點(diǎn),計(jì)算二階導(dǎo)數(shù)需要計(jì)算36個(gè)鄰點(diǎn)。用Tex指令取出鄰點(diǎn)上的標(biāo)量值。
圖1 C.Sigg 方法
計(jì)算導(dǎo)數(shù)和曲率需要大量的計(jì)算,本文按需計(jì)算導(dǎo)數(shù),曲率和光照。這樣可以大大加速計(jì)算。在GPU光線投射算法中,用C.Sigg方法計(jì)算標(biāo)量值,需要用8個(gè)Tex指令取出當(dāng)前采樣的(8個(gè))鄰點(diǎn)上標(biāo)量值,結(jié)合前一采樣點(diǎn)的標(biāo)量值,得到當(dāng)前積分步的不透明度,如果不透明度為非0,用 C.Sigg方法計(jì)算計(jì)算導(dǎo)數(shù),需要用 84個(gè)Tex指令取出(84個(gè))鄰點(diǎn)上標(biāo)量值,然后計(jì)算曲率以及光照效果。這樣有效地加速繪制過(guò)程。
目前的 GPU已經(jīng)實(shí)現(xiàn)基于線性濾波的紋理映射,沒(méi)有實(shí)現(xiàn)基于高次濾波的紋理映射。在線性濾波條件下,近似計(jì)算各種導(dǎo)數(shù)。
圖2 領(lǐng)域和鄰點(diǎn)示意圖
考慮以采樣點(diǎn)為中心的2×2×2軸向長(zhǎng)方形網(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)的值。因?yàn)楣饩€段上的采樣點(diǎn)一般不在網(wǎng)格節(jié)點(diǎn)上,所以,D的節(jié)點(diǎn)一般情況下,位于原網(wǎng)格單元的內(nèi)部,而不是在節(jié)點(diǎn)上。采用中心差分計(jì)算一階導(dǎo)數(shù)和二階導(dǎo)數(shù)需要18個(gè)鄰點(diǎn),以光線采樣點(diǎn)為中心,取18個(gè)鄰點(diǎn),如圖2(b)。鄰點(diǎn)的物理值通過(guò)三線性插值計(jì)算,通過(guò)紋理映射實(shí)現(xiàn)。采用向量指令計(jì)算梯度和Hessian矩陣,用一個(gè)向量保存對(duì)角元素( fxx,fyy,fzz),另一個(gè)向量保存( fxy,fyz,fzx)
在光線投射算法中,先用1個(gè)指令取出當(dāng)前的標(biāo)量值,結(jié)合前一采樣點(diǎn)的標(biāo)量值,得到當(dāng)前積分步的不透明度,如果不透明度非0,采用18鄰點(diǎn)方法,計(jì)算梯度,Hessian矩陣和曲率,是一個(gè)實(shí)用的方法。
稱為總曲率. 將 trace(P HP(P HP )T)直接展開(kāi)比較麻煩。采用矩陣計(jì)算會(huì)簡(jiǎn)捷一些,但仍然需要采用較多步數(shù),因此,不同于文獻(xiàn)[1],本文不直接計(jì)算,而是先計(jì)算高斯曲率,下面給出高斯曲率公式的簡(jiǎn)單推導(dǎo)過(guò)程
這樣,可以用向量乘法和DP3指令實(shí)現(xiàn)以上公式[11]。導(dǎo)數(shù)、二階導(dǎo)數(shù)和曲率均在物理空間計(jì)算,而不是在圖像空間計(jì)算。在計(jì)算出b和c后,就可以計(jì)算主曲率。
主曲率可以描述曲面的局部特征,根據(jù)主曲率,曲面上的點(diǎn)可以分為橢圓點(diǎn),雙曲點(diǎn),和拋物點(diǎn)。曲面的幾何特征如“溝谷”,“山脊”,“山峰”,“洼坑”,線狀特征和粒狀特征等與曲率密切相關(guān)。根據(jù)需要設(shè)計(jì)關(guān)于主曲率的函數(shù)F(κ1,κ2),來(lái)突出感興趣的特征。
實(shí)現(xiàn)基于代理面的 GPU光線投射體繪制的基本步驟是:
1)把數(shù)據(jù)場(chǎng)裝入GPU內(nèi)存;
2)裝入寫(xiě)緩沖區(qū)和光線積分片元程序;3)連接寫(xiě)緩沖區(qū)的片元程序;
4)繪制數(shù)據(jù)體的后表面到FBO;
5)連接GPU光線投射的片元程序;
6)設(shè)置GPU光線投射程序的局部參數(shù);7)繪制數(shù)據(jù)體的前表面。
開(kāi)始時(shí),實(shí)現(xiàn)基于B樣條濾波的曲率GPU光線投射,未考慮按需計(jì)算導(dǎo)數(shù),繪制速度較慢。為了加速繪制,采用按需計(jì)算導(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);
計(jì)算光線方向;
修正采樣步長(zhǎng);
計(jì)算當(dāng)前光線總的積分步數(shù);
計(jì)算單步增量;
LOOP (255,0,1){
LOOP (255,0,1){
將當(dāng)前采樣點(diǎn)的坐標(biāo)變換到基本紋理空間;
計(jì)算當(dāng)前采樣點(diǎn)的標(biāo)量值scal;
if(當(dāng)前采樣點(diǎn)是第一個(gè)采樣點(diǎn))v.x = scal;
否則{
當(dāng)前積分步終點(diǎn)標(biāo)量值v.y = scal;
用v查預(yù)積分表得到不透明度;
v.x = scal;
如果當(dāng)前積分步不透明度非0{
用SIMD方法計(jì)算標(biāo)量、法向量、
偏導(dǎo)數(shù)和二階導(dǎo)數(shù);
計(jì)算法向量的模;
采用SIMD方法計(jì)算高斯曲率和主曲率之和;
計(jì)算F;
規(guī)格化梯度量;
根據(jù)梯度量查調(diào)制因子;
將F轉(zhuǎn)換為顏色;
用當(dāng)前積分步的不透明度預(yù)乘顏色;
計(jì)算光照;
用當(dāng)前積分步的不透明度修正鏡面光;
用梯度量調(diào)制環(huán)境光、漫反射光和鏡面光;
累積顏色;
調(diào)制當(dāng)前積分步的不透明度;
累積光線段的透明度;
如果累積透明度小于0.02,退出循環(huán);
}
}
前進(jìn)一步;
如果累積步數(shù)等于總步數(shù),退出循環(huán);
}
}
計(jì)算光線不透明度;
輸出顏色和不透明度;
注:LOOP (255,0,1)表示循環(huán)255次。
在GPU光線投射的片元程序中實(shí)現(xiàn)算法1,算法采用按需計(jì)算導(dǎo)數(shù)的方法。根據(jù)當(dāng)前采樣點(diǎn)的標(biāo)量和前一采樣點(diǎn)的標(biāo)量,查預(yù)積分表取得當(dāng)前積分步的不透明度,只有不透明度非0時(shí),才計(jì)算梯度、Hessian矩陣、曲率和光照。計(jì)算標(biāo)量時(shí),如果采用線性濾波,用一個(gè)Tex指令取出當(dāng)前采樣點(diǎn)的標(biāo)量值[11];若采用B樣條濾波,用8個(gè)Tex指令取出8個(gè)鄰點(diǎn)。在計(jì)算導(dǎo)數(shù)時(shí),若采用線性濾波需要18個(gè)鄰點(diǎn);若采用B樣條濾波,需要 84個(gè)鄰點(diǎn)。顏色轉(zhuǎn)換函數(shù)把關(guān)于主曲率的函數(shù)F(κ1,κ2)轉(zhuǎn)換成顏色。顏色轉(zhuǎn)換函數(shù)和預(yù)積分表用紋理實(shí)現(xiàn)。用 GPU快速計(jì)算預(yù)積分[12],這樣可以交互修改轉(zhuǎn)換函數(shù)。函數(shù)F的值緊縮到[0.02,0.98],如果它大于 0.98,置為 0.98;如果小于0.02,置為0.02。
本算法充分利用SIMD技術(shù)進(jìn)行計(jì)算。采用光線早中止加速繪制過(guò)程。在本文中,采用SIMD方法計(jì)算高斯曲率和主曲率之和,只需要 17條指令。計(jì)算高斯曲率的一些中間結(jié)果可以用來(lái)計(jì)算主曲率之和。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;// 計(jì)算高斯曲率;
DP2 temp10,temp9,{1,-2,0,0};
MUL temp10,temp10,temp1.x;;// 計(jì)算主曲率之和;
我們已經(jīng)基于 VTK[13]實(shí)現(xiàn)了一個(gè)可視化原型系統(tǒng),將本文算法集成到原型系統(tǒng)中進(jìn)行測(cè)試。在原型系統(tǒng)中,實(shí)現(xiàn)了修改轉(zhuǎn)換函數(shù)以及保存、恢復(fù)交互參數(shù)的功能。
在GPU光線投射中,采用高次B-樣條濾波可以起到光順作用,消除走樣現(xiàn)象,繪制質(zhì)量很好。但是由于采用B樣條濾波,計(jì)算量較大,繪制速度較慢。本文同時(shí)實(shí)現(xiàn)了基于線性濾波和B樣條濾波的按需計(jì)算導(dǎo)數(shù)的GPU光線投射方法。在原型系統(tǒng)中,增加了選擇曲率函數(shù)的交互功能;修改保存和恢復(fù)功能,使之包括曲率函數(shù)。
在交互過(guò)程中,先采用基于線性濾波的繪制方法交互選擇必要的相機(jī)參數(shù)、轉(zhuǎn)換函數(shù)和曲率函數(shù),再用基于B樣條濾波的方法繪制高質(zhì)量圖像。
本文在GPU上實(shí)現(xiàn)了基于線性濾波和B樣條濾波的按需計(jì)算導(dǎo)數(shù)、曲率和光照效果的光線投射方法,分別記為L(zhǎng)RT和BRT,為了便于比較,還實(shí)現(xiàn)了基于B樣條濾波的非按需計(jì)算導(dǎo)數(shù)(曲率和光照效果按需計(jì)算)的光線投射方法,記為BRTN。本文在裝有NVIDIA NVS 4200M 顯卡的intel i5機(jī)器(便攜式電腦)上測(cè)試算法。以下測(cè)試中,除非有特殊聲明,采樣步長(zhǎng)為最小網(wǎng)格間隔的1/4光線積分計(jì)算中,采用光線早中止技術(shù),不透明度決定了實(shí)際的采樣步數(shù),顏色轉(zhuǎn)換函數(shù)不影響采樣步數(shù)。在測(cè)試?yán)L制速度時(shí),對(duì)一種數(shù)據(jù),選擇一個(gè)曲率函數(shù)即可。在測(cè)試過(guò)程中,本文利用基于線性濾波的光線投射方法 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ù),采樣步長(zhǎng)為最小網(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可以看出,按需計(jì)算可以大大加快繪制速度,線性濾波快于B樣條濾波。采用線性濾波可以進(jìn)行交互繪制。
LRT先用1個(gè)Tex指令取出當(dāng)前的標(biāo)量值,結(jié)合前一采樣點(diǎn)的標(biāo)量值,得到當(dāng)前積分步的不透明度,如果不透明度非0,采用18鄰點(diǎn)方法,計(jì)算梯度,Hessian矩陣和曲率,是一個(gè)實(shí)用的方法。BRT用C.Sigg方法計(jì)算標(biāo)量值,需要用8個(gè)Tex指令取出當(dāng)前采樣的(8個(gè))鄰點(diǎn)上標(biāo)量值,結(jié)合前一采樣點(diǎn)的標(biāo)量值,得到當(dāng)前積分步的不透明度,如果不透明度為非0,用C.Sigg方法計(jì)算計(jì)算導(dǎo)數(shù),需要用84個(gè)Tex指令取出(84個(gè))鄰點(diǎn)上標(biāo)量值,然后計(jì)算曲率以及光照效果,從而有效地加速繪制過(guò)程。BRT繪制質(zhì)量比基于LRT好,LRT速度比BRT快。本文的方法可以交互選擇曲率函數(shù),突出感興趣的特征,從而幫助科學(xué)家有效地分析數(shù)據(jù)。
圖7 平均曲率和總曲率的繪制效果(BRT繪制)
表1 繪制時(shí)間比較
[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ì)算機(jī)學(xué)報(bào),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] 鄭 杰,姬紅兵. 基于動(dòng)態(tài)紋理載入的大規(guī)模數(shù)據(jù)場(chǎng)體繪制[J]. 中國(guó)圖象圖形學(xué)報(bào),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ù)場(chǎng) GPU 光線投射[J]. 中國(guó)圖象圖形學(xué)報(bào),2011,16(7): 1269-1275.
[13] Martin K,Schroeder W,Lorensen B. The Visualization ToolKit(VTK)[EB/OL]. 2008. http://www.vtk.org/