Julián Quiroga

Scene Flow Estimation from RGBD Images

Keywords: motion, depth, scene flow, rgbd, variational, semi-rigid

This thesis addresses the problem of reliably recovering a 3D motion field, or scene flow, from a temporal pair of RGBD images. We propose a semi-rigid estimation framework for the robust computation of scene flow, taking advantage of color and depth information, and an alternating variational minimization framework for recovering rigid and non-rigid components of the 3D motion field. Previous attempts to estimate scene flow from RGBD images have extended optical flow approaches without fully exploiting depth data or have formulated the estimation in 3D space disregarding the semi-rigidity of real scenes. We demonstrate that scene flow can be robustly and accurately computed in the image domain by solving for 3D motions consistent with color and depth, encouraging an adjustable combination between local and piecewise rigidity. Additionally, we show that solving for the 3D motion field can be seen as a specific case of a more general estimation problem of a 6D field of rigid motions. Accordingly, we formulate scene flow estimation as the search of an optimal field of twist motions achieving state-of-the-art results.

Estimation du Flot de Scène à partir des Images RGBD Résumé:

Cette thèse aborde le problème du calcul de manière fiable d'un champ de mouvement 3D, appelé flot de scène, à partir d'une paire d'images RGBD prises à des instants différents. Nous proposons un schéma d'estimation semi-rigide pour le calcul robuste du flot de scène, en prenant compte de l'information de couleur et de profondeur, et un cadre de minimisation alternée variationnelle pour récupérer les composantes rigides et non rigides du champ de mouvement 3D. Les tentatives précédentes pour estimer le flot de scène à partir des images RGBD étaient des extensions des approches de flux optique, et n'exploitaient pas totalement les données de profondeur, ou bien elles formulaient l'estimation dans l'espace 3D sans tenir compte de la semi-rigidité des scènes réelles. Nous démontrons que le flot de scène peut être calculé de manière robuste et précise dans le domaine de l'image en reconstruisant un mouvement 3D cohérent avec la couleur et la profondeur, en encourageant une combinaison réglable entre rigidité locale et par morceaux. En outre, nous montrons que le calcul du champ de mouvement 3D peut être considéré comme un cas particulier d'un problème d'estimation plus général d'un champ de mouvements rigides à 6 dimensions. L'estimation du flot de scène est donc formulée comme la recherche d'un champ optimal de mouvements rigides. Nous montrons finalement que notre méthode permet d'obtenir des résultats comparables à l'état de l'art. In the real world everything changes. Changes occur for many reasons but nothing remains the same. Some of these changes are sudden while others are subtle. Some are induced while others are provoked, but always these changes reveal information that was hidden from our initial perception. Visual changes brought by relative motions of the environment are particularly meaningful, allowing the perception of, and interaction with, the surrounding world. For instance, a chaotic motion alerts us to possible dangers. A peculiar way of motion tells us that someone known is approaching. A subtle gesture of the hand is enough to communicate our dislike.

Mots

Being able to understand visual changes is a wonderful ability of human beings, which inspires the development of computer vision systems aiming to improve the way we perceive, and we are perceived, in our environment.

A computer vision system sees the world through a camera, which usually corresponds to a color-sensitive retina providing a planar representation of the scene. Visual changes in the world are perceived as displacements of the color patterns on the retina, and the computation of such those motions enables understanding of a changing environment. Image motion estimate provides powerful cues for visual systems, and has been exploited in a wide variety of applications including autonomous navigation, action recognition, human-machine interaction, threedimensional (3D) reconstruction. For example, by estimating and learning motion patterns, visual systems are able to recognize actions, activities and gestures. By using a single moving camera and tracking a set of key points it is possible to recover the 3D structure of a static scene. Even, the shape of a time-varying 3D surface can be estimated, after a sufficiently long observation. Nevertheless, this planar representation of the scene is only a projection of the 3D world, and therefore the perception and estimation of some motions may become a challenge. For instance, the lack of texture hinders motion estimation, because there is not enough information to disambiguate between all of the possible motions that can explain the observations. Moreover, it is complicated to handle with partial occlusions, which can severely confound motion estimation. For such of challenges, motion estimation using monocular color images is still a very active area of research. There is a physical limitation due to the use of a single view, which prevents the reliable estimation of the full 3D motion field of the scene. If a fully calibrated stereo or multi-view camera system is available, it is possible to estimate both the motion and the structure of the scene. Evidently, when the number of views increases, there is more information available to deal with occlusions and motion ambiguities, but at the same time, such system becomes less affordable and portable.

Recently with the arrival of depth cameras, based either on time-of-flight (ToF) or structured light sensing, it has been possible to have direct measurements of the 3D structure of the scene. Such cameras are also known as RGBD (RGB and Depth) cameras, since they simultaneously provide color and depth information. Access to structure information opens the door to the estimation of motions directly in the 3D scene, using a single view. For this reason, current depth sensors have decreased the requirements of the system needed for this kind of tasks. Moreover, the availability of RGBD representation of the scene demands the reformulation of the 3D motion estimation problem. Usually, when only color information is available, motion is solved to be consistent with the displacement of color patterns. However, for RGBD cameras, both color and 3D structure can be simultaneously exploited. For instance, action recognition is usually based on features obtained from 2D trajectories computed on the image. The simple addition of the depth information to every trajectory significantly improves the recognition performance. However, it is clear that this way depth data is not fully exploited. The combined use of color and depth can yield more accurate trajectories, incorporating at the same time, the 3D structure information. In general, the use of depth data in addition to color can yield better results for computer vision tasks, such as recognition, interaction and reconstruction. Both color and depth data complement each other and can be combined as whole to obtain the most benefit from RGBD images.

Scene flow is defined as the motion field in 3D space, and can be computed from a single view when using an RGBD sensor. Being able to confidently estimate the scene flow makes it possible to directly take advantage of real world motion. In this thesis, we address the problem of scene flow estimation from RGBD images; by rethinking the way color and depth are exploited. The two main questions for scene flow estimation from RGBD images are: how to fully exploit the data?, and which motion model should be used?. We go beyond previous works and propose a scene flow framework that exploits the local and piecewise rigidity of real world scenes. By modeling the motion as a field of twists, our method encourages piecewise smooth solutions of rigid body motions. We give a general formulation to solve for local and global rigid motions by jointly using intensity and depth data. In order to deal efficiently with a moving camera, we model motion as a rigid component plus a non-rigid residual and propose an alternating solver. The evaluation demonstrates that use of the proposed framework achieves the best results in the most commonly used scene flow benchmark. Moreover, through additional experiments we indicate the general applicability of our approach in a variety of different scenarios.

Thesis outline

The following six chapters present the basics, methods and results for scene flow estimation from RGBD images.

Chapter 2. The first part of this chapter introduces two important concepts about the scene: the RGBD representation, based on color and depth data, and the scene flow, which models the motion observed in the scene and is the main concern of this thesis. The second part presents three useful tools for motion representation and estimation. Fundamentals of rigid body motion are reviewed, particularly, the twist motion representation that we use to model the motion in the scene. We summarize the concept of total variation, describing its role and properties as regularizer and presenting some basic solvers. Finally, we present the basics of robust least-squares estimation.

Chapter 3. The chapter analyses the estimation of 3D motion when color and structure of the scene are available. We first define the scene flow estimation problem, pointing out the main challenges of this task and analyzing properties of real scenes that can be exploited to aid the estimation. We then present a survey of the most relevant work done in scene flow estimation, where the advantages and limitations of most known techniques are studied. Simultaneously, we explain how the ideas supported by this thesis go beyond previous methods. Having presented the state of the art for scene flow estimation, we finally give a brief exposition of the main components of our approach.

Chapter 4. In this chapter, we present alternative representations for 3D motion in the image domain. First, we introduce the parameterization of the scene flow as a field of rigid motions, using a twist representation. Here we define a warping function to constrain every twist motion on the image, and which is also provided with a linear version allowing an incremental solution of rigid motion. Alternatively, we present the direct representation of the scene flow as a 3D motion field, using a warping function to exploit RGBD data. Finally, we review other possibilities for motion representation and we conclude with a table containing the three cases considered in this thesis.

Chapter 5. The framework for scene flow estimation from RGBD images is presented. We begin by motivating the semi-rigid constraint of the 3D motion field on the image and then we formulate the scene flow energy based on the twist motion representation. Subsequently, each of the components of this energy is described. By considering the 3D motion representation, we formulate an alternative scene flow energy, which is also described in detail. We then present the minimization of both energies and a coarse-to-fine procedure. Finally, we introduce an approach to estimate the scene flow with a moving camera using the proposed framework.

Chapter 6. This chapter presents the experiments performed using the proposed framework. The first part contains the motion estimation results on the Middlebury stereo dataset. Using this benchmark we compare the accuracy of our approach with the most relevant scene flow methods. Also, we analyze the different components our framework and the alternatives that it supports. The second part of the experimentation considers images from RGBD sensors. We investigate the performance of our method in different and interesting scenarios presenting nonrigid and rigid motions.

Chapter 7. We conclude by reviewing the formulated questions and describing how this thesis contributes to their resolution. We summarize the limitations of this approach and propose some directions for future work. Finally, we conclude by discussing the impact of our work, some of the lessons learned and perspectives on scene flow estimation from RGBD images.

Background

RGBD representation of a scene

Let Ω be the image domain defined as a regular L x ⇥ L y grid of pixels, indexed as x =( x, y), for x =1 , 2,...,L x , y =1 , 2,...,L y . The use of a RGBD camera, e.g., Microsoft Kinect for Xbox or Asus Xtion Pro Live, enables color (RGB) and depth (D) discrete representations of the observed scene in the image domain. For every time t the scene is registered as a RGBD image S t (x)={I c (x),Z t (x)} where images I(x) and Z(x) provide color and depth, respectively, for every pixel x 2 Ω. This kind of camera provides both RGB and depth with VGA resolution (640⇥480 pixels). The RGB image uses 8-bit pixel representation with a Bayer color filter. The depth image uses 11-bit pixels to encode the 3D structure of the scene, using Light Coding by PrimeSense. This technology codes the scene with near-IR light and with an infrared camera that captures distorted light to measure the relative Fig. 2.2: The RGBD image is displayed by showing the gray valued of every pixel having a valid depth measure. We also incorporate information about 3D structure as contour lines, which are defined as a function of the point-wise changes in depth data.

distance of the surface. As is shown in Figure 2.1 b), some pixels do not have a valid depth measure for one of these reasons: i ) a depth out of the sensor range (80cm to 400cm for the Kinect for Xbox ), ii ) an occluded point in the triangulation between IR emitter and IR camera, iii ) a point out of the field of view of the depth sensor but visible for the RGB camera, or iv ) a poor or non reflective surface to IR light.

In this RGBD representation every pixel x 2 Ω, having valid depth measure, represents a 3D point X of the scene, or scene point. For illustration, in some cases we display the RGBD image as a single image, by blending gray values and depth data, as is described in Figure 2.2. Below we present the camera model used to define a bijective mapping between 3D space and image domain.

Camera model

The camera frame is defined with its origin at the optical center and its (X, Y, Z) axes respectively in the direction of the image axes and the optical axis, with positive Z pointing in the direction of sight, as is shown in Figure 2.3.

In the RGBD representation, every 3D point that is visible for the sensor, is parametrized into the image domain Ω. The projection π : R 3 ! R 2 maps a 3D point X =(X, Y, Z) T onto a pixel x =(x, y) T on Ω by:

π(X)= ✓ f x X Z + c x ,f y Y Z + c y ◆ T , (2.1)
where f x and f y are the focal lengths of the camera and (c x ,c y ) its principal point.

The inverse projection π -1 : R 2 ⇥ R ! R 3 back-projects an image point to the 3D space for a given depth z, as follows: Given a RGBD image S(x)={I c (x),Z(x)}, the notations X or π -1 (x,Z(x)) are equivalently used as the 3D representation of image point x.

π -1 (x,z)= ✓ z x -c x f x ,z y -c y f y ,z ◆ T . ( 2

Scene flow definition

Let X t =( X t ,Y t ,Z t ) T be coordinates at time t of a scene point in the camera reference frame. The scene flow, v(x):Ω ! R 3 , is defined as the 3D motion field describing the motion of every visible 3D point between two time steps:

X t+1 = X t + v (x t ) .
(2.

3)

The 3D motion of every visible scene point induces a motion in the image domain, which is defined as image flow and corresponds to the 2D motion field u(x):Ω ! R 2 , given by the projection of the scene flow. Figure 2.3 shows the motion and projection of a 3D point belonging to a rigid surface. An illustration of 2D and 3D motion fields is present in Figure 3.1. Note that the image flow should not be confused with the optical flow, which is defined as the apparent motion of the brightness patterns in the image, since they both may not coincide. For instance, not every 3D motion in the scene yields to a change in the brightness images, e.g., the motion of an untextured surface, and conversely, not every change in brightness is brought by a the motion of a 3D surface, e.g., a change in illumination.

Rigid body motion

On a rigid surface the distance between any two given points remains constant in time regardless of external forces exerted on it. Accordingly, a rigid body motion, or rigid motion, is defined as a motion which preserves distance between points of the moving surface. The group action of a rigid body transformation can be written as T (X)=RX + t,w h e r et 2 R 3 is a translation, and R is a 3 ⇥ 3 matrix belonging to SO(3), the special orthogonal group of rotation matrices in R 3 , defined as:

SO(3) = R 2 R 3⇥3 : RR T = I 3⇥3 , det(R)=±1 .
(2.4)

A rigid body motion can be seen a rotation R followed by a translation t, as is illustrated in Figure 2.5. Using homogeneous coordinates, under the action of rigid motion G, a 3D point X =(X, 1) T is transformed into X0 , according to:

X0 = G X, with G = Rt 0 1⇥3 1 ! 2 SE(3), (2.5)
with SE(3) the special Euclidean group representing rigid transformations in R 3 :

SE(3) = (R, t):t 2 R 3 , R 2 SO(3) = R ⇥ SO(3). (2.6)
The group of matrices, SE (3), is a group with the multiplication as operation, and the identity matrix, I 4⇥4 , as the identity. Therefore, the composition of two rigid motion, G 1 and G 2 , is done by a matrix multiplication, yielding to the rigid motion G = G 1 G 2 . Moreover, every rigid motion G has inverse an inverse given by G -1 . It is important to remark that SE(3) and SO(3) are Lie groups, i.e., groups with a smooth manifold structure. Any Lie group G defines a Lie algebra g through an exponential function. Below we define the exponential function relating the Lie group SE(3) with its Lie algebra se (3). Unlike Lie groups, Lie algebras are linear spaces allowing to understand and solve problems on the Lie group in a easier way.

Twist representation

Since any G 2 SE(3) has only 6 degrees of freedom, a more convenient and compact representation is the 6-parameter twist. Every rigid motion can be described as a rotation around a 3D axis ω =( ω X , ω Y , ω Z ) T and a translation along this axis, giving as a function of a 3D translation τ =(τ x , τ y , τ z ) T , see Figure 2.6. Therefore it can be shown that for any arbitrary G 2 SE(3) there exists an equivalent ξ 2 R 6 twist representation. A twist ξ =(ω, τ ) can be converted into the G representation with the following exponential map:

G = e ξ = I 4⇥4 + ξ + ( ξ) 2 2! + ( ξ) 3 3! + ••• , (2.7)
where ξ = ωτ

0 1⇥3 1 ! , with ω = 0 B B @ 0 -ω Z ω Y ω Z 0 -ω X -ω Y ω X 0 1 C C A .
(2.8)

Every twist vector ξ corresponds to a 4⇥4 matrix ξ 2 se(3), which is the Lie algebra of SE(3), via the hat operator given by the left side of equation (2.8). Similarly, every rotation vector ω is related to a 3 ⇥ 3 skew-symmetric matrix ω 2 so(3), the Lie algebra of SO(3), via the right side of equation (2.8). Alternatively, the exponential map (2.7) can be expressed in a compact form as follows: where the rotation matrix e ω is given by Rodrigues's formula:

e ξ = e ω Aτ 0 1⇥3 1 ! , (2.9) θ k ω τ ξ =(ω, τ ) ←-ω O 3D rotation 3D traslation
e ω = I 3⇥3 + ω |ω| sin |ω| + ω2 |ω| 2 (1 -cos |ω|), (2.10)
and

A = I 3⇥3 + ω |ω| 2 (1 -cos |ω|)+ ω2 |ω| 3 (|ω| -sin |ω|), for ω 6 = 0.
(2.11)

The exponential map (2.7), or equivalently (2.10), goes from se(3) to SE(3).

Correspondingly, for each G 2 SE(3) there exists a twist representation given by the logarithmic map:

ξ = log(R, t)=(ω,A -1 t), (2.12)
where

ω = 1 2 sin θ (R -R T ), with θ = cos -1 ✓ trace(R) -1 2 ◆ , for R 6 = I 3⇥3 , (2.13)
and

A -1 = I 3⇥3 - ω 2 + 2 sin |ω| -|ω|(1 + cos |ω|) 2|ω| 2 sin |ω| ω2 , for ω 6 = 0. (2.14)

Small-rotation model

Any given rotation vector ω encodes the amount, or angle, of rotation in its magnitude, |ω|, and the axis of rotation as the unitary vector ω/|ω|. In the case of the angle of rotation approaches zero, Rodrigues's formula (2.10) yields to:

lim |ω|!0 e ω = I 3⇥3 +ω . (2.15)
Therefore, for a small-rotation rigid motion, the rotation matrix is well approximated as e ω ⇡ I 3⇥3 +ω [START_REF] Fang | Solving three-dimensional small-rotation motion equations[END_REF]. Similarly, if |ω| ! 0w eh a v e A ⇡ I 3⇥3 in equation (2.11) and the full rigid motion can be approximated by: (2.16) which implies that t ⇡ τ . This linear approximation is useful to analyze a small rigid motion in the 3D space and in the image domain. Moreover, equation (2.16) allows to write the scene flow induced by a small-rotation rigid motion as a linear function of the twist. Let ξ be a twist motion for a small rotation, under the rigid motion, the 3D point X is transformed into X 0 inducing the scene flow given by:

e ξ ⇡ I 4⇥4 + ξ,
X 0 -X =([X] ⇥ | I 3⇥3 ) ξ, (2.17)
where the operator [•] ⇥ takes the elements of X to construct a skew-symmetric matrix (the cross product matrix), as is done in the right side of equation (2.8).

Total variation

Total variation was introduce to Computer Vision by Rudin, Osher and Fatemi (ROF) in their seminal work [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] for image denoising. The ROF denoising model can be written as:

min u Z Ω |ru(x)|dx + η 2 ku -f k 2 2 , (2.18)
where Ω ⇢ R 2 denotes the image domain, and f : Ω ! R, the degraded observed image. The right side of the model (2.18)i st h efidelity term, which encourages the restored image u to be close to the input image f in L 2 , the space of squareintegrable functions. The left term is the total variation (TV), which is used a Tikhonov regularization and has the capability to preserve discontinuities (edges) while removing noise and other undesired details, see Figure 2.7. The constant η balances the tradeoff between a good fit to the data and a smooth solution. In the TV term, ru(x) is the gradient of u evaluated at x and its magnitude is given by: The solution of (2.18), which is also known as the primal formulation,i sa challenge, since is non-differentiable at |ru(x)| and is highly nonlinear. For this reason the dual formulation of the ROF-model has also been studied [START_REF] Chan | A nonlinear primal-dual method for total variation-based image restoration[END_REF] and has been shown that the desired denoised image u ⇤ can be solved by:

|ru(x)| = s ∂u ∂x 2 (x)+ ∂u ∂y 2 (x) . (2.19)
u ⇤ = f + 1 η r • w ⇤ , (2.20)
where w ⇤ : Ω ! R 2 is the solution of the dual problem:

min |u|1 kr • w + λf k 2 2 .
(2.21)

The solution of the dual formulation presents its own challenges since it brings extra constraints and its solution is not unique, see [START_REF] Zhu | Duality-based algorithms for totalvariation-regularized image restoration[END_REF] for more details. However, the dual problem was the first to allow an exact solution of the ROF model by using a projection algorithm [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF]. The ROF model (2.18) has been successfully used in other computer vision tasks. Particularly, for motion estimation, TV is preferred to other strategies, since TV: i ) does not penalize discontinuities of the motion, ii ) does not penalize smooth functions and iii ) preserves exactly the location of discontinuous edges, see Figure 2.8.

Discrete ROF model

The discrete total variation of a scalar function u : Ω ! R is defined by: (2.22) where the components of the gradient ru =( u x ,u y ), are computed using forward differences, see [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF]. Accordingly, the discrete ROF model is given by: min et al. (2007). From left to right: (a) ground truth of the motion, (b) motion estimate using l2-regularization [START_REF] Horn | Determining optical flow[END_REF], and (c) motion estimate using TV-regularization [START_REF] Wedel | Statistical and Geometrical Approaches to Visual Motion Analysis[END_REF]. Unlike l2-regularization that over-smooths the solution because its quadratic penalty, TV respects motion discontinuities thanks to its l1-based penalty. (Images courtesy of Mark Werlberger)

TV(u)= X x2Ω |ru(x)| = X x2Ω q u 2 x (x)+u 2 y (x),
u TV S (u)+ η 2 X x2Ω |u(x) -f (x)| 2 , (2.23)

Chambolle's projection algorithm

Problem (2.23) can be solved using the projection algorithm by [START_REF] Chambolle | Total variation minimization and a class of binary mrf models[END_REF], which works on the dual formulation and provides the solution:

u(x)=f (x)+ 1 η r • p(x) (2.24)
where the dual variable p : Ω ! R 2 is computed iteratively as:

p n+1 (x)= p n (x)+(ην)ru n (x) max {1, |p n (x)+(ην)ru n (x)|} (2.25)
with p 0 = 0 and ν < 1/4. However, this approach is only a choice between a large set of alternatives to solve this kind of problems, see [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF].

Arrow-Hurwicz primal-dual algorithm

Algorithms that simultaneously solve the primal and dual formulation converge significantly faster than methods based solely on the primal or dual formulation.

An efficient and simple primal-dual formulation is presented by [START_REF] Zhu | An efficient primal-dual hybrid gradient algorithm for total variation image restoration[END_REF], where the dual variable is iteratively computed as:

p n+1 (x)= p n (x)+(ν 1 )ru n (x) max {1, |p n (x)+(ν 1 )ru n (x)|} (2.26)
and the solution u is iteratively given by:

u n+1 (x)= u n (x) -ν 2 r • p(x)+(ην 2 )f 1+ην 2 (2.27)
with p 0 = 0 and u 0 = 0, and the parameters ν 1 > 0, ν 2 > 0 and ν 1 ν 2 < 1/2.

Staircaising Effect

When total variation is used to regularize, it can suffer from the staircaising problem. This effect is characterized by artificial constant areas in the solution.

It is possible to correct these artifacts by replacing the l 1 norm in the TV term (2.22) by the Huber norm [START_REF] Werlberger | Anisotropic Huber-L1 Optical Flow[END_REF]:

|u| α = ( |u| 2 2α |u|  α |u| -α 2 |u| > α , (2.28)
where α > 0 is a small parameter balancing between the l 1 and the l 2 penalization. Quadratic regularization is applied on small values while TV regularization is used with larger values, thus diminishing the staircase effect. The Huber norm can be easily used with any of the algorithms presented above.

The starcaising problem can also be solved by penalizing the variation of higher derivatives of u, as is propesed by [START_REF] Bredies | Total generalized variation[END_REF]w i t ht h etotal generalized variation (TGV). Particularly, the second order TGV is defined by:

TGV 2 (u)=min v α 0 X x2Ω |ru(x) -v(x)| + α 1 X x2Ω |rv(x)|, (2.29)
where α 0 and α 1 balance between the first and second order derivative terms.

Penalizing the variation of the second order derivative benefits smooth solutions, eliminating the staircase effect.

Vectorial total variation

The idea of total variation can be extended to a vector-valued function u : Ω ! R N . [START_REF] Goldluecke | The natural vectorial total variation which arises from geometric measure theory[END_REF] analyze some of the options to define a vectorial total variation (VTV). Such approaches can be divided in two classes: those that compute the TV channel by channel and those that consider the Riemann geometry of the image manifold. An intuitive way to define the VTV is to sum up the contribution of every channel. The channel-by-channel VTV is defined by: .30) This definition allows a fast and simple regularization separately by channel, where any of the two algorithms shown above can be used. However, since there is no coupling between channels, the correlation between the different components is not exploited. For instance, depending of the nature of u, its discontinuities may appear at the same position in multiples channels and this information can benefit the regularization to preserve the multidimensional edge. This is the case when TV regularization is applied to color images or a motion field, the latter is the main concern of this thesis.

TV S (u)= N X i=1 TV(u i ). ( 2
By considering u as a two-dimensional Riemann manifold, the derivative matrix, Du := (ru 1 , ru 2 ..., ru N ) T : Ω ! R N ⇥2 , can be used to compute the metric tensor g(u)=( D u) T Du. In this case, the eigenvector λ + corresponding to the largest eigenvalue, gives the direction of the vectorial edge. [START_REF] Goldluecke | An approach to vectorial total variation based on geometric measure theory[END_REF] propose the use of the Jacobian J 1 from geometric measure theory, which corresponds to p λ + and can be computed as the largest singular value of matrix Du. Accordingly, the VTV can be defined by:

TV σ 1 (u)= X x2Ω σ 1 (Du(x)) , (2.31)
where σ 1 is the largest singular value of the derivative matrix. The metric TV σ 1 has two main properties: it weighs the contribution of every channel, unlike TV S , and it supports a common edge direction for all channels. This latter property is the main advantage of TV σ 1 with respect other metrics. For example, the Frobenius norm of Du can be used to define the metric: (2.32) allowing a coupling between all channels, but since it does not support a common edge, the vectorial edge is strongly penalized. It can be shown that [START_REF] Goldluecke | An approach to vectorial total variation based on geometric measure theory[END_REF], reflecting the softer penalization performed by TV σ 1 which protects the vectorial edge.

TV F (u)= X x2Ω kDu(x)k ,
TV σ 1 (u)  TV F (u)  TV S (u)

Robust Lucas-Kanade framework

The goal of Lucas-Kanade algorithm is to align a template image, T (x), to an input image, I(x), therefore, it allows to compute optical flow or to track an image patch between two frames. Following [START_REF] Baker | Lucas-Kanade 20 years on: A unifying framework[END_REF], this problem can be stated as finding the parameter vector θ =(θ 1 ,...,θ n ) T which minimizes a squared sum of differences between the template T an the current image I: (2.33) where the warp W(x, θ) maps each template pixel x to a pixel on the image. The warping function can be defined using an image motion model, e.g., one of those presented in Chapter 4. This approach was introduced by [START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF] for image registration, assuming that the optical flow corresponds to a single (continuously varying) motion locally; that is, the variation inside the patch is assumed to be Gaussian. Formulation (2.33) is a least-squares problem, which can be seen as a particular case of maximum a posteriori estimation as is shown below.

θ = arg min θ X x [I (W(x, θ)) -T (x)] 2 ,
Bayesian reasoning can provide an alternative solution.

Bayesian derivation of the Lucas-Kanade framework

Assume that we want to estimate parameters θ = {θ 1 ,...,θ n } given a set of observations ρ = {⇢(x 1 , θ),...,⇢(x L , θ)}, where each observation is defined as the residual between the warped image and the template according to ⇢(x i , θ)= I (W(x i , θ)) -T (x i ). Let P (θ)b et h eap r i o r idistribution of the parameters. Baye's law converts the prior belief about parameters θ,i n t oa posterior probability P (θ|ρ)b yu s i n gt h elikelihood function as follows: .34) thus, the maximum a posteriori (MAP) estimate is defined as:

P (θ|ρ)= P (ρ|θ)P (θ) P (ρ) , ( 2
θ = arg max θ P (θ|ρ). (2.35)
Because of noise, residuals are not zero, but follow a particular distribution P (ρ|θ). Assuming that the noise of all observations is independent and identically distributed, and since P (ρ) does not depend on θ, the MAP estimate becomes:

θ = arg max θ Y x P (⇢(x, θ)|θ)P (θ).
(2.36)

The prior P (θ) models the prior distribution of the parameter. In the Lucas-Kanade framework this belief corresponds to the prior distribution of image motions, which is assumed uniform, so that every possible parameter, or motion, is equiprobable and P (θ) can be dropped from (2.36). Therefore, the MAP estimate is equivalent to minimizing the negative log of the likelihood function:

θ = arg min θ - X x log (P (⇢(x, θ)) . (2.37)
If the noise is Gaussian, e.g., P (⇢(x, θ)|θ) ⇠ N (⇢(x, θ), 0, σ), the MAP estimate applied to the Lucas-Kanade framework, gives: 2.38) where the constant related to σ can be removed.

θ = arg min θ 1 2σ 2 X x ⇢ 2 (x, θ), (

Robust least-squares

If the noise distribution is not Gaussian, the least-squares estimate may be skewed from the real (or desirable) solution. This is a very common case for motion estimation using the Lucas-Kanade framework. The Gaussian assumption is commonly violated due to occlusion or multiple motions, yielding undesired solutions. Particularly, in a least-squares approach, outliers contribute too much to the overall solution, because they are not expected in the prior distribution. In order to improve robustness, the estimator should be able to recover the structure that best fits the majority of the data while identifying and rejecting outliers.

The goal of a robust estimator is to find the values for the parameters, θ, that best fit model W(x i , θ) for the observed measurements ρ, in cases where the data differs statistically from the model assumptions. For this reason the quadratic, or l 2 -, penalty, brought by the Gaussian assumption, is replaced by a function ,whic his called a M -estimator since it is related to the maximum-likelihood estimate, shown in Equation 2.37. Robust least-squares estimation is formulated as:

θ = arg min θ X x ⇢ 2 (x, θ) , (2.39)
The specific M -estimator, (s 2 ), characterizes the insensitivity to outliers or deviations (robustness) of the estimation. The case (s 2 )=s 2 corresponds to the quadratic penalty, where all square residues are treated the same and outliers may skew the solution. If Equation (2.39)i sd i fferentiable, the derivative of the Mestimator weighs every observation in the solution, and can be used to analyse the influence of outliers in the estimate. For instance, the Lorentzian is a M -estimator defined as:

a (s 2 ) = log ✓ 1+ s 2 2a 2 ◆ , with 0 a (s 2 )= 1 s 2 +2a 2 , (2.40)
has a derivative which decreases as s 2 moves away from zero, reducing the influence of outliers. A less robust M -estimator, i.e., with a slower decreasing derivative, is the Charbonnier penalty, is given by: .41) corresponding to a differentiable approximation of the l 1 norm. This penalty has shown to perform better than other M -estimators for optical flow estimation [START_REF] Sun | Secrets of optical flow estimation and their principles[END_REF]. Figure 2.9 shows these three penalties and their derivatives.

ε (s 2 )= p s 2 + " 2 , with 0 ε (s 2 )= 0.5 p s 2 + " 2 , ( 2

Iterative solution

The robust Lukas-Kanade estimate is defined by: (2.42) which is a nonlinear optimization problem, which can be solved by iterative reweighted least squares (IRLS) with the Gauss-Newton algorithm. Assuming that an initial estimate θ is known, the input image I can be approximated around it, by using the first-order Taylor expansion. Accordingly, (2.42) of finding an incremental ∆θ that minimizes

θ = arg min θ X x [I (W(x, p)) -T (x)] 2 ,
X x  ⇢(x, θ)+I x @W @θ ∆θ 2 ! , (2.43)
with ⇢(x, θ)=I (W(x, θ))-T (x) the image residue, and where I x =(@I/@x, @I/@y) denotes the input image gradient and @W/@θ the Jacobian of the warp, at evaluated at W (x, θ). Finding the minimum of expression (2.43) requires an iterative solution. Taking the partial derivative with respect to ∆θ gives: .44) and at the minimum of (2.43), equation (2.44) is zero, so that the incremental scene flow ∆θ can be computed by the following iterative procedure:

2 X x 0 (⇢ 2 (x, θ + ∆θ)) ✓ I x @W @θ ◆ T ✓ ⇢(x, θ)+I x @W @θ ∆θ ◆ , ( 2
1. Initialize: ∆θ = 0.

2. Update weights with the current estimation, ✓, and compute the increment:

∆✓ = -H -1 X x 0 (⇢ 2 (x, θ)) ✓ I x @W @θ ◆ T ⇢(x, θ), (2.45)
where H is the Gauss-Newton approximation of the Hessian matrix:

H = X x 0 (⇢ 2 (x, θ)) ✓ I x @W @θ ◆ T ✓ I x @W @θ ◆ . ( 2

Problem Analysis

Problem definition

The 3D motion field of a scene, or scene flow, is useful for several computer vision applications, such as action recognition, interaction, or 3D modeling of nonrigid objects. Scene flow methods can be classified into two main groups depending on if the 3D structure of the scene is estimated or provided. The first group of methods uses a stereo or multi-view camera system, where both scene flow and depth are estimated from images. The second group uses RGBD images as input. In this case, the depth given by the sensor may be used directly for scene flow estimation. This thesis proposes a new formulation for this second group.

The two main questions for scene flow estimation from RGBD images are: (a) how to fully exploit both sources of data, and (b) which motion model should be used to compute a confident scene flow. Regarding the first question, 2D motion estimation has made a great progress in the last years and powerful tools are available to estimate optical flow from color images [START_REF] Sun | Secrets of optical flow estimation and their principles[END_REF]. Optical flow is defined as the apparent motion of brightness pattern in the image and is closely related to image flow, which is the projection of the 3D motion field on the image. However, optical and image flow may not coincide, as is show in Figure 3.1. Moreover, changes on the motion field are not always reflected in the optical flow or can be hard to estimate, such as the motion of untextured surfaces. Also, changes in depth and occlusions are a challenge for current optical flow methods since there is no 3D knowledge to disambiguate the possible 2D motions. Thanks to the RGBD sensor, information about the 3D structure of the scene is available. By using depth data alone it is possible to estimate a global rigid motion, e.g., the camera motion, as is proposed by [START_REF] Newcombe | Kinectfusion: Real-time dense surface mapping and tracking[END_REF]u s i n gIterative Point Cloud (ICP) [START_REF] Besl | A method for registration of 3-d shapes[END_REF]. However, solving for a nonrigid motion field is a more challenging problem having a wider scope. Although it is possible to perform nonrigid registration of depth data [START_REF] Li | Robust single-view geometry and motion reconstruction[END_REF], low-noise and high-resolution depth data is required, which prevents the use of current RGBD sensors. The availability of RGBD data opens the door to a new scenario for motion estimation, where some of these drawbacks can be overcome.

When both color and depth are provided, there is no consensus on how to most effectively combine these two sources of information. A straightforward approach would be to compute the optical flow from RGB images and to infer the scene flow from the depth data. Alternatively, the texture of depth images can be used to help the 2D motion estimation and simultaneously used to compute the changes in depth. It is also possible to generate colored 3D point clouds or 3D patches by exploiting the depth data and try to estimate a 3D transformation consistent with the input images. Nevertheless, it is not necessary to explicitly represent points in 3D: the scene structure can be represented in the image domain, where color and depth data can be coupled using a projective function. This way, depth changes influence the motion in the image domain and consistency constraints can be formulated jointly over the color and depth images. We survey some of these alternatives by reviewing the related work to this thesis in the next Section.

The second question, concerning the representation of the motion, has been addressed less. Scene flow estimation using intensity and depth is an ill-posed problem and regularization is needed. Instead of solving for 2D motions, as in optical flow estimation, the 3D motion vector of each point can be directly solved to minimize color and intensity constraints in a data term. Moreover, the whole 3D motion field can be regularized to get spatially smooth solutions while preserving discontinuities. Since depth data is available, it can be used in the regularization to preserve motion discontinuities along depth edges, where independent motions are probable to appear, and also can be exploited to define a prior of the local rigidity. However, solving for a piecewise smooth solution of the 3D motion field may not be the best choice for some motions of interest. For example, a 3D rotation of a rigid surface induces a great variety of 3D motions that are hardly well regularized by such an approach. A similar issue occurs when the RGBD sensor is moving with respect to the scene. For this reason the motion representation should be able to model rigid and nonrigid components of the scene flow, while favoring the regularization procedure. Moreover, scene flow estimation can be formulated in 3D space or in the image domain and is not clear in which domain the parametrization of the scene yields more accurate results. While 3D space is the natural domain of the scene flow, RGBD images are sensed and provided in the image domain and this 2D parametrization can be exploited. This thesis aims to make a step forward towards a better understanding of scene flow estimation from RGBD images. For this reason, in order to cover the two questions previously formulated, this study focuses in four main aspects: i )t h e coupling of intensity and depth, ii ) the representation of motion, iii ) the constraints of motion in RGBD data and iv ) the regularization. Temporal consistency of the scene flow is out of the scope of this thesis and the estimation problem is formulated on a frame-to-frame scenario.

Related work

Scene flow was first introduced by [START_REF] Vedula | Three-dimensional scene flow[END_REF] as the 3D motion field of the scene. Since this seminal work several approaches have been proposed to compute the 3D motion field. If a stereo or multi-view (SoMV) camera system is available, scene flow can be computed by enforcing consistency with the observed optical flows [START_REF] Vedula | Threedimensional scene flow[END_REF], by a decoupled [START_REF] Wedel | Efficient dense scene flow from sparse of dense stereo data[END_REF] or joint [START_REF] Huguet | A variational method for scene flow estimation from stereo sequences[END_REF][START_REF] Basha | Multi-view scene flow estimation: A view centered variational approach[END_REF] estimation of structure and motion, or by assuming a local rigidity of the scene [START_REF] Vogel | 3D scene flow estimation with a rigid motion prior[END_REF]. On the other hand, if a depth sensor is available the estimation of the structure is not needed and color and depth information can be simultaneously used. In this case scene flow can be computed by considering the depth data as an additional image channel [START_REF] Spies | Dense range flow from depth and intensity data[END_REF], by coupling color and depth in the image domain with a projective function [START_REF] Quiroga | Scene flow by tracking in intensity and depth data[END_REF][START_REF] Quiroga | Local scene flow by tracking in intensity and depth[END_REF] or by reasoning directly in 3D space [START_REF] Hadfield | Scene particles: Unregularized particle-based scene flow estimation[END_REF][START_REF] Hornacek | Sphereflow: 6 dof scene flow from RGB-D pairs[END_REF]. Although in this thesis we assume that a depth sensor is available and the estimation of the structure is not needed, below we summarize the most relevant previous methods based on SoMV systems and RGBD images. Scene flow estimations from SoMV and RGBD images are very related, as well as with optical flow estimation, so that their jointly survey could help to a better understanding of this motion estimation problem. A more detailed survey is done for RGBD-based methods in order to clarify the main contributions of this thesis.

Scene flow from stereo or multi-view

Most scene flow methods assume a stereo, or multi-view camera system, in which the motion and the geometry of the scene are jointly estimated, in some cases, under a known scene structure. Inspired by optical flow methods, the most common approach for estimating scene flow is to perform an optimization on a global energy, including photometric constraints and some regularization. Since optical flow is an approximation of the projection of the 3D motion field on the image plane, an intuitive way to compute scene flow is to reconstruct it from the optical flow measured in a multi-view camera system, as proposed by [START_REF] Vedula | Three-dimensional scene flow[END_REF], or including a simultaneously structure estimation as by [START_REF] Zhang | On 3D scene flow and structure estimation[END_REF]. Although it is possible to estimate optical dealing with large displacements [START_REF] Brox | Large displacement optical flow: Descriptor matching in variational motion estimation[END_REF] and preserving small details [START_REF] Xu | Motion detail preserving optical flow estimation[END_REF], it is difficult to recover a scene flow compatible with several observed optical flows that may be contradictory. Moreover occlusion handling is a challenge since it could affect both the per-view optical flow computation and the posterior scene flow estimation using multiple views. Some authors introduce constraints from a fully calibrated stereo or multiview camera system. [START_REF] Huguet | A variational method for scene flow estimation from stereo sequences[END_REF] simultaneously compute the optical flow field and two disparities maps by extending the variational optical flow framework by [START_REF] Brox | High accuracy optical flow estimation based on a theory for warping[END_REF] to a stereo system: optical flow constraints for each camera (left, right) and stereo constraints for each time (t, t + 1). This way a temporal consistency is enforced on the disparity map estimation between t and t+1. Nevertheless, this approach may be impractical due to its high computational time. [START_REF] Wedel | Efficient dense scene flow from sparse of dense stereo data[END_REF] decouple disparity and motion estimation to achieve real-time computation, since each problem can be separately solved more efficiently. Although the spatio-temporal information is not fully exploited, the optical flow is enforced to be consistent with the computed disparities. [START_REF] Valgaerts | Joint estimation of motion, structure and geometry from stereo sequences[END_REF] generalize the stereo-based scene flow formulation by embedding epipolar constraints in the energy functional, showing that scene flow and the stereo structure can be simultaneously estimated. Since motion and depth discontinuities do not necessarily coincide, regularization is performed separately. All these methods suffer from the smoothness constraints brought by 2D parametrization, especially due to the regularization on the optical flow rather than directly on the scene flow. [START_REF] Basha | Multi-view scene flow estimation: A view centered variational approach[END_REF] improve estimation by formulating the problem as a point cloud in 3D space supporting multiple views. In this method scene flow and the estimated disparity are regularized independently using total variation (TV). In a closely related approach, [START_REF] Vogel | 3D scene flow estimation with a rigid motion prior[END_REF] regularize the scene flow by encouraging a locally rigid 3D motion field, outperforming TV regularization. In a later work, [START_REF] Vogel | Piecewise rigid scene flow[END_REF] present a more constrained assumption and represent the 3D scene by a collection of rigidly moving planes. The scene flow is computed by iteratively alternating three procedures: plane fitting, pixel-to-plane assignation and motion estimation for every plane. This methods achieves state-of-the-art results in the KITTI stereo flow benchmark [START_REF] Geiger | Are we ready for autonomous driving? the kitti vision benchmark suite[END_REF].

Alternatively, other methods formulate scene flow estimation as a 3D tracking problem. On one hand, it is possible to use an explicit model of the scene and compute a 3D motion consistent with the pre-frame reconstruction and the photometric information. [START_REF] Neumann | Spatio-temporal stereo using multi-resolution division surfaces[END_REF] use a subdivision surface model, where the shape and motion of an object are initialized independently, frame to frame, and then refined simultaneously. [START_REF] Fukurawa | Dense 3D motion capture from synchronized video streams[END_REF] track an initial 3D mesh model of the scene by a local estimation and a posterior full deformation of the mesh. The motion of each vertex is estimated using photometric information across different views, afterwards the full mesh is deformed enforcing smoothness and local rigidity. Another possibility is to track a sparse set of 3D points or surface elements [START_REF] Carceroni | Multi-view scene captured by surfer sampling: From video streams to non-rigid 3D motion[END_REF][START_REF] Devernay | Multi-camera scene flow by tracking 3D points and surfels[END_REF]. [START_REF] Carceroni | Multi-view scene captured by surfer sampling: From video streams to non-rigid 3D motion[END_REF] model the scene as a set of planar surfels (surface elements) and solve for their pose and motion using photometric constraints across different views. This approach requires a well-controlled lighting and acquisition setup, and because its complexity, scene flow solution is only suitable in a limited volume. [START_REF] Devernay | Multi-camera scene flow by tracking 3D points and surfels[END_REF] directly get a set of 3D trajectories for 3D points from which the scene flow is derived. However, this method suffers from drift problems and its proposed method for point visibility handling is difficult to perform.

Scene flow from RGBD images

The first work using intensity and depth was by [START_REF] Spies | Dense range flow from depth and intensity data[END_REF], where the optical flow formulation by [START_REF] Horn | Determining optical flow[END_REF]i se x t e n d e dt oi n c l u d ed e p t h data. In this approach, depth data is used simply as an additional channel in the variational formulation of the optical flow trough the definition of a range flow equation encouraging the consistency between 2D motion and changes in depth. This way both intensity and depth image gradients are combined to jointly estimate the optical and the range flow, which are enforced to be smooth. [START_REF] Lukins | Colour constrained 4d flow[END_REF] extend this approach to multiple color channels and one aligned depth image. Nevertheless, in both methods, the scene is assumed to be captured by an orthographic camera and there is no an optimal coupling between optical and range flows. Moreover, due to the Horn-Schunck framework, these methods suffer from the early linearization of the constancy constraints and from over-smoothing along motion boundaries because of the L 2 -regularization. The coupling issue can be solved by using a projective function that allows to constrain the 3D motion in the image domain. Because depth information is available it is possible to define a projective function for every pixel using its depth measure. This way changes in depth are used to support the constraint of motion in the image domain. Following this idea, [START_REF] Letouzey | Scene flow from depth and color images[END_REF]p r o j e c tt h es c e n efl o w in the image domain using a projective matrix to relate 3D motion and image flow. Thus scene flow can be directly solved by using photometric consistency constraints. However, this method uses a projective function that only supports smalls changes in depth, also it suffers for the well known drawbacks of the Horn-Schunck framework. Furthermore the depth data is not fully exploited since there is no range flow constraint and is only used to support the projective function.

Rather than computing a dense scene flow, we define a 2D warping function to couple image motion and 3D motion [START_REF] Quiroga | Scene flow by tracking in intensity and depth data[END_REF], allowing for a joint local constraint of scene flow on intensity and depth data. The local 3D motion is modeled as a 3D translation and solved as a tracking problem in RGBD data, using a Lucas-Kanade based framework [START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF][START_REF] Baker | Lucas-Kanade 20 years on: A unifying framework[END_REF]. Unlike the projective matrix [START_REF] Letouzey | Scene flow from depth and color images[END_REF], this warping function allows to project the true 3D motion on the image and can be approximated by a linear model, to iteratively solve for the scene flow. Although the method is able to deal with large displacements, it fails on untextured regions and more complex motions, such as rotations. In order to solve for dense scene flow, a regularization procedure that preserves motion discontinuities is required. Inspired by optical flow methods, the 3D motion field can be assumed to be piecewise smooth, and total variation (TV) can be applied as regularizer.

The work by [START_REF] Herbst | RGB-D flow: Dense 3-D motion estimation using color and depth[END_REF] follows this idea, but as [START_REF] Spies | Dense range flow from depth and intensity data[END_REF], it lacks a coupling between optical and range flows, and the regularization is done on the optical flow rather than on the scene flow. Such that approach lacks of the accuracy to recover a confident scene flow, and it suffers of the same drawbacks that optical flow methods. The regularization of 3D motion field is more reliable than encouraging smoothness on the apparent motion in the image domain. In [START_REF] Quiroga | Local/global scene flow estimation[END_REF] we present a variational extension of the local method for scene flow [START_REF] Quiroga | Scene flow by tracking in intensity and depth data[END_REF]. In this work, the tracking approach in intensity and depth is regularized to compute a dense scene flow enforcing piecewise smoothness. In order to preserve motion discontinuities along depth edges, a weighted TV is applied on each component of the 3D motion field.

All these methods assume spatial smoothness of the scene flow, which is a reasonable assumption for translational motions but not the optimal choice for rotations. Under a rotation, even close scene points present different 3D motions. Particularly, in the case of a moving camera the regularization of the motion field can be a challenge. For this reason, we present an over-parametrization of the scene flow (Quiroga et al., 2014a), where each scene point is allowed to follow a rigid body motion. This way, the regularization can be done on a field of rigid motions, favoring piecewise solutions, which is a better choice for real scenes. A similar idea is presented by [START_REF] Rosman | Group-valued regularization for analysis of articulated motion[END_REF], where a regularization of a field of rigid motions is proposed. However, our work (Quiroga et al., 2014a)d i ffers from such that by [START_REF] Rosman | Group-valued regularization for analysis of articulated motion[END_REF], in three ways. First, a more compact representation is used of the rigid-body motion via the 6-parameter twist representation, instead of a R 12 embedding. Second, this approach solves and regularizes the rigid motion field at the same time. Finally, it is proposed a decoupled regularization of the rotational and translational fields, which simplifies the optimization and allows the use of different TV strategies on each field. Also, as we present in [START_REF] Quiroga | Local/global scene flow estimation[END_REF], a depth-based weighting function is used to avoid penalization along surface discontinuities and a warping function is defined to couple the twist motion and the optical flow. A very similar warp is presented by [START_REF] Kerl | Robust odometry estimation for RGB-D cameras[END_REF] to solve for a global rigid motion from RGBD images. However, we use the warping function to locally constrain the rigid motion field in the image domain, while the global rigid motion estimation can be seen as a particular case of this formulation. Moreover, unlike [START_REF] Kerl | Robust odometry estimation for RGB-D cameras[END_REF], a depth consistency constraint is defined to fully exploit both sources of data.

An alternative approach is the representation of the scene in 3D space. [START_REF] Hadfield | Scene particles: Unregularized particle-based scene flow estimation[END_REF] estimate the scene flow using particle filtering in a 3D colored point cloud representation of the scene. In that approach, a large set of 3D motion hypotheses must be generated and tested for each 3D point, leading to high computational costs. Also, although the propagation of multiple hypotheses could avoid to get stuck at local minima, it is not clear the benefit brought by formulating the problem on an unstructured 3D point cloud. [START_REF] Hornacek | Sphereflow: 6 dof scene flow from RGB-D pairs[END_REF] compute scene flow by matching rigid patches of 3D points identified as inliers of spheres, where RGB and depth information are simultaneously used to define a matching cost. Rather than solve for 3D motions, this approach estimates an optimal rigid motion for every patch, by generating and propagating a set of hypotheses. The scene flow is Iteratively regularized by enforcing local rigidity, which is done as a labeling optimization problem. This approach is very related to our work (Quiroga et al., 2014a), in the sense that both solve for a dense field of rigid motions and encourage local rigidity. However, these two approaches differ in several aspects, such as the kind of RGBD consistency constraints, and the optimization strategy: while [START_REF] Hornacek | Sphereflow: 6 dof scene flow from RGB-D pairs[END_REF] use discrete optimization via a modified PatchMatch [START_REF] Hornacek | Depth super resolution by rigid body self-similarity in 3d[END_REF], Quiroga et al. (2014a) use a Gauss-Newton algorithm to perform continuos optimization. Nevertheless, the main difference concerns our scene flow framework that minimizes a local/global energy, which enables an adjustable combination between local and piecewise rigidity.

Proposed approach

In this thesis, we take advantage of the fact that most real scenes can be well modeled as locally or piecewise rigid, i.e., the scene is composed of 3D independently rigid components. We define a semi-rigid scene as one that its 3D motion field is well approximated as locally or piece-wise rigid, or a combination of both, see Figure 3.2. This way scene flow estimation can benefit from the rigid properties of the scene while allowing motion discontinuities and outliers, being able to model complex motion patterns. In order to exploit this assumption, we define a framework for scene flow estimation from RGBD images, using an over-parametrization of the motion. We represent the scene flow as a vector field of rigid body motions, which is expressed as a field of twist motions and solved directly in the image domain. This representation helps the regularization process, which, instead of directly penalizing variations in the 3D motion field, encourages piecewise smooth solutions of rigid motions by using a weighted TV regularization to preserve motion discontinuities along surface edges. Moreover, this semi-rigid approach can be constrained in the image domain to fully exploit intensity and depth data. For this purpose, we define a projective function that makes it possible to constrain every rigid motion on the image and to define an incremental solver, which is embedded in a robust least-square estimation to exploit local rigidity. Thus, we can solve for the local twist motion that best explains the observed intensity and depth data, gaining robustness under noise. The local solver, in conjunction with the TV regularization of the twist field, provides an adjustable combination between local and piecewise rigidity. This formulation is flexible enough to support different data constraints and regularization strategies, and it can be adapted to more specialized problems. By using the same general framework, it is possible to model the scene flow as a global rigid motion plus a nonrigid residual, which is particularly useful when estimating the motion of deformable objects in conjunction with a moving camera.

As we present in the experiments, this framework is flexible enough to estimate a global rigid body motion or to solve for a general 3D motion field in challenging setups. This way, we are able to model the motion of the scene as a global rigid motion and a nonrigid residual.

Scene Motion Representation

Real world scenes are composed of a set of 3D surfaces, which in the more general case, are free to move and deform independently. Using a depth sensor we are provided with a finite set of samples of the scene. This is assumed to be enough to represent the most important details of the 3D structure and to observe the motions of interest in the scene. In the RGBD representation, every valid sample of the scene, or scene point, is characterized by a RGB color and a depth measure.

If between two time steps there exists a motion of any 3D surface or of the RGBD camera, or both simultaneously, the RGBD representation of the scene changes, enabling the estimation of the 3D motion field.

Perhaps the most intuitive approach to solve for the scene flow is to use the RGBD data without imposing any rule about the 3D motion. If any scene point moves freely there are no limitations for the motions that can be modeled. However, this independence between scene points, which may favor the richness of motions modeled, leads to a very ill posed problem for two main reasons. First, the color of every scene element is not discriminative since several scene points may have the same color. Moreover, because of the presence of noise and illumination changes, color-based consistency assumptions are commonly violated. Second, using each depth measure alone does not make it possible to constrain 3D motion since depth measures change with the motions in the scene. Therefore, any scene point with a similar color would be a potential candidate for the new position of the scene point and the 3D structure information cannot be exploited.

An appropriate motion model is required to reliably estimate the 3D motion of every scene point . This model should be general enough to cover a wide variety of interesting 3D motions but also it should allows constraints on the set of motions, making it possible to state the scene flow formulation as a well-posed problem. Most real world scenes are composed of 3D surfaces that can be well approximated as locally or piece-wisely rigid, even very complex scenes. We define a semi-rigid surface as a surface that any of its feasible motions can be well approximated as a locally or piece-wise rigid motion, or a combination of both. This assumption is fundamental because allows the use of a smoothness prior on the motion of scene points. Note that the assumption of semi-rigidity does not imply planar surfaces. On a locally rigid surface, the 3D distance between scene points remains constant and the 3D motions of these points are closely related, as well as the motion of their projections on the image. Usually this assumption is exploited by assuming local solutions of the optical flow or the 3D motion field. However, points belonging to a rigid surface only have the same 3D motion if the rotational component is zero. Otherwise, every point performs a different 3D displacement, as is observed from Equation (2.17). Although that kind of regularity yields a better statement of the problem, it does not totally exploit the properties of a locally rigid scene. Similarly, the properties of a piece-wise rigid scene explain the benefits of using a TV-based regularization, as is done in optical flow methods, but it is the chosen representation of motion which determines the degree of utilization of the semi-rigid assumption in the motion solution.

In this thesis we go beyond previous approaches by considering a semi-rigid scene. For this purpose, scene flow is represented by a field of twist motions, where local and piece-wise rigidity can be fully exploited. This over-parametrization makes it possible to model scenes containing rigid and nonrigid motions, and as is presented in Section 4.4, it includes and generalizes other motion representations.

Twist-motion representation

Every 3D point X, or scene point, is assumed to belong to a rigid surface S,w h i c h is free to move between two time steps. Let ⇠ be the 6-vector twist representing the rigid motion of the surface S. Under the twist action, every surface point performs a 3D motion and the assemblage of motions of all surfaces generates the observed scene flow. Instead of directly representing the elements of the motion field as 3D vectors, we use an over-parametrized model to describe the motion of every point as the rigid transformation of the surface. The 3D motion field of the scene is embedded in a twist motion field, ⇠(x):Ω ! R 6 , where the motion of a 3D point X,b e t w e e nt w ot i m es t e p s ,i sm o d e l e db y :

Xt+1 = e ξ(xt)
Xt .

(4.1)

In the twist-motion representation, every scene point follows a rigid motion, which is represented by 6 parameters. This way is possible to encourage piecewise smooth or constant solutions on this field of twist motions. This representation exploits the local rigidity of the scene in order to better approximate complex 3D motions. Nevertheless, this representation does not expect any a priori semi-rigid configuration enabling a flexible formulation of the scene flow estimation, as is presented in Chapter 5.

Between any two time steps, the scene flow, v(x), can be directly recovered from a given twist motion field, ⇠(x), as follows:

v(x) 1 ! = ⇣ e ξ(x) -I 4⇥4 ⌘ X . (4.2)

Twist motion on the image

Under the action of a twist ⇠ between t and t + 1, the motion of every 3D point induces an image flow (u, v) in its RGBD representation. Changes and invariances of the RGBD data can be used to solve for the scene flow, as is described in Section 5.2. For this purpose, it is useful to define a warping function that allows the constraint of a rigid motion in the image domain, enabling the exploitation of color and depth images. Particularly, because the twist representation, this warp is required to map the twist action to the image domain. Let W ξ (x, ⇠):R 2 ⇥R 6 ! R 2 be the warping function that maps each (non-occluded) scene point to its new location after the rigid motion ⇠. The warping function is defined as:

W ξ (x, ⇠)=⇡ ⇣ e ξ X⌘ , with X = ⇡ -1 (x,Z(x)) 1 ! . (4.3)
Accordingly, under the action of twist motion ⇠, a 3D point X t becomes X t+1 , and its projection on the image is modeled by:

x t+1 = W ξ (x t , ⇠). (4.4)
Let S k (x) be a component of the RGBD image, e.g., brightness or depth, evaluated at x 2 Ω, and let ⇢ S k (x, ⇠)b ead i fferentiable consistency function between S k t (x) and S k t+1 (W ξ (x, ⇠)). Every twist motion vector ⇠ is solved aiming to minimize one or multiple consistency functions like ⇢ S k (x, ⇠). Optimization problems on manifolds, such as SE(3), can be solved by calculating incremental steps in the tangent space to the manifold. Particularly, the tangent space of the Lie group SE(3) corresponds to its Lie algebra se (3). If a gradient-based procedure is used, it requires computation of the derivatives of the warp with respect to the twist motion ⇠,w h e r e ⇠ 2 se(3), see Equation (2.8). The warp (4.3) is nonlinear in ⇠ yielding to a complicated expression for its Jacobian. However, this can be well approximated by a linear approximation around any small neighborhood on the manifold, as is presented below. The definition of this linearized approximation simplifies the optimization procedure providing a simple expression of the Jacobian, which benefits the efficiency of any solver. Moreover, the linearization makes it possible to analyze how each component of the twist motion affects the motion in the image domain.

Linearization of the warping function

Let x be the projection on the image of a given 3D point X. Under the action of twist ⇠,t h ep i x e lx becomes x ξ = W ξ (x, ⇠). Moving away from ⇠, on the manifold, induces an additional image flow on the projection of the 3D point. Let ∆⇠ be a small twist increment, i.e, close to the identity in SE(3), applied from ⇠,yieldingto the equivalent twist motion ⇠ ∆ = log(e ∆ξ e ξ ). The projection of 3D point X,u n d e r the action of ⇠ ∆ , can be expressed using the warping function as follows:

W ξ (x, log(e ∆ ξ e ξ )) = W ξ (x ξ , ∆⇠)=x ξ + δx(x ξ , ∆⇠), (4.5)
where δx(x ξ , ∆⇠) is the image flow induced by the twist increment ∆⇠.K n o w i n g that ∆⇠ is close to the identity, the twist induces a small rotation and so that the exponential function can be well approximated by the linear expression e ∆ξ ⇡ I + ∆⇠, see Section 2.2.2. Therefore, the image flow can be written as follows:

δx(x ξ , ∆⇠)=⇡ ⇣ e ∆ξ Xξ ⌘ -⇡ ⇣ Xξ ⌘ ⇡ ⇡ ⇣ (I + ∆⇠) Xξ ⌘ -⇡ ⇣ Xξ ⌘ , (4.6)
with

X ξ = ⇡ -1 (x ξ ,Z(x ξ ))
, and where we abuse of notation of ⇡ to project 3D points directly in homogeneous coordinates. Applying the twist and the projective function in (4.6), we have:

δx(x ξ , ∆⇠) ⇡ 0 @ f x X ξ +∆X Z ξ +∆Z + c x f y Y ξ +∆Y Z ξ +∆Z + c y 1 A - 0 @ f x X ξ Z ξ + c x f y Y ξ Z ξ + c y 1 A , (4.7)
with (∆X, ∆Y, ∆Z) T the vector that contains the components of the 3D displacement induced by the twist ∆⇠, which are given by:

0 B B @ ∆X ∆Y ∆Z 1 C C A = 0 B B @ 0 -Z ξ Y ξ Z ξ 0 -X ξ -Y ξ X ξ 0 1 C C A 0 B B @ ∆! X ∆! Y ∆! Z 1 C C A + 0 B B @ ∆⌧ X ∆⌧ Y ∆⌧ Z 1 C C A . (4.8)
Equation ( 4.9) can be simplified if the Z-component of the 3D motion is very small with respect to the depth of the 3D point. By assuming that |∆Z/Z| ⌧ 1 then

Z ξ (Z ξ + ∆Z) ⇡ Z ξ 2
, and the induced flow becomes:

δx(x ξ , ∆⇠) ⇡ 0 B @ f x X ξ ∆Z+Z ξ ∆X Z 2 ξ f y Y ξ ∆Z+Z ξ ∆Y Z 2 ξ 1 C A . (4.9)
Therefore, by defining J(x ξ ) as the Jacobian matrix evaluated at x ξ , as follows:

J(x ξ )= 0 B @ fx Z ξ 0 - x ξ Z ξ - x ξ y ξ fy fx+x 2 ξ fy - y ξ fx fy 0 fy Z ξ - y ξ Z ξ - fy+y 2 ξ fx x ξ y ξ fx x ξ fy fx 1 C A , (4.10)
and given a small twist increment ∆⇠, the warping function can be well approximated around the twist ⇠, by the following linear version:

W ξ (x, log(e ∆ξ e ξ )) = W ξ (x, ⇠)+J(x ξ )∆⇠. (4.11)

3D-motion representation

The scene flow can be represented directly as a 3D motion field, which is its original definition in Section 2.1.2. In this representation, every 3D point X is also assumed to belong to a moving rigid surface S, but instead of an explicit representation of the rigid motion, the 3D displacement of every surface points is considered. Let x be the projection of X. Under the motion of the surface, the 3D displacement of X can be written as

v(x)=(v X ,v Y ,v Z ).
The assemblage of all motions forms the scene flow, v(x):Ω ! R 3 , where the motion of every 3D point, between two time steps, is modeled by:

X t+1 = X t + v(x t ) . (4.12)
In the 3D-motion representation, every scene point performs a 3D displacement that is directly represented by 3 parameters. Under this representation, the local rigidity assumption cannot be directly exploited. However it is still possible to encourage smooth piecewise solutions, favoring semi-rigid motions. Moreover, it is also possible to easily fit a local rigid motion, which can be useful to encourage semi-rigid solutions, as is presented in Section 5.3.2.

3D motion on the image

The 3D motion of a scene point induces an image flow,( u, v), in its projection on the image. The definition of a warping function makes it possible to constrain every 3D motion vector on the RGBD images. In order to define the 3D-motion warp, we use the twist warping function (4.3), by considering the 3D displacement as the particular case of the twist motion, where the rotational component is zero. Let W v (x, v):R 2 ⇥ R 3 ! R 2 be the warping function that maps every scene point into its new location, after the action of 3D motion v, as follows:

W v (x, v)=⇡ (X + v) , with X = ⇡ -1 (x,Z(x)).
(4.13) Therefore, if between two time steps, the point X t performs a 3D motion v, becoming X t+1 , its projection on the image is modeled by:

x t+1 = x t + W v (x t , v). (4.14)

Linearization

The warping function (4.13) can be linearized around a given 3D motion v. Under 3D motion v, scene point x becomes

x v = W v (x, v). Let ∆v =(∆v X , ∆v Y , ∆v Z ) T
be a small 3D motion increment, then the projection of the 3D point X, after the motion given by v + ∆v, can be expressed using the warping function as follows:

W v (x, v + ∆v)=x v + δx(x v , ∆v), (4.15)
with δx(x v , ∆v) the image flow induced by ∆v. Knowing that v can be seen as a twist with a zero rotational component, we use the result (4.10), from Section 4.1.2. By defining J(x v ) as the Jacobian matrix evaluated at x v as follows:

J(x v )= 1 Z v 0 @ fx 0 -x v 0 fy -y v 1 A , (4.16)
and given a small increment ∆v, the warping function can be well approximated around v, by the following linear version:

W v (x, v + ∆v)=W v (x, v)+J(x v )∆v. ( 4
.17)

Other motion representations 4.3.1 Rigid body model

The inter-frame image flow induced by the 3D motion of a rigid surface has been modeled in the literature using different assumptions. For brevity, let the focal lengths, {f x ,f y }, be the unit and the optical center, (c x ,c y ), be at the image origin.

The instantaneous velocity of a retinal point (x, y), corresponding to a given 3D point with known depth Z,i sm od e l e db yLonguet-Higgins and Prazdny (1980) as:

0 @ u v 1 A = 1 Z 0 @ 10-x 01-y 1 A 0 B B @ t X t Y t Z 1 C C A - 0 @ xy -(1 + x 2 ) y 1+y 2 -xy -x 1 A 0 B B @ ! X ! Y ! Z 1 C C A , (4.18)
with (u, v)=( @x/@t, @y/@t), and where t

=( t X ,t Y ,t Z ) T and ! =( ! X , ! Y , ! Z ) T
are the instantaneous translation and angular velocity vectors, respectively. [START_REF] Horn | Direct methods for recovering motion[END_REF] use this instantaneous velocity model to describe the inter-frame pixel motion, where velocities become displacements between t and t + 1, and the optical flow is defined as (u, v)=( x t+1x t ,y t+1y t ). [START_REF] Adiv | Determining three-dimensional motion and structure from optical flow generated by several moving objects[END_REF] shows that equation (4.18) approximates the image flow induced by a rigid motion under the assumptions that the view of field of the camera is not very large, and V Z /Z ratio and rotational velocities are very small. This model is named rigid body model in the hierarchy classification of motions proposed by [START_REF] Bergen | Hierarchical model-based motion estimation[END_REF], and it is characterized by the 6 degrees of freedom (DOF) of the rigid motion. Also, this model coincides with the linearized warping function of the twist motion (4.11) and is only valid to capture the image motion for rigid motions close to the identity.

Affine flow model

When the distance between the background surfaces and the camera is large, it is possible to approximate the image motion as an affine transformation [START_REF] Bergen | Hierarchical model-based motion estimation[END_REF] using the following model:

0 @ u v 1 A = 0 @ a 1 + a 2 x + a 3 y a 4 + a 5 x + a 6 y 1 A . (4.19)
This affine flow model can be seen as the particular case of the rigid body model with only 4 DOF, where the rotational components ! x and ! y are assumed to be negligible. Following equation ( 4.18) we can write the image flow as:

0 @ u v 1 A = 1 Z 0 @ 10-x 01-y 1 A 0 B B @ t X t Y t Z 1 C C A + 0 @ -y x 1 A ! Z . (4.20)
This way the induced flow is modeled by the coefficients {a 1 ,...,a 6 }, which are functions of the 3D displacement (t X ,t Y ,t Z ) and the rotational component ! Z .

Planar surface flow model

Equation ( 4.18) assumes that the depth of every pixel is known, which is not always the case for many applications. If the optical flow corresponds to the projection of a planar surface,

k X X + k Y Y + k Z Z = 1, the plane parameters (k X ,k Y ,k Z ) constrain
the flow induced by the 3D motion of the surface. The depth variable in (4.18) can be removed by using the plane equation, thus the image flow is expressed as a function of 8 parameters [START_REF] Adiv | Determining three-dimensional motion and structure from optical flow generated by several moving objects[END_REF]:

0 @ u v 1 A = 0 @ a 1 + a 2 x + a 3 y + a 7 x 2 + a 8 xy a 4 + a 5 x + a 6 y + a 7 xy + a 8 y 2 1 A , (4.21)
where coefficients {a 1 ,...,a 8 } depend on the rigid motion parameters (t X ,t Y ,t Z ) T and (! X , ! Y , ! Z ) T and the surface parameters (k X ,k Y ,k Z ) T . This 8 DOF optical flow model is called planar surface flow model by [START_REF] Bergen | Hierarchical model-based motion estimation[END_REF]. Also, this model is known as the homography model by [START_REF] Baker | Lucas-Kanade 20 years on: A unifying framework[END_REF].

Orthographic camera models

Rigid body and affine models depend on depth data to compute the induced motion. Therefore they make it possible to use depth data jointly with RGB information to constrain motion in the image domain. However, it is also possible to consider a simpler camera model, as an orthographic camera, which is the most widely used assumption in optical flow methods. In this case, the scene flow is inferred from the computed optical flow and the provided depth data. Considering an orthographic camera model allows to assess the contribution of incorporating depth data in the 2D motion constraint.

The action of a rigid transformation G = {R, t} can be written as X t+1 = RX t + t, so that the orthographic projection of the resulting 3D point X t+1 is given by:

⇡ o (X t+1 )= 0 @ R 11 X t + R 12 Y t + R 13 Z t R 21 X t + R 22 Y t + R 23 Z t 1 A + 0 @ t X t Y 1 A , (4.22)
where the orthographic projection

⇡ o : R 3 ! R 2 maps a 3D point X =( X, Y, Z) T to a pixel x =(x, y) T on Ω,b y : ⇡ o (X)=(k x X + c x ,k y Y + c y ) T , (4.23)
with k x and k y parameters of the camera, and (c x ,c y ) its principal point. For the sake of simplicity, we assume that k x = k y = 1 and c x = c y = 0. Therefore, the induced image flow, u =(u, v) T ,b e t w e e nt and t + 1, is given by:

u = ⇡ o (X t+1 ) -⇡ o (X t ), (4.24)
and can be represented by the following model:

0 @ u v 1 A = 0 @ (R 11 -1) x + R 12 y +(R 13 Z + t X ) R 21 x +(R 22 -1) y +(R 23 Z + t Y ) 1 A , (4.25)
which corresponds to the affine motion model. In this case, the depth information Z is irrelevant and is combined with the independent coefficients a 1 and a 4 of the affine model, see equations (4.19) and (4.20). This 4 DOF model is used by [START_REF] Nir | Over-parameterized variational optical flow[END_REF] to compute an over-parametrized optical flow. A simpler model is defined by assuming a negligible rotation with an orthographic camera model. In this simplified model the image flow is given by:

0 @ u v 1 A = 0 @ t X t Y 1 A , (4.26)
where the 3D surfaces are assumed to only perform 3D translations and be captured by an orthographic camera. This 2 DOF model, which we name 2D-motion representation, is the most used in the optical flow literature. The warping function, W u (x, u):R 2 ⇥R 2 ! R 2 , is defined to map a pixel to its new location as a function of the perceived optical flow u, according to:

W o (x, u)=x + u. (4.27)
The scene flow, v(x), can be inferred using the optical flow u and the depth image Z(x), as follows:

v(x)=⇡ -1 (x + u,Z(x + u)) -⇡ -1 (x,Z(x)). (4.28)

Summary of motion representations

In this thesis, we consider three different motion representations, as is shown in Table 4.1. We are mainly concerned in the proposed twist-motion representation, which enables the direct use of the semi-rigid properties of the scene for the 3D motion field estimation. We also use the 3D-motion representation, which allows the direct solution of the scene flow and can be seen as particular case of the twist-motion representation. Both twist-motion and 3D-motion representations exploit depth data through a projective warping function. On the other hand, the 2D-motion representation considers an orthographic camera, neglecting the depth changes in the constraint of motion in the image, and yielding to a simpler formulation of the motion estimation problem. Although every representation have its own warping function, for the sake of simplicity, the generic expression W(x, •) is used, and the parameter ⇠, v, or u is used to identify the respective motion representation.

Representation Warping function

Small-increment model Scene flow

6-Twist W ξ (x, ξ)=π ⇣ e ξ X⌘ W ξ (x, ∆ξ)=x + J ξ (x)∆ξ ⇣ e ξ -I4⇥4 ⌘ X 3D-motion Wv(x, v)=π (X + v) Wv(x, ∆v)=x + Jv(x)∆v v 2D-motion Wo(x, u)=x + u Wo(x, ∆u)=x + ∆u π -1 (x + u) -π -1 (x)

Scene Flow Model

The RGBD representation of a scene changes when the observed 3D surfaces move or the camera undergoes motion, or both simultaneously. These changes in RGBD data, together with the definition of a set of invariants of the representation, allow the estimation of the 3D motion field presented in the scene. For example, the color patterns of a given surface are expected to remain unchanged when it moves in the scene, but also its color must differ from that of other surfaces, in order to perceive its motion in the RGB data. Similarly, when a rigid surface moves in the scene, it is expected that its 3D structure does not change and its intrinsic properties remain consistent in the observed depth data. However, at the same time, the representation should be able to capture the relative motion between the surface and the sensor. Figures 5.1 and 5.2 present examples of changes and invariances in the RGBD representation when surfaces move relative to the sensor. The definition of an appropriate RGBD-based framework for scene flow estimation makes it possible to take advantage of both sources of data: color and depth. However, there are several ways to define the scene flow model and different assumption can be used on the scene, in order to exploit the RGBD data. We base our framework on the assumption of a semi-rigid scene and on the use of constancy constraints directly on the RGBD images. Below we describe below our approach to exploit color and depth data for scene flow estimation, using the twist-and 3D-motion representations.

The scene flow estimation problem is formulated as the determination of the interframe 3D motion between two time steps. Accordingly, given two pairs of RGBD images {I c1 ,Z 1 } and {I c2 ,Z 2 } the goal is to solve for the 3D motion field that best explains the observed data, aiming to determine the real 3D motion executed in the scene. The use of an appropriate representation of motion makes possible to project 3D motions into the image domain and define consistency assumptions on color and depth data. This way plausible 3D-motion vectors can be solved to satisfy or minimize one or several RGBD-based constraints.

A motion constraint example

Let I k (x) be a color component of a scene point X, after an inter-frame 3D motion v, it is intuitive to assume that component I k remains unchanged, yielding to the constraint:

I k 1 (x)=I k 2 (⇡(X + v)). (5.1)
This constancy constraint is not linear on v, but can be linearized, using a first-order Taylor series expansion, to incrementally solve for the 3D motion. This linearization assumes that the image changes linearly with the projection of X + v, which is not always the case, especially under large motions in the image domain. This effect can be reduce by solving iteratively for small increments of v, where the linearization is more reliable. However, it is clear that only one equation as (5.1) does not uniquely determine the 3D motion v for a given scene point X: one equation with three unknowns. Since the full RGB information is available it seems reasonable to define similar constancy constraints to exploit the three color components. Nevertheless, the gradient information of the RGB component is highly correlated yielding in most cases to a ill conditioned system of equations that do not allow a reliable estimation of motion. Such limitation is known as the aperture problem in optical flow estimation and implies that only the component of the 2D motion, normal to The movement of the camera, as wel l as that of the person, produce a 3D motion of every scene point. Particularly, the head turning produces more changes in intensity than in depth, while the camera motion is better perceived on depth data in some regions, such as neck and t-shirt.

the gradient direction, can be determined from a single constraint. This disability is not exactly the same in the scene flow estimation problem for two main reasons: i) a 3D representation of the scene is provided by the depth data and ii) we aim to solve for a 3D motion rather than simply the 2D motion on the image. Therefore, it is evident that at least a second assumption is needed. The depth data can be used to formulate an additional constraint to ensure the consistency of the 3D motion with the 3D structure of the scene, as is presented in Section 5.2.2. Solving for the 3D motion using single constraints on depth and a color images, yields a kind of 3D aperture problem. Motion in the image domain can be determined only if is normal to the gradient direction in color or depth images. Moreover, color and depth gradient information could coincide yielding to the 2D aperture problem. Conversely, any change in the Z-component of the motion can be determined since it is directly reflected on the depth data, only limited by physics restrictions of the depth sensor. Besides, vanishing gradients, noise and outliers are also very common in color and depth images, making of the 3D motion computation an ill-posed problem that cannot be solved independently for each point.

We present below a variational framework to estimate scene flow by exploiting local and global properties of semi-rigid scenes.

Scene Flow Energy

In order to get a reliable estimation of the scene flow, it is required to impose a set of assumptions on the 3D motion field. In the case of optical flow estimation, mostly a local or global regularity is used to model motion. On one hand, a local method assumes that the motion in a local neighborhood on the image is constant. This way it is possible to locally solve for the motion that best explains the local data, ignoring the surrounding information. This idea dates from the seminal work by [START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF] for stereo registration, which have been widely used in optical flow. On the other hand, a global method assumes that the motion field is smooth. In their seminal work [START_REF] Horn | Determining optical flow[END_REF] introduce a global variational approach, where a (global) dependency is defined on the motion at every pixel with those at all the other pixels. Although local methods are robust to noise and can handle larger displacements, current global variational methods are by far the top performer in motion estimation tasks [START_REF] Sun | Secrets of optical flow estimation and their principles[END_REF]. Nevertheless, as is proposed by [START_REF] Bruhn | Lucas/kanade meets horn/schunck: Combining local and global optic flow methods[END_REF], it is also possible to benefit from local and global approaches in the same formulation. We follow this belief to exploit the semi-rigidity properties of real scenes from RGBD images. For this purpose, we use local methods to take advantage of the local rigidity of the scene while global methods are used to encourage piecewise smooth solutions, which are expected to be observed due to the assumed piecewise rigidity of the scene. Accordingly, we formulate the scene flow computation as the minimization problem of the energy:

E(⇠)=E D (⇠)+↵E S (⇠), (5.2)
where the twist-motion representation is used, aiming to better exploit the semirigidity of the scene and which, at the same time, being able to cover other motions models. The first term, E D (⇠), is the data term, which measures how consistent is the estimated scene flow with the observed intensity and depth data. The smoothness term, E S (⇠), is defined to favor smooth motion fields but preserving motion discontinuities. The local rigidity assumption can be exploited both in the data and smoothness terms. By using this formulation is possible to get an adjustable combination between local and piecewise rigidity in the solution. We describe below the definition of each of components of the scene flow energy, aiming to exploit the RGBD representation of the scene.

Data term

Our goal is to solve for the scene flow that best explains the observed color and depth data. The definition of the data term is based in two main assumptions: i ) the color of the scene points remains (approximately) unchanged and ii )t h ed e p t h data is able to reliably capture the 3D structure of the scene at every time step.

In order to constrain motion in the RGBD data the scene flow has to be projected to the image domain. Considering the twist-motion representation, the warping function given by Equation ( 4.3), models the motion of each scene point, enabling the formulation of RGBD constraints for elements of this twist field on the image.

Color-based constraint

As in most optical flow methods, we exploit the color invariance through the gray value image by defining a brightness constancy assumption (BCA):

I 2 (W(x, ⇠)) = I 1 (x), (5.3)
where I(x) represent the gray value of pixel x. Although it is also possible to define constraints by using each of the RGB components or another color space, the BCA has proved to be simple and effective. The main drawback of the BCA is its sensitivity to illumination changes, which are common in real scenes. Moreover, current RGBD sensors perform automatic white balance and exposure control by default, which can alter gray values when the composition of the scene change suddenly or when surfaces approach or move away from the camera. Therefore the invariance expected in (5.3) is often violated. In order to gain robustness against slight changes in brightness, [START_REF] Brox | High accuracy optical flow estimation based on a theory for warping[END_REF] define a gradient invariance of the gray value image. However this assumption fails under rotations, especially around the Z-axis since the gradient orientation changes. Nevertheless, the magnitude of the gradient, I g (x)=|rI(x)|, is expected to remain unchanged under a rotation, and in general under any rigid motion. Accordingly, we define a gradient constancy assumption (GCA) as:

I g 2 (W(x, ⇠)) = I g 1 (x).
(5.4) Figure 5.3 shows an example of how changes in the gray value image, not related with the motion in the scene, affects brightness and gradient images.

Depth-based constraint

Variations of the 3D structure, induced by the scene flow, must be consistent with the observed depth data at any time step. Therefore, it is possible to define a depth variation constraint (DVC ) as follows:

Z 2 (W(x, ⇠)) = Z 1 (x)+δ Z (x, ⇠), (5.5)
where δ Z (x, ⇠) denotes the variation in depth induced on the 3D point ⇡ -1 (x,Z 1 (x)) by the twist motion ⇠. Let δ 3D (x, ⇠)=( δ X , δ Y , δ Z ) T represent the 3D variation induced on x by twist motion ⇠, defined as:

δ 3D (x, ⇠)= ⇣ e ξ -I 4⇥4 ⌘ X, with X = 0 @ ⇡ -1 (x,Z 1 (x)) 1 1
A .

(5.6)

Accordingly, the variation in depth, δ Z , is obtained from the third component of vector δ 3D . This constraint enforces the consistency between the motion captured by the depth sensor and the estimated motion.

3D-motion case. If the 3D-motion representation is used for the scene flow, v =(v X ,v Y ,v Z ) T , the DVC takes a simpler form, which is given by:

Z 2 (W(x, v)) = Z 1 (x)+v Z , (5.7)
reflecting that the third component of the 3D motion field must be consistent with the observed change in depth.

Consistency measure

Different factors make the measurement of these consistency constraints a challenge. First, the presence of noise on color and depth data produces perturbations in the RGBD data, frequently violating the consistency terms. Second, lens distortions and the inaccuracy of the depth measures may bring errors in the projection of the 3D motion on the image. Third, partial or total occlusions of the 3D surfaces have to be handle since consistency constraints are not valid at these regions and could strongly affect the motion estimation. Finally, external changes uncorrelated to the 3D motion, as illumination changes, specular surfaces and shadows, produce inconsistency between the RGBD representation of the scene and the scene flow. For these reasons, RGBD based constraints are not obeyed in general and it is required to measure how consistent is a given estimation with the observed data.

Let ⇢ I , ⇢ g and ⇢ Z be the residuals of BCA, GCA and DVC constraints, respectively, given by:

⇢ I (x, ⇠)=I 2 (W(x, ⇠)) -I 1 (x) , (5.8) ⇢ g (x, ⇠)=I g 2 (W(x, ⇠)) -I g 1 (x) ,
(5.9)

⇢ Z (x, ⇠)=Z 2 (W(x, ⇠)) -(Z 1 (x)+δ Z (x, ⇠)).
(5.10)

In order to cope with outliers in brightness and depth brought by noise, occlusions or motion inconsistencies, a robust norm is required. In contrast to quadratic penalties that give too much influence to outliers in the motion estimation, l 1 norm has proved to be very effective to deal with outliers in the data term in global optical flow methods [START_REF] Wedel | Efficient dense scene flow from sparse of dense stereo data[END_REF]. We use the Charbonnier penalty Ψ s 2 = p s 2 +0.001 2 ,w h i c hi sad i fferentiable approximation of the absolute value, see Figure 2.9. This penalty is applied jointly to BCA and GCA, and separately to DVC, to get the total RGBD consistency measure: 5.11) where the constant γ balances brightness and gradient residuals, while λ weights intensity and depth terms.

⇢ data (x, ⇠)=Ψ ⇢ 2 I (x, ⇠)+γ⇢ 2 g (x, ⇠) + λΨ ⇢ 2 Z (x, ⇠) , (

Local consistency encouragement

It is well known that local methods are often more robust under noise than global methods, allowing reliable motion estimations even facing contradictory or missing data. Although local methods allow fast sparse solutions, they fail to provide dense estimations since the surrounding information is not exploited to disambiguate for motion in untextured regions. Local methods properties can benefit and complement the global motion estimation, as is proposed by [START_REF] Bruhn | Lucas/kanade meets horn/schunck: Combining local and global optic flow methods[END_REF].

Most scenes of interest can be well modeled as locally rigid scenes, i.e. they can be seen as scenes composed by independent 3D rigid parts. Under this context, the local assumption can benefit the motion estimation. However, the local assumption can bias the estimation if the representation of motion does not properly model these local properties. For example, under the rotation of a small rigid surface, neither the optical flow nor the 3D motion field should be straightforward encouraged to be locally consistent since motion patterns differ between them. Conversely, the twist representation is able to reflect the local rigidity of the scene on the image and a local estimation allows to directly exploit this properties. Therefore we define a local consistency error ⇢ l data , which measures how well a given twist ⇠ locally satisfies the RGBD constraints, according to:

⇢ l data (x, ⇠)= X x 0 2N (x)
w(x 0 ) ⇢ data x 0 , ⇠ (5.12)

where N (x) is an image neighborhood centering on x, and w is a specific weighting function. The local neighborhood N (x) can be setting to a fix size or adjusted as a function of the depth in order to cover a constant size in the scene. It is clear that if N (x)={x} that local measure ⇢ l data becomes the point-wise error ⇢ data . The weighting function can be chosen to follow a given distribution, e.g. uniform or Gaussian, or to be defined using the observed data. In this thesis we consider a fixed-size square function, for two cases: the uniform distribution, where weights are the unity, i.e., w(x 0 ) = 1, and a depth-driven weighting function, where weights decay with the difference in depth relative to the window center, defined as: (5.13) reflecting the belief that close scene points are more likely to follow the same rigid motion, and avoiding the inclusion of scene points from both sides of a depth edge.

w(x 0 )=e β|Z(x)-Z(x 0 )| ,
Finally, we modify the robust norm, which is able to deal with a tractable percentage of outliers, smaller than 50%, however, because total or partial occlusion of 3D surfaces, this percentage can be larger and skew motion estimation. In order to deal with of outliers brought by occlusion, the robust norm is truncated from a given value, s 2 tr , as is shown in Figure 5.4, cropping these scene points from the estimation to reduce the effect of partial occlusions.

Full data term

The full data term can be written as:

E D (⇠)= X x X x 0 2N (x)
w(x 0 ) ⇢ data x 0 , ⇠ .

(5.14)

Smoothness term

Without regularization, equation ( 5.11) alone is not sufficient to constrain the twist motion for a given point since there is an infinite number of twists that satisfy these constraints. On the other hand, equation ( 5.12) allows to solve for a local motion on well textured regions but in most cases fails to capture the full motion field. We exploit the semi-rigidity of the scene to define a regularization term on the field of twist motions, which aims to favor piecewise smooth solutions. Particularly, we use TV as regularizer on this field of rigid motions for three reasons: i ) does not penalizes discontinuities, or edges, ii ) does not penalizes smooth functions and iii ) preserves exactly the location of discontinuous edges. Moreover there exist efficient solvers for TV-based regularizers.

Regularization on the twist field

A twist motion field ⇠ can be decomposed into a rotational field !,w i t he ω : Ω ! SO(3), and a 3D motion field ⌧ : Ω ! R 3 . Since they are decorrelated by nature, we regularize each field independently using a weighted TV, which allows piecewise smooth solutions while preserving motion discontinuities. Although discontinuities of these two fields are expected to coincide, a decoupled regularization scheme is preferred for three reasons: i ) it allows to benefit from available, efficient TV regularizers, simplifying the optimization of the energy, ii ) it yields a general framework for the twist-motion and 3D-motion representations of the scene flow and iii ) it let to apply different TV regularization strategies to each component. Moreover, thanks to the depth data is possible to define a weighting function to preserve discontinuities.

Depth-based discontinuities preserving

In a semi-rigid scene it is expected that close 3D points move under the action of the same rigid body motion, while distant points may or may not follow the same rigid motion. However, in most cases, discontinuities of the motion field coincide with the boundary of the observed 3D surface. This fact can be exploited using the 2D parametrization of the 3D surface given by the depth image Z(x). Particularly, the magnitude of the gradient |rZ(x)| encodes the 3D distance between neighbor scene points and expose the edges of the surface as is shown in Figure 5.5.W e define the decreasing positive function c(x)=e -β|rZ(x)| 2 (5.15) to quantify the continuity of the surface in a given point x, as is shown in Figure 5.5. The use of the weighting function c(x) prevent the regularization across the discontinuities of the surface, as is presented below.

Translational field regularization

Regularization of the fields ! and ⌧ poses different challenges. Elements of ⌧ lie in the Euclidean space R 3 , and the problem corresponds to a vector-valued function regularization. Different vectorial TV approaches can be used to regularize ⌧ ,s u c h as those described in [START_REF] Goldluecke | An approach to vectorial total variation based on geometric measure theory[END_REF]. Particularly, the channelby-channel l 1 norm has been successfully used for optical flow [START_REF] Wedel | Statistical and Geometrical Approaches to Visual Motion Analysis[END_REF].

We follow this approach to define the weighted TV of ⌧ as: 5.16) where kr⌧ k

TV c (⌧ )= X x c(x) kr⌧ (x)k , (
:= |r⌧ X | + |r⌧ Y | + |r⌧ Z | and c(x)
is the weighting function (5.15), which helps to preserve motion discontinuities along edges of the 3D surface. The absolute value can be replaced by the Huber norm [START_REF] Werlberger | Anisotropic Huber-L1 Optical Flow[END_REF], which is given by (2.28), to reduce the staircasing effect. Efficient solvers for both norms exist and some of them are presented in [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF]. The staircasing effect can also be reduced if TV is replaced by the total generalized variation (TGV) [START_REF] Bredies | Total generalized variation[END_REF], where higher order derivatives of ⌧ are not penalized.

Rotational field regularization

Elements of ! are rotations in the Lie group SO(3), embedded in R 3 through the exponential map (2.7), and the regularization has to be done on this manifold. In order to apply a TV regularization, a notion of variation should be used for elements of SO(3), which is not a vector space. Given two rotations e ω1 ,e ω2 2 SO(3), the residual rotation can be defined as e ω2 e -ω 1 , representing the rotation required to go from e ω1 to e ω2 in SO (3). The product in logarithmic coordinates can be expressed as e ω1 e ω2 = e µ(ω 1 ,ω 2 ) , where the mapping µ can be expanded in a Taylor series around the identity, using the Baker-Campbell-Hausdorff formula: .17) where [•, •] is the Lie bracket in so (3). Close to the identity, equation (5.17)i sw e l l approximated by its first-order terms, yielding to e ω2 e -ω 1 ⇡ !2 -!1 . Therefore, for small rotations the variation measure can be defined as the matrix subtraction in so(3), or equivalently, as a vector difference for the embedding in R 3 . Accordingly, with r being the image gradient operator, the derivative image matrix, D! := (r! X , r! Z , r! Z ) T : Ω ! R 3⇥2 , approximates the horizontal and vertical pointwise variations of ! on the image. Following [START_REF] Goldluecke | An approach to vectorial total variation based on geometric measure theory[END_REF], we define the TV as the sum over the largest singular value σ 1 of the derivative matrix: .18) This TV approach supports a common edge direction for the three components, which is a desirable property for regularization of the field of rotations. Moreover, deviations are less penalized with respect to other measures, e.g., the Frobenius norm [START_REF] Rosman | Group-valued regularization for analysis of articulated motion[END_REF], and efficient solvers are available. However, this TV definition approximates the real structure of the manifold, yielding biased measures when it is far from the identity. The more the rotation is away from the identity, the more its variations in SO( 3) are penalized as is shown hereinafter. Given two rotations

µ(! 1 , !2 )=! 1 +! 2 + 1 2 [! 1 , !2 ]+O(|(! 1 , !2 )| 3 ), ( 5
TV σ (!)= X x c(x) σ 1 (D!(x)) . ( 5
! 1 = ✓ 1 - ! 1 and ! 2 = ✓ 2 - ! 2
,w i t h✓ the angle and -! the unitary axis vector of the rotation, and writing

✓ 2 = ✓ 1 +δ✓, leads to ! 2 -! 1 = ✓ 1 ( - ! 2 -- ! 1 )-δ✓ - ! 2
.T h i s linearly dependent penalization usually is not a problem, since large rotations imply larger motion on the image. Thus, a stronger regularization can be reasonable. Moreover, large rotation caused by a global motion of the scene or the camera can be optimized separately and compensated, as we present in Section 5.6. Optionally, this over-penalization can be removed by expressing each rotation as ! = ✓ -! and applying a vectorial TV on -! and a scalar TV on ✓. This alternative definition does not increase the computational cost much, since the solver for (5.18) computes the norm of every element of ! and the field generated by ✓ is efficiently regularized using a scalar TV solver.

Full regularization term

The regularization on the twist field is formulated as:

E S (⇠)=TV c (⌧ )+TV σ (!).
(5.19)

Regularization on the 3D motion field

When a rigid surface moves in the scene, the 3D distance between any two of its points remains constant. Therefore, it is intuitive to expect that close points perform the same, or a very similar, 3D motion, and encourage piecewise smooth solutions through the regularization. This assumption has been successfully used in the regularization of the optical flow. However, although 3D distances remain unchanged, 3D motions performed for points of the moving rigid surface are not the same in general. Using the small-rotation model for a rigid motion, given by Equation (2.16), the scene flow v induced on a 3D point X, is given by: 0

B B @ v X v Y v Z 1 C C A = 0 B B @ 0 -ZY Z 0 -X -YX 0 1 C C A 0 B B @ ! X ! Y ! Z 1 C C A + 0 B B @ t X t Y t Z 1 C C A , (5.20)
from which is clear that if the rotational component is not zero, the scene flow is a function of the 3D point, and so that, it differs from one point to another on the rigid surface. Only a zero rotation component of the rigid motion, i.e., a pure translational motion, guarantees a constant scene flow. Using the 3D-motion representation of the scene flow is not possible to directly exploit the piecewise rigidity of the scene. However, close 3D points that follows the same rigid motion, are expected to present similar 3D motions and piecewise smoothness of the 3D motion field is a valid approximation. Therefore, we propose a regularization approach using this piecewise smoothness assumption combined with a rigid motion prior [START_REF] Vogel | 3D scene flow estimation with a rigid motion prior[END_REF], which encourages a local rigidity without an explicit computation of the rigid motion field. This way local and piece-wise rigidity are combined into the regularization of the 3D motion field.

TV-based regularization on the 3D motion field

In order to favor smooth piecewise solution of the scene flow, a weighted TV is applied. The regularization is done channel by channel according to:

TV c (v)= X x c(x) krv(x)k , ( 5

Rigid-motion prior on the 3D motion field

Aiming to exploit the semi-rigidity of the scene, the 3D motion field is encouraged to be locally consistent with a rigid motion. Let NR(v) be the nonrigid residual functional, defined as:

NR (v)= X x Ψ ⇣ Rig N v (x) -v(x) 2 ⌘ (5.22)
where Rig N v (x) is the local N ⇥ N projection of the scene flow v onto the closest rigid motion subspace. The functional NR (v) measures the total non-local rigidity of the scene flow and it is provided with the robust norm Ψ s 2 = p s 2 +0.001 2 , to deal with outliers. The definition of the local rigid projection Rig N v (x) is based on the local fitting of a rigid motion (!, t), using the small-rotation approximation (5.20). With N (x) being a N ⇥ N neighborhood centering at x, the local rigid motion (!, t) is defined via the following weighted least-squares problem:

min (ω,t) X x 0 2N (x) w(x, x 0 ) [ X 0 ] ⇥ | I 3⇥3 0 @ ! t 1 A -v(x 0 ) 2 ,
(5.23) with X 0 = ⇡ -1 (x 0 ,Z(x 0 )) and [•] ⇥ the cross product matrix. The weighting function w x (x, x 0 ) aim to favor points belonging to the same surface of x. For this purpose, each pixel x 0 is weighted by e -kd(x,x 0 ) 2 ,w h e r ed(x, x 0 ) is the 3D Euclidean distance between x and x 0 , and k an adjustable parameter. We follow the solution of problem (5.23)p r e s e n t e db yVogel et al. (2011). Let M 2 R 3N 2 ⇥6 and P 2 R 3N 2 ⇥1 be the concatenation of matrices ([ X 0 ] ⇥ | I 3⇥3 ) 2 R 3⇥6 and scene flow vectors v(x 0 ) 2 R 3 , respectively, whose elements are taken from N x . Let µ 2 R 3N 2 ⇥3N 2 be the diagonal matrix, whose elements correspond to the weights given by w x (x, x 0 ). Using this notation, the rigid motion that solves (5.23) is given by: (5.24) therefor, the projection of v(x) onto the local rigid motion subspace is given by:

0 @ ! t 1 A =(M T µM ) -1 M T µP ,
Rig N v (x)=([X] ⇥ | I 3⇥3 )(M T µM ) -1 M T µP .
(5.25)

The rigid projection requires the inversion of the matrix (M T µM ) 2 R 6⇥6 ,w h i c h can be done efficiently. Although equation ( 5.25) allows the projection of a local scene flow without the explicit computation of the rigid motion, due to the linear approximation it assumes a small rotation. It is also possible to iteratively compute the full rigid motion that best fit the local scene flow, in which case, is required to explicitly compute (!, t).

Full regularization term

The regularization term of the scene flow is defined as: 5.26) where ⇣ balances the total variation and the locally-rigid terms.

E S (v)=TV c (v)+⇣ NR(v), (

Optimization of the energy

The minimization of the proposed energy (5.2) is a challenge due to the l 1 -norm used in the regularization term, which is non-differentiable. However, it is possible to introduce an auxiliary variable [START_REF] Wang | A new alternating minimization algorithm for total variation image reconstruction[END_REF] to decompose the optimization into two simpler problems. In particular, the auxiliary twist field χ is used to reformulate the minimization problem as follows:

min ξ,χ E D (⇠)+ 1 2 X x |⇠(x) -χ(x)| 2 + ↵E S (χ) (5.27)
where  is a small numerical variable. Note that the linking term between ⇠ and χ is the distance on the tangent space at the identity in SE(3). The solution of (5.27) converges to that of (5.2) as  ! 0. Minimization of this energy is performed by alternating the following two optimization problems:

i. For fixed χ, estimate ⇠ that minimizes:

X x ⇢ l data (x, ⇠(x)) + 1 2 |⇠(x) -χ(x)| 2 .
(5.28)

ii. For fixed ⇠ =(!, ⌧ ), compute χ =($, ⇡) that minimizes:

( TV c (⇡)+ ⌘ 2 X x |⇡(x) -⌧ (x)| 2 ) + ( TV σ ($)+ ⌘ 2 X x |$(x) -!(x)| 2 ) , (5.29)
where ⌘ =(↵) -1 . Below we present how to solve for data-based and smoothnessbased problems, which are defined as the minimization of (5.28) and ( 5.29), respectively. Afterwards, the minimization of the scene flow energy for the 3Dmotion representation is presented.

Data-based energy minimization

The optimization of (5.28) can be done pointwise by minimizing: .30) This minimization corresponds to a nonlinear optimization problem, which can be solved by IRLS [START_REF] Green | Iteratively Reweighted Least Squares for Maximum Likelihood Estimation, and some Robust and Resistant Alternatives[END_REF] with a Gauss-Newton numerical algorithm. For this purpose, the energy is linearized at a given estimate ⇠, to calculate a rigid motion increment, ∆⇠, on the tangent space to the manifold. This procedure is applied iteratively and the twist motion updating is done via ⇠ log(e ∆ξ e ξ ).

1 2 |⇠ -χ| 2 + X x2N (x 0 ) Ψ ⇢ 2 I (x, ⇠)+γ⇢ 2 g (x, ⇠) + λΨ ⇢ 2 Z (x, ⇠) . ( 5

Linearization of the data-based energy

Energy (5.30) can be linearized around the current estimate ⇠, using a first-order Taylor series expansion. Particularly, for the brightness residual we have: 5.31) where I x =( @I 2 /@x, @I 2 /@y) is the image gradient and @W/@⇠ is the Jacobian J of the warp, which is given by (4.10), and both evaluated at the current estimate of the scene point position, x ξ = W(x, ⇠). Therefore, the residual becomes:

⇢ I ⇣ x, log(e ∆ ξ e ξ ) ⌘ ⇡ I 2 (x ξ )+I x @W @⇠ ∆⇠ -I 1 (x) , (
⇢ I ⇣ x, log(e ∆ ξ e ξ ) ⌘ ⇡ ⇢ I (x, ⇠)+I x J∆⇠, (5.32)
and, similarly, for the gradient residual we have:

⇢ g ⇣ x, log(e ∆ ξ e ξ ) ⌘ ⇡ ⇢ g (x, ⇠)+I g x J∆⇠.
(5.33)

Using the first-order expansion for the depth residual yields: 5.34) where Z x =( @Z 2 /@x, @Z 2 /@y) is the depth image gradient evaluated at x ξ .T h e last term, δ Z , is the variation in depth which can be linearized by approximating the 3D motion, δ 3D ,i n d u c e db y∆⇠, as follows:

⇢ Z ⇣ x, log(e ∆ ξ e ξ ) ⌘ ⇡ Z 2 (x ξ )+Z x J∆⇠ - ⇣ Z 1 (x)+δ Z ⇣ x, log(e ∆ ξ e ξ ) ⌘⌘ , (
δ 3D (x, log(e ∆ ξ e ξ )) = ⇣ e ∆ξ e ξ -I ⌘ X ⇡ ⇣ (I + ∆⇠)e ξ -I ⌘ X ⇡ ⇣ e ξ -I ⌘ X + ∆⇠ ⇣ e ξ X⌘ = δ 3D (x, ⇠)+ ∆⇠ Xξ . (5.35)
Therefore, the depth residual can be written as:

⇢ Z ⇣ x, log(e ∆ ξ e ξ ) ⌘ ⇡ ⇢ Z (x, ⇠)+(Z x J -K) ∆⇠, (5.36)
where the 1 ⇥ 6 vector K is defined as K

= D ([X ξ ] ⇥ |I 3⇥3 ), with D =( 0 , 0, 1)
isolating the third component and [•] ⇥ denoting the cross product matrix. Finally, the linearized version of the data term can be expressed as: 5.37) where the parameters (x, ⇠) are considered implicit.

1 2 log(e ∆ ξ e ξ ) -χ 2 + X x2N (x 0 ) Ψ ⇣ |⇢ I + I x J∆⇠| 2 + γ |⇢ g + I g x J∆⇠| 2 ⌘ + λΨ ⇣ |⇢ Z +(Z x J -K) ∆⇠| 2 ⌘ , (

Twist motion increment

Taking the partial derivative of (5.37)withrespectto∆⇠ and setting it to zero, the increment ∆⇠ can be computed as: 5.38) where the 6 ⇥ 1 gradient vector g is given by:

∆⇠ = -H -1 g, (
g = X x2N (x 0 ) Ψ 0 ⇢ 2 I + γ⇢ 2 g ⇥ (I x J) T ⇢ I + γ(I g x J) T ⇢ g ⇤ +λ Ψ 0 ⇢ 2 Z (Z x J -K) T ⇢ Z o + 1  log ⇣ (e ξ ) -1 e χ⌘ , (5.39)
with Ψ 0 the derivative of the robust norm, which is always evaluated at the current estimate ⇠.T h e6⇥ 6 matrix H is the Gauss-Newton approximation of the Hessian matrix, which is given by:

H = X x2N (x 0 ) Ψ 0 ⇢ 2 I + γ⇢ 2 g h (I x J) T (I x J)+γ (I g x J) T (I g x J) i + λ Ψ 0 ⇢ 2 Z (Z x J -K) T (Z x J -K)+ 1 2 I 6⇥6 . (5.40)
A reliable computation of (5.38) requires that the Hessian matrix H is invertible and well conditioned. The construction of H is done as a weighted linear combination of 6 ⇥ 6 matrices, which are computed using gradient images I x , I g x and Z x , and the Jacobian J. The main drawback of local methods following the Lucas-Kanade framework [START_REF] Baker | Lucas-Kanade 20 years on: A unifying framework[END_REF], occurs on regions presenting vanishing gradients, in both intensity and depth images for RGBD data. In such that case, the Hessian matrix H is non-invertible and/or ill-conditioned, and so that, being unable to perform dense motion estimation. However, in our approach, the regularization applied includes the auxiliary field χ, which provides an interesting stability to the optimization procedure. The identity matrix I 6⇥6 allows to smoothly switch between the Gauss-Newton algorithm and a gradient descent algorithm. On one hand, if the gradients vanish, H tends to a constant diagonal matrix and the twist increment equation (5.38) leans to a gradient descent updating where the initial estimation is favored. On the other hand, if the gradient information is meaningful, H -1 guides the twist motion increment according to the gradient images and the Jacobian, in a Gauss-Newton method. This way, equation (5.38) avoids to estimate motions on flat regions of the image, giving to the smoothness procedure the task of completing the estimation. Unlike the 2D motion on the image plane, depth changes can always be computed as is demonstrated in Section 5.4.3, where a close look to the Hessian Matrix is done for the 3D-motion representation case.

Smoothness-based energy minimization

Each side of the energy (5.29) can be seen as vectorial image denoising problem with a TV-L 2 model. The left side corresponds to a channel-by-channel TV. This regularization is done independently per channel and coincides with the ROF-model [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF]. This problem is solved using the Chambolle's projection algorithm [START_REF] Chambolle | Total variation minimization and a class of binary mrf models[END_REF], which is based on gradient descent and re-projection on the dual-ROF model. The right side is a vectorial TV regularization proposed by [START_REF] Goldluecke | An approach to vectorial total variation based on geometric measure theory[END_REF]. This problem can be solved using the Bermudez-Moreno algorithm [START_REF] Bermudez | Duality methods for solving variational inequalities[END_REF][START_REF] Goldluecke | The natural vectorial total variation which arises from geometric measure theory[END_REF], which is a generalization of the Chambolle's projection algorithm [START_REF] Aujol | Some first-order algorithms for total variation based image restoration[END_REF]. Below we describe the solution of both TV-based problems.

Channel-by-channel TV

Given a fixed 3D motion field ⌧ =( ⌧ X , ⌧ Y , ⌧ Z ), the goal is to compute ⇡ = (⇡ X , ⇡ Y , ⇡ Z ) that minimizes the left side of the energy (5.29). For every channel i = {X, Y, Z}, this problem corresponds to the following weighted formulation of the ROF image denoising model:

X x c(x) |r⇡ i (x)| + ⌘ 2 |⇡ i (x) -⌧ i (x)| 2 .
(5.41)

An efficient solution, using the dual ROF-model, is given by the Chambolle's projection algorithm [START_REF] Chambolle | Total variation minimization and a class of binary mrf models[END_REF] according to:

⇡ i (x)=⌧ i (x)+ c(x) ⌘ r • p(x) , (5.42)
where p =( p 1 ,p 2 ):Ω ! R 2 is the solution of the dual problem, which is defined (almost everywhere) as the unite vector r⇡ i /|r⇡ i |.S e t t i n gp 0 = 0, ⇡ 0 i = ⌧ i and choosing ⌫  1/4, the dual variable p can be computed recursively following:

p n+1 (x)=Π S (p n (x)+(⌘⌫) r⇡ n i (x)) , (5.43)
with the projection Π S defined as:

Π S (p(x)) = p(x) max {1, |p(x)|} .
(5.44)

Vectorial TV

Given a rotational field ! =(! X , ! Y , ! Z ) one wants to solve for $ =($ X , $ Y , $ Z ) that minimizes: (5.45) This problem is solved using the Bermudez-Moreno's framework [START_REF] Aujol | Some first-order algorithms for total variation based image restoration[END_REF]. This is a general projection-based algorithm that can be used to solve the dual problem of a vectorial TV image denoising, as is presented by [START_REF] Goldluecke | The natural vectorial total variation which arises from geometric measure theory[END_REF]. Using the Bermudez-Moreno's algorithm the minimizer of (5.45) is given by: 5.46) where P =(p 1 , p 2 , p 3 ):Ω ! R 3⇥2 is the solution of the dual problem, and Div P, the divergence matrix defined as Div P := (r • p 1 , r • p 2 , r • p 3 ) T : Ω ! R 3 .S e t t i n g P 0 = 0, $ 0 = ! and choosing ⌫  1/4, the dual variable P can be computed recursively following:

X x c(x) σ 1 (D!(x)) + ⌘ 2 |$(x) -!(x)| 2 ,
$(x)=!(x)+ c(x) ⌘ Div P(x) , (
P n+1 (x)=Π K (P n (x)+(⌘⌫)D$ n (x)) , (5.47)
with D! := (r! X , r! Z , r! Z ) T : Ω ! R 3 , the derivative matrix of !, and the projection, Π K , defined as:

Π K (P(x)) = PVΣ + Σ l 1 V T .
(5.48)

Being P = UΣV the singular value decomposition of P, and Σ, the diagonal matrix containing the singular values {σ 1 , σ 2 }, matrix Σ + denotes the pseudo inverse of Σ and matrix Σ l 1 the projection of the elements of Σ onto the l 1 -unit ball.

Minimization on the 3D motion field

Using the 3D-motion representation of the scene flow, the problem is formulated as the minimization of the following energy:

E(v)=E D (v)+↵ {TV c (v)+λ NR(v)} , (5.49)
where the data term, E D , is given by (5.14), with ⇠ =(0, v), i.e., a zero rotational component, and the smoothness term is formed by the total variation and nonrigid terms. The minimization of the energy (5.49) is done by introducing two auxiliary motion fields to decompose the optimization into three simpler problems. Let u and w be 3D motion fields, then the scene flow v is solved to minimize:

E D (v)+ 1 2 X x |v(x) -u(x)| 2 + ↵TV c (u)+ 1 2 X x |v(x) -w(x)| 2 + ↵λ NR(w).
(5.50) The minimum of (5.50) converges to the solution of (5.49) as  ! 0, and can be computed by alternating the three following optimization problems: i. For fixed u and w, estimate v that minimizes: .51) ii. For fixed v, compute u that minimizes: (5.52) where ⌘ =(↵) -1 .

X x ⇢ l data (x, v(x)) + 1 2 n |v(x) -u(x)| 2 + |v(x) -w(x)| 2 o . ( 5
TV c (u)+ ⌘ 2 X x |u(x) -v(x)| 2 ,
iii. For fixed v, compute w that minimizes:

NR(w)+ ⌘ 2 X x |w(x) -v(x)| 2 , (5.53) with ⌘ =(↵λ) -1 .
Problem ii. is solved using the channel-by-channel TV approach presented in Section 5.4.2. The minimization of problem i., or data-based energy, and of problem iii., or nonrigid-based energy, is presented below.

Data-based energy minimization

The optimization of (5.51) can be done pointwise by minimizing:

1 2 n |v -u| 2 + |v -w| 2 o + X x2N (x 0 ) Ψ ⇢ 2 I (x, v)+γ⇢ 2 g (x, v) + λΨ ⇢ 2 Z (x, v) .
(5.54) This problem can be solved using IRLS with a Gauss-Newton algorithm, where the scene flow increment, ∆v, is iteratively calculated by linearizing the energy at the current estimate v. Thus, the scene flow solution is updating via v v + ∆v. Following (5.37), energy (5.54) is linearized around the current estimate v as follows: 5.55) where the Jacobian J of the warp, is given by Equation ( 4.16) and the 3D vector D =(0, 0, 1) isolates the Z-component. Taking the partial derivative of (5.55)w i t h respect to ∆v, and setting it to zero, the increment ∆v can be computed as: 5.56) where g is the 3 ⇥ 1 gradient vector defined as: 5.57) with the derivative of the robust norm Ψ 0 evaluated at the current estimate v, and the 3 ⇥ 3 approximation of the Hessian matrix H given by:

1 2 n |v -u| 2 + |v -w| 2 o + X x2N (x 0 ) Ψ ⇣ |⇢ I +(I x J)∆v| 2 + γ |⇢ g +(I g x J)∆v| 2 ⌘ + λΨ ⇣ |⇢ Z +(Z x J -D) ∆v| 2 ⌘ , (
∆v = -H -1 g, (
g = X x2N (x 0 ) Ψ 0 ⇢ 2 I + γ⇢ 2 g ⇥ (I x J) T ⇢ I + γ(I g x J) T ⇢ g ⇤ +λ Ψ 0 ⇢ 2 Z (Z x J -D) T ⇢ Z o + 2  ✓ v - u + w 2 ◆ , (
H = X x2N (x 0 ) Ψ 0 ⇢ 2 I + γ⇢ 2 g h (I x J) T (I x J)+γ (I g x J) T (I g x J) i + λ Ψ 0 ⇢ 2 Z (Z x J -D) T (Z x J -D)+ 1  I 3⇥3 . (5.58)
The construction of H is a weighted linear combination of 3 ⇥ 3 matrices, which are computed from gradient images I x , I g

x and Z x , and the Jacobian, J, of the warp. For this reason, the gradient information determines the ability to estimate motion in the image domain for a given region. For instance, each of the matrices which depends on the intensity gradient, can be written as: 5.59) with I Σ = -(xI x + yI y ), combining horizontal and vertical components of the gradient as a function of the pixel location. The 2 ⇥ 2 inner matrix in (5.59), I T x I x , characterizes the gradient information pointwise, and when summed with all matrices in the neighborhood, produces the covariance matrix of the gradient,which has been used to determine regions that can be tracked well [START_REF] Shi | Good features to track[END_REF]. The full 3 ⇥ 3 matrix takes into account the projective transformation depending on the pixel location, and the inverse square of its depth, which weighs the contribution of every matrix. If the point-wise gradient vanishes or Ψ 0 ! 0, due to the presence of an outlier, there is a null contribution of this pixel into the overall construction of H. Gradient and depth images also contribute in the solution, particularly, the gradient of the depth image brings a decorrelated information, capturing the texture of the 3D surface which benefits the characterization of the region. The matrices depending on the depth gradient take the following form:

(I x J) T (I x J)=J T 0 @ I 2 x I x I y I x I y I 2 y 1 A J = 1 Z(x) 2 0 B B @ I 2 x I x I y I x I Σ I x I y I 2 y I y I Σ I x I Σ I y I Σ I 2 Σ 1 C C A (
(Z x J) T (I x J)= 1 Z(x) 2 0 B B @ Z 2 x Z x Z y Z x (Z Σ -1) Z x Z y Z 2 y Z y (Z Σ -1) Z x (Z Σ -1) Z y (Z Σ -1) (Z Σ -1) 2 1 C C A (5.60)
with Z Σ = -(xZ x + yZ y ). This way gradient information from intensity and depth images complement each other, allowing motion estimation on regions where the use alone of each source of data cannot. However, when both intensity and depth gradients vanish, it is not possible to estimate motion in the image since I x J, I g x J and Z x J in equation ( 5.56) go to zero. Observe that even on this untextured region the matrix H is still invertible thanks to the regularization. This is an important properties because it prevents the estimation of the 2D motion on regions with no enough gradient information, conversely, it is still able to compute the Z-component of the motion. Observe that when gradients vanish the 2D motion estimation is done in the regularization step, which corrects and completes the estimation using the global information. However, at the same time, if intensity and depth gradients vanish the increment equation (5.56) 5.61) allowing to compute the changes in depth by minimizing the depth residual ⇢ Z .

becomes: ∆v = H -1 2 4 X x2N (x 0 ) λ Ψ 0 ⇢ 2 Z D T ⇢ Z + 2  ✓ v - u + w 2 ◆ 3 5 , (

nonrigid-based energy minimization

The minimization of (5.53) is done per-pixel by solving: (5.62) where the parameter x is considered implicit. By setting the initial estimate as w = v, the scene flow can be iteratively computed via w w + ∆w,w h e r et h e 3D motion increment, ∆w, is given by: .63) with ⇢ NR = Rig N ww, the nonrigid residual, and Ψ 0 the robust norm derivative.

min w Ψ ⇣ Rig N w -w 2 ⌘ + ⌘ 2 |w -v| 2 ,
∆w = 2Ψ 0 ⇣ |⇢ NR | 2 ⌘ ⇢ NR + ⌘(v -w) 2+⌘ , ( 5

Coarse-to-fine estimation

Current RGBD sensors are able to provide more than 15 registered images per second, so that, inter frame 3D motions are expected in the range of few centimeters. However, the displacement on the image brought by the projection of the motion of a scene point becomes larger, the closer the 3D surface is from the camera, yielding quiet often to motions of more than one pixel of magnitude. Similarly, when the sensor itself moves rigidly relative to the scene, further scene points can induce large displacements in the image domain, even for a subtle camera motion. In such those cases and many others, images displacements are larger than one pixel between two frames and the scene flow energy is expected to be multi-modal, i.e. having multiple local minima. For this reason, the proposed minimization can be easily stuck in a local minimum. In order to find the global minimum, it is possible to apply multi-scale strategy [START_REF] Anandan | A computational framework and an algorithm for the measurement of visual motion[END_REF][START_REF] Brox | High accuracy optical flow estimation based on a theory for warping[END_REF]: solving coarser and smoothed versions of the problem, each of which may have a unique minimum, hopefully close to the global minimum of the original problem. In this strategy, every coarser version of the problem is formulated by downsampling the original images and its solution is used to initialize the solution of a finer version.

This multi-scale strategy is adopted in order to deal with larger motions. We construct a RGBD image pyramid with a downsampling factor of 2. A Gaussian anti-aliasing filter is applied to the intensity image at the original resolution and the intensity image pyramid is built using bicubic downsampling. For the depth image, a 5 ⇥ 5 median filter is used at the original resolution and the depth image pyramid is constructed by averaging pixels in non-overlapped neighborhoods of 2⇥2, where only pixels with a valid depth measure are used. Averaging on depth images is avoided to prevent the propagation of pixels having no valid depth measures. Having a pyramid with levels l = {0, 1,...,L}, with 0 the original resolution, the computation is started at level L an the estimated twist field is directly propagated to the next lower level, as is shown in Figure 5.6, for a pyramid of 2 levels. The camera matrix is scaled at each level by the factor 2 l , compensating the changes in the projection of scene points due to downsampling. At each level M loops are performed, consisting of M GN iterations of the Gauss-Newton procedure followed by M TV iterations of the TV solver. The linking constant  is styled at each scale. Algorithm 1 presents the general scheme of the scene flow estimation using the twist-motion representation. In order to reduce computational time, the domain of ⇠ is set to correspond to the domain of RGBD images at each level of the pyramid. For this reason, before moving to a finer level, a linear interpolation is applied on ⇠ to match the appropriate image domain.

The same multi-scale strategy can be applied if the 3D-motion representation is used. The scene flow v is filtering by component with a 5 ⇥ 5 median filter, before being propagated from one level to the finer one, robustly integratating flow estimates over large spatial neighborhoods [START_REF] Sun | Secrets of optical flow estimation and their principles[END_REF]. Algorithm 2 presents the general scheme of the scene flow estimation using the 3D-motion representation.

Rigid plus nonrigid model

In many applications, the sensor itself moves relative to the observed scene and causes a dominant global motion in the overall motion field. In this situation, compensating for the motion of the camera can simplify the estimation and regularization of the scene flow. Figure 5.7 presents an example of how global rigid compensation modifies the motion estimation problem. Moreover, for some applications, such as SLAM and 3D reconstruction, the camera ego-motion is required. Particularly, for 3D reconstruction of deformable objects, the camera motion is needed to register partial 3D reconstructions, which are constructed by compensating the nonrigid motions. Therefore, we consider splitting the motion of the scene into a globally rigid component, motion relative to the dominant object/background, plus a nonrigid residual field, ⇠ =(⌧, !). We assume that a large part of the scene follows the same rigid motion. Accordingly, the scene flow can be recovered from the composition of twist fields given by χ = log(e ξ + e ξR -I 4⇥4 ), and the estimation problem is formulated as: 5.64) with E Rig (χ) and E Res (χ) the rigid and nonrigid energies, respectively. It is worth noting that the separation of the camera motion is in addition to the framework presented above, i.e., the nonrigid part can still deal with local motion.

⇠ R =(⌧ R , ! R ) 2 R 6 , capturing the camera w(x, ξ) w(x, ξ) w(x, ξ) ξ ξ ξ Cam = {f /2, c/2} Cam = {f /4, c/4} Cam = {f , c} S 0 1 = {I 1 ,Z 1 } S 0 2 = {I 2 ,Z 2 } S 1 1 S 2 1 S 2
min χ E Rig (χ)+E Res (χ) , (

Rigid energy

The estimation of a global rigid motion can be seen as a particular case of the framework given by (5.2), where every twist of the motion field follows the same rigid motion and so that, no regularization is needed. The camera motion, or the dominant rigid motion, can be estimated using the data term (5.14), by considering every pixel (or a subset of Ω) to solve for a unique twist ⇠ R in a robust least-squares estimation. Accordingly, the rigid component of the energy is defined as: 5.65) where w I (x) and w Z (x) weigh the per-pixel contribution of intensity-and depthbased terms. Setting w I (x) = 1 and w Z (x)=λ yields the same formulation given by (5.12), but in general they can be defined as per-pixel functions of the confidence in intensity and depth data. Also it is possible to use only intensity or depth constraints. For example, setting w I (x) = 0 and w Z (x) = 1 corresponds to estimate a rigid transformation consistent only with the observed 3D surface.

E Rig (χ)= X x w I (x)Ψ ⇢ 2 I (x, χ)+γ⇢ 2 g (x, χ) + w Z (x)Ψ ⇢ 2 Z (x, χ) , (

Nonrigid energy

The residual nonrigid motion can be computed following the framework (5.2) for scene flow estimation. Accordingly, the nonrigid energy is given by:

E Res (χ)=E D (χ)+↵E S (χ).
(5.66)

Rigid plus nonrigid estimation

The estimation of rigid and nonrigid components of the motion is a challenge since both are superimposed in the overall motion field affecting its separately estimation.

If the rigid body motion is extremely dominant, covering by far most of the observed scene, it can be estimated independently as is done in the example of Figure 5.7. However, if the nonrigid motion is dominant, the independent estimation of the rigid motion is strongly biased by the nonrigid component. Similarly, the rigid component, especially when is brought by the camera motion, can produce large displacements on the image making harder the separately estimation of the nonrigid component. In order to reliably estimate the rigid and nonrigid components of the scene flow, we alternately estimate both components into the coarse-to-fine procedure. Accordingly, and knowing that χ = log(e ξ + e ξR -I 4⇥4 ), energy (5.64) is minimized by the following iterative, alternating estimation of ⇠ R and ⇠:

a. Given a fixed ⇠, solve for the twist motion ⇠ R that minimizes E Rig (χ).

b. Given a fixed ⇠ R , solve for the twist field ⇠ that minimizes E Res (χ).

Step a. is equivalent to minimize E Rig (⇠ R ) as is done in Section 5.4.1, but adding to the warp of every scene point x, the 3D motion offset given by (e ξ -I 4⇥4 ) X,which compensates the nonrigid component. This minimization is done by iteratively applying (5.38), with a zero auxiliary flow, i.e., cropping the last term of equations (5.39) and (5.40). Accordingly, the increment ∆⇠ R is computed as follows: (5.67) where gradient vector g R is given by:

∆⇠ R = -H -1 R g R ,
g R = X x2Ω Ψ 0 I,G ⇥ (I x J) T ⇢ I + γ(I g x J) T ⇢ g ⇤ + λ Ψ 0 Z (Z x J -K) T ⇢ Z ,
and the approximation of the Hessian matrix is given by:

H R = X x2Ω Ψ 0 I,G h (I x J) T (I x J)+γ (I g x J) T (I g x J) i + λ Ψ 0 Z (Z x J -K) T (Z x J -K).
On the other hand, step b. is equivalent to minimize E Res (⇠) as in Section 5.4, by iterating steps i and ii. The 3D motion offset, given by (e ξR -I 4⇥4 ) X, should be added into the warping function of point x, to compensate for the rigid component.

Algorithm 3 presents the general scheme of the scene flow estimation by modeling motion as a global rigid component plus a nonrigid residual.

Experiments

Middlebury datasets

The Middlebury stereo dataset [START_REF] Scharstein | High-accuracy stereo depth maps using structured light[END_REF] is commonly used as benchmark to compare scene flow methods [START_REF] Huguet | A variational method for scene flow estimation from stereo sequences[END_REF][START_REF] Basha | Multi-view scene flow estimation: A view centered variational approach[END_REF]Hadfield andBowden, 2011, 2014;[START_REF] Quiroga | Scene flow by tracking in intensity and depth data[END_REF][START_REF] Quiroga | Local/global scene flow estimation[END_REF]Quiroga et al., , 2014a;;[START_REF] Hornacek | Sphereflow: 6 dof scene flow from RGB-D pairs[END_REF]. Using images of one of these datasets is equivalent to a fixed camera observing an object moving in X direction. As in previous methods, we take images 2 and 6 as the first and second RGBD images, respectively. Gray images are obtained from the RGB images, while the ground truth disparity map of each image is used as depth channel, see Figure 6.1. These datasets present a simple 3D motion field, where the ground truth for the scene corresponds to the camera motion along the X axis, which coincides with the baseline of the stereo setup. For this reason, Y -and Z-components of the scene flow are zero while the X-component is constant. On the other hand, the optical flow is not constant and is given by the disparity map of the first image, as is shown in Figure 6.2. For the experiments we consider three Middlebury datasets: Teddy, Cones and Venus. Particularly, Teddy and Cones datasets present a large displacement optical flow, whose magnitude in pixels is in the range (12.5, 52.75) and (5.5, 55), respectively.

Error measures

Error metrics for optical flow

In order to compare with previous methods, the scene flow error is firstly measured in the image domain using the root mean squared error (RMSE) and the average angle error (AAE) of the optical flow. Let (u, v) and (u GT ,v GT ) denote the image flow vector and the ground-truth flow vector, respectively. The angular error (AE) is the angle in 3D space between (u, v, 1) and (u GT ,v GT , 1), which is computed as:

AE = cos -1 0 @ u • u GT + v • v GT +1 p u 2 + v 2 +1 q u 2 GT + v 2 GT +1 1 A .
Note that the AE corresponds to an angle in the range [0, ⇡]. Accordingly, the AEE is defined by:

AAE = 1 |Ω v | X x2Ωv AE(x), (6.1)
with |Ω v | the cardinality of set Ω v , which is defined as the subset of Ω containing the non-occluded pixels having a valid depth measure.

The endpoint error (EE) is defined as the Euclidean distance on the image between (u, v) and (u GT ,v GT ), and is given by: Using the EE, the RMSE of the optical flow is defined by:

EE = q (u -u GT ) 2 +(v -v GT ) 2 .
RMSE = s 1 |Ω v | X x2Ωv EE(x) 2 . (6.2)
Error metric for changes in depth

In order to measure how consistent is Z-component of the scene flow with the changes in depth, we use the root mean squared error of the Z-motion (RMSE Z ). This measure has been used to compare scene flow methods (Hadfield andBowden, 2011, 2014;[START_REF] Hornacek | Sphereflow: 6 dof scene flow from RGB-D pairs[END_REF]. Being u GT the ground-truth flow vector and v Z the Z-component of the 3D motion field v, and with {Z 1 ,Z 2 }, the first and second depth images, the RMSE-Z is defined by:

RMSE Z = s 1 |Ω v | X x2Ωv [(Z 2 (x + u GT ) -Z 1 (x)) -v Z (x)] 2 , (6.3)
Error metrics for scene flow

The computation of 2D motion errors does not necessarily reflect the errors in the scene since even small 2D displacements can produce large 3D motions. Moreover, it is preferable to compute 3D motions errors directly rather than measure independently errors of the optical flow and changes in depth, which can partially hide 3D motion errors. For instance, in Middlebury datasets the Z-component of the motion is zero, so that any scene flow estimation having no changes in depth will achieve a perfect score in RMSE Z . In order to measure the overall scene flow error we define the average normalized error of the scene flow (ANE V ). The point-wise scene flow error is defined as the Euclidean distance between the estimated 3D motion vector, v =( v X ,v Y ,v Z ), and the ground-truth 3D motion,

v GT =( v GT X ,v GT Y ,v GT Z )
. Subsequently, every error measure is normalized by the magnitude of the ground-truth scene flow, which is assumed to be non-zero, to present the error as a percentage of the motion magnitude. Accordingly, the ANE V is defined by:

ANE V = 1 |Ω v | X x2Ωv q (v X -v GT X ) 2 +(v Y -v GT Y ) 2 +(v Z -v GT Z ) 2 |v GT | . (6.4)
Particularly, for the Middlebury datasets the magnitude of the scene flow is constant and is given by the baseline of the stereo setup. Additional to the ANE V ,w e compute the statistic Rp% to show the percentage of pixels having a normalized scene flow error smaller or equal to p% of the ground truth magnitude.

Baseline methods

We consider two baseline methods that follows the proposed framework presented in Chapter 5: 6TwistFlow and 3DmotionFlow. We named 6TwistFlow the method that uses the twist-motion representation to minimize the scene flow energy (5.2).

The algorithm using the 3D-motion representation to minimize energy (5.49), is called 3DmotionFlow. The main parameters are listed in Table 6.1.

For both baseline methods we set the following parameters the same: L = 5, ↵ = N 2 β = 0, γ = 0, s 2 tr = max (off), M TV = 50, M GN = 5. Now, for the 6TwistFlow method we independently set: λ =0 .1, N = 5, M = 5, and (l)=N -2 10 -3 10 l . Finally, for the 3DmotionFlow we set: λ = 1, N = 11, M (l)={1, 2, 3, 4, 5, 5}, ⇣ = 0 (the local-rigid prior set off in E S ), N ζ = 5, and (l)= N -2 ⇥ {2 -2 , 10 -1 , 10 0 , 10 2 , 10 4 , 10 6 }. For both methods, the weighting function in the data term (5.12) follows a uniform distribution.

Comparison with other methods

We consider three groups of methods in the comparison: i ) optical flow methods, ii ) stereo-based methods and iii ) RGBD-based methods.

Optical flow methods

Optical flow methods use only the RGB information of first and second images to estimate the 2D motion. The scene flow is computed by inferring the 3D motion with the provided depth data, as is shown in Table 4.1 for the 2D-motion representation. We consider the large displacement optical flow algorithm by Brox and Malik (2011) (OpticalFlow1 ) (OpticalFlow2 ). Comparison with optical flow methods allow to asses the relative contribution of using depth data in addition.

Stereo-based methods

Stereo

Results

Table 6.2 presents the results on the Middlebury datasets for the eight algorithms considered. The proposed algorithm 6TwistFlow achieves the best results for RMSE and RMSE Z , confirming the benefit brought by the twist representation. Observe that in most cases this method outperforms SceneFlow2, which also use a field of rigid motion, but in a discrete optimization scheme, which could benefit from the propagation of the correct rigid motion. AEE results of both methods are very similar, differing only in about 1/4 of degree. On the other hand, the alternative approach 3DmotionFlow outperforms most of the considered methods, except SceneFlow2. Particularly, it is interesting to remark that 3DmotionFlow outperforms SceneFlow2, which directly represent the scene in 3D and solve for the scene flow testing and propagating 3D motion hypothesis. Although the 3D motion of every scene point is the same and few correct hypothesis could favor the final solution through the particle filtering, the incremental 3D motion estimation done by 3DmotionFlow performs better. Moreover, the 2D parametrization of the scene, used both by the 6TwistFlow and by the 3DmotionFlow algorithms, allows to accurately solve for the 3D motion field. It is important to remark that SceneFlow1 and SceneFlow2 use a 3D representation of the scene, as a 3D point cloud and a set of 3D patches, respectively. Therefore, there is no any advantage of using a 3D representation for the motion estimation on these datasets.

Regarding stereo-based methods, the measured errors are higher in most cases since the depth information is estimated at the same time. However, thanks to the two additional views these methods benefits from smaller occluded regions, and the constraints brought by the stereo setup favor the estimation of the purely horizontal 2D motion. Finally, optical flow methods are the worst performer, which is expected since they only use the brightness of two views. Particularly, Teddy and Cones datasets are challenging motion estimation problems because the observed optical flow presents a wide range, around 40 and 50 pixels, respectively. Unlike their RMSE is between the highest errors, for both methods, the AEE is very low. This accurate angular estimation is explained by the strong TV regularization done directly on the optical flow and because all pixels move in the same direction. Finally, it is observed from Table 6.2 that the critical measure for optical flow methods is the estimation of changes in depth, although there is no Z-motion in these datasets. Due to the lack of depth information in the motion estimation, even small errors on the image motion yield to large error in the 3D motion inferred using the provided depth data. Accordingly, the jointly use of brightness and depth data allows to a more accurate estimation of motion on the image domain, and simultaneously, changes in depth.

Analysis of the scene flow framework

The proposed scene flow framework allows the implementation of different approaches through the setting of parameters. For example, local and piecewise rigidity can be adjusted by varying the size of the local neighborhood, N (x), and the weight of the smoothness term, ↵, respectively. Similarly, it is possible to configure the balance between several terms, such as brightness and gradient constraints, brightness-and depth-based constraints, and data en smoothness energies. We analyze below how the configuration of some of these parameters affects the performance of the scene flow framework on Teddy and Cones datasets. Additional to RMSE, AAE and RMSE Z , we compute ANE V and R5% to directly compare 3D motion errors.

Local rigidity

The scene flow is assumed to be semi-rigid, i.e., local and piecewise rigid. While 6TwistFlow directly encourages local rigidity in data term, 3DmotionFlow allows to incorporate the local rigidity in both terms: in the data term it is done by encouraging a constant 3D translation, while in the smoothness term a local rigid motion is fitted and the residual motion minimized. Using baseline methods, we perform an experiment varying the size of N (x), thereby controlling the amount of local rigidity assumed in the data term. We start from N (x)={x}, i.e., no local rigidity, to N (x) being a 11 ⇥ 11 window, as is presented in Table 6.3. In the 3Dmotion case, the bigger the neighborhood the more accurate the results, showing that a strong local rigidity encouragement favors the scene flow estimation. On the other hand, the best performance is obtained for middle size windows, 5 ⇥ 5 and 7 ⇥ 7. Although RMSE subtly decreases as the size of the window increases, ANE v and AAE present an inflection point around 5 ⇥ 5, as is shown in Figure 6.3.T h i si s a meaningful result for 6TwistFlow, showing that a bigger local neighborhood does not imply a better estimation of the 3D motion, even for this constant scene flow. Moreover, unlike 3DmotionFlow, accurate results for 6TwistFlow start from with a small window, 3 ⇥ 3, therefore it can be inferred that the amount of local rigidity assumption has less significance in the final results.

Particularly for 3DmotionFlow, local rigidity can be encouraged in data and smoothness terms. In order to asses the benefit of each strategy, we consider Table 6.3: Variation of the local rigidity in the data term. Top results in red and blue for 3DmotionFlow and 6TwistFlow, respectively. Baseline methods are highlighted in green.

0" 0,5"

1" 3" 5" 7" 9" 11"

6TwistFlow:+errors+vs.+N"

ANEz" RMSE" AAE" 0" 0,5" two variations of the baseline method: 3DmotionE S and 3DmotionALL.T h efi r s t , 3DmotionE S , removes the local encouragement in the data term by setting N (x)= x, i.e., one pixel neighborhood, and includes the local rigid prior in the smoothness term. We modify the linking parameter to improve the accuracy by setting (l)={10 -2 , 10 -1 , 10 0 , 10 0 , 10 1 , 10 1 }, giving more weight to the regularization in order to improve the stability in the minimization of the data term. The second approach, 3DmotionALL, applies the local-rigid prior both in data and smoothness term, which is equivalent to add the local-rigid prior in the regularization of the baseline method. In both cases, the local fitting is done on a 11 ⇥ 11 neighborhood, i.e., N ζ = 11, and the local-rigid parameter is set as ⇣ = 3. The use of the localrigid prior in addition to the local estimation in the data term, yields to a more accurate estimation of the motion, as is presented in Table 6.4. On the other hand, the encouragement alone of the local motion seems to be more effective in the data term than in the regularization procedure.

Piecewise rigidity

Piecewise rigidity is encouraged in the smoothness term by TV regularization. We investigate the influence of this assumption in the performance of 6TwistFlow and 3DmotionFlow, by varying the number of iterations of the TV solver, M TV . As is shown in Table 6.5, a stronger piecewise rigidity yields to a better results since the global rigid motion of the dataset. Also observe that for M TV = 0, there is no TV regularization favoring piecewise rigidity but there exists a Tikhonov regularization on the energies from the linking term controlled by . This kind of regularization favors small norm solutions. In this case, unlike 6TwistFlow that provides totally wrong estimations, the performance of 3DmotionFlow is still acceptable. This is explained by the median filtering on 3DmotionF low, which is directly related with l 1 minimization [START_REF] Sun | Secrets of optical flow estimation and their principles[END_REF] and therefore favor piecewise smooth solutions.

Motion representation

This experiment compares the different motion representation. Alternatively to 6TwistFlow and 3DmotionFlow, it is also possible to use a simpler motion representation considering a orthographic camera, as in most optical flow methods. As is explained in Section 4.3.4, in this case, the depth data is no directly used to model the image motion. The proposed framework can be straightforwardly used with the 2D-motion representation. We named 2DmotionFlow the method based on energy (5.49), with the warping function given by (4.27). This baseline method is set with

L = 5, ↵ = N 2 , β = 0, γ = 0, s 2 tr = max (off) λ =0 .1, N = 7, M = 5, M GN = 5, M TV = 100, ⇣ = 0, N ζ = 5 and (l)=N -2 ⇥ {2, 10 1 ,
20, 10 2 , 10 3 , 10 3 }. Table 6.6 shows errors for the three representations, and includes GroundTruth, which is obtained by applying the ground truth 3D motion given by the baseline of the stereo setup. Note that for GroundTruth the RMSE is not zero due to the dataset disparities, which are used to compute the ground truth optical flow, have quarter-pixel accuracy. As it was concluded from previous experiments, the method based on the twist representation performs better that the method using the 3D representation. Observe that the use of the 2D-motion representation does not bring any advantage over optical flow methods for the estimation of the image motion, see Table 6.2. For this reason, the use of depth data only as an additional brightness channel, wastes useful information for the constraint of the image motion. On the other, the estimation of changes in depth by 2DmotionFlow outperforms optical flow methods; disparity errors are 0.06 and 0.12, for Teddy and Cones, respectively. This is thanks to the use of a depth constancy assumption, which allows to regularize the estimation of the Z-component of the scene flow. Conversely, optical flow methods only infer the 3D motion field from the provided depth data. Finally, optical flow results are presented in Figures 6.4 and 6.5.

Brightness and depth balance

The use of a warping function allows to constrain the scene flow on the image and therefore to exploit both intensity and depth data. We experiment by varying the balancing term λ to modify the contribution of brightness and data terms in the scene flow energy. We consider λ = {0, 0.1, 1, 1}. Notation λ = 1 stands for the case where only depth constancy constraints are used in the data term. Table 6.7 presents the results. It can be noted that the simultaneous constraint of the motion in brightness and depth achieves the best performance for each representation. Particularly, for the twist representation there is subtle difference between the baseline method and the constraint only on the brightness image. This can be explained by the zero Z-component of the scene flow and by the ability of the twist representation to constrain this simple motion on the brightness image. However, the same difference becomes larger for the 3D-motion representation. On the other hand, for the 2D-motion representation the addition of depth constraints brings a very slight advantages respect of the optical flow case (λ = 0), even this latter performs better in some cases. The constraint of the motion only on depth data is the worst performer in each case, since the texture of depth images is by far less meaningful that such of brightness images. Finally, it is important to remark that the best balance varies from one motion representation to another, and therefore it is no possible to conclude about the optimal combination.

Processing time and resolution

We measured the runtime of non-optimized implementations of baseline methods 3DmotionFlow and 6TwistFlow, using the single-core 64-bit Intel Xeon E5520 @ 2.27GHz with 8 Gb of RAM. Table 6.8 summarizes the computational times of both approaches, and includes the runtimes reported of SceneFlow1 [START_REF] Hadfield | Scene particles: Unregularized particle-based scene flow estimation[END_REF], which also uses a single-core machine. Proposed methods 3DmotionFlow and 6TwistFlow requires less than 1/10 and 1/5 of runtime than SceneFlow1, respectively, and both achieve better estimations. Also observe that thanks to its simpler parametrization, 3DmotionFlow runs almost twice faster than 6TwistFlow. Comparison with SceneFlow2 was not possible because its runtime is not reported [START_REF] Hornacek | Sphereflow: 6 dof scene flow from RGB-D pairs[END_REF], but authors claim it takes about 300 s. Computational speed is not our primary concern but there are several ways to speed up the scene flow methods. Particularly, we focus in a simple modification using our implementations in the same machine. For this purpose, we modify the resolution of the scene representation to get faster estimations, and at the same time, we measure the accuracy of the obtained estimate. First, we experiment with an early stop of the coarse-to-fine procedure as follows. Being at level l = PyrStop of the RGBD pyramid, the estimation is stopped and the final motion field at the original resolution is obtained by linear interpolation. Table 6.9 presents the average results on Teddy and Cones datasets. Observe that with PyrStop = 1, the computation is stopped at the first octave (half of the initial scale) and runtimes are reduced by about 2 and 4, for 3DmotionFlow and 6TwistFlow, respectively, while the accuracy is almost equivalent to that obtained at the original resolution.

Finally, we consider only every mth pixel of the motion field for the Gauss-Newton procedure in the data term. The result of each operation is copied to the surrounding pixels. For example, setting m(l)={4, 4, 2, 2, 1, 1} as a function of the pyramid level l, for 3DmotionFlow, the runtime is reduced to 5.9 s, i.e., about 7 times faster, while the accuracy stands almost unchanged. Similarly, setting m(l)={6, 6, 4, 4, 2, 2} for 6TwistFlow, the runtime becomes 13.8 s, i.e., about 6 times faster, with nearly the same accuracy. Table 6.9: Early stopping in the coarse-to-fine procedure. PyrStop stands for the level of the pyramid where the computation is stopped. The resulting motion field is linearly interpolated to compute the final motion field.

Scene flow from RGBD images

We performed further experiments on more complex scenes using two RGBD sensors: the Microsoft Kinect for Xbox and the Asus Xtion Pro Live. We consider three different setups: i ) a fixed or moving camera assuming a nonrigid motion, ii )a fixed or moving camera observing a rigid scene and iii ) a moving camera capturing deformable objects. In each case, we show the input images, the optical flow, and one or more components of the scene flow. Also, depending of the experiment we present the warped images, the resulting brightness and depth residuals and the 3D reconstruction. The following example presents a scene flow estimation example.

A fixed camera is observing a non rigid scene, as is showed in Figure 6.6.T h e scene is segmented in the image domain using the observed depth measures of the first depth image. For this example, the estimation is done only for scene points within 50 cm and 200 cm from the RGBD sensor. The bottom images in Figure 6.6 give an idea of the segmentation. It is assumed that the 3D surfaces of interest are within this given range. No motion estimation is done for the remaining points. Nevertheless, all points of the second image, having a valid depth measure, are considered for the scene flow estimation. Results are presented as is shown Figure 6.7. The optical flow is visualized using the Middlebury color code [START_REF] Scharstein | High-accuracy stereo depth maps using structured light[END_REF]. For the scene flow, we show each component using a coldto-warm code, where green color is zero motion, and warmer and colder colors represent positive and negative velocities, respectively. Finally, Figure 6.8 shows the brightness and depth residuals, which correspond to the differences between the original and warped RGBD representations of the scene. The bottom row presents the warped image and differences with the estimated scene flow. Brightness differences are scaled by 8 for visualization, and depth differences go from black to white, for 0 cm to 1 cm, and which for larger differences. Cyan color stands for an invalid depth difference.

Nonrigid motion estimation

In this section we perform some nonrigid estimation experiments. Firstly, we consider a large Z-motion case presented by [START_REF] Hornacek | Sphereflow: 6 dof scene flow from RGB-D pairs[END_REF], where it is possible to compare some of the current scene flow methods. Secondly, four further experiments are performed using the proposed approaches.

Large Z-motion experiment [START_REF] Hornacek | Sphereflow: 6 dof scene flow from RGB-D pairs[END_REF] propose to estimate in both directions, the motion represented by the pair of RGBD images shown in Figure 6.9. This is a challenging setup due to the large displacement in Z-direction, illumination changes, partial occlusions and the low resolution of the RGB images, 320 ⇥ 240 pixels. Results of four scene flow methods are presented in Figure 6.10, including our approach [START_REF] Quiroga | Local/global scene flow estimation[END_REF], which corresponds to the baseline method 3DmotionFlow in the previous section. Observe that [START_REF] Herbst | RGB-D flow: Dense 3-D motion estimation using color and depth[END_REF] and [START_REF] Hadfield | Scene particles: Unregularized particle-based scene flow estimation[END_REF] fail capturing the coarse motion, while [START_REF] Hornacek | Sphereflow: 6 dof scene flow from RGB-D pairs[END_REF] and [START_REF] Quiroga | Local/global scene flow estimation[END_REF] are able to give a rough idea of the observed motion. However, it is not possible to measure the accuracy using only the estimated optical flow. We replicated the experiments using the methods 6TwistFlow and 3DmotionALL, and also include the estimation using the large displacement method (LDOF) by [START_REF] Brox | Large displacement optical flow: Descriptor matching in variational motion estimation[END_REF], to have an idea of the performance of current optical flow methods. Results are presented in Figure 6.11, where it can be seen that all methods are able to capture the rough motion with some small differences. In order to have an idea of the accuracy, we warp second images using the estimated motion. The resulting warped images are shown in Figure 6.12. It can be seen that 6TwistFlow models the best the large Z-motion performed by the person on the right, especially visible in the warped texture of the face. However, there exist some problems estimating the motion of one of the hands in the second case, due to the low resolution of depth data and the splitting of the surface from frame 2 to 1. On the other hand, 3DmotionALL is able to model the motion of both hands in both cases but it deforms the texture of the face due to the strong regularization that domains the estimation. The optical flow method only gives a roughly estimation of the observed motion. Note that none of the methods is able to estimate the head motion of the woman on the right, which presents a large rotation between both frames. Figure 6.13 shows the Z-component of the scene flow, where is evident the benefit of using the depth constraints in conjunction with a projective warp that models the image motion. The optical flow method is unable to determine the motion boundaries and suffers from the depth data noise, since changes in depth are directly estimated using the raw depth data.

Hand rotation

This experiment presents a hand rotating in the middle of the scene, moving quasirigidly, as is shown in Figure 6.14. The depth range is set to include only the surface that belongs to the arm and forearm. Figure 6.15 presents the resulting optical flow and two components of the scene flow, for each method. Regarding the optical flow, it can be noted that 3DmotionALL subtly distorts the motion field due to smoothness performed on the 3D motion field. The distortion of the motion is more evident for the optical flow method, since the regularization is done on the 2D motion field. On the other hand, thanks to the regularization done on the twist motion field, 6TwistFlow performs the most accurate estimation.

The advantages of the proposed methods are more evident when observing the components of the scene flow. Unlike LDOF that estimates the X-component of the 3D motion as curve levels, 6TwistFlow and 3DmotionALL are able to accurately reflect the smooth transition brought by the 3D rotation. Similarly, the estimation of the Z-component from both methods, reveal the smooth changes in depth of the 3D surface while rotating. Contrary to the results obtained by using the optical flow estimate, which is affected by the 2D errors and the noise in depth data.

Two arms rotation

In this experiment a fixed camera observes a motion performed with arms and hands. While hands are rotating inwards, the elbows lift, as is shown in Figure 6.16. This composite motion generates a varied optical flow, which is well estimated by our two methods, as is presented in the top images of Figure 6.17. Moreover, it can be seen that small rotations and articulated motions are well described. Observe that the optical flow method provides a smoother solution, showing be less sensitive to the specular highlights appearing on the skin. This is because its simpler 2D motion model, which is stronger regularized on the image. A higher parametrization seems to be more sensitive to this specular effect and a stronger regularization may be required. However, as is observed in the bottom images of Figure 6.17,t h i s smoother optical flow does not imply a more accurate estimate of the scene flow. The 6TwistFlow method performs the best estimation of changes in depth, being able to model the motion of thumbs, palm hands and elbows.

Poster stretching

This experiment presents the stretching of a poster, as is shown in Figure 6.18. This is a challenging motion for the semi-rigid assumption due to the non-uniform deformation of the 3D surface. Figure 6.19 presents the resulting optical flow and Z-component of the scene flow. The two proposed methods are able to capture the poster deformation, thus it is possible to accurately estimate the changes in depth when the poster is folded. The gradient constancy constraint plays an important role here, since the sensor applies an automatic white balancing. Observe that the smoothness of the optical flow, inside the poster, varies according to the motion representation where the regularization is applied. Also, as can be seen on the left side of the poster, the partial occlusion brought by the stretching is estimated in a different way for each approach. Particularly, 6TwistFlow and 3DmotionALL suffer from the misalignment between the poster edges in brightness and depth data, and from the depth uncertainty near to the object boundaries. On other hand, LDOF estimates a sharp motion discontinuity since it does not have to explain the ambiguous observed data and thanks to the good contrast between the poster and the background. However, depth changes are better captured by our two methods, as is shown in the bottom images of Figure 6.19. Particularly, 3DmotionALL achieves the best performance in this experiments, accurately estimating the Zmotion and dealing in a better way with the partial occlusion, as is shown in the resulting warped images in Figure 6.20. Note that we do not use explicit occlusion awareness, but the robust norm is able to reduce its effect. Fig. 6.20: Warped images for the poster stretching experiment. Unlike LDOF, 6TwistFlow and 3DmotionALL are able to correctly warp the poster. However, observe that these two methods expand the shirt sleeve close the poster, due to the partial occlusion. In order to reduce these artifacts is required to include an explicit occlusion awareness in the estimation.

Teddy bear

This experiment considers a pair of images of the RGBD dataset [START_REF] Sturm | A benchmark for the evaluation of RGB-D slam systems[END_REF], where a moving camera observes a rigid scene. Figure 6.21 shows the input images. Unlike Middlebury stereo datasets, the observed motion in this experiment has a nonzero rotation component and therefore the scene flow is not constant. Optical flow results are presented in Figure 6.22, and as is expected because its motion representation, 6TwistFlow performs the most accurate estimation.

Rigid motion estimation

The estimation of a single rigid motion is a particular case of the proposed framework (5.2) and can be done by minimizing the rigid energy defined by 5.65. Using this formulation is possible to solve for the rigid motion that best explains both intensity and depth. For example, for the camera motion shown by the images in Figure 6.23,w es e tw I (x) = 1 and w Z (x)=λ =0 .1, to obtain the optical flow and Z-motion presented in Figure 6.24. In this case, parameter λ allows to control the balance between brightness and depth constraints in the rigid energy. Also, it is possible to use only brightness or depth constraints, by setting λ = 0 or λ = 1, respectively. A comparison of the final residuals of brightness and depth, for different values of λ, is presented in Figure 6.25.

3D reconstruction

The optimal balance between brightness and depth is an open question and should be set according to the confidence in the observed measures. Also, this balance can be defined point-wise, as we proposed in our rigid energy. Throughout experimentation we found out that setting λ =0 .1, for all pixels, allows accurate estimation using the Kinect and Xtion sensors. In order to illustrate the precision of the rigid estimation we perform a 3D reconstruction using both sensors. Getting a 3D model from a single depth image is not accurate due to noise and occlusions, but adding a couple of warped views is possible to improve the quality of the reconstruction. Accordingly, to complete the 3D model, we choose a reference frame and warp a number of close depth images using the rigid estimations. The reconstruction is done using the Volumetric Range Image Processing Package (VripPack) by [START_REF] Curless | A volumetric method for building complex models from range images[END_REF], which merges range images into a compressed volumetric grid. Experiments using images from the Xtion and Kinect sensor are presented in Figures 6.26 and 6.27, respectively. Finally, we compare the quality of the reconstruction while varying the balance between brightness and depth. In this experiment a fixed camera observes a person seated on a swivel chair, whom is requested for rotating in front of the camera, as rigid as possible.

The reference frames and limit images are shown in Figure 6.28, giving an idea of the visible surface. Observe that the considered images belong to a large rotation and therefore some parts of the reference surface result partially occluded. Moreover, the illumination is not uniform and there exist specular highlights, especially on the skin face. In Figure 6.29 the reconstructions for different values of λ are compared.

In this challenging experiment, the best results are achieved using only the depth constraints. The jointly use of brightness and depth constraints, outperforms the only-brightness estimation but is less accurate than the results by the only-depth approach.

Nonrigid plus rigid motion estimation

In many applications, the sensor itself moves relative to the observed scene and compensating for the motion of the camera can simplify the estimation and regularization of the scene flow, as is shown in Figure 5.7 . Moreover, for 3D reconstruction of deformable objects, the camera motion is needed to register partial 3D reconstructions. Using the energy 5.64 we are able to split the motion of the scene into a globally rigid component, ⇠ R =( ! R , ⌧ R ) 2 R 6 , capturing the camera motion relative to the dominant object/background, plus a nonrigid residual field, ⇠ =(!, ⌧ ), as we demonstrate in the two following two examples.

In the first experiment, we use the Xtion sensor to capture the motion of head and hands performed by a person, as is shown Figure 6.30. We compute the relative motion from a given reference frame, to a previous and a posterior frames. Components of the estimated motion are presented in Figure 6.31. Particularly, the motion of head and right hand are well described by the nonrigid component. On the other hand, the estimation of the left hand motion is a challenge, especially from the reference frame to the first frame, due to the sudden deformation and partial occlusion. Using the estimated scene flow, the first and last frames can be warped to the reference frame, as we done in Figure 6.32. Observe that despite the limited resolution of the RGBD images, it is possible to register the observations even under this nonrigid motion. This can be useful for constructing partial 3D models of deformable objects using a moving camera. Moreover, occlusions across Observe that the rigid motion correctly compensates for the body motion. The nonrigid estimation is able to roughly register the rotating head, as is noted by observing the right location of glasses, the beard and the ear. The motion estimation of the left hand is not well modeled.

views can be estimated observing the depth and brightness residuals, as we illustrate in the following experiment.

In the second experiment, a Kinect sensor is subtly moved to observe the motion of the head of a person, as is shown in Figure 6.33. Components of the estimated motion are presented in Figure 6.34, where can be observed that the motion of the camera and head are well estimated. Figure 6.35 presents the resulting warped images using the estimated scene flow. Note that the head rotation is correctly compensated except for the occluded regions. By observing the brightness and depth residuals is possible to determine which regions of the reference frame are occluded in a given view, as we present in Figure 6.36. Particularly, the depth data provides confident measures to detect occlusion, unlike the brightness images, which could suffer from the specular highlights and illumination changes. Brightness and depth residuals are presented in top and bottom images, respectively. Observe that high and undetermined residuals can be used to reliably estimate the partial occlusions. On the other hand, occlusions are also visible in the brightness residual, but can be confusable with the specular highlights which also present high residual.

Conclusion

There have been two questions guiding the development of this thesis. The first question concerns the usage of both sources of data: "How to fully exploit color and depth to compute a reliably scene flow?" The second question addresses the way the motion is parametrized: "Which motion representation should be used to compute a confident scene flow?" These two questions defined the starting point of our exploration and inspired us to develop this semi-rigid framework to compute dense scene flow from RGBD images. Below we summarize the main contributions of this thesis and conclude with our vision of researching perspectives on this field.

Contributions

Color and depth exploitation 1. A warping function that allows to constrain the scene flow on RGBD images is presented. With help from the depth channel, this warp models the image motion as a fully projective function of the 3D motion. This way is possible to define consistency constraints directly on RGBD images.

2. In order to take advantage of the depth channel, we define constraints directly on the depth image, by measuring how consistent is the scene flow estimate with the observed 3D surface.

3. We simultaneously constrain the scene flow to be consistent with brightness and depth data. This allows us to exploit both sources of measurements, which are uncorrelated by nature. This way we are able to solve for the scene flow that best explains the RGBD observations.

4. The proposed framework is flexible enough to support different data terms. The set of constraints can be directly expanded, or modified, to work with color components in any color space, and also to include constraints from invariants properties of 3D surfaces, on the depth data.

5. We have intro duced different strategies to benefit from the depth channel for scene flow estimation. Depth data has been used to: i ) model the image motion via the warping function, ii ) directly constrain the scene flow in depth images, iii ) guide the TV regularization, preventing smoothness across surface edges, and iv ) estimate a prior of the local rigidity.

Motion representation

1. A semi-rigid framework for dense scene flow estimation is presented, taking advantage of the local and piecewise rigidity of most real scenes. The proposed framework supports different motion representation, data terms and regularization strategies. The evaluation demonstrates that the proposed method achieves the best results in the most commonly used scene flow benchmark.

2. We introduce an over-parametrization of the scene flow by using a field of rigid motions, with a twist motion representation. This way semi-rigid properties of the scene are directly exploited. This motion representation is provided with a warping function, allowing the formulation of RGBD constraints on the image plane.

3. Using the twist-motion representation, a semi-rigid scene flow energy is defined. Locally-rigid motions are encouraged in the data term, while piecewise rigid solutions are favor in the smoothness term. The framework provides an adjustable combination between local and piecewise rigidity.

4. We also present an alternative scene flow energy by directly using a 3D-motion representation, which is provided of a corresponding warping function. In this case, we show that the locally rigidity can be exploited in both data and smoothness terms.

5. We derive a linear version of the warp for each motion representation, enabling the solution of the scene flow by an incremental minimization. We fully describe the optimization of the two proposed energies.

6. By using the same general framework, we model the scene flow as a global rigid motion plus a nonrigid residual. This is particularly useful when estimating the motion of deformable objects in conjunction with a moving camera.

7. We have performed different experiments to evaluate some of the possibilities of the proposed framework, such as the motion representation, the brightness and depth terms balance, and the local and piecewise rigidity assumptions.

Through additional experiments we indicate the general applicability of our approach in a variety of different scenarios.

Open questions

The study done in this thesis moves forward the state of the art in scene flow by formulating a general framework for scene flow estimation. This is only a first step in a long journey of 3D motion estimation from RGBD images and our attempt opens the door to future directions for exploration.

Accuracy

Nowadays there is not an appropriate, available benchmark for scene flow. Because the 3D motion field resulting from camera translation is constant, Middlebury stereo datasets are not well suited to fully evaluate the performance of scene flow methods. Although the quality of the results on the sequences from RGBD sensors is motivating, and important remarks are done, it was not possible to measure the 3D errors for interesting nonrigid motions. Therefore, the optimal setting of parameters or the right balance between assumptions is still a query, as well as the best configuration of the proposed approach for a given specific task. For example, it is not clear the correct balance between brightness and depth terms, which could be controlling by defining pointwise confidence measures of the observed data. Moreover, further exploration is required on the appropriate amount of local and piecewise rigidity assumptions. The availability of such benchmark will allow the comparison with current methods and support the development of the new ones.

It is expected that for the MPI Sintel flow dataset [START_REF] Butler | A naturalistic open source movie for optical flow evaluation[END_REF], which is derived from the open source 3D animated short film Sintel, the ground truth of the depth channel be released and can be used as scene flow benchmark. Unlike Middlebury images, this dataset contains long sequences and nonrigid motions.

Depth data

Current RGBD sensors provide depth measures with a reasonable precision. For the Kinect for Xbox errors are smaller than 2mm at 1m from sensor [START_REF] Khoshelham | Accuracy and resolution of kinect depth data for indoor mapping applications[END_REF], but increase quadratically with depth to about 4cm at the maximum distance of 5m. We have used the Charbonnier penalty to deal with the error in depth measures. However, a better modeling could be beneficial since errors in depth have a two-side influence in our approach: on the image motion modeling and on depth constancy constraints. Moreover, the quality of the registration between color and depth data is a function of errors in the depth measurement.

For this reason, in some cases is evident the misalignment between color and depth, especially close to the borders where there exists a large uncertainty on depth measures. We have used the automatic registration given by each sensor, which is based on the factory parameters, but also is possible to calibrate each sensor to improve the registration [START_REF] Herrera | Accurate and practical calibration of a depth and color camera pair[END_REF]. Also, a preprocessing step can improve the quality of depth channel. For example, a hole-filling algorithm and a bilateral filtering jointly color images to enhace and correct the edges of the 3D surface [START_REF] Richardt | Coherent spatiotemporal filtering, upsampling and rendering of rgbz videos[END_REF].

It is important to remark that the precision, resolution and working range of depth sensors have improved in last years and is possible that some of these current issues on depth data fade away soon.

Large displacements

How to deal with large displacement is still an open question. In the rigid case, as in the Middlebury dataset, the semi-rigid assumption benefit from the constant scene flow. Nevertheless, in a nonrigid case, large displacements are harder to estimate. If the motion of smaller structures is similar to the motion of larger ones, the coarseto-fine approach works well and improves global convergence. However, if the scale of a structure is smaller than its displacement or its motion is no coherent with larger structures, the scene flow is not well estimated. Larger structures dominate the estimation in lower scales, and then local minima in higher scales prevent the correct estimation. It is possible to include a set of sparse 2D matches in the variational formulation [START_REF] Brox | Large displacement optical flow: Descriptor matching in variational motion estimation[END_REF], e.g., using SIFT or SURF features. We follow this approach in [START_REF] Quiroga | Local/global scene flow estimation[END_REF], where the scene flow is skewed to satisfy a set of 3D motion hypothesis to deal with large displacements. Nevertheless, the low quality of brightness images and the lack of texture, as presented on the skin and uniform clothes, make the matching of 2D features a challenge. Also, most of these features are located at edges of the brightness image, coinciding with edges of the 3D surface where depth measures are not reliable or available, and therefore no valid 3D motion hypothesis. For this reason, a further exploration has to be done. An interesting alternative can be to combine the proposed framework with a discrete minimization procedure to avoid get trapped in poor local optima, as is done by [START_REF] Lempitsky | Fusionflow: Discrete-continuous optimization for optical flow estimation[END_REF] using graph cuts. Also, it would possible to directly perturb the field of rigid motions to deal with larger motions, as is proposed by [START_REF] Hornacek | Sphereflow: 6 dof scene flow from RGB-D pairs[END_REF]. Finally, according to the specific application, a priori information of the scene of the deformable object can be exploited to achieve reliable estimations.

Regularization of the motion field

Total variation is of the most popular regularizers for variational motion estimation because its nice properties on discontinuities. However, its extension to functions with values in a manifold is an open problem. This is precisely the regularization case on the field of rigid motions. We have decoupled the regularization procedure for the rotational and translational part and proposed some approximations to simplify the optimization. The decoupled regularization may not be the optimal choice from the modeling point of view since the discontinuities of both components are expected to coincide. Moreover, it is only an approximation of the real structure of the manifold. However, the decoupling allows for a simple and effective solution, and also allows the framework to support other regularization strategies. Recently, [START_REF] Lellmann | Total variation regularization for functions with values in a manifold[END_REF] reformulate TV as a multi-label optimization problem with an infinite number of labels. This hard optimization problem is approximately solved using convex relaxation and is applied to regularize three-dimensional rotations. Future advances on manifold regularization may provide even more accurate and faster solvers that can be used with our parameterization. On the other hand, the regularization of the 3D motion field can be exactly formulated using TV. Nevertheless, there exist different alternatives that can be adopted and a further exploration is required. Results using channel-by-channel and vectorial TV are very similar in our experimentation. Also, it is possible to use alternative norms, as the Huber norm, or consider the total generalized variation [START_REF] Bredies | Total generalized variation[END_REF]. An appropriate scene flow benchmark is required to analyze in detail advantages and differences among these strategies.

Future work

Temporal consistency

Most scene flow applications use RGBD sequences as input and therefore temporal information can be exploited. Enforcing temporal consistency has proved to be useful for optical flow estimation [START_REF] Volz | Modeling temporal coherence for optical flow[END_REF], and also recently, it has been used for scene flow estimation from stereo [START_REF] Vogel | View-consistent 3d scene flow estimation over multiple frames[END_REF], achieving stateof-the-art results. Particularly, 3D motion estimation can benefit from temporal information thanks to the reluctance of objects to change their way of moving (inertia). For this reason, using a sufficiently high frame rate, real scenes present temporal consistency over more than two frames due to the inertia of moving surfaces. Accordingly, locally rigid motions can be assumed to be approximately constant on a short time interval. Also, for the 3D motion representation, smooth 3D trajectories can be encouraged in the solution. Using more than two frames yields to a more robust formulation, where the influence of noise, occlusion and missing data can be reduced.

Real-time implementation

Computational speed was not our primary concern in this work, but there are several ways to speed up a scene flow method that follows our framework. First of all, it is highly parallelizable since computations are done independently per pixel and only considering few nearby local values. Moreover, the local rigidity is configurable and can be even reduced to one pixel to accelerate the method. Also, the image warping can be applied every n iterations instead of every iteration, an early stopping criterion can be used in the Gauss-Newton and TV solvers. Finally, depending on the application, a subsampled version of the field of rigid motions can be computed, i.e., considering only every nth pixel in the motion field, or alternatively stopping the computation in a coarser scale in the RGBD pyramid. Time benefit of some of these strategies is analyzed during our experimentation.

Scene flow descriptors

There is no work that directly computes scene flow to perform tasks such as action recognition, gesture classification or interaction. Probably, this is due to the fact that most existing methods require fully calibrated stereo or multi-view camera systems, which are not always available. Besides, most of these methods require a lot of processing time, becoming not suitable for real time applications, and their performance is limited facing nonrigid motions, which are frequent in such tasks. On the other hand, the optical flow, which is related with the scene flow projection on the image, has been successfully used in motion analysis. For example, histograms of optical flow are commonly used in state-of-the-art techniques in action recognition to construct descriptors over spatio-temporal interest points [START_REF] Laptev | Space-time interest points[END_REF][START_REF] Niebles | Modelling temporal structure of decomposable motion segments for activity classification[END_REF] and to extract 2D trajectories by tracking key-points [START_REF] Messing | Activity recognition using the velocity histories of tracked keypoints[END_REF][START_REF] Matikainen | Trajectons: Action recognition through the motion analysis of tracked features[END_REF]. Furthermore, since trajectory based methods outperform other state-of-the-art approaches for action classification [START_REF] Wang | Action recognition by dense trajectories[END_REF], it is promising to use scene flow to capture motion information by extracting accurate 3D trajectories, which can be used to model in a better way a wider variety of motions.

Using the proposed framework we are able to solve for an accurate dense scene flow, but is also possible to solve for a set of sparse 3D trajectories, both in a more precise way than previous methods. Currently methods performing motion analysis from RBGB data, such as the pose estimation approach by [START_REF] Shotton | Real-time human pose recognition in parts from a single depth image[END_REF], are based on a human body model, which requires a huge database for training and are constrained to the specific application. Moreover, color information is not simultaneously exploited. We differ from such approaches and aim to compute a model-less scene flow. Since we do not use any restriction on the shape of moving objects or the type of movements, our approach can be used in different tasks where 3D motion information is useful. The definition of appropriate scene flow descriptors, encoding the estimated 3D motion field, can yield better and novel applications on recognition, classification and interaction using RGBD sensors. For instance, the distribution of rigid motions, histograms of 3D trajectories, displacements or orientation, can be used as descriptors.

3D reconstruction of deformable objects

Using the proposed semi-rigid energy is possible to model the motion of the scene as a rigid component plus a non-rigid residual. This motion split is particularly useful for modeling non-rigid objects. Usually, the object of interest is captured using a moving camera, in order to observe regions that are not visible from a single view. Accordingly, the camera is moved around the 3D surface, capturing new observations to complement the desired model. However, the deformable object itself moves in a non-rigid way, and this motion has to be compensated to correctly create the 3D model. Note that this deformation of the object makes the registration task a harder problem, but at the same time, can be useful since may reveal occluded regions of the object that can be used to complete the model. By estimating the rigid and non-rigid components of the scene flow, as we propose, is possible to register close RGBD frames into a reference view to create a partial 3D model. A priori information of the object is required in order exploit the observations from other frames that are not visible for the reference view, allowing the completion of the partial 3D model. This process can be repeated for a set of reference views and partials model can be combined to have a variable 3D model of the object that varies according to the given view.

Discussion

The two questions studied in this work are inspired by a common goal, the confident estimation of scene flow for a variety of scenarios. The real-time ability of determining the 3D motion field of a scene can yield computer vision applications to take full advantage of the changing 3D world. This thesis rethinks the way 3D motion is estimated, to take the most possible advantage from current RGBD sensors. We strongly believe that the ideas developed in this work will take us a step closer to reaching this goal.

We have presented a new framework to compute dense scene flow from RGBD images by exploiting the semi-rigid properties of real world scenes. By modeling motion by a dense field of twist motions, we are able to directly encourage local and piecewise rigid solutions. We demonstrate that this motion representation allows very accurate estimations in a set of interesting experiments. Current computational tools provide alternatives to develop an implementation that will be suited for practical applications. Moreover, recent advances of optimization on manifolds will enable the exact formulation of the problem and its solution.

In order to exploit RGBD images we have developed a projective function supporting different motion representations. Using this warp is possible to directly exploit color and depth images provided by a RGBD sensor. In addition to the classical brightness-based constancy constraints, we introduce a depth constraint which allows the fully exploitation of depth data. By jointly minimizing the brightness and depth constraints we incrementally solve for the scene flow that best explains the provided images. The variational solution is directly formulated in the image domain, showing be able to accurately estimate rigid and nonrigid motions. Alternative representations of the scene can be required, particularly for 3D reconstruction applications, and discrete optimization can be used in addition to the variational formulation, to deal with larger displacement. Future advances of RGBD sensors will benefit the proposed framework, enabling access to a better quality data.

In this thesis we go beyond previous RGBD scene flow methods, by generalizing some of them and reformulating the way physical properties are used for 3D motion estimation. However, the problem scene flow estimation is not solved and there are several challenges waiting for being resolved, as some of those we have described above. Particularly, interesting challenges are posed by the requirement of measuring the confidence of components of the RGBD data and reliably detecting partial occlusions, which can take us to a closer point from a confident scene flow estimation.

  -Clés: mouvement, profondeur, flot de scène, rgbd, variationnelle, semiv X ,v Y ,v Z ) T 3D motion vector ξ =(ω, τ ) T twist motion xii ω =(ω X , ω Y , ω Z ) T 3D rotation vector τ =(τ X , τ Y , τ Z ) T

Figure 2 .

 2 1 shows a pair of color and depth images captured by a Kinect for Xbox.

Fig. 2 . 1 :

 21 Fig. 2.1: RGBD representation. (a) RGB image and (b) depth image. The depth image is encoded in gray levels, where the closer a point is at the brighter is drawn. Pixels having no valid depth measures are drawn in black. Note that field of views of RGB and depth cameras are different.

Fig. 2 . 3 :

 23 Fig. 2.3: Motion of a rigid surface point in 3D space and in the image domain

Fig. 2 . 4 :

 24 Fig. 2.4: From left to right: (a) 3D and (b) 2D motion fields represented in the image domain.

Fig. 2 . 5 :

 25 Fig. 2.5: Rigid body motion represented as a rotation followed by a rotation.

Fig. 2 . 6 :

 26 Fig. 2.6: Motion induced by twist ξ =( ω, τ ).T h e 3 D r o t a t i o n a l c o m p o n e n t i s e x p r e s s e d a s ω = θω ,w h e r es c a l a rθ = |ω| is the amount of rotation, and the unitary vector denoted by the over-arrow,ω = ω/|ω|,i st h ea x i so fr o t a t i o n . T h et r a n s l a t i o n a lc o m p o n e n t ,τ ,g i v e st h em o t i o n along the 3D axis by kωτ ,w h e r ekω =ω/|ω| 2 .T h e3 Da x i si sal i n ei nω direction going through the point ωτ/|ω| 2 .

Fig. 2 .

 2 Fig. 2.7: ROF-model denoising. From left to right: (a) original image, (b) noisy image and (c) denoised image. (Images courtesy of Mark Werlberger).

Fig. 2 . 8 :

 28 Fig. 2.8: Comparison of regularization strategies using the Army optical flow dataset by Bakeret al. (2007). From left to right: (a) ground truth of the motion, (b) motion estimate using l2-regularization[START_REF] Horn | Determining optical flow[END_REF], and (c) motion estimate using TV-regularization[START_REF] Wedel | Statistical and Geometrical Approaches to Visual Motion Analysis[END_REF]. Unlike l2-regularization that over-smooths the solution because its quadratic penalty, TV respects motion discontinuities thanks to its l1-based penalty. (Images courtesy of Mark Werlberger)

Fig. 3 . 1 :

 31 Fig. 3.1: Difference between apparent and real motion. From left to right: (a) a barber's pole rotating to the right, (b) optical flow and (c) 2D motion field. While the pole is rotating, the color patterns make perceive an upward motion, which does not correspond with the motion field.

Fig. 3 . 2 :

 32 Fig. 3.2: Local and piecewise rigidity assumptions. (a) Input image, (b) local rigidity and (c) piecewise rigidity. Every color patch indicates that inner scene points present (approximated) the same rigid motion. Local rigidity favor the best rigid motion explaining the motion inside every patch. Piece-wise rigidity encourage a set of rigid moving segments.

Fig. 5 . 1 :

 51 Fig. 5.1: RGBD changes using a fixed camera. Form left to right, top to bottom: (a) first image, (b) second image, and absolute difference in (c) intensity and (d) depth. Moving surfaces produce changes in intensity and depth, while static scene points remain unchanged, except by the noise.

Fig. 5 . 2 :

 52 Fig. 5.2: RGBD changes with a moving camera. From left to right, top to bottom: (a) first image, (b) second image, and absolute difference in (c) intensity and (d) depth.The movement of the camera, as wel l as that of the person, produce a 3D motion of every scene point. Particularly, the head turning produces more changes in intensity than in depth, while the camera motion is better perceived on depth data in some regions, such as neck and t-shirt.

Fig. 5 . 3 :

 53 Fig. 5.3: Brightness constancy violation. From left to right, top to bottom: (a) first image, (b) second image, (c) image difference, (d) first and (e) second images of the magnitude of the gradient, and (f ) their difference. While the camera is moving, there exists an automatic white balance that modifies the gray value image, as can be noted in the image difference (c), changes that are not related to the motion. On the other hand, the magnitude of the gradient is too little affected since it measures the image structure, or texture, regardless uniform changes in luminance.

Fig. 5 . 4 :

 54 Fig. 5.4: Truncated Charbonnier penalty. From left to right: (a) truncated version of the robust norm Ψ s 2 = p s 2 +0.001 2 ,f o rag i v e nS 2 tr ,a n d( b )i t sd e r i v a t i v e . T h i st r u n c a t e dp e n a l t yd r o p s all square residues larger than S 2 tr from the local estimation.

Fig. 5 . 5 :

 55 Fig. 5.5: Depth-based discontinuities preserving. From left to right: (a) depth image, Z(x),( b ) magnitude of the gradient, |rZ(x)| ,a n d( c )d e p t h -b a s e dw e i g h t i n gf u n c t i o n ,c(x).

  .21) with krvk := |rv X | + |rv Y | + |rv Z | and c(x) the weighting function (5.15).

2 S 1 2Fig. 5 . 6 :

 2156 Fig. 5.6: RGBD image pyramid of 2 levels. Intensity and depth pyramids are built with a downsampling factor of 2. Camera parameters f = {fx,cx} and c =( cx,xy) are scaled at each level l,w i t hl = {0, 1, 2},b yt h ef a c t o r2 l .T h efi e l do fr i g i dm o t i o nξ is directly propagate from one level to the next, starting from the coarsest level, given by images {S 2 1 ,S 2 2 },t ot h efi n e s t ,i m a g e sa t the original resolution {S 0 1 ,S 0 2 }.

Fig. 5 .

 5 Fig. 5.7: Rigid motion compensation. From left to right, top to bottom: (a) first image, (b) second image, (c) difference between input images, (d) second image warped to compensate rigid motion, (e) image blended using first image and the compensated, and (f ) difference between first image and the compensated. Once the rigid motion has been compensated the motion estimation becomes the computation of the nonrigid residual of the motion.

Fig. 6 . 1 :

 61 Fig. 6.1: Teddy stereo dataset. From left to right: (a) first and (b) second RGBD image pairs. Every depth image is computed from the disparity map of the corresponding view.

Fig. 6 . 2 :

 62 Fig. 6.2: Optical flow for the Teddy dataset. From left to right: (a) Blend of first and second brightness images, (b) ground truth of the optical flow, and (c) optical flow color coding. Observe that every pixel performs a purely horizontal motion, whose magnitude is a function of its depth.

  and the motion detail preserving algorithm by Xu et al. (2012) local neighborhood, N ⇥ N centering window M number of alternations in each level of the pyramid M TV iterations of the TV solver M GN iterations of the Gauss-Newton algorithm ⇣ local-rigidity parameter (3DmotionFlow) N ζ size of the local neighborhood for the rigid fitting (3DmotionFlow)

  -based methods do not assume RGBD images and simultaneously estimate the optical flow and disparity maps by considering RGB images 2, 4, 6 and 8 of each dataset. Images 2 and 4 form the first stereo pair, and images 6 and 8 the second one. We consider the variational method from stereo by Huguet and Devernay (2007)( StereoFlow1 ), which was the first to use these datasets for scene flow, and the multi-view variational algorithm by Basha et al. (2010)( StereoFlow2 ). Scene flow methods RGBD-based methods, as those proposed in this thesis, use simultaneously RGB and depth images for the estimation. We consider two scene flow methods: the particle filter based algorithm by Hadfield and Bowden (2014)( SceneFlow1 ) and the sphere flow algorithm by Hornacek et al. (2014)( SceneFlow2 ).

Fig. 6 . 3 :

 63 Fig. 6.3: Averages of Teddy and Cones datasets for the motion errors RMSE, AAE and ANEZ , as functions of the size of the local neighborhood.

Fig. 6 . 4 :

 64 Fig. 6.4: Optical flow estimates for the Teddy stereo dataset. From left to right, from top to bottom: (a) Ground-truth, and (b) twist-motion, (c) 3D-motion and (d) 2D-motion representations.

Fig. 6 . 5 :

 65 Fig. 6.5: Optical flow estimates for the Cones stereo dataset. From left to right, from top to bottom: (a) Ground-truth, and (b) twist-motion, (c) 3D-motion and (d) 2D-motion representations.

Fig. 6 . 6 :

 66 Fig. 6.6: Inputs for scene flow estimation. From left to right: (a) brightness and (b) 3D surface of the first and second frames. The 3D surface is illustrated using directly the depth image, where ad e c r e a s i n gf u n c t i o no ft h em a g n i t u d eo ft h eg r a d i e n tm o d u l a t e st h ev i s u a l i z a t i o n ,g i v i n ga ni d e a of the 3D structure. Only the 3D surface within 50 cm and 200 cm from the sensor is shown.

Fig. 6 . 7 :

 67 Fig. 6.7: Scene flow results. From left to right, and from top to bottom: (a) optical flow, and (b) X-, (c) Y -a n d( d )Z-components of the scene flow, respectively.

Fig. 6 . 8 :

 68 Fig. 6.8: Brightness and depth residuals. From left to right: (a) warped image, and (b) brightness and (c) depth residuals, respectively. The top row shows results assuming a null scene flow.The bottom row presents the warped image and differences with the estimated scene flow. Brightness differences are scaled by 8 for visualization, and depth differences go from black to white, for 0 cm to 1 cm, and which for larger differences. Cyan color stands for an invalid depth difference.

Fig. 6 .

 6 Fig. 6.9: Large Z-motion experiment proposed by Hornacek et al. (2014). From left to right: (a) color image, (b) 3D surface and (c) RGBD image. First and second images are presented in the top and bottom rows, respectively. The experiment estimates the motion in both directions.

Fig. 6 .

 6 Fig.6.10: Optical flow results of the large Z-motion experiment presented by[START_REF] Hornacek | Sphereflow: 6 dof scene flow from RGB-D pairs[END_REF]. The top images present the estimation from frame 1 to 2, while the bottom images for the contrary sense. Observe, that only (a) and (b) are able to estimate the motion in both cases.

Fig. 6 .

 6 Fig. 6.11: Optical flow results of the large Z-motion experiment. All methods are able to overcome the large change in depth and give a roughly estimation of the performed motion.

Fig. 6 .

 6 Fig.6.12: Warped images for the large Z-motion experiment. The use of depth data allows a better modeled of the observed motion. Particularly, 6TwistFlow performs the most accurate estimate of the Z-motion. None of the method is able to model the large head rotation on the left side.

Fig. 6 .

 6 Fig.6.13: Changes in depth for the large Z-motion experiment. 6TwistFlow and 3DmotionALL are able to accurately estimate the changes in depth. Conversely, the optical flow method suffers from the 2D regularization and the lack of depth data, failing to estimate the motion boundaries.

Fig. 6 .

 6 Fig. 6.14: Inputs for the hand rotation experiment. From left to right: (a) first and (b) second pairs, of brightness and 3D structure images, respectively.

Fig. 6 .

 6 Fig. 6.15: Results of the hand rotation experiment. From top to bottom: (a) optical flow, and (b) X-component and (c) Z-component of the scene flow.

Fig. 6 .

 6 Fig. 6.16: Inputs for the two arms rotation experiment. From left to right: (a) first and (b) second image pairs, of brightness and 3D structure.

Fig. 6 .

 6 Fig.6.17: Results of the two arms rotation experiment. From top to bottom: optical flow and Z-component of the scene flow. LDOF provides a smoother optical flow estimation, but it does not correspond with the most accurate scene flow estimate. On the other hand, 6TwistFlow is more sensitive to the small rotation and deformation, and is able to model the motion of thumbs, palm hands and elbows. 3DmotionALL also provides an accurate estimation.

Fig. 6 .

 6 Fig. 6.18: Inputs for the poster stretching experiment. From left to right: (a) first and (b) second image pairs, of brightness and 3D structure.

Fig. 6 .

 6 Fig. 6.19: Results of the poster stretching experiment. From top to bottom: optical flow and Z-component of the scene flow. 3DmotionALL achieves the best performance in this experiments.

Fig. 6 .

 6 Fig. 6.21: Inputs for the teddy bear experiment. From left to right: (a) first and (b) second image pairs, of brightness and 3D structure.

Fig. 6 .

 6 Fig. 6.22: Optical flow results of the Teddy bear experiment. From Left to right, top to bottom: (a) rigid estimation, (b) 6TwistFlow, (c) 3DmotionALL and (d) LDOF.

Fig. 6 .

 6 Fig. 6.23: Inputs for the rigid motion experiment. From left to right: (a) first and (b) second image pairs, of brightness and 3D structure.

Fig. 6 .

 6 Fig. 6.24: Results of the rigid motion experiment. From left to right: (a) optical flow and (b) Z-component of the scene flow, respectively.

1 Fig. 6 .

 16 Fig. 6.25: Brightness and depth residual while varying λ.B r i g h t n e s s a n d d e p t h r e s i d u a l s a r e presented in top and bottom images, respectively. For λ =0,onlythebrightnessc onstr aintsar euse d and the depth residual presents higher values, as is expected. On the other hand, if λ = 1,o n l yt h e depth constraints are used, and the depth residuals achieve their minimum while brightness residuals becomes large. Setting λ =0.1 balances the contribution of brightness and depth constraints, jointly minimizing both brightness and depth residuals.

Fig. 6 .

 6 Fig. 6.26: Teddy bear 3D reconstruction, images from the RGBD dataset (Sturm et al., 2012). From left to right: (a) reference image, and 3D reconstruction using (b) the reference depth image and (c) by warping 12 close views in addition to the reference depth image.

Fig. 6 .

 6 Fig. 6.27: Rotating hand 3D reconstruction. From left to right: (a) reference image, and 3D reconstruction using (b) the reference depth image and (c) by warping 10 close views in addition to the reference depth image.

Fig. 6 .

 6 Fig. 6.28: Reference and limit frames for the 3D reconstruction experiment. Brightness and 3D structure images are presented in top and bottom images, respectively.

1 Fig. 6 .

 16 Fig.6.29: Varying λ in the 3D reconstruction experiment. The only use of brightness constraints is highly affected by specular highlights, occlusions and the large rotation. On the other hand, the formulation of constraints only on the depth images, improves the quality of the 3D reconstruction. The jointly use of brightness and depth constraints, outperforms the only-brightness estimation but is less accurate than the results by the only-depth approach.

Fig. 6 .

 6 Fig.6.30: Inputs for the rigid plus nonrigid scene flow estimation using the Xtion sensor. Brightness and 3D structure images are presented in top and bottom images, respectively. While the camera is rotating to the right, the head rotates to the left and hands moves in opposite directions, one becoming partially occluded.

Fig. 6 .

 6 Fig.6.32: Warped images for the rigid plus nonrigid scene flow estimation using the Xtion sensor. Observe that the rigid motion correctly compensates for the body motion. The nonrigid estimation is able to roughly register the rotating head, as is noted by observing the right location of glasses, the beard and the ear. The motion estimation of the left hand is not well modeled.

Fig. 6 .

 6 Fig.6.33: Inputs for the rigid plus nonrigid scene flow estimation using the Kinect sensor. Brightness and 3D structure images are presented in top and bottom images, respectively. While the camera is subtly moving, the head rotates and lifts.

Fig. 6 .

 6 Fig. 6.34: Results of the rigid plus nonrigid scene flow estimation using the Kinect sensor.

Fig. 6 .

 6 Fig. 6.35: Warped images for the rigid plus nonrigid scene flow estimation using the Kinect sensor.

Fig. 6 .

 6 Fig.6.36: Residuals of the rigid plus nonrigid scene flow estimation using the Kinect sensor. Brightness and depth residuals are presented in top and bottom images, respectively. Observe that high and undetermined residuals can be used to reliably estimate the partial occlusions. On the other hand, occlusions are also visible in the brightness residual, but can be confusable with the specular highlights which also present high residual.

  .2)

		S t		S t+1
	X t	v(x t )		X t+1
	Y
	Z
	X	x t	(u, v)	x t+1
	Camera Frame
	y
				Image Plane
		x
		Optical Center

  becomes the problem

	ψ(s 2 )	Quadratic	ψ '(s 2 )	Quadratic
		Lorentzian		Lorentzian
		Charbonnier		Charbonnier
		s 2		s 2
	Fig. 2.9: Comparison of quadratic, Lorentzian (a =1)a n dC h a r b o n n i e r( ε =0.001)p e n a l t i e s .

Table 4 .

 4

1: Summary of motion representations.

Table 6 . 1 :

 61 Parameters of baseline scene flow methods.

Table 6 . 2 :

 62 Comparison with other methods using Middlebury datasets. First and second best performances are highlighted in red and orange, respectively.

			Te d d y			Cones			Ve nu s
		RMSE	AAE	RMSE Z	RMSE	AAE	RMSE Z	RMSE	AAE	RMSE Z
	OpticalFlow1	2.83	0.39	1.75	3.20	0.39	0.47	0.72	1.28	0.14
	OpticalFlow2	1.66	0.21	1.15	1.70	0.28	0.50	0.3	1.43	0.22
	StereoFlow1	1.10	0.69	n/a	1.25	0.51	n/a	0.31	0.98	n/a
	StereoFlow2	0.58	0.39	n/a	0.57	1.01	n/a	0.16	1.58	n/a
	SceneFlow1	1.24	1.01	0.06	0.83	0.83	0.03	0.36	1.03	0.02
	SceneFlow2	0.54	0.52	0.02	0.35	0.15	0.01	0.26	0.53	0.02
	6TwistFlow	0.35	0.35	0.02	0.40	0.34	0.02	0.15	0.84	0.01
	3DmotionFlow	0.71	0.67	0.04	0.51	0.31	0.04	0.2	1.02	0.00

Table 6 . 4 :

 64 Local rigidity in the smoothness term. Top results are shown in red. The best results are achieved when the local-rigid prior in the smoothness term is applied in addition to the local rigidity assumption in the data term.

			Te d d y			Cones
		RMSE	AAE	ANEv	P5%	RMSE	AAE	ANEv	P5%
	3DmotionFlow	0.71	0.67	2.01	89.7	0.51	0.34	1.37	98.1
	3DMotionE S	1.41	1.27	3.80	78.2	1.50	1.49	4.26	77.7
	3DMotionALL	0.68	0.64	1.98	89.9	0.48	0.32	1.30	98.2

Table 6 . 5 :

 65 Variation of the piecewise rigidity assumption. Top results in red and blue for 3DmotionFlow and 6TwistFlow, respectively. Baseline methods are highlighted in green.

			Te d d y			Cones
		RMSE	AAE	ANEv	P5%	RMSE	AAE	ANEv	P5%
	GroundTruth	0.13	0.00	0.00	100.	0.24	0.00	0.00	100.
	6TwistMotion	0.35	0.35	1.63	100.	0.40	0.34	1.22	100.
	3DmotionFlow	0.71	0.67	2.01	89.7	0.51	0.34	1.37	98.1
	2DmotionFlow	3.44	0.89	10.4	56.6	3.04	0.74	10.7	68.2

Table 6 . 6 :

 66 Comparison of motion representations. Groundtruth results are highlighted in red. The twist-motion representation achieves the best performance.

Table 6 . 7 :

 67 Variation of the balance between brightness and depth. Top results in orange, red and blue for 2DmotionFlow, 3DmotionFlow and 6TwistFlow, respectively. Baseline methods are highlighted in green. The simultaneous constraint in brightness and depth achieves the best results.

				Te d d y			Cones
		λ	RMSE	AAE	ANEv	P5%	RMSE	AAE	ANEv	P5%
		0	3.41	0.94	14.2	57.0	3.12	0.79	26.1	66.6
	2DmotionFlow	0.1 1	3.44 3.63	0.89 1.42	10.4 11.6	56.6 50.9	3.04 4.06	0.74 0.78	10.7 12.3	68.2 64.8
		1	9.25	9.78	28.0	28.4	7.07	2.92	46.7	25.0
		0	1.12	0.91	13.2	28.2	1.19	0.94	15.7	22.1
	3DmotionFlow	0.1 1	0.98 0.71	0.76 0.67	2.38 2.01	89.1 89.7	0.72 0.51	0.60 0.34	1.85 1.37	93.3 98.1
		1	3.26	4.33	8.47	67.7	1.18	0.82	3.13	86.3
		0	0.38	0.36	1.66	100.	0.47	0.42	4.09	70.2
	6TwistFlow	0.1 1	0.35 0.54	0.35 0.57	1.63 2.21	100. 93.17	0.40 0.45	0.34 0.25	1.22 1.18	100. 99.2
		1	0.88	0.59	3.5	70.3	1.28	1.19	3.85	67.5

Acknowledgments

I would like to thank my advisors Frederic Dévérnay and James Crowley for the opportunity of doing this thesis under their supervision. I also would like to thank Thomas Brox for the interesting discussions and for providing a nice research environment at University of Freiburg.

I would like to especially thank Sergi Pujades and Laurent Boiron for their unconditional support and encouragement while this rough road trip. I must also specially thank INRIA Grenoble for providing a wonderful research environment and for supporting me in every research activity while this period.

This project would have been impossible without the support of the Pontificia Universidad Javeriana and Colciencias. Also I am grateful for the support of Campus France, Region Rhones-Alpes and GdR ISIS.

Algorithm 1 Scene flow estimation using the twist-motion representation.

Super-scripted l denotes the pyramid level.

Input: two RGBD images S 1 = {I 1 ,Z 1 } and S 2 = {I 2 ,Z 2 } Output: scene flow v from S 1 to S 2

Construct RGBD pyramids S l 1 and S l 2 , l =0,...,L; Initialize ⇠ =(!, ⌧ )=0, χ =($, ⇡)=0, and l = L; while l ≥ 0 do for w =1toM do for n =1toM GN do Warp images I l 2 and Z l 2 using ⇠; Eq. ( 4. Algorithm 2 Scene flow estimation using the 3D-motion representation.

Super-scripted l denotes the pyramid level.

Construct RGBD pyramids S l 1 and S l 2 , l =0,...,L; Initialize v = 0, w = 0, u = 0, and l = L; while l ≥ 0 do for w =1toM do for n =1toM GN do Warp images I l 2 and Z l 2 using ⇠; Eq. ( 4.13) Calculate residues ⇢ l I , ⇢ l g and ⇢ l Z ; Update v; Eq. ( 5.56) end for for n =1toM TV do Iterate once for solving u; Eq. ( 5.42) end for Encourage local rigidity on w; Eq. ( 5.63) end for Median filter on v; end while Algorithm 3 Rigid plus nonrigid scene flow estimation.

Construct RGBD pyramids S l 1 and S l 2 , l =0,...,L; Initialize ⇠ R = 0, ⇠ = 0, and l = L; while l ≥ 0 do for k =1toK do for w =1toM do for n =1toM GN do Warp images I l 2 and Z l 2 using χ = log(e ξ + e ξR -I 4⇥4 ); Eq. ( 4.3) Calculate residues ⇢ l I , ⇢ l g and ⇢ l Z ; Update ⇠; Eq. ( 5.38) end for for n =1toM TV do Iterate once for solving $; Eq. ( 5.42) Iterate once for solving ⇡; Eq. ( 5.46) end for end for for n =1toM GN do Warp images I l 2 and Z l 2 using χ = log(e ξ + e ξR -I 4⇥4 ); Eq. ( 4.3) Calculate residues ⇢ l I , ⇢ l g and ⇢ l Z ; Update ⇠ R ; Eq. ( 5.67) end for end for end while Compute the scene flow v from ⇠; Eq. ( 4.2)

List of Tables

List of Algorithms