Next: 11.1.3.2 Theoretical analysis of Up: 11.1.3 Medial axis Previous: 11.1.3 Medial axis   Contents   Index


11.1.3.1 Algorithms for determining the MAT or related sets

The MAT was introduced and explored by Blum [30] and further explored by Blum [29] and Blum and Nagel [31] to describe biological shape. Soon after it was introduced, various algorithms for computing the MAT were developed for special planar regions. Montanari [272] developed an algorithm to compute the MAT of a multiply-connected polygonal region. His algorithm proceeds by identifying significant branch points and propagating the boundary contour inward, while connecting the branch points with appropriate linear or parabolic segments. A more efficient algorithm for computing the MAT of a convex polygonal region in was presented by Preparata [334] along with an algorithm for a non-convex polygonal region. Lee [229] developed an algorithm for polygonal regions with non-convex corners. Srinivasan and Nackman [405] presented an algorithm for multiply connected polygonal regions with holes. Gursoy and Patrikalakis [140,141,142] developed an algorithm to compute the MAT of a multiply connected planar region bounded by line segments and circular arcs, and used this algorithm to generate finite element meshes automatically and determine global shape characteristics. Guibas and Stolfi [137] investigate the relationship between the Voronoi Diagram and the Delaunay Triangulation, and develop the quad-edge data structure to represent them. Sugihara [416] investigates the use of Voronoi Diagrams to approximate various types of generalized Voronoi Diagrams. Rosenfeld [350] considers different representations of shapes based on an axis and a generation rule. Held's book [157] contains a comprehensive review of Voronoi Diagram algorithms, which he uses in the context of pocket machining. Another comprehensive review of the state of the art in Voronoi Diagram algorithms has been compiled by Aurenhammer [14].

Other work has concentrated on discrete and approximate approaches to determine the MAT or its related sets. Nackman [282] proposes a 3-D algorithm to use a polyhedral approximation of a smooth boundary and produce a polyhedral approximation to the skeleton. The algorithm is an extension of Bookstein's line skeleton approach [36] to 3-D. It takes as input a polyhedral surface made up of convex polygons and generates a connected graph of convex polygons approximating the MA of the original object; since the input polyhedron is assumed to be an approximation to a smooth curved object, the output is not the skeleton of the polyhedron itself but rather a collection of polygons approximately tangent to the skeleton of the true object. Lavender et al. [226] use an octree-based approach to determine the Voronoi Diagram. Their algorithm works on set-theoretic solid models, composed of unions, intersections, and differences of primitive regions represented by a collection of polynomial inequalities, and produces an octree (or quadtree in two dimensions) which divides space into Voronoi regions at some specified resolution. Scott et al. [370] discuss a method for determining the Symmetric Axis (a superset of the MA) which is based on a combined wave/diffusion process in the plane. Their algorithm proceeds by assigning each boundary pixel a unit displacement above the plane and every other pixel a zero displacement, and then numerically propagating a wave from the boundary. The wave is attenuated by a diffusion process to reduce numerical error, and local maxima in the wave are declared to lie on the Symmetric Axis. Although useful for binary images at low resolutions, the error may be large for higher resolutions. Memory and processing requirements for this method tend to be quite high as well.

Brandt, and Brandt and Algazi [41,40,39] find a continuous approximation to the skeleton in both the planar and the 3-D case by first discretizing the boundary. The boundary is sampled at a given sampling density, yielding a set of discrete points which form a pixelized or voxelized approximation to the boundary. The next step is to run an efficient discrete-point Voronoi diagram on the set of points. Finally, portions of the skeleton which result from the effects of quantization are pruned away [39]. This approach attempts to classify each of the vertices in the interior Voronoi diagram according to how many foot points the vertex has. The number of footpoints is determined by taking the associated maximal sphere at the vertex, increasing the radius slightly, and intersecting the dilated sphere with the boundary. This intersection partitions the surface of the sphere into areas which lie either inside or outside the region. Each area lying outside the boundary is assumed to correspond to a footpoint; since the most commonly occurring type of skeleton point has two footpoints, only these points are kept, and the rest are pruned away.

Chiang [54] takes a planar region bounded by piecewise curves and performs a cellular decomposition of the plane in a neighborhood of the region. Each cell is assigned an approximate distance to the nearest point on the boundary of the region using an algorithm due to Danielsson [70] which computes the Euclidean distance transform. In order for us to explain the Euclidean distance transform, we consider a given binary image (where each pixel is assigned 0 or 1); we call the set of pixels with value 0 . The distance map is defined as a scalar function on

    (11.1)

where is a distance function. If is the Euclidean distance between two pixels
    (11.2)

then the distance map is called the Euclidean distance map. This information is later used to find a starting point for tracing axis branches in two dimensions and for recognizing when the tracing has passed the end of a branch. The tracing itself uses the distance information to determine on which boundary elements the footpoints of the current Medial Axis point lie. Once these elements are known, a set of simultaneous equations describing the local structure of the MA near the given point is formed. Using these equations to determine the tangent to the MA at the given point, a short distance along the tangent is traversed, the point is refined with Newton iteration, and another tracing step is taken. At each step, the distance information is used to determine whether or not the current branch has become inactive. If so, a branch point or an end point has been hit, and the tracing either proceeds along another branch or stops. Although the tracing is not extended to three dimensions, Chiang [54] notes that the same Euclidean distance transform in 3-D may be used to determine an approximation to the skeleton. One simple way of using the distance transform in this way is to identify those points which have locally maximal distance values after the distance transform is carried out. Such points are clearly close to centers of maximal disks, and so can be considered to provide an approximation to the skeleton.

Sudhalkar et al. [415] introduce a set called the box skeleton which they argue has the properties which make the MAT desirable as an alternate representation of shape. In particular, the box skeleton exhibits dimensional reduction, homotopic equivalence and invertibility. However, their skeleton is defined using the norm (the box norm) instead of the Euclidean norm and thus may be quite different from the Medial Axis. Their algorithm for determining the box skeleton operates on discrete objects made of unit squares (or cubes, in 3-D) and proceeds by thinning the object while maintaining homotopic equivalence to the original object. In order to perform the thinning, the object is transformed into a graph; in the planar case, the boundaries between adjacent pixels are considered to be edges of the graph, and the intersections of these edges are the vertices of the graph. The first thinning step proceeds by replacing the graph by that portion of the dual graph which is interior to the original (primal) graph. Since this procedure alone may result in a disconnected skeleton, the boundary ``shrink wraps'' around the skeleton as it thins. Procedures based on these concepts are developed for both 2-D and 3-D discrete objects.

Most of the 3-D algorithms in existence (such as the ones above) are fundamentally discrete algorithms. To our knowledge, few continuous approaches have been proposed, due largely to the computational complexity involved. One of the few such techniques is developed by Hoffmann [170], who proposed a method for assembling the skeleton of a CSG object. His method proceeds by determining points of closest approach between pairs of boundary elements and checking these points to make sure they are in fact on the Medial Axis. (Each of these points are on the MA if and only if the distance to the pair of elements is less than or equal to the distances to the other boundary elements.) The points are then sorted in order of increasing distance from the boundary, and then a local analysis around each point is performed, in an attempt to identify whether the point lies on a face, edge, or vertex of the Medial Axis. This determination is made by identifying all boundary elements which lie at the same minimum distance from the point and forming a set of simultaneous equations in variables which describe the equidistantial set that the point belongs to. Based on the rank of the Jacobian of this set of equations, the point is predicted to lie on a face, edge, or vertex of the skeleton. Then neighboring faces and edges are traced out in order of increasing distance from the boundary. The method requires intersection of equidistantial sets with one another in order to trim away portions which do not belong to the Medial Axis. Related papers by Dutta and Hoffmann [81,82] consider the exact representation of the bisectors which appear as skeleton branches in the skeleton of CSG objects bounded by planes, natural quadrics, and torii.

Reddy and Turkiyyah [341] propose an algorithm for determining the skeleton of a 3-D polyhedron based on a generalization of the Voronoi Diagram. They compute an abstract Delaunay Triangulation of the polyhedron and use the result to obtain the dual, the generalized Voronoi Diagram. The Delaunay triangulation computed is a generalization of the usual Delaunay triangulation, which connects isolated nodes together with line segments. In the generalization, the nodes represent parts of a polyhedron (specifically, either a face or a non-convex edge or a non-convex vertex of the polyhedron) and therefore the triangulation is an abstract graph. It still maintains duality with the generalized Voronoi Diagram, and is easier to compute since the Voronoi Diagram may contain trimmed quadric surfaces. The skeleton of the polyhedron is obtained by trimming away certain elements of the generalized Voronoi Diagram which are in the Voronoi Diagram but not the skeleton (for example, the equidistantial point set between a face and a non-convex vertex bounding it). The algorithm can explicitly determine certain critical points of the skeleton, but does not contain highly accurate explicit representations of the curves and surfaces making up the skeleton. Under the same approach, Turkiyyah et al. [428] developed an accelerated algorithm using Delaunay triangulation with a local optimization scheme for the generation of accurate skeletons of 3D solid models. The accelerated algorithm has linear complexity in terms of the number of points used for skeleton approximation.

Sheehy et al. [385,384] investigate the use of a domain Delaunay triangulation on a distribution of points on the boundary (in a manner similar to Brandt [39]) to attempt to determine the topological features of the Medial Axis of a 3-D solid. The steps of an algorithm to compute the Medial Axis from these features are outlined.

Sherbrooke et al. [393] developed an algorithm for determining the MAT of a general polyhedral solid with simply connected faces and without cavities. In [394] and [391], that algorithm was refined, simplified and extended to work for polyhedral solids of arbitrary genus without cavities, with nonconvex vertices and edges. The algorithm is based on a classification scheme that relates different pieces of the medial axis to one another even in the presence of degenerate MA points. Vertices of the MA are connected to one another by tracing along adjacent edges, and finally the faces of the axis are found by traversing closed loops of vertices and edges. The completeness, complexity and stability of the algorithm were also analyzed.

Etzion and Rappoport [91] presented an algorithm for computing the Voronoi Diagram of a 3-D polyhedron based on subdivision of space. This method enables local and partial computation of the Voronoi diagram. By separating the computation of the symbolic (Voronoi graph) and geometric parts of the diagram, the algorithm tends to be more robust. In addition, a tracing algorithm for the MAT of a 3-D polyhedron is developed by Culver et al. [68] to improve accuracy using exact arithmetic and exact geometric representations. In order to improve the efficiency in searching MA elements, spatial decomposition and linear programming are utilized. Sheets in this methods are in the form of quadrics. The method also includes a new algorithm for analysis of the topology of an algebraic plane curve and a fast numerical method for implicit geometric computation.

Wolter and his associates [340,451] compute medial curves on a surface, which is the locus of points being equidistant from two given curves on the surface, utilizing the geodesic offset function. Their method is also applicable to the plane curve case. Also Wolter and his associates [214] applied the above method to compute a Voronoi diagram on a parametric surface instead of the usual Voronoi diagram in Euclidean space.

In addition to the work described above which discusses the determination of the MAT, there has been additional work on the inverse problem of reconstructing the original solid from the Medial Axis and the radius function. Reconstruction of boundary surfaces from curves and surfaces of the Medial Axis in 3-D is discussed by Gelston and Dutta [125]. Vermeer [431] considers the problem of boundary reconstruction from the MAT in 2-D and 3-D. Recently, Verroust and Lazarus [432] developed a method to extract skeletal curves from a set of scattered points that are sampled from a surface using a geodesic graph and distance map with various levels of decomposition. The skeletal curves extracted can be applied in surface reconstruction.

Work on the MAT and its applications has been somewhat limited by the difficulty of developing an algorithm which is robust, accurate, and efficient to carry out the transformation especially for curved solids. Most recent work on the 3-D problem has tried to generate a discrete approximation to the actual MAT; while such algorithms may be well-suited for some meshing applications, particularly if they involve some degree of user interaction, they are less satisfactory for modeling since they do not capture the topological structure of the shape accurately. However, because the skeleton is typically path connected, and because branch points on the skeleton may be expressed as solutions of a set of simultaneous nonlinear polynomial equations, we believe that a continuous approach to the problem is promising especially if combined with interval or even exact arithmetic methods, as for example in the recent work by Culver et al. [68] for 3-D polyhedra.



Next: 11.1.3.2 Theoretical analysis of Up: 11.1.3 Medial axis Previous: 11.1.3 Medial axis   Contents   Index
December 2009