An Efficient Probabilistic Registration Based on Shape Descriptor for Heritage Field Inspection

: Heritage documentation is implemented by digitally recording historical artifacts for the conservation and protection of these cultural heritage objects. As e ﬃ cient spatial data acquisition tools, laser scanners have been widely used to collect highly accurate three-dimensional (3D) point clouds without damaging the original structure and the environment. To ensure the integrity and quality of the collected data, ﬁeld inspection (i.e., on-spot checking the data quality) should be carried out to determine the need for additional measurements (i.e., extra laser scanning for areas with quality issues such as data missing and quality degradation). To facilitate inspection of all collected point clouds, especially checking the quality issues in overlaps between adjacent scans, all scans should be registered together. Thus, a point cloud registration method that is able to register scans fast and robustly is required. To fulﬁll the aim, this study proposes an e ﬃ cient probabilistic registration for free-form cultural heritage objects by integrating the proposed principal direction descriptor and curve constraints. We developed a novel shape descriptor based on a local frame of principal directions. Within the frame, its density and distance feature images were generated to describe the shape of the local surface. We then embedded the descriptor into a probabilistic framework to reject ambiguous matches. Spatial curves were integrated as constraints to delimit the solution space. Finally, a multi-view registration was used to reﬁne the position and orientation of each scan for the ﬁeld inspection. Comprehensive experiments show that the proposed method was able to perform well in terms of rotation error, translation error, robustness, and runtime and outperformed some commonly used approaches. indicating that the accuracy and robustness of the method are not dependent on the curves. The curves can be used to improve efficiency but do not affect robustness. Considering the results, the method shows good robustness, as it is able to exploit the advantages of the probabilistic method.


Introduction
Precise documentation of cultural heritage objects is essential for their protection and scientific studies during the restoration and renovation process. To meet the requirements for these applications, highly accurate 3D surface information is required [1]. Lots of measurement systems and technologies are used to collect the information [2][3][4][5]. Among them, laser scanner is an efficient tool to collect accurate 3D point clouds, widely used for cultural heritage objects without damaging their original match the surface geometry and intensity information. Although image-integrated methods perform well, they are susceptible to the degradation of images (e.g., distortion, illumination).
The laser point-based methods can be further divided into coarse registration and fine registration according to the registration errors. The aim of coarse registration is to estimate an initial transformation between two point clouds. Most coarse registration methods are based on geometric primitives (e.g., feature point, linear and planar feature). In this category, the feature point-based methods are more popular due to their flexibility and broad applicability. Generally, they use 3D detectors (e.g., 3D Harris [31], 3D SIFT [32], 3D SURF [33], curvature and curvature change-based [34,35], heat kernel signature (HKS)-based [36], local surface patches (LSP)-based [37], and Laplace Beltrami scale space (LBSS)-based [38] detectors) to extract key points from the original point clouds, and use the above 3D descriptors to measure the similarity to determine the correspondences. For example, Hussnain et al. used an adaptive variant of Harris-operator to detect corner points and used LATCH binary descriptor to describe them. Finally, the relative Euclidean-distance and angles between sets of points were exploited to match them [39]. Petricek et al. [40] proposed two types of detectors based on the covariance matrix of points and normal respectively. They computed the saliency based on the eigenvalues to extract key-points, and constructed the descriptor based on local reference frames to determine the correspondences for registration. He et al. [41] estimated the curvature of each point based on its eigenvalues to extract key-points. Then Tsallis entropy of spin image and the reflection intensity of laser sensor were combined to match key points to realize the registration. There are many detector and descriptor-based methods proposed to register TLS point clouds. For instance, Dong et al. calculated a local descriptor-binary shape context (BSC)-and a global descriptor-vector of locally aggregated descriptor (VLAD)-to describe the local and global characteristics of each point. Then an adjacent graph is formed to register multiple point clouds iteratively [42]. Bueno et al. calculated an eigenvector-based descriptor that describes the linear, planar and scattered characteristics of local surface. The change of curvature is another point descriptor used. These two 3D descriptors are used to detect key points and then a 4 Points Congruent Sets registration is performed [43]. Many methods have been used to register cultural heritage. For example, Shanoer et al. realized the registration by the manual selection of at least three corresponding points [44]. Tournas et al. proposed a target point-based registration for cultural heritage documentation. These methods have better accuracy, but they are labor intensive and time consuming [45]. In contrast, without any targets, Shao et al. used Super4PCS (i.e., randomly selecting four coplanar points and matching all the 4-point configurations with similar geometric properties) algorithm to realize coarse registration [46]. Bae et al. proposed to use the change of geometric curvature and approximate normal vector to search for possible corresponding points for registration. This method was applied for deformation monitoring of large Buddha [47]. These methods work well for specific applications, but for the filed inspection, a more stable and efficient registration method is still required.
Linear and planar feature-based methods have also attracted interest from researchers [48]. Among this category, canny edge detection [49,50], Hough transform [51,52], line segment detector [53] were extended to detect linear features in 3D point clouds. For example, Date et al. [54] used straight lines to register point clouds for bridges. They extracted horizontal and vertical lines based on planar regions and used two pairs of nonparallel lines to solve the transformation. Similarly, Cui et al. [55] applied linear features to register the panoramic images and LiDAR point clouds collected by a mobile mapping system (MMS). Moreover, many methods have been proposed to detect planar features. For instance, Vosselman and Dijkman used a 3D version of the Hough transform to extract roof planes for 3D reconstruction [56]. Vosselman et al. proposed to group neighborhood points that share the same property (e.g., surface normal) to detect smooth surfaces from scan line datasets and improved the surfaces by extending the split-and-merge method in image [57]. Pu and Vosselman proposed knowledge-based feature constraints to detect linear and planar features of building façade (e.g., walls, doors, roofs) from laser point clouds. Then concave polygons were fitted to reconstruct these components [58]. Khoshelham [59] used a normalization algorithm to fit planes. A closed-form solution based on point-plane correspondences was proposed to determine the transformation. Similar works can also be found in [60,61]. Linear and planar feature-based methods perform well for urban and industrial scenes with lots of regular structures but operate poorly for free-form objects with less line/plane features.
Most recently, deep learning-based methods have been extensively explored. These techniques learn high-level feature information from a large volume of data for registration. For example, Zhang et al. [62] divided the point clouds into regular grids and used a 3D convolution network to extract features for matching. Pujol-Miró et al. [63] proposed a multiview-based method, which generates multiple views from different angles using a single point cloud. Image convolutional neural networks have been used to determine the correspondences between views. PointNet [64] and PointNet++ [65] are pioneering methods that directly process unordered point clouds. T-Net and different feeding orders have also been proposed to make the process transform-invariant. Similarly, researchers have proposed PPFNet [66], PointNetLK [67], and CorsNet [68] for point cloud registration. Other commonly used methods include the random sample consensus (RANSAC) [69], the lifting optimization method [70], the branch-and-bound (BnB)-based method [71], and the simultaneous localization and mapping (SLAM) method [72].
Within the category of coarse registration methods, the probabilistic registration method is a well-known and promising approach. Different from the traditional one-to-one matching strategy, one-to-many matching patterns are allowed. More combinations of correspondences are considered to improve the flexibility and robustness of the technique. Most of the previous studies were based on the Expectation-Maximization (EM) algorithm. For example, Evangelidis et al. [73] proposed a joint registration of multiple point clouds (JR-MPC). They used a Gaussian mixture model (GMM) to describe the distribution of one point cloud and assigned points from the other point cloud to the GMM. The EM algorithm was then used to estimate both the GMM and the transformation parameters. Myronenko et al. [74] introduced the coherent point drift (CPD) algorithm for the point cloud registration, using GMM centroids to describe one point cloud and fitting it to the other point cloud by maximizing a likelihood function. The use of the CPD algorithm was able to generate accurate and robust results. The variants of CPD include the automatic estimation of outliers [75] and computational efficiency improvements [76].
After the initial transformation, fine registration can then be used to further refine the transformation. The iterative closest point (ICP) algorithm [77,78] is a well-known fine registration method that matches correspondence by searching the nearest point and solves the transformation parameters by minimizing the distance errors iteratively. Because the original ICP algorithm has many limitations, many variants have been proposed, mainly focusing on improving computational expense, correspondence establishment, and robustness [79,80]. The normal distribution transform (NDT) algorithm is another widely used fine registration method. It was first used for robotic localization in 2D space before it was extended into 3D space [81]. NDT uses multiple Gaussian distributions to describe different point clouds and estimate the transformation parameters by maximizing the similarity between Gaussian distributions. The variants of NDT mainly improve the computational time and the convergence basin [82,83].

Processing Pipeline of Digital Documentation
To better show the research content of this article. The whole processing pipeline of our digital documental is shown in Figure 1. The pipeline consists of the following steps, namely planning of scanning, preparation of protective measures, measures for scanning preparation, field inspection, accurate registration, accurate modeling and digital documentation. maximizing the similarity between Gaussian distributions. The variants of NDT mainly improve the computational time and the convergence basin [82,83]. (1) Planning of scanning. In this step, the number of stations and the scanning area of each station were planned. Take the data acquisition in Mogao Grottoes for example, the following steps were employed. (a) Take photos covering the operation area. According to the characteristics of the scene, a scanning order was determined (as Figure 2 shows). (b) For each Buddha, its shape and surroundings were taken into account to determine the number of stations and the scanning area of each station. Generally, for the Buddha with complex shape, more stations were planned to ensure the data quality, and vise verse. Moreover, the space needed for the target net and its fixed support should also be taken into account. (c) To facilitate recording and operation, the planned stations were named after the parts of body (as Figure 3 shows). To better show the research content of this article. The whole processing pipeline of our digital documental is shown in Figure 1. The pipeline consists of the following steps, namely planning of scanning, preparation of protective measures, measures for scanning preparation, field inspection, accurate registration, accurate modeling and digital documentation.

Processing Pipeline of Digital Documentation
(1) Planning of scanning. In this step, the number of stations and the scanning area of each station were planned. Take the data acquisition in Mogao Grottoes for example, the following steps were employed. (a) Take photos covering the operation area. According to the characteristics of the scene, a scanning order was determined (as Figure 2 shows). (b) For each Buddha, its shape and surroundings were taken into account to determine the number of stations and the scanning area of each station. Generally, for the Buddha with complex shape, more stations were planned to ensure the data quality, and vise verse. Moreover, the space needed for the target net and its fixed support should also be taken into account. (c) To facilitate recording and operation, the planned stations were named after the parts of body (as Figure 3 shows).  (2) Preparation of protective measures. For the effective protection of cultural heritage objects, many necessary facilities (e.g., glass fender, air monitor, protective blanket) were prepared and distributed.
(3) Measures for scanning preparation. This step prepares the facilities (e.g., power unit, lighting equipment, work shelf) for scanning operation. For each station, a special target net was made for scanning, which should be close to the shape of scanning area. To protect the historical artifacts, the target nets were not allowed to touch them. Thus, a fixed support used for hanging and fixing them  To better show the research content of this article. The whole processing pipeline of our digital documental is shown in Figure 1. The pipeline consists of the following steps, namely planning of scanning, preparation of protective measures, measures for scanning preparation, field inspection, accurate registration, accurate modeling and digital documentation.
(1) Planning of scanning. In this step, the number of stations and the scanning area of each station were planned. Take the data acquisition in Mogao Grottoes for example, the following steps were employed. (a) Take photos covering the operation area. According to the characteristics of the scene, a scanning order was determined (as Figure 2 shows). (b) For each Buddha, its shape and surroundings were taken into account to determine the number of stations and the scanning area of each station. Generally, for the Buddha with complex shape, more stations were planned to ensure the data quality, and vise verse. Moreover, the space needed for the target net and its fixed support should also be taken into account. (c) To facilitate recording and operation, the planned stations were named after the parts of body (as Figure 3 shows).  (2) Preparation of protective measures. For the effective protection of cultural heritage objects, many necessary facilities (e.g., glass fender, air monitor, protective blanket) were prepared and distributed.
(3) Measures for scanning preparation. This step prepares the facilities (e.g., power unit, lighting equipment, work shelf) for scanning operation. For each station, a special target net was made for scanning, which should be close to the shape of scanning area. To protect the historical artifacts, the target nets were not allowed to touch them. Thus, a fixed support used for hanging and fixing them (2) Preparation of protective measures. For the effective protection of cultural heritage objects, many necessary facilities (e.g., glass fender, air monitor, protective blanket) were prepared and distributed.
(3) Measures for scanning preparation. This step prepares the facilities (e.g., power unit, lighting equipment, work shelf) for scanning operation. For each station, a special target net was made for scanning, which should be close to the shape of scanning area. To protect the historical artifacts, the target nets were not allowed to touch them. Thus, a fixed support used for hanging and fixing them is needed and designed for each station (as Figure 4 shows).
(4) Field inspection. To ensure the data quality and save on labor costs, this step uses a rapid registration for field inspection to determine the need for extra scanning. Even though detailed plans were designed in advance, long hours of scanning and limited working space challenge the data collector's endurance, energy and concentration. In this case, signal loss, target net movement and resolution variations caused by manual errors are difficult to avoid, which will further lead to data degradation issues. Therefore, field inspection is required in the practical applications.
(5) Accurate post-registration in interior work. In this step, data preprocessing (e.g., noise removal, down sampling, duplicates removal) is carried out. Then a manual registration and a fine registration is used to get high accurate alignment (submillimeter level) for accurate modeling.
(6) Accurate modeling. Based on the accurate alignment, this step is used to construct tin-based meshes, which includes many model optimization procedures (e.g., holes filling, bridging, rendering by normal vectors, nails removal).
(7) Digital documentation. This step stores the original data, data acquisition diary, accurate 3D models, description files, and the final report.   (4) Field inspection. To ensure the data quality and save on labor costs, this step uses a rapid registration for field inspection to determine the need for extra scanning. Even though detailed plans were designed in advance, long hours of scanning and limited working space challenge the data collector's endurance, energy and concentration. In this case, signal loss, target net movement and resolution variations caused by manual errors are difficult to avoid, which will further lead to data degradation issues. Therefore, field inspection is required in the practical applications.

Experimental Point Cloud Datasets and Work Flow
(5) Accurate post-registration in interior work. In this step, data preprocessing (e.g., noise removal, down sampling, duplicates removal) is carried out. Then a manual registration and a fine registration is used to get high accurate alignment (submillimeter level) for accurate modeling.
(6) Accurate modeling. Based on the accurate alignment, this step is used to construct tin-based meshes, which includes many model optimization procedures (e.g., holes filling, bridging, rendering by normal vectors, nails removal).
(7) Digital documentation. This step stores the original data, data acquisition diary, accurate 3D models, description files, and the final report.

Experimental Point Cloud Datasets and Work Flow
The proposed method was used on scan data from the Mo Kao Grotto at Dunhuang and Wuwei museum in evaluating its performance. All the point clouds were collected using HandySCAN 3D (a hand-held scanner from Creaform, Lévis, Quebec, Canada). Before scanning, a target net was set up to help position the laser beam cross of each moment. Adjacent scans had different overlaps and initial positions. Scans over the General from the Mo Kao Grotto and the Buddha from the Wuwei museum were used for field inspection. We collected the General data in July 2012 and the Buddha data in October 2014. Six scans were set up for each, and 0.67 and 1.44 million points were collected separately. The average point density for the General was about 1.0 mm, and 0.8 mm for the Buddha. The General was 1.6 m in height and 0.9 m in width, while the Buddha in Wuwei museum was 1.4 m in height and 1.2 m in width. The illustration of the scanned point clouds for the General and the Buddha is shown in Figure 5.   The proposed method mainly includes three parts: the descriptor construction and similarity measurement, the construction of probabilistic framework encoding curve constraints, and the multi-view adjustment method and inspection. The workflow is shown in Figure 6. The proposed method was used on scan data from the Mo Kao Grotto at Dunhuang and Wuwei museum in evaluating its performance. All the point clouds were collected using HandySCAN 3D (a hand-held scanner from Creaform, Lévis, Quebec, Canada). Before scanning, a target net was set up to help position the laser beam cross of each moment. Adjacent scans had different overlaps and initial positions. Scans over the General from the Mo Kao Grotto and the Buddha from the Wuwei museum were used for field inspection. We collected the General data in July 2012 and the Buddha data in October 2014. Six scans were set up for each, and 0.67 and 1.44 million points were collected separately. The average point density for the General was about 1.0 mm, and 0.8 mm for the Buddha. The General was 1.6 m in height and 0.9 m in width, while the Buddha in Wuwei museum was 1.4 m in height and 1.2 m in width. The illustration of the scanned point clouds for the General and the Buddha is shown in Figure 5.

Experimental Point Cloud Datasets and Work Flow
The proposed method mainly includes three parts: the descriptor construction and similarity measurement, the construction of probabilistic framework encoding curve constraints, and the multiview adjustment method and inspection. The workflow is shown in Figure 6.

Problem Formulation
To obtain the complete data of cultural heritage objects for field inspection, pairwise registration was employed as a fundamental kernel problem. To make the method easy to understand, we formulate the general problem of pairwise registration as follows, given two overlapping point clouds P and Q. Assume there are M correspondences between them given by The expression to estimate the rigid transformation parameters R and T is as follows:

Problem Formulation
To obtain the complete data of cultural heritage objects for field inspection, pairwise registration was employed as a fundamental kernel problem. To make the method easy to understand, we formulate the general problem of pairwise registration as follows, given two overlapping point clouds P and Q. Assume there are M correspondences between them given by C = (p i , q i ) i, (1, 2, . . . , M) , where (p i , q i ) denotes one correspondence. The expression to estimate the rigid transformation parameters R and T is as follows: where Tr(·) denotes the transformation; R(θ) is a 3 × 3 rotation matrix; and, T is a 3 × 1 translation vector. To estimate the transformation parameters, this study constructed an objective function that encodes the candidate correspondences assigned by the matching probability matrix between P and Q: The transformation parameters and covariance are solved in an iterative way within a probabilistic framework.

Principal Direction Descriptor of Local Surface
Correspondence establishment between different scans is a precondition of registration. Considering the requirements for field inspection, a fast and robust 3D descriptor tailored for cultural heritage objects (e.g., a free-form object like Buddha) is needed. This section introduces the proposed principal direction descriptor of the local surface.

Construction of Local Frame
Curvature quantifies curvedness and indicates the undulation of the surface, which is important in various applications (e.g., shape analysis, object recognition, anisotropic texture mapping). Curvature is characterized by two perpendicular principal directions (i.e., maximum and minimum principal directions), encoding abundant geometric information of the local surface. These geometric data are intrinsic properties that are invariant to rigid body transformation. Along the principal directions, the surface bends the most. These characteristics make them suitable as the local frame for the descriptors. However, the complex estimation algorithm limits their application.
The principal directions indicate the courses with the steepest and slowest surface change within the local surface. Instead of using complex parametric quadratic surface fitting, we propose a fast and simple estimation method. Assume that p 0 is one point from the point cloud, N P indicates its neighborhood point set with N neighbors within a radius r, and p i is one of the neighbors. Let − − → n 0 be the normal vector of p 0 , estimated by the principal components analysis (PCA) algorithm, and the direction is unified by the consistent propagation algorithm. To ensure robustness, the point density and distance of each neighbor are also considered. T p is the tangent plane of p 0 . The maximum principal direction is estimated as: Here, − −−− → p 0 p i is the projected unit vector from p 0 to p i within T p , which can be written as: w d i is the weight of point density of p i , w s i is the weight of slope of p i , written as: Here, N i is the number of neighborhood points of p i within radius r. The point density weight implies that the sparser neighborhood point set is, the larger the weight, and vice versa. HD i is the L2 λ is the coefficient that balances the point density weight and the slope weight (suggest λ = 0.5). w s i implies that the steeper the value of p i , the larger the weight, and vice versa.
The minimum principal direction is estimated as: Here × is the cross product. Then, a local frame, with − − → ψ v as u and v axis and p i as the origin, is constructed.

Generation of Descriptor Images and Similarity Measurement
To improve the efficiency of probabilistic registration, a sampling method [84] is used to sample each scan to about 2000 points. For each sampled point, its neighborhood points of raw point cloud (about 400 points are suggested) is used to generate the descriptor images: density and distance feature images. We transform the neighborhood point set of p i into the local frame to obtain their local coordinates, written as: These transformed points are divided into M × M grids (M = 27 is suggested). To better describe the shape of local surface, the size of each grid is set about the mean point span. Each grid records the projection distance (from the neighbor to the tangent plane) of each neighbor and the number of points falling within. Next, we calculate the distance and density features of each grid to form shape images. To reduce the noise sensitivity, only points whose distance is within three times the Root Mean Square Error (RMSE) are kept. The distance feature f dis tan ce Here, max(density) and min(density) are the maximum and minimum density values, while max(dis tan ce) and min(dis tan ce) are the maximum and minimum distance values from all grids, as shown in Figure 7.
These transformed points are divided into M×M grids (M = 27 is suggested). To better describe the shape of local surface, the size of each grid is set about the mean point span. Each grid records the projection distance (from the neighbor to the tangent plane) of each neighbor and the number of points falling within. Next, we calculate the distance and density features of each grid to form shape images. To reduce the noise sensitivity, only points whose distance is within three times the Root Mean Square Error (RMSE) are kept. The distance feature Here, λ and γ are the weights that balance the density term and the distance term (suggest:  For the two density images a and b, the similarity S density (a, b) between them is measured according to the correlation coefficient. A similar procedure is used for the similarity between two distance images S dis tan ce (a, b). The final similarity between the two descriptors is calculated by combining the similarities: S(a, b) = λS density (i, j) + γS dis tan ce (i, j), s.t. λ + γ = 1.0 (9) Here, λ and γ are the weights that balance the density term and the distance term (suggest: λ = 0.3).

Spatial Curves Extraction from a Free-Form Object
Spatial curves are common and important geometric primitives in free-form artifact objects (e.g., Buddha). To improve robustness, our developed method [23] was used to extract spatial curves. In the algorithm, curvature and several geometric constraints are combined to extract complete and accurate curves. A brief diagram is shown in Figure 8, illustrating how we made our method self-contained.

Spatial Curves Extraction from a Free-Form Object
Spatial curves are common and important geometric primitives in free-form artifact objects (e.g., Buddha). To improve robustness, our developed method [23] was used to extract spatial curves. In the algorithm, curvature and several geometric constraints are combined to extract complete and accurate curves. A brief diagram is shown in Figure 8, illustrating how we made our method selfcontained.

Improved Pairwise Probabilistic Registration
Probabilistic registration treats the registration problem as a GMM's (Gaussian Mixture Model) probability density estimation problem. Based on the original probabilistic framework of registration [74], this section embeds the principal direction descriptor into the posterior matching probability and integrates the curve constraints to improve performance. To make the discussion concise, we only present the improvements in this section.

Probabilistic Registration Based on Principal Direction Descriptor
Take one point cloud Y3×M as the centroid of the component of a GMM, and the other point cloud X3×N is assumed to be generated by these components. For each point pair, a matching probability is assigned. If the two-point clouds are optimally aligned, the sum of the probabilities should be maximum. Registration can be realized by minimizing the objective function, written as: whereθ represents the transformation, and 2 δ is the covariance controlling the motion of each iteration. These variables are unknowns that will be determined by the EM algorithm. E-step

Improved Pairwise Probabilistic Registration
Probabilistic registration treats the registration problem as a GMM's (Gaussian Mixture Model) probability density estimation problem. Based on the original probabilistic framework of registration [74], this section embeds the principal direction descriptor into the posterior matching probability and integrates the curve constraints to improve performance. To make the discussion concise, we only present the improvements in this section.

Probabilistic Registration Based on Principal Direction Descriptor
Take one point cloud Y 3×M as the centroid of the component of a GMM, and the other point cloud X 3×N is assumed to be generated by these components. For each point pair, a matching probability is assigned. If the two-point clouds are optimally aligned, the sum of the probabilities should be maximum. Registration can be realized by minimizing the objective function, written as: where θ represents the transformation, and δ 2 is the covariance controlling the motion of each iteration. These variables are unknowns that will be determined by the EM algorithm. E-step estimates the transformation parameters and uses Bayes' theorem to calculate the posterior matching probability, while M-step maximizes the lower bound of Equation 10 to estimate θ and δ 2 . Tr(·) indicates the transformation of the current iteration. N p is the summation of matching probability P entries, D = 3 is the dimension of a 3D point. P old (m x n ) is the posterior probability that indicates the correspondence determination. We embed the principal direction descriptor in the posterior matching probability as: Here, the Euclidean distance and the similarity of the descriptor are combined to calculate the ρ(x n , y m ), written as: where w 1 is the weight coefficient that balances the Euclidean distance term and the descriptor term, and w indicates the amount of noise and outliers (w = 0.2 is suggested). All posterior probabilities form a M × N matching probability matrix P. The E-step and M-step are iterated to obtain the unknowns.

Spatial Curve Constraints
This section embeds the curve constraints into the probabilistic framework to further improve the robustness and effectiveness. Assume that a candidate matching pair is x n and y m . If only one of them (x n or y m ) belongs to a spatial curve, the matching probability is determined by considering the Euclidean distance (shown in Equation 13). If both of them belong to a spatial curve, the similarity of neighborhood distribution is also considered within the matching probability and is written as: 2δ 2 x n − Tr(y m ) 2 , i f x n or y m is curve point − 1 2δ 2 x n − Tr(y m ) 2 + w 1 · S(m, n) + w 2 · Q(m, n), i f both x n and y m are curve points (13) Here, S(m, n) represents the similarity between two descriptors of x n and y m , and Q(m, n) is the integrated similarity of curve points within three neighborhood scales. If this candidate matching pair (x n and y m ) are true correspondence, this term will be large. w 1 and w 2 are the weight coefficients that balance the three terms in Equation (13). If neither x n and y m are curve points, we can calculate their probability according to Equation (11). Equation (13) implies that if this candidate matching pair (x n and y m ) are not true correspondence, their matching probability will be low, and vice versa. Q(m, n) is written as: Here, N 1 , N 2 , and N 3 are the number of curve points from three scales (with different neighborhood radii); x i is one neighborhood curve point of x n ; and, y j is its nearest point after the current transformation.
We then calculate ρ x i , y j according to Equation (12). Equation (14) demonstrates that for one true correspondence, their neighborhood curve points should also be similar within three scales (as shown in Figure 9). This correspondence would be given a large matching probability.
Here, N1, N2, and N3 are the number of curve points from three scales (with different neighborhood radii); xi is one neighborhood curve point of xn; and, yj is its nearest point after the current transformation. We then calculate ( ) , i j x y ρ according to Equation (12). Equation (14) demonstrates that for one true correspondence, their neighborhood curve points should also be similar within three scales (as shown in Figure 9). This correspondence would be given a large matching probability. Figure 9. Schematic diagram of curve constraints: red circles indicate neighborhood areas of three scales from two point clouds.

Multi-Scan Registration and Field Inspection
Based on the pairwise registration results, we constructed an undirected graph to describe the topology between scans. Each node indicates one scan, while the edge connecting two nodes indicates the two scans with overlap. For one pair of scans S 1 and S 2 , we determined each point's (in S 1 ) nearest point in S 2 to calculate the overlapping ratio. For example, if the overlapping ratio between S 1 and S 2 is 30%, which is larger than the 5% threshold, a new edge is generated to connect them. The weight of this edge is given by the overlapping ratio. A graph was constructed based on all scans. We then used the adjustment method in [85] to refine the position of each scan.
After obtaining a complete point cloud of one artifact object, field inspection was carried out to determine if there are places with missing or quality degraded data. Extra scanning work was executed based on the inspection results. An image of actual fieldwork activity and the software interface for inspection are shown in Figure 10.

Multi-Scan Registration and Field Inspection
Based on the pairwise registration results, we constructed an undirected graph to describe the topology between scans. Each node indicates one scan, while the edge connecting two nodes indicates the two scans with overlap. For one pair of scans S1 and S2, we determined each point's (in S1) nearest point in S2 to calculate the overlapping ratio. For example, if the overlapping ratio between S1 and S2 is 30%, which is larger than the 5% threshold, a new edge is generated to connect them. The weight of this edge is given by the overlapping ratio. A graph was constructed based on all scans. We then used the adjustment method in [85] to refine the position of each scan. After obtaining a complete point cloud of one artifact object, field inspection was carried out to determine if there are places with missing or quality degraded data. Extra scanning work was executed based on the inspection results. An image of actual fieldwork activity and the software interface for inspection are shown in Figure 10.

Results
To demonstrate the performance of the proposed method, we carried out experiments based on the datasets introduced in Section 2. The implementation details of the experiments are described in this section.

Correspondence Establishment
For a probabilistic registration method, the transformation solved after the first iteration is

Results
To demonstrate the performance of the proposed method, we carried out experiments based on the datasets introduced in Section 2. The implementation details of the experiments are described in this section.

Correspondence Establishment
For a probabilistic registration method, the transformation solved after the first iteration is important. Good first transformation can efficiently reduce the solution space and improve registration efficiency. Otherwise, it requires more iterations and easily falls into a local minimum. To determine a good first transformation, well matched correspondences should be input for the first iteration. For this purpose, we developed a matching mechanism that combines the descriptor-based matching strategy and curve constraint based rejection scheme for the suppression of incorrect matches. To demonstrate the matched correspondences for the first iteration, we selected the top 2000 pairs of correspondences based on the matching probability matrix. According to the number of incorrect correspondences, 40 matches were randomly selected and are shown in Figure 11. As shown in Figure 11, incorrect ratios of matches from the first iteration were measured at 18.5%, 20%, 21.4%, and 15.8%. For a probabilistic method, this amount of incorrect matches is considered acceptable, which means a good transformation can be obtained. Moreover, these experimental data represent the common appearances of free-form cultural heritage, some even have obvious curved surface (i.e., Figure 11 (c)). The mean incorrect ratios of matches is controlled within 20%, this is because the local frame of the principal direction descriptor is suitable for describing the local geometric properties. These results demonstrate the distinctiveness and descriptiveness of the proposed descriptor and the effectiveness of the matching mechanism. Note that during the iterations, the transformed point clouds gradually shift towards the correct position, and more and more correct correspondences can be matched.

Registration Results
To improve the effectiveness of registration, we used the Geomagic Studio software (Geomagic, Chapelhill, NC, USA) to sample each scan uniformly to a 3.0 mm average point density (about 6000 points). The registration accuracy of the General and the Buddha data is shown in Figure 12. The color is rendered using the distance from one point in the target point cloud to its nearest point of the transformed template point cloud in the overlap areas. To qualitatively demonstrate the registration errors, the results of manual registration were used to evaluate the registration accuracy. Manual registration was performed by manually selecting the corresponding points from the adjacent scans, followed by the ICP refinement. Table 1 provides the evaluation results of the registration performance. As shown in Figure 11, incorrect ratios of matches from the first iteration were measured at 18.5%, 20%, 21.4%, and 15.8%. For a probabilistic method, this amount of incorrect matches is considered acceptable, which means a good transformation can be obtained. Moreover, these experimental data represent the common appearances of free-form cultural heritage, some even have obvious curved surface (i.e., Figure 11c). The mean incorrect ratios of matches is controlled within 20%, this is because the local frame of the principal direction descriptor is suitable for describing the local geometric properties. These results demonstrate the distinctiveness and descriptiveness of the proposed descriptor and the effectiveness of the matching mechanism. Note that during the iterations, the transformed point clouds gradually shift towards the correct position, and more and more correct correspondences can be matched.

Registration Results
To improve the effectiveness of registration, we used the Geomagic Studio software (Geomagic, Chapelhill, NC, USA) to sample each scan uniformly to a 3.0 mm average point density (about 6000 points). The registration accuracy of the General and the Buddha data is shown in Figure 12. The color is rendered using the distance from one point in the target point cloud to its nearest point of the transformed template point cloud in the overlap areas. To qualitatively demonstrate the registration errors, the results of manual registration were used to evaluate the registration accuracy. Manual registration was performed by manually selecting the corresponding points from the adjacent scans, followed by the ICP refinement. Table 1 provides the evaluation results of the registration performance. color is rendered using the distance from one point in the target point cloud to its nearest point of the transformed template point cloud in the overlap areas. To qualitatively demonstrate the registration errors, the results of manual registration were used to evaluate the registration accuracy. Manual registration was performed by manually selecting the corresponding points from the adjacent scans, followed by the ICP refinement. Table 1 provides the evaluation results of the registration performance.   These scans have different overlaps, varying levels of noise, and different initial positions and orientations. Figure 12 and Table 1 indicate that all the pairs of scans do not fall into local convergence and are registered successfully. Column Std. from Table 1 shows that their mean registration errors are within 2.0 mm, demonstrating a good distribution of registration errors. Although several scans have few curves (for example, General data scan 3&4, and 4&5 in Figure 7a) or limited overlaps (for example, 21% for General scan 1&2, and 23% for Buddha scan 2&5), they were still able to get good results (about 0.5 mm mean error for the General data and about 1.5 mm mean error for Buddha data). The above results suggest that the proposed method has good stability. By statistical analysis, considering the actual size of cultural heritage object (about 1.5 m in height and 1.0m in width), the mean errors and std. values are sufficient for the application of field inspection.
Note that the scan 5&6 for the General and scan 2&7 for the Buddha have relatively large errors (1.21 and 2.43 mm, respectively), which were probably caused by repetitive structures resulting in ambiguous correspondences and affecting registration accuracy. The registration errors of the multi-view method are 0.38 mm and 0.69 mm, demonstrating the method's effectiveness. In addition, Table 1 shows that the number of iterations needed for convergence is within 10, and the runtime is about 1 minute for each pair of scans. These results suggest that our approach method provides a very fast and efficient technique that is tailored and suitable for field inspection.

Field Inspection Result
To ensure the reliability of the inspection, we performed a multi-view method to further refine the position of each scan. Figure 13 shows the generated images for the General and the Buddha after undergoing the multi-view method. The red circles indicate the problem places identified by inspection. As shown in the figure, many geometric details (e.g., defects) have been preserved well. Through field inspection, we are able to identify the places with missing data (for example, areas A and B in General data, and areas C and D in Buddha data) and those with poor quality (for example, area E in Buddha data). To satisfy the requirements for high accurate heritage documentation, these areas would require extra field scanning.

Pairwise Registration Comparison
To further evaluate the method's performance, we compared the proposed method to other pairwise registration methods, including the popular RANSAC method [69], the BnB method [71], and the LM method (lifting method) [70]. The codes can be obtained from the link (https://github.com/ZhipengCai/Demo---Practical-optimal-registration-of-terrestrial-LiDAR-scanpairs). The following parameters of the compared methods are tuned to get best performance. Specifically, in RANSAC method, the probability of searching a valid sample set is 0.98; In BnB method, 300×300×300 grids are used for 3D Distance Transform (DT) computation, and the convergence threshold is set to 0.01; In LM method, the annealing rate is set to 1.2. The transformation parameters of manual registration were used to evaluate the proposed method. To ensure the reliability, manual registration was performed by carefully selecting some corresponding points from adjacent scans, and the registration result was further refined by the ICP method. The rotation and translation errors were calculated by the deviations from the transformation parameters of manual

Pairwise Registration Comparison
To further evaluate the method's performance, we compared the proposed method to other pairwise registration methods, including the popular RANSAC method [69], the BnB method [71], and the LM method (lifting method) [70]. The codes can be obtained from the link (https://github.com/ ZhipengCai/Demo---Practical-optimal-registration-of-terrestrial-LiDAR-scan-pairs). The following parameters of the compared methods are tuned to get best performance. Specifically, in RANSAC method, the probability of searching a valid sample set is 0.98; In BnB method, 300×300×300 grids are used for 3D Distance Transform (DT) computation, and the convergence threshold is set to 0.01; In LM method, the annealing rate is set to 1.2. The transformation parameters of manual registration were used to evaluate the proposed method. To ensure the reliability, manual registration was performed by carefully selecting some corresponding points from adjacent scans, and the registration result was further refined by the ICP method. The rotation and translation errors were calculated by the deviations from the transformation parameters of manual registration using: ∆A = ∆ϕ 2 + ∆ω 2 + ∆γ 2 and Figure 14 shows the translation and rotation errors and the runtime for each method. The experiments were implemented on a computer with 16 GB RAM and an Intel Core i7-6700HQ @2.60 GHz CPU. As shown in Figures 14(a) and (b), the proposed method performed the best. LM method had a competitive performance both in rotation and translation estimation. This is because LM is also a softassign-based optimization method, but our method is specialized for the free-form cultural heritage objects that geometric properties and curve constraints are combined efficiently, leading to a better result. The RANSAC method could sometimes generate better results than the proposed method (for example, in scan 2&4 for the Buddha data in Figure 14(a); in scan 2&3 and 4&5 for the General data in Figure 14(b)). This result was probably dependent on the randomly selected primitives. The BnB method had the worst performance. Given that the BnB method's main advantage is obtaining global optimization in a large solution space, the tested scans had small initial positions and orientations, which are more suited for the proposed method. Figure 14(c) shows that the LM method performs the best among the compared methods, and the proposed approach performed competitively. The runtime for the RANSAC method was not stable (for example, scan 2&4, 2&7 of Buddha data needs the least runtime). In conclusion, the proposed method is efficient both in accuracy and in runtime. As shown in Figure 14a,b, the proposed method performed the best. LM method had a competitive performance both in rotation and translation estimation. This is because LM is also a softassign-based optimization method, but our method is specialized for the free-form cultural heritage objects that geometric properties and curve constraints are combined efficiently, leading to a better result. The RANSAC method could sometimes generate better results than the proposed method (for example, in scan 2&4 for the Buddha data in Figure 14a; in scan 2&3 and 4&5 for the General data in Figure 14b). This result was probably dependent on the randomly selected primitives. The BnB method had the worst performance. Given that the BnB method's main advantage is obtaining global optimization in a large solution space, the tested scans had small initial positions and orientations, which are more suited for the proposed method. Figure 14c shows that the LM method performs the best among the compared methods, and the proposed approach performed competitively. The runtime for the RANSAC method was not stable (for example, scan 2&4, 2&7 of Buddha data needs the least runtime). In conclusion, the proposed method is efficient both in accuracy and in runtime.

Evaluation of Robustness Performance
To test the robustness of the proposed method comprehensively, we simulated various datasets for testing. The pairwise registration results are shown in Figure 15. Figure 15a shows the registration with different amounts of noise (8%, 4%, and 2% Gaussian noise are added to the original point clouds). Figure 15b shows the registration for different overlaps (i.e., 30%, 55%, and 80%) of the original point clouds. Figure 15c shows the registration results on objects with few curves. To test the robustness of the proposed method comprehensively, we simulated various datasets for testing. The pairwise registration results are shown in Figure 15. Figure 15(a) shows the registration with different amounts of noise (8%, 4%, and 2% Gaussian noise are added to the original point clouds). Figure 15(b) shows the registration for different overlaps (i.e., 30%, 55%, and 80%) of the original point clouds. Figure 15(c) shows the registration results on objects with few curves.
As shown in Figure 15, the point clouds for all the situations were registered successfully. Figure  15(a) shows the proposed method is robust to noise and outliers. Although the middle column gets relatively poor results (i.e., about 2.0 mm), the point clouds can still be aligned together. Figure 15(b) shows the proposed method performs robustly for different overlap ratios. However, for the 30% overlap ratio (the left column in Figure 15(a)), the proposed method generated comparatively poor results, indicating that the overlap value can affect the method's accuracy. In Figure 15(c), the method performs well, indicating that the accuracy and robustness of the method are not dependent on the curves. The curves can be used to improve efficiency but do not affect robustness. Considering the results, the method shows good robustness, as it is able to exploit the advantages of the probabilistic method.

Ablation Study
To directly show the contributions of the proposed method, an ablation study was implemented, consisting of two components: the descriptor ablation and the curve constraint ablation. Descriptor ablation is used to check the role of the proposed descriptor by removing the principal direction descriptor from the method and preserving the Euclidean term and the curve constraint term. For the As shown in Figure 15, the point clouds for all the situations were registered successfully. Figure 15a shows the proposed method is robust to noise and outliers. Although the middle column gets relatively poor results (i.e., about 2.0 mm), the point clouds can still be aligned together. Figure 15b shows the proposed method performs robustly for different overlap ratios. However, for the 30% overlap ratio (the left column in Figure 15a), the proposed method generated comparatively poor results, indicating that the overlap value can affect the method's accuracy. In Figure 15c, the method performs well, indicating that the accuracy and robustness of the method are not dependent on the curves. The curves can be used to improve efficiency but do not affect robustness. Considering the results, the method shows good robustness, as it is able to exploit the advantages of the probabilistic method.

Ablation Study
To directly show the contributions of the proposed method, an ablation study was implemented, consisting of two components: the descriptor ablation and the curve constraint ablation. Descriptor ablation is used to check the role of the proposed descriptor by removing the principal direction descriptor from the method and preserving the Euclidean term and the curve constraint term. For the curve constraint ablation, the role of the curve constraints is tested by removing them. Three pairs of scans were randomly selected from the General and the Buddha data for testing. The translation and rotation errors were then calculated based on the true value of the transformation, as discussed in Section 4.1.
The results of the ablation study are shown in Figure 16.

Ablation Study
To directly show the contributions of the proposed method, an ablation study was implemented, consisting of two components: the descriptor ablation and the curve constraint ablation. Descriptor ablation is used to check the role of the proposed descriptor by removing the principal direction descriptor from the method and preserving the Euclidean term and the curve constraint term. For the curve constraint ablation, the role of the curve constraints is tested by removing them. Three pairs of scans were randomly selected from the General and the Buddha data for testing. The translation and rotation errors were then calculated based on the true value of the transformation, as discussed in Section 4.1. The results of the ablation study are shown in Figure 16. As shown in Figure 16, the results indicate that the descriptor and the curve constraints contribute significantly to the effectiveness and registration accuracy. Specifically, without the descriptor term, the method would result in considerably larger registration errors (almost three times than the proposed method) and would require more time (almost twice than the proposed method). For scans 2&4 and 2&5 of the Buddha data in Figure 16(a) and (b), even the correct convergence cannot be ensured, which could suggest that the descriptor rejects much of the incorrect matches, thereby avoiding local convergence and improving efficiency. On the other hand, without the curve constraints, the registration errors are similar to the proposed method. However, the runtime needed for convergence is significantly increased. This suggests that the curve constraints As shown in Figure 16, the results indicate that the descriptor and the curve constraints contribute significantly to the effectiveness and registration accuracy. Specifically, without the descriptor term, the method would result in considerably larger registration errors (almost three times than the proposed method) and would require more time (almost twice than the proposed method). For scans 2&4 and 2&5 of the Buddha data in Figure 16a,b, even the correct convergence cannot be ensured, which could suggest that the descriptor rejects much of the incorrect matches, thereby avoiding local convergence and improving efficiency. On the other hand, without the curve constraints, the registration errors are similar to the proposed method. However, the runtime needed for convergence is significantly increased. This suggests that the curve constraints efficiently limit the solution space, leading to much faster convergence. In conclusion, the above experiments demonstrate that the descriptor and curve constraints improve the accuracy, robustness, and effectiveness of the alignment.

Conclusions and Future Work
Heritage documentation, aimed at the conservation and protection of invaluable cultural heritage objects, is implemented by recording historical monuments and artifacts in digital forms. As a commonly used approach, laser scanning is able to collect highly accurate 3D data without damaging these historical artifacts. During data collection, inspection is important in order to check the integrity and quality of the scanned data. To fulfill this aim, this paper proposed a fast and robust probabilistic registration by combining shape descriptor and curve constraints for free-form objects (e.g., Buddha) and validated its performance using real-world artifact datasets. Comprehensive experiments showed that the proposed method performed well in terms of robustness and runtime and outperformed commonly used approaches. Although the method was able to provide satisfactory results, it still had difficulty for ambiguous objects with repetitive structures, which led to incorrect correspondences and transformation.
For future work, we are planning to design a more robust shape descriptor with good descriptiveness particularly suited for large-scale heritage monuments. We are also planning on improving the probabilistic framework to make it more general and scientific and extend its usefulness to other applications (e.g., train track matching, oil tank deformation monitoring, or reconstruction of cultural heritage). Color information can also be used in combination with our probabilistic method to improve its performance in the future research.