低劑量照射條件下CT影像快速重建算法研究

董 建，陳 昊，孟慶寬，秦轉萍，楊耿煌

（天津職業技術師范大學天津市信息傳感與智能控制重點實驗室，天津300222）

隨著計算機斷層成像（computed tomography，CT）技術的發展，圖像質量的明顯改善，CT檢查越來越多地被用于一些復雜疾病的診斷，并已經成為臨床工作中不可缺少的影像學診斷工具[1-2]。然而，CT是一種放射性檢查，具有電離輻射效應，患者過多地接受此類檢查會導致電離輻射在人體內的積累，有引起癌變的風險[3]。因此，加速低輻射劑量的醫用CT器械的研發已經成為醫生與患者的共同期待。降低CT設備的輻射劑量可從改進硬件設備和開發重建算法兩方面取得突破。近年來，CT設備制造商對其關鍵部件，如探測器、準直器、適形濾波器、自動曝光控制器等的技術更新，已經為降低輻射作出了重要貢獻[4-5]。在圖像重建方面，利用稀疏角度或受限角度投影數據、低管電壓下的投影數據和非全局投影數據進行圖像重建是在算法層面降低輻射的3個有效措施。近年來在該領域的論文數量增長迅速，也預示著超低輻射劑量CT的成像算法開發已成為研究熱點之一[6-9]。本文聚焦利用低管電壓下的投影數據進行圖像重建的問題，對含有大量泊松噪聲的投影數據進行重建。

投影數據的不完整性給圖像重建帶來了巨大挑戰，搭載于現行CT機中的濾波反投影（filtered backprojection，FBP）算法重建出的圖像含有大量噪聲或放射狀偽影，嚴重影響了正常組織結構和病變區域的刻畫，甚至基本的迭代型算法如代數式重建（algebraic reconstruction technique，ART）和同步迭代（simultaneous iterative technique，SIRT）算法也無法獲得滿意的重建圖像。2006年，Donoho等[10-11]提出了壓縮感知（compressed sensing，CS）理論，為此類欠定問題的求解提供了理論依據。CS理論的求解機制一般要求先設計一個包含數據保真項和正則化項的目標函數，運用最優化理論實現目標函數的最小化，在此過程中得到所求數據的最優解。正則化項的設計對重建圖像的質量起著至關重要的作用。全變分（total variation，TV）正則化在去除圖像噪聲、保存物體邊緣方面能發揮有效作用，已經被視為該領域的優秀方案[12]。然而，越來越多的研究發現TV模型并不能完美重建所有類型的圖像，尤其在重建含有較多高頻信息的圖像時表現出明顯缺陷，如出現補丁狀偽影、鋸齒狀物體邊緣、丟失圖像細節等[13]。因此，設計更優秀的正則化項已經變得迫在眉睫。近年來，深度學習技術也越來越多地被應用于低輻射劑量CT圖像重建領域[14-15]。

本文創新了目標函數的設計方式，在TV正則化項的基礎上引入基于非線性濾波的補充正則項。新方案使得非線性濾波操作與CS框架有機結合，進一步增強了算法對圖像的重建能力。具體體現在TV模型只考慮中心像素與2個鄰域像素的灰度相似度，而非線性濾波能將考慮范圍擴展至包圍中心像素的更大鄰域范圍。研究驗證了4種非線性濾波分別嵌入CS框架后的新算法的重建效果，確認了新算法的有效性，利用凸優化領域的臨近點算法對目標函數進行了最小化處理，構建出以行結構為基礎的快速迭代算法，并利用腹部圖像的重建驗證了算法的收斂性。

1 腹部CT圖像與投影數據

在人體各部分組織的CT影像中，腹部影像涵蓋了人體的眾多器官，像素的灰度跨度較大，而且圖像中既包含了具有區域灰度一致性的較大面積器官，又包含了血管、胃腸容物、軟組織褶皺等圖像的細節信息。腹部影像可謂同時具備了眾多的低頻和高頻信息，是用來驗證圖像重建算法性能的優質圖像樣本。本實驗選取了腹部斷層作為重建對象，對實際臨床圖像進行了去噪處理的實驗用圖像如圖1所示。圖像的像素數為512×512，像素灰度值反映CT值，其跨度范圍為[-1024HU，2269HU]。

圖1 人體腹部CT圖像

人體腹部CT圖像的投影數據如圖2所示。圖2（a）是圖1所示腹部CT圖像的投影數據，是在[0°，180°]范圍內以180°/270的角度間隔進行的數據采集，因此該投影數據的像素數為512×270。此投影數據是對圖1中圖像進行再投影運算取得的數據，并非CT儀器中導出的圖1圖像重建前的原始投影數據。圖2（a）的投影數據是在以平行束結構建模下取得的。本文聚焦低管電壓下的投影數據的圖像重建，向圖2（a）的投影數據中加入泊松噪聲，即可再現低管電壓下的數據采集結果。圖2（b）為光子數為4.0×106時引入了泊松噪聲的投影數據。低管電壓即代表了低劑量照射條件下的CT檢查。

圖2 人體腹部CT圖像的投影數據

本研究選用C編程語言在Visual Studio 2017運行平臺上進行了實現。所使用的PC機參數為Intel（R）Xeon（R）2.40 GHz 2.40 GHz的CPU處理器，32.0 GB內存，64位Windows 10教育版操作系統。

2 算法介紹

2.1 問題定義

CT圖像重建即利用重建算法將投影數據轉化為CT圖像的過程。依據CS理論，本文設計了如式（1）所示的目標函數

式中：x=（x1，x2，…，xJ）T是一個J維向量，表示待求解圖像；b=（b1，b2，…，bI）T是一個I維向量，代表已知投影數據；矩陣A代表I×J的系統矩陣，由CT裝置決定；是L2范數型數據保真項，其保證待求解向真值逐步逼近是對圖像進行濾波后又進行了稀疏約束；‖x‖TV是使重建圖像具有分片光滑特性。兩項操作均在去除了部分噪聲后保留了圖像的高頻信息。2種范數的展開式為

TV范數定義為

xm，n代表圖像x在（m，n）處的像素值，本文僅在此處將圖像看作一個二維矩陣以期對TV范數做直觀說明。根據定義，式（3）可看作灰度梯度的L1范數，即因此TV范數也是一種特殊的L1范數。中的M操作代表對圖像的非線性低通濾波，該正則化項使非線性濾波操作被引入CS框架，使非線性濾波的優勢得以在CS中發揮重要作用，這也是本研究設計該正則化項的出發點。目標函數中的β和γ代表控制各正則化項作用強弱的超參數。為了使非線性濾波發揮優勢以及抑制TV范數的缺陷，設定β取較大值，確定其主體地位，而γ取較小值，使TV項起輔助作用。

2.2 正則化項中的非線性濾波

本文重點驗證了中值濾波、雙邊濾波、非局部均值濾波和非局部中值濾波分別嵌入正則化項后圖像重建的效果。

中值濾波是一種基本的非線性濾波，它將中心像素點及其鄰域的像素值由小到大排列，用位于序列的中間值代替中心像素的原有像素值。以同樣的模式對圖像中的每一個像素點進行遍歷，完成中值濾波對圖像的處理。中值濾波能有效去除脈沖性噪聲，且在圖像處理任務中能有效地保存圖像中物體的邊緣。

雙邊濾波是一種進化了的均值濾波，其首先計算中心像素及其鄰域像素的權重，再用求出的加權平均值代替中心像素的原有像素值。以同樣的模式對每一個像素點進行遍歷。和傳統的圖像平滑算法不同，雙邊濾波在計算權重過程中除了考慮像素間幾何距離上的接近程度外，還考慮像素間灰度的一致性。這使得雙邊濾波不僅能有效地去除噪聲，還能保存圖像上的邊緣信息。

非局部均值濾波在雙邊濾波的基礎上更強調像素塊間的灰度一致性，非局部均值濾波的搜索窗口與像素塊的示意圖如圖3所示。

圖3 非局部均值濾波的搜索窗口與像素塊的示意圖

黃色部分是以像素xj0為中心的搜索窗口，橙色和藍色部分是分別以xj0和xj0′為中心的像素塊。搜索窗口的中心像素的像素值應被窗口中每個像素的加權平均值替代，因此計算搜索窗口中每個像素的權重成為關鍵環節。權重的大小由2個因素決定，分別為中心像素xj0與搜索窗口中某一像素如xj0′間的高斯距離，另一因素為分別以xj0和xj0′為中心的像素塊的灰度相似性。像素xj0′的權重wj0′的計算公式為

式中：D2（xj0，xj0′）為像素xj0與xj0′間的高斯距離；N為像素塊中的像素數；參數δ1和δ2用來控制距離因素與灰度相似度因素對wj0′的貢獻度。

非局部中值濾波中搜索窗口與像素塊的選取策略以及搜索窗口中像素權重的計算方法與非局部均值濾波一致。當式（6）中的不等式成立時，以像素xjs′代替中心像素xj0，式中S代表搜索窗口中的像素數。

2.3 算法推導

為了構建快速的行處理型重建算法，本研究將目標函數中的數據保真項按系統矩陣的行進行拆分，見式（7）和式（8），并且用函數g（x）和h（x）分別代表2個正則化項。

其中：向量ai代表系統矩陣A的第i個行向量，aiT是向量ai的轉置。目標函數f（x）的最小化過程將用到凸優化理論中的重要工具，即臨近點算法。文獻[16]對該理論進行了詳細介紹，本文不再贅述。由臨近點算法推導得出算法執行過程如下

其中：k代表總循環次數，式中關于x的各種符號代表各中間過程圖像。由式（9）可知，在變量i的循環過程中，任何一個投影數據元素bi均能促成一次圖像的更新。只要程序的執行遵循特定的數據訪問順序，該算法就能像ART算法一樣具有快速收斂性質，文獻[17]對數據訪問順序的設計做了詳細說明，本實驗亦采用該類數據訪問順序。prox算子中的α(k)代表步長參數，為了避免數據陷入極限環，令α(k)=α0/（1+ε·k），α0和ε取值依具體任務確定。為了在圖像平滑和去噪的同時達到程序快速收斂的效果，算法引入了跨度參數P，使得在數據保真項作用的間隙讓2個正則化項依次起作用。式（9）還不具備算法的顯式表達，以下將對算法中各個prox算子作出說明。

根據臨近點算法，關于fi（x）的prox算子展開如下

對式（10）求解，可得到如下迭代式

函數g（x）的prox算子展開如下

其中，Mx的存在（x是待求未知量，無法對其進行濾波操作）阻礙了式（12）的最小化。當步長參數α(k)取值足夠小時，x無限接近已知量x(k，i+1，0，0)。因此，在α(k)取值足夠小的前提下，Mx可由當前圖像的濾波結果Mx(k，i+1，0，0)代替。至此，式（12）的求解可由封閉形式的軟閾值處理實現，可得到迭代式

函數h（x）的prox算子展開如下

式（14）要實現TV項與二次項和的最小化，此問題形式是一個在CS領域常見的最優化問題。可將其看作是以x(k，i+1，1，0)為輸入圖像的TV項約束的圖像去噪問題，被稱作Chambolle投影算法的標準化算法可用于該問題的求解，此處引用了該算法[18－19]。

至此，低劑量照射條件下CT影像的行處理型重建算法已經構建完成。

2.4 算法展示

初始化 輸入像素數為N×M的投影數據b。設置步長參數的控制變量（α0，ε）。設置β>0，γ>0的超參數以及P>0的跨度參數。

步驟1 賦值總循環次數k=0，像素數為N×N的初始圖像x(0)=0。

步驟2 計算步長參數α(k)=α0/（1+ε·k），令x(k，0，0，0)=x(k)。

步驟3 圖像更新程序：

步驟4 令x(k+1)=x(k，i+1，1，1)，返回步驟2。

輸出 當循環條件被滿足時，輸出圖像x(k+1)。

3 實驗結果及分析

3.1 傳統算法重建結果

首先，本文用搭載于現行CT機中的FBP算法對低管電壓采集條件下的投影數據（圖2（b））進行重建，重建結果如圖4所示。結果顯示，有大量噪聲傳遞到腹部的斷層影像，各組織中的細節信息被覆蓋，嚴重損害了該斷層影像的臨床診斷價值。

圖4 FBP算法的重建結果

其次，本文用CS領域的經典算法TV最小化算法進行了圖2（b）的重建實驗，即只允許式（1）中目標函數的數據保真項和TV項起作用，令β=0。步長參數的控制變量α0=10，ε=1 000，跨度參數P=9×512，總循環次數k=50。本文對TV項的超參數γ的取值進行了探索性實驗，代表性結果如圖5所示。圖5（a）為γ=1.0時的重建圖像，圖5（b）為γ=100時的重建圖像。結果表明：γ取值小則TV項無法達到去噪效果，γ取值大則TV項對圖像過度平滑，使圖像細節消失以及出現鋸齒狀邊緣。γ的其他取值亦未帶來滿意的重建效果。

圖5 TV算法的重建結果

3.2 新算法重建結果

圖6為本文提出的基于非線性壓縮感知4種算法的重建結果。

圖6 新算法的重建結果

實驗所選定的4種情況的共同參數：α0=10，ε=1 000，β=1 500，γ=20，P=9×512。圖6展示的是各實驗的第50次的運行結果。圖6（a）的實驗中，所設置的中值濾波器的搜索窗口為5×5。由圖6（a）可知，Median-CS算法對圖像進行了較好的去噪，然而其導致了部分圖像細節的丟失，如肝臟中的血管信息。同時，結果圖像類似于TV算法的結果，出現了部分鋸齒狀邊緣。圖6（b）的實驗中，雙邊濾波的搜索窗口為5×5，設置δ1=10 000，δ2=120。由圖6（b）可知，Bilateral-CS算法使圖像獲得了良好的平滑效果，然而圖像的細節被過度平滑而失真。圖6（c）的實驗中，非局部均值濾波的搜索窗口和像素塊的大小分別為7×7，5×5，設置δ1=δ2=45。由圖6（c）可知，NLWmeans-CS算法使圖像獲得了較好的平滑效果，保存了圖像的細節信息，但在圖像的高頻區域有密集性噪聲殘留，在圖像的低頻區域有分散性噪聲點殘留。圖6（d）的實驗中，非局部中值濾波的搜索窗口和像素塊的大小分別為7×7，5×5，設置δ1=5，δ2=1 000。由圖6（d）知，NLWmedian-CS算法獲得了最好的圖像重建效果，不僅對圖像的噪聲去除較完善，還較好地保留了圖像的細節信息。

圖6（b）—（d）的重建結果中，對圖像的高頻部分出現密集性噪聲殘留的機理闡述如下：在非局部濾波中，設定的搜索窗口中每個像素的權重由2部分因素決定，除了與中心像素的距離因素相關外，還與中心像素塊的灰度相似度因素有關。距中心像素越遠、與中心像素塊的灰度值差距越大，該像素獲得的權重越小。當中心像素處于灰度變化緩慢的低頻區域時，其周圍像素容易獲得較大權重，因此中心像素能夠達到良好的平滑效果。而當中心像素位于灰度變化急促的高頻區域時，其周圍像素獲得的權重過小，故不能使中心像素達到有效的平滑效果，導致在高頻區域出現密集性噪聲殘留。也正是得益于此項特性，圖像高頻區域的細節信息才得以較好地保存。而圖6（a）的中值濾波不屬于非局部濾波，不計算搜索窗口中像素的權重，因此不能依據圖像的性質進行選擇性的圖像平滑。

3.3 圖像與算法評價

為了更好地評價圖像重建算法的效能，本文引入了均方根誤差（root mean square error，RMSE）和峰值信噪比（peak signal-to-noise ratio，PSNR）2種圖像的客觀評價指標[20]，其定義為

式（15）中xj和分別代表重建圖像和參考圖像的第j個像素的像素值。RMSE反映的是圖像間的差異程度，是一種基于像素誤差的圖像質量客觀評價指標，用于衡量被處理圖像與理想參考圖像之間的差異。RMSE值越小表示被處理圖像的質量越好。式（16）中max（x）表示圖像x中的最大像素值。PSNR是一種基于誤差敏感的圖像質量評價指標，其值越大表示被處理圖像的質量越好。

表1為3.1與3.2中各實驗結果圖像的RMSE值和PSNR值以及迭代型算法的運行時間（設定總迭代次數k=50）。由表1可知，NLWmedian-CS算法的結果圖像獲得了最小RMSE值和最大PSNR值，表明該算法對圖像的重建效果最佳，這與3.2中的主觀觀察結果一致。算法的運行時間方面，基于非線性壓縮感知的圖像重建算法的時間成本較高，且運行時間與算法所選擇的非線性濾波器種類有關。

表1 結果圖像評價指標與算法運行時間

另外，本文驗證了基于非線性壓縮感知的圖像重建算法的收斂性，記錄了各次迭代結果圖像的RMSE值和PSNR值隨迭代次數的變化情況。圖7為NLW means-CS算法的運行過程中RMSE值的變化軌跡，圖8為該算法的運行過程中PSNR值的變化軌跡。由圖7和圖8可知，算法運行到第10次前后便已經收斂到特定區間。因此，圖示結果支持本文選取各個實驗的第50次迭代的運行結果作為結果圖像。

圖7 RMSE值隨迭代次數的變化

圖8 PSNR值隨迭代次數的變化

4 結語

本文研究了低劑量照射條件下基于非線性壓縮感知的CT影像重建算法的開發問題。對低管電壓條件下所采集的腹部投影數據進行了影像重建，對比了傳統的FBP算法、TV最小化算法和本文提出的基于非線性壓縮感知的算法在腹部影像重建中的效果，把對圖像去噪和保存圖像細節信息具有良好作用的非線性濾波操作引入了壓縮感知的理論框架，有效地強化了壓縮感知在圖像重建領域的作用效果。通過實驗分別驗證了中值濾波、雙邊濾波、非局部均值濾波和非局部中值濾波在壓縮感知算法中的作用效果。實驗結果表明，新算法在CT影像重建任務中比傳統算法在去除圖像噪音、保存圖像細節、刻畫精確邊緣方面具有明顯優勢。然而，重建效果會根據引入的非線性濾波種類的不同而有所差異，主觀觀察和客觀評價指標都表明基于非局部中值濾波的壓縮感知算法獲得了最優的圖像重建效果。新算法的運行需要多個變量參數的共同作用，尤其是引入非局部濾波的情況下，不同的參數組合作用出不同的重建效果，因此參數的設置與訓練是本研究面臨的挑戰。制定能夠依據不同重建任務而進行自適應設置的參數決策系統是下一步研究的方向。針對新方案運行耗時長的問題，課題組將結合頻率濾波來提高算法的運行效率。本文所提出的基于非線性壓縮感知的CT影像重建算法不僅對低劑量照射條件下的圖像重建具有良好的應用效果，而且具有應用在稀疏角度或有限角度掃描的CT圖像重建中的潛在價值。