Flexible point cloud matching method based on three-dimensional image feature points

Abstract

Flexible and robust point cloud matching is important for three-dimensional surface measurement. This article proposes a new matching method based on three-dimensional image feature points. First, an intrinsic shape signature algorithm is used to detect the key shape feature points, using a weighted three-dimensional occupational histogram of the data points within the angular space, which is a view-independent representation of the three-dimensional shape. Then, the point feature histogram is used to represent the underlying surface model properties at a point whose computation is based on the combination of certain geometrical relations between the point’s nearest k-neighbors. The two-view point clouds are robustly matched using the proposed double neighborhood constraint of minimizing the sum of the Euclidean distances between the local neighbors of the point and feature point. The proposed optimization method is immune to noise, reduces the search range for matching points, and improves the correct feature point matching rate for a weak surface texture. The matching accuracy and stability of the proposed method are verified using experiments. This method can be used for a flat surface with weak features and in other applications. The method has a larger application range than the traditional methods.

Keywords

Point cloud match point feature histogram three-dimensional feature points double neighborhood

Introduction

Three-dimensional (3D) point cloud matching is a key task in inverse engineering.¹ Multi-view measurement point clouds are registered and integrated into the same coordinate system for a complete digital model. The solution of matching methods for point clouds is generally to establish descriptors for data association. These descriptors are constructed using a 3D pattern recognition algorithm, which can be used to achieve partial or entire surface matching. The solution is widely used in computer vision research for model construction, visualization, identification, classification, scene understanding, and so on. The method has the advantages of easy operation and extensive application.

At present, there are many research methods for point cloud matching.^2–9 The iterative closest point (ICP) algorithm³ is one of the most popular methods. The optimal transformation matrix of a two-view point cloud is acquired through continuously reducing the Euclidean distance between the point clouds until the best match between the two views is achieved. However, this method has a restriction that limits the initial position of the point cloud. A good initial position will greatly reduce the computational burden and improve the matching accuracy; otherwise, the algorithm can easily fall into local convergence, resulting in matching and registration errors. In order to precisely estimate the initial rotation and translation vector, one study⁴ developed a point cloud matching method that used topological features, including the border points and other salient points. Another⁵ introduced a mapping method based on a least squares conformal map, which transformed the 3D matching to two-dimensional (2D) matching. A third study⁷ derived a matching model using the overlapping relationship of different planes in the point cloud and then designed a point cloud registration model based on point and plane features. Tuladhar et al.⁸ developed a two-step registration approach for aligning partially overlapped 3D scanned data, and the fine registration is the optimized step for the coarse step. The method improved the registration efficiency and accuracy. Wei and Pu⁹ proposed a coarse matching method using irregular surface points and the central point of their neighborhood as matching points.

The above methods have two disadvantages. First, these methods are easily affected by noise, with a low robustness. Second, these low-dimensional metrics that include the curvature and distance are not very discriminative, but are very susceptible to matching error. For these reasons, many researchers have focused on 3D feature points with prolific details to describe the point cloud.¹⁰ Johnson¹¹ introduced a spin image feature descriptor. The point normal and tangent plane distance between one point and another are used as 2D feature points that describe the characteristic information around the point. Tombari et al.¹² used direction histogram signature descriptors to construct a topological structure for the spherical neighborhood of each point. Frome¹³ proposed a 3D shape context descriptor and established a spherical neighborhood structure for each point. The descriptor was divided into sub-regions. The invariant features were shown in the weighting descriptors. Rusu et al.^14,15 developed persistent feature histograms (PFH) and a fast point feature histogram (FPFH) method. This method has been used widely in computer vision research. Rusu et al.¹⁶ did not compute the feature for each point, but instead used the global descriptor viewpoint feature histogram (VFH), which joined viewpoint characteristic information based on the FPFH. Drost et al.¹⁷ used global description information based on directed feature matching points. A rapid voting mechanism was used in the local area for matching.

This study developed a new robust matching method based on double neighborhood constraint optimization. The proposed double constraint optimization method is not only immune to noise but also reduces the search range for matching points and improves the correct feature point matching rate for a weak surface texture. The remainder of this article is organized as follows. The next section introduces the principle of our method, which includes point cloud pretreatment, the intrinsic shape signature (ISS) algorithm, the FPFH algorithm, and our double constraint optimization method. Section “Experiments and analysis” presents the applications for the 3D registration problem using our double constraint optimization method and discusses how the efficiency of the method was tested on noisy scanned data using several experiments. Finally, we conclude and provide insight on our future work in section “Conclusion.”

Principle of our method

Point cloud pre-processing

Point cloud noise is inevitable in the measurement process. The matching results will have errors if the noisy points are not filtered. Thus, a filtering and interpolation operation is needed before matching. We use two simple methods to filter and interpolate noisy points. First, the distances between the point and its neighboring points are computed. If we suppose that these distances have a Gaussian distribution, a threshold value can be estimated, along with its mean and variance value. This threshold is used to separate the noisy points. Then, the curvature-based interpolation algorithm is used for the scanned data. The point cloud after pre-processing has a structured format and is acceptable to use for point cloud matching.

Key feature point detection based on ISS

The scanned point cloud has not only a huge quantity of data but also a different signature with a different view. The matching efficiency and accuracy are low with a huge point cloud. Thus, it is a good solution to extract the key feature points using a view-independent method, which are then used to compute the transform relation between different views. The ISS¹⁸ is a stable and discriminative method for characterizing the local region of a point cloud.

This method uses the popular surface normal reference to provide a full 3D coordinate system for local shape descriptors. The z-axis of the frame is the surface normal. The x- and y-axes are computed using the right-hand principle. A highly discriminative shape feature is computed as a weighted 3D occupational histogram of data points in its spherical neighborhood using a 3D partition that evenly divides the angular space at a 3D point. This choice of feature extraction ensures robustness in relation to the data noise, as well as any errors in computed reference frames.

First, the K-D tree is built for a 3D point cloud. Compute a weight for each point p_i inversely related to the number of points in its spherical neighborhood as follows

w_{ij} = 1 / ‖ p_{i} - p_{j} ‖, p_{j} \in D_{m}

(1)

where p_i is the center point in the pixel window mxm and p_j is the computed 3D point at the square D_m of the mxm pixel window.

Then, a weighted scatter matrix COV(p_i) can be determined for using all points p_j within a square D_m

COV (p_{i}) = \sum_{P_{j} \in D_{m}} w_{ij} (p_{j} - p_{i}) (p_{j} - p_{i})^{T} / \sum_{P_{j} \in D_{m}} w_{ij}

(2)

Finally, we can compute its eigen values ${λ_{i}^{1}, λ_{i}^{2}, λ_{i}^{3}}$ in the order of decreasing magnitude and their eigen vectors. The point with the relations of $λ_{i}^{2} / λ_{i}^{1} \leq ε_{1}$ and $λ_{i}^{3} / λ_{i}^{2} \leq ε_{2}$ is the key feature point, where $ε_{1}$ and $ε_{2}$ are the thresholds.

Computation of FPFHs

The computation of an FPFH at a point p relies on the presence of 3D coordinates and the estimated surface normal. For each point p, all of p’s neighbors are enclosed in a sphere within the k-neighborhood. Then, for every pair of points p and q, their estimated normals are n_p and n_q, respectively. Taking p as the coordinate origin, a Darboux uvn frame is built as follows: (1) $u = n_{p}$ , $v = u \times (p - q) / d$ , and $w = u \times v$ , where $d = ‖ p - q ‖_{2}$ is the Euclidean distance from p to q, as shown in Figure 1. The difference between $p (n_{p})$ and $q (n_{q})$ is shown as follows: (1) $α = \arccos (v \cdot n_{q})$ , (2) $ϕ = \arccos (u \cdot (p - q) / d)$ , and (3) $θ = \arctan (w \cdot n_{p}, u \cdot n_{p})$ . Then, the quad of $(α, ϕ, θ, d)$ , which is composed of p and its k-neighborhood point, is shown in the histogram in some statistical way. The PFH of point p is acquired.

Figure 1.

Reference frame.

Figure 2(a) presents an influence region diagram of the PFH computation for a query point (p). The p is marked in red and placed in the middle of a circle (sphere in 3D), and all its k-neighbors are fully interconnected in a mesh. It shows the angular characteristic relationship between p and pki. P6–P10 are the outside points of p’s k-neighbors. The FPFH is a simplified version of the PFH, which reduces the computational complexity of the algorithm, while retaining most of the discriminative power of the PFH.

Figure 2.

Computation model of descriptor: (a) descriptor of PFH and (b) descriptor of FPFH.

To simplify the histogram feature computation, the relationships between the point p and its neighbors, as shown in equation (3) and Figure 2(b), are computed for each query point p, which is called the simplified point feature histogram (SPFH).¹⁵ Then, each point’s k-neighbors are re-determined. The colorful lines in Figure 2(b) are the computation range of the SPFH. The neighboring SPFH values are used to weight the final histogram of point p

F P F H (D_{q}) = S P F H (D_{q}) + \frac{1}{k} \sum_{i = 1}^{k} \frac{1}{w_{k}} \cdot S P F H (D_{k})

(3)

Double neighborhood constraint descriptor

To find two matching point clouds, the sample consensus initial alignment (SAC-IA) algorithm is used to compute the feature histogram of the local neighborhood points.¹³ The point cloud is shown at the left in Figure 3, with the black dots depicting the local point clouds and the red dots showing the key feature points. A local neighbor is shown as $Ψ_{1}$ . The point feature histogram with SAC-IA cannot be used to sufficiently determine the characteristics of point p, because the range of the neighborhood is not large enough.¹⁹ Therefore, there are many wrong matches when it is used for a weak surface texture.

Figure 3.

Double neighborhood constraint model.

We developed a new method to solve this problem. As shown in Figure 3, $Ψ_{2}$ is a feature point neighborhood that presents the distribution of key feature points. Each key point uses its FPFH to show a wider neighborhood range because the calculation of the histogram descriptors and matching relies on the neighborhood point cloud distribution. Because the feature point has a view-independent characteristic, the key point’s FPFH for $Ψ_{2}$ has the ability to distinguish a weak feature texture. The neighborhood $Ψ_{2}$ covers the key point distribution of a wide local range, whereas neighborhood $Ψ_{1}$ covers a point cloud distribution for a small local range. Both neighborhoods are used for matching, which can greatly improve the robustness of point cloud matching.

The double neighborhood constraint steps consist of the following:

If points P and Q are key feature points detected from cloud points P₁ and P₂, respectively, we can compute every key feature point’s FPFH within its local neighborhood:

Three points are sampled randomly from point set P to form set $s_{p}$ .

The matched key feature point list of $s_{p}$ is then obtained from point set Q using its FPFH, forming set $M$ . The Euclidean distance between them is set at $d_{1}$ . Note that because there are many similar feature points on the object’s surface, the match relationship of list $M$ is not one to one.

For each pair of matched key feature points, the FPFH of neighborhood $Ψ_{2}$ is computed, and the Euclidean distance $d_{2}$ is obtained.

We then find three key feature points whose distance $d_{1} + d_{2}$ is the minimum for point set $M$ . These three points are set $s_{q}$ , and the transform matrix $T$ between $s_{p}$ and $s_{q}$ is estimated.

The match points of $M$ are eliminated because their match error is the maximum. Then, the match errors of the remaining points of set M are computed and set as $E_{R}$ .

If $E_{R} < E_{\min}$ , set $E_{R} = E_{\min}$ , and $T = T_{opt}$ .

Repeat steps 1 to 7. The algorithm will stop when any one of the following conditions is met:

(a)The maximum iteration number Max_Iterations is reached.

(b) $E_{\min} < E_{cov}$ , where $E_{cov}$ is the minimum permissible error estimated based on the model.

Experiments and analysis

Some experiments were performed to verify the stability and applicability of the method. Its immunity to noise, applicability to a weak texture surface, and matching accuracy were verified using the following experiments.

Noise immunity test and analysis

The rabbit model data were used in this experiment. Two point clouds P and Q were selected from a registered model. Point cloud Q was rotated 30° along a direction vector (10, 20, 0) and then the data were moved along the vector (0.05, 0.1, 0). The two point clouds with different views are shown as Figure 4(a). Then, different level noises were added to the model data. The SAC-IA-based algorithm and our double neighborhood constraint method were used for matching, as shown in Figure 4(b) and (c), respectively. We can clearly see that our matching method is better than the SAC-IA algorithm. There are some obviously gaps in the matching results of Figure 4(b), whereas Figure 4(c) shows the good matching results found with our method.

Figure 4.

Matching result comparison of different methods: (a) two-view point clouds, (b) matching result with SAC-IA method, and (c) matching result with double neighborhood constraint method.

Some Gaussian noise was added to the two resource point clouds. The matching results are shown in Figure 5 with different noise levels $(σ)$ . The number of correctly matched points decreases as the noise level $σ$ increases. The number of correctly matched points with the SAC-IA algorithm rapidly decreases from 108 to 30, but the number of correctly matched points with our method remained close to 75. Many correctly matched points were acquired even with the highest noise level $σ = 0.9$ , as shown in Figure 5.

Figure 5.

Correct match points comparison under different levels of noise.

Then, we compared the transform matrices with the two matching methods, as listed in Table 1. The transform matrices with the different methods were nearly same when the noise level was very low. The error increased with the SAC-IA method. Our method had good error convergence with an increase in the noise level.

Table 1.

Motion parameter comparison under different levels of noise (in mm).

Standarddeviation $σ$	SAC-IA						Our method
	$θ_{x}$	$θ_{y}$	$θ_{z}$	$x$	$y$	$z$	$θ_{x}$	$θ_{y}$	$θ_{z}$	$x$	$y$	$z$
0.1	0.2200	0.4500	0.5200	0.0018	0.0019	0.0008	0.2200	0.4500	0.3100	0.0007	0.0010	0.0008
0.3	0.2500	0.6500	0.7100	0.0018	0.0019	0.0008	0.2300	0.6500	0.4200	0.0007	0.0010	0.0008
0.5	0.3100	0.7100	0.7500	0.0019	0.0020	0.0017	0.2500	0.7100	0.4500	0.0009	0.0013	0.0020
0.7	0.3600	0.8100	1.5000	0.0019	0.0021	0.0038	0.3400	0.8100	0.4500	0.0012	0.0018	0.0020
0.9	0.5200	1.0400	1.8900	0.0020	0.0034	0.0040	0.3400	0.8110	0.4800	0.0012	0.0018	0.0021

Matching experiments for weak texture surface

In order to verify that our method would work for a weak texture surface, the smooth tooth model was used as a measurement object. Figure 6(a) and (b) shows the point cloud and 3D model surface of the tooth. It is obvious that the surface of the tooth is very smooth, with almost no features.

Figure 6.

Measurement data of tooth with weak texture surface: (a) point cloud and (b) 3D model.

The maximum iteration was set at $Max_Iterations = 1000$ , and various registration results were found with different values for the minimum registration error $E_{cov}$ . The matching results are shown in Figure 7(a)–(d). We can find that the false matches increased with the value of $E_{cov}$ . Conversely, the convergence condition for registration became stronger as the minimum registration error $E_{cov}$ decreased. Thus, the best match was acquired when the number of false matches was small. The minimum registration error $E_{cov}$ was one for the tooth data in this study.

Figure 7.

Matching results with different values for $E_{cov}$ : (a) $E_{cov} = 3$ , (b) $E_{cov} = 2$ , (c) $E_{cov} = 1.5$ , and (d) $E_{cov} = 1$ .

Verification of matching accuracy

Many experiments were used to verify the matching accuracy of the proposed method, with the results listed in Table 2. The advantage of our method is seen in the comparison of different matching parameters. The quantities of the correctly matched points were 97 for the gypsum jaw with the traditional method (SAC-IA), 105 for the original method without sample consensus (PFH), and 110 with our method. The random sample consensus (RANSAC) method was used to filter the matching results because there were many incorrectly matched points. The numbers of correctly matched points after filtering were approximately 61, 72, and 106. The average distances between two overlapping point clouds were 0.434, 0.432, and 0.216 mm, and the standard deviations were approximately 0.198, 0.191, and 0.106 mm. Thus, our double neighborhood matching method was found to be better than the traditional methods for a weak surface texture (Figure 8).

Table 2.

Experimental results.

Test model	Our method (RANSAC filter)	SAC-IA method(RANSAC filter)	PFH method (RANSAC filter)	Our method (average distance)	SAC-IA method(average distance)	PFH method(average distance)
Jaw	110 (106)	97 (61)	105 (72)	0.216000	0.434000	0.432000
Hand	33 (31)	31 (20)	32 (23)	0.003500	0.00970	0.009400
Cat	93 (93)	35 (10)	68 (13)	0.206000	0.557000	0.547000
Bunny	46 (40)	50 (28)	58 (29)	0.000513	0.000817	0.000817

RANSAC: random sample consensus; SAC-IA: sample consensus initial alignment; PFH: persistent feature histogram.

Figure 8.

Matching results with our double neighborhood method: (a) cat model, (b) correct match, and (c) matched point cloud.

Table 2 lists the detailed results with different methods. The correct matching rate was improved by an average of 33% with our method. The experimental data were analyzed using the Geomagic Studio 12 software.

Conclusion

This study developed a double neighborhood matching method with novel 3D feature points for difficult matching problems with a complex point cloud, random space position, high noise level, and weak surface texture. Both of the advantages of a local neighborhood and key feature points were used to eliminate incorrectly matched points. The correct matches were increased by 30% even if there was some noise. The experiment results showed that our method is superior, with a better correct matching rate and accuracy even when used for a weak feature surface. This method provides new opportunities for random point cloud matching, with a performance that is superior to the traditional matching method.

Footnotes

Handling Editor: Michal Kuciej

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the Fundamental Research Funds for the Central Universities of China (No. NS2017033).

ORCID iD

Haihua Cui

References

Tam

GKL

Cheng

Lai

et al . Registration of 3D point clouds and meshes: a survey from rigid to nonrigid. IEEE T Vis Comput Gr 2013; 19: 1199–1217.

Sang

Zhang

et al . Robust non-rigid point registration based on feature-dependant finite mixture model. Pattern Recogn Lett 2013; 34: 1557–1565.

Besl

McKay

. A Method for registration of 3-D shapes. IEEE T Pattern Anal 1992; 14: 239–256.

Kuang

. Application of neighborhood feature in point cloud registration. J Comput Appl 2012; 32: 762–765.

Wang

Jin

et al . Conformal geometry and its applications on 3D shape matching, recognition, and stitching. IEEE T Pattern Anal 2007; 29: 1209–1220.

Wan

Ruan

et al . 3D SMoSIFT: three-dimensional sparse motion scale invariant feature transform for activity recognition from RGB-D videos. J Electron Imag 2014; 23: 023017.

Zhang

Yao

Chen

. Research on building scanning point cloud registration based on geometrical feature constraint. J Geodesy Geodynam 2015; 35: 416–419.

Tuladhar

Lee

Ahn

. A novel coarse-to-fine registration approach for aligning partially overlapped 3D scanned data. Int J Comput Integ M 2017: 31: 718–726.

Wei

. Point cloud integration base on distances between points and their neighborhood centroids. J Image Graph 2011; 16: 886–891.

10.

Kim

Hilton

. Evaluation of 3D feature descriptors for multi-modal data registration. In: Proceedings of the international conference on 3D Vision, Seattle, WA, 29 June–1 July 2013, pp.119–126. New York: IEEE.

11.

Johnson

. Spin-images: a representation for 3-D surface matching. PhD thesis, Carnegie Mellon University, Pittsburgh, PA, 1997.

12.

Tombari

Salti

Di Stefano

. Unique signatures of histograms for local surface description. In: Proceedings of European conference on computer vision, Heraklion, 5–11 September 2010, pp.356–369. Berlin: Springer.

13.

Frome

Huber

Kolluri

et al . Recognizing objects in range data using regional point descriptors. In: Proceedings of European conference on computer vision, Prague, 11–14 May, pp.224–237. Berlin: Springer.

14.

Rusu

Blodow

Marton

et al . Aligning point cloud views using persistent feature histograms. In: Proceedings of IEEE/RSJ international conference on intelligent robots and systems, Nice, 22–26 September 2008, pp.3384–3391. New York: IEEE.

15.

Rusu

Blodow

Beetz

. Fast point feature histograms (FPFH) for 3D registration. In: Proceedings of 2009 IEEE international conference on robotics and automation, Kobe, Japan, 12–17 May 2009, pp.3212–3217. New York: IEEE.

16.

Rusu

Bradski

Thibaux

et al . Fast 3D recognition and pose using the viewp oint feature histogram. In: Proceedings of intelligent robots and systems (IROS), Taipei, Taiwan, 18–22October 2010, pp.2155–2162. New York: IEEE.

17.

Drost

Ulrich

Navab

et al . Model globally, match locally: efficient and robust 3D object recognition. In: Proceedings of computer vision and pattern recognition (CVPR), San Francisco, CA, 13–18 June 2010, pp.998–1005. New York: IEEE.

18.

Zhong

. Intrinsic shape signatures: a shape descriptor for 3D object recognition. In: IEEE 12th international conference on computer vision workshops, ICCV workshops, Kyoto, Japan, 27 September–4 October 2009, pp.689–696. New York: IEEE.

19.

Tombari

Salti

Di Stefano

. Performance evaluation of 3D keypoint detectors. Int J Comput Vision 2013; 102: 198–220.