L Alberti
email: lalberti@sophia.inria.fr

B Mourrain
email: mourrain@sophia.inria.fr

GALAAD J Wintz
email: jwintz@sophia.inria.fr

Topology and arrangement computation of semi-algebraic planar curves

Keywords: topology, arrangement, implicit, parametric, curves

We describe a new subdivision method to efficiently compute the topology and the arrangement of implicit planar curves. We emphasize that the output topology and arrangement are guaranteed to be correct. Although we focus on the implicit case, the algorithm can also treat parametric or piecewise linear curves without much additional work and no theoretical difficulties.

The method isolates singular points from regular parts and deals with them independently. The topology near singular points is guaranteed through topological degree computation. In either case the topology inside regions is recovered from information on the boundary of a cell of the subdivision.

Obtained regions are segmented to provide an efficient insertion operation while dynamically maintaining an arrangement structure.

We use enveloping techniques of the polynomial represented in the Bernstein basis to achieve both efficiency and certification. It is finally shown on examples that this algorithm is able to handle curves defined by high degree polynomials with large coefficients, to identify regions of interest and use the resulting structure for either efficient rendering of implicit curves, point localization or boolean operation computation.

Introduction

When manipulating a geometric object there are several operations that one wants to be able to carry out, such as visualizing it in an accurate way, have a topologically correct discretization of the object, be able to individuate the connected components of the object (or of its complement), or to test whether a point is in a specific connected component. Furthermore, one wants to be able to dynamically manipulate these geometric objects and have them interact together through operations such as intersection, complement and union.

A very flexible and theoretically robust setting for doing so is semi-algebraic geometry. Our algorithm, that handles all possible boolean operations on semi-algebraic objects in the plane, makes it possible to consider practical examples such as those which come from CAD (Computer Aided Design) or CSG (Constructive Solid Geometry). In order to do so, we study curves and "regions" endowed with the connected component predicate in a rectangular domain

D 0 = [a, b] × [c, d] ⊂ R 2 .
Curves are simply sets defined by one polynomial equality. Regions are the connected components of the complement of a curve.

This article, which corresponds to an expanded version of [START_REF] Alberti | Visualisation of implicit algebraic curves[END_REF] and [START_REF] Wintz | A subdivision arrangement algorithm for semi-algebraic curves: an overview[END_REF], presents a new algorithm to answer the following question: how to interactively and accurately compute the correct topology and regions defined by one or several curves. We first describe the general subdivision framework that we propose to compute these arrangement of curves. Then, we move to the problem of representing region groups and implementing localization, intersection and update operations on them. The solution we present to this problem is a new combinatorial approach that speeds up computations. Finally, we describe the regularity criterion used in the subdivision algorithm and how to make it effective for algebraic curves.

To get a nice visual representation of the curve C we're interested in, we produce a piecewise linear approximation of it (within a given Hausdorff distance ǫ > 0). Besides, we guarantee that this piecewise linear complex is actually topologically equivalent to C. What we guarantee precisely is that there is an injective continuous semi-algebraic deformation of C into its piecewise linear approximation. This is made possible by controlling what happens close to singular points by means of the topological degree which is a quantity that can be computed on the boundary of an isolating box for the singular point. This is a new idea that allows us to control the neighborhood of a singular point only relying on 1-dimensional information around the point. The approach we choose is subdivision. One strong point of this approach regarding rendering is that we can easily control the Hausdorff distance to the actual curve and that we can use the subdivision incrementally to create an interactive visualization environment. In other words, should the user ask to zoom in on any specific part of the curve, we can always refine the approximation by subdividing further to reduce the approximation error to a size smaller than a pixel. This interactivity feature is hardly achievable with physically based methods like the one described in [START_REF] Rösch | Interactive visualization of implicit surfaces with singularities[END_REF]. On the efficiency side, this approach makes as much use of floating number computation as possible which speeds up the computation and enables us to treat curves defined by polynomials with large coefficients and large degree. Such curves appear, for instance, when applying computer algebra techniques on exact representations of geometric objects. We illustrate it on an example of a selfintersection curve of a bicubic parametrized surfaces and on a discriminant curve related to a conjecture on bivariate systems. As for certification, the output correctness is proved under the assumption that we have a certified multivariate solver that can isolate the roots of a zero-dimensional system. Multivariate solving is a difficult problem in itself that can be solved efficiently in most cases by sleeve methods or subdivision methods for example [START_REF] Sherbrooke | Computation of the solutions of nonlinear polynomial systems[END_REF][START_REF] Elber | Geometric constraint solver using multivariate rational spline functions[END_REF][START_REF] Mourrain | Subdivision methods for solving polynomial equations[END_REF], but no matter the approach it seems necessary, in some cases, to resort to purely algebraic techniques such as rational univariate representations of roots [START_REF] Basu | Algorithms in Real Algebraic Geometry[END_REF][START_REF] Elkadi | Introduction à la résolution des systèmes polynomiaux[END_REF]. We will not go into further details about multivariate solving and refer to previous citations and references therein instead.

Arrangements of geometric objects is a field of computational geometry which has been studied for years [START_REF] Agarwal | Arrangements and their applications[END_REF][START_REF] Halperin | Arrangements and their applications in robotics: recent developments[END_REF]. They allow to perform boolean operations which are fundamental tools in Geometric Modeling.

The new method that we describe in this paper for computing arrangements of semi-algebraic curves, provides an efficient way to localize intersection points of region boundaries by storing geo-metric information on the vanishing of the functions which define the algebraic curves in hierarchical structures. Its originality is to prevent useless computation by stopping the subdivision as soon as the topology of the object is known in a cell of subdivision.

Associated topology computation has lead us to a generic subdivision arrangement algorithm where input objects considered in a given domain are subdivided until being regular in a cell of subdivision, building a quadtree of cells. From this quadtree, we easily obtain regions which are organized within an augmented influence graph to describe the arrangement of objects. It is dynamic in the sense that we can maintain this structure while new objects are inserted or existing objects are removed. These insertions and deletions involve to be able to compute respectively intersections and unions of regions.

Our algorithm brings several contributions. First, the nature of a subdivision scheme allows us to focus, if desired, on regions of interest, leading to a multi-resolution arrangement algorithm, parametrized by a level of approximation. Another contribution is that thanks to its generic approach, it can be used to compute an heterogeneous arrangement, i.e. an arrangement of objects having different representations. The resulting regions are also equipped with data structures which lead to efficient operations and are directly usable in a CSG context.

We now give a description of the content of the subsequent subsections: Section 2 places back our method in its scientifical and historical context by explaining what is the current state of the art. Section 3 explains how the generic arrangement computation works and how it is implemented. Section 4 contains the description of the topology algorithm that is used for implicit curves. We isolate the roots of a bivariate polynomial system, using either a Bernstein subdivision solver to approximate efficiently C or algebraic techniques to certify the result. Section 5 shows some experimental results.

Previous work

There are two main types of algorithms that compares to ours concerning topology computation.

The first type is inspired by the Cylindrical Algebraic Decomposition [START_REF] Collins | Quantifier elimination for real closed fields by cylindrical algebraic decomposition[END_REF] algorithm. They use projection techniques based on a conceptual sweeping line perpendicular to some axis that detects the critical topological events, such as tangents to the sweeping planes and singularities. They involve the exact computation of critical points and genericity condition tests and adjacency tests. The approach has been applied successfully to curves in 2D, and even in 3D, 4D [START_REF] Grandine | A new approach to the surface intersection problem[END_REF][START_REF] Keyser | Efficient and exact manipulation of algebraic points and curves[END_REF][START_REF] González-Vega | Efficient topology determination of implicitly defined algebraic plane curves[END_REF][START_REF] Gatellier | Computing the topology of 3dimensional algebraic curves[END_REF][START_REF] Tao | Computing the topology of real algebraic space curves[END_REF] and extended to surfaces [START_REF] Cheng | Determining the topology of real algebraic surfaces[END_REF][START_REF] Mourrain | Isotopic meshing of a real algebraic surface[END_REF].

However, they assume exact input equations and rely on the analysis of the curve at the critical values of its projection. From an algebraic point of view, they involve the computation of (sub)resultants polynomial and of their roots which are algebraic numbers. This can be a bottleneck in many examples with large degree and large coefficients, for which the resultant is difficult to compute, and its real roots even harder to manipulate. Moreover, as these algorithms work by projection, they have to compute every point in the fibers above the points in the projection. In other words, most points that they compute are actually useless for the computation of the final topological description. The complexity of the algorithm can also vary wildly, depending on the direction of projection we choose. And non-degeneracy conditions have to be checked (which can be difficult by itself) to ensure the correctness of the algorithm. The problem is that the choice of projection is not at all related to the geometry of the curve. This is why the Cylindrical Algebraic Decomposition methods are hardly efficient in practice, and are facing complexity problems in higher dimension. They are also intrinsically delicate to apply using approximate computation.

The other type of methods relies on subdivision techniques of the original domain. This process is most commonly used to get approximations of the curve in terms of Hausdorff distance. The most famous family of algorithm using this approach is the marching cube algorithms family [START_REF] Lorensen | Marching cubes: a high resolution 3d surface construction algor ithm[END_REF]. It does not give any guarantee on the topological correctness of its output, but it has inspired some algorithms that do certify that their output has the same topology as the curve (usually in the smooth case). They have already been used for solving several complicated equations. See [START_REF] Sherbrooke | Computation of the solutions of nonlinear polynomial systems[END_REF][START_REF] Elber | Geometric constraint solver using multivariate rational spline functions[END_REF] and the recent improvements proposed in [START_REF] Mourrain | Subdivision methods for solving polynomial equations[END_REF], exploiting preconditioning techniques. Extensions of this approach to higher dimensional objects have also been considered [START_REF] Sederberg | Algorithm for algebraic curve intersection[END_REF][START_REF] Keyser | Efficient and exact manipulation of algebraic points and curves[END_REF][START_REF] Hass | Guaranteed Consistency of Surface Intersections and Trimmed Surfaces Using a Coupled Topology Resolution and Domain Decomposition Scheme[END_REF][START_REF] Joon-Kyung | Contouring 1-and 2-Manifolds in Arbitrary Dimensions[END_REF][START_REF] Liang | Subdivision methods for 2d and 3d implicit curves[END_REF][START_REF] Alberti | Meshing implicit algebraic surfaces: the smooth case[END_REF][START_REF] Plantinga | Isotopic approximation of implicit curves and surfaces[END_REF]. These subdivision methods usually fail when singular points exist in the domain. If a threshold on the minimal size for boxes is not set, the algorithm would run forever. Indeed at singularities, no matter the scale of approximation, the shape and topology of the algebraic objects remain similar.

Our algorithm, like the one in [START_REF] Ratschek | Scci-hybrid methods for 2d curve tracing[END_REF], is hybrid. It combines approximation properties with certification and adaptivity. It subdivides the domain D 0 into regular regions in which the curve is smooth and regions that contain singular points. In the regular regions, we can approximate the curve as precisely as we want and the "singular" regions can be made as small as required. The algorithm computes the topology inside the regions by using what happens on their boundary and we use enveloping techniques to efficiently treat large input equations. This scheme is refined into two concrete algorithms, one being purely numerical and the other one using some algebraic computations. This method combines the advantages of subdivision and Cylindrical Algebraic Decomposition like methods. Its complexity is intrinsic to the geometry of the curve (like the subdivision methods) and it avoids the main drawback of projections methods because it does not need to lift points.

While current methods proposed to compute an arrangement of algebraic curves allow to handle objects of degree 1 such as line segments, objects of degree 2 such as circular arcs and curves of higher degree are still investigated [START_REF] Milenkovic | Calculating approximate curve arrangements using rounded arithmetic[END_REF][START_REF] Wolpert | Jacobi curves: Computing the exact topology of arrangements of non-singular algebraic curves[END_REF][START_REF] Halperin | Arrangements[END_REF][START_REF] Fogel | Arrangements[END_REF][START_REF] Milenkovic | An approximate arrangement algorithm for semi-algebraic curves[END_REF][START_REF] Hanniel | An exact, complete and efficient computation of arrangements of bézier curves[END_REF] and can be used for computing an arrangement of surfaces [START_REF] Mourrain | On the computation of an arrangement of quadrics in 3d[END_REF][START_REF] Mourrain | Isotopic meshing of a real algebraic surface[END_REF]. However, their multiple representations is a major difficulty which leads to as many algorithms as representations exist.

While current methods using a sweep approach [START_REF] Bentley | Algorithms for reporting and counting geometric intersections[END_REF] focus on events, which are critical points for a projection direction, our method using a subdivision approach focuses on regularity criteria and regular regions. When using sweep methods, events are treated when the sweep line encounters points of interest where a projection on a line becomes critical, reducing the dimension of the problem but increasing its computational difficulty (for instance by computing resultants and by lifting points in the case of implicit curves).

On the contrary, a subdivision method avoids the analysis at critical values by enclosing the singular points in a domain from which the problem in hand can be solved. Such methods are less sensitive to numerical approximations of objects and of their intersection points. They have already been used in geometric modelers such as [START_REF] Bowyer | Svlis: Introduction and user manual[END_REF]. Their application to arrangement computation have recently emerged such as in [START_REF] Bowyer | Interval methods in geometric modeling[END_REF][START_REF] Hijazi | Computing arrangements using subdivision and interval arithmetic[END_REF] where interval arithmetic is used to classify cells in the subdivision process. Subdivision methods are also very efficient for isolating the roots of polynomial equations, which appear in geometric problems [START_REF] Sherbrooke | Computation of the solutions of nonlinear polynomials systems[END_REF][START_REF] Elber | Geometric constraint solver using multivariate rational spline functions[END_REF][START_REF] Garloff | Investigation of a subdivision based algorithm for solving systems of polynomial equations[END_REF][START_REF] Mourrain | Subdivision methods for solving polynomial equations[END_REF]. They have also been extended for the approximation of one or two dimensional objects [START_REF] Joon-Kyung | Contouring 1-and 2-Manifolds in Arbitrary Dimensions[END_REF][START_REF] Alberti | Meshing implicit algebraic surfaces: the smooth case[END_REF][START_REF] Liang | Subdivision methods for 2d and 3d implicit curves[END_REF].

Arrangement computation

In this first section we set the basic notations and definitions we use throughout the paper. Then the following sections present the specifics of the method for topology computation, and for manipulating arrangements.

Notations and definitions

The objects that we consider in this framework are semi-algebraic curves whose representation can either be discrete, parametric or implicit. We will denote by O = {o 1 , . . . , o n } the set of input objects of the arrangement algorithm.

For a subset S ⊂ R 2 , we denote by S • its interior, by S its closure, and by ∂S its boundary. We call domain any closed set D such that D • = D and D is simply connected. And we call region any open set R which is a connected component of the complement of an algebraic curve. The boundary of a domain D is denoted ∂D.

We call branch of a curve C (relative to a domain D), any smooth closed connected subset of C ∩ D (i.e. C ∞ diffeomorphic to [0, 1]) and maximal for the inclusion.

We call half branch of a curve C at a point p ∈ D • or half branch originating from p ∈ D • , any smooth closed connected subset of C ∩ D which has one endpoint on ∂D and which is maximal for the inclusion.

A tangent to C is a line, which intersects C with multiplicity ≥ 2. We extend this definition to discrete curves as follows. For a point on one segment of the curve, the line supporting the segment is a tangent line. For a vertex belonging to two segments, a tangent line is any line through this vertex and which defines a half-plane containing the two segments. For a vertex belonging to more segments, any line through this vertex is a tangent line.

A curve C in a domain D is x-critical (resp. y-critical) if it contains a point with a vertical (resp. horizontal) tangent.

A point of C which has a tangent line in any direction is called singular.

The arrangement computation will produce regions which are connected components of the plane R 2 , whose interior does not intersect the objects of O. They might also be called faces of the arrangement. These regions are constructed such that their boundary is a set of edges (a segment of a curved object and two vertices) on objects of O.

In addition to this information, we will associate to a region a bounding box, the set of objects involved in its edges and which determine it. Sign conditions can also be associated to a region with regard of the type of objects which determine it.

The generic arrangement algorithm works incrementally by introducing objects one by one, as follows. For each new object o k , 1. the regions defined by this new object o k independently of the others are computed, 2. these regions are inserted in the arrangement data-structure A k-1 and yields the new data-structure A k .

When a region is inserted in the arrangement data-structure, conflicts are detected and "resolved". We say that two regions conflict together if their intersection is non empty. In what follows, to make the discussion easier, we may say that an object conflicts with a region, which means that a conflict exists between regions determined by this object and another region. We may also say that two objects conflict together, which means that a conflict exists between regions defined by these objects.

In order to be able to dynamically add or remove objects, we use an augmented influence graph I a (see figure 11), which is an influence graph connected together with a conflict graph, that we describe now. An influence graph is a directed, acyclic and connected graph. It possesses a single root, and its nodes correspond to the regions created by an algorithm during its execution. Therefore, a node corresponds to a region defined over the current set of objects at some point during the execution of the algorithm. The influence graph possesses two essential properties: at each step of the algorithm, a region defined over the current set of objects is associated with a leaf of the influence graph, and, the domain of influence of a region associated with a node of the influence graph is contained in the union of the domains of influence of the regions associated with the parents of that node.

In addition to the usual information stored in the influence graph, the augmented influence graph stores a conflict graph between the objects in the current set O and the regions stored in the nodes of the influence graph. This conflict graph is a system of interconnected lists: to each region stored in a node of the influence graph, corresponds a list of objects of O with which it conflicts, and, to each object in the current set O corresponds a list of regions stored in the entire influence graph that conflict with it.

The arrangement is represented by the set of leaves in the augmented influence graph, but, the latter still contains the information required to be able to remove an object from it, or to add an object to it. Indeed, leaf nodes constitute the current arrangement A k , where k is the number of objects in the arrangement, while other nodes (non-root and non-leaf nodes) allow to keep track of the incremental construction of the arrangement A 1 . . . A k-1 . We remind that a root node has no parent, an internal node has both parents and children and a leaf node has no children.

To compute an arrangement A n of a collection O of n objects o 1 . . . o n , we randomly insert the elements of the collection into the arrangement structure A. This structure can be maintained while removing for example the object o k , leading to the arrangement A Σ where Σ = {o 1 . . . o k-1 , o k+1 . . . o n } is the corresponding chronological sequence.

Removing an object from an arrangement can be achieved following the general design explained in [START_REF] Boissonnat | Algorithmic geometry[END_REF]: after the deletion of an object o k , the algorithm reinserts the objects o l of higher chronological rank l > k to create new nodes and re-parent unhooked nodes in the augmented influence graph. We do not detail the removal any further but describe in great details the insertion operation in the next section.

Generic algorithm

In this section, we present the generic insertion operation of the algorithm, i.e. the parts which are not related to the type of input objects. The latter makes the assumption that type related functions will be found in a specialization of it. Inserting an object o k into a structure built for a set O k-1 can be handled in four phases.

Computing regions. In this phase, regions are computed from the topology of o k independently of objects of O k-1 .

Segmenting the boundary. In this phase, computed regions are equipped with additional data structures to help the introduction in the current arrangement.

Locating conflicts. In this phase, newly computed regions are checked for conflict with regions defining the current arrangement.

Updating regions. In this phase, conflicts are dealt with, possibly leading to new regions which are inserted in the data structure.

Computing regions

Regions defined by an object within a given domain D 0 (generally, its bounding box) cannot be computed directly from its representation. To be able to compute a set of regions, from an object in a given domain, we have to ensure that we are in a configuration where we are able to deduce its topology from information on the border of the bounding domain. For the sake of simplicity from here on, we will consider the configurations shown in figure 1, from which we can deduce the topology of regions simply by connecting points of interest of the domain, with points on its border.

When an object is not in one of these configurations, it is subdivided into smaller parts. This process is iterated until one of these configurations is detected. An object in such a configuration is said to be regular. The subdivision is then driven by a regularity test.

There is a strong link between determining if a cell has been subdivided enough and computing the regions defined by the object within the cell. These two type-related functionalities being dependent one from another, they are left for the specialization of an algorithm computing an arrangement of objects of type t. There is however a general scheme to get the regions defined by an object in a given domain.

To each object, we associate a quadtree used to keep track of the subdivision process which allows to deduce regions. The root of the quadtree stores the bounding box of o k . A list of cells is initialized with the root node. While this list is not empty, we check its first element for regularity: if o k is deemed regular in the cell, regions are computed from its topology within the cell and stored in the corresponding node of the quadtree. Else, the cell is subdivided into four children which are appended to the list of cells to be checked for regularity and the quadtree is updated accordingly.

Once all cells have been processed, the leaves of the quadtree contain sub-regions whose union constitutes the regions defined by o k . To compute this union, these regions are merged traversing the tree from its leaves to its root in a process called fusion. The subdivision algorithm, summarized in 3.1 ends up with the root node representing the bounding box of o k containing the regions determined Output: a list of regions. Create an quadtree Q and set its root to B 0 ; Create a list of cells C and initialize it with

[B 0 ]; while C = ∅ do c = pop(C) ; o = objects(O, c) ; if regular(o, c) then Q ← topology(o, c) ; else C ← subdivide(o, c) ; end end return fusion(Q) ;
The following operations remain to be clarified: regularity: the specific operation which checks if regions can be computed from an object within a cell of the subdivision, i.e. if the object is regular in the cell.

subdivide: the generic operation which subdivides a cell into four children, saving computation effort.

topology: the specific operation which computes regions in a regular cell. fusion: the generic operation which merges regions stored in each node of the tree. The regularity test allows to determine if regions can be computed in a cell from information on the boundary of this cell. It is a representation dependent operation and therefore provided in one specialization of this generic algorithm (see section 3.3).

The computation of regions determined by an object in a regular subdivision cell is also provided in a specialization of the generic algorithm, since it strongly depends on the regularity test, which corresponds to different configurations of input objects within a cell (see section 3.3).

When a cell is subdivided into four child cells, we can save computation effort by inheriting the information from a parent to its children , so that we ensure this information is computed only once.

The inheritance is performed in two steps as shown in figure 2: 1. The information contained in the parent cell is inherited by all its children, from a level l to a level l + 1 in the tree (vertical inheritance) 2. Each time the information is computed for one child, it is inherited by its neighbors (horizontal inheritance). Beyond the performance improvement, this inheritance ensures the information shared by different cells to be consistent, especially if solvers used to compute the intersections of an object with a cell are approximate, as it may be the case when using subdivision solvers (see section 3.3). Also, by using a subdivision policy that determines the location of the subdivision depending of the representation of the object, it is possible to significantly decrease the depth of the subdivision.

The fusion of a quadtree consists in merging at each level and across the levels of the tree, the regions stored in the nodes, from its leaves (only nodes where regions have been computed) to the root (where final regions determined by an object will end). The algorithm takes as input a quadtree obtained from the subdivision, in which, only leaves contain regions, and propagates these regions from the leaves to the root. At each step, the set of regions contained in the children of a non-leaf node are merged together.

Segmenting the boundary of a region

Finding all the regions that conflict with a new object is a frequent operation in the arrangement computation so we need to pay special attention to its complexity. Indeed, a naive tree traversal of I a where each region (stored in a node) is checked for intersection with the regions (stored in newly created nodes) to be inserted would lead to a O(p 2 ) complexity, where p is the number of nodes in I a . Also, the respective quadtrees of each object are only used to determine the topology of regions and do not subdivide the same domain. Besides, this structure is related to an object whereas the concept of conflict is relative to a region. To find out the set of regions of I a which conflict with regions created by o k we build a structure called the region segmentation. A region segmentation is a balanced tree data structure in which each node has two children (internal nodes), or none (leaves). Leaf nodes contain the edges defining the boundary of the region as well as their bounding boxes, and internal nodes are associated to the union of the boxes of their children. This way, the root of the tree defining the segmentation is associated to a domain of influence, a subset of the bounding box from which the region has been computed, containing the whole region.

For more efficiency when querying the tree, we propose a simple procedure to build a balanced tree from the list of edges defining a region. As shown in figure 3, to build the segmentation, the list of edges is "flattened", ordered and traversed to build the tree. First, a node is created for each edge and the bounding box of the edge is associated to it, yielding a list of nodes. Second, this list is traversed to build the tree. Since we want the tree to be balanced, the traversal is performed from left to right at step k, then from right to left at step k + 1. Each time a couple of nodes is found in the list, we create a node containing the union of their bounding boxes and reparent the couple of nodes to it, until no more nodes exist in the list. This structure leads to an efficient algorithm to check whether two regions conflict together and helps dealing with them.

A s u b d i v i s i o n c e l l i n a q u a d t

Locating conflicts

To get the list of conflicting zones, we simply compare the segmentations of the regions stored in the nodes of the augmented influence graph with the segmentation of the region currently inserted. First, we compare their domain of influence, i.e. the roots of their respective segmentations. Second, we compare the boxes associated to internal nodes and proceed to child nodes as long as internal nodes boxes do intersect.

If the resulting list is non empty, the intersection of the curve segments are computed using type specific intersection methods according to the type of the objects defining the edges in question. If several intersection points are found, we refine the boundary segmentations in order to have at most one intersection per box.

If the resulting list of pairs of intersecting boxes is empty and if there is no inclusion of a region in the other, there is no conflict and the region can be inserted to I a . The inclusion test simply consists in locating one point of the region to be inserted. If the latter falls inside the other region, and since we already know there is no intersection of edges, we can conclude that the inserted region lies inside the other one. It is then inserted in the graph as a child node of the one representing the other region. To locate a point in a region, we only have to count the number of intersections of a half-line stemming out from the point in an arbitrary direction. The point is inside the region if this number of intersections is odd, the point is outside the region otherwise. The intersection test can be carried out in a reduced complexity using the region segmentation. Indeed, we only test the half-line for intersection with an edge of the region if the latter intersects its bounding box from the segmentation.

f 1 f 2 (a) Computing navigation information f 1 f 2 f 1 ∩ f 2 (b) Walk-about border of regions

Updating regions

When two nodes of the region segmentations intersect, a conflict is detected, yielding a conflict zone (or box). The intersection of the conflicting regions is computed in three phases. First, intersection points between regions are computed in the conflict box. Second, navigation information are computed for each intersection point, as in [START_REF] Zalik | A quick intersection algorithm for arbitrary polygons[END_REF] (see figure 4a). Finally, a walkabout from an intersection point to the next along the edges of the intersecting regions is performed to construct the intersection (figure 4b).

Computing navigation information for each intersection point consists in examining intersections of the edges with the conflict box, to know which edge to follow from an intersection point to the next during the walkabout. To do so, we compute intersection points of each edge of each region (which are oriented counter-clockwise for their outside border and clockwise for their inside border if holes exist) with the bounding box.

When computing the intersection of regions, the walkabout proceeds along the edge whose intersection point with the box is to the left of the other, considering the orientation of edges. This can be easily obtained by sorting intersection points on the boundary of the conflict box, see figure 4a and4b.

The resulting intersection regions are segmented and then inserted in the augmented influence graph as a child node of the regions from which they have been computed.

Specialization for specific curves

When computing an arrangement of curves, we only have to specialize object-type specific functions to meet the generic arrangement algorithm's requirements. We have to provide the following tools:

• Regularity criteria of a curve in a domain.

• Localisation of characteritic points: points with vertical or horizontal tangents, singular points, intersection points, extremal points.

In the following sections, we will detail how to specialize this two points for implicit, paremetric and piecewise curves.

Implicit curves. We denote by f k (x, y) ∈ R[x, y] the polynomial defining the implicit curve corresponding to the object o k . We assume that f k is a squarefree polynomial in Q[x, y]. That is, f k has no square factors in its factorization over CC[x, y], or more intrinsically the ideal

(f k ) generated by f k is radical. The curve o k is the zero set Z(f k ) = {(x, y) ∈ D 0 |f k (x, y) = 0} where D 0 := [a, b]×[c, d] ⊂ R 2 .
is the rectangular domain in which we carry out all of our computations.

The assumption on the f k being squarefree is easily enforceable as it suffices to divide f k by its greatest common divisors with its derivatives along the coordinate axes (or with only one of its greatest common divisor with its derivative along a generic direction). For more details on computing the squarefree part of a polynomial, we refer the reader to the literature on this classical matter [START_REF] Zur Gathen | Modern computer algebra[END_REF]. Of course, once having taken the squarefree part of f k we have not modified Z(f k ) which is the object we are interested in.

The set of singular points of C is

S := Z(f, ∂ x f, ∂ y f ) = {(x, y) ∈ R 2 |f (x, y) = ∂ x f (x, y) = ∂ y f (x, y) = 0}. The set of extremal points of f is Z e (f ) := Z(∂ x f, ∂ y f ) = {(x, y) ∈ R 2 |∂ x f (x, y) = ∂ y f (x, y) = 0}.
We recall that a tangent to the curve C is a line, which intersects C with multiplicity ≥ 2. In particular, any line through a singular point of C is tangent to C.

The specific operations for the arrangement will be performed on the Bernstein representation of f k , starting with the initial domain

D = [a, b] × [c, d]: f k (x, y) = d x,k i=0 d y,k j=0 b k i,j B i d x,k (x; a, b) B j d y,k (y; c, d),
where d x,k is the degree in x of f k and d y,k its degree in y and

B i d (x; u, v) = d i (x -u) i (v -x) d-i (v - u) -d (for 0 ≤ i ≤ d, u < v).
Using the method described in section 4, we obtain the topology of the curve within a regular cell of subdivision, applying the connection algorithm 4.1. The topology of local regions within the regular cell is easily deduced from the topology of the curves lying in the cell together with its corner points. The edges of regions are oriented counterclockwise to ease up following steps of the arrangement algorithm and specify on which side of the border lies the interior of a region.

After the segmentation step, the detection of conflicts reduces either to intersecting regular segments of two different objects or to testing that an endpoint of a regular segment of an object belongs to a given region. To this end, our favorite polynomial solver (see Section 4.4) is used to solve the bivariate polynomial equations f k , f l associate to the objects in potential conflict in the region of interest.

Parametric curves. The computation of an arrangement of parametric curves only requires to be able to compute information on the boundary of a subdivision cell, the subdivision criteria itself, to check if a curve is regular or not within a cell, and a way to compute the intersection of two segments of parametric curves, for detecting region conflicts.

We denote by x(t) and y(t) ∈ R(t) the rational functions representing the parametric curve corresponding to the object o k . For sake of simplicity, we will assume that o k is the image of [0, 1] by the map σ k : t → (x(t), y(t)).

In order to meet the configurations presented in figure 1, we first partition the interval [0, 1] in intervals where x ′ (t) = 0 and y ′ (t) = 0 on the interior. On each one of these intervals, the curve is x and y monotonic. The corresponding bounding box of this curve segment is defined by image by σ k of the end points of the interval. This ensures that the regularity test is realized. Next, we localize which pairs of bounding boxes of two non-consecutive segments intersect, using the segment structure described in Section 3.2.2. If we find two such boxes, which are the images of the intervals I and I ′ , we test if there exists (t, s) ∈ I × I ′ with t = s such that x(t) = x(s) and y(t) = y(s). This reduces to solving a bivariate polynomial system of polynomial equations. Using our polynomial solver (Section 4.4), we isolate the solution and refine the boundary segmentation by inserting isolating boxes around the image of these points by σ k . This allows to test simply regular cells.

To determine the topology of the regions within a cell we have to compute the intersection points of the object it contains with its border. In the case of parametric curves, this can be done easily by solving univariate equations of the form x(t) = x 0 or y(t) = y 0 using solvers as [START_REF] Ioannis | Real algebraic numbers: Complexity analysis and experimentations[END_REF] and by checking that the image by σ k of these roots is on the border of the cell.

Again, we finally turn around these points on the border of the cell in clockwise order and connect them to points of intersect inside the cell, such as self-intersection points or points where a vertical or horizontal tangent exist.

To compute the intersection points of two parametric curves σ k , σ l (with k = l), we solve the bivariate system σ k (t) = σ l (s), with t, s in the intervals ⊂ [0, 1] corresponding to the curve segments. Here, we can use again our polynomial solver (Section 4.4) for this purpose. This allows to detect conflicts between regions which edges are made up of parametric curve segments.

Piecewise-linear curves. When computing an arrangement of piecewise linear curves, the only challenge is to be able to check if the line segments defining a curve are regular in a cell of subdivision. Indeed, computing the intersection of a linear curve segment with a cell border poses no problem since they have the same representation and detecting self-intersection points is not difficult either since it consists in intersecting two line segments.

To check a piecewise linear curve for regularity in a cell, we gather the curve segments lying in the cell having the same derivative sign in the directions x and y. To each curve segment s, we associate a monotony code M (s) given by the tests x i+1 ≥ x i and y i+1 ≥ y i . We traverse the list of curve segments looking for transitions where the monotony code changes, i.e. positions i such that M (s i ) = M (s i+1 ). The transition points are used to define bounding boxes in the the boundary segmentation and then to detect self-intersecting points, following the same approach as for parametric curves.

Regularity and topology of an algebraic curve

The main ingredients of our iterative approach, are to check the regularity criteria of one or several curves in a given domain and to deduce their topology inside this domain. The subdivision algorithm find a partition of D 0 into what we call simple domains D i for which we can compute the topology. Then we can piece together the topologies of the simple domain by gluing them on their boundaries. For each kind of simple domains, we have a connection algorithm that computes a piecewise linear approximation of the curve inside simple domains of that type.

In this section, we describe these operations in details for algebraic curves. The structure of this exposition is as follows: We first define the types of simple domains we will use in subsection 4.1. Then we devote subsections 4.2 and 4.3 to each type of simple domains to explain how to test whether a given domain is indeed of a given type and how to compute the topology in the domain once its type is known. Finally in subsection 4.4, we explain how to actually obtain a partition of the initial Figure 5: From left to right: A simply singular domain, a domain that fails to be simply singular because there are too many points on its boundary, another one that fails to be simply singular because there are several extremal points of f inside it, an x-regular domain, and a domain that fails to be x-regular because there are vertical tangents to C inside it.

bounding box into simple domains with a subdivision algorithm that makes use of the tests introduced before, together with known root isolation techniques which are shortly presented in the same section.

Regularity criteria

We distinguish three different types of simple domains: x-regular domains, y-regular domains and simple singular domains. Definition 4.1 A domain D is x-regular (resp. y-regular) for C if C is smooth in D and it has no vertical (resp. horizontal) tangents. This is algebraically formulated as the following condition:

Z(f, ∂ y f ) ∩ D = ∅ (resp. Z(f, ∂ x f ) ∩ D = ∅).
We might equivalently say that the curve C is x-regular (resp. y-regular) in D instead of saying that D is x-regular (resp. y-regular) in D. Remark 4.2 Pay attention to the fact that x-regularity is a condition on the partial derivative along y. It ensures that the orthogonal projection on the x-axis is a submersion. The same remark applies to y regularity.

Finally we say for short that a curve is regular in D, or equivalently that D is regular for C if C is x-regular or y-regular in D. One may at this point be slightly confused by the fact that we speak freely of number of points on the boundary although it could happen that the curve runs along the boundary and there are thus infinitely many points on the boundary. As a matter of fact, the domains are rectangles, and if this happened this would mean that either (yα) or (xα) could be factored out of the defining equation of C (for some α ∈ Q). It makes it very easy to test for this condition and it can be remediated in several ways: effectively factoring out the linear factor and considering it as a new object that can be added back to the curve afterwards by virtue of its simplicity, or avoid creating domains whose boundary are along a line of the curve by tweaking the final subdvision algorithm. For the sake of clarity we will assume the intersection of C with every domain considered is 0-dimensional in the rest of the discussion as we feel it would pointlessly burden the exposition.

Finally, these definitions will allow us to prove that the topology in regular domains is that of a number of lines stacked up as a function graph, and to prove that the topology inside simply singular domains is that of a cone over the trace of the curve on the domain's boundary. This will be explained in the next two subsections.

Regular domains

In this section, we consider a curve C in R 2 , defined by the equation

f (x, y) = 0 with f ∈ Q[x, y] and a domain D = [a, b] × [c, d] ⊂ R 2 .
We are going to show that if C is x-regular in D, then its topology can be deduced from its intersection with the boundary ∂D. By symmetry the same applies when C is y-regular. We only require that ∂D ∩ C be 0-dimensional. This is a again a mild requirement that can be easily taken care of when choosing a partition of the initial domain or factoring out the corresponding linear component out of C. And we necessarily fall into one of these cases because ∂D ∩ C is 0-dimensional.

These conditions can be effectively tested using the sign s y of ∂ y f , the order k of the first x derivative of f that does not vanish, and the sign s x of ∂ k f . k is well defined because if all this partial derivatives were 0, the whole horizontal line would be included in C which would mean C ∩ ∂D is not 0-dimensional. In the following the points with double index (+-or -+) are considered as double points, one with "smaller x component" than the other (although they correspond to a single point that has only one x component). The one with smaller x component gets the left part of the double index, and the one to its right (bigger x component) gets the right part.

Lemma 4.8 If C is x-regular in D, then a branch of C ∩ D connects a point p of x-index + to a point q of x-index -, such that x p < x q .
Proof. As the curve is x-regular, it has no vertical tangent and thus no closed loop in D. Consequently, each of the interior connected components of C ∩ D intersects ∂D in two distinct points p, q ∈ C ∩ ∂D. Without loss of generality we can assume that x p ≤ x q .

Let s : t ∈ [0, 1] → s t = (x(t), y(t)) be an analytic parameterization of the branch [p, q] with s 0 = (x(0), y(0)) = p, s 1 = (x(1), y(1)) = q. Assume that the x-index of q is +, then by definition of the x-index, there is a point q ′ = s(T ) on the branch [p, q] such that x q ′ > x q . This means x(T ) > x(1)

Location of p on ∂D

Parity of k

Sign determinant Resulting x-index p ∈](a, c), (a, d)[

I p = + p ∈](b, c), (b, d)[ I p = - p ∈](a, c), (b, c)[ k odd s y s x > 0 I p = + p ∈](a, c), (b, c)[ k odd s y s x < 0 I p = - p ∈](a, c), (b, c)[ k even s y s x > 0 I p = +- p ∈](a, c), (b, c)[ k even s y s x < 0 I p = -+ p ∈](a, d), (b, d)[ k odd s y s x > 0 I p = - p ∈](a, d), (b, d)[ k odd s y s x < 0 I p = + p ∈](a, d), (b, d)[ k even s y s x > 0 I p = -+ p ∈](a, d), (b, d)[ k even s y s x < 0 I p = +- p = (a, c) s y s x (-1) k < 0 I p = + p = (a, c) s y s x (-1) k > 0 I p = +- p = (a, d) s y s x (-1) k > 0 I p = + p = (a, d) s y s x (-1) k < 0 I p = +- p = (b, c) s y s x (-1) k < 0 I p = - p = (b, c) s y s x (-1) k > 0 I p = +- p = (b, d) s y s x (-1) k > 0 I p = - p = (b, d) s y s x (-1) k < 0 I p = +- Table 4.2. Definition of the x-index in a box D = [a, b] × [c, d].
and thus ∂ t x(t) is negative on the whole segment [T, 1] as it does not vanish by x-regularity. Similarly for p, as x p ≤ x q we also have x(T ) > x(0) and thus ∂ t x(t) is positive on the whole of [0, T ]. Therefore ∂ t x(T ) is both negative and positive as pictured on figure 7a, this is a contradiction and thus the x-index of q is -. The symmetrical reasoning shows that the x-index of p is +. Finally, since ∂ t x(t) > 0 for s ∈ [0, 1], we have x p < x q , which proves the lemma. Lemma 4.9 Suppose that C is x-regular in D and let p, q be two consecutive points of C ∩ ∂D with: q such that x q is minimal among the points with x-index= -, and x p < x q , then p, q belong to the same branch of C ∩ D.

Proof. Suppose that p and q are not on the same branch. Let p ′ be the other endpoint of the branch going to q. Let q ′ be the other endpoint of the branch starting from p. By lemma 4.8, x-index(p ′ ) = + and x p ′ < x q . By that same lemma, x-index(q ′ ) =and x p < x q ′ . The branch (p ′ , q) separates D in two connected components. We call C r the one whose boundary B r = ∂C r contains the point p (see figure 7b).

Because (p ′ , q) and (p, q ′ ) do not intersect, p and q ′ are in the same connected component of D -(p ′ , q) and on B r .

Consider the sub-boundary {x ≥ x q } ∩ B r . It must be connected. Otherwise the branch (p ′ , q) would intersect x = x q in two distinct points and the curve would have a x-critical point in between. We denote by q, q, the endpoints of {x ≥ x q } ∩ B r (with possibily q = q). We decompose B r as the union of arcs B r = (p ′ , q) ∪ (q, q) ∪ (q, p ′ ) with (q, q) ⊂ ∂D, (q, p ′ ) ⊂ ∂D. By minimality of x q , we have x q ′ ≥ x q so that q ′ ∈ {x ≥ x q } ∩ B r = (q, q). Because x p < x q and p ∈ B r and p ∈ (p ′ , q), we have p ∈ (q, p ′ ) ⊂ ∂D.

This proves that p is in-between p ′ and q ′ and q ′ in-between p and q on ∂D. Therefore, p and q cannot be consecutive points of C on ∂D. By way of contradiction, we conclude that p and q must be on the same branch of C. Proof. We prove the proposition by induction on the number N (C) of points on C ∩ ∂D. We denote this set of points by L. Since the curve has no vertical tangent in D and has no closed loop, each of the connected components of C ∩ D • have exactly two distinct endpoints on ∂D. Thus if N (C) = 0, then there is no branch of C in D. Assume now that N (C) > 0, and let us find two consecutive points p, q of L with x-index(p) = +, x-index(q) = -, x p < x q and x q minimal. By lemma 4.9, the points p, q are the endpoints of the branch of C. Removing this branch from C, we obtain a new curve C ′ which is still x-regular and such that N (C ′ ) < N (C). We conclude by induction hypothesis, that the topology of C ′ and thus of C is uniquely determined. Proposition 4.11 If C has at most one x-critical or y-critical point ∈ D, which is also smooth, then its topology in D is uniquely determined by its intersection with the boundary of D.

Proof. Suppose C has at most one x-critical point in D, which is smooth, then the curve is smooth in D and has no closed loop inside D (otherwise the number of x-critical points would be at least 2). Therefore, the branches are intersecting ∂D in two points. If there is no x-critical point on a branch, by Lemma 4.8 their x-index ∈ {-, +} are distinct. If the branch has a x-critical point of even multiplicity, then the x-index of the end-points of the branch in C are the same. If there are only two points of C on ∂D, then this branch is connecting the two points. As the curve is smooth, the branches are not intersecting. If there are more points, and thus more than 2 branches, the branch with the even x-critical point is separating the set of branches into two disjoints subsets of branches with no x-critical points. Changing the orientation of the x-axis if necessary, we can find consecutive points p, q on ∂D which satisfies the hypothesis of lemma 4.9. By this lemma, they are necessarily on the same branch of one of these two subsets. Removing this branch from C and processing recursively in this way, we end up either with no point on ∂D or two points on ∂D with the same x-index. These points are necessarily connected by the branch containing the x-critical point of C in D.

This leads to the following algorithm: Take a point q such that x q is minimal among the points in L with x-index= -; Take the point p that follows or precedes q in L such that x p < x q (thus x-index(p) = +) ; Add the arc [p, q] to the set B of branches and remove p, q from L ; end Notice that a sufficient condition for the x (resp. y) regularity of f in a domain D is that the coefficients of ∂ y (resp. ∂ x f ) in the Bernstein basis on D are all > 0 or < 0. In this case the connection algorithm can be simplified even further. See [START_REF] Alberti | Meshing implicit algebraic surfaces: the smooth case[END_REF] for more details.

Simply singular domains

In this section we deal with simple singular domains. Recall that those domains D are such that they contain a unique critical point p of f where f also vanishes (i.e. p is a singular point of C) and the number of half branches at p is the same as the number of points of C on their boundary. As mentioned earlier we will see in section 4.4, how to compute such a domain. This section will focus on explaining how to test whether a given domain is simply singular, as well as how to compute its topology.

The next two subsections explain how it is possible to effectively count the number of half branches of C at p. The first one introduces the essential theoretical tool to do this which is the topological degree, [START_REF] Lloyd | Degree theory[END_REF]. Then the subsection that follows deals with Khimshiashvili theorem, [START_REF] Khimšiašvili | The local degree of a smooth mapping[END_REF][START_REF] Aoki | On the number of branches of a plane curve germ[END_REF][START_REF] Szafraniec | On the number of branches of a 1-dimensional semianalytic set[END_REF] which actually relates the number of half branches and the topological degree. This enables us to check effectively that the number of points in ∂D ∩ C is the same as the number of half branches at p.

The third and last subsection for simply singular domains show that the topology inside a simply singular domain is conic. This fact gives rise a very simple connection algorithm for them.

Topological Degree

In this section, we recall the definition of the topological degree in two dimensions and how it can be computed. See [START_REF] Lloyd | Degree theory[END_REF][START_REF] Stenger | Computing the topological degree of a mapping in R n[END_REF] for more details.

Let D be a bounded open domain of R 2 and F = (f 1 , f 2 ) : D → R 2 a bivariate function which is two times continuously differentiable in D. A point p ∈ R 2 is said to be a regular value of F on D if the roots of the equation F (x, y) = p in D are simple roots, i.e. the determinant of the Jacobian J F of F at these roots is nonzero). It can be proved that this construction does not depend on the regular value q in the same connected component of R 2 -F (∂D) as p [START_REF] Lloyd | Degree theory[END_REF]. If p is a regular value of F on D, we can take q = p. Remark 4. [START_REF] Dickenstein | Extremal real algebraic geometry and a-discriminants[END_REF] The topological degree has a geometric interpretation known as the degree of the "Gauss map". It is the number of times F (p) goes around F (D) when p goes around D one time. And it is negative when F reverses the orientation of D. The red arrows in fig. 8 picture the F (p) on the boundary. This viewpoint allows to use the strong geometric intuition behind the gradient field when F is the gradient map of f .

Let us now give a more explicit formula for computing this topological degree, which involves only information on the boundary of D.

Proposition 4.14 [START_REF] Stenger | Computing the topological degree of a mapping in R n[END_REF] Assume here that the boundary D is a polygon and that it is decomposed in reverse clock-wise order into the union of segments

∂D = ∪ g i=1 [p i , p i+1 ], p g+1 = p 1 ,
in such a way that one of the component

f σi (σ i ∈ {1, 2}) of F = (f 1 , f 2 ) has a constant sign ( = 0) on [p i , p i+1 ]. Then deg[F, D, (0, 0)] = 1 8 g i=1 (-1) σi-1 sg(f σi (p i )) sg(f σi (p i+1 )) sg(f σi+1 (p i )) sg(f σi+1 (p i+1 )) (1)
where f 1 = f 3 and sg(x) denotes the sign of x. Thus in order to compute the topological degree of F on a domain D bounded by a polygon, we need to separate the roots of f 1 from the roots of f 2 on ∂D by points p 1 , . . . , p g+1 at which we compute the sign of f 1 and f 2 . This will be performed on each segment of the boundary of D, by a univariate root isolation method working simultaneously on f 1 and f 2 , that we will described in the next section.

Figure 8 shows a sequence of points p 1 , . . . , p 9 , which decomposes ∂D into segments on which one of the two functions (f 1 = 0 and f 2 = 0 are represented by the plain and dash curves) has a constant sign. Computing the sign of these functions and applying formula (1) yields the topological degree of F = (f 1 , f 2 ) on D at (0, 0).

Counting the number of branches

Let us consider a curve C in a domain D ⊂ R 2 , defined by the equation f (x, y) = 0 with f (x, y) ∈ R[x, y]. Let ∇f = (∂ x f, ∂ y f ) be the gradient of f . A point p ∈ C is singular if ∇f (p) = 0. We defined a real half branches of C at p, as a connected component of C -{p} ∩ D(p, ǫ) for ǫ > 0 small enough.

The topological degree of ∇f can be used to count the number of half branches at a singular point, based on the following theorem: Theorem 4.15 (Khimshiashvili [START_REF] Khimšiašvili | The local degree of a smooth mapping[END_REF][START_REF] Aoki | On the number of branches of a plane curve germ[END_REF][START_REF] Szafraniec | On the number of branches of a 1-dimensional semianalytic set[END_REF]) Suppose that p is the only root of ∇f = 0 in D. Then the number N of real half branches at p of the curve defined by f (x, y) = f (p) is

N = 2 (1 -deg[∇f, D, (0, 0)]).
(2)

We will denote by N (f, D) the number given by Formula (2).

In order to count the number of branches of C at a singular point p ∈ C, first we isolate the singular point p in a domain D, so that ∇f does not vanishes elsewhere in D. Then we compute the topological degree deg[∇f, D, (0, 0)], as described previously, by isolating the roots of ∂ x f and ∂ y f on ∂D.

The subdivision procedure

Let D 0 = [a, b] × [c, d] be a domain of R 2 .
The goal of this section is to describe effective methods to partition D 0 into simple domains. The difficult step of this approach is to isolate the roots of

Z e (f ) = {(x, y) ∈ D 0 , ∂ x f (x, y) = 0, ∂ y (f )(x, y) = 0}.
which are on C, with the following property:

• There is only one point p of Z e (f ) in each isolating domain D (and it is on C, that is singular)

• The number of points in C ∩ ∂D 0 is the number of half-branches at the singular point p (that is

N (f, D) = 2 (1 -deg[∇f, D, 0])).
We present two approaches. The first one exploits the Bernstein representation of f and subdivision techniques to isolate the roots of Z e (f ), while identifying domains where the curve is regular. It outputs an approximation of C to a precision that is given as input to the algorithm. We prove that, for a sufficiently high precision, the algorithm output has the same topology as C. The second algorithm is based on algebraic techniques (namely Rational Univariate Representation) and is guaranteed to output the correct topology.

The two following methods do the isolation work in a different way but they share the test described in section 4.3 to count the number of half branches at a singular point.

Recursive Bernstein subdivision method

We describe here the subdivision method used to obtain such isolating domains, with a specialization of the approach used in [START_REF] Mourrain | Subdivision methods for solving polynomial equations[END_REF]. See also [START_REF] Sherbrooke | Computation of the solutions of nonlinear polynomial systems[END_REF][START_REF] Elber | Geometric constraint solver using multivariate rational spline functions[END_REF]. This method which we recall for polynomials in Q[x, y] applies for general multivariate polynomials. We are going to consider the system f (x, y) = 0,

∂ x f (x, y) = 0, ∂ y f (x, y) = 0 in the domain D 0 = [a, b] × [c, d].
Each of these polynomials is expressed in the Bernstein basis on where h ∈ {f, ∂ x f, ∂ y f } and d x is the degree of h in x, d y the degree of h in y. By using a method described in [START_REF] Mourrain | Subdivision methods for solving polynomial equations[END_REF] we can quickly generate a set of boxes where the curve is x or y-regular and small set of boxes of size smaller than a given precision ǫ > 0 that isolates the part of the curve where we don't yet know what is happening. The principle of this method is to either reduce a box by using convexity inequalities on Bernstein bases or to split the boxes if the inequalities do not apply. This is the main loop of the subdivision algorithm, which is combined with preconditioning techniques to improve the performance of the solver. The computation is iterated until the size of the box is smaller than ǫ.

It is not the purpose of this article to go into the details of the method expounded in [START_REF] Mourrain | Subdivision methods for solving polynomial equations[END_REF]. In this method the domain are boxes [a, b] × [c, d] and when they are reduced, they are reduced in only one. That is, either the [a, b] or the [c, d] part of the box is shrinked. When that happens, the convexity inequalities imply that one of the functions f , ∂ x f , ∂ y f does not vanish in the regions which are removed from the box. Thus the curve C in these regions is regular and according to section 4.2, its topology can be deduced from the intersection of the curve with the boundary of the region.

This method can be adapted to our implicit curve problem, and yields the following algorithm:

This algorithm decomposes the initial domain into regions where the topology is known and a set of non-intersecting boxes of size ≤ ǫ where ♯(∂D ∩ C) = 2 (1deg[∇f, D, 0]) (this is (2)). If ǫ corresponds to the size of a pixel, the visualization of the curve will be correct, except in these pixel boxes, which we call singular regions. Inside them equation ( 2) holds, and if in addition there is a unique critical point of f , which is also on C, then the computed topology is correct.

During the subdivision process we have to zoom on domains or equivalently to scale the variables (x := λ x, y := λ y). In order to handle the numerical instability problems, which may happen in this scaling step or when we have to deal polynomials with large coefficients and degrees, we use the following enveloping techniques, which allows us to compute with fixed precision numbers: To analyze the curve C defined by the polynomial f ∈ Q[x, y] on a domain D = I × J,

• we convert f to the Bernstein basis on the domain D using exact arithmetic: f (x, y) = i,j γ i,j B i dx (x; I) B j dy (y; J) • we round up and down to the nearest machine precision number γ i,j ≤ γ i,j ≤ γ i,j , so that we have f (x, y) ≤ f (x, y) ≤ f (x, y) on D.

• We use the interval coefficients [γ i,j , γ i,j ] to test the sign conditions and to remove the regular regions.

It can be proved that if ǫ is small enough, then this algorithm compute the topology of C (but for space limitation reasons, we do not include the proof here).

Remark that if Z(f ) is smooth in a domain D, this algorithm can be run with ǫ = 0 and will terminate (and output the correct topology) as every subdomain will ultimately be x-regular or yregular.

One could think that "scaling" f to have its coefficient lie inside a "good" range could be useful to enhance rounding. This is to a very large extent a fallacy as IEEE standard machine floating numbers are encoded in the format 0.m10 e where m (the significand) and e (the exponent) are integers with a fixed number of digits in base 2 (23 for m and 8 for e). The operation of scaling would therefore affect primarily e without improving in any way the rounding of the rational number which primarily takes place in m. Scaling would prove useful in situations where the amplitude of the exponents of the coefficients of the polynomial is less than 2 8 (= 256) and the lower exponents are below -128. This situation hardly ever arises in practice where the exponents typically lay between -30 and 30 even in the most degenerate cases encountered after a great deal of subdivisions. And this is by no means close to 10 -128 .

Subdivision based on rational univariate representation

Choosing the precision parameter ǫ smaller than some bound was enough to certify the output of the previous algorithm. The drawback is that the bounds are difficult to compute and are bad because uniform. The algebraic technique we present hereafter, namely RUR (rational univariate representation), is guaranteed to yield the correct topology. It allows the algorithm to use coarser approximations of roots (when the critical points of f are far away from each others).

We explain in short what RUR are in the bivariate case (see [START_REF] Basu | Algorithms in Real Algebraic Geometry[END_REF] for more details). When given a system of equations E = {f 1 = 0, f 2 = 0} in R2 with 0-dimensional solution space, it is possible to find polynomials P, P 1 , P 2 ∈ R[u] so that we have Z(E) = P1 P ′ (α), P2 P ′ (α) | α ∈ R, P (α) = 0 where P is squarefree and P ′ is its derivative. In other words, the roots of E are the image of the roots of P by a rational map. A RUR of the roots of E can be computed by finding a separating linear function and using resultant or Groebner basis techniques.

In our case the following problem arises: Z e (f ) can have 1-dimensional components. Because we are dealing with curves in R 2 , we can easily separate the 1-dimensional part from the 0-dimensional part by computing g := gcd(∂ x f, ∂ y f ). We define

Z 1 e (f ) = Z(g), Z 0 e (f ) = Z ∂ x f g , ∂ y f g .
Among the points in Z 0 e (f ) we want to be able to tell those that are in C, that is those which are singular points of Z(f ). This way we can isolate the singular points of C from the rest of N e (f ). Since f is square-free, the singular locus S of f is 0-dimensional and Z 0 e (f ) ∩ C = Z1 e (f ) ∩ C. Therefore we compute (P, F 1 , F 2 ) a RUR for Z 0 e (f ) instead of Z e (f ) to isolate the critical points of f . And to tell which points are on C we compute Q = gcd(P, num f (F 1 , F 2 )) where num takes the numerator an irreducible rational fraction. It can be check easily that (Q, F 1 , F 2 ) is a RUR for Z 0 e (f ) ∩ C by using the fact that P ′ and P have no common roots. Now, we use this RUR to isolate the roots of squarefree polynomial P using a univariate solver (see e.g.. [START_REF] Ioannis | Real algebraic numbers: Complexity analysis and experimentations[END_REF]). By using interval arithmetic one can find isolating intervals for the roots of Z 0 e (f ) by computing the images of the isolating intervals of the roots of P by F 1 := P1 P ′ and F 2 := P2 P ′ . This generates boxes containing these roots. If the boxes intersect we refine the isolating intervals of the roots of P until the boxes do not intersect anymore. Finally, using again interval arithmetic, we check that g does not vanish in these isolating boxes. Otherwise we refine them until it doesn't.

Keeping the boxes which correspond to roots of Q, we obtain isolating boxes which contain a single singular point. For each isolating box D, we compute the topological degree. If N (f, D) is not the number of points of C ∩ ∂D, we refine the isolating box.

This yields isolating boxes for the singular points of C, which are simply singular. The complement of the isolating boxes is divided into boxes on which we apply the previous subdivision algorithm for smooth curves.

Examples

Algorithms presented in this paper have been completely implemented with the algebraic geometric modeling software Axel 1 . The efficiency of the topology algorithm presented here allows a real time manipulation of algebraic objects within the software, whereas current solutions usually only propose ray tracing algorithm for visualization. Here are some significant illustrations. More examples 2 and videos3 can be found on the software's website.

The curve in figure 9 is the preimage in the parameter space of a self-intersection point of a bicubic surface. Its equation has been obtained by resultant computation. It is of total degree 76 and of degree 44 in each parameter. Its coefficients are of maximal bit size 590. It takes 7s to visualize this curve.

Figure 10 shows the discriminant curve of a bivariate system with few monomials that gives a counter-example to Kushnirenko's conjecture [START_REF] Dickenstein | Extremal real algebraic geometry and a-discriminants[END_REF]. It is of degree 47 in x and y, and the maximal bit size of its coefficient is of order 300. It takes less that 10 seconds to visualize it within the modeler. The area that is emphasized looks like a cusp point, but when we blow it up in the thumbnail shown in bottom right corner of figure 10, we see that it is actually made of 3 cusp points and 3 crossings. The counter-example comes from this area.

The dynamic part of the arrangement algorithm has lead to a daemon in the modeler which looks for insertions and removal of objects to maintain its data structure. The generic part has lead to an abstract algorithm implemented following the template method and visitor design patterns [START_REF] Gamma | Design patterns: elements of reusable objectoriented software[END_REF]. We have put a special emphasis on keeping the structure accessible to the user. Using the model-viewcontroller pattern, data structures such as the augmented influence graph and the various quadtrees can be interactively displayed and queried for point location (see figure 11).

Figure 11 illustrates how the algorithm behaves with a large amount of implicit curves of degree up to 9 in degenerate cases, that is to say with tangential intersection points and coincident singular points. In this example, the subdivision process takes no longer than 7.224 seconds, the corresponding quadtree depth is 10.

Summary and outlook

Topology computation methods are usually compared in terms of exactness and efficiency. The method we propose certifies the topology of singular curves even for very high degrees. Moreover, using the subdivision scheme together with the Bernstein representation of input object coupled with enveloping techniques, we achieve very fast computations. To the best of our knowledge, our method is more than 20 times faster than best existing other implementations.

Arrangement computation is a new research field in computational geometry in which sweep methods are systematically employed, as a natural evolution of vertical decomposition algorithms. Such methods are very sensitive to numerical errors whereas ours benefits from the topology algorithm efficiency and correctness. Keeping track of the subdivision process also permits faster queries on the arrangement result.

Both the topology and the arrangement algorithm remain to be extended to treat surfaces. The technique has already proved reliable for computing the topology of implicit space curves [START_REF] Liang | Subdivision methods for 2d and 3d implicit curves[END_REF] and the subdivision scheme naturally extends to higher dimensions.

Figure 1 :

 1 Figure 1: Allowed configurations within a cell.

Figure 2 :

 2 Figure 2: Inheritance.

Figure 3 :

 3 Figure 3: Building the region segmentation.

Figure 4 :

 4 Figure 4: A generic boolean operation

Definition 4 . 3

 43 A domain D is simply singular for C if S ∩ D = {p} and if the number n of half branches of C at the singular point p is equal to ♯(∂D ∩ C), the number of points of C on the boundary of D.

Remark 4 . 4 Definition 4 . 5 Figure 6 :

 44456 Figure 6: x-indices of an x-regular domain

Table 4 .

 4 2 summarizes how to obtain the x-index I p of a point p ∈ D from these 3 numbers for a box D = [a, b] × [c, d].

  (a) p and q cannot have same x-index.(b) Illustrative figure for proof of lemma 4.9

Figure 7 :

 7 Figure 7: Connection algorithm.

Proposition 4 .

 4 10 Let C = Z(f ). If D is a x-regular domain, the topology of C in D is uniquely determined by its intersection C ∩ ∂D with the boundary of D.

Algorithm 4 . 1 :

 41 Connection for a x-regular domain Input: an algebraic curve C and a domain D = [a, b] × [c, d] ⊂ R 2 such that C has no vertical tangent in D Output: the set B of branches of C in D Isolate the points C ∩ ∂D and compute their x-index ; Order the points of C ∩ ∂D with non-zero x-indices clockwise and store them in the circular list L ; while L = ∅ do

Definition 4 . 12

 412 Let p ∈ R 2 and suppose further that the roots of the equation F (x, y) = p, are not located on the boundary ∂D.Then the topological degree of F at p relative to D, denoted by deg[F, D, p], is defined bydeg[F, D, p] = x∈D:F (x)=q sign det J F (x),for q a regular value of F on D in the connected component of R 2 -F (∂D) containing p.

Figure 8 :

 8 Figure 8: Computing the topological degree

  j B i dx (x; a, b) B j dy (y; c, d),

Figure 9 :

 9 Figure 9: Topology of a high degree curve in Axel.

Figure 10 :

 10 Figure 10: Topology of a curve with hidden cusp points.

Figure 11 :

 11 Figure 11: Computing an arrangement in Axel.

http://axel.inria.fr

http://axel.inria.fr/user/screenshots

http://axel.inria.fr/user/screencasts

Let us describe now the algorithm used to compute the topological degree of ∇f in a domain D = [a, b] × [c, d]. According to formula [START_REF] Agarwal | Arrangements and their applications[END_REF], this reduces to separating the roots of ∂ x f ∂ y f on the boundary of D, which consists in 2 horizontal and vertical segments. The problem can thus be transformed into isolating the roots of univariate polynomials on a given interval. Hereafter, these polynomials will be called g 1 (t), g 2 (t) and the interval [u, v] ⊂ R. For instance, one of the 4 cases to consider will be g 1 (t) = ∂ x f (t, c), g 2 (t) = ∂ y f (t, c), u = a, v = b. We recall briefly the subdivision method described in [START_REF] Mourrain | On the complexity of isolating real roots and computing with certainty the topological degree[END_REF][START_REF] Mourrain | Bernstein's basis and real root isolation[END_REF][START_REF] Ioannis | Real algebraic numbers: Complexity analysis and experimentations[END_REF], which can be used for this purpose. First we express our polynomials

where

The number of sign variations of the sequence

. By a variant of Descartes rule [START_REF] Basu | Algorithms in Real Algebraic Geometry[END_REF], it bounds the number of roots of g k on the interval [u, v] and is equal modulo 2 to it. Thus if

. This is the main ingredient of the subdivision algorithm [START_REF] Ioannis | Real algebraic numbers: Complexity analysis and experimentations[END_REF], which splits the interval using de Casteljau algorithm [START_REF] Farin | Curves and surfaces for computer aided geometric design : a practical guide[END_REF]

and remove it otherwise. It iterates the process on each subintervals until the number of sign variation is 0 or 1. The complexity analysis of the algorithm is described in [START_REF] Ioannis | Real algebraic numbers: Complexity analysis and experimentations[END_REF]. See also [START_REF] Eigenwillig | Almost tight recursion tree bounds for the descartes method[END_REF].

In our case, we need to compute intervals on which one of the polynomial g 1 or g 2 has a constant sign. Thus we replace the subdivision test by the following:

• otherwise we split it and compute the Bernstein representation of g k (k = 1, 2) on the two subintervals using de Casteljau algorithm and repeat the process. This yields the following algorithm for computing the topological degree of ∇f = (f 1 (x, y), f 2 (x, y)) on D:

If we assume that ∂ x f and ∂ y f have no common root on the boundary of D, it can be proved (by the same arguments as those used in [START_REF] Basu | Algorithms in Real Algebraic Geometry[END_REF][START_REF] Mourrain | Bernstein's basis and real root isolation[END_REF][START_REF] Ioannis | Real algebraic numbers: Complexity analysis and experimentations[END_REF]) that this algorithm terminate and output a sequence of intervals on which one of the functions g 1 , g 2 has no sign variation. The complexity analysis of this method is described in [START_REF] Mourrain | On the complexity of isolating real roots and computing with certainty the topological degree[END_REF]. This analysis can be improved by exploiting the recent results in [START_REF] Ioannis | Real algebraic numbers: Complexity analysis and experimentations[END_REF].

Conic structure and connection algorithm

Finally we prove that the topology in a simple singular domain D is conic and write a connection algorithm for theses domains.

Let A ⊂ R n and p ∈ R n . We call cone over A with center p the set p ⋆ A := q∈A [p, q]. Proof. s is the unique critical point of f in D. If the endpoint of a half branch at s is not on ∂D, the half branch has to be a closed loop inside D. In that case, f would be extremal at some point p ( = s) Compute the restriction g 1 (t) (resp. g 2 (t)) of f 1 (resp. f 2 ) on this side segment I and its representation in the corresponding Bernstein basis ; L := {I}; while L = ∅ do pop up an interval [p, q] from L ; if V (g 1 ; p, q) = 0 or V (g 2 ; p, q) = 0 then insert p, q clockwise in the circular list B ; else split [p, q] in half an insert the two subintervals in L; end end end Compute N given by formula (1) for the points in the circular list B ; inside the loop, and p would be another critical point of f inside D. Thus, by way of contradiction, the endpoints of half branches at s have to be on ∂D.

The number of half branches at s is exactly ♯(∂D ∩ C). As no two half branches can have the same endpoint on ∂D (that would be another singular point in D), all points on ∂D are endpoints of half branches at s. Thus, at this point, we know that the connected component of s inside D is conic.

But in fact, there is no other connected component: Suppose we have another connected component α of C intersecting D. As all points of ∂D ∩ C are connected to s, we have α ⊂ D. α is a smooth 1-dimensional manifold because s is the only singular point. Therefore α is a closed loop inside D (s might be inside it). We look at the complement of C in R 2 , it has a bounded connected component because one of them is inside the loop α. As f vanishes on the boundary of this component, f has an extremum inside it. This extremum cannot be s as it is in the complement of f , which is impossible. Thus, C ∩ D is connected. This concludes our argument as we have proved that C ∩ D is equal to the connected component of s inside D and that it has the topology of a cone over ∂D ∩ C which is what we claimed. Remark 4.17 We do not have to suppose that D is convex, simply connected would suffice. But we only work with convex sets (boxes) and the denomination "conic topology" originates from the convex case.

In the end the connection algorithm is extremely simple. We just proved that the topology inside these domains is conic, that is C ∩ D can be deformed into a cone over C ∩ ∂D. Therefore the connection algorithm for (convex) simply singular domains is to first compute the points q i of C ∩ ∂D, then choose an arbitrary point p inside D and finally for every q i , connect q i and p by a half branch segment b i = [p, q i ].