UvA-DARE

: When fitting archaeological artifacts, one would like to have a representation that simplifies fragments while preserving their complementarity. In this paper, we propose to employ the scale-spaces of mathematical morphology to hierarchically simplify potentially fitting fracture surfaces. We study the masking effect when morphological operations are applied to selected subsets of objects. Since fitting locally depends on the complementarity of fractures only, we introduce ‘Boundary Morphology’ on surfaces rather than volumes. Moreover, demonstrating the Lipschitz nature of the terracotta fractures informs our novel extrusion method to compute both closing and opening operations simultaneously. We also show that in this proposed representation the effects of abrasion and uncertainty are naturally bounded, justifying the morphological approach. This work is an extension of our contribution earlier published in the proceedings of ISMM2019 [10]


Introduction
The GRAVITATE H2020 project aims at providing archaeologists with the virtual tools to analyse digital artifacts, distributed across several collections (https://gravitate-project.eu/). The project involves archaeological institutes and research partners with experience in 3D geometry processing, semantic analysis and ICT integration. The artifacts 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 properties. Reassembly of broken artifacts 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 [9] which partitions its surface into significant sub-parts called 'facets'. Each facet is characterized by its own geometrical properties of the roughness of its surface and the 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 on the object facets).
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. We could then explore the possible fit and approximate alignment at the coarse level, and refine this when descending to finer resolutions, all at approximately the same computational cost per level. But to be relevant to alignment, such a simplification should be contact preserving and sufficiently specific to structure and guide the pair-wise alignment.
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 [16]. Structure-aware mesh decimation, to preserve the structure of the mesh by representing it as a set of pre-computed canonical proxies, has been studied by [22]. Such simplifications are strongly related to characterizing a surface by differential geometric features, such as those computed by some means of discrete differential geometry [7]. Indeed most methods to reassemble broken objects [12,14,18,20,26] are 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 of corresponding points is not appropriate; some extensions of ICP that include distance transforms [11] 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. Contrary to linear methods, it is not affected by the depth of missing parts, but merely by their spatial scale. Such scale spaces have been studied in [2,3,15], and it has been demonstrated that they simplify the shape, when measured in terms of certain descriptive features such as local maxima [15], or zero-crossings of the curvature [6]. 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 morphology of complementary objects (rather than of 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 new 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. This paper does not perform the actual fitting, but Section 6 outlines how the proposed morphological representation of the fracture surfaces can be employed in optimal pair-wise alignment.
This paper extends and supersedes our work published in [10]. We now give full technical detail in the form of lemmas and proofs of the validity of complementarity properties within a masked region in Section 3.1 and 3.2, preparing this by precise notation in Section 2. Section 6 is also new contribution relative to [10].

Basics of Mathematical Morphology
Mathematical morphology (MM) was developed in the 1960s as a non-linear framework to describe shape in its qualitative (morphological) aspects, with structuring elements acting as probes to describe the shape hierarchically at different scales. MM flourished mostly in image (pre)processing, in which its set of non-linear operations produces results that cannot be obtained by linear, convolutional filtering. This effect-oriented use of MM has obscured the fact that it is actually a powerful language for describing shapes in contact [8].
In this paper, we employ it as such.
Since it is always clear which structuring element we use, we adopt the following notational shorthand. Let X be a connected set in R 3 and S be our structuring element. Then: Reflected opening:ˇ(X) =δε(X) = (X ⊖Š) ⊕Š With this notation, the underlying principles of MM [23] that we will employ can be summarized as follows: (B 1 ) Duality: Erosion (opening) and dilation (closing) are dual operations relative to complementation (denoted X → X c ) of the sets they act on. This means that they satisfy: The dilation is distributive over union and the erosion is distributive over intersection: but the converse statements are weaker: Local Knowledge Principle : Morphological operations satisfy the "local knowledge principle". Shapes are usually seen through a window (e.g, image frame), within which MM is restricted: In words, for any bounded window M in which we want to define a morphological transformation Ψ(X), there exists a bounded window M * , in which the knowledge of X is sufficient to locally perform the transformation Ψ.

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 volume occupies), the two fragments X and Y, viewed as sets, 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.
No Gaps ⇐⇒ Completeness: In our work, one does not acquire the fragments in alignment, or in perfect condition, or even together with their counterpart. We choose MM scale spaces to process the fracture region of each of the fragments separately, producing a simplified representation for each fracture facet of each object at a range of scales ρ. 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 5). 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). The complementarity of two exactly fitting objects (within a Mask M) and their morphological opening ( ) and closing ( ) may thus be summarized with the following schema: Here c M * refers to complementarity within a more restricted 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., ρ(XM ) ≠ ( ρ(X))M ). We refine this procedural sketch in the following section.

Masked MM
Before we move to the scale space complementarity, we need intermediate results on the effects of masking morphological operations when applied to one object only (independent of its complementary counterpart). The setup will be that within a mask M we have only X ∩ M available. We are interested in the largest possible region M * in which the results of our operations are identical to these operations applied to the unmasked X. So our results will have the form: for each operation α, with the post-masking set M * to be determined. Therefore the lemmas will all be specific instances of the principle of local knowledge (B 3 ). For many specific examples of X and M, the masking by M does not really interact with operations on X. For instance, if X is a point at the origin, and M is a large set enclosing the origin, then small dilations will satisfy δ(X ∩ M) = δ(X). In this case M * = M would suffice. However, if just outside M there are elements of the set X, these will dilate into the mask M, and the identity then only holds in a shrunken version of M; in fact, optimally shrunken toε(M), as we will show below in Lemma 3.2. In our typical use of masking, M will select a part of X (it is meant to focus on a fracture facet of an artifact), so the precise delimiting of the validity matters. We are therefore interested in lemmas that hold universally, for all X and M. To formulate this more precisely, let us introduce the weak notion of 'necessity' required to formulate the lemmas.
We will call a parametrized set a 'generally necessary lower bound' for the validity of a theorem involving sets X and M, if for all X and all M, there is no smaller set for which the theorem holds. In order to demonstrate that a particular set is a generally necessary lower bound, it suffices to show a specific choice of X and M in which masking of the result by a larger set would make the equality invalid.

Lemma 3.1. For masked erosion, the generally necessary lower bound is ε(M), i.e. ε(X∩M)∩ε(M) = ε(X)∩ε(M), for all X and M.
Proof. Sufficiency is straightforward from the commutation and distributivity properties of the operations of erosion and intersection: We give a critical example to demonstrate the general necessity in 1D. Consider a mask M = (−∞, 0), the half line left of zero, and an object X = (−∞, 0], an interval that includes the point at zero. Picking as our structuring element the asymmetric 1-sphere B = [−ρ, 0], we find that ε(X) equals (−∞, −ρ], so it includes the point at −ρ. Yet ε(X ∩ M) equals (∞, −ρ), and therefore does not include this point. We can precisely eliminate this point by limiting the post-masking of the results by M * = ε(M), and no more severe limitation than this is needed.
Our counterexample to prove the general necessity of the lower bound was taken from 1D, for clarity. This is more universal than it may seem. It has been shown that the boundary points that can possibly be present in the result of erosions (or dilations) emanate along straight lines from the original set boundaries, as the structuring element scales (see [27], Section IV). In morphology in n dimensions, there will always be at least one critical 1-dimensional propagation direction that determines the lower bound of validity, through the effect most purely exposed in this 1D counterexample. Proof. We first show sufficiency of the bound: The second term is empty by duality between dilation and erosion: The general necessity of the masking byε(M) can be established by giving an example in which no larger mask (i.e., smaller eroded mask) would work to preserve the identity. In  Proof. Sufficiency: Proof. Sufficiency: Actually, the post-masking validity regions of opening and closing are identical, since we have: For notation convenience, we will denote both εε andεε by ε 2 , which better shows their equality. In fact the resulting structuring element has interesting properties: Lemma 3.6. The structuring element S ⊕Š involved in ε 2 is symmetrical and contains the origin. This implies that the post-masking region for opening and closing only shrinks the sets symmetrically but does not displace them, even for general S.

Masked Complementarity
In this section we prove that when originally there was exact complementarity between X and Y within a mask M (i.e, X c M = Y M ), then exact complementarity still holds between (X) andˇ(Y), or between (X) anď (Y), but is only guaranteed within a doubly eroded mask M * = ε 2 (M). This is important to establish before we move into MM scale space, where the operations will be performed with scaled version of the structuring element S. It shows precisely how the validity region depends on S and its scale.  Proof.
In contrast to this universal validity of non-intersection within M, the completeness property needs to be restricted, as shown by the following lemmas.
(by Lemmas 3.2 and 3.1) Necessity in the sense of 'generally necessary lower bound' follows from the general necessity of the erosion and dilation masking of Lemmas 3.1 and 3.2. Proof. Sufficiency of limiting by ε 2 (M):

Ball Scale Spaces for Archaeology
We will not need MM in its full generality, which would allow one to choose a structuring element for all basic operations to affect the directions and extent of dilation and erosion (and hence of closing and opening). In our archaeological context, we assume that the objects are placed in arbitrary orientation, and that abrasion (by wear or wind) has no preferential direction. As a consequence, we should perform all our morphological operations with structuring elements that are balls (filled spheres), of a specified radius ρ as the effective scale. We now refine our notation to be specific on the radius and use ρ as a subscript to our δ as in δρ. The balls we use are convex, symmetric and contain the origin; this simplifies our treatment since we no longer need to distinguish between S andŠ, and consequently ε andε etc. At a scale ρ within a suitably chosen mask ε 2 ρ (M), the opened version ρ(X) will be complementary to the closed version ρ(Y ), and vice versa. The closing ρ(X) simplifies the local geometry of the fracture X at its valleys, and the opening ρ(Y ) does the same to Y at its peaks. Therefore, the complementarity of ρ(X) with ρ(Y ) becomes less specific in those areas, and vice versa. To maintain specificity of both peaks and valleys of the common fracture surface, we must therefore compute both opening ρ(X) and closing ρ(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.
For different scales R and ρ with R > ρ, the complementarity and containment relationships illustrated by schema 3 can now be extended as follows: Here cρ is a shorthand for c ε 2 ρ (M) and ε 2 ρ (M) = ερ(ερ(M)) =ερ(ερ(M)), which is equivalent to erosion of M by a ball of twice the radius (ε 2 ρ = ε 2ρ ). 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 ρ of more than half this thickness would results in empty sets, (see Figure 2), 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 the fragment volume). This should suffice for our puzzle, since whether two broken fragments can be refit locally depends on the complementary shape of their fractures only.

Boundary Morphology
In Section 3, we had masked objects X M and Y M , with exact complementarity within a chosen mask 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 mathematical morphology 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 ∂A (and ignore the masking effects for the moment). As we have seen, for true mathematical set 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 ∂ to define our boundary operator (which we denote by ∂ ): Thus ∂ A returns the set of points of A that have neighbours in 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. We assume that the result of the boundary operator ∂ includes this allocation of normal directions (and hence so does ∂ ). This orientation should change to its opposite when we consider the boundary of the complement. To more easily denote this, we write the complement as the prefix operator c . We then have: which defines the complementation c of the boundary. We define what it means to dilate and erode a boundary, defining operators δ and ε in terms of classical volumetric operators by: δ ∂ A ≡ ∂ δA and ε ∂ A ≡ ∂ εA.
It now follows that we can perform a reflected boundary erosion as a boundary dilation on the complemented (oppositely oriented) boundary: which implies thatε = c δ c andδ = c ε c .
In fact, this is merely the duality between dilation and erosion, extended to their boundary versions (with c playing the role of complementation). So once we can dilate a boundary, we can use that capability to produce both the dilation δ and the reflected erosionε (by doing dilation on c ∂ A, the oppositely oriented boundary, and complementing the result): This can be easily extended to boundary closing and opening defined as: for it follows from eqs. (9) and (10) that we can also perform a reflected boundary opening as a boundary closing on the complemented boundary: which implies thatˇ= c c andˇ= c c .
For symmetric structuring elements like balls, we can ignore the checkmarks within the above MM notations. Now let us relate this to each of the objects X and Y within a mask M. The mask is chosen to focus on (part of) the common fracture F. With X and Y in perfect contact, the set F is simply F = X ∩ Y ∩ M (the overbar denoting set closure). Considering only the object X in isolation, for the contact we focus on the part of its boundary masked by M. We turn this point set into an oriented boundary facet for X, borrowing the outward pointing normals of X at the chosen locations. This defines a facet F X : Some care is needed: this definition produces only points of the desired F if M itself does not contribute to the boundary. Therefore, we have to take the mask M as an open set. Similarly, one defines F Y ≡ ∂(Y ∩ M). One would then like to define the boundary dilation δ through: or, denoting the masking by M as an operator, which looks like the natural extension of the unmasked case of eq. (8). But after dilation, the masking by a general M could crop parts of the dilated boundary facet that should have been included, since they were caused by points of the original boundary within the mask. This issue can be resolved simply, by constructing masks M for boundary MM only as the union of normal line segments emanating from the boundary facet F X , since the differential equations of morphology transport points along straight lines in the direction of the normals [27]. With F X defined by eq. (14), we have produced a proper oriented boundary facet, which can take the place of ∂ A of eq. (6) in the derivations above. We can then define boundary morphology in this manner, with the complementarity relationships like F Y = c F X indicating that F Y and F X are the same set of points, but with opposite orientation. This permits us to translate the volumetric schema 5 to describe the complementarity of two exactly fitting fracture boundaries F X and F Y (at scales R and ρ with R > ρ), as follows: Here cρ is the boundary complementation within a mask ε 2 ρ (M). Terracotta is uniform and brittle. We assert that 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) within an orthogonal coordinate frame taking a visibility direction as function axis (while noting that another fracture facet of the same fragment may require a different coordinate frame to be seen as a function). By choosing our frames well, the visibility assumption of the surface can be related to the 'Lipschitz condition' of the corresponding function [24].

The Lipschitz Condition
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 maximum Lipschitz slope s, we first collect the local normal vectors of the fracture facet on a unit sphere of directions, and then fit the smallest possible cone that contains this set (see Figure 3). 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. The corresponding Lipschitz principal direction acts as the 'average breakage direction' for the whole fracture facet. The 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 fracture facets with a Lipschitz slope of approximately 1. Even if it does not hold, we can always artificially separate the fracture regions by the Faceting preprocessing procedure [9] that delivers the piecewise Lipschitz facets of the fracture region. The fracture surfaces can thus be guaranteed to be locally Lipschitz. 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 ∂ F and one with inverted normals c ∂ F (see Figure 4). The outward propagation of the copy with non-inverted normals is the fracture surface dilation (δρ), while the propagation of the one with inverted normals is the erosion (ερ), by eq. (10). By a subsequent inward propagation of the resulting expanded surfaces, by the same amount ρ, one then acquires the closing ( ρ ) and the opening ( ρ ) of the fracture.

Extrusion
Focusing our morphological simplification on the extruded fracture volume has the following advantages: (a) Extending the functional morphological scale (ball radius ρ) beyond the minimum thickness of the archaeological object by only needing to perform the closing operation.  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 not straightforward (the dilation and erosion operations lead to many additional vertices as the mesh faces start intersecting [13], and the ball structuring elements add many more). At this point of our narrative, it would seem that Point Morphology [5] is also a good candidate; however we will see in Section 5 that we need the structure of the distance field to determine sensible morphological features for pair-wise alignment, rather than just the simplified surfaces by themselves.

Distance-Transform Based MM Implementation
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 one more dimension [25] 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. This map specifies the original boundary point(s) to which each point of the MM surface can be traced back. Our implementation relies on distance transforms in a volumetric representation [28], 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. The maximum outward and inward displacement that could take place due to discretization effects is no more than √ 3 g/2, since each point is displaced by maximum half the grid size along each of its dimensions. 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. To generate our distance field, the method proposed by [17] is employed, which calculates the euclidean distance transform in linear time on a binary voxelized representation of the object as shown in Figure 5. The closing of the fracture volume with a ball of radius ρ (closing is all we need by virtue of eq. (13)) is performed by first extracting the outward level set at distance ρ, then re-computing the distance field of the background and extracting the inward level set at distance ρ. This generates a maximum error of √ 3 g for the extracted surface relative to the true closing. Intuitively, morphologically closing the entire extruded volume is equivalent to rolling a ball (our chosen structuring element) over 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. At increasingly coarse scales, we compute the closed and the opened versions of the fracture facet. Moving up the scale space, the generated MM surfaces are simpler in shape than the original fracture they represent while maintaining the fitting (i.e., complementarity) essentials. Figure 6 shows the closed and opened simplified MM surfaces at 6 different scales of a given fracture facet with 5.8K vertices and 11K faces. The grid resolution is g = 0.2 mm of size 119 × 323 × 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 × 625 × 331.

MM Surfaces
The fracture facet 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 O(g −3 ): computing the distance transform for 500 × 500 × 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 since the morphological characterization of the fracture facets of each artifact is an o ine process at ingestion time. We performed our experiments on an Intel(R) Xeon(R) CPU 2.9 GHz × 8 computer with 256 GB RAM.

Inexact Complementarity
In practice, there are many aspects to our data that prevent complementarity from being exact. Fortunately, 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 how we expect 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 surface undergoes a morphological opening α of a small size α.
Assuming this, the opening scale space will be unaffected for ρ > α, since ρ ( α (X)) = ρ (X) for ρ ≥ α. However, the closing is affected by the abrasion. Therefore, exact complementarity of ρ (X ′ ) and ρ (Y ′ ) of the schema in eq. (17) no longer holds for the abraded versions X ′ = α (X) and Y ′ = α (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 α capped by a Lipschitz cone of slope s shows that the maximal deviation of a possible original surface and its α-opening is ( √ 1 + s 2 − 1) ≃ 1 2 αs 2 . Since X ′ does not differ by more than 1 2 αs 2 from X, also ρ (X ′ ) will not differ from ρ (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 subfragments, 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 then 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 further 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 (defined in Section 4.4). 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 [11].

Towards Fitting
The classical way of fitting fracture surfaces is based on matching significant corresponding features characterizing the fracture surfaces of counter-fitting fragments [4,14,21]. These features are mostly detected and extracted using differential geometric techniques which are linearly sensitive to missing information, as has been clarified before in Section 1. Such features are not expected to be stable under abrasion and many of them are not compact enough to allow the cheap hierarchical fitting we are anticipating. In this section, we introduce our intended morphological features and how they can be structured in a compact scale space representation that can robustly and cheaply guide our next phase of optimal pairwise alignment.
The key observation is that the closing process results in a surface of which large parts are spherical caps (see Figure 6, especially the coarser resolutions, and Figure 7). Moreover, as the scale increases, the surface tends to be hierarchically composed of fewer, larger and more clearly separated spherical caps. The centers of the balls to which these spherical caps belong can be extracted from the provenance map of the distance transform computation (which specifies which point on the original fracture generated a given point on the MM scale space surface). In fact, these centers are medial axis points of the complement of the extruded fracture mesh volume. Medial axes are notoriously sensitive to discretization, but they can be made robust [1]. With increasing scale, the ongoing morphological simplification can be characterized in terms of fewer yet more distinctive medial axis feature points, which can therefore serve as 3D morphological features. Moreover, these points can be tagged by their significance for the MM reconstruction (see Figure 7(a)). Figure 7(b) is an exploration of this idea. It shows the 30 most distinctive MM feature points which are capable of representing 30% of the simplified (closed) fracture facet at scale 30 mm. The fracture facet is the same as in Figure 6, which is about 6 × 2.5 cm size. Moving up in the hierarchical scale space, these MM feature points are increasingly capable of describing more surface area of the overall simplified fracture surface. It seems natural to use each cap's area for weighting the significance of the corresponding feature point. We thus arrive at the representation of a fracture surface by a set of 3D weighted points, at well-defined positions above it (for closing) and below it (for opening), at a certain scale ρ. Such 'landmarks' may then be used to characterize the morphologically simplified versions of the fracture surface at different scales. We could further reduce the morphological characterization by only keeping the medial axis points at critical scales, as is usual in linear scale spaces. We would then consider replacing our linear time Euclidean distance transform (based on [17]) by the linear time medial axis chamfer approximation of [19].
According to our simplification framework, a morphological graph characterizing the closing scale space of a given fracture surface is expected to be very similar to the morphological graph characterizing the opening scale space of its counter-fitting fracture surface. Alignment of fracture surfaces using morphological graphs can be done by means of using the well-known Procrustes method. One should then use an extended version of this algorithm which allows for missing point correspondences. But before we perform our alignment in this way, we are currently studying the reproducibility, stability and robustness properties of the morphological feature points and the feasibility of computing them directly from the distance field.

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 archaeological fragments. We showed what the effect of masking is on the interpretability of MM results in Sections 3.1 and 3.2 and also how to perform morphology on boundaries in Section 4.1. The assumption that each fracture facet is a Lipschitz function allowed us to set up the extrusion method in Section 4.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 4.5. The resulting surfaces are simpler, but still contain the essentials of complementarity for identification and fitting of counterparts. In Section 5, we discussed the sources of inexact complementarity and how their effects are naturally bounded in our method.
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) along the lines sketched in Section 6. 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.