- Email: [email protected]
Similarity measures for blending polygonal shapes
Computers & Graphics 25 (2001) 13}27
Shape Blending
Similarity measures for blending polygonal shapes Michela Mortara, Michela Spagnuolo* Istituto per la Matematica Applicata, Consiglio Nazionale delle Ricerche, Via De Marini 6, 16149 Genova, Italy
Abstract This paper presents a solution to the correspondence problem for polygon blending, based on a morphological characterization of the shape. The approximate skeleton is suggested for coding the shape as a graph of meaningful areas, and the correspondence problem is therefore solved as a skeleton matching process. The correspondences are found taking into account the global structure of the shapes, and also local similarity criteria are used to further re"ne the matching. The approximated skeleton resembles the medial axis, but it is much simpler to extract and may be easily pruned to reduce the e!ects of small scale details of the shapes' boundary. The "rst implementation results con"rm the validity of the approach for driving the blending both in surface reconstruction or morphing contexts. 2001 Elsevier Science Ltd. All rights reserved. Keywords: Skeleton; Shape similarity; Similarity measures; Feature recognition; Polygon blending
1. Introduction The problem of smoothly transforming, or blending, one shape into another is an intriguing problem with several di!erent applications, ranging from the entertainment industry to engineering applications. In particular, object-space blending deals with the transformation of an object representation into another one. When shapes are modelled as polygons or polyhedra, the blending process is generally seen as composed of two basic steps: "rstly, correspondences among vertices of the source and target models are set, and then an interpolation function is de"ned to produce as many intermediate shapes as required by the application. One of the major problems is to produce intuitive and pleasant intermediate shapes, with a blending which naturally transforms features of the source into similar features of the target shape, avoiding self-intersections of the boundaries or sudden shape changes in intermediate steps [1,2]. Moreover, the blending should be accomplished with minimal interaction of the user.
* Corresponding author. Tel.: #39-01-0647-5677; fax: #3901-0647-5660. E-mail address: [email protected] (M. Spagnuolo).
In a slightly di!erent sense, blending polygonal shapes is necessary for reconstructing surfaces from cross sections. In this case, cross-sections are considered pair-wise to de"ne a tiling among them, that is, a triangular mesh whose edges join vertices belonging to the considered sections. The result of the reconstruction obviously depends on the sampling density, that is, the distance between sections. Since many surfaces could give rise to the observed cross sections, an optimisation is performed with respect to some objective function in order to choose the best surface among the possible ones. Many di!erent goal functions have been used, mainly based on local criteria of the geometric entities involved [3,4]. Blending and tiling polygonal shapes, therefore, exhibit similar problems especially when the shapes to blend are rather di!erent, either in a topological sense (e.g. polygons with di!erent genus) or in a morphological sense. This situation occurs quite frequently in morphing, while in surface reconstruction it generally corresponds to a severely underconstrained tiling problem. In both cases, it is important to have methods for automatically setting correspondences among features of the source and target polygons. Several solutions have been proposed for this problem. For example, physically based methods are aimed at minimising some deformation functions, as in [5] where the minimization concerns the work needed to
0097-8493/01/$ - see front matter 2001 Elsevier Science Ltd. All rights reserved. PII: S 0 0 9 7 - 8 4 9 3 ( 0 0 ) 0 0 1 0 4 - 7
14
M. Mortara, M. Spagnuolo / Computers & Graphics 25 (2001) 13}27
stretch and bend edges and angles along the contours, determining one-to-one correspondences between the vertices. A similar concept has been used in [6] for minimising a curvature-based deformation work function. In industrial design context, a technique has been proposed for locating signi"cant vertices along contours, which represent shape features [8]. Then, local similarity criteria are established for vertices and elements, based on vertex location, curvature, incident and emitting angles and element orientation. For the shape characterization of curves, curvature primitives have been proposed in [7], which could lead to a feature-based and multi resolution contour matching. All these methods fail to take into account the global structure of the shapes for the location of correspondent vertices, and also do not consider the inner part of the shape whose transformation highly in#uences the blending result, as pointed out in [2]. For this aim, it is important to de"ne suitable descriptions of the shape, which may be used to characterize both the inner and boundary of polygons in a computationally usable form. Some works following this principle use structures like modi"ed Voronoi diagrams [10] or a global matching between smooth curves [11] for this purpose. In particular, Shapira and Rappoport [12] use a planar graph obtained by joining the star points of each star-shaped polygon in which the original shape can be divided. Their approach interpolates shapes by "rst interpolating between the skeletons and then unfolding the star pieces from the skeleton. The star skeleton is appropriate for the shape blending problem for two reasons: "rst, it represents all points in the shape's interior as well as on its boundary. Second, it de"nes an explicit dependence between interpolated polygon vertices relatively to a common structure. In any case, the method requires the polygons to have compatible skeletons, that is to say, they must be isomorphic. We have initially considered the problem of shape description for the reconstruction of surfaces from pro"les, for digital terrain modelling and reverse engineering applications. In these contexts, surface reconstruction have been accomplished in two steps: "rstly, the characterization of pro"le shapes is performed, using a multiresolution approach to detect global and local curvature features; then, surface characteristic lines are extracted by joining similar pro"le features, and are used as constraints in the mesh "tting step [13,14]. With regard to polygons, we have used a similar method to e!ectively render generalized surfaces represented by a topological graph, which codes the only adjacencies between critical sections of the original surface. Therefore, the polygons to tile are generally rather distant and may be arbitrarily shaped [15]. In this paper, we suggest the use of an approximate skeleton to describe the shape of polygons and to drive the vertex correspondence phase. The idea we would like
to pursue is to "nd a reasonable match among the approximate skeletons of the polygons to blend, possibly even if the skeletons are quite di!erent. The approximate skeleton of the polygons, which resembles the medial axis, is much simpler to extract and less sensible to convexities on the boundary. With respect to the star skeleton representation, moreover, the approximate skeleton induces a shape decomposition into morphologically meaningful components and is really e$cient to compute. By using global shape structures, it is possible to roughly align the shapes according to their relative direction of best similarity, set global correspondences among shape, and further re"ne the matching, at a lower level, by adding correspondences due to local geometric comparison. The reminder of the paper is organized as follows. In Section 2, the approximate skeleton is introduced and compared with the medial axis transformation; a pruning method for discarding irrelevant details in the resulting shape description is also sketched. Then, the algorithm used to match approximate skeletons is described. In Section 4, some examples are shown of vertex correspondences extracted with our method, and simple blending using linear interpolation of intermediate shapes are given. Conclusions and on-going developments are "nally discussed.
2. Shape characterisation: the approximate skeleton Two-dimensional shapes have been intensively studied in many disciplines, ranging from pure mathematics to psychology, with the aim of de"ning high-level descriptions for the recognition, classi"cation and processing of objects. Characterising a shape means constructing a computational description of the most representative features of the shape, and among the various approaches taken to solve this problem the most suitable one for our purpose is the structural decomposition. According to this paradigm, an object shape is coded as a collection of parts, usually of few basic types, along with their spatial relationships. Therefore, structural decomposition may be seen as a global approach to shape description and the comparison among shapes is generally handled as a labelled graph matching process. Generally, graph matching is a di$cult computational problem, which can be simpli"ed by embedding the problem in its geometric domain. The basic shape elements, or features, to be used for the decomposition should show some properties [16]: shape elements should not change under translation, rotation or scaling (invariance); the description should be locally insensitive to modi"cation of the shape occurring far from the current focus (rich local support); and a small perturbation of the shape should produce only a small perturbation of the description (tolerance to noise).
M. Mortara, M. Spagnuolo / Computers & Graphics 25 (2001) 13}27
Unfortunately, it is rather di$cult to de"ne shape descriptors which ful"l all these requirements. 2.1. The medial axis Among the existing techniques for describing polygonal shapes, the medial axis (MA) is generally considered the more elegant and e!ective one. An intuitive de"nition of the skeleton in the continuum was given by Blum [17], who described the skeleton by analogy with a "re front which starts at the boundary of the shape and propagates isotropically towards the interior. The medial axis is de"ned by the locations at which the "re fronts collide. More formally, the medial axis of a shape S is the locus of centres of all maximal discs of S, that is, those discs contained in S which are not contained in any other disc in S. Equivalently, if B(S) is the boundary of S, then the medial axis of B(S) is the set of points in S having at least two nearest neighbours on B(S). The medial axis, together with the radius function, i.e. the distance from each point on the axis to the nearest point on the boundary, de"ne the medial axis transformation (MAT). The power of this representation is that the shape's boundary and its MAT are equivalent and one can be computed from the other, therefore a two-dimensional object is e!ectively transformed into a one-dimensional graph-like structure. The radius function and the graph-like structure of the MAT induce a decomposition of the shape into subparts which correspond to an intuitive notion of shape protrusion and branching. The points on the axis may indeed be classi"ed according to the number of corresponding nearest neighbours on the shape boundary (see Fig. 1): end points have a single contact point with the boundary, normal points have two contact points and branch points have three or more contact points. A branch point indi-
15
cates a split of the shape into diverging protrusions, while end points indicate a local maximum of curvature of the shape boundary. Therefore, the medial axis can be seen as a graph whose nodes are the end and branch points, and whose arcs are de"ned by normal points. If the shape is a polygon, the MAT is a tree-like planar graph whose arcs are composed by straight-line segments and portions of parabolic curves, and algorithms have been de"ned to compute the MAT from the Voronoi diagram of the polygon elements in O(n log n) time, where n is the number of edges in the polygon [18]. Finally, if the shape is simply connected, then its medial axis graph is a tree, while cycles will appear in the MAT around each hole in the shape. 2.2. The approximate skeleton The idea of approximate skeleton has been used in the "eld of surface reconstruction as a tool for localizing characteristic lines of the surface from a set of its contours [19,20]. Intuitively, an approximate skeleton is generated using the constrained Delaunay triangulation of the polygon as support structure and then by tracing medial segments of the inner triangles, according to rules which take into account the number of constrained edges of each triangle. In this way, a kind of skeleton is produced which globally has the same descriptive e!ect of the medial axis, in the sense that it induces a polygon decomposition in terms of prutrusion-like features. Moreover, the extraction of the approximate skeleton is computationally e$cient. More precisely, let S be a simply-connected polygonal shape having as boundary the polygon P"+P " i"1,2, n, P 3R,; let also ¹ (P) be the DeG G " launay triangulation of S constrained to its boundary
Fig. 1. The medial axis of two polygonal shapes.
16
M. Mortara, M. Spagnuolo / Computers & Graphics 25 (2001) 13}27
Fig. 2. The Delaunay triangulation of a polygon, the triangles having two constrained edges are depicted in dark grey, and the triangles with three inner edges are coloured in light grey.
Fig. 3. The delineation of the AS: some steps of the construction (a) and the "nal AS (b).
[21]. The triangles of ¹ (P) may have either one or two " constrained edges, i.e. edges lying on the boundary P, or may have three non-constrained internal edges (see Fig. 2). Triangles with three inner edges de"ne branch sites of the skeleton, triangles with two edges on the boundary identify terminal nodes of the skeleton, and "nally triangles having only one edge on the boundary contribute to an arc of the skeleton. The AS can be easily traced using medial lines, starting either from branch sites or terminal nodes as follows. With reference to Fig. 3, where the polygon P and its constrained Delaunay triangulation are depicted, let us start the process with any triangle ¹ identifying a terminal node of AS. If v is the vertex of G ¹ shared by the two boundary edges e and e , then the G "rst part of the AS arc ending in v is de"ned by (v, p) where p is the middle point of the opposite edge e . Now, the delineation of the AS continues in the adjacent triangle ¹ according to its type: H E ¹ has no boundary edges: a branching in the AS occurs, H which causes a branching node and two new arcs being
de"ned; the branching node is represented by the centre of mass of ¹ , c, which is joined to the middle H points of its three edges (points q and r in Fig. 3); the tracing continues recursively on the two newly generated arcs, (c, p) and (c, r). E ¹ has 1 boundary edge: the current arc goes on, by H joining the previous one with the middle point of the next non-boundary edge as, for example, (r, t) in Fig. 3. E ¹ has 2 boundary edges: this condition implies that the H current AS arc is terminated, by drawing the segment joining the current arc portion to the vertex common to the two boundary edges. Therefore, the AS can be formally de"ned as follows. De5nition 1. Given a polygonal shape S, with boundary de"ned by the polygon P"+P " i"1,2, n, P 3R,, G G and given the Delaunay triangulation of S constrained to its boundary ¹ (P), the approximate skeleton of S " is the graph AS(S)"(N, A) de"ned by the node set N"+¹ " ¹ 3¹ (P) and ¹ has either no constrained G G " G edges or two constrained edges, and by the arc set
M. Mortara, M. Spagnuolo / Computers & Graphics 25 (2001) 13}27
17
Fig. 4. The MAT of a polygon (a) and the approximate skeleton (b) of the same polygon.
A"+(¹ , ¹ ) " there is a medial path in ¹ (P) between G H " ¹ and ¹ ,. G H The algorithm for extracting AS has to be slightly modi"ed if we have to compute the skeleton of a nonsimply connected polygon. In this case, a starting triangle with two boundary edges may not even exist. If a triangle exists with no edges on the boundary, then the process can start by generating a branching node with degree 3 in the skeleton, and by proceeding recursively on the three generated directions; otherwise, the construction starts from any triangle with one edge on the boundary, with the creation of a dummy node with degree 2, which will be deleted once the skeleton delineation is completed. The computational complexity of the algorithm for the AS extraction is optimal: the choice of the start triangle requires at most ¹ tests, where ¹ is the number of triangles; then, the graph is constructed performing simple operations for each visited triangle, and each triangle is visited exactly one time. Therefore, the algorithm requires O(¹) time. 2.3. Shape description using the approximate skeleton Let us now discuss on the use of the approximate skeleton for describing the shape of a polygon. First of all, the AS is a planar graph having at most degree 3, and this property will be used in the similarity matching step. Compared to the medial axis, the AS is combinatorially simpler in the sense that it generally has fewer arcs and nodes than the MAT. Also, the geometric arcs of the AS are made of line segments and therefore queries like edge length are very easy to implement and fast to compute. In
Fig. 4, an example is shown of the medial axis and the approximate skeleton computed on the same polygon. The decomposition induced by the AS nicely identi"es branching and protrusions of the polygon shape, similarly to that induced by the medial axis. With reference to Fig. 5(a), each arc ending in a terminal node identi"es a protrusion-like feature of the shape (see the portions of shape coloured in grey); each arc between two branching nodes de"nes an inner part of the shape. Therefore, the AS can be e!ectively used to describe the shape of a polygon by attaching to its elements appropriate measures of the components they represent. More precisely, each arc is linked to the area of its component, to the length of the corresponding medial path and to an `average orientationa of the component de"ned as the segment joining the two nodes de"ning the arc (see Fig. 5(b)). Moreover, the decomposition induced by the AS can be easily organised hierarchically according to the scale of the resulting components, following a technique similar to the one used for the medial axis [16,22,23]. Obviously, the approximate skeleton is invariant with respect to shape rotation, and provides a rich local support. However, the MAT and the AS behave similarly with respect to small ondulations or perturbations of the boundary, that is, small modi"cations of the boundary may produce a change in the skeleton structure. This e!ect is shown in Fig. 6 with a simple example. Even if this behaviour is negative for the sake of shape description, the metrics used for characterizing shape components can be used to discard irrelevant details kept in the skeleton. It is straightforward, indeed, to de"ne algorithms for pruning the AS so that only signi"cant arcs are
18
M. Mortara, M. Spagnuolo / Computers & Graphics 25 (2001) 13}27
Fig. 5. The decomposition induced by the AS: protrusions are coloured in light blue and inner components in light orange (a); the dashed segments represents the average orientation associated to arcs (b).
Fig. 6. Small changes in the polygon boundary a!ect the skeleton structure: the original polygon, the medial axis and the approximate skeleton, and the pruned approximate skeleton.
kept in the shape description [23]. An example of this pruning method is depicted in Fig. 6. Intuitively, at the "rst step, protrusions composed by a single triangle with relative small area (depending on the scale and on the average triangle area) are eliminated, and their parent branching node disappears in the AS. The process can be iterated and the size of arcs to prune can be tuned to match the required tolerance. In Fig. 7, the pruning results are shown, where the same shape is represented with a di!erent sampling density, which produces small modi"cations in the AS. As the number of points increases on the boundary, the number of small and irrelevant arcs increases as well, and the pruning algorithm is used to shave away non-signi"cant arcs. In conclusion, the approximate skeleton provides a description of the polygon shape which represents not only boundary features but also takes into account the inner of the polygon.
3. Shape similarity evaluation: structure and geometry Having de"ned a method for describing a polygonal shape, the idea is to establish a set of correspondences between the polygon's boundaries by reasoning on their structural similarity. Eventually, this process will support blending between a source and a target polygon, so that the shape features of the source polygon evolve, if possible, in corresponding similar features of the target polygon. First of all, given two polygons P and P and their respective approximate skeletons AS and AS , it is quite natural to consider that P and P are structurally equiv alent if AS and AS are isomorphic as graphs. There fore, "nding correspondences among shape components is related to the problem of establishing not only a global isomorphism, but also partial matchings among subgraphs. If the problem of establishing isomorphisms among graphs is a complex problem, subgraph matching is even more complicated to be handled in the general
M. Mortara, M. Spagnuolo / Computers & Graphics 25 (2001) 13}27
19
Fig. 7. In the "rst column, the same shape is coded with di!erent sampling; in the second, the computed approximate skeletons and their pruning in the third column.
case [24]. Here, we present an heuristics which strongly uses the geometric embedding of the graphs. For reasoning on skeleton similarity, however, we require that P and P are at least skeleton-comparable, that is, they have at least one node of degree 3; otherwise, the correspondence among features is established using a simpler point-wise approach, as described in the following subsections. 3.1. Skeleton similarity The shape similarity among polygons is evaluated using the following heuristics. First of all, we observe that the branching sites of the approximate skeleton can be seen as the kernel components of the shape, and therefore the comparison should primarily consider branching nodes. A branching node in#uences the shape of a polygon according to the `weighta of the feature components starting at that node, features which are characterized by a geometric signature. Therefore, we start with the de"nition of a measure of similarity for arcs, based on a combination of geometric quantities which take into account the size and the relative orientation of the components. Then, the similarity of branching nodes is de"ned in terms of the combination of their incident arcs which maximizes the sum of the related similarity contributes. Finally, the shapes are compared by searching the &&most similara branching nodes in the two skeletons, and then visiting the graphs from these matched nodes, following each pair of similar arcs until nodes belonging to the
boundary are reached. In this way, correspondences between nodes with degree one, those lying on the boundary, can be established. More precisely, the similarity between arcs is evaluated as a composition of some geometric properties, as follows. De5nition 2 (Similarity between arcs). Let e and e be two arcs of two approximate skeletons; the similarity between arcs, Sim (e , e ), is de"ned as ? Sim (e , e )"aS (e , e )#bS (e , e ) ? ?LEJC JCLERF #dS (e , e )#cS (e , e ), (1) ?PC? RWNC where E S (e , e )"(1803!h)/1803, where h is the angle ?LEJC between the average orientation of e and e . E S (e , e )"l /l if l (l , l /l otherwise, where JCLERF l and l are the length of medial paths corresponding to e and e . E S (e , e )"A /A if A (A , A /A otherwise, ?PC? where A and A are the areas of the two shape components corresponding to e and e. E S (e , e )"1 i! the vertices of e and e have the RWNC same degree, 0 otherwise (see Fig. 8). The parameters a b d c are weights which can be changed according to the requirements of the speci"c input. In particular, it is interesting to consider the e!ect of setting to zero the orientation weight, a. In surface
20
M. Mortara, M. Spagnuolo / Computers & Graphics 25 (2001) 13}27
Fig. 8. Entities used to evaluate the similarity among arcs of the approximate skeleton: in (a) and (b) the average orientation and the medial paths are depicted; in (c) the orientation di!erence between the protrusion alignment.
reconstruction problem, for example, the position of the polygons is "xed and therefore it is reasonable to consider the relative alignment of the shape components to judge their similarity. Conversely, in morphing, similar features should be matched even if the shape is rotated about the branching node. In this latter case, it is useful to discard the contribute due to the relative orientation. These aspects will be discussed also in Section 4. The similarity between branching nodes is computed by evaluating the similarity between each pair of arcs incident in the nodes. De5nition 3 (Similarity of branching nodes). Let A and B be two nodes of degree three of two approximate skeletons; let a , a , a (resp. b , b , b ) be the arcs inci dent in A (resp. B) ordered around A (resp. B). The similarity between A and B, Sim (A, B) is de"ned as L Sim (A, B)"max+(Sim (a , b )#Sim (a , b ) L ? ? #Sim (a , b )), (Sim (a , b ) ? ? #Sim (a , b )#Sim (a , b )), ? ? (Sim (a , b )#Sim (a , b )#Sim (a , b )), ? ? ? (Sim (a , b )#Sim (a , b )#Sim (a , b )), ? ? ? (Sim (a , b )#Sim (a , b )#Sim (a , b )), ? ? ? (Sim (a , b )#Sim (a , b ) ? ? # Sim (a , b )),/distance(A, B). (2) ? Thus, similarity is evaluated by taking the maximum of the sum of the values of incident arc similarity for all the
possible combinations of arcs couples, taken consistently with the ordering around their node. Again, the approach to assess the similarity should take into account the context in which it is performed. If the polygon position is "xed, than the match has to be pursued without allowing any rotation of the shapes. In other words, the orientation of the incident arcs has to be considered and has to be related to a common reference system, in order to avoid torsions in the reconstructed surface. Conversely, if the polygons can rotate during the transformation, the best matching pair can be searched in a slightly di!erent way: for each permutation, the "rst arc couple is aligned, which de"nes a rotation of one shape with respect to the other; then, the relative alignment of the other arcs is evaluated in the rotated reference frame. The parameter distance(A, B) represents the Euclidean distance between points A and B in the plane. If the polygons' position is not "xed a priori, imagine that you move the polygons till their centres of mass coincide, and then the distance between A and B can be computed. The "rst step of the shape comparison process consists in "nding a starting pair of similar branching nodes of the two polygons. Since this pair will be used as the starting point for visiting the two graphs in parallel, it is reasonable to search for the most similar branching nodes in the two shapes. Therefore, the branching nodes which maximise (2) will be chosen among all the possible node combinations. In Fig. 9, an example of this "rst step is shown. Having found the "rst pair of similar branching nodes and having established, at the same time, the three pairs of similar incident arcs, the graphs are recursively visited
M. Mortara, M. Spagnuolo / Computers & Graphics 25 (2001) 13}27
21
Fig. 9. The "rst pair of similar branching nodes.
Fig. 10. The elements of the approximate skeleton involved in the matching process.
following the most similar arcs, until terminal nodes are reached. Such nodes are said to be similar by skeleton and de"ne correspondences among the polygon boundaries. With reference to Fig. 10, the pseudo-code of the matching algorithm is the following:
S "(N , E ) and S "(N , E ) skeletons ? ? ? @ @ @ A" set of correspondences ∀n 3N with degree(n )"3 ? ? ? ∀n 3N with degree(n )"3 @ @ @ if Sim (n , n ) is maximum L ? @ + p Qn ? ? p Qn @ @ , let (a , b ) (a , b ) and (a , b ) be the combination which maximizes Sim (n , n ) L ? @
arc
with a , a , a 3E arcs incident in p ? ? and b , b , b 3E arcs incident in p @ @ match(p , a , p , b ); ? @ match(p , a , p , b ); ? @ match(p , a , p , b ); ? @ The Function match is implemented as follows: match(p , a , p , b ) ? @ + let q "target node of a ? let q "target node of b @ if degree(q )"degree(q )"3 ? @ + let a , a 3E be the other two arcs incident ? in q ? and b , b 3E the other two arcs incident @ in q @ if Sim (a , b );Sim (a , b )'Sim (a , b ); ? ? ? Sim (a , b ) ?
22
M. Mortara, M. Spagnuolo / Computers & Graphics 25 (2001) 13}27
+
3.2. Local similarity match(q , a , q , b ); ? @ match(q , a , q , b ); ? @
, else + match(q , a , q , b ); ? @ match(q , a , q , b ); ? @ ,
, else if degree(q )"degree(q )"1 ? @ AQ(q , q ) ? @ else if degree(q )"1 and degree(q )"3 ? @ + match(p , a , q , b ); ? @ match(p , a , q , b ); ? @ , else if degree(q )"3 and degree(q )"1 ? @ + match(p , b , q , a ); @ ? match(p , b , q , a );, @ ? , With regard to the computational complexity of the similarity evaluation process, if N is the maximum of the number of branching nodes in the two skeletons, then the computation of the "rst couple of similar branching nodes takes at most O(N). Then, in the worst case, visiting the most similar arcs in the two skeletons until similar terminal nodes are reached implies to make as many steps as the number of nodes of the skeletons, which is can be expressed as a linear function of the number of branching nodes. Therefore, the complexity of the whole matching algorithm requires O(N) time.
The correspondences found using the skeleton matching approach divide the shapes, and particularly the boundaries, into correspondent zones: inside each pair of corresponding zones, further correspondences can be set between vertices by local similarity criteria, based on the geometric characteristics of each single vertex. Local similarity criterions are based on proximity information, concavity or convexity, and curvature at each vertex [15]. Proximity is used in the sense that two vertices on two di!erent polygons are judged similar if they are `quite closea (a similar approach has been adopted in [9]). More precisely, we de"ne similarity by neighbourhood as follows: P+ Q 0 d(P, Q)"min+d(P, Q ) " Q 3C ,#d , L C
(3)
where C , C are the two polygons, P3C and Q, Q 3C , d(P, Q) is the Euclidean distance between P and Q and d is an optional user de"ned threshold C (see Fig. 11). Two points on two di!erent polygons are similar by concavity/convexity if they are both concave or convex. Obviously, this criterion is too general to point out useful correspondences, but it is however valid to "nd points which can not be similar, reducing the research area. Finally, if the concept of curvature at a vertex is approximated with the radius of the circle passing for the vertex P, its predecessor and its successor along the boundary, we de"ne similarity by curvature as follows: P+ B 0 "R(P , P, P )!R(B , B, B )")R , (4) ASPT NPCA QCA NPCA QCA C
Fig. 11. Delaunay triangulation and the ASs computed for two polygons; correspondences de"ned by skeleton similarity are depicted in red and those similar by neighbourhood are dashed; correspondences between points similar by curvature are detected among points similar by neighbourhood.
M. Mortara, M. Spagnuolo / Computers & Graphics 25 (2001) 13}27
23
Fig. 12. Two shapes, the structural and local correspondences found, and the intermediate shape.
Fig. 13. Intermediate polygons generated between two shapes.
where P, P , P 3C , B, B , B 3C , P , P and NPCA QCA NPCA QCA NPCA QCA B , B are the previous and the following vertices with NPCA QCA respect to P and B on the same polygon, R(P , P , P ) is the radius of the circle connecting P , P , P , and R is C a user-de"ned threshold. An example of the similarity-based matching is shown in Fig. 11, where the correspondences found using the skeleton matching as well as the local criteria are shown.
4. Shape matching for blending Since the contribution of this paper regards the automatic extraction of global and local correspondences, the choice of vertex interpolation schemes is not discussed here, and examples are given of the generation of intermediate shapes using a simple linear interpolation.
24
M. Mortara, M. Spagnuolo / Computers & Graphics 25 (2001) 13}27
Fig. 14. Correspondences between two di!erent polygons and one intermediate shape.
Fig. 15. Intermediate shapes obtained using the approximate skeleton.
First of all, the integration of skeleton-based and local shape matching seems to give satisfactory results also when the shapes are rather di!erent. In Fig. 12, the shape of an apple is matched against the shape of a pear and reasonable correspondences are found. In all the examples, the correspondence problem is solved trying the skeleton similarity at "rst, if the polygons are structurally compatible, which will produce a set of principal attaches. Then, local shape similarity through curvature and neighbourhood criterion is considered to determine local attaches, by processing only the points belonging to the corresponding zones. The intermediate polygons shown in the examples are computed as a linear interpolation between the polygons coordinates, constrained to the correspondences found. This technique has been developed for surface reconstruction, and it assumes that the two polygons are
positioned on two di!erent planes. As previously stated, however, the approach can be easily extended to more general situations by assigning a di!erent elevation to the two polygons and positioning them so that their centres of mass coincide with respect to the X> plane [15]. The intermediate shape is then constructed by intersecting a plane with the segments joining the corresponding points, and the distance assigned to the polygon planes, if not "xed, together with the number of intermediate sections computed determines the `speeda of the transformation (see Fig. 13). In Fig. 14, another example is shown of correspondences found between shapes with very di!erent skeletons. In this situation, only few points are detected using the skeleton, while the local strategy will "nd more points. If the source polygon must be transformed into more than one target polygon, a change in the topology of the
M. Mortara, M. Spagnuolo / Computers & Graphics 25 (2001) 13}27
25
Fig. 16. Shape transformation involving a branching hole, with a di!erent number intermediate sections in (a) and (b).
polygon during the transformation occurs, which can be handled as a branching instance in surface reconstruction. Let us call the source shape the pre-branch polygon and the n target shapes the post-branch polygons. Our approach is to generate a kind of pseudo-polygon between the pre-branch and the post-branch ones, as in [15]. The pseudo-polygon is obtained by computing the approximate skeleton of the non-simply connected polygon obtained by projecting, or scaling and moving, the post-branch shape onto the pre-branch one, and by pruning the AS of all terminal arcs, independently of their size. The resulting pseudo-polygon is made of several cycles connected together, one for each post-branch polygon. Then, vertex correspondences can be extracted for the pre-branch polygon with respect to the external boundary of the computed pseudo-polygon, and each
single cycle of the shaved skeleton will correspond to exactly one post-branch contour (see Fig. 15). The same method can be applied for branching inner components, as depicted in Fig. 16. If at least one of the polygons to blend is non-simply connected, di!erent strategies have to be used. Let us "rst consider the case in which the source polygon is simply connected, while the other has one (or more) hole. This means that the polygon with holes is described by one external contour and several internal contours. In this situation, the contour of the source polygon is matched with the external contour of the target, using the same algorithm discussed in the previous section. Then, the internal contours of the target are projected onto the source polygon (eventually scaling and moving them so that they fall inside the polygon),
26
M. Mortara, M. Spagnuolo / Computers & Graphics 25 (2001) 13}27
Fig. 17. Transformation between polygons having di!erent topology.
and a point is chosen inside each projections, either the barycentre or any point on the polygon skeleton. These points, which are the seeds for creating the holes, are linked with the more signi"cant points of the corresponding internal contour (e.g. points with highest curvature), giving a set of correspondences that will be interpolated to form intermediate holes that are created and grow from the simple polygon towards the "nal shape, as sketched in Fig. 17. The same idea can be followed if we have to blend polygons with a di!erent number of holes. In this case, some holes must be transformed from the source polygon to the target, while others must be created or must disappear. The holes which have a correspondence are one to one warped, following the algorithm previously discussed for matching simply connected polygons; the ones which have to be generated (or removed) are created from their seeds (or collapse to them). The transformation strategy between polygons with di!erent topology is currently being developed.
5. Conclusions and future work In this paper, a top-down approach for "nding corresponding features of polygonal shapes is described. The approximate skeleton is proposed to describe the shape as a decomposition into morphologically meaningful areas. A strategy for matching the approximate skeletons has been described, which gives a solution to the correspondence problem by taking into account the global structure of the shapes. Also, local similarity criteria based on geometric properties of vertices are used to further re"ne the shape matching. The "rst implementation results con"rm the validity of the approach, which is useful for driving the blending both in surface reconstruction or morphing contexts. Another feature of the method presented here is the use of a pruning step which
makes the approximate skeleton less sensitive to small perturbations and noise on the polygon boundary. If polygons have very di!erent structure and/or topology, the shapes may be judged too dissimilar to produce a meaningful segmentation of the boundaries from the skeleton matching. If no structural decomposition is found, intermediate polygons may be generated following the local similarity analysis on the whole boundaries. Moreover, the skeleton matching can be used also for evaluating the rotation which best align the shapes, before a more precise blending is applied. Beside improving the matching strategy, we are currently working on an extension of the shape description so that a geometric signature can be assigned not only to convex-like features (i.e. protrusions), but also to concave areas of the polygon. Therefore, we are extending the approximate skeleton to the `externala shape, using the polygon boundary as a constraint set for the Delaunay triangulation, but without forcing the boundary of the triangulation to be the boundary of the polygon. Therefore, tracing medial paths of triangles belonging to the convex hull but not to the interior of the polygon, a compact description is obtained of both the inner and the external shape of the polygon. Moreover, since the approximate skeleton can be nicely coded in a hierarchical fashion according to the size of the shape components, we are also studying the possibility of extending the graph matching at di!erent scales. Finally, improvements of the topological operations needed to create/delete or merge/split shapes with an arbitrary number of holes are also planned.
Acknowledgements This work has been partially supported by the Programma Nazionale di Ricerche in Antartide, within the Project `3D sea#oor modelling from irregularly distributed data pointsa. Thanks are given to all the people of the Computer Graphics group at IMA, especially to Silvia Biasotti and Bianca Falcidieno, for their helpful comments.
References [1] Cohen-Or D, Levin D, Solomovici A. Three dimensional distance "eld metamorphosis. ACM Transaction of Graphics, 1998;17:116}41. [2] Alexa M, Cohen-Or D, Levin D. As-rigid-as-possible shape interpolation. SIGGRAPH Proceedings, 2000. [3] Fuchs H, Kedem ZM, Uselton SP. Optimal surface reconstruction from planar contours. Communications of the ACM 1977;20(10):693}702. [4] Meyers D. Reconstruction of surfaces from planar contours. PhD dissertation, University of Washington, 1991.
M. Mortara, M. Spagnuolo / Computers & Graphics 25 (2001) 13}27 [5] Sederberg T, Greenwood E. A physically based approach to 2d shape blending. Computer Graphics (Proc. SIGGRAPH) 1992;26:25}34. [6] Basri R, Costa L, Geiger D, Jacobs D. Determining the similarity of deformable shapes. Vision Research 1998;38:2365}85. [7] Asada H, Brady M. The curvature primal sketch. IEEE Transaction on Pattern Analysis and Machine Intelligence 1986;PAMI 8(1):2}14. [8] Hui KC, Yadong-Li X. A feature based shape blending technique for industrial design. Computer Aided Design 1998;30(10):823}34. [9] Barequet G, Shapiro D. Multilevel sensitive reconstruction of polyhedral surfaces from parallel slices. The Visual Computer 2000;16:116}33. [10] Oliva J-M, Perrin M, Coquillart S. 3D reconstruction of complex polyhedral shapes from contours using a sempli"ed generalised Voronoi diagram. In: Rossignac J, Sillion F (guest editors), Eurographics 96, vol. 15. Blackwell Publishers, 1996. [11] Cohen S, Elber G, Bar-Yehuda R. Matching of freeform curves. CAD 1997;29(5):369}78. [12] Shapira M, Rappoport A. Shape blending using the starskeleton representation. IEEE Computer Graphics and Applications 1995;15:44}50. [13] Spagnuolo M. Shape-based Reconstruction of Natural Surfaces: an Application to the Antarctic Sea Floor. In: Craglia M, Couclelis H (editors). Geographic Information Research: Bridging the Atlantic. London: Taylor and Francis, 1997. [14] Raviola A, Spagnuolo M. Shape-based Surface Reconstruction from pro"les for rapid/virtual prototyping. International Conference on Numerisation3D, Paris, May 1999.
27
[15] Biasotti S, Mortara M, Spagnuolo M. Surface compression and reconstruction using Reeb graphs and shape analysis. Proceedings of the Spring Conference on Computer Graphics 2000, Bratislava, May 2000. [16] Ogniewicz RL. Discrete Voronoi skeletons. PhD thesis, Swiss Federal Institute of Technology, ZuK rich, 1992. [17] Blum H. A transformation for extracting new descriptors of shape. In: Whaten-Dunn W, Editor. Proceedings of the Symposium Models for perception of Speech and Visual form. Cambridge, MA: MIT Press, 1967. 362}80. [18] Lee DT. Medial axis transformation of a planar shape. IEEE Transactions On Pattern Analysis and Machine Intelligence 1982;PAMI-4(4):363}9. [19] Aumann G, Ebner H, Tang L. Automatic derivation of skeleton lines from digitized contours. ISPRS Journal of Photogrammetry and Remote Sensing 1991;46:259}68. [20] De Martino M, Ferrino M. An example of automated shape analysis to solve human perception problems in anthropology. International Journal of Shape Modeling 1996;2(1):69}84. [21] De Floriani L, Falcidieno B, Pienovi C. Delaunay-based representation of surfaces de"ned over arbitrarily shaped domains. Computer Vision, Graphics and Image Processing 1985;32. [22] Dill AR, Levine MD, Noble PB. Multiple resolution skeletons. IEEE Transactions on Pattern Analysis and Machine Intelligence 1987;PAMI-9(4):495}504. [23] Mortara M, Spagnuolo M. Hierarchical representation of 2D polygons based on approximate skeletons, IMA, Technical Report No. 8/00, 2000. [24] Pelillo A, Siddiqi K, Zucker SW. Matching hierarchical structures using association graphs. IEEE Transactions on Pattern Analysis and Machine Intelligence 1999;PAMI21(11):1105}19.
Shape Blending
Similarity measures for blending polygonal shapes Michela Mortara, Michela Spagnuolo* Istituto per la Matematica Applicata, Consiglio Nazionale delle Ricerche, Via De Marini 6, 16149 Genova, Italy
Abstract This paper presents a solution to the correspondence problem for polygon blending, based on a morphological characterization of the shape. The approximate skeleton is suggested for coding the shape as a graph of meaningful areas, and the correspondence problem is therefore solved as a skeleton matching process. The correspondences are found taking into account the global structure of the shapes, and also local similarity criteria are used to further re"ne the matching. The approximated skeleton resembles the medial axis, but it is much simpler to extract and may be easily pruned to reduce the e!ects of small scale details of the shapes' boundary. The "rst implementation results con"rm the validity of the approach for driving the blending both in surface reconstruction or morphing contexts. 2001 Elsevier Science Ltd. All rights reserved. Keywords: Skeleton; Shape similarity; Similarity measures; Feature recognition; Polygon blending
1. Introduction The problem of smoothly transforming, or blending, one shape into another is an intriguing problem with several di!erent applications, ranging from the entertainment industry to engineering applications. In particular, object-space blending deals with the transformation of an object representation into another one. When shapes are modelled as polygons or polyhedra, the blending process is generally seen as composed of two basic steps: "rstly, correspondences among vertices of the source and target models are set, and then an interpolation function is de"ned to produce as many intermediate shapes as required by the application. One of the major problems is to produce intuitive and pleasant intermediate shapes, with a blending which naturally transforms features of the source into similar features of the target shape, avoiding self-intersections of the boundaries or sudden shape changes in intermediate steps [1,2]. Moreover, the blending should be accomplished with minimal interaction of the user.
* Corresponding author. Tel.: #39-01-0647-5677; fax: #3901-0647-5660. E-mail address: [email protected] (M. Spagnuolo).
In a slightly di!erent sense, blending polygonal shapes is necessary for reconstructing surfaces from cross sections. In this case, cross-sections are considered pair-wise to de"ne a tiling among them, that is, a triangular mesh whose edges join vertices belonging to the considered sections. The result of the reconstruction obviously depends on the sampling density, that is, the distance between sections. Since many surfaces could give rise to the observed cross sections, an optimisation is performed with respect to some objective function in order to choose the best surface among the possible ones. Many di!erent goal functions have been used, mainly based on local criteria of the geometric entities involved [3,4]. Blending and tiling polygonal shapes, therefore, exhibit similar problems especially when the shapes to blend are rather di!erent, either in a topological sense (e.g. polygons with di!erent genus) or in a morphological sense. This situation occurs quite frequently in morphing, while in surface reconstruction it generally corresponds to a severely underconstrained tiling problem. In both cases, it is important to have methods for automatically setting correspondences among features of the source and target polygons. Several solutions have been proposed for this problem. For example, physically based methods are aimed at minimising some deformation functions, as in [5] where the minimization concerns the work needed to
0097-8493/01/$ - see front matter 2001 Elsevier Science Ltd. All rights reserved. PII: S 0 0 9 7 - 8 4 9 3 ( 0 0 ) 0 0 1 0 4 - 7
14
M. Mortara, M. Spagnuolo / Computers & Graphics 25 (2001) 13}27
stretch and bend edges and angles along the contours, determining one-to-one correspondences between the vertices. A similar concept has been used in [6] for minimising a curvature-based deformation work function. In industrial design context, a technique has been proposed for locating signi"cant vertices along contours, which represent shape features [8]. Then, local similarity criteria are established for vertices and elements, based on vertex location, curvature, incident and emitting angles and element orientation. For the shape characterization of curves, curvature primitives have been proposed in [7], which could lead to a feature-based and multi resolution contour matching. All these methods fail to take into account the global structure of the shapes for the location of correspondent vertices, and also do not consider the inner part of the shape whose transformation highly in#uences the blending result, as pointed out in [2]. For this aim, it is important to de"ne suitable descriptions of the shape, which may be used to characterize both the inner and boundary of polygons in a computationally usable form. Some works following this principle use structures like modi"ed Voronoi diagrams [10] or a global matching between smooth curves [11] for this purpose. In particular, Shapira and Rappoport [12] use a planar graph obtained by joining the star points of each star-shaped polygon in which the original shape can be divided. Their approach interpolates shapes by "rst interpolating between the skeletons and then unfolding the star pieces from the skeleton. The star skeleton is appropriate for the shape blending problem for two reasons: "rst, it represents all points in the shape's interior as well as on its boundary. Second, it de"nes an explicit dependence between interpolated polygon vertices relatively to a common structure. In any case, the method requires the polygons to have compatible skeletons, that is to say, they must be isomorphic. We have initially considered the problem of shape description for the reconstruction of surfaces from pro"les, for digital terrain modelling and reverse engineering applications. In these contexts, surface reconstruction have been accomplished in two steps: "rstly, the characterization of pro"le shapes is performed, using a multiresolution approach to detect global and local curvature features; then, surface characteristic lines are extracted by joining similar pro"le features, and are used as constraints in the mesh "tting step [13,14]. With regard to polygons, we have used a similar method to e!ectively render generalized surfaces represented by a topological graph, which codes the only adjacencies between critical sections of the original surface. Therefore, the polygons to tile are generally rather distant and may be arbitrarily shaped [15]. In this paper, we suggest the use of an approximate skeleton to describe the shape of polygons and to drive the vertex correspondence phase. The idea we would like
to pursue is to "nd a reasonable match among the approximate skeletons of the polygons to blend, possibly even if the skeletons are quite di!erent. The approximate skeleton of the polygons, which resembles the medial axis, is much simpler to extract and less sensible to convexities on the boundary. With respect to the star skeleton representation, moreover, the approximate skeleton induces a shape decomposition into morphologically meaningful components and is really e$cient to compute. By using global shape structures, it is possible to roughly align the shapes according to their relative direction of best similarity, set global correspondences among shape, and further re"ne the matching, at a lower level, by adding correspondences due to local geometric comparison. The reminder of the paper is organized as follows. In Section 2, the approximate skeleton is introduced and compared with the medial axis transformation; a pruning method for discarding irrelevant details in the resulting shape description is also sketched. Then, the algorithm used to match approximate skeletons is described. In Section 4, some examples are shown of vertex correspondences extracted with our method, and simple blending using linear interpolation of intermediate shapes are given. Conclusions and on-going developments are "nally discussed.
2. Shape characterisation: the approximate skeleton Two-dimensional shapes have been intensively studied in many disciplines, ranging from pure mathematics to psychology, with the aim of de"ning high-level descriptions for the recognition, classi"cation and processing of objects. Characterising a shape means constructing a computational description of the most representative features of the shape, and among the various approaches taken to solve this problem the most suitable one for our purpose is the structural decomposition. According to this paradigm, an object shape is coded as a collection of parts, usually of few basic types, along with their spatial relationships. Therefore, structural decomposition may be seen as a global approach to shape description and the comparison among shapes is generally handled as a labelled graph matching process. Generally, graph matching is a di$cult computational problem, which can be simpli"ed by embedding the problem in its geometric domain. The basic shape elements, or features, to be used for the decomposition should show some properties [16]: shape elements should not change under translation, rotation or scaling (invariance); the description should be locally insensitive to modi"cation of the shape occurring far from the current focus (rich local support); and a small perturbation of the shape should produce only a small perturbation of the description (tolerance to noise).
M. Mortara, M. Spagnuolo / Computers & Graphics 25 (2001) 13}27
Unfortunately, it is rather di$cult to de"ne shape descriptors which ful"l all these requirements. 2.1. The medial axis Among the existing techniques for describing polygonal shapes, the medial axis (MA) is generally considered the more elegant and e!ective one. An intuitive de"nition of the skeleton in the continuum was given by Blum [17], who described the skeleton by analogy with a "re front which starts at the boundary of the shape and propagates isotropically towards the interior. The medial axis is de"ned by the locations at which the "re fronts collide. More formally, the medial axis of a shape S is the locus of centres of all maximal discs of S, that is, those discs contained in S which are not contained in any other disc in S. Equivalently, if B(S) is the boundary of S, then the medial axis of B(S) is the set of points in S having at least two nearest neighbours on B(S). The medial axis, together with the radius function, i.e. the distance from each point on the axis to the nearest point on the boundary, de"ne the medial axis transformation (MAT). The power of this representation is that the shape's boundary and its MAT are equivalent and one can be computed from the other, therefore a two-dimensional object is e!ectively transformed into a one-dimensional graph-like structure. The radius function and the graph-like structure of the MAT induce a decomposition of the shape into subparts which correspond to an intuitive notion of shape protrusion and branching. The points on the axis may indeed be classi"ed according to the number of corresponding nearest neighbours on the shape boundary (see Fig. 1): end points have a single contact point with the boundary, normal points have two contact points and branch points have three or more contact points. A branch point indi-
15
cates a split of the shape into diverging protrusions, while end points indicate a local maximum of curvature of the shape boundary. Therefore, the medial axis can be seen as a graph whose nodes are the end and branch points, and whose arcs are de"ned by normal points. If the shape is a polygon, the MAT is a tree-like planar graph whose arcs are composed by straight-line segments and portions of parabolic curves, and algorithms have been de"ned to compute the MAT from the Voronoi diagram of the polygon elements in O(n log n) time, where n is the number of edges in the polygon [18]. Finally, if the shape is simply connected, then its medial axis graph is a tree, while cycles will appear in the MAT around each hole in the shape. 2.2. The approximate skeleton The idea of approximate skeleton has been used in the "eld of surface reconstruction as a tool for localizing characteristic lines of the surface from a set of its contours [19,20]. Intuitively, an approximate skeleton is generated using the constrained Delaunay triangulation of the polygon as support structure and then by tracing medial segments of the inner triangles, according to rules which take into account the number of constrained edges of each triangle. In this way, a kind of skeleton is produced which globally has the same descriptive e!ect of the medial axis, in the sense that it induces a polygon decomposition in terms of prutrusion-like features. Moreover, the extraction of the approximate skeleton is computationally e$cient. More precisely, let S be a simply-connected polygonal shape having as boundary the polygon P"+P " i"1,2, n, P 3R,; let also ¹ (P) be the DeG G " launay triangulation of S constrained to its boundary
Fig. 1. The medial axis of two polygonal shapes.
16
M. Mortara, M. Spagnuolo / Computers & Graphics 25 (2001) 13}27
Fig. 2. The Delaunay triangulation of a polygon, the triangles having two constrained edges are depicted in dark grey, and the triangles with three inner edges are coloured in light grey.
Fig. 3. The delineation of the AS: some steps of the construction (a) and the "nal AS (b).
[21]. The triangles of ¹ (P) may have either one or two " constrained edges, i.e. edges lying on the boundary P, or may have three non-constrained internal edges (see Fig. 2). Triangles with three inner edges de"ne branch sites of the skeleton, triangles with two edges on the boundary identify terminal nodes of the skeleton, and "nally triangles having only one edge on the boundary contribute to an arc of the skeleton. The AS can be easily traced using medial lines, starting either from branch sites or terminal nodes as follows. With reference to Fig. 3, where the polygon P and its constrained Delaunay triangulation are depicted, let us start the process with any triangle ¹ identifying a terminal node of AS. If v is the vertex of G ¹ shared by the two boundary edges e and e , then the G "rst part of the AS arc ending in v is de"ned by (v, p) where p is the middle point of the opposite edge e . Now, the delineation of the AS continues in the adjacent triangle ¹ according to its type: H E ¹ has no boundary edges: a branching in the AS occurs, H which causes a branching node and two new arcs being
de"ned; the branching node is represented by the centre of mass of ¹ , c, which is joined to the middle H points of its three edges (points q and r in Fig. 3); the tracing continues recursively on the two newly generated arcs, (c, p) and (c, r). E ¹ has 1 boundary edge: the current arc goes on, by H joining the previous one with the middle point of the next non-boundary edge as, for example, (r, t) in Fig. 3. E ¹ has 2 boundary edges: this condition implies that the H current AS arc is terminated, by drawing the segment joining the current arc portion to the vertex common to the two boundary edges. Therefore, the AS can be formally de"ned as follows. De5nition 1. Given a polygonal shape S, with boundary de"ned by the polygon P"+P " i"1,2, n, P 3R,, G G and given the Delaunay triangulation of S constrained to its boundary ¹ (P), the approximate skeleton of S " is the graph AS(S)"(N, A) de"ned by the node set N"+¹ " ¹ 3¹ (P) and ¹ has either no constrained G G " G edges or two constrained edges, and by the arc set
M. Mortara, M. Spagnuolo / Computers & Graphics 25 (2001) 13}27
17
Fig. 4. The MAT of a polygon (a) and the approximate skeleton (b) of the same polygon.
A"+(¹ , ¹ ) " there is a medial path in ¹ (P) between G H " ¹ and ¹ ,. G H The algorithm for extracting AS has to be slightly modi"ed if we have to compute the skeleton of a nonsimply connected polygon. In this case, a starting triangle with two boundary edges may not even exist. If a triangle exists with no edges on the boundary, then the process can start by generating a branching node with degree 3 in the skeleton, and by proceeding recursively on the three generated directions; otherwise, the construction starts from any triangle with one edge on the boundary, with the creation of a dummy node with degree 2, which will be deleted once the skeleton delineation is completed. The computational complexity of the algorithm for the AS extraction is optimal: the choice of the start triangle requires at most ¹ tests, where ¹ is the number of triangles; then, the graph is constructed performing simple operations for each visited triangle, and each triangle is visited exactly one time. Therefore, the algorithm requires O(¹) time. 2.3. Shape description using the approximate skeleton Let us now discuss on the use of the approximate skeleton for describing the shape of a polygon. First of all, the AS is a planar graph having at most degree 3, and this property will be used in the similarity matching step. Compared to the medial axis, the AS is combinatorially simpler in the sense that it generally has fewer arcs and nodes than the MAT. Also, the geometric arcs of the AS are made of line segments and therefore queries like edge length are very easy to implement and fast to compute. In
Fig. 4, an example is shown of the medial axis and the approximate skeleton computed on the same polygon. The decomposition induced by the AS nicely identi"es branching and protrusions of the polygon shape, similarly to that induced by the medial axis. With reference to Fig. 5(a), each arc ending in a terminal node identi"es a protrusion-like feature of the shape (see the portions of shape coloured in grey); each arc between two branching nodes de"nes an inner part of the shape. Therefore, the AS can be e!ectively used to describe the shape of a polygon by attaching to its elements appropriate measures of the components they represent. More precisely, each arc is linked to the area of its component, to the length of the corresponding medial path and to an `average orientationa of the component de"ned as the segment joining the two nodes de"ning the arc (see Fig. 5(b)). Moreover, the decomposition induced by the AS can be easily organised hierarchically according to the scale of the resulting components, following a technique similar to the one used for the medial axis [16,22,23]. Obviously, the approximate skeleton is invariant with respect to shape rotation, and provides a rich local support. However, the MAT and the AS behave similarly with respect to small ondulations or perturbations of the boundary, that is, small modi"cations of the boundary may produce a change in the skeleton structure. This e!ect is shown in Fig. 6 with a simple example. Even if this behaviour is negative for the sake of shape description, the metrics used for characterizing shape components can be used to discard irrelevant details kept in the skeleton. It is straightforward, indeed, to de"ne algorithms for pruning the AS so that only signi"cant arcs are
18
M. Mortara, M. Spagnuolo / Computers & Graphics 25 (2001) 13}27
Fig. 5. The decomposition induced by the AS: protrusions are coloured in light blue and inner components in light orange (a); the dashed segments represents the average orientation associated to arcs (b).
Fig. 6. Small changes in the polygon boundary a!ect the skeleton structure: the original polygon, the medial axis and the approximate skeleton, and the pruned approximate skeleton.
kept in the shape description [23]. An example of this pruning method is depicted in Fig. 6. Intuitively, at the "rst step, protrusions composed by a single triangle with relative small area (depending on the scale and on the average triangle area) are eliminated, and their parent branching node disappears in the AS. The process can be iterated and the size of arcs to prune can be tuned to match the required tolerance. In Fig. 7, the pruning results are shown, where the same shape is represented with a di!erent sampling density, which produces small modi"cations in the AS. As the number of points increases on the boundary, the number of small and irrelevant arcs increases as well, and the pruning algorithm is used to shave away non-signi"cant arcs. In conclusion, the approximate skeleton provides a description of the polygon shape which represents not only boundary features but also takes into account the inner of the polygon.
3. Shape similarity evaluation: structure and geometry Having de"ned a method for describing a polygonal shape, the idea is to establish a set of correspondences between the polygon's boundaries by reasoning on their structural similarity. Eventually, this process will support blending between a source and a target polygon, so that the shape features of the source polygon evolve, if possible, in corresponding similar features of the target polygon. First of all, given two polygons P and P and their respective approximate skeletons AS and AS , it is quite natural to consider that P and P are structurally equiv alent if AS and AS are isomorphic as graphs. There fore, "nding correspondences among shape components is related to the problem of establishing not only a global isomorphism, but also partial matchings among subgraphs. If the problem of establishing isomorphisms among graphs is a complex problem, subgraph matching is even more complicated to be handled in the general
M. Mortara, M. Spagnuolo / Computers & Graphics 25 (2001) 13}27
19
Fig. 7. In the "rst column, the same shape is coded with di!erent sampling; in the second, the computed approximate skeletons and their pruning in the third column.
case [24]. Here, we present an heuristics which strongly uses the geometric embedding of the graphs. For reasoning on skeleton similarity, however, we require that P and P are at least skeleton-comparable, that is, they have at least one node of degree 3; otherwise, the correspondence among features is established using a simpler point-wise approach, as described in the following subsections. 3.1. Skeleton similarity The shape similarity among polygons is evaluated using the following heuristics. First of all, we observe that the branching sites of the approximate skeleton can be seen as the kernel components of the shape, and therefore the comparison should primarily consider branching nodes. A branching node in#uences the shape of a polygon according to the `weighta of the feature components starting at that node, features which are characterized by a geometric signature. Therefore, we start with the de"nition of a measure of similarity for arcs, based on a combination of geometric quantities which take into account the size and the relative orientation of the components. Then, the similarity of branching nodes is de"ned in terms of the combination of their incident arcs which maximizes the sum of the related similarity contributes. Finally, the shapes are compared by searching the &&most similara branching nodes in the two skeletons, and then visiting the graphs from these matched nodes, following each pair of similar arcs until nodes belonging to the
boundary are reached. In this way, correspondences between nodes with degree one, those lying on the boundary, can be established. More precisely, the similarity between arcs is evaluated as a composition of some geometric properties, as follows. De5nition 2 (Similarity between arcs). Let e and e be two arcs of two approximate skeletons; the similarity between arcs, Sim (e , e ), is de"ned as ? Sim (e , e )"aS (e , e )#bS (e , e ) ? ?LEJC JCLERF #dS (e , e )#cS (e , e ), (1) ?PC? RWNC where E S (e , e )"(1803!h)/1803, where h is the angle ?LEJC between the average orientation of e and e . E S (e , e )"l /l if l (l , l /l otherwise, where JCLERF l and l are the length of medial paths corresponding to e and e . E S (e , e )"A /A if A (A , A /A otherwise, ?PC? where A and A are the areas of the two shape components corresponding to e and e. E S (e , e )"1 i! the vertices of e and e have the RWNC same degree, 0 otherwise (see Fig. 8). The parameters a b d c are weights which can be changed according to the requirements of the speci"c input. In particular, it is interesting to consider the e!ect of setting to zero the orientation weight, a. In surface
20
M. Mortara, M. Spagnuolo / Computers & Graphics 25 (2001) 13}27
Fig. 8. Entities used to evaluate the similarity among arcs of the approximate skeleton: in (a) and (b) the average orientation and the medial paths are depicted; in (c) the orientation di!erence between the protrusion alignment.
reconstruction problem, for example, the position of the polygons is "xed and therefore it is reasonable to consider the relative alignment of the shape components to judge their similarity. Conversely, in morphing, similar features should be matched even if the shape is rotated about the branching node. In this latter case, it is useful to discard the contribute due to the relative orientation. These aspects will be discussed also in Section 4. The similarity between branching nodes is computed by evaluating the similarity between each pair of arcs incident in the nodes. De5nition 3 (Similarity of branching nodes). Let A and B be two nodes of degree three of two approximate skeletons; let a , a , a (resp. b , b , b ) be the arcs inci dent in A (resp. B) ordered around A (resp. B). The similarity between A and B, Sim (A, B) is de"ned as L Sim (A, B)"max+(Sim (a , b )#Sim (a , b ) L ? ? #Sim (a , b )), (Sim (a , b ) ? ? #Sim (a , b )#Sim (a , b )), ? ? (Sim (a , b )#Sim (a , b )#Sim (a , b )), ? ? ? (Sim (a , b )#Sim (a , b )#Sim (a , b )), ? ? ? (Sim (a , b )#Sim (a , b )#Sim (a , b )), ? ? ? (Sim (a , b )#Sim (a , b ) ? ? # Sim (a , b )),/distance(A, B). (2) ? Thus, similarity is evaluated by taking the maximum of the sum of the values of incident arc similarity for all the
possible combinations of arcs couples, taken consistently with the ordering around their node. Again, the approach to assess the similarity should take into account the context in which it is performed. If the polygon position is "xed, than the match has to be pursued without allowing any rotation of the shapes. In other words, the orientation of the incident arcs has to be considered and has to be related to a common reference system, in order to avoid torsions in the reconstructed surface. Conversely, if the polygons can rotate during the transformation, the best matching pair can be searched in a slightly di!erent way: for each permutation, the "rst arc couple is aligned, which de"nes a rotation of one shape with respect to the other; then, the relative alignment of the other arcs is evaluated in the rotated reference frame. The parameter distance(A, B) represents the Euclidean distance between points A and B in the plane. If the polygons' position is not "xed a priori, imagine that you move the polygons till their centres of mass coincide, and then the distance between A and B can be computed. The "rst step of the shape comparison process consists in "nding a starting pair of similar branching nodes of the two polygons. Since this pair will be used as the starting point for visiting the two graphs in parallel, it is reasonable to search for the most similar branching nodes in the two shapes. Therefore, the branching nodes which maximise (2) will be chosen among all the possible node combinations. In Fig. 9, an example of this "rst step is shown. Having found the "rst pair of similar branching nodes and having established, at the same time, the three pairs of similar incident arcs, the graphs are recursively visited
M. Mortara, M. Spagnuolo / Computers & Graphics 25 (2001) 13}27
21
Fig. 9. The "rst pair of similar branching nodes.
Fig. 10. The elements of the approximate skeleton involved in the matching process.
following the most similar arcs, until terminal nodes are reached. Such nodes are said to be similar by skeleton and de"ne correspondences among the polygon boundaries. With reference to Fig. 10, the pseudo-code of the matching algorithm is the following:
S "(N , E ) and S "(N , E ) skeletons ? ? ? @ @ @ A" set of correspondences ∀n 3N with degree(n )"3 ? ? ? ∀n 3N with degree(n )"3 @ @ @ if Sim (n , n ) is maximum L ? @ + p Qn ? ? p Qn @ @ , let (a , b ) (a , b ) and (a , b ) be the combination which maximizes Sim (n , n ) L ? @
arc
with a , a , a 3E arcs incident in p ? ? and b , b , b 3E arcs incident in p @ @ match(p , a , p , b ); ? @ match(p , a , p , b ); ? @ match(p , a , p , b ); ? @ The Function match is implemented as follows: match(p , a , p , b ) ? @ + let q "target node of a ? let q "target node of b @ if degree(q )"degree(q )"3 ? @ + let a , a 3E be the other two arcs incident ? in q ? and b , b 3E the other two arcs incident @ in q @ if Sim (a , b );Sim (a , b )'Sim (a , b ); ? ? ? Sim (a , b ) ?
22
M. Mortara, M. Spagnuolo / Computers & Graphics 25 (2001) 13}27
+
3.2. Local similarity match(q , a , q , b ); ? @ match(q , a , q , b ); ? @
, else + match(q , a , q , b ); ? @ match(q , a , q , b ); ? @ ,
, else if degree(q )"degree(q )"1 ? @ AQ(q , q ) ? @ else if degree(q )"1 and degree(q )"3 ? @ + match(p , a , q , b ); ? @ match(p , a , q , b ); ? @ , else if degree(q )"3 and degree(q )"1 ? @ + match(p , b , q , a ); @ ? match(p , b , q , a );, @ ? , With regard to the computational complexity of the similarity evaluation process, if N is the maximum of the number of branching nodes in the two skeletons, then the computation of the "rst couple of similar branching nodes takes at most O(N). Then, in the worst case, visiting the most similar arcs in the two skeletons until similar terminal nodes are reached implies to make as many steps as the number of nodes of the skeletons, which is can be expressed as a linear function of the number of branching nodes. Therefore, the complexity of the whole matching algorithm requires O(N) time.
The correspondences found using the skeleton matching approach divide the shapes, and particularly the boundaries, into correspondent zones: inside each pair of corresponding zones, further correspondences can be set between vertices by local similarity criteria, based on the geometric characteristics of each single vertex. Local similarity criterions are based on proximity information, concavity or convexity, and curvature at each vertex [15]. Proximity is used in the sense that two vertices on two di!erent polygons are judged similar if they are `quite closea (a similar approach has been adopted in [9]). More precisely, we de"ne similarity by neighbourhood as follows: P+ Q 0 d(P, Q)"min+d(P, Q ) " Q 3C ,#d , L C
(3)
where C , C are the two polygons, P3C and Q, Q 3C , d(P, Q) is the Euclidean distance between P and Q and d is an optional user de"ned threshold C (see Fig. 11). Two points on two di!erent polygons are similar by concavity/convexity if they are both concave or convex. Obviously, this criterion is too general to point out useful correspondences, but it is however valid to "nd points which can not be similar, reducing the research area. Finally, if the concept of curvature at a vertex is approximated with the radius of the circle passing for the vertex P, its predecessor and its successor along the boundary, we de"ne similarity by curvature as follows: P+ B 0 "R(P , P, P )!R(B , B, B )")R , (4) ASPT NPCA QCA NPCA QCA C
Fig. 11. Delaunay triangulation and the ASs computed for two polygons; correspondences de"ned by skeleton similarity are depicted in red and those similar by neighbourhood are dashed; correspondences between points similar by curvature are detected among points similar by neighbourhood.
M. Mortara, M. Spagnuolo / Computers & Graphics 25 (2001) 13}27
23
Fig. 12. Two shapes, the structural and local correspondences found, and the intermediate shape.
Fig. 13. Intermediate polygons generated between two shapes.
where P, P , P 3C , B, B , B 3C , P , P and NPCA QCA NPCA QCA NPCA QCA B , B are the previous and the following vertices with NPCA QCA respect to P and B on the same polygon, R(P , P , P ) is the radius of the circle connecting P , P , P , and R is C a user-de"ned threshold. An example of the similarity-based matching is shown in Fig. 11, where the correspondences found using the skeleton matching as well as the local criteria are shown.
4. Shape matching for blending Since the contribution of this paper regards the automatic extraction of global and local correspondences, the choice of vertex interpolation schemes is not discussed here, and examples are given of the generation of intermediate shapes using a simple linear interpolation.
24
M. Mortara, M. Spagnuolo / Computers & Graphics 25 (2001) 13}27
Fig. 14. Correspondences between two di!erent polygons and one intermediate shape.
Fig. 15. Intermediate shapes obtained using the approximate skeleton.
First of all, the integration of skeleton-based and local shape matching seems to give satisfactory results also when the shapes are rather di!erent. In Fig. 12, the shape of an apple is matched against the shape of a pear and reasonable correspondences are found. In all the examples, the correspondence problem is solved trying the skeleton similarity at "rst, if the polygons are structurally compatible, which will produce a set of principal attaches. Then, local shape similarity through curvature and neighbourhood criterion is considered to determine local attaches, by processing only the points belonging to the corresponding zones. The intermediate polygons shown in the examples are computed as a linear interpolation between the polygons coordinates, constrained to the correspondences found. This technique has been developed for surface reconstruction, and it assumes that the two polygons are
positioned on two di!erent planes. As previously stated, however, the approach can be easily extended to more general situations by assigning a di!erent elevation to the two polygons and positioning them so that their centres of mass coincide with respect to the X> plane [15]. The intermediate shape is then constructed by intersecting a plane with the segments joining the corresponding points, and the distance assigned to the polygon planes, if not "xed, together with the number of intermediate sections computed determines the `speeda of the transformation (see Fig. 13). In Fig. 14, another example is shown of correspondences found between shapes with very di!erent skeletons. In this situation, only few points are detected using the skeleton, while the local strategy will "nd more points. If the source polygon must be transformed into more than one target polygon, a change in the topology of the
M. Mortara, M. Spagnuolo / Computers & Graphics 25 (2001) 13}27
25
Fig. 16. Shape transformation involving a branching hole, with a di!erent number intermediate sections in (a) and (b).
polygon during the transformation occurs, which can be handled as a branching instance in surface reconstruction. Let us call the source shape the pre-branch polygon and the n target shapes the post-branch polygons. Our approach is to generate a kind of pseudo-polygon between the pre-branch and the post-branch ones, as in [15]. The pseudo-polygon is obtained by computing the approximate skeleton of the non-simply connected polygon obtained by projecting, or scaling and moving, the post-branch shape onto the pre-branch one, and by pruning the AS of all terminal arcs, independently of their size. The resulting pseudo-polygon is made of several cycles connected together, one for each post-branch polygon. Then, vertex correspondences can be extracted for the pre-branch polygon with respect to the external boundary of the computed pseudo-polygon, and each
single cycle of the shaved skeleton will correspond to exactly one post-branch contour (see Fig. 15). The same method can be applied for branching inner components, as depicted in Fig. 16. If at least one of the polygons to blend is non-simply connected, di!erent strategies have to be used. Let us "rst consider the case in which the source polygon is simply connected, while the other has one (or more) hole. This means that the polygon with holes is described by one external contour and several internal contours. In this situation, the contour of the source polygon is matched with the external contour of the target, using the same algorithm discussed in the previous section. Then, the internal contours of the target are projected onto the source polygon (eventually scaling and moving them so that they fall inside the polygon),
26
M. Mortara, M. Spagnuolo / Computers & Graphics 25 (2001) 13}27
Fig. 17. Transformation between polygons having di!erent topology.
and a point is chosen inside each projections, either the barycentre or any point on the polygon skeleton. These points, which are the seeds for creating the holes, are linked with the more signi"cant points of the corresponding internal contour (e.g. points with highest curvature), giving a set of correspondences that will be interpolated to form intermediate holes that are created and grow from the simple polygon towards the "nal shape, as sketched in Fig. 17. The same idea can be followed if we have to blend polygons with a di!erent number of holes. In this case, some holes must be transformed from the source polygon to the target, while others must be created or must disappear. The holes which have a correspondence are one to one warped, following the algorithm previously discussed for matching simply connected polygons; the ones which have to be generated (or removed) are created from their seeds (or collapse to them). The transformation strategy between polygons with di!erent topology is currently being developed.
5. Conclusions and future work In this paper, a top-down approach for "nding corresponding features of polygonal shapes is described. The approximate skeleton is proposed to describe the shape as a decomposition into morphologically meaningful areas. A strategy for matching the approximate skeletons has been described, which gives a solution to the correspondence problem by taking into account the global structure of the shapes. Also, local similarity criteria based on geometric properties of vertices are used to further re"ne the shape matching. The "rst implementation results con"rm the validity of the approach, which is useful for driving the blending both in surface reconstruction or morphing contexts. Another feature of the method presented here is the use of a pruning step which
makes the approximate skeleton less sensitive to small perturbations and noise on the polygon boundary. If polygons have very di!erent structure and/or topology, the shapes may be judged too dissimilar to produce a meaningful segmentation of the boundaries from the skeleton matching. If no structural decomposition is found, intermediate polygons may be generated following the local similarity analysis on the whole boundaries. Moreover, the skeleton matching can be used also for evaluating the rotation which best align the shapes, before a more precise blending is applied. Beside improving the matching strategy, we are currently working on an extension of the shape description so that a geometric signature can be assigned not only to convex-like features (i.e. protrusions), but also to concave areas of the polygon. Therefore, we are extending the approximate skeleton to the `externala shape, using the polygon boundary as a constraint set for the Delaunay triangulation, but without forcing the boundary of the triangulation to be the boundary of the polygon. Therefore, tracing medial paths of triangles belonging to the convex hull but not to the interior of the polygon, a compact description is obtained of both the inner and the external shape of the polygon. Moreover, since the approximate skeleton can be nicely coded in a hierarchical fashion according to the size of the shape components, we are also studying the possibility of extending the graph matching at di!erent scales. Finally, improvements of the topological operations needed to create/delete or merge/split shapes with an arbitrary number of holes are also planned.
Acknowledgements This work has been partially supported by the Programma Nazionale di Ricerche in Antartide, within the Project `3D sea#oor modelling from irregularly distributed data pointsa. Thanks are given to all the people of the Computer Graphics group at IMA, especially to Silvia Biasotti and Bianca Falcidieno, for their helpful comments.
References [1] Cohen-Or D, Levin D, Solomovici A. Three dimensional distance "eld metamorphosis. ACM Transaction of Graphics, 1998;17:116}41. [2] Alexa M, Cohen-Or D, Levin D. As-rigid-as-possible shape interpolation. SIGGRAPH Proceedings, 2000. [3] Fuchs H, Kedem ZM, Uselton SP. Optimal surface reconstruction from planar contours. Communications of the ACM 1977;20(10):693}702. [4] Meyers D. Reconstruction of surfaces from planar contours. PhD dissertation, University of Washington, 1991.
M. Mortara, M. Spagnuolo / Computers & Graphics 25 (2001) 13}27 [5] Sederberg T, Greenwood E. A physically based approach to 2d shape blending. Computer Graphics (Proc. SIGGRAPH) 1992;26:25}34. [6] Basri R, Costa L, Geiger D, Jacobs D. Determining the similarity of deformable shapes. Vision Research 1998;38:2365}85. [7] Asada H, Brady M. The curvature primal sketch. IEEE Transaction on Pattern Analysis and Machine Intelligence 1986;PAMI 8(1):2}14. [8] Hui KC, Yadong-Li X. A feature based shape blending technique for industrial design. Computer Aided Design 1998;30(10):823}34. [9] Barequet G, Shapiro D. Multilevel sensitive reconstruction of polyhedral surfaces from parallel slices. The Visual Computer 2000;16:116}33. [10] Oliva J-M, Perrin M, Coquillart S. 3D reconstruction of complex polyhedral shapes from contours using a sempli"ed generalised Voronoi diagram. In: Rossignac J, Sillion F (guest editors), Eurographics 96, vol. 15. Blackwell Publishers, 1996. [11] Cohen S, Elber G, Bar-Yehuda R. Matching of freeform curves. CAD 1997;29(5):369}78. [12] Shapira M, Rappoport A. Shape blending using the starskeleton representation. IEEE Computer Graphics and Applications 1995;15:44}50. [13] Spagnuolo M. Shape-based Reconstruction of Natural Surfaces: an Application to the Antarctic Sea Floor. In: Craglia M, Couclelis H (editors). Geographic Information Research: Bridging the Atlantic. London: Taylor and Francis, 1997. [14] Raviola A, Spagnuolo M. Shape-based Surface Reconstruction from pro"les for rapid/virtual prototyping. International Conference on Numerisation3D, Paris, May 1999.
27
[15] Biasotti S, Mortara M, Spagnuolo M. Surface compression and reconstruction using Reeb graphs and shape analysis. Proceedings of the Spring Conference on Computer Graphics 2000, Bratislava, May 2000. [16] Ogniewicz RL. Discrete Voronoi skeletons. PhD thesis, Swiss Federal Institute of Technology, ZuK rich, 1992. [17] Blum H. A transformation for extracting new descriptors of shape. In: Whaten-Dunn W, Editor. Proceedings of the Symposium Models for perception of Speech and Visual form. Cambridge, MA: MIT Press, 1967. 362}80. [18] Lee DT. Medial axis transformation of a planar shape. IEEE Transactions On Pattern Analysis and Machine Intelligence 1982;PAMI-4(4):363}9. [19] Aumann G, Ebner H, Tang L. Automatic derivation of skeleton lines from digitized contours. ISPRS Journal of Photogrammetry and Remote Sensing 1991;46:259}68. [20] De Martino M, Ferrino M. An example of automated shape analysis to solve human perception problems in anthropology. International Journal of Shape Modeling 1996;2(1):69}84. [21] De Floriani L, Falcidieno B, Pienovi C. Delaunay-based representation of surfaces de"ned over arbitrarily shaped domains. Computer Vision, Graphics and Image Processing 1985;32. [22] Dill AR, Levine MD, Noble PB. Multiple resolution skeletons. IEEE Transactions on Pattern Analysis and Machine Intelligence 1987;PAMI-9(4):495}504. [23] Mortara M, Spagnuolo M. Hierarchical representation of 2D polygons based on approximate skeletons, IMA, Technical Report No. 8/00, 2000. [24] Pelillo A, Siddiqi K, Zucker SW. Matching hierarchical structures using association graphs. IEEE Transactions on Pattern Analysis and Machine Intelligence 1999;PAMI21(11):1105}19.