基于衰減補償的CBCT短掃描實時重建方法

【作 者】楊坪坪，馮漢升,，許繼偉，楊洋，宋云濤,

1 中國科學技術大學，合肥市，230026

2 中科離子醫學技術裝備有限公司，合肥市，230088

3 中國科學院等離子體物理研究所，合肥市，230031

0 引言

錐形束CT（cone beam computed tomography,CBCT）[1-2]因錐形束射線利用率高，采集效率高，能有效減少輻射劑量，提高成像效率，提高層間分辨率而應用于口腔、頭部、胸腔以及盆腔成像中。在臨床應用中，每日就診人數較多，每次患者掃描進行圖像重建的時間有限，因此對CBCT系統的重建算法運行速度要求較高。基于短掃描的重建方式所需投影幀數更少，重建時間減短，且掃描劑量較少，因此成為當前研究熱點。由于數據冗余，在掃描過程中沒有必要讓探測器做360o全掃描，掃描范圍可以小于360o，這種掃描方式叫做短掃描。基于短掃描加權的FDK（Feldkamp）圖像重建算法[3]既能保證重建速度，掃描劑量也更少，因此商業應用潛力更加巨大。1982年Parker首先提出了圓軌跡短掃描扇形束概念[4]，指出適當地對短掃描原始數據進行加權，重建得到的圖像質量基本上等價于使用全掃描投影數據進行圖像重建得到的圖像質量，并提出了一種二維加權方式，獲得了較好的短掃描重建圖像質量。同時由于錐形束成像采用圓掃描軌跡，投影數據不完全，因此FDK重建圖像存在強度衰減現象[5-6]，研究普遍采用反投影加權進行強度衰減補償，其中倒高斯加權和倒余弦加權公式重建效果最佳，但都各有缺點，倒高斯加權重建水平向衰減補償以及倒余弦加權重建軸向衰減補償不佳。根據已有的衰減補償算法的不足，本研究提出新的衰減補償算法，在提高倒高斯加權軸向衰減補償的同時提高水平向衰減補償，有效地提高了重建圖像質量；并針對基于三維重建開源庫RTK庫（Reconstruction Toolkit，醫學圖像重建庫）實現的短掃描加權FDK算法進行流程改進，克服了基于RTK庫一次讀入所有投影圖像的缺點，實時性地每次讀入一幀投影圖像進行重建，同步并行執行患者采集和三維重建過程，大大減少從圖像開始采集到生成圖像重建之間的總耗時，同時應用新的衰減補償算法改善圖像質量。

1 FDK算法原理

三維圖像重建的定義是從物體的投影數據得到物體的內部斷層圖像的過程，FDK算法是其中應用最廣泛的算法，需要利用投影圖像數據與重建幾何位置信息，實施加權、濾波、反投影三個步驟即可完成重建。

1.1 錐形束圖像重建FDK算法

FDK算法是由Feldkamps等[7-8]在1984年提出，針對錐形束投影數據進行近似重建。該算法可以由扇形束的FBP（濾波反投影）算法推導得到。

針對圓軌跡的FDK算法公式如下:

式中，βmin為初始視角，βmax為最終視角，R為物體中心到射線源距離，Z為探測器平面上探測器像素對應探測器坐標系的距離。p(t,l,β)⊕h(t)行濾波后，經加權β)，即可反投影得到重建圖像。

其中，p(t,l,β)為投影數據，t是探測器上的行坐標，l為探測器上的列坐標，h(t)為一維斜坡濾波核，參數幾何位置關系見圖1。

圖1 錐形束示意圖Fig.1 Schematic diagram of cone beam

1.2 短掃描加權FDK算法

根據1982年Parker首先提出的圓軌跡扇形束短掃描的概念以及加權概念，短掃描重建需要在加權基礎上乘上額外的Parker加權因子。

針對第i幀的Parker加權因子wi(α,β)公式如下：

其中：α是投影數據位置相對于重建物體中心角度，β是X射線源對應角度。△是扇形束角度，其參數幾何位置關系，如圖2所示。

圖2 Parker加權因子中各參數對應幾何關系Fig.2 Geometric relationship in the Parker weighting factor

2 衰減補償算法改進

2.1 算法介紹及改進

FDK重建域中，中心斷層物體能被精確重建，非中心斷層物體的重建精度隨與中心斷層距離迅速下降，其形狀貌似山頂函數，稱為強度衰減現象，如圖3所示。

圖3 均勻密度球模型橫截面512層與FDK重建斷層圖像橫截面對比Fig.3 Comparison of the cross section of the uniform density sphere model (512 layers) and the cross section of the FDK reconstruction central tomographic image

由此，學者們提出許多函數進行強度衰減補償，其中最有效的是倒高斯加權函數[9]和倒余弦加權函數[5]。其一般形式如下:

倒高斯加權函數：

其中：z∈[-zmax,zmax]

倒余弦加權函數：

其中，x、y、z分別為重建點坐標，k、c1、c2分別為倒高斯加權函數和倒余弦加權函數的加權因子，D為射線源到重建中心點距離。根據文獻[9]表明，倒高斯函數最優加權因子k=20.5；倒余弦加權函數最優加權因子c1=1.4，c2=1。倒高斯函數相比較倒余弦函數有更好的軸向衰減補償效果；而倒余弦函數因為考慮了重建點到原點距離r因素，因此有更好的水平向衰減補償效果。

由上，本研究提出新的加權公式，兼顧倒高斯加權函數與倒余弦加權函數的優勢，在倒高斯函數的基礎上加上距離r因素，稱為改進倒高斯加權函數，如下：

2.2 改進算法比較

實驗采用均勻球密度模型，利用Matlab進行仿真重建。重建角度從0o到209o，每隔1o取一張投影，錐角取28.72o。重建體素512×512×512個，體素尺寸0.08 mm×0.08 mm×0.08 mm，探測器像素512×512個，探測器尺寸51.2 mm×51.2 mm，源點到重建物體中心距離80 mm，到探測器中心距離100 mm。在醫學成像的仿真實驗中，常用歸一化均方誤差評價重建圖像質量，定義為:

其中：xi為物體模型真實強度值，為重建圖像像素值。

共設置5組數據，即原始球數據、FDK重建數據、倒高斯加權數據、倒余弦加權數據、改進倒高斯加權數據，其中倒高斯加權、倒余弦加權均采用最優加權因子，改進倒高斯加權采用加權因子k1=20.5，k2=1.4。對比如圖4～圖7所示。

圖4 X軸第256層（均勻球模型)、X軸第244層（SheepLogan[10]模型）模型重建橫截面圖像對比（窗口[0.9 1.25]）Fig.4 X-axis layer 256 (uniform sphere model),X-axis layer 244 (SheepLogan[10] model) model reconstruction cross-sectional image comparison (window [0.9 1.25])

圖5 軸向x=256,y=256掃描線強度對比Fig.5 Axial x=256,y=256 scan line intensity comparison

圖6 水平向y=256,z=128掃描線強度對比Fig.6 Scanning line intensity comparison of horizontal y=256,z=128

圖7 水平向x=256, z=128掃描線強度對比Fig.7 Scanning line intensity comparison of horizontal x=256,z=128

根據掃描線強度對比圖及歸一化均方誤差表（見表1）比較可知，經過衰減強度補償的FDK算法圖像質量都有明顯的提升，其中改進倒高斯加權重建軸向衰減補償優勢更大，比倒余弦重建軸向衰減補償誤差減小了11.019%，比倒高斯重建水平向128層衰減補償誤差減小了1.291%，軸向重建質量幾乎趨近于原始圖像，且水平向衰減補償較倒高斯加權也有較大的提升，接近于倒余弦加權重建，同時豎直面方向上改進倒高斯加權重建衰減補償誤差更小，而水平面上重建衰減補償誤差接近于倒余弦加權重建。因此，改進倒高斯加權重建是一種更有效的衰減補償重建方法。

表1 不同掃描位置下的歸一化均方誤差比較Tab.1 Comparison of normalized mean square error at different scanning positions

3 短掃描FDK算法實施流程及實驗結果對比

利用改進倒高斯加權公式，本研究針對原有基于RTK的短掃描加權FDK重建方法進行了改進。首先，原有方法在開始重建之前即讀入所有投影數據，無法滿足實時性，改進方法通過同步患者采集和重建過程，實現了實時重建；其次，在改進方法反投影的過程中進行改進倒高斯加權，改善圖像質量。基于以上兩點，做了5組性能對比實驗，分別是本地Parker加權重建（原有方法）、未加權實時重建、倒高斯加權實時重建、倒余弦加權實時重建、改進倒高斯加權實時重建（改進方法），濾波方式采用Hamming濾波，截止頻率為0.5。其中本地Parker加權重建過程模擬基于RTK一次讀入所有投影數據方法。

3.1 算法實施流程圖

具體實施流程如圖8所示。

圖8 算法實施流程Fig.8 Algorithm implementation flowchart

3.1.1 本地Parker加權重建流程

（1）等待探測器采集投影圖像時間26.4 s；

（2）一次性加載投影數據和幾何位置信息；

（3）經過加權、濾波、Parker乘積遍歷加權、反投影進行FDK重建；

（4）重建結果輸出，計算重建時間。

3.1.2 未加權實時重建流程

（1）一次性加載幾何位置信息；

（2）模擬患者采集過程，每間隔70 ms將投影數據發送給重建線程；

（3）經過加權、濾波、反投影進行FDK重建；

（4）重建結果輸出，計算重建時間。

3.1.3 衰減補償實時加權實時重建流程

（1）一次性加載幾何位置信息；

（2）模擬患者采集，每間隔70 ms將投影數據發送給重建線程；

（3）經過加權、濾波、Parker乘積遍歷加權、衰減補償函數加權反投影進行FDK重建；

（4）重建結果輸出，計算時間。

3.2 實驗對比

本次測試以CatPhan500性能測試體模（見圖9）代替患者，同時結合投影圖像實例（見圖10）進一步說明重建對比結果。

圖9 CatPhan500性能測試體模Fig.9 CatPhan500 performance test body mask

3.2.1 測試采用的設備性能

所使用的實驗平臺為Intel Core i5-8400 CPU，內存8 GB，GPU為NVIDIA GeForce GTX 1050 Ti，設備上可用的全局內存總量為4 096 MB。

3.2.2 測試所用圖像重建參數

射線源到重建物體中心距離為1 967 mm，射線源到探測器中心距離2 967 mm，錐角8.32o。投影圖像角度從0.5o到189o，間隔0.5o取一幀，其中一幀投影圖像見圖10，每幀像素個數為1 440×1 440，像素大小為0.3 mm×0.3 mm，每個像素占2個字節。重建圖像像素個數為512×150×512，像素大小為0.5 mm×2 mm×0.5 mm，每個像素占2個字節。測試數據采用短掃描數據，重建結果經過射野外視場消除、線性變換從圖像像素值轉換為CT值。

圖10 投影圖像示例Fig.10 Example of projected image

3.2.3 重建結果對比

重建結果對比如圖11所示，可以看出，相比較未加權實時重建結果，加權重建圖像變得更加清晰，圖像偽影明顯變少，成像質量提高了許多，這是由于重建方式為短掃描重建，且有冗余數據，當不進行Parker加權時，有冗余數據的探測器角度會多次反投影，而沒有冗余數據的探測器角度只會進行一次反投影，造成重建圖像質量明顯不一。此外，圖層2在不同重建方法中都可以看到比較明顯的偽影，初步考慮為由于散射造成的邊緣模糊現象。錐形束CT成像中，X射線一次與物體接觸的體積更大，散射發生概率也會更大，因此可以考慮散射矯正進一步改善重建圖像質量。

圖11 短掃描未加權實時重建、本地Parker加權重建、倒高斯加權實時重建、倒余弦加權實時重建、改進倒高斯加權實時重建結果對比 Fig.11 Results comparison of short scan unweighted real-time reconstruction,local Parker weighted reconstruction,Gaussian weighted real-time reconstruction,inverse cosine weighted real-time reconstruction,improved inverse Gaussian weighted real-time reconstruction

從圖11及表2可以看出，5種重建方式中實時重建耗費時間減小接近一半，更好地滿足臨床上對時效性的要求。本地Parker加權重建方法采用一次性加載所有投影圖像數據，缺點是不能實時重建，無法滿足實際需求。未加權實時重建方法重建過程中不計算Parker權重，所以重建時間最短。改進倒高斯加權實時重建實時同步計算Parker加權因子并作用于投影圖像數據，只需要一次遍歷，且在反投影過程中實施改進倒高斯加權，Parker加權實時重建方法重建圖像質量更佳，時間更短。通過模擬患者掃描時間（26.4 s）和改進倒高斯加權實時重建方法重建時間（27.932 s）可以看出，改進倒高斯加權實時重建方法在患者采集完成后1～2 s內即可完成整個重建過程，達到了實時重建、改進重建圖像質量的目的。

表2 五種重建方法耗費時間對比Tab.2 Time consumption comparison of five reconstruction methods

4 結論

通過同步患者采集和三維圖像重建過程，在載入每幀投影數據的同時進行重建，且在反投影過程中進行新的衰減補償加權，克服了基于RTK庫的短掃描加權FDK算法一次載入所有投影圖像無法滿足實時重建的問題，同時能有效改善圖像質量，在患者采集完成后1～2 s內完成重建，實現了短掃描Parker加權FDK算法實時重建。針對基于短掃描的Parker加權FDK算法實時重建問題進行了研究，但仍有許多不足的地方，需要進一步考慮三維圖像質量的改進方法，對圖像進行散射矯正以及設計濾波函數改善圖像質量，更好滿足臨床需求。