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.