Modeling Multi-Object Configurations via Medial/Skeletal Linking Structures

TL;DR

This paper presents a mathematical framework using skeletal linking structures to model and analyze configurations of multiple objects in images, capturing both individual shapes and their spatial relationships.

Contribution

It introduces a novel skeletal linking structure that extends existing models to handle multiple objects and their interactions in 2D and 3D images.

Findings

Classification of Blum linking structures for generic configurations

Extension of skeletal structures to multi-object configurations

Enhanced modeling of internal shapes and external spatial relations

Abstract

We introduce a method for modeling a configuration of objects in 2D or 3D images using a mathematical "skeletal linking structure" which will simultaneously capture the individual shape features of the objects and their positional information relative to one another. The objects may either have smooth boundaries and be disjoint from the others or share common portions of their boundaries with other objects in a piecewise smooth manner. These structures include a special class of "Blum medial linking structures," which are intrinsically associated to the configuration and build upon the Blum medial axes of the individual objects. We give a classification of the properties of Blum linking structures for generic configurations. The skeletal linking structures add increased flexibility for modeling configurations of objects by relaxing the Blum conditions and they extend in a minimal way the individual "skeletal structures" which have been previously used for modeling individual objects and capturing their geometric properties. This allows for the mathematical methods introduced for single objects to be significantly extended to the entire configuration of objects. These methods not only capture the internal shape structures of the individual objects but also the external structure of the neighboring regions of the objects.

Tables

Table 1. Table 1. Classification of Linking Types for Blum Medial Linking Structures in ℝ 3 superscript ℝ 3 {\mathbb{R}}^{3}

	Linking Type	Dimension	Description of Linking
	$A_1^2$ Linking
i)	$\left(A_1^2:A_1^2,A_1^2\right)$	2	between $2$ smooth points
ii)	$\left(A_1^2:A_1^3,A_1^2\right)$	1	between a smooth point
			and a Y-junction point
iii)	$\left(A_1^2:A_3,A_1^2\right)$	1	between a smooth point
			and an edge point
iv)	$\left(A_1^2:A_3{A}_1,A_1^2\right)$	0	between a fin point
			and a smooth point
v)	$\left(A_1^2:A_1{A}_3,A_1^2\right)$	0	between a smooth point
			associated to a fin point and
			another smooth point
vi)	$\left(A_1^2:A_1^4,A_1^2\right)$	0	between a smooth point
			and a $6$-junction point
vii)	$\left(A_1^2:A_1^3,A_1^3\right)$	0	between $2$ Y-junction points
viii)	$\left(A_1^2:A_3,A_3\right)$	0	between $2$ edge points
ix)	$\left(A_1^2:A_1^3,A_3\right)$	0	between a Y-junction point
			and an edge point
	$A_1^3$, $A_1^4$ and $A_1{A}_3$ Linking
x)	$\left(A_1^3:A_1^2,A_1^2,A_1^2\right)$	1	between $3$ smooth points
xi)	$\left(A_1^3:A_1^3,A_1^2,A_1^2\right)$	0	between $2$ smooth points
			and a Y-junction point
xii)	$\left(A_1^3:A_3,A_1^2,A_1^2\right)$	0	between $2$ smooth points
			and an edge point
xiii)	$\left(A_1^4:A_1^2,A_1^2,A_1^2,A_1^2\right)$	0	between $4$ smooth points
xiv)	$\left(A_1{A}_3:A_1^2,A_1^2\right)$	0	$A_1{A}_3$ linking between $2$
			smooth points
	Pure Self-Linking
xv)	$\left(A_3:A_1^2\right)$	1	edge-type self-linking with
			a smooth point
xvi)	$\left(A_3:A_1^3\right)$	0	edge-type self-linking with
			a Y-junction point
xvii)	$\left(A_3:A_3\right)$	0	edge-type self-linking with
			an edge point

Equations

$${\lambda }_i(x,t)$$

$${\chi }_i(x,t)$$

Full text

The skeletal linking structures add increased flexibility for modeling configurations of objects by relaxing the Blum conditions and they extend in a minimal way the individual “skeletal structures”which have been previously used for modeling individual objects and capturing their geometric properties. This allows for the mathematical methods introduced for single objects to be significantly extended to the entire configuration of objects. These methods not only capture the internal shape structures of the individual objects but also the external structure of the neighboring regions of the objects.

In our subsequent second paper [DG2] we use these structures to identify specific external regions which capture positional information about neighboring objects, and we develop numerical measures for closeness of portions of objects and their significance for the configuration. This allows us to use the same mathematical structures to simultaneously analyze both the shape properties of the individual objects and positional properties of the configuration. This provides a framework for analyzing the statistical properties of collections of similar configurations such as for applications to medical imaging.

(1) Partially supported by the Simons Foundation grant 230298, the National Science Foundation grant DMS-1105470 and DARPA grant HR0011-09-1-0055. (2) This paper contains work from this author’s Ph. D. dissertation at Univ. of North Carolina

1. Introduction

Given a collection of objects, a basic question is how we may simultaneously model both the shapes of the individual objects and their positional information relative to one another. Motivating examples are provided by 2D and 3D medical images, in which we encounter collections of objects which might be organs, glands, arteries, bones, etc (see e.g. Figure 1). A number of researchers have already begun to recognize the importance of using the relative positions of objects in medical images to aid in analyzing physical features for diagnosis and treatment (see especially the work of Stephen Pizer and coworkers in MIDAG at UNC for both time series of a single patient and for populations of patients [LPJ], [GSJ], [JSM], [JPR], [Jg], and [CP]).

In these papers, significant use is made of the mathematical methods which are applied when modeling a single object using a “skeletal structure”[D1], [D3]. Such a structure generally consists of a “skeletal set”, which is a stratified set within the object, together with a multi-valued “radial vector field”  from points on the set to the boundary of the object. These form a flexible class of structures obtained by relaxing conditions on Blum medial axes, and these allow one to overcome various shortcomings of the Blum medial axis for modeling purposes. They have played an important role by providing mathematical tools and numerical criteria for fitting discrete deformable templates as models for the individual objects in 3D medial images. These discrete templates, called “M-reps” or in their more general forms “S-reps”, are discrete versions of skeletal structures where the skeletal set is a polyhedral surface with boundary edge, which is topologically a $2$-disk (versus the Blum medial axis which would have singularities), and with vectors defined at the vertices. In turn the fitting of the discrete deformable templates to individual objects in the medical images combine the numerical criteria with statistical priors formed from training sets of images. In carrying out the stages of fitting of the models to the objects in the images, the mathematics of skeletal structures is used to guarantee nonsingularity of the models, for the interpolation of the discrete models to yield smooth surfaces, and for ensuring the regularity of resulting object model boundaries.

This approach has been successfully used in numerous cases for modeling objects in medical images, including e.g. the hippocampus, brain ventricles, corpus callosum, kidney, liver, and pelvic area involving the prostate, bladder and rectum. For example, many references to papers that use M-reps and S-reps for the analysis of medical images of single objects are given in Chapters 1, 8 and 9 in the book edited and partially written by S. Pizer and K. Siddiqi, see reference [PS].

For certain collection of objects they have have been combined with user chosen, rather adhoc ways, to capture the relative positions between certain objects. While considerable work has been devoted to the application to individual objects, there has not been developed an effective mathematical approach for an entire collection. It is the goal of this paper to describe a mathematical framework for modeling an entire configuration of objects in 2D or 3D, which builds upon the success obtained for individual objects.

We will concentrate on configurations of objects in 2D and 3D images. These can be modeled by a collection of distinct compact regions $\{\Omega_{i}\}$ in ${\mathbb{R}}^{2}$ or ${\mathbb{R}}^{3}$ with piecewise smooth generic boundaries $\mathcal{B}_{i}$, where each pair of regions are either disjoint or only meet along portions of their boundaries in specified generic ways (see e.g. Figure 2). The mathematical results apply more generally to regions in ${\mathbb{R}}^{n}$ (see [DG]) even though we concentrate here on ${\mathbb{R}}^{2}$ and ${\mathbb{R}}^{3}$.

The geometric properties of the configuration are determined by both the shapes of the individual objects (regions) and their positions in the overall configuration. The “shapes”of the regions capture both the local and global geometry (and topology) of the regions. The overall “positional geometry”of the configuration involves such information as: the measure of relative closeness of portions of regions, characterization of “neighboring regions,”and the “relative significance”of an individual region within the configuration.

Single numerical measures such as the Gromov-Hausdorff distance between objects in such configurations measure overall differences without identifying specific feature variations responsible for the differences. Also invariants that would be appropriate for a collection of points fail to incorporate how shape differences between objects contribute to differences in positional structure of such configurations. We will introduce for such configurations “medial and skeletal linking structures,”which allow us to simultaneously capture local and global shape properties of the individual objects and their “positional geometry.”

In doing so we have three goals:

1. (1)

minimizing the redundancy of the geometric information provided by the structure while allowing existing mathematical methods for single regions to be extended to configurations; 2. (2)

robustness of the structure’s properties under small perturbations of the object’s shapes and positions in the configuration; and 3. (3)

obtaining quantitative measures from the structure for use in statistical analysis of populations of “configurations of the same type”.

To limit the amount of redundancy, the linking structures extend in a minimal way the notion of a “skeletal structure”for a single compact region $\Omega$ with smooth boundary $\mathcal{B}$, developed by the first author [D1] (or see [D3] for regions in ${\mathbb{R}}^{2}$ and ${\mathbb{R}}^{3}$). It consists of a pair $(M,U)$, where the “skeletal set”$M$ is a “Whitney stratified set”in the region and $U$ is a multivalued “radial vector field”defined on $M$. Skeletal structures are a more flexible class of structures which relax several of the conditions satisfied by the Blum medial axis of a region with smooth boundary [BN] The Blum medial axis is a special type of skeletal structure (with $U$ consisting of the vectors from points of $M$ to the points of tangency). The resulting skeletal structure still retains the mathematical tools for analyzing the shape and geometric properties of the associated region.

The added linking structure consists of a multi-valued “linking function”$\ell_{i}$ defined on the skeletal set $M_{i}$ for each region $\Omega_{i}$ and a refinement of the Whitney stratification of $M_{i}$ on which $\ell_{i}$ is stratawise smooth. The linking functions $\ell_{i}$ together with the radial vector fields $U_{i}$ yield multivalued “linking vector fields”$L_{i}$, which satisfy certain linking conditions. Even though the structures are defined on the skeletal sets within the regions, the linking vector fields allow us to capture geometric properties of the external region as well.

A special type of skeletal linking structure is a Blum medial linking structure, which builds upon the Blum medial axes of the individual regions. We classify the properties of the Blum linking structures for generic configurations. In obtaining this structure, we identify the regions which remain unlinked and classify their local generic structure of their singular boundaries (which is also the singular set of the convex hull of the configuration). We do this by introducing the spherical axis, which is the analog of the Blum medial axis but for directions in ${\mathbb{R}}^{n}$, instead of distances; and it is defined using the family of height functions on the region boundaries.

Just as for skeletal structures for single regions, the skeletal linking structure is a more flexible structure obtained by relaxing a number of the conditions for the Blum linking structure. This allows more flexibility in applying skeletal structures to model objects. Then, as we have already explained, this flexibility is what allows simplified skeletal structures to be used as discrete deformable templates for modeling objects [P], overcoming the lack of $C^{1}$-stability of the Blum medial axis, and providing discrete models to which statistical analysis can be applied [P2]. Likewise the skeletal linking structures have analogous properties, exhibiting stability under $C^{1}$ perturbations and allowing discrete models for applying statistical analysis. This motivates their potential usefulness for modeling and computer imaging questions for medicine and biology.

A skeletal structure enables both the local, relative, and global geometric properties of individual objects and their boundaries to be computed from the “medial geometry”of the radial vector field on the skeletal structure [D2] and [D4] (or see [D3] for regions in ${\mathbb{R}}^{2}$ and ${\mathbb{R}}^{3}$). This includes conditions ensuring the nonsingularity of a “radial flow”from the skeletal set to the boundary, providing a parametrization of the object interior by the level sets of the flow and implying the smoothness of the boundary [D1, Thm 2.5]. The medial/skeletal linking structure then extends these results to apply to the entire collection of objects, and includes the global geometry of the external region.

In Section 2, we explain what we mean by multi-object configurations and how we mathematically model them using multi-region configurations. Next, in Section 3 we introduce medial/skeletal linking structures for multi-region configurations and give their basic properties. In Section 4, we introduce two versions of a “Blum linking structure”for a general configuration. We first give the generic properties of the Blum medial axis for a region with singular boundary having corners and edge curves (in the case of ${\mathbb{R}}^{3}$). The Blum medial axis will now contain the singular points of the boundaries in its closure and we give an “edge-corner normal form”that the Blum medial axis will exhibit near such singular points. We further give the local form of the stratification of the boundary by points associated to the singular points of the medial axis. It is the generic interplay between these stratifications for adjoining regions and the linking medial axis that gives the generic properties of the linking structure.

For a generic multi-object configuration, which allows shared boundaries that have singular points, we establish the existence of a generic “full Blum linking structure”for the configuration (Theorem 4.13); and later list the generic linking types for the 3D case in Section 6. In the special case where all of the regions are disjoint with smooth generic boundaries, this directly yields a “Blum medial linking structure”(Theorem 4.12). This special case for disjoint regions was obtained in the thesis of the second author [Ga]. In Section 5 we explain how to modify for general configurations the resulting full Blum linking structure near the singular points of the boundaries to obtain a skeletal linking structure. Lastly in Section 7 we explain how the method involving M-reps and S-reps used as deformable templates for single regions can be expanded to obtain deformable templates for an entire configuration based on skeletal linking structures.

In a second paper [DG2], we will use the linking structure we have introduced to determine properties of the “positional geometry”of the configuration. This will include: identifying distinguishing external regions capturing positional geometry; identifying neighboring regions via linking between these regions; introducing and computing numerical invariants of the positional geometry for measures of closeness and significance of regions in terms of volumetric measurements; computing volumetric invariants (which involve regions outside the configuration) as “skeletal linking integrals”; and combining these results to obtain a tiered graph structure which provides a hierarchical ordering of the regions. This will provide for the comparison and statistical analysis of collections of objects in ${\mathbb{R}}^{2}$ and ${\mathbb{R}}^{3}$.

The authors would like to thank Stephen Pizer for sharing with us his early work with his coworkers involving multiple objects in medical images. This led us to seek a completely mathematical approach to these problems.

2. Modeling Multi-Object Configurations in ${\mathbb{R}}^{2}$ and ${\mathbb{R}}^{3}$

Local Models for Objects at Singular Points on Boundary

We begin by defining what exactly we mean by a “multi-object configuration.”We consider objects which have smooth boundaries if they are separated from the other objects. However, if two objects meet along their boundaries then we allow two different situations, where objects either have flexible boundaries or are rigid. The resulting possible configurations depend upon which combinations occur. To model such meetings along boundaries, we first describe the form that the object’s boundaries take at the edges of such shared boundary regions.

First, we model the individual objects by compact connected 2D or 3D regions $\Omega$ with boundaries $\mathcal{B}$ which are smooth manifolds with boundaries and corners. We say that $\Omega$ is a manifold with boundaries and corners if each point $x\in\mathcal{B}$ is either a smooth point of the boundary or is modeled, via a diffeomorphic mapping of a neighborhood of the origin of ${\mathbb{R}}^{2}$, resp. ${\mathbb{R}}^{3}$, sending either a closed quadrant or half space of ${\mathbb{R}}^{2}$, resp. a closed octant, quarter-space, or half-space of ${\mathbb{R}}^{3}$ to a neighborhood of $x$ in $\Omega$, with the boundaries of these regions in ${\mathbb{R}}^{2}$, resp. ${\mathbb{R}}^{3}$, mapping to the boundary $\mathcal{B}$, see e.g. Figures 3 and 4.

Then, the boundary $\mathcal{B}$ is stratified by the corner points, the edge curves (for ${\mathbb{R}}^{3}$), and the smooth boundary regions (referred to as regular points of the boundary). The union of the corner points and edge curves (in ${\mathbb{R}}^{3}$) form the singular points of the boundary.

Second, we describe how we model objects sharing common boundary regions. If two such regions intersect it will only be on their boundaries, and to describe the common boundary regions, we consider two cases depending on whether the regions are all flexible, or one of them is rigid. In the first case, for 3D objects, we model the edges of regions of common boundary regions by either a) $P_{2}$; b) $P_{3}$ in Figure 6. In these figures, one of the regions in the figures may denote the external region complementary to the objects. If one of the regions is rigid then instead the regions meet as in c) $Q_{1}$ or d) $Q_{2}$ in Figure 6, where the region with smooth boundary is the rigid one, and again one of the other regions may denote the external region complementary to the objects. For 2D objects, the local models for meeting are more simply given by either points of type $P_{2}$ in a) and b) in Figure 5 or type $Q_{2}$ in c) and d) in the same figure. We refer to a configuration of regions satisfying the above two conditions regarding their boundary structure and their common boundary regions as satisfying the combined boundary intersection condition.

Remark 2.1**.**

Physically such local configurations of type a) or b) in Figures 5 or 6 are generic when objects with flexible boundaries have physical contact. For example, these are typical of the singularities occurring where multiple soap bubbles are joined. In general Maxwell’s Principle from physics implies that for boundaries of “conflicting regions”  the singular boundaries have such specific forms (determined by a singularity theoretic analysis of the energy potential functions). For c) and d), the region below the curve in Figure 5 or to the left of the plane in Figure 6 would represent a rigid region with which one or more flexible regions have contact.

We are going to model a multi-object configuration mathematically with a multi-region configuration, which we define next.

Definition 2.2**.**

A multi-region configuration consists of a collection of compact regions $\Omega_{i}\subset{\mathbb{R}}^{2}$ or ${\mathbb{R}}^{3}$, $i=1,...,m$, with smooth boundaries and corners $\mathcal{B}_{i}$ and satisfying the combined boundary intersection condition.

Then, a multi-region configuration will encode as data : i) the topological structure of each region; ii) the regions which share common boundaries and their individual types (flexible or rigid); and iii) the topological structure of the shared boundary regions between two regions. Different regions having the same data will be viewed as having the same configuration type, and changes in the type will be viewed as resulting from a transition of types, see e.g. §5.

The configuration type can change due to: regions that were disjoint coming in contact, or the reverse; boundary regions of contact undergoing topological changes; regions undergoing changes in their topological type (e.g. a hole being created or destroyed); or more subtle positional changes. These are modeled by transitions of configuration type. The generic transitions are the most common ones, but we will not attempt to classify or analyze them here. If we wish to study variations and modifications within a single configuration type, then these can be most effectively described via an embedding mapping of a configuration which we consider next.

The Space of Equivalent Configurations via Mappings of a

Model

We introduce a model configuration $\boldsymbol{\Delta}$ for $\boldsymbol{\Omega}$ in ${\mathbb{R}}^{n+1}$, $n=1,2$ which is a configuration of multi-regions $\{\Delta_{i}\}$ satisfying Definition 2.2. A configuration of the same type as $\boldsymbol{\Omega}$ will be obtained by an embedding $\Phi:\boldsymbol{\Delta}\to{\mathbb{R}}^{n+1}$. This means that $\Phi:\cup_{i}\Delta_{i}\to{\mathbb{R}}^{n+1}$ extends to a smooth embedding in some neighborhood, and it restricts to diffeomorphisms $\Delta_{i}\simeq\Omega_{i}$ for each $i$, with $\mathcal{B}_{i}$ denoting the boundary of $\Omega_{i}$. Even though the configuration varies with $\Phi$ we still use the notation $\boldsymbol{\Omega}$ for the resulting (varying) image configuration (with a specific $\Phi$ understood). In particular, a configuration of model type $\boldsymbol{\Delta}$ has all of the data of $\boldsymbol{\Delta}$ built into it, so different regions $\Omega_{i}$ and $\Omega_{j}$ will meet along a shared boundary only when the corresponding $\Delta_{i}$ and $\Delta_{j}$ already do so and otherwise will not meet.

The space of configurations of type $\boldsymbol{\Delta}$ is given by the infinite dimensional space of embeddings ${\rm Emb}\,(\boldsymbol{\Delta},{\mathbb{R}}^{n+1})$. Then, by generic properties of a configuration we will mean properties satisfied for an open dense set of embeddings in ${\rm Emb}\,(\boldsymbol{\Delta},{\mathbb{R}}^{n+1})$. This means that for a configuration that exhibits generic properties, a sufficiently small perturbation will not destroy these properties, and if a configuration does not exhibit these properties, “almost any arbitrarily small perturbation”which is performed will ensure that it does. The generic properties which we give are shown in [DG] to satisfy this property.

We say the regions $\Omega_{i}$ and $\Omega_{j}$ are adjoining regions if $\mathcal{B}_{i}\cap\mathcal{B}_{j}\neq\emptyset$. Also, the external boundary region of $\mathcal{B}_{i}$ will denote that portion of the boundary which is shared by the external region. For example, in the multi-region configuration $\boldsymbol{\Omega}$ in Figure 2 the regions $\Omega_{1}$, $\Omega_{5}$, $\Omega_{6}$, are adjoined to $\Omega_{2}$ (and $\Omega_{1}$ and $\Omega_{5}$ with each other), and all of the regions have external boundary regions. These features will persist for any configuration of the same type as this particular $\boldsymbol{\Omega}$.

Remark 2.3**.**

In the special case of a configuration consisting of disjoint regions with smooth boundaries, each boundary is entirely external. In such a case, we shall see that the geometric relations between the regions are entirely captured via “linking behavior”  in the external region.

Because the external region extends indefinitely, we will frequently find for computational reasons that it is preferable to have the configuration contained in a “bounded region”. By this we mean an ambient region $\tilde{\Omega}$ so that $\Omega_{i}\subset\tilde{\Omega}$ for each $i$, then we say that $\boldsymbol{\Omega}$ is a configuration bounded by $\tilde{\Omega}$. Such a $\tilde{\Omega}$ might be a bounding box or disk or an intrinsic region containing the configuration. Then we will also consider bounded configurations either given by an embedding $\tilde{\Phi}:\tilde{\Delta}\to{\mathbb{R}}^{2}$ or ${\mathbb{R}}^{3}$, with $\tilde{\Phi}(\tilde{\Delta})$ denoting $\tilde{\Omega}$, and $\Phi=\tilde{\Phi}|\boldsymbol{\Delta}$; or we fix $\tilde{\Omega}$ and consider embeddings $\Phi:\boldsymbol{\Delta}\to{\rm int}\,(\tilde{\Omega})$.

Configurations Allowing Containment of Regions

In our definition of a multi-region configuration, we have explicitly excluded one region being contained in another. However, given a configuration which allows this, we can easily identify such a configuration with the type we have already given. To do so, if one region is contained in another $\Omega_{i}\subset\Omega_{j}$, then we may represent $\Omega_{j}$ as a union of two regions $\Omega_{i}$ and the closure of $\Omega_{j}\backslash({\rm int}\,(\Omega_{i})\cup(\mathcal{B}_{i}\cap\mathcal{B}_{j}))$, which we refer to as the region complement to $\Omega_{i}$ in $\Omega_{j}$. By repeating this process a number of times we arrive at a representation of the configuration as a multi-region configuration in the sense of Definition 2.2. See Figure 7.

3. Skeletal Linking Structures for Multi-Region Configurations in

${\mathbb{R}}^{2}$ and ${\mathbb{R}}^{3}$

The skeletal linking structures for multi-object configurations will be constructed for the multi-region configurations modeling them. They will build upon the skeletal structures for individual regions. We begin by recalling their basic definitions and simplest properties.

Skeletal Structures for Single Regions

We begin by recalling [D1] (or also see [D3]) that a “skeletal structure”$(M,U)$ in ${\mathbb{R}}^{n+1}$ [D1, Def. 1.13], consists of: a Whitney stratified set $M$ which satisfies the conditions for being a “skeletal set”  [D1, Def. 1.2] and a multivalued “radial vector field”  $U$ on $M$ which satisfies the conditions of [D1, Def. 1.5] and the “local initial conditions”[D1, Def. 1.7]. $M$ consists of smooth strata of dimension $n$, and the set of singular strata $M_{sing}$, with $\partial M$ denoting the singular strata where $M$ is locally a manifold with boundary (for which we use special “boundary coordinates”). For images in ${\mathbb{R}}^{2}$ and ${\mathbb{R}}^{3}$, we are interested in the cases for $n=1,2$; however, the general description is independent of dimension.

Without restating the conditions, we remark that these conditions allow each value of $U$ on a smooth point $x$ of $M$ to extend to a smooth vector field $U$ on a neighborhood of $x$, and various mathematical constructions on the smooth strata to be extended to the singular strata $M_{sing}$, see [D1, §2]. Using the multivalued radial vector field $U$, we define a stratified set $\tilde{M}$, called the “double of $M$”  and a finite-to-one stratified mapping $\pi:\tilde{M}\to M$, see [D1, §3]. Points of $\tilde{M}$ consist of all pairs $\tilde{x}=(x,U)$ where $U$ is a value of the radial vector field at $x$, and $\pi(x,U)=x$ (see Figure 8). $\tilde{M}$ provides a mathematical method to keep track of both sides of each stratum of $M$. It also allows multivalued objects on $M$ to be viewed as single-valued objects on $\tilde{M}$.

Likewise we can keep track of the radial lines from each side of $M$ with the space $N_{+}=\tilde{M}\times[0,\infty)$, where for $\tilde{x}=(x,U)\in\tilde{M}$, the radial line (actually a half-line) is parametrized by $[0,\infty)$ by $c\mapsto c\cdot U$. Then, using $N_{+}$, we can define the “radial flow”. In a neighborhood $W$ of a point $x_{0}\in M$ with a smooth single-valued choice for $U$, we define a local representation of the radial flow by $\psi_{t}(x)=x+t\cdot U(x)$. Together the local definitions yield a global radial flow as a mapping $\Psi:N_{+}\to{\mathbb{R}}^{n+1}$ defined for $\tilde{x}=(x,U)$ by $\Psi(\tilde{x},t)=\psi_{t}(x)$. Beginning with a skeletal structure $(M,U)$, we can associate a “region”  $\Omega$ which is the image $\Psi(N_{1})$, where $N_{1}=\tilde{M}\times[0,1]$, and its “boundary”  $\mathcal{B}=\{x+U(x):x\in M\,\mbox{ all values of }U\}$.

Provided certain curvature and compatibility conditions are satisfied, then by [D1, Thm 2.5] the radial flow defines a stratawise smooth diffeomorphism between $N_{1}\backslash\tilde{M}$ and $\Omega\backslash M$; however the level sets of the flow $\mathcal{B}_{t}$ are stratified rather than being smooth. From the radial flow we define the radial map $\psi_{1}(x)=x+U(x)$ from $\tilde{M}$ to $\mathcal{B}$. We may then relate the boundary and skeletal set via the radial flow and the radial map.

A standard example we consider will be the Blum medial axis $M$ of a region $\Omega$ with generic smooth boundary $\mathcal{B}$ and its associated (multivalued) radial vector field $U$. Then, the associated boundary $\mathcal{B}$ will be the original boundary of the object.

The compatibility condition involves the compatibility $1$-form $\eta_{U}=dr+\omega_{U}$, where $\|U\|=r$ is the radial function and $\omega_{U}(v)=v\cdot\mathbf{u}$ where $U=r\cdot\mathbf{u}$, with $\mathbf{u}$ the multivalued unit vector field along $U$. Then, $\eta_{U}$ is also multivalued with one value for each value of $U$. If $\eta_{U}$ vanishes on a neighborhood of $x$ in a smooth region of $M$, then $U(x)$ is orthogonal to the boundary $\mathcal{B}$ at the point $\psi_{1}(x)$, and we refer to $U$ as being “partially Blum”at $x$. Also, if $\eta_{U}$ vanishes on $M_{sing}$, then the boundary $\mathcal{B}$ will be weakly $C^{1}$ on the image of the singular set. However, if it does not, then the boundary will have corners and edges. Hence, skeletal structures can also be used for regions with corners and edges.

Skeletal Linking Structures for Multi-Region Configurations

We are ready to introduce skeletal linking structures for multi-region configurations. These structures will accomplish multiple goals. The most significant are the following.

• i)

Extend the skeletal structures for the individual regions in a minimal way to obtain a unified structure which also incorporates the positional information of the objects.

• ii)

For generic configurations of disjoint regions with smooth boundaries, provide a Blum medial linking structure which incorporates the Blum medial axes of the individual objects.

• iii)

For general multi-region configurations with common boundaries, provide for a modification of the resulting Blum structure to give a skeletal linking structure.

• iv)

Handle both unbounded and bounded multi-region configurations.

In Part II [DG2] we shall also see that the skeletal linking structure has a second important function allowing us to answer various questions involving the “positional geometry”of the regions in the configuration.

Remark 3.1**.**

Even though a configuration itself is bounded, the complement is unbounded. The properties we shall give will be valid for the unbounded external region; however, for both computational and measurement reasons, it is desirable to consider the configuration as lying in a bounding region. This is the “bounded case”. There are a number of possibilities for such bounding regions: intrinsic bounding region, a bounding box, the convex hull, or a bounded region resulting from imposing a threshold, etc. We shall consider further the relation of these with a skeletal linking structure modeling a given configuration at the end of this section.

We begin by giving versions of the definition for both the bounded and unbounded cases.

Definition 3.2**.**

A skeletal linking structure for a multi-region configuration $\{\Omega_{i}\}$ in ${\mathbb{R}}^{2}$ or ${\mathbb{R}}^{3}$ consists of a triple $(M_{i},U_{i},\ell_{i})$ for each region $\Omega_{i}$ with the following properties.

• S1)

$(M_{i},U_{i})$ is a skeletal structure for $\Omega_{i}$ for each $i$ with $U_{i}=r_{i}\cdot\mathbf{u}_{i}$ a scalar multiple of $\mathbf{u}_{i}$ a (multivalued) unit vector field and $r_{i}$ the multivalued radial function on $M_{i}$.

• S2)

$\ell_{i}$ is a (multivalued) linking function defined on $M_{i}$ (excluding the strata $M_{i\,\infty}$, see L4 below), with one value for each value of $U_{i}$, for which the corresponding values satisfy $\ell_{i}\geq r_{i}$, and it yields a (multivalued) linking vector field $L_{i}=\ell_{i}\cdot\mathbf{u}_{i}$ .

• S3)

The canonical stratification of $\tilde{M}_{i}$ has a stratified refinement $\mathcal{S}_{i}$, which we refer to as the labeled refinement.

By $\mathcal{S}_{i}$ being a “labeled refinement”of the stratification $\tilde{M}_{i}$ we mean it is a refinement in the usual sense of stratifications in that each stratum of $\tilde{M}_{i}$ is a union of strata of $\mathcal{S}_{i}$; and they are labeled by the linking types which occur on the strata.

In addition, they satisfy the following four linking conditions.

Conditions for a skeletal linking structure

• L1)

$\ell_{i}$ and $L_{i}$ are continuous where defined on $M_{i}$ and are smooth on strata of $\mathcal{S}_{i}$.

• L2)

The “linking flow”(see (3.1) below) obtained by extending the radial flow is nonsingular and for the strata $S_{i\,j}$ of $\mathcal{S}_{i}$, the images of the linking flow are disjoint and each $W_{i\,j}=\{x+L_{i}(x):x\in S_{i\,j}\}$ is smooth.

• L3)

The strata $\{W_{i\,j}\}$ from the distinct regions either agree or are disjoint and together they form a stratified set $M_{0}$, which we shall refer to as the (external) linking axis.

• L4)

There are strata $M_{i\,\infty}\subset\tilde{M}_{i}$ on which there is no linking so the linking function $\ell_{i}$ is undefined. On the union of these strata $M_{\infty}=\cup_{i}M_{i\,\infty}$, the global radial flow restricted to $N_{+}|M_{\infty}$ is a diffeomorphism with image the complement of the image of the linking flow. In the bounded case, with $\tilde{\Omega}$ the enclosing region of the configuration, it is required that the boundary of $\tilde{\Omega}$ is transverse to the stratification of $M_{0}$ and where the linking vector field extends beyond $\tilde{\Omega}$, it is truncated at the boundary of $\tilde{\Omega}$ (this includes $M_{\infty}$).

We denote the region on the boundary corresponding to $M_{\infty}$ by $\mathcal{B}_{\infty}$ and that corresponding to $M_{i\,\infty}$ by $\mathcal{B}_{i\,\infty}$.

Remark 3.3**.**

By property L4), $M_{i\,\infty}$ does not exhibit any linking with any other region. We will view it as either the unlinked region or alternately as being linked to $\infty$, where we may view the linking function as being $\infty$ on $M_{i\,\infty}$. In the bounded case we modify the linking vector field so it is truncated at the boundary of $\tilde{\Omega}$. We can also introduce a “linking vector field on $M_{\infty}$”  by extending the radial vector field until it meets the boundary of $\tilde{\Omega}$ (see e.g. Figure 12).

For this definition, we must define the “linking flow”  which is an extension of the radial flow. We define the linking flow from $M_{i}$ by

[TABLE]

As with the radial flow it is actually defined from $\tilde{M}_{i}$ (or $\tilde{M}_{i}\backslash M_{\infty}$). The combined linking flows $\lambda_{i}$ from all of the $M_{i}$ will be jointly referred to as the linking flow $\lambda$. For fixed $t$ we will frequently denote $\lambda(\cdot,t)$ by $\lambda_{t}$.

**Convention: ** It is convenient to view the collection of objects for the linking structure as together forming a single object, so we will adopt notation for the entire collection. This includes $M$ for the union of the $M_{i}$ for $i>0$, and similarly for $\tilde{M}$. On $M$ (or $\tilde{M}$) we have the radial vector field $U$ and radial function $r$ formed from the individual $U_{i}$ and $r_{i}$, the linking function $\ell$ and linking vector field $L$ formed from the individual $\ell_{i}$ and $L_{i}$; as well as the linking flow $\lambda$ and $M_{\infty}$ already defined.

We see that for $0\leq t\leq\frac{1}{2}$ the flow is the radial flow at twice its speed; hence, the level sets of the linking flow, $\mathcal{B}_{i\,t}$, for time $0\leq t\leq\frac{1}{2}$ will be those of the radial flow. For $\frac{1}{2}\leq t\leq 1$ the linking flow is from the boundary $\mathcal{B}_{i}$ to the linking strata of the external medial linking axis.

By the linking flow being nonsingular we mean it is a piecewise smooth homeomorphism, which for each stratum $S_{i\,j}\subset\tilde{M}_{i}$, is smooth and nonsingular on $S_{i\,j}\times\left[0,\frac{1}{2}\right]$ and either: $S_{i\,j}\times\left[\frac{1}{2},1\right]$ is smooth and nonsingular if $S_{i\,j}$ is a stratum associated to strata in $\mathcal{B}_{i\,0}$. For $S_{i\,j}$ not associated to strata in $\mathcal{B}_{i\,0}$, $\ell_{i}=r_{i}$ on $S_{i\,j}$, so the linking flow on $S_{i\,j}\times\left[\frac{1}{2},1\right]$ is constant as a function of $t$. That the linking flow is nonsingular will follow from the analogue of the conditions given in [D1, §3] for the nonsingularity of the radial flow. These will be given in Part II [DG2], when we use the linking flow to establish geometric properties of the configuration.

Linking between Regions and between Skeletal Sets

A skeletal linking structure allows us to introduce the notion of linking of points in different (or the same) regions and of regions themselves being linked. We say that two points $x\in\tilde{M}_{i}$ and $x^{\prime}\in\tilde{M}_{j}$ are linked if the linking flows satisfy $\lambda_{i}(x,1)=\lambda_{j}(x^{\prime},1)$. This is equivalent to saying that for the values of the linking vector fields $L_{i}(x)$ and $L_{j}(x^{\prime})$, $x+L_{i}(x)=x^{\prime}+L_{j}(x^{\prime})$. Then, by linking property $L3)$, the set of points in $\tilde{M}_{i}$ and $\tilde{M}_{j}$ which are linked consist of a union of strata of the stratifications $\mathcal{S}_{i}$ and $\mathcal{S}_{j}$. Furthermore, if the linking flows on strata from $S_{i\,k}\subset\tilde{M}_{i}$ and $S_{j\,k^{\prime}}\subset\tilde{M}_{j}$ yield the same stratum $W\subset M_{0}$, then we refer to the strata as being linked via the linking stratum $W$. Then, $\mu_{i\,j}=\lambda_{j}(\cdot,1)^{-1}\circ\lambda_{i}(\cdot,1)$ defines a diffeomorphic linking correspondence between $S_{i\,k}$ and $S_{j\,k^{\prime}}$.

In Part II [DG2] we will introduce a collection of regions which capture geometrically the linking relations between the different regions. For now we concentrate on understanding the types of linking that can occur. There are several possible different kinds of linking. More than two points may be linked at a given point in $M_{0}$. Of these more than one may be from the same region. If all of the points are from a single region, then we call the linking self-linking, which occurs at indentations of regions. If there is a mixture of self-linking and linking involving other regions then we refer to the linking as partial linking, see Figure 11.

Remark 3.4**.**

If $\Omega_{i}$ and $\Omega_{j}$ share a common boundary region, then certain strata in $M_{i}$ and $M_{j}$ are linked via points in this boundary region, and for those $x\in M_{i}$, $\ell_{i}(x)=r_{i}(x)$; see Figure 10.

We next consider in more detail the bounded case.

Skeletal Linking Structures in the Bounded Case

For the bounded case we suppose that the configuration lies in the interior of a “bounding region”$\tilde{\Omega}$ so that the boundary $\partial\tilde{\Omega}$ is transverse to: i) the strata of the external medial axis $M_{0}$, ii) the extension of the radial lines from $M_{\infty}$, and iii) the linking line segments. If the boundary $\partial\tilde{\Omega}$ has singular points, then for ii) and iii) by the lines being transverse to $\partial\tilde{\Omega}$, we mean that at a singular point, the line is transverse to the limits of tangent spaces from the smooth points of $\partial\tilde{\Omega}$. For any convex region $\tilde{\Omega}$ with piecewise smooth boundary, the limiting tangent planes of $\partial\tilde{\Omega}$ are supporting hyperplanes for $\tilde{\Omega}$. Hence any line in the tangent plane lies outside ${\rm int}\,(\tilde{\Omega})$. Thus, the radial lines from regions in the configuration will always meet the boundary transversely.

We can alter the linking vector field either by truncating it where it meets the boundary $\partial\tilde{\Omega}$ or defining it on $M_{\infty}$ and then refining the stratification so that on appropriate strata the linking vector field ends at $\partial\tilde{\Omega}$. It is now defined on all of $M_{i}$ for the individual regions $\Omega_{i}$. Because we are either reducing $\ell_{i}$, or defining $L_{i}$ on $M_{i\,\infty}$, the linking flow is still nonsingular.

Thus, we have compact versions of the regions defined for the unbounded case. There are a number of different possibilities for such bounding regions.

Possibilities for a Bounding Region $\tilde{\Omega}$:

*Bounding Box or Bounding Convex Region: * The box requires a center and directions and sizes for the edges of the box. For this, we would need to first normalize the center and directions for the sides of the box and then normalize the sizes of the edges either using a fixed size or one based on feature sizes of the configuration (see e.g. Figure 12 (a)).

• 2)

*Convex Hull: * The smallest convex region which contains a configuration is the convex hull of the configuration. For a generic configuration, the convex hull consists of the regions $\mathcal{B}_{i\,\infty}$ together with the truncated envelope of the family of degenerate supporting hyperplanes which meet the configuration with a degenerate tangency or at multiple points. a) of Figure 17 in § 4, which is a configuration in ${\mathbb{R}}^{2}$, the envelope consists of line segments joining the doubly tangent points corresponding to the points in the spherical axis (b) of Figure 17. In ${\mathbb{R}}^{3}$, the envelope consists of a triangular portion of a triply tangent plane, with line segments joining pairs of points with a bitangent supporting plane, and a decreasing family of segments ending at a degenerate point.

• 3)

*Intrinsic Bounding Region: * If the configuration is naturally contained in an (a priori given) intrinsic region, which is modeled by $\tilde{\Omega}$, then provided that the extensions of the radial lines intersect $\partial\tilde{\Omega}$ transversely, then we can modify the linking structure as in the convex case to have it defined on all of the $M_{i}$, and terminating at $\partial\tilde{\Omega}$, if linking has not already occurred (see e.g. Figure 12 (b)).

• 4)

*Threshold for Linking: * The external region can be bounded by placing a threshold $\tau$ on the $\ell_{i}$ so that the $L_{i}$ remain in a bounded region. This can be done in two different ways. An absolute threshold restricts the external regions to those arising from linking vector fields with $\ell_{i}\leq\tau$; and a truncated threshold restricts to a bounded region formed by replacing $\ell_{i}$ by $\ell_{i}^{\prime}=\min\{\ell_{i},\tau\}$. For the first type, only part of the region would have an external linking neighborhood; while in the second, the entire region would. As for the convex case, since we are replacing $\ell_{i}$ by smaller values, the linking flow will remain nonsingular. We obtain modified versions of the regions lying in a bounded region.

4. Blum Medial Linking Structure for a Generic Multi-Region

Configuration

In this section we consider two types of linking structures extending the Blum medial axes of the individual objects (i.e. regions). To do so we give a number of results. These include: giving the generic structure of the Blum medial axis for a region with singular boundary allowing corners and edges (in ${\mathbb{R}}^{3}$); introducing the “spherical axis”for the configuration, which plays the analogous role for directions as the medial axis plays for distances; and giving the generic properties for the “full Blum medial linking structure”for a general configuration, and its special form in the case of a configuration of disjoint regions with smooth boundaries.

Blum Medial Axis for a Single Region with Generic Singular

Boundary

For a single region we will extend the classification of the generic local structure of the Blum medial axis for regions with smooth boundaries in ${\mathbb{R}}^{2}$ and ${\mathbb{R}}^{3}$ due to Yomdin [Y], Mather [M], and Giblin-Kimia [GK] (see also [P] and [PS]). For ${\mathbb{R}}^{2}$ the singular points are either $Y$-branch points where $3$ smooth curves meet non-tangentially at a point, or end points of a smooth curve. For ${\mathbb{R}}^{3}$ the singular points are either edge points ($A_{3}$), fin points ($A_{1}A_{3}$), $Y$-branch curve points ($A_{1}^{3}$), or $6$-junction points ($A_{1}^{4}$), where six surface segments meet at a point along four $Y$-branch curves. These are illustrated in Figure 13. The notation $A_{1}A_{3}$ or $A_{1}^{k}$ denotes the singular behavior of the “distance to the boundary”function at the points on the boundary associated to the indicated points on the medial axis.

Stratification of the Boundary

The associated points on the boundary form a stratification of the boundary classified by the types of singular points. The boundary classification is shown in [DG, Thm 4.4] to have the generic forms given in Figure 14.

Stratification of the Blum Medial Axis at Singular

Boundary Points

Unlike the case of a region with smooth boundary, the Blum medial axis of a region with singular boundary may extend up to the singular portion. However, in the case of singular points which are either corners or edges (in ${\mathbb{R}}^{3}$) we can give a precise characterization of the Blum medial axis near these points.

Definition 4.1**.**

The edge-corner normal form for the Blum medial axis of an edge or corner point $x$ in ${\mathbb{R}}^{2}$ or ${\mathbb{R}}^{3}$ is the image via a smooth diffeomorphism, sending the origin to $x$, of one of the models given in Figure 15. In ${\mathbb{R}}^{2}$, the model for the medial axis is given for the corner model by the set of points $(x,y)$ with $x=y\geq 0$. In ${\mathbb{R}}^{3}$, the models are given by either: the set of points $(x,y,z)$ with $x=y\geq 0$ for the edge model, or by the set of points $(x,y,z)$ with $x,y,z\geq 0$ and one of $x=y\geq 0$, $x=z\geq 0$, or $y=z\geq 0$.

We give the following strengthened form for the structure of the Blum medial axis allowing singular points, see [DG, Thm 4.4]

Theorem 4.2** (Generic Properties of the Blum Medial Axis).**

For a generic region $\Omega$ in ${\mathbb{R}}^{2}$ or ${\mathbb{R}}^{3}$ with singular boundary allowing corners and edges (for ${\mathbb{R}}^{3}$), the Blum medial axis $M$ has the following properties:

• i)

it has the same local properties in the interior of the region ${\rm int}\,(\Omega)$ as given above for smooth regions (including the local singular structure and properties of the radial vector fields as illustrated for 3D in Figure 13);

• ii)

the corresponding stratifications of the boundary lie in the smooth strata (i.e. they miss the edges and corners) and are the same as for the smooth case; and

• iii)

the closure of the Blum medial axis in the boundary consists of the singular points of the boundary; at these points, the Blum medial axis has the corresponding edge-corner normal form in Definition 4.1.

Remark 4.3**.**

The Blum medial structure $M$ for a region with boundary and corners contains the singular points of the boundary in its closure, and at such points the radial vector field $U=0$. Hence, $(M,U)$ does not define a skeletal structure in the strict sense. However, it can still be used, just as with skeletal structures, to compute the local, relative, and global geometry and topology of both the region and its boundary. Hence, we can view it as a “relaxed skeletal structure,”where “relaxed”means that $M$ includes the singular boundary points and $U=0$ on these points. Alternatively, we can modify the Blum medial axis near the singular points of the boundary so it becomes a skeletal structure, see §5.

Spherical Axis of a Configuration

The generic local properties of the Blum medial axis for individual regions described in the last section will also apply to the external linking axis for the Blum linking structure which will be intrinsic for a configuration. There are still two remaining contributions to the linking structure. The one involves identifying the types of linking which may occur generically between different regions. This leads to part of the refinement of the stratification of the Blum medial axes of the individual regions and the external linking axis. The remaining contribution concerns the portions of the Blum medial axes, denote $M_{\infty}$, and their corresponding boundary points $\mathcal{B}_{\infty}$ where no linking occurs. These are points on the boundaries where no other regions (or different portions of the same region) would be visible. We characterize the boundary of $\mathcal{B}_{\infty}$ using the “spherical axis”and the associated “spherical structure”, which we now proceed to define.

Along with the Blum medial axis, we introduce its analog for directions in place of distances. Directions are given by vectors $\mathbf{v}\in S^{n}$, where $S^{n}$ is the unit sphere in ${\mathbb{R}}^{n+1}$. Given such a direction and either a region $\Omega$ or a configuration $\boldsymbol{\Omega}$, the supporting hyperplane in the direction $\mathbf{v}$ is given by an equation $x\cdot\mathbf{v}=c$, where all points $x$ of $\Omega$ or $\boldsymbol{\Omega}$ satisfy $x\cdot v\leq c$ (i.e. it is contained in the half-space defined by $x\cdot\mathbf{v}\leq c$), with equality at one or more points. The intersection of the supporting hyperplanes defines the convex hull of the region or configuration.

We define the spherical axis $\mathcal{Z}\subset S^{n}$ of $\Omega$, or the configuration $\boldsymbol{\Omega}$, to consist of directions $\mathbf{v}\in S^{n}$ for which the supporting hyperplanes $x\cdot\mathbf{v}=c$ for the convex hull of $\Omega$ or $\boldsymbol{\Omega}$ have two or more tangencies with $\mathcal{B}$ or a degenerate tangency, see e.g. Figure 16. If we let $\tau:\mathcal{B}\times S^{n}\to{\mathbb{R}}$ be defined by $\tau(x,\mathbf{v})=x\cdot\mathbf{v}$, then $\tau$ is the family of “height functions” on $\mathcal{B}$. The spherical axis is the set of $v\in S^{n}$ at which the absolute maximum of $\tau(\cdot,v)$ occurs at multiple points or is a degenerate maximum (this is the “Maxwell set”of $-\tau$ ). Then, similar methods to those used by Mather can be applied to determine the generic structure of $\mathcal{Z}$ (see [DG, Thm 4.6]). In the special cases of ${\mathbb{R}}^{2}$ and ${\mathbb{R}}^{3}$ they give the following.

Theorem 4.4** (Generic Properties of the Spherical Axis).**

For generic regions (with smooth boundaries) or generic configurations, the spherical axis has the following local structure. In ${\mathbb{R}}^{2}$, $\mathcal{Z}$ consists of a discrete set of points where supporting lines have two tangent points (see e.g. Figure 17). In ${\mathbb{R}}^{3}$, it consists of a collection of branched curves with endpoints. At smooth points of the curves there is a double tangency of the supporting plane. The branching occurs at $Y$-branch points which correspond to triple tangencies. The end points correspond to points where the height function $\tau$ has an $A_{3}$ singularity. These are illustrated in Figure 16.

Spherical Structure for a Configuration

Just as for the Blum medial axis, we may associate to the spherical axis $\mathcal{Z}$, both a height function $h$ and a multivalued vector field $V$. To do this we need to initialize the origin. Then, the height for a point $\mathbf{u}$ on the spherical axis, which is a unit vector, defined for a point $x\in\mathcal{B}$ is just the dot product $x\cdot\mathbf{u}$. Furthermore, there may be multiple points associated to $\mathbf{u}$, all of which lie in the supporting hyperplane $H$ defined by $x\cdot\mathbf{u}=h(\mathbf{u})$, for the maximum value $h=h(\mathbf{u})$ for $\mathcal{B}$. For each point $x_{i}\in H\cap\mathcal{B}$, there is a vector $V=x_{i}-h\mathbf{u}$ orthogonal to the line spanned by $\mathbf{u}$. This defines a multivalued vector field $V$.

Definition 4.5**.**

The full spherical structure for the spherical axis of the multi-region configuration is the triple $(\mathcal{Z},h,V)$, consisting of the spherical axis $\mathcal{Z}$, the height function $h$, and the multivalued vector field $V$. This depends upon the choice of an origin on which all height functions are [math].

From the spherical structure, we can reconstruct the boundary of $\mathcal{B}_{\infty}$ by $x=V(\mathbf{u})+h(\mathbf{u})\mathbf{u}$ for $\mathbf{u}\in\mathcal{Z}$ and the multiple values of $V$ at $\mathbf{u}$. Here $x$ denotes a collection of points corresponding to the values of $V$.

Then, the regions in ${\rm int}\,(\mathcal{B}_{\infty})$ are the regions in the complement of the boundary of $\mathcal{B}_{\infty}$ which have supporting hyperplanes for at least one point in one of the corresponding complementary regions to the spherical axis. If we have in addition the height function for the configuration defined on all of $S^{n}$, then we can construct the supporting hyperplanes for all $\mathbf{u}\in S^{n}$, and the envelope of these hyperplanes yields $\mathcal{B}_{\infty}$.

Blum Medial Linking Structure

We now consider the analogous Blum medial linking structure for a generic multi-region configuration. First, we note that if the configuration has regions with boundaries and corners, then the Blum medial axes of the individual regions will not define skeletal structures. This is because the Blum medial axis will actually meet the boundary at the edges and corners. However, in the case of disjoint regions $\{\Omega_{i}\}$ in ${\mathbb{R}}^{n+1}$ with smooth generic boundaries (which do not intersect on their boundaries) there is a natural Blum version of a linking structure, which we introduce.

Definition 4.6**.**

Given a multi-region configuration of disjoint regions $\boldsymbol{\Omega}=\{\Omega_{i}\}$ in ${\mathbb{R}}^{n+1}$, for $i=1,\dots,m$, with smooth generic boundaries (which do not intersect on their boundaries), a Blum medial linking structure is a skeletal linking structure for which:

• B1)

the $M_{i}$ are the Blum medial axes of the regions $\Omega_{i}$ with $U_{i}$ the corresponding radial vector fields;

• B2)

the linking axis $M_{0}$ is the Blum medial axis of the exterior region $\Omega_{0}$ (and we refer to it as the medial linking axis); and

• B3)

the $M_{i\,\infty}$ are the points in $\tilde{M}_{i}$ corresponding to points on $\mathcal{B}_{i}$ for which a height function has an absolute maximum (or minimum for the height function for the opposite direction).

Remark 4.7**.**

It follows from B2) that if $x\in M_{i}$ and $x^{\prime}\in M_{j}$ are linked, the corresponding values of the radial and linking functions satisfy $\ell_{i}(x)-r_{i}(x)=\ell_{j}(x^{\prime})-r_{j}(x^{\prime})$.

Generic Linking Properties

We consider a generic configuration of regions $\{\Omega_{i}\}$ with $\mathcal{B}_{i}$ the boundary of $\Omega_{i}$, and $M_{i}$ the Blum medial axis of $\Omega_{i}$. If each region is generic, then by Theorem 4.2 the Blum medial axis of each region has generic local structure. For 3D, the points on the boundary are of types $A_{3}$, $A_{1}^{3}$, $A_{1}^{4}$, $A_{3}A_{1}$ and $A_{1}A_{3}$, as shown in Figure 14, along with the remaining points of type $A_{1}^{2}$ corresponding to smooth points on the medial axis. For 2D, there are only isolated points of type $A_{3}$, corresponding to end points of the medial axis and $A_{1}^{3}$ points corresponding to $Y$-branch points, with the remaining points of type $A_{1}^{2}$.

A point of $\mathcal{B}_{i}$ may simultaneously have a role as both a point on the boundary of $\Omega_{i}$ and as a point on the boundary of the external region, see Figure 18. As such it has two distinct Blum descriptions, and hence simultaneously belongs to strata on the boundaries corresponding to given Blum behavior exhibited for the individual regions and for the external region. Generic linking between regions for a configuration satisfies special genericity properties which is explained by how these strata intersect.

We consider a collection of smooth boundary points $S=\{x_{1},\dots,x_{k}\}$ with $x_{i}\in\mathcal{B}_{j_{i}}$ associated to points $u_{i}$ on the medial axes $M_{j_{i}}$. Suppose the $u_{i}$ are linked at the point $u_{0}$ on the external medial axis $M_{0}$. In the generic case, this means $k\leq 3$ for ${\mathbb{R}}^{2}$ and $k\leq 4$ for ${\mathbb{R}}^{3}$. We use the notation $A_{\boldsymbol{\alpha}_{i}}$ to denote the Blum type of $x_{i}$ for the point $u_{i}$ for the individual region $\Omega_{j_{i}}$ and denote the corresponding strata of the boundary consisting of points of this type by $\Sigma^{\boldsymbol{\alpha}_{i}}_{\mathcal{B}_{j_{i}}}$. Likewise, we let $A_{\boldsymbol{\beta}}$ denote the Blum type for the point $u_{0}$ on the external medial axis, and denote the strata in the external medial axis by $\Sigma_{M_{0}}^{(\boldsymbol{\beta})}$ and the corresponding strata on the boundary by $\Sigma^{\boldsymbol{\beta}}_{\mathcal{B}}$. Then, the Blum linking type is denoted by $(A_{\boldsymbol{\beta}}:A_{\boldsymbol{\alpha}_{i_{1}}},\dots,A_{\boldsymbol{\alpha}_{i_{k}}})$.

Definition 4.8**.**

The set of points $S=\{x_{1},\dots,x_{k}\}$ as above exhibits generic Blum linking of type $(A_{\boldsymbol{\beta}}:A_{\boldsymbol{\alpha}_{i_{1}}},\dots,A_{\boldsymbol{\alpha}_{i_{k}}})$ if:

• i)

the strata $\Sigma_{\mathcal{B}}^{(\boldsymbol{\beta})}\subset\mathcal{B}$ and $\Sigma_{\mathcal{B}_{j_{i}}}^{(\boldsymbol{\alpha}_{i})}\subset\mathcal{B}_{j_{i}}$ intersect transversely in $\mathcal{B}_{j_{i}}$; and

• ii)

the images in $M_{0}$ of the strata $\Sigma_{\mathcal{B}}^{(\boldsymbol{\beta})}\cap\Sigma_{\mathcal{B}_{j_{i}}}^{(\boldsymbol{\alpha}_{i})}$ under the linking flow intersect transversely in the stratum $\Sigma_{M_{0}}^{(\boldsymbol{\beta})}$ (see e.g. Figure 19).

For example, generic linking in ${\mathbb{R}}^{2}$ means that zero-dimensional (isolated point) strata in one region can only be linked to points in the $A_{1}^{2}$ strata in another. Also, in ${\mathbb{R}}^{3}$ a one-dimensional stratum can intersect other one-dimensional stratum (transversely in a single point) in a boundary $\mathcal{B}_{i}$ and this can then also be linked to $A_{1}^{2}$ strata in another boundary as in Figure 19.

Example 4.9** (Generic Linking in ${\mathbb{R}}^{2}$).**

There are only five types of generic linking for configurations in ${\mathbb{R}}^{2}$ (see Figure 11):

• i)

$(A_{1}^{2}:A_{1}^{2},A_{1}^{2})$:  linking smooth points of two medial axes at a smooth point of the external medial axis, forming one-dimensional strata;

• ii)

$(A_{1}^{2}:A_{3},A_{1}^{2})$:  linking an end point of one medial axis with a smooth point of another at an isolated smooth point of the external medial axis, forming a zero-dimensional stratum;

• iii)

$(A_{1}^{2}:A_{1}^{3},A_{1}^{2})$:  linking a $Y$-branch point of one medial axis with a smooth point of another at an isolated smooth point of the external medial axis, forming a zero-dimensional stratum;

• iv)

$(A_{1}^{3}:A_{1}^{2},A_{1}^{2},A_{1}^{2})$:  linking smooth points of three medial axes at a $Y$ branch point of the external medial axis, forming a zero-dimensional stratum;

• v)

$(A_{3}:A_{1}^{2})$:  a single smooth point of a medial axis is self-linked at an end point of the external medial axis.

We note that all five types can occur for self-linking; but only the first four can occur for linking between distinct regions; and partial linking can only occur via iv).

By these being the only five generically occurring linking types, we mean that any other possibility, such as e.g. two end points of different regions being linked, may be arranged for a specially constructed configuration; however, almost all slight perturbations of this special configuration will destroy this special linking feature and only involve the five basic types.

For 3D, the situation becomes more complicated and there are seventeen generic linking types listed in Table LABEL:Table1 in §6. However, as explained in §7, by replacing the Blum linking structure with a simpler skeletal linking structure, we can significantly reduce the number of linking types.

Remark 4.10**.**

For a general multi-region configuration, we can substitute in place of $\Omega_{0}$ a region $\Omega_{i}$ which has multiple adjoining regions (including possibly the complement $\Omega_{0}$) and the definition of “generic linking”of the adjoining regions relative to $\Omega_{i}$ has the same form as in Definition 4.8. Generically they have the same properties as for $\Omega_{0}$.

Generic Structure for $\mathcal{B}_{\infty}$ and $M_{\infty}$

For a region $\Omega_{i}$ of a configuration, the points in $M_{i\,\infty}$ (and $\mathcal{B}_{i\,\infty}$) are not involved in linking; however, they do involve points which have a type $A_{\boldsymbol{\alpha}}$ based on the internal medial structure of $\Omega_{i}$. By the generic structure of $M_{i\,\infty}$ and $\mathcal{B}_{i\,\infty}$ we are interested in both the generic properties of the stratifications of $\mathcal{B}_{i\,\infty}$ and $M_{i\,\infty}$ resulting from the spherical axis and their relation with the strata $\Sigma^{\boldsymbol{\alpha}}_{\mathcal{B}_{i}}$.

The interior points of $\mathcal{B}_{i\,\infty}$ are those points where the supporting hyperplane meets the configuration at a single nondegenerate tangency. In addition, the boundary of $\mathcal{B}_{i\,\infty}$ consists of strata $\Sigma^{(\boldsymbol{\alpha})}_{\infty}$ corresponding to the strata $\Sigma^{(\boldsymbol{\alpha})}_{\mathcal{Z}}$ of the spherical axis $\mathcal{Z}$ and are of the same dimensions. The strata $\Sigma^{(\boldsymbol{\alpha})}_{\infty}$ lie in the smooth strata of the $\mathcal{B}_{i}$.

Definition 4.11**.**

By $M_{\infty}$ and $\mathcal{B}_{\infty}$ having generic structure we mean:

the spherical axis of the configuration is generic and each $\mathcal{B}_{i\,\infty}$ has strata with the resulting generic local structure given by Theorem 4.4;

• 2)

the strata of this stratification of $\mathcal{B}_{i\,\infty}$ intersect tranversally the strata $\Sigma^{(\boldsymbol{\beta})}_{\mathcal{B}_{i}}$ for the Blum medial axis of $\Omega_{i}$; and

• 3)

the strata of $M_{i\,\infty}$ are the images in $\tilde{M}_{i}$ of the transverse intersections $\Sigma^{(\boldsymbol{\alpha})}_{\infty}\cap\Sigma^{(\boldsymbol{\beta})}_{\mathcal{B}_{i}}$.

Examples of the generic local structure for the boundary strata of $M_{\infty}$ for a multi-region configuration in ${\mathbb{R}}^{3}$ are shown in Figure 16: a) of type $A_{1}^{2}$, b) of type $A_{1}^{3}$, and c) of type $A_{3}$. The darker curves (including the darker dashed curves) denote the boundary strata bounding regions of $M_{\infty}$ (consisting of points whose outward pointing normals point away from the other regions).

Existence of a Blum Medial Linking Structure

Before stating the general form of a full Blum linking structure for a general configuration, we first give a special case where the configuration $\boldsymbol{\Omega}$ consists of disjoint regions with smooth boundaries. Then, the existence of Blum medial linking structures is guaranteed by the following, which in addition ensures generic linking (see [Ga, Thms 3.4.3 and 3.4.4] and [DG, Thm 4.12]).

Theorem 4.12** (Existence of Blum Medial Linking Structure).**

For a generic configuration $\boldsymbol{\Omega}=\{\Omega_{i}\}$ in ${\mathbb{R}}^{2}$ or ${\mathbb{R}}^{3}$ consisting of disjoint regions with smooth boundaries (i.e. they do not intersect on their boundaries), and with $\boldsymbol{\Omega}$ contained in the interior of a given compact region $\tilde{\Omega}$,

1. (1)

the configuration has a Blum medial linking structure such that each $M_{i}$ (including $M_{0}$) has generic local properties given by Theorem 4.2; 2. (2)

the linking structure exhibits generic linking as in Definition 4.8; 3. (3)

$M_{\infty}$* and $\mathcal{B}_{\infty}$ have generic structure as given in Definition 4.11; and* 4. (4)

in the case that $\tilde{\Omega}$ is convex, the properties for a linking structure in the bounded case hold.

This result is a special case of the following general result for a general configuration allowing adjoined regions in [DG, Thm 4.13].

Theorem 4.13** (Full Blum Linking Structures).**

A generic multi-region configuration $\boldsymbol{\Omega}=\{\Omega_{i}\}$ in ${\mathbb{R}}^{2}$ or ${\mathbb{R}}^{3}$ contained in the interior of a compact region $\tilde{\Omega}$ has the following generic properties:

• i)

each $\Omega_{i}$ has a Blum medial axis $M_{i}$ exhibiting the generic local properties at interior points of $\Omega_{i}$ given by Theorem 4.12;

• ii)

the external region in $\tilde{\Omega}$ has a medial linking axis $M_{0}$ which exhibits the generic local Blum properties;

• iii)

the local structure of the Blum axes $M_{i}$ (including $i=0$) near a singular boundary point has the local generic edge-corner normal form given by Definition 4.1;

• iv)

at a smooth boundary point of a region $\Omega_{i}$ of type $Q_{k}$, $k=2,3$, the strata of $Q_{k}$-points transversally intersects the strata of Blum type points $\Sigma_{\mathcal{B}_{i}}^{(\boldsymbol{\alpha})}$;

• v)

generic linking occurs between the smooth points* of the regions and no linking occurs at edge-corner points; and*

• vi)

$\mathcal{B}_{i\,\infty}$* is contained in the smooth strata of the $\mathcal{B}_{i}$, and $\mathcal{B}_{i\,\infty}$ and $M_{i\,\infty}$ exhibit the generic properties given in Definition 4.11.*

Note: Property v) holds as well for generic linking between adjoining regions of a given region $\Omega_{i}$ relative to the region $\Omega_{i}$.

5. Modifying the Full Blum Medial Structure to a Skeletal Linking

Structure

We know by Theorem 4.13 that a generic multi-region configuration has a Blum medial structure. If the regions are disjoint with smooth boundaries, then the Blum linking structure is a skeletal linking structure. However, if the configuration contains regions which adjoin, then the Blum linking structure does not satisfy all of the conditions for being a skeletal linking structure. Specifically the individual Blum medial axes of both the regions and the complement will extend to the singular points of the boundaries.

There are two perspectives on this. On the one hand, as mentioned in Remark 4.3, we may view this as a “relaxed form of a skeletal linking structure.”We shall see that from this structure we still obtain all of the local, relative, and global geometry of the individual regions and the positional geometry of the configuration. However, if we consider the stability and deformation properties, such a structure does not behave well.

We describe two approaches to modifying the full Blum linking structure to obtain a skeletal linking structure. One approach is when the configuration with adjoined regions can be viewed as a deformation of a configuration with disjoint regions. The second is to modify the full Blum linking structure by a process of “smoothing the corners”of the regions. We describe each of these.

Example of Evolving Skeletal Linking Structure for Simple

Generic Transition

We illustrate the method for a simple generic transition under which two initially disjoint regions $\Omega_{i\,0}$, $i=1,2$, undergo a deformation $\Omega_{i\,t}$ to become adjoined at $t=t_{0}$ and form a common boundary region $Z_{t}$ for $t>t_{0}$. This is illustrated in Figure 20.

However, this causes a transition in the full Blum medial linking structure as in changing from a) to b) in Figure 21. There is a discontinuous change in the Blum medial axis somewhat analogous to the introduction of a $C^{1}$ bump in the boundary of a region which forces the Blum medial axis to introduce a new branch.

An alternate approach is to deform the Blum medial linking structure of the regions before the transition to a skeletal linking structure as the transition occurs and continues for the adjoined regions to evolve from a) to c) in Figure 21. This is achieved by keeping the skeletal axes for each region but altering the lengths of both the radial vectors so they extend to the common boundary and altering the linking vectors to the external linking axis as shown.

This gives a family of skeletal linking structures for the varying configuration $\boldsymbol{\Omega}_{t}^{\prime}=\{\Omega_{1\,t},\Omega_{2\,t}\}$, which evolve continuously (and stratawise differentiably on the added strata $S_{i\,t}$), see Figure 21.

Remark 5.1**.**

This example illustrates one significant advantage of skeletal linking structures over the full Blum linking structure for general multi-object configurations; namely, certain changes in configuration type can be incorporated as continuous variation in the skeletal linking structure.

Modifying the Full Blum Linking Structure via Smoothing

A second approach for a configuration with adjoining regions is to modify the full Blum linking structure to a skeletal linking structure. One way to accomplish this is by “smoothing of the corners of the regions”in a small neighborhood of the edge/corner set as in Figure 22. The goal is to modify the regions $\Omega_{i}$ in a small neighborhood $W$ of the edge corner points so that the smooth region boundaries are smooth, interior to $\Omega_{i}$ and agree with the original boundaries outside of $W$. This is done in such a way that the radial vector field of the Blum medial axes of the smooth regions can be extended to be transverse to the original boundary and can be extended to meet the external medial axis transversely. These extended vectors are then the linking vectors for a skeletal structure. The details of this can be found in [DG, Chap. 5].

6. Classification of Linking Types for Blum Medial

Linking Structures in ${\mathbb{R}}^{3}$

For a generic configuration $\boldsymbol{\Omega}$ in ${\mathbb{R}}^{3}$, the Blum medial linking structure exhibits generic linking properties given in Table LABEL:Table1. We first briefly explain the features of the table. Dimension refers to the dimension in ${\mathbb{R}}^{3}$ of the strata where the given linking type occurs. There are three types of linking: 1) linking between points on distinct medial axes; 2) “partial linking”  involving more than one point from one medial axis and point(s) from another; and 3) “self-linking”  where linking is between points from a single medial axis. “Pure linking type”  refers to cases only occurring for self linking. The $A_{1}^{2}$ and $A_{1}A_{3}$ linking can occur for either linking or self-linking; while $A_{1}^{3}$ and $A_{1}^{4}$ linking can occur for any of the three linking types.

7. Simplifying the Blum Linking Structure for Medical Imaging

For modeling configurations in 2D and 3D medical images, the strategy we propose builds upon that used in modeling single objects using M-reps and S-reps as developed by Pizer and coworkers at UNC MIDAG. The Blum medial axis is replaced by a simpler skeletal structure where the skeletal set (also called the medial set) is replaced by a surface with boundary diffeomorphic to a $2$-dimensional disk, which is modeled by a rectangular grid with a multi-valued vector field at the vertices of the grid. Then, the grid together with vectors is the template which is fitted to a specific object in each of a training set of medical images. The size of the grid is experimentally chosen based on the object. The fitting involves an optimization process based on several criteria involving nonsingularity of the structure, regularity of the interpolated vector fields, nonsingularity of the constructed object boundary, and closeness of fit to an expert segmentation of the object. These criteria use in several ways the underlying mathematics of skeletal structures to numerically measure the criterion and contribute to the interpolation technique.

The fitting of the templates provides a discrete set of data for each object in the training set; but the data are points on a product manifold. A form of nonlinear PCA is applied to identify the most significant feature directions. Then an iterative process is applied to further refine the fits and ultimately to yield a fitting process which is applied to new images using the preceding criteria and the statistical priors. This method provides extremely high quality fit yielding high segmentation accuracy for a number of different cases that have been studied, and led to commercial software for imaging the male pelvic region.

In several cases mentioned earlier where several objects in the same image are considered, additional user chosen positional information is attached which improves the positioning of the templates to the individual objects. What the work here provides for a configuration of objects is an analysis of how the entire collection can be simultaneously modeled beginning with a Blum linking structure. Because of several weak features and disadvantages of the Blum medial axis, this linking structure would be replaced by a more general but simplified skeletal linking structure. This linking structure would still restrict to give an S-rep representation for each object. By using skeletal sets without singular points (except the edge points), the number of linking types are reduced from $17$ to $9$. Moreover, if the external linking axis can be simplified by, for example, avoiding self-linking, then this can be further reduce to $7$ types in Table LABEL:Table1: i), iii), viii), x), xii), xiii), and xiv). These concern the singular structure of the external linking axis and the relation of the edges of the skeletal sets with each other and the singular points of the linking axis. The added data for a fitted linking structure, beyond that for the S-rep data for each object, would include the values of the linking functions at the vertex points, and the points on the skeletal sets corresponding to the linking types. The linking function values are easily included in the data, while for the linking types more investigation is needed to effectively include these so the statistical analysis includes all information.

Right now this has begun on a region in the neck which includes the three objects illustrated in Figure 1 b). There is still needed the detailed fitting for a training set with experimental decisions made on the size of the grids for the S-reps to also incorporate the extra data. This will involve the use of additional mathematical structure built upon that for skeletal structures, which is explained in the related paper [DG2]. Then, the process will proceed to arrive at a coherent template for the entire configuration and a specific optimization procedure for the fitting.