====================== Registration framework ====================== Deformation modules =================== Mathematical definition ----------------------- A deformation module is a structure :math:`(\mathcal O, H, \zeta, \xi, c)` where :math:`\mathcal O` is the space of **geometrical descriptors**, :math:`H` is the space of **controls**, :math:`\zeta : \mathcal O \times H \mapsto C^\ell_0 (\mathbb R^d, \mathbb R^d)` is the **field generator**, :math:`\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 :math:`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 :math:`M^i = (\mathcal O^i, H^i, \zeta^i, \xi^i, c^i)`, :math:`1 \leq i \leq N`, to build the compound module :math:`(\mathcal O, H, \zeta, \xi, c)` defined by :math:`\mathcal O = \prod_{i=1}^N \mathcal O^i`, :math:`H = \prod_{i=1}^N H^i`, :math:`\zeta: ( q,h) \in \mathcal O \times H \mapsto \sum_i \zeta^i_{q^i} h^i`, :math:`\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) )`and :math:`c: ( q,h) \in \mathcal O \times H \mapsto \sum_i c^i_{q^i} h^i` with :math:`q= (q_1, \dots, q_N)`, :math:`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 :math:`(\mathcal O, H, \zeta, \xi, c)`. It is then parametrized by a trajectory :math:`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 :math:`t \in [0,1] \mapsto v_t \doteq \zeta_{q_t} (h_t)` can be integrated: the flow :math:`\varphi^{\zeta_q (h)}` defined by :math:`\varphi_{t=0}^{\zeta_q(h)} = Id` and :math:`\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 :math:`q_0 \in \mathcal O` and *momentum* :math:`p_0 \in T_{q_0} ^\ast \mathcal O` and satisfy the following *shooting equation* .. math:: \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} where .. math:: 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 :math:`Z_q : H \mapsto H^\ast` is such that :math:`(Z_q (h), h) = c_q (h)` for all :math:`(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 :math:`q` and an attribute **cotan** implementing the momentum :math:`p` which are used as initial conditions for the shooting equation. The geodesic control :math:`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 :math:`(f_S, f_T)` in a shape space :math:`S` and a deformation module :math:`M = (\mathcal O, H, \zeta, \xi, c)`, the goal is to estimate the *best* (geodesic) modular large deformation generated by :math:`M` transporting :math:`f_S` as close as possible to :math:`f_T`. This amounts to minimizing the registration functional :math:`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 :math:`((q, q^S), (p p^S))` satisfying the shooting equation for the combination of :math:`M` and the silent deformation module induced by :math:`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 :math:`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 :math:`K`. We only consider for the moment scalar Gaussian kernel :math:`K_\sigma` with :math:`\sigma >0`. Let :math:`P` be an fixed integer, the deformation module generating a sum of :math:`P` translations is defined by: * :math:`\mathcal O = (\mathbb R^d)^P` * :math:`H = (\mathbb R^d)^P` * :math:`\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 :math:`q = (x_1, \dots, x_P)` and :math:`h = (h_1, \dots, h_P)` * :math:`\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 :math:`q = (x_1, \dots, x_P)` * :math:`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 :math:`\sigma` (scale of the scalar Gaussian kernel), an integer :math:`d` (dimension of the ambient space) and an integer :math:`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 :math:`N` be an integer and :math:`f : (\mathbb R^d)^N \mapsto (\mathbb R^d)^P`, :math:`g: (\mathbb R^d)^N \mapsto (\mathbb R^d)^P` two functions, we define the corresponding deformation module by * :math:`\mathcal O = (\mathbb R^d)^N` * :math:`H = \mathbb R` * :math:`\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 :math:`f(q) = (f(q)_1, \dots, f(q)_P)` and :math:`g(q) = (g(q)_1, \dots, g(q)_P)` * :math:`\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 :math:`q = (x_1, \dots, x_N)` * :math:`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 :math:`\sigma` (scale of the scalar Gaussian kernel), an integer :math:`d` (dimension of the ambient space), an integer :math:`N` (number of points for the geometrical descriptor) and two functions :math:`f\_support`, :math:`f\_vectors` implementing respectively the functions :math:`f` and :math:`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 :math:`\sigma` (scale of the scalar Gaussian kernel) and an integer :math:`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 :math:`\mathcal S` be a shape space, the silent deformation module induced by :math:`\mathcal S` is defined by: * :math:`\mathcal O = \mathcal S` * :math:`H = \{ 0 \}` * :math:`\zeta: (q,h) \in \mathcal O \times H \mapsto 0` * :math:`\xi: (q,v) \in \mathcal O \times C^\ell (\mathbb R^d, \mathbb R^d) \mapsto \xi_q (v)` (defined by the shape space :math:`\mathcal S`) * :math:`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 :math:`d` (dimension of the ambient space) and an integer :math:`N` (number of points for the geometrical descriptor). Compound deformation module --------------------------- Mathematical definition ^^^^^^^^^^^^^^^^^^^^^^^ This module :math:`(\mathcal O, H, \zeta, \xi, c)` is built from a collection of :math:`1 \leq i \leq N` deformation modules :math:`M^i = (\mathcal O^i, H^i, \zeta^i, \xi^i, c^i)` and is defined by: * :math:`\mathcal O = \prod_{i=1}^N \mathcal O^i` * :math:`H = \prod_{i=1}^N H^i` * :math:`\zeta: ( q,h) \in \mathcal O \times H \mapsto \sum_i \zeta^i_{q^i} h^i` * :math:`\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) )` * :math:`c: ( q,h) \in \mathcal O \times H \mapsto \sum_i c^i_{q^i} h^i` with :math:`q= (q_1, \dots, q_N)`, :math:`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 :math:`\zeta` is defined as the solution of a minimization problem: .. math:: \zeta_q (h) = argmin \{ |v|_V^2 + \frac{1}{\nu}|S_q (v) - A_q (h)|^2 \} with - :math:`V` a RKHS of vector field, - :math:`\nu` a positive scalar - :math:`S_q : V \to Y` a linear surjective constraint operator on vector fields that takes values in the space of constraints :math:`Y` (vector space of finite dimension) - :math:`A_q : H \to Y` is a linear operator which defines the value that one wants to observe. The space of geometrical descriptors :math:`\mathcal O` and controls :math:`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 :math:`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 .. math:: \zeta_q (h) = KS_q^*\lambda with .. math:: \lambda=(S_qKS_q^*+\nu I)^{-1}A_q(h) and :math:`K` the kernel of the RKHS :math:`V`. Order 0 ------- In this first example, the geometrical descriptor :math:`q` is made of points and the constraint operator :math:`S_q` returns the values of the input vector field on these points. Mathematical definition ^^^^^^^^^^^^^^^^^^^^^^^ Let :math:`K_\sigma` be a scalar Gaussian kernel with :math:`\sigma >0`, :math:`V` be the corresponding RKHS, :math:`P` be an fixed integer and :math:`\nu >0`. We define: * :math:`\mathcal O = (\mathbb R^d)^P` * :math:`H = (\mathbb R^d)^P` * :math:`Y = (\mathbb R^d)^P` * :math:`S_q: v \in V \mapsto (v(x_1), \dots, v(x_P))` with :math:`q = (x_1, \dots, x_P)` * :math:`A_q : h \in H \mapsto (K_q + \nu Id ) h` The field generator is then given by .. math:: \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 :math:`q = (x_1, \dots, x_P)` and :math:`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 :math:`\sigma` (scale of the scalar Gaussian kernel), an integer :math:`d` (dimension of the ambient space), an integer :math:`P` (number of local translations) and a scalar :math:`\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 :math:`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 :math:`\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 :math:`x_i` this infinitesimal strain tensor to be equal to the *growth factor* defined by a choice of eigenvectors :math:`R_i` and eigenvalues :math:`\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 :math:`\alpha_i` for each :math:`i` is the way to define the imposed deformation structure. Structural relations among the different :math:`\alpha_i` s can be captured saying that :math:`\alpha_i = \alpha_i (h)` where :math:`h` is a control parameter. The operator :math:`C : h \in H \mapsto (\alpha_1 (h), \dots, \alpha_N (h) ) \in (\mathbb R^d)^N` is called the *growth model operator*. Let :math:`K_\sigma` be a scalar Gaussian kernel with :math:`\sigma >0`, :math:`V` be the corresponding RKHS, :math:`N` be an fixed integer, :math:`P` be an fixed integer, :math:`\nu >0` and :math:`N` linear operators :math:`\alpha_i : \mathbb R^P \mapsto \mathbb R^d`. The corresponding implicit deformation module of order 1 is defined by: * :math:`\mathcal O = \{ q = ( (x_i, R_i) )_{1 \leq i \leq N} \in (\mathbb R^d \times SO_d)^N \}` * :math:`H = \mathbb R^P` * :math:`Y = S (\mathbb R^d)^N` * :math:`S_q: v \in V \mapsto (\epsilon_{x_1}(v), \dots, \epsilon_{x_N}(v))` with :math:`q = ( (x_i, R_i) )_{1 \leq i \leq N}` * :math:`A_q : h \in H \mapsto (R_1 D_1(h) R_1^T, \dots, R_N D_N(h) R_N^T)` with :math:`q = ( (x_i, R_i) )_{1 \leq i \leq N}` and for each :math:`i` :math:`D_i (h) = diag (\alpha_i (h) )` (diagonal matrix). Implementation ^^^^^^^^^^^^^^ This module is implemented in the class **ImplicitModule1** which is initialized by a scalar :math:`\sigma` (scale of the scalar Gaussian kernel), an integer :math:`d` (dimension of the ambient space), an integer :math:`N` (number of points on which the constraints are imposed), an integer :math:`p` (dimension of the control), a scalar :math:`\nu >0` and a **growth model tensor** :math:`C` of size :math:`N \times d \times P` (implementing the growth model operator). .. toctree:: :maxdepth: 2 :caption: Modules