Reliable estimation of dense optical flow fields with large displacements
Abstract
56
Full text
Reliable Estimation of Dense Optical Flow Fields with Large Displacements Luis Alvarez 1 , Joachim Weickert 2 , and Javier Sanchez 1 1 Departamento de Informatica y Sistemas, Universidad de Las Palmas, Campus de Tara, SP-35017 Las Palmas, Spain. E-mail: f lalvarez,jsanchez g @dis.ulpgc.es WWW: http://serdis.dis.ulpgc.es/ f lalvarez,jsanchez g 2 Computer Vision, Graphics, and Pattern Recognition Group Department of Mathematics and Computer Science University of Mannheim D-68131 Mannheim, Germany E-mail: [email protected]-mannheim.de WWW: http://www.ti.uni-mannheim.de/ bmg/weickert Abstract In this pap er we show that a classic optical owtechnique by Nagel and Enkelmann (1986) can b e regarded as an early anisotropic diusion metho d with a diusion tensor. Weintro duce three improvements into the mo del form ulation that (i) avoid inconsistencies caused bycentering the brightness term and the smo othness term in dierent images, (ii) use a linear scale-space fo cusing strategy from coarse to ne scales for avoiding convergence to physically irrelevant lo cal minima, and (iii) create an energy functional that is invariant under linear brightness changes. Applying a gradient descent metho d to the resulting energy functional leads to a system of diusion{reaction equations. We prove that this system has a unique solution under realistic assumptions on the initial data, and we present an ecient linear implicit numerical scheme in detail. Our metho d creates ow elds with 100 % density over the entire image domain, it is robust under a large range of parameter variations, and it can recover displacement elds that are far b eyond the typical one-pixel limits which are characteristic for many dierential metho ds for determining optical o w. We show that it p erforms b etter than the classic optical ow metho ds with 100 % density that are evaluated by Barron et al. (1994). Our software is available from the Internet. Keywords: image sequences, optical ow, dierential metho ds, anisotropic diusion, linear scale-space, regularization, nite dierence metho ds, p erformance evaluation 1
1 Intro duction Optical ow computation consists of nding the apparent motion of ob jects in a sequence of images. Recovering this displacement eld is a key problem in computer vision and much research has b een devoted to this eld during the last two decades. For a survey of these activities we refer to Mitiche and Bouthemy 37], and p erformance evaluations of some of the most p opular algorithms include pap ers of Barron et al. 7], Jahne and Haussecker 32], and Galvin et al. 22]. One imp ortant class of optical ow metho ds consists of so-called dierential metho ds. Often they are considered as useful only in the case of small displacement elds. The goal of the present pap er is to show that acombination of linear and nonlinear scale-space ideas may leadtoawell-p osed dierential metho d that allows to recover the optical ow between two images with high accuracy, even in the case of large displacement elds. We consider two images I 1 ( x y )and I 2 ( x y ) (dened on R 2 to simplify the discussion) which represent two consecutive views in a sequence of images. Under the assumption that corresp onding pixels have equal grey values, the determination of the optical ow from I 1 to I 2 comes down to nding a function h ( x y )=( u ( x y ) v ( x y )) suchthat I 1 ( x y ) = I 2 ( x + u ( x y ) y + v ( x y )), 8 ( x y ) 2 R 2 : (1) To compute h ( x y ) the preceding equality is usually linearized yielding the so-called linearized optical ow constraint I 1 ( x ) ; I 2 ( x ) r I 2 ( x ) h ( x ) 8 x (2) where x := ( x y ). The linearized optical ow constraint is based on the assumption that the ob ject displacements h ( x ) are small or that the image is slowly varying in space. In other cases, this linearization is no longer valid. Frequently, instead of equation (1), the alternative equality I 1 ( x ; u ( x y ) y ; v ( x y )) = I 2 ( x y ), 8 ( x y ) 2 R 2 (3) is used. In this case the displacement h ( x y )is centred in the image I 2 ( x y ) : The determination of optical owis a classic ill-p osed problem in computer vision 10], and it requires to be supplemented with additional regularizing assumptions. The regularization by Horn and Schunck 28] reects the assumption that the optical ow eld varies smo othly in space. However, since many natural image sequences are b etter describ ed in terms of piecewise smo oth ow elds separated by discontinuities, much research has b een done to mo dify the Horn and Schunck approach in order to p ermit such discontinuous ow elds see 6, 11,12,14,17,18,24,26,34,36,42,43,46,49,55] and the references therein. An imp ortant improvement in this direction has b een achieved by Nagel and Enkelmann 42] in 1986 (see also 39]). They consider the following minimization problem: E NE ( h )= Z R 2 ( I 1 ( x ; u ( x y ) y ; v ( x y )) ; I 2 ( x y )) 2 dx (4) + C Z R 2 tr ace ; r h T D ( r I 1 ) ; r h dx 2
where C is a p ositive constant and D ( r I 1 )is a regularized pro jection matrix in the direction p erp endicular of r I 1 : D ( r I 1 )= 1 jr I 1 j 2 +2 2 8 > < > : 0 @ @I 1 @y ; @I 1 @x 1 A 0 @ @I 1 @y ; @I 1 @x 1 A T + 2 Id 9 > = > : (5) In this formulation, Id denotes the identity matrix. The advantage of this metho d is that it inhibits blurring of the ow across b oundaries of I 1 at lo cations where jr I 1 j . In spite of its merits, however, this metho d still leaves ro om for improvements: (i) The Nagel{Enkelmann mo del uses an optical ow constraint whichis centred in I 2 , while the pro jection matrix D in the smo othness term dep ends on I 1 .This inconsistency may create artifacts for large displacementelds. (ii) Refraining from a linearization of the optical ow constraint has the consequence that the energy functional (6) may b e nonconvex. In this case p opular algorithms such as gradient descent metho ds may get trapp ed in physically irrelevant lo cal minima. (iii) Minimizers of the energy functional (4) are not invariant under linear brightness changes of the images I 1 and I 2 . In the present pap er we will address these p oints byintro ducing three improvements into the Nagel{Enkelmann mo del: (i) We design an energy functional that consistently centers both the optical ow constraint and the smo othness constraint in the same image. (ii) We encourage convergence to the global energy minimum byemb edding the metho d into a linear scale-space framework that allows to fo cus down from coarse to ne scales in small steps. (iii) We intro duce an adaptation of the parameters C and to the dynamic range of the images such that the resulting energy functional is invariant under linear brightness rescalings. This adaptation is particularly useful in the context of our scale-space fo cusing which alters the dynamic range of the images. Applying the gradient descent metho d to our mo del leads to a coupled system of two diusion{reaction equations, for which we establish the existence of aunique solution. Interestingly, these equations can be related to anisotropic diusion ltering with a diusion tensor. We present an ecient numerical scheme that is based on a linear implicit nite dierence discretization. Afterwards, we discuss the role of the mo del parameters and demonstrate that our mo del allows very accurate recovery of optic ow elds for a large range of parameters. This is done by considering b oth synthetic image sequences, for which ground truth ow elds exist, as well as a real-world test sequence. 3
Owing to the scale-space fo cusing, our mo del is particularly suited for recovering large displacement elds. The pap er is organized as follows: In Section 2 we describ e our optical ow metho d that incorp orates the three improvements, and we showthat the Nagel{Enkelmann metho d and its mo dications are closely related to anisotropic diusion ltering. In Section 3 we present existence and uniqueness results for the nonlinear parab olic system that arises from using the gradient descent metho d for minimizing the energy functionals. Section 4 describ es an ecientnumerical discretization of this system based on a linear implicit nite dierence scheme. Section 5 claries the role of the mo del parameters, and in Section 6 we present exp erimental results on synthetic and real-world image sequences. Finally, in Section 7 we conclude with a summary. Related work. Pro esmans et al. 46, 45] studied a related approach that also disp enses with a linearization of the optical ow constraint in order to allowfor larger displacements. Their metho d, however, requires six coupled partial dierential equations and its nonlinear diusion pro cess uses a scalar-valued diusivity instead of a diusion tensor. Their discontinuity-preserving smo othing is ow-driven while ours is imagedriven. Another PDE technique that is similar in vein to the work of Pro esmans et al. is a stereo metho d by Shah 50]. Other ow-driven regularizations with discontinuitypreserving prop erties include the work of Aub ert et al. 6], Cohen 17], Deriche et al. 18], Hinterb erger 27], Kumar et al. 34], Schnorr 49], Weickert 55], and Weickert and Schnorr 57]. Related sto chastic regularization approaches have been studied by Black and Anandan 11, 12], Blanc{Feraud et al. 14], Heitz and Bouthemy 26], and Memin and Perez 36]. The image-driven anisotropic Nagel{Enkelmann approach has b een sub ject to many subsequent studies. Examples include later work by Nagel 40, 41] as well as research by Schnorr 47, 48] and Snyder 51]. A multigrid realization of this metho d has been describ ed by Enkelmann 19], and a related pyramid framework is studied by Anandan 5]. An isotropic image-driven optic ow regularization is investigated by Alvarez et al. 1]. With resp ect to emb eddings into a linear scale-space framework our metho d can b e also b e related to the optical ow approach of Florack et al. 21]. Their metho d diers from ours in that it is purely linear, applies scale selection mechanisms and do es not use discontinuity-preserving nonlinear smo othness terms. Our fo cusing strategy for avoiding to end up in irrelevantlo cal minima also resembles the graduated non-convexity (GNC) algorithms of Blake and Zisserman 13]. A preliminary version of our work has been presented at a conference 4], and a related optical ow metho d has been used by Hinterb erger 27] to generate amovie b etween two images. 2The Mo del In this section we consider three mo dications of the Nagel{Enkelmann mo del in order to improve its p erformance in the case of large displacement elds. We also discuss relations b etween this metho d and anisotropic diusion ltering. 4
2.1 Consistent Centering Wehave seen that the energy functional (4) uses an optical ow constraint and a smo othness term that are is centred in dierent images. Our exp eriments showed that this inconsistency may lead to artifacts when the displacement eld is large. As a remedy, we consider a mo died energy functional where b oth the optical ow constraint and the smo othness constraint are related to I 1 : E ( h )= Z R 2 ( I 1 ( x y ) ; I 2 ( x + u ( x y ) y + v ( x y ))) 2 dx + C Z R 2 tr ace ; r h T D ( r I 1 ) ; r h dx: (6) The asso ciated Euler-Lagrange equations are given by the PDE system C div ( D ( r I 1 ) r u )+ ; I 1 ( x ) ; I 2 ( x + h ( x )) @I 2 @x ( x + h ( x )) = 0 (7) C div ( D ( r I 1 ) r v )+ ; I 1 ( x ) ; I 2 ( x + h ( x )) @I 2 @y ( x + h ( x )) = 0 : (8) In this pap er, we are interested in solutions of the equations (7)-(8) in the case of large displacement elds and images that are not necessarily slowly varying in space. Therefore, we do not use the linearized optic ow constraint (2) in the ab ove system. 2.2 Relations to Anisotropic Diusion Filtering We obtain the solutions of the Euler{Lagrange equations (7)-(8) by calculating the asymptotic state ( t !1 )ofthe parab olic system @u @t = C div ( D ( r I 1 ) r u )+ ; I 1 ( x ) ; I 2 ( x + h ( x )) @I 2 @x ( x + h ( x )) (9) @v @t = C div ( D ( r I 1 ) r v )+ ; I 1 ( x ) ; I 2 ( x + h ( x )) @I 2 @y ( x + h ( x )) : (10) These equations do also arise when the steep est descent metho d is applied in order to minimize the energy (6). Interestingly, this coupled system of diusion{reaction equations reveals adiusion tensor which resembles the one used for edge-enhancing anisotropic diusion ltering. Indeed, D ( r I 1 ) has the eigenvectors v 1 := r I 1 and v 2 := r I ? 1 . The corresp onding eigenvalues are given by 1 ( jr I 1 j )= 2 jr I 1 j 2 +2 2 (11) 2 ( jr I 1 j )= jr I 1 j 2 + 2 jr I 1 j 2 +2 2 : (12) We observe, that 1 + 2 =1 holds indep endently of r I 1 .In the interior of ob jects we have jr I 1 j ! 0, and therefore 1 ! 1 = 2 and 2 ! 1 = 2. At ideal edges where 5
jr I 1 j ! 1 , we obtain 1 ! 0and 2 ! 1. Thus, we have isotropic b ehaviour within regions, and at image b oundaries the pro cess smo othes anisotropically along the edge. This b ehaviour is very similar to edge-enhancing anisotropic diusion ltering 53], and it is also close in spirit to the mo died mean-curvature motion considered in 3]. In this sense, one may regard the Nagel{Enkelmann metho d as an early predecessor of mo dern PDE techniques for image restoration. One structural dierence, however, should be observed: the optical ow equations (9){(10) use a temp orally constant diusion tensor, while the nonlinear diusion tensor of anisotropic diusion ltering is a function of the evolving image itself. Hence, the Nagel{Enkelmann mo del is anisotropic and space-variant, but it remains linear in its diusion part. Related linear anisotropic diusion lters have b een pioneered by Iijima in the sixties and seventies in the context of optical character recognition see 56] and the references therein. For a detailed treatment of anisotropic diusion ltering we refer to 54], an axiomatic classication of mean-curvature motion and related morphological PDEs for image analysis is presented in 2], and recent collections of pap ers on PDEbased image smo othing metho ds include 8, 16, 25, 44]. 2.3 Recovering Large Displacements by Scale-Space Fo cusing The energy functional (6) may b e nonconvex due to its data term without linearization. In this case we cannot exp ect the uniqueness of solutions of the elliptic system (7)- (8). As a consequence, the asymptotic state of the parab olic system (9)-(10), whichwe use for approximating the optical ow, dep ends on the initial data. Typically, we may exp ect that the algorithm converges to a lo cal minimizer of the energy functional (6) that is lo cated in the vicinity of the initial data. When we have small displacements in the scene, the natural choice is to take u v 0 as initialization of the ow. For large displacement elds, however, this may not work, and we need b etter initial data. To this end, we embed our metho d into a linear scale-space framework 29, 56]. Considering the problem at a coarse scale avoids that the algorithm gets trapp ed in physically irrelevant lo cal minima. The coarse-scale solution serves then as initial data for solving the problem at a ner scale. Scale fo cusing has along tradition in linear scale-space theory (see e.g. Bergholm 9] for an early approach), and in spite of the fact that some theoretical questions remain op en, it has not lost its p opularity. For more details on linear scale-space theory we refer to 20, 30, 31, 35, 52]. Using a scale-space approach enables us also to p erform a ner and more reliable scale fo cusing as it would be the case for related pyramid or multigrid approaches. We pro ceed as follows. First, weintro duce a linear scale factor in the parab olic PDE system in order to end up with @u @t = C div ( D ( r I 1 ) r u )+ ; I 1 ( x ) ; I 2 ( x + h ( x )) @I 2 @x ( x + h ( x )) (13) @v @t = C div ( D ( r I 1 ) r v )+ ; I 1 ( x ) ; I 2 ( x + h ( x )) @I 2 @y ( x + h ( x )) (14) where I 1 := G I 1 , I 2 := G I 2 , h ( x ):= ( u ( x ) v ( x )), and G I j represents the convolution of I j with a Gaussian of standard deviation . 6
The convolution with a Gaussian blends the information in the images and allows us to recover a connection b etween the ob jects in I 1 and I 2 . In our application, this global supp ort prop ertythatischaracteristic for linear diusion scale-spaces is very imp ortant. It makes them favourable over morphological scale-spaces in the sense of 2], since the latter ones cannot transp ort information b etween top ologically disconnected ob jects. We start with a large initial scale 0 . Then we compute the optical ow ( u 0 v 0 ) at scale 0 as the asymptotic state of the solution of the ab ove PDE system using as initial data u v 0 : Next, we cho ose a number of scales n < n ; 1 <:::: < 0 ,and for each scale i we compute the optical ow ( u i v i )as the asymptotic state of the ab ove PDE system with initial data ( u i ; 1 v i ; 1 ). The nal computed ow corresp onds to the smallest scale n . In accordance with the logarithmic sampling strategy in linear scale-space theory 33], we cho ose i := i 0 with some decay rate 2 (0 1). 2.4 Invariance Under Linear Greyvalue Transformations A remaining shortcoming of the mo died mo del is that the energy E ( h ) is not invariant under grey level transformation of the form ( I 1 I 2 ) ! ( kI 1 kI 2 ), where k is a constant. Therefore, the choice of the parameters dep ends strongly on the image contrast. This is esp ecially problematic when the metho d is emb edded in the scale-space fo cusing strategy, since the amount of smo othing inuences the contrast range in the regularized images G I 1 and G I 2 . We address this problem by normalizing the parameters C and in such a way that the energy E ( h ) b ecomes invariant under grey level transformation of the form ( I 1 I 2 ) ! ( kI 1 kI 2 ). We compute C and by means of two parameters and s 2 (0 1) that are calculated via C = max x ( j ( r G I 1 )( x ) j 2 ) s = Z 0 H jr G I 1 j ( z ) dz where H jr G I 1 j ( z ) represents the normalized histogram of jr G I 1 j . We name s the isotropy fraction . When s ! 0 the diusion op erator b ecomes anisotropic at all lo cations, and when s ! 1, it leads to isotropic diusion everywhere. So now C = C ( r G I 1 ), and = ( s r G I 1 ) : With this normalization of C and , the energy E ( h ) is invariant under grey level transformation of the form ( I 1 I 2 ) ! ( kI 1 kI 2 ). In practical applications of our metho d it is thus sucient to sp ecify the parameters and s instead of C and . The parameters C and are then automatically adjusted to the dynamic image range in each step of the fo cusing pro cedure. 7
3 Existence and Uniqueness of the Parab olic System In this section we show the existence and uniqueness of solutions of the parab olic system (13)-(14) where D ( r I 1 )isgiven by(5). The parameters C and can b e arbitrary p ositive real numbers. In particular, they may be determined as describ ed in the previous section. First we intro duce an abstract framework where we study the ab ove system. This framework is used to show the existence and uniqueness of the solutions afterwards. 3.1 Abstract Framework For simplicitywe assume that the images are dened on the entire space R 2 : We assume that the input images I 1 and I 2 b elong to the functional space L 2 ( R 2 ) : Let H = L 2 ( R 2 ) L 2 ( R 2 ), and let us denote by A : D ( A ) H ! H the dierential op erator dened by A ( h )= ; C 0 @ div ( D ( r I 1 ) r u ) div ( D ( r I 1 ) r v ) 1 A : If I 1 2 L 2 ( R 2 ) then I 1 2 W 1 1 ( R 2 ), so r I 1 is b ounded and the eigenvalues of the diusion tensor D ( r I 1 ) are strictly p ositive. Therefore, as C > 0, the op erator A ( h ) is a maximal monotone op erator. For more details ab out maximal monotone op erators we refer to Brezis 15]. Next, let us intro duce the function F : H ! H dened by F ( h )= ; I 1 ; I 2 ( Id + h ) r I 2 ( Id + h ) : Then the abstract evolution problem can b e written as 8 < : dh dt + Ah = F ( h ) in H , 8 t 2 0 T ] h (0) = h 0 in H: (15) Any classical solution h 2 C 1 (0 T ] H ) \ C (0 T ] D ( A )) of (15) is given by h ( t )= S ( t ) h 0 + Z t 0 S ( t ; s ) F ( h ( s )) ds (16) where f S ( t ) g t> 0 is the contraction semi-group asso ciated to the homogeneous problem. Denition. We say that h 2 C (0 T ] H ) is a generalized solution of (15) if it satises (16). 3.2 Existence and Uniqueness Result In order to prove existence and uniqueness, we have to establish a lemma rst. Lemma 1 Suppose that I 1 I 2 2 L 2 ( R 2 ) ,then F is Lipschitz-continuous, and the Lipschitz constant L depends on the functions I 1 and I 2 and on . 8
Pro of: First we note that if I 1 I 2 2 L 2 ( R 2 ), then we have in particular that I 2 2 W 1 1 ( R 2 ) and I 1 2 L 1 ( R 2 ). Let h 1 h 2 2 H . For the i -th comp onent of F ( h 1 ) ; F ( h 2 ), i =1 2, we have the following p ointwise estimate. j F i ( h 1 ) ; F i ( h 2 ) j = j ( I 1 ; I 2 ( Id + h 1 )) @ i I 2 ( Id + h 1 ) ; ( I 1 ; I 2 ( Id + h 2 )) @ i I 2 ( Id + h 2 ) j j I 2 ( Id + h 1 ) @ i I 2 ( Id + h 1 ) ; I 2 ( Id + h 2 ) @ i I 2 ( Id + h 2 ) j + j I 1 j : j @ i I 2 ( Id + h 1 ) ; @ i I 2 ( Id + h 2 ) j 1 2 j @ i ( j I 2 j 2 )( Id + h 1 ) ; @ i ( j I 2 j 2 )( Id + h 2 ) j + k I 1 k 1 : j @ i I 2 ( Id + h 1 ) ; @ i I 2 ( Id + h 2 ) j 1 2 C Lip ( @ i ( j I 2 j 2 )) : j h 1 ; h 2 j + k I 1 k 1 :C Lip ( @ i I 2 ) : j h 1 ; h 2 j 1 2 C Lip ( @ i ( j I 2 j 2 )) + k I 1 k 1 :C Lip ( @ i I 2 ) : j h 1 ; h 2 j where C Lip ( f ) denotes the Lipschitz constant of the function f . We nally deduce that k F ( h 1 ) ; F ( h 2 ) k H = k F 1 ( h 1 ) ; F 1 ( h 2 ) k L 2 + k F 2 ( h 1 ) ; F 2 ( h 2 ) k L 2 2 X i =1 1 2 C Lip ( @ i ( j I 2 j 2 )) + k I 1 k 1 :C Lip ( @ i I 2 ) : k h 1 ; h 2 k H : We conclude the proofofthe lemma by setting L = 2 X i =1 1 2 L ; @ i ( j I 2 j 2 ) + k I 1 k 1 :C Lip ( @ i I 2 ) : This shows the assertion. Now we can state the existence and uniqueness result for problem (13)-(14). Theorem 1 Suppose that I 1 I 2 2 L 2 ( R 2 ) then, for al l h 0 2 H ,there exists a unique generalized solution h ( t ) 2 C (0 1 H ) of (13)-(14). Pro of: The assumptions on I 1 and I 2 allowus to apply Lemma 1. Assume that h 1 ( t )and h 2 ( t )are solutions of (16) for initial conditions h 1 (0) and h 2 (0), then we have, using the fact that ; A is dissipative (which yields k S ( t ) f k H k f k H ), and the Lipschitz continuity of F the following estimate. k h 1 ( t ) ; h 2 ( t ) k H k h 1 (0) ; h 2 (0) k H + L Z t 0 k h 1 ( s ) ; h 2 ( s ) k H ds: 9
Figure 4: Computation of the optical ow for the taxi sequence (frames 15 and 19) with =0 : 6, s =0 : 1, 0 =10, n =0 : 8, and =0 : 95. 16
Figure 5: Vector plot of the optical ow from Figure 4. 17
Figure 6: Computation of the optical owfor the Square2 sequence with = 0 : 6, s =0 : 1, 0 = 10, n =1, and =0 : 95. heavily on this p ostpro cessing. Variational approaches with smo othness terms do not require such a p ostpro cessing step as they automatically yield ow elds with 100 % density. In Figures 6 and 7 weshow the computed optical ow for the Square2 sequence that depicts a square moving with velo city (4 = 3 4 = 3). Table 1gives a comparison with the results of Barron et al. for some classic optic ow techniques that create ow elds with 100 % density. It can be seen that our prop osed technique reveals smaller errors than these metho ds. In particular, this also shows that our three mo dications improve Nagel's metho d substantially. While the implementation of Nagel's metho d in 7] gives an angular error of 34 : 57 , our metho d reveals an error of 10 : 97 . In this example Barron et al. assume that the background moves in the same direction as the square. However, as the background is constant the displacement is not well dened in this area. If we fo cus our attention on the error of the computed ow within the interior of the square we obtain an average angular error of 0 : 85. This shows that the computed owis very accurate in the interior of the square. Next we draw our attention to the most complex synthetic test sequence from 7], the Yosemite sequence with cloudy sky.It contains displacements of up to ve pixels. Our optical ow results are shown in Figures 8 and 9, and a juxtap osition with other metho ds can be found in Table 2. Again our technique outp erforms all metho ds from 7] which yield ow elds with 100 %density. With an angular error of 5 : 53 it even 18
Figure 7: Vector plot of the optical ow from Figure 6. Table 1: Comparison between the results from 7] with 100 % density and our metho d for the Square2 sequence. Technique Aver. Error Stand. Deviat. Density Horn and Schunck (original) 47.21 14.60 100 % Horn and Schunck (mo died) 32.81 13.67 100 % Nagel 34.57 14.38 100 % Anandan (unthresholded) 31.46 18.31 100 % Singh (step 1) 49.03 21.38 100 % Singh (step 2) 46.12 18.64 100 % our metho d 10.97 9.60 100 % 19
Figure 8: Computation of the optical ow for the Yosemite sequence with = 0 : 6, s =0 : 1, 0 =5, n =1, and =0 : 95. Table 2: Comparison between the results from 7] with 100 % density and our metho d for the Yosemite sequence. Technique Aver. Error Stand. Deviat. Density Horn and Schunck (original) 31.69 31.18 100 % Horn and Schunck (mo died) 9.78 16.19 100 % Nagel 10.22 16.51 100 % Anandan (unthresholded) 13.36 15.64 100 % Uras et al. (unthresholded) 8.94 15.61 100 % Singh (step 2) 10.03 13.13 100 % our metho d 5.53 7.40 100 % 20
Figure 9: Vector plot of the optical ow from Figure 8. 21
reaches the estimation quality of typical metho ds with 30 % density, and the standard deviation of 7 : 40 is lower than the standard deviation of all metho ds that have been evaluated in 7]: the best metho d (Lucas and Kanade with 2 5 : 0) had an average angular error of 3 : 22 with a standard deviation of 8 : 92 and a density of only 8 : 7 %. In order to evaluate the robustness of our algorithm with resp ect to the choice of parameters we presen t in Table 3 the errors for the Yosemite sequence taking dierent values of the parameters. To simplify the presentation, we xed the nest scale to n := 1, and as numerical parameters we used := 10 and T := 500. These parameters are almost indep endent of the image and can therefore be set to default values. Hence, we vary only the parameters , s , and 0 in Table 3. First of all it can be seen that our metho d outp erforms all metho ds in 7] with 100 % density not only in case of optimized parameters, but also for a rather large range of parameter settings. Let us now study the parameter inuence in more detail. One imp ortant observation from Table 3 is that the decay parameter has an imp ortant inuence of the result: values around 0 : 5, as are implicitely used bytypical pyramidbased fo cusing algorithms, are by far not optimal. Aslow fo cusing with =0 : 95 gives signicantly better results. Our exp erience with other images suggests that may be xed to this value for all applications. Cho osing to o a small value for the isotropy fraction s do es hardly worsen the results, while for larger values the smo othness term b ecomes isotropic almost everywhere and approximates the Horn and Schunckscheme 28]. In order to avoid the resulting deteriorations, we prop ose to x s := 0 : 1, which means that the metho d smo othes anisotropically at 90% of all lo cations. Regarding the smo othness parameter ,our metho d app eared to be rather robust with resp ect to overand underestimations. We have thus used axed value of 0 : 6for all exp eriments in the present pap er. As already mentioned, the initial scale 0 should be chosen suchthat it covers the largest exp ected displacements. We found that overestimations are less critical than underestimations. This also conrms the use of the fo cusing strategy. Too small values increase the danger of ending up in a physically irrelevant lo cal minimum. Actually, 0 was basically the only parameter that we had to adapt in order to analyse dierent image sequences. Since it has a clear physical interpretation, this adaptation was simple. Remark. More detailed information ab out the exp eriments in this section can be found at the web site http://serdis.dis.ulpgc.es/ lalva rez/ rese arc h/de mos .In particular, some movies to illustrate the fo cusing strategy are presented. At this site we also provide a window oriented image pro cessing software named XMegaWave (see 23]) which includes the algorithm that we have develop ed in this pap er. 7 Conclusions Usually, when computer vision researchers deal with variational metho ds for optical ow calculations, they linearize the optical ow constraint. Except for those cases where the images a suciently slowly varying in space, linearization, however, do es only work for 22
Table 3: Errors for the Yosemite sequence, using dierent parameters of the algorithm smo othness init. scale isotr. fract. decay rate angul. error stand. dev. 0 s 0.4 5 0.1 0.90 5.61 7.46 0.5 " " " 5.57 7.41 0.6 " " " 5.55 7.37 0.7 " " " 5.56 7.33 1.0 " " " 5.69 7.24 0.6 1 0.1 0.90 16.83 15.23 " 2.5 " " 5.92 7.31 " 5 " " 5.55 7.37 "10 " " 5.54 7.37 "15 " " 5.81 8.45 0.6 5 0.01 0.90 5.70 7.92 " " 0.1 " 5.55 7.37 " " 0.2 " 5.70 7.31 " " 0.5 " 6.38 8.14 " " 0.8 " 7.31 9.76 " " 0.9 " 7.64 10.37 " " 0.99 " 8.04 11.21 0.6 5 0.1 0.50 7.25 7.58 " " " 0.70 6.14 7.36 " " " 0.80 5.75 7.33 " " " 0.95 5.53 7.40 " " " 0.99 5.56 7.45 23
small displacements. In this pap er we intro duced three improvements into a classical metho d by Nagel and Enkelmann where no linearization is used. We identied this metho d as two coupled linear anisotropic diusion lters with a nonlinear reaction term. Weshowed that this parab olic system is well-p osed from a mathematical viewp oint, and we presented a linear implicit nite dierence scheme for its ecientnumerical solution. In order to avoid that the algorithms converges to physically irrelevant lo cal minima, we emb edded it into a linear scale-space approach for fo cusing the solution from a coarse to anescale. A detailed quantitative analysis using test sequences with ground truth data showed the following results. The metho d can recover displacements of more than 10 pixels will go o d accuracy. It p erforms signicantly b etter than Nagel's original metho d and all other metho ds with 100 % density that are evaluated by Barron et al. 7]. The p erformance hardly deteriorated for quite a large range of parameters. This allows to use default parameter settings for many applications. We are currently investigating the use of our metho d for related matching problems such as stereo reconstruction. It is our hop e that our metho d that combines anisotropic diusion{reaction equations with linear scale-space techniques may serve as a motivation to study other combinations of linear and nonlinear scale-space approaches for solving computer vision problems. Acknowledgement. This work has b een supp orted by the Europ ean TMR network Viscosity Solutions and their Applications . References 1] L. Alvarez, J. Esclar in, M. Lefebure and J. Sanchez, A PDE model for computing the optical ow ,Pro c. XVI Congreso de Ecuaciones Diferenciales y Aplicaciones (C.E.D.Y.A. XVI, Las Palmas de Gran Canaria, Sept. 21{24, 1999), 1349{1356, 1999. 2] L. Alvarez, F. Guichard, P.-L. Lions, J.-M. Morel, Axioms and fundamental equations in image processing ,Arch. Rational Mech. Anal., Vol. 123, 199{257, 1993. 3] L. Alvarez, P.-L. Lions, J.-M. Morel, Image selective smoothing and edge detection by nonlinear diusion. II , SIAM J. Numer. Anal., Vol. 29, 845{866, 1992. 4] L. Alvarez, J. Weickert, J. Sanchez, Ascale-space approach to nonlocal optical ow calculations ,M. Nielsen, P. Johansen, O.F. Olsen, J. Weickert (Eds.), Scale-space theories in computer vision, Lecture Notes in Computer Science, Springer, Berlin, Vol. 1682, 235{246, 1999. 5] P. Anandan, A computational framework and an algorithm for the measurement of visual motion ,Int. J. Comput. Vision, Vol. 2, 283{310, 1989. 24
6] G. Aub ert, R. Deriche, P. Kornprobst, Computing optical ow via variational techniques , to app ear in SIAM J. Math. Anal. 7] J.L. Barron, D.J. Fleet, S.S. Beauchemin, Performance of optical ow techniques , Int. J. Comput. Vision, Vol. 12, 43{77, 1994. 8] M.-O. Berger, R. Deriche, I. Herlin, J. Jare, J.-M. Morel (Eds.), ICAOS '96: Images, wavelets and PDEs , Lecture Notes in Control and Information Sciences, Vol. 219, Springer, London, 1996. 9] F. Bergholm, Edge focusing ,IEEETrans. Pattern Anal. Mach. Intell., Vol. 9, 726{ 741, 1987. 10] M. Bertero, T.A. Poggio, V. Torre, Il l-posed problems in early vision ,Pro c. IEEE, Vol. 76, 869{889, 1988. 11] M.J. Black, P. Anandan, Robust dynamic motion estimation over time , Pro c. IEEE Comp. So c. Conf. on Computer Vision and Pattern Recognition (CVPR '91, Maui, June 3{6, 1991), IEEE Computer So ciety Press, Los Alamitos, 292{302, 1991. 12] M.J. Black, P. Anandan, The robust estimation of multiple motions: Parametric and piecewise smooth ow elds , Computer Vision and Image Understanding, Vol. 63, 75{104, 1996. 13] A. Blake, A. Zisserman, Visual reconstruction , MIT Press, Cambridge (Mass.), 1987. 14] L. Blanc{Feraud, M. Barlaud, T. Gaidon, Motion estimation involving discontinuities in a multiresolution scheme , Optical Engineering, Vol. 32, No. 7, 1475{1482, 1993. 15] H. Brezis, Operateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert , North Holland, Amsterdam, 1973. 16] V. Caselles, J.M. Morel, G. Sapiro, A. Tannenbaum (Eds.), Special issue on partial dierential equations and geometry-driven diusion in image processing and analysis , IEEE Trans. Image Pro c, Vol. 7, No. 3, March 1998. 17] I. Cohen, Nonlinear variational method for optical ow computation ,Pro c. Eighth Scandinavian Conf. on Image Analysis (SCIA '93, Troms$, May 25{28, 1993), Vol. 1, 523{530, 1993. 18] R. Deriche, P. Kornprobst, G. Aub ert, Optical-ow estimation while preserving its discontinuities: A variational approach , Pro c. Second Asian Conf. Computer Vision (ACCV '95, Singap ore, December 5{8, 1995), Vol. 2, 290{295, 1995. 19] W. Enkelmann, Investigation of multigrid algorithms for the estimation of optical ow elds in image sequences , Computer Vision, Graphics and Image Pro cessing, Vol. 43, 150{177, 1988. 25