Asynchronous Blob Tracker for Event Cameras

Event cameras report the relative log intensity change of brightness for each pixel asynchronously. We consider the likelihood $\ell(\xi,t,\rho|e_{k})$ of an event $e_{k}$ occurring at pixel location $\xi\in\mathbb{R}^{2}$ at time $t\geq 0$ and with polarity $\rho\in\{\pm 1\}$. We will term such a likelihood model an event blob if the spatial distribution of the conditional likelihood $\ell(\xi,\rho|t,e_{k})$ for a fixed time $t$ is blob-like. We propose a spatio-temporal Gaussian event blob likelihood function

where $p(t)\in\mathbb{R}^{2}$ is the pixel location of the centre of the object and $\Lambda(t)\in\mathbb{R}^{2\times 2}$ is the shape of the object (encoded as the principal square root of the second-order moment of the event spatio-temporal intensity). The positive-definite matrix $\Lambda(t)>0$ is written

for principal correlations $\lambda^{1}(t),\lambda^{2}(t)>0$ and orientation angle $\theta(t)$ where $R(\theta(t))$ is the associated rotation matrix. The scalar $\gamma(t)$ encodes the temporal dependence of the likelihood associated with changing event rate and $B(\rho|\xi,t)$ denotes a binomial distribution associated with whether the event has polarity $\{\pm 1\}$ depending on position and time. If $\gamma$ is integrable on a time interval of interest, then the spatio-temporal Gaussian event blob likelihood could be normalised to produce a probability density on this time interval. In this paper, however, we will use the conditional likelihood $\ell(\xi,\rho|t,e_{k})$, a spatial Gaussian distribution with time-varying parameters. In the sequel, we will often omit the time index from the state variables $p=p(t)$, $\theta=\theta(t)$, etc., to make the notation more concise.

The archetypal examples of event blob targets are flickering objects such as an LED or fluorescent light [44]. A flickering target generates events at pixel locations in proportion to the frequency of the flicker and intensity of the source at that pixel. The likelihood of an event occurring at a given moment in time is proportional to the rate of events and the binomial probability of polarity is independent of position. Thus, for a target with constant flicker, the spatial distribution of the unipolar (ignoring event polarity) event blob likelihood is directly related to the intensity of the source and is analogous to image intensity. An example of the event arrival histogram for LED tail lights of cars is shown in the top row of Figure 2. Here the image blobs in Figure 2b-c are roughly Gaussian.

In the case of a non-flickering target, the situation is more complex since the motion of the target in the image plane is required to generate events. In this case, the events are asymmetrically arranged around the target centre depending on whether the target is bright or dark with respect to the background and the binomial distribution $B(\rho|\xi,t)$ is not independent of the spatial parameter $\xi$. For a bright target against a dark background, the events in the direction of motion are positive while the trailing events are negative, and vice-versa if the background is brighter than the target. If both positive and negative events are equally considered, then the resulting density is symmetric around the centre point of the target (for an idealised sensor). The rate of events at the target centre should be zero since the intensity of the target blob at this point is a maximum (or minimum) with respect to the direction of motion of the target. The second row of Figure 2 provides a good example of a non-flickering moving target. The density of events in Figure 2e clearly shows the leading and trailing edge effects where the target is moving to the left and down in the image. The effective likelihood over the whole target, however, appears to be a single blob as seen in Figure 2f.

In practice, we have found that a simple event blob model and the associated spatio-temporal Gaussian likelihood model proposed in (1) works well for a wide range of blob-like intensity targets. For non-flickering targets (e.g., Figure 4), we introduce a minor modification to the proposed algorithm (cf. §V-D) to compensate for the bimodal distribution discussed above.

Consider a sequence of events ${e_{k}:=(\xi_{k},t_{k},\rho_{k})}$ for k $\geq$ 0, at pixel locations $\{\xi_{k}\}$, times $\{t_{k}\}$ and with polarities $\{\rho_{k}\}$ associated with a target event blob modelled by (1). We will consider the unipolar model for event blobs where the polarity $\rho_{k}$ is ignored and all events are considered equally. The associated likelihood is the marginal

In addition, we will ignore noise in the time-stamp and consider the conditional likelihood $\ell(\xi|t=t_{k},e_{k})$ fixing time at $t=t_{k}$. With modern event cameras, where the time resolution is as low as 1$\mu$s resolution this assumption is justified in a wide range of scenarios. Based on these assumptions then a generative noise model for the event location is given by

where coordinate $p(t_{k})$ is the true location of the object and we emphasise that the times $\{t_{1},t_{2},\ldots,t_{k},\ldots\}$ are not periodic and depend on the asynchronous time-stamp of each event. The measurement noise $\eta_{k}$ is a zero-mean independent and identically distributed (i.i.d.) Gaussian process with known covariance $I_{2}$. The noise process is scaled through the square-root covariance $\Lambda(t_{k})$ given by (2) with parameters $\lambda^{1}(t_{k}),\lambda^{2}(t_{k})>0$ and $\theta(t_{k})$.

For a 2D object in the image plane, the proposed model considers the states $p=(p_{x},p_{y})^{\top}$ for the object centre, and the shape parameters $(\lambda^{1},\lambda^{2})^{\top}$ and $\theta$ at pixel-level. The parameter $\theta$ can be thought of as a spatial state parameter describing the orientation of the object, while the principal correlations $\lambda=(\lambda^{1},\lambda^{2})^{\top}$ encode the eigenvalues of the second-order ‘visual’ moment of the object. We add a linear velocity state $v=(v_{x},v_{y})^{\top}$ and angular velocity $q$ for the spatial states $(p_{x},p_{y},\theta)$. The state vector $x$ is

If an Inertial Measurement Unit (IMU) is available on the camera then the gyroscope measurement $\Omega=(\Omega_{x},\Omega_{y},\Omega_{z})^{\top}$, can be used to provide feed-forward prediction for the state evolution in $p$ and $\theta$. Define the matrix

and the partial state function [45, Chapter 15.2.1]

where $\bar{p}_{x}=(p_{x}-p_{x}^{0}$) and $\bar{p}_{y}=(p_{y}-p_{y}^{0}$) are the pixel coordinates relative to the principal point and $f$ is the camera focal length. The proposed state model is

where we model uncertainty in the state evolution using continuous Wiener processes $w_{p}$, $w_{v}$, $w_{\theta}$, $w_{q}$ and $w_{\lambda}$ with $Q^{p},Q^{v},Q^{\lambda}\in\mathbb{R}^{2\times 2}$ and $Q^{\theta},Q^{q}\in\mathbb{R}^{1\times 1}$ symmetric positive definite matrices.

Image blob velocity is modelled as a filter state, while change in velocity (acceleration) is compensated for in the process noise. Camera angular velocity data $\Omega$ is the major source of ego-motion induced optical flow for image tracking tasks. Ego-motion associated with camera IMU accelerometer measurements is not modelled in the present work as it is typically insignificant compared to the rotational ego-motion and independent motion of a target. In addition, since there is often no reliable linear velocity sensor available on a real-world robotic system, the feature velocity estimation provided by the AEB tracker may be the best available sensor modality to estimate linear motion of the camera. Due to its low latency and high temporal resolution, the feature velocity estimate of the AEB tracker is ideally suited to be used in this way. Even if the cameras linear velocity is available, it is necessary to know or estimate the range of the feature to include feature kinematics in the filter, introducing additional complexity. The benefit of such an approach is unclear since the AEB tracker can easily track all but the most extreme motion without including feature kinematics.

The predicted state $\hat{x}$ is computed by integrating the dynamics (8) without noise. We linearise the dynamics (8) about the predicted state $\hat{x}$ in order to predict the covariance $\hat{\Sigma}$.

Linearising $f_{p}$ in Equation (7) about $\hat{p}$ (the position components of $\hat{x}$) yields

It follows that the dynamics of the true state $x$, linearised about the estimated state $\hat{x}$, are

where $w$ is Wiener process with positive definite matrix

and the state matrix $A=A(t)$ is a block matrix corresponding to the state partition of (5),

with $J$ as defined in (6). Note that, if $\Omega$ is unavailable or zero (the camera is stationary), then $A_{\hat{p}}=0_{2\times 2}$. In this case, the dynamics (9) are exactly linear and are thus independent of the linearisation point $\hat{x}$; that is, $\mathrm{d}x=A(x-\hat{x})\mathrm{d}t+\mathrm{d}w$ and $A$ does not depend on $\hat{x}$.

The state covariance $\hat{\Sigma}$, in the absence of measurements, is predicted as the solution of a continuous-time linear Gaussian diffusion process given by

where $A$ is computed around $\hat{x}$.

Let $\{t_{k}\}$ be a sequence of times associated with an event stream. For each given time $t_{k}$, define the asynchronous system-state to be

Analogously, define asynchronous estimated state parameters $(\hat{x}_{k},\hat{\Sigma}_{k})$ at time $t_{k}$. Given estimated state parameters at time $t_{k}$ the predicted state $\hat{x}_{k+1}^{-}:=x(t_{k+1})$ at time $t_{k+1}$ is the solution of (8) without noise for initial condition $x(t_{k})=\hat{x}_{k}$. The predicted covariance $\hat{\Sigma}_{k+1}^{-}:=\Sigma(t_{k+1})$ at time $t_{k+1}$ is the solution of (12) for initial condition $\Sigma(t_{k})=\hat{\Sigma}_{k}$. The actual computation of the prediction step used for the tracker algorithm is detailed in §V-E.

Let $\{e_{k}\}$ be a sequence of events. The natural measurement available for event blob target tracking is the location $\xi_{k}$ of each event $e_{k}=(\xi_{k},t_{k},\rho_{k})$, available at discrete asynchronous times $t_{k}$. The associated generative noise model for the raw event location measurement is given by (4). However, this model cannot be used directly in the proposed event-blob target tracking filter since the square-root uncertainty $\Lambda_{k}$ in the generative noise model is itself a state to be estimated by the filter (i.e., determined by variables $\lambda_{k},\theta_{k}$ in $x_{k}$). Instead, we will derive a pseudo measurement with known measurement covariance and use the theory of Kalman filtering with constraints [48].

Define a measurement function by

where $H$ is a function of the state $x_{k}$ parametrised by the measurement $\xi_{k}$. Recalling (4), then by construction

where recall $\eta_{k}\sim\mathcal{N}(0,I_{2})$ is an independent and identically distributed (i.i.d.) Gaussian process with known covariance. In particular, the expected value $E[\Lambda_{k}^{-1}(\xi_{k}-p_{k})]=0$.

Define a new measurement model

with pseudo measurements $y_{k}\equiv 0\in\mathbb{R}^{2}$. The generative noise model (15) with measurements $y_{k}=0$ has known stochastic parameters $\eta_{k}\sim\mathcal{N}(0,I_{2})$ and can be used in a Kalman filter construction.

The measurement (15) is insufficient to provide observability of the full state $x_{k}$. Intuitively, this can be seen by noting that a decreasing measurement error $\|y_{k}-\Lambda_{k}^{-1}(\xi_{k}-p_{k})\|=\|\Lambda_{k}^{-1}(\xi_{k}-p_{k})\|$ can be modelled either as indicating $p_{k}$ should be moved towards $\xi_{k}$ or that $\Lambda_{k}$ should be increased (see Figure 16 in the ablation study Section XI-B).

To overcome this issue, we introduce an additional pseudo measurement specifically designed to observe the shape parameter $\Lambda_{k}$. The approach is to construct a chi-squared statistic from a buffer of prior state estimates that allows estimation of the shape parameters $\lambda_{k}=(\lambda^{1}_{k},\lambda^{2}_{k})$.

Assume that the extended Kalman filter has been operating for $n$ asynchronous iterations, with $n$ events received for the target. During initialisation we will run a bootstrap filter to generate an initial $n$ time-stamps ${t_{k}}$ to enable a warm start for the full filter. For a fixed index $k$ and window length $n$, let $\hat{p}^{-}_{k-j}$ and $\hat{\theta}^{-}_{k-j}$ denote the state-estimate filter predictions based on measurements $\xi_{k-j}$, where $j=(1,\ldots,n)$. Define $\hat{P}_{k}^{-}$ to be the buffer of state prediction estimates

and $\Xi_{k}$ to be the buffer of the associated events

For $j=(1,\ldots,n)$, define

where $(\lambda^{1}_{k},\lambda^{2}_{k})$ are the principal correlations at time $k$, while the prior filter states are used to rotate the image moment to the best estimate of its orientation. Note that we do not use all the available data at time $k$ for the estimation of $\hat{\theta}_{k-j}^{-}$ since doing so would introduce undesirable stochastic dependencies.

Let $\beta\geq\|\Sigma^{-}_{\hat{p}_{k-j}}\|_{2}$ be a constant estimate of an over-bound for the uncertainty in $\hat{p}_{k-j}$. Define

Note that $E[\hat{p}_{k-j}^{-}]=p_{k-j}$ and $E[\hat{\Lambda}_{k-j}]=\Lambda_{k-j}$. Replacing the estimates in (19) by their expected values yields a scaled version of the event generation model (4). Thus, the primary contribution to the uncertainty in $\chi$ will be a Gaussian distribution $\mathcal{N}(0,\frac{1}{1+\beta}I_{2})$. Variance in $\hat{p}_{k-j}^{-}$ is bounded by the estimated parameter $\beta>0$. Since this variance has a similar structure to the variance in $\xi_{k}$ it can be approximated as contributing additional uncertainty $\frac{\beta}{1+\beta}I_{2}$ to $\chi_{k-j}$. The additional uncertainty in $\hat{\theta}_{k-j}^{-}$ is only present in the rotation matrix $\mathcal{R}(\hat{\theta}_{k-j}^{-})$ and its contribution to uncertainty in $\chi$ is negligible. Thus, the uncertainty in $\chi$ can be modelled by

Define a new virtual output function $G$

with dependence on $x_{k}$ through $\hat{\Lambda}_{k-j}^{-1}$ in $\chi_{k-j}$ (cf. Equation (19)). Since $\chi_{k-j}$ are normally distributed with unit covariance, the output function $G(x_{k};\hat{P}_{k}^{-},\Xi_{k})$ is the sum of the squares of $n$ independent 2-dimensional Gaussian random variables. Consequently, it follows a chi-squared distribution with an order of $2n$. According to the central limit theorem [49], as $n$ becomes sufficiently large, the chi-squared distribution converges to a normal distribution with a mean of $2n$ and a variance of $4n$. Based on this, we propose a pseudo measurement model

with pseudo measurements $z_{k}\equiv 2n$.

This development depends on the bound $\beta\geq\|\Sigma^{-}_{\hat{p}_{k-j}}\|_{2}$. In practice, the uncertainty in $\hat{p}_{k-j}$ is often much smaller compared to the uncertainty in the event position $\xi_{k}$ associated with the actual shape of the object. Typically $\Lambda_{k}$ is 100 to 1000 larger than $\Sigma^{-}_{\hat{p}_{k}}$, corresponding to a $\beta\in[10^{-2},10^{-3}]$. Although this bound makes little difference in practice, it provides theoretical certainty that the filter will not become overconfident.

Define the combined virtual measurement

The expected value is $m_{k}=(0,0,2n)^{\top}$ for all $k$ by the definitions of $y_{k}$ and $z_{k}$. Note that $y_{k}$ is a function of the (current) event $\xi_{k}$ and (current) state $x_{k}$, while $z_{k}$ is a function only of the (past) events and the (past) filter estimate before the update step in the Kalman filter. It follows that the random variables in $y_{k}$ (15) and $z_{k}$ (21) are independent. This leads to the generative noise model

where $H(x_{k};\xi_{k})$ and $G(x_{k};\hat{P}_{k}^{-},\Xi_{k})$ are independent.

We linearise the two non-linear measurement models $y_{k}$ and $z_{k}$ (15)-(21) by computing the Jacobian matrices of the $H(x_{k};\xi_{k})$ and $G(x_{k};\hat{P}_{k}^{-},\Xi_{k})$ measurement functions. The linearised observation model $C_{k}$ for the state $x_{k}$ is written

The detailed partial derivatives are provided in Appendix XI-A.

From here, we follow the standard extended Kalman filter algorithm to compute the pre-fit residual Kalman filter gain $K_{k}$ to be

Finally, the updated state estimate and covariance estimate are given by

Following the precedent set by our primary comparison algorithms (ACE [50] and HASTE [23]), we use manual selection (click on screen) or predefined locations for specific experimental studies to initialise tracks in the suite of experimental studies documented in Section VI. In practice, the AEB tracker would be initialised using detection algorithms methodologies tailored to specific applications, such as high speed particle tracking [12] or eye tracking [14], which require specialised detection methodologies. There are a range of potential algorithms in the literature that could be used for this purpose [20, 16]. Further discussion of automated track initialisation is beyond the scope of the present paper.

We initialise targets with zero initial linear and angular velocity, and the initial target size is set to a value at least two times larger than the maximum expected blob size to make the initial transient of the filter more robust.

In addition to the initial states, the proposed filter requires estimates of process noise covariance and an initial prior for the state covariance. The tuning of these parameters follows standard practice for Kalman filter algorithms, allowing for data-based methods like Expectation Maximisation [51, 52]. Alternatively, manual tuning of process noise and prior covariances can be performed to align with the specific application and anticipated motion. The covariance prior for blob location and size is determined based on an estimation of the error in the blob detection methodology. The covariance for linear and angular velocity is set in an order of magnitude larger than the position estimate to allow the filter to quickly converge for these states. In cases of highly erratic motion, the value of $Q$ can be increased to model stochastic variation in the linear optical flow of the target. For less erratic motions, $Q$ can be decreased to enhance the algorithm’s resilience to outliers and noisy data. The better the predictive model of the motion of the target is, the smaller the value of $Q$ can be made. In particular, when an IMU is present and ego-motion of the target for camera rotation can be estimated, this will tend to allow for significantly smaller $Q$ and improve outlier rejection and robustness of the tracking performance. This is important since the ego-rotation of the camera generates the most spurious outlier events from the background texture. Default parameters for the covariance prior and $Q$ that we have found to work well across a wide range of scenarios is documented in the companion code.

Data association in visual target tracking is arguably one of the most difficult problems. In this paper, we take a simple approach that considers all events in a neighbourhood of the event blob to be inliers. As long as the signal-to-noise ratio of the number of events generated by an event blob target to the background noise events is high enough, the disturbance from background noise events will be small and the filter will maintain track. Furthermore, any background noise events that are spatially homogeneous, such as an overall change of illumination, will be distributed in the developed pseudo measurements and will have only a marginal effect on the filter tracking performance. It is only when two objects generating a large number of events cross in the image that the algorithm may lose track. This may be due to two targets crossing in the image, such as when the tail lights from one car occlude the tail lights of another car. Or when the ego-motion of the camera causes a high-contrast background feature to cross behind a low-intensity target. Measuring camera ego-motion by an IMU and predicting the target’s position in the corresponding direction helps to address the second case.

There remains the question of choosing the neighbourhood in which to associate inliers and outliers. The size of this neighbourhood must be adapted dynamically to adjust for the changing size of the target in the image. Our data association approach uses the nearest neighbour classifier with a dynamic threshold $\sigma$. That is, any event $\xi_{k}$ that lies less than $\sigma$ pixels from the predicted state estimate $\hat{p}^{-}_{k}$ is associated with the target and used in the filter. The threshold $\sigma$ is chosen as a low-pass version of $\max(\hat{\lambda}^{1}_{k},\hat{\lambda}^{2}_{k})$. In continuous-time, this low-pass filer is written

where $\alpha$ is the filter gain and $b$ is the desired ratio between the distance threshold and the estimated target size. Integrating (29) over the time interval $(t_{k}-t_{k-1})$ for $\max(\hat{\lambda}^{1}_{k},\hat{\lambda}^{2}_{k})$ constant yields

The gain $\alpha$ is chosen based on the expected continuous-time dynamics of the target in the image.

Our AEB tracker tracks multiple targets independently. Each target is updated as an individual state asynchronously as events arrive. This approach allows us to maintain separate and accurate tracking for multiple blob-like targets. When a target is detected a separate state model $x=(p,v,\theta,q,\lambda)\in\mathbb{R}^{8}$ and the corresponding matrices (e.g., $Q$, $\Xi$, etc.,) is initialised for each new target as discussed in Section V-A. Events are processed asynchronously and assigned to separate targets using our data association method as discussed in Section V-B. These events are used to update the state of the associated target through our proposed algorithm. Any event not associated with a target is discarded.

Flickering targets such as LED lights or non-flickering drones with spinning rotors trigger events that can be approximated to a spatial-temporal Gaussian event blob (Figure 2). However, a moving non-flickering object triggers a bimodal distribution of event data where the leading edge events have one polarity while the trailing edge events have the opposite polarity, depending on the intensity contrast of the blob to the background (see Figure 4). In this case, an unmodified version of the AEB tracker algorithm tends to track either the leading or trailing edge of the target and may sometimes switch between the two edges causing inconsistent velocity estimation. This behaviour only happens for non-flickering blob targets and arises because of the bimodal distribution of events.

To overcome this, we add a simple polarity offset parameter $\Delta=(\Delta_{x},\Delta_{y})^{\top}$ that models the pixel offset, with respect to a central point, between the leading and trailing edge polarity of the observed blob. We expand the state model $x_{k}=(p_{k},v_{k},\theta_{k},q_{k},\lambda_{k},\Delta_{k})$ to include the offset parameter. We add rotation dynamics with stochastic diffusion for the offset parameter to the state dynamics (8)

recalling $J$ from (6) and $\Delta_{k}=\Delta(t_{k})$. The offset is used to compensate the event position $\xi_{k}$ in the first pseudo measurement function

where $\rho_{k}$ is the event polarity at event timestamp $t_{k}$. Note that the polarity offset will converge to the corresponding direction as the velocity estimate if the blob is a light blob on a dark background and to the opposite direction to the velocity if it is a dark blob on a light background in order to balance the event polarities correctly. An analogous change is applied to the second pseudo measurement function (21).

Since events are triggered in microsecond time resolution, the continuous-time prediction step in Equation (9)-(12) are well-approximated by a first-order Euler integration scheme. For each timestamp $t_{k}$ drawn from the events associated with a specific blob, define the time interval

The predicted estimate of the prior state $\hat{x}_{k}^{-}$ and state covariance $\hat{\Sigma}^{-}_{k}$ are computed by

The AEB tracker includes two main steps: data association and the extended Kalman filter. The data association step has a linear event complexity of $\mathcal{O}(N)$ for $N$ input events, effectively passing only a subset of $M$ events to the extended Kalman filter, where $M\ll N$. The extended Kalman filter operates with a time complexity of $\mathcal{O}(nM)$, where $n$ denotes the buffer length in (16)-(17). In practice, the buffer length is typically set to $n<10$, and it is worth noting that $N$ remains significantly larger than $nM$. As a result, our algorithm maintains an overall linear event complexity of $\mathcal{O}(N)$.

Since the data association step only scales linearly with the event rate and this process significantly reduces the number of events used in the Extended Kalman filter, the algorithm is highly efficient for real-time implementation. The efficiency and asynchronous processing of the AEB tracker make it well-suited for implementation on FPGA. Such systems have the potential to be particularly advantageous for high-performance embedded robotic systems.

For the ‘Fast-Spinning’ indoor experiment, we record high-resolution events using a Prophesee Gen4 pure event camera (720$\times$1280 pixels). We evaluate the tracking performance of different event-based corner and blob tracking algorithms at varying speeds. To achieve this, we have constructed a spinning disk test-bed (Figure 6) with adjustable spin rates, and we use a black square on a white disk as the target to provide corners for the corner detection algorithms to operate.

For the ‘Fast Moving Camera with Ego-motion’ experiment, we record data using a hybrid event-frame DAVIS 240C camera (180$\times$240 pixels) that outputs synchronised events, reference frames and IMU data. Note that we do not record frame data because the hybrid DAVIS cameras usually generate a large number of noisy events at the shutter time of each frame when operating in hybrid mode [55].

In this section, we evaluate our method against the most relevant state-of-the-art object-tracking algorithms that use only events as input, operate efficiently (preferably asynchronously), and require no training process. The comparisons include ACE from [50], four HASTE methods from [23], the generic object tracking method from the industrial Prophesee software [53] and the Rectangular Cluster Tracker (RCT) function from the jAER software [54]. In addition, we broaden our scope by comparing with event-based corner detection methods, including the popular asynchronous corner tracking algorithm Arc* [20] and the widely used benchmark method EOF [16]. Note that we limit our comparisons to event-based methods due to the extreme nature of our targeted scenarios in this paper. Frame-based cameras lack the capability to capture information for tracking fast motion or in poor lighting due to motion blur and jumps in blob location between frames (Fig. 6). All algorithms (where velocity is estimated) use a constant velocity model, that is, linear trajectories for the blobs. We do not provide algorithms with prior knowledge that the target trajectory is circular. This mismatch makes the dataset challenging at sufficiently high optical flow.

The main panel of Figure 5 plots the $x$-component of the blob position for each of the different algorithms during an experiment where the optical flow increases from $\sim$100 pixels/s to more than 11000 pixels/s over 90 seconds. The graphic above the main panel plots the period of time for which each algorithm provides reliable tracking of the blob. The plotted line indicates the period up to which the algorithm loses track the first time. The solid dot terminations indicate the end-of-life of the target track while the open circle terminations to the dotted line indicates that future tracks last less than one cycle. The bottom sub-figures highlight zooms of short windows of data to demonstrate key behaviour of the algorithms. The first window, at $\sim$1000 pixels/s, captures the moment when the ACE [50] and the four HASTE algorithms [23] begin to lose track. These algorithms have enduring state estimates and since the blob is travelling in a circular trajectory they often recapture the blob at a later moment as it recrosses the location of the state estimate. This intermittent tracking persists for some time but tracking failures become more regular and eventually the estimates drift from the circular route of the blob and tracking is lost. We also note that ACE [50] and HASTE [23] depend heavily on good identification of the feature template and this is only reliably possible at low speeds. The second figure captures the moment that the Prophesee [53] method loses track, while the third window captures the moment that the jAER [54] loses track. The Prophesee [53] was run with default parameters while the jAER [54] was hand-tuned working with the author (jAER [54]) to achieve the best performance on the test dataset possible. Both these algorithms delete target hypotheses when they lose track. This is indicated with a terminal point on the track in the zoomed windows. Both these algorithms also include detection routines and initialise new targets, which (if they happen to be the correct blob) are able to maintain track for a few additional cycles. However, the first tracking failure indicates the onset of non-robustness in the tracking. Both algorithms regularly lose track for optical flow higher than that indicated in the graphic above the main plot. Our proposed method continues to track stably up to the physical limitations of the experimental platform ($>$ 6 rev/s). Table I presents the maximum optical flow that each algorithm could reliably maintain track.

In addition to the position estimation, jAER [54] and our proposed method also estimate the target velocity. In Figure 7, we plot the velocity estimates $V_{x}=\dot{x}$ and $V_{y}=\dot{y}$ of the target by the jAER [54] and our method. The figures illustrate the linear increase in peak velocity of the oscillation to around 11000 pixels/s in both the $x$ and $y$ directions. The two zoomed-in plots show the detailed behaviour of the velocity estimates of both algorithms for $\sim$9000 pixels/s and $\sim$11000 pixels/s and demonstrate the stability of our algorithm in comparison to the much noisier estimates provided by the jAER algorithm [54].

Apart from the blob tracking methods, we have also provided qualitative comparison studies with event-based corner detection methods, including the popular asynchronous corner tracking algorithm Arc* [20] and the widely used benchmark method EOF [16] at $\sim$1000 pixels/s and $\sim$9000 pixels/s is shown in Figure 8. We mark the corners detected by EOF [16] by green filled dots, Arc* [20] by magenta circles and the position estimated by our proposed method by an orange star. The two compared algorithms mostly detect corners around the target at low speed but fail at high speed and mainly detect the centre of the spinning test bed. Neither algorithm is able to provide reliable target tracking at high optical flow for this simple scenario.

This section evaluates the effectiveness of compensating for camera ego-motion using IMU data in our tracker. The experiment involved tracking the black shapes on a white background in an image where the camera was rapidly shaken by hand. The dataset is similar to, and we used the same shapes image, as the popular dataset proposed by [56], but with more aggressive camera motions, and all targets remain in the camera field-of-view. Maximum optical flow rates are $\sim$1000 pixels/s with abrupt changes in direction.

In Figures 9 and 10, we compare the tracking trajectory of our AEB tracker with and without including camera ego-motion using IMU data. Brighter colours represent trajectories where no IMU data is available, while darker colours correspond to the case where gyroscope measurements from a camera-mounted IMU are available and used to predict ego-motion in the filter as outlined in Section IV. Figure 10 demonstrates that the ‘no IMU’ trajectories lag slightly behind those of the ego-motion compensated model. This effect is also visible in Fig. 9 where the position estimates of the ‘no IMU’ filter are lagging the latest event data slightly more than the estimates obtained using IMU to compensate ego-motion. Compensating camera ego-motion also reduces overshoot when the camera motion changes direction at high speed (Fig. 10). The combination of these effects makes the ‘no IMU’ case slightly less robust and we were able to induce a tracking failure when the camera was shaken at nearly 6Hz with optical flow rates between $\pm$1000 pixels/s. The authors note that although including camera ego-motion is clearly beneficial, the performance of the algorithm without ego-motion compensation is better than expected and is sufficient for many real-world applications.

Apart from recording event data using a Prophesee camera, we also use a FLIR RGB frame camera (Chameleon3USB3, 2048$\times$1536 pixels) to capture high quality image frames for the real-world case study. The two cameras were placed side-by-side and time-synchronised by an external trigger (see Figure 13). The RGB camera frames were re-projected to the event camera image plane using stereo camera calibration. The re-projection mismatch due to the disparity between the two cameras is minor for far-away outdoor scenes. Both datasets include high-speed targets and challenging lighting conditions, allowing a comprehensive evaluation of tracking and shape estimation, as well as for downstream tasks such as time-to-contact estimation and range estimation. For the flying quadrotor dataset, we used a small quadrotor equipped with a Real-Time Kinematic (RTK) positioning GPS (see Figure 13).

In this experiment, the event camera system was mounted on a tripod (held down by hand in the foot well) inside a car driving on a public road at night. The proposed algorithm is used to track the tail lights of other vehicles, and we demonstrate its performance by using the output states to estimate time-to-contact (TTC) at more than 100kHz.

We demonstrate the performance of the proposed AEB tracker by tracking a flying quadrotor under various camera motions and lighting conditions shown in Figure 12. The average speed of the drone in the image is around 150 pixels/s, with some fast motions, such as fast flips, reaching speeds of 400 pixels/s. Figure 12a)-b) demonstrates the tracking performance in well-lit conditions, while data shown in Figure 12c)-d) demonstrate the performance in poorly lit conditions. The RGB images and event reconstruction images are provided for reference only and the filter is implemented on pure asynchronous event data. The red ellipses superimposed on the event reconstruction images illustrate the position, size, and orientation of the quadrotor as estimated by our AEB tracker. The minor axis of the ellipse aligns with the drone’s heading direction. Please also refer to the supplementary video to see the tracking performance.

The first two columns in Figure 12 show the reference frames and event reconstruction recorded with a stationary camera. In these conditions, the majority of the events are triggered by the movement of the quadrotor since it is the only object in the scene that changes brightness or motion. Therefore, the event reconstruction frames mainly consist of a gray background with an event blob triggered by the target. With good lighting conditions and a simple background, the target is easily distinguishable in both frames and event data (see Figure 12a). The advantage of our proposed algorithm over frame-based algorithms becomes more apparent in low-light conditions or environments with complex backgrounds. As shown in Figure 12b), the complex background in the frames makes tracking using RGB frames extremely challenging when the quadrotor crosses trees. In contrast, our AEB tracker can still track the quadrotor successfully. In the poorly lit sequences in Figure 12c)-d), the RGB frames only capture a faint green light emitted by an LED on the quadrotor, making tracking nearly impossible. On the other hand, the event camera can still capture unique event blobs of the quadrotor even in low light or in front of a complex background, providing sufficient information for our AEB tracker to achieve reliable tracking performance even in these challenging scenarios.

In the third column in Figure 12, when the camera is also moving, the background triggers many events, posing a challenge for tracking. The dynamic nearest neighbour data association approach (cf. Section V-B) selectively feeds events to the filter based on their proximity to the centre of the target and estimation of the target size. Only events within a dynamic threshold from the target’s centre are considered, effectively rejecting the vast majority of background noise. Although overall the tracking performance is excellent, the present algorithm can still lose track if substantial ego-motion of the camera occurs just as the target crosses a high contrast feature in the scene. This is particularly the case if the camera ego-motion tracks the target, since this reduces the optical flow of the target and increases the optical flow of the scene. In such situations, a scene feature blob can capture the filter state and draw it away from the target. Interestingly, this scenario is almost always related to the rotational ego-motion of the camera, since such motion generates the majority of optical flow in a scene. The corresponding filter velocity of the incorrect blob track could then be computed a-priori from the angular velocity of the camera, and this prior information could be exploited to reject false background tracks. Implementation of this concept is beyond the scope of the present paper but is relevant in the context of the future potential of the algorithm.

The partial derivatives of $H(x_{k};\xi_{k})$ with respect to $p_{k}$, $v_{k}$ and $\lambda_{k}$ are given by

The partial derivatives of $G(x_{k};\hat{P}_{k}^{-},\Xi_{k})$ with respect to $p_{k}$, $v_{k}$ and $\lambda_{k}$ are given by

where

In this section we present an ablation study to evaluate the effectiveness of our second measurement function $z_{k}$ in the combined measurement function $m_{k}$ (22) and experimentally demonstrate that adding $z_{k}$ enhances the observability of the shape parameter $\lambda$, as discussed in Section IV-B.

We use the ‘Fast-Spinning’ dataset for this comparison. Figure 16 illustrates the size estimation and tracking performance of using the combined measurement function $m_{k}$ in (a) and using only the first measurement function $y_{k}$ in (b). The background intensity is generated using high-pass filter [57] for display only. The estimated target position and size are marked by red ellipses.

We initialise the target position and size at timestamp $t_{0}$. After 0.2 milliseconds, the ellipse computed by the combined measurement function in (a) quickly shrinks to the corresponding target size, while the ellipse of using the single measurement $y_{k}$ in (b) remains almost unchanged. Subsequently, the ellipse in (a) further shrinks to the actual target size and the filter tracks the moving target accurately. However, the ellipse in (b) using only the first measurement function $y_{k}$, tracks the target but demonstrates a gradual increase in estimated size and eventually loses track.