Voronoi diagrams of semi-algebraic sets by Frangois Anton M . A . T . D . R , Universite Laval, 1990 Diplome de l'lnstitut de Topometrie du Conservatoire National des Arts et Metiers, France, 1992 M . Sc., Universite Laval, 1995 A THESIS S U B M I T T E D I N P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F D o c t o r of P h i l o s o p h y in T H E F A C U L T Y O F G R A D U A T E STUDIES (Department of Computer Science) We accept this thesis as conforming to the required standard The University of British Columbia January 2004 © Frangois Anton, 2004 Abstract Most of the curves and surfaces encountered in geometric modelling are denned as the set of solutions of a system of algebraic equations and inequalities (semi-algebraic sets). Many problems from different fields involve proximity queries like finding the (nearest) neighbours or quantifying the neighbourliness of two objects. The Voronoi diagram of a set of sites is a decomposition of space into prox-imal regions. The proximal region of a site is the locus of points closer to that site than to any other one. Voronoi diagrams allow one to answer proximity queries after locating a query point in the Voronoi zone it belongs to. The dual graph of the Voronoi diagram is called the Delaunay graph. Only approximations by conies can guarantee a proper order of continuity at contact points, which is necessary for guaranteeing the exactness of the Delaunay graph. The theoretical purpose of this thesis is to elucidate the basic algebraic and geometric properties of the offset to an algebraic curve and to reduce the semi-algebraic computation of the Delaunay graph to eigenvalues computations. The practical objective of this thesis is the certified computation of the Delaunay graph for low degree semi-algebraic sets embedded in the Euclidean plane. The methodology combines interval analysis and computational algebraic geometry. The central idea of this thesis is that a (one time) symbolic preprocessing may accelerate the certified numerical evaluation of the Delaunay graph conflict locator. The symbolic preprocessing is the computation of the implicit equation of the generalised offset to conies. The reduction of the Delaunay graph conflict loc-ator for conies from a semi-algebraic problem to a linear algebra problem has been possible through the use of the generalised Voronoi vertex (a concept introduced in this thesis). The certified numerical computation of the Delaunay graph has been possible ii by using an interval analysis based library for solving zero-dimensional systems of equations and inequalities (ALIAS). The certified computation of the Delaunay graph relies on theorems on the uniqueness of a root in given intervals (Kantorovitch, Moore-Krawczyk). For conies, the computations get much faster by considering only the implicit equations of the generalised offsets. iii Contents Abstract ii Contents iv List of Tables viii List of Figures ix List of Symbols, Nomenclature and Abbreviations xiii Acknowledgements xvii Dedication xix 1 Introduction 1 1.1 The problem 10 1.2 The motivation . . . ; 13 1.3 Outline of the research 17 2 Generalisations of the Voronoi diagram for curved objects 24 2.1 Voronoi diagrams of general manifolds 24 2.2 Voronoi diagrams for curved objects 28 2.3 The Voronoi diagram and medial axis transform for planar domains with curved boundaries 29 3 Combining symbolic computation and scientific computation 38 iv 3.1 The exact symbolic Delaunay graph conflict locator for circles . . . . 40 3.1.1 Preliminaries 40 3.1.2 The Delaunay graph conflict locator for additively weighted points 42 3.1.3 The Delaunay graph conflict locator for circles 53 3.2 Formulating the Delaunay graph conflict locator for curves from those curves 62 3.3 A n hybrid approach linking symbolic computation and scientific com-putation . 66 3.3.1 The sparse resultant 67 3.3.2 Numerical methods for computing exactly the signs of the eigenvalues of large sparse matrices 74 3.3.3 A L I A S 76 4 The offset to an algebraic curve 81 4.1 Equations defining the offsets 82 4.2 The degree of the generalised offset curve 89 4.3 The degree of the generalised offset to a conic 102 4.4 A n implicit equation of the generalised offset to a conic 105 4.5 Conclusion's 118 5 The Delaunay graph conflict locator for conies 120 5.1 Preliminaries 120 5.2 Formalisation of the Delaunay graph conflict locator 121 5.3 The simplification of the Delaunay graph conflict locator 126 5.4 The algebraic precomputations 131 5.5 The numerical computation of the Delaunay graph conflict locator . 132 6 The Delaunay graph conflict locator for semi-algebraic sets 135 v 6.1 The algebraic equations and inequalities of the Delaunay graph con-flict locator 136 6.2 The formulation of the Delaunay graph conflict locator with ALIAS 138 6.3 The hybrid symbolic/scientific computation of the Delaunay graph conflict locator for conies 142 7 Conclusions 146 7.1 Limitations of this research 149 7.2 Future research 150 Bibliography 152 Appendix A The Macaulay 2 program for the exact Delaunay graph conflict locator for the Additively Weighted Voronoi diagram 164 Appendix B A n implicit equation of the generalised offset to a conic 167 B. l The Maple program for calling Resin- for obtaining the matrix of the ' sparse resultant 167 B.2 The Maxima program for computing the sparse resultant 172 B.3 An implicit equation for parabolas 173 B. 4 An implicit equation for ellipses, circles or hyperbolas 174 Appendix C The Singular program for obtaining the matrix of the sparse resultant of the Delaunay graph conflict locator for conies 175 C l The case of parabolas 175 C. 2 The case of ellipses and/or hyperbolas 177 Appendix D A Maple program for the ALIAS based Delaunay graph conflict locator for semi-algebraic sets 181 D. l The system without the implicit equations of the generalised offsets 182 D.2 The system with the implicit equations of the generalised offsets . . 183 vi Index 187 v i i List of Tables 3.1 The comparison of different algebraic methods for computing the Delaunay graph conflict locator for a circle, a parabola, a hyperbola and a circle 66 4.1 The degree of the generalised offset to conies 105 5.1 The mixed volumes and sparse resultant degrees of the different con-figurations of conies 132 6.1 Some running time results with different ALIAS parameters on the system with the generalised offsets 140 6.2 Some running time results for ellipses with the equations of the ori-ginal curves 140 6.3 Some running time results for ellipses with the equations of the ori-ginal curves for the partial system 141 6.4 Some running time results for ellipses with the equations of the gen-eralised offsets 142 6.5 The ratios of the running time results without/with the equations of the generalised offsets 142 6.6 Some running time results for hyperbolas with the original equations of the hyperbolas 144 6.7 Some running time results for hyperbolas with the equations of the generalised offsets 144 viii List of Figures 1.0.1 The ordinary Voronoi diagram of points, and its topology expressed by the Delaunay triangulation 4 1.0.2 The bisector of a point and a line segment 5 1.0.3 The influence zone of a point with respect to a line 5 1.0.4 The Voronoi zone of a cubic 5 1.0.5 The Voronoi diagram of a circle, an ellipse and a hyperbola . . . . 6 1.0.6 The Delaunay graph of the sites of Figure 1.0.5 '. . 7 1.0.7 The 3 empty circles for the sites of Figure 1.0.5 8 1.2.1 The relative position with respect to the bisector must be constant 18 2.1.1 The space of generalised spheres 26 2.1.2 The circles tangent to a given circle in the space of generalised spheres 27 2.3.1 The difference between (a) the true Voronoi diagram of a planar domain bounded by curves, and the computed from piecewise-linear boundary approximations with (b) 15 segments and (c) 60 segments 31 2.3.2 The internal part of the Voronoi diagram and the medial axis are identical for the exact boundary while they differ for the approxim-ate boundary 32 2.3.3 The bisector of two curves as the envelope of a family of point/curve bisectors 33 2.3.4 The left and right candidate points 36 3.1.1 The Additively Weighted Voronoi diagram . 41 ix 3.1.2 The Apollonius Tenth problem 42 3.1.3 The starting configuration 45 3.1.4 The real configuration after addition of the fourth weighted point . 46 3.1.5 The configuration computed by an approximate algorithm 47 3.1.6 The Delaunay graph conflict locator for the Additively Weighted Voronoi diagram 48 3.1.7 The Additively Weighted Voronoi vertex as the common intersection of three expanding circles 49 3.1.8 The Additively Weighted Voronoi vertex as the common intersection of three shrinking circles 49 3.1.9 There is no such triangle in the old Delaunay graph because of the absence of a real Voronoi vertex 52 3.1.10 Two triangles can possibly disappear simultaneously by the addition of a single weighted point 53 3.1.11 Seven Apollonius circles centres that are true Voronoi vertices . . . 54 3.1.12 Four Apollonius circles centres that are true Voronoi vertices . . . 57 3.1.13 Two Apollonius circles centres that are true Voronoi vertices (first case) 58 3.1.14 Two Apollonius circles centres that are true Voronoi vertices (second case) 58 3.2.1 The Delaunay graph conflict locator that certifies whether the ad-dition of a curve C 4 changes or does not change each one of the triangles induced by three curves C\, C2, and C 3 63 3.3.1 The Minkowski sum of the point sets A and B 69 3.3.2 A mixed subdivision of the Minkowski sum shown in Figure 3.3.1 and the corresponding mixed volume 70 4.1.1 The strophoid (C) and its true offset (thick lines) 83 4.1.2 The strophoid and its generalised offset 84 x 4.1.3 Relationships between the generalised offset to the strophoid and the strophoid 85 4.1.4 The construction of a point of the generalised offset . 88 4.4.1 The Newton polytope of / 107 4.4.2 The Newton polytopes of n and of an — (3f 108 4.4.3 The Newton polytopes of d and of / — ad 108 4.4.4 The mixed volume of / and an — (3f 109 4.4.5 The mixed volume of an — f3f and / — ad 109 4.4.6 The mixed volume of / and / — ad 109 4.4.7 The Newton polytope of / for ellipses and hyperbolas I l l 4.4.8 The Newton polytope of n for ellipses and hyperbolas I l l 4.4.9 The Newton polytopes of d and of ad — f for ellipses and hyperbolas. 112 4.4.10 The mixed volume of / and n for ellipses and hyperbolas 112 4.4.11 The mixed volume of n and / — d for ellipses and hyperbolas. . . . 113 4.4.12 The mixed volume of / and / — d for ellipses and hyperbolas. . . . 113 4.4.13 A n ellipse and its 0.5—generalised offset 1.14 4.4.14 The same ellipse and its 3—generalised offset 115 4.4.15 A n hyperbola and its 0.5—generalised offset 115 4.4.16 The same hyperbola and its 3—generalised offset 115 4.4.17 The Newton polytope of / for parabolas . 116 4.4.18 The Newton polytope of n for parabolas 116 4.4.19 The Newton polytopes of d and of / — d for parabolas 117 4.4.20 The mixed volume of / and n for parabolas 117 4.4.21 The mixed volume of n and / — d for parabolas 118 4.4.22 The mixed volume of / and / — ci for parabolas 118 4.4.23 A parabola and its 0.5—generalised offset 119 4.4.24 The same parabola and its 3—generalised offset 119 5.2.1 A generalised Voronoi vertex of three conies 122 xi 5.2.2 A true Voronoi vertex of three conies 122 5.2.3 The insertion of C4 induces a conflict with the triangle C1C2C3 . . 123 5.2.4 The new Delaunay graph after insertion of C 4 124 5.2.5 The Delaunay graph conflict locator for conies . 124 5.3.1 The Newton polytope of the implicit equation of the r—generalised offset to a parabola in a system centred on its summit and with x-axis the axis of symmetry of the parabola 127 5.3.2 The Newton polytope of the implicit equation of the r—generalised offset to a parabola in an arbitrary system 128 5.3.3 The Newton polytope of the difference of the implicit equations of the generalised offset to a parabola in a nice system and in an arbitrary system 129 5.3.4 The Newton polytope of the implicit equation of the r—generalised offset to an ellipse or an hyperbola in a system centred on the centre of the conic and with axes the axes of symmetry of the conic . . . 129 5.3.5 The Newton polytope of the implicit equation of the r—generalised offset to an ellipse or an hyperbola in an arbitrary system 130 5.3.6 The Newton polytope of a linear combination of the implicit equa-tions of the generalised offset to an ellipse or an hyperbola in a nice system and in an arbitrary system 130 5.5.1 The test case with four ellipses 133 6.3.1 The test case with four hyperbolas 143 xi i List of Symbols, Nomenclature and Abbreviations A + B, 68 R, 62 Au 68 i? = P ! + ... + P m , 69 B = V(f)c K4, 90 ^es/,,,...,^, 71 B (si,Sj), 3 5 C R*, 69 C", 46 Us, 25 Cl, 45 ^(5), 3 C = V(f), 85 F ( / i , . . . , / s ) c £ , 83 C 5 , 99 ^ ( s i , 5 ) , 3 . CQO, 98 Vh 136 ' Ci, 40, 62 f o Z ( Q ) , 68 D = V(d) CKA, 90 VolN (Q), 68 ZX7(«S), 3 91 D(si,Sj), 3 Ws,91 F ( ^ ) , 79 . . Woo, 91 F i , 79 Wa, 91 7 = c { n C £ n C ^ , 4 6 X i , 136 K[xf1,...,x±1]=K[x, x " 1 ] , 68 Z , 25 K[xi,...,xN], 68, 82 V, 86 K [Xl, . . . ,X m ]q j , 98 i z , 25 MV(QI,...,QJV) , 69 M, 86 MI/_ i = M y (Q0,...,Qi. M , 137 72 P N , 71 M s , 51 6 , 86 P i , 40 5, 3 Q, 68 © i , 137 Q = Q i + ... + Qm C M 7 v , 69 5, 3 xii i S(M,Z), 25 4 , 4 7 ^ 9 3 Ci(xi,yi),62 (h, -Js), 93 d ( M ) C i ) , 5 4 (M,d) , 3 d(M,Pi),40 fez), 7 7 <* 3 (ari,2/i),44 d ( s ,y ,u ,v ) , 87 (a*, 2/0. 6 2 di, 137 C = {CI,...,CJV},'53 di (xi,yi,x,y,r), 62 •F<(a:i,...,:rm) {=,>}(), 79 / f e y ) ; 85 2:, 79 / r 9 i M -V {n) c K4, 90 m = («,/?), 86 136 mg, 51 V ( C i , C ) , 54 n(s,2/,u,u), 87 V(Pi),V,41 n i M 136 V(C) , 54 mix^yux^), 62 V(7>),41 p = f e y ) , 8 7 *4> 7 2 Pi , 40 93,98 g = (u,v),85 / ^ ( 7 ( F ) , K ( G ) ) , 9 8 q = (yi,...,yN), 136 r* P(/,ff), 98 r, 46 0 , 87 r i ; 53 V 7, 91 ^ , 40 7,91 a;, 2/, 46 Tr, 87 x a , 68 v/> 9 3 79 1 3n — (-^il i -^ i;v )> 136 a = (au...,aN) e ZN, 68 sf, 79 xiv abstract Voronoi diagrams, 11 ideal, 49 additively weighted Voronoi diagram, inclusion monotonic, 77 41 influence zone, 3 additively weighted Voronoi vertex,/45 intersection multiplicity, 98 affine variety, 83 interval arithmetic rules, 77 algebraic variety, 1 interval function, 77 algebraically closed field, 82 interval number, 77 A L I A S / 3 B , 80 irreducible, 92 A L I A S / M a x 3 B , 80 localisation, 98 ALIAS/mixed_bisection, 78 ALIAS/single.bisection, 78 manifold, 25 metric, 3 Bezout number, 70 metric space, 3 bisector, 3 Minkowski sum, 68 boxes, 79 mixed bisection, 78 centre of an ellipse or hyperbola, 110 mixed cell, 69 component at infinity, 90 mixed volume, 69 multiplication operator, 51 degree, 83 Delaunay graph, 3 polyhedral subdivision, 69 Delaunay graph conflict locator, 8 prime ideal, 98 Delaunay graph predicate, 8 projective completion, 90 dimension, 50 projective, hypersurface, 92 projective space, 71 full bisection, 78 projective variety, 83 generalised offset, 84 quasi-projective variety, 71 Grobner basis, 50 quotient algebra, 50 homogenisation, 91 rational function, 98 XV reducible, 92 ring of polynomials, 50 semi-algebraic set, 1 single bisection, 78 singular component, 90 singular point, 86 sites, 3 sparse resultant, 71 square-free, 90 summit of a parabola, 110 support of a polynomial, 68 tangent space, 85 to touch, 13 toric variety, 71 true offset, 83 Voronoi diagram, 3 Voronoi region, 3 weight circle, 40 weighted point, 40 Zarisi closed set, 71 Zariski closure, 71 Zariski topology, 71 Acknowledgements I would like to express my most profound gratitude to my supervisor Dr. David Kirkpatrick for accepting to supervise me in a research field that was far from his research interests, for his continuous support and availability, and for his invaluable help with most useful comments, questions and suggestions as well as constant in-terest in my research. I would never have been able to achieve my dream of working on Voronoi diagrams for curves without his invaluable supervision. I would like to show my deepest gratitude to my co-supervisor Dr. James Carrell, and the other members of my thesis committee Drs. Anne Condon and Joel Friedman for their interest in my research and their support, interesting questions, comments and sug-gestions. I would like to acknowledge and express my deep gratitude for the support of Dr. Jean-Pierre Merlet during and after my visit at INRIA (Institut National de Recherche en Informatique et en Automatique) Sophia Antipolis. He introduced me to the A L I A S library for analysing and solving zero-dimensional systems of equa-tions and inequalities and to its technicalities. I would like to express my gratitude to Drs. Ioannis Emiris, Monique Teillaud, and Bernard Mourrain for their comments and corrections, and their help both during and after my visit at INRIA Sophia A n -tipolis. I would like to express my gratitude to Dr. Jean-Daniel Boissonnat for his comments and interest in my research during my visit at INRIA Sophia Antipolis Computational Geometry group ( P R I S M E project). I would like to thank the support staff from the "Unite Mixte de Service" CNRS/Ecole Polytechnique Medicis (France) and the Department of Computer Sci-ence at U B C for their kind technical support. I am particularly grateful to Dave Brent, Michael Sanderson, Teresa Gomez-Diaz, Chris Majewski, Eric Schost, and Pierre Lafon. I would like to thank my fellow graduate students at the Theory group and Beta Lab at U B C and at the P R I S M E and G A L A A D projects at INRIA Sophia Antipolis for their interest for and support to my research. I am particularly grate-xvii ful to Maja Dimitrijevic, Sanja Rogic, Julia Flototto, Philippe Guigue, Frangois Rebufat, David Cohen-Steiner, Drs. Ellen Gethner and Alex Brodsky, Stephane Durocher, Dan Tulpan, Mirela Andronescu and Mark McCann. M y Ph.D. thesis research has been made possible by a Natural Sciences and Engineering Research Council of Canada (NSERC) Post Graduate Student-ship (PGS B) , a British Columbia Advanced Systems Institute (BC ASI) Graduate Scholarship, a "bourse doctorale" of the Fonds Chercheurs et Aide a la Recherche (FCAR) , a "Bourse du Gouvernement Frangais" (BGF) , a scholarship from the In-stitut National de Recherche en Informatique et en Automatique (INRIA), teaching assistant contracts from the Department of Computer Science at U B C and research assistant contracts from my supervisor Dr. David Kirkpatrick. Finally, I would like to express my strong gratitude to my wife Dr. Darka Mioc for accepting my visit at INRIA Sophia Antipolis as well as my long research and endless computations. FRANCOIS ANTON The University of British Columbia January 2004 xviii A Marie-Christine et a Darka. A mes parents Roger Anton et Francisca Castro. Chapter 1 Introduction We can encounter almost everywhere in the real world curved objects in a three dimensional Euclidean space. Most of the curves and surfaces encountered in geo-metric modelling are defined as the set of common zeroes of a set of polynomials (algebraic varieties) or subsets of algebraic varieties defined by one or more algeb-raic inequalities (semi-algebraic sets). More formally, let us recall the following definition. Def in i t ion 1.0.1. (Semi-algebraic set, adapted from [BCR98, Definition 2.1.4]) A semi-algebraic set of R ^ is a subset of the form \Jf]{xeRN\fitjMtjo}, i=\j=\ where the fcj are polynomials in the variables x\, ...,XN with coefficients in R and * j j is either < or =, for i = 1, ...,s and j = 1, ...,rj. See [BR90, BCR98] for an introduction on semi-algebraic sets. Examples of semi-algebraic sets include Bezier, Spline curves [Baj94] in geometric modelling, Geographic Information Systems, Computer Graphics and Aeronautics; the coupler curve [Mer96] in mechanism theory, workspace and singular configurations [Mou96] in robotics, etc. Maps are unthinkable without curved objects. Many problems from 1 different fields involve proximity queries like finding the nearest neighbour, finding all the neighbours, or checking or quantifying the neighbourliness of two objects. The potential applications of proximity queries on semi-algebraic sets em-bedded in two or three dimensional spaces include the problem of motion planning in a real environment for robots [Mou96] and Geographic Information Systems with curved objects being spatial primitives [Mio02]. The retraction planning problem [OY85, OSY86, OSY87] in robotics is strongly linked to questions of proximity among the projections of real-world objects in real-world environments onto the plane. The optimal path of a robot (considered as a disk) to transport the widest load while avoiding obstacles is a subset of the Voronoi diagram of those obstacles. Even if we assume that the obstacles a robot might encounter will not be curved in most of practical situations, the trajectories of the other robots or moving pieces that the robot may encounter, are curves that are expressed as semi-algebraic sets. The growth models used in several natural sciences are also strongly linked to prox-imity: the growth of crystal aggregates (the Johnson-Mehl model [OBS92]), the growth of trees in a forest (Voronoi diagrams [MB97]), etc. In geography, the study of influence zones and spatial analysis is also strongly linked to proximity queries on curved objects in the plane [VC90]. In all these applications, the qualitative know-ledge of the neighbourliness is more critical than the quantitative knowledge of the Voronoi diagram. In crystallography, the exact knowledge of the neighbourliness is necessary for the prediction of the crystallisation process. Voronoi diagrams are irregular tessellations of the space, where space is con-tinuous and structured by discrete objects [AKOO, OBS92]. The Voronoi diagram [Vor07, Vor08, VorlO] (see Figure 1.0.1) of a set of sites is a decomposition of the space into proximal regions (one for each site). Sites were points for the first his-torical Voronoi diagrams [Vor07, Vor08, VorlO], but in this thesis we will explore sets of circles, conies and more generally semi-algebraic sets. The proximal region corresponding to one site (i.e. its Voronoi region) is the set of points of the space 2 that are closer to that site than to any other site of the set of sites [OBS92]. We will recall now the formal definitions of the Voronoi diagram and of the Delaunay graph. For this purpose, we need to recall some basic definitions. Definition 1.0.2. (Metric) Let M be an arbitrary set. A metric on M is a mapping d : MxM —> R + such that for any elements o, b, and cof M , the following conditions are fulfilled: d (a, 6) = 0 a = b, d(a,b) = d(b,a), and d{a1c) < d(a,b) + d(b,c). (M,d) is then called a metric space, and d(a,b) is the distance between a and b. Remark 1.0.3. The Euclidean space RN with the Euclidean distance 5 is a metric space (RN,6). Let M = M.N, and 6 denote a distance between points. Let S = { s i , s m } C M,m > 2 be a set of m different subsets of M, which we call sites. The distance between a point x and a site Sj C M is defined as d (x, Sj) = iniyeSi {5 (x, y)}. Definition 1.0.4. (Bisector) For Si,Sj E S,Si^ Sj, the bisector B (SJ, SJ) of Sj with respect to Sj is: B(si,Sj) = {x G M\d(x,Si) = d(x,Sj)} (see example on Figure 1.0.2). Definition 1.0.5. (Influence zone) For Si,Sj € S,Si ^ Sj, the influence zone D(si,Sj) of Si with respect to Sj is: D(si,Sj) = {x G M\d(x,Si) < d(x,Sj)} (see example on Figure 1.0.3). Definition 1.0.6. (Voronoi region) The Voronoi region V(si,S) of Sj € <S with respect to the set S is: V(SJ,«S) = H^es S J - = / S J s j ) ( s e e e x a m p l e o n Figure 1.0.4). Definition 1.0.7. (Voronoi diagram) The Voronoi diagram of S is the union V (S) = ( J S i e S dV (si,S) of all region boundaries (see example on Figure 1.0.5). Definition 1.0.8. (Delaunay graph) The Delaunay graph DG (S) of S is the dual graph of V (S) defined as follows: 3 Figure 1.0.1: The ordinary Voronoi diagram (plain lines) of points (squares), and its topology expressed by the Delaunay triangulation (dashed lines) 4 Figure 1.0.2: The bisector (parabola) of a point and a line segment Figure 1.0.3: The influence zone (hashed) of a point with respect to a line Figure 1.0.4: The Voronoi zone (dark grey) of a cubic (light grey). The two objects are a circle and a cubic. 5 Figure 1.0.5: The Voronoi diagram (light lines) of a circle, an ellipse and a hyperbola (dark lines) • the set of vertices of DG (S) is <S, • for each (N — 1)—dimensional facet of V (S) that belongs to the common boundary of V (SJ,<S) and of V (SJ,S) with s», Sj € <S and Sj / Sj, there is an edge of DG (S) between Sj and Sj and reciprocally, and • for each vertex of V (S) that belongs to the common boundary of V(piltS),... ,V (siN+2,S), with VA; <E {1,...,N + 2} ,sik € S all distinct, there exists a complete graph between the Sik,k e { 1 , J V + 2}, and reciprocally (see example on Figure 1.0.6). The one-dimensional elements of the Voronoi diagram are called Voronoi edges. The points of intersection of the Voronoi edges are called Voronoi vertices. The Voronoi vertices are points that have at least N + 1 nearest neighbours among the sites of S. In the plane, the Voronoi diagram forms a network of vertices and edges. In the plane, when sites are points in general position, the Delaunay graph is a triangulation known as the Delaunay triangulation. In the plane, the Delaunay graph satisfies the following empty circle criterion: no site intersects the interior of 0 HYPERBOLA . •# '^5— v - -— — N * . ^^CTRCLE "x *" - ^ •» Figure 1.0.6: The Delaunay graph of the sites of Figure 1.0.5 the circles touching (tangent to without intersecting the interior of) the sites that are the vertices of any triangle of the Delaunay graph (see Figure 1.0.7). Once the Voronoi region a query point belongs to has been identified, it is easy to answer proximity queries. The closest site from the query point is the site whose Voronoi region is the Voronoi region that has been identified. The Voronoi diagram defines a neighbourhood relationship'among sites: two sites are neighbours if, and only if, their Voronoi regions are adjacent, or alternatively, there exists an edge between them in the Delaunay graph. The certified computation of the Delaunay graph is important for two reas-ons. B y certified computation, we mean a computation whose output is correct. First, unlike the Voronoi diagram, the Delaunay graph is a discrete structure, and thus it does not lend itself to approximations. Second, the inaccurate computation of this Delaunay graph can induce inconsistencies within this graph (see Sections 2.3 and 3.1.2), which may cause a program that updates this graph to crash. This is particularly true for the randomised incremental algorithm for the construction of the Voronoi diagram of semi-algebraic sets. The algorithm that certifies whether the facets of the Delaunay graph whose vertices are N + 1 given sites would remain Figure 1.0.7: The 3 empty circles for the sites of Figure 1.0.5 or not after the addition of a given site is central to the design of a semi-dynamic (allowing additions of sites but not deletions) algorithm for the construction of the Voronoi diagram for semi-algebraic sets. This algorithm is called the "Delaunay graph conflict locator'''' in the reminder of this thesis. Its input is a (N + 2)-tuple of sites, and its output is the list of all the Voronoi vertices corresponding to the (N — 1)—dimensional facets of the Delaunay graph having the first N + 1 sites as vertices that would not remain in the Delaunay graph after the addition of the (N + 2)^ n site, and a value that certifies the presence in that list (for each Voronoi vertex). The fact the addition of the site would imply the disappearance of a Delaunay graph facet is called a conflict. Thus, it justifies the name of "Delaunay graph conflict locator". In the context of the ordinary Voronoi diagram of points in the plane, the concept that is analogous to the Delaunay graph conflict locator is the Delaunay graph predicate, which certifies whether a triangle of the Delaunay triangulation would remain or not after the addition of a point. The development of the Delaunay graph conflict locator is the main objective of this thesis. The certified knowledge of the Delaunay graph for curved objects may sound like a purely theoretical knowledge that is not central in practical applications. 8 This is not the case in some applications. These applications include crystallo-graphy, metallography and V L S I layout. The Johnson-Mehl tessellations (which-generalise several weighted Voronoi diagrams) [OBS92] play a central role in the Kolmogorov-Johnson-Mehl-Avrami [Kol37] nucleation and growth kinetics theory. The Kolmogorov, theory provides an exact description of the kinetics during the heat-ing and cooling processes in crystallography (the Kolmogorov equation [Kol37]). The certified knowledge of the neighbourliness among molecules is central to the prediction of the formation of crystal aggregates. In metallography, the analysis of precipitate sizes in aluminium alloys through Transmission Electronic Microscopy [Des03, Section 1.2.2] provides an exact measurement of the cross sections of these precipitates when they are rodes with a fixed number of orientations [Des03, Section 1.2.2]. In V L S I design, the second order Voronoi diagram of the layout is used in the computation of the critical area, a measure of a circuit layout's sensitivity to spot defects [CPX02, Section 1], A n important concern on critical area computation is the robustness [CPX02, Section 1]. In the context of the application to robotics, one might object that the real-world objects a robot may collide with are approximate because of their man-ufacturing process. This is indeed true, but their specification is exact and can be expressed as semi-algebraic sets owed to the possibility of algebraic translation of the geometric specification of the methods of manufacturing process. Indeed, the mech-anical manufacturing process methods such as turning, countersinking, sharpening, drilling, bending, roll bending, or pressing can be expressed as geometric transform-ations, which in turn can be expressed algebraically. The finished result (mechanical component) is specified as a semi-algebraic set. Another limitation of approximative algorithms for the computation of the Delaunay graph is that when approximate computations are performed on objects defined approximately (within some geometric tolerance), the propagation of the errors can be critical, especially if the final computation involves approximate in-9 termediary computations. Interval analysis allows one to certify computations on objects defined approximately. Indeed, interval analysis allows one to consider ob-jects denned approximately since these objects can be specified by intervals. Interval analysis allows one to certify computations on intervals by providing bounds on the results. Finally, on another hand, the certified computation of the Delaunay graph participates to the recent move in the development of numerical and simulation software as well as computer algebra systems to exact systems [BCSS98]. 1.1 The problem The proximity queries stated above could be effectively answered if the Delaunay graph for sets of geometric objects could be computed in an efficient and certified way. This would require the embedding of the Delaunay graph and the location of the query point in that embedded graph. The embedded Delaunay graph and the Voronoi diagram are dual subdivisions of space, which can be stored in a quad-edge data structure [GS85]. The first and most explored Voronoi diagram is the Voronoi diagram for a set of points [Vor07, Vor08, VorlO] in the Euclidean plane or in the three-dimensional Euclidean space (see Figure 1.0.1). Voronoi diagrams have been generalised in many different ways including by modifying the space in which they are embedded (see [Aur87, OBS92] for a general survey of Voronoi diagrams): higher dimensional Eu-clidean spaces, non Euclidean geometries (e.g. Laguerre geometry, hyperbolic geo-metry, etc.). Fewer generalisations of Voronoi diagrams correspond to extending the possible sites from points to geometric objects. The only generalisations of this kind that have been explored are "abstract Voronoi diagrams" [Kle89] and the Voronoi diagram for lines [OBS92]: (the sites are points, line segments, cir-cular arcs or piecewise analytic curves). The Voronoi diagram for curved objects 10 does not necessarily satisfy the property of the abstract Voronoi diagrams. Ab-stract Voronoi diagrams are two-dimensional Voronoi diagrams defined by topolo-gical properties [Kle89] (which is that the bisector of any pair of sites is an unboun-ded simple curve). Therefore, the Voronoi diagram of N curved objects cannot be computed with the randomised O (N log N) algorithm for the construction of ab-stract Voronoi diagrams [Kle89], The Voronoi diagrams for lines that have already been studied are the Voronoi diagram for circles (set of sites comprising circles) [KKSOlb, K K S O l a , KKSOO], the Voronoi diagram for sets of points and straight line segments, the Voronoi diagram for sets of points, line segments and circular arcs [Yap87], and the Voronoi diagram for planar domains with curved boundaries (piecewise analytic) [RF99a, RF99b]. In principle, the Voronoi diagram can be gen-eralised to sets of sites comprising more general geometric objects (especially general curved objects). We will refer to these generalisations of the Voronoi diagram as the Voronoi diagram of sets of geometric objects (points, curves, surfaces; see Figure 1.0.5). These Voronoi diagrams of sets of geometric objects have been far less ex-plored, and they have been computed only approximately: as approximation al-gorithms decomposing the objects by points sets [OBS92] or approximating the computation of Voronoi vertices by a Newton-Raphson scheme for curves with ra-tional parameterisation [RF99a, RF99b]. B y computation of the Voronoi diagram, we mean computation of the coordinates of the Voronoi vertices, of the equations of the Voronoi edges, and of the network formed by theses vertices and edges. The first type of approximation algorithms for constructing Voronoi diagrams is not guaran-teed to give topologically correct results, because the Voronoi diagram is very sens-itive to the order of continuity at contact points (see Section 2.3 and [RF99a]). The second type of approximation algorithms for constructing Voronoi diagrams does not directly address the exactness of the Delaunay graph, because the basic compu-11 tation is the computation of the Voronoi vertex instead of the computation of the Delaunay graph predicate, and this predicate was not addressed in [RF99a, RF99b]. The neighbourhood relationship among sites is addressed indirectly through the identification of the Voronoi vertices and their classification. The curves that they addressed are parametric curves admitting rational parameterisations (i.e. ratios of polynomials), and therefore, it excludes conies (see Section 2.3). Moreover, their computation of the Voronoi vertices uses basic techniques (projective resultants) of higher complexity because these techniques consider common zeroes in the pro-jective space instead of in the affine space or in the algebraic torus (C*)N. This difference in complexity is explained later in the section on sparse elimination (see Section 3.3.1). This higher complexity of the algebraic computation techniques used is particularly significant because those projective resultant computations are not al-gebraic precomputations done only once, but computations done each time a bisector is computed. The algorithm for constructing the Voronoi diagram for points, line segments and circular arcs proposed in [Yap87] proceeds by a divide-and-conquer paradigm using vertical slabs, which excludes an incremental construction of the Delaunay graph. The computation of the Delaunay graph conflict locator for conies and more generally for semi-algebraic sets is the main problem that is being addressed in this thesis. The Delaunay graph conflict locator is the basic tool for maintaining the Delaunay graph when the curved objects are introduced sequentially. A direct application of this Delaunay graph conflict locator is a randomised incremental algorithm for the construction of the Voronoi diagram of semi-algebraic sets. This algorithm could be used for maintaining a semi-dynamic (allowing only insertions but not deletions) Voronoi diagram for semi-algebraic sets. Such a Voronoi diagram could be stored with the embedded Delaunay graph in a quad-edge data structure [GS85]. We will see how such an algorithm can be designed from the Delaunay 12 graph conflict locator in the following section. 1.2 T h e motivation In this section, we will examine how the Delaunay graph conflict locator can be used to maintain the Voronoi diagram of semi-algebraic sets in the plane as those semi-algebraic sets are introduced one by one. Finally, we will ennounce a necessary and sufficient condition for the connectivity of the Voronoi diagram of semi-algebraic sets in the projective plane, that has a direct application in the representation of spatial data at different resolutions. Knowing the Voronoi diagram V (<S) of a set S = { s i , . . . , sm} C R 2 of at least two semi-algebraic sets (m > 1) and its embedded Delaunay graph DG (<S) stored in a quad-edge data structure, we would like to get the Voronoi diagram V (S U {sm+i}), where sm+i is a semi-algebraic set of M 2 . In all this section, we will say that a circle C touches a semi-algebraic set s; if, and only if, C is tangent to S j and no point of s; is contained in the interior of C. The Voronoi edges and vertices of V (S) may or may not be present in V (5 U { s m + i } ) . Each new Voronoi vertex w induced by the addition of s m + i necessarily belongs to two Voronoi edges of V (S), because two of the three closest sites to w necessarily belong to <S. The new Voronoi edges induced by the addition of s m + i will clearly connect Voronoi vertices of V (S) to new Voronoi vertices induced by the addition of sm+i or new Voronoi vertices between themselves. Any of these later Voronoi edges e' must be incident with one of the former Voronoi edges at each extremity of e' (because the Voronoi vertex at each extremity of e' belongs to only one new Voronoi edge, i.e. e'). Any of the former Voronoi edges e must be a subset of a Voronoi edge of V (<S), since e must be a new Voronoi edge between sites of S (otherwise the Voronoi vertex belonging to V (S) at one of the extremities of e by the definition of e would be a new Voronoi vertex). Thus, to get V (S U {sm+i}), we need to know which Voronoi vertices and edges of V (5) will not be present in V (S U { s m + i } ) , which Voronoi edges of V (<S) 13 will be shortened in V (S U { s m + i } ) and which new Voronoi edges will connect new Voronoi vertices between themselves. We can test whether each Voronoi vertex v of V (S) will be present in V ( 5 U { s m + i } ) , .Let us suppose that v is a Voronoi vertex of Sj, Sj and s^- v will remain in V (S U {sm+i}) if, and only if, no point of s m + i is contained in the interior of the circle centered on v that touches Sj, Sj and s^. This is a sub-problem of the Delaunay graph conflict locator that can be tested by giving S{, Sj, Sk and s m + i as input to the Delaunay graph conflict locator, and then retain only the solutions where the Voronoi vertex is v. We can test whether each Voronoi edge e of V (<S) will be present in V (S U {sm+i}). Let us suppose that e is a locus of points having Sj and Sj as closest sites, e will disappear entirely from V (S U {sm+i}) if, and only if, a point of s m + i is contained in the interior of each circle centered on e and touching Si, Sj and each common neighbour Sk to Sj and Sj in DG(S) in turn. This can be tested by giving Si, Sj, Sk and s m +i as input to the Delaunay graph conflict locator and then retaining only the solutions where the Voronoi vertex belongs to e. e will be shortened (possibly inducing one or more new Voronoi edges) in V (<S U { s m + i } ) if, and only if, there exists Voronoi vertices of S{, Sj and sm+i on e and there is no point of any common neighbour s& to Sj and Sj in DG (S) in the interior of a circle centered on e and touching Sj, Sj and s m +i- The centre of each one of such circles will be a new Voronoi vertex in V (S U { s m + i } ) . This can be tested by giving s;, Sj, sm+i and s& as input to the Delaunay graph conflict locator and then retaining only the solutions where the Voronoi vertex belongs to e. This may not be the best way to proceed, but the Delaunay graph conflict locator is sufficient to maintain the Voronoi diagram of semi-algebraic sets. Tests might be limited to edges and vertices on the boundaries of the Voronoi regions V (si,S) ,Si G S that intersect sm+\ and of the Voronoi regions V (SJ, <S), Sj G S 14 adjacent to a Voronoi region V (SJ,<S). Indeed, a point (and thus a semi-algebraic set) can steal its Voronoi region only from the Voronoi region it belongs to and the adjacent Voronoi regions. We will finish this section with a necessary and sufficient condition for the connectivity of the Voronoi diagram of connected semi-algebraic sets in the project-ive plane. This result allows the characterisation of dangling edges in the Delaunay graph corresponding to the presence of closed edges in the Voronoi diagram. In or-der to proceed, let us recall some notations used in point set topology: let s denote the closure of s, and s denote the interior of s in the sense of the point set topology in R 2 . Note that if s bounds a closed domain then the interior of s is meant to be the interior of the closed domain bounded by s. P r o p o s i t i o n 1.2.1. (Connectivity of the Voronoi diagram in the plane) The Voro-noi diagram V (S) of a set S — {s\,..., s m } C M2 of at least two connected semi-algebraic sets (m > 1) considered in P 2 is not connected if, and only if, there exist a subset I of [1,... ,m] and one index j of [1,... ,m] such that Vi G / , Sj CSJ and VA; 6 [1, . . . ,m] \ I, si fl s% = s~j 0 s£ = 0. Proof. If: Assume there exist a subset / of [1 , . . . , m] and one index j of [ 1 , . . . , m] such that Vi G / , Sj Cs°- and Vfc € [1,...,m)\I, s~ir\s~k — SjCisJ = 0. Let si € S with I € [1, . . . , m] \ I. Let S = \JieI S{. Since S CSj, any circle touching both a S{, i G / and Sj must be contained in sj. Since SDJj — s~jHsJ = 0, no circle can touch each of a Si,i G /, Sj and s;. Thus, there is no point that has a Sj , i G /, Sj and s/ as nearest neighbours. Thus, there is no Voronoi vertex of a Si,i G /, Sj and s;. Since there is no Voronoi vertex of a Sj, i G /, Sj and an si with I G [1 , . . . , m] \ I, there are no Voronoi vertices on the bisector of 5 and Sj. Since S C\Jl = S C\sl = $, any circle centred on the bisector of S and Sj and touching both S and Sj does not intersect any site with k G [ 1 , . . . , rn] \ I. Thus, the bisector of S and Sj is contained in V (iS). Since Sj is connected and S Cs°, the bisector of S and Sj is a closed curve. 15 Thus, the Voronoi diagram of S is not connected in P2. Only if: Assume the Voronoi diagram of S is not connected in P2. Then, V (<S) has at least two connected components. Thus, at least one of these connected components does not have points at infinity. Let us consider the connected com-ponent (let us call it C\) that does not have points at infinity. Since C\ is composed of Voronoi edges1, each edge in C\ must end at either a Voronoi vertex or a point at infinity. Since C\ does not have any point at infinity, all Voronoi edges in' C\ connect Voronoi vertices. Thus C\ is a network of vertices and edges linking those vertices. The regions that this network defines are Voronoi regions. Let T> be the union of the closure of those Voronoi regions. V is a closed set set by its defini-tion. Let us consider now the semi-algebraic sets si,l 6 L whose Voronoi regions are contained in V. Let S = UzeL s z- From the definition of a semi-algebraic set, its is straightforward that the union of two semi-algebraic sets is a semi-algebraic set. Thus S is a semi-algebraic set. We will now consider S as a site instead of each o one of the si,l € L. The influence zone of S = (J;e£ s ' ^ s dearly V, because the influence zone of a union of semi-algebraic sets is clearly the closure of the union of the Voronoi regions of those semi-algebraic sets. Let e = dV. It is a portion of the bisector of S and another semi-algebraic set. Let us call it Sj. If not all the bisector of S and Sj was contained in V (S), then e would end at Voronoi vertices (a point on the Voronoi diagram has at least two closest sites) or the point at infinity, a contradiction with e not being connected. Thus, the bisector of S and of Sj. is contained in V (S), and it is equal to e. By the definition of e, e must be a closed curve. Assume the positions of S and Sj with respect to e are not always the same. Then, S and Sj must intersect. The bisector of S and Sj must have two branches near the intersection points (see Figure 1.2.1). Since e is a closed curve and S is contained in the interior of e, Sj must be closed, and the other branches must be unbounded (a contradiction with e not being connected in P2). Thus, the positions x a one-dimensional component of the Voronoi diagram, which is also the locus of points having two nearest sites 16 of S and Sj with respect to e are always the same along e. Since Sj is connected, S is contained in the interior of e and the positions of S and Sj with respect to e are always the same along e, S CSJ. Since e is the bisector of 5 and Sj and belongs to V (S), any circle centred on e and touching both S and s^ - does not intersect any site s/c with A; € [1, . . . ,m] \ I . Thus, V7c € [1, . . . ,m] \ n = s j n s]t = 0. • The only cases of disconnected (considered in P2) Voronoi diagrams corres-pond to one or more sites (semi-algebraic sets) contained in the interior of another site. This property has a direct application in Geographic Information Systems. When the same region 7Z bounded by a semi-algebraic set S is represented at dif-ferent scales, the representation of the details inside 1Z does not change the Voronoi diagram outside TZ. The edges of the Delaunay graph corresponding to a discon-nected Voronoi diagram (considered in P2) are respectively dangling edges or cut edges (the Delaunay graph is not bi-connected and removing a cut edge induces two connected components). It is possible to detect if there exists one or more sites Si,i £ I contained in the interior of another site Sj by checking that there exists no Voronoi vertex of S j , Sj and any € S distinct from Si and Sj. This is a subproblem of the Delaunay graph conflict locator. 1.3 Outline of the research The main theoretical objectives of this thesis are the determination of a general for-mula for degree of the offset to (i.e., the locus of points at a given distance from) an algebraic curve, and the reduction of the Delaunay graph conflict locator for algeb-raic curves from a semi-algebraic problem to a linear algebra problem (computing the eigenvalues of a matrix). The main practical objectives of this thesis are the computation of an implicit equation of the generalised offset (i.e., the locus of centres of circles of a given radius tangent) to a conic defined by a formal polynomial, the 17 Figure 1.2.1: The relative position with respect to the bisector must be constant 18 exact symbolic computation of the sparse resultant matrix for the Delaunay graph conflict locator for conies, and the certified computation of the Delaunay graph conflict locator for conies and for semi-algebraic sets. In all the computations, the central ideas were to use the lightest computational techniques and to simplify the formalisation of the solutions. The original contributions in this thesis are as follows. 1. A general formula for the degree of the generalised offset to an algebraic curve has been determined by studying the algebraic properties of these offsets. The knowledge of the degree of the generalised offset to a conic is used in the identification of the implicit equation of the generalised offset to a conic as a factor of a sparse resultant; 2. The number of points on which the Delaunay graph conflict locator for conies is evaluated has been computed. It corresponds to the number of lines of the matrix whose eigenvalues need to be computed for the algebraic computation of the Delaunay graph conflict locator for conies; 3. The Voronoi diagram and the Delaunay graph for circles have been computed exactly and symbolically through a completely symbolic conflict locator; 4. The computation of the Delaunay graph of conies has been reduced from a semi-algebraic problem to a linear algebra problem. The matrix for which the eigenvalues of a Schur complement of one of its submatrices give the answer to the Delaunay graph conflict locator for conies has been computed symbolically; 5. The computation of the Delaunay graph of semi-algebraic sets embedded in a two-dimensional space has been done using A L I A S . Although there is no known lower nor upper bound for this problem, the running time is satisfactory for exploration purposes. 19 We have also outlined how the Delaunay graph conflict locator could be used for the incremental construction of the Voronoi diagram of semi-algebraic sets. We have also proved a necessary and sufficient condition of connectivity of the Voronoi diagram for semi-algebraic sets in the plane, which has a direct application in Geographic Information Systems. The originality of this contribution resides in the formalisms that have been used (the generalised offset, see Chapter 4) and introduced (the generalised Voronoi vertex, see Chapter 5), and in the fact this work represents the first computation of the Delaunay graph conflict locator for general semi-algebraic sets. The thesis contents are organised as follows. Chapter 2 reviews the gener-alisations of the Voronoi diagram for curved objects. Voronoi diagrams for general semi-algebraic sets have not been studied previously. The problems closest to the Voronoi diagram for semi-algebraic sets are the Voronoi diagrams for manifolds, studied by Devillers et al. [DMT92] in 1992, the Voronoi diagram for curved ob-jects, studied by Al t and Schwarzkopf [AS95] in 1995, and the Voronoi diagram for planar domains with curved boundaries [RF99a, RF99b]. Chapter 3 introduces the approach used to design the computation of the Delaunay graph conflict locator for semi-algebraic sets. The approach used in this thesis involves combining symbolic precomputations and numerical computations to find the fastest computation of that conflict locator. The central tools for formalising the Delaunay graph conflict locator are the generalised offset (the locus of points that are locally at a given distance from a given geometric object [ASS99]) and the generalised Voronoi vertex (a concept introduced in this thesis: see Section 5.2). Wi th regard to symbolic computing, the approach includes working with the "geometry" of the monomials (see Section 3.3.1) composing a polynomial (i.e. representing the exponents of the monomials appearing in a polynomial as points in 20 an AT—dimensional space, where ./V is the number of variables we wish to eliminate) and the sparse resultant algorithms of Emiris [EC95, CEOO] and of Singular [GPS01], This is used to get an implicit equation of the generalised offset to a conic (see Section 4.4), and the matrix for which the eigenvalues of a Schur complement give the answer to the Delaunay graph conflict locator for conies (see Chapter 5). Wi th respect to scientific computing, the interval analysis and consistency methods implemented in A L I A S [MerOO] have allowed us to design and implement the Delaunay graph conflict locator for semi-algebraic sets and to obtain another Delaunay graph conflict locator for conies (by applying the interval analysis and consistency methods to the implicit equation of the generalised offset to a conic). Matlab eigs function has been used to compute the eigenvalues of the matrix ob-tained through symbolic computing for the Delaunay graph conflict locator for conies Chapter 4 studies the algebraic properties of the offset to an algebraic curve. From these algebraic properties, we have obtained a general formula for the degree of the offset to an algebraic curve. We applied this formula to conies in order to get the degree of the offset to conies. We used the sparse resultant algorithm of Emiris [EC95, CEOO] and the geometry of monomials to get an implicit equation of the generalised offset to a conic. Chapter 5 presents the algebraic computation of the Delaunay graph conflict locator for conies. Recent theoretical achievements from Algebraic Geometry and Computational Algebraic Geometry allow one to solve systems of algebraic equa-tions by solving a linear algebra problem (computing eigenvalues [CL098]). The eigenvalue problem has been studied for a long time and it is computationally tract-able [CDOO]. By computing in the quotient algebra [Lan02, Section 3, Chapter II] of the ring of polynomials by the ideal corresponding to a zero-dimensional algebraic variety, it is theoretically possible to transform a complex problem of resolution of 21 a system of algebraic equations into the more tractable and explored linear algebra problem of finding eigenvalues. Indeed, the quotient of the ring of polynomials in the variables of the problem by the ideal corresponding to the zero-dimensional variety object of the study is a finite dimensional vector space. Thus, every linear mapping can be expressed by a matrix in any (finite) basis of the quotient algebra. Moreover, the values of a given polynomial g at the points of the zero-dimensional variety correspond to the eigenvalues of the matrix of the map corresponding to the multiplication by g (see Theorem 4.5 on page 54 in [CL098]). This is the basis of the symbolic part of the Delaunay graph conflict locator for conies. The sparse resultant matrix corresponding to that conflict locator is a sparse matrix, and sparse methods for eigenvalues [Pin02, NikOO, vdV99, CDOO, Kra92] can be applied on it. Chapter 6 presents the interval analysis based computation of the Delaunay graph conflict locator for semi-algebraic sets. Interval analysis provides a more gen-eral approach for solving systems of algebraic equations and inequalities. Some new tools have appeared recently for solving systems of equations and inequalities with real coefficients based on interval analysis (see A L I A S [MerOO] and Section 3.3.3). A L I A S is a library developed at INRIA Sophia Antipolis, by Dr. Jean-Pierre Mer-let. A L I A S considers the real roots of zero-dimensional systems of equations and inequalities. A L I A S uses the P R O F I L / B I A S (Programmer's Runtime Optimized Fast Interval Library [Knii94]) library for evaluating intervals. It uses different, theorems from Real Algebraic Geometry [BCR98] for analysing as well as solving zero-dimensional semi-algebraic systems. The certified computation of the Delaunay graph conflict locator relies on theorems on the uniqueness of a root in given inter-vals (Kantorovitch, Moore-Krawczyk) and on the certified computation of function intervals by the P R O F I L / B I A S library. This computation uses a bisection process on one or all the variables using either only the equations of the system, or using the Jacobian of the system (Moore-Krawczyk test for finding "exactly" the solutions), 22 or using the Jacobian and the Hessian of the system (with Kantorovitch, Moore-Krawczyk tests). We first used A L I A S on the original system of algebraic equations and inequalities that specifies the Delaunay graph conflict locator for semi-algebraic sets. Then, we have obtained faster computations for conies by replacing the ori-ginal system by a system simplified by introduction of the implicit equation of the generalised offset to a conic. Chapter 7 summarises the results obtained in this thesis and presents the avenues of future work. 23 Chapter 2 Generalisations of the Voronoi diagram for curved objects In this chapter, we present an overview of the generalisations of the Voronoi diagram for curved objects. These are the generalised Voronoi diagrams- which are most closely related to the problem we are addressing in this thesis. These generalisations will be presented in increasing order of relevance for the present thesis. In Section 2.1, we will briefly introduce the Voronoi diagrams for manifolds, studied by Devillers et al. [DMT92] in 1992. In Section 2.2, we will introduce the Voronoi diagram for curved objects, studied by Alt and Schwarzkopf [AS95] in 1995. Finally, in Section 2.3, we will present in more detail the Voronoi diagram for planar domains with curved boundaries [RF99a, RF99b]. 2.1 V o r o n o i diagrams of general manifolds The Voronoi diagrams for manifolds were analysed in the context of the space of generalised spheres by Devillers et al. [DMT92] in 1992 (see Figure 2.1.1). In RN, an hypersphere S centred on $ € M.N and of radius r is mapped to a point S £ MN+1 (at the distance r 2 below the intersection of a vertical line passing through $ and the paraboloid of equation K — <f>2 — $r> = 0) of the space of generalised spheres (see 24 the upper part of Figure 2.1.1). The image of a concentric pencil of circles is the vertical line passing through its centre (see Figure 2.1.1). The sites of the Voronoi diagram for manifolds are manifolds of any dimen-sions embedded in RN [DMT92]. The distance induced by the metric in is denoted by S. An iV-dimensional topological manifold is a Hausdorff (or separated) space i.e., such that every point in it has an open neighbourhood homeomorphic to the open ball in RN. The distance is the usual distance from a point M to a manifold Z: 6 (M, Z) = mmP&z \MP\. Therefore, 6 (M, Z) is the radius of the min-imum sphere centred at M and "tangent" to Z (such that Z intersects the sphere but not its interior). Let Tz be the set of all the spheres tangent to Z. An example of this set Tz with Z being a circle S is shown on Figure 2.1.2 (the corresponding set is denoted Ts). In the space of generalised spheres S(E), if Z is an analytic manifold (i.e. such as all the connecting maps are infinitely often differentiable), then Tz is an analytic manifold (see [DMT92, page 20]): it is the upper envelope of the polar1 planes of the point-spheres (points) of Z with respect to the earlier mentioned paraboloid. Let S be a set of sites (manifolds). Let Us be the upper envelope of the manifolds Tz for all Z e S. The intersection of the Us with a vertical line <5 = M in the space of spheres gives the sphere with centre M whose radius is the distance to the nearest neighbour of M in S. The projection of the Us on the n-dimensional Euclidean space is the Voronoi diagram of <S. However, by bounding the upper envelope Us, what corres-ponds to considering only the spheres that are contained in a big sphere enclosing all the sites, we get a compact convex set. The work of Devillers et al. [DMT92] did not include any algorithm for the construction of the Voronoi diagram of general manifolds. 1The polar plane of a point M with respect to a quadric Q is the locus of the harmonic conjugates of M with respect to the two intersections of a line through M intersecting Q with Q. 25 26 Figure 2.1.2: The circles tangent to a given circle in the space of generalised spheres (taken from [DMT92]) 27 2.2 Voronoi diagrams for curved objects A n incremental randomised algorithm for the construction of the Voronoi diagram for curved objects in the Euclidean plane has been proposed by Al t and Schwarzkopf [AS95], and it is motivated by the applications in motion planning "which lead to the so-called retraction method" [AS95]. The approach of A l t and Schwarzkopf is to decompose the complex curves into open curves that are characterised by the property that no circle touches them in more than one point. Such curves are called "harmless". This decomposition is motivated by the fact that the Voronoi diagram for connected curved objects is not necessarily connected and there may be Voro-noi edges between two different parts of the same site (these edges are called self Voronoi edges) and Voronoi vertices between three different parts of the same site (these Voronoi vertices are called self Voronoi vertices). The self Voronoi edges are loci of centres of circles tangent to without containing different parts of the same site, while self Voronoi vertices are centres of circles tangent to without containing different parts of the same site. The curves are broken up into "harmless" pieces to prevent the computation of self Voronoi vertices and self Voronoi edges by the ran-domised incremental algorithm for the construction of the Voronoi diagram. They prove that this decomposition into harmless pieces assures that the Voronoi diagram is connected, no region is empty, and each region is simply connected. The first basic assumption of this approach is that curves are abstract ob-jects, and "certain elementary operations are available as black boxes" [AS95]. These operations are the following constructions: "finding the points having the same dis-tance from three given points, finding all points of a given slope, finding points where the curvature has a local maximum, finding the representation of a bisector given the representation of two curves, and finding intersection points of given curves" [AS95]. The curves are encoded as their parametric representation (7 : / —» R 2 where / C R is some closed interval). The second basic assumption is that curves 28 are supposed to be regular and simple. They define a harmless site as either a point, an open circular arc, or a harmless curve (which is also open by definition); and they define a harmless site collection as a finite set S of pair-wise disjoint harmless sites with the condition that for every circular arc and harmless curve of S, its endpoints are also members of S. This definition of harmless sites excludes closed straight line segments and circular arcs as well as circles (a circular arc should be open, and a circle has to be partitioned into a point and an open circular arc). The points responsible for the non-harmlessness of curves are the local maxima of the absolute value of the curvature (see proof at the beginning of section 3 in [AS95]). In the incremental randomised algorithm for constructing the Voronoi dia-gram for curved objects in the plane, the insertion of the sites is done in two steps. First, a point acting as place holder is computed for each one-dimensional harmless site. The point objects and these place holders are inserted first in a random order. Then, the one-dimensional harmless sites are inserted in a random order. The pre-dicates and constructions needed for this incremental randomised algorithm have been considered as "black boxes", and no implementation results have been presen-ted. Algebraic curves are harmful, and the decomposition of an algebraic curve into harmless sites may require a number of cuts bounded by its degree minus one. 2.3 The Voronoi diagram and medial axis transform for planar domains with curved boundaries We will review the work of Rajesh Ramammurthy and Rida T. Farouki [RF99b, RF99a]; of Rida T. Farouki and Rajesh Ramammurthy [FR98b, FR98a]; and of Rida T. Farouki and John K . Johnstone [FJ94b, FJ94a] on Voronoi diagrams for planar domains with curved boundaries, and point-curve and curve-curve bisectors. This review will also consider a result from [CCM97] stating the sensitivity of the 29 Voronoi diagram to the order of continuity of the approximations at contact points that is cited in some of the earlier mentioned papers. The Voronoi diagram and medial axis of a closed bounded planar domain are basic geometric entities associated with that domain. The Voronoi diagram of a domain bounded by N curve segments is a network specifying a partition of the plane into N regions such that each point within a given region is at least as close to its associated curve segment as to all other curve segments. The me-dial axis is the locus of centres of maximum-radius circles that may be inscribed within the domain. This locus forms also a network. The computation of these geo-metric entities involves the computation of the nodes and edges of these networks. The edges of these networks are part of point/curve, curve/curve, or self bisectors. The nodes of these networks are the centres of circles touching the boundary at at least two distinct points. The basic assumption (and focus of the work) of Rida T. Farouki et al. [RF99b, RF99a, FR98b, FR98a, FJ94b, FJ94a] is that planar domains are considered with piecewise analytic boundaries. For domains with poly-gonal or piecewise linear/circular boundaries, the bisectors are conies, and efficient algorithms that yield essentially (i.e. topologically) exact Voronoi diagrams have been developed. These bisectors admit rational parameterisations. However, for planar domains bounded by curves having a polynomial or rational parameterisa-tion, the corresponding curve/curve bisectors do not necessarily admit such simple rational parameterisations. One possible approach for computing the bisectors of those curves having a rational parameterisation is to use approximations by the earlier mentioned piece-wise linear/circular curves. However, this yields results "that are not even qualitat-ively (topologically) correct" [RF99a]. The argument presented in [RF99a] involves two counter examples. In one of them, "the discrepancy between the true Voronoi 30 Figure 2.3.1: The difference between (a) the true Voronoi diagram of a planar domain bounded by curves, and the computed from piecewise-linear boundary ap-proximations with (b) 15 segments and (c) 60 segments (taken from [RF99a]) diagram and the approximated one grows as the tolerance on the approximation is tightened by introducing further linear/circular approximating segments" (see Fig-ure 2.3.1). In the other one, "the Voronoi diagram and the medial axis are identical for the exact boundary while they differ for the approximate boundary" (see Figure 2.3.2). The explanation that is given to justify this strange behaviour is that both the Voronoi diagram and the medial axis are very sensitive to the order of continuity of the boundary curve segments at their contact points [CCM97]. The approach used by Ramammurthy and Farouki [RF99b, RF99a] is: • to determine the rational parameterisations of the Voronoi edges that admit them, and • to provide piecewise-rational approximations that satisfy a prescribed geomet-rical tolerance for the remaining Voronoi edges. The main results that are at the basis of their algorithm are: • the bisector of a point and a rational curve segment can be described exactly (e.g. in the customary Bezier form [FJ94b]); 31 . Figure 2.3.2: The internal part of the Voronoi diagram and the medial axis are identical for the exact boundary while they differ for the approximate boundary (taken from [RF99a]) • the bisector of two rational (or even polynomial) curves is not necessarily a rational locus, but by expressing the curve/curve bisector as the envelope of a family of point/curve bisectors (see Figure 2.3.3), the generation of the sequences of singularities (tangent discontinuities) of the curve/curve bisector can be reduced to a family of univariate polynomial root-finding problems, • the true curve/curve bisector can be approximated to any given geometric tolerance by means of adaptive subdivision and error measures for geometric Hermite interpolants, and its "singularities (tangent discontinuities) can be captured in an essentially exact manner" [FR98b]. This work allows one to only construct the Voronoi diagram approximately since some Voronoi vertex computations are approximate. Moreover, the compu-tation of the bisectors is very heavy since it involves discretising both curves, and then computing the envelopes of the two families of point-curve bisectors. The only algebraic curves that can be processed with the proposed projective resultant based techniques are the algebraic curves admitting rational parameterisations. Finally, 32 Figure 2.3.3: The bisector of two curves as the envelope of a family of point/curve bisectors (taken from [FR98b]) 33 this work does not directly address the problem of computing exactly the Delaunay graph conflict locator. We will review the results presented above as well as the algorithm for com-puting the Voronoi diagram. We have presented this work in further detail hereafter because it is the only implementation of an algorithm for constructing the Voronoi diagram of quite general curved objects. The algorithm for computing the Voronoi diagram of a planar domain bounded by a rational curve is incremental: one boundary segment is introduced at a time. To perform the Boolean operations involved in the incremental construction of the Voronoi diagram, a complete description of the oriented, boundaries of the involved regions is necessary and sufficient. Such a description involves a classification of the Voronoi vertices. While some kinds of Voronoi vertices involve only rational bisectors, and they can be computed by standard curve intersection algorithms, which are essen-tially exact; other kinds of Voronoi vertices involve non-rational bisector segments, which must be approximated. Therefore, their location computed as intersection of non-rational bisectors is inherently approximate. Such non-rational bisectors are approximated using a polynomial interpolation method that guarantees the desired order of continuity at contact points (the Hermite interpolant). Moreover, the com-puted Voronoi vertices are refined through a Newton-Raphson scheme. The computation and classification of the Voronoi vertices are done by in-termediate computations of curve/curve bisectors. These curve/curve bisectors are generally portions of high-order algebraic curves that do not admit rational para-meterisations. The approach of Farouki and Ramammurthy is to approach these segments to a prescribed geometrical tolerance. In order to do so, they provide a means of recognising transition points between segments belonging to different types and singularities. 34 The philosophy of the method to computing curve/curve bisectors is to break the curve/curve bisector into a sequence of segments between singular points, and then to transform the problem of generating ordered sets of points along the bisector to a sequence of univariate polynomial root-finding problems, using point/curve bisectors as an intermediate tool. Each curve is discretised in turn. The curve/curve bisector is computed as the curves formed by the left and right candidate points (see Figure 2.3.4). The left and right candidate points are the points on the left and right side of the discretised curve respectively that belong to the envelope of the family of point/curve bisectors. The left and right candidate points are the centres of circles passing through the point on the discretised curve and tangent to the other curve or passing through an extremity of that other curve. The locus of the left and right candidate points is a superset of the curve/curve bisector, called the untrimmed bisector. The untrimmed curve/curve bisectors are trimmed by identifying the values between discrete samples on the untrimmed bi-sector at which there is a corresponding change in the status for the candidates: re-tained/discarded, i.e., belonging/not belonging to the trimmed bisector. The Voro-noi vertices and the singular points on the curve/curve bisector are also identified. For each retained point, the unit tangent and curvature of the curve/curve bisector are computed. Hermite interpolants are used to construct the curve/curve bisectors between Voronoi vertices and/or singular points. If the Hermite inter-polants do not satisfy the specified tolerance, then additional intermediary points may be added. Finally, the point/curve bisectors are computed by trimming in a similar way the point/curve bisectors (which admit rational parameterisations). The problem we address in this thesis is the certified computation of the 35 Figure 2.3.4: The left and right candidate points (taken from [FR98b]) 36 Delaunay graph of sets of semi-algebraic sets. The fact that we address the Voro-noi diagram of semi-algebraic sets makes of our work an extension of the work on Voronoi diagrams for planar domains with curved boundaries. Indeed, we consider semi-algebraic sets which include algebraic varieties, which in turn strictly include algebraic curves with rational parameterisations. The proper conies (circle, ellipse, parabola and hyperbola) are the simplest example of algebraic curves that do not ad-mit rational parameterisations, and that cannot be processed by the work presented in this section. This results from that the resultant techniques used on the rational parameterisations exclude curves without an algebraic rational parameterisation. 37 Chapter 3 Combining symbolic computation and scientific computation In the first section, I will present the Delaunay graph conflict locator for additively weighted points and for circles. This presentation will allow me to introduce in an easier way the reader to the formalisms of the generalised offset (which will be studied in Chapter 4), and of the generalised Voronoi vertex (which will be introduced in Chapter 5). In the second section, I will present different attempts made to compute the Delaunay graph conflict locator by using a formulation of the conflict locator that was based on the original curves, both using symbolic computational techniques and numerical computational techniques. These were the first attempts at computing the Delaunay graph conflict locator. These experiments suggested that a purely algebraic solution starting from the original curves was not tractable (except for the simple case of the Voronoi diagram of circles). In the third section, I will briefly introduce the key computational techniques used in this research that constitute a hybrid approach integrating symbolic algeb-38 raic precomputations with numerical computational techniques for finding eigen-values and for solving systems of equations and inequalities. The different aspects of this hybrid symbolic-numerical approach will be introduced in Chapters 5 (for conies) and 6 (for arbitrary semi-algebraic sets). Generally, one performs first al-gebraic precomputations, and then take over with numerical computations. The main challenge is to identify where the symbolic precomputations should stop, or alternatively, where the numerical analysis techniques should start. Two central ideas have driven this thesis. The first one is that knowing the structure of the set of solutions may help finding the solutions. For the structure of the solution set, algebraic geometry plays a central role. The second one is that some (one time) symbolic preprocessing may accelerate the certified numerical evaluation of the Delaunay graph conflict locator to be performed several times. The main computational challenges that have been encountered during this thesis are identifying which computational techniques might work and which com-putational techniques will probably not work, and finding the optimum between the legitimate wish to compute everything exactly, symbolically and algebraically, and the universal scope of the numerical computational methods for solving systems of equations and inequalities. We will see that although the Delaunay graph of sets of circles can be com-puted symbolically and exactly, and the Delaunay graph of sets of algebraic curves can be specified theoretically, the exact computation of the Delaunay graph of sets of algebraic curves is beyond the present limits of exact computational methods such as Grobner bases [Gr639, Buc92, Buc70, Buc79, Buc88, B C K 8 8 , Buc98] and classic projective resultants. Indeed, the complexity of the computation of Grobner bases is doubly exponential in the number of variables (see the doubly exponen-tial lower bound in [MM84, Huy86] and the doubly exponential upper bound in [MM84, Giu84]), and the complexity of the computation of the sparse resultant 39 is exponential in the number of variables [Emi96]. This exponential complexity of the computation of Grobner bases or resultants does not affect the complexity of the evaluation of the conflict locator because those computations are algebraic precomputations done only once before implementing the algorithm for evaluating the Delaunay graph conflict locator. So, it makes sense to try to attempt such a huge complexity algebraic precomputations. However, the most limitative resource for those algebraic precomputations is the amount of Random Access Memory and of Virtual Memory accessible. Also, the size of the matrix of the sparse resultant determines the complexity of the numerical computations for finding eigenvalues performed each time the Delaunay graph conflict locator is computed. 3.1 The exact symbolic Delaunay graph conflict locator for circles We will first present the exact symbolic Delaunay graph conflict locator for additively weighted points when weighted points are introduced one by one, and then introduce what changes for circles. For this purpose, we will present some preliminaries about Additively Weighted Voronoi diagrams. 3.1.1 Pre l iminar ies Let N be the set of integers, E be the set of real numbers, and R 2 be the Euclidean plane. Let V = { P i , P J V } be the set of generators or sites, where P; is the weighted point located at pi € R 2 and of weight Wi € R. Let Cj be the circle centred at pi and of radius Wi, which we call weight circle hereafter. The definitions of bisector, influence zone, Voronoi region and Voronoi dia-gram presented in Chapter 1 generalise to the case where the set of sites <S is a set of weighted points V, and the distance d(M, Pi) (called additive distance) between a point M and a site Pi is d ( M , Pj) = S (M,pi) — Wi, where 8 is the Euclidean distance 40 between points. The Voronoi region of Pi with respect to the set V is defined by: V (PuV) = {M e M 2 | V j ^ i : S (M,pi) - wt < 6 (M,Pj) - W j } . The Additively Weighted Voronoi diagram of V is defined by: V (V) = \JP.ePdV (Pi,V). The Additively Weighted Voronoi diagram is illustrated on Figure 3.1.1: the weight circles are drawn as plain disks with a small hole at their centres, the Additively Weighted Voronoi diagram is drawn in plain thick hyperbola segments, and the Delaunay graph is drawn in dashed lines. The Additively Weighted Voronoi diagram defines a network composed of edges (loci of points having two nearest neighbours), and vertices (loci of points having three nearest neighbours). The Additively Weighted Voronoi diagram is related to the Apollonius Tenth problem. The Apollonius Tenth problem is to find a circle T tangent to three given circles C\, C2, C3 (see Figure 3.1.2). For additively weighted points; we will see later in this section that only the circles that are either externally tangent to each of three given circles C\, C2, C3 or internally tangent to each of C\, C2, C3, are relevant to the Delaunay graph conflict locator. The centres of the circles that are solutions to the Apollonius Tenth problem are the first example encountered in this Figure 3.1.1: The Additively Weighted Voronoi diagram 41 thesis of generalised Voronoi vertices (a concept that we will introduce in Section 5.2). Informally, generalised Voronoi vertices are the centres of circles tangent to N + 1 sites, where N is the dimension of the Euclidean space. Hereafter we will call the solutions of the Apollonius Tenth problem Apol-lonius circles. The centres of the Apollonius circles that are either externally tangent to each of three given circles G\, C% C3 or internally tangent to each of C\, C2, C3 are the first example encountered in this thesis of true Voronoi vertices (i.e. centres of circles that are touch N + 1 sites where N is the dimension of the Euclidean space). 3.1.2 The Delaunay graph conflict locator for addi t ively weighted In this subsection, we present an exact algebraic conflict locator for the Delaunay graph of additively weighted points (i.e. the dual graph of the Additively Weighted Voronoi diagram). The maximum degree of the polynomials which need to be eval-uated to compute this Delaunay conflict locator is 16 (thus, we say that the degree of the conflict locator is 16). This Delaunay graph conflict locator would be the core of a randomised incremental algorithm for constructing the Additively Weighted Figure 3.1.2: The Apollonius Tenth problem points 42 Voronoi diagram since the Additively Weighted Voronoi diagram is an abstract Voronoi diagram [Kle89], and thus, it can be constructed with the randomised in-cremental algorithm of Klein [Kle89]. This work solves the robustness issue in the work of Anton, Mioc and Gold [AMG98] on dynamic Additively Weighted Voro-noi diagrams by providing an exact conflict locator. The exact computation of the Additively Weighted Voronoi diagram has not been addressed until Anton et al. [ABMY02]. That paper addressed the exact predicate for the off-line construction of the dual graph of the Additively Weighted Voronoi diagram from the dual of the Power Voronoi diagram of spheres by using the relationship between the Additively Weighted Voronoi diagram in the plane and the Power Voronoi diagram 1 of spheres in the three-dimensional space. In their independent work, Karavelas and Emiris [KE02] provided several exact predicates of maximum degree 16 for achieving the same "in-circle/orientation/edge-conflict-type/difference of radii" test as we do in a single conflict locator. They reduced the degree of their predicate from 28 to 20 and then to 16 using Sturm sequences and invariants. The motivation for an exact conflict locator lies in the fact that without an exact computation of the Delaunay graph of additively weighted points, some geometric and topologic inconsistencies may appear. This is illustrated with an example. The starting configuration is shown on Figure 3.1.3. There are three weighted points (whose corresponding weight circles are drawn). The Delaunay graph is drawn in dashed lines. The Apollonius circles tangent to the weight circles have been drawn in dotted lines. The real configuration after addition of a fourth weighted point is shown on Figure 3.1.4. The configuration that might have been computed by an approximate algorithm is shown on Figure 3.1.5: the difference between real and perceived situations has been exaggerated to show the difference. 1The Power Voronoi diagram is a generalised Voronoi diagram where sites are hyper-spheres and the distance between a point and a site is the power of that point with respect to that site [Aur87]. 43 The old Apollonius circles have been adequately perceived to be invalid with respect to the newly inserted weighted point. About the new Voronoi vertices, while on the right of the figure two new Voronoi vertices have been identified as valid with respect to their potential neighbours, on the left of the figure, only one Voronoi vertex has been identified as being valid with respect to its potential neighbours. While the new Voronoi edge between the middle and bottom weighted points can be drawn between the two new Voronoi vertices of the new, middle and bottom weighted points; the Voronoi edge between the top and new weighted points cannot be drawn, because there is no valid Voronoi vertex on the left. There is an inconsistency within the topology: there is one new Voronoi vertex (the Voronoi vertex of the top new and middle weighted points) that cannot be linked by a new Voronoi edge to any other new Voronoi vertex and thus, that Voronoi vertex is incident to only two Voronoi edges. That additively weighted Voronoi diagram that might have been computed by an approximative algorithm is not an additively weighted Voronoi diagram. Thus, even if we perturbate the input weighted points, we will never get that additively weighted Voronoi diagram. We consider the maintenance of the Delaunay graph of additively weighted points in an incremental way: we check the validity of all the triangles of the Delaunay graph whose vertices are P i , P2, P3 with respect to a newly inserted weighted point P 4 [AKM02]. Thus, the input of the conflict locator is constituted by four points: the first three are supposed to define a triangle in the Delaunay graph, and the last one is the newly inserted weighted point. Let (xi,yi) be the co-ordinates of pi, for i. = 1,2,3,4. There are two possible outcomes to the above test of validity: either the triangles are valid with respect to the newly inserted weighted point and the triangles remain in the Delaunay graph, or one or two triangles are not valid with respect to the newly inserted weighted point and those triangles will not be present in the Delaunay graph any longer. We can see an example of the later case in Figure 3.1.6. A triangle having P1P2P3 as vertices is not valid with respect 44 Figure 3.1.3: The starting configuration to the weighted point P4. Thus, it will not belong any longer to the Delaunay graph after the insertion of P4. The conflict locator consists of determining which ones of the additively weighted Voronoi vertices of P\, P2 and P3 will not remain after the insertion of P4. This is equivalent in turn to the additive distance from which ones of the additively weighted Voronoi vertices of P i , P2 and P3 to P4 is smaller than the additive distance of that Voronoi vertex to P i (or P2 or P3, see Figure 3.1.6). Any additively weighted Voronoi vertex I oi Pi, P2, and P3 with coordinates (x, y) can be obtained algebraically by computing the common intersection of the three circles C[, C 2 and C 3 expanding (see Figure 3.1.7), or shrinking (see Figure 3.1.8) from the three first circles C\, C2 and C3 all at the same rate. The common 45 Figure 3.1.4: The real configuration after addition of the fourth weighted point (bold weight circle) signed expansion of the first three circles is denoted by r. Each circle C" centred on (x, y) and of radius r is either externally tangent to the first three circles (if the expansion r is positive) or internally tangent to the first three circles (if the expansion r is negative). The centres coordinates x,y and radii r of the circles C" centred on the intersections I = C[ D C'2 n C'3 and either externally or internally tangent to each of Ci, C2, and C 3 can be computed algebraically as the solutions of the following system of three quadratic equations in the variables x, y and r: 46 Figure 3.1.5: The configuration computed by an approximate algorithm cfT (x,y,r) = (x- x i ) 2 + {y- y\f - (w1 + rf = 0 < c'2 (x, y, r) = (x - x2)2 + (y - J^)2 - («>2 + r) 2 = 0 4 (x, y, r) = (x - x 3 ) 2 + (y - y 3) 2 - ( ^ 3 + r) 2 = 0 Subtracting one of the equations (say dx (x, y, r) = 0) from the remaining two (c2 (x,y,r) = 0 and c3 (x,y, r) = 0) results in a system of 2 linear equations, from which x and y may be expressed as linear functions of r. Substitution in the first equation dx (x,y,r) = 0 then leads to a quadratic equation in r. This means that the unknown quantities x, y, r can be expressed with quadratic radicals as functions of the given centres and radii. Though the simplest thing to do now would be to compute the two Voronoi vertices and use their computed coordinates and corresponding signed expansion in 47 Figure 3.1.6: The Delaunay graph conflict locator for the Additively Weighted Voro-noi diagram the computation of the values certifying the output of the Delaunay graph conflict locator, it is not desirable because this method would not be generalisable to conies or higher degree algebraic curves. We will detail hereafter only the computation of the values certifying the presence in the output list. To get the exact Delaunay graph conflict locator in a more elegant and generalisable way, we evaluated the values certifying the conflict locator output without relying on the computation of the Voronoi vertices as an intermediary computation. This is done by evaluating the values taken by the polynomial function expressing the relative position of C\ with respect to C" on the set of solutions of the system (i.e. the common zeroes of 48 Figure 3.1.7: The Additively Weighted Voronoi vertex as the common intersection of three expanding circles Figure 3.1.8: The Additively Weighted Voronoi vertex as the common intersection of three shrinking circles the three polynomials c[, c 2 and c 3 ). This is possible thanks to the translation that exists between geometry and algebra. More specifically, to the geometric set X of the set of common zeroes of the three polynomials c'x, c 2 and c 3 in Kz, where K is an algebraically closed field [Lan02, Definition before Theorem 1, Section 2, Chapter VII], we can associate the set of all polynomials vanishing on the points of X, i.e., the set of polynomials fic[ + /2C2 + /3C3 where the fc, i = 1,2,3 are polynomials in the three variables x,y,r with coefficients in K. This set is the ideal [GP02, Defin-ition 1.3.1] {c i , c 2 , c 3 ) . The set of polynomials with coefficients in K forms with 49 the addition and the multiplication of polynomials a ring: the ring of polynomials [GP02, Definition 1.1.3]. It is easy to see that a polynomial function g(x,y,r) on K3 is mapped to a polynomial function on X if we recursively subtract from g any polynomial in g belonging to (c^c^Cg) until no monomial in g can be divided by each one of the lexicographically highest monomials in c'x, c 2 and c 3 . The result of this mapping gives a canonic representative of the reminder of the Euclidean division of the polynomial g by the polynomials c[, c 2 and c 3 . The image of the ring of polynomials by this mapping is called the quotient algebra [Lan02, Section 3, Chapter II] of the ring of polynomials by the ideal (c^c^c^) . It is also easy to see that (c'i,c'2 — c^c^ — c'1)=(c'1, c 2 , c 3 ). Moreover, if we recursively subtract from g any polynomial in g belonging to (dx, c 2 — c\, c 3 — c[) t i l l the only monomials in g are 1 and r, we get the same result as the preceding mapping. The polynomials c[, c 2 — c[, c'3 — Ci constitute what is called a Grobner basis [GP02, Definition 1.6.1] of the ideal (c'^c^c^). The monomials 1 and r are standa-rd monomials. Grobner bases are used in Computational Algebraic Geometry in order to compute a ca-nonic representative of the remainder of the division of one polynomial by several polynomials generating a given ideal / . This canonic representative belongs to the quotient algebra of the ring of polynomials by the ideal / . The Grobner basis for this system provides a set of polynomials that define uniquely the algebraic rela-tionships between variables for the solutions of the system. The initial (largest with respect to some monomial order [CL098]) monomials of each one of the polynomi-als of the Grobner basis form an ideal. The monomials that do not pertain to this ideal form a basis for the representatives of the equivalence class of the remainders of the division of a polynomial by the polynomials of the system in the quotient algebra. These monomials are called standard monomials. The size of this basis equals the dimension [GP02, see definition on page 414] of the quotient algebra and the number of solutions of the system counted with their multiplicity [Lan02]. In the case of the conflict locator for the additively weighted Voronoi diagram, there 50 are two solutions. The polynomial g = (£4 — x)2 + (2/4 — y)2 — (r + r^)2 expresses the relative position of C\ with respect to C". Indeed C" is tangent to C 4 if, and only if, the Euclidean distance between the centres of C" and of C4 (i.e., (x,y) and ^4) equals the sum of the radii r and r^, i.e. (2:4 — x)2 + (y 4 — y)2 — (r + T 4 ) 2 = 0. The open balls bounded by C" and C4 intersect if, and only if, the Euclidean distance between the centres of C" and of C4 is smaller than the sum of the radii r and r^, i.e. (X4 — x)2 + (2/4 — y)2 — (r + r 4 ) 2 < 0. The circles C" and C\ are disjoint if, and only if, the Euclidean distance between the centres of C" and of C4 is greater than the sum of the radii r and r\, i.e. (X4 — x)2 + (^4 — y)2 — (r + r 4 ) 2 > 0. We considered the operation of multiplication of polynomials by the polynomial g. This multiplication operator is a linear mapping. The operation of this mapping on the canonic representative of the reminder of the division of a polynomial by c[, c'2 and c'3 is also a linear mapping that can be expressed by a matrix since the quotient algebra has a finite dimension. ( ^00 moi \ of the following multi-7nio rnn J plication operator on the quotient algebra: mg : [/] [gf]. The eigenvalues of Mg are the values of g taken on X (see Theorem 4.5, page 54 in [CL098]). The eigenvalues of Mg are the solutions of det (Mg — XI) = 0, where / denotes the 2 x 2 identity matrix, i.e. the roots of A 2 - A (m0o + m i l ) + (moomii) - (m 0 imio) = 0 (3.1.1) The values certifying the presence in the list output by the Delaunay graph conflict locator are the signs of the values taken by g, and they are determined by the sign of the roots of Equation 3.1.1 (which are the eigenvalues of Mg). If there is only one eigenvalue and it is 0 then the fourth circle is tangent to the circle externally tangent to the first three circles. The sign of A (where A = 51 ("loo + m l l ) 2 — 4(moomoi — moimio) ) cannot be negative because this would be equivalent to the fact there would be no triangle with vertices C\, C% and C 3 in the old Delaunay graph (because of the absence of real Voronoi vertex, see Figure 3.1.9). Thus, sign (A) is 0 or positive, and we have to evaluate the sign of the roots of the quadratic equation. Figure 3.1.9: There is no such triangle in the old Delaunay graph because of the absence of a real Voronoi vertex When there is only one double root of Equation 3.1.1 then we have the following two possibilities. Either the value of the root of Equation 3.1.1 is positive or 0 and the triangle will remain in the new Delaunay graph, or the value of the root of Equation 3.1.1 is negative and the triangle will disappear in the new Delaunay graph (see Figure 3.1.6). When there are two real roots of Equation 3.1.1, we have two triangles to consider (see Figure 3.1.10). The triangles that correspond to the roots with a negative value will disappear in the new Delaunay graph (see Figure 3.1.10). There is not much interest in showing the elements of the matrix of the multiplication operator here, but the Macaulay 2 [GS] code is presented in Appendix A . The exact algebraic computation of the Delaunay graph conflict locator we have presented in the previous paragraph is not generalisable to the other proper conies or 52 Figure 3.1.10: Two triangles can possibly disappear simultaneously by the addition of a single weighted point higher degree algebraic curves. Indeed, the size of the multiplication operator matrix is greater than 4 for the other proper conies and for higher degree algebraic curves (see Section 5.3), and an algebraic equation of degree 5 or more is not necessarily solvable by radicals (see [BB96, Theorem 8.4.8]). Even if we can obtain the matrix of the multiplication operator symbolically, we will need numerical methods for computing the eigenvalues of that matrix, which give the answer to the Delaunay graph conflict locator. We will now present the Delaunay graph conflict locator for circles, emphasising on the changes with respect to the Delaunay graph of additively weighted points presented in this subsection. 3.1.3 T h e Delaunay graph conflict locator for circles Let C = { C i , C / v } be the set of generators or sites, with all the Cj being circles in R 2 . Let pi be the centre of Cj and rj be the radius of Cj . The definitions of bisector, influence zone, Voronoi region and Voronoi dia-gram presented in Chapter 1 generalise to the case where the set of sites S is a set of circles C, and the distance d ( M , Cj) between a point M and a site Cj is the Euclidean distance between M and the closest point on Cj from M , i.e. 53 Figure 3.1.11: Seven Apollonius circles centres that are true Voronoi vertices d(M, Ci) = \8 (M,pi) — rj | , where 6 is the Euclidean distance between points. Ob-serve that assuming Q is centred on pi and r» = u>i for i = 1, ..,N, this distance is the absolute value of the additive distance used in the previous subsection. The Voronoi region of Ci with respect to the set C is thus defined by: V (CitC) = {M G M 2 | V j ^ i : \S (M,Pi) - n\ < \S (M,Pj) - rj\}. The Voronoi dia-gram of C is defined by: V (C) = [}c.&c dV {ChC). In the previous subsection, we observed that two Apollonius circles centres are true Voronoi vertices of the Additively Weighted Voronoi diagram (the circles that are either externally or internally tangent to three given circles). When the sites are circles, up to seven of the eight Apollonius circles may be relevant to the Delaunay graph conflict locator (see Figure 3.1.11). We consider the maintenance of the Delaunay graph of circles in an incre-mental way: we check the validity of all the triangles of the Delaunay graph whose vertices are a given triple of circles with respect to a given newly inserted circle. Thus, four circles C\, C%, Cz and C 4 are given: the first three are supposed to define one or more triangles in the Delaunay graph, and the last one is the newly inserted 54 circle. Let (xi,yi) be the coordinates of pi for i = 1,2,3,4. There are two possible outcomes to the above test of validity. Either the triangles are valid with respect to the newly inserted weighted point and the triangles remain in the new Delaunay graph, or there is at least one triangle that is not valid with respect to the newly in-serted weighted point and these triangles will not be present in the Delaunay graph any longer. The Apollonius circles of C\, C2 and C 3 can be obtained algebraically by computing the common intersection of the three circles C[, C2 and C'3 (see Figure 3.1.7) expanding or shrinking from the three first circles C\, C2 and C 3 all with the same absolute value of the rate. The common unsigned expansion of the first three circles is denoted by r. The coordinates of the intersection I of C[, C2 and C3 are denoted (x, y). The circle C" centred on (x, y) and of radius r is tangent to the first three circles. Thus, the Apollonius circles are the solutions of one of the eight following systems (I) of three quadratic equations in three unknowns x,y,r: ( ( x _ 2 ; i ) 2 + ( 2 / _ J / 1 ) 2 _ ( ; i ± r ) 2 = 0 < (x - x2)2 + {y - y2)2 - (r 2 ± rf = 0 • ^ (x - xzf + {y- y3)2 - (r 3 ± rf = 0 By replacing r by — r in one of the preceding systems of equations, we still get another one of the preceding systems of equations. Thus, let us suppose r is the signed expansion of C\. Then, we can reformulate the preceding systems of equations as the following systems (II) of equations: f (x - xt)2 + (y- yx)2 - {n + r)2 = 0 < {x - x2)2 + (y- V2? ~ (r2 ± rf = 0 ^ (x - x3f + (y- yzf - (r3 ± rf = 0 Now let us consider for each system (II) the set X of solutions of the system (II) of equations in K3, where K is an algebraically closed field. Subtracting one of the equations from the remaining two results in a system 55 of 2 linear equations, from which x and y may be expressed as linear functions of r. Substitution in the first equation then leads to a quadratic equation in r. This means that the unknown quantities x,y,r can be expressed with quadratic radicals as functions of the given centres and radii for each one of the systems of equations above. As before, though the simplest thing to do now would be to compute the two Voronoi vertices and use their computed coordinates and corresponding signed ex-pansion in the computation of the values certifying the output of the Delaunay graph conflict locator, it is not desirable because this method would not be generalisable to conies or higher degree curves. For the Delaunay graph of additively weighted points, the true Voronoi ver-tices are the solutions of one system of algebraic equations. Unlike the previous case, for the Delaunay graph of circles, the true Voronoi vertices are not all the solutions of one system of algebraic equations, but a subset of the solutions of four systems of algebraic equations. The solutions of the algebraic equations are the Apollonius circles, whose centres are generalised Voronoi vertices (a concept that we will intro-duce in Section 5.2). We thus need to determine which Apollonius circles centres are potentially true Voronoi vertices (only the real Apollonius circles centres can be true Voronoi vertices). There are four possible determinations of the true Voronoi vertices from Apollonius circles centres of C\, C2 and C3: first case if C\, C2 and C 3 mutually intersect, then the real circles among the seven Apollonius circles that are not internally tangent to each of Ci , C2 and C 3 correspond to true Voronoi vertices (their centres are true Voronoi vertices, see Figure 3.1.11), and reciprocally. second case if one circle (say C\) intersects the two others (C2 and C 3 ) which do not intersect, then only the real Apollonius circles that are either externally tangent to each of C\, C 2 and C 3 , or internally tangent to C\ and externally 56 Figure 3.1.12: Four Apollonius circles centres that are true Voronoi vertices tangent to C2 and C 3 correspond to true Voronoi vertices (their centres are true Voronoi vertices, see Figure 3.1.12). t h i r d case if two circles (say C\ and C2) intersect the interior of the third one (C3) and at least one of them (say C\) is contained in the interior of C 3 , then only the real Apollonius circles that are externally tangent to C\ and C2 and internally tangent to C 3 correspond to true Voronoi vertices (their centres are true Voronoi vertices, see Figure 3.1.13). four th case otherwise (if none of the three situations above apply), only the real Apollonius circles that are externally tangent to C\, C2 and C 3 correspond to true Voronoi vertices (their centres are true Voronoi vertices, see Figure 3.1.14). The case where one circle (say C\) lies in the interior of a second circle (say C2), which lies in the interior of the third circle (C3), or only one circle (say C{) lies within the interior of one of the other ones (say C2) cannot happen because then, there would be no Voronoi vertices and the triangle C\C2C3 would not exist in the 57 Figure 3.1.13: Two Apollonius circles centres that are true Voronoi vertices (first case) Delaunay graph. Now that we have seen the different cases of true Voronoi vertices, we will see how we can test in which case we are and which solutions of the systems of equations (II) described above correspond to true Voronoi vertices. first case C\, C2 and C 3 mutually intersect if, and only if, d{p\,p2) — r\ — r 2 < 0 and d(pi,ps) — r\ — 7-3 < 0 and d(p2,ps) — T2 — ^3 < 0. The computation of this test can be done exactly, since the only variables that are not input to the Delaunay graph conflict locator are the distances, and these distances are expressed by radicals. Indeed, we need to test the sign of the difference of a radical and a number which do not depend on intermediary computations. The true Voronoi vertices are the real solutions of all the systems of equations (II) such that r > 0. second case C\ intersects C2 and C 3 , and C2 and C 3 have no point of intersection if, and only if, d (pi,P2)-n-r2 < 0 and d (p1,p3) - r : - r 3 < 0 and d (p2,P3) -T2 — f2, > 0. The computation of this test can be done exactly for the same reasons as the previous case. The true Voronoi vertices are the real solutions of the system of equations: with r < 0. third case C\ lies in the interior of C 3 and C 2 intersects the interior of C 3 if, and . only if, d(pi,pz) + r\ - r 3 < 0 and d(p2,P3) - r 2 - r 3 < 0 and (x\ - x 3 ) 2 + the same reasons as the previous case. The true Voronoi vertices are the real solutions of the system of equations: (yi—yz)2 — 7-3 < 0. The computation of this test can be done exactly for 59 I (x - Xl)2 + (y - y i ) 2 - (r-i + r)2 = 0 {x - x2)2 + (y - y2)2 - (r 2 + r)2 = 0 (x - x3)2 + {y- y3)2 - (r 3 - r ) 2 = 0 such that r > 0. fourth case this is the case if all the previous three tests failed. The true Voronoi vertices are the real solutions of the system of equations: f (x - xx) 2 + (y - yx)2 - (n + r ) 2 = 0 < (x - x 2 ) 2 + (y - y2)2 - (r 2 + r ) 2 = 0 .. ^ (x - x 3 ) 2 + (y - y3)2 - (r 3 + r ) 2 = 0 with r > 0. As before, we used the same algebraic machinery to compute the values of polynomials that are taken by the true Voronoi vertices without solving any intermediate system of equations. We computed the Grobner basis of the ideal of X for each one of the systems (II) encountered. Each one of these Grobner basis consists of the earlier mentioned quadratic equation in r and linear equations in x, y and r. For the Delaunay graph of additively weighted points, we observed that eval-uating the signs of a single polynomial (g = ( X 4 — x ) 2 + (3/4 — y)2 — (r + r 4 ) 2 ) taken on the real points of X was enough to provide the values certifying the presence in the list output by the conflict locator. As before, we can check for the existence of real solutions by evaluating the sign of the discriminant of the characteristic polyno-mial. We will suppose the real solutions to the systems (II) have been tested. Unlike in the previous case, here we need to evaluate the signs taken by both g and r on each one of the points of X. Indeed, we need not only to check the relative position of C4 with respect to the Apollonius circles, but we need for each Apollonius circle, to check the relative position of C4 with respect to that Apollonius circle, and to check whether that Apollonius circle corresponds to a true Voronoi vertex. As before, we considered the operation of multiplication of polynomials by 60 the polynomial g, whose sign expresses the relative position of C 4 with respect to C". We also considered the operation of multiplication of polynomials by the polynomial r, whose sign allows one to check whether the solutions correspond to true Voronoi vertices. These operations are linear mappings. The operations of these mappings on the canonic representative of the remainder of the Euclidean division of a polynomial by the three polynomials of the system are also linear mappings that can be expressed by a matrix. We need to be able to associate the signs of the values of g with the signs of the values of r taken on the (real) solutions of each system (II). For a given system (II), let Mg and Mr be the matrices of the result of the multiplication by g and by r respectively on the canonic representative of the remainder of the division of a polynomial by the three polynomials of the system. Since these multiplication maps commute, it is possible to use the transformation matrix obtained during the computation of the Jordan form of one of these matrices to triangularise the other matrix by a simple multiplication of matrices [CL098]. Indeed, the computation of the Jordan form for Mg gives the triangular matrix P~lMgP of the Schur form of that matrix where P is a unitary matrix called the transformation matrix; and P~lMTP is triangular. Finally, we can obtain the solutions by reading the diagonal entries in turn in each one of the Jordan forms of these matrices (the diagonal entries of the Jordan form of a matrix are its eigenvalues). The row number on each one of these matrices corresponds to the index of the solution. B y evaluating the signs of the diagonal entries in the Jordan forms of Mg and of Mr on the same line, we associate the signs of the values of g with the signs of the values of r taken on the solutions of each system (II). We will see that though this method for computing several polynomials in the quotient algebra can be generalised to algebras of higher dimension, it will not be possible in practice to use this method to compute the Delaunay graph for conies. 61 3.2 Formulating the Delaunay graph conflict locator for curves from those curves In this section, we describe different experiences with different Computer Algebra Systems and the interval analysis based solver A L I A S in trying to compute the Delaunay graph conflict locator for algebraic curves. To compute that conflict locator, we need first to write the algebraic equa-tions (and inequalities) that should be satisfied by all the geometric loci that will intervene in the evaluation of the conflict locator. The input of the Delaunay graph conflict locator is a quadruple ( C i , C 2 , C 3 , C 4 ) of semi-algebraic one-dimensional sets in the plane. These geometric loci that are considered in the conflict locator are one point (xi,yi) on each one of the four curves Ci, i = 1,.., 4, and the circle (whose centre will be denoted (x,y) and radius r) touching C i , C 2 and C3 (see Figure 3.2.1). The distance between (x,y) and (£4,2/4) will be denoted R. The formulation of the Delaunay graph conflict locator will be described in more detail in Sections 4.1 (generalised offset), 5.3 (conflict locator for conies) and 6.1 (conflict locator for semi-algebraic sets). For each one of the four curves Ci,i = 1,..,4, there are three equations: the implicit equation of the curve Cf. a (xi,yi), the equation of the nor-mal to the curve Cj at the point (XJ,2/j): nj (x{,yi,xty), and the equation expressing the distance between the point (x'i,yi) and the point (x,y): di(xi,yi,x,y,r) for i = 1,2,3 and d<± (£4,2/4,x,y, R). The equality of the distances between (x,y) arid (#i,2/i), between (x,y) and (x 2,2/ 2), and between (x,y) and (£3,2/3) and r expresses the fact that the circle centred at (x,y) of radius r is tangent to the three curves C\, C 2 , and C3. The equation d\ expresses that the distance between (x,y) and (£4,2/4) is R. The equation of the normal is half of the differential of the square distance between the point (xi,2/i) and the point (x,y). It vanishes at the local extrema of that square distance and expresses the necessary condition for (xi, y{) to be a closest point on Ci from the point (x,y). These equations form a system of 12 equations 62 Figure 3.2.1: The Delaunay graph conflict locator that certifies whether the addition of a curve C4 changes or does not change each one of the triangles induced by three curves C\, C2, and C 3 in 12 unknowns. As in Section 3.1, we need to write programs for solving the corresponding system of algebraic equations. Such systems have a finite number of solutions (the solutions constitute a finite set of points i.e. a zero-dimensional variety). As be-fore we consider the quotient algebra of the ring of polynomials in the variables Of the system (x\,XN) by the zero-dimensional ideal I generated by the algebraic equations of the system to be solved. Computing in the quotient algebra involves generalising Euclid's algorithm for the division of one polynomial in one variable by another polynomial in one variable to the division of one polynomial / in the vari-ables XI,...,XN by the generators fi,...,fs of I and having a canonic representative for the coset (i.e. the equivalence class) of / . This can be done by computing a 63 Grobner basis of the ideal I, and then dividing the polynomial / by the generators of the Grobner basis. If the system has a finite number of solutions (points), the quotient algebra dimension is finite and therefore, its bases are finite. We can define a multiplication map mg in this quotient algebra that maps the coset of p to the coset of gp. Let Mg be the matrix of this linear map in the basis B of the quotient algebra. The key property of multiplication operators for evaluating a polynomial g at the points of a zero-dimensional algebraic variety is that the eigenvalues of Mg are the results of the evaluation of the polynomial g on the points of the zero-dimensional variety (see Section 3.3.1). As before we can compute the matrices MXi of the multiplication maps by each one of the variables (x\, . . . , X N ) . Since these multiplication maps commute, it is possible to use the transformation matrix obtained during the computation of the Jordan form of one of these matrices to triangularise other matrices by a simple multiplication of matrices [CL098] (see Section 3.1). Grobner bases can be computed in most Computer Algebra systems. Singular [GPS01] gave incorrect normal forms (i.e. the remainder of the Euclidean division of a polynomial by the polynomials of the Grobner basis), while CoCoA [CNROO] did not give any result within a month of computation. The tests with G B / R S [Fau, Rou] have been done on the L E O N machines of the U M S Medicis [CNR] at Ecole Polytechnique (1 proc alpha 500 Mhz, 640 M o R A M ) from U M S Medicis [CNR]. The tests with Macaulay 2 have been done on the L E O N machines also. The tests with Maple [CGGL92] have been done on a Pentium III 550 Mhz, 1024 Mo R A M . The tests with Maxima [GG82] have been done on an Ultra 5-10 station 440 Mhz, 768 Mo R A M , because the affme module (used for Grobner basis computations) was not available from U M S Medicis machines. The tests with Maxima were done on G I U L I A machines (2 proc P H I 933 Mhz, 1 Go R A M ) from U M S Medicis. The tests with toric resultants involve a preprocessing on Maple in order to produce some files 64 necessary for calling Emiris's [Emi97] "Resin" C program of computation of toric resultants, and a post-processing with Maxima to compute the determinant of the matrix returned by Resin where the symbolic coefficients have been replaced by their formal values. The time corresponding to the Toric Resultant row includes all these three steps. The precomputations with Maple for "Resin" were done on the same machine as the one used for the earlier mentioned tests on Maple. Similarly, the computations with Maxima after "Resin" were done on the same machine as the one used for the earlier mentioned tests on Maxima. The column titled "Canonic conic" corresponds to an equation of the conic of the form ax2 + by2 + c = 0. The column titled "Generic conic" corresponds to a generic equation of a conic of the form ax2 + bxy + cy2 + dx + ey + f — 0. The column titled "+ monomial" corresponds to a preprocessing step on the system to remove the leading monomial by linear combination of the current polynomial with the other polynomials. In the case of Grobner basis, new variables have been introduced and these variables correspond to invariants. Invariants can be useful for rewriting the polynomials as polynomials with fewer monomials and lower degree to make them more manageable in Computer Algebra Systems. Particularly the following property of the Voronoi diagram is useful: the image of the Voronoi diagram of a set of curves by a motion is the Voronoi diagram of the images of the original curves by that motion. Even a partial rewriting of polynomials in the variables x\, ...,x^ as polynomials in x\, ...,x^ and some invariants might simplify the polynomials. The invariants that are relevant here are the fundamental invariants of the group of motions. These are the invariants of the special orthogonal group for points (scalar products and vector products of the vectors involved in the problem: the vectors between points on the original curves and the centre of the circle touching the first three curves). The following tests deal with the Delaunay graph conflict locator for an input constituted by a circle, a parabola, a hyperbola and a circle defined with numeric parameters. The results are summarised in Table 3.1. When there is no numerical 65 without mon + monomial G B / R S 12 hours segmentation fault Toric Resultant "not enough memory" "not enough memory" Maxima bad normal forms bad normal forms Bezoutian-Maple "object too large" "object too large" Macaulay 2 "not zero-dimensional ideal" "not zero-dimensional ideal" A L I A S 11 min. 11 min. Table 3.1: The comparison of different algebraic methods for computing the Delaunay graph conflict locator for a circle, a parabola, a hyperbola and a circle defined with numeric coefficients value, the table reports either the error reported by the Computer Algebra System between quotation marks, or the problem that induced a wrong computation (bad normal forms). The normal forms are canonic representatives of the equivalence class of the reminder of the division of one polynomial by the polynomials of an ideal. We can observe that using invariants to simplify the writing of the polyno-mials before the computation does not bring any acceleration in the computation. We can also observe that approaches based on Grobner bases induce much slower computations ( G B / R S [Fau, Rou]) than interval analysis based methods (ALIAS) . Finally, projective resultants based approaches do not induce any result in the case of Maple [CGGL92]. 3.3 An hybrid approach linking symbolic computation and scientific computation These first experiments with Computer Algebra Systems forced us to distance ourselves from a pure algebraic approach oriented towards Grobner bases and to adopt a hybrid symbolic/scientific computing approach to try to find the optimum running time combination with the guarantee of a certified Delaunay graph con-flict locator. This optimum corresponds to a tradeoff point between what can be 66 achieved with algebraic precomputations (possibly giving rise to huge matrices), and what can be achieved with numerical computations. We will present in this section the key characteristics of the main tools that have been used to find this optimum, and to obtain the Delaunay graph conflict locator for semi-algebraic sets. 3.3.1 T h e sparse resultant In this subsection, we review the main concepts and results about mixed subdivi-sions, mixed volumes, and sparse resultants that will be used in Section 4.4 and Chapter 5. We put the emphasis on the properties of the sparse resultant that have allowed us to get the results presented in Section 4.4 and Chapter 5. Indeed, the sparse resultant of the equations presented in the last section cannot be computed because of a lack of memory (see Table 3.1) even on machines with 4Gb of R A M and more than 6Gb of virtual memory [CNR]. We could compute the sparse resultant needed for the Delaunay graph conflict locator only after simplifying the polynomi-als involved to reduce the sparse complexity of the system of equations. For this purpose, we have used a key property of the Newton polytope and the mixed volume as well as a specific usage of a sparse resultant. The key property of the mixed volume for our work with sparse resultants is that if we replace a polytope (say Qi) by a polytope Q\ whose Newton polytope is strictly included in the Newton polytope of Qi, the mixed volume of Qi, • • • , Q N decreases. We will present this property later in this section. This key property allowed us to get the equation of the generalised offset to a conic (see Section 4.4) and the sparse resultant matrix needed for the algebraic computation of the conflict locator (see Chapter 5). Let us now introduce the concepts of Newton polytope and mixed volume and the key property mentioned in this paragraph. The classical projective resultant of N + 1 polynomials in N + k affine vari-ables is a polynomial in k variables, which characterises the solvability of a system 67 [CEOO, CL098]. Thus, it allows the elimination of N variables, and is therefore also called eliminant. The degree (i.e., the algebraic complexity) of the projective result-ant of several polynomials relatively to the coefficients of one of these polynomials is the product of the degrees of the other polynomials. Sparse resultants generalise the classical (projective) resultant and exploit the monomial structure of the polynomi-als of the system [CEOO]. The Newton polytope expresses the monomial structure of a polynomial, i.e. the monomials appearing in the polynomial. The degree of the sparse resultant is determined only by information about the exponent vectors of the polynomials and it is generally lower than the complexity of the classical projective resultant [CL098]. The sparse resultant has been extensively treated in Canny and Emiris [CEOO, EC95]. Let K be an algebraically closed field (see [Lan02] for a definition of "algeb-raically closed field"). Let K[xi,XJV] denote the ring of polynomials in the vari-ables x\, . . . , X N over the field K. Let K [xf1, ...yxff] = K [ x , £ _ 1 ] denote the field of Laurent polynomials in the variables xi,...,xpj over K. If a= (au...,aN) G ZN, then let xa denote the monomial If Q is a polytope of RN, then let Vol (Q) and Volpj (Q)denote the ./V—dimensional volume of Q. Def in i t i on 3.3.1. (Support of a polynomial [CEOO, Definition 3.1, page 420]) The support Ai of a polynomial fi^Klxf1, . . . j X ^ 1 ] is the set of exponent vectors in ZN corresponding to non-zero coefficients, i.e. fi = YlaeAi cax°'> c<* ^ ^ n e Newton polytope Qi of fi in RN is the convex hull of A{. For example, for the strophoid of equation y2 — x2 — x 3 = 0, the exponent vectors are: (0,2), (2,0), and (3,0). Def in i t ion 3.3.2. (Minkowski sum [EC95, Definition 3.2, page 121]) The Minkowski sum A+B of point sets A and B in is the point set A+B = {a + b\a € A, b G B} (see example on Figure 3.3.1). 68 Figure 3.3.1: The Minkowski sum of the point sets A and B Definition 3.3.3. (Mixed volume [CEOO, Definition 3.4, page 421]) Let N polytopes QI,...,QN C M.N whose vertices belong to Z N . The mixed volume MV (Qi,Qyv) is the coefficient of the monomial in Vol (X\Qi + ... + XNQN)-Definition 3.3.4. (Polyhedral subdivision [CEOO, Definition 4.3, page 421]) A poly-hedral subdivision of a compact set S C M.N is a collection of polyhedra whose union equals S, such that each intersection of two polyhedra of the same dimension is an-other polyhedron in the subdivision of lower dimension. The polyhedra of maximal dimension are called maximal cells or facets. Definition 3.3.5. (Mixed subdivision [CL098, Definition 6.5, page 344]) Let Q = Qi + ... + Qm C be a Minkowski sum of polytopes, and assume that Q has dimension N. Then a subdivision Ri,...,Rs of Q is a mixed subdivision if each cell Ri can be written as a Minkowski sum R{ = F\ + ... + Fm where each Fi is a face of Qi and N = dim (F\) + ... + dim (Fm) (see example on Figure 3.3.2). Definition 3.3.6. (Mixed cell [CL098, Definition 6.6]) Suppose that R = Fi + ...+ Fm is a cell in a mixed subdivision of Q = Q\ + ... + Qm. Then R is called a mixed cell if dim (Fi) < 1 for all i. T h e o r e m 3.3.7. ([CL098, Theorem 6.7]) Given polytopes QI,...,QN C RN and a mixed subdivision of Q = Q\ + ... + QN, the mixed volume MVJv (Qi, ...,QN) *5 computed by the formula MVjv (Qi, ...,QN) = 'ERVOIN (R), where the sum is over all mixed cells R of the mixed subdivision. 69 Figure 3.3.2: A mixed subdivision of the Minkowski sum shown in Figure 3.3.1 and the corresponding mixed volume (the mixed cells are not filled and the filled cells are copies of the Newton polytopes) P r o p o s i t i o n 3.3.8. (Bernstein's Theorem [EC95, Theorem 3.6, page 122]) For the system pi, ...,PN 6 K [xf1,x^1] = K [ x , x _ 1 ] , the sum of the multiplicities (see Definition 4-1-5) of solutions in (K*)N is infinite, or bounded by the mixed volume of the Newton polytopes. If the coefficients are generic, then this bound is exact. Generally, this bound (also known in the literature as the B K K bound be-cause of the existence of closely related papers by Kushnirenko and Khovanskii) is lower than the Bezout number (i.e., the product of the degrees of the polynomials), which expresses the sum of multiplicities of isolated solutions in PN. Now let us introduce the sparse resultant and its use for the computation of the implicit equation of the generalised offset to a conic (see Section 4.4) and the multiplication operator matrix needed for the algebraic computation of the conflict locator (see Chapter 5). The main interest of the sparse resultant in the computation of the implicit equation of the generalised offset to a conic is that it allows one to do the elimination of N variables from N + 1 polynomials like projective multi-polynomial resultants with an algebraic complexity that is lower than the Bezout number, and therefore, the matrices obtained are smaller, and the computations can be performed faster. In order to introduce the sparse resultant, let us recall some basic definitions: Def in i t ion 3.3.9. (Zariski topology, adapted from [GP02, Definition A.2.1, para-70 graph following Lemma A.2.4 and Lemma A.4.3]) In the Zariski topology, a closed set is a subset consisting of all common zeroes of finitely many polynomials with coefficients in K. A quasi-projective variety is an open subset of a closed projective set. The Zariski closure of a set X is the smallest closed set containing X. The sparse resultant is a necessary condition of existence of solutions of a system of algebraic equations in (C*)^: Given a set of exponent vectors A={ct\,a/}cZ^, let C (A) — {c\xai + ... + cixai : Cj € C} be the set of polynomials whose terms all have exponents in A. Consider N + l Laurent polynomials fi£C(A{). Let Qi = Convex hull ( A ) and MV^MV (Q0,..., Q j _ i , Q i + U ..., QN). Let Z (Ao,AN)C£ ( ^ 4 O ) X . . . X £ (AN) be the Zariski closure of the set of all (/o, - , / j v ) for which / 0 (xi, ...,xN)=...-fN (xi, ...,xN) has a solution in (C*)^ . Definition 3.3.10. (Sparse resultant, [CL098, Theorem 6.2 p. 342]) Assume that Qi is an N-dimensional polytope for all i. Then, there is an irreducible polyno-mial called the sparse resultant RCSA0,...,AN in the coefficients of the / j such that (/o,-- , / iv) e Z(Ao,...,AN) & ResAot...:AN =0 . However, it is not sufficient generally. On the toric variety constructed from the Minkowski sum of the Newton polytopes of the polynomials of a system of algebraic equations, the sparse resultant is a necessary and sufficient condition of existence of solutions of the system. Let us define first the toric variety constructed from the Minkowski sum of the Newton polytopes of the polynomials of a system of algebraic equations. For this purpose, we recall that denotes the N—dimensional projective space over K, i.e., the set of lines of K N + 1 going through the origin O of the coordinate system of K N + 1 (adapted from [GP02, Definition A.4.1]). Definition 3.3.11. (Toric variety, [CL098, p. 307]) Let A = {mu...,mi} C ZN, and suppose that / j = antmi + ... + autmi, i = 0 , N are N + l Laurent polynomials 71 in L (A). Assume that the convex hull of A has dimension N. Then consider the map </>A : (C*)N ^ P1'1 defined by <f>A ( t u t N ) = (tmi, tm'). Then, the toric variety XA is the Zariski closure of the image of <f)A. The necessary and sufficient condition of existence of solutions of the system in the preceding toric variety is: Theorem 3.3.12. (adapted from [CL098, Theorem 3.4]) ResA{f0, ...,fN) = 0 if, and only if, fo — ... = fiy = 0 has a solution in The sparse resultant is a factor of the determinant of the Newton matrix, which can be computed using a mixed subdivision [CEOO, CL098] . The sparse elimination using the sparse resultant allowed us to obtain an implicit equation of the generalised offset to a conic (see Section 4.4). The main interest of the sparse resultant of the polynomials fo,---,fN hi the variables x\, . . . ,XJV in the algebraic computation of the Delaunay graph is that it allows one to compute the values taken by the polynomial fo at V ( / i , / # ) . The sparse resultant matrix allows one to do computations in the quotient algebra A = K [x\,..., XN] I {fi,//y). We will see this now. We need now to intro-duce some notations for this purpose. Let fo, . . . /JV be N + 1 Laurent polynomials. Let QO,---,QN be their Newton polytopes. Let R be the sparse resultant of the polynomials fo,---fN- Let degj.R denote the total degree of the resultant R in the coefficients of the polynomial / j . Let MV_i — MV (Qo, ...,Qi-i,Qi+i, . . . , Q N ) . Theorem 3.3.13. ([EC95, Theorem 3.10]) The sparse resultant is separately ho-mogeneous in the coefficients Ci of each fi and its degree in these coefficients equals the mixed volume of the other N Newton polytopes, i.e., deg^R = MV-i. From the last theorem, it follows that the total degree of the sparse res-ultant equals the sum of the mixed volumes of N Newton polytopes: deg (R) = YliLo MV-i. Let Q = Qo + ... + QN- For an infinitesimal vector 5 € Q^, we define e = (Q + 5)nZN. 72 T h e o r e m 3.3.14. [CL098, Theorem 6.17 p. 354] For the set e described above, the cosets [x@] for (3 £ i form a basis of the quotient ring C [xf1, . . j X ^ 1 ] / (fi,...,/#). We define p, : i -* \J{Ai : p aij € Ai p € a = Fo + ... + aij + ... + Fyv, dim Fj > 0, Vj > i, and Bi = {p — p (p) : p 6 e, p (p) € Ai}. We decompose the (M 0 0 Moi\ so that the M 1 0 MuJ rows and columns of Moo correspond to elements of Bo. The elements of BQ are in one-to-one correspondence with the integer coordinates points of the mixed cells of any mixed subdivision of f\,f/y. Thus, the cardinality of BQ is the mixed volume M V - o . Let A = C [xf1, . ^x^ 1 ] / ( / i , ...,/#) be the quotient algebra of the ring of polynomials C [xf 1 , . . .x^ 1 ] by the ideal /JV) generated by / I , . . . , / J V -T h e o r e m 3.3.15. [CL098, Theorem 6.21, p. 356] Let fi € L (Ai) be generic Laurent polynomials, and let fo = UQ + u\X\ + ... + UJVXN. Using the basis from Theorem 3.3.14, the matrix Mf0 of the multiplication map mj0 : A —> A which maps [g] to [fog] is the transpose of the Schur complement of the matrix M\\: M = Moo - M o i M f ^ M i o . Thus the size of the matrix of the multiplication map mj0 is the mixed volume of Q\, ...,QN, which is the sparse bound for the number of common zeroes °f /i> -"J/JV (counted with their multiplicities). This is coherent with the fact the dimension of the quotient algebra is the number of common zeroes of / I , . . - , / J V (counted with their multiplicities). The values taken by /o on V (fi, . . . , / iv) are the eigenvalues of the matrix of the multiplication map my 0 : A —> A (see [CL098, Theorem 4.5]) Moreover, since M/ 0 = UQI + uiMXl + ... + UNMXN, where MXi is the matrix of the map of the multiplication by X{, by [CL098, Corollary 4.3, p. 53], M = UQI + u\M\ + ... + U^MJV , where each Mi is obtained as in Theorem 3.3.15. Then Mf0 = MT implies that MXi = Mi . 73 As observed in Section 3.1, we can compute the matrices MXi of the multiplic-ation maps by each one of the variables {x\,..., xpi) by a simultaneous diagonalisation of all these matrices [CL098]. 3.3.2 N u m e r i c a l methods for comput ing exactly the signs of the eigenvalues of large sparse matrices In this subsection, we present numerical methods for computing exactly the signs of the eigenvalues of large sparse matrices used in the numerical computation of the conflict locator that follows the algebraic precomputation of the matrix of the sparse resultant for the Delaunay graph conflict locator (see Section 5.5). The results re-ported in this section can be found in [Hea02], For computing all the eigenvalues of an arbitrary real matrix, the standard approach is to reduce the matrix to the Hessenberg form (lower Hessenberg: a^ - = 0 for i > j + 1, or upper Hessenberg: aij = 0 for i < j — 1), and then apply QR iteration on that Hessenberg matrix [Hea02]. However, for very large sparse matrices, like the matrix of the sparse res-ultant or the multiplication operator matrix for the Delaunay graph conflict locator computation, standard algorithms for reduction to the Hessenberg form are prohib-itively time and memory consuming [Hea02]. The current method of choice for general sparse square matrices of size N is the Arnoldi method [Hea02]. The Arnoldi method for large sparse non-symmetric eigenvalue problems is implemented in ARPACK [LSY98]. It is the basis for the MATLAB [Hea02] function eigs for computing a few (the six having highest modulus by default) to almost all (up to the iV — 1 having the highest modulus) eigenvalues, and eigenvectors of a matrix. The Arnoldi method is a Krylov subspace method, whose main assump-tion is that the input matrix is best considered as a linear operator, with which one can form matrix-vector products. Krylov subspace methods are based on a 74 simple method (known as the power iteration) for computing a single eigenvalue of an N x N matrix A, which multiplies an arbitrary nonzero vector by successively higher powers of the matrix A. Assuming A has a unique eigenvalue X\ of max-imum modulus and v\ is its corresponding eigenvector, power iteration converges to a multiple of v\. Subspace iteration methods (also known as simultaneous iter-ation methods in the literature) use power iteration with several different starting vectors to compute several eigenvalues of a matrix simultaneously. The sequence of vector spaces spanned by the product of these starting vectors by successively higher powers of the matrix A will converge towards the invariant space spanned by the eigenvectors v\,...,vp corresponding to the p largest eigenvalues of A in modulus A i , . . . , A p . Krylov subspace methods provide similarity reduction to the Hessenberg form using only matrix-vector multiplication. The application of the mathematical setup presented above for computing eigenvalues is difficult because the columns of the matrices computed by the Krylov subspace methods converge towards multiples of the dominant eigenvector of A and thus, they become exceedingly ill-conditioned. In order to remedy to this an orthonormalisation (othogonalising the new vector with respect to all the previous ones and normalising it) of the vector used for matrix-vector multiplication is done at each iteration. This algorithm is owed to Arnoldi [Hea02]. Since the Arnoldi method requires at iteration k a matrix-vector multiplication by A plus O (kN) for the orthonormalisation, plus O (A;3) for the com-putation of the eigenvalues, it is run for a few iterations and then restarted with a new starting vector that is relatively rich in components of the desired eigenvectors. A few repetitions of the restarted Arnoldi process produce excellent approximations to the extreme eigenvalues of A. The certified computation of the sign of eigenvalues of sparse matrices can be done by computing tight bounds on the intervals taken by eigenvalues [Kra92]. This method gives an interval bound for the exact eigenvalue without the user providing error estimations. The existence of the true eigenvalue within the computed bounds 75 is "simultaneously proven by the method" [Kra92]. 3.3.3 ALIAS In this subsection, we present the key properties of the A L I A S library that have been used for the computation of the Delaunay graph conflict locator for semi-algebraic sets presented in Chapter 6. Interval analysis is a well-known method for computing bounds of a function, being given bounds on the variables of that function [Han92, MerOl]. The basic mathematical object in interval analysis is the interval instead of the variable. The operators need to be redefined to operate on the intervals instead of the variables. This leads to an interval arithmetic. In the same way, most usual mathematical functions are redefined by an interval equival-ent. Many packages implement basic interval operations [Han92, MerOl]. However, interval analysis has two main drawbacks. The first drawback is that the bounds of a function depend heavily on the way the function is written [MerOl]. This can lead to or be accompanied with an over-estimation of the bounds [MerOl]. The second drawback is that it is difficult to test rapidly the efficiency of an interval analysis based algorithm [MerOl]. The implementation of such an algorithm involves three different levels of software: a basic level where the basic operations of interval arith-metic are performed, an end-user level where the bounds for the function to be evaluated will be computed using the functions of the first level, and an algorithm level where the function evaluation of the second level will be used to solve a prob-lem [MerOl]. A L I A S [MerOO] is a library of algorithms enabling to analyse and solve zero-dimensional systems of equations and inequalities with real coefficients. A L I A S allows "the user to focus on the algorithmic part of the problem, while offering a convenient way to use interval analysis and, furthermore, enabling one to easily change the analytic form of the function that will be evaluated" (from [MerOl]). A L I A S is composed of a C++ library (the kernel of A L I A S ) , a Maple [CGGL92] interface (which enables to produce the code that will solve a system of equations 76 and inequalities directly from M A P L E by using the C + + library), and a parser (which enables to perform the interval evaluation of a function whose analytical formulation is written in a file). The analysis of the zero-dimensional system of equations and inequalities encountered in the Delaunay graph conflict locator computation (see Chapter 6) relies on results about the uniqueness of a root in an interval (Kantorovitch [Kan57] and Moore-Krawczyk [Moo77]). The mathematical setup is as follows [Han92]. A n interval number is a real, closed interval (x, x). Arithmetic rules are replaced by interval arithmetic rules, e.g. let two interval numbers X = (x, x), Y = (y_, y), then X + Y = (x + y,x + y) and X — Y = (x_ — y,x — y). A n interval function is an interval-valued function of one or more interval arguments. A n interval function F is said to be inclusion monotonic if, Xi C Yi for i e [1,N] implies F(Xi,...,XN) C F ( Y i , Y N ) . A fundamental theorem is that any rational interval function evaluated with a fixed sequence of op-erations involving only addition, subtraction, multiplication and division is inclusion monotonic [Han92]. Theorem 3.3.16. Moore theorem [Moo77j Let a system of N equations in N un-knowns: F = {Fi (xi, ...,a;jv) = 0,i £ [l,iV]} each Fi being at least Cl. Let X be an interval vector for {x\,XN}, y a point inside X and Y an arbitrary nonsin-gular real matrix. Define K as K (X) = y - YF (y) + {I - YF' (X)} (X-y). If ' K (X) C X and \\I — YF' (X) || < 1 then there is a unique solution of F in X . This unique solution can be found using the Krawczyk solving method [Moo77]. Theorem 3.3.17. Kantorovitch theorem [Kan57] There exists a unique solution x* of the functional equation P (x) = 0, where P is an operator twice continuously differentiable, which maps the normed space X onto Y, if the following conditions are satisfied: 77 1. There exists a real valued majoring function Q (x) on an interval (ZQ, z') (i.e., ||P (x0) || < Q (z0) and \\P' (x) \\ < Q' (z) if, \\x - x0\\ < z - z0 < z' - z0) for which the relation Q (z) = 0 has real roots z\, z2 (zo < zi < zi < z'); 2. There exists an inverse operator To = —\P' (xo)]_1, B = —[Q' (^o)]_1 > 0; 3. | | r 0 P(x 0 ) | | <BQ(z0); 4. | | r °P" (x) || < SQ" (z) for ||x - soil < z - z0 < z' - zQ. The solution x* is bounded by \\x* — xo\\ < z\ — ZQ and furthermore it is unique in j|cc — Xo|| £ z2~ZQ. The approximations x/y obtained by the Newton method (ir/v+i — Xjv —[P'(x/v)]_1P (xjv)j and its modification (x^+i=x^—[P'(xo)]~lP(x^)) con-verge to x*. A l l the general purpose solving procedures for zero-dimensional systems have been tested in the computation of the Delaunay graph conflict locator for semi-algebraic sets. These general purpose solving procedures are based on a bisection process on one (single bisection), several (mixed bisection) or all the variables (full bisection) using either: • only the equations and inequalities of the system, • the equations and inequalities of the system and the Jacobian of the system (Moore-Krawczyk test for finding "exactly" the solutions), • the equations and inequalities of the system and the Jacobian and Hessian of the system (with Kantorovitch and Moore-Krawczyk tests). The bisection process (single, mixed or full) can be set by setting the "ALIAS/single_ bisection" parameter (to 2, 1 or 0) and in the case of mixed bisection, the number of bisected variables can be set by setting the "ALIAS/mixed_bisection" parameter. The variables that will be bisected will be the ones having the largest interval width. 78 Let x\,..., XJJI be the set of unknowns and let Z\ = I (x\,x\) [x^xj, be the set of m intervals in which we are searching the solutions of the N equations or inequalities F\ (x\,xm) { = , >} 0 , F N O I , xm) { = , >} 0. The indices in the iteration indices. Let Fi be the interval value of Ti when this equation is evaluated for the interval value (xi, x~j") , ( x m , x^) of the unknowns while F (Tj) be the N—dimensional interval vector constituted of the Fi when the unknowns have the interval value defined by the set Tj. The algorithms use a list of interval vectors (or boxes) T whose maximal size M is an input of the program. The general purpose solving algorithm bisects the input intervals until either their width is lower than an accuracy on the variables e or the width of the equation interval is lower than an accuracy on the equations ep (provided there is enough storage space in the list to store the intervals) [MerOO]. Then, if all the equations and inequalities intervals are acceptable (they contain 0 for an equation or they contain positive or negative values according to the sense for an inequality), we get a new solution, if one of them is not acceptable, there is no solution of the system within the current variable intervals. The new interval vectors are added to the list of T with an ordering which aims at considering first the input intervals having the highest probability of containing a solution. There are two ordering criteria: maximum equation ordering and maximum middle point equation ordering. In the maximum equation ordering, the boxes are ordered along the value of C = Max \Fk (l),Fk(l)j for all k in [1,N] (the first box will have lowest C). In the maximum middle-point ordering, the boxes are ordered along the value of C = Max [Fk (Ci), Fk (Ci)) where Ci is the vector whose components are the middle points of the intervals I. This full bisection process on all the variables simultaneously may induce a combinatorial explosion (at each iteration, 2^ new interval vectors or boxes are produced). Instead of bisecting on all the variables simultaneously, it is possible to bisect only one variable at each iteration. This may xf are unknowns indices, while the exponents in the x\ (and indices in the Tj) are 79 reduce the computation time as the number of function evaluations may be reduced [MerOO]. The variable that will be bisected is the one for which the bisection will produce the intervals with the lowest criteria except if for at least one variable the intervals that will result from the bisection cannot possibly contain a solution. Moreover, to avoid bisecting always the same variable, another test is used: let di be the width of the interval [XJ, x7j and dmax be the maximum of all the df if - r ^ — < 0.1, Xi is not considered as a possible bisection variable. "max . Aside from these bisection processes, it is possible to use another bisection method called the 3B approach (by setting the " A L I A S / 3 B " parameter to 1 and the maximal range " A L I A S / M a x 3 B " and minimal range "ALIAS/Del ta3B" para-meters). Each variable Xi and its range [XJ, xi] are considered in turn. Let x™ be the middle point of this range. First, the interval evaluations for the equations and inequalities in the system with the full ranges of the variables except for the variable i where the range is [x,, x™] are computed. Clearly, if one of the equations or inequalities is not satisfied, it is possible to reduce the range of the variable i to [x™,x7j- If this is not the case, let's define a new x™ as the middle point of the interval [XJ, x™] and repeat the process until either we have found an equation or an inequality that is not satisfied (inducing a new reduction of a variable interval) or the width of the interval [XJ, x™] is lower than a given threshold S. A similar pro-cedure can be used to reduce the input interval on the right. We may additionally select a subset of equations and/or inequalities whose intervals will be evaluated. This can be done by setting the "ALIAS/SubEq3B" variable [MerOO]. 80 Chapter 4 The offset to an algebraic curve While the offset (i.e. locus of points at a given distance, see an example on Figure 4.1.1) to a curve or a surface and the bisector of two curves or surfaces have been studied for their applications (see [HV91]), nothing has been written about the degree of the polynomials defining these objects, and no implicit equation has been given, even in the simple case of conies. In the reminder of this thesis, we will call the offset true offset (to contrast it with the generalised offset, that we will review in this chapter). What Hoffmann and Vermeer [HV91] define as offset curves, Arrondo, Sendra and Sendra [ASS99] define as generalised offset curves (see Figure 4.1.2) . The extraneous solutions (corresponding to a singular point of the curve or surface, for example, the circle centred on the self intersection of the strophoid on Figure 4.1.2 and of radius the offset parameter) have been addressed in [HV91]. Hoffmann and Vermeer [HV91] did not address the computations, but they gave some examples computed using Grobner bases (see also [Hof90]). Arrondo, Sendra and Sendra [ASS99] computed the genus of the generalised offset curve when the field has characteristic [Lan02, end of Section 2 of Chapter II] zero. Farouki and Neff studied the analytic properties [FN90b] as well as the algebraic properties [FN90a] of the true offset to a planar parametric curve. In this chapter, we will address the degree of the true offset to an algebraic 81 plane curve in its most general setting, i.e. in an algebraically closed field of zero characteristic. Knowing the degree of the generalised offset to conies allowed us to identify the factor of the sparse resultant that correspond to the implicit equation of the generalised offset to a conic (see Section 4.4). This implicit equation is central to the formalisation and computation of the Delaunay conflict locator for conies (see Chapter 5 and Section 6.3). Moreover, the degree of the generalised offset to conies determines the Bezout bound 1 on the degree of the algebraic variety on which we will evaluate the Delaunay graph conflict locator. Our main contributions are a general formula for the degree of the true offset curve and its application for the determination of an implicit equation of the generalised offset to a conic (which is the Zariski closure of the true offset, see Section 4.1). The conic is defined implicitly by a formal polynomial (i.e. a polynomial whose coefficients are formal constants). We used the general formula for the degree of the true offset curve to eliminate the extraneous factors from the sparse resultant [CEOO, CL098] to get an implicit equation defining the generalised offset to a conic (see Section 4.4). This chapter is organised as follows: in section 4.1, we study the equation of the offset. In section 4.2, we study the algebraic properties of the true offset to an algebraic curve in order to determine its degree. In section 4.3, we apply the results of section 4.2 to the conies. In section 4.4 we use the results of sections 4.3 and 3.3.1 to compute an implicit equation of the generalised offset to a conic. 4.1 Equations defining the offsets Let us first recall some basic definitions about algebraically closed fields, rings of polynomials and algebraic varieties. Let i f be a zero characteristic algebraic-ally closed field, i.e., a field such that all polynomials with coefficients in K have a root in K [Lan02, Definition before Theorem 1, Section 2, Chapter VII] and Vx <E K, N 6 N : Nx ^ 0. Let K [xi,x/v] be the ring of polynomials [GP02, Defin-*In this case, the Bezout number is the degree of the generalised offset to the power 4 82 Figure 4.1.1: The strophoid (C) and its true offset (thick lines) ition 1.1.3] in the variables x\, ...,xpj with coefficients in K. In all this chapter, we assume,that V ( / i , f s ) C E denotes the algebraic variety embedded in E defined by the polynomials fi, . . . , / s , i.e., the set of all'the points of E whose coordinates are common zeroes of all the polynomials fi, ...,fs [Sha94]. If E = KN is an affine space, then the variety is an affine variety [GP02, Definition A.2.1]. If E = P ^ is a projective space over K, then the variety is a projective variety [GP02, Definition A.4.2]. We can now recall the notion of degree. The degree of a projective variety X C P ^ is the maximum number of points of intersection of X with a projective linear subspace p ^ - d i m X Q f c o m p l e m e n t a r y dimension in general position with re-spect to X (see page 234 in [Sha94]). Thus, the degree of a projective variety is the degree of its maximal dimensional component. In this chapter, we focus on algebraic curves in the affine space K2 — C 2 . We thus suppose K = C in the remaining of this thesis. We will now introduce the true offset curve. Definition 4.1.1. (True offset) Let C = V (/) C K2, for / G K[x,y], be an algebraic curve and R € M + be the offset parameter. The true R—offset curve to C 83 o C o 0 c o Figure 4.1.2: The strophoid and its generalised offset is the locus of points of R 2 being at the distance R from C (see Figure 4.1.1). Remark 4.1.2. This is equivalent to saying that each point q = (u, v) of the true offset curve is the centre of a circle T> of radius R that is tangent to C, and does not contain any point of C in its interior. We will suppose that / has positive degree (i.e., / is not constant), which implies that C is not empty and not equal to K2. We will also suppose that R ^ 0 unless stated otherwise. Let us now introduce the generalised offset and emphasize its differences with the true offset. There exists a superset of the true .R-offset curve, called the generalised P-offset curve and denoted by O that is defined as the locus of points that are locally at the distance R from the given curve (see example on Figure 4.1.2): Definition 4.1.3. (Generalised offset, adapted from [ASS99]) The generalised offset to a hypersurface v at distance R is the Zariski closure of the set of intersection points of the spheres with centre on a non-singular point of v and radius R, and the normal lines to v at the centre of the spheres. This is equivalent to saying that each point q = (u, v) of the generalised offset 84 1.5-y c © . . . <Sfiy / o o / / / \ v ^ \ -1 .5-\\\° o \ \ \ c Figure 4.1.3: Relationships between the generalised offset to the strophoid and the strophoid curve O is the centre of a circle V of radius R that is tangent to C, but may contain points of C in its interior (see Figure 4.1.3). We will now establish the systems of equations and inequalities that define the generalised B-offset O and the true .R-offset to an algebraic curve C. For this purpose, we will describe which points have to be removed from the generalised offset in order to get the true offset. Let us first introduce the following notation and definitions. Let q = (u, v) be an arbitrary point on the generalised P-offset (see Figure 4.1.3). If C is the affine variety defined by / € K[x,y] (i.e. C — V (/)), then the normal to C at a given point m = (a, (3) £ C is the variety defined by n (i.e. V (n)), for n (u, v) = -fy (m) • (u - a) + fx (m) • (v - @) where fx = d^gx'y'> and fy = 9^y'y^ denote the partial derivatives of / . Definition 4.1.4. (Tangent space [Sha94]) The tangent space to an algebraic variety V fs) C E at m 6 V (fx,fs) is the locus of points on lines tangent to V(fi,...,fa) at m. 85 Definition 4.1.5. (Singular point) A point m of an algebraic variety V (fx,fs) C E where fi G K[xx,x2, ...,xrf] is called a singular point if, and only if, the tangent space at m is E, or equivalently, ^ (m) = 0 for all i = 1,.., s and j = 1,.., n. We have already seen an example of singular point: the self intersection of the strophoid. We are ready to describe the points that belong to both the generalised offset and the true offset to an algebraic curve. For a given q = {u,v) G R 2 , let M be the set of points m = (a, (5) G C P | R 2 such that q G V (n). This condition is achieved whenever the normal to C at m passes through q, or m is singular. In the general case, M is finite. However, if C is a circle centred on q, then M. — C f )R 2 . To get in all cases a finite set of points m of C P ] R 2 such that q G V (n), we use & = M when M. is finite, and & = {w} for an arbitrary point w of C f ) R 2 when C is a circle centred on g. Lemma 4.1.6. The set of all the closest points on C P)R 2 from q is contained in M. Proof. The polynomial n defining M expresses half of the differential of the scalar product qm • qm (which is equal to the square of the Euclidean distance 5 (q,m)) with respect to m. The closest points on C f] R 2 from q are global minima of the Euclidean distance 5(q,m), so, the differential (and thus, n) vanishes on them. • Recall the definition of the power of a point p, with respect to a circle centred at c of radius r: it is equal to cp • cp — r 2. The power is positive, zero or negative if p is outside, on, or inside the circle respectively. Let V be the circle of radius R centred on q (see Figure 4.1.3). Lemma 4.1.7. The minimum power of the points of C f] R 2 with respect to T> is at least the minimum power of the points of (3 with respect to V. Proof. The points of C f ) R 2 having minimum power with respect to V are the closest points on C Q R 2 from q, because V is centred on q. In Lemma 4.1.6, we 86 have seen that the set of all closest points on C P| R 2 from q is contained in M. Thus, the minimum power of the points of C f] M 2 with respect to V is the power of a point of M. with respect to T>. If C is not a circle centred on q, M = 6 , and we are done. If C is a circle centred on q, all the points of M are at the same distance from q: the radius of C. Thus, the minimum power of the points of C P ] K 2 with respect to T> is the power of the point w ol &. • A direct consequence of Remark 4.1.2 is that the power of any point of Cf]R2 with respect to any circle V centred on a point q of the true i?-offset and of radius R is positive or zero. Lemma 4.1.7 allows us to restate this last condition as the power of any of the points of the sets & with respect to any of the circles V must be positive or zero. B y Definition 4.1.3, the point q on the generalised i?-offset curve O can be constructed from a non-singular point p = (x,y) on C as the intersection of the normal to C at p and the circle centred on p, and of radius R (see Figure 4.1.4). This circle is the variety V (d), where d (x, y, u, v) = (u- x)2 + (v~ y)2 - R2, whereas the normal is the variety V (n), where n(x,y,u,v) = -fy-(u-x) + fx- {v - y). We are ready to write the equations of the offset. Let us consider the map 7r : K6 —> K2 defined by 7r ((x,y,u, v,a,(3)) = (u,v). The generalised i?-offset O is the Zariski closure of the image by 7r of the variety of K6 defined by the following system of equations and inequalities: / (x,y) = n (x,y,u,v) = d (x,y,u, v) = 0 fx(x,y) # 0 or fy(x,y) # 0 The first line in the preceding system of equations and inequalities contains the algebraic equations defining the point p on C and the point q on the generalised 87 Figure 4.1.4: The construction of a point of the generalised offset offset to C. The second line is a necessary and sufficient condition for p not being singular. The true .R-offset is obtained as the difference of the generalised .R-offset O and the union of each one of the images by n of the sets defined by the following system of equations and inequalities for each point m = (a, 3) of 6 : ( f{a,3) = n(a,B,u,v) = 0 ( T J - a) 2 + (v - Bf - B? < 0 We will now ennounce some fundamental properties of the true offset that are central to the design of the Delaunay graph conflict locator. Proposition 4.1.8. The true offset to an algebraic curve is a semi-algebraic set. Proof. The true .R-offset is defined as the difference of two sets. The second set is a finite union of sets (since & is finite), each one of them being the projection of a set defined by a system of equations and inequalities. Such a set is a semi-algebraic set. The projection of a semi-algebraic set is semi-algebraic [BCR98, Thm.2.2.1]. A finite union or difference of semi-algebraic sets is semi-algebraic (see Note following 88 Definition 2.1.4 in [BCR98]). Thus, the second set is a semi-algebraic set. The first set is also semi-algebraic, as the image of an algebraic set by a projection. • Proposition 4.1.9. The true offset to an algebraic curve is not necessarily algeb-raic. Proof. We have proved that the true offset is a difference of two sets, where the first one is the projection of an algebraic set, and the second one is a semi-algebraic set. If the second set is different from the empty set, then the true offset can-not be algebraic. Let's consider the strophoid, which is the affine algebraic variety V (y2 — x2 — a;3) (see Figure 4.1.2). The true offset curve to a strophoid differs from its generalised offset because around the origin of the coordinate system, the gener-alised offset has two branches which intersect the strophoid. Indeed, the intersection points of the strophoid and its generalised offset cannot be part of a true offset (with positive offset value) to the strophoid, because their distance to the strophoid is zero. Thus, any positive true offset to the strophoid is not algebraic. • The fact the true offset is not an algebraic set implies that the true Voronoi vertices (as intersections of 3 true offsets) are not algebraic, and thus, they cannot be defined by a system of algebraic equations. The fact the true offset is a semi-algebraic set implies that the true Voronoi vertices are semi-algebraic sets (as finite intersection of semi-algebraic sets). The algebraic variety closest to the true Voronoi vertex in the plane is the intersections of three generalised offsets, that we will call a generalised Voronoi vertex. This notion will be formally defined in Chapter 5, and it will be used as a central tool for the Delaunay graph conflict locator for conies. 4.2 The degree of the generalised offset curve In this section, we will prove a general formula for the degree of the generalised offset to an algebraic curve. Let us start with some notations. Consider polyno-mials in K[x,y,u,v] and the projective space P 4, in which K4 is embedded. The 89 homogenisation variable will be denoted as t. We consider the point q = (u, v) on the generalised offset curve O constructed from an arbitrary point p = (x, y) on C = V (/) C K2. In this section, the varieties will be considered in different under-lying spaces. We will use the notation V (/) extensively hereafter. When necessary, we will precise the underlying space by an inclusion, e.g. V (/) C K4 or V (/) C K2. Now let us define the different affine varieties that allow one to define the generalised offset. We will consider these affine varieties in K4 since the polynomials involved belong to K[x,y,u,v]. Let B — V (/) C K4 be the image of C by the inclusion map from K2 to K4 defined by (x,y) t-> (x,y,u,v). Let N = V (n), D = V (d) and V(fx,fy) be considered in K4 in the same way. Now, V (n) is a proper subset of K4 provided that p is not singular. If / is not square-free (i.e. / admits a factorisation with square factors), there is an infinity of singular points. The generalised offset curve O is the image of (V (/) n V (n) n V (d)) \ V (fx, fy) by the canonical projection 7r : K4 —> K2 onto the (u, u)-plane. We will determine the degree of the generalised offset considered in K4. Re-call that the definition of degree of a projective variety we gave at the beginning of Section 4.1 depends on the dimension of the projective variety. Thus, we need to determine the dimension of the generalised offset. For this purpose, we will con-sider the projective completion (i.e. the smallest projective variety containing it) of V (/, n, d) C K4. The motivation for the consideration of the projective completion instead of the affine variety lies in the conditions of the theorem (Theorem 4.2.4) we will use in order to determine the dimension of the generalised offset. Then, we will decompose this projective completion as the union of its component at infinity (i.e. the points with homogeneous coordinate equal to zero), its singular component (points induced by singular points p on C), and the generalised offset considered in K4 (i.e. the affine variety V(f,n,d) \ V (fx,fy) C K4). We will thus obtain the generalised offset considered in K4 as a difference of projective varieties, and we will determine their dimensions and degrees in order to determine the degree 90 of the generalised offset. Now, we need some notations related to homogenisation, projective completion and component at infinity. Notation 4.2.1. Let • / denote the homogenisation (i.e. the replacement of each monomial m in / by 7 n £ d e s ^ - d e s ( m ) where t is the homogenisation variable and deg denotes the degree) of the polynomial / , • V denote the projective completion of the variety V, • fT denote the polynomial defining the component at infinity of V (/). The homogenisation of a polynomial / defines the projective completion of the variety V (f): VTj) = V fj)\ see [CL097, Sec.8.4]. Clearly, V (f) , V (n) and V (d) are all subsets of P 4. Let W := Ylj^nV\^nV{dj. Thus, W = V (J,n,d) C P 4 . Let Wa be the subset of W defined as W \ (v ( /„ / , )U7 ( t ) ) C P 4 . Wa is the generalised offset considered in KA. Wa is a quasi-projective variety since V(f) D V (n) n V(d) and V(fx,fy) U V (t) C P 4 are projective varieties. Let us define the projective variety Ws := W n V (fx, fy) and the projective variety Woo '•— W CiV (t). Thus, the generalised offset considered in K4 is the following affine variety: Wa = W\ (W^ (J Ws). We will now determine the dimension of the component at infinity W^ of W. Since W = V (f,n,d) C P 4 and l^oo := W n V (t), W^ = V (fT,nT,dT,t). Thus, in order to determine the dimension of Woo, we will examine how the di-mension is changed when going from V ( / r ) to V (fT,dT), and from V (fT,dT) to V (fT,dT,t), and finally from V (fT,dT,t) to V (fT,dT,t,nT) = W^. For this purpose we use the following theorem on the dimension of an hypersurface. We will now introduce regular functions, hypersurfaces and irreducible closed sets, in order to recall the theorem on the dimension of a hypersurface. Hypersurfaces will also be 91 encountered in the proofsof Lemma 4.2.15 on the dimension of Woo and of Theorem 4.2.18 on the degree of the generalised offset. Def in i t ion 4.2.2. (Regular function on an affine closed set) Let X be a closed set in the affine space KN. A function / given on X is called regular if there exists a polynomial F with coefficients in k such that / (x) = F (x) for all points x 6 X [Sha94], If X is closed in FN and F ^ 0 is a complex valued regular function on X, then we denote by XF the closed subset of X, known as a projective hypersurface, defined by F = 0. Def in i t ion 4.2.3. (Irreducible closed set [Sha94]) A closed set X is'called reducible if there exist closed subsets X\ C X, X2 C X, X\ ^ X, X2 ^ X, such that X = X\ U X2. Otherwise, X is called irreducible. Thus, a closed set X = V (/) is irreducible if, and only if, / cannot be written as / = h\h2 with h\,h2 S K[x,y,u,v] and h± and h2 are not constant polynomials. In such case, we call the polynomial / irreducible. The following Theorem on the dimension of a hypersurface can be found as [Sha94, Thm.4,Ch.l,Sec.6] or [CL097, Cor.4,p.459]. T h e o r e m 4.2.4. (Theorem on the dimension of a hypersurface) If a complex valued regular function F does not vanish on an irreducible projective variety X, then dimXp = dimX — 1. In order to apply the theorem on the dimension of an hypersurface on a pro-jective variety X = V (/) whose irreducibility is not known, we determine if none of the irreducible components of X is contained in an irreducible component of V (F). Indeed, if an irreducible component V (hi) is contained in an irreducible component V (Fi) of V (F), then Fi and F vanish on V (hi), and the dimension of V (hi) f] Xp is the same as the dimension of V (hi). In Lemma 4.2.6, we determine that none of 92 the irreducible components of V (fT) c P 4 is contained in an irreducible component of V (dT) C P 4 . In Lemma 4.2.10, we determine that none of the irreducible com-ponents of V (fT, dT) C P 4 is contained in an irreducible component of V (t) C P 4 . Finally in Lemma 4.2.11, we determine when none of the irreducible components of V (fT, dT, t) C P 4 is contained in an irreducible component of V (nT) C P 4 . Lemma 4.2.15 will conclude the determination of the degree of the one-dimensional component of W^. Here are the notations and definitions needed to ennounce and prove Lemma 4.2.6. Let V (fT) = V (h{) U V (h2) U ••• U V (hs) be a minimal (i.e. such as V (hi) <f_ V (hj) for any j ^ i) decomposition of V (fT) into irreducible closed sets. Such a minimal decomposition of V (dT) into irreducible closed sets contains at most two components since the degree of dT — (u — x)2 + (v — y)2 is 2. V (dT) = V (di) U V (d2), where di := —t (u — x) + (v — y), d2 :— i (u — x) + (v — y) and L is a root of the equation x2 + 1 = 0 in the algebraically closed field K, is a minimal decomposition of V (dT) into irreducible closed sets. Notation 4.2.5. For fi, -. • ,fs € K[xi, ...,XN], let (fi, ...,fs) denote the ideal gen-erated by fi,...,fs in K[xi,XJV]. Let y denote the radical of the ideal [Sha94, CL097] . The radical yfl of an ideal / of a ring A is \fl = {a 6 A\3N e N : aN e./}. We can now state and prove Lemma 4.2.6. Lemma 4.2.6. None of the irreducible components of V (fT) C P 4 is contained in an irreducible component of V (dT) C P 4 . Proof. We will prove it by contradiction. Let us assume that V (hi) C V (dT), i.e., dT E I(V(hi)) C K[x,y,u,v,t] for some i 6 { l ,2 , . . . , s} . B y Hilbert's Nullstel-lensatz [Sha94], I (V (In)) = vW-Since hi is irreducible, y/(hi) = (hi). Thus, there exists g £ K [x,y,u,v,t] such that dT = g • hi. The sum of the degrees of g 93 and of hi equals 2. Now, either the degree of g is 0 and the degree of hi is 2, or the degrees of g and of hi are both 1, or the degree of g is 2 and the degree of hi is 0. The first case is not possible, because if the degree of hi is 2, then hi is reducible. In the second case, dT — g • (ax + By + 7) , where a, 3 and 7 are coefficients of K. This would imply that the terms in v? and the terms in v2 of dT must come from g. These terms induce terms in u2x, u2y, v2x, and v2y which cannot be cancelled. In the last case, hi reduces to a constant, which is impossible since V (hi) / 0. • In order to present Lemmas 4.2.10 and 4.2.11, we need to recall some defin-itions and results about regular functions and regular mappings on projective vari-eties. Definition 4 .2 .7 . (Regular mapping of affine closed sets) Let X be a closed set of KN and Y be a closed set of KN. A mapping / : X —> Y is called regular if there exists N regular functions fx, . . . ,/JV on X such that / (x) = (fx (x), ..,/jv (x)) for all x G X [Sha94]. Let FN denote the N—dimensional projective space, so that a point £ 6 fN is given by N + 1 elements (£Q : ... : £N) of K and not all the & are 0 [Sha94]. Let be the subset of FN consisting of all the points for which & ^ 0. Definition 4 .2 .8 . (Regular mapping of quasi-projective varieties) Let / : X —• Y be a mapping of quasi projective varieties and Y C P ^ . This mapping is called regular if for every point x € X and every open affine set containing the point / (x) there exists a neighbourhood U of x such that / (U) C , and the mapping / : U is regular [Sha94]. Theorem 4 .2 .9 . ([Sha94, Thm.8, Sect.1.6]) Let f : X -+ Y be a regular map between projective varieties with f (X) — Y. Suppose that Y is irreducible and that all the fibres / _ 1 (y) for y G Y are irreducible and of the same dimension, then X is irreducible. 94 We are now ready to state and prove Lemmas 4.2.10 and 4.2.11. L e m m a 4.2.10. None of the irreducible components ofV (fT,dT) C P 4 is contained in an irreducible component of V (t) C P 4 . Proof. Notice that V (fT) and V (hi) are not contained in the hyperplane at infinity V (t) and that the irreducible components of V(fT,dT) are the V(hi,dj) for all i = 1,.. . , s, j = 1,2. Consider the following sequence of projections: V (hi, dj) C P 4 -> V (hi) C P 2 TTij • (t : x : y : u : v) >—> (t : x : y) Clearly, these Wij are regular mappings of projective varieties. Clearly also, each fibre of these mappings is irreducible. The fibres TT^1 (LO) for LO G V (hi) have dimension 1 since the points on these fibres have fixed x, y and t coordinates (those of LO), and their other coordinates u and v are related by the equation of dj. Thus, all the fibres have the same dimension. Then, we can apply Theorem 4.2.9, and conclude that the V(hi,dj) are irreducible. We will show that none of these V (hi,dj) is contained in an irreducible component of V (t). Let's suppose t G / (V (hi, dj)) for some i G {1,2, ...,s} and j G {1,2}. B y Hilbert's Nullstellensatz [Sha94], I(V(hi,dj)) = y/(hi,dj). Since the V (hi, dj) are irreducible, y/ (hi, dj) = (hi,dj). Then there exists a,b G K [x,y,u,v,t\ such that t = ahi + bdj. Since hi and dj don't have monomials with t nor constant terms, t in ahi + bdj must come from a or b or both. Since hi and dj don't have constant terms, the monomial of least total degree containing t in ahi+bdj must have a degree greater than or equal to 1 in the other variables. This is a contradiction. • L e m m a 4.2.11. None of the irreducible components of V (fT,dT,t) C P 4 is con-tained in an irreducible component of V (nT) C P 4 if, and only if, ( / J ) + L {fy) £ (hi) and (fT) - t (fT) i (hi) for i = 1,2,..., s. Proof. Now, consider the following sequence of mappings: V (hi,dj,t) C P 4 - • V (hi) C P 2 TTij • (t : x : y '. u : v) i—• (t : x : y) 95 Notice that V (fT,(F,t) is contained in the hyperplane at infinity V(t), and thus does V (hi,dj,t). Clearly, these mappings are regular mappings of projective varieties and TT^ (V (hi,dj,t)) = V (hi, t). Clearly also, each fibre of these mappings is irreducible, and V (hi,t) is also irreducible. The fibres T T " . 1 (U>) for u> G V (hi,t) have dimension 1 since the points on these fibres have fixed x, y coordinates (those of u), their t coordinate equals 0, and their other coordinate v is related to u by the equation of dj, and is therefore also fixed. Thus, all the fibres have the same dimension. Then, we can apply Theorem 4.2.9, and conclude that the V (hi, dj, t) are irreducible. We will prove the contrapositive of the conditional statements. Let's sup-pose ( / J ) + L (fy) G (hi) for some i G {1, 2 , s } . We have to show that nT — 0 on V (hi,d\,t). Since ( / J ) + i(fy) G (hi), there exists g G K[x,y] such as {fl) +<>{fy) = 9- hi- But, hi = 0 on V (hudx,t). Thus, ( / J ) + L ( / J ) = 0, i.e., fy — i / J . Then, replacing in nT, we get nT = — t / J (u — x) + / J (t> - y) = (—L (u — x) + (v — y)) / J on V (hi,di,t). Since cii = 0 on V (/ij, d\,t), —i (u — x) + (?j — y) = 0, and therefore, nT = 0 on V (hi,d\,t). We can prove in the same way that if ( / J ) - t ( / J ) G (hi) then n r = 0 on V (hu d2). Reciprocally, let us suppose nT = 0 on V (hi,d\,t) for some % G { 1 , 2 , . . . , s}. Since none of nT,hi,d\ depend on t and nT — 0 on V (/ i i , dx, t), nT = 0 on V ( / i i , c i i ) . Since di = 0 on V (hi, d±), then —i (u — x) + (v — y) — 0, and ( n r = - / J ( U - x ) + l / J ( t i - x ) = 0 \ nT = ifyr{v-y) + fZ(v-y) = 0 on V (hi, d\). Since ( — / J + t / J ) = /- ( t / J + / J ) , we can rewrite the last system of equations as f (-fyr+^)(u-x)=0 1 (-/J + /^J) (--?/) = o on V (hi,di). 96 The subvariety V (hi, d\,u — x, v — y) of V (hi,d\) is a one-dimensional variety in the two-dimensional variety V (hi,d\). The open set V (hi,di)\V (hi,d\,u — x, v — y) is a dense open subset of V (hi, d\). From the last system, we know that {—fy + tfT) vanishes on this dense open subset. Now, we consider V (hi, d\, — fT + tf-jT). This is a closed projective set. We know that it contains the open set V (hi,d\) \ V (hi,d\,u — x, v — y). Thus, it contains also its Zariski closure, i.e., V (hi,d\). Thus, (—fT + ifT) vanishes on V (hi,d\). Therefore, there exists a, b e K [x, y,u, v] such that —fT + tfT = ahi + bdi, i.e., -fy + <-/J = ahi + b (-L (u-x) + v- y). Let auv be the sum of all the terms of a containing the variables u or v. Since {—fy + tfx) does not depend on any of the variables u and v, auvhi + b(—iu + v) — 0. Thus, auvhi = b(iu — v). Since the left hand side of this equality has terms in x or y, b must also have terms in x or y. The last two facts imply that there exists a polynomial e € K [x, y, u, v] such that b = ehi and auv = e (tu — v). Since b = ehi, b e (hi). Thus, -fT + ifl e (hi). Thus, / J + t / J = -L ( - / J + tfT) G (hi). In the same way, we can prove that if nT = 0 on V (hi, d2, t) c f 4 then / J — tfT G (hi). • We are now going to analyse the dimension of W^. N o t a t i o n 4.2.12. Let I = 1 if / J + t / J G (hi) or / J - ifT G (hi) for some i G { 1 , 2 , s } , and / = 0 otherwise. L e m m a 4.2.13. The dimension ofW^ is equal to I. Proof. We can see W^x, in the following two equivalent ways We*, = B nl?n D n V(t) C P 4 , or Woo = V (fT, nT, dT, t) C P 4 . Considering the last expression, by Lemmas 4.2.6, 4.2.10 and three repeated applications of Theorem 4.2.4, we get that V (fT,dT,t) C P 4 has dimension 4 — 3 = 1. Thus, the dimension of Woo is 0 or 1. B y Theorem 4.2.4, the dimension of Woo is 0 if, and only if, I = 0 by Lemma 4.2.11. • 97 The following Lemmas 4.2.15 and 4.2.16 give the degrees of the one-dimen-sional component at infinity and of the one dimensional singular component of the projective variety W. They will allow us to conclude with the degree of the generalised offset in Theorem 4.2.18. We recall the definition of localisation (see [GP02, Definition 1.4.4]) at a prime ideal (see [GP02, Definition 1.3.10]) 3^ C K [t\,tm], where i f is a field. We denote by K [x\, ...,xm\^ the set of all rational functions f /g such that f,g& K[xi, ...,xm] where g 0 «p. Def in i t ion 4.2.14. (Intersection multiplicity, adapted from [GP02, Def.A.8.16, p.480]) Let f,ge K[x,y], p = (pi,p 2) G V (f)C]V (g) C K2, and let Mv = (x -Pi,y~P2) be the maximal ideal of p. We define fip (/, g) := dim f c (K[x, V\MVI (/, 9)), and call it the intersection multiplicity of / and g at p. Now let F,G G K[z, x, y] be homogeneous polynomials, let p = (po : p\ : p2) G V(F)f]V(G) C P 2 , and let Mp — (pox — p\z,poy — p2z) be the homogeneous ideal of p. Assume that po ^ 0 then, for the intersection multiplicity \xv (F, G) := dimfc (K[z, x, y]Mp/(F, G)), it is easy to see that it equals pp (/, g), where / = F | z = i and5 = G | 2 = i . If C, D are projective curves such as I (C) — (F), and / (D) = (G), then HP(C,D) := up(F,G) = np (V (F), V (G)) is the intersection multiplicity of C, D at p. Let CQO be the component at infinity of the projective completion of C , i.e. C = V{fT,t)cK*. ' L e m m a 4.2.15. The degree of the one-dimensional component O/WQO is: 21 /v ( F ( / T ) , y ( n T ) ) , peCoo where p' is an arbitrary point ofWoo whose projection on the projective (t, x, y) -plane (or equivalently on Coo) is p. 98 Proof. Lemma 4.2.13 implies that Woo has a one-dimensional component if, and only if, 1 = 1. In this case, there exists an irreducible component V (hi) such that nT = 0 on V (hi,dj,t). Thus, the points of are the same as the ones of ,fTtt)^ji but their multiplicities differ. The one-dimensional component of V ^ (dT, fT'^t)^j is constituted of all the points p' whose projection on the projective (t, x, y)-plane is a point p of C^o C P 2 and whose projection to the affine (u,u)-plane is the circle V (dT). There are also the two isolated cyclic points ( 0 : 0 : 0 : 1 : ± t ) in Woo , but they do not belong to the one-dimensional component of Wx). The degree of the one-dimensional component of V (^J(dT, fT, t)^j equals the product of the degree of dT (which is 2) by the number of isolated points in Coo- The degree of the one-dimensional component of Woo is twice the sum of the multiplicities of V (fT,nT) C P 4 at the points p' for all the points p of Coo- • Let Cs be the affine subvariety of C composed of all its singular points, i.e. Cs = V(fxJyJ)c.K\ L e m m a 4 .2 .16 . The degree of the one-dimensional component o /Ws is: 2 £ M g , ( F (/) , V-(n)) , qecs where q' is an arbitrary point of Ws whose projection on the affine (x,y)-plane (or equivalently on Cs) is q. Proof. Each point q of Cs induces a trivial equation n. Thus, at the level of Ws, it induces a one-dimensional variety that consists of all the points q' whose projection on the affine (x,y)-plane is q and whose projection on the projective (t, u, v)-plane is a projective circle centred at q and of radius R. It follows that Ws does not have a one-dimensional component at infinity. The only component at infinity of Ws are the points (0 : XQ : yo : 1 : yo ± i (1 — XQ)), where (xo, yo) is a common root of fT, fZ(x,y), and (x,y)). 99 The points of Ws are the same as the points of V \J(fx,fyJ,d)j, but their multiplicities differ. The degree of the one-dimensional component of V ^yj(fx, fy, f,d) J equals the product of the degree of d (which is 2) by the number of isolated points in Cs- The degree of the one-dimensional component of Ws is twice the sum of the multiplicities of V (/, n) C P 4 at the points q' for all the points q oi Cs. • Remark 4.2.17. We could think that since the only variables in common to / and n are x and y, and these are the only variables of / , the intersection multiplicity °f V (/) C P 4 and V (n) C P 4 at q' should equal the intersection multiplicity of V (/) C K2 and V (n) C K2 at q. In fact this is not necessarily true. Indeed, by the projection from the four-dimensional projective space to a projective plane, the intersection multiplicity may be lowered. A n example is / = x2y — y3. The intersection multiplicity of V (/) C P 4 and V (n) C P 4 at q' = (1 : 0 : 0 : 0 : 0) is 9 while the intersection multiplicity of V (/) C K2 and V (n) C K2 at q = (0,0) is 6. The degree of Ws is 18, which corresponds to twice the intersection multiplicity of V (/) C P 4 and V (n) C P 4 at The point q' was taken on the 0—generalised offset to the curve defined by / = x2y -y3: it satisfies d = (u - x)2 + (v - y)2 - R2t2 = 0 with R — 0. However, when R = 0, Ws is 0-dimensional, but its degree corresponds to that announced by Lemma 4.2.16. In this case, Woo is 0-dimensional (the three points ( 0 : 1 : 0 : 1 : 0 ) , (0 : 1 : 1 : 1 : 1) and (0 : 1 : - 1 : 1 : -1)) , and I = 0. Finally, the degree of W is 18. However, in this case the degree of the canonical projection TT : P 4 \{ (1 : 0 : 0 : 0 : 0)}—>K2 is not any more 1, but 6. Indeed, any point (x,y) of the generalised offset is the image of possibly 3 double points (t : x : y : u : v) of Wa by TT. This justifies that the degree of the 0—generalised offset is 3. T h e o r e m 4 .2 .18 . The degree of the generalised offset to an algebraic curve V ( / ) C 100 K2 of degree m, such as f is square-free is : 2m 2 - 21 ] T AV (V (fT) > V ("T)) - 2 E / V (/) > ^ (« 0) where p' and q' are defined as in Lemmas 4-2.15 and 4-2.16. Proof. Bezout's Theorem [Sha94] tells us that for projective varieties, the degree of the intersection of two such varieties is generically equal to the product of the degrees of these varieties. Therefore, the degree of W = B tlJ7 D D is generically 2m 2 . The quasi-projective varieties W, Wa, Woo, and Ws are related by Wa = W \ ( W O Q U Ws). Since W is defined by three polynomials, its dimension is at least 1 according to Theorem 4.2.4. Since n is a polynomial in two more variables than / , n cannot identically vanish on B. Since d is a polynomial whose coefficients are polynomials in R, R ^ 0 and the coefficients of / and d do not depend on R, d cannot identically vanish on BriJ7, and the dimension of W is exactly one. Since W is one-dimensional, the degree of W is the sum of the degrees of its one-dimensional components. Thus, the degree of Wa equals the degree of W minus the degree of the one-dimensional component of Woo U Ws. Since Ws has no one-dimensional component at infinity (see proof of Lemma 4.2.16), the one-dimensional components of W^ and of Ws are disjoint, and the degree of the one-dimensional component of W^ U Ws is the sum of the degrees of the one-dimensional components of Woo and of Ws- By Lemmas 4.2.15 and 4.2.16, the degree of Wa is 2m 2 -21 J2 AV (V (f) , V (nT)) - 2 £ ^ (V (/) , V (n)) . (4.2.1) By definition, t never vanishes on Wa. Finally, the generalised offset O is the image of Wa by the canonical projection it : P 4 \ { (1 : 0 : 0 : 0 : 0)} -» K2 : (t : x : y : u : v) ^ The multiplicity of the extraneous variety corresponding to a singular point is at least the product of the valuations of the polynomials / and n at a generic point peCoo (f , f) with degree 1. Thus, the degree of O is (4.2.1). • 101 of Ws. The valuation of / corresponds to the multiplicity mp of the singular point of C. The valuation of n is the valuation of / minus 1 since n involves the partial derivatives of / . Thus, the multiplicity of the extraneous variety corresponding to a singular point is at least mp(mp — 1), and the degree of the corresponding extraneous factor is at least 2mp (mp — 1). 4.3 The degree of the generalised offset to a conic In this section, we use the general formula for the degree of the generalised offset to a conic developed at the preceding section in order to compute the degree of the generalised offset to the different conies. In this section, we deal with real affine curves. P r o p o s i t i o n 4.3.1. The degree of the generalised offset to a circle is 4-Proof. The projective completion of the circle is the projective variety V [[x-atf+ {y-btf-rH2^. Replacing t by 0 in the previous polynomial, we get fT = x2 + y2 = (x + ty) (x — ty) = h\h% By taking the partial derivatives, we get fl = 2x, / J = 2y. Thus, fT + tfT = 2(x + ty) e (hi). Thus, the dimension of the component at infinity of the projective completion of the generalised offset is 1. Since the component at infinity of a circle is the two cyclic points ( 0 : 1 : ±t)» and V (fT) and V (nT) do not meet tangentially, ]CpeCoo Mp' {V (/T) , V ( n T ) ) = 2, and the degree of W^ is 4. Since a circle does not have singular points, Cs = 0- Thus, the degree of the generalised offset to the circle is 2 -2 2 — 4 — 0 = 4. • P r o p o s i t i o n 4.3.2. The degree of the generalised offset to an ellipse or a hyperbola is 8. Proof. The ellipse can be considered in some coordinate system as the affine variety V (j^ + ^ — lj. The projective completion is the projective variety V (b2x2 + o?y2 — a2b2t2). Replacing t by 0 in the previous polynomial, we get 102 fT = b2x2 + a2y2 = hi • h2, where h\ = bx + Lay and h2 = bx — cay. By taking the partial derivatives we get /J = 2b2x, fT = 2a2y. Thus, fi + tfy = 2 {b2x + ta2y) ^ fa), and /J - tfT = 2 (b2x - ta2y) £ fa) unless a2 = b2. Thus, the dimension of the component at infinity of the projective completion of the generalised offset equals 0 unless a = b, in which case the conic is a circle. Thus, I = 0. The ellipse (or the hyperbola) does not have singular points, thus Cs = 0. Thus, the degree of the generalised offset to the ellipse is 2 • 2 2 = 8. The same reasoning allows one to conclude with the same degree for the hyperbola: the difference with the case of the ellipse is that the polynomial ^ + ^ - 1 is replaced by £ - - 1. • Proposition 4.3.3. The degree of the generalised offset to a parabola is 6. Proof. The parabola can be considered in some coordinate system as the affine variety V (y2 — 2px). The projective completion is the projective variety V (y2 — 2pxt). The projective completion of the normal at the point (x,y) is the projective variety V (—2yu + 2xy — 2ptv + 2pty). Replacing t by 0, we get fT — y2 = h\, where h\ — y. By taking the partial derivatives, we get /J = 0, fT — 2y. Thus, /J + tfT = 2iy E fa). Thus, the dimension of the component at infinity of the projective completion of the generalised offset equals 1. The component at infinity is given by fT — y2 — 0, hence ( 0 : 1 : 0 ) . Since V (fT) and V (nT) do not meet tangentially, X^peCoc / V (fT) > ^ (nT)) = 1' a n o - ^ e d e § r e e of is 2. Since a parabola does not have singular points, Cs — 0. Thus, the degree of the generalised offset to the parabola is: 2 • 2 2 — 2 • 1 • 1 — 8 — 2 = 6. • Lastly, we consider the generalised offset curve of two straight lines. The two straight lines are the following affine variety V ((ax + by + c) (dx + ey + /)). It is the union of V (ax + by + c) and of V (dx + ey + / ) . The projective completion of the two straight lines variety is the affine variety V ((ax + by + ct) (dx + ey + ft)). Replacing t by 0 in the previous polynomial, we 103 get fT = (ax + by) (dx + ey) = hi • h2, where h\ — ax + by and h2 = dx + ey. By taking the partial derivatives of the previous polynomial, we get / J = a (dx + ey) + d (ax + by), fT = b (dx + ey) + e (ax + by). Thus, fx + Lfy = (a (dx.+ ey) + d (ax + by)) + i (b (dx + ey) + e (ax + by)) and fx - Lfy=(a (dx + ey) +d (ax + by))-t(b (dx + ey) + e (ax + by)). L e m m a 4.3.4. I = 1 44> ae — bd = 0 or a ± ib — 0 or d ± te = 0. Proof. =>: There must exist a polynomial # of K[x,y] such as (a (dx + ey) + d (ax + by)) ± i(b (dx + ey) + e (ax + by)) = g (ax + by) or (a (dx + ey) + d (ax + by)) ± L (b (dx + ey) + e (ax + by)) = g(dx + ey). Thus, /• * (ax + fey) (d ± ie) + (dx + ey) (a ± ib) = g (ax + 6y) < or (ax + 6y) (d ± te) + (dx + ey) (a ± cb) — g (dx + ey) The first equality implies that ae — bd = 0 or a ± ib — 0. The second equality implies that ae — bd = 0 or d ± te = 0. For conies with real coefficients, this implies that ae — bd = 0. : If ae — bd = 0 then there exists z E k such that dx + ey = z (ax + by). Thus, / J ± ifT = (ax + by) (d ± te) + (dx + ey) (a ± ib) = (ax + by) (d±ie + z(a± ib)), a n d / J ± 6 / J e (M-If a ± i& = 0, / J ± t = (ax + by) (d ± te) + (dx + ey) (a ± ib) = (d ± ie)(ax + by). Thus, / J ± t / J E (hi). If d ± te = 0, then / J ± t / J = (ax + by) (d ± ie) + (dx + ey) (a ± t6) = (a±ib) (dx + ey). Thus, fT±ifT£{h2). In all these cases, 2 = 1. • Proposition 4.3.5. The degree of the generalised offset curve of two distinct straight lines is A. Proof. According to Lemma 4.3.4, I = 1 O ae—bd = 0 o r a ± t 6 = 0ord±ie = 0. How-ever, for real straight lines (i.e. defined by polynomials with real coefficients), the 104 dimension of the component at infinity of the projective completion of the general-ised offset equals 1 if, and only if ae—bd = 0, since a±ib = 0ord±ie = 0 is impossible with a, b, d, e being all real coefficients. Thus, the dimension of the component at infinity of the projective completion of the generalised offset equals 1 if, and only if, the two straight lines are parallel and distinct (if they were not / would not be square-free). The component at infinity of two parallel and distinct straight lines is one-dimensional and it is the zero set of fT = ^(a + d)2 + (b + e) 2^ (ax + by)2, i.e. of (ax + by)2. Thus, since V (fT) and V (nT) do not meet tangentially, Y^peCoo Mp' , y ( n T)) = 2 and the degree of is 4. If the two straight lines are parallel, the variety does not have singular points, and Cs = 0. Thus, the degree of the generalised offset to two parallel and distinct straight lines is 2 • 2 2 — 4 — 0 = 4. Otherwise, the two straight lines have a single real intersection, I = 0, and the variety has exactly one singular point (the intersection point) of multiplicity 2. Thus, since V (f) and V (fi) do not meet tangentially, ]C geCs Mg' {V (/) > ^ (™)) = 2, and the degree of the generalised offset to two non-parallel straight lines is 2-22—0—4 = 4. • The results of this section are summarised in the following Table 4.1. Conic ellipse/hyperbola parabola circle two lines Offset degree 8 6 4 4 Table 4.1: The degree of the generalised offset to conies 4.4 An implicit equation of the generalised offset to a conic In this section, we use the results of the preceding section in order to compute an implicit equation of the generalised offset to a conic as a factor of a sparse resultant. This implicit equation of the generalised offset to a conic will be used in both the algebraic formalisation (in Chapter 5) and the numerical computation (in Section 105 6.3) of the Delaunay graph conflict locator for conies. A conic C can be defined implicitly as the variety defined by a second de-gree formal polynomial: / (x, y) = ax2 + fixy + ^y2 + 5x + ey + £ = 0. We shall compute an implicit equation of its generalised offset. The partial derivat-ives are fx = 2ax + (3y + 5 and fy = (3x + 2jy + e. A n equation of the nor-mal J\f = V (n) C KA to the original conic at the point (x,y) is: n (x,y,u,v) = — (f3x + 2jy + e) (u — x) + (2ax + j3y + 5) (v — y) — 0. If a conic is not degenerate (proper conic different from the union of two lines), then it has no singular points, and CS = WS = 0. The generalised offset to a conic is the Zariski closure of the projection of the affine variety V(f,n,d) \ V (fx,fy) onto the (u,v) plane. Its implicit equation is thus a factor of the resultant expressing the elimination of the variables x and y from the three polynomials f,n,d. However, if the conic is not degenerate, an implicit equation of the generalised offset is the sparse resultant itself. Since we are using the sparse resultant instead of the "normal" projective resultant, we should pay attention to the fact the sparse resultant is a necessary condition of existence of common solutions in (C*)N instead of PN. The axes of equation x = 0 and y = 0 can be part of the generalised offset to a conic only if the conic is the degenerate union of two straight lines with one of them being the line of equation x = ± r or y = ± r , where r is the offset parameter. Since the conic is defined generically, the generalised offset will not generically contain one of the axes of equation x = 0 and y = 0. The main objective that has been sought in the computations is to simplify the polynomials. This simplification has to induce a system of equations equivalent to the original system of equations in order to generate the same set of zeroes (i.e. the same variety). In the case of the sparse resultant computation, this has been 106 Figure 4.4.1: The Newton polytope of / . achieved by replacing the original system of algebraic equations by an equivalent system of algebraic equations, where one polynomial is replaced by a linear com-bination of the polynomials in the system of equations having a Newton polytope strictly inscribed in the Newton polytope of the original polynomial. The sparse resultants have been computed thanks to the sparse resultant software developed by Emiris [EC95, Emi97]. This software allowed us to compute the sparse resultant matrix. We computed the determinant of this matrix and its factorisation. The degree of the generalised offset to conies has been determined in Section 4.3 (see Table 4.1). This allowed us to identify the factor that corresponds to an implicit equation of the generalised offset. In the case where a and B are different from 0, we call the conic "generic". If not stated otherwise, we will suppose the conic is generic in all this subsection. The Newton polytope of the polynomial defining B = V (/) c K4 is illustrated in Figure 4.4.1. The polynomial defining J\f = V (n) C K4 can be rewritten in the following way exhibiting its monomials in the variables x and y, which need to be eliminated in order to get an equation of the generalised offset: n (x,y,u,v) = Bx1 — By2 + 2 (7 — a) xy + (~Bu + e) x + (Bv — 5) y + (—eu + 5v) = 0. The monomial in x2 can be eliminated if we replace n(x,y,u,v) by an (x,y,u,v) — Bf(x,y). The Newton 107 Figure 4.4.2: The Newton polytopes of n and of an — Bf. Figure 4.4.3: The Newton polytopes of d and of / — ad polytopes of n (x, y, u, v) and of an (x, y, u, v) — Bf (x, y) are shown in Figure 4.4.2. It is easy to see from Figures 4.4.1, 4.4.2 and 4.4.3, that no other monomial can be eliminated in addition to x 2 if we replace n by a linear combination of n, / and d. If we replace d(x,y,u,v) by / (x,y) — ad(x,y,u,v), the monomial in x 2 disappears. The Newton polytopes of d (x, y, u, v) and of / (x, y) — ad (x, y, u, v) are shown in Figure 4.4.3. It is easy to see from Figures 4.4.1, 4.4.2 and 4.4.3, that no other monomial can be eliminated in addition to x 2 if we replace d by a linear combination of d, f and n. The mixed volumes of / and an — Bf and of / and f — ad are 4 (see Figures 4.4.4 and 4.4.6), and the mixed volume of an — Bf and f — ad is 3 (see Figure 4.4.5). In our search for an equivalent system of algebraic equations, we could have supposed that n or d will be unchanged instead of / . In both cases, we would arrive to the same mixed volume. We get an equation for the generalised offset as the 108 Figure 4.4.4: The mixed volume of / and an — y •'•"/V MV=3 .'*••- • X Figure 4.4.5: The mixed volume of an — (3f and / Figure 4.4.6: The mixed volume of / and / — 109 sparse resultant of the three polynomials / , an — Bf and / — ad. The code written to get this sparse resultant has been placed in Appendices B . l (Maple code generat-ing the files that will be input to the "Resin" program [Emi97]), and B.2 (Maxima code used to compute the sparse resultant as a factor of the determinant of the mat-rix returned by the "Resin" program). This implicit equation is publicly available at: http:// www.cs.ubc.ca/~fanton. However, it is possible to further simplify the equation of a non-degenerate conic by using a different coordinate system. For an ellipse or an hyperbola, a nice coordinate system has its origin at its centre (the intersection of its axes of symmetry), and axes the axes of symmetry of the conic. For a parabola, a nice coordinate system has its origin at its summit (intersection of the parabola with its axis of symmetry), and one of the axes is the axis of symmetry of the parabola. In the case of an ellipse or an hyperbola, the equation of the conic in a co-ordinate system with origin at the centre of the conic, and axes the axes of symmetry 2 2 of the conic, simplifies to ^ ± |? — 1 = 0, assuming both a and b are different from 0. B y multiplying this equation by a 2 , we get an equivalent equation of the form x 2 ± ^y2 — a? — 0. Then, by replacing ± ^ by c, and - a 2 by e, we get an equivalent equation of the form x 2 + cy2 + e = 0 where e ± 6 2c = 0. Let / = x 2 + cy2 + e. The Newton polytope of the polynomial defining B = V (f) C K4 is illus-trated in Figure 4.4.7. The polynomial defining M — V (n) C K4 can be rewritten in the following way exhibiting its monomials in the variables x and y, which need to be eliminated in order to get an equation of the generalised offset: n (x, y, u, v) = —2cy (u — x) + 2x (v — y) = 0. The Newton polytope of n (x,y,u, v) is shown in Figure 4.4.8. It is easy to see from Figures 4.4.7 and 4.4.9 that it is not possible to remove monomials from n without adding new monomials by replacing n by any linear combination of /, n and d. 110 Figure 4.4.7: The Newton polytope of / for ellipses and hyperbolas. 2 0 y Figure 4.4.8: The Newton polytope of n for ellipses and hyperbolas. i i r Figure 4.4.9: The Newton polytopes of d and of ad - f for ellipses and hyperbolas. y , M V = 4 \ \ \ X Figure 4.4.10: The mixed volume of / and n for ellipses and hyperbolas. If we replace d(x,y,u,v) = (u — x)2 + (v — y)2 — R2 = 0 by f (x,y) — d(x,y,u,v), the monomial in x2 disappears. The Newton polytopes of d(x,y,u,v) and of f(x,y) — d(x,y,u,v) are shown in Figure 4.4.9. As with that equation in the case of a generic conic, no other monomial can be eliminated in addition to x2 if we replace d by a linear combination of d, f and n (see Figures 4.4.7, 4.4.8 and 4.4.9. The mixed volumes of / and n and of / and f — d are 4 (see Figures 4.4.10 and 4.4.12), and the mixed volume of n and / — d is 3 (see Figure 4.4.11). As before, in our search for an equivalent system of algebraic equations, we could have supposed that n or d will be unchanged instead of / . In both cases, we would arrive to greater mixed volumes. The sparse resultant of / , n, / — d in the variables u and v in the case where the equation / of the ellipse or hyperbola is expressed in a coordinate system centred on its centre and with axes its axis of 112 Figure 4.4.11: The mixed volume of n and / - d for ellipses and hyperbolas. Figure 4.4.12: The mixed volume of / and / - d for ellipses and hyperbolas. 113 Figure 4.4.13: A n ellipse and its 0.5—generalised offset symmetry has been placed in Appendix B.4. The degree (8) of this equation matches the degree of the generalised offset to an ellipse or hyperbola obtained in Section 4.3. By a simple change of coordinate system, we can obtain easily the general equation of the generalised offset to an ellipse or an hyperbola. As a means of verification of these symbolic computations, if we replace the formal coefficients of the equation of an ellipse or an hyperbola by their numerical values, we get exactly the same equation of the generalised offset when we use the general equation of the generalised offset to a conic or the general equation of the generalised offset to an ellipse or an hyperbola. Two examples of generalised offsets to an ellipse computed using the preceding equation are shown in Figures 4.4.13 and 4.4.14. Two examples of generalised offsets to a hyperbola computed using the preceding equation are shown in Figures 4.4.15 and 4.4.16. In the case of a parabola, the equation of the conic in a coordinate system with origin at the summit of the parabola, and one of the axes being the axis of the parabola, simplifies to y2 — 2px = 0. The Newton polytope of the polynomial defining B = V (/) c K4 is illus-trated in Figure 4.4.17. The polynomial defining Af = V (n) C K4 can be rewritten in the following way exhibiting its monomials in the variables x and y, which need to be eliminated in 114 Figure 4.4.14: The same ellipse and its 3—generalised offset Figure 4.4.15: A n hyperbola and its 0.5—generalised offset Figure 4.4.16: The same hyperbola and its 3—generalised offset 115 Figure 4.4.17: The Newton polytope of / for parabolas. Figure 4.4.18: The Newton polytope of n for parabolas. order to get an equation of the generalised offset: n(x,y,u,v) = —2yu + 2xy — 2pv + 2py = 0.. This polynomial has monomials in xy, y, and a constant term. It is easy to see from Figures 4.4.17 and 4.4.19 that it is not possible to remove monomials from n without adding new monomials by replacing n by any linear combination of /, n and d. The Newton polytope of n(x,y,u,v) is shown in Figure 4.4.18. If we replace d(x,y,u,v) — (u — x)2 + (v — y)2 — R2 = 0 by f (x,y) — d(x,y,u,v), the monomial in y2 disappears. The Newton polytopes of d(x, y,u,v) and of / (x,y) — d(x,y,u,v) are shown in Figure 4.4.19. It is easy to see from Figures 4.4.17, 4.4.18 and 4.4.19 that no other monomial can be eliminated in addition to y1 if we replace d by a linear combination of d, f and n. The mixed volumes of / and n and of n and / — d are 3 (see Figures 4.4.20 116 Figure 4.4.19: The Newton polytopes of d and of / - d for parabolas. Figure 4.4.20: The mixed volume of / and n for parabolas. and 4.4.21), and the mixed volume of / and / — d is 4 (see Figure 4.4.22). The sparse resultant of / , n, / — d in the variables u and v in the case'where the equation / of the parabola is expressed in a coordinate system centred on its summit, and one axis of the coordinate system is its axis of symmetry has been placed in Appendix B.3. The degree (6) of this equation matches the degree of the generalised offset to a parabola obtained in Section 4.3. By a simple change of coordinate system, we can obtain easily the general equation of the generalised offset to a parabola. As before, as a means of verification of these symbolic computations, if we replace the formal coefficients of the equation of a parabola by their numerical values, we get exactly the same equation of the generalised offset when we use the general equation of the generalised offset to a conic or the general equation of the generalised offset to a parabola. Two examples of generalised offsets to a parabola computed using 117 Figure 4.4.21: The mixed volume of-n and / - d for parabolas. MV=4 \ Figure 4.4.22: The mixed volume of / and / — d for parabolas, the preceding equation are shown in Figures 4.4.23 and 4.4.24. 4.5 Conclusions We have obtained a general formula for the degree of the generalised offset to an algebraic curve defined in its most general setting: an implicit equation with coef-ficients in a zero characteristic algebraically closed field. We have applied it to compute the degree of the generalised offset to conies. We have obtained an implicit equation of the generalised offset to a conic defined by a formal polynomial. We have also obtained simplified equations in two cases: in the first case, the conic is a circle, an ellipse or an hyperbola defined by a formal polynomial; and in the second case, the conic is a parabola defined by a formal polynomial. The same simplific-118 Figure 4.4.23: A parabola and its 0.5—generalised offset ation can be obtained in the case of a degenerate conic (i.e. two straight lines). This implicit equation of the generalised offset to a conic has been used in order to get algebraic descriptions for the generalised Voronoi vertex of three conies and the algebraic Delaunay graph conflict locator for conies. Figure 4.4.24: The same parabola and its 3—generalised offset 119 Chapter 5 The Delaunay graph conflict locator for conies In this chapter, we present an algebraic approach to the computation of the Delaunay graph conflict locator in the case of conies. In Section 5.1, we will show how the definitions related to generalised Voronoi diagrams presented in Chapter 1 adapt to the case where the sites are conies. In Section 5.2, we will present the formalisation of the Delaunay graph conflict locator for conies. In Section 5.3, we present the simplification of the equations defining the Delaunay graph conflict locator for conies. In Section 5.4, we present how these simplications allowed us to compute the matrix of the sparse resultant used for the computation of the Delaunay graph conflict locator. Finally, in Section 5.5, we present the numerical computation of the Delaunay graph conflict locator. 5 . 1 Preliminaries We consider now M = R 2 . Let <S = {C\,Cm} c M, m > 2 be a set of m different conies in M. Let us define the distance d (p, Si) between a point p and a conic Cj as d (p, = inf 6 (p, q) \q 6 S j , where 5 denote the Euclidean distance between points. The definitions of influence zone, bisector, VoronOi region and Voronoi diagram we 120 presented in Chapter 1 hold, and in this context they are the definitions of influence zone, bisector, Voronoi region and Voronoi diagram of conies. 5.2 Formalisation of the Delaunay graph conflict loc-ator In this section, we consider the maintenance of the Delaunay graph for conies in an incremental way: we check the validity of the triangles of the Delaunay graph of the set of sites S formed by a given triple of conies with respect to a newly in-serted conic. The reason of the simultaneous treatment of all the triangles of the Delaunay graph formed by a given triple is algebraic: we cannot treat portions of curves because they are semi-algebraic sets instead of being algebraic sets. We need to treat whole algebraic curves. Thus, four conies Ci, C2, C 3 and C4 are given: the first three are supposed to be the vertices of one or more triangles in the Delaunay graph, and the last one is the newly inserted conic. Our approach is similar to the one we used in [AKM02]. We consider now the notion of generalised offset of a conic (see Section 4.1), which can be seen as an expansion/shrinking of conies. We introduce the notion of generalised Voronoi vertex: Def in i t ion 5.2.1. (generalised Voronoi vertex) A generalised Voronoi vertex of three conies C\, C2, and C 3 is a point of intersection of the generalised offsets of C i , C2, and C 3 with the same offset parameter (see Example on Figure 5.2.1). The generalised Voronoi vertex shown in Figure 5.2.1 is not a true Voronoi vertex because the circle centred on that vertex and tangent to the three conies has points of one of the conies in its interior. In contrast, the generalised Voronoi vertex shown in Figure 5.2.2 is a true Voronoi vertex: the circle centred on that vertex and tangent to the three conies has no points of any of the conies in its interior. 121 Figure 5.2.1: A generalised Voronoi vertex (dot) of three conies (thick lines) Figure 5.2.2: A true Voronoi vertex of three conies (thick lines) 122 Figure 5.2.3: The insertion of C4 induces a conflict with the triangle CiC2C3 The objective of the Delaunay graph conflict locator is to determine whether or not the insertion of a new conic C4 would change each one of the triangles of the Delaunay graph DG (S) of S corresponding to the triple Ci, C2, C 3 . There are two possible outcomes to the conflict locator: either the triangles CiC2C% are valid with respect to C4 and the triangles remain in the new Delaunay graph, or some triangles are not valid with respect to C4 and these triangles will not be present in the Delaunay graph (and in the quad-edge data structure storing it) any longer. We can see an example of the later case in Figure 5.2.3. One of the two triangles C\C2C-i is not valid with respect to the conic C4, thus that triangle will not belong any more to the Delaunay graph (see Figure 5.2.4). The Delaunay graph conflict locator consists of determining which of the true Voronoi vertices of the conies C\, C2, and C 3 are at a distance (denoted as R) with respect to C4 lower than the distance (denoted as r) with respect to C\, C2, and C 3 (see Figure 5.2.5). Let (x,y) be the coordinates of a generalised Voronoi vertex of Ci, C2, and C 3 . Let r be the local distance between the generalised Voronoi vertex of C\, C2, and C 3 and C\. Let R be the local distance between the generalised Voronoi vertex of Ci, C2, and C 3 and C4. Let X be the set of all possible values of (x,y,r,R) that 123 Figure 5.2.4: The new Delaunay graph after insertion of C 4 : the edges in plain thin lines remain, those in dashed lines disappear, and those in plain thick lines appear Figure 5.2.5: The Delaunay graph conflict locator for conies 124 are solutions to the system of four equations composed of the implicit equations of the three r—generalised offsets to C i , C 2 , and C3 and the R—generalised offset to C 4 in the four unknowns x,y,r,R, such that (x,y) is a true Voronoi vertex of C\, C 2 , and C3. The Delaunay graph conflict locator is computed in two phases by evaluating the sign of the polynomial R — r among the solutions (x, y, r, R)oi the system com-posed of the three r—generalised offsets to C\, C 2 , and C3 and the R—generalised offset to C 4 in the four unknowns x,y,r,R, and computing the x,y,r coordinates of the solutions for which R — r has the desired sign, in order to isolate them. The computation of the x,y,r coordinates is not used as an intermediary computation of the Delaunay graph conflict locator computation. It is used in order to check that the coordinates of a solution of the preceding system of equations is a true Voronoi vertex (this is done in a second phase which is explained below). If one of the points of X is such that R — r < 0, then there is a change in the Delaunay graph. The true Voronoi vertices are the generalised Voronoi vertices whose empty circle does not conflict with any of the first three (defining) conies. This is checked in the second phase by evaluating the sign of R-r on the solutions of each one of the systems composed of the three r—generalised offsets to C\, C2, and C3 and the R—generalised offset to a fourth conic, where the fourth conic is alternatively each one of the first three conies, and isolating the x,y,r coordinates of the solutions for which R — r> 0. The true Voronoi vertices (x, y) are those for which R — r > 0 for any solution (x, y, r, R) of any system composed of the three r—generalised offsets to C i , C 2 , and C3 and the R—generalised offset to C4 in the four unknowns x,y,r,R, where C4 is alternatively C i , C 2 , and C3. We can summarise this paragraph with the following theorem: T h e o r e m 5.2.2. The Delaunay graph conflict locator output list is empty if, and only if, R — r does not take a negative real value on the points of X. 125 5.3 The simplification of the Delaunay graph conflict locator In this section, we will present the simplification of the algebraic equations specify-ing the generalised Voronoi vertices of C\, C2 and C 3 and their local distances r to C i or C 2 or C 3 , and R to C 4 . This simplification allowed us to compute the sparse resultant needed for the algebraic computation of the Delaunay graph conflict loc-ator. Since the offsets to ellipses or hyperbolas have degree 8 while the offsets to parabolas have degree 6 (see Table 4.1), we will call C i a conic whose offset has maximum degree among C i , C 2 , and C 3 : one of the ellipses or hyperbolas among C\, C 2 , and C 3 if there is one, or a parabola if there are no ellipses nor hyperbolas among C i , C 2 , and C 3 . We will suppose without loss of generality that we can make a change of coordinate system such that the x-axis of our new coordinate system is parallel to the main axis of symmetry of CL , and that the origin of the new coordin-ate system is the summit of C i if C i is a parabola or its centre if C i is a conic with centre (ellipse or hyperbola). Thus in this new system of coordinates, the equation 2 2 rt of C i is either ^±p- — 1 = 0 or y — 2pix — 0. The equations of the conic Cj , i ^ 1 can be obtained from the generic equation of a conic in a system centered on the summit/centre of Cj and with x-axis the main axis of symmetry of Ci by applying an affine transformation given by: x = a ; X — BiY + 7, and y = fyX + a j Y + <5j. In the same way we can obtain an implicit equation of the r—generalised offset to Cj from the generic equation of the r—generalised offset to Cj in a system centered on the summit/centre of Cj and with x-axis the main axis of symmetry of Cj by applying the affine transformation given by: x = ctiX — BiY + 7J and y = BiX + aiY + <5j. We have already studied and obtained an implicit equation of the generalised offset to the different types of conies in Section 4.4. Observe that the variable r appears 126 Figure 5.3.1: The Newton polytope of the implicit equation of the r—generalised offset to a parabola in a system centred on its summit and with x-axis the axis of symmetry of the parabola always as powers of r2 since the only term in which it occurs is the term r2 of the polynomials di (x,y,u,v) expressing the distance between the point on the curve Cj and the point on the generalised offset. Therefore, we can consider the following change of variable r2— > r for all the conies. Although, we have to eliminate the four variables x, y, r, R, and thus we need to consider the Newton polytopes in K4, we represented all the Newton polytopes of this section in K3 since the variable R only appears in the implicit equation of the R—generalised offset to the fourth conic. A n implicit equation of the r—generalised offset to a parabola of parameter p in a system centred on its summit and with x-axis its axis of symmetry has been placed in Appendix B.3. The Newton polytope of the implicit equation of the r—generalised offset to a parabola in a system centred on its summit and with x-axis the axis of symmetry of the parabola is shown on Figure 5.3.1. The Newton polytope of the implicit equation of the r—generalised offset to a parabola in an arbitrary system is shown on Figure 5.3.2. However, if we subtract the implicit equation of the r—generalised offset 127 Figure 5.3.2: The Newton polytope of> the implicit equation of the i—generalised . offset to a parabola in an arbitrary system to a parabola in a nice system centred on its summit and with x-axis its axis of symmetry from this equation, the term in r 3 disappears, and the corresponding Newton polytope is shown on Figure 5.3.3. A n implicit equation of the generalised offset to an ellipse or hyperbola of big axis a and small axis 6 in a system centred on the centre of the conic and with axes the axes of symmetry of the conic has been placed in Appendix B.4. The Newton polytope of the implicit equation of the r—generalised offset to an ellipse or hyperbola in a system centred on the centre of the conic and with axes the axes of symmetry of the conic is shown on Figure 5.3.4. The Newton polytope of the implicit equation of the i—generalised offset to an ellipse or an hyperbola in an arbitrary system is shown on Figure 5.3.5. However, in a linear combination of the implicit equations of the r—generalised offset to a parabola in a nice system centred on the centre of the conic and with axes the axes of symmetry of the conic and in an arbitrary system, the term in r 4 disappears, and the corresponding Newton polytope is shown on Figure 5.3.6. 128 Figure 5.3.3: The Newton polytope of the difference of the implicit equations of the generalised offset to a parabola in a nice system and in an arbitrary system Figure 5.3.4: The Newton polytope of the implicit equation of the r-generalised offset to an ellipse or an hyperbola in a system centred on the centre of the conic and with axes the axes of symmetry of the conic 129 Figure 5.3.5: The Newton polytope of the implicit equation of the r-generalised offset to an ellipse or an hyperbola in an arbitrary system Figure 5.3.6: The Newton polytope of a linear combination of the implicit equations of the generalised offset to an ellipse or an hyperbola in a nice system and in an arbitrary system 130 5.4 The algebraic precomputations In this section, we present how the simplifications described in the preceding section allowed us to compute the matrix of the sparse resultant needed for the algebraic computation of the Delaunay graph conflict locator of conies. Let Oj denote the polynomial expressing the implicit equation of the gener-alised offset of the conic Ci in the coordinate system based on C\. In order to evaluate the sign of the polynomial R — r among the solutions (x,y,r, R) of the system composed of the three i—generalised offsets to C\, C2, and C 3 and the R—generalised offset to C 4 in the four unknowns x,y,r,R, we evaluate the sparse resultant matrix (i.e. the Newton matrix) of the polynomials R — r, °i> ° 2 , ° 3 , a n d 0 4 . Indeed, the matrix of the multiplication map by R — r in the quotient algebra K[x, y, r, R]/ < 0 1 , 0 2 , 0 3 , 0 4 > is the transpose of the Schur com-plement of the submatrix M\\ of the sparse resultant matrix of the polynomials R — r, o\, o2, 0 3 , and 0 4 (see Section 3.3.1). The sparse resultant matrix is a square matrix of size up to 7995 (in the case of four ellipses/hyperbolas). The com-putation of this Schur complement involves the computation of the inverse of the matrix M\\, which is a square matrix of up to 7995 — 1024 = 6971 rows (in the case of four ellipses/hyperbolas). Even with a U N I X workstation equipped with 4 Gb of R A M (cosimo or ginevra machines at Medicis [CNR]), that computation on a symbolic matrix halts due to a lack of memory. Thus, the Schur complement computation must be done after having replaced formal coefficients by their numer-ical values (fractions). The different possible configurations and the corresponding mixed volumes and the degree of the sparse resultant are shown in Table 5.1. The sparse resultant computations were done using Singular [GPS01], on both an Apple PowerBook G4 with 1Gb of R A M running under Mac OS X and the cosimo ma-chine at Medicis [CNR] (a Compaq DS20E - Alpha E V 6 with 4Gb of R A M running under OSF/1 4.Of). We have given the Singular [GPS01] code for the cases of four 131 C i c2 Cz c 4 MV-2 MV-3 M V 1 4 MV-o deg (RA) p P P P 108 108 108 108 324 756 p P P E h 144 144 144 108 432 972 E h P P P 108 144 144 144 432 972 E h E h P P 144 144 192 192 576 1248 E h E h E h P 192 192 192 256 768 1600 E h E h P P 144 192 192 144 576 1248 E h E h P E h 192 192 256 192 768 1600 E h E h E h E h 256 256 256 256 1024 2048 Table 5.1: The mixed volumes and sparse resultant degrees of the different config-urations of conies (Eh stands for ellipse or hyperbola, P stands for parabola) parabolas and four ellipses/circles/hyperbolas (see Appendix C). 5 . 5 The numerical computation of the Delaunay graph conflict locator Observe that the sparse resultant matrix of the Delaunay graph conflict locator is very sparse: for the case of four ellipses/circles/hyperbolas, the sparse resultant matrix is a 7995 x 7995 matrix with only 589610 non-zero entries, i.e. a density of 0.92%! The computation of the Schur complement and of the eigenvalues should be done using methods for sparse matrices. Such computations using standard matrix methods can take several days on Singular [GPS01]. The computation of eigenvalues with A R P A C K takes O (N) time where N is the number of lines of the square matrix. The computation of the x,y,r coordinates can be done in theory by simultaneous orthogonalisation of the matrix of the multiplication operator (see Section 3.3.1). The positions (row or column number) of solutions for which R — r has the desired sign are first stored. Then, the coordinates corresponding to R—r having the desired sign can be read on the diagonal of the corresponding orthogonalised multiplication operator matrix. In practice however, the computation of the Schur complement of the matrix M\\ involves the computation of the inverse of the matrix M\\. In the 132 Figure 5.5.1: The test case with four ellipses case of four ellipses/circles/hyperbolas, the size of M\\ is 7995 — 1024 = 6971. We have chosen a common case of four ellipses for the computation of the Delaunay graph conflict locator (see Figure 5.5.1). The computation of the Schur complement takes 1 minute and half, and the computation of the eigenvalues takes also 1 minute and half. The main problem in the computation of the Schur complement is that the matrix M\\ is badly conditioned due to the fact the numerical coefficients of this matrix are very huge. This induces a precision problem in Matlab [LM90]. This is due to the fact Matlab does all the computations considering matrix coefficients as doubles. The alternative is to do the computation of the Schur complement of the Newton matrix using a Computer Algebra System (like Maple [CGGL92], Singular [GPS01], or Maxima [GG82]) that does exact computations on fractional numbers, and then use a method of certification (like [Kra92]) that computes tight bounds on the intervals taken by the eigenvalues. The problem with this computation is that it takes far too much time for being a practical implementation within the context of a randomised incremental algorithm for the construction of the Voronoi diagram of conies. Another alternative is to compute the intervals taken by the entries of the Schur complement, and then certify the computation of the eigenvalues. However, we can't apply the certification method for computing eigenvalues of matrices with exact entries [Kra92]. This is due to the fact the Schur complement is not a continuous 133 function of the matrix entries. We already observed this problem in the simultaneous orthogonalisation of the matrices of the multiplication operators in Section 3.2. Thus, we can conclude that even though we could compute the Delaunay graph conflict locator "exactly" by computing the Schur complement of a submat-rix of the sparse resultant (or Newton) matrix with a Computer Algebra System that does exact computations on fractions, this method would not lead to a prac-tical implementation of the conflict locator that is fast enough for the randomised incremental construction of the Voronoi diagram of conies. The main problem if we compute the Schur complement already mentioned above using floating point computations is that we must certify the results using interval analysis on the Schur complement computation, and then we should do the same thing (using interval analysis) on the computation of the eigenvalues. 134 Chapter 6 The Delaunay graph conflict locator for semi-algebraic sets In this chapter, we will first present a Delaunay graph conflict locator for sites being semi-algebraic sets. This conflict locator checks whether the addition of a given semi-algebraic set would remove or not each facet of the Delaunay graph DG (S) of a set S of semi-algebraic sets whose vertices are a given (N + 1)—tuple of semi-algebraic sets. Thus, the input is an (iV + 2)—tuple of semi-algebraic sets, where the first N + 1 semi-algebraic sets define one or more facets of the Delaunay graph, and the (N + 2)*^ semi-algebraic set is the semi-algebraic set being added. This conflict locator will be based on the A L I A S library [MerOO], whose important features used here have been presented in Section 3.3.3. We will first present the system of equations and inequalities that must be satisfied if the conflict locator outcome is positive (the addition of the semi-algebraic set would remove one or more triangles) in Section 1. We will then present in Section 2 how we check the existence of solutions to the preceding system of algebraic equations and inequalities using A L I A S (with different solving techniques and corresponding parameters). We will then present in Section 3 how the implicit equation of the generalised offset to a conic (see Section 4.4) can be used in order to accelerate the computation of the 135 conflict locator in the case the semi-algebraic sets are conies. This involves both symbolic algebraic precomputations and scientific computations. 6.1 The algebraic equations and inequalities of the De-launay graph conflict locator The definition of a semi-algebraic set has already been presented in Definition 1.0.1 in Chapter 1. Let X\,Xjv+2, be semi-algebraic sets. Let us suppose that each semi-algebraic set Xi is defined as UjLi flfc=i {x e ^N\fi,j,k *i,j,k 0}, where fijtk € R , xiN\ and *ijtk is either < or =, for i = 1,2,3,4, j — l , . . . , S j and k = 1 , r ^ j . Let us assume without loss of generality that each P|fe=i {x e ^-N\fi,j,k *i,j,k 0} for each Xi is defined by at least one non-trivial algebraic equation (i.e. different from the zero polynomial). This is equivalent to each component of each semi-algebraic set being defined as the restriction of a proper closed set in the Zariski topology sense. Thus each component of each semi-algebraic set has a codimension greater or equal to 1. If our starting assumption is not valid in the case we treat, we can make it valid by adding the equations corresponding to fijtk — 0 for each k) such that j is the index of a component that is not defined as in the assumption and i is the index of the semi-algebraic set to which the component belongs. This realises the (real) topological closure of the set of solutions of the fij^ *ij,k 0 such that j is the index of a component that is not defined as in the assumption and i is the index of the semi-algebraic set to which the component belongs. Let us denote Vi as the intersection of all the V (/ij,fc) such that *ijtk is = for each i — 1,2,3,4. Let Mi be the normal space to Vi at the point Xi = (x^, ...,XiN). Each / j j ^ defining Vi induces N — 1 polynomials ^i,j,fc,; with I = 1,...,N — 1 that are the equations defining the normal to V (fij,k) at Xi. A point q = (yi, ...,yyv) belongs to Mi if its coordinates satisfy each one of the equations of the normal spaces to V (fij^k) at Xi 136 such that *ij,fc is =. For agiveng = (yi, ...,yjv), let Aft be the the set of points rrii = (z^, ...,ZiN) € Xi such that g belongs to the normal space to Vi at the point rrii. In the general case, each set Mi is a finite set of points. However, if Vi contains a portion of hypersphere PHS (q,p) centered on q, then Aii contains that portion of hypersphere; To get in all cases a finite set of points rrii of Vi, we use &i = Aii when Aii is finite, and © j P| PHS (q, p) = {wi} for an arbitrary point Wi of PHS (q,p) when Vi contains a portion of hypersphere PHS (q,p) centered on q. We are now able to write the system of algebraic equations and inequalities that define the outcome of the Delaunay graph conflict locator. Let us consider the map 7r : K3N —• KN defined by TT (xi,q, rrii) = Q-The point q is at the distance r from the point X j if, and only if, the distance between q and X j is r. This is expressed algebraically by the equation di (q, X j ) = (yi ~ xh)2 + ... + (yN - xiNf - r2 = 0. The generalised r-offset Oi to Xi is the image by TT of the points of K2N defined by the following system of equations and inequalities: 3j£ [ l , S i ] , V / c € [l,rij], fi,j,k (pi) *ijtk 0 if •,• ~ • i 1 S « = " ,j,k 1 S , di (xi,q) = 0 V/ = 1,.., N-l, nijM (xi, q) = 0 ( fXil (xi) ^ 0 or . . . or fXiN (xt) # 0 The true r-offset to Xi is obtained as the difference of the generalised r-offset Cj to Xi and the union of each one of the images by TT of the semi-algebraic sets defined 137 by the following system of equations and inequalities for each point rrii of 6i: 3j£ [ l , S i ] , V f c € [l,rid], fi,j,k (m-i) *i,j,k 0 i,j>k l h — i < f / (rrn) = 0 < \/l = 1,..,N -l,nitjtkti(mi,q) = 0 [ [ d(mhq) < 0 It is obvious that a true Voronoi vertex of X\, ...,XN+\ is a point of inter-section of the true r—offsets to X-y,.XJV+I respectively. Now, what is left to write is first, that the true Voronoi vertices of XI,...,XN+I are at the distance R from XN+2, or alternatively, that the true Voronoi vertices of X\, ...,XN+I belong to the true R—offset to XN+2, and finally to evaluate the signs of the (real) values of R — r. Consider the (N + 2)—dimensional points whose first N coordinates are the coordinates of a true Voronoi vertex of X \ , X N + \ , and the remaining two are the distancesr between that true Voronoi vertex and X\, ...,XN+I, and R between that true Voronoi vertex and X ^ + 2 - The Delaunay graph conflict locator output list will not be empty if, and only if, R — r < 0 for one of these (N + 2)—dimensional points. 6.2 The formulation of the Delaunay graph conflict loc-ator with ALIAS A l l these equations and the final inequality can be written in a Maple [CGGL92] file (see Maple program in Appendix D). We have tested different solving pro-cedures as well as different A L I A S parameters for computing the Delaunay graph conflict locator for semi-algebraic sets on conies. We first observed better results with the 3B consistency method (ALIAS/3B=1) than without (ALIAS/3B=0) , both in terms of running time and in terms of number of searched boxes. Consequently, we tested only 3B based searching techniques. We can either apply the 3B consist-138 ency method on all the equations and/or inequalities, or we can specify a subset of equations and/or inequalities on which the 3B consistency method will be applied (the ALIAS/subeq3B list variable). Typically, one chooses the algebraic equations and inequalities with lowest degree for the SubEq3B subset. We can also select a maximal interval range (ALIAS/Max3B) , that is the maximal range that the vari-ables intervals should have in order for the 3B method to start being applied. B y setting this maximum interval range to the maximum variable interval range, we can force the 3B method to be applied from the starting variables intervals. We can also select the tolerance interval range for the 3B method, by setting the ALIAS/Del ta3B variable. Some testing of these A L I A S bisection modes was done on the same ex-ample of ellipses as in Chapter 5 (see Figure 5.5.1). The results are summarised in Table 6.1. The executable was run always on the same machine at off-peak hours (in the evening). The significant results is that single bisection is much faster than mixed bisection, which in turn is faster than full bisection. These running times illustrate the impact of the exponential growth of the number of boxes with the number of variables on the running times. The 3B consistency method improves the results obtained by single bisection. There is a trade-off between the additional time required by the additional 3B interval evaluations and the reduction of the variable intervals by the 3B method. The optimum in the example of Table 6.1 is to start the 3B method after the variable intervals have been reduced by 4 by the normal single bisection process. We have tested the full system of equations and inequalities corresponding to the Delaunay graph conflict locator as well as the partial system for identifying the generalised Voronoi vertices that are true or computing the conflict locator assuming we know which generalised Voronoi vertices are true Voronoi vertices. Note that we need to execute the program N + 2 times instead of one if we use the partial system (N + l times for identifying the true Voronoi vertices and once for computing the Delaunay graph conflict locator assuming the true Voronoi vertices have been 139 Bisection process Running time for optimised GradientSolve Full bisection 2 h 34 min 42 s Mixed bisection (half the variables) 26 min 26 s Single bisection + 3B 2 min 26 s + 3B with half Max3B 2 min 8 s ' . + 3B with one quarter Max3B 2 min 2 s + 3B with one sixth Max3B 4 min 49 s + 3B with one eight Max3B 4 min 49 s Table 6.1: Some running time results with different A L I A S parameters on the system with the generalised offsets (see Table 6.4) Running time without subEq3B with subEq3B optimised General Solve 6 min 38 s 6 min 26 s non-optimised General Solve optimised Gradient Solve 2 h 56 min 10 s 2 h 55 min 4 s non-optimised Gradient Solve optimised Hessian Solve 20 h 17 min 42 s 20 h 19 min 33 s non-optimised Hessian Solve Table 6.2: Some running time results for ellipses with the equations of the original curves identified). The results are summarised in Tables 6.2 and 6.3. While there is a time disadvantage in running the solver on the full system with respect to running the solver (four times) on the partial systems for the Gradi-ent and Hessian based solvers, running the General solver on the full system is less time consuming. What is important to observe at this point is that if the full system is not used, in addition to running four times the solver, we need to check that the coordinates and local distance to the defining semi-algebraic sets of each generalised Voronoi vertex corresponding to the solutions correspond to a true Voronoi vertex. The full system has all the constraints (inequalities in it). This is likely why the full system can be solved faster than the four partial systems. Moreover, the General solver is the fastest one on the full and partial systems based on the equations of the original curves. Finally, specifying a subset of equations and/or inequalities 140 Running time without subEq3B with subEq3B optimised General Solve 2 min 51 s 2 min 38 s non-optimised General Solve optimised Gradient Solve 26 min 52 s 27 min 8 s non-optimised Gradient Solve optimised Hessian Solve 1 h 43 min 42 s 1 h 45 min 38 s non-optimised Hessian Solve 1 h 49 min 9 s Table 6.3: Some running time results for ellipses with the equations of the original curves for the partial system (ALIAS/subeq3B) on which the 3B bisection process will be applied improves the running time of the general solver while increases those of the gradient and hessian based solvers. The general solver seems to be more efficient on the full system based on the equations specifying the semi-algebraic sets because the computations of the gradient and of the hessian are more time consuming when the number of variables increases. A way to improve the running time is to improve the computations of the intervals. What we have done is to use the parser to convert the Maple [CGGL92] expressions into C++ code that uses the A L I A S C + + library in order to do interval computations and identify intervals with possible solutions. The interval computations resulting from this automatised process may not be optimal. This is true because the computation of the interval taken by a function depends on the way this function is written. We have tried to optimise the way functions are written for the interval computations on the general solver code. However, this did not lead to any significant improvement on the running time of the optimised code. It looks like there is a trade-off between the running time improvement owed to the change of expression of the function and the optimisation done by the C + + compiler (g++ version 2.95.2). 141 Running time without subEq3B with subEq3B optimised General Solve 12 min 37 s 12 min 56 s non-optimised General Solve optimised Gradient Solve 2 min 26 s 4 min 57 s non-optimised Gradient Solve optimised Hessian Solve 3 min 41 s 3 min 42 s non-optimised Hessian Solve Table 6.4: Some running time results for ellipses with the equations of the generalised offsets Ratios (without/with generalised offset) without subEq3B with subEq3B optimised General Solve 0.526 0.49.7 non-optimised General Solve optimised Gradient Solve 72.40 35.37 non-optimised Gradient Solve optimised Hessian Solve 330.6 329.6 non-optimised Hessian Solve Table 6.5: The ratios of the running time results without/with the equations of the generalised offsets 6.3 The hybrid symbolic/scientific computation of the Delaunay graph conflict locator for conies In the previous section, we have presented the computation of the Delaunay graph conflict locator for semi-algebraic sets. In this section, we will present how we can use the implicit equation of the generalised offset to a conic in the A L I A S compu-tations, and show some results that illustrate the considerable reduction of running time realised by using the equations of the generalised offsets instead of the equa-tions of the original curves. Table 6.4 shows some results on the ellipses that were used for the tests of the previous section. These results show a decrease of running time except for the general solver with ratios varying between 35 and 331 (see Table 6.5). 142 Figure 6.3.1: The test case with four hyperbolas In order to see if the same result could be achieved with hyperbolas, we also tested the Delaunay graph conflict locator on hyperbolas. The test case we used is illustrated on Figure 6.3.1. The running times obtained using the original implicit equations of the hyperbolas are summarised in Table 6.6, while the running times obtained using the implicit equations of the generalised offsets to the hyperbolas are summarised in Table 6.7. The general solver and 3B consistency method with spe-cification of a subset of equations to be evaluated in the 3B method on the system with the implicit equations of the generalised offsets gives the best running time with 6 minutes and 20 seconds. Another possibility of expression of the equations of conies is the replacement of the implicit equations of conies by their parametric equations in polar coordin-ates1. This did not lead to any running time improvement with any of the solving :The general polar equation of a non degenerate conic with respect to one of its foci is: r = 1_epcosg, where e is the eccentricity of the conic and p is that eccentricity times the distance from the foci to the directrix. 143 Running time without subEq3B with subEq3B optimised General Solve 17 min 23 s . 15 min 40 s non-optimised General Solve optimised Gradient Solve l h 10 min 33 s l h 07 min 28 s non-optimised Gradient Solve optimised Hessian Solve l h 45 min 11 s non-optimised Hessian Solve l h 36 min 0 s Table 6.6: Some running time results for hyperbolas with the original equations of the hyperbolas Running time without subEq3B with subEq3B optimised General Solve 10 min 8 s 6 min 20 s non-optimised General Solve optimised Gradient Solve 1 h 2 min 2 s 35 min 28 s non-optimised Gradient Solve optimised Hessian Solve 2 h 12 min 17 s non-optimised Hessian Solve 1 h 49 min 9 s Table 6.7: Some running time results for hyperbolas with the equations of the generalised offsets 144 techniques used earlier. Chapter 7 Conclusions The theoretical purpose of this thesis as defined in Chapter 1 is the elucidation of the basic algebraic and geometric properties of the generalised offset to a curve that are central to the Delaunay graph conflict locator. One of the practical objectives (defined in Chapter 1) is the computation of the Delaunay graph conflict locator for algebraic varieties in the case of conies. The other practical purpose of this thesis is the computation of the Delaunay graph conflict locator for low degree (semi-) algebraic sets embedded in the Euclidean plane. We will briefly review the achievements obtained in this thesis. Wi th respect to the theoretical objectives, we have obtained what we fixed as objectives. A general formula for the degree of the generalised offset to an algebraic curve has been presented in Section 4.2. The degree of the generalised offset to a conic and the degree of the Delaunay graph conflict locator for conies have been computed (see Sections 4.3 and 5.4 respectively). We have also reduced the problem of the computation of the Delaunay graph conflict locator from a semi-algebraic problem to a linear algebra problem (computing the eigenvalues of a matrix, as described in Chapter 5). Wi th respect to the practical objectives of this thesis, we have obtained what we fixed as objectives with algebraic precomputations limited to the im-146 plicit equation of the generalised offset to a conic and the numerical computa-tions done with the interval analysis based solver A L I A S . The implicit equation of the generalised offset of a conic defined by a polynomial with formal coefficients (ax 2 + by2 + cxy + dx + ey + / = 0) has been computed symbolically (see Section 4.4). The exact Delaunay graph of the Voronoi diagram for circles has been computed symbolically (see Section 3.1.3). This allows one to construct the Voronoi diagram of circles exactly. The matrix of the sparse resultant of the polynomials specifying the Delaunay graph conflict locator for conies was computed symbolically (see Sec-tion 5.4). The Schur complement of a submatrix of this matrix is the multiplication operator matrix whose eigenvalues allow one to answer the Delaunay graph conflict locator. Its computation has not been done symbolically, but numerically. The main problem with this computation is that it takes far too much time to have a practical interest in the randomised incremental construction of the Voro-noi diagram for conies. The multiplication operator matrix.and of its eigenvalues were computed using Matlab [LM90]. However, this computation cannot be done with enough precision due to the bad conditioning of the matrix of the sparse res-ultant. The certified Delaunay graph of algebraic curves and of semi-algebraic sets defined by numeric polynomials has been computed using the interval analysis and consistency based solver " A L I A S " (see Chapter 6). The certified computation of the Delaunay graph of conies using interval analysis gradient and hessian based solvers can benefit from the use of the implicit equation of the generalised offset to a conic (from 35 to 331 times faster, 2 min 26 s for ellipses, 6 min 20 s for hyperbolas). This result confirms the idea that knowing the structure of the set of solutions may help finding the solutions. Even though there is no known theoretical lower nor upper bound for the interval analysis based solvers, in practice those solvers can compute the Delaunay graph conflict locator much faster than the computation of the eigenvalues of the Schur complement of a submatrix of the sparse resultant matrix. The first time 147 we encountered a result similar to this one was when we found in Section 3.2 that the same computation of the Delaunay graph conflict locator was taking 11 minutes with A L I A S and 12 hours using G B / R S [Fau, Rou] (by computing the Grobner basis and the eigenvalues of the multiplication operator matrices). Moreover, the computation of the Delaunay graph conflict locator by interval analysis based solvers is more direct than the computation of the conflict locator through eigenvalues of the Schur complement of a submatrix of the sparse resultant matrix. Indeed, we can test whether the generalised Voronoi vertices are true Voro-noi vertices and the Delaunay graph conflict locator simultaneously using A L I A S . We do not have to evaluate first the conflict locator assuming the generalised Voronoi vertices are true Voronoi vertices and then test which generalised Voronoi vertices invalid with respect to the fourth conic were true Voronoi vertices of the first three conies. Finally the computation of the Delaunay graph conflict locator by interval analysis based solvers can be used for semi-algebraic sets of arbitrary degree, which is not possible for algebraic techniques because we have reached with conies the limit on the memory that can be used in current systems (4 Gb of R A M and 6 Gb of virtual memory on Medicis machines [CNR]). This is due to the exponential complexity of the size of the matrices involved in the computations. Moreover, the computation of the Delaunay graph conflict locator by interval analysis based solvers can be easily generalised to general regular curves, which is absolutely impossible for algebraic techniques. This thesis has presented what we believe is the first computation of the Delaunay graph conflict locator for semi-algebraic sets (and in particular conies), and its application to a semi-dynamic algorithm for the construction of the Voronoi diagram of semi-algebraic sets. The computations involved in the certified Delaunay graph conflict locator may be required only in almost degenerate cases where the semi-algebraic set being 148 added touches or almost touches one or more Delaunay graph empty circles. This situation is not the most probable one, at least statistically. The algorithms presen-ted in this thesis could be combined with simpler algorithms when the degenerate cases can be filtered. 7.1 L i m i t a t i o n s of this research The implicit equation of the generalised offset to a conic is quite long and complex. We could have tried to use invariants in order to try to simplify it after we obtained the implicit equation. A simplification of this implicit equation of the generalised offset to a conic would imply simpler polynomials (i.e. with fewer monomials and maybe a lower degree). Simpler implicit equations for the generalised offsets would imply that the sparse resultant matrix of the Delaunay graph conflict locator would be smaller, provided the new variables (the invariants used) allow one to identify unambiguously the generalised Voronoi vertices and their local distances to the first three conies and to the fourth conic. The problem here is to find such a set of new variables (the invariants used). The invariants of the special orthogonal group for the generalised Voronoi vertex look like the only invariants that allow one to identify unambiguously the generalised Voronoi vertices and their local distances to the first three conies and to the fourth conic. Wi th respect to the Computer Algebra Systems, I tested a large number of them: all the public domain ones (CoCoA [CNROO], Macaulay 2 [GS], Maxima [GG82], and Singular [GPS01]), and the most commonly commercial one used for Grobner bases computations (Magma), as well as general Computer Algebra Sys-tems such as Maple [CGGL92] or Mathematica [W6199]. Wi th regard to numerical computations, I have tried only interval analysis based methods. Other methods such as homotopy continuation [Li03] allow one to compute the solutions of systems of polynomial equations. Wi th respect to the interval analysis based solvers, I tested only A L I A S for two reasons. The first one is 149 the availability of the designers of A L I A S at INRIA Sophia-Antipolis, where I made a visit of one year and half. The second one is that A L I A S is an attempt to bring together in a single library tools from interval analysis and consistency methods. There are a lot of interval analysis solvers available either on the public domain or commercially. Other alternatives could have been tested. The libraries used at the lower levels in order to perform interval evaluations for functions are however more limited in number. A L I A S uses the P R O F I L / B I A S [Knii94] library for performing interval evaluations. 7.2 Future research The future research should try to go beyond the limitations of this research as out-lined in the previous section. A n interesting direct generalisation worth to explore is the case of regular curves rather than semi-algebraic sets. However, the prac-tical applications of such a generalisation does not seem to be more important than the case of semi-algebraic sets, because regular curves can be approximated to any geometrical tolerance by semi-algebraic sets. What is critical is that the Voronoi diagram and the Delaunay graph are very sensitive to the continuity of the first order and second order derivatives at contact points. This is the reason for which approximation of curves by line segments does not guarantee the exactness of the Delaunay graph. It would be useful to explore the differences in running times of the Delaunay graph conflict locator computation using regular curves or such an approximation of regular curves by semi-algebraic sets (and especially conies). If the implicit equation of the generalised offset to a conic could be simplified by using new variables such as invariants, it would be worth trying to compute the Grobner basis of the ideal generated by the three r—generalised offsets to three conies C i , C2 and C3 and the R—generalised offset to a fourth conic C4. Knowing the Grobner basis of that ideal, two alternative approaches may be used in order to evaluate the conflict locator: either the linear algebra approach in the quotient 150 algebra (as developed in Section 3.1.3 or in Chapter 5), or the triangular set approach (as implemented in the triang library of Singular [GPS01]). 151 Bibliography [ABMY02] Frahcgois Anton, Jean-Daniel Boissonnat, Darka Mioc, and Mariette Yvinec. A n exact predicate for the optimal construction of the A d - 1 ditively Weighted Voronoi diagram. In Proceedings of the European Workshop on Computational Geometry 2002, Warsaw, Poland, 2002. [AK00] Franz Aurenhammer and Rolf Klein. Voronoi diagrams. In Handbook of computational geometry, pages 201-290. North-Holland, Amsterdam, 2000. [AKM02] Francois Anton, David Kirkpatrick, and Darka Mioc. A n exact algebraic predicate for the maintenance of the topology of the additively weighted voronoi diagram. In The Fourteenth Canadian Conference on Compu-tational Geometry, pages 72-76, Lethbridge, Alberta, Canada, 2002. [AMG98] Francois Anton, Darka Mioc, and Christopher M . Gold. Dynamic A d -ditively Weighted Voronoi diagram made easy. In Proceedings of the 10th Canadian Conference on Computational Geometry, August 1998, Montreal, Canada, pages 92-93, 1998. [AS95] Helmut Al t and Otfried Schwarzkopf. The Voronoi diagram of curved objects. In Proc. 11th ACM Symp. Computational Geometry, pages 89-97. A C M Press, 5-7 June 1995. 152 [ASS99] Enrique Arrondo, Juana Sendra, and J . Rafael Sendra. Genus formula for generalized offset curves. J. Pure Appl. Algebra, 136(3):199-209, 1999. [Aur87] Franz Aurenhammer. Power diagrams: properties, algorithms and ap-plications. SIAM J. Comput., 16(l):78-96, 1987. [Baj94] Chandrajit L . Bajaj. Some applications of constructive real algebraic geometry. In Algebraic geometry and its applications (West Lafayette, IN, 1990), pages 393-405. Springer, New York, 1994. [BB96] John A . Beachy and William D. Blair. Abstract Algebra. Waveland Press Inc., 1996. [BCK88] Bruno Buchberger, George E . Collins, and Bernhard Kutzler. Algebraic methods for geometric reasoning. In Annual review of computer science, Vol. 3, pages 85-119. Annual Reviews, Palo Alto, C A , 1988. [BCR98] Jacek Bochnak, Michel Coste, and Marie-Frangoise Roy. Real algeb-raic geometry. Springer-Verlag, Berlin, 1998. Translated from the 1987 French original, Revised by the authors. [BCSS98] Lenore Blum, Felipe Cucker, Michael Shub, and Steve Smale. Com-plexity and real computation. Springer-Verlag, New York, 1998. Wi th a foreword by Richard M . Karp. [BR90] Riccardo Benedetti and Jean-Jacques Risler. Real algebraic and semi-algebraic sets. Hermann, Paris, 1990. [Buc70] Bruno Buchberger. E in algorithmisches Kriterium fur die Losbarkeit eines algebraischen Gleichungssystems. Aequationes Math., 4:374-383, 1970. 153 [Buc79] Bruno Buchberger. A criterion for detecting unnecessary reductions in the construction of Grobner-bases. In Symbolic and algebraic computa-tion (EUROSAM '79, Internal Sympos., Marseille, 1979), volume 72 of Lecture Notes in Comput. Sci., pages 3-21. Springer, Berlin, 1979. [Buc88] . Bruno Buchberger. Applications of Grobner bases in nonlinear compu-tational geometry. In Mathematical aspects of scientific software (Min-neapolis, Minn., 1986/87), volume 14 of IMA Vol. Math. Appi, pages 59-87. Springer, New York, 1988. [Buc92] Bruno Buchberger. Grobner bases: an introduction. In Automata, lan-guages and programming (Vienna, 1992), volume 623 of Lecture Notes in Comput. Sci., pages 378-379. Springer, Berlin, 1992. [Buc98] Bruno Buchberger. Introduction to Grobner bases. In Grobner bases and applications (Linz, 1998), volume 251 of London Math. Soc. Lecture Note Ser., pages 3-31. Cambridge Univ. Press, Cambridge, 1998. [CCM97] Hyeong In Choi, Sung Woo Choi, and Hwan Pyo Moon. Mathematical theory of medial axis transform. Pacific J. Math., 181(l):57-88, 1997. [CDOO] Tzu-Yi Chen and James W . Demmel. Balancing sparse matrices for computing eigenvalues. In Proceedings of the International Workshop on Accurate Solution of Eigenvalue Problems (University Park, PA, 1998), volume 309, pages 261-287, 2000. [CEOO] John F. Canny and Ioannis Z. Emiris. A subdivision-based algorithm for the sparse resultant. J. ACM, 47(3):417-451, 2000. [CGGL92] Bruce W . Char, Keith O. Geddes, Gaston H . Gonnet, and Benton L . Leoung. Maple V. Springer, Berlin, 1992. 154 [CL097] David Cox, John Little, and Donal O'Shea. Ideals, varieties, and al-gorithms. Springer-Verlag, New York, second edition, 1997. A n intro-duction to computational algebraic geometry and commutative algebra. [CL098] David Cox, John Little, and Donal O'Shea. Using algebraic geometry. Springer-Verlag, New York, 1998. [CNR] CNRS/Ecole Poly technique, http://www.medicis.polytechnique.fr. Unite Mixte de Services Medicis. [CNROO] Antonio Capani, Gianfranco Niesi, and Lorenzo Robbiano. CoCoA, a system for doing Computations in Commutative Algebra. Available via anonymous ftp from cocoa .d ima.un ige . i t , 4.0 edition, 2000. [CPX02] Zhenming Chen, Evanthia Papadopoulou, and Jinhui X u . Robust al-gorithm for k-gon voronoi diagram construction. In Abstracts for the Fourteenth Canadian Conference on Computational Geometry CCCG '02, pages 77-81, Lethbridge, Alberta, Canada, August 2002. Univer-sity of Lethbridge. [Des03] Alexis Deschamps. Handbook of Aluminum, chapter Analytical Tech-niques .for Aluminium Alloys, vol. 2, Alloy Production and Materials manufacturing, pages 155-192. Marcel Dekker, Inc., New York, USA, 2003. [DMT92] Olivier Devillers, Stefan Meiser, and Monique Teillaud. The space of spheres, a geometric tool to unify duality results on Voronoi diagrams. Rapport de recherche 1620, INRIA, 1992. [EC95] Ioannis Z. Emiris and John F . Canny. Efficient incremental algorithms for the sparse resultant and the mixed volume. J. Symbolic Comput., 20(2):117-149, 1995. 155 [Emi96] Ioannis Z. Emiris. On the complexity of sparse elimination. J. Com-plexity, 12(2):134-166, 1996. [Emi97] Ioannis Z. Emiris. A general solver based on sparse resultants: Numer-ical issues and kinematic applications. Technical Report 3110, SAFIR, 1997. [Fau] Jean-Charles Faugere. Gb. Available from http://www-calfo r . I i p 6.fr/ j cf/index.html. [FJ94a] R. T. Farouki and J . K . Johnstone. Computing point/curve and curve/curve bisectors. In Design and application of curves and sur-faces (Edinburgh, 1992), pages 327-354. Oxford Univ. Press, New York, 1994. [FJ94b] Rida T. Farouki and John K . Johnstone. The bisector of a point and a plane parametric curve. Comput. Aided Geom. Design, 11(2):117—151, 1994. [FN90a] Rida T. Farouki and C. Andrew Neff. Algebraic properties of plane offset curves. Comput. Aided Geom. Design, 7(1-4): 101-127, 1990. Curves and surfaces in C A G D '89 (Oberwolfach, 1989). [FN90b] Rida T. Farouki and C. Andrew Neff. Analytic properties of plane offset curves. Comput. Aided Geom. Design, 7(1-4) :83-99, 1990. Curves and surfaces in C A G D '89 (Oberwolfach, 1989). [FR98a] Rida T. Farouki and Rajesh Ramamurthy. Degenerate point/curve and curve/curve bisectors arising in medial axis computations for planar domains with curved boundaries. Comput. Aided Geom. Design, 15(6):615-635, 1998. 156 \ [FR98b] Rida T. Farouki and Rajesh Ramamurthy. Specified-precision compu-tation of curve/curve bisectors. Internat. J. Comput. Geom. Appl., 8(5-6):599-617, 1998. [GG82] V . Ellen Golden and Mathlab Group. Introductory MACSYMA docu-mentation: a collection of papers: (a) An introduction to. ITS for the Macsyma user, (b) ITS easy once ITS explained, and (c) Macsyma primer. Massachusetts Institute of Technology, Laboratory for Com-puter Science, Cambridge, M A , USA, revised edition, 1982. [Giu84] Marc Giusti. Some effectivity problems in polynomial ideal theory. In EUROS AM 84 (Cambridge, 1984), volume 174 of Lecture Notes in Com-put. Sci., pages 159-171. Springer, Berlin, 1984. [GP02] Gert-Martin Greuel and Gerhard Pfister. A Singular introduction to commutative algebra. Springer-Verlag, Berlin, 2002. Wi th contributions by Olaf Bachmann, Christoph Lossen and Hans Schonemann, Wi th 1 C D - R O M (Windows, Macintosh, and UNIX) . [GPS01] Gert-Martin Greuel, Gerhard Pfister, and Hans Schonemann. SIN-GULAR 2.0. A Computer Algebra System for Polynomial Computa-tions, Centre for Computer Algebra, University of Kaiserslautern, 2001. h t t p : / / w w w . s i n g u l a r . u n i - k l . d e . [Gr639] Wolfgang Grobner. Uber die algebraischen eigenschaften der integrale von linearen differentialgleichungen mit konstanten koeffizienten. Mon-atsh. der Math., 1939. [GS] Daniel R. Grayson and Michael E . Stillman. Macaulay 2, a software system for research in algebraic geometry. Available at http://www.math.uiuc.edu/Macaulay2/. 157 [GS85] Leonidas Guibas and Jorge Stolfi. Primitives for the manipulation of general subdivisions and the computation of voronoi. ACM Transactions on Graphics (TOG), 4(2):74-123, 1985. [Han92] Eldon Hansen. Global optimization using interval analysis, volume 165 of Monographs and Textbooks in Pure and Applied Mathematics. Marcel Dekker Inc., New York, 1992. [Hea02] Michael T. Heath. Scientific Computing : An Introductory Survey. McGraw-Hill Higher Education, second edition edition, 2002. [Hof90] Christoph M . Hoffmann. Algebraic and numerical techniques for offsets and blends. In Computation of curves and surfaces (Puerto de la Cruz, 1989), pages 499-528. Kluwer Acad. P u b l , Dordrecht, 1990. [Huy86] . Dung T. Huynh. A superexponential lower bound for Grobner bases and Church-Rosser commutative Time systems. Inform, and Control, 68(l-3):196-206, 1986. [HV91] Christoph M . Hoffmann and Pamela J . Vermeer. Eliminating extraneous solutions in curve and surface operations. Intemat. J. Comput. Geom. Appl, l(l):47-66, 1991. [Kan57] Leonid Vitalyevich Kantorovitch. On some further applications of the Newton approximation method. Vestnik Leningrad. Univ. Ser. Mat. Meh. Astr., 12(7).-68-103, 1957. [KE02] Menelaos I. Karayelas and Ioannis Z. Emiris. Predicates for the Planar Additively Weighted Voronoi Diagram. E C G Technical Report E C G -TR-122201-01, INRIA, 2002. [KKS00] Deok-Soo K i m , Donguk K i m , and Kokichi Sugihara. Voronoi diagram of a circle set constructed from Voronoi diagram of a point set. In Al-158 gorithms and computation (Taipei, 2000), volume 1969 of Lecture Notes in Comput. Sci., pages 432-443. Springer, Berlin, 2000. [KKSOla] Deok-Soo K i m , Donguk K i m , and Kokichi Sugihara. Voronoi diagram of a circle set from Voronoi diagram of a point set. I. Topology. Comput. Aided Geom. Design, 18(6):541-562, 2001. [KKSOlb] Deok-Soo K i m , Donguk K i m , and Kokichi Sugihara. Voronoi diagram of a circle set from Voronoi diagram of a point set. II. Geometry. Comput. Aided Geom. Design, 18(6):563-585, 2001. [Kle89] Rolf Klein. Concrete and abstract Voronoi diagrams. Springer-Verlag, Berlin, 1989. [Knii94] Olaf Kniippel. P R O F I L / B I A S — a fast interval library. Computing, 53(3-4):277-287, 1994. International Symposium on Scientific Computing, Computer Arithmetic and Validated Numerics (Vienna, 1993). [Kol37] Andrei Nikolaevich Kolmogorov. A statistical theory for the recrystal-lization of metals. Akad. nauk SSSR, Izv., Ser. Matem., l(3):355-359, 1937. [Kra92] Walter Kramer. Verified solution of eigenvalue problems with sparse matrices. In Computational and applied mathematics, I (Dublin, 1991), pages 277-287. North-Holland, Amsterdam, 1992. [Lan02] Serge Lang. Algebra, volume 211 of Graduate Texts in Mathematics. Springer-Verlag, New York, third edition, 2002. [Li03] Tien-Yien L i . Numerical solution of polynomial systems by homotopy continuation methods. In Handbook of numerical analysis, Vol. XI, Handb. Numer. Anal., X I , pages 209-304. North-Holland, Amsterdam, 2003. 159 [LM90] Jack N . Little and Cleve B. Moler. Matlab — User's Guide. MathWorks, Inc., Cochinate Place, 24 Prime Park Way, Natick, M A 01760, January 1990. [LSY98] Rich B . Lehoucq, Danny C. Sorensen, and Chao Yang. ARPACK users' guide. Software, Environments, and Tools. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, P A , 1998. Solution of large-scale eigenvalue problems with implicitly restarted Arnoldi methods. [MB97] Frangois Mercier and Olivier Baujard. Voronoi diagrams to model forest dynamics in french guiana. In Proceedings of Geo Computation '97 & Sire '97, pages 161-171, University of Otago, New Zeland, August 1997. [Mer96] Jean-Pierre Merlet. Some algebraic geometry problems arising in the field of mechanism theory. In Algorithms in algebraic geometry and applications (Santander, 1994), pages 271-283. Birkhauser, Basel, 1996. [MerOO] Jean-Pierre Merlet. Alias: an interval analysis based library for solv-ing and analyzing system of equations. In SEA, Toulouse, France, 14-16 June 2000. [MerOl] Jean-Pierre Merlet. A parser for the interval evaluation of analytical functions and its application to engineering problems. J. Symbolic Corn-put, 31(4):475-486, 2001. [Mio02] Darka Mioc. The Voronoi spatio-temporal data structure. PhD thesis, Universite Laval, 2002. [MM84] H . Michael Moller and Ferdinando Mora. Upper and lower bounds for the degree of Groebner bases. In EUROS AM 84 (Cambridge, 1984), volume 174 of Lecture Notes in Comput. Sci., pages 172-183. Springer, Berlin, 1984. 160 [Moo77] Ramon E. Moore. A test for existence of solutions to nonlinear systems. SIAM J. Numer. Anal, 14(4):611-615, 1977. [Mou96] Bernard Mourrain. Enumeration problems in geometry, robotics and vision. In Algorithms in algebraic geometry and applications (Santander, 1994), pages 285-306. Birkhauser, Basel, 1996. [NikOO] Jorgen L. Nikolajsen. An improved Laguerre eigensolver for unsymmet-ric matrices. SIAM J. Sci. Comput, 22(3).-822-834 (electronic), 2000. [OBS92] Atsuyuki Okabe, Barry Boots, and Kokichi Sugihara. Spatial tessella-tions: concepts and applications of Voronoi diagrams. John Wiley & Sons Ltd., Chichester, 1992. With a foreword by D. G. Kendall. [OSY86] Colm 6'Dunlaing, Micha Sharir, and Chee-K. Yap. Generalized Voronoi diagrams for moving a ladder. I. Topological analysis. Comm. Pure Appl. Math., 39(4):423-483, 1986. [OSY87] Colm 6'Dunlaing, Micha Sharir, and Chee Yap. Generalized Voronoi diagrams for a ladder. II. Efficient construction of the diagram. Algorith-mica, 2(l):27-59, 1987. [OY85] Colm 6'Dunlaing and Chee-K. Yap. A "retraction" method for planning the motion of a disc. J. Algorithms, 6(1):104—111, 1985. [Pin02] Giorgio Pini. Leftmost eigenvalue of real and complex sparse matrices on parallel computer using approximate inverse preconditioning. Parallel Algorithms Appl., 17(l):41-58, 2002. [RF99a] Rajesh Ramamurthy and Rida T. Farouki. Voronoi diagram and medial axis algorithm for planar domains with curved boundaries. I. Theoretical foundations. J. Comput. Appl. Math., 102(1):119-141, 1999. Special issue: computational methods in computer graphics. 161 [RF99b] Rajesh Ramamurthy and Rida T. Farouki. Voronoi diagram and medial axis algorithm for planar domains with curved boundaries. II. Detailed algorithm description. J. Comput. Appl. Math., 102(2):253-277, 1999. [Rou] Fabrice Rouillier. RS RealSolving. Available from h t t p : / / s p a c e s . I i p 6 . f r / r o u i l l i e / S o f t s / R S / i n d e x . p h p . [Sha94] Igor R. Shafarevich. Basic algebraic geometry. 1. Springer-Verlag, Ber-lin, second edition, 1994. Varieties in projective space, Translated from the 1988 Russian edition and with notes by Miles Reid. [VC90] Christine Voiron-Canicio. Analyse spatiale et analyse d'images par la morphologie mathematique. R E C L U S , 1990. [vdV99] Henk van der Vorst. Subspace iteration methods with preconditioning for eigenvalues of very large matrices. In Recent advances in numerical methods and applications, II (Sofia, 1998), pages 124-134. World Sci. Publishing, River Edge, N J , 1999. [Vor07] Georgii Feodosevich Voronoi. Nouvelles applications des parametres continus a la theorie des formes quadratiques. premier memoire. sur quelques proprietes des formes quadratiques positives parfaites. Journal fur die reine und angewandte Mathematik, 133:97-178, 1907. [Vor08] Georgii Feodosevich Voronoi. Nouvelles applications des parametres continus a la theorie des formes quadratiques. deuxieme memoire. recherches sur les paralleloedres primitifs. premiere partie. partition uni-forme de l'espace analytique a n dimensions a l'aide des translations d'un meme polyedre convexe. Journal fur die reine und angewandte Math-ematik, 134:198-287, 1908. [VorlO] Georgii Feodosevich Voronoi. Nouvelles applications des parametres continus a la theorie des formes quadratiques. deuxieme memoire. 162 recherches sur les paralleloedres primitifs. seconde partie. domaines de formes quadratiques correspondant aux differents types de paralleloedres primitifs. Journal fur die reine und angewandte Mathematik, 136:67-181, 1910. [Wol99] Stephen Wolfram. The Mathematica book. Cambridge University Press and Wolfram Research, Inc., New York, NY, USA and 100 Trade Center Drive, Champaign, IL 61820-7237, USA, fourth edition, 1999. [Yap87] Chee-K. Yap. An O(nlogn) algorithm for the Voronoi diagram of a set of simple curve segments. Discrete Comput. Geom., 2(4):365-393, 1987. 163 Appendix A The Macaulay 2 program for the exact Delaunay graph conflict locator for the Additively Weighted Voronoi diagram gbTrace 4 dim FractionField := F -> 0 P = frac(QQ[a,b,c,d,e,f,g,h,i,j,k,l]) R = P[x,y,t] c e r c l e l = (x-a)~2+(y-b)~2-(c+t)~2 cercle2 = (x-d)~2+(y-e)~2-(f+t)~2 cercle3 = (x-g)~2+(y-h)~2-(i+t)~2 emptycircle = ideal(cerclel,cercle2,cercle3) ecgb = gb emptycircle print ecgb eckb = basis cokernel gens ecgb print eckb 164 k l = sort(flatten(entries(eckb))) kmind = splice {0..#kl - 1} scan(kl,entry->print ring entry); hashlist = pack(2,mingle(kl,kmind)); feetmon = applyKeys(hashTable hashlist, key->toString(key)); compmat = f -> (htl=apply(kl,be-> hashTable(pack (2 ,mingle(apply(flatten(entries((coefficients((f*be) %ecgb))#0)), item -> feetmon#(toString(item))),flatten(entries((coefficients ((f*be)y„ecgb))#l)))))); matrix(table(#kl,#kl,(i,j)->if (htl#i)#?j then (htl#i)#j else 0 ) ) ) ; matp2 = compmat((x-j)~2+(y-k)~2-(t+l)"2); mOO = matp2_(0,0) mOl = matp2_(0,l) mlO = matp2_(l ,0) mil = matp2_(l,l) cmOO = coefficients mOO cmOOO = cm00#0 cmOOl = cm00#l cmOl = coefficients mOl cmOlO = cm01#0 cmOll = cm01#l cmlO = coefficients mlO cmlOO = cmlO#0 cmlOl = cmlO#l cmll = coefficients mil cmllO = cmll#0 165 c m l l l = c m l l # l 166 Appendix B A n implicit equation of the generalised offset to a conic .1 The Maple program for calling Resin for obtaining the matrix of the sparse resultant > res tar t ;conique:=a*x"2+b*x*y+c*y~2+d*x+e*y+f; conique := ax2 + bxy + cy2 + dx + ey + f . > n o r m a l e : = e x p a n d ( - d i f f ( c o n i q u e , y ) * ( u - x ) + d i f f ( c o n i q u e , x ) * ( v - y ) ) ; normale :=— bxu + bx2 — 2cyu + 2cyx — eu + ex + 2axv — 2axy + byv — by2 + dv — dy > d i s t a n c e : = ( u - x )~2+ ( v - y )~2 - r~2; distance := (u — x)2 + (v — y)2 — r2 167 > collect(normale,[x,y]); > collect(expand(X*conique+Y*distance),[x,y]); > coeff s(collect(expand(X*conique+Y*distance) , [x,y]) ,x, 'al'') ;al; > alpha:=coeffs(collect(expand(X*conique+Y*distance),[x,y]),x)[3] > jcoeffs(collect(expand(X*conique+Y*distance),[y,x]),y,'a2');a2; > beta:=coeffs(collect(expand(X*conique+Y*distance),[y,x]),y)[3]; > solve({alpha=b,beta=-b},{X,Y});simplify(expand(subs(solve( > {alpha=b,beta=-b},{X,Y}),X*conique+Y*distance))); bl:= > denom(simplify(expand(subs(solve({alpha=b,beta=-b}, > {X,Y}),X*conique+Y*distance)))) > numer(simplify(expand(subs(solve({alpha=b,beta=-b}, > {X,Y}),X*conique+Y*distance))))\ > simplify(expand(bl*normale-numer(simplify(expand(subs( > solve({alpha=b,beta=-b},{X,Y}),X*conique+Y*distance)))))); bx2 + ((2 c — 2a)y + e — bu + 2av)x — by2 + (—2cu — d + bv)y — eu + dv (Xa + Y)x2 + (Xd + Xby-2Yu)x + (Xc + Y)y2 + (-2Yv + Xe)y-Yr2 + X f + Yu2 + Yv2 .(X c + Y)y2 + (-2Y v + X e)y -Yr2 + X f + Y u2 + Y v2, Xd + Xby-2Yu, Xa + Y a:=Xa + Y {Xa + Y)x2 + (X d - 2Y u)x -Y r2 + X f + Y u2 + Y v2, -2Yv + Xbx + Xe, Xc + Y i . 'y, y2 (3:=Xc + Y {Y'=J±±A,X = 2J^} a — c a — c b(ax2 + 2bxy + cy2 + 2dx + 2ey + 2 / - au2 + 2aux - av2 + 2avy - ay2 + ar2 — cu2 + 2cux — cx2 — cv2 + 2cv y + cr2)/(a — c) hi := a — c 168 b(ax2 + 2bxy + cy2 + 2dx + 2ey + 2 f — av? + 2aux — av2 + 2avy — ay2 + ar2 — cu2 + 2cux — cx2 — cv2 + 2cvy + cr 2 ) —3bcvy — bcux + bcv2 — bcr2 — bavy — 3baux — 2bey + bau2 + bav2 — bar2 + bcu2 — 2acyu + 4acyx — 2caxv + 2a2xv — 2b2 xy — 2b f — cex — cdv + cdy + 2c2yu — aeu — 2c2yx + aex — 2a2xy + adv — ady + ceu -2bdx > collect(distance,[x,y]); > collect(expand(X2*conique+Y2*normale),[x,y]) ; > coeffs(collect(expand(X2*conique+Y2*normale),[x,y]),x,'al');al; > alpha:=coeffs(collect(expand(X2*conique+Y2*normale),[x,y]),x)[3 > ];coeffs(collect(expand(X2*conique+Y2*normale),[y,x]),y,'a2'); > a2;beta:=coeffs(collect(expand(X2*conique+Y2*normale),[y,x]),y) > [3];solve({alpha=l,beta=l},{X2,Y2});simplify(expand(subs(solve( > {alpha=l,beta=l},{X2,Y2}),X2*conique+Y2*normale))); > b2:=denom(simplify(expand(subs(solve({alpha=l,beta=l}, > {X2,Y2}),X2*conique+Y2*normale)))); > numer(simplify(expand(subs(solve({alpha=l,beta=l},{X2.Y2}), > X2*conique+Y2*normale)))); > simplify(expand(b2*distance-numer(simplify(expand(subs(solve( > {alpha=l,beta=l},{X2,Y2}),X2*conique+Y2*normale)))))); u2 — 2 u x + x2 + v2 — 2 v y + y2 — r2 (Y2 b + X2 a) x2 + ( ( - 2 Y2 a + 2 Y2 c + X2 b) y - Y2 bu + Y2 e + 2 Y2 a v+ X2d)x + [X2 c - Y2 b)y2 + (X2 e - 2 Y2cu + Y2bv - Y2d)y + X2 f - Y2 eu + Y2dv (X2c - Y2 b)y2 + (X2 e - 2 Y2cu + Y2bv - Y2d)y + X2 f - Y2eu + Y2 dv, (-2Y2a + 2Y2c + X2b)y - Y2bu + Y2e + 2Y2av + X2d, Y2b + X2a a:= Y2b + X2a 169 (Y2 b + X2 a) x2 + (-Y2 bu + Y2 e + 2 Y2 av + X2 d) x + X2 f - Y2 eu + Y2dv, (-2Y2a + 2Y2c + X2b)x + X2e-2Y2cu+Y2bv-Y2d,X2c-Y2b f3:= X2c- Y2b {Y2 = - " ~ \ , X2 = 2 — } 1 b(a + c) a + cJ (bevy — bcux + bcx2 — bavy + baux + 2bey + bay2 + 2acyu — Aacyx i + 2caxv - 2a2xv + bcy2 + 2b2 xy -\-2b f + cex + cdv — cdy — 2c2yu + aeu + 2c2yx — aex + 2a2xy — adv + ady — ceu + bax2 + 2bdx)/(b (a + c)) • b2 := b (a + c) bevy — bcux + bcx2 — bavy + baux + 2bey + bay2 + 2acyu — Aacyx + 2.caxv — 2a2xv + bcy2+2b2xy + 2bf + cex + cdv — cdy — 2 c2 yu + aeu + 2c2yx — aex + 2 a2 xy — adv + ady — ceu + bax2 + 2bdx —3bcvy — bcux + bcv2 — bcr2 — bavy — 3baux — 2bey + bau2 + bav2 — bar2 + bcu2 — 2acyu + Aacy x — 2caxv + 2a2 xv — 2b2 xy — 2b f — cex — cdv + cdy + 2c2yu — aeu — 2c2yx + aex — 2a2xy + adv — ady + ceu -2bdx > read("/cs/beta/People/Anton/ag/toric/mapl2form"); > read("/cs/beta/People/Anton/ag/toric/maplib"); Warning, the p ro tec ted names norm and t r ace have been rede f ined and unprotected ., typelArray := proc(tarray, ttype) local i , size; i f not evalb(type(tarray, array)) t h e n R E T U R N (false) e n d i f ; i f not typ e(tarray, vector) t h e n RETURN( i r ae ) end i f ; size :— nops(convert(tarray, list)); for i to size do i f not evalb(type(iarray i, ttype)) t h e n R E T U R N (false) end i f e n d do; R E T U R N (true) end p r o c 170 typeMatrix := proc(tmatrix, ttype) l o c a H ; i f not evalb (type (tmatrix, matrix)) t h e n R E T U R N (false) e n d i f ; for i to rowdim(tmatrix) do i f not evalb(typelArray(row(imaira, i), ttype)) t h e n R E T U R N (false) e n d i f end do; R E T U R N (£rue) end p r o c > PList: = [expand(conique), > simplify(expand(bl*normale-numer(simplify(expand(subs( > solve({alpha=b,beta=-b},{X,Y}),X*conique+Y*distance)))))), > simplify(expand(b2*distance-numer(simplify(expand(subs( > solve({alpha=l,beta=l},{X2,Y2}), > X2*conique+Y2*normale))'))))] ;VList: = [x,y] ; PList := [ax2 + bxy + cy2 + dx + ey + f,—3bcvy — bcux + bcv2 — bcr2— bavy — Sbaux — 2bey + bau2 + bav2 — bar2 + bcu2 — 2acy u + 4acy x — 2caxv + 2a2 xv — 2b2 xy — 2b f — cex — cdv + cdy + 2c2 yu — aeu — 2c2yx + aex — 2a2xy + adv — ady + ceu — 2bdx, —2bcvy — 2bcux + bcv2 — bcr2 — 2bavy — 2baux + bau2 + bav2 — bar2 + bcu2 — Y2 ex + abx2 - X2bxy+ Y2bxu + cby2 + cbx2 + aby2 - X2 f + 2 Y2 cy u -2Y2cyx-2Y2axv + 2Y2axy - Y2byv- X2 ey - X2dx- X2 cy2 — X2 ax2 + Y2 dy - Y2 dv + Y2 by2 + Y2 eu - Y2bx2} VList := [x, y] > sys_forC(PList,VList,u, > Vcs/beta/People/Anton/ag/toric/offsetgeneralemtplus',2); writing for C program: 5 arguments; returns symb.coeffs Writing exponents of 3 polys in 2 vars to '/cs/beta/People/Anton/ag/toric/offsetgeneralemtplus.exps'. written all exponents Writing coeffile '/cs/beta/People/Anton/ag/toric/offsetgeneralemtplus.coef' with hidden 171 degree=2 3 polyns i n 2 v a r i a b l e s [x, y] v a r s l i s t and hidden u wrote monoms, c o e f f i l e ; r e t u r n symbolic coeffs {cSxl = -X2 f + bcu2 +bcv2 - bcr2 - bar2 + Y2 eu + bau2 + bav2 -Y2dv, c2x4 = -2a2 +Aac- 2c2 - 2b2, c2x3 — —Sbcv — 2be + 2c2u — bav — ad — 2acu + cd, clx6 = b, c2xl — bav2 + ceu + bcv2 — bcr2 + adv — cdv + bcu2 + bau2 — 2 b f — bar2 — aeu, c2x2 = — bcu — 2bd — 2cav — 3bau + ae — ce + 2a2 v, clx4 = c, clx5 — a, clx2 = d, clxS — e, clxl — f, c3x5 = ab - Y2b + cb - X2 a, c3x6 = 2Y2a-X2b-2Y2c, c3x4 = -X2c + ab+ Y2b + cb, c3x2 = Y2bu-2bcu - Y2 e - 2 Y2 av - 2bau - X2 d, c3x3 = Y2d-2bcv + 2 Y2cu- X2e- Y2bv - 2bav) B.2 The M a x i m a program for computing the sparse res-ultant c3x2 : a - c ; c3x3 : -a*r~2+a*u~2+a*v~2-f; c3x4 : - 2 * a * u - d ; c3x5 : - 2 * a * v - e ; c l x 2 : a ; c l x 3 : c ; c l x 4 : f ; c l x 5 : d ; c l x 6 : e; c 2 x l : b~2-2*a*c+2*a"2; c2x2 : b*c+a*b; c2x3 : -a*d*v+b*f+a*e*u; c2x4 : b*d+a*b*u-2*a"2*v-a*e; 172 c2x5 : b*e-a*b*v+2*a*c*u+a*d; c 3 x l : - b ; c l x l : b ; M : m a t r i x ( c l x l , c l x 2 , c l x 3 , c l x 4 , c l x 5 , c l x 6 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ] , c l x 2 , 0 , c l x l , 0 , 0 , c l x 5 , 0 , 0 , 0 , c l x 3 , c l x 4 , c l x 6 , 0 ] , c 2 x l , 0 , c 2 x 2 , c 2 x 3 , c 2 x 4 , c 2 x 5 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ] , 0 , 0 , 0 , c l x 6 , c l x l , c l x 3 , c l x 2 , c l x 4 , c l x 5 , 0 , 0 , 0 , 0 ] , c 2 x 2 , c 2 x l , 0 , 0 , c 2 x 5 , 0 , c 2 x 4 , 0 , c 2 x 3 , 0 , 0 , 0 , 0 ] , 0 , 0 , 0 , c 2 x 5 , c 2 x l , c 2 x 2 , 0 , c 2 x 3 , c 2 x 4 , 0 , 0 , 0 , 0 ] , 0 , 0 , c 2 x l , 0 , 0 , c 2 x 4 , 0 , 0 , 0 , c 2 x 2 , c 2 x 3 , c 2 x 5 , 0 ] , c 3 x l , 0 , c 3 x 2 , c 3 x 3 , c 3 x 4 , c 3 x 5 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ] , c 3 x 2 , c 3 x l , 0 , 0 , c 3 x 5 , 0 , c 3 x 4 , 0 , c 3 x 3 , 0 , 0 , 0 , 0 ] , 0 , 0 , 0 , c 3 x 5 , c 3 x l , c 3 x 2 , 0 , c 3 x 3 , c 3 x 4 , 0 , 0 , 0 , 0 ] , 0 , 0 , c 3 x l , 0 , 0 , c 3 x 4 , 0 , 0 , 0 , c 3 x 2 , c 3 x 3 , c 3 x 5 , 0 ] , 0 , 0 , 0 , c l x 5 , c l x 2 , c l x l , 0 , 0 , 0 , 0 , c l x 6 , c l x 3 , c l x 4 ] , 0 , 0 , 0 , c 2 x 4 , 0 , c 2 x l , 0 , 0 , 0 , 0 , c 2 x 5 , c 2 x 2 , c2x3] p : d e t e r m i n a n t ( M ) ; r : f a c t o r ( p ) ; B.3 A n implicit equation for parabolas A n implicit equation of a parabola in a system centred on its summit (the intersec-tion of its axis of symmetry with itself) and the x axis being its axis of symmetry is: y2 — 2px = 0. A n implicit equation of the r—generalised offset to a parabola is: 4 * y 6 + 4 * x2 * y 4 — 20*p*x * y 4 — (12 * r 2 — p 2 ) * y 4 - 16*p*x 3 * y 2 — 8* (r 2 — 4*p 2 ) * x 2 *y2 + 173 4* (p* r2 — p 3 ) *x * y2 + 2 * (6 * r 4 — 10 *p 2 * r 2 ) * y 2 + 16 *p 2 * x 4 — 16 * (p* r 2 + p 3 ) * x 3 + 4 * ( r 4 — 2 * p 2 * r 2 + p 4 ) * x 2 + 16*(p*r 4 +p 3 *r 2 )*x — (4*r 6 + 8*p 2 *r 4 + 4*p 4 *r 2 ) . B.4 A n implicit equation for ellipses, circles or hyper-bolas 2 2 A n implicit equation of the conic is: ^ + ±p — 1 = 0 where ± stands for + in the case of an ellipse or a circle and — in the case of an hyperbola. A n implicit equation of the r—generalised offset to an ellipse, a circle or an hyperbola is: c 4 * u 8 + 2 * ( c 4 + c 3 ) * u 2 * 7 j 6 - 2 * ( ( 2 * c 4 - c 3 ) * r 2 + ( c 4 - 2 * c 3 ) * e ) * 7 j 6 + (c 4 + 4 * c 3 + c 2 ) * u 4 * u 4 - 2 * ( ( 3 * c 4 - c 3 + c 2 ) * r 2 - ( c 4 - c 3 + 3 * c 2 ) * e ) * t i 2 * 7 j 4 + ( (6*c 4 -6*c 3 + c 2 )*r 4 +2*(3*c 4 -5*c 3 +3*c 2 )*e*r 2 +(c 4 -6*c 3 +6*c 2 )*e 2 )*v 4 +2*(c 3 +c 2 )*u 6 *ij 2 -2* ( ( c 4 - c 3 +3*c 2 )* r 2 - (3*c 3 - c 2 +c)*e )*u 4 * t ; 2 +2*( (3*c 4 -5*c 3 +3*c 2 )* r 4 - (2*c 4 -3* c 3 -3*c 2 +2*c)*e*r 2 +(3*c 3 -5*c 2 +3*c)*e 2 )*u 2 *7j 2 -2*((2*c 4 -3*c 3 +c 2 )*r 6 +(3*c 4 -4*c 3+2*c 2 —c)*e*r4+(c4—2*c3+4*c2 —3*c)*e2*r2 — (c3—3*c2+2*c)*e3)*t>2+c2*'u8+ 2*((c 3 -2*c 2 )*r 2 + (2*c 2 -c)*e)*u 6 +((c 4 -6*c 3 +6*c 2 )*r 4 +2*(3*c 3 -5*c 2 +3*c)*e* r 2 +(6*c 2 -6*c+l )*e 2 )*u 4 -2*( (c 4 -3*c 3 +2*c 2 )*r 6 - (c 4 -2*c 3 +4*c 2 -3*c)*e*r 4 - (3* c 3 -4*c 2 +2*c - l )*e 2 *r 2 - (2*c 2 -3*c+ l )*e 3 )*u 2 +( (c 4 -2*c 3 +c 2 )* r 8 +2*(c 4 - c 3 - c 2 + c)*e*r 6 + ( c 4 +2*c 3 -6*c 2 +2*c+l )*e 2 *r 4 +2*(c 3 - c 2 - c+ l )*e 3 *r 2 + (c 2 -2*c+l )*e 4 ) where c = ± | £ and e = —a2 and ± stands for + in the case of an ellipse or a circle, and — in the case of an hyperbola. 174 Appendix C The Singular program for obtaining the matrix of the sparse resultant of the Delaunay graph conflict locator for conies C l The case of parabolas r i n g p r e d i c a t r = ( 0 , a , b ) c , d , e , f , g , h ) i , j , k , l , m , n , o , p ) u 0 , u l , u 2 , u 3 , u 4 ) , ( x . y . r . s ) , d p ; p o l y offse t l=4*y~6+4*x~2*y~4-20*p*x*y~4-12*r*y~4+p~2*y~4 -16*p*x~3*y~2-8*r*x~2*y~2+32*p~2*x~2*y~2+4*p*r*x*y~2 -4*p"3*x*y~2+12*r~2*y~2-20*p~2*r*y~2+16*p~2*x~4-16* p*r*x~3-16*p~3*x~3+4*r"2*x~2^-8*p"2*r*x~2+4*p~4*x~2 +16*p*r~2*x+16*p~3*r*x-4*r~3-8*p"2*r~2-4*p~4*r; p o l y offset2=4*(b*x+a*y+d)"6+4*(a*x-b*y+c)"2*(b*x+a*y+d)"4 -20*o*(a*x-b*y+c)*(b*x+a*y+d)~4-12*r*(b*x+a*y+d)~4+o~2* (b*x+a*y+d)~4-16*o*(a*x-b*y+c)"3*(b*x+a*y+d)"2 175 -8*r*(a*x-b*y+c)~2*(b*x+a*y+d)~2+32*o~2*(a*x-b*y+c)"2* (b*x+a*y+d)~2+4*o*r*(a*x-b*y+c)*(b*x+a*y+d)"2 -4*o~3*(a*x-b*y+c)*(b*x+a*y+d)~2+12*r~2*(b*x+a*y+d)"2 -20*o~2*r*(b*x+a*y+d)~2+16*o~2*(a*x-b*y+c)"4 -16*o*r*(a*x-b*y+c)~3-16*o~3*(a*x-b*y+c)"3 +4* r ~2*(a*x-b*y+c)~2-8*o~2*r*(a*x-b*y+c)"2 +4*o~4*(a*x-b*y+c)~2+16*o*r~2*(a*x-b*y+c)+ 16*o~3*r*(a*x-b*y+c) -4*r~3-8*o~2*r~2-4*o~4*r ; p o l y of fse t3=4*(f*x+e*y+h)"6+4*(e*x-f*y+g)"2*(f*x+e*y+h)"4 -20*n*(e*x-f*y+g)*( f*x+e*y+h)~4-12*r*( f*x+e*y+h)"4 +n~2*(f*x+e*y+h)~4-16*n*(e*x-f*y+g)"3*(f*x+e*y+h)"2 -8*r*(e*x-f*y+g)~2*(f*x+e*y+h)~2+32*n~2*(e*x-f*y+g)"2* (f*x+e*y+h)~2+4*n*r*(e*x-f*y+g)*(f*x+e*y+h)"2 -4*n~3*(e*x-f*y+g)*(f*x+e*y+h)~2+12*r~2*(f*x+e*y+h)"2 -20*rT2*r*( f*x+e*y+h)~2+16*iT2*(e*x- f*y+g)"4 -16*n*r* (e*x- f*y+g)~3-16*n~3*(e*x- f*y+g)"3 +4*r~2*(e*x- f*y+g)~2-8* rT2*r* (e*x- f*y+g)"2 +4*n~4*(e*x-f*y+g)"2+16*n*r~2*(e*x-f*y+g)+ 16*n~3*r* (e*x - f*y+g) -4* r~3 -8* rT2*r~2 -4*n~4*r ; p o l y o f f s e t 4 = 4 * ( j * x + i * y + l ) " 6 + 4 * ( i * x - j * y + k ) ~ 2 * ( j * x + i * y + l ) " 4 - 2 0 * m * ( i * x - j *y+k)*( j * x + i * y + l ) ~ 4 - 1 2 * s * ( j * x + i * y + l ) " 4 +m"2*(j * x + i * y + l ) ~ 4 - 1 6 * m * ( i * x - j *y+k)"3*( j * x + i * y + l ) " 2 - 8 * s * ( i * x - j *y+k)"2*( j *x+ i*y+ l )~2+32*nT2*( i*x - j*y+k)~2* ( j *x+ i*y+ l ) "2+4*m*s* ( i*x - j *y+k)*( j * x + i * y + l ) " 2 - 4 * m ~ 3 * ( i * x - j * y + k ) * ( j *x+i*y+l )"2+12*s"2*( j * x + i * y + l ) " 2 -20*m~2*s*( j*x+i*y+l )"2+16*m"2*( i*x- j*y+k)~4 - 1 6 * m * s * ( i * x - j * y + k ) ~ 3 - 1 6 * n T 3 * ( i * x - j * y + k ) " 3 + 4 * s ~ 2 * ( i * x - j * y + k ) " 2 - 8 * m " 2 * s * ( i * x - j *y+k)~2 176 +4*nT4*(i*x-j*y+k)"2+16*m*s~2*(i*x-j*y+k)+ 16*m~3*s*(i*x-j*y+k)-4*s~3-8*m~2*s~2-4*nT4*s; ideal predicati=u0+ul*x+u2*y+u3*r+u4*s,offsetl,offset2-offsetl, offset3-offsetl, offset4-(i~6+j~2*i~4)*offsetl; module predicatm=mpresmat(predicati,0); C.2 The case of ellipses and/or hyperbolas ring predicatr = (0,a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,t,u,u0,ul,u2, u3,u4,v), (x.y.r.s),(dp(4),C); poly offsetl=a~4*y~8+2*(a~4+a~3)*x~2*y~6-2*((2*a~4-a~3)*r +(a~4-2*a~3)*b)*y~6+(a~4+4*a~3+a~2)*x"4*y~4 -2*((3*a~4-a~3+a~2)*r-(a~4-a~3+3*a~2)*b)*x~2*y~4 +((6*a~4-6*a~3+a~2)*r~2+2*(3*a~4-5*a~3+3*a~2)*b*r +(a~4-6*a~3+6*a~2)*b"2)*y~4+2*(a~3+a~2)*x~6*y"2 -2*((a~4-a~3+3*a"2)*r-(3*a~3-a~2+a)*b)*x~4*y~2 +2*((3*a~4-5*a~3+3*a~2)*r~2-(2*a~4-3*a~3-3*a~2+2*a)*b*r +(3*a~3-5*a~2+3*a)*b~2)*x~2*y~2-2*((2*a~4-3*a~3+a~2)*r~3 +(3*a~4-4*a~3+2*a~2-a)*b*r~2+(a~4-2*a~3+4*a~2-3*a)*b~2*r -(a~3-3*a~2+2*a)*b~3)*y~2+a~2*x~8+2*((a~3-2*a~2)*r +(2*a"2-a)*b)*x"6+((a~4-6*a~3+6*a~2)*r~2 +2*(3*a~3-5*a~2+3*a)*b*r+(6*a~2-6*a+l)*b"2)*x~4 -2*((a~4-3*a"3+2*a"2)*r"3-(a~4-2*a*3+4*a"2-3*a)*b*r"2 - (3*a~3-4*a~2+2*a-l) *b~2*r- (2*a~2-3*a+l)*b~3) *x~2 +((a"4-2*a"3+a"2)*r~4+2*(a~4-a~3-a~2+a)*b*r~3 +(a~4+2*a~3-6*a~2+2*a+l)*b~2*r~2+2*(a~3-a~2-a+l)*b~3*r +(a~2-2*a+l)*b~4); 177 p o l y o f f se t2=c~4*(f*x+e*y+h)"8+2*(c~4+c"3)*(e*x- f*y+g)"2* ( f*x+e*y+h)"6 -2*( (2*c"4 -c"3 )* r+(c"4 -2*c"3 )*d)* ( f*x+e*y+h)"6+ (c"4+4*c~3+c"2)* (e*x- f*y+g)"4*( f*x+e*y+h)"4 -2*( (3*c"4-c"3+c~2)*r -(c"4-c~3+3*c"2)*d)*(e*x- f*y+g)"2*( f*x+e*y+h)~4+( (6*c"4-6*c"3+c"2)*r"2 +2*(3*c"4-5*c"3+3*c"2)*d*r+(c~4-6*c"3+6*c"2)*d"2)*(f*x+e*y+h)"4 +2*(c"3+c"2 )* (e*x- f*y+g)"6*( f*x+e*y+h)"2 -2* ( (c"4 -c"3+3*c"2 )* r - (3*c"3 -c"2+c )*d )* (e*x - f*y+g)"4* ( f*x+e*y+h)"2 +2* ( (3*c" 4 -5*c" 3+ 3*c"2 )* r "2 - (2*c"4 -3*c"3 -3*c"2+2*c )*d* r + (3*c"3-5*c~2+3*c)*d~2)*(e*x-f*y+g)"2*(f*x+e*y+h)"2 - 2 * ( ( 2 * c " 4 - 3 * c " 3 + c " 2 ) * r " 3 + ( 3 * c " 4 - 4 * c ~ 3 + 2 * c " 2 - c ) * d * r " 2 + (c"4 -2*c"3+4*c"2 -3*c )*d"2*r - ( c"3 -3*c"2+2*c)*d"3 )* ( f*x+e*y+h)"2 + c " 2 * ( e * x - f * y + g ) " 8 + 2 * ( ( c " 3 - 2 * c " 2 ) * r + ( 2 * c ~ 2 - c ) * d ) * ( e * x - f * y + g ) " 6 +( (c~4-6*c"3+6*c"2)*r~2+2*(3*c"3-5*c"2+3*c)*d*r+(6*c"2-6*c+l )*d"2)* ( e * x - f * y + g ) " 4 - 2 * ( ( c " 4 - 3 * c " 3 + 2 * c " 2 ) * r " 3 - ( c " 4 - 2 * c " 3 + 4 * c " 2 - 3 * c ) * d * r " 2 - ( 3 * c ~ 3 - 4 * c " 2 + 2 * c - l ) * d ~ 2 * r - ( 2 * c " 2 - 3 * c + l ) * d " 3 ) * ( e * x - f * y + g ) " 2 +( (c"4-2*c~3+c"2)*r~4+2*(c~4-c"3-c~2+c)*d*r"3+ (c"4+2*c"3-6*c"2+2*c+ l )*d"2*r"2+2*(c"3 -c"2 -c+ l )*d"3*r+ ( c " 2 - 2 * c + l ) * d " 4 ) ; p o l y o f f s e t 3 = i " 4 * ( l * x + k * y + n ) " 8 + 2 * ( i ~ 4 + i " 3 ) * ( k * x - l * y + m ) " 2 * ( l * x + k * y + n ) ~ 6 - 2 * ( ( 2 * i " 4 - i " 3 ) * r + ( i " 4 - 2 * i " 3 ) * j ) * ( l * x + k * y + n ) " 6 + ( i " 4 + 4 * i " 3 + i " 2 ) * ( k * x - l * y + m ) " 4 * ( l * x + k * y + n ) " 4 - 2 * ( ( 3 * i " 4 - i " 3 + i ~ 2 ) * r - ( i " 4 - i " 3 + 3 * i " 2 ) * j ) * ( k * x - l * y + m ) " 2 * ( l * x + k * y + n ) " 4 + ( ( 6 * i " 4 - 6 * i " 3 + i " 2 ) * r " 2 + 2 * ( 3 * i " 4 - 5 * i " 3 + 3 * i ~ 2 ) * j * r + ( i ~ 4 - 6 * i " 3 + 6 * i " 2 ) * j " 2 ) * ( l * x + k * y + n ) " 4 + 2 * ( i " 3 + i " 2 ) * ( k * x - l * y + m ) " 6 * ( l * x + k * y + n ) " 2 - 2 * ( ( i " 4 - i " 3 + 3 * i " 2 ) * r - ( 3 * i " 3 - i " 2 + i ) * j ) * ( k * x - l * y + m ) " 4 * ( l * x + k * y + n ) " 2 + 2 * ( ( 3 * i " 4 - 5 * i " 3 + 3 * i ~ 2 ) * r " 2 - ( 2 * i " 4 - 3 * i " 3 - 3 * i " 2 + 2 * i ) * j * r + ( 3 * i ~ 3 - 5 * i ~ 2 + 3 * i ) * j " 2 ) * ( k * x - l * y + m ) " 2 * ( l * x + k * y + n ) " 2 - 2 * ( ( 2 * i " 4 - 3 * i " 3 + i " 2 ) * r " 3 + ( 3 * i " 4 - 4 * i " 3 + 2 * i " 2 - i ) * j * r ~ 2 178 + ( i " 4 - 2 * i " 3 + 4 * i " 2 - 3 * i ) * j ~ 2 * r - ( i " 3 - 3 * i " 2 + 2 * i ) * j " 3 ) * ( l * x + k * y + n ) " 2 + i " 2 * ( k * x - l * y + m ) " 8 + 2 * ( ( i " 3 - 2 * i " 2 ) * r + ( 2 * i " 2 - i ) * j ) * ( k * x - l * y + m ) " 6 + ( ( i " 4 - 6 * i " 3 + 6 * i " 2 ) * r " 2 + 2 * ( 3 * i " 3 - 5 * i " 2 + 3 * i ) * j * r + ( 6 * i " 2 - 6 * i + l ) * j ~ 2 ) ( k * x - l * y + m ) " 4 - 2 * ( ( i " 4 - 3 * i " 3 + 2 * i " 2 ) * r " 3 - ( i ~ 4 - 2 * i " 3 + 4 * i " 2 - 3 * i ) * j * r " 2 - ( 3 * i " 3 - 4 * i ~ 2 + 2 * i - l ) * j " 2 * r - ( 2 * i " 2 - 3 * i + l ) * j " 3 ) * ( k * x - l * y + m ) " 2 + ( ( i ~ 4 - 2 * i " 3 + i ~ 2 ) * r " 4 + 2 * ( i " 4 - i ~ 3 - i " 2 + i ) * j * r " 3 + ( i " 4 + 2 * i " 3 - 6 * i ~ 2 + 2 * i + l ) * j " 2 * r ~ 2 + 2 * ( i " 3 - i " 2 - i + l ) * j " 3 * r + ( i " 2 - 2 * i + l ) * j " 4 ) ; p o l y of fse t4=o"4*( t*x+q*y+v)"8+2*(o~4+o"3)*(q*x- t*y+u)"2* ( t*x+q*y+v)"6-2*( (2*o"4-o"3)*s+(o"4-2*o"3)*p)*( t*x+q*y+v)"6 +(o"4+4*o"3+o"2)*(q*x- t*y+u)"4*( t*x+q*y+v)"4-2*((3*o~4-o"3+o"2)*s - (o~4-o~3+3*o~2)*p)*(q*x- t*y+u)"2*( t*x+q*y+v)"4+((6*o~4-6*o~3+o"2) s"2+2*(3*o~4-5*o"3+3*o"2)*p*s+(o"4-6*o"3+6*o~2)*p"2)*(t*x+q*y+v)~4 +2*(o"3+o"2)*(q*x- t*y+u)"6*( t*x+q*y+v)"2-2*( (o"4-o"3+3*o"2)*s - (3*o"3-o"2+o)*p)*(q*x- t*y+u)"4*( t*x+q*y+v)"2+ 2*((3*o"4-5*o~3+3*o"2)*s"2-(2*o"4-3*o"3-3*o~2+2*o)*p*s +(3*o"3-5*o"2+3*o)*p"2)*(q*x- t*y+u)"2*( t*x+q*y+v)"2 -2*( (2*o"4-3*o"3+o"2)*s"3+(3*o"4-4*o"3+2*o"2-o)*p*s"2 +(o"4-2*o"3+4*o~2-3*o)*p"2*s-(o"3-3*o"2+2*o)*p"3)*( t*x+q*y+v)~2 +o"2*(q*x- t*y+u)"8+2*( (o"3-2*o"2)*s+(2*o"2-o )*p)* (q*x- t*y+u)"6 +((o"4-6*o"3+6*o"2)*s"2+2*(3*o"3-5*o"2+3*o)*p*s + (6*o"2 -6*o+ l )*p"2 )* (q*x- t*y+u)"4 -2* ( (o"4 -3*o"3+2*o"2)*s"3 - (o"4 -2*o"3+4*o"2 -3*o )*p*s"2 - (3*o"3 -4*o"2+2*o- l )*p"2*s - (2*o"2-3*o+l )*p"3)* (q*x- t*y+u)"2+( (o"4-2*o~3+o"2) *s"4 +2*(o"4-o"3-o"2+o)*p*s"3+(o"4+2*o"3-6*o"2+2*o+l)*p"2*s"2 +2*(o~3-o"2-o+l )*p"3*s+(o"2-2*o+l )*p~4) ; i d e a l p i d e a l = u 0 + u l * x + u 2 * y + u 3 * r + u 4 * s , o f f s e t l , ( a " 4 - 2 * a " 3 + a " 2 ) * o f f s e t 2 - ( c " 4 - 2 * c " 3 + c " 2 ) * o f f s e t l , 179 ( a ~ 4 - 2 * a ~ 3 + a ~ 2 ) * o f f s e t 3 - ( i ~ 4 - 2 * i ~ 3 + i ~ 2 ) * o f f s e t l , (a"2)*offse t4- (o"4*q~4*t"4+2*o~4*q"2*t"6+o"4*t"8+2*o~3*q~6*t~2 +4*o~3*q~4*t"4+2*o~3*q"2*t'*6+o"2*q~8+2*o~2*q""6*t"2+o"2*q~4*t~4)* o f f s e t l ; module pmodu le=mpresma t (p idea l , 0 ) ; 180 Appendix D A Maple program for the ALIAS based Delaunay graph conflict locator for semi-algebraic sets After having specified the equations and inequality, we need to specify the set of equations (EQ) and the set of unknowns (VAR) of the system of equations and inequalities. Then, we need to specify the search intervals for all the variables (IN). The Maple program must specify where the A L I A S Maple library is located (by set-ting the Iibname variable if it is not located in the standard Maple library location), where the A L I A S C + + library is located (by setting the A L I A S / L i b variable) as well as where the BIAS/Prof i l C + + library is located (by setting the A L I A S / P r o f i l vari-able). Then, some A L I A S parameters (prefixed with ALIAS) can be specified (de-bugging: ALIAS/debug, bisection behaviour: ALIAS/single_bisection, optimisation: ALIAS/optimised, 3B parameters: A L I A S / 3 B , A L I A S / M a x 3 B , ALIAS/Del ta3B) . Then, we issue the parser command that will convert the previous code into C++ code that uses the C + + A L I A S library (for example GradientSolve in the case of 181 the Maple program shown here). D.l The system without the implicit equations of the generalised offsets OA := subs(alpha=-3,beta=2,xl"2/beta+yl"2/(beta/alpha)-l); OB := subs(u=5*x2-3*y2+7,v=2*x2+3*y2+4,alpha=5,beta=4, u~2/beta+v~2/(beta/alpha)-1); OC := subs(u=3*x3-2*y3+l,v=3*x3+2*y3+4,alpha=2,beta=l, u~2/beta+v~2/(beta/alpha)-1); NA := -diff(OA,yl)*(x-xl)+diff(OA.xl)*(y-yl); NB := -diff(0B,y2)*(x-x2)+diff(0B,x2)*(y-y2); NC := -diff(0C,y3)*(x-x3)+diff(0C,x3)*(y-y3); DA := (x-xl)"2+(y-yl)~2-r"2; DB := (x-x2)~2+(y-y2)~2 -r~2; DC ':= (x-x3)~2+(y-y3)~2 - r"2; OD := subs(u=5*x4-2*y4+l,v=x4+4*y4+2,alpha=3,beta=2, u"2/beta+v"2/(beta/alpha)-1); ND := -diff(0D,y4)*(x-x4)+diff(0D,x4)*(y-y4); PVD:=(x-x4)"2+(y-y4)"2-r"2<0; TOA:=subs(xl=xtl,yl=ytl,OA); TOB:=subs(x2=xt2,y2=yt2,OB); TOC:=subs(x3=xt3,y3=yt3,OC); PVA:=(x-xtl)"2+(y-ytl)"2-r"2>=or TNA := -diff(TOA,ytl)*(x-xtl)+diff(TOA,xtl)*(y-ytl); TNB := -diff(TOB,yt2)*(x-xt2)+diff(TOB,xt2)*(y-yt2); TNC := -diff(TOC,yt3)*(x-xt3)+diff(TOC,xt3)*(y-yt3); PVB:=(x-xt2)"2+(y-yt2)-2-r"2>=0; 182 PVC:=(x-xt3)~2+(y-yt3)"2-r"2>=0; EQ:=[OA,OB,OC,NA,NB,NC,DA,DB,DC,TOA,TOB,TOC,TNA,TNB,TNC, PVA,PVB,PVC,OD,ND,PVD]; VAR:=[xl,yl,x2,y2,x3,y3,x,y,xtl,ytl,xt2,yt2,xt3,yt3,x4,y4,r]; libname:="/cs/beta/People/Anton/ag/ALIAS/ALIAS/maple",libname; with(ALIAS); 'ALIAS/lib':="/cs/beta/People/Anton/ag/ALIAS/ALIAS/Lib-solaris"; 'ALIAS/profil':="/cs/beta/People/Anton/ag/ALIAS/ALIAS/Profil"; IN:-[[-1000,1000]]; for i from 1 to 15 do IN: = [op(IN), [-1000,1000]]: od; IN:=[op(IN),[0,1000]]; 'ALIAS/debug' := 1; 'ALIAS/single_bisection':=2; 'ALIAS/optimized':=0; 'ALIAS/3B':=0; 'ALIAS/SubEq3B': = [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1]; 'ALIAS/Max3B':=2000; 'ALIAS/Delta3B':=0.1; 'ALIAS/maxkraw':=30; GradientSolve(EQ,VAR,IN); D.2 The system with the implicit equations of the gen-eralised offsets offsete:=alpha~4*v~8+2*alpha~4*u~2*v~6+ 2*alpha~3*iT2*v~6-4*alpha~4*r*v~6 +2*alpha"3*r*v"6-2*alpha~4*beta*v"6+4*alpb.a"3*beta*v"6 183 +alpha~4*u~4*v~4+4*alpha~3*u~4*v~4+alpha~2*u~4*vM -6*alpha~4*r*u"2*v~4+2*alpha~3*r*u~2*v~4 -2*alpha~2*r*u~2*v~4+2*alpha~4*beta*u"2*v~4 -2*alpha~3*beta*u~2*v~4+6*alpha~2*beta*u~2*v~4 +6*alpha~4*r~2*v~4-6*alpha~3*r"2*v~4+alpha~2*r~2*v~4 +6*alpha~4*beta*r*v~4-10*alpha~3*beta*r*v~4 +6*alpha~2*beta*r*v"4+alpha~4*beta~2*v~4 -6*alpha~3*beta~2*v~4+6*alpha~2*beta~2*v~4 +2*alpha~3*u."6*v"2+2*alpha"2*u~6*v~2 -2*a lpha"4*r*u"4*v"2+2*a lpha"3*r*u~4*v~2 -6*alpha~2*r*u~4*v~2+6*alpha~3*beta*u~4*v~2 -2*alpha~2*beta*u~4*v~2+2*alpha*beta*u~4*v~2 +6*alpha~4*r~2*u~2*v~2-10*alpha~3*r~2*u~2*v~2 +6*alpha~2*r~2*u~2*v~2-4*alpha~4*beta*r*u~2*v~2 +6*alpha~3*beta*r*u"2*v"2+6*alpha~2*beta*r*u"2*v"2 -4*alpha*beta*r*u~2*v~2+6*alpha~3*beta~2*u~2*v~2 -10*alpha~2*beta~2*u~2*v~2+6*alpha*beta~2*u~2*v~2 -4*a lpha"4* r "3*v"2+6*a lpha~3*r~3*v"2 -2*a lpha~2*r"3*v~2 -6*alpha~4*beta*r~2*v~2+8*alpha~3*beta*r~2*v~2 -4*alpha~2*beta*r~2*v~2+2*alpha*beta*r~2*v~2 -2*alpha~4*beta~2*r*v~2+4*alpha~3*beta~2*r*v~2 -8*alpha~2*beta~2*r*v~2+6*alpha*beta~2*r*v~2 +2*alpha~3*beta~3*v~2-6*alpha~2*beta~3*v~2 +4*alpha*beta~3*v~2+alpha~2*u~8+2*alpha~3*r*iT6 -4*a lpha"2*r*u~6+4*a lpha"2*be ta*u"6-2*a lpha*be ta*u"6 +alpha~4*r"2*u~4-6*alpha"3*r~2*u~4+6*alpha"2*r"2*u~4 +6*alpha~3*beta*r*u~4-10*alpha~2*beta*r*u~4 +6*alpha*beta*r*u~4+6*alpha~2*beta~2*u~4 184 -6*alpha*beta~2*u~4+beta~2*u~4-2*alpha~4*r~3*u~2 +6*alpha~3*r~3*u~2-4*alpha~2*r~3*u~2 +2*alpha~4*beta*r~2*u~2-4*alpha~3*beta*r~2*u~2 +8*alpha~2*beta*r~2*u~2-6*alpha*beta*r~2*u~2 +6*alpha~3*beta~2*r*u~2-8*alpha~2*beta~2*r*u~2 +4*alpha*beta"2*r*u"2-2*beta"2*r*u"2+4*alpha"2*beta"3*u"2 -6*alpha*beta"3*u~2+2*beta"3*u"2+alpha"4*r"4-2*alpha"3*r~4 +alpha~2*r~4+2*alpha~4*beta*r~3-2*alpha~3*beta*r~3 -2*alpha"2*beta*r"3+2*alpha*beta*r"3+alpha"4*beta~2*r"2 +2*alpha"3*beta~2*r"2-6*alpha"2*beta"2*r"2+2*alpha*beta~2*r"2 +beta~2*r~2+2*alpha~3*beta~3*r-2*alpha~2*beta~3*r -2*alpha*beta~3*r+2*beta~3*r+alpha~2*beta~4 -2*alpha*beta"4+beta"4; offsetl:=simplify(subs(alpha=3,beta=-2,u=x,v=y,offsete)); offset2:=simplify(subs(u=5*x-3*y+7,v=2*x+3*y+4,alpha=5, beta=-4,offsete)); offset3:=simplify(subs(u=3*x-2*y+l,v=3*x+2*y+4,alpha=2, beta=-l,offsete)); offset4:=simplify(subs(u=5*x-2*y+l,v=x+4*y+2,alpha=3, beta=-2,r=R4,offsete)); offsetlv:=simplify(subs(r=Rl,offsetl)); offset2v:=subs(r=R2,offset2); offset3v:=subs(r=R3,offset3); t1:=Rl-r>=0;t2:=R2-r>=0;t3:=R3-r>=0;t4:=R4-r<0; EQ: = [offset1,offset2,offset3,offset4,offsetlv,offset2v, offset3v,tl,t2,t3 )t4]; VAR: = [x)y,r,Rl,R2,R3,R4] ; 185 libname:="/cs/beta/People/Anton/ag/ALIAS/ALIAS/maple".libname; with(ALIAS): 'ALIAS/lib': = M/cs/beta/People/Anton/ag/ALIAS/ALIAS/Lib-solaris": 'ALIAS/profil':=M/cs/beta/People/Anton/ag/ALIAS/ALIAS/Profil"; 'ALIAS/lib': ="/cs/beta/People/Anton/ag/ALIAS/ALIAS/Lib-solaris»; IN: = [[-1000,1000], [-1000,1000],[0,1000],[0,1000],[0,1000], [0,1000], [0,1000]]; •ALIAS/debug':=1; 'ALIAS/single_bisection':=2; 'ALIAS/optimized':=1; 'ALIAS/3B':=1; <ALIAS/Max3B':=2000; 'ALIAS/Delta3B':=0.1; 'ALIAS/maxkraw':=30; 'ALIAS/SubEq3B':=[0,0,0,0,0,0,0,1,1,1,1]; GradientSolve(EQ,VAR,IN); 186 Index algebraic variety, 1, 37, 146 algebraically closed field, 49, 68, 82, 93, 118 A L I A S , xvii , 21-23, 62, 76, 80, 135, 138, 139, 141, 142, 148-150, 181 Apollonius circles, 42, 44, 54-58 Apollonius tenth problem, 41-43 Delaunay graph, 4, 6-8, 10, 12, 14, 15, 17, 19, 22, 37-45, 52-56, 60, 61, 66, 72, 82, 88, 120, 121, 123-126, 135, 137, 147, 150 Delaunay graph conflict locator, 8, 12-14, 17, 19-23, 38-40, 42, 43, 45, 48, 51-54, 56, 59, 62, 63, 65-67, 74, 76-78, 106, 119-121, 123-126, 131-135, 138, 139, 142, 143, 146-150, 164, 175, 181 Delaunay graph predicate, 8, 12, 43 field, 68, 81 generalised offset, 17, 19-21, 23, 38, 72, 81, 82, 84, 85, 89, 90, 98, 100, 102-108, 114, 117-119, 121, 125-131, 135, 142, 146, 147, 149, 150, 167, 174 generalised Voronoi diagrams, 10 generalised Voronoi vertex, 20, 38, 119, 121-123, 140, 149 Grobner basis, 39, 40, 50, 60, 64-66, 81, 148-150 Hermite, 32, 34, 35 Newton-Raphson, 11, 34, 78 offset, 17 proximal region, 2 semi-algebraic set, 1, 2, 9, 12, 13, 15, 19-23, 37, 76, 88, 121, 135, 136, 138, 140, 141, 147, 148, 150 sparse resultant, 19, 21, 22, 39, 67, 68, 70-74, 82, 106, 107, 110, 112, 117, 131, 132, 134, 147-149, 167, 172, 175 true offset, 81, 82, 84, 88, 89 true Voronoi vertex, 15, 121, 125, 138 .Voronoi diagram, 2-4, 6, 10-12, 19, 20, 22, 24, 25, 28-32, 34, 37, 41, 119-121 Voronoi diagram for planar domains with curved boundaries, 29 187 Voronoi diagrams for curved objects, 28 Voronoi diagrams of general manifolds, 24 188
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Voronoi diagrams of semi-algebraic sets
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Voronoi diagrams of semi-algebraic sets Anton, François 2004
pdf
Page Metadata
Item Metadata
Title | Voronoi diagrams of semi-algebraic sets |
Creator |
Anton, François |
Date Issued | 2004 |
Description | Most of the curves and surfaces encountered in geometric modelling are denned as the set of solutions of a system of algebraic equations and inequalities (semialgebraic sets). Many problems from different fields involve proximity queries like finding the (nearest) neighbours or quantifying the neighbourliness of two objects. The Voronoi diagram of a set of sites is a decomposition of space into proximal regions. The proximal region of a site is the locus of points closer to that site than to any other one. Voronoi diagrams allow one to answer proximity queries after locating a query point in the Voronoi zone it belongs to. The dual graph of the Voronoi diagram is called the Delaunay graph. Only approximations by conies can guarantee a proper order of continuity at contact points, which is necessary for guaranteeing the exactness of the Delaunay graph. The theoretical purpose of this thesis is to elucidate the basic algebraic and geometric properties of the offset to an algebraic curve and to reduce the semialgebraic computation of the Delaunay graph to eigenvalues computations. The practical objective of this thesis is the certified computation of the Delaunay graph for low degree semi-algebraic sets embedded in the Euclidean plane. The methodology combines interval analysis and computational algebraic geometry. The central idea of this thesis is that a (one time) symbolic preprocessing may accelerate the certified numerical evaluation of the Delaunay graph conflict locator. The symbolic preprocessing is the computation of the implicit equation of the generalised offset to conies. The reduction of the Delaunay graph conflict locator for conies from a semi-algebraic problem to a linear algebra problem has been possible through the use of the generalised Voronoi vertex (a concept introduced in this thesis). The certified numerical computation of the Delaunay graph has been possible by using an interval analysis based library for solving zero-dimensional systems of equations and inequalities (ALIAS). The certified computation of the Delaunay graph relies on theorems on the uniqueness of a root in given intervals (Kantorovitch, Moore-Krawczyk). For conies, the computations get much faster by considering only the implicit equations of the generalised offsets. |
Extent | 11654252 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-11-27 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0051543 |
URI | http://hdl.handle.net/2429/15860 |
Degree |
Doctor of Philosophy - PhD |
Program |
Computer Science |
Affiliation |
Science, Faculty of Computer Science, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2004-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_2004-901432.pdf [ 11.11MB ]
- Metadata
- JSON: 831-1.0051543.json
- JSON-LD: 831-1.0051543-ld.json
- RDF/XML (Pretty): 831-1.0051543-rdf.xml
- RDF/JSON: 831-1.0051543-rdf.json
- Turtle: 831-1.0051543-turtle.txt
- N-Triples: 831-1.0051543-rdf-ntriples.txt
- Original Record: 831-1.0051543-source.json
- Full Text
- 831-1.0051543-fulltext.txt
- Citation
- 831-1.0051543.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0051543/manifest