Medial Axis Isoperimetric Profiles

While there is no known efficient algorithm to compute the full isoperimetric profile, much is known about $\partial E_{\Omega}(t)$. In particular, for any fixed value of $t$, the free boundary of $E_{\Omega}(t)$ , $\partial E_{\Omega}(t)\setminus\partial\Omega$ must have constant and nonsingular mean curvature [GMT80]. In addition, $\partial E_{\Omega}(t)$ and $\partial\Omega$ must meet tangentially, and any free boundary of $E_{\Omega}(t)$ must be a circular arc [BES49]. For a small subset of values of $t$, $E_{\Omega}(t)$ and $\mathrm{IP}_{\Omega}(t)$ are known in closed form [LS19]. We will expand on this in §3.3 after establishing more notation and geometric tools.

Literature on the isoperimetric profile has typically focused on suitably regular Jordan domains.

Beyond restricting attention to Jordan domains, there is typically an additional regularity constraint that $\mathrm{area}(\partial\Omega)=0$, where $\mathrm{area}(\cdot)$ denotes the Lebesgue 2-dimensional measure. Domains satisfying this condition include any finite polygon and Jordan domains with smooth boundary. For contrast, Osgood curves are Jordan curves that can have nonzero area [Kno17], but for our purposes these will not be considered. The set of suitably regular Jordan domains will be denoted $\mathbb{J}$.

We will further refine these domains $\Omega$ by the characteristics of their necks.

Let $r_{n}$ denote the smallest value of $\rho$ for which $\Omega$ has a neck. This definition is built intuitively on the idea of a bottleneck where objects of a certain size are unable to pass through. A domain that does not satisfy Definition 3.3 for any $\rho$ has no neck, i.e., $r_{n}=\infty$. For example, star domains have no necks. In contrast, it is clear that the shape in Figure 3 has $r_{n}<\infty$. This definition does not have a concept of multiple necks nor does it take into account the locations of necks. We will provide a new definition in §6.1 that generalizes necks to incorporate these concepts.

Natural variations of a shape are given by morphological operations. Furthermore $E_{\Omega}(t)$ for limited values of $t$ can be characterized by the morphological opening of $\Omega$. We will also use morphological closing in §7.2 to clean up our results. For these reasons, we provide their definitions; see [Ser00] for more discussion.

Let $B_{r}(x)$ be a ball of radius $r$ centered at $x\in{\mathbb{R}}^{2}$ and $B_{r}\coloneqq B_{r}(0)$. Given a starting domain $\Omega\subset{\mathbb{R}}^{2}$, its erosion by $B_{r}$ is

Its dilation by $B_{r}$ is

Its opening by $B_{r}$ is

Finally, its closing by $B_{r}$ is

These operations and their results are depicted in Figure 3. Intuitively, the opening can be thought of as a blurring of the convex corners of a shape. The closing can be thought of as a way of filling in small holes or smoothing out concave corners of a shape.

With this notation, we recall a property of the isoperimetric profile, namely that $E_{\Omega}(t)=\Omega\circ B_{r}$ when $r<r_{n}$ with corresponding area $t=\mathrm{area}(\Omega\circ B_{r})$ [LS19]. It follows that when restricted to no neck domains, the isoperimetric profile can be computed by opening [LS19, Theorem 2.3]. For domains with necks, however, the opening is not generally equal to $E_{\Omega}(t)$. This is exemplified in Figure 3 where the opening has singular curvature near the neck. Furthermore, we find that in practice the range $t\in[\mathrm{area}(\Omega\circ B_{r_{n}}),\mathrm{area}(\Omega)]$ accounts for a very small portion of the isoperimetric profile on domains with necks.

One of our key ideas is to connect the widely used medial axis to the isoperimetric profile. The medial axis allows us to formulate a new definition of neck that will be helpful for our analysis.

The opening of a 2D domain is simple to compute from modern geometry processing tools and generates feasible solutions to problem (IP) on general domains. Unfortunately, the IP upper bound generated by measuring the perimeter of $\Omega\circ B_{r}$ exhibits several undesirable traits on domains with necks. Firstly, $\Omega\circ B_{r}$ varies discontinuously with $r$, leaving large gaps in the profile; we repair this issue in §7.2. Second, this IP bound is not monotonic, and we know the true isoperimetric profile should always monotonically increase. Finally, $\partial(\Omega\circ B_{r})$ can have singular curvature on free boundaries, while $\partial E_{\Omega}(t)$ cannot. For these reasons direct computation of $\Omega\circ B_{r}$ does not extend well to domains with necks. To fix these shortcomings, we propose a novel medial axis based approach.

Recall from §3.4 that the medial axis transform is invertible i.e. we can reconstruct $\Omega$ from $\mathrm{MA}(\Omega)$ and its radius function $r(x)$. Unsurprisingly, performing the reconstruction on a subset of $\mathrm{MA}(\Omega)$ will produce a subset of $\Omega$. We will refer to these as limited reconstructions.

When $g=\mathrm{MA}(\Omega)$, $\Omega_{g}=\Omega$. It turns out that for a particular choice of $g$, we can relate the morphological opening of a domain to its limited reconstruction.

For proof, see §LABEL:proof:medaxrecon. The relationship between morphological operations and the morphological skeleton has been studied in image processing [Lan77], where a similar statement is made in the discrete setting. Recall that the isoperimetric profile of no neck domains in $\mathbb{J}$ can be constructed by morphological opening. By Proposition 1 these can now be constructed from subsets $g$ of the medial axis. Furthermore, $g$ can be built iteratively due to the following result:

This motivates our algorithm in §5 for constructing the isoperimetric profile by starting with $g_{\mathrm{inr}(\Omega)}$ and iteratively increasing $g$ by greedily absorbing neighboring medial axis points of highest radius. We will see in §7 why this medial axis based approach far outperforms direct computation of $\Omega\circ B_{r}$.

Since our algorithm will traverse the medial axis, we consider here how much the isoperimetric profile of $\Omega_{g}$ will change as $g$ grows differentially. We can simplify this problem by looking at only a single medial axis edge segment $\gamma:[0,1]\rightarrow{\mathbb{R}}^{2}$, and let $g(\tau_{2})=\bigcup_{\tau\in[\tau_{1},\tau_{2}]}\gamma(\tau)$, for $\tau_{1},\tau_{2}\in[0,1]$ and $\tau_{1}<\tau_{2}$. We derive that

For derivation see “differentialGrowth.nb” in supplemental materials. Equation (7) tells us that if we infinitesimally increase $g$ by absorbing vertices of the medial axis near $x$, the slope of the isoperimetric profile will be $\nicefrac{{1}}{{r(x)}}$. A consequence of Equation (7) is that the differential growth of the isoperimetric profile when growing a circle of fixed center from radius 0 is infinite, i.e., differentially, it is costly for $F$ in problem (IP) to be disconnected.

For comparison with the direct computation of $\Omega\circ B_{r}$ suggested in §4.1, we pick 100 evenly spaced values of $r_{i}\in[0,\mathrm{inr}(\Omega)]$, compute $\Omega\circ B_{r_{i}}$ and its corresponding area and perimeter. This is easily implemented with Matlab’s polybuffer methods and forms a second upper bound that in the majority of test cases performs worse than our medial axis based approach. In very few cases, it can perform better, in which case our final IP upper bound is the minimum of these two upper bounds.

Since our algorithm is derived jointly from the medial axis and the concept of necks, we create a generalized definition for necks based on the medial axis transform. From this definition, we can quantify the size of multiple necks as well as their locations (see Figure 6). This will allow us to prove results that are not possible to state with Definition 3.3.

Algorithmically, the computation of the necks $\eta(\Omega)$ is done by first computing the medial axis, and then iterating over its edges and nodes while checking the neck conditions. We denote the maximum neck radius by $r_{m}=\max\{r(x):x\in\eta(\Omega)\}$. Similarly we denote its minimum neck radius by $r_{n}=\min\{r(x):x\in\eta(\Omega)\}$. If a domain has no necks then we say $r_{n}=\infty=r_{m}$. The maximum neck radius $r_{m}$ is not part of the previous definition of necks but will be used in Proposition 4 to give theoretical guarantees about our algorithm. A nice property of our neck definition is that it is compatible with the previous neck definition, i.e., the minimal neck radius of both are identical.

For proof see §LABEL:proof:neckneck. Note that this equivalence does not hold for domains outside $\mathbb{J}$. Consider the annulus with a pinch in it, illustrated in Figure 6. This does not have a neck by Definition 3.3, but by Definition 6.1 it has a single neck marked by the black point.

While much is known about the isoperimetric profile for shapes with no necks ($r_{n}=\infty$), much less is known when $r_{n}$ is finite. One might hope that if $r_{n}$ were in some sense large enough, we could still compute the isoperimetric profile with some confidence.

Here we define the thick neck condition which will be used in Proposition 5 to analyze the behavior of our algorithm on domains whose necks are large enough. Recall that $\mathrm{inB}(\Omega)$ is the largest circle inscribed in $\Omega$.

For examples of domains that satisfy the thick neck condition and their profiles, see Figure 12.

By construction, Algorithm 1 produces an upper bound to the isoperimetric profile. This is because $\overline{F_{i}}$ is always a feasible candidate to (IP) with $t_{i}=\mathrm{area}(\overline{F_{i}})$. By applying Equation (7) to domains in $\mathbb{J}$, our bound is strictly monotonically increasing, similarly to the actual IP. Our construction also satisfies the condition in [GMT80], which states that $\partial E_{\Omega}(t)$ must lie tangent to $\partial\Omega$ and its free boundary must have no points of singular curvature. This is because the free boundary of $\overline{F_{i}}$ is made of circular arcs that start and end on $\partial\Omega$ and whose centers are on different edges of the medial axis. If these circular arcs intersect to form a point of singular curvature, then a non-contractable loop can be drawn in $\Omega$ by tracing from the intersection through both edges of the medial axis, back to $\mathrm{inx}(\Omega)$, implying $\Omega\notin\mathbb{J}$.

We provide guarantees for when our bound is tight for any domain in $\mathbb{J}$. While we cannot guarantee that our bound is tight for all $t_{i}$, by Proposition 4, we know our bound is tight for $t_{i}$ in a limited range.

While Proposition 4 does not cover the entire IP, it applies to any domain in $\mathbb{J}$ regardless of the presence or size of necks. To achieve this result, this proposition uses $r_{m}$, a quantity that is not defined in prior work.

Application of Proposition 4 is demonstrated in Figure 7. We compute our upper bound to the IP by Algorithm 1 on the Candy domain. Its profile is depicted between two events: the area of the “max inscribed circle” $\mathrm{area}(\mathrm{inB}(\Omega))$, and the “max area” $\mathrm{area}(\Omega)$. We apply Proposition 4, which tells us that our medial axis bound for the IP is tight to the left of the “Conservatively tight” purple event line and to the right of the “Minimal Neck” purple event line. For values of $t_{i}$ between the two purple lines, we upper-bound the IP.

Given some assumptions on the domain $\Omega$, we show that our IP bound exactly recovers the IP.

Intuitively, this is because the smallest perimeter cost for having a disconnected $E_{\Omega}(t)$ is larger than the cost for continuing to grow $E_{\Omega}(t)$ via limited reconstructions. For proof, see §LABEL:proof:bigneck. Proposition 5 allows us to compute the exact isoperimetric profile for domains satisfying the thick neck condition. We show in §8.3 examples of domains satisfying this condition and that while previous methods of computing the isoperimetric profile are not tight, our method provides the exact profile for this larger class of domains.

Following Equation (8) we compute $\Omega_{g}$ by the union of floating precision polygons. This sometimes introduces artifacts like holes, or jagged boundary at a scale that is invisible. While seemingly negligible, many holes or jagged regions add up to significantly increase the perimeter thus loosening our upper bound. This is solved easily by some basic post processing. First we remove all holes in $F$, then we perform morphological closing with a ball of small radius, and lastly we intersect it with $\Omega$ to make it feasible again. These steps are summarized in Algorithm 2. Figure 8 shows the change in $\overline{F}_{i}$ before and after post processing. Figure 9 demonstrates the difference in the IP with and without post processing.

We also perform a repair step to the upper bound obtained by direct computation of $\Omega\circ B_{r}$. Since $\Omega\circ B_{r}$ changes discontinuously with $r$ for domains with necks, this leaves gaps in the IP. An example of this discontinuous transition is visualized in Figure 9. The discontinuous change happens due to the creation of a new disconnected point $x^{-}$ in $\Omega\ominus B_{r}$. $\Omega\circ B_{r}$ will instantly include a circle centered on $x^{-}$ of radius $r$, which abruptly increases area and perimeter. We detect when the number of disconnected regions in $\Omega\ominus B_{r}$ changes, and compute $x^{-}$ as the center of the bounding box of the smallest area component of $\Omega\ominus B_{r}$. We then inflate a disk centered at $x^{-}$ from radius $0$ to $r$. This procedure continuously fills in the IP bound and is illustrated in Figure 9.

We compare our algorithm against the few existing computational approaches to the isoperimetric profile. First we compare Algorithm 1 with direct morphological opening. While both generate upper bounds, direct opening has much fewer guarantees and generally produces a looser bound. Next we compare with the TV approach [DLSS19] which generates lower bounds. In theory, it generates the convex lower envelope of the IP.

Code for our algorithm, figure generation, and all parameters used can be found in supplementary materials. We show profiles starting from the area of the maximum inscribed circle because the profile before then is uninformative.

On no neck domains both the morphological opening procedure and our medial axis traversal will produce the exact isoperimetric profile. This is depicted in Figure 11. We test on a star domain Star as well as a non-star, no neck domain no-neck-J. We choose $r_{l}=.1$ which computes most of, but not the entire, profile. This is why our profile stops before the maximum area is reached. The total variation lower bound roughly follows the convex lower envelope of our profile, but overshoots in a few regions. As a lower bound, it leaves slack room, while our algorithm is exact.

Note that the choice of $r_{l}=.1$ is specific to the domains we tested on. As $r_{l}$ roughly corresponds to the smallest length scale of a domain our IP bound will capture, a smaller domain requires a smaller $r_{l}$.

On thick neck domains our algorithm computes the exact IP, while morphological opening computes a much looser upper bound, and total variation computes the convex lower envelope of the profile. This is demonstrated in Figure 12 on the Jellyfish, Hourglass, and Worm domains. Note that the “Minimal Neck” area, $\mathrm{area}(\Omega\circ r_{n})$, is often very close to the maximum area. Prior methods would only be able to compute the exact profile in the limited region from $\mathrm{area}(\Omega\circ r_{n})$ to $\mathrm{area}(\Omega)$, while our method computes the entire profile. Due to the prescence of necks, direct morphological opening produces much worse upper bounds than our method.

We test on a set of geographic and hand drawn data. Maps of Alabama, California, Delaware, and District 1 of Alabama were obtained from [FIP20]. A map of Mozambique was obtained from [Vem20]. In addition, we create a drawing, Boxes, to test on. Alabama, California, and Delaware satisfied the thick neck condition while District 1, Mozambique, and Boxes do not. Our profiles for these domains are shown in Figure 13. Our bound computes the exact IP for all $t$ on Alabama, California, and Delaware. [DLSS19] continues to provide convex lower envelopes with slack and an occasional intersection with the upper bound due to imperfect discretization. This lower bound is particularly loose, on Alabama, and California revealing a significant gap for most of the profile. Morphological opening produces upper bounds as well, but are loose compared to our medial axis based approach and has no theoretical guarantees except for a small region on the tail of the profile.

On the District 1, Mozambique, and Boxes domains, we use Proposition 4 to compute “Conservatively tight” and “Minimal Neck” events that allow us to guarantee our bound is tight on a larger range of $t$ values. For both District 1, and Mozambique, our algorithm produces much tighter upper bounds than morphological opening. An exception occurs on Boxes, where the morphological opening produces a tighter bound. This is explained by the fact that the optimal $E_{\Omega}(t)$ is disconnected for much of the profile, while our algorithm by construction only produces connected polygons. For this reason it would be worthwhile to explore using disconnected medial axis reconstructions based on a criteria like the thick neck condition. Further exploration of disconnected limited reconstructions is discussed in §LABEL:sec:discuss and left to future work.

By combining the bounds and guarantees provided in this paper and those of [DLSS19, LS19], we compute the envelope inside which the isoperimetric profile must lie. This is visualized in Figure 14. We find that this produces a relatively narrow region of uncertainty where the IP is unknown. This visualization also makes clear the rich space of multi scale behaviors that different domains can have.

For example, consider the Jellyfish (brown). It consists of a large cap region which has mostly low resolution features, followed by several tendrils that require high resolution data to capture. This is reflected in its isoperimetric profile where it starts shallow for most of the area, before sharply veering up to have one of the largest perimeters in this figure. Delaware (bright green) demonstrates similar behavior to a lesser extent. It starts with a shallow profile indicating its coarse lower portion, before its profile curves upwards to capture the high resolution details of its upper quarter. Delaware requires less high resolution data to capture than the Jellyfish, however, as indicated by its profile sitting below that of the Jellyfish. Our results with the isoperimetric profile make it possible to computationally and numerically compare these domains’ multi-scale properties.

A list of runtimes for our experiments is provided in Table LABEL:table:timings. The runtime of our algorithm depends most strongly on number of vertices $|V|$, and loosely on other parameters $d_{c}$, $r_{l}$, and $r_{n}$. For morphological opening, the runtime depends primarily on the number of samples of the IP that are computed. For all times listed, we sample 100 values of the IP. For TV, runtime depends on the number of samples of the IP that are computed as well as the pixel resolution used to rasterize the polygon. We compute 20 samples of the TV IP and use 100 vertical pixels for the rasterization.