Feature detection and tracking in optical flow on non-flat manifolds

Feature detection and tracking in optical flow on non-flat manifoldsSheraz Khan a ,c ,d ,⇑,Julien Lefevre b ,Habib Ammari c ,Sylvain Baillet a ,eaDepartment of Neurology,Medical College of Wisconsin,Milwaukee,USAbAix-Marseille Univ,Département d’Informatique de Luminy /CNRS,LSIS,UMR 6168,13397Marseille,France cCenter of Applied Mathematics,Ecole Polytechnique,France dDepartment of Neurology,MGH/Harvard Medical School,Boston,USA eMontreal Neurological Institute,McGill University,Montreal,Canadaa r t i c l e i n f o Article history:Received 22September 2010Available online 19September 2011Communicated by A.Shokoufandeh Keywords:Optical flowHelmholtz–Hodge decomposition Feature detection Image processingRiemannian formalism Vector fieldsa b s t r a c tOptical flow is a classical approach to estimating the velocity vector fields associated to illuminated objects traveling onto manifolds.The extraction of rotational (vortices)or curl-free (sources or sinks)fea-tures of interest from these vector fields can be obtained from their Helmholtz–Hodge decomposition (HHD).However,the applications of existing HHD techniques are limited to flat,2D domains.Here we demonstrate the extension of the HHD to vector fields defined over arbitrary surface manifolds.We pro-pose a Riemannian variational formalism,and illustrate the proposed methodology with synthetic and empirical examples of optical-flow vector field decompositions obtained on a variety of surface objects.Ó2011Elsevier B.V.All rights reserved.1.IntroductionThe optical flow in a scene is the apparent motion of objects,induced by the observed temporal variations of local brightness.Under certain conditions,it provides an adequate model for objects’displacement fields,Horn and Schunck (1981).In most applications,however,the definition of the optical flow is restricted to two-dimensional (2D),i.e.flat,surface manifolds.This hinders some potentially high-impact applications of the optical flow,e.g.the study of turbulent flow over general surfaces.We recently intro-duced a new variational framework to estimate optical-flow vector fields on non-flat surfaces using a Riemannian formulation,Lefèvre and Baillet (2008).In most applications,features of interest such as sources,sinks and traveling waves,need to be extracted from time-varying optical-flow vector fields,and be interpreted as meaningful components of the kinematics of the observed phenomena.The Helmholtz–Hodge decomposition (HHD)of vector fields has proven to be a useful technique in that respect,Chorin and Mars-den (1993).HHD proceeds by decomposing vector fields defined in 2D flat domains,Guo et al.(2005)or full 3D geometry,Tong et al.(2003),into the sum of:a curl-free component,deriving from the gradient of a scalar potential U ;a non-diverging component,deriving from the curl of a scalar (in 2D)or vector (in 3D)potential A ;an harmonic component H ,with zero Laplacian.Singularities or vanishing points in the optical-flow vector field can be detected using these components,as shown in e.g.,Guo et al.(2005),and as we shall detail in Section 3.3of this article.As already mentioned,a large number of new fields of application would ben-efit from a generic and principled approach to the HHD of vector fields on surface manifolds:e.g.,experimental fluid dynamics and turbulence (Corpetti et al.,2003;Palit et al.,2005);physiological modeling (Guo et al.,2006);structural and functional neuroimag-ing (Lefevre et al.,2009),and the compressed representation of large,distributed vector fields (Scheuermann and Tricoche,2005).It has been suggested that approximations to surface-based HHD could be achieved using polyhedral surfaces,in the local Euclidean geometry of the manifold (Polthier and Preuss,2003).This however proves to be problematic,when the vector fields are defined over highly-curved surface supports.Indeed,we have shown previously that the surface curvature plays a crucial role in the robust estimation of vector fields within the tangent bundle of the surface support (Lefèvre and Baillet,2008).Our results also showed that the non-flatness properties of the surface may cause convergence issues in the estimation process of the vector field.We note however that new results in the Euclidian domain that do not require a meshed tessellation of the surface were published recently and indicate a promising alternative to our approach (Petronetto et al.,2009).0167-8655/$-see front matter Ó2011Elsevier B.V.All rights reserved.doi:10.1016/j.patrec.2011.09.017⇑Corresponding author at:Department of Neurology,Massachusetts General Hospital,Harvard Medical School,Boston,MA 02129,USA.Tel.:+16176435634;fax:+16179485966.E-mail address:sheraz@ (S.Khan).We show in this article how the HHD can be performed onto non-flat domains and we extended the study to the extraction of basic kinetic features of interest from the HHD components.In line with our previous developments(Lefèvre and Baillet,2008),the present results were worked in the framework of differential geometry on2-Riemannian manifolds and a numerical implemen-tation using the Finite Element Method(FEM).Important concepts and notations from this framework are revisited in Section2.Sec-tion3introduces the proposed framework for HHD.The method is then featured with a variety of results in Section4.We emphasize that the software implementation of the meth-ods are available for download,as a plugin to the Brainstorm aca-demic software for electromagnetic brain imaging(http:// /brainstorm,requires Matlab and has been optimized for multithreaded computation).2.Technical backgroundThis section is a summary of the concepts and notations that are appropriate to the2-Riemannian manifold framework and its FEM implementation.More details can be found in(Do Carmo,1993; Lefèvre and Baillet,2008).2.1.Differential geometryLet us consider M,a surface(or2-Riemannian manifold) equipped by a metric or scalar product g(Á,Á)to measure distances and angles on the surface.M can be parameterized with local charts(x1,x2).Thus,it is possible to obtain a normal vector at each surface point:n p¼@@x1Â@@x2:Note that n p does not depend on the(x1,x2)parametrization.Let U:M!R be a function defined on the surface,and its dif-ferential dU:T M!R,acting on T M,the space of vectorfields.We then can write:dUðv1;v2Þ¼@U1v1þ@U2v2:Let us then define the gradient and divergence operators through duality:d UðVÞ¼gð$M U;VÞ;ZM U div M H¼ÀZMgðH;$M UÞ:We emphasize that the gradient and divergence measures are also independent of the surface parametrization.Lastly,let us define two additional functional spaces that will be useful for a theorem proposed in Section4:H1ðMÞis the space of differentiable functions on the manifold M,and C1ðMÞis the space of vectorfields,which satisfy good properties for regularity(as dis-cussed in(Druet et al.,2004;Lefèvre and Baillet,2008)).2.2.OpticalflowUnder the seminal hypothesis of the conservation of a scalar field I along streamlines,the opticalflow V is a vectorfield which satisfies:@t IþgðV;r M IÞ¼0:ð1ÞNote that the scalar product g(Á,Á)is modified by the local curva-ture of the surface domain M.The solution to Eq.(1)is not unique as long as the components of V(p,t)orthogonal to r M I are left unconstrained.This so-called‘aperture problem’has been addressed by a large number of methods using e.g.,regularization approaches.These latter may be formalized as the minimization problem of an energy functional,which both includes the regular-ity of the solution and the agreement with the model:EðVÞ¼ZM@I@tþgðV;r M IÞ2d lþkZMCðVÞd l;ð2Þwhere d l is a volume form of the manifold M.In this study,we used a regularity operator that is quadratic in the generalized gradient of the expected vectorfield:CðVÞ¼Trðt r VÁr VÞ;ð3Þwe note that the gradient of a vectorfield must be defined as the covariant derivative associated to the manifold M,so that it can be considered as an intrinsic tensor.We refer to the general con-cepts of differential geometry for more information(Do Carmo, 1993).3.Helmholtz–Hodge decomposition on a2-Riemannian manifoldWe now introduce an extended framework to perform HHD on surfaces and show that it can be applied to any vectorfield defined on a2-Riemannian manifold,M.3.1.Definitions and theoremLet usfirst define the scalar and vectorial curls by:Curl M:H1ðMÞ!C1ðMÞcurl M:C1ðMÞ!H1ðMÞCurl M A¼$M AÂn;curl M H¼div MðHÂnÞ:Note that these intrinsic expressions are not dependent on the parametrization of the surface.As established in(Polthier and Preuss,2003),we assume that M is a closed manifold,i.e.it has no boundaries.Given V,a vectorfield in C1ðMÞ,there exists functions U and A(up to an additive con-stant)in H1ðMÞand a vectorfield H in C1ðMÞsuch that:V¼$M UþCurl M AþH;ð4Þwherecurl Mð$M UÞ¼0;div MðCurl M AÞ¼0;div M H¼0;curl M H¼0:To prove this theorem,we show how U and A may be con-structed.Following classical constructions,we consider U and A as minimizers of the two following convex functionals,which are not necessarily unique:ZMk VÀ$M U k2d l;ð5ÞZMk VÀCurl M A k2d l;ð6Þwhere kÁk is the norm associated to the Riemannian metric g(Á,Á).These two functionals carry a global minimum on H1ðMÞ,which satisfies:8/2H1ðMÞ;ZMgðV;$M/Þd l¼ZMgð$M U;$M/Þd l;ð7Þ8/2H1ðMÞ;ZMgðV;Curl M/Þd l¼ZMgðCurl M A;Curl M/Þd l:ð8ÞThese two equations will play a key role in the subsequent numerical computations,as detailed below.The rest of the proof is provided in the Appendix.2048S.Khan et al./Pattern Recognition Letters32(2011)2047–20523.2.DiscretizationEqs.(7)and (8)are crucial since they pave the way to numerical implementation,when H 1ðMÞis approximated by a subspace of fi-nite elements (e.g.,continuous linear piecewise functions).We now detail the numerical implementation of Eqs.(9)and(10),which are defined over a surface tessellation c Mthat approx-imates the ideal manifold.Let N be the number of nodes in the tes-sellation (Fig.1).In the FEM implementation,we define N continuous piecewise affine functions /i ,decreasing from a value of 1at node i ,to 0at the other triangle nodes.Using these basis functions (/1,...,/n ),we can write U =(U 1,...,U n )T ,A =(A 1,...,A n )T ,and Eqs.(7)and (8)become,using array notations:Z Mg ð$M /i ;$M /j Þ i ;jU ¼Z Mg ðV ;$M /i Þi;ð9ÞZMg ðCurl M /i ;Curl M /j Þ i ;jA ¼Z Mg ðV ;Curl M /i Þ i;ð10Þwhere [Á]i ,j is a N ÂN square matrix and [Á]i is a column vector of length N .The harmonic component H of the vector field V is obtained as:H ¼V À$M U ÀCurl M A :ð11ÞThe gradient of each basis function is constant on each triangle and can be simply computed with geometrical quantities (see Lefè-vre and Baillet (2008)).Hence Eq.(9)reads:3.3.HHD feature detection as critical points of potentialsThe critical points of a vector field are often obtained from the eigenvalues of their Jacobian matrix.However,it has been shown that feature detection from critical points in global potential fields is less sensitive to noise in the data (Tong et al.,2003)and there-fore less prone to false positive detections than with methods based on the extraction of local Jacobian eigenvalues (Mann and Rockwood,2002).We therefore define critical points in the flow as local extrema of the divergence-free potential A (representing rotation)and curl-free potential U (representing divergence).A sink (respectively,a source)is thus defined as a local maxi-mum (respectively,minimum)of U .Counterclockwise (respec-tively,clockwise)vortices are featured by local maximum (respectively,minimum)of A .Further,since H is a curl-free vector field with zero divergence,traveling objects may be detected from vector elements bearing the largest magnitudes in the H vector field.4.ResultsThe proposed approach was first evaluated by detecting sources and sinks on a spherical surface bearing a rotating vector field,as shown in Fig.2.A possible real-life situation with such vector field patterns would be the air flow over a rotating spherical projectile.The HHD decomposition was applied on this vector field and vortices were identified from the critical points in A ,which was textured over the surface,as shown in Fig.2.In the case of a coun-ter-clockwise rotating vector field,the vortex was found as the maxima in the distribution of A ;for a clockwise rotating vector S.Khan et al./Pattern Recognition Letters 32(2011)2047–20522049shows a snapshot of a diverging vector field with two sources,overlaid with a texture of the values of U ,thereby revealing the correct locations of the sources in the vector field.Fig.3(c)and (d)extend further these results using vector fields with both rotating and diverging components.It is shown that features of different categories can be extracted at once from a vec-tor field.Our colormap shows sources and clockwise vortices in blue (local minima);sinks and counter-clockwise vortices in red (local maxima).In another set of simulations,a traveling source and a traveling vortex were tracked on the surface of the bunny manifold,by iden-tifying the critical points of the U and A HHD scalar fields (Fig.4(a)).Additionally,a surface patch of constant illumination was dis-placed onto the surface using the advection equation along a pre-defined vector field (Lefèvre and Baillet,2008)that defined a trajectory.This object was tracked with the vectors of largest norm in H vector field (Fig.4(b)).Fig.4shows snapshots of the tracking simulation.A black arrow has been used to trace path of sources and vortex in Fig.4(a)and moving the object in Fig.4(b).We further tested the HHD using experimental functional brain imaging data obtained from magnetoencephalography (MEG)(Baillet et al.,2001).Brain responses were recorded over 151MEG sensors following repeated somatosensory stimulation of the median nerve of the right hand (200trials),sampled at 1024Hz.Cortical currents were constrained on a surface tessella-tion of the subject’s brain containing about 50,000triangle nodes.The brain surface tessellation was obtained from the tissue seg-mentation of T1weighted magnetic resonance image data using brainVISA ( ).The cortical currents were esti-mated with a regularized minimum-norm inverse model as de-tailed in (Baillet et al.,2001).The vector field of the optical flow of these currents were computed using the method described in (Lefèvre and Baillet,2008),as implemented by us in the Brain-Storm software package (/brainstorm ).HHD was then applied to detect sources and sinks in the vector fields of the cortical currents.The source pattern at about 30ms after stimulus onset obtained from the event-related average of all trials revealed the expected primary somatosensory brain re-sponse.Fig.5shows a strong divergent pattern of the vector field over the central sulcus contralateral to the side of the stimulation.Computations to obtain all HHD components on the 50,000-node brain surface tessellation took about 10s using a conven-tional workstation with 8cores.Additional results are featured as Supplementary material .5.ConclusionIn this contribution,we developed a framework for the HHD of optical flow vector fields over 2-Riemannian manifolds.We sug-gest that features of interest regarding the kinetics supported by vector fields can be extracted and tracked through the potential components of the HHD.The required computations are simple and are made available to the academic community as part as the distribution of BrainStorm.Evaluation of this framework under real and simulated environments may open promising applications in the emerging field of multi-dimensional imaging.Appendix AWe show that the scalar potentials U and A can be obtained up to a constant.If we consider two minimizers U 1and U 2for the functional in (5)then from Eq.(7),we get:8/2H 1ðMÞ;ZMg ð$M ðU 1ÀU 2Þ;$M /Þd l ¼0:Green’s formula allows to transform this equation (given that M has no boundaries)as follows:8/2H 1ðMÞ;ZM/div M $M ðU 1ÀU 2Þd l ¼0;which corresponds to the Laplace equation:D M ðU 1ÀU 2Þ¼0;where D M ¼div M $M is the Laplace–Beltrami operator.The same approach can be followed for two minimizers A 1and A 2of the functional in Eq.(6).Hence,writing R =U 1ÀU 2and R 0=A 1ÀA 2,we obtained the following two conditions:D M R ¼0;D M R 0¼0:Multiplying by R the first equation and integrating on the man-ifold weget:rotating vector fields over a spherical surface object (shown with green arrows).(a)Counter-clockwise rotating vector field and divergence-free potentials A obtained from the HHD of the vector fields were mapped onto the surface as a colored texture.In both cases,extrema in the A map (maximum in (a),minimum in (b)).Values of A are shown in arbitrary,scaled units.(For interpretation of the reader is referred to the web version of this article.)ZMR D M R d l¼0:Applying Green’s formula yields:ZMgð$M R;$M RÞd l¼0:It follows that$M R must equal zero and R must be a constant.The same approach can be applied to R0,which must be constant.We then define the conditions on the vectorfield H.Wefirst write H¼VÀ$M UþCurl M A,and obtain:div M H¼div M VÀdiv M$M U and curl M H¼curl M VÀcurl M Curl M A;since div M Curl M A¼0and curl M$M U¼0.It follows that:8/2H1ðMÞ;ZM /div M H d l¼ZM/ðdiv M VÀdiv M$M UÞd l:Using Green’s formula and Eq.(7),we obtain: 8/2H1ðMÞ;ZM/div M H d l¼ZM gðV;$M/ÞÀZMgð$M U;$M/Þd l¼0:We have then proved that div M H¼0.The same strategy applies to obtain curl M H¼0.Appendix B.Supplementary dataSupplementary data associated with this article can be found,in the online version,at doi:10.1016/j.patrec.2011.09.017.ReferencesBaillet,S.,Mosher,J.,Leahy,R.,2001.Electromagnetic brain mapping.IEEE Signal Process.Mag.18(6),14–30.Chorin,A.,Marsden,J.,1993.A Mathematical Introduction to Fluid Mechanics.Springer.Corpetti,T.,MTmin, E.,PTrez,P.,2003.Extraction of singular points from dense motionfields:an analytic approach.J.Math.Imaging Vision19, 175–198.Do Carmo,M.,1993.Riemannian Geometry.Birkhauser.Druet,O.,Hebey, E.,Robert, F.,2004.Blow-Up Theory for Elliptic PDEs in Riemannian Geometry.Princeton University Press,Princeton,N.J(Chapter.Background Material,pp.1–12).Guo,Q.,Mandal,M.,Li,M.,2005.Efficient Hodge–Helmholtz decomposition of motionfields.Pattern Recognition Lett.26(4),493–501.Guo,Q.,Mandal,M.,Liu,G.,Kavanagh,K.,2006.Cardiac video analysis using Hodge–Helmholtzfield put.Biol.Med.36(1),1–20.Horn, B.,Schunck, B.,1981.Determining Optical Flow.Artif.Intell.17(1-3), 185–203.Lefèvre,J.,Baillet,S.,2008.Opticalflow and advection on2-Riemannian manifolds:A common framework.IEEE Trans.Pattern Anal.Machine Intell.30(6),1081–1092rs.fr/publications/2008/LB08.Lefevre,J.,Leroy,F.,Khan,S.,Dubois,J.,Huppi,P.,Baillet,S.,Mangin,J.,2009.Identification of Growth Seeds in the Neonate Brain through Surfacic Helmholtz Decomposition.In:21st Information Processing in Medical Imaging, Williamsburg.Mann,S.,Rockwood,A.,puting singularities of3D vectorfields with geometric algebra.Proc.IEEE Visual.,283–289.Palit,B.,Basu,A.,Mandal,M.,2005.Applications of the discrete Hodge Helmholtz decomposition to image and video processing.Lect.Notes Comput.Sci.3776, 497.Petronetto,F.,Paiva,A.,Lage,M.,Tavares,G.,Lopes,H.,Lewiner,T.,2009.Meshless Helmholtz–Hodge decomposition.IEEE put.Graphics, 338–349.Polthier,K.,Preuss,E.,2003.Identifying vectorfields singularities using a discrete Hodge decomposition.Visual.Math.3,113–134.Scheuermann,G.,Tricoche,X.,2005.Topological Methods for Flow Visualization.The Visualization Handbook,341–356.Tong,Y.,Lombeyda,S.,Hirani,A.,Desbrun,M.,2003.Discrete multiscale vectorfield decomposition.ACM Trans.Graphics22(3),445–452.。