A Connected Component Labeling Algorithm for Implicitly-Defined Domains

In this paper, we develop a connected component labeling algorithm for implicitly-defined domains specified by multivariate polynomials. In particular, we consider domains of the form $\Omega:=\mathcal{U}\setminus\{\phi=0\}$, where $\phi:\mathbb{R}^{d}\to\mathbb{R}$ is a given polynomial and $\mathcal{U}\subset\mathbb{R}^{d}$ is a given bounded $d$-dimensional hyperrectangle, and consider the problem of implementing a labeling function $\chi:\Omega\to\mathbb{N}$ such that $\chi(x)=\chi(y)$ if and only if the two points $x,y\in\Omega$ are path-connected.

Our interest in this labeling problem stems from the author’s prior work on quadrature algorithms for multi-component implicitly-defined domains [18]. For a general, piecewise-smooth integrand function $f$, these algorithms output a quadrature scheme of the form $\int_{\mathcal{U}}f\approx\sum_{i}w_{i}f(x_{i})$ and have the following key property: if $V$ is a connected component of $\mathcal{U}\setminus\{\phi=0\}$ and $f$ is sufficiently smooth on $V$, then $\int_{V}f\approx\sum_{x_{i}\in V}w_{i}f(x_{i})$ represents a high-order accurate quadrature scheme on $V$. In other words, to build a quadrature scheme on $V$, one may simply discard the quadrature points outside of $V$, leaving the weights of the remaining points unmodified. Consequently, one may think of the bulk quadrature scheme over $\mathcal{U}$ as an agglomeration of smaller quadrature schemes over the individual connected components of $\mathcal{U}\setminus\{\phi=0\}$. The algorithms developed in [18], however, do not automatically determine this grouping; one of the aims of the present work is to design a simple and efficient algorithm for this purpose. Besides this application, connected component analysis arises in a wide variety of different settings, including in computer-aided geometric design, shape modeling and pattern matching, and in building topological feature descriptors of physical data sets [12].

A number of different approaches could be taken to solve the connected component labeling problem:

• •

  One approach is to apply the workhorse tools of real algebraic geometry, such as the cylindrical algebraic decomposition (CAD) algorithm pioneered by Collins [6, 5], or by computing roadmaps [2, 3]. In exact or arbitrary-precision arithmetic, CAD and roadmap methods can be used to identify and analyze the connected components of general, arbitrarily complex semi-algebraic sets (of which $\Omega$ is one particular instance). However, these methods can be computationally expensive; moreover, their implementation in fixed-precision arithmetic requires considerable care, even for relatively simple geometry.

• •

  A second approach is to apply the computational methods of Morse theory, Reeb graphs, and contour trees [11, 12]. Roughly speaking, these methods compute the topology of the set of level sets of $\phi$ via the location of its extrema and saddle points, combined with gradient path tracing algorithms. These methods become increasingly more intricate to implement as the dimension $d$ increases.

• •

  A third approach might be to isolate individual components of $\Omega^{+}:=\mathcal{U}\cap\{\phi>0\}$ by finding a polynomial which can separate them (followed by a similar procedure for $\Omega^{-}:=\mathcal{U}\cap\{\phi<0\}$). For example, if we could find a polynomial $u:\mathbb{R}^{d}\to\mathbb{R}$ such that $\Omega^{+}\cap\{u=0\}$ is empty and $u(x)u(y)<0$ for two sample points $x,y\in\Omega^{+}$, it would follow that the zero level set of $u$ forms a kind of divider that separates the two connected components associated with $x$ and $y$. To find this polynomial, one could use polynomial positivstellensatz or sums-of-squares methods [13, 14] to compute, if possible, a polynomial $w$ such that: (i) $w>0$ on $\Omega^{+}$; (ii) $w=u^{2}$; and (iii) $u(x)u(y)<0$. Part (i), (ignoring conditions (ii) and (iii) and various other subtleties,) can be solved efficiently via linear or semidefinite programming techniques, see, e.g., [14, 1]; however, conditions (ii) and (iii) end up creating a quadratically-constrained quadratic program for the coefficients of $w$, whose associated feasibility condition involves a challenging non-convex constraint. These ideas, among various others for computing separating polynomials, were investigated as part of the present work; however, no sufficiently elegant or efficient approach was found.

• •

  Finally, a fourth and considerably simpler approach is to subdivide the constraint domain $\mathcal{U}$ into small enough cells for which the topology of $\Omega$ is easy to determine. Depending on the tessellation method, this approach can be quite effective (fast, simple to implement, high resolution power), but in some cases the subdivision process may need to be stopped, potentially leading to some components being incorrectly merged or broken.

The preceding discussion serves to highlight that, depending on the final application, one must ultimately decide on a suitable compromise between: (i) absolute correctness/certifiability of the connected component labeling algorithm, potentially needing costly arbitrary-precision computation (as exemplified by CAD-based methods); (ii) robustness, e.g., in suitably handling fixed-precision roundoff errors or uncertainty in the input polynomial coefficients; (iii) computational complexity of constructing the labeling function as well as its subsequent evaluation; and (iv) implementation simplicity. The target application of this work concerns the development of numerical methods for multi-physics simulations involving highly complex geometry such as liquid atomization dynamics; in particular, the input polynomial $\phi$ represents the fluid interface and is computed dynamically. This application necessitates prioritizing aspects (ii), (iii), and (iv), which, in turn, necessitates a design choice for how to handle edge cases such as nearly touching components or interfacial self-intersections. Our choice here is to require that components are never broken, i.e., if $x,y$ are path-connected in $\Omega$ and $\chi(x),\chi(y)$ is the output of the labeling algorithm, then $\chi(x)=\chi(y)$ is an absolute certainty. It follows that the labeling algorithm must be permitted to merge or “glue” together distinct components in the (presumably rare) cases of uncertain topology. The algorithm developed in this work glues components only if they are exactly or nearly touching, relative to a user-defined gap threshold, typically orders of magnitude smaller than the length scale of $\mathcal{U}$.

With this design objective in mind, the connected component labeling algorithm developed here follows the fourth approach mentioned earlier: it operates by recursively subdividing the constraint hyperrectangle $\mathcal{U}\subset\mathbb{R}^{d}$ into hyperrectangular subcells until the topology of $\Omega$ thereon is sufficiently simple or has reached a smallest-permitted size. This is achieved through a combination of quadtrees (in 2D) or octrees (in 3D), along with simple yet effective topology tests using properties of Bernstein polynomials. Once the tree is constructed, its leaf cells form the vertices of a graph whose edges are determined by sign attainability tests on the faces of adjacent leaf cells. The connected components of this graph ultimately decide the overall labeling: $\chi(x)$ equals the label of the leaf cell containing $x$, the latter found by fast tree traversal methods. We show in this paper that the overall approach is efficient in handling various degrees of geometric complexity. In particular, if the recursive subdivision process terminates without reaching the user-defined smallest subcell size, then the labeling is provably correct. Applied to a fixed class of randomly generated geometry, the percentage of cases in which the algorithm is uncertain about the topology decreases exponentially as a function of the maximum tree depth; moreover, only a small fraction of these cases yield artificially-glued components.

An outline for the rest of the paper is as follows. In §​ 2, we establish some preliminaries on Bernstein polynomials and define the simply connected topology test. §​ 3 presents the main connected component labeling algorithm, followed by a discussion of its features. Numerical tests are presented in §​ 4, analyzing the algorithm’s success and failure rates on randomly generated multi-component geometry in 2D and 3D as well as its behavior on particular examples exhibiting cusps, self-intersections, and other kinds of singularities. Concluding remarks with a brief discussion of possible extensions are given in §​ 5.

A key part of the connected component labeling algorithm is an effective means for evaluating the range of attainable values of a polynomial. One of the most accurate and straightforward methods for doing so is through the use of the Bernstein basis [8, 15, 10, 4]. These methods make use of a convex hull property and guarantee that a polynomial’s value, at any point in its rectangular reference domain, is no larger (or smaller) than its maximum (or minimum) coefficient in the Bernstein basis. Especially useful in the present setting, these bounds become monotonically more accurate under the stable operations of Bernstein subdivision. In addition, owing to the hyperrectangular constraint domain $\mathcal{U}$ and its subdivision into hyperrectangular subcells, it is particularly natural to use a tensor-product basis. Accordingly, a tensor-product Bernstein basis is adopted throughout this work.

In $d$ dimensions, let $\phi$ be a tensor-product Bernstein polynomial of degree $n=(n_{1},\ldots,n_{d})$, defined relative to the hyperrectangle $U=[\alpha_{1},\beta_{1}]\times\cdots\times[\alpha_{d},\beta_{d}]$; $\phi$ takes the form

where

are the one-dimensional Bernstein basis functions of degree $\eta$ relative to the interval $[\alpha,\beta]$, $c_{i}=c_{i_{1},\ldots,i_{d}}$ denotes the $i$th Bernstein coefficient of $\phi$ for a multi-index $i\in{\mathbb{N}}^{n}:=[0,n_{1}]\times\cdots\times[0,n_{d}]$, and $\smash{b_{i}^{n,U}}(x)=\smash{\prod_{j=1}^{d}}\smash{b_{i_{j}}^{n_{j},[\alpha_{j},\beta_{j}]}(x_{i_{j}})}$ denotes the product of basis functions. This representation gives the Bernstein coefficients of $\phi$ relative to the given hyperrectangle $U$; subdivision refers to transforming these coefficients relative to a subset hyperrectangle and can be stably computed using the de Casteljau algorithm [7].

In the connected component labeling algorithm, we need a method to assess whether the topology of $U\cap\{\phi>0\}$ and $U\cap\{\phi<0\}$ is simple enough. This is achieved through the following three concepts.¹¹1Definition 2.1 has, in various guises, appeared before in the literature; the remaining two concepts (Definitions 2.2 and 2.3, along with their implications) are perhaps new, and target specifically the objectives of the present work.

An equivalent viewpoint comes from a convenient property of differentiation in the Bernstein basis: up to a multiplicative factor, the coefficients of the derivative are formed via first-order divided differences of the input polynomial’s coefficients. Therefore, $\phi$ is coefficient monotone if and only if the Bernstein coefficients of its derivative in the corresponding direction are either all non-negative or all non-positive. It follows that coefficient monotone polynomials are monotone in the conventional sense and therefore attain their extrema on the corresponding faces of the hyperrectangle. Tests of Bernstein coefficient monotonicity have a variety of applications, including, e.g., in polynomial range evaluation [16]. In the present application, an important property is that, for a coefficient monotone polynomial $\phi$, for any $x\in U\cap\{\phi>0\}$, $x$ can be path-connected to one of the corresponding faces of the hyperrectangle without leaving the region $U\cap\{\phi>0\}$, and analogously for points in $U\cap\{\phi<0\}$.

We next establish a sufficiently accurate means to determine whether a polynomial attains positive or negative values on a given hyperrectangle. The following definition sets the requirements on an algorithm implementing this task.

Note that a method which exactly computes the minimum and maximum value of $\phi$ on $U$ would trivially yield a sign evaluation algorithm. However, computing the extrema of arbitrary polynomials is a non-trivial and computationally expensive task, perhaps as difficult as the connected component labeling problem itself. The conditions of a sign evaluation algorithm are weaker than an exact computation, thereby allowing for simpler and more efficient algorithms; nevertheless, some degree of certitude is required. In particular, condition (iii) requires an implementation be at least as accurate as conventional Bernstein range evaluation; this ensures it inherits at least the same level of accuracy under the action of subdivision. Condition (iv) says that if $\phi$ happens to be coefficient monotone, the sign evaluation must be at least as accurate as what can be determined solely from the corresponding faces where $\phi$’s extrema are known to occur. One particularly important property is that if a sign evaluation algorithm outputs either $\{-1\}$ or $\{+1\}$, then we may conclude, with absolute certainty, $\phi$ is uniformly signed throughout $U$. An example of a simple sign evaluation algorithm is given in the next section.

Finally, a notion of whether the topology of $U\setminus\{\phi=0\}$ is “simple enough” is recursively defined as follows. This definition implicitly depends on the presence of a suitable sign evaluation algorithm $\sigma_{d}$, so it is assumed that one has been implemented and fixed ahead of time.

A simply connected polynomial guarantees a simple topology of both $\Omega^{-}:=U\cap\{\phi<0\}$ and $\Omega^{+}:=U\cap\{\phi>0\}$, as follows. If conditions (i) or (ii) in Definition 2.3 are met, then $\phi$ is non-zero throughout the hyperrectangle and so exactly one of $\Omega^{\pm}$ is empty and the other the whole hyperrectangle. Otherwise, there is a coordinate axis $k$ such that any point in $x\in\Omega^{\pm}$ can be path-connected to one of the corresponding faces. These faces are themselves simply connected; inductively it follows that if $\Omega^{+}$ (resp., $\Omega^{-}$) is nonempty, then $\Omega^{+}$ (resp., $\Omega^{-}$) has exactly one connected component.³³3In fact, the sets $\Omega^{+}$ or $\Omega^{-}$ are simply connected in the conventional topological sense, but it is unnecessary to prove this fact in the present work.

To prove the correctness of the connected component labeling algorithm, it is useful to establish the following two properties of simply connected polynomials:

• •

  The first property is that if $\phi$ is simply connected on a hyperrectangle $U$, then $\phi$ is simply connected on every face of $U$. This can be shown via induction on the dimension; we illustrate here with a three-dimensional example. Suppose $\phi$ is simply connected on a rectangular 3D prism. If $\sigma_{d}(\phi)\in\bigl{\{}\{-1\},\{+1\}\bigr{\}}$, then it is uniformly signed and thus trivially simply connected on every face of the prism. Otherwise, $\phi$ is coefficient monotone in the up direction, say, such that the restriction of $\phi$ to the bottom and top faces is simply connected; in particular, by induction, we have that $\phi$ is simply connected on every edge of the bottom and top face. Now, on each vertical face of the prism (those faces excluding the bottom and top), the restriction of $\phi$ is coefficient monotone in the up direction, and, as was just concluded, $\phi$ is simply connected on the corresponding lower and upper edges. Therefore, $\phi$ is simply connected on every vertical face. The preceding arguments can be extended to any dimension, and the one-dimensional base case of the inductive argument trivially holds.

• •

  The second property is that if $\phi$ is simply connected on a hyperrectangle $U$, then a sign evaluation algorithm yields exact results, i.e., the output of $\sigma_{d}(\phi,U)$ is precisely $\{\operatorname{sign}(\phi(x)):x\in U\}$. Similar to before, this can be shown via induction. Suppose $\phi$ is simply connected on $U$. If $\sigma_{d}(\phi,U)\in\smash{\bigl{\{}\{-1\},\{+1\}\bigr{\}}}$, then $\phi$ is with absolute certainty uniformly signed and $\sigma_{d}$ is exact. Otherwise, if $\phi$ is coefficient monotone on $U$ in the up direction, say, then $\sigma_{d}(\phi,U)$ is at least as accurate as $\sigma_{d-1}$ applied to the restriction of $\phi$ on the lower and upper faces of $U$; by induction, the latter calculation is exact, and therefore $\sigma_{d}$ is, too.

Combining these two properties, we observe that if $\phi$ is simply connected on a hyperrectangle $U$, then a sign evaluation algorithm yields exact results on every face of $U$.

With these preliminaries established, we are now in a position to describe the connected component labeling algorithm. First, we describe a sufficiently-accurate Bernstein polynomial range evaluation algorithm, which approximately (and in various cases, exactly) evaluates the range of $\phi$ on a given hyperrectangle $U$. Our particular implementation is given in Algorithm 1 and has the following properties:

• •

  The output is an interval such that $[\inf_{x\in U}\phi(x),\sup_{x\in U}\phi(x)]\subseteq\operatorname{range}(\phi,U)$ always holds.

• •

  The algorithm is at least as accurate as conventional Bernstein range evaluation, i.e., $\operatorname{range}(\phi,U)\subseteq[\min_{i}c_{i},\max_{i}c_{i}]$, the latter interval bounding the minimum and maximum value of the input’s Bernstein coefficients.

• •

  For the simple cases of linear, quadratic, or cubic polynomials in $d=1$ dimensions, it is a particularly simple and efficient task to exactly evaluate the range of $\phi$; this is implemented on line 2 and is an easy tweak benefiting overall accuracy.

• •

  If $\phi$ happens to be coefficient monotone on $U$, then the range evaluation algorithm is as least as accurate as what can be determined by evaluating the range on the corresponding faces of $U$ where the extrema of $\phi$ is known to occur (line 8). Further, if $\phi$ happens to be coefficient monotone in multiple directions, then the tightest possible range is returned (lines 9–10).

The main purpose of the range evaluation algorithm is to implement a sign evaluation algorithm. In particular, we define

In other words, $\sigma_{d}(\phi)$ is defined by the possible signs of $\phi$ as computed by the range evaluation algorithm. By construction, it satisfies all the required properties set out in Definition 2.2. In particular, it is exact for low-degree one-dimensional polynomials, whenever $0\notin\operatorname{range}(\phi)$, as well as various other cases. Using the sign evaluation algorithm, it is also straightforward to implement an algorithm that tests whether $\phi$ is simply connected on a given hyperrectangle $U$; our implementation carries out the algorithmic process naturally suggested by Definition 2.3.

With these ingredients, we can now describe the connected component labeling algorithm. Given a Bernstein polynomial $\phi$ and a hyperrectangular constraint domain $\mathcal{U}$, we first test whether $\phi$ is simply connected on $\mathcal{U}$. If it is, then $\Omega^{+}:=\mathcal{U}\cap\{\phi>0\}$ and $\Omega^{-}:=\mathcal{U}\cap\{\phi<0\}$ are either empty or have exactly one component, in which case the labeling problem is trivial. Otherwise, we recursively subdivide $\mathcal{U}$ into smaller and smaller subcells until either $\phi$ is simply connected on the subcell or a maximum recursion depth is met. In this work, we have opted for a quadtree/octree-style subdivision wherein the subcells maintain the aspect ratio of $\mathcal{U}$, though other subdivision algorithms are certainly possible. Upon conclusion of the subdivision process, we have a grid in which most (and in many cases, all) leaf cells have simple topology. In the second phase of the algorithm, two graphs ${\mathcal{G}}^{\pm}$ are constructed, each representing the connectivity of the cells overlapping $\Omega^{+}$ and $\Omega^{-}$. For example, if $\sigma_{d-1}$ applied to a shared face $\mathcal{F}=U_{i}\cap U_{j}$ of the grid yields $s\subseteq\{-1,0,+1\}$ and $-1\in s$ (resp., $+1\in s$), then an edge $(U_{i},U_{j})$ is inserted into ${\mathcal{G}}^{-}$ (resp., ${\mathcal{G}}^{+}$). These edges reflect the fact that there is likely (and in most cases, definitely) a path connecting $U_{i}\cap\Omega^{\pm}$ and $U_{j}\cap\Omega^{\pm}$. Ultimately, it is the connected components of ${\mathcal{G}}^{-}$ and ${\mathcal{G}}^{+}$ which define the possible labels of the overall algorithm; computing this labeling can be done with negligible cost via efficient disjoint-set/merge-find data structures. In essence, these graphs establish the globally topology of $\Omega$, whereas each leaf cell of the tree handles the local topology. Finally, to determine which component an arbitrary point $x\in\Omega$ belongs, we apply standard (and fast) tree traversal algorithms to find a leaf cell containing $x$, and then adopt its corresponding label from ${\mathcal{G}}^{+}$ (if $\phi(x)>0$) or ${\mathcal{G}}^{-}$ (if $\phi(x)<0$).

Algorithm 2 and Algorithm 3 sketch the implementation of the connected component labeling algorithm. A few complementary remarks are in order:

• •

  In Algorithm 2, whenever the restriction of $\phi$ to a subcell $U$ or its face is required, the corresponding Bernstein coefficients can be efficiently and accurately computed via the de Casteljau algorithm.

• •

  On line 7, $\Lambda$ specifies the user-defined maximum recursion depth. For example, if the constraint domain is the unit cube $[0,1]^{d}$, the parameter corresponds to the smallest-possible subcell having width $2^{-\Lambda}$. The effect of $\Lambda$ on the overall labeling is studied in detail in the next section.

• •

  Algorithm 2 only assigns labels to cells that are path-connected to other cells; as such, if a cell contains a part of $\Omega^{+}$, say, that has not been path-connected to any other cell overlapping with $\Omega^{+}$, then that cell does not (yet) have a $\mathchoice{\mathbin{\ooalign{$\displaystyle+$\crcr$\displaystyle\bigcirc$}}}{\mathbin{\ooalign{$\textstyle+$\crcr$\textstyle\bigcirc$}}}{\mathbin{\ooalign{$\scriptstyle+$\crcr$\scriptstyle\bigcirc$}}}{\mathbin{\ooalign{$\scriptscriptstyle+$\crcr$\scriptscriptstyle\bigcirc$}}}$ label defined. Our main motivation for this mainly concerns the scenarios for which the smallest-possible cell size is reached whereon the topology of $\Omega^{+}$ is still uncertain: in these cases, we do not want to unnecessarily create a new label when $\Omega^{+}$ thereon may actually be empty. In a sense, if a sign evaluation algorithm is unable to certify the presence of $\Omega^{+}$, then doing so will be up to the subsequent evaluation of $\chi$; indeed, Algorithm 3 creates a label only when a point $x\in\Omega^{+}$ is given which proves the nonemptiness of $\Omega^{+}$ in that cell.

• •

  By caching the results of the simply connected tests on the leaf cells, as computed by the first phase of Algorithm 2, various steps in the second and last phase can be accelerated. For example, it is unnecessary to compute $s:=\sigma_{d-1}(\phi|_{\mathcal{F}},\mathcal{F})$ on a shared face $\mathcal{F}$ if it is already known that $\phi$ is uniformly signed on one of the corresponding cells. On a related note, if a path connection already exists between $U_{i}$ and $U_{j}$, it is unnecessary to execute lines 18–23 because these lines would not alter that connectivity; thus, some faces can be skipped during the second phase of the algorithm.

• •

  Our particular implementation of the subdivision tree represents each node by a bit-field of width 16 bits. One bit specifies whether the node is a leaf or non-leaf node; the remaining 15 bits of a non-leaf node points to the children of that node; the remaining 15 bits of a leaf node store the type of the leaf cell (2 bits, flagging whether it is uniformly negative, uniformly positive, mixed sign and simply connected, or not simply connected) and 6 bits each to cache the $\mathchoice{\mathbin{\ooalign{$\displaystyle+$\crcr$\displaystyle\bigcirc$}}}{\mathbin{\ooalign{$\textstyle+$\crcr$\textstyle\bigcirc$}}}{\mathbin{\ooalign{$\scriptstyle+$\crcr$\scriptstyle\bigcirc$}}}{\mathbin{\ooalign{$\scriptscriptstyle+$\crcr$\scriptscriptstyle\bigcirc$}}}$ and $\mathchoice{\mathbin{\ooalign{$\displaystyle-$\crcr$\displaystyle\bigcirc$}}}{\mathbin{\ooalign{$\textstyle-$\crcr$\textstyle\bigcirc$}}}{\mathbin{\ooalign{$\scriptstyle-$\crcr$\scriptstyle\bigcirc$}}}{\mathbin{\ooalign{$\scriptscriptstyle-$\crcr$\scriptscriptstyle\bigcirc$}}}$ labels created by Algorithms 2 and 3. This implies an upper bound of $2^{6}-1=63$ components for each phase, ample for our purposes. In addition, our implementation of Algorithm 2 bypasses the explicit creation of the graphs ${\mathcal{G}}^{\pm}$; instead, an efficient disjoint-set/merge-find data structure is used to represent the edge formations and subsequent label creation and lookup.

• •

  Finally, our implementation also makes use of a fuzzy threshold to robustly handle roundoff error in fixed-precision arithmetic. A tolerance $\epsilon$ is incorporated into the sign evaluation algorithm as follows: if $R=\operatorname{range}(\phi,U)$, then $\sigma_{d}(\phi,U)=\{+1\}$ iff $\inf R\geq\epsilon$, $\sigma_{d}(\phi,U)=\{-1\}$ iff $\sup R\leq-\epsilon$, and $\sigma_{d}(\phi)=\{-1,0,+1\}$ in all other cases. This tolerance helps to prevent inconsistencies arising from roundoff error, such as what might occur when the zero level set $\{\phi=0\}$ grazes the boundary of a cell. Here, the tolerance is chosen to scale relative to the Bernstein coefficients of $\phi$ on $\mathcal{U}$; in particular, we have set $\epsilon=10^{3}\,\epsilon_{0}\max_{i}|c_{i}|$, where $\epsilon_{0}$ is machine epsilon (approximately $10^{-16}$ in double-precision arithmetic).

To conclude this section, we prove the correctness of the connected component labeling algorithm in the case the subdivision process finishes without hitting the maximum-imposed recursion depth. In other words, every leaf cell of the tree satisfies the simply connected property. In the prior section we saw that the sign evaluation algorithm is exact in this case, and so two adjacent cells $U_{i}$ and $U_{j}$ are $\mathchoice{\mathbin{\ooalign{$\displaystyle+$\crcr$\displaystyle\bigcirc$}}}{\mathbin{\ooalign{$\textstyle+$\crcr$\textstyle\bigcirc$}}}{\mathbin{\ooalign{$\scriptstyle+$\crcr$\scriptstyle\bigcirc$}}}{\mathbin{\ooalign{$\scriptscriptstyle+$\crcr$\scriptscriptstyle\bigcirc$}}}$-connected if and only if $(U_{i}\cap U_{j})\cap\{\phi>0\}$ is nonempty. It follows that if $x,y\in\Omega^{+}:={\mathcal{U}}\cap\{\phi>0\}$ are path connected, then a corresponding path would go through zero or more faces overlapping with $\Omega^{+}$, all of which must have been assigned the same label by Algorithm 2. Conversely, if for two points $x,y\in\Omega^{+}$ the evaluation of Algorithm 3 yields $\chi(x)=\chi(y)$, then there must be a sequence of leaf cells $U_{1},U_{2},\ldots,U_{m}$ with $x\in U_{1}$ and $y\in U_{m}$ such that, for each pair, $(U_{i}\cap U_{i+1})\cap\{\phi>0\}$ is nonempty. Therefore, there is a sequence of points $z_{i}\in(U_{i}\cap U_{i+1})\cap\{\phi>0\}$, and, because each $U_{i}\cap\Omega^{+}$ has exactly one component, there exists a path connecting $x\rightarrow z_{1}\rightarrow z_{2}\rightarrow\cdots\rightarrow z_{m-1}\rightarrow y$, and so $x$ and $y$ are path connected. The same arguments apply analogously to the negative region, and it is trivial to see that $\chi(x)\neq\chi(y)$ whenever $\phi(x)\phi(y)<0$, since the negative and positive regions do not have overlapping label identifiers. In summary, when all the leaf cells satisfy the simply connected property, the output of Algorithm 3 is exact in the sense that $\chi(x)=\chi(y)$ if and only if $x$ and $y$ are path connected. In fact, we will observe in the next section that in many instances, the labeling algorithm continues to be exact even when the maximum recursion depth is reached.

In this section, we assess the effectiveness of the connected component labeling algorithm on randomly generated geometry in 2D and 3D as well as on particular examples exhibiting cusps, self-intersections, and other kinds of singularities.

The connected component labeling algorithm developed here operates by recursively subdividing the constraint hyperrectangular domain $\mathcal{U}$ into progressively smaller subcells until the topology of $\Omega:=\mathcal{U}\setminus\{\phi=0\}$ thereon is sufficiently simple. This is achieved through a “simply connected” test exploiting some useful properties of the Bernstein polynomial basis, namely its effective methods for polynomial range evaluation. If the subdivision process succeeds, then the labeling is certifiably exact, i.e., $\chi(x)=\chi(y)$ if and only if $x,y\in\Omega$ are path-connected. At its core, however, the connected component labeling problem is ill-posed in the sense that tiny changes in the coefficients of the input polynomial $\phi$ could alter the number and topology of the connected components. We treated that here by allowing distinct components to be glued together, but only when the topology is uncertain: gluing occurs mainly when two components are nearly touching relative to the length scale of the smallest-permitted subcell; gluing also occurs for various edge cases such as interfacial self-intersections or junctions, or grazing-type singularities such as non-negative polynomials for which $\{\phi>0\}$ has multiple components.

We have yet to discuss the speed of the algorithm. In general, the computational complexity is difficult to precisely quantify, being intricately dependent on the input polynomial $\phi$. For a fixed spatial dimension $d$ and bounded polynomial degree, the construction of the subdivision tree (Algorithm 2) has a worst-case complexity $o(2^{d\Lambda})$; in practice, however, it is usually significantly faster. Computation of $\chi(x)$ for some $x\in\Omega$ (Algorithm 3) has a worst-case complexity ${\mathcal{O}}(\Lambda)$; in practice, however, it is usually significantly faster. Indeed, for a given fixed polynomial, as soon as the subdivision process resolves its implicitly-defined geometry, the algorithmic complexity becomes independent of $\Lambda$.⁸⁸8In pathological cases, such as those illustrated in Fig.​ 7, the computational complexity depends on the measure of the input polynomial’s set of singularities; in particular, if its set of singularities (inside $\mathcal{U}$) is nonempty and has maximal manifold dimension $d_{s}<d$, the worst-case complexity of Algorithm 2 is ${\mathcal{O}}(\lambda)$ when $d_{s}=0$, or ${\mathcal{O}}(2^{d_{s}\Lambda})$ when $d_{s}>0$, as $\Lambda\to\infty$. A comparison to the high-order quadrature algorithms developed in [18] is perhaps useful: solving the connected component labeling problem is about as fast as building a moderate order quadrature scheme. These quadrature algorithms are already “fast” in some sense, so the speed of the labeling algorithm is sufficient for our intended usage. As an additional indication, the cost of tree construction for the randomly generated geometry problems in §​ 4.1 is about $8$ microseconds per instance in 2D and $0.2$ milliseconds per instance in 3D; about two thirds of that time is spent in the first phase of Algorithm 2, about one third in the second phase, and negligible time in third phase. After tree construction and on the same set of test problems, the cost of evaluating $\chi(x)$ for a random point $x$ is about $0.05$ (resp., $0.1$) microseconds per evaluation in 2D (resp., 3D).⁹⁹9Timing measurements obtained on an Intel Xeon E3-1535m v6 laptop, single core, operating at approximately 3.5 Ghz.

Finally, we mention here some possibilities for extending the algorithmic approach. One immediate possibility is to develop a more sophisticated subcell topology test, e.g., $\Omega^{\pm}\cap U$ being star-connected or via an analysis of convexity [9]. Accompanying this change, a more sophisticated sign evaluation algorithm might also be required, so as to appropriately path-connect neighboring cells, consistent with the updated topology test. Depending on the end application, however, these alterations might come with an increased computational cost. Another possibility is to replace the quadtree/octree-style subdivision with a k-d tree subdivision, or, more generally, a binary space partitioning (BSP) tree. In this setting, the most efficient subdivision process would likely occur by orienting to the features of the polynomial, e.g., by choosing hyperplanes which are locally parallel to its level sets. One would then have to generalize the range evaluation algorithms to polytopes and contend with more complex algorithms for building the cell connectivity. On the other hand, the quadtree/octree approach allows for a considerably simpler implementation, partly due to the use of tensor-product polynomials as well as dimension reduction and recursion techniques. As a final possibility, we note that the basic idea of the subdivision algorithm could be extended to non-polynomial level set functions; for example, one approach might be to apply automatic differentiation and interval arithmetic techniques (see, e.g., [17]) to devise a suitable subcell topology test for arbitrary level set functions.

This research was supported in part by the Applied Mathematics Program of the U.S. Department of Energy (DOE) Office of Advanced Scientific Computing Research under contract number DE-AC02-05CH11231, and by a DOE Office of Science Early Career Research Program award.