Rémi Imbach

Guillaume Moroz

Marc Pouget

Reliable Location with Respect to the Projection of a Smooth Space Curve *

Keywords: Singular curve topology, point location algorithm, geometric approximation AMS subject classifications: 65D99

Consider a plane curve B defined as the projection of the intersection of two surfaces in R 3 or as the apparent contour of a surface. In general, B has node or cusp singular points and thus is a singular curve. Our main contribution is the computation of a data structure answering point location queries with respect to the subdivision of the plane induced by B. This data structure is composed of an approximation of the space curve together with a topological representation of its projection B. Since B is a singular curve, it is challenging to design a method based only on reliable numerical algorithms.

In recent work, the authors show how to describe the set of singularities of B as regular solutions of a so-called ball system suitable for a numerical subdivision solver. Here, the space curve is first enclosed in a set of boxes with a certified path-tracker to restrict the domain where the ball system is solved. Boxes around singular points are then computed such that the correct topology of the curve inside these boxes can be deduced from the intersections of the curve with their boundaries. The tracking of the space curve is then used to connect the smooth branches to the singular points. The subdivision of the plane induced by B is encoded as an extended planar combinatorial map allowing point location. We experimented with our method, and we show that our reliable numerical approach can handle classes of examples that symbolic methods cannot.

Introduction

Let C be a smooth space curve defined as the intersection of two surfaces P = Q = 0, with P , Q real smooth functions in x, y, z. We aim at computing, in a compact domain B0 ⊂ R 2 , the geometry and the topology of B = π (x,y) (C) where π (x,y) : R 3 → R 2 denotes projection into the (x, y)-plane. In general, B is not smooth and has singular points, i.e. points where B has no well defined tangent direction. Generically, the only singular points of a projected curve are transversal crossings of two branches of the curve, called nodes. A special occurrence of our problem is the case where Q = Pz, the derivative of P with respect to z. B is then called the apparent contour of the surface P = 0, and generic singular points of B are nodes and cusps, i.e. projections of points where C has a vertical tangent. Figure 1 shows, for a torus P = 0, its intersection with the surface Pz = 0 in bold, and its apparent contour with cusp and node singularities.

By computing the geometry of B in B0, we mean being able to draw B ∩ B0 with an arbitrary precision. To reach this goal, we use a certified interval path tracker to compute a sequence of boxes (i.e. multi-dimensional extensions of intervals) of width as small as desired, such that each box intersects C, and their union encloses C. The projection of these boxes is thus a geometric approximation of the plane curve B.

We want our topological encoding to be able to answer location queries. The curve B decomposes the box B0 into connected components of dimension two that we call faces. For a point p in B0 but not on B, one wants to find the face to which p belongs. Such queries can be answered by considering a combinatorial encoding of the embedding of the projected curve, i.e. the subdivision of B0 induced by the curve B. We use Combinatorial Maps (CMaps defined in Definition 1.2) and their extensions to the non-connected case, the eXtended Planar Maps (XPMaps defined in Definition 1.3), to encode embeddings of plane curves. In a first step, we focus on singularities and isolate the singular points of B in boxes of width as small as desired.

In the second step we compute the topology in these isolating boxes. Finally, we use the tracking of the space curve to connect the singular points and construct the representation of the embedding.

The work presented here is a first step toward the computation of the topology with reliable approximated geometry of the apparent contour of a smooth manifold of R n . Such a manifold arises naturally in the design of parallel or cable mechanisms, and its apparent contour represents the boundary of the workspace of such a mechanism [START_REF] Merlet | Parallel Robots[END_REF]. The encoding of the topology and the geometry we propose can thus be seen as a reliable tool to validate a robot configuration, to check if the clearance with respect to special configurations is large enough, or to check whether or not the robot passes through a singularity during a motion.

The paper is organized as follows. Section 2 describes how the curve C is enclosed by tracking. Section 3 reviews the encoding of singularities of B by the ball system and shows how the enclosure of C is used to restrict its solution domain. For an apparent contour, an algorithm is presented to determine the type, node or cusp, of a singularity isolated with the ball system. Section 4 is dedicated to the computation of the local topology at special points, i.e. singular points and x-extreme points. Section 5 explains the construction of the XPMap representing the embedding of the curve B. Section 6 reports experiments on the implementation of our numerical approach. The remainder of Section 6 presents previous work, details our contributions, formally defines our geometric and topological representations, and reviews some basics of reliable numerical interval solvers.

Figure 1: Left: a torus P = 0, the curve P = Pz = 0 in bold, its apparent contour, and a zoom zone. The vertical dashed lines show fibers over three singularities. Right: a zoom, with preimages near the projection of cusp and node singularities.

Previous work

State-of-the-art symbolic methods that compute the topology of plane real curves defined by polynomials are based on the Cylindrical Algebraic Decomposition, and use resultant and sub-resultant theory to isolate critical points [START_REF] Hong | An efficient method for analyzing the topology of plane real algebraic curves[END_REF][START_REF] Mourrain | Algebraic issues in Computational Geometry[END_REF]. One advantage of these methods is that they can handle any type of singularity of the curve. The drawbacks are their high complexity as a function of the degree of the curve and the global aspect of the approach: computing the topology in the whole plane or in a small box have almost the same cost.

Numerical methods based on interval arithmetic are able to compute and certify the topology of a non-singular curve [START_REF] Kearfott | An interval step control for continuation methods[END_REF][START_REF] Martin | Certified parallelotope continuation for one-manifolds[END_REF][START_REF] Plantinga | Isotopic approximation of implicit curves and surfaces[END_REF]. One advantage is the local aspect of the approach: the topology can be computed in a small box, and the smaller the box the faster the computation. However, they fail near any singular point of the curve. Isolating singularities of a planar curve f (x, y) = 0 with a numerical method is a challenge, since the set of singular points is described by the non-square system f = ∂f ∂x = ∂f ∂y = 0, and singularities are not necessarily regular solutions of this system. The latter system can be translated into a square system using combinations of its equations with first derivatives [START_REF] Dedieu | Points fixes, zéros et la méthode de Newton[END_REF], and non-regular solutions can be handled through deflation systems [START_REF] Leykin | Newton's method with deflation for isolated singularities of polynomial systems[END_REF][START_REF] Ojika | Deflation algorithm for the multiple roots of a system of nonlinear equations[END_REF], but the resulting systems are usually still overdetermined.

We are not aware of a numerical algorithm that can certify in practice the computation of the topology of any singular curve, but several promising approaches have been presented. The subdivision approach presented in [START_REF] Burr | Complete subdivision algorithms ii: Isotopic meshing of singular algebraic curves[END_REF] is an extension of the singular case of the one in [START_REF] Plantinga | Isotopic approximation of implicit curves and surfaces[END_REF], relying on global non-adaptive separation bounds for algebraic systems. This approach can theoretically certify the topology of any singular curve, but due to the worst-case bounds, this algorithm cannot be practical. A numerical algebraic geometric approach is presented by [START_REF] Lu | Finding all real points of a complex curve[END_REF] using irreducible decomposition, generic projection and plane sweep, deflation and homotopy to compute the topology of a singular curve in any codimension. Even though this work has been implemented in [START_REF] Daniel | Bertini: Software for numerical algebraic geometry[END_REF] 1 , certification of all the algorithm steps appears to be challenge. The numerical approach in [START_REF] Corless | Computing the topology of a real algebraic plane curve whose defining equations are available only "by values[END_REF], based on Bezoutian and eigenvalue computation, can handle singular curves, but even if multiprecision gives accurate results, no certification is provided.

Instead of designing a general numerical method able to handle any singular curve, an alternative is to focus on restricted classes of singular curves. A natural example is when the plane curve to be studied is the projection of a smooth space curve living in a higher dimension. According to the classification of singularities of mappings (see [START_REF] Arnold | Singularities of Differentiable Maps: Volume I: The Classification of Critical Points Caustics and Wave Fronts[END_REF][START_REF] Demazure | Bifurcations and catastrophes: geometry of solutions to nonlinear problems[END_REF][START_REF] Whitney | On singularities of mappings of euclidean spaces. i. mappings of the plane into the plane[END_REF] for example), it appears that the generic singularities of such a projected curve are only transversal intersections of two branches, and in the case of the apparent contour of a surface, ordinary cusps also occur. From an algorithmic point of view, the authors of [START_REF] Delanoue | A numerical approach to compute the topology of the apparent contour of a smooth mapping from R 2 to R 2[END_REF] use these elements to derive an algorithm isolating the singularities arising in generic mappings from R 2 to R 2 . Our problem of isolating the singularities of the projection of a generic algebraic space curve was investigated in [START_REF] Imbach | A certified numerical algorithm for the topology of resultant and discriminant curves[END_REF]. The authors use resultant and sub-resultant theory to represent the singularities as the solutions of a regular bivariate system suited to a branch and bound solving approach. To overcome the drawbacks of resultant and sub-resultant, [START_REF] Imbach | Numeric and certified isolation of the singularities of the projection of a smooth space curve[END_REF] contains studies of the geometric configurations of the space curve that induce singularities on the projected curve, and descriptions of them as the regular solutions of a four dimensional system.

Detailed overview

Our main contribution is the computation of a data structure allowing location of a point with respect to the projection of a smooth space curve. Our data structure is the union of a geometric approximation of the space curve and a topological representation of its projection.

For the geometric approximation of the space curve C and thus its projection B, we compute a sequence of 3-dimensional boxes (Ci) m i=1 enclosing all the connected components of C, that is C ⊂ m i=1 Ci. We use a reliable numerical solver (Algorithm 1) to find intersections of the curve with the boundary of the input box and at least one point on each connected component of C using a critical point method. We then use these points as starting points for a tracking algorithm (Specified in Algorithm 2 and detailed in Appendix A) which is an adapted version of the one in [START_REF] Martin | Certified parallelotope continuation for one-manifolds[END_REF]. Since we want to only use numerical algorithms, we must avoid some degenerate configurations for the systems of critical points and boundary points. Such assumptions, stated in Section 2.1, are satisfied for generic curves and projections.

For the topological representation, the basic numerical tools are the reliable numerical solver and tracker already used for the geometric approximation. It is worth noting that while the topology of the space curve C is directly given by the connected components of boxes enclosing it, the topology of the plane curve B has little to do with its enclosure by the projection of the box enclosure of C. The general idea is first to isolate the singularities and critical points of the plane curve B, then to refine the boxes around those special points until the topology inside these boxes is trivial. Finally, we use smooth path tracking to connect the singularities and the critical points. Since we use numerical algorithms, we require assumptions in Sections 3.1 and 4.2, that are generically satisfied.

The first step is to isolate the singularities of the plane curve B. In previous work [START_REF] Imbach | Numeric and certified isolation of the singularities of the projection of a smooth space curve[END_REF], we have shown how the singularities can be described as the regular solutions of a so-called ball system involving 4 equations in 4 unknowns. The ball system could be solved in R 4 with a reliable numerical solver; however, in four dimensions, this global subdivision approach becomes costly. To overcome this issue, we use the box enclosure of the space curve to restrict the solution domain of the ball system.

The second step is to compute witness boxes for the singularities and critical points of B, that is isolating boxes such that the topology of the curve inside these boxes can be deduced from the intersections of the curve with their boundaries. Special care is devoted to the refinement of boxes such that the intersections of the curve with their boundaries do not eventually occur at the corners.

The last step uses the tracking of the space curve to connect its smooth branches to the witness boxes of the singular and critical points and thus compute a combinatorial map for each connected component. The construction of the extended planar map encoding the embedding of B is incremental on the connected components. The point location algorithm is based on a vertical ray shooting principle but eventually needs to compute intersections of the space curve C with planes.

We implemented the presented algorithms and tested them to compute the topology of apparent contours of algebraic surfaces of degrees up to 15. Our experiments show that our method can handle classes of examples symbolic methods cannot handle, and that multi-precision arithmetic is needed for such difficult examples. More specifically, for isolation of the singularities, the efficiency of our approach to restrict the solving domain of the ball system is demonstrated.

Notation and definitions

In this paper, real intervals are connected sets [a, b] with a, b ∈ R ∪ {±∞} and a ≤ b. Lowercase boldface letters denote real intervals and uppercase boldface letters denote boxes, that is, vectors of intervals. Let x be a real interval, let l(x) denote its lower bound, u(x) its upper bound and w (x) its width, defined as u

(x)-l(x) if x is bounded and ∞ otherwise. If x is bounded, m(x) denotes its midpoint. If ∈ R + , x holds for [m(x) -w (x) 2 , m(x) + w (x) 2 ]
. Let X = (x1, . . . , xn) be a box, ∂X denotes the boundary of X and i(X) = X \ ∂X its interior, where A \ B is the set of elements of A not in B. The width w (X) is defined as max 1≤i≤n w (xi) and the midpoint m(x) as (m(x1), . . . , m(xn)) if X is bounded. If ∈ R + , X is defined to mean ( x1, . . . , xn).

Uppercase letters denote sets of boxes. A domain of R n is a set defined as X \ ( Y∈Y i(Y)) where X is a box in R n and Y a possibly empty set of boxes, of R n . If X is a domain of R n , ∂X denotes its boundary.

For a differentiable real function P of the variables x1, . . . , xn, Px i denotes its partial derivative with respect to xi, and Px i x j its derivative with respect to xi and xj. Let P1, . . . , Pn be n real functions of x1, . . . , xn, a solution x of P1 = . . . = Pn = 0 is regular if the Jacobian matrix A = [(Pi)x j ] evaluated at x has full rank; otherwise x is singular.

Cursive letters denote sets of points. We mostly work with points, boxes and curves in R 2 , R 3 or R 4 . We use the following naming scheme: objects in R 2 are named with the letter B, in R 3 with the letter C and in R 4 with the letter D.

A graph is a triplet G = (V, E, I) where V, E are two finite sets whose elements are called vertices and edges, respectively, and

I : E → V × V is called the incidence relation of G. G is directed (resp. non-directed ) if
images of edges by I are seen as ordered sets (respectively non-ordered pairs).

Geometrical and topological representations

Our first goal is to compute a geometrical approximation of the projected curve B in a box B0. Let us formalize this notion as follows.

v 7 v 5 v 3 v 1 v 2 v 4 v 8 v 6 v 11 v 10 v 9 h - 1 h + 13 h - 10 h - 5 h + 4 h + 9 h - 7 h - 11 h + 10 h + 11 h - 13 h + 14 h - 14 h + 1 h + 12 h - 12 h - 9 h - 8 h + 7 h + 6 h + 8 h + 5 h - 3 h - 2 h + 2 h + 3 h - 4 h - 6 H0 = ∅, H1 = {h + 1 , h - 1 , h + 2 , h - 2 , h + 13 , h - 13 , h + 14 , h - 14 }, H2 = {h + 3 , h - 3 , . . . , h + 12 , h - 12 } edges: (h + i , h - i ) for 1 ≤ i ≤ 14, vertices: v1 = (h - 1 , h + 14 , h - 13 ),. . . faces: f0 = ∅, f1 = (h + 13 , h + 14 ), f2 = (h + 1 , h + 2 , h - 14 ), f3 = (h + 1 , h - 13 , h - 2 ), f4 = (h + 3 , h - 4 , . . . , h - 12 ), f5 = (h - 3 , h + 12 , h - 11 , h + 4 ), f6 = (h + 5 , h - 10 ), f7 = (h + 6 , h - 9 , h + 8 , h - 7 ). ext: (H 1 , f1), (H 2 , f4) cont: (f1, f0), (f4, f3)
1.1 (δ-approximation) Let X be a subset of R n . A sequence of boxes (Xi) m i=1 is a δ-approximation of X if X ⊂ m i=1
Xi, and for all 1 ≤ i ≤ m, w (Xi) ≤ δ and X ∩ Xi = ∅.

We will consider in particular δ-approximations of C in C0 = B0 × R or in a domain C0. Obviously the projection of a δ-approximation of C in C0 is a δ-approximation of B in B0, in particular it yields a drawing of B with the guarantee that any point of the drawing is at most at distance δ from the original curve. The left part of Figure 2 represents a δ-approximation of the curve C in C0 and its projection, for the torus example.

We also want to encode the topology of the curve B ∩ B0 as an embedding of a graph G, to be able to answer location queries for points of B0. This is achieved by computing a set V of vertices representing points on the curve with at least one vertex per connected component, and such that (B ∩ B0) \ V is a set of smooth curves identified as the set E of edges. The right part of Figure 2 shows an embedding of such a graph G when B ∩ B0 is the apparent contour of the torus. The embedding of G defined by B ∩ B0 is encoded by an extended planar map. More precisely, each connected component of (B ∩ B0) ∪ ∂B0 is encoded by a combinatorial map, and the inclusion of a component in a face of another is encoded by the extended map. We thus recall the definition of these combinatorial structures. Definition 1.2 ( [START_REF] Sergei | Graphs on surfaces and their applications[END_REF][START_REF] Köthe | Xpmaps and topological segmentation-a unified approach to finite topologies in the plane[END_REF]) A Combinatorial Map (CMap) is a triplet (H, σ, α) where H is a set of half-edges, σ is a permutation on H and α an involution on H.

An edge of G is associated with a cycle of α, and a vertex with a cycle of σ. Cycles of σ encode counter-clockwise orderings of outgoing half-edges around vertices. The cycles of the permutation ϕ = σ -1 • α describe the faces of the combinatorial map. All such face cycles are in counter-clockwise order, except for one which is called the exterior face of the combinatorial map.

When (B ∩ B0) ∪ ∂B0 has several connected components, it remains to encode the containment relationship between the exterior face of each CMap within a non-exterior face of another CMap. In [START_REF] Köthe | Xpmaps and topological segmentation-a unified approach to finite topologies in the plane[END_REF] In our case, the infinite face is R 2 , Figure 3 shows an XPMap characterizing (B ∩ B0) ∪ ∂B0 as an embedding of the apparent contour of the torus. The cycle f1 = (h + 13 , h + 14 ) is the exterior face of the connected component of (B ∩ B0) ∪ ∂B0 containing ∂B0. This exterior face f1 is contained in the infinite face f0 creating a hole in it. More gererally, any non-infinite face of an XPMap can be described by exactly one clockwise cycle of half-edges and possibly several inner counter-clockwise cycles of half-edges corresponding to the exterior face of the CMaps it contains.

Interval arithmetic tools

The certification of our algorithms is based on interval arithmetic (see [START_REF] Kearfott | Rigorous global search: continuous problems. Nonconvex optimization and its applications[END_REF][START_REF] Ramon E Moore | Introduction to interval analysis[END_REF][START_REF] Neumaier | Interval methods for systems of equations[END_REF][START_REF] Stahl | Interval Methods for Bounding the Range of Polynomials and Solving Systems of Nonlinear Equations[END_REF]), that is a way of computing with intervals (whose endpoints are floating numbers) instead of computing with floating numbers, while carefully handling rounding to overcome numerical approximations that naturally occur with floating number arithmetic.

Convergent interval functions

If F be a function, an inclusion function for F is a function defined on boxes such that the image of a box X is a box Y ⊇ {F (x)|x ∈ X}. As a consequence, if 0 / ∈ Y, one can certify that F does not vanish in X. An inclusion function for F is convergent if for any decreasing sequence of boxes converging to a point {p}, the image of the sequence converges to {F (p)}.

In this article, we consider C k -smooth functions, that is functions that are k ≥ 1 times differentiable and with continuous partial differentials. In addition, we also assume that we are given convergent inclusion functions for such functions and all their differentials up to order k.

Criteria for existence and uniqueness of solutions

Let P1, . . . , Pn be n functions C 1 in n unknowns and S = {P1 = 0, . . . , Pn = 0} the associated system of equations. A box X ⊂ R n isolates a solution of S if there exists a unique x ∈ X such that P1(x) = . . . = Pn(x) = 0.

Several criteria can be found in the interval arithmetic literature that certify existence and uniqueness of a solution of S in a box, see for instance [START_REF] Kearfott | Rigorous global search: continuous problems. Nonconvex optimization and its applications[END_REF][START_REF] Ramon E Moore | Introduction to interval analysis[END_REF][START_REF] Neumaier | Interval methods for systems of equations[END_REF][START_REF] Stahl | Interval Methods for Bounding the Range of Polynomials and Solving Systems of Nonlinear Equations[END_REF]. Most of them are based on the Brouwer fixed point theorem, and use interval Newton operators that contract a box around a solution.

Letting F be the multivariate function with components P1, . . . , Pn and X an interval of R n , interval Newton operators are of the form N (X) = y -V , where y ∈ X, and V is a box containing the set {v | ∃x ∈ X such that J(x)v = F (y)}, where J is the Jacobian of F . Among other interval Newton operators is the Krawczyk operator [START_REF] Krawczyk | Newton-Algorithmen zur Bestimmung von Nullstellen mit Fehlerschranken[END_REF][START_REF] Ramon E Moore | Introduction to interval analysis[END_REF], that takes y as the midpoint of X and an approximate inverse of the Jacobian of F in y to determine the box V . Let us denote by KS the Krawczyk operator for the system S, and recall that i(X) is the interior of X. Namely,

KS(X) = y -YF (y) + {I -YJ(X)}(X -y),
where Y is a nonsingular real matrix approximating the inverse of the Jacobian matrix of F at m(X). The Krawczyk operator satisfies:

• KS(X) ⊂ i(X) ⇒ X contains a unique solution of S and the sequence X (0) = X, X (i+1) = KS(X (i) ) converges asymptotically quadratically to this solution.

• KS(X) ∩ X = ∅ ⇒ there is no solution of S in X.

Reliable numerical isolation

Interval evaluation, the Krawczyk operator and bisection of boxes can be used together to design a simple reliable numerical method to isolate all solutions of S in a given initial box X0. Such methods are described for instance in [START_REF] Kearfott | Rigorous global search: continuous problems. Nonconvex optimization and its applications[END_REF][START_REF] Ramon E Moore | Introduction to interval analysis[END_REF][START_REF] Neumaier | Interval methods for systems of equations[END_REF][START_REF] Stahl | Interval Methods for Bounding the Range of Polynomials and Solving Systems of Nonlinear Equations[END_REF], and called interval branch and bound algorithms or subdivision methods. These algorithms successfully isolate all the solutions of S in a bounded domain, provided they are Algorithm 1 IsolateSols(S,X 0 ,δ)

Input: A bounded box X0 of R n , a system S of n functions C 1 in n unknowns with only regular solutions in X0, the Krawczyk operator KS, and a non-negative real number δ. Output: Two sets X sol , X ind of boxes such that:

• boxes of X sol ∪ X ind are pairwise disjoint,

• x ∈ X0 is a solution of S ⇒ ∃X ∈ X sol ∪ X ind s.t. x ∈ X, • X ∈ X sol ⇒ X ⊂ i(X0) and KS(X) ⊂ i(X),
• X ∈ X ind ⇒ X ∩ ∂X0 = ∅ and KS(X) ⊂ i(X) and w (X) < δ.

regular and in the interior of the domain. Note that in the case of a polynomial system it is possible to extend such branch and bound methods to unbounded initial domains; see [START_REF] Neumaier | Interval methods for systems of equations[END_REF]Section 5.6] or [START_REF] Stahl | Interval Methods for Bounding the Range of Polynomials and Solving Systems of Nonlinear Equations[END_REF]Section 5.10].

In the following, we consider the procedure IsolateSols with the specifications given in Algorithm 1. We briefly explain the output of this algorithm, and refer to [START_REF] Imbach | A Subdivision Solver for Systems of Large Dense Polynomials[END_REF] for details. The algorithm does not identify solutions that are exactly on the boundary of the input domain X0, and uses δ-inflation. A consequence is that a solution on ∂X0 can only be isolated in a box containing a part of the boundary in its interior; the refinement of such a box is stopped due to the threshold on the width δ, and the box is output in X ind . Similarly, for a regular solution not on, but near the boundary (even outside X0), that is at distance less than δ from ∂X0, the algorithm may return a box in X ind . On the other hand, with the additional assumption that there is no solution on the boundary, setting δ = 0 certifies the isolation of all solutions in the domain X0, that is, that X ind is empty.

Enclosing the Space Curve C

In this section, we introduce Algorithm 3, computing a δ-approximation of C. This goal is achieved by first finding at least one point on each connected component of C, then using these points as initial points for a certified path tracker. Let us first introduce the notions of x-or y-critical points and boundary points.

= Q = R = 0, where R = PyQz -PzQy. Similarly, p ∈ C is y-critical if it is the solution of the system P = Q = R = 0 where R = PxQz -PzQx. Definition 2.2 If C is a box or a domain in R 3 , a point of C ∩∂C is called a boundary point.
Isolating boxes for boundary points and x-critical points of C are computed with the procedure IsolateSols defined in Algorithm 1. We introduce the notion of implicit points to manipulate a point known as the unique solution of a system S in a box C = x × y × z. Definition 2.3 Let P, Q define the curve C, let E ∈ R[x, y, z], and let S be the system

P = Q = E = 0. The ordered pair (C, S) is an implicit point of C if KS(C) ⊂ i(C).
If p is the unique solution of S in C, we say that p is implicitly defined by (C, S).

In the special case where E = x -x0 with x0 ∈ R, we say that

(x0 × y × z, S) is an implicit point of C if KS 0 (y×z) ⊂ i(y×z), where S0 = {P (x0, y, z) = Q(x0, y, z) = 0}.
Similarly, in the special case where E = y-y0 with y0 ∈ R, we say that

(x×y0 ×z, S) is an implicit point of C if KS 0 (x×z) ⊂ i(x×z), where S0 = {P (x, y0, z) = Q(x, y0, z) = 0}.
In Sec. 2.2, we characterize a connected component of C ∩ C0: it is diffeomorphic either to a circle and contains at least two x-critical points, or to [0, 1] and its extremities are boundary points. We also show how to obtain these points as implicit points. In Sec. 2.3, we specify Algorithm 2, our certified path tracker, and in Sec. 2.4, we specify Algorithm 3 for computing a δ-approximation of C. In Sec. 2.1, we state the assumptions allowing our approach to be correct and terminating.

Assumptions

Recall that B0 = (x0, y0) and C0 = B0 × R. We assume the functions P and Q are C3 functions2 and that convergent inclusion functions are given for P, Q and all their derivatives up to order 3. We make the following assumptions on P, Q and C: (A5) Over a point of ∂B0, P = Q = 0 has only one solution and no solutions above its corners.

(A1) The curve C is smooth above the box B0. (A2) C is compact over B0 and z0 is a bounded interval such that C ⊂ B0 × z0.

Connected components of C ∩ C 0 and initial boxes

The following proposition characterizes the topology of the connected components of C ∩ C0. Proof: The first part of the proposition is straightforward. One dimensional manifolds are diffeomorphic either to ]0, 1[, or to ]0, 1], or to [0, 1], or to a circle. Let C k be a connected component of C. From assumption (A2), it is compact, hence it is diffeomorphic either to [0, 1], or to a circle. Suppose first C k has an intersection with ∂C0. From (A3), this intersection is a point, hence C k is diffeomorphic to [0, 1] and has a second intersection with ∂C0. Suppose now that C k does not intersect ∂C0. Hence it is diffeomorphic to a circle and, since it is compact, it has minimum and maximum x-coordinates. Assertion (b) follows.

As a direct consequence of Prop. 2.1, the following corollary gives a constructive characterization of a set containing at least one point on each connected component of C ∩ C0.

Corollary 2.1 Consider the following systems of equations:

(S1) P (l(x0), y, z) = Q(l(x0), y, z) = 0, for y ∈ y0 and z ∈ R. (S2) P (u(x0), y, z) = Q(u(x0), y, z) = 0, for y ∈ y0 and z ∈ R. (S3) P (x, l(y0), z) = Q(x, l(y0), z) = 0, for x ∈ x0 and z ∈ R. (S4) P (x, u(y0), z) = Q(x, u(y0), z) = 0, for x ∈ x0 and z ∈ R. (S5) P (x, y, z) = Q(x, y, z) = R(x, y, z) = 0, for (x, y, z) ∈ C0.
Assuming (A3), (A4), the set of solutions of S1, . . . , S4 is finite and is the set of boundary points of C ∩ C0, the set of solutions of S5 is finite and is the set of x-critical points of C ∩ C0. Then assuming (A2), the set of solutions of S1, . . . , S5 is a finite set of points in B0 × z0 containing at least one point on each connected component of C.

The solutions of S5 are obtained by calling IsolateSols(S5,B0 × z0,0); from assumption (A4) this process terminates and the outputs C sol , C ind are such that C ind = ∅ (since S5 has no solution on ∂C0) and C sol contains boxes isolating the solutions of S5. For each solution C ∈ C sol , the implicit point (C, S5) defines an x-critical point of C.

The solutions of S1 are obtained by calling IsolateSols(S1,y0 × z0,0). From assumptions (A3) and (A5), S1 has a finite number of solutions on y0 × R that are regular and has no solutions on ∂y0 × R. Hence the procedure terminates, and the two obtained sets B sol , B ind are such that B ind = ∅ and B sol contains two dimensional boxes isolating the solutions of S1. A solution B = (y, z) in B sol defines an implicit point (C, {P, Q, E}) where C = l(x0) × y × z and E is the polynomial x -l(x0), as described in Definition 2.3. When needed, C is contracted with l(x0)×KS 1 (π (y,z) (C)). The solutions of S2, S3 and S4 are obtained similarly.

Certified numerical path-tracking

Our path-tracking procedure Track is specified in Algorithm 2, while its description is given in Appendix A. The tracker is used in Algorithm 3 to compute the connected components of C ∩ C0, and later in Section 5.1.2 to deduce the topology of B.

Computing a δ-approximation of C

In the following, we assume that P, Q satisfy assumptions (A1) and (A2) of Section 2.1: C is smooth and compact in C0. Based on the procedure Track, Algorithm 3 computes a δ-approximation of C ∩ C0.

Algorithm 2 Track(< P, Q >,C 0 ,(C 0 , S 0 ),{(C j , S j )} j ,δ)

Input: A system P = Q = 0 defining a smooth compact curve C in the domain C0,
an implicit point (C 0 , S 0 ) of C, a finite set {(C j , S j )}j of implicit points of C containing (but not restricted to) the boundary points of C in C0 and (C 0 , S 0 ), δ > 0. Output: A δ-approximation of the connected component of C∩C0 containing (C 0 , S 0 ) and the set C on of implicit points of {(C j , S j )}j that are on the same connected component as (C 0 , S 0 ).

Algorithm 3 Compute a δ-approximation of C in C 0 Input: A system P = Q = 0 and an initial domain C0 = B0 × R defining a smooth curve C ∩ C0, the set C b of implicit points defining the boundary points of C ∩ C0, and the set C x of implicit points defining the x-critical points of C ∩ C0, δ > 0. When L b is empty, each connected component that is diffeomorphic to [0, 1] has been approximated, and the implicit points in L x define the x-critical points that belong to the connected components that are diffeomorphic to circles.

Output: A sequence of boxes (Ci) m i=1 that is a δ-approximation of C ∩ C0. 1: Let L b (resp. L x ) be a list containing elements of C b (resp. C x ) 2: k ← 0, m k ← 0 3: while L b = ∅ do 4: k = k + 1, (C, S) ← pop front(L b ) 5: ((Ci) m k i=m k-1 +1 ,C on ) ← Track(< P, Q >,C0,(C, S),C b ∪ C x ,δ) 6: remove C on ∩ C b from L b and C on ∩ C x from L x 7: while L x = ∅ do 8: k = k + 1, (C, S) ← pop front(L x ) 9: ((Ci) m k i=m k-1 +1 ,C on ) ← Track(< P, Q >,C0,(C, S),C b ∪ C x ,δ)
The connected components of C that are diffeomorphic to a circle are approximated in the while loop beginning at Step 7 of Algorithm 3. Recall that each connected component of C ∩ C0 that is diffeomorphic to a circle contains at least two x-critical points. Each time the loop is performed, the size of L x decreases, and when L x is empty, each connected component of C has been approximated.

Isolating Singularities of B

When C is defined by two analytic maps P, Q, [START_REF] Imbach | Numeric and certified isolation of the singularities of the projection of a smooth space curve[END_REF] describes, under genericity conditions on P, Q, the type of singularities arising in the projection B: they are only nodes (two branches of C induce a self intersection in B), or cusps (C has a vertical tangent).

[12] also introduces a system called ball system whose solutions are in one-to-one correspondence with the singularities of B, and shows that the ball system only has regular solutions if and only if singularities of B are either nodes or ordinary cusps.

We first restate the assumptions and the main results of [START_REF] Imbach | Numeric and certified isolation of the singularities of the projection of a smooth space curve[END_REF]. We then show how an enclosure of C helps to restrict the domain where the ball system is solved while ensuring that all cusps and nodes are obtained. Then, we present Algorithm 4, that decides, for a given solution of the ball system, if the corresponding singularity is an ordinary cusp or a node when B is an apparent contour. Within this section, Σ denotes the set of singular points of B ∩ B0.

Assumptions

Consider the following assumptions:

(A6) For any (α, β) in B0, the system P (α, β, z) = Q(α, β, z) = 0 has at most 2 real roots, counting multiplicities.

(A7) There are finitely many points (α, β) in B0 such that P (α, β, z) = Q(α, β, z) = 0 has 2 real roots, counting multiplicities.

(A8) The singularities of the curve B in B0 are either nodes or ordinary cusps.

Notice that the Thom Transversality Theorem implies that (A1) . . . , (A8) hold for generic smooth maps P, Q defining C (see [8, Theorem 3.9.7 and §4.7]).

Ball system

Following a geometric modelling, [START_REF] Imbach | Numeric and certified isolation of the singularities of the projection of a smooth space curve[END_REF] defines a 4 dimensional system whose solutions map to the singularities of B. In this modelling, two solutions (x, y, z1) and (x, y, z2) of P = Q = 0 (or P = Pz = 0) are mapped to the point (x, y, c, r2) with c = (z1 + z2)/2 and r2 = (z1 -z2) 2 . Figure 4 illustrates this mapping for singularities of the apparent contour of a torus. We review the main results of [START_REF] Imbach | Numeric and certified isolation of the singularities of the projection of a smooth space curve[END_REF]. Lemma 3.1 ([12, Lemma 4]) Let P, Q be two analytic functions in x, y, z satisfying the Assumptions (A1), (A2), (A6) and (A7), and let S be the set of solutions of the socalled ball system:

                   1 2 (P (x, y, c + √ r2) + P (x, y, c - √ r2)) = 0 1 2 √ r2 (P (x, y, c + √ r2) -P (x, y, c - √ r2)) = 0 1 2 (Q(x, y, c + √ r2) + Q(x, y, c - √ r2)) = 0 1 2 √ r2 (Q(x, y, c + √ r2) -Q(x, y, c - √ r2)) = 0 (1) in B0 × R × R + . Then π (x,y) (S) = Σ
, where π (x,y) is the projection from R 4 to the (x, y)-plane. In addition, the solutions of the ball system with r2 = 0 map to cusps in 2d and the curve in 3d has a vertical tangent (collinear to the z-axis); whereas the solutions with r2 = 0 map to nodes.

Solution domain

In [START_REF] Imbach | Numeric and certified isolation of the singularities of the projection of a smooth space curve[END_REF], the ball system is solved within the box B0 × R × R + with a subdivision solver. Using the δ-approximation (Ci) 1≤i≤m of C computed by Algorithm 3, we propose to reduce significantly the search domain for the singularities. Indeed, given a singular point σ of B, there exists 1 ≤ i ≤ m such that σ ∈ Bi, where Bi = π (x,y) (Ci). Hence it is possible to isolate all singularities by solving the ball system within Bi × R × R + , for 1 ≤ i ≤ m. In addition, Proposition 3.1 shows how to bound the solution domain in the c and r2 components.

Proposition 3.1 Let (Ci) 1≤i≤m be a δ-approximation of C. For 1 ≤ i ≤ m, let Ci = (xi, yi, zi), Bi = π (x,y) (Ci) = (xi, yi), and for 1 ≤ i < j ≤ m, let Bij = (xij, yij) = Bi ∩ Bj
and consider the sets:

• Di = (xi, yi, zi, [0, ( w (z i ) 2 ) 2 ]), • Dij = (xij, yij, (z i +z j ) 2
, (

(z i -z j ) 2 ) 2 )
Then, under the Assumptions (A1), (A2), (A6) and (A7), all solutions of the ball system lie in ( 1≤i≤m Di) ∪ ( 1≤i<j≤m Dij).

Proof: Let p = (xp, yp, cp, rp) ∈ B0 × R × R + be a solution of the ball system and σ = π (x,y) (p) the corresponding singularity in Σ. From Assumption (A8), σ is either an ordinary cusp or a node. Suppose first it is an ordinary cusp. Then rp = 0, and σ is the projection of a single point p = (xp, yp, cp) of C. Hence there exists 1 ≤ i ≤ m such that p ∈ Ci. Thus, we have cp ∈ zi and p ∈ Di (see Figure 5).

Suppose now σ is a node. Then rp > 0, and σ is the projection of two points

p-= (xp, yp, cp - √ rp) and p+ = (xp, yp, cp + √ rp) of C. Hence there exist 1 ≤ i ≤ m and 1 ≤ j ≤ m such that p-∈ Ci and p+ ∈ Cj. If i = j, we have cp ∈ zi and rp ∈ [0, ( w (z i ) 2
) 2 ], and finally p ∈ Di. If i = j (this case is illustrated in Figure 5), cp lies in

z i +z j 2
that is the center of the two intervals zi and zj, and rp lies in (

z i -z j 2
) 2 , that is the square of the corresponding radius.

Singularities of an apparent contour

For a generic 3-dimensional curve, the singularities of its projection B are only nodes. When the curve B is the apparent contour of a surface, its singularities generically also include cusps. We introduce Algorithm 4 to distinguish these two types of singularities.

Let D = (x, y, c, r) be a box isolating a solution p = (x, y, c, r) ∈ D of the ball system and let σ = π (x,y) (p) be the associated singularity of B. If 0 / ∈ r then r = 0 and σ is a node. Otherwise, σ can be either a cusp or a node. Recall that σ is an ordinary cusp of B only if it is the projection of a point of CP ∩Pz that has a vertical tangent (collinear to the z-axis). In other words, if σ = (α, β) is a cusp, there exists a unique

γ ∈ R such that σ = π (x,y) (α, β, γ) and P (α, β, γ) = Pz(α, β, γ) = Pzz(α, β, γ) = 0. According to assumption (A6), γ is a triple root of P (α, β, z), that is, Pzzz(α, β, γ) = 0.
The regularity of the curve CP ∩Pz thus implies that (α, β, γ) is a regular solution of the system P = Pz = Pzz = 0 (see [START_REF] Imbach | A certified numerical algorithm for the topology of resultant and discriminant curves[END_REF]Lemma 10]). Noting that K (P,Pz ,Pzz ) is the Krawczyk operator for the latter system, the test on Line 2 of Algorithm 4 is eventually true for a small enough box whose projection contains a cusp. Algorithm 4 Singularity types for an apparent contour

Input: P in Q[x, y, z], a box D = (x, y, c, r) s.t K b (D) ⊂ i(D)
, where K b is the Krawczyk operator for the ball system. Output: The type of the singularity contained in π (x,y) (D) = (x, y).

1: while 0 ∈ r do 2:

if K (P,Pz ,Pzz ) ((x, y, c)) ⊂ i((x, y, c)) then return cusp 3:

D = K b (D)
4: return node and π (x,y) (p) is an ordinary cusp.

Topology at Special Points of B

In this section, we compute witness boxes for the singularities and critical points of B, that is, isolating boxes such that the topology of the curve inside these boxes can be deduced from the intersections of the curve with their boundaries. Special care is taken to refine the boxes such that the intersections of the curve with their boundaries do not eventually occur at the corners. The case of a node is detailed; the other cases are only sketched, since they are similar.

Boundary points

We

Preprocessing of nodes, cusps and x-extreme points

We add the following assumptions.

(A9) The two points of C above a node of B are neither x-critical nor y-critical points of C.

(A10) P = Q = R = 0 has finitely many regular solutions in C0 and none in ∂C0, and they are the y-critical points of the curve C.

Note that a cusp of B corresponds to a point of C with a vertical tangent (collinear to the z-axis), that is both an x and y-critical point of C. With all our assumptions, the x or y-critical points of C are in i(C0), and • x and y-critical points of C are in one-to-one correspondence with cusps of B.

• x-critical points of C that are not y-critical are in one-to-one correspondence with smooth x-critical points of B.

• y-critical points of C that are not x-critical are in one-to-one correspondence with smooth y-critical points of B. The first step is to compute implicit points for y-critical points of C, by solving the system P = Q = R = 0 with IsolateSols on B0 ×z0. Assumption (A10) ensures the termination of the process with input = 0. Then, the 4-dimensional boxes isolating nodes and cusps (the solutions of the ball system) and the 3-dimensional boxes isolating x and y-critical points (the solutions of P = Q = R = 0 and P = Q = R = 0) are refined until their projections in the (xy)-plane satisfy

• a box of a cusp overlaps exactly one x-critical and one y-critical point,

• a box of a node overlaps no x nor y-critical point,

• a box of an x-critical point that is not overlapping a cusp, does not overlap any y-critical point.

An x-critical point of B that is smooth, i.e. is not also y-critical, is called x-extreme. A point of B a points that is an x-extreme, a node or a cusp is called a special point. Similarly, as in [START_REF] Imbach | A certified numerical algorithm for the topology of resultant and discriminant curves[END_REF], we define a witness box for a special point of the curve B. Definition 4.1 A witness box for a special point (x-extreme, node or cusp) of the curve B is a box containing this point, such that the topology of the curve inside the box is the one of the graph connecting its center to the crossings of the curve on its boundary.

This definition implies that a witness box has 4 crossings of the curve on its boundary for a node and 2 crossings for a cusp or an x-extreme point. The right part of Figure 6 shows a witness box for a node.

Topology at a node singularity

1: Let B = (xB, yB) 2: (X sol 1 , X ind 1 ) = IsolateSols(P (x = l(xB), y, z) = Q(x = l(xB), y, z) = 0,yB × z0,w (yB)/4) 3: (X sol 2 , X ind 2 ) = IsolateSols(P (x = u(xB), y, z) = Q(x = u(xB), y, z) = 0,yB × z0,w (yB)/4) 4: if |X ind 1 |> 0 or |X ind 2 |> 0 then 5: return False 6: (X sol 3 , X ind 3 ) = IsolateSols(P (x, y = l(yB), z) = Q(x, y = l(yB), z) = 0,xB × z0,0) 7: (X sol 4 , X ind 4 ) = IsolateSols(P (x, y = u(yB), z) = Q(x, y = u(yB), z) = 0,xB × z0,0) 8: if not |X sol 1 |+|X sol 2 |+|X sol 3 |+|X sol 4 |= 4 then 9:
return False 10: return True witness box when the curve B crosses its boundary 4 times, or equivalently the curve C crosses ∂B × R 4 times. To avoid the problem of crossings on the corners of B, the refinement is performed such that its x-coordinate is exponentially smaller than its y-coordinate, so that B will eventually cross only the left and right sides of B and far from the corners (see the right part of Figure 6). Suppose Algorithm 5 returns True. This implies that B crosses the boundary of B in exactly four distinct points. Since C has neither x nor y-critical point above B, the two branches of C, those projections pass through the node, are x and y-monotone above B. The projections of these two branches thus cross the boundary of B in exactly four distinct points. In addition, there cannot be any other branch of C, since it would either generate more crossings on the boundary or the existence of critical points in the box; thus (ii) follows.

Proof of (iii). If there exists a box X ∈ X ind 1 , it should include one of the points in l(xB) × ∂yB and its width should be less than w (yB)/4. Since the curve B does not intersect N W (∂B) ∪ SW (∂B) such a box cannot contain any solution; thus X ind 1 is empty. Similarly, X ind 2 is empty and the 2 solutions on E(∂B) (resp. W (∂B)) are reported in X sol 1 (resp. X sol 2 ). In addition, the curve B intersects neither S(∂B) nor N (∂B), and X sol 3 and X sol 4 are empty. Hence the number of reported solutions is 4, and the algorithm returns True.

Note that Algorithm 5 may return a false negative: a box that is a witness box may not be classified as such. This can happen when the curve crosses the box near its corners. The idea of Algorithm 6 is to refine the box of a node to avoid such a case, that is, such that property (iii) of Proposition 4.1 eventually holds. The sequence B n of boxes constructed in Algorithm 6 is illustrated in the right part of Figure 6.

m > n, B m ⊂ B n σ . Note also that N W (∂B n )∪N E(∂B n )∪ SW (∂B n ) ∪ SE(∂B n ) ⊂ B n σ .
The boxes B n and B both contain the node σ in their interior, so the condition B n ⊆ B of Line 9 will be true for any n large enough.

It remains to show that for n large enough Algorithm 6 will report B n as witness. We will show that B n satisfies the sufficient condition (iii) of Proposition 4.1.

Since B × R does not contain any x and y-critical points of C, for n large enough, one can assume that C ∩ (B n σ × R) has two connected components C 1 and C 2 that are x and y-monotone. These two branches project onto two curves B 1 = π (x,y) (C 1 ) and B 2 = π (x,y) (C 2 ) that are x-monotone, cross at σ and have bounded slopes. On the other hand, the aspect ratio bn an = 2 n b 0 a 0 of the box B n is increasing. For n large enough, B 1 (resp. B 2 ) thus crosses E(∂B) and W (∂B) exactly once, and there is no other

intersection on N W (∂B n ) ∪ N (∂B n ) ∪ N E(∂B n ) ∪ SW (∂B n ) ∪ S(∂B n ) ∪ SE(∂B n ).
Condition (iii) of Proposition 4.1 is thus satisfied, and the returned box B n is a witness box.

Topology at a cusp singularity

Algorithms 5 and 6 are easily adapted for a cusp. The input of Algorithm 6 should then be a box D containing a unique solution of the ball system that projects into a cusp σ and such that π (x,y) (D) × R does not contain any other x or y-critical point of C than the point projecting onto the cusp. The input of Algorithm 5 should then be a box B containing a unique cusp of B, such that B × R does not contain any other x or y-critical point of C other than the point projecting onto the cusp. The test for the number of reported solutions at Line 8 should be with the value 2 instead of 4. The output will be true when B intersects E(∂B) exactly once, W (∂B) exactly once, and

does not intersect N W (∂B) ∪ N (∂B) ∪ N E(∂B) ∪ SW (∂B) ∪ S(∂B) ∪ SE(∂B).
For the proof of correctness, the same arguments hold when the limit of the tangent to B at the cusp is not collinear with the y-axis, since in this case the slope of the branch is bounded. When this limit is the y-axis, the same algorithm with the variables x and y swapped will behave as above. In other words, the box B n is elongated in the x-direction instead of the y-direction, and intersections of the curve with its boundary will eventually appear on the north or south sides and far from the corners. The solution is thus to run the two algorithms in parallel, and stop as soon as one has identified a witness box.

Topology at an x-extreme point

For an x-extreme point, the method is similar to the one for a cusp with a limit of the tangents collinear with the y-axis. Indeed, the tangent to the curve at an x-extreme point is collinear with the y-axis, so that the curve is locally y-monotone. The box B n is thus elongated in the x-direction, and intersections of the curve with its boundary will eventually appear on the north or south sides and far from the corners.

Global Topology of B as an Embedded Graph

In this section, the certified tracking of the curve C ∩ C0 together with the local topology at the special points of its projection B ∩ B0 are combined to compute the global topology of B ∩ B0. We use an XPMap to encode this topology and design a point location algorithm, that is, given p ∈ B0 \ B, find the connected component of B0 \ B to which p belongs.

In Section 5.1, we compute a graph G such that B ∩ B0 is an embedding of G. The vertices of G are the special points computed in Section 4, and the edges are computed by tracking between these points. In addition, to restrict the point location algorithm to the box B0, we consider the curve B ∂ = (B ∩ B0) ∪ ∂B0 as the embedding of a graph G ∂ . The graph G ∂ and its embedding B ∂ have in general several connected components. Section 5.2 shows how to compute a CMap encoding the topology of one component. The only geometric task for this is to order branches of B around node singularities. Section 5. [START_REF] Burr | Complete subdivision algorithms ii: Isotopic meshing of singular algebraic curves[END_REF] shows how to answer location queries with respect to one connected component of B ∂ . In Section 5.4, the XPMap encoding the topology of B ∂ is constructed, and the point location algorithm is generalized to this structure.

Computing a graph of which B ∩ B 0 is an embedding

Let B x and B n be the sets of witness boxes for the x-extreme points or cusps, and let the nodes of B ∩ B0 be computed in Section 4. We define V x and V n as the cylinders above witness boxes: Note that implicit points of C b define connections of C ∩ C0 in C0. During the computation of witness boxes as described in Section 4, the connections are also computed. We thus assume that for a box C ∈ V x ∪ V n , the connections of C ∩ C0 in C are given by the set connect(C) = {. . . , (Ci, Si), . . .} of disjoint implicit points (i.e. Ci are pairwise disjoint). If C ∈ V b , we let connect(C) = {(C, S)} where S is the system defining the boundary containing C (see Corollary 2.1).

V x = {B × R | B ∈ B x } and V n = {B × R | B ∈ B n }. We note that C0 (resp. B0) is the domain in R 3 (resp. R 2 ) defined as C0 \ C∈V x ∪V n i(C) (resp. π (x,y) (C0)).
Let V be the set V x ∪V n ∪V b . Let E be the set of connected components of C ∩C0. We define the incidence relation I : E → V × V such that for e ∈ E, I(e) = (C, C ) if and only if the extremities of e are an out-connection of C ∩ C0 in C and an inconnection of C ∩ C0 in C . By construction, the projections of edges do not cross, and the incidence of four edges at a node is correctly encoded. This yields the following proposition.

Proposition 5.1 Let G be the graph (V, E, I). The curve B ∩ B0 is an embedding of G, seen as a non-directed graph.

We denote by v1, v2, . . . the vertices of G (i.e. the boxes of V = V x ∪ V n ∪ V b ) and by e1, e2, . . . the edges of G, (i.e. the connected components of C ∩ C0). The graph G representing the apparent contour of the torus is shown in the left part of Figure 7.

We show in Section 5.1.1 how to distinguish in and out connections of C ∩ C0 in boxes of V , and, in Section 5.1.2, how to compute the incidence relation I by computing δ-approximations of the connected components of C ∩C0. We consider G to be equipped with an application approx that maps each e ∈ E to its δ-approximation. Furthermore, approx satisfies the three following properties: In the following, we consider the faces of B ∩ B0 in B0 that are the connected components of B0 \ B. It is thus convenient to represent the boundary of B0 explicitly in a graph

G ∂ = (V ∂ , E ∂ , I ∂ ) of which the curve B ∂ = (B ∩ B0) ∪ ∂B0 is an embedding. G ∂ is defined as follows. If V b is empty, then V ∂ = V ∪ {B ∂ × R} where B ∂ is a box reduced to the left bottom corner of B0 and E ∂ = E ∪ {e ∂ } where e ∂ is ∂B0 \ B ∂ . Otherwise (V b is not empty), V ∂ = V and E ∂ = E ∪ E where E is the set of connected components of ∂B0 \ C∈V b i(π (x,y) (C)). Then I ∂ (e)

In and out connections

We characterize in and out connections of C ∩ C0 in a box of V x ∪ V n or in C0 using the direction of the tangent to the curve. For this, let E * (∂B) (resp. W * (∂B)) denote the set (E(∂B) ∪ SE(∂B) ∪ N E(∂B)) ∩ B (resp. (W (∂B) ∪ SW (∂B) ∪ N W (∂B)) ∩ B). We first consider the case where C ∈ V x ∪ V n . Proposition 5.2 Let c be a connection of C ∩ C0 in a box C ∈ V x ∪ V n , and B = π (x,y) (C). Let t = (tx, ty, tz) be a tangent vector of C at c. Then Proof:According to Section 4, the boxes in V x ∪V n do not contain any x or y-critical points on their boundary, so property (i) holds. Claim (ii) (resp. (iii)) rephrases the conditions such that the tangent vector pointing out of the box is oriented to the left (resp. right) thus yielding an out-connection, (resp. in-connection).

Proposition 5.2 is easily adapted for boundary points, that is, when c is a connection of C ∩ C0 in C0. From assumptions (A4), (A10), ∂C0 contains neither x nor y-critical points, and (i) follows. Then one can easily show that properties (ii ) and (iii ), obtained from (ii) and (iii) by swapping E and W , and N and S, hold, since in this case the tangent vector pointing inside the box is relevant.

Proposition 5.2 is used to decide if a connection c implicitly defined by the ordered pair (Ci, Si) is an out or an in connection, as follows. Ci can be contracted with the Krawczyk operator KS i until 0 / ∈ R(Ci) and 0 / ∈ R (Ci), where R and R as in Definition 2.1 are the x and y component of a tangent vector. This process terminates from (i) of Proposition 5.2, and (ii), (iii), (ii ), (iii ) are used to finish the argument. The right part of Figure 7 illustrates this in the case of a node singularity. The cones containing the tangent vector (R(Ci), R (Ci)) are drawn with dashed lines for two of the four boxes isolating the connections in its witness box.

Computing the incidence relation I and δ-approximations of edges

Let C + (resp. C -) be the set of implicit points defining the out (resp. in) connections of C ∩ C0 in boxes of V x ∪ V n and in C0. For each out-connection (Ci, Si) ∈ C + , the process Track(< P, Q >,C0,(Ci, Si),C -∪{(Ci, Si)},δ) defined in Algorithm 2 is performed. This process terminates and returns a δ-approximation (C e i ) me 1 of the connected component e of C ∩ C0 having (Ci, Si) as left-most point and the set C on containing exactly the implicit point (Cj, Sj) defining the right-most point of e. This defines the edge e ∈ E with I(e) = (v, v ), where v (resp. v ) is the vertex of G having (Ci, Si) (resp. (Cj, Sj)) as one of its out (resp. in) connections. We let approx(e, δ) = (C e i ) me i=1 . To ensure that the application approx satisfies the properties (a1), (a2) and (a3) listed above, the tracking of edges is refined as follows. Each time a new approximation approx(e, δ) is computed, (a2) is checked and a δ 2 -approximation of e is computed while it does not hold. This process terminates due to the properties of witness boxes. Then for each e for which approx(e , δ) is already known, (a1) is checked for e, e . While it does not hold, δ 2 -approximations of e and e are computed. This process terminates since e and e are, by construction, non-intersecting smooth curves. Finally (a3) is checked for e, e . Let approx(e, δ) = (C e i ) me i=1 and approx(e , δ) = (C e i ) m e i=1 . Since (a2)

Algorithm 7 CMap construction

Input: A graph G k = (V k , E k , I k ) Output: A CMap (H k , σ k , α k ) encoding the embedding B k ∂ of G k 1: Let V k = {i | vi ∈ V k }, E k = {i | ei ∈ E k }
and (Oi) i∈V k be empty sets. 2: Let H k = ∅, α k be an involution on H k and σ k be a permutation on H k 3: for i ∈ E k do 4:

Let H k = H k ∪ {h + i , h - i }, α k = α k ∪ {(h + i , h - i )} 5: Suppose I k (ei) = (vj, v l ) with j, l ∈ V k and let Oj = Oj ∪{h + i }, O l = O l ∪{h - i } 6: for i ∈ V k do 7:
Order counter-clockwise the connections of Oi around ∂(π (x,y) (vi))

8: Let σ k be the permutation on H k which orbits are Oi, for i ∈ V k . 9: return (H k , σ k , α k ) holds, the ball system can have solutions only in the 4-dimensional boxes constructed from C e 1 , C e 1 , C e me , C e m e as in Proposition 3.1. Then, δ 2 -approximations of e and e are computed while the ball system has a solution in one of the latter boxes. This process terminates since each singularity of B has a strictly positive distance from any point of the boundary of a witness box around it. To avoid x-extreme points, one has to check that boxes C e 1 and C e me do not contain solutions of the system (S5) as defined in Corollary 2.1, and refine the boxes if needed.

Computing the CMap of one connected component

Let B 1 ∂ , . . . , B n ∂ be the connected components of B ∂ . Then G ∂ has n connected com- ponents G 1 = (V 1 , E 1 , I 1 ), . . . , G n = (V n , E n , I n ) such
that the vertices of G i are the cusps, nodes, x-extreme and boundary points of B i ∂ , and the projections of the edges of G i are the connected components of B i ∂ \ i(B0). We suppose that G 1 is the connected component containing the boundary points of B ∂ . We aim here at computing, for each 1

≤ k ≤ n , a CMap (H k , σ k , α k ) representing the embedding B k ∂ of G k = (V k , E k , I k ).
As emphasized in [START_REF] Sergei | Graphs on surfaces and their applications[END_REF][ §1. 3.3], or in [START_REF] Köthe | Xpmaps and topological segmentation-a unified approach to finite topologies in the plane[END_REF], computing the faces of an embedding reduces to ordering counter-clockwise the edges around each vertex. The CMap (H k , σ k , α k ) representing the embedding

B k ∂ of G k = (V k , E k , I k
) encodes this order in the permutation σ k , and each face of B k ∂ is an orbit of ϕ k = (σ k ) -1 • α k which is an ordered sequence of half-edges describing a counter-clockwise walk around it. A vertex of V k corresponds either to a node, a cusp, a boundary point or a x-extreme point of B k . Around a cusp or an x-extreme point p, B k \ p has only two branches; thus there is no need to order them.

Let p be a node of B k and B its associated witness box. Since the connected components of (B k \ p) ∩ B are non intersecting curves linking p to the projections of the connections of C in B × R, a counter-clockwise ordering of the branches of B k \ p around p is given by the counter-clockwise order of the connections of B k in B (see the right part of Figure 7). Now let p be a boundary point of B k ∩ B0. A counter-clockwise ordering of the branches of B k ∂ \ p around p can be directly deduced from the part of the boundary to which p belongs. Consider the left part of Figure 7, and let p be the boundary point of which v1 is a witness box. Since p is on E * (∂B), a counter-clockwise ordering is necessarily (e1, e14, e13).

Algorithm 8 Point-face location on the boundary of a witness box

Input: The CMap (H k , σ k , α k ) of B k ∂ , a box C ∈ C x ∪ C n ,
+ i and h - i , such that h + i is oriented as ei, i.e. if I k (ei) = (vj, v l ), h + i leaves vj and h - i leaves v l .
In other words, h + i represents a walk from left to right along π (x,y) (ei). In addition, half-edges leaving a common vertex are collected. The second loop aims at ordering the half-edges around vertices to define faces.

Faces of a CMap

To define the relative positions of the CMaps encoding the different connected components and thus computing an XPMap, the exterior face of each CMap is first identified. Then, we propose a procedure to identify to which face of B k a point p / ∈ B k belongs.

Exterior face

Point-face location on the boundary of a witness box. As an intermediate step, Algorithm 8 identifies the face containing a point located on the boundary of a witness box. In the algorithm, h * stands for a half-edge h + or h -. In Step 3, the connections of C in C are contracted with the appropriated Krawczyk operator until their projections are pairwise disjoint and do not contain p. This step terminates, since p / ∈ B. Then, it suffices to order counter-clockwise p and the projections of the connections of C in C on ∂(π (x,y) (C)) to conclude. The right (resp. left) part of Figure 8 illustrates this procedure when π (x,y) (C) is the witness box of a node (resp. cusp). The orbit of the node vertex is (h

- 1 , h + 2 , h + 3 , h - 4
) and (c1, p, c2, . . .) is a counterclockwise order on the boundary of the box. The face containing p is given by the orbit of ϕ k containing h - 1 .

Exterior face and leftmost box. We now explain Algorithm 9, that computes the orbit of ϕ k describing the exterior face of B k ∂ . Note that by construction all faces but the exterior face are described by a counter-clockwise cycle of half-edges. The case where k = 1 is directly addressed, since B 1 ∂ contains ∂B0. The exterior face of B 1 ∂ is exactly R 2 \ B0 and is described by the orbit of ϕ 1 containing the half-edges associated

0 0 1 1 0 1 h - 1 v c 1 c 2 h - 2 h + 2 h + 1 h + 4 h - 4 h + 1 h - 1 c 4 c 1 c 3 h + 3 h - 3 h + 2 h - 2 c 2 p p v

Algorithm 9 Exterior face

Input: The CMap (H k , σ k , α k ) of B k ∂ Output: The orbit of ϕ k describing the exterior face of B k ∂ and a leftmost box C enclosing B k ∂ if k > 1, C = ∅ if k = 1. 1: if k = 1 then 2:
Let e be an edge of E ∂ \ E and h + be the half-edge associated with e

3:

Let f be the orbit of ϕ k containing h + and C = ∅ 4: else 5:

Find a leftmost box C enclosing B k ∂ and let p be any point on W * (∂(

Let f be the result of Algorithm 8 with input (H k , σ k , α k ), C and p

7: return f , C
with edges of E ∂ \ E (see Steps 1, 2 and 3 of Algorithm 9). As an illustration, the exterior face of B 1 ∂ on Figure 3 is the orbit (h + 13 , h + 14 ). Suppose now k > 1, and recall that B k ∂ is an embedding of G k = (V k , E k , I k ). Recall also that for B = (x, y), W * (∂B) = (l(x), y). We state the following remark to identify the exterior face of B k ∂ .

Remark and Definition 5.1 Let k > 1.

(i) B k ∂ has at least two cusps or x-critical points, that is V k ∩ V x contains at least two elements.

We call a leftmost box enclosing

B k ∂ a box C = (x, y, R) of V k ∩ V x minimizing l(x) over all boxes of V k ∩ V x .
Let C be one of the leftmost boxes enclosing B k ∂ , and

B = π (x,y) (C). Then (ii) any point of W * (∂B) lies in the exterior face of B k ∂ , (iii) if p ∈ W * (∂B) then p / ∈ B.
Point (i) is a direct consequence of Proposition 2.1. In order to prove point (ii), consider a point p ∈ W * (∂B), and suppose it does not belong to the exterior face of B k ∂ . Hence B k ∂ has necessarily a cusp or x-critical point lying to the left of p. As a consequence, there is a box C = (x , y , R) in V k ∩ V x such that l(x ) < px, where px is the x-coordinate of p, and C is not a leftmost box enclosing B k ∂ . Consider now a point p ∈ W * (∂B) such that p ∈ B. Since C is the witness box, one has necessarily p ∈ B k ∂ and a contradiction follows. Hence (iii) holds. Algorithm 9 uses Remark 5.1 in Step 5. From point (iii), the call to the procedure described in Algorithm 8 terminates and returns the orbit of ϕ k describing the exterior face of B k ∂ .

Algorithm 10 Point-face location in a CMap

Input:

A point p = (xp, yp) ∈ i(B0). The CMap (H k , σ k , α k ) encoding the em- bedding B k ∂ of G k = (V k , E k , I k ). The exterior face f ext of (H k , σ k , α k ). The boundary edges E = E ∂ \ E of G ∂ . Preconditions: (c1) p is not in B k ∂ , (c2)
p is not in the projection of any box of V k and (c3) ∀ej ∈ E k , p is not in the projection of any box of approx(ej, δ). Output: The orbit of ϕ k describing the face of B k ∂ that contains p.

1: Let S = {P (xp, y, z) = Q(xp, y, z) = 0}, B * = (y * , z * ) = ∅ and f * = ∅ 2: for ej ∈ E k \ E do //ej is not a boundary edge 3:
Let X sol be the result of Algorithm 11 with input ej, p and S

4:

if X sol = ∅ then 5:

Let B * j = (y * j , z * j ) be the unique element in X sol 6: if l(y * j ) < l(y * ) then Let f * be the orbit of ϕ k containing h - j and let

if B * = ∅ then B * = B *
B * = B * j 10: for C ∈ (V x ∪ V n ) ∩ V k where C = (x, y, z) do 11: if xp ∈ i(x) and yp < l(y) < l(y * ) then // when B * = ∅ we set l(y * ) = ∞ 12:
Let y * = y 13:

Let f * be the result of Algorithm 8 with input (H k , σ k , α k ), C and the point (xp, l(y))

14: if f * = ∅ then 15: if k = 1 then 16:
Let e ∈ E be the boundary edge above p and f * be the orbit of ϕ k containing h -

17:

else Let f * be the exterior face f ext of (H k , σ k , α k ) given by Algorithm 9

18: return f *

Point-face location in a CMap

We now describe Algorithm 10 for the point-face location in a CMap. The preconditions (c1), (c2), (c3) ensure its termination. In order to show the correctness of Algorithm 10, we consider the segment B seg = xp × [yp, u(y0)] and the set C seg = B seg × R, where B0 = (x0, y0). The idea of the proposed approach is to find the closest intersection

q of B k ∂ with B seg . If q belongs to a smooth component π (x,y) (ej) of B k ∂ (i.e. a connected component of B k ∂ ∩ B0
) the face is given by the orbit containing the half-edge h

- j . If q is in the witness box C ∈ (V x ∪ V n ) ∩ V k ,

Algorithm 11 Edge intersection

Input: An edge e ∈ E k , a point p = (xp, yp), the system S = {P (xp, y, z) = Q(xp, y, z) = 0}. Preconditions: (c1) p is not in B k ∂ and (c3) ∀ej ∈ E k , p is not in the projection of any box of approx(ej, δ). Output: A set X sol containing the unique intersection of e with C seg if it exists, X sol = ∅ otherwise.

1: Let (Ci) m i=1 = approx(e, δ) and ∀1 ≤ i ≤ m, Ci = (xi, yi, zi) and Xi = (yi, zi) 2: for 1 ≤ i ≤ m s.t. xp ∈ xi and l(yi) > yp do //Ci is located above p 3: Let (X sol , X ind ) = IsolateSols(S,1.1Xi,min( 0.1w (x i ) 2 , 0.1w (y i ) 2
))

4:

if X sol = ∅ then 5: break 6: return X sol Consider now the case where C seg does not intersect boxes of V k but intersects some edges. The intersections of C seg with the smooth and x-monotone components ej ∈ E k \ E are computed in the loop beginning in Step 2. Algorithm 11 computes, if it exists, the unique intersection q = (xp, yq, zq) of ej with C seg and returns a nonempty set containing B * j ⊂ C seg with (yq, zq) ∈ i(B * j ). B * j is made disjoint in the y-coordinate with the current intersection B * in Step 8, and B * j is updated in Step 9 to contain the intersection of C with C seg which projection is the closest above p. Due to the orientation of half-edges, the face below an edge ej is described by the orbit of ϕ k containing h - j . Thus f * as updated in Step 9 is the face containing p. Before describing the last case (i.e. when C seg intersects some boxes of (V x ∪ V c ∪ V b ) ∩ V k ) we remark that C seg can intersect only one boundary box in V b ∩ V k corresponding to the implicit point (B, S4) where S4 is defined as in Corollary 2.1 and B = (x, z) is such that xp ∈ x. Letting ej ∈ E k be such that (B, S4) is a connection of ej and approx(ej, δ) = (Ci) m j i=1 , it follows from property (a2) of approx that the boundary point implicitly defined by (B, S4) is either in C1 or in Cm j . If the intersection (xp, u(y0)) is in C1 or in Cm j , the appropriated face has been determined above. Otherwise, (B, S4) can be refined with KS 4 until (xp, u(y0)) / ∈ B, and this case is handled in Steps 14 to 17 of the algorithm.

Suppose now that C seg intersects some boxes of (V x ∪ V n ) ∩ V k . We can assume that there is a box C ∈ (V x ∪ V n ) ∩ V k with l(y) < l(y * ), otherwise the output of the algorithm is determined by the loop of Step 2 or Steps 14 to 17. Such a box does not contain p from precondition (c2). The loop beginning in Step 10 finds the box (x, y, z) ∈ (V x ∪ V n ) ∩ V k intersecting C seg that minimizes l(y), thus p lies in the same face than the point (xp, l(y)). Note that (xp, l(y)) / ∈ B k ∂ , otherwise (xp, l(y)) is the projection of a connection of C in C, and B * found in the loop beginning at Step 2 necessarily satisfy l(y * ) ≤ l(y). As a consequence, the input arguments given to Algorithm 8 in Step 13 satisfy its preconditions.

Proof:[Proof of termination of Algorithm 10] It is shown below that if its input arguments satisfy preconditions (c1) and (c3), Algorithm 11 terminates. We have already shown that Algorithm 8 terminates due to condition (c1). It only remains to prove the termination of the while loop in Step 8. The points of the curve represented by the boxes B * and B * k belong to two different edges which are disjoint in projection on the (x, y)-plane, so after a finite number of contractions the y-coordinates will be disjoint intervals.

Proof:[Proof of correctness and termination of Algorithm 11] Let approx(e, δ) = (Ci) m i=1 and Ci = (xi, yi, zi) and Xi = (yi, zi). From properties (a2) and (a3) of approx, e = C ∩ m i=1 Ci contains no x-extreme point of C and is a smooth xmonotone curve. Hence its intersections with C seg are regular solutions of S. From property (a1) of approx, for any 1 ≤ i ≤ m, 1.1Ci does not intersect any box of any other approximation. Hence for any 1 ≤ i ≤ m, 1.1Xi contains at most one solution of S. If 1.1Xi contains a solution of S it is non-singular as a point of e . As a first consequence of this, the call IsolateSols(S,1.1Xi,min( 0.1w (x i ) 2 , 0.1w (y i )

Constructing the XPMap

The relation cont that assigns the exterior face of a CMap to a non-exterior face of another CMap (or to the infinite face) can be computed iteratively thanks to the following remark, that can be proved similarly to Remark 5.1. Let (H k , σ k , α k ) be the result of Algorithm 7 with input G k

4:

Let f k , C k be the result of Algorithm 9 with input (H k , σ k , α k )

5: Let H = H 1 , σ = σ 1 , α = α 1 , ext = {(H 1 , f1)} and cont = {(f0, f1)} 6: Suppose C k = (x k , y k , z k ) for 1 < k ≤ n 7: Re-index G k , (H k , σ k , α k ), f k , C k wrt increasing l(x k ) for 1 < k ≤ n 8: for 1 < k ≤ n do 9:
Let p be any point on W * (∂(π (x,y) (C k )))

10:

Let i = 0 and f = f k 11:

while f = f k-i do 12:

Let i = i + 1

13:

Let f be the result of Algorithm 10 with input p, (H k-i , σ k-i , α k-i ) and f k-i 14:

Let H = H ∪ H k , σ = σ ∪ σ k , α = α ∪ α k , ext = ext ∪ {(H k , f k )}, cont = cont ∪ {f, f k } 15: return (H, H 0 , σ, α, ext, cont)
box of V 1 ∂ ∪ . . . ∪ V k-1 ∂ . Then, ∀ei ∈ E 1 ∂ ∪ . . . ∪ E k-1 ∂
, p is not in any box of approx(ei, δ), as a consequence of the property (ii) of the application approx.

we have several points that project on the same point and step 6 will not terminate, since it cannot separate them. Finally, if (A3), (A4), (A5), (A7), (A8) or (A10) is not satisfied, then some of the steps using IsolateSols will not terminate, since they will be applied to a system with singular solutions.

Our implementation

The two cornerstones of our implementation are the procedure IsolateSols specified in Algorithm 1 and the procedure Track specified in Algorithm 3 and described in Appendix A.

Our implementation of the procedure IsolateSols with subdivision, centeredform interval evaluation and Krawczyk operator is exhaustively described in [START_REF] Imbach | A Subdivision Solver for Systems of Large Dense Polynomials[END_REF]. It uses adaptive multi-precision arithmetic. It is available in the mathematical software SageMath4 as the package subdivision solver5 . Solving a system of polynomial equations in an unbounded box reduces to solving several transformed systems in bounded boxes (see. [START_REF] Neumaier | Interval methods for systems of equations[END_REF]Section 5.6] or [START_REF] Stahl | Interval Methods for Bounding the Range of Polynomials and Solving Systems of Nonlinear Equations[END_REF]Section 5.10]). In particular, we don't need to know in advance the bounded interval z0 for the z-component. Solving a system on a box of the shape x×y×R reduces to solving two systems on the bounded boxes x×y×[-1, 1] and hence requires two calls to the solver provided by subdivision solver. Notice that a δ-approximation of C ∩ C0 yields bounds for the z-coordinates of C ∩ C0. Once Step ( 3) is performed, we use these bounds to avoid the second call in each algorithm involving solving a system in a box with unbounded z-component.

The procedure Track is described in detail in Appendix A. In this description, we assume that computations are carried out with arbitrary precision. We implemented this procedure in Python within SageMath only for machine precision; our tracker stops when the width of the boxes of the enclosure reaches the machine precision. However, we never encountered this case in the experiments reported below. In Steps (3) and [START_REF] Delanoue | A numerical approach to compute the topology of the apparent contour of a smooth mapping from R 2 to R 2[END_REF], the initial value δ = 1 is used as input for Track. In step [START_REF] Delanoue | A numerical approach to compute the topology of the apparent contour of a smooth mapping from R 2 to R 2[END_REF], the conditions (a2) and (a3) are enforced during the tracking process. The condition (a1) is checked a posteriori on the δ-approximations. Each time the latter condition does not hold, i.e. each time two boxes C e i and C e j of δ-approximations of edges e and e have a non-empty intersection, the δ-approximations of e and e are refined within C e i and C e j while enforcing (a2) and (a3) until (a1) holds. The other algorithms used by our approach have been implemented in Python within SageMath.

Experimental data

The surfaces are defined by random dense polynomials P in Z[x, y, z] with odd total degrees from 5 to 15 and integer coefficients chosen uniformly in -2 8 , 2 8 . We isolate the singularities and compute the topology of the apparent contour of the surface defined by P with two state-of-the-art methods and the approach described here. For each degree, five instances are considered, and we give averages of sequential times in seconds and standard deviations for each method.

State-of-the-art methods

In [START_REF] Imbach | A certified numerical algorithm for the topology of resultant and discriminant curves[END_REF], the singularities of an apparent contour are characterized as the real solutions of a system of two polynomials that are coefficients of the sub-resultant chain. This system can be solved with subdivision solver that has been designed as a solver dedicated to large dense polynomials. We use this approach to isolate the singularities in

B0 = [-1, 1] × [-1, 1].
The package Isotop6 for Maple computes the topology of a plane curve in R 2 , see [START_REF] Cheng | On the topology of real algebraic plane curves[END_REF]. Here, Isotop is used to compute the topology of the resultant of P and Pz with respect to z.

Results

We report results obtained with these three methods for isolating the singularities and computing the topology of the apparent contour of algebraic surfaces. We first give details for one of the surfaces of degree 13 we tested before giving synthetic data for surfaces of degrees from 5 to 15. Running times given below are sequential running times in seconds on an Intel(R) Core(TM) i5-3317U CPU @ 1.70GHz machine with Linux.

Details for a surface of degree 13. We detail the computation of the topology of the apparent contour of one of the surfaces of degree 13 we considered. Figure 9 displays its apparent contour and Table 1 details each step of the computation. Columns t give the running times. For Step [START_REF] Dedieu | Points fixes, zéros et la méthode de Newton[END_REF], the column i gives the maximum number of times a box isolating a singularity or an x-critical point is contracted before it can be certified as a witness box. The symbol * signifies that multi-precision was required to address this step. Table 2. The first group of columns reports on the isolation of the singularities of the apparent contour in B0 with the approach using the sub-resultant system [START_REF] Imbach | A certified numerical algorithm for the topology of resultant and discriminant curves[END_REF] and Steps (1) to (4) of our approach. The former approach suffers from the size in term of degree, number of monomials and bit-size of the coefficients of the equations of the subresultant system. For instance, the first equation of this system for the first polynomial of degree 9 we tested has degree 57, more than 1700 monomials and its coefficient bitsize is more than 130. Furthermore, the running times have a high standard deviation, and machine precision was not sufficient to carry out the computations for examples of high degree.

The group of columns "Computing topology" refers to the computation of the topology of the apparent contour, including the isolation of singularities. The column Isotop reports the running times of Isotop applied to the resultant of P and Pz with respect to z. As expected, the size of the resultant polynomial excludes this approach for surfaces of high degree. We tried Isotop for d up to 8. For d = 8, the running time was 1924 seconds.

In contrast, our method does not consider any resultant or sub-resultant polynomial; hence we deal with systems having almost the same degree and bit size as the input. On the other hand, we have up to four variables instead of two. The machine precision was sufficient to isolate the singularities of all examples, and only one example, the one detailed above, required use of more precision for the computation of the topology. Singularities are characterized here as the solutions of the ball system that involves four equations in four unknowns. The results in [START_REF] Imbach | Numeric and certified isolation of the singularities of the projection of a smooth space curve[END_REF] already showed the advantage of our subdivision solver over the symbolic approach, but they also showed the limitation of our approach when solving a four dimensional system on a large domain instead of a two dimensional one even with higher degree and bitsize. Here, the δ-approximation of C enables us to filter the domain where the ball system is solved, and the results of Table 2 show the efficiency of this strategy. For the computation of the singularities in Steps (1) to (4): Step (2) solves systems in 2 dimensions and is the least time consuming; Step (1) solves systems in 3 dimensions in the large domain B0 × R which is more time consuming; Step (4) solves systems in 4 dimensions in m + m∩ small boxes. For the particular case of degree 13 detailed above, m + m∩ = 3490. Comparing with the results of [START_REF] Imbach | Numeric and certified isolation of the singularities of the projection of a smooth space curve[END_REF] that performed isolation on a unique large domain, the benefit the δ-approximation to reduce the solution domain is dramatic.

In Step (5), the y-critical points are found by solving a system in 3 dimensions in the m boxes of the δ-approximation of C. Comparing the running times to address this step with the ones required for Step (1) illustrates again how using the δ-approximation of C speeds up the computation.

In Step (6), witness boxes for x-extreme points and singularities are computed. It is adressed by combining Algorithm 6, that contracts a box containing a node, and Algorithm 5, that checks if the resulting box is a witness box, as well as their equivalents for x-extreme points and cusps. In almost all cases we tried, Algorithm 5 succeeded with one iteration of Algorithm 6. The time for Step [START_REF] Dedieu | Points fixes, zéros et la méthode de Newton[END_REF] approximately corresponds to |C x |+|D n | times half the time for Step (2), except when i is high (see details given in Appendix B). It appeared in our experiments that most of the time spent in the procedure associated with Algorithm 5 is used to construct the systems defining the boundary points (that requires partial substitutions in large polynomials) rather than to solve them. For the example of degree 13 detailed above, approximately 95% of the time required by the latter procedure is spent defining these systems. For this example, multi-precision arithmetic was required to carry out Step [START_REF] Dedieu | Points fixes, zéros et la méthode de Newton[END_REF]. Determining a witness box for an x-extreme point for this particular example required 5 iterations in the equivalent of Algorithm 6 for x-extreme points before the equivalent of Algorithm 5 for x-extreme points certifies it to witness the x-extreme point. It should be considered as a proof of robustness of our implementation more than as a drawback of the approach. Similarly to Step (2), Algorithm 5 involves isolating the intersections of C on the boundaries of a cylinder. Step (2) requires two calls to the numerical 0-dimensional solver whereas in Algorithm 5 we use the bounds for the z-component given by the δ-approximation of C ∩ C0. Hence Algorithm 5 requires roughly half the time of Step (2) to execute.

Among Steps (1) to (8), Step [START_REF] Delanoue | A numerical approach to compute the topology of the apparent contour of a smooth mapping from R 2 to R 2[END_REF] is the most time consuming. It consists of two substeps. First, δ-approximations of connected components of C ∩ C0 are computed while ensuring that conditions (a2) and (a3) hold. This results in approximations with smaller boxes than in Step (3); see the details in Appendix B. Then, the condition (a1) is checked for each pair of δ-approximation. In most examples we tested, (a1) was satisfied without refining the approximations. For the example of degree 13 detailed above, computing the approximations required 18.8s and checking (a1) required 4.6s.

Most of the apparent contours we computed have only one connected component; this explains the running times for Step [START_REF] Demazure | Bifurcations and catastrophes: geometry of solutions to nonlinear problems[END_REF]. See also column n in the table in Appendix B.

Conclusion

In this article, we provided an algorithm computing the XPMap of the projection of a space curve. The space curve is defined by the intersection of two surfaces represented implicitly by convergent interval functions. Compared to the state of the art, our method is the first to solve this problem reliably for functions that are not polynomial. Moreover, in the polynomial case, our experiments show that our method handles polynomials of degree up to 15 that are not yet reachable by sub-resultant approaches. As future work, we plan to generalize our method to the projections of curves from R n to R 2 . In this case, we would need to change the systems modelling the node and cusp singularities, but the tracking of the curve extends easily. We would also like to address the problem of computing the projection of 2-manifolds from R 4 to R 3 . Several certified numerical path-tracking algorithms can be found in the literature [START_REF] Faudot | A new robust algorithm to trace curves[END_REF][START_REF] Kearfott | An interval step control for continuation methods[END_REF][START_REF] Martin | Certified parallelotope continuation for one-manifolds[END_REF]. They approximate a smooth connected curve in a bounded box with a sequence of sets (Ci) m k i=1 . The sets Ci are in general boxes as in [START_REF] Faudot | A new robust algorithm to trace curves[END_REF][START_REF] Kearfott | An interval step control for continuation methods[END_REF]. It has recently been proposed to use parallelotopes instead of boxes [START_REF] Martin | Certified parallelotope continuation for one-manifolds[END_REF]; a parallelotope C is the image by an affine transformation of a box C. Aligning parallelotopes along the tangent to the curve yields a more efficient approximation.

For the sake of generality, we consider here approximating sets Ci that are either boxes or parallelotopes. If C is such a set, h(C) denotes its box hull, that is the smallest box containing C. The operators ∂, i() are directly generalized for parallelotopes. Let C be a parallelotope and C, f be the box and the affine map such that C = f (C). The width w (C) of C is the width of C and if ∈ R, the -inflation C is defined as f ( C).

Section A.1 introduces Algorithm 13 using -inflation to decide if an implicit point lies in a set C. Section A.2 recalls the parallelotope path-tracking algorithm of [START_REF] Martin | Certified parallelotope continuation for one-manifolds[END_REF], and Section A.3 shows how it is adapted to meet the specifications of Algorithm 2.

A.1 Deciding if an implicit point lies in a set of the approximation

Let (C, S) define implicitly a point p. To check if p belongs to a set C, a naive approach consists of contracting C on p until C ⊂ i(C) or C ∩ C = ∅. If p lies on ∂C, this leads to a non-terminating process. To tackle this pitfall, we use the -inflation approach described in Algorithm 13, that takes as input C and a real number > 1 and returns true if p ∈ C and false if p / ∈ C. Notice that when p ∈ C \ C, it may return either true or false.

Proof:[Proof of termination and correctness of Algorithm 13] Assume > 1 and let C be C. After applying Step 3 of Algorithm 13 a finite number of times, C is strictly included in a ball of diameter -1 2 w (C). Hence C can not intersect both ∂C and ∂C; thus either C ⊂ i(C ) or C ∩ C = ∅, and the while loop of Algorithm 13 terminates, so does Algorithm 13.

B Details for Results

Isolating singularities

Step (

Figure 2 :

 2 Figure 2: Left: a δ-approximation of C in B0 × R, and a δ-approximation of B in B0.Right: an embedding of the graph G for the apparent contour of a torus. Black circles represent the points associated with the vertices of G, and thick curves the curves associated with its edges. Thin lines represent the boundary of B0.

Figure 3 :

 3 Figure 3: An XPMap encoding the topology of the apparent contour of the torus in B0. Each edge is represented by a pair of half-edges, which is a cycle of the permutation α. Curved arrows around vertices represent cycles of the permutation ϕ.

Definition

Definition 2 . 1 A

 21 point p ∈ C is x-critical if the x component of the tangent of C at p vanishes. The x-critical points of C are the solutions of the system P

(

  A3) P (x * , y, z) = Q(x * , y, z) = 0 has finitely many regular solutions when x * = l(x0) or u(x0), y ∈ y0 and z ∈ R. P (x, y * , z) = Q(x, y * , z) = 0 has finitely many regular solutions when y * = l(y0) or u(y0), x ∈ x0 and z ∈ R. (A4) P = Q = R = 0 has finitely many regular solutions in C0, and no solution in ∂C0.

Proposition 2 . 1

 21 Assuming (A1), the connected components of C ∩ C0 are smooth one dimensional manifolds possibly with boundary. In addition, assuming (A2) and (A3), any connected component C k of C ∩ C0 satisfies at least one of the following statements: (a) C k has exactly two boundary points, (b) C k has at least two x-critical points.

10 :

 10 remove C on ∩ C x from L x 11: Let m = m k 12: return (Ci) m i=1 Algorithm 3 first computes δ-approximations for every connected component of C that is diffeomorphic to [0, 1]. This is addressed by calling the procedure Track with boundary points as initial points. The boundary points are implicitly defined by elements of C b . Thus letting C k be a connected component containing a boundary point defined by (C, S) ∈ C b , the call to the procedure Track in Step 5 terminates, returns a δ-approximation of C k and identifies the two extremities and the x-critical points of C k . When entering the while loop in Step 3 for the first time, the list L b contains the two extremities of each connected component of C ∩ C0 that is diffeomorphic to [0, 1]. Each time the instruction in Step 6 is performed, the two extremities of the connected component that has been tracked are removed from L b , and the size of L b decreases.

Figure 4 :

 4 Figure4: Singularities of the apparent contour B of the torus. For node and cusp singularities of B, their preimages on the space curve C as well as corresponding centers c and radii r2 = r 2 for the ball system are represented.

Figure 5 :

 5 Figure 5: Some boxes and their projections containing singularities of B. Cusps singularities are in boxes Bi, nodes in boxes Bi ∩ Bj.

  While refining the box D, Algorithm 4 thus terminates for a box D * = (x * , y * , c * , r * ) ⊂ D, containing p such that either 0 / ∈ r * and π (x,y) (p) is a node, or K (P,Pz ,Pzz ) ((x * , y * , c * )) ⊂ i((x * , y * , c * )),

  define the boundary points of C as the set C ∩ (∂B0 × R), and the boundary points of B as the set B ∩ ∂B0. With Assumptions (A3) -(A5), the boundary points of C and B are in one-to-one correspondence, are neither x nor y-critical points of C nor of B, and are not the corners of B0. The isolating boxes of boundary points of C can thus be refined until their projections onto ∂B0 are disjoint.

Figure 6 :

 6 Figure 6: Left (resp. middle): the decomposition of B (resp. ∂B). Right: a witness box B n for a node.

Algorithm 6

 6 computes a witness box for a node. The idea is that if the curve C does not contain an x or y-critical point above an isolating box B of the node, B is a Algorithm 5 IsWitnessNodeBox(B) Input: A box B containing a unique node of B, such that B × R does not contain any x or y-critical point of C. Output: False when the box is not witness. True when B intersects exactly twice E(∂B), twice W (∂B) and does not intersect N W (∂B) ∪ N (∂B) ∪ N E(∂B) ∪ SW (∂B) ∪ S(∂B) ∪ SE(∂B).

Algorithm 6 1 : 2 : 2 5: repeat 6 : 7 :

 612267 Let B = (xB, yB) be a 2-dimensional box in the (x, y)-plane and let a = w (xB) and b = w (yB). Assuming a < b, we decompose B into three boxes and ∂B into eight closed segments as illustrated in the left and middle parts of Figure 6, and we name the boxes and segments with respect to their cardinal directions around the center of B. In particular, M (B) is (xB, m(yB) + [-a/2, a/2]), and E(∂B) and W (∂B) are the respective segments (u(xB), m(yB) + [-b/4, b/4]) and (l(xB), m(yB) + [-b/4, b/4]).Proposition 4.1 Let B be a box containing a unique node of B such that B × R does not contain any x or y-critical point of C. Then (i) Algorithm 5 terminates, (ii) if Algorithm 5 returns True then B is a witness box for the node, (iii) if B intersects exactly twice E(∂B), twice W (∂B) and does not intersect N W (∂B) ∪ N (∂B) ∪ N E(∂B) ∪ SW (∂B) ∪ S(∂B) ∪ SE(∂B), then Algorithm 5 returns True.Proof: Since there is no x-critical point of C in B×R, the system P (x = l(xB), y, z) = Q(x = l(xB), y, z) = 0 (resp. P (x = u(xB), y, z) = Q(x = u(xB), y, z) = 0) has only regular solutions on yB × R. In addition, since w (yB)/4 > 0, the procedure IsolateSols called in Step 2 (resp.Step 3) of Algorithm 5 terminates even if solutions lie on ∂yB × R. Suppose now the sets of undetermined boxes X ind 1 , X ind 2 are empty. Then P = Q = 0 has no solution above the corners of B, and in particular the systems P (x, y = l(yB), z) = Q(x, y = l(yB), z) = 0 and P (x, y = u(yB), z) = Q(x, y = Witness box for a node Input: A box D containing a unique solution of the ball system that projects into a node σ and such that π (x,y) (D) × R does not contain any x or y-critical point of C. Output: A witness box for σ. Let D = D and B = π (x,y) (D) Let n = 0, an = w (B), bn = w (B) and B n = B 3: repeat 4: n = n + 1, an = an-1/4 and bn = bn-1/D = K b (D ), B = π (x,y) (D ) B n = m(B ) + ([-an/2, an/2], [-bn/2, bn/2]) 8: until B ⊆ M (B n ) 9: until B n ⊆ B and IsWitnessNodeBox(B n ) 10: return B n u(yB), z) = 0 have no solution on ∂xB × R. Hence the procedures IsolateSols called in Step 6 and 7 of Algorithm 5 terminate if and only if the two latter systems have only regular solutions on xB × R, and this is the case since C has no y-critical point over B; (i) follows.

Proposition 4 . 2

 42 Algorithm 6 correctly computes a witness box.Proof: For n fixed, the boxes M (B n ) and B have the same center, and the width of M (B n ) is an. During the repeat loop of line 5, the width of B is strictly decreasing, and the condition B ⊆ M (B n ) of Line 8 will be true after a finite number of loops.If B n σ = σ + [-an, an] × [-2bn, 2bn], this is a strictly decreasing sequence with respect to inclusion, and for

Definition 5 . 1

 51 Recall that C b is the set of implicit points for the boundary points. We define V b as {C | (C, S) ∈ C b }. Since x-critical points of C ∩ C0 are isolated in boxes of V x , the connected components of C ∩ C0 are diffeomorphic to [0, 1] and x-monotone. Their endpoints are C ∩ C∈V x ∪V n ∪V b ∂C 3 . Hence a point of the latter set is either the left-most or the right-most point of a connected component of C ∩ C0, and we state the following definition. Let C be a box of V x ∪ V n or the box C0. We call a point c of C ∩ ∂C a connection of C ∩ C0 in C. If c is the left-most point of a connected component of C ∩ C0, we say that c is an out-connection. Otherwise we call c an in-connection.

Figure 7 :

 7 Figure 7: Left: The directed graph G ∂ for the apparent contour of the torus. Its edges are represented by thin arrows linking out-connections to in-connections of C ∩ C0 in boxes of V . Thick lines represent the curve B in witness boxes. Dashed arrows represent the clockwise walk on ∂B0. Right: A witness box B for a node, and, in bold, the four boxes isolating the connections of B in B.

( a1 )

 a1 If e, e are two edges of G, then approx(e, δ) = (C e i ) me i=1 and approx(e , δ) = (C e i ) m e i=1 satisfy 1.1C e i ∩ C e i = ∅ for all 1 ≤ i ≤ me, for all 1 ≤ i ≤ m e . (a2) If e is an edge of G, C a box of V and approx(e, δ) = (C e i ) me i=1 , then π (x,y) (C e i ) ∩ π (x,y) (C) = ∅ if and only if i = 1 or i = me and e has a connection in C. (a3) If e is an edge of G and approx(e, δ) = (C e i ) me i=1 then π (x,y) (C e 1 ), π (x,y) (C e m ) contain no node, cusp or x-extreme point of B.

  is defined as I(e) if e ∈ E. Otherwise e ∈ E ∂ \E and I ∂ (e) = (v, v ) if and only if the extremities of e are extremities of the segments π (x,y) (v) and π (x,y) (v ), and the walk from v to v around i(B0) is a clockwise walk. It is clear that the curve B ∂ is an embedding of G ∂ . The left part of Figure 7 shows, in dashed arrows, the edges of E ∂ \ E.

  (i) tx = 0 and ty = 0;(ii) c is an out-connection of C ∩ C0 in C if and only if π (x,y) (c) ∈ E * (∂B) or π (x,y) (c) ∈ N (∂B)and txty > 0 or π (x,y) (c) ∈ S(∂B) and txty < 0;(iii) c is an in-connection of C ∩ C0 in C if and only if π (x,y) (c) ∈ W * (∂B) or π (x,y) (c) ∈ N (∂B)and txty < 0 or π (x,y) (c) ∈ S(∂B) and txty > 0.

1 : 3 : 4 : 5 :

 1345 and a point p ∈ ∂B where B = π (x,y) (C). Preconditions: p is not in B Output: The orbit of ϕ k describing the face of B k ∂ to which p belongs. Let v = (. . . , h * i , . . .) be the orbit of σ k corresponding to C 2: Let {. . . , ci = (Ci, Si), . . .} be the connections of C in C associated to v Contract each ci with Ci = KS i (Ci) until their projections are disjoint and p ∈ π (x,y) (Ci) Order counter-clockwise connections ci and p on ∂B and let (c1, p, c2, . . .) be such an ordering Let f be the orbit of ϕ k containing h * 1 6: return f Algorithm 7 performs the construction of the CMap of a connected component B k ∂ of B ∂ . In a first loop, pairs of half-edges are created from edges: to each edge ei ∈ E k are associated two half-edges h

Figure 8 :

 8 Figure 8: Finding the face containing a point p on the boundary of the witness box of a cusp (left) and of a node (right).

8 :

 8 while y * j ∩ y * = ∅ do let B * j = KS(B * j ) and B * = KS(B * ) 9:

  we use Algorithm 8 to determine the face containing the intersection of B seg with S(∂π (x,y) (C)). Let us now give further details. Proof:[Proof of correctness of Algorithm 10] Suppose first that C seg does not intersect any box of approx(ei, δ) where ei ∈ E k nor any box of V k . Then B seg does not intersect B k ∂ \ ∂B0. If k = 1, i.e. B k ∂ is the component of B ∂ containing ∂B0, we let p be the point (xp, u(y0)) ∈ ∂B0. Then p lies on a boundary edge ei of E = E ∂ \ E and the face of B k ∂ containing p is described by the orbit of ϕ k containing h - i . If k > 1, then p belongs to the exterior face of B k ∂ . The algorithm handles these cases in Steps 14 to 17.

2 )

 2 ) in Step 3 of Algorithm 11 terminates, and since any δ-approximation contains a finite number of boxes, Algorithm 11 terminates for any input. Recall now that e ⊂ e . Since both e and e are smooth curves, a point of e belongs to the same connected component of C \ {c ∈ C|π (x,y) (c) is a singularity or a x-extreme point of B ∂ }. Hence for any 1 ≤ i ≤ m, a solution of S in 1.1Xi is a point of the edge represented by e. Reciprocally since e ⊂ m i=1 Ci (approx(e, δ) is a δ-approximation of e), there exists 1 ≤ i ≤ m such that e ∩ C seg is a solution of S in Xi satisfying xp ∈ xi and l(yi) > yp from preconditions (c1) and (c3) and X sol contains the corresponding solution. The correctness of Algorithm 11 follows.

5. 4

 4 Embedding G ∂ In this part we show how to compute an XPMap (H, H 0 , σ, α, ext, cont) representing the topology of B ∂ , as defined in Section 1. It is shown above how to compute a CMap for each connected component B k ∂ of B ∂ , how to identify its exterior face and how to perform point-face location in B k ∂ . Thus computing the XPMap reduces to computing the relation cont that assigns each exterior face to a non-exterior face of another CMap or to the infinite face. We describe in Section 5.4.1 the construction of the relation cont. Point-face location in an XPMap is discussed in Section 5.4.2.

Remark 5 . 1 1 :

 511 Let B k ∂ and B k ∂ with k, k > 1 be two components of B ∂ , let C = (x, y, R) be a leftmost box enclosing B k ∂ and let C = (x , y , R) be a leftmost box enclosing B k ∂ . If l(x) ≤ l(x ), then B k lies in the exterior face of B k . Suppose B 1 ∂ , . . . , B n ∂ are indexed with respect to increasing l(x k ) whereC k = (x k , y k , z k ) is a leftmost box enclosing B k ∂ . As a consequence of Remark 5.1, given k, k , k such that 1 ≤ k < k < k ≤ n , the exterior face of B k ∂ is not contained in any non-exterior face of B k ∂ and is contained in at least one non-exterior face of B k ∂ (since it is at least contained in a non-exterior face of B 1 ∂ ). As a second consequence, if the exterior face f k of B k ∂ is contained in two non-exterior faces f , f of respective components B k ∂ , B k ∂ with 1 ≤ k < k , then the exterior face f k of B k ∂ is contained in f and f k is contained in the non-exterior face f of B k ∂ ∪ B k∂ . Algorithm 12 uses these arguments in an iterative process to construct the XPMap (H, H 0 , σ, α, ext, cont) representing the topology of B ∂ . Recall that (H 0 , σ 0 , α 0 ) is an empty CMap representing the infinite face R 2 \B0. We denote by f0 the empty orbit of ϕ 0 . The loop beginning in Step 2 computes CMaps (H k , σ k , α k ) for B k for 1 ≤ k ≤ n , Algorithm 12 XPMap construction Input: The graph G ∂ and its connected components G k = (V k , E k , I k ) for 1 ≤ k ≤ n Output: An XPMap (H, H 0 , σ, α, ext, cont) representing the topology of B ∂ Let H 0 = ∅ and f0 represent the infinite face 2: for 1 ≤ k ≤ n do 3:

1 ∂ 1 ∂ 1 ∂..Remark 5 . 2

 11152 together with their exterior faces f k and the leftmost box enclosing B k . Since B is the connected component of B ∂ containing ∂B0, its exterior face is necessarily contained in the infinite face. In Step 5, an XPMap (H, H 0 , σ, α, ext, cont) encoding the embedding of B 1 is initialized with H = H 1 , σ = σ 1 , α = α 1 , ext = {(H 1 , f1)} and cont = {(f0, f1)}. In Step 7, leftmost boxes C k are sorted with respect to increasing l(x k ) and objects are re-indexed to apply the arguments detailed above. Let k > 1, and consider the k-th iteration of the for loop beginning in Step 8. Suppose the relation cont has been properly constructed for B 1 ∂ , . . . , B k-, i.e. (H, H 0 , σ, α, ext, cont) encodes the faces of B <k ∂ = B 1 ∂ ∪ . . . ∪ B k-One has to find the non-exterior face of B <k ∂ containing the exterior face f k of B k ∂ . One first finds a point p in f k by taking any point on W * (∂(π (x,y) (C k ))) as justified in Remark 5.1. Due to the property of a witness box, p lies in the same face of B <k ∂ as B k ∂ . The while loop in Step 11 finds this face by calling Algorithm 10 with input p, (H k-i , σ k-i , α k-i ) and f k-i for increasing i. Remark 5.2 below shows that this input meets the preconditions of Algorithm 10 for any i and thus ensures the termination of these calls. An argument stated above shows that the first non-exterior face f obtained is the inner face of B <k ∂ containing f k , i.e. any other non-exterior face f of B k ∂ containing f k contains f . As a consequence, after Step 14, the relation cont is properly constructed for B 1 ∂ , . . . , B k ∂ and (H, H 0 , σ, α, ext, cont) encodes the faces of B <k+1 ∂ Let C be a leftmost box enclosing B k ∂ , and p a point on W * (∂(π (x,y) (C))). From point (iii) of Remark 5.1, p / ∈ B. Since witness boxes of V are pairwise disjoint, p is not in any

Figure 9 :

 9 Figure 9: Left: the apparent contour of an algebraic surface of degree 13. Center (resp. Right): a detailed view of the same curve in the dashed box of the leftmost (resp. central) part. Cusps and nodes lying outside zooming areas are marked with squares.

Algorithm 13

 13 isSolInSet((C, S),C, ) Input: An implicit point (C, S) defining a point p, a set C and a real number > 1. Output: true if p ∈ C, false if p / ∈ C. 1: C = C 2: while not ( (C ⊂ i(C )) or (C ∩ C = ∅) ) do 3: C = KS(C) 4: if C ⊂ i(C ) then return true 5: return false A Implementation of Algorithm 2

Remark A. 1

 1 Let (C, S) be a point of C, and C and > 1 be such that both C ∩ C and C ∩ C are diffeomorphic to [0, 1]. If isSolInSet((C, S),C, ) returns true, then (C, S) defines a point on C k ∩ i( C), where C k is such that C ∩ C = C k ∩ C. Proof:[Proof of Remark A.1] Suppose isSolInSet((C, S),C, ) returns true, and (C, S) defines a point on a connected component C k = C k of C in C. Then C ∩ C contains two connected components of C, leading to a contradiction.

a combinatorial structure, called eXtended Planar Map (XPMap) is proposed to represent such a relationship. Let σ| H i (resp. α| H i ) denote the restriction of the permutation σ (resp. α) defined on H to elements of a subset H i ⊆ H.

  where H = H 1 ∪ . . . ∪ H n is the union of pairwise disjoint non empty sets of half-edges, H0 is an empty set of half-edges representing the infinite face, (H i , σ| H i , α| H i ) are CMaps for all 1 ≤ i ≤ n , ext is a relation that labels one face of each CMap (H i , σ| H i , α| H i ) as the exterior face, and cont is a relation that assigns each exterior face to one non-exterior face of a CMap (H i , σ| H i , α| H i ) or to the infinite face.

	Definition 1.3 ([16]) An eXtended Planar Map (XPMap) is a tuple (H, H0, σ, α,
	ext, cont)

Table 1 :

 1 Computation of the apparent contour of an algebraic surface of degree 13.

					Isolating singularities
	Step	(1)		(2)			(3)			(4)
	d	t	|C x |	t	|C b |	t		m, m∩	t	|D n |, |D c |
	13	13.7	14	3.44	16	9.28	1012, 2478	7.87	4, 7
					Computing topology
	Step	(5)		(6)			(7)			(8)
	d	t	|C y |	t	i	t		m	t	n
	13	1.57	8	33*	5	23.4		1940	0.34	2
			Isolating singularities			Computing topology
	method	Sub-resultant approach	Our approach	Isotop	Our approach
	domain	B0			B0		R 2		B0
	d	t ± σ		n	t ± σ	n	t ± σ	t ± σ	n
	5	0.06 ± 0.04	0	1.60 ± 0.46	0	4.78 ± 0.26	3.17 ± 1.01	0
	7	3.03 ± 3.30	0	3.73 ± 0.40	0	251 ± 16.1	8.13 ± 1.86	0
	9	304 ± 478	1	10.2 ± 5.45	0	-		24.1 ± 16.2	0
	11	>3600		4	36.0 ± 8.53	0	-		75.5 ± 13.5	0
	13	-			43.5 ± 18.1	0	-		90.6 ± 37.9	1
	15	-			97.9 ± 47.1	0	-		169 ± 71.0	0

Table 2 :

 2

Sequential running times t in seconds (averaged over five runs), standard deviations σ, number n of runs requiring multi-precision for isolating the singularities and computing the topology of apparent contours of algebraic surfaces of degree d. B0 is [-1, 1] × [-1, 1] and -means that the process has not been run.

Table 3 .

 3 This table details how the times reported in Table2are distributed among the main steps that our method performs.

		Isolating singularities	Computing topology	total
	Step	(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)
	d	t	t	t	t	t	t	t	t	t
	5	0.14	0.36	1.03	0.06	0.02	0.31	1.22	0.00	3.17
	7	0.64	0.70	2.21	0.17	0.05	1.47	2.84	0.02	8.13
	9	2.85	1.41	4.66	1.27	0.24	5.78	7.91	0.02	24.1
	11	8.08	2.12	16.4	9.44	0.90	10.8	27.6	0.08	75.5
	13	15.6	3.36	16.3	8.14	1.22	17.6*	28.1	0.06	90.6
	15	30.5	4.87	29.1	33.3	2.12	24.7	45.2	0.00	169

Table 3 :

 3 Distribution of times given in Table2between the main steps of our approach.

[START_REF] Arnold | Singularities of Differentiable Maps: Volume I: The Classification of Critical Points Caustics and Wave Fronts[END_REF]

x-critical points; (2) boundary points; (3) δ-approximation of C; (4) solving the ball system; (5) y-critical points; (6) witness boxes singularities; (7) connecting special points; (8) computing the XPMap; Symbol *: for one example, multi-precision was required.

See also http://www.bertinireal.com/.

In Algorithm 4, we need to isolate the solutions of {P = Pz = Pzz = 0}, the Jacobian of which involves derivatives of P of order

Recall that boxes of V b have exactly one coordinate reduced to one point; thus∂C = C if C ∈ V b .

http://www.sagemath.org/

http://subdiv-solver.gforge.inria.fr

http://vegas.loria.fr/isotop/

This research was supported by the ANR JCJC SingCAST (ANR-13-JS02-0006).

Point-face location in a XPMap

We finally propose here a procedure that, given any point p ∈ i(B0) and the XPMap encoding the embedding B ∂ , finds the face of B ∂ to which p belongs. Such a procedure is easily obtained by rewriting Algorithm 10 while generalizing it for XPMaps.

We restate, in this general setting, the preconditions that ensure the modified algorithm will terminate: (c1 ) p is not in B, (c2 ) p is not in the projection of any box of V ∂ and (c3 ) ∀e k ∈ E, p is not in the projection of any box of approx(e k , δ).

Conditions (c2 ) and (c3 ) can be satisfied by refining boxes of V x ∪ V n and δapproximations of elements of E, provided that (c1 ) holds. Checking (c1 ) with a fully numeric method is challenging since it goes back to the problem of deciding zero. In the special case where the coordinates of p are rational numbers and P, Q are polynomials, this can be addressed with symbolic computation, by computing g(z) = gcd(P (x p , y p , z), Q(x p , y p , z)). If degz(g) ≥ 1 then p ∈ B, in particular if degz(g) = 2 and g has two distinct roots, p is a node of B and if g has a double root, p is a cusp of B. Otherwise degz(g) = 0 and p / ∈ B.

Implementation and Results

We implemented the method presented in this paper and tested it to compute the topology of apparent contours of algebraic surfaces of degrees up to 15. We briefly describe this implementation, then we present the experiments we carried out and compare our approach with state-of-the-art methods. Let us recall the main steps that are performed to compute the topology of an apparent contour:

(1) Computing the set C x of x-critical points by calling IsolateSols on S5.

(2) Computing the set C b of boundary points by calling IsolateSols on S1, S2, S3, S4.

(3) Computing a δ-approximation (Ci) m i=1 of C with Algorithm 3. (4) Computing the sets D n of nodes and D c of cusps. This is performed by solving the ball system S b by calling IsolateSols(S b ,.,0) for each Di and each Dij as defined in Proposition 3.1. We denote by m∩ the number of non-empty Dij.

Then nodes and cusps are distinguished with Algorithm 4.

(5) Computing the set C y of y-critical points with IsolateSols({P, Q, R },Ci,0) for each 1 ≤ i ≤ m, then separating the special points (i.e smooth critical points, nodes and cusps of B) in the projection, as described in Section 4.2. Note that all these steps terminate when the assumptions (A1), . . . , (A10) are satisfied. Moreover, the assumptions are generically satisfied and the algorithm terminates if and only if the assumptions are satisfied. Indeed, if (A1), (A2) is not satisfied then the tracking steps in (3) and ( 7) will not terminate. If (A6) or (A9) is not satisfied then A.2 The path-tracking algorithm of [START_REF] Martin | Certified parallelotope continuation for one-manifolds[END_REF] The authors of [START_REF] Martin | Certified parallelotope continuation for one-manifolds[END_REF] propose a parameterized version of the Krawczyk operator to certify that in a given set, in practice a parallelotope, a curve is diffeomorphic to [0, 1]. This operator is the cornerstone of an algorithm that uses a classical adaptive step-length control and constructs iteratively a certified approximation of a connected component of a curve within a bounded initial domain. We state here the specifications of this algorithm in our context of a curve of R 3 defined by two polynomials P and Q.

Ci is diffeomorphic to a one-dimensional manifold. Furthermore, one has either C1 ∩Cm = ∅ and in this case C∩ m i=1 Ci is diffeomorphic to a circle, or C1∩∂C0 = ∅ or Cm∩∂C0 = ∅ and in this case C ∩ m i=1 Ci is diffeomorphic [0, 1]. The initial point p0 is used to construct an initial parallelotope C0, and this step succeeds if p0 is close enough to C. Otherwise the flag failure is returned. Then, d determines the direction in which C is followed. At each step, a new parallelotope Ci with a step-length γ is constructed, with α < γ < β. Several properties are checked on Ci that guarantee the correctness of the algorithm, and γ is decreased until either γ ≤ α or Ci satisfies the latter properties. When γ ≤ α, the algorithm stops and returns the flag failure. In case of success, γ is increased until it reaches β.

When α > 0, the algorithm terminates. When α = 0, it terminates with the flag success, provided that C is smooth in C0 and p0 is sufficiently close to C.

A.3 Meeting the specifications of Algorithm 2

We now show how to modify the algorithm described above to meet the specifications of the procedure Track described in Algorithm 2. These specifications are reproduced here for the sake of readability. Input: A system P = Q = 0 defining a smooth curve C, a domain C0, an implicit point (C 0 , S 0 ) of C, a finite set {(C j , S j )}j of implicit points containing the boundary points of C in C0 and (C 0 , S 0 ), δ > 0. Output: A δ-approximation of the connected component of C∩C0 containing (C 0 , S 0 ) and the set C on of implicit points of {(C j , S j )}j that are on the same connected component than (C 0 , S 0 ). First, the procedure Track takes as input a domain C0 instead of a box C0; this leads in the algorithm to testing inclusion of parallelotopes within a domain rather than within a box. Recall that in our case the domain is possibly unbounded in z. We also fix α to 0. Provided (A2) and (A3) are satisfied (C is smooth and bounded above B0), this modified algorithm terminates with flag success if p0 is sufficiently close to C.

Then, the procedure Track takes as input a finite set {(C j , S j )}j and returns a set C on of implicit points of {(C j , S j )}j that are on the approximated connected component. To achieve this goal, we modify the algorithm so it constructs a pair (Ci, i) with i > 1 such that both C ∩ Ci and C ∩ iCi are diffeomorphic to [0, 1]. According to Remark A.1, if the procedure isSolInSet((C j , S j ),Ci, i) returns true, it is a guarantee that (C j , S j ) defines a point on C in iCi. We also require that at Algorithm 14 TrackFromPoint(< P, Q >,C 0 ,p 0 ,{(C j , S j )} j ,δ) Input: A system P = Q = 0 defining C, an initial domain C0, an initial point p0 ∈ i(C0), a finite set {(C j , S j )}j of implicit points containing the boundary points of C in C0, δ > 0. Output: A flag in {success, failure}, a sequence C enc and a set C on . If the flag is failure, C enc and C on are empty. Otherwise, C enc is the sequence of sets (Ci) m i=1 with m > 1 such that (h(Ci)) m i=1 is a δ-approximation of a connected component C k of C ∩ C0 and C on contains implicit points of {(C j , S j )}j that are on C k . most one point of {(C j , S j )}j is in each Ci; this can be tested in the same way. The latter property guarantees in particular that each Ci intersecting the boundary of the domain contains no more than one boundary point and that the extremities of the approximated connected component of C ∩ C0 are properly identified and reported in C on .

Finally, the procedure Track takes as input δ > 0 and returns a δ-approximation of a connected component. To achieve this, the condition w (h(Ci)) < δ is enforced at each iteration, and this replaces the condition γ ≤ β. Then, when the algorithm described in Section A.2 with d = 1 as input terminates with the flag success and if it is detected that the tracked component is diffeomorphic to [0, 1], another call with d = -1 allows approximation of the whole connected component.

These modifications give rise to the procedure TrackFromPoint(,,,,) specified in Algorithm 14. When assuming (A2) and (A3), this algorithm terminates with the flag success if p0 is sufficiently close to C.

In contrast to the latter procedure, Track takes as input an implicit point (C 0 , S 0 ) on C rather than a point p0 and guarantees that the returned δ-approximation is a δ-approximation of the connected component containing (C 0 , S 0 ). One can use p0 = m C 0 and call the procedure TrackFromPoint(< P, Q >,C0,p0,{(C j , S j )}j,δ). When the returned flag is success, the obtained δ-approximation is not necessarily the one of the connected component containing (C 0 , S 0 ), in particular when two connected components are close to m C 0 . On the other hand, since (C 0 , S 0 ) ∈ {(C j , S j )}j, the returned δ-approximation is the one of the connected component containing (C 0 , S 0 ) if and only if (C 0 , S 0 ) ∈ C on . If it is not the case, C 0 is contracted with K S 0 (C 0 ) and TrackFromPoint(< P, Q >,C0,m C 0 ,{(C j , S j )}j,δ) is called again. This is performed until the returned flag is success and (C 0 , S 0 ) ∈ C on . The termination of this recursion in ensured, since m C 0 becomes arbitrarily close to C and the initial set of the approximation will contain the implicit point (C 0 , S 0 ). Thus, after a finite number of iterations, TrackFromPoint(,,,,) returns the flag success, and (h(Ci)) m i=1 and C on are a suitable output for Track with arguments (< P, Q >, C0, (C 0 , S 0 ), {(C j , S j )}j, δ).