Second Order Variational Optic Flow Estimation

Second Order Variational Optic Flow Estimation L. Alvarez, C.A. Castaño, M. Garcı́a, K. Krissian, L. Mazorra, A. Salgado, and J. Sánchez Departamento de Informática y Sistemas, Universidad de Las Palmas de Gran Canaria, Campus de Tafira s/n, 35017 Las Palmas de Gran Canaria, Spain Tel.: +34 928 458708, Fax: +34 928 458711 {lalvarez,ccastano,mgarcia,krissian,lmazorra,asalgado, jsanchez}@dis.ulpgc.es Abstract. In this paper we present a variational approach to accurately estimate the motion vector field in a image sequence introducing a second order Taylor expansion of the flow in the energy function to be minimized. This feature allows us to simultaneously obtain, in addition, an estimation of the partial derivatives of the motion vector field. The performance of our approach is illustrated with the estimation of the displacement vector field on the well known Yosemite sequence and compared to other techniques from the state of the art. 1 Introduction Optic flow estimation is a problem that has focused the attention of many researchers in the domain of image processing, probably due to the huge amount of applications in which the estimated displacement between consecutive frames is an important source of information. In this sense, many different approaches have been presented in the literature, starting from classical methods like the ones proposed by Horn and Schunck [1] or by Lucas and Kanade [2], trying either to overcome some limitations of the existing methods or to exploit some a priori information in order to improve the accuracy on the estimation. For instance, different regularization terms have been proposed in variational approaches like in [3,4], pyramidal decompositions have also been proposed in order to detect large displacements [5,6] or spatio-temporal regularization constraints have also been taken into account [7,8]. In this paper, we propose a variational formulation to accurately estimate the motion vector field in a image sequence introducing a second order Taylor expansion of the flow in the energy function to be minimized. This feature allows us to simultaneously obtain, in addition, an estimation of the partial derivatives of the motion vector field. This idea is quite interesting in the field of fluid dynamics, since partial derivatives is an important source of information in this field as they appear in the Navier-Stokes equations and other intersting parameter such as divergence, vorticity, strain rate tensor and dissipation rate. Moreover, since a regularity constraint is also imposed on the second order flow parameters, the estimated flow will better preserve the continuity behavior assumed in fluid dynamics. R. Moreno-Dı́az et al. (Eds.): EUROCAST 2007, LNCS 4739, pp. 646–653, 2007. c Springer-Verlag Berlin Heidelberg 2007  Second Order Variational Optic Flow Estimation 647 The paper is organized as follows. In section 2, we describe the details of our variational approach and the way we have adapted the energy function to be minimized in order to directly estimate significative fluid parameters such as vorticity and strain rate tensor components. In section 3, we present numerical experiments on synthetic and real data and we analyze the behavior of the proposed method in comparison with other standard techniques described in the literature. Finally, in section 4, we present the main conclusions of the paper. 2 A Variational Approach for Second Order Motion Estimation To estimate the optic flow of a given sequence we propose a variational approach based on the minimization of an energy function E(ũ) defined for each point x0 with the special fact that it depends not only on the displacement vector components u = (u, v)T , but also on their partial derivatives ũ = (u, v, ux , uy , vx , vy )T , as shown in the following equation: E(ũ) = E(u, v, ux , uy , vx , vy ) = 2     ux uy T (x − x0 )) dx (1) Ω(x0 ) Kσ (x − x0 ) I1 (x) − I2 (x + (u, v) + vx vy where I1 (x) is the first image of the sequence and I2 (x) is the following frame, where the time increment from frame to frame is assumed to be normalized (∆t = 1), Kσ (x − x0 ) is a Gaussian kernel with standard deviation σ which weighs the pixels in the domain Ω(x0 ) centered at point x0 . Hence, the goal is to find the components of vector ũ, that is, the displacement and its partial derivatives such that minimize the error between I1 (x) and I2 (x) displaced by the unknown flow. At a first glance, the dependence on the partial derivatives (ux , uy , vx , vy ) might seem a redundant operation since they can be computed from the obtained motion vector field. However, it can be shown that numerical differentiation of the estimated motion vector field yields to inaccurate results mainly due to undesired numerical error amplification [9]. 2.1 Energy Minimization In order to be able to obtain the six components vector of unknowns ũ = (u, v, ux , uy , vx , vy )T that minimize Eq. 1 we formulate the solution at step n + 1 as a function of the solution at step n and the six components vector of residuals h̃ = (hu , hv , hux , huy , hvx , hvy )T computed at each step, as shows the following expression: ũn+1 = ũn + h̃ . (2) Then, introducing Eq. 2 in Eq. 1 we are able to obtain an iterative operator towards the local minimum of the energy function similar to a gradient descent 648 L. Alvarez et al. algorithm [10]. But first, the following approximations are still necessary to simplify the term that depends on I2 (x) in order to obtain our iterative operator:  n+1 n+1  ux uy n+1 n+1 T (x − x0 )) = I2 (x + (u ,v ) + vxn+1 vyn+1    n n ux uy hux huy I2 (x + (un , v n )T + (x − x0 ) + (hu , hv )T + (x − x0 )) . hvx hvy vxn vyn To simplify notation, let us define I2n as:   n n ux uy x − x0 ) I2n (x) = I2 (x + (un , v n )T + y − y0 vxn vyn (3) where the vector x − x0 has been decomposed into their components (x− x0 , y − y0 )T . The partial derivatives are denoted as:  n n  ∂I2 ux uy x − x0 n (x) = I2,x (x + (un , v n )T + ). (4) vxn vyn y − y0 ∂x n I2,y (x) = ∂I2 (x + (un , v n )T + ∂y  unx uny vxn vyn   x − x0 ). y − y0 (5) Then, using the new notation, the solution at step n + 1 can be approximated by its Taylor expansion over ũn to provide: I2n+1 (x) ≃ I2n (x)     T  n I2,x (x) hux huy x − x0 T . + (hu , hv ) + n I2,y (x) hvx hvy y − y0 (6) Finally, by using vectorial notation I2 (x)n+1 can be expressed as: T I2n+1 (x) ≃ I2n (x) + În 2 (x) · h̃ (7) where n n n n În 2 (x) = (I2,x (x), I2,y (x), (x − x0 )I2,x (x), (y − y0 )I2,x (x), n n (x − x0 )I2,y (x), (y − y0 )I2,y (x))T . (8) Hence, introducing Eq. 7 in Eq. 1 we obtain the approximation for the energy function shown in Eq. 9, where, in addition, we introduce a regularization term in the energy function weighted by the parameter α. The role of this regularization term is to provide at every point a smooth vector field, by means of an additional constraint on the norm of the vector h̃, which is forced to be as small as possible. In this sense, the parameter α determines the importance of this additional constraint on the vector field.    2 T dx+α  h̃ 2 dx. Kσ (x−x0 ) I1 (x) − I2n (x) − În (x) · h̃ Ẽ(h̃) = 2 Ω(x0 ) Ω(x0 ) (9) Second Order Variational Optic Flow Estimation 649 This formulation allows us to easily obtain an analytical expression to compute the local minimum of the energy as a function of h̃:    T ∇Ẽ(h̃) = −2 Ω(x0 ) Kσ (x − x0 ) I1 (x) − I2n (x) − În În (x) · h̃ 2 (x)dx + 2  (10) 2α Ω(x0 ) h̃dx = 0 which is equivalent to:  Ω(x0 )  Ω(x0 ) Kσ (x − x0 ) (I1 (x) − I2n (x)) În 2 (x)dx =   T n Kσ (x − x0 )(În 2 (x) · Î2 (x) ) + αI h̃dx . (11) This system of equations can be expressed with the standard matrix notation Ax = b, where the vector of unknowns in this case is the vector h̃, while the system matrix and the independent term can be computed as:    n T A= Kσ (x − x0 )(În (x) · Î (x) ) + αI dx . (12) 2 2 Ω(x0 )  b= (13) Kσ (x − x0 ) (I1 (x) − I2n (x)) În 2 (x)dx . Ω(x0 ) The solution of the system of equations is given by h̃ = A−1 b, so we only have to invert the 6 × 6 matrix A or use any other algorithm to solve the system of equations. We observe that if α > 0, matrix A is positive definite and therefore the system of equations is well-posed (it has a unique solution). It means that the parameter α avoids instabilities in the solution of the linear system of equations. 2.2 Numerical Implementation Next, we describe the main steps we followed to derive an efficient algorithm of the method proposed in section 2, including the numerical considerations involved on the computation of integrals and spatial derivatives in the proposed method. As it can be seen in Algorithm 1, our variational approach starts from an initial estimation of the motion vector field u0 which can be obtained using any other optic flow estimation method. Then, we start the iterative procedure towards the local minimum of the energy function given in Eq. 1. At each iteration we check the convergence of our algorithm in order to discard the solutions that do not provide the optimal response. This algorithm is then applied for every point in the image (or a grid at a given scale) to obtain the desired first and second order flow parameters estimation, using information within a neighborhood of that point determined by the σ, parameter of the Gaussian kernel. The initialization of the energy at the first step requires the computation of the image derivatives. In our implementation, we have used central finite differences because it is a good compromise between an easy implementation, low time of computation and low error propagation [9]. 650 L. Alvarez et al. Algorithm 1. Implementation of the second order flow estimation method Initialization of the parameters: σ (window size), α, ũ0 = 0, E(ũ0 ), N iter, F ailures = 0, M axF ailures Initial Estimation of the motion field u0 = (u, v) with any estimation flow technique for n = 0 to N iter do n Update I2n (x) and În 2 using ũ Computation of the energy E(ũn ), A and b Solve system of equations h̃ = A−1 b if E(ũn ) ≤ E(ũn−1 ) then ũoptimized = ũn + h̃ {The method converges towards the solution} Update ũn+1 F ailures = 0 else if E(ũn ) ≥ E(ũn−1 ) and F ailures ≤ M axF ailures then Update ũn+1 Reject this solution as ũoptimized {The method diverges from the solution} F ailures + + else Exit end if end for 3 Results In order to show the performance of the proposed approach, we present here the numerical experiments we have performed to evaluate the accuracy of the method, working on both synthetic and real data. It is interesting to remark that the initial estimation of the flow can be obtained by any desired method, which must be, at least, a rough approximation of the underlying motion vector field so that the algorithm is able to converge. In our work, we use the method described in [11] to achieve such initialization. 3.1 Experiment 1: Yosemite Sequence In our first experiment we use the well-known Yosemite sequence to evaluate the performance of our approach. From a qualitative point of view, it is interesting to notice that the vector field obtained with our approach is usually smoother than the initial estimation. This is more clearly seen on Fig. 1, where on the left side we represent the angular error between the groundtruth and the initial estimation, while on the right we show the angular error for our response. As it can be seen, this error feature is lower and more regular with our method. From a quantitative point of view, table 1 shows some statistics on the euclidean error (Eq. 14) and the angular error (Eq. 15) for that sequence according to the ideas in [12]. The error was computed removing the clouds in order to be able to provide a meaningful error measure, following the recommendation in [13]. As it Second Order Variational Optic Flow Estimation 651 Fig. 1. On the left, angular error between the groundtruth and the initial estimation of the vector field. On the right, the vector field obtained with our method. Table 1. Quantitative error measures obtained for the initial estimation and for our method Initialization Our Method Angular Error Std. Dev. Euclidean Error Std. Dev. 1.6696 1.8765 1.6598 0.01 1.5432 1.8467 1.6598 0.01 can be seen in the referred table, the error rate is lower once we use our method to process the initial vector field. N 1  |ui − uref i | . N i=1 (14) 180  arccos(ui · uref i ) . N π i=1 (15) RM SV D = N ψE = 3.2 Experiment 2: Satellite Sequence Finally, we present the results we obtained with our approach using real satellite image sequence provided by the Laboratoire de Meteorologie Dynamique (LMD) from Paris (France). This data tracks the hurricane Vince that became the first known tropical cyclone to reach the Iberian Peninsula between October 8th and 11th, 2005. The main interest in processing this kind of data is to estimate the cloud motion between consecutive frames. In Fig. 2 we compare the vector field obtained with our approach (represented in yellow, light arrows) with that used to initialize our method (represented in red, dark arrows). As it can be seen, the vector field obtained is more regular than the initial estimated flow. 652 L. Alvarez et al. Fig. 2. Details of the motion vector field in the satellite sequence. In red (dark arrows), we represent the method used to initialize our algorithm and in yellow (light arrows) we display the vector field obtained with our method. 4 Conclusions In this paper, we have presented a new variational method to estimate the motion vector field in a image sequence. In this sense, our variational method improves the quality of the motion vector field used as initial estimation introducing a regularity constraint on the first and second order moments of the motion vector field to be estimated. Our quantitative results on synthetic data show the encouraging results we obtain using our method as a post-processing step for the regularization process. The performance of our method is also tested with real satellite image sequences, where the vector fields we obtain also present a regular behavior. Acknowledgements This work has been funded by the European Commission under the Specific Targeted Research Project FLUID (contract no. FP6-513663). References 1. Horn, B., Schunck, B.: Determining optical flow. MIT Aritificial Intelligence Laboratory (1980) 2. Lucas, B., Kanade, T.: An iterative image registration technique with an application to stereo vision. In: Proc. Seventh International Joint Conference on Artificial Intelligence, Vancouver (1981) 674 – 679 3. Nagel, H.H., Enkelmann, W.: An investigation of smoothness constraints for the estimation of displacement vector fields from image sequences. IEEE Trans. Pat. Anal. and Mach. Intell. 8 (1986) 565–593 4. Aubert, G., Deriche, D., Kornprobst, P.: Computing optical flow via variational techniques. SIAM Journal on Applied Mathematics 60(1) (1999) 152–182 5. Mémin, E., Pérez, P.: Dense estimation and object-based segmentation of the optical flow with robust techniques. IEEE Transactions on Image Processing 7(5) (1998) Second Order Variational Optic Flow Estimation 653 6. Alvarez, L., Weickert, J., Sánchez, J.: Reliable estimation of dense optical flow fields with large displacements. International Journal of Computer Vision 39(1) (2000) 41–56 7. Weickert, J., Schnorr, C.: A theoretical framework for convex regularizers in pdebased computation of image motion. Draft, Departament of Mathematics and Computer Science (2000) 8. Weickert, J., Papenberg, N., Bruhn, A., Brox, T.: High accuracy optical flow estimation based on a theory for warping. Volume 4. (2004) 25–36 9. Raffel, M., Willert, C., Kompenhans, J.: Particle Image Velocimetry. A Practical Guide. Springer Verlag (1998) 10. Fletcher, R., Reeves, C.: Function minimization by conjugate gradients. Computer Journal 7 (1964) 149–154 11. Alvarez, L., Castaño, C., Garcı́a, M., Krissian, K., Mazorra, L., Salgado, A., Sánchez, J.: Symmetric optical flow. In: Proceeding of EuroCast, Las Palmas, Spain (2007) 12. Barron, J., Fleet, D., Beauchemin, S.: Performance of optical flow techniques. IJCV 12(1) (1994) 43–77 13. Black, M.: http://www.cs.brown.edu/people/black/sequences/yosfaq.html (2006) View publication stats