一种基于变分光流估计肺部4D‑CT图像超分辨率处理方法与流程

本发明涉及肺部图像处理领域，尤其涉及一种基于变分光流估计肺部4d-ct图像超分辨率处理方法。

背景技术：

放射治疗是医学上治疗肺癌最常用的手段。肺部4d-ct图像因其能够提供必要的呼吸运动信息来引导医生进行精确地放射治疗，从而得到越来越多的关注。肺部4d-ct图像是由呼吸运动周期中的多个相位各自对应的3d-ct图像按照获取时间的先后顺序排序得到的，能够包含呼吸周期内肺部器官以及肿瘤区域的空间结构，以及运动信息。这些信息可以辅助医生精确定位靶区，有助于病人的放射治疗。但是，众所周知ct图像的获取伴随着高剂量照射，为了降低对病人的辐射量，往往只能降低沿纵向的采样，从而导致肺部4d-ct图像的层间分辨率远低于层内分辨率。在观察3d肺部图像的冠状图时，为了还原真实肺部形状，通常会将3d图像沿z轴方向进行插值放大。传统插值算法会使图像变模糊，因此，如何提高4d-ct图像超分辨率重建图像质量是亟待解决的问题。

技术实现要素：

为了克服上述现有技术中的不足，本发明提供一种基于变分光流估计肺部4d-ct图像超分辨率处理方法，处理方法包括：

s1：采用clg变分光流估计算法得到浮动图像，浮动图像与其余相位图像形成光流场；

s2：结合光流场以及图像放大倍数，得到图像间的仿射变换；

s3：采用改进的非局部迭代反投影算法重建出高分辨率图像。

优选地，步骤s1包括：

基于光流场基本理论，将肺部4d-ct图像看作肺部图像i随时间发生形变而生成的图像序列；即呼吸运动过程中肺部图像各个像素位置的亮度沿着运动轨迹保持不变，即：

i(x(t),y(t),t)＝c(1)

对光流的l1范数约束和各项异性扩散约束，以及对数据项的双边滤波操作，得到肺部4d-ct图像的clg变分光流估计模型，如下式所示：

其中λ为数据残差项与数据平滑项之间的权重系数，bfw表示对数据项进行双边滤波，u＝(u,v)是待求解的光流场；模型引入了各向异性扩散因子以实现光流的各向异性扩散；f(u)为配准后图像i1与参考帧图像i0的差，定义如下：

其中x＝(x,y)为像素点坐标，u0是光流场的初始估计。

优选地，步骤s2结合光流场以及图像放大倍数，得到图像间的仿射变换具体包括：

定义ψi(x)为张量x在第i维度上的融合算子，光流模型可转化成如下等价：

其中，vi,1≤i≤2是松弛变量。

优选地，步骤s3采用改进的非局部迭代反投影算法重建出高分辨率图像具体包括：

设置高分辨率图像的初始估计对初始估计按照图像退化模型生成与集合相对应的低分辨率图像估计集通过反向投影重建误差{i-i⁽⁰⁾}到高分辨率估计来提高超分辨率重建效果；迭代此过程以最小化误差函数；

第n次迭代中，低分辨率图像估计的生成可通过下面的模型模拟实现：

其中，是指图像的仿射变换；h是点扩散函数；↓s为下采样算子；反投影重建过程可表示为：

式中，是tk的逆；p是反投影算子，p的取值将影响算法迭代时间；↑s表示上采样算子；

对ibp算法中的投影重建误差进行非局部均值滤波，保留图像中的高频细节；

利用数据加权平均法进行融合，根据低分辨率图像光流估计误差确定该图像的重建误差在反投影过程中的权重，抑制光流估计误差对重建结果的影响，反投影过程可表示为：

将ωk设置为ik与间欧式距离dk的指数函数，即：

其中t是控制核衰变速度的参数；

将图像所包含的所有相位图像作为已知的低分辨率图像集合，从中选出某一相位图像作为高分辨率图像。

优选地，为了求解式(4)中的条件极值问题，将该约束问题转换成无约束的增广拉格朗日函数，如下式所示：

应用admm算法得到下面的迭代求解步骤：

step1：固定vi和wi，更新u^k+1：

转换为：

step2：固定u和wi，求解vi,i＝1,2

其中，{vi}表示集合{v1,v2}；计算vi,i＝1,2，且等价写成下式

式中的式(9)中分解成规模更小的、并行求解的lasso问题，即

ν是张量vi所包含的一维张量，t是γi中与ν相对应的一维张量，k是ν的维数，di是d中与ν相对应一维张量的第i个元素，式(10)是一维全变分问题求解，并行计算的引入使得vi的计算可以在毫秒时间内完成；

step3：计算wi,i＝1,2

该迭代求解过程在原始残差和对偶残差满足特定条件时结束。

从以上技术方案可以看出，本发明具有以下优点：

基于变分光流估计肺部4d-ct图像超分辨率处理方法构建了一个用于求解肺部4d-ct不同相位图像之间的光流场的变分光流模型；然后，利用快速交替方向乘子法求解该模型，得到不同相位图像之间的光流场；最后，基于光流场，并利用非局部迭代反投影超分辨率重建算法，实现了高分辨率肺部图像的重建.实验结果表明，与已有算法相比，本方法在增强图像纹理结构的同时更好地保留了图像的轮廓。

为了提高重建图像的质量和速度，采用迭代反投影方法，充分利用4d-ct不同相位间的互补信息实现快速的ct图像超分辨率重建，并在以下两个方面进行了创新：其一，鉴于影响超分辨率重建算法效果的主要因素是图像配准的精度，采用了局部和全局相结合clg的变分光流估计模型获得更加精准稠密的光流场，并采用快速交替方向乘子法实现光流模型的快速求解；其二，采用基于边缘的非局部迭代反投影算法，有效增强超分辨率图像中的边缘和纹理细。

附图说明

为了更清楚地说明本发明的技术方案，下面将对描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为滤波器大小对光流精度的影响示意图；

图2为clg变分光流估计算法得到的光流场和根据光流场对基准图像运动补偿后得到补偿图像以及补偿图像与其对应的浮动图像间的差图；

图3为不同算法得到的冠状图超分辨率图像视觉效果对比.前三行图片分别对应以第二份4d-ct数据中第一相位、第五相位以及第八相位为基准图像得到的超分辨率结果，第四行图像依次为第三行图像中红框区域放大后的图像；

图4为不同算法得到的矢状面超分辨率图像视觉效果对比，图片排列方式与图3相一致。

具体实施方式

为使得本发明的发明目的、特征、优点能够更加的明显和易懂，下面将运用具体的实施例及附图，对本发明保护的技术方案进行清楚、完整地描述，显然，下面所描述的实施例仅仅是本发明一部分实施例，而非全部的实施例。基于本专利中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其它实施例，都属于本专利保护的范围。

在本发明中，clg变分光流估计模型是用于计算视频中相邻帧间光流场的模型，4d-ct数据中不同相位间光流场亦可用此模型求解。本发明针对肺部4d-ct图像的clg变分光流估计模型，以及利用快速交替方向乘子算法求解该光流模型的过程。

光流估计是利用图像序列中参考帧与当前帧之间像素值数据的时域相关性来确定当前帧中各像素位置相对于参考帧的“运动”。光流计算模型都是基于光流基本等式，通过引入光流约束条件来实现光流求解的。clg变分光流估计是一种对噪声具有鲁棒性且更加准确的光流计算模型。该模型中引入了对光流的l1范数约束，从而保证了模型对噪声的鲁棒性；模型引入了对数据项的双边滤波操作，以减小相邻像素之间的相互影响；该模型基于图像局部特性实现了光流的各向异性传播，有效地抑制光流过传播带来的负影响。本发明构建适用于肺部4d-ct数据不同相位图像间光流估计的计算模型。

本发明基于光流场基本理论，将肺部4d-ct图像看作肺部图像i随时间发生形变而生成的图像序列，也称为图像一致性假设，即呼吸运动过程中肺部图像各个像素位置的亮度沿着运动轨迹保持不变，即：

i(x(t),y(t),t)＝c(1)

引入对光流的l1范数约束和各项异性扩散约束，以及对数据项的双边滤波操作，得到肺部4d-ct图像的clg变分光流估计模型，如下式所示：

其中λ为数据残差项与数据平滑项之间的权重系数，bfw表示对数据项进行双边滤波，u＝(u,v)是待求解的光流场；模型引入了各向异性扩散因子以实现光流的各向异性扩散；f(u)为配准后图像i1(即除参考帧外的所有帧图像)与参考帧图像i0的差，定义如下：

其中x＝(x,y)为像素点坐标，u0是光流场的初始估计。数据残差项的局部约束使得邻域像素的光流会影响中心像素点光流的计算，为了确保在得到稠密光流场的同时弱化邻域像素的影响，引入对数据残差项的滤波操作.滤波算法有双边滤波以及基于奇异值分解的滤波算法，其中双边滤波算法实现简单且能满足需求，故本发明采用双边滤波对数据残差项进行约束。此外，在clg变分光流估计模型的求解过程中，光流场的各向同性传播会造成光流的过传播，从而降低光流场的精度。本发明中clg变分光流模型通过加入光流的各向异性滤波操作，有效地抑制了光流的过传播。

在本发明中，为了实现3d-ct序列图像间光流场的快速求解，采用快速交替方向乘子法(alternatingdirectionmethodofmultipliers,admm)算法求解本文中的变分光流模型。admm是一种求解优化问题的计算框架，通过将全局问题分解成多个相互独立、较小且易于求解的局部子问题，之后协调子问题的解得到全局问题的解.admm算法结合了增广拉格朗日方法以及求解凸优化问题的对偶分解策略的优势，实现了对凸优化问题的快速求解.定义ψi(x)为张量x在第i维度上的融合算子，光流模型可转化成如下等价问题：

其中，vi,1≤i≤2是松弛变量.为了求解式(4)中的条件极值问题，将该约束问题转换成无约束的增广拉格朗日函数，如下式所示：

应用admm算法得到下面的迭代求解步骤：

step1：固定vi和wi，更新u^k+1：

最优解为：

step2：固定u和wi，求解vi,i＝1,2

其中，{vi}表示集合{v1,v2}。式(8)是可分解的，独立地计算vi,i＝1,2，且可等价写成下式

式中的式(9)中的问题可以分解成规模更小的、可并行求解的lasso问题，即

ν是张量vi所包含的一维张量，t是γi中与ν相对应的一维张量，k是ν的维数，di是d中与ν相对应一维张量的第i个元素。式(10)是一维全变分问题有很多算法可以快速求解，并行计算的引入使得vi的计算可以在毫秒时间内完成。

step3：计算wi,i＝1,2

该迭代求解过程在原始残差和对偶残差满足特定条件时结束。

由于迭代反投影重建过程中并未考虑边缘的方向和强度信息，导致最终收敛于边缘模糊的高分辨率图像。为了提高图像中边缘区域的质量，非局部迭代反投影算法提出借助非局部冗余信息减小反投影重建误差的单幅图像超分辨率重建算法。

在本发明中，通过迭代地最小化重建误差实现多幅低分辨率图像序列重建得到高分辨率图像的迭代反投影算法。对于输入的低分辨率图像集合首先给出高分辨率图像的初始估计然后对初始估计按照图像退化模型生成与集合相对应的低分辨率图像估计集通过反向投影重建误差{i-i⁽⁰⁾}到高分辨率估计来提高超分辨率重建效果.迭代此过程以最小化误差函数。

第n次迭代中，低分辨率图像估计的生成可通过下面的模型模拟实现:

其中，是指图像的仿射变换；h是点扩散函数；↓s为下采样算子，反投影重建过程可表示为：

式中，是tk的逆；p是反投影算子，p的取值将影响算法迭代时间；↑s表示上采样算子。

在本发明中，为了更好地利用不同帧图像间的互补信息，结合迭代反投影算法和非局部迭代反投影算法提出了改进的非局部迭代反投影，实现4d-ct图像超分辨率重建。本发明算法在以下两个方面具有创新性：一方面，对ibp算法中的投影重建误差进行非局部均值滤波，来保留图像中的高频细节；另一方面，对反投影重建过程中的图像融合策略进行了改进。

ibp算法在反投影重建过程中采用了灰度平均的像素级图像数据融合策略，该策略简单且便于实现，但在输入低分辨率图像间光流场估计误差较大的情况下会严重影响超分辨率重建结果的质量。基于此，本发明利用数据加权平均法进行融合，根据低分辨率图像光流估计误差确定该图像的重建误差在反投影过程中的权重。这样便可以很好地抑制光流估计误差对重建结果的影响，从而得到高质量的高分辨率图像.修改后的反投影过程可表示为：

将ωk设置为ik与间欧式距离dk的指数函数，即：

其中t是控制核衰变速度的参数。

本发明将4d-ct图像所包含的所有相位图像作为已知的低分辨率图像集合，从中选出某一相位图像作为高分辨率图像。

在本发明中，通过将多维变分问题转换成标准的l1范数问题，并证明了转换后的问题满足收敛的约束条件，且时间复杂度为o(1/k)，而光流求解算法仅是该算法中的特例(二维变分问题)，所以收敛且时间复杂度为o(1/2)。

当待配准图像与基准间存在较大位移时，式(1)图像灰度一致性假设会造成估计结果误差较大，为此，在求解过程中，利用由粗到细的图像金字塔方法提高光流估计精度。式(2)中各项异性扩散因子d的参数设置为典型值α＝5,β＝1/2，权重系数λ＝1000。增广拉格朗日函数中的惩罚项参数设置为ρ＝10，由光流求解过程可知，双边滤波器的邻域大小的选取对光流估计结果影响最大。因此，采用几种不同大小滤波器并分别对两幅图像间存在大位移和小位移的情况进行了对比实验，如图1所示。

图1中，横轴表示滤波器邻域大小，纵轴表示图像插值误差.插值误差是一种光流估计效果评价标准，其值为真实图像与基于光流补偿得到估计图像间的均方差。由图1可知，当滤波器大小为3×3，5×5时，对于大位移图像间的光流估计效果较差；当滤波器大小为9×9,11×11,13×13时，虽然对图像间存在大位移的情况下光流估计结果较精确，但计算时间较长且对图像间存在小位移情况估计偏差变大了，这是由于滤波器的邻域变大，使得影响当前中心像素的邻域像素变多，从而影响估计结果；当滤波器大小为7×7时，能够很好地处理大位移光流估计问题，且处理时间也能达到要求。因此，实验中将双边滤波器的大小设置为7×7。

图2显示了使用基于clg变分光流估计算法对不同相位中的三幅冠状面图像进行光流估计的两个实例。这两幅冠状面图像是从同一幅4d-ct图像中选取不同相位上的冠状面图像(a)、(b)和(f)，其中(a)作为参考帧图像，(b)和(f)分别作为图像，使用clg变分光流估计算法分别得到(a)和(b)之间的光流场图(c)以及(a)和(f)之间的光流场图(g)。利用得到的光流场(c)对图(a)进行运动补偿后得到图(d)；利用得到的光流场(g)对图(a)进行运动补偿后得到图(h)。图(e)是图(b)和图(d)之间的差图像，图(i)是图(f)和图(h)之间的差图像。观察可知，图(b)和图(f)中肺部相对于(a)整体向下运动且(a)和(b)之间的位移较小、(a)和(f)之间的位移较大，这与得到光流场一致.结果可见，(b)和(d)以及(f)和(h)的差异很小，说明clg变分光流估计能准确估计图像间的位移场，有助于更好地进行超分辨率重建。

实验所用数据是从可形变图像配准库获取的肺部4d-ct数据集，从中选取了5组4d-ct数据(数据1-数据5).每组数据含有呼吸运动过程中的10个不同相位，包括呼气末端相位和吸气末端相位.数据的层内像素尺寸范围从(0.97×0.97)mm²到(1.116×1.16)mm²，层厚为2.5mm.下面将分别展示clg变分光流估计和最终的超分辨率重建结果。

非局部迭代反投影算法可以很好地抑制高对比度边缘的“振铃”现象。本发明将算法与非局部反投影算法进行了对比，以评估算法在边界保持方面的效果。同时，将基于全搜索的运动估计超分辨算法加入对比，论证光流估计对超分辨重建效果的影响。

图3展示了采用不同算法对数据2进行冠状面图像重建的结果图像.由于数据2中层内像素尺寸为(1.16×1.16)mm²且层厚2.5mm，为了使轴向采样率和层内采样率相同，本发明将数据2冠状图上采样倍数设置为2.15，图3中列图像分别是利用双三次插值(bicubic)算法、基于全搜索运动估计(fullsearch)的迭代反投影重建算法、非局部迭代(nlibp)反投影重建算法以及算法得到的冠状面图像超分辨率结果。前三行图像是分别对应数据2中不同相位(相位5、相位8和相位1)的超分辨率重建结果。为了更好地对比不同算法重建效果，在第四行给出了第三行图像中红框区域细节的放大图。对比前三行图像可以发现，当以不同相位图像作为基准图像时得到的超分辨率结果图像间存在差异，但在图像质量方面相差无几。对比放大后的肺实质周围的血管与组织边缘处的细节可以看出，本发明算法重建出的图像具有更清晰的结构，同时边缘和细节信息也得到加强。

图4则给出采用不同算法对数据2进行矢状面图像重建得到的结果图像。图4中图片排列方式与图3一致.结果表明算法的超分辨率重建效果最好，与其他算法相比，算法在重建图像的边缘处理上做得最好。

为了更好地评价重建结果图像的清晰度，采用平均梯度作为量化标准。图像边缘附近的灰度值变化率较大，这种变化率大小可以用来表示图像的清晰度。图像平均梯度就是用图像多维方向上灰度变化率之和来表征图像的相对清晰度，值越大图像越清晰，可用公式表示：

其中，i(x,y)是图像在该位置上的灰度值；和分别表示图像在x方向和y方向上的梯度；m和n分别是图像的行数以及列数。

表(1)列出了分别利用双三次插值算法、全搜索运动估计、clg光流估计以及算法对5幅4d-ct图像选择不同相位图像为浮动图像，对不同平面进行超分辨率重建，然后计算得到所有重建结果图像平均梯度值的平均值。对比表(1)中数据可知，较前三个方法，方法得到的图像平均梯度更高，图像清晰度明显增强。

表1对以下五组数据，用四种不同算法得到的平均梯度

从图3和图4中的细节放大图中可以看出，全搜索算法由于得到的光流场不够精确，导致重建图像中出现一些噪声。本发明算法利用基于clg变分光流估计算法计算得到比全搜索估计更加精确和稠密的光流场，所以超分辨结果要优于全搜索算法。由于非局部迭代反投影只利用低分辨率图像自身的冗余信息，导致超分辨率结果图像在具有较少信息冗余的器官组织边缘出现了模糊。而本发明算法则是很好地利用不同帧图像间的互补信息，在保留细节信息的同时提高了重建图像的清晰度。对比表1中5份4d-ct数据重建结果图像的平均梯度值可知，本发明算法得到的平均梯度值都是最大的，这表明本发明算法更好地保留了图像的细节信息。

实验用c++语言实现了本发明中光流求解算法，并在配备英特尔i5处理器(主频为3.2ghz)的主机上进行了测试。对于分辨率为256×99的浮动图像和基准图像，光流计算时间大约为300ms，极大地缩短了光流计算时间，得益于更短的光流计算时间，对于10幅图像组成的实验图像序列，使用本发明超分辨率算法可以在30s内得到超分辨结果。

本发明提出了一种基于clg变分光流估计的超分辨率重建算法，用于提高4d-ct冠状面图像分辨率。将不同相位的3d-ct图像看作不同“帧”图像，基于clg变分光流模型，并利用快速admm算法迅速收敛的性质，实现了浮动“帧”图像与基准“帧”图像间光流场的快速计算；然后，基于得到的光流场，改进了非局部迭代反投影算法，并利用不同相位间的互补信息增强重建了图像的结构细节信息.实验结果表明，本发明算法不仅有效增强了图像的纹理结构，而且能够更好地保留图像中的轮廓信息。

对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本发明。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。