Invariant Feature Point based ICP with the RANSAC for 3D Registration

INTRODUCTION

The task of registering three dimensional data sets is a fundamental problem in reverse engineering, facility mapping, computer vision, arising whenever two or more 3D data sets must be aligned in a common coordinate system. This study addresses the usage of IFP-ICPR to reduce the computation time and improve the convergence of ICP in registering 3D data. We present a new surface curvature based method for invariant feature point extraction, a robust ICP process with outlier rejection and an experimental evaluation demonstrating reduction in computation time and improvement in convergence.

3D registration is typically accomplished by a variant of the ICP algorithm (Besl and McKay, 1992). ICP is an iterative descent procedure which seeks to minimize the sum of the squared distances between all points in the reference data and their closest points in the input data. Chen and Medioni (1992) independently developed an approach similar to ICP, which minimized the sum of squared distance between reference data points and a local planar approximation of the input data points. However, there are two main drawbacks for the above two classic ICP algorithms: (1) Convergence to local optimal; (2) high computation time for large datasets.

To make ICP converge to good solution, spin image (Johnson, 1997), principal component analysis (Kim et al., 2003), RANSAC-based darces (Chen et al., 1999) are used as coarse registration for good initial transformation parameter estimation; rejection of the worst pairs as outliers based on point-to-point or point-to-surface distance (Rusinkiewicz and Levoy, 2001) are used to ensure the correct direction of convergence; spherical harmonics invariant features (Sharp et al., 2002) are added into ICP iterations. As coarse registration only finds the initial parameters, then the convergence of ICP is mainly determined by the matching during ICP iterations, where the outlier removal is important. To reduce the computation cost, different sampling strategies are used: uniformly sampling (Turk and Levoy, 1994), random sampling, normal space sampling (Rusinkiewicz and Levoy, 2001). For good sampling strategies, the shape of the object is well described by small amount of subset of the 3D data. In this study, a new surface curvature based invariant feature point extraction is used as the sampling method to cut down the computation cost effectively. As the extracted invariant feature points are more accurate to describe the shape of 3D data and the embedded RANSAC removes outliers effectively at each iteration of ICP, the proposed IFP-ICPR ensures the convergence to good solution.

BACKGROUND STUDY

Iterative closest point registration: To describe 3D registration with ICP, a general statement is described as follow. Given two point sets in U³, one denotes the reference 3D data:

and the other is the input 3D data:

The standard ICP algorithm achieves rigid registration by iteratively registering the input data to reference data with rigid transformation and it has two basic steps.

Firstly, find correspondence:

between data P and X the rigid transformation (R^k-1, t^k-1) after the (k-1)th iteration:

(1)

Secondly, compute the kth rotation and translation transformations (R^k, t^k) based on the current correspondence:

	(2)

Surface normal vector estimating using covariance analysis: For the estimation of the surface normal vector in a 3D data, the covariance matrix, sometimes called the variance-covariance matrix, will be utilized since the first order plane fitting is equivalent to the eigenvalue problem of the covariance matrix (Bae, 2006). Figure 1 presents a point p_i and its neighbor points q_j(j = 1,2,...,N) used for estimating surface normal vector. The covariance matrix COV(p_i)∈U³xU³ is expressed as:

(3)

Fig. 1:	First order plane fitted by neighbor points

where, c_i is the center of q_j and e_m is the eigenvector of the (m+1)th smallest eigenvalue λ_m. Since, COV(p_i) is a real, positive, semi-definite matrix, its eigenvalues are always greater or equal to zero (Bae, 2006). The eigenvector of the minimum eigenvalue is the estimated normal vector of the surface formed by p_i and its k neighbor points.

Change of surface curvature estimation for invariant feature point extraction: The common way for invariant feature point extraction is to utilize the change of surface curvature. There are many ways to define curvature, e.g., through Gaussian and mean curvature matrix. Hoppe et al. (1992) proposed a covariance analysis method for the estimation of the normal vector with consistent orientation. The covariance analysis method has been also utilized for the estimation of local curvature estimation using the ratio between the minimum eigenvalue and the sum of the eigenvalues. It has been called as the surface variance by Pauly et al. (2002).

Each eigenvalue of the covariance matrix represents the spatial variation along the direction of the corresponding eigenvector. The ratio of the minimum eigenvalue and the sum of the eigenvalues approximates the change of curvature, M_c(p_i), as follows:

(4)

INVARIANT FEATURE POINT BASED ICP WITH THE RANSAC ALGORITHM

Modified surface curvature estimation for invariant feature point extraction: The classic surface curvature estimation in Eq. 4 depends on the three eigenvalues in Eq. 3, if there is position noise in 3D data, the eigenvalue λ_i(i = 0,1,2) will be biased, then the curvature is incorrect to describe the surface. One way to solve this problem is to use some eigenvalue in λ_i(i = 0,1,2).

Fig. 2:	3D curvature computation with sphere radius

Wang et al. (2009) only used e₀ to construct a sphere and then compute the curvature. Here, we modify the construction of sphere by e₀ and the fitted plane and the radius of sphere is set as the surface curvature.

Figure 2 shows the point p_i with its estimated surface curvature M(p_i) = 1/r. The distance of point p_i to its tangent plane is:

(5)

The surface is represented by the spherical patch cutting by the fitted plane. According to s² =μ²-d² and s² = r²-(r-d)², the curvature is computed as to be a good criterion for detect invariant feature point of crease and corner points,

(6)

By computing the maximum curvature estimation M_max(p) of all points, the curvature is defined as:

(7)

where, M_n(p_i) is in the range of [0,1]. When M_n(p_i) near to 1, the possibility to be crease or corner point is high.

Outlier removal procedure with the RANSAC: When the registration error is small, it is expected to be a set of very reliable corresponding points. However, it does not mean that there are no outliers, which precludes the registration algorithm from estimating accurate relative transformation parameters between points. Therefore, a method for removing outliers is required.

Many procedures of outlier detection use as large number of data point as possible and then detect outlier among them. However, the RANSAC procedure tries to use as small number of data as possible if the problem allows:

In the case of RANSAC, the process is repeated for K iterations and the candidate solution with the highest number of inliers based on squared residual error is selected as the best transformation model. In the case of RANSAC, the number of iterations can be determined as:

(8)

where, p is the desired probability of selecting at least one transformation model that is free from outliers within K iterations, n_inliers is the number of candidates that fit the current estimated model, N is the total number of candidates and S is the minimum number of candidates needed to fit the transformation model.

A description of the proposed method, invariant feature point based ICP with the RANSAC(IFP-ICPR), is presented in this section. As ICP converge to local optimal, initialization of the transformation model parameters is need. However, as the proposed IFP-ICPR has a good converge range, a common coarse registration method is enough. So, we select PCA (Kim et al., 2003), which is very fast without identifying point or curve correspondences compared to other coarse registration algorithms.

The pseudo code of the proposed IFP-ICPR is presented as follows:

Algorithm: IFP-ICPR

The flowchart of the algorithm is shown in Fig. 3.

Fig. 3:	Flowchart of IFP-ICPR

EXPERIMENTS AND PERFORMANCE EVALUATION

To evaluate the performance of IFP-ICPR, we conducted two groups of experiments: (1) point extraction experiments to validate the modified invariant feature point extraction; (2) overall registration experiments with simulated and real reconstructed 3D data to show the registration performance.

Point extraction experiments: The 3D data used for feature point extraction is the simulated cube data, the total point number is 2646 and number of the crease and corner points is 447. We use 3D rigid and affine model to simulate two different transformations, then the invariance of the feature point extractor is tested. To do comparisons, the surface variance method (Pauly et al., 2002) is selected, the number of local neighbor N and th_M are set to different values. The point extraction results are given in Fig. 4. To evaluate objectively, Fig. 5 give the correct ratio of the extracted 3D feature point under different curvature threshold th_M and N.

Comparing the number of extracted points in Fig. 4, we can see that surface variance based method extracted more points than ours with original, rigid transformed and affine transformed data. However, if we look into the extracted points, especially between Fig. 4e and f, our extracted points are more correctly and sparsely located at the creases and corners than surface variance based method. As the shape of object can be more easily described by creases and corners, our modified algorithm is better to describe the shape of object. Considering Fig. 4b, d and f, the correct/entire numbers of extracted points almost remain the same, but the numbers change from 210/261 to 92/180 in Fig. 4c and e, then we can draw the conclusion that our modified algorithm is much less conditioned by transformation, that is to say more invariant than surface variance based method.

To evaluate the performance under different curvature threshold th_M and N objectively, Fig. 5 gives the correct ratio of the extracted points. Considering the five curves (th_M =0.2, 0.3, 0.4, 0.5, 0.6, respectively) in Fig. 5, the range of correct ratio for surface variance based method is between 0.5 and 1, while the range of our modified method is between 0.7 and 1, so we can make the conclusion that the correct ratio of our modified algorithm is higher than the surface variance algorithm (approximately 20%). To compare every curve with the changing of N, the curves of our modified method is more flat than surface variance algorithm, thus our modified algorithm is more stable than the surface variance algorithm, which means we can select smaller N (we set N = 9 in the following experiments).

3D registration experiments: We select original ICP (Besl and McKay, 1992), uniformly sampled ICP (US-ICP, Turk and Levoy, 1994), invariant feature point based ICP (IFP-ICP, remove the RANSAC out of the proposed IFP-ICPR algorithm) and the proposed IFP-ICPR for comparisons.

To make the results more comparable, PCA is used to initialize R⁰ and t⁰ for the 4 algorithms. The maximum iteration number N is set to 30, ICP convergence tolerance τ is set to 10^-1. For the proposed algorithm, th_M and Δth_M are set to 0.6 and 0.1, respectively, the stopping threshold th_ε is set to the distance of the nearest of points in reference or input points.

To compare the 3D registration results, subjective and objective evaluation are used. For subjective evaluation, the registration results of input data are shown in the same coordinates with the reference data. If the registration results match the reference data smoothly, the registration is good, otherwise it is poor. For objective evaluation, RMSE (Root Mean Square Error) is used and is defined as follows:

Fig. 4:	Invariant feature point extraction from simulated cube data under N =15 and curvature threshold thM = 0.5 (correct numbers/entire numbers). (a) Surface variance (215/273), (b) Our modified algorithm (135/135), (c) Surface variance (210/261), (d) Our modified algorithm (134/134), (e) Surface variance (92/180) and (f) Our modified algorithm (128/132)

Fig. 5:	Correct ratio of the extracted points with the variation of thM and N

Table 1:	Parameters estimation/RMSE/computing time of simulated cube data registration

(9)

where, (x^R_i, y^R_i, z^R_i) and (x^I_i, y^I_i, z^I_i) are the ith pair of test point, the total number M is set to 10 and T(•) is the estimated transformation.

Simulated cube data registration: For the cube data, we use rigid transformation to simulate the input data, the truth 6 parameters (ω, φ, κ) are the rotated degrees in X, Y, Z axis, respectively and ΔX, ΔY, ΔZ are the translations are shown in Table 1.

From the registration results in Fig. 6, we can see obvious mismatched part in Fig. 6a and b, which means that original ICP and US-ICP fail to register the input data and the registered results of IFP-ICP and IFP-ICPR match the reference cube data well in Fig. 6c and d. However, according to the estimated parameters in Table 1, the error of estimated 6 parameters (ω, φ, κ, ΔX, ΔY, ΔZ) of our IFP-ICPR achieves the minimum among the four algorithms and the RMSE is 0.33, which is much better ICP, US-ICP and IFP-ICP. Considering the RMSE curves of the four algorithms over iterations in Fig. 7 , IFP-ICP and the proposed IFP-ICPR converge dramatically, however, only our IFP-ICPR converge to smaller RMSE.

The time required is very important, all the methods have been programmed using Matlab in a Pentium Dual Core 2.2 GHz because MATLAB guarantees an easy implementation. As ICP takes all data points into consideration, the computation time is 74.24 sec. The time cost of US-ICP, IFP-ICP are one order of magnitude smaller than ICP and our IFP-ICPR achieve one more orders magnitude smaller at 0.12 sec.

Fig. 6:	Registration results of simulated cube data. (a) Original ICP, (b) US-ICP, (c) IFP-ICP and (d) IFP-ICPR

Fig. 7:	Registration RMSE of simulated cube data over iterations

Real reconstructed data registration: The real reconstructed 3D data are shown in Fig. 8. Compared to the simulated cube data above, the real reconstructed data are different: (1) The true parameters of rigid transformation model are unknown, (2) the data size is almost 11 times larger than the cube data and (3) 3D changes exist between reference and input data, which result in more outliers during ICP iterations.

According to the registration results in Fig. 9, small translation error exists in the registration results of ICP, US-ICP and IFP-ICP, only our IFP-ICPR register the input data to the reference data correctly, where the RMSE in Table 2 shows the same results. As there are changes of the 4th and 6th building of the input 3D data in Fig. 6b, which induces the existence of outliers in real reconstructed data, IFP-ICP does not converge to good solution as in Fig. 10 so RMSE curve of IFP-ICP is much above IFP-ICPR RMSE curve.

Fig. 8:	Real reconstructed data with multiple objectswith 3D changes. (a) Reference data and (b) input data

Fig. 9:	Registration results of real reconstructed data. (a) Original ICP, (b) US-ICP, (c) IFP-ICP and (d) IFP-ICPR

Table 2:	Registration RMSE and computation time of multiple whole objects data

Fig. 10:	Registration RMSE of real reconstructed data over iterations

As the data size is larger, US-ICP and IFP-ICP cut down the computation time by two orders magnitude and the time for our IFP-ICPR is only 3.81 sec, which is one order magnitude less than US-ICP and IFP-ICP.