提出三维凸集投影(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算法实现三维图像的超分辨率重建性能优于其他方法；混合显示方式能够满足肺结节高分辨率三维可视化的需要。
