稀疏掃描型醫用CT的快速成像算法開發

董 建，陳 昊，楊耿煌，丁 洋，楊曉霞，車艷秋

（1.天津職業技術師范大學天津市信息傳感與智能控制重點實驗室，天津 300222；2.天津職業技術師范大學工程實訓中心，天津 300222）

此前的研究已經表明電離輻射存在致癌風險，且計算機斷層掃描（computed tomography，CT）輻射已經被證明與癌癥發病率增加有關[1]。臨床上常通過減少CT掃描時設備旋轉一周的投影視圖數量或降低X射線管功率的方法來減少CT檢查時的輻射暴露[2]。對于低劑量CT的圖像重建可使用投影域預處理、迭代重建算法等方案[3]。對于投影數量不足的圖像重建，常用的濾波反投影重建算法會不可避免地產生條狀偽影[4]。而得益于壓縮感知（compressed sensing，CS）理論的發展，針對稀疏投影或含噪聲投影數據，采用壓縮感知的方法進行圖像重建，逐漸得到廣泛的研究與論證[5]。該類算法可更準確地重建出細節保留完整、偽影分布較少的CT圖像，能夠滿足人們對CT成像質量日漸提高的要求[6-7]。因此，基于壓縮感知理論的CT圖像重建算法開始被廣泛研究和應用[8]。

壓縮感知理論在圖像重建領域中的應用框架是以由數據保真項和正則化項（懲罰因子）共同組成的代價函數為基礎，利用凸函數最優化理論，通過對代價函數進行最小化而得到最優解。正則化項的設計對于求解該類欠定問題具有重要作用。壓縮感知理論在此領域應用時與CT成像的相關性取決于CT圖像是否稀疏，在壓縮感知理論的正則化項設計中常用稀疏變換來實現圖像的稀疏性。2006年Donoho在文獻[5]中指出離散梯度變換和小波變換經常被用于對圖像進行稀疏變換。目前常用的稀疏變換是對圖像的全變分（total variation，TV）變換[9]。基于TV正則項的壓縮感知算法被證明在圖像的稀疏重建中具有良好的效果，且在CT圖像重建領域表現出較強的噪聲抑制特性[10]。文獻[11]提出了全變分最小化的PICS（prior image compressed sensing）算法，利用先驗圖像約束的壓縮感知來提升稀疏投影角度下的圖像重建質量，但此類基于TV正則項的壓縮感知算法在重建后仍存在偽影殘留[12]。

研究發現導致在基于TV正則項的重建圖像中出現潛在模糊結構和補丁狀偽影的主要原因是TV正則項通過線性變換提高預恢復圖像的稀疏性，此過程易過多考慮圖像各像素之間的位置相關性而忽略圖像邊緣、圖像紋理部分等灰度不相關像素間的過渡。為解決上述問題，本團隊在2016年提出了基于非線性濾波的壓縮感知重建算法理論。在進行稀疏變換時，原有的TV線性濾波方法被非線性濾波代替，經過實驗得出基于此法的壓縮感知算法在保留圖像邊緣和細節紋理及消除大塊模糊上具有良好效果[13-14]。

此非線性壓縮感知算法在圖像重建中常進行上千次迭代以得到收斂圖像，且單次迭代運算量大，運算時間長，導致此算法的收斂速度慢。本團隊在2017年提出了基于非線性濾波器的壓縮感知的加速算法，將代價函數拆分成若干個子函數并使用凸函數優化中近端優先理論，每次以行優先準則進行迭代[15]，此算法提高了收斂速度，使算法的可用性大幅提高。

非線性濾波壓縮感知算法在對于恢復圖像交錯處細節和邊緣信息以及抑制大塊模糊方面具有良好效果，但經加速后單次迭代時間仍較長，收斂速度慢。研究發現，基于壓縮感知算法的CT重建技術大多使用空間域濾波方法實現對重建圖像稀疏化的過程。而濾波器的設計對圖像重建質量和重建速度起著關鍵作用[16]。相對的頻域濾波器特別是低通頻域濾波器可去除圖像的高頻噪聲等干擾，且對于平滑圖像具有良好效果[17]。同時低通頻域濾波過程中總體時間復雜度為O（n），遠低于先前研究中所選取的非線性濾波器的時間復雜度O（n2）。在低通濾波器中，巴特沃斯低通濾波器對圖像邊緣的模糊程度大量減小，且沒有振鈴效應，可以在對圖像噪聲及偽影平滑過程中相對較多地保存圖像邊緣細節[18]。指數型低通濾波器具有更快的衰減特性且同樣沒有振鈴效應[17]。

結合以上2種低通頻域濾波器的特點，在原研究的基礎上，用新的低通頻域濾波器替代改進后的壓縮感知重建加速算法中的非線性空間域濾波器進行實驗。本文將通過使用基于低通濾波的壓縮感知重建算法實現在相對較高地保留原可用等級的圖像邊緣和紋理細節的前提下，加快重建圖像收斂速度，使基于此算法的CT圖像重建技術可用性進一步提高。

1 本文算法

1.1 模型的建立

本研究的目的是利用稀疏角度的投影數據重建出臨床可用的CT圖像。此處對圖像重建過程建立模型，將利用投影數據進行圖像重建的問題轉化為對方程組Ax=d的求解問題。其中，向量d=（d1，d2，…，dI）T為掃描得到的投影數據，x=（x1，x2，…，xJ）T為待重建的圖像，A為一個I×J維度的系統矩陣，由CT器械決定。本文研究稀疏角度投影的重建問題，即d的維度I遠小于x的維度J，因此基于該問題的數學模型可轉化為求解欠定方程組的問題[19]。基于壓縮感知理論建立代價函數

式中：L（x）為數據保真項；U（x）為正則化項。

L（x）由每次重建時投影數據的估計值與CT掃描時投影數據真實值的最小二乘差來限定，表示為。系統矩陣A在任務執行中表示對待重建圖像x的取投影操作。正則化項U（x）代表對圖像的稀疏變換。將代價函數進行展開表示

式中：β為用來控制正則化項作用強弱的超參數，亦起到限定運算過程中不發生過擬合的作用；W為對向量x的稀疏化操作，即提取圖像的高頻信息，將x轉換為稀疏向量。

式中：I為單位算子；L操作表示對圖像的空間域平滑濾波。

因此，可將代價函數進一步表達為

1.2 正則化項中的頻域濾波

基于壓縮感知算法的CT重建技術大多使用空間域濾波方法來實現圖像稀疏化的過程，且團隊以往的研究發現基于空間域的NL-median（NonLocal-median）濾波器在重建圖像過程中可以高質量地恢復待重建圖像，但收斂速度過于緩慢[21]。本文的創新點在于將既往研究的空間域濾波升級為低時間復雜度的低通頻域濾波，主要驗證了巴特沃斯低通濾波器和指數型低通濾波器被嵌入正則化項后的工作機制。巴特沃斯低通頻率濾波器和指數型低通頻域濾波器的傳遞函數分別由式（4）、（5）定義，為了更好地描述頻率域處理機制，使用二維矩陣來代表圖像。

式中：u，v為頻率變量；D（u，v）為坐標到中心點的歐式距離；D0為巴特沃斯低通濾波器的濾波半徑；n為濾波器的階數，視實驗條件確定。

由于改用低通濾波器后算法濾波過程是對圖像頻域信息進行操作，在濾波開始前需將由數據保真項限定后的圖像信息通過傅里葉變換從空間域變換為頻率域，離散圖像f（x，y）用式（6）實現傅里葉變換。

式中：M為圖像的行像素數；N為圖像的列像素數。

對于變換后的圖像在頻域中實部和虛部的值分別進行基于式（4）和式（5）的濾波運算，將得到的濾波后的圖像的頻域值進行如式（7）所示的傅里葉逆變換。由此完成正則化中的稀疏變換過程。

綜上所述，基于頻域濾波的壓縮感知在稀疏角度CT圖像重建中的代價函數為

式中：F和F-1分別為圖像的傅里葉正變換與逆變換；H為低通頻域濾波操作。

1.3 近端分裂理論

實施算法推導前，有必要先對其基礎內容的近端分裂理論及其對算法推導的重大貢獻做簡要介紹。首先考慮以下凸優化問題

假定式中f（x）為可能不可微分的下半連續凸函數。近端算子法（proximal operator）為求解此類最優化問題提供了解決方案。對應于f（x）的近端算子表達式為

式中：λ為步長參數。

近端算子具備2種優秀屬性：①其可針對任何一個下半連續凸函數進行定義，即使該函數不可微（如L1范數）；②在步長參數λ＞0的前提下，滿足x=proxλf（g）的定點x與使凸函數f（x）取最小值時的變量x相一致。這2個性質可以應用迭代式求得凸函數f（x）取得最小值時的變量x

式中：k為迭代次數。

該迭代算法叫做近端最小化算法，其為不可微凸函數的最優化問題提供了強有力的解決方案。

凸函數f（x）可拆分成若干分量函數fi（x），即，其中，fi（x）（i=1，2，…，I）均為可能不可微，但滿足下半連續的凸函數。此時，假定凸函數的近端算子求解難度巨大，而其各個分量函數fi（x）的近端算子能夠由低計算成本求得。近端分裂理論能夠為該條件下的凸優化問題提供良好解決方案，其處理方式為在求解由x(k)到x(k+1)的過程中，進一步將問題分解為按照i=1，2，…，I的順序進行各分量函數所對應的近端算子的計算，即

式中：當i=I時，x(k+1)=x(k，i+1)；k為主循環變量；i為子循環變量。步長參數λ(k)值隨k值改變。為了使所求解序列x(k)收斂于函數f（x）取最小值時的變量x，λ(k)需滿足如下條件

近端分裂理論為本文L1范數重建算法的構建提供了理論依據。

1.4 算法實現過程

為了進一步加快算法的運算速度，基于近端算子和近端分裂理論將系統矩陣按行進行拆分。依據近端分裂理論，代價函數式（8）可表示為

式中：g（x）為正則化項；fi（x）為數據保真項。

式中：I為所測得投影數據的維數；αi為對應的系統矩陣A的第i個向量；αiT為向量αi的轉置。

對于數據保真項，由于此算法在迭代時以行更新，可參考ART（algebraic reconstruction technique）算法加速中使用隨機順序行更新法實現對本算法數據保真項的求解加速[22]。對于正則化項，此處引入跨度參數S，使在數據保真項求解的主循環中以S為跨度進行I/S次濾波，S的最優取值視實驗條件確定。數據保真項的行迭代過程如下

式中：proxλ(k)fi（x(k，i)）為近端算子運算；λ(k)為迭代時的步長參數；k為主迭代因子。為了避免數據陷入極限環，令λ(k)=λ(0)/（1+ε·k），λ(0)和ε取值依具體任務確定；i為每次計算數據保真項時的行迭代因子。

通過對保真項求解可得

對正則化項g（x）進行prox算子展開如下

在迭代求解過程中，當λ(k)取值足夠小時，x無限接近已知量x(k，i+1)，因此在λ(k)取值足夠小的前提下，H·F（x）可由當前圖像的濾波結果H·F（x(k，i+1)）代替。因此，式（18）的求解過程可由式（19）的固定形式的軟閾值更新實現。

至此完成了本文中提出的重建算法的構建。基于頻域濾波的壓縮感知應用于稀疏掃描CT圖像重建的算法總結如下：

輸 入 像素數為N×M的投影數據d。

初始化 設置步長參數的控制變量λ(0)，ε，設置β＞0的超參數以及S＞0的跨度參數。

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

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

步驟3 圖像更新程序

傅里葉正變換，求得F（x(k，i+1)）；低通頻域濾波處理，求得H·F（x(k，i+1)）；傅里葉逆變換，求得F-1（H·F（x(k，i+1)）），并賦值x(k，i+1)←F-1（H·F（x(k，i+1)））。

for（j=0；j＜N*N；j++）do

式（19）軟閾值處理。

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

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

實驗中發現將非線性濾波器用低通頻域率濾波器替代后，重建圖像的收斂速度顯著提升，但重建圖像中有放射狀偽影未被去除。由此根據此前重建算法的研究，在上述算法的基礎上，加入微弱的TV正則化項用于去除重建圖像中的放射狀偽影。代價函數式升級為式（20）。諸如式（20）雙正則化項的代價函數的最優化及算法推導與前述推導方法類似，在文獻[21]中有詳細論述。

式中：γ為用來控制TV項作用強弱的超參數；‖x‖TV為對圖像進行TV最小化。

本研究中限定取值足夠小，即讓TV項起到微小的輔助作用，旨在有效利用TV項在圖像處理中對噪聲的抑制作用。被定義為TV范數，其可看作是一種特殊的L1范數，即灰度梯度的L1范數，用于去除重建圖像中的放射狀偽影。

2 實驗結果與分析

在本文實驗設計中選用口腔CT平掃影像，該影像牙齒、牙槽骨面積大且與其他組織的射線吸收系數差異較大，在圖像重建過程中易在高亮部位產生放射狀偽影，此外牙齒形狀及排列方式的細節豐富，因此可采用此影像來對比新舊算法重建效果。本實驗采用圖像如圖1所示，像素為512×512，圖像CT值（hounsfield unit，HU）最小處為-1 000，最大處3 071，由此可較準確地反映圖像中不同部位的對比度差異。

圖1 人體口腔CT平掃圖像

本研究分別對圖1從180°中以間隔為2取90個角度的投影數據，間隔為3取60個角度的投影數據以及間隔為4取45個角度的投影數據模擬真實情況下的稀疏角度的CT掃描過程。以模擬投影數據為輸入進行基于NL-median和巴特沃斯低通頻域濾波器、指數型低通頻域濾波器壓縮感知圖像重建對比實驗。

2.1 重建結果對比展示

針對新算法中使用的巴特沃斯低通濾波器和指數型低通濾波器的參數D0和n的選擇進行如下實驗：以90個角度的投影數據以及巴特沃斯低通濾波器為例，固定參數n=10，分別設置D0=1、D0=400與D0=500進行重建實驗，對比發現D0的最佳取值位于1～500，經過實驗得到其最佳取值為400。固定D0=400，分別設置n=1、n=10與n=100進行重建實驗，對比發現n的最佳取值位于1～100，經實驗得到其最佳取值為10。90個角度投影數據不同頻域濾波參數重建對比圖像如圖2所示，基于45、60、90個投影數量的不同正則化項的壓縮感知重建算法的重建圖像結果如圖3所示。

圖2 90個角度投影數據不同頻域濾波參數重建對比圖像

圖3 基于不用濾波器不同投影數量的壓縮感知重建圖像

圖3（a）、（b）、（c）從左至右依次是依據45個、60個、90個投影角度正弦圖的重建結果。

在本實驗的重建過程中，對于不同的正則化項選取同樣的共同參數：λ(0)=0.2，ε=1，β=1 200，S=4。圖3所展示的重建圖像是分別迭代100次后得到的重建結果。在圖3（a）的實驗中，對于NL-median濾波器設置的搜索窗口和濾波窗口分別為7×7和5×5，設置δ1=10，δ2=10 000。基于此，重建的圖像邊緣比較清晰，細節恢復比較完整。在圖3（b）的實驗中，對于巴特沃斯濾波器設定濾波半徑D0=400，傳遞函數中參數n=10，同時在迭代100次的前60次為了消除重建圖像的放射狀偽影加入γ=2的TV正則化項。基于此，重建的圖像同樣可以觀測到，重建圖像邊緣較清晰，細節較完整。在圖3（c）的實驗中，將指數型低通濾波器濾波半徑設為D0=400，傳遞函數中參數n=10，在迭代100次的前60次為了消除重建圖像的放射狀偽影加入γ=2的TV正則化項。基于此，重建的圖像可以觀測到同上述巴特沃斯頻域濾波器有類似的效果。

2.2 重建圖像及算法的客觀評價

為了更好地評價基于不同濾波器的CS算法在重建后得到的收斂圖像的質量優劣，引入了平均絕對誤差（mean absolute deviation，MAD）作為評價指標。MAD的運算公式為

式中：J為圖像像素個數，在本實驗中取值512×512；xj為每次運算后得到的圖像的各個像素的CT值；Xj為原始圖像，即圖1中各像素的CT值。在此評價體系中MAD的值越小，最終的收斂圖像相較于原圖的重建質量越好，同時用最終收斂時間來評價各個算法執行圖像重建任務的可用性。上述各算法分別以45個、60個、90個投影角度迭代100次后的收斂圖像的各評價指標數值如表1所示。

表1 以45、60、90個投影角度重建圖像的評價指標

為了更精確地比較重建圖像與原始圖像的CT值，重建圖像與原圖第256行像素CT值的擬合曲線如圖4所示。

圖4 重建圖像與原圖第256行像素CT值的擬合曲線

圖4（a）、（b）、（c）分別給出了基于45個、60個、90個角度的投影數據重建后的第256行的各個像素的CT值與原始圖像CT值的擬合情況。從圖4（a）、（b）可知，在以45個、60個角度的投影數據重建后，基于巴特沃斯和指數型濾波的CS重建算法與NLmedian-CS在對圖像恢復上均在圖像邊緣跳變處有沖激，在圖像過渡平緩的區域擬合效果良好。從圖4（b）、（c）可知，紅色標注處所示的當投影數據達到60個、90個角度時，基于巴特沃斯和指數型濾波的CS重建算法的圖像在圖片CT值突變沖激處比原圖像的擬合效果有很大改善。

3 結語

本研究為了進一步提高基于稀疏投影數據的圖像重建算法的可用性，對之前提出的NL-median-CS算法中的正則化項在稀疏化時的濾波器進行了重新選擇，基于新型頻域濾波器的壓縮感知算法在稀疏投影重建時亦可保持較高圖像質量。由實驗結果可知，圖像細節紋理、邊緣信息可相對完整地保留，且基于新濾波器的CS算法在迭代速度上有了顯著提高，大幅降低了收斂時間。實驗結果顯示，相同迭代次數下基于新濾波器的CS算法的迭代時間約是先前的基于NL-median濾波器的1/10，顯著增強了該算法的可用性。在本文提出的基于新濾波器的CS算法中仍有較多參數需要最優化，且為了獲得較好的重建效果常需要繁雜的參數訓練過程，該過程不利于算法魯棒性的提升。同時，從整體上看，非線性壓縮感知重建算法在重建中獲得收斂圖像的速度仍與正在CT中使用的FBP算法有較大差距。使算法中各參數可根據不同部位的CT圖像自適應改變是本團隊下一步研究的主要方向。本文所提出新的基于頻域濾波的壓縮感知CT圖像重建算法顯著提高了圖像的重建速度，且較好地保證了圖像質量，具有一定的實際應用價值。