The LDDMM framework proposes to compare shapes (meshes or images) using diffeomorphic transformations of the 2D or 3D ambient space between these objects. To use it, we need to be able to parametrize a large family of transformations of the 2D or 3D space as well as to be able to compute distances between the objects.

[Bône et al. 2018] provides a concise reference summarizing the theory behind Deformetrica applications. For more details, see [Durrleman et al. 2014].

## Parametric family of deformations

To generate deformations, deformetrica uses a set of `n`

control points `(q_i)_{i=1,\ldots,n}`

and of vectors `(\mu_i)_{i=1,\ldots,n}`

attached to these control points, called momenta. These are in turn used to generate a vector field `X`

on the full 2D or 3D space along the formula:

` X(x) = \sum_{i=1}^p K(x,q_i) \cdot \mu_i`

where `X(x)`

is the vector field evaluated at the position `x`

, and where `K(x,y) = \exp\left(-\frac{\|x-y\|^2}{\sigma^2}\right)`

is a Gaussian kernel of width `\sigma`

. This width `\sigma`

is a critical hyperparameter of any deformetrica run, to be specified in the `deformation-parameters`

of the `model.xml`

. It is to be understood as the typical width of the described deformation. With a large `\sigma`

, only very smooth and global deformations are generated whereas with a small `\sigma`

, the deformation can have a very thin level of precision. In any case, it is advised to **start with a low value for \sigma** and to increase it if need be.

A time dynamic is prescribed for the control points and the momenta, according to the Hamiltonian equations:

```
\begin{cases}
&\dot q(t) = K[q(t), q(t)] \cdot \mu(t)\\
&\dot \mu(t) = -\frac{1}{2}\nabla_{q}\left\{ K[q(t), q(t)] \cdot \mu(t)^T \mu(t)\right\}
\end{cases}
```

which correspond to the geodesic equation of the considered space of deformations. The integration of these equations is done using a simple Euler or a Runge-Kutta scheme. Each iteration requires to compute convolutions between the control points and the momenta, and is therefore quadratic in the number of control points. Below is an example of time-varying vector field generated this way:

Overall, a deformation of the ambient space is **fully parametrized by an initial set of control points and an initial set of momenta**. We denote this transformation `\Phi_{q,\mu}`

where `q,\mu`

are the initial control points and momenta. Consequently, each model will be formulated as an optimization problem on these variables.

Such a deformation **acts on meshes** by direct application to the vertices of the meshes, as shown below:

On images, seen as function `I`

from `\mathbb R^2`

(or `\mathbb R^3`

) to `\mathbb R`

, it acts via:

`\Phi_{q,\mu}(I) = I \circ \Phi_{q,\mu}^{-1}.`

The deformation of a mesh or an image under such a deformation is computed by convolution between the control points and their momenta and the vertices of the meshes or the voxels of the images. It has a complexity which is in `n \times d`

where `n`

is the number of control points and `d`

is the number of vertices/voxels of the shape. Consequently, on large objects, the use of GPU computations provide a sensible speed-up.

## Distances between shapes

Now that we can deform shapes, we need a criterion indicating how close two shapes are.

### Distances between meshes

A mesh is described by its vertices `(v_p)_{p=1,\ldots, d}`

and by its connectivity. The faces of `PolyLine`

objects are segments (see the `skulls`

examples). The faces of `SurfaceMesh`

objects are triangles (see the `brain_structures`

examples). In each case, it i possible to compute the centers `(c_p)_{p=1,\ldots, f}`

and the normals `(n_p)_{p=1,\ldots, f}`

of the edges.
Two kind of distances are considered between meshes:

#### Current distance

The current distance between two meshes of the same type is then given by:

`d\left((n_p^\alpha, c_p^\alpha)_{p=1,\ldots,f^\alpha}, (n_q^\beta, c_q^\beta)_{q=1,\ldots,f^\beta}\right)^2 = \sum_p \sum_q K^W(c_p^\alpha, c_q^\beta) ((n_p^\alpha)^\top n_q^\beta).`

where `K^W`

is a gaussian kernel with width `\sigma^W`

.

#### Varifold distance

The varifold distance between two meshes of the same type is then given by:

`d\left((n_p^\alpha, c_p^\alpha)_{p=1,\ldots,f^\alpha}, (n_q^\beta, c_q^\beta)_{q=1,\ldots,f^\beta}\right)^2 = \sum_p \sum_q K^W(c_p^\alpha, c_q^\beta) \frac{[(n_p^\alpha)^\top n_q^\beta]^2}{\|n_p^\alpha\| \|n_q^\beta\|}.`

where `K^W`

is a gaussian kernel with width `\sigma^W`

.

In both cases, `\sigma^W`

is another hyper-parameter to be set for each mesh object in the `model.xml`

file. It depends on the level of fit that is desired in the computation. A small value will demand the fit to be very close while a large value will only demand the fit to be globally good. It is advised to start with large values for this width, typically a fifth of the length of the meshes, before re-running with lower values of needed.

### Distances between images

Only the `\ell^2`

distance is implemented in deformetrica: it is the sum of the squares of the differences of two images voxel-wise.

### Multi-objects distances

When the template object is composed of several objects, the distances between two observations is a weighted sum of the distances between each pair of corresponding objects. For instance, if each multi-object is composed of an hippocampus with noise-std `\varepsilon_h`

and of an amygdala with noise-std `\varepsilon_a`

, then the distance between two multi-objects `(h_\alpha, a_\alpha)`

, `(h_\beta, a_\beta)`

is:

`d((h_\alpha, a_\alpha), (h_\beta, a_\beta))^2 = \frac{1}{\varepsilon_h^2}d(h_\alpha, h_\beta)^2 + \frac{1}{\varepsilon_a^2} d(a_\alpha, a_\beta)^2`

## Loss functions for the models

Each model can be formulated as an optimization problem on a loss function.

### Registration

For a simple registration between two shapes `C_1, C_2`

, the loss function is:

`f(q,\mu) = d\left(\Phi_{q,\mu}(C_1), C_2)\right)^2 + R(q, \mu)`

where the first term measures data-attachment i.e. how well the deformation of the first shape corresponds to the second shape, whereas the second term is the norm of the deformation. This second term acts as a regularizer requiring the deformation to have low energy. The trade-off between attachment and regularization is managed through the **noise-std** f each object, as detailed in the model.xml file page.

### Deterministic Atlas

This model computes the mean of a set of observations, and finds momenta so that each observation is well reconstructed by the mean deformed along these momenta. If we denote `(C_i)_{i=1,\ldots,n}`

the observations, the loss function is:

`f(T, q,(\mu_i)_{i=1,\ldots,n}) = \sum_{i=1}^n \left( d\left(\Phi_{q,\mu_i}(T), C_i)\right)^2 + R(q, \mu_i) \right)`

The optimization is done on `T, q`

and `(\mu_i)_{i=1,\ldots,n}`

, although `q`

is fixed by default in deformetrica.

### Geodesic regression

Geodesic regressions seeks the geodesic which best fits a collection of observations `(C_i)_{i=1,\ldots,n}`

observed at times `(t_i)_{i=1,\ldots,n}`

. The loss function is:

`f(T,q, \mu) = \sum_{i=1}^n d\left(\Phi_{q, t_i \mu}(T), C_i \right)^2 + R(q, \mu)`

### Principal geodesic analysis

As for the other models, we have here a template shape `T`

and a set of control points `c`

. We suppose that we have a collection of shapes `(C_i)_{i=1,\ldots,n}`

.

Principal geodesic analysis estimated `d`

principal directions of shape variability. These are stored in the columns of a matrix `A`

: each column contains the momenta describing a principal direction of variability in the diffeomorphism space.

Then, each observed shape `C_i`

is given a latent coordinate `l_i`

which is such that `S_i`

is best reconstructed via diffeomorphic transformation of a template `T`

along the diffeorphism with initial momenta `Al`

: a linear combination of the principal directions described in `A`

. Finally, the optimized cost function is the sum of the reconstruction errors for each shape added to a regularization term encouraging the latent coordinates to be distributed along a zero-centered normal distribution, as for PCA:

`f(T, c, A, (l_i)_i) = \sum_i \left( d\left(\Phi_{c, Al_i}\right)^2 + R(l_i) \right).`

where `R(l_i)`

is a log likelihood of the unit-variance zero-mean normal distribution. Once again, the noise-variance parameters included in the distance `d`

control the trade-off regularization/reconstruction and the trade-off between the different objects in each observation.

Finally, note that the principal geodesic analysis model in deformetrica starts with a deterministic atlas whose output is used to initialize `T, c C_i, l_i`

by tangent PCA.