序列圖像的高精度面繪制方法

劉宏,閔曙輝

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

1 引言

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

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

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

2 數學基礎

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

有理函數插值是在工程技術﹑計算機算法中采用的方法,它的數學模型為:

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

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

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

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

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

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

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

(舍去最后一項得)

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

表1

3 有理插值的程序實現

熟悉了有理插值的數學原理后,可以看出用遞歸方法編程實現該算法是比較可行的辦法。

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

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));

}

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

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

有了有理插值的程序實現方法,下面就通過一個簡單的函數來簡單測試該算法的性能:

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

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

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

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

4 高精度的曲面繪制方法

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

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

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

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

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

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

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

5 總結

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

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

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

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

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

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

[6] 李慶揚,關治,白峰杉.數值計算原理[M].北京:清華大學出版社,2000.