Registration framework

Deformation modules

Mathematical definition

A deformation module is a structure \((\mathcal O, H, \zeta, \xi, c)\) where \(\mathcal O\) is the space of geometrical descriptors, \(H\) is the space of controls, \(\zeta : \mathcal O \times H \mapsto C^\ell_0 (\mathbb R^d, \mathbb R^d)\) is the field generator, \(\xi : (q,v) \in \mathcal O \times C^\ell (\mathbb R^d, \mathbb R^d) \mapsto \xi_q (v) = v \cdot q \in T \mathcal O\) is the infinitesimal action and \(c : \mathcal O \times H \mapsto \mathbb R^+\) is the cost.

The field generator specifies the vector fields that can be generated by the deformation module, i.e. the structure that is incorporated in the deformation model.

It is possible to combine several deformation modules \(M^i = (\mathcal O^i, H^i, \zeta^i, \xi^i, c^i)\), \(1 \leq i \leq N\), to build the compound module \((\mathcal O, H, \zeta, \xi, c)\) defined by \(\mathcal O = \prod_{i=1}^N \mathcal O^i\), \(H = \prod_{i=1}^N H^i\), \(\zeta: ( q,h) \in \mathcal O \times H \mapsto \sum_i \zeta^i_{q^i} h^i\), \(\xi : (q,v) \in \mathcal O \times C^\ell (\mathbb R^d, \mathbb R^d) \mapsto (\xi_q^1 (v), \dots, \xi_q^1 (v) )\) with \(q= (q_1, \dots, q_N)\), \(h=(h_1, \dots, h_N)\).

Implementation

A deformation module is defined as a class whose main method is a field generator which specifies, for each module, the vector fields that can be generated. It uses the manifold and the control of the deformation modules and returns a structured vector field. The manifold is itself a class implementing the notion of geometrical descriptor and, more generally, of all geometrical objects which can be displaced through the action of a vector field. This action is defined by the infinitesimal action method. The control variable is a tensor. Structured vector fields implement the mathematical notion of vector field: they can be applied to a tensor of points and return a tensor of speeds. The second key method for a deformation module is the cost which associates a scalar cost to a couple of a manifold and a control.

The combination of deformation module is implemented via the class CompoundModule which is built from a list of deformation modules.

Modular large deformations

Mathematical definition

A modular large deformation is a diffeomorphism built as a flow of modular vector fields, i.e. a trajectory of vector fields generated by a given deformation module \((\mathcal O, H, \zeta, \xi, c)\). It is then parametrized by a trajectory \(t \in [0, 1] \mapsto (q_t, h_t) \in \mathcal O \times H\) of geometrical descriptor and control. Under some conditions, the corresponding trajectory of field \(t \in [0,1] \mapsto v_t \doteq \zeta_{q_t} (h_t)\) can be integrated: the flow \(\varphi^{\zeta_q (h)}\) defined by \(\varphi_{t=0}^{\zeta_q(h)} = Id\) and \(\dot{\varphi}_t^{\zeta_q(h)} = \zeta_q(h) \circ \varphi^{\zeta_q(h)}\) exists.

We are interested in particular trajectories: the geodesic ones. They are parametrized by an initial value of geometrical descriptor \(q_0 \in \mathcal O\) and momentum \(p_0 \in T_{q_0} ^\ast \mathcal O\) and satisfy the following shooting equation

\[\begin{split}\begin{cases} \dot{q}_t & = & \frac{\partial H}{\partial p} (q_t, p_t, h_t) = \xi_{q_t} \circ \zeta_{q_t} (h_t) \\ \dot{p}_t & = & -\frac{\partial H}{\partial q} (q_t, p_t, h_t) \\ h_t & = & Z_{q_t}^{-1} (\xi_{q_t} \circ \zeta_{q_t})^\ast p_t \end{cases}\end{split}\]

where

\[H : (q, p, h) \in T^\ast \mathcal O \times H \mapsto (p | \xi_{q_t} \circ \zeta_q (h)) - \frac{1}{2} c_q(h)\]

is the Hamiltonian and \(Z_q : H \mapsto H^\ast\) is such that \((Z_q (h), h) = c_q (h)\) for all \((q,h) \in \mathcal O \times H\).

Implementation

The shoot class implements the shooting equations. It takes as input an Hamiltonian implementing the mathematical Hamiltonian which is itself defined from a deformation module. The manifold of the deformation module contains an attribute gd implementing the geometrical descriptor \(q\) and an attribute cotan implementing the momentum \(p\) which are used as initial conditions for the shooting equation. The geodesic control \(h_t = Z_{q_t}^{-1} (\xi_{q_t} \circ \zeta_{q_t}) ^\ast p_t\) is computed at each time using the method compute_geodesic_control of the deformation module.

The shoot class can also take in input a list of controls (implementing a trajectory of controls) which are used at each time instead of the geodesic ones. This can be used for instance to study the influence of one specific component of the controls.

Modular registration

Mathematical definition

Given a couple of source and target shapes \((f_S, f_T)\) in a shape space \(S\) and a deformation module \(M = (\mathcal O, H, \zeta, \xi, c)\), the goal is to estimate the best (geodesic) modular large deformation generated by \(M\) transporting \(f_S\) as close as possible to \(f_T\). This amounts to minimizing the registration functional \(J(q_0, p_0, p^S_0) \in T^\ast \mathcal O \times T^\ast_{f_S} \mathcal F \mapsto \int_0^1 c_{q_t} (h_t) d t + D (\varphi^{\zeta_q (h)}_{t=1} \cdot f_S , f_T)\) with \(((q, q^S), (p p^S))\) satisfying the shooting equation for the combination of \(M\) and the silent deformation module induced by \(S\).

Implementation

The class RegistrationModel is the keystone allowing to perform a registration using a chosen deformation model. It is initialized with lists of Deformable, Attachment (defining the attachment functions for these data) and DeformationModules. A Deformable is a class implementing the notion of geometric data: it gathers point clouds, meshes and images. To each element of the Deformable list, is associated an Attachment implementing a distance between to Deformables of the same class. Several Attachments are implemented such as euclidean distance for point clouds, the varifold distance for curves or meshes and \(L^2\) distance for images. The method compute_deformed allows to compute the shooting equation with the method compute_deformed, the initial value of geometrical descriptor and momentum being stored in the attribute init_manifold. It also enables to compute the registration function with the evaluate method which takes as input a list of Deformable targets. In order to minimize this functional, a Fitter class is created from RegistrationModel and runs the optimisation. It interfaces most PyTorch and Scipy optimizers along with a gradient descent with linear search optimization algorithm. It is possible to optimize the parameters of the deformation modules with a key-word other_parameters in RegistrationModel. It is also possible to add a callback function that is evaluated before model evaluation. This allows us to estimate some parameters of the model as functions of meta-parameters which can be estimated.

Implemented deformation modules

There are two categories of deformation modules: the explicit ones where the field generator is explicitly given in the definition, and the implicit ones

Explicit deformation modules

These deformation modules are said to be explicit in the sense that their field generator is explicitly given in the definition.

Local translations

Mathematical definition

This module generates a sum of local translations, localized by a chosen kernel \(K\). We only consider for the moment scalar Gaussian kernel \(K_\sigma\) with \(\sigma >0\).

Let \(P\) be an fixed integer, the deformation module generating a sum of \(P\) translations is defined by:

  • \(\mathcal O = (\mathbb R^d)^P\)

  • \(H = (\mathbb R^d)^P\)

  • \(\zeta: (q,h) \in \mathcal O \times H \mapsto \sum_{i=1}^P K(x_i, \cdot) h_i \in C^\ell (\mathbb R^d, \mathbb R^d)\) where \(q = (x_1, \dots, x_P)\) and \(h = (h_1, \dots, h_P)\)

  • \(\xi: (q,v) \in \mathcal O \times C^\ell (\mathbb R^d, \mathbb R^d) \mapsto (v(x_1), \dots, v(x_P)) \in T_q \mathcal O\) where \(q = (x_1, \dots, x_P)\)

  • \(c: (q,h) \in \mathcal O \times H \mapsto |\zeta_q (h)|^2_{V_\sigma} \in \mathbb R^+\)

Implementation

This module is implemented in the class Translations which is initialized by a scalar \(\sigma\) (scale of the scalar Gaussian kernel), an integer \(d\) (dimension of the ambient space) and an integer \(P\) (number of local translations).

Local constrained translations

Mathematical definition

This deformation module builds also sum of local translations but imposes some links between them to constrain the generated field. Let \(N\) be an integer and \(f : (\mathbb R^d)^N \mapsto (\mathbb R^d)^P\), \(g: (\mathbb R^d)^N \mapsto (\mathbb R^d)^P\) two functions, we define the corresponding deformation module by

  • \(\mathcal O = (\mathbb R^d)^N\)

  • \(H = \mathbb R\)

  • \(\zeta: (q,h) \in \mathcal O \times H \mapsto h \sum_{i=1}^P K(f(q)_i, \cdot) g(x)_i \in C^\ell (\mathbb R^d, \mathbb R^d)\) where \(f(q) = (f(q)_1, \dots, f(q)_P)\) and \(g(q) = (g(q)_1, \dots, g(q)_P)\)

  • \(\xi: (q,v) \in \mathcal O \times C^\ell (\mathbb R^d, \mathbb R^d) \mapsto (v(x_1), \dots, v(x_P)) \in T_q \mathcal O\) where \(q = (x_1, \dots, x_N)\)

  • \(c: (q,h) \in \mathcal O \times H \mapsto |\zeta_q (h)|^2_{V_\sigma} \in \mathbb R^+\)

Implementation

This module is implemented in the class LocalConstrainedTranslations which is initialized by a scalar \(\sigma\) (scale of the scalar Gaussian kernel), an integer \(d\) (dimension of the ambient space), an integer \(N\) (number of points for the geometrical descriptor) and two functions \(f\_support\), \(f\_vectors\) implementing respectively the functions \(f\) and \(g\). Two particular deformation modules generating local constrained translations are implemented: LocalScaling (generating one local scaling) and LocalRotation (generating one local rotation). They are initialized by a scalar \(\sigma\) (scale of the scalar Gaussian kernel) and an integer \(d\) (dimension of the ambient space). Their geometrical descriptor is made of one point.

Silent

Mathematical definition

This module generates a null vector field: when it is combined with other ones, it does not contribute to the generated vector field but its geometrical descriptors are deformed by it. Let \(\mathcal S\) be a shape space, the silent deformation module induced by \(\mathcal S\) is defined by:

  • \(\mathcal O = \mathcal S\)

  • \(H = \{ 0 \}\)

  • \(\zeta: (q,h) \in \mathcal O \times H \mapsto 0\)

  • \(\xi: (q,v) \in \mathcal O \times C^\ell (\mathbb R^d, \mathbb R^d) \mapsto \xi_q (v)\) (defined by the shape space \(\mathcal S\))

  • \(c: (q,h) \in \mathcal O \times H \mapsto 0\)

Implementation

The silent module induced by a landmark shape space is implemented in the class Silent. It is initialized by an integer \(d\) (dimension of the ambient space) and an integer \(N\) (number of points for the geometrical descriptor).

Compound deformation module

Mathematical definition

This module \((\mathcal O, H, \zeta, \xi, c)\) is built from a collection of \(1 \leq i \leq N\) deformation modules \(M^i = (\mathcal O^i, H^i, \zeta^i, \xi^i, c^i)\) and is defined by:

  • \(\mathcal O = \prod_{i=1}^N \mathcal O^i\)

  • \(H = \prod_{i=1}^N H^i\)

  • \(\zeta: ( q,h) \in \mathcal O \times H \mapsto \sum_i \zeta^i_{q^i} h^i\)

  • \(\xi : (q,v) \in \mathcal O \times C^\ell (\mathbb R^d, \mathbb R^d) \mapsto (\xi_q^1 (v), \dots, \xi_q^1 (v) )\)

  • \(c: ( q,h) \in \mathcal O \times H \mapsto \sum_i c^i_{q^i} h^i\) with \(q= (q_1, \dots, q_N)\), \(h=(h_1, \dots, h_N)\).

Implementation

The combination of deformation modules is implemented in the class CompoundModule initialized by a list of DeformationModule.

Implicit deformation modules

In many cases, the deformation prior may not be given in the form of an explicit type of vector field but more in some properties that the vector fields should satisfy. In order to address this problem, we define implicit deformation modules, where the field generator \(\zeta\) is defined as the solution of a minimization problem:

\[\zeta_q (h) = argmin \{ |v|_V^2 + \frac{1}{\nu}|S_q (v) - A_q (h)|^2 \}\]

with

  • \(V\) a RKHS of vector field,

  • \(\nu\) a positive scalar

  • \(S_q : V \to Y\) a linear surjective constraint operator on vector fields that takes values in the space of constraints \(Y\) (vector space of finite dimension)

  • \(A_q : H \to Y\) is a linear operator which defines the value that one wants to observe.

The space of geometrical descriptors \(\mathcal O\) and controls \(H\) will be specific to each type of implicit deformation module, we implemented two case: order 0 and oreder 1. The cost of the deformation module is then given by \(c_q (h) =|\zeta_q (h)|_V^2 + \frac{1}{\nu}|S_q (\zeta_q (h)) - A_q (h)|^2\).

From this definition, it is possible to compute the explicit expression for the field generator

\[\zeta_q (h) = KS_q^*\lambda\]

with

\[\lambda=(S_qKS_q^*+\nu I)^{-1}A_q(h)\]

and \(K\) the kernel of the RKHS \(V\).

Order 0

In this first example, the geometrical descriptor \(q\) is made of points and the constraint operator \(S_q\) returns the values of the input vector field on these points.

Mathematical definition

Let \(K_\sigma\) be a scalar Gaussian kernel with \(\sigma >0\), \(V\) be the corresponding RKHS, \(P\) be an fixed integer and \(\nu >0\). We define:

  • \(\mathcal O = (\mathbb R^d)^P\)

  • \(H = (\mathbb R^d)^P\)

  • \(Y = (\mathbb R^d)^P\)

  • \(S_q: v \in V \mapsto (v(x_1), \dots, v(x_P))\) with \(q = (x_1, \dots, x_P)\)

  • \(A_q : h \in H \mapsto (K_q + \nu Id ) h\)

The field generator is then given by

\[\zeta: (q,h) \in \mathcal O \times H \mapsto \sum_{i=1}^P K(x_i, \cdot) h_i \in C^\ell (\mathbb R^d, \mathbb R^d)\]

where \(q = (x_1, \dots, x_P)\) and \(h = (h_1, \dots, h_P)\). This implicit deformation module of order 0 is therefore a regularized version of the Local translations one.

Implementation

This module is implemented in the class ImplicitModule0 which is initialized by a scalar \(\sigma\) (scale of the scalar Gaussian kernel), an integer \(d\) (dimension of the ambient space), an integer \(P\) (number of local translations) and a scalar \(\nu >0\) (regularization parameter).

Order 1

The intuitive idea behind this module is to be able to incorporate objects related elastic properties in the deformation model. More specifically, we constrain the local changes of lengths induced by the action of the vector field at some specified locations \(x_i\) and along some specified directions attached to the object.

Mathematical definition

Constraining such local changes of lengths amounts to constraining the infinitesimal strain tensor \(\epsilon_{x_i}(v) \doteq \frac{D v (x_i) + Dv (x_i)^\ast}{2}\) which is a symmetric tensor capturing the local metric changes (expansion or dilation along given directions). We constrain at each point \(x_i\) this infinitesimal strain tensor to be equal to the growth factor defined by a choice of eigenvectors \(R_i\) and eigenvalues \(\alpha_i\). The eigenvectors correspond to geometric variables (directions) attached to the object and are then part of the geometrical descriptor. The choice of the eigenvalues \(\alpha_i\) for each \(i\) is the way to define the imposed deformation structure. Structural relations among the different \(\alpha_i\) s can be captured saying that \(\alpha_i = \alpha_i (h)\) where \(h\) is a control parameter. The operator \(C : h \in H \mapsto (\alpha_1 (h), \dots, \alpha_N (h) ) \in (\mathbb R^d)^N\) is called the growth model operator.

Let \(K_\sigma\) be a scalar Gaussian kernel with \(\sigma >0\), \(V\) be the corresponding RKHS, \(N\) be an fixed integer, \(P\) be an fixed integer, \(\nu >0\) and \(N\) linear operators \(\alpha_i : \mathbb R^P \mapsto \mathbb R^d\). The corresponding implicit deformation module of order 1 is defined by:

  • \(\mathcal O = \{ q = ( (x_i, R_i) )_{1 \leq i \leq N} \in (\mathbb R^d \times SO_d)^N \}\)

  • \(H = \mathbb R^P\)

  • \(Y = S (\mathbb R^d)^N\)

  • \(S_q: v \in V \mapsto (\epsilon_{x_1}(v), \dots, \epsilon_{x_N}(v))\) with \(q = ( (x_i, R_i) )_{1 \leq i \leq N}\)

  • \(A_q : h \in H \mapsto (R_1 D_1(h) R_1^T, \dots, R_N D_N(h) R_N^T)\) with \(q = ( (x_i, R_i) )_{1 \leq i \leq N}\) and for each \(i\) \(D_i (h) = diag (\alpha_i (h) )\) (diagonal matrix).

Implementation

This module is implemented in the class ImplicitModule1 which is initialized by a scalar \(\sigma\) (scale of the scalar Gaussian kernel), an integer \(d\) (dimension of the ambient space), an integer \(N\) (number of points on which the constraints are imposed), an integer \(p\) (dimension of the control), a scalar \(\nu >0\) and a growth model tensor \(C\) of size \(N \times d \times P\) (implementing the growth model operator).