Nonrigid Points Alignment with Soft-weighted Selection

1 Introduction

As a crucial technique of computer vision, point set registration, is the process of recovering the correspondences and determining the spatial transformation between the point sets that contain the same context sampled at different time, from different viewpoints or by different sensors. In most applications, the matching results greatly influence the accuracy of the subsequent processing tasks like dynamic target detection, simultaneous localization and mapping (SLAM)[Hong and Kim, 2016] and object reconstruction [Zhou et al., 2015]. Therefore, it is necessary to develop a robust method that can automatically align point sets with sub-pixel accuracy.

Recently, lots of PSR methods have been proposed, which mainly include two categories. One uses rigid transformation to model the spatial deformation between two input point sets, which is simple but has a lower accuracy for complex distortions, and the other models the spatial deformation as nonrigid transformation by a set of base functions. Since spatial deformation between point sets has multi-degree of freedom in real scenario, non-rigid registration is closer to real deformation. So, it is more widely used than the rigid ones. Owing to following difficulties, PSR becomes pretty challenging: (1) degradations: the distribution of point set would be more complex in the presence of large degradations, e.g., noise, occlusion, outlier, scale, and multi-view changes. The noisy data often means the points cannot be matched precisely. Moreover, the data with occlusion and outliers could cause some points fail to find the true correspondences. Consequently, the registration accuracy is unstable for these factors; (2) the numerical optimization often falls into local minima; (3) the computational complexity is higher when aligning the point set with great number.

∗Corresponding author: crabwq@nwpu.edu.cn

For these challenges, earlier methods mainly applied distance ratio to recover point correspondences by comparing the distance of the closest neighbor with that of the secondclosest neighbor [Lowe, 2004]. However, some outliers or missing matches always exist in this manner. Particularly, the outliers or missing matches become more serious as the distractions, e.g., noise, similar structure or improper threshold, increase. Note that the outliers have a negative effect on the estimation of spatial transformation while less matches mean more unreliable parameter calculation. To remove outliers, geometric structures [Zhang et al., 2014] and spatial constraints [Kahaki et al., 2016], are often used. Although these methods performwell, it is infeasible to align non-rigid point sets for the unknown underlaying transformation.

Unlike above matching, another often-used approach is the one based on iterative optimization. Such methods first formulate point matching as an optimization problem, and then solve it via numerical algorithms. The typical methods based on Gaussian fields criterion [Myronenko and Song, 2010] used a Gaussian mixture model to maximize the overlap between the features in both images. Whereas this criterion is generally differentiable, some gradient-based optimization algorithms can be directly utilized to solve it. Besides, it does not need explicit point correspondences and is more robust to noise as well. However, most of them are only suitable for rigid deformations. For the nonrigid case, Ma et al. [Ma et al., 2015c] improved it via a regularized constraint and sparse approximation. Wang et al. [Wang et al., 2017] used a laplacian regularized term to preserve the intrinsic geometry of the transformed set. Qu et al. [Qu et al., 2017] solved nonrigid point matching problem by regression and clustering under a Bayesian framework. However, these methods only aim to minimize the global error. Since the position accuracies of different points are usually different, this means they should have different contributions to the accumulation of global error. In this sense, only minimizing the global error could make parameter estimation imprecise.

There are often some points (called confused samples) that should be outliers or noises in the transformed set, and the global error would be smaller when using them than without them as shown in Fig. 1. Such confused points easily mislead the iterative model, even make it struck in local extremum. In this paper, we introduce a self-selected mechanism to remove the confused points by guiding the convergence of the model and propose a robust point matching criterion. Our main contributions are:

• We formulate the point matching problem as a sparse approximation with low rank constraint in RKHS, which makes our model more efficient.
• An adaptive self-selected mechanism is presented to remove the confused points, which adaptively select reliable points as iteration and make our the model focus more on the faithful samples;
• An equivalent matching number constraint is embedded into the non-rigid criterion, which enhances the reliability of point matching.

2 RelatedWork

Point matching has been widely applied to computer vision. Over past few years, many point matching methods have been proposed, which can be divided into two aspects. The first is the ones based on RANSAC [Lowe, 2004]. As a robust estimator, RANSAC tries to find a feasible subset to estimate the transformation by multiple sampling. Based on it, several improvements are presented [Chou and Wang, 2015].For instance, Litman et al. [Litman et al., 2015] proposed an inverting RANSAC to detect inliers in the given set. To get more matches, Wu et al. [Wu et al., 2015] first applied a fast sample consensus (FSC) to obtain initial correct matches, and then tried to find the precise matches from the initial matches as many as possible. Although these methods perform well in most cases, they are designed for rigid matching and unsuitable for nonrigid alignment.

Another popular strategy is to iteratively match the points and solve the optimal spatial transformation, such as coherence point drift (CPD) [Myronenko and Song, 2010], progressive vector field consensus (PVFC) [Ma et al., 2015a], robust point matching with L2E estimator (PRM-L2E) [Ma et al., 2015c], image gradient orientations (IGO) [Zheng et al., 2017]. Besides, Tsin [Tsin and Kanade, 2004] employed the kernel correlation (KC) to align the points. Afterwards, Bing [Jian and Vemuri, 2011] introduced the Gaussian mixture models (GMMs) and expanded the method [Tsin and Kanade, 2004]. To preserve local or global structures better, Zheng and Doermann [Zheng and Doermann, 2006] constructed a simple graph matching initialized with the shape context distance, and proposed a robust point matching by preserving the local structures (RPM-LNS). Ma et al. [Ma et al., 2016] formulated point alignment as the estimation of a mixture of densities. Bai et al. [Bai et al., 2018] introduced local connectivity constraint by k-connected neighbors and connectivity matrix. In these works, PSR is often formulated as the optimal problems based on GMMs or graph matching [Pinheiro et al., 2017], which has a limited calculation efficiency in real applications.

3 Methodology

3.1 ProblemFormulation

Given two point sets, the model set X = {xi}ni=1 and the data set Y = {yr}mr=1 , the goal of registration is to recover the correspondence between {xi}ni=1 and {yr}mr=1 . This process is moving the model set X onto the data set Y by a series of transformation τ. In general, the best spatial transformation τ(·) : Rd ↦→ Rd is just the one that has the maximum pointto-point overlapping between X and Y by τ , where d is the dimension of one point. When we assume that the noise on the inlier correspondence (xi, yr) is Gaussian on each component with zero mean and uniform standard deviation σ, the correspondence (xi, yr) satisfies yr - τ(xi) ∼ N(0, σ2 I), where I is a d×d identity matrix. Thus, the point correspondence problem can be regarded as finding the largest samples (i.e., the inlier points) that obey the normal density model.

Defining displacement function f(x) = τ(x) - x, above matching problem can be further written as follows

where

Φ(f) is a regularization term controlled by constant λ for trading-off the two terms, and the stabilizer Φ(f) is an inner product Φ(f) = ‖f‖2Γ = 〈f , f〉Γ , where ‖ · ‖2Γ is the Hilbertian norm of the vector-valued RKHS, d(x, y) denotes the Euclidean distance between vector x and y, and wir is the weight coefficient.

In Eq. 1, Gaussian distance is used as measure function. And if model point xi and data point yr are matched, the Gaussian distance is equal to 1 for d(τ (xi), yr) = 0, whereas if xi, yr are not matched or outliers, d(τ (xi), yr) is generally larger and the Gaussian distance is closer to 0. So, Eq. 1 aims to find the largest matches. However, due to the disturbances in the data, e.g., noises, local similar structures and minor outliers, the model is more easily confused so that converges to local extremum. To avoid above problems, a self-selected mechanism is introduced to guide model convergence by dynamic sample selection and incoordinate treatment in this paper. Specifically, sample label is assigned for each model point to select suitable sample points first. A self-selected function w.r.t. the label is defined to measure how much contribution each sample point devotes to the model, which can restrain the matching number and work for dynamically updating the label. As a result, we formulate the point alignment as the following self-selected regularized criterion:

where W is the weight matrix, written as W = [w1, · · · ,wm] ∈ Rn×m with wr = [wr1, · · · ,wrn]T ; υi ∈ [0, 1] is the sample label for the model point xi that shows whether it is selected. κ(:, ς) is a function controlled by the parameter ς > 0, which specifies how xi is assigned. ℓ(x) = 1 - exp{-x} and ‖ · ‖ denotes the L2 norm.

3.2 Self-selectedMechanism

Here, label vi is used as a measure of selection or deselection, which can be performed by a hard-value within {0,1}, i.e., setting selection with vi = 1 while deselection with vi = 0. But such hard-value cannot reflect the importance of different samples to the model. In this paper, we treat the selected samples with unequal weights. vi ∈ [0, 1]n is set in a mixture of a hard 0-1 and a soft-assignment manner. More importantly, each vi is dynamically updated in order to select the best sample points for current model during each iteration. For this purpose, we introduce a self-selected function κ(υi; ς) w.r.t. vi, just as in [Jiang et al., 2014], and it is crystallized as follows

where the extra parameter ζ is uesd to control the strength of the weights assigned to the selected points.

Let ℓi = ∑r wirℓ(ρir) be the accumulated error of model point xi. To discuss vi better, omit the items without vi , and then for each model point xi, Eq. 3 can be simplified as

Note that Eq. 5 is continuous and differentiable with repect to vi. Thus, it is easy to derive its optimal solution υ∗ i by taking a derivative with respect to vi, then we have

It is clear that each vi is related to its accumulated error in Eq. 6. Smaller accumulated error means more contributions to the model for the sample point with larger vi, whereas the contribution of the sample point varies from small to nothing as ℓi increasing. Moreover, when ℓi is less than a small constant, the corresponding vi ∈ [0, 1] becomes 1, which enhances the influence of the less but faithful sample points. Conversely, vi ∈ [0, 1] is assigned as 0 and it works for reducing the effect of those unstable samples on the model. In other words, the accumulated error actually reflects the sample accuracy in position, which can be used to measure its importance to the model. Generally, higher accuracy means the point is more important and the sample point would be more suitable for the model. From this view, such strategy can select the faithful points in each iteration. As iteration, it makes our model preferentially select some reliable points and then gradually increase the points with larger influence on the model.

From Eq.6, each label υi is updated according to ℓi, which makes our model focus more on the reliable and accurate points in each iteration. By iterative point selection, our method gradually selects the points from reliable ones to less reliable ones, which are continually added into the final point set. The final points can be directly applied to estimate the displacement function f(x) better and more reliably. In other words, this is an accuracy-aware matching process. Meanwhile, the sum of the function κ(υi, ς) with respect to the label υi can be actually regarded as the matching number constraint and larger value of the sum means more matches.

3.3 Transformation Estimation

The Gaussian kernel is symmetric and positive define. It also makes the regularization term easy to rewrite in form under the Representer theorem. According to the theory of RKMS, the displacement function f can be expressed as a linear combination of multiple Gaussian kernels in RKHS. More details can be seen in [Ma et al., 2015b]. Then the function f takes the following specific form of Gaussian redial basis function:

where the coefficient ai is a column vector, and the Gaussian kernel Γ(xi, xj) = exp {-β‖xi -xj‖2} · I . Moreover, the regularization item Φ(f) is specified as the L2 norm, i.e.,

Φ(f ) = ‖f‖2Γ = 〈f , f 〉Γ . Then Eq. 3 can be further described as follows:

where Γ is the kernel matrix with Γi,j = π(xi, xj) = exp {-β‖xi -xj‖2}, and Γi represents the ith row of Γ. C = (a1, · · · , an)T is the coefficient matrix , and tr(·) denotes the trace of matrix, and the constant ℑ satisfies ℑ = ∑i υi ∑r wir .

3.4 Alternating Optimization

Since Eq.8 is not convex, it is almost impossible to find the global minimum by existing algorithms. Whereas Gaussian fields is differentiable and preferably convex in local area around the optimal position, a stable local extremum is enough for most applications. So, the gradient-based numerical approaches (e.g., quasi-Newton method, conjugate gradient method) can be used to solve it by deterministic annealing (DA). As iteration, a stable local minimum can be reached with a better chance [Wang et al., 2017]. There are two variables to optimize in Eq. 8. Here, an alternating optimization is used, namely, optimizing one variable with the other fixed in an alternating way. For υ, it is easy to obtain its optimal value. We only need to update it by Eq. 6 during iteration.

To optimize f , first ignore the items without the f because the derivatives of these items w.r.t f are 0. Then the equation becomes the following objective function Eσ(C) :

Since Eq. 9 is always local continuously differentiable w.r.t C, take its derivative w.r.t parameter C, then we have

where ρir = xTi +ΓiC-yT r σ . Here, Eq. 9 is convex only in the local neighborhood around the optimal solution. Thus, it is almost unlikely to directly get the global minimum. Using the derivative for Eq. 10, we solve the equation by the gradientbased quasi-Newton algorithm with DA.

During the optimization, a rigid-to-nonrigid technique is used to enhance the model convergence by performing DA on scale parameter σ2 . In other words, by varying σ from a larger initial value, it gradually makes Eq. 9 approach the local area around the real minimum. More specially, we initialize σ with a large value, and gradually update it with σ ↦→ νσ in each iteration, where ν is the annealing rate. Note that the criterion is convex over a large area around the global minimum at a large σ. As σ decreases, the position of the global minimum will tend towards smoothness. Considering that Eq. 9

Algorithm 1 Self-selected point matching

Input: Reference points with descriptors {xi, F(xi)}ni=1 , and sensed points with descriptors {yr, F(yr)}mr=1 , parameters β, λ, η, ν, n0, iterations Titer, µ Output: Optimal transformation τ

• 1: Calculate the Grammatrix Γ and matrix U;
• 2: Initialize parameter σ and the coefficient matrix Cs
• 3: Compute the distance between {F(xi)}ni=1 and {F(yi)}mr=1 as the weights {wir}n,m i,r=1
• 4: repeat (DA)
• 5: t ← t +1
• 6: Using the quasi-Newton algorithm (e.g. BFGS ) to solve Eq. 12 based on υ(t-1)
• 7: Update the coefficient matrix Cs = argmin Cs Eσ(Cs)
• 8: Update υ: calculate the υ(t) via the Eq. 6
• 9: Using the Cs, υ(t) to calculate the root mean square error rmse(t) of the selected points pairs
• 10: if ς < ςmax then
• 11: ς ←µς; ζ ←µζ
• 12: end if
• 13: Update annealing rate σ = νσ
• 14: until rmse(t) > min 1≤˜ τ≤t-1 rmse(˜ τ) or t ≥ Titer
• 15: Cs, υ∗ = argmin ˜ τ rmse(˜ τ)
• 16: Calculate the transformation τ (x) = x + f (x) by Eq.7
• 17: Return τ;

is convex in a local small area near its minimum, it is possible to converge to a new global minimum using the global minimum of last iteration as initial value of next iteration. So, the global minimum is more easier to obtain as iteration.

It is complex to directly use all the points in both time and space. A suboptimal but simpler method may be more effective. Kernel approximation is efficient for selecting an RKHS. Low-rank kernel approximation can yield a large speed improvements with little accuracy loss and constrains both the nonrigid transformation and space [Wang et al., 2017]. Specifically, low-rank kernel matrix Γs is the closest n0-rank matrix Γ, which satisfies the Frobenius norm,

Using the eigenvalue decomposition of Γ, the approximated kennel matrix can be taken as Γs = PΛP, where Λ is an n0 × n0 diagonal matrix with n0 largest eigenvalue and P is an n × n0 matrix with the eigenvectors. Finally, the regularized Gaussian criterion (9) and its derivative changes as

where ˜ ρir = xTi +ΓiCs-yT r σ , Cs denotes the new n0 × d coefficient matrix with elements (a1, . . . , an0 )T, Un×n0 = PΛ and Ui represents the ith row of U. We summarise the optimization procedure in Algorithm 1.

Note that the time complexity of Gaussian fields criterion is about O(n2m+n3 ). However, the time complexity reduces to O(nm) by the low-rank kernel matrix approximation.

4 Experiments

In this section, we assess our approach by comparing with several methods on different datasets, i.e., the synthesized 2D shapes dataset [Zheng and Doermann, 2006], the benchmark IMM Face Databaset 1 and the Oxford dataset 2 . Since different methods are proposed for different applications, we select the targeted methods on different datasets for a fair comparison. For the 2D shape and 3D face datasets, we make comparisons with the typical methods, mainly including CPD [Myronenko and Song, 2010], L2E [Ma et al., 2015b], RGF ([Ma et al., 2015c]), PR-GLS [Ma et al., 2016] and context-aware gaussian fields (CA-GF) [Wang et al., 2017]. On the Oxford dataset, some latest algorithms based on RANSAC [Lowe, 2004] are made comparisons, such as FSC [Wu et al., 2015] and the imprecise points removing algorithm (IPRA) [Wu et al., 2015]. All methods are implemented in Matlab, and tested on the same environment.

Parameter Setting In the experiments, a deterministic annealing technique strategy is applied on the scale parameter σ2 to improve the algorithm convergence. More specially, give a large initial value of σ2 , and reduce them with a fixed annealing rate ν by σ = νσ. We empirically set σ = 2, and ν = 0.93 throughout the whole experiments. The parameter β identifies the range width of the interaction between points. Since the points have been normalized to 0 mean and unit variance, β would be similar to different samples, and we set β = 0.2 in the experiments. The parameter of regularization term includes λ which works for trade-off the smoothness, and λ is affected by the degree of data degradation, which is set by λ = 0.1. Note that parameter η reflects the effect of feature descriptor on the model. So, it can be set as 1 for the same kind of data. Besides, we initialize Cs as zero vector.

Whereas the matching error is large at the beginning of iteration, the tolerance of the model to the error would be large for selecting more points. So, we dynamically update ζ, ς that decide the tolerance of the model to the error. As analyzed in Fig.1, we give a large initial value of ζ, ς and reduce them with a fixed rate µ. Since the optimization would be terminated, we choose a lower bound of ζ and set ζfinal = 1.0, and also set a lower bound ςfinal = 1.2 to control the weight of confused points. ς is fixed when ς ≤ ςfinal.

Low-rank kernel matrix approximation is applied to reduce the computational cost. Parameter n0 is the number of the selected eigenvalues, and as a trade-off, it controls the balance between runtime and matching accuracy. The proposed method is best for n0 ∈ [15, 20], where eigenvector P and eigenvalue are calculated by the fast Gauss transform (FGT)

1 Available at: http://www.imm.dtu.dk/ aam/datasets/datasets.html

2 Available at: http://www.robots.ox.ac.uk/vgg/data/dataaff.html

and its experimental analysis is shown in Fig.3. From it, we set n0 = 15 in the experiments.

4.1 Results on Non-rigid Synthesized Data

Fig.4 gives the matching results of our method on point sets with deformation, noise, outliers and rotation. In each group, four degradation levels are designed to test the robustness of our method. The intuitive experimental results in Fig.4 show that the model point sets are well aligned onto the data sets except for the data sets distorted by noises.

Moreover, for the data with different deformations, noises, outliers and occlusion, we list the visual results of the typical methods, e.g., CPD, L2E, RGF, PR-GLS and the latest CAGF (see Fig.5). From it, the shapes badly eroded by noises or outliers, are still well-aligned by our method. Especially for the occlusion, other methods almost fail or only align part of the shape while our method nearly gets a complete alignment. It owns to our soft selection. Initially, the points with large error are blocked by the soft labels, which makes the model sequential to select points, not to disturb other points. More reliable points are first selected to gradually enhance the model until it converges. Even if later points are noises or outliers, they are not enough to affect the estimated model. As presented in Table.1, our method has lower errors.

4.2 Registration on IMMFace Dataset

In this section, we test our method on the IMM Face dataset. It contains 240 images with 640 × 480 resolution on 40 human faces. For this data, each face has 6 samples and each sample contains 58 landmarks. In the experiment, our goal is to align the model points to the data landmarks. Due to complex facial expression, the deformation is nonrigid and also serious in local area. Aligning such points is difficult. However, Gaussian kernel can approach the real deformation better as shown in Fig.7. From the details of Fig.7, RGF can align the edges well, but it must be correct to single point for PSR. Since some landmark pairs are not the actual correspondences in real face, it is almost unlikely to align all the points. Such points indeed affect CPD, RGF and PR-GLS, and only obtains a rough alignment. As expected, our method selects more reliable points, and achieves an accurate alignment.

4.3 Registration on Real Image

The Oxford dataset is a widely used dataset in image registration. It consists of more than 40 image pairs with different variations. Fig.6 shows the results on the Oxford dataset. To test the validity of our method in real scenario, we further employ our approach to match the real images with complex deformations, which are captured at different time from different viewpoints and contain large scale changes, local selfsimilar texture even nonrigid distortions (see Fig.8). In this experiment, RGF, K nearest neighbor with the triangle-area representation (KNN-TAR) [Zhang et al., 2014], RANSAC, FSC and IPRA, is used to make comparisons.

5 Conclusion

In this paper, a robust nonrigid Gaussian fields criterion with self-selection is proposed for point matching. The main idea of the proposed method is a novel non-rigid transformation estimation with adaptive self-selected mechanism, and the biggest difference with traditional density estimation methods is that our method alters the convergence manner of model by adaptive sample selection. Moreover, an equivalent matching number constraint is embedded into the nonrigid criterion, which enhances the reliability of point matching. Experiments on three matching datasets show that our method has higher registration accuracy and correct matches than stateof-the-art approaches.

Acknowledgments

This work was supported by the National Key R&D Program of China under Grant 2017YFB1002202, National Natural Science Foundation of China under Grant 61773316 and 61379094, and the Open Research Fund of Key Laboratory of Spectral Imaging Technology, Chinese Academy of Sciences.