基于图像灰度引导的非局部立体像对密集匹配方法与流程

本发明涉及一种影像匹配方法，尤其是涉及一种基于图像灰度引导的非局部立体像对密集匹配方法。

背景技术：

立体像对密集匹配是摄影测量和计算机视觉领域经久不衰的研究热点之一，其主要任务是逐像素地找到两张影像之间的同名点。密集匹配在DEM\DSM生产、数字城市三维重建、图像渲染、无人自动驾驶、机器人导航、增强现实、虚拟现实等方面有着很广泛的应用前景。至今为止，绝大多数的立体像对密集匹配算法都可以总结为4个步骤：1.代价计算；2.代价积聚；3.视差计算；4.视差修正。

尽管密集匹配已经经过数十年的发展，技术日趋成熟，其在代价计算、无纹理区域代价积聚等方面仍有很大的研究空间。当前已经出现了双边滤波局部匹配算法、半全局密集匹配算法(SGM)、基于图割/置信传播的全局密集匹配算法等等，但是无论哪一种密集匹配算法，均不能很好地兼顾匹配精度和匹配鲁棒性。局部匹配算法实现简单，计算速度快，在视差边缘具有较好的匹配精度，但是整体匹配鲁棒性差；半全局匹配算法SGM同样能够实现快速匹配，在纹理丰富区域具有较好的匹配效果，但是在纹理贫乏区域匹配鲁棒性差；全局匹配算法匹配鲁棒性好，但是在视差边缘容易出现过度平滑的问题。因此，需要研究一种新的密集匹配算法，能够兼具匹配精度和匹配鲁棒性，生成密集的三维点云。

技术实现要素：

本发明主要是解决传统密集匹配方法在视差边缘匹配精度差、整体匹配结果不鲁棒的问题，提出了一种基于图像灰度引导的立体像对密集匹配方法，能够采用改进的HOG算子，合理描述同名像素之间的相似性，构建代价矩阵；根据设计的二次函数来约束代价传递，采用八方向迭代方式进行非局部代价积聚，获得稳定的代价积聚结果；最后对初始视差图进行修正，并生成密集的高精度三维点云。

本发明的上述技术问题主要是通过下述技术方案得以解决的：

一种基于图像灰度引导的非局部立体像对密集匹配方法，其特征在于，包括以下步骤：

步骤1，建立代价矩阵：采用改进HOG算子，作为描述同名点之间相似性的测度，建立代价矩阵，其中，改进的HOG算子基于以下定义：

定义在一个窗口范围内(一般取5x5或者7x7像素大小的窗口)，立体像对之间的辐射畸变呈线性变化，即满足以下公式：

g_r(p_r)＝c·g_l(p_l)+t (1)

式中，p_l、p_r分别表示左右影像同名点；当立体像对为核线影像时，满足关系式p_rx＝p_l-d，其中d表示视差；g_l、g_r分别表示左右影像同名点的灰度；c、t是线性方程的系数，c表示比例因子，t表示偏移因子；

对立体像对同时计算梯度，可以在局部小窗口范围内消除偏移因子t；可采用Sobel算子来计算影像梯度：

Gx_r(p_r)＝c·Gx_l(p_l)Gy_r(p_r)＝c·Gy_l(p_l) (2)

式中，Gx_l(p_l)、Gx_r(p_r)分别表示点p_l、p_r在水平方向的灰度梯度；Gy_l(p_l)、Gy_r(p_r)分别表示点p_l、p_r在垂直方向的灰度梯度；

为了进一步消除比例因子c，根据式(2)，计算梯度的方向角：

式中，θ_l(p_l)、θ_r(p_r)分别表示点p_l、p_r的梯度方向角；按照象限的不同，值域是[0,360°)；由式(3)可见，梯度方向具有良好的抗线性辐射畸变的能力；

以影像中任意一个像素p为中心，开辟一个WxW大小的窗口，作为一个基本的HOG描述单元(cell)；将梯度方向的值域范围划分为12个区间，每个区间对应30°的范围；将每个区间的初始计数均设为0；统计单元内所有像素的梯度方向角；判断单元内像素的梯度方向落在哪个区间，则该区间对应的计数加1；依次遍历单元内所有的像素，最后，将每个区间的计数除以窗口内像素的总数；根据每个区间的统计情况，构造一个12维的向量，作为描述像素p的特征描述子，如式(4)所示；

V_HOG(p)＝(b₀,b₁,b₂,b₃,b₄,b₅,b₆,b₇,b₈,b₉,b₁₀,b₁₁)^T (4)

式中，b_i(i＝0～11)表示每个区间的值；V_HOG(p)表示像素p的特征描述算子；

步骤2，基于图像灰度引导的非局部代价积聚：根据图像的灰度信息，设计一个二次函数来约束代价传递，采用八方向迭代的方式进行非局部代价积聚，建立可靠的代价积聚矩阵；

步骤3，三维点云生成：根据代价积聚矩阵，生成初始视差图，采用左右一致性检测的方法剔除粗差，修正视差图，生成稠密的高精度三维点云。

在上述的一种基于图像灰度引导的非局部立体像对密集匹配方法，所述的步骤1中建立代价矩阵的具体方法如下：将同名点HOG特征描述向量之间的距离作为代价，基于如下公式；

式中，epl(p_l,d)表示点p_l,在视差为d情况下的同名点，即p_r＝epl(p_l,d)；C_HOG(p_l,d)表示点p_l、p_r之间用HOG计算的代价；表示左影像点p_l的HOG特征描述符；表示右影像点p_r的特征描述符；

在上述的一种基于图像灰度引导的非局部立体像对密集匹配方法，所述的步骤2中，基于图像灰度引导的非局部代价积聚方式：

定义影像上灰度差绝对值在15个像素以内，且中心像素八邻域内的像素，其视差应该是满足一致性条件；根据假设条件，列出如式(6)所示的代价传递方程：

式中，L(p,d)表示像素p在当前路径对应视差d的累积代价；r表示路径的方向；C_HOG(p,d)表示当前像素p对应视差d的HOG代价；p–1表示在当前路径上，像素p的前一个像素；P₁、P₂表示惩罚项，在灰度同质区域内，当相邻像素视差平滑变化时，给予微小的惩罚P₁；当相邻像素视差发生阶跃性的变化，给予较大的惩罚P₂；这样做的好处是：在尽量遵守灰度同质区域邻近像素视差一致的原则下，灵活地允许像素视差发生变化；由于灰度同质区域视差平滑变化的情况比较常见，所以设定的参数P₂要远大于P₁；参数T表示像素p和像素p–1灰度/颜色的接近程度，用于约束邻域像素之间代价的传递。

在上述的一种基于图像灰度引导的非局部立体像对密集匹配方法，所述建立代价矩阵步骤中，用于约束邻域像素之间代价的传递采用一个二次函数来表示；具体是基于二次函数的核函数，在[0,2σ]范围内，保证斜率是递增的，如式(7)所示；

式中，T_G表示采用高斯函数计算出来的T值；a表示二次项系数；σ是平滑因子。

在上述的一种基于图像灰度引导的非局部立体像对密集匹配方法，所述建立代价矩阵步骤中，代价积聚结果与代价累积路径密切相关，具体基于8方向的代价积聚策略；不仅能够包含水平/垂直方向，还包括了4个对角线的方向；算法需要两次迭代，在第一次迭代，需要对8个方向的代价累积结果进行相加，如式(8)所示；

式中，L表示根据式(6)算出来的代价累积结果；r表示扫描线方向，包括0°方向、45°方向、90°方向、135°方向、180°方向、225°方向、270°方向以及315°方向；S表示8个方向扫描线总的代价累积结果；

经过第一次迭代，对每个像素来说，仅仅在扫描线方向上存在路径，扫描线之外的像素是无路径相连的；因此，将第一次迭代的累积结果S作为新的代价，参与第二次迭代计算，有：

式中，表示第二次迭代中，r方向扫描线上的代价累积结果；

最后，对8个方向扫描线的代价累积结果相加，如式(10)所示；经过两次迭代，影像上的每个像素，与整幅影像其它所有的像素之间，都存在路径，或者连通，或者非连通，从而实现了非局部代价累积，增强代价累积的可靠性；

式中，S²表示第二次迭代中，8个方向总的代价累积结果；5根据权利要求2所述的一种基于图像灰度引导的非局部立体像对密集匹配方法，其特征在于，所述的步骤3中，三维点云生成方式如下：

根据步骤2获得的代价累积矩阵，采用WTA策略，计算每个像素的视差，如式(11)所：

其中，d(p)表示像素p在视差图上的视差；表示在视差范围内，最小S²所对应的视差d；采用WTA策略计算影像上所有像素的视差后，即可获得初始视差图；

以左影像为参考影像，可以获得左影像的初始视差图；以右影像作为参考影像，可以获得右影像的初始视差图，采用左右一致性检测的方法，可以剔除误匹配点，获得修正的视差图，如式(12)所示：

式中，d(p_l)表示像素p在左影像的视差；d(p_r)表示像素p在右影像的视差；Invalid表示无效视差；如果同一个像素在左影像和右影像的视差相差大于1个像素，则认为是误匹配点，予以剔除；

根据修正的视差图，可以得到每个有效像素的同名点，根据同名点坐标和影像的内外方位元素，可以生成对应的三维点云，如式(13)所示：

式中，X、Y、Z表示物方点的三维坐标；x_i、y_i(i＝1,2)表示同名像点坐标；f表示焦距；R_i(i＝1,2)表示旋转矩阵；X_si、Y_si、Z_si(i＝1,2)表示外方位线元素。

因此，本发明具有如下优点：在影像之间存在辐射畸变的情况下，仍然能得到精确可靠的匹配结果；采用图像灰度引导的代价传递方式，能够保证视差边缘的匹配精度；采用八方向迭代式的非局部代价积聚方式，能够保证纹理贫乏区域的匹配鲁棒性，最终生成密集的高精度三维点云。本发明在航天摄影测量、航空摄影测量、低空摄影测量和近景摄影测量、数字城市三维重建、无人车自动驾驶等领域有较好的应用前景。

附图说明

图1a为代价计算窗口；

图1b为梯度方向直方图；

图2为八方向迭代式非局部代价积聚路径。

具体实施方式

下面通过实施例，并结合附图，对本发明的技术方案作进一步具体的说明。

实施例：

本发明提供的技术方案是，根据核线立体像对，采用改进HOG算子计算代价，建立代价矩阵；根据图像灰度信息，设计二次函数来约束代价的传递，并采用八方向迭代式的非局部代价积聚策略，获得可靠的代价积聚结果；采用WTA策略计算视差，根据左右一致性检测剔除误匹配，生成稠密的高精度三维点云，包括以下步骤：

步骤1.构建HOG代价矩阵。

代价是描述同名点之间相似性的测度。一个好的代价，应该能精确描述像素及其邻域的灰度特征。本发明将HOG算子作为密集匹配的代价。如果对影像上每个像素都计算完整的HOG特征，那么计算量将非常巨大。同时，HOG算子也不具备抗线性辐射畸变的能力。因此，本发明对HOG算子进行了改进，不仅能够加快代价计算速度，而且使得改进的HOG特征具有抗线性辐射畸变的能力。

影像的辐射情况受到传感器本身、时相、光照条件等因素的影响。实际生产中，立体像对之间存在复杂的辐射畸变。本文假设在很小的窗口范围内，立体像对之间的辐射畸变呈线性变化，即满足以下公式：

g_r(p_r)＝c·g_l(p_l)+t

式中，p_l、p_r分别表示左右影像同名点。当立体像对为核线影像时，满足关系式p_rx＝p_l-d，其中d表示视差。g_l、g_r分别表示左右影像同名点的灰度。c、t是线性方程的系数，c表示比例因子，t表示偏移因子。

对立体像对同时计算梯度，可以在局部小窗口范围内消除偏移因子t。可采用Sobel算子来计算影像梯度：

Gx_r(p_r)＝c·Gx_l(p_l)Gy_r(p_r)＝c·Gy_l(p_l)

式中，Gx_l(p_l)、Gx_r(p_r)分别表示点p_l、p_r在水平方向的灰度梯度；Gy_l(p_l)、Gy_r(p_r)分别表示点p_l、p_r在垂直方向的灰度梯度。

首先，确定左影像参考影像，采用Sobel算子，分别计算左右两张影像的水平梯度图和垂直梯度图。根据左右影像的梯度图，采用下式，计算梯度方向图，获得每个像素的梯度方向：

θ_r(p_r)＝θ_l(p_l)

θ_r(p_r)＝arctan(Gy_r(p_r)/Gx_r(p_r))

θ_l(p_l)＝arctan(Gy_l(p_l)/Gx_l(p_l))

式中，p_l表示左影像的像素；p_r表示右影像的像素；θ表示梯度方向；Gx表示水平方向梯度；Gy表示垂直方向梯度。

我们以影像中任意一个像素p为中心，开辟一个WxW(W一般取5-9)大小的窗口，作为一个基本的HOG描述单元(cell)。将梯度方向的值域范围划分为12个区间，每个区间对应30°的范围。将每个区间的初始计数均设为0。统计单元内所有像素的梯度方向角。判断单元内像素的梯度方向落在哪个区间，则该区间对应的计数加1，如图1所示。图1(a)表示描述单元，箭头方向代表每个像素的梯度方向；像素的背景色与区间是一一对应的。图1(b)表示单元内的区间，区间内的数字表示每个区间的计数。依次遍历单元内所有的像素，最后，将每个区间的计数除以窗口内像素的总数。根据每个区间的统计情况，构造一个12维的向量，作为描述像素p的特征描述子，如下式所示。

V_HOG(p)＝(b₀,b₁,b₂,b₃,b₄,b₅,b₆,b₇,b₈,b₉,b₁₀,b₁₁)^T

式中，b_i(i＝0～11)表示每个区间的值；V_HOG(p)表示像素p的特征描述算子。

与传统HOG算法相比，本发明在构建HOG描述子的过程中，忽略了梯度模的信息，从而使得构造的特征描述子，具有线性辐射畸变不变性。此外，传统的HOG算法需要将若干个基本描述单元(Cell)组合成一个块(Block)，来构成更大的特征描述符。由于密集匹配的代价积聚环节能够将邻域像素的代价积聚在一起，因此为了避免重复计算，本发明对像素特征的描述仅停留在单元(Cell)层面上。与传统HOG算法相比，无疑大大加快了计算速度。

左影像像素p_l可获得对应视差d的HOG特征描述子；右影像像素p_r同样可获得对应视差d的HOG特征描述子。本发明将同名点HOG特征描述向量之间的距离作为代价测度，如下式所示。

式中，epl(p_l,d)表示点p_l,在视差为d情况下的同名点，即p_r＝epl(p_l,d)；C_HOG(p_l,d)表示点p_l、p_r之间用HOG计算的代价；表示左影像点p_l的HOG特征描述符；表示右影像点p_r的特征描述符。

根据上式，可计算左影像上所有像素p_l对应视差d的HOG代价。假设影像宽为W，高为H，视差范围是DR，则需要构建一个WxHxDR大小的代价矩阵，来存储左影像上所有像素的HOG代价。

步骤2.八方向迭代式的非局部代价积聚。

在构建代价矩阵后，沿着扫描线方向进行代价传递。扫描线的方向包括水平方向、垂直方向和对角线方向，如图2中绿色扫描线所示。本发明假设：影像上灰度接近，且距离较近的像素，其视差应该是尽量满足一致性条件。根据假设条件，代价传递方程如下式所示：

式中，L(p,d)表示像素p在当前路径对应视差d的累积代价；r表示路径的方向；C_HOG(p,d)表示当前像素p对应视差d的HOG代价；p–1表示在当前路径上，像素p的前一个像素；P₁、P₂表示惩罚项，一般设定的参数P₂要远大于P₁；参数T表示像素p和像素p–1灰度/颜色的接近程度，用于约束邻域像素之间代价的传递，本发明采用一个二次函数来表示，

传统的密集匹配算法采用高斯核函数来计算参数T。高斯核函数是一个减函数，其斜率也是由大到小逐次递减。这意味着当灰度差在0附近产生微小的变化时，高斯函数值T的下降幅度是最大的！在影像上，即使是灰度同质区域，区域内部像素的灰度也不可能完全一样。我们希望在灰度差不等于0，但仍然较小的情况下，得到较大的T值，从而能够增强同质区域内的代价传递。传统算法取较大的平滑因子σ(σ＝15-25)来解决这个问题，但是在非同质区域像素间代价的传递时，较大的σ也会计算出较大的权值。造成这个问题的关键在于高斯核函数的斜率是递减的。因此，本发明设计了一种基于二次函数的核函数，在[0,2σ]范围内，保证斜率是递增的，如下式所示：

Δg＝|g(p+1)-g(p)|a＝(e^-2-1)/4σ²

式中，T_G表示采用高斯函数计算出来的T值；a表示二次项系数；σ是平滑因子，一般取5。

代价积聚结果与代价累积路径密切相关。本发明定义从起点到终点上所有像素均属于同一区域的路径为连通路径；否则为非连通路径。一个好的代价累积路径应该尽可能的覆盖整个同质区域，即灰度同质区域内像素之间的路径是连通的。大多数匹配算法采用水平/垂直扫描线的方法来定义路径。但是仅仅采用水平方向和垂直方向的扫描线，容易造成灰度同质区域内像素之间的路径不连通。

单单一条扫描线的代价积聚结果是不稳定的，很容易产生“条纹效应”。为了解决这个问题，本发明提出了一种基于8方向的非局部代价积聚策略。不仅能够包含水平/垂直方向，还包括了4个对角线的方向。将八个方向的代价积聚结果进行累加，以增强代价积聚的鲁棒性。算法需要两次迭代，在第一次迭代，需要对8个方向的代价累积结果进行相加，如下式所示。

式中，L表示单一扫描线方向计算出来的代价累积结果；r表示扫描线方向，包括0°方向、45°方向、90°方向、135°方向、180°方向、225°方向、270°方向以及315°方向；S表示8个方向扫描线总的代价累积结果。

经过第一次迭代，对每个像素来说，仅仅在扫描线方向上存在路径，扫描线之外的像素是无路径相连的。即每个像素仅仅受到扫描线方向像素的影响，不受扫描线之外像素的影响。这种代价积聚方式本质上是局部的，代价积聚结果并不鲁棒，尤其在纹理贫乏区域的代价积聚。因此，为了进一步增强代价积聚的鲁棒性，将第一次迭代的累积结果S作为新的代价，参与第二次迭代计算，每条扫描线分别进行代价传递，有：

式中，表示第二次迭代中，r方向扫描线上的代价累积结果。

最后，对8个方向扫描线的代价累积结果相加，如下式所示。经过两次迭代，影像上的每个像素，与整幅影像其它所有的像素之间，都存在路径，从而实现了非局部代价累积，增强代价累积的可靠性，如图2所示。图2表示把代价从像素p₁传递到像素p₃的路径，圆形表示影像像素。虚线直线表示像素p₁的8条扫描线方向；实线直线表示像素p₃的8条扫描线方向。这两组扫描线之间的交点，定义了代价传递的路径。交点用黑色的圆圈表示。箭头表示代价传递的方向。从图2可以看出，基于8方向的非局部代价累积策略，能够为p₁到p₃的代价传递，提供多条路径，有水平方向的路径、垂直方向的路径以及对角线方向的路径。

式中，S²表示第二次迭代中，8个方向总的非局部代价累积结果。

每个像素p对应视差d，有一个代价累积结果S²。遍历整张影像，会存在非常多的代价累积结果。假设影像宽为W，高为H，视差范围是DR，则需要构建一个WxHxDR大小的累积矩阵，来存储左影像上所有像素的累积代价S²。

步骤3.三维点云生成。

在视差范围内，参考影像上每个像素都会对应多个累积代价。采用WTA(Winner Takes ALL)策略，将多个累积代价中最小代价对应的视差作为该像素的最终视差，如下式所示：

其中，d(p)表示像素p在视差图上的视差；表示在视差范围内，最小S²所对应的视差d。

采用WTA策略计算影像上所有像素的视差后，即可获得初始视差图。由于遮挡、影像辐射等原因，不可避免地存在误匹配问题，需要通过一定的手段来检测和剔除误匹配点。本发明采用左右一致性的方法来检测和剔除误匹配。首先，以左影像为参考影像，可以获得左影像的初始视差图；然后，以右影像作为参考影像，按照步骤1和步骤2，同样可以获得右影像的初始视差图。如果匹配点为正确点，则左影像像素和右影像同名像素之间的视差应该是一致的；否则，应该是不一致的。因此，采用左右一致性检测的方法，可以剔除误匹配点，获得修正的视差图，如下式所示：

式中，d(p_l)表示像素p在左影像的视差；d(p_r)表示像素p在右影像的视差；Invalid表示无效视差。如果同一个像素在左影像和右影像的视差相差大于1个像素，则认为是误匹配点，予以剔除。

根据修正的视差图，可以得到每个有效像素的同名点，如下式所示：

x_r＝x_l-d y_r＝y_l

式中，(x_l,y_l)、(x_r,y_r)表示影像的同名点坐标。

根据同名点坐标和影像的内外方位元素，可以生成对应的三维点云，如下式所示：

B_u＝X_S2-X_S1B_v＝Y_S2-Y_S1B_w＝Z_S2-Z_S1

X＝X_S1+U₁＝X_S2+U₂

Y＝Y_S1+V₁＝Y_S2+V₂

Z＝Z_S1+W₁＝Z_S2+W₂

式中，X、Y、Z表示物方点的三维坐标；f表示焦距；R_i(i＝1,2)表示旋转矩阵；X_si、Y_si、Z_si(i＝1,2)表示外方位线元素。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。