Keywords

1 Introduction

The template plays a pivotal role in template-based fragments reassembly tasks. The premise that the template can work lies in that it has the same geometric scale as the fragments. However, the fragments and the template often have scale variances in practical applications. The usual solution to this problem is to conduct scale normalization for the fragments and the template. Existing scale normalization algorithms could be roughly divided into two categories: global information-based methods [1, 2] and feature point-based methods [3, 4]. However, the former is only applicable for full shapes and the latter is strongly dependent on feature points. It is often difficult to obtain correct normalization results when the fragments are small or the features of the shape are not discriminative.

This paper presents a scale normalization algorithm for template-based fragments reassembly, which effectively solves the reassembly problem of fragments and template with scale variances. The Multi-Radii Geodesic Disk Spectrum (MR-GDS) descriptor is firstly proposed which explicitly fuses the geometric scale. The scale variances are computed by comparing the GDS of the fragment and the MR-GDS of the template. The proposed algorithm needs only one pair of matched points and is suitable for fragments with few extremities. Pipeline of our algorithm is illustrated in Fig. 1, which involves the anchor point selection stage, the scale variance estimation stage, and the post verification stage.

Fig. 1.
figure 1

Pipeline of the proposed algorithm. The anchor point selection stage is on the left, the scale variance estimation stage based on MR-GDS is on the middle, and the post verification stage is on the right.

In summary, main contributions of our work is as follows:

  1. (1)

    A scale normalization algorithm is proposed, which solves the reassembly problem of fragments and template with scale variances. The Multi-Radii Geodesic Disk Spectrum (MR-GDS) is firstly proposed to conduct scale normalization.

  2. (2)

    A post-verification strategy based on Interactive Closest Points (ICP), which conducts ICP for geodesic disks with different scale normalization, is proposed to optimize the scale estimation results.

The remainder of the paper is arranged as follows. The related work will be introduced in Sect. 2. In Sect. 3, pipelines of our algorithm are discussed in detail. Next, experimental results and discussions are presented in Sect. 4. Finally, conclusions are summarized in Sect. 5.

2 Related Work

In this section, methods for 3D shape scale normalization are generally summarized. Some methods transform the problem from Euclidean space into another domain and find the scale integrally. Savva et al. [2] address the problem of recovering reliable sizes for a collection of models defined using scales with unknown correspondence to physical units. Their approach provides quite correct size estimates for 3D models by combining category-based size priors and size observations from 3D scenes.

Zaharescu et al. [3] proposed a scale-independent feature point extractor that called MeshDOG and a scale-independent feature descriptor that called MeshHOG, which extract and match feature points on rigid or non-rigid shapes. Bronstein et al. [4] proposed the Scale Invariant Heat Kernel Signature (SI-HKS) descriptor and applied it to the identification of non-rigid shapes, which solves the feature extracting and feature matching problem of different scales.

Some methods find corresponding feature points firstly and then conduct matching. Sahillio\(\breve{g}\)lu and Yemez [5] presents a scale normalization method for isometric shape matching in a combinatorial matching framework, which resolves the scale ambiguity by finding a coarse matching between extremities of shapes based on a novel scale-invariant isometric distortion measure. However, their method needs at least three pairs of points with accurate correspondences. In many cases, it is difficult to find enough corresponding points, especially when the area of the fragment is less than half of the template.

3 Pipeline of Our Algorithm

In this section, pipeline of our algorithm is presented, which contains the anchor points selection, the scale variance estimation based on MR-GDS, and the post verification.

3.1 Anchor Points Selection

A pair of matched points on the fragment and the template is regarded as anchor points. The anchor points are significant since the scale normalization and fragments reassembly are all based on correct anchor points. In order to ensure the accuracy of the corresponding anchor points, we proposed a strategy based on Gaussian curvature to select the anchor points interactively.

Let the fragment and the template be F and T, respectively. The non-boundary points with the top k largest Gaussian curvature on F are selected firstly, and an interval \([{{g}_{\min }},{{g}_{\max }}]\) is obtained according to the sorting k Gaussian curvature. Secondly, the points whose Gaussian curvature are within \([{{g}_{\min }},{{g}_{\max }}]\) are selected on T. Thirdly, a pair of corresponding points is selected interactively as the anchor points. In the experiments, we generally set \(k\le 10\) to reduce the difficulties of interactive selection. Figure 2 illustrates a pair of anchor points (red dots) selected on the fragment and the template.

Fig. 2.
figure 2

A pair of anchor points (red dots) on the template (a) and the fragment (b). (Color figure online)

3.2 Scale Variance Estimation Based on MR-GDS

In this section, we firstly recap the background of Geodesic Disk Spectrum (GDS) and then propose the MR-GDS descriptor. The scale variance are estimated finally.

GDS. Shape Index was initially proposed for the graphical visualization of the surface [6]. Dorai et al. [7] employed a modified definition for surface point identifying, which is defined as follows:

$$\begin{aligned} {{S}_{I}}(p)=\frac{1}{2}-\frac{1}{\pi }{{\tan }^{-1}}(\frac{{{\kappa }_{1}}(p)+{{\kappa }_{2}}(p)}{{{\kappa }_{1}}(p)-{{\kappa }_{2}}(p)}) \end{aligned}$$
(1)

where \({{S}_{I}}(p)\) ranges during [0, 1], \({\kappa _1}(p)\) is the maximum normal curvature, \({\kappa _2}(p)\) is the minimum normal curvature and \({\kappa _1}(p) > {\kappa _2}(p)\). Shape Index provides a continuous gradation between salient shapes, such as convex, saddle, concave, and it owns a large vocabulary to describe subtle shape variations. Shape Index is appropriate for the description of the distribution of the disk since the value ranges from [0, 1].

Du et al. [8] proposed Geodesic Disk Spectrum (GDS) for the part-in-whole matching task. Let P be a 3D manifold surface of n vertices, a geodesic disk GD(pr) with center p and radius r is defined as follows:

$$\begin{aligned} GD(p,r)=\{v,\text { }gd(v,p)\le r\} \end{aligned}$$
(2)

where gd(vp) is the geodesic distance between any point v to the center p. GDS is obtained by counting the distribution of Shape Index for all the points within a geodesic disk. GDS of one geodesic disk is given by:

$$\begin{aligned} GDS(p,r)=(\frac{{{n}_{1}}}{n},\frac{{{n}_{2}}}{n},...,\frac{{{n}_{N}}}{n}) \end{aligned}$$
(3)

where N is the number of bins range in [0,1] and \({{n}_{1}},{{n}_{2}},...,{{n}_{N}}\) are the number of points that fall into respective bins.

MR-GDS. The MR-GDS descriptor is proposed in this section. The MR-GDS descriptor of the template T depends on the multi-radii geodesic disks, and the multi-radii geodesic disks center on the anchor point q need to be extracted.

Procedures of computing the MR-GDS descriptor are listed as follows. Firstly, the farthest point s away from q on T is found and the distance between distance s and q is set as \({{d}_{\max }}\). Secondly, divide \({{d}_{\max }}\) into \(n-1\) discrete intervals equally as \(\left\{ {{d_1},\mathrm{{ }}{d_2}, \ldots ,{d_n}} \right\} \), and all the geodesic disks with the center q and the radius \({{d}_{i}}\) are extracted. Figure 3 shows a set of geodesic disks with different radii on T. Note that n is usually set as \({{d}_{\max }}/mesh\_resolution\), where \(mesh\_resolution\) is the resolution of the template, which is the average point-to-point distance. This configuration ensures enough samplings of geodesic disks and contributes to accurate scale estimation results. Finally, the GDS descriptor is computed for all geodesic disks with different radii, and they combine into the MR-GDS descriptor together.

Fig. 3.
figure 3

A set of geodesic disks with different radii on T. The red regions are geodesic disks and the radii of the geodesic disks are is increasing from left to right. (Color figure online)

Scale Variance Estimation. Scale normalization is actually the calculation of the scale variance between the fragments and the template. Based on the anchor points obtained in Sect. 3.1, the largest enclosed geodesic disk on F with the center of anchor point is extracted firstly to compute the GDS descriptor. Then, the scale variance will be estimated by comparing the GDS of the fragment and the MR-GDS of the template. The specific procedures of scale variance estimation are listed as follows.

Firstly, extract the largest geodesic disk on F and compute the GDS descriptor of the fragment. Aiming at F, the largest enclosed geodesic disk \(GD\left( p,r \right) \) centered on the anchor point p is needed. All geodesic distances between boundary points and p are computed firstly, and the shortest one is set as the radius r of \(GD\left( p,r \right) \). Note that principle curvatures of boundary points are abnormal from the ones on the template, so the two-ring neighbor points of the boundary are deleted, which eliminate the affections of GDS descriptor calculating. Figure 4 illustrates the extracted geodesic disk on the remained fragment regions. The GDS descriptor is then computed.

Fig. 4.
figure 4

(a) The gray region is two-ring neighbor points of boundary on fragment, which is deleted in order to eliminate the affections of GDS descriptor calculating. (b) The red region is the largest geodesic disk on the remained fragment. (Color figure online)

Secondly, comparing GDS of the fragment and the MR-GDS of the template are conducted. The similarity between \(GDS\left( p,r \right) \) and \(GDS\left( q,{{d}_{i}} \right) \) is given as follows:

$$\begin{aligned} dis=E(GDS(p,r)-GDS(q,{{d}_{i}})) \end{aligned}$$
(4)

Note that the similarities are measured as the variance of differences between each value of GDS. The most similar m geodesic disks are selected as candidate disks and their radii \({{d}_{i}}\) are recorded. The scale variance is computed as the ratio between \({{d}_{i}}\) and r, where \({{d}_{i}}\) is one of the candidate disks. At last, m scale variances are obtained. Figure 5 shows the candidate GDS descriptors of the most similar m geodesic disks, where \(m=3\) in this case. We can see that the GDS descriptors between \(GD\left( p,r \right) \) and the candidate disks are quite similar.

Fig. 5.
figure 5

The three similar GDS compare with \(GDS\left( p,r \right) \), which will be used for the next post-verification strategy.

3.3 Post-verification Strategy Based on ICP

A post-verification strategy based on ICP [9] is proposed to find the most accurate scale variance. The largest geodesic disk of the fragment is geometrically modified according to the m scale variances. Then, the largest geodesic disk of the fragment is registered with m corresponding geodesic disks on the template, and the registration error is achieved using ICP. Among them, the one with the lowest registration error is selected and the corresponding scale variance is selected as the final scale estimation result. Detailed descriptions of this algorithm is shown in Algorithm 1. The transformations between the fragment and the template are also obtained through ICP and the fragment can be matched to the template. The fragments are thus reassembled when multiple fragments are all matched to the template.

figure a

4 Experiments

In the section, datasets used in our experiments are firstly introduced. Experimental results and discussions on scale estimation and fragments reassembly are then presented.

4.1 Datasets

We build two datasets to verify the proposed algorithm. Examples of the two datasets are shown in Fig. 6 on the same perspective. Note that, the Ground Truth is obtained manually and all fragments have different translations compared with the template. The first dataset is used to evaluate the effectiveness of the proposed scale normalization algorithm and the second dataset is used to evaluate the effectiveness of the proposed algorithm for fragments reassembly.

Dataset 1: Since there is no existing benchmarks for 3D shape scale normalization, we constructed dataset 1 to verify the effectiveness of the proposed algorithm. The scanned Terracotta Warriors template is randomly segmented. Geomagic Studio 2012 is used to perform random scale transformation on each fragment and the scale transformation are recorded as Ground Truth.

Dataset 2: In order to verify the practicability of the scale normalization algorithm in the fragment reassembly, we use the artificial Terracotta Warriors as a template for the experiment, while the fragments are still in the original sizes.

Fig. 6.
figure 6

The datasets on the same perspective. (a) shows dataset 1, where the template Terracotta Warriors 3D model is segmented into four fragments randomly and the fragments are modified in random scales. These scales are regarded as the Ground Truth. (b) shows dataset 2, where the artificial Terracotta Warriors 3D model is set as the template and five fragments are with real scales.

4.2 Experimental Results

Scale Normalization. Experiments for scale normalization are conducted on dataset 1. The accuracy of scale estimation is evaluated by the error \(\varepsilon \) and the error ratio \(\delta \), which are defined in the follow. Scale is the experimental results and GT is Ground Truth.

$$\begin{aligned} \varepsilon =scale-GT,\text { }\delta =\frac{scale-GT}{GT} \end{aligned}$$
(5)

The scale estimation results of our algorithm are shown in Table 1, where \(mesh\_resolution\) represents the average point to point distance of the template, scale represents our result, respectively. \(\varepsilon \) and \(\delta \) in Table 1 are both quite small, which illustrates the effectiveness of the proposed algorithm.

Table 1. The scale estimation results. Mesh resolution represents the average point-to-point distance of the template, scale represents our result, respectively.

The time costs are shown in Table 2 for the four fragments. The main time costs can be divided into three parts, which are the computation of GDS for fragments, the computation of MR-GDS of the template, and the post verification using ICP. The time costs are all with the milliseconds level, which is quite fast.

Table 2. The time costs of the experiments. The main time costs can be divided into three parts, which are the computation of GDS for fragments t1, the computation of MR-GDS of the template t2, and the post verification using ICP t3. The time costs are all with the milliseconds level, which is quite fast.

Our algorithm is compared with the scale normalization algorithm proposed by Sahillio\(\breve{g}\)lu and Yemez [5]. Their algorithm all failed in estimating the scales of the two datasets we used. Their method failed in finding three pairs of corresponding points, and can not compute the scale variances through proportional constraints. This illustrates the advantage of our algorithm in dealing with relative small fragments.

Fragment Reassembly. Experiments of fragment reassembly are conducted on dataset 2. The scales of these fragments compared with the template are estimated firstly. Then, the fragments are modified according to the estimated scales in order to conduct the ICP algorithm. At last, using the transformations obtained by ICP, the fragment could be registered to the template. To numerically evaluate reassembly accuracy, the reference error \({{\varepsilon }_{R}}\) and the reference error ratio \({{\delta }_{R}}\) are defined. Let S, \({{S}^{'}}\) be the template model and the reassembled model, respectively. The reference error \({{\varepsilon }_{R}}\) is defined as follow

$$\begin{aligned} {{\varepsilon }_{R}}=\sqrt{\sum \limits _{v\in S,{{v}^{'}}\in \widehat{S}}{{{d}^{2}}(v,{{v}^{'}})}/|S|} \end{aligned}$$
(6)

where |S| represents the number of points of S, d(; ) represents the Euclidean distance. The reference error ratio \({{\delta }_{R}}\) is measured as a percentage of the diagonal length of bounding box.

The results of fragments reassembly by our algorithm are shown in Fig. 7. The scanned fragments are shown in Fig. 7(a) and (b) shows the reassembly results. The error distribution of each point is shown in Fig. 7(c), where the red color represents a large error and the blue color represents the opposite. It can be intuitively seen that the error is relatively small.

Fig. 7.
figure 7

The fragments reassembly results of Terracotta Warriors. (a) The scanned fragments. (b) Reassembly results. (c) The error distribution of each point, red represents a large error while blue represents the opposite. (Color figure online)

Quantitative results are shown in Table 3, including reference error \({{\varepsilon }_{R}}\) and the reference error ratio \({{\delta }_{R}}\). It can be seen that our algorithm obtains a high accuracy of reassembly.

Table 3. Quantitative results of the fragment reassembly. \({{\varepsilon }_{R}}\) is the reference error and \({{\delta }_{R}}\) is the reference error ratio.

5 Conclusions

In this paper, we propose a scale normalization algorithm based on MR-GDS for archaeological fragments reassembly. The MR-GDS descriptor is firstly proposed which integrate the scale information into the descriptor. With only one pair of matched points, the GDS descriptor of the fragment and the MR-GDS descriptor of the template could be computed. The scale variance can thus be estimated by comparing the GDS descriptor and the MR-GDS descriptor. In order to improve the accuracy of scale estimation, a post-verification strategy based on ICP is also proposed. Experiments are conducted for both scale normalization tasks and template-based fragments reassembly tasks, which all achieved accurate results. Compared with other algorithms, our algorithm needs only one pair of matched points and is suitable for fragments with few extremities, which expands the application scenes greatly for template-based archaeological fragments reassembly tasks.