Feature-Driven Warps

In this section, we specialize the generic Feature-Driven parameterization presented in section A.3.1 for two types of warps: the TPS and the FFD warps. Since the representational power of the TPS warp and of the FFD warp are equivalent (see experiments in section A.6.1), we focused our experiments on the TPS warp. However, it is important to show how the FFD warp can actually be used in the Feature-Driven framework. In particular, we show how the standard FFD model can be extended in order to be compatible with the warp reversion operation.

The Feature-Driven Thin-Plate Spline Warp

Definition

Ignoring the parameters, a TPS  $\bar{\omega}_{\scriptscriptstyle \mathsf{T}}$ is an  $\mathbb{R}^2 \rightarrow \mathbb{R}$ function. It is the Radial Basis Function that minimizes the integral bending energy. In its natural parameterization, a TPS is driven by a set of $l+3$ weights  $\bar{p}_k \in \mathbb{R}$. These weights are grouped in a vector of parameters  $\bar{\mathbf{p}} \in \mathbb{R}^{l+3}$. The evaluation of a TPS at the point  $\mathbf{q}^\mathsf{T}= (x \; y)$ is given by:

$$\bar{\omega}_{\scriptscriptstyle \mathsf{T}}(\mathbf{q} ; \bar{\m... ...{q} , \mathbf{c}_i)\right) + \bar{p}_{l+1} x + \bar{p}_{l+2}y + \bar{p}_{l+3}.$$ (A.23)

The $l$ 2D points  $\mathbf{c}_k$ are called the centers. They are also the driving features in the texture image. They can be located at any place but, in practice, we place them on a regular grid. The function $d^2$ gives the squared euclidean distance between its two arguments. The function $\rho$ is the TPS basis function and is defined by  $\rho(r) = r^2 \log(r)$ for $r > 0$ and  $\rho(0) = 0$. In matrix form, equation (A.23) is equivalent to:

$$\bar{\omega}_{\scriptscriptstyle \mathsf{T}}(\mathbf{q} ; \bar{\mathbf{p}}) = \boldsymbol{\ell}_{\mathbf{q}}^\mathsf{T}\bar{\mathbf{p}},$$ (A.24)

with $\boldsymbol{\ell}_{\mathbf{q}}^\mathsf{T}= \begin{pmatrix}\rho(d^2(\mathbf{q},... ...}, \mathbf{c}_l)) & \mathbf{q}^\mathsf{T}& 1 \end{pmatrix} \in \mathbb{R}^{l+3}$.

Standard $\mathbb{R}^2 \rightarrow \mathbb{R}^2$ TPS warps are obtained by replacing the scalar weights $\bar{p}_i$ by the control points  $\mathbf{p}_k = (p_k^x \; p_k^y) \in \mathbb{R}^2$. The control points are grouped in a single matrix of parameters  $\mathsf{P} \in \mathbb{R}^{(l+3) \times 2}$ defined by  $\mathsf{P}^\mathsf{T}= \left( \mathbf{p}_1 \; \ldots \; \mathbf{p}_{l+3} \right)$. The TPS warp is thus defined by:

$$\omega_{\scriptscriptstyle \mathsf{T}}(\mathbf{q} ; \mathbf{p}) =... ...{T}\mathsf{P}, \qquad \textrm{with } \mathbf{p} = \boldsymbol{\nu}(\mathsf{P}).$$ (A.25)

Feature-Driven Parameterization

The Feature-Driven parameterization of the TPS warp consists in replacing the control points by some features (i.e. points) in the current image. A point  $\mathbf{a}_k^\mathsf{T}= (a_k^x \; a_k^y) \in \mathbb{R}^2$ is assigned to each center  $\mathbf{c}_k$ defined in the texture image. The features  $\mathbf{a}_k$ are grouped in a single matrix $\mathsf{A} \in \mathbb{R}^{l \times 2}$. Similarly, the centers  $\mathbf{c}_k$ are grouped in a matrix  $\mathsf{C} \in \mathbb{R}^{l \times 2}$. Following (24), the control points of a TPS can be determined from the correspondences  $\mathbf{c}_k \leftrightarrow \mathbf{a}_k$:

$$\omega_{\scriptscriptstyle \mathsf{T}}(\mathbf{c}_k ; \mathbf{p}) = \mathbf{a}_k \qquad \forall k \in \{1, \ldots, l\},$$ (A.26)

while enforcing the 3 `side-conditions' ensuring that the TPS has square integrable second derivatives (more details can be found in (193)):

$$\sum_{i=1}^l x_i \mathbf{p}_i = \mathbf{0}_{2 \times 1}, \quad \... ...hbf{0}_{2 \times 1}, \quad \sum_{i=1}^l \mathbf{p}_i = \mathbf 1_{2 \times 1}.$$ (A.27)

Combining these $l+3$ conditions in a single matrix gives the following exactly determined linear system:

$$\mathsf{M}_\lambda \mathsf{P} = \begin{pmatrix} \mathsf{A} \mathbf{0}_{\scriptscriptstyle 3 \times 1} \end{pmatrix},$$ (A.28)

with  $\mathsf{M}_\lambda \in \mathbb{R}^{(l+3) \times (l+3)}$ the matrix defined by:

$$\mathsf{M}_\lambda = \begin{pmatrix} \mathsf{N}_\lambda & \math... ...athsf{Q}^\mathsf{T}& \mathbf{0}_{\scriptscriptstyle 3 \times 3} \end{pmatrix},$$ (A.29)

with $\mathsf{N}_\lambda = \mathsf{N} + \lambda \mathsf{I}_l$, $\mathsf{N}^\mathsf{T}= \left( \boldsymbol{\ell}_{\mathbf{c}_1}^\mathsf{T}\quad \cdots \quad \boldsymbol{\ell}_{\mathbf{c}_l}^\mathsf{T}\right)$ and  $\mathsf{Q} = \begin{pmatrix}\mathsf{C} & \mathbf{1}_{\scriptscriptstyle l\times 1} \end{pmatrix} \in \mathbb{R}^{l \times 3}$. Adding  $\lambda \mathsf{I}_l$ to $\mathsf{N}$ acts as a regularizer. Determining the control points  $\mathsf{P}$ from the equation (A.28) can be done in a straightforward manner as the solution of an exactly determined linear system. The resulting matrix of control points, denoted  $\mathsf{P}_\lambda$, is a nonlinear function of the regularization parameter $\lambda$ and a linear function of the features  $\mathsf{A}$:

$$\mathsf{P}_\lambda = \mathsf{M}_\lambda^{-1} \begin{pmatrix} \mathsf{A} \mathbf{0}_{\scriptscriptstyle 3 \times 1} \end{pmatrix}.$$ (A.30)

$\mathsf{P}_\lambda$ is a linear `back-projection' of the feature matrix  $\mathsf{A}$. It can be computed efficiently using the blockwise matrix inversion formulas:

$$\mathsf{P}_\lambda = \mathsf{E}_\lambda \mathsf{A}$$ (A.31)

with:

$$\mathsf{E}_\lambda = \begin{pmatrix} \mathsf{N}_\lambda^{-1}\le... ...hsf{Q} \right)^{-1}\mathsf{Q}^\mathsf{T}\mathsf{N}_\lambda^{-1} \end{pmatrix}.$$ (A.32)

This expression has the advantages of separating $\lambda$ and  $\mathsf{A}$ and introduces units: while  $\mathsf{P}_\lambda$ has no obvious unit, $\mathsf{A}$ in general has (e.g. pixels, meters). Finally, if we replace the natural parameters  $\mathbf{p}$ in the definition of the TPS warp  $\omega_{\scriptscriptstyle \mathsf{T}}$ (equation (A.25)) by their expression given in the equation (A.32), we get the Feature-Driven parameterization of the TPS warp, denoted  $\mathcal{W}_{\scriptscriptstyle \mathsf{T}}$:

$$\mathcal{W}_{\scriptscriptstyle \mathsf{T}}(\mathbf{q} ; \mathbf{... ...da \mathsf{A}, \qquad \textrm{with } \mathbf{a} = \boldsymbol{\nu}(\mathsf{A}).$$ (A.33)

We use the notation  $\mathcal{W}_{\scriptscriptstyle \mathsf{T}}(\mathbf{q} ; \mathbf{a})$ for $\mathcal{W}_{\scriptscriptstyle \mathsf{T}}(\mathbf{q} ; \mathbf{a}, 10^{-4})$. We choose $\lambda = 10^{-4}$ to ensure good numerical conditioning of the matrix  $\mathsf{N}_\lambda$.

Jacobian matrix of the warp.

The Jacobian matrix of the warp is needed by the Gauss-Newton based algorithms for local registration (see e.g. , section A.2.1 or section A.4.1). We denote  $\mathsf{K}_{\scriptscriptstyle \mathsf{T}}$ the Jacobian matrix of the TPS warp evaluated at the point  $\mathbf{q}$. It is defined by  $\mathsf{K}_{\scriptscriptstyle \mathsf{T}}= \frac{\partial \mathcal{W}_{\scrip... ...{\partial \mathbf{a}}(\mathbf{q} ; \mathbf{a}) \in \mathbb{R}^{2 \times 2(l+3)}$ and is given by:

$$\mathsf{K}_{\scriptscriptstyle \mathsf{T}} = \begin{pmatrix} \frac{\partial \mathcal{W}_{\scriptscriptstyl... ...athsf{T}}^y}{\partial \mathbf{a}^y}(\mathbf{q} ; \mathbf{a}) \end{pmatrix}$$

$$= \begin{pmatrix} \boldsymbol{\ell}_{\mathbf{q}}^\mathsf{T}\mat... ...} & \boldsymbol{\ell}_{\mathbf{q}}^\mathsf{T}\mathsf{E}_\lambda \end{pmatrix},$$ (A.34)

where $\mathcal{W}_{\scriptscriptstyle \mathsf{T}}^x$ and  $\mathcal{W}_{\scriptscriptstyle \mathsf{T}}^y$ are the first and the second coordinates of the warp  $\mathcal{W}_{\scriptscriptstyle \mathsf{T}}$ and  $\mathsf{A} = \left(\mathbf{a}^x \; \mathbf{a}^y\right)$.

The Feature-Driven Free-Form Deformation

Tensor-product B-Splines are a particular model of Free-Form Deformations. They are a general model of polynomial functions which have been proved to be useful for image registration (162). Even if there is a wide variety of B-Splines (with various degrees for the polynomial basis or by choosing exotic knot sequences), we limit our study to the case of the Uniform Cubic B-Splines since it best matches the needs of image registration. For the sake of simplicity, we will abbreviate it FFD.

Definition

Monodimensional case.

Ignoring the parameters, a monodimensional FFD  $\bar{\omega}_{\scriptscriptstyle \mathsf{F}}$ is an  $\mathbb{R}\rightarrow \mathbb{R}$ function defined as a linear combination of the basis functions $N_i$ weighted by the scalars $\bar{p}_k$ called the weights:

$$\bar{\omega}_{\scriptscriptstyle \mathsf{F}}(x ; \bar{\mathbf{p}}) = \sum_{i=1}^m \bar{p}_i N_i(x),$$ (A.35)

where  $\bar{\mathbf{p}} \in \mathbb{R}^m$ is the vector that contains all the weights $\bar{p}_k$. The basis functions are defined using a knot sequence, i.e. a non-decreasing sequence  $k_1 < \ldots < k_{m+4}$. The FFD is said to be uniform when the knot sequence is uniform, i.e. all the knot intervals  $[k_i,k_{i+1}]$ have the same length $s$. In this case, the basis functions $N_i$ are defined by using four polynomials of degree three, the blending functions (see figure A.8 for an illustration):

$$N_i(x) = \left\{ \begin{array}{l} b_1(x) = \frac{1}{6} \hat x... ...} 0 \qquad\qquad\qquad\qquad\qquad\textrm{otherwise} \end{array} \right.$$ (A.36)

where $\hat x$ is the normalized abscissa of $x$ defined as  $\hat x = \frac{x - k_I}{s}$ for $x \in [k_I,k_{I+1}]$.

Domain.

We can see from equation (A.35) that an FFD is non-zero only over the interval  $[k_1, k_{m+4}]$. However, it is common practice to reduce the domain to  $[k_4,k_{m+1}]$. By doing so, there are always exactly 4 non-zero basis functions on each knot interval, as figure A.8 illustrates.

Figure A.8: (a) The basis functions of an FFD are bell-shaped curves with bounded support. They are defined using 4 polynomial pieces of degree three: $b_1$, $b_2$, $b_3$ and $b_4$. (b) The usual (or natural) definition domain of an FFD is represented by the non-grayed part.

(a)	(b)

FFD warp.

The standard $\mathbb{R}^2 \rightarrow \mathbb{R}^2$ FFD warp is obtained as the two-way tensor-product of monodimensional FFD s. Using its natural parameterization, the evaluation of the FFD warp  $\omega_{\scriptscriptstyle \mathsf{F}}$ at the point  $\mathbf{q} = \left( x \; y \right)^\mathsf{T}$ is given by:

$$\omega_{\scriptscriptstyle \mathsf{F}}(\mathbf{q} ; \mathbf{p}) =... ...}^n \sum_{i=1}^m \mathbf{p}_k N_i(x) N_j(y), \quad \textrm{with } k = (j-1)m+i.$$ (A.37)

The $mn$ control points  $\mathbf{p}_k$ are grouped in the vector  $\mathbf{p} \in \mathbb{R}^{2mn}$ that is defined as  $\mathbf{p} = \boldsymbol{\nu}(\mathsf{P})$ where  $\mathsf{P} \in \mathbb{R}^{mn \times 2}$ is the matrix given by  $\mathsf{P}^\mathsf{T}= \left( \mathbf{p}_1 \; \ldots \; \mathbf{p}_{mn} \right)$. The control points of an FFD warp are not more meaningful than the ones of a TPS warp. They are not interpolated: they just act as `attractors' to the warp. Equation (A.37) can be rewritten in matrix form:

$$\omega_{\scriptscriptstyle \mathsf{F}}(\mathbf{q} ; \mathbf{p}) = \boldsymbol{\ell}_{\mathbf{q}}^\mathsf{T}\mathsf{P},$$ (A.38)

where $\boldsymbol{\ell}_{\mathbf{q}} \in \mathbb{R}^{mn}$ is the vector defined by:

$$\boldsymbol{\ell}_{\mathbf{q}}^\mathsf{T}= \left( N_1(x)N_1(y) \:\ldots\: N_m(x)N_1(y) \:\ldots\: N_m(x) N_n(y) \right).$$ (A.39)

Feature-Driven Parameterization

The Feature-Driven parameterization of the FFD warp is similar to the one of the TPS warp in the sense that it makes the warp driven by features expressed in pixels in both the texture and the current images. The centers of the TPS warp were used as features in the texture image. Such centers do not exist for FFD warps. We thus introduce a set of points $\mathbf{c}_k$ that will be used as features in the texture image. We call these points centers for consistency with the TPS warps. We use $l=mn$ centers located on a regular grid. A feature  $\mathbf{a}_k$ in the current image is associated to every center  $\mathbf{c}_k$. The control points  $\mathbf{p}$ of the FFD warp can be determined from the correspondences  $\mathbf{c}_k \leftrightarrow \mathbf{a}_k$ by enforcing the following constraints:

$$\omega_{\scriptscriptstyle \mathsf{F}}(\mathbf{c}_k ; \mathbf{p}) = \mathbf{a}_k, \qquad \forall k \in \{1,\ldots,l\}.$$ (A.40)

Since the number of features is equal to the number of degrees of freedom of the FFD warp, the determination of the parameters from the features can be carried out with an exactly determined linear system:

$$\mathsf{M} \mathsf{P} = \mathsf{A},$$ (A.41)

with $\mathsf{M}^\mathsf{T}= \begin{pmatrix}\boldsymbol{\ell}_{\mathbf{c}_1} & \ldots & \boldsymbol{\ell}_{\mathbf{c}_l}\end{pmatrix} \in \mathbb{R}^{l\times l}$, $\mathsf{P} = \zeta(\mathbf{p}) \in \mathbb{R}^{l \times 2}$ and  $\mathsf{A}^\mathsf{T}= \left( \mathbf{a}_1 \; \ldots \; \mathbf{a}_{mn} \right) \in \mathbb{R}^{2 \times l}$. The solution of the linear system of equation (A.41) can be written $\mathsf{P} = \mathsf{E} \mathsf{A}$ where $\mathsf{E} = \mathsf{M}^{-1}$. The existence of the matrix  $\mathsf{E}$ is guaranteed if the Schoenberg-Whitney conditions are satisfied (see (57)) as it is the case when the centers are located on a regular grid. Note that the matrix  $\mathsf{E}$ can be pre-computed. Finally, the Feature-Driven parameterization of the FFD warp, denoted  $\mathcal{W}_{\scriptscriptstyle \mathsf{F}}$, is given by replacing the natural parameters  $\mathbf{p}$ in equation (A.38) with their expression in function of the features  $\mathbf{a} = \boldsymbol{\nu}(\mathsf{A}) \in \mathbb{R}^{2l}$:

$$\mathcal{W}_{\scriptscriptstyle \mathsf{F}}(\mathbf{q} ; \mathbf{a}) = \boldsymbol{\ell}_{\mathbf{q}}^\mathsf{T}\mathsf{E} \mathsf{A}.$$ (A.42)

Jacobian matrix of the warp.

The Jacobian matrix  $\mathsf{K}_{\scriptscriptstyle \mathsf{F}}\in \mathbb{R}^{2\times 2l}$ for FFD warps can be computed following exactly the same reasoning as for the TPS warp:

$$\mathsf{K}_{\scriptscriptstyle \mathsf{F}} = \begin{pmatrix} \frac{\partial \mathcal{W}_{\scriptscriptstyl... ...athsf{F}}^y}{\partial \mathbf{a}^y}(\mathbf{q} ; \mathbf{a}) \end{pmatrix}$$

$$= \begin{pmatrix} \boldsymbol{\ell}_{\mathbf{q}}^\mathsf{T}\mat... ...\times l} & \boldsymbol{\ell}_{\mathbf{q}}^\mathsf{T}\mathsf{E} \end{pmatrix}.$$ (A.43)

Extrapolation

The computations involved in the warp reversion operation (see section A.3.3) can lead to evaluate a warp outside of its natural definition domain. More precisely, in equation (A.11), nothing ensures that the features of the vector  $\mathbf{v}$ lies in the domain of the warp. While this is not a problem with the TPS warp whose domain is infinite, extra work need to be done with the FFD warp. Indeed, with the previous definition, it is possible to evaluate an FFD warp outside of its natural domain but it is meaningless since it collapses to 0. In this section, we propose a new method to extrapolate an FFD warp outside of its domain making it virtually infinite.

The principle of the method is simple: a linear extension is added to the basis that crosses the boundaries of the domain (with some extra conditions of continuity and differentiability). While this seems almost trivial in the monodimensional case, it is less simple in two dimensions, i.e. for warps. Our strategy consists in defining the extension in 1D and, then, propagate it to the 2D case using the usual tensor-product.

We present the extrapolation approach in the monodimensional case and for the leftmost boundary of the domain (i.e. , the knot $k_4$). The four non-zero bases that cross this boundary are $N_1$, $N_2$, $N_3$ and $N_4$. Our idea is to drop the part of these bases that are outside the domain and to replace them with a linear extension. We call $N_1^e$, $N_2^e$, $N_3^e$ and $N_4^e$ the bases resulting from this process. In addition to be linear, we enforce the following constraints in order to preserve continuity and differentiability:

$$\left. \begin{array}{l} N_i^e(k_4) = N_i(k_4) \displaystyl... ...rtial N_i}{\partial x}(k_4) \end{array} \right\} \forall i \in {1,\ldots,4}.$$ (A.44)

For the sake of simplicity and without loss of generality, we consider that the leftmost boundary coincides with zero ($k_4=0$) and that the length of the knot intervals is consistently one ($s=1$). Under all these constraints, it follows that:

$$N_1^e(x) = \begin{cases} \displaystyle -\frac{x}{2}+\frac{1}{6} & \textrm{if $x \in (-\infty, 0]$} N_1(x) & \textrm{otherwise} \end{cases}$$

$$N_2^e(x) = \begin{cases} \displaystyle \frac{2}{3} & \textrm{if $x \in (-\infty, 0]$} N_2(x) & \textrm{otherwise} \end{cases}$$

$$N_3^e(x) = \begin{cases} \displaystyle \frac{x}{2}+\frac{1}{6} & \textrm{if $x \in (-\infty, 0]$} N_3(x) & \textrm{otherwise} \end{cases}$$

$$N_4^e(x) = \begin{cases} \displaystyle 0 & \textrm{if $x \in (-\infty, 0]$} N_4(x) & \textrm{otherwise} \end{cases}$$

The extended basis for the rightmost boundary are obtained by symmetry. Figure A.9 illustrates the resulting extended basis functions. The twodimensional counterparts of these newly defined extended basis functions are obtained using the tensor-product.

Figure A.9: (a) Standard basis functions. (b) Extended basis functions that allow one to extrapolate outside the natural domain  $[k_4,k_{m+1}]$.

(a)		(b)

The proposed extension gives a remarkably good behavior to the extrapolating functions. See figure A.10 and figure A.11 for an illustration in 1D and 2D respectively. Besides, the fact that the basis functions form a partition of unity remains true ( $\sum_{i=1}^4 N_i^e(x) = 1, \forall x \in (-\infty,0]$).

Figure A.10: Examples of our extrapolating FFD in the monodimensional case. The extrapolating parts are represented with dashed lines.

Figure A.11: Examples of our extrapolating FFD warp. The dark part of the meshes represents the warp over its initial domain while the light part is extrapolated.