The motion of 3D structures is an important information to extract from a scene to build geometric and dynamic descriptions of its content. Such information is also essential to a large set of vision applications such as: virtual reality synthesis, scene segmentation and autonomous navigation. Scene flows are vector fields that map points of a 3D structure to their instantaneous velocity vectors. Because of this close relationship, estimating the scene flow usually implies to estimating the structure and vice-versa.

The Structure From Motion (SFM) is one of the classical computer vision problems that have been largely studied during the past few decades by the machine vision community (Maybank, 1993). However, SFM's high vulnerability to images' noise and to camera calibration errors raised questions regarding its applicability in real-world scenarios (Tomasi and Zhang, 1995). Currently, with the increasing demand for realistic and high definition 3D content, many ready-to-use sensors are now able to provide dense 3D points clouds in real-time (such as: laser range-finders, structured light vision sensors,…). These devices allow to decouple the structure reconstruction from the motion estimation and to focus the effort on motion extraction and its characterization.

To achieve dense scene flow estimation, state-of-the-art techniques estimate depth maps and compute optical flows for each camera separately. In a second stage they combine both to estimate the 3D flow. This approach parametrizes the motion problem on the image plane, i.e., in 2D and is the most commonly found in the existing literature (Vedula et al., 1999; Zhang et al., 2001; Isard and MacCormick, 2006; Wedel et al., 2011). A 2D parametrization is however more prone to discontinuities since a smooth 3D signal may be projected into a discontinuous 2D signal due to occlusions.

In Basha et al. (2013), the depth map and the optical flow are solved simultaneously rather than in a sequential manner, as authors argue, for a better coupling between spatial and temporal information. In Hadfield and Bowden (2014) and Park et al. (2012), the motion flow is extracted and refined directly from the 3D points clouds by using particle filtering or tensor voting techniques. Optical flows are only estimated for comparison purposes or for initial scene flow estimation.

A second requirement for obtaining dense flow estimation is to introduce some form of regularization. For that purpose, one recurrent hypothesis is to assume local rigid body motion and therefore induce local constant velocity, i.e., points on a non-deformable surface will have the same velocity. Regularization is often performed by minimizing an energy function with variational formulations (Zhang et al., 2001; Min and Sohn, 2006; Huguet and Devernay, 2007). Energy minimization has proven to be a successful technique for both 2D and 3D flow parametrization. It is however computationally greedy and it makes it difficult to achieve real-time estimation without embedding a dedicated powerful computational unit (e.g., GPU). Scene flow can also be computed from local descriptors of reconstructed surfaces such as surfel that encodes the local geometry and the reflectance information of the shapes (Carceroni and Kutulakos, 2002). Motion is then estimated in an integrative manner by matching descriptors over time. Several authors adopted the same idea of addressing the scene flow as a problem of characterization and tracking 3D surfaces over time. Varanasi et al. (2008) proposes to describe and track the surfaces by sparse features matching and extend this to a dense estimation using smoothing operations based on the Laplacian diffusion. Patch based techniques have also been used in Popham et al. (2010) and Cagniart et al. (2010) to split complex surfaces into simpler ones. Their matching and relative pose estimation for each patch allows to estimate the scene flow densely.

This paper introduces a new solution to estimate scene flow using properties of 4D (3D space+time) spaces without the need to use luminance. We will show that the use of the time allows to go beyond the conventional framework that relies on the combined use of luminance and depth information (Herbst et al., 2013). The paper is initially intended to operate on high temporal resolution 3D depth information output from a binocular neuromorphic event-based camera stereo rig. As introduced in Rogister et al. (2012) and Carneiro et al. (2013), event-based cameras allow to estimate depth and produce 3D point clouds at unprecedented accuracy (>1 kHz in real-time) at very low computational and energy cost using conventional processing hardware. We will show that the method can be used even in the case of lower temporal resolution and it can be applied to any 3D data such as the ones output from: RGB-D cameras (Khoshelham and Elberink, 2012), time-of-flight range-imaging sensors (Hansard et al., 2012), laser range finder and even conventional camera based systems that are also able to provide robust 3D reconstructions with a reasonable accuracy using optimized implementations. It is however important to notice that beyond the heavy computational and energy requirement, all these techniques rarely exceed frame rates beyond 90 Hz.

We will then show that the use of timed 4D spaces (3D space + time) allow to derive more efficient techniques than state of art techniques. The method assumes locally non-deformable spatiotemporal surfaces swept by 3D moving structures. We show that under such hypothesis, the velocity estimation is reduced to a one dimensional search over ℝ, the set of real numbers, and the dense estimation is directly achieved using local spatiotemporal planes. An additional advantage is its ability to determine velocities collinear to moving edges assuming it is possible to identify local 3D structures across the trajectory. This work can be seen as a generalization of the previous work on the event-based estimation of 2D visual motion flow (Benosman et al., 2014) to higher dimensional spaces.

Biological retinas do not encode visual scenes as collection of static frames, but rather as a continuous stream of asynchronous spikes. Neuromorphic vision sensors replicate partially this mechanism by encoding visual information with high temporal resolution asynchronous streams of events. Since the pioneering work of Mahowald (1992) that built the first retina on silicon, several major improvements have been made for what is now refered to as the “neuromorphic silicon retinas.” One of the most important achievements is the Dynamic Vision Sensor (DVS) (Lichtsteiner et al., 2008), a 128 × 128 pixel resolution sensor which encodes light intensity changes into a stream of asynchronous events. Each pixel responds independently to contrast changes producing ON and OFF events (respectively to increase or decrease in light intensity) at microsecond resolution.

Posch et al. designed the Asynchronous Time-based Imaging Sensor (ATIS) (Posch et al., 2011), a 302 × 240 pixel resolution sensor which measures absolute luminance information when a contrast change event occurs. The sensor provides a 143 dB dynamic range gray-level information asynchronously encoded as the temporal difference of two exposure measurement events. Its typical temporal accuracy is around 1μs. The reader can refer to Delbrück et al. (2010) for a complete review of the existing neuromorphic visual sensors.

We define a 3D event as a 4-components vector (x, y, z, t)^T. It can be increased to 5 components if the luminance information is available. Let us consider a smooth edge C which can be assumed planar within a small enough spatial neighborhood. If the velocity of C is constant, then as time increases, the edge generates a ruled surface S in the direction of the velocity v. The surface can be algebraically defined by the equation:

where p∈C. Figure 1 shows an illustration of such ruled surface.

The velocity vector is according to Equation (1) the directrix of the ruled surface swept by the edge, hence the estimation of v is equivalent to determining the surface's directrix. In addition to Equation (1), if the surface is smooth (i.e., of class C1 at least), we get a second equation satisfied by v:

because the directrix v is contained in the tangent plane Tp (Sommerville, 1934). ∇S refers to the gradient of S. Only the direction of v can be deduced from the two scalar equations since v has 3 components. Its norm can be set arbitrarily to 1. To determine the exact amplitude, additional constraints are required. A possible way to estimate amplitudes is to apply a shape registration technique, since the velocity vector is the vector joining the two consecutive positions of the shape when it moves. We then propose to estimate the velocity in a two steps operation:

To get enough equations to estimate v, we propose to study three surfaces derived from S. Let S1, S2 and S3 be respectively the surfaces built from Equation (1) in each coordinate frame (X, Y, T), (Y, Z, T), and (Z, X, T). Because of the constant velocity hypothesis, we get three surfaces with implicit equations of the form: (3)Sk(i, j, t)=Sk((pipj0)+ t (vivj1))=0, where (i, j) is any pair of elements in {(x, y), (y, z), (z, x)} and k indexes the k^th element of this list e.g., if k = 1, (i, j) = (x, y). This means we are working with the x, y and t components of S.

These surfaces are also ruled surfaces of respective directrices (v_x,v_y, 1)^T, (v_y,v_z, 1)^T and (v_z,v_x, 1)^T and their generatrices are the restrictions of C to (X, Y, T), (Y, Z, T) and (Z, X, T). For the same reason, given Equation (2), we can establish for each Sk the equation:

As illustrated by Figure 2, we now have three geometric constraints, which can be rearranged into a matrix form:

with the convention that Sk,x (respectively y, z) is the partial derivative with respect to x (respectively y,z). To determine v, the ideal case would be to have M invertible i.e., it is full ranked. There is no obvious way to tell from the general expression of M.

Solving Equation (5) for v cannot be done without knowing the analytic equations of Sk, so we propose to apply a local plane fitting to establish the matrix M. The choice of a plane instead of a more complex surface is motivated by the fitting simplicity and its computational cost even though planes give rise to rank-2 matrices M, as it will be shown further.

Let Π₁, Π₂, and Π₃ be the planes that are fitted locally to the surfaces S1, S2, and S3 respectively. They then can be locally expressed using the plane's implicit equation as: (6)Sk(i,j,t)=Π1T(ijt1)=0, where ΠkT=(ak, bk, ck, dk), for 1 ≤ k ≤ 3.

If we derive Equation (6) with respect to each of the spatial and temporal components and for each Sk, then Equation (5) becomes

Under the local plane hypothesis we previously made, it is possible to determine the rank of M. For that purpose, we assume the hypothesis that the edge C is a straight line segment defined by a point p₀, a direction vector u, and parametrized by a real α:

and the equation of S is changed into:

Figure 3 depicts the case where C is a line and the resulting ruled surface S, obtained by sweeping lines in the direction of v is a plane. The vector (u_i, u_j, 0)^T is by construction parallel to Π_k, then: (10)nkT(uiuj0)=0, where n_k = (a_k, b_k, c_k)^T is the normal to Π_k. The three similar equations for the three possible k lead to: (11)Mu=0. This shows u as an element of the kernel of M. u is not the null vector because C is not reduced to a point, thus M is non-invertible and the rank of M is not larger than 2. The rank deficiency of M means we only have two linearly independent scalar equations from Equation (7), however we can still express two of the velocity components as functions of the last one, e.g., v_x:

where Mq=(0det(M)a3b10)=(000).

This last equation shows that q is collinear to u if M is rank 2, hence we deduce from Equation (7) that Mr = (c₁, c₂, c₃)^T.

Remark 1. M's rank can only be reduced to one if there is no spatial translation. The swept structures in the subspaces defined by (ijT) are vertical lines. Such a case is a particular case which is detected when no ruled surface is generated. It does not concern structures undergoing rotations since points not on the rotation axis will have a non null tangential velocity.

Remark 2. Expressing v as a one parameter vector fails if and only if the rank of M is less than 2 i.e., if edges do not generate planes. However, some plane configurations require larger effort to achieve the closed form of v such as the case where the plane is perpendicular to one of the spatial frame axes. For example, when the X-axis is normal to the plane, Equation (12) is not valid as b₁ and a₃ are equal to zero. This problem can be solved by expressing v either as a function of v_y or v_z. In that case, we can see that v_x = −c₃/b₃ and v_z is a function of v_y. The problem of finding v is again reduced to the search for the correct value of one of its component.

As shown in the previous section, from Equation (12), the assumption of local constant velocity motion of straight edges allows to establish a simple linear relation between the velocity vector and the surface swept by the edge points. Estimating the velocity becomes equivalent to identifying the correct real value v_x. This is a registration problem for which we need to initiate the point cloud within a spatio-temporal neighborhood as a given structure. We then translate it according to vectors v, parametrized by v_x. A matching operation is then performed for several sampled values of v_x, the correct v_x is the one producing the smallest matching error at the time and location given by the velocity vector (see Figure 4). The procedure to estimate the velocity via the shape registration is explained in detail by Algorithm 1: the search for v is now a minimization problem of a error cost function E, which is built as explained in the next section.

A local point cloud centered on the event (p1, t1)^T, is temporally consistent in the sense that any of its element will be captured at closely the same time. If in addition, the luminance L of the events is available, then the cloud local rigidity also ensures that L is consistent independently of time. We can therefore state that when the point cloud that moves from p₁ at t₁ to p₂ at t₂, the local geometric structure and the luminance should be preserved. We can formalize the structure matching operation as a minimization of the energy E problem and stated as follows:

where E_S, E_T, and E_L are respectively the geometric, the temporal and the luminance energies. E_T and E_S are minimal as long as the cloud is not deforming when it moves from p₁ to p₂. If the events' brightness is also preserved during this motion then E_L is also minimal.

We define the 3D events cloud S(p_i, t_i) as:

This set contains all 3D points spatiotemporally close to p_i i.e., points within a neighborhood of p_i of radius Δ_s in space and length Δ_t in time. The energy cost associated to each sampled velocity vector for a given point p₀ is computed according to Algorithm 2.

E_L is the sum of the smallest luminance difference between all pairs of (p_i, q_j) and E_S is the mean value of the smallest distances of each p_i to each q_j. It is also called the mean closest point between both points clouds and is a dissimilarity measure often used for example in the Iterative Closest Point (ICP) problem (Besl and McKay, 1992). The correct v is given by the value v_x which minimizes the energy function E:

More elaborate registration techniques to track deformable 3D surfaces may be used for this matching operation. We can mention the most notable ones, Starck and Hilton (2007), Ahmed et al. (2008), and Zeng et al. (2010), that are not using any shapes prior. Accurate registrations are achieved by combined use of several surface features, followed by a coarse to fine scheme. These techniques are however not suitable in their actual form for processing textureless and event-based inputs.

To minimize E with respect to v_x, we also applied a dichotomic search strategy to sample possible values of v_x and match local 3D structure accordingly. Let R = [R₁, R_r] be a real interval that is set large enough at the beginning of the search to make sure it contains ṽ_x. Fixing R large enough is only necessary when no recent past estimations of the velocity have been calculated at p_i, otherwise the length of R is defined from the previous estimation of v_x. To determine precisely ṽ_x, R is subdivided into r equal length intervals and the centers of all intervals give a set of possible values for v_x. The error cost function is computed for each v_x and the interval producing the smallest E is used to update R.

This operation is iterated until E is below a preset threshold and after a minimum number of iterations. This threshold is defined experimentally with the purpose of optimizing the structure matching process by limiting the search to an acceptable matching error. This threshold can be related to the point cloud density and if it is set to zero, then the maximum number of matching iterations is always performed. r is usually set to 5, however it can be larger. Estimation accuracy increases with r but at the cost of longer processing time.

The correct estimation of the velocity is conditioned by the spatiotemporal neighborhood, defined as the spatiotemporal volume of dimensions (Δ_x×Δ_y×Δ_z×Δ_t), in which the 3D point cloud has moved from time t to t + dt. A large neighborhood will allow to find the correct match, but at the cost of processing a large set of data, on the contrary, a too small one will not allow to match the local structures. The spatiotemporal neighborhood must also be resized automatically and dynamically in accordance to the 3D points' velocity. In our implementation, we deal with the problem by adjusting a linear function on the neighborhood size e.g., s_k = (Δ_x, Δ_y, Δ_z, Δ_t)^T is a linear combination of the m previous values s_k−1, …, s_k−m: (19)sk=∑i=1maisk−i, where the coefficients a_i are estimated with a standard linear prediction coding scheme (Durbin, 1959). The value of m is usually set to 5 according to experimental results while the initial value s₀ is deduced from the coarse estimation of the initial velocity i.e., the mean translation between the first two frames. Thus, we have s0=(v0dt,dt)T, assuming v₀ is the initial estimate of the velocity.

The asynchronous 3D flow extraction from points clouds can be achieve by implementing Algorithms 1, 2 and improved if necessary with the optimal neighborhood estimation. This proposed approach does not require clusters of 3D points captured at the same t as one uses to have with frame-based reconstructions, yet it can still be applied if the inputs are frame-based.

Four simulated scenes are synthetized: (Figure 5) a smooth translation of a wire cube at constant amplitude; (Figure 6) a 3D car model undergoing a straight translation at 10 m per second; (Figure 7) the same 3D car model describing a circular motion; (Figure 8) a pure rotation of a sphere at constant angular speed.

For each scene, the velocity flow is computed using the geometric structure information alone (only the 3D points' positions and timestamps are given in the simulation). The energy cost function in Algorithm 2 is reduced to E_s. The flow performance is measured by two quantities conventionally used to validate optical flow, the angular error, which is the angle defined by the estimated normalized velocity vector v~ and the ground-truth v. The angle is given by the inverse cosine of the scalar product of 2 vectors: (20)arccos(v˜Tv/|v˜||v|). This measure has been introduced in Fleet and Jepson (1990) to assess the accuracy of the flow direction. However, the angular error can be biased by large velocity vectors for which the differences in amplitude can be significant and in the same time, the angular errors are small. In that sense the angular error is favoring the large motion over the smaller ones. To compensate for that bias, a second performance measure, the endpoint error, introduced by Otte and Nagel (1994) is conjointly used. This endpoint error is the norm of the difference between the estimated velocity and the real one: (21)|v−v˜|. Both estimated angular error and endpoint error are represented with a color scaled representation (Figures 5–15). For the moving cube, the maximal error occurs at the beginning of the motion and is due to the fitting spatio-temporal neighborhood, chosen as the best compromise for the entire motion.

The results on synthetic data, summarized in Table 1, show the ability of the method to estimate densely and smoothly the velocity field. The rotating sphere is a challenging because the 3D points composing the surface are not spatially uniform. The non uniform acceleration on the sphere also implies non uniform tangential velocity of 3D points on which we fit the local planes. This explains why the velocity estimation is slightly less accurate for rotations. For translations (cube and car), the velocity is estimated with higher accuracy since the direction has a mean angular error of 0.04 rad (with a peak value of 0.1 rad) and a mean endpoint error of 0.8%, with a peak value of 1.2% when normalized by the ground-truth objects sizes (respectively the cube edge length, the sphere diameter, and the length of the car). For rotations, the accuracy has the same order of magnitude: around 0.15 rad and 2.2% for respectively the angular and the endpoint error.

The second set of results is obtained from real scenes showing a moving person in the scene. The 3D point clouds are provided by a Kinect sensor that also measures the RGB intensity. The Kinect provides depth information for every detected pixel. The background pixels representing the room's walls are removed via depth segmentation. In these sequences, the person is a nice example of a deformable target with limbs moving at different non-constant velocities. However, the local constant speed hypothesis holds. It is sufficient to allow a smooth estimation of the scene flow. Scene flows estimations are given as two sets of results. The first one uses only geometric constraints, when the scene luminance is not available for the structure registration operation. The second set uses the additional information brought by the luminance in addition to the geometry.

The flow estimation for each sequence is assessed in two ways:

This subsection provides the 3D scene flow using event-based cameras (DVS) as described in Carneiro et al. (2013). Computed 3D data have a high resolution of 1μs. The input to the scene flow estimation are asynchronous 3D point clouds of a hand closing and opening in front of the stereo rig, while the second sequence is a moving face captured by the same stereo rig. The hand speed is of the order of one meter per second while the face moved slower (several cm per second). The first sequence's results are shown in Figures 12, 13. The events generated by the hand's contours are sufficient to estimate 3D flow estimation. The direction and amplitude are consistent with the motion. The mean morphing error, (8.7%), is at the same order of magnitude than the previous experiments.

The second 3D scene flow estimation from the event-based cameras are shown without providing morphing error for readability reason. The second sequence shows a face moving in front of the event-based cameras. The color-coded flow are shown in Figure 14 with row (a) color-coding the flow's directions while (b), shows the motion amplitudes.

This third set of results obtained with the event-based plane fitting technique (without luminance information) is compared to the method published in Hadfield and Bowden (2014). This paper models a set of moving 3D points using a particle filter that supports multiple motion hypotheses to estimate the 3D scene flow from the 3D points provided by a Kinect.

The sequence and its estimated scene flow presented in Hadfield and Bowden (2014) were kindly provided by the authors and are shown as output without additional processing. Figure 15 shows samples of this sequence along with the estimated 3D scene flow: the velocity amplitude and direction are shown in two separate color-coded representations for the two methods. One can notice that the event-based plane fitting method produces smoother results that are consistant with the scene content, especially the velocity is expected to be maximal at the foot when the kick is accomplished.

We provide the morphing error in Table 3. It shows the error computed for 5 frames taken at some arbitrary regular time interval, using both method. This is consistant with the previous analysis showing that the event-based plane fitting method's performance is higher than the state-of-the-art frame-based technique. Unsurprisingly, the overall estimation accuracy for both methods is lower than the previous experiments mainly because the Kinect's inability to accurately capture fast motions that give rise to blurry images.

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.