CN105844630A - 一种双目视觉的图像超分辨率融合去噪方法 - Google Patents

Abstract

本发明公开了一种双目视觉的图像超分辨率融合去噪方法。其步骤为：(1)输入待超分辨率融合的双目图像；(2)全局位置配准图像；(3)局部位置匹配图像；(4)更新双目图像；(5)图像超分辨率融合；(6)图像去噪；(7)输出最终高分辨率图像。本发明在图像配准中加入局部位置配准方法，并采用拉普拉斯算子应用于图像超分辨率融合，再对融合图像使用非局部均值滤波方法去噪，克服了现有技术无法解决图像局部位置不匹配，细节增强不足，以及不具备抑制噪声的缺陷，得到图像细节增强，噪声减小的超分辨率融合图像。

Description

一种双目视觉的图像超分辨率融合去噪方法

技术领域

本发明属于图像处理技术领域，更进一步涉及图像超分辨率技术领域中的一种双目视觉的图像超分辨率融合去噪方法。本发明可以应用在智能终端设备或专业摄像设备的图像后处理。

背景技术

双目视觉的图像超分辨率融合方法是指从双目摄像机获得的同一场景的左灰度图和右灰度图，利用它们在时空上的相关性及信息上的互补性，采用重建技术进行图像超分辨率融合。重建技术假定超分辨率图像在适当的变形、平移和子采样及噪声干扰下，利用多帧低分辨率图像作为数据一致性约束，并结合图像先验知识进行求解。该方法一般包括两个重要部分：配准和重建。配准是获得多帧低分辨率图像间的亚像素精度的相对运动；重建是利用先验知识，对目标图像进行优化求解。

华南理工大学申请的专利“一种基于SIFT的场景物体实时配准方法”(申请日：2015年9月30日，申请号：CN201510646525.7，公开号：CN105279522A)公开了一种图像配准的方法。该方法是对图像提取SIFT特征，然后利用RANSAC优化算法和轮廓排除法优化匹配结果，得到匹配的特征点对，并采用SIFT算子来提取特征点，然后利用RANSAC优化算法和轮廓排除法来得到优化的匹配结果。该方法能够较好的匹配双目图像的全局位置，但是仍然存在的不足之处是，双目图像的左灰度图和右灰度图的不仅存在全局位置的仿射变换，还存在局部位置的平移、旋转、遮挡等变换，该方法不能解决局部位置的变换问题。

Farsiu S,Robinson M D,Elad M,et al.“Fast and robust multiframe superresolution”([J].Image processing,IEEE Transactions on,2004,13(10):1327-1344.)中提出一种多帧图像的超分辨率方法。该方法采用BTV算子应用于图像超分辨率融合中。该方法对融合后的图像的边缘保持较好，对匹配的精确度有一定鲁棒性。但是，该方法仍然存在的不足之处是，该方法对于图像的细节增强不足，对于噪声较大的图像，不具备去噪能力。

发明内容

本发明的目的在于克服上述已有技术的缺点，提出一种双目图像超分辨率的方法，以解决局部位置无法匹配的问题，同时增强图像细节，减小图像噪声。

为实现上述目的，本发明的技术方案是首先将左灰度图与右灰度图做SIFT全局位置配准，然后再对两图做局部位置匹配，得到完全匹配的左灰度图与右灰度图，接着利用拉普拉斯算子对两图进行超分辨率融合，最后，利用非局部均值滤波方法对融合后的图像进行去噪，得到最终高分辨率图像。

本发明的具体步骤如下：

(1)输入待超分辨率融合的双目图像：

输入待超分辨率融合的的双目图像，其中，左摄像头拍摄的图像为左灰度图，右摄像头拍摄的图像为右灰度图；

(2)全局位置配准图像：

(2a)采用尺度不变特征转换SIFT方法，对左灰度图和右灰度图处理，得到左灰度图特征点和右灰度图的特征点；

(2b)利用欧氏距离，从左灰度图中选取一个最终特征点，并在右灰度图找与该最终特征点最近的两个特征点，计算这两个特征点中距离中最近的距离与次近的距离的比值，判断该比值是否小于0.4，若是，则将选取的左灰度图的最终特征点与这两个特征点中距离最近的特征点匹配，将匹配好的特征点匹配对加入特征点匹配对集合；否则，丢弃选取的左灰度图的最终特征点；

(2c)判断左灰度图中的最终特征点是否选取完，若是，执行步骤(2d)，否则，执行步骤(2b)；

(2d)利用随机抽样一致RANSAC算法，剔除误特征点匹配对，得到特征点匹配对；

(2e)利用得到特征点匹配对的坐标对应关系，列出坐标对应方程组，求解方程组，得到仿射变换矩阵；

(2f)利用仿射变换矩阵，配准左灰度图和右灰度图，得到全局位置配准后的左灰度图和右灰度图；

(3)局部位置匹配图像：

(3a)将左灰度图的像素值拷贝至一个空白的图像集中，得到参考图像；将右灰度图的像素值拷贝至另一个空白的图像集中，得到目标图像；创建一个空白的重构图像集；

(3b)在参考图像和重构图像中，以40个像素为宽度，等距离的将参考图像和重构图像均分为多个正方形的宏块，得到均分宏块后的参考图像和重构图像；选取均分宏块后的参考图像左上角的第一个宏块，作为参考宏块；选取均分宏块后的重构图像左上角的第一个宏块，作为重构宏块；将运动向量初始化为零向量；

(3c)判断参考宏块是否位于参考图像左侧第一列，若是，则将参考宏块的臂长取值为20个像素；否则，按照下式，计算参考宏块的臂长：

Γ＝MAX{C_X,C_Y}

其中，Γ表示参考宏块的臂长，MAX表示取最大值操作，C_X表示运动向量的水平分量，C_Y表示运动向量的垂直分量；

(3d)以参考宏块在参考图像中的坐标，在目标图像的相同坐标处，划定大小相同的目标宏块；

(3e)取目标宏块的上下左右方向、距离为臂长的四个宏块，再取目标宏块用运动向量指向位置的宏块，将这五个宏块和目标宏块，共6个宏块作为检测宏块；

(3f)按照下式，分别计算6个检测宏块与参考宏块的代价值：

D＝Σ|E-O|

其中，D表示检测宏块与参考宏块的代价值，Σ表示求和操作，|·|表示取绝对值操作，E表示参考宏块，O表示检测宏块；

(3g)从6个代价值中找出代价值最小的检测宏块，用代价值最小的检测宏块替换目标宏块，得到十字宏块；

(3h)取十字宏块的上下左右四个方向、距离为1的四个宏块，将这四个宏块和十字宏块，共5个宏块作为比较宏块；

(3i)按照下式，分别计算比较宏块与参考宏块的5个代价值：

G＝∑|N-I|

其中，G表示比较宏块与参考宏块的代价值，∑表示求和操作，|·|表示取绝对值操作，N表示参考宏块，I表示比较宏块；

(3j)从5个代价值中找出代价值最小的比较宏块，判断该代价值最小的比较宏块是否为十字宏块，若是，则执行步骤(3k)，否则，用该代价值最小的比较宏块替换十字宏块，执行步骤(3h)；

(3k)将运动向量更新为十字宏块相对于参考宏块的矢量距离，将十字宏块的像素值拷贝到重构宏块；

(3l)判断参考宏块是否为参考图像的最后一个宏块，若是，则得到重构图像，执行步骤(4)；否则，选取参考图像的下一个宏块，作为参考宏块；选取重构图像的下一个宏块，作为重构宏块，执行步骤(3c)；

(4)更新双目图像：

(4a)将重构图像与目标图像相减，得到差值图像；

(4b)创建一个空白的新重构图像集；

(4c)按照长宽均为40像素，分别对差值图像和新重构图像等距离的划分大小相同的宏块，得到划分宏块后的差值图像；

(4d)从划分宏块后的差值图像中，每次选取一个宏块，判断该宏块的像素值总和是否大于差值阈值1600，若是，则用参考图像对应位置的宏块替换新重构图像对应位置的宏块，否则，按照下式，计算新重构图像对应的宏块：

J＝(1-K)×L+K×M

其中，J表示新重构图像对应的宏块，K表示宏块边界权值矩阵，L表示重构宏块，M表示参考宏块；

(4e)判断选取的宏块是否为差值图像的最后一个宏块，若是，则得到新重构图像，执行步骤(4f)，否则，执行步骤(4d)；

(4f)将参考图像的像素值拷贝至左灰度图，将新重构图像的像素值拷贝至右灰度图，完成左灰度图和右灰度图的更新；

(5)图像超分辨率融合：

(5a)利用光流法，计算左灰度图与右灰度图的估计位移；

(5b)按照下式，计算低分辨率帧间位移和高分辨率帧间位移：

u＝FLO(RND(a*2)/2)

v＝MOD(RND(a*2),2)

其中，u表示低分辨率帧间位移，FLO表示向下取整操作，RND表示四舍五入操作，a表示左灰度图与右灰度图的估计位移，/表示除法操作，v表示高分辨率帧间位移，MOD表示取余操作；

(5c)按照下式，计算新右灰度图：

Gn＝VEC(u,Gr)

其中，Gn表示新右灰度图，VEC表示位移变换操作，u表示低分辨率帧间位移；Gr表示右灰度图；

(5d)将新右灰度图的像素值拷贝到右灰度图中；

(5e)按照下式，计算高分辨率估计图：

Z＝INT(Gl,2)

其中，Z表示高分辨率估计图，INT表示插值操作，Gl表示左灰度图，2表示插值的倍数；

(5f)将超分辨率迭代次数初始化为1；

(5g)按照下式，计算高分辨率估计图的反向图像：

Zr＝DSP(COV(VEC(v,Z),h))

其中，Zr表示高分辨率估计图的反向图像，DSP表示下采样操作，COV表示卷积操作，VEC表示位移变换操作，v表示高分辨率帧间位移，Z表示高分辨率估计图，h表示高分辨率估计图的降质矩阵；

(5h)按照下式，计算高分辨率估计图的反向梯度：

Zg＝VEC(COV(USP(Zr),h^T),v^T)

其中，Zg表示高分辨率估计图的反向梯度，VEC表示位移变换操作，COV表示卷积操作，USP表示上采样操作，^T表示转置操作，Zr表示高分辨率估计图的反向图像，h表示高分辨率估计图的降质矩阵，v表示高分辨率帧间位移；

(5i)按照下式，计算高分辨率估计图的反向正则：

Ze＝COV(COV(Lp,Z),Lp^T)

其中，Ze表示高分辨率估计图的反向正则，COV表示卷积操作，Lp表示拉普拉斯算子，Z表示高分辨率估计图，^T表示转置操作；

(5j)按照下式，计算高分辨率估计图的迭代差值：

Zn＝be*(Zg+lm*Ze)

其中，Zn表示计算高分辨率估计图的迭代差值，be表示高分辨率估计图的迭代系数，Zg表示高分辨率估计图的反向梯度，lm表示高分辨率估计图的正则系数，Ze表示高分辨率估计图的反向正则；

(5k)将高分辨率估计图的像素值减去高分辨率估计图的迭代差值，得到更新后的高分辨率估计图；

(5k)判断超分辨率迭代次数是否大于等于迭代阈值20，若是，则执行步骤(6)，否则，将超分辨率迭代次数加上1，执行步骤(5g)；

(6)图像去噪：

(6a)创建与高分辨率估计图相同大小的空白的去噪图像集，创建值为0的计数矩阵；

(6b)选取高分辨率估计图的一个像素，将该像素作为处理像素；

(6c)在高分辨率估计图中，以处理像素为中心点，划定宽为40的正方形块，将该正方形块作为处理宏块；

(6d)在高分辨率估计图中，以处理像素为中心点，划定宽为200的正方形块，将该正方形块作为处理搜索窗；

(6e)创建宽为40的正方形框，将该正方形框作为移动窗口；

(6f)选取处理搜索窗内的一个宽为40的正方形块，将该正方形块的像素值拷贝到移动窗口；

(6g)按照下式，计算移动窗口与处理宏块的相关权值数组：

$\mathrm{\ω}\[e,f,p,q\]=\exp \{-\frac{{||Wm-Wc||}_2^2}{2{\mathrm{\σ}}^2}\}\×\sqrt{{(e-p)}^2+{(f-q)}^2}\×\mathrm{\χ}$

其中，ω表示在移动窗口与处理宏块的相关权值数组，e表示移动窗口中心点的横坐标值，f表示移动窗口中心点的纵坐标值，p表示处理宏块中心点的横坐标值，q表示处理宏块中心点的纵坐标值，exp表示以自然常数为底的指数运算操作，‖·‖表示取模值操作，Wm表示移动窗口，Wc表示处理宏块，σ表示相关权值的控制系数，表示开方操作，χ表示相关权值的距离系数；

(6h)判断选取的处理搜索窗内的正方形块是否为处理搜索窗内的最后一个正方形块，若是，则执行步骤(6f)，否则，执行步骤(6g)；

(6f)按照下式，计算去噪宏块：

$Wd=\frac{\mathrm{\λ}\·Wc+\Σ(\mathrm{\ω}\[e,f,p,q\]\·Wm)}{\mathrm{\λ}+\Σ\mathrm{\ω}\[e,f,p,q\]}$

其中，Wd表示去噪宏块，λ表示处理宏块的融合权值，Wc表示处理宏块，∑表示求和操作，ω表示在移动窗口与处理宏块的相关权值数组，e表示移动窗口中心点的横坐标值，f表示移动窗口中心点的纵坐标值，p表示处理宏块中心点的横坐标值，q表示处理宏块中心点的纵坐标值，Wm表示移动窗口；

(6g)将去噪宏块累加到去噪图像中，将计数矩阵中对应与去噪宏块位置的值加上1；

(6h)判断处理像素是否为高分辨率估计图的最后一个像素，若是，则执行步骤(6i)，否则，执行步骤(6b)；

(6i)将去噪图像除以计数矩阵，得到最终高分辨率图像；

(7)输出最终高分辨率图像。

本发明与现有的技术相比具有以下优点：

第一，由于本发明对双目图像的左灰度图与右灰度图做全局位置配准后，再对全局位置配准后的双目图像的左灰度图和右灰度图两图做局部位置配准，克服了现有技术无法解决局部位置不匹配的问题，使得本发明可以获得全局与局部都完全匹配的双目图像，有利于后期对图像的超分辨率融合处理。

第二，由于本发明对双目图像进行配准处理之后，采用拉普拉斯算子应用于双目图像的超分辨率融合处理中，克服了现有技术中融合图像细节增强不足的缺陷，使得本发明得到的融合图像分辨率比现有技术的图像细节明显提高。

第三，由于本发明对双目图像进行超分辨率融合处理之后，采用非局部均值滤波方法对融合图像进行去噪，克服了现有技术中超分辨率融合图像不具备抑制噪声的问题，使得本发明比现有技术得到的融合图像信噪比大大提高。

附图说明

图1为本发明的流程图。

图2为本发明的仿真图。

具体实施方式

下面结合附图对本发明做进一步的描述。

参照附图1，对本发明的具体实施方式描述如下。

步骤1，输入待超分辨率融合的双目图像。

输入待超分辨率融合的的双目图像，其中，左摄像头拍摄的图像为左灰度图，右摄像头拍摄的图像为右灰度图。

步骤2，全局位置配准图像；

(2a)采用尺度不变特征转换SIFT方法，对左灰度图和右灰度图处理，得到左灰度图特征点和右灰度图的特征点；

尺度不变特征转换SIFT方法的具体步骤如下：

第1步，对左灰度图和右灰度图，分别用尺度递增的高斯滤波器进行滤波，得到滤波后的图像，将滤波后的图像按照尺度递增顺序每8个划分为一组，组成一个子八度；将所有子八度组成高斯金字塔图层；将高斯金字塔中相邻的图层相减，得到差分图层，将所有差分图层组成差分高斯金字塔；

第2步，在差分高斯金字塔中，逐一选取每个图层中的像素点，将所选取的像素点与选取像素所在图层的其余8个相邻像素点和选取像素点的上下相邻图层的18个相邻像素点的值逐一进行比较，若选取像素的值是所有比较的像素点的最大值或者是最小值，则将该像素点作为候选的特征点；

第3步，去除对噪声敏感的低对比度的候选特征点和具有不稳定的边缘响应的候选特征点，剩余的是最终的特征点；

第4步，按照下式，计算以最终特征点为中心的邻域像素的梯度的模值和方向：

$m(x,y)=\sqrt{{\left(A\right(x+1,y)-A(x-1,y\left)\right)}^2+{\left(A\right(x,y+1)-A(x,y-1\left)\right)}^2)}$

$\mathrm{\θ}(x,y)=arctan\left(\frac{A(x,y+1)-A(x,y-1)}{A(x+1,y)-A(x-1,y)}\right)$

其中，m(x,y)表示以最终特征点为中心的邻域像素的梯度的模值，A(x+1,y)表示以最终特征点为中心的邻域像素的右边像素在高斯空间的值，A(x-1,y)表示以最终特征点为中心的邻域像素的左边像素在高斯空间的值，A(x,y+1)表示以最终特征点为中心的邻域像素的下边像素在高斯空间的值，A(x,y-1)表示以最终特征点为中心的邻域像素的上边像素在高斯空间的值，x表示以最终特征点为中心的邻域像素的纵坐标值，y表示以最终特征点为中心的邻域像素的横坐标值，θ(x,y)表示以最终特征点为中心的邻域像素的梯度的方向，arctan表示反正切操作；

第5步，以最终特征点为中心，选择16×16的邻域，并将该邻域划分为16个4×4的子区域，按照最终特征点为中心的邻域像素的梯度方向，将以最终特征点为中心的邻域像素在每个子区域上按照0°、45°、135°、180°、225°、270°、315°、360°共8个方向划分，并在每个方向将所有以(2b)利用欧氏距离，从左灰度图中选取一个最终特征点，并在右灰度图找与该最终特征点最近的两个特征点，计算这两个特征点中距离中最近的距离与次近的距离的比值，判断该比值是否小于0.4，若是，则将选取的左灰度图的最终特征点与这两个特征点中距离最近的特征点匹配，将匹配好的特征点匹配对加入特征点匹配对集合；否则，丢弃选取的左灰度图的最终特征点；

(2c)判断左灰度图中的最终特征点是否选取完，若是，执行步骤(2d)，否则，执行步骤(2b)；

(2d)利用随机抽样一致RANSAC算法，剔除误特征点匹配对，得到特征点匹配对；

(2e)利用得到特征点匹配对的坐标对应关系，列出坐标对应方程组，求解方程组，得到仿射变换矩阵；

(2f)利用仿射变换矩阵，配准左灰度图和右灰度图，得到全局位置配准后的左灰度图和右灰度图；

(2b)利用欧氏距离，从左灰度图中选取一个最终特征点，并在右灰度图找与该最终特征点最近的两个特征点，计算这两个特征点中距离中最近的距离与次近的距离的比值，判断该比值是否小于0.4，若是，则将选取的左灰度图的最终特征点与这两个特征点中距离最近的特征点匹配，将匹配好的特征点匹配对加入特征点匹配对集合；否则，丢弃选取的左灰度图的最终特征点；

(2c)判断左灰度图中的最终特征点是否选取完，若是，执行步骤(2d)，否则，执行步骤(2b)；

(2d)利用随机抽样一致RANSAC算法，剔除误特征点匹配对，得到特征点匹配对；

(2e)利用得到特征点匹配对的坐标对应关系，列出坐标对应方程组，求解方程组，得到仿射变换矩阵；

(2f)利用仿射变换矩阵，配准左灰度图和右灰度图，得到全局位置配准后的左灰度图和右灰度图。

步骤3，局部位置匹配图像；

(3a)将左灰度图的像素值拷贝至一个空白的图像集中，得到参考图像；将右灰度图的像素值拷贝至另一个空白的图像集中，得到目标图像；创建一个空白的重构图像集；

(3b)在参考图像和重构图像中，以40个像素为宽度，等距离的将参考图像和重构图像均分为多个正方形的宏块，得到均分宏块后的参考图像和重构图像；选取均分宏块后的参考图像左上角的第一个宏块，作为参考宏块；选取均分宏块后的重构图像左上角的第一个宏块，作为重构宏块；将运动向量初始化为零向量；

(3c)判断参考宏块是否位于参考图像左侧第一列，若是，则将参考宏块的臂长取值为20个像素；否则，按照下式，计算参考宏块的臂长：

Γ＝MAX{C_X,C_Y}

其中，Γ表示参考宏块的臂长，MAX表示取最大值操作，C_X表示运动向量的水平分量，C_Y表示运动向量的垂直分量；

(3d)以参考宏块在参考图像中的坐标，在目标图像的相同坐标处，划定大小相同的目标宏块；

(3e)取目标宏块的上下左右方向、距离为臂长的四个宏块，再取目标宏块用运动向量指向位置的宏块，将这五个宏块和目标宏块，共6个宏块作为检测宏块；

(3f)按照下式，分别计算6个检测宏块与参考宏块的代价值：

D＝Σ|E-O|

其中，D表示检测宏块与参考宏块的代价值，Σ表示求和操作，|·|表示取绝对值操作，E表示参考宏块，O表示检测宏块；

(3g)从6个代价值中找出代价值最小的检测宏块，用代价值最小的检测宏块替换目标宏块，得到十字宏块；

(3h)取十字宏块的上下左右四个方向、距离为1的四个宏块，将这四个宏块和十字宏块，共5个宏块作为比较宏块；

(3i)按照下式，分别计算比较宏块与参考宏块的5个代价值：

G＝Σ|N-I|

其中，G表示比较宏块与参考宏块的代价值，∑表示求和操作，|·|表示取绝对值操作，N表示参考宏块，I表示比较宏块；

(3j)从5个代价值中找出代价值最小的比较宏块，判断该代价值最小的比较宏块是否为十字宏块，若是，则执行步骤(3k)，否则，用该代价值最小的比较宏块替换十字宏块，执行步骤(3h)；

(3k)将运动向量更新为十字宏块相对于参考宏块的矢量距离，将十字宏块的像素值拷贝到重构宏块；

(3l)判断参考宏块是否为参考图像的最后一个宏块，若是，则得到重构图像，执行步骤4；否则，选取参考图像的下一个宏块，作为参考宏块；选取重构图像的下一个宏块，作为重构宏块，执行步骤(3c)。

步骤4，更新双目图像；

(4a)将重构图像与目标图像相减，得到差值图像；

(4b)创建一个空白的新重构图像集；

(4c)按照长宽均为40像素，分别对差值图像和新重构图像等距离的划分大小相同的宏块，得到划分宏块后的差值图像；

(4d)从划分宏块后的差值图像中，每次选取一个宏块，判断该宏块的像素值总和是否大于差值阈值1600，若是，则用参考图像对应位置的宏块替换新重构图像对应位置的宏块，否则，按照下式，计算新重构图像对应的宏块：

J＝(1-K)×L+K×M

其中，J表示新重构图像对应的宏块，K表示宏块边界权值矩阵，L表示重构宏块，M表示参考宏块；

(4e)判断选取的宏块是否为差值图像的最后一个宏块，若是，则得到新重构图像，执行步骤(4f)，否则，执行步骤(4d)；

(4f)将参考图像的像素值拷贝至左灰度图，将新重构图像的像素值拷贝至右灰度图，完成左灰度图和右灰度图的更新。

步骤5，图像超分辨率融合；

(5a)利用光流法，计算左灰度图与右灰度图的估计位移；

光流法的具体步骤如下：

第1步：对左灰度图和右灰度图，分别进行第一次下采样，得到左灰度图和右灰度图的第一次下采样图像，并将该图像作为左灰度图和右灰度图的第一层图层；分别对左灰度图和右灰度图的第一次下采样图像进行第二次下采样，得到左灰度图和右灰度图的第二次下采样图像，并将该图像作为左灰度图和右灰度图的第二层图层；分别对左灰度图和右灰度图的第二次下采样图像进行第三次下采样，得到左灰度图和右灰度图的第三次下采样图像，并将该图像作为左灰度图和右灰度图的第三层图层，得到左灰度图和右灰度图各自的三层下采样图层；对左灰度图和右灰度图各自的三层下采样图层，分别进行高斯滤波操作，得到滤波后的图层，将滤波后的图层组成左灰度图和右灰度图的高斯金字塔；

第2步：将左灰度图的高斯金字塔的第三层图层拷贝至的左处理图层，将右灰度图的高斯金字塔的第三层图层拷贝至右处理图层；

第3步：按照下式，计算左处理图层的水平方向梯度和垂直方向梯度：

Iy＝COV(V,Ry)

Ix＝COV(V,Rx)

其中，Iy表示左处理图层的水平方向梯度，COV表示卷积操作，V表示左处理图层，Ry表示sobel水平算子，Ix表示左处理图层的垂直方向梯度，Rx表示sobel垂直算子；

第4步：按照下式，计算左处理图层的偏导向量：

R＝[Ix Iy]

其中，R表示左处理图层的偏导向量，[·]表示组合向量操作，Ix表示左处理图层的垂直方向梯度，Iy表示左处理图层的水平方向梯度；

第5步：按照下式，计算左处理图层的偏导矩阵：

U＝R^T*R

其中，U表示左处理图层的偏导矩阵，R表示左处理图层的偏导向量，^T表示转置操作，*表示矩阵相乘操作；

第6步：将左灰度图和右灰度图的估计位移初始化为0；

第7步：按照下式，更新右处理图层：

Cl＝VEC(a,Cl)

其中，Cl表示右处理图层，VEC表示位移变换操作，a表示左灰度图与右灰度图的估计位移；

第8步：按照下式，更新左灰度图与右灰度图的估计位移：

b＝U\(R^T×(V-W))

其中，b表示左灰度图与右灰度图的估计位移增量，U表示左处理图层的偏导矩阵，\表示矩阵左除操作，R表示左处理图层的偏导向量，^T表示转置操作，V表示左处理图层，W表示右处理图层；

第9步：将左灰度图与右灰度图的估计位移加上左灰度图与右灰度图的估计位移增量；

第10步：判断左灰度图与右灰度图的估计位移增量的模值是否大于终止阈值0.2，若是，则执行第11步，否则，将左灰度图与右灰度图的估计位移乘以2，执行第7步；

第11步：判断左处理图层是否为左灰度图高斯金字塔的第一层图层，若是，则得到左灰度图与右灰度图的估计位移；否则，将左处理图层在左灰度图高斯金字塔的上一层图层拷贝至左处理图层，将右处理图层在右灰度图高斯金字塔的上一层图层拷贝至右处理图层，执行第3步。

(5b)按照下式，计算低分辨率帧间位移和高分辨率帧间位移：

u＝FLO(RND(a*2)/2)

v＝MOD(RND(a*2),2)

其中，u表示低分辨率帧间位移，FLO表示向下取整操作，RND表示四舍五入操作，a表示左灰度图与右灰度图的估计位移，/表示除法操作，v表示高分辨率帧间位移，MOD表示取余操作；

(5c)按照下式，计算新右灰度图：

Gn＝VEC(u,Gr)

其中，Gn表示新右灰度图，VEC表示位移变换操作，u表示低分辨率帧间位移；

Gr表示右灰度图；

(5d)将新右灰度图的像素值拷贝到右灰度图中；

(5e)按照下式，计算高分辨率估计图：

Z＝INT(Gl,2)

其中，Z表示高分辨率估计图，INT表示插值操作，Gl表示左灰度图，2表示插值的倍数；

(5f)将超分辨率迭代次数初始化为1；

(5g)按照下式，计算高分辨率估计图的反向图像：

Zr＝DSP(COV(VEC(v,Z),h))

其中，Zr表示高分辨率估计图的反向图像，DSP表示下采样操作，COV表示卷积操作，VEC表示位移变换操作，v表示高分辨率帧间位移，Z表示高分辨率估计图，h表示高分辨率估计图的降质矩阵；

(5h)按照下式，计算高分辨率估计图的反向梯度：

Zg＝VEC(COV(USP(Zr),h^T),v^T)

其中，Zg表示高分辨率估计图的反向梯度，VEC表示位移变换操作，COV表示卷积操作，USP表示上采样操作，^T表示转置操作，Zr表示高分辨率估计图的反向图像，h表示高分辨率估计图的降质矩阵，v表示高分辨率帧间位移；

(5i)按照下式，计算高分辨率估计图的反向正则：

Ze＝COV(COV(Lp,Z),Lp^T)

其中，Ze表示高分辨率估计图的反向正则，COV表示卷积操作，Lp表示拉普拉斯算子，Z表示高分辨率估计图，^T表示转置操作；

(5j)按照下式，计算高分辨率估计图的迭代差值：

Zn＝be*(Zg+lm*Ze)

其中，Zn表示计算高分辨率估计图的迭代差值，be表示高分辨率估计图的迭代系数，Zg表示高分辨率估计图的反向梯度，lm表示高分辨率估计图的正则系数，Ze表示高分辨率估计图的反向正则；

(5k)将高分辨率估计图的像素值减去高分辨率估计图的迭代差值，得到更新后的高分辨率估计图；

(5k)判断超分辨率迭代次数是否大于等于迭代阈值20，若是，则执行步骤6，否则，将超分辨率迭代次数加上1，执行步骤(5g)。

步骤6，图像去噪；

(6a)创建与高分辨率估计图相同大小的空白的去噪图像集，创建值为0的计数矩阵；

(6b)选取高分辨率估计图的一个像素，将该像素作为处理像素；

(6c)在高分辨率估计图中，以处理像素为中心点，划定宽为40的正方形块，将该正方形块作为处理宏块；

(6d)在高分辨率估计图中，以处理像素为中心点，划定宽为200的正方形块，将该正方形块作为处理搜索窗；

(6e)创建宽为40的正方形框，将该正方形框作为移动窗口；

(6f)选取处理搜索窗内的一个宽为40的正方形块，将该正方形块的像素值拷贝到移动窗口；

(6g)按照下式，计算移动窗口与处理宏块的相关权值数组：

$\mathrm{\ω}\[e,f,p,q\]=\exp \{-\frac{{||Wm-Wc||}_2^2}{2{\mathrm{\σ}}^2}\}\×\sqrt{{(e-p)}^2+{(f-q)}^2}\×\mathrm{\χ}$

其中，ω表示在移动窗口与处理宏块的相关权值数组，e表示移动窗口中心点的横坐标值，f表示移动窗口中心点的纵坐标值，p表示处理宏块中心点的横坐标值，q表示处理宏块中心点的纵坐标值，exp表示以自然常数为底的指数运算操作，‖·‖表示取模值操作，Wm表示移动窗口，Wc表示处理宏块，σ表示相关权值的控制系数，表示开方操作，χ表示相关权值的距离系数；

(6h)判断选取的处理搜索窗内的正方形块是否为处理搜索窗内的最后一个正方形块，若是，则执行步骤(6f)，否则，执行步骤(6g)；

(6f)按照下式，计算去噪宏块：

$Wd=\frac{\mathrm{\λ}\·Wc+\Σ(\mathrm{\ω}\[e,f,p,q\]\·Wm)}{\mathrm{\λ}+\Σ\mathrm{\ω}\[e,f,p,q\]}$

其中，Wd表示去噪宏块，λ表示处理宏块的融合权值，Wc表示处理宏块，∑表示求和操作，ω表示在移动窗口与处理宏块的相关权值数组，e表示移动窗口中心点的横坐标值，f表示移动窗口中心点的纵坐标值，p表示处理宏块中心点的横坐标值，q表示处理宏块中心点的纵坐标值，Wm表示移动窗口；

(6g)将去噪宏块累加到去噪图像中，将计数矩阵中对应与去噪宏块位置的值加上1；

(6h)判断处理像素是否为高分辨率估计图的最后一个像素，若是，则执行步骤(6i)，否则，执行步骤(6b)；

(6i)将去噪图像除以计数矩阵，得到最终高分辨率图像。

步骤7，输出最终高分辨率图像。

下面结合附图2对本发明的仿真效果做进一步说明。

1.仿真条件：

本发明的仿真实验所使用的待处理图像通过固定在支架的两台三星Galaxy S5手机在同一时刻拍摄的，图像分辨率为3984*2988。截取拍摄的两幅图像中显著区域分辨率为600*400的图像，作为双目图像的左灰度图与右灰度图。实验仿真环境为MATLAB 7.10.0(R2010a)Service Pack 3,Intel(R)Pentium(R)4CPU 3.00GHz,WindowXP Professional。

2.仿真结果与分析：

附图2是本发明的仿真结果图，其中，附图2(a)为左灰度图；附图2(b)为右灰度图；附图2(c)为对左灰度图与右灰度图使用双三次插值处理后的效果图；附图2(d)为对左灰度图与右灰度图使用本发明处理后的效果图。

对比附图2中的四幅子图，可以看到，使用双三次插值处理后的效果图中，文字图案的模糊边缘较大，图像噪声明显，图像细节不清晰。而可以看到，使用本发明处理后的效果图中，文字图案的模糊边缘较小，图像噪声明显减少，图像细节清晰度提高，从而有效地解决了上述的问题。由此可以看出，本发明可以有效的增强图像细节，减少图像噪声，提高图像超分辨率的质量，比已有技术有较好的客观评价效果。

Claims (3)

1.一种双目视觉的图像超分辨率融合去噪方法，包括如下步骤：

(1)输入待超分辨率融合的双目图像：

输入待超分辨率融合的的双目图像，其中，左摄像头拍摄的图像为左灰度图，右摄像头拍摄的图像为右灰度图；

(2)全局位置配准图像：

(2a)采用尺度不变特征转换SIFT方法，对左灰度图和右灰度图处理，得到左灰度图特征点和右灰度图的特征点；

(2b)利用欧氏距离，从左灰度图中选取一个最终特征点，并在右灰度图找与该最终特征点最近的两个特征点，计算这两个特征点中距离中最近的距离与次近的距离的比值，判断该比值是否小于0.4，若是，则将选取的左灰度图的最终特征点与这两个特征点中距离最近的特征点匹配，将匹配好的特征点匹配对加入特征点匹配对集合；否则，丢弃选取的左灰度图的最终特征点；

(2c)判断左灰度图中的最终特征点是否选取完，若是，执行步骤(2d)，否则，执行步骤(2b)；

(2d)利用随机抽样一致RANSAC算法，剔除误特征点匹配对，得到特征点匹配对；

(2e)利用得到特征点匹配对的坐标对应关系，列出坐标对应方程组，求解方程组，得到仿射变换矩阵；

(2f)利用仿射变换矩阵，配准左灰度图和右灰度图，得到全局位置配准后的左灰度图和右灰度图；

(3)局部位置匹配图像：

(3a)将左灰度图的像素值拷贝至一个空白的图像集中，得到参考图像；将右灰度图的像素值拷贝至另一个空白的图像集中，得到目标图像；创建一个空白的重构图像集；

(3b)在参考图像和重构图像中，以40个像素为宽度，等距离的将参考图像和重构图像均分为多个正方形的宏块，得到均分宏块后的参考图像和重构图像；选取均分宏块后的参考图像左上角的第一个宏块，作为参考宏块；选取均分宏块后的重构图像左上角的第一个宏块，作为重构宏块；将运动向量初始化为零向量；

(3c)判断参考宏块是否位于参考图像左侧第一列，若是，则将参考宏块的臂长取值为20个像素；否则，按照下式，计算参考宏块的臂长：

Γ＝MAX{C_X,C_Y}

其中，Γ表示参考宏块的臂长，MAX表示取最大值操作，C_X表示运动向量的水平分量，C_Y表示运动向量的垂直分量；

(3d)以参考宏块在参考图像中的坐标，在目标图像的相同坐标处，划定大小相同的目标宏块；

(3e)取目标宏块的上下左右方向、距离为臂长的四个宏块，再取目标宏块用运动向量指向位置的宏块，将这五个宏块和目标宏块，共6个宏块作为检测宏块；

(3f)按照下式，分别计算6个检测宏块与参考宏块的代价值：

D＝∑|E-O|

其中，D表示检测宏块与参考宏块的代价值，∑表示求和操作，|·|表示取绝对值操作，E表示参考宏块，O表示检测宏块；

(3g)从6个代价值中找出代价值最小的检测宏块，用代价值最小的检测宏块替换目标宏块，得到十字宏块；

(3h)取十字宏块的上下左右四个方向、距离为1的四个宏块，将这四个宏块和十字宏块，共5个宏块作为比较宏块；

(3i)按照下式，分别计算比较宏块与参考宏块的5个代价值：

G＝∑|N-I|

其中，G表示比较宏块与参考宏块的代价值，∑表示求和操作，|·|表示取绝对值操作，N表示参考宏块，I表示比较宏块；

(3j)从5个代价值中找出代价值最小的比较宏块，判断该代价值最小的比较宏块是否为十字宏块，若是，则执行步骤(3k)，否则，用该代价值最小的比较宏块替换十字宏块，执行步骤(3h)；

(3k)将运动向量更新为十字宏块相对于参考宏块的矢量距离，将十字宏块的像素值拷贝到重构宏块；

(3l)判断参考宏块是否为参考图像的最后一个宏块，若是，则得到重构图像，执行步骤(4)；否则，选取参考图像的下一个宏块，作为参考宏块；选取重构图像的下一个宏块，作为重构宏块，执行步骤(3c)；

(4)更新双目图像：

(4a)将重构图像与目标图像相减，得到差值图像；

(4b)创建一个空白的新重构图像集；

(4c)按照长宽均为40像素，分别对差值图像和新重构图像等距离的划分大小相同的宏块，得到划分宏块后的差值图像；

(4d)从划分宏块后的差值图像中，每次选取一个宏块，判断该宏块的像素值总和是否大于差值阈值1600，若是，则用参考图像对应位置的宏块替换新重构图像对应位置的宏块，否则，按照下式，计算新重构图像对应的宏块：

J＝(1-K)×L+K×M

其中，J表示新重构图像对应的宏块，K表示宏块边界权值矩阵，L表示重构宏块，M表示参考宏块；

(4e)判断选取的宏块是否为差值图像的最后一个宏块，若是，则得到新重构图像，执行步骤(4f)，否则，执行步骤(4d)；

(4f)将参考图像的像素值拷贝至左灰度图，将新重构图像的像素值拷贝至右灰度图，完成左灰度图和右灰度图的更新；

(5)图像超分辨率融合：

(5a)利用光流法，计算左灰度图与右灰度图的估计位移；

(5b)按照下式，计算低分辨率帧间位移和高分辨率帧间位移：

u＝FLO(RND(a*2)/2)

v＝MOD(RND(a*2),2)

其中，u表示低分辨率帧间位移，FLO表示向下取整操作，RND表示四舍五入操作，a表示左灰度图与右灰度图的估计位移，/表示除法操作，v表示高分辨率帧间位移，MOD表示取余操作；

(5c)按照下式，计算新右灰度图：

Gn＝VEC(u,Gr)

其中，Gn表示新右灰度图，VEC表示位移变换操作，u表示低分辨率帧间位移；Gr表示右灰度图；

(5d)将新右灰度图的像素值拷贝到右灰度图中；

(5e)按照下式，计算高分辨率估计图：

Z＝INT(Gl,2)

其中，Z表示高分辨率估计图，INT表示插值操作，Gl表示左灰度图，2表示插值的倍数；

(5f)将超分辨率迭代次数初始化为1；

(5g)按照下式，计算高分辨率估计图的反向图像：

Zr＝DSP(COV(VEC(v,Z),h))

其中，Zr表示高分辨率估计图的反向图像，DSP表示下采样操作，COV表示卷积操作，VEC表示位移变换操作，v表示高分辨率帧间位移，Z表示高分辨率估计图，h表示高分辨率估计图的降质矩阵；

(5h)按照下式，计算高分辨率估计图的反向梯度：

Zg＝VEC(COV(USP(Zr),h^T),v^T)

其中，Zg表示高分辨率估计图的反向梯度，VEC表示位移变换操作，COV表示卷积操作，USP表示上采样操作，^T表示转置操作，Zr表示高分辨率估计图的反向图像，h表示高分辨率估计图的降质矩阵，v表示高分辨率帧间位移；

(5i)按照下式，计算高分辨率估计图的反向正则：

Ze＝COV(COV(Lp,Z),Lp^T)

其中，Ze表示高分辨率估计图的反向正则，COV表示卷积操作，Lp表示拉普拉斯算子，Z表示高分辨率估计图，^T表示转置操作；

(5j)按照下式，计算高分辨率估计图的迭代差值：

Zn＝be*(Zg+lm*Ze)

其中，Zn表示计算高分辨率估计图的迭代差值，be表示高分辨率估计图的迭代系数，Zg表示高分辨率估计图的反向梯度，lm表示高分辨率估计图的正则系数，Ze表示高分辨率估计图的反向正则；

(5k)将高分辨率估计图的像素值减去高分辨率估计图的迭代差值，得到更新后的高分辨率估计图；

(5k)判断超分辨率迭代次数是否大于等于迭代阈值20，若是，则执行步骤(6)，否则，将超分辨率迭代次数加上1，执行步骤(5g)；

(6)图像去噪：

(6a)创建与高分辨率估计图相同大小的空白的去噪图像集，创建值为0的计数矩阵；

(6b)选取高分辨率估计图的一个像素，将该像素作为处理像素；

(6c)在高分辨率估计图中，以处理像素为中心点，划定宽为40的正方形块，将该正方形块作为处理宏块；

(6d)在高分辨率估计图中，以处理像素为中心点，划定宽为200的正方形块，将该正方形块作为处理搜索窗；

(6e)创建宽为40的正方形框，将该正方形框作为移动窗口；

(6f)选取处理搜索窗内的一个宽为40的正方形块，将该正方形块的像素值拷贝到移动窗口；

(6g)按照下式，计算移动窗口与处理宏块的相关权值数组：

$\mathrm{\ω}\[e,f,p,q\]=\exp \{-\frac{\left|\right|Wm-Wc|{|}_2^2}{2{\mathrm{\σ}}^2}\}\×\sqrt{{(e-p)}^2+{(f-q)}^2}\×\mathrm{\χ}$

其中，ω表示在移动窗口与处理宏块的相关权值数组，e表示移动窗口中心点的横坐标值，f表示移动窗口中心点的纵坐标值，p表示处理宏块中心点的横坐标值，q表示处理宏块中心点的纵坐标值，exp表示以自然常数为底的指数运算操作，||·||表示取模值操作，Wm表示移动窗口，Wc表示处理宏块，σ表示相关权值的控制系数，表示开方操作，χ表示相关权值的距离系数；

(6h)判断选取的处理搜索窗内的正方形块是否为处理搜索窗内的最后一个正方形块，若是，则执行步骤(6f)，否则，执行步骤(6g)；

(6f)按照下式，计算去噪宏块：

$Wd=\frac{\mathrm{\λ}\·Wc+\mathrm{\Σ}(\mathrm{\ω}\[e,f,p,q\]\·Wm)}{\mathrm{\λ}+\mathrm{\Σ}\mathrm{\ω}\[e,f,p,q\]}$

其中，Wd表示去噪宏块，λ表示处理宏块的融合权值，Wc表示处理宏块，∑表示求和操作，ω表示在移动窗口与处理宏块的相关权值数组，e表示移动窗口中心点的横坐标值，f表示移动窗口中心点的纵坐标值，p表示处理宏块中心点的横坐标值，q表示处理宏块中心点的纵坐标值，Wm表示移动窗口；

(6g)将去噪宏块累加到去噪图像中，将计数矩阵中对应与去噪宏块位置的值加上1；

(6h)判断处理像素是否为高分辨率估计图的最后一个像素，若是，则执行步骤(6i)，否则，执行步骤(6b)；

(6i)将去噪图像除以计数矩阵，得到最终高分辨率图像；

(7)输出最终高分辨率图像。

2.根据权利要求1所述的一种双目视觉的图像超分辨率融合去噪方法，其特征在于：步骤(2a)中所述的尺度不变特征转换SIFT方法的具体步骤如下：

第1步，对左灰度图和右灰度图，分别用尺度递增的高斯滤波器进行滤波，得到滤波后的图像，将滤波后的图像按照尺度递增顺序每8个划分为一组，组成一个子八度；将所有子八度组成高斯金字塔图层；将高斯金字塔中相邻的图层相减，得到差分图层，将所有差分图层组成差分高斯金字塔；

第2步，在差分高斯金字塔中，逐一选取每个图层中的像素点，将所选取的像素点与选取像素所在图层的其余8个相邻像素点和选取像素点的上下相邻图层的18个相邻像素点的值逐一进行比较，若选取像素的值是所有比较的像素点的最大值或者是最小值，则将该像素点作为候选的特征点；

第3步，去除对噪声敏感的低对比度的候选特征点和具有不稳定的边缘响应的候选特征点，剩余的作为最终的特征点；

第4步，按照下式，计算以最终特征点为中心的邻域像素的梯度的模值和方向：

$m(x,y)=\sqrt{{\left(A\right(x+1,y)-A(x-1,y\left)\right)}^2+{\left(A\right(x,y+1)-A(x,y-1\left)\right)}^2)}$

$\mathrm{\θ}(x,y)=\arctan \left(\frac{A(x,y+1)-A(x,y-1)}{A(x+1,y)-A(x-1,y)}\right)$

其中，m(x,y)表示以最终特征点为中心的邻域像素的梯度的模值，A(x+1,y)表示以最终特征点为中心的邻域像素的右边像素在高斯空间的值，A(x-1,y)表示以最终特征点为中心的邻域像素的左边像素在高斯空间的值，A(x,y+1)表示以最终特征点为中心的邻域像素的下边像素在高斯空间的值，A(x,y-1)表示以最终特征点为中心的邻域像素的上边像素在高斯空间的值，x表示以最终特征点为中心的邻域像素的纵坐标值，y表示以最终特征点为中心的邻域像素的横坐标值，θ(x,y)表示以最终特征点为中心的邻域像素的梯度的方向，arctan表示反正切操作；

第5步，以最终特征点为中心，选择16×16的邻域，并将该邻域划分为16个4×4的子区域，按照最终特征点为中心的邻域像素的梯度方向，将以最终特征点为中心的邻域像素在每个子区域上按照0°、45°、135°、180°、225°、270°、315°、360°共8个方向划分，并在每个方向将所有以最终特征点为中心的邻域像素的梯度模值累加，最终生成128维的特征向量。

3.根据权利要求1所述的一种双目视觉的图像超分辨率融合去噪方法，其特征在于：步骤(5a)中所述的光流法的具体步骤如下：

第1步：对左灰度图和右灰度图，分别进行第一次下采样，得到左灰度图和右灰度图的第一次下采样图像，并将该图像作为左灰度图和右灰度图的第一层图层；分别对左灰度图和右灰度图的第一次下采样图像进行第二次下采样，得到左灰度图和右灰度图的第二次下采样图像，并将该图像作为左灰度图和右灰度图的第二层图层；分别对左灰度图和右灰度图的第二次下采样图像进行第三次下采样，得到左灰度图和右灰度图的第三次下采样图像，并将该图像作为左灰度图和右灰度图的第三层图层，得到左灰度图和右灰度图各自的三层下采样图层；对左灰度图和右灰度图各自的三层下采样图层，分别进行高斯滤波操作，得到滤波后的图层，将滤波后的图层组成左灰度图和右灰度图的高斯金字塔；

第2步：将左灰度图的高斯金字塔的第三层图层拷贝至的左处理图层，将右灰度图的高斯金字塔的第三层图层拷贝至右处理图层；

第3步：按照下式，计算左处理图层的水平方向梯度值和垂直方向梯度值：

Iy＝COV(V,Ry)

Ix＝COV(V,Rx)

其中，Iy表示左处理图层的水平方向梯度值，COV表示卷积操作，V表示左处理图层，Ry表示sobel水平算子，Ix表示左处理图层的垂直方向梯度值，Rx表示sobel垂直算子；

第4步：按照下式，计算左处理图层的偏导向量：

R＝[Ix Iy]

其中，R表示左处理图层的偏导向量，[·]表示组合向量操作，Ix表示左处理图层的垂直方向梯度值，Iy表示左处理图层的水平方向梯度值；

第5步：按照下式，计算左处理图层的偏导矩阵：

U＝R^T*R

其中，U表示左处理图层的偏导矩阵，R表示左处理图层的偏导向量，^T表示转置操作，*表示矩阵相乘操作；

第6步：将左灰度图和右灰度图的估计位移初始化为0；

第7步：按照下式，更新右处理图层：

Cl＝VEC(a,Cl)

其中，Cl表示右处理图层，VEC表示位移变换操作，a表示左灰度图与右灰度图的估计位移；

第8步：按照下式，更新左灰度图与右灰度图的估计位移：

b＝U\(R^T×(V-W))

其中，b表示左灰度图与右灰度图的估计位移增量，U表示左处理图层的偏导矩阵，\表示矩阵左除操作，R表示左处理图层的偏导向量，^T表示转置操作，V表示左处理图层，W表示右处理图层；

第9步：将左灰度图与右灰度图的估计位移加上左灰度图与右灰度图的估计位移增量；

第10步：判断左灰度图与右灰度图的估计位移增量的模值是否大于终止阈值0.2，若是，则执行第11步，否则，将左灰度图与右灰度图的估计位移乘以2，执行第7步；

第11步：判断左处理图层是否为左灰度图高斯金字塔的第一层图层，若是，则得到左灰度图与右灰度图的估计位移；否则，将左处理图层在左灰度图高斯金字塔的上一层图层拷贝至左处理图层，将右处理图层在右灰度图高斯金字塔的上一层图层拷贝至右处理图层，执行第3步。