【GPU編程】體繪製：預積分傳輸函數(pre-integrated transfer function)

體繪製中的傳輸函數除了pre-classification和post-classification以外，另外一種是pre-integrated transfer function。這種方法相比前面的兩種方法，在效能和繪製品質上，都有明顯的提高。下面進行此方法的介紹。

一、原理簡介

先分類和後分類都沒有解決正確重構三維資料場所需的採用頻率問題。對於離散形式的體繪製積分方程，當採樣間隔d趨於0時，才能產生完全正確的映像。這是因為，根據從採用定律，只有滿足奈奎斯特定律才能正確重構。然後由於傳輸函數的非線性特徵，會造成體繪製積分正確賦值所需的採樣頻率的急劇增加。因為沿射線方向採樣時 ，場c (S (x) ) 和 k( s ( x ) ) 的奈魁斯特頻率為兩個傳遞
函數c( S ) 和  k ( S ) 的奈魁斯特頻率 中的最大值與標量場 S ( x ) 的 奈魁斯特頻率的乘積。進行體繪製積分時，往往由於採樣頻率不能滿足該奈魁斯 特頻率．而 造成映像品質下降，而 且過高的樣頻率會增加運算量，降低繪製速度。

分類(classification)是用傳輸函數把原始體資料中的標量值s= s(x)映射為RGBA值（顏色值和阻光度係數）。pre-integrated分類的目標是把數學積分分解為兩部分：一是沿視線方向對連續標量場s(x)採樣，此時採樣所需的奈奎斯特頻率與傳遞函數無關。二是傳輸函數顏色映射q(s)和阻光度映射k(s)。

第一步是沿著某條光線進行採樣（標量值）。注意到這次採樣的奈奎斯特頻率值不受到傳輸函數的影響。為了達到pre-integrated分類的目的，採樣值定義了一個一維的基於片段的線性標量值地區。對這種片段線性標量值地區的積分可以通過為每一個線性片段產生一個尋找表來高效實現。這個尋找表的三個參數是：前面的採樣點的值、後面的那個採樣點的值和片段的長度，一所示：

圖一 光線上第i個片段的顏色和阻光度計算

第i個片段的阻光度計算公式（公式1）：

它是前採樣點的值Sf、後採樣點的值Sb、片段長度d的函數。

第i個片段的顏色值的計算公式（公式2）：

這樣，pre-integrated分類通過計算下列公式來獲得體繪製一條光線上的的顏色積分值（公式3）：

這個公式的推導過程如下：

一條光線的體繪製積分如下：

通過對光線上第i段的分析，上式中的阻光度係數可以估算為：

ai是第i段內的累計阻光度，約等於：

有時候可以進一步簡化為：

1-ai就是第i段內的累計透明度。顏色部分也可以簡化計算為：

那麼，對於一條光線的繪製積分可以估算簡化為：

這樣，back-to-front合成演算法為：

其中，是第i段的累計顏色。

經常可以寫為：，這樣：

那麼I可以改寫為：

那麼，back-to-front合成演算法：

預積分分類的體繪製積分方法計算射線 上相鄰兩個採樣點之間線段的顏色和不透明度，該顏色和不透明度值由參數Sf、Sb和d共 同決定；而其他分類方法只計算每個採樣點位置的顏色和不透明度，由該採樣點位置的標量值 S和採樣間隔 d決定 。此與一般的體繪製方法相 比， 預積分分類體繪製方法計算映像象素強度I時 ，首先對射上每個小線段進行積分 ， 得 到厚度為 d的薄板層
的顏色和不透明度值，再進行體繪製積分運算；而先分類和後分類的方法只計算採樣點處的顏色和不透明度值作體繪製積分運算，因此也稱前者為slab-by-slab的積分方式，後者為slice-by-slice積分方式。其二所示 ：

  

圖二 兩種積分方式

預積分分類的一個缺點就是每當傳輸函數發生變化時，需要重新計算大量的光線片段的顏色和阻光度的值。這對於遊戲和娛樂應用程式來說，不是一個需要擔心的問題。因此，如果需要互動修改傳輸函數，比如科學可視化領域，則此方法不適合。

二、演算法最佳化

我們可以採取一些措施來對演算法進行一些最佳化。

1、可以通過把片段的長度設定為常量，來把映射尋找表的維度從3降為2.很明顯，這適合於等間距採樣。它也適合於3D紋理的正投影，對於透視投影也是比較好的估算。到那時對於axis-aligned的2D紋理來說，不是一個好的方法。

2、對於傳輸函數的局部改變，沒必要重新計算整個尋找表。比如傳輸函數中標量值s發生了改變，只需要重新計算尋找表中的這一一些值：C(Sf,Sb,d)，a(Sf,Sb,d)，其中Sf<=S<=Sb，或者Sb<=S<=Sf。最壞情況下，尋找表的一半內容需要重新計算。

3、引進積分函數，能夠極大地提升效率。

積分函數可以加速上面的幾個公式的計算。

例如，函數f(x)：

的積分運算可以通過引入積分函數G(s)來計算：

對於我們的預積分方法，我們只需要為所有的標量值引入一個積分函數。對於標量值Sf和Sf的所有組合的積分可以通過積分函數的計算來獲得。

這樣，阻光度的計算可以變形為：

採用得積分函數是：

要對顏色值的計算採用積分函數需要忽略光線在slab內部的衰減：

積分函數是：

三、程式實現

產生預積分尋找表的C/C++代碼：

//compute a 256 * 256 pre-integration map tableGLuint algo_preIntegrate(unsigned char * transferFuncTable){GLuint colorTexture;double r=0.f,g=0.f,b=0.f,a=0.f;int rcol,gcol,bcol,acol;double rInt[256],gInt[256],bInt[256],aInt[256];GLubyte lookupImg[256*256*4];int smin,smax;double factor,tauc;int lookupindex = 0;// compute integral functionsfor (int i1 = 1; i1 < 256; i1++){tauc = ( transferFuncTable[ ( i1 - 1 ) * 4 + 3 ] + transferFuncTable[ i1 * 4 + 3 ] ) / 2.0;r = r + ( transferFuncTable[ ( i1 - 1 ) * 4 + 0 ] + transferFuncTable[ i1 * 4 + 0 ] ) / 2.0 * tauc / 255.0;g = g + ( transferFuncTable[ ( i1 - 1 ) * 4 + 1 ] + transferFuncTable[ i1 * 4 + 1 ] ) / 2.0 * tauc / 255.0;b = b + ( transferFuncTable[ ( i1 - 1 ) * 4 + 2 ] + transferFuncTable[ i1 * 4 + 2 ] ) / 2.0 * tauc / 255.0;a = a + tauc;rInt[i1] = r; gInt[i1] = g;bInt[i1] = b;aInt[i1] = a;}rInt[0] = gInt[0] = bInt[0] = aInt[0] = 0;// compute look-up table from integral functionsfor ( int sb = 0; sb < 256; sb++ ){for ( int sf = 0; sf < 256; sf++ ){if ( sb < sf ){smin = sb;smax = sf;}else{smin = sf;smax = sb;}if ( smin != smax){factor = 1.0 / (double)(smax -smin);rcol = ( rInt[smax] - rInt[smin] ) * factor;gcol = ( gInt[smax] - gInt[smin] ) * factor;bcol = ( bInt[smax] - bInt[smin] ) * factor;acol = 256 * (1.0 - exp(-(aInt[smax] - aInt[smin])*factor/255.0));} else{if ( sb==0 && sf==0){rcol = 0;gcol = 0;bcol = 0;acol = 256 * (1.0 - exp(-(aInt[smin])*factor/255.0));}else{rcol = (lookupImg[(sb-1)*256*4+sf*4+0] + lookupImg[sb*256*4+(sf-1)*4+0]) * 0.5;gcol = (lookupImg[(sb-1)*256*4+sf*4+1] + lookupImg[sb*256*4+(sf-1)*4+1]) * 0.5;bcol = (lookupImg[(sb-1)*256*4+sf*4+2] + lookupImg[sb*256*4+(sf-1)*4+2]) * 0.5;acol = (lookupImg[(sb-1)*256*4+sf*4+3] + lookupImg[sb*256*4+(sf-1)*4+3]) * 0.5;}}CLAMP(0,rcol,255);CLAMP(0,gcol,255);CLAMP(0,bcol,255);CLAMP(0,acol,255);lookupImg[lookupindex++] = rcol;lookupImg[lookupindex++] = gcol;lookupImg[lookupindex++] = bcol;lookupImg[lookupindex++] = acol;}}//create textureglGenTextures(1, &colorTexture);glBindTexture(GL_TEXTURE_2D,colorTexture);glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP);glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP);glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);glTexImage2D(GL_TEXTURE_2D,0,GL_RGBA,256,256,0,GL_RGBA,GL_UNSIGNED_BYTE,lookupImg);return colorTexture;}

Cg程式：

頂點程式：

// Pre-integration vertex program in Cgstruct appin {float4 Position : POSITION;float4 TCoords0 : TEXCOORD0;};struct v2f {float4 HPosition : POSITION;float4 TCoords0 : TEXCOORD0;float4 TCoords1 : TEXCOORD1;};v2f main(appin IN,uniform float4x4 ModelViewProj,uniform float4x4 ModelView,uniform float4x4 ModelViewI,uniform float4x4 TexMatrix,uniform float SliceDistance){v2f OUT;// compute texture coordinate for sFOUT.TCoords0 = mul(TexMatrix, IN.TCoords0);// transform view pos and view dir to obj spacefloat4 vPosition = float4(0,0,0,1);vPosition = mul(ModelViewI, vPosition);float4 vDir = float4(0.f,0.f,-1.f,1.f);vDir = normalize(mul(ModelViewI, vDir));// compute position of sBfloat4 eyeToVert = normalize(IN.Position - vPosition);float4 sB = IN.Position- eyeToVert * (SliceDistance / dot(vDir, eyeToVert));// compute texture coordinate for sBOUT.TCoords1 = mul(TexMatrix, sB);// transform vertex position into homogeneous clip spaceOUT.HPosition = mul(ModelViewProj, IN.Position);return OUT;}

片段程式：

// Pre-integration fragment program in Cgstruct v2f {float4 TexCoord0 : TEXCOORD0;float4 TexCoord1 : TEXCOORD1;};float4 main(v2f IN,uniform sampler3D Volume,uniform sampler2D PreIntegrationTable) : COLOR{float4 lookup;// sample front scalarlookup.x = tex3D(Volume, IN.TexCoord0.xyz).x;// sample back scalarlookup.y = tex3D(Volume, IN.TexCoord1.xyz).x;// lookup and return pre-integrated valuereturn tex2D(PreIntegrationTable, lookup.xy);}

最後給出一個程式分別實現Pre-Integrated演算法和Post-Classification演算法的效果對比圖，是一個256*256*225的.raw體資料三維重建出來的結果，二者採用同樣的傳輸函數和相等的採樣步長。可以很清楚的看出前者效果要好得多~

Pre-Integrated演算法：

Post-Classification演算法：