Matthew Trager

Martial Hebert

Jean Ponce

Consistency of silhouettes and their duals

Silhouettes provide rich information on threedimensional shape, since the intersection of the associated visual cones generates the "visual hull", which encloses and approximates the original shape. However, not all silhouettes can actually be projections of the same object in space: this simple observation has implications in object recognition and multi-view segmentation, and has been (often implicitly) used as a basis for camera calibration. In this paper, we investigate the conditions for multiple silhouettes, or more generally arbitrary closed image sets, to be geometrically "consistent". We present this notion as a natural generalization of traditional multi-view geometry, which deals with consistency for points. After discussing some general results, we present a "dual" formulation for consistency, that gives conditions for a family of planar sets to be sections of the same object. Finally, we introduce a more general notion of silhouette "compatibility" under partial knowledge of the camera projections, and point out some possible directions for future research.

Introduction

When can a set of 2D silhouettes be projections of the same object in space? This seemingly simple question is related to a variety of practical problems in computer vision, such as multi-view segmentation [10], object recognition [START_REF] Zitnick | Object instance recognition using triplets of feature symbols[END_REF], and multi-view stereo [START_REF] Furukawa | Carved visual hulls for imagebased modeling[END_REF]. Geometric consistency is sometimes taken for granted, when appearance-based features give reasonable evidence that the silhouettes are associated to the same object. However, it is clear that incorporating geometric constraints can be important, either in the process of inference, or for correcting the effects of noisy data.

In this paper, we analyze the notion of "consistency" for silhouettes and for more general closed image sets (Figure 1). We consider opaque objects projected in different images and assume, in our initial setting, the knowledge of all camera parameters. The theory can be seen as a natural extension of classical multi-view geometry, which provides conditions for points (and sometimes lines) to be consistent (i.e., to correspond) in terms of given camera projections [START_REF] Faugeras | The geometry of multiple images: the laws that govern the formation of multiple images of a scene and some of their applications[END_REF][START_REF] Hartley | Multiple view geometry in computer vision[END_REF]. This also expands an analogy initiated by the generalized epipolar constraint introduced in [START_REF] Åström | Generalised epipolar constraints[END_REF]. Throughout our discussion, we point out several results that hold for arbitrary closed image sets (or sometimes convex image sets), that have identical counterparts in the theory of point correspondences.

As consistency is clearly not a metric property of silhouettes, it is natural to adopt the framework of projective geometry [START_REF] Coxeter | Projective geometry[END_REF], since this eliminates various degenerate situations and allows, for example, to unify the cases of orthographic and perspective projections. This is also typical in multi-view geometry. 1 Another advantage of the projective language is that it provides a homogeneous formulation of duality. In particular, we exploit the fact that perspective projections are related to planar sections in the dual space [27], to define a very natural "dual" notion of consistency, expressing conditions for a family of planar sets to be sections of the same object.

Finally, we also consider the case of having only partial knowledge of the camera parameters, and discuss a more general concept of "compatibility" for silhouettes. This extends a setting first considered in [START_REF] Bottino | Introducing a new problem: Shape-from-silhouette when the relative positions of the viewpoints is unknown. Pattern Analysis and Machine Intelligence[END_REF].

The overall goal of the paper is not to give the "final" answer to the complex problem of silhouette consistency, but rather to make a first formal foray in that field, with precise definitions that have been missing so far in a general setting, spelling out rigorously what is known in this area and adding a set of new results. For example, we do not deal with algorithmic issues here, and we assume throughout the presentation an ideal setting with no noise. We believe this to be a necessary first step, much in the same way as multi-view geometry initially characterizes exact point correspondences, and then makes use of the theory to infer camera parameters from real world data.

Previous work.

The most widespread application of silhouette consistency has been for designing alternatives to point-based methods for camera calibration, required for dealing with smooth and textureless surfaces. Indeed, ever since the seminal work of Rieger [START_REF] Rieger | Three-dimensional motion from fixed points of a deforming profile curve[END_REF], the problem of estimating camera motion or calibration parameters using only silhouettes has received considerable attention: see, e.g., [4,[START_REF] Furukawa | Structure and motion from images of smooth textureless objects[END_REF][START_REF] Hernández | Silhouette coherence for camera calibration under circular motion[END_REF][START_REF] Joshi | Structure and motion estimation from dynamic silhouettes under perspective projection[END_REF][START_REF] Mendonc ¸a | Epipolar geometry from profiles under circular motion[END_REF][START_REF] Wong | Reconstruction of sculpture from its profiles with unknown camera positions[END_REF]. Albeit with some variations, all these methods exploit (more or less directly) the geometric constraints provided by the epipolar tangencies [START_REF] Åström | Generalised epipolar constraints[END_REF].

In addition to camera calibration, silhouette consistency has been enforced explicitly for other tasks; for example multi-view segmentation [START_REF] Campbell | Automatic 3d object segmentation in multiple views using volumetric graph-cuts[END_REF], or 3D-reconstruction [START_REF] Cremers | Multiview stereo and silhouette consistency via convex functionals over convex domains[END_REF] and recognition [START_REF] Lazebnik | On pencils of tangent planes and the recognition of smooth 3d shapes from silhouettes[END_REF]. Another interesting "artistic" application is discussed in [26].

There exists limited theoretical work on silhouette consistency, and it is always restricted to special situations. In particular, the problem of determining whether a family of silhouettes can correspond to a real object is considered in [START_REF] Bottino | Introducing a new problem: Shape-from-silhouette when the relative positions of the viewpoints is unknown. Pattern Analysis and Machine Intelligence[END_REF][START_REF] Jacobs | Judging whether multiple silhouettes can come from the same object[END_REF][START_REF] Yi | Efficiently determining silhouette consistency[END_REF], but results are only given for the case of orthographic projections and somewhat restricted camera motion. Some theoretical facts, e.g., the fact that epipolar tangency conditions do not imply global consistency, can be found in [4,6].

Finally, the duality between projections and planar sections is well known in both the euclidean (orthographic) setting [START_REF] Van Hove | Silhouette-slice theorems[END_REF]30], as well as in the projective (perspective) case [27,[START_REF] Tevelev | Projective duality and homogeneous spaces[END_REF]. However, we are not aware of work on consistency for planar sections.

Main contributions.

• We formally introduce a notion of geometric consistency for arbitrary closed sets in the case of general projective (perspective) cameras. We present some new results (Propositions 1, 3), and collect others which are scattered in more applied work (Propositions 2, 4).

• We investigate in detail the relationship between pairwise consistency, epipolar tangencies, and our more general notion of consistency (Sections 2.2 and 2.3). We also observe the existence of special "tangential points", which play an important role in constraining the consistency for silhouettes.

• We restate the notion of consistency in terms of duality, expressing the condition for planar sets to be sections of the same object. For convex silhouettes, we show that the dual of the visual hull coincides with the convex hull of the dual image of the silhouettes (Proposition 6).

• We define the notion of compatibility for silhouettes, which characterizes silhouettes that may be geometrically consistent for appropriate camera parameters. This setting generalizes a viewpoint first introduced in [START_REF] Bottino | Introducing a new problem: Shape-from-silhouette when the relative positions of the viewpoints is unknown. Pattern Analysis and Machine Intelligence[END_REF].

Notation. Our analysis will be coordinate free, so we consider projective cameras as linear maps M : P 3 \{c} → P 2 , where c is the camera pinhole or center. The action of a camera will be indicated with M(p) = u. For any camera M and set of points T in an image, we define the associated visual cone as M -1 (T ), where M -1 denotes the pre-image set (this definition excludes the pinhole from the visual cone, however this will be irrelevant for our discussion).

Consistency of image sets

In this section, we introduce a notion of geometric consistency for arbitrary closed sets in different images. Our definition is very natural, and similar concepts have previously been used to introduce "incoherence" measures for silhouettes [4,[START_REF] Hernández | Silhouette coherence for camera calibration under circular motion[END_REF]. Compared to these works, we focus on analyzing some theoretical properties of consistency, rather than finding strategies for putting it into practice. After a general discussion, we consider the important case of two silhouettes (Section 2.2), pointing out how the consistency condition is basically equivalent to the popular "epipolar tangency" constraint, but only applied to extremal tangents. We then investigate in more detail the difference between pairwise consistency, and general consistency, and argue the importance of special "tangential points", which constrain consistent silhouettes (Section 2.3).

Basic definitions

Let M 1 , . . . , M n be n perspective cameras with distinct centers c 1 , . . . , c n , and let T 1 , . . . , T n be a family of closed sets, one in each image. For example, the sets T i could be a finite collection of points, curves or closed regions. For each i = 1, . . . , n, we let C i = M -1 i (T i ) be the visual cone associated to T i . Definition 1. The sets T 1 , . . . , T n are said to be consistent (relative to the cameras M 1 , . . . , M n ) if there exists a nonempty set R ⊆ P 3 \ {c 1 , . . . , c n } such that M i (R) = T i for all i = 1, . . . , n.

When T 1 , . . . , T n are consistent for M 1 , . . . , M n , the visual hull associated with T 1 , . . . , T n is given by H = i C i , and it is the largest set that projects onto T 1 , . . . , T n .

If all T 1 , . . . , T n are singletons, then consistency reduces to the classical notion of point correspondence [START_REF] Faugeras | The geometry of multiple images: the laws that govern the formation of multiple images of a scene and some of their applications[END_REF][START_REF] Hartley | Multiple view geometry in computer vision[END_REF]; in this case, the visual hull is simply the triangulated 3D-point. Extending this analogy, consistent image sets can be seen as n-tuples of projections of all possible objects in space. Moreover, it is clear that there is a one-to-one correspondence between n-tuples of consistent sets, and the (exact) visual hulls associated with a fixed set of n cameras.

In principle, the concept of visual hull is not well defined if the original silhouettes (or image sets) are not geometrically consistent. This is rarely taken into consideration, and it is customary to define the visual hull simply as the intersection of the cones i C i for arbitrary (non necessarily consistent) silhouettes: this operation can be justified by noting that if i C i is not empty, then it is the visual hull associated with the subsets Ti = M i ( j C j ) ⊆ T i , which will always be consistent. In fact, consistency is clearly equivalent to the fact that Ti = M i ( j C j ) = T i for all i = 1, . . . , n, or to T i ⊆ M i ( j C j ), since the opposite inclusion is always true. We collect a few other simple but useful properties: Proposition 1. Let T 1 , . . . , T n be arbitrary closed image sets.

1. T 1 , . . . , T n are consistent if and only if for each i = 1, . . . , n, and for all u i ∈ T i , the visual ray M -1 (u i ) intersects j =i C j .

2. T 1 , . . . , T n are consistent if and only if

T i ⊆ M i   j =i C j   , ∀i ∈ {1, . . . , n}. (1)
See Figure 2.

3. If T 1 , . . . , T n are consistent, then any subfamily T i1 , . . . , T is is consistent (for the associated cameras M i1 , . . . , M is ).

Proof. The first property follows from the fact that

T i ⊆ M i ( j C j ) can be expressed as M -1 i (u i ) ∩ j C j = ∅ for all u i ∈ T i , which in turn is equivalent to M -1 i (u i ) ∩ j =i C j = ∅, since M -1 i (u i ) ⊆ C i .
The second and third properties are consequences of the first one.

This might be a good moment to point out that the notion of geometric consistency discussed in this paper is somewhat independent from a more intuitive (but less formal) concept of "similarity" of appearance. Indeed, consistent silhouettes may actually look completely different (as in the example in Figure 1), while, on the other hand, almost identical silhouettes may be geometrically inconsistent. Thus, the concept might be well suited for being used alongside more traditional feature-based methods for recognition.

Finally, for the rest of the paper we will make the following assumption for all n-tuples of closed image sets T 1 , . . . , T n and cameras M 1 , . . . , M n :

(A) For each camera center c i and visual cone C j , with i = j, c i does not belong to C j . This condition is useful for excluding various degenerate situations; for example, it guarantees that the visual hulls associated with all subfamilies of T 1 , . . . , T n are closed sets. 2

Pairwise consistency

Let us assume that we are given only two image sets T 1 , T 2 and two cameras M 1 , M 2 . According to Proposition 1, we know that T 1 , T 2 are consistent for M 1 , M 2 if and only if

T 1 ⊆ M 1 (C 2 ) , and T 2 ⊆ M 2 (C 1 ) (2)
In other words, we require for each set to be contained in the projection of the visual cone associated to the other one.

Pairwise consistency is closely related to the popular epipolar tangency constraint [START_REF] Åström | Generalised epipolar constraints[END_REF][START_REF] Wong | Reconstruction of sculpture from its profiles with unknown camera positions[END_REF]. Indeed, we can restate the condition (2) in terms of epipolar geometry as follows.

Proposition 2. Two image sets T 1 , T 2 are consistent if and only if the set of epipolar lines in the first image which intersect T 1 is in epipolar correspondence with the set of epipolar lines in the second image which intersect T 2 .

Proof. The statement can be seen as a consequence of the first property in Proposition 1. In fact, the epipolar correspondence condition guarantees that for all i = 1, 2, and for every point u i ∈ T i there exists at least one corresponding point u j ∈ T j (j = i), so that triangulating all pairs of associated points (i.e., intersecting the cones C 1 , C 2 ) we obtain a set R ⊆ P 3 \ {c 1 , c 2 } that projects exactly onto T 1 and T 2 . Note that assumption (A) is equivalent to the fact that, in each image, the epipole lies outside of the given set. This result can also be stated in terms of visual cones: if we consider the pencil of planes passing through the two camera centers, then pairwise consistency is equivalent to the fact that each visual cone meets the same set of planes in the pencil.

If we assume that T 1 , T 2 are connected regions bounded by smooth curves, then pairwise consistency basically reduces to the fact that extremal epipolar tangents (i.e., "outermost" epipolar lines that are tangent to the contours) are in epipolar correspondence [4] (Figure 3). On the other hand, pairwise consistency does not require non-extremal epipolar tangents to be matched. In fact, perhaps somewhat unintuitively, the visual hull generated by two silhouettes with non-extremal epipolar tangencies (matched or unmatched) will always present a complicated topology because it must take into account for all possible self-occlusions (Figure 4). The difference between the two examples shown is just that in the case (top) of unmatched internal epipolar tangencies, any object which is consistent with the silhouettes must project with occlusions (i.e., the inconsistent "branches" are necessarily projections of different parts of the observed object); on the other hand, in the case (bottom) where internal epipolar tangencies are matched, there exists an object (different from the visual hull) that has a simple topology and that projects onto the silhouettes without occlusions.

Returning to the case of an arbitrary number of sets, we will say that T 1 , . . . , T n are pairwise consistent if each pair T i , T j , i = j is consistent. From Proposition 1, we know that this holds whenever T 1 , . . . , T n are consistent. The converse is not true, as pointed out in [4, 6] (see also Figure 6). However, it would be useful to clarify the practical distinction between these two notions: much of our discussion in the following will be aimed at a better understanding of this issue. For example, in the case of three sets, pairwise consistency often implies that there is at least an "approximate" consistency: this property directly generalizes the fact that three non-coplanar visual rays that converge pairwise will always converge [START_REF] Hartley | Multiple view geometry in computer vision[END_REF]. not intersect the plane spanned by c 1 , c 2 , c 3 , then i C i is not empty.

Proof. It is sufficient to prove that, say,

T 1 ∩M 1 (C 2 ∩C 3 ) is not empty. Let R = M 1 (C 2 ) ∩ M 1 (C 3
). The assumptions on the centers guarantee that R is a quadrilateral (it is the intersection of two connected projected cones; see Figure 5). From pairwise consistency (2), we know that T 1 is tightly "inscribed" in R, meaning that T 1 ⊆ R and T 1 intersects of the four edges of R. The same holds for

M 1 (C 2 ∩ C 3 ): indeed M 1 (C 2 ∩ C 3 ) ⊆ M 1 (C 2 ) ∩ M 1 (C 3 ) holds

From pairwise to general consistency

We have already noted that pairwise consistency does not imply general consistency [4, 6]: in fact, as shown in the example in Figure 6, the two notions are not equivalent even in Figure 6: Using orthogonal orthographic projections (that may be thought as "side", "top", and "front" views), we observe that two appropriately "clipped" squares and one full square are pairwise consistent but not globally consistent: indeed, the reprojection of the visual hull associated to the first two silhouettes into the third view does not contain the corresponding silhouette (cf. Proposition 1), since the top right corner is missing. In fact, the volume intersection of the three visual cones (shown on the right) does not project onto the third silhouette.

the restricted case of convex image sets (contrary to a claim in [4]). Interestingly, however, we can show that for convex sets (i.e., sets that are convex in an affine chart: we refer to the next section for a discussion on convexity in the projective setting) general consistency is actually guaranteed by "triplet-wise" consistency. A similar fact has been pointed out in [START_REF] Yi | Efficiently determining silhouette consistency[END_REF].

Proposition 4. Let T 1 , . . . , T n be closed convex sets. If T i , T j , T k are consistent for every {i, j, k} ⊆ {1, . . . , n} then T 1 , . . . , T n are consistent.

Proof. From Proposition 1, it is sufficient to prove that for every i = 1, . . . , n, and for all u i ∈ T i , we have M -1 (u i ) ∩ j =i C j = ∅. Because of convexity, the sets U i j = M -1 i (u i ) ∩ C j are intervals; moreover, they intersect pairwise because of the assumption of "triplet-wise" consistency: this implies that they all intersect. 3 This statement closely resembles "Helly-type" theorems in computational geometry [START_REF] Toth | Handbook of discrete and computational geometry[END_REF], and is a generalization of the fact that point correspondence is always implied by triplet-wise point correspondence. An inspection of the proof also shows that it is not actually necessary for the sets T i to be convex, but only for their intersections with epipolar lines (rather than all lines) to be intervals.

Returning to the case of arbitrary image sets, we see from (2) that pairwise consistency is equivalent to

T i ⊆ j =i M i (C j ), i = 1, . . . , n. (3)
Moreover, if the sets T i are consistent there is a chain of inclusions T i ⊆ M i j =i C j ⊆ j =i M i (C j ). In fact, as pointed out in [4], pairwise and general consistency can 3 To be precise, it should be noted that, although in projective space, none of the cones C j contain c i (because of our assumption (A)), so we may treat the sets U i j as intervals on a real line. be seen as the same formal condition except for the inverted order of the projection and intersection of the visual cones.

It is also natural to ask ourselves how the consistency condition actually constrains the shape of the image sets. In this regard, we observe that if T 1 , . . . , T n are consistent for M 1 , . . . , M n , then a set, say T 1 , may be replaced by any strictly larger one T1 without affecting consistency, provided that T1 ⊆ M 1 j =1 C j , because this guarantees that condition (1) remains satisfied in each image. For example, as shown in Figure 8, we are allowed to modify the shape of a set T i with a small protrusion, in the neighborhood of any point that does not belong to the boundary of M i j =i C j . On the other hands, points belonging to the intersection of T i with the boundary of M i j =i C j are more constrained: we will refer to these points as "tangential points", since the associated visual rays are "tangent" to the visual hull associated with the remaining silhouettes j =i C j . Tangential points play a crucial role in characterizing consistency, and we intend to investigate their properties in more detail in the future. For the moment, we simply point out that the visual ray corresponding to a tangential point will typically graze the surface of one of the remaining visual cones, so that the ray will correspond to an epipolar tangent (not necessarily extremal) for the associated image (Figure 9, left). In general, it is also possible for the ray to intersect the set j =i C j only at an "intersection curve" [23], so that it is not tangent to any cone, and gives rise to a "tangential triple point" on the visual hull surface (Figure 9, right): however, it is not hard to realize that this case will not occur for generic projections of smooth solids, essentially because points on a smooth surface are never visible from three generic viewing directions (one also has to note tangential triple points are "hard points" [22], i.e., points of the visual hull that are guaranteed to belong to the actual object being observed). that the new silhouette is still contained in the reprojection of the visuall hull associated to the other views, consistency is preserved. Note that the reprojection onto the third view is not the same as in the previous case, however condition (1) remains satisfied since the new reprojected set is larger.

A dual view of consistency

In this section we revisit the notion of consistency from the viewpoint of duality. In particular, we discuss the relationship between consistency for projections, as presented in the previous section, and a different notion of consistency for planar sections of solids. We focus in our presentation on the case of convex sets, since the duality is much simpler in this setting. We defer a more general analysis to future work.

Duality

The basis of many similar notions of duality is the fact that points and hyperplanes in some n-dimensional space can play symmetric roles. For example, in R n any hyperplane through the origin can be represented with its orthogonal vector. In projective space P n , all hyperplanes correspond to points of the dual projective space (P n ) * : in fact, any k-dimensional linear subspace H in P n , can be associated with an (n -k)-dimensional linear subspace H * in (P n ) * , and this duality interchanges the role of "join" and "meet" [START_REF] Coxeter | Projective geometry[END_REF].

This identification of hyperplanes in P n with points in (P n ) * can be used to introduce notions of "duality" for more general objects in P n . For example, if S is a smooth or algebraic hypersurface in P n , (the closure of) the set of tangent hyperplanes at points of S forms a dual hypersurface in (P n ) * [START_REF] Tevelev | Projective duality and homogeneous spaces[END_REF]. However, dual hypersurfaces will typically have self-intersections (e.g., in the case of curves, crossings correspond to bitangents of the original curve). For this reason, we will consider instead duality for (properly) convex sets in projective space [START_REF] Goldman | Projective geometry on manifolds[END_REF]. In the following, we will say that a set K in P n is convex when there exists an affine chart (i.e., an appropriate dehomogenization of projective coordinates) so that K is compact and convex in the usual affine sense. This is equivalent to a more analytical definition, for which K is convex if it is the projectivization of a cone K in R n+1 that is convex, closed, and pointed (i.e., it does not contain any line) [START_REF] Boyd | Convex optimization[END_REF]. To any convex set K in P n , we can now associate a dual set K • in (P n ) * , that can be Right: a tangential point may also be the projection of a "tangential triple point": in this case, the associated viewing ray is not tangent to any visual cone. characterized geometrically as the closure of the set of all hyperplanes H such that H ∩ K is empty. Alternatively, the same set K • is the projectivization of the polar cone ( K) • of a convex cone K in R n+1 that corresponds to K.

For an appropriate choice of coordinates, convex duality in projective space essentially coincides with the usual concept of polarity for affine sets [START_REF] Boyd | Convex optimization[END_REF]. However, while many properties of affine polarity only apply for convex sets containing the origin (e.g., biduality: K •• = K), the projective framework is completely "homogeneous", and does not require similar conditions. Moreover, this setting is useful for dealing with general perspective projections.

All the results in the next section can be understood geometrically (and proven using purely synthetic arguments), by exploiting the characterization of duality as the set of "complementary" hyperplanes. More formal analytical proofs can be given using cones in R n+1 : see the supplementary material for details.

Duality and visual hulls

A perspective projection M : P 3 \{c} → P 2 with center c also defines a dual map M * : (P 2 ) * → (P 3 ) * that associates lines in P 2 with planes in P 3 through c. In particular, the image of M * is a plane in the dual space (P 3 ) * , namely the dual of the center c * (i.e., the set of planes containing c). The following important result holds: Proposition 5. Let M be a perspective projection with center c, and let K ⊆ P 3 \ {c} be a convex set. Then

M(K) = L is equivalent to M * (L • ) = K • ∩ c * . (4)
Indeed, the geometric intuition for condition (4) is that a line does not meet L if and only if its preimage does not meet K. This result basically states fact that projecting a convex object is equivalent to taking a planar section of its dual. We also refer to [27], where a similar idea is used in a "smooth" setting to investigate the qualitative relationship between image contours and projective shapes.

Motivated by Proposition 5, we introduce a notion of "sectional consistency". We assume that we are given projective maps P 1 , . . . , P n where each P i : P 2 → P 3 identifies P 2 with a plane π i in P 3 . The planes π 1 , . . . , π n are distinct. Definition 2. A family L 1 , . . . , L n of sets in P 2 is said to be sectionally consistent (relative to the embeddings P 1 , . . . , P n ) if there exists K ⊆ P 3 such that K ∩ π i = P i (L i ) for all i = 1, . . . , n. This notion of sectional consistency is "dual" to geometric consistency from Definition 1. Indeed: Proposition 6. A family L 1 , . . . , L n of convex sets in P 2 is consistent for a set of projections M 1 , .

= K • , where K is the convex hull of M * 1 (L • 1 ), . . . , M * n (L • n ). 4
Proof. The first claim is a direct consequence of (4). For the second part, we can use the fact that for arbitrary closed convex cones K 1 , . . . , K n in R N :

[Conv(K 1 ∪ . . . ∪ K n )] • = K • 1 ∩ . . . ∩ K • n . (5)
See for example [30]. In fact, it is sufficient to take

K i = M * i (L • i ), and observe that (M * i (L • i ))
• is in fact the visual cone C i (see the supplementary material for details). Consistency is arguably more intuitive in its dual formulation, since planar sets in space are conceptually easier to grasp than families of cones. For example, the pairwise consistency constraint for two convex silhouettes T 1 , T 2 dualizes to the fact that M * 1 (T • 1 ) and M * 2 (T • 2 ) have the same intersection with the "dual baseline" (c 1 ∧ c 2 ) * (here ∧ denotes the "join" of two points): see Figure 10.

We also point out a duality between the "combinatorial" structure of the visual hull H and of its dual K = H • . Indeed, by interpreting the supporting planes of H as points on the boundary of K, we see that a plane π supports a cone (patch) C i if and and only if, in the dual space, it represents a point belonging to the planar curve that is the boundary of M * i (L • i ); on the other hand, if a supporting plane π meets the visual hull at the intersection of two cones C i , C j (i.e., at an intersection curve [23]), then it is dual to a point on a ruled patch, joining the two boundaries of M * i (L • i ) and M * j (L • j ); finally, if π supports H at the intersection of three cones C i , C j , C k (i.e., a triple point), then it represents a point of the dual hull K that belongs to a planar patch, spanning across the three boundaries of

M * i (L • i ), M * j (L • j ), M * k (L • k )
(note that ruled and planar patches are typical for convex hulls of curves in space [START_REF] Koenderink | Solid shape[END_REF]). See Figure 11.

The notion of sectional consistency given in Definition 2 is reminiscent of questions in geometric tomography [START_REF] Gardner | Geometric tomography[END_REF], or stereology [START_REF] Howard | Unbiased stereology: threedimensional measurement in microscopy[END_REF]. In tomography, for example, the duality between projections and sections is well studied, but typically in an affine setting that considers only orthographic projections. Nevertheless, it is quite possible that tools from these related fields could provide interesting new insight for problems in computer vision.

Compatible silhouettes

Throughout the paper, we have always considered families of silhouettes (or image sets) T 1 , . . . , T n together with known camera projections M 1 , . . . , M n . However, we can introduce a more general notion of "geometric compatibility", that can be applied in the case of incomplete knowl-edge of the camera parameters. More precisely, we define an arbitrary family of silhouettes (or image sets) T 1 , . . . , T n to be compatible with some partial knowledge P of the camera parameters when there exist projections M 1 , . . . , M n that agree with P and for which T 1 , . . . , T n are consistent according to Definition 1.

Simply put, this notion characterizes the most general condition for a family of silhouettes (or sets) to be feasible projections of a single object. The study of similar issues was initiated in [START_REF] Bottino | Introducing a new problem: Shape-from-silhouette when the relative positions of the viewpoints is unknown. Pattern Analysis and Machine Intelligence[END_REF], where the authors assume the external parameters of the cameras to be unknown, and analyze in detail the geometric constraints for compatibility in a particular case (orthographic viewing directions parallel to the same plane).

The first problem raised by the definition of compatibility is to understand whether it actually constrains the silhouettes and the camera parameters, or if instead it is just a vacuous condition. For example, in [START_REF] Bottino | Introducing a new problem: Shape-from-silhouette when the relative positions of the viewpoints is unknown. Pattern Analysis and Machine Intelligence[END_REF] it is pointed out that if one considers internal as well as external parameters to be unknown, then any family of silhouettes will be compatible (they note that considering a convex object, and applying local protrusions with appropriate shapes, one is able to produce arbitrary silhouettes by placing cameras near the surface). This interesting observation, however, would seem to imply that SfM methods can never exploit the geometry of the silhouettes in order to recover camera parameters, since a particular family of silhouettes would provide no information on the viewing conditions. On the other hand, we note that the construction proposed in [START_REF] Bottino | Introducing a new problem: Shape-from-silhouette when the relative positions of the viewpoints is unknown. Pattern Analysis and Machine Intelligence[END_REF] violates our assumption (A) for all pairs of cameras and silhouettes (viewing cones must be extremely "wide", and thus contain all other centers which lie near the surface of the convex object). We argue instead that by considering only certain regions in the space of parameters, the compatibility of silhouettes does provide constraints on the camera projections, even without knowledge of internal parameters. For example, consider two pictures with two arbitrary silhouettes (closed connected regions with non empty interior) T 1 , T 2 , one in each image: in this case, it is easy to realize that the fundamental matrices for which T 1 , T 2 are consistent and such that each epipole lies outside the corresponding silhouette (as in condition (A)), form a space that always has 5 degrees of freedom, compared to 7 in the unconstrained case. Indeed, for any choice of epipoles e 1 / ∈ T 1 and e 2 / ∈ T 2 , there is one-parameter family of feasible fundamental matrices, corresponding to the homographies between the pencil of lines through e 1 and e 2 that associate the projections of T 1 and T 2 (as segments of the pencil). On the other hand, assuming e 1 ∈ Int(T 1 ) and e 2 ∈ Int(T 2 ), the space of fundamental matrices for which T 1 , T 2 are consistent has 7 degrees of freedom, since the homography of epipolar lines is unconstrained. The fact that a pair of silhouettes imposes (with certain assumptions) two conditions on the epipolar geometry, justifies methods that have exploited two silhouettes to recover camera parameters when the ambiguity of camera motion had only two degrees of freedom, namely in the case of pure translation motion [START_REF] Cross | Parallax geometry of smooth surfaces in multiple views[END_REF], or when a visible plane provides a homography between the views [9].

In our opinion, a better understanding of the theory of "silhouette compatibility" can be important for several practical reasons, such as spelling out conditions for when silhouettes may or may not be used to determine camera geometry (and possibly help design better algorithms), or similarly to give conditions for a family of silhouettes alone to provide a unique representation of (the visual hull of) an object. Duality might also prove to be a useful tool for investigating these kinds of questions: much in the same way as in Section 3, we realize that the notion of compatibility of silhouettes can be expressed in terms of the compatibility of planar regions, which need to be "assembled" consistently in order to be feasible sections of a single object.

Conclusions

We have analyzed in detail the notion of "geometric consistency" for arbitrary image sets, in a setting that can be seen as an extension of traditional multi-view geometry. In the case of convex silhouettes we have also discussed a "dual" interpretation of consistency, expressing conditions for planar sets to be sections of a single object. These concepts lead to a more general notion of silhouette "compatibility", that does not require (complete) knowledge of the camera parameters.

We plan to extend this work in various directions. On the practical side, our results need to be revisited for dealing with noisy data, and the theory may be used for comparing different measures of "inconsistency" such as the ones considered in [4, [START_REF] Hernández | Silhouette coherence for camera calibration under circular motion[END_REF][START_REF] Mendonc ¸a | Epipolar geometry from profiles under circular motion[END_REF]. On the theoretical side, in addition to the questions on "compatibility" discussed above, there remains to gain a better understanding of dual consistency for non-convex silhouettes. In geometric tomography, for example, it is typical to study sections of convex bodies [START_REF] Gardner | Geometric tomography[END_REF]; however, for applications in vision, this assumption is probably restrictive. Solutions to all of these problems would be useful in many practical settings, and would help us clarify the fundamental relationship between two-dimensional projections, and the natural concept of "shape".

Figure 1 :

 1 Figure 1: Geometrically consistent silhouettes are feasible projections of a single object.

Figure 2 :

 2 Figure 2: The silhouettes from Figure 1 are geometrically consistent (relative to three orthogonal orthographic projections), since each silhouette is contained in the reprojection (in gray) of the intersected visual cones associated to the other views. See Proposition 1.

Figure 3 :

 3 Figure 3: Parwise consistency: the sets of epipolar lines intersecting each silhouette must be in epipolar correspondence. In particular, the extremal lines lines A 1 , A 2 and B 1 , B 2 must correspond.

Proposition 3 .

 3 Figure4: The visual hull generated by two "hook-shaped" silhouettes always presents self-occlusions, whether the internal epipolar tangents are unmatched (top) or matched (bottom). In these examples we considered orthogonal orthographic projections, but the behavior is completely general.

Figure 5 :

 5 Figure 5: Proof of Proposition 3. The set T1 (dark gray) and M1(C2 ∩C3) are both "inscribed" in the quadrilateral M1(C2)∩ M1(C3) (light gray), and must thus intersect.

  , and M 1 (C 2 ∩ C 3 ) intersects all edges of R since C 2 , C 3 are themselves pairwise consistent (an extremal epipolar line in M 1 (C 2 ) is the projection of a line in C 2 which must intersect C 2 ∩ C 3 ). The claim follows from continuity arguments: for example, consider paths in T 1 and M 1 (C 2 ∩ C 3 ) connecting different pairs of opposite edges.

Figure 7 :

 7 Figure 7: Proof of Proposition 4: the visual ray corresponding to p meets the other visual cones in intervals that intersect pairwise, and thus must all intersect.

Figure 8 :

 8 Figure 8: If we modify the second silhouette from Example 1 so

Figure 9 :

 9 Figure 9: Tangential points are the intersections of a silhouette Ti with the boundary of Mi j =i Cj . Left: the silhouette of a torus and the reprojection of the visual hull generated by three different views. The six tangential points are epipolar tangencies (colors indicate projections of cone surfaces and the associated epipolar tangency points).Right: a tangential point may also be the projection of a "tangential triple point": in this case, the associated viewing ray is not tangent to any visual cone.

Figure 10 :

 10 Figure 10: Dual-pairwise-consistency. The dual images of the silhouettes M * 1 (T • 1 ) and M * 2 (T • 2 ) are sectionally consistent if and only if they have the same intersection with the "dual baseline" (c1 ∧ c2) * .

Figure 11 :

 11 Figure 11: Structure of the "dual visual hull". The boundary of the dual visual hull is composed of planar curves, ruled surface patches and planar triangular patches: each of these are associated to specific components of the visual hull. See text for details.

  . . , M n if and only if L • 1 , . . . , L • n are sectionally consistent for the embeddings M * 1 , . . . , M * n . Moreover, if consistency holds, and H is the visual hull associated with L 1 , . . . L n , then H

The most significant difference between the euclidean and projective frameworks in our setting is that in the latter case visual cones are twosided. An alternative (but perhaps less natural) approach would have been the use of oriented projective geometry[START_REF] Stolfi | Oriented projective geometry: A framework for geometric computations[END_REF].

This condition is not actually necessary for all of our results, and weaker assumptions may often be considered. However, for the sake of simplicity, we give a single condition that is valid throughout the paper.

The convex hull of a closed set contained in an affine chart in P n is independent of the choice of projective coordinates only when the set is connected: this is true forM * 1 (L • 1 ) ∪ . . . ∪ M * n (L • n ) since assumption (A) guarantees that M * i (L • i ) ∪ M * j (L • j )is connected for all i, j.

Acknowledgments. This work was supported in part by the ERC grant VideoWorld, the Institut Universitaire de France, the Inria -CMU associate team GAYA, and ONR MURI N000141010934.


Chen Liang

B Mourrain

J P Pavone

Subdivision Methods for the Topology of 2d and 3d Implicit Curves

In this paper, we describe a subdivision method for handling algebraic implicit curves in 2d and 3d. We use the representation of polynomials in the Bernstein basis associated with a given box, to check if the topology of the curve is determined inside this box, from its points on the border of the box. Subdivision solvers are used for computing these points on the faces of the box, and segments joining these points are deduced to get a graph isotopic to the curve. Using envelop of polynomials, we show how this method allow to handle efficiently and accurately implicit curves with large coefficients. We report on implementation aspects and experimentations on 2d curves such as ridge curves or self intersection curves of parameterized surfaces, and on silhouette curves of implicit surfaces, showing the interesting practical behavior of this approach.

Introduction

In this paper, we address the problem of computing the topology of 3D curves resulting from the intersection of two algebraic surfaces. Algebraic curves and surfaces are compact representations of shapes, which can be complex and have numerous advantages over parametric ones, such as easy determination of inside/outside of the surface. This is particularly useful when we have to apply logical operations (union, subtraction, etc.) between two solid objects, defined implicitly. In such problems, computing the intersection of two surfaces is a critical operation, which has to be performed efficiently and accurately. Implicit curves and surfaces have also disadvantages such as difficulty in performing graphical display, but the method that we propose in this paper is step towards handling such problems, since it allows fast display of this implicitly 2d and 3d curves. On the other hand, dealing with parameterized surfaces naturally leads to the computation of implicit curves. Let us mention in particular, the computation of the intersection curve of two surfaces, self-intersection curves, plane sections and ridge curves (which are defined implicitly on these parameterised surfaces, though they are usually approximated by parameterised curves). Such problems reduce to the analysis of a curve defined by n-1 polynomial equations, in a space of dimension n (here n = 2, 3, 4).

One major obstacle for adopting implicit representations instead of parametric representations concerns the piecewise linear approximation of such curves or surfaces for visualization purposes. A brute force approach would be an exhaustive evaluation for approximating the zero level set, which is obviously very inefficient. A typical alternative scenario is to adopt a divide-and-conquer approach. Larger undetermined domains are broken down to smaller predictable domains in which the topological feature and eventually, the curve/surface itself can be inferred efficiently. An objective of this paper is to describe an efficient method, which allows us to capture the topology of an implicit curve, when this curve is smooth 1 , but also to localize the singular points if they exist.

The problem of computing the topology of curves has been approach in different ways. A first family of methods is based on a sweeping approach. For 2D planar algebraic curves, such approach has been studied in [START_REF] Grandine | A new approach to the surface intersection problem[END_REF] and [START_REF] González-Vega | Efficient topology determination of implicitly defined algebraic plane curves[END_REF]. It was later extended by Gatellier et al. in [START_REF] Gatellier | Computing the topology of 3-dimensional algebraic curves[END_REF] to the 3D spatial curves resulting from the intersection of two algebraic surfaces. See also [START_REF] Morken | Computing the topology of real algebraic space curves[END_REF]. These methods use a conceptual sweeping line/plane perpendicular to some projection axis, and detect the critical topological events, such as tangents to the sweeping planes and singularities. The final output of these methods are a graph of connected vertices complying to the topology of the original curves. A notable problem of aforementioned approaches is that they relies of the computation of sub-resultant sequences, which can be a bottleneck in many examples with large degree and large coefficients (see Section 4.1).

Another family of methods are the subdivision based techniques, which uses a simple criterion to remove domains which do not contain the roots. A crucial problem involved here is how to efficiently and reliably deduce the root information in a given interval (or a bounding box). In these methods, instead of using monomial representation, we represent the equations using Bernstein basis [START_REF] Farin | Curves and Surfaces for Computer Aided Geometric Design: A Practical Guide[END_REF]. Among early attempts, Sederberg [START_REF] Sederberg | Algorithm for algebraic curve intersection[END_REF] converted an algebraic curve in to piecewise triangular Bernstein basis. See also [START_REF] Keyser | Efficient and exact manipulation of algebraic points and curves[END_REF] combining symbolic and numeric techniques to compute the topology of 2D curves. The approach of [HFH + 05] for computing the curves of intersection of two parameterised surfaces is also combining subdivision techniques with regularity criterion, exploiting the properties of the intersection curve in the 2D parameter domains.

The first problem of computing roots of univariate polynomials has been analyzed for instance in [START_REF] Mourrain | Bernstein's basis and real root isolation[END_REF], where root information tests are by based on Descartes' Law of Sign and its variant in the Bernstein basis. This approach has been extended to the approximation of isolated roots of multivariate systems. In [START_REF] Sherbrooke | Computation of the solutions of nonlinear polynomial systems[END_REF], the author used tensor product version of Bernstein basis and integrated domain reduction techniques to speed up the convergence and reduce the number of subdivisions. In [START_REF] Elber | Geometric constraint solver using multivariate rational spline functions[END_REF], the emphasis is put on the subdivision process, and stopping criterion based on the normal cone to the surface patch. In [START_REF] Mourrain | Subdivision methods for solving polynomial equations[END_REF], this approach has been improved by introducing pre-conditioning and univariatesolver steps. The complexity of the method is also analyzed in terms of intrinsic differential invariants.

The application of subdivision methods for handling higher dimensional objects is not so well developed. In [START_REF] Joon-Kyung | Contouring 1and 2-Manifolds in Arbitrary Dimensions[END_REF] a method which subdivides up to some precision level, and applies dual marching cube approach to connect points on the curve or to mesh a surface is described. The variety is covered by boxes of a given size, and the connectivity of these cells is used to deduce the piecewise linear approximation. In [START_REF] Alberti | Meshing implicit algebraic surfaces: the smooth case[END_REF], a subdivision approach exploiting the sign variation of the coefficients in the Bernstein basis in order to certify the topology of the surface in a cell, is used for the purpose of polygonalizing an implicit algebraic surface.

The work of this paper is in the spirit of this former approach. We apply a subdivision approach also exploiting the properties of the Bernstein polynomial representation. We describe a simple regularity test extending the criterion of [START_REF] Alberti | Meshing implicit algebraic surfaces: the smooth case[END_REF] to curves, which allows us to detect easily when the topology of the curve in a cell is uniquely determined from its intersection with the border of the cell. This provides an efficient test for stopping earlier the subdivision process and branching to path following methods if we are interested in a good geometrical approximation of the curve.

We address the same question as in [START_REF] Gatellier | Computing the topology of 3-dimensional algebraic curves[END_REF], but with this new methods, we are able to solve the following problems already identified in this paper: -To achieve higher numerical stability by operating on Bernstein basis instead of monomial basis; -Through subdivision on three principle directions, i.e. x, y, z (or x, y, z), to isolate the domain containing the singularities from those containing regular curve segments. This divide-and-conquer approach, in principle, should simplify the graph building algorithm adopted in [START_REF] Gatellier | Computing the topology of 3-dimensional algebraic curves[END_REF] where the whole domain has to be considered.

However, for the treatment of singular points, we have to introduced a threshold to stop the subdivision. Contrarily to [START_REF] Gatellier | Computing the topology of 3-dimensional algebraic curves[END_REF], we do not certify the topology at singular points, but computed boxes of size , containing these singularities.

On the contrary, we show that our approach is able to handle implicit curves with large equation (of total degree about 80 with coefficients of bit-size 200), which resultant-based techniques are not able to treat.

This paper is organized as follows: in Section 2, we will review some of the relevant concepts and theorem required by our proposed algorithm; Section 3 is devoted to outlining our proposed algorithms and the details about how the essential steps in our algorithm are handled. We will show the experiment results in Section 4 and conclude in Section 5 with the problems and possible improvements over the currently proposed algorithm.

Fundamental ingredients

This section introduces the theoretical background of Bernstein polynomial representation and how it is related to the problem we want to solve. For a domain D ⊂ n , we denote by

• D its interior, by D its closure. For a box D = [a 0 , b 0 ] × [a 1 , b 1 ] × [a 2 , b 2 ] ⊂ 3 ,
its x-face (resp. y-face, z-face) are its faces normal to the direction x (resp. y, z).

Univariate Bernstein Basis

Given an arbitrary univariate polynomial function f (x) ∈ ¡ , we can convert it into the representation of Bernstein basis of degree d, which is defined by:

f (x) = i b i B d i (x), and (1)
B d i (x) = d i x i (1 -x) d-i (2)
where b i is usually referred as controlling coefficients. Such conversion is done through a basis conversion [START_REF] Farin | Curves and Surfaces for Computer Aided Geometric Design: A Practical Guide[END_REF]. The above formula can be generalized to an arbitrary interval [a, b] by a variable substitution x = (ba)x + a. We denote by

B i d (x; a, b) d i (x -a) i (b -x) d-i (b -a) -d the corresponding Bernstein basis on [a, b].
There are several useful properties regarding Bernstein basis given as follows:

-Convex-Hull Properties:

Since i B i d (x; a, b) ≡ 1 and ∀x ∈ [a, b], B i d (x;
a, b) ≥ 0 where i = 0, ..., d, the graph of f (x) = 0, which is given by (x, f (x)), should always lie within the convex-hull defined by the control coefficients [START_REF] Farin | Curves and surfaces for computer aided geometric design : a practical guide[END_REF].

-Subdivision (de Casteljau): Given t 0 ∈ [0, 1], f (x) can be represented piecewisely by:

f (x) = d i=0 b (i) 0 B i d (x; a, c) = d i=0 b (d-i) i B i d (x; c, b), where (3) b
(k) i = (1 -t 0 )b (k-1) i + t 0 b (k-1) i+1 and c = (1 -t 0 )a + t 0 b. (4)
Another interesting property of this representation is related to Descartes's Law of signs. The definition of Descartes's Law for a sequence of coefficients b k = b i |i = 1, ..., k is defined recursively:

V (b k+1 ) = V (b k ) + 1, if b i b i+1 < 0 0, else (5)
With this definition, we have:

Theorem 1. Given a polynomial f (x) = n i b i B d i (x; a, b), the number N of real roots of f on ]a, b[ is less than or equal to V (b), where b = (b i ), i = 1, ..., n and N ≡ V (b) mod 2.
The theorem 1 enables a simple yet efficient test of the existence of real roots in a given domain. This test is essential to our algorithm, as it serves as a key criterion to classify whether a domain has certified topology, without actually computing the curve. This allow the our algorithm to execute in reasonably short time, as demonstrated in our experiments.

Generalization to Multivariate Case

The univariate Bernstein basis representation can be generalized to multivariate ones. Briefly speaking, we can rewrite the definition (Eq. ( 1)) in the form of tensor products. Suppose for x = (x 0 , ..., x n-1 ) ∈ n , f = (x) ∈ ¡ [x] having the maximum degree d = (d 0 , ..., d n-1 ) has the form:

f (x) = d0 k0=0 ... dn kn=0 b k1,...,kn B d0 k0 (x 0 )...B d0 kn (x n ) (6)
For a polynomial of n variables, the coefficients can be viewed as a tensor of dimension n. The de Casteljau subdivision for the multivariate case proceeds similarly to the univariate one, since the subdivision can be done independently with regards to a particular variable x i .

Based on these properties, a subdivision solver which can be seen as an improvement of the Interval Projected Polyhedron algorithm in [START_REF] Sherbrooke | Computation of the solutions of nonlinear polynomial systems[END_REF], is described in [START_REF] Mourrain | Subdivision methods for solving polynomial equations[END_REF]. It uses the following operations: The multivariate functions to be solved are enclosed in-between two univariate functions, for each variable. For this purpose, the Bernstein control points of the functions are projected in each direction and the upper and lower envelop are used to define these enveloping univariate polynomials. A lower and upper approximation of the roots of these univariate polynomials are used to reduce the domain. If the reduction is not sufficient, the domain is split. These reduction operations are improved by pre-conditioning steps. See [START_REF] Mourrain | Subdivision methods for solving polynomial equations[END_REF] for more details.

Algorithmic ingredients

We consider the problem of computing the topology of the curve, denoted hereafter as C, resulting from the intersection of two known algebraic surfaces, namely, f (x) = 0 and g(x) = 0 defined in 3 , with f, g ∈ [x, y, z]. Our discussion is confined to the case where f and g has no common divisor other than 1, so that their intersection has dimension of 1. We assume moreover that (f, g) is radical or equivalently that the resultant of f (x, y, z), g(x, y, z) with respect to z after a generic change of coordinates, is square free.

Tangent Vector Field

The tangent vector on C serves as the key to our analysis of topology of the curve. It serves as an important indicator of topological feature of C. While it is computationally prohibitive to compute the tangent vector at each point on C, we can reach some useful conclusion about the topology of the curve by looking into the tangent vector field defined below:

t = t x (x)e x + t y (x)e y + t z (x) e z = f ∧ g = e x e y e z ∂ x f ∂ y f ∂ z f ∂ x g ∂ y g ∂ z g (7)
where e x , e y and e z are the unit vectors along the principle axis x, y and z, respectively; t x , t y and t z are functions of x = (x, y, z).

Singularities on the curve can be easily characterized, as t vanishes at those points. In [START_REF] Gatellier | Computing the topology of 3-dimensional algebraic curves[END_REF], the author also tried to localize the point having a tangent parallel to a virtual sweeping plane. They are connected together with the singularities to form the final topological graph. In order to do this, the whole curve is projected onto some principle projection planes. However, the projected planar curve in many cases has a very different topology as C. In our proposed algorithm, we exploit the subdivision along all three principle axes simultaneously and the critical events are either reduced to regular case (such as for tangents) or localized (such as for intersections). The topology graph can be built without explicitly computing the exact position of the singularities.

Regularity Test

In this section, we are going now to describe how to detect boxes, for which the topology of the curve can be determined. We will use the following notions: Definition 1. We say that a curve C ∈ n is regular in a compact domain D ⊂ n , if its topology is uniquely determined from its intersection with the boundary D.

The aim of the method is to give a simple criterion for the regularity of a curve in a box.

To form the topological graph for this domain, we only need to compute the intersections between the curve and the boundary of this domain, and there exists a unique graph to link these intersections so that this graph complies to the true topology of the original curve.

2D case: For 2D planar algebraic curve C defined by a polynomial equation f (x, y) = 0, and denoting the partial derivative of f w.r.t x by ∂ x f , we have the following direct property:

Proposition 1. If ∂ y f (x, y) = 0 (resp. ∂ x f (x, y) = 0) in a domain D = [a 0 , b 0 ] × [a 1 , b 1 ] ⊂ 2 , the curve C is regular on D.
Proof. Suppose that ∂ y f (x, y) = 0 in D. Then C is smooth, since its normal vector is defined everywhere, and has no vertical tangents in D. By the implicit function theorem, the connected components of C∩ Moreover, as there is no vertical tangent, each of these branches starts and ends at a point on the border ∂D of D. Notice that two branches may share a starting or ending point, when the curve is tangent (with even multiplicity) to ∂D.

Thus, computing the points of C ∩ ∂D, repeating a point if its multiplicity is even, sorting them by lexicographic order such that x > y ((x 0 , y 0 ) > (x 1 , y 1 ) if x 0 > x 1 or x 0 = x 1 and y 0 > y 1 ), we obtain a sequence of points p 1 , p 2 , . . . , p 2 s-1 , p 2 s such that the curve C in D is isotopic to the union of the non-intersecting segments [p 1 , p 2 ], . . . , [p 2s-1 , p 2s ]. In other words, the topology of C is uniquely determined from its intersection points with ∂D and C is regular on D. If ∂ y f = 0 on D (resp. ∂ x f = 0), we will say that C is x-regular (resp. y-regular).

A sufficient condition for f to be x-regular (respectively y-regular) is that the Bernstein coefficients of the first derivative of f against y (respectively x) maintains a constant sign (see also [START_REF] Alberti | Meshing implicit algebraic surfaces: the smooth case[END_REF]). By Descartes' law, this statement implies that the sign variation in this direction should be at most 1.

To put it in another way, by solely studying the sign variations of the tangential gradient vector of the curve (represented in Bernstein basis), i.e. (∂ y f (x), -∂ x f (x)), we are able to detect when the curve is regular on D and to determine uniquely the topological graph.

3D case: The 2D approach can be generalized to the 3D case where the tangential gradient vector of the curve C defined by the intersection of two algebraic surfaces, namely f (x, y, z) = 0 and g(x, y, z) = 0, is given by t = (f ) ∧ (g) (see Eh. ( 7)). Similar to the 2D case, we can represent each component of t in the Bernstein basis for a given domain (in cube shape

) D = [a 0 , b 0 ]×[a 1 , b 1 ]×[a 2 , b 2 ].
The sign change of the resulting Bernstein coefficients enables a simple regularity test with minimal computation effort.

We describe a first and simple regularity criterion:

Proposition 2. The 3D spatial curve C defined by f = 0 and g = 0 is regular on D, if t x (x) = 0 on D, and -∂ y h = 0 (or ∂ z h = 0) on D, for h = f or h = g.

Proof. Suppose that t x (x) = 0 and ∂ z (f ) = 0 on D. It implies that C is smooth in D. Consider two branches of C in D and project them by π z onto a (x, y)-plane.

Their projection cannot intersect at an interior point. Otherwise, there would be two points p 1 , p 2 ∈ D, such that f (p 1 ) = 0, f (p 2 ) = 0 and π z (p 1 ) = π z (p 2 ), which implies that ∂ z f (p) vanishes for an intermediate point ∈]p 1 , p 2 [ in D. This is impossible by hypothesis. Consequently, the branches of C project bijectively onto the branches of π z (C). Their tangent vector is the projection (t x (x), t y (x)) of the tangent vector of C. By proposition 1, π z (C) is regular, so that the topology of π z (C) , and thus of C, is uniquely determined by the intersection points of C with the border of D.

A similar criterion applies by symmetry, exchanging the roles of the x, y, z coordinates.

Let us give now a finer regularity criterion, which is computationally less expensive: Proposition 3. If C is smooth in D and if for all x 0 ∈ , the plane x = x 0 plane has at most one intersection point with the curve C in D, then C is regular on D.

Proof. Consider the projection π z (C) of the curve C in D along the z direction. Then the components of C in D projects bijectively on the (y, z) plane. Otherwise, there exist two points p 0 and p 1 lying on C such that π z (p 0 ) = π z (p 1 ) = (x 0 , y 0 ), then p 0 and p 1 belong to x = x 0 which are functions of the form y = Φ(x). Otherwise, there exist two points on π z (C) and (and on C ∩ D) with the same x-coordinate. Consequently, for x ∈ [a 0 , b 0 ] there is at most one branch of π z (C) in D above x, and the connected components of C∩

Proof. Let us fix x 0 ∈ [a 0 , b 0 ] where D = [a 0 , b 0 ] × [a 1 , b 1 ] × [a 2 , b 2 ], let U = {x 0 } × [a 1 , b 1 ] × [a 2 , b 2 ] and let Φ x0 : (x 0 , y, z) ∈ U → (f (x 0 , y, z), g(x 0 , y, z)).
We are going to prove that under our hypothesis, Φ x0 is injective. The Jacobian t x (x 0 , y, z) of Φ x0 does not vanish on U , so that Φ x0 is locally injective. We consider the level-set f (x) = f 0 for some f 0 ∈ f (U ). It cannot contain a closed loop in U , otherwise we would have (∂ y f, ∂ z f ) = 0 (and thus t x = 0) in U ⊂ D.

We deduce that each connected component of f (x) = f 0 in U intersects ∂U in two points. Now suppose that Φ x0 is not injective on U , so that we have two points

p 1 , p 2 ∈ U such that Φ x0 (p 1 ) = Φ x0 (p 2 ).
If p 1 and p 2 are on the same connected component of the level set f (x) = f 0 (where f 0 = f (p 1 ) = f (p 2 )) in U , then g reaches the same value at p 1 and p 2 on this level set, so that by Role's theorem, there exists a point p ∈ U in-between p 1 and p 2 , such that Jac(Φ x0 )(p) = t x (p) = 0. By hypothesis, this is impossible.

Thus p 1 and p 2 belongs to two different connected components of f (x) = f 0 in U . Consequently the value f 0 is reached at 4 distinct points of ∂U , which implies that f has at least 4 extrema on ∂U . Now note that up to a change of variable z = a 2z, we can assume that

∂ z f > 0 on both y = a 1 , y = b 1 faces. Then if ∂ y f < 0 on z = a 2 , we have f (x 0 , a 1 , b 2 ) > f (x 0 , a 1 , a 2 ) > f (x 0 , b 1 , a 2 ) and (a 2 , a 3 ) is not a local extrema.
Otherwise ∂ y f > 0 and (b 1 , a 2 ) is not a local extrema. In both cases, we do not have 4 extrema, which proves that φ x0 is injective and that the intersection of C with the plane x = x 0 in D is at most one point. So by proposition 3, we deduce that C is regular in D.

For more details on the injectivity properties, see [START_REF] Pavone | Auto-intersection de surfaces pamatrées réelles[END_REF]. Here also, a similar criterion applies by symmetry, exchanging the roles of the x, y, z coordinates.

If one of these criteria applies with t i (x) = 0 on D (for i = x, y, z), we will say that C is i-regular on D.

From a practical point of view, the test that t i (x) = 0 or ∂ i (h) for i = x, y or z, h = f or g, is replaced by the stronger condition that their coefficients on the Bernstein basis of D have a constant sign, which is straightforward to check. Similarly, such a property on the faces of D is also direct, since the coefficients of a polynomial, with a minimal (resp. maximal) x-indices (resp. y-indices, zindices) are its Bernstein coefficients on the corresponding face.

In addition to these tests, we also test whether both surfaces penetrate the cell, since a point on the curve must lie on both surfaces. This test could be done by looking at the sign change of the Bernstein coefficients of the surfaces with regards to that cell. If no sign change occurs, we can rule out the possibility that the cell contains any portion of the curve C, hence terminate the subdivision early. In this case, we will also say that the cell is regular.

The regularity criterion is sufficient for us to uniquely construct the topological graph g of C within D. Without loss of generality, we suppose that the curve C is x-regular in D. Hence, there is no singularity of C in D. Furthermore, this also guarantees that there is no 'turning-back' of the curve tangent along x-direction, so the mapping of C onto the x axis is injective. Intuitively, the mapped curve should be a series of non-overlapping line segments, of which the ends correspond to the intersections between the curve C and the cell, and such mapping is injective.

This property leads to a unique way to connect those intersection points, once they are computed (see section 3.3), in order to obtain a graph representing the topology of C. Here is how this graph is computed in practice: suppose C is i-regular in the domain D, and that we have computed the set of intersection points V = {v j } of the curve with the boundary of D. First, we sort the elements in V comparing vectors by their i-th coordinates. Assuming the sorted points v j are indexed by j = 0, 1, 2, ..., we form the edges v k , v k+1 , for k = 0, 2, 4, ... However, a special case has to be taken into account, that is when v j has a multiplicity m i > 1, for instance, when C is tangent to the bounding domain D at v i . In this case, we can treat v j conceptually as a multiple point which plays the role of m i points. In this way, we proceed the connecting process in the same manner as we do for the general case. To determine the multiplicity of a point v j , we only have to evaluate the derivatives of C at this point.

Hierarchical subdivision

We adopted a hierarchical octree to partition the 3 space, for several reasons:

each cell of the octree is equivalent to a cube-shaped domain D; which stores the coefficients of the polynomials in the Bernstein basis of the corresponding domain.

we can take cares of faces shared by cells, to minimize the number of calls to solvers; -the hierarchical structure of octree allows us to terminate (stop further subdivision) early when a cell is deemed regular or irrelevant.

We begin by setting a initial bounding domain D 0 to a root cell. A cell is subdivided if the curve C defined in the correspondent domain fails the regularity test. For each subdivision, we result in several smaller domains in form of subcells. For each of them, we repeat the regularity test and, if necessary, further subdivides. The subdivision of a cell will terminate either when the curve within is deemed regular, or the size of the cell is beyond a predefined precision .

There are several techniques to save computation efforts. As the sub-cells share certain faces with their parent cell, the earlier computed intersections on the parent cell's faces are inherited directly by the sub-cells. In addition, subcells split from the same parent cell do share some faces as well. Once again, the shared faces should be computed exactly once.

Once a new face is introduced in the octree decomposition, the bivariate solver described in section 2.2, is called directly with the Bernstein coefficients of the polynomials on this face. The points we found are shared by neighbor cells, connected to this face in the octree.

Symbolic-numeric approach

Some geometric operations such as computing the self-intersection curve or the ridge curve of a parameterized surface leads to the computation of implicit curves of high degree with coefficients of large size. This is either due to projection techniques (see [START_REF] Galligo | Self-intersections of a Bézier bicubic surface[END_REF]), or to their definition through composed operations (see [START_REF] Cazals | Topologically certified approximation of umbilics and ridges on polynomial parametric surface[END_REF]). In order to be able to handle such curve, the main difficulty is to control the result, using approximate computation, since exact computation though possible, would be prohibitive. We describe here the symbolic-numeric approach that we have developed for this purpose.

We assume that the input equations are given with exact (large) rational (or integer) numbers (even if the input is given with floating point numbers, we will consider it as an exact input). In order to compute the topology of C in a domain D, we convert its representation in the Bernstein basis of D, using exact rational arithmetic.

Once this conversion is done, we normalize the equation, by dividing by the coefficient of maximal norm. For each resulting rational coefficient c, we compute the smallest interval [c, c] represented with floating point numbers and containing c.

Then, the subdivision process is performed, using interval arithmetic. The regularity criterion, which reduces to sign evaluations, is applied on these interval coefficients. We use the following convention: a interval is < 0 (resp. > 0) if all its elements are < 0 (resp. > 0). If the interval contains 0, we say that its sign is indeterminate.

If the regularity test fails, -either the sign of all the coefficients of the polynomial are indeterminate, and we re-convert the exact polynomial to its representation on the corresponding sub-domain and restart the approximation process. -or we subdivide the domain, as in the usual case.

Algorithm Outline

The proposed algorithm for 3D curves is outlined as following:

Algorithm 31 Computing the topology of the curve C: Input: f(x) and g(x) polynomials ∈ ¢ [x, y, z], a tolerance and a list of bounding domain Output: The graph represented by a set of vertices V , which are either 3D points or boxes (with size less than ) bounding the singularities, and a set of connections E that are representing the edges of the resulting graph.

D 0 ← [a 0 , b 0 ] × [a 1 , b 1 ] × [a 2 , b 2 ] (a i,bi ∈ ).
We do not describe the algorithm for 2D curves, which is basically a specialization of this one.

Experiments

Our proposed algorithm is implemented as a part of synaps (SYmbolic Numeric APplicationS) library2 . The experiments have been carried out on a 3.4GHz PC, under Linux.

Planar curves of high degree with large coefficients

In this section, we report on the application of the 2d algorithm, in the case of large integer coefficients. The first example is about ridge curve. Ridge curves correspond to local extrema of curvature taken in the principal direction of the surface, after some algebraic manipulations they can be obtained as implicit curve (see [START_REF] Cazals | Topologically certified approximation of umbilics and ridges on polynomial parametric surface[END_REF]). See also [START_REF] Thirion | The 3D Marching Lines Algorithm and its Applications to Crest Lines Extraction[END_REF] and [START_REF] Thirion | Computing the Differential Characteristics of Isointensity Surfaces[END_REF] for other related approaches.

In the example it corresponds to a bicubic surface, the input polynomial is of total degree 84, of multidegree (43, 43) with 1907 monomials. The coefficients are integers encoded on at most 65 bits. For the precision = 10 -3 which controls the singularity localization, it takes 30 seconds. The topology is certified except in tiny boxes (which contains the singularity points). Notice that a pure algebraic approach, exploiting the specificity of problem and with a very efficient Gröbner engine takes about 10 minutes to certify the topology (see [START_REF] Cazals | Topologically certified approximation of umbilics and ridges on polynomial parametric surface[END_REF]).

The second example is a projection of a self-intersection curve of bicubic patch, computed by resultant techniques (see [START_REF] Galligo | Self-intersections of a Bézier bicubic surface[END_REF]). The input polynomial is of total degree 76, of multidegree (44, 44) with 1905 monomials. The coefficients are integers of at most 288 bits. It takes 5 seconds, for this example with the same precision = 10 -3 . 1) f (x) = 0.85934x 2 + 0.259387xy + 0.880419y 2 + 0.524937xz -0.484008yz + 0.510242z 2 -1 g(x) = 0.95309x 2 + 0.303149xy + 0.510242y 2 -0.200075xz + 0.64647yz + 0.786669z 2 -1 time: 80 msec 2) f (x) = -0.125x 2 -0.0583493xy + 0.493569y 2 + 0.966682xz -1.5073yz -0.368569z 2 -0.865971x -0.433067y -0.250095z g

(x) = x 2 + y 2 + z 2 -2 time: 20 msec 3) f (x) = 2x 2 + y 2 + z 2 -4 g(x) = x 2 + 2xy + y 2 -2yz -2z 2 + 2zx time: 30 msec 4) f (x) = x 4 + y 4 + 2x 2 y 2 + 2x 2 + 2y 2 -x -y -z g(x) = x 4 + 2x 2 y 2 + y 4 + 3x 2 y -y 3 + z 2 time: 130 msec (

Silhouette curves of implicit surfaces

The following samples are taken from http://www-sop.inria.fr/galaad/surfaces/.

We intersect the surface with its polar variety in one direction (here the x direction). In other words, we intersect the surface with the surface defined by one of its first order derivative (here ∂ x f ), to extract its silhouette. The surfaces that we used are called respectively Tetrahedral, Q3, Q1 and Barth Sextic (see Fig. 3): 5) f (x) = x 4 + 2x 2 y 2 + 2x 2 z 2 + y 4 + 2y 2 z 2 + z 4 + 8xyz -10x 2 -10y 2 -10z 2 + 25 g(x) = 4x 3 + 4xy 2 + 4xz 2 + 8yz -20x time: 510 msec 6) f (x) = 5.229914547374508y 2 z 2 + 3.597883597883598x 2 y 2 + y 4 + z 4x 4 -19.49816368932737xyz + 5.229914547374508x 2 -7.43880040039534y 2 g(x) = -3.597883597883598z 2 +7.43880040039534z 2 x 2 -110.45982909yz 2 + 7.195767196x 2 y + 4y 3 -19.49816368932737xz -14.87760080y time: 330 msec 7) f (x) = x 4 +y 4 +z 4 -4x 2 -4y 2 z 2 -4y 2 -4z 2 x 2 -4z 2 -4x 2 y 2 +20.7846xyz +1 g(x) = 4x 3 -8x -8xz 2 -8xy 2 + 20.7846yz time: 730 msec 8) f (x) = 67.77708776x 2 y 2 z 2 -27.41640789x 4 y 2 -27.41640789x 2 z 4 +10.47213596x 4 z 2 -27.41640789y 4 z 2 + 10.47213596y 4 x 2 + 10.47213596y

Discussion

The algorithm proposed in this paper offers a generic method for computing the topological graph of spatial curves resulting from the intersection of two algebraic surfaces. As demonstrated in the experiments, it is rather robust despite the increase of the complexity of the curve. The major weakness of this approach, however, is that certain apparently simply situation could result in a lot of subdivisions, such as the curve with parallel structures which are very close to each other.

As specified, in this method, the singular points are only isolated into small boxes of size , and we do not certify the connection of the branches at these points. An additional work would be necessary, to certify the singularity type. We are currently investigating on this problem.

Another weakness of this approach is the memory consumption (as the storage requirement is approximately cubic to the depth of the subdivision). This problem however, can be mitigated by a divide-and-conquer approach and pure programming techniques.

•D

  are the graph of functions of the form y = φ(x). The closure of such a connected component is called hereafter a branch of C in D. As ∂ y f (x, y) = 0 in D, for a given x ∈ [a 0 , b 0 ] there is at most one branch of C in D above x. Consequently, the connected components of C∩ • D project bijectively onto non-overlapping open intervals of [a 0 , b 0 ].

•

  D project bijectively onto non-overlapping open intervals of [a 0 , b 0 ] as π z (C) does. We conclude as in the 2D case (proposition 1), by sorting the points of C ∩ ∂D according to their x-coordinates, and by gathering them by consecutive pairs corresponding to the starting and ending points of branches of C ∩ D. Proposition 4. The 3D spatial curve C defined by f = 0 and g = 0 is regular on D, if t x (x) = 0 on D, and -∂ y h = 0 on z-faces, and ∂ z h = 0 and its has the same sign on both y-faces of D, for h = f or h = g.

-

  Step 0: (initialization step) domain list D ← D 0 ; vertex list V ← NIL; connectivity list E ← NIL; -Step 1: compute t ← (f ) ∧ (g) given by Eq. (7); -Step 2: convert f , g and t into Bernstein basis representation; -Step 3: while D is not empty, pick a D in D:• Step 3.1: compute V the set of intersection points between the boundary of the domain D and the curve C; • Step 3.2: if the size of D is larger than :* if the curve C within the domain D is regular (see section 3.2):• sort and connect the points v ∈ V ; the connectivities are stored in E; * else if the domain D is not regular:• subdivide D and append the subdivided domains into the domain list D • else if the size of D is not larger than :* add the domain D as a 'box' vertex into V ; * this vertex is connected with all intersections v ∈ V of D; these connectivities are also appended to E; • Step 3.4: remove D from D and repeat Step 3;

Fig. 1 .

 1 Fig. 1. Topological descriptions of high complexity curves

Fig. 2 .

 2 Fig. 2. Topological descriptions of the intersection curve for 4 pairs of low-order algebraic surfaces.

Fig. 3 .

 3 Fig. 3. Topological descriptions of the silhouette curves of algebraic surfaces.

  2 z 4 -4.236067978x 4 -8.472135956x 2 y 2 -8.472135956x 2 z 2 +8.472135956x 2 -4.236067978y 4 -8.472135956y 2 z 2 + 8.472135956y 2 -4.236067978z 4 + 8.472135956z 2 -4.236067978 g(x) = 135.5541755xy 2 z 2 -109.6656316x 3 y 2 -54.83281578xz 4 +41.88854384x 3 z 2 + 20.94427192y 4 x-16.94427191x 3 -16.94427191xy 2 -16.94427191xz 2 +16.94427191x time: 4010 msec

http://www-sop.inria.fr/galaad/software/synaps/

http://www-sop.inria.fr/galaad/software/axel/

Acknowledgements: We would like to thanks J. Wintz, for his very nice software axel 3 for the visualisation and manipulation of algebraic objects, that we used to produce the pictures of the curves. We acknowledge the partial support of Aim@Shape (IST NoE 506766) and ACS (IST Fet Open 006413).