PRIORITY CLAIM

This application claims the benefit of U.S. Provisional Patent Application Serial No. 62/459,039 filed February 14, 2017, the disclosure of which is incorporated herein by reference in its entirety.

TECHNICAL FIELD

This specification relates generally to computer systems for event- based feature tracking using asynchronous event-based sensors.

BACKGROUND

Asynchronous event-based sensors present new challenges in basic robot vision problems like feature tracking. The few existing approaches rely on grouping events into models and computing optical flow after assigning future events to those models. Such a hard commitment in data association attenuates the optical flow quality and causes shorter flow tracks.

SUMMARY

This specification describes a soft data association modeled with probabilities. The association probabilities are computed in an intertwined expectation maximization (EM) scheme with the optical flow computation that maximizes the expectation (marginalization) over all associations. In addition, to enable longer tracks we compute the affine deformation with respect to the initial point and use the resulting residual as a measure of persistence. The computed optical flow enables a varying temporal integration different for every feature and sized inversely proportional to the length of the flow. The results can be seen in egomotion and very fast vehicle sequences and we show the superiority over standard frame-based cameras.

The subject matter described herein may be implemented in hardware, software, firmware, or any combination thereof. As such, the terms "function" or "node" as used herein refer to hardware, which may also include software and/or firmware components, for implementing the feature(s) being described. In some exemplary implementations, the subject matter described herein may be implemented using a computer readable medium having stored thereon computer executable instructions that when executed by the processor of a computer control the computer to perform steps. Exemplary computer readable media suitable for implementing the subject matter described herein include non-transitory computer readable media, such as disk memory devices, chip memory devices, programmable logic devices, and application specific integrated circuits. In addition, a computer readable medium that implements the subject matter described herein may be located on a single device or computing platform or may be distributed across multiple devices or computing platforms.

BRIEF DESCRIPTION OF THE DRAWINGS

Figure 1 illustrates selected features tracked on a semi truck;

Figure 2 is a graphical outline of the algorithm;

Figure 3 shows the norm of feature position error between the method and Kanade-Lucas-Tomasi (KLT);

Figure 4 shows a comparison between frame-based and integrated- event images;

Figure 5 shows images of a truck driving on a highway;

Figure 6 includes a chart illustrating optical flow estimates and KLT tracking;

Figure 7 is a flow diagram of an example method for feature tracking in asynchronous event-based sensors that relies on probabilistic data association; and

Figure 8 shows an overview of the EVIO algorithm. DETAILED DESCRIPTION

I. Introduction

Robot vision continues to primarily rely on the main paradigm of frame-based cameras that acquire whole frames with fixed time exposure and frame rate. An alternative camera paradigm has emerged recently, that captures at almost unlimited frame rate changes in intensity and records events at specific time-points and image locations. Such sensors like the DVS [1 ] or DAVIS [2] cameras enable tasks entailing very fast motions and high dynamic range. However, to facilitate these tasks we have to redefine basic vision problems like optical flow or feature tracking because of the lack of representations that could make use of fundamental assumptions like the brightness change constraint [3].

The fundamental challenge underlying event-based tracking is the lack of any data association between event and established features. Grouping the events and searching for the most similar event group is impossible because events are arriving asynchronously and grouping would require a time window specification. Without knowing the optical flow we are unable to assign a new event to the closest feature unless the flow is under one pixel. In this document, we introduce a novel approach for feature definition by addressing the data association for each event to a feature as a hidden soft random variable. We regard the associations as probabilities because we do not need to make a hard commitment of an event to a feature. We apply an expectation-maximization (EM) scheme, where given optical flow we compute probabilities (weights) for data association and then we take the expectation over these probabilities in order to compute the optical flow. Inspired by the Good Features to Track approach [4] we model the alignment between events and features as an affine transformation and end the duration of the feature based on the quality of the alignment as well as the convergence of the EM iteration. Grouping of the events into a feature is not by a fixed spatiotemporal window but rather by a lifetime defined by a fixed length of the optical flow computed in the previous steps. Figure 1 illustrates selected features tracked on a semi truck driving at 60 miles/hr, 3 meters from a camera. Intermediate images are generated by integrating events for a period equal to three times their lifetimes.

We show in very fast motion scenarios that we can track robustly features over long trajectories. In this document, we make the following novel contributions to the state of the art of event-based tracking:

• Events are grouped into features based on lifespan defined by the length of the optical flow.

• Assignment of events to existing features is soft and computed as probability based on a predicted flow.

• Flow is computed as a maximization of the expectation over all data associations.

• Deformation of the feature is modeled as affine and the residual of the affine fit serves as a termination criterion.

II. Related Work

Litzenberger et al. [5], inspired by mean-shift tracking [6], create clusters of events by assigning events to the closest centroid. Each cluster is weighted by the mean frequency of the events and inactive clusters are eliminated. Kim et al. [7] estimate a 3D-rotation for the purpose of mosaicking by updating a particle filter with the likelihood of each new event given the current pose estimate. [8] propose an approach where an event is assigned to the spatially closest feature model and its Euclidean transformation and scaling with respect to the model is computed. Initialization is achieved by fitting spatiotemporal planes to the event space. Lagorce et al. [9] define features using the Hough transform and then assign events using the ICP principle. Tschechne et al. [10] introduced the notion of a motion streak using biological vision principles where event tracks are detected by tuning spatiotemporal orientation over a longer temporal support. [1 1 ] and [12] were the first to combine a frame based camera and event-based sensor on the same pixel array for tracking. Using a corner and an edge detector this approach initializes a feature patch which is enhanced by new events that are registered using a 2D Euclidean transformation. The commitment of an event to a feature is hard and hence the registration is prone to false event associations.

A common characteristic of the above approaches is the hard commitment, usually via ICP, to the assignment of an event to a model/feature with a subsequent estimate of a transformation given that commitment. Our approach integrates both data association and transformation estimation into a sound probabilistic scheme that makes it less prone to wrong correspondences. It does not make use of grayscale feature initializations. It is tested in very fast sequences where we show superiority over standard frame-based techniques. Barranco et al. [13] have created an evaluation dataset which offers ground truth optical flow but not longer feature trajectories. It provides self-motion estimates and we plan to use the dataset on the application of our feature tracking in visual odometry. III. Problem Formulation

Let be the projection of F onto the image plane at time t

where K is the camera calibration matrix and [R(t) 7(f)] is the camera pose. In the remainder, we refer to the projections f(t) as features and consider a set of features F(t). Given a feature define a spatial window be a set of 3-D

points, whose projections onto the image plane at time 0 are

contained within the window B(s). Let P ^f (t) denote the set of point projections associated with feature at time t. At discrete times fi , ... , tn, the sensor generates a set of events where

and j _{S an un} k _n own many-to-one function

representing the data association between the events {e,} and projections {ø,}. Problem (Event-based Feature Tracking). Given a set of events ε generated by the point projections estimate the feature

projections F(t) in the image plane over time.

IV. Method

In Section IV-A, we introduce an optical flow based constraint that groups events projected from the same point within a spatiotemporal window. Section IV-B then shows that we can optimize this constraint over the optical flow using the Expectation Maximization algorithm. The resulting flow can then be used to reconstruct the set of point projections within the spatial window, which we then use in Section IV-C to refine the feature position using the expectation maximization-iterative closest point (EM-ICP) algorithm [14]. Our tracking method then iterates between Section IV-B and Section IV-C to track a given set of features over time in the event stream. Section IV-D outlines our technique to select the size of the temporal windows in an asynchronous fashion, Section IVE details our method for initializing the feature positions, and Section IV-F summarizes our method for estimating the image features within each window.

A Spatiotemporal Optical Flow Constraint

The motion of a feature in the image plane can be described using its optical flow f(t) as follows:

Where v(t) is the average flow of f(0) over time. If t is sufficiently small, we can assume that the average flow v is constant and equal to the average flows of all point projections P(0) associated with f(0). We can define a spatiotemporal window around f(0) as the collection of events up to time t that propagate backwards onto S(s): Thus, provided that t is small, events corresponding to the same point in P(0) should propagate backwards onto the same image location. In other words, the following equality should hold for any pair of events:

However, since the data association π between events and 3D points is unknown, we can hope to satisfy the above requirement only in expectation:

where ^and we assume that ττ(/) is independent of TT(K ).

Given an affine transformation (A,b) and the flow v of feature f(0), we model the noise in the event generation process by defining the probability that event e,- was generated from point p _s as proportional to the pdf of a Normal distribution with mean p _j and covariance ∑, i.e.,

Where the argument {pj} is the set of points over which the means are defined. From here on, we will assume that rij with no argument implies that the set is the point projections {pj}. Note also that∑ is a parameter.

One can see that, by minimizing the distance between propagated events that correspond to the same point, we are effectively enforcing a sparsity constraint over the spatial domain into which the events propagate.

We propose an iterative approach to estimate the data association probabilities rij between the events and points the affine transformation A, b, and the optical flow v of feature f.

B. EM Optical Flow Estimation

In this section, we propose an Expectation Maximization algorithm for solving (5) over a spatiotemporal window W(s, t) with a set of events

Within this window, our optical flow constraint becomes

In the E step, we update the r _tj and rkj, given v using (6). Initially, the set of point positions {pj} is unknown, and so we first approximate the {pj} by the set of propagated events

where v' is the true optical flow. In addition, as A and b are unknown, we initialize them as A = I and b = 0. The full update, then, is rij({e _/ }).

The M step now involves solving for v given the As we assumed that the average optical flow v is constant, (7) is a linear least squares problem in v, which corresponds to the general overdetermined system:

To get the normal equations, we multiply both sides on the right

We iterate equations (6) and (8) until convergence of the error (4). As in [14], we reject outlier matches by thresholding the likelihood w _ik when computing (8) by setting all the w _ik higher than some threshold ε to 0.

C. Feature Alignment

The flow estimate from Section IV-B can then be used to propagate the events within the window to a common time to- Give a correct flow, this set of events is then the approximation to the projection of the points P _j at time t ₀ , up to an affine transformation. As a result, given the set of point projections at time to, _ _{we can} align the events with their corresponding points using a similar EM formulation as in Section IV-B. The cost function for this alignment is that of the EM-ICP algorithm [14]: We can minimize this cost function using exactly the steps from Section IV-B. In the E step, we can use (6) to update

The M step is also similar:

Where:

As in Section IV-B, we iterate (6) and (1 1 ) until the error function (10) converges. We then use the new estimate for b to refine the prior estimate for the image feature positions, and propagate them to the next window as:

If the EM fails to converge, or if the value of (10) is too high after convergence, we consider the tracker lost and abandon the feature.

D. Temporal Window Selection

Many event based techniques have worked on temporal windows of fixed size. However, these techniques have many similar drawbacks to traditional cameras, as the amount of information within the windows is highly variable depending on the optical flow within the image. Due to the quantization of the spatial domain, no information about the optical flow can be gained from the event stream until the projected points have moved at least one pixel within the image. On the other hand, too large of a window may violate the constant optical flow assumption made in Section IV-B. To ensure that the temporal windows are of an appropriate size, we dynamically adjust them using the concept of event 'lifetimes' [15]. Given the optical flow of a feature within a prior spatiotemporal window, we can estimate its lifetime T as the expected time for the feature to move one pixel in the spatial domain: For robustness against erroneous flow estimates, we estimate the lifetimes of several windows, and set the next temporal window size as k times the median lifetime. In our experiments, we observed that k = 3 was a reasonable value to avoid large optical flow deviations while still capturing a sufficient number of events. This technique can be extended to have separate temporal window sizes for each tracked feature for fully asynchronous tracking between features. However, we found that, in our testing, the variation in optical flow of our features was not large enough to require this.

E. Feature Selection

As this method relies on the assumption that the projected points are sparse, it will fail on spatial windows with dense points throughout. In addition, the matching scheme in Section IV-C suffers from the same aperture problem as traditional feature matching techniques. To avoid selecting such windows for our tracker, we propagate all events within a temporal window onto the image plane with zero flow to generate an image. We then use the Harris corner detector [16] to select spatial windows with large variations, i.e., variations that exceed a configurable threshold number of variations, in image intensity (in this case, event density) in multiple directions.

F. Point Set Generation

On the first iteration of this algorithm, we run the flow estimation in Section IV-B, and propagate events within the first spatiotemporal window to the time of the last event in the window, T. In the next step, we use these propagated events as the projected point positions, and match the events in the next window propagated to time T. We then concatenate the next set of propagated, aligned events to our point set, up to some fixed maximum number of points N, and continue to propagate this set to the beginning of the next temporal window at each iteration. To reduce the computation time for matching against this potentially large feature set, we perform the sphere decimation algorithm in [14] to reduce the cardinality of this set once N points are collected.

Figure 2 is a graphical outline of the algorithm. A first image 200 shows an event stream within a spatiotemporal window. Diagonal lines are formed by the linear optical flow. A second image 202 shows events integrated directly onto the image with no flow correction. A third image 204 shows propagated events with estimated flow. Motion blur has been removed. A fourth image 206 shows a later set of propagated events before affine warping. The size of the lighter circles are the weights of each point after decimation. A fifth image 208 shows a later set of propagated events after affine warping.

V. Experiments

We present the results of our approach using a DAVIS- 240C sensor [2] in two situations. First, we compare the tracking accuracy of our tracking algorithm on a structured, textured area at normal speeds against traditional image based tracking on the frame-based intensity values from the DAVIS. We then demonstrate qualitative results of our algorithm on tracking a vehicle on a highway travelling at roughly 60 miles/hr, which we qualitatively compare to the tracking results on the 240 frames per second (FPS) output of an iPhone 6.

In each experiment, we used 31x31 pixel patches, with∑, set to 2 x / ₂ . At the beginning of each sequence, a manually picked integration time is selected to start the algorithm to guarantee that the first temporal window contained significant apparent motion in the spatial domain. However, this time is very robust, and we found that any integration time where points moved at least 2 pixels was sufficient. In both experiments, 20 features were generated, with new features initialized if fewer than 12 remained. A. Comparison with Frame-Based Tracking

To quantitatively analyze the performance of our algorithm, we recorded a sequence where the DAVIS camera was moved in front of a textured surface (Figure 3), and compare our results to the KLT tracker [17]. Due to the relatively low frame rate of 25Hz for the frame based images on the DAVIS, this motion was restricted to relatively low speeds. Features were initialized from the integrated event image, and tracked in both the event stream as well as the frame based images until the majority of features were lost in both trackers. During the one second tracking period, the features moved on average 100 pixels.

Figure 3 shows the norm of feature position error between the method and KLT. Figure 4 shows a comparison between frame-based and integrated-event images. Figure 5 shows images of a truck driving on a highway recorded from the 240 FPS video.

We show the mean tracking error for features that had not been discarded in Figure 3, where the black line is the mean tracking error over all the features, and the lighter region is one standard deviation around the mean error. As the event based measurements arrive much faster than the frame based ones, we interpolate the event based feature position estimates to the nearest frame based position estimate in time using the event based optical flow estimate. The overall mean error from this technique is 0.9492 pixels, which is comparable to the state of the art in this topic. B. Tracking on Scenes with High Apparent Motion

To test the algorithm on scenes with very high apparent motion, the camera was placed on the side of a highway with a speed limit of 60 miles per hour. Cars passed the camera at a distance between 3-4 meters, and passed the field of view in under 0.5s. We present here the results of tracking on a semi truck driving by at the posted speed limit. In this sequence, the average flow magnitude was 4000 pixels/s, and the (roughly) 15m long truck passed the camera's field of view in 800ms. The frame based images from the DAVIS sensor for these vehicles were almost completely blurred out. For comparison, we also recorded the scene with an iPhone 6 at 240 FPS (Fig. 5), on which we also ran a KLT tracker. The 240 FPS video is sufficient to capture the motion in this sequence, but is beginning to show motion blur on the order of one or two pixels. The two cameras' extrinsic parameters were estimated using stereo calibration. Unfortunately, due to the relatively close distance of the vehicles to the camera, we were unable to accurately warp the images onto one another for a quantitative comparison, and so we will instead give qualitative comparisons for our flow estimation based on a warp of the iPhone images assuming that points all have a depth of 3 meters.

We visualize a subset of the feature tracks in Fig. 1 . It is interesting to note that, while the first integrated event image (superimposed over the iPhone image) has a significant amount motion blur, the subsequent images have structures only a few pixels thick, due to the lifetime estimation from the optical flow.

In Fig. 6, we analyze the distribution of the direction of the optical flow vectors estimated by our method and by the KLT tracker. Figure 6 includes a chart 600 illustrating optical flow estimates (light areas) and KLT tracking (dark areas). Figure 6 also includes a polar histogram 602 (20 bins) of optical flow directions estimated by our method. Figure 6 also includes a polar histogram 604 (20 bins) of optical flow directions estimated by KLT.

We can see that the majority of flow vectors lie between 0 and 20 degrees. This can also be seen in the chart 600 in Fig. 6, which shows individual flow vectors, with optical flow calculated within tracks shorter than 20ms removed. From these plots, we can see that both the direction and magnitude of the KLT flow vectors are very similar, although they should not perfectly correspond.

VI. Conclusion

This specification describes a novel approach for feature tracking in asynchronous event-based sensors that relies on probabilistic data association. Estimating optical flow becomes, thus, not sensitive to erroneous associations of new events and is computed from the expectation over all associations. To increase persistence of our tracks we compute the affine transformation for each feature with respect to the starting time. Existing approaches use a hard correspondence commitment and usually compute a similitude transformation. The spatiotemporal support of our features is adaptive and defined by the size of the flow rather than a fixed time or a number of events. We show that it outperforms classic KLT trackers when they are applied to 240 FPS cameras capturing very fast motions of the order of one field of view per half a second. We plan a realtime implementation so that we can apply it in real robot scenarios of visual odometry and moving object tracking.

Figure 7 is a flow diagram of an example method 700 for feature tracking in asynchronous event-based sensors that relies on probabilistic data association. The method 700 can be performed, for example, by a computer system programmed to perform aspects of the algorithm described above with reference to Figures 1 -6. The method 700 includes receiving asynchronous events from asynchronous event-based sensor (702). The method 700 includes initializing one or more features by integrating events from an initial temporal window into an initial image and detecting the features in the initial image (704). The method 700 includes tracking the features after the initial temporal window, including assigning events corresponding to later temporal windows to features by determining, for each event and each feature, a probability that the event is associated with the feature (706). The method 700 can end when the tracking indicates that the feature is lost. Based on the tracking, motion of one or more physical objects can be identified.

The following sections present an algorithm to fuse a purely event- based tracking algorithm with an inertial measurement unit, to provide accurate metric tracking of a camera's full 6dof pose. The algorithm is asynchronous, and provides measurement updates at a rate proportional to the camera velocity. The algorithm selects features in the image plane, and tracks spatiotemporal windows around these features within the event stream. An Extended Kalman Filter with a structureless measurement model then fuses the feature tracks with the output of the IMU. The camera poses from the filter are then used to initialize the next step of the tracker and reject failed tracks.

Problem Formulation

Consider a sensor package consisting of an inertial measurement unit (IMU) and an event-based camera. Without loss of generality, assume that the camera and IMU frames are coincident. The state of the sensor package:

consists of its position p ε R ³ , its velocity v ε R ³ , the orientation of the global frame in the sensor frame represented by a unit quaternion q ε SO(3), and the accelerometer and gyroscope measurement biases, ba and bg, respectively.

At discrete times the IMU provides acceleration and

angular velocity measurements The environment, in

which the sensor operates, is modeled as a collection of landmarks L :=

At discrete times U , t ₂ , . . ., the event-based camera generates events E := {(Xj, tj)} which measure the perspective projection of the landmark positions as follows:

where _a n unknown function representing

the data association between the events E and the landmarks L and R(q) is the rotation matrix corresponding to the rotation generated by quaternion q.

Problem 1 (Event-based Visual Inertial Odometry). Given inertial measurements I and event measurements E, estimate the sensor state s(t) over time.

Overview

The visual tracker uses the sensor state and event information to track the projections of sets of landmarks, collectively called features, within the image plane over time, while the filtering algorithm fuses these visual feature tracks with the IMU data to update the sensor state. In order to fully utilize the asynchronous nature of event-based cameras, the temporal window at each step adaptively changes according to the optical flow. The full outline of our algorithm can be found in Fig. 8 and Alg. 2.

Algorithm 2 EVIO

Figure 8 shows an overview of the EVIO algorithm. Data from the event-based camera and IMU is processed in temporal windows determined by our algorithm. For each temporal window, the IMU values are used to propagate the state, and features are tracked using two expectation maximization steps that estimate the optical flow of the features and their alignment with respect to a template. Outlier correspondences are removed, and the results are stored in a feature track server. As features are lost, their feature tracks are parsed through a second random sample consensus (RANSAC) step, and the resulting tracks are used to update the sensor state. The estimated optical flows for all of the features are then used to determine the size of the next temporal window.

The algorithm leverages the property that all events generated by the same landmark lie on a curve in the spatiotemporal domain and, given the parameters of the curve, can be propagated along the curve in time to a single point in space (barring error from noise). In addition, the gradient along this curve at any point in time represents the optical flow of the landmark projection, at that time. Therefore, we can reduce the tracking problem to one of finding the data association between events and landmark projections, and estimating the gradient along events generated by the same point.

To simplify the problem, we make the assumption that optical flow within small spatiotemporal windows is constant. That is, the curves within these windows become lines, and estimating the flow is equivalent to estimating the slope of these lines. To impose this constraint, we dynamically update the size of the temporal window, dt, based on the time for each projected landmark position to travel k pixels in the image, estimated using the computed optical flow. This way, we are assuming constant optical flow for small displacements, which has been shown to hold in practice. Given these temporal windows, we compute feature position information within a discrete set of nonoverlapping time windows {[T-i , T ₂ ], [T ₂ , T ₃ ], . . . } where T _i+1 - 7} = dt,.

Problem 1a (Event-based Feature Tracking). Given event measurements E, the camera state s,- := s(T,) at time and a temporal window size dt _h estimate the feature projections F(t) for t ε [Ti·, T _/ +1 ] in the image plane and the next temporal window size dt _i+1 .

Given the solution to Problem 1 a and a set of readings from the IMU, our state estimation algorithm then employs an Extended Kalman Filter with a structureless vision model. Structured extended Kalman filter (EKF) schemes impose vision constraints on all the features between each two subsequent camera poses, and as a result optimize over both the camera poses and the feature positions. A structureless model, on the other hand, allows us to impose constraints between all camera poses that observe each feature, resulting in a state vector containing only the IMU and camera poses. This drastically reduces the size of the state vector, and allows for marginalization over long feature tracks. We pose the state estimation problem as the sub problem below:

Problem 1 b (Visual Inertial Odometry). Given inertial measurements / at times {t _k }, a camera state s,- at time 7, a temporal window size dt _h and feature tracks F(t) for t ε [Ti·, T _i+1 ], estimate the camera state S _/ +i.

Event-based Feature Tracking

Given a camera state Sj and a temporal window [Ti·, T _i+1 ], the goal of Problem 1 a is to track a collection of features F(t) in the image domain for t ε [Ti, T _i+1 ], whose positions are initialized using traditional image techniques. Our approach performs two expectation-maximization (EM) optimizations to align a 3-D spatiotemporal window of events with a prior 2-D template in the image plane. Optical Flow Estimation

The motion of a feature f(t) in the image plane can be described using its optical flow f(t) as follows: where the average flow of f(s) for j _s sufficiently small, we can assume that the average

flow u is constant and equal to the average flows of all landmarks L _j whose projections are close to f in the image plane. To formalize this observation, we define a spatiotemporal window around f(Ti) as the collection of events, that propagate backwards to and landmarks whose projections at time T,- are close to the vicinity of f(Ti) : where is the window size in pixels,

as defined above.

Provided that [Ti·, T _i+1 ] is small, the following equality should hold for any event e _k and landmark L _j in Wf.

where the indicator requires that event k associates with landmark j. However, since the data association between events and landmarks in unknown, we can hope to satisfy the above requirement only in expectation: where Given the flow u of feature f(T _/ ), due

to the measurement model in (2), we model the probability that event e _k was generated from landmark L _j as proportional to the probability density function of a Gaussian distribution with mean /,· and covariance

be an enumeration of the events in Wj. The optical flow constraints in (5) and the data association probabilities r _kJ allow us to estimate the optical flow u of feature f(Tj) as follows:

Unfortunately, the landmark projections

compute the data association probabilities are unknown approximate in the above optimization with the set

of forward-propagated events from the previous iteration to the current time. This set provides a close approximation to because,

as Uj-i approaches the true optical flow, (4) requires that every point in approaches some point in l _j at time Τ,·, due to continuity in the

projected feature trajectories. This leads to the following result for estimating the optical flow of feature f(Ti).

Proposition 1

Given a sufficiently small temporal window so that the average optical flow u of a feature IS

approximately constant and a spatiotemporal window W _i of n _i events associated with f, the optical flow u can be estimated by iterating the following EM steps:

Proof. Given an optical flow estimate, the E step computes the likelihood /¾· of the data association between events e _k and approximate landmark projections by propagating the events backwards in time and applying the measurement model in (2). Given r _kJ , the M step is a weighted linear least squares problem in u, which corresponds to the overdetermined system ud ^T = Y, where:

To get the normal equations, we multiply both sides on the right by d and obtain u = (Y d)/(d ^T d).

During the initialization of each feature, no prior flow is known, and so we further substitute the prior events in with the current events and flow, O ~nce aga .in, t..h.is approx .imat ..ion

approaches the true projected landmark positions as u approaches the true flow. The M step, in this case, becomes:

where r _kJ is computed as in (8). This method is significantly slower, as the distance computation in the E step uses two different sets of points every iteration. This is detrimental for most neighbor search data structures such as k-d trees, as the data structure must be reconstructed at every iteration. However, as this step is only performed once per feature, it is effectively marginalized out for long feature tracks.

Template Alignment and RANSAC

While Prop. 1 is sufficient to solve Problem 1 a, the feature position estimates may drift as a feature is being tracked over time, due to noise in each flow estimate. To correct this drift, we estimate it as the affine warping that warps (7) to the template in the first camera pose that observed the feature. We assume that the corresponding landmarks {/} are planar in 3-D, so that we can alternatively represent the affine wrapping as a 3-D rotation and scaling. The 3-D rotation '*Rj from the current camera pose at Tj to the first camera pose at T, _* can be obtained from the filter used to solve Problem 1 b. Hence, in this section we focus on estimating only a scaling ^σ and translation b between First, we rotate each

point to the first camera frame and center at the rotated feature position as follows: where is the projection function defined in (2). Note that ^,¾r propagates events to time T _i+1 , and so we substitute f(7}) + Ujdtj as an estimate for f(T, ₊₁ ). Then, using the same idea as above, we seek the scaling and translation b that minimize the mismatch between and

This optimization problem has a similar form to problem (6) and, as before, can be solved via the following EM steps: where the M step is solved as in scaled ICP. Algorithm 3 Event-based Feature Tracking

Outlier Rejection

In order to remove outliers from the above optimizations, only pairs of points and approximate projected landmarks, with Mahalanobis distance below a set threshold are used in

the optimizations. This also serves the purpose of heavily reducing the number of computation.

After all of the features have been updated, two-point RANSAC is performed given the feature correspondences and the rotation between the frames from the state to remove features whose tracking have failed. Given two correspondences and the rotation, we estimate the essential matrix, and evaluate the Sampson error on the set of correspondences to determine the inlier set. The complete feature tracking procedure is illustrated in Fig. 2 and summarized in Alg. 3.

Temporal Window Size Selection

To set the temporal window size such that each feature moves k pixels within the window, we leverage the concept of 'event lifetimes'. Given an estimate of the optical flow, u, of a feature f at a given time, T), the expected time for the feature to move 1 pixel in the image is Therefore, to estimate the time for a feature to move k pixels, the time is simply dt{f) = Given a set of tracked features, we set the next

, . , .

temporal window size as:

Assuming that the differences in depth between the features is not large, this window size should ensure that each feature will travel roughly k pixels in the next temporal window. For all of our experiments, k is set to 3.

Feature Detection

Like traditional image based tracking techniques, our event-based feature tracker suffers from the aperture problem. Given a feature window with only a straight line, we can only estimate the component of the optical flow that is normal to the slope of this line. As a result, features must be carefully selected to avoid selecting windows with a single strong edge direction. To do so, we find 'corners' in the image that have strong edges in multiple directions. In order to produce an image from the event stream, we simply take the orthogonal projection of the events onto the image plane. As we constrain each temporal window such that features only travel k pixels, this projection should reconstruct the edge map of the underlying image, with up to k pixels of motion blur, which should not be enough to corrupt the corner detection. The actual corner detection is performed with FAST corners, with the image split into cells of fixed size, and the corner with the highest Shi-Tomasi score within each cell being selected. State Estimation

To estimate the 3D pose of the camera over time, we employ an Extended Kalman Filter with a structureless vision model. For compactness, we do not expand on the fine details of the filter. At time 7}, the filter tracks the current sensor state (1 ) as well as all past camera poses that observed a feature that is currently being tracked. The full state, then, is:

where n is the length of the oldest tracked feature.

Between update steps, the prediction for the sensor state is propagated using the IMU measurements that fall in between the update steps. Note that, due to the high temporal resolution of the event based camera, there may be multiple update steps in between each IMU measurement. In that case, we use the last IMU measurement to perform the propagation.

Given linear acceleration a _k and angular velocity measurements, the sensor state is propagated using 5th order Runge-Kutta integration:

To perform the covariance propagation, we adopt a discrete time model and covariance prediction update. When an update from the tracker arrives, we augment the state with a new camera pose at the current time, and update the covariance using the Jacobian that maps the IMU state to the camera state.

We then process any discarded features that need to be marginalized. For any such feature ή , the 3D position of the feature can be estimated using its past observations and camera poses by Gauss Newton optimization, assuming the camera poses are known. The projection of this estimate into a given camera pose can then be computed using (2). The residual, r®, for each feature at each camera pose is the difference between the observed and estimated feature positions. We then left multiply r® by the left null space, A, of the feature Jacobian, H _F , to eliminate the feature position up to a first order approximation:

The elimination procedure is performed for all features, and the remaining uncorrelated residuals, are stacked to obtain the final residual r ₀ . We perform one final step to reduce the dimensionality of the above residual. Taking the QR decomposition of the matrix Ho, we can eliminate a large part of the residual. The EKF update step is then where K is the Kalman gain.

When a feature track is to be marginalized, we apply a second RANSAC step to find the largest set of inliers that project to the same point in space, based on reprojection error. This removes moving objects and other erroneous measurements from the track.

Algorithm 4 State Estimation

Although specific examples and features have been described above, these examples and features are not intended to limit the scope of the present disclosure, even where only a single example is described with respect to a particular feature. Examples of features provided in the disclosure are intended to be illustrative rather than restrictive unless stated otherwise. The above description is intended to cover such alternatives, modifications, and equivalents as would be apparent to a person skilled in the art having the benefit of this disclosure.

The scope of the present disclosure includes any feature or combination of features disclosed in this specification (either explicitly or implicitly), or any generalization of features disclosed, whether or not such features or generalizations mitigate any or all of the problems described in this specification. Accordingly, new claims may be formulated during prosecution of this application (or an application claiming priority to this application) to any such combination of features. In particular, with reference to the appended claims, features from dependent claims may be combined with those of the independent claims and features from respective independent claims may be combined in any appropriate manner and not merely in the specific combinations enumerated in the appended claims. References

The following documents are hereby incorporated by reference:

[I ] P. Lichtsteiner, C. Posch, and T. Delbruck, "A 128_ 128 120 db 15 ps latency asynchronous temporal contrast vision sensor," IEEE journal of solid-state circuits, vol. 43, no. 2, pp. 566-576, 2008.

[2] C. Brandii, R. Berner, M. Yang, S.-C. Liu, and T. Delbruck, "A 240x 180 130 db 3 μs latency global shutter spatiotemporal vision sensor,"

IEEE Journal of Solid-State Circuits, vol. 49, no. 10, pp. 2333-2341 ,

2014.

[3] B. K. Horn and B. G. Schunck, "Determining optical flow," Artificial intelligence, vol. 17, no. 1 -3, pp. 185-203, 1981 .

[4] J. Shi and C. Tomasi, "Good features to track," in Computer Vision and Pattern Recognition, 1994. Proceedings CVPR'94., 1994 IEEE

Computer Society Conference on, pp. 593-600, IEEE, 1994.

[5] M. Litzenberger, C. Posch, D. Bauer, A. Belbachir, P. Schon, B. Kohn, and H. Garn, "Embedded vision system for real-time object tracking using an asynchronous transient vision sensor," in 2006 IEEE 12th

Digital Signal Processing Workshop & 4th IEEE Signal Processing

Education Workshop, pp. 173-178, IEEE, 2006.

[6] D. Comaniciu, V. Ramesh, and P. Meer, "Real-time tracking of nonrigid objects using mean shift," in Computer Vision and Pattern

Recognition, 2000. Proceedings. IEEE Conference on, vol. 2, pp. 142- 149, IEEE, 2000.

[7] H. Kim, A. Handa, R. Benosman, S. H. Ieng, and A. J. Davison,

"Simultaneous mosaicing and tracking with an event camera," J. Solid State Circ, vol. 43, pp. 566-576, 2008.

[8] Z. Ni, S. H. Ieng, C. Posch, S. R ^' egnier, and R. Benosman, "Visual tracking using neuromorphic asynchronous event-based cameras,"

Neural computation, 2015.

[9] X. Lagorce, C. Meyer, S. H. Ieng, D. Filliat, and R. Benosman,

"Asynchronous event-based multikernel algorithm for high-speed visual features tracking," IEEE transactions on neural networks and

learning systems, vol. 26, no. 8, pp. 1710-1720, 2015.

[10] S. Tschechne, T. Brosch, R. Sailer, N. von Egloffstein, L. I. Abdul- Kreem, and H. Neumann, "On event-based motion detection and

integration," in Proceedings of the 8th International Conference on

Bioinspired Information and Communications Technologies, BICT '14,

2014.

[I I ] D. Tedaldi, G. Gallego, E. Mueggler, and D. Scaramuzza, "Feature detection and tracking with the dynamic and active-pixel vision sensor

(davis)," in Int. Conf. on Event-Based Control, Comm. and Signal

Proc.(EBCCSP), Krakow, Poland, 2016.

[12] B. Kueng, E. Mueggler, G. Gallego, and D. Scaramuzza, "Low-latency visual odometry using event-based feature tracks," in IEEE/RSJ International Conference on Robotics and Intelligent Systems, IEEE/RSJ,

2016.

[13] F. Barranco, C. Fermuller, Y. Aloimonos, and T. Delbruck, "A dataset for visual navigation with neuromorphic methods," Frontiers in neuroscience, vol. 10, 2016.

[14] S. Granger and X. Pennec, "Multi-scale em-icp: A fast and robust approach for surface registration," in European Conference on Computer Vision, pp. 418-432, Springer, 2002.

[15] E. Mueggler, C. Forster, N. Baumli, G. Gallego, and D. Scaramuzza, "Lifetime estimation of events from dynamic vision sensors," in 2015

IEEE International Conference on Robotics and Automation (ICRA), pp. 4874-4881 , IEEE, 2015.

[16] C. Harris and M. Stephens, "A combined corner and edge detector.," Citeseer, 1988.

[17] B. Lucas and T. Kanade, "An iterative image registration technique with an application to stereo vision," in Int. Joint Conf. on Artificial

Intelligence (IJCAI), vol. 81 , pp. 674-679, 1981 .