We use the event-based OF computation method proposed in Benosman et al. (2014) which is known for its robustness and its algorithmic efficiency (see Clady et al., 2014, 2015; Mueggler et al., 2015b). More bio-inspired event-based OF computation methods such as Brosch et al. (2015) and Orchard and Etienne-Cummings (2014) can be used but they are computationally more expensive.

A function Σ_e that maps to each p the time t is defined locally:

Applying the inverse function theorem of calculus, the vector ∇Σ_e measures the rate and the direction of change of time with respect to space: it is the normal optical flow, noted v=(vx,vy)T, such as:

This equation could be defined assuming that the surface described by the visual events (generated by a moving edge) in the space-time reference frame (XYt)^T is continuous. This assumption is validated through a regularization process proposed in order to locally estimate this surface as a spatiotemporal plane (fitted directly on the local event stream). In this work the implementation proposed in Clady et al. (2015) has been chosen because it proposes mechanisms to automatically adapt the temporal dimension of the local neighborhood to the edge's dynamics, and to reject estimations of optical flow probably wrong and due to noise. This algorithm allows us to consider a function that associates for each valid visual event e∈E, a so-called visual motion event, noted v_e, such as:

where (v, θ)^T corresponds to the intensity (i.e., speed) and the direction of the normal visual flow.

Remark 1. Note that the polarity of visual events is not conserved by the function (Equation 2). Indeed, in the applications proposed in this article, it is not useful to “memorize” if the visual flow has been computed on a positive or negative event stream. If required, the feature can be augmented in order to distinguish the distribution along “positive contours” from the one along “negative contours.”

As we said, the feature corresponds to the estimated distribution of the optical flow along the (local or global) contours in the visual scene. This distribution is evaluated on a grid-based sampling in the polar reference frame of the visual flow, such as it is subdivided into an interval set {vl¯}l={(θl¯,vl¯)T}l where θl¯ is an angle based interval and vl¯ is an intensity based interval. Such a discretization of the velocity subspace is consistent with biologic observations about orientation (cf. Hubel and Wiesel, 1962, 1968) and speed (cf. Priebe et al., 2006) selectivity in V1 cells and human psychophysical experiments about speed discrimination as in Orban et al. (1984) and Kime et al. (2014, 2016). Here, we parametrize the grid sampling mostly according to these biologic observations and human psychophysical experiments. However, its ranges and precisions could be set in relation with targeted tasks, optimizing them according to given performance criteria. We define the centers {θl}l of the angle intervals such as: θl∈[0,…,2πiNθ,…,2πNθ−1Nθ], with i ∈ [0, N_θ − 1]; 2πNθ is the length of the interval and thus the angular precision of the grid. With N_θ = 36, we barely reach the precision (~10°) observed for V1 simple cells (see Hubel and Wiesel, 1962, 1968). For the velocity intensity, we propose a non-regular speed-based sampling, where {v^l} are the centers of the speed based intervals on a logarithmic scale. The sampling is then operated such that vl∈[vmin,…,vminγi,…,vminγNv−1], with i ∈ [0, N_v − 1] and γ = 1 + ϵ_v (ϵ_v > 0). This discretization strategy ensures an a priori constant relative precision in speed estimation: Δvlvl≈ϵv. Setting ϵ_v to 0.1 will barely correspond to the relative speed-discrimination threshold (10%) observed in human psychophysical experiments (see Orban et al., 1984; Kime et al., 2014, 2016). v_min has been fixed to 1pixel.s⁻¹ and N_v to 73 in order that vmax=vminγNv−1 is close to 1000pixels.s⁻¹, i.e., inversely close to the temporal precision of the visual events, estimated over 1 ms (cf. Akolkar et al., 2015). Motions with intensities less than v_min are then discarded: they are assumed as belonging to static or faraway objects in the background visual scene. Motions with intensities higher than v_max are also discarded because noise associated to their computation can a priori be considered as too high.

Finally, the feature, noted F∈F, is defined as a matrix corresponding to this grid, and associated to a spatiotemporal point (p, t)^T of the retina (or to the entire visual scene for a global approach), and computed as:

where:

w_t is a temporal exponentially decay function (or kernel), inspired by the synchrony measure of spike trains proposed in van Rossum (2001), such that: (4)wtl(t−tj)=H(t−tj)exp(−αvl(t−tj))

where H(·) is the Heaviside step function and α parametrizes the global decreasing dynamic. In our experiments (see Sections 3 and 4), we fixed α to 0.8, i.e., close to 1 in order to mostly take into account the current edges while slightly smoothing them in order to make F less sensitive to both noise and missing data. This kernel gives indeed more weight (or a higher probability value) to events generated by current edges, i.e., the events with timings close to t, while also respecting an isoprobabilistic representation of the edges whatever their dynamics, as we will discuss below (see Section 2.3). Of course, other temporal kernels [Gaussian-based in Schreiber et al. (2003),…] can be envisioned, but this one has the advantage of being causal and of leading to an incremental computation of the feature (see Equation 7).

As we said previously, the feature can be incrementally updated at each occurring visual motion event v_i, considering that Fp,0(vl,θl)=1NθNv for all (v^l, θ^l)^T (in order to consider, at time t = 0, an uniform distribution for the considered velocity-space), such as:

Remark 2. The feature works like a voting matrix, i.e., each visual motion event votes for the speed and direction interval it belongs (and its neighboring intervals through the weighting kernel w_v, Equation 6). More visual events there are, more robust the feature will be. Conversely, the feature will be more sensitive to noise in low light or low contrast situations.

In addition the feature F can be related to a probabilistic distribution while normalizing it to sum up to 1, i.e., to divide it with ∑lF(vl,θl).

In the global approach, F_p,t is independent of p; it can then be noted F_t. Note that the feature is noted F (without sub-index) in this article when the application context (local or global approach) is not relevant or obvious.

Another important point to highlight is that the temporal decreasing function w_t (Equation 4) is related to the speed v^l. Indeed, τl=1vl is the time during which an edge travels through a pixel or in other words, the estimated lifetime of its observation at a given location p, as already remarked in Clady et al. (2015) and Mueggler et al. (2015b). Including it as decay factor in the temporal kernel (Equations 4 and 7) provides a more isoprobabilistic representation of the moving edges in F, i.e., depending only of their contrasts whatever their respective dynamics.

In order to concretely illustrate this point, Figure 3 represents two synchrony images I built integrating a visual event stream and with two different strategies for decay factor τ (related to the speed or not), such as, for each occurring visual event e_i, I(p,ti)=I(p,ti−1)exp(−ti−ti−1τ)+δ(||p−pi||) where δ(·) is the Dirac function. The left image (Figure 3A) results from this equation with a constant τ = cst (whatever the dynamics of the edges), and the middle image (Figure 3B) with a speed-tuned τ=1v. As shown in the right image (Figure 3C), which is the subtraction of both previous images without a speed-tuned factor the high-velocity edges (resulting from the moving and forward leg) are over-represented and the low-velocity edges (resulting from the backward leg) are under-represented in the corresponding synchrony image (Figure 3A). The moving edges are more equitably represented in the second synchrony image (Figure 3B) with a speed-tuned temporal kernel and, by extension, in feature F. Results in Section 3.2 show this equitable representation is very important to obtain accurate results.

The proposed strategy is also consistent with biological observations. Indeed Bair and Movshon (2004) showed that the effective integration time of the computations in direction-selective cells changes with stimulus speed; the integration time for slow motions is longer than that for fast motions. This is modeled in Equation (4) as a decay factor inversely proportional to the speed intensity.

The organization of the feature in a polar coordinate frame based grid, greatly facilitates its computation and its update. The representation of the visual motion information into speed and direction coordinates grants that each speed-tuned decay factor can be associated to an element of the grid, and not directly to the velocity associated to the occurring visual motion event. The latter indicates only which elements in the grid have to be incremented. A bio-inspired implementation can be envisioned where visual motion events are conveyed by selective lines (each line conveying only the motion events v_e included in its associated interval, (v,θ)T∈(vl¯,θl¯)T) from a neuron layer computing the optical flow to a leaky integrate-and-fire (LIF) neural layer (cf. Gerstner and Kistler, 2002), in which each neuron could be assimilated with an element of the feature; this selectivity of lines could result from the selectivity of neurons in the first neuron layer.

Indeed the following model (notations are inspired by Lee et al., 2016) can be used to update the membrane potential of a LIF neuron for a given input event (or spike):

where τ_mp is the membrane time constant, w_k is the synaptic weight of the k-th synapse (through which the input event or spike arrives) and w_dyn is a dynamic weight controlling a refractory period (see Gerstner and Kistler, 2002; Lee et al., 2016 for more details). This model is very similar to the incremental updating equation of our feature, Equation (7). The only things missing are the dynamic weight w_dyn and a firing threshold V_th in order to output approximatively the value of the corresponding feature's element as an event stream (or spike train), and then approximatively following a rate-coding model. Here, the refractory period should be set close to 0 (probably as a small fraction of the integration time τl=1vl), in order to allow (quasi-)simultaneous visual events in the neighborhood (i.e., the events generate by the same contour moving across several pixels in the neighborhood) to contribute equitably to the neuron's potential, i.e., the value of the corresponding element of the feature.

For the local approach, a leaky integrate-and-fire neural layer has to be implemented for each pixel; this neural layer collects the visual motion events from the receptive field, Ω_pi (defined as ∣∣ p − p_i ∣∣< 2σ_s) defined by the corresponding bi-variate spatial kernel (Equation 5). This local computation is detailed in Algorithm 1. For the global approach, only one neural layer is required, collecting the visual motion events estimated over the entire retina.

Finally, Figure 4 shows that the distribution of optical flow representation in the global approach (Figure 4C) summarizes the principal motions observed in the visual scene. This property will allow us to propose a machine learning based approach to recognize gestures in Section 4. In the next Section, we will demonstrate that the local version can be also used to detect particular interest points, i.e., corners.

Remark 3. If the photodiode of the retina's pixel is not square as for the ATIS's one (see Posch et al., 2010 and Figure 1A), the frequency of a set of events emitted by a pixel will be not the same when a contour moves horizontally or vertically in the pixel's field of view (contour's speed and contrast are considered equal in both cases), because the contour travels the same surface of the photodiode during different time periods. In this case, keeping a decay factor invariant whatever the direction of the motion will introduce a bias, favoring one direction over another, in F. To avoid this bias, a cone-pixel with an ellipse-based basis (and not a disk-based basis as illustrated in Figure 2) can be implicitly considered in a correcting function α_θ(·) introduced in Equations 4 and 7 (instead of the constant smoothing parameter α); it is depending on the direction θ^l of the visual motion and defined as:

where α ∈ [0, 1] and e=1−(ab)2 is the eccentricity of the ellipse, with a and b the width and the length of the photodiode, respectively. The second term of this equation increases the decay factor in the direction of the principal axis of the ellipse, rebalancing the representation of the moving edges in F.

In the local approach, normalized F_p,t is the distribution of the normal velocities along the contours around the spatiotemporal location (p, t)^T. In an ideal case illustrated in Figure 7, if this location corresponds to a corner location, F_p,t is null execpt around two velocity coordinates, (vⁿ, θⁿ)^T and (v^m, θ^m)^T, corresponding to both normal visual motions of the intersecting edges.

In order to evaluate the detectors, we reproduced one of the experiments proposed in Clady et al. (2015), the one with the most quantitative evaluations. It consists into a swinging wired 3D cube shown to a neuromorphic camera (DVS, see Figure 8).

A complete accuracy evaluation, comparing the results obtained with the geometric-based method given in Clady et al. (2015) and the methods proposed in this article, is provided in Figure 9. The corner events parameters (spatial location and velocity) and the 11 corners ones (obtained with the ground-truth) are compared using different measures of errors. Each corner event is associated to the spatially closest ground-truth corner's trajectory.

In order to propose a fair evaluation, the thresholds used in the different methods have been set in order to detect the same number of corner events (1500) and other algorithms' parameters have been set as the ones proposed in Clady et al. (2015) (see Remark 3.6). The distribution of the corner events per corner's trajectory is shown in Figure 11A. We can observe that the distributions using the geometric-based and the 2-maxima decision based methods are closely similar. However, the one obtained with the velocity-constraint decision based method is unbalanced, with a great number (close to the third of the corner events) of detections around a particular corner, corner number 5. This can be explained by the fact that the proposed method is less spatially precise than the geometric-based one (cf. the curves in Figure 9A and Remark 3.6) and, as we can see in Figure 10, the edges around this corner generated more events than the others because they are generated by “clean” intersecting edges, see Figure 8, and then verifying well the ideal conditions for the optical flow estimation, and because it is a X-junction. It is not the case for the corners number 1, 7, and 11, for example; the high speed of the cube (close to 500pix.s⁻¹, i.e., inversely close to the precision of the event timings) and their badly shaped structures (they correspond to connections between the different wires constituting the cube) make their detection very hard due to the local bad quality of the event streams (in particular, there are numerous missing events as we can see in Figure 11).

Remark 7. Note that accuracy results in Figure 9 concern median evaluations over the 11 ground-truth corners. Each corner is associated to the spatially closest ground-truth corners trajectory. Each set of corner events (associated to a ground-truth corner) is sorted according to one of the evaluation criteria (type of errors). The Y%-most accurate corner events are then selected. Finally, the accuracy median value for this evaluation criterion is computed over all ground-truths corners. So these evaluations are a priori not (or weakly) biased by these differences in distributions.

We can observe that the detectors proposed in this article are influenced by the quantification of the grid; especially in the Figure 9C representing the angular precision of the estimated speed direction. Indeed a lot of corner events have a direction-related precision close to 5°, the half of the direction-related interval length. The velocity-constrainst based decision method is less clearly influenced because it takes into account more elements in the feature (not only the elements with the maximal values, but also their neighboring elements) to estimate the speed.

In addition, Figure 11B shows the detections distribution for both feature-based methods, with or without speed-tuned temporal kernels. In the approaches without speed-tuning, the temporal decreasing factor τ has been fixed as τ=1vmean, where v_mean is the mean velocity computed over all corners and the stream duration (150 ms). Without speed-tuning, some corners are not or not often detected, in particular corners number 6 and 8. They correspond to X-junctions between two intersecting edges with quite different dynamics, because generated by front and back wires. Furthermore, the accuracy performances for the approaches without speed-tuned temporal kernels are significantly lower than the ones with speed-tuned kernels, as shown in Figure 9.

Finally, if we consider that a corner event detection is valid if the distance error is <3pixels, the geometric-based method generates only 2% of false alarms (with a median velocity error around 10% and a median direction error around 3° for the positive detections), while this rate rises to 8% and to 18% for the velocity-constraint decision and 2-maxima decision based methods, respectively (with a median velocity error around 10% and a median direction error around 8°, for both).

We have demonstrated that the proposed feature can be used (in its local approach) to detect corners in event streams. Even if the detectors are slightly less precise and more sensitive to the quality of the event streams than the other method proposed in the literature, our feature-based approaches are more efficient in terms of memory and computation loads.

Indeed the method in Clady et al. (2015) requires to memorize the stream of the visual motion events (see Equation 2) and spatiotemporal extrapolations of them (called “normal events”) and operates quite complex computations between them. In the approach presented in this article, the visual motion events are integrated directly in the neighboring features, and corner detection related computations are operated only using the feature at the spatiotemporal location of the current event. We have measured important differences in terms of computation time between their different implementations; e.g., for the event stream used for the above evaluations, the feature-based approaches are ~10 times faster. Table 1 presents the distribution of mean computation times obtained with the different approaches and over 10 repetitions (for 1500 detections). But as the method in Clady et al. (2015) has been only implemented on Matlab (Matlab2015b), they should be taken with caution; it is indeed known that memory can be poorly managed on Matlab. Measuring the computation time without code lines dedicated to memory management (which is a crucial part of the method in Clady et al., 2015), the gain is still around 40%. While the geometric-based method is only envisioned in Clady et al. (2015) for a real time implementation on massively parallel computers such as the SpiNNaker board (see Furber et al., 2013; Orchard et al., 2015a), the feature-based approaches run in real-time on a standard computer (in C++ on a Intel Core i7-4790K @ 4GHz, using only one core and without any optimization such as integer arithmetic instead of floating point based computations, e.g., Schraudolph, 1999; Cawley, 2000) for weakly complex visual scenes such as the one presented in this study.

Beyond this operational asset, the greatest strength of the proposed feature-based approaches lies in fact that they lead to a solution of the corner detection issue on event streams based on classical event-based neural network models (leaky integrate-and-fire neural network, coincidence detectors, etc.) as it is highlighted in Section 2.3 and Remark 5.

In order to reduce the dimensionality of the feature (it is often required in machine learning, in order to address the “curse of dimensionality” issue) and to provide (global speed- and) scale-invariance property to the gesture representation, F can be transformed into a more compact representation, noted h (h_p,t or h_t, in local or global approaches, respectively) and named hoof-like in reference to the Histogram of Oriented Optical Flow (HOOF) introduced by Chaudhry et al. (2009) in frame-based vision. This transformation consists in summing the intensities of the optical flow vectors with respect to their directions.

From the feature F, hp,t=[hp,t;1,….,hp,t;Nθ]T can be easily obtained:

In the global approach, normalization (to sum to 1) makes the hoof-like feature globally speed- and scale-invariant. Figure 4D represents the histogram of oriented optical flows computed globally on an event stream capturing a walking human (Figure 4A).

We propose a classification architecture where the problem is framed as a Bayes filter, that is estimating the probabilities of gestures recursively over time using incoming measurements, given as the hoof-like features ht0:tk∈H computed globally from every visual events [e₀, e_k].

Then we note the state gi∈G, the gesture (numerated i, i ∈ [1, K]) that the user is currently performing. A state g⁰ is added in G, in order to consider the not-considered gestures or the instants while the user is not performing a hand gesture.

The camera observes the user's action and at each occurring feature estimates a distribution over the current state gtki:

where htk∈H is the observation of the gesture occurring at time t_k.

To estimate this probability, a time update and a measurement update are performed alternately. The time update updates the belief that the user is performing a specific gesture given previous information:

The time update includes a transition probability from the previous state to the current state. As no-contextual information is available here, we assume that an user is likely to perform the same gesture, and at each timestamp has a large probability of transitioning to the same state:

with τ_g set to 150 ms, less than the half duration of shorter gestures. This assumption means that the gesture's certainty slowly decays over time, in the absence of corroborating information, converging to a uniform distribution (even if no event is observed).

The measurements update combines the previous belief with the newest observation to update each belief state, such as:

In order to estimate p(htk∣gtki), we propose a machine learning based approach to compute and select generative models for gesture. It is decomposed into two steps:

For the first step, we collect hoof-like features computed while the users (included in the training database, see Section 4.3.1) performed a gesture gⁱ, i ∈ [1, K]. Then a k-means algorithm is applied on them in order to compute N candidate models, noted m^gi.

Adaboost (Freund and Schapire, 1996) is an iterative algorithm that finds, from a feature set, some weak but discriminative classification functions and combines them in a strong classification function:

where B and b are the strong and weak classification functions, respectively, and λ is a weight coefficient for each b. T is the threshold of the strong classifier B. The principle of the Adaboost algorithm is to select, at each iteration, a new weak classifier in favor of the instances (or features) misclassified by previous classifiers, through a weighting process attributing more influence to misclassified instances.

Note that a threshold value, noted th_B, can be defined (such as the condition in Equation 18 can be written: 2∑s=1Sλs∑s=1Sλsbs≥thB) in order to optimize a particular classification performance. During the learning step, its default value is 1; this means a classification frontier at the middle of the margin (see Schapire et al., 1998). Increasing or reducing its value correspond to moving the frontier closer or further to the positive class, respectively.

In literature, discriminative training of generative models, as we propose here, has been shown as efficient learning methods in numerous applications as object or human detection (Holub et al., 2005; Negri et al., 2008; Wang et al., 2011), face or character recognition (Prevost et al., 2005; Grabner et al., 2007) or for medical purposes (Deselaers et al., 2008; Wang et al., 2015). The proposed classifier based on the training and the selection of generative models in a discriminative way, combines indeed the main characteristics of discriminative and generative approaches: discriminative power and generalization ability, respectively. The latter is in particular very important in our application, when a weak amount of labeled training data is available, see Section 4.3.1.

Following the framework described in Jing et al. (2008), we propose to design weak classifiers as generative ones, associated to each candidate models msgi (s ∈ [1, N]):

where d(·, ·) is the Euclidean distance and θsgi parametrizes the likelihood function f and is computed at each iteration of the algorithm through a maximum-likelihood estimation (taking into account the weights attributed to features).

During training, Adaboost based algorithm tends to select iteratively the most discriminative and complementary models for each gesture. We limit the number of selected models, such as the relative difference between F-measure (computed on training database, see Section 4.3.1) obtained at the corresponding iteration is superior or equal to 95% of its maximum (obtained with a greater number of iterations of Adaboost algorithm). Let us remind that F-measure is defined as 2×precision×recallprecision+recall. Optimizing it means also to determine a number of models for which an acceptable compromise between precision (the ratio of positive detections to instances belonging to performed gestures) and recall (the ratio of positive detections to all instances detected as belonging to gestures) is reached.

The probability p(htk∣gtki) is then estimated as proportional to a measure (∈ [0, 1]) operated between the hoof-like feature and the set of selected models (applying a sigmoidal function to the output of the strong classifier):

with i ∈ [1, K] and thBi is the threshold obtained optimizing the F-measure. The probability associated to not-considered gesture (or no-gesture), noted g⁰, is then defined as:

Figure 12 presents the obtained classification architecture. Finally a gesture's class G_tk at each time is attributed from the distribution of probabilities, defined as:

Remark 8. Even if our implementation is based on a learning process not directly related to neural approaches (essentially due to the limited size of the database), we can observe that the resulting classification architecture could be fully implemented in an event-based framework. Through a rate-coding model, hoof-like features could be computed and transmitted from the leaky integrate-and-fire neural network, corresponding to the feature computation, as evoked in Section 2.3, to neural networks performing their comparison with gesture models (considering maybe another distance than the Euclidean one used here) and outputting positive events when they match; these positive events corresponding to the weak classifier responses (bsi). The coefficients λsi would be then assimilated to synaptic weights. The other operations, in particular involved in Bayes filters, would correspond to feedback lines and basic mathematical operations that can be modeled using precise timing and event-based paradigms as demonstrated in Lagorce and Benosman (2015).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.