Rigid Shape Registration Based on Extended Hamiltonian Learning

Shape registration, finding the correct alignment of two sets of data, plays a significant role in computer vision such as objection recognition and image analysis. The iterative closest point (ICP) algorithm is one of well known and widely used algorithms in this area. The main purpose of this paper is to incorporate ICP with the fast convergent extended Hamiltonian learning (EHL), so called EHL-ICP algorithm, to perform planar and spatial rigid shape registration. By treating the registration error as the potential for the extended Hamiltonian system, the rigid shape registration is modelled as an optimization problem on the special Euclidean group SE(n) (n=2,3). Our method is robust to initial values and parameters. Compared with some state-of-art methods, our approach shows better efficiency and accuracy by simulation experiments.


Introduction
Point set registration is to find an optimal transformation that aligns one set to the other. This problem occurs in many applications, such as motion tracking [1,2], target localization [3], super-resolution [4], mosaicing [5] and medical image analysis [6][7][8][9][10]. Depending on different goals, the methods can be categorized as coarse registration and fine registration [11,12]. The former is to find an initial estimation between two point sets while the latter is to obtain a more accurate solution. In general, the popular method used in coarse registration is point-to-point and iterative, such as point signature and spin image. Some other methods may include principle component analysis (PCA), algebraic surface model and principle curvature. While methods like iterative closed point (ICP), Chen's method, signed distance fields and genetic algorithms, are more commonly used in fine registration.
Literately, the ICP algorithm is efficient and accurate among all point set registration techniques. However, it requires that the initial configuration of the two point sets is sufficiently close and it is sensitive to noise and outliers. Moreover, the ICP algorithm is unable to handle affine case. To improve the robustness and cope with more situations, many researchers extended the ICP method. In particular, Ying et al. applied the theory of differential geometry and Lie group to ICP algorithm and formed a unified framework to design fast and robust shape registration algorithms [13][14][15][16][17][18][19][20][21].
From the geometric point of view, given the correspondence between two sets of points, the affine registration is to find an element in the general linear group that transfers one set to the other optimally. In this paper, we deal with the rigid registration, which uses the geometry of the special Euclidean group. From this perspective, differential geometry is a powerful tool not only in registration problems but also in computer vision extensively. There is a vast community dealing with shape registration/analysis/comparison from geometrical viewpoints. For example, E. Celledoni etc. applied the theory of Lie groups and homogeneous manifolds to the problem of shape analysis [22,23]. A.L. Brigant introduced an algorithm that finds an optimal matching between two curves by computing the geodesic of the infinite-dimensional manifold of curves [24]. X. Pennec etc. concentrated on feature-based approaches for rigid registrations using differential geometry of surfaces [25]. Other methods, using harmonic analysis and statistical optimizations, can be found in [26][27][28][29]. Experiments showed that these geometry-based approaches performed better than other state-of-art methods.
There is also an immense literature in applying Hamiltonian dynamical systems to learning theory. F. Barbaresco studied the symplectic extension of Souriau Lie groups thermodynamics and used this model for data analysis and machine learning on Lie groups [30,31]. S. Fiori proposed extended Hamiltonian learning (EHL) on Riemannian manifolds motivated by Hamilton variational principle [32,33]. Compared with the classical gradient-based learning, EHL has the advantages of averaging out the oscillations and mitigating the plateau effect. In this paper, we investigate the 2D/3D rigid registration problem from the view of Hamiltonian learning. The innovation of our proposed method is treating the registration error as the potential of the extended Hamiltonian system on the underlying space of transformation, the special Euclidean group. Under the system of extended Hamiltonian, an EHL-based algorithm is designed to achieve the optimal transformation corresponding to the minimum of the registration error. Some experiments are carried out to validate the efficiency and robustness of our algorithm.
The structure of the paper is organized as follows. In Section 2, some basic background about geometric structures of the special Euclidean group SE(2) and SE(3) is reviewed, including the Riemannian metric, geodesic, and Riemannian gradient. In Section 3, extended Hamiltonian learning on a general Riemannian manifold is discussed. And in the following section, we formulate the method for 2D/3D rigid registration problem by extended Hamiltonian learning on SE(2) and SE(3). In Section 5, some numerical experiments prove the efficiency of this proposed method, comparing with Du's SVD method (SVD) [16], Ying's Iwasawa decomposition method (ID) [21], Ying's Lie group optimization method (LGO) [18]. At last, conclusions are presented and possible improvements of our method for future work are discussed.

Geometry of Special Euclidean Groups
In this section, we review some basics about the special Euclidean group SE(n), which is the underlying geometric space for rigid registration. The special Euclidean group consists of rotations and translations in the Euclidean space. It is the semi-direct product of the special orthogonal group SO(n) and R n , Being a Lie group, SE(n) is equipped with the smooth structure and group structure simultaneously. Any element g in SE(n) can be represented as a pair (r, t), where r ∈ SO(n) and t ∈ R n stands for the rotation and translation respectively. In the matrix form, g can be written as which is an one-to-one correspondence with the pair (r, t). The action of the special group SE(n) on the Euclidean space is defined as The Lie algebra of SE(n), denoted by se(n), contains infinitesimal rotations and translations. A general element h in se(n) is a pair (J, v), where J is a n × n skew-symmetric matrix and v is a vector in R n . In matrix form, h can be rewritten as The exponential mapping and the logarithm mapping provide methods to transfer information between the nonlinear manifold SE(n) and the linear space se(n). Concretely, the exponential map of a general Lie algebra element is defined by In low dimensions the exponential map can be written in a compact form, as presented by J.M. Selig [34]. Denote ω = ω 2 For any h ∈ se (3), direct computation shows that it satisfies the quartic equation h 4 + ω 2 h 2 = 0. Thus the higher order terms in the taylor expansion (5) can be simplified. Explicitly, we have Especially, for pure rotations, where v = 0 in the expression (4) and h satisfies the cubic equation h 3 + ω 2 h = 0, the exponential map can be written as while for pure translations, where J = 0 in the expression (4) and h 2 = 0, the exponential map satisfies The exponential map for pure rotation degenerates on se(2) to be as the following, An inner product on the Lie algebra g can be extended to a Riemannian metric on the Lie group G by left translation. Specifically, the inner product at se(n) is defined by where m > 0 is a fixed parameter, (J 1 , v 1 ), (J 2 , v 2 ) ∈ se(n). The left-invariant metric, for two tangent vectors α 1 and α 2 at an arbitrary group element g ∈ SE(n) is defined as The right-invariant metric is defined similarly. In physics, the left-invariance and right-invariance correspond to the independence of the choices of the inertial frame and the body-fixed frame respectively [35]. A bi-invariant metric means it is both left-invariant and right-invariant. However, there is no bi-invariant metric on SE(n) [36]. Thus, we adopt the left-invariant metric on SE(n) throughout this paper.
Let f be a function defined on Riemannian manifold M, ∂ x f and ∇ x f be Euclidean gradient and Riemannian gradient, respectively. The relation between ∂ x f and ∇ x f is governed by [32] where ·, · E and ·, · R denote Euclidean inner product and Riemannian metric respectively. For the Riemannian metric defined in (12) on SE(n), we have

Extended Hamiltonian Learning on SE(n)
We first introduce the extended Hamilton learning on general Riemannian manifolds [32,33]. Then we specify on special Euclidean group SE(n).
On a Riemannian manifold M, equipped with a metric ·, · x at T x M, the extended Hamiltonian principle is In Equation (15), x(t) ∈ M stands for the trajectory of a particle moving along M andẋ(t) is the corresponding instantaneous velocity. The function K : T x M × T x M → R denotes the kinetic energy and V represents associated potential energy. The symbol δ represents the variation of the action in the dynamical system, namely, the particle slides from a point to an infinitely close point. An element f x ∈ T x M indicates the dissipation force at the point x. This system degenerates to conservative if the dissipation force disappears everywhere. Following [32,33], the kinetic energy is adopted as the symmetric bilinear form K x (ẋ,ẋ) = 1 2 ẋ,ẋ x under the assumption of unit mass. From viscosity theory the dissipation term is assumed to be f x = −µẋ with µ ≥ 0 in the paper. Then, the equations for such dynamical system read where Γ x is the Christoffel symmetric form, and ∇ x V denotes the Riemannian gradient of V.
In order to implement simulations, we turn the continuous version (16) into a discrete one, which yields where R x k : T x k M → M is the exponential map and η > 0 is the selected step size.
When it comes to the special Euclidean group SE(n), we need to compute the Christoffel symmetric form. However, in the discrete case, it suffices to know the form in a subspace so(n) of se(n). Specifically, the Christoffel symmetric form for g = (r, t) ∈ SE(n), h = (J, 0) ∈ se(n) reads Γ g (gh, gh) = (−rJ 2 , 0).
Then Equation (17) can be rewritten as Note that from the iteration (19) one cannot ensure that J k+1 is skew-symmetric. Consequently r k+1 will not be an element in SE(n). Thus, J k+1 should be modified to keep the validity of the iteration. To do this we make an orthogonal projection from the Euclidean space R n 2 to the subspace consisting of skew-symmetric matrices. i.e., For any matrix J k+1 , it can be decomposed as the sum of a skew-symmetric part and a symmetric part (19), the symmetric part 1 2 (J k+1 − J T k+1 ) would be a negligible error when the step is small enough. Thus we apply the skew-symmetric part in the tangent space T r k+1 SE(n). See Figure 1. We remark that this step is nothing but the to compute the covariate derivative on Euclidean submanifolds by definition [37]. Hence, the iteration becomes For details about the convergence of extended Hamiltonian learning, refer to [32]. The necessary condition for convergence is that µ satisfies √ 2λ max < µ < η −1 , where λ max is the maximum eigenvalue of the Hessian matrix of V.

2D/3D Rigid Shape Registration Based on Extended Hamiltonian Learning
Given two n-D (n = 2,3) data sets X = {x i } N x i=1 and Y = {y j } Ny j=1 and a correspondence π : {1, · · · , N x } → {1, · · · , N y }, denoting z i = y π(i) , find an element g = (r, t) ∈ SE(n) such that the cost function achieves it minimum. The fundamental steps based on ICP for rigid registration are: First, for the current fixed transformation g k = (r k , t k ) ∈ SE(n), find Z k ⊆ Y with N x elements such that the subset minimizes Second, once the correspondent data set Z k is obtained, update the transformation g k as g k+1 for minimizing In fact, the translation t can eliminated by coordinating the centers of the two data sets. Let x c and y c be centers of X and Y, respectively. For the obtained set Z k , the least squares solution (r k+1 , t k+1 ) of (23) indicates t k+1 = y c − r k+1 x c .
Thus, with the centralized datax i = x i − x c , optimization problems (22) and (23) can be simplified as Here, we regard the registration error as potential of extended Hamiltonian system on the special Euclidean group SE(n). The Euclidean gradient of V(r) is given as from which we can compute the Riemannian gradient of V(r) by (14).
Therefore, for two given data sets , we summarize the method based on extended Hamiltonian learning as Algorithm 1.

Numerical Results
All data samples used in this section is from the MPEG-7 shape B database and all programs are written in Matlab 2018a and run by PC with AMD Athlon II P340 Dual-Core processor, 2.20 GHz CPU, and 2 GB RAM.

2D Rigid Shape Registration
In 2D case, our method appears to be robust and insensitive to initial values. Moreover, a near-optimal registration can be obtained within a few steps. To give a visualization of our method, we test chicken-2 and chicken-3 in MPEG-7 shape B database as model data set and test data set. The initial rotation is set to be the identity and the initial translation is set to be the difference of two means of data. The numerical results are shown in Figure 2. We select nine groups of rigid data in MPEG-7 shape B database to run the experiments and compare with Du's SVD method, Ying's Iwasawa decomposition method and Ying's Lie group optimization method. To give a quantitive comparison we define the root mean square (RMS) error to be Though the geodesic distance on SE(2) is more reasonable in theory, practically it is difficult to compute since we do not know the true rigid transformation.
The precision is set to be = 0.00001. Comparison results are displayed in Figure 3. The resulted RMS errors are displayed in Table 1. Other methods require the careful choice of initial values, whereas for our method we simply choose identity as the initial rotation and difference of means as initial translation. We can find that EHL-ICP algorithm is more robust to the size and shape of point cloud data. Other algorithms may have good performance on small point sets with simple shapes but lose precision when point sets are complex.

3D Shape Registration
Similar to the 2D case, we selected a group of 3D models including bunny, chair, cactus, dinosaur, elephant and block to verify the validity of our algorithm. The initial rotation is set to be identity and the initial translation is set to be the difference of means of model data and test data. The visualized results are represented in Figure 4. Note that we do not require a subtle choice of initial values and parameters. The results demonstrate the efficacy of our method.

Conclusions and Future Works
Shape registration plays a significant role in computer vision, where the task is to transfer one set of points to the other. Since the iterative closest point method is widely used in registration problems yet having some shortcomings, this paper proposes the EHL-ICP method, which incorporating the extended Hamiltonian learning with the ICP algorithm, to deal with the 2D and 3D rigid shape registration problem. By regarding the registration error as the potential of the extended Hamiltonian system, we formulate rigid registration as an optimization problem on the special Euclidean group SE(n) (n = 2, 3). Numerical results show that our method is more effective and accurate when compared with other methods. Moreover, our method is robust with respect to initial values in both dimensions, which provides a good choice for rough registration.
For future work, we may generalize the extended Hamiltonian learning method to different registration problems. There are two hot topics worth mentioning. The first is to use optimal transportation for data set registration [38,39]. The basic idea is that shape data can be viewed as a sum of Dirac measures in a given space and the difference of shapes is taken to be the Wasserstein distance. To find a (non)rigid transformation is to find a measure-preserving map. Another possible application is affine registration where shapes are distorted and there is no rigid transformation [40,41]. From this viewpoint, we should consider extended Hamiltonian learning on the general linear group GL(n), where more techniques should be developed.