The Journal of the Society for Art and Science, Vol. 23, No. 2, pp. 4:1 – 4:17 (2024)

A New Matching Algorithm for Stone Tool Reassembly Based on
Contour Points of Flake Surface
Amgalan Altansukh1) Mengbo You1) Enkhbayar Altantsetseg 2 )Oyundolgor Khorloo 2 ) Fumito Chiba 3 ) Kouichi Konno 1)
1)Graduate School of Science and Engineering, Iwate University
2) National University of Mongolia3) LANG CO.,LTD.
amgaa0809(at)gmail.com, {ymb, konno21}(at)iwate-u.ac.jp

Abstract
Reassembling fragmented stone tools from the Paleolithic era remains a signiﬁcant archaeological challenge.
Our study addresses this challenge by introducing a matching algorithm with a primary emphasis on enhancing
computational eﬃciency and reﬁning the ﬂake surface matching process. Our method builds upon previous
study, speciﬁcally targeting reductions in the computation time. The algorithm was tested on a dataset comprising 43 stone models. A critical aspect of this study is ﬂake surface matching which is a fundamental aspectof stone tool reassembly. By optimizing the computational cost, we aim to provide archaeologists with a more
eﬃcient and accurate tool for reconstructing archaeological artifacts.
Keywords: Flake surface, Matching algorithm, Point cloud, Reassembly, Stone tool restoration– 4:1 –

The Journal of the Society for Art and Science, Vol. 23, No. 2, pp. 4:1 – 4:17 (2024)

1 Introduction
The Paleolithic era witnessed the production of awide range of stone tools, including cutting implements and weapons. Several archaeologicalsites have been excavated and analyzed. Thesesites often contain remnants of debris generated
during the production and utilization of stone
tools. Considering their durability compared toorganic materials, such as bones, antlers, andwood, stone artifacts provide signiﬁcant evidenceof the locations and periods of early human activities because of their geographical distribution,
and adaptive capacities in various environments[1].
Stone tools oﬀer valuable insights for early human toolmakers, including their production techniques, life patterns, and transfer of things [2].
For instance, archaeologists analyze the assembly order for investigating creation of stone tools.
Figure 1 shows an example of excavated stones,
which were restored. The restored mother rockis referred to as the “Joining Material.” Joining materials have played a signiﬁcant role in theanalysis of human activities during this period.
Moreover, these joining materials have a significant educational value and can be displayed inhistorical museums and educational institutions[1].
However, the completion and analysis of relics[1] is challenging. Even though relics are subject to measurement through scanning methods[3] and several issues can be solved swiftly andaccurately by computational algorithms, the taskof reassembling stone tools remains intricate anddemanding.

Figure 1: Example of excavated stones.

Archaeologists manually reassembled stone
tools using traditional restoration techniques, as
shown in Figure 2(a). However, this approachhas two major limitations. First, the assembly
process is time-consuming because it requires theidentiﬁcation of correctly matched surfaces bytrial and error. Second, the assembly and disassembly of stone tools require special archaeological knowledge. These are challenging tasks,
because not all stone tools constructed from themother rock can be excavated.

(a)

(b)

Figure 2: a) Traditional reassembly.
b) Making stone tools.

In recent years, several reassembly approaches
have been developed for the reconstruction offragmented archaeological artifacts such as pottery and fresco wall paintings. However, the successful application of these methods to reassemble
stone tools is limited for owing to the following
reasons: (1) Irregular shapes of stone fragments:
stone tools exhibit highly irregular shapes, making it challenging for descriptor-based methods
to extract crucial features. The lack of distinct regional features on ﬂake surfaces further
complicates the feature extraction process. (2)Global matching algorithms have proven inadequate, considering the requirement for partial
matching between pairs of ﬂake surfaces. The
previous study [16] introduced a method of stone
matching. [16] extracts the boundary of the ﬂake
surface to determine the initial position. However, [16] the excessive computation time and difﬁculties with partial matching.
Therefore, our study introduces a new ﬂake
surface matching algorithm designed to signiﬁcantly improve the computational time while accurately accommodating partial matches. This
study aims to address these issues. The proposed method involves aligning two ﬂake surfaces– 4:2 –

The Journal of the Society for Art and Science, Vol. 23, No. 2, pp. 4:1 – 4:17 (2024)

based on the maximum number of matched contour (edge) points. By introducing a tunable
matching parameter, a user can control the balance between the computation time and accuracy,
thus addressing a critical bottleneck in stone tool
reassembly.
The contribution of this study is a new ﬁnetuned approach to ﬂake surface matching thatutilizes the contour points of ﬂake surfaces. Thisapproach not only enhances the eﬃciency of thematching process but also proves eﬀective inhandling the unique challenges posed by partial matching scenarios. We conﬁrmed that our
method was eﬀective for reassembling stone tools.

2 Related Work
2.1 Stone Tools
Numerous stone tools have been manufactured
prehistorically through rock striking process.
Crafting a stone tool involves the repetitive action of striking the edge of a rock with a pebble,
resulting in the peeling of ﬂakes to modify andreﬁne the tool shape, as shown in Figure 2(b).
Within this archaeological context, an assemblyof stone tools, known as joining materials, represents a reconstruction in which the core andﬂakes originate from the same rock. In addition,
the joining materials are dispersed across the excavation site. Restoration involves reassemblingthese fragments to recreate the original composition and structure of a stone tool [1].

2.2 Reassembly Methods
Extensive studies have been conducted on the reassembly fractured objects. Huang et al. [4] introduced a feature-based alignment method thateﬀectively addressed these scenarios and achieved
remarkable results. However, their method iscomplex and, comprises various specialized algorithms for tasks such as segmentation, multiscale feature extraction, correspondence determination, registration, collision detection, and supervised learning. Consequently, implementation
and adoption of these methods pose signiﬁcantchallenges. Brown et al. [5] proposed a method

speciﬁcally designed to reassemble fragments ofwall paintings. Willis et al. [6] presented a system tailored to pottery shard reassembly.
In [7], an innovative method was proposed
for reassembling axially symmetric ceramic potsfrom three-dimensional (3D) scanned fragments.
Their method addressed speciﬁc challenges withceramic artifacts, particularly axial symmetry.
Their approach is noteworthy because it incrementally adds fragments using a beam searchtechnique. This signiﬁcantly mitigates the false
positive matches. Their method not only improves the match accuracy through geometricdescriptors but also facilitates the simultaneous
reassembly of multiple pots from mixed collections. However, the application of their techniqueto the reassembly of stone tools is challenging.
Reassembly stone tools lack the axial symmetry
present in ceramic pots and require considerationof the ﬂake removal sequence from the core stone.
In [8], a new methodology for generating synthetic 3D fragmented data was proposed to facilitate the evaluation of object restoration, particularly in the context of cultural heritage. Theirapproach was designed to overcome the challengesof limited availability of real test data. Theirmethod produces artiﬁcial fracturing from an input object without physical simulation. In addition, their method uses a cutter object to create
new breaking edges and thereafter fragments theobject. It generates a large-scale fragment test
dataset from existing cultural heritage models.
These datasets have advantage of ground truth(the input object before fracturing), which is often missed. However, it does not directly relate
to the reassembly stone tool task or reassemblytask. It focuses on the reverse process of reassembling models.
In [9], a mesh-based approach to create restorations integrated with broken objects was introduced. Utilizing the 3D scanned meshes ofboth broken objects and their intact counterparts, their method generates a restoration pieceby smooth. It focuses on objects for which themissing parts are predictable, and its goal is torestore the original form of the objects. However,its application to the reassembly of stone tools
poses several challenges. First, it relies on the– 4:3 –

The Journal of the Society for Art and Science, Vol. 23, No. 2, pp. 4:1 – 4:17 (2024)

predictability of the shape of the missing part.
This condition is not satisﬁed in stone tool reassembly in which the shape and size of the ﬂakes
are unpredictable. Second, the methodology ismesh-based, whereas stone data were represented
in a point cloud. Finally, their proposed method
does not address partial matching or the sequential order of assembly, both of which are crucial
for stone tool reassembly tasks.
The Fantastic Breaks dataset [10] was designed to facilitate machine learning study in automated reassembly by providing a comprehensive collection of 3D scans of real-world broken
objects and their counterparts. It focuses on acomplete object reconstruction using predeﬁned
models. The use of predeﬁned models may notalign with the stone tool reassembly tasks. This
often requires addressing the unpredictable fragment shapes and sizes. Machine learning approaches for stone tool reassembly may face challenges owing to the requirement for extensive,well-prepared datasets and the computational demands for processing large datasets. Stone tool
datasets are often characterized by a lack of extensive pretested data. A highly detailed point
cloud representation can limit the applicabilityof machine learning methods. Speciﬁcally, ourdatasets with 280,000 points per stone can strain
the GPU resources [11]. This renders the learning
process more challenging.

2.3 Matching Methods for Flake Surface
Iterative Closest Point (ICP) [12], Super4PCS[13], fast point feature histograms (FPFH) [14]and random sample consensus (RANSAC) [15] algorithms are commonly employed in point cloud
registration tasks. Furthermore, these methods
have been adapted and applied to surface matching scenarios such as ﬂake surface matching.
The ICP algorithm [12], which is one of the cornerstones of point cloud registration, has demonstrated wide-ranging eﬀectiveness in aligning 3D
surfaces. However, it may be diﬃcult to obtain an accurate initial alignment, particularly
in cases involving partial ﬂake surfaces in stone
tools, where achieving a precise alignment can beparticularly challenging.

Super4PCS [13], recognized for its proﬁciency
in matching partial 3D shapes, eﬃciently identiﬁes local surface correspondences. Considering this capacity, it is a compelling candidate foradapting to the task of matching partial ﬂake surfaces. The inherent ability of the algorithm tohandle partial matching aligns well with the result of the present study.
The FPFH [14] calculates the correspondences
and correspondence probabilities, providing valuable guidance for subsequent RANSAC [15]matching algorithms. In [15], a robust method for
estimating transformation models between data
points is presented, particularly in the presenceof outliers, and is a crucial step in aligning partial
ﬂake surfaces. Their approach, which utilizes the
FPFH [14] for the initial correlation, addresses
challenges related to partial matching and aligning surfaces with potentially irregular shapes.
Our previous study [16] focused on the sametask, considering the shape of the ﬂake surface
by utilizing contour points, which is similar to
our approach. However, it was characterized byprolonged computation times and limitations ineﬀectively matching partial ﬂake surfaces.

3 Proposed Method
3.1 Overview
The input of our method comprised a collectionof point cloud of stone tools, acquired throughmeasurements, as described in [3]. Our proposed
reassembly method builds upon the pipeline outlined in [16], with a primary focus on reducing
the computation time and addressing limitations.
The pipeline is executed using the following procedure:
1. Extracting ﬂake surface and identifying contour points.
2. Finding the best ﬂake surface among the candidate ﬂakes. For each ﬂake surface of thecore stone, a matching algorithm is run todetermine the optimal matching ﬂake surface. This operation calculates a transformation matrix and iterates the process until– 4:4 –

The Journal of the Society for Art and Science, Vol. 23, No. 2, pp. 4:1 – 4:17 (2024)

the best matching surface is identiﬁed within
the candidates.
3. Transformation of matched ﬂake. All ﬂake
surfaces associated with an unique ﬂake
stone of the matched surface are transformedby the transformation matrix.
4. Reconstruction.
After integrating thematched ﬂake stone with the core stone, theﬂake surfaces of the core stone are thereafterreconstructed.
5. Iterative matching. The matching process isiterated until all stone tools are successfully
matched.
In contrast to a previous study [16], our study
introduces several pivotal enhancements. Initially, the region growing algorithm [17] is employed to extract ﬂake surfaces. However, thealgorithm challenges in capturing the boundary
points of ﬂake surfaces owing to its sensitivityof normal vector variation and curvatures near
sharp edges. To address this shortcoming, we
subsequently perform a boundary correction stepthat has not been used in the previous study. In
addition, our matching algorithm features a tunable matching parameter that enables the optimization of ﬂake surface matching for enhanced
eﬃciency. This parameterized approach diﬀers
signiﬁcantly from that used in the previous study.
Although the D2 algorithm [18] is employed to select candidate part in the previous study, we decide to employ our matching algorithm by tuning
the parameter for candidate ordering, as our candidate ordering improved without D2. Finally,
our method performs matching and reconstruction operations directly on the point cloud, as
opposed to the mesh-based approach used in theprevious study. This shift enhances the computational eﬃciency.

3.2 Extracting ﬂake surface andidentifying its contour points
The region growing algorithm segments a point
cloud starting from a seed point and adding
neighboring points that satisfy the criteria of geometric features such as diﬀerences in point nor-mals and curvatures. Near sharp edges, the large
amount of shape variation in these local geometric features can result in over-segmentation nearthem. This is largely owing to the algorithm’s
sensitivity to abrupt changes in normal vectors
and curvatures.
In the region growing segmentation, the angle
threshold α limits the addition of point to thosewith similar normal vectors in a region, ensuring the smoothness of the region. The curvature threshold c excludes points from the region
with abrupt changes in curvatures. The minimum number of points l is the required numberof points for a region to be considered valid.
The algorithm categorizes segments as valid orinvalid based on the number of points they contain. Regions that are over-segmented, particularly near sharp edges, may be classiﬁed as invalid owing to their small size. If the numberof points belonging to a segment is less than l,points in these invalid segments are consideredunsegmented points, as shown in the gray points
in Figure 3. These unsegmented points are crucial, because they contain essential information
regarding the true boundaries of the ﬂake surfaces. Therefore, our method corrects the boundary of the region by merging unsegmented points
as a re-segmented operation.

(a)

(b)

Figure 3: a) Result of the region growing with anangle threshold α set to 1.5◦ , a curvaturethreshold c to 0.1 and the minimum number of point in a valid segment l to 500.
b) A case of the unsegmented point withnormal vector.

To correct the boundaries of the ﬂake surfaces,
The points were re-segmented using ﬁve steps.
The ﬁrst two steps are introduced in Point Cloud
Library [17]. Steps from 3 to 5 are additional
steps to correct boundary points of the segmentation process. The detail of steps is as follows:– 4:5 –

The Journal of the Society for Art and Science, Vol. 23, No. 2, pp. 4:1 – 4:17 (2024)

1. Running the region growing segmentation
with the angle threshold α and curvaturethreshold c.
2. The segments are categorized into twogroups based on the number of points. Regions with a point count exceeding l are considered as valid segments, represented by the
purple, blue, green, and yellow points in Figure 3(a). Conversely, those with fewer points
are identiﬁed as over-segmented regions that
included unsegmented points, which are depicted as gray points in Figure 3(a).
3. Find the k nearest neighbors of each unsegmented point. The k is determined that atleast one point included in the valid segment
is selected as the nearest neighbor point.
4. For each unsegmented point:(a) Determine if the neighboring points
belong to a valid segment. If a neighboring point belongs to a valid segment, thatsegment is considered a bordered segmentwith the unsegmented point. Thus, the unsegmented point is a candidate of merging
to bordered segment. This scenario is visualized with the unsegmented point as alight gray point and the neighboring point
as a light purple point in Figure 3(b). The
neighboring point is considered a border
point. When two or more neighboring points
are found in the same valid segment, selectthe nearest neighboring point as the border
point.
(b) Calculate the angle ϕ between thenormal vectors of border point invalid segment and the unsegmented point, illustrated
as m and n in Figure 3(b), respectively. This
calculation is computed by Equation (1).


n·m
(1)
ϕ = arccos|n||m|(c) Merge the unsegmented point into
the bordering valid segment with the minimum degree of ϕ.
5. Repeat step 4 until all unsegmented point
have been merged into valid segments.

After extracting the ﬂake surface using the region growing segmentation with the boundary
correction step, the next step involves identifying the contour points of the ﬂake surface. First,
a ﬁtting plane is computed for the ﬂake surface
using Principal Component Analysis [19] to ensure an accurate representation. Subsequently,
all the points within the ﬂake are mapped onto
the ﬁtting plane, thereby eﬀectively aligning themwith the surface morphology. Finally, a twodimensional concave hull algorithm [20] is executed. This crucial step aids in the precise identiﬁcation of contour points in the clockwise direction and oﬀers valuable insights into the characteristics of the ﬂake surface.

3.3 Reassembly
Lithic materials exhibits distinct characteristics
representing the joining surface shape and separation from the reassembly of other fractured
objects [16], as shown in Figure 4. The ﬁrst distinguishing property is the sequence of ﬂake generation from a single core [21]. Unlike arbitrary
ﬂake matching, stone tool reassembly requires asequential order. For instance, if ﬂake A precedes
ﬂake B in the peeling process, then, considering
the sequential extraction of ﬂakes from the corestone, the matching procedure entails a reverse
order.

Figure 4: Case of making a stone tool (ﬂake A ispeeled ﬁrst followed by ﬂake B).

The second property pertains to the ﬂat andsmooth of most ﬂake surfaces in the matching.
This is because of the capacity of the motherstone to divide sharply [1]. Consequently, conventional matching algorithms that rely on surface features may not be suitable.
The third property arises when a single ﬂake
surface can fracture into several fragments. As
shown in Figure 4, the matching of ﬂake A re-– 4:6 –

The Journal of the Society for Art and Science, Vol. 23, No. 2, pp. 4:1 – 4:17 (2024)

quires prior matching of ﬂake B with the corestone to create ﬂake surface Fc , which is a composite derived from ﬂake surfaces F1 and F2 .
Consequently, the matched ﬂake surfaces requireda reconstruction process to identify the subsequent surface.
In our approach, stone cores are designatedmanually, initiating the matching process on the
ﬂake surfaces of the core stones. As our dataset
comprised multiple core stones, they are reassembled in succession. Each ﬂake surface of the corestone is matched with every ﬂake surface of theﬂakes to identify the best match. To enhancethe eﬃciency and reduce the number of matchingtasks, the candidate order for each ﬂake surface
is calculated. This order is sorted based on thecorresponding candidate score.
3.3.1 New ﬁne-tuned ﬂake surface matchingalgorithm
A ﬁne-tuned matching algorithm is developed fortwo purposes during reassembly. It is designed to
calculate both the candidate and matching scoresin pairwise matching. The matching algorithm
calculates a matrix Mst , by which a source surface in the candidate stone is mapped onto a target surface in the core stone. The algorithm calculates both the candidate and matching scoresby tuning the matching parameter n, which determines the level of detail during the matchingprocess.
Figure 5 shows the diﬀerence between the candidate and matching scores for diﬀerent values ofn. For instance, when the matching parametern is set to 6, the algorithm works on all combinations between every 6 points in the target surface and source contours. When it is set to 12
as an example, the algorithm works on all combinations between every 12 points in the target
and source contours. In this case, the number of
combinations is reduced, and the working process
faster than n is set to 6. In contrast, the accuracy tends to decrease depending on the contour
shape. Therefore, n is determined by considering
the balance of speed and accuracy.
Our reassembly method requires less computational time to calculate the candidate score. How-

Figure 5: Tuning the precision parameter.

ever, the matching accuracy must be high whenreﬁtting a candidate for the core stone and calculating the matching score.
A ﬁve-point set {pi , pi+n , pi−n , di , ct } is constructed, where pi is i-th contour point of thetarget surface t, n is the matching parameter,
pi+n and pi−n are separated from pi by n, ct
denotes the center point of the target surface,and di denotes the centroid point of the triangle
whose vertices are pi+n , pi+n , and ct . Similarly,
another ﬁve-point set {pj , pj+n , pj−n , dj , cs }is constructed in the same manner on the source
surface s, as shown in Figure 6.
The normal vector ki of the triangle, control
vector mi , and centroid di of the triangle are calculated using Equations (2), (3), and (4), respectively. All the points and vectors are constructedin the same manner on the source surface. The
3D rotation matrix Rij is determined by satisfying Equation (6) owing to ki ⊥ mi and kj ⊥ mj .
The 3D translation matrix Tij is calculated using
Equation (5). The overall transformation matrix
Mij is thereafter computed using Equation (7),integrating both rotation and translation to facilitate transformation in a 3D space.

– 4:7 –

ki = (ct − pi−n ) × (pi+n − pi−n )(2)di = (pi+n + pi−n + ct )/3(3)

mi = d i − p i(4)di = Tij · dj

(5)(ki = Rij · kj
mi = Rij · mj

Mij = Tij · Rij

(6)

(7)

The Journal of the Society for Art and Science, Vol. 23, No. 2, pp. 4:1 – 4:17 (2024)

Figure 6: Mapping a source into a target surface.

Equation (8) computes error metric e at point
pi , which is the contour point of the target surface and determines the minimum distance between pi and any contour point p′j on the transformed source surface by Mij . The number ofmatched target contour points rij is calculated using Equation (9), where Es denotes the maximum
edge length in the concave hull [20]. The matching transformation matrix Mst , which aligns thesource surface s to the target surface t, is selected based on max rst , which is the highest
value of rij as shown in Equation (10). Algorithm 1 outlines our ﬂake surface matching algorithm; nt and ns denote the respective numbersof contour points; ct and cs denote the center
points; n denotes the matching parameter in the
input; Mst denotes the matching transformation;
max rst denotes the highest number of matchedcontour points in the output.
e(pi ) = min(|pi , p′j |);

rij =

X
i+=n

(

j, i+ = n1,

if e(pi ) < n · Es ,0,otherwise

max rst = max(rij )

(8)

(9)

(10)

The uniform point density of the dataset enables a straightforward representation of thematched area using the number of corresponding points. In this context, qi denotes any pointon the target surface t, whereas qj is a point on
the source surface s that is nearest point to thetangent plane of qi .

The two orthogonal distances can be calculated
by surface s and t. One is the distance from thetangent plane of qi to qj , and the other is the
distance from the tangent plane of qj to qi . In
our method, if both distances are less than dc andthe distance between qi and qj is also less thandc , qi is identiﬁed as a correspondence point of
qj . dc is determined 1.5 by experiment.
Subsequently, matched point pairs of surface s
and t are derived and the matched surface percentages are determined using Equation (11). In
this equation, Pu denotes the matched area percentage of the surface u, Nu denotes the numberof corresponding points of the surface u, Tu denotes the total number of points on surface u.
Speciﬁcally, Pt and Ps are the instances of Pu forthe target surface t and s, respectively. Pt and Ps
are indicated the percentage of each matched surface. Moreover, the matching score, denoted as
Sst , or candidate score, is computed using Equation (12), where max rst denotes the number ofmatched contour points, and Pt and Ps represent
the matched surface percentage of the target andsource surfaces, respectively.
A higher Sst value indicates a greater level ofmatching between the target t and source s surfaces, implying that there are more matched contour points and a larger matched area between
them. Conversely, a lower Sst value implies alower matching level, indicating fewer matchedcontour points and a smaller matched area between the surfaces.

– 4:8 –

Pu =

Nu
Tu

(11)

The Journal of the Society for Art and Science, Vol. 23, No. 2, pp. 4:1 – 4:17 (2024)

Sst = max rst

p

Pt Ps

(12)

3.3.2 Flake surface reconstruction
Algorithm 1: Flake Surface Matching
Algorithm
Input : contour points of the target,
nt , ct , contour points of thesource, ns , cs , n
Output: Mst , max rst
M′ ← ∅;
for i ← 0 to nt dopi ← i-th contour point of the target;
for j ← 0 to ns do
pj ← j-th point of the source;
Mij ←getTransformMatrix(pi , ct , pj , cs , n);// by Equation
7′′S
M ← M Mij ;j + = n;end
i + = n;end
rst ← 0;foreach Mij in M′ do
S ′ ← {p′j | p′j = Mij · pj , pj is any
contour point of the source };rij ← 0;for i ← 0 to nt dopi ← i-th contour point of thetarget;e = min ∥pi − p′j ∥, p′j ∈ S ′ ;if e < Es · n thenrij + = 1;
end
i + = n;end
if rij > max rst thenmax rst ← rij ;Mst ← Mij ;end
end

The process of reconstructing the original ﬂake
surface involves the detection and integration ofdivided ﬂake surfaces. Figure 7 shows a visual representation of the reconstruction, with
the ﬂake surfaces Fa in blue and Fb in green.
The reconstruction proceeds through the following steps:

Figure 7: Flake surface reconstruction.

1. Each pair of ﬂake surfaces Fa and Fb frommatched ﬂakes undergoes a search for theclosest contour point pbj of Fb for each contour point pai . If the distance between pai
and pbj is less than the threshold named dr ,the pair is designated as a corresponding contour point pair, red points in Figure 7.
2. Calculate the angle between na and nb , average normal vector of Fa and Fb , respectively. In addition, the number of corresponding pair is counted.
3. The corresponding pair derived by Step 2 ishigher than the threshold w and the angle
between na and nb is less than the thresholdwθ , then the ﬂake surface pair Fa and Fb ismerged into a single ﬂake surface Fc .

4 Results and limitation
The implementation was performed on a PC withan Intel Core i7-10700 CPU and 16 GB memory. The experiment involved testing data from– 4:9 –

The Journal of the Society for Art and Science, Vol. 23, No. 2, pp. 4:1 – 4:17 (2024)

43 stone models, as shown in Figure 8. The data
could be reassembled into three groups. N o.01
and N o.20, shown in thick red box are the corestones and Group 1 and 2, respectively. The stone
tools belongs to Group 3 are role of dummy data
for assembling Group 1 and Group 2. There were12,300,222 points, with an average of 286,052
points per stone.

Figure 8: Experiment data.

4.1 Experimental Results
This section describes the evaluation and experimental results. First, we examined an evaluation
score and execution time according to the tunable matching parameter n, when n is changed
from 1 to 12. Second, we compared our methodwith related methods. Thereafter, we presented
experimental results. Finally, we compared our
results with the previous study [16].
To assess the accuracy of the resulting shape,
we calculated the evaluation score which is how
the matching surface shape coincides. To achievethis, Equation (12) which indicates the surface
coincident equation between surface s and surfacet, was used as follows:
Se = Sst n=1(13)

where Se means score of the surface coincident.
To evaluate how to match the surface shape
Equation (13) was applied to compare the resulting shape of other methods. In Equation (13), we
set n to 1 to be fair to the shape evaluation when
we evaluate the resulting shape. In other words,

all points belonging to surfaces s and t are used
to evaluate the resulting shape. If Se is high, itmeans that the shape is more match.
Figure 9(a) illustrates that the evaluation score
Se for the four matching pairs when the matching parameter n varied from 1 to 12. Figure 9(b)shows the computation time for the same parameter variations. Figure 9(a) shows the highest
scores for each pair, indicated by red circles. As
shown in Figure 9(a), three of the four had highest scores when parameter n was set to 6. In contrast, Figure 9(b) shows the computation timeof the four matching pairs. The average computation time of each matching pair for parameter n =6 was 1.52, and that for n =4 was 6.87
s. This shows that when n =6, the computation time is approximately 4.5 times faster thanwhen n =4. Therefore, to balance computation
time with accuracy—as reﬂected by the evaluation score—in the performance of the algorithm,
we set n =6 for the comparative analysis.
We compare our method and previous works.
The comparison involved benchmarking ourmatching algorithm against Super4PCS [13] andFPFH-SACIA [14] with 100 iterations employedfor partial matching tasks. For a fair comparison
across all methods, we applied Se , the evaluation
score, calculated by Equation (13).
The ground truth was manually aligned to imitate the real reassembly results achieved by archaeologists, thereby providing a benchmark forevaluation. Table 1 summarizes the results ofthis comparison. In this table, column “t” denotes the ID of the target surface, formed by
combining the stone identiﬁcation and ﬂake surface ID(e.g.,“01 0” indicates core stone ID “01”and ﬂake surface ID “0”). Similarly, column “s”
refers to the source ﬂake surface ID, following
the same identiﬁcation pattern. “S4PCS” represents the absolute diﬀerence between the evaluation score of Super4PCS and the evaluation
score of the ground truth. “FPFH” and “Our”
indicate the same manner. These diﬀerences arecalculated using Equation(14), where D denotes
the absolute diﬀerence, Sgt denotes the evaluation score of the ground truth, and Se denotes
the evaluation score of the respective algorithm.
In this table, the proposed ﬂake surface matching– 4:10 –The Journal of the Society for Art and Science, Vol. 23, No. 2, pp. 4:1 – 4:17 (2024)

algorithm yields the closest evaluation score to
the ground truth across all the target and source
pairs. The performance is graphically presented
in Figure 10. These comparison results highlightthat our algorithm performed better than existing
methods for ﬂake surface matching, particularly
partial matching.
D = |Sgt − Se |(a) 01 1 vs 06 1.

(14)(b) 01 5 vs 10 3.

Figure 10: Visual representation of comparison result.

(a) Relationship between n and the evaluation score Se .

(b) Relationship between n and the computation time.

Figure 9: Parameter sensitivity test.

Table 1: The absolute diﬀerence between the
ground truth and the evaluation score.
t
01 14
01 1
01 0
01 5

s
12 1
61
07 0
10 3

S4PCS
331.4
142.9
248.3
162.8

FPFH
374.6
115.8
157.6
157.2

Our
1.7
19.1
21.5
6.5

We describe the experimental results of our
method. The segmentation process for extracting the ﬂake surfaces requires the estimation ofpoint normal vectors and point curvatures. Previous studies [22] and [19] are employed for normal vector estimation and point curvature calculation, respectively. During segmentation described by Section 3.2, the parameters are tuned
as follows: the angle threshold α is set to 1.5◦ ,curvature threshold c is set to 0.1, minimum number of points for a valid segment l is set to 500.
The number of nearest neighbor k is set to 4 in theboundary correction step. The dataset in Figure
8 is segmented, extracting 311 ﬂake surfaces andaveraging 7 surfaces per model for the matchingprocess. Figure 11 shows the segmentation results for point cloud N o.10. Figure 11(a) displays
the point cloud with point normal vectors, while
Figure 11(b) shows the outcome of region growing segmentation using [17], where unsegmented
points are shown in gray color, resulting in 67962
unsegmented points of 276336. In contrast, Figure 11(c) shows our boundary correction result,which extracts 8 ﬂake surfaces, and all points aresegmented. In Figure 11(b), (c), and (d), the various colors represent diﬀerent ﬂake surfaces, andthe contour points are highlighted in red. Figure
12 (a) and (b) show the results of the segmentation for stones N o.06 and N o.07, respectively.
The left side of each ﬁgure shows the result of the– 4:11 –The Journal of the Society for Art and Science, Vol. 23, No. 2, pp. 4:1 – 4:17 (2024)

region growing, and the right shows the result of
the boundary correction step. As shown in Figure
11 and 12, our method enables to correct boundary.

(a)

(b)

(c)

(a)

(b)

(c)

(d)

(d)

Figure 11: Result of the boundary correction.
a) Point cloud.
b) Result of region growing segmentation.
c) Result of our boundary correction step.
d) Result of extracting contour points.

Figure 13: Result of ﬂake surface matching.
a) Stone models N o.01 and N o.06.
b) Correspondence between ﬂake surfaces
01 1 and 06 1.
c), d) Result of matching in diﬀerent
views.

(a)(b)

Figure 12: Result of the segmentation.
a) N o.06.
b) N o.07.

Figure 13 shows the matching results between
the two ﬂake surfaces 01 1 and 06 1. The correlation between the two surfaces is shown in
Figure 13(b). Green points represent correlation
points, red points represent unmatched surface
points and yellow points represents matched contour points, and blue points indicate unmatchedcontour points on the target surface.
Figure 14 shows a matching scenario in which
ﬂake surfaces require reconstruction for alignment. Three ﬂake surfaces of stone tools N o.10,
N o.06 and N o.16 were reconstructed to achievealignment with the surface of stone N o.17.
Tables 2 and 3 summarize the speciﬁc detailsof matching Groups 1 and 2. In these tables,

columns “t” and “s” maintain the same format
as that in Table 1, denoting the respective target
and source ﬂake surfaces. The “Order” column
signiﬁes the order of ﬂake surface for the best
matching, determined by the candidate score.
“Time” represents the duration required for computing the ﬁtting transformation matrix for eachpair. For instance, ﬁnding the best match for
ﬂake surface 01 1 in Group 1 requires 44.7 s.
On average, each matching process uses approximately 7.45 s per pair, considering 6 (candidate
order) matching times. All ﬂake surfaces initiatedmatching from stone core N o.01 and N o.20. During the experiment, the maximum edge length inconcave hull Es was set to 0.5 and the matching tune parameter n was set to 6 for matching,
16 for ordering the candidates. Additionally, wedeﬁned two distance thresholds, dc and dr , both
set to 1.5. The threshold w is set to 30, and theangular threshold wθ is set to 15 degrees.
Group 1 resulted in the reassembly of 16 stone– 4:12 –The Journal of the Society for Art and Science, Vol. 23, No. 2, pp. 4:1 – 4:17 (2024)

(a)

(b)

(c)

Figure 14: Reconstruction result.
a) Reconstructed ﬂake surfaces of N o.6,
N o.10 and N o.16.
b) Source ﬂake N o.17.
c) Result of matching.

tools, whereas Group 2 was reassembled of 16
stone tools. Flake stone N o.14 remained unmatched using our method, as explained in detailin Section 4.2. Figure 15 shows the ﬁnal matching
results for both groups, accompanied by imagesof the manually matched imitations.
Table 4 summarizes the comparison result ofour method against our previous study [16] fortwo groups of stones. In the table, “N.Dis” indicates the normalized distance measurement inmillimeters, introduced in previous studies [16],[23]. This measurement method is utilized to
identify the most compatible ﬂake surfaces intheir studies. It provides a standard to measure
the diﬀerence between two ﬂake surfaces on a unit
area. Essentially, a smaller normalized distance
between the surfaces indicates a superior match,signifying a closer similarity or better alignment
between two ﬂake surfaces. The column “Time”
records the computation time for each method inseconds. The computation time for the previous
method is normalized to account for the diﬀerence in CPU. Our study was tested on the Intel
Core i7-10700, while the previous study utilized
an Intel Core i7-4790. The normalization factorof 2.68, derived from the relative ﬂoating-point
math speed of the two CPUs [24], is applied tothe computation times of the previous method toutilize a fair comparison.
In Group 1, our method completed the task

in 400 seconds, whereas the previous study required 5912 seconds. In Group 2, our method
required 290 seconds, compared with the previous 1641 seconds. In the normalized distancemeasurements, our method demonstrated an average of 0.012 mm for Group 1, in contrast to
0.033 mm reported in the previous study. Similarly, for Group 2, our method achieved 0.016
mm, compared with 0.027 mm in the previous
study. Moreover, we successfully addressed thelimitation related to partial matching, allowing
us to match 4 more stones in Group 2 and ultimately reassemble all the stones in that group.
Table 2: Reassembly of Group 1
t
01 14
01 1
01 4
01 0
01 49
01 43
01 55
01 39
01 57
01 53
01 59
01 66
01 63
01 64
01 68

s
12 0
06 0
10 6
07 0
05 1
16 0
11 3
17 0
02 0
18 1
13 1
08 3
19 0
15 0
04 1

Order Score
1/260 417.71
6/250 188.40
2/243 306.41
1/234 459.82
1/225 378.01
2/218 219.58
2/209 175.32
1/205 578.64
4/203 140.52
1/193 281.97
1/187 435.26
1/185 554.52
4/178 346.89
3/172 315.48
6/166 147.44
Total time:Time(s)
27.22
44.70
20.96
37.23
10.37
29.09
19.08
29.75
15.20
37.58
25.39
11.05
57.80
23.83
11.11
400.36

4.2 Limitation
Figure 16 shows a scenario in which proposed
method could not correctly match for stone
N o.14. Red circles indicate the correspondenceof the correct matching. Matched contour pointsof the two ﬂake surfaces are indicated in yellow,
whereas the unmatched contour points are presented in blue. Our matching method is designed
for alignment based on the maximum number ofmatched contour points max rst . However, in
this case, max rst did not achieved the correct
matching. Therefore, our method encountered– 4:13 –The Journal of the Society for Art and Science, Vol. 23, No. 2, pp. 4:1 – 4:17 (2024)

Table 3: Reassembly Group 2
t
20 10
20 9
20 10
20 5
20 2
20 7
20 12
20 18
20 33
20 5
20 38
20 43
20 0
20 44
20 54

s
34 1
27 2
36 0
33 6
25 2
31 1
26 0
29 0
31 0
30 2
28 1
35 2
22 0
23 3
24 1

Order Score
1/162 442.66
1/153 360.03
2/146 229.89
1/138 331.59
1/132 219.45
2/127 217.92
4/117 196.15
2/108 234.27
5/102 302.02
3/99 278.34
2/91 305.22
1/83 408.63
3/79 166.33
2/74
60.67
1/64
60.12
Total time:Time(s)
7.01
14.29
19.12
6.72
18.61
18.76
51.69
29.87
33.77
17.40
13.58
7.76
18.52
26.32
6.57
289.98

(a)

(b)

(c)

(d)

Figure 15: Result of reassembly.
a), b) Reassembled results by our method.
c), d) Reassembled ﬂakes of two groups.

Table 4: Comparison table with previous study

Group 1
Group 2

Our
Previous
Our
Previous

N.Dis (mm)
0.012
0.033
0.016
0.027

Time (s)
400
5912
290
1641

diﬃculties in this scenario. In the experiments,
this case occurred only once.

5 Conclusion and future work
In this study, we introduced a novel matching algorithm for stone tool reassembly based on contour points derived from ﬂake surfaces. Our approach addressed the critical challenges encountered in a previous study [16], resulting in signiﬁcant reductions in computation time and improvements in partial matching.
However, it is essential to acknowledge the limitations of our method, particularly in cases inwhich maximum contour point matching does notyield a correct match. Although these instances
are rare, they should be highlighted in futurestudies.

Acknowledgments
A part of this study was supported by JSPS
KAKENHI (Grant Number JP22K00998).

References[1] K. Matsufuji, S. Monta, Yoku wakaru koukogaku (Understand archaeology), Minerva
Shobo, Kyoto, pp. 18, 2010.
[2] K. Suzuki, Koukogaku Nyuumon (Archaeology Introduction), University of Tokyo
Press, Kyoto, JP, 1988.
[3] E. Altantsetseg, Y. Muraki, F. Chiba, K.
Konno, 3D Surface Reconstruction of Stone
Tools by Using Four-Directional Measurement Machine, International Journal of Virtual Reality, Vol. 10, No. 1, pp. 37-43, 2011.
[4] Q. Huang, S. Flory, N. Gelfand, M. Hofer,
Reassembling Fractured Objects by Geometric Matching, ACM Transactions on Graphics(TOG), vol. 25, Issue 3, pp. 569-578, 2006.

– 4:14 –The Journal of the Society for Art and Science, Vol. 23, No. 2, pp. 4:1 – 4:17 (2024)[10] N. Lamb, C. Palmer, B. Molloy, S. Banerjee,
N. Banerjee, Fantastic Breaks: A Datasetof Paired 3D Scans of Real-World Broken
Objects and Their Complete Counterparts,
2023 IEEE/CVF Conference on Computer
Vision and Pattern Recognition (CVPR),
Vancouver, BC, Canada, pp. 4681-4691,
2023.

Figure 16: Limitation of our method: Failed matching scenario for stone N o.14.

[5] B. Bronn, C. Toler-Franklin, A System
for High-Volume Acquisition and Matching of Fresco Fragments: Reassembling
Theran Wall Paintings, ACM Transactions
on Graphics(TOG), vol. 27, Issue 3, pp. 1-9,
2008.
[6] A. Wills, Stochastic 3D Geometric Models
for Classiﬁcation, Deformation, and Estimation. Ph.D. thesis, Brown Univ. Press., 2004.
[7] J. H. Hong, S. J. Yoo, M. A. Zeeshan, Y.
M. Kim and J. Kim, Structure-from-Sherds:
Incremental 3D Reassembly of Axially Symmetric Pots from Unordered and Mixed
Fragment Collections, 2021 IEEE/CVF International Conference on Computer Vision(ICCV), Montreal, QC, Canada, pp. 54235431, 2021.
[8] R. Gregor, D. Bauer, I. Sipiran, P. Perakis,
T. Schreck, Automatic 3D Object Fracturing for Evaluation of Partial Retrieval andObject Restoration Tasks - Benchmark andApplication to 3D Cultural Heritage Data,
The Eurographics Association, 2015.
[9] N. Lamb, S. Banerjee, N. Banerjee, Automated reconstruction of smoothly joining 3D
printed restorations to ﬁx broken objects, In
Proceedings of the 3rd Annual ACM Symposium on Computational Fabrication (SCF
’19), Association for Computing Machinery,
New York, NY, USA, Article 3, pp. 1 ‒ 12,
2019.

[11] M. Zhang, H. You, P. Kadam, S. Liu, C.
Kuo, PointHop: An Explainable Machine
Learning Method for Point Cloud Classiﬁcation, IEEE Transactions on Multimedia, vol.
22, no. 7, pp. 1744-1755, 2020.
[12] P. Besl, N. McKay, A method for Registration of 3-d Shapes, IEEE Trans Pattern Anal
Mach Intel, vol. 14, no. 2, pp. 239 ‒ 256,
1992.
[13] N. Mellodo, N. Mitra, D. Aiger, Super 4PCS
Fast Global Point cloud Registration via
Smart Indexing, Computer Graphics Forum,
vol. 33, no. 5, pp. 205 ‒ 2015, 2014.
[14] R.B. Rusu, N. Blodow, M. Beetz, Fast Point
Feature Histograms (FPFH) for 3D registration, 2009 IEEE International Conference on
Robotics and Automation, pp. 3212 ‒ 3217,
2009.
[15] M.A. Fischler, R.C. Bolles, Random Sample Consensus: A Paradigm for Model Fitting with Applications to Image Analysis
and Automated Cartography, Readings in
Computer Vision, pp. 726 ‒ 740, 1987.
[16] X. Yang, K. Matsuyama, K. Konno, A New
Method of Reﬁtting Mixture Lithic Materials by Geometric Matching of Flake Surfaces,
The Journal of the Society for Art and Science, Vol. 15, No. 4, pp. 167-176, 2016.
[17] Point Cloud Library, Region growing segmentation.
[18] R. Osade, T. Funkhouser, B. Chazelle, D.
Dobkin, Matching 3D models with shape
distributions, Shape Modeling and Applications, SMI 2001 International Conference,
pp. 154-166, 2001.

– 4:15 –The Journal of the Society for Art and Science, Vol. 23, No. 2, pp. 4:1 – 4:17 (2024)[19] Mark Pauly, Richard Keiser, Markus Gross,Multi-scale Feature Extraction on PointSampled Surfaces, Eurographics, Vol. 22,
No. 3, 2003.

interest include geometric modeling, virtual reality, and simulation.

[20] H. Edelsbrunner, D. Kirkpatrick, R. Seidel,
On The Shape of a Set of Points in The
Plane, Information Theory, IEEE Transactions on Information Theory, vol. 29, no. 4,
pp. 551- 559, 1983.

Mengbo You[21] A. Igarashi, Hakuhen hakuri genri: Seisei no
zengo kankei (Flakes peeling principle: The
context of generation), Sekki zukuri no ziken
koukogaku (Experimental archaeology of theStone tool Making), Lithic Technology Research Society, Gakuseisha, Tokyo, JP, pp.
22-35, 2004.
[22] A. Altansukh, M. You, E. Altantsetseg, O.
Khorloo, K. Konno, A Study on Automatic
Flake Surface Segmentation of Stone Tools
by Calculating Shape Features, IWAIT,
SPIE Digital Library, Vol. 12592, pp. 123128, 2023.
[23] K. Yamahara, K. Konno, F. Chiba, M.
Satoh, A Method of Detecting Adjacent
Flakes in Stone Tool Restoration by Extracting Peeling Surfaces, Japan Society for Archaeological Information., Vol. 17, No. 1-2,
pp.23-31, 2011.

is an assistant professor in the Faculty of Science & Engineering, Iwate University. He received the B.S. degree in 2012 from the Computer
Science Department, Northwest A&F University,
the M.S. degree in 2015 and the Dr.Eng. in 2018,both from Iwate University, Japan. His researchinterests include machine learning, image processing, object detection and deep learning. He is amember of the Society for Art and Science, andthe Institute of Image Electronics Engineers ofJapan.

Enkhbayar Altantsetseg[24] Inter i7-10700 vs. Intel i7-4790, CPU Benchmark. Accessed on: Mar.17, 2024, Available:
https://www.cpubenchmark.net/compare/3747vs2226/Intel-i7-10700-vs-Intel-i7-4790.
received the B.S. and M.S. in mathematics from
National University of Mongolia in 1995 and1997, respectively. He earned his Dr. Eng. indesign and media technology from Iwate University in 2013. He is currently a professor of the
School of Information Technology and Electronics
of National University of Mongolia. His researchinterests include geometric modeling, virtual reality, augmented reality, and simulation.

Amgalan Altansukh

is a Ph.D candidate of Graduate School of Scienceand Engineering at Iwate University. He received
B.S. and M.S degrees from the National University of Mongolia in 2018 and 2019. His research

Oyundolgor Khorloo– 4:16 –The Journal of the Society for Art and Science, Vol. 23, No. 2, pp. 4:1 – 4:17 (2024)

is currently an associate professor in the Department of Information & Computer Science at
School of Information Technology and Electronics
of National University of Mongolia. Her researchinterests include computer graphics, virtual reality, noise-based animation, and simulation of natural phenomena. She received a B.S. degree inapplied mathematics from National University ofMongolia in 1995 and an M.S. degree in computerscience from University of Colorado at Denver in
2001. She received her Dr. Eng. in design andmedia technology from Iwate University in 2012.

is a professor of Faculty of Science and Engineering at Iwate University. He received a BS in Information Science in 1985 from the University ofTsukuba. He earned his Dr.Eng. in precision machinery engineering from the University of Tokyoin 1996. He joined the solid modeling project at
RICOH from 1985 to 1999, and the XVL projectat Lattice Technology in 2000. He worked on
an associate professor of Faculty of Engineering
at Iwate University from 2001 to 2009. He has
written a book Introduction to 3D shape processing. His research interests include 3D modeling,
3D surface data compression, archaeological relics
restoration. He is a member of The Society for
Art and Science, The Institute of Image Information and Television Engineers , Japan Association for Archaeoinformatics, Information Processing Society of Japan, and EuroGraphics.

Fumito Chiba

received the Dr.Eng. degree from the Department of Electronic Engineering and Computer
Science, Graduate School of Engineering, Iwate
University. He worked as a research associate in
the Department of Computer Science, Faculty ofEngineering, Iwate University, and now he is theManaging Director of LANG CO., LTD. He is engaged in research and development of 3D shapemeasurement and processing of archaeological artifacts. He is a member of Japan Association for
Archaeoinformatics.
Kouichi Konno– 4:17 –