Providence, April 27, 1998


On 3D Medial "Axes" Computations

Presented by Frederic Leymarie

Brown University

From 1980 -> today

ToC:

Mathematical Morphology : 3D Thinning

Computational Geometry: Polyhedra & Max. Inscribed Spheres

Wave Propagation: Cellular automata & Distance Transforms

Initial assumptions:

In other words: "We have a bounded section of an equidistantly sampled Euclidean space of dimension D=3."


Mathematical Morphology : 3D Thinning

This is also called:

Early work by Lobreght in 1980 who used the Euler number to determine when to remove a voxel :

Sufficiency condition:

Toriwaki showed (1982) that a an additional condition was required: only one "cluster" of voxels (a single and simply connected set) should be attached to the central voxel, to be removed.

Near Euclidean metric

Verwer (1991) introduces "bucket sorting" to minimize voxels' status update. This is done by ordering voxels in terms of distance to the object's surface. It optimizes numerical complexity and permits to use near Euclidean metrics (integer or vector-based).

Formalism for Thinning in (N-D) sampled Euclidean Spaces

Jonker (1995) answers the following questions:

Specificity of Thinning:

"Thinning algorithms operate within an image on the objects, to which connectivity can be associated."

Definition of skeletonization in 3D - notion of Dimensional Reduction :

Connectivity paradox (of squared/cubical grids):

An object (foreground or FG) can pass through the background (BG) when it shouldn't and vice-versa (in 2D this leads to the 4-connected - 8-connected distinction between FG and BG). In 3D, we get the following choices to resolve the "paradox":

Euclidean metric:

Firstly, perform an EDT (Euclidean Distance Transform) to assign a distance value to each object voxel. The EDT is a mapping from a two-valued set (object - background) to a multi-valued set, where distance values (from the object's boundary) are assigned to the background (or to the interior of the object).

Secondly, do the thinning in order of increasing distance.

Characteristics of thinning

Good points:

  1. Simple process (numerically speaking and in terms of implementation).
  2. Attempt at simulating a Grassfire propagation to elicit medial symmetries.
  3. Dimensional reduction is available.

Weak points:

  1. Assumes a compact object/bounding surface (cannot easily deal with partial data).
  2. Choice of connectivity criterion (variability in results).
  3. All "medial" points are treated equally (except end-points --> ribs -->ridges).
  4. No obvious generalization to other symmetry computations (such as Mid-chord loci or the full Symmetry Set).
  5. Fundamentally "discrete" process: no obvious generalization to continuous thinning processes.

Computational Geometry:
Polyhedra & Max. Inscribed Spheres

In 2D, this is equivalent to representing the object by sets of (convex) polygons and deriving an analytic locus of symmetries (e.g. cf. Preparata in 1977). Thus, most of the perspiration is spent on obtaining a valid representation of an object in terms of polygons (in 2D) or polyhedra (in 3D).

This approach is strongly influenced by CSG (Computational Solid Geometry) representations in use within CAD-systems. In other words, the usefulness of the method is driven by applications to "human-made" objects, as modeled today in computer systems. Many references in the recent literature, through the 80's and 90's.


Fig.1 : Sheets, Seams and Rims (from [Sherbrooke96]).

N.B.: Connected components of end points are called rims (our "ribs").

Delaunay triangulation and Voronoï Diagrams

The properties of the Delaunay triangulation of point sets has been exploited as a means of skeletal construction in both 2D and 3D. E.g. Yu et al. (1991) relies on the 3D Delaunay triangulation of a dense distribution of points lying on the boundary of the solid. It is demonstrated that the circumspheres of the resulting tetrahedra approximate maximal spheres.

DEFINITION: Delaunay Triangulation

Given a set S of N distinct points p1, p2, ... pN, in R3, each point pi in S can be associated with a unique polyhedron vi, 1 <= i <= N, where vi is the locus of all points x in R3 that are closer to pi than to any other point of S i.e.,

vi = { x : d(x , pi) < d(x , pj) , for all i not j , pi, pj in S }

where d is the Euclidean distance. vi is an open convex polyhedron sometimes referred to as the Voronoi polyhedron or Voronoi region.

The Delaunay triangulation of S consists of the set of tetrahedra formed by connecting the neighboring points pi and pj between the adjacent polyhedra vi, vj.

A fundamental property of the Delaunay triangulation of S is the empty circumsphere property. The circumsphere associated with each Delaunay tetrahedron does not contain any points of S in its interior. This is a direct consequence of the Voronoi construction since the vertices of the Voronoi diagram are the circumcenters of the associated tetrahedra.


Fig.2 In 2D, the circumcircle of a Delaunay triangle approximates an inscribed circle (from [Sheehy96]).

Given a set of points S lying on the boundary of a solid T and its corresponding Delaunay triangulation, it can be demonstrated that as S becomes infinitely dense the circumsphere associated with each tetrahedron in the interior of T tends to a maximal sphere and its circumcenter to a point on the medial surface.

Generating and triangulating an infinitely dense distribution of points is not possible in practice. From sparse sets however, the approach is practical and leads to the correct topology.


Fig.3 Boundary & MA of a 30-element shell (from [Sherbrooke96]).

Characteristics of Delaunay/Voronoi-based methods

Good points:

  1. It becomes possible to generate "analytic" representations of medial axes.
  2. Medial points can be distinguished.
  3. Chord midpoints can be retrieved.

Weak points:

  1. Assumes a sparse representation of the object/bounding surface (numerical complexity blows up with dense data).
  2. Sparseness may lead to a poor "shape" approximation.
  3. Does not deal with curved objects very well: complexity rises sharply.
  4. No obvious generalization to the full symmetry set elicitation
  5. Fundamentally "global" process: cannot deal easily with open-ends or missing information (e.g. due to occlusion).
  6. Characteristics of medial axes (seems, end-points, rims, sheets, ...) obtained in a post-processing stage.

Wave Propagation :
Cellular automata & Distance Transforms

In an attempt to address problems of previous approaches as well as going further ahead in designing a powerful process for shape description and representation we propose a "new" approach for dealing with symmetry elicitation in 3D.

What we are aiming at:

The solution: Huygen's principle revisited

Let Phi(q0,t) be the wavefront of the point q0 after time t. For every point q of this front, consider the wavefront after time s: Phi(q,s).

Then, the wavefront of the point q0 after time s+t, Phi(q0, s+t), will be the envelope of the fronts Phi(q,s), such that q in Phi(q0,t).


Fig.4 Illustration of Huygens' principle.

In this formulation, one can replace the "point" q0 by a curve, surface, or in general, by a closed set.

Furthermore, the 3D space can be replaced by any smooth manifold, and the propagation of light can be replaced by the propagation of any disturbance transmitting itself "locally" [Arnold89].

We can also replace the metric from Euclidean to a Riemannian, and model anisotropic "shape" properties (of the medium, or space).

The difficulty, is to translate this principle in a discrete world relying upon a lattice, such as the X-Y-Z grid made of voxels, traditionally used in image processing.

We have developed such a "discrete" method relying upon a 3D grid. We call it 3D-CADT for:

In this approach, we simulate propagation of information by the use of a class of cellular automata "oriented" in the direction of the wavefront. Equivalently, these automata rely upon the metric characteristic of the underlying space (the "indicatrix"). For Euclidean 3D space, this is a sphere or directions.

The resultant ray system, is the system of normals to a surface in Euclidean space. In the ngb. of a smooth surface, its normals form a smooth fibration, until at some distance from the surface, various normals may begin to intersect one another.

Historical note: These "singularities" were first investigated by Archimedes (caustics of light rays, around 300 B.C.).

Some examples

1 - A point and a plane generate a paraboloid surface


Fig.5 Initial level sets (d=3) from a sphere and a plane together with
resultant medial axis (a parabola).


Fig.6 Initial level sets (d=10) from a sphere and a plane together with
resultant medial axis (a parabola).

2 - An open cylinder generates a stick axis


Fig.7 Initial level sets (d=5) from a cylinder together with
resultant medial axis (a "stick").

3 - Intersection of a paraboloid & a plane : medial axis with a A1^3 medial curve


Fig.8 "Peter's" object : z1 = (ay+b)x^2 (b>0) , z2 = cx + dy + e (e > 1/2b)


Fig.9 The computed skeleton for Peter's object.

Computation of 3D shocks

3D shocks have been proposed by Kimia.These are based on curvature properties of a level hypersurface, £, in a 4D space which represents the evolution of a 3D object under a wave propagation (or under a motion defined by a so-called "reaction-diffusion" PDE).

Shock type

Gradient

Principal curvatures

1-n

£' disc.

Undefined

2-1

£' = 0

K1.K2 < 0 & K3 != 0

2-2

£' = 0

K1.K2 < 0 & K3 = 0

3-1

£' = 0

K1.K2 > 0 & K3 = 0

3-2

£' = 0

K1=K2=0 & K3!=0

4

£' = 0

K1.K2 > 0 & K2.K3 > 0

Table 1: The classification of 3D shocks based on £


Figure 10: The 4 classes of shocks


Next ....


Some pointers


Page created by Frederic Leymarie, 1998.
Comments, suggestions, etc., mail to: leymarie@lems.brown.edu