基于运动估计的肺4D-CT图像冠矢状面超分辨率重建

肖 珊 王婷婷 张 煜

(南方医科大学生物医学工程学院，广州 510515)

引言

肺癌是常见的恶性肿瘤之一，也是确诊后存活率最低的癌症之一，特别是近几十年以来，各国肺癌的发病率和死亡率都在急剧上升，是名副其实的健康杀手。放射治疗是治疗肺癌的有效手段之一，而指导放射治疗计划和剂量投射的信息通常来自于肺自由呼吸3D-CT的图像。然而，3D数据没有呼吸相位信息，这意味着在确定肺组织的结构和形状时，缺乏足够的运动信息。因此，为了在肺肿瘤放射治疗时，确保有足够的肿瘤覆盖面，肿瘤估计区域一般是根据目标运动向各个方向扩展，再进行放射治疗，而对靶体积不必要的扩大治疗将伤害正常组织。

基于肺4D-CT图像的精确放射治疗是当前治疗肺癌的有效手段，肺4D-CT提供了一个全面的高精度放射治疗呼吸运动表征。在肺4D-CT数据中，由于有多个相位的图像，各相位图像可以有助于获取呼吸运动信息，是放射治疗目标精确定位的关键，因此肺4D-CT在肺肿瘤精确放射治疗中发挥着越来越重要的作用。肺4D-CT数据的获取通常是根据床位和肺体积，将多个自由呼吸3D-CT数据段排序而得。然而，由于CT固有的高剂量照射[1-2]，沿纵向(Z轴方向)的密集采样往往是不实际的，从而导致层内分辨率远低于层间分辨率，造成数据显著的各向异性。因此，在对各相位3D数据进行冠矢状面观察时，为了获得正确比例的图像，需要根据3D数据的层间分辨率和层内分辨率的比例，沿Z轴方向进行插值放大。常用的插值方法为最近邻或双线性插值法，这些方法都会导致图像模糊，尤其当层间分辨率与层内分辨率的比例差别比较大时。笔者针对此问题，根据肺4D-CT图像数据的特点，采用超分辨率技术重建清晰的肺冠矢状面图像。

图像超分辨率技术是指利用多帧关于同一场景的有相互位移的低分辨率降质图像来重建高分辨率高质量图像的技术[3]。肺4D-CT图像，是由10-20个不同相位的3D图像所组成，其每一相位对应不同的肺运动时刻。因此对于某一3D图像的某一位置的低分辨率(Z轴方向)冠(矢)状面图像，其他相位对应位置的冠(矢)状面图像可认为是相关的、有运动位移的多“帧”同一场景低分辨率图像，故可采用图像超分辨率重建技术来获得高分辨率肺4D-CT冠(矢)状面图像，这是本项工作的基本思想。

针对肺4D-CT冠矢状面超分辨率重建问题，本文首先分析了图像退化基本模型;随后，采用基于完全搜索块匹配算法[4]估计不同“帧”肺冠矢状面图像之间的运动场;最后，以此运动场为基础，采用迭代反投影(IBP)法[5]重建高分辨率肺4D-CT冠矢状面图像。实验结果表明，所提出的方法可获得更清晰的结果图像。

1 理论与方法

1.1 图像退化模型和超分辨率重建原理

图像退化基本模型[6-7]如图1所示。

图1 图像的退化模型Fig.1 The image degradation model

这个模型可用公式表示为

式中:gk表示K幅低分辨率图像中的第k幅，f表示需要重建的高分辨率图像;Mk表示几何变换，由运动估计求得;Bk表示模糊矩阵，是由光学系统本身、成像系统与原始场景的相对运动以及低分辨率传感器的点扩散函数(point spread function，PSF)形成的;D表示下采样矩阵;ηk是系统加性噪声。

本研究的超分辨率图像重建是这个退化过程的逆问题。也就是说，利用多帧关于同一场景(不同相位)的有相互位移的低分辨率肺4D-CT冠矢状面图像，可以重建高分辨率清晰肺冠矢状面图像。图2显示了一例肺4D-CT图像中的3个相位(分别为相位7、2、0)对应的同一位置肺冠状面图像，其中黑色带箭头实线标示了相位间的呼吸运动趋势。

图2 不同相位对应的冠状面图像Fig.2 The corresponding corona images of different phases

1.2 不同相位的肺4D-CT冠矢状面图像运动估计

图像块匹配算法是目前一些国际标准组织推荐的运动估计方案，也是图像序列稳定中最常用的一种运动估计算法。笔者采用完全搜索块匹配算法，对肺4D-CT不同相位对应的冠矢状面图像之间进行运动场估计。

完全搜索块匹配算法原理如图3所示。假设块内各像素做相同运动，根据一定的匹配准则，在前一帧(即参考帧)的某一给定搜索区域内，找出与当前帧中某一块最相似的块，即匹配块[8]。由匹配块与当前块的相对位置，计算出运动位移。

图3 块匹配运动估计原理Fig.3 Schematic of block-matching motion estimation

完全搜索块匹配算法有两种常用的匹配准则:最小均方误差匹配准则和最小绝对误差匹配准则。笔者采用的是常用的最小绝对误差匹配准则[9]，定义为

式中:图像块的大小为N×N，左上角的坐标为(p，q)，本研究选取的块大小为16×16;运动矢量为(vx，vy);It(i，j)和 It-Δt(i，j)分别为当前帧和参考帧在像素(i，j)处的值。在搜索区域内使S(vx，vy)值最小的运动矢量即为当前块的最优运动矢量。

实验前，从某组肺数据中分别提取10个相位某位置对应的冠状面和矢状面图像，作为实验初始低分辨率图像，大小均为256像素×99像素，并将其插值放大到256像素×216像素。如图4所示，(a)为一例肺部冠状面初始低分辨率图像，(b)为插值到256像素×216像素的图像。从图中可看出，插值放大后的图像较模糊。

图4 肺4D-CT数据的冠状面图像。(a)原始低分辨率图像;(b)插值后的图像Fig.4 Example coronal view of 4D-CT lung data.(a)The original low-resolution image;(b)Interpolated image

由此，取其中需重建的某相位插值图像作为初始假设的高分辨率图;然后，将其与其他相位插值放大的图像通过运动估计分别得到对应各相位图像之间的运动矢量场，进而由迭代反投影算法重建出高分辨率图像。

1.3 肺4D-CT冠矢状面图像迭代反投影超分辨率重建

迭代反投影(iterative back-projection，IBP)算法是由Irani和Peleg率先提出[10]，主要过程为:首先，假设一个初始高分辨率图像f(0)，将需要重建的某一原始低分辨率图像插值放大为一幅初始高分辨率图像;然后，根据退化模型模拟成像过程，得到一个低分辨图像的集合(k=1，...，K)，K 表示原始序列低分辨率图像的数量，n为迭代次数。如果f(0)是正确的高分辨率图像，则得到的这一集合就应该和已知的低分辨率图像序列集合{gk}相等。如果f(0)不是理想的高分辨率图像，则{gk}与将产生差值{gk-}。将此差值返回给f(0)并将其更新，得到一个相对较接近的正确结果。反复迭代该过程，使误差函数达到最小值，即

式中:Tk表示从f到gk的二维几何变换，且可逆，即为运动估计获得的运动矢量(vx，vy);↓s是下采样算子。

重建的高分辨率图像的更新过程可表示为

式中，↑s表示上采样算子，p表示背投影算子。

为使算法收敛，p满足条件‖δ-h*p‖2≤1，其中δ为单一脉冲;基于稳定性和收敛速度的考虑[11]，设定p=h-1。h是模糊算子，采用高斯模型，标准差设为3。此标准差值通过实验获得，因为肺图像较模糊，故高斯退化模型的标准差亦较大。

1.4 超分辨率重建结果量化评价

本实验采用一个公共可用的数据集[11]，该数据集由10组肺4D-CT数据组成，每组数据包含10个相位。每组数据中的各相位3D图像大小分别为256像素×256像素×94像素～256像素×256像素×112像素，图像层内分辨率为0.97～1.16 mm，层间分辨率均为2.5 mm。实验针对10组数据，在不同相位中选取冠矢状面低分辨率图像，然后分别采用最近邻插值、双线性插值和本算法进行图像重建。由于剂量限制，临床上尚无法采集高分辨率的4D-CT原始数据。

本研究采用边缘宽度[12]和图像平均梯度[13]这两个量化评价指标来评价重建效果。

边缘宽度(edge widths)是衡量图像空间分辨率的量化参数，其值越小，边界越清晰，分辨率越高，文献[12]给出了它的计算方法。

首先，沿图像边缘采集一系列像素点，并拟合一双曲线函数，有

式中，x为像素点坐标，y(x)为像素点灰度，c为采集的像素点中心，a为与边缘宽度成反比的参数。

根据点集，可以拟合出式(6)，并解得a。而后，边缘宽度可计算为[11]

此外，采用图像平均梯度来评价重建图像[13]。平均梯度能够反映出图像细微反差的程度，其值越大，表明影像越清晰，反差越好，细节保持得越好，可定义为

式中，f(i，j)、▽if(i，j)和▽jf(i，j)分别为像素点灰度及其在行、列方向上的梯度，M和N分别为图像的行数和列数。

2 结果

2.1 运动估计实验结果

图5显示了某组数据中典型的两不同相位间冠状面和矢状面图像运动估计结果，可以明显看出肺的整体运动趋势。同时，由于心脏运动和膈肌运动的影响，在靠近心脏和膈肌的地方，存在一些横向收缩等运动。

图6是对运动估计准确性进行评价的实验结果。从原始实验数据中，任意选取两个相位的冠状面图像(a)和(b);利用完全搜索块匹配运动估计算法，得出两图像之间的运动矢量场;根据该运动矢量场，对图(a)进行运动补偿后得到图像(c)。可以看出，模拟图像(c)与原始图像(b)相似度很高。图像(d)为原图(b)与模拟图(c)之间的差值图像，可看出两图像误差较小，说明了全搜索块匹配运动估计方法能较准确地估计出不同相位图像间的运动场，有助于下一步迭代反投影(IBP)超分辨率重建算法的实现。

图5 冠矢状面运动矢量场估计。(a)冠状面运动矢量场估计结果;(b)矢状面运动矢量场估计结果Fig.5 Motion estimation of coronal and sagittal views.(a)Motion vector field of coronal images;(b)Motion vector field of sagittal images

图6 运动估计结果。(a)和(b)不同相位的两幅图像;(c)对(a)进行运动补偿后的模拟图像;(d)图像(c)与(b)之间的误差Fig.6 The result of motion estimation.(a)and(b)The images of different phases;(c)The motion compensation image from(a);(d)The difference of(c)and(b)

2.2 超分辨率重建实验结果

2.2.1 量化评价

根据本文1.4节定义的边缘宽度和图像平均梯度两个量化指标，分别针对10组数据重建的图像进行了量化评价。对某原始低分辨率图像采用最近邻插值、双线性插值和本算法，分别重建肺部图像，如图7所示。在相同边界采集相同数量连续的像素点，把所取得的边界像素点灰度在图7(a)中分别标记出来，并且根据式(6)画出相应的拟合曲线。图7(b)为肺部矢状面图，其中的白色矩形框表示所取边缘的位置。从图7(a)的拟合曲线可看出，利用本方法重建图像的边缘处像素灰度下降最快，而最近邻插值和双线性插值的变化相对平缓。这说明，利用本方法重建的图像边界更清晰，分辨率更高。同时，在10组数据中，针对不同方法重建的图像结果上采集了多组边界像素点，并根据式(7)计算出边缘宽度及总平均值，如表1所示。可见，本方法与传统的插值方法(如最近邻插值、双线性插值法)相比，图像边缘宽度分别显著降低(P＜0.001)，可见边界质量有明显提升。

表1 各算法重建图像的边缘宽度比较Tab.1 The edge widths value for image reconstructed by different algorithms

图7 不同图像相同边缘像素点灰度随像素坐标的变化。(a)不同方法重建的图像边缘灰度变化的拟合结果。(b)肺部矢状面图中的白色矩形框表示所取边缘的位置Fig.7 A sample edge from different image and the corresponding fitting curve.(a)The fitting curve corresponding to nearest interpolation，linear interpolation and our method，respectively.(b)The white rectangle of lung sagittal view show the taken edge position.

现将表1中3组不同方法的边缘宽度值进行Levene方差齐性检验结果(F=1.047，P=0.365)显示，三组方差齐。对3组数据进行单向方差分析，F=134.943，P＜0.001，按 P=0.05的检验水准，差异有显著性意义，可认为最近邻插值、双线性插值、本方法的边缘宽度值总体均数不全相等。其中，本算法的边缘宽度值最小，双线性插值的边缘宽度值次之，最近邻插值的边缘宽度值最大。Bonferroni多重比较结果见表2，可认为最近邻插值的边缘宽度大于双线性插值的边缘宽度值，最近邻插值的边缘宽度大于本算法的边缘宽度，双线性插值的边缘宽度大于本算法的边缘宽度。

表2 Bonferroni多重比较结果Tab.2 Multiple comparison results of Bonferroni

平均梯度既反映图像中的微小细节反差与纹理变化特征，也反映图像的清晰度。同样是在10组数据中，根据式(8)计算双线性插值和本算法重建图像的平均梯度，如表3所示。可见，本算法较双线性插值重建的图像平均梯度显著提高(P ＜0.001)，图像更清晰，细节保持得更好。

将表3中两种不同方法的平均梯度值进行Levene方差齐性检验结果(F=0.752，P=0.397)显示，两组方差齐。对两组数据进行两独立样本t检验，t=-6.885，P ＜0.001，按 α=0.05的检验水准，差异有显著性意义，可认为双线性插值的平均梯度低于本算法的平均梯度。

表3 重建图像的平均梯度值Tab.3 The average gradient value for image reconstructed by different algorithms

2.2.2 视觉评价

图8显示了某组数据中一相位(相位0)的典型冠矢状面图像重建效果，从左至右分别为最近邻插值、双线性插值和本方法的重建结果，其中(a)和(c)分别为冠状面和矢状面图像。为了更清楚地进行效果对比，分别将黑色方框内的区域进行放大，置于原图下方显示。从实验结果可看出，本算法获得的超分辨率重建结果比传统的最近邻插值、双线性插值方法有较大的改进，图像的分辨率有较明显的提高;从局部放大图像来看，肺实质中的血管和周边组织的亮度及清晰度有明显增强。

3 讨论和结论

图8 实验结果(每行从左至右，分别为最近邻插值、双线性插值及本算法的重建结果)。(a)冠状面;(b)(a)中各图像方框区域的局部放大;(c)矢状面;(d)(c)中各图像方框区域的局部放大Fig.8 The experimental results(from left to right，each row are the reconstruction results of nearest interpolation，bilinear interpolation and our method，respectively).(a)coronal view;(b)The enlarged image in black areas of(a);(c)sagittal view;(d)The enlarged image in black areas of(c)

为了提高肺4D-CT数据层间分辨率，本研究采用基于运动估计的多帧图像超分辨率重建技术。该技术融合了所有原始图像信息，即充分利用了肺部各相位图像的呼吸运动信息，避免了重建过程信息单一的不足，进而重建出优于传统插值方法的高分辨率肺部图像。

在本研究中，运动估计算法的选取尤为重要。为保留原始数据的全部信息，本研究采用块匹配运动估计算法;在估计两个不同相位图像之间运动矢量场的过程中，为保证运动矢量的精度，采用完全搜索块匹配方法，但是计算速度缓慢、时间较长。在图像块搜索过程中，常用图像块大小有8像素×8像素、16像素×16像素和32像素×32像素等。若选取的图像块过大，将会明显增加搜索计算;而如果图像块过小，由于肺纹理比较稀疏，会导致匹配搜索精度不足。因此，本研究设定的图像块大小为16像素×16像素，同时考虑肺的运动幅度不大(像素位移一般为10个像素以内)，故将搜索范围设定为16像素×16像素，也就是在目标像素点周围16像素×16像素的范围内搜索。

本算法的不足之处是运算时间较长，重建图像大约耗时40 min;主要时间消耗在图像运动估计过程中，占总体时间的98%。在实验中，由于要将需重建的特定相位图像与其他9个相位图像通过运动估计分别得到相应的运动矢量Tk，且在迭代过程中要重新进行运动估计，因此导致运算量比较大。在下一步的工作中，拟采用文献[14]所提的快速算法进行运动估计。从文献[14]报道的结果看，其运动估计速度可比完全块搜索提高20余倍，这有望大幅提高本方法的速度。

分辨率提高技术可以克服由采集时间和放射剂量造成的图像分辨率限制。4D-CT在肺癌放射治疗中发挥着重要的作用，它能提供一个全面的呼吸运动信息，有助于跟踪肿瘤运动，实施精确放射治疗，并减少对正常组织的损伤。在本研究中，笔者提出了一个新型的基于运动估计的肺4D-CT图像冠矢状面超分辨率重建算法，根据肺4D-CT图像数据特点，将不同相位的对应冠矢状面图像视为不同的“帧”图像，采用完全搜索块匹配的运动估计算法来估计图像之间的运动场，再利用迭代反投影(IBP)超分辨率重建算法，重建出清晰的肺冠矢状面图像。实验结果表明，笔者提出的超分辨率重建算法在视觉评价方面及定量评价方面，均优于传统的最近邻插值和线性插值法。

[1]Khan F，Bell G，Antony J，et al.The use of 4DCT to reduce lung dose:A dosimetric analysis[J].Medical Dosimetry，2010，34(4):273-278.

[2]Li T，Schreibmann E，Thorndyke B，et al.Radiation dose reduction in four-dimensional computed tomography[J].Medical Physics，2005，32(12):3650-3660.

[3]浦剑，张军平，黄华.超分辨率算法研究综述[J].山东大学学报 (工学版)，2009，39(1):27-32.

[4]De Vos L，Stegherr M.Parameterizable VLSI architectures for the full-search block-matching algorithm[J].IEEE Transactions on Circuits and Systems，1989，36(10):1309-1316.

[5]郭伟伟，章品正.基于迭代反投影的超分辨率图像重建[J].计算机科学与探索，2009，3(3):321-329.

[6]李桐.超分辨率图像重建技术[J].哈尔滨师范大学自然科学学报，2006，5:69-71.

[7]Rousseau F.A non-local approach for image super-resolution using intermodality priors[J].Medical Image Analysis，2010，14(4):594-605.

[8]黄新生，杨庆伟，王亦平，等.图像序列运动估计技术综述[J].计算机仿真，2008，5:180-183.

[9]朱长征，沈振康.一种改进的完全搜索块匹配算法[J].红外与激光工程，2004，33(4):389-391.

[10]Irani M，Peleg S.Motion analysis for image enhancement:Resolution，occlusion，and transparency[J].Journal of Visual Communication and Image Representation，1993，4(4):324-335.

[11]Castillo R，Castillo E，Guerra R，et al.A framework for evaluation of deformable image registration spatial accuracy using large landmark point sets[J].Physics in Medicine and Biology，2009，54(7):1849-1870.

[12]Greenspan H， Oz G， Kiryati N， et al. MRI inter-slice reconstruction using super-resolution[J].Magnetic Resonance Imaging，2002，20(5):437-446.

[13]徐美芳，刘晶红.基于边缘保持的航拍图像凸集投影超分辨率重建算法[J].液晶与显示，2010，25(6):873-878.

[14]Chan SH，Võ DT，Nguyen TQ.Subpixel motion estimation without interpolation[C]//Acoustics Speech and Signal Processing(ICASSP)，Dallas:IEEE，2010:722-725.