俞龍江，戴道清，鄒魯民

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

0 引 言

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

1 各向異性擴(kuò)散濾波偽影消除算法

1.1 濾波后投影重建算法

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

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

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

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

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

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

由圖1與圖2對(duì)比可見，圖1所示的錐束CT 成像結(jié)構(gòu)，成像物體繞軸旋轉(zhuǎn)，而從成像物體坐標(biāo)系來(lái)看，可把X 線光源與探測(cè)器視為繞軸旋轉(zhuǎn)，旋轉(zhuǎn)軸與運(yùn)動(dòng)軌道平面垂直。圖2 （a）所示的是一種常見的體層合成的成像結(jié)構(gòu)，成像物體靜止，球管與探測(cè)器繞物體中心軸旋轉(zhuǎn)，旋轉(zhuǎn)范圍從θmax到－θmax，從圖2 （b）所示的成像物體坐標(biāo)系來(lái)看，旋轉(zhuǎn)軸與運(yùn)動(dòng)軌道平面也是垂直的。因此，從成像幾何結(jié)構(gòu)角度來(lái)看，可以將體層合成視為有限投影角度范圍（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)進(jìn)行后投影操作。

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

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

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

1.2 各向異性擴(kuò)散濾波算法

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

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

式中：u——三維體素空間，t——擴(kuò)散時(shí)間，D——擴(kuò)散張量。對(duì)式 （1）的擴(kuò)散方程進(jìn)行迭代前向差分近似得到

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

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

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

式中：D——對(duì)稱矩陣，滿足

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

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

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

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

2 仿真實(shí)驗(yàn)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié)束語(yǔ)

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

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

［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.