體層合成的各向異性擴散濾波重建算法

俞龍江，戴道清，鄒魯民

（1.中山大學(xué) 數(shù)學(xué)與計算科學(xué)學(xué)院，廣東 廣州510275；2.珠海友通科技有限公司 醫(yī)療信息系統(tǒng)部，廣東 珠海519085）

0 引 言

體層合成 （tomosynthesis）［1］是一種有限角度圖像重建技術(shù)，其成像角度通常在60度以內(nèi)，與傳統(tǒng)CT 相比劑量明顯減小，能夠極大減少病人所受的輻射劑量和檢查所需的時間，降低醫(yī)療設(shè)備成本，提高了醫(yī)生的診斷效率，具有廣泛的應(yīng)用前景，能產(chǎn)生巨大的經(jīng)濟效益并具有良好的社會效益。目前體層合成重建技術(shù)已經(jīng)開始用于國外的高端X 線設(shè)備上，該技術(shù)僅被少數(shù)公司所掌握，尚未普遍應(yīng)用，但是這項技術(shù)可以應(yīng)用的范圍非常廣泛，在牙科檢查［2］、關(guān)節(jié) 檢 查［3］、乳 腺 癌 早 期 篩 查［4］、胸 部 檢 查［5］、放療定位［6］等方面都可發(fā)揮CT 無法替代的作用。然而，這種有限角度的成像結(jié)構(gòu)會產(chǎn)生切片間的偽影，這種偽影會影響閱片時對病灶的判斷，從而影響診斷的正確率。傳統(tǒng)的濾波后投影重建算法無法消除切片間偽影，而目前迭代類型的重建算法難以保證實時性，尚未得到普遍應(yīng)用。

1 各向異性擴散濾波偽影消除算法

1.1 濾波后投影重建算法

自從1972年體層合成由Grant提出以來，針對這種成像結(jié)構(gòu)提出的各種重建算法，主要分為濾波反投影算法和迭代算法兩大類［7］。濾波反投影算法主要是在平移－疊加的后投影操作基礎(chǔ)上，再利用頻率域斜坡濾波來提高重建圖像的對比度［8］。然而，頻率域斜坡濾波在放大高頻信號以達到提高圖像對比度的同時，對有限角度投影產(chǎn)生的偽影也進行了放大，特別是在體層合成應(yīng)用中，由于投影角度的數(shù)目與360度相比非常小 （通常在60度以內(nèi)），經(jīng)濾波后投影重建后產(chǎn)生的偽影尤為明顯，嚴重影響了重建圖像的質(zhì)量。而迭代重建算法收斂速度較慢，需要上百次迭代才能達到較好的重建效果，不能滿足體層合成重建應(yīng)用的實時性要求，故目前濾波后投影算法依然是提供體層合成功能的醫(yī)療設(shè)備主流廠商的選擇，如Siemens［9］、GE［10］、Philips［11］、島津［12］等。如何抑制體層合成產(chǎn)生的偽影，是各大廠家秘而不宣的申請專利保護的技術(shù)，這固然保護了廠家的利益，但也阻礙了體層合成技術(shù)的應(yīng)用與推廣。

從實際應(yīng)用角度考慮，本文選擇濾波后投影類型的算法進行圖像重建，通過將頻率域斜坡濾波器替換為可抑制體層合成產(chǎn)生的偽影的濾波器，從而滿足體層合成重建的圖像質(zhì)量要求。接下來對本文提出的重建算法進行具體闡述。

重建算法的具體構(gòu)造取決于成像系統(tǒng)幾何結(jié)構(gòu)，下面分別給出錐束CT 與體層合成的成像幾何結(jié)構(gòu)圖進行對比，通過對比來演示體層合成的成像特點及與錐束CT 的異同點。圖1是錐束CT 的成像幾何結(jié)構(gòu)。

圖1 錐束CT 成像幾何結(jié)構(gòu)

圖2是體層合成的成像幾何結(jié)構(gòu)。

圖2 （a）是體層合成系統(tǒng)圖像采集方式的原理演示圖，圖2 （b）是對應(yīng)的運動坐標系演示圖。

由圖1與圖2對比可見，圖1所示的錐束CT 成像結(jié)構(gòu)，成像物體繞軸旋轉(zhuǎn)，而從成像物體坐標系來看，可把X 線光源與探測器視為繞軸旋轉(zhuǎn)，旋轉(zhuǎn)軸與運動軌道平面垂直。圖2 （a）所示的是一種常見的體層合成的成像結(jié)構(gòu)，成像物體靜止，球管與探測器繞物體中心軸旋轉(zhuǎn)，旋轉(zhuǎn)范圍從θmax到－θmax，從圖2 （b）所示的成像物體坐標系來看，旋轉(zhuǎn)軸與運動軌道平面也是垂直的。因此，從成像幾何結(jié)構(gòu)角度來看，可以將體層合成視為有限投影角度范圍（2倍θmax）錐束CT，且旋轉(zhuǎn)軸方向從圖1中沿紙面向上變換到圖2 （b）中從紙內(nèi)指向讀者。

圖2 一種常見的體層合成的成像幾何結(jié)構(gòu)

根據(jù)上述分析得出的結(jié)論，即體層合成可看作有限投影角度范圍的錐束CT，故錐束CT 上的濾波后投影算法同樣適用于體層合成圖像重建。根據(jù)體層合成的成像結(jié)構(gòu)稍作調(diào)整，得到的濾波后投影算法具體如下：

（1）旋轉(zhuǎn)軸變換。旋轉(zhuǎn)軸方向從圖1中沿紙面向上變換到圖2中從紙內(nèi)指向讀者。

（2）后投影。采用后投影算法所遵循的錐束CT 幾何結(jié)構(gòu)進行后投影操作。

（3）濾波操作。由于濾波后投影算法中濾波與后投影操作的順序是可以調(diào)換的，故可采用先進行后投影再濾波的重建算法，這里將采用2個濾波器，分別是各向異性擴散濾波器和切片厚度濾波器，前者用于抑制體層合成由于投影角度不足而產(chǎn)生的偽影；后者用于調(diào)整切片厚度，使切片切換時呈現(xiàn)自然的視覺過渡。

其中，切片厚度濾波器操作在垂直于切片平面的坐標軸方向，目的是對切片厚度進行調(diào)整，使感興趣區(qū)域清晰呈現(xiàn)于特定切片內(nèi)，且切片間的過渡要足夠平滑，不至于使人眼在閱片時感覺變化太突然。切片厚度濾波器可以選擇窗函數(shù)進行，例如Siemens方法的漢寧 （Hanning）窗函數(shù)［9］，以及島津方法的高斯窗函數(shù)［12］。其中，高斯窗函數(shù)參數(shù)設(shè)置簡單，且能提供切片間的平滑過渡，故本文選擇高斯窗函數(shù)進行切片厚度調(diào)整。

各向異性擴散濾波器是本文提出的算法核心，下面進行具體闡述。

1.2 各向異性擴散濾波算法

在體層合成重建的三維圖像中，感興趣的對象應(yīng)能在其對應(yīng)的某一個或某幾個切片圖像中清晰可見，而在其它切片圖像中模糊不清甚至是完全不可見。然而，體層合成具有的有限角度投影的特點，由其本身的圖像采集幾何結(jié)構(gòu)決定了切片間的偽影不可能完全消除，即在某一切片圖像中應(yīng)該不可見的對象，在該切片圖像中產(chǎn)生偽影。這種偽影會遮擋該切片圖像中本應(yīng)清晰呈現(xiàn)的對象，造成模糊，還會使本應(yīng)在其它切片圖像中的對象的輪廓呈現(xiàn)在該切片圖像中，影響對該對象的準確位置的判斷。切片間偽影來源于體層合成的圖像采集結(jié)構(gòu)，無論采用濾波反投影重建算法還是迭代重建算法，都會在重建結(jié)果中出現(xiàn)這種偽影。

為消除這種切片間偽影，本文提出采用三維各向異性擴散濾波器對重建體數(shù)據(jù)進行偽影消除。使用該濾波器的目的是消除體層合成的切片間偽影，并在濾波操作后保持對象的對比度不會降低。各向異性擴散濾波器的擴散方程通 常 為［13］

式中：u——三維體素空間，t——擴散時間，D——擴散張量。對式 （1）的擴散方程進行迭代前向差分近似得到

式中：k——迭代次數(shù)。擴散張量D 是對圖像進行結(jié)構(gòu)張量的特征分解得到的，結(jié)構(gòu)張量J 由下式得到

式中：Kρ——一個高斯加權(quán)函數(shù)，ρ——該高斯函數(shù)的參數(shù)σ 值，在使用中一般將其設(shè)置為1。對結(jié)構(gòu)張量J 的特征分解得到

式中：λ——特征分解得到的特征值，v——特征值λ 對應(yīng)的特征向量。從而根據(jù)式 （4）得到擴散張量D 為

式中：D——對稱矩陣，滿足

將式 （6）代入式 （1），可得到等式右邊的散度運算為

式中：ix、iy、iz——x、y、z方向上的通量

體層合成圖像重建的要求是消除體層合成設(shè)備在圖像采集過程中產(chǎn)生的切片間偽影，這意味著濾波器要有一定平滑能力，而同時要對濾波操作造成的平滑進行增強，這又意味著濾波器要有一定的增強能力。綜上，選擇邊緣增強型各向異性擴散濾波器，其特征值滿足［14］

由于切片間偽影在切片圖像中表現(xiàn)為一種模糊的偽影，并隨著切片圖像數(shù)目增加而緩慢衰減，因此各向異性擴散濾波器根據(jù)其濾波器原理可通過擴散方式，將這種模糊偽影通過切片圖像內(nèi)和切片圖像間的擴散，將其逐漸減弱，從而達到消除切片間偽影的目的。同時由于選擇的各向異性擴散濾波器是邊緣增強型的，故不會因為濾波操作而造成對象邊緣的平滑，保持了濾波后圖像的對比度。

2 仿真實驗

在X 線體層合成圖像重建系統(tǒng)設(shè)計和開發(fā)過程中，通常要進行一系列試驗來測量系統(tǒng)各項性能指標，衡量各種影響圖像重建質(zhì)量的因素并驗證設(shè)計的有效性。由于影響X 線體層合成圖像重建的因素較多，整個試驗過程非常耗時且成本極高，因此圖像重建算法的開發(fā)階段采用計算機仿真技術(shù)，可避免各種擾動因素的疊加，從而降低處理難度，加速研究工作的進展。這種仿真技術(shù)是根據(jù)體層合成成像系統(tǒng)的物理模型，并對其加入各種影響成像質(zhì)量的因素模型進行數(shù)字仿真，例如掃描機構(gòu)的機械運動不穩(wěn)定、采集系統(tǒng)幾何結(jié)構(gòu)的偏差、散射噪聲、射束硬化偽影、金屬偽影等。本文在實驗中采用珠海友通公司開發(fā)的體層合成數(shù)字仿真平臺，生成符合體層合成圖像采集結(jié)構(gòu)的投影數(shù)據(jù)，進行重建算法的設(shè)計。仿真實驗中的體層合成系統(tǒng)的圖像采集及重建參數(shù)見表1。

仿真實驗采用一個數(shù)字生成的體模，模擬人胸腔及其內(nèi)部的肺部結(jié)節(jié)，如圖3所示。

表1 體層合成仿真實驗參數(shù)設(shè)置

圖3 人胸腔及肺部結(jié)節(jié)的數(shù)字仿真體模

對人胸腔及肺部結(jié)節(jié)的數(shù)字仿真體模按表1所示的采集及重建參數(shù)進行投影，對生成的投影數(shù)據(jù)進行圖像重建，重建算法選擇常用的濾波后投影算法，得到如圖4所示的重建結(jié)果。

圖4 濾波后投影重建得到的切片圖像

圖4是重建結(jié)果中的一張切片圖像。從圖4 中可見，該切片圖像中呈現(xiàn)出3個圓拱形的偽影，是體模中的肋骨在該切片圖像上形成的切片間偽影。在沒有先驗知識的情況下，無法知道這是偽影，還是組織本來就應(yīng)呈現(xiàn)的影像。

在使用本文提出的算法進行實驗時，需要對算法參數(shù)對性能的影響進行分析，以便正確設(shè)置算法參數(shù)。

各向異性擴散濾波算法的參數(shù)主要包括式 （1）中的擴散時間t和式 （9）中的控制增強對比度的參數(shù)λe，參數(shù)設(shè)置及原則如下：

擴散時間t是控制各向異性擴散濾波器的平滑程度的參數(shù)，隨著擴散時間的增加，圖像將變得越來越平滑，對象與背景間的對比度會越來越低。根據(jù)圖像總體對比度選擇一個合適的擴散時間，使擴散操作后保證組織對比度滿足正常的閱片要求。

參數(shù)λe是用于平滑切片間偽影的門限值，梯度小于這個門限值的圖像區(qū)域，各向異性擴散濾波器的平滑程度隨梯度值減小而增大，即這些區(qū)域的梯度值越小，平滑程度就越大，而梯度大于這個門限值的圖像區(qū)域，各向異性擴散濾波器的平滑程度隨梯度值增大而減小，即這些區(qū)域的梯度值越大，平滑程度就越小，通常存在切片間偽影區(qū)域的梯度小于組織區(qū)域的梯度，為消除切片間偽影，設(shè)置平滑梯度門限值參數(shù)λe在切片間偽影區(qū)域的梯度和組織區(qū)域的梯度數(shù)值之間。如果設(shè)置該門限值越接近切片間偽影區(qū)域的梯度，組織區(qū)域受濾波器平滑操作的影響就越小，即正常區(qū)域的對比度保持的就越好；如果設(shè)置該門限值越接近組織區(qū)域的梯度，切片間偽影受濾波器平滑操作的影響就越大，即切片間偽影的消除程度就越大，圖像就越平滑。一般選擇以保持圖像內(nèi)組織的對比度為前提的情況下進行切片間偽影的平滑，例如在切片圖像總體梯度值范圍歸一到0到1之間后，切片間偽影區(qū)域的梯度數(shù)值處于0.02以下的范圍，則可把該門限值設(shè)置為0.05，略大于切片間偽影區(qū)域的梯度數(shù)值即可。

圖5是采用本文提出的算法重建得到的結(jié)果。圖5是與圖4所示的相同位置的切片圖像。在實驗中，取擴散時間t＝50，梯度門限值λe＝0.02。從圖5中可見，本文提出的算法處理后得到的切片圖像，與圖4所示的傳統(tǒng)的濾波后投影重建結(jié)果相比，切片間偽影得到顯著的消除，圖像中的對象與周圍背景過渡自然，從而提高了重建圖像質(zhì)量。

圖5 各向異性濾波重建結(jié)果

為定量地對比切片間偽影的消除性能，這里定義目標可見度這一指標來衡量切片間偽影的消除結(jié)果。設(shè)切片中某一目標區(qū)域內(nèi)的平均灰度值為O，其周圍背景的平均灰度值為B，且滿足0＜O＜B＜1，則目標可見度D 表示為

根據(jù)體層合成重建原理，在理想情況下，某一目標應(yīng)在某一特定切片最為清晰，可設(shè)在這一切片內(nèi)的目標可見度為1，隨著切片數(shù)目增加，該目標在其它切片內(nèi)應(yīng)迅速變得模糊，即O 顯著變大，B 與O 在數(shù)值上迅速接近，使該目標可見度D 迅速衰減為0。而實際重建中，由于存在切片間偽影，該目標可見度D 衰減緩慢，在其它切片中仍可見到該目標的影子，D 的數(shù)值越低，衰減越慢，對應(yīng)的切片間偽影就越強。

下面通過實例來演示目標可見度是如何計算得到的。圖6是濾波后投影重建算法和本文提出的各向異性濾波重建算法得到的第60張切片：圖6 （a）濾波后投影重建得到的切片，圖6 （b）是各向異性濾波重建得到的切片。圖6（a）和圖6 （b）中位于圖像下方用矩形框標出的圓形黑點，是體模中對人肺部結(jié)節(jié)的模擬，選擇該目標進行目標可見度計算，根據(jù)體模的預(yù)先設(shè)定，該目標應(yīng)只存在于第60張切片內(nèi)，故設(shè)在第60張切片內(nèi)的該目標可見度為1。

圖6 2個算法重建得到的第60張切片

圖7是濾波后投影重建算法和本文提出的各向異性濾波重建算法得到的第100張切片：圖7 （a）濾波后投影重建得到的切片，圖7 （b）是各向異性濾波重建得到的切片。圖7 （a）和圖7 （b）中位于圖像下方用矩形框標出的區(qū)域，是選擇的肺結(jié)節(jié)目標對應(yīng)的偽影區(qū)域，在這2個矩形區(qū)域中分別按式 （10）計算目標可見度D，同理，分別對2個算法重建得到的其它切片計算目標可見度，得到該目標從第60張切片到第100張切片對應(yīng)的目標可見度結(jié)果，如圖8所示。

圖7 2個算法重建得到的第100張切片

圖8是按式 （10）分別計算濾波后投影算法與本文提出的各向異性擴散濾波算法的目標可見度，以對比2個算法對切片間偽影的消除結(jié)果。

圖8 偽影可見性對比結(jié)果

從圖8中可見，橫坐標是從第60張切片到第100張切片的切片序數(shù)，縱坐標是目標的可見度，該目標在重建結(jié)果中應(yīng)只存在于第60張切片中，因此在第60張切片的目標可見度為1，隨著切片數(shù)的增加，該目標的可見度應(yīng)迅速衰減為0。圖6中實線為濾波后投影重建算法得到的結(jié)果，從圖6中可見，隨著切片數(shù)目的增加，目標可見度的衰減緩慢，直到第100張切片時，該目標的可見度仍有0.3，從而造成該目標在其它切片內(nèi)仍依稀可見，存在明顯的切片間偽影。相比而言，虛線所示的本文提出的算法結(jié)果，可隨切片數(shù)目增加顯著地降低該目標的可見度，到100張切片時，目標可見度已衰減到接近0，從而顯著減少了切片間偽影。

3 結(jié)束語

本文提出用于體層合成的各向異性擴散濾波重建算法，經(jīng)過實驗驗證，可顯著消除切片間偽影，與傳統(tǒng)的濾波后投影重建算法相比，圖像質(zhì)量得到明顯改善，從而確保了更好的閱片效果。

目前實驗在數(shù)字仿真的體模上進行，有利于加入各種影響偽影消除的因素比如噪聲、運動模糊、高衰減對象等，來驗證算法性能，數(shù)字仿真實驗同時可對重建算法進行不斷完善，使其應(yīng)用于實際采集圖像成為可能。

［1］Karellas A，Lo JY，Orton CG.Point／counterpoint：Cone beam x－ray CT will be superior to digital X－ray tomosynthesis in imaging the breast and delineating cancer ［J］.Medical Physics，2008，35 （2）：409－411.

［2］Cho MK，Kim HK，Kim SS，et al.Development of dental tomosynthesis system ［C］／／IEEE Nuclear Science Symposium，2007：3788－3791.

［3］Kalinosky B，Sabol JM，Piacsek K，et al.Quantifying the tibiofemoral joint space using X－ray tomosynthesis［J］.Medical Physics，2011，38 （12）：6672－6682.

［4］Sahiner B，Chan HP，Hadjiiski LM，et al.Computer－aided detection of clustered microcalcifications in digital breast tomosynthesis：A 3D approach ［J］.Medical Physics，2012，39（1）：28－39.

［5］Wang J，Dobbins JT III，Li Q.Automated lung segmentation in digital chest tomosynthesis［J］.Medical Physics，2012，39（2）：732－741.

［6］Brunet－Benkhoucha M，Verhaegen F，Lassalle S，et al.Clinical implementation of a digital tomosynthesis－based seed reconstruction algorithm for intraoperative postimplant dose evaluation in low dose rate prostate brachytherapy ［J］.Medical Physics，2009，36 （11）：5235－5244.

［7］Baker JA，Lo J.Y.Breast tomosynthesis：State－of－the－art and review of the literature ［J］.Academic Radiology，2011，18 （10）：1298－1310.

［8］Gengsheng Lawrence Zeng.Medical Image Reconstruction：A conceptual tutorial［M］.NY：Springer，2010.

［9］Zhao B，Zhao W.Three－dimensional linear system analysis for breast tomosynthesis［J］.Medical Physics，2008，35 （12）：5219－5232.

［10］Deller T，Jabri KN，Sabol JM，et al.Effect of acquisition parameters on image quality in digital tomosynthesis ［C］／／Proc SPIE，2007：1－11.

［11］Erhard K，Grass M，Hitziger S，et al.Generalized filtered back－projection for digital breast tomosynthesis reconstruction［C］／／Proc SPIE，2012：1－7.

［12］Gomi T，Hirano H，Umeda T.Evaluation of the X－ray digital linear tomosynthesis reconstruction processing method for metal artifact reduction ［J］.Computerized Medical Imaging and Graphics，2009，33 （4）：267－274.

［13］WANG Dakai，HOU Yuqing.Image processing based on partial differential equations［M］.Science Press，2008 （in Chinese）. ［王大凱，侯榆青.圖像處理的偏微分方程方法［M］.科學(xué)出版社，2008.］

［14］Mendrik A，Vonken E，Rutten A，et al.Noise reduction in computed tomography scans using 3－D anisotropic hybrid diffusion with continuous switch ［J］.IEEE Transactions on Medical Imaging，2009，28 （10）：1585－1594.