肺4D断层扫描（4D-CT）是当今的肺癌放射治疗中必不可少的元素。肺4D-CT图像在肺癌靶向治疗中捕获肺部随呼吸运动的信息，能够指导精确放射治疗。然而，由于存在一些因素，例如人体能承受的放射剂量有限，成像硬件以及时间限制等导致肺4D-CT图像在Z轴方向上的分辨率较低，并且图像常伴有伪影和噪声。而在大多数临床工作中，常需要高分辨率的肺4D-CT图像来引导放射治疗，因此提高肺4D-CT图像的工作十分重要^[1]。

本文的主要目的是通过超分辨率重建技术来提高肺4D-CT图像Z轴的分辨率。超分辨率重建技术常用的方法是结合多幅低分辨率图像提供的不同空间信息来恢复高分辨率图像。例如：迭代反投影算法^[2]，POCS 算法^[3]等，这些算法的结果严重依赖配准结果的精度，并且需要不断地迭代以重建出高分辨率图像。超分辨率重建的另一类方法是基于学习的方法。该方法主要是通过训练图像学习高分辨率图像与低分辨率图像之间的关系，来恢复低分辨率图像的细节信息。例如Yang 等^[4]人提出基于稀疏表达的方法来重建高分辨率图像。 基于学习的方法可以得到较好的结果，但是依然存在一个弊端就是需要高分辨率数据作为训练集，因此基于学习的方法在肺4D-CT数据上并不适用。目前也有一部分超分辨率重建方法是基于MAP-MRF框架来构建超分辨率重建的模型^[5]。例如，Wallach等^[6]提出基于最大后验的超分辨率重建方法来提高正电子断层成像质量。基于MAP-MRF方法的优势在于能较灵活地结合多种形式的正则化项来得到好的重建结果，但是传统的估计最大后验的方法求解常需要结合迭代过程，而且多数情况下会导致非凸集优化问题^[7]。

本文主要工作是提出一种基于图割方法的超分辨率重建方法来提高肺4D-CT图像在Z轴方向上的分辨率。本文的方法分为两个步骤。首先，在MAP-MRF框架下构建肺4D-CT图像超分辨率重建的模型，并构建成全局图模式，然后，用图割方法求解能量函数的最优值。为了提高求解速度和精度，我们采用α-β swap 方法^[8]进行图割优化。实验结果显示，与经典的线性插值和POCS超分辨率重建算法相比，本文的方法在视觉和定量方面都获得了更优的结果。

1 方法 1.1 图割

能量函数优化技术在计算机视觉领域已有长时间发展，图割是图像能量函数优化常用方法之一^[9]。假设 g=（V,E）表示带权图，顶点称为终端，V是顶点的集合，E 是边的集合，图中有两个与其他顶点相连的终端，该图的一个切割就是E的一个子集，使得两个终端之间没有通路连接。图割问题就是要寻找这样的一个切割，使得中包含的边的权值最小。本文将超分辨率重建问题转化成图割优化问题。假设给图中所有的像素点都赋有标签，找到像素点最佳的标签匹配值，使图像能量函数值最小，最佳的标签匹配即为重建出高分辨率图像的像素值。

1.2 图像退化模型

图像退化原理^[10]如图 1所示。

图像退化模型可以用公式（1）表示为：

${{g}_{k}}=D{{H}_{k}}{{M}_{k}}f+{{n}_{k}},$

(1)

其中g_k表示第k 幅图像，f 是要重建出的高分辨率图像。M_k是f 相对于g_k的几何变换矩阵，H_k为光学模糊矩阵，是由光学系统本身、成像系统与原始场景的相对运动以及低分辨率传感器的点扩散函数（point spread function,PSF）形成的；D为降采样矩阵；n_k则表示系统加性噪声。

超分辨率重建过程是图像退化过程的逆问题，其原理就是利用多帧关于同一场景（不同相位）的有相互位移的低分辨率图像，重建高分辨率清晰的图像。在本文中，我们在MAP-MRF框架下构建图像超分辨率重建能量函数，如公式（2）所示：

$E\left( f/g \right)=\sum\limits_{k}{{{\left\| D{{H}_{k}}{{M}_{k}}f-{{g}_{k}} \right\|}^{2}}+\sum\limits_{p,q\in N}{{{V}_{p,q}}\left( {{f}_{p}},{{f}_{q}} \right)}}$

(2)

其中V_p,q( f_p,f_q )为平滑项，N为邻域，f_p，f_q为像素点p，q 的灰度值。

1.3 基于图割的肺4D-CT超分辨率重建 1.3.1 全局图构建

2001年Kolmogorov等^[11]提出图割方法用于能量函数优化后，图割方法就常用于求取多维图像能量函数的全局最小值。考虑到4D-CT数据的特点，一个数据包含多相位低分辨率图像，且每相位图像之间有着相互的联系。因此本文提出了一个全局图割方法来优化针对4D-CT数据建立的图像能量函数。首先将4D-CT图像用一个全局图表示，然后在MAP-MRF 框架下构建多相位图像能量函数，优化求解像素点最佳的标签匹配值。全局图可以用图 2表示。

1.3.2 能量函数构建

图割方法能优化能量函数，前提是函数必须是图的表达形式^[11]，所以我们构建包含全部 相位的全局图形式的能量函数。公式（3）可以写为：

$E\left( {{f}_{k}} \right)=\sum\limits_{p\in S}{\sum\limits_{k=0}^{T}{{{\left( h*{{f}_{k}}\left( p \right)-F_{k}^{B}\left( p \right) \right)}^{2}}+\lambda }}\sum{_{k}^{T}=0\sum{_{p,q\in N}{{V}_{p,q}}\left( {{f}_{k}}\left( p \right),{{f}_{k}}\left( q \right) \right)}}$

(3)

其中f_k为待求解的各相位高分辨率重建图像。S为初始高分辨率图像像素集，T表示相位数，p为高分辨率图像上的像素点，h 表示PSF，λ 为平滑项系数。 V_p,q( f_k (p),f_k (q))表示平滑项，可以保证重建出的每一相位高分辨率图像保持平滑，N表示邻域。我们设定

${{V}_{p,q}}\left( {{f}_{k}}\left( p \right),{{f}_{k}}\left( q \right) \right)=\min \left( \tau ,\left| {{f}_{k}}\left( p \right)-\left. {{f}_{k}}\left( q \right) \right| \right. \right)$

(4)

其中τ为阈值大小。

F^B_k 表示相位K通过三次样条插值得到的初始高分辨率图像fk0，采用Demons配准方法分别向其他相位的低分辨率图像上投影，然后根据投影位置插值重建的高分辨率观察图像^[7]。

为了使式（4）可以采用图割的方法优化求解，我们对h*f_k ( p)进行简化^[12]，卷积h*f_k ( p)可以写成：

$h*{{f}_{k}}\left( p \right)={{w}_{pp}}{{f}_{k}}\left( q \right)+\sum\limits_{r\in {{N}_{r}}}{{{w}_{pr}}{{f}_{k}}\left( r \right)}$

(5)

其中r是像素p的邻域N_r内的点。这里假设邻域像素点对中心点的影响较小，我们设定w_pp值为1，w_pr值为0，则简化后全局图能量函数如下：

$E\left( {{f}_{k}} \right)\sum\limits_{p\in S}{\sum\limits_{k=0}^{T}{{{\left( {{f}_{k}}\left( p \right)-F_{k}^{B}\left( p \right) \right)}^{2}}+\lambda }}\sum\limits_{k=0}^{T}{\sum\limits_{q\in N}{{{V}_{p,q}}}}\left( {{f}_{k}}\left( p \right),{{f}_{k}}\left( q \right) \right)$

(6)

1.3.3 全局图能量函数优化

图割方法有两种最小化算法，swap算法和expansion算法^[8]。本文中我们选择α-β swap算法来优化能量函数（6）。α-β swap算法主要是通过计算数据项$\sum\limits_{p\in S}{\sum\limits_{k=0}^{T}{{{\left( {{f}_{k}}\left( p \right)-F_{k}^{B}\left( p \right) \right)}^{2}}}}$值，如果数据项满足条件，那么像素点的标签由α变成β，反之保持原来的标签，通过这样的策略求解最佳的标签与像素点的配，即最优的高分辨率图像。

1.4 基于全局图割方法的肺4D-CT图像超分辨率重建

本文基于全局图割方法的肺4D-CT图像超分辨率重建步骤如下：

输入：肺4D-CT的低分辨率图像；（1）由肺4D-CT 数据获得不同相位的低分辨率图像序列；（2）将不同相位的低分辨率图像均进行插值，得到初始的高分辨率图像；（3）根据图像退化模型建立能量函数式（2）；（4）选取任意一相位图像插值得到的初始高分辨率图像的像素值集作为所有像素点匹配的初始标签集；（5）计算各相位的投影重建高分辨率观察图像F B K，构建全局图的能量函数式（6）；（6）通过图割方法和α-β swap算法优化求解建立的全局图能量函数式（6）；输出：肺4D-CT的高分辨率图像。

2 结果

本文采用DIR 实验室提供的一个公共可用的肺 4D-CT数据集进行方法评估^[13]。该数据集由10 组肺 4D-CT数据组成。每组肺4D-CT包含10个相位，包括极端的吸气和呼气的末端。对于数据1到数据5，图像大小分别为256´256，层内的体素尺寸范围从0.97 mm´0.97mm~1.16mm´1.16 mm之间。数据6到数据10，图像的大小为512´512，层内的体素尺寸均为0.97 mm´0.97mm。对于所有的数据，层间间距均为2.5mm。

2.1 视觉评价

图 3显示了典型的冠状面和矢状面用不同方法重建的结果。第一行和第三行从左到右分别为冠状面和矢状面由线性插值，POCS算法和本文方法重建结果， 第一行和第三行中矩形里的内容被放大显示在第二行和第四行。我们可以看出，与线性插值和POCS算法结果相比，本文方法重建出来的结果要更加清晰。

2.2 量化评价

我们采用平均梯度^[14]和边缘宽度^[15]作为量化评价指标。

平均梯度能反映图像微小细节反差变化的速率，能表征图像的相对清晰程度。值越大，清晰度就越高，可用公式表示为：

$\nabla f=\frac{1}{\left( M-1 \right)\left( N-1 \right)}\sum\limits_{i=1}^{M-1}{\sum\limits_{j=1}^{N-1}{\sqrt{\frac{\nabla _{i}^{2}f\left( i,j \right)+\nabla _{i}^{2}f\left( i,j \right)}{2}}}}$

(7)

其中，f (i,j)，▽f_i(i,j)和▽f_j(i,j)分别是像素点灰度以及其在行、列方向上的梯度；M和N分别为图像的行数和列数。我们分别对10组肺4D-CT数据用不同的方法重建出来的高分辨率图像进行平均梯度计算。结果如表 1 所示。为了使结果对比的更明显，我们还列出了平均梯度的增加百分比。较线性插值和POCS算法，我们的方法重建出来图像的平均梯度明显提高，说明图像更清晰，细节保持更好。

边缘宽度是衡量空间分辨率的量化参数。边缘宽度主要反映所选边缘‘锐利’程度，其计算公式（8）如下：

$width\left[ edge \right]=\frac{4.4}{a}$

(8)

其中a又双曲线函数$y\left( x \right)=\frac{1}{1+\exp \left( -a\left( x-c \right) \right)}$拟合得到。函数中y(x)为像素点灰度，x为像素点坐标，c为采集的像素点中心。

我们分别对10组数据进行量化评价，在线性插值， POCS算法和文本方法重建出来的图像上局部选取3 个边缘区域A，B，C（图 4）。

根据公式（8），我们对三种方法重建得到高分辨率图像上取3个边缘区域进行边缘宽度计算，结果如下表 2所示。很明显我们方法与线性插值和POCS算法相比较，边缘宽度都有了一定的减小，边缘细节有所增强。

3 结论

高分辨率图像在当今的医疗和科研中起着重要的作用。如可以提高分割精度、改进可视化效果等^[16]。图像处理技术可以提高图像的分辨率，例如基本的插值方法可以改变图像大小。但是这些插值方法会使图像产生一些伪影，而且图像中没有增加额外的信息。超分辨率重建技术可以有效地提高图像分辨率。肺4D-CT图像在肺癌的方式在治疗中发挥着重要的作用，它可以提供肺部随呼吸运动的信息，有利用跟踪肿瘤运动。提高肺4D-CT图像的质量，有利于实施更精确的肿瘤放射治疗。

本文采用了一种基于全局图割方法的超分辨率重建方法来提高肺4D-CT图像的Z轴分辨率。该方法主要特点是以MAP-MRF框架建立一个肺4D-CT多相位高分辨率图像重建的全局能量函数，然后用图割方法和 α-β swap 算法优化能量函数从而重建高分辨率图像。 本方法的优势在于不需要迭代求解，同时由于图割方法本身的特性能获得全局最优解。另外，采用的马尔科夫正则化项，有助于保持图像的平滑，降低重建图像的噪声。实验结果证明，本文方法相比与传统的线性插值和 POCS算法，该方法无论是在视觉上还是在量化评估都优于以上的两种方法。