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

Consider a plane curve B defined as the projection of the intersection of two analytic 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 only based on reliable numerical algorithms.

In a 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 our method and show that our reliable numerical approach can handle classes of examples that are not reachable by symbolic methods.

Introduction

Let C be a smooth space curve defined as the intersection of two analytic surfaces P = Q = 0, with P, Q real analytic functions in x, y, z. We aim at computing, in a compact domain B 0 ⊂ R 2 , the geometry and the topology of B = π (x,y) (C) where π (x,y) : R 3 → R 2 denotes the projection in 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 = P z , 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. Fig. 1 shows, for a torus P = 0, its intersection with the surface P z = 0 in bold line and its apparent contour that has cusp and node singularities.

Figure 1: Left: a torus P = 0, in bold line the curve P = P z = 0, its apparent contour, and a zoom zone.

Right: a zoom, with preimages by the projection of cusp and node singularities.

By computing the geometry of B in B 0 , we mean being able to draw B ∩ B 0 with an arbitrary precision. To reach this goal, we use a certified interval path tracker to compute a sequence of boxes (i.e. multidimensional 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, that is, given p ∈ B 0 \ B, find the face of B 0 \ B 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 B 0 induced by the curve B. We use Combinatorial Maps (CMaps) and their extensions to the non-connected case, the eXtended Planar Maps (XPMaps), 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. Then the second step is to compute the topology in these isolating boxes. Finally the tracking of the space curve is used 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 algebraic or analytic variety of R n . Such a variety 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 during a motion the robot passes through a singularity or not.

The paper is organized as follows. Section 2 describes how the curve C is enclosed by tracking. Section 3 recalls the encoding of singularities of B by the ball system and shows how the enclosure of C is used to restrict its solving 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 remaining of this section presents previous work, details our contributions, defines formally our geometric and topological representations and recalls basics about reliable numerical interval solvers.

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 [Hon96, MPS + 06]. One advantage of these methods is that they can handle any type of singularities 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 together with interval arithmetic are able to compute and certify the topology of a non-singular curve [MGGJ13, [START_REF] Kearfott | An interval step control for continuation methods[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] Ojika | Deflation algorithm for the multiple roots of a system of nonlinear equations[END_REF][START_REF] Leykin | Newton's method with deflation for isolated singularities of polynomial systems[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 curves, but several promising approaches have been presented. The subdivision approach presented by [START_REF] Burr | Complete subdivision algorithms ii: Isotopic meshing of singular algebraic curves[END_REF] is an extension to the singular case of the one by [START_REF] Plantinga | Isotopic approximation of implicit curves and surfaces[END_REF] relying on global nonadaptive separation bounds for algebraic systems. This approach can theoretically certify the topology of any singular curve, but due to these 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 if this work has been implemented by [START_REF] Daniel | Bertini: Software for numerical algebraic geometry[END_REF] 1 , the certification of all the algorithm steps appears as a challenge. The numerical approach by [CDTF + 13], 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 higher dimension. According to the classification of singularities of mappings (see [START_REF] Whitney | On singularities of mappings of euclidean spaces. i. mappings of the plane into the plane[END_REF][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] 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] studies the geometric configurations of the space curve that induce singularities on the projected curve, and describes them as the regular solutions of a four dimensional system.

Detailed overview

Our main contribution is the computation of a data structure allowing to locate 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 (C i ) m i=1 enclosing all the connected components of C, that is C ⊂ m i=1 C i . We use a reliable numerical solver (Algo. 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 Algo. 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, and 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 a 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 solving 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. A 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 not reachable by symbolic methods and that multi-precision arithmetic is needed for such difficult examples. More specifically, for the isolation of the singularities, the efficiency of our approach to restrict the solving domain of the ball system is demonstrated.

Notations and definitions

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

(x) -l(x) if x is bounded, ∞ 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 = (x 1 , . . . , x n ) be a box, ∂X denotes the boundary of X and i(X) = X \ ∂X its interior. The width w(X) is defined as max 1≤i≤n w(x i ) and the midpoint m(x) as (m(x 1 ), . . . , m(

x n )) if X is bounded. If ∈ R + , X holds for ( x 1 , . . . , x n ).
Uppercase letters denotes sets of boxes. A domain of R n is a set defined as X \ ( Y∈Y i(Y)) where X is a box of 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 real analytic function P in the variables x 1 , . . . , x n , P xi denotes its partial derivative with respect to x i , and P xixj its derivative with respect to x i and x j . Let P 1 , . . . , P n be n real analytic functions in x 1 , . . . , x n , a solution x of P 1 = . . . = P n = 0 is regular if the jacobian matrix A = [(P i ) xj ] evaluated in x has full rank, otherwise x is singular.

Cursive letters denotes 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 triple G = (V, E, I) where V, E are two finite sets which elements are respectively called vertices and edges 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 couples (resp. pairs).

Geometrical and topological representations

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

0 10 0 1 1 0 0 1 1 00 11

0 1 0 0 1 1 0 1 0 0 1 1 0 1 0 0 1 1 0 1 ∂B 0 Figure 2: Left: a δ-approximation of C in B 0 × R,
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 H 0 = ∅, H 1 = {h + 1 , h - 1 , h + 2 , h - 2 , h + 13 , h - 13 , h + 14 , h - 14 }, H 2 = {h + 3 , h - 3 , . . . , h + 12 , h - 12 } edges: (h + i , h - i ) for 1 ≤ i ≤ 14, vertices: v 1 = (h - 1 , h + 14 , h - 13 ),. . . faces: f 0 = ∅, f 1 = (h + 13 , h + 14 ), f 2 = (h + 1 , h + 2 , h - 14 ), f 3 = (h + 1 , h - 13 , h - 2 ), f 4 = (h + 3 , h - 4 , . . . , h - 12 ), f 5 = (h - 3 , h + 12 , h - 11 , h + 4 ), f 6 = (h + 5 , h - 10 ), f 7 = (h + 6 , h - 9 , h + 8 , h - 7 ). ext: (H 1 , f 1 ), (H 2 , f 4 ) cont: (f 1 , f 0 ), (f 4 , f 3 )
Figure 3: An XPMap encoding the topology of the apparent contour of the torus in B 0 . 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 1 (δ-approximation). Let X be a subset of R n . A sequence of boxes

(X i ) m i=1 is a δ-approximation of X if X ⊂ m i=1 X i and for all 1 ≤ i ≤ m, w(X i ) ≤ δ and X ∩ X i = ∅.
We will consider in particular δ-approximations of C in C 0 = B 0 × R or in a domain C 0 . Obviously the projection of a δ-approximation of C in C 0 is a δ-approximation of B in B 0 , 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. Left part of Fig. 2 represents a δ-approximation of the curve C in C 0 for the example of the torus, and its projection.

We also want to encode the topology of the curve B ∩ B 0 as an embedding of a graph G to be able to answer location queries for points of B 0 . This is achieved by computing a set V of vertices representing points of the curve with at least one vertex per connected component and such that (B ∩ B 0 ) \ V is a set of smooth curves identified as the set E of edges. Right part of Fig. 2 shows an embedding of such a graph G when B ∩ B 0 is the apparent contour of the torus. The embedding of G defined by B ∩ B 0 is encoded by an extended planar map. More precisely, each connected component of (B ∩ B 0 ) ∪ ∂B 0 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 2 ( [LZ13, Köt02]). A Combinatorial Map (CMap) is a triple (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 to a cycle of σ. Cycles of σ encode counterclockwise 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 but one which is called the exterior face of the combinatorial map.

When (B∩B 0 )∪∂B 0 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. [START_REF] Köthe | Xpmaps and topological segmentation-a unified approach to finite topologies in the plane[END_REF] proposes a combinatorial structure, called eXtended Planar Map (XPMap) to represent such a relationship. Let σ| H i (resp. α| H i ) denotes the restriction of the permutation σ (resp. α) defined on H to elements of a subset H i ⊆ H.

Definition 3 ( [Köt02]

). An eXtended Planar Map (XPMap) is a tuple (H, H 0 , σ, α, ext, cont) where H = H 1 ∪ . . . ∪ H n is the union of pairwise disjoint non empty sets of half-edge, H 0 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.

In our case, the infinite face is R 2 , Fig. 3 shows an XPMap characterizing (B ∩B 0 )∪∂B 0 as an embedding of the apparent contour of the torus. The cycle f 1 = (h + 13 , h + 14 ) is the exterior face of the connected component of (B ∩ B 0 ) ∪ ∂B 0 containing ∂B 0 . This exterior face f 1 is contained in the infinite face f 0 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] Neumaier | Interval methods for systems of equations[END_REF][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] 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 (which 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.

Interval arithmetic

Usual arithmetic operations such as additions, multiplications and so on can be extended to intervals and boxes. Letting * being an operator and its interval extension, one has X Y ⊇ {x * y|x ∈ X, y ∈ Y}.

Algorithm 1 IsolateSols(S,X 0 ,δ) Input: A bounded box X 0 of R n , a system S of n analytic functions in n unknowns with only regular solutions in X 0 , the Krawczyk operator K S , 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 ∈ X 0 is a solution of S ⇒ ∃X ∈ X sol ∪ X ind s.t. x ∈ X, • X ∈ X sol ⇒ X ⊂ i(X 0 ) and K S (X) ⊂ i(X),
• X ∈ X ind ⇒ X ∩ ∂X 0 = ∅ and K S (X) ⊂ i(X) and w(X) < δ.

If P 1 is an analytic function and X a box, the evaluation of P 1 over X results in an interval Y that satisfies Y ⊇ {P 1 (x)|x ∈ X}, and in general Y strictly includes the range of values, that is Y over-estimates {P 1 (x)|x ∈ X}. As a consequence, if 0 / ∈ Y then one can certify that P does not vanish in X.

Criteria for existence and uniqueness of solutions

Let P 1 , . . . , P n be n analytic functions in n unknowns, S = {P 1 = 0, . . . , P n = 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 P 1 (x) = . . . = P n (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] Neumaier | Interval methods for systems of equations[END_REF][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] 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 P 1 , . . . , P n 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 solutions of the linear system J(X)v = F (y) where J(X) is the interval evaluation of the jacobian matrix 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] Kearfott | Rigorous global search: continuous problems. Nonconvex optimization and its applications[END_REF] that takes y as the middle of X and an approximate inverse of the derivative of F in y to determine the box V. Let us note K S the Krawczyk operator for the system S, and recall that i(X) is the interior of X. Important results about the Krawczyk operator are:

• K S (X) ⊂ i(X) ⇒ X contains a unique solution of S and the sequence X (0) = X, X (i+1) = K S (X (i) )

converges asymptotically quadratically to this solution.

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

Reliable numerical isolation

Interval evaluation, 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 X 0 . Such methods are described for instance in [Neu90, Kea96, MKC09, Sta95] 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 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 [Neu90, Section 5.6] or [Sta95, 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 more details. The algorithm does not identify solutions that are exactly on the boundary of the input domain X 0 and uses δ-inflation. A consequence is that a solution on ∂X 0 can only be isolated in a box containing a part of the boundary in its interior; the refinement of such a box is stopped thanks 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 X 0 ), that is at distance less than δ from ∂X 0 , 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 X 0 or in other words that X ind is empty.

Enclosing the space curve C

In this section, we introduce Algo. 3 computing a δ-approximation of C. This goal is achieved by first finding at least one point on each connected component of C and 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 = P y Q z -P z Q y . Similarly, p ∈ C is y-critical if it is the solution of the system P = Q = R = 0 where R = P x Q z -P z Q x .
Definition 5. Let C be a box or a domain of 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 Algo. 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 6. Let P, Q define the curve C, E ∈ R[x, y, z] and S be the system

P = Q = E = 0. The couple (C, S) is an implicit point of C if K S (C) ⊂ i(C).
Let p be the unique solution of S in C, we say that p is implicitly defined by (C, S).

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

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

(x × y 0 × z, S) is an implicit point of C if K S0 (x × z) ⊂ i(x × z), where S 0 = {P (x, y 0 , z) = Q(x, y 0 , z) = 0}.
In Sec. 2.2, we characterize a connected component of C ∩ C 0 : 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 present the specifications of Algo. 2, our certified path tracker and in Sec. 2.4, Algo. 3 computing a δ-approximation of C. In Sec. 2.1, we state the assumptions allowing our approach being correct and terminating.

Assumptions

Recall that B 0 = (x 0 , y 0 ) and C 0 = B 0 × R. We define the following assumptions on P, Q and C:

(A 1 ) The curve C is smooth above the box B 0 . (A 2 ) C is compact over B 0 and z 0 is a bounded interval such that C ⊂ B 0 × z 0 .
(A 3 ) P (x * , y, z) = Q(x * , y, z) = 0 has finitely many regular solutions when x * = l(x 0 ) or u(x 0 ), y ∈ y 0 and z ∈ R.

P (x, y * , z) = Q(x, y * , z) = 0 has finitely many regular solutions when y * = l(y 0 ) or u(y 0 ), x ∈ x 0 and z ∈ R.

(A 4 ) P = Q = R = 0 has finitely many regular solutions in C 0 , and no solution in ∂C 0 .

(A 5 ) Over a point of ∂B 0 , P = Q = 0 has only one solution and no solution above its corners.

Connected components of C ∩ C 0 and initial boxes

The following proposition characterizes the topology of the connected components of C ∩ C 0 .

Proposition 7. Assuming (A 1 ), the connected components of C ∩ C 0 are smooth one dimensional manifolds possibly with boundary. In addition, assuming (A 2 ) and (A 3 ), any connected component C k of C ∩C 0 satisfies at least one of the following statements: 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 (A 2 ), it is compact, hence it is diffeomorphic either to [0, 1], or to a circle. Suppose first C k has an intersection with ∂C 0 . From (A 3 ), this intersection is a point, hence C k is diffeomorphic to [0, 1] and has a second intersection with ∂C 0 . Suppose now that C k does not intersect ∂C 0 . 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. 7, the following corollary gives a constructive characterization of a set containing at least one point on each connected component of C ∩ C 0 .

Corollary 8. Consider the following systems of equations:

(S 1 ) P (l(x 0 ), y, z) = Q(l(x 0 ), y, z) = 0, for y ∈ y 0 and z ∈ R. (S 2 ) P (u(x 0 ), y, z) = Q(u(x 0 ), y, z) = 0, for y ∈ y 0 and z ∈ R. (S 3 ) P (x, l(y 0 ), z) = Q(x, l(y 0 ), z) = 0, for x ∈ x 0 and z ∈ R. (S 4 ) P (x, u(y 0 ), z) = Q(x, u(y 0 ), z) = 0, for x ∈ x 0 and z ∈ R. (S 5 ) P (x, y, z) = Q(x, y, z) = R(x, y, z) = 0, for (x, y, z) ∈ C 0 .
Assuming (A 3 ), (A 4 ), the set of solutions of S 1 , . . . , S 4 is finite and is the set of boundary points of C ∩ C 0 , the set of solutions of S 5 is finite and is the set of x-critical points of C ∩ C 0 . Then assuming (A 2 ), the set of solutions of S 1 , . . . , S 5 is a finite set of points in B 0 × z 0 containing at least one point on each connected component of C.

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

The solutions of S 1 are obtained by calling IsolateSols(S 1 ,y 0 × z 0 ,0). From assumptions (A 3 ) and (A 5 ), S 1 has a finite number of solutions on y 0 × R that are regular and has no solution on ∂y 0 × 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 S 1 . A solution B = (y, z) in B sol defines an implicit point (C, < P, Q, E >) where C = l(x 0 ) × y × z and E is the polynomial x -l(x 0 ), as described in Def. 6. When needed, C is contracted with l(x 0 ) × K S1 (π (y,z) (C)). The solutions of S 2 , S 3 and S 4 are obtained similarly.

Certified numerical path-tracking

Our path-tracking procedure Track is specified in Algo. 2, its description is postponed to appendix A. The tracker is used in Algo. 3 to compute the connected components of C ∩ C 0 and later in Sec. 5.1.2 to deduce the topology of B.

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 C 0 , 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 C 0 and (C 0 , S 0 ), δ > 0. Output: A δ-approximation of the connected component of C ∩ C 0 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 ).

Algorithm 3 Compute a δ-approximation of C in C 0 Input: A system P = Q = 0 and an initial domain

C 0 = B 0 × R defining a smooth curve C ∩ C 0 , the set C b of implicit points defining the boundary points of C ∩ C 0 , the set C x of implicit points defining the x-critical points of C ∩ C 0 , δ > 0. Output: A sequence of boxes (C i ) m i=1 that is a δ-approximation of C ∩ C 0 . 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: ((C i ) m k i=m k-1 +1 ,C on ) ← Track(< P, Q >,C 0 ,(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: ((C i ) m k i=m k-1 +1 ,C on ) ← Track(< P, Q >,C 0 ,(C, S),C b ∪ C x ,δ) 10: remove C on ∩ C x from L x 11: Let m = m k 12: return (C i ) m i=1

Computing a δ-approximation of C

In the following, we assume that P, Q satisfy assumptions (A 1 ) and (A 2 ) of Sec. 2.1: C is smooth and compact in C 0 . Based on the procedure Track, Algo. 3 computes a δ-approximation of C ∩ C 0 .

Algo. 3 first computes δ-approximations for every connected components of C that are 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 for the first time in the while loop in Step 3, the list L b contains the two extremities of each connected component of C ∩ C 0 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. 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.

The connected components of C that are diffeomorphic to a circle are approximated in the while loop beginning at Step 7 of Algo. 3. Recall that each connected component of C ∩ C 0 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). [IMP15] also introduces a system called ball system which 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]. Then we 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 Algo. 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 ∩ B 0 .

Assumptions

Consider the following assumptions:

(A 6 ) For any (α, β) in B 0 , the system P (α, β, z) = Q(α, β, z) = 0 has at most 2 real roots counted with multiplicities.

(A 7 ) There is finitely many points (α, β) in B 0 such that P (α, β, z) = Q(α, β, z) = 0 has 2 real roots counted with multiplicities.

(A 8 ) π (x,y) restricted to the curve C is a proper map, that is the inverse image of a compact is compact.

(A 9 ) The singularities of the curve B in B 0 are either nodes or ordinary cusps.

Assumption (A 2 ) given in Sec. 2.1 is a consequence of (A 8 ). Notice that Thom Transversality Theorem implies that (A 1 ), (A 6 ), . . . , (A 9 ) hold for generic analytic maps P, Q defining C (see [Dem00, Th. 3.9.7 and §4.7]).

Ball system

Following a geometric modelling, [IMP15] defines a 4 dimensional system which solutions maps to the singularities of B. In this modelling, two solutions (x, y, z 1 ) and (x, y, z 2 ) of P = Q = 0 (or P = P z = 0) are mapped to the point (x, y, c, r 2 ) with c = (z 1 + z 2 )/2 and r 2 = (z 1 -z 2 ) 2 . Fig. 4 illustrates this mapping for singularities of the apparent contour of a torus. We recall 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 9 ( [IMP15, Lemma 4]). Let P, Q be two analytic functions in x, y, z satisfying the Assumptions (A 1 ), (A 6 ) -(A 8 ), and S be the set of solutions of the so-called ball system:

       1 2 (P (x, y, c + √ r 2 ) + P (x, y, c - √ r 2 )) = 0 1 2 √ r2 (P (x, y, c + √ r 2 ) -P (x, y, c - √ r 2 )) = 0 1 2 (Q(x, y, c + √ r 2 ) + Q(x, y, c - √ r 2 )) = 0 1 2 √ r2 (Q(x, y, c + √ r 2 ) -Q(x, y, c - √ r 2 )) = 0 (1) in B 0 × R × R + . Then π (x,y) (S) = Σ
, where π (x,y) is the projection from R 4 to the (x, y)-plane.

Lemma 10 ( [IMP15, Lemma 5]). Under the Assumptions (A 1 ), (A 6 ) -(A 8 ), all the solutions of the ball system in B 0 × R × R + are regular if and only if (A 9 ) is satisfied.

In addition, the solutions of the ball system with r 2 = 0 map to cusps in 2d and the curve in 3d has a vertical tangent (collinear to the z-axis); whereas the solutions with r 2 = 0 map to nodes.

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

Solving 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 B 0 × R × R + with a subdivision solver. Using the δ-approximation (C i ) 1≤i≤m of C computed by Algo. 3, we propose to reduce significantly the domain of research for the singularities. Indeed, given a singular point σ of B, there exists 1

≤ i ≤ m such that σ ∈ B i , where B i = π (x,y) (C i ).
Hence it is possible to isolate all singularities by solving the ball system within

B i × R × R + , for 1 ≤ i ≤ m.
In addition, Prop. 11 shows how to bound the solving domain in the c and r 2 components.

Proposition 11. Let (C i ) 1≤i≤m be a δ-approximation of C. For 1 ≤ i ≤ m, let C i = (x i , y i , z i ), B i = π (x,y) (C i ) = (x i , y i ), and for 1 ≤ i < j ≤ m, let B ij = (x ij , y ij ) = B i ∩ B j
and consider the sets:

• D i = (x i , y i , z i , [0, ( w(zi)
2 ) 2 ]),

• D ij = (x ij , y ij , (zi+zj ) 2
, (

(zi-zj ) 2 ) 2 )
Then, under the Assumptions (A 1 ), (A 6 ) -(A 9 ), all solutions of the ball system lie in

( 1≤i≤m D i ) ∪ ( 1≤i<j≤m D ij ).
Proof. Let p = (x p , y p , c p , r p ) ∈ B 0 × R × R + be a solution of the ball system and σ = π (x,y) (p) the corresponding singularity in Σ. From Assumption (A 9 ), σ is either an ordinary cusp or a node. Suppose first it is an ordinary cusp. Then r p = 0, and σ is the projection of a single point p = (x p , y p , c p ) of C. Hence there exists 1 ≤ i ≤ m such that p ∈ C i . As a consequence we have c p ∈ z i and p ∈ D i (see Fig. 5).

Suppose now σ is a node. Then r p > 0, and σ is the projection of two points p -= (x p , y p , c p -√ r p ) and

p + = (x p , y p , c p + √ r p ) of C. Hence there exist 1 ≤ i ≤ m and 1 ≤ j ≤ m such that p -∈ C i and p + ∈ C j .
If i = j, we have c p ∈ z i and r p ∈ [0, ( w(zi) 2 ) 2 ], and finally p ∈ D i . If i = j (this case is illustrated in Fig. 5), c p lies in

zi+zj 2
that is the center of the two intervals z i and z j , and r p lies in ( zi-zj 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 smooth analytic surface, then its singularities generically also include cusps. We introduce Algo. 4 to distinguish these two types of singularities.

B i , nodes in boxes B i ∩ B j .
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 ∈ 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 C P ∩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 (α, β, γ) = P z (α, β, γ) = P zz (α, β, γ) = 0. According to assumption (A 6 ), γ is a triple root of P (α, β, z), that is P zzz (α, β, γ) = 0. The regularity of the curve C P ∩Pz thus implies that (α, β, γ) is a regular solution of the system P = P z = P zz = 0 (see [IMP17, Lemma 10]). Noting K (P,Pz,Pzz) the Krawczyk operator for the latter system, the test on Line 2 of Algo. 4 is eventually true for a small enough box those projection contains a cusp.

While refining the box D, Algo. 4 hence 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 * )) 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. A 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 case of a node is detailed and the other cases are only sketched since they are similar.

Boundary points

We define the boundary points of C as the set C ∩ (∂B 0 × R) and the boundary points of B as the set B ∩ ∂B 0 . With Assumptions (A 3 ) -(A 5 ), the boundary points of C and B are in one-to-one correspondence, are not x nor y-critical of C nor B and are not the corners of B 0 . The isolating boxes of boundary points of C can thus be refined until their projection on ∂B 0 are disjoint.

Preprocessing of nodes, cusps and x-extreme points

We add the following assumptions. 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(C 0 ) and • x and y-critical points of C are in one-to-one correspondence with cusps of B.

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

• y-critical not x-critical points of C 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 B 0 × z 0 . Assumption (A 11 ) 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 cusps 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. We call special points of B a point that is either x-extreme, node or cusp. 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 12. A witness box for a special point (x-extreme, node or cusp) of the curve B is a box containing this point and 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. Right part of Fig. 6 shows a witness box for a node.

Topology at a node singularity

Algo. 6 computes a witness box for a node. The idea is that if the curve C does not contain x or y-critical point above an isolating box B of the node, B is a 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 right part of Fig. 6). Proof. Since there is no x-critical point of C in B×R, the system

P (x = l(x B ), y, z) = Q(x = l(x B ), y, z) = 0 (resp. P (x = u(x B ), y, z) = Q(x = u(x B ), y, z) = 0
) has only regular solutions on y B × R. In addition, since w(y B )/4 > 0, the procedure IsolateSols called in Step 2 (resp.

Step 3) of Algo. 5 terminates even if solutions lie on ∂y B × 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(y B ), z) = Q(x, y = l(y B ), z) = 0 and P (x, y = u(y B ), z) = Q(x, y = u(y B ), z) = 0 have no solution on ∂x B × R. Hence the procedures IsolateSols called in Step 6 and 7 of Algo. 5 terminates if and only if the two later systems have only regular solutions on x B × R, and this is the case since C has no y-critical point over B. (i) follows.

Suppose Algo. 5 returns True. This implies that B crosses the boundary of B in exactly four distinct points. Since C does not have 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(x B ) × ∂y B and its width should be less than w(y B )/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 does not intersect S(∂B) nor until B ⊆ M (B n ) 9: until B n ⊆ B and IsWitnessNodeBox(B n ) 10: return B n 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.

1: Let B = (x B , y B ) 2: (X sol 1 , X ind 1 ) = IsolateSols(P (x = l(x B ), y, z) = Q(x = l(x B ), y, z) = 0,y B × z 0 ,w(y B )/4) 3: (X sol 2 , X ind 2 ) = IsolateSols(P (x = u(x B ), y, z) = Q(x = u(x B ), y, z) = 0,y B × z 0 ,w(y B )/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(y B ), z) = Q(x, y = l(y B ), z) = 0,x B × z 0 ,0) 7: (X sol 4 , X ind 4 ) = IsolateSols(P (x, y = u(y B ), z) = Q(x, y = u(y B ), z) = 0,x B × z 0 ,0) 8: if not |X sol 1 |+|X sol 2 |+|X sol 3 |+|X sol 4 |=
Note that Algo. 5 may return false negative: a box that is witness may not be classified as such. This can happen when the curve crosses the box near its corners. The idea of Algo. 6 is to refine the box of a node to avoid such a case, that is such that property (iii) of Prop. 13 eventually holds. The sequence B n of boxes constructed in Algo. 6 is illustrated in the right part of Fig. 6. Proposition 14. Algo. 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 a n . 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.

Let

B n σ = σ + [-a n , a n ] × [-2b n , 2b n ],
this is a strictly decreasing sequence for the inclusion and for

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, B n will be reported as witness by Algo. 6. We will show that B n satisfies the sufficient condition (iii) of Prop. 13.

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 b0 a0 of the box B n is increasing. For n large enough, B 1 (resp. B 2 ) is thus crossing 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 Prop. 13 is hence satisfied and the returned box B n is a witness box.

Topology at a cusp singularity

Algo. 5 and 6 are easily adapted for a cusp. The input of Algo. 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 Algo. 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 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 exactly once E(∂B), once W (∂B) 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. Or 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 in parallel the two algorithms 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 the 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 ∩ C 0 together with the local topology at the special points of its projection B ∩ B 0 are combined to compute the global topology of B ∩ B 0 . We use a XPMap to encode this topology and design a point location algorithm, that is given p ∈ B 0 \ B, find the connected component of B 0 \ B to which p belongs.

In Sec. 5.1, we compute a graph G such that B ∩ B 0 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 B 0 , we consider the curve B ∂ = (B ∩ B 0 ) ∪ ∂B 0 as the embedding of a graph G ∂ . The graph G ∂ and its embedding B ∂ have in general several connected components. Sec. 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. Sec. 5.3 shows how to answer location queries with respect to one connected component of B ∂ . In Sec. 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 the nodes of B ∩ B 0 computed in Section 4. We define V x and V n as the cylinders above witness boxes: Remark that implicit points of C b define connections of C ∩ C 0 in C 0 . 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 ∩ C 0 in C are given by the set connect(C) = {. . . , (C i , S i ), . . .} of disjoint implicit points (i.e. C i 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 8).

V x = {B × R | B ∈ B x } and V n = {B × R | B ∈ B n }. We note C 0 (resp. B 0 ) the domain of R 3 (resp. R 2 ) defined as C 0 \ C∈V x ∪V n i(C) (resp. π (x,y) (C 0 )). Recall that
Let V be the set V x ∪ V n ∪ V b . Let E be the set of connected components of C ∩ C 0 . 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 ∩ C 0 in C and an in-connection of C ∩ C 0 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 16. Let G be the graph (V, E, I). The curve B∩B 0 is an embedding of G, seen as a non-directed graph.

We note v 1 , v 2 , . . . the vertices of G (i.e. the boxes of V = V x ∪ V n ∪ V b ) and e 1 , e 2 , . . . the edges of G, (i.e. the connected components of C ∩ C 0 ). The graph G representing the apparent contour of the torus is shown in the left part of Fig. 7.

We show in Sec. 5.1.1 how to distinguish in and out connections of C ∩ C 0 in boxes of V and in Sec. 5.1.2 how to compute the incidence relation I by computing δ-approximations of the connected components of C ∩ C 0 . We consider that G is equipped with an application approx that maps each e ∈ E to its δapproximation. Furthermore, approx satisfies the three following properties: 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(B 0 ) is a clockwise walk.

G ∂ = (V ∂ , E ∂ , I ∂ ) of which the curve B ∂ = (B ∩ B 0 ) ∪ ∂B 0 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 B 0 and E ∂ = E ∪ {e ∂ } where e ∂ is ∂B 0 \ B ∂ . Otherwise (V b is not empty), V ∂ = V and E ∂ = E ∪ E where E is the set of connected components of ∂B 0 \ C∈V b i(π (x,y) (C)). Then I ∂ (e) is defined as I(e) if e ∈ E.
It is clear that the curve B ∂ is an embedding of G ∂ . The left part of Fig. 7 shows, in dashed arrows, the edges of E ∂ \ E.

In and out connections

We characterize in and out connections of C ∩ C 0 in a box of V x ∪ V n or in C 0 using the direction of the tangent of 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 17. Let c be a connection of C ∩C 0 in a box C ∈ V x ∪V n , and B = π (x,y) (C). Let t = (t x , t y , t z ) be a tangent vector of C at c. Then (i) t x = 0 and t y = 0;

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

Proof. According to Section 4, the boxes in V x ∪ V n do not contain any x or y-critical point on their boundary so that 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).

Prop. 17 is easily adapted for boundary points that is when c is a connection of C ∩ C 0 in C 0 . From assumptions (A 4 ), (A 11 ), ∂C 0 contains neither x nor y-critical point 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 this is now the tangent vector pointing inside the box that is relevant.

Prop. 17 is used to decide if a connection c implicitly defined by the couple (C i , S i ) is an out or an in connection as follows. C i can be contracted with the Krawczyk operator K Si until 0 / ∈ R(C i ) and 0 / ∈ R (C i ), where R and R as in Def. 4 are the x and y component of a tangent vector. This process terminates from (i) of Prop. 17, and (ii), (iii), (ii ), (iii ) allows to conclude. Right part of Fig. 7 illustrates this in the case of a node singularity. For two among the four boxes isolating the connections in its witness box, the cones containing the tangent vector (R(C i ), R (C i )) are drawn with dashed lines.

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∩C 0 in boxes of V x ∪ V n and in C 0 . For each out-connection (C i , S i ) ∈ C + , the process Track(< P, Q >,C 0 ,(C i , S i ),C -∪{(C i , S i )},δ) defined in Algo. 2 is performed. This process terminates and returns a δ-approximation (C e i ) me 1 of the connected component e of C ∩ C 0 having (C i , S i ) as left-most point and the set C on containing exactly the implicit point (C j , S j ) 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 (C i , S i ) (resp. (C j , S j )) 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) holds, the ball system can have solutions only in the 4-dimensional boxes constructed as in Prop. 11 from C e 1 , C e 1 , C e me , C e m e . 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 with 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 solution of the system (S 5 ) as defined in Corollary 8 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 components 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(B 0 ). 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 [LZ13][ §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 order 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 be 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 Fig. 7).

Let now p be a boundary point of B k ∩ B 0 . 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 Fig. 7,

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 | v i ∈ V k }, E k = {i | e i ∈ E k } and (O i ) 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 (e i ) = (v j , v l ) with j, l ∈ V k and let O j = O j ∪ {h + i }, O l = O l ∪ {h - i } 6: for i ∈ V k do 7:
Order counter-clockwise the connections of O i around ∂(π (x,y) (v i ))

8: Let σ k be the permutation on H k which orbits are O i , for i ∈ V k . 9: return (H k , σ k , α k ) 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 , 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. 1: Let v = (. . . , h * i , . . .) be the orbit of σ k corresponding to C 2: Let {. . . , c i = (C i , S i ), . . .} be the connections of C in C associated to v 3: Contract each c i with C i = K S i (C i )
until their projections are disjoint and p ∈ π (x,y) (C i ) 4: Order counter-clockwise connections c i and p ∂B and let (c 1 , p, c 2 , . . .) be such an ordering 5: Let f be the orbit of ϕ k containing h * 1 6: return f and let p be the boundary point of which v 1 is a witness box. Since p is on E * (∂B), a counter-clockwise ordering is necessarily (e 1 , e 14 , e 13 ).

Algo. 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 e i ∈ E k are associated two half-edges h + i and h - i , such that h + i is oriented as e i , i.e. if I k (e i ) = (v j , v l ), h + i leaves v j and h - i leaves v l . In other words, h + i represents a walk from left to right along π (x,y) (e i ). 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, Algo. 8 identifies the face containing a point located on the boundary of a witness box. In the algorithm, h * stands for an 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 Fig. 8 illustrates this procedure when π (x,y) (C)

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
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 Let f be the result of Algo. 8 with input (H k , σ k , α k ), C and p

7: return f , 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 (c 1 , p, c 2 , . . .) is a counter-clockwise 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 Algo. 9 computing 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 counterclockwise cycle of half-edges. The case where k = 1 is directly adressed since B 1 ∂ contains ∂B 0 . The exterior face of B 1 ∂ is exactly R 2 \ B 0 and is described by the orbit of ϕ 1 containing the half-edges associated with edges of E ∂ \ E (see Steps 1, 2 and 3 of Algo. 9). As an illustration, the exterior face of B 1 ∂ on Fig. 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 18. 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 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 box 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 Prop. 7. 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 ) < p x , where p x is the x-coordinate of p, and C is not a leftmost box enclosing B k ∂ . Consider now a point p ∈ W * (∂B)

Algorithm 10 Point-face location in a CMap Input:

A point p = (x p , y p ) ∈ i(B 0 ). The CMap (H k , σ k , α k ) encoding the embedding 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) ∀e j ∈ E k , p is not in the projection of any box of approx(e j , δ). Output: The orbit of ϕ k describing the face of B k ∂ that contains p.

1: Let S = {P (x p , y, z) = Q(x p , y, z) = 0}, B * = (y * , z * ) = ∅ and f * = ∅ 2: for e j ∈ E k \ E do
//e j is not a boundary edge 3:

Let X sol be the result of Algo. 11 with input e j , p and S 4:

if X sol = ∅ then 5:

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

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

12:

Let y * = y 13:

Let f * be the result of Algo. 8 with input (H k , σ k , α k ), C and the point (x p , 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 Algo. 9

18: return f * such that p ∈ B. Since C is the witness box, one has necessarily p ∈ B k ∂ and a contradiction follows. Hence (iii) holds.

Algo. 9 uses Rem. 18 in Step 5. From point (iii), the call to the procedure described in Algo. 8 terminates and returns the orbit of ϕ k describing the exterior face of B k ∂ .

Point-face location in a CMap

We now describe Algo. 10 for the point-face location in a CMap. The preconditions (c1), (c2), (c3) ensure its termination. In order to show the correctness of Algo. 10, we consider the segment B seg = x p × [y p , u(y 0 )] and the set C seg = B seg × R where B 0 = (x 0 , y 0 ). 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) (e j ) of B k ∂ (i.e. a connected component of B k ∂ ∩ B 0 ) 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 ,
we use Algo. 8 to determine the face containing the intersection of B seg with S(∂π (x,y) (C)). Let us now give further details.

Proof of correctness of Algo. 10. Suppose first that C seg does not intersect any box of approx(e i , δ) where Algorithm 11 Edge intersection Input: An edge e ∈ E k , a point p = (x p , y p ), the system S = {P (x p , y, z) = Q(x p , y, z) = 0}. Preconditions: (c1) p is not in B k ∂ and (c3) ∀e j ∈ E k , p is not in the projection of any box of approx(e j , δ). Output: A set X sol containing the unique intersection of e with C seg if it exists, X sol = ∅ otherwise.

e i ∈ E k nor any box of V k . Then B seg does not intersect B k ∂ \ ∂B 0 . If k = 1, i.e. B k ∂ is the component of B ∂ containing ∂B 0 ,
1: Let (C i ) m i=1 = approx(e, δ) and ∀1 ≤ i ≤ m, C i = (x i , y i , z i ) and X i = (y i , z i ) 2: for 1 ≤ i ≤ m s.t. x p ∈ x i and l(y i ) > y p do //C i is located above p 3: Let (X sol , X ind ) = IsolateSols(S,1.1X i ,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 e j ∈ E k \ E are computed in the loop beginning in Step 2. Algo. 11 computes, if it exists, the unique intersection q = (x p , y q , z q ) of e j with C seg and returns a non-empty set containing B * j ⊂ C seg with (y q , z q ) ∈ 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 e j 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, S 4 ) where S 4 is defined as in Corollary 8 and B = (x, z) is such that x p ∈ x. Letting e j ∈ E k be such that (B, S 4 ) is a connexion of e j and approx(e j , δ) = (C i ) mj i=1 , it follows from property (a2) of approx that the boundary point implicitly defined by (B, S 4 ) is either in C 1 or in C mj . If the intersection (x p , u(y 0 )) is in C 1 or in C mj , the appropriated face has been determined above. Otherwise, (B, S 4 ) can be refined with K S4 until (x p , u(y 0 )) / ∈ 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 (x p , l(y)). Remark that (x p , l(y)) / ∈ B k ∂ , otherwise (x p , 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 Algo. 8 in Step 13 satisfy its preconditions.

Proof of termination of Algo. 10. It is shown below that if its input arguments satisfy preconditions (c1) and (c3), Algo. 11 terminates. We have already shown that Algo. 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 of correctness and termination of Algo. 11. Let approx(e, δ) = (C i ) m i=1 and C i = (x i , y i , z i ) and X i = (y i , z i ). From properties (a2) and (a3) of approx, e = C ∩ m i=1 C i contains no x-extreme point of C and is a smooth x-monotone curve. Hence its intersections with C seg are regular solutions of S. From property (a1) of approx, for any 1 ≤ i ≤ m, 1.1C i does not intersect any box of any other approximation. Hence for any 1 ≤ i ≤ m, 1.1X i contains at most one solution of S. If 1.1X i 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.1X i ,min( 0.1w(xi) 2 , 0.1w(yi)

2

)) in

Step 3 of Algo. 11 terminates, and since any δ-approximation contains a finite number of boxes, Algo. 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.1X i is a point of the edge represented by e. Reciprocally since e ⊂ m i=1 C i (approx(e, δ) is a δ-approximation of e), there exists 1 ≤ i ≤ m such that e ∩ C seg is a solution of S in X i satisfying x p ∈ x i and l(y i ) > y p from preconditions (c1) and (c3) and X sol contains the corresponding solution. The correctness of Algo. 11 follows.

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 compute the relation cont that assigns each exterior face to a non-exterior face of another CMap or to the infinite face. We describe in Sec. 5.4.1 the construction of the relation cont. Point-face location in an XPMap is discussed in Sec. 5.4.2.

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 as Rem. 18.

Remark 19. 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 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 ) where C k = (x k , y k , z k ) is a leftmost box enclosing B k ∂ .
As a consequence of Rem. 19, 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 ∂ .
Algo. 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 \ B 0 . We name f 0 the empty orbit of ϕ 0 . The loop beginning in Step 2 computes CMaps (H k , σ k , α k ) for B k for 1 ≤ k ≤ n , together with their exterior faces f k and the leftmost box enclosing B k . Since B 1 ∂ is the connected component of B ∂ containing ∂B 0 , 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 , f 1 )} and cont = {(f 0 , f 1 )}.
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 Let (H k , σ k , α k ) be the result of Algo. 7 with input G k 4:

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

5: Let H = H 1 , σ = σ 1 , α = α 1 , ext = {(H 1 , f 1 )} and cont = {(f 0 , f 1 )} 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 Algo. 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) of B <k ∂ containing f k , i.e.
box of V 1 ∂ ∪ . . . ∪ V k-1 ∂ . Then ∀e i ∈ E 1 ∂ ∪ . . . ∪ E k-1 ∂
, p is not in any box of approx(e i , δ) as a consequence of the property (ii) of the application approx.

Point-face location in a XPMap

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

We restate, in this general setting, the preconditions yielding the modified algorithm to 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 , δ).

(c2 ) and (c3 ) can be satisfied by refining boxes of V x ∪V n and δ-approximations of elements of E, providing 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 deg z (g) ≥ 1 then p ∈ B, in particular if deg z (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 deg z (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 S 5 .

(2) Computing the set C b of boundary points by calling IsolateSols on S 1 , S 2 , S 3 , S 4 .

(3) Computing a δ-approximation (C i ) m i=1 of C with Algo. 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 D i and each D ij as defined in Prop. 11. We note m ∩ the number of non-empty D ij . Then nodes and cusps are distinguished with Algo. 4.

(5) Computing the set C y of y-critical points with IsolateSols(< P, Q, R >,C i ,0) for each 1 ≤ i ≤ m, then separating the special points (i.e smooth critical points, nodes and cusps of B) in projection as described in Sec. 4.2.

(6) Computing the sets B x and B n of witness boxes for x-critical (possibly cusp) points and nodes of B.

(7) Computing the graph G ∂ : the projections of connections of x-critical points and nodes of B are made disjoint, then a δ-approximation for each connected components of C ∩ C 0 is computed with calls to Track while ensuring that conditions (a1), (a2) and (a3) described in Sec. 5.1 are satisfied. Let us call m the sum of the number of boxes of each such δ-approximation.

(8) Embedding the graph with Algo. 12.

Our implementation

The two angular stones of our implementation are the procedure IsolateSols specified in Algo. 1 and the procedure Track specified in Algo. 3 and described in appendix A.

Our implementation of the procedure IsolateSols with subdivision, centered-form 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 as the package subdivision solver3 for the mathematical software SageMath4 . Solving a system of polynomial equations in an unbounded box reduces to solve several transformed systems in bounded boxes (see. [START_REF] Neumaier | Interval methods for systems of equations[END_REF]Section 5.6] or [Sta95, Section 5.10]). In particular, we don't need to know in advance the bounded interval z 0 for the z-component. Solving a system on a box of the shape x × y × R reduces to solve 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 ∩ C 0 yields bounds for the z-coordinates of C ∩ C 0 . 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 (7), the initial value δ = 1 is used as input for Track. In step (7), 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

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

B 0 = [-1, 1] × [-1, 1].
The package Isotop5 for Maple computes the topology of a plane curve in R 2 , see [CLP + 10]. Here, Isotop is used to compute the topology of the resultant of P and P z 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. Fig. 9 displays its apparent contour and Tab. 1 details each step of the computation. Columns t give the running times. For Step (6), 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 * notifies that multi-precision was required to address this step.

Isolating singularities

Computing topology method Sub-resultant approach Our approach Isotop Our approach domain Table 2. The first group of columns reports on the isolation of the singularities of the apparent contour in B 0 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 sub-resultant 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 bit-size 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.

B 0 B 0 R 2 B 0 d t ± σ n t ± σ n t ± σ t ±
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 P z 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 bitsize 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 to use 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 upon 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 to filter the domain where the ball system is solved, and the results of Table 2 shows the efficiency of this strategy. For the computation of the singularities in Steps to (4): Step (2) solves systems in 2d and is the least time consuming; Step (1) solves systems in 3d in the large domain B 0 × R which is more time consuming;

Step (4) solves systems in 4d 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 solving domain is dramatic.

In Step (5), the y-critical points are found by solving a system in 3d 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 Algo. 6 that contracts a box containing a node and Algo. 5 that checks if the resulting box is a witness box, and their equivalents for x-extreme points and cusps. In almost all cases we experimented, Algo. 5 succeeded with one iteration of Algo. 6. The time for Step (6) 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 Algo. 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 for defining these systems. For this example, multi-precision arithmetic was required to carry out Step (6). Determining a witness box for an x-extreme point for this particular example required 5 iterations in the equivalent of Algo. 6 for x-extreme points before the equivalent of Algo. 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), Algo. 5 involves isolating the intersections of C on the boundaries of a cylinder. Step (2) requires two calls to the numerical 0-dim solver whereas in Algo. 5 we use the bounds for the z-component given by the δ-approximation of C ∩ C 0 . Hence Algo. 5 requires roughly speaking half the time of Step (2) to execute.

Among Steps (1) to (8), Step (7) is the most time consuming. It consists in two substeps. First, δapproximations of connected components of C ∩ C 0 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 (8). See also column n in the table in appendix B.

Algorithm 13 isSolInSet((C, S),C, ) Input: An implicit point (C, S) defining a point p, a set C and a real number > 1. For the sake of generality, we consider here approximating sets C i 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).

Output: true if p ∈ C, false if p / ∈ C. 1: C = C 2: while not ( (C ⊂ i(C )) or (C ∩ C = ∅) ) do 3: C = K S (C) 4: if C ⊂ i(C )
Sec. A.1 introduces Algo. 13 using -inflation to decide if an implicit point lies in a set C. Sec. A.2 recalls the parallelotope path-tracking algorithm of [START_REF] Martin | Certified parallelotope continuation for one-manifolds[END_REF], and Sec. A.3 shows how it is adapted to meet the specifications of Algo. 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 in 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 Algo. 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 of termination and correctness of Algo. 13. Assume > 1 and let C be C. After applying a finite number of times Step 3 of Algo. 13, 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 Algo. 13 terminates, so does Algo. 13.

A.2 The path-tracking algorithm of [MGGJ13]

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 corner-stone 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. The initial point p 0 is used to construct an initial parallelotope C 0 and this step succeeds if p 0 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 C i with a step-length γ is constructed, with α < γ < β. Several properties are checked on C i that guarantee the correctness of the algorithm and γ is decreased until either γ ≤ α or C i satisfies the latter properties. When γ ≤ α, the algorithm stops and returns the flag failure. γ is increased in case of success until it reaches β.

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

A.3 Meeting the specifications of Algo. 2

We now show how to modify the algorithm described above to meet the specifications of the procedure Track described in Algo. 2 that are reproduced here for the sake of readability. Input: A system P = Q = 0 defining a smooth curve C, a domain C 0 , 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 C 0 and (C 0 , S 0 ), δ > 0. Output: A δ-approximation of the connected component of C ∩ C 0 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 C 0 instead of a box C 0 , this leads to test in the algorithm 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. Providing that (A 2 ) and (A 3 ) are satisfied (C is smooth and bounded above B 0 ), this modified algorithm terminates with flag success if p 0 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 that it constructs a pair (C i , i ) with i > 1 such that both C ∩ C i and C ∩ i C i are diffeomorphic to [0, 1]. According to Rem. 21, if the procedure isSolInSet((C j , S j ),C i , i ) returns true, it is a guarantee that (C j , S j ) defines a point on C in i C i . We also require that at most one point of {(C j , S j )} j is in each C i , this can be tested in the same way. The latter property guarantees in particular that each C i intersecting the boundary of the domain contains no more than one boundary point and that the extremities of the approximated connected component of C ∩ C 0 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(C i )) < δ is enforced at each iteration, and this replaces the condition γ ≤ β. Then, when the algorithm described in Sec. 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], an other call with d = -1 allows to approximate the whole connected component.

These modifications give rise to the procedure TrackFromPoint(,,,,) specified in Algo. 14. When assuming (A 2 ) and (A 3 ), this algorithm terminates with the flag success if p 0 is sufficiently close to C.

In contrast with the latter procedure, Track takes as input an implicit point (C 0 , S 0 ) on C rather than a point p 0 and guarantees that the returned δ-approximation is a δ-approximation of the connected

B Details for results

Isolating singularities

Computing

Definition 4 .

 4 A 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

(

  a) C k has exactly two boundary points, (b) C k has at least two x-critical points.

Figure 5 :

 5 Figure 5: Some boxes and their projections containing singularities of B. Cusps singularities are in boxes

2 :

 2 if K (P,Pz,Pzz) ((x, y, c)) ⊂ i((x, y, c)) then return cusp 3: D = K b (D) 4: return node 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) the associated singularity of B. If 0 /

(A 10 )

 10 The two points of C above a node of B are not x-critical nor y-critical points of C.(A 11 ) P = Q = R = 0 has finitely many regular solutions in C 0 and none in ∂C 0 , they are the y-critical points of the curve C.

Figure 6 :

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

Algorithm 5

 5 IsWitnessNodeBox(B) Input: A box B containing a unique node of B, s.t. 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).

D

  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 σ. 1: Let D = D and B = π (x,y) (D) 2: Let n = 0, a n = w(B), b n = w(B) and B n = B 3: repeat 4: n = n + 1, a n = a n-1 /4 and b n = b n-1 /2 = K b (D ), B = π (x,y) (D ) 7: B n = m(B ) + ([-a n /2, a n /2], [-b n /2, b n /2]) 8:

  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 ∩ C 0 are isolated in boxes of V x , the connected components of C ∩ C 0 are diffeomorphic to [0, 1] and x-monotone. Their endpoints are C ∩ C∈V x ∪V n ∪V b ∂C 2 . Hence a point of the latter set is either the left-most or the right-most point of a connected component of C ∩ C 0 , and we state the following definition. Definition 15. Let C be a box of V x ∪ V n or the box C 0 . We call a point c of C ∩ ∂C a connection of C ∩ C 0 in C. If c is the left-most point of a connected component of C ∩ C 0 , we say that c is an out-connection. Otherwise we call it an in-connection.

(

  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 ≤ m e , 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 = m e 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 In the following, we consider the faces of B ∩ B 0 in B 0 , that are the connected components of B 0 \ B. It is thus convenient to represent explicitly the boundary of B 0 in a graph

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 ∩ C 0 in boxes of V . Thick lines represent the curve B in witness boxes. Dashed arrows represent the clockwise walk on ∂B 0 . Right: A witness box B for a node, and in bold the four boxes isolating the connections of B in B.

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).

3 :

 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 * (∂(π (x,y) (C)))6:

  j ∩ y * = ∅ do let B * j = K S (B * j ) and B * = K S (B * ) 9:

  we let p be the point (x p , u(y 0 )) ∈ ∂B 0 . Then p lies on a boundary edge e i 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.

∂.

  Step 8. Suppose the relation cont has been properly constructed for B 1 ∂ , . . . , B k-1 ∂ , i.e. (H, H 0 , σ, α, ext, cont) encodes the faces of B <k ∂ = B 1 ∂ ∪ . . . ∪ B k-1 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 Rem. 18. Due to the property of a witness box, p lies in the same face of B <k ∂ than B k ∂ . The while loop inStep 11 finds this face by calling Algo. 10 with input p, (H k-i , σ k-i , α k-i ) and f k-i for increasing i. Rem. 20 below shows that this input meets the preconditions of Algo. 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 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 ∂ 1: Let H 0 = ∅ and f 0 represent the infinite face 2: for 1 ≤ k ≤ n do 3:

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.

  then return true 5: return false A Implementation of Algo. 2 Several certified numerical path-tracking algorithms can be found in the literature [KX94, FM07, MGGJ13]. They approximate a smooth connected curve in a bounded box with a sequence of sets (C i ) m k i=1 . The sets C i are in general boxes as in [KX94, FM07]. It has recently been proposed to use parallelotopes instead of boxes [MGGJ13]; 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.

Remark 21 .

 21 Let (C, S) be a point of C, and C and > 1 be such that both C ∩C andC ∩ 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 of Rem. 21. 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, which rises a contradiction.

  Input: A system P = Q = 0 defining C, an initial bounded box C 0 , an initial point p 0 ∈ C 0 , two strictly positive real numbers α, β, an integer d ∈ {-1, 1}. Output: A flag in {success, failure} and a sequence C enc . If the flag is failure, C enc is empty. Otherwise, C enc is a sequence of sets(C i ) m i=1 with m > 1 such that C ∩ C i is diffeomorphic to [0, 1] for each i, C i ⊂ i(C 0 ) for i / ∈ {1, m} and C ∩ m i=1 C i is diffeomorphic to a one-dimensional manifold. Furthermore, one has either C 1 ∩ C m = ∅ and in this case C ∩ m i=1 C i is diffeomorphic to a circle, or C 1 ∩ ∂C 0 = ∅ or C m ∩ ∂C 0 = ∅ and in this case C ∩ m i=1 C i is diffeomorphic [0, 1].

  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 Since witness boxes of V are pairwise disjoint, p is not in any

	∂	.
	Remark 20. Let C be a leftmost box enclosing B k ∂ , and p be a point on W * (∂(π (x,y) (C))). From point (iii)
	of Rem. 18, p / ∈ B.

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. B 0 is [-1, 1] × [-1, 1] and -means that the process has not been run.

	σ	n

Table 1 :

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

Table 3 .

 3 This table details how the times reported inTable 2 are distributed among the main steps that our method performs.

	Computing topology	total

Table 3 :

 3 Distribution of times given in Table2between the main steps of our approach. (1) 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 www.bertinireal.com.

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

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

http://www.sagemath.org/

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

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

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 C 0 , an initial point p 0 ∈ i(C 0 ), a finite set {(C j , S j )} j of implicit points containing the boundary points of C in C 0 , δ > 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

component containing (C 0 , S 0 ). One can use p 0 = m(C 0 ) and call the procedure TrackFromPoint(< P, Q >,C 0 ,p 0 ,{(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 >,C 0 ,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(C i )) m i=1 and C on are a suitable output for Track(< P, Q >,C 0 ,(C 0 , S 0 ),{(C j , S j )} j ,δ).