Contrast Maximization-Based Feature Tracking for Visual Odometry with an Event Camera

Abstract

As a new type of vision sensor, the dynamic and active-pixel vision sensor (DAVIS) outputs image intensity and asynchronous event streams in the same pixel array. We present a novel visual odometry algorithm based on the DAVIS in this paper. The Harris detector and the Canny detector are utilized to extract an initialized tracking template from the image sequence. The spatio-temporal window is selected by determining the life cycle of the asynchronous event streams. The alignment on timestamps is achieved by tracking the motion relationship between the template and events within the window. A contrast maximization algorithm is adopted for the estimation of the optical flow. The IMU data are used to calibrate the position of the templates during the update process that is exploited to estimate camera trajectories via the ICP algorithm. In the end, the proposed visual odometry algorithm is evaluated in several public object tracking scenarios and compared with several other algorithms. The tracking results show that our visual odometry algorithm can achieve better accuracy and lower latency tracking trajectory than other methods.

1. Introduction

Visual odometry has played an important role in robot navigation [1], intelligent transportation [2], and intelligent inspection [3]. Although several decades of active research have led to a certain level of maturity, we still face challenges in scenes with high dynamics, low texture, or harsh lighting conditions [4]. Conventional visual odometry generally acquires images from frame-based cameras and estimates camera motion from several adjacent images [5]. However, for objects with a high dynamic range or high-speed motion, frame-based cameras cannot obtain clear images. Therefore, it is difficult for a frame-based camera to extract the feature points of the image and estimate the camera pose in this case. Furthermore, during the blind period between frames, the frame-based camera does not capture the precise motion of the feature information. Frame-based cameras capture redundant information in static scenes, resulting in a waste of storage and computational resources [6].

Bio-inspired event cameras, such as the dynamic vision sensor (DVS) [7], overcome the above-mentioned limitations of frame-based cameras [8]. As a new type of bio-heuristic vision sensor, the event camera has a completely different mode from the conventional camera [9], and event cameras only output the change in pixel-level brightness. For each pixel, when the intensity change reaches a certain threshold, an event is triggered; the event carries information on pixel coordinates, timestamps, and polarity [10]. Compared with the low frame rate, large delay, and insufficient dynamic response range in conventional visual cameras, an event camera has the characteristics of fast response speed, high dynamic range, and low power consumption [11]. Although traditional frame-based tracking algorithms have been vigorously developed, asynchronous event streams cannot be directly handled by current frame-based pose estimation methods [12]. As a result, event-based tracking algorithms are required. The dynamic and active-pixel vision sensor (DAVIS) is a sensor that combines a frame-based camera with an asynchronous event-based camera in the same pixel array. Based on the assumption that events are mainly triggered by high gradient edges in the image, the optimal motion parameters for events can be computed by maximizing contrast in the Image of Warped Events (IWE) [13]. Various reward functions for evaluating contrast were proposed and analyzed in recent work by Gallego et al. [14]. Contrast maximization algorithms for events have also been successfully used to solve various problems of event cameras, such as optical flow estimation [15,16], motion segmentation [17,18], 3D reconstruction [19], and motion estimation [20]. Since event streams cannot provide absolute brightness values and synchronized output with image frames, a contrast maximization algorithm is utilized to resolve event and image frame data associations in this paper, while the event-based feature tracking can be further simplified.

In this paper, we present a novel visual odometry algorithm based on template edges using a DAVIS camera. In our tracking approach, we leverage a combination of event- and frame-based camera measurements. The tracked initial tracking template is extracted by the feature detector in the image, and the spatio-temporal windows are selected by determining the life cycle of the asynchronous event streams. The calibrated tracking templates are utilized to calculate the event camera trajectory. The main contributions are summarized as follows:

• Compared with traditional event-by-event tracking methods [20,21,22], a new tracking mechanism is presented to resolve data associations. A contrast maximization method is adopted to calculate the displacement parameters of the events, and the IMU data are used to calibrate the rotation parameters of the events, which greatly enhances the calculation speed and accuracy of the event stream.
• Since the ICP algorithm is highly dependent on the depth of the scene, a robust Beta-Gaussian distribution depth filter is presented to obtain a more accurate depth of the tracking template than depth estimation with only triangulation [23,24].
• We successfully apply our method to evaluation experiments in several different scenarios on public event camera datasets. Compared with the visual odometry algorithm [7,9], the proposed algorithm can achieve better performance and obtain a lower-latency camera trajectory.

The rest of this article is organized as follows. The related work on event-based camera visual odometry is presented in Section 2. The realization steps of fusion events and image frame tracking are described in Section 3. The effectiveness of the algorithm is verified on several DAVIS datasets in Section 4. Finally, conclusions are given in Section 5.

2. Related Work

Feature detection and tracking methods for frame-based cameras are well established. However, these frame-based methods cannot track the blind time between two adjacent frames, and frame-based cameras still capture information in all pixel arrays, even in scenes without motion [17]. In contrast, an event-based camera only acquires information about the areas of the scene where the intensity value has changed, and it fills the blind spots between adjacent frames with a higher asynchronous response rate. The advantages of event-based cameras are better suited for applications such as driverless vehicles and motion tracking [22,25].

Event-based visual odometry can be divided into two main methods: one is to use traditional methods to accumulate event information for tracking, and the other is to track directly based on asynchronous event streams. Gallego et al. [26] presented the pose tracking of a 6-DoF event camera from an existing photometric depth map (intensity + depth information). Stoffregen et al. introduced event accumulation frames with edge contrast maximization for motion segmentation [27]. Alzugaray and Chli [28] proposed a purely event-based corner detector and a new corner tracker, proving that it is possible to directly detect and track corners in event streams. Although there is a great deal of research devoted to event-based feature detection, very little work has been done to consider the problem of event tracking. Some approaches for localization and mapping with event cameras have several similarities with our method. Gallego et al. [29] detected independent moving objects by tracking corners detected in event images integrated over short time windows. Mueggler et al. [30] proposed a DVS-based ego-motion estimation method that uses a continuous-time framework to directly integrate the information transmitted by the camera. The pose trajectory of the DVS is estimated based on the observed events. Kim et al. [31] estimated 6-DoF camera motion, log intensity gradient, and inverse depth and used three decoupled probabilistic filters in real time.

Recently, Mueggler et al. [32] extended a frame and event tracker for DAVIS cameras and integrated events from event cameras for high-speed object tracking. In contrast to the above-mentioned approach, we consider the connection between event stream and images in our modeling process, and higher-rate trajectory tracking is provided.

3. Main Methods

This paper proposes a visual odometry algorithm based on DAVIS. The main methods are divided into two parts: feature detection and feature tracking. As shown in Figure 1, the entire process was divided into six steps. We detect features in the image sequence. The contrast maximization algorithm is then utilized to match the event stream in the corresponding spatio-temporal windows and calculate the optical flow of the event. The estimated motion parameters and the IMU measurement are exploited to calibrate the tracking template. We then detect the depth values of the tracking template using a depth filter. Finally, the 6-DoF event camera pose can be estimated by the ICP algorithm.

3.1. Feature Detection

Since events in edge regions of the scene are triggered more frequently than those in low-texture regions, we devise suitable features for tracking. As shown in Figure 2, the method first extracts feature points and edge maps from image frames by the Harris detector [9] and the Canny detector [20] in the feature detection stage. Then, it selects the edge in the specified area around each feature point as the template edge of the feature point. All regions are square in the same size, which is a tunable parameter. Our method does not need to provide frames at a constant rate, as they are only used to initialize features. Keyframes can be added to replace missing features.

3.2. Feature Tracking

Assuming that feature point $x$ is detected from the image frame at time $t_0$, the motion of feature point $x$ can be described as follows:

$$x\left(t\right)=x\left(t_0\right)+{\displaystyle {\int }_{t_0}^t\dot{x}\left(s\right)}ds$$

(1)

where $\dot{x}\left(s\right)$ represents the position differential of the feature point $x$ at time $s$.

A set of events is selected from the event stream, the time of feature point $x$ is selected as the spatio-temporal window corresponding to initial time $t_0$, and $W$ is the set of events in the spatio-temporal window.

$$W={\left\{e_i|t_0<t_{e_i}<t_1\right\}}_{i=1}^n$$

(2)

Here, $e_i$ represents the $i-th$ event in the spatio-temporal window, and $t_{e_i}$ represents the time when event $e_{\mathrm{i}}$ occurs. $n$ represents the number of events. Since $\left[t_0,t_1\right]$ is the first sub-time interval, the value of $t_1$ can be calculated by setting the size of the spatio-temporal window. For example, let $W$ contain 10,000 events, that is, $n=$ 10,000.

3.2.1. Choice of Spatio-Temporal Windows

In order to achieve asynchrony of the feature point tracking method, the size of the sub-time interval is determined by the method during real-time operation. The faster the scene moves, the smaller the sub-time interval and the faster the update frequency of feature points.

The specific calculation process is as follows. After the optical flow of all feature points in this iteration is obtained, the size of the next sub-time interval is calculated from the optical flow. We define $x_{\mathrm{i}}$ to represent the $i-th$ feature point, $i=\left\{1,\cdots ,m\right\}$, and the number of feature points is $m$. ${\theta }_i^n$ is the optical flow of feature point $x_{\mathrm{i}}$ in the $n-th$ sub-time interval $\left[t_{n-1},t_n\right]$. Given the optical flow of all feature points ${\left\{{\theta }_i^n\right\}}_{i=1}^m$ in the sub-time interval $\left[t_{n-1},t_n\right]$, the next sub-time interval $\left[t_n,t_{n+1}\right]$ can be calculated:

$$\left\{\begin{array}{l}{\theta }_{average}^n=({\displaystyle \sum _{i=1}^m}{\theta }_i^n)/m\\ t_{n+1}=t_n+3/{\theta }_{average}^n\end{array}\right.$$

(3)

where the unit of the number 3 is pixels, and ${\mathbf{\theta }}_{average}^n$ represents the average optical flow of all feature points in the sub-time interval $\left[t_{n-1},t_n\right]$. The $t_{n+1}$ value is calculated through the formula, and the time required for the feature points in the previous time interval to move by 3 pixels on average is used as the estimated value of the current interval $\left[t_n,t_{n+1}\right]$.

3.2.2. Maximizing IWE Contrast

After the spatio-temporal window is determined by corresponding to the feature point, we use the contrast maximization algorithm to match the event set $W$ around the feature point $x$ with the template point set. It is assumed that all template points have the same optical flow as feature point $x$ (the optical flow $\theta$ of all pixels in the region is the same), and the optical flow of feature point $x$ is constant in the sub-time interval $\left[t_0,t_1\right]$. Let the optical flow $\theta$ of the feature point $x$ in the time interval $\left[t_0,t_1\right]$ be defined as $v$. For event $e_{\mathrm{i}}$ in $W$, as shown in Figure 3, Image of Warped Events (IWE) is used to calculate its position $X_k^{\prime }$ at time $t_0$. The formula is as follows:

$$X_k^{\prime }=\left[\begin{array}{c}{x}_k^{\prime }\\ y_k^{\prime }\end{array}\right]=\left[\begin{array}{c}{x}_k\\ y_k\end{array}\right]+\left[\begin{array}{c}{v}_1\\ v_2\end{array}\right]•\left(t_0-t_{\mathit{ref}}\right)$$

(4)

The weighted IWE is defined as:

$$I_j\left(x\right)=\sum _{k=1}^{N_e}{P}_{kj}\delta (X-X_{kj}^{\prime })$$

(5)

where $X_{kj}^{\prime }$ is the position of the $k-th$ event after it is distorted along the $j-th$ optical flow $\theta$, $N_e$ represents the number of events, $\delta$ represents the Dirac function, $P_{kj}$ represents the probability that the $k-th$ event belongs to the $j-th$ optical flow, and $I_j\left(x\right)$ represents the IWE corresponding to the $j-th$ optical flow.

Events are aligned by image contrast, which is defined by a sharpness/dispersion metric, such as the variance:

$$Var\left(I_j\right)={\displaystyle {\int }_{\Omega }{\left(I_j\left(x\right)-{\mu }_j\right)}^{2}dx}$$

(6)

where $\Omega$ is the image plane and ${\mu }_j$ is the mean of the $j-th$ Image of Warped Event.

$$\theta \leftarrow \theta +\mu {\nabla }_{\theta }\left({\displaystyle \sum _{j=1}^{N_l}Var\left(I_j\right)}\right)$$

(7)

where, $\mu$ represents the step size, and $N_l$ represents the number of clusters.

3.2.3. Template Edge Update

After the optical flow of the feature points is obtained, the feature points and template edges are updated. However, when the camera is rotated, the template edges and template points move at significantly different speeds, and the points farther from the center of rotation move faster. Therefore, if the template points are updated by only using optical flow, the position of feature points will quickly deviate from the true value due to the existence of rotation factors. The process of updating the stencil edge is divided into two steps. First, the position of the stencil edge is updated by using the optical flow, and then the position of the stencil edge is corrected by using the IMU data.

The optical flow $\theta$ is used to update the position of feature point $x$ and the corresponding position of the template edge $x_j$. Assuming that the optical flow of the feature point is constant in the sub-time interval $\left[t_0,t_1\right]$, the optical flow $\theta$ is then utilized to update the position of feature point $x$; the formula is as follows:

$$x\left(t_1\right)=x\left(t_0\right)+\theta ·\left(t_j-t_0\right)$$

(8)

$x\left(t_0\right)$ is the position of the initially extracted feature point, $x\left(t_j\right)$ is the position of the updated feature point, and $t_j$ is the time of the tracking template in the sub-time interval $\left[t_0,t_1\right]$.

In order to eliminate the influence of rotation on the template edge update, we introduce IMU data to correct the template edge position. As a result, the corrected position is closer to the true value of the template edge position. The relative position of the template points relative to the feature point is calculated at time $t_j$, and then the IMU data are used to correct the relative position. For the template point $x_j$, we define its relative position:

$$x_j^{relative}=x_j\left(t\right)-x\left(t\right)$$

(9)

Since the data in the IMU coordinate system and the tracking template from the camera coordinate system are not in the same system, we need to transform the data. The IMU measurements contain the accelerometer measures and gyroscope measures for each axis, but what we need during the update process is the rotation matrix between the templates. So, we convert linear acceleration and angular velocity into Euler angles and then convert the Euler angles into a rotation matrix. The accelerometer can calculate the roll angle ${\alpha }_{acc}$ and pitch angle ${\beta }_{acc}$ at the time of rest. The gyroscope is the angular velocity integral in the time interval and can calculate three angles: the roll angle ${\alpha }_{gyro}$, pitch angle ${\beta }_{gyro}$, and yaw angle ${\gamma }_{gyro}$. The Euler angle fusion formula can be expressed as:

$$\left\{\begin{array}{c}\alpha ={\alpha }_{\mathrm{g}yro}+\left({\alpha }_{acc}-\right.{\alpha }_{\mathrm{g}yro})\cdot k\hfill \\ \beta ={\beta }_{\mathrm{g}yro}+\left({\beta }_{acc}-{\beta }_{\mathrm{g}yro}\right)\cdot k\hfill \\ \gamma ={\gamma }_{\mathrm{g}yro}\hfill \end{array}\right.$$

(10)

where $\alpha$, $\beta$, and $\gamma$ are obtained after the complementary fusion of the accelerometer and gyroscope pose; $k$ represents the proportional coefficient, which needs to be adjusted according to the actual situation, such as by choosing 0.3. The conversion of Euler angles to a rotation matrix is expressed as:

$$R_{wi}=\left[\begin{array}{ccc}\cos \gamma & -\sin \gamma & 0\\ \sin \gamma & \cos \gamma & 0\\ 0& 0& 1\end{array}\right]\ast \left[\begin{array}{ccc}\cos \beta & 0& \sin \beta \\ 0& 1& 0\\ -\sin \beta & 0& \cos \beta \end{array}\right]\ast \left[\begin{array}{ccc}1& 0& 0\\ 0& \cos \alpha & -\sin \alpha \\ 0& \sin \alpha & \cos \alpha \end{array}\right]$$

(11)

where $R_{wi}$ is the rotation matrix in the IMU coordinate system.

Then, we solve the transformation of the rotation matrix between the camera coordinate system and the IMU coordinate system. The transformation relation is denoted as:

$$R_{wc}=R_{wi}\cdot R_{ic}$$

(12)

where $R_{wc}$ is the rotation matrix in the camera coordinate system, and $R_{ic}$ represents the transformation matrix between the IMU and camera coordinate systems during camera calibration. The IJRR dataset takes the first frame of the camera as the world coordinate system. We assume that the pose transformation of the camera is linear. The rotation matrix of the tracking template can be obtained by linear interpolation of the time scale. The process is represented as follows:

$$R_j=inv\left(R_{wc}^{t_0}\right)\cdot R_{wc}^{t_1}\cdot \frac{t_j-t_1}{t_0-t_1}$$

(13)

where $R_{wc}^{t_0}$ is the rotation matrix in the camera coordinate system at time $t_0$. $R_{wc}^{t_1}$ is the rotation matrix in the camera coordinate system at time $t_1$, and it is calculated from the rotation matrix $R_{wi}^{t_1}$ in the IMU coordinate system at time $t_1$ and the transformation matrix $R_{ic}$ between the camera and IMU; $R_j$ is the interpolated rotation matrix between $t_0$ and $t_j$.

The symbols $X_j$ and $X$ respectively represent the 3D coordinates corresponding to the template point $x_j$ and the feature point $x$ in the camera coordinate system at time $t_0$.

At time $t_0$:

$$\left\{\begin{array}{c}{s}_{x_j}^0{x}_j\left(t_0\right)=K{X}_j\\ s_x^{0}x\left(t_0\right)=KX\end{array}\right.$$

(14)

At time $t_j$:

$$\left\{\begin{array}{l}{s}_{x_j}^1{x}_j\left(t_j\right)=K\left(R_j{X}_j+t\right)\hfill \\ s_x^{1}x\left(t_j\right)=K\left(R_{j}X+t\right)\hfill \end{array}\right.$$

(15)

where $t$ is the translation vector between $t_0$ and $t_j$, $R_j$ is the interpolated rotation matrix between $t_0$ and $t_j$, and $K$ represents the internal parameter matrix of the camera.

Substituting the above formula into (9) and normalizing the 3D coordinates, we obtain:

$$x_j^{relative}\left(t_j\right)=Nor\left(K{R}_j{K}^{-1}{x}_j\left(t_0\right)\right)-Nor\left(K{R}_j{K}^{-1}x\left(t_0\right)\right)$$

(16)

At this time, the position of $x_j$ at time $t_j$ is obtained by adding the position of the feature point $x$ and the relative position of $x_j^{relative}$ at time $t_j$:

$$x_j\left(t_j\right)=x_j^{relative}\left(t_j\right)+x\left(t_j\right)$$

(17)

Finally, the updating of all template points can be completed through the formula.

3.2.4. Depth Estimation

After the corresponding position of the edge map of each frame is obtained, the 3D position coordinates of the edge in the space can be recovered through the matching relationship of the pixels on the edge between different edge maps. The triangulation method is used, and the coordinates of the 3D point in space are recovered by using the pixel positions of the 3D points observed from different viewing angles. The same 3D point can be observed in multiple frames, and the corresponding spatial 3D point coordinates can be calculated for any two frames. The strategy of our method is using the properties of the Gaussian distribution to fuse the current observations, assuming that the depth values follow a Gaussian distribution. In the final estimation of the same 3D point, the strategy is based on a uniform Gaussian mixture distribution.

First, a triangulation method is used to recover the depth of pixels in the edge map. Assuming that a certain pixel $p_0$ of the edge graph and pixel point $p_1$ of the edge graph are a pair of matching points, and the two pixels correspond to the same 3D point $P$ in the space, the formula holds as follows:

$$s_0{K}^{-1}{p}_0=s_{1}R{K}^{-1}{p}_1+t$$

(18)

where $R$ and $t$ are the rotation matrix and translation vector of the event camera from the edge graph to the edge graph, respectively. $s_0$ and $s_1$ represent the depths of the 3D point $P$ in the coordinate system of the event camera at times $t_i$ and $t_j$, respectively. Equation (18) is further simplified to obtain:

$$s_1=\frac{-\left(K^{-1}{p}_0\right)^t}{\left(K^{-1}{p}_0\right)^R{K}^{-1}{p}_1}$$

(19)

Then, a depth value $y_i$ can be calculated from the corresponding points in the two images. The distribution of $y$ can be jointly represented by a Gaussian distribution and a uniform distribution:

$$p\left(y_i|\widehat{y},\pi \right)=\pi N\left(y_i|\widehat{y},{\tau }_i^2\right)+\left(1-\pi \right)U\left(y_i|y_{min},y_{max}\right)$$

(20)

where $N\left(y_i|\widehat{y},{\tau }_i^2\right)$ is a Gaussian distribution centered on the true value $\widehat{y}$, and ${\tau }_i^2$ is its variance. $\pi$ is the probability of estimating the correct depth; the closer the depth measurement is to the true value, the closer $\pi$ is to 1. As shown in Figure 4, $y_{min}$ and $y_{\mathit{max}}$ are the upper and lower limits of the uniform distribution.

Equation (20) can be approximated by a Beta-Gaussian distribution. The step 4 of Algorithm 1 summarizes the depth estimation steps of our algorithm.

$$p(Z,\pi |a_i,b_i,{\mu }_i,{\sigma }_i)=\mathrm{Beta}(\pi |a_i,b_i\left)N\right(Z|u_i,{\sigma }_i^2)$$

(21)

3.2.5. Pose Estimation

According to Algorithm 1, we can calculate corresponding 3D point sets $P=\left\{p_1,p_2,\cdots ,p_n\right\}$ and $Q=\left\{q_1,q_2,\cdots ,q_n\right\}$ through Section 3.2.3 and Section 3.2.4, where the number of corresponding point pairs is $n$. Then, the optimal coordinate transformation is calculated by the ICP algorithm, that is, the rotation matrix $R$ and the translation vector $t$; the problem can be described as follows:

$$(R,t)=\underset{R\in SO\left(3\right),t\in {ℝ}^3}{\arg \min }{\displaystyle \sum _{i=1}^n{w}_i}\parallel (R{p}_i+t)-q_i{\parallel }^2$$

(22)

where $w_i$ represents the weight of each point. $R$ and $t$ are our required rotation matrix and translation vector. To reduce the reprojection error during tracking, the Tukey weight function is used as follows:

$$w_i=\left\{\begin{array}{cc}{\left(1-\frac{x^2}{b^2}\right)}^2& \left|x\right|\le \left|b\right|\\ 0& otherwise\end{array}\right.$$

(23)

where $x=\Vert p_w^{-}-p_w\Vert$ and $b=5$ pixels. Additionally, we find that with $w_i$ multiplied by the internal probability $\pi$, the estimated results are significantly better. Template points that are well tracked get the highest weight, while template points whose matching does not converge are usually removed due to their errors being too large.

Algorithm 1: Pseudo-code of pose estimation.

Step 1: Extracting feature points and edge points

1: Extract feature points by Harris detector

2: Extract edge points by Canny detector

3: Pick edge points around feature points as template edges

Step 2: Computing the optical flow and the spatio-temporal windows of events

4: Input: template edges P and event point set W

Output: optical flow $\theta$

Initialize optical flow $\theta =$ 0

$\mathrm{for} i=1:m$ do

for $j=1:n$ do

$\hspace{1em}I\left(x\right)={\displaystyle \sum _{j=1}^n}\delta (X_i-X_j^{\prime });\mathit{Var}\left(I\right)={\displaystyle {\int }_{\Omega }{\left(I\left(x\right)-\mu \right)}^2\mathit{dx};}\theta \leftarrow \theta +\mu {\nabla }_{\theta }\left({\displaystyle \sum _{k=1}^{N_l}{\mathit{Var}(I}_k)}\right)$

end for

end for

5: Update the spatio-temporal window: $t_{n+1}=t_n+3/{\theta }_{average}^n$

Step 3: Updating the template edges using the optical flow and IMU data

6: Update the postion of feature points: $x_i\left(t_1\right)=x_i\left(t_0\right)+v.\left(t_1-t_0\right)$

7: Update the relative position: $x_j^{relative}\left(t_1\right)=Nor\left(K{R}_j{K}^{-1}{x}_j\left(t_0\right)\right)-Nor\left(K{R}_j{K}^{-1}{x}_i\left(t_0\right)\right)$

8: Update the postion of edge points: $x_j\left(t_1\right)=x_j^{relative}\left(t_1\right)+x\left(t_1\right)$

Step 4: Computing the depth value of template edges using a depth filter

9: Triangulation depth: $y_i$

10: Depth filter: $q(y,\pi |a_i,b_i,{\mu }_i,{\sigma }_i)=\mathrm{Beta}(\pi {|a}_i,b_i\left)N\right(y\mid u_i,{\sigma }_i^2)$

Input: Triangulation depth: $y_i$

Output: depth value $y$

$\hspace{1em}\mathrm{for} i=1:n$ do

$\hspace{1em}\hspace{1em}{\mu }_i^{\prime }=C_1^{\prime }{m}_i+C_2^{\prime }{\mu }_i$; ${{\sigma }_i^{\prime }}^2=C_1^{\prime }(y_i^2+m_i^2)+C_2^{\prime }({\sigma }_i^2+{\mu }_i^2)$

$\hspace{1em}\hspace{1em}{f}_i=C_1^{\prime }\frac{a_i + 1}{a_i + b_i + 1} + C_2^{\prime }\frac{a_i}{a_i + b_i + 1}$

$\hspace{1em}\hspace{1em}{e}_i=C_1^{\prime }\frac{\left(a_i + 1\right)\left(a_i + 2\right)}{\left(a_i + b_i + 1\right)\left(a + b_i + 2\right)} + C_2^{\prime }\frac{a_i(a_i + 1)}{\left(a_i + b_i + 1\right)\left(a_i + b_i + 2\right)}$

$\hspace{1em}\hspace{1em}{a}_i^{\prime }=\frac{e_i-f_i}{f_i-\frac{e_i}{f_i}}$; $b_i^{\prime }=\frac{1-f_i}{f_i}\cdot a_i^{\prime }$; $\frac{1}{y_i^2}=\frac{1}{{\sigma }_i^2}+\frac{1}{{\tau }_i^2}$; $m_i=y_i^2(\frac{{\mu }_i}{{\sigma }_i^2}+\frac{x_i}{{\tau }_i^2})$

end for

Step 5: Pose estimation using the ICP agorithm

11: ICP agorithm: $(R,t)=\underset{R\in SO\left(3\right),t\in {ℝ}^3}{\arg \min }{\displaystyle \sum _{i=1}^n{w}_i}\parallel (R{p}_i+t)-q_i{\parallel }^2$

Output

$6-\mathrm{DoF} \mathrm{event} \mathrm{camera} \mathrm{pose} \mathrm{translation}\in R^3$, $\mathrm{rotation} \in SO\left(3\right)$.

4. Experiments

In order to demonstrate the performance of our proposed algorithm, some widely used event datasets were utilized to estimate camera trajectories. These datasets were generated using the DAVIS240C from iniLabs, and they contain events, images, and IMU measurements. The algorithm was evaluated by selecting some challenging sequences from the IJRR event camera datasets [33]. The sequences are from several indoor and outdoor scenes with different lighting conditions. Experiment 1 and Experiment 2 were evaluated indoors. Compared with Experiment 1, Experiment 2 had richer scene textures. In Experiment 3, the camera trajectory was evaluated outdoors. In several experiments, the camera trajectories obtained by our method were compared with frame-tracking-based (ORB) visual odometry and event-tracking-based (EVO) visual odometry, and we also compared the trajectories of our method with and without IMU data. Among the methods tested, only ours uses IMU data to calibrate rotation. The absolute pose error was calculated between the poses estimated by other algorithms and our method and the ground truth.

In the IJRR datasets [33], we selected representative scenes for trajectory estimation. We tested different visual odometry methods on three sequences: shapes_6dof, boxes_6dof, and outdoors_6dof. Figure 5, Figure 6 and Figure 7 show the test results of the different algorithms on the three sequences. These figures show the trajectories produced by visual odometry such as Our_method, ORB, and EVO.

Table 1. Absolute pose error (m/s), position error (m/s), and rotation error (deg/s) test results in different scenarios with different methods [33].

	Sequences	Our_Method	ORB	Our_Method (+IMU)	EVO
Absolute pose error	shapes_6dof	0.0815	0.0780	0.0315	0.0435
	boxes_6dof	0.0434	0.0369	0.0204	0.0344
	outdoors_6dof	0.1416	0.2357	0.0461	0.0227
Position error	shapes_6dof	0.1562	0.0845	0.0391	0.0496
	boxes_6dof	0.2255	0.0956	0.0436	0.0586
	outdoors_6dof	0.1728	0.0735	0.0325	0.0265
Rotation error	shapes_6dof	2.4562	1.8342	0.6443	0.9093
	boxes_6dof	3.5432	1.6901	0.5426	0.5084
	outdoors_6dof	3.4734	1.6042	0.4453	0.7253

shows the position error, orientation error, and absolute pose error of different methods compared to the ground truth in different scenarios. The position error was calculated as the Euclidean distance between the estimated value of the camera trajectory and the true value. Rotation error was computed as the angular error based on the estimated and true values of the rotation matrix. The absolute pose error compares the estimated trajectory and the reference trajectory and calculates the statistics of the entire trajectory, which is suitable for testing the global consistency of the trajectory. Furthermore, to reflect the contribution of the IMU to the estimation calibration, we obtained the results of our method on the shapes_6dof, boxes_6dof, and outdoors_6dof sequences without fusing the IMU data.

In Experiment 1 (Figure 5), we tracked simple shapes using different visual odometry methods. As the camera movement speed increases, the resolution of the captured image gradually decreases. The ORB method does not extract enough feature points, which leads to tracking failure. Thus, we chose a tracking duration of only 35 s. The position errors of the different methods were 15.62 cm, 8.45 cm, 3.91 cm, and 4.96 cm for Our_method, ORB, Our_method(+IMU), and EVO, respectively. The rotation errors of the different methods were 2.46 deg, 1.83 deg, 0.64 deg, and 0.91 deg for Our_method, ORB, Our_method(+IMU), and EVO, respectively. In Experiment 2 (Figure 6), the natural textures in the box scene are richer and also generate more event information. Since the EVO algorithm tracks the camera trajectory through pure event information, too many event streams cause the tracking time to be longer. The position errors of the different methods were 22.55 cm, 9.56 cm, 4.36 cm, and 5.86 cm for Our_method, ORB, Our_method (+IMU), and EVO, respectively. The rotation errors of the different methods were 3.54 deg, 1.69 deg, 0.54 deg, and 0.51 deg for Our_method, ORB, Our_method (+IMU), and EVO, respectively. In Experiment 3 (Figure 7), we evaluated several methods on outdoor datasets. The position errors of the different methods were 17.28 cm, 7.35 cm, 3.25 cm, and 2.65 cm for Our_method, ORB, Our_method (+IMU), and EVO, respectively. The rotation errors of the different methods were 3.47 deg, 1.60 deg, 0.45 deg, and 0.73 deg for Our_method, ORB, Our_method (+IMU), and EVO, respectively. The mean scene depth was 3 m. Overall, our algorithm performed slightly better than EVO in tracking performance and far outperformed ORB, but our tracking time was much faster than that of EVO in complex textured scenes. Our algorithm provides low-latency pose updates and preserves the nature of event-based data. These results represent the success of fusing frames and events for tracking in natural scenes and 6-DoF motion.

5. Conclusions

In this paper, we presented a novel visual odometry algorithm based on the DAVIS. The initial tracking template was extracted from the image sequence. A contrast maximization algorithm was presented to estimate the optical flow estimated by minimizing the distance between the event and the templates in the spatio-temporal windows. Then, IMU data were used to calibrate the rotational position of the tracking template. A Beta-Gaussian distribution depth filter was presented to update the depth of each pixel on the edge of templates. These templates were used to achieve lower-latency camera trajectories by the ICP algorithm. We tested our method on several scenes with different textures in the DAVIS dataset. Compared with the visual odometry algorithm ORB and EVO, the proposed algorithm showed more advantageous performance in terms of accuracy and robustness.