基于圖像自相似性的多尺度稀疏表示肺4D-CT圖像超分辨率重建

王婷婷，張煜

1.南方醫科大學南方醫院設備器材科，廣東廣州510515；2.南方醫科大學生物醫學工程學院，廣東廣州510515

前言

四維計算機斷層攝影（Four-Dimensional Computed Tomography,4D-CT）能夠真實準確地反映靶區的運動范圍及動態特征，被廣泛地應用于肺癌的精確放射治療之中并取得良好的效果［1］。然而，CT 的高劑量照射特性以及持久的掃描時間，致使數據沿縱向（Z軸方向）無法實現密集采樣，導致數據層間分辨率遠低于層內分辨率，造成各項數據異性，不僅影響圖像可視化而且妨礙數據分析［2-3］。超分辨率（Super Resolution,SR）重建算法是一種提高肺4D-CT 圖像分辨率的有效方法。將數據中不同相位、對應位置的圖像作為多“幀”低分辨率（Low Resolution,LR）圖像輸入，經過SR重建處理，獲得高分辨率（High Resolution,HR）輸出圖像。

本研究提出一種基于圖像自相似性的多尺度稀疏表示SR重建方法，它是一種基于機器學習的重建算法［4-6］，主要解決以下兩大難點：（1）基于稀疏表示的SR重建方法，一般需要相對應的HR和LR圖像塊兩兩組對作為訓練集，以訓練得到反映HR 和LR 圖像之間關系的字典。但是對于肺4D-CT 圖像而言，各項數據異性致使無法獲取冠矢狀面的HR 圖像。在此種情形下，該如何構建訓練集數據？（2）考慮到肺4D-CT 圖像的解剖結構特征存在于不同尺度之下，應該如何選取圖像塊的大小，才能做到同時捕獲不同尺度下的解剖信息？針對第一個問題，筆者引入自相似性的概念。圖像自相似性［7］是源自于對圖像的觀察，發現一些細小的結構在整幅圖像中重復出現。此概念已經被廣泛應用到圖像壓縮［8］、圖像復原［9］、非局部均值去噪［10］以及SR 重建［11］中。我們探索發現，肺4D-CT 數據冠矢狀面圖像和橫斷面圖像不同尺度下的組織結構均存在一定的相似性。因此，筆者利用橫斷面圖像來代替冠矢狀面圖像，構成HR 和LR 圖像塊對，作為訓練集。具體內容將在第1.2 節中詳細闡述。針對第二個問題，引入多尺度策略。圖像多尺度分析最早可以追溯到1980年，源自于高斯和拉普拉斯金字塔的提出［12］。鑒于多尺度分析的優越性，大量算法引入此策略構建不同的多尺度 字 典，典 型 的 包 括Wavelets［13］、Contourlets［14］、Curvelets［15］等。筆者根據四叉樹原理將“根”圖像塊劃分為各個尺度的子圖像塊，并建立全局多尺度字典，詳細內容將在第1.3節中具體介紹。

不同于基于模型的SR 重建算法，本研究算法不僅能避免運動估計，更能夠改善HR圖像細微結構的顯示質量。對比單一尺度稀疏表示的SR算法，本研究的多尺度算法能有效地捕獲肺部不同尺度下的解剖結構特征，重建出細節信息更為豐富的HR 圖像。同時，量化評價指標也體現了本研究算法的優越性。

1 材料與方法

1.1 算法綜述

圖1 基于圖像自相似性的多尺度稀疏表示肺4D-CT圖像SR重建算法示意圖Fig.1 Overview of self-similarity-based multiscale sparse representation for lung 4D-CT image super-resolution reconstruction

第一步：利用橫斷面圖像組成HR 和LR 圖像塊對構成訓練集，訓練得出成對的全局多尺度HR字典DH和LR 字典DL；第二步：輸入的每一個LR 圖像塊PLi，尋求其經過LR字典DL表示后的稀疏系數α，由該系數及HR字典DH重建出此LR圖像塊對應的HR圖像塊PHi。最終，通過拼接所有的HR圖像塊以輸出完整的HR圖像。

1.2 肺4D-CT圖像自相似性

本節中，結合后續依照四叉樹原理所構造的多尺度分析，分別探究了圖像塊尺寸為16×16、8×8和4×4 的冠矢狀面圖像同橫斷面圖像之間的自相似性［16］。利用結構相似性（Structure Similarity,SSIM）這一指標來度量圖像塊之間的相似性，其計算公式可以表示為：

其中，Pt表示橫斷面圖像塊，Pc表示冠矢狀面圖像塊，l(Pt,Pc)、c(Pt,Pc)、s(Pt,Pc)分別為：

其中，μ(Pt)、μ(Pc)分別表示圖像塊Pt和Pc的均值，σ(Pt)和σ(Pc)分別為Pt和Pc方差，σ(Pt,Pc)為Pt和Pc的協方差，C1、C2、C3均為常數。

圖2給出冠矢狀面圖像同橫斷面圖像之間自相似性的驗證結果。其中，圖2a 顯示的是針對肺4D-CT 數據自相似性的量化評價結果。對于大小為16×16的圖像塊來說，74%的冠矢狀面圖像塊能夠在橫斷面上找到至少一個相似塊；而能夠找到50 個相似塊的比例占到了55%。隨著圖像塊尺寸的減小，該比例也在顯著提高。當圖像塊大小為4×4 時，95%的冠矢狀面圖像塊能夠在橫斷面上找到20 個相似塊。圖2b給出冠矢狀面圖像同橫斷面圖像之間自相似性的圖像塊示例。其中，冠狀面與橫斷面的相似結構用紅色方框標出；矢狀面與橫斷面的相似結構用黃色方框標出。

圖2 冠矢狀面圖像同橫斷面圖像之間自相似性的驗證結果Fig.2 Verification result of cross-scale self-similarity between coronal/sagittal images and transverse images

由此，筆者通過實驗從定量和定性兩個方面對冠矢狀面圖像同橫斷面圖像之間的自相似性進行驗證。基于此相似性，面對缺失冠矢狀面HR圖像的情形，筆者采用橫斷面HR和LR圖像塊作為訓練集，來聯合訓練HR和LR字典。

1.3 多尺度稀疏表示SR重建

考慮到肺部圖像解剖結構特征存在于不同尺度之下的數據特性，引入多尺度策略，采用四叉樹劃分原則［17］，構建多尺度圖像塊，利用基于稀疏表示的SR重建算法［18］，以橫斷面HR和LR圖像塊作為訓練集，聯合訓練得出一對全局多尺度字典：HR 字典DH和LR 字典DL，并使得每一對HR 和LR 圖像塊具有相同的稀疏系數，而后對輸入的LR圖像塊尋求其經LR 字典DL表示后的稀疏系數，直接由該系數及HR字典DH來估計得到相應的HR 圖像塊，實現基于圖像自相似性的多尺度稀疏表示肺4D-CT 圖像SR重建。

1.3.1 多尺度模型 如圖3所示，根據四叉樹原理將“根”圖像塊依次向下劃分，得到各個尺度下的子圖像塊。假設“根”圖像塊的大小是n，則依照四叉樹原則其子圖像塊的大小為其中s=0,1,…,S-1表示四叉樹的深度，S表示尺度的個數。

對于尺度s來說，四叉樹上共存在4s個位置。Ds∈Rns×ks表示尺度為s下任意一個位置的字典，該字典由ks個原子構成，每個原子是維數為的列向量。為了能夠同時抓取多尺度數據特性，需要構建一個全局含有所有尺度下各個位置信息的字典D∈Rn×k。因此，對于多尺度字典D而言，它共含有個原子。

圖3 四叉樹分塊方式示意圖Fig.3 Quad-tree model selected for the proposed multiscale framework

由于尺度s的不同，構成單個尺度字典Ds的原子維數ns會發生相應變化。在全局多尺度字典D中如何設定原子的維數呢？筆者設定原子的維數為n，對于維數為ns=的原子，用零將其填滿成n維。圖4給出尺度個數S=3的全局多尺度字典中不同尺度下的原子示例。

圖4 尺度個數S=3的全局多尺度字典中不同尺度下的原子示例Fig.4 Example of multiscale global dictionary atoms with 3 scales

1.3.2 SR重建 通過聯合訓練階段，學習得到成對的全局多尺度HR 字典DH和LR 字典DL。它們可以使得訓練集中成對的HR和LR圖像塊具有相同的稀疏系數，稀疏編碼問題表示如下：

其中，XH={x1,x2,...,}xm表示從訓練集取樣的一系列HR圖像塊，YL={y1,y2,...,}ym表示與之對應的LR圖像塊，XH和YL構成了成對的訓練集圖像塊為l0范數，它表示非零項的個數。Z分別為圖像塊XH和YL的稀疏表示系數矩陣。N和M分別表示HR和LR圖像塊向量形式下的維數?？梢詫⑹剑?）改寫為：

其中，

由于目標函數式（6）關于自變量Dc和Z為非凸函數，無法采用一般優化算法有效求解。但是二者中若一項固定，那么對于另一項來說，上述目標函數是凸函數。因此，一般采用兩步迭代的方法求解。初始化字典Dc之后，循環進行稀疏編碼階段和字典學習階段，直到目標函數收斂或達到預設的迭代次數為止。

（1）稀疏編碼：首先固定全局多尺度字典Dc，通過OMP 算法計算獲得訓練集中“根”圖像塊Xc的稀疏表示系數矩陣Z：

OMP 算法在計算“根”圖像塊Xc的稀疏系數時，是針對“根”圖像塊下的每一個尺度依次計算的，從s=0 到s=S-1，在同一尺度下再逐一計算每一個位置，共計4s位置。

（2）更新字典：與稀疏編碼的方式相同，更新字典中的每一個原子dsl(1≤l≤ks) 也是逐一尺度，逐一位置進行的。在尺度s下，選出所有使用到第l個原子的子圖像塊：

其中，[s,p表]示位于“根”圖像塊Xc尺度為s、位置為p的子圖像塊；z(s,l,p)是對應子圖像塊其系數矩陣的第l行。

對于每一個子圖像塊[s,p]∈wsl，重新計算殘差，此時參與計算的是字典中非dsl的其他原子：

其中，Tsp表示二元矩陣用以從“根”圖像塊中提出子圖像塊[s,p。]以wsl中所有子圖像塊的殘差elsp為列向量，就構成矩陣對Esl進行SVD 分解，即：

用矩陣U的第一列更新字典原子，用矩陣V的第一列乘以Δ(1,1)來更新系數向量z(s,l,p)中集合wsl對應的元素。得到全局多尺度字典Dc即得到成對的多尺度HR 字典DH和LR 字典DL。重建階段，給出LR圖像Y的圖像塊y，需要通過LR字典DL得到它的稀疏系數α：

其中，λ為均衡系數，以權衡稀疏近似誤差和非零項個數的比重。同樣地，采用OMP 算法來求解此稀疏表示最優化問題。

得到LR圖像塊的稀疏系數α之后，HR圖像X的圖像塊x就可表示為x=DHα，從而重建出HR圖像塊，將所有的HR圖像塊拼接后即可輸出最終HR圖像。

2 實驗結果

分別采用仿真數據和真實數據對所提算法進行實驗驗證，同時從視覺和量化兩方面評價實驗結果。

2.1 訓練數據

由于肺4D-CT 數據在軸向是低分辨率的，導致無法獲取冠矢狀面的HR圖像。由此，通過實驗對冠矢狀面圖像同橫斷面圖像之間自相似性進行定性及定量的分析，驗證其有效性?；诖讼嗨菩裕捎脵M斷面HR和LR圖像塊作為訓練集，來聯合訓練HR和LR字典。

按照圖像退化模型式（13）對真實的HR 圖像進行降質處理，得到與之對應的LR圖像：

其中，gk表示第k幅LR 觀測圖像，F表示原始HR 圖像，Mk表示運動變形矩陣，Bk表示模糊矩陣，D表示降采樣矩陣，nk表示加性噪聲。

首先進行降采樣，其系數為2；其次，通過疊加肺部的運動變形場來模擬圖像變形；最后采用高斯函數使圖像模糊。圖5顯示了一組橫斷面的HR 和LR圖像。將所有HR和LR圖像分塊后，每個HR圖像塊和與其對應的LR 圖像塊組對，構成圖像塊對Pair={XH,YL}，以此作為訓練集，為下一步字典訓練奠定基礎。

2.2 多尺度字典

圖5 橫斷面HR和LR圖像舉例Fig.5 Example of high-resolution and low-resolution transverse images

圖6顯示了一組3 個尺度下（S=3）的全局HR 字典DH={Dh0,D1h,D2h}。該全局字典DH由3 個尺度下的字典構成，分別是：圖6a 中s=0 時，字典Dh0；圖6b 中s=1時，字典Dh1；圖6c中s=2時，字典Dh2。實驗中，設定每個尺度下的原子個數均為256。與圖像自相似性驗證實驗中圖像塊的尺寸一致，這里設定“根”圖像塊的大小為16×16，子圖像塊的大小分別是8×8和4×4。因此，字典Dh0中每個原子維數n0=256，字典Dh1中每個原子維數n1=64，字典Dh2中每個原子維數n2=16。特別說明的是，聯合訓練的方式致使與全局HR 字典DH一同得到的還有LR 字典DL，其組成結構與DH相同。

2.3 仿真實驗

本節同樣采用橫斷面圖像作為仿真數據，從視覺和量化評價兩方面對實驗結果進行比較。仿真數據的構造與2.1節中的方式相一致，此處不再贅述。

2.3.1 視覺評價 圖7給出不同算法下橫斷面圖像的重建結果圖。左側第1列顯示的是真實HR圖像。第2、3、4列給出基于單一尺度稀疏表示的SR算法重建結果，其圖像塊大小依次為4×4、8×8、16×16。本章提出基于多尺度的算法結果顯示在第5列，實驗中尺度個數設定為3，“根”圖像塊的大小為16×16，子圖像塊依照四叉樹原則依次為8×8和4×4。各圖中的紅色方框區域圖像均被放大后放置在其下方，以便加強對比。得益于多尺度的分析及處理，本章提出的算法與單一尺度算法相比，可以有效地捕獲肺部不同尺度下的解剖結構特征，從而得到細節信息更為豐富的HR圖像。

圖8給出本章算法與雙線性插值算法及POCS算法的橫斷面重建結果對比圖。同樣左側第1 列顯示的是真實HR圖像。雙線性插值算法和POCS算法的重建結果分別顯示在第2、3 列。第4 列展示的是本章算法結果。將各圖中的紅色方框區域圖像放大后置于第2行便于直觀比較。對比雙線性插值算法，本章算法重建出的HR 圖像清晰度顯著提升；對比POCS算法，本章算法不僅重建得到邊緣有所增強的HR 圖像，更避免了POCS 算法中對結果速度及精度造成限制的圖像配準步驟。

圖6 3個尺度下的全局HR字典DH= {Dh0,D1h,Dh2}Fig.6 Learned 3-scale global high-resolution dictionaryDH= {D0h,D1h,Dh2}

圖7 不同算法下橫斷面圖像的重建結果圖Fig.7 Reconstruction results obtained by different algorithms

圖8 本章算法與雙線性插值算法及POCS算法的橫斷面重建結果對比圖Fig.8 Transverse views of reconstructed images for comparison of the proposed approach with bilinear interpolation and POCS algorithm

2.3.2 量化評價 利用均方根誤差（RMSE）對仿真數據結果進行量化評價。RMSE 值表征兩組數據之間的離散程度，可以衡量真實圖像和重建圖像之間的偏差，其值越小代表兩幅圖像越接近。用公式表示如下：

其中，I表示真實圖像，I表示重建圖像。

根據式（14），分別計算5組數據對不同分塊大小下單一尺度稀疏表示算法、雙線性插值算法、POCS算法和本章算法重建結果的RMSE值，結果如表1所示。雙線性插值算法所得的RMSE 值最高。雖然基于單一尺度稀疏表示的算法分別采用了與多尺度算法中“根”圖像塊及子圖像塊一致的分塊大小，但是單一的處理方式使得大塊中較小結構無法細化，小塊與其所在的局部較大結構也無法保持良好的一致連續性，致使重建效果不佳。而與之相反的多尺度分析，不僅可以把握圖像不同尺度下的結構特征，同時還是一個由大到小、由粗到精的優化處理過程。因此，它所重建出圖像的RMSE值更低，也更為準確。本章算法盡管與POCS 算法重建結果的RMSE 值相差并不顯著，但是本章算法有效避免了運動估計的過程。

2.4 多平面重建結果

表1 不同分塊大小下單一尺度稀疏表示算法、雙線性插值算法、POCS算法和本章算法重建圖像的RMSE值比較結果Tab.1 Comparison of root mean square errors among single-scale sparse representation with different patch sizes,bilinear interpolation,POCS algorithm and the proposed approach

圖9給出單一尺度算法和本章算法下冠矢狀面圖像重建結果。第1、2、3 列是不同分塊大?。?×4、8×8、16×16）下單一尺度稀疏表示算法的重建結果。本章算法的結果顯示在第4 列。為了更加直觀地對比，各圖中的紅色方框區域圖像均被放大后放置在其下方。與基于單一尺度的算法相比，本章算法在捕獲肺4D-CT 數據不同尺度的解剖信息方面更有效，重建出視覺效果更好的HR圖像。

圖10顯示雙線性插值算法、POCS算法和本章算法下冠矢狀面圖像重建結果。第1、2 列分別是雙線性插值算法和POCS 算法的重建結果。第3 列是本章算法的重建結果。各圖中的紅色方框區域圖像均被放大后放置在其下方，以便于比較。對比雙線性插值算法，本章算法重建出的HR圖像清晰度顯著提升；較于POCS算法，本章算法能夠改善HR圖像細微結構的顯示質量。

3 討論

基于模型方法需要采用配準或塊匹配方法進行運動估計，這導致重建速度和精度往往受到配準速度和精度的制約。而基于學習的SR重建方法可以避免這一過程。因此，本研究提出一種基于圖像自相似性的多尺度稀疏表示肺4D-CT圖像超分辨重建方法。通過驗證冠矢狀面圖像和橫斷面圖像中不同尺度下組織結構的自相似性，利用橫斷面圖像來代替冠矢狀面圖像，將對應的HR和LR圖像塊兩兩組對，作為訓練集。其次，面對肺4D-CT 圖像解剖結構特征存在于不同尺度之下的數據特點，引入多尺度策略，根據四叉樹原理將“根”圖像塊劃分為各個尺度的子圖像塊，并建立全局多尺度字典。本研究提出的算法不僅有效地解決了在缺失原始HR 圖像情形下構建訓練集的難題，還通過多尺度分析的手段同時捕獲了肺部不同尺度下的解剖結構特征。

雖然筆者通過實驗從定量和定性兩方面均對冠矢狀面圖像同橫斷面圖像之間的自相似性進行驗證，但是重建結果顯示還是存在一定的干擾，接下來將引入金字塔模型［19］，直接由輸入圖像本身生成訓練集，來對提出的算法進行進一步改善。另外，還將進一步研究如何使該算法能夠用于放大倍數更大的圖像重建，以滿足臨床對于圖像質量及圖像放大的雙重需求。

圖9 單一尺度算法和本章算法下冠矢狀面圖像重建結果Fig.9 Coronal and sagittal views of image reconstructed with single-scale approach and the proposed approach

圖10 雙線性插值算法、POCS算法和本章算法下冠矢狀面圖像重建結果Fig.10 Coronal and sagittal views of images reconstructed with bilinear interpolation,POCS algorithm and the proposed approach

總的來說，本研究提出的算法不僅能夠重建出細節信息更為豐富清晰的HR圖像，還能避免對結果速度及精度產生影響的運動估計，同時也可為4D心臟圖像、4D超聲圖像的SR重建提供新思路。