A three-dimensional point cloud registration based on entropy and particle swarm optimization

A three-dimensional (3D) point cloud registration based on entropy and particle swarm algorithm (EPSA) is proposed in the paper. The algorithm can effectively suppress noise and improve registration accuracy. Firstly, in order to find the k-nearest neighbor of point, the relationship of points is established by k-d tree. The noise is suppressed by the mean of neighbor points. Secondly, the gravity center of two point clouds is calculated to find the translation matrix T. Thirdly, the rotation matrix R is gotten through particle swarm optimization (PSO). While performing the PSO, the entropy information is selected as the fitness function. Lastly, the experiment results are presented. They demonstrate that the algorithm is valuable and robust. It can effectively improve the accuracy of rigid registration.


Introduction
3D point cloud registration is widely applied in computer vision, computer graphics and so on. It is always a research hotspot in the computer field. In recent years, information science and technology have developed rapidly. [1][2][3][4][5][6][7][8][9][10][11] The field represented by artificial intelligence has achieved fruitful scientific research results. 12,13 Numerous descriptions have been proposed from rigid registration to nonrigid registration. 14 Rigid registration is a challenging work. It is because there is noise and unwanted points in the data and the original positions affect algorithm performance. Therefore, there is still much work to perform, such as removing the noise [15][16][17] and improving the registration accuracy and so on.

Rigid registration
For the original point cloud P and the target one Q, where there are a lot of overlapping between the two, the rigid registration is to find the rotation matrix R and the translation matrix T to transform the original point cloud to the target one. The equation is defined as 1 School of Information Engineering, Southwest University of Science and Technology, Mianyang, People's Republic of China where q i and p i are the corresponding points, R is the rotation matrix, and T is the translation matrix. For all kinds of rigid registration methods, most of them build registration algorithm models through various constraints.

Transform constraint
The nearest point method can be solved by establishing constraint potential consistency correspondence. It is to complete the final registration by selecting the nearest point as consistent correspondence point. The iterative closest point (ICP) is solved by establishing the following constraint condition Due to the good performance of ICP, many scholars have proposed many improved algorithms to improve the computation speed and robustness. [18][19][20] Feature constraint The geometric properties of point cloud, such as curvature and normal vectors, remain unchanged in rigid transformation. 21,22 Vary features can form an eigenvector. For the higher the eigenvector dimension, the lower the probability of all the features which can be matched, the high-dimensional eigenvectors can simplify data. Many surveys propose registration methods of feature constraints. [23][24][25] However, such method needs that the features of point cloud are obvious and easy to be extracted, and more time will be spent in the process of feature extraction. Kase et al. 26 proposed the extended Gaussian curvature, and a matching rate equation is used to determine the difference between the corresponding point sets. The extended Gaussian curvature is defined as where q i and p i are the corresponding points, and k 1 and k 2 are the main curvatures. The data can be greatly simplified and the registration efficiency is improved by extracting features.

Significance constraint
The significant areas are different from the surrounding areas. Significance can be used to measure local information to detect the key points or key areas. [27][28][29] The significance is usually used to reduce the potential mismatch points. The significance methods are usually geometric scale spatial analysis, significant scale and curvature based on vision, multi-scale sliding, maximum stability region extremum, and so on.

Regularization constraint
Regularization is constrained by adding penalty terms to the target function. Regularization constraint contains prior information, which avoids the occurrence of local minimum value and improves the search efficiency. Gold and Rangarajan 30 proposed that the rigid registration is regarded as a continuous optimization problem and dealt with problem between rigid transformation and consistency correspondence. The regularization term based on entropy is defined as where M is a consistency correspondence matrix.
According to the definition of entropy, when all points achieve registration, the entropy reaches the maximum.

Search constraint
Search constraint is mainly aimed at improving registration efficiency, including localization method and hierarchical search method. Jost and Hugli 31 proposed a method that speeds up the iteration of ICP with the above coarse to fine search technology and is refined gradually to obtain a more reliable consistency correspondence.

Contribution of this article
From the descriptions above, it is shown that 3D point cloud registration is the fundamental in computer field and still has many problems to further study and discussion. Our survey aims to solve the two rigid problems: (1) getting rid of the noise and (2) improving the accuracy. The challenge work in our article is how to improve the accuracy using the search constraint. The contributions of the article are introduced in the following: The k-d tree is used to build the relationship in point cloud to find the k-nearest neighbor (KNN) of the point, where the mean filter is employed to optimize the neighbor of the point. The noise of the point cloud is removed. The method of calculating the entropy of the point cloud is introduced in this article. It is used as the fitness function in particle swarm optimization (PSO) to search the best R. The results of the experiment show that our method effectively improved the accuracy.

Problem formulation
Point cloud registration is used to seek consistent correspondence between different datasets and to transform the different coordinate systems into the same coordinate system to gather the full data of object. Due to the limitation of 3D scanning technology, different datasets are usually obtained from different observation points. Each observation point is in a different coordinate system, so point cloud registration is an important aspect in 3D data acquisition. Point cloud registration is divided into rigid registration and nonrigid registration. Our method is the rigid registration and starts with formula (1) to find R and T. How to efficiently find R and T and improve the accuracy are the main problems solved by us. However, how to remove the noise is also a challenge work which we have to face.

Definition of entropy and particle swarm algorithm
We propose entropy and particle swarm algorithm (EPSA) to achieve rigid registration. The process of our algorithm is shown in Figure 1. First, the original and target point clouds are the input data. Then, the EPSA descriptor uses k-d tree to find the KNN of the point, where the mean filter is used to get rid of the noise. Second, the EPSA descriptor calculates the center of two point clouds to find the translation matrix T; the information entropy is defined to calculate the entropy of the two point clouds; and PSO, which is the search constraint, is employed to search the rotation matrix R at random. The stopping criterion in PSO is that the entropy of the original point cloud is the closest to the target one. Once R and T are found, we can achieve the registration.
k-d tree k-d tree is also named as k-dimensional tree. It is a highdimensional data structure and applied in searching the KNN such as the KNN matching in high-dimensional image feature vectors in the image retrieval, that is, optimization of KNN algorithm. The k-d tree is a highdimensional binary search tree. Unlike the common binary search tree, the tree stores the k-dimensional data.
In the searching process, the k-dimensional data need to decide from which dimension the data are selected. Then, the data are compared with the root node in the selected dimension. How to select the dimension and ensure the number of nodes in the left subtrees is as equal as possible to the right one after selecting dimension are two challenge works to face. In order to select the dimension, the variance in each dimension is calculated and the dimension of the largest variance is selected; for the larger the variance, the more the data disperse. The right and left subtrees are divided according to the pivot value to ensure the number of nodes in the left subtrees is as equal as possible to the right one. Table 1 shows the data structures in each node of the k-d tree.
From Table 1, it can be seen that building k-d tree is a hierarchical recursive process. The pseudo-code is shown as follows: Name: CreateKDTree. Input: Data set. Output: K d tree.
Step 1. If Data set is empty, returns an empty K d tree.
Step 2. Call the node generating function.
(a) Split is determined: the data variance in each dimension is calculated. Then, the corresponding dimension in which data variance is the maximum is the value of Split. The large data variance indicates that the data along the coordinate axis are scattered widely, and the data segmentation in this direction has a better resolution. (b) Node data is determined: the data set is sorted by Split. The data in the middle is selected as Node_ data, Data_set' = Data_set\Node_ data, where it means getting riding of Node_data in Data_set.
Step 8. right = CreateKDTree(dataright,Right_Rang), the parent of the right is set to K_d tree.
In this article, k-d tree is used to establish the relationship between points to find the KNN of the point. Once the KNN of the point is found, we use mean filter to get rid of the noise. For KNNs of point p i1 , p i2 , . . . , p ik , k = 1, 2, . . ., of p i , the value q i by mean filtering is defined as Therefore, p i = q i .

Information entropy
The definition of information entropy is proposed by Shannon. 32 It is applied in calculating the average amount of the information. The information associated with the probability of data. The lower the probability, the more is the information which the event carries. The amount of the information conveyed by each event depends on the random variable's value. The value is the information entropy.
In this article, the information entropy is developed to calculate the entropy of point cloud. First, we need to project the point cloud to the X , Y , Z axis, respectively. The information entropy of X , Y , Z is the average amount of information. The formula is defined in the following  where p x is the probability of the 3D point clouds which are projected to X-axis, and it is related to frequency. We need to standardize data in order to calculate the frequency. The standardizing formula is defined in the following 3 1000 mod 1000 S(y i ) = y i À min(y) max(y) À min(y) 3 1000 mod 1000 where S(x i ) is the standardized value within a region (0-1000) in X-axis when the point cloud is projected to the X-axis. max() is a function to find the maximum, and min() is a function to find the minimum. Similarly, we have S(y i ) and S(z i ). The frequency is obtained (See Figure 2), and it is to count the number of times which the value appears between 0 and 1000. Then, p x is obtained, which is related to frequency and calculated as where n xi , n yi , n zi is the number of point cloud which is projected to the X , Y , Z axis and N is the total number of point clouds. The same procedure may be easily adapted to get H(y) and H(z). H sum is defined as the sum of information entropy when the 3D point clouds are projected to the X , Y , Z axis PSO To find the rotation angle to generate point-to-point correspondence, we study the optimization algorithm [33][34][35] and we research PSO. PSO is a populate algorithm to solve the optimization problems. 36 The particles move in the whole search space. A candidate solution is represented by a particle. According to the rules, each particle searches the better position. The basic formula of PSO is defined in the following where D is the spatial dimension; X is the particle i's position; V is the particle i's velocity; p i is the particle i's best position in history; p g is the best position for all particles; c 1 , c 2 are the learning factors or accelerators, which are usually set to 2; and r 1 , r 2 are the pseudorandom numbers with a region (0-1). The particle of PSO in this article starts by generating random rotation angle within an initialization region (0°-360°). The particle's velocity is usually set to zero or to small random values in order to make the particle move in the search space during the first iterations. Once the stopping criterion is met, the velocities and rotation angles of the particles are no more updated in the algorithm. The stopping criterion in the algorithm is that the H sum from the source point cloud by rotating is closest to the H sum from the target one. The flow diagram is shown in Figure 3, where R x is a rotation angle around the X-axis. Similarly, we have R y and R z .

Experimental simulation
This experiment is performed in MATLAB R2012b and the datasets are from Stanford's experimental database.

Robustness
In order to assess the robustness of the noise, the experimental dataset contains two models (''cow'' and ''feet of man''). In the test, we use eight nearest neighbors of each point and add the random noise to the point cloud. The result is shown in Figure 4. From the results, we find that the descriptor is workable when facing such a challenging case.

Registration accuracy test
In order to test the registration accuracy, the experimental dataset contains three models (''bunny,''''cow,'' and ''man''). The dataset used in complete overlapping is shown in Figure 5 and in partial overlapping is shown in Figure 6. The red one is the original point cloud and the blue one is the target point cloud.
Performance of EPSA. We calculate the mean square error (MSE) to further test the performance of EPSA. MSE is defined in the following where MSE is the mean square error; x ok , y ok , z ok is the source point cloud; x tk , y tk , z tk is the target one; and n is the smallest size of the point clouds.
In order to make the data in experiment, first, we calculate the center of gravity of two point clouds. For the point cloud P = p i x i , y i , z i ð Þji = 1, 2, 3, . . . , N f g , the gravity is defined as where mean() is the mean function. Then, T, which is the difference in the gravity center between two clouds, is obtained. Second, the entropy of the two point clouds is calculated, and PSO is used to search the R. The experimental flow is shown in Figure 7 and the steps are shown as follows: Input: The original point cloud and the target point cloud.
Output: The rotation matrix R g .
Step 1. Set the number of iterations, M, and the number of randomly generated rotation angles, N.   Step 2. Set c 1 = 2, c 2 = 2, where c 1 , c 2 are the learning factors or accelerators.
Step 4. The original point cloud is transformed in R i , i = 1, 2, . . . , N , to the target one, respectively.
Step 5. The differences e i , i = 1, 2, . . . , N , between the entropy of the two point clouds are calculated after rotation, respectively. Step 6. Find the minimum e in e i , i = 1, 2, . . . , N , then the corresponding R in the minimum e is the best R in this round.
Step 7. Set R p = R, R g = R, where R p is the particle's best position in history, and R g is the best position for all particles.  where r 1i , r 2i are the pseudo-random numbers and rand() is the random function, which generates a random number between 0 and 1 from a uniform distribution.  Table 2 shows the size of the experimental datasets and the MSE by comparing it with the descriptor coherent point drift (CPD). 37 Case of the point cloud partially overlapped. Case of the point cloud partially overlapped is that the data of the point cloud are only part of the object. The results in partial overlapping are shown in Figure 11, and Table   3 shows the MSE of the experimental datasets which are compared with the descriptor CPD.     Tables 2 and 3, we find that the experimental datasets achieved the high correct rate followed by the proposed descriptor.

Conclusion
In this paper, a 3D point cloud registration based on entropy and particle swarm algorithm is proposed. Firstly, to find the k-nearest neighbor, the k-d tree is employed in the point cloud to establish the relationship in the points. The noise is suppressed by the mean of neighbor points. Secondly, in order to improve the registration accuracy, the gravity center of two point clouds is calculated to find T and the search constraint is used to find the best R. Once the R and T are found, the original point cloud can be transformed to the target one. Lastly, the experiment results are presented. They demonstrate that the algorithm is workable.

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 jointly supported by the National Natural Science Foundation of China (grant nos 11705122   MSE: mean square error; EPSA: entropy and particle swarm algorithm; CPD: coherent point drift.