Complementarity-Preserving Fracture Morphology for Archaeological Fragments

TL;DR

This paper introduces a morphological scale space approach to simplify archaeological fracture surfaces, preserving their complementarity and robustness against abrasion, aiding in fragment fitting analysis.

Contribution

It presents a novel hierarchical simplification method for fracture surfaces using Lipschitz properties and a new embedding for morphological operations.

Findings

Preserves complementarity of fracture surfaces.

Insensitive to abrasion effects.

Effective hierarchical simplification pipeline.

Abstract

We propose to employ scale spaces of mathematical morphology to hierarchically simplify fracture surfaces of complementarity fitting archaeological fragments. This representation preserves complementarity and is insensitive to different kinds of abrasion affecting the exact fitting of the original fragments. We present a pipeline for morphologically simplifying fracture surfaces, based on their Lipschitz nature; its core is a new embedding of fracture surfaces to simultaneously compute both closing and opening morphological operations, using distance transforms.

Figures

Equations

$$⟺$$

$${\gamma }_R(X_M)\subseteq {\gamma }_{\rho }(X_M)\subseteq X_M\subseteq {\varphi }_{\rho }(X_M)\subseteq {\varphi }_R(X_M)$$

$$\Updownarrow c\Updownarrow c\Updownarrow c\Updownarrow c\Updownarrow c$$

$${\varphi }_R(Y_M)\supseteq {\varphi }_{\rho }(Y_M)\supseteq Y_M\supseteq {\gamma }_{\rho }(Y_M)\supseteq {\gamma }_R(Y_M)$$

$$\overline{\partial }A=\partial \overline{A}$$

$$\overline{c}\overline{\partial }A\equiv \overline{\partial }cA(=\overline{\partial }{A}^c),$$

$$\overline{\delta }\overline{\partial }A\equiv \overline{\partial }\delta A\text{and}\overline{\epsilon }\overline{\partial }A\equiv \overline{\partial }\epsilon A.$$

$$\overline{\epsilon }\overline{\partial }A=\overline{\partial }\epsilon A=\overline{\partial }c\delta cA=\overline{c}\overline{\partial }\delta cA=\overline{c}\overline{\delta }\overline{\partial }cA=\overline{c}\overline{\delta }\overline{c}\overline{\partial }A.$$

$$\overline{\epsilon }=\overline{c}\overline{\delta }\overline{c}\text{and}\overline{\delta }=\overline{c}\overline{\epsilon }\overline{c}.$$

$$\overline{\varphi }\overline{\partial }A\equiv \overline{\epsilon }\overline{\delta }\overline{\partial }A\text{and}\overline{\gamma }\overline{\partial }A\equiv \overline{\delta }\overline{\epsilon }\overline{\partial }A,$$

$$\overline{\gamma }\overline{\partial }A=\overline{\delta }\overline{\epsilon }\overline{\partial }A=\overline{c}\overline{\epsilon }\overline{c}\overline{\epsilon }\overline{\partial }A=\overline{c}\overline{\epsilon }\overline{\delta }\overline{c}\overline{\partial }A=\overline{c}\overline{\varphi }\overline{c}\overline{\partial }A.$$

$$\overline{\gamma }=\overline{c}\overline{\varphi }\overline{c}\text{and}\overline{\varphi }=\overline{c}\overline{\gamma }\overline{c}.$$

$${\overline{\gamma }}_R\overline{\partial }(X_M)\subseteq {\overline{\gamma }}_{\rho }\overline{\partial }(X_M)\subseteq \overline{\partial }(X_M)\subseteq {\overline{\varphi }}_{\rho }\overline{\partial }(X_M)\subseteq {\overline{\varphi }}_R\overline{\partial }(X_M)$$

$$\Updownarrow \overline{c}\Updownarrow \overline{c}\Updownarrow \overline{c}\Updownarrow \overline{c}\Updownarrow \overline{c}$$

$${\overline{\varphi }}_R\overline{\partial }(Y_M)\supseteq {\overline{\varphi }}_{\rho }\overline{\partial }(Y_M)\supseteq \overline{\partial }(Y_M)\supseteq {\overline{\gamma }}_{\rho }\overline{\partial }(Y_M)\supseteq {\overline{\gamma }}_R\overline{\partial }(Y_M)$$

Full text

1 Introduction

The GRAVITATE H2020 project aims at providing archaeologists with the virtual tools to analyse digital artefacts, distributed across several collections (https://gravitate-project.eu/). The artefacts are digitally scanned by standard scanning techniques in the form of 3D meshes, capturing the geometrical properties of each object and some of its photometric proprieties. Reassembly of broken artefacts is one of the core objectives of the project.

The test case is on terracotta, which does not deform; but typical fragments are lacking material through abrasion and chipping. Each fragment undergoes a preliminary preprocessing step [5] which partitions its surface into significant sub-parts called ‘facets’. Each facet is characterized by its own geometrical properties of roughness of its surface and sharpness of its boundary, which in turn guides its categorization as either belonging to the fracture region or outside skin region of the fragment.

Our contribution to the project focuses on structuring the pair-wise alignment of promising fragments nominated for fitting by other selection modules. The computational approach is based purely on the geometrical properties of the fracture facets, ignoring other clues (such as possible continuity of decorative patterns).

Optimal pair-wise alignment of potentially counter-fitting fragments is computationally expensive. We therefore take a hierarchical approach, considering the fragments in increasing resolution from simplified to detailed.

The common way of simplifying shapes represented as meshes is by various forms of linear filtering of their vertices; letting them move with a differentiable flow (such as based on the Poisson equation) to produce a smoother model [11]. Structure-aware mesh decimation to preserve the structure of the mesh by representing it as a set of pre-computed canonical proxies has been also studied by [16]. Such simplifications are strongly related to characterizing a surface by differential geometric features, such as those computed by some means of discrete differential geometry [12]. Indeed most methods to reassemble broken objects [9, 19, 15, 7, 14] are all based on such features. However, a small chip, resulting in a missing part of an object, will affect such differential measures in a linear fashion: the deeper the missing bit, the stronger the simplified shape is affected. This is undesirable: the actual match between part and counterpart is simply lacking some local evidence of a more binary nature. The effect of this gap should not depend on the depth of the hole (a gap with the same surface outline but twice as deep is not twice as bad); so our representation method should not be too sensitive to this either. For the same reason, applying ICP (Iterative Closest Point) algorithm to find a matching counterpart based on squared errors is not appropriate; some extensions of ICP that include distance transforms [6] are more in line of what we might employ.

In our hierarchical approach, we need a representation where a coarser, simplified version of the fracture surface maintains its potential fit with a coarser, simplified version of the counterpart, under abrasion-type noise. Mathematical morphology with its opening and closing scale spaces, of increasing size, is the natural choice for just that kind of representation. Such scale spaces have been studied in [10, 1, 2], and it has been demonstrated that they simplify the shape, when measured in terms of certain descriptive features such as local maxima [10], or zero-crossings of the curvature [4]. But our intended use for alignment and complementarity checking appears to be new.

In this paper, we therefore hone standard Mathematical Morphology (MM) to this novel application. First, we define scale-based complementarity of objects (rather than one object and its background). Then, we define morphological operations working only on the fracture region of the fragments. We implement this ‘Boundary Morphology’ by taking advantage of the Lipschitz nature of the individual fracture facets. Our extrusion method uses a single distance transform to simultaneously produce both scaled dilations and scaled erosions of a fracture facet (as in Figure 1), leading to a simplified representation at every scale. We briefly discuss how the fracture Lipschitz property also naturally gives quantitative bounds on the detectable complementarity and collision-free alignment at each scale of our representation, even in the presence of abrasion.

2 Complementarity Preserving Scale Space

Consider a fracture between two fragments $X$ and $Y$ in perfect alignment (or, if you prefer, consider a newly developed crack in an object). Within a well-chosen mask $M$ (such as the part of space the object occupies), the two fragments $X$ and $Y$ are almost complementary in the usual MM sense, with the fracture area being ambiguous (does it belong to $X$ or $Y$?). In an implementation such an infinitely thin layer is not going to make a difference, so we arbitrarily choose $X$ to be a closed set and $Y$ to be open. Then the $M$-restricted parts of $X$ and $Y$, $X_{M}(=X\cap M)$ and $Y_{M}(=Y\cap M)$ are two relatively complementary sets $(X_{M}^{c}=Y_{M})$. The properties of exact complementarity (no overlap, no gaps) can then be algebraically formulated as:

[TABLE]

In archaeology, we do not acquire the fragments to be in alignment, or in perfect condition. We choose MM scale spaces to process the fracture region of each of the fragments separately, producing a simplified representation for each fracture facet at a range of scales $\rho$. Morphological scale spaces are well suited to this, since they can simplify objects while closely maintaining local complementarity, even when objects are slightly damaged (as we will discuss in Section 4).

Performing erosions or dilations would change the objects considerably, so we prefer to use openings and closings to process the fragments (even though much essential structure is already contained in eroded and dilated versions). Because of the arbitrary orientation of presented objects, and an assumed isotropy of fracturing, our morphological scale space is built using structuring elements that are balls of radius $\rho$.

At a scale $\rho$, the opened version $\gamma_{\rho}(X)$ will be complementary (in the above sense) to the closed version $\phi_{\rho}(Y)$, and vice versa. Since a closing $\phi_{\rho}(X)$ simplifies the local geometry of the fracture at valleys, the complementarity with $\gamma_{\rho}(Y)$ becomes less specific in those areas (one can easily construct non-complementary counterparts $Y^{\prime}$ with the same opening $\gamma_{\rho}(Y^{\prime})=\gamma_{\rho}(Y)$, that do not fit the original $X$). To maintain specificity of peaks and valleys of the common fracture surface, we must therefore compute both $\gamma_{\rho}(X)$ and $\phi_{\rho}(X)$ for each fragment $X$. Since both opening and closing are increasing, coarser levels of the scale space are guaranteed to contain less detail, thus enabling the hierarchical approach (as well as being a justification for calling the process a ‘simplification’).

In a schema, if $R>\rho$, we have the following complementarity and containment relationships:

[TABLE]

Here $c$ refers to complementarity within a well-chosen mask $M$ containing the common fracture of $X$ and $Y$. This ‘well-chosen’ actually hides some essential details, since masking does not commute with morphology (e.g., $\phi_{\rho}(X_{M})\neq(\phi_{\rho}(X))_{M})$. It is not hard to show that when originally there was exact complementarity between $X$ and $Y$ within a mask $M$, then in the scale space exact complementarity still holds between $\phi_{\rho}(X)$ and $\gamma_{\rho}(Y)$, or between $\gamma_{\rho}(X)$ and $\phi_{\rho}(Y)$, but is only guaranteed within a doubly eroded mask $\varepsilon_{2\rho}(M)$.

3 Fracture Morphology

Consider a broken vase producing thin sherds as fragments. Volumetric morphological simplification of the fragments would be limited in scale by their thickness: an opening by a ball with a radius $\rho$ of more than half this thickness would results in empty sets, which would only be complementary to their counterpart in a trivial way. Since the fractures themselves have informative morphological structure transcending this scale, we design in this section a way to focus the morphological processing purely on the fracture surface (rather than on fragment volumes). This should suffice for our puzzle, since whether two broken fragments can be refit locally depends on the complementary shape of their fractures only.

3.1 Boundary Morphology

In Section 2, we had masked objects $X_{M}$ and $Y_{M}$, with complementarity in $M$ (and a slight issue of how we treat their common boundary fracture $F$). Now consider only one of them, no longer in contact, and investigate how to compute its openings and closings. We are only interested in the effect on $F$, not on the remainder of the object volume within $M$.

$F$ is a boundary, not a volume, and the classical MM does not apply immediately. We first define what we mean by applying mathematical morphology to a boundary. Let us call it ‘Boundary Morphology’ and denote its operations by an over-bar. Consider a general object $A$ of which we take the boundary $\partial A$ (and ignore the masking effects for the moment). As we have seen, for true complementarity some objects may be open sets, others closed sets. We need to take the boundary of either, so we will always close the set and then apply the usual boundary operator $\partial$ to define our boundary operator (which we denote by $\overline{\partial}\,$):

[TABLE]

Thus $\overline{\partial}\,A$ returns the set of points of $\overline{A}$ that have neighbours in ${\overline{A}}^{c}$.

In order to distinguish erosion and dilation, we need to make the boundary oriented, for instance by denoting the outward pointing normal at each location. Then this orientation changes to its opposite if we consider the boundary of the complement. To more easily denote this, write the complement as the prefix operator $c\,$. We then have:

[TABLE]

which defines the complementation $\overline{c}\,$ of the boundary.

We define what it means to dilate and erode a boundary, defining operators $\overline{\delta}\,$ and $\overline{\varepsilon}\,$ in terms of classical volumetric operators by:

[TABLE]

It now follows that we can perform an erosion as a dilation on the complemented (oppositely oriented) boundary:

[TABLE]

and accordingly,

[TABLE]

In fact, this is merely the duality between dilation and erosion, extended to their boundary versions (with $\overline{c}$ playing the role of complementation). So once we can dilate a boundary, we can use that both to produce the dilation $\overline{\delta}\,$ and the erosion $\overline{\varepsilon}\,$ (by doing dilation on $\overline{c}\,\overline{\partial}\,A$, the oppositely oriented boundary):

This can be easily extended to boundary closing $\overline{\phi}\,$ and opening $\overline{\gamma}\,$ defined as:

[TABLE]

for it follows from Eqs. (7) and (8) that we can also perform an opening as a closing on the complemented boundary:

[TABLE]

and accordingly,

[TABLE]

Now, the schema of eq. (3) can be rewritten to describe the complementarity of two exactly fitting fracture boundaries (at scales $R$ and $\rho$ with $R>\rho$), as follows:

[TABLE]

where $\overline{c}\,$ (complementarity) is taken within $M$.

3.2 Lipschitz Condition

Terracotta is uniform and brittle. It fractures rather simply such that there is at least one 3D direction from which the whole fracture facet is completely visible, i.e. a ray in that direction emanating from each point of the fracture will not encounter any other point of that local fracture facet. Accordingly, a given fracture facet can be represented as a function (a Monge patch); noting that another fracture facet of the same fragment may require a different visibility direction to be seen as a function. The visibility assumption can be efficiently characterized in terms of the ‘Lipschitz condition’ which frequently occurs in quantitative mathematical morphology [17].

For any pair of points on the graph of a slope-limited Lipschitz function with slope $s$, the absolute value of the slope of the line connecting them is always less than $s$ (i.e., the Lipschitz constant). As a consequence, there exists a double cone whose vertex can be moved along any such Lipschitz-continuous function, so that it always remains entirely visible from the principal direction of the cone while the rest of the function is completely outside the cone. For a local fracture to meet the visibility condition, it should meet the Lipschitz condition over its entire domain.

To compute the Lipschitz slope $s$, we first collect the local normal vectors of the fracture facet on a unit sphere of directions, and then fit the largest possible cone inside this set (see Figure 2). The axis of this cone we call the cone direction, and the opening angle we call the cone angle. Let us use the cone direction as the principal direction to align the fracture facet, so that it can be described as a Lipschitz function. This fracture function is then bounded in Lipschitz slope $s$, which is the cotangent of the cone angle.

In practice, the Lipschitz condition appears to intrinsically hold for each of our terracotta fractures. Even if it does not hold, we can always artificially separate the fracture regions by the Faceting preprocessing procedure [5] that delivers the piecewise Lipschitz facets of the fracture region.

3.3 Extrusion

The fracture surfaces, extracted using the Faceting preprocessing operation [5], are locally Lipschitz. The corresponding Lipschitz principal direction acts as the ‘average breakage direction’ for the whole fracture facet. Therefore, we can generate a thickened fracture volume by extruding the facet along such visibility direction without self-intersection. Such an extrusion effect is equivalent to making two copies of the fracture surface bounded by a generalized cylinder: one with original (non-inverted) normals $\overline{\partial}\,F$ and one with inverted normals $\overline{c}\,\overline{\partial}\,F$(see Figure 3). The outward propagation of the copy with non-inverted normals is the fracture surface dilation ($\overline{\delta}_{\rho}$), while the propagation of the one with inverted normals is the erosion ($\overline{\varepsilon}_{\rho}$), by eq. (8). By a subsequent inward propagation of the resulting expanded surfaces, by the same amount $\rho$, one then acquires the closing ($\overline{\phi}_{\rho}$) and the opening ($\overline{\gamma}_{\rho}$) of the fracture.

Focusing our morphological simplification on the extruded fracture volume has the following advantages:

1. (a)

Extending the functional morphological scale (ball radius $\rho$) beyond the minimum thickness of the archaeological object by only needing to do closing. 2. (b)

Permitting simultaneous propagation of the fracture surface to both sides, thus allowing the production of MM opened and closed surfaces in one go. 3. (c)

Avoiding the dilation of elements outside the mask (which could affect the outcome within the mask) by having the extruded surface surrounded by sufficient empty space completely isolated from outside influences. 4. (d)

Avoiding needless processing of the non-fracture regions of the object.

3.4 Distance-Transform based MM implementation

We obtain our data as 3D meshes, so it would be natural to apply the operations of boundary morphology to fracture facets represented directly as extruded 3D meshes. However, applying morphological operations directly on a mesh representation is notoriously hard (the dilation and erosion operations lead to many additional vertices as the mesh faces start intersecting [8], and the ball structuring elements add many more). At this point of our narrative, it would seem that Point Morphology [3] is also a good candidate; we will see in Section 4 that we need more than just the simplified surfaces by themselves to determine our morphological features.

An alternative implementation of morphology, especially when done in a scale space context, is by means of the isosurfaces of distance functions. Since we do want to perform our operations at different scales of a ball-shaped structuring element, this is especially attractive. Two basic approaches exist: distance transform considered as a numerical level set in a space of 1 more dimension [18] or classic distance transform on a 3D grid. We adopted the latter approach which: a) suffices as a straightforward demonstration of the validity of our concepts and b) enables us to keep track of the original fracture points propagation in the distance field and their contribution to the generated MM surfaces through what we call the provenance map.

Our implementation relies on distance transforms in a volumetric representation [20], oriented with the Lipschitz principal direction along its z-direction. We convert the extruded mesh to binary representation by embedding the fracture surface in a 3D voxelized grid of a well chosen resolution, of step size $g$. We derived that the maximum outward and inward displacement that could take place due to discretization effects is no more than $\sqrt{3}g/2$. This implies that distances (and hence MM surfaces) are affected by no more than this amount. The generated grid is padded with an empty region sufficient to contain the maximally dilated versions of the fracture surfaces. For generating our distance field, the method proposed by [13] is employed, which calculates the euclidean distance transform in linear time on a binary voxelized representation of the object as shown in Figure 4. The closing of the fracture volume with a ball of radius $\rho$ (closing is all we need by virtue of eq. (11)) is performed by first extracting the outward level set at distance $\rho$, then re-computing the distance field of the background and extracting the inward level set at distance ${\rho}$.

3.5 MM Surfaces

Intuitively, morphologically closing the entire extruded volume is equivalent to rolling a ball on the fracture surface from both sides. The upper rolling (constrained by the fracture peaks) is equivalent to the closing effect, while the lower rolling (constrained by the fracture valleys) is equivalent to the desired opening effect.

Figure 5 shows the closed and opened simplified MM surfaces at 6 different scales of a given fracture surface with $5.8K$ vertices and $11K$ faces. The grid resolution is $g=0.2$ mm of size $119\times 323\times 29$ and is padded with $151$ voxels from each side to permit the outward propagation of the fracture surface up to $30$ mm without losing information, thus producing a grid of size $421\times 625\times 331$.

The fracture Lipschitz computations together with the extrusion and embedding take less than $2$ seconds for a mesh of size $10K$. The distance field computation time is more affected by the resolution of the grid: computing the distance transform for $500\times 500\times 500$ grid takes at most 1 minute. Extracting the closed and opened MM surfaces takes less than $50$ seconds for the same grid size. In the GRAVITATE system, such times are considered acceptable. We performed our experiments on an Intel(R) Xeon(R) CPU $2.9GHZ$ $\times$ $8$ computer with $256GB$ RAM.

4 Inexact Complementarity

In practice, there are many aspects to our data that prevents complementarity from being exact. However, using the Lipschitz condition for fracture surface characterization together with MM for hierarchical fracture simplification enables us to establish the bounds on the amount of inexactness. In this section, we briefly discuss the sources causing inexactness and and how to treat their associated deviations.

• •

Abrasion*.* When a fractal surface as brittle as terracotta is rubbed, this will tend to take off the sharp local peaks without affecting the valleys. We propose that a reasonable model of such effect is that the object undergoes a morphological opening $\gamma_{\alpha}$ of a small size $\alpha$.

Assuming this, the opening scale space will be unaffected for $\rho>\alpha$, since $\gamma_{\rho}(\gamma_{\alpha}(X))=\gamma_{\rho}(X)$ for $\rho\geq\alpha$. However, the closing is affected by the abrasion. Therefore, exact complementarity of $\gamma_{\rho}(X^{\prime})$ and $\phi_{\rho}(Y^{\prime})$ of eq. (3) schema no longer holds for the abraded version $X^{\prime}=\gamma_{\alpha}(X)$ and $Y^{\prime}=\gamma_{\alpha}(Y)$, even if it did for the original $X$ and $Y$. Nevertheless, the difference is bounded for the Lipschitz functions we are considering as the fracture facets. A simple sketch of a local sphere of radius $\alpha$ capped by a Lipschitz cone of slope $s$ shows that the maximal deviation of a possible original surface and its $\alpha$-opening is (\sqrt{1+s^{2}}-1)\simeq{\mbox{\frac{1}{2}}}\alpha s^{2}. Since $X^{\prime}$ does not differ by more than {\mbox{\frac{1}{2}}}\alpha s^{2} from $X$, also $\phi_{\rho}(X^{\prime})$ will not differ from $\phi_{\rho}(X)$ by more than this amount. With that in mind, we can still test the possibility of complementarity, to within this bound, of the opened and closed abraded fractures.

• •

Fracture Adjacency Effects*.* Each facet is delineated by an outer border (contour). That border defines the complementarity zone of its processed MM surfaces. In practice, there are two different kinds of facet borders: fracture-skin border and fracture-fracture border. The former is the contour’s segment that is adjacent to a facet which belongs to the exterior skin of the original object. By contrast, a fracture-fracture border occurs when the original fragment is further broken into sub-fragments, causing the adjacency of more than one fracture facet.

If the facet has only fracture-skin borders, then complementarity with its counterpart still holds over the entire surfaces of its simplified MM versions across all scales. For facets with fracture-fracture borders, one should expect complementarity to hold within regions away from the border. But there is a scale-dependent zone adjacent to the fracture-fracture borders where the MM surfaces are no longer guaranteed to be complementary to those of its counterpart. In that zone, the representation of either part may have been affected by the secondary breaking. The bounds we can obtain on this zone under the general assumptions we made (such as Lipschitz) are weak. We have found that we can much more specifically delineate it through analysing the provenance map of the distance transform. Our current work is on efficient representation of the simplified scale space surfaces, taking this effect into account.

• •

Misalignment*.* Complementarity only exists in perfect alignment of the two counterparts; any offset in translation or orientation will destroy it. But this effect is not as boolean as Eqs. (1) and (2) makes it appear: the complementarity is quantifiable by means of the MM scale spaces. The distance transform method of producing the openings and closings at different scales can be employed to guide the alignment, from coarser to finer scales, with the separation between one surface’s opening and closing at each scale giving an indication of reasonable bounds on the fine alignment at that scale.

5 Conclusion and Future Work

We have presented morphological simplification of the fracture facets of archaeological fragments based on their scanned 3D mesh representation, as a preparatory phase for optimal pair-wise alignment. This is a problem that lends itself very well to treatment by mathematical morphology, since that framework provides complementarity-preserving simplification of shapes in a manner that is insensitive to the kind of missing information and abrasion that we expect in the archaeological fragments. We showed how to perform morphology on boundaries in Section 3.1. The assumption that each fracture facet is a Lipschitz function allowed us to set up the extrusion method in Section 3.3. The duplicate internal and external copies of the local fracture facet surface enable computing the opening and closing scale spaces simultaneously, in a volumetric representation in Section 3.5.

The resulting surfaces are simpler, but still contain the essentials of complementarity. We are currently investigating the compact characterization of the MM scale space surfaces by means of characteristic scale space medial axis points, computed directly from the distance transform and its provenance map (specifying which points are influential at each location). This should allow the use of standard registration algorithms on those morphologically relevant feature points only. Moreover, we plan to use the provenance map of the distance transform to automatically control the reliability of those representations when parts of the original fracture facet are missing.