DEVELOPMENT OF AN ALGORITHMIC PROCEDURE FOR THE DETECTION OF CONJUGATE FRAGMENTS

: The rapid development of Computer Vision has contributed to the widening of the techniques and methods utilized by archaeologists for the digitization and reconstruction of historic objects by automating the matching of fragments, small or large. This paper proposes a novel method for the detection of conjugate fragments, based mainly on their geometry. Subsequently the application of the Fragmatch algorithm is presented, with an extensive analysis of both of its parts; the global and the partial matching of surfaces. The method proposed is based on the comparison of vectors and surfaces, performed linearly, for simplicity and speed. A series of simulations have been performed in order to test the limits of the algorithm for the noise and the accuracy of scanning, for the number of scan points, as well as for the wear of the surfaces and the diversity of shapes. Problems that have been encountered during the application of these examples are interpreted and ways of dealing with them are being proposed. In addition a practical application is presented to test the algorithm in real conditions. Finally, the key points of this work are being mentioned, followed by an analysis of the advantages and disadvantages of the proposed Fragmatch algorithm along with proposals for future work.


INTRODUCTION
Over the last few years, archaeological science increasingly utilizes the possibilities offered by the use of computers and especially by the innovations implemented in the field of Computer Vision.A significant reduction in processing time required for the reconstruction of ancient monuments with the application of automatic or semi-automatic methods and the 3D digitization of fragments has been witnessed.Within this framework, this paper focuses on exploring an innovative algorithm, suitable for three-dimensional matching of broken objects.Specifically an algorithmic process, which is based mainly on the 3D geometry of the fragments, is being proposed and developed.In this paper the limits, the reliability and the problems of this process are being examined through simulations and a practical application.

RELATED PREVIOUS WORK
Reconstructing broken objects from fragments with computeraided techniques represents a challenge that has been studied for many years (Willis and Cooper 2008).The first approach in this discipline faced the problem of solving jigsaw puzzles (Freeman andGarder 1964, Goldberg et al. 2004).Jigsaw puzzles have a number of important constraints, which aid matching and represent the most basic scenario for automatic reconstruction in 2D space.Contour matching techniques provide solutions to more general problems, without distinguishing between specific edges or features.Applied to 2D space, some of the most outstanding techniques can be found mainly in (Da Gama Leitao and Stolfi, 2002) and in (Papaodysseus, et al., 2002).A very common extension of contour matching techniques in archaeological fragment reconstruction is pottery re-assembly techniques.These are applied to revolution surfaces and can take advantage on the additional constraint of axial symmetry (Kampel and Sablatnig, 2004, Karasik and Smilansky, 2008, Papaodysseus et al., 2002, Ucoluk and Toroslu, 1999).Surface matching techniques are applied to more general problems, considering 3D input data, and provide solutions to problems with six degrees of freedom (3 translations and 3 rotations).The first approach was presented by Papaioannou and Karabassi (2003), Papaioannou et al. (2002) and Papaioannou et al. (2001) with the underlying assumption that the fractured faces were nearly planar and they matched each other completely.Using a projective space, GPU (Graphics Processing Unit) depth maps where analyzed to reconstruct the original object.An optimization of this technique is further presented by Belenguer and Vidal (2012).Stanford's Digital Forma Urbis Romae project (Koller and Levoy 2006) dealt with heavily eroded fragments, whose fractured surfaces sometimes did not even touched each other.In this approach, instead of using the geometry, reconstruction is done by matching incisions on the fragment's top surfaces.
The technique presented by Huang et al. (2006) reassembles solid objects by first identifying fractured regions and then generating clusters of feature patches for alignment-based matching.This approach also introduces an effective technique for multi piece global matching of fragments.Finally, Brown et al. (2010) exploit the orientation constraints of flat fragments to achieve a simple, fast matcher based on edge geometry.The proposed technique analyzes exhaustively every possible alignment of a pair of fragments in a few seconds.

FRAGMATCH ALGORITHM
The idea that spawned the Fragmatch algorithm is that the comparison between the broken surfaces should be accomplished by a linear and quick method.That is the main innovation in relation to all previous investigations.

Data input
The algorithm accepts as input data point clouds, i.e.X, Y, and Z for each point of the broken surfaces and not of the whole fragment.It is assumed that while scanning the broken surface should be located in front of the scanner, so that the Z axis defines the depth of the broken surface.The program asks the user to manually enter three pieces of information, i.e. i) the number of pieces ii) the number of broken surfaces for each piece and iii) the name of the file for each broken surface, which corresponds the to the data point cloud.

Comparison based on ICP algorithm
The use of ICP has been employed in order to accelerate the software application.The most common use of the ICP is to match two groups of data points into a common coordinate system (Gelfand et al., 2003).However in this case it is used as a method of comparing the similarity between the broken surfaces.During the application of ICP, two tables (A and B) are being created.The results from the pair wise comparisons of all the surfaces are being placed in table A. The first two columns of table A contain the compared surfaces; the third column takes values 0 or 1 depending on the success of the ICP matching.The fourth column contains the value of the error.Table B contains the compared surfaces, which have been evaluated with an acceptable error.

Creating surfaces and vectors
The point clouds are converted to surfaces through the process of Delaunay triangulation; this is the only non-linear procedure and is necessary for the implementation of Fragmatch algorithm.The triangles formed for each broken surface are sorted according to their size, while triangles with sides exceeding a manually set limit are excluded.
After the formation of the surfaces, the normal vector of each triangle is determined.The application point is the centre of gravity of each triangular surface.Each vector's norm is equal to the area of the applied triangle.This characteristic ensures the uniqueness of triangles and vectors.

Global Matching
For each broken surface a table K is being created.The dimensions of this table are mx7, m being the number of the created triangles.The first three columns of table K define the coordinates of the triangle's centre of gravity, which is also the normal vector application point.The next three columns define the three vector's components and the last column contains the norm of the vertical vector, which is the area of the triangle.Table K is sorted in descending order of the seventh column.Then the angles between the vectors are being calculated.The calculation of all the angles is being made by the dot product.
Each angle is being formed by the respective vector with the first vector of the table.The angles are compared based on the order they were registered.Then the subtractions between two angles are being calculated.The final success rate depends on the number of subtractions that were below the tolerance limit.The area is calculated according to the seventh column of table K. Two sums are being calculated for each comparison, one sum for every surface.The matching percentage is estimated based on the following formula: If a 1 and a 2 are sums of the triangle areas and a 1 <a 2 then Area (%) = a 1 /a 2 *100 (1) else Area (%) = a 2 /a 1 *100 (2)

Angle Tolerance:
The tolerance of the subtraction between two angles is defined by the following formula: d is the surface's wear and α is the triangle's side.The following figure illustrates two conjugate surfaces; each of them is composed by two triangles.Theoretically the two compared angles are θ 1 and θ 2 , however because of the appearing wear in the second surface, the formed angle, between the face 3 and face 4 is the θ 2 ' (Figure 2).The subtraction between angles θ 1 and θ 2 is the requested angle φ' (Figure 3).

Partial Matching
In partial matching the same methodology is being applied, i.e. comparison of vector angles and areas.By partial matching the correspondence between fragments of different dimensions is implied.As in the case of global matching, for each broken surface a table K is again created.Assuming that the smaller surface is subsurface of a bigger one, common elements are being identified based on the seventh column of table K.The elements, whose area's subtraction is less than a defined threshold, are considered as common.Then all accepted elements are placed in a new table F. Having two matrices of the same length, matching is being performed in the same way as in the case of global matching.In this case except of the angles' and areas' success rate, a percentage of the greater surface covered by the smaller one, is being determined and displayed.K.

Steps of the algorithm
The input data include (i) the number of pieces, (ii) the number of broken surfaces of each piece and (iii) the name of the file for each broken surface, which corresponds the to the point cloud.The output data is the table with possible conjugate couples and their matching percentages both for the area and angles criterion.
Then the algorithm follows the following procedure:

Introduction
The first simulation is going to determine the limits of tolerable noise.The second factor, the accuracy of the scan, is going to be studied with the same methodology as the first factor.Then the next simulation ensures whether the algorithm is affected not only by the number of scanning points, but also by the number of created triangles.One of the most common effects for the surfaces of fragments is corrosion, mainly due to time.Hence, the next simulation studies the performance of the algorithm, handling damaged surfaces.The last simulation analyses whether the algorithm's performance is being affected by the size and shape of the surfaces.The algorithm was created to examine objects whose sides range from 2-3 cm to 1 meter and therefore object dimensions vary between these limits.
Comparisons in all simulations were made between two broken surfaces.

Noise and accuracy of scanning
In this simulation the size of the compared surfaces is shown in the following The algorithm is being performed for three surface pairs.In the first case, the two surfaces are invariant.In the second case the second surface has undergone alterations up to 0.1mm.In the last case the second surface has undergone alterations up to 1mm.The results are presented below: The results of the first simulation reveal that with a fixed number of triangles, the success rate of the algorithm is reduced, as the distortions increase, which is more or less expected.Since the distortion in the z-axis (depth) increases, different triangles are being created in two different surfaces.However the final result in the above cases is that the three pairs constitute possible conjugate fragments.

Number of scan points
The second part of simulations evaluates the results as far as the number of scan points and the created triangles.The algorithm was applied four times in this experiment.In each application, the input data were two point clouds.The first does not have any distortion, while the second has distortion up to 0.1mm in all three axes x, y and z.The studied subjects had the same dimensions so as to be comparable.The results (Table 4) indicate that objects of the same size, as the number of points and triangles increases, the percentage of acceptable angles is being decreased.This implies that the algorithm has a greater tolerance when dealing with the area section, compared to angles section.In

Wear of the surfaces
In these simulations the input data are again two point clouds.
The first point cloud is not supposed to have suffered any damage, while the second had damages introduced in all three axes.Dimensions, surfaces elements and results of the comparisons are presented in the following tables: Both results of these simulations can be characterized as successful.In the second simulation several different triangles have been formed.The direction of triangles and vectors has changed and therefore the percentage of angles decreased so much.

Diversity of shapes
A shattered object generates completely random surfaces.Limits of the algorithm attributed to the shape diversity are detected by this simulation.Four applications of the algorithm were performed, each one with a completely different shape.
Once again the input data were two point clouds; the first one with no distortions and the second one with distortions up to 0.1mm along the direction of all three axes.Results of the comparisons are presented in table 9: In these four applications the percentage of angles is constant and close to 80%, while the percentage of the area is also stable and close to 100%.The algorithm does not depend on the shape of the surface, as the results in different cases are considered similar.

PRACTICAL APPLICATION
In this section, results of a practical implementation of the algorithm are presented.The fragments used for these tests come from the Laboratory of Geotechnics and Rock Mechanics and are actually cylindrical pieces from test drillings through various rock layers.These specimens have then been tested for shear in the laboratory until they fracture.It is these fractured, eroded and shattered surfaces that were put to test by the algorithm.
Scans were performed with the "Next Engine" 3D laser scanner.Three fragments were scanned: a pair of conjugate fragments 1 & 2 (Figure 10) and one single broken piece, i.e. fragment 3 Figure 11).A point cloud of the "conjugate" piece of fragment 3 (fragment 4) was created using Matlab.Fragment Fragment Area(%) Angles(%) The results of Table 10 are considered successful only for simulated conjugate fragments 3 and 4. On the contrary in cases of physical fragments (1 & 2), the number of triangles created for the comparisons, proves to be very high for the proper operation of the algorithm, as concluded during the simulations.
For the second test, the scanned points were decreased to approximately 6000 using Geomagic ® software.Despite this decimation the areas of interest maintained the essential level of detail.The average area of triangles increased to 0.5 mm 2 .Moreover angle φ tolerance could be reduced to 13 grad.
Results from test 2 are being presented in Table 11: Fragment Fragment Area(%) Angles(%) From the successful results in Table 11 of the practical application using really dubious conjugate surfaces prove that even under adverse conditions the algorithm performs well and that the only critical parameter is the number of 3D points used to describe the surface and form the triangles.

Advantages and disadvantages
Summarizing the most important advantages and disadvantages of Fragmatch algorithm are being presented.First of all, the algorithm performs a linear comparison between the elements that ensures fast comparisons.Using Fragmatch algorithm it is not necessary to use any another software to create triangles and surfaces.The strong point of the algorithm is the use of the normal vectors whose norm is the actual area of the triangle.The algorithm performs fast surface comparison regardless of shape or form.The necessary input data are easy to acquire and are considered standard nowadays, hence the required time for that is very limited.Last but not least, not only global, but also partial matching is being performed.
On the other hand, Fragmatch algorithm identifies correspondences between conjugate fragments, but it does not unite the broken pieces.Another disadvantage of the software is that receives as input data only scanned fractured surfaces and not the entire scanned piece.Algorithm's sensitivity is in the number of scan points; the user should be very careful in choosing the right scan density.However, this density may always be decreased later.Finally the size of the input data is limited due to the limitations of the software Matlab.

Proposals
Fragmatch algorithm has been created having as primary goals, the linear comparison and the rapid implementation.The algorithm can be extended based on the following proposals.Determine a correlation coefficient to deal with the multitude of points with the angle section.That factor could eliminate the sensitivity at many scan points.An extension could be made so that broken surfaces could be found automatically by the software.In addition another goal is the unification of pieces and the update of the entire object.Furthermore the transcription of the algorithm in a high level programming language would eliminate the problem with limited size input data.Matching of the pieces could be made more easily by the identification of characteristic forms or by exploiting information about the shape of the final product.Finally an automatic virtual completion of missing pieces could also be added.
Figures 4 & 5 the surfaces are visualized and in Figures 6 & 7 the two graphs show the variation of the matching percentages with the variation of the number of scan points.

Figure 4 :
Figure 4: The first pair of surfaces, with 50 points and 78 triangles each.

Figure 5 :
Figure 5: The third pair of surfaces, with 5000 points and 9055 triangles each.

Figure 6 :
Figure 6: The percentages of the angles in the four phases of the simulation.

Figure 7 :
Figure 7: The percentages of the area during the four phases of the simulation.

Figure 11 .
Figure 11.Fragment 3Two tests were performed.For the first implementation of the algorithm, each cloud consisted of approximately 21000 points.At this high point density, small triangles with an average area of 0.14 mm 2 , were created.Consequently the tolerance angle (φ) is set at 50 grad.The comparison of both area and angles for fragments 3 and 4 should derive results that reach 100%, since they are mathematically created as conjugate.The next higher results of compatibility are expected from the comparison of fragments 1 and 2, as they are physically conjugate fragments.The results of the test are presented below: 1.It compares point clouds with ICP algorithm and creates tables A and B. 2. It creates surfaces and vectors for each fragment.3. It performs the comparison of the pieces.If table B is empty, the pieces are being compared for partial matching.Tables A1 and B1 are being created with the final results.If table B is not empty, its data are being compared for global matching and a table E1 with the final results is being created.Remaining elements are compared for partial matching and matrices A1 and B1 with the final results are being created.

Table 1 .
Size of the surfaces

Table 2 .
Characteristics of the Surfaces.

Table 3 .
Results of the three cases

Table 4 .
Results of the four applications of the algorithm

Table 5 .
Dimensions of the simulation surfaces

Table 6 .
Characteristics of the simulation surfaces

Table 7 .
Results from the first pair of compared surfaces

Table 8 .
Results from the second pair of compared surfaces

Table 9 .
Results of the four applications of the algorithm

Table 10 .
Results from the first test (conjugate pairs are in the shades lines)

Table 11 .
Results from the second test (conjugated pairs are in the shaded lines)