An Algorithm for Fitting Sphere Target of Terrestrial LiDAR

The sphere target played a vital role in terrestrial LiDAR applications, and solving its geometrical center based on point cloud was a widely concerned problem. In this study, we proposed a newly finite random search algorithm for sphere target fitting. Based on the point cloud data and the geometric characteristics of the sphere target, the algorithm realized the target sphere fitting from the perspective of probability and statistics with the help of parameter estimation. Firstly, an initial constraint space was constructed, and the initial center and radius were determined by finite random search. Then, the optimal spherical center and radius were determined gradually through continuous iterative optimization. We tested the algorithm with the simulated and realistic point cloud. Experimental results showed that the proposed algorithm could be effectively applied to all kinds of point cloud fitting. When the coverage rate was bigger than 30%, the fitting accuracy could reach within 0.01 mm for all kinds of point clouds. When the coverage rate was less than 20%, the fitting accuracy can reach ±1 mm, although it was reduced to a certain extent.


Introduction
Terrestrial light detection and ranging (LiDAR), also known as terrestrial laser scanning (TLS), could quickly acquire the high-resolution point cloud on the target surface by high-speed laser scanning and had brought the traditional single point measurement into the era of surface measurement. A point cloud was a set of data points in space, which was a collection of a large number of discrete measuring points on the external surface of an object. Each point position had its set of Cartesian coordinates (X, Y, Z). TLS technology had the characteristics of active, non-contact, high resolution, high precision, and rapid and flexible data acquisition. It could go deep into the complex field environment and realize the complete collection of the various large, irregular and non-standard entity or real scene 3D data. It has been wildly used in many works, such as engineering surveys, cultural relic protection, disaster monitoring, reverse 3D reconstruction et al. [1][2][3][4][5]. The sphere target (ST) had a typical spatial rotation symmetry and standard parameterized form. Part of its contour information could be obtained from any angle of view, with which the spherical center and radius could be effectively solved. Therefore, it was widely used in the multi-class application research of terrestrial LiDAR, such as the calibration and check of a terrestrial laser scanner, scanning accuracy evaluation, registration, and georeferencing of point clouds et al. [6][7][8][9][10][11]. The geometric center of the sphere target was inside the sphere and could not be obtained directly through measurement. Usually, what we utilized TLS to collect directly was the point cloud data on the surface of the target ball. Based on such data to determine its internal geometric center, this involved the fitting problem of the sphere target. At the same time, the sphere fitting problem was also a common problem to be solved in object tracking, pattern recognition, robotics, camera calibration, and other research work [12][13][14][15][16][17].
Judging from the related literature, some scholars have conducted related research on the problem of sphere fitting. Forbes took the center and the radius of the sphere as the parameters to be sought and analyzed the fitting algorithms of several types of spheres and other geometric bodies. These algorithms were mainly suitable for noise-free point cloud data with a high coverage rate (CR) [18]. Nievergelt used a least-squares method based on algebraic distances to calculate the center of the sphere. Although his method had advantages in computational efficiency, it usually did not provide satisfactory results [19]. Späth, Shakarji, and Ahn, et al. used improved least-squares methods to perform sphere fitting [20][21][22][23]. Clouse used conjugate gradient descent to calculate the sphere's center, which used both cost function evaluations, and evaluations of the derivative to find a set of parameters that produce a local minimum cost [24]. Witzgall respectively used algebraic fitting and geometric fitting to perform sphere fitting. With the help of the concept of deviation between data point and sphere, the arithmetic fitting was solved by leastsquare through linear regression. The geometric fitting used the orthogonal least-squares solution [25]. Sumith used a fast geometric method to fit the center and radius of the sphere, and the fitting accuracy was better than the ordinary least squares estimator (OLS) [26]. Liu used a nonlinear least-squares method to achieve sphere fitting [27]. Fei used a constrained nonlinear least-squares fitting (CNLSF) algorithm to realize the fitting of spheres with a small segment angles strategy [28]. Lesouple used an expectation-maximization method to achieve the fitting of spheres [29].
At present, most of the sphere fitting algorithms mainly rely on least-squares minimization methods to obtain their centers, such as linear least-squares, nonlinear least-squares, the total least squares method as well as the weighted total least squares method to eliminate the error of the coefficient matrix [30][31][32]. From the theory of least squares, the least-squares estimation assumed that the mean of data noise was zero, resulting in an unbiased parameter estimation. If the noise variance was known, the minimum variance parameter estimation could be obtained by selecting appropriate weights on the data. In addition, least squares estimation implicitly assumed that the entire data set could only be explained by one parameter vector of a given model [33,34]. Numerous studies have clearly shown that least-squares estimation could easily violate these assumptions. Sometimes, even if the data contained only one "bad" datum, the least-squares estimate may be seriously disturbed, resulting in low calculation accuracy. In addition to the least-squares method, there were also some other methods, such as a minimum zone sphere, maximum inscribed sphere, minimum circumscribed sphere [35][36][37]. These methods mainly take advantage of linearization to fit the sphere with the help of mathematics or geometry. The sphere target fitting itself was a nonlinear problem, which inevitably led to the loss of accuracy in the linearization process. At the same time, the number of points in a sphere target point cloud was usually more than thousands, which would cause a large calculation matrix and low computational efficiency.
As we all know, in TLS work, no matter what type of sphere target we used, it had a specific geometric size, that is to say, the spatial distribution of the point cloud of any sphere target had a particular range that we call bounding box, and this bounding box contains the whole point cloud of the sphere target, including the noises. From the geometric characteristics, the geometric center and radius of the sphere target must be in this bounding box. Thus, we could adopt a search strategy to look for the optical center and radius in the bounding box that satisfies the specific error criteria. In this study, combining the point cloud and geometric characteristics of point cloud, we developed a finite random search alogorithm for the sphere target fitting. Our proposed algorithm mainly aimed to achieve a better sphere target fitting after the point cloud extraction of a singular sphere target has been completed. Its main objective is to calculate the geometric center accurately based on the point cloud data of a single target sphere. The detailed design of the algorithm is described in Section 2. In this paper, we did not discuss how to extract point cloud data of an individual target sphere from a complex point cloud, but there were many solutions for this problem [38].

Methods and Data
Given a point cloud of a sphere target T = {(x i , y i , z i )|i = 1, 2, · · · , n ≥ 3} obtained by TLS, let (X, Y, Z) be the unknown center, and let R be the unknown radius of the sphere target. In a specific scanning coordinate system, both the geometric center (X, Y, Z) and radius R of the sphere target were determined. During the data acquisition process, affected by factors such as the instrument itself and the external environment, a point cloud was inevitably mixed with noise [39,40]. Sphere target fitting was to extract the center and radius of the sphere target from the point cloud with unknown distribution and outliers. This could be viewed as an optimal parameter estimation problem. In this problem, we regarded the geometric center (X, Y, Z) and radius R of the sphere target as the parameters to be solved and took the target point cloud as the observation value. Using the point cloud to fit the geometric center and radius could be regarded as finding the optimal parameters that meet the specific decision rules.
We took the centroid of the sphere target point cloud as the center and took more than two times the radius length as the constraint to construct an initial bounding box. According to the geometric characteristics of the sphere target, its geometric center and radius must be within the bounding box. Based on this feature, we could solve the problem of the sphere target fitting by using the idea of probability theory and parameter estimation.
n} be composed of four characteristic quantities, where (X i , Y i , Z i ) was the potential geometric center of the target sphere and R i was the potential geometrical radius of the target sphere. The four characteristic quantities (X, Y, Z, R) were continuous variables, and their values should be infinite in theory. From the perspective of probability and statistics, in the process of finite random search, the probability of obtaining the optimal value was related to the size of the sample space. The larger the sample space, the lower the probability of finding the optimal value. Conversely, the smaller the sample space, the higher the probability of finding the optimal value [41,42]. In this study, we proposed a finite random search algorithm suitable for sphere target fitting combined with the point cloud and geometric characteristics of the sphere target. The primary technical process of the algorithm is shown in Figure 1.

Initial Parameters
The initial parameters were mainly composed of the random search times N loop , the iterative optimization times N opt , the estimated radius of the target sphere R set , the total error threshold E min and the scale scaling factor α.

•
N loop referred to the number of random searches for the optimal value in a set-limited search space. • N opt referred to the predetermined number of times to update the search space in the iterative search process. • R set referred to the pre-estimated geometric radius of the sphere target. In TLS work, the geometric radius of the sphere target used was usually known, and if the radius could not be determined, R set could be set to a relatively sizeable rough value. • E min referred to the preset value of the total error that determined the end of the fitting prematurely. Setting the threshold here could realize that after the predetermined fitting accuracy was reached, the fitting process was ended early to improve the execution efficiency of the algorithm. • α referred to the adjustment factor of the constrained space scaling when the constrained space was updated.

Error Criteria
The fitting problem of the sphere target point cloud could be regarded as the problem of estimating parameters from noisy data, where the judgment criterion played an important role. That would influence the accuracy of the estimated parameters, computation efficiency, and the robustness to predictable or unpredictable errors. To find the optimal center and radius of the target sphere, we must first determine the optimal center and radius measurement criteria. The error criteria of the existing sphere target fitting algorithm were mainly divided into geometric error and arithmetic error, each of which has its advantages and disadvantages [43]. In our proposed algorithm, an arithmetic error was chosen as the criterion.
Let the current fitting center of the target sphere be (X, Y, Z) and the radius R. The current total error E total could be calculated by Equation (1).
where (x i , y i , z i ) were the coordinates of the measured points of the point cloud, and n was the number of the measuring points contained in the point cloud.

Initial Parameters
The initial parameters were mainly composed of the random search times loop , the iterative optimization times opt , the estimated radius of the target sphere set , the total error threshold and the scale scaling factor . In the proposed algorithm, a limited sample of possible solutions was generated with the help of the constraint space, and their total error was calculated by using Equation (1) one by one. The sample with the smallest total error was taken as the optimal solution under the current constraint space. More detailed usage is described in Sections 2.3 and 2.4. Theoretically, the smaller the total error E total was, the higher the coincidence degree between the point cloud and simulated sphere surface was, the more accurate the current center and radius were; otherwise, the deviation of the center or radius was more significant.

Initial Constraint Space and Sample Construction
In terms of the geometric characteristics of the sphere target, given the point cloud of a sphere target, its centroid could be solved, and this centroid could be used as the geometric center of the bounding box of the point cloud. By adjusting the side length of the bounding box, we could construct a space cube containing the point cloud of the sphere target and the complete range of the sphere target to be fitted. Then, the center of the sphere target must be in this space cube. We called this cube the constrained space S, and the cube constructed for the first time was the initial constraint space S (1) . In the algorithm in this paper, the initial constraint space was mainly constructed based on the point cloud and the preset parameter R set , and the construction principle is shown in Figure 2. (1) The centroid ( ̅ , , ̅ ) of the sphere target point cloud could be calculated and used as the center of the initial constraint space.
(2) With O as the geometric center and as the constraint, the space ranges of the X, Y, and Z axes of the initial constraint space and radius range were determined by Equation (2). In general, was set to be slightly larger than the real radius of the sphere target to ensure that the constructed initial constraint space covered the entire target sphere. In practical scanning operations, the sphere target used usually had a clear geometric size, so the estimated radius of the sphere target should be easy to determine.
(3) With the initial constraint space and radius as constraints, a sample space ( ) = ( ) , ( ) , ( ) , ( ) | = 1,2, ⋯ , composed of four characteristic quantities was constructed, and it contained randomly generated samples.  (1) The centroid O(x, y, z) of the sphere target point cloud could be calculated and used as the center of the initial constraint space.
(2) With O as the geometric center and R set as the constraint, the space ranges of the X, Y, and Z axes of the initial constraint space and radius range were determined by Equation (2). In general, R set was set to be slightly larger than the real radius of the sphere target to ensure that the constructed initial constraint space covered the entire target sphere. In practical scanning operations, the sphere target used usually had a clear geometric size, so the estimated radius of the sphere target should be easy to determine. (3) With the initial constraint space and radius as constraints, a sample space (1) i | i = 1, 2, · · · , N loop composed of four characteristic quantities was constructed, and it contained N loop randomly generated samples.
(4)After determining the value ranges of the four characteristic quantities in the sample space, the current scale of the characteristic quantities L (1) The E total of each sample in U (1) was calculated in sequence by the error criterion of Section 2.2, and compared with the preset total error threshold E min . If E total was less than E min , it indicated that the current sample had reached the predetermined accuracy of parameter estimation. At this point, this sample could be considered as the optimal solution of parameters. Otherwise, the search should continue until all samples in U (1) were tested. After traversing all the samples in U (1) , if no sample satisfying E total < E min was found, the sample with the smallest E total was selected as the optimal solution of the parameter in the current constraint space.

Constraint Space and Sample Update
From the perspective of probability and statistics, in the process of finite random search, the probability of obtaining the optimal value was related to the size of the constraint space. The smaller the constraint space, the higher the probability of finding the optimal value [44,45]. In the last constraint space, the center and radius determined by finite random search did not meet the set accuracy, so the sample space needed to be further optimized to accurately determine the center and radius.
Suppose that in the constraint space S (i) , the optimal solution of center and radius obtained from N loop samples was X (i) , Y (i) , Z (i) , R (i) , and the scale of the characteristic quantity was L (i) (1) Each feature quantity could be scaled by using Equation (4). Here α ∈ (0, 1) was the preset scale scaling factor, which was used to adjust the scaling speed of each feature quantity. The smaller the value was, the faster the scaling speed was.
as the center and the updated feature quantity scale as the constraint, a new constraint space S (i+1) was constructed. The value ranges of the four characteristic quantities of each sample in sample Equation (5). In this case, the value of R (i+1) min should be paid attention. The radius of the sphere target should not be negative, so R (i+1) min must always be greater than 0.
(3) With the updated constraint space S (i+1) and radius as constraints, a new sample space U (i+1) containing N loop samples was randomly generated.
According to the same search strategy, the E total of each sample in the sample space U (i+1) was calculated one by one and detected whether it was less than the preset threshold E min of the total error. After traversing all the samples in U (i+1) , if no sample satisfying E total < E min was found, the sample with the minor total error was selected as the optimal solution of the parameters in the current constraint space.

Iterative Optimization
Based on the constraint space constructed, the finite random sample space was generated, and the optimal sample in the sample space was determined as the optimal solution in the current constraint space. The optimal solution satisfying the predetermined accuracy was gradually found by constantly updating the search and sample spaces. In general, after more than ten iterations of optimization, the optimal value could be determined. In the iterative optimization process, the search process may end when E total was less than E min , or it may end after the predetermined N opt iterative optimization. At the same time, considering the noise contained in the point cloud, the E total obtained may always fail to reach E min during the next finite iteration of N opt , but after a finite iteration search, E total might always stop decreasing. At this time, it was necessary to detect a situation that if the E total did not decrease during two consecutive iterations, the iterative process should end.

Accuracy Estimation
The fitting accuracy of the sphere target center directly determined the scope of application of the algorithm. When the center of the sphere target was known, the fitting accuracy could be evaluated by the root mean square error (RMSE) of the center [46,47].
Let the real center of the sphere target be (X o , Y o , Z o ) and the fitted center X f , Y f , Z f . The RMSE of the fitted center could be obtained by Equation (6).
where in this study, N = 1. Because the center of the sphere target was located inside the sphere, it could not be measured directly, so it was difficult to accurately evaluate the fitting accuracy of the center of the sphere target in practical work. To effectively test the fitting effect of the proposed algorithm, simulation data with known center and radius were used to test the algorithm.

Experiment
To test the feasibility of the algorithm proposed in this paper, we designed two groups of simulated data without noise point cloud and noisy point cloud for testing. The noisefree point cloud mainly simulates the data processed by various denoising operations. The noisy point cloud mainly simulates the data obtained by TLS directly and without previous noise processing. The test platform uses Windows 10 system, Intel I5-6500 CPU, 16G memory, and Intel Graphics 530 Graphics card.

Point Cloud Simulation
A target sphere was a standard spherical geometry whose point cloud consisted of several measuring points on the surface of the target sphere. The spherical coordinates (x, y, z) of any measurement point p in the sphere target point cloud could be determined by the spherical radius r, the zenith angle θ ∈ 0, 180 • and the plane projection angle ϕ ∈ 0, 360 • , as shown in Figure 3. In the scanning coordinate system, the scanning coordinate (X, Y, Z) of measuring point p was affected not only by spherical coordinates but also by the position of spherical coordinates (x 0 , y 0 , z 0 ) in the scanning coordinate system. At this point, the scanning coordinate of measuring point p should be calculated by Equation (7). The sphere target point cloud was simulated by equally dividing the zenith angle θ and plane projection angle ϕ.
where in this study, N = 1. Because the center of the sphere target was located inside the sphere, it could not be measured directly, so it was difficult to accurately evaluate the fitting accuracy of the center of the sphere target in practical work. To effectively test the fitting effect of the proposed algorithm, simulation data with known center and radius were used to test the algorithm.

Experiment
To test the feasibility of the algorithm proposed in this paper, we designed two groups of simulated data without noise point cloud and noisy point cloud for testing. The noise-free point cloud mainly simulates the data processed by various denoising operations. The noisy point cloud mainly simulates the data obtained by TLS directly and without previous noise processing. The test platform uses Windows 10 system, Intel I5-6500 CPU, 16G memory, and Intel Graphics 530 Graphics card.

Point Cloud Simulation
A target sphere was a standard spherical geometry whose point cloud consisted of several measuring points on the surface of the target sphere. The spherical coordinates ( , , ) of any measurement point p in the sphere target point cloud could be determined by the spherical radius r, the zenith angle ∈ [0, 180°] and the plane projection angle ∈ [0, 360°], as shown in Figure 3. In the scanning coordinate system, the scanning coordinate (X, Y, Z) of measuring point p was affected not only by spherical coordinates but also by the position of spherical coordinates ( , , ) in the scanning coordinate system. At this point, the scanning coordinate of measuring point p should be calculated by Equation (7). The sphere target point cloud was simulated by equally dividing the zenith angle θ and plane projection angle φ.  Due to the influence of TLS scanning field of view, occlusion of the target sphere itself, scanning distance between TLS and sphere target, single station scanning could only obtain target sphere point cloud data with maximum coverage of 50% [48,49]. Coverage rate (CR) here was the percentage of the surface area covered by the pointing cloud to the total area of the target sphere. It was mainly calculated based on the proportion of the surface area S occupied by the point cloud to the sphere's surface area S, as shown in Figure 4.
Due to the influence of TLS scanning field of view, occlusion of the target sphere itself, scanning distance between TLS and sphere target, single station scanning could only obtain target sphere point cloud data with maximum coverage of 50% [48,49]. Coverage rate (CR) here was the percentage of the surface area covered by the pointing cloud to the total area of the target sphere. It was mainly calculated based on the proportion of the surface area S′ occupied by the point cloud to the sphere's surface area S, as shown in Figure 4. The CR simulation of experimental data was realized by adjusting the range of zenith angle θ. According to the surface area calculation formula of the spherical crown and sphere, the coverage rate c of the simulated point cloud could be calculated by Equation  The CR simulation of experimental data was realized by adjusting the range of zenith angle θ. According to the surface area calculation formula of the spherical crown and sphere, the coverage rate c of the simulated point cloud could be calculated by Equation (8), and the meanings of each parameter are shown in Figure 4 [50,51].
where S was the surface area of the spherical crown, S was the surface area of the sphere, h was the height of the spherical crown, and r was the radius of the sphere.
In the simulation of a sphere target point cloud, the coverage rate c of the point cloud was given, and the value range of the zenith angle θ ∈ (0, θ max ) should be determined according to the CR c, where θ max could be calculated by Equation (9).
where OO was the distance from the sphere's center to the underside of the spherical crown, and other symbols were the same as Equation (8).
Influenced by the instrument performance and external environment (wind, air humidity, illumination, etc.), noise would inevitably be mixed in the point cloud obtained by TLS. According to the error theory of geodesy, the noise in the measured data usually satisfied the Gaussian distribution N µ, σ 2 , where µ and σ 2 were the expected and variance of the Gaussian distribution, respectively. In the simulation of experimental data, Gaussian noise was added to the point cloud to simulate the real noise in the scanning process. In all noise point cloud simulations, µ was set to 0, and σ was set to the maximum deviation value of the point cloud.

Noise-Free Point Cloud
For noise-free point clouds of sphere target with different coverage rates, five simulated data were generated by using method 3.1, as shown in Figure 5. Among them, coverage rates of (a)~(e) were 50%, 40%, 30%, 20% and 10%, respectively. The center and radius (X, Y, Z, R) of all the simulated point clouds were (1000, 1000, 100, 0.0725), and the unit was the meter. The sampling interval of both the zenith angle θ and the plane projection angle ϕ was 3 • . We use the proposed algorithm to fit the five simulated point clouds, and the fitting process is shown in Figure 5. For the initial parameter setting in the fitting process, N loop , N opt , R set , E min , and α were 1000, 30, 0.08 m, 0.001, and 0.25, respectively.

Noise-Free Point Cloud
For noise-free point clouds of sphere target with different coverage rates, five simulated data were generated by using method 3.1, as shown in Figure 5. Among them, coverage rates of (a)~(e) were 50%, 40%, 30%, 20% and 10%, respectively. The center and radius (X, Y, Z, R) of all the simulated point clouds were (1000, 1000, 100, 0.0725), and the unit was the meter. The sampling interval of both the zenith angle θ and the plane projection angle φ was 3°. We use the proposed algorithm to fit the five simulated point clouds, and the fitting process is shown in Figure 5. For the initial parameter setting in the fitting process, loop , opt , set , , and α were 1000, 30, 0.08 m, 0.001, and 0.25, respectively. fitting process of (a), (g) fitting process of (b), (h) fitting process of (c), (i) fitting process of (d), (j) fitting process of (e).
After all sphere target fitting was completed, we calculated statistics on the number of points in the point cloud, the RMSE of the fitting center, the total error at the end of the fitting, iteration times, and running time, as shown in Table 1. According to the statistical results, for the noise-free point cloud, when the coverage rate reaches more than 30%, the total error threshold of 0.001 could be achieved, and the fitting center's RMSE was less fitting process of (a), (g) fitting process of (b), (h) fitting process of (c), (i) fitting process of (d), (j) fitting process of (e).
After all sphere target fitting was completed, we calculated statistics on the number of points in the point cloud, the RMSE of the fitting center, the total error at the end of the fitting, iteration times, and running time, as shown in Table 1. According to the statistical results, for the noise-free point cloud, when the coverage rate reaches more than 30%, the total error threshold of 0.001 could be achieved, and the fitting center's RMSE was less than 0.01 mm after less than 15 iterations of optimization. When the coverage was below 20%, the iterative optimization times would increase, but generally, the fitting would be completed after about 25 iterations of optimization. At this time, the total error did not reach the predetermined total error threshold of 0.001, and fitting should be the end of the iteration process when the adjacent total error did not decrease, and both the deviations of X, Y, and Z axes and the fitting center's RMSE were all less than 1 mm. Experimental results showed that the iterative optimization times were mainly affected by the coverage rate in the absence of noise. When the coverage rate was greater than 30%, ideal fitting results could be obtained after 15 iterations of optimization. When the coverage was less than 20%, the number of iterations would increase, but generally, the fitting could be completed after less than 25 iterations. The running time was mainly affected by the number of measuring points in the point cloud. When the coverage rate was high, the number of measuring points was large, and the running time was extended. When the coverage rate was low, the measuring point data was small, and the running time was short. From the calculation time of 10 simulated data, the fitting could be completed in less than 0.5 s.

Noisy Point Cloud
Affected by instrument performance, scanning distance, incident angle, reflection intensity, surface roughness, the external environment (e.g., ambient vibration, wind, temperature), and other factors, the TLS point cloud cannot avoid mixed noise. From various scanners' existing nominal technical parameters, the ranging error within the scanning distance of 100 m was basically within ±2.0 mm. Considering the influence of the surrounding environment and other factors, we add ±5.0 mm noise to the simulation data. We used method 3.1 to simulate five sphere target point clouds with different coverage rates of ±5 mm noise, as shown in Figure 6. Among them, the coverage rates of (a)~(e) were 50%, 40%, 30%, 20% and 10%, respectively. The center and radius (X, Y, Z, R) of all the simulated point clouds were (1000, 1000, 100, 0.0725), and the unit was the meter. The sampling interval of both the zenith angle θ and the plane projection angle ϕ was 3 • . We use the proposed algorithm to fit the five simulated point clouds, respectively, and the fitting process is shown in Figure 6. The initial parameter setting of the algorithm was the same as that of the noise-free point cloud fitting, N loop , N opt , R set , E min , and α were 1000, 30, 0.08 m, 0.001, and 0.25, respectively.
were 50%, 40%, 30%, 20% and 10%, respectively. The center and radius (X, Y, Z, R) of all the simulated point clouds were (1000, 1000, 100, 0.0725), and the unit was the meter. The sampling interval of both the zenith angle θ and the plane projection angle φ was 3°. We use the proposed algorithm to fit the five simulated point clouds, respectively, and the fitting process is shown in Figure 6. The initial parameter setting of the algorithm was the same as that of the noise-free point cloud fitting, loop , opt , set , , and α were 1000, 30, 0.08 m, 0.001, and 0.25, respectively. process of (a), (g) fitting process of (b), (h) fitting process of (c), (i) fitting process of (d), (j) fitting process of (e).
After every sphere target fitting was completed, we made statistics on the number of points in the point cloud, the RMSE of the fitting center, the total error at the end of the fitting, iteration times, and running time, as shown in Table 2. According to the statistical results, for the noisy point clouds with different coverage rates, the fitting process ends when the total errors stop decreasing during the iterative optimization process, and the total errors at the end of the process do not reach the preset total error threshold of 0.001. When the coverage was above 30%, the fitting center with RMSE less than 0.001 mm could be obtained after about 25 iterations of optimization. When the coverage rate was 20%, the fitting center with RMSE less than 1 mm could be obtained after about 20 iterations.  CR 10, (f) fitting process of (a), (g) fitting process of (b), (h) fitting process of (c), (i) fitting process of (d), (j) fitting process of (e).
After every sphere target fitting was completed, we made statistics on the number of points in the point cloud, the RMSE of the fitting center, the total error at the end of the fitting, iteration times, and running time, as shown in Table 2. According to the statistical results, for the noisy point clouds with different coverage rates, the fitting process ends when the total errors stop decreasing during the iterative optimization process, and the total errors at the end of the process do not reach the preset total error threshold of 0.001. When the coverage was above 30%, the fitting center with RMSE less than 0.001 mm could be obtained after about 25 iterations of optimization. When the coverage rate was 20%, the fitting center with RMSE less than 1 mm could be obtained after about 20 iterations. When the coverage rate was 10%, the fitting center's RMSE was about 2 mm after about 20 iterations of optimization. At this time, the error of X, Y, and Z axes primarily occurs in the Z-axis, because the point cloud and noise were mainly concentrated in the Z-axis, which may affect the fitting accuracy of the sphere target center in this direction to some extent. It could be seen from the comparison of the fitting results of noise-free point cloud and noisy point cloud when the point cloud was mixed with noise, the number of iterative optimization of point cloud with coverage of more than 30% increases obviously, while the number of iterative optimization of point cloud with coverage of less than 20% does not change significantly. The accuracy of the fitting center was mainly affected by the coverage rate. Whether there was noise mixed in the point cloud or not, the RMSE of the fitting center was less than 0.01 mm when the coverage rate reached more than 30%. The fitting accuracy of the point cloud with a 20% coverage rate was about 1 mm. The point cloud with a 10% coverage rate was susceptible to noise, and the fitting accuracy would be reduced to a certain extent after mixing with noise. The running time was mainly restricted by the number of measuring points in the point cloud, but noise and coverage rate was not significantly affected. The more the number of measuring points, the longer the runtime would be.

Realistic Point Cloud
In order to test the applicability of the proposed algorithm to real target balls, the point clouds of five real sphere targets were acquired by TLS. The experimental site was selected in a small square on the campus of Nanjing Tech University, where five sphere targets with different distances were arranged, and the Faro Laser Scanner Focus 3D X330 was used to collect their point clouds. The arrangement of the Target balls is shown in Figure 7a. According to their distance from the scanner, they were named Target 1~5 respectively from near to far. The sphere target selected was Faro's standard sphere target, whose real geometric radius was 0.0725 m. The FARO Focus 3D X330 was a high-speed 3D scanner with extra-long range, which could scan objects up to 330 m away even in direct sunlight. Its nominal ranging error was ±2 mm, and ranging noise was 0.3 mm at a distance of 10 and 25 m when the reflectivity was 90%. The point cloud of the whole scene was obtained through single-site scanning according to the point spacing setting of 7.67 mm@10 m, as shown in Figure 7b. According to the measurement principle of a 3D laser scanner, the spatial distribution of scanned point cloud showed divergence, that is, the farther away from the scanner, the larger the point spacing and the smaller the spatial density, which could be clearly shown in the orthographic projection of the scanned point cloud in Figure 7b. From the original point clouds obtained without any denoising process, point clouds of five sphere targets were manually extracted, as shown in Figure 7c. In practical scanning work, the real geometric center of the sphere target could not be measured directly, so its true value cannot be accurately obtained. In order to effectively compare the fitted results of the proposed algorithm with other software or algorithms, the geometric centers of five target spheres were extracted by the point cloud processing software of Faro's SCENE and regarded as reference values. In SCENE software, the fitting of the sphere target needed to set the radius of the sphere target in advance, which was set as the real radius of the target ball (0.0725 m). Meanwhile, our algorithm and the least square (LS) algorithm were used to fit the five point clouds, respectively, and the fitting results were compared with the reference values, as shown in Table 3. Here, dX, dY, dZ and dR were the absolute values of deviations between the fitted In practical scanning work, the real geometric center of the sphere target could not be measured directly, so its true value cannot be accurately obtained. In order to effectively compare the fitted results of the proposed algorithm with other software or algorithms, the geometric centers of five target spheres were extracted by the point cloud processing software of Faro's SCENE and regarded as reference values. In SCENE software, the fitting of the sphere target needed to set the radius of the sphere target in advance, which was set as the real radius of the target ball (0.0725 m). Meanwhile, our algorithm and the least square (LS) algorithm were used to fit the five point clouds, respectively, and the fitting results were compared with the reference values, as shown in Table 3. Here, dX, dY, dZ and dR were the absolute values of deviations between the fitted center and radius of the sphere target and the reference value, respectively, and RMSE was the error value of the fitted center, which could be calculated by Equation (5). The preset parameters required in our algorithm were the same as those in Section 3.2. When LS was applied to fit, the initial values of center and radius needed to be given, which were respectively set as the centroid of the point cloud and the real radius of the sphere target. From the scanning and fitting experiments of real sphere targets, the following conclusions could be drawn.
(1) The points and coverage rate of the point cloud were directly affected by the distance between the sphere target and the scanner. It could be seen from Figure 8a that as the distance between the sphere target and the scanner increased, both the number of measuring points and the coverage rate decreased accordingly, which was determined by the performance of the instrument. In actual scanning work, the coverage rate was usually less than 40%. For example, Target 1, which was 3.316 m away from the scanner, had a coverage rate of only 35%.
(2) Our algorithm was efficient in real sphere target fitting. From the iterative optimization times and runtime, our algorithm could complete the fitting after less than 20 iterative optimizations, and the runtime was less than 0.5 s, as shown in Figure 8b.
(3) The fitting accuracy of our algorithm was comparable to that of commercial software SCENE. Under the assumption that the centers of the sphere targets by SCENE were the true value, the deviation of X, Y and Z and RMSE of the fitted center of our algorithm were all less than 1 mm, as shown in Table 3. From another perspective, the applicability of our algorithm was better than that of commercial software SCENE. The reason was that in SCENE's fitting work, the true radius of the sphere target must first be accurately set, but our algorithm only needed a rough estimate, and it would then be automatically optimized. From the experiments we conducted, setting the radius in our algorithm to a known value would improve the efficiency and fitting accuracy of the algorithm to a certain extent. However, considering the versatility of the algorithm, it was still chosen as an unknown parameter to be solved here.
(4) The fitting accuracy and noise immunity of our algorithm were better than that of the least squares algorithm. It can be seen from Figure 7c that Target 1 and Target 2 had no obvious noise. At this point, the fitting accuracy of the two methods was equivalent. Target 3~4 all contained obvious noises. Especially in the case of obvious outliers in Target 3, our algorithm could still achieve a fitting accuracy of RMSE less than 1 mm, while LS had an obvious large deviation, as shown in Figure 8c. The radius of the actual sphere target used in the experiment was known. From the fitting error of radius, the fitting error of the two algorithms was less than 1 mm when there was no obvious noise influence. However, when there was obvious noise, our algorithm could still be applied stably, while LS was greatly disturbed and had serious deviation, as shown in Figure 8d.
Target 3~4 all contained obvious noises. Especially in the case of obvious outliers in Target 3, our algorithm could still achieve a fitting accuracy of RMSE less than 1 mm, while LS had an obvious large deviation, as shown in Figure 8c. The radius of the actual sphere target used in the experiment was known. From the fitting error of radius, the fitting error of the two algorithms was less than 1 mm when there was no obvious noise influence. However, when there was obvious noise, our algorithm could still be applied stably, while LS was greatly disturbed and had serious deviation, as shown in Figure 8d.

Parameters and Efficiency
The preset radius set of the sphere target directly affects the size of the initial space constructed. In TLS work, the sphere target used usually has a fixed geometric radius, in which case set could be easily determined. If the radius of the sphere target cannot be determined, a fairly large value could be set. After several iterations of optimization, the range of the radius could be basically determined. According to the experimental results, set should be controlled within six times of the actual radius of the sphere target, so that the comparatively accurate radius could be obtained after five iterations. If the value was too large, the fitting accuracy would be low.
When updating the constraint space, the scaling factor α directly determines the contraction speed of the constraint space, and the smaller the value is, the faster the contraction speed is. The experimental results show that α was generally set to about 0.25 to ensure a good fitting effect. At the same time, the setting of α should take into account the point cloud coverage. When the CR was more excellent than 20%, α could be reduced to speed up the search. When the point cloud coverage was less than 20%, the α could be increased to improve the accuracy of the solution.
The random search times loop directly determine the number of samples generated in the finite constraint space, and its value was recommended to be between 500 and 1000. If loop was too large, the number of sample points was too large, resulting in low fitting efficiency; if loop was too small, the number of samples would be insufficient, resulting in insufficient fitting accuracy. According to the experimental results, when loop was set to 1000, all kinds of point clouds could be fitted within 1 s, which was acceptable to us, as shown in Figure 9a. The value of opt should be adjusted to that of loop . If the value of loop was relatively small, such as 500, the value of opt should be increased. If the loop was correspondingly large, such as 1000, it should be possible to reduce the opt appropriately. According to the experimental results, when loop was set to 1000, an ideal fitting result could generally be achieved within 30 times, as shown in Figure 9b.

Parameters and Efficiency
The preset radius R set of the sphere target directly affects the size of the initial space constructed. In TLS work, the sphere target used usually has a fixed geometric radius, in which case R set could be easily determined. If the radius of the sphere target cannot be determined, a fairly large value could be set. After several iterations of optimization, the range of the radius could be basically determined. According to the experimental results, R set should be controlled within six times of the actual radius of the sphere target, so that the comparatively accurate radius could be obtained after five iterations. If the value was too large, the fitting accuracy would be low.
When updating the constraint space, the scaling factor α directly determines the contraction speed of the constraint space, and the smaller the value is, the faster the contraction speed is. The experimental results show that α was generally set to about 0.25 to ensure a good fitting effect. At the same time, the setting of α should take into account the point cloud coverage. When the CR was more excellent than 20%, α could be reduced to speed up the search. When the point cloud coverage was less than 20%, the α could be increased to improve the accuracy of the solution.
The random search times N loop directly determine the number of samples generated in the finite constraint space, and its value was recommended to be between 500 and 1000. If N loop was too large, the number of sample points was too large, resulting in low fitting efficiency; if N loop was too small, the number of samples would be insufficient, resulting in insufficient fitting accuracy. According to the experimental results, when N loop was set to 1000, all kinds of point clouds could be fitted within 1 s, which was acceptable to us, as shown in Figure 9a. The value of N opt should be adjusted to that of N loop . If the value of N loop was relatively small, such as 500, the value of N opt should be increased. If the N loop was correspondingly large, such as 1000, it should be possible to reduce the N opt appropriately. According to the experimental results, when N loop was set to 1000, an ideal fitting result could generally be achieved within 30 times, as shown in Figure 9b. The total error threshold was the preset condition for the end of the iteration. In the iterative optimization process, if the of the current sample was detected to be less than , the iteration would end in advance, and this sample would be taken as the optimal solution of the sphere target fitting. The value should not be too large, generally less than 0.001. In the algorithm, the check judgment of iterative optimization accuracy was set. If the fitting accuracy of two adjacent iterations was no longer improved, the judgment would be terminated in advance. This judgment was very effective in the fitting process of the noisy point cloud. As shown in Table 2, for the noisy point cloud, the at the end of fitting did not reach , but due to the influence of noise, the precision of in the adjacent iterative process did not improve, so the fitting was ended before reaching the preset iteration optimization number opt .
The time complexity of the proposed algorithm was O loop * opt [52], and the efficiency of the algorithm was directly restricted by loop and opt . It could be seen from the experimental results that several parameters were not absolutely and completely independent, but have some influence on each other. Therefore, the parameter reference value given here was an empirical value, which should be adjusted appropriately according to the actual situation in the practical application of the algorithm.

Robustness and Accuracy
Robustness was an important feature of an algorithm. It directly determines its application scope. The proposed algorithm not only makes use of the point cloud data itself but also considers the geometric characteristics of the target sphere, and realizes the target sphere fitting, employing probability and statistic theory. As shown in Figure 10, from the fitting process of 10 simulated sphere targets, regardless of whether there was noise in the sphere target point cloud or not, with the optimization of the constraint space, the total error would become smaller and smaller, that is, the fitting accuracy would gradually improve. Generally, it would be basically stable after five iterations, and further optimization would follow.   The total error threshold E min was the preset condition for the end of the iteration. In the iterative optimization process, if the E total of the current sample was detected to be less than E min , the iteration would end in advance, and this sample would be taken as the optimal solution of the sphere target fitting. The value should not be too large, generally less than 0.001. In the algorithm, the check judgment of iterative optimization accuracy was set. If the fitting accuracy of two adjacent iterations was no longer improved, the judgment would be terminated in advance. This judgment was very effective in the fitting process of the noisy point cloud. As shown in Table 2, for the noisy point cloud, the E total at the end of fitting did not reach E min , but due to the influence of noise, the precision of E total in the adjacent iterative process did not improve, so the fitting was ended before reaching the preset iteration optimization number N opt .
The time complexity of the proposed algorithm was O N loop * N opt [52], and the efficiency of the algorithm was directly restricted by N loop and N opt . It could be seen from the experimental results that several parameters were not absolutely and completely independent, but have some influence on each other. Therefore, the parameter reference value given here was an empirical value, which should be adjusted appropriately according to the actual situation in the practical application of the algorithm.

Robustness and Accuracy
Robustness was an important feature of an algorithm. It directly determines its application scope. The proposed algorithm not only makes use of the point cloud data itself but also considers the geometric characteristics of the target sphere, and realizes the target sphere fitting, employing probability and statistic theory. As shown in Figure 10, from the fitting process of 10 simulated sphere targets, regardless of whether there was noise in the sphere target point cloud or not, with the optimization of the constraint space, the total error E total would become smaller and smaller, that is, the fitting accuracy would gradually improve. Generally, it would be basically stable after five iterations, and further optimization would follow. Figure 10. Detailed fitting process of experimental data: (a-e) were the total error changes during the fitting process of five simulated noise-free point clouds with different coverage rates, (f-j) were the total error changes during the fitting process of five simulated noisy point clouds with different coverage rates.
The sphere target fitting method based on the least-squares completely depends on the target sphere point cloud data, and the noise in the data was easy to cause the illconditioned coefficient matrix, which leads to the decrease in fitting accuracy, especially in the case of low coverage rate, it was easy to cause the fitting failure. Although the fitting effect could be improved through error in variables (EIV), observation value weighting, and other strategies, many prior assumptions of such methods were often untenable. In practice, such as equal precision of all measurement points, weight determined by reflection intensity or incident angle, etc. [53][54][55][56][57]. We find the solution from the perspective of probability theory not only relies on the point cloud data but is based on considering the geometric characteristics of the target sphere, through the global optimal parameter estimation to find the best fitting results. This method not only avoids the defects of the least square fitting method but also overcomes the influence of various noises on the fitting accuracy. From the experimental results, when the coverage rate reaches more than 30%, regardless of whether there was noise in the point cloud data, the proposed algorithm could achieve a fitting accuracy of more than 0.01 mm, which was beyond the reach of all current least-square fitting methods. When the coverage rate was less than 20%, the fitting accuracy of the proposed algorithm was reduced to some extent, but it could reach the accuracy of about 1 mm in point cloud fitting with mixed noise. In the verification process, we also conducted fitting tests on point clouds with large noise such as 10 and 15 mm, and found that the proposed algorithm could achieve fitting accuracy of more than 1 mm as long as the coverage rate remained above 20%. This fully proves that the proposed algorithm was practicable and could be used for fitting different sphere target point clouds. Figure 10. Detailed fitting process of experimental data: (a-e) were the total error changes during the fitting process of five simulated noise-free point clouds with different coverage rates, (f-j) were the total error changes during the fitting process of five simulated noisy point clouds with different coverage rates.
The sphere target fitting method based on the least-squares completely depends on the target sphere point cloud data, and the noise in the data was easy to cause the illconditioned coefficient matrix, which leads to the decrease in fitting accuracy, especially in the case of low coverage rate, it was easy to cause the fitting failure. Although the fitting effect could be improved through error in variables (EIV), observation value weighting, and other strategies, many prior assumptions of such methods were often untenable. In practice, such as equal precision of all measurement points, weight determined by reflection intensity or incident angle, etc. [53][54][55][56][57]. We find the solution from the perspective of probability theory not only relies on the point cloud data but is based on considering the geometric characteristics of the target sphere, through the global optimal parameter estimation to find the best fitting results. This method not only avoids the defects of the least square fitting method but also overcomes the influence of various noises on the fitting accuracy. From the experimental results, when the coverage rate reaches more than 30%, regardless of whether there was noise in the point cloud data, the proposed algorithm could achieve a fitting accuracy of more than 0.01 mm, which was beyond the reach of all current least-square fitting methods. When the coverage rate was less than 20%, the fitting accuracy of the proposed algorithm was reduced to some extent, but it could reach the accuracy of about 1 mm in point cloud fitting with mixed noise. In the verification process, we also conducted fitting tests on point clouds with large noise such as 10 and 15 mm, and found that the proposed algorithm could achieve fitting accuracy of more than 1 mm as long as the coverage rate remained above 20%. This fully proves that the proposed algorithm was practicable and could be used for fitting different sphere target point clouds.

Conclusions
Combined with the target sphere's point cloud and geometric characteristics, a finite random search algorithm was proposed for fast calculation of the center of the target sphere. This algorithm was suitable for all types of sphere target point cloud data and had high fitting accuracy and operation efficiency. According to the experimental results, when the coverage rate reaches more than 20%, regardless of whether the point cloud contains noise, the rapid fitting could be completed within 1 s, and the fitting error was less than 1 mm. When the coverage was less than 20%, the noise-free point cloud could also reach the fitting error of less than 1 mm. When the coverage rate was less than 20%, the fitting error of the noise-free point cloud could also reach within 1mm, while the fitting accuracy of the noisy point cloud was basically within 2 mm, although the accuracy of noisy point cloud decreases relatively.
Future research would focus on the problem of accurate fitting of noise point clouds with low coverage. The noise was removed continuously in the iterative optimization process through spatial clustering or error sorting methods to improve the fitting accuracy further.