李長春, 李元金

(1.滁州職業(yè)技術(shù)學(xué)院 信息工程學(xué)院，安徽 滁州 239000; 2.滁州學(xué)院 計(jì)算機(jī)與信息工程學(xué)院， 安徽 滁州 239000)

近年來，隨著科技的發(fā)展，CT[1]成像系統(tǒng)越來越廣泛地應(yīng)用于臨床醫(yī)學(xué)領(lǐng)域.正常情況下，CT成像系統(tǒng)能夠提供清晰的CT圖像.但是，當(dāng)投射區(qū)域含有金屬等高密度物質(zhì)時，重建之后的CT圖像將出現(xiàn)從金屬等高密度物質(zhì)出發(fā)向四周發(fā)射的條狀偽影，這就是金屬偽影.這些條狀偽影不僅降低了CT圖像的對比度，而且使得金屬等高密度物質(zhì)與人體組織之間變得模糊.有時，這種偽影會阻礙醫(yī)務(wù)工作者對病理的診斷，甚至阻礙了CT成像系統(tǒng)的應(yīng)用.

為了提高CT圖像的清晰度和對比度，以便提高對疾病診斷的正確率，人們提出了多種CT圖像金屬偽影校正方法.總體上來說，從已發(fā)表的論文可以看出，CT圖像金屬偽影校正方法可以分為基于硬件的CT圖像金屬偽影校正方法(hardware-based correction approach)和基于軟件的CT圖像金屬偽影校正方法(software-based correction method).基于硬件的CT圖像金屬偽影校正方法，是CT圖像出現(xiàn)金屬偽影后最早的校正方法.在這些方法中，可以將投影前的金屬等高密度物質(zhì)從人體中取走，或者增加投影時的電壓，但是基于硬件的校正方法不易實(shí)現(xiàn)，并且對人體有害，所以人們更多地求助于基于軟件的CT圖像金屬影響校正方法.基于軟件的CT圖像金屬偽影校正方法(software-based correction method)又可以分為基于插值的CT圖像金屬偽影校正方法(interpolation-based correction method)[2-10]和基于迭代的CT圖像金屬偽影校正方法[11-13](iteration-based correction method).受限于算法運(yùn)行速度和現(xiàn)有CT成像系統(tǒng)體系結(jié)構(gòu)等方面的限制，雖然有一些學(xué)者研究了基于迭代的CT圖像金屬偽影校正方法(iteration-based correction method)，但實(shí)踐中應(yīng)用甚少.相反，人們更多地求助于基于插值的CT圖像金屬偽影校正方法(interpolation-based correction method).

1987年，Kalender Willi A在文獻(xiàn)[2]中首次提出了使用簡單的線性插值方法校正CT圖像中出現(xiàn)的金屬偽影，這是最早的基于插值的CT圖像金屬偽影校正方法(interpolation-based correction method).在這之后，學(xué)者們在線性插值CT圖像金屬偽影校正方法的基礎(chǔ)上對多項(xiàng)插值在CT圖像金屬偽影校正方法進(jìn)行了研究.本文將在前人結(jié)果的基礎(chǔ)上，研究多項(xiàng)式次數(shù)對CT圖像金屬偽影校正方法性能影響比較研究.

1 多項(xiàng)式插值校正方法

圖1(A)、(B)分別顯示了基于插值的CT圖像金屬偽影校正算法流程圖.其中，圖1(A)的CT圖像金屬偽影校正算法主要步驟如下：

1)投影數(shù)據(jù)輸入；

2)CT圖像重建；

3)CT圖像分割和投影；

4)插值；

5)CT圖像重建和金屬對象補(bǔ)償.

圖1(B)的CT圖像金屬偽影校正算法主要步驟如下：

1)CT圖像重建；

2)CT圖像分割和投影；

3)插值；

4)CT圖像重建和金屬對象補(bǔ)償.

圖1(A)與(B)主要區(qū)別在于圖1(A)起源于投影數(shù)據(jù)，而圖1(B)起源于CT圖像.

1.1 投影數(shù)據(jù)輸入

投影數(shù)據(jù)來源于CT成像系統(tǒng)，是CT成像系統(tǒng)通過掃描投射體后得到的結(jié)果.在基于插值的CT圖像金屬偽影校正算法流程圖中，圖1(A)通過輸入獲取到的投影數(shù)據(jù)，得到含有金屬偽影的CT圖像，當(dāng)CT成像系統(tǒng)中投影數(shù)據(jù)獲取不到或者很難獲取到時，直接利用含有金屬偽影CT圖像，如圖1(B).

圖1 基于插值的CT圖像金屬偽影校正算法流程圖Figure 1 Metal artifact correction algorithm flow chart for CT image based on interpolation

1.2 CT圖像重建

CT圖像重建是由投影數(shù)據(jù)重建成CT圖像的過程，為醫(yī)學(xué)診斷提供了手段.在這個過程中，使用的方法主要包括解析法和迭代重建算法.其中迭代法又包括代數(shù)重建法(Algebraic reconstruction technique，ART)、聯(lián)合迭代重建法(Simultaneous Iterative Reconstruction Technique，SIRT)和基于統(tǒng)計(jì)學(xué)的優(yōu)化方法(Optimization Method based on Statistics，OMS)；解析法直接傅立葉重建法和波反抽影法.

迭代法重建中首先假設(shè)斷層截面是由一個未知的數(shù)字矩陣組成的，然后由測量投影數(shù)據(jù)建立一組未知向量的代數(shù)方程式，通過方程組求解圖像向量.迭代重建算法由于計(jì)算代價大、普適性較差，僅在少數(shù)場合應(yīng)用.

相反，由于解析法具有重建速度快而且符合當(dāng)代CT成像系統(tǒng)體系結(jié)構(gòu)，因此，實(shí)際中使用較多.在解析法中，濾波反投影重建算法(Filter Back Projection，F(xiàn)BP)是目前常用的CT圖像重建算法，在實(shí)際之中使用最多.

1.3 圖像分割和重投影

圖像分割是圖像處理領(lǐng)域常用的名詞之一.它的意思是將預(yù)處理后的圖像根據(jù)圖像處理目的將圖像劃分成若干個特定的、具有獨(dú)特性質(zhì)的區(qū)域.圖像分割是圖像處理的關(guān)鍵步驟和圖像處理后續(xù)步驟的前提和基礎(chǔ)，也是計(jì)算機(jī)視覺研究中的一個經(jīng)典難題，已經(jīng)成為圖像理解領(lǐng)域關(guān)注的一個熱點(diǎn).由于不同圖像具有自己的特點(diǎn)，因此，時至今日沒有統(tǒng)一的圖像分割方法. 一般來說，基于閾值法的方法是最早的圖像分割方法.根據(jù)分割時獲取閾值方法的不同，基于閾值法的圖像分割方法又可以分為全局閾值法、局部閾值法、自動閾值法等等.除了基于閾值的圖像分割方法外，圖像分割方法還包括基于區(qū)域的圖像分割、基于邊緣信息的圖像分割方法、基于形態(tài)學(xué)的圖像分割方法、基于先驗(yàn)信息的圖像分割方法、基于特定理論的圖像分割方法等等.

為了實(shí)驗(yàn)簡單起見，本文實(shí)驗(yàn)時使用了基于閾值的圖像分割方法對采集來的CT圖像進(jìn)行分割.閾值分割是最常見、最簡單，也最適用的圖像分割算法.該算法主要根據(jù)圖像像素的灰度值來確定對應(yīng)的分割閾值.根據(jù)選閾值方式的不同，可以分為以下兩種形式.1)如果圖像目標(biāo)單一，只需要選取一個確定的閾值，就可以將圖像劃分為目標(biāo)和背景兩大類，這就是單閾值分割；2)如果圖像目標(biāo)比較復(fù)雜，需要根據(jù)目標(biāo)區(qū)域和背景區(qū)域的特點(diǎn)選擇不同的閾值對圖像進(jìn)行分割，這就是多閾值圖像分割.在使用閾值對圖像分割過程中，最關(guān)鍵是如何取得一個合適的閾值.式(1)是使用閾值對圖像進(jìn)行分割公式：

(1)

其中：T是閾值，f(x,y)是待分割的圖像在坐標(biāo)(x,y)處的灰度值，g(x,y)是分割后的圖像灰度值.公式(1)的表示圖像在(x,y)處的灰度值大于T時，分割之后圖像的灰度值為a，表示圖像在(x,y)處的灰度值小于等于T時，分割之后圖像的灰度值為b.通常情況下，a和b分別為1和0，分別表示前景與背景.從公式(1)中可以看出，T的確定對f(x,y)分割起著關(guān)鍵性的作用，只有確定精確的T才能精確的把圖像分割開；相反，如果T確定不好，將很難將目標(biāo)圖像和背景圖像分割開.

1.4 插值

通過確定閾值并對圖像進(jìn)行分割后，原圖像被分割成金屬圖像和非金屬圖像.為了消除原CT圖像中的金屬偽影，算法將獲取到的金屬圖像與非金屬圖像進(jìn)行投影，并得到對應(yīng)的金屬圖像投影正弦圖，非金屬圖像投影正弦圖.

為了消除CT圖像中的金屬偽影，首先，根據(jù)金屬對象投影正弦圖確定其在原始CT圖像投影正弦圖中的位置與形狀，然后使用非金屬投影正弦圖數(shù)據(jù)對金屬對象投影正弦圖中的數(shù)據(jù)進(jìn)行插值并填充對應(yīng)的位置.

實(shí)驗(yàn)過程中，為了加快算法的運(yùn)行速度，這里使用了插值方式去除金屬偽影.插值法有內(nèi)插值法和外插值法，這里主要用了內(nèi)插值法.內(nèi)插值法是一種根據(jù)已經(jīng)二個或以上端點(diǎn)坐標(biāo)值，求出端點(diǎn)之間某點(diǎn)坐標(biāo)及其對應(yīng)值.

根據(jù)插值時所使用基的不同，又可以分為拉格朗日插值、牛頓插值、埃特肯插值、埃爾米特插值和樣條插值等.這里重點(diǎn)介紹一下拉格朗日插值，以滿足實(shí)驗(yàn)的需要.

拉格朗日插值問題可表述為：求作次數(shù)≤n多項(xiàng)式Ln(x)使?jié)M足插值條件：

Ln(x)=yi(i=0,…,n)即為拉格朗日(Lagrange)插值.

拉格朗日(Lagrange)插值公式的基本思想是：把構(gòu)造插值多項(xiàng)式的問題轉(zhuǎn)化為構(gòu)造(n+1)個插值基函數(shù)li(x)(i=0,1,…,n).

為了比較次數(shù)對CT圖像校正結(jié)果的影響，這里選擇了拉格朗日線性插值、拉格朗日二次多項(xiàng)式插值和拉格朗日四次多項(xiàng)式插值.

1.4.1 拉格朗日線性插值[14-15]

拉格朗日線性插值的問題：已知函數(shù)y=f(x)在點(diǎn)x0，x1上的值為y0，y1，求作一次式，使?jié)M足條件：

L1(x0)=y0,L1(x1)=y1

(2)

如圖2所示可知.由直線兩點(diǎn)式可知，通過A，B的直線方程為，

圖2 線性插值示意圖Figure 2 Linear interpolation diagram

(3)

式(3)稱為線性插值公式.該式又可以被寫成：

L1(x)=y0l0(x)+y1l1(x)

(4)

(5)

1.4.2 二次多項(xiàng)式插值

拉格朗日(Lagrange)二次多項(xiàng)式插值問題可以表述為：求作二次多項(xiàng)式L2(x)，使?jié)M足條件：

L2(xj)=yj(j=0,1,2)

(6)

拉格朗日二次多項(xiàng)式插值是拉格朗日線性插值的延伸和擴(kuò)展，在幾何解釋是用通過三個點(diǎn)的拋物線來近似考察曲線.類似于拉格朗日線性插值，構(gòu)造基函數(shù)，要求滿足：

L2(x)=y0l0(x)+y1l1(x)+y2l2(x)

(7)

即拉格朗日二次多項(xiàng)式插值表達(dá)式為：

(8)

1.4.3 四次多項(xiàng)式插值

拉格朗日(Lagrange)四次多項(xiàng)式插值問題可以表述為：求作四次多項(xiàng)式L4(x)，使?jié)M足條件：

L4(xj)=yj(j=0,…,4)

(9)

同樣地，拉格朗日四次多項(xiàng)式插值是拉格朗日線性插值的延伸和擴(kuò)展，是拉格朗日多項(xiàng)式插值的特例.類似于拉格朗日線性插值，構(gòu)造基函數(shù)，要求滿足：

L4(x)=y0l0(x)+y1l1(x)+

y2l2(x)+y3l3(x)+y4l4(x)

(10)

即拉格朗日四次多項(xiàng)式插值表達(dá)式為：

(11)

1.5 CT圖像重建和金屬對象填補(bǔ)

在非金屬投影正弦圖中，金屬對象所在區(qū)域使用插值方法獲取數(shù)據(jù)并填充后(不含金屬部分)，使用濾波反投影算法重建校正后的正先弦圖，得校正后的圖像.由于重建后的CT圖像中不含有對應(yīng)的金屬對象補(bǔ)充，必須將金屬對象填補(bǔ)到重建后的CT圖像中，得到最終CT圖像.

2 實(shí)驗(yàn)與結(jié)果分析

為了比較多項(xiàng)式次數(shù)對CT圖像校正結(jié)果的影響，使用了臨床CT圖像.臨床CT圖像是從國外某醫(yī)療機(jī)構(gòu)網(wǎng)站上下載.實(shí)驗(yàn)環(huán)境是Matlab R2012a和Win7操作系統(tǒng).

實(shí)驗(yàn)時，參加比較的量主要是均方根差(root mean square error，RMSE)和 以及系統(tǒng)運(yùn)行時間.其中均方根差，也稱為方均根偏移是一種常用的測量數(shù)值之間差異的量度，計(jì)算公式如(12)所示：

(12)

其中：M和N分別表示圖像的寬度和高度.xij是圖像中第i行和第j列的灰度值.X是圖像灰度值的平均值.

圖3顯示了理想圖像、模擬圖像和校正后圖像，表1是其對應(yīng)參數(shù)表.

表1 模擬模型參數(shù)Table 1 simulation model parameters

表2顯示了模型理想圖像與含有金屬偽影的原始圖像、線性插值校正后圖像、二次多項(xiàng)式插值校正后圖像和四次多項(xiàng)式插值校正圖像.其中“L-MAR圖像”、“Q1-MAR圖像”和“Q2-MAR圖像”分別表示使用線性插值校正結(jié)果、二次多項(xiàng)式插值校正結(jié)果和四次多項(xiàng)式校正結(jié)果.

表2 各種圖像與參考圖像之間的均方根差Table 2 Root mean square difference between various images and Reference image

從表2可以看出，原始圖像的均方根差是0.4270，L-MAR圖像的均方根差是0.3465，Q1-MAR圖像的均方根差是0.3209，Q2-MAR圖像的均方根差是0.3102.

圖4表示線性插值校正法運(yùn)行時間、二次多項(xiàng)式插值校正法運(yùn)行時間和四次多項(xiàng)式插值校正法運(yùn)行時間.

從圖4可以看出，線性插值校正算法運(yùn)行時間為0.297 704，二次多項(xiàng)式插值校正算法運(yùn)行時間為0.302 814，四次多項(xiàng)式插值校正算法運(yùn)行時間為0.325 467.

圖4 線性插值、二次多項(xiàng)式插值和四次多項(xiàng)式插值運(yùn)行時間Figure 4 Running time of linear interpolation,quadratic polynomial interpolation and quartic polynomial interpolation

3 結(jié) 語

放射狀的金屬偽影是影響CT圖像質(zhì)量的重要因素之一.因此，在基于CT圖像的臨床應(yīng)用領(lǐng)域，在使用CT圖像之前必須要進(jìn)行CT圖像清晰化.線性插值作為最早使用的校正方法具有簡單且運(yùn)行快的特點(diǎn).在線性插值的基礎(chǔ)上，人們又提出的多項(xiàng)式插值.通過實(shí)驗(yàn)，本文比較了多項(xiàng)式的次數(shù)對CT圖像校正結(jié)果的影響.實(shí)驗(yàn)結(jié)果表明，一次線性插值校正方法運(yùn)行速度最快、二次線性插值校正方法運(yùn)行速度次之，四次線性插值校正方法運(yùn)行速度最慢.均方根差(RMSE)方面，原始含有金屬偽影的CT圖像最大，一次線性插值校正方法次之、二次線性插值校正方法再次之，四次線性插值校正方法最小.

總而言之，在使用多項(xiàng)式插值校正CT圖像金屬偽影的過程中，隨著多項(xiàng)式次數(shù)的增加程序運(yùn)行的越來越慢，運(yùn)行時間越來越長；圖像質(zhì)量方面，均方根差(RMSE)越來越小.