The variational particle-mesh method for matching curves
aa r X i v : . [ m a t h . NA ] J u l The variational particle-mesh method for matchingcurves
C J Cotter
Department of Aeronautics, Imperial College London, SW7 2AZE-mail: [email protected]
Abstract.
Diffeomorphic matching (only one of several names for this technique)is a technique for non-rigid registration of curves and surfaces in which the curve orsurface is embedded in the flow of a time-series of vector fields. One seeks the flowbetween two topologically-equivalent curves or surfaces which minimises some metricdefined on the vector fields, i.e. the flow closest to the identity in some sense.In this paper, we describe a new particle-mesh discretisation for the evolution ofthe geodesic flow and the embedded shape. Particle-mesh algorithms are very naturalfor this problem because Lagrangian particles (particles moving with the flow) canrepresent the movement of the shape whereas the vector field is Eulerian and hencebest represented on a static mesh. We explain the derivation of the method, andprove conservation properties: the discrete method has a set of conserved momentacorresponding to the particle-relabelling symmetry which converge to conservedquantities in the continuous problem. We also introduce a new discretisation for thegeometric current matching condition of (Vaillant and Glaunes, 2005). We illustratethe method and the derived properties with numerical examples.
Keywords : symplectic integrators, diffeomorphic shape matching, geodesics on thediffeomorphism group, EPDiffAMS classification scheme numbers: 65P10,65N21PACS numbers: 87.57.N-
Submitted to:
J. Phys. A: Math. Gen. he variational particle-mesh method For Darryl Holm on the occasion of his 60th birthday
1. Introduction
Diffeomorphic matching is a numerical framework for quantifying the differencesbetween geometric information (such as curves, surfaces, images or vector fields) usingdeformations from one geometric object to another. The geometric object is embeddedin the flow generated by a time-series of vector fields with a chosen norm (such asthe H -norm) defining the distance along the path between geometric objects. Thecomputational task is to calculate the velocity fields which minimise this distancesuch that one geometric object is mapped to another. This framework was originallyintroduced in a series of papers including [GM98, CY01, MY01, MTY03] being extendedto match distributions (which can model curves and surfaces) in [GTY04], and tomatch geometric currents (also for modelling curves and surfaces) in [VG05]. Variousnumerical approaches have been proposed for solving the optimisation problem, eitherby optimising the functional directly (with an extra term to penalise flows which do notmap close to the target shape) as described in [BMTY05], or by solving the equations ofmotion and shooting for a match between shapes by adjusting the initial conditions as in[TMT02, MM06a, MMS06]. The main challenge remains to find a numerical approachwhich is accurate and efficient, since the problem of computing the shortest path is ahigh-dimensional optimisation problem.In this paper we introduce a new numerical discretisation for the diffeomorphicmatching problem in the context of matching (although it can also be used for matchingsurfaces, images and vector fields). This method uses a similar approach to theHamiltonian Particle-Mesh (HPM) method [FGR02], with the difference being thatHPM uses the particle-mesh discretisation to interpolate density from the particles to themesh, whereas in this application the particle-mesh discretisation is used to interpolatemomentum (which takes a central role in the diffeomorphic matching framework).In section 2 we give a review of the diffeomorphic matching approach applied tocurves in the plane, and establish the notation which will be used in the other sections.We also discuss the role of momentum, the conditions used to establish whether thecurves have been matched, and the implications of the particle-relabelling symmetrysatisfied by the equations of motion, as well as the connection with EPDiff and theCamassa-Holm equation. In section 3 we introduce the particle-mesh discretisation anddiscuss the discrete symmetries and conservation laws, as well as a discretisation ofthe current matching condition and a description of solution methods. In section 5 weillustrate the properties of the numerical method applied to computing the shortest pathbetween two test shapes. Finally, in section 6 we give a summary and outlook. he variational particle-mesh method
2. Diffeomorphic matching of embedded curves
In this section we describe the problem of matching one embedded curve onto anotherusing diffeomorphisms. For simplicity we shall focus on simple closed curves in the planealthough the approach is easily generalised to other structures.
We take two simple curves embedded in R , C A and C B , with parameterisations Q A ( s ) , Q B ( s ) , s ∈ [0 , π ) . Note that the curves can equally well be represented by reparameterisations of the curves i.e. by Q A ( η ( s )) , Q B ( η ( s )) , s ∈ [0 , π ) , where η ( s ) is a diffeomorphism of the circle. Our aim is to find a matching processwhich is independent of such reparameterisations. If we take a vector field which defines a fluid flow u ( x , t ), and embed a curve in thatflow, then the curve satisfies ∂∂t Q ( t ; s ) = u ( Q ( t ; s ) , t ) , (1) i.e. each point of the curve moves with the vector evaluated at that point.The aim of the calculation is to search amongst time-series of vector fields u ( x , t ), t ∈ [0 ,
1] such that (1) is satisfied, with the boundary conditions Q (0; s ) = Q A ( s ) , Q (1; s ) = Q B ( η ( s )) , for some (unspecified) reparameterisation η . If these conditions are satisfied then wesay that Q A is matched onto Q B by the vector field time series u ( x , t ). We choose a norm for vector fields, such as the H nα norm defined by k u k H nα = Z R u ( x ) · (1 − α ∇ ) n u ( x )dVol( x )Given a time series u ( x , t ), we can measure the total amount of deformation in the flowgenerated by u by the integral S ( u ) = Z k u k d t. (2)We wish to find the flow that maps C A to C B which is “nearest to the identity” withrespect to the choice of norm. This leads to the following optimal control problem[GTY06]: he variational particle-mesh method Definition 1 (Optimisation problem for curve matching)
Let Q A , Q B : S → R be parameterisations of two curves C A , C B in the plane, and let k · k be a norm forvector fields in the plane. Then the distance between curve C A and C B is defined to bethe minimum over all vector fields u of the functional Z k u k d t subject to the following constraints: • Dynamic constraint : ∂∂t Q ( t ; s ) = u ( Q ( t ; s ) , t ) , ∀ t ∈ [0 , . • Matching conditions : Q (0; s ) = Q A ( s ) , Q (1; s ) = Q B ( η ( s )) , where η is somereparameterisation of the circle.2.4. Momentum We can enforce the dynamic constraint by introducing Lagrange multipliers P ( t ; s ), sothat the optimisation problem becomes δ Z k u k + Z πs =0 P · ( ˙ Q − u ( Q ))d s d t = 0 , (3)subject to the matching conditions give above. The Euler-Lagrange equations then give δlδ u = Z πs =0 P δ ( x − Q )d s, (4)˙ Q = u ( Q ) , (5)˙ P = − P · ∇ u ( Q ) , (6)where l = k u k /
2. For example, if we choose the H nα -norm then δlδ u = (1 − α ∇ ) n u . We see that optimal velocity fields take the form Z πs =0 P G ( x − Q )d s, (7)where G is the Green’s function associated with the chosen norm, e.g. for the H nα -normit is the Green’s function for the operator (1 − α ∇ ) n .Equation (4) allows us to write u as a function of P and Q , leading us to noticethat equations (5-6) are canonically Hamiltonian, with Hamiltonian function given by H = l ( u ( P , Q )) = 12 k u ( P , Q ) k , i.e. , half the square of the norm of the velocity field written as a function of P and Q . It is for this reason that we refer to the Lagrange multiplier P as the momentum associated with the curve. This Hamiltonian structure arises from the fact that theequations have been derived from a variational principle defined by the extreme pointsof the functional. Since the Hamiltonian for this system is time-independent, it isconserved along the trajectory, and hence the norm of the velocity is also conserved. he variational particle-mesh method There are a number of different possible ways to pose the matching conditionmathematically. One matching condition which is invariant under reparameterisationsis based on defining a singular vector field v Q ( x ) = Z πs =0 ∂ Q ∂s δ ( x − Q )d s. (8)We define a functional for these singular vector fields f (cid:2) v Q (cid:3) = Z v Q · K ∗ v Q dVol( x )= Z πs =0 ∂ Q ∂s · Z πs ′ =0 ∂ Q ∂s ′ K ( Q ( s ) − Q ( s ′ ))d s d s ′ . (9)where K is some smooth kernel function. When Q matches Q B then f [ v Q − v Q B ]vanishes. This is called the current matching condition [VG05].This is a weaker condition that setting the value of Q ( s ) for each s at time t = 1, andhence we need to add another boundary condition to get a unique solution to equations(4-6). In [MTY03] it was shown that the solution which minimises the functional is theone with P initially normal to the curve i.e. P · ∂ Q ∂s = 0 , at t = 0 , and this extra boundary condition means that equations (4-6) have a unique solution. The optimisation problem in definition 1 is invariant under symmetries given byreparameterisations of the curve η . Noether’s theorem tells us that the equations ofmotion (4-6) have conserved quantities which are generated by these symmetries.To compute the conserved quantities we compute the infinitesimal generators ofthese symmetries δ Q ( s ; t ) = ∂ Q ∂s · ξ ( s )where ξ ( s ) is a vector field on the circle. The cotangent lift of this infinitesimal symmetryis δ ( Q ( s ; t ) , P ( s ; t )) = (cid:18) ∂ Q ∂s · ξ ( s ) , − P ( s ; t ) · ∂ Q ∂s · ξ ( s ) (cid:19) which is a Hamiltonian flow generated by the Hamiltonian functional h ξ = Z πs =0 P ( s ) · ∂ Q ∂s · ξ ( s )d s. By Noether’s theorem, h ξ is conserved along solutions of (4)-(6), and since ξ ( s ) isarbitrary, we see that P ( s ) · ∂ Q ∂s he variational particle-mesh method s . In particular, this means that if themomentum is initially normal to the curve i.e. P ( s ; 0) · ∂ Q ∂s ( s ; 0) = 0 , then the momentum is normal to the curve for all values of t along the solution.Therefore, optimal velocity fields take the form of equation (7) with the added constraintthat the momentum P is normal to the curve. As described in [CH, CHH07], because the dynamic constraint is given by a Lie algebraaction of velocity fields on embedded curves, it is possible to eliminate P and Q bytaking the time derivative of equation (4) and making use of equations (5-6). This leadsto a PDE defined on the whole of R given by m t + ∇ · ( um ) + ( ∇ u ) T m = 0 , m = δlδ u . This is the Euler-Poinc´are equation on the diffeomorphism group, abbreviated as EPDiff[HM04, HMR98], which is the equation for geodesic flow on the diffeomorphism group.Although we do not explicitly solve EPDiff during the computation of the optimal flow,it is useful to understand the computed solutions as singular solutions of EPDiff. In onedimension and with the H α -norm, EPDiff becomes the Camassa-Holm equation [CH93]which is completely integrable with singular soliton solutions.
3. Particle-mesh discretisation
We take a discrete mechanics and optimal control [JMOB05] approach by applying thediscretisation directly to the functional (3) and deriving the resulting equations. Thediscretisation used is a particle-mesh method with: • the vector fields discretised on a fixed, finite mesh , and • the curve discretised as a set of moving points.The principle benefits of this discretisation are that different norms can be defined onthe mesh without the need to calculate Green’s functions. With large numbers of pointsthe method becomes very efficient as particle momenta can be interpolated to the mesh,then the norm operator can be inverted, then the velocity values can be interpolatedback to the particle positions. We take a fixed set of points on a mesh { x k } n g k =1 (for the numerical examples in thispaper we used an equispaced square mesh) and give each mesh point x k a vector u k . he variational particle-mesh method { u k } n g k =1 of vectors to a general point x in the plane usinga linear interpolation u ( x ) = n g X k =1 u k ψ k ( x ) . (10)For the examples in this paper we set ψ k to be a tensor product of cubic B-splinefunctions centred on x = x k . Remark 2
The set of vectors on the mesh generate a finite dimensional subspace ofthe infinite dimensional space of vector fields. However, the finite dimensional subspaceis not closed under the Lie bracket for vector fields which means we will not be able toobtain a discrete EPDiff equation from the discrete equations of motion by eliminating Q and P as in section (2.7). Once we have this discrete representation of the vector fields we can define adiscretised norm using standard mesh methods (such as finite difference, finite volume etc. ). For the the examples in this paper we took periodic boundary conditions, and useda spectral discretisation of the (1 − α ∇ ) n operator using discrete Fourier transforms,and a simple Riemann sum for the integration (which gives spectral accuracy in thiscase). Since the Green’s functions decay exponentially over the lengthscale α , theboundary conditions should not affect the solution as long as the curve is sufficientlyfar from the boundary. Other boundary conditions are possible if, for example, a finite-difference method is used to discretise the norm operator. We replace the parameterised curve Q ( s ) by a finite set of points in the plane { Q β } n p β =1 .The equation (1) gets replaced by the semi-discrete (continuous time/discrete space)equation ˙ Q β = n g X k =1 u k ψ k ( Q β ) . (11)We also have to choose a discretisation of integration around the loop; again for theexamples in this paper we use a Riemann sum. After the particle-mesh discretisation we obtain a semi-discrete functional by integratingthe discrete vector field norm from t = 0 to t = 1 and introducing Lagrange multipliersto enforce the constraint (11) to get Z k u k g + X β P β · ˙ Q β − n g X k =1 u k ψ k ( Q β ) ! d t. (12) he variational particle-mesh method ∂∂ u k k u k g = X β P β ψ k ( Q β ) , (13)˙ Q β = n g X k =1 u k ψ k ( Q β ) , (14)˙ P β = − P β · X k u k ∇ ψ k ( Q β ) . (15)Once again, equation (13) allows us to write u as a function of P and Q , leading us tonotice that equations (14-15) are canonically Hamiltonian, with Hamiltonian functiongiven by k u ( P , Q ) k g / Equation (11) can be discretised using any one-step method, such as a Runge-Kuttamethod. For simplicity, in this paper we use the forward Euler discretisation, but thederivation of the equations is very similar for other methods. The equation becomes Q n +1 β = Q nβ + ∆ t X k u n +1 ψ k ( Q n ) . The discrete functional becomes N X n =1 ∆ t k u k g + X β P n +1 β · Q n +1 β − Q nβ − ∆ t n g X k =1 u n +1 k ψ k ( Q nβ ) !! , (16)where ∆ t = 1 /N , and the discrete Euler-Lagrange equations are ∂∂ u k k u n +1 k g = X β P n +1 β ψ k ( Q nβ ) , (17) Q n +1 β = Q nβ + ∆ t X k u n +1 k ψ k ( Q nβ ) , (18) P n +1 β = P nβ − ∆ t P n +1 · X k u n +1 k ∇ ψ k ( Q nβ ) . (19)These equations provide a variational integrator [LMOW03] for the semi-discreteequations (13-15). Hence, they give a symplectic integrator (see [LR05] for a survey ofthese methods) for the semi-discrete equations in Hamiltonian form. In this particularcase we obtain the first-order symplectic Euler method; higher-order partitioned Runge-Kutta methods can be obtained by discretising (11) using any Runge-Kutta method. In this section we discuss the properties of applying a variational integrator to equations(13)-(15), and their possible benefits for the problem of matching curves.The properties that make variational integrators the best choice for longintegrations are: he variational particle-mesh method • Modified Hamiltonian : if the Hamiltonian H ( P , Q ) is analytic then it is possibleto use backward error analysis to find a modified Hamiltonian˜ H = H + ∆ t p ∆ H ( P , Q ; ∆ t )which is conserved over times | t | ≤ c exp c / ∆ t . See [LR05] for a summary and references. • Conserved momenta : if the continuous-time Euler-Lagrange equations haveconserved momenta associated with symmetries of the Lagrangian, then discreteEuler-Lagrange equations also conserve these momenta [LMOW03].
In the problem of matching curves, the Hamiltonian is half the squared norm for vectorfields, which is conserved along trajectories as described in section 2.4. If a symplecticintegrator is used then backward error analysis guarantees that the Hamiltonian willbe approximately conserved for long times (although it should be noted that no theoryexists for the piecewise-cubic functions used in the examples of this paper). It remainsan open question as to whether the conservation of the Hamiltonian is important forthese problems on finite-time intervals, since they are on short time intervals and sothe value of the Hamiltonian will also be approximately preserved by variable step-sizemethods using error control.The main cost of using variational integrators is that they are implicit for this typeof system where the Hamiltonian is a function of P and Q , but during the optimisationalgorithm we are able to store values along the trajectory from previous integrationswhich can be used as initial guesses for solving the implicit equations using Newtoniteration. The variational integrator allows large timesteps to be taken in that case(although of course this is the case for any implicit method, symplectic or otherwise). As described in section 2.6, the equations (4-6) have a symmetry underreparameterisation of the curve which has an associated conserved momentum P · ∂ Q /∂s . This means that for the optimal solution, P remains normal to the curve alongthe trajectory. Since we have discretised the curve we have broken that symmetry, butas we shall show in this section, it is still possible to recover a sense in which P is normalto the curve along the solution.Equation (10) defines a vector field over the whole space parameterised by thevectors { u k } n g k =1 associated with the mesh points. This means that we can choose aparameterised curve Q ( s ) which passes through our discrete points { Q β } n p β =1 in sequenceand follow its evolution along the flow by solving ∂∂t Q ( s ) = X k u k ψ k ( Q ( s )) . (20) he variational particle-mesh method ∂∂t g ( x , t ) = X k u k ψ k ( g ( x , t )) , where g ( x , t ) is the flow map taking points from their position at time 0 to their positionat time t ( x plays the role of a coordinate on “label space” here), and we can follow theJacobian of this map ∂∂t ∂ g ( x , t ) ∂ x = X k u k ∇ ψ k ( g ( x , t )) · ∂ g ( x , t ) ∂ x . In particular, we can evaluate this equation at each of our discrete points Q β :˙ J β = X k u k ∇ ψ k ( Q β ) · J β , J β = ∂ g ∂ x ( Q β (0) , t ) , J β (0) = Id . (21)A variation δ Q β (0) in the initial conditions Q β (0) leads to a variation in the entiretrajectory given by δ Q β ( t ) = J β ( t ) δ Q β (0) . (22)This variation generates a symmetry of the equations, as shown in the following lemma. Lemma 3
The infinitesimal transformation given by equations (21)-(22) together with δ P β = 0 , β = 1 , . . . , n p , δ u k = 0 , k = 1 , . . . , n k , is a symmetry of equations (13)-(15).Proof. Since the equations have been derived from an action principle, we simply needto show that the infinitesimal transformation causes the action (12) to vanish. If weapply the transformation to the integrand (the Lagrangian) then we obtain δL = δ k u k g + X β P β · ˙ Q β − n g X k =1 u k ψ k ( Q β ) !! , = X β P β · δ ˙ Q β − n g X k =1 u k ∇ ψ k ( Q β ) · δ Q β ! , = X β P β · ˙ J β · δ Q β (0) − n g X k =1 u k ∇ ψ k ( Q β ) · J β δ Q β (0) ! , = X β P β · X k u k ∇ ψ k ( Q β ) · J β · δ Q β (0) − n g X k =1 u k ∇ ψ k ( Q β ) · J β δ Q β (0) ! , = 0 , and hence the result. (cid:3) This symmetry of the equations has an associated conserved momenta followingNoether’s theorem, as described in the following proposition. he variational particle-mesh method Proposition 4
For each β = 1 , . . . , n p , the quantity J Tβ P β is conserved.Proof. Applying the symmetry to the action principle and substituting the equations ofmotion (13)-(15) gives0 = δ Z k u k g + X β P β · ˙ Q β − n g X k =1 u k ψ k ( Q β ) ! d t = Z X β ˙ Q β − n g X k =1 u k ψ k ( Q β ) ! · δ Q β + dd t X β P β J β · δ Q β (0) ! d t = X β ( P β ( t ) J β ( t ) − P β ( t ) J β ( t )) · δ Q β (0) . We obtain the result since δ Q β (0) is an arbitrary vector. (cid:3) As described in [LMOW03], these conservation laws satisfied by the semi-discreteequations will be preserved by numerical methods provided that they are derived froma discrete variational principle i.e. following the framework described in section 3.4.If we choose a continuous curve through our set of points that moves with the flowaccording to equation (20), then an approximation to ∂ Q ( s β ) /∂s · d s is given at time t by J β ∆ Q β (0) , β = 1 , . . . , n p , where ∆ Q β (0) is an approximation to ∂ Q ( s β ) /∂s · d s at time t = 0. This leads to thefollowing corollary: Corollary 5 If P β is chosen to be normal to ∆ Q β (0) at time t = 0 for all β i.e. P isinitially normal to the curve, then P is normal to J β ∆ Q β (0) for all t along the flow, i.e. P stays normal to the curve.Proof. For each β = 1 , . . . , n p , the component of P tangential to the shape is P β ( t ) · J β ( t )∆ Q β (0) = P β (0) · J β (0)∆ Q β (0)= P β (0) · Id∆ Q β (0)= P β (0)∆ Q β (0) = 0 . (cid:3) If the equations of motion for Q and P are discretised using a variational integratorthen this property is preserved, provided that the discrete equation for the evolution of J is defined as the gradient of the evolution equation for Q , since then it generates asymmetry of the discrete variational principle just as in the continuous time case.For example, the discrete equation for the gradient of the time-discrete flow for Q obtained from the symplectic Euler method is J n +1 β = J nβ + ∆ t X k u n +1 k ∇ ψ ( Q nβ ) J nβ . (23) he variational particle-mesh method Proposition 6
For each β = 1 , . . . , n p , the quantity J Tβ P β is conserved when theequations are integrated using the symplectic Euler method, and J is obtained fromequation 23.Proof. δA = δ N X n =1
12 ∆ t k u k g + X β P n +1 β · Q n +1 β − Q nβ − ∆ t n g X k =1 u n +1 k ψ k ( Q nβ ) !! , = N X n =1 X β P n +1 β · δ Q n +1 β − δ Q nβ − n g X k =1 u n +1 k ∇ ψ k ( Q nβ ) · δ Q nβ !! , = N X n =1 X β P n +1 β · ( J n +1 β − J nβ ) · δ Q β − n g X k =1 u n +1 k ∇ ψ k ( Q nβ ) · J nβ δ Q β ! , = X β P n +1 β · X k u n +1 k ∇ ψ k ( Q nβ ) · J nβ · δ Q β − n g X k =1 u n +1 k ∇ ψ k ( Q nβ ) · J nβ δ Q β ! , = 0 , Applying the symmetry to the discrete action principle 16 and substituting the equationsof motion (13)-(15) gives0 = X β (cid:0) P n +1 β J n +1 β − P nβ J nβ (cid:1) · δ Q β . We obtain the result since δ Q β is an arbitrary vector. (cid:3) We can verify the discrete conservation of P β J β directly since P nβ J nβ = P n +1 β (Id + X k u n +1 k ∇ ψ k ( Q nβ )) J nβ = P n +1 β J n +1 β . Similar results can be obtained when higher-order variational integrators are used bydiscretising the Q equation using a Runge-Kutta method.As discussed in section 2.5, the norm-minimising solution has momentum normalto the shape along the whole trajectory. Good preservation of the tangential componentof momentum is important because it means that it is only necessary to constrain themomentum to be normal to the shape in the initial conditions when a shooting algorithmis used (discussed in section 4). We can also apply a particle-mesh discretisation to the matching condition described insection 2.5. The particle-mesh representation of equation (8) is v Q k = X β ∆ Q β ψ k ( Q β ) , he variational particle-mesh method Q β is an approximation to ∂ Q ( s ) /∂s d s at s = s β . For example, a simple finitedifference approximation gives v Q k = X β ( Q β − Q β − ) ψ k ( Q β ) . This is the approximation which we will use in the computed numerical examples usedin this paper.Once the singular vector field is evaluated on the mesh we can apply standard meshdiscretisation methods to compute the functional (9):ˆ f [ v Q ] = X kl K kl v Q k v Q l where K kl is a discretisation of the kernel operator K . For the computed numericalexamples used in this paper, the kernel operator used was the inverse operator(1 − α ∇ ) − discretised using discrete Fourier transforms.This results in a numerical discretisation of the functional used for the currentmatching condition. After numerical discretisation the functional will not have aminimum at zero any more, and we must aim to minimise this functional rather thanfind the zero to achieve matching. In [VG05], a mesh-free method was applied to the current matching problem. Themain cost in this type of method is in summing up Green’s functions on each of thepoints on the curve to calculate their velocity, and the value of the matching condition.They employed a fast multipoles algorithm [GS91] which has a computational cost ofsize O ( N log N ) (where N is the total number of particles), although the multiplicativeconstant in the scaling can be relatively big, requiring larger N before benefits are seen.For a particle-mesh method, the cost of evaluating the momentum on the grid is O ( N ), and the cost of the FFT is O ( M log M ) (where M is the number of rows ofgridpoints i.e. the total number of grid points is M ), although one could reduce thisby discretising the operator used in the norm for velocity ( e.g. the Helmholtz operator inthe examples given here) and just performing a few iterations of a method such as Jacobi,resulting in an O ( M ) cost. In the case where one is matching dense information (imagesor measures, for example), then typically N = cM (with c >
1) and the particle-meshapproach produces a method which is competitive with fast multipoles methods. Forthe case of matching curves, typically N = cM (with c >
1) and a good fast multipoleimplementation will be more efficient. However, any mesh-based operator inversion iseasily parallelised using standard methods (as is the particle-mesh operation), and sothis could be a strength of the particle-mesh approach in the curve case. There are alsoother advantages to using a mesh, for example it is very easy to modify the operator toinvestigate the behaviour with different kernels. he variational particle-mesh method
4. Solution methods
As discussed in section 3.8, numerical discretisation of the matching condition meansthat we cannot require that the matching functional vanishes, and so we must take thisinto account in the solution approach. There are two main approaches to obtaining anumerical approximation to the optimal flow between embedded curves:(i)
Inexact matching : In this approach, suggested in [MTY03], we “softly” enforcethe matching condition by adding a penalty term to the functional (2): S ( u ) = Z k u k d t + 1 σ f [ v Q ] . (24)One then attempts to minimise this functional directly by applying a descentalgorithm (such as nonlinear conjugate gradients) by varying the time series ofvector fields and computing the implicit gradient of the functional f . An alternativemethod is to introduce the dynamical constraint using the Lagrange multipliers P as in section 2.4, and to solve the resulting Euler-Lagrange equations (4-6)(modified to accomodate the penalty functional) using Newton iteration. Thisapproach becomes attractive when it is possible to find a good initial guess at thesolution e.g. by specifying a path of curves between C A and C B and approximatelysolving equation (6) to get an initial guess for P given the initial guess for Q .It is worth mentioning that the minimisation of the functional 24 becomes ill-conditioned as σ → method ofmultipliers , described in [Ber82], introduces Lagrange multiplier variables alongwith the penalty parameter σ and allows a reduction of the error in matchingwithout making σ arbitrarily large.(ii) Minimisation by shooting : In this approach, advocated in [MM06b], we seekinitial conditions for P such that the functional 9 is minimised. The gradient of thefunctional with respect to the initial conditions for P is computed by solving theEuler-Lagrange equations (4-6) from t = 0 to t = 1 with those initial conditions,computing the gradients of f [ u Q N ] with respect to Q Nβ , and propagating thosegradients back to P β using the adjoint equations (see [Gun03], for example). Theproblem can then be solved using a nonlinear gradient algorithm such as nonlinearconjugate gradients. As described in section 2.5, P must be constrained to benormal to the curve in order to obtain the optimal path.This is the approach which we used in computing the numerical examples.
5. Numerical examples
In this section we show a computed curve-matching calculation of two curves in theplane. The curve was parameterised using 420 points, with a square 128 ×
128 mesh ofsize 2 π × π . The velocity norm used was the H α -norm with α = 0 .
4, discretised on he variational particle-mesh method Figure 1.
Test curves used for the example matching computation.
Figure 2.
Snapshots of the optimal path between two test curves computed using theparticle-mesh discretisation. The computed curve is plotted with a continuous line,superimposed on the target curve plotted with a dashed line. Top row: the curve at t = [0 , . , . t = [0 . , . , the mesh using discrete Fourier transform, and the kernel used for the current matchingwas the Green’s function of the (1 − α ∇ ) − operator with α = 0 . he variational particle-mesh method Figure 3.
Plots illustrating the optimal deformation map which maps between ourtwo test curves. The numerical solution specifies a vector field which is definedeverywhere, which can be used to transport other points which are not on the curve.This calculation was performed on a set of points on equispaced grid lines to showhow space is being deformed around the curve. Top: the grid lines and curve beforedeformation. Bottom: the grid lines and curve after deformation. Note that the flowvector fields are far from divergence-free, as can be seen by inspecting the areas of thesquares on the deformed grid. Note also that the optimal flow only deforms space nearto the curve. he variational particle-mesh method Figure 4.
Plots showing the evolution of the momentum P along the optimal flow.The plots are taken from the same snapshots as in figure 2, with vectors showing thedirection and magnitude of P on the curve. Note that P remains normal to the curvethroughout. Figure 5.
Plots showing J β ∆ Q β around the curve which is an approximation to ∂ Q /∂s d s as described in section 3.5. This approximation remains normal to the curvethroughout, with the increase in magnitudes showing regions of stretching in the flow. he variational particle-mesh method + ] optimize.fmin ncg routine which applies the Newton conjugate gradients algorithmwith the hessian computed by applying finite differences to the gradient. The algorithmwas run out until the functional was reduced to a value of 1 . × − (with an initialvalue of 0 . J β ∆ Q β (0) which isan approximation to ∂ Q /∂s d s . The quantity P β P β · J β ∆ Q β remains within round-offerror of zero throughout, confirming the results of section 3.5.
6. Summary and outlook
In this paper we introduced a new particle-mesh discretisation for diffeomorphicmatching of curves, which can also be used for matching images. In this method thevector fields used to transport the curves are represented on a fixed mesh, whilst thecurves themselves are represented by a finite set of moving points. Since the discreteequations arise from discretising an action principle, they are variational integratorswhich have many favourable properties. As discussed in [MM06a], whilst the benefitsof variational integrators have been established for long integrations such as those forcelestial mechanics or molecular dynamics, the benefits for short time optimal controlproblems such as the curve matching problem discussed in this paper are not so clear.There are a number of drawbacks and benefits which need to be investigated withfurther testing. In this paper we showed that the variational integrators that arise fromthe particle-mesh discretisation have a discrete form of the momentum conservationlaw which leads to the curve momentum remaining normal to the curve throughoutthe computed trajectory between shapes. We also noted that the integrators have amodified Hamiltonian which is conserved over long times (but not exponentially longsince the compactly supported basis functions ψ k ( x ) are not analytic) which can beinterpreted as a modified metric for the discrete equations. We then showed illustrativeexamples obtained from the calculation of the optimal trajectory between two closedcurves in the plane.The main computational challenge for diffeomorphic matching remains the design ofefficient algorithms to obtain optimal trajectories for large datasets. Of the two solutionmethods described in section 4, the inexact matching method applied to the particle-mesh discretisation may be best applied in parallel by using a Newton-Krylov method he variational particle-mesh method [Ber82] D. P. Bertsekas. Constrained Optimization and Lagrange Multiplier Methods . AcademicPress, 1982.[BMTY05] M. F. Beg, M. I. Miller, A. Trouv´e, and L. Younes. Computing large deformation metricmappings via geodesics flows of diffeomorphisms.
International Journal Of ComputerVision , 61(2):139–157, 2005.[CH] C. J. Cotter and D. D. Holm. Continuous and discrete Clebsch variational principles. Toappear in Foundations of Computational Mechanics, 2008.[CH93] R. Camassa and D. D. Holm. An integrable shallow water equation with peaked solitons.
Phys. Rev. Lett. , 71:1661–1664, 1993.[CH08] C. J. Cotter and D. D. Holm. Discrete momentum maps for lattice EPDiff. To appear inHandbook of Numerical Analysis, 2008.[CHH07] C. J. Cotter, D. D. Holm, and P. E. Hydon. Multisymplectic formulation of fluid dynamicsusing the inverse map.
Proc. Roy. Soc. A , 463(2086):2671–2687, 2007.[Cot05] C. J. Cotter. A general approach for producing Hamiltonian numerical schemes for fluidequations. arXiv:math.NA/0501468, 2005.[CY01] V. Camion and L. Younes. Geodesic interpolating splines. In M Figueiredo, J Zerubia,and K Jain, A, editors,
EMMCVPR 2001 , volume 2134 of
Lecture notes in computersciences . Springer, 2001.[FGR02] J. Frank, G. Gottwald, and S. Reich. A Hamiltonian particle-mesh method for the rotatingshallow-water equations. In
Lecture Notes in Computational Science and Engineering ,volume 26, pages 131–142. Springer-Verlag, 2002.[GM98] U. Grenander and M. I. Miller. Computational anatomy: An emerging discipline.
Quarterly of Applied Mathematics , LVI(4):617–694, 1998.[GS91] L. Greengard and J. Strain. The fast Gauss transform.
SIAM Journal of Scientific he variational particle-mesh method Statistical Computing , 12(79–94), 1991.[GTY04] J. Glaunes, A. Trouv´e, and L. Younes. Diffeomorphic matching of distributions: A newapproach for unlabelled point-sets and sub-manifolds matching. In
IEEE ComputerSociety Conference on Computer Vision and Pattern Recognition , volume 2, pages 712–718, 2004.[GTY06] J. Glaunes, A. Trouve, and L. Younes. Modeling planar shape variation via hamiltonianflows of curves. In H. Krim and A. Yezzi Jr., editors,
Statistics and Analysis of Shapes .Birkh¨auser, 2006.[Gun03] M. D. Gunzburger.
Perspectives in Flow Control and Optimization . Advances in designand control. SIAM, Philadelphia, USA, 2003.[HM04] D. D. Holm and J. E. Marsden. Momentum maps & measure valued solutions of theEuler-Poincar´e equations for the diffeomorphism group.
Progr. Math. , 232:203–235,2004. http://arxiv.org/abs/nlin.CD/0312048.[HMR98] D. D. Holm, J. E. Marsden, and T. S. Ratiu. The Euler–Poincar´e equations and semidirectproducts with applications to continuum theories.
Adv. in Math. , 137:1–81, 1998.http://xxx.lanl.gov/abs/chao-dyn/9801015.[JMOB05] O. Junge, J. Marsden, and S. Ober-Blobaum. Discrete mechanics and optimal control. In
Proceedings of the 16th IFAC World Congress , 2005.[JOP + Finite Element Methods: 1970s and Beyond , pages 85–146. CIMNE,Barcelona, Spain, 2003.[LR05] B. Leimkuhler and S. Reich.
Simulating Hamiltonian Dynamics . CUP, 2005.[MM06a] R. McLachlan and S. Marsland. Discrete mechanics and optimal control for imageregistration.
ANZIAM Journal , 48, 2006.[MM06b] R. I. McLachlan and S. Marsland. The Kelvin-Helmholtz instability of momentum sheetsin the Euler equations for planar diffeomorphisms.
SIAM Journal on Applied DynamicalSystems , 5(4):726–758, 2006.[MMS06] A Mills, S Marsland, and T Shardlow.
Biomedical Image Registration , chapter Computingthe Geodesic Interpolating Spline. Number 169–177. Springer, 2006.[MTY03] M. I. Miller, A. Trouv´e, and L. Younes. Geodesic shooting in computational anatomy.Technical report, Center for Imaging Science, Johns Hopkins University, 2003.[MY01] M. I. Miller and L. Younes. Group action, diffeomorphism and matching: a generalframework.
Int. J. Comp. Vis. , 41:61–84, 2001.[TMT02] C. Twinings, S. Marsland, and C. Taylor. Measuring geodesic distances on the space ofbounded diffeomorphisms. In
British Macine Vision Conference , 2002.[TY04] A Trouv´e and L Younes. Metamorphoses through lie group action. Technical report,Center for Imaging Science, Johns Hopkins University, 2004.[VG05] Marc Vaillant and Joan Glaunes. Surface matching via currents. In