Efficient Gauss-Newton-Krylov momentum conservation constrained PDE-LDDMM using the band-limited vector field parameterization
Monica Hernandez






Efficient Jacobi PDE-EPDiff LDDMM     





Computer Sciences Department 
 Aragon Institute on Engineering Research 
 University of Zaragoza 
 mhg@unizar.es

 Efficient Gauss-Newton-Krylov momentum conservation constrained PDE-LDDMM using the band-limited vector field parameterization 
    Monica Hernandez
    Received ...; accepted ...
================================================================================================================================





The class of non-rigid registration methods proposed in the framework of PDE-constrained 
Large Deformation Diffeomorphic Metric Mapping is a particularly interesting family of physically 
meaningful diffeomorphic registration methods. 
PDE-constrained LDDMM methods are formulated as constrained variational problems,
where the different physical models are imposed using the associated partial differential equations 
as hard constraints.
Inexact Newton-Krylov optimization has shown an excellent numerical accuracy and an extraordinarily 
fast convergence rate in this framework.
However, the Galerkin representation of the non-stationary velocity fields does not provide proper geodesic paths.
In a previous work, we proposed a method for PDE-constrained LDDMM parameterized in the space of initial velocity 
fields under the EPDiff equation.
The proposed method provided geodesics in the framework of PDE-constrained LDDMM, and 
it showed performance competitive to benchmark PDE-constrained LDDMM and EPDiff-LDDMM methods.
However, the major drawback of this method was the large memory load inherent to PDE-constrained LDDMM methods and 
the increased computational time with respect to the benchmark methods.
In this work we optimize the computational complexity of the method using the band-limited vector field parameterization
closing the loop with our previous works.







§ INTRODUCTION


Deformable image registration is the process of computing spatial transformations between different
images so that corresponding points represent the same anatomical location. 
There exists a vast literature on deformable image registration methods with differences on the transformation
characterization, regularizers, image similarity metrics, optimization methods, and additional 
constraints <cit.>.
In the last two decades, diffeomorphic registration has arisen as a powerful paradigm for deformable image
registration <cit.>. 
Diffeomorphisms (i.e., smooth and invertible transformations) have become fundamental inputs in Computational
Anatomy.
There exist different big families of diffeomorphic registration methods.
Our attention in the last years has been focused to PDE-constrained diffeomorphic registration due to its relevance 
in the last decade.

PDE-constrained diffeomorphic registration augments the original variational formulation with Partial 
Differential Equations (PDEs) of interest. 
The framework seems to be very appropriate for the computation of physically meaningful transformations.
The first method was proposed by Hart et al. <cit.>. 
In that work, the problem was formulated as a PDE-constrained control problem subject to the state PDE
and the relationship with Beg et al. LDDMM <cit.> was stated.
Later on, Vialard et al. proposed a PDE-constrained method parameterized on the initial 
momentum <cit.>.
More recently, Mang et al. have proposed a PDE-constrained method that extended the gradient-descent 
optimization in Hart et al. approach to inexact Newton-Krylov optimization <cit.>.
In a previous ArXiv publication, we have bridged the gap between Vialard et al. and Mang et al. work
by proposing a novel PDE-constrained LDDMM method parameterized on the initial momentum <cit.>.
In addition we have faced the huge computational complexity of PDE-constrained LDDMM using the band-limited vector
field parameterization <cit.>.

This work closes the loop between our previous works <cit.>
by formulating the method in <cit.> in the space of band-limited vector fields.
This document is intended to be a self contained equation guide of all the related methods and provide the equations of 
the closing loop methods. The results section shows, as a proof of concept, the potential of our efficient method 
for Computational Anatomy applications.



§ RELATED METHODS




 §.§ LDDMM


Let I_0, and I_1 be the source and the target images defined on the image domain Ω⊆ℝ^d.
We denote with Diff(Ω) to the Riemannian manifold of diffeomorphisms on Ω.
V is the tangent space of the Riemannian structure at the identity diffeomorphism, id.
V is made of smooth vector fields on Ω.
The Riemannian metric is defined from the scalar product in V 

    ⟨ v, w ⟩_V = ⟨ v, w ⟩_L^2 = ∫_Ω⟨ v(x), w(x) ⟩ dΩ,

where = (Id - αΔ)^s, α >0, s ∈ℕ
is the invertible self-adjoint differential operator associated with the differential structure of Diff(Ω).
We denote with K to the inverse of L.

The LDDMM variational problem is given by the minimization of the energy functional

    E(v) =  1/2∫_0^1 ⟨ v_t, v_t ⟩_L^2 dt + 1/σ^2‖ I_0 ∘ (ϕ^v_1)^-1 - I_1‖_L^2^2.


The problem is posed in the space of time-varying smooth flows of velocity fields in V, v∈ L^2([0,1],V).
Given the smooth flow v:[0,1] → V, v_t:Ω→^d∈ V, 
the diffeomorphism ϕ_1^v is defined as the solution at time 1 to the transport equation 
d_t ϕ_t^v = v_t ∘ϕ_t^v with initial condition ϕ_0^v = id.
The transformation (ϕ^v_1)^-1 computed from the minimum of E(v) is the diffeomorphism that 
solves the LDDMM registration problem between I_0 and I_1.
The optimization of Equation <ref> was originally approached in <cit.> using gradient-descent in L^2([0,1],V), 
yielding the update equation 


    v_t^n+1 = v_t^n - ϵ (∇_v E(v))_t.




 §.§ EPDiff-LDDMM


The geodesics of Diff(Ω) under the right-invariant Riemannian metric are uniquely determined by 
the time-varying flows of velocity fields that satisfy the Euler-Poincaré equation (EPDiff) <cit.>

    ∂_t v = -ad_v^† v = - K ad_v^* L v = - K [ (Dv)^T Lv + D(Lv) v + L v ∇· v ].

with initial condition v_0 ∈ V.

LDDMM can be posed in the space of initial velocity fields
 

    E(v_0) = 1/2⟨ v_0, v_0 ⟩_L^2 + 1/σ^2‖ I_0 ∘ (ϕ^v_1)^-1 - I_1 ‖_L^2^2,


where (ϕ^v_1)^-1 is the solution at time 1 to the transport equation of the flow v_t that satisfies 
the EPDiff equation for v_0. 
The optimization of Equation <ref> was originally approached using gradient-descent in V <cit.>

    v_0^n+1 = v_0^n - ϵ∇_v_0 E(v_0).


More recently, it has been proposed in <cit.> to compute the gradient at t = 1 and to integrate backward
the reduced adjoint Jacobi field equations <cit.> 

    ∂_t U_t + ad_v_t^† U_t = 0  in Ω× [0,1) 
    ∂_t δ v_t + U_t - ad_v_tδ v_t + ad_δ v_t^† v_t = 0  in Ω× [0,1)

with initial conditions U(1) = ∇_v_1 E(v_0) and δ v(1) = 0,
to get the gradient update at t = 0, 


    v_0^n+1 = v_0^n - ϵδ v(0).




 §.§ PDE-LDDMM subject to the state equation


The PDE-constrained LDDMM variational problem is given by the minimization of 

    E(v) = 1/2∫_0^1 ⟨ v_t, v_t ⟩_L^2 dt + 1/σ^2‖ m(1) - I_1‖_L^2^2,

subject to the state equation

    ∂_t m(t) + ∇ m(t) · v_t = 0  in Ω× (0,1],

with initial condition m(0) = I_0. 
The compressible PDE-constrained problem was proposed by Hart et al. with gradient-descent 
optimization <cit.>.
Mag et al. introduced the incompressibility constraint and solved the problem using inexact 
Newton-Krylov optimization <cit.>.
 
In PDE- LDDMM, the gradient and the Hessian are computed using the method of Lagrange multipliers. 
Thus, we define the Lagrange multiplier λ:Ω× [0,1] →ℝ associated with 
the state equation, and we build the augmented Lagrangian

    E_aug(v) =  E(v) +  
    ∫_0^1 ⟨λ(t), ∂_t m(t) + ∇ m(t) · v_t ⟩_L^2 dt.


The first-order variation of the augmented Lagrangian yields the expression of the gradient 

    ∂_t m(t) + ∇ m(t) · v_t = 0  in Ω× (0,1] 
    
    -∂_t λ(t) - ∇· ( λ(t) · v_t ) = 0  in Ω× [0,1) 
    
    (∇_v E_aug(v))_t =  v_t + λ(t) ·∇ m(t)  in Ω× [0,1]


subject to the initial and final conditions 
m(0) = I_0 and λ(1) = -2/σ^2(m(1) - I_1)  in Ω. 
In the following, we will recall m as the state variable and λ as the adjoint variable.
Equations <ref> and <ref> will be recalled as the state and adjoint equations, respectively.

The second-order variation of the augmented Lagrangian yields the expression of the Hessian-vector product

    (H_v E_aug(v))_t ·δ v(t) = 
    δ v(t) + δλ(t) ·∇ m(t) - λ(t) ·∇δ m(t)  in Ω× [0,1]

where

    ∂_t δ m (t) + ∇δ m(t) · v_t + ∇ m(t) ·δ v(t) = 0  in Ω× (0,1] 
    
    -∂_t δλ(t) - ∇· ( δλ(t) · v_t ) + ∇· ( λ(t) ·δ v(t) ) = 0  in Ω× [0,1)

subject to 
δ m(0) = 0 and δλ(1) = -2/σ^2δ m(1)  in Ω.
Equation <ref> corresponds with the incremental state equation.
Equation <ref> corresponds with the incremental adjoint equation.



 §.§ PDE-LDDMM subject to the deformation state equation


This variant of PDE-constrained LDDMM is formulated from the minimization of Equation <ref>
subject to the deformation state equation

    ∂_t ϕ(t) + D ϕ(t) · v_t = 0  in Ω× (0,1],

and the incompressibility constraint

    γ∇· v_t = 0  in Ω× [0, 1].
 
The compressible PDE-constrained problem was proposed by Polzin et al. with gradient-descent 
optimization <cit.>.

The Lagrange multipliers are ρ:Ω× [0,1] →ℝ^d, associated with 
the deformation state equation, and p: Ω× [0,1] →ℝ^d, associated with the incompressibility 
constraint.
The augmented Lagrangian is given by

    E_aug(v) =  E(v) +  
    ∫_0^1 ⟨ρ(t), ∂_t ϕ(t) + D ϕ(t) · v_t ⟩_L^2 dt.


The expression of the gradient is given by the first-order variation of the augmented Lagrangian

    ∂_t ϕ(t) + D ϕ(t) · v_t = 0  in Ω× (0,1] 
    
    -∂_t ρ(t) - ∇· (ρ(t) · v_t ) = 0  in Ω× [0,1) 
    
    (∇_v E_aug(v))_t =  v_t + D ϕ(t) ·ρ(t)  in Ω× [0,1]

subject to the initial and final conditions 
ϕ(0) = id, and ρ(1) = λ(1) ·∇ m(1). 

From the second-order variation of the augmented Lagrangian, we obtain the expression of the Hessian-vector product 

    (H_v E_aug(v))_t ·δ v(t) = 
    δ v(t) + D δϕ(t) ·ρ(t) -  D ϕ(t) ·δρ(t)  in Ω× [0,1]


where


    ∂_t δϕ(t) + D δϕ(t) · v_t + D ϕ(t) ·δ v(t) = 0  in Ω× (0,1] 
    
    -∂_t δρ(t) - ∇· ( δρ(t) · v_t ) + ∇· (ρ(t) ·δ v(t)) = 0  in Ω× [0,1)


subject to δϕ(0) = 0, δρ(1) = δλ(1) ·∇ m(1) - λ(1) ·δ∇ m(1).



 §.§ Jacobi PDE-EPDiff LDDMM subject to the state equation


In <cit.> we proposed bridging the gap between momentum conservation constrained LDDMM and 
the PDE-constrained LDDMM method in <cit.> using the adjoint Jacobi field equations.
With this approach, the integration of the adjoint equation is not needed.
During the derivation of the equations we first explored the idea of transporting the gradient and the Hessian-vector
product using the adjoint Jacobi equations.
However, we found that the resulting method did not converge.
We found that transporting the vectors differently (i.e. using the adjoint Jacobi equations for the gradient, 
and the incremental adjoint Jacobi equations for the Hessian-vector products) yields the desired convergence behavior.

Thus, the PDE-constrained problem is given by the minimization of the energy functional

    E(v_0) = 1/2⟨ L v_0, v_0 ⟩_L^2 + 1/σ^2‖ m(1) - I_1 ‖_L^2^2,

subject to the EPDiff and the state equations

    ∂_t v_t + ad_v_t^† v_t = 0  in Ω× (0,1] 
    ∂_t m(t) + ∇ m(t) · v_t = 0  in Ω× (0,1],


with initial conditions v(0) = v_0 and m(0) = I_0, respectively.

Optimization is performed combining the method of Lagrange multipliers with inexact Gauss-Newton-Krylov methods
in the following way. Let
w: Ω× [0,1] →ℝ^d and
λ: Ω× [0,1] →ℝ 
be the Lagrange multipliers associated with the EPDiff and the state equations.
We build the augmented Lagrangian

    E_aug(v_0) =  E(v_0) + ∫_0^1 ⟨ w(t), ∂_t v(t) + ad_v_t^† v_t ⟩_L^2 dt 
    
    + ∫_0^1 ⟨λ(t), ∂_t m(t) + ∇ m(t) · v_t ⟩_L^2 dt.


Similarly to <cit.>, the gradient is computed at t = 1, ∇_v_1 E(v_0) = λ(1) ·∇ m(1)



and integrated backward using the reduced adjoint Jacobi field equations (Equation <ref>)
to obtain ∇_v_0 E(v_0).









The second-order variations of the augmented Lagrangian on w and λ yield the incremental EPDiff and 
incremental state equations, needed for the computation of the Hessian-vector product. Thus,

    ∂_t δ v_t + ad_δv_t^† v_t + ad_v_t^†δ v_t= 0  in Ω× (0,1] 
    ∂_t δ m(t) + ∇δ m(t) · v_t + ∇ m(t) ·δ v_t = 0  in Ω× (0,1]

with initial conditions δ v(0) = 0 and δ m(0) = 0.

The Hessian-vector product H_v_0 E(v_0) ·δv_0 is computed from the Hessian-vector product at t = 1, which is integrated backward
using the reduced incremental adjoint Jacobi field equations

    ∂_t δ U + ad_δ v^† U + ad_v^†δ U= 0  in Ω× [0,1) 
    ∂_t δ w + δ U - ad_δ v w - ad_v δ w + ad_δ w^† v + ad_w^†δ v = 0  in Ω× [0,1)


with initial conditions δ U(1) = K( δλ(1) ·∇ m(1) - λ(1) ·∇δ m(1) ), 
and δ w(1) = 0. 





 §.§ Jacobi PDE-EPDiff LDDMM subject the deformation state equation


In <cit.> we explored the behavior of different variants of the PDE-constrained LDDMM problem with the band-limited vector field parameterization.
The best performing method was the Newton-Krylov extension of the method proposed in <cit.>.
In this work we provide the equations of the Jacobi PDE-EPDiff LDDMM version of the method.

The PDE-constrained problem is given by the minimization of the energy functional

    E(v_0) = 1/2⟨ L v_0, v_0 ⟩_L^2 + 1/σ^2‖ m(1) - I_1 ‖_L^2^2,

subject to the EPDiff and the deformation state equations

    ∂_t v_t + ad_v_t^† v_t = 0  in Ω× (0,1] 
    ∂_t ϕ(t) + D ϕ(t) · v_t = 0  in Ω× (0,1],


with initial conditions v(0) = v_0 and ϕ(0) = id, respectively.
The state variable m is computed from m(t) = I_0 ∘ϕ(t).

The augmented Lagrangian is given by

    E_aug(v_0) =  E(v_0) + ∫_0^1 ⟨ w(t), ∂_t v(t) + ad_v_t^† v_t ⟩_L^2 dt 
    
    + ∫_0^1 ⟨ρ(t), ∂_t ϕ(t) + D ϕ(t) · v_t ⟩_L^2 dt.


The gradient is computed at t = 1, ∇_v_1 E(v_0) = D ϕ(1) ·ρ(1)
and integrated backward using the reduced adjoint Jacobi field equations (Equation <ref>)
to obtain ∇_v_0 E(v_0).
Thus, 


    ∂_t ϕ(t) + D ϕ(t) · v_t = 0  in Ω× (0,1], 
    ∂_t U_t + ad_v_t^† U_t = 0  in Ω× [0,1) 
    ∂_t δ v_t + U_t - ad_v_tδ v_t + ad_δ v_t^† v_t = 0  in Ω× [0,1)


where U(1) = K( D ϕ(1) ·ρ(1) ).

The Hessian-vector product H_v_0 E( v_0 ) ·δv_0 is computed from the Hessian-vector product at t = 1, which is integrated backward
using the reduced incremental adjoint Jacobi field equations


    ∂_t δϕ(t) + D δϕ(t) · v_t + D ϕ(t) ·δ v_t = 0  in Ω× (0,1], 
    ∂_t δ U_t + ad_δ v_t^† U_t + ad_v_t^†δ U_t = 0  in Ω× [0,1) 
    ∂_t δ w_t + δ U_t - ad_δ v_t w_t - ad_v_tδ w_t + ad_δ w_t^† v_t + ad_w_t^†δ v_t = 0  in Ω× [0,1)


with initial conditions δ U(1) = K( D δϕ(1) ·ρ(1) - D ϕ(1) ·δρ(1) ), and δ w(1) = 0.



 §.§ Gauss-Newton-Krylov optimization


By construction, the Hessian is positive definite in the proximity of a local minimum.
However, it can be indefinite or singular far away from the solution.
In this case, the search directions obtained with PCG are not guaranteed to be descent directions. 
In order to overcome this problem, one can use a Gauss-Newton approximation dropping expressions of H_v E(v_0) ·δ  v_0
to guarantee that the matrix is definite positive.

The minimization using a second-order inexact Gauss-Newton-Krylov method yields to the update equation

    v_0^n+1 = v_0^n - ϵδ v_0^n,

where δ v_0^n is computed from PCG on the system 

    H_v_0 E( v_0^n) ·δv_0^n = - ∇_v_0 E(v_0^n).

In this work, we consider CG with the gradient and the Hessian computed on V instead of L^2.



§ METHODS PARAMETERIZED IN THE SPACE OF BAND-LIMITED VECTOR FIELDS




 §.§ Background on the space of band-limited vector fields


Let Ω be the discrete Fourier domain truncated with frequency bounds K_1, …, K_d.
We denote with V the space of discretized band-limited vector fields on Ω with these frequency bounds. 
The elements in V are represented in the Fourier domain as ṽ: Ω→ℂ^d,
ṽ(k_1, …, k_d), and in the spatial domain as ι(ṽ):Ω→ℝ^d,

    ι(ṽ)(x_1, …, x_d) = ∑_k_1 = 0^K_1…∑_k_d = 0^K_dṽ(k_1, …, k_d) e^2 π i k_1 x_1… e^2 π i k_d x_d.


The application ι:V→ V denotes the natural inclusion mapping of V in V.
The aplication π: V →V denotes the projection of V onto V.

The space of band-limited vector fields V has a finite-dimensional Lie algebra structure using the truncated 
convolution in the definition of the Lie bracket <cit.>.
We denote with Diff(Ω) to the finite-dimensional Riemannian manifold of diffeomorphisms on Ω 
with corresponding Lie algebra V.
The Riemannian metric in Diff(Ω) is defined from the scalar product

    ⟨ṽ, w̃⟩_Ṽ = ⟨ṽ, w̃⟩_l^2,

where  is the projection of operator  in the truncated Fourier domain. 
Similarly, we will denote with K̃, ∇, and ∇· to 
the projection of operators K, ∇, and ∇· in the truncated Fourier domain.
In addition, we will denote with ⋆ to the truncated convolution.

The EPDiff-equation in the space of band-limited vector fields is given by

    ∂_t ṽ_t = ad^†_ṽ_tṽ_t = -K̃ [ (D̃ṽ)^T ⋆L̃ṽ + D̃(L̃ṽ) ⋆ṽ + L̃ṽ∇⋆ṽ ].

The adjoint operator is given by

    ad_ṽw̃ = D̃ṽ⋆w̃ - D̃w̃⋆ṽ.




 §.§ BL Jacobi PDE-EPDiff LDDMM subject to the state equation


The variational problem is given by the minimization of 

    E(ṽ_0) = 1/2⟨L̃ṽ_0, ṽ_0 ⟩_l^2 + 1/σ^2‖ m(1) - I_1 ‖_L^2^2

subject to 

    ∂_t ṽ_t + ad^†_ṽ_tṽ_t = 0  in Ω× (0,1] 
    ∂_t m(t) + ∇ m(t) ·ι(ṽ_t) = 0  in Ω× (0,1],

with initial conditions ṽ(0) = ṽ_0 and m(0) = I_0. 

The expression of the gradient is computed from the reduced adjoint Jacobi field equations in the space of band-limited vectors yielding

    ∂_t Ũ_t + ad^†_ṽ_tŨ_t = 0  in Ω× [0,1) 
    ∂_t δṽ_t + Ũ - ad^†_ṽ_tδṽ_t + ad^†_δṽ_tṽ_t = 0  in Ω× [0,1)

where  

    Ũ(1) = K̃( π( λ(1) ·∇ m(1) ) ).


The expression of the Hessian-vector product is computed from the reduced incremental adjoint Jacobi equations in the space of band-limited vector fields

    ∂_t δŨ_t + ad^†_δ̃ ̃ṽ_tŨ_t + ad^†_ṽ_tδŨ_t = 0  in Ω× [0,1) 
    ∂_t δw̃_t + δŨ_t - ad_δṽ_tw̃_t - ad_ṽ_tδw̃_t + ad_δw̃_t^†ṽ_t + ad_w̃_t^†δṽ_t = 0  in Ω× [0,1)

where  

    δŨ(1) = π( K̃(δλ(1) ·∇ m(1) - λ(1) ·∇δ m(1)) ).




 §.§ BL Jacobi PDE-EPDiff LDDMM subject to the deformation state equation


The variational problem is given by the minimization of Equation <ref> subject to 

    ∂_t ṽ_t + ad^†_ṽ_tṽ_t = 0  in Ω× [0,1) 
    ∂_t ϕ̃(t) + Dϕ̃(t) ⋆ṽ_t = 0  in Ω× [0,1).

The expression of the gradient is given by

    ∂_t ṽ_t + ad^†_ṽ_tṽ_t = 0  in Ω× [0,1) 
    ∂_t ϕ̃(t) + Dϕ̃(t) ⋆ṽ_t = 0  in Ω× (0,1]  
    ∂_t Ũ_t + ad^†_ṽ_tŨ_t = 0  in Ω× (0,1] 
    ∂_t δṽ + Ũ - ad_ṽδṽ + ad^†_δṽṽ = 0  in Ω× [0,1)


with initial conditions Ũ(1) = K̃(D̃ϕ̃(1) ⋆ρ̃(1)).

On the other hand, the expression of the Hessian-vector product is given by

    ∂_t δṽ_t + ad^†_δ̃ ̃ṽ_tṽ_t + ad^†_ṽ_tδ̃ ̃ṽ_t = 0  in Ω× [0,1) 
    ∂_t δϕ̃(t) + Dδϕ̃(t) ⋆ṽ_t + Dϕ̃(t) ⋆δṽ(t) = 0  
    ∂_t δŨ_t + ad^†_δ̃ ̃ṽ_tŨ_t + ad^†_ṽ_tδŨ_t = 0  in Ω× [0,1) 
    ∂_t δw̃_t + δŨ_t - ad_δṽ_tw̃_t - ad_ṽ_tδw̃_t + ad_δw̃_t^†ṽ_t + ad_w̃_t^†δṽ_t = 0  in Ω× [0,1)

where  

    δŨ(1) = K̃(D̃δϕ̃(1) ⋆ρ̃(1) - D̃ϕ̃(1) ⋆δρ̃(1)).




§ RESULTS


The experiments have been conducted on the Non-rigid Image Registration Evaluation Project database (NIREP) with volumes of size 180 × 210 × 180. 
Figure <ref> shows the source and target images and the differences before registration.



Figure <ref> shows the MSE_rel and ‖ g‖_∞,rel convergence curves obtained during the optimization for band sizes
of 16, 32, 40, 48, 56 and 64. The figure shows that both methods converge to similar MSE_rel values. 
However, Jacobi PDE-LDDMM subject to the deformation state equation shows smaller ‖ g‖_∞,rel, which indicates a better convergence behavior.
Table <ref> shows the numeric values after 10 iterations. 
For BL sizes of 32 both methods achieve acceptable MSE_rel values.
Jacobi PDE-LDDMM subject to the deformation state equation slightly outperformed the method subject to the state equation.





Figure <ref> shows the deformed images and the differences after registration.
In the figure, it can be appreciated the accuracy achieved by the proposed methods.
The difference between the deformed images is hardly perceptible.



Finally, Table <ref> shows the VRAM memory load and the computation time 
exhibit by the proposed methods. It should be noticed that Jacobi PDE-LDDMM subject to the state equation
in the spatial domain did not fit the memory of our available graphics card (11 GBs) and the computation
time for a downsampled example of our data was 1825.50 seconds. Therefore, the band limited vector field
parameterization definitively shows up a considerable computational saving.





§ ACKNOWLEDGEMENTS

This work was partially supported by Spanish research grant TIN2016-80347-R.

splncs


