Adaptive grid generation for discretizing implicit complexes

We present a method for generating a simplicial (e

We present a method for generating a simplicial (e.g., triangular or tetrahedral) grid to enable adaptive discretization of implicit shapes defined by a vector function.Such shapes, which we call implicit complexes, are generalizations of implicit surfaces and useful for representing non-smooth and non-manifold structures.While adaptive grid generation has been extensively studied for polygonizing implicit surfaces, few methods are designed for implicit complexes.Our method can generate adaptive grids for several implicit complexes, including arrangements of implicit surfaces, CSG shapes, material interfaces, and curve networks.Importantly, our method adapts

INTRODUCTION
Implicit representations are widely used in computer graphics.A smooth and manifold surface in 3D can be represented as the level set of a scalar function  : R 3 → R, also known as an implicit surface.The implicit surface enjoys a number of benefits, including a simple definition, easy modification (in both geometry and topology), and convenience for operations such as offsets and boolean.Shapes beyond smooth manifolds can also be represented implicitly, by making use of a vector function  : R 3 → R  where  > 1.For example, a CSG shape [Requicha and Voelcker 1977], often used to represent piecewise smooth surfaces, consists of subsets of multiple implicit surfaces, each can be defined by the level set of a scalar component of  .A non-manifold network of surfaces can be represented either as the arrangement of multiple implicit surfaces [Bagley et al. 2016;Guo et al. 2021], or as the interface between regions where one component of  is greater than the rest [Bertram et al. 2005;Zhang et al. 2008].A vector function can also define lower-dimensional shapes, such as spatial curves defined by the intersection of two implicit surfaces [Burns et al. 2005;Edelsbrunner and Harer 2002;Kohlbrenner et al. 2023].We refer to a shape implicitly defined by some vector function as an implicit complex.Some examples can be found in Figure 1 (bottom).
To be useful for downstream applications, an implicit surface or complex must be discretized.Many discretization methods, such as Marching Cubes [Lorensen and Cline 1987] and Marching Tetrahedra [Bloomenthal 1988], operate on a spatial grid consisting of cubic or simplicial cells (reviewed in Section 2.1).The grid structure plays a key role in the performance of these methods.While a finer grid generally leads to higher discretization accuracy, it comes at the cost of slower grid generation and discretization, a larger memory footprint, and a larger output size.An ideal grid for discretization should be adaptive in that the finer grid cells are located only where higher accuracy is needed.These locations are typically where the implicit shape has a non-trivial geometry or topology.
While adaptive grid generation has been extensively studied for polygonizing implicit surfaces, works on implicit complexes have been scarce (see Section 2.2).A key challenge in the latter is adapting the grid structure to accurately discretize the intersection of multiple implicit surfaces.Such intersections form the lower-dimensional elements in an implicit complex, such as the sharp edges and corners in a CSG shape or the junction graph in a non-manifold surface network.Note that the intersection of implicit surfaces may have a more complex geometry than the surfaces themselves, as shown in the example of Figure 2 (a).As a result, grids adapted only to the implicit surfaces could suffer from either insufficient resolution at their intersections (Figure 2 (b)) or excessive refinement elsewhere (Figure 2 (c)).Existing methods that can capture such intersections either are limited to low-degree polynomials [Dupont et al. 2007;Mourrain et al. 2005;Schömer and Wolpert 2006] or rely on interval analysis [Allgower and Georg 1989;Boissonnat and Wintraecken 2022], and hence they are not suited for general functions where tight intervals can be difficult, if not impossible, to obtain.
This paper presents a new method for generating adaptive simplicial grids for discretizing a variety of implicit complexes (such as those in Figure 1).A key feature of the method is that it adapts the grid structure to the geometry of both the implicit surfaces and their intersections.This allows elements at all dimensions in an implicit complex to be accurately discretized without excessive refinement of the grid (see Figure 2 (d)).Unlike existing intersectionaware methods, our method does not require intervals, and it can be applied to any function whose values and gradients can be queried.
Our main contribution is a set of criteria that decide whether a simplicial cell needs to be refined, given a particular definition of the implicit complex.The core of these criteria are novel tests that (approximately) check if an -simplex contains the intersection of  implicit surfaces, and if such intersection can be well approximated by linear interpolation.Our tests can be efficiently and robustly implemented in any dimension  and for any  ≤ , without the need for interval analysis.These tests become the building blocks of the refinement criteria for an implicit complex, which include additional checks to avoid unnecessary refinement along parts of an implicit surface (or intersection) that does not lie on the complex (e.g., the trimmed portion of a primitive in CSG).
Our second contribution is a new method for top-down refinement of a simplicial grid, given some refinement criteria.Our method is a variant of the classic longest edge bisection (LEB) method [Plaza and Rivara 2003;Rivara 1991], which bisects a simplex at the midpoint of its longest edge.While LEB generates high quality cells with smoothly varying sizes, the grid can be overly fine for the purpose of discretizing a single implicit shape (see Figure 3 (a)).Our variant produces fewer cells than LEB by allowing cells away from the implicit shape, which are irrelevant to discretization, to have worse shapes and hence higher adaptivity (see Figure 3 (b)).
The rest of the paper is organized as follows.After discussing the previous works in Section 2, we review the implicit representations in Section 3. The technical discussion starts with our tests on the implicit surface intersections in Section 4, followed by our criteria for implicit complexes in Section 5, and it ends with our refinement algorithm in Section 6.We present our experimental results in Section 7 and conclude with a discussion in Section 8. Our method produces much fewer triangles away from the curve, while maintaining the shape and adaptivity of triangles around the curve (gray).

PREVIOUS WORK
There is an extensive literature on discretizing implicit shapes.We refer the reader to the excellent and comprehensive survey of De Araújo et al. [2015].Here we selectively review works that are most relevant to ours, namely methods that discretize the shape using an embedded grid and methods to generate such an embedding.Note that these methods are different from works that produce grids that conform to (rather than surround) a given surface [Dey and Slatton 2013;Oudot et al. 2010].

Grid-based discretization
A common approach to discretize an implicit surface is to divide the domain into cells and approximate the surface locally within each cell using polygons [Whitney 1957].Marching Cubes [Lorensen and Cline 1987] and its many variants adopt uniform cubic cells.The lack of adaptivity of a cubic grid motivates methods that can work on an adaptive octree [Ju et al. 2002;Schaefer and Warren 2004;Varadhan et al. 2004].Alternatively, adaptivity can be achieved by using simplicial cells (e.g., tetrahedra in 3D), which can tile space with varying sizes.Another advantage of simplicial cells is that the values at the simplex's vertices can be interpolated by a linear function.The linearity enables simple and unambiguous discretization, which is exploited by the Marching Tetrahedra method [ Bloomenthal 1988].

Grid generation for discretization
Grid generation is important for fields such as finite element analysis (FEA) [Brandts et al. 2020;George et al. 2017] and data visualization [Borgo et al. 2004].As the criteria of a "good" grid differ by how the grid is used, we focus on methods designed for discretizing implicit shapes.Such methods (ours included) typically take a topdown refinement approach, where an initial coarse grid is iteratively refined where needed.The key questions that need to be answered are (1) which cells need to be refined, and (2) how to refine them.We review existing solutions for each problem.
2.2.1 Refinement criteria.For the purpose of discretization, a cell is generally regarded as refinable if it contains some part of the implicit shape (often known as the zero-crossing test), and if that part is insufficiently approximated by the discretization.
Many refinement criteria have been developed for implicit surfaces.Criteria based on interval analysis allow the refinement algorithm to offer strong topological guarantees [Boissonnat et al. 2008;Chattopadhyay et al. 2012;de Figueiredo et al. 2006;Plantinga and Vegter 2006;Stander and Hart 1997] and/or accurately capture the geometric details [de Figueiredo et al. 2006].However, tight intervals are difficult to obtain for many practical functions, and impossible when the function is provided as a "black-box" (e.g., a deep neural network), while loose intervals may lead to excessive grid refinement [Chattopadhyay et al. 2012].Criteria that do not rely on intervals resort to approximate means to measure the curvature of the function [Azernikov and Fischer 2005;Bloomenthal 1988;Hui and Jiang 1999;Molino et al. 2003;Schmidt 1993] or the deviation of the discretization from the implicit surface [Hall and Warren 1990;Liang and Zhang 2014;Petersen et al. 1987;Zhang et al. 2005].
In contrast, refinement criteria for implicit complexes are few and far between.Criteria have been developed to exactly discretize the IA of low-degree polynomials [Dupont et al. 2007;Mourrain et al. 2005;Schömer and Wolpert 2006].Interval-based criteria can offer topological and/or geometric guarantees in discretizing the intersection of multiple level sets [Allgower and Georg 1989;Boissonnat and Wintraecken 2022] or singularities in the projection of surfaces from R 4 [Diatta et al. 2019].However, these criteria are not suited for general functions or functions without tight intervals.A few other criteria have been proposed for IA [Kim et al. 2000], CSG [de Miras and Feito 2002;Tobler et al. 1995], and MI [Zhang and Qian 2012], but all of them are concerned with only the surface approximation quality.To the best of our knowledge, there are no interval-free refinement criteria for any implicit complex that consider the approximation quality of the low-dimensional elements in the complex (i.e., intersection curves and points).

Refinement method.
For discretization purposes, we consider refinement methods that produce a conformal simplicial grid, where adjacent cells share complete faces (i.e., no "hanging" nodes).Another desirable feature of a simplicial grid is well-shaped cells.In the context of discretization, poorly shaped (i.e., near-degenerate) cells are more likely to yield poorly shaped triangles in the discretization.Although such elements can be removed in a post-process, for example by relocating grid points [Hall and Warren 1990;Raman and Wenger 2008]) or remeshing [de Figueiredo et al. 2006;Dillard et al. 2007], excessive amount of near-degenerate elements can be challenging to remove.
One way to produce well-shaped simplices adaptively is regularly subdividing each simplex that needs to be refined, followed by a local refinement of adjacent simplices to maintain a conformal grid [Bey 1995;Grosso et al. 1997;Hall and Warren 1990;Hui and Jiang 1999;Kim et al. 2000;Liu and Joe 1995;Zhang 1995;Zhao et al. 2021].Note that subdividing a simplex regularly creates many new cells.A more economic approach is bisection, which splits a simplex into two halves at the midpoint of a chosen edge.A common choice of the bisecting edge is the longest edge [Plaza and Rivara 2003;Rivara 1991], although other choices have also been investigated [Arnold et al. 2000;Bänsch 1991;Belda Ferrín et al. 2022;Kossaczký 1994;Mitchell 2016].Like regular subdivision, bisecting a simplex is followed by local refinement of surrounding simplices to maintain conformity.When applied to the Kuhn-subdivision of a square or cube, the longest edge bisection (LEB) method creates simplices that belong to 2 or 3 similarity classes [Maubach 1995].Efficient data structures have also been developed [Gerstner 2003;Weiss andDe Floriani 2009, 2011].
The excellent cell quality, combined with improved compactness over regular simplex subdivision, make LEB ideal for applications such as FEA [Anderson et al. 2021] and data visualization [Gerstner and Pajarola 2000;Gregorski et al. 2002;Pascucci 2004;Zhou et al. 1997].However, the small number of similarity classes of cells limits the adaptivity of LEB for discretizing implicit shapes.This results in unnecessarily high refinement away from the discretization (see Figure 3 (a)), which we try to avoid in our new method.

PRELIMINARIES: IMPLICIT REPRESENTATIONS
Before introducing our method, we briefly review several implicit shape definitions that are considered in this paper.

Level sets and zero sets
The level set of a scalar function  : R  → R at level  ∈ R, which we denote by  −1 (), is the loci of points where  evaluates to : For a smooth function  , the level set is generally a smooth ( − 1)dimensional manifold.The level set is better known in computer graphics as the implicit curve or implicit surface for  = 2, 3. We call the level set  −1 (0) the zero set, which we abbreviate as  −1 .The definition in Equation 1 naturally generalizes to the level set of a vector-valued function  : R  → R  , where  > 1, as the loci of points where  evaluates to a vector level .We shall refer to  −1 () the vector level set when  > 1 and scalar level set when  = 1.Note that  can be considered as a set of  scalar functions { 1 , . . .,   }, which we call the components of  .Geometrically,  −1 () is the intersection of the  scalar level sets of its components, where  = { 1 , . . .,   }.As a result,  −1 () is a ( − )-dimensional manifold, and it generally only exists for  ≤ .The vector level set is also known as the implicit manifold [Allgower and Georg 1989;Kohlbrenner et al. 2023] or iso-manifold [Boissonnat and Wintraecken 2022].For  = 3, the vector level set is the intersection curve ( = 2) or point ( = 3) of  implicit surfaces.We continue to use  −1 to denote the vector zero set,  −1 (0), which is the intersection of  scalar zero sets  −1  .

Implicit complexes
We use the term implicit complex to refer generally to any shape composed of zero sets of one or multiple (scalar or vector) functions.
Our method is designed for the implicit complexes described below, which have found many applications in computer graphics.Each of these implicit complexes is defined by a single vector function  : R  → R  for some  > 1.Specifically, the complex consists of zero sets of either subfunctions of  or their difference bases.We call a function  ′ a subfunction of  , if the components of  ′ are a subset of components of  .Given a function  ′ with  components, the difference basis,  ( ′ ), is a linearly independent set of  − 1 pairwise differences between the components of  ′ .
(1) The implicit arrangement (IA) of  is the loci of points where at least one of its component evaluates to zero: Geometrically, the IA is the union of zero sets of all subfunctions of  with no more than  components.In R 3 , the IA consists of multiple implicit surfaces together with their intersection curves and points.An 2D illustration is given in Figure 4 (right).The IA has been used for representing complex objects with interior partitioning, such as geological structures [Bagley et al. 2016;Guo et al. 2021].It also includes other implicit complexes as its subsets, such as the ones described below.
(2) Constructive solid geometry (CSG) [Requicha and Voelcker 1977], a popular modeling approach, constructs a solid shape from a collection of primitive shapes (e.g., planes, spheres, cylinders, etc.) using set operations (union, intersection, and complement).The boundary of an CSG shape can be implicitly defined by a vector function  , whose components encode the primitives as scalar zero sets, and an CSG expression  consisting of max, min and negation operators that encode the set operations.The CSG boundary is the scalar zero set of the composition  •  : Since the composition  •  is typically only  0 continuous, the CSG boundary is a piecewise-smooth ( − 1)-dimensional manifold.Each smooth ( − 1)-dimensional piece is a subset of the scalar zero set of a component of  , and adjacent smooth pieces meet at vector zero sets of subfunctions of  .By definition,  ( ) ⊆ ( ) (see illustration in Figure 4 middle).
(3) The material interface (MI) of  partitions a domain into regions where one component of  dominates; that is, some   is strictly greater than all other   ( ≠ ) in that region.Equivalently, MI is the loci where multiple components of  attain the maximum, The MI is a network of ( − 1)-dimensional manifolds meeting at non-manifold junctions.Each ( − 1)-dimensional piece separates regions dominated by two components of  (say   ,   ), and it lies on the scalar zero set of their difference (  −   ).The junction where more than two regions meet lies on the vector zero set of the difference basis of the subfunction of  whose components dominate those regions (see illustration in Figure 4 left).The MI is a subset of the IA of all pairwise differences of the components of  .MI has been widely used for representing complex objects with interior partitions, such as multi-phase fluids [Kim 2010;Losasso et al. 2006], composite materials [Dillard et al. 2007;Shammaa et al. 2010], and anatomical structures [Bertram et al. 2005;Zhang et al. 2008].
(4) Finally, the lower-dimensional elements of an implicit complex form an implicit complex themselves, which consist exclusively of vector zero sets.This includes the intersections of two or more scalar zero sets in IA, the locations where two or more smooth pieces meet in CSG, and the junction where more than two regions meet in MI.In R 3 , the intersection curves of implicit surfaces have been used for visualizing line features in volume data [Burns et al. 2005;Edelsbrunner and Harer 2002] and reconstructing lower-dimensional shapes [Kohlbrenner et al. 2023].

REFINEMENT CRITERIA FOR ZERO SETS
We start by presenting our refinement criteria for zero sets of scalar and, more importantly, vector functions, which are the building blocks of implicit complexes.Our criteria answer the following general question: given a function  : R  → R  for any dimension  and  ∈ [1, ] and an -simplex , does  need to be refined to better discretize the zero set  −1 ?In R 3 , our criteria consider a single implicit surface ( = 1), the intersection curve of two surfaces ( = 2), or the intersection point of three surfaces ( = 3).
Our criteria follow the same principles as existing ones for implicit surfaces.That is, refinement is not needed if either  does not contain any part of  −1 , or if  −1 is already close enough to the discretization.Unlike existing criteria, our criteria are formulated so that they are generalizable to vector zero sets and also can be efficiently computed.Specifically, let  be the linear approximation of  inside  by barycentric interpolation of values of  at the vertices of .We consider the zero set of  ,  −1 , as the discretization of  −1 in .We deem  refinable for  −1 if it passes two tests: • Zero-crossing test: • Distance test: where   (,  ) = sup  ∈  (,  ) is the one-sided Hausdorff distance from  to  , and  is a user-defined threshold.
The distance test checks if the discretization is too far from the zero set.Our specific choices in the test, including the use of the one-sided Hausdorff distance and the complete zero set  −1 (instead of its portion in ), are made to enable efficient computations, as we shall see below.
While these tests can be performed using interval analysis, intervals on the values and derivatives are not always available for arbitrary functions.Instead, we perform the tests on local approximations of  , called proxies, that can be constructed by sampling the functions at spatial locations (Section 4.1).Our design of the proxy leads to efficient implementation of both tests (Sections 4.2 and 4.3); see Equations 9 and 15.

Proxy construction
We consider the class of functions  represented as a convex combination with linear precision as our proxy.Such a function is defined by control points  1 , . . .,   that lie inside or on the boundary of the simplex , where each   is associated with a control value (a length- vector) where the weights   () satisfy, for all  ∈ : (1) Convexity:   () ≥ 0 for  = 1, . . ., , and  =1   () = 1.(2) Linear precision:  =1   ()  = .
The convexity of the weights ensures that  () lies in the convex hull of  = { 1 , . . .,   }, where each   is treated as a point in R  .With linear precision,  can reproduce the linear interpolation function  by giving   the linearly interpolated value  (  ), which we denote as   .These properties will be exploited later in implementing the two tests.An example of such  is the Bezier simplex [Farin 1986], where the control points lie on a regular grid in  and the weights   are the Berstein polynomials.One could also place the control points in arbitrary locations on the boundary of  and use one of the generalized barycentric coordinates [Floater 2015] as the weights (as long as the coordinates are positive).
The refinement criteria to be introduced in Sections 4.2 and 4.3 apply to any proxy satisfying the properties above.In our implementation, we adopt the cubic Bezier simplex as the proxy.The control points   in a cubic Bezier simplex consist of all vertices of , the two trisectors on each edge of , and the centroid of each triangle face of  (see Figure 5).Although such a proxy can only reproduce cubic polynomials, we have found it to be effective in approximating a variety of non-cubic, and even non-polynomial functions (see Section 7).We obtain the control values   , known as the Bezier ordinates, by evaluating  and its gradient at the simplex vertices (see details in Appendix A).

Zero-crossing test
We now examine the refinement criteria on a proxy function  .The zero-crossing test checks if  −1 ∩ ≠ ∅.As discussed in the previous section, the value  () at any point  ∈  lies in the -dimensional convex hull of , the control values of  .As a result, if there exists some point  such that  () = 0, then 0 lies in this convex hull.So a necessary condition to pass the zero-crossing test is where  is the origin of R  and  () is the convex hull of .
The test in Equation 9 is a generalization of the zero-crossing test for implicit surfaces (i.e.,  = 1) that has been commonly used in previous works [Hall and Warren 1990;Petersen 1984], which checks if the range of control values  encloses 0. For  > 1, the test amounts to checking if  is not an extreme point in the union set  ∪ .The extremity check can be written as a linear programming (LP) problem with the same number of constraints as the number of points ( in our case) [Ottmann et al. 1995], which can be solved in time linear to  if the space dimension ( in our case) is fixed [Megiddo 1984].For lower dimensions (e.g.,  = 2, 3), we use a faster, brute-force implementation that will be described in Section 7.4.

Distance test
We next turn to the distance test on the proxy  , which checks if where  is the linear approximation of  in .While the exact (onesided) Hausdorff distance between the two zero sets can be difficult to compute, we observe that it has the following upper bound: (The second equality holds because  () = 0 for all  ∈  −1 ∩ ).
In other words, the upper bound is the maximal distance from any point  ∈  to the level set of the linear function  at the level  ().In the following, we will show that this maximum distance is attained at a control point of  , and hence the test in Equation 10can be expressed solely by the control values .
We start by deriving an explicit form for the upper bound.It helps to first consider the simple case of  = 1.Since  is a scalar function, its level sets are ( − 1)-dimensional hyperplane in R  orthogonal to the gradient  = ∇ .The distance from  to the level set of  at level  () (i.e., Note that this distance is signed, and a positive distance implies that  is on the opposite side of the level set as .
Furthermore,  must lie in the subspace spanned by the basis  = { 1 , . . .,   }.We can therefore find  as where  = (  ) −1 .The distance  (,  −1 (  ())) can be written as To find the maximal distance, observe that both  (),  () are convex combinations with the same weights.This is because  can be reproduced by  when the control values  are the barycentric interpolants , as discussed in Section 4.1.As a result, the difference  () −  () is a convex combination of values  − , and hence it lies inside the convex hull of  −  (treated as points in R  ).Since linear transformations (such as ) preserve convexity,  lies in the convex hull of the transformed points  ( − ).Finally, utilizing the fact that the maximum distance from the origin to a convex hull is realized at its vertices, we arrive at the bound: are scalar values and  is the gradient vector of  .This test is similar to tests proposed by other researchers [Gregorski et al. 2002;Petersen et al. 1987;Zhang et al. 2005], who also consider bounding the difference between a function and its linear approximation after a gradient-based normalization.Unlike previous works, our test can be applied generally to any vector zero set ( > 1).
One issue that arises in implementing the test in Equation 15 is that it involves inverting matrices that could become singular.For example, computing the gradient  needs to invert a matrix composed of the edge vectors of simplex , which is singular if  is degenerate (i.e., has no -dimensional measure).Computing  requires inverting   , which is singular if rank() <  (e.g., a component of  has zero gradient, or two components have colinear gradients, etc.).To address the issue, we can rewrite the test in an inversion-free form (see Appendix B).The re-formulation not only is robust against singular matrices, but also exhibits improved numerical accuracy when the matrices are almost singular.

REFINEMENT CRITERIA FOR IMPLICIT COMPLEXES
Building on the refinement criteria for zero sets, we now describe criteria for the implicit complexes discussed in Section 3. Recall that such a complex  is defined by a vector function  : R  → R  where  can be any positive integer (possibly greater than ), and it consists of zero sets of either the subfunctions of  or their difference bases.We start by presenting a general structure of the criteria that is shared by all implicit complexes being considered (Section 5.1), and then discuss details specific to each complex (Section 5.2).

General criteria
Similar to the refinability for zero sets, a simplex  is considered refinable for  if it contains some portion of  and if this portion is too far from its discretization (i.e., the implicit complex defined by the linear approximation  ).
A straight-forward way to check refinability is to see if  is refinable for any zero set that makes up .However, this would lead to over-refinement if  contains partial zero sets, as in the cases of CSG and MI.For example, a 3D CSG shape consists of patches of implicit surfaces that are trimmed by other implicit surfaces.To avoid refining around the part of a zero set that does not belong to , we identify component functions of  that contribute to the portion of  in  as active components.The definition of an active component differs by the type of the implicit complex, which will be introduced in the next section.We deem  refinable if it is refinable for any zero set making up  that involves only active components of  in .
To efficiently check refinability, we observe that the zero set of a function ℎ (or its difference basis,  (ℎ)) intersects  only if the zero set of every subfunction of ℎ also intersects .We therefore examine the subfunctions of ℎ in the order of increasing number of components, and we ignore a function ℎ with  components if  fails the zero-crossing test for any subfunction (or its difference basis) of ℎ with  − 1 components.The pseudo-code of our implementation is given in Algorithm 1.Here, gather(  ) produces a list of ( + 1)-component subfunctions ℎ of  such that every -component subfunction of ℎ is in   .
Our criteria accommodate a separate distance threshold   for checking zero sets of a -component function, for each  ∈ [1, ].This makes it possible to refine the grid only for discretizing the lower-dimensional elements of an implicit complex (e.g., the sharp features in CSG or the junction graph in MI).Specifically, setting   = ∞ for all  <  −  will prevent the grid to be refined around zero sets whose dimension is greater than .

Active components
We now detail, for each type of implicit complex , the criteria for a component   of  to be active in .To avoid missing geometric details, our criteria are intentionally conservative, in that we want to identify all   that possibly contribute to the portion of  in .Our criteria will again utilize the proxy function  , as constructed in Section 4.1, and particularly its control values .Note that, since  is a convex combination of , the value of a component   () at any point  ∈  is bounded within the range of the control values   , which we denote by   .Our criteria for active components will be expressed by this range.
• IA: Since ( ) (Equation 3) contains only complete zero sets, we deem a component   active if its scalar zero set intersects .Using the zero-crossing test (Equation 15), we only need to check if 0 ∈   .• CSG: Recall that  ( ) (Equation 4) is the scalar zero set of the composition  •  , where  is some CSG expression consisting of min, max, and negation operations on the components of  .We consider a component active if its scalar zero set might lie on ( •  ) −1 ∩ , if it is not empty.Leveraging the hierarchical structure of , we collect the list of active components in a recursive manner.We first estimate the range of  •  , denoted by  (), recursively on  following [Duff 1992;Hijazi et al. 2010], where △ can be either min or max, ¬ is negation, and   ,  ℎ denotes the lower and higher ends of a range  .We then recursively collect the list of active components for an expression , denoted by (), as the union of the active lists for the operands of , if  () contains zero, and empty otherwise: if  =   and 0 ∈   ∅, otherwise • MI: Recall that  ( ) (Equation 5) partitions the space into regions within which one component of  dominates the others.We therefore consider a component   active if it dominates all other components at some point in .A necessary condition is that there is no range   ( ≠ ) that is strictly higher than   .Here, we say a range [ 1 ,  1 ] is strictly higher (or lower) than another range It is easy to verify that the ranges in the set  = { 1 , . . .,   } satisfying this condition are those whose higher end is no smaller than the highest lower end among all ranges in .
These ranges, which belong to active components, can be efficiently found by two sweeps over , first identifying the highest lower end and then checking it against each range.

SIMPLICIAL REFINEMENT
We now consider how to refine a simplicial grid using the criteria described in the previous sections.We adopt the top-down refinement regime that is common in grid generation.Given an initial grid and some criteria (e.g., Algorithm 1), our goal is to produce a maximally refined grid where no cell is deemed refinable by the criteria.While our implementation is specialized for the practically useful cases of  = 2, 3, our methodology is not specific to dimensionality.Our method is a variation of the classical longest edge bisection (LEB) method.Recall from Section 2.2.2 that this scheme refines a simplicial cell by splitting it into two at the mid-point of its longest edge, followed by a local refinement of surrounding cells to restore a conformal grid.To avoid the unnecessary refinement in LEB away from the implicit shape, as seen in Figure 3 (a), our key idea is to relax the requirement that every cell has to be split by its longest edge.Instead, for any chosen edge , we bisect all cells incident to  at the mid-point of , regardless of whether  is the longest edge of that cell (see illustration in Figure 7).Note that the operation maintains the conformity of the grid, which avoids the local refinement in LEB.Our choice of the bisecting edge  is motivated by the need to maintain cell quality.Specifically, we call an edge refinable if any of its incident cells is refinable according to the given refinement criteria.At each iteration of the algorithm, we bisect the longest refinable edge in the current grid.We call our method the longest refinable edge bisection (LREB) method.
We have observed that LREB produces more adaptive and compact grids than LEB, as seen in Figure 3 (b).The adaptivity is made possible by larger and more irregularly-shaped cells away from the implicit geometry.Although lacking a theoretical bound, we provide experimental evidence in Section 7.3 that the worst quality among cells enclosing the implicit shape produced by LREB appears to be lower-bounded for a variety of implicit surfaces and complexes.
We give the pseudo-code of our complete refinement algorithm in Algorithm 2. We use a priority queue to maintain all refinable cells, prioritized by the length of each cell's longest edge.Each iteration splits the cell at the head of the queue by its longest edge (which is the longest refinable edge in the grid), as well as all cells (refinable or not) sharing that edge.New cells are checked for refinability and the refinable ones are placed into the queue.The algorithm stops when the queue is empty or some termination criteria is met (e.g., the longest refinable edge is shorter than a given threshold ).

RESULTS
We show results of our refinement method on implicit complexes in both 2D and 3D.In all examples, the domain is a unit square (in 2D) or cube (in 3D).To create the initial grid (  in Algorithm 2), we adopt the Kuhn-subdivision, which splits a square into two triangles and a cube into six tetrahedra, although our algorithm is not specific to this initialization.Unless stated otherwise, the minimum edge length threshold ( in Algorithm 2) is set to 0, and hence the algorithm terminates only when the grid is fully refined (i.e., no cell is deemed refinable).We adopt the discretization method of [Du et al. 2022], which computes the exact IA and MI of the piecewise linear functions in 2D and 3D.CSG is discretized by first computing the IA and extracting the relevant subset.

Zero-crossing and distance tests
We start by evaluating the effectiveness of our two tests, zerocrossing test and distance test (Section 4), on zero sets of functions with different number of components.
We first examine the tests on a scalar zero set.Although they are similar to previous tests in the literature, we include them not only for completeness, but also since their behavior echoes those we shall see on vector zero sets.Figure 8 shows the refined grids for the same implicit curve used in Figure 3, by only performing the zero-crossing test (with some  > 0 to avoid infinite refinement) or the distance test.Observe that the zero-crossing test alone results in highly adaptive triangles away from the curve, as desired, but overrefines the region along the curve.On the contrary, the distance test alone creates desirable cell sizes and shapes along the curve, but over-refines the rest of the domain.The latter is because the distance test considers the level sets at all levels (not just the zero set).Combining both tests yields well-shaped and adaptive cells along the implicit curve while avoiding over-refinement away from the curve, as seen in Figure 3 (b).We now examine the tests on vector zero sets.We start with a 2-component function  = { 1 ,  2 } in 2D in Figure 9.The implicit curves  −1 1 ,  −1 2 are circles with large radii and are almost tangent at their intersection points ( −1 ).As the curves are nearly flat, checking only the criteria on each curve leads to a coarse grid, where the intersection points are missed by the discretization (see (a)).If we additionally perform the distance test on the intersection points ( −1 ), we obtain a more refined grid that recovers the intersection points (see the arrows in (b)), but the refinement extends to everywhere that the two curves are close by.By also performing the zero-crossing test on the intersections, our algorithm localizes the refinement to the vicinity of the intersection points (see (c)).We show a 3D example of 2-component function in Figure 10.The implicit surfaces  −1 1 ,  −1 2 consist of a sphere and a cylinder, whose intersection curve ( −1 ) has a region of high curvature that is poorly discretized if we only check the refinement criteria on each surface.Adding the zero-crossing test on the intersection curve (with some  > 0) results in densely refined cells everywhere along the curve, whereas performing only the distance test can adapt the level of refinement to the curvature of the curve but at the cost of over-refining nearby regions where the two surfaces are close.By combining both tests on the intersection curve, our algorithm captures the intersection curve without unnecessary refinement.
We next consider a 3-component function example in Figure 11.The implicit surfaces  −1 1 ,  −1 2 ,  −1 3 are three intersecting cylinders whose axes are almost parallel.These cylinders intersect at several near-coinciding curves (zero sets of 2-component subfunctions), which in turn intersect at two points  −1 (red dots in the top view).However, if we only check the refinement criteria on the surfaces and their intersection curves, the intersection points are missed on the discretization (see (a)).This is because both the surfaces and the intersection curves are nearly "flat" and hence a coarse grid resolution is sufficient to approximate them well.The two  intersection points are recovered only after we additionally perform the tests on them, which refines the grid in their vicinity (see (b)).
To further validate our distance test, we measure the Hausdorff distance between a zero set and its discretization on a refined grid.We consider zero sets with analytical forms, including the zero set of a scalar function (a sphere with radius 0.3), a 2-component function (the intersection curve in Figure 10), and a 3-component function (the two intersection points in Figure 11).Although our distance test is based on the one-sided Hausdorff distance, we measure the two-sided (symmetric) Hausdorff distance for validation, and we do so at a range of threshold values.As shown in Figure 12, the distance between the actual zero set and the discretization remains consistently lower than the threshold value used.

Active components
We evaluate the effectiveness of active components, which are used in our refinement criteria to avoid unnecessary refinement for implicit complexes that contain partial zero sets, such as CSG and MI (Section 5.2).
Figure 13 shows a simple MI in 3D, as the Voronoi diagram of three spheres with different radii (also known as the Apollonius Diagram).Without restricting the refinement criteria to only active components (i.e., by treating every component of  as active in every cell), the grid is refined both over the MI patches and along their extensions (see arrows in (a)).These extensions lie on the zero set of the difference between a pair of components, but they are dominated (i.e., having lower values than) by another component.The overrefinement is avoided by considering only active components.
An alternative approach for avoiding unnecessary refinement in discretizing CSG is to ignore empty grid cells that do not intersect with the CSG surface [Duff 1992;Hijazi et al. 2010].However, this approach still considers all components of  in a non-empty cell, even though they do not all contribute to the shape.This can result in over-refinement where a primitive (e.g., the cube faces in Figure 14 (a)) intersects the trimmed-away part of another primitive (e.g., the sphere in Figure 14 (a)).In contrast, our refinement criteria can reduce such interference by considering only active components in a non-empty cell, as shown in Figure 14 (b).

Simplicial refinement
We compare the compactness of the grid between our refinement method (LREB) and the classical longest edge bisection (LEB) method (as implemented in the MFEM package [Anderson et al. 2021]) on an implicit surface example in Figure 15.Observe that LREB consistently produces around a third of the total number of cells produced by LEB at the same threshold.
The compactness of LREB comes at the cost of reduced cell quality, particularly those away from the implicit shape.To evaluate the quality of cells around the shape, we took all 3D shapes that we have used so far in the paper and measured the worst quality of the cells around the shape as the distance threshold  decreases to a very small amount (0.01% of the bounding box size).Here the quality of  Comparing the longest edge bisection (LEB) method and our longest refinable edge bisection (LREB) method on the tear-drop surface (a surface of revolution whose profile is defined by the tear-drop curve in Figure 3), showing a cross-section of the grid refined at  = 0.005 (top; cells containing the geometry are highlighted), the implicit surface extracted from the grid (middle), and the graph of number of cells as a function of the distance threshold  (bottom).
a tetrahedron is measured as the ratio between its in-radius over its out-radius, normalized by the maximum of such ratio (1/3).As shown in Figure 16, the worst quality across all thresholds in all examples is around 0.03.Observe that the graphs for many examples (except for CSG) oscillate rather than trending downward, which suggests the possible existence of a lower bound.

Implementation and performance
Our algorithm is implemented in C++ for  = 2, 3 dimensions without using any parallelization or third-party packages.We use an alternative approach for checking extremity in our zero-crossing test (Equation 9), which we found to be faster than LP due to the small problem size.Recall that the zero-crossing test checks whether the origin of R  is an extremal point in a set of  + 1 points, where  ≤  and  is the number of control points in the cubic Bezier simplex (10 for  = 2 and 20 for  = 3).We enumerate all -tuples of points.For each tuple, we check if the origin lies on a different side of the line ( = 2) or plane ( = 3) formed by the tuple from the remaining points that are not in the tuple.If so, the test returns false.The test returns true when the enumeration successfully completes.All timings are recorded on a PC with a 10-core Intel i9-10900X 3.70Hz CPU and 64GB of RAM.
We evaluate the performance of our algorithm on the IA of 18 implicit spheres with the same radius that are located in a regular pattern (see Figure 17 top-left).Figure 17 (top-right) plots the running time of each part of our algorithm as the distance threshold  decreases.Observe that the running time is dominated by evaluating the function values and gradients at the grid points (red) and checking the criteria (blue), and the latter in turn is dominated by computing the proxy (solid line) and performing the tests on implicit surfaces (short dashes).Note that while the tests on intersection curves and points (medium and long dashes) are more expensive to perform than those on surfaces, those tests are only invoked for cells where there are more than one (two for MI) active components, which are a small fraction of total cells.Our algorithm scales roughly linearly with both the number of function components and the total number of cells, as shown in the two plots in Figure 17 bottom.

More examples
Figure 1 shows examples of all four types of implicit complexes that are considered in this paper, including an IA of 10 implicit surfaces (a), a CSG shape defined by the same set of surfaces (b), an MI as the Voronoi diagram of five lines (c), and the network of junction curves of the MI (d).Observe from the visualization of the grids that our method can effectively adapt the grid resolution to the geometric complexity in each example without unnecessary refinement.
We next show a few challenging examples of discretizing implicit surfaces and complexes to further demonstrate the benefit of our algorithm, particularly its adaptivity and refinement around intersections.
An adaptive grid is particularly useful when the function is expensive to evaluate.As the runtime of grid generation is dominated by function evaluations, reducing the number of grid points (where evaluation takes place) is the key to improving efficiency.Figure 18 gives an example of such a function, which is an Hermite RBF defined by around 1000 points sampled from the Max Planck model (obtained from [Huang et al. 2019]).Note that the evaluation time of a Hermite RBF scales linearly with the number of points.Evaluating this function on a uniform cubic grid of size 100 3 (over 1 million evaluations) took 27.76 seconds, whereas our algorithm took 1.19 seconds to generated a grid with only 21, 286 points (at  = 0.00125).Despite a much more compact grid, the implicit surface computed on our adaptive grid better captures geometric details, such as eyes and ears, than that computed from the uniform grid.
Another scenario that demands adaptivity is when the implicit shape contains features at different scales.While a uniform grid would need an excessive number of cells to capture the smallest feature, an adaptive grid would need much fewer cells as it can  adapt the cell size to the feature's scale.To this end, we created a stress test in Figure 19 where the input function is the SDF (signed distance function) from the surface of three nested "key" shapes at decreasing sizes.Our algorithm at  = 0.0005 creates a grid of 2, 225, 796 cells that enables the discretization to capture all three keys (observe the variable triangle sizes on different part of the keys in the zoom-in view).In comparison, the surface extracted from a uniform tetrahedral grid with 5 × 10 6 cells (obtained from a 100 3 cubic grid) can barely reconstruct the largest key.
Finally, we demonstrate the intersection-awareness of our algorithm on several CSG examples.The shape in Figure 20 is created by subtracting three cylinders from a torus.Since the cylinders are almost tangent to the torus, the shape contains a few very narrow spots (e.g., in the boxed region).As our algorithm adapts the grid to the geometry of both the surface patches and their intersections, it can faithfully capture the shape of each narrow band (see the zoomed-in view in (a)).In contrast, adapting the grid to only the surfaces (i.e., only checking the refinement criteria for scalar zero sets) leads to artifacts at the narrow spots (see the zoomed-in view in (b)).Similar observations can be made in Figure 21, which takes the union of 10 randomly located tori and subtracts it from the union of another 10 tori, and Figure 22, which performs subtraction (a) and union (b) of two copies of the same dog head surface (each defined by an Hermite RBF function) that are slightly shifted from each other.By adapting the grid to intersection curves in addition to the surfaces, our algorithm better captures fine features and produces smoother intersection curves.

DISCUSSIONS
We presented a new method for generating adaptively refined simplicial grids for discretizing implicit surfaces and complexes.Our method adapts the grid resolution to not only surfaces but also their lower-dimensional intersections, thus enabling more accurate discretizations without excessive refinement.Our method does not require interval analysis, and hence can be applied to arbitrary functions, and it can produce grids for several types of commonly used implicit complexes.
The effectiveness of our method largely depends on how well the input function can be approximated by the proxy within each simplex.Our choice of the proxy, the cubic Bezier simplex, may fail to adequately capture the variation of non-trivial functions and thus lead to either insufficient or redundant refinement of the grid.We demonstrate these two failure modes in Figures 23 and 24.In Figure 23, we discretize a high-order algebraic curve in R 2 (the rose  curve) using two domains with different sizes.Our method produces an adequately refined grid in the larger domain (top) but fails to refine beyond the initial grid in the smaller domain (bottom).This is because the proxies in both triangles of the initial grid do not pass the zero-crossing test.In Figure 24, our method over-refines the grid in a region far from the zero-set (near the center).This is because the function there has a high amount of oscillation (see the insert), which causes our proxy to over-estimate the function range.Using  a higher-order proxy may alleviate these issues, but ultimately an interval analysis may be necessary.
While our refinement method (LREB) produces more compact grids than the classical LEB method, it comes at the cost of significantly worsened simplex quality away from the implicit geometry.This makes LREB unsuited for applications that require well-shaped simplices over the entire domain, such as FEA and visualization.Another drawback of LREB is that it needs to explicitly store the vertex coordinates and simplex connectivity.In contrast, the regularity of simplices produced by LEB enables space-efficient data structures that can implicitly encode the geometry and connectivity of the grid (see review in [Weiss and De Floriani 2011]).
There are several venues of future work that we wish to pursue.Our current construction of the proxy function requires evaluating gradients.This can be avoided by first constructing the triangular Lagrange interpolant from values of  at the control points [Farin 1986] and then converting to the Bezier form.We would like to explore refinement criteria for other implicitly defined shapes, such as F-reps [Pasko et al. 1995], extremal features (e.g., ridges and valleys) of unsigned functions, and topological features in scalar or vector functions (e.g., separatrices).Another interesting direction is to extend our method to higher-order grids [Jiang et al. 2021;Khanteimouri and Campen 2023].The extension could leverage the common representation of a typical curved simplex and our proxy function -both as a Bezier simplex.

A CUBIC BEZIER SIMPLEX
We give details of the cubic Bezier simplex and its construction from sampled values and gradients at the simplex vertices.A thorough discussion can be found in [Farin 1986].
We obtain the control values (or the Bezier ordinates)  1 , . . .,   following the "nine parameter interpolant" method described in [Farin 1986], which interpolates given values and gradient vectors at each vertex of .Specifically, let  () and ∇ () be the value and gradient at a vertex , • If   is a vertex of , then   =  (  ).
• If   is a trisector of edge     and closer to   , then

B INVERSION-FREE DISTANCE TEST
We show how the distance test of Equation 15 can be re-formulated without matrix inversions.To further reduce numerical errors when implementing on floating-point numbers, we also remove divisions and square-roots in the calculations.
To demonstrate the improved numerical accuracy of the reformulated test (Equation 16) over the original (Equation 15), we compared both in discretizing the IA of a two-component function in Figure 25.The components  1 ,  2 are identical spherical functions slightly shifted from each other, and hence their gradients  1 ,  2 are almost identical (and hence    is almost singular) along the circle where the two spheres intersect.Compared to the reformulated test (shown in (c)), the original test results in significantly more refinement along the intersection curve (see (b)).This indicates that the original test over-estimates the distance error (left-hand side of Equation 15), which we attribute to inverting the almost-singular   .

Fig. 1 .
Fig. 1.Our method generates simplicial grids (top) that enable adaptive discretization of a variety of implicit shapes defined by vector functions (bottom), such as arrangements of multiple implicit surfaces (a), CSG shapes (b), non-manifold material interfaces (e), and curve networks (d).The grids are adapted to the geometric complexity of not only the surface components but also their intersections (shown as wires on top).

Fig. 2 .
Fig. 2. (a): Two views of a plane (blue, implicit surface of  (, ,  ) =  ) and a nearly flat surface (yellow, implicit surface of  (, ,  ) = 0.1 * ( 2 +  3 − 0.005) − ), whose intersection curve (red) is far from being flat.(b,c): Discretization on a grid adapted only to the two implicit surfaces either cannot capture well the intersection curve (b) or over-refines the surface (c).(d): Discretization on a grid adapted to both the surfaces and their intersection curve (generated by our method) accurately captures the curve without over-refinement.

Fig. 3 .
Fig. 3. Refining a triangular grid using the longest edge bisection (LEB) method (a) and our method (b) for discretizing the tear-drop curve [de Figueiredo et al. 2006] (red, implicit curve of  (, ) =  5 +  4 − 2 2 ).Our method produces much fewer triangles away from the curve, while maintaining the shape and adaptivity of triangles around the curve (gray).

Fig. 4 .
Fig. 4. Illustration of several types of implicit complexes in R 2 .

Fig. 6 .
Fig. 6.Illustration for computing the distance from a point  to the level set  −1 (  ( ) ) for a scalar function in 1D (left) and a two-component function in 3D (right).
Figure 6 (left) gives an illustration for the 1-dimensional case.Now consider  > 1.The vector level set  −1 (  ()) is the intersection of  scalar level sets  −1  (   ()), each being a hyperplane orthogonal to a gradient vector   = ∇  (see illustration in Figure 6 (right)).Let  be the vector from  to its nearest point on this intersection.Since the signed distance from each scalar level set to  is (   () −   ())/|  |, we have the identities 1 , of scalar zero sets ( = 1), the inside part of the max operator in Equation 15 reduces to |(  −   )|/||, where the   and

Fig. 8 .
Fig. 8. Refined grids and resulting discretization for the tear-drop curve by performing only the zero-crossing test (a) or only the distance test (b).The result of using both tests is shown in Figure 3 (b).

Fig. 9 .
Fig. 9.The refined grids and resulting discretizations of the IA of a 2component function  = {  1 ,  2 } in 2D created using different settings: checking only the criteria on each scalar level set  −1  (a), also performing the distance test on the vector level set  −1 (b), and performing both zero-crossing and distance tests on  −1 (c).Each curve segment in the IA is assigned a different color, and arrows point to  −1 .

Fig. 10 .
Fig. 10.Discretization of the IA of a 2-component function  = {  1 ,  2 } in 3D (top) on grids refined under three settings (showing zoomed-in views): checking only the criteria on each  −1  (a), adding only the zero-crossing test on  −1 (b), adding only the distance test on  −1 (c), and performing both tests on  −1 (d).Each patch in the IA is assigned a different color.

Fig. 11 .
Fig. 11.Discretization of the IA of a 3-component function  in 3D (top;showing two views) on grids refined without checking the refinement criteria on the 3-vector zero set  −1 (but checking all scalar and 2-vector zero sets) (a), and with the criteria turned on (b).

Fig. 12 .
Fig. 12.The two-sided Hausdorff distance between the analytical and discrete zero sets of a scalar function (a sphere of radius 0.3), a 2-component function (Figure 10), and a 3-component function (Figure 11), over decreasing refinement threshold .

Fig. 13 .
Fig. 13.Comparing discretization of an MI (the Voronoi diagram of three spheres) on a grid refined by applying the refinement criteria to all components (a) and only active components (b).Observe the over-refinement in (a), pointed by arrows.

Fig. 14 .
Fig. 14.Comparing discretization of a CSG shape, consisting of a cube and a truncated sphere (the complete sphere is shown in transparency), on a grid refined by only avoiding empty cells (a) and by considering only active components in each cell (right).Observe the over-refinement in (a) where the sphere intersects the cube.

Fig
Fig.15.Comparing the longest edge bisection (LEB) method and our longest refinable edge bisection (LREB) method on the tear-drop surface (a surface of revolution whose profile is defined by the tear-drop curve in Figure3), showing a cross-section of the grid refined at  = 0.005 (top; cells containing the geometry are highlighted), the implicit surface extracted from the grid (middle), and the graph of number of cells as a function of the distance threshold  (bottom).

Fig. 16 .
Fig. 16.Worst quality (measured by normalized ratio between in-radius and out-radius) among tetrahedra containing the implicit shape generated by our LREB method as a function of decreasing , for all 3D examples shown so far in the paper.

Fig. 17 .
Fig.17.Timing of our algorithm on an IA consisting of intersecting spheres (top-left; using  = 0.0001), plotted against decreasing threshold  (topright; for 18 spheres), increasing number of spheres (bottom-left; at  = 0.0001), and increasing number of tetrahedra (bottom-right; for 18 spheres by varying ).The timing is broken down into evaluation of function values and gradients (red), cell splitting (green), and checking the refinement criteria, including computing the proxy and quantities needed for robust evaluation (Appendix B; solid blue), and performing the zero-crossing and distance tests for subfunctions with different number of components  (dashed blue).

Fig. 18 .
Fig. 18.Comparing the implicit surfaces of a Hermite RBF function with 1000 sites (left) extracted on a uniform tetrahedral grid (tessellation of a 100 3 cubic grid) and on our adaptively refined grid ( = 0.00125).

Fig. 19 .
Fig. 19.Comparing the implicit surface of an SDF extracted on our adaptively refined grid ( = 0.0005) and on a uniform tetrahedral grid (tessellation of a 100 3 cubic grid).The SDF is created by Shadertoy user Flopine [2020], released under CC BY-NC-SA 3.0 License.

Fig. 20
Fig. 20.A CSG shape defined by subtracting three cylinders from a torus, discretized on a grid that is adapted to both the surfaces and their intersection curves (a) or adapted only to surfaces (b).

Fig. 21
Fig. 21.(a): A CSG shape defined by random subtraction and addition of 20 tori, discretized on a grid that is adapted to both the surfaces and their intersection curves (a) or adapted only to surfaces (b).

Fig. 22 .
Fig. 22. (a,b): The result of subtraction (a) or addition (b) between two slightly shifted dog head surfaces (each defined by an Hermite RBF function), discretized on our adaptively refined grid.(c,d): Comparing discretizations with or without adaptive refinement for the intersection curves.

Fig. 23 .
Fig. 23.Discretizing the 4-pedal rose curve (zero set of  (, ) = 4 2  2 − ( 2 +  2 ) 3 ) using a larger domain (top) and a smaller domain (bottom).The curves are visualized on the left as the intersection (red) between the height surface of  (, ) (blue) and the zero plane (gray).The dotted curve in lower-right indicates where the rose curve should be, but missed by the discretization, due to an under-refined grid (made up of two triangles).

Fig. 24 .
Fig. 24.Discretizing a wavy circle (zero set of  (,  ) = 0.2 0.1 sin(16 ) +  2 − 1, where ,  are polar coordinates).Left: the curve visualized as the intersection (red) between the height surface of the function (blue) and the zero plane (gray).The zoom-in shows the region around the origin.Right: the refined grid and discretization (red).The arrow points to an over-refined region at the origin, far away from the curve.
the centroid of a triangle with vertices  and edge trisectors , then control values   at the triangle centroids do not affect the interpolation of the given values and gradients.The choices made above have the property that the resulting interpolant reproduces all polynomial  up to quadratics.

Fig. 25 .
Fig. 25.IA of a two-component function, whose components are slighted shifted spherical distance functions, discretized on a grid refined by the original distance test in Equation 15 (b) and the reformulated test in Equation 16 (c).An exploded view of the IA is shown in (a). .