劉宏,閔曙輝

(中國傳媒大學(xué)信息工程學(xué)院,北京 100024)

1 引言

科學(xué)技術(shù)﹑科學(xué)儀器設(shè)備的發(fā)展,產(chǎn)生了大量的序列圖像,如醫(yī)用的 CT圖像﹑ MRI圖像等。在序列圖像的應(yīng)用初期,人們是通過閱讀單個的圖像,在大腦中形成所關(guān)心部分的空間位置、形狀、大小信息的,這一過程需要有一定的經(jīng)驗(yàn),否則人腦中形成的圖形和實(shí)際情況相差很大,這樣就不能完整體現(xiàn)序列圖像的價值,比如在外科手術(shù)、牙齒修復(fù)等方面的應(yīng)用就受到制約。

人們應(yīng)用計算機(jī)﹑數(shù)學(xué)﹑圖像處理等科學(xué)技術(shù),把所關(guān)心的等值面(灰度值相等的面)繪制出來,這樣就可以形成一個三維空間的曲面,并且可以交互,極大地方便了應(yīng)用。這方面最典型的繪制曲面的方法有:等高線﹑ MC算法﹑ MT算法﹑剖分立方體方法(Dividing Cubes)。在以上方法中,人們是以三角面片代替等值面,又由于是采用線性插值的方法,所以造成的誤差較大。在形成三角面片時,MC算法有連接二義性問題﹑拓?fù)洳灰恢滦詥栴}。

分析造成上述問題的根本原因就是在做等值面時采用線性插值的方法。體素是 MC算法中的一個重要概念.它指的是一張 CT片上某方格上的四個像素點(diǎn),與另外一張 CT片上相對應(yīng)方格上的四個像素點(diǎn)所組成的長方體。在做等值面時,首先給出等值面的灰度值,然后判斷等值面與體素的每一條邊是否相交,如果相交,用線性插值的方法求出交點(diǎn),最后把這些交點(diǎn)按照一定的規(guī)則連接成三角面片,最終形成等值面。有些算法如 SMC,如果等值面與體素的某條邊相交,甚至采用了取線段最大灰度值的端點(diǎn)為等值面所經(jīng)過的點(diǎn)或取線段的最小灰度值的端點(diǎn)為等值面所經(jīng)過的點(diǎn),這樣的方法誤差會更大。本文采用在工程技術(shù)﹑計算機(jī)技術(shù)中采用的有理函數(shù)插值方法,得到了滿意的結(jié)果。

2 數(shù)學(xué)基礎(chǔ)

在數(shù)值計算里,常用的插值方法有線性插值方法﹑拉格朗日(lagrange)插值﹑哈米特(hermite)插值等。就實(shí)用的情況看,采用線性插值在三維重建中會造成偽曲面,如 CT片中的骨頭的灰度值為250,而空氣的灰度值為 0。給一個灰度值介于 0到250的灰度值 C,它就會在空氣、骨頭之間形成一個等值面。采用拉格朗日插值,由于多項(xiàng)式的階數(shù)較高,插值函數(shù) f(x)在插值點(diǎn)之間會產(chǎn)生偽擾動,這樣拉格朗日插值也不能采用。哈米特插值方法除了插值多項(xiàng)式的階數(shù)較高外,還必須知道插值點(diǎn)的導(dǎo)數(shù)值,這在實(shí)際應(yīng)用中難以辦到。

有理函數(shù)插值是在工程技術(shù)﹑計算機(jī)算法中采用的方法,它的數(shù)學(xué)模型為:

當(dāng) x=xi時, yi=f(xi) i=0,1,2,…,n+m.要找一個有理函數(shù)=f(xi), i=0,1,2,…n+m

在構(gòu)造有理函數(shù)插值時,用到連分式的概念以及一些特殊記法,下面給出例子加以說明。給出有理分式用輾轉(zhuǎn)相除化為連分式。

由上面把有理函數(shù)寫成連分式的過程及記號方法,可以把滿足 R(xi)=f(xi) i=0,1,2,…,n的有理插值函數(shù)記為:

其中 ck(k=0,1,…,n)由插值條件 R(xi)=f(xi) i=0,1,…,n確定。

下面用反差商的方法得到 ck(k=1,2,…,n+1),同時也得到了 R(x)。

定義 1 給定一組點(diǎn)集{xi, i=0,1,2,…},如果函數(shù)序列{vk(x)}滿足如下關(guān)系

就稱 vk(x)為函數(shù) f(x)在點(diǎn)集{xi, i=1,2,…}上的 k階反差商。

(舍去最后一項(xiàng)得)

由上可以看出,若是求出了 vi(xi),i=0,1,2,……,也就求出了 R(x)。即 f(x)的有理插值函數(shù)。求 vi(xi) i=0,1,2,……,n可以采用構(gòu)造反差商表來得到。

表1

3 有理插值的程序?qū)崿F(xiàn)

熟悉了有理插值的數(shù)學(xué)原理后,可以看出用遞歸方法編程實(shí)現(xiàn)該算法是比較可行的辦法。

首先,我們需要一個函數(shù)來計算反差商表中各個值。從表 1的規(guī)律總結(jié)出以下遞歸函數(shù):

double value(double x[],double y[],int k,int i){

if(k==0)return y[i];else

return(x[i]-x[k-1])/(value(x,y,k-1,i)-value(x,y,k-1,k-1));

}

下面說明一下函數(shù)的幾個參數(shù):x[]和 y[]這 2個數(shù)組分別代表進(jìn)行有理插值的一組抽樣點(diǎn)的 x坐標(biāo)和 y坐標(biāo)數(shù)列。k和 i與 vk(xi)=

這個函數(shù)用來計算有理插值算法產(chǎn)生的連分式,同樣使用了遞歸方法,并調(diào)用了第一個函數(shù)取得反差商表中的數(shù)據(jù)。與 value函數(shù)一樣,參數(shù)中的 x[]和 y[]這 2個數(shù)組分別代表進(jìn)行有理插值的一組抽樣點(diǎn)的x坐標(biāo)和y坐標(biāo)數(shù)列。n是指有理插值所用的抽樣點(diǎn)個數(shù)。x0表示我們要由長度為 n的一組采樣點(diǎn) x[]和y[]通過有理插值得到橫坐標(biāo)為 x0處的估計縱坐標(biāo)。這時,有理插值得到的縱坐標(biāo)用 y0=value(x,y,0,0)+rational(x0,x,y,1,n)來表示。

有了有理插值的程序?qū)崿F(xiàn)方法,下面就通過一個簡單的函數(shù)來簡單測試該算法的性能:

對于對數(shù)函數(shù) ln(x+1),在區(qū)間[0,1]內(nèi)等距離插入四個點(diǎn),分別求出 ln(x+1)、線性插值、有理函數(shù)插值在各個插值點(diǎn)中間的值,可以看出有理函數(shù)插值的精度是線性插值精度的一萬多倍。比如 ln(x+1)的精確值為 0.0953098,線性插值的函數(shù)值為0.0911608,有理函數(shù)插值的函數(shù)值為0.0953098.線性插值的絕對誤差是 0.0041494,有理函數(shù)插值的絕對誤差是 3.53e-0.07。有理函數(shù)插值的精度是線性插值精度的 11753.6倍。

在同一個坐標(biāo)中畫出 ln(x+1)的圖像以及 ln(x+1)的有理插值﹑線性插值函數(shù)的圖像,可以看出,有理插值函數(shù)的圖像是和 ln(x+1)的圖像是重合的。如圖 1所示。中的 k和 i這 2個下標(biāo)是對應(yīng)的。從反差分表中可以看到,當(dāng) k為 0時,也就是第一列數(shù)值,就等于 y[]數(shù)組相應(yīng)的值。

另外還需要一個函數(shù)來計算有理插值計算出來的值。

圖 1 ln(x+1)和它的線性插值函數(shù)、有理函數(shù)插值函數(shù)的圖像

4 高精度的曲面繪制方法

選取一系列 CT斷層圖像,輸入這些圖像數(shù)據(jù),由于原始數(shù)據(jù)可能有噪聲,所以先將原始數(shù)據(jù)經(jīng)過2次中值濾波,并把這些處理后的數(shù)據(jù)存放在一個數(shù)組里。CT斷層圖像的灰度值是標(biāo)量數(shù)據(jù),數(shù)據(jù)類型是浮點(diǎn)型的。

以往的曲面繪制方法是在給出一個灰度值后判斷等值面與每一個體素的每一條邊是否相交,如果相交,用線性插值的方法找到等值面與體素的交點(diǎn),然后用三角面片來近似代替等值面。高精度的曲面繪制方法是:給出一個等值面的灰度值 c0,先在每一張 CT片中找到灰度值為 c0的若干個點(diǎn),用點(diǎn)列{(xi,yj,zij)}i=1,2,…,n;j=1,2,…,m表示,其中{(xi,yj,zi1)}ni=1表示從第一張 CT片中取出的灰度值為 c0的點(diǎn)列,其它的可類推。

對于每一個點(diǎn)列{(xi,yj,zij)}ni=1,過這個點(diǎn)列,用有理插值函數(shù) R(x)來逼近這條等高線,用 Jj(j=1,2,…,m)表示。同樣地,對于點(diǎn)列 {(xi,yj,zij)}mi=1,過這個點(diǎn)列,用有理插值函數(shù) R(x)來逼近這一等高線,用 wj(j=1,2,…,n)表示。在三維空間中,把這兩組等高線畫出來,就可以得到等值面的近似曲面。

下面兩個圖像,分別是用高精度面繪制方法得到的皮膚表面和用 MC算法得到的等值面。從圖像的效果看,高精度面繪制的曲面網(wǎng)能夠很好地貼合皮膚,視覺效果也十分滿意。如圖 2和圖 3所示。

圖2 高精度面繪制的頭部皮膚的效果圖

圖3 MC算法繪制的頭部皮膚的效果圖

需要特別強(qiáng)調(diào)的是在尋找灰度值為 c0的點(diǎn)(x0,y0,z0)時,程序判斷[f(x0,y,z0)-c0]×[f(x0,y+1,z0)-c0]≦ 0是否成立,其中 f(x0,y,z0)是圖像在點(diǎn)(x0,y,z0)的灰度值,如果成立,說明灰度值為c0的點(diǎn)就落在(x0,y,z0)與(x0,y+1,z0)之間。由于這兩點(diǎn)之間的距離很小,所以用線性插值的方法尋找 y0。令y+k。

5 總結(jié)

有理函數(shù)插值算法是一種非常精確的曲線擬合算法,尤其當(dāng)實(shí)際曲線比較接近一些理想數(shù)學(xué)函數(shù)時,通過有理函數(shù)插值可以由較少的抽樣點(diǎn)擬合出精確度很高的曲線。將其用于醫(yī)學(xué)圖像三維重建不會產(chǎn)生傳統(tǒng) MC算法容易引起的連接二義性問題,同時又能繪制出精確度較高的表面,其誤差也遠(yuǎn)遠(yuǎn)小于傳統(tǒng)使用的線性插值算法,而且擬合出的曲線更加美觀﹑光滑、自然。根據(jù)有理函數(shù)插值算法的以上特點(diǎn),該方法不僅僅可以用于醫(yī)學(xué)圖像,而且也可以在其他研究方面得到廣泛的應(yīng)用,如表面設(shè)計,工業(yè)器件三維圖像重建等。

[1] Burden R L,Faires J D.數(shù)值分析(第七版影印版)[M].北京:高等教育出版社,2001.

[2] 朱恒軍,楊莘元.邊緣檢測技術(shù)在 CT圖像預(yù)處理中的應(yīng)用[J].哈爾濱大學(xué)學(xué)報,2007,23(1).

[3] 陳秦玉 .人體三維重建的實(shí)踐和技術(shù)研究[C].浙江大學(xué)博士學(xué)位論文,2004.

[4] 羅朝東.面向有封閉內(nèi)腔和外形的復(fù)雜工件快速原型制造的 ICT圖像處理及反求工程研究[C].重慶大學(xué)碩士學(xué)位論文,2006.

[5] 安洪振,劉金匯.用于工業(yè) CT圖像三維重建的邊緣信息提取算法研究[J].核電子學(xué)與探測技術(shù),2007,27(5).

[6] 李慶揚(yáng),關(guān)治,白峰杉.數(shù)值計算原理[M].北京:清華大學(xué)出版社,2000.