TPS-RPM review
We first review the mathematical framework of TPS-RPM proposed by Chui et al. [4]. Given the source point-set X = {x
i
, i = 1, 2, . . . , K} and the target point-set Y = {y
j
, j = 1, 2, . . . , N} in a region Ω, the goal of TPS-RPM is to find the optimal transformation h : Ω → Ω that maps the source point-set X to the target point-set Y , as well as estimating the corresponding relationship between X and Y . In TPS-RPM, TPS is employed to model the transformation with parameters (a, W), which maps points in X as,
(1)
where a and W are affine transform matrix and warp coefficient matrix respectively, w
j
is an element of matrix W, r
ij
= ||x
i
− x
j
|| is the distance norm between point x
i
and x
j
, ϕ(r
ij
) is the basis function of TPS.
A fuzzy correspondence matrix M with dimension (K + 1) × (N + 1) is defined to describe the correspondence between points. Since the one-to-one correspondence relationship between point sets X and Y will probably not always exist, outlier point is defined as corresponding point of the isolated point. Therefore, each row and each column of matrix M has an extra outlier point. The fuzzy correspondence of point x
i
and y
j
is defined as follows:
(2)
where T is the temperature in the anneal procedure of TPS-RPM. The fuzzy correspondence matrix is subject to for j ∈ {1, 2, . . . , N}, for i ∈ {1, 2, . . . , K}, and m
ij
∈ [0, 1]. The nearer the distance between the mapped x
i
and y
j
is, the more likely a corresponding relationship exists between x
i
and y
j
.
RPM employed soft assign and deterministic annealing technique to estimate the fuzzy correspondence matrix M and the transformation h simultaneously that minimize the following cost function:
(3)
subject to for j ∈ {1, 2, . . . , N}, for i ∈ {1, 2, . . . , K}, and m
ij
∈ [0, 1]. Here, h is the transformation maps two point-sets with components h
d
, d = 1, . . . , D, where D is the dimension of Ω. L is the linear elastic operator, and ||Lh||2 is used to measure the smoothness of the transformation h, i.e.
(4)
The cost function is derived from a statistical physics model. The term is a barrier function, which is used to push the minimum of the cost function away from the discrete points. The temperature T contorls the degree of convexity of the cost function [3]. When T is sufficiently small, the cost function is ensured to be convex. λ and ζ are regularization parameters. In the TPS-RPM algorithm, Expectation-maximization (EM) algorithm is adopted to solve M and h iteratively, the detailed process can be found in literature [4].
When TPS-RPM is used to register the source image I and target image J , the source point set X and the target point set Y are extracted from I and J respectively. Next, TPS-RPM is employed to estimate the forward transformation h : X → Y , which is the transformation to map the source image I to the target image J so that I(h(x)) = J . When image J is registered to image I, the backward transformation g : Y → X maps the image J to image I so that J (g(x)) = I. As previously mentioned, it is required that the forward transformation and the backward transformation are inversely consistent, i.e. g ○ h = id and h ○ g = id, where id is the identity map, to ensure the correspondence between the two images to be consistent. However, the forward transformation h and the backward transformation g is not dependent to each other for TPS-RPM, since TPS is an unidirectional function which results in a non-consistent correspondence between the two images except at the control points, that is, g ○ h ≠ id, h ○ g ≠ id and g ○ h ≠ h ○ g. Furthermore, the value of the fuzzy correspondence matrix M is computed based on the mapping errors in the forward transformation only, the mapping errors from Y to X will not be penalized, which leads to a bias matching result.
Inverse consistent robust point matching
Firstly, we introduce several notations used in this paper. The forward transformation from the source point set X to the target point set Y is denoted as h(x), the displacement field is u(x) = h(x) − x; the backward transformation from Y to X is denoted as g(x), the displacement field is w(x) = g(x) − x. The inverse of the forward transformation is h− 1, the corresponding displacement field is , and the inverse of the backward transformation is g− 1, the corresponding displacement field is . An inverse consistent registration is required to satisfy g ○ h = id and h ○ g = id. In other words, the forward and backward transformations estimated by an inverse consistent registration should satisfy g = h− 1and h = g− 1in the region Ω. Johnson et al. [24] defined the inverse consistency constraint as ||h − g− 1||2 + ||g − h− 1||2, which makes sure that the function of forward transformation is similar with the inverse function of backward transformation. Correspondingly, the function of backward transformation is as similar with the inverse function of forward transformation as possible. We impose the inverse consistency constraint on the RPM optimization problem by minimizing the cost function given by
(5)
In (5), the mapping errors between two point sets are extended as the combination of distance between the target point and the mapped position of the source point using the forward transformation, and the distance between the source point and the mapped position of the target point using the backward transformation, instead of only using the forward mapping errors. Both the smoothness of the forward and backward transformations ||Lh||2 + ||Lg||2 are included in the cost function. Χ is the weighting parameters to make a trade-off between the inverse consistent error and other terms.
The goal of the inverse consistent robust point matching is to estimate the inversely consistent forward and backward transformations for two sets of points concurrently, as well as making clear the correspondence between X and Y bidirectionally. The correspondence matrix in traditional RPM is only based on the unidirectional transformation between the target point set and the source point set, while the value of correspondence m
ij
for two points in our algorithm is inversely proportional to the mapping errors of points bi-directionally, i.e.
(6)
Furthermore, to register images, the similarity of local image is introduced to the correspondence as,
(7)
where, I(x
i
) and J (y
j
) are two local regions centred at x
i
in image I and y
j
in image J . I(h(x)) and J (g(x)) are deformed images of I and J using the forward transformation and the backward transformation respectively. corr is the correlation coefficient used to measure the similarity between two local regions. T
s
is the temperature parameters of image similarity. By introducing image information to the fuzzy correspondence matrix, improvement of image matching is achieved for the inverse consistent RPM.
To find M, h and g that used to optimize formula (5), we still use the iterative strategy proposed in [4]. The iterative process includes the E step to calculate the fuzzy correspondence matrix according to the current estimated forward and backward transformations. Next, it performs the M step to estimate the forward and backward transformations on the basis of the current estimated fuzzy correspondence matrix. By dropping the terms independent of h and g, it is needed to minimize the following objective function:
(8)
where
, i = 1, 2, . . . , K and , i = 1, 2, . . . , N are the virtual points computed in the forward and backward directions respectively. Moreover, v
i
is expected to be corresponding to x
i
, and z
j
is expected to be corresponding to y
j
also. v
i
and z
j
are held fixed during the procedure of the M step. Then, the optimization problem is to find the optimal forward and backward transformations h and g given four point sets {x
i
}, {y
j
}, {v
i
} and {z
j
}, where {x
i
} are corresponding to {v
i
}, and {y
j
} is corresponding to {z
j
}. The iterative process continuously alternates the E step with the M step until it converges. Next, we will discuss how to calculate transformations h and g at the same time by optimizing Ec(h, g).
Modified consistent landmark thin-plate spline registration
Given two point sets with known correspondence relationship, the Consistent Landmark Thin-Plate Spline (CLTPS) registration algorithm [24] was originally proposed to solve the inversely consistent transformations between these two point sets. During the procedure of CLTPS, only two point sets are used to estimate the forward and backward transformations simultaneously. However, there are four point sets {x
i
}, {y
j
}, {v
i
} and {z
j
} in RPM. Based on the correspondence between {x
i
} and {v
i
}, and the correspondence between {y
j
} and {z
j
}, an intuitive approach to estimate the forward transformation h is to let {x
i
} be the source point set and {v
i
} be the target point set. Conversely, let {y
j
} be the source point set and {z
j
} be the target point set to estimate the backward transformation g. Details of CLTPS can be referred in [24].
However, there several existed problems in CLTPS: (1) the mapped positions of control points are oscillated near their target positions, instead of mapping exactly to the target positions [31]; (2) topology of the forward and backward transformations cannot be ensured to be preserved.
Firstly, there is a minor oscillation problem in CLTPS algorithm. In CLTPS, the forward and backward displacements are updated iteratively using the temporary forward and temporary backward transformations f
1(x) and f
2(x), where f
1(x) is estimated by considering the current mapped position of {x
i
} and {v
i
} as the source and target control point sets respectively, and f
2(x) is estimated by considering the current mapped position of {y
j
} and {z
j
} as the source and target control point sets respectively. However, in CLPTS, x
i
can be mapped to a location near to v
i
, but cannot be mapped to v
i
exactly. The same goes for y
j
also. To tackle the oscillation problem of CLTPS, we propose a new approach to update the forward and backward displacements iteratively. Denote r
i
and s
j
as the temporary mapped positions of x
i
and y
j
respectively. After the k th iteration, x
i
is mapped to r
i
using the current forward displacement u
k
(x), and y
j
is mapped to s
j
using the current backward displacement w
k
(x). We update the forward and backward displacements iteratively as follows:
(9)
We use the forward displacement to demonstrate the improvement of the update. α = 1 is assumed so as to simplify the analysis, then, at the k + 1th iteration, the displacement of x
i
is,
(10)
It implies that x
i
is mapped to v
i
exactly using our approach. Similarly, we can prove that y
j
is mapped exactly to z
j
using the backward displacement.
Secondly, the forward and backward transformations estimated by CLTPS cannot be ensured to be topology-preserving, since the temporary transformations f
1(x) and f
2(x) are estimated by TPS, which does not enforce one-to-one mapping. Topology preservation is an important property of a deformation, which ensures that connected structures remain connected, and that the neighborhood relationships between structures are maintained before and after warping [33]. In image registration, topology preservation of deformation fields can prevent disappearing of existing structures or introducing new artificial structures after image warping. However, transformations estimated by TPS are not constrained to be topology-preserving as they are motivated by small deformation kinematics [34], and they do not allow for large deformations that maintain the topology of the template [35]. To preserve the topology of the deformation field, the regularized TPS proposed by Chui et al. [4] is employed to estimate the temporary forward transformation f
1(x) and the temporary backward transformation f
2(x), which preserves topology of deformation fields better than TPS. As shown in Figure 1, the source points (circle points) are expected to be mapped to the target points (star points). The regularized TPS produces a smooth and topology-preserving deformation field, while TPS makes the deformation field folding, which is non-topology-preserving. The parameter used in the regularization procedure is decreased gradually to preserve the correspondence between points. Moreover, h and g are required to be topology-preserving, so, after each adjustment of u(x) and w(x), the Jacobian values of h and g are computed, when either of the minimum Jacobian values of h and g is negative, the iteration is stop.
Finally, r
i
and s
j
are required to be updated as the newest mapped positions of x
i
and y
j
for each iteration. So, after the update of the forward and backward transformations, r
i
and s
j
are updated correspondingly using the latest transformations respectively. More importantly, r
i
is closer and closer to v
i
with the increase in the number of iterations, rather than swinging nearby v
i
as CLTPS. Similarly, s
j
is closer and closer to z
j
in the iteration process. All these ensures that x
i
is mapped exactly to its target position v
i
, and y
j
is mapped exactly to its target position z
j
using the modified consistent landmark thin-plate spline registration algorithm.
Details of the modified consistent landmark thin-plate spline registration are described in algorithm 1.
Algorithm 1 Modified Consistent Landmark Thin-Plate Spline (CLTPS) registration algorithm using four points sets.
1: Let r
i
= x
i
, s
j
= y
j
; u(x) = 0, w(x) = 0, the steps α and β, the mapping error threshold ξ of control point, and the maximum number of iteration m
iter
, k = 1.
2: Regularized TPS is performed to estimate the temporary forward transformation f
1(x) based on the correspondence between r
i
and v
i
, and the temporary backward transformation f
2(x) based on the correspondence between s
j
and z
j
.
3: Update the displacements,
u(x) = u(x) + αu∗(x), u∗(x) = u
t
(u(x) + x), u
t
(x) = f
1(x) − x,
w(x) = w(x) + αw∗(x), w∗(x) = w
t
(w(x) + x), w
t
(x) = f
2(x) − x.
4: Compute the Jacobian values of the forward and backward transformations, if the minimum Jacobian values of the forward or backward transformations are negative, the iteration is terminated.
5: Get h− 1(x), the inverse function of forward transformation and g− 1(x), the inverse function of backward transformation.
6: Update displacement field of forward and backward transformation. u(x) = u(x) − β[u(x) − g− 1(x) + x], meanwhile, w(x) = w(x) − β[w(x) − h− 1(x) + x].
7: r
i
and s
j
are updated as r
i
= x
i
+ u(x
i
), s
j
= y
j
+ w(y
j
).
8: Check whether the termination condition is met. If k > m
iter
or |u(x
i
) − (v
i
− x
i
)| < ξ or |w(y
j
) − (z
j
− y
j
)| < ξ, the iteration is terminated; otherwise, k = k + 1, go to step 2.