A Closed-Form Solution to Planar Feature-Based Registration of LiDAR Point Clouds

: Since pairwise registration is a necessary step for the seamless fusion of point clouds from neighboring stations, a closed-form solution to planar feature-based registration of LiDAR (Light Detection and Ranging) point clouds is proposed in this paper. Based on the Plücker coordinate-based representation of linear features in three-dimensional space, a quad tuple-based representation of planar features is introduced, which makes it possible to directly determine the difference between any two planar features. Dual quaternions are employed to represent spatial transformation and operations between dual quaternions and the quad tuple-based representation of planar features are given, with which an error norm is constructed. Based on L2-norm-minimization, detailed derivations of the proposed solution are explained step by step. Two experiments were designed in which simulated data and real data were both used to verify the correctness and the feasibility of the proposed solution. With the simulated data, the calculated registration results were consistent with the pre-established parameters, which veriﬁes the correctness of the presented solution. With the real data, the calculated registration results were consistent with the results calculated by iterative methods. Conclusions can be drawn from the two experiments: (1) The proposed solution does not require any initial estimates of the unknown parameters in advance, which assures the stability and robustness of the solution; (2) Using dual quaternions to represent spatial transformation greatly reduces the additional constraints in the estimation process.


Introduction
With the fast development of Light Detection and Ranging (LiDAR) techniques and their successful application in three-dimensional data acquisition, point cloud registration has attracted significant attention for its role in the fusion of LiDAR point clouds from two neighboring stations. The essence of point cloud registration is to estimate the transformation parameters between the two neighboring stations, which is also known as spatial transformation. As is known, a spatial transformation can be explained as a rotation around the x, y, and z axes, a translation along the three axes, and a scale factor based on the centroid of the coordinate system.
Based on the different registration primitives used for the estimation of unknown transformation parameters, the available methods can be categorized into point feature-based methods [1,2], linear feature-based methods [3][4][5], planar feature-based methods [6][7][8], and hybrid feature-based methods [9][10][11]. Until now, point features are still the most popular and widely used registration primitives for simple mathematical expressions. However, affected by the characteristics of LiDAR technology, the extraction of point features from point clouds often has low accuracy without pre-established man-made reflectors. Except for point features, linear features are another popular registration primitive. Compared Table 1. Advantages and disadvantages of the iterative methods and the closed-form methods.

Registration Methods Advantages Disadvantages
Iterative Methods (1) The most popular and widely used; (2) The derivation of the formulas is simple.
(1) Initial approximate estimates of those unknown parameters must be determined in advance; (2) The number of iterations is closely related to the choice of those initial estimates.

Closed-form Methods
(1) No initial estimates of the unknown transformation parameters are needed in advance; (2) Registration results can be obtained in only one step.
(1) The derivation of the formulas is complex.
Based on the abovementioned analysis, a closed-form solution to planar feature-based registration of point clouds is proposed. Firstly, a quad tuple-based representation of planar features is given and dual quaternions are then employed to represent the spatial transformation. After the operations between the dual quaternions and the quad tuple are explained, an error norm (error function) is constructed by supposing that the two conjugate planar features are equivalent after registration. Based on L2-norm-minimization, detailed derivations of the proposed solution are given step by step. Lastly, two experiments are designed to verify the correctness and feasibility of the proposed method, in which simulated and real data are both incorporated.
The remainder of the paper is organized as follows. Section 2 reviews some related work. Section 3 gives the operations between dual quaternions and the quad tuple-based representation of planar features in three-dimensional space. Section 4 explains the proposed solution, in which the detailed derivations of all formulas are given step by step. Section 5 shows the experiments and the results. Section 6 discusses the proposed solution and gives suggestions for future work. Section 7 concludes the paper.

Quaternion and Its Application in Point Cloud Registration
Quaternion was first proposed by Hamilton W. R., which has been proved to be a convenient and effective way to describe rotation in three-dimensional space [17]. Due to its compactness and high efficiency, quaternion has attracted considerable attention from researchers working in different fields. The representative work can be summarized as follows: Horn introduced unit quaternion to solve the absolute orientation problem in photogrammetry [13]; Shen et al. used unit quaternion in the transformation of two sets of three-dimensional coordinates [12]; Zeng et al. presented a unit quaternion-based, iterative solution to coordinate transformation in geodesy [18]; Joseph et al. introduced unit quaternion in robot arm manipulation and presented an extended Kalman filter-based algorithm for the estimation of human motion [19]; Kim et al. presented a similar algorithm to Joseph et al., which employed unit quaternion for the real-time estimation of orientation in robot arm manipulation [20].
As is known, spatial transformation mainly consists of a rotation and a translation. However, unit quaternion can only represent the spatial rotation in three-dimensional space, as shown in Figure 1. Later, dual quaternion was presented to represent spatial rotation and translation simultaneously, as shown in Figure 2, in which two quaternions are integrated with the aid of a dual number. The first successful application of dual quaternion in estimating the unknown spatial transformation parameters was introduced by Walker et al. [15]. Based on the analysis of correspondences between dual quaternionbased and matrix-based representations, a single cost function was formulated, which enabled the simultaneous calculation of six parameters in point feature-based registration. Instead of estimating only the six transformation parameters, Wang et al. added the scale parameter to a single cost function, which enabled the simultaneous derivation of rotation, translation, and scale parameters [16]. Prošková et al. also introduced dual quaternion to represent spatial transformation; the presented approach was successfully applied for deriving the seven parameter-based transformation between two sets of three-dimensional coordinates [21,22]. of three-dimensional coordinates [12]; Zeng et al. presented a unit qua ative solution to coordinate transformation in geodesy [18]; Joseph et quaternion in robot arm manipulation and presented an extended K algorithm for the estimation of human motion [19]; Kim et al. prese rithm to Joseph et al., which employed unit quaternion for the real-tim entation in robot arm manipulation [20].
As is known, spatial transformation mainly consists of a rotatio However, unit quaternion can only represent the spatial rotation in space, as shown in Figure 1. Later, dual quaternion was presented rotation and translation simultaneously, as shown in Figure 2, in whi are integrated with the aid of a dual number. The first successful appl ternion in estimating the unknown spatial transformation parameters Walker et al. [15]. Based on the analysis of correspondences betwee based and matrix-based representations, a single cost function was for abled the simultaneous calculation of six parameters in point feature Instead of estimating only the six transformation parameters, Wang e parameter to a single cost function, which enabled the simultaneous tion, translation, and scale parameters [16]. Prošková et al. also intro nion to represent spatial transformation; the presented approach was s for deriving the seven parameter-based transformation between two s sional coordinates [21,22].
n , then the rotation can be expressed as Based on the above analysis, conclusions can be drawn that unit quaternion can only represent rotation in three-dimensional space; when it is applied in spatial transformation, such as the seven parameter-based Helmert transform, the rotation parameters must be estimated first, on which basis the translation parameters and the scale factor will later be estimated. In the case of errors being introduced into the estimation of the rotation parameters, the accuracy of the translation parameters and the scale factor will certainly be affected, which makes it a tendency that dual quaternions are introduced to represent spatial transformation.

Planar Feature-Based Registration Methods
Compared to point and linear features, planar features exist widely in those point clouds acquired from man-made buildings. On the condition that no specific man-made reflectors are used, point and linear features are often extracted by fitting and the intersection of adjacent planar features. Therefore, if planar features can be used directly in point cloud registration, the amount of calculation will be reduced greatly. In other words, the employment of planar features in point cloud registration will provide more conditions for estimating the transformation parameters than using point and linear features only. Of all existing planar feature-based registration methods, some use the distance between the point and its corresponding planar feature to construct error norms [6,23], others use the parallelism between the two normal vectors and the distance between the two conjugate planar features to construct error norms [7,[24][25][26][27][28]. Detailed explanations are as follows: Grant et al. presented an iterative solution to pairwise point cloud registration, which is based on the correspondences between point and planar features [23]. Pavan et al. presented a closed-form solution to planar feature-based registration of terrestrial LiDAR point clouds [27], which was later applied in the global refinement for the terrestrial laser scanner (TLS) data registration [7]. Moreover, Wang et al. [24], Zhang et al. [25], and Previtali et al. [26] separately presented a planar feature-based registration algorithm based on the correspondences between each pair of planar features in which the Rodriguez matrix, Euler angles, and quaternions were, respectively, employed to represent spatial rotation. In addition, Khoshelham presented a closed-form solution to planar feature-based pairwise point cloud registration [6] in which the Kronecker product was employed to linearize the nonlinear equations; rigid transformation and similarity transformation are both given in detail; the parallelism between normal vectors of the two conjugate planar features was not fully utilized. Föstner and Khoshelham presented an optimal solution and three direct solutions for efficient motion estimation from plane-to-plane correspondences and provided an analysis of the accuracy of the solutions, comparing their performance with the classical iterative closest point (ICP) algorithm [28].
In general, few of the available methods incorporate the scale parameter in planar feature-based registration, which makes it difficult to apply them in the fusion of data from different sources. Moreover, some methods use point and planar features as registration primitives; the role of normal directions is neglected in the estimation of transformation parameters.
To summarize: (1) When unit quaternion is introduced to represent spatial rotation in point cloud registration, the accuracy of the translation parameters and the scale factor will surely be affected in the case of the introduction of errors in the estimation of rotation parameters; (2) Planar feature is a good alternative in point cloud registration, however, the determination of the difference between each corresponding planar feature is difficult with traditional mathematical expressions; (3) Few methods incorporate the scale parameters in the planar feature-based registration, which makes it difficult to apply them in the fusion of data from different sources; (4) Most available methods estimate the registration parameters by iterative calculation, the initial approximate estimates of the unknown parameters must be determined in advance to assist the linearization. Based on the above analysis, our goal was to develop a closed-form solution to a planar feature-based registration algorithm in which a quad tuple-based representation of planar features in three-dimensional space is given, which made it possible to directly determine the difference between each pair of corresponding planar features; dual quaternions were employed to represent the spatial transformation in three-dimensional space which made it possible to estimate all the registration parameters simultaneously. The scale factor was also considered in our registration methods, which made it possible to apply our method in the fusion of data from different sources.

Representation of a Plane in Three-Dimensional Space
Traditionally, a plane is represented by the normal vector → n and any one point → p lying on it. Since the selection of the point → p is random, there is more than one representation for a single plane, which makes it difficult to compare and determine the difference between two planes. To ensure the unique representation of a plane in three-dimensional space, a quad tuple-based representation method was employed, as shown in Figure 3, in which a plane can be represented as a quad tupleΓ, as shown in Equation (1): where → l represents the normalized normal vector that forms n , and m represents the plane moment, that is, the distance between the origin and the plane that forms fusion of data from different sources; (4) Most available methods estimate the registration parameters by iterative calculation, the initial approximate estimates of the unknown parameters must be determined in advance to assist the linearization. Based on the above analysis, our goal was to develop a closed-form solution to a planar feature-based registration algorithm in which a quad tuple-based representation of planar features in threedimensional space is given, which made it possible to directly determine the difference between each pair of corresponding planar features; dual quaternions were employed to represent the spatial transformation in three-dimensional space which made it possible to estimate all the registration parameters simultaneously. The scale factor was also considered in our registration methods, which made it possible to apply our method in the fusion of data from different sources.

Representation of a Plane in Three-Dimensional Space
Traditionally, a plane is represented by the normal vector n and any one point p lying on it. Since the selection of the point p is random, there is more than one representation for a single plane, which makes it difficult to compare and determine the difference between two planes. To ensure the unique representation of a plane in three-dimensional space, a quad tuple-based representation method was employed, as shown in

Definition of a Dual Quaternion
A quaternion is a composite of four real numbers, which is generally represented in the form of a quad tuple r , as shown in Equation (2): Figure 3. The quad tuple-based representation of a plane in three-dimensional space.

Definition of a Dual Quaternion
A quaternion is a composite of four real numbers, which is generally represented in the form of a quad tuple . r, as shown in Equation (2): r is usually called a unit quaternion. A dual quaternion is a composite of two quaternions and a dual number ε, as shown in Equation (3) where . r and .
s are both quaternions, each of which corresponds to the real part and the dual part ofq separately; when . r T . r = 1,q is a unit dual quaternion, which is the default definition of a dual quaternion.

Operation Rules for Dual Quaternions
Given two dual quaternionsq 1 = .
s 2 , the addition, subtraction, and multiplication between them can be expressed as shown in Equation (4): Furthermore, multiplication between . r 1 and . r 2 can also be represented as shown in Equation (5):

Dual Quaternion-Based Transformation of Planar Features
According to the definition of a dual quaternion, the spatial transformation of a plane in three-dimensional space can be expressed as shown in Equation (6): whereq is the dual quaternion corresponding to the spatial transformation;q * is the conjugate ofq that formsq * = . r * + ε . s * ;Γ a andΓ b represent a pair of conjugate planar features, Γ a represents the one before transformation andΓ b represents the one after transformation.
To realize the dual quaternion-based similarity transformation, a plane should be represented as shown in Equation (7) According to Chasles' theorem [29], when the scale parameter is not considered, six independent parameters are needed to represent the transformation between two different coordinate frames. However, when a dual quaternion is used, there are eight parameters in total. Two constraints should be added to estimate the unknown parameters for the spatial similarity transformation, as shown in Equation (8): r is the quaternion that represents the spatial rotation and .
s is the quaternion that represents the spatial translation. The explanation in Equation (8) is that . r is a unit quaternion and it is orthogonal to . s.

Construction of the Objective Functions
In three-dimensional space, the transformation of a plane can be expressed as shown in Equation (9): where → l b and m b represent the normal vector and the moment of the plane before transformation, respectively; → l a and m a represent the normal vector and the moment after transformation, respectively; R, → t , and µ separately represent the transformation matrix, the translation vector, and the scale factor between the two coordinate systems.
With a dual quaternion, Equation (9) can be further expressed as Equation (10): where .
r is the unit quaternion corresponding to the rotation matrix R, r, Equation (10) can be rewritten as Equation (11): where . s * is the conjugate of . s. Based on Equation (7), Equation (11) can be further rewritten as Equation (12): Considering the existence of random errors, the planar feature-based registration approach aims to minimize the difference between  (13) and (14): The transformation parameters can be obtained when the expression f = f 1 2 + f 2 2 reaches its minimum. Considering that f 1 2 and f 2 2 are both positive, when each term reaches its minimum, f will be minimal. The optimal value of . r will be obtained by the minimization of Equation (13); . s and µ will be obtained by the minimization of Equation (14). (13) can be decomposed as Equation (15):

Solution of Rotation Quaternions
Using Equation (5) as the restriction, the error function can be expressed as Equation (16): where λ 1 is a Lagrange multiplier constant.
Taking the partial derivative of Equation (16) yields Equation (17): Making A = − 1 2 C l + C T l , Equation (17) can be further expressed as Equation (18): According to the definition of the eigenvalues and eigenvectors of a matrix, .
r represents one of the four eigenvectors of A, and λ 1 represents the corresponding eigenvalue of . r. Among the four eigenvectors that satisfy Equation (18), the optimal solution can be determined by referring back to Equation (17).
Multiplying Equation (17) with . r T yields Equation (19): Substituting Equation (19) into Equation (16) yields Equation (20): When λ 1 reaches its maximum, Equation (20) will be minimized. The conclusion can be drawn that the optimal solution of . r equals to the eigenvector corresponding to maximum eigenvalue of A.

Solution of the Translation Quaternions
l b , and decomposing Equation (14), we can obtain Equation (21): Using Equation (5) as the restriction, the best quaternion . s to represent the translation can be obtained when Equation (22) is minimized: where λ 2 is a Lagrange multiplier constant. Taking the partial derivative of Equation (22), with respect to . s and µ, and making them equal to zero, we obtain Equations (23) and (24): Based on Equation (24), the scale factor µ can be expressed as Equation (25): Substituting Equation (25) into Equation (23) yields Equation (26): Multiplying Equation (26) with . r T , we obtain Equation (27): With Equation (27), λ 2 can be expressed as Equation (28): Since C m2 T + C m3 T is a skew-symmetric matrix, the first item of Equation (28) will be zero, that is, T . r = 0. Equation (28) can be simplified, as shown in Equation (29): The quaternion corresponding to the translation between the two neighboring LiDAR stations can be obtained using Equation (30):

Implementation of the Proposed Solution
As in Figure 4, the given point clouds from the two neighboring LiDAR stations, namely the reference station and the unregistered station, suppose that (1) Based on the normal vectors of each pair of conjugate planar features, construct matrix A and calculate the maximum eigenvalue and its corresponding eigenvector r  of A; (2) Calculate the quaternion s  using Equation (26); (3) Calculate the scale factor  using Equation (25); (4) Calculate the quaternion t using Equation (30). Using r , s , and  , a plane in three-dimensional space can be transformed from one coordinate system to another. More importantly, since each pair of conjugate planar features is extracted from the two neighboring LiDAR stations, point clouds will be merged using r , s , and  .

Experiments and Results
The proposed planar feature-based registration algorithm was implemented using Matlab. Two experiments were designed to verify the correctness and the feasibility in which simulated data and real data were both incorporated.

Simulated Data
The first experiment was designed to verify the correctness of the proposed algorithm. Five pairs of simulated planar features were designed as shown in Figure 5. Mathematical descriptions of the simulated planar features were as shown in Table 2. The preestablished spatial similarity transformation parameters were as shown in Table 3.

Solution of Unknown Parameters
Least Squares Criteria Dual Quaternion s, and µ, a plane in three-dimensional space can be transformed from one coordinate system to another. More importantly, since each pair of conjugate planar features is extracted from the two neighboring LiDAR stations, point clouds will be merged using . r, . s, and µ.

Experiments and Results
The proposed planar feature-based registration algorithm was implemented using Matlab. Two experiments were designed to verify the correctness and the feasibility in which simulated data and real data were both incorporated.

Simulated Data
The first experiment was designed to verify the correctness of the proposed algorithm. Five pairs of simulated planar features were designed as shown in Figure 5. Mathematical descriptions of the simulated planar features were as shown in Table 2. The pre-established spatial similarity transformation parameters were as shown in Table 3.     Table 3. Pre-established spatial similarity transformation parameters. The obtained transformation parameters were as shown in Table 4; the residuals between each pair of conjugate planar features after registration are also given. Based on the results given in Table 4, the calculated rotation matrix and the scale factor were both consistent with the pre-established ones. The small deviation between the calculated translation vector and the pre-established one can be attributed to the rounding errors in the calculation, which can largely be ignored in practical applications. The conclusion can be drawn that the proposed solution is correct and the obtained result is as expected.

Real Data
The second experiment was designed to verify the feasibility of the proposed algorithm. The point clouds were collected by a model Riegl LMS-Z420i terrestrial laser scanner. In order to collect the complete point cloud data of the building, a total of eight observation stations were set up around it. The average distance between each observation station and the building was about 100 m, and the average sampling interval was set to 4 cm. Seven pairs of conjugate planar features were extracted from the two neighboring point clouds as shown in Figure 6, and the extracted planar features were as shown in Table 5.
rithm. The point clouds were collected by a model Riegl LMS-Z420i terrestrial laser scanner. In order to collect the complete point cloud data of the building, a total of eight observation stations were set up around it. The average distance between each observation station and the building was about 100 m, and the average sampling interval was set to 4 cm. Seven pairs of conjugate planar features were extracted from the two neighboring point clouds as shown in Figure 6, and the extracted planar features were as shown in Table 5.  The obtained transformation parameters and residuals between each pair of conjugate planar features are both given in Table 6. To further verify the correctness and the feasibility of the proposed solution, a unit quaternion-based, iterative method [8] was employed to estimate the transformation parameters; the results are also shown in Table 6. The obtained transformation parameters and residuals between each pair of conjugate planar features are both given in Table 6. To further verify the correctness and the feasibility of the proposed solution, a unit quaternion-based, iterative method [8] was employed to estimate the transformation parameters; the results are also shown in Table 6. The visual effects of the point clouds from the two neighboring LiDAR stations before and after registration were as shown in Figure 5.
As shown in Table 6 and Figure 7, by using the obtained transformation parameters to register the two neighboring stations, there were no significant residuals between each pair of conjugate planar features after registration. More importantly, there were also no significant differences between the results obtained by our method and by the unit quaternion-based iterative method [8]; the root mean square error (RMSE) of the normal and moment's differences between conjugate planar features were exactly the same. The conclusion can be drawn that the newly proposed closed-form solution to pairwise registration of LiDAR point clouds is correct and feasible. The visual effects of the point clouds from the two neighboring LiDAR stations before and after registration were as shown in Figure 5.
As shown in Table 6 and Figure 7, by using the obtained transformation parameters to register the two neighboring stations, there were no significant residuals between each pair of conjugate planar features after registration. More importantly, there were also no significant differences between the results obtained by our method and by the unit quaternion-based iterative method [8]; the root mean square error (RMSE) of the normal and moment's differences between conjugate planar features were exactly the same. The conclusion can be drawn that the newly proposed closed-form solution to pairwise registration of LiDAR point clouds is correct and feasible.

Discussion
On the assumption that the normal directions of each pair of conjugate planar features are the same, the proposed quad tuple-based representation method can provide a unique mathematical expression of any planar feature in three-dimensional space; the differences between two planar features can be determined by direct comparison, which makes it convenient to propose and implement a closed-form solution to planar featurebased point cloud registration. The results of the two designed experiments were both correct, which proves that the proposed closed-form solution is capable of dealing with point cloud registration problems on the condition that sufficient corresponding planar features are provided as registration primitives.

Discussion
On the assumption that the normal directions of each pair of conjugate planar features are the same, the proposed quad tuple-based representation method can provide a unique mathematical expression of any planar feature in three-dimensional space; the differences between two planar features can be determined by direct comparison, which makes it convenient to propose and implement a closed-form solution to planar feature-based point cloud registration. The results of the two designed experiments were both correct, which proves that the proposed closed-form solution is capable of dealing with point cloud registration problems on the condition that sufficient corresponding planar features are provided as registration primitives.
It should be noted that any vector in three-dimensional space can be expressed in another form with the same module and an opposite direction; however, the prerequisite of our solution is that normal directions of each pair of conjugate planar features should be exactly the same after registration. On the condition that the prerequisite is not fulfilled, ensuring the proper run of the proposed solution will be one of our objectives in the future.

Conclusions
On the condition that point and linear features are not sufficient, planar features are good alternatives to ensure the proper run of the point cloud registration algorithm. Compared to point and linear features, extracting planar features from point clouds is more convenient. Furthermore, the impact of random errors can be significantly reduced by least squares fitting. Based on the quad tuple-based representation method of planar features in three-dimensional space, a closed-form solution to planar feature-based registration of point clouds is proposed. Two experiments were designed, in which both simulated data and real data were used to verify the correctness and effectiveness of the proposed solution. With the simulated data, the calculated registration results were consistent with the preestablished parameters and with the real data, there were also no significant differences between the results obtained by our method and by the available iterative method [8]; the root mean square error (RMSE) of the normal and moment's differences between conjugate planar features were the same.
The conclusion can be drawn that the newly proposed closed-form solution to pairwise registration of LiDAR point clouds is correct and feasible. Moreover, the advantages of the proposed solution can be summarized as follows: (1) By using eigenvalue decomposition to replace the linearization of the objective function, the presented solution does not require any initial estimates of unknown transformation parameters in advance, which assures the stability and robustness of the algorithm; (2) On the condition that no man-made reflectors are used, planar features are directly extracted from point clouds using least squares fitting, the impact of random errors can be reduced significantly and the reliability of the estimated transformation parameters is higher than point and linear feature-based methods; (3) In contrast to Euler angle-based and rotation matrix-based methods, using dual quaternions to represent spatial rotation greatly reduces additional constraints in the estimation of similarity transformation parameters.