一種投影數據恢復的扇束CT局部重建算法

梁亞星，張 權,，張 萌，陳 燕，桂志國,

（1.中北大學 電子測試技術國家重點實驗室，山西 太原 030051; 2.中北大學 生物醫學成像與影像大數據山西省重點實驗室，山西 太原 030051; 3.中北大學信息與通信工程學院，山西 太原 030051）

0 引 言

X-ray計算機斷層成像技術（CT）是一種重要的現代醫療診斷手段，在疾病的早期篩查和后續的診斷治療等方面發揮著重要作用。新型冠狀病毒肺炎（COVID-19）傳染性強、潛伏期長。而高分辨率胸部CT敏感度高，能檢測出毫米級病灶[1]，目前國家衛生健康委員的診療方案中已將CT作為診斷COVID-19的常規首選影像學手段。

鑒于同一患者在診斷過程中往往需要進行多次CT復查，且大多數患者免疫力降低，不宜接受過多的X線輻射[2-3]，故在掃描過程中降低X射線劑量是有必要的。基于感興趣區域的局部斷層掃描成像技術由于可以較好地控制輻射劑量，因而在臨床治療中具有應用推廣價值。該技術在掃描過程中僅覆蓋可疑的病灶部位，避免了對正常區域的不必要射線照射，同時硬件設備的運行損耗也得到降低。但由于采用局部掃描的方式，投影數據被截斷，在重建的局部圖像邊界存在高亮的圓環，且含有灰度移位偽影，從而降低了重建圖像的質量。

為有效消除截斷偽影，傳統方法在擴充缺失數據時多將投影數據進行外插延拓，基本思想是依據某一投影角度測量的已知投影數據對截斷投影數據進行擴充。其中，陳云斌[4]、王克軍[5]將某一投影角度下處于邊緣的兩個投影數據作為常數，分別向兩側擴充，直至填充整個數據區域；薛少峰[6]利用截斷數據邊界內已知的若干連續投影數據進行曲線擬合，并根據該擬合曲線值向外插值，獲取缺失的投影數據。此外，王浩[7]則是根據截斷投影數據相鄰部分的已知投影數據，利用級比生成法擴充缺失數據。

此外，近年來深度學習由于其強大的圖像特征提取能力，逐漸被應用于提高低劑量CT圖像質量[8-11]。其中，Li等[10]對于密集和稀疏角度采樣的截斷投影數據，提出采用神經網絡擬合FBP過程，實現了CT的精確重建；王蕾[11]則利用深度生成網絡預測不完全投影數據的缺失部分，進一步重建圖像。同時，深度學習處理圖像的效果往往也受限制于低劑量CT圖像質量退化程度、訓練樣本數目以及算法的計算復雜度，相比較之下，傳統的局部重建方法重建速度較快，應用范圍較廣，仍然具有較大的研究價值。

為此，受文獻[12]基于線過程模型的平滑性約束集投影重建圖像，文獻[13]根據變分偏微分方程實現圖像修復啟發。基于截斷投影數據的特點，以及線過程模型描述圖像的邊緣結構特性，本文將線過程模型和全變分模型引入截斷投影數據恢復過程中，提出了一種新的局部重建算法。

1 算法的引出

圖1 局部重建的投影數據

故對截斷數據進行恢復時，需要依據已知的A部分投影數據外推截斷邊界以外的B部分投影數據。全變分模型最初應用于圖像去噪，后被逐漸應用于修復不同程度的破損圖像，因截斷數據恢復過程本質上也是一個借助待恢復區B的鄰域信息A恢復較為完整的投影數據問題，本文研究了將全變分模型用于恢復截斷投影數據的有效性。

由于待恢復投影數據點存在鄰域數據缺失的問題，利用全變分模型恢復截斷數據時存在不足，需要合理利用已知的投影數據值對待恢復數據值的鄰域進行預插值。

2 本文算法

算法的基本思想包括兩部分：根據已知的截斷投影數據預估待修復點鄰域的投影數據；將線過程模型檢測出的截斷數據邊緣信息作為恢復截斷數據過程中的擴充邊界，并利用全變分模型恢復缺失的投影數據。基于恢復后的投影數據，進一步濾波反投影重建得到CT局部圖像。

2.1 聯合線過程和全變分模型的截斷投影數據恢復

2.1.1 邊緣檢測

其中閾值Th決定線過程對于投影數據邊緣的敏感度，閾值小，則被判斷存在線過程的區域越大，對截斷投影數據的邊緣信息保存就越多。

故在利用全變分模型對截斷數據進行擴充過程中，需要對判定線過程區域的閾值Th進行自適應調節。在擴充過程的起始一段時間，閾值設置較大，以便快速獲取截斷數據的邊緣位置，此時依據全變分模型擴充后的投影數據值較小；自適應調節閾值大小，經過一定迭代次數后，使得經由全變分模型擴散的投影數據與截斷數據邊界處的投影數據值逐漸接近，從而達到截斷數據恢復的目的。

2.1.2 全變分模型用于投影數據恢復的存在問題

根據偏微分方程理論中的熱傳導數學模型，將投影數據的灰度值視為平面物體的溫度，則投影數據恢復過程，可視為由截斷數據鄰域已知的投影數據作為固定點提供熱量的熱傳遞過程。由于全變分圖像修復模型采用各向異性的信息擴散方式，可以較好地保護邊緣特征和細節信息。因此，在獲取截斷投影數據的邊界位置后，采用全變分模型擴充缺失的投影數據。

全變分模型修復圖像的代價函數為

運用Language乘子法將有約束的極值問題轉化為

式中：r——正的實函數；

u——圖像的像素值；

λ——Language乘子。

引入時間量t，采用梯度下降法求解，則得到最速下降方程為

其中?·、?分別代表散度和梯度算子。

根據Euler-Lagrange方程，將偏微分采用差分代替，引入半點中心插值法求解函數u。若采用3×3像素鄰域，則獲得求解公式如下：

其中已知投影數據待恢復點o周圍的四鄰域像素點、四鄰域半像素點，求解過程中循環次數、平滑系數和步長采用控制變量法，經過多次實驗確定投影數據恢復的參數最適值。因局部掃描下投影數據不完整，待修復的像素點處于邊界位置時，其3×3像素鄰域的部分像素點是不完整的，無法直接利用公式（6）進行求解，故應合理利用已知的投影數據值，對待修復點的鄰域點進行預插值。

2.2 待恢復點鄰域預插值

由于待恢復的投影數據灰度值僅與其鄰域內像素有關，與較遠像素位置無關。全變分恢復過程采用3×3像素鄰域時，需充分利用鄰近的截斷投影數據灰度值，對待恢復點像素的鄰域點進行預插值。

假設ui,j為一個待恢復點，m為待恢復投影數據所在的探測器位置，如圖2所示。由于第m、m-1行探測器投影值為零，采用3×3像素鄰域恢復中心像素點ui,j時不免代入零值，降低數據恢復的效率，故第m-1行探測器的數據值將通過第 m+1行已知的3個投影灰度值加權獲得，即

圖2 待恢復點ui,j及其鄰域

其中權重wi根據待插值點ui與已知點u0的歐氏距離比確定，即

對待恢復點ui,j的3×3像素鄰域仍存在的未知投影數據值ui,j+1，可通過上次全變分模型恢復的投影數據值ui,j-1等值估計，即

到此，待恢復點 ui,j的3×3像素鄰域點已完成預插值，則可根據全變分模型的公式（6），求解待恢復點ui,j的投影數據值。則基于全變分框架的截斷投影數據恢復算法框圖如圖3所示。

圖3 基于全變分恢復截斷投影數據的局部重建算法

綜上所述，所提算法具體步驟如下：

1） 截斷數據獲取：將Shepp-Logan體模進行360°全局掃描，得到圖4（a）所示完整的投影數據，對完整數據截斷得到局部的掃描數據，如圖4（b）所示；

2） 投影數據邊緣檢測：根據公式（2），利用線過程模型檢測出的截斷投影數據的邊緣信息作為擴充條件，圖4（c）為截斷投影數據垂直方向的線過程，即檢測到的待恢復投影數據的邊界位置；

3） 待恢復點預插值：在此基礎上應用公式（7）、（9）對待恢復點像素的3×3像素鄰域點進行預插值；

4） 應用公式（6）對截斷投影數據進行恢復，恢復后的投影數據圖4（d）所示；

圖4 仿真模型投影數據算法過程圖

5） 依據恢復后的投影數據，經濾波反投影重建得到CT局部重建圖像。

3 實驗結果與分析

實驗中，采用Shepp-Logan頭部體模和實際骨盆數據，選取以圖像中心為圓心、半徑為50像素圓形區域為待重建感興趣區域，重建圖像為256像素×256像素，驗證算法進行局部重建的可行性。

為了進行對比驗證，將截斷數據分別進行常數延拓、對稱延拓，并重建圖像。其中，常數延拓將截斷投影數據上、下邊緣處的投影數據作為常數向兩側擴充；對稱延拓則是以截斷數據邊界處的數據作為對稱點，以對稱的方式，分別對上、下邊界之外的缺失數據進行擴充。

此外，為了定量評價算法的有效性，采用歸一化均方距離判據NMSD、歸一化平均絕對距離判據NAAD、均方誤差MSE、峰值信噪比PSNR和結構相似性SSIM進行進一步驗證：

式中：ui,j和vi,j——原始圖像和重建圖像第i行、j列的像素灰度值；

μu和μv——原始圖像和重建圖像的平均灰度值；

σu和σv——原始圖像和重建圖像的標準差，圖像的像素為N×N個。

3.1 仿真數據實驗

實驗采用256像素×256像素的Shepp-Logan仿真模型。采用等角度扇束掃描，投影角度數為360個，探元數目246個，選取中間68個探元數據作為截斷投影數據。圖5（a）、（b）分別是原始體模及其ROI放大顯示圖。

圖5 仿真模型

圖6左列為截斷數據，以及截斷投影數據經常數延拓、對稱延拓和本文算法三種算法擴充后FBP重建結果，右列分別為局部重建結果分別與原圖相應的感興趣區域圖5（b）的差值圖。顯然，截斷數據局部重建圖6（a）和對應差值圖6（b）的邊界處都存在環狀的高亮截斷偽影，且重建圖像灰度值產生畸變，邊界區域像素的灰度值偏高。與其他截斷數據擴充后的FBP局部重建圖比較，本文算法重建圖6（g）抑制了邊界處的環狀截斷偽影，且圖像的對比度提高、均勻性較好。比較重建圖像差值圖，可以看出圖6（b）、圖6（d）和圖6（f）的左下角方框區域內存在部分程度損失的邊緣信息，相比之下圖6（h）中的邊緣殘留更少，進一步說明本文算法重建局部圖像的邊緣結構信息保存的更加完整。

圖6 左列為不同算法FBP局部重建圖，右列為ROI差值圖

為了更加客觀評價本文算法，表1列出了各局部重建算法定量評價參數。分析比較可見，本文算法的NMSD、NAAD和MSE值較低，表明重建圖像與ROI差異程度小，PSNR值較高表明失真程度小，SSIM數值表明本文算法局部重建圖像與ROI的結構相似程度高，從而驗證了本文算法恢復截斷投影數據的有效性。

表1 仿真數據經不同算法恢復后的局部重建圖像定量評價參數

3.2 實際數據實驗

為了進一步驗證算法的實用性，實驗采用實際骨盆數據。投影角度數360個，探元數目512個，選取中間180個探元數據為截斷投影數據，如圖7所示。

圖7 實際骨盆

圖8（a）、（b）分別為完整投影數據和截斷投影數據，圖8（c）為采用本文算法恢復后的投影數據。

圖8 實際骨盆投影數據

圖9左列為截斷投影數據經不同算法擴充后FBP重建結果，右列為對應差值圖。圖9（a）、（b）的邊界處都存在環狀亮偽影，這與仿真模型截斷投影數據的局部重建結果邊界處存在的高亮環狀偽影相一致。將局部重建結果圖與感興趣區域圖7（b）比較，圖9（e）的矩形框標注內存在明顯的圖像信息缺失，上下邊界處重建效果較差，這在定量評價參數中也得到體現。而本文算法局部重建圖像圖9（g）的邊界處環狀截斷偽影得到抑制，視覺效果比較好，對比差值圖9（h）也存在更少的細節信息。

圖9 重建圖與插值圖（左列為不同算法FBP局部重建圖，右列為ROI差值圖）

圖10為對應算法局部重建結果圖與完整數據重建圖像ROI的結構相似索引圖。比較圖10（a）～（d），本文算法處理的圖10（d）矩形框標注區域內的像素亮的程度比較高，進一步表明本文算法重建的局部圖像質量較高。

圖10 不同算法重建圖的結構相似索引圖

表2列出了局部重建算法定量評價參數，與其他截斷數據擴充方法比較，本文算法NMSD、NAAD和MSE值降低，PSNR值和SSIM值有明顯的提高。本文算法局部重建圖像的有效性進一步得到了驗證

表2 實際數據經不同算法恢復后的局部重建圖像定量評價參數

4 結束語

本文提出了一種基于投影數據恢復的扇束CT局部重建方法。利用線過程模型檢測的截斷投影數據邊緣信息作為擴充邊界；在恢復過程中，針對待恢復點像素的鄰域點的數據缺失問題，利用已知投影數據的距離比加權進行預插值，然后依據全變分模型對中心投影數據值進行恢復，改善了直接使用截斷數據重建圖像造成的截斷偽影問題。對比其他截斷投影數據擴充方法，從仿真數據和實際數據的局部重建圖像的視覺效果及客觀評價參數等方面，評估了算法的可行性。