Pairwise Matching for 3D Fragment Reassembly Based on Boundary Curves and Concave-Convex Patches

We present a simple pairwise matching method for 3D fragment reassembly that uses boundary curves and concave-convex patches to accelerate and optimize the matching. Given the boundary curves of two fracture surfaces that completely or partially coincide, we can quickly exclude the fracture surface pairs with low boundary curve similarities, which can substantially reduce the computational cost of the subsequent patch matching, where we extract and delineate concave-convex patches of the selected fracture surfaces. A modified iterative closest point algorithm is applied on these concave-convex patches to refine the alignment. Finally, we determine the matched 3D fragments according to the overlap ratio of their fracture surfaces. The results of experiments on real-world examples demonstrate that our proposed algorithm is both accurate and efficient.


I. INTRODUCTION
Computer-aided 3D fragment reassembly is highly important in many fields, such as archaeology, paleontology and medicine. In many cases, 3D reconstruction focuses mainly on the geometry of the fragments, since information such as color and texture is not reliable for assembly [1]. Hence, shape matching and 3D scanned digital fragment alignment are challenging problems in computer graphics and computer vision [1], [2].
Geometry-matching-based reconstruction consists of two main successive steps: pairwise matching and multi-fragment matching. The former focuses on pairwise matching of the fracture surfaces among all fragments, while the latter aims at reconstructing the entire original object according to a set of candidate matching pairs of fragments. Pairwise matching is crucial to multi-fragment matching and has been a very important technique for 3D fragment alignment [3], [4].
Recently, studies [3]- [13] matched fragments based on features of the fracture surface, which used feature patches [4], [5], feature curves [3], [6], [8], and feature points [9]. For example, Huang et al. [4] clustered points to generate The associate editor coordinating the review of this manuscript and approving it for publication was Zahid Akhtar . feature patches of a fracture surface and used a forward search algorithm to identify similar feature patches for pairwise matching. Wu and Wang [6] used the Hausdorff distance and the modified 4-point congruent set algorithms to identify potential simple chordless cycles for matching fragments for the reassembly of fractured sand particles. In their follow-up work [7], for matching individual fragments to their mother particle, the iterative closest point (ICP) algorithm was used to attempt to match each face of each child particle to the mother particle and to evaluate the matching results. Altantsetseg et al. [8] guided the pairwise matching of fragments by applying a descriptor that contained a cluster of feature points to the fracture surfaces and curves along the principal directions of the cluster. Others proposed non-feature-based methods for the pairwise matching of fragments. Winkelbach and Wahl [10] used a binary tree structure in a depth-first approach. Mavridis et al. [11] found an optimal match according to all vertices on the fracture surfaces and used the centroid and the normal distribution of the fracture surfaces or a standard RANSAC (Random sample consensus) for coarse alignment. Stanford's Digital Forma Urbis Romae project [12] dealt with heavily eroded fragments. Instead of using the geometry of the fracture surface, reconstruction was conducted by matching annotated incisions on the fragments' top surfaces. In medicine, Paulano et al. [13] presented an automatic method for calculating the contact zone between two bone fragments from computed tomography (CT) images and applied the ICP algorithm to the fragment to reduce the fracture.
In this work, we combine curve matching and surface matching to improve the accuracy and efficiency of fragment matching. We directly extract the 3D boundary curve and surface features from the fracture surfaces, which are used for coarse matching and fine alignment, respectively. Our method can handle arbitrarily shaped fracture surfaces and partial matches. Our work is inspired by the work of Huang et al. [4]; however, there are several differences between our method and theirs in the aspects of fragment preprocessing, surface segmentation, feature extraction and the process of pairwise matching.
Our algorithm makes two main contributions: First, fast and simple curve matching can exclude most pairs with low surface similarities while keeping only a few candidates with high similarities. This can massively reduce the computational cost of the subsequent patch matching. Second, after extracting concave-convex patches of fracture surface pairs, a modified iterative closest point algorithm is proposed for application on the patches to optimize the assembly results.
In the experiments, we evaluate the performance of the proposed algorithm on real-world examples.
In the following, we introduce two geometric features, namely, boundary curves and concave-convex patches, in Section II and Section III, respectively. In Section IV, we introduce the proposed pairwise matching approach by using the extracted features. The experimental results are presented in Section V. Finally, we present the conclusions of the paper in Section VI.

II. BOUNDARY CURVE EXTRACTION AND COMPARISON A. IDENTIFYING THE FRACTURE SURFACE AND EXTRACTING THE BOUNDARY CURVE
Our target object is typically broken into several fragments, whose fracture surfaces are often bumpy; these bumps are defined as noise. Therefore, we first smooth the fragments by using adaptive smoothing [14] without shrinking, which has the advantage of controlling the degree of noise by adjusting a smoothing radius. Fig. 1 shows an original brick fragment and its smoothed version. The degree of noise of the smoothed fragment is obviously reduced, which facilitates the extraction of more accurate concave and convex patches. After smoothing, to segment the fragments into a set of surfaces that are bounded by sharp edges, we distinguish edge points and non-edge points using multi-scale curvedness [15]. The curvedness at a point p on a surface can be estimated as: where k 1 and k 2 are the maximum and minimum principal curvatures, respectively. The curvedness can be used to indicate how sharp or gentle a surface is. We use the multiscale method that is proposed in [15] to robustly estimate the curvedness by fitting a quadratic surface to neighborhood points at various scales [16]. Then, surface segmentation is realized by using a simple patch-growing algorithm [3]. Afterwards, we use the surface roughness [4] to recognize coarse fracture surfaces created by the broken object, which typically differ from the smooth and fine surfaces of the unbroken object. If the original and fracture surfaces do not differ substantially in terms of roughness, our algorithm cannot distinguish them. In this case, we must match all surfaces, including all original and fracture surfaces. However, the original surfaces can be easily excluded in the stage of pairwise matching because it is difficult to find proper surfaces that match with them.
We now have identified all fracture surfaces of the fragment; however, their boundary points are unordered and inaccurate. Thus, we use the boundary tracking algorithm [17] to obtain ordered and more accurate boundary points, which constitute the final boundary curves. Fig. 2 illustrates the extraction procedure of the fracture surface and the boundary curve of a fragment of a stone brick. Fig. 2(a) presents the multi-scale curvedness of the key points and Fig. 2(b) presents the segmentation of the fracture surface. After boundary tracking, we obtain the boundary curves, as shown in Fig. 2(c).

B. BOUNDARY CURVE COMPARISON
When two fracture surfaces can be matched, their boundary curves should completely or partially coincide. Therefore, we can compare two boundary curves to judge their similarity. As illustrated in Fig. 2(c), because the extracted boundary curve includes a few noise points, the boundary curve is smoothed using a Gaussian filter to remove the noise prior to the similarity comparison. The points p(i) on the extracted boundary curve and the Gaussian filter f G for p(i) are defined as follows: A filter of width w = 6 performs well, and the performance at a specified filter width is related to the resolution of the boundary curve. As points with a distance of greater than 3σ from the central point make a negligible contribution when applying a Gaussian filter with a standard deviation of σ , the filter is designed using σ = w/3. The details are available in [3]. Fig. 2(d) shows the smoothed boundary curve after the application of the Gaussian filter.
We compare the boundary curves C 1 and C 2 , which are described by strings . . , (κ 2 n , τ 2 n )}, based on the discrete curvature κ and torsion τ . Here, m and n are the numbers of points of C 1 and C 2 , respectively. This problem is addressed as a circular substring matching problem. This substring matching problem is based on an m × n similarity matrix , in which element (i, j) is the difference between points p 1 i ∈ C 1 and p 2 j ∈ C 2 , which is expressed as the mean Euclidean distance of p 1 i and p 2 j in the vicinity of the two points [3]. It is defined as follows: We assume that the two points are similar if (i, j) < ε . A similar curve segment is a sequence of consecutive diagonal elements of : The boundary curves of fragments that are matched partially coincide only on a small segment. We require that two similar surfaces share a boundary of at least 1/8th the length of the shortest boundary curve to avoid failing to detect the correct match.

III. CONCAVE-CONVEX PATCH EXTRACTION AND COMPARISON A. CONCAVE-CONVEX PATCH EXTRACTION
Based on the mean curvature and the Gaussian curvature at each point on a surface, local surfaces that surround the point can be subdivided into eight basic types: peak, pit, ridge, saddle ridge, valley, saddle valley, flat and minimal surface [18]. After analyzing various material fragments, we find that concave patches and convex patches mainly include pit points and peak points, respectively. The main point types between concave patches and convex patches are saddle ridges and saddle valleys.
We cluster peak points to form convex patches and cluster pit points to form concave patches. There may be several non-peak points within a group of peak points or several non-pit points within a group of pit points; we merge them to large adjacent groups. Because fracture surfaces are complicated and irregular, an ideal surface type is rarely expected. We set two thresholds, namely, ε K and ε H : if the Gaussian curvature satisfies K ≥ ε K and the mean curvature satisfies H ≤ ε H , the point is a peak point; if K ≥ ε K and H > ε H , the point is a pit point.

B. CONCAVE-CONVEX PATCH COMPARISON
We define each patch R as follows: where µ(R) and σ (R) are the mean and the standard deviation, respectively, of the curvature values of all points on patch R, which are used to represent the concavity or convexity degree of the patch R. The size signature S(R) and anisotropy signature A(R) that are associated with R are computed via principal component analysis (PCA) [19] as follows.
Let m be the total number of points on the patch R and H i be the mean curvature of a point p i ∈ R. Then, µ(R) and σ (R) are respectively defined as Let p be the barycenter of patch R = {p 1 , . . . , p m }, and let M be a 3 × 3 covariance matrix: where λ R 0 ≥ λ R 1 ≥ λ R 2 and n R 0 ,n R 1 ,n R 2 are the eigenvalues and the corresponding eigenvectors, respectively, of M . λ R i and n R i are also known as the principal components and the principal directions of R. S(R) = (λ R 0 + λ R 1 + λ R 2 ) 1/2 and A(R) = λ R 1 /λ R 2 1/2 are the size signature and the anisotropy signature, respectively. The size similarity SS(R 1 , R 2 ) and the VOLUME 8, 2020 anisotropy similarity AS(R 1 , R 2 ) are defined as follows: If two patches R 1 and R 2 satisfy the following: we claim that (R 1 , R 2 ) is a similar patch pair.

IV. PROPOSED PAIRWISE FRAGMENT MATCHING
If two fracture surfaces match with each other, their boundary curves will completely or partially coincide, and the convex/concave patches of one fracture surface should match with the concave/convex patches of the paired fracture surface. In addition, since curve matching is much faster than surface matching, it can directly exclude many redundant surface pairs and reduce the high computational complexity of the subsequent concave-convex patch matching. Our pairwise matching method consists of three steps: (a) Two fragments are roughly aligned quickly via boundary curve matching. (b) A modified ICP algorithm that is based on concave-convex patches is used for fine alignment, which is more effective than the ICP algorithm that is based on all points. (c) A validity testing strategy that is based on the overlap ratio of two fracture surfaces is applied to determine whether the two corresponding fragments match well. In our algorithm, when two fragments match successfully, both their fracture surfaces and their boundary curves match. Therefore, our algorithm yields fewer mismatches than methods that only use curve matching or surface matching.

A. BOUNDARY-CURVE-BASED COARSE ALIGNMENT
For boundary curves of two fracture surfaces, we employ the boundary curve comparison algorithm that was proposed in Section II and obtain all similar curve segment pairs. For each pair of similar curve segments, we calculate a transformation via the dual quaternion method [20] to coarsely align the two fragments. Next, we assign scores to all coarse alignments according to the lengths of the overlapping curve segments and store the corresponding coarse alignments in a priority queue. Starting from the head of the queue, we apply a modified ICP algorithm that integrates concave-convex patches to obtain a fine alignment.

B. CONCAVE-CONVEX-PATCH-BASED FINE ALIGNMENT
For fracture surface matching, in most cases, part of one surface coincides with part of the other surface. We present a modified ICP algorithm that is based on the probability iterative closest point algorithm (PICP) [21] and apply it to the concave-convex patches to obtain the fine alignment.
First, we use points on concave-convex patches instead of all surface points to find the closest points. For each point x i ∈ R 1 j ⊂ S 1 , we find its closest point y i ∈ R 2 j ⊂ S 2 , where patch R 2 j is the similar patch of R 1 j . If R 1 j is not similar to patch R 2 j on S 2 , we discard point x i . The similarity between patch R 1 j and R 2 j is computed according to formula (10). In this case, S 1 and S 2 either partially coincide or do not coincide at all. Moreover, we only use some of the points on the concaveconvex patches instead of all points on the fracture surfaces for the ICP algorithm. This is sufficient for the convergence of the algorithm to the correct solution. Hence, this strategy reduces the search space for reasonable closest points and accelerates the algorithm.
Second, we use three constraints to discard outliers. The closest point pair (x i , y i ) is identified as an outlier if one of the following constraints is satisfied: (a) the distance constraint: x i − y i > ε d , namely, the distance of points x i and y i is larger than a threshold ε d ; (b) the normal constraint: n(x i ) − n(y i ) > ε n , namely, the difference between the normals at points x i and y i exceeds a threshold ε n ; and (c) the concavity-convexity degree constraint: where H (·) is the mean curvature, namely, the concavityconvexity degree differs substantially between points x i and y i .
The method of removing outliers enables us to realize part-part surface matching. Using the three constraints to remove outliers, we can obtain the closest point pairs with fewer outliers. The closest point pairs are used to compute the transformation, which yields an accurate estimate and accelerates the convergence.

C. VALIDITY TESTING STRATEGY
The validity testing strategy is based on the overlapping area ratio of two fracture surfaces. The area of a surface is described by the number of points on it. According to our experimental analysis, when the fracture surfaces matched correctly, the minimum overlapping percentage between the surfaces was approximately 20%. Hence, 20% is set as the threshold for successful matching.
Therefore, we consider two fracture surfaces S 1 and S 2 to be correctly matched or to fit well if the size of their overlapping area is larger than 20% of the area of one of the fracture surfaces [19].
To define the overlapping area between surface S 1 and S 2 , we consider points from both S 1 and S 2 and their closest points on another surface. Considering a point x i from S 1 as an example, we define y i in the overlapping area between S 1 and S 2 if it satisfies: where y i is the closest point of x i on S 1 , n(·) denotes the normal and H (·) is the mean curvature. In the implementation, we find that even if a pair of incorrectly matched fracture surfaces have several consistent point correspondences, their overlapping fracture surface is small. Therefore, most incorrect fracture surface matches are eliminated by the validity test.

V. EXPERIMENTAL RESULTS
We applied the proposed algorithm to several real-world broken objects of various materials, which included brick (stone) and cake (mortar) models from Vienna University of Technology [4]. The brick and cake fragments are presented in [22]. The triangular meshes of the head (terracotta) fragments are obtained using a 3D laser scanner, namely, Handy3dscan. All fragments are scanned from real fractured models. The experiment was conducted on a computer with a Pentium(R)/3.4 GHz CPU and 8.0 GB of RAM.

A. THRESHOLD PARAMETER SETTING
Parts of our algorithm are guided by thresholds. Following references [4] and [19], we use statistical methods to estimate the thresholds. According to the attributes of the thresholds, we classify the thresholds into two categories.
The thresholds in the first category are used in the similarity comparison of surfaces. They include mainly ε for boundary curve similarity; ε µ , ε σ , ε a , ε s in Eq. (10) for concave-convex patch similarity; and ε dis , ε θ , ε curv in Eq. (11) for validity testing. To prevent failure to detect the correct match, these thresholds are set to large values. For false matches, the corresponding error is much larger than these thresholds. We set ε = 0.3, ε µ = 0.3, ε σ = 0.2, ε a = 0.1, ε s = 0.1,ε dis = 0.01,ε θ = 0.01 and ε curv = 0.03.  The thresholds in the second category include ε K and ε H for concave-convex patch extraction and ε d , ε n and ε c for discarding outliers when finding the closest point pairs. We use all fragments to estimate and analyze these thresholds. When ε K and ε H are set in the interval [0.001, 0.01], the difference in the results of concave and convex region extraction is small. Thus, we set ε K = 0.005 and ε H = 0.005. We found that large values of ε d , ε n and ε c , which are used to prune the false point pairs, do not greatly influence the final results but do affect the running time of the ICP algorithm. Therefore, in our experiments, ε d , ε n and ε c are respectively set as where N is the total number of closest point pairs of the two fracture surfaces in the iteration.   our ICP algorithm. Fig. 3(a) shows two similar fractures and their boundary curves and concave-convex patches. The coarse alignment results of two fragments, boundary curves and concave-convex patches via similar curve segment pair (C 1 , C 2 ) in Fig. 3(a) are present in Fig. 3(b), and the fine alignment results that are obtained using our ICP algorithm based on similar concave-convex patches are presented in Fig. 3(c). Fig. 4 presents the convergence processes of PICP [21] and our modified ICP on two fragments that are shown in Fig. 3. The root mean square (RMS) error [21] is used to measure the residual error in registration algorithms. Our modified ICP method converges faster than PICP under the same convergence condition. Fig. 5 shows six brick fragments (Fig. 5(a)), their pairwise matching results (Fig. 5(b)), and the global reassembly result (Fig. 5(c)) that is based on pairwise matching. Similarly, Fig. 6 shows a matching example of a head model with five fragments.

C. PAIRWISE MATCHING AND GLOBAL ASSEMBLY VIA THE PROPOSED METHOD FOR REAL FRAGMENTS
In this work, we focus mainly on the pairwise matching of fragments. Based on our pairwise matching, we use the simple and greedy method in [5] to globally reassemble fragments. It is used to evaluate the performance of our method. Based on the proposed pairwise matching method, we can obtain all fracture surface pairs of all fragments for the subsequent complete reconstruction. The global assembly process is as follows: First, the fracture surface pair with the maximum overlap region is selected and assembled after calculating the related transformation matrix. This process is repeated for the remaining fracture surface pairs of unassembled fragments until all fragment pairs have been assembled.
For the brick and head fragments, we can find all potential matches. A complete match is a match in which two fracture surfaces coincide entirely. The remaining matches are partial matches, in which only part of the surface coincides with all or part of the other surface. The related experimental results in Fig. 5(b) and Fig. 6(c) demonstrate that our method can realize the complete or partial matching of broken objects.
Tables 1-2 compare the performances of the PICP algorithm [21] and our algorithm on the brick and head fragments that are presented above in terms of the RMS error, number of iterations and speed. The proposed method outperforms the PICP algorithm [21] with fewer iterations and higher execution speed under the same convergence condition. The accuracy of our algorithm is also slightly higher than that of the PICP algorithm [21]. Fig. 7 and Table 3 present an experimental comparison between the original surfaces and the smoothed surfaces of the brick fragments. From the original fragment in Fig. 7(a), 131 concave-convex patches ( Fig. 7(b)) were extracted. Most of these patches are too small to accurately reflect the concavity and convexity characteristics of the fracture surface. Then, 34 concave-convex patches (Fig. 7(d)) were extracted from the smoothed fragment in Fig. 7(c).

D. ABLATION STUDY OF THE PROPOSED METHOD
Compared with the original fragments, the speeds of the algorithms on the smoothed fragments increased by approximately 18.5%. For the smoothed fragments, we can find all potential matches, whereas original brick fragments 2-5 were not matched successfully. In Table 3, for those successfully matched fragment pairs, the average RMS errors of smooth fragments is approximately 0.078 of the error of the original fragments with the same number of iterations.   Although fewer concave-convex patches were extracted based on the smooth fragments, these patches can more robustly describe the concave and convex shapes and characteristics of the fracture surface. Furthermore, because the search space can be reduced and the matching speed can be increased, higher accuracy is realized.
In both archeological and medical applications, it is very common for some fragments to be unavailable. If a piece of a fragment is missing, our method can still find all potential fragment pairs in the remaining fragments. However, the global reassembly results are affected. For example, Fig. 8(a) presents the reassembly results with brick 5 missing. In this case, the absence of brick 5 will not affect the assembly result of the remaining fragments. Fig. 8(b) presents the reassembly results with brick 1 missing, in which brick 3 is isolated because it only pairs with brick 1.

VI. CONCLUSION
In this paper, we present a simple but effective method for the pairwise matching of 3D fragments, which consists of boundary-curve-based coarse matching and concave-convexpatch-based fine alignment. First, we compare the fracture surfaces' boundary curves of 3D fragments that have high similarities, which are regarded as matched surfaces. This simple boundary-curve-based comparison can facilitate fast coarse matching to quickly exclude low-similarity boundary curves from surface pairing. Next, based on the results of the previous step, we conduct fine alignment according to concave-convex patches of the selected fracture surfaces. Afterwards, we determine whether two fragments fit well based on their overlap ratio. The experimental results demonstrate that our method realizes efficient and accurate matching for various types of fragments. MINGQUAN ZHOU is currently a Professor with the School of Information Science and Technology, Northwest University, and the College of Artificial Intelligence, Beijing Normal University, a Ph.D. Supervisor, and an Outstanding Member of the Chinese Computer Society. His main research interests include virtual reality, image processing, visualization technology, and so on. VOLUME 8, 2020