CN110390685B - 一种基于事件相机的特征点跟踪方法 - Google Patents

一种基于事件相机的特征点跟踪方法 Download PDF

Info

Publication number
CN110390685B
CN110390685B CN201910672162.2A CN201910672162A CN110390685B CN 110390685 B CN110390685 B CN 110390685B CN 201910672162 A CN201910672162 A CN 201910672162A CN 110390685 B CN110390685 B CN 110390685B
Authority
CN
China
Prior art keywords
event
module
time
template
feature points
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201910672162.2A
Other languages
English (en)
Other versions
CN110390685A (zh
Inventor
史殿习
李凯月
李睿豪
伽晗
王明坤
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
National University of Defense Technology
Original Assignee
National University of Defense Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by National University of Defense Technology filed Critical National University of Defense Technology
Priority to CN201910672162.2A priority Critical patent/CN110390685B/zh
Publication of CN110390685A publication Critical patent/CN110390685A/zh
Application granted granted Critical
Publication of CN110390685B publication Critical patent/CN110390685B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C21/00—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
    • G01C21/10—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration
    • G01C21/12—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning
    • G01C21/16—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning by integrating acceleration or speed, i.e. inertial navigation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/20—Analysis of motion
    • G06T7/246—Analysis of motion using feature-based methods, e.g. the tracking of corners or segments
    • G06T7/248—Analysis of motion using feature-based methods, e.g. the tracking of corners or segments involving reference images or patches
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/20—Analysis of motion
    • G06T7/269—Analysis of motion using gradient-based methods
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/70—Determining position or orientation of objects or cameras
    • G06T7/73—Determining position or orientation of objects or cameras using feature-based methods
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10016—Video; Image sequence

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Theoretical Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Multimedia (AREA)
  • Automation & Control Theory (AREA)
  • Image Analysis (AREA)

Abstract

本发明公开了一种基于事件相机的特征点跟踪方法，目的是提高特征点跟踪精度。技术方案是构建由数据获取模块、初始化模块、事件集选取模块、匹配模块、特征点更新模块、模板边更新模块组成的基于事件相机的特征点跟踪***。初始化模块从图像帧中提取特征点和边图；事件集选取模块在特征点周围从事件流中选取特征点的事件集S；匹配模块将S与特征点周围的模板边进行匹配，求出tk时刻n个特征点的光流集合Gk，特征点更新模块根据Gk计算tk+1时刻n个特征点的位置集合FDk+1，模板边更新模块使用IMU数据对PBDk中的位置进行更新，得到tk+1时刻n个特征点对应的模板边的位置集合PBDk+1。采用本发明既可提高事件流上跟踪特征点的精度，又可延长特征点平均跟踪时间。

Description

一种基于事件相机的特征点跟踪方法
技术领域
本发明涉及计算机图像处理领域，具体涉及运用事件相机完成对图像中特征点的跟踪。
背景技术
SLAM(全称同时定位与建图)作为机器人领域一个重要分支，近年来得到了大家的广泛研究。SLAM试图解决这样的问题：一个机器人在未知的环境中运动，如何通过对环境的观测确定自身的运动轨迹，同时构建出环境的地图。SLAM技术正是为了实现这个目标涉及到的诸多技术的总和。一个完整的SLAM***主要包括前端视觉里程计部分和后端优化部分。视觉里程计部分估计机器人状态，估计方法主要分为两种方法：特征点法和直接法。特征点法是目前机器人状态估计的主流方法，即首先在图像中提取特征点，对不同帧之间的特征点进行匹配，进而对匹配的特征点对做相关操作，估算相机的位姿。常用的点特征包括Harris角点、SIFT、SURF、ORB、HOG特征。不同于特征点法，直接法可以省去提取特征点的过程，直接利用图像中的灰度信息来估计机器人状态，但该方法尚不成熟，鲁棒性较差。
然而，无论是特征点法或直接法，使用标准相机在应对极端环境时仍存在精度和鲁棒性的问题。极端情况主要包括两种：相机高速运动的情况，采用标准相机获取图像时，若相机移动过快，则获取到的图像会出现运动模糊现象；高动态范围场景的情况，场景中光强变化强烈，前后帧明暗变化明显。在这些极端情况下，使用标准相机会严重影响直接法和特征点法的算法效果。此外，标准相机无法提供帧与帧之间更精确的特征点运动轨迹。并且标准相机在静态场景中会产生冗余信息，这不仅导致存储资源的浪费，而且会在后续的图像处理过程中消耗大量额外的计算资源。
一种生物启发的事件相机的出现，克服了标准相机的上述限制。事件相机的典型代表为DVS(Dynamic Vision Sensor，中文名为动态视觉传感器，2011年由PatrickLichtsteiner等人在期刊IEEE Journal of Solid-State Circuits第43卷，第2期，第566-576页上发表的文章《A 128×128 120dB 15μs Latency Asynchronous TemporalContrast Vision Sensor》即《一个128×128像素，120dB动态范围，15微秒延迟的异步时间对比视觉传感器》提出，由瑞士的iviVation AG公司生产。与标准相机不同，事件相机仅输出像素级别的亮度变化。在像素阵列中，每当强度变化超过阈值时，对应的像素便会独立地产生一个输出，输出被称为“事件”。所以，不同于标准相机，事件相机输出的数据为时空当中的事件流。由于低延迟，事件相机在快速移动场景中具有很大的优势。此外，事件相机还可以避免在缓慢变化的场景中记录冗余信息。2014年，一种新的事件相机——DAVIS(Dynamic and Active-pixel Vision Sensor，中文名为动态和主动式像素传感器)由Christian Brandli等人在期刊IEEE Journal of Solid-State Circuits第49卷，第10期，第2333-2344页上发表的文章《A 240×180130dB 3μs Latency Global ShutterSpatiotemporal Vision Sensor》即《一个240×180像素，130dB动态范围，3微秒延迟的全局快门时空视觉传感器》提出，后由瑞士的iviVation AG公司生产。DAVIS将标准相机和事件相机DVS结合在了一起，可以同时输出图像帧和事件流。
对于使用事件相机进行特征点跟踪，Zhu等人在文章《Event-based featuretracking with probabilistic data association》即《使用概率数据关联实现基于事件的特征点跟踪》(2017年，发表于会议International Conference on Robotics andAutomation(ICRA)，第4465–4470页上)中提出了一种直接在事件流中提取特征点，并使用事件流计算特征点光流，从而实现特征点跟踪的方法。由于该方法只使用事件流计算特征点光流，导致跟踪精度不高。Kueng等人在文章《Low-Latency Visual Odometry usingEvent-based Feature Tracks》即《使用基于事件的特征点跟踪的低延迟视觉里程计》(2016年，发表于会议2016 International Conference on Intelligent Robots andSystems(IROS)，Inspec检索号为16504091)中提出了一种同时结合图像帧和事件流的特征点跟踪方法，该方法使用图像帧中的信息与事件流进行几何配准，从而跟踪特征点。使用该方法，每接收一个事件便更新一次特征点的位置，这便使得计算量变大，并且此方法会引入更多的误差，导致特征点跟踪精度变差，无法实现长时间的跟踪。
因此，现有的使用事件相机进行特征点跟踪的方法仍然具有跟踪精度不高，跟踪时长较短的缺点。
发明内容
本发明要解决的技术问题是提供一种基于事件相机的特征点跟踪方法，使用的事件相机为DAVIS，既提高事件流上跟踪特征点的精度，又延长特征点平均跟踪时间。
为解决此问题，本发明提出了一种基于事件相机的特征点跟踪方法，使用的事件相机为DAVIS，本发明同时使用图像帧中的边与事件流中的事件进行匹配，求特征点的光流，由于同时利用了两种形式的信息，使得计算得到的光流更为准确，从而提高特征点跟踪精度。且本发明引入IMU(Inertial Measurement Unit，中文名为惯性测量单元)来更新边的位置，使得跟踪过程中边的位置更为准确，从而延长特征点平均跟踪时间。
具体技术方案是：
第一步，构建基于事件相机的特征点跟踪***。事件相机的特征点跟踪***由数据获取模块、初始化模块、事件集选取模块、匹配模块、特征点更新模块、模板边更新模块组成。
数据获取模块与初始化模块、事件集选取模块、模板边更新模块相连。数据集获取模块从苏黎世大学的公开事件相机数据集“The Event-Camera Dataset and Simulator”(中文名为“事件相机数据集和模拟器”，该数据集由DAVIS采集，数据集中包括图像帧、事件流和IMU数据)上下载数据，并将获取到的图像帧发送给初始化模块，将事件流发送给事件集选取模块，将IMU数据发送给模板边更新模块。
初始化模块与数据获取模块、事件集选取模块、匹配模块相连。初始化模块从数据获取模块接收图像帧，从图像帧中提取特征点和边图，得到特征点的位置和特征点周围的模板边的位置(边图中特征点周围的边称作对应特征点的模板边)，将特征点的位置发送给事件集选取模块，将特征点周围的模板边的位置发送给匹配模块。
事件集选取模块与数据获取模块、初始化模块、特征点更新模块、匹配模块相连，事件集选取模块从数据获取模块接收事件流，从初始化模块(第一次循环)或特征点更新模块(从第二次循环开始)接收特征点的位置，从特征点更新模块接收特征点的光流(从第二次循环开始)，在每个特征点周围从事件流中选取特征点的事件集，将每个特征点的位置以及对应的事件集发送给匹配模块。
匹配模块与初始化模块、事件集选取模块、特征点更新模块、模板边更新模块相连，匹配模块从事件集选取模块接收特征点的位置以及每个特征点对应的事件集，从初始化模块(第一次循环)或模板边更新模块(从第二次循环开始)接收特征点周围的模板边的位置，将特征点的事件集与特征点周围的模板边进行匹配，求出每个特征点的光流，将每个特征点的位置、特征点周围的模板边的位置和特征点的光流发送给特征点更新模块。
特征点更新模块与匹配模块、模板边更新模块、事件集选取模块相连，特征点更新模块从匹配模块接收特征点的位置、特征点周围的模板边的位置、特征点的光流，利用特征点的光流计算特征点新的位置，将特征点周围的模板边的位置以及特征点新的位置发送给模板边更新模块，将特征点的光流和特征点新的位置发送给事件集选取模块，并将特征点新的位置(即特征点跟踪结果)输出，供用户查看。
模板边更新模块与数据获取模块、特征点更新模块、匹配模块相连，模板边更新模块从特征点更新模块接收特征点周围的模板边的位置和特征点新的位置，从数据获取模块接收IMU数据，通过IMU数据来更新特征点周围的模板边的位置，将模板边更新后的位置发送给匹配模块。
第二步，数据获取模块从事件相机数据集“The Event-Camera Dataset andSimulator”中获取图像帧、事件流和IMU数据。
第三步，采用基于事件相机的特征点跟踪***对从t0时刻开始到tN时刻结束的时间段内数据获取模块得到的事件流中的特征点进行跟踪。跟踪过程中将时间区间[t0,tN]分成一系列的子时间区间[t0,t1],...,[tk,tk+1],...,[tN-1,tN]，第k+1个时间子区间表示为[tk,tk+1]，N表示子时间区间个数，N的大小由事件流时间长度决定，取值范围为N≥1。跟踪过程如下：
3.1，令时间序号k＝0；
3.2，初始化模块进行初始化操作，具体步骤如下：
3.2.1初始化模块采用Harris角点检测法(1988年，由C.G.Harris等人在会议Alvey Vision Conference，第15卷，第50期上发表的文章《A combined corner and edgedetector》即《组合式角点和边缘检测器》提出)从数据获取模块得到的图像帧中提取特征点，将提取到的特征点放到集合FD中，令FD为{f1,...,fi,...,fn}，fi代表检测到的第i个特征点，n为特征点个数。将特征点的位置看做时间的函数，将t0时刻FD中n个特征点的位置放在特征点位置集合FD0中，FD0＝{f1(t0),...,fi(t0),...,fn(t0)}，fi(t0)代表特征点fi在t0时刻的位置，将FD0发送给事件集选取模块。
3.2.2初始化模块使用Canny边检测法(1987年，由John Canny等人在会议Readings in Computer Vision，第184-203页上发表的文章《A computational approachto edge detection》即《边缘检测的计算方法》提出)从数据获取模块得到的图像帧中提取边图，每个图像帧对应一个边图。
3.2.3初始化模块选取FD中n个特征点对应的模板边(由模板边上所有模板点的集合表示)，将FD中n个特征点对应的t0时刻模板边的位置集合PBD0发送给匹配模块。方法如下：
3.2.3.1令i＝1；
3.2.3.2以特征点fi在t0时刻的位置fi(t0)为中心，在fi(t0)周围选取一个矩形区域
Figure BDA0002142074750000041
大小为s×s，即
Figure BDA0002142074750000042
的长和宽均为s，s的范围为20到30个像素。将3.2.2检测到的边图中在
Figure BDA0002142074750000043
内的边作为fi的模板边。模板边上的像素点即为fi的模板点。定义
Figure BDA0002142074750000044
内的模板点的集合
Figure BDA0002142074750000045
Figure BDA0002142074750000051
也即fi对应的模板边，将
Figure BDA0002142074750000052
放到模板边集合PB中：
Figure BDA0002142074750000053
pj为fi的第j个模板点。将pj的位置视作时间的函数，pj(t0)表示pj在t0时刻的位置。
Figure BDA0002142074750000054
表示在t0时刻pj在矩形区域
Figure BDA0002142074750000055
内。mi表示集合
Figure BDA0002142074750000056
中模板点的个数。
3.2.3.3令i＝i+1；
3.2.3.4若i≤n，转3.2.3.2；否则，表示得到了FD中n个特征点对应的模板边构成的模板边集合PB，
Figure BDA0002142074750000057
定义
Figure BDA0002142074750000058
中模板点在t0时刻的位置的集合
Figure BDA0002142074750000059
令PB中n个模板边中模板点在t0时刻的位置集合
Figure BDA00021420747500000510
将PBD0发送给匹配模块，转第四步。
第四步，事件集选取模块从数据获取模块接收事件流，根据不同的k值分别从初始化模块或特征点更新模块接收不同的数据，在n个特征点周围从事件流中选取特征点的事件集S，将FDk和特征点的事件集S发送给匹配模块，方法是：
4.1若k＝0，事件集选取模块从数据获取模块接收事件流，从初始化模块接收tk时刻n个特征点的位置集合FDk，FDk＝{f1(tk),...,fi(tk),...,fn(tk)}(此时即FD0＝{f1(t0),...,fi(t0),...,fn(t0)})，令t1＝t0+1，单位为秒，转4.3；
4.2若k≥1，事件集选取模块从数据获取模块接收事件流，从特征点更新模块接收tk时刻n个特征点的位置集合FDk和tk-1时刻n个特征点的光流集合Gk-1，
Figure BDA00021420747500000511
其中
Figure BDA00021420747500000512
表示tk-1时刻特征点fi的光流，预估子时间区间[tk,tk+1]的大小，根据公式(2)计算tk+1：
Figure BDA00021420747500000513
其中
Figure BDA00021420747500000514
表示在子时间区间[tk-1,tk]上特征点fi的光流，从特征点更新模块得到。
Figure BDA00021420747500000515
表示在子时间区间[tk-1,tk]上所有特征点的平均光流。通过公式(2)计算tk+1的物理意义在于，上一时间区间内的特征点平均运动3个像素需要的时间作为当前区间[tk,tk+1]大小的预估值。转4.3；
4.3在tk时刻n个特征点的位置集合FDk中每个位置的周围选取每个特征点对应的事件集。
4.3.1令i＝1；
4.3.2从事件流中选出符合公式(3)要求的事件，将这些事件放到特征点fi在tk时刻对应的事件集
Figure BDA0002142074750000061
中，并将
Figure BDA0002142074750000062
放到事件集S中：
Figure BDA0002142074750000063
Figure BDA0002142074750000064
表示在一个三维时空区域内的事件的集合，空间范围为特征点fi在tk时刻的位置fi(tk)周围的矩形区域Hfi，时间范围为区间[tk,tk+1]。ed代表
Figure BDA0002142074750000065
中的第d个事件，ed来自于事件流，并且在公式(3)中指定的三维时空区域内，d＝{1,2,...,zi}，zi表示
Figure BDA0002142074750000066
内事件的个数。事件ed的表示形式为
Figure BDA00021420747500000613
xd表示事件ed在像素坐标系中的坐标，
Figure BDA0002142074750000067
表示事件ed发生的时间，
Figure BDA0002142074750000068
表示ed的像素坐标在
Figure BDA0002142074750000069
内。
4.3.3令i＝i+1；
4.3.4若i≤n，转4.3.2；否则，说明得到了tk时刻FD中n个特征点对应的事件集S,
Figure BDA00021420747500000610
并将FDk和S发送给匹配模块，转第五步。
第五步，匹配模块从事件集选取模块接收FDk和S，若k＝0，从初始化模块接收n个特征点对应的模板边的位置集合PBDk，转第六步；若k≥1，从模板边更新模块接收n个特征点对应的模板边的位置集合PBDk，转第六步；
第六步，匹配模块将S与特征点周围的模板边进行匹配，求出tk时刻n个特征点的光流集合Gk，将n个特征点的位置集合FDk、n个特征点对应的模板边的位置集合PBDk和光流集合Gk发送给特征点更新模块。具体方法为：
6.1令i＝1；
6.2对特征点fi构造匹配误差函数：
对于
Figure BDA00021420747500000611
内的事件ed，修正事件ed在tk时刻的位置，公式如下：
Figure BDA00021420747500000612
x′d表示计算得到的事件ed在tk时刻的位置，简称为事件ed运动修正的位置，
Figure BDA0002142074750000071
表示在时间区间[tk,tk+1]上特征点fi的光流。符号·表示点乘，下文中出现该符号的地方均表示点乘，在不引起歧义的情况下，下文中部分公式对该符号做了省略。为了表示方便，定义符号
Figure BDA0002142074750000072
构造匹配误差函数如下：
Figure BDA0002142074750000073
ε表示误差，pj(tk)代表tk时刻模板点pj的位置。rdj表示ed由模板点pj生成的概率，即ed与pj匹配的概率。这里将rdj作为权重，rdj为事件ed与模板点pj相匹配的概率，即rdj越大，事件ed与模板点pj之间在tk时刻的距离差占总误差的比例越大。公式中双竖线表示对双竖线内的向量进行取模操作，下同。rdj的计算公式如下：
Figure BDA0002142074750000074
该公式的物理意义在于，分子为事件ed运动修正的位置与模板点pj距离的平方，分母为事件ed运动修正的位置与特征点fi对应的mi个模板点的距离的平方和，将两者相除作为事件ed与模板点pj相匹配的概率。
6.3采用EM-ICP算法(是一种匹配算法，2002年，由Sébastien Granger等人在会议European Conference on Computer Vision，第418-432页发表的文章《Multi-scale EM-ICP:A fast and robust approach for surface registration》即《多尺度EM-ICP：一种快速而稳健的表面配准方法》提出)。最小化6.2步构造得到的匹配误差函数，求解得到最优的光流
Figure BDA0002142074750000075
求解过程如下：
6.3.1初始化
Figure BDA0002142074750000076
6.3.2通过公式(6)计算rdj；
6.3.3更新光流：
Figure BDA0002142074750000081
6.3.4计算光流的变化量
Figure BDA0002142074750000082
并令
Figure BDA0002142074750000083
将
Figure BDA0002142074750000084
放到tk时刻光流集合Gk中。
6.3.5若
Figure BDA0002142074750000085
表示得到最终的优化结果
Figure BDA0002142074750000086
转6.4；若
Figure BDA0002142074750000087
转6.3.2。σ为光流的变化量阈值，取值范围为0.01到0.1，单位为像素每秒(pixel/s)。
6.4令i＝i+1；
6.5若i≤n，转6.2；否则说明得到了tk时刻n个特征点的光流，即得到了tk时刻光流集合Gk，
Figure BDA0002142074750000088
将FDk、PBDk和Gk发送给特征点更新模块，转第七步。
第七步，特征点更新模块从匹配模块接收tk时刻n个特征点的位置集合FDk、tk时刻n个特征点对应的模板边的位置集合PBDk和n个特征点的光流集合Gk，通过光流计算tk+1时刻n个特征点的位置集合FDk+1，将Gk和FDk+1发送给事件集选取模块，将FDk+1和PBDk发送给模板边更新模块。方法是：
7.1令i＝1；
7.2计算tk+1时刻fi的位置fi(tk+1)，将fi(tk+1)放到集合FDk+1中，
Figure BDA0002142074750000089
Figure BDA00021420747500000810
为光流与时间相乘，表示从tk时刻到tk+1时刻特征点fi运动距离。
7.3令i＝i+1；
7.4若i≤n，转7.2；否则说明得到了tk+1时刻n个特征点的位置，并得到集合FDk+1，FDk+1＝{f1(tk+1),...,fi(tk+1),...,fn(tk+1)}。将Gk和FDk+1发送给事件集选取模块，将FDk+1和PBDk发送给模板边更新模块；并将tk+1时刻n个特征点的位置集合FDk+1进行显示或存储到结果文件中供用户查看，转第八步。
第八步，模板边更新模块从数据获取模块接收IMU数据，从特征点更新模块接收tk+1时刻n个特征点的位置集合FDk+1和tk时刻特征点周围模板边的位置集合PBDk，使用IMU数据对PBDk中的位置进行更新，得到tk+1时刻n个特征点对应的模板边的位置集合PBDk+1，将PBDk+1发送给匹配模块，。具体方法如下：
8.1令i＝1；
8.2更新特征点fi对应的模板边在tk+1时刻的位置，方法是：
8.2.1令j＝1；
8.2.2定义符号F为与fi相对应的三维空间中的点，Pj为与模板点pj相对应的三维空间中的点，F和Pj均用三维坐标(x,y,z)的形式表示。将模板点pj使用齐次坐标形式表示，格式为(x,y,1)，前两维分别表示pj在像素坐标系下的横坐标和纵坐标，于是如下方程成立：
在tk时刻，
Figure BDA0002142074750000091
在tk+1时刻，
Figure BDA0002142074750000092
其中K为事件相机的内参矩阵，为事件相机自带参数。R为从tk时刻到tk+1时刻事件相机的旋转矩阵，t为从tk时刻到tk+1时刻事件相机的平移向量，均通过获取到的IMU数据计算得到。
Figure BDA0002142074750000093
为Pj在tk时刻相机坐标系下的深度、
Figure BDA0002142074750000094
为Pj在tk+1时刻相机坐标系下的深度，
Figure BDA0002142074750000095
为F在tk时刻相机坐标系下的深度、
Figure BDA0002142074750000096
为F在tk+1时刻相机坐标系下的深度。
8.2.3将公式(10)中的两个等式相减，得到tk+1时刻模板点pj相对于fi的相对位置
Figure BDA0002142074750000097
计算公式如下：
Figure BDA0002142074750000098
由于在像素坐标系下模板点pj在特征点fi的临近区域内，所以在空间中对应的点Pj和F同样相近，在该领域中，此类情况即可认为Pj和F在tk+1时刻的相机坐标系中具有相同的深度，即
Figure BDA0002142074750000101
将公式(9)带入公式(11)，并对公式(11)进行化简得到：
Figure BDA0002142074750000102
考虑pj和fi均为齐次坐标表示，进一步对公式(12)化简得到tk+1时刻模板点相对位置的公式为：
Figure BDA0002142074750000103
符号Nor()表示齐次化操作，即将括号内的坐标转化为齐次坐标。通过公式(13)得到更新后的tk+1时刻模板点pj相对于fi的相对位置
Figure BDA0002142074750000104
按公式(14)得到tk+1时刻模板点pj的位置pj(tk+1)，将pj(tk+1)放到tk+1时刻特征点fi周围模板边的位置集合
Figure BDA00021420747500001010
Figure BDA0002142074750000105
8.2.4令j＝j+1；
8.2.5如果j≤mi，转8.2.2；否则表示得到了
Figure BDA0002142074750000106
Figure BDA0002142074750000107
将
Figure BDA0002142074750000108
放到tk+1时刻n个特征点周围模板边的位置集合PBDk+1中，转8.3；
8.3令i＝i+1；
8.4若i≤n，转8.2；否则说明得到PBDk+1，
Figure BDA0002142074750000109
将PBDk+1发送给匹配模块，转第九步。
第九步，令k＝k+1。
第十步，若k<N，转第四步；否则结束。
采用本发明可以达到以下技术效果：
1.本发明同时利用了DAVIS输出的图像帧和事件流两种数据信息，首先使用Harris角点检测法提取特征点并使用Canny边检测法提取边图，然后从事件流中选择空间时间窗口,并使用EM-ICP算法将窗口内的事件与边图上的模板点进行匹配，估计特征点的光流，进而通过光流跟踪特征点，这样的跟踪充分利用了图像帧中的边缘信息和又利用了事件流的异步性，在保证异步跟踪的同时，提高了特征点跟踪的精度。
2.在模板更新时，本发明使用IMU更新模板点的位置，在跟踪过程中提高了模板边的跟踪精度，从而提高了计算光流的准确性，进一步提高了特征点的跟踪精度。
本发明在苏黎世大学发布的事件相机数据集“The Event-Camera Dataset andSimulator”(事件相机数据集和模拟器)上进行了实验验证，并与Zhu等人提出的特征点跟踪方法、Kueng等人提出的方法做了对比实验，实验结果表明本发明不仅提高了特征点跟踪精度，且延长了特征点跟踪时间。
附图说明
图1为本发明总体流程图；
图2为本发明第一步构建的基于事件相机的特征点跟踪***逻辑结构图；
图3为本发明与现有跟踪方法平均跟踪精度误差对比实验结果；
图4为本发明与现有跟踪方法平均跟踪时间对比实验结果。
具体实施方式
图1为本发明总体流程图；如图1所示，本发明包括以下步骤：
第一步，构建基于事件相机的特征点跟踪***。事件相机的特征点跟踪***如图2所示，由数据获取模块、初始化模块、事件集选取模块、匹配模块、特征点更新模块、模板边更新模块组成。
数据获取模块与初始化模块、事件集选取模块、模板边更新模块相连。数据集获取模块从苏黎世大学的公开事件相机数据集“The Event-Camera Dataset and Simulator”上下载数据，并将获取到的图像帧发送给初始化模块，将事件流发送给事件集选取模块，将IMU数据发送给模板边更新模块。
初始化模块与数据获取模块、事件集选取模块、匹配模块相连。初始化模块从数据获取模块接收图像帧，从图像帧中提取特征点和边图，得到特征点的位置和特征点周围的模板边的位置，将特征点的位置发送给事件集选取模块，将特征点周围的模板边的位置发送给匹配模块。
事件集选取模块与数据获取模块、初始化模块、特征点更新模块、匹配模块相连，事件集选取模块从数据获取模块接收事件流，从初始化模块(第一次循环)或特征点更新模块(从第二次循环开始)接收特征点的位置，从特征点更新模块接收特征点的光流(从第二次循环开始)，在每个特征点周围从事件流中选取特征点的事件集，将每个特征点的位置以及对应的事件集发送给匹配模块。
匹配模块与初始化模块、事件集选取模块、特征点更新模块、模板边更新模块相连，匹配模块从事件集选取模块接收特征点的位置以及每个特征点对应的事件集，从初始化模块(第一次循环)或模板边更新模块(从第二次循环开始)接收特征点周围的模板边的位置，将特征点的事件集与特征点周围的模板边进行匹配，求出每个特征点的光流，将每个特征点的位置、特征点周围的模板边的位置和特征点的光流发送给特征点更新模块。
特征点更新模块与匹配模块、模板边更新模块、事件集选取模块相连，特征点更新模块从匹配模块接收特征点的位置、特征点周围的模板边的位置、特征点的光流，利用特征点的光流计算特征点新的位置，将特征点周围的模板边的位置以及特征点新的位置发送给模板边更新模块，将特征点的光流和特征点新的位置发送给事件集选取模块，并将特征点新的位置(即特征点跟踪结果)输出，供用户查看。
模板边更新模块与数据获取模块、特征点更新模块、匹配模块相连，模板边更新模块从特征点更新模块接收特征点周围的模板边的位置和特征点新的位置，从数据获取模块接收IMU数据，通过IMU数据来更新特征点周围的模板边的位置，将模板边更新后的位置发送给匹配模块。
第二步，数据获取模块从事件相机数据集“The Event-Camera Dataset andSimulator”中获取图像帧、事件流和IMU数据。
第三步，采用基于事件相机的特征点跟踪***对从t0时刻开始到tN时刻结束的时间段内数据获取模块得到的事件流中的特征点进行跟踪。跟踪过程中将时间区间[t0,tN]分成一系列的子时间区间[t0,t1],...,[tk,tk+1],...,[tN-1,tN]，第k+1个时间子区间表示为[tk,tk+1]，N表示子时间区间个数，N的大小由事件流时间长度决定，取值范围为N≥1。跟踪过程如下：
3.1，令时间序号k＝0；
3.2，初始化模块进行初始化操作，具体步骤如下：
3.2.1初始化模块采用Harris角点检测法从数据获取模块得到的图像帧中提取特征点，将提取到的特征点放到集合FD中，令FD为{f1,...,fi,...,fn}，fi代表检测到的第i个特征点，n为特征点个数。将特征点的位置看做时间的函数，将t0时刻FD中n个特征点的位置放在特征点位置集合FD0中，FD0＝{f1(t0),...,fi(t0),...,fn(t0)}，fi(t0)代表特征点fi在t0时刻的位置，将FD0发送给事件集选取模块。
3.2.2初始化模块使用Canny边检测法从数据获取模块得到的图像帧中提取边图，每个图像帧对应一个边图。
3.2.3初始化模块选取FD中n个特征点对应的模板边(由模板边上所有模板点的集合表示)，将FD中n个特征点对应的t0时刻模板边的位置集合PBD0发送给匹配模块。方法如下：
3.2.3.1令i＝1；
3.2.3.2以特征点fi在t0时刻的位置fi(t0)为中心，在fi(t0)周围选取一个矩形区域
Figure BDA0002142074750000131
大小为s×s，即
Figure BDA0002142074750000132
的长和宽均为s，s的范围为20到30个像素。将3.2.2检测到的边图中在
Figure BDA0002142074750000133
内的边作为fi的模板边。模板边上的像素点即为fi的模板点。定义
Figure BDA0002142074750000134
内的模板点的集合
Figure BDA0002142074750000135
Figure BDA0002142074750000136
也即fi对应的模板边，将
Figure BDA0002142074750000137
放到模板边集合PB中：
Figure BDA0002142074750000138
pj为fi的第j个模板点。将pj的位置视作时间的函数，pj(t0)表示pj在t0时刻的位置。
Figure BDA0002142074750000139
表示在t0时刻pj在矩形区域
Figure BDA00021420747500001310
内。mi表示集合
Figure BDA00021420747500001311
中模板点的个数。
3.2.3.3令i＝i+1；
3.2.3.4若i≤n，转3.2.3.2；否则，表示得到了FD中n个特征点对应的模板边构成的模板边集合PB，
Figure BDA00021420747500001312
定义
Figure BDA00021420747500001313
中模板点在t0时刻的位置的集合
Figure BDA00021420747500001314
令PB中n个模板边中模板点在t0时刻的位置集合
Figure BDA00021420747500001315
将PBD0发送给匹配模块，转第四步。
第四步，事件集选取模块从数据获取模块接收事件流，根据不同的k值分别从初始化模块或特征点更新模块接收不同的数据，在n个特征点周围从事件流中选取特征点的事件集S，将FDk和特征点的事件集S发送给匹配模块，方法是：
4.1若k＝0，事件集选取模块从数据获取模块接收事件流，从初始化模块接收tk时刻n个特征点的位置集合FD0，令t1＝t0+1，单位为秒，转4.3；
4.2若k≥1，事件集选取模块从数据获取模块接收事件流，从特征点更新模块接收tk时刻n个特征点的位置集合FDk和tk-1时刻n个特征点的光流集合Gk-1，
Figure BDA00021420747500001316
其中
Figure BDA0002142074750000141
表示tk-1时刻特征点fi的光流，预估子时间区间[tk,tk+1]的大小，根据公式(2)计算tk+1：
Figure BDA0002142074750000142
其中
Figure BDA0002142074750000143
表示在子时间区间[tk-1,tk]上特征点fi的光流，从特征点更新模块得到。
Figure BDA0002142074750000144
表示在子时间区间[tk-1,tk]上所有特征点的平均光流。通过公式(2)计算tk+1的物理意义在于，上一时间区间内的特征点平均运动3个像素需要的时间作为当前区间[tk,tk+1]大小的预估值。转4.3；
4.3在tk时刻n个特征点的位置集合FDk中每个位置的周围选取每个特征点对应的事件集。
4.3.1令i＝1；
4.3.2从事件流中选出符合公式(3)要求的事件，将这些事件放到特征点fi在tk时刻对应的事件集
Figure BDA0002142074750000145
中，并将
Figure BDA0002142074750000146
放到事件集S中：
Figure BDA0002142074750000147
Figure BDA0002142074750000148
表示在一个三维时空区域内的事件的集合，空间范围为特征点fi在tk时刻的位置fi(tk)周围的矩形区域
Figure BDA0002142074750000149
时间范围为区间[tk,tk+1]。ed代表
Figure BDA00021420747500001410
中的第d个事件，d＝{1,2,...,zi}，zi表示
Figure BDA00021420747500001411
内事件的个数。事件ed的表示形式为
Figure BDA00021420747500001412
xd表示事件ed在像素坐标系中的坐标，
Figure BDA00021420747500001413
表示事件ed发生的时间，
Figure BDA00021420747500001414
表示ed的像素坐标在
Figure BDA00021420747500001415
内。
4.3.3令i＝i+1；
4.3.4若i≤n，转4.3.2；否则，说明得到了tk时刻FD中n个特征点对应的事件集S,
Figure BDA00021420747500001416
并将FDk和S发送给匹配模块，转第五步。
第五步，匹配模块从事件集选取模块接收FDk和S，若k＝0，从初始化模块接收n个特征点对应的模板边的位置集合PBDk，转第六步；若k≥1，从模板边更新模块接收n个特征点对应的模板边的位置集合PBDk，转第六步；
第六步，匹配模块将S与特征点周围的模板边进行匹配，求出tk时刻n个特征点的光流集合Gk，将n个特征点的位置集合FDk、n个特征点对应的模板边的位置集合PBDk和光流集合Gk发送给特征点更新模块。具体方法为：
6.1令i＝1；
6.2对特征点fi构造匹配误差函数：
对于
Figure BDA0002142074750000151
内的事件ed，修正事件ed在tk时刻的位置，公式如下：
Figure BDA0002142074750000152
x′d表示计算得到的事件ed在tk时刻的位置，简称为事件ed运动修正的位置，
Figure BDA0002142074750000153
表示在时间区间[tk,tk+1]上特征点fi的光流。符号·表示点乘，下文中出现该符号的地方均表示点乘，在不引起歧义的情况下，下文中部分公式对该符号做了省略。为了表示方便，定义符号
Figure BDA0002142074750000154
构造匹配误差函数如下：
Figure BDA0002142074750000155
ε表示误差，pj(tk)代表tk时刻模板点pj的位置。rdj表示ed由模板点pj生成的概率，即ed与pj匹配的概率。这里将rdj作为权重，rdj为事件ed与模板点pj相匹配的概率，rdj的计算公式如下：
Figure BDA0002142074750000156
6.3采用EM-ICP算法最小化6.2步构造得到的匹配误差函数，求解得到最优的光流
Figure BDA0002142074750000157
求解过程如下：
6.3.1初始化
Figure BDA0002142074750000158
6.3.2通过公式(6)计算rdj；
6.3.3更新光流：
Figure BDA0002142074750000161
6.3.4计算光流的变化量
Figure BDA0002142074750000162
并令
Figure BDA0002142074750000163
将
Figure BDA0002142074750000164
放到tk时刻光流集合Gk中。
6.3.5若
Figure BDA0002142074750000165
表示得到最终的优化结果
Figure BDA0002142074750000166
转6.4；若
Figure BDA0002142074750000167
转6.3.2。σ为光流的变化量阈值，取值范围为0.01到0.1，单位为像素每秒。
6.4令i＝i+1；
6.5若i≤n，转6.2；否则说明得到了tk时刻n个特征点的光流，即得到了tk时刻光流集合Gk，
Figure BDA0002142074750000168
将FDk、PBDk和Gk发送给特征点更新模块，转第七步。
第七步，特征点更新模块从匹配模块接收tk时刻n个特征点的位置集合FDk、tk时刻n个特征点对应的模板边的位置集合PBDk和n个特征点的光流集合Gk，通过光流计算tk+1时刻n个特征点的位置集合FDk+1，将Gk和FDk+1发送给事件集选取模块，将FDk+1和PBDk发送给模板边更新模块。方法是：
7.1令i＝1；
7.2计算tk+1时刻fi的位置fi(tk+1)，将fi(tk+1)放到集合FDk+1中，
Figure BDA0002142074750000169
Figure BDA00021420747500001610
为光流与时间相乘，表示从tk时刻到tk+1时刻特征点fi运动距离。
7.3令i＝i+1；
7.4若i≤n，转7.2；否则说明得到了tk+1时刻n个特征点的位置，并得到集合FDk+1，FDk+1＝{f1(tk+1),...,fi(tk+1),...,fn(tk+1)}。将Gk和FDk+1发送给事件集选取模块，将FDk+1和PBDk发送给模板边更新模块；并将tk+1时刻n个特征点的位置集合FDk+1进行显示或存储到结果文件中供用户查看，转第八步。
第八步，模板边更新模块从数据获取模块接收IMU数据，从特征点更新模块接收tk+1时刻n个特征点的位置集合FDk+1和tk时刻特征点周围模板边的位置集合PBDk，使用IMU数据对PBDk中的位置进行更新，得到tk+1时刻n个特征点对应的模板边的位置集合PBDk+1，将PBDk+1发送给匹配模块，。具体方法如下：
8.1令i＝1；
8.2更新特征点fi对应的模板边在tk+1时刻的位置，方法是：
8.2.1令j＝1；
8.2.2定义符号F为与fi相对应的三维空间中的点，Pj为与模板点pj相对应的三维空间中的点，F和Pj均用三维坐标(x,y,z)的形式表示。将模板点pj使用齐次坐标形式表示，格式为(x,y,1)，前两维分别表示pj在像素坐标系下的横坐标和纵坐标，于是如下方程成立：
在tk时刻，
Figure BDA0002142074750000171
在tk+1时刻，
Figure BDA0002142074750000172
其中K为事件相机的内参矩阵，为事件相机自带参数。R为从tk时刻到tk+1时刻事件相机的旋转矩阵，t为从tk时刻到tk+1时刻事件相机的平移向量，均通过获取到的IMU数据计算得到。
Figure BDA0002142074750000173
为Pj在tk时刻相机坐标系下的深度、
Figure BDA0002142074750000174
为Pj在tk+1时刻相机坐标系下的深度，
Figure BDA0002142074750000175
为F在tk时刻相机坐标系下的深度、
Figure BDA0002142074750000176
为F在tk+1时刻相机坐标系下的深度。
8.2.3将公式(10)中的两个等式相减，得到tk+1时刻模板点pj相对于fi的相对位置
Figure BDA0002142074750000177
计算公式如下：
Figure BDA0002142074750000178
由于在像素坐标系下模板点pj在特征点fi的临近区域内，所以在空间中对应的点Pj和F同样相近，在该领域中，此类情况即可认为Pj和F在tk+1时刻的相机坐标系中具有相同的深度，即
Figure BDA0002142074750000181
将公式(9)带入公式(11)，并对公式(11)进行化简得到：
Figure BDA0002142074750000182
考虑pj和fi均为齐次坐标表示，进一步对公式(12)化简得到tk+1时刻模板点相对位置的公式为：
Figure BDA0002142074750000183
符号Nor()表示齐次化操作，即将括号内的坐标转化为齐次坐标。通过公式(13)得到更新后的tk+1时刻模板点pj相对于fi的相对位置
Figure BDA0002142074750000184
按公式(14)得到tk+1时刻模板点pj的位置pj(tk+1)，将pj(tk+1)放到tk+1时刻特征点fi周围模板边的位置集合
Figure BDA0002142074750000185
Figure BDA0002142074750000186
8.2.4令j＝j+1；
8.2.5如果j≤mi，转8.2.2；否则表示得到了
Figure BDA0002142074750000187
Figure BDA0002142074750000188
将
Figure BDA0002142074750000189
放到tk+1时刻n个特征点周围模板边的位置集合PBDk+1中，转8.3；
8.3令i＝i+1；
8.4若i≤n，转8.2；否则说明得到PBDk+1，
Figure BDA00021420747500001810
将PBDk+1发送给匹配模块，转第九步。
第九步，令k＝k+1。
第十步，若k<N，转第四步；否则结束。
图3为本发明与现有跟踪方法平均跟踪精度误差对比实验结果；该实验结果为在“The Event-Camera Dataset and Simulator”数据集上采用本发明进行测试得到的结果。实验环境为一台配置为i7 2.8GHz CPU、8G RAM的笔记本。该实验的评价指标为特征点的平均跟踪误差，单位为像素。图中左侧为该数据集中数据序列的名称，右侧上方为特征点的平均跟踪误差。图中三列实验数据分别为本发明、Zhu等人的方法、Kueng等人的方法在相同的测试数据序列、同样的实验环境下测试的结果。实验结果显示，与其它两个方法相比，本发明在所有测试数据序列上都具有更低的平均跟踪误差。图中“×”表示无数据。
图4为本发明与现有跟踪方法平均跟踪时间对比实验结果。该实验与图3对应的实验的测试数据集和实验环境相同。该实验的评价指标为特征点的平均跟踪时间，单位为秒。图中左侧为该数据集中数据序列的名称，右侧上方为特征点的平均跟踪时间。图中三列实验数据分别为本***、Zhu等人的***、Kueng等人的***在相同的测试数据系列、同样的实验环境下测试的结果。实验结果显示，与其它两个方法相比，本发明在除数据序列“boxes_translation”和“boxes_6dof”之外的其他数据序列上，均能够实现更长时间的跟踪。

Claims (5)

1.一种基于事件相机的特征点跟踪方法，其特征在于包括以下步骤：
第一步，构建基于事件相机的特征点跟踪***，事件相机的特征点跟踪***由数据获取模块、初始化模块、事件集选取模块、匹配模块、特征点更新模块、模板边更新模块组成；
数据获取模块与初始化模块、事件集选取模块、模板边更新模块相连；数据集获取模块从事件相机数据集上下载数据，将获取到的图像帧发送给初始化模块，将事件流发送给事件集选取模块，将IMU数据即惯性测量单元数据发送给模板边更新模块；
初始化模块与数据获取模块、事件集选取模块、匹配模块相连；初始化模块从数据获取模块接收图像帧，从图像帧中提取特征点和边图，得到特征点的位置和特征点周围的模板边的位置，将特征点的位置发送给事件集选取模块，将特征点周围的模板边的位置发送给匹配模块；
事件集选取模块与数据获取模块、初始化模块、特征点更新模块、匹配模块相连，事件集选取模块从数据获取模块接收事件流，从初始化模块或特征点更新模块接收特征点的位置，从特征点更新模块接收特征点的光流，在每个特征点周围从事件流中选取特征点的事件集，将每个特征点的位置以及对应的事件集发送给匹配模块；
匹配模块与初始化模块、事件集选取模块、特征点更新模块、模板边更新模块相连，匹配模块从事件集选取模块接收特征点的位置以及每个特征点对应的事件集，从初始化模块或模板边更新模块接收特征点周围的模板边的位置，将特征点的事件集与特征点周围的模板边进行匹配，求出每个特征点的光流，将每个特征点的位置、特征点周围的模板边的位置和特征点的光流发送给特征点更新模块；
特征点更新模块与匹配模块、模板边更新模块、事件集选取模块相连，特征点更新模块从匹配模块接收特征点的位置、特征点周围的模板边的位置、特征点的光流，利用特征点的光流计算特征点新的位置，将特征点周围的模板边的位置以及特征点新的位置发送给模板边更新模块，将特征点的光流和特征点新的位置发送给事件集选取模块，并将特征点新的位置输出；
模板边更新模块与数据获取模块、特征点更新模块、匹配模块相连，模板边更新模块从特征点更新模块接收特征点周围的模板边的位置和特征点新的位置，从数据获取模块接收IMU数据，通过IMU数据来更新特征点周围的模板边的位置，将模板边更新后的位置发送给匹配模块；
第二步，数据获取模块从事件相机数据集中获取图像帧、事件流和IMU数据；
第三步，采用基于事件相机的特征点跟踪***对从t0时刻开始到tN时刻结束的时间段内数据获取模块得到的事件流中的特征点进行跟踪，将时间区间[t0,tN]分成一系列的子时间区间[t0,t1],...,[tk,tk+1],...,[tN-1,tN]，N表示子时间区间个数，N的大小由事件流时间长度决定，取值范围为N≥1；第k+1个时间子区间表示为[tk,tk+1]；跟踪过程如下：
3.1，令时间序号k＝0；
3.2，若k＝0，初始化模块进行初始化操作，方法为：
3.2.1初始化模块采用Harris角点检测法从数据获取模块得到的图像帧中提取特征点，将提取到的特征点放到集合FD中，令FD为{f0,...,fi,...,fn}，fi代表检测到的第i个特征点，n为特征点个数；将特征点的位置看做时间的函数，将t0时刻FD中n个特征点的位置放在特征点位置集合FD0中，FD0＝{f0(t0),...,fi(t0),...,fn(t0)}，fi(t0)代表特征点fi在t0时刻的位置，将FD0发送给事件集选取模块；
3.2.2初始化模块使用Canny边检测法从数据获取模块得到的图像帧中提取边图，每个图像帧对应一个边图；
3.2.3初始化模块选取FD中n个特征点对应的模板边，将FD中n个特征点对应的t0时刻模板边的位置集合PBD0发送给匹配模块，方法为：
3.2.3.1令i＝1；
3.2.3.2以特征点fi在t0时刻的位置fi(t0)为中心，在fi(t0)周围选取一个矩形区域
Figure FDA0002904542740000021
大小为s×s，即该矩形的长和宽均为s；将3.2.2检测到的边图中在
Figure FDA0002904542740000022
内的边作为fi的模板边，模板边上的像素点即为fi的模板点；定义
Figure FDA0002904542740000023
内的模板点的集合
Figure FDA0002904542740000024
Figure FDA0002904542740000025
也即fi对应的模板边，将
Figure FDA0002904542740000026
放到模板边集合PB中,
Figure FDA0002904542740000027
pj为fi的第j个模板点，将pj的位置视作时间的函数，pj(t0)表示pj在t0时刻的位置，
Figure FDA0002904542740000028
表示在t0时刻pj在矩形区域
Figure FDA0002904542740000029
内，mi表示
Figure FDA00029045427400000210
中模板点的个数；
3.2.3.3令i＝i+1；
3.2.3.4若i≤n，转3.2.3.2；否则，表示得到了FD中n个特征点对应的模板边构成的模板边集合PB，
Figure FDA0002904542740000031
定义
Figure FDA0002904542740000032
中模板点在t0时刻的位置的集合
Figure FDA0002904542740000033
令PB中n个模板边中模板点在t0时刻的位置集合
Figure FDA0002904542740000034
将PBD0发送给匹配模块，转第四步；
第四步，事件集选取模块从数据获取模块接收事件流，根据不同的k值分别从初始化模块或特征点更新模块接收不同的数据，在n个特征点周围从事件流中选取特征点的事件集S，将FDk和特征点的事件集S发送给匹配模块，方法是：
4.1若k＝0，事件集选取模块从数据获取模块接收事件流，从初始化模块接收tk时刻n个特征点的位置集合FD0，令t1＝t0+1，单位为秒，转4.3；
4.2若k≥1，事件集选取模块从数据获取模块接收事件流，从特征点更新模块接收tk时刻n个特征点的位置集合FDk和tk-1时刻n个特征点的光流集合Gk-1，
Figure FDA0002904542740000035
预估子时间区间[tk,tk+1]的大小，根据公式(2)计算tk+1：
Figure FDA0002904542740000036
其中
Figure FDA0002904542740000037
表示在子时间区间[tk-1,tk]上特征点fi的光流；
Figure FDA0002904542740000038
表示在子时间区间[tk-1,tk]上所有特征点的平均光流，转4.3；
4.3在tk时刻n个特征点的位置集合FDk中每个位置的周围选取每个特征点对应的事件集，方法是：
4.3.1令i＝1；
4.3.2从事件流中选出符合公式(3)要求的事件，将这些事件放到特征点fi在tk时刻对应的事件集
Figure FDA0002904542740000039
中，并将
Figure FDA00029045427400000310
放到事件集S中：
Figure FDA00029045427400000311
Figure FDA00029045427400000312
表示在一个三维时空区域内的事件的集合，空间范围为特征点fi在tk时刻的位置fi(tk)周围的矩形区域
Figure FDA00029045427400000313
时间范围为区间[tk,tk+1]；ed代表
Figure FDA00029045427400000314
中的第d个事件，d＝{1,2,...,zi}，zi表示
Figure FDA0002904542740000041
内事件的个数；
Figure FDA0002904542740000042
xd表示事件ed在像素坐标系中的坐标，
Figure FDA0002904542740000043
表示事件ed发生的时间，
Figure FDA0002904542740000044
表示ed的像素坐标在
Figure FDA0002904542740000045
内；
4.3.3令i＝i+1；
4.3.4若i≤n，转4.3.2；否则，说明得到了tk时刻FD中n个特征点对应的事件集S，
Figure FDA0002904542740000046
将FDk和S发送给匹配模块，转第五步；
第五步，匹配模块从事件集选取模块接收FDk和S，若k＝0，从初始化模块接收n个特征点对应的模板边的位置集合PBDk，转第六步；若k≥1，从模板边更新模块接收n个特征点对应的模板边的位置集合PBDk，转第六步；
第六步，匹配模块将S与特征点周围的模板边进行匹配，求出tk时刻n个特征点的光流集合Gk，将n个特征点的位置集合FDk、n个特征点对应的模板边的位置集合PBDk和光流集合Gk发送给特征点更新模块，具体方法为：
6.1令i＝1；
6.2对特征点fi构造匹配误差函数，方法是：
按公式(4)修正事件ed在tk时刻的位置：
Figure FDA0002904542740000047
x′d表示计算得到的事件ed在tk时刻的位置，简称为事件ed运动修正的位置，
Figure FDA0002904542740000048
表示在时间区间[tk,tk+1]上特征点fi的光流；符号·表示点乘，定义符号
Figure FDA0002904542740000049
构造匹配误差函数如下：
Figure FDA00029045427400000410
e表示误差，pj(tk)代表tk时刻模板点pj的位置；rdj表示ed由模板点pj生成的概率，即ed与pj匹配的概率；公式中双竖线表示对双竖线内的向量进行取模操作，rdj的计算公式如下：
Figure FDA00029045427400000411
6.3采用EM-ICP算法最小化匹配误差函数，求解得到最优的光流
Figure FDA0002904542740000051
方法是：
6.3.1初始化
Figure FDA0002904542740000052
6.3.2通过公式(6)计算rdj；
6.3.3更新光流：
Figure FDA0002904542740000053
6.3.4计算光流的变化量
Figure FDA0002904542740000054
并令
Figure FDA0002904542740000055
将
Figure FDA0002904542740000056
放到tk时刻光流集合Gk中；
6.3.5若
Figure FDA0002904542740000057
表示得到最终的优化结果
Figure FDA0002904542740000058
转6.4；若
Figure FDA0002904542740000059
转6.3.2；σ为光流的变化量阈值；
6.4令i＝i+1；
6.5若i≤n，转6.2；否则说明得到了tk时刻光流集合Gk，
Figure FDA00029045427400000510
将FDk、PBDk和Gk发送给特征点更新模块，转第七步；
第七步，特征点更新模块从匹配模块接收FDk、PBDk和Gk，通过光流计算tk+1时刻n个特征点的位置集合FDk+1，将Gk和FDk+1发送给事件集选取模块，将FDk+1和PBDk发送给模板边更新模块，方法是：
7.1令i＝1；
7.2计算tk+1时刻fi的位置fi(tk+1)，将fi(tk+1)放到集合FDk+1中，
Figure FDA00029045427400000511
Figure FDA00029045427400000512
为光流与时间相乘，表示从tk时刻到tk+1时刻特征点fi运动距离；
7.3令i＝i+1；
7.4若i≤n，转7.2；否则说明得到了tk+1时刻n个特征点的位置，并得到集合FDk+1，FDk+1＝{f1(tk+1),...,fi(tk+1),...,fn(tk+1)}；将Gk和FDk+1发送给事件集选取模块，将FDk+1和PBDk发送给模板边更新模块；并将tk+1时刻n个特征点的位置集合FDk+1进行显示或存储到结果文件中，转第八步；
第八步，模板边更新模块从数据获取模块接收IMU数据，从特征点更新模块接收FDk+1和PBDk，使用IMU数据对PBDk中的位置进行更新，得到tk+1时刻n个特征点对应的模板边的位置集合PBDk+1，将PBDk+1发送给匹配模块，方法如下：
8.1令i＝1；
8.2更新特征点fi对应的模板边在tk+1时刻的位置，方法是：
8.2.1令j＝1；
8.2.2定义符号F为与fi相对应的三维空间中的点，Pj为与模板点pj相对应的三维空间中的点，F和Pj均用三维坐标(x,y,z)的形式表示；将模板点pj使用齐次坐标形式表示，格式为(x,y,1)，前两维分别表示pj在像素坐标系下的横坐标和纵坐标，得到如下方程：
在tk时刻，
Figure FDA0002904542740000061
在tk+1时刻，
Figure FDA0002904542740000062
其中K为事件相机的内参矩阵，为事件相机自带参数；R为从tk时刻到tk+1时刻事件相机的旋转矩阵，t为从tk时刻到tk+1时刻事件相机的平移向量，均通过获取到的IMU数据计算得到；
Figure FDA0002904542740000063
为Pj在tk时刻相机坐标系下的深度、
Figure FDA0002904542740000064
为Pj在tk+1时刻相机坐标系下的深度，
Figure FDA0002904542740000065
为F在tk时刻相机坐标系下的深度、
Figure FDA0002904542740000066
为F在tk+1时刻相机坐标系下的深度；
8.2.3将公式(10)中的两个等式相减，得到tk+1时刻模板点pj相对于fi的相对位置
Figure FDA0002904542740000067
Figure FDA0002904542740000068
Figure FDA0002904542740000071
Pj和F在tk+1时刻的相机坐标系中具有相同的深度，即
Figure FDA0002904542740000072
将公式(9)带入公式(11)，并对公式(11)进行化简得到：
Figure FDA0002904542740000073
对公式(12)化简得到tk+1时刻模板点相对位置的公式为：
Figure FDA0002904542740000074
符号Nor()表示齐次化操作，即将括号内的坐标转化为齐次坐标；按公式(14)得到tk+1时刻模板点pj的位置pj(tk+1)，将pj(tk+1)放到tk+1时刻特征点fi周围模板边的位置集合
Figure FDA0002904542740000075
Figure FDA0002904542740000076
8.2.4令j＝j+1；
8.2.5如果j≤mi，转8.2.2；否则表示得到了
Figure FDA0002904542740000077
Figure FDA0002904542740000078
将
Figure FDA0002904542740000079
放到tk+1时刻n个特征点周围模板边的位置集合PBDk+1中，转8.3；
8.3令i＝i+1；
8.4若i≤n，转8.2；否则说明得到PBDk+1，
Figure FDA00029045427400000710
将PBDk+1发送给匹配模块，转第九步；
第九步，令k＝k+1；
第十步，若k<N，转第四步；否则结束。
2.如权利要求1所述的一种基于事件相机的特征点跟踪方法，其特征在于所述事件相机数据集指由DAVIS采集的“事件相机数据集和模拟器”即“The Event-Camera Datasetand Simulator”，事件相机数据集中包括图像帧、事件流和IMU数据；DAVIS指动态和主动式像素传感器。
3.如权利要求1所述的一种基于事件相机的特征点跟踪方法，其特征在于所述模板边由模板边上所有模板点的集合表示。
4.如权利要求1所述的一种基于事件相机的特征点跟踪方法，其特征在于3.2.3.2步所述s的范围为20到30个像素。
5.如权利要求1所述的一种基于事件相机的特征点跟踪方法，其特征在于6.3.5步所述σ的取值范围为0.01到0.1，单位为像素每秒。
CN201910672162.2A 2019-07-24 2019-07-24 一种基于事件相机的特征点跟踪方法 Active CN110390685B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910672162.2A CN110390685B (zh) 2019-07-24 2019-07-24 一种基于事件相机的特征点跟踪方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910672162.2A CN110390685B (zh) 2019-07-24 2019-07-24 一种基于事件相机的特征点跟踪方法

Publications (2)

Publication Number Publication Date
CN110390685A CN110390685A (zh) 2019-10-29
CN110390685B true CN110390685B (zh) 2021-03-09

Family

ID=68287327

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910672162.2A Active CN110390685B (zh) 2019-07-24 2019-07-24 一种基于事件相机的特征点跟踪方法

Country Status (1)

Country Link
CN (1) CN110390685B (zh)

Families Citing this family (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111062337B (zh) * 2019-12-19 2023-08-04 北京迈格威科技有限公司 人流方向检测方法及装置、存储介质和电子设备
CN113066104B (zh) * 2021-03-25 2024-04-19 三星（中国）半导体有限公司 角点检测方法和角点检测装置
CN113160218B (zh) * 2021-05-12 2023-06-20 深圳龙岗智能视听研究院 一种基于事件相机的检测物体运动强度的方法
CN113506321A (zh) * 2021-07-15 2021-10-15 清华大学 图像处理方法及装置、电子设备和存储介质
CN115984336A (zh) * 2021-10-14 2023-04-18 华为技术有限公司 一种光流估计方法和装置
CN116188533B (zh) * 2023-04-23 2023-08-08 深圳时识科技有限公司 特征点跟踪方法与装置、电子设备

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109697726A (zh) * 2019-01-09 2019-04-30 厦门大学 一种基于事件相机的端对端目标运动估计方法
WO2019099337A1 (en) * 2017-11-14 2019-05-23 Kaban Technologies Llc Event camera-based deformable object tracking

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109426782B (zh) * 2017-08-29 2023-09-19 北京三星通信技术研究有限公司 对象检测方法和用于对象检测的神经网络***
CN107808407B (zh) * 2017-10-16 2020-12-18 亿航智能设备（广州）有限公司 基于双目相机的无人机视觉slam方法、无人机及存储介质
CN109934862A (zh) * 2019-02-22 2019-06-25 上海大学 一种点线特征结合的双目视觉slam方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2019099337A1 (en) * 2017-11-14 2019-05-23 Kaban Technologies Llc Event camera-based deformable object tracking
CN109697726A (zh) * 2019-01-09 2019-04-30 厦门大学 一种基于事件相机的端对端目标运动估计方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
Asynchronous Corner Detection and Tracking for Event Cameras in Real Time;Ignacio Alzugaray et al.;《IEEE Robotics and Automation Letters 》;20180622;第3卷(第4期);第3177-3184页 *
Event-based feature tracking with probabilistic data association;Alex Zihao Zhu et al.;《2017 IEEE International Conference on Robotics and Automation (ICRA)》;20170603;第4465-4470页 *
Low-latency visual odometry using event-based feature tracks;Beat Kueng et al.;《2016 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)》;20161014;第16-23页 *
基于无人机视觉的场景感知方法研究;谢榛;《中国博士学位论文全文数据库 工程科技Ⅱ辑》;20190115(第1期);第C031-73页 *

Also Published As

Publication number Publication date
CN110390685A (zh) 2019-10-29

Similar Documents

Publication Publication Date Title
CN110390685B (zh) 一种基于事件相机的特征点跟踪方法
CN111462200B (zh) 一种跨视频行人定位追踪方法、***及设备
CN109387204B (zh) 面向室内动态环境的移动机器人同步定位与构图方法
Harville Stereo person tracking with adaptive plan-view templates of height and occupancy statistics
Clipp et al. Parallel, real-time visual SLAM
Irani et al. Computing occluding and transparent motions
Saeedi et al. Vision-based 3-D trajectory tracking for unknown environments
CN108492316A (zh) 一种终端的定位方法和装置
CN103854283B (zh) 一种基于在线学习的移动增强现实跟踪注册方法
CN113286194A (zh) 视频处理方法、装置、电子设备及可读存储介质
US11138742B2 (en) Event-based feature tracking
JP4668360B2 (ja) 移動体検出方法及び移動体検出装置
CN109472820B (zh) 单目rgb-d相机实时人脸重建方法及装置
WO2007056711A2 (en) Tracking using an elastic cluster of trackers
US10636190B2 (en) Methods and systems for exploiting per-pixel motion conflicts to extract primary and secondary motions in augmented reality systems
CN106530407A (zh) 一种用于虚拟现实的三维全景拼接方法、装置和***
CN112861808B (zh) 动态手势识别方法、装置、计算机设备及可读存储介质
WO2023134114A1 (zh) 运动目标的检测方法、检测设备以及存储介质
CN114494150A (zh) 一种基于半直接法的单目视觉里程计的设计方法
Dardelet et al. An event-by-event feature detection and tracking invariant to motion direction and velocity
Zhang et al. Target tracking for mobile robot platforms via object matching and background anti-matching
Xue et al. Event-based non-rigid reconstruction from contours
CN103345762B (zh) 基于流形学习的贝叶斯视觉跟踪方法
CN108534797A (zh) 一种实时高精度视觉里程计方法
CN114612545A (zh) 图像分析方法及相关模型的训练方法、装置、设备和介质

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant