基于低維流形先驗的低劑量CT重建方法

牛善洲，梁禮境，李 碩，張夢真，邱 洋，劉漢明，李 楠

1.贛南師范大學數學與計算機科學學院，江西贛州341000

2.贛南師范大學贛州市計算成像重點實驗室，江西贛州341000

3.贛南師范大學經濟管理學院，江西贛州341000

X 射線計算機斷層（computed tomography，CT）成像在疾病防治與診斷方面取得了巨大成就，其作為一種無痛、無創傷性的檢查方法，可用于全身多臟器檢查，是現代影像學的杰出代表[1]。CT 圖像的質量與X 射線的輻射劑量緊密相關，劑量越高圖像質量越好，但是過高劑量的X 射線照射會造成患者脫發、皮膚灼傷，甚至可能誘發癌癥或者遺傳性疾病[2-4]。針對CT檢查過程中的X 射線照射劑量問題，世界衛生組織、國際放射委員會以及國際醫學物理組織制定了X 射線照射劑量保證和劑量控制標準，極力主張X 線CT 檢查應遵循實踐正當性、防護最優化的原則，希望以最小的代價和劑量獲取最好的CT 影像診斷效果。因此，在保證圖像質量的前提下最大限度地減少X射線的輻射劑量已成為CT成像領域亟待解決的關鍵問題。

降低管電流或者管電壓是實現低劑量CT成像最簡單有效的方法，但是投影數據會受到量子噪聲和電子噪聲的污染，從而使得FBP 算法重建圖像的質量嚴重退化，無法滿足臨床診斷的需求[5-6]。低劑量CT 圖像重建方法主要分為基于圖像域重建和基于投影域重建兩大類。基于圖像域的重建方法是在圖像域中對待重建的圖像進行系統建模，然后通過迭代算法求解目標函數進行圖像重建[7-8]。但是，由于迭代重建需反復進行投影與反投影運算，導致運算量大，重建速度慢，無法滿足快速診斷的需要。基于投影域的方法利用投影數據的噪聲統計特性建立投影數據恢復模型，得到高質量的投影數據后利用FBP 算法重建最終CT 圖像[9-15]。該類方法與圖像域迭代重建方法相比計算時間大幅減少，在去除噪聲的同時可以很好地保持圖像的空間分辨率。Lu 等[16]通過大量的臨床低劑量CT 實驗以及相應的理論分析，提出經對數變換后的投影數據噪聲近似滿足均值和方差呈非線性關系的高斯分布。基于投影數據這一重要統計特性，Wang等[12-13，17]提出多種基于吉波斯隨機場先驗的懲罰加權最小二乘（penalized weighted least-squares，PWLS）低劑量CT 投影數據恢復方法，如基于KL 變換的懲罰加權最小二乘方法[13]、多尺度加權懲罰加權最小二乘方法[12]等。但Wang等提出的基于吉波斯隨機場先驗的懲罰加權最小二乘方法的目標函數使用的是二次懲罰函數，而在投影數據恢復過程中，二次懲罰函數對投影數據的均勻區域和邊緣區域使用相同的懲罰權重會造成圖像的邊緣和結構細節信息的丟失，從而會導致重建圖像的空間分辨率下降。

圖像塊流形的一個重要特性是其具有低維結構，并且包含豐富的結構信息。圖像塊可以視為嵌在高維空間中低維流形的點云，極小化圖像塊流形的維數可以作為先驗信息進行圖像恢復。低維流形先驗將線性對象的低秩先驗擴展到具有更復雜結構的數據對象上，其可以更好地恢復圖像的結構細節信息。因此，將低維流形的維數作為先驗信息引入到低劑量CT 投影數據恢復過程中，提出了一個基于低維流形懲罰加權最小二乘（penalized weighted least-squares based on low-dimensional manifold，PWLS-LDMM）的低劑量CT 投影數據恢復模型，通過極小化相應的目標函數可以得到理想的投影數據，再使用FBP 算法對恢復后的投影數據進行CT 重建。實驗結果表明該方法在有效抑制低劑量CT圖像噪聲和偽影的同時可以很好地保持圖像的邊緣和結構細節特征。

1 模型與方法

1.1 PWLS重建模型

先前的研究工作表明[18]，在低劑量CT 掃描協議下經過系統校準和對數轉換的投影數據近似遵循高斯分布，數據樣本均值和方差之間滿足如下的公式：

其中，I0是入射X線強度，是第i個探測器單元上投影數據（對數變換后）的均值，是電子噪聲的方差。設y=(y1,y2,…,yM)T為系統校正和對數變換后的投影數據測量值，q=(q1,q2,…,qM)T為待估計投影數據的理想值，則基于PWLS的低劑量CT重建模型為：

1.2 低維流形先驗

設投影數據q為m×n的矩陣，選取q中的任意一個像素點(i,j)，在像素點(i,j)處選取大小為s1×s2的圖像塊Ψx(q)，其中，s1和s2分別為圖像塊Ψx(q)的行數和列數，為s1×s2矩陣左上角的像素。對q中所有像素點的圖像塊Ψx(q)組合成q的點云Ψ(q)：

其中，d=s1×s2,Rd為d維的歐幾里德空間。

基于微分流形理論，點云Ψ(q)可以模擬為一個嵌在Rd空間中的低維光滑流形M=∪l Ml，其中Ml是對應q中不同區域的流形，M為q的圖像塊流形。投影數據q的低維流形先驗可以寫為[19]：

其中，dim(M)為流形M的維數，αi(w)=wi為指標函數，且w=(w1,w2,…,wd)∈M,?Mαi(w)為指標函數αi(w)在流形M上的梯度且為歐幾里德范數。

1.3 PWLS-LDMM重建模型

PWLS-LDMM的投影數據恢復模型可以寫為：

其中，p和q都為待恢復投影數據，T 為轉置運算，Σ為對角矩陣，對角矩陣對角線上的元素為方差β2>0 為正則化參數為q的低維流形先驗。

使用如下的交替優化算法求解PWLS-LDMM投影數據恢復模型：

其中，k為迭代次數。

由于式（6）的目標函數為光滑的二次凸函數，可以對其求導并令導數為零來進行求解。得到其解為：

其中，I為單位矩陣，Σ-1為對角矩陣Σ的逆矩陣，β1為正則化參數，(Σ-1+β1I)為對角矩陣。

使用交替優化法求解式（7），先固定流形更新投影數據，再固定投影數據更新流形。具體求解步驟如下：

（1）先固定流形初始值M(k)，通過求解下式得到指標函數α(k+1)i(i=1,2,…,d)和投影數據q(k+1)：

其中，Ψ*為Ψ的伴隨算子，Ψ*Ψ為對角矩陣。

（2）通過指標函數對流形進行更新，公式如下：交替求解公式（6）、（7）得到恢復后的投影數據后，再使用FBP算法，即可重建出最終CT圖像。

綜上所述，基于低維流形先驗加權最小二乘的低劑量CT重建方法的計算步驟可以寫成如下所示：

步驟1 初始化令k=0；

步驟2 將qk代入式（6），通過式（8）得到pk+1；

步驟3 將pk+1代入式（7），通過式（9）~（13）得到αk+1和qk+1，并通過式（14）對流形進行更新；

步驟4 重復步驟2~步驟3，直至滿足迭代終止條件；

步驟5 利用FBP算法將步驟3得到的投影數據qk+1重建出CT圖像。

1.4 對比方法

為了驗證PWLS-LDMM 方法的有效性，將其與傳統的FBP方法，基于二次模先驗的懲罰加權最小二乘方法[14]（penalized weighted least-squares via quadratic membrane，PWLS-QM）和基于字典學習先驗的懲罰加權最小二乘（penalized weighted least-squares via dictionary learning，PWLS-DL）方法[20]進行比較。

PWLS-QM方法的目標函數為：

其中，Ni是投影數據中第i個像素的一階鄰域，wim在水平方向取值為1，在垂直方向取值為0.25。

PWLS-DL方法的目標函數為：

其中，D∈RL×K為訓練字典，αs∈RL×1是含有少量非零元的向量，Es∈RL×K是從圖像f提取的塊矩陣，vs是拉格朗日乘子，γ是超參數。

1.5 實驗設置

分別使用Shepp-Logan 體膜（如圖1（a）所示）和XCAT 體膜（如圖1（b）所示）以及臨床數據實驗來對所提出的PWLS-LDMM方法進行驗證，并使用文獻[21]中的方法仿真生成低劑量CT 投影數據。數值體膜實驗中，CT成像幾何采用扇形束和弧形探測器，X射線源到旋轉中心和探測器的距離分別為570 mm 和1 040 mm，旋轉角在[]0,2π 間采樣值為1 160，每個采樣角對應672個探測器單元，探測器單元的大小為1.407 mm。圖像維數為512×512，像素大小為1 mm×1 mm。Shepp-Logan和XCAT 體膜室驗的入射光子強度I0分別為5.0×104、3.0×104，電子噪聲的方差σ2e都為10.0。在臨床數據實驗中的腹部投影數據，由商業單排CT 設備采集。X 射線源到旋轉中心和探測器的距離分別為570 mm 和1 040 mm，旋轉角在間采樣值為1 160，每個采樣角對應672個探測器單元，探測器單元的大小為1.407 mm。在球管電壓為120 kV 下，分別獲取400 mAs（高劑量）和50 mAs（低劑量）的CT 投影數據。重建圖像的維數為512×512，像素大小為0.6 mm×0.6 mm。分別使用相對均方根誤差（relative root mean square error，RRMSE），結構相似性指標[22]（structural similarity，SSIM）和特征相似性指標[23]（feature similarity，FSIM）對重建圖像進行定量分析。使用MatlabR2020b（Ubuntu）編程，所有程序都在Intel?Reon?Gold 6234 2.33 GHz，16 核處理器，64 GB內存PC機上實現。

圖1 數值體膜圖像Fig.1 Digital phantom

2 實驗結果與分析

2.1 Shepp-Logan體膜實驗

圖2為不同方法重建的Shepp-Logan 體膜圖像。圖2（a）為FBP算法重建的圖像；圖2（b）為PWLS-QM方法重建的圖像；圖2（c）為PWLS-DL 方法重建的圖像；圖2（d）為PWLS-LDMM方法重建的圖像。從圖2可以看出，FBP 方法重建的圖像質量嚴重退化，存在大量的噪聲和條形偽影。PWLS-GS重建的圖像相比于FBP重建圖像、噪聲和偽影得到抑制，但仍含有一些噪聲和條形偽影。PWLS-DL方法重建的圖像噪聲得到了一定程度上的抑制，但仍含有條形偽影。PWLS-LDMM 方法重建的圖像中噪聲和條形偽影都被有效抑制，同時可以保持邊緣和細節結構特征。因此，本文提出的PWLSLDMM方法在抑制噪聲偽影和保持圖像分辨率方面都有良好表現，可以保持良好的圖像一致性。

圖2 不同方法重建的Shepp-Logan體膜圖像Fig.2 Shepp-Logan images reconstructed by different methods

表1為Shepp-Logan 體膜重建圖像的結構相似性、特征相似性和相對均方根誤差。與FPB、PWLS和PWLSDL 方法相比，PWLS-LDMM 方法的相對均方根誤差分別降低了64.87%、54.81%和7.02%。與FBP、PWLS-QM和PWLS-DL 方法相比，PWLS-LDMM 方法的SSIM 值分別提高了16.78%、1.88%和1.91%；FSIM 值分別提高了12.68%、2.73%和1.65%。PWLS-LDMM 方法重建的圖像具有最高的SSIM 值和FSIM 值，即最接近于真實的體膜圖像。

表1 Shepp-Logan體膜圖像的RRMSE、SSIM和FSIMTable 1 RRMSE，SSIM and FSIM of Shepp-Logan

為了更清晰地比較各種方法對噪聲和偽影的抑制效果，圖3給出了圖1中感興趣區域（如圖1（a）中紅色方框所示）的局部放大圖。與FBP、PWLS-QM 和PWLSDL 三種方法相比，PWLS-LDMM 方法能在有效去除噪聲和偽影的同時可以保持圖像的邊緣信息和結構細節特征。進一步，使用SSIM 和FSIM 值對局部放大圖進行定量分析。圖4 給出了圖3 中各重建結果的SSIM 值和FSIM 值。與FBP、PWLS-QM、PWLS-DL 方法相比，PWLS-LDMM 方法的SSIM 值分別提高了329.71%、11.74%和23.41%；FSIM 值分別提高了111.84%、8.31%和13.48%。PWLS-LDMM方法的SSIM和FSIM值均高于其他方法的SSIM 值和FSIM 值，即PWLS-LDMM 方法重建的圖像最接近真實體膜。

圖3 Shepp-Logan體膜的局部放大圖Fig.3 Zoomed-in views of Shepp-Logan phantom

圖4 圖3中圖像的SSIM值和FSIM值Fig.4 SSIM and FSIM values of images in Figure 3

2.2 XCAT體膜實驗

圖5為不同方法重建的XCAT 體膜圖像。圖5（a）為FBP方法重建的圖像；圖5（b）為PWLS-QM方法重建的圖像；圖5（c）為PWLS-DL 方法重建的圖像；圖5（d）為PWLS-LDMM 方法重建的圖像。從圖5 可以看出，FBP 方法重建的圖像存在較多的噪聲和條形偽影；PWLS-GS 重建的圖像中仍含有一些噪聲和條形偽影；PWLS-DL方法重建的圖像噪聲得到了一定程度上的抑制，但仍含有條形偽影；PWLS-LDMM方法重建的圖像中噪聲和條形偽影都被有效抑制。因此，本文提出的PWLS-LDMM方法在抑制噪聲偽影和保持圖像分辨率方面都有良好表現，對比其他三種方法得到了更優質的重建圖像。

圖5 不同方法重建的XCAT體膜圖像Fig.5 XCAT images reconstructed by different methods

表2為XCAT 體膜重建圖像的結構相似性、特征相似性和相對均方根誤差。如表2所示，PWLS-LDMM方法與其他方法相比，在降低相對均方根誤差方面有很好的表現。與FBP、PWLS-QM 和PWLS-DL 方法相比，PWLS-LDMM方法的相對均方根誤差分別降低了37.46%、22.17%和11.48%。與FBP、PWLS-QM 和PWLS-DL 方法相比，PWLS-LDMM方法的SSIM值分別提高了3.33%、0.73%和1.22%；FSIM 值分別提高了0.83%、0.68%和0.71%。PWLS-LDMM 方法的SSIM 和FSIM 值均高于其他方法的SSIM 值和FSIM 值，即PWLS-LDMM 方法重建的圖像最接近真實體膜。

表2 XCAT體膜圖像的RRMSE、SSIM和FSIM Table 2 RRMSE，SSIM and FSIM of XCAT images

同樣，為了更清晰地比較各種方法對噪聲和偽影的抑制效果，圖6給出了圖1中感興趣區域（如圖1（b）中大小為60×100 的紅色方框所示）的局部放大圖，并使用SSIM 和FSIM 值來對結果的感興趣區域進行定量分析。圖7 給出了圖6 中各重建結果的SSIM 值和FSIM值。與FBP、PWLS-QM 和PWLS-DL 方法相比，PWLSLDMM 方法的SSIM 值分別提高了80.17%、8.97%和20.46%；FSIM 值分別提高了26.17%、4.98%和8.30%。PWLS-LDMM方法的SSIM和FSIM值均高于其他方法的SSIM 值和FSIM 值，即PWLS-LDMM 方法重建的圖像最接近真實體膜。

圖6 XCAT體膜圖像的局部放大Fig.6 Zoomed-in views of XCAT phantom

圖7 圖6中圖像的SSIM值和FSIM值Fig.7 SSIM and FSIM values of images in Figure 6

2.3 臨床數據實驗

圖8為臨床腹部CT投影數據重建圖像。圖8（a）為400 mAs條件下FBP算法重建得到的高劑量圖像，其作為金標準圖像來對實驗結果進行定量分析；圖8（b）為50 mAs 條件下FBP 方法重建的圖像；圖8（c）為50 mAs條件下PWLS-QM 方法重建的圖像；圖8（d）為50 mAs條件下PWLS-DL方法重建的圖像；圖8（e）為50 mAs條件下PWLS-LDMM 方法重建的圖像。從圖中可以看出，本文提出方法在噪聲濾除和邊緣細節以及圖像分辨率保持方面表現更為突出，優于對比方法。

圖8 臨床數據重建圖像Fig.8 Images reconstructed by clinical data

表3為低劑量臨床數據重建的CT圖像的結構相似性、特征相似性和相對均方根誤差。如表3所示，PWLSLDMM 方法與其他方法相比，在降低相對均方根誤差方面有很好的表現。與FBP、PWLS-QM 和PWLS-DL方法相比，PWLS-LDMM 方法的相對均方根誤差分別降低了45.96%、25.61%和15.87%。與FBP、PWLS-QM和PWLS-DL 方法相比，PWLS-LDMM 方法的SSIM 值分別提高了19.12%、7.46%和8.63%；FSIM 值分別提高了4.37%、1.83%和1.14%。PWLS-LDMM 方法的SSIM和FSIM 值均高于其他方法的SSIM 值和FSIM 值，即PWLS-LDMM方法重建的圖像最接近高劑量圖像。

表3 低劑量臨床數據重建的CT圖像的RRMSE、SSIM和FSIMTable 3 RRMSE，SSIM and FSIM of low-dose CT images reconstructed by clinical data

進一步，為了更清晰地比較各種方法對噪聲和偽影的抑制效果，圖9給出了重建圖像感興趣區域（如圖8（a）中紅色方框所示）的局部放大圖。圖10給出了圖9中各重建結果的SSIM值和FSIM值。與FBP、PWLS-QM和PWLS-DL 方法相比，PWLS-LDMM 方法的SSIM 值分別提高了46.88%、1.02%和26.71%；FSIM值分別提高了8.34%、2.59%和3.78%。

圖9 臨床數據重建CT圖像的局部放大圖Fig.9 Zoomed-in views of CT images reconstructed by clinical data

圖10 圖9中圖像的SSIM值和FSIM值Fig.10 SSIM and FSIM values of images in Figure 9

3 結論

為了提高低劑量CT 重建圖像的質量，將低維流形作為先驗信息引入投影數據的恢復過程中，提出了一種基于低維流形先驗的低劑量CT重建方法。數值仿真和臨床數據實驗結果表明，本文提出的低劑量CT 投影數據恢復模型在噪聲和條形偽影抑制方面都取得了很好效果。為了驗證PWLS-LDMM 方法的有效性，使用多種定量指標對重建結果進行分析與比較。與FBP、PWLS-QM 和PWLS-DL 對比方法相比，PWLS-LDMM方法在相對均方根誤差、結構相似性和特征相似性等方面均有上佳的表現。PWLS-LDMM方法在噪聲去除，條形偽影抑制已經結構細節保持方面明顯優于對比方法。