Ridges and Umbilics of Polynomial Parametric Surfaces Frédéric Cazals, Jean-Charles Faugère, Marc Pouget, Fabrice Rouillier To cite this version: Frédéric Cazals, Jean-Charles Faugère, Marc Pouget, Fabrice Rouillier. Ridges and Umbilics of Polynomial Parametric Surfaces. B. Juttler and R. Piene. Geometric Modeling and Algebraic Geometry, Springer, pp.141–159, 2008, 978-3-540-72184-0. <10.1007/978-3-540-72185-7 8>. <inria-00329762> HAL Id: inria-00329762 https://hal.inria.fr/inria-00329762 Submitted on 13 Oct 2008 HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers. L’archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d’enseignement et de recherche français ou étrangers, des laboratoires publics ou privés. Ridges and umbilics of polynomial parametric surfaces Frédéric Cazals, Jean-Charles Faugère, Marc Pouget, and Fabrice Rouillier INRIA Sophia-Antipolis, Geometrica project, 2004 route des Lucioles, BP 93, F-06902 Sophia-Antipolis, FRANCE. Frederic.Cazals@sophia.inria.fr Marc.Pouget@gmail.com and INRIA Rocquencourt and Université Pierre et Marie Curie-Paris6, UMR 7606, LIP6, Salsa project, Domaine de Voluceau, BP 105, F-78153 Le Chesnay Cedex, FRANCE. Jean-Charles.Faugere@inria.fr Fabrice.Rouillier@inria.fr Abstract. Given a smooth surface, a blue (red) ridge is a curve along which the maximum (minimum) principal curvature has an extremum along its curvature line. Ridges are curves of extremal curvature and therefore encode important informations used in segmentation, registration, matching and surface analysis. State of the art methods for ridge extraction either report red and blue ridges simultaneously or separately —in which case a local orientation procedure of principal directions is needed, but no method developed so far topologically certifies the curves reported. In this context, we make two contributions. First, for any smooth parametric surface, we exhibit the implicit equation P = 0 of the singular curve P encoding all ridges of the surface (blue and red), we analyze its singularities and we explain how colors can be recovered. Second, we instantiate to the algebraic setting the implicit equation P = 0. For a polynomial surface, this equation defines an algebraic curve, and we develop the first certified algorithm to produce a topologically certified approximation of it. The algorithm exploits the singular structure of P —umbilics and purple points, and reduces the problem to solving zero dimensional systems using Rational Univariate Representations and isolate roots of univariate rational polynomials. An experimental section illustrates the efficiency of the algorithm on a Bezier patch. 1 Introduction 1.1 Ridges Differential properties of surfaces embedded in R3 are a fascinating topic per se, and have long been of interest for artists and mathematicians, as illustrated by the parabolic lines drawn by Felix Klein on the Apollo of Belvedere [HCV52], and also by the developments reported in [Koe90]. Beyond these noble considerations, the recent development of laser range scanners and medical images shed light on the importance of being 2 F. Cazals, J.C. Faugère, M. Pouget, F. Rouillier able to analyze discrete datasets consisting of point clouds in 3D or medical images — grids of 3D voxels. Whenever the datasets processed model piecewise smooth surfaces, a precise description of the models naturally calls for differential properties. In particular, applications such as shape matching [HGY+ 99], surface analysis [HGY+ 99], or registration [PAT00] require the characterization of high order properties and in particular the characterization of curves of extremal curvatures, which are precisely the so-called ridges. Interestingly, (selected) ridges are also central in the analysis of Delaunay based surface meshing algorithms [ABL03]. A comprehensive description of ridges can be found in [Por71,HGY+ 99,Por01], and in the sequel, we just introduce the basic notions so as to discuss our contributions. Consider a smooth embedded surface, and denote k1 and k2 the principal curvatures, with k1 ≥ k2 . Umbilics are the points where k1 = k2 . For any non umbilical point, the corresponding principal directions of curvature are well defined, and we denote them d1 and d2 . In local coordinates, we denote h, i the inner product induced by the ambient Euclidean space, and dk1 , dk2 the gradients of the principal curvatures. We then define ridges (see also fig. 1 and 2 for an illustration) : Definition 1. A non umbilical point is called – a blue ridge point if the extremality coefficient b0 = hdk1 , d1 i vanishes, i.e. b0 = 0. – a red ridge point if the extremality coefficient b3 = hdk2 , d2 i vanishes, i.e. b3 = 0. Notice that, as the principal curvatures are not differentiable at umbilics, the extremality coefficients are not defined at such points. In addition, the sign of the extremality coefficients is not defined since each principal direction can be oriented by two opposite unit vectors. Apart from umbilics, special points on ridges are purple points, which correspond to intersections between red and a blue ridges. The calculation of ridges poses difficult problems, which are of three kinds. Topological difficulties. Ridges of a smooth surface form a singular curve on the surface, with self-intersections at umbilics (more precisely at so-called 3-ridges umbilics), and purple points. From a topological viewpoint, reporting ridges of a surface faithfully requires reporting umbilics and purple points. Numerical difficulties. As pointed out above, ridges are characterized through third order derivatives of the surface. Estimating them depends on the particular type of surface processed —implicitly defined, parameterized, discretized by a mesh— and is numerically a difficult task. Orientation difficulties. Since coefficients b0 and b3 depend on a given orientation of the principal directions, their sign is not well defined. Practically, tracking the sign change of functions whose sign depends on the particular orientation of the frame in which they are expressed poses a problem. In particular, tracking a zero-crossing of b0 or b3 from sign changes along a curve segment on the surface imposes to find a coherent orientation of the principal frame at the endpoints. Given two principal directions at these endpoints, one way to find a local orientation consists of choosing two vectors so that they make an acute angle, whence the name Acute Rule. This rule has been Ridges and umbilics of polynomial parametric surfaces 3 used since the very beginning of computer examination of ridges [Mor90,Mor96], and is implicitly used in almost all algorithms. But the question of specifying conditions guaranteeing the decisions made are correct has only been addressed recently [CP05]. Another approach is to extract the zero level set of the Gaussian extremality Eg = b0 b3 defined in [Thi96]. This function has a well defined sign independent from the orientation, but it is still not defined at umbilics and one cannot distinguish blue from red ridges. Fig. 1: Umbilics, ridges, and principal blue foliation on the ellipsoid Fig. 2: Schematic view of the ridges with umbilics (black points) and purple points (purple). Max of k1 : blue; Min of k1 : green; Min of k2 : red; Max of k2 : yellow 1.2 Previous work Given the previous difficulties, no algorithm reporting ridges in a certified fashion has been developed as of today. Most contributions deal with sampled surfaces known 4 F. Cazals, J.C. Faugère, M. Pouget, F. Rouillier through a mesh, and a complete review of these contributions can be found in [CP05]. In the following, we focus on contributions related to parametric surfaces. Reporting umbilics. Umbilics of a surface are always traversed by ridges, so that reporting ridges faithfully requires reporting umbilics. To do so, Morris [Mor90] minimizes the function k1 −k2 , which vanishes exactly at umbilics. Meakawa et al. [MWP96] define a polynomial system whose roots are the umbilics. This system is solved with the rounded interval arithmetic projected polyhedron method. This algorithm uses specific properties of the Bernstein basis of polynomials and interval arithmetic. The domain is recursively subdivided and a set of boxes containing the umbilics is output, but neither existence nor uniqueness of an umbilic in a box is guaranteed. Reporting ridges. The only method dedicated to parametric surfaces we are aware of is that of Morris [Mor90,Mor96]. The parametric domain is triangulated and zero crossings are sought on edges. Local orientation of the principal directions are needed but only provided with a heuristic. This enables to detect crossings assuming (i)there is at most one such crossing on an edge (ii)the orientation of the principal directions is correct. As this simple algorithm fails near umbilics, these points are located first and crossings are found on a circle around the umbilic. Equation of the ridge curve. Ridges can be characterized either as extrema of principal curvatures along their curvature lines as in definition 1, or by analyzing the contact between the surface and spheres [HGY+ 99]. For parametric surfaces, this later approach allows a global characterization of ridges [Por01, Chapter 11] as a 1 dimensional smooth submanifold in a 7 dimensional space. For any point c ∈ R3 , the contact function between the surface and the sphere centered at c is defined by: V (c) = 1 h c, c i − h s − c, s − c i = h c, s i − h s, s i. 2 2 (1) The surface s is parameterized by w ∈ R2 . To see the connection with definition 1, let V (c)k denote the k-times linear form of the derivative of order k w.r.t. w associated to point c. Ridges are essentially characterized in the following theorem [Por01, Thm 11.10 page 191]: Theorem 2. Let s be a regular surface in R3 with only ordinary umbilics. Then these points (c, w, u) of R3 × R2 × R2 where V (c)1 (w) = 0, V (c)2 (w)u = 0, V (c)3 (w)u3 = 0, || u ||= 1 (2) form a smooth submanifold of R3 × R2 × R2 of dimension 1, with certain exceptions that will appear in the course of the proof. Shifting from this seven-dimensional space to the parametric space, the theory of algebraic invariants is used to derive the equation of the ridge curve as the zero set of an invariant function in [Gra05]. Reporting the topology of an algebraic curve. In the case of a polynomial parametric surface, we recast the problem of approximating ridges into the field of algebraic geometry. We recall that the standard tool to compute a graph encoding the topology of a 2-D or 3-D curve is the Cylindrical Algebraic Decomposition (CAD) [GVN02,GLMT04]. Ridges and umbilics of polynomial parametric surfaces 5 1.3 Contributions and paper overview Given the difficulties recalled in section 1.1, the ultimate goal is the development of a certified algorithm reporting the topology of ridges. Let Φ(u, v) be a smooth parameterized surface over a domain D ⊂ R2 . We meet this goal in two steps. Implicitization of the ridge curve. First, we exhibit the implicit equation P = 0 of the singular curve P encoding all ridges of the surface (blue and red), and show how to recover the colors from factors of P . We also derive zero dimensional systems coding the different singularities of P . These contributions are fundamentally different from previous work. With respect to [Por01], we have a global implicit formulation as a curve in a 2D domain. With respect to [HGY+ 99] —in particular the analysis carried out using extremality coefficients, we provide a global rather than local analysis of the ridge curve. With respect to [Mor90,Mor96], we do not resort to local orientation procedures. Finally with respect to [Gra05], we are able to recover ridge colors —which are not algebraic invariants. The point of view of [Gra05] is to define invariants as functions of the fundamental forms and their derivatives. The equation of ridges is given in this setting. If one further specializes this equation for a surface given by a parametrisation, the result will be the equation we gave (up to a constant factor). The point of view of our approach is to work from the beginning on a parametrized surface. The definition of ridges envolves principal curvatures and principal directions of curvature which are independant of the given parametrization, but we explicit all these invariants wrt the parametrization and its derivatives. Hence we end with a polynomial with integer coefficients whose variables are the partial derivatives of the parametrization up to the third order. This polynomial is the same for any other parametrization. Certifying the approximation of P = 0. Second, exploiting as far as possible the geometry of P, we develop the first algorithm able to compute the ridges topology of a polynomial parametric surface. Our algorithm avoids the main difficulties of CAD methods: (i) singular and critical points are sequentially computed directly in 2D; (ii) no generic assumption is required, i.e. several critical or singular points may have the same horizontal projection; (iii) no computation with algebraic numbers is involved. Because algorithms based on the Cylindrical Algebraic Decomposition are not effective for our high degree curves such as P = 0, our algorithm is to the best of our knowledge the only one able to certify properties of the curve P = 0. The paper is organized as follows. Preliminary differential lemmas are proved in section 2. The implicit equation for ridges is derived in section 3. The algorithm to compute the topology of the ridge curve is described in section 4. Section 5 provides results of our algorithm on a degree 4 Bezier surface. 1.4 Notations Ridges and umbilics. At any non umbilical point of the surface, the maximal (minimal) principal curvature is denoted k1 (k2 ), and its associated direction d1 (d2 ). Anything related to the maximal (minimal) curvature is qualified blue (red), for example we shall 6 F. Cazals, J.C. Faugère, M. Pouget, F. Rouillier speak of the blue curvature for k1 or the red direction for d2 . Since we shall make precise statements about ridges, it should be recalled that, according to definition 1, umbilics are not ridge points. Differential calculus. For a bivariate function f (u, v), the partial derivatives are de3 f noted with indices, for example fuuv = ∂∂2 u∂v . The gradient of f is denoted f1 or df = (fu , fv ). The quadratic form induced by the second derivatives is denoted f2 (u, v) = 2 fuu u2 +2fuv uv+fvv v 2 . The discriminant of this form is denoted δ(f2 ) = fuv −fuu fvv . The cubic form induced by the third derivatives in denoted f3 (u, v) = fuuu u3 + 3fuuv u2 v + 3fuvv uv 2 + fvvv v 3 . The discriminant of this form is denoted δ(f3 ) = 2 2 4(fuuu fuvv − fuuv )(fuuv fvvv − fuvv ) − (fuuu fvvv − fuuv fuvv )2 . Let f be a real bivariate polynomial and F the real algebraic curve defined by f . A point (u, v) ∈ C2 is called – a singular point of F if f (u, v) = 0, fu (u, v) = 0 and fv (u, v) = 0; – a critical point of F if f (u, v) = 0, fu (u, v) = 0 and fv (u, v) 6= 0 (such a point has an horizontal tangent, we call it critical because if one fixes the v coordinate, then the restricted function is critical wrt the u coordinate, this notion will be usefull in section 4); – a regular point of F if f (u, v) = 0 and it is neither singular nor critical. If the domain D of study is a subset of R2 , one calls fiber a cross section of this domain at a given ordinate or abscissa. Misc. The inner product of two vectors x, y is denoted h x, y i, the norm of x is ||x|| = h x, x i1/2 and the exterior product is x ∧ext y. 2 Manipulations involving the Weingarten map of the surface Let Φ be the parameterization of class C k for k ≥ 4. Principal directions and curvatures of the surface are expressed in terms of second order derivatives of Φ. More precisely, the matrices of the first and second fundamental forms in the basis (Φu , Φv ) of the tangent space are     ef h Φu , Φu i h Φu , Φv i , I= = f g h Φu , Φv i h Φv , Φv i     N l m h Nn , Φuu i h Nn , Φuv i . II = = , with N = Φu ∧ext Φv , Nn = m n h Nn , Φuv i h Nn , Φvv i ||N || To compute the principal directions and curvatures, one resorts to the Weingarten map, whose matrix in the basis (Φu , Φv ) is given by W = (wij ) = I −1 II. The Weingarten map is a self-adjoint operator 1 of the tangent space [dC76]. The principal directions 1 A self-adjoint map L over a vector space V with a bilinear form < ., . > is a linear map such that h Lu, v i = h u, Lv i, for all u, v ∈ V . Such a map can be diagonalized in an orthonormal basis of V . Ridges and umbilics of polynomial parametric surfaces 7 di and principal curvatures ki are the eigenvectors and eigenvalues of the matrix W . Observing that ||N ||2 = det I, matrix W can also be written as follows:   1 AB W = . (3) (det I)3/2 C D This equation defines coefficients A, B, C and D as polynomials wrt the derivative of the parameterization Φ up to the second order, for example : gl − f m A = (det I)3/2 w11 = (det I)3/2 det I √ √ √ = det I(gh N/ det I, Φuu i − f h N/ det I, Φuv i) = gh N, Φuu i − f h N, Φuv i. As a general rule, in the following calculations, we will be interested to derive quantities which are polynomials wrt the derivatives of the parameterization. These calculations are based on quantities (principal curvatures and directions) which are independant of a given parameterization, hence the derived formula are valid for any parameterization. In order for normals and curvatures to be well defined, we assume the surface is regular i.e. det(I) 6= 0. 2.1 Principal curvatures. The characteristic polynomial of W is PW (k) = k 2 − tr(W )k + det(W ) = k 2 − (w11 + w22 )k + w11 w22 − w12 w21 . Its discriminant is ∆(k) = (tr(W ))2 − 4 det(W ) = (w11 − w22 )2 + 4w12 w21 . A simplification of this discriminant leads to the definition of the following function, denoted p2 : p2 = (det I)3 ∆(k) = (A − D)2 + 4BC The principal curvatures ki , with the convention k1 ≥ k2 , are the eigenvalues of W , that is: √ √ A + D − p2 A + D + p2 k = . (4) k1 = 2 2(det I)3/2 2(det I)3/2 A point is called an umbilic if the principal curvatures are equal. One has: Lemma 3. The two following equivalent conditions characterize umbilics: 1. p2 = 0 2. A = D ∧ B = C = 0. Proof. Since det(I) 6= 0, one has p2 = 0 ⇔ ∆(k) = 0, hence 1. characterizes umbilics. Condition 2. trivially implies 1. To prove the converse, assume that p2 = 0 i.e. the Weingarten map has a single eigenvalue k. This linear map is self-adjoint hence diagonalizable in an orthogonal basis, and the diagonal form is a multiple of the identity. It is easily checked that the matrix remains a multiple of the identity in any basis of the tangent space, in particular in the basis (Φu , Φv ), which implies condition 2. 8 F. Cazals, J.C. Faugère, M. Pouget, F. Rouillier 2.2 Principal directions. Let us focus on the maximum principal direction d1 . A vector of direction d1 is an eigenvector of W for the eigenvalue k1 . Denote: W − k1 Id =   w11 − k1 w12 . w21 w22 − k1 (5) At non umbilic points, the matrix W − k1 Id has rank one, hence either (−w12 , w11 − k1 ) 6= (0, 0) or (−w22 + k1 , w21 ) 6= (0, 0). Using the expression of W given by Eq. 3, up to a normalization factor of (det I)3/2 , a non zero maximal principal vector can be chosen as either √ v1 = 2(det I)3/2 (−w12 , w11 − k1 ) = (−2B, A − D − p2 ) √ or w1 = 2(det I)3/2 (−w22 + k1 , w21 ) = (A − D + p2 , 2C). For the minimal principal direction d2 one chooses v2 = (−2B, A − D + √ w2 = (A − D − p2 , 2C). √ p2 ) and Lemma 4. One has the following relations: v1 = (0, 0) ⇔ (B = 0 ∧ A ≥ D), v2 = (0, 0) ⇔ (B = 0 ∧ A ≤ D), w1 = (0, 0) ⇔ (C = 0 ∧ A ≤ D), w2 = (0, 0) ⇔ (C = 0 ∧ A ≥ D). Proof. The proofs being equivalent, we focus on the first one. One has: ( B=0 p v1 = (0, 0) ⇔ A − D = (A − D)2 + 4BC ( B=0 ⇔ A−D ≥0 ( B=0 p ⇔ A − D = (A − D)2 A direct consequence of lemma 4 is the following: Observation. 1 1. The two vector fields v1 and w1 vanish simultaneously exactly at umbilics. The same holds for v2 and w2 . 2. The equation {v1 = (0, 0) ∨ v2 = (0, 0)} is equivalent to B = 0. 3 Implicitly defining ridges In this section, we prove the implicit expression P = 0 of ridges where P is a polynomial wrt the derivatives of the parameterization up to the third order. Before diving into the technicalities, we first outline the method. Ridges and umbilics of polynomial parametric surfaces 9 3.1 Method outline Principal curvatures and directions read from the Weingarten map of the surface. At each point which is not an umbilic, we have just defined two vector fields v1 or w1 which are collinear with d1 , with the additional property that one (at least) of these two vectors is non vanishing. Let z stand for one of these non vanishing vectors and denote y1 one of the two opposite unit vectors orienting d1 . The nullity of b0 = h dk1 , y1 i is equivalent to that of h dk1 , z i —that is the normalization of the vector along which the directional derivative is computed does not matter. Using v1 and w1 we obtain two independent equations of blue ridges. Each has the drawback of encoding, in addition to blue ridge points, the points where v1 (or w1 ) vanishes. As a consequence of observation 1, the conjunction of these two equations defines the set of blue ridges union the set of umbilics. The same holds for red ridges and the minimal principal vector fields v2 ad w2 . One has to note the symmetry between the equations for blue and red ridges as stated in lemma 5. Eventually, combining the equation for blue ridges with v1 and the equation for red ridges with v2 gives the set of blue ridges union the set of red ridges union the set of zeros of v1 = 0 or v2 = 0. This last set is also B = 0 (observation 1), hence dividing by B allows to eradicate these spurious points and yields the equation P = 0 of blue ridges together with red ridges (see theorem 8). One can think of this equation as an improved version of the Gaussian extremality Eg = b0 b3 defined in [Thi96]. Our strategy cumulates several advantages: (i)blue and red ridges are processed at once, and the information is encoded in a single equation (ii)orientation issues arising when one is tracking the zero crossings of b0 or b3 disappear. The only drawback is that one looses the color of the ridge. But this color is recovered with the evaluation of the sign of factors of the expression P . 3.2 Precisions of vocabulary In the statement of the results, we shall use the following terminology. Umbilic points are points where both principal curvatures are equal. A ridge point is a point which is not an umbilic, and is an extremum of a principal curvature along its curvature line. A ridge point is further called a blue (red) ridge point if it is an extremum of the blue (red) curvature along its line. A ridge point may be both blue and red, in which case it is called a purple point. 3.3 Implicit equation of ridges Lemma 5. For a regular surface, there exists functions a, a′ , b, b′ which are polynomials wrt A, B, C, D, det I as well as their first derivatives and hence are polynomials wrt the derivatives of the parameterization up to the third order, such that: √ 1. the set of blue ridges union {v1 = 0} has equation a p2 + b = 0, √ 2. the set of blue ridges union {w1 = 0} has equation a′ p2 +(b′ = 0, √ a p2 + b = 0 3. the set of blue ridges union the set of umbilics has equation √ a′ p2 + b′ = 0 10 F. Cazals, J.C. Faugère, M. Pouget, F. Rouillier √ 4. the set of red ridges union {v2 = 0} has equation a p2 − b = 0, √ 5. the set of red ridges union {w2 = 0} has equation a′ p2 −(b′ = 0, √ a p2 − b = 0 6. the set of red ridges union the set of umbilics has equation √ a′ p2 − b′ = 0 Moreover, a, a′ , b, b′ are defined by the equations: √ √ √ √ a p2 +b = p2 (det I)5/2 h dk1 , v1 i ; a′ p2 +b′ = p2 (det I)5/2 h dk1 , w1 i. (6) Proof. The principal curvatures are not differentiable at umbilics, hence the equation √ h dk1 , v1 i is not defined at umbilics. On the other hand, since p2 only vanishes at um√ bilics and det I never vanishes, the equation p2 (det I)5/2 h dk1 , v1 i defines the same √ set. In addition, this equation is also well defined at umbilics. Rewritten as a p2 + b, the functions a and b are then polynomial functions of the derivatives of the parameterization Φ up to the third order. The explicit expressions of a and b are given in [CFPR05b]. This equation describes the set of blue ridge points union the set where v1 vanishes. A similar derivation yields the second claim. Finally, the third claim follows from observation 1. Results for red ridges are similar and the reader is referred to [CFPR05b] for the details. Lemma 6. 1. If p2 = 0 then a = b = a′ ( = b′ = 0. a = b = a′ = b′ = 0 2. The set of purple points has equation p2 6= 0 Proof. 1. If p2 = 0, one has A = D and B = C = 0. Substituting these conditions in the expressions of a, a′ , b, b′ gives the result, computations are sketched in [CFPR05b]. 2. Let p be a purple point, it is a ridge point and hence not an umbilic, then p2 6= 0. The point p is a blue and a red ridge point, hence it satisfies all equations of lemma 5. √ √ If a 6= 0 then equations 1. and 4. imply p2 = −b/a = b/a hence b = 0 and p2 = 0 which is a contradiction. Consequently, a = 0 and again equation 1. implies b = 0. A similar argument with equation 2. and 5. gives a′ = b′ = 0. The converse is trivial: if a = b = a′ = b′ = 0 then equations 3. and 6. imply that the point is a purple point or an umbilic. The additional condition p2 6= 0 excludes umbilics. The following definition is a technical tool to state the next theorem in a simple way. The meaning of the function Signridge introduced here will be clear from the proof of the theorem. Essentially, this function describes all the possible sign configurations for ab and a′ b′ at a ridge point. Definition 7. The function Signridge takes the values ( ( ab < 0 ab ≤ 0 -1 if or , ′ ′ ab ≤0 a′ b′ < 0 ( ( ab ≥ 0 ab > 0 or , +1 if a′ b′ > 0 a′ b′ ≥ 0 Ridges and umbilics of polynomial parametric surfaces 11 0 if ab = a′ b′ = 0. Theorem 8. The set of blue ridges union the set of red ridges union the set of umbilics has equation P = 0 where P = (a2 p2 − b2 )/B is a polynomial wrt A, B, C, D, det I as well as their first derivatives and hence is a polynomial wrt the derivatives of the parameterization up to the third order. For a point of this set P, one has: – If p2 = 0, the point is an umbilic. – If p2 6= 0 then: • if Signridge = −1 then the point is a blue ridge point, • if Signridge = +1 then the point is a red ridge point, • if Signridge = 0 then the point is a purple point. Proof. To form the equation of P, following the characterization of red and blue ridges in lemma 5, and the vanishing of the vector fields v1 and v2 in lemma 4, we take the product of equations 1. and 4. of lemma 5. This gives an equation encoding ridges and the zeros of v1 or v2 which are precisely those of B. Dividing by B allows to eradicates these points and yields the equation of ridges and umbilics. One can also use equations 2. and 5. of lemma 5 and end with the expression P = −(a′2 p2 − b′2 )/C. The equivalence between the two equations of P and another alternative one P = 2(a′ b − ab′ ) is proved formally with the help of Maple, see [CFPR05b]. To qualify points on P, first observe that the case p2 = 0 has already been considered in lemma 6. Therefore, assume p2 6= 0, and first notice the following two simple facts: – The equation (a2 p2 − b2 )/B = 0 for P implies that a = 0 ⇔ b = 0 ⇔ ab = 0. Similarly, the equation −(a′2 p2 − b′2 )/C = 0 for P implies that a′ = 0 ⇔ b′ = 0 ⇔ a′ b′ = 0. – If ab 6= 0 and a′ b′ 6= 0, the equation ab′ − a′ b = 0 for P implies b/a = b′ /a′ , that is the signs of ab and a′ b′ agree. These two facts explain the introduction of the function Signridge of definition 7. This function enumerates all disjoint possible configurations of signs for ab and a′ b′ for a point on P. One can now study the different cases w.r.t. the signs of ab and a′ b′ or equivalently the values of the function Signridge . Assume Signridge = −1. √ √ –First case: ab < 0. The equation (a2 p2 −b2 )/2B = 0 implies that (a p2 +b)(a p2 − √ √ b) = 0. Since p2 > 0, one must have a p2 + b = 0 which is equation 1 of lemma 5. From the second simple fact, either a′ b′ < 0 or a′ b′ = 0. √ – For the first sub-case a′ b′ < 0, equation (a′2 p2 − b′2 )/C = 0 implies (a′ p2 + √ ′ ′√ ′ ′ ′√ b )(a p2 − b ) = 0. Since p2 > 0, one must have a p2 + b = 0 which is equation 2 of lemma 5. – For the second sub-case a′ b′ = 0, one has a′ = b′ = 0 and the equation 2 of lemma 5 is also satisfied. (Moreover, equation 5 is also satisfied which implies that w2 = 0). 12 F. Cazals, J.C. Faugère, M. Pouget, F. Rouillier In both cases, equations 1 and 2 or equivalently equation 3 are satisfied. Since one has excluded umbilics, the point is a blue ridge point. –Second case: ab = 0. One has a′ b′ < 0 the study is similar to the above. ab = 0 implies equation 1 and a′ b′ < 0 implies equation 2 of lemma 5. The point is a blue ridge point. Assume Signridge = 1. This case is the exact symmetric of the previous, one only has to exchange the roles of a, b and a′ , b′ . Assume Signridge = 0. The first simple fact implies a = b = a′ = b′ = 0 and lemma 6 identifies a purple point. As shown along the proof, the conjunctions <, = and =, < in the definition of Signridge = −1 correspond to the blue ridge points where the vector fields w2 and v2 vanish. The same holds for Signridge = 1 and w1 and v1 . One can also observe that the basic ingredient of the previous proof is to transform an equation with a square root into a system with an inequality. More formally: Observation. 2 For x, y, z real numbers and z ≥ 0, one has: ( √ x2 z − y 2 = 0 x z + y = 0 ⇐⇒ xy ≤ 0 (7) 3.4 Singular points of P Having characterized umbilics, purple points and ridges in the domain D with implicit equations, an interesting question is to relate the properties of these equations to the classical differential geometric properties of these points. In particular, recall that generically (with the description of surfaces with Monge patches and contact theory [HGY+ 99]), umbilics of a surface are either 1-ridge umbilics or 3-ridge umbilics. This means that there are either 1 or 3 non-singular ridge branches passing through an umbilic. The later are obviously singular points of P since three branches of the curve are crossing at the umbilic. For the former ones, it is appealing to believe they are regular points since the tangent space to the ridge curve on the surface at such points is well defined and can be derived from the cubic of the Monge form [HGY+ 99]. Unfortunately, one has: Proposition 9. Umbilics are singular points of multiplicity at least 3 of the function P (i.e. the gradient and the Hessian of P vanish). Proof. Following the notations of Porteous [Por01], denote Pk , k = 1, . . . , 3 the kth times linear form associated with P , that is Pk = [∂P/(∂uk−i ∂v i )]i=0,...,k . Phrased differently, P1 is the gradient, P2 is the vector whose three entries encodes the Hessian of P , etc. To show that the multiplicity of an umbilic of coordinates (u0 , v0 ) is at least three, we need to show that P1 (u0 , v0 ) = [0, 0], P2 (u0 , v0 ) = [0, 0, 0]. We Ridges and umbilics of polynomial parametric surfaces 13 naturally do not know the coordinates of umbilics, but lemma 3 provides the umbilical conditions. The proof consists of computing derivatives and performing the appropriate substitutions under Maple, and is given in [CFPR05b]. We can go one step further so as to relate the type of the cubic P3 —the third derivative of P — to the number of non-singular ridge branches at the umbilic. Classical arguments of singularity theory give the following classification (the proof can be found in [CFPR05b]). Proposition 10. The classification of an umbilic as 1-ridge or 3-ridges from P3 goes as follows: – If P3 is elliptic, that is the discriminant of P3 is positive (δ(P3 ) > 0), then the umbilic is a 3-ridge umbilic and the 3 tangent lines to the ridges at the umbilic are distinct. – If P3 is hyperbolic (δ(P3 ) < 0) then the umbilic is a 1-ridge umbilic. Note that the classifications of umbilics with the Monge cubic CM and the cubic P3 do not coincide. The Monge cubic CM is the third order part of the Taylor expansion of the surface parametrized by its tangent plane. Indeed if CM is elliptic, it may occur that two ridges have the same tangent. In such a case, the cubic P3 is not elliptic since δ(P3 ) = 0. Since purple points correspond to the intersection of two ridges, one has: Proposition 11. Purple points are singular points of multiplicity at least 2 of the function P (i.e. the gradient of P vanishes). Proof. It follows from the equation P = 2(a′ b − ab′ ) that dP = 2(d(a′ )b + a′ d(b) − d(a)b − ad(b)). At purple points one has a = a′ = b = b′ = 0 hence dP = 0. 3.5 Polynomial surfaces A fundamental class of surface used in Computer Aided Geometric Design consists of Bezier surfaces and splines. In this section, we state some elementary observations on the objects studied so far, for the particular case of polynomial parametric surfaces. Notice that the parameterization can be general, in which case Φ(u, v) = (x(u, v), y(u, v), z(u, v)), or can be a height function Φ(u, v) = (u, v, z(u, v)). We first observe that if Φ is a polynomial, all its derivatives are also polynomials. Thus in the polynomial case the equation of ridges, which is a polynomial wrt to these derivatives, is algebraic. Hence the set of all ridges and umbilics is globally described by an algebraic curve. As a corollary of Thm. 8, one can give upper bounds for the total degree of the polynomial P wrt that of the parameterization. Distinguishing the cases where Φ is a general parameterization or a height function (that is Φ(u, v) = (u, v, h(u, v))) with h(u, v) and denoting d the total degree of Φ, P has total degree 33d − 40 or 15d − 22 for a height function. Note that in the case of a height function, P is divided by its factor (det I)2 (cf. [CFPR05b]). 14 F. Cazals, J.C. Faugère, M. Pouget, F. Rouillier In the more general case where the parameterization is given by rational fractions of polynomials, P is a rational function of the surface parameters too. The denominator of P codes the points where the surface is not defined and away from these points, the numerator codes the ridges and umbilics. 4 Certified topological approximation In this section, we circumvent the difficulties of the Cylindrical Algebraic Decomposition (CAD) and develop a certified algorithm to compute the topology of P. Consider a parameterized surface Φ(u, v), the parameterization being polynomial with rational coefficients. Let P be the curve encoding the ridges of Φ(u, v). We aim at studying P on the compact box domain D = [a, b] × [c, d]. Given a real algebraic curve, the standard way to approximate it consists of resorting to the CAD. Running the CAD requires computing singular points and critical points of the curve —points with a horizontal tangent. Theoretically, these points are defined by zero-dimensional systems. Practically, because of the high degree of the polynomials involved, the calculations may not go through. Replacing the bottlenecks of the CAD by a resolution method adapted to the singular structure of P, we develop an algorithm producing a graph G embedded in the domain D, which is isotopic to the curve P of ridges in D. Key points are that: 1. no generic assumption is required, i.e. several critical or singular points may have the same horizontal projection; 2. no computation with algebraic numbers is involved. 4.1 Algebraic tools Two algebraic methods are ubiquitously called by our algorithm: univariate root isolation and rational univariate representation. We briefly present these tools and give references for the interested reader. Univariate root isolation. This tool enables to isolate roots of univariate polynomials whose coefficients are rational numbers, by means of intervals with rational bounds. The method uses the Descartes rule and is fully explained in [RZ03]. Rational univariate representation [Rou99]. The Rational Univariate Representation is, with the end-user point of view, the simplest way for representing symbolically the roots of a zero-dimensional system without loosing information (multiplicities or real roots) since one can get all the information on the roots of the system by solving univariate polynomials. Given a zero-dimensional system I =< p1 , . . . , ps > where the pi ∈ Q[X1 , . . . , Xn ], a Rational Univariate Representation of V(I) has the following shape : ft (T ) = 0, X1 = gt,X1 (T ) gt,Xn (T ) gt,1 (T ) , . . . , Xn = gt,1 (T ) , where ft , gt,1 , gt,X1 , . . . , gt,Xn ∈ Q[T ] (T is a new variable). It is uniquely defined w.r.t. a given polynomial t which separates V (I) (injective on V (I)), the polynomial ft being necessarily the characteristic polynomial of mt in Q[X1 , . . . , Xn ]/I. The RUR defines a bijection between the roots of I and those of ft preserving the multiplicities and the real roots : Ridges and umbilics of polynomial parametric surfaces 15 V(I)(∩R) ≈ V(ft )(∩R) α = (α1 , . . . , αn ) → t(α) g (t(α)) gt,Xn (t(α)) 1 , . . . , ) ← t(α) ( gt,X gt,1 (t(α)) t,1 (t(α)) The RUR also enables efficient evaluation of the sign of polynomials at the roots of a system. 4.2 Assumptions on the ridge curve and study points Following the theoretical study performed in section 3.4, the only assumption made is that the surface admits generic ridges in the sense that real singularities of P satisfy the following conditions: – Real singularities of P are of multiplicity at most 3. – Real singularities of multiplicity 2 are called purple points. They satisfy the system Sp = {a = b = a′ = b′ = 0, δ(P2 ) > 0, p2 6= 0}. In addition, this implies that two real branches of P are passing through a purple point. – Real singularities of multiplicity 3 are called umbilics and they satisfy the system Su = {p2 = 0} = {p2 = 0, P = 0, Pu = 0, Pv = 0}. In addition, if δ(P3 ) denote the discriminant of the cubic of the third derivatives of P at an umbilic, one has: • if δ(P3 ) > 0, then the umbilic is called a 3-ridge umbilic and three real branches of P are passing through the umbilic with three distinct tangents; • if δ(P3 ) < 0, then the umbilic is called a 1-ridge umbilic and one real branch of P is passing through the umbilic. As we shall see in section 4.5, these conditions are checked during the processing of the algorithm. Given this structure of singular points, the algorithm successively isolate umbilics, purple points and critical points. As a system defining one set of these points also includes the points of the previous system, we use a localization method to simplify the calculations. The points reported at each stage are characterized as roots of a zerodimensional system —a system with a finite number of complex solutions, together with the number of half-branches of the curve connected to each point. In addition, points on the border of the domain of study need a special care. This setting leads to the definition of study points: Definition 12. Study points are points in D which are – real singularities of P, that is Ss = Su ∪ Sp , with Su = S1R ∪ S3R and • S1R = {p2 = P = Pu = Pv = 0, δ(P3 ) < 0} • S3R = {p2 = P = Pu = Pv = 0, δ(P3 ) > 0} • Sp = {a = b = a′ = b′ = 0, δ(P2 ) > 0, p2 6= 0} = {a = b = a′ = b′ = 0, δ(P2 ) > 0} \ Su – real critical points of P in the v-direction (i.e. points with a horizontal tangent which are not singularities of P) defined by the system Sc = {P = Pu = 0, Pv 6= 0}; – intersections of P with the left and right sides of the box D satisfying the system Sb = {P (a, v) = 0, v ∈ [c, d]} ∪ {P (b, v) = 0, v ∈ [c, d]}. Such a point may also be critical or singular. 16 F. Cazals, J.C. Faugère, M. Pouget, F. Rouillier 4.3 Output specification Definition 13. Let G be a graph whose vertices are points of D and edges are nonintersecting straight line-segments between vertices. Let the topology on G be induced by that of D. We say that G is a topological approximation of the ridge curve P on the domain D if G is ambient isotopic to P ∩ D in D. More formally, there exists a function F : D × [0, 1] −→ D such that: – F is continuous; – ∀t ∈ [0, 1], Ft = F (., t) is an homeomorphism of D onto itself; – F0 = IdD and F1 (P ∩ D) = G. Note that homeomorphic approximation is weaker and our algorithm actually gives isotopy. In addition, our construction allows to identify singularities of P to a subset of vertices of G while controlling the error on the geometric positions. We can also color edges of G with the color of the ridge curve it is isotopic to. Once this topological sketch is given, one can easily compute a more accurate geometrical picture. 4.4 Method outline Taking the square free part of P , we can assume P is square free. We can also assume P has no part which is a horizontal segment —parallel to the u-axis. Otherwise this means that a whole horizontal line is a component of P . In other words, the content of P wrt u is a polynomial in v and we can study this factor separately and divide P by this factor. Eventually, to get the whole topology of the curve, one has to merge the components. Our algorithms consists of the following five stages: 1. Isolating study points. Study point are isolated in 2D with rational univariate representations (RUR). Study points within a common fiber are identified. 2. Regularization of the study boxes. We know the number of branches of the curve going through each study point. The boxes of study points are reduced so as to be able to define the number of branches coming from the bottom and from the top. 3. Computing regular points in study fibers. In each fiber of a study point, the ucoordinates of intersection points with P other than study points are computed. 4. Adding intermediate rational fibers. Add rational fibers between study points fibers and isolate the u-coordinates of intersection points with P. 5. Performing connections. This information is enough to perform the connections. Consider the cylinder between two consecutive fibers, the number of branches connected from above the lower fiber is the same than the number of branches connected from below the higher fiber. Hence there is only one way to perform connections with non-intersecting straight segments. 4.5 Step 1. Isolating study points The method to identify these study points is to compute a RUR of the system defining them. More precisely, we sequentially solve the following systems: Ridges and umbilics of polynomial parametric surfaces 17 1. The system Su from which the sets S1R and S3R are distinguished by evaluating the sign of δ(P3 ). 2. The system Sp for purple points. 3. The system Sc for critical points. 4. The system Sb for border points, that is intersections of P with the left and right sides of the box D. Solving this system together with one of the previous identifies border points which are also singular or critical. Selecting only points belonging to D reduces to adding inequalities to the systems and is well managed by the RUR. According to [Rou99], solving such systems is equivalent to solving zero-dimensional systems without inequalities when the number of inequations remains small compared to the number of variables. The RUR of the study points provides a way to compute a box around each study point qi which is a product of two intervals [u1i ; u2i ] × [vi1 ; vi2 ]. The intervals can be as small as desired. Until now, we only have separate informations on the different systems. In order to identify study points having the same v-coordinate, we need to cross these informations. First we compute isolation intervals for all the v-coordinates of all the study points together, denote I this list of intervals. If two study points with the same v-coordinate are solutions of two different systems, the gcd of polynomials enable to identify them: – Initialize the list I with all the isolation intervals of all the v-coordinates of the different systems. – Let A and B be the square free polynomials defining the v-coordinates of two different systems, and IA , IB the lists of isolation intervals of their roots. Let C = gcd(A, B) and IC the list of isolation intervals of its roots. One can refine the elements of IC until they intersect only one element of IA and one element of IB . Then replace these two intervals in I by the single interval which is the intersection of the three intervals. Do the same for every pair of systems. – I then contains intervals defining different real numbers in one-to-one correspondence with the v-coordinates of the study points. It remains to refine these intervals until they are all disjoint. Second, we compare the intervals of I and those of the 2d boxes of the study points. Let two study points qi and qj be represented by [u1i ; u2i ]×[vi1 ; vi2 ] and [u1j ; u2j ]×[vj1 ; vj2 ] with [vi1 ; vi2 ] ∩ [vj1 ; vj2 ] 6= ∅. One cannot, a priori, decide if these two points have the same v-coordinate or if a refinement of the boxes will end with disjoint v-intervals. On the other hand, with the list I, such a decision is straightforward. The boxes of the study points are refined until each [vi1 ; vi2 ] intersects only one interval [wi1 ; wi2 ] of the list I. Then two study points intersecting the same interval [wi1 ; wi2 ] are in the same fiber. Finally, one can refine the u-coordinates of the study points with the same v coordinate until they are represented with disjoint intervals since, thanks to localizations, all the computed points are distinct. Checking genericity conditions of section 4.2. First, real singularities shall be the union of purple and umbilical points, this reduces to compare the systems for singular points and for purple and umbilical points. Second, showing that δ(P3 ) 6= 0 for umbilics and δ(P2 ) > 0 for purple points reduces to sign evaluation of polynomials at the roots of a system (see section 4.1). 18 F. Cazals, J.C. Faugère, M. Pouget, F. Rouillier v Number of branches above n+ i,j vi2 αi Number of branches below − ni,j vi1 1 βi,1 2 βi,1 u1i,1 u2i,1 1 βi,l i 2 βi,l i u1i,mi u2i,mi u 1(2) Fig. 3: Notations for a fiber involving several critical/singular points: ui,j are used for 1(2) study points, βi,j for simple points. v αi+1 δi vi2 αi vi1 u Fig. 4: Performing connections between the study point fiber αi and the intermediate fiber δi 4.6 Step 2. Regularization of the study boxes At this stage, we have computed isolating boxes of all study points {qi,j , i = 1 . . . s, j = 1 . . . mi } : the v-coordinates α1 , ..., αs are isolated by intervals [vi1 ; vi2 ], i = 1..s and the u-coordinates of the mi study points in each fiber αi are isolated by intervals [u1i,j ; u2i,j ], j = 1..mi . We know the number of branches of the curve passing through each study point : it is 6 for a 3-ridge umbilic, 4 for a purple and 2 for others. We want to compute the number branches coming from the bottom and from the top. We first reduce the box Ridges and umbilics of polynomial parametric surfaces 19 until the number of intersections between the curve and the border of the box matches the known number of branches connected to the study point. Then the intersections are obviously in one-to-one correspondance with the branches. Second, as in [SW05] for example, we reduce the hight of the box again if necessary so that intersections only occur on the top or the bottom of the box. Compting the number of intersections reduces to solve 4 univariate polynomials with rational coefficients. Reducing a box means refining its representation with the RUR. 4.7 Step 3. Computing regular points in study fibers We now compute the regular points in each fiber P (u, αi ) = 0. Computing the regular points of each fiber is now equivalent to computing the roots of the polynomials P (u, αi ) outside the intervals representing the u-coordinates of the study points (which contain all the multiple roots of P (u, αi )). Denote [u1i,j ; u2i,j ], j = 1..mi the intervals representing the u-coordinates of the study points on the fiber of αi and [vi1 , vi2 ] an interval containing (strictly) αi and no other αj , j 6= i. Substituting v by any rational value q ∈ [vi1 , vi2 ] in P (u, v) gives a univariate polynomial with rational coefficients P (u, q). We then isolate the (simple) 2 1 i roots of this polynomial P (u, q) on the domain [a, b] \ ∪m j=1 [ui,j ; ui,j ] : the algorithm 2 1 returns intervals [βi,j ; βi,j ], j = 1 . . . li representing these roots. To summurize the information up to this point : we have, along each fiber, a collection of points si,j , i = 1 . . . s, j = 1, . . . , mi + li , which are either study points or regular points of P. Each such point is isolated in a box i.e. a product of intervals and − comes with two integers (n+ i,j , ni,j ) denoting the number of branches in D connected from above and from below. 4.8 Step 4. Adding intermediate rational fibers Consider now an intermediate fiber, i.e. a fiber associated with v = δi i = 1 . . . s − 1, with δi a rational number in-between the intervals of isolation of two consecutive values αi and αi+1 . If the fibers v = c or v = d are not fibers of study points, then they are added as fibers δ0 or δs . Getting the structure of such fibers amounts to solving a univariate polynomial with rational coefficients, which is done using the algorithm described in section 4.1. Thus, each such fiber also comes with a collection of points, isolated in boxes, for which one − knows that n+ i,j = ni,j = 1. 4.9 Step 5. Performing connections We thus obtain a full and certified description of the fibers: all the intersection points with P and their number of branches connected. We know, by construction, that the branches of P between fibers have empty intersection. The number of branches connected from above a fiber is the same than the number of branches connected from 20 F. Cazals, J.C. Faugère, M. Pouget, F. Rouillier below the next fiber. Hence there is only one way to perform connections with nonintersecting straight segments. More precisely, vertices of the graph are the centers of isolation boxes, and edges are line-segments joining them. Notice that using the intermediate fibers v = δi is compulsory if one wishes to get a graph G isotopic to P. If not, whenever two branches have common starting points and endpoints, the embedding of the graph G obtained is not valid since two arcs are identified. The algorithm is illustrated on Fig. 4. In addition – If a singular point box have width δ, then the distance between the singular point and the vertex representing it is less than δ. – One can compute the sign of the function Signridge defined in 4.2 for each regular point of each intermediate fiber. This defines the color of the ridge branch it belongs to. Then one can assign to each edge of the graph the color of its end point which is on an intermediate fiber. 5 Illustration We provide the topology of ridges for two Bezier surfaces defined over the domain D = [0, 1] × [0, 1]. The first surface has control points   [0, 0, 0] [1/4, 0, 0] [2/4, 0, 0] [3/4, 0, 0] [4/4, 0, 0] [0, 1/4, 0] [1/4, 1/4, 1] [2/4, 1/4, −1] [3/4, 1/4, −1] [4/4, 1/4, 0]   [0, 2/4, 0] [1/4, 2/4, −1] [2/4, 2/4, 1] [3/4, 2/4, 1] [4/4, 2/4, 0]   [0, 3/4, 0] [1/4, 3/4, 1] [2/4, 3/4, −1] [3/4, 3/4, 1] [4/4, 3/4, 0] [0, 4/4, 0] [1/4, 4/4, 0] [2/4, 4/4, 0] [3/4, 4/4, 0] [4/4, 4/4, 0] Alternatively, this surface can be expressed as the graph of the total degree 8 polynomial h(u, v) for (u, v) ∈ [0, 1]2 : h(u, v) = 116u4 v 4 − 200u4 v 3 + 108u4 v 2 − 24u4 v − 312u3 v 4 + 592u3 v 3 − 360u3 v 2 + 80u3 v + 252u2 v 4 − 504u2 v 3 + 324u2 v 2 − 72u2 v − 56uv 4 + 112uv 3 − 72uv 2 + 16uv. The computation of the implicit curve has been performed using Maple 9.5 and requires less than one minute (see [CFPR05b]). It is a bivariate polynomial P (u, v) of total degree 84, of degree 43 in u, degree 43 in v with 1907 terms and coefficients with up to 53 digits. Figure 5 displays the topological approximation graph of the ridge curve in the parametric domain D computed with the algorithm of section 4. There are 19 critical points (black dots), 17 purple points (pink dots) and 8 umbilics, 3 of which are 3-ridge (green) and 5 are 1-ridge (yellow). We have computed the subsets Su , Sp and Sc by using the software FG B and RS (http://fgbrs.lip6.fr). The RUR can be computed as shown in [Rou99] or alternatively, Gröbner basis can be computed first using [Fau99] or [Fau02]. We tested both methods and the computation time for the biggest system Sc does not exceed 10 minutes with a Pentium M 1.6 Ghz. The following table gives the main characteristics of these systems : Ridges and umbilics of polynomial parametric surfaces 21 System # of roots ∈ C # of roots ∈ R # of real roots ∈ D Su 160 16 8 Sp 749 47 17 Sc 1432 44 19 In order to have more insight of the geometric meaning of the ridge curve, the surface and its ridges are displayed on Fig. 6. This plot is computed without topological certification with the rs tci points function (from RS software, see also [CFPR05a]) from the polynomial P and then lifted on the surface. Fig. 5: Bi-quartic Bézier example : isotopic approximation of the ridge curve with 3-ridge umbilics (green), 1-ridge umbilics (yellow), purple points (pink) and critical points (points with an horizontal tangent displayed in black). 22 F. Cazals, J.C. Faugère, M. Pouget, F. Rouillier Fig. 6: Plot of the bi-quartic Bezier surface with ridges and umbilics (boxed in green) The second surface is a bi-quadratic Bézier of equation Φ (u, v) = [2/3 v + 2/3 uv − 1/3 u2 v + 1/3 v 2 − 2/3 v 2 u + 1/3 u2 v 2 , 1/2 u + 1/2 u2 + uv − u2 v − 1/2 v 2 u + 1/2 u2 v 2 , 1 + 3 v 2 − u − 4 v + 5 uv + u2 v − 7/2 v 2 u − 5/2 u2 v 2 ] The ridge curve has total degree 56 and partial degrees 33 with 1078 terms and coefficients with up to 15 digits. The computation of the biggest system of study points Sc takes 4.5 minutes. On this example, study point boxes have to be refined up to a size of less than 2−255 to compute the topology. The following table gives the main characteristics of the study point systems : System # of roots ∈ C # of real roots ∈ D Su 70 1 Sp 293 6 Sc 695 5 Ridges and umbilics of polynomial parametric surfaces 23 Figure 7 diplays the topology of the ridges. In addition to study points, the regular points of all fibers are displayed as small black dots. Fig. 7: Bi-quadratic Bézier example : isotopic approximation of the ridge curve with the same color coding as in Fig. 5, in addition small black dots are regular points of all fibers. 6 Conclusion This paper sets the implicit equation P = 0 of the singular curve encoding all ridges of a smooth parametric surface. For parametric algebraic surfaces, we developed an algorithm to report a topologically certified approximation of the ridges. This algorithm is computationally demanding in terms of algebra. It is in a sense complementary to the heuristic one developed in a companion paper [CP05], which is working directly on a triangulation of the surface, and provide a fast way to report non certified results. 24 F. Cazals, J.C. Faugère, M. Pouget, F. Rouillier The method developed for the computation of the topology of the ridges can be generalized for other algebraic curves. It gives an alternative to usual algorithms based on the CAD provided one knows the geometry of curve branches at singularities. Acknowledgements F. Cazals and M. Pouget acknowledge the support of the AIM@Shape and ACS European projects. Jean-Pierre Merlet is acknowledged for fruitful discussions. References [ABL03] D. Attali, J.-D. Boissonnat, and A. Lieutier. Complexity of the delaunay triangulation of points on surfaces the smooth case. In ACM SoCG, San Diego, 2003. [CFPR05a] F. Cazals, J.-C. Faugère, M. Pouget, and F. Rouillier. Topologically certified approximation of umbilics and ridges on polynomial parametric surface. Technical Report 5674, INRIA, 2005. [CFPR05b] F. Cazals, J.-C. Faugères, M. Pouget, and F. Rouillier. The implicit structure of ridges of a smooth parametric surface. Technical Report 5608, INRIA, 2005. To appear in CAGD. [CP05] F. Cazals and M. Pouget. Topology driven algorithms for ridge extraction on meshes. Technical Report RR-5526, INRIA, 2005. [dC76] M. de Carmo. Differential Geometry of Curves and Surfaces. Prentice Hall, Englewood Cliffs, NJ, 1976. [Fau99] J.-C. Faugère. A new efficient algorithm for computing gröbner bases (f4 ). Journal of Pure and Applied Algebra, 139(1-3):61–88, June 1999. [Fau02] Jean-Charles Faugère. A new efficient algorithm for computing gröbner bases without reduction to zero f5 . In International Symposium on Symbolic and Algebraic Computation Symposium - ISSAC 2002, Villeneuve d’Ascq, France, Jul 2002. [GLMT04] G. Gatellier, A. Labrouzy, B. Mourrain, and J.-P. Técourt. Computing the topology of 3-dimensional algebraic curves. In Computational Methods for Algebraic Spline Surfaces, pages 27–44. Springer-Verlag, 2004. [Gra05] J. Gravesen. Third order invariants of surfaces. In T. Dokken and B. Juttler, editors, Computational methods for algebraic spline surfaces. Springer, 2005. [GVN02] L. Gonzalez-Vega and I. Necula. Efficient topology determination of implicitly defined algebraic plane curves. Computer Aided Geometric Design, 19(9), 2002. [HCV52] D. Hilbert and S. Cohn-Vossen. Geometry and the Imagination. Chelsea, 1952. [HGY+ 99] P. W. Hallinan, G. Gordon, A.L. Yuille, P. Giblin, and D. Mumford. Two-and ThreeDimensional Patterns of the Face. A.K.Peters, 1999. [Koe90] J.J. Koenderink. Solid Shape. MIT, 1990. [Mor90] R. Morris. Symmetry of Curves and the Geometry of Surfaces: two Explorations with the aid of Computer Graphics. Phd Thesis, 1990. [Mor96] R. Morris. The sub-parabolic lines of a surface. In Glen Mullineux, editor, Mathematics of Surfaces VI, IMA new series 58, pages 79–102. Clarendon Press, Oxford, 1996. [MWP96] T. Maekawa, F. Wolter, and N. Patrikalakis. Umbilics and lines of curvature for shape interrogation. Computer Aided Geometric Design, 13:133–161, 1996. Ridges and umbilics of polynomial parametric surfaces [PAT00] [Por71] [Por01] [Rou99] [RZ03] [SW05] [Thi96] 25 X. Pennec, N. Ayache, and J.-P. Thirion. Landmark-based registration using features identified through differential geometry. In I. Bankman, editor, Handbook of Medical Imaging. Academic Press, 2000. I. Porteous. The normal singularities of a submanifold. J. Diff. Geom., 5, 1971. I. Porteous. Geometric Differentiation (2nd Edition). Cambridge University Press, 2001. F. Rouillier. Solving zero-dimensional systems through the rational univariate representation. Journal of Applicable Algebra in Engineering, Communication and Computing, 9(5):433–461, 1999. F. Rouillier and P. Zimmermann. Efficient isolation of polynomial real roots. Journal of Computational and Applied Mathematics, 162(1):33–50, 2003. Raimund Seidel and Nicola Wolpert. On the exact computation of the topology of real algebraic curves. In SCG ’05: Proceedings of the twenty-first annual symposium on Computational geometry, pages 107–115, New York, NY, USA, 2005. ACM Press. J.-P. Thirion. The extremal mesh and the understanding of 3d surfaces. International Journal of Computer Vision, 19(2):115–128, August 1996.