Application of computer vision techniques for 3D matching and retrieval of archaeological objects

Basic requirements

A system of shape recognition must be able to discover shape similarities between partially isometric objects, that is, between objects that share shape characteristics even if they are not identical. For example, a researcher might need 3D models of all the anthropomorphic figurines in a museum. In this case, the query should not be affected by the fact that one object lacks a head, while another is missing a leg. Also, it must retrieve several complete models of the same class, regardless of whether they differ in certain details (Gal & Cohen-Or, 2006). The three objects shown in Figure 1 illustrate this situation; they belong to the same class of anthropomorphic figures, but their heads and arms show some morphological differences.

Figure 1. Three figurines found in the Sacred Precinct of Tenochtitlan.

They are similar in their overall geometry, but differ in some details, such as the design of the head and the position of the arms. A recognition system must be able to find the 3D models of these objects, despite their partial differences in shape.

The second requirement is that the system be able to detect similarities without being affected by affine variations such as rotation, translation, reflection and scale of the 3D models. Ideally, different combinations of translation, rotation, and scale applied to equal objects during the digitisation process should not affect the system’s capacity to recognise their similarity.

To fulfil those requirements, it is necessary to apply a specialized set of methods of computer vision specifically designed to identify objects’ similarities efficiently, within the processing limits of today’s computers and that are sufficiently discriminatory to resolve the requirements of cultural heritage institutions.

Development of a search-engine module

The development of shape recognition systems has been the subject of research since the 1980s, especially in the fields of computer vision, geometric modelling, and machine learning (Besl & Jain, 1985; Bustos et al., 2004, 2005; Campbell & Flynn 2001; Jain & Mishra 2014; Lara López et al., 2017; Loncaric 1998; Tangelder & Veltkamp 2008; Theologou et al., 2014; Veltkamp & Hagedoorn, 1999). In the early 2000’s. Funkhouser et al. (2003) developed a system for retrieving 3D models applying the shape descriptor of spherical harmonics, implementing an interface that allowed queries by example and also sketch-based searches. Another early system was developed by Paquet and Rioux (2000) with shape descriptors based on tensors of inertia, distribution of normal vectors, distribution of cords and multiresolution analysis. Its interface allows entering a combination of parameters such as scale, shape or colour to search the database. Other innovations were due to Suzuki (2001) who included the processing of material data (colour and texture) as search criteria; although its descriptors require normalization of 3D models via PCA.

In the field of cultural heritage, Rowe et al. (2001) and Rowe and Razdan (2002) implemented a shape-based search engine for analysis and retrieval of native American ceramic vessels. Objects were modelled as parametric surfaces and the interface allows query by example and sketch-based query, and links to descriptive data. Another interesting system was designed by Schurmans et al. (2002) according to specific archaeological research objectives. For example, indicators of craft specialization can be gathered from the morphology of ceramic vessels. This involves matching shapes, as well as text, numeric, and vessel data calculated with the system tools.

Those projects demonstrate that the effectiveness of a recognition system depends above all on implementing efficiently two basic procedures. The first one consists in calculating a “shape descriptor”, that is a numerical representation of the form of each 3D model. Such descriptor can represent the global geometry of the object or a sample of its local features. The computation of shape descriptors involves a combination of mathematical, statistical, and more recently Machine Learning methods to represent shape in a numerical array or feature vector (Tangelder & Veltkamp, 2008). Some examples of characteristics encoded by shape descriptors are the curvature or orientation of a certain quantity of patches drawn around points chosen in a random manner (Jiménez-Badillo et al., 2013), or alternately, signals calculated with spherical functions (Kazhdan & Funkhouser, 2002; Kazhdan et al. 2003), reflective symmetry (Kazhdan et al., 2002, 2004a, 2004b), spin-images (Johnson and Hebert, 1999), shape-contexts (Mori et al., 2001), histograms of spherical orientations (Roman-Rangel et al., 2016), and many others as described in bibliographic surveys by Bustos et al. (2004, 2005), Lara López et al. (2017), and Rostami et al. (2019).

More recently, deep-learning models have been developed to compute shape descriptors for the recognition of generic objects. The spectrum of techniques includes supervised and self-supervised approaches (Herrewegen et al., 2023), such as Neural Networks (Krestenitis et al., 2020; Xie et al., 2017; Yang et al., 2021), Convolutional Neural Networks (Feng et al., 2018; Kim & Chae, 2022), Autoencoders (Furuya & Ohbuchi, 2020), Point Cloud Networks (Wang et al., 2020), Graph Neural Networks (He et al., 2020), and Geometric Deep Learning techniques (Bronstein et al., 2017). These methods, however, require large data sets for training the models, a condition rarely met by cultural heritage collections, including the one used for this project.

In any case, the final objective is that the shape of the object is characterized in the best possible manner, as to constitute a “signature” (numerical representation) of the object readable by a computer. As mentioned above, the numerical descriptor must represent the shape regardless the object’s position, orientation and scale.

The second procedure consists in creating an index of the numerical representations of all the objects (i.e. their shape descriptors), to facilitate the matching operation. The comparison between 3D models is done by measuring their degree of similarity with a mathematical function, such as Euclidean distance, that indicates the degree of their resemblance, so that when a query is implemented the system can rank objects from the more to the less similar (Bustos et al., 2004; Tangelder & Veltkamp, 2008).

The search-engine module developed over the course of this project is based on the implementation of three different global descriptors, namely shape distributions (Osada et al., 2002), reflective symmetry (Kazhdan et al., 2002, 2004a, 2004b) and spherical harmonic functions (Kazhdan and Funkhouser, 2002; Kazhdan et al., 2003).

As for determining the degree of dissimilarity between objects, four measures have been implemented: Euclidean distance, City block (Manhattan) distance, Chebychev distance,¹ and Minimum Coordinate distance (Table 1). In the case of shape distributions, each dissimilarity measure has been implemented in two norms: the probability density function (pdf), and the cumulative distribution function.

The following sections describe, in layman’s terms, the methods to compute the three shape descriptors selected for the implementation of the search-engine module.

Shape distributions

The simplest shape descriptors included in the search-engine module are five probability distributions proposed by Osada et al. (2002). The names given to the descriptors depend on the type of calculation, “A” stands for angle and “D” for distance:

As Osada et al. (2002) recommend, we compute those variables for a very large sample of points selected from the surface-mesh of each 3D model, specifically 1,048,576 points (i.e. 1024 × 1024 points). The measurements were then transformed into a frequency histogram (probability distribution), which could then be used as the global signature of the object’s shape. Once the shape histograms for all the objects had been computed, a normalization step was necessary to standardize the scales of all the histograms in order to avoid matching errors due to variations in the size of the objects. The objective is finding the scale that produces the minimal dissimilarity measure during the comparison of two object’s histograms. To achieve this, one of the methods proposed by Osada et al. (2002) involves the following steps: align both shape distributions (i.e. histograms) so that the mean sample in each distribution equals 1; then find the minimum value D(f (x), sg (sx)) for values of log s from -10, 10, in 100 equally spaced intervals; where f and g represent the shape distributions of two models; x correspond to the 3D model; and D corresponds to the distance function chose from A3, D1, D2, D3, or D4. Finally, select the minimum value among the results and use it as the dissimilarity measure for the two normalised shape distributions. This guarantees that two objects of the same shape but different sizes are recognized as similar, and vice versa, two objects of the same size but different shape are recognised as different.

The resemblance between any pair of objects can be determined by applying a function that measures dissimilarity between distributions (i.e. histograms), for example Euclidean distance or any of the other measures mentioned in Table 1.

Figure 2 shows histograms resulting from the descriptor A3 (angle between three random points), representing the shape of four archaeological objects. Notice the probability distributions of the two models on the left, reflecting the differences between the long, wavy form of the serpentiform sceptre and the flat, wide anthropomorphic figurine. For an elongated figure, the angles between vectors tend to concentrate around the mode, while for flat objects the histogram would have a flat distribution because there would not be a predominant value of angle between vectors. In contrast, the images on the right correspond to two vessels whose histograms are quite similar because their shapes are also alike. Through this kind of comparison, the recognition system manages to identify similarities or differences between archaeological objects.

Figure 2. Frequency probability histograms resulting from obtaining measurements of angles between vectors for four objects from the Templo Mayor museum collection.

Notice the great difference between the histogram of the serpentiform scepter (first object on the left), which shows the mode close to zero, and the histogram of the flat anthropomorphic figure (second from left to right). On the other hand, the great similarity of the histograms corresponding to the so-called Tlaloc vessels located on the right side can be appreciated. Comparing these histograms allows the recognition system to determine whether or not two objects belong to the same class.

Reflective symmetry descriptor

The second representation of shape, more complex but at the same time more effective, is the reflective symmetry descriptor proposed by Kazhdan et al. (2002, 2004a, 2004b). As these authors point out, symmetry —or the lack of it — is one of the most distinctive characteristics of any object.

Given a 3D model, denoted by function g, the concept of reflective symmetry implies that there is a reflection function γ, such that $g=\gamma \left(g\right)$. This means that the pointwise distance between the points of surface g and the points of surface $\gamma \left(g\right)$ is zero.

Kazhdan et al. (2002) propose quantifying reflective symmetry with respect to several cutting planes, oriented on perpendicular axes that pass through the centroid of the model. For any given plane P cutting the shape f, the method consists in finding the function g such that $g=\gamma$ with $⟦f-g⟧$ as small as possible. Mathematically this is expressed as:

where SD stands for Symmetry Descriptor. The more symmetric the shape f with respect to plane P, the smaller the value of $‖f-g‖$. Large values of $‖f-g‖$ indicate that the surface is less symmetric.

The calculation of reflective symmetry can be performed quicker and more efficiently by transforming the description of the surface mesh (3D model) into a discrete volumetric representation (i.e. voxel grid). The process starts by immersing the triangular surface mesh inside a regular 3D grid. When a triangle of the mesh intersects a voxel of the 3D grid, such voxel is assigned a value of 1. Such rasterization process allows determining where the points and triangles of the mesh are located on the 3D grid. Notice that there are voxels in the grid that do not intersect the 3D mesh and therefore lack any information. For calculating the reflective symmetry descriptor, it is convenient to add to these voxels information related to how far they are from the surface of the model, for which the distance transform is used. This transform consists of assigning each voxel the distance to the nearest voxel belonging to the model. Then, the distance is transformed to a measure of similarity with the Gaussian function. Additional voxelization methods available can be found in Aleksandrov et al. (2021) and Huang et al. (1998). The resulting discrete representation consists of a 3D set of voxels, which appears like the archaeological model shown in Figure 3.²

Figure 3. A 3D model of an archaeological figure and its voxel representation.

Once the voxel grid has been labeled in this way, the descriptor can be calculated. Broadly speaking, what it is done is to assume that when passing a cutting plane through the voxel grid there is perfect reflective symmetry between the two halves, so that one of the two halves can be replaced with the other if that property were fulfilled. Then, the reflective symmetry distance is calculated in the real model and compare it with the assumed model to measure any difference. If both representations are equal (zero distance), then there is perfect symmetry and a radius of 1 is assigned to the corresponding plane. If not, a value less than 1 is assigned according to how different these figures are.

Finally, the measures of symmetry obtained from a number of cutting planes (i.e. axes of symmetry) are concatenated to generate a 3D graph, describing the model’s global symmetries. Values near to 1 indicate perfect symmetry, while those near zero indicate that the two halves of a model are too asymmetrical. This graph is used to compare an object with any other for the matching and retrieval application. Visual representations of the descriptors obtained from three archaeological models are shown in Figure 4.

Figure 4. Graphs of the reflexive symmetry descriptors obtained for three objects from the Templo Mayor collection.

As with the histograms shown in Figure 2, reflexive symmetry descriptors offer another method for calculating dissimilarity between 3D models.

Spherical harmonics descriptor

A third way of describing shapes on a computer is to consider them as outcomes of mathematical functions. Each stroke of a drawing, for example, can be regarded as a mix of 2D functions. The numerical representation of the whole drawing would be the sum of many functions. In practice, the functions are unknown, but they can be calculated by applying standard mathematical procedures such as the Fourier Transform, which would find the specific mix of simple functions that represent the complete drawing.

Something similar happens with 3D objects, but in that case the function describing the shape is defined on the surface of the sphere. One way of describing the shape of a surface is calculating the so-called spherical harmonic functions (Figure 5). Intuitively, we can think of the spherical harmonic functions as “Lego” pieces that, together, help built the shape of complex 3D objects. This is possible, because mathematically speaking, spherical harmonics constitute a complete set of orthogonal functions and therefore form an orthonormal basis, upon which any function defined on the sphere (like the shape of a 3D model) can be expressed as the sum of these spherical harmonics.

Figure 5. Graphical representation of harmonic functions calculated to obtain shape descriptors from 3D models.

The exact combination of spherical harmonic functions needed to describe a particular object can be found though a harmonic analysis. This method divides the complex surface of a 3D model into sums of relatively simple components.

Harmonic analysis is the branch of mathematics that deals with the problem of representing functions as the combination of basic elements called waves or harmonics. Typically, the term harmonic refers to functions with sinusoidal variations, but more strictly, it indicates any solution of Laplace’s Equation. The Fourier series is an example of a complete set of harmonics, which consists of sine and cosine waves of different frequencies.

In this work, we adopted the Spherical Harmonic Transform (SHT) to obtain a reduced representation of a 3D mesh. This method is a powerful tool for describing data on a sphere using spherical harmonics as basic functions. Given a function $f\left(\theta ,\phi \right)$ in the spherical coordinates $\theta$ and $\phi$ the decomposition of $f\left(\theta ,\phi \right)$ in spherical harmonics $Y_l^m\left(\theta ,\phi \right)$ is written as:

Here, l ≥ 0 and m are integers such that $\left|m\right|\le l$, $c_{\mathit{lm}}$ is the coefficient of the harmonic $Y_l^m\left(\theta ,\phi \right)$, and the general form of $Y_l^m\left(\theta ,\phi \right)$ is:

where $P_l^m\left(x\right)$ is a Legendre polynomial:

The problem in the Spherical Harmonic Transform is to calculate the coefficients $c_{\mathit{lm}}$.

In practice, it is not possible to calculate the coefficients of all the spherical harmonic functions. For this reason, we limit the order of the harmonics to a fixed value b (for instance 16 or 32) so that:

Finally, the coefficients $c_{\mathit{lm}}$ are estimated by finding the least-squares solution to equation (5). That is to say, for a set of n points $\left\{\left({\theta }_1,{\phi }_1\right),\left({\theta }_2,{\phi }_2\right),\dots ,\left({\theta }_n,{\phi }_n\right)\right\}$ where $f\left(\theta ,\phi \right)$ is evaluated, we calculate the values of the coefficients $c_{\mathit{lm}}$ that minimize:

To describe a 3D mesh using the Spherical Harmonic Transform, we define the function $f^r\left(\theta ,\phi \right)$ as the intersection between the voxelized version of the 3D mesh and the sphere of radix r, both centered at the origin. The function $f^r\left(\theta ,\phi \right)$ takes the value 1 only if the sphere intersects a voxel of the mesh at the point $(\theta ,\phi$), otherwise, this function is 0. The Spherical Harmonic Transform is applied to different radii so that the functions $f^r\left(\theta ,\phi \right)$ are characterized by their corresponding harmonic coefficients.

The simplest harmonic function is the sphere, so if the object resembles a balloon, only one harmonic component of degree zero is enough for describing it. However, if the model has a more complex shape, then it is essential to calculate several higher order harmonic functions.

There are several methods to compute shape descriptors based on spherical harmonics. Some require a priori registration of the model along principal axes (Saupe and Vranic 2001; Vranic & Saupe, 2001; Vranic et al., 2001), but these are not good to process 3D models of the same class digitised with different orientations (Funkhouser et al., 2003). A method that solves that limitation is the one proposed by Kazhdan and Funkhouser (2002), and Kazhdan et al. (2003) and it is the one implemented during this project. In practice, the descriptor is computed as follows:

Figure 6. Graphs of the decomposition coefficients obtained by calculating harmonic functions for three objects of the collection of the Templo Mayor collection.

To compare two objects using their harmonic representations, it is simply necessary to compute the Euclidean distance between them: “Thus, finding the K closest models to a query is equivalent to solving the nearest-neighbour problem” (Kazhdan & Funkhouser, 2002). Figure 6 shows the spherical harmonics descriptor for three objects.

User interface

Some early systems were tested with collections of 3D models produced with computer-aided design (parametric models) software and acquired on the internet. In contrast, we have developed a first module (MeshAnalyser) of a search engine called ArcheoShape that uses real archaeological objects. At this stage, the objective is to assess how good are the shape descriptors described in the previous sections to match and retrieve real archaeological objects before continuing the development of the entire system and deploying it within a museum environment. The module developed over the course of this project is freely available in the following GitHub repository (Mendoza-Montoya, 2023b), which contains the source code, written in C++, ready to be compiled in Windows and Linux, as well as an executable file. Instructions to compile are included in the GitHub repository. A sample of ten 3D models of archaeological artefacts are also provided under the Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International. These resources will allow any user to test the MeshAnalyser module in a local computer, as well as replicate or customized the module for his own purposes.

To test the implementation of the shape descriptors, we used a sample of nearly 500 archaeological artefacts from the Museo del Templo Mayor. The collection is available for research purposes through specific agreements with Instituto Nacional de Antropología e Historia.³

The Museo del Templo Mayor preserves objects discovered between 1978 and 1982 within the area occupied by the Sacred Precinct of Tenochtitlan, the most important religious centre of the Aztecs and nowadays a famous archaeological site adjacent to Zócalo square in Mexico City (Matos Moctezuma, 1988). The core collection includes more than 8000 objects from ritual offerings found in the main pyramid temple (i.e. Templo Mayor) of the site and its surroundings, and include ritual artefacts, flora, fauna, and human remains (López Luján, 1994; Nagao, 1985). The collection has increased considerably in recent years thanks to the excavations conducted in the same site by different research teams led by archaeologists Barrera Rivera & Islas Domínguez (2018), Barrera Rodríguez, and López Luján (2012, ), 2019. Indeed, between 2012 and 2019, 43 new offerings, containing 13,925 artefacts and 35,648 samples of organic material have been reported. Digitization of this collection is still at a very early stage, but we have been able to acquire a sample of 495 objects, including stone-masks, anthropomorphic and zoomorphic figures, clay vessels (bowls, pots, jars, braziers), religious paraphernalia like sceptres, earplugs, ritual pendants, as well as flint sacrificial knifes, flutes, and models of drums made of clay or stone.

The implementation of the prototype was divided into two independent offline and online jobs. The main offline task consists in computing the shape distributions, reflective symmetry and spherical harmonics descriptors of all the archaeological 3D models available (“Database of descriptors” in Figure 7). The actual matching and retrieval operations are done online.

Figure 7. Architecture of MeshAnalyser.

The module performs four basic operations as follows:

The system’s architecture is crafted to reflect the sequential execution of these operations (Figure 7). As mentioned above, we have made an executable file of the module available through the link provided below in the section Software Availability, as well as the source code. Every section of the source code includes detailed comments that will be updated to reflect changes in its functionality. This would facilitate the replication of the software by others as the software evolves.

Next, we explain how these operations interact when a new 3D mesh is introduced into the system.

Figure 8. Screenshot of the window to open a file in MeshAnalizer.

Figure 9. Screenshot of the window that displays the properties of a 3D mesh opened with MeshAnalizer.

Figure 10. Rendering a triangular mesh as a voxelised model.

Figure 11. Screenshot of the window to apply the Shape Distribution algorithm to perform a matching and retrieval operation with MeshAnalizer.

Figure 12. Screenshot of the window to apply the Relective Simmetry algorithm to perform a matching and retrieval operation with MeshAnalizer.

Figure 13. Screenshot of the window to apply the Spherical Harmonics algorithm to perform a matching and retrieval operation with MeshAnalizer.

Figure 14. Illustration of the software developed for the search and recovery of 3D models of archaeological objects.

Above, on the left, the consultation model (i.e. anthropomorphic sculpture) is shown; below left illustrates obtaining the reflexive symmetry descriptor for that query model; on the right are the search results obtained when the user requests to compare the query model with the models stored in the repository. It can be seen that the system retrieves all objects similar to the query model.

Figure 15. Another example of searching and retrieving 3D models.

In this case, all copies of anthropomorphic figures were requested from the system applying the descriptor of harmonic functions.

Figure 16. A third example of how the software works.

In this case stone masks were recovered. It should be noted that, despite the fact that the query model lacks a fragment, the system was able to produce the expected results, even recovering a mask fragment that clearly belongs to the class of the objects that the user expected.