Global Registration of Terrestrial Laser Scanner Point Clouds Using Plane-to-Plane Correspondences

: Registration of point clouds is a central problem in many mapping and monitoring applications, such as outdoor and indoor mapping, high-speed railway track inspection, heritage documentation, building information modeling, and others. However, ensuring the global consistency of the registration is still a challenging task when there are multiple point clouds because the di ﬀ erent scans should be transformed into a common coordinate frame. The aim of this paper is the registration of multiple terrestrial laser scanner point clouds. We present a plane-based matching algorithm to ﬁnd plane-to-plane correspondences using a new parametrization based on complex numbers. The multiplication of complex numbers is based on analysis of the quadrants to avoid the ambiguity in the calculation of the rotation angle formed between normal vectors of adjacent planes. As a matching step may contain several matrix operations, our strategy is applied to reduce the number of mathematical operations. We also design a novel method for global reﬁnement of terrestrial laser scanner data based on plane-to-plane correspondences. The rotation parameters are globally reﬁned using operations of quaternion multiplication, while the translation parameters are reﬁned using the parameters of planes. The global reﬁnement is done non-iteratively. The experimental results show that the proposed plane-based matching algorithm e ﬃ ciently ﬁnds plane correspondences in partial overlapping scans providing approximate values for the global registration, and indicate that an accuracy better than 8 cm can be achieved by using our global ﬁne plane-to-plane registration method.


Introduction
In recent years, the rising role of building information modeling (BIM) has driven the market for static terrestrial laser scanners (TLS). TLS is an efficient and cost-effective technology for rapid and accurate collection of 3D point clouds. TLS point clouds are the most suitable data source for the generation of as-built BIMs, which are a crucial tool for many construction and architecture professionals. To derive globally consistent 3D point cloud models from TLS with high positional accuracy, registration is a mandatory task. Typically, in the existing frameworks, the registration of multiple point clouds is divided into a pairwise registration step and a global fine registration step. Pairwise registration involves finding feature correspondences between pairs of point clouds and minimizing the sum of residuals over all such feature correspondences for the estimation of transformation parameters (3D rotation matrix and 3D translation vector), which establish the relative orientation for each pair of scans in a common coordinate system. In practice, pairwise registration using free-form correspondences (e.g., iterative closest point algorithm [1]), feature point-based (e.g., keypoints) methods [2][3][4][5][6][7], or primitive-based (e.g., lines or planar surfaces) approaches [8][9][10][11][12][13] should be applied first to obtain the transformation parameters. However, a problem which arises in the registration of multiple point clouds is that corresponding scan features may still present significant residual errors from pairwise registration task. As a direct consequence, most of the pairs of registered point clouds will contain loop closing errors, thus reducing their local metric accuracy. Then, the global refinement should be used to minimize the influence of pairwise alignment errors on the transformation parameters, by distributing the residual errors evenly across the scan project [14]. Our present contribution is twofold. First, a plane matching strategy by decoupling rotation and translation parameters is proposed. The proposed matching algorithm makes use of a new parametrization based on complex numbers to find the correspondence between pairs of segmented planes. Second, a new global closed-form solution is proposed via graph-based formulation adapted for plane-to-plane correspondences. The proposed solution refines the sensor positions by treating the 3D points and their corresponding surface normal as observations. To the best of our knowledge, the proposed method is the first using this approach.
The rest of the paper is organized as follows. Section 2 provides the related works in point cloud registration, helping the reader to gain an insight into the registration problem. The studied area, our proposed plane-based matching approach, and our global fine registration solution based on plane correspondences are presented in Section 3. Experimental evaluation of the proposed solution on real datasets and a discussion about the potential and limitation of the proposed method are given in Section 4. Finally, the paper is concluded in Section 5.

Related Work
Typically, standard global refinement solutions first find the transformation parameters using [1], then evenly redistribute loop closing errors with the help of a graph-based optimization such as the one proposed by [15], which has been extended to 3D by [16]. This task is formulated as a scan graph, in which each scan denotes a node, and each edge highlights the spatial constraint between the pairs of nodes. A globally consistent registration of multiple point clouds via graph optimization was proposed by [14]. The authors find the coarse alignments using an unambiguous geometric keypoint configuration able to reduce the number of candidates that are sampled in the iteration process, called on K-4PCS [6], and the pairwise candidate solutions are filtered using a discrete graphical model of the scan network. Then, least-squares optimization [15] is used to evenly distribute the residual error. Ji et al. [13] generated a globally consistent 3D map of high-speed viaduct point clouds using closing condition and external geometric constraints. Their proposed method uses artificial targets for the iterative pairwise registration step. The pairwise transformations are also refined with [1], followed by global refinement to evenly spread the residual errors, as described in [15]. Yang et al. [17] use a branch-and-bound strategy to globally solve the objective function. It first uses the iterative closest point algorithm for a coarse registration task and refines the transformation parameters with the algorithm of [15]. The approach proposed by Huber and Hebert [2] uses surfaces and a Bayesian filter for alignment of multi-view 3D point clouds. The main disadvantage of this method is its high computational cost. Other global refinement approaches rely on general graph optimization [18], bundle adjustment [19], low-rank sparse decomposition strategy [20], kernel-based energy function [21] and visual bag of words [22], which are computationally less attractive.
The aforementioned existing methods for registration of multiple point clouds are prominently iterative free-form or point-based correspondences. Despite their popularity, iterative solutions have several limitations. For instance, they are not effective for sparse point clouds, require good approximate values to avoid convergence to weak local minima, they are sensitive to noise, and involve a matching step which incurs a high computational cost. Furthermore, when a downsampling step is applied, as in [6], the details of objects in the scene may be lost [12]. In contrast, closed-form solutions estimate the transformation parameters in one-step and they do not require approximate values.
Additionally, plane-based approaches can achieve alignments with higher accuracy, are less influenced by the presence of outliers [23], and are more robust in identifying corresponding feature pairs [24]. In [25] a closed-form solution for pose-graph relaxation is introduced for enhancing the consistency of 3D maps. First, planar surfaces are extracted from the point clouds and matched. In the subsequent step, the pairwise registration of scans is performed. Assuming that the rotation parameters are accurately estimated, the sensor pose vectors are then refined using the least-squares solution of [15]. Pavan and dos Santos [26] introduce a global closed-form refinement of multiple point clouds with local consistency of planes. They use the similarity of plane properties and the geometric constraint formed by planar surfaces to identify correspondences and place all the rotation parameters into a common coordinate system by exploiting a property of quaternions. Moreover, the global refinement procedure is done using a combination of [1] and the global refinement algorithm of [15]. Yan et al. [27] presented a framework for global registration of building point clouds using portals (windows and doors) of rooms to connect scans with limited overlap, while the global refinement step is formulated based on linear integer programming. First, they classify all the points into horizontal and vertical groups according to the normals of the points. Then, portals are extracted by detecting boundary points in the vertical segments. After, pairwise matching are encoded in a connection graph. The globally consistent registration of the point clouds is obtained by choosing an optimal subset of the pairwise matchings from the graph. To avoid enumerate an extremely large number of constraints, an acceleration scheme by iterative adding constraints is proposed. During the registration step, local alignments remain conflicting after each optimization step. Thus, the set of constraints are added and the optimization is realized again. This process is repeated until the conflict free registration is obtained. The registration is refined using [1].
In general, closed-form solutions are more efficient because they provide the best transformation in one-step, without requiring initial approximations for each pair of scans. Closed-form solutions avoid convergence issues of iterative methods and have been shown to achieve more efficient and accurate registration results [28]. Since urban environments are characterized by an abundance of planar surfaces and current point cloud registration methods are inefficient for large-scale data [29], a closed-form global registration method based on planes have clear advantages for urban mapping applications [28,30]. While closed-form solutions are more suitable for pairwise registration of point clouds, there is a lack of global closed-form solutions in the literature, especially for plane-based approaches. In this article, we propose a new closed-form solution for global registration of point clouds from plane-to-plane correspondences.

Studied Area
We evaluate our proposed method on three different TLS datasets. We use one outdoor and two indoor scan projects for our experiments. Detailed information of the experiments is listed in Table 1. The Faro 880 TLS offers a near-spherical field of view made possible by a 320 • vertical angle scanning range and a 180 • horizontal field with a linearity error in the rangefinder less than 3 mm at a range of 10 m, a ranging error about 1 mm and the measurement speed of 122,000 pts/sec. The Faro Focus S120 offers a distance accuracy up to ± 2 mm, a range from 0.6 m up to 120 m, a measurement rate up to 976,000 pts/sec, vertical field of view (vertical/horizontal) of 305 • /360 • , and a maximal vertical scan speed of 97 Hz. Figure 1 shows the marker position for each dataset.
Remote Sens. 2020, 11, x FOR PEER REVIEW 19 of 19 a range of 10 m, a ranging error about 1 mm and the measurement speed of 122,000 pts/sec. The Faro Focus S120 offers a distance accuracy up to ± 2 mm, a range from 0.6 m up to 120 m, a measurement rate up to 976,000 pts/sec, vertical field of view (vertical/horizontal) of 305°/360°, and a maximal vertical scan speed of 97 Hz. Figure 1 shows the marker position for each dataset. Our outdoor scan project "Patio Batel" (located in a urban area in the south of Brazil, -30°04′09"S, -51°24′20"E) is part of a densely urban environment, in which there are cars, pedestrians, trees, poles, streets and building façades, as shown in Figure 1a. The "Lape" (located in a urban area in the south of Brazil, -30°04′09"S, -51°24′20"E) and the "Royal exhibition building" " (located in an urban scene in the south of Australia, -38°34′06"S, 145°18′07"E) experiments were conducted to evaluate the performance of the method in an indoor environment. The first dataset contains features as tables, computers, chairs and façade constructions, as shown in Figure 1b. The "Royal exhibition building" experiment was carried out to evaluate the performance of the method with scans of a historical building with complex interior structures. The indoor area features as chairs and façade constructions, as shown in Figure 1c. Figure 2 shows some detailed views of unregistered point clouds for each datasets. Our outdoor scan project "Patio Batel" (located in a urban area in the south of Brazil, -30 • 04 09"S, -51 • 24 20"E) is part of a densely urban environment, in which there are cars, pedestrians, trees, poles, streets and building façades, as shown in Figure 1a. The "Lape" (located in a urban area in the south of Brazil, -30 • 04 09"S, -51 • 24 20"E) and the "Royal exhibition building" " (located in an urban scene in the south of Australia, -38 • 34 06"S, 145 • 18 07"E) experiments were conducted to evaluate the performance of the method in an indoor environment. The first dataset contains features as tables, computers, chairs and façade constructions, as shown in Figure 1b. The "Royal exhibition building" experiment was carried out to evaluate the performance of the method with scans of a historical building with complex interior structures. The indoor area features as chairs and façade constructions, as shown in Figure 1c. Figure 2 shows some detailed views of unregistered point clouds for each datasets.

Method
In this section, we first describe the general prerequisites for TLS scans acquisition and preprocessing with a brief description of plane segmentation using the random sample consensus (RANSAC) algorithm [31]. Then, we present the closed-form registration method that first finds a set of plane correspondences with complex numbers and minimizes the sum of residuals over all such plane correspondences. Finally, the closed-form solution, which uses a set of plane-to-plane correspondences to consistently redistribute the residual errors across the whole project, is presented. In Figure 3, the proposed framework is shown, and the involved steps are illustrated. The details of these steps are explained in the subsequent sections.

Method
In this section, we first describe the general prerequisites for TLS scans acquisition and preprocessing with a brief description of plane segmentation using the random sample consensus (RANSAC) algorithm [31]. Then, we present the closed-form registration method that first finds a set of plane correspondences with complex numbers and minimizes the sum of residuals over all such plane correspondences. Finally, the closed-form solution, which uses a set of plane-to-plane correspondences to consistently redistribute the residual errors across the whole project, is presented. In Figure 3, the proposed framework is shown, and the involved steps are illustrated. The details of these steps are explained in the subsequent sections.

Plane Segmentation with the RANSAC Algorithm
Planar surfaces are abundant in human-made environments, and their extraction requires less storage than voxel-grid representations. However, due to the varying point density of TLS data, measurement errors can lead to sparse outliers. These can influence the plane extraction task, and reduce the registration accuracy. To avoid these weaknesses, we start by removing outliers from the raw point cloud (ℵ) using the statistical outlier removal technique proposed by [32]. The filtered point cloud ℵ is obtained based on a statistical analysis of a query point with respect to its surrounding neighbors ∈ ℵ. Thus, given a filtered reference point cloud ℵ and a filtered target point cloud ℵ with both inliers and outliers, the RANSAC algorithm is used for plane extraction task. This is done by randomly selecting three non-collinear points in ℵ , followed by a voting scheme to find the optimal fitting plane model. The best fit plane is estimated by finding the eigenvector corresponding to the smallest eigenvalue of the covariance of weighted point coordinates, as described in [17].

Estimation of Transformation Parameters from Plane Correspondences
Typically, pairwise registration of point clouds follows a two-step strategy. First, corresponding features are required in overlapping scans, and second, transformation parameters are estimated to align a pair of point clouds relative to each other in a common coordinate system. Our proposed strategy is based on plane correspondences and estimates the rotation and translation parameters separately. Given two sets of segmented planar surfaces ⊂ ℵ and ′ ⊂ ℵ , the transformation from ∈ ℵ to ′ ∈ ℵ consists of a 3D rigid motion, as follows:

Plane Segmentation with the RANSAC Algorithm
Planar surfaces are abundant in human-made environments, and their extraction requires less storage than voxel-grid representations. However, due to the varying point density of TLS data, measurement errors can lead to sparse outliers. These can influence the plane extraction task, and reduce the registration accuracy. To avoid these weaknesses, we start by removing outliers from the raw point cloud (ℵ) using the statistical outlier removal technique proposed by [32]. The filtered point cloud ℵ is obtained based on a statistical analysis of a query point p i with respect to its surrounding neighbors k ∈ ℵ. Thus, given a filtered reference point cloud ℵr and a filtered target point cloud ℵt with both inliers and outliers, the RANSAC algorithm is used for plane extraction task. This is done by randomly selecting three non-collinear points in ℵ , followed by a voting scheme to find the optimal fitting plane model. The best fit plane is estimated by finding the eigenvector corresponding to the smallest eigenvalue of the covariance of weighted point coordinates, as described in [17].

Estimation of Transformation Parameters from Plane Correspondences
Typically, pairwise registration of point clouds follows a two-step strategy. First, corresponding features are required in overlapping scans, and second, transformation parameters are estimated to align a pair of point clouds relative to each other in a common coordinate system. Our proposed strategy is based on plane correspondences and estimates the rotation and translation parameters separately. Given two sets of segmented planar surfaces S i ⊂ ℵr and S i ⊂ ℵt , the transformation from p i ∈ ℵr to p i ∈ ℵt consists of a 3D rigid motion, as follows: where R denotes the 3D rotation, t represents a 3D translation and The condition that a point p i lies on a plane π = u T , d T in ℵr , where u = [a , b , c ] T is the normal vector of plane π and d its distance to the origin, can be written as [26]: where u is the unit normal vector of the plane π and d its distance to the origin. By substituting Equation (1) into (2), we can obtain the following expression: from Equation (3), we have: note that the rotation R from vector u to the vector u is expressed as: We now have two separate sets of equations per plane correspondence, one for t, as presented in Equation (4), and one for R, as can be seen in Equation (5). Given a set of putative correspondences, we can estimate the transformation parameters R and t using a closed-form solution. As described by [26], a system of equations At = l + v can be formed to estimate the parameter t, where A and l are obtained by stacking u and d − d for all plane-to-plane correspondences and v contains the residual values. The least-squares solution for t can be obtained as: The estimate of R is obtained using Horn's solution [33]. In the subsequent section, we describe our proposed method to find corresponding planes using complex numbers.

Plane Matching with Complex Numbers
In its basic form, the estimation of the transformation parameters relies on the presence of pairs of corresponding planes. Dold and Brenner [34] give a plane matching method for TLS point clouds. They investigated the combinatorial complexity of the search for corresponding planes. Our plane matching algorithm extends this method by adding a false positive detection task. This is done by inferring the relative angle between patch plane vectors using complex numbers. It works in three steps: (1) Classification of the segmented planes; (2) Determination of the rotation between pairs of candidate planes; and (3) Calculation of approximate translation.
• Classifying the segmented planes First, we assume that the TLS is mounted on a tripod and is levelled using a bubble level. Consequently, planes extracted from TLS data have small deviations from vertical and horizontal planes [35]. Thus, given two sets of segmented planar surfaces S i ⊂ ℵr and S i ⊂ ℵt , our algorithm separates these in different classes, i.e., horizontal planes (representing the planes parallel to the floor) and vertical planes (representing the wall planes of buildings) using the third component of the normal vector of S i or S i . The horizontal planes are those whose normal vectors and the z axis form an angle smaller than 3 degrees, while the vertical planes are those whose normal vectors and the z axis form an angle larger than 89 degrees. Such classification procedure is essential to reduce the search space, and for the optimization of the matching performance.
• Calculating the rotation between normal vectors of pairs of candidate planes Second, given a set of segmented planes in both the reference point cloud (ℵr ) and the target point cloud (ℵt ), as presented in Figure 4a, in order to calculate the rotation angle (θ) between the normal vectors (Figure 4b) of each pair of candidate planes for correspondence (Figure 4c), we propose the use of complex numbers. Figure 4 shows an example of segmented planes, the normal vectors and the relative roation angle between two segmented planes in ℵr and ℵt . Remote Sens. 2020, 11, x FOR PEER REVIEW 19 of 19 As shown in Figure 4c, and represent the rotation angles of u to u ′ and u to u ′, respectively, where u and u ′ are the normals of the pair of segmented planes in ℵ and ′ in ℵ ′, and u and u ′ are the normals of the pair of segmented planes in ℵ and ′ in ℵ ′.
For vertical planes, we can assume that is only applied around the axis; then, for every pair of adjacent non-parallel vertical planes, we have: Equation (6) can be rewritten as: where S( ) is the rotation matrix in the plane. For all angles , the matrices S( ) form a special orthogonal group (2). All rotation S( ) in the plane can be represented by a complex number. Then, the rotation in the plane of a vector ( , ) ∈ ℝ can be rewritten by the multiplication of complex numbers, as follows [36]: where = cos + sin . Thus, the complex number can be calculated by multiplying Equation (8) by the inverse of the complex number ( + ), as follows: Then, the rotation angle between ( , ) and ( , ′) can be calculated using the following expression: Subsequently, the algorithm checks the correspondence for every two pairs of non-parallel vertical planes calculating the difference between their relative angle ( ), as follows: The pseudo-correspondences between the two pairs of non-parallel vertical planes are established if | | ≤ 1 degree, otherwise, the correspondence hypothesis is rejected. Please note that, to avoid the singularity problem, the planes must be non-parallel vertical planes, since vertical planes are influenced by the rotation of the TLS sensor.
• Selection of pairwise plane correspondences In the third step, the components and are calculated using the pairs of corresponding non-parallel vertical planes. For this step, we can reduce the search space by initially assuming that As shown in Figure 4c, θ 1 and θ 2 represent the rotation angles of u 1 to u 1 and u 2 to u 2 , respectively, where u 1 and u 1 are the normals of the pair of segmented planes π 1 in ℵr and π 1 in ℵr , and u 2 and u 2 are the normals of the pair of segmented planes π 2 in ℵt and π 2 in ℵt . For vertical planes, we can assume that R is only applied around the z axis; then, for every pair of adjacent non-parallel vertical planes, we have: Equation (6) can be rewritten as: where S(θ) is the rotation matrix in the xy plane. For all angles θ, the matrices S(θ) form a special orthogonal group M(2). All rotation S(θ) in the plane can be represented by a complex number. Then, the rotation in the plane of a vector (a, b) ∈ R 2 can be rewritten by the multiplication of complex numbers, as follows [36]: where e iθ = cos θ + i sin θ. Thus, the complex number e iθ can be calculated by multiplying Equation (8) by the inverse of the complex number (a + ib), as follows: Then, the rotation angle between (a, b) and (a , b ) can be calculated using the following expression: Subsequently, the algorithm checks the correspondence for every two pairs of non-parallel vertical planes calculating the difference between their relative angle (ε), as follows: The pseudo-correspondences between the two pairs of non-parallel vertical planes are established if |ε| ≤ 1 degree, otherwise, the correspondence hypothesis is rejected. Please note that, to avoid the singularity problem, the planes must be non-parallel vertical planes, since vertical planes are influenced by the rotation of the TLS sensor.
• Selection of pairwise plane correspondences In the third step, the components t x and t y are calculated using the pairs of corresponding non-parallel vertical planes. For this step, we can reduce the search space by initially assuming that t z = 0. Thus, we can use two corresponding planes to estimate an approximation of t x and t y , as follows: which can be further decomposed to: where The solution for Equation (13) is given as follows: An issue which requires attention in this solution is the correspondence between non-parallel vertical planes. The choice of parallel vertical planes can lead to an inconsistent system of equations with no solution. Similarly, the solution for component t z is obtained using the third pair of corresponding non-parallel vertical planes, as follows: where u 3 = 0 0 1 T and u 3 = 0 0 1 T are parallels to the z axis. Next, all pairs of vertical planes considered as pseudo-corresponding are combined to find the correct correspondences. This is done by using the approximated translation values and the rotation angle, as follows: where β is the rotation angle from normal vector u in S i to the combined normal vector u in S i , d and d are the perpendicular distances between the origin of the coordinate system and the combined planes, respectively.
For δ values less than a predefined threshold and for the smallest error xy value, the pair of corresponding planes is added to the set of pseudo-correspondences with respect to the planes π 1 ↔ π 1 and π 2 ↔ π 2 . Then, the two steps of rotation angle estimation and the approximated estimation of t x and t y components are repeated using other planes and a new set of pseudo-corresponding planes is found. The combinations inserted in this set are regarded as true correspondences. In addition, the third component of translation parameters (t z ) is recalculated using corresponding horizontal planes. Thus, pairs of pseudo-corresponding horizontal planes are combined to find the correct correspondences using the following expression: When error z is less than a predefined threshold, the pair of corresponding horizontal planes is added to the set of pseudo-correspondences with respect to the planes π 3 ↔ π 3 . The combinations inserted in this set are also regarded as true correspondences. After all verifications, the potential corresponding planes are used to refine the translation parameters.

Proposed Global Plane-to-Plane Refinement Solution
Due to the accumulated residual errors resulting from the consecutive pairwise registration steps, a global refinement task should be used to minimize the loop closure error across all scan pairs in the project [14]. Given a discrete set of putative pairwise registrations, we develop a global refinement solution based on a plane-to-plane approach. The global refinement is done with the help of a pose-graph structure.
According to [37], a graph G(X, T) is formed by nodes (X) and edges (T), where nodes represent the pose of each point cloud and edges denote each pairwise registration (transformation of pairs points cloud) with sufficient mutual overlap. An example of a graph structure is shown in Figure 5. The variables X i ∈ X denote point clouds and the each pairwise registration T ij ∈ T are represented by transformation parameters M ij (rotation) and t ij (translation) between the point clouds X i and X j , respectively, with j, i ∈ {0, 1, 2, . . . , n}. The variables ∈ denote point clouds and the each pairwise registration ∈ are represented by transformation parameters (rotation) and (translation) between the point clouds and , respectively, with , ∈ 0,1, 2, … , . • Global refinement of rotation parameters The rotation parameters are globally refined ( ) using our previous work [26], as follows: where represents the 3D rotation matrix from node to the node (reference) and denotes the 3D rotation matrix from to with and = 0, 1, 2,…, (see Figure 4). Using quaternions, Equation (18) can be rewritten as follows: where represents the 3D rotation matrix from node to the node (reference) and M denotes the 3D rotation matrix from to and Δ is the residual error, with , ∈ 0,1, 2, … , . Thus, Equation (19) is rewritten as follows: where is a 4×4 matrix associated with left product of the quaternion , , and Ε are 4×1 vectors that represent the quaternions , and , respectively. Equation (20) can be rewritten as follows: • Global refinement of rotation parameters The rotation parameters are globally refined (R j ) using our previous work [26], as follows: where R i represents the 3D rotation matrix from node X i to the node X 0 (reference) and M ij denotes the 3D rotation matrix from X i to X j with j and i = 0, 1, 2,. . ., n (see Figure 4). Using quaternions, Equation (18) can be rewritten as follows: where R i represents the 3D rotation matrix from node X i to the node X 0 (reference) and M ij denotes the 3D rotation matrix from X i to X j and ∆ ij is the residual error, with j, i ∈ {0, 1, 2, . . . , n}. Thus, Equation (19) is rewritten as follows: where L q ij is a 4 × 4 matrix associated with left product of the quaternion . q ij , Q i , Q j and E ij are 4 × 1 vectors that represent the quaternionsˆ. q i ,ˆ. q j and ij , respectively. Equation (20) can be rewritten as follows: where Ω represents a 4n × 4m matrix formed by partial derivatives with respect to all components of Q j , Q i is the column vector 4n × 1 formed by concatenation of all matrices Q i and Q j , n is the number of the nodes and m denotes the number of the edges on the graph. Then, the total least square method [38] can be used to estimate the matrix Ψ (or R r ), as follows: Thus, the matrix Ψ that minimizes the Equation (22) is the eigenvector corresponding to the smallest eigenvalue of the matrix Ω T Ω. To obtainˆ. q i it is essential to normalize each Q i . As this solution is direct, it does not require initial approximations.
• Proposed global refinement of translation parameters By combining a pair of point clouds, a difference R r t i between two sensor poses x i and x j can be obtained, as follows: where R r represents the globally refined rotation parameters, t ij is the 3D translation vector from node X i to X j , x i is a column matrix that represents the 3D position of the node X i with x 0 = 0 0 0 T and x j is a column matrix that represents the 3D position of the node X j . As described, the translation parameters can be estimated from plane correspondences: where u T jk denotes the normal vector from segmented plane π jk ∈ X j and d jk its perpendicular distance from the origin, and d ik the distance from the origin to plane π ik ∈ X i . Multiplying Equation (24) by u T ijk R T r and substituting it into Equation (23) yields: For each pair of point clouds within the network, the proposed method calculates the residual errors based on plane-to-plane correspondences. As the Equation (25) is linear with respect to the sensor positions x 1 , x 2 ,. . ., x n and assuming that d ik − d jk is the vector of observations and u T jk is the normal vector for all plane-to-plane correspondences, respectively, the solution x i − x j can be globally refined using the least squares method (LSM). Thus, the set of liner equations Ax + v = l can be formed to estimate the sensor positions, where A is the coefficient matrix containing the partial derivatives in Equation (25), l is the vector of constants containing d − d for all plane-to-plane correspondences and v is the residual vector. The least-squares solution for x i − x j can be obtained as: Please note that the proposed method has been adapted to use plane parameters to estimate x 2 ,. . ., x n , instead of feature points as used by [15]. To obtain a solution, at least four pairs of corresponding planes are needed. Fortunately, this is not a major issue in urban environments, since a large number of planar surfaces are available.

Experiments and Results
First, the plane matching method is evaluated in terms of matching success rate (SR), by simply counting the number of correctly matched planes over all pairs of point clouds and the matching time (MT). Its performance also is compared to the K-4PCS registration algorithm [6], which is available in the point cloud library (PCL) [39]. Secondly, the global refinement method is evaluated in terms of mean registration accuracy. The residual root means square errors (RMSE) before and after global refinement are calculated based on the plane-to-plane distances between the corresponding planes. The performance of the proposed global refinement method is evaluated by comparison with three other approaches. The first is a combination of iterative closest point algorithm [1] and the global refinement method proposed by Lu and Milios [15]. The second approach is the global registration using our pipeline but omitting the global refinement of the rotation parameters. The third approach is the global registration proposed by Theiler et al. [14]. The first approach is available in the point cloud library (PCL) [39], while the source code from the third approach can be found in [40].

Pre-Processing Data
For all the datasets, the outliers were detected and removed using the statistical outlier removal algorithm. With the removal of the outliers, the number of points contained in the point clouds was reduced by 25%. Basically, this step of the method removes from the data sample all points outside the interval µ k ± α·σ k . In this paper, the values assumed for the variables µ k = 0.050 m (mean distance between neighbouring points) and α = 0.10 (restriction factor) were determined empirically and were those that best represented the expected sampling of the object on the surface. The surface planes were extracted with RANSAC algorithm using a distance threshold 0.01 m, and 100 iterations. Figure 6 shows an example of the segmented planes for pairs of point clouds obtained by using the RANSAC algorithm. As can be seen, ground planes could not be correctly classified because small inclinations exist due to the different ground levels in all scanned regions. However, this issue does not affect the performance of the matching algorithm.   The planes were classified and segmented into horizontal and vertical planes, with their normals and perpendicular distances calculated with respect to the origin. The error xy and error z threshold values were respectively 0.15 and 0.01 m for Patio batel and Lape datasets, while for the Royal exhibition building dataset, the value 0.01 m was used for error xy and error z .

Plane Matching Evaluation
As previously mentioned, we tested our plane-based matching algorithm against the K-4PCS algorithm proposed by [6]. For the tests conducted with K-4PCS, we used the following matching algorithm and parameter settings: 3D scale feature transform (SIFT) algorithm for keypoints detection, and voxel size τ = 50 mm, respectively. Table 2 shows the results obtained with both our plane-based matching algorithm and the K-4PCS method. For the "Patio batel" dataset, a sufficiently large number of corresponding planes was obtained for the pairs of point clouds containing salient structures. For these pairs of point clouds, the method achieves high SR (>80%). The SR is generally high for point clouds with large overlap (about 40%). The lower SR for a few of the pairs is due to the smaller overlap (30%) and the high degree of symmetry of planar surfaces in the scene, which reduces the chance for finding correct matches. Table 2 also shows the metric accuracy of the estimated transformation parameters represented by the root mean square error between true correspondences (RMSE R ). The RMSE R obtained with the K-4PCS algorithm is about 0.60 m, while our method achieved a RMSE R around 0.40 m. Probable reason for a large RMSE R value found is the lower overlap between the pairs of point clouds. Note also that both algorithms use plane correspondences, instead of targets. Thus, to improve the pairwise registration accuracy is necessary to apply the global refinement. For the "Lape" dataset, the proposed plane matching has been carried out with the optimal parameters, error xy and error z , set to 0.15 m and 0.01 m, respectively. From the test area, nine evenly distributed corresponding pairs of planes were found, which indicates a matching success rate above 75% (see Table 2). The reason for the moderate success rate is the high number of parallel planes discarded by the proposed plane matching algorithm to avoid shift errors. For the "Royal building exhibition" dataset, the optimal parameters are set to 0.01 m for both the planimetric and altimetric values. This setting is based on the overlap area between the point clouds: where for larger overlaps smaller values for the parameters are used. As can be observed, the reference and target point clouds are well registered. The mean SR for both the experiments is around 95%, as shown in Table 2. However, some failed cases with a large degree of scene symmetry occur. The results of the successful cases reveal that our proposed plane matching algorithm works well. For plane pairwise registration methods, the translational errors are larger than rotational errors due the estimate of normal vectors. Errors of orientation from normal vectors should produce large translation errors. For these experiments, the results reveal that the pairwise registration obtains an RMSE higher than 0.3 degrees and 0.2 m for the rotation error and translation error, respectively. As can be observed in Table 2, the K-4PCS achieve highest matching SR; however, the matching time is considerably higher. Examples of the results obtained with our proposed pairwise registration method can be seen in Figure 7.

Global Refinement Evaluation
As previously mentioned, for comparisons with state-of-the-art research, the renowned global refinement approach presented by [15] and the approach of Theiler et al. [14] are introduced. The results of the 3D globally consistent point clouds obtained with the proposed method can be seen in Figure 8.

Global Refinement Evaluation
As previously mentioned, for comparisons with state-of-the-art research, the renowned global refinement approach presented by [15] and the approach of Theiler et al. [14] are introduced. The results of the 3D globally consistent point clouds obtained with the proposed method can be seen in Figure 8.

Global Refinement Evaluation
As previously mentioned, for comparisons with state-of-the-art research, the renowned global refinement approach presented by [15] and the approach of Theiler et al. [14] are introduced. The results of the 3D globally consistent point clouds obtained with the proposed method can be seen in Figure 8.   Figure 9 shows the mean residual RMSE obtained for each dataset using the proposed method and the baseline approaches. The RMSE values shown that the proposed global refinement method achieves more accurate results than both the method of [15] and omitting the global refinement of the rotation parameters. It can also be seen that scans with more geometric constraints (plane correspondences with other scans) are registered more accurately after global refinement. The lowest registration accuracy is obtained for sensor poses III and VIII which have fewer corresponding planes. Figure 9 shows the mean residual RMSE obtained for each dataset using the proposed method and the baseline approaches. The RMSE values shown that the proposed global refinement method achieves more accurate results than both the method of [15] and omitting the global refinement of the rotation parameters. It can also be seen that scans with more geometric constraints (plane correspondences with other scans) are registered more accurately after global refinement. The lowest registration accuracy is obtained for sensor poses III and VIII which have fewer corresponding planes.  Figure 10 shows the RMSE value of the registration before and after the global refinement. As expected, the proposed global refinement method significantly reduces the misalignments. After the global refinement we achieved accuracies better than 8 cm, that is close to the expected measurement accuracy of the scanning device.   Figure 10 shows the RMSE value of the registration before and after the global refinement. As expected, the proposed global refinement method significantly reduces the misalignments. After the global refinement we achieved accuracies better than 8 cm, that is close to the expected measurement accuracy of the scanning device. Figure 9 shows the mean residual RMSE obtained for each dataset using the proposed method and the baseline approaches. The RMSE values shown that the proposed global refinement method achieves more accurate results than both the method of [15] and omitting the global refinement of the rotation parameters. It can also be seen that scans with more geometric constraints (plane correspondences with other scans) are registered more accurately after global refinement. The lowest registration accuracy is obtained for sensor poses III and VIII which have fewer corresponding planes.  Figure 10 shows the RMSE value of the registration before and after the global refinement. As expected, the proposed global refinement method significantly reduces the misalignments. After the global refinement we achieved accuracies better than 8 cm, that is close to the expected measurement accuracy of the scanning device.

Discussion
The proposed point cloud registration algorithm involves three prerequisite algorithms, including outlier removal, plane segmentation, and plane-based matching. The complexity of the proposed algorithm is centred on the plane-based matching algorithm. In practice, by calculating θ using the dot product of the two normal, only internal angles are verified, because such a process is limited to the interval from 0 to 180 degrees, producing ambiguity between the angles. However, a normal vector can be rotated from −180 to 180 degrees, that is, internal angles and external angles should be verified. To overcome this limitation, we use complex numbers to find the angle between the normal vectors of each pair of candidate planes referencing to the point at which this plane is rotated. Using complex numbers helps to avoid ambiguity between the angles through an analysis of the quadrants (see Equation (11)). Note also that we constraint the sign of the parameter d to be positive to orient all u towards the viewpoint. Using complex numbers to calculate the rotation angles between adjacent non-parallel vertical planes, we can decrease the number of mathematic operations comparisons with matrix operations. For the plane-based matching algorithm, plane extraction and manual setting of parameters is mandatory, as these arguments reduce the number of false positives and the number of candidate features. A close inspection of false-positive detection step revealed that the proposed plane-based mathing algorithm can reduce the number of false positives about 79%. As expected, failure cases are caused by the symmetry of the scenes, which are often in the indoor areas. As limitations of the plane-based scheme, we can point out the requirement for horizontal and vertical planes. Although surface planes are abundant in urban scenes, the behavior of the proposed strategy is not consistent for forest environments and under occlusion or very low overlap. The pairwise registration performance is better with high valid overlapping area such as for indoor environments. However, indoor environments remind us that parallel planes will increase the risk of failure for the proposed plane matching algorithm. In light of the proposed plane-based matching algorithm a minimum of three non-parallel pairs of correct plane correspondences must be found at the matching step. Fortunately, this is not a major issue in man-made environments (e.g., cities) and indoor environments that can be used in several applications, such as indoor navigation, infrastructure inspection, façade 3D modelling, cultural heritage documentation, mapping, augmented reality, 3D modelling of industrial objects, and others. With respect to the matching performance, as the rotation of the TLS is around z axis, the horizontal planes are less influenced by rotation of the TLS, while vertical planes are influenced by the rotation of the TLS. Thus, we assume that horizontal planes are those whose normal vectors and the z axis form an angle smaller than 3 degree and vertical planes are those whose normal vectors and the z axis form an angle larger than 89 degree. This is essential for reducing the search space, and for the optimization of the matching performance. As a direct consequence of the TLS scans, different plane parameters may be obtained from different views. As a result, the plane-based registration is influenced by the plane extraction accuracy. The d parameter is defined as a function of the u (normal plane) and point coordinates, while the precision of plane fitting is dependent on both the noise level of the points and the choice of the coordinate system, as presented by [28]. Since the estimation of translation parameter (t) in the plane-based solution is dependent on the d parameter of the noisy planes, a solution is normalized point coordinates to improve the estimate of d parameter in plane fitting, and consequently uses the variance of d in plane fitting to improve the estimate of translation parameters (t) during the pairwise registration solution. However, this has not been tested, and is left for future work.
Compared with other global refinement methods, our method is a closed-form solution and refines both rotational and translational parameters. From our point of view, the refinement of the rotations is essential for the correct refinement of the sensor positions. Thus, as previously mentioned, we first refine the rotations using the approach indicated by [26], while the TLS positions are corrected with our plane-to-plane refinement method. In practice, the pairs of point clouds with lowest overlapping point clouds present worst results for the pairwise registration task and the alignment is not good enough to be refined by our method. To disambiguate the pairwise transformations, geometric constraints should be added in the pose-graph optimization. Moreover, pairs of point clouds that lack salient structures, and contain more occlusions and vegetation present the largest errors. Although our global refinement method is compatible in terms of accuracy compared to the state-of-the-art, it is much faster, corresponding to 55% less computation time, as can be seen in see Figure 11.
Remote Sens. 2020, 11, x FOR PEER REVIEW 19 of 19 the largest errors. Although our global refinement method is compatible in terms of accuracy compared to the state-of-the-art, it is much faster, corresponding to 55% less computation time, as can be seen in see Figure 11. Figure 11. Results of the global registration time for each tested method.

Conclusions
In this paper, we presented a framework for automatic marker-less registration of multiple terrestrial laser scanner point clouds. First, we introduced a plane-based matching scheme that relies on a new parametrization using complex numbers to find the correspondence between pairs of segmented planes. The RANSAC algorithm is used to segment plane surfaces in TLS data. The segmented planes are then classified into vertical and horizontal planes, and their correspondences are obtained using our proposed plane-based matching algorithm. The key novel aspect of our planebased matching algorithm is the multiplication of complex numbers based on analysis of the quadrants to avoid the ambiguity in the calculation of the rotation angle formed between normal vectors of adjacent planes. It also is able to reduce the number of mathematic matrix operations during the correspondence task. Secondly, we formulated the global fine registration as a graphbased formulation adapted for plane-to-plane correspondences. The main characteristic of this proposed solution is that the refinement of TLS positions by treating the 3D points and their corresponding surface normal as observations. By globally refining the rotation parameters and the translation parameters, TLS position can be accurately obtained. Since our global registration method is non-iterative, multiple point clouds can be quickly registered. In our results, we demonstrated the potential of our method in registering point clouds of outdoor and indoor urban environments with reasonable overlapping. Future work involves conducting additional experiments using dataset of different sources, i.e., photogrammetry and RGB-D data. In addition, the loop closing detection can be executed automatically.

Conclusions
In this paper, we presented a framework for automatic marker-less registration of multiple terrestrial laser scanner point clouds. First, we introduced a plane-based matching scheme that relies on a new parametrization using complex numbers to find the correspondence between pairs of segmented planes. The RANSAC algorithm is used to segment plane surfaces in TLS data. The segmented planes are then classified into vertical and horizontal planes, and their correspondences are obtained using our proposed plane-based matching algorithm. The key novel aspect of our plane-based matching algorithm is the multiplication of complex numbers based on analysis of the quadrants to avoid the ambiguity in the calculation of the rotation angle formed between normal vectors of adjacent planes. It also is able to reduce the number of mathematic matrix operations during the correspondence task. Secondly, we formulated the global fine registration as a graph-based formulation adapted for plane-to-plane correspondences. The main characteristic of this proposed solution is that the refinement of TLS positions by treating the 3D points and their corresponding surface normal as observations. By globally refining the rotation parameters and the translation parameters, TLS position can be accurately obtained. Since our global registration method is non-iterative, multiple point clouds can be quickly registered. In our results, we demonstrated the potential of our method in registering point clouds of outdoor and indoor urban environments with reasonable overlapping. Future work involves conducting additional experiments using dataset of different sources, i.e., photogrammetry and RGB-D data. In addition, the loop closing detection can be executed automatically.