提出三維凸集投影(3D POCS)算法,用于實現三維肺部CT圖像的超分辨率重建;并采用多分辨率混合顯示方式實現肺結節的三維可視化。首先,構建多個有亞像素級位移的低分辨率三維圖像,并生成參考圖像;其次,利用三維運動估計方法,將低分辨率圖像映射到高分辨率參考圖像上;利用一致性約束凸集對三維參考圖像進行修正,迭代重建出高分辨率三維圖像;最后,混合顯示不同分辨率下繪制的圖像。實驗選取5組圖集做性能評價,并與3種插值方法進行比較。主客觀兩方面的評價顯示,3D POCS算法實現三維圖像的超分辨率重建性能優于其他方法;混合顯示方式能夠滿足肺結節的高分辨率三維可視化的需要。
引用本文： 王兵, 樊星, 楊穎, 田學東, 顧力栩. CT圖像肺結節的三維超分辨率重建及顯示. 生物醫學工程學雜志, 2015, 32(4): 788-794. doi: 10.7507/1001-5515.20150143 復制
引言
根據世界衛生組織2014年5月發布的最新資料顯示，肺癌是全球前十位死亡原因的第五位，癌癥死亡的首位[1]。目前，醫學上對肺癌的早期檢測主要利用計算機斷層掃描技術(computer tomography,CT)。國內外研究表明，計算機輔助檢測/診斷(computer assisted detection/diagnosis，CAD)系統對肺部CT圖像進行自動分析，不僅減少了醫生的工作量，而且提高了檢測和診斷的準確性，為肺癌的早期檢測和診斷提供了有力的支持[2]。在CAD系統中，像病變部位及周圍的器官組織等醫生感興趣區域(regions of interest，ROI)的三維逼真可視化具有重要的意義，醫生可以交互地從多個角度觀察ROI大小、形態特征，獲得病變部位更為準確和客觀的信息以做出更為準確的診斷[3]。然而，由于CT固有的高劑量照射，沿Z軸方向的密集釆樣往往是不實際的，導致層間分辨率遠低于層內分辨率，造成數據顯著的各項異性。另外，由于結節體積相對較小(小結節通常直徑不超過5 mm)，對應的切片數量少，直接進行三維繪制難以達到醫學診斷所需的真實感效果。因此，需要三維圖像沿Z軸增加層與層之間的信息量，以獲得肺結節高清晰度可視化效果。采用插值方法增加層間切片數量，如線性插值、最近鄰插值和立體卷積插值等，雖然獲得的三維圖像體素數有所提高,但圖像所包含的頻率信息并沒有什么變化,因此圖像并沒有增加細節信息，三維繪制之后物體表面仍會存在明顯的梯田狀粗紋理效果[4]。增加圖像細節信息的一種有效方法是超分辨率重建技術，即由同一場景的多幀具有互補信息的低分辨率圖像序列，按照相應的算法獲得一幅或多幅清晰的高分辨率圖像的過程。Harris和Goodman于20世紀60年代初提出了針對單張圖像重建的超分辨率重建技術的概念和方法，之后許多研究者在此基礎上相繼提出橢圓球波函數法、線性外推法以及疊加正弦模板法等；1984年Tsai和Huang引入多幀增強的概念，并且首次在頻域中提出基于序列圖像的超分辨率重建的概念和方法[5]。利用序列和多幅圖像進行超分辨率重建成為研究熱點。20世紀90年代以來，各種空間域算法的提出，使得超分辨率重建技術得到了更大的發展。空間域算法主要有迭代反投影方法[6]、凸集投影法(projection onto convex sets,POCS)[7]、統計復原方法[8](最大后驗概念估計器和最大似然估計器)等。基于學習的超分辨率重建方法從大量的訓練樣本集中獲取先驗知識作為超分辨率的依據，在不增加輸入圖像樣本數量的前提下，產生新的高頻細節，提高重建圖像的質量。由于層析成像技術的特殊機制，使得超分辨率重建技術在醫學圖像處理領域獲得廣泛的研究[9]。基于集合理論重建的理論體系是超分辨率重建領域中的一種主流方法，其實現簡單，能方便地加入先驗信息，可很好地保持高分辨率圖像上的邊緣和細節，適用于醫學圖像中局部精確信息的復原或重建[7]。然而傳統的POCS算法是針對二維圖像的超分辨率重建，只考慮圖像內部信息而忽略了圖像之間的關聯信息。對于醫學斷層掃描圖像，由于病灶特征信息不僅存在于各個切片中(層內)，也存在于多個連續的切片之間(層間)，因此，需要一種基于更高維的超分辨率重建方法。本文將二維圖像超分辨率重建的POCS方法擴展至三維數據場，提出三維POCS方法，即利用多個有亞像素級位移的低分辨率三維圖像間的信息重建出高分辨率的三維圖像，用于三維肺部CT圖像的超分辨率重建；并采用多分辨率混合顯示方式實現肺結節的高清晰度三維可視化。
1 三維凸集投影方法
POCS方法是使用一系列凸約束集合來描述超分辨率重建方法中的一些特性，如正定性、能量有界性、平滑性等，這些約束集對高分辨率圖像的解進行約束，用于保證與觀測數據(低分辨率圖像)的一致性，通常稱之為一致性約束凸集，超分辨率圖像解空間與這一系列的凸約束集合相交。POCS方法是一個迭代過程，預先估計一幀高分辨率圖像作為成像空間中的初始參考圖像，然后從參考圖像的某一點開始，將圖像當前估計值投影到一致性約束凸集上，判斷是否滿足所有的凸形約束，如果不滿足，則將其殘差反投影到高分辨率估計圖像上對其進行修正，經過有限次的迭代，直到找到滿足所有凸約束集的解，重建出高分辨率圖像。因此，POCS重建算法主要分為三個具體步驟：構造參考圖像、進行運動估計、對參考圖像進行修正。在一個病例僅有一個三維肺部CT圖像的情況下，重建之前需要構造有亞像素級位移的多個低分辨率三維圖像，模擬多個相關的、有運動位移的、反映同一場景的低分辨率三維圖像，由多個這樣的三維圖像重建一個高分辨率的三維圖像。
1.1 構造有亞像素級位移的低分辨率三維圖像
在超分辨率重建中，普遍接受的退化成像過程模型為“變形-模糊”過程，結合肺部橫截面CT圖像特點，假定圖像之間的相對運動(圖像變形)采用全局平移運動模擬。設理想高分辨率三維圖像為g(x,y,z)，其空間分辨率為qN×qN×bP，其中qN×qN為層內像素數，bP為層數，q和b為降采樣因子或放大倍數。對于三維肺部CT圖像，通常bP>qN。為了盡可能獲得正確比例的三維圖像，圖像層內和層間使用不同的降采樣間隔。為了模擬多個相關的、有運動位移的、反映同一場景的低分辨率三維圖像，需要對g(x,y,z)進行多次降采樣，并且兩次降采樣之間移動小于一個像素的距離，稱為亞像素級位移，各圖像之間的相對位移保持一致。
圖 1為4個互有1/2像素位移的降采樣過程示意圖，移動路線分別為不移動、X軸正向移動、Y軸正向移動、X軸負向移動。虛線為高分辨率數據，實線為降采樣數據，降采樣因子q=2，b=1，得到分辨率為N×N×bP的多個低分辨率圖像。即，通過每次整體錯開一個高分辨率像素大小的位移，經過降采樣后獲取的低分辨率圖像之間的相對位移即為亞像素級位移，從而獲得有亞像素級位移的一組低分辨率數據。
圖1
4個有亞像素級位移的三維低分辨率數據的構建過程
Figure1.
Construction of 4 sub-pixel 3D low-resolution data
1.2 構造參考圖像
POCS算法需要選取一個降質的低分辨率圖像作為初始參考圖像，即迭代空間的初始解。理論上來說初始參考圖像可以是圖像空間中的任意一點，最后根據先驗條件都可以定位到相交的凸集。然而，為了減少迭代次數和保證收斂，希望初始解越接近最后的重建圖像越好。一般采用對降質低分辨率圖像進行插值，得到與預估計的高分辨率圖像大小相同的初始參考圖像，作為之后運動估計中其他低分辨率圖像的參考。常用的插值方法包含三線性插值、最近鄰插值以及三立方插值等。綜合考慮計算量和插值效果基礎上，本文引入兩次線性插值，先對圖像的層內采用雙線性插值，再對層間采用線性插值。線性插值是以相鄰兩層切片與待插切片的距離為權重進行線性變換，通常使待插切片位于相鄰兩切片的中間位置。
1.3 運動估計
運動估計就是求同一物體在兩幀低分辨率維圖像中的位置差，即這個物體的兩幀圖像中的“運動”。在POCS方法中，要把低分辨率圖像序列投影到參考圖像上進行修正，為保證低分辨率圖像中的點投影到參考圖像中的正確位置上就需要運動估計[10]。醫學圖像的運動估計中，解剖結構可視為剛體或是近似剛體，最常用到的變換是全局剛體變換，平移變換和旋轉變換屬于剛體變換[11]。
設點(x,y,z)在三維空間經過全局剛體變換(平移、旋轉)到點(x′,y′,z′)可表示為：
| $\left[ {\begin{array}{*{20}{c}} {x'}\\ {y'}\\ {z'} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\cos \theta }&{ - \sin \theta }&0\\ {\sin \theta }&{\cos \theta }&0\\ 0&0&1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} 1&0&0\\ 0&{\cos \omega }&{ - \sin \omega }\\ 0&{\sin \omega }&{\cos \omega } \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\cos \rho }&0&{ - \sin \rho }\\ 0&1&0\\ {\sin \rho }&0&{\cos \rho } \end{array}} \right]\left[ {\begin{array}{*{20}{c}} x\\ y\\ z \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {{t_x}}\\ {{t_y}}\\ {{t_z}} \end{array}} \right]$ |
式中：tx、ty和tz分別為X軸、Y軸和Z軸方向的平移量；θ、ω和ρ分別為XY、YZ和XZ三個平面上的旋轉角，這里只考慮圖像圍繞Z軸有很小旋轉的情況，因此變換模型可簡化為：
| $\left[ {\begin{array}{*{20}{c}} {x'}\\ {y'}\\ {z'} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\cos \theta }&{ - \sin \theta }&0\\ {\sin \theta }&{\cos \theta }&0\\ 0&0&1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} x\\ y\\ z \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {{t_x}}\\ {{t_y}}\\ {{t_z}} \end{array}} \right]$ |
令g1(x,y,z)和g2(x′,y′,z′)分別是圖像g1和g2在點(x,y,z)和(x′,y′,z′)的灰度值，成像點(x,y,z)運動到(x′,y′,z′)，則根據亮度守恒假設(即假設物體點的亮度在其運動軌跡上保持不變)有：
| ${g_1}\left( {x,y,z} \right) = {g_2}\left( {x',y',z'} \right) = {g_2}\left( {x\cos \theta - y\sin \theta + {t_x},y\cos \theta + x\sin \theta + {t_y},z + {t_z}} \right)$ |
在θ很小的情況下，cosθ≈1,sinθ≈θ，經三元泰勒展開保留前三項后則有：
| ${g_1}\left( {x,y,z} \right) \approx {g_2}\left( {x,y,z} \right) + \left( {{t_x} - y\theta } \right){g_{2x}}\left( {x,y,z} \right) + \left( {{t_y},x\theta } \right){g_{2y}}\left( {x,y,z} \right) + {t_z}{g_{2z}}\left( {x,y,z} \right)$ |
式中
| ${g_{2x}}\left( {x,y,z} \right) = \frac{{\partial {g_2}\left( {x,y,z} \right)}}{{\partial x}};{g_{2y}}\left( {x,y,z} \right) = \frac{{\partial {g_2}\left( {x,y,z} \right)}}{{\partial y}};{g_{2z}}\left( {x,y,z} \right) = \frac{{\partial {g_2}\left( {x,y,z} \right)}}{{\partial z}}$ |
為了估計運動參數()，建立最小化誤差函數：
| $\left( {{{\hat t}_x},{{\hat t}_y},{{\hat t}_z},\hat \theta } \right) = \mathop {ArgMinE}\limits_{{t_x},{t_y},{t_z},\theta } \left( {{t_x},{t_y},{t_z},\theta } \right) = \mathop {ArgMinE}\limits_{{t_x},{t_y},{t_z},\theta } {\sum\limits_{x,y,z} {\left[ {{g_2}\left( {x,y,z} \right) + \left( {{t_x} - y\theta } \right){g_{2x}}\left( {x,y,z} \right) + \left( {{t_y} + x\theta } \right){g_{2y}}\left( {x,y,z} \right) + {t_z}{g_{2z}}\left( {x,y,z} \right) - {g_1}\left( {x,y,z} \right)} \right]} ^2}$ |
分別計算E(tx,ty,tz,θ)對3個待求參數的偏導數，并令結果為0，最終求得全局運動參數為
| ${\left[ {\theta ,{t_x},{t_y},{t_z}} \right]^T} = {\left[ {\begin{array}{*{20}{c}} {\sum {{R^2}} }&{\sum {{R_{2x}}} }&{\sum {R{g_{2y}}} }&{\sum {R{g_{2z}}} }\\ {\sum {R{g_{2x}}} }&{\sum {R_{2x}^2} }&{\sum {{g_{2x}}{g_{2y}}} }&{\sum {{g_{2x}}{g_{2z}}} }\\ {\sum {R{g_{2y}}} }&{\sum {{g_{2x}}{g_{2y}}} }&{\sum {g_{2y}^2} }&{\sum {{g_{2y}}{g_{2z}}} }\\ {\sum {R{g_{2z}}} }&{\sum {{g_{2x}}{g_{2z}}} }&{\sum {{g_{2y}}{g_{2z}}} }&{\sum {g_{2z}^2} } \end{array}} \right]^{ - 1}}\left[ {\begin{array}{*{20}{c}} {\sum {R{g_t}} }\\ {\sum {{g_{2x}}{g_t}} }\\ {\sum {{g_{2y}}{g_t}} }\\ {\sum {{g_{2z}}{g_t}} } \end{array}} \right]$ |
式中：R=xg2y-yg2x;gt=g1-g2。
1.4 修正參考幀
參照二維空間中的POCS算法，利用數據一致性約束凸集對三維參考數據進行修正，重建出高分辨率數據。三維數據一致性約束凸集為：
| ${C_{{m_1},{m_2},{m_3}}} = \left\{ {\hat f\left( {{i_1},{i_2},{i_3}} \right):\left| {{r^{\left( g \right)}}\left( {{m_1},{m_2},{m_3}} \right)} \right| \le {\delta _0}} \right\}$ |
式中表示高分辨率圖像的估計值；r(g)(m1,m2,m3)表示殘差，δ0表示修正閾值。殘差的計算過程如下：
| ${r^{\left( g \right)}}\left( {{m_1},{m_2},{m_3}} \right) = g\left( {{m_1},{m_2},{m_3}} \right) - \left\{ {\sum\limits_{{i_1} = 0}^{{M_1} - 1} {\sum\limits_{{i_2} = 0}^{{M_2} - 1} {\sum\limits_{{i_3} = 0}^{{M_3} - 1} {f\left( {{i_1},{i_2},{i_3}} \right)} } } H\left( {{m_1},{m_2},{m_3};{i_1},{i_2},{i_3}} \right)} \right\}$ |
其反映了重建圖像與真實圖像之間的差別，而根據數據一致性約束的目的是要減少估計圖像與真實圖像之間的差異，因此將點x(m1,m2,m3)在Cm1,m2,m3上的投影Pi1,i2,i3[x(m1,m2,m3)]定義為：
| ${P_{{i_1},{i_2},{i_3}}}\left[ {x\left( {{m_1},{m_2},{m_3}} \right)} \right] = \left\{ {\begin{array}{*{20}{c}} {x\left( {{m_1},{m_2},{m_3}} \right) + \frac{{{r^{\left( g \right)}}\left( {{m_1},{m_2},{m_3}} \right) - {\delta _0}}}{{\sum {\sum {{H^2}\left( {{m_1},{m_2},{m_3};{i_1},{i_2},{i_3}} \right)} } }}}&{{r^{\left( g \right)}}\left( {{m_1},{m_2},{m_3}} \right) > {\delta _0}}\\ {x\left( {{m_1},{m_2},{m_3}} \right)}&{ - {\delta _0} < {r^{\left( g \right)}}\left( {{m_1},{m_2},{m_3}} \right) < {\delta _0}}\\ {x\left( {{m_1},{m_2},{m_3}} \right) + \frac{{{r^{\left( g \right)}}\left( {{m_1},{m_2},{m_3}} \right) + {\delta _0}}}{{\sum {\sum {{H^2}\left( {{m_1},{m_2},{m_3};{i_1},{i_2},{i_3}} \right)} } }}}&{{r^{\left( g \right)}}\left( {{m_1},{m_2},{m_3}} \right) < - {\delta _0}} \end{array}} \right.$ |
修正閾值δ0的選取會影響重建結果:δ0的值越大，實際修正的點越少，則引入的噪聲越少，但圖像修正質量不夠細致；而δ0的值越小，實際修正的點增多，殘差不大的一部分像素也會進行修正，引入的噪聲增多，導致重建圖像質量下降[12]。
2 結節的混合顯示方式
目前，醫學圖像三維可視化通常使用的方法為面繪制方法和體繪制方法[13]。通過對一系列二維圖像進行邊界識別等分割處理，拼接擬合物體表面來描述物體的三維結構的方法稱為面繪制方法，也稱為間接繪制方法。通過將每個體素直接投影到顯示屏來描述物體三維結構的方法稱為體繪制方法，也稱為直接繪制方法。面繪制方法在算法效率、實時交互性、圖像品質等方面有好的性能，廣泛應用于醫學圖像的可視化。在肺部CT圖像的三維可視化中，雖然相對于整個肺實質部分，結節體積小，占用的切片少，卻是醫生最關注的部分，因此，需要更高的可視化質量。采用上一節所述的超分辨率重構技術可以提高圖像的顯示質量，但圖像切片數量成倍增加，占用內存量大，以損失顯示的實時性為代價。本文使用多分辨率混合顯示方式實現結節的三維可視化。即結節與周圍組織分別繪制，并混合顯示。首先，采用圖像分割方法將結節從肺實質中提取出來，分割后剩余部分稱為背景。這樣做可以使結節與周圍的具有相似灰度級的粘連組織(血管或胸膜)分開，有利于結節的單獨顯示。然后對包含結節的圖像采用3D POCS超分辨率重建方法獲得高分辨率三維圖像，背景仍然保持低分辨率。最后，對結節和背景分別進行面繪制，并混合顯示。多分辨率混合顯示的目的是，既可保留結節足夠的細節信息，獲得高質量的肺結節顯示效果，又能在整體上降低圖像處理量，提高顯示的實時性，同時保留了結節與周圍組織之間的相對位置信息。
3 實驗結果及分析
3.1 實驗數據及環境
實驗數據來自于美國國家癌癥研究所(National Cancer Institute，NCI)發布的肺部圖像數據庫聯盟(Lung Image Database Consortium，LIDC)，該影像數據庫可以通過互聯網下載得到(https://imaging.nci.nih.gov/ncia/login.jsf/)。圖像層間分辨率為0.625~3.0 mm，層內分辨率為0.488~0.946 mm(512像素×512像素)。圖像格式為16位DICOM。仿真實驗針對5例不同患者的肺部CT圖像，按CT層厚分組，兩組為1.0 mm，圖集編號分別為190和195；兩組為1.25 mm，圖集編號分別為275和306；一組為1.5 mm，圖集編號為034。實驗的硬件平臺為3.2 GHz Pentium(R) Dual-Core CPU，內存為2.0 GB，以及NVIDIA GeForce GT 430顯卡。操作系統為WINDOWS XP，開發平臺為Matlab 7.1以及Visual Studio 2008,開發工具包VTK5.4。
3.2 實驗步驟
針對肺結節的三維超分辨率重建及混合顯示步驟：
(1)預處理：對每層圖像進行高斯濾波去噪；用Fast Marching算法提取結節。
(2)構造低分辨率圖像：按照1.1節所述方法，對提取的結節圖像構造有亞像素級位移的低分辨率三維圖像，低分辨率三維圖像數量為4，q取值2。
(3)超分辨率重建：對提取的結節圖像重建出高分辨率圖像。選擇位移為零的三維圖像作為初始參考圖像；按照1.3節給出的方法，其他3個低分辨率三維圖像分別與參考圖像進行運動估計；并根據1.4節的方法對參考圖像進行迭代修正，獲得結節的三維高分辨率圖像。
(4)對結節和背景分別用Marching Cube算法進行面繪制，并混合顯示。
針對整個肺部CT 圖像的三維超分辨率重建步驟：
(1)構造低分辨率圖像：按照1.1節所述方法，對整個CT圖像構造有亞像素級位移的低分辨率三維圖像，低分辨率三維圖像數量為4，q取值2。
(2)超分辨率重建：選擇位移為零的三維圖像作為初始參考圖像；按照1.3節給出的方法，其他3個低分辨率三維圖像分別與參考圖像進行運動估計；并根據1.4節的方法對參考圖像進行迭代修正，獲得結節的三維高分辨率圖像。
(3)用Marching Cube算法進行面繪制。可以有選擇地繪制不同的組織，如血管、骨骼。
3.3 實驗結果及評價
將本文提出的3D POCS方法與3個傳統插值方法(三立方插值、三線性插值、最近鄰插值)[4]的重建效果進行比較。重建效果評價包括主觀評價和客觀評價。選擇峰值信噪比(peak signal to noise ratio,PSNR)作為客觀評價的指標，PSNR(單位為dB)表示重建圖像與原始圖像對比的失真程度，即
| ${\rm{PSNR = }}10 \times \lg \left[ {f_{\max }^2/{\rm{MSE}}} \right]$ |
其中fmax表示圖像灰度值函數的峰值，MSE為均方誤差(mean square error)，表示重建圖像與原始圖像灰度差的平方期望值，即
| ${\rm{MSE = }}\frac{1}{{{\rm{KML}}}}\sum\limits_{{i_3} = 1}^K {\sum\limits_{{i_2} = 1}^M {\sum\limits_{{i_1} = 1}^L {{{\left[ {\hat f\left( {{i_1},{i_2},{i_3}} \right) - f\left( {{i_1},{i_2},{i_3}} \right)} \right]}^2}} } } $ |
其中K×M表示每個切片的像素數，L表示三維圖像的切片數。f(i1,i2,i3) 表示原始圖像f在點(i1,i2,i3)的灰度值，表示該點重建后的灰度值。PSNR的值越大，代表失真越少，通常，PSNR>40 dB表示效果比較好。
3.3.1 客觀評價
POCS算法中，修正閾值δ0的選取影響著重建獲得的圖像質量，需要評估δ0的有效取值范圍。針對5組結節圖像分別對δ0從0~300取16個不同的值，取值間隔為20，并計算相應的PSNR值，用于評價δ0取不同值情況下圖像重建的性能，比較結果如圖 2所示。
圖2
不同修正閾值的選取對重建性能的影響
Figure2.
Influences of different threshold selection on reconstruction of performance
圖 2表明，當δ0取值為20時，圖集034的PSNR達到最大值；然而，圖集190和195的PSNR曲線從此開始下降；對于圖集275和306，δ0取80之前的值，PSNR曲線處于遞增狀態，但增加得不明顯。總體上看，δ0在10~40之間取值時，各個圖集的PSNR值較大，表示重建效果好。
相對于肺實質部分，結節的體積較小，切片數量少，一個最大直徑30 mm的結節平均有40張切片，為了獲得更逼真的結節圖像，重建時可以給δ0賦相對較小的值。對這5組結節圖像，δ0取值為12，迭代次數為40，4種重建方法的PSNR評價指標對比如表 1所示。而整個肺部CT圖像重建時，為了提高其重建的實時性，δ0值取15，迭代次數為20。
表1
PSNR評價結節圖像的重建效果(dB)
Table1.
Evaluation of reconstruction of pulmonary nodule images using PSNR (dB)
表 1表明，本文提出的3D POCS超分辨率重建算法對于結節圖像進行重建，能夠取得比3種插值方法更好的性能。
3.3.2 主觀評價
為了更直觀地展示3D POCS算法重建肺部CT圖像的效果，針對三種情況做主觀評價。第一種情況是結節圖像重建及混合顯示。為了能在重建圖像與原圖之間作性能對比，將含有結節的圖像作為高分辨率原圖，降采樣時q=2，b=2。然后，按照3.2節所述的針對肺結節的三維重建步驟重建高分辨率結節圖像，重建圖像的分辨率應與原圖的分辨率保持一致。第二種情況是整個肺部CT圖像重建。第三種情況是整個圖像重建后對骨骼做面繪制。三種情況分別與插值方法作對比，對比結果如圖 3、圖 4和圖 5所示。
圖 3顯示以3D POCS方法獲得的高分辨率圖像，在進行三維繪制時，更接近原始高分辨率圖像，結節和骨骼組織的三維混合顯示可以展示出結節在胸部的位置信息。從主觀上判斷該方法適用于三維肺部CT圖像的高分辨率重建。
圖3
結節圖像重建及顯示效果對比
Figure3.
Comparison of reconstruction and 3D visualization of pulmonary nodules by different methods
從圖 4可以看出，與其他插值方法比較，3D POCS算法重建的圖像整體更清晰；高灰度級的血管、結節與低灰度級的肺實質之間的對比度更高；圓形區域中結節的邊緣細節(如毛刺特征)更明顯；矩形區域中肺裂線更容易辨識。
圖4
整個肺部CT圖像重建效果對比(其中一層)
圓形區域中大的亮色塊是結節；矩形區域中亮線是肺裂線。(a)三立方插值算法；(b)三線性插值算法；(c)最近鄰插值算法；(d)3D POCS算法
Figure4. Comparison of reconstruction of the whole lung CT images by different methods (one of the given layers)nodule in circles and fissure in rectangles shown as high gray level in each sub image. (a) cubic interpolation; (b) tri-interpolation; (c) nearest interpolation; (d) 3D POCS
圖 5顯示了圖像重建后骨骼面繪制的效果。三線性插值算法和最鄰近插值算法重建后圖像的質量沒能得到更好的改善，骨骼表面不夠光滑，仍然有明顯的梯田狀粗紋理現象。三立方插值法重建的效果優于其他兩種插值法，但局部區域仍然不夠細膩，如圖 5(a)中的兩處橢圓形區域內，可看到骨骼表面粗糙的紋理。3D POCS算法重建效果優于其他三種插值方法，有效地降低了骨骼表面的梯田狀紋理，改善了圖像整體質量。
圖5
骨骼圖像重建及面繪制效果對比
(a)三立方插值算法；(b)三線性插值算法；(c)最近鄰插值算法；(d)3D POCS算法
Figure5. Comparison of bone reconstruction and surface rendering by different methods(a) cubic interpolation; (b) tri-interpolation; (c) nearest interpolation; (d) 3D POCS
4 結論
提出三維凸集投影(3D POCS)算法，用于實現三維肺部CT圖像的超分辨率重建；并采用多分辨率混合顯示方式實現肺結節的三維可視化。實驗選取5組圖集做性能評價，兩組切片層厚為1.0 mm，兩組為1.25 mm，一組為1.5 mm，圖集中切片內像素數為512×512。與3種插值方法進行比較，主客觀兩方面的評價顯示，3D POCS算法實現三維圖像的超分辨率重建性能優于其他方法；混合顯示方式能夠滿足肺結節高分辨率三維可視化的需要。
引言
根據世界衛生組織2014年5月發布的最新資料顯示，肺癌是全球前十位死亡原因的第五位，癌癥死亡的首位[1]。目前，醫學上對肺癌的早期檢測主要利用計算機斷層掃描技術(computer tomography,CT)。國內外研究表明，計算機輔助檢測/診斷(computer assisted detection/diagnosis，CAD)系統對肺部CT圖像進行自動分析，不僅減少了醫生的工作量，而且提高了檢測和診斷的準確性，為肺癌的早期檢測和診斷提供了有力的支持[2]。在CAD系統中，像病變部位及周圍的器官組織等醫生感興趣區域(regions of interest，ROI)的三維逼真可視化具有重要的意義，醫生可以交互地從多個角度觀察ROI大小、形態特征，獲得病變部位更為準確和客觀的信息以做出更為準確的診斷[3]。然而，由于CT固有的高劑量照射，沿Z軸方向的密集釆樣往往是不實際的，導致層間分辨率遠低于層內分辨率，造成數據顯著的各項異性。另外，由于結節體積相對較小(小結節通常直徑不超過5 mm)，對應的切片數量少，直接進行三維繪制難以達到醫學診斷所需的真實感效果。因此，需要三維圖像沿Z軸增加層與層之間的信息量，以獲得肺結節高清晰度可視化效果。采用插值方法增加層間切片數量，如線性插值、最近鄰插值和立體卷積插值等，雖然獲得的三維圖像體素數有所提高,但圖像所包含的頻率信息并沒有什么變化,因此圖像并沒有增加細節信息，三維繪制之后物體表面仍會存在明顯的梯田狀粗紋理效果[4]。增加圖像細節信息的一種有效方法是超分辨率重建技術，即由同一場景的多幀具有互補信息的低分辨率圖像序列，按照相應的算法獲得一幅或多幅清晰的高分辨率圖像的過程。Harris和Goodman于20世紀60年代初提出了針對單張圖像重建的超分辨率重建技術的概念和方法，之后許多研究者在此基礎上相繼提出橢圓球波函數法、線性外推法以及疊加正弦模板法等；1984年Tsai和Huang引入多幀增強的概念，并且首次在頻域中提出基于序列圖像的超分辨率重建的概念和方法[5]。利用序列和多幅圖像進行超分辨率重建成為研究熱點。20世紀90年代以來，各種空間域算法的提出，使得超分辨率重建技術得到了更大的發展。空間域算法主要有迭代反投影方法[6]、凸集投影法(projection onto convex sets,POCS)[7]、統計復原方法[8](最大后驗概念估計器和最大似然估計器)等。基于學習的超分辨率重建方法從大量的訓練樣本集中獲取先驗知識作為超分辨率的依據，在不增加輸入圖像樣本數量的前提下，產生新的高頻細節，提高重建圖像的質量。由于層析成像技術的特殊機制，使得超分辨率重建技術在醫學圖像處理領域獲得廣泛的研究[9]。基于集合理論重建的理論體系是超分辨率重建領域中的一種主流方法，其實現簡單，能方便地加入先驗信息，可很好地保持高分辨率圖像上的邊緣和細節，適用于醫學圖像中局部精確信息的復原或重建[7]。然而傳統的POCS算法是針對二維圖像的超分辨率重建，只考慮圖像內部信息而忽略了圖像之間的關聯信息。對于醫學斷層掃描圖像，由于病灶特征信息不僅存在于各個切片中(層內)，也存在于多個連續的切片之間(層間)，因此，需要一種基于更高維的超分辨率重建方法。本文將二維圖像超分辨率重建的POCS方法擴展至三維數據場，提出三維POCS方法，即利用多個有亞像素級位移的低分辨率三維圖像間的信息重建出高分辨率的三維圖像，用于三維肺部CT圖像的超分辨率重建；并采用多分辨率混合顯示方式實現肺結節的高清晰度三維可視化。
1 三維凸集投影方法
POCS方法是使用一系列凸約束集合來描述超分辨率重建方法中的一些特性，如正定性、能量有界性、平滑性等，這些約束集對高分辨率圖像的解進行約束，用于保證與觀測數據(低分辨率圖像)的一致性，通常稱之為一致性約束凸集，超分辨率圖像解空間與這一系列的凸約束集合相交。POCS方法是一個迭代過程，預先估計一幀高分辨率圖像作為成像空間中的初始參考圖像，然后從參考圖像的某一點開始，將圖像當前估計值投影到一致性約束凸集上，判斷是否滿足所有的凸形約束，如果不滿足，則將其殘差反投影到高分辨率估計圖像上對其進行修正，經過有限次的迭代，直到找到滿足所有凸約束集的解，重建出高分辨率圖像。因此，POCS重建算法主要分為三個具體步驟：構造參考圖像、進行運動估計、對參考圖像進行修正。在一個病例僅有一個三維肺部CT圖像的情況下，重建之前需要構造有亞像素級位移的多個低分辨率三維圖像，模擬多個相關的、有運動位移的、反映同一場景的低分辨率三維圖像，由多個這樣的三維圖像重建一個高分辨率的三維圖像。
1.1 構造有亞像素級位移的低分辨率三維圖像
在超分辨率重建中，普遍接受的退化成像過程模型為“變形-模糊”過程，結合肺部橫截面CT圖像特點，假定圖像之間的相對運動(圖像變形)采用全局平移運動模擬。設理想高分辨率三維圖像為g(x,y,z)，其空間分辨率為qN×qN×bP，其中qN×qN為層內像素數，bP為層數，q和b為降采樣因子或放大倍數。對于三維肺部CT圖像，通常bP>qN。為了盡可能獲得正確比例的三維圖像，圖像層內和層間使用不同的降采樣間隔。為了模擬多個相關的、有運動位移的、反映同一場景的低分辨率三維圖像，需要對g(x,y,z)進行多次降采樣，并且兩次降采樣之間移動小于一個像素的距離，稱為亞像素級位移，各圖像之間的相對位移保持一致。
圖 1為4個互有1/2像素位移的降采樣過程示意圖，移動路線分別為不移動、X軸正向移動、Y軸正向移動、X軸負向移動。虛線為高分辨率數據，實線為降采樣數據，降采樣因子q=2，b=1，得到分辨率為N×N×bP的多個低分辨率圖像。即，通過每次整體錯開一個高分辨率像素大小的位移，經過降采樣后獲取的低分辨率圖像之間的相對位移即為亞像素級位移，從而獲得有亞像素級位移的一組低分辨率數據。
圖1
4個有亞像素級位移的三維低分辨率數據的構建過程
Figure1.
Construction of 4 sub-pixel 3D low-resolution data
1.2 構造參考圖像
POCS算法需要選取一個降質的低分辨率圖像作為初始參考圖像，即迭代空間的初始解。理論上來說初始參考圖像可以是圖像空間中的任意一點，最后根據先驗條件都可以定位到相交的凸集。然而，為了減少迭代次數和保證收斂，希望初始解越接近最后的重建圖像越好。一般采用對降質低分辨率圖像進行插值，得到與預估計的高分辨率圖像大小相同的初始參考圖像，作為之后運動估計中其他低分辨率圖像的參考。常用的插值方法包含三線性插值、最近鄰插值以及三立方插值等。綜合考慮計算量和插值效果基礎上，本文引入兩次線性插值，先對圖像的層內采用雙線性插值，再對層間采用線性插值。線性插值是以相鄰兩層切片與待插切片的距離為權重進行線性變換，通常使待插切片位于相鄰兩切片的中間位置。
1.3 運動估計
運動估計就是求同一物體在兩幀低分辨率維圖像中的位置差，即這個物體的兩幀圖像中的“運動”。在POCS方法中，要把低分辨率圖像序列投影到參考圖像上進行修正，為保證低分辨率圖像中的點投影到參考圖像中的正確位置上就需要運動估計[10]。醫學圖像的運動估計中，解剖結構可視為剛體或是近似剛體，最常用到的變換是全局剛體變換，平移變換和旋轉變換屬于剛體變換[11]。
設點(x,y,z)在三維空間經過全局剛體變換(平移、旋轉)到點(x′,y′,z′)可表示為：
| $\left[ {\begin{array}{*{20}{c}} {x'}\\ {y'}\\ {z'} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\cos \theta }&{ - \sin \theta }&0\\ {\sin \theta }&{\cos \theta }&0\\ 0&0&1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} 1&0&0\\ 0&{\cos \omega }&{ - \sin \omega }\\ 0&{\sin \omega }&{\cos \omega } \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\cos \rho }&0&{ - \sin \rho }\\ 0&1&0\\ {\sin \rho }&0&{\cos \rho } \end{array}} \right]\left[ {\begin{array}{*{20}{c}} x\\ y\\ z \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {{t_x}}\\ {{t_y}}\\ {{t_z}} \end{array}} \right]$ |
式中：tx、ty和tz分別為X軸、Y軸和Z軸方向的平移量；θ、ω和ρ分別為XY、YZ和XZ三個平面上的旋轉角，這里只考慮圖像圍繞Z軸有很小旋轉的情況，因此變換模型可簡化為：
| $\left[ {\begin{array}{*{20}{c}} {x'}\\ {y'}\\ {z'} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\cos \theta }&{ - \sin \theta }&0\\ {\sin \theta }&{\cos \theta }&0\\ 0&0&1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} x\\ y\\ z \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {{t_x}}\\ {{t_y}}\\ {{t_z}} \end{array}} \right]$ |
令g1(x,y,z)和g2(x′,y′,z′)分別是圖像g1和g2在點(x,y,z)和(x′,y′,z′)的灰度值，成像點(x,y,z)運動到(x′,y′,z′)，則根據亮度守恒假設(即假設物體點的亮度在其運動軌跡上保持不變)有：
| ${g_1}\left( {x,y,z} \right) = {g_2}\left( {x',y',z'} \right) = {g_2}\left( {x\cos \theta - y\sin \theta + {t_x},y\cos \theta + x\sin \theta + {t_y},z + {t_z}} \right)$ |
在θ很小的情況下，cosθ≈1,sinθ≈θ，經三元泰勒展開保留前三項后則有：
| ${g_1}\left( {x,y,z} \right) \approx {g_2}\left( {x,y,z} \right) + \left( {{t_x} - y\theta } \right){g_{2x}}\left( {x,y,z} \right) + \left( {{t_y},x\theta } \right){g_{2y}}\left( {x,y,z} \right) + {t_z}{g_{2z}}\left( {x,y,z} \right)$ |
式中
| ${g_{2x}}\left( {x,y,z} \right) = \frac{{\partial {g_2}\left( {x,y,z} \right)}}{{\partial x}};{g_{2y}}\left( {x,y,z} \right) = \frac{{\partial {g_2}\left( {x,y,z} \right)}}{{\partial y}};{g_{2z}}\left( {x,y,z} \right) = \frac{{\partial {g_2}\left( {x,y,z} \right)}}{{\partial z}}$ |
為了估計運動參數()，建立最小化誤差函數：
| $\left( {{{\hat t}_x},{{\hat t}_y},{{\hat t}_z},\hat \theta } \right) = \mathop {ArgMinE}\limits_{{t_x},{t_y},{t_z},\theta } \left( {{t_x},{t_y},{t_z},\theta } \right) = \mathop {ArgMinE}\limits_{{t_x},{t_y},{t_z},\theta } {\sum\limits_{x,y,z} {\left[ {{g_2}\left( {x,y,z} \right) + \left( {{t_x} - y\theta } \right){g_{2x}}\left( {x,y,z} \right) + \left( {{t_y} + x\theta } \right){g_{2y}}\left( {x,y,z} \right) + {t_z}{g_{2z}}\left( {x,y,z} \right) - {g_1}\left( {x,y,z} \right)} \right]} ^2}$ |
分別計算E(tx,ty,tz,θ)對3個待求參數的偏導數，并令結果為0，最終求得全局運動參數為
| ${\left[ {\theta ,{t_x},{t_y},{t_z}} \right]^T} = {\left[ {\begin{array}{*{20}{c}} {\sum {{R^2}} }&{\sum {{R_{2x}}} }&{\sum {R{g_{2y}}} }&{\sum {R{g_{2z}}} }\\ {\sum {R{g_{2x}}} }&{\sum {R_{2x}^2} }&{\sum {{g_{2x}}{g_{2y}}} }&{\sum {{g_{2x}}{g_{2z}}} }\\ {\sum {R{g_{2y}}} }&{\sum {{g_{2x}}{g_{2y}}} }&{\sum {g_{2y}^2} }&{\sum {{g_{2y}}{g_{2z}}} }\\ {\sum {R{g_{2z}}} }&{\sum {{g_{2x}}{g_{2z}}} }&{\sum {{g_{2y}}{g_{2z}}} }&{\sum {g_{2z}^2} } \end{array}} \right]^{ - 1}}\left[ {\begin{array}{*{20}{c}} {\sum {R{g_t}} }\\ {\sum {{g_{2x}}{g_t}} }\\ {\sum {{g_{2y}}{g_t}} }\\ {\sum {{g_{2z}}{g_t}} } \end{array}} \right]$ |
式中：R=xg2y-yg2x;gt=g1-g2。
1.4 修正參考幀
參照二維空間中的POCS算法，利用數據一致性約束凸集對三維參考數據進行修正，重建出高分辨率數據。三維數據一致性約束凸集為：
| ${C_{{m_1},{m_2},{m_3}}} = \left\{ {\hat f\left( {{i_1},{i_2},{i_3}} \right):\left| {{r^{\left( g \right)}}\left( {{m_1},{m_2},{m_3}} \right)} \right| \le {\delta _0}} \right\}$ |
式中表示高分辨率圖像的估計值；r(g)(m1,m2,m3)表示殘差，δ0表示修正閾值。殘差的計算過程如下：
| ${r^{\left( g \right)}}\left( {{m_1},{m_2},{m_3}} \right) = g\left( {{m_1},{m_2},{m_3}} \right) - \left\{ {\sum\limits_{{i_1} = 0}^{{M_1} - 1} {\sum\limits_{{i_2} = 0}^{{M_2} - 1} {\sum\limits_{{i_3} = 0}^{{M_3} - 1} {f\left( {{i_1},{i_2},{i_3}} \right)} } } H\left( {{m_1},{m_2},{m_3};{i_1},{i_2},{i_3}} \right)} \right\}$ |
其反映了重建圖像與真實圖像之間的差別，而根據數據一致性約束的目的是要減少估計圖像與真實圖像之間的差異，因此將點x(m1,m2,m3)在Cm1,m2,m3上的投影Pi1,i2,i3[x(m1,m2,m3)]定義為：
| ${P_{{i_1},{i_2},{i_3}}}\left[ {x\left( {{m_1},{m_2},{m_3}} \right)} \right] = \left\{ {\begin{array}{*{20}{c}} {x\left( {{m_1},{m_2},{m_3}} \right) + \frac{{{r^{\left( g \right)}}\left( {{m_1},{m_2},{m_3}} \right) - {\delta _0}}}{{\sum {\sum {{H^2}\left( {{m_1},{m_2},{m_3};{i_1},{i_2},{i_3}} \right)} } }}}&{{r^{\left( g \right)}}\left( {{m_1},{m_2},{m_3}} \right) > {\delta _0}}\\ {x\left( {{m_1},{m_2},{m_3}} \right)}&{ - {\delta _0} < {r^{\left( g \right)}}\left( {{m_1},{m_2},{m_3}} \right) < {\delta _0}}\\ {x\left( {{m_1},{m_2},{m_3}} \right) + \frac{{{r^{\left( g \right)}}\left( {{m_1},{m_2},{m_3}} \right) + {\delta _0}}}{{\sum {\sum {{H^2}\left( {{m_1},{m_2},{m_3};{i_1},{i_2},{i_3}} \right)} } }}}&{{r^{\left( g \right)}}\left( {{m_1},{m_2},{m_3}} \right) < - {\delta _0}} \end{array}} \right.$ |
修正閾值δ0的選取會影響重建結果:δ0的值越大，實際修正的點越少，則引入的噪聲越少，但圖像修正質量不夠細致；而δ0的值越小，實際修正的點增多，殘差不大的一部分像素也會進行修正，引入的噪聲增多，導致重建圖像質量下降[12]。
2 結節的混合顯示方式
目前，醫學圖像三維可視化通常使用的方法為面繪制方法和體繪制方法[13]。通過對一系列二維圖像進行邊界識別等分割處理，拼接擬合物體表面來描述物體的三維結構的方法稱為面繪制方法，也稱為間接繪制方法。通過將每個體素直接投影到顯示屏來描述物體三維結構的方法稱為體繪制方法，也稱為直接繪制方法。面繪制方法在算法效率、實時交互性、圖像品質等方面有好的性能，廣泛應用于醫學圖像的可視化。在肺部CT圖像的三維可視化中，雖然相對于整個肺實質部分，結節體積小，占用的切片少，卻是醫生最關注的部分，因此，需要更高的可視化質量。采用上一節所述的超分辨率重構技術可以提高圖像的顯示質量，但圖像切片數量成倍增加，占用內存量大，以損失顯示的實時性為代價。本文使用多分辨率混合顯示方式實現結節的三維可視化。即結節與周圍組織分別繪制，并混合顯示。首先，采用圖像分割方法將結節從肺實質中提取出來，分割后剩余部分稱為背景。這樣做可以使結節與周圍的具有相似灰度級的粘連組織(血管或胸膜)分開，有利于結節的單獨顯示。然后對包含結節的圖像采用3D POCS超分辨率重建方法獲得高分辨率三維圖像，背景仍然保持低分辨率。最后，對結節和背景分別進行面繪制，并混合顯示。多分辨率混合顯示的目的是，既可保留結節足夠的細節信息，獲得高質量的肺結節顯示效果，又能在整體上降低圖像處理量，提高顯示的實時性，同時保留了結節與周圍組織之間的相對位置信息。
3 實驗結果及分析
3.1 實驗數據及環境
實驗數據來自于美國國家癌癥研究所(National Cancer Institute，NCI)發布的肺部圖像數據庫聯盟(Lung Image Database Consortium，LIDC)，該影像數據庫可以通過互聯網下載得到(https://imaging.nci.nih.gov/ncia/login.jsf/)。圖像層間分辨率為0.625~3.0 mm，層內分辨率為0.488~0.946 mm(512像素×512像素)。圖像格式為16位DICOM。仿真實驗針對5例不同患者的肺部CT圖像，按CT層厚分組，兩組為1.0 mm，圖集編號分別為190和195；兩組為1.25 mm，圖集編號分別為275和306；一組為1.5 mm，圖集編號為034。實驗的硬件平臺為3.2 GHz Pentium(R) Dual-Core CPU，內存為2.0 GB，以及NVIDIA GeForce GT 430顯卡。操作系統為WINDOWS XP，開發平臺為Matlab 7.1以及Visual Studio 2008,開發工具包VTK5.4。
3.2 實驗步驟
針對肺結節的三維超分辨率重建及混合顯示步驟：
(1)預處理：對每層圖像進行高斯濾波去噪；用Fast Marching算法提取結節。
(2)構造低分辨率圖像：按照1.1節所述方法，對提取的結節圖像構造有亞像素級位移的低分辨率三維圖像，低分辨率三維圖像數量為4，q取值2。
(3)超分辨率重建：對提取的結節圖像重建出高分辨率圖像。選擇位移為零的三維圖像作為初始參考圖像；按照1.3節給出的方法，其他3個低分辨率三維圖像分別與參考圖像進行運動估計；并根據1.4節的方法對參考圖像進行迭代修正，獲得結節的三維高分辨率圖像。
(4)對結節和背景分別用Marching Cube算法進行面繪制，并混合顯示。
針對整個肺部CT 圖像的三維超分辨率重建步驟：
(1)構造低分辨率圖像：按照1.1節所述方法，對整個CT圖像構造有亞像素級位移的低分辨率三維圖像，低分辨率三維圖像數量為4，q取值2。
(2)超分辨率重建：選擇位移為零的三維圖像作為初始參考圖像；按照1.3節給出的方法，其他3個低分辨率三維圖像分別與參考圖像進行運動估計；并根據1.4節的方法對參考圖像進行迭代修正，獲得結節的三維高分辨率圖像。
(3)用Marching Cube算法進行面繪制。可以有選擇地繪制不同的組織，如血管、骨骼。
3.3 實驗結果及評價
將本文提出的3D POCS方法與3個傳統插值方法(三立方插值、三線性插值、最近鄰插值)[4]的重建效果進行比較。重建效果評價包括主觀評價和客觀評價。選擇峰值信噪比(peak signal to noise ratio,PSNR)作為客觀評價的指標，PSNR(單位為dB)表示重建圖像與原始圖像對比的失真程度，即
| ${\rm{PSNR = }}10 \times \lg \left[ {f_{\max }^2/{\rm{MSE}}} \right]$ |
其中fmax表示圖像灰度值函數的峰值，MSE為均方誤差(mean square error)，表示重建圖像與原始圖像灰度差的平方期望值，即
| ${\rm{MSE = }}\frac{1}{{{\rm{KML}}}}\sum\limits_{{i_3} = 1}^K {\sum\limits_{{i_2} = 1}^M {\sum\limits_{{i_1} = 1}^L {{{\left[ {\hat f\left( {{i_1},{i_2},{i_3}} \right) - f\left( {{i_1},{i_2},{i_3}} \right)} \right]}^2}} } } $ |
其中K×M表示每個切片的像素數，L表示三維圖像的切片數。f(i1,i2,i3) 表示原始圖像f在點(i1,i2,i3)的灰度值，表示該點重建后的灰度值。PSNR的值越大，代表失真越少，通常，PSNR>40 dB表示效果比較好。
3.3.1 客觀評價
POCS算法中，修正閾值δ0的選取影響著重建獲得的圖像質量，需要評估δ0的有效取值范圍。針對5組結節圖像分別對δ0從0~300取16個不同的值，取值間隔為20，并計算相應的PSNR值，用于評價δ0取不同值情況下圖像重建的性能，比較結果如圖 2所示。
圖2
不同修正閾值的選取對重建性能的影響
Figure2.
Influences of different threshold selection on reconstruction of performance
圖 2表明，當δ0取值為20時，圖集034的PSNR達到最大值；然而，圖集190和195的PSNR曲線從此開始下降；對于圖集275和306，δ0取80之前的值，PSNR曲線處于遞增狀態，但增加得不明顯。總體上看，δ0在10~40之間取值時，各個圖集的PSNR值較大，表示重建效果好。
相對于肺實質部分，結節的體積較小，切片數量少，一個最大直徑30 mm的結節平均有40張切片，為了獲得更逼真的結節圖像，重建時可以給δ0賦相對較小的值。對這5組結節圖像，δ0取值為12，迭代次數為40，4種重建方法的PSNR評價指標對比如表 1所示。而整個肺部CT圖像重建時，為了提高其重建的實時性，δ0值取15，迭代次數為20。
表1
PSNR評價結節圖像的重建效果(dB)
Table1.
Evaluation of reconstruction of pulmonary nodule images using PSNR (dB)
表 1表明，本文提出的3D POCS超分辨率重建算法對于結節圖像進行重建，能夠取得比3種插值方法更好的性能。
3.3.2 主觀評價
為了更直觀地展示3D POCS算法重建肺部CT圖像的效果，針對三種情況做主觀評價。第一種情況是結節圖像重建及混合顯示。為了能在重建圖像與原圖之間作性能對比，將含有結節的圖像作為高分辨率原圖，降采樣時q=2，b=2。然后，按照3.2節所述的針對肺結節的三維重建步驟重建高分辨率結節圖像，重建圖像的分辨率應與原圖的分辨率保持一致。第二種情況是整個肺部CT圖像重建。第三種情況是整個圖像重建后對骨骼做面繪制。三種情況分別與插值方法作對比，對比結果如圖 3、圖 4和圖 5所示。
圖 3顯示以3D POCS方法獲得的高分辨率圖像，在進行三維繪制時，更接近原始高分辨率圖像，結節和骨骼組織的三維混合顯示可以展示出結節在胸部的位置信息。從主觀上判斷該方法適用于三維肺部CT圖像的高分辨率重建。
圖3
結節圖像重建及顯示效果對比
Figure3.
Comparison of reconstruction and 3D visualization of pulmonary nodules by different methods
從圖 4可以看出，與其他插值方法比較，3D POCS算法重建的圖像整體更清晰；高灰度級的血管、結節與低灰度級的肺實質之間的對比度更高；圓形區域中結節的邊緣細節(如毛刺特征)更明顯；矩形區域中肺裂線更容易辨識。
圖4
整個肺部CT圖像重建效果對比(其中一層)
圓形區域中大的亮色塊是結節；矩形區域中亮線是肺裂線。(a)三立方插值算法；(b)三線性插值算法；(c)最近鄰插值算法；(d)3D POCS算法
Figure4. Comparison of reconstruction of the whole lung CT images by different methods (one of the given layers)nodule in circles and fissure in rectangles shown as high gray level in each sub image. (a) cubic interpolation; (b) tri-interpolation; (c) nearest interpolation; (d) 3D POCS
圖 5顯示了圖像重建后骨骼面繪制的效果。三線性插值算法和最鄰近插值算法重建后圖像的質量沒能得到更好的改善，骨骼表面不夠光滑，仍然有明顯的梯田狀粗紋理現象。三立方插值法重建的效果優于其他兩種插值法，但局部區域仍然不夠細膩，如圖 5(a)中的兩處橢圓形區域內，可看到骨骼表面粗糙的紋理。3D POCS算法重建效果優于其他三種插值方法，有效地降低了骨骼表面的梯田狀紋理，改善了圖像整體質量。
圖5
骨骼圖像重建及面繪制效果對比
(a)三立方插值算法；(b)三線性插值算法；(c)最近鄰插值算法；(d)3D POCS算法
Figure5. Comparison of bone reconstruction and surface rendering by different methods(a) cubic interpolation; (b) tri-interpolation; (c) nearest interpolation; (d) 3D POCS
4 結論
提出三維凸集投影(3D POCS)算法，用于實現三維肺部CT圖像的超分辨率重建；并采用多分辨率混合顯示方式實現肺結節的三維可視化。實驗選取5組圖集做性能評價，兩組切片層厚為1.0 mm，兩組為1.25 mm，一組為1.5 mm，圖集中切片內像素數為512×512。與3種插值方法進行比較，主客觀兩方面的評價顯示，3D POCS算法實現三維圖像的超分辨率重建性能優于其他方法；混合顯示方式能夠滿足肺結節高分辨率三維可視化的需要。
