Isotopic Arrangement of Simple Curves: an Exact Numerical Approach based on Subdivision
IIsotopic Arrangement of Simple Curves:an Exact Numerical Approach based on Subdivision
Jyh-Ming Lien
Department of Computer ScienceGeorge Mason UniversityFairfax, VA 22030, [email protected]
Vikram Sharma
The Institute of Mathematical Sciences, HBNI, Chennai, [email protected]
Gert Vegter
Johann Bernoulli Institute for Mathematics and Computer ScienceNijenborgh 9, 9747 AG GroningenThe [email protected]
Chee Yap
Department of Computer ScienceCourant Institute of Mathematical SciencesNew York University. New York, [email protected]
Abstract
This paper presents the first purely numerical (i.e., non-algebraic) subdivision algorithm forthe isotopic approximation of a simple arrangement of curves. The arrangement is “simple” inthe sense that any three curves have no common intersection, any two curves intersect transver-sally, and each curve is non-singular. A curve is given as the zero set of an analytic function f : R → R , and effective interval forms of f, ∂f∂x , ∂f∂y are available. Our solution generalizes theisotopic curve approximation algorithms of Plantinga-Vegter (2004) and Lin-Yap (2009).We use certified numerical primitives based on interval methods. Such algorithms havemany favorable properties: they are practical, easy to implement, suffer no implementationgaps, integrate topological with geometric computation, and have adaptive as well as localcomplexity.A version of this paper without the appendices appeared in [9]. a r X i v : . [ c s . C G ] S e p Introduction
We address problems in computing approximations to curves and surfaces. Most algebraic algo-rithms for curve approximation begin by computing a combinatorial object K first. To compute K ,we typically use algebraic projection (i.e., resultant computation), followed by root isolation andlifting. But most applications will also require the geometric realization G . Thus we will need aseparate (numerical) algorithm to compute G . This aspect is typically not considered by algebraicalgorithms.In this paper, we describe a new approach for computing curve arrangements based on purelynumerical (i.e., non-algebraic) primitives. Our approach will integrate the computation of the com-binatorial ( K ) and geometric ( G ) parts. This leads to simpler implementation. Our numericalprimitives are designed to work directly with arbitrary precision dyadic (BigFloat) numbers, avoid-ing any “implementation gap” that may mar abstract algorithms. Furthermore, machine arithmeticcan be used as long as no over-/underflow occurs, and thus they can serve as efficient filters [3].We now explain our specific problem, and illustrate the preceding notions of K and G . By a simple curve arrangement we mean a collection of non-singular curves such that no three ofthem intersect, and any two of them intersect transversally. The simple arrangement of three ormore curves can, in some sense, be reduced to the case of two curves (see the Final Remarks).Let F : R → R , where F ( x, y ) = ( f ( x, y ) , g ( x, y )) is a pair of analytic functions. It genericallydefines two planar curves S = f − (0) ⊆ R and T = g − (0). We call F = 0 a simple system of equations if { S, T } is a simple curve arrangement. Throughout this paper, F = ( f, g ) will befixed unless otherwise indicated. Figure 1 illustrates such an arrangement for the curves definedby f ( x, y ) = y − x and g ( x, y ) = x + y −
1. The concept of hyperplane arrangement is highlyclassical in computational geometry [5]. Recent interest focuses on nonlinear arrangements [2]. (a) (
S, T ) arrangement y = x x + y = 1 (b) ( S, T )-decomposition K ∗ x + y = 1 y = x (c) cell complex K Figure 1: Arrangement of two curves, y = x and x + y = 1Our basic problem is the following: suppose we are given an (cid:15) > B ⊆ R ,called the region-of-interest or ROI, which is usually in the shape of an axes-aligned box. We1ant to compute an (cid:15) -approximation to the arrangement of the pair ( S, T ) of curves restricted to B . This will be a planar straightline graph G = ( V, E ) where V is a finite set of points in B and E is a set of polygonal paths in B . Each path e ∈ V connects a pair of points in V , and nopath intersects another path or any point in V (except at endpoints). Moreover, E is partitionedinto two sets E = E S ∪ E T such that ∪ E T (resp., ∪ E S ) is an approximation of T (resp., S ).The correctness of this graph G has two aspects: (A) topological correctness, and (B) geometriccorrectness. Geometric correctness (B) is easy to formulate: it requires that the set ∪ E S ⊆ B is (cid:15) -close to S in the sense of Hausdorff distance: d H ( S, ∪ E S ) ≤ (cid:15) . Similarly, the ∪ E T is (cid:15) -close to T .If we specify (cid:15) = ∞ , then we are basically unconcerned about geometric closeness.Topological correctness (A) is harder to capture. One definition is based on the notion of “celldecomposition”. A (cell) decomposition of B is a partition K ∗ of B into a collection of setscalled cells, each c ∗ ∈ K ∗ homeomorphic to a closed i -dimensional ball ( i ∈ { , , } ); we call c ∗ an i -cell and its dimension is dim( c ∗ ) = i . If b ∗ is an i -cell and c ∗ an ( i + 1)-cell, we say b ∗ bounds c ∗ if b ∗ is contained in the boundary ∂c ∗ of c ∗ . Call K ∗ an ( S, T )-decomposition of B if the set( S ∪ T ) ∩ B is a union of some subset of 0- and 1-cells of K ∗ . A ( S, T )-decomposition is illustratedin Figure 1(b).A cell complex K is an (abstract) set such that each c ∈ K has a specified dim( c ) ∈ { , , } together with a binary relation B ⊆ K × K such that ( b, c ) ∈ B implies dim( b ) + 1 = dim( c ). Wesay that the decomposition K ∗ is a realization of K , or K is an abstraction of K ∗ , if there is a1-1 correspondence between the cells c ∗ of K ∗ with the elements c ∈ K such that dim( c ∗ ) = dim( c ),and moreover the relation ( b, c ) ∈ B iff b ∗ bounds c ∗ in K ∗ . Figure 1(c) shows the abstraction K of the decomposition in Figure 1(b).Our algorithmic goal is to compute a planar straightline graph (PSLG for short [19]) G = ( V, E )which approximates (
S, T ) in a box B . Such a graph G naturally determines a decomposition K ∗ ( G ) of B as follows: the set of 0-cells is V , the set of 1-cells is E and the set of 2-cells is simplythe connected components of B \ ( V ∪ ( (cid:83) E )). Finally, we say G is topologically correct if thereexists an ( S, T )-decomposition K ∗ such that K ∗ and K ∗ ( G ) are realizations of the same abstractcell complex. ¶
1. Towards Numerical Computational Geometry.
The overall agenda in this line ofresearch is to explore new modalities for designing geometric algorithms. We are interested inexploiting weaker numerical primitives that are only complete in a certain limiting sense. Unliketraditional exact algorithms, our algorithms must strongly interact with these weaker primitives,and exploit adaptivity. The key challenge is to achieve the kind of exactness and guarantees thatis typically missing in numerical algorithms. See [26] for a discussion of “numerical computationalgeometry”.In the algebraic approach, one must compute the abstract complex K before the approximateembedded graph G . Indeed, most algebraic algorithms do not fully address the computation of G . In contrast to such a “decoupled” approach, our algorithm provides an integrated approachwhereby we can commence to compute G (incrementally) even before we know K in its entirety.Ultimately, we would be able to determine K exactly — this can be done using zero bounds as in[25, 4]. The advantage here is that our integrated approach can cut off this computation at anydesired resolution, without fully resolving all aspects of the topology. This is useful in applicationslike visualization.Unlike exact algebraic primitives, our use of analytic (numerical) primitives means that ourapproach is applicable to the much larger class of analytic curves. Numerical algorithms are rel-atively easy to implement and have adaptive as well as “local” complexity. Adaptive means thatthe worst case complexity does not characterize the complexity for most inputs, and local means2he computational effort is restricted to ROI.One disadvantage of our current method is that it places some strong restrictions on the classof curve arrangements: the curves must be non-singular with pairwise transversal intersections inthe ROI. In practice, these restrictions can be ameliorated in different ways. The complete removalof such restrictions is a topic of great research interest.The algorithms in this paper fall under the popular literature on Marching-cube type algorithms[16]. There are many heuristic algorithms here which are widely used. The input for these algo-rithms can vary considerably. E.g., Varadhan et al. [24, 23] discuss input functions F : R → R that might be a discretized function, or a CSG model or some polygonal model – each assumptionhas its own exactness challenge. All current exact algorithms for curve arrangements are based on algebraic projection, i.e., theyneed some resultant computation. The disadvantage of projection is the large number of cells:even in relatively simple examples, the graph can be large as seen as Figure 1(c). For manyapplications, the 2-cells may be omitted, but the graph remains large. There are several knowntechniques to reduce this (double-exponential in dimension) explosion in the number of cells. In thispaper, we avoid cell decomposition, but base our topological correctness on the concept of isotopy.Our algorithm uses the well-known subdivision paradigm, and produces a subdivision of the inputdomain into boxes. Figure ?? illustrates the form of output from our subdivision algorithm usingour previous example of y = x and x + y = 1. The number of subdivision boxes tend to beeven more numerous than cells in the decomposition approach. But these numbers are not directlycomparable to number of cells for three reasons: (1) Subdivision boxes are very cheap to generate.(2) Most of these boxes can be instantly discarded as inessential for the final output (we keep themfor visualization purposes). (3) Unlike cells, our subdivision boxes play a double role: they are usedfor (A) topological determination as well as (B) in determining geometric accuracy.The approach of this paper has previously been successfully applied to the isotopic approxima-tion of a single non-singular curve or surface by Plantinga and Vegter [18, 17] and Lin and Yap[11, 10]. The current paper is a non-trivial extension of these previous works.We now define the notion of isotopy for arrangements. For our problem on arrangements, weneed to extend the standard definitions of isotopy. Suppose
S, T ⊆ R are two closed sets and (cid:15) >
0. First recall that S and T are (ambient) isotopic if there exists a continuous mapping γ : [0 , × R → R (1)such that for each t ∈ [0 , γ t : R → R (with γ t ( x, y ) = γ ( t, x, y )) is a homeo-morphism, γ is the identity map, and γ ( S ) = T . If, in addition, d H ( S, T ) ≤ (cid:15) (where d H is theHausdorff distance on closed sets) we say that they are (cid:15) -isotopic . We will write S (cid:15) (cid:39) T (via γ )in this case. Note that we may omit mention of (cid:15) , in which case it is assumed that (cid:15) = ∞ . The figure is not produced by the algorithm of this paper because the implementation is currently underway.Instead, it is produced by the Cxy Algorithm for approximating a non-singular curve [11], using the input curve fg = 0. Thus the intersection points are singularities which the Cxy algorithm cannot resolve, but this does notprevent its computation to some cut-off bound. Also, the Cxy algorithm does not know which part of the arrangementis the f -curve and which is the g -curve.
3e now generalize this to arrangement of sets. Let S = ( S , . . . , S m ) and T = ( T , . . . , T m ) betwo sequences of m closed sets. For each non-empty subset J ⊆ { , , . . . , m } , let S J denote theintersection ∩ i ∈ J S i . Similarly for T J . We say that S and T are isotopic if there exists a continuousmapping γ as in (1) such that for each non-empty subset J ⊆ { , , . . . , m } , we have S J (cid:15) (cid:39) T J (via γ ) . We also call γ an isotopy from S to T . For simple curve arrangements, the critical problem tosolve is the case m = 2. We assume the two curves S , S are restricted to a region or box B . Ourbasic problem is to compute a pair of curves ( T , T ) such that( T , T ) (cid:15) (cid:39) ( S ∩ B, S ∩ B ) . (2)The approximations ( T , T ) produced by our algorithms will be piecewise linear curves. See [1] fora general discussion of isotopy of the case m = 1. In Appendix A, we provide the necessary definitions; these are consistent with the terminology inthe related work [11]. For now, we rely on common terms that are mostly self-explanatory. ¶
2. Box Complexes and Subdivision Trees.
Our fundamental data structure is a sub-division tree T rooted in some box B . In 2-D, T is the well-known quad-tree and B is arectangle. Each internal node of T has four congruent children. The boxes of a subdivision treeare non-degenerate (i.e., 2-dimensional). They need not be squares, but for the correctness of ouralgorithm, their aspect ratios must be ≤
2. For any region R ⊆ R , we define a subdivision of R to be a set S = { R , . . . , R n } of subregions such that R = ∪ ni =1 R i and the interiors of R i ’s arepairwise disjoint. If each R i is a box, we call S a box subdivision . The box subdivision is a box complex if for any two adjacent boxes B, B (cid:48) ∈ S , their intersection ∂ ( B ) ∩ ∂ ( B (cid:48) ) is side ofeither B or B (cid:48) . Clearly, the set S of leaf boxes of T forms a box complex of B . But in this paper,we need to consider a more general subdivision of B that is obtained as the leaf boxes of a finitenumber of subdivision trees. A segment of a box complex S is the side of a box of S that does notproperly contain the side of an adjacent box. Therefore every side of a box of S is a finite unionof segments. We say the box complex S is balanced if every side is either a segment or the unionof two segments. A segment is called bichromatic w.r.t. a curve S if S has different signs on theendpoints of the segment; otherwise call it monochromatic .Although ( S, T ) is simple, we need to consider degeneracies induced by a subdivision S : wesay ( S, T ) is S -regular if S ∪ T does not intersect any corner of a box in S . This can be effectivelyachieved by an infinitesimal perturbation of S and T using a trick in [18]: when we evaluate thesign of f at a box corner, we simply regard a 0 sign to be +1. ¶
3. Normalization.
Consider an isotopy of the arrangement (
S, T ) into another arrangement( S (cid:48) , T (cid:48) ). Let us write ( S, T ) t for the arrangement at time t ∈ [0 ,
1] during this transformation. Thus(
S, T ) = ( S, T ) and (
S, T ) = ( S (cid:48) , T (cid:48) ). The isotopy is said to S -regular provided, for all t ∈ [0 , S, T ) t is S -regular. We say that ( S, T ) is S -normalized if:(N0) ( S, T ) is S -regular.(N1) Each subdivision box B of S contains at most one point of S ∩ T .(N2) Let X ∈ { S, T } . Then X intersects each segment of S at most onceCall ( S (cid:48) , T (cid:48) ) a S -normalization of ( S, T ) if there exists a S -regular isotopy from ( S, T ) to( S (cid:48) , T (cid:48) ) such that ( S (cid:48) , T (cid:48) ) is S -normalized. Our algorithm will construct an S -normalization ( S (cid:48) , T (cid:48) )of ( S, T ). 4
4. Box Predicates.
We will use a variety of box predicates. These predicates will determinethe subdivision process. Typically, we will keep subdividing boxes until some Boolean combinationof some box predicates hold.Let h : R → R be any real function. Recall (Appendix A) that we assume an interval formu-lation of h denoted h : R → R where R denotes the set of closed intervals and R can beviewed as the set of boxes. We introduce a pair of box predicates denoted C h and C h , defined as C h ( B ) ≡ (cid:54)∈ h ( B ) ,C h ( B ) ≡ (cid:54)∈ ( h x ( B )) + ( h y ( B )) . (cid:27) (3)Note that C h as taken from Plantinga-Vegter, where the interval operation I is defined as { xy : x, y ∈ I } and not { x : x ∈ I } . An alternative to C h would be the weaker C hxy predicate fromLin-Yap [11], but the corresponding algorithm would would be more involved. So for now, we focuson the C h predicate. We classify boxes using these predicates: • Box B is h -excluded if it satisfies C h ( B ). • Box B is h -included if it fails C h ( B ) but satisfies C h ( B ). • Box B is resolved if it satisfies the predicate( C f ∨ C f ) ∧ ( C g ∨ C g ) . (4) • Box B is excluded if it satisfies C f ∧ C g . Note that excluded boxes are resolved. • Box B is a candidate if it is resolved but not excluded. • Candidate boxes can be further classified into three subtypes: f -candidates are those thatare f -included but g -excluded, g -candidates is similarly defined, and f g -candidates arethose that are f - and g -included. ¶
5. Root Boxes.
We define a root box to be any box B where B ∩ S ∩ T has exactly one point.We next consider two predicates that will allow us to detect root boxes. One is the Jacobiancondition , JC( B ) ≡ / ∈ det( J F ( B ))where J F ( B ) is the Jacobian of F = ( f, g ) evaluated on B . If JC( B ) holds, then B has at mostone root of f = g = 0, The other is the Moore-Kioustelidis condition
MK( B ) [14] which can beviewed as a preconditioned form of the famous Miranda Test [8]; for other existence tests based oninterval arithmetic see [6]. If MK( B ) holds, then B has at least one root of f = g = 0. We providethe details for this predicate in Appendix B; see (9). Therefore, when JC( B ) and MK( B ) holds,we know that B is a root box. The use of Miranda’s test combined with the Jacobian conditionhas been used earlier to isolate the common roots [12]. What is new in this paper is its applicationto the simple curve arrangement problem. Our algorithm will produce a graph G = ( V, E ) where vertices v ∈ V are points in R and edgesare line segments connecting pairs of vertices. Moreover, each edge E will be labeled as an S -edgeor a T -edge. The union of these edges will provide a polygonal (cid:15) -approximation of ( S, T ). We nowgive an overview of the issues and solution.First, we describe how the vertices of V are introduced.(V0) We introduce a vertex in the center of a root box B .5V1) We evaluate f, g at the endpoints of segments of B . If h ∈ { f, g } is bichromatic on asegment of B , then we must introduce an h -vertex somewhere in the segment. In a balancedsubdivision, an S -normalized pair ( S (cid:48) , T (cid:48) ) of curves has at most two h -vertices on an edge ofa box B .(V2) Introducing vertices on the edges of a box B is straightforward if B is an f -candidate or a g -candidate. When B is a f g -candidate, we may have an edge e containing both a f -vertexand a g -vertex. In the next section we will show how to find the relative order of these twovertices.Next we discuss how to introduce the edges E , which are line segments completely containedin a box. • If B is a root box, we just connect the vertex at its midpoint cen( B ) to each of the verticeson the edges of B . There will be exactly two f -vertices and two g -vertices. • If B is a f -candidate or g -candidate, then the connection is trivial in the regular case. Inthe balanced case, the rules from the previous work of Plantinga-Vegter [18] assures us of thecorrect connection. • If B is a f g -candidate, but not a root box, we know that the f -segment and g -segment willnot intersect. Some f g -candidates need global information to resolve them: when there aretwo edges where each edge contains both an f - and a g -vertex. Their relative order must bedetermined globally from root boxes or from boxes where their relative order is known. Wewill show how to propagate this information in § Suppose ( S (cid:48) , T (cid:48) ) is the normalization of ( S, T ) relative to the box B , i.e., ( S (cid:48) , T (cid:48) ) is an isotopictransformation of ( S, T ) which respects the four corners of B . We now determine the isotopy typeof ( S (cid:48) , T (cid:48) ) in a root box B . The possible combinatorial types fall under one of the 8 patterns asshown in Figure 2. We put them in three groups (I, II, III) for our analysis. (IIIc)(Ib) (Ic) (IIb) (cid:9) (cid:9)⊕⊕(cid:9) (cid:9)⊕(cid:9) f = 0 : g = 0 :KEY:(Ia) (IIIa) (IIIb) (IIa) Figure 2: Local intersection patterns of the normalized curves ( S (cid:48) , T (cid:48) )Following the standard Marching Cube technique, we evaluate the sign of the functions f, g atthe four corners of B . If f has different signs at the endpoints of an edge e of B , then we mustintroduce an f -vertex somewhere in the interior of e . Our normalization assumptions imply thatthere are either zero or two f -vertices on the boundary of B . We treat g similarly. Our aim is toconnect the two f -vertices, the two g -vertices, and a point in the center of the box which representsthe common root with line segments such that the graph G obtained is an isotopic approximationof ( S (cid:48) ∩ B , T (cid:48) ∩ B ). There is a subtlety: the method exploits “local non-isotopy” [18, 11], meaningthat we do not guarantee that S ∩ B is isotopic to the segment introduced to connect two f -vertices.6owever, the graph G will be locally isotopic to the normalized curves ( S (cid:48) , T (cid:48) ), i.e., G ∩ B is isotopicto ( S (cid:48) ∩ B, T (cid:48) ∩ B ) in each subdivision box B .The issue before us is the relative placements of an f -vertex and g -vertex in case they bothoccur in e ; e.g., the patterns in group II in Figure 2. The main result of this section is the following. Theorem . Let B be a root box that satisfies MK( B ) . Then the signs of f and g at each of thefour corners of B determine the combinatorial type of the normalized curves S (cid:48) , T (cid:48) in B . Moreover,these combinatorial types fall under one of the five types in Groups II and III in Figure 2. The main idea of the proof is that if MK( B ) holds for a box B then there exists an edge e of B such that either f ( e ) > cg ( e ), or g ( e ) > cf ( e ), for some c >
0. Given such an e , we can find therelative order of the f -vertex and g -vertex on e . See Appendix C for details of the proof. By an aligned box we mean one that can be obtained as a node of a subdivision tree rooted at theregion-of-interest (ROI) B ; otherwise, it is said to be non-aligned . For instance, in Figure 3(a),let the box with corners p, q, r, s be B . Then the figure shows the four children of B , which arealigned, as well as the non-aligned box (1 / B whose corners are p (cid:48) , q (cid:48) , r (cid:48) , s (cid:48) . Note that (1 / B can be obtained as the union of aligned boxes. We are interested in non-aligned boxes that can beobtained as a finite union of aligned boxes. In the simplest case of non-alignment, a box B is saidto be half-aligned if it is equal to the union of congruent aligned boxes of size w ( B ) /
2. Thus if B is aligned then both (1 / B and 2 B are half-aligned. B B B B B (b) Standard Subdivision of 8 Bpsr qq p s r (a) Half-aligned (1 / B = ( p q r s ) Figure 3: (a) B = ( pqrs ) is aligned, (b) 2 B is a root box.In most subdivision algorithms, it is enough to work with aligned boxes. But to treat rootboxes, we see an essential need to work with non-aligned boxes. The reason is that if we apply theMoore-Kioustelidis predicate to aligned boxes, non-termination may occur when a root of F lieson the boundary of an aligned boxes. But such roots can be detected in the interior of non-alignedboxes. This issue is often ignored in the literature, but it needs to be properly treated in exactalgorithms. Some discussions may be found in Stahl [22] and Kamath [7]; in the univariate case, asolution is suggested by Rote [20] for splines.Therefore, given an aligned box B , we provide a procedure to detect if 2 B is a root box. Weconsider the nested sequence of boxes B ⊂ B ⊂ B ⊂ B as illustrated in figure 3(b). Our goal isto detect 2 B as a root box, but because of alignment issues, we must also treat the larger box 8 B which is called the extended root box corresponding to B .We construct the following standard subdivision of 8 B , denoted Std( B ), into sub-boxes:7 Subdivide 6 B into 9 boxes, each congruent to 2 B (indeed, 2 B is one of these 9 boxes). • The annular region 8 B \ B is partitioned into 28 boxes, each congruent to B . These arecalled the ring boxes .See Figure 3(b) for illustration. Note that Std( B ) is balanced. None of the subdivision boxes arealigned, but the ring boxes are half-aligned. ¶
6. Conforming Subdivisions.
Let Π be a subdivision of a region R . A box B (cid:48) ∈ Π is a boundary box of the subdivision if ∂B (cid:48) intersects ∂R . In the following definitions, we fix a region R ⊆ B and fix a box B such that 8 B ⊆ R . Also let k ≥ for R \ B is called externally k -conforming for B if it has three properties:Π is balanced, the union Π ∪ { B } is a box complex, and for each box B (cid:48) ∈ Π , if B (cid:48) is adjacentto 8 B then w ( B (cid:48) ) = w ( B ) / k . A subdivision Π of 8 B is called internally k -conforming for B ifΠ is balanced, and for every boundary box B (cid:48) of Π , w ( B (cid:48) ) = w ( B ) / k − . Note for instance thatif Π is the standard subdivision of 8 B , then it is internally 1-conforming for B . Below we showhow to achieve subdivisions of 8 B that is internally k -conforming for B for k ≥
2. The followingis immediate: If Π is externally k -conforming for B , and Π is internally k -conforming for B ,then their union Π ∪ Π is a balanced subdivision of R . Note that if k > ∪ Π may cause the edges of a root box 2 B to split into two segments(but not more); see Figure 4. This can be handled by a case analysis similar to Theorem 1 basedon Lemma 7. An alternative approach is to replace 8 B by 10 B which would have an extra ring ofboxes congruent to B . In this case, we can handle any k > B . This gives a simple and effective solution. ¶
7. Strong Root Isolation.
Suppose 2 B is a root box. We say 2 B is strongly isolated if thefollowing conditions hold • (P1) The following four predicates hold: C f (8 B ) , C g (8 B ) , JC(6 B ) , MK(2 B ). • (P2) F = ( f, g ) has no roots in the annulus 8 B \ B .The predicates in (P1) ensures that 2 B is a root box. It is not hard to see that if 2 B contains aroot of F and is sufficiently small, then properties (P1) and (P2) will hold. The reason for MK(2 B )(not just MK( B ) is to ensure that we test the Moore-Kioustelidis predicate on overlapping boxes,so that roots on the boundary of an aligned box B will appear in the interior of 2 B . The reasonfor JC(6 B ) instead of JC(2 B ) is that there can be two boxes 2 B and 2 B (cid:48) such that both of themsatisfy MK-test and they overlap. The test JC(6 B ) ensures that if there are two such boxes thenthey correspond to the same root, and so discard one of them. ¶
8. Root Refinement:
Let B be an aligned box from the subdivision queue such that 2 B is aroot box. We give a subroutine to refine such a root box 2 B . It it important that in our refinementmethod all the sub-boxes remain dyadic boxes, assuming the input boxes are dyadic. The idea isto cover 2 B with a covering of aligned boxes, which must be of size w ( B ) /
2, and check whetherMK-test holds for the doubling of any of these 16 boxes. If not, then subdivide these boxes andcontinue recursively with the f g -candidates. See Appendix A for more details.
Our overall algorithm begins with the (trivial) subdivision tree T rooted at the ROI B but withno other nodes. The algorithm amounts to repeatedly expansion of the candidate leafs in T until avariety of global properties hold. We given an overview of the algorithm in a sequence of 9 stages ;see Appendix C. 8
9. Stage I: Resolution Subdivision
The high level description of this stage is easy: keepexpanding any leaf B of T that is not resolved (see (4)). Recall that resolved boxes are eitherexcluded or candidates. As each box is resolved, it is placed in one of the following four queues: Q for excluded boxes, Q f for f -candidates, Q g for g -candidates, and Q fg for f g -candidates Besidesthese four global queues, we also use these additional queues: Q JC , Q MK , Q Root correspondingroughly to boxes that satisfies the JC and MK predicates, or are found to be root boxes.
The boxesin all the queues are always aligned boxes . ¶
10. Stage II: Jacobian Stage.
Remove a box B from Q fg and do the following: If JC(6 B )holds then put B into Q JC , otherwise, subdivide B and distribute the children into Q , Q f , Q g , Q fg . ¶
11. Stage III: MK Stage.
For every box B ∈ Q JC we subdivide it until either we find asub-box B (cid:48) such that MK(2 B (cid:48) ) holds, or we have identified all sub-boxes as one of Q , Q f , Q g , Q fg . ¶
12. Stage IV: Strong Root Isolation Stage
We assume that Q MK is a priority queue, whereboxes are popped starting from the largest size. For each such box B check whether 8 B is disjointfrom 8 B (cid:48) , for all its neighbors B (cid:48) ; if not then replace B with RefineRoot( B ). We now have obtaineda queue Q Root containing root boxes for all the roots in ROI. The next step is to externally conformStd( B ) with the rest of the subdivision tree T . ¶
13. Stage V: Pruning T In this stage we will turn OFF some leaf boxes in On ( T ) dependingon how they interact with the extended root boxes 8 B . The aim is to “blackout” the 8 B regions fromROI, and ensure that the boxes abutting it are all aligned boxes. Let B (cid:48) be the great-grandparentof B in T . Then we get the list of leaf boxes that cover the interior of B (cid:48) and another list ofboxes that are its neighbors. For each box B tmp in these lists, we turn it OFF if it is contained in8 B ; if it overlaps 8 B then we subdivided it and proceed with its children. Let T (cid:48) be the resultingsubdivision tree. ¶
14. Stage VI: Balancing and Externally Conforming
Recall the standard balancingprocedure for a subdivision T of a region B from the appendix. We will construct a balancedand externally conformal subdivision of B \ ∪ i B i , where 8 B i ’s are pairwise disjoint extended rootboxes. For each box 8 B i , we add a conceptual box to T (cid:48) , with depth either one more than itssmallest neighbor, or if all the neighbors of 8 B are larger than w ( B ) then one more than the depthof B in T . Call the standard balancing procedure on the modified T (cid:48) . By Lemma 3, we will getthe desired subdivision; after balancing the boxes bordering 8 B will all be of the same size, namely w ( B ) / k , for some k ≥ ¶
15. Stage VII: Internally Conforming Extended Root Boxes
Consider any extendedroot box 8 B and its standard subdivision Std( B ). Given a k > B ). Note that since k > B ). To get a balanced conformal subdivisionof Std( B ), we initialize a priority queue Q with all the boxes on the exterior of 8 B (all of them areof the same size) and the 37 boxes in Std( B ). Then we initiate the standard balancing procedureon Q . See Figure 4(c) for an illustration of this procedure; the box B (cid:48) has width w ( B ) /
8. We dothis balancing step for each of the extended root boxes 8 B . The union of these subdivisions withthe balanced subdivision of B \ ∪ i B i gives us a balanced subdivision of B , our ROI. ¶
16. Stage VIII: PV-Construction
For each box in Q f , connect its two f -vertices with aline segment; do the same for boxes in Q g . For each box in Q Root place a vertex at its center andconnect the two f -vertices and the two g -vertices with this vertex according to the cases shown inGroups II and III. of Figure 2. At the end of this stage, the only queue that remains unprocessedis Q fg . The next stage resolves these boxes. 9
17. Stage IX: Resolving Ambiguous f g -candidates
We call an f g -candidate box ambigu-ous if they have the same set of bichromatic segments; otherwise, call the box unambiguous . Bydefinition, boxes where f and g do not share a bichromatic segment are unambiguous. However,some ambiguous boxes can be made unambiguous locally. From Theorem 1 we know that ambigu-ous root boxes can be made unambiguous. Also, boxes where the two shared bichromatic segmentsare on adjacent edges can be made unambiguous by repeated subdivisions of the edges until wereach a segment in one of the edges that is bichromatic for one curve and monochromatic for theother; this will happen along one of the edges since both C f and C g hold. A similar approachworks to resolve ambiguous boxes that share an edge with B and a common bichromatic segmentis on this edge, because by assumption boundary of B does not contain a root of f, g . From theseunambiguous boxes, we propagate the ordering of the f -vertex and g -vertex on the shared edge totheir ambiguous neighbors. B B B B B Figure 4: An internally conformal subdivision of Std( B ). ¶
18. Correctness of Algorithm
We must prove that our graph G = ( V, E ) is isotopic to thearrangement (
S, T ) in box B . Suppose there are k roots, | S ∩ T | = k . Our correctness requiresthat none of these roots lie in ∂B . Our algorithm produces the following data: we have “wellisolated” the roots in this sense: we have found k aligned boxes, B , . . . , B k such that 2 B i is aroot box, 8 B i ⊆ B , and the interiors of the 8 B i ’s are pairwise disjoint. Next, we have constructedsubdivisions, S , S , . . . , S k where S i is a subdivision of 8 B i ( i = 1 , . . . , k ) and S is a subdivision of B \ ∪ ki =1 B i . Moreover,the union of all these subdivisions, denoted S ∗ , constitutes a balanced box complex of B . Theorem . The PSLG G computed by the algorithm is a S ∗ -normalization of the curves ( S, T ) . We sketch the arguments here: let ( S (cid:48) , T (cid:48) ) be a S ∗ -normalization of ( S, T ). The graph G willbe obtained as the union of G B for all B ∈ S ∗ , where each G B is a PSLG contained in box B . Weknow from Theorem 1 how to construct a PSLG G B ⊆ B that is isotopic to ( S (cid:48) , T (cid:48) ) in each rootbox B . We know from Plantinga-Vegter how to construct PSLG G SB that are isotopic to S (cid:48) in eachnon-root box B . Similarly we have G TB . But we need to form their ”union”, which is the PSLG G B that is isotopic to ( S (cid:48) , T (cid:48) ) in B . For this purpose, we need to know the relative ordering of the10 -vertex and g -vertex on each segment of B that is bichromatic for both curves. This informationis resolved by Stage IX of our construction. We have presented a complete numerical algorithm for the isotopic arrangement of two simplecurves. The underlying paradigm is Domain Subdivision, coupled with box predicates and effectiveforms of the Miranda Test. Moreover, we crucially exploit the previous isotopic approximationalgorithms of Plantinga-Vegter [18] for a single curve.The algorithm is very implementable: despite the many stages, each stage involves iterationusing well-known data structures. A full implementation and comparisons with other methods isplanned; we have currently implemented the root isolation part.The extension of this work to the simple arrangement of multiple curves is of great interest.Many of the techniques we have developed for 2 curves will obviously extend. One possible way touse our work for multiple curves is as follows: first compute the root boxes 2 B i of all the pairwiseintersections, and make them “well isolated” in the sense that 8 B i boxes are pairwise disjoint, asbefore. Then we compute a balanced, conforming subdivision S of complement of the union ofthese 8 B boxes. Moreover, we need to resolve ambiguities, i.e., relative ordering of curves on acommon segment. Some of this can be resolved by propagation, but there will be need for recursivesubdivision in general. In the full paper, we will provide such a description.A general open problem is to prove polynomial complexity bounds for such subdivision algo-rithms. As a first step, we would like to prove that the root isolation part is polynomial-time. Thiswould be a generalization of our recent work on continuous amortization for real and complex roots[21]. References [1] J.-D. Boissonnat, D. Cohen-Steiner, B. Mourrain, G. Rote, and G. Vegter. Meshing of surfaces.In Boissonnat and Teillaud [2]. Chapter 5.[2] J.-D. Boissonnat and M. Teillaud, editors.
Effective Computational Geometry for Curves andSurfaces . Springer, 2006.[3] H. Br¨onnimann, C. Burnikel, and S. Pion. Interval arithmetic yields efficient dynamic filtersfor computational geometry.
Discrete Applied Math. , 109(1-2):25–47, 2001.[4] M. Burr, S. Choi, B. Galehouse, and C. Yap. Complete subdivision algorithms, II: Isotopicmeshing of singular algebraic curves.
J. Symbolic Computation , 47(2):131–152, 2012. SpecialIssue for ISSAC 2008.[5] M. de Berg, M. van Kreveld, M. Overmars, and O. Schwarzkopf.
Computational Geometry:Algorithms and Applications . Springer-Verlag, Berlin, revised 3rd edition edition, 2008.[6] A. Frommer and B. Lang. Existence Tests for Solutions of Nonlinear Equations Using Borsuk’sTheorem.
SIAM J. Numer. Anal. , 43(3):1348–1361, 2005.[7] N. Kamath. Subdivision algorithms for complex root isolation: Empirical comparisons. Mscthesis, Oxford University, Oxford Computing Laboratory, Aug. 2010.118] W. Kulpa. The Poincar´e-Miranda theorem.
The American Mathematical Monthly , 104(6):545–550, Jun–Jul 1997.[9] J.-M. Lien, V. Sharma, G. Vegter, and C. Yap. Isotopic arrangement of simple curves: Anexact numerical approach based on subdivision. In H. Hong and C. Yap, editors,
MathematicalSoftware – ICMS 2014 , volume LNCS 8592, pages 277–282. Springer, 2014. Seoul, Korea, Aug5-9, 2014.[10] L. Lin.
Adaptive Isotopic Approximation of Nonsingular Curves and Surfaces . Ph.D. thesis,New York University, Sept. 2011.[11] L. Lin and C. Yap. Adaptive isotopic approximation of nonsingular curves: the parameteriz-ability and nonlocal isotopy approach.
Discrete and Comp. Geom. , 45(4):760–795, 2011.[12] A. Mantzaflaris, B. Mourrain, and E. P. Tsigaridas. On continued fraction expansion of realroots of polynomial systems, complexity and condition numbers.
Theoretical Computer Science ,412:2312–2330, 2011.[13] R. E. Moore.
Interval Analysis . Prentice Hall, Englewood Cliffs, NJ, 1966.[14] R. E. Moore and J. B. Kioustelidis. A simple test for accuracy of approximate solutions tononlinear (or linear) systems. In
SIAM J. Numer. Anal. [15], pages 521–529.[15] R. E. Moore and J. B. Kioustelidis. A simple test for accuracy of approximate solutions tononlinear (or linear) systems.
SIAM J. Numer. Anal. , 17(4):521–529, 1980.[16] T. S. Newman and H. Yi. A survey of the marching cubes algorithm.
Computers & Graphics ,30:854879, 2006.[17] S. Plantinga.
Certified Algorithms for Implicit Surfaces . Ph.D. thesis, Groningen University,Institute for Mathematics and Computing Science, Groningen , Netherlands, Dec. 2006.[18] S. Plantinga and G. Vegter. Isotopic approximation of implicit curves and surfaces. In
Proc.Eurographics Symposium on Geometry Processing , pages 245–254, New York, 2004. ACMPress.[19] F. P. Preparata and M. I. Shamos.
Computational Geometry . Springer-Verlag, 1985.[20] G. Rote. Extension of geometric filtering techniques to higher-degree parametric curves – curveintersection by the subdivision-supercomposition method. Technical report, Freie Universit¨atBerlin, Institute of Computer Science, 2008. ACS Technical Report No.: ACS-TR-361503-01.[21] M. Sagraloff and C. K. Yap. A simple but exact and efficient algorithm for complex rootisolation. In I. Z. Emiris, editor, , pages 353–360,2011. June 8-11, San Jose, California.[22] V. Stahl.
Interval Methods for Bounding the Range of Polynomials and Solving Systems ofNonlinear Equations . Ph.D. thesis, Johannes Kepler University, Linz, 1995.[23] G. Varadhan, S. Krishnan, Y. J. Kim, S. Diggavi, and D. Manocha. Efficient max-norm dis-tance computation and reliable voxelization. In
Proc. Symp. on Geometry Processing (SGP’03) ,pages 116–126, 2003. 1224] G. Varadhan, S. Krishnan, T. Sriram, and D. Manocha. Topology preserving surface extractionusing adaptive subdivision. In
Proc. Symp. on Geometry Processing (SGP’04) , pages 235–244,2004.[25] C. K. Yap. Complete subdivision algorithms, I: Intersection of Bezier curves. In , pages 217–226, July 2006.[26] C. K. Yap. In praise of numerical computation. In S. Albers, H. Alt, and S. N¨aher, editors,
Efficient Algorithms , volume 5760 of
Lect. Notes in C.S. , pages 308–407. Springer-Verlag, 2009.13 ppendices
Appendix A Basic Concepts
We fix the terminology for well-known concepts in boxes, interval arithmetic and subdivision trees.We define these concepts in d -dimensions. Of course, the algorithms in this paper work in d = 2. ¶
19. Boxes.
Let R denote the set of closed intervals. We may identify R with degenerateintervals [ a, a ] ∈ R . Also R d is the d -fold Cartesian product of R . Elements of R d are called d -boxes . The width of B is ( w ( I ) , . . . , w ( I d )) where the width of an interval I = [ a, b ] is w ( I ) = b − a .the same (resp., differ by at most 1). If B, B (cid:48) are two boxes in R d , we say they are k -neighbors if B ∩ B (cid:48) has dimension k . So k ∈ {− , , , d − } , where the empty set has dimension −
1. Wesay B and B (cid:48) are adjacent if they are ( d − d sides (sometimes called edges ) and 2 d corners . The boundary of a box B is denoted ∂B . ¶
20. Box Functions.
Interval arithmetic [13] is central to our computational toolkit. If f : R d → R is a real function, then we call a function of the form f : R d → R an inclusionfunction for f if for all B ∈ R d , f ( B ) contains f ( B ) = { f ( p ) : p ∈ B } . Call f a box function for f if it is an inclusion function for f and for all { B i ∈ R d : i ∈ N } , if B i converges monotonicallyto a point p ∈ R then f ( B i ) converges monotonically to f ( p ). Note that box functions are easyto construct for polynomials and common real functions. ¶
21. Subdivision Trees.
Our fundamental data structure is a quad-tree or subdivision tree T : the nodes of T are boxes in R d , and each internal node B has 2 d children which are congruentsub-boxes, with pairwise disjoint interiors, and whose union is B . In order to use T to representregions of complex geometry, we assume that each leaf of T is (arbitrarily) either turned ON orturned OFF. The union of all the ON-leaves is denoted R ( T ), called the region-of-interest (ROI).Let On ( T ) denote the set of ON-leaves of T . We call On ( T ) a subdivision of R ( T ). In general,a subdivision of a set X ⊆ R d is a collection C of sets in R d such that ∪ C = X and the relativeinterior of the sets in C are pairwise disjoint. One of the basic operations on subdivision trees isto take an ON-leaf B of T and to “expand it”, i.e., to split B into 2 d congruent sub-boxes andattach them as children of B . Thus B becomes an internal node and its children become leaves ofthe expanded T . By definition, the children of B remain ON-leaves. Thus the ROI is not affectedby expansion.A segment of T is a line segment of the form B ∩ B (cid:48) where B, B (cid:48) are adjacent boxes in T . Notethat a segment is always an edge of some box, but some box edges are not segments. In general,an edge is subdivided into a finite number of segments.The boxes of a subdivision tree are assumed to be non-degenerate, i.e., they are d -dimensional.In our algorithms, certain ON-leaves are called “candidates box”. Unless otherwise noted, we couldassume every ON-leaf is a candidate box. We then say T is balanced if, for any two candidateboxes, if they are adjacent then their depths differ by at most one. Traversing neighbors in a subdivision of ROI:
Given a subdivision tree T partitioningthe ROI, a crucial sub-procedure required by the algorithm is the ability to get the neighbors of aleaf-box in T . One way to achieve this is to associate two pointers with every edge of a leaf boxof T , namely the pointers that point to the extreme neighbors along the edge (there may be onlyone such neighbor, in which the two pointers point to the same box). Thus we associate 8 pointerswith every leaf-box. We will often say the “eight neighbors” of a box to refer to the boxes pointedby these eight pointers, where we count the boxes with multiplicity. We can list all the neighborsof a leaf-box B in T using these eight pointers. 14 tandard Balancing Procedure: Let Q tmp be a priority queue of all the leaves in T ;the deeper the level the higher the priority.While Q tmp is non-empty do B ← Q tmp .pop ().For each neighbor B tmp of B doIf B tmp is not balanced w.r.t. B subdivide B tmp and add its children to Q tmp .There can be at most two neighbors of B that need to be subdivided, because B shares two edgeswith its siblings and so the boxes neighboring B along those edges are balanced w.r.t. B ; theunbalanced boxes can occur on the remaining two edges. Moreover, for any neighbor B tmp that issubdivided only one of its children neighbors B . Balancing also has the following nice property,which intuitively says that the boxes produced in the ensuing subdivision cannot all be very small. Lemma . Suppose we are balancing a box B , and let B (cid:48) be its violating larger neighbor. Let e bethe edge of B (cid:48) shared with B and e (cid:48) be the opposite edge. Then the subdivision of B (cid:48) caused by B while balancing will split the edge e (cid:48) only once. In the subdivision tree of B (cid:48) , the two children that share e (cid:48) are in a different subdivision treecompared to the child of B (cid:48) that is adjacent to B and shares e ; see Figure 5. Balancing produces asubdivision tree of B (cid:48) that has only one path, with leaves hanging from it, that ends in a box whosesize is double the size of B . The number of leaves in this tree are 3 · (log w ( B (cid:48) ) − log w ( B ) − B e e B Figure 5: A subdivision caused balancing.15 ppendix B The Moore-Kioustelidis Test for Roots
Although our paper is focused on arrangement of curves, we shall temporarily consider a moregeneral setting of a continuous function F : R n → R n in n -space. Let the coordinate functions of F be denoted ( f , . . . , f n ). If B = (cid:81) ni =1 I i ⊆ R n is a box, we write B + i and B − i for the pair offaces of B whose outward normal are (respectively) the positive and negative i th semi-axis. Thus,if I i = [ a i , b i ] then B − i = I × · · · × I i − × a i × I i +1 × · · · × I n , and B + i is similar, but with b i inplace of a i . The center of a box B , cen( B ), is defined as the vector (( a + b ) / , , . . . , ( a n + b n ) / λ , define the scaled box λB := { λ ( x − cen( B ) + cen( B )) | x ∈ B } . For X ∈ R , define the magnitude of X , mag( X ) := max x ∈ X | x | .Miranda’s theorem [8] gives us a sufficient condition for the existence of roots of F in the interiorof box B : Proposition . Let F = ( f , . . . , f n ) : R n → R n be a continuous function,and B a box. A sufficient condition that F has a root in the interior of B is that f i ( B + i ) > , f i ( B − i ) < holds for each i = 1 , . . . , n . Remark: we have stated Miranda’s theorem in the simplest possible form. For instance, oursimple form could be generalized by replacing (5) with the following condition: f i takes a definitesign s + i ∈ {− , +1 } on B + i , takes a definite sign s − i on B − i , and s + i s − i = − . But the simplifiedform implies this more general form since we can replace the system F = ( f , . . . , f n ) by (cid:101) F = ( s +1 f , . . . , s + n f n ) , since the systems F and (cid:101) F have exactly the same set of roots. The usual statement of Miranda’stheorem is even general, where (5) is replaced by: there exists a permutation π of the indices { , . . . , n } with this property: for each i , f i has definite signs s + i and s − i on B + π ( i ) and B − π ( i ) (re-spectively), where s + i s − i = − . We shall see that there is no need to find such a permutation, ifwe transform F appropriately. Moore and Kioustelidis [15] give the following effective form of theMiranda test: Proposition . Let F :=( f , . . . , f n ) : R n → R n be a continuous func-tion with appropriate box functions. Write f i,j := ∂f i /∂x j . For any box B with width w ( B ) =( w , . . . , w n ) , if for all i = 1 , . . . , nf i (cen( B + i )) · f i (cen( B − i )) < , (6) | f i (cen( B + i )) | > n (cid:88) j =1 ,j (cid:54) = i mag( f i,j ( B + i )) w j , and (7) | f i (cen( B − i )) | > n (cid:88) j =1 ,j (cid:54) = i mag( f i,j ( B − i )) w j , (8) then F has a zero in the interior of B . roof. Using the mean-value interval extension of f , we know that f i ( B + i ) ⊆ f i (cen( B + i )) + ∇ f i ( B + i ) · ( B + i − cen( B + i ));note the dot-product on the RHS is the inner-product of interval vectors. But ∇ f i ( B + i ) · ( B + i − cen( B + i ))= n (cid:88) j =1 f i,j ( B + i )([ x j , x j ] − ( x j + x j ) / . Since x i = x i , the i th entry in the summation vanishes on the RHS and hence we obtain ∇ f i ( B + i ) · ( B + i − cen( B + i ))= n (cid:88) j =1 ,j (cid:54) = i f i,j ( B + i )([ x j , x j ] − ( x j + x j ) / n (cid:88) j =1 ,j (cid:54) = i f i,j ( B + i ) ( x j − x j )2 [ − , n (cid:88) j =1 ,j (cid:54) = i mag( f i,j ( B + i )) ( x j − x j )2 [ − , n (cid:88) j =1 ,j (cid:54) = i mag( f i,j ( B + i )) ( x j − x j )2 [ − , n (cid:88) j =1 ,j (cid:54) = i mag( f i,j ( B + i ))( w j / [ − , . Thus w ( ∇ f i ( B + i ) · ( B + i − cen( B + i ))) = n (cid:88) j =1 ,j (cid:54) = i mag( f i,j ( B + i )) w j . Therefore, (7) implies that 0 (cid:54)∈ f i ( B + i ). Similarly, (8) implies that 0 (cid:54)∈ f i ( B − i ). By (6), f i takesopposite signs on the faces B + i and B − i , and so Miranda’s theorem implies B contains a root in itsinterior.Miranda’s test is not a “complete” method for detecting roots in the following sense: there aresystems F = 0 whose roots cannot be detected by Miranda’s test, even in the general form thatallows permutation π . For instance, let F = ( f, g ) where f = x + y and g = x − y . Then norectangle B ⊆ R containing the root (0 ,
0) will pass the generalized Miranda test.The solution is a “preconditioning” trick. Consider a transformation of F to G := Y F , where Y is a suitable non-singular matrix in the box B . Note that G and F have the same sets of roots. Toperform the Miranda Test on a box B , we choose Y to be the inverse of any non-singular Jacobian J F ( m ) where m ∈ B . More precisely, MK-test for a system F on a box B is the effective Miranda-test ap-plied to the system J F ( m ) − F , where m := cen( B ), and the Jacobian is non-singular. (9)17his idea was first mentioned by Kioustelidis and its completeness was shown by Moore-Kioustelidis [15]. We reproduce their result, but to do that we need some notation and the MeanValue Theorem in higher dimensions.Given x, y ∈ R , the notation x ± y denotes a number of the form x + θy , where θ is such that0 ≤ | θ | ≤
1; thus “ ± ” hides the θ implicit in the definition. We further extend this notation tomatrices in the following sense: for two matrices A, B , the matrix A ± B :=[ a ij ± b ij ]; also, for ascalar λ , the matrix A ± λ :=[ a ij ± λ ]. We now recall the Mean Value Theorem for F : R n → R n :Given two points x , y ∈ R n , there exists a matrix K with non-negative entries such that F ( x ) − F ( y ) = ( J F ( y ) ± K (cid:107) x − y (cid:107) ) · ( x − y ) . (10)To see this claim, we apply the mean value theorem twice in each of the components of F to obtain f i ( x ) − f i ( y )= ( f i, ( y ) ± K i, (cid:107) x − y (cid:107) , · · · , f i,n ( y ) ± K i,n (cid:107) x − y (cid:107) ) · ( x − y )= ∇ f i ( y ) · ( x − y ) ± ( K i, , . . . , K i,n ) · ( x − y ) (cid:107) x − y (cid:107) for i = 1 , . . . , n . Lemma . Let F be a zero-dimensional system of polynomials. For all sufficiently small open boxes B containing a single root α of F the modified system G := J F ( m ( X )) − F , if well defined, satisfiesthe conditions in Miranda’s theorem, namely for i = 1 , . . . , n , g i ( B + i ) ≥ and g i ( B − i ) ≤ .Proof. Let x be a point on the boundary of the box B . From the definition of G and from themean value theorem (10) we know that G ( x ) = J F ( m ) − ( F ( α ) + ( J F ( m ) ± K (cid:107) x − α (cid:107) ) · ( x − α ))= J F ( m ) − ( J F ( m ) + K (cid:107) x − α (cid:107) ) · ( x − α ))= (1 ± (cid:107) J F ( m ) − K (cid:107) ∞ (cid:107) x − α (cid:107) ) · ( x − α ) . The i th component in the vector(1 ± (cid:107) J F ( m ) − K (cid:107) ∞ (cid:107) x − α (cid:107) ) · ( x − α ) (11)is the polynomial g i ( B ), so we obtain | g i ( x ) − ( x i − α i ) | ≤ (cid:107) x − α (cid:107)(cid:107) J F ( m ) − K (cid:107) ∞ n (cid:88) j =1 | x j − α j | . (12)The term on the RHS (cid:107) x − α (cid:107)(cid:107) J F ( m ) − K (cid:107) ∞ n (cid:88) j =1 | x j − α j |≤ (cid:107) ˆ w ( B ) (cid:107) (cid:107) J F ( m ) − K (cid:107) ∞ , because (cid:107) x − α (cid:107) ≤ (cid:107) ˆ w ( B ) (cid:107) ≤ (cid:107) ˆ w ( B ) (cid:107) and (cid:80) nj =1 | x j − α j | ≤ (cid:107) ˆ w ( B ) (cid:107) . Suppose the box B is suchthat 2 (cid:107) ˆ w ( B ) (cid:107) (cid:107) J F ( m ) − K (cid:107) ∞ < min i =1 ,...,n (cid:107) α − B ± i (cid:107) then we claim that for all i = 1 , . . . , n , g i ( B + i ) ≥ g i ( B − i ) ≤
0. This is because for all x ∈ B + i , | x i − α i | = | x i − α i | = (cid:107) α − B + i (cid:107) , since the projection of α on B − i is ( α , . . . , α i − , x i , α i +1 , . . . , α n );similar argument applies for x ∈ B − i . Thus the term on the RHS in (12) is smaller than | x i − α i | / g i ( B + i ) ≥ g i ( B − i ) ≤ G ( x ) has the same sign pattern as x − α on the boundary of the box B .18his “orthogonalization” around the zero by the pre-conditioning step helps us avoid findingthe permutation matrix in the general Miranda’s test. Note, however, that if the root is on theboundary of the box then the above proof breaks down.19 ppendix C Proofs and Details Proof of Theorem 1:
We will need the following lemma for the proof.
Lemma . If a box B satisfies MK( B ) and an f -vertex and a g -vertex share an edge e of B thenwe can determine the relative order of the normalized curves ( S (cid:48) , T (cid:48) ) along e .Proof. Since the MK( B ) test passed along e , we know that there are real numbers a, b such thateither a · f ( e ) > b · g ( e ) or a · f ( e ) < b · g ( e ). To see this, recall that MK( B ) test replaces the system F = ( f, g ) T by the system (cid:98) F = J · F , where J is the inverse of the Jacobian of F evaluated atcen( B ), and performs the Miranda test, Proposition 5, for (cid:98) F . If J = (cid:20) a − bc d (cid:21) and (cid:98) F = ( (cid:98) f , (cid:98) g ) T then (cid:98) f = a · f − b · g . The Miranda test on (cid:98) F asserts that there is an edge e for which either (cid:98) f ( e ) > (cid:98) f ( e ) <
0. The first inequality is equivalent to a · f ( e ) > b · g ( e ), and the second inequality isequivalent to a · f ( e ) < b · g ( e ). In the rest of the proof we assume that a · f ( e ) > b · g ( e ); theanalysis in the other case is same.Neither a nor b can vanish, since that would imply that either f or g has a constant sign on e , which is a contradiction as both f and g have a vertex on e . Let e ( t ) be a parametrization of e with endpoints e (0) and e (1). Let T f ⊆ (0 ,
1) be such that f ( e ( t )) = 0 for all t ∈ T f , and let t f be the smallest element in T f ; similarly define T g and t g . Since both f and g change sign across e , we know that the cardinality of T f and T g is odd. Any normalization ( S (cid:48) , T (cid:48) ) of ( S, T ) relativeto B will remove all but one element from both T f and T g , while maintaining the relative order ofthe remaining element. That order is the same as the order of t f and t g along e . Thus we want todetermine whether t f < t g or t g < t f . Suppose ab >
0. Then f ( e ) > c · g ( e ) for some c >
0. Thereare two cases to consider: • f ( e (0)) >
0: then f ( e ( t g )) > cg ( e ( t g )) = 0, which implies that f is positive at e ([0 , t g ]) andso t f > t g ; • f ( e (0)) <
0: this similarly implies t f < t g .If ab < g ( e ) > c · f ( e ), for some c >
0, and the claim follows from similar arguments. ¶
22. Group I Patterns.
Notice that using the sign of f, g at the corners of B , we can neverdetect these patterns. For instance, for Figure 2(Ia), we will not detect the presence of the curve S (cid:48) because f has the same sign on every corner of the box. So we first show that they cannot arise. Lemma . Suppose box B satisfies MK( B ) . Then the patterns in Group I of Figure 2 cannot occur.Proof. Let e be an edge of B and suppose S (cid:48) ∪ T (cid:48) intersect e in three consecutive points e ( t ) , e ( t ) , e ( t )( t < t < t ) where e ( t ) is a parametrization of e . The “pattern” of these intersections is the triple( p , p , p ) where p i ∈ { f, g } . For instance, if e is the top edge of the box in Figure 2(Ia), thenthe pattern is either ( f, g, f ) or ( g, f, g ). Our claim is equivalent to showing that the intersectionpattern of any three consecutive intersections of S (cid:48) ∪ T (cid:48) on any edge e of B cannot be ( f, g, f ) or( g, f, g ).From Lemma 7 we know that f ( e ) > c · g ( e ), for some c ∈ R (cid:54) =0 ; let us assume c >
0. Considerthe ( f, g, f ) pattern (the other pattern is similar). Consider the sign of g at the point e ( t − ε ) and e ( t + ε ) for sufficiently small ε >
0. Then g must have different signs at these points — this isbecause as we move from e ( t − ε ) to e ( t + ε ), the function g changes sign exactly once, at e ( t ).Likewise, we see that f must have the same sign at e ( t − ε ) and e ( t − ε ), because as we movefrom e ( t − ε ) to e ( t + ε ), the function f changes sign exactly twice, at e ( t ) and e ( t ). Thus f ( e ( t − ε )) > g ( e ( t − ε )) iff f ( e ( t − ε )) < g ( e ( t − ε )). This is a contradiction.20
23. Group II Patterns.
Suppose f, g have sign agreement on B . We can determine fromthese signs the two edges that contains f - and g -vertices. Suppose e is such an edge. So there isan f -vertex and a g -vertex on e , and from Lemma 7 we know their relative ordering. ¶
24. Group III Patterns.
Let us say that f, g have sign agreement on B if there is a sign s ∈ { +1 , − } such that sign ( f ( c ) g ( c )) = s for each corner c of B . Observe that Group II patternsarise precisely because f, g have sign agreement; likewise Group III patterns arise precisely because f, g do not have sign agreement. We claim that the patterns in Group III can be determined bysigns of f and g at the corners of B . First of all, by evaluating the signs of f and g on the cornersof B , we can determine whether or not f, g have sign agreement of B . If not then we can determinewhether the pattern is (IIIa), (IIIb) or (IIIc). If (IIIa), the pattern is completely determined. If(IIIb), there is an edge e containing both an f - and a g -vertex, and we need to know their relativeorder on e . This is determined by the positions of the other f -vertex and other g -vertex: this isbecause the order of the four f - and g -vertices on the boundary of B must be alternating: f, g, f, g .A similar remark applies in case (IIIc).To summarize the proof of Theorem 1: Lemma 8 implies that Group I patterns cannot occur;for Group II patterns we can determine the relative order from Lemma 7 and for Group III patternsthe ordering is immediate. ¶
25. The RefineRoot Procedure:
RefineRoot( B ) (cid:47) Assume that
JC(6 B ) holds. (cid:47) Thus no neighbor of B can be an MK-box. Input: an aligned box B with 2 B as the root box.Output: an aligned box B ∗ ⊂ B with 2 B ∗ as the root box. Algorithm (cid:46)
Remove B from Q MK .Subdivide the neighbors of B until the size of theneighborhood of B is w ( B ) / Q , Q f , Q g , Q fg .Initialize Q tmp with the neighbors of B and its children.While Q tmp is non-empty do B tmp ← Q tmp .pop ().If MK(2 B tmp ) holds thenEmpty Q tmp into Q fg .Return B tmp and add it to Q MK .Else Subdivide B tmp and add its children to Q , Q f , Q g and Q tmp if they satisfyrespective predicates..Correctness: The subdivision of B and its neighborhood of size w ( B ) / B , the root boxcorresponding to B . Let B (cid:48) be any of these 16 boxes. Since JC(6 B ) holds, if MK(2 B (cid:48) ) holds for abox B (cid:48) then the root in 2 B (cid:48) is exactly the root in 2 B .We now give the details of various stages mentioned in § ¶
26. Details of Stage III: Q JC is non-empty B ← Q JC .pop (). Q tmp ← { B } While Q tmp is non-empty do B tmp ← Q tmp .pop ().If MK(2 B tmp ) holds.Push B tmp into Q MK . Empty Q tmp into Q fg .Else subdivide B tmp and distribute the childreninto Q , Q f , Q g , Q tmp (after testing for the corresponding predicates).For each box B ∈ Q MK doIf there is another box B (cid:48) ∈ Q MK such that2 B ∩ B (cid:48) (cid:54) = ∅ then remove B (cid:48) from Q MK .Note that we only search for a root in f g -candidate boxes. This is justified by Lemma 6 andthe observation that eventually the root will be contained in the interior of the doubling of an f g -candidate box. At the end, Q JC is empty and Q MK contains a set of root boxes. Moreover, thelast loop ensures no two boxes B, B (cid:48) ∈ Q MK correspond to the same root, i.e., 2 B ∩ B (cid:48) = ∅ . Theboxes in Q fg do not contain any root. ¶
27. Details of Stage V:
For each box B in Q Root do the following steps.Initialize Q tmp with all the neighbors of B in T .While Q tmp is non-empty do B tmp ← Q tmp .pop ().If B tmp ⊂ B then turn it OFF andadd its neighbors to Q tmp .If the interior of B tmp intersects the interior of 8 B thensubdivide it and add its children to Q tmp . (cid:47) NOTE: Whenever we subdivide a box B tmp (cid:47) we remove it from one of the queues Q f , Q g , or Q fg (cid:47) and add its children to the appropriate queue. Since 8 B is half-aligned, there is a refinement of T such that every box in this refinement iseither contained in 8 B or does not intersect its interior. Thus the procedure described above willterminate. Let T (cid:48) be the refinement of T with blacked-out regions corresponding to extended rootboxes. ¶
28. Details of Stage VI: B ∈ Q Root doLet m be the largest depth amongst all theneighbors B tmp of 8 B in T (cid:48) .Let (cid:96) be the depth of B in the subdivision tree T . (cid:47) Thus w ( B tmp ) = w ( B )2 (cid:96) − m If m > (cid:96) then k ← m ; else k ← (cid:96) + 1.Add a conceptual leaf box to T (cid:48) that represents 8 B .Set the depth of this box to k + 1 and initializeits 8 pointers to the 8 neighbors of 8 B in T (cid:48) .Let T (cid:48)(cid:48) be the resulting subdivision tree.Let Q be the priority queue of all the leaves in T (cid:48)(cid:48) ;the deeper the level the higher the priority.Initiate the standard balancing procedureon Q with one difference: whenever we pop a conceptualbox 8 B we check the depth of its neighbors and ifnecessary reset the depth of 8 B to one more thanthe depth of its deepest neighbor. (cid:47) NOTE: Whenever we subdivide a box B tmp we (cid:47) remove it from one of the queues Q f , Q g , or Q fg (cid:47) and add its children to the appropriate queue. We claim that at the end of this procedure the tree T (cid:48) is balanced, and all the neighbors ofextended root boxes 8 B i in T (cid:48) are of the same size, namely w ( B i ) / k , for some k ≥
1. Thebalancing of B \ ∪ i (8 B i ) follows from the proof of correctness for standard balancing procedure.The conformity follows because a conceptual box is always deeper in T (cid:48)(cid:48) than its neighbors, soit will never be subdivided, and its neighbors will always be twice its size. The modification tothe standard balancing is required, because a smallest neighbor B tmp of 8 B in T (cid:48) could have beensubdivided by a box that is adjacent to B tmp along the edge that is not abutting 8 B or any of theneighbors of 8 B . However, this can only happen once because of the balancing property, Lemma 3. ¶
29. Details of Stage IX: Q tmp with all the root boxes. (cid:47) Q tmp will contain unambiguous boxes. For each box B ∈ Q fg doIf there is pair of f -vertex and g -vertex that do not sharea segment of B thenConnect the two f -vertices with an edge;connect the two g -vertices with an edge;ensure that the two edges do not intersect.Add B to Q tmp and remove it from Q fg . (cid:47) In the remaining boxes, the two pairs of (cid:47) ( f, g ) -vertices share the same segments. If the two pairs of ( f, g )-vertices are on edges e , e (cid:48) that share a vertex then (cid:47) Call such a box a
Transition Box (cid:47)
These boxes definitely appear in a covering (cid:47) of nested f g -loops; they can appear otherwise.
Subdivide both e and e (cid:48) until we reacha subset e (cid:48)(cid:48) in one of the edges such thatonly one of the curves f or g changes sign on e (cid:48)(cid:48) ;say e (cid:48)(cid:48) ⊂ e and f changes sign on it.Check which side of e \ e (cid:48)(cid:48) does g change sign;order the f -vertex and g -vertexalong e accordingly; connect the f -verticesand g -vertices respecting this order;add B to Q tmp and remove it from Q fg .If B shares an edge e with B thenSubdivide e until we reach a subset e (cid:48)(cid:48) ⊂ e such that only one of the curves f , g changessign on e (cid:48)(cid:48) . Check which side of e \ e (cid:48)(cid:48) containsthe other curve. Order the vertices accordinglyand connect the f -vertices and g -vertices.Add B to Q tmp and remove it from Q fg . (cid:47) The boxes in Q tmp are all unambiguous boxes. While Q tmp is non-empty do B ← Q tmp .pop ()For each ambiguous f g -candidate B tmp of B doOrder the f -vertices and g -vertices on theshared segment between B tmp and B accordingto their ordering in B ; connect the pair of f -verticesand g -vertices in B tmp respecting this ordering.Add B tmp to Q tmp and remove it from Q fg . (cid:47) Thus all the f g -neighbors are unambiguous.
In practice, we should first resolve boxes that can be traced to root boxes. Then we should resolvetransition boxes and propagate their ordering. Finally, in the remaining ambiguous boxes, weshould resolve the boundary boxes and propagate their ordering. At the end of this stage Q fgfg