Estimating Diffeomorphic Mappings between Templates and Noisy Data: Variance Bounds on the Estimated Canonical Volume Form
Daniel J. Tward Partha Mitra Michael I. Miller


AIP/123-QED







 [Also at ]Physics Department, XYZ University.

 Second.Author@institution.edu

 Authors' institution and/or address

 This line break forced with \\

MUSO Collaboration


 http://www.Second.institution.edu/ Charlie.Author

 Second institution and/or address

 This line break forced


 Third institution, the second for Charlie Author



 Authors' institution and/or address

 This line break forced with \\


CLEO Collaboration



dtward@cis.jhu.edu
Center for Imaging Science, Johns Hopkins University, Baltimore, MD, 21218

Cold Spring Harbor Laboratory, Cold Spring Harbor, NY, 11724

Department of Biomedical Engineering, Johns Hopkins University, Baltimore, MD, 21218



             



An article usually includes an abstract, a concise summary of the work
covered at length in the main body of the article. 


  Usage
Secondary publications and information retrieval purposes.

  PACS numbers
May be entered using the  command.

  Structure
You may use the  environment to structure your abstract;
use the optional argument of the  command to give the category of each item. 



Valid PACS appear here
                             

                              

Estimating Diffeomorphic Mappings between Templates and Noisy Data: Variance Bounds on the Estimated Canonical Volume Form


    Michael I. Miller
    December 30, 2023
=============================================================================================================================










Anatomy is undergoing a renaissance driven by availability of large digital data sets generated by light microscopy. A central computational task is to map individual data volumes to standardized templates. This is accomplished by regularized estimation of a diffeomorphic transformation between the coordinate systems of the individual data and the template, building the transformation incrementally by integrating a smooth flow field. The canonical volume form of this transformation is used to quantify local growth, atrophy, or cell density. While multiple implementations exist for this estimation, less attention has been paid to the variance of the estimated diffeomorphism for noisy data. Notably, there is an infinite dimensional un-observable space defined by those diffeomorphisms which leave the template invariant. These form the stabilizer subgroup of the diffeomorphic group acting on the template. The corresponding flat directions in the energy landscape are expected to lead to increased estimation variance.  Here we show that a least-action principle used to generate geodesics in the space of diffeomorphisms connecting the subject brain to the template removes the stabilizer. This provides reduced-variance estimates of the volume form. Using simulations we demonstrate that the asymmetric large deformation diffeomorphic mapping methods (LDDMM), which explicitly incorporate the asymmetry between idealized template images and noisy empirical images, provide lower variance estimators than their symmetrized counterparts (cf. ANTs). We derive Cramer-Rao bounds for the variances in the limit of small deformations. Analytical results are shown for the Jacobian in terms of perturbations of the vector fields and divergence of the vector field.













Estimating Diffeomorphic Mappings between Templates and Noisy Data: Variance Bounds on the Estimated Canonical Volume Form


    Michael I. Miller
    December 30, 2023
=============================================================================================================================



















§ INTRODUCTION

Computational Anatomy (CA) is a growing discipline. Starting with initial work <cit.> directed towards the
the study of transformations between anatomical coordinate systems suitable for volumetric images of the sub-compartments of the human brain acquired largely using MRI, contemporary applications have been extended to much larger data volumes acquired using light microscopy.
Infinite dimensional diffeomorphisms constitute the central transformation group for studying shape and form.
<cit.>
The diffeomorphism model underlying this analysis assumes that the space of measured MRI and optical imagery can be generated from exemplars or templates via diffeomorphic changes of coordinates.

At the mesoscopic scale the variation of the diffeomorphic change in coordinates from one brain to another can represent transverse individual variation, including pathological conditions, or longitudinal developmental variation.



Of particular interest in mesoscale neuroanatomy are quantities such as the spatial densities of cellular somata or neuronal processes.
Mapping estimates of these quantities to a template or reference space requires estimation of the change in the local scale as captured by the determinant of the metric tensor (the canonical volume form determined by the Jacobian determinant of the mapping). These applications point to the importance of uncertainty estimation for the diffeomorphic transformations involved. This is the subject area of the current paper.

Dense mapping between coordinate systems began with
the low-dimensional matrix Lie groups forming the basis of Kendall's shape theory.<cit.>
Their infinite dimensional analogue, the diffeomorphisms between coordinate systems, were first introduced by Christensen et al. <cit.>
and have occupied a central role in CA.
<cit.> 
The small deformation methods have been associated with Bookstein's landmark matching <cit.> and subsequently for image matching.<cit.> 
Large deformations were studied as simply topology preserving transformations without a metric structure.<cit.>
The symmetric approaches
for large deformations were variants of these methods.<cit.>
The large deformation diffeomorphic metric mapping algorithms (LDDMM) emerged corresponding to Lagrangian and Hamilton's principles applied to the flow fields incrementally generating the diffeomorphic transformations involved
<cit.>
and provide a metric between images and diffeomorphisms. 
The survey article by <cit.> places these works in the greater context of deformable registration.



Uncertainty of nonrigid registration algorithms has been investigated for several applications such as spatially adaptive smoothing of population level data or uncertainty visualization for surgical planning.<cit.> These approaches to uncertainty have typically used resampling techniques and linear elastic or spline based deformation models as opposed to diffeomorphisms.  Theoretical bounds on uncertainty have been investigated for the case of translation-only registration.<cit.>

In this paper we study problems associated with the estimation of the (local) canonical volume form (ie √(det(g)) where g is the metric tensor of the diffeomorphism) via methods based on large deformation diffeomorphic metric mapping (LDDMM) <cit.>.
As an example of an application where knowledge of this local scale change is important, consider mapping estimated cell densities from an individual brain to a template. Investigators may wish to distinguish differential changes in cell density of different types, from an overall change in scale given by the diffeomorphism relating the individual brain to the template. This requires an estimation of the volume form.  In human neuroanatomy, this local change in scale is commonly used to quantify patterns of tissue atrophy in neurodegenerative disease  <cit.>.

In the notation followed in the paper, the coordinate transformation ϕ has  Jacobian [D ϕ]_ij  = ( ∂ϕ_i /∂ x_j) and first fundamental form (metric tensor) given by 

    g(ϕ) =  (D ϕ )^T (D ϕ )  .

The canonical volume form is given by the square-root of the determinant √( |  g |).

We show that the geodesic equations associated with the LDDMM method provide a crucial reduction, removing the nuisance dimensions of the  underlying symmetries associated with the stabilizer of the template, i.e. the subgroup of the diffeomorphic transformations that leave the deforming template unchanged.
 This template-centered reduction of the stabilizer gives rise to the asymmetric mapping properties central to LDDMM, and provides robustness when the imaging targets suffer from an incomplete and or noisy measurement process. We examine the mean square error of estimation of the
canonical volume form.
We demonstrate that the asymmetry of LDDMM coupled to geodesic reduction of indeterminate dimensions of the flow leads to favorable performance of this estimate in the presence of noise and variability when compared to symmetric methods for image matching originally proposed by
<cit.>, and <cit.>.  
We also calculate the Cramer-Rao bound for the variance of the volume form in the case of small deformations.

Explicitly studying the behavior of the canonical volume form under uncertainty, as contrasted to registration accuracy which has been addressed by the community in some situations, is essential for drawing meaningful conclusions about cell or process density and brain morphometry <cit.>.  The current work provides a theoretical basis for the observations presented and for related observations by other authors. It clarifies important differences between symmetric and asymmetric methods as they relate to the uncertainty of the volume form estimates.

The paper is organized as follows:


  * We first review the theoretical methods including two optimal control problems that are stated for retrieving the diffeomorphism and the fundamental forms describing changes of coordinate systems between images. 

  * The necessary condition for the solution to the variational problem is stated in terms of the Euler-Lagrange equation on the conjugate momentum. It is shown that this equation incorporates a Hamiltonian reduction of the infinite-dimensional symmetry group corresponding to the stabilizer of the template. 

  * We then derive a new analytical result for the Cramer-Rao bound on the variance of the fundamental form in the presence of measurement noise. 

  * Following this we show results from large deformation simulations showing the decrease in variance for the asymmetric LDDMM which privileges the template coordinate system as ground truth and removes the stabilizer with respect to the template. 

  * These results are compared to symmetric methods derived by Christensen and Johnson<cit.> and Avants and Gee et al.<cit.> 
It is shown that the Hamiltonian-reduced LDDMM and symmetric methods behave similarly at low-noise but LDDMM outperforms with increasing noise. 


The significant new results in this paper thus address uncertainty estimation of the canonical volume form for estimated diffeomorphic transformations for both large deformation LDDMM and symmetric methods as well as small deformation methods.





§ METHODS

Ethics approval is not required for this work.


 §.§ Theoretical Methods



  §.§.§ Geodesic Flows of Diffeomorphisms for Dense Transformation of Coordinate Systems














The diffeomorphism model in Computational Anatomy posits that the diffeomorphism group acts on templates I_temp(x), x ∈ X ⊂ℝ^3 via group action to generate the space of observed anatomical data from individual subjects ℐ,

    ℐ = { I = I_temp∘ϕ_t^-1 , ϕ_t ∈} ,

with diffeomorphisms ϕ_t:X  → X, t ∈ [0,1], ϕ_t ∈ generated via flows:

    ϕ̇_t = v_t ∘ϕ_t  ,  ϕ_0 =  id .

where id is the identity transformation.  This is also termed the random orbit model: in group theoretic terms, ℐ is the orbit of I_temp under .

The Eulerian vector fields v_t: X →ℝ^3
are constrained to be spatially smooth, supporting at least 1-continuous derivative in space, ensuring that the flows are well defined and with smooth inverse.<cit.>
They are modeled to be smooth with a finite
norm (V,·_V )  defined by a differential operator A: V → V^*,
v_V^2 = ∫_X Av · v dx.
It is conventional to use powers of the Laplacian for the differential operator A with a sufficient number of generalized derivatives such that the flow fields are guaranteed to have at least 1 continuous derivative in space.<cit.>
The kernel K of the associated reproducing kernel Hilbert space is at least 1-time continuously differentiable in the spatial variables, and is given by the Green's kernel of A.






Mapping individual brains to reference atlases such as the Allen mouse atlas <cit.> at the micron scale or the Mori human atlas <cit.> at millimeter scale is performed via bijective coordinate transformation between anatomical coordinate systems ϕ.

The diffeomorphic change in coordinates is not directly observable and must be inferred from observed brain image data subjected to measurement noise and technical variations. The coordinate system transformations may be estimated by solving the diffeomorphic matching problems as a solution of an optimal control problem. Different optimal control problems are obtained for large and small deformations. 

To set up the optimal control problem, we define the transformation ϕ_t as a time dependent state, t →ϕ_t ∈. The velocity field
t → v_t is taken to be the control variable. The state satisfies the dynamical equations ϕ̇_t = v_t ∘ϕ_t, with initial condition ϕ_0=id. 
The goal of the optimal control problem is to drive the state from an initial condition of identity (corresponding to the template), to a state that matches template coordinates to target coordinates. This can be enforced by a cost function for the final state,

    U(ϕ_1) = 1/2σ^2∫_X | J - I ∘ϕ_1^-1 |^2 dx  .

where J is an observed target image. We assume that any differences in position or orientation have already been accounted for by an appropriate similarity transform before deforming the coordinate grid using the diffeomorphic procedure.  A scalar parameter σ controls the importance of this data attachment term relative to the regularization term defined below.

Since there are an infinite number of possible flows we use the principle of least-action to minimize a running cost given by the integrated kinetic energy of the vector field  v = ϕ̇∘ϕ^-1 of the flow (LDDMM <cit.>).
This regularization cost term can be interpreted as a kinetic energy Lagrangian:

    L(ϕ_t,ϕ̇_t)   =1/2∫_X A (ϕ̇_t ∘ϕ_t^-1)·ϕ̇_t ∘ϕ_t^-1 dx 
       =1/2∫_X Av_t · v_t dx  .


The cost function for the optimal control problem is obtained by combining the kinetic energy term (the regularization term) with an end-point cost. Thus, the optimal control problem involves Lagrangian mechanics of the infinite dimensional state defined by the coordinate transformation ϕ_t. We study two versions of this problem, corresponding to large and small deformations (control problems 1 and 2).  
















 




 


















[Large Deformation]



    [                  ϕ̇_t= v_t ∘ϕ_t  , ϕ_0 = id; min_ϕ C(ϕ) ≐∫_0^1 L(ϕ_t,ϕ̇_t) dt   + U(ϕ_1) ]


For small deformations
<cit.>
there is no time index, and v represents a displacement field rather than a flow field, ϕ = id +v. 


[Small Deformation]


    [                                ϕ=id+v; min_v C(ϕ) ≐1/2∫_X A v · v dx  + U(ϕ) ]





The solution to the large deformation Control Problem 1 has classical conjugate momentum  p_t = Av_t ∘ϕ_t |D ϕ_t |  <cit.> which satisfies geodesic equations


    p_t = α_t ∇ I ,    α_t = 1/σ^2 (I-J∘ϕ_1) |D ϕ_1| (D ϕ_t^T)^-1 .


For small deformations, Control Problem 2,
p=Av satisfying

    p= α∇ I  ,  α=
    1/σ^2 (I-J∘ϕ )|D ϕ |(D ϕ^T)^-1 .


The conjugate momentum satisfying LDDMM 
equation (<ref>)
is perpendicular to the level lines of the template following the image gradient.
This normal condition is satisfied over the entire path of deformation.  This is observed most easily by considering the Eulerian momentum, related to the conjugate momentum by a change of coordinates, Av_t = 1/σ^2|Dϕ_1t^-1|(J(ϕ_1t^-1) - I(ϕ_t^-1))∇ [I(ϕ^-1_t)], where ϕ_1t = ϕ_t(ϕ_1^-1).  Notice that Av_t is parallel to ∇ [I(φ_t^-1)] at every point.
See Appendix <ref> for a proof.

The fact that the conjugate momentum is in the range space at every point of the gradient of the template is a Hamiltonian reduction that removes the symmetries of the template given by the stabilizer subgroup. The non-identifiable motions that are tangent to level lines are thus removed.
This kills off the nuisance dimensions, essentially suppressing components in the null-space of level lines of the template as it flows.<cit.>

For small deformation matching the normal condition of the stabilizer corresponds to a single vector condition that Av is in the span of ∇ I the gradient of the image template.
Small deformations allow us to both calculate the Cramer-Rao bound (see below) as well as perform
a direct perturbation argument. 


Via a perturbation J → J +δ  J we can calculate
the accuracy of the divergence of v and how it determines the Jacobian determinant and the canonical volume measure related to the gradient of the template.


For deformations close to the identity, a perturbation J ↦ J + δ J results in a perturbation to the canonical volume form:

    δ v    = [A + 1/σ^2 (∇ I)  (∇ I)^T]^-1∇ I δ J
    
    |Dϕ|    ↦ |Dϕ| + |Dϕ|div[δ v]

where ∇ I is a column vector, the outer-product giving a 3× 3 matrix.  

See Appendix <ref> and <ref> for small deformation proofs.  
Equation
(<ref>) clearly shows the crucial matrix involving the gradient of the image and the prior represented via the differential operator A for determining the 
perturbation on the canonical volume form through the divergence.
The Cramer-Rao bound reflects this theme.



  §.§.§ Cramer-Rao Bound (CRB) for Small Deformations


We now examine the variational estimator and the variance bound for small deformations in the finite dimensional setting of n-dimensional vector fields v^n(x) = ∑_i=1^n θ^i ψ^i(x). Here ψ^i is some suitable family of expansion functions. For computing the CRB, 
we take the observed data to be a conditionally Gaussian random field, conditioned on the mean I ∘ϕ^n-1 with additive noise

    J= I ∘ (ϕ^n)^-1 + noise   ,   ϕ^n = id + v^n  .

The noise is taken to be zero-mean with non-white inverse covariance Q. 
The log-likelihood on the n-dimensional cylinder Θ^n = (θ^1,…,θ^n) is

    ℓ_n (J;Θ_n) 
    
       =
    -
    1/2∫_X∫_X (J(x)-I ∘ϕ^n -1 (x)) 
            Q(x,y)(J(y)-I ∘ϕ^n -1 (y)) dx dy   .

For
white-noise the inverse-covariance is Q(x,y) = 1/σ^2 (x)δ(x-y) the n-dimensional log-likelihood (<ref>), for σ^2(x) a variance at location x.
Adding the finite-dimensional prior term - 1/2∫_X Av^n(x) · v^n(x) dx
gives a proper maximum a-posterior estimator (MAP) on n-dimensions:

    max_θ^1, … , θ^n logπ^n (Θ^n  | J ) 
       = - 1/2∫_X Av^n · v^n dx + ℓ^n (J; θ^1,…, θ^n)
     .

The Fisher-information is the n × n matrix, i,j=1,…,n:

    ℐ_ij(Θ_n)    = E_J|Θ^n{∂ℓ^n(J;Θ^n)/∂θ^i∂ℓ^n (J;Θ^n)/∂θ^j} .

Taking expectation over Θ^n gives the Bayesian version of the Fisher-information: 

    ℐ^B = ∫_R^nℐ(J;Θ^n) π(d Θ_n)   .
   

Defining for all (x,y),
Q^ϕ (x, y ) = Q( ϕ(x), ϕ(y) |D ϕ (x) | | D ϕ (y) |,

    ℐ_ij(Θ_n)
    = ∫_X ∫_X   ψ^iT(x)[Dϕ(x)]^-T∇ I(x) 
            Q^ϕ (x,y)  ∇  I^T(y)[Dϕ(y)]^-1ψ^j(y)dxdy  
     .

Neglecting quadratic terms using small deformations, noting that the linear terms have an expected value of 0, and taking white noise variance σ^2  gives

    ℐ^B_ij   = ∫_X ψ_i^T(x) (  A  +  1/σ^2(x)∇ I(x) ∇  I^T(x)  ) ψ_j(x) dx 
            + E[o(v^2)]  ,

giving the lower bound on the estimator:

    Cov[Θ̂_n] ≥  [ℐ^B]^-1


Noteworthy is
the CR bound (<ref>) is at the core the same form as direct perturbation
(<ref>).
We are particularly interested in the variance of div(v̂), because it directly relates to variance of the Jacobian.  This can be written as a linear functional of v, using a test function w which is nonzero in a small neighborhood of the point y (i.e. the sensitivity of an image voxel), namely ∫_X w(x)divv(x)dx.
In this case the information inequality, which is concerned with estimating functions of the parameter v, becomes

    Var[div[v̂](y)]    ≥∑_i,j=1^n ( ∫_X w(x) divψ_i(x) dx) 
            [ℐ^B]^-1_ij( ∫_X w(x) divψ_j(x)dx)








































  §.§.§ The Stabilizer: The normal condition for the nonlinear group action

Unfortunately there is an infinite dimensional
nuisance parameter, termed the stabilizer, which is not uniquely determined in estimating the coordinates of the bijections between coordinate systems.
The geodesic motion kills the nuisance parameter, suppressing components in the instantaneous null-space of level lines of the template as it flows <cit.>;
this is Theorem 4 of <cit.>.
Notice there is no momentum tangent to the level lines of the flow of the template image, p_t = α_t ∇ I in 
(<ref>).
This is not true for other diffeomorphic mapping methods.
It results from the metric property of the geodesic equation
with it's associated conservation law. 


To understand flows which are normal to the level lines and in the span of the gradient of the template,
examine the stabilizer. 


Define the subgroup
𝒮⊂ as the stabilizer of template I if for all ϕ∈𝒮,

    I ∘ϕ = I  .




Figure <ref> shows examples of mappings from the stabilizer of an image of a human hippocampus, to a close numerical approximation. The left column shows the identity mapping on the grid; the right two columns show two mappings from the stabilizer.


Vector fields that are normal to the level lines in the span ∇ I do not generate flow through the stabilizer group.
The stabilizer is generated from flows ϕ̇_t = w ∘ϕ_t , ϕ_0=id for which the vector fields are

tangent to level lines of the template:

    V_I = { w ∈ V : ∇ I · w = 0 } .

That being normal to the level lines is a necessary condition to null out the stabilizer, examine 

ϕ_ϵ = id + ϵ w, then

    I ∘ϕ_ϵ = I  + ϵ∇ I ·  w + o(ϵ)  .

For I ∘ϕ_ϵ = I it must be the case that ∇ I · w = 0, implying w must be normal
to the level lines.
This is the group action version of the pseudo-inverse condition for inverting a matrix with non-zero null-space.





 §.§ Numerical Experiments



  §.§.§ Simulating Variance Bounds on Large Deformations

The geodesic properties of LDDMM imply lower variance estimates of the first fundamental form and the fundamental volume measure.
To illustrate this we create simulated images by generating large deformation diffeomorphisms under the random orbit model,  deforming a template image under these known transformations, and applying Gaussian noise.  Our template corresponds to the binary segmentations of the anatomical section in Fig. <ref> (left column),  a sagittal section of an ex vivo image of the human medial temporal lobe, where hippocampal subfields and surrounding areas are visible.  


We generate random Gaussian vector fields v_0(x) = v^n(x), x ∈ X as the initial condition of geodesic solutions for deformation.
Since geodesics satisfy a conservation law 

    d/dt D ϕ_t^T p_t = 0

proved in Appendix <ref> equation <ref>, we can synthesize the initial condition p_0 = Av^n and generate a random spray of deformations ϕ_1.
The Gaussian vector field is
generated from an n-dimensional expansion ψ_i:  X →ℝ^3,i=1,2,…,

    v^n(x)    = ∑_i = 1^n θ^i ψ^i(x)  ,

coefficients θ^i ∈ℝ, i = 1,…,n distributed as multivariate Gaussian, zero-mean θ̅^i = 0 and covariance E[θ^i θ^j] = Σ^ij.

We choose the covariance of vector fields v^n to be well defined for large n,  sampled at pairs of points x,y ∈ℝ^3, 

    K^n(x,y)    ≐ E[v^n(x)v^nT(y)] = ∑_i,jψ^i(x) E[θ^i θ^j] ψ^jT(x)  .

For large n this tends to covariance specified by the Green's kernel of the operator A originally stated by Beg
<cit.>,
 the inverse of A = (1 -a^2Δ)^4 for Δ the Laplacian at spatial scale of a = 0.25  mm, or 2 pixels).  To achieve this, we choose an expansion ψ^i that corresponds to a superposition of Green's kernels located on each voxel where the image gradient is nonzero (3432 expansion functions).
  
For the statistical characterization, 100 realizations of random deformation are generated running equation (<ref>) from v^n as initial condition.  100 randomly generated images were constructed from these deformations applied to our template (Fig. <ref> left), embedded in additive Gaussian noise with standard deviation from 0 to 0.5, spatial correlation of 0 (white noise) or 1.5 pixels (obtained by convolving with a Gaussian of standard deviation 1.5 pixels).  Correlated noise is common in medical imaging systems such as radiography <cit.>, CT <cit.>, or MRI <cit.>.  Similarly, 100 randomly generated images were constructed using the identity transformation instead of a random transformation.


Applying diffeomorphic mapping for 100 images produced using the identity transformation (θ_i=0 ∀ i), and 100 images produced using these known transformations then gives variance estimates of the diffeomorphic mapping methods.  For contrast we compare LDDMM to the symmetrized version discussed in <cit.> ANTs which does not exploit the reduced representation as symmetry supports momentum on both the template and target:

    max_v: ϕ^v = ∫_0^1 v_t ∘ϕ_t dt + id
    -1/2∫_0^1 v_t^2_V dx - 1/2σ^2I∘ϕ^-1 - J^2_L_2
                -1/2∫_0^1 v_t^2_V dx  - 1/2σ^2I - J∘ϕ^2_L_2 .

Gradient descent is used to determine the optimal v, giving estimates of  ϕ̂ and the canonical form |D ϕ̂|. We use exactly the same parameters for both algorithms.

Parameters were fixed or varied experimentally as summarized in Table <ref> in the Appendix <ref>.  In our computational implementation we use two constants in our objective function to be minimized:  1/σ^2_V multiplying the regularization (kinetic energy) term, and 1/σ^2_I multiplying the data attachment term.  Since we are estimating the optimizer, but not the optimum, this is equivalent to choosing a single 1/σ^2 = σ^2_V/σ^2_I.




  §.§.§ Calculating the Cramer Rao Bound

Finite dimensional cylinders were used to illustrate the Cramer Rao lower bound on the divergence of v̂_n = ∑_i=1^n θ_i ψ_i from  (<ref>).  Expansion functions are Green's kernels of A located on image boundaries, where the image is downsampled by a factor of (from left to right) 8 (76 expansion functions), 4 (318 expansion functions), 
or 1 (not downsampled, 3432 expansion functions).  Explicitly inverting the ill-conditioned matrix ℐ^B was avoided by solving the linear system implied by (<ref>), which was performed using Matlab's  for positive definite symmetric matrices.



§ RESULTS



 §.§ Variance bounds

Figure <ref> shows sections (column 1) through the high-field medial temporal lobe MRI phantom used for numerical experiments, with the segmentation into the entorhinal cortex and it substructures subiculum and Cornu Ammonis (CA) compartments.
Figure <ref> (columns 2,3) show results for simulating known transformations of the template with additive white noise.
Column 2 shows the results of solving LDDMM equation (<ref>)
to estimate each random deformation ϕ. The top row shows identity transformation of coordinates; the bottom row shows random deformations ϕ.
LDDMM in homogeneous regions has RMSE error which drops to nearly zero across the entire image, with more
error incurred at the contrast boundaries. Column 3 shows the symmetric case which has significantly higher RMSE across the image than LDDMM.  






To compare homogeneous regions to boundary regions, we show RMSE in Figs. <ref>, <ref> as a function of noise at several regions identified in Fig. <ref> (bottom left). 
Results for each of these regions are shown; 
Fig. <ref> examines white noise; Fig. <ref> examines correlated noise.


The symmetric method degrades at higher noise levels, as well as with correlated noise.  

In the homogeneous regions the LDDMM has no uncertainty, and as one moves away from gradients the uncertainty quickly drops to zero.
Location (a) illustrates a place of high gradient on the boundary in which we see at reasonable noise levels equivalent performance.
The effect of noise correlation is to increase RMSE, the symmetric method being particularly harshly affected.







 §.§ Cramer-Rao bound

Examples of the Cramer-Rao bound on the divergence of v̂_n for different n are shown in the top row of Fig. <ref>.  One observers an asymmetry between uncertainty tangent and normal to level lines, which has been noted by other authors <cit.> and is implied by our consideration of the stabilizer of the diffeomorphism group.  In the bottom row  of Fig. <ref> we show the lower bound as it varies with changing balance between the image gradient (<ref>) term to the quadratic prior term (<ref>).  These results have an important implication for image mapping parameter selection in the presence of noise.  Note that when both terms are multiplied by the same constant, the bound increases linearly and the displayed image will look the same.  The behavior as we transition from “prior dominant” to “image dominant” can be seen as we move from left to right.


























§ DISCUSSION
















In this work we demonstrated that image registration algorithms which give comparable accuracy in terms of alignment of observable anatomical boundaries, can lead to widely different estimates of the canonical volume form.  In particular we showed that by employing the non-symmetric procedure associated to geodesic matching of LDDMM, we have explicitly accounted for the stabilizer of the diffeomorphism group, reducing the dimensionality of mappings to that of the observable image boundaries, results in favorable performance in the presence of noise.
Interestingly the small perturbation argument of
Equation (<ref>) also illustrates the importance of the asymmetry, and shows that it is instructive to consider the effect of image variability on estimation of the canonical volume form in 2 scenarios.  
First, because of the null-space of identifiability of motions along level lines, the homogeneous regions that are distant from image gradients by an amount larger than the spatial scale of K have clear identifiability issues in the symmetric methods. This leads to favorable performance of the asymmetric LDDMM method.
Second, boundaries between anatomical structures, such as gray white matter interfaces, have high image gradients.  In these cases the δ v is approximately related to δ J through the pseudoinverse of (∇ I) (∇ I)^T.  

The simulation of noise in the images clearly illustrates that
the departure between source and channel in Shannon's classic model is at the heart of the random orbit model, and implies the importance of asymmetry.
In this setting, the targets are not in the homogeneous orbit of the template under diffeomorphisms. The prior distribution and the space associated to it is separate from the targets as outputs of the noise channel. 
Noisy targets have gradients which are nonzero everywhere, which explains the clear superiority of the asymmetric approach of LDDMM in this situation.

This fills an important gap in our knowledge of how image registration algorithms can or should be used to quantify biological properties of tissue.



In brain morphomety <cit.>, the canonical volume form at each voxel in a 3D image is used to quantify structural differences between populations, for example measuring atrophy due to disease or aging.  In microscopy it is convenient to idealize cellular locations in terms of spatial point process models <cit.>.  Here the canonical volume used to define a relationship between counts (at the micron scale) and densities (at the mm scale) in a standard atlas coordinate system. During the earliest phases of preclinical Alzheimer's disease for example, when cells have not yet died, cell density changes due to dissolution of the neuropil.


When choosing a mapping procedure, one must consider that algorithms or parameters that may be optimal for registration accuracy, may be considerably sub optimal for quantifying local expansion or contraction.

















































One challenge here is in considering the effect of the many possible parameter choices that define mapping algorithms. The parameters considered here represent realistic choices for neuroimaging applications, and were chosen based on experience in studying brain morphometry.  An important parameter not explicitly varied was the spatial scale of the transformation's smoothness.  Claims made for “homogeneous regions” or “regions near a boundary”, should be understood as related to this spatial scale.  The MR image considered for our experiment was chosen because it has features at many different spatial scales, from the closely packed structures of the dentate gyrus, to the broadly uniform white matter of the angular bundle.  We have observed experimentally that as the spatial scale of the deformation grows very large, every location in the image becomes “near a boundary”, and the differences between the two methods considered become less pronounced.











One limitation of this work is that the experiments considered only two dimensional images.  This was for reasons of computational speed when considering statistical ensembles of mapping results, as well as for ease of display and interpretation.  The theoretical developments here do hold in three dimensions, and we expect our findings to generalize.




















This work is important for two important biomedical applications.  The first is the effect of data noise on atlas registration procedures.  When considering areas such as surgical planning, the quantification of “average performance” may be insufficient.  The consequence of deviations from average in this context can incredibly costly. The second is the extension of these atlas mapping technologies to morphometry, and the study of cell and process densities. 

In the new era of computational anatomy enabled by large volumes of light microscopic data, quantification of the uncertainty in the coordinate mapping between individual data sets and templates is important both for fundamental science applications, and for applications to the study of pathological conditions or computer aided diagnoses.   






§ ACKNOWLEDGEMENTS

Funding: This work was supported by the National Institutes of Health [P41-EB015909, 
R01-EB020062,
R01-NS102670, 
U19-MH114821, and U01 MH114824 
];  the Kavli Neuroscience Discovery Institute, the Crick-Clay Professorship, CSHL, and the H N Mahabala Chair, IIT Madras.









§ EXPERIMENTAL PARAMETERS

Parameters used in our simulations are listed in Table <ref>
.  











§ PROOF OF LDDMM INEXACT MATCHING



We will need the perturbation of the inverse. 
Let ϕ^ϵ = ϕ + ϵδϕ, computing the variation
uses the fact that (ϕ + ϵδϕ) ∘ (ϕ^-1 + ϵδϕ^-1) = id giving

    δϕ^-1 = - (D ϕ)|_ϕ^-1^-1δϕ |_ϕ^-1 .


Define the total cost C(ϕ) for LDDMM to be minimized:

    C(ϕ) = ∫_0^1∫_X L(ϕ,ϕ̇)dx dt + U(ϕ_1)  , 
    
      U(ϕ_1) = 1/2σ^2∫_X |J-I∘ϕ_1^-1|^2 dx

and L(ϕ,ϕ̇) is given by
(<ref>).
The Euler-Lagrange equation follows from first order perturbation, ϕ→ϕ + ϵδϕ, ϕ̇↦ϕ̇+ ϵd/dtδϕ:

    d/d ϵ C(ϕ+ϵδϕ)|_ϵ=0 =
    ∫_0^1 ∫_X ∂ L (ϕ,ϕ̇)/∂ϕ·δϕ dx dt 
       + ∫_0^1 ∫_X ∂ L(ϕ,ϕ̇)/∂ϕ̇·d/dtδϕ dx dt  + d/dϵ U(ϕ_1+ϵϕ_1)|_ϵ=0 .

Integrating by parts gives two equations using zero-boundary δϕ_0=0:

    ∫_0^1 ∫_X ( ∂ L (ϕ,ϕ̇)/∂ϕ - d/dt∂ L(ϕ,ϕ̇)/∂ϕ̇) · 
    (δϕ) dx dt  = 0;
    
       ∫_X ∂ L(ϕ_1,ϕ̇_1)/∂ϕ̇_1·δϕ_1 dx + d/dϵ U(ϕ_1 + ϵϕ_1) |_ϵ=0  = 0 .

To calculate the variation of the endpoint term gives

    d/dϵ|_ϵ = 0 U(ϕ_1+ϵδϕ_1) 
       =  ∫_X 1/σ^2 (J - I ∘ϕ_1^-1) (D ϕ_1)_ϕ_1^-1^-1T (∇ I)_ϕ_1^-1·δϕ_1 |_ϕ_1^-1 dx
    
       =
    ∫_X1/σ^2 (J ∘ϕ_1 - I )  (D ϕ_1)^-1T∇ I  |D ϕ_1 | ·δϕ_1  dx
     .

Using the fact that the conjugate momentum is p= ∂ L/∂ϕ̇= [Av] ∘ϕ |Dϕ|, then
Equation (<ref>)
is the Euler-Lagrange equation;

    d/dt p_t   +(Dv_t)^T∘ϕ_t  p_t =0  ;

(<ref>) is the boundary matching term
reducing to

    ∫ (p_1 +1/σ^2 (J ∘ϕ_1 - I )  (D ϕ_1)^-1T∇ I )  |D ϕ_1 |·δϕ_1  dx = 0
     .

giving

    p_1 +1/σ^2 (J ∘ϕ_1 - I )  (D ϕ_1)^-1T∇ I   |D ϕ_1| = 0
     .


The result for p_t in terms of p_1 is obtained by applying the result of Appendix <ref> twice. First calculate at time 0: p_0 = Dϕ_1^T p_1, then at time t: p_t = Dϕ_t^-TDϕ_1^T p_1, giving the result

    p_t + 1/σ^2(J∘ϕ_1 - I)Dϕ_t^-T∇ I |Dϕ_1| = 0






§ EULER-LAGRANGE AND CONSERVATION ARE EQUIVALENT




    ṗ_t   +(Dv_t)^T∘ϕ_t  p_t  =0  ; 
     
           Dϕ_t^T p_t  =   p_0  .


To see this,
take the derivative

    d/dt D ϕ_t^T p_t    =
    D(v_t ∘ϕ_t)^T p_t + D ϕ_t^T d/dt p_t 
       = D ϕ_t^T Dv_t^T ∘ϕ_t p_t +D ϕ_t^T ṗ_t  .

The last equality is zero by Euler-Lagrange equation (<ref>).






















§ SMALL DEFORMATION MOMENTUM



Let ϕ^ϵ = ϕ + ϵδϕ (ϕ), then v^ϵ = v+ ϵδϕ (ϕ).
Computing the variation requires perturbation of the inverse, δϕ^-1 = - (D ϕ)|_ϕ^-1^-1δϕ |_ϕ^-1 determined in (<ref>).

Then we have

    d/d ϵ C(ϕ+ϵδϕ) |_ϵ = 0
       =  d/dϵ1/2 ∫_X A v^ϵ· v^ϵ dx +1/2 σ^2 J- I ∘ϕ^ϵ -1^2 |_ϵ=0
       =
    ∫_X Av ·δϕ (ϕ)   dx 
            +  ∫_X 1/σ^2 (J - I ∘ϕ^-1 )  D (ϕ^-1)^T ( ∇ I)_ϕ^-1·δϕ  dx
     .

Substituting into the second integral x=ϕ(y) with dx =|dϕ| dy gives
us (<ref>).





§ EFFECT OF SMALL PERTURBATIONS



We consider the perturbation J ↦ J + ϵδ J and v ↦ v + ϵδ v on  (<ref>), seeking a relationship between the two to first order in ϵ∈ℝ.

    d/dϵ[Av + ϵ Aδ v] |_ϵ = 0 = d/dϵ1/σ^2(I - (J + ϵδ J)∘(ϕ + ϵδ v)) 
            D[ϕ + ϵδ v]^-T∇ I |Dϕ + ϵδ dv|_ϵ = 0

The left hand side is simply Aδ v.  The right hand side includes a product of three terms.  The variation of the first gives - δ J(ϕ) - DJ(ϕ)δ v.  That of the second gives -[Dϕ]^-TDδ v^T [D ϕ]^-T, while the third gives div[δ v]|Dϕ|.  Combining these gives

    Aδ v    = -1/σ^2[DJ(ϕ)δ v + δ J(ϕ)][Dϕ]^-T∇ I |dϕ|
       - 1/σ^2[I - J(ϕ)][Dϕ]^-TDδ v^T [dϕ]^-T∇ I |dϕ|
        +1/σ^2[I - J(ϕ)][Dϕ]^-T|Dϕ|div[ϕ]


We consider one of the cases examined experimentally, with I = J and therefore ϕ = id (identity).  Only the first term is nonzero, giving

    Aδ v = 1/σ^2[DI δ v + δ J]∇ I

Rearranging gives the relationship (<ref>).







