Consistent Landmark Thin-Plate Spline Registration

The averaged unidirectional landmark-based thin-plate spline (AUL-TPS) image registration algorithm produces consistent correspondence only at the landmark locations. The consistent landmark-based, thin-plate spline (CL-TPS) image registration algorithm is designed to align the landmark points and minimize the consistency errors across the entire image domain.

The CL-TPS algorithm is solved by minimizing the cost function given by

$$C = \rho \int_{\Omega} \vert\vert{\mathcal L} u(x)\vert\vert^2 + \vert\vert{\mathcal L} w(x)\vert\vert^2 dx$$

$$+ \chi \int_{\Omega} \vert\vert u(x)-\tilde{w}(x)\vert\vert^2 + \vert\vert w(x)-\tilde{u}(x)\vert\vert^2 dx$$

$$p_i +u(p_i) = q_i \text { and } q_i+w(q_i) = p_i \text{ for } i=1, \dots, M .$$ (4)

The first integral of the cost function defines the bending energy of the thin-plate spline for the displacement fields $u$ and $w$ associated with the forward and reverse transformations, respectively. This term penalizes large derivatives of the displacement fields and provides the smooth interpolation away from the landmarks. The second integral is called the inverse consistency constraint (ICC) and is minimized when the forward and reverse transformations are inverses of one another. This integral couples the estimation of the forward and reverse transformations together and penalizes transformations that are not inverses of one another. The constants $\rho$ and $\chi$ define the relative importance of the bending energy minimization and the inverse consistency terms of the cost function. Notice that this problem is a nonlinear minimization problem since the inverse consistency constraint is a function of the inverse-forward $h^{-1}(x) = x+\tilde{u}(x)$ and inverse-reverse $g^{-1}(x) = x+\tilde{w}(x)$ transformations.

Equation 4 is minimized numerically using the CL-TPS algorithm described in Figure 3. The algorithm is initialized with the forward and reverse displacement fields $u$ and $w$ either set to zero as in Figure 3 or with the result of a previous registration algorithm. The temporary variables $r_l$ and $s_l$ are initially set equal to the landmark locations $q_l$ and $p_l$, respectively, for $l=1, \ldots M$. The value of $r_l$ converges from $q_l$ to $p_l$ as the algorithm converges, and in similar fashion, the value of $s_l$ converges from $p_l$ to $q_l$.

Figure 3: The CL-TPS algorithm registers two images by matching corresponding landmarks in the images while minimizing the inverse consistency error between the forward and reverse transformations.

Consistent Landmark Thin-plate Spline (CL-TPS) Registration Algorithm

1. Initialization: Set $u(x)=0$, $w(x)=0$, $r_l=q_l$, and $s_l=p_l$.
2. Compute $f_1(x)$ that satisfies $\nabla^4 f_1(x) = 0$ subject to $f_1(r_l)=p_l-r_l \quad \forall l$ using the periodic boundary UL-TPS algorithm.
3. Compute $f_2(x)$ that satisfies $\nabla^4 f_2(x) = 0$ subject to $f_2(s_l)=q_l-s_l \quad \forall l$ using the periodic boundary UL-TPS algorithm.
4. Set $u(x)=u(x)+\alpha f_1(x)$ and $w(x)=w(x)+\alpha f_2(x)$.
5. Set $r_l=q_l+u(r_l)$ and $s_l=p_l+w(s_l)$.
6. Compute $h^{-1}$ and $g^{-1}$ using procedure described in [34].
7. Set $\tilde{u}(x)= h^{-1}(x)-x$ and $\tilde{w}(x)=g^{-1}(x)-x$.
8. Set $u(x)=u(x)+\beta [u(x) - \tilde{w}(x)]$ and $w(x)=w(x)+\beta [w(x) - \tilde{u}(x)]$.
9. If the maximum landmark error $\vert u(q_l)-(p_l-q_l)\vert$ or $\vert w(p_l)-(q_l-p_l)\vert$ is greater than a threshold $\epsilon_1$ or the maximum inverse consistency error $\vert u(x) - \tilde{w}(x)\vert$ or $\vert w(x) - \tilde{u}(x)\vert$ is greater than a threshold $\epsilon_2$ then Goto 2.

At each iteration of the algorithm, the unidirectional landmark thin-plate spline (UL-TPS) algorithm with periodic boundary conditions is used to solve for the perturbation field $f_1$ that minimizes the distance between the current position of $r_l$ and its final position $p_l$. The perturbation field $f_1$ times the step size $\alpha$ is added to the current estimate of the forward displacement field $u$ where $\alpha$ is a positive number less than one. This procedure is repeated to update the reverse displacement field $w$. Next, the forward displacement field $u$ is updated with the step size $\beta$ times the gradient of the inverse consistency constraint with respect to $u$ assuming that $\tilde{w}$ is constant. The displacement field $\tilde{w}$ is computed by taking the inverse of the transformation $g(x)=x+w(x)$ as described in our previous paper describing the consistent intensity registration algorithm [34]. This step is repeated in the reverse direction to update the displacement field $w(x)$. These steps are repeated until the landmark error and the inverse consistency error fall below problem specific thresholds or until a specified number of iterations are reached. In practice, this algorithm converges to an acceptable solution within five to ten iterations and therefore we use a maximum number of iterations as our stopping criteria.

Copyright © 2002 • The University of Iowa. All rights reserved. Iowa City, Iowa 52242
Questions or Comments: gary-christensen@uiowa.edu