混合高斯/泊松最大似然函數(shù)下的CBCT圖像重建

鄭蓉珍，趙 芳，李 波，田 昕*

(1.武漢大學(xué) 電子信息學(xué)院，湖北 武漢 430072;2.武漢大學(xué) 中南醫(yī)院心血管內(nèi)科，湖北 武漢 430071;3.武漢大學(xué) 口腔醫(yī)院放射科，湖北 武漢 430079)

1 引 言

錐形束計算機(jī)斷層成像技術(shù)(Cone Beam Computed Tomography，CBCT)因其X射線輻射劑量低、分辨率高、數(shù)據(jù)采集時間短等優(yōu)點，近年來逐漸成為CT領(lǐng)域的研究熱點。在低輻射劑量下，CBCT圖像重建質(zhì)量容易受到噪聲等因素的干擾。加大輻射劑量可以提升圖像質(zhì)量，但同時也會增加患者患癌和遺傳病變的風(fēng)險[1]。所以，如何在低輻射劑量下實現(xiàn)高質(zhì)量CBCT圖像重建是一個非常有意義的研究問題。

CBCT三維重建算法包括解析算法和迭代算法。以Feldkamp-Davis-Kress (FDK)算法三維圖像重建算法，一種易于實現(xiàn)的基于圓軌跡掃描的近似三維圖像重建算法為代表的解析算法，重建速度快，在小錐角時具有良好的重建效果，但如果投影數(shù)據(jù)采樣不充分或是測量噪聲較大時，解析法重建出的圖像效果不甚理想，且對噪聲敏感。而各種迭代算法，例如同時型迭代重建算法(Simultaneous Algebraic Reconstruction Technique，SART)，可以恢復(fù)出高質(zhì)量的圖像，但其缺點同樣是噪聲會在迭代過程中被逐漸放大，嚴(yán)重影響重建圖像質(zhì)量。目前對CBCT圖像去噪的主要思路是首先對投影數(shù)據(jù)去噪，然后通過去噪后的投影數(shù)據(jù)實現(xiàn)CBCT圖像的三維重建。在去噪方法中，非局部均值算法具有良好的性能，但算法計算量大，去噪時間比較長[2]，因此，在CBCT去噪中應(yīng)用較為受限。2010年，Yong Yin等對小波算法提出改進(jìn)，提出了適合CBCT圖像特征的去噪算法[3]，但該方法不具有自適應(yīng)性。近年來，三維塊匹配濾波(Block-Matching and 3D filtering, BM3D)算法具有較好的去噪性能，被廣泛應(yīng)用于CBCT的圖像去噪中[4]。

綜上所述，傳統(tǒng)方法中對于低輻射劑量CBCT圖像重建中的噪聲去除是先去噪，再重建。但是，這種處理方式存在如下問題：(1)雖然可以去除投影數(shù)據(jù)中的噪聲，但是同時也降低了投影數(shù)據(jù)中圖像邊緣的質(zhì)量，從而使得由投影數(shù)據(jù)重建的三維圖像往往會出現(xiàn)過平滑的現(xiàn)象；(2)在去噪過程中，未充分考慮低輻射劑量下CBCT投影數(shù)據(jù)中真實數(shù)據(jù)和噪聲的分布特性。一方面，探測器在工作過程中會不可避免的產(chǎn)生服從高斯分布的加性系統(tǒng)噪聲；另一方面由于輻射劑量的減少，所以在X射線穿過被測人體時會產(chǎn)生更嚴(yán)重的光電效應(yīng)、光子散射等現(xiàn)象，在光電轉(zhuǎn)換過程中光子數(shù)的減少會產(chǎn)生嚴(yán)重的泊松噪聲[5]，因此，在低輻射劑量情況下CBCT系統(tǒng)投影數(shù)據(jù)噪聲分布近似為復(fù)雜的高斯/泊松混合分布模型，而不是傳統(tǒng)的高斯噪聲模型。

本文將噪聲去除與圖像重建過程建立在統(tǒng)一的重建模型中，它包含一個基于混合高斯/泊松最大似然函數(shù)的保真項和一個基于三維全變分正則化方法的約束項。保真項用于描述混合高斯/泊松噪聲環(huán)境下重建值與觀測值盡可能的相似，約束項要求在重建過程中去除噪聲的同時盡可能地保證重建圖像的邊緣信息。我們通過可分離近似法和擴(kuò)展拉格朗日法對上述模型進(jìn)行了求解，并通過仿真數(shù)據(jù)和真實數(shù)據(jù)驗證了算法的有效性。

2 圖像重建算法

2.1 CBCT三維圖像重建原理及數(shù)學(xué)描述

典型的基于等距算法原理的CBCT三維圖像重建過程描述如下：

如圖1(a)所表示的笛卡爾坐標(biāo)系xyz旋轉(zhuǎn)β角度可得圖1(b)tsz表示的坐標(biāo)系圖，其中p′和ξ分別表示經(jīng)過待重建點P(t,s,z)映射到虛擬探測器上的橫縱坐標(biāo)，D表示射線源S到到虛擬探測器坐標(biāo)原點O的距離。

圖1 FDK掃描投影結(jié)構(gòu)示意圖Fig.1 FDK scan projection structure diagram

FDK算法的重建過程可以通過式(1)進(jìn)行表示：

f(t,s,ξ)=

(1)

其中：U表示加權(quán)因子，Rβ(p,ξ)表示在角度β下經(jīng)過被測物體的投影，具體過程可以參見文獻(xiàn)[6]。

2.2 基于混合高斯泊松最大似然函數(shù)的CBCT三維圖像重建算法

2.2.1 混合高斯/泊松噪聲環(huán)境下的CBCT圖像重建模型

基于以上推導(dǎo)本文所提出的CBCT重建模型可以表示為：

(2)

其中：x代表需要重建的三維CBCT圖像；f(x)是最大似然函數(shù)的保真項，用于約束重建值與觀測值的相似程度；g(x)是三維全變分正則化項，用于去除重建過程中的噪聲；λ是正則化參數(shù)。

對于高斯噪聲而言，假設(shè)其分布服從標(biāo)準(zhǔn)正態(tài)高斯分布而言，最大化似然函數(shù)f(x)可以得到：

(3)

其中：Ai代表服從式(1)的第i個角度下的投影矩陣，yi代表第i個角度的投影數(shù)據(jù)(觀測值)，M代表總的投影數(shù)目。

(a)混合高斯/泊松最大似然函數(shù)的保真項

(4)

其中：[bi]k和[yi]k分別代表向量bi和yi中第k個數(shù)據(jù)。式(4)較為復(fù)雜，不宜直接求解其最大似然函數(shù)。采用GAT算法的類似思路[7]，對輸入投影數(shù)據(jù)進(jìn)行Anscombe變換，其過程如下：

(5)

在此情況下yi的概率密度函數(shù)可以表示為[8]：

p(yi|bi)=

(6)

已知bi=Aix，那么關(guān)于x的最大似然函數(shù)可以表示為[9]：

(7)

(b)三維全變分正則化項

g(x)是正則化項，考慮到重建過程中噪聲的影響，這里可以引入三維全變分正則化項(3D Total Variation，3DTV)[10]來對重建圖像進(jìn)行約束，使得重建過程在去除噪聲的同時能夠較好地保留圖像的邊緣和細(xì)節(jié)信息。這里采用各項同性3D全變分正則化，此時g(x)可以表示為：

g(x)=‖x‖TV=

(8)

其中：βx,βy和βz是沿著不同方向的權(quán)重系數(shù),Dx,Dy和Dz分別代表x,y和z方向的差分計算符號。假設(shè)x是三維重建圖像x(x,y,z)對應(yīng)的一維矢量(將三維圖像按照首尾鏈接的方式變成一維矢量)，該轉(zhuǎn)化過程用vec(·)表示，則Dxx，Dyx和Dzx的計算過程如下：

(9)

邊界處的值可以采用補(bǔ)零的方式進(jìn)行填充。

2.2.2 模型優(yōu)化求解

為了描述方便，在下文中變量均轉(zhuǎn)化為矢量進(jìn)行描述?；诳煞蛛x近似的方法[11]，式(2)可以轉(zhuǎn)化為如下兩個迭代步驟：

(10)

其中：

(11)

f(x)=

(12)

(13)

其中：ρr是正則化參數(shù)，u是中間變量，y是拉格朗日乘子。與D類似，u和y可以表示為：

分別對x,u和y進(jìn)行求解，從而可以將式(13)轉(zhuǎn)化為求解如式(14)所示的三個子問題：

(14)

可以得到：

這里擬采用Barzilar-Borwein方法來更新公式(11)中的迭代步長a(k)[14]。

令g(k)為f(x)在第k次迭代時的梯度值，則：

(16)

其中：s(k-1)=x(k)-x(k-1),y(k-1)=g(k)-g(k-1)。

2.2.3 算法描述

總結(jié)前文步驟，所提出算法的流程如下：

算法:混合高斯/泊松最大似然函數(shù)下的CBCT圖像重建輸入:投影數(shù)據(jù)yi,投影矩陣Ai,循環(huán)次數(shù)TO和TI,正則化參數(shù)λ,其他參數(shù)ρr,βx,βy,βz,σ2變量初始賦值:x(1),a(1),y(1),u(1) 1.根據(jù)公式(5),計算yi;2.For k=1: TO3.根據(jù)公式(12),計算f(x);4.根據(jù)公式(11),計算z(k); For t=1: TI5.根據(jù)公式(15),計算x(t);6.根據(jù)公式(16),計算u(t);7.根據(jù)公式(17),計算y(t); End8.令x(k)=x(TI);9.根據(jù)公式(18),計算a(k); End輸出:x

3 實驗結(jié)果與分析

3.1 仿真數(shù)據(jù)實驗

仿真實驗數(shù)據(jù)為head phantom模體，模體大小為128×128×128，每個體素大小為4×4×4 mm3,投影規(guī)格為256×200×211，投影的角度為0°～210°，射線源到探測器的距離為1 500 mm，射線源到模體的距離為1 100 mm。通過調(diào)整仿真數(shù)據(jù)中平均入射光子數(shù)I0仿真不同劑量的X射線下得到的不含噪聲的投影數(shù)據(jù)，對投影數(shù)據(jù)添加不同程度的混合高斯/泊松噪聲進(jìn)行后續(xù)實驗，對其添加泊松噪聲，泊松噪聲的大小與I0大小成負(fù)相關(guān)關(guān)系，添加的高斯噪聲通過調(diào)整期均方差σ調(diào)整其大小，σ越大添加的高斯噪聲越大。

(a)主觀結(jié)果

對模體在I0=50 000的輻射劑量情況下得到的投影數(shù)據(jù)添加混合高斯/泊松噪聲后各算法重建head phantom模體取第45個切片，其中泊松噪聲的大小與I0大小成負(fù)相關(guān)關(guān)系，高斯噪聲：σ=50。重建圖像中紅色方框區(qū)域為目標(biāo)放大區(qū)域，簡稱ROI區(qū)域。圖2第1行分別表示FDK，SART，預(yù)先用BM3D對投影數(shù)據(jù)去噪然后用FDK算法重建和本文算法的重建結(jié)果，圖2第2行分別表示各算法重建圖像相同ROI區(qū)域的放大效果圖。之所擇圖2(a)中標(biāo)紅的區(qū)域為目標(biāo)放大區(qū)域，因為此區(qū)域包含了豐富的圖像細(xì)節(jié)和邊緣信息，放大之后便于區(qū)別各算法在保持重建圖像邊緣和細(xì)節(jié)信息方面的能力(彩圖見期刊電子版)。對比其他算法的重建圖像和本文算法的重建圖像可以看本文算法在去除噪聲的同時很好地保留了圖像的細(xì)節(jié)和邊緣信息。

圖2 仿真數(shù)據(jù)主觀實驗結(jié)果.Fig.2 Subjective experimental result of simulation data

(b) 客觀結(jié)果

表1，表2利用各種客觀圖像質(zhì)量評價指標(biāo)RMSE，CC，SSIM，UQI，PSNR及各算法的收斂次數(shù)ITE來對重建圖像性能進(jìn)行輔助評價。

表1 仿真數(shù)據(jù)客觀實驗結(jié)果(I0=100 000)

Tab.1 Objective experimental results of simulation data.(I0=100 000)

算法RMSECCUQIPSNRITEFDK0.067 80.9560.95471.501 41SART0.076 10.9470.93870.501 734BM3D+FDK0.068 00.9560.95371.476 71文中算法0.052 90.9730.97273.664 37

表2 仿真數(shù)據(jù)客觀實驗結(jié)果(I0=50 000)

Tab.2 Objective experimental result of simulation data.(I0=50 000)

算法RMSECCUQIPSNRITEFDK0.076 80.9430.94270.426 41SART0.084 30.9360.92169.612 425BM3D+FDK0.071 10.9520.94971.095 91文中算法0.061 70.9650.96172.322 16

表1，表2分別為在I0=100 000和I0=50 000的輻射劑量情況下的投影數(shù)據(jù)添加混合高斯/泊松噪聲后各算法重建head phantom模體中取第45個切片的各項客觀評價指標(biāo)。其中兩組數(shù)據(jù)所添加混合噪聲中泊松噪聲的大小與I0大小成負(fù)相關(guān)關(guān)系，高斯噪聲相同：σ=50。對比各算法重建圖像的各客觀評價指標(biāo)可以看出，不同輻射劑量不同強(qiáng)度混合高斯/泊松噪聲情況下本文算法重建圖像的同一切片誤差最小，例如，相對于其他方法而言，PSNR最高可以提升2.1 dB；同時和原圖的相關(guān)系數(shù)最高，結(jié)構(gòu)相似性最高，UQI指標(biāo)和信噪比也都是最高的。且圖像質(zhì)量隨著噪聲強(qiáng)度的增強(qiáng)退化越來越嚴(yán)重的情況下本文算法仍然能保持較高的重建圖像質(zhì)量。從計算復(fù)雜度而言，由于BM3D需要進(jìn)行圖像塊匹配，因此，BM3D+FDK具有最高的計算復(fù)雜度；本文算法的計算復(fù)雜度次之，F(xiàn)DK具有最低的計算復(fù)雜度。從算法收斂性而言，對比迭代算法SART達(dá)到收斂時需要幾十次的迭代次數(shù)，本文算法收斂的迭代次數(shù)均在10次以內(nèi)，有很高的重建效率。雖然本文算法相對于經(jīng)典的FDK算法而言，算法的計算復(fù)雜度有所提升，迭代次數(shù)有所增加，會對快速重建造成一定的影響；但是，快速重建問題可以通過基于GPU的并行加速算法來解決，因此并不會對算法的實際使用造成太大的影響。

3.2 真實數(shù)據(jù)實驗

實驗數(shù)據(jù)采用美國瓦利安(Varian)公司新一代250×200 mm2非晶硅平板探測器PaxScan 2520DX所獲取的投影數(shù)據(jù)，其分辨率為768×960，共有360張投影數(shù)據(jù)。重建三維圖像后取其不同的切片如圖3所示。對比重建方法包括FDK，SART，BM3D+FDK及本文算法。對比可以看出，傳統(tǒng)的FDK及SART算法在重建過程中容易受到噪聲的干擾，且隨著噪聲的增強(qiáng)重建圖像質(zhì)量退化嚴(yán)重，因此，重建圖像質(zhì)量較差。BM3D+FDK算法對投影數(shù)據(jù)進(jìn)行去噪，可以有效地提升重建圖像的質(zhì)量，但是同時也會使得重建圖像邊緣模糊。而本文算法在重建過程中有效地去除了噪聲，并在很大程度上保留了圖像的細(xì)節(jié)信息，具有較高重建圖像質(zhì)量。

圖3 真實數(shù)據(jù)實驗結(jié)果Fig.3 Reconstructed data of real projection data

4 結(jié) 論

本文針對低劑量下CBCT圖像重建過程中的噪聲有效去除問題，提出了一種基于高斯/泊松最大似然函數(shù)的CBCT三維圖像重建算法。首先提出了一種結(jié)合圖像去噪與三維重建的統(tǒng)一模型，包括一個保真項和一個約束項。在該模型中，考慮到實際存在的高斯/泊松混合噪聲，保真項通過最大化高斯/泊松似然函數(shù)，來約束重建圖像與觀測數(shù)據(jù)盡可能的相似。由于重建圖像是三維圖像，因此，可以通過三維全變分作為約束項，來減少噪聲對重建圖像質(zhì)量的影響。本文通過head phantom仿真數(shù)據(jù)和真實投影數(shù)據(jù)進(jìn)行三維重建，并和其他的比較經(jīng)典的重建算法FDK，SART和BM3D+FDK算法進(jìn)行了主客觀對比分析，結(jié)果證明本文算法具有較好的去噪效果，例如，相對于其他方法而言，PSNR最高可以提升2.1 dB，同時本文方法能在很大程度上保留重建圖像的細(xì)節(jié)信息，收斂速度快，重建效率高。本文算法可以在低劑量下重建出較好質(zhì)量的CBCT圖像，可以被廣泛應(yīng)用到各種醫(yī)學(xué)研究領(lǐng)域，例如牙科診斷。下一步的主要工作方向是研究基于GPU的并行加速算法，從而更進(jìn)一步提升本文算法的運行速度。