[发明专利]一种基于事件相机的特征点跟踪方法

权利要求书

1.一种基于事件相机的特征点跟踪方法，其特征在于包括以下步骤：

第一步，构建基于事件相机的特征点跟踪系统，事件相机的特征点跟踪系统由数据获取模块、初始化模块、事件集选取模块、匹配模块、特征点更新模块、模板边更新模块组成；

数据获取模块与初始化模块、事件集选取模块、模板边更新模块相连；数据集获取模块从事件相机数据集上下载数据，将获取到的图像帧发送给初始化模块，将事件流发送给事件集选取模块，将IMU数据即惯性测量单元数据发送给模板边更新模块；

初始化模块与数据获取模块、事件集选取模块、匹配模块相连；初始化模块从数据获取模块接收图像帧，从图像帧中提取特征点和边图，得到特征点的位置和特征点周围的模板边的位置，将特征点的位置发送给事件集选取模块，将特征点周围的模板边的位置发送给匹配模块；

事件集选取模块与数据获取模块、初始化模块、特征点更新模块、匹配模块相连，事件集选取模块从数据获取模块接收事件流，从初始化模块或特征点更新模块接收特征点的位置，从特征点更新模块接收特征点的光流，在每个特征点周围从事件流中选取特征点的事件集，将每个特征点的位置以及对应的事件集发送给匹配模块；

匹配模块与初始化模块、事件集选取模块、特征点更新模块、模板边更新模块相连，匹配模块从事件集选取模块接收特征点的位置以及每个特征点对应的事件集，从初始化模块或模板边更新模块接收特征点周围的模板边的位置，将特征点的事件集与特征点周围的模板边进行匹配，求出每个特征点的光流，将每个特征点的位置、特征点周围的模板边的位置和特征点的光流发送给特征点更新模块；

特征点更新模块与匹配模块、模板边更新模块、事件集选取模块相连，特征点更新模块从匹配模块接收特征点的位置、特征点周围的模板边的位置、特征点的光流，利用特征点的光流计算特征点新的位置，将特征点周围的模板边的位置以及特征点新的位置发送给模板边更新模块，将特征点的光流和特征点新的位置发送给事件集选取模块，并将特征点新的位置输出；

模板边更新模块与数据获取模块、特征点更新模块、匹配模块相连，模板边更新模块从特征点更新模块接收特征点周围的模板边的位置和特征点新的位置，从数据获取模块接收IMU数据，通过IMU数据来更新特征点周围的模板边的位置，将模板边更新后的位置发送给匹配模块；

第二步，数据获取模块从事件相机数据集中获取图像帧、事件流和IMU数据；

第三步，采用基于事件相机的特征点跟踪系统对从t₀时刻开始到t_N时刻结束的时间段内数据获取模块得到的事件流中的特征点进行跟踪，将时间区间[t₀,t_N]分成一系列的子时间区间[t₀,t₁],...,[t_k,t_k+1],...,[t_N-1,t_N]，N表示子时间区间个数，N的大小由事件流时间长度决定，取值范围为N≥1；第k+1个时间子区间表示为[t_k,t_k+1]；跟踪过程如下：

3.1，令时间序号k＝0；

3.2，若k＝0，初始化模块进行初始化操作，方法为：

3.2.1初始化模块采用Harris角点检测法从数据获取模块得到的图像帧中提取特征点，将提取到的特征点放到集合FD中，令FD为{f₀,...,f_i,...,f_n}，f_i代表检测到的第i个特征点，n为特征点个数；将特征点的位置看做时间的函数，将t₀时刻FD中n个特征点的位置放在特征点位置集合FD₀中，FD₀＝{f₀(t₀),...,f_i(t₀),...,f_n(t₀)}，f_i(t₀)代表特征点f_i在t₀时刻的位置，将FD₀发送给事件集选取模块；

3.2.2初始化模块使用Canny边检测法从数据获取模块得到的图像帧中提取边图，每个图像帧对应一个边图；

3.2.3初始化模块选取FD中n个特征点对应的模板边，将FD中n个特征点对应的t₀时刻模板边的位置集合PBD₀发送给匹配模块，方法为：

3.2.3.1令i＝1；

3.2.3.2以特征点f_i在t₀时刻的位置f_i(t₀)为中心，在f_i(t₀)周围选取一个矩形区域大小为s×s，即该矩形的长和宽均为s；将3.2.2检测到的边图中在内的边作为f_i的模板边，模板边上的像素点即为f_i的模板点；定义内的模板点的集合也即f_i对应的模板边，将放到模板边集合PB中,

p_j为f_i的第j个模板点，将p_j的位置视作时间的函数，p_j(t₀)表示p_j在t₀时刻的位置，表示在t₀时刻p_j在矩形区域内，m_i表示中模板点的个数；

3.2.3.3令i＝i+1；

3.2.3.4若i≤n，转3.2.3.2；否则，表示得到了FD中n个特征点对应的模板边构成的模板边集合PB，定义中模板点在t₀时刻的位置的集合令PB中n个模板边中模板点在t₀时刻的位置集合将PBD₀发送给匹配模块，转第四步；

第四步，事件集选取模块从数据获取模块接收事件流，根据不同的k值分别从初始化模块或特征点更新模块接收不同的数据，在n个特征点周围从事件流中选取特征点的事件集S，将FD_k和特征点的事件集S发送给匹配模块，方法是：

4.1若k＝0，事件集选取模块从数据获取模块接收事件流，从初始化模块接收t_k时刻n个特征点的位置集合FD₀，令t₁＝t₀+1，单位为秒，转4.3；

4.2若k≥1，事件集选取模块从数据获取模块接收事件流，从特征点更新模块接收t_k时刻n个特征点的位置集合FD_k和t_k-1时刻n个特征点的光流集合G_k-1，预估子时间区间[t_k,t_k+1]的大小，根据公式(2)计算t_k+1：

其中表示在子时间区间[t_k-1,t_k]上特征点f_i的光流；表示在子时间区间[t_k-1,t_k]上所有特征点的平均光流，转4.3；

4.3在t_k时刻n个特征点的位置集合FD_k中每个位置的周围选取每个特征点对应的事件集，方法是：

4.3.1令i＝1；

4.3.2从事件流中选出符合公式(3)要求的事件，将这些事件放到特征点f_i在t_k时刻对应的事件集中，并将放到事件集S中：

表示在一个三维时空区域内的事件的集合，空间范围为特征点f_i在t_k时刻的位置f_i(t_k)周围的矩形区域时间范围为区间[t_k,t_k+1]；e_d代表中的第d个事件，d＝{1,2,...,z_i}，z_i表示内事件的个数；x_d表示事件e_d在像素坐标系中的坐标，表示事件e_d发生的时间，表示e_d的像素坐标在内；

4.3.3令i＝i+1；

4.3.4若i≤n，转4.3.2；否则，说明得到了t_k时刻FD中n个特征点对应的事件集S，将FD_k和S发送给匹配模块，转第五步；

第五步，匹配模块从事件集选取模块接收FD_k和S，若k＝0，从初始化模块接收n个特征点对应的模板边的位置集合PBD_k，转第六步；若k≥1，从模板边更新模块接收n个特征点对应的模板边的位置集合PBD_k，转第六步；

第六步，匹配模块将S与特征点周围的模板边进行匹配，求出t_k时刻n个特征点的光流集合G_k，将n个特征点的位置集合FD_k、n个特征点对应的模板边的位置集合PBD_k和光流集合G_k发送给特征点更新模块，具体方法为：

6.1令i＝1；

6.2对特征点f_i构造匹配误差函数，方法是：

按公式(4)修正事件e_d在t_k时刻的位置：

x′_d表示计算得到的事件e_d在t_k时刻的位置，简称为事件e_d运动修正的位置，表示在时间区间[t_k,t_k+1]上特征点f_i的光流；符号·表示点乘，定义符号

构造匹配误差函数如下：

e表示误差，p_j(t_k)代表t_k时刻模板点p_j的位置；r_dj表示e_d由模板点p_j生成的概率，即e_d与p_j匹配的概率；公式中双竖线表示对双竖线内的向量进行取模操作，r_dj的计算公式如下：

6.3采用EM-ICP算法最小化匹配误差函数，求解得到最优的光流方法是：

6.3.1初始化

6.3.2通过公式(6)计算r_dj；

6.3.3更新光流：

6.3.4计算光流的变化量并令将放到t_k时刻光流集合G_k中；

6.3.5若表示得到最终的优化结果转6.4；若转6.3.2；σ为光流的变化量阈值；

6.4令i＝i+1；

6.5若i≤n，转6.2；否则说明得到了t_k时刻光流集合G_k，将FD_k、PBD_k和G_k发送给特征点更新模块，转第七步；

第七步，特征点更新模块从匹配模块接收FD_k、PBD_k和G_k，通过光流计算t_k+1时刻n个特征点的位置集合FD_k+1，将G_k和FD_k+1发送给事件集选取模块，将FD_k+1和PBD_k发送给模板边更新模块，方法是：

7.1令i＝1；

7.2计算t_k+1时刻f_i的位置f_i(t_k+1)，将f_i(t_k+1)放到集合FD_k+1中，

为光流与时间相乘，表示从t_k时刻到t_k+1时刻特征点f_i运动距离；

7.3令i＝i+1；

7.4若i≤n，转7.2；否则说明得到了t_k+1时刻n个特征点的位置，并得到集合FD_k+1，FD_k+1＝{f₁(t_k+1),...,f_i(t_k+1),...,f_n(t_k+1)}；将G_k和FD_k+1发送给事件集选取模块，将FD_k+1和PBD_k发送给模板边更新模块；并将t_k+1时刻n个特征点的位置集合FD_k+1进行显示或存储到结果文件中，转第八步；

第八步，模板边更新模块从数据获取模块接收IMU数据，从特征点更新模块接收FD_k+1和PBD_k，使用IMU数据对PBD_k中的位置进行更新，得到t_k+1时刻n个特征点对应的模板边的位置集合PBD_k+1，将PBD_k+1发送给匹配模块，方法如下：

8.1令i＝1；

8.2更新特征点f_i对应的模板边在t_k+1时刻的位置，方法是：

8.2.1令j＝1；

8.2.2定义符号F为与f_i相对应的三维空间中的点，P_j为与模板点p_j相对应的三维空间中的点，F和P_j均用三维坐标(x,y,z)的形式表示；将模板点p_j使用齐次坐标形式表示，格式为(x,y,1)，前两维分别表示p_j在像素坐标系下的横坐标和纵坐标，得到如下方程：

在t_k时刻，

在t_k+1时刻，

其中K为事件相机的内参矩阵，为事件相机自带参数；R为从t_k时刻到t_k+1时刻事件相机的旋转矩阵，t为从t_k时刻到t_k+1时刻事件相机的平移向量，均通过获取到的IMU数据计算得到；为P_j在t_k时刻相机坐标系下的深度、为P_j在t_k+1时刻相机坐标系下的深度，为F在t_k时刻相机坐标系下的深度、为F在t_k+1时刻相机坐标系下的深度；

8.2.3将公式(10)中的两个等式相减，得到t_k+1时刻模板点p_j相对于f_i的相对位置

P_j和F在t_k+1时刻的相机坐标系中具有相同的深度，即将公式(9)带入公式(11)，并对公式(11)进行化简得到：

对公式(12)化简得到t_k+1时刻模板点相对位置的公式为：

符号Nor()表示齐次化操作，即将括号内的坐标转化为齐次坐标；按公式(14)得到t_k+1时刻模板点p_j的位置p_j(t_k+1)，将p_j(t_k+1)放到t_k+1时刻特征点f_i周围模板边的位置集合

8.2.4令j＝j+1；

8.2.5如果j≤m_i，转8.2.2；否则表示得到了将放到t_k+1时刻n个特征点周围模板边的位置集合PBD_k+1中，转8.3；

8.3令i＝i+1；

8.4若i≤n，转8.2；否则说明得到PBD_k+1，将PBD_k+1发送给匹配模块，转第九步；

第九步，令k＝k+1；

第十步，若kN，转第四步；否则结束。