SPPD: A Novel Reassembly Method for 3D Terracotta Warrior Fragments Based on Fracture Surface Information

Abstract

As one of China′s most precious cultural relics, the excavation and protection of the Terracotta Warriors pose significant challenges to archaeologists. A fairly common situation in the excavation is that the Terracotta Warriors are mostly found in the form of fragments, and manual reassembly among numerous fragments is laborious and time-consuming. This work presents a fracture-surface-based reassembling method, which is composed of SiamesePointNet, principal component analysis (PCA), and deep closest point (DCP), and is named SPPD. Firstly, SiamesePointNet is proposed to determine whether a pair of point clouds of 3D Terracotta Warrior fragments can be reassembled. Then, a coarse-to-fine registration method based on PCA and DCP is proposed to register the two fragments into a reassembled one. The above two steps iterate until the termination condition is met. A series of experiments on real-world examples are conducted, and the results demonstrate that the proposed method performs better than the conventional reassembling methods. We hope this work can provide a valuable tool for the virtual restoration of three-dimension cultural heritage artifacts.

1. Introduction

Terracotta Warriors are the burial objects of China′s first emperor. They are known for their large numbers, vivid appearances, and diverse postures, and they are of great value in history, archaeology, art, and the tourist industry. However, owing to thousands of years of weathering erosion and earth crust movement, a large number of the Terracotta Warriors have been damaged and gathered in piles of fragments [1]. To restore the original appearance, the reassembly task is essential. Nevertheless, the traditional manual reassembly is tedious and laborious for archaeologists and may induce secondary damage to the fragments [2]. Advances in computer technology have made it possible to replace manual operations with virtual-reassembly techniques. Due to the simple structure and pure information, the point cloud is utilized to represent virtual fragments. Therefore, virtual reassembly is a complex problem in the field of the point cloud, involving feature extraction, point cloud segmentation, sampling, optimization, etc. To date, some methods have been exploited to address this problem. In general, these methods mainly consist of two stages: pairwise matching of fragments and the registration of matching pairs.

In the pairwise matching stage, a general opinion is to extract one or more geometric features of a pair of fragments and determine whether they match through features similarity. According to the geometric features, methods can be roughly classified into point-based [3,4], curve-based [5,6,7], and surface-based [8,9,10]. For example, Pan et al. [3] detected key points of point cloud and encoded their local features with the binary shape contest (BSC) descriptor. An energy function that models the similarity of critical points in the hybrid space of BSC feature and Euclidean geometry was introduced to formulate the corresponding matching task. Altantsetseg et al. [5] proposed a descriptor containing curves along the main direction of the feature point cluster. The main idea was to compare the curves by Fourier transform and compute the respective Fourier series. Son et al. [9] proposed the surface signature, a surface-based descriptor that describes the geometric characteristics based on the concave/convex information to measure the similarity of different fracture surfaces. Intuitively, from points to lines to surfaces, the extracted features contain more information and are less affected by noise. In addition, some researchers have fused multiple features to enhance the robustness of the algorithm. Li et al. [11] used boundary curves of two fracture surfaces to determine preliminary matching and concave–convex patches to refine the alignment, respectively.

As for the stage of point cloud registration, methods can be roughly classified into three categories [12]: distance-based [13,14,15,16,17], filtering-based [18,19,20], and probability-based methods [21,22,23,24]. The iterative closest point (ICP) [13] algorithm represents the first category, and many subsequent algorithms are based on it. Given a pair of point clouds, ICP takes Euclidian distance between them as the optimization objective and updates the correspondences of points and spatial pose iteratively. To solve more specific and complex problems, a series of ICP-based methods have been developed [14,15,16,17]. However, ICP-style methods still have some shortcomings, such as computational cost, convergence time, local optimization, noise sensitivity, and so on [25]. Li et al. [18] view the point set registration problem as a non-linear state-space model and adopt the cubature split covariance intersection filter (CSCIF) for this problem. The experiment demonstrates the precision and robustness of this method to noise and outliers. Jian et al. [21] represent point clouds as a Gaussian mixture model and reformulate the problem of point cloud registration to align two Gaussian mixtures. In addition to the three categories mentioned above, a coarse-to-fine registration strategy is commonly used for 3D point cloud registration [26]. Kim et al. [27] combined PCA with ICP. They used PCA-based global registration to roughly align the 3D data, which yields the initial coarse estimation of the position. The local registration technique, which uses ICP and the Levenberg–Marquardt algorithm, was applied to the coarse estimate to obtain a fine registration.

Recently, deep learning has also been widely applied to the field of 3D point cloud recognition, classification, segmentation, retrieval, detection, reconstruction, etc., and has achieved outstanding performance [28,29,30,31,32]. As a pioneer, PointNet [29] provides a deep learning method for 3D point cloud tasks. It takes unstructured point clouds as an input for the neural network, extracts the semantic information through point-level convolution layers, and acquires the global information through pooling layers, to encode the deep expression of the whole point cloud. Its excellent performance in classification and segmentation tasks shows the convictive effectiveness. It inspires both the matching and registration stages of the reassembly task. Various deep-learning-based methods have been proposed to solve the problem of point cloud registration [33,34,35,36]. The deep closest point (DCP) [35] model utilizes a point cloud embedding network, an attention-based module, and a differentiable singular value decomposition (SVD) layer to extract the rigid transformation. In the experimental part of this paper, the data is ideal. We can form a target object by performing a rigid transformation on any source object, which means that the two are only different in spatial position and pose. However, things are different in reality. Apart from external damage, the fragments are usually scanned separately, which means different densities, shapes, and surface structures of point clouds to be registered, so DCP does not perform as well as it should.

This work proposes a novel fracture-surface-based reassembly method for 3D Terracotta Warrior fragments, which consist of SiamesePointNet (proposed based on Siamese network [37], PointNet [29]), PCA [27], and DCP [37]. The flowchart of this method is shown in Figure 1.

Firstly, we build a dataset of all the fragments to be reassembled. If the number of fragments in this dataset is more than 1, there may be a matching relationship. Secondly, we enter the matching stage where we preprocess the data, extract the fracture surfaces of fragments, and judge the similarity between fracture surfaces according to the score calculated with SiamesePointNet. Thirdly, we select the most similar fracture surfaces in the registration stage and apply them to a coarse-to-fine registration method combining PCA and DCP. Then, the dataset is updated by replacing the original matching fragments with the registration result. The above process is repeated until any one of the two termination conditions is met: either the number of fragments in the dataset is equal to 1, or the maximum similarity score is not greater than the threshold ζ, which means that there is no potential matching relationship anymore. This is a greedy strategy, taking full account of the intermediate registration products, ensuring that each round of reassembly objects is the optimal choice, and ultimately guarantees global optimization and avoids mismatch.

Our key contributions can be summarized as follows:

• A novel reassembly method for 3D Terracotta Warrior fragments is proposed;

• A deep network called SiamesePointNet to evaluate the similarity of fracture surfaces is proposed;

• A coarse-to-fine registration method combining PCA and DCP is proposed;

• A series of real-world-based experiments have been conducted to prove the validity of the proposed method.

The rest of this work is structured below. In Section 2, data preprocessing and fracture surface extraction, SiamesePointNet network, and PCA-DCP based registration method are introduced. In Section 3, we conduct a series of real-world-data-based experiments to demonstrate the effectiveness of our approach. Finally, we present a conclusion and outlook in Section 4.

2. Materials and Methods

2.1. Data Preprocessing and Fracture Surfaces Extraction

As real-world objects, Terracotta Warriors should be represented in a suitable form, and fracture surfaces should be extracted for matching and registration stages. The fragment data came from the Emperor Qinshihuang′s Mausoleum Site Museum and was collected by Northwestern University with an Artec Eva handy scanner. As there is a large amount of noise in the preliminary scanned data, it is necessary to denoise it with Artec Studio Professional and Geomagic Wrap 2017 first. The fracture surface is selected as the data source because of its more extensive area, richer information, and continuous shape. This work follows Li′s approach [11] to accomplish region segmentation to strip fracture surfaces from the entire fragment. Additionally, random sampling is implemented to prepare normalized data as input to the following network. These steps are shown in Figure 2.

2.2. Network Architecture for Matching Stage

Inspired by the Siamese network [37], which is used to judge the similarity between 2D images, we propose a modified network to calculate the similarity score between 3D point clouds, named SiamesePointNet. The network architecture is depicted in Figure 3. The main intention of SiamesePointNet is to measure similarity in a quantitative way to provide a reliable basis for determining the matching relationship and selecting registration objects.

Our network takes a Siamese structure as the overall framework. The feature is extracted through two parallel PointNet networks. The weights are adjusted by matching labels so that SiamesePointNet can migrate from the original classification and segmentation task of PointNet to evaluate similarity. Intuitively, this network mainly includes two parts: a feature extractor module and a scoring module. The two inputs (point cloud A and B) are in pairs. Specifically, they are denoted as $A=\left\{x_1,\cdots ,x_i,\cdots ,x_N\right\}\subset {ℝ}^3$, $B=\left\{y_1,\cdots ,y_j,\cdots ,y_N\right\}\subset {ℝ}^3$, respectively. The feature extractor module is obliged to obtain a correct mapping $F:{ℝ}^{\mathrm{N}\times 3}\to {ℝ}^{\mathrm{K}}$ from a collection of scattered points to a K-dimensional normalized feature vector, which is easy to compare. On account of the feature transformation, point-level convolution, and pooling layer, PointNet [29] can map point cloud to a meaningful feature space, regardless of the pose and coordinate system, which disturbances the judgment of the similarity. Therefore, PointNet is embedded in the feature extractor module, where A is mapped to a feature vector $F\left(A\right)=\left(a_1,a_2,\cdots ,a_k\right)$, and B is the same. In the scoring module, an energy function $E\left(A,B\right)=\Vert F\left(A\right)-F{\left(B\right)\Vert }_2$ is used to describe the distance between two feature vectors, and the final similarity score is negatively correlated with this $E\left(A,B\right)$. Therefore, point clouds with similar shapes will be mapped to feature vectors close to each other, and finally, lead to a high score. In other words, a high score suggests that the two fragments may be sub-parts of an original fragment and break on the fracture surface, so they can be reassembled to recover the complete fragment, i.e., they match.

To constrain the network to associate the matching relationship with the similarity score correctly, the contrastive loss is used and shown as below:

(1)$loss=\frac{1}{n}{\displaystyle \sum }f{g}_i\sqrt{f{p}_i}+\left(1-f{g}_i\right)\sqrt{\{max(1-f{p}_i,0)\}}$

where $fg$ refers to the label in the ground truth (1 for the match, 0 for not), and $fp$ refers to the output similarity score. This score, roughly distributed between 0–1, is a manifestation of the possibility of matching. After the network training is completed, the output score guides selecting registration objects for the registration stage, the subsequent stage of the proposed reassembly method.

As shown in the example in Figure 4, given two fragments, the similarity scores among all fracture surfaces are calculated; the matching pair with the highest score is selected as the objects to be registered.

2.3. PCA-DCP Based Method for Registration Stage

The problem in the registration stage can be stated as:

Given $A=\left\{x_1,x_2,\cdots ,x_N\right\}\subset {ℝ}^3$, $B=\left\{y_1,y_2,\cdots ,y_N\right\}\subset {ℝ}^3$, our purpose is to find apposite rigid transformation parameters $\left[R_{xy},T_{xy}\right]$ where $R_{xy}\in \mathrm{SO}\left(3\right)$ and $T_{xy}\in {ℝ}^3$ which transform $A$ to approach $B$. Transform result $A^{*}$ is obtained by applying the parameters to $A$. If transformation parameters correct, $A^{*}$ and $B$ will completely coincide in position in an ideal situation.

As shown in Figure 5, our registration method follows the coarse-to-fine principle, where the PCA method obtains the initial position estimation, and DCP [35] brings the fine adjustment result on this basis.

Inspired by [38], the main idea of the PCA coarse registration method is utilizing PCA to obtain the main axes system (MAS, a general direction descriptor, consists of three mutually orthogonal unit vectors, each called an axis) of the fracture surface, and then registering MASs instead of the whole surfaces. Since MAS is attached to the fracture surface, when MASs are registered, fracture surfaces will also be registered by applying transformation parameters for the entire fracture surfaces. Furthermore, because the fracture surface is a subset of the fragment, they share the same coordinate system. In the same way, the registration parameters of the fracture surfaces can also be applied to the whole fragments. Therefore, the registration of the two fragments is completed by registering two MASs of fracture surfaces; Figure 6 shows this process.

However, the PCA registration method is ambiguous. Since each axis in a MAS has two corresponding potential choices among all the axes in another MAS, there are eight possible registration schemes between them. Under the right-hand rule, half of the potential situations are excluded so that only four strategies remain. For lack of effective priori methods to determine which of the four possible schemes is correct, we can implement them and choose the best one by measuring their performances. For fracture surfaces A and B, the MAS matrix is denoted as $MA=\left[\begin{array}{ccc}{a}_1& a_2& a_3\end{array}\right]$ and $MB=\left[\begin{array}{ccc}{b}_1& b_2& b_3\end{array}\right]$, and every column vector stands for the direction of an axis. The details of 4 potential schemes are described as follows and shown in Figure 7.

$MA=\left[\begin{array}{ccc}{a}_1& a_2& a_3\end{array}\right],MB=\left[\begin{array}{ccc}{b}_1& b_2& b_3\end{array}\right]$

$MA=\left[\begin{array}{ccc}-a_1& -a_2& a_3\end{array}\right],MB=\left[\begin{array}{ccc}{b}_1& b_2& b_3\end{array}\right]$

$MA=\left[\begin{array}{ccc}{a}_1& -a_2& -a_3\end{array}\right],MB=\left[\begin{array}{ccc}{b}_1& b_2& b_3\end{array}\right]$

$MA=\left[\begin{array}{ccc}-a_1& a_2& -a_3\end{array}\right],MB=\left[\begin{array}{ccc}{b}_1& b_2& b_3\end{array}\right]$

In any scheme, the rotation matrix is given by:

(2)$R=MB\ast M{A}^{-1}$

The mean of coordinates should be acquired to obtain the translation vector.

(3)$\overline{A}=\frac{1}{N}{\displaystyle \sum }_{i=1}^N{x}_i$

(4)$\overline{B}=\frac{1}{N}{\displaystyle \sum }_{j=1}^N{y}_j$

The translation vector $T=\overline{B}-\overline{A}$, and corresponding registration results can be obtained:

(5)$A^{*}=A\ast R+T$

The essence of the evaluating registration result is to evaluate the degree of overlap between $A^{*}$ and $B$. The shortest distance indicates the degree of overlap from point to surface.

(6)$D_{AB}=\frac{1}{n}{\displaystyle \sum }\Vert p^{*}{}_{y_j}-p_{x_i}\Vert ,$

(7)$p^{*}{}_{y_j}=argmin(\Vert p_{y_j}-p_{x_i}\Vert )$

(8)$D_{BA}=\frac{1}{n}{\displaystyle \sum }\Vert p^{*}{}_{x_j}-p_{y_i}\Vert ,$

(9)$p^{*}{}_{x_j}=argmin(\Vert p_{x_j}-p_{y_i}\Vert ).$

This $D$ with a small value indicates a good registration performance. All these four schemes are implemented, and the scheme corresponding to the smallest $D$ is elected as the final decision. Based on PCA coarse estimation, the result of fine registration can be further obtained through DCP [35]. Avoiding the iterative strategy, point-level optimization, and risk of non-convergence, our PCA-DCP method has an excellent performance in efficiency and robustness.

3. Results

All experiment data is digitized by the Visualization Laboratory of Northwest University, Shaanxi Province, China. SiamesePointNet network for the matching stage, the PCA-DCP-based method for the registration stage, and reassembly of multiple real-world Terracotta Warrior fragments are discussed in this section. The experiment is conducted on a computer with an Intel(R) Core [email protected] CPU and two NVIDIA GeForce RTX 2080Ti GPUs.

3.1. SiamesePointNET Network for Matching Stage

A small dataset containing 1800 matched pairs and 1560 non-matched pairs is collected and divided into training, validation, and test sets according to the ratio of 6:2:2. Throughout the experiment, we keep the feature dimension K = 1024. We follow [29] and experiment with an input point number of 1024. The accuracy achieves 96.77 and 93.69 in matched samples and non-matched samples, respectively. The average accuracy between two classes is 95.39, and the overall accuracy of all samples is 95.23. It demonstrates the effectiveness and applicability of SiamesePointNet in a small dataset. After that, the sampling rate is changed in the course of data processing, and the impact of the scale of input points is discussed on the proposed network. Intuitively, it is learned from Table 1 that the trend of accuracy gradually increases as the number of input points increases. When the number of points reaches 2048, accuracy achieves a state of saturation. According to the principle of Occam′s razor, 2048 is determined as the number of input points of our network.

For comparison, in contrast to our deep learning approach, two traditional descriptors are applied to the matching stage: FPFH [39] and SHOT [40] features. Their performance on the same dataset is recorded in Table 2. The default number of points is kept at 2048. The experimental results show that the proposed method is superior to traditional methods on the same dataset.

3.2. PCA-DCP Based Method for Registration Stage

3.2.1. Comparison of Registration Results of Different Methods

The data for the registration method is 360 matching fracture surfaces, which can be correctly registered if adequately handled. Table 3 demonstrates the effectiveness of our method and its peers. Firstly, the success rate of each method is recorded as a rough metric. Registration results that differ too much with ground truth in direction and distance are defined as failures, a contrast to success. After eliminating these failures, the MSE, RMSE, and MAE of the corresponding points of the rest are calculated to measure their performance in detail. Finally, the average number of iterations is recorded to reflect consumption.

The first part is the vertical and quick comparison. As shown in Table 3, the PCA-DCP method in this paper has the highest success rate, and its MSE, RMSE, and MAE are slightly higher than the RANSAC method. Between iterative methods, the participation of the PCA method in this paper reduces the resource overhead.

The second part is the horizontal and detailed comparison and analysis. Due to the iterative computation and dependence on the initialization, traditional ICP performs the worst in success rate, error, and consumption. The RANSAC [41] method in row 3 is based on the FPFH feature, and its error performance is outstanding. However, the error performance is mainly due to the low success rate. We cannot deny that it performs very well on several successful samples, but such a success rate is unacceptable. The success rate of NDT [42] is too low to cope with real application scenarios. The results of methods using only PCA or DCP are in rows 5 and 6, respectively. The PCA method achieves a high success rate but is not excellent enough in terms of error, while the DCP method performance garners a poor success rate, but the error is smaller. In the 8th row, our PCA-DCP combination method achieves the highest success rate and slight error close to RANSAC, outperforming both PCA and DCP, which illustrates the mutual positive effects of the two ways.

Since the PCA method provides a reasonable initial estimate, it can also serve as a preprocess for other methods, such as ICP. Compared with traditional ICP, PCA-ICP has improved a lot. However, it is still inferior to PCA. Through the experimental results, we believe that PCA can provide good initialization for many cases where ICP would fail, which is reflected in the improvement of the success rate compared to ICP. Nevertheless, as a global coarse registration method, the initial position estimation provided by PCA is not accurate enough, so it is easy to make the PCA + ICP method fall into a local optimum in the process of iteratively calculating the corresponding points, which eventually leads to a decrease in performance compared to PCA. The comparison in rows 7 and 8 shows that PCA and DCP are much better than the combination of PCA and ICP.

PCA is outstanding in initial estimation but not accurate enough in detail, while DCP performs well in detail but depends on initialization. Combining them can absorb their advantages and remedy their shortcomings, showing the effect of “1 + 1>2”.

To examine the performance of these methods on each sample, we save all the successful registration results and sort them according to MAE, as shown in Figure 8. It can be concluded that our method maintains the best performance in terms of robustness and details on different samples compared to other methods.

Figure 9 is a comparison of different methods on a specific pair of real-world data. In this example, PCA obtains a coarse initial pose estimation but is not good enough in detail. When no initial estimate is given, the DCP method leads to failure: there is a significant deviation between the two fracture surfaces in spatial position. The PCA-DCP method proposed in this work achieves the best result. Figure 10 shows the registration visualization of our method on several samples.

3.2.2. Robustness Experiments

The robustness of our registration method is verified concerning the initial angle, initial distance, and noise, and the average MAE over all samples is shown in Figure 11. In the angle-robustness experiment, different upper limits for the initial angle (15°, 30°, 45°, 60°, 75°, 90°, 105°, and 120°) are set, and a specific fracture surface from ground truth is randomly rotated within the upper limit. The distance-robustness experiment is similar, and the maximum initial distance (10 mm, 20 mm, 30 mm, 40 mm, 50 mm, 60 mm, 70 mm, 80 mm, 90 mm, and 100 mm) specifies the translation distance between two surfaces. In the noise-robustness experiment, Gaussian noise is added to ground truth with gradually attenuated signal-to-noise ratios (SNR). The initial value of SNR is set to 100 db. As the attenuation rate increases from 0 to 70%, samples in ground truth gradually become noisier and noisier. The average MAE of the RANSAC method over all samples is significant, reaching more than 30, although its performance is stable on datasets with different angles, distances, and noise. Similarly, unlike distance-based methods, the results of NDT often show serious deviations in the registration direction, resulting in a higher mean MAE. Therefore, their MAE are not shown in the figure.

It is worth noting that all non-convergent samples are filtered. In ICP-participated methods (traditional ICP and PCA + ICP), there is no statistical result for MAE if it does not converge. A threshold with a value of 200 is set, and once the number of iterations exceeds this threshold, it is considered non-convergence.

Due to the randomness in data manufacturing, we should pay more attention to the overall trend rather than the value of each sample. Our method exhibits the most stable and robust properties no matter the initial angle, distance, or noise. Figure 12 shows the registration results of several noisy fracture surfaces.

3.3. Reassembling Feasibility on Multiple Real-World Terracotta Warrior Fragments

To investigate the effectiveness of our framework, several sets of experiments are conducted on multiple real-world Terracotta Warrior fragments. Due to the diversity and uncertainty of the real-world Terracotta Warrior fragments reassembling scene, the performances of the proposed framework are different in specific cases. This section is arranged to demonstrate the feasibility of this framework through a complete reassembly process of a typical sample. A further discussion of the limitations of the proposed method will be mentioned in Section 4.

The results are satisfactory when fracture surfaces are sufficiently similar, as shown in Figure 13. Figure 13a shows the five pieces we are faced with. In this scenario, the largest fragments reach a volume of about 2000 cm², and the smallest fragments are less than a tenth of that. However, since the data concerned by the proposed method is the fracture surface, the area difference becomes small after the fracture surfaces are extracted. In the beginning, a dataset for five fragments is created. Then, we calculate the similarity among fracture surfaces, select the most similar pair “3–4” to register and update the dataset, as shown in Figure 13b. At this stage, only four fragments are left. The original fragments “3” and “4” are replaced by a merged fragment “3–4”, and the fracture surface where the registration was implemented was blocked and will not be considered in the next round of matching and registration steps. As with the above process, after several steps, we obtain the result shown in Figure 13e. At this point, there is only one fragment left, which suggests that all the fragments have been reassembled together, and our termination condition is met so that the reassembly process ends.

4. Discussion

This article proposes a 3D fragments reassembling method named SPPD based on fracture surface similarity, which mainly includes a similarity measurement network called SiamesePointNet and a PCA-DCP registration method. PointNet is embedded into a Siamese structure to quantitatively analyze the similarity of fracture surfaces in the matching stage. PCA is used to obtain a rough initial position estimate in the registration stage, and DCP is used to refine the results further. Extensive experiments prove the accuracy of SiamesePointNet on small-scale samples, the efficiency and robustness of the PCA-DCP method, and the effectiveness of the proposed method on real-world Terracotta Warriors.

Given the richness and reliability of the information on the fracture surface, the core of our method lies in the fracture surface. Once the information on the fracture surface is unreliable, for example, if the fracture surfaces are severely damaged or partially similar, problems will arise in our method. Figure 14 reflects the situation of severe damage, and our SPPD method judges that the two fragments do not match. This is a wrong judgment contrary to the ground truth.

Although our registration method can obtain excellent registration results, there is a premise that the similarity of the fracture surfaces is sufficiently high. As shown in Figure 15, when given fragments with low similarity, our registration result will be unsatisfactory in detail. Therefore, how to complete the reassembling task with insufficient information on the fracture surface is a problem to be overcome in our later work.

In the reassembly process of entire Terracotta warriors, severe damage to the fracture surfaces inevitably occurs. At present, the proposed method undertakes the reassembly of parts with high integrity. After the reassembly of these parts is completed, the rest are difficult parts requiring manual guidance from experts. Figure 16 shows the reassembling result of an entire Terracotta warrior. The number of fragments to be processed is 49, and seven intermediate results were successfully obtained, leaving nine individual fragments without matching objects. These are then artificially reassembled together to obtain the final result. It is believed that the problem of incomplete fracture surface will be overcome in the future, and the dependence on manual operation will be gradually reduced.

In this framework, the computational overhead caused by the rising number of input fragments cannot be ignored. Fortunately, however, this number will not keep increasing endlessly. In the excavation and preservation of the fragments, there is a weak spatial continuity. That is to say, the fragments from the same Terracotta Warrior are close to each other in space. Therefore, the maximum number of fragments that our system needs to deal with at one time is usually not more than 100, which is an acceptable range. Even so, the search for a faster solution is a potential future research direction.

5. Conclusions

In conclusion, the purpose of the current study is to reassemble the 3D fragments correctly and effectively. The experiment demonstrates the practicability of the real-world Terracotta Warriors of the proposed method. We hope this work can provide a valuable tool for the virtual restoration of three-dimension cultural heritage artifacts.

Author Contributions

Conceptualization, Wenmin Yao and Xin Cao; methodology, Xin Cao and Kang Li; software, Wenmin Yao; validation, Xin Cao and Kang Li; investigation, Wenmin Yao and Wenlong Tang; resources, Wenmin Yao; data curation, Tong Chu; writing—original draft preparation, Wenmin Yao; writing—review and editing, Xin Cao; visualization, Fengjun Zhao and Jingyu Wang; supervision, Guohua Geng and Mingquan Zhou. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Key Research and Development Program of China (No. 2019YFC1521102); the Project funded by China Post-doctoral Science Foundation, grant number No. 2018M643719, the Young Talent Support Program of the Shaanxi Association for Science and Technology, grant numberNo.20190107; The Key Research and Development Program of Shaanxi Province (No. 2019GY-215), The Major research and development project of Qinghai(No. 2020-SF-143).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available on request from the corresponding author. The data are not publicly available due to that the Terracotta Warriors involve a policy of secrecy over cultural heritage.

Acknowledgments

We thank the Emperor Qinshihuang′s Mausoleum Site Museum for providing the Terracotta Warriors data.

Conflicts of Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. The flowchart of the proposed method.

Figure 2. Flow chart of data processing. (a) Fragment denoising; (b) fracture surface extraction; (c) random sampling result.

Figure 3. SiamesePointNet network architecture.

Figure 4. An example of calculating the similarity score of the fracture surface and select objects to be registered. (a) Two fragments and their fracture surfaces; (b) all potential matching pairs; (c) “2–2” match with the highest similarity score; (d) registration of pair “2–2”.

Figure 5. Coarse-to-fine registration structure.

Figure 6. Schematic diagram illustrating the PCA method. (a) MASs of two fracture surfaces; (b) registration of two MASs; (c) fracture surfaces registration; (d) fragments registration.

Figure 7. Potential registration schemes between 2 fracture surfaces. F1 and F2 represent different fracture surfaces, and their direction of the text reflects the order of registration. (a) The potential registration scheme where MA=[a1a2a3],MB=[b1b2b3]; (b) The potential registration scheme where MA=[−a1−a2a3],MB=[b1b2b3]; (c) The potential registration scheme where MA=[a1−a2−a3],MB=[b1b2b3]; (d) The potential registration scheme where MA=[−a1a2−a3],MB=[b1b2b3].

Figure 8. Sorted MAE comparison.

Figure 9. Comparison of different methods on a specific pair of real-world data.

Figure 10. Registration visualization on several samples.

Figure 11. Experiment results on robustness. (a) Mean MAE on different initial angles; (b) mean MAE on the different initial distances; (c) mean MAE on noisy data.

Figure 12. Examples of registering several noisy point clouds with our method.

Figure 13. The reassembly result of multiple fragments from a Terracotta arm. (a) Five fragments; (b), (c), and (d) steps to reassembly; (e) reassembly result.

Figure 14. Judgment in severely damaged condition. (a) Two fragments; (b) judgment obtained by our SPPD method that the two do not match; (c) ground truth.

Figure 15. Reassembly results in partially similar conditions. (a) Two fragments; (b) the reassembly result of our SPPD method; (c) ground truth.

Figure 16. Reassembling result of an entire Terracotta warrior (the head and feet can no longer be found; they may have been devastated).

Table 1

Accuracy on different numbers of input points.

Points Number	256	512	1024	2048	4096
Acc avg class	92.85%	93.10%	95.23%	96.53%	95.60%
Acc overall	92.86%	93.15%	95.39%	96.58%	95.83%

Table 2

Accuracy of different methods.

Method	FPFH	SHOT	OURS
Acc avg class	77.93%	77.04%	95.60%
Acc overall	77.38%	77.08%	95.83%

Table 3

Registration results of different methods.

Method	Suc_Rate	MSE	RMSE	MAE	Iteration
ICP	20.56%	9.48	3.08	1.73	23.74
RANSAC	23.33%	1.17	1.08	0.64	------
NDT	4.72%	13.83	3.72	1.95	------
PCA	51.67%	5.28	2.30	1.18	------
DCP	21.38%	4.26	2.06	1.01	------
PCA + ICP	35.83%	4.77	2.19	1.21	9.03
PCA + DCP	53.89%	3.11	1.76	0.92	------