一种基于四叉树的光流跟踪方法

本发明属于计算机视觉跟踪技术领域，具体涉及一种基于四叉树的光流跟踪方法。

背景技术：

视觉目标跟踪是计算机视觉领域的一个重要研究问题，在军民领域都有着广泛应用的场景，包括安防领域中的车辆跟踪、巡检领域中的机器人导航以及智慧生活领域的人机交互等方面。

在视觉目标跟踪方法中，第一类方法是特征点匹配法。该方法是基于特征点检测和描述子匹配来实现目标跟踪，具有精度高的优势。但该方法随着时间增加，计算复杂度骤增，且易跟丢，导致算法实时性不高、鲁棒性不强，所以有较大的场景约束限制。

在视觉目标跟踪方法中，第二类方法是光流跟踪方法。该方法主要采用基于角点的光流跟踪，具有鲁棒性强、实时性好的优势。但该方法对噪声、光源以及尺度变化过于敏感，且角点的检测存在随机性，导致目标的跟踪精度不高。

技术实现要素：

(1)技术问题

为了满足在大范围真实环境下兼具快速性、准确性、鲁棒性的视觉目标跟踪需求，本发明提出一种基于四叉树的光流跟踪方法。该方法融合了特征点匹配法和光流跟踪方法的各自优势，既可以解决特征点匹配法特征点易跟丢等问题，又可以解决光流法定位精度不高等问题。

(2)技术方案

一种基于四叉树的光流跟踪方法，包括如下步骤：

步骤1：实时接收图像序列，对图像进行自适应直方图均衡化预处理；

步骤2：构建图像的非线性尺度空间,完成对特征点的实时提取；

步骤3：对图像基节点进行分裂，基于四叉树算法，均匀提取各图像块特征点；

步骤4：采用稀疏光流法对特征点进行跟踪；

步骤5：输出图像平移向量。

根据本发明的示例性实施例，所述步骤1中，将图像划分为边角区域、边非角区域和对角区域；如果是边角区域，则采用直方图均衡化处理；如果是边非角区域，则采用线性插值处理；如果是对角区域，则采用双线性插值处理，由此完成对图像的自适应直方图均衡化预处理。

根据本发明的示例性实施例，所述步骤2中，构建图像的非线性尺度空间,完成对特征点的实时提取包括如下步骤：

步骤2.1：通过各向异性扩散滤波构建图像的非线性尺度空间，并利用快速显示扩散算法对图像进行过滤，计算不同尺度下图像的海森矩阵的特征值，求取其极值点，设定判断阈值，即实现对图像中特征点的识别；

步骤2.2：对识别的特征点的响应极值进行降序排列。

根据本发明的示例性实施例，所述步骤3中的对图像基节点进行分裂，基于四叉树算法，均匀提取各图像块特征点包括如下步骤：

步骤3.1：对图像进行根节点划分，定义图像根节点为w/h个，其中w为图像像素宽度，h为图像像素高度，如果w/h计算值为小数，则取整为n；

步骤3.2：一个根节点代表一个图像区域，此时图像被分为n个区域，每个区域大小为w/n*h；

步骤3.3：计算当前根节点区域内的特征点数目，如果根节点区域内特征点数目大于一，则把当前根节点分成四个子节点，子节点的区域是根节点图像区域的四等分，如果根节点区域内特征点数不大于一，则终止根节点分裂；

步骤3.4：将分出来的子节点作为新的根节点重复步骤3.3，而将原先的根节点从列表中删除，不再分裂，如此循环，直到图像中节点数目大于na，或者当前节点无法再继续分裂，则节点分裂结束，最后从每个节点区域选取质量最好的特征点，其中，na为系统需要用来光流跟踪的特征点数目。

根据本发明的示例性实施例，所述步骤4中采用的稀疏光流法是lucas-kanade稀疏光流法。

根据本发明的示例性实施例，步骤2.1中，通过各向异性扩散滤波构建图像的非线性尺度空间可通过非线性偏微分方程来描述:

其中，l是图像亮度矩阵，div与分别表示散度算子与梯度算子，c(·)是传导函数，时间t对应的尺度因子，x,y为像素点在图像帧x,y方向上的坐标；

其中，lσ(·)是图像亮度矩阵l经过高斯平滑后的图像，是图像lσ的梯度算子；

其中，k是控制扩散级别的对比度因子。

根据本发明的示例性实施例，k是梯度图像中亮度累计概率分布为70％所对应灰度值。

根据本发明的示例性实施例，步骤2.1中，利用快速显示扩散算法对图像进行过滤时，扩散步长τj可表示为：

其中，τmax是不破坏迭代过程稳定性的最大步长，n是最小内循环长度，j是第j次快速显示扩散循环。

根据本发明的示例性实施例，步骤2.1中，计算不同尺度下图像的海森矩阵的特征值，求取其极值点，设定判断阈值，即实现对图像中特征点的识别包括：采用海森矩阵求取特征值，判断该特征值是否大于预设的阈值t，通过对比特征值的极值点完成图像中特征点的识别。

根据本发明的示例性实施例，：步骤2.1中，计算不同尺度下图像的海森矩阵的特征值，求取其极值点，设定判断阈值，即实现对图像中特征点的识别时，利用空间尺度的泰勒展开式来求特征点坐标：

l(x)是图像亮度的空间尺度函数，x是特征点坐标，l是已知极值点的图像亮度空间尺度函数。

(3)有益效果

本发明所提出的一种基于四叉树的光流跟踪方法，该方法融合了传统的光流跟踪方法和特征点匹配法的各自优势，既可以克服光流跟踪方法在弱纹理场景下精度不高的缺点，也可以解决特征点匹配法在长时间运行环境下计算复杂度高、鲁棒性低的问题。因此，本发明所提出的方法是一种适用于大范围真实环境的目标跟踪算法，兼具快速性、准确性和鲁棒性的优点。

附图说明

图1为根据本发明实施例的基于四叉树的光流跟踪方法流程图；

图2(a)-(d)为harris角点检测算法、orb算法、brisk算法和本发明的算法特征点检测对比试例图；

图3(a)-(b)为根据本发明实施例的四叉树示意图；

图4(a)-(b)为本发明特征点检测算法以及本发明融合四叉树的特征点检测算法进行了对比试验图；

图5为在公开数据集mh_01_easy上的本发明的方法、vins-mono(关闭回环)方法的对比试验图。

具体实施方式

下面详细描述本发明的实施例，所述实施例是示例性的，旨在用于解释本说明，而不能理解为对本发明的限制。

为了使本发明的目的、技术方案和优点更加清楚明白，将结合具体实施例，并参考附图，对本发明进一步详细说明。

图1为本发明实施例的流程图。如图1所示，本发明所提出的一种基于四叉树的光流跟踪方法，包括以下步骤：

步骤1：实时接收视觉传感器所发送的图像序列，对图像数据进行预处理。

本申请中图像数据也可称作图像信息或者图像序列。

在一实施例中，接收从公开数据集mh_01_easy中的单目图像信息，对图像进行灰度化处理，利用自适应直方图算法对图像进行图像增强处理；

图像灰度化处理用到了opencv中cvtcolor函数，可用下式得到：

gray＝(r+g+b)/3(1)

式中gray是处理后的图像像素点值，r、g、b分别表示图像中像素点的红、绿、蓝三色亮度值。

自适应直方图算法是为了将图像的像素灰度尽可能占有全部的灰度级,并使之分布均匀，以提高图像对比度，使得一些不明显的特征点凸显。对于图像i(x,y)，灰度范围为[0，l]，l通常取值256，统计图像中所有灰度数出现的概率密度函数为：

其中px(i)是像素灰度为i的概率密度函数，ni是图像中像素灰度值为i的个数，n是所有像素点个数总和，l是灰度级数，代表像素点的亮暗差别。

然后计算图像像素灰度的概率分布函数，定义为：

其中p(·)是累计概率分布函数，得到累计概率密度函数后，开始计算均衡化后的灰度值，直方图均衡化算式如下：

其中h(·)是灰度值变换函数，round(·)是四舍五入取整函数。

上述直方图均衡化是指对全局图像进行均衡化，而自适应直方图均衡化是对图像局部区域的直方图进行均衡化，具体操作流程就是将图像分为多个矩形块大小，对于每个矩形块子图，分别进行直方图均衡化处理，为了避免矩阵块之间像素值分布不连续，产生明显的分割线，需要限制对比度，在本方法中采用双线性插值方法，双线性插值公式如下所示：

f(m)＝(1-u)(1-v)*f(i,j)+(1-u)*v*f(i,j+1)+u*(1-v)*f(i+1,j)+u*v*f(i+1,j+1)(5)

其中f(·)是亮度变换函数，m是待求像素点的绝对坐标，(u,v)是待求像素点相对于特征点的相对坐标，(i,j)是左上角的特征点绝对坐标。

将原始图像区域分成三种情况处理，首先需要判断当前图像待处理区域是边角区域、边非角区域和对角区域，边角区域指图像中四个角区域，边非角区域指图像中四条边但不包含角的区域，对角区域指图像中不是边角以及边非角的其他区域。如果是边角区域矩阵块，则用直方图均衡化处理；如果是边非角区域，则用线性插值的方法；如果是对角区域，则用双线性插值的方法，完成对原始图像的自适应直方图均衡化处理，以达到对图像增强的效果。

需要说明的是，线性插值算法属于图像处理领域常用算法，此处不再赘述。

步骤2：构建图像的非线性尺度空间,完成对特征点的实时提取，包括以下步骤：

步骤2.1：通过各向异性扩散滤波构建图像的非线性尺度空间，并利用快速显示扩散算法对图像进行过滤，在此基础上，计算不同尺度下图像的海森(hessian)矩阵的特征值，求取其极值点，设定判断阈值，即实现对图像中特征点的识别；

非线性扩散滤波方法(即各向异性扩散滤波方法)是将图像亮度(l)在不同尺度上的变化视为某种形式的流动函数的散度，可以通过非线性偏微分方程来描述:

式中，l是图像亮度矩阵，div与分别表示散度算子与梯度算子，c(·)是传导函数，时间t对应尺度因子。通过设置合适的传导函数c(x,y,t)，可以使得图像的局部结构特性进行自适应扩散，沿着图像边缘方向扩散的快，而图像跨边缘方向扩散的慢，以此保留图像的局部精度。

在实施例中传导函数的构建方式为：

式中，lσ(·)是图像亮度矩阵l经过高斯平滑后的图像，是图像lσ的梯度算子。函数g(x)表达形式如下：

式中，参数k是控制扩散级别的对比度因子，决定图像平滑过程中被加强或被削弱的边缘信息程度，如果k值越大，保留的边缘信息越少，图像越平滑；反之，则边缘信息越多。

在一实施例中，k是梯度图像中亮度累计概率分布为70％所对应灰度值。

由于非线性偏微分方程并没有解析解，需要通过离散偏微分方程来求取近似解，故采用快速显示扩散(fed)进行方程求解，fed结合了显式和半隐式求解的优点，主要思想是运行m个周期的n阶显式扩散步骤，扩散步长τj如下式所示：

式中，τmax是不破坏迭代过程稳定性的最大步长，每一个fed周期所对应的时长θn为：

其中(9)式中的某些τj可能违反稳定性条件，因此将(6)式用向量矩阵表示，如下式所示：

式中，a(lⁱ)是图像信息的传导矩阵，τ为一个恒定的时间迭代步长，上式经过变换后可转换为：

lⁱ⁺¹＝(i+τa(lⁱ))lⁱ(12)

其中i是单位矩阵，考虑到先验信息估计l^i+1,0＝lⁱ，经过迭代之后可以得到一个n阶变步长的fed周期为：

l^(i+1,j+1)＝(i+τja(lⁱ))l^(i+1,j)j＝0,1,...,n-1(13)

在整个fed周期中a(lⁱ)保持不变，当每个fed周期完成，就需要计算新的矩阵值a(lⁱ)。

在此基础上构建非线性尺度空间，尺度空间构建和经典算法sift算法类似，通过设置图像金字塔层数o和组数s后，将层数和组数通过下式映射到对应的尺度σi：

σi(o,s)＝2^ο+s/sο∈[0,o-1],s∈[0,s-1],i∈[0,m](14)

式中，σi是高斯滤波参数标准差，ο是某一层数，s是某一组数，m是fed循环次数，也等于图像金字塔中层数o与组数s的乘积，即得到滤波后的图像总数。

在本实施例中，最大步长τmax为0.25，图像金字塔层数o为4，组数s为4，m为16。

由于非线性扩散不是线性扩散，需要将高斯滤波参数标准差σi转换为时间参数，则需要通过下式实现转换：

式中ti是时间参数，当结束某一层滤波之后，在图像金字塔中下一层进行降采样处理，随后修改对比度因子k，重复上述操作，一个o层(每层s组)图像金字塔构建完成。

在本实施例中，σ0为1.6。

待非线性尺度空间构建完成后，采用hessian矩阵求取特征值，先判断该特征值是否大于预设定的阈值t，再通过对比特征值的极值点完成局部特征点的提取。hessian矩阵行列式值计算公式如下：

式中，σ为当前层尺度参数σi的整数值；lxx，lyy分别是图像水平和垂直方向的二阶偏导；lxy为交叉偏导数。由于特征点检测在不同尺度下，还需要尺度进行归一化，归一化尺度因子oⁱ指图像金字塔中第i层。

在本实施例中，设定的阈值t为0.001。

计算各像素点hessian矩阵特征值，并在所有3×3像素区域中找到中心点的特征值为极大值的像素点，将该极大值点再与金字塔上下两层对应位置σi×σi个像素点对比，验证其是否为极值点，如果是，则定义该点为特征点。在此基础上，利用空间尺度的泰勒展开式来求特征点坐标：

式中，l(x)是图像亮度的空间尺度函数；x是特征点坐标，其特征点亚像素坐标解为：

由于图像金字塔层数越高，其像素间隔距离相对于底层图像越大，将特征点坐标精确到亚像素级别，可以较好的定位其在更高层金字塔中特征点的位置，可以提高特征点检测精度。

在本实施例中，分别使用harris角点检测算法、orb算法、brisk算法和本发明的算法对特征点检测进行对比试验，如图2中的左上、右上、左下以及右下四个区域所示。从图2中可以看出，本发明的算法检测出来的特征点最多，且分布较为平均，多达1500点。

步骤2.2：在识别特征点的基础上，对特征点的响应极值进行降序排列。

步骤3：对图像基节点进行分裂，基于四叉树算法，均匀提取各图像块特征点，包括以下步骤：

步骤3.1：首先对图像进行根节点划分，具体划分准则是，定义图像根节点为w/h个，其中w为图像像素宽度，h为图像像素高度，由于w/h计算值为小数，取整为n；

在本实施例中，w为752，h为480，n为1。

步骤3.2：一个根节点代表一个图像区域，此时图像被分为n个区域，区域大小为w/n*h，即按图像像素宽度n等分；

步骤3.3：计算当前根节点区域内的特征点数目，如果根节点内特征点数目大于一，把当前根节点分成四个子节点，子节点的区域是根节点图像区域的四等分，分别为左上、右上、左下、右下，如果根节点内特征点数不大于一(即小于或者等于1)，则终止分裂。

在本实施例中，对根节点采用如图3所示的四叉树方法进行分解。即将图3中a、b、c、d区域按照步骤3.3的规则进行分裂。

步骤3.4：将分出来的子节点作为新的根节点继续分裂，而将原先的根节点从列表中删除，不再分裂，如此循环分裂，直到图像中节点数目大于na，或者当前节点无法再继续分裂，则节点分裂结束。最后从每个节点区域选取质量最好的特征点，na为系统需要用来光流跟踪的特征点数目。从每个节点区域选取质量最好的特征点时，可基于特征点的响应极值降序排列中选择一个或者多个特征点。

在本实施例中，na为200。分别使用原始的特征点检测算法以及本发明融合四叉树的特征点检测算法进行了对比试验，试验结果如图4中左图、右图所示。从图4中可以看出，采用四叉树的特征点检测算法后特征点明显稀疏了很多，也没有出现特征点密集的情况。

步骤4：采用稀疏光流法对特征点进行跟踪。

稀疏光流法核心思想是基于相邻像素的空间一致性和亮度恒定，由此得出：

i(x,y,t)＝i(x+dx,y+dy,t+dt)(19)

式中i(x,y,t)表示t时刻图像中坐标为(x,y)像素点光强度，dx,dy,dt分别表示在经过dt时刻后，像素点在图像帧的x,y轴方向的移动距离；上式用泰勒展开后，可变换为：

式中ε可忽略不计，由式(19)(20)联立可得：

令u,v是光流沿着图像x,y轴的速度矢量

由式(21)两边同除dt,可得到最终光流法数学表达式：

ixu+iyv+it＝0(23)

式中ix，iy分别表示图像中的像素点灰度值沿x,y方向的图像梯度，it表示图像中的像素点灰度值沿时间参数t的时间梯度。

为求解光流方程中u,ν参数，在本实施例中通过构建5×5图像区域，可得到下式：

式中，pi指图像中的像素点，即

变形求解可得到下式：

其中：

使用最小二乘算法去拟合，可以得到光流方程中的u,v参数，即能得出图像的位移变化，实现视觉跟踪。

步骤5：输出图像平移向量。

为了验证本发明的正确性以及合理性，将本发明整体移植到vins-mono框架下，与原vins-mono算法进行对比。本实施例中采用开源数据集euroc中mh_01_easy对本算法进行测试，实验平台：cpu为i5-7300hq，主频为2.5ghz，内存为16gb，gpu为nvidiageforcegtx1050ti，搭载ubuntu16.04。

euroc数据集搭载了双目相机、惯性测量单元、激光雷达，以及维肯动作捕捉系统(用于计算真值)，其中mh_01_easy是拍摄于室内环境，总长度80.626米，时长187秒，相机采样频率为20hz，imu频率为200hz，激光雷达为20hz,维肯动作捕捉系统可以提供六自由度的位姿测量，频率为100hz。

在本实施例中，对视觉惯性里程计系统进行测试，采用绝对位姿估计(ape)作为评价指标，ape误差是计算估计的位姿和维肯动作捕捉系统中得到的真值位姿之间距离得到。将本发明的方法、vins-mono(关闭回环)方法进行了对比试验，如图5所示。其中黑虚线为动作捕捉的真值轨迹，右上角第一实线为本方法实验轨迹，右上角第二实线为vins-mono(关闭回环)的实验轨迹。为了更好地比较试验结果，两种方法的ape结果如下表所示：

表中max,min,mean以及rmse分别代表定位误差中的最大值、最小值、平均值以及均方根误差，其中均方根误差是评价定位系统精确性的关键指标。

通过上表对比发现，本发明的ape中最大定位误差为0.386241m,相对于vins-mono减小了6.69％；最小定位误差为0.009391m，相对于vins-mono减小了67.62％；平均定位误差为0.183887m，相对于vins-mono减小了10.16％；均方根误差为0.210039m，相比于vins-mono，定位精度提高了6.08％。

由上述数据可以表明本发明在四个定位误差指标上，都有不同程度的提高。可见，基于本发明的视觉惯性里程计，相对于vins-mono算法，提高了定位精度，可以较好的满足机器人导航需求。

需要注意的是，具体实施方式仅仅是对本发明技术方案的解释和说明，不能以此限定权利保护范围。凡根据本发明权利要求书和说明书所做的仅仅是局部改变的，仍应落入本发明的保护范围。