Stefanie Hahmann
email: stefanie.hahmann@imag.fr

Alexander Belyaev
email: belyaev@mpi-sb.mpg.de

Laurent Busé

Gershon Elber
email: gershon@cs.technion.ac.il

Bernard Mourrain
email: bernard.mourrain@sophia.inria.fr

Christian Roessl
email: christian.roessl@sophia.inria.fr

Christian Rössl

Shape Interrogation

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.

Introduction

Shape interrogation is the process of extraction of information from a geometric model. Surface interrogation is of central importance in modern Computer Graphics and Computer Aided Design (CAD) systems. Wherever geometrical models are used, they often need to be analyzed with respect to different aspects like, for example, visual pleasantness, technical smoothness, geometric constraints or surface intrinsic properties. The various methods, which are presented in this survey can be used to detect surface imperfections, to analyze shapes or to visualize different forms. We not only restrict the shapes to be investigated to free-form surfaces, but include polygonal meshes as well as algebraic surfaces. Artefacts of subdivision surfaces are subject of Chapter 4 of this book [START_REF] Catalano | Subdivision surfaces[END_REF]. Particular attention is paid to stable numerical and symbolic solving of algebraic systems of equations, a problem that arises in most shape interrogation methods.

In Section 2, fundamental notions of differential geometry are briefly recalled. Interrogation methods for polygonal meshes are discussed in Section 3. First and second order shape interrogation and visualization techniques are discussed in Sections 4, 5, focusing mainly on free-form curves and surfaces. The computation and visualization of characteristic curves on surfaces is subject of Section 6. Section 7 discusses the use of robust symbolic computation methods for shape interrogation. Interrogation of algebraic curves and surfaces is finally discussed in Section 8, in particular the transversal problem of solving of algebraic systems of equations is described.

Differential Geometry of curves and surfaces

Fundamental notions of differential geometry of curves and surfaces that are needed in the following of the paper will briefly be reviewed in this section. For a complete bibliography on differential geometry the reader is referred to standard literature [START_REF] Do | Differential Geometry of curves and surfaces[END_REF]102,[START_REF] Hilbert | Geometry and the imagination[END_REF]176].

Curves

A parametric curve is a mapping x from I = [a, b] ⊂ IR into IR n of class C r (r ≥ 1). x is called regular, if dx dt (t) = 0 for all t ∈ I. If L is the length of x([a, b]), there exists a unique parameter transformation s from I into [0, L] such that for all t 0 , t 1 ∈ [0, L] the length of the arc x([t 0 , t 1 ]) is equal to s(t 1 ) -s(t 0 ). For all t ∈ [a, b] s(t) = t a dx dt dt. s is called the arc length parameterization. It is a geometric invariant of a curve and is therefore also called natural parameterization. Let x : [0, L] → IR 3 , s → x(s) be a regular and naturally parameterized curve of class C 3 , such that x ′′ (s) = 0 for all s ∈]0, L[, then

• v 1 (s) := x ′ (s) is called tangent vector of x in s.

• v 2 (s) := x ′′
x ′′ is called unit normal vector of x in s. • v 3 (s) := v 1 (s) × v 2 (s) is called binormal vector of x in s, where × denotes the vector product (cross product) in IR 3 . {v 1 (s), v 2 (s), v 3 (s)} form an orthonormal basis of IR 3 called the Frenet frame of x in s. The following holds: v 1 , v 2 , v 3 are mappings of class C 1 , and

v ′ 1 = κ 1 v 2 v ′ 2 = -κ 1 v 1 + κ 2 v 3 v ′ 3 = -κ 2 v 2 where κ(s) = x ′′ , τ (s) = |x ′ , x ′′ , x ′′′ | x ′′
are mappings of class C 1 and C 0 respectively. |•, •, •| denotes the determinant of the matrix formed by the three vector arguments. κ and τ are called curvature and torsion of the curve x. The curvature measures the deviation of a curve from a straight line, and the torsion measures the deviation of a curve from being planar.

Surfaces

A parametric surface is a mapping X from Ω ⊂ IR 2 into IR 3 of class C r (r ≥ 1). X is called regular if for all u = (u, v) ∈ Ω, dX u is an invertible linear mapping. The two partial derivatives of X in u are denoted by X u (u) and X v (u). The affine subspace T u X := {X(u) + λX u (u) + µX v (u) | (λ, µ) ∈ IR 2 } is called tangent plane to X in u.

The unit normal vector field N is given by

N := X u × X v X u × X v .
The moving frame {X u , X v , N } is the Gauss frame. The Gauss frame is in general not an orthogonal frame.

The bilinear form on T u X induced by the inner product of IR 3 is called the first fundamental form of the surface. The matrix representation of the first fundamental form I u with respect to the basis {X u , X v } of T u X is given by G = (g ij ) with i, j = 1, 2:

g 11 g 12 g 21 g 22 = X u , X u X u , X v X v , X u X v , X v
where <, > denotes the scalar product. The first fundamental form I u is symmetric, positive definite and geometrically invariant. The first fundamental form allows measurements on the surface (length of curves, angles of tangent vectors, areas of regions) without referring back to the space IR 3 , in which the surface lies. The linear mapping L u L u : T u X → T u X x → dN u • dX -1 u (x) is called the Weingarten map.

The bilinear symmetric form II u defined on T u X by

II u (x, y) = L u (x), y
is called the second fundamental form of the surface X.

Its matrix in the basis {X u , X u } of T u X is denoted H = (h ij ) with i, j = 1, 2:

h 11 h 12 h 21 h 22 = N, X uu N, X uv N, X vu N, X vv .

The matrix HG -1 of the Weingarten map L u is symmetric and real and therefore it has two real eigenvalues κ 1 , κ 2 with corresponding orthogonal eigenvectors. κ 1 , κ 2 are called principle curvatures of the surface X, also labeled as κ max , κ min . The product of the principle curvatures

K = κ 1 • κ 2 = det(L u ) = det(H) det(G)
is called the Gaussian curvature and its mean M = 1 2 (κ 1 + κ 2 ) = trace(L u ) is called the mean curvature.

Another approach for the principle curvatures is the following: Let A := ∆u • X u + ∆v • X v be a tangent vector with A = 1. If we intersect the surface with the plane given by N and A, we get an intersection curve y with the following properties: ẏ(s) = A and e 2 = ±N where e 2 is the principal normal vector of the space curve y. The implicit function theorem implies the existence of this normal section curve. To calculate the extreme values of the curvature of a normal section curve (the normal section curvature) we can use the method of Lagrange multipliers because we are looking for the extreme values of the normal section curvature κ N with the condition ẏ(s) = 1.

As a result of these considerations we obtain the following. Unless the normal section curvature is the same for all directions there are two perpendicular directions A 1 and A 2 in which κ N attains its absolute maximum and its absolute minimum values. These directions are the principal directions with the corresponding normal section curvatures κ 1 and κ 2 .

For A = A 1 cos ϕ + A 2 sin ϕ we get Euler's formula:

κ N = κ 1 cos 2 ϕ + κ 2 sin 2 ϕ,
If the principal directions are taken as coordinate axes, Euler's formula implies the so-called Dupin indicatrix:

κ 1 (u) 2 + κ 2 (u) 2 = ±1. (1)
We use the Dupin indicatrices as a tool to visualize curvature situations on surfaces. The Dupin indicatrices at elliptic points (K > 0) are ellipses, at hyperbolic points (K < 0) pairs of hyperbolas, and at parabolic points (K = 0) pairs of parallel lines. Flat points (κ 1 = κ 2 = 0) are degenerated parabolic cases. Points with κ 1 = κ 2 are called umbilical points.

Interrogation of discrete shapes

Polygonal meshes constitute the primary tool for 3D surface representation and are frequently used in a wide range of scientific applications, including computer graphics, visualization, and numerical simulations. Two fundamental questions of surface approximation by polygonal meshes concern approximation quality (accuracy) [START_REF] Frey | Surface mesh quality evaluation[END_REF] and the relation between the accuracy and size of the approximation [START_REF] Friedel | Variational normal meshes[END_REF]. Recently both of these questions were also addressed in [START_REF] Cohen-Steiner | Variational shape approximation[END_REF] where a variational approach for surface approximation by polygonal meshes was developed. Shape approximation with polygonal meshes is discussed in more detail in Chapter 2 of this book [START_REF] Alliez | Recent advances in remeshing surfaces[END_REF]. Accurate estimation of geometric properties of a surface from its discrete approximation is important for many applications. Nevertheless there is no consensus on how to achieve accurate estimations of simple surface attributes such as the normal vector and curvatures [122]. An accurate polygonal approximation of surface geometry in a least-squares sense [START_REF] Frey | Surface mesh quality evaluation[END_REF][START_REF] Cohen-Steiner | Variational shape approximation[END_REF] does not guarantee accurate approximations of surface normals and curvatures by their discrete counterparts [121,132,119,[START_REF] Borrelli | On the angular defect of triangulations and the pointwise approximation of curvatures[END_REF]. Thus, deriving accurate, consistent, and numerically robust estimates for the surface normal vector and curvature tensor remains an area of active and creative research today.

Surface Normal Estimation

Given a smooth surface approximated by a dense triangle mesh, an accurate and robust estimation of vertex normals is important for a number of tasks including smooth shading [START_REF] Gouraud | Continuous shading of curved surfaces[END_REF]156], curvature estimation (see, e.g., [180]), and feature extraction (see, e.g., [START_REF] Hoschek | A geometric concept of reverse engineering of shape: Approximation and feature lines[END_REF]).

Usually the normal vector at a vertex of a triangle mesh is estimated as the normalized weighted sum of normals of the incident facets (triangles). A survey of various methods to estimate the normal vector can be found in [174]. Uniform (equal) weights are justified in [START_REF] Glassner | Computing surface normals for 3D models[END_REF] via finite difference approximations. In [180] the weights are chosen to be equal to the areas of the incident triangles. Weighting by the inverse areas was considered in [174,[START_REF] Hoschek | A geometric concept of reverse engineering of shape: Approximation and feature lines[END_REF], and weights equal to the facet angles at the vertex are proposed in [185]. A weighting scheme assuming that the mesh locally approximates a sphere was developed in [120]. The vertex normal vector can be also obtained from the mean curvature vector and, therefore, mean curvature vector estimates proposed in [START_REF] Desbrun | Implicit fairing of irregular meshes using diffusion and curvature flow[END_REF]122] lead to approximations of the vertex normals. A standard approach for testing and comparing various methods to estimate surface normals and curvatures consists of tessellating known (analytical) surfaces and comparing the estimates from the resulting mesh and from the original surface [START_REF] Hamann | Curvature approximation for triangulated surfaces[END_REF]180,104,[START_REF] Csákány | Computation of local differential parameters on irregular meshes[END_REF]122]. An interesting statistical approach was recently proposed in [125,126]. First steps towards a rigorous mathematical analysis and comprehensive comparison of various weighting schemes are made in [106].

Curvature Tensor Estimation

Estimates of the curvature tensor on polygonal meshes are applied in a variety of applications ranging from the detection of surface defects to the detection of features. Many techniques have been proposed (see, e.g., [153] for a recent survey), in this section we provide an overview of different approaches.

In order to estimate the curvature tensor at a vertex a certain neighborhood of this vertex is considered, typically its 1-ring. A common approach is to first discretize the normal curvature along edges. Given is an edge (i, j), vertex positions X i , X j , and the normal N i , then

κ ij = 2 (X j -X i ), N i X j -X i 2 (2)
provides an approximation of the normal curvature at X i in the tangent direction which results from projecting X i and X j into the tangent plane defined by N i . This expression can be interpreted geometrically as fitting the osculating circle interpolating X i and X j with normal N i at X i (cf. [130]). Alternatively, the equation can be derived from discretizing the curvature of a smooth planar curve (cf. [180]). With estimates κ ij of the normal curvature for all edges incident to vertex i, Euler's formula can be applied to relate the κ ij to the unknown principal curvatures (and principal directions). Then approximates to the principal curvatures can be obtained either directly as functions of the eigenvalues of a symmetric matrix ([180, 147]) or from solving a least-squares problem ([130, 122]). Alternatively, the trapezoid rule is applied in [188] to get a discrete approximation of the mean curvature M expressed as the integral over the normal curvatures κ N , the Gaussian curvature K is obtained from a similar integral over κ 2 N , then M and K define the principal curvatures. Exact quadrature formulas for curvature estimation are provided in [107].

Another class of techniques for curvature tensor estimation locally fits a smooth parametric surface patch and then derives the differential quantities from that. This leaves the choice for the surface -typically polynomials of low degree -the geometric quantities to interpolate or approximate -e.g., the vertex positions in a 1-ring neighborhood -and a projection operator to obtain a parameterization -in general the projection into the tangent plane. A straightforward choice is to consider the quadratic height surface

z(x, y) = 1 2 a 20 x 2 + a 11 xy + 1 2 a 02 y 2 ,
for a local coordinate system spanned by the normal N i (in z-direction) and two orthogonal tangent vectors (in x-and y-direction) and with origin X i = 0 [START_REF] Goldfeather | A novel cubic-order algorithm for approximating principal directions vectors[END_REF]. Then the parameters a 20 , a 11 , and a 02 obtained as a leastsquares solution are the elements of the symmetric matrix defining the Weingarten map. This can be interpreted as estimating the normal curvature from parabolas rather than circles (as with ( 2)) and then solving a least-squares system like in [122].

In [189] a quadratic Taylor polynomial of different form is applied, namely

X(u, v) = X u u + X v v + 1 2 u 2 X uu + X uv uv + 1 2 v 2 X vv .
The coefficients of the local least-squares approximating polynomial are the first and second order partials and hence define the fundamental forms. For robustness reasons, an exponential map is used as projection operator rather than a simple projection to the tangent plane. The use of a cubic approximation scheme which takes into account vertex normals in the 1-ring is proposed in [START_REF] Goldfeather | A novel cubic-order algorithm for approximating principal directions vectors[END_REF]. As the normals themselves are local estimates, this effectively enlarges the neighborhood. Again, a least-squares problem is solved to find the coefficients of a cubic height surface, where the Weingarten matrix is obtained entirely from the quadratic terms in the same way as before.

In general, least-squares methods may suffer from degenerate cases -even for reasonable geometric configurations -which lead to ill-conditioned system matrices. In [189] the polynomial basis is successively reduced in such cases. An alternative is to provide more samples e.g. from linear interpolation. In [START_REF] Cazals | Estimating differential quantities using polynomial fitting of osculating jets[END_REF] the patch fitting approach is discussed from an approximation theory point of view including robustness and numerical issues. For high-quality and consistent estimation of curvatures and their derivatives, [145] applies a (rather expensive) global fitting of an implicit surface to the surface mesh.

In contrast to the previously mentioned techniques, tensor averaging methods estimate the curvature tensor as an average over a certain region of a polyhedral mesh. In [START_REF] Cohen-Steiner | Restricted delaunay triangulations and normal cycle[END_REF] the curvature tensor is derived building upon the theory of normal cycles. This work includes a proof of convergence under certain sampling conditions based on geometric measure theory. The curvature tensor is defined at each point along an edge, and all contributions are integrated over a small region, see also [START_REF] Alliez | Anisotropic polygonal remeshing[END_REF]. A similar discrete curvature measure is applied in [START_REF] Hildebrandt | Anisotropic filtering of non-linear surface features[END_REF].

Alternative approaches locally consider a triangle with given vertex normals. In [167], the directional derivatives of the normal are expressed as finite differences for every edge of a triangle. The resulting system of six equations is set up from the vertex positions (in parameter space) and normals and then solved for the three unknowns of the Weingarten matrix in leastsquares sense. The tensors which are obtained per triangle are transformed to a common coordinate system to get a per-vertex average over the 1-ring. The algorithm can be applied with only slight modifications to compute curvature derivatives from the prior result.

In [181] the curvature tensor is estimated as smooth function (rather than a constant value) per triangle. This technique is inspired by Phong shading [156], where the vertex normals are linearly interpolated over the triangle. These interpolated normals are used to define the first and second order partials of the unit normal. This yields a piecewise smooth function defining the curvature tensor and elegant expressions for the Gaussian and mean curvature. Although this function is in general not continuous over edges of the triangulation, the approximation error is comparable to other approaches. For the estimation at vertices, the error is reduced by taking averages from all incident triangles.

Applications to Discrete Shape Analysis

The techniques reviewed in the previous section enable the estimation of curvature on discrete shapes: curvature estimates such as principal curvatures, Gaussian curvature and mean curvature are available at every vertex. These values can then be linearly interpolated in triangles. This is illustrated in define a vector field on the surface. Figure 2(c) and (d) shows lines of curvature obtained from stream line integration. In addition to these examples, many surface interrogation methods which were initially developed for smooth surfaces can be adapted easily to work in the discrete setting. This applies to first order analysis (Section 4) using estimates of the surface normal: reflection lines can be simulated by environment mapping techniques, highlight lines and isophotes can be emulated similarly. With curvature estimates being available, second order analysis (Section 5) can be applied. For the computation of discrete characteristic lines (Section 6), curvature derivatives are approximated by appropriate differences.

The following sections discuss shape analysis of smooth surfaces. Interrogation of discrete shapes follows the general ideas closely and applies estimates of surface normals and curvature.

First-Order Shape Analysis

First-order surface interrogation methods make generally use of the surface normal vector by simulation of particular light reflecting behavior of the surface. The light reflection methods all simulate the special reflection behavior of light sources or light lines on the surface. Due to the intuitive understanding that everybody has when he observes light reflections, these methods are very effective in detecting surface irregularities. They are therefore very well suitable for testing the fairness of surfaces. Because the surface normals are involved in the computation of these lines, they also can be used to visualize first order discontinuities, like tangent discontinuities.

Reflection lines

The reflection line method determines unwanted dents by emphasizing irregularities in the reflection line pattern of parallel light lines. Let X(u, v) be a representation of the surface to investigate, and let N (u, v) be the unit normal vector of the surface. Furthermore a light line L is given in parameter form:

L(t) = L 0 + t • s
where L 0 is a point on L, s is a vector defining the direction of L, t ∈ IR. The reflection line is the projection of the line L on the surface X, which can be seen from the fixed eye point A, if the light line L is reflected on the surface, see Figure 3(a). From geometric dependencies the following reflection condition is derived:

b + λa = 2 N (u, v) • b N (u, v) with λ := b a , (3)
where a = P -A, b = L -P . Equation (3) has to be solved for the unknown parameters u and v of the reflection point P . These three non-linear equations can be reduced to two equations by eliminating λ; they can then be solved by numerical methods, but the existence and uniqueness of solutions has to be ensured by an appropriate choice of the eye point A [START_REF] Kaufmann | Smoothing surfaces using reflection lines for families of splines[END_REF][START_REF] Klass | Correction of local irregularities using reflection lines[END_REF]. To analyze visually the surface one uses a set of parallel reflection lines with direction s, a fixed eye point A, and one steps along each curve of the set. Figure 4(a) shows a reflection line pattern on a part of a hair dryer and visualizes some surface irregularities.

Highlight lines

A highlight line is defined as the loci of all points on the surface where the distance between the surface normal and the light line is zero. The linear light source idealized by a straight line with an infinite extension

L(t) = L 0 + Bt
(L 0 is a point on L, B is a vector defining the direction of L, t ∈ IR), is positioned above the surface under consideration, see Figure 3(b). The highlight line method also detects surface irregularities and tangent discontinuities by visualizing special light reflections on the surface. In comparison with the reflection line method, the highlight lines are calculated independently from any observers view point. For a given surface point X(u, v) let N (u, v) be the unit normal vector. The surface point X(u, v) belongs to the highlight line if both lines, L(t) and the extended surface normal

E(s) = X(u, v) + s • N (u, v) , s ∈ IR intersect, i.e. if the perpendicular distance d = [B × N ] • [L 0 -X] [B × N ]
between these lines is zero, see Figure 3(b). This method can be extended to highlight bands, lines where d ≤ r (r fixed) is verified. For details on the algorithms to compute highlight lines see [START_REF] Beier | Highlight-line algorithm for realtime surface quality assessment[END_REF].

Isophotes

Isophotes are lines of equal light intensity. If X(u, v) is a parameterization of the surface and L the direction of a parallel lighting, then the isophote condition is given by:

N (u, v) • L = c ,
where c ∈ IR is fixed, see Figure 3(c). Note that silhouettes are special isophotes (c = 0) with respect to the light source. Isoclines are lines of equal normal inclination with respect to some direction V . If X(u, v) is a parameterization of the surface, then the isocline condition is given by:

N (u, v) • V = c
where N (u, v) is the unit normal field of X and c ∈ IR is fixed. In other words, isophotes are isoclines with respect to the light source direction. Similar to reflection lines and highlight lines, the isophotes provide a powerful tool to visualize small surface irregularities, which can not be seen with a simple wireframe or a shaded surface image. In Figure 4(b) we use 20 different values for c in order to get an isophote pattern on a NURBS test surface. Now, as stated out in the introduction of this section, the light reflection methods can be used to visualize first and second order discontinuities, because the surface normal vector is always involved in the line definitions. In fact, if the surface is C r -continuous, then the isophotes are C r-1 -continuous curves (see [157] for more details). A curvature discontinuity can be recognized, where the isophotes possess tangent discontinuities (breaks). One should nevertheless be careful by using isophotes for this purpose, because sometimes the break points of the isophotes at curvature discontinuities may not be clearly recognized, because of an ill-conditioned light direction. This special case occurs if the orthogonal projection of the light direction L in the tangent plane at a boundary point X(u, v) is parallel to the tangent of the isophote at this point.

Isophotes for curvature discontinuity:

There is another isophote method, which on one hand is an automatic method (independent of a special light direction), but which on the other hand only visualizes curvature discontinuities across the boundaries of a patch work. It makes use of the fact that along a common boundary curve y between two surface patches that join only with tangent plane continuity the Dupin indicatrices i 1 and i 2 on both sides are different. In general there are two conjugate diameters of the Dupin indicatrix. This relation degenerates at parabolic points, because the asymptotic direction (the direction in which the normal section curvature vanishes) is the conjugate to itself, but also conjugate to all other directions. At planar points, we have this degeneration for each (tangent) direction. Since both patches have a common boundary curve, and the tangent planes along that curve are unique, the Dupin indicatrices i 1 , i 2 have a common diameter, but differ in the other. We now consider an isophote c passing through P . The tangent t i of c at P with respect to X i is conjugate to the orthogonal projection f of the light ray onto the tangent plane (i = 1, 2), see Figure 5. In general the isophote c shows a tangent discontinuity at P if the Dupin indicatrices of X 1 and X 2 are not equal, but we have to avoid the situations f = ẏ = t and f = t ′ . More details can be found in [161].

Detection of inflections

Orthotomics and the polarity method are both interactive interrogation tools capable to detect only one particular type of surface "imperfection": the change of the sign in the Gaussian curvature. For example, surface with only convex iso-parameter lines are not necessarily convex, i.e. their Gaussian cur-vature is not required to be positive at all surface points. Such surface imperfections are difficult to detect visually in this case and therefore a curvature based surface analysis is needed like color maps or generalized focal surfaces, see Section 5. The following methods in contrast can visualize a change of sign in the Gaussian curvature without computing second order derivatives of the surface.

Orthotomics

In [START_REF] Hoschek | Detecting regions with undesirable curvature[END_REF] it has been shown that for a regular surface X(u, v) and for a point P that does not lie on the surface or on any tangential plane of the surface the k-orthotomic surface Y k (u, v) with respect to P defined by

Y k (u, v) = P + k (X(u, v) -P ) • N (u, v) N (u, v) ,
where N (u, v) is the unit normal vector of the surface has a singularity in (u 0 , w 0 ), if and only if the Gaussian curvature of X vanishes, or changes its sign at this point. To illustrate this method we consider a Bézier surface with completely convex parameter lines, see in Figure 6(left). But this surface is not convex: as shown in Figure 6(right), the orthotomic analysis emphasizes the change of sign of the Gaussian curvature in the corner region.

Polarity method

The polarity method is a further method able to detect unwanted changes in the sign of the Gaussian curvature without computing second order derivatives of the surface. It works for curves as well. It uses the polar image of a curve or surface, where the singularities (cusps, edge of regression) of this image indicate the existence of points with vanishing Gaussian curvature. The polar surface looks similar to the orthotomic surface, because the center of polarity is chosen to be equal to the projection point of the orthotomic analysis. For more information about the polarity method and on how removing the inflections see [START_REF] Hoschek | Smoothing of curves and surfaces[END_REF].

Geodesic paths on surfaces and meshes

Geodesic paths, or simply geodesics, on a surface are surface curves which connect two surface points with minimum path length. A thorough study of geodesics and their role of in surface interrogation requires much more attention than the present overview can provide. So below we give only a brief literature survey. Geodesics deliver rich information about surface geometry and, therefore, have various theoretical and practical applications. In particular, detecting geodesic paths on surfaces approximated by triangles meshes is a common operation for many graphics and modeling tasks such as mesh parameterization [103], mesh segmentation [START_REF] Katz | Hierarchical mesh decomposition using fuzzy clustering and cuts[END_REF], skinning [175], mesh watermarking [162], and mesh editing [100].

A rigorous mathematical treatment of geodesics can be found in [102,[START_REF] Eisenhart | An introduction to differential geometry[END_REF]. Some numerical aspects are presented in [START_REF] Farouki | Graphical methods for surface differential geometry[END_REF]. An algorithm for approximate computation of geodesic paths on smooth parametric surfaces has been explored in [155,154]. Various algorithms exit for computing geodesic paths and distances. The so-called MMP algorithm [124] computes an exact solution for the "single point, all distances" shortest path problem by partitioning each mesh edge into a set of intervals over which the exact distance can be computed. In [179] an accelerated implementation of this algorithm is presented. An algorithm to solve the "single source, single distance" geodesic problem is given in [START_REF] Kapoor | Efficient computation of geodesic shortest paths[END_REF]. See also [123] for a broad survey of algorithms for computing shortest paths on graphs.

Second-order shape analysis

Surface curvature is of central importance for surface design. Often the result is required to be mathematically smooth (continuous in the 2nd derivative) and aesthetically pleasing, i.e. have smooth flowing highlights and shadows. To obtain an aesthetically pleasing shape, the designer works with the curvature. A color map (see Section 5.5) can be used to visualize curvature (Gaussian, principal curvatures) over the surface. The problem is the good choice of the color scale, which depends on the curvature function and therefore on the underlying surface.

The surface interrogation methods presented in this section are therefore curvature analysis tools which are able to detect all surface imperfections related to curvature, like bumps, curvature discontinuity, convexity, and so on.

Local shape analysis with Gaussian curvature

Let us look at a smooth surface in a neighborhood of one of its point.

The simplest classification of local surface shapes is given by the the sign of the Gaussian curvature The Gaussian curvature of a surface can be expressed through the coefficients of the first fundamental form. Thus we arrive at the following famous result called Gauss's Theorema Egregium: the Gaussian curvature of a surface is a bending invariant. Now let us consider a simple geometrical interpretation of the Gaussian curvature, by means of which Gauss originally introduced it.

K = κ 1 • κ 2 . K > 0. The normal curvatures κ N (ϕ)
Consider a two-sided surface in three-dimensional space. Let us transport the positive unit normal vector from each point of the surface to the origin. The ends of these vectors lie on the unit sphere. We obtain the mapping of the surface into the unit sphere, see Figure 9. It is called the Gauss map.

The Gauss mapping takes areas on surfaces to areas on the unit sphere. Consider the unit surface normals at the surface points within the area ∆S on the surface. Let us denote the area on the unit sphere (solid angle) corresponding to ∆S by ∆A. It turns out that the Gaussian curvature at the point is the limit of the ratio of these areas:

K = lim ∆S→0 ∆A ∆S .
This remarkable formula resembles the definition of the curvature of the plane curves: κ = dϕ/ds. The Gauss map can be used for detecting spherical, cylindrical, and conical regions on a surface [START_REF] Besl | Invariant surface characteristics for 3D object recognition in range images[END_REF].

Gaussian Sphere Surface Gauss map

S A A S K = lim

Focal Surface and Corresponding Surface Features

For a smooth surface X = X(u, v) its focal surface is given by

X F (u, v) = X(u, v) + N (u, v) κ(u, v) , κ = κ 1 , κ 2 ,
where N (u, v) is the oriented normal. The focal surface is formed by the principal centers of curvature and consists of two sheets corresponding to the maximal and minimal principal curvatures κ 1 and κ 2 . One can show that the focal surface is the envelope of the surface normals. In geometrical optics [START_REF] Herman | A Treatise on Geometrical Optics[END_REF], a caustic generated by a family of rays is defined as the envelope of the family. Thus the focal surface is the caustic of the family of surface normals.

The focal surface can be also defined as a surface swept by the singularities of the offset surfaces

O d (u, v) = X(u, v) + d N (u, v).
The focal surface is the 3D analogue of the evolute of a planar curve and has singularities. The singularities of the focal surface consist of space curves called focal ribs.

Ridges, the surface curves corresponding to the focal ribs are natural generalization of the curve vertices for surfaces. The ridges can be defined as sets of surface points where the principal curvatures have extremes along their associated principal directions and points where the principal curvatures are equal (umbilics). A thorough study of the ridges and their properties is conducted by Porteous [158]. See also [START_REF] Hallinan | Twoand Tree-Dimensional Patterns of the Face[END_REF] where a detail classification of the ridges is presented. Below we briefly discuss the ridges from a singularity theory point of view.

Near a point on a focal rib the focal surface can be locally represented in the parametric form (c 1 t 3 , c 2 t 2 , s), where c 1 = 0 and c 2 = 0, in well chosen coordinates (s, t). The focal ribs themselves have singularities at points corresponding to the umbilics and those ridge points where one of the principal curvatures has an inflection along its corresponding curvature line. Generic (typical) singularities of the focal surface are shown in Figure 10. The umbilics and ridge points can be also characterized as surface points where the osculating spheres (spheres of curvature) have high-order contacts with the surface. Therefore the umbilics and ridges are invariant under inversion of the surface with respect to any sphere.

The focal surface points can be also described in terms of degenerate singular points of distance functions. Given a surface and a point in 3D, let us consider the distance function from the point and restrict the function onto the surface. This gives a three-dimensional family of distance functions defined on the surface and parameterized by points in 3D. Now the focal surface is generated by those point-parameters for which the distance function has degenerate critical points. A typical degenerate critical point has on of the following two forms ±s 2 + t 3 in proper coordinates s and t. If the pointparameter is a typical point on a focal rib, the distance function has a critical point in one of the following four forms: ±s 2 ± t 4 . More degenerate critical points occur when the point-parameter is located either at a swallowtail singularity of the focal surface or at an umbilical points. It is interesting that the cut locus of the surface [190] (skeleton or medial axis of a figure bounded by the surface) consists of those point-parameters which define the distance functions with two equal global minima. Thus, as illustrated in Figure 11, the edges of the skeleton are located at focal ribs. The focal surface possesses many interesting properties. For example, for each line of curvature on a surface there is a corresponding line on the corresponding sheet of the focal surface. It can be shown that those raised lines of curvature are geodesics on the focal surface [159,131].

In [118,[START_REF] Ko | Shape intrinsic fingerprints for free-form object matching[END_REF] umbilics are used for shape interrogation and shape matching purposes. Statistics of various types of umbilics on random surfaces computed and analyzed in [START_REF] Berry | Umbilic points on gaussian random surfaces[END_REF] may have have many potential applications for for inspecting and interrogating surface properties.

Hedgehog diagrams and curvature plots

The hedgehog diagrams and curvature plots are well known interrogation tools for planar curves [START_REF] Beier | The porcupine technique: principles, applications, and algorithms[END_REF][START_REF] Farin | Curvature and the fairness of curves and surfaces[END_REF]. A hedgehog diagram for planar curves visualizes the curve normals proportional to the curvature value at some curve points. A new curve is obtained by Xhedgehog (t) = X(t) + κN (t) thus visualizing curvature distribution and discontinuity. The inspection of surfaces with these methods can be done by applying them to planar curves on the surface (intersections of the surface with planes). [START_REF] Kjellander | Smoothing of bicubic parametric surfaces[END_REF] shows an example of application. Hedgehog diagrams for entire surfaces are nevertheless difficult to interpret and are therefore not to be recommended.

Generalized focal surfaces

Although the idea of generalized focal surfaces is quite similar to hedgehog diagrams, their application area is much larger. Instead of drawing surface normals proportional to a function value, only the point on the surface normal proportional to the function is drawn. The loci of all these points is the generalized focal surface. This method was introduced by [START_REF] Hagen | Generalized focal surfaces : A new method for surface interrogation[END_REF], and is based on the concept of focal surfaces which are known from line geometry, introduced in Section 5.2. The generalization of this classical concept leads to the generalized focal surfaces:

F (u, v) = X(u, v) + s • f (κ 1 , κ 2 ) • N (u, v) , with s ∈ IR
where N is the unit normal vector of the surface X. f is a real valued function of the parameter values (u, v). The variable offset function f can be any arbitrary scalar function, but in the context of surface interrogation it is quite natural to take f as a function depending on the principal curvatures κ 1 , κ 2 of X, f.ex.

f = κ 1 κ 2 Gaussian curvature, f = 1 2 (κ 1 + κ 2 ) mean curvature, f = (κ 2 1 + κ 2 2 ) energy functional, f = |κ 1 | + |κ 2 | absolute curvature, f = κ i principal curvatures, f = 1 κi
focal points, f = const offset surfaces. This not only enables to visualize a particular curvature behavior, but it can interrogate and visualize surfaces with respect to various criteria: A convexity test can be performed using the Gaussian curvature offset f = κ 1 • κ 2 = K. A surface is locally convex at X(u, v), if the Gaussian curvature is positive at this point. Often a surface is called non-convex, if there is a change in the sign of the Gaussian curvature. the two surfaces X(u, v) and F (u, v) intersect at the parabolic points, see Figure 12(a). The generalized focal surface therefore pin points directly on the area where the sign of K changes in contrast to orthotomics (Section 4) which are also used to test the convexity. Flat points which are special umbilic points with κ 1 = κ 2 = 0 can be detected using

f = |κ 1 | + |κ 2 | as well as f = κ 2 1 +κ 2 2 .
Flat points are undesired surface points because they make the surface bumpy. Curvature discontinuity can be visualized through gaps in the surface F with f = κ 2 1 + κ 2 2 since it is a second order surface analysis tool, see Figure 12(b). Visualizing surface irregularities: Surfaces are aesthetically pleasing if they have "nice" light reflections. Thus similar to reflection lines the generalized focal surfaces are also a tool for visualizing such surface imperfections because they are very sensitive to small irregularities in the shape. In Figure 12(b) part of a hair dryer is shown. It consists of biquintic C 1 -continuous patches. The iso-parametric lines do not reflect the bump in the surface, which is however emphasized by the focal analysis. Another aspect of surface analysis is the visualization of technical aspects. A surface which should be treated by a spherical cutter is not allowed to have a curvature radius smaller than the radius of the cutter R cutter . The generalized focal surfaces are able to detect such undesired regions by intersection with the surface X. The offset function to choose in this special case, is f = 1 Rcutter -κ max . Figure 12(c) shows such a surface which is not allowed to be cut. Generalized focal surfaces not only visualize surface imperfections, they also give the user a 3D impression of the relative amount of the offset function over the surface, what color maps can't do.

Color mappings

Color is used to emphasize features on the surface. Texturing can emphasize the spatial perception of an 2D image of the surface. A color-coded map is an application, which associates to a scalar function value a specific color. The color scale presents an even gradation of color corresponding to the range of function values. Colors are principally used to visualize either continuously or discontinuously any scalar function over a surface [START_REF] Dill | An application of color graphics to the display of surface curvature[END_REF][START_REF] Beck | Surface analysis methods[END_REF][START_REF] Barnhill | Twists, curvature and surface interrogation[END_REF][START_REF] Forrest | On the rendering of surfaces[END_REF], like pressure, temperature, or curvature, see Figure 13. Colors are used as a fourth dimension and show the user immediately and quantitatively how the function varies over the surface. An even gradation of the linear or cyclic color coding is important to visualize the rapid curvature variation by the presence of color fringes. Beck et al. [START_REF] Beck | Surface analysis methods[END_REF] propose to use the HSI (hue, saturation, intensity) model and to perform transformations between this space and the three primary colors RGB. See [START_REF] Foley | Computer Graphics. Principles and Practice[END_REF] for more details on color spaces and transformations. An example of discrete color-coding of the interval [0,1] is the following one:

Interval

Red Green Blue Color 0.0 -0. The main difficult of this simple interrogation method is the choice of a convenient color scale, which obviously depends on the function values to be visualized.

Pseudo texture

The use of colors for displaying a surface helps to emphasize the 3D understanding of an 2D image by simulating shadows, perspective and depth of the object. An artificial texturing is an aid for visualizing rendered surfaces.

Isoparametric lines are commonly used, but they are in some situations ambiguous. Schweitzer [170] projects equally spaced dots of equal size over the surface in order to increase the visual perception of the form.

Characteristic lines

Drawing lines on surfaces is a powerful and widely used tool for analysis and visualization of surface features. The techniques of isolines, lines of curvature, geodesic paths and ridges are presented. Numerous graphical examples are illustrated in [159,[START_REF] Farouki | Graphical methods for surface differential geometry[END_REF]. In the last three cases a set of lines on the surface can be created, and should be interpreted with the knowledge of differential geometry. They are the most sophisticated tools from the mathematician's point of view. The user should interpret the lines of curvature or the geodesic paths.

Isolines

Isolines are lines of a constant characteristic value on the surface. They provide an interrogation tool with a wide variety of applications. They help analyzing surface characteristics, and they are used to visualize the distribution of scalar quantities over the surface. The visualization of a certain number of isolines, with respect to an even distribution of the characteristic values allows to study the behavior of these values.

Contour lines are planar lines on the surface which are all parallel to a fixed reference plane. Closed contour lines indicate maxima and minima of the surface with respect to the direction given by the plane's normal vector [START_REF] Hartwig | Isolinien und schnitte in coonschen flächen[END_REF][START_REF] Beck | Surface analysis methods[END_REF] Parabolic lines are isolines of zero Gaussian curvature on the surface. They are of particular interest for intrinsic surface interrogation, since they divide the surface into elliptic and hyperbolic regions and they reflect therefore the local curvature behavior of a surface. Parabolic lines are special Gaussian curvature lines, see Figure 14. In [START_REF] Hilbert | Geometry and the imagination[END_REF] a more complex example with the statue Apollo Belvedere is drawn.

Lines of curvature, umbilics

Lines of curvature are curves whose tangent directions coincide with those of the principal directions, which are orthogonal. They form therefore an orthogonal net on the surface.

The net of lines of curvature becomes singular at an umbilical point where κ 1 and κ 2 are identical and the principal directions are indeterminate. Some numerical integration method is used to calculate the lines of curvature. But the integration process becomes unstable near an umbilic. Unfortunately umbilics appear frequently on free-form surfaces. A recent work about umbilics [117], destined for use in CAGD (Computer Aided Geometric Design), presents a procedure to compute the lines of curvature near an umbilic. And in [116] a computational method to locate all isolated umbilics on parametric polynomial surfaces is described. The discrete field of principle curvature directions computed on a surface mesh has been used for remeshing [START_REF] Alliez | Anisotropic polygonal remeshing[END_REF]. More details about umbilics and lines of curvature figures are found in classical differential geometry literature [START_REF] Darboux | Leçons sur la théorie générale des surfaces[END_REF], or in a more recent book [159].

Curvature Extrema for Shape Interrogation

Surface features invariant under rotations, translations, and scaling are important for studying shapes of 3D objects. The ridge curves discussed briefly in Section 5.2 are among the most important view-and scale-invariant features of a smooth surface.

The ridges are defined as the extremes of the principal curvatures along their corresponding curvature lines and constitute powerful surface descriptors. They have been intensively studied in connection with research on the accommodation of the eye lens [START_REF] Gullstrand | Zur Kenntnis der Kreispunkte[END_REF], structural geology [163] and geomorphology [109], human perception [START_REF] Hoffman | Parts of recognition[END_REF], image analysis [191,129,127,[START_REF] Eberly | Ridges in Image and Data Analysis[END_REF]110], quality control of free-form surfaces [START_REF] Hosaka | Modeling of Curves and Surfaces in CAD/CAM[END_REF], reverse engineering [START_REF] Hoschek | A geometric concept of reverse engineering of shape: Approximation and feature lines[END_REF], analysis and registration of anatomical structures [START_REF] Guéziec | Medical image registration using geometric hashing[END_REF][START_REF] Grenader | Computational anatomy: An emerging discipline[END_REF]151], face recognition [START_REF] Hallinan | Twoand Tree-Dimensional Patterns of the Face[END_REF], and non-photorealistic surface rendering [START_REF] Interrante | Enhancing transparent skin surfaces with ridge and valley lines[END_REF]114,[START_REF] Decarlo | Suggestive contours for conveying shape[END_REF]. (See also references therein.)

An explanation of why some ridges are good for sketching complex 3D shapes can be found in [191]: given a grey-scale image of an illuminated 3D object, under general illumination and reflection conditions, the zero-crossings of the second directional derivative of the image intensity along the direction of the image intensity gradient occur near the extremes of the principal curvature along their principal directions. Thus the projections of ridges onto the image plane are usually located near edges, the most salient image features. Some subsets of ridges play an important role in perceptual shape organization. Human perception experiments suggest the so-called minima rule [START_REF] Hoffman | Parts of recognition[END_REF] which sets region boundaries along lines divides shapes into parts at negative minima of the principal curvatures along their lines of curvature. The minima rule was employed in [146] for mesh segmentation purposes.

The ridges on a surface have interesting relations with the skeleton (medial axis) of a figure bounded by the surface and can be described via high-order contacts between the surface and its osculating spheres. See [158,101,192,159,[START_REF] Belyaev | Ridges, ravines, and singularities[END_REF], [START_REF] Hallinan | Twoand Tree-Dimensional Patterns of the Face[END_REF]Chapter 6], and recent reviews in [START_REF] Cazals | Smooth surfaces, umbilics, lines of curvatures, foliations, ridges and the medial axis: a concise overview[END_REF][START_REF] Cazals | Ridges and umbilics of a sampled smooth surface: a complete picture gearing toward topological coherence[END_REF] for rigorous mathematical treatments revealing beautiful properties of these curvature features. Surface landmarks associated with the ridges were considered in [101,131,160]. Bifurcations of the ridges on dynamic shapes were studied in [159, [START_REF] Bruce | Ridges, crests and sub-parabolic lines of evolving surfaces[END_REF][START_REF] Bruce | Families of surfaces: focal sets, ridges and umbilics[END_REF]160,112].

Recently the so-called crest lines, a subset of the ridges consisting of the extremes of the principal curvature maximal in absolute value along its corresponding curvature line, draw much attention because of their ability to represent surface creases [184,127,151,177,145,[START_REF] Hildebrandt | Smooth feature lines on surface meshes[END_REF]. See also references therein. One motivation for describing surface creases as the crest lines is based upon the following analogy with edges of grey-scale images [145]. Consider a surface and its Gauss map which associates with every point p of the surface the oriented normal vector n(p). The derivative ∇n(p) (Jacobian matrix) of the Gauss map measures the variation of the normal vector near p, i.e., how the surface bends near p. It is easy to see that the eigenvalues and eigenvectors of ∇n(p) are the principal curvatures and principal directions of the surface at p, respectively. Thus the maximal variation the surface normal is achieved in the principal direction of the principal curvature maximal in absolute value. So it is natural to define surface creases as loci of points where the positive (negative) variation of the surface normal in the direction of its maximal change attains a local maximum (minimum). Figure 16 shows the crest lines detected on various models represented by dense triangle meshes.

Practical detection of the ridges and their subsets is a difficult computational task since it involves estimating of high-order surface derivatives. Various techniques were proposed for detecting the ridge lines and their subsets on For shape interrogation purposes (shape quality control and analysis of aesthetic free-form surfaces), the ridges were used in [START_REF] Hosaka | Modeling of Curves and Surfaces in CAD/CAM[END_REF][START_REF] Higashi | Analysis of aesthetic free-form surfaces by surface edges[END_REF]. Moreton and Sequin [130] used the sum of the squared derivatives of the principal curvatures along their corresponding curvature lines as a measure of surface fairness.

Often, instead of the ridges and their subsets defined via extremes of the principal curvatures, simpler surface features are detected. In geometric modeling, there has been considerable effort to develop robust methods for detecting surface creases, curves on a surface where the surface bends sharply. Interesting methods for crease detection on dense triangle meshes and pointsampled surfaces were proposed in [START_REF] Hoschek | A geometric concept of reverse engineering of shape: Approximation and feature lines[END_REF]166,[START_REF] Hubeli | Multiresolution feature extraction from unstructured meshes[END_REF][START_REF] Gumhold | Feature extraction from point clouds[END_REF]178,148,150] Whereas the ridges were first studied one hundred years ago [START_REF] Gullstrand | Zur Kenntnis der Kreispunkte[END_REF] and have rich history [159], the so-called sub-parabolic lines, the loci of points where one of the principal curvatures has an extreme value when moving along the curvature line corresponding to another principal curvature. The subparabolic lines were introduced in [START_REF] Bruce | Folding maps and focal sets[END_REF] and studied in [159, 131, 160]. They possess many remarkable properties: the sub-parabolic lines correspond to the parabolic lines on the focal surface, hence the name, and consist of geodesic inflections of the lines of curvature [131]. The sub-parabolic lines can be also detected by examining the profiles of surfaces [131].

Special Surface Points

In this section, following [131] we consider special surface points which lie on the ridges and sub-parabolic lines. We adapt the color scheme proposed by Porteous [158,159]. Let us give the principal curvatures and corresponding principal directions, parabolic lines, and sheets of the focal surface a color (red or blue) in order to distinguish between them. The red (blue) sub-parabolic line consists of the extremes of the red (blue) principal curvature along the blue (red) curvature line. The following surface landmarks are useful for surface interrogation purposes:

• Umbilic points. See [118,149] for application of umbilics in surface matching and shape interrogation. • A ridge and sub-parabolic line of the same color cross. The principal curvature of the same color takes an extreme value there (maximum, minimum, or saddle). and the shape index

S = - 2 π arctan κ 1 + κ 2 κ 1 -κ 2 .
These measures are often more convenient for practical purposes then the standard curvature descriptors {κ 1 , κ 2 } and {M, K}, where K and M are the Gaussian and mean curvatures, respectively. In [142] it was suggested to use local maxima of the curvedness to define surface corner points.

Robust Symbolic based Shape Interrogation and Analysis

Interrogation of polynomial and rational surfaces could be made with the aid of symbolic processing. The advantage of the symbolic approach over sampling of properties, like curvature, at a discrete set of point stems from the ability to analyze the properties globally and provide global (error) bounds. Many properties of free-form geometry are differential and can be derived after executing a few basic operations over the polynomial or rational representation of the original interrogated curve C(t) or surface S(u, v), namely: differentiations, summations and products. We also assume the availability of a zero set finding tool, an operation that is equivalent to intersecting a polynomial or a rational function with a line in R 2 (a plane in R 3 ). As a simple example, consider the curvature field of a planar regular curve C(t) = (x(t), y(t)) that is equal to:

κ(t) = x ′ (t)y ′′ (t) -y ′ (t)x ′′ (t) (x ′ 2 (t) + y ′ 2 (t)) 2/3 .
κ(t) is not rational due to the fractional power in the denominator, in the normalization factor. Nonetheless, if one only seeks the inflection points of C(t), only the numerator of κ needs to be considered. Then, the solution of the constraint of

x ′ (t)y ′′ (t) -y ′ (t)x ′′ (t) = 0 (4) finds all the inflection points in the regular planar curve C(t), if any. In Equation ( 4), the problem of finding all the inflection points of a planar regular curve was reduced to that of finding a zero set. Differentiation and products were used to compute the inflection points' constraints. Differentiation of piecewise polynomials and rationals is well known [START_REF] Cohen | Geometric Modeling with Splines: An Introduction[END_REF][START_REF] Farin | Curves and Surfaces for Computer Aided Geometric Design[END_REF]. Similarly, the addition of two (piecewise) polynomials that share a function space (same order and knot sequence) is realized by simply adding the corresponding coefficients. Two polynomials could be elevated to the same function space via knot insertion and degree elevations; see [START_REF] Cohen | Geometric Modeling with Splines: An Introduction[END_REF][START_REF] Farin | Curves and Surfaces for Computer Aided Geometric Design[END_REF] for more details. Products are the last operator we seek, an operation also required because of the quotient rules over addition and differentiation of rationals. Products are more complex to compute (see [START_REF] Elber | Free form surface analysis using a hybrid of symbolic and numerical computation[END_REF][START_REF] Farin | Curves and Surfaces for Computer Aided Geometric Design[END_REF]) but, clearly, products of piecewise polynomials and/or rationals are piecewise polynomials and/or rationals as well.

In summary, the ability to form a closure and compute a differential property in the piecewise polynomial and/or rational domains, makes it far simpler and robust to analyze that property. While κ is not rational, its numerator is and so inflection points could be detected as a zero set of x ′ (t)y ′′ (t)-y ′ (t)x ′′ (t). For similar reasons, the unit normal N (t) of C(t) is not rational but both κ(t)N (t) and N (t)/κ(t) are rational. Hence, x-extreme points and y-extreme points on C(t) can be identified as κ(t)N (t), (0, 1) = 0, and κ(t)N (t), (1, 0) = 0, and the local maximum curvature locations in C(t) are detectable [START_REF] Elber | Symbolic and numeric computation in curve interrogation[END_REF] as the zeros of d κ(t)N (t), κ(t)N (t) dt , yet another rational function.

In [START_REF] Elber | Symbolic and numeric computation in curve interrogation[END_REF], points of extreme curvature, or alternatively, inflection points are detected using these schemes. In addition, a scheme to approximate an arclength reparametrizations for piecewise polynomial and/or rational curves is presented.

In the next section, Section 7.1, we will demonstrate the power of symbolic based interrogation in geometric design, for curvature analysis. In Section 7.2, silhouette curves, isoclines and isophotes curves, and reflection curves are all shown to be reducible to zero set finding. Then, in Section 7.3, we consider the problem of symbolic recognition of simple primitive surface shapes.

Curvature Analysis

Reexamining the second order differential analysis of parametric surfaces (recall Section 2), it turns out that given a rational surface S(u, v), the Gaussian curvature K is rational whereas the mean curvature M is not (while M 2 is). In [START_REF] Elber | Second-order surface analysis using hybrid symbolic and numeric operators[END_REF], a rational form of (the numerator of) K is symbolically computed and its zeros are used to robustly extract the parabolic lines of the surface. Figure 17 presents one example of computing the parabolic curves for a bicubic surface patch as the zeros of K. While M is not rational, one can compute M 2 as a rational form. Similarly, the form of κ 2 1 + κ 2 2 , where κ i , i = 1, 2, are the two principle curvatures, is rational and can capture regions that are highly curved. By subdividing the original surface into regions that prescribe different values of κ 2 1 + κ 2 2 , one can separate the surface into regions that could be NC-machined more efficiently with different sizes of ball-and flat-end cutters [START_REF] Elber | Freeform surface region optimization for three-and five-axis milling[END_REF]. Let K 0 = κ 2 1 + κ 2 2 at S 0 = S(u 0 , v 0 ). Then, the normal curvature at S 0 is bounded from above by √ K or an NC ball end cutter of radius 1/ √ K could be locally fitted to S 0 without (local) gouging. Figure 18 shows one such example where a surface is divided into regions of different values of extreme curvature, K = κ 2 1 + κ 2 2 . See also Equation (1).

Silhouette, Isoclines/Isophotes and Reflection lines

The extraction of silhouettes of a free-form surface could be easily reduced to a zero set finding problem. Looking at a rational surface S(u, v) from direction vector V , the silhouettes of S are characterized as the rational constraints of where N (u, v) = ∂S ∂u × ∂S ∂v . If the view is a perspective view through point P (the eye), the silhouettes could be derived as the rational form of

N (u, v), V = 0,
N (u, v), S(u, v) -P = 0.
Interestingly enough, highlight lines [START_REF] Beier | Highlight-line algorithm for realtime surface quality assessment[END_REF] (see Section 4.2), isoclines and isophotes (see Section 4.3) could be similarly reduced to a zero set finding, using symbolic manipulation. Let the unit view direction vector for which isoclines are sought be V . Then, positions on surface S(u, v) that present a normal with a constant inclination angle of α degrees could be characterized as

N (u, v) N (u, v) , V = cos(α),
which is not a rational but could be made into one by squaring both sides as,

N (u, v), V 2 -N (u, v) 2 cos 2 (α) = 0, (5)
at the cost of extraction both the + cos(α) and the -cos(α) isoclines, simultaneously. Figure 19 shows an example of subdividing a free form surface into regions of steep slopes (more than 45 degrees) and shallow slopes, using isoclines' analysis. Such a dichotomy might be desired, for example, in layered manufacturing processing where support is to be added to the geometry only below a certain slope. Reflection lines (see Section 4.1) can also be reduced to rational zero set constraints as follows. An incoming ray V that hits surface S(u, v) will be reflected in direction r(u, v),

r(u, v) = 2N (u, v) -V N (u, v), N (u, v) N (u, v), V . (6)
In practice, Equation ( 6) might be difficult to work with near silhouettes (where N (u, v), V vanish) and so, in [START_REF] Elber | Curve evaluation and interrogation on surfaces[END_REF], 2N (u, v) N (u, v), V -Fig. 19. Isoclines at 45 degrees from the vertical direction V . Left: the function whose zero set (see Equation ( 5)) prescribes the isoclines of the surface shown in the right figure is presented. Right: Isoclines also serve to split the surface into regions of slopes (normals) with more than 45 degrees (in thin lines) and regions of less than 45 degrees (in thick lines) with respect to the vertical direction V .

V N (u, v), N (u, v) was proposed as a better alternative. In [START_REF] Elber | Curve evaluation and interrogation on surfaces[END_REF], reflection ovals, or reflections of circular curves, were similarly reduced to zero set finding problems.

Surface Recognition

A fundamental question when analyzing free-form geometry is whether the given curve or surface is of a basic primitive nature. That is, a curve is tested if it is a circle, or a surface is tested if it is a cylinder, or alternatively, a surface of revolution. In [START_REF] Elber | Geometric shape recognition of freeform curves and surfaces[END_REF], these questions are answered using symbolic differential analysis. A rational curve is a circle if its rational squared curvature field, κ 2 (t) = κ(t)N (t), κ(t)N (t) , is constant. In other words, given a B-spline curve C(t), all its coefficient of the B-spline representation of κ 2 (t) should be the same and in fact equal to the square of the reciprocal of the radius of the circle. Alternatively, the evolute curve,

E(t) = C(t) + N (t)/κ(t),
which is also rational, should vanish (along with all its control points) at the circle's center locations.

A surface called the mean evolute surface,

E(u, v) = S(u, v) + N (u, v) 2M (u, v) ,
where M is the mean curvature (see Section 2.2) is also defined in [START_REF] Elber | Geometric shape recognition of freeform curves and surfaces[END_REF] and was shown to be rational for rational surface S(u, v). If S(u, v) is a circular cone, E(u, v) is reduced to a line, the cone's center axis. Figure 20 presents two such examples. In [START_REF] Elber | Geometric shape recognition of freeform curves and surfaces[END_REF], the connection is made between rational surfaces of revolution and rational pseudo-focal surfaces (see Section 5.4)

F u (u, v) = S(u, v) + N (u,v) κu(u,v)
, where κ u is the normal curvature of S(u, v) in the u isoparametric direction. If the u iso-parametric directions are the latitude lines of the surface of revolution, then F u reduces to the center axis line of the surface of revolution.

(a) (b) S(u, v) E(u, v) S(u, v) E(u, v)
Fig. 20. The mean evolute surface reduces to the center axis line of a circular cone or cylinder. In (a), a polynomial approximation of a cylinder surface S(u, v) with unconventional parameterization is presented along with its mean evolute E(u, v).

(b) presents a similar view of a portion of a polynomial approximation of a region of a circular cone along with its mean evolute.

For more information, see the recent book on shape interrogation in geometric design and manufacturing [149] that discusses many of the above topics as well as intersection problems, distance queries, curvature properties, and geodesics and offsets curves and surfaces.

Interrogation of algebraic curves and surfaces

In this section we will focus on particular geometric models: the algebraic curves and surfaces. We will show how to solve in this context some important shape interrogation problems as singularity detection, intersection problems and offset computation. It turns out that all these problems require at one point to solve an algebraic system of equations, this step being the crucial one. We thus articulate this section mainly on methods that can be applied on these algebraic systems.

Most of the curves and surfaces used in CAGD are given by parametric equations, as defined in Section 2. Planar rational curves in CAGD are typically defined as

x(t) = a(t) c(t) , y(t) = b(t) c(t)
where a(t), b(t) and c(t) are polynomials in the Bernstein basis for rational Bézier curves or in the B-spline basis for NURBS. Note that the algebraic methods most commonly use polynomials in the power basis and polynomials can be converted from Bernstein basis to power basis. Parametric rational surfaces in CAGD are defined by

x(u, v) = a(u, v) d(u, v) , y(u, v) = b(u, v) d(u, v) , z(u, v) = c(u, v) d(u, v)
where a(u, v), b(u, v), c(u, v) and d(u, v) are polynomials. Most of the shape interrogation problems for algebraic curves and surfaces can be translated in terms of a system of polynomial equations, as this has been widely illustrated in the previous sections (see also the extensive work of Thomas Sederberg on this topic, e.g. [172]). Consequently, methods for solving polynomial systems are required. The aim of this section is to give a quick overview of such methods. In order to be as much concrete as possible we mention the following two typical problems of shape interrogation that can easily be reduced to polynomial system solving: Singularity detection. An important problem in CAGD is the detection of singularities of a 3D-surface. If an algebraic surface S is given implicitly by P (x, y, z) = 0 (that is

S = {(x, y, z) ∈ R 3 ; P (x, y, z) = 0}), a point (a, b, c) on S is singular if ∂P ∂x (a, b, c) = ∂P ∂y (a, b, c) = ∂P ∂z (a, b, c) = 0.
It is then clear that the singular points of S are the common roots of the polynomials P, ∂P ∂x , ∂P ∂y , ∂P ∂z .

Computation of intersection points. Given two parameteric curves, one would like to compute their intersection points. By implicitizing one of the two curves this problem is reduced to finding the real roots of a univariate polynomial which is obtained by substituting the parameterization of a curve into the implicit equation of the second one. Similar approaches can be used to compute curve/surface intersection points and more generally to compute a parameterization of an intersection surface/surface curve. Though we are manipulating objects in dimension 3, the polynomial systems that we consider might involve more that 3 variables. For instance, the intersection of 2 parametric surfaces involve 4 variables. Therefore, we are not going to restrict the number of variables in the methods that we are going to describe. Hereafter, the variables will be denoted x 1 , . . . , x n . However, since these systems come from real geometric modeling problems, we will consider only polynomials with real coefficients.

Since the problem of solving polynomial equations goes back to the ancient Greeks and Chinese, it is not surprising to see that a large number of methods exists to handle this task. Several families of solvers can however be identified:

• Analytic solvers exploit the value of a functional f = (f 1 , . . . , f m ) and its derivatives in order to converge to a solution or all the solutions of f = 0. Typical examples are Newton like methods, minimization methods, etc. • Subdivision methods use an exclusion criterion to remove a domain if it does not contain a root of f = 0 or refine the search in sub-domains otherwise. These solvers are often used to isolate the real roots in a given domain.

• Algebraic solvers exploit the known relations between the unknowns. They are based on polynomial manipulations and involve effective algebraic geometry tools.

We are going to focus essentially on the two last families, which yield information on all the roots (resp. in a fixed domain).

Subdivision methods

The methods that we describe in this section, exploit the properties of Bernstein's basis for representing univariate and multivariate polynomials. The Bernstein polynomials are ubiquitous in geometric modeling. The representation of a polynomial in the Bernstein basis is known to be numerically more stable than the monomial basis representation [START_REF] Farouki | On the optimal stability of the bernstein basis[END_REF][START_REF] Farouki | On the numerical condition of polynomials in bernstein form[END_REF]. It has a direct geometric meaning, in terms of control points and useful properties such that the convex hull and the variation diminishing property. These properties in conjunction with the subdivision nature of Bernstein's polynomials explain the large variety of algorithms proposed until today for solving univariate polynomials, starting with Lane and Riesenfeld The situation in the multivariate case has not been studied so extensively. Two main sub-families coexist: a first family which are based on subdivision techniques like [START_REF] Elber | Geometric constraint solver using multivariate rational spline functions[END_REF]171]; a second family of solvers which are based on reduction techniques as [173]. We briefly describe these approaches, starting with univariate polynomials. For more details, see [136].

Univariate polynomials

Any polynomial f (x) ∈ IR[x], of degree d, can be represented as f (x) = d i=0 b i B d i (x) where B d i (x) = d i (1 -x) d-i x i , i = 0, . . . , d
is the Bernstein basis associated to the interval [0, 1]. Similarly, we will say that a sequence b represents the polynomial f on the interval [r, s] if: A first property of this representation is that the derivative f ′ of f , is represented by the control coefficients:

f (x) = d i=0 b i d i 1 (s -r) n (x -r) i (s -x) d-i . The polynomials B i d (x; r, s) := d i 1 (s-r) n (x -r) i (s -x) n-i form
d∆b := d(b i+1 -b i ) 0 i d-1 .
Another fundamental algorithm that we will use on such a representation is the de Casteljau algorithm [START_REF] Farin | Curves and Surfaces for Computer Aided Geometric Design[END_REF]. It allows us to subdivide the representation of f into the two sub-representations on the intervals [r, (1 -x)r + xs] and [(1 -x)r + xs, s]. It requires at most 2d(d + 1) arithmetic operations. For a more detailed list of properties of this representation, we refer for instance to [START_REF] Farin | Curves and Surfaces for Computer Aided Geometric Design[END_REF]. A simple but interesting property that we are going to use is the following: Theorem 1 (Descartes rule). The number of real roots of the polynomial

f (x) = b i B i d (x; r, s) in ]r,

Multivariate polynomials

By a direct extension to the multivariate case, any polynomial f (x 1 , . . . , x n ) ∈ IR[x 1 , . . . , x n ] = IR[x] of degree d i in the variable x i , can be decomposed as:

f (x 1 , . . . , x n ) = d1 i1=0 • • • dn in=0 b i1,...,in B i1 d1 (x 1 ; r 1 , s 1 ) • • • B in dn x(x n ; r n , s n ) where ( B i1 d1 (x 1 ; r 1 , s 1 ) • • • B in dn (x n ; r n , s n )) 0≤i1≤d1,.
..,0≤in≤dn is the tensor product Bernstein basis on the domain De Casteljau's algorithm applies in each of the direction x i , , i = 1, . . . , n so that we can split this representation in these directions. This algorithm can be used either to split the domain or to restrict the representation to a sub-domain. For a univariate polynomial of degree d, this costs 2 (d + 1)d arithmetic operations. For a multivariate polynomial of degree d i in x i , we check that this restriction operation costs 2

D := [a 1 , b 1 ] × • • • × [r n , s n ] ⊂ IR
n i=1 d i n i=1 (d i + 1) = O(d n+1
), where d = max{d 1 , . . . , d n }. Thus, as the dimension and the degree increase, a good method to isolate the roots, should consider carefully when to apply this reduction operation, in order to save the computation time.

Notice that the univariate Bernstein representation also extends to the socalled triangular Bernstein basis representation. This representation can also be used in our approach, but we will concentrate on the tensor product one. For any f ∈ IR[x] and j = 1, . . . , n, let

m j (f ; x j ) = dj ij =0 min {0≤i k ≤d k ,k =j} b i1,...,in B ij dj (x j ; r j , s j ) M j (f ; x j ) = dj ij =0 max {0≤i k ≤d k ,k =j} b i1,...,in B ij dj (x j ; r j , s j ).
We have the following property: for any u = (u 1 , . . . , u n ) ∈ D, and any j = 1, . . . , n, we have

m(f ; u j ) ≤ f (u) ≤ M (f ; u j ).
As a direct consequence, for any root u = (u 1 , . . . , u n ) of the equation f (x) = 0 in the domain D, we have µ j ≤ u j ≤ µ j where • µ j (resp. µ j ) is either a root of m j (f ; x j ) = 0 or M j (f ; x j ) = 0 in [r j , s j ] or r j (resp. s j ) if m j (f ; x j ) = 0 (resp. M j (f ; x j ) = 0) has no root on [r j , s j ],

• m j (f ; u) ≤ 0 ≤ M j (f ; u) for u ∈ [µ j , µ j ].
This transforms the problem of approximating the real roots of multivariate polynomials into problems on univariate polynomials.

Univariate Root Solver

Descartes rule (see theorem 1) yields a simple subdivision algorithm, which splits the domain when the number of sign variation of the control coefficients is bigger than 2. In the presence of a multiple root, the number of sign changes of a representation on any interval containing a multiple root is bigger than 2, and the algorithm splits the box until its size is smaller than a given ε. A detailled analysis of the behavior of the algorithm, has been carried out, using a partial inverse of Descartes rule given by the two circles theorem. See [140],

[137], [START_REF] Ioannis | Real algebraic numbers: Complexity analysis and experimentations[END_REF], [START_REF] Eigenwillig | Almost tight recursion tree bounds for the descartes method[END_REF].

This algorithm yields, in the presence of simple roots, an interval isolating the roots. But usually in practice, we are interested in approximating this root within a given precision ε. In order to approximate the isolated roots within ε, further steps of bisection may be required, either using de Casteljau's algorithm, Newton-like methods, or variants such as in [165].

Multivariate root finding

In this section, we consider a system of s equations in n variables f 1 (x 1 , . . . , x n ) = 0, . . . , f s (x 1 , . . . , x n ) = 0 with coefficients in IR, that we will also denote by f(x) = 0. We are looking for an approximation of the real roots of f(x) = 0 in the domain

D = [r 1 , s 1 ]× • • • × [r n , s n ],
within a precision ǫ. The general framework of the families of algorithms that we will consider consists in (1) applying a preconditioning step on the equations, (2) in reducing the domain, and (3) if the reduction ratio is too small, in splitting the domain, until the size of the domain is smaller than a given epsilon. The solvers that we will consider are parameterized by the • Preconditioner, that is, a transformation of the initial system f = 0 into a system M f = 0 (with an M invertible matrix), which has a better numerical behavior. We consider -Global transformation which aims at increasing the distance between the equations, considered as vectors of coefficients, -Local straightening which multiplies by the inverse of the Jacobian matrix at the center of the box, if it exists (it applies only for square systems).

• Reduction strategy, that is, the technique used to reduce the initial domain, for searching the roots of the system. We consider -Convex hull reduction as described in [173].

-Extreme root reduction, which consists in computing the first (resp. last) root of the polynomial m j (f k ; u j ), (resp. M j (f k ; u j )), in the interval [r j , s j ]. The improvement compared with the previous approach can be substantial (see Figure 8.1). • Subdivision strategy, that is, the technique used to subdivide the domain, in order to simplify the forthcoming steps, for searching of roots of the system. Here some simple rules that can be used to subdivide a domain, either in the parameter domain, or in the image.

This family of algorithms has been implemented in the C++ library synaps5 and compared on several benchmarks. It appears that the strategy of local preconditionning with an emphasis on reduction is significantly better than the other strategies. The performances of such a solver are very good.

Algebraic methods

We now turn to algebraic methods for solving polynomial systems. We will mainly discuss the ones based on resultant matrix constructions, but also mention a method based on normal form computations which generalizes the well-known concept of Gröbner basis.

Resultant-based methods

Resultant theory. The theory of resultant is devoted to the study of conditions on the coefficients of an over-determined system, to have a solution in a fixed variety. The typical situation is the case of a system of n + 1 equations in a space of dimension n of the form:

f c :=        f 0 (x) = k0 j=0 c 0,j ψ 0,j (x) . . . f n (x) = kn j=0 c n,j ψ n,j (x)
where c = (c i,j ) are parameters, x is a point of a variety X of dimension n, and the vector L i = (ψ i,j ) j=0,...,ki is a regular map [START_REF] Harris | Algebraic geometry[END_REF] from X to the projective space P ki independent of c. The elimination problem consists, in this case, in finding necessary (and sufficient) conditions on c such that the system f c = 0 has a solution in X. From a geometric point of view, we look for the values of parameters c = (c i,j ) such that there exists x ∈ X with ki j=0 c i,j ψ i,j (x) = 0 for i = 0, . . . , n. It turns out that these parameters are exactly the zero locus of a unique polynomial equation (defined up to the multiplication by a non-zero constant) in c which is called the resultant of f 0 , . . . , f n and is denoted by Res X (f c ). It is a quite attractive object because one can compute it through some matrix constructions.

Construction of resultant matrices. In order to construct a non-trivial multiple of Res X (f c ), we apply the following strategy. A vector L(c, x) of polynomials in Z[c][x], where x denotes a coordinate system of a projective space containing X, is constructed in such a way that (1) the polynomial entries of L(c, x) are generically independent, (2) the set v(x) of monomials (or polynomials) in x needed to decompose all the polynomials of L(c, x) has the same size than L(c, x), and

(3) the polynomials in L(c, x) vanish when the input system has a common root in the variety X.

The coefficient matrix of the polynomials entries of L(c, x) with respect to the set v(x) yields a matrix S(c) whose entries are in Z[c]. Its determinant is nonzero, according to the point (1).

Above and hereafter, we use the term generic, which means that the property that we are considering is true on an open subset of the coefficient space for c. Our aim is to construct matrices S(c), which can be used for generic systems of a certain class of polynomial equations. In practice, the problem is not posed in these terms. We are given a system of equations and it may happen that the construction we are considering yields a degenerate matrix S(c). In this case, the system is not generic for the resultant formulation and we have to chose another class of systems for which we are in a generic position. This explains why a lot of different types of resultant formulations have been studied; we will give a list in a moment.

By construction, we have v(x) t S(c) = L(c, x) t . Thus, according to the point (2), if ζ is a common root of a specialized system f c0 = 0, we have

L(c 0 , ζ) = 0 and v(ζ) t S(c 0 ) = 0. If (generically) v(ζ)
is not zero at a common root ζ ∈ X of f c0 = 0, we deduce that det S(c) vanishes when the system has a common root in X. Therefore det S(c) is a non-trivial multiple of its equation, that is of the resultant Res X (f c ).

A usual way to construct these resultant matrices (which extends Sylvester's construction for the classical resultant of two univariate polynomial), consists in choosing for L, a list of monomial multiples of the polynomials f i . In this case, the matrix S(c) is the matrix of a map of the form

S : x E0 × • • • × x En → x F (q 0 , . . . , q n ) → g = n i=0
q i f i where x Ei is the vector subspace generated by a specific set of monomials x Ei . The entries of S(c) are filled as follows: every column of S is indexed by an element of some {E i } i=0,...,n and every row by an element of F ; equivalently, the columns and rows are indexed by the monomials of q i and the monomials of g, respectively. The coefficient in the row corresponding to β ∈ F and in the column corresponding to α ∈ E j is the coefficient of x β in x α f j . The coefficients of monomials which do not explicitly appear in x α f j have a zero entry. Thus the matrix S(c) is divided into blocks (S 0 , . . . , S n ), where each block S i depends linearly on the coefficients of f i . Notice that such resultant matrices have a quasi-Toeplitz structure which can be exploited to accelerate many computations with them, by almost one order of magnitude in terms of the matrix dimension. This relies on FFT for fast multivariate polynomial and dense structured matrix arithmetic; a smaller acceleration is achieved by applying Karatsuba's divide-and-conquer arithmetic [135,[START_REF] Emiris | Symbolic and numeric methods for exploiting structure in constructing resultant matrices[END_REF].

Different resultant formulations. We recall briefly different resultant formulations which can be used in geometric problems. A detailed description is beyond the scope of this book. The formulation will be chosen according to the geometric properties of the problem to solve.

• multivariate resultants: They correspond to the classical case studied in [115,187]. Here X is the projective space P n , L i is the vector of all monomials of a fixed degree d i , and the function f i is a generic homogeneous polynomial of degree d i . The most used resultant matrix in this context is the Macaulay's construction [115] which can be seen as an extension of the Sylvester's method to the multivariate case. However, some other multivariate resultant matrices have been developed (see e.g. [START_REF] Jouanolou | Formes d'inertie et résultant: un formulaire[END_REF][START_REF] Dickenstein | Explicit formulas for the multivariate resultant[END_REF]).

• toric (or sparse) resultants: They have been introduced in [START_REF] Kapranov | Chow polytopes and general resultants[END_REF], then developed in [START_REF] Gelfand | Discriminants, resultants, and multidimensional determinants. Mathematics: Theory & Applications[END_REF]. It takes into account the monomial support of the input polynomials of a system, and not only their respective degree. Thus it is possible to work with polynomials having negative exponents, that is Laurent polynomials. Methods for constructing toric resultant matrices can be found in [START_REF] Canny | An algorithm for the Newton resultant[END_REF][START_REF] Emiris | A practical method for the sparse resultant[END_REF][START_REF] Martin | Macaulay style formulas for sparse resultants[END_REF].

• Residual resultants: In many situations coming from practical problems, the equations have common zeroes which are independent of the parameters of the problems, and which we are not interested in. The residual resultant constructions has been designed to take into account these degenerate cases. It is described, as well as matrix construction, in [START_REF] Busé | Resultant over the residual of a complete intersection[END_REF][START_REF] Busé | Residual resultant over the projective plane and the implicitization problem[END_REF][START_REF] Busé | Using projection operators in computer aided geometric design[END_REF]. A more general construction has been developed in [START_REF] Busé | Generalized resultants over unirational algebraic varieties[END_REF] whose associated matrix construction is the so-called called Bezoutian matrix.

Solving polynomial systems via eigenvalues computations. Let f 0 (x), f 1 (x), . . ., f n (x) be polynomials in n variables x = (x 1 , . . . , x n ). By choosing an adapted resultant formulation one can construct a resultant matrix S associated to this system. It turns out that this matrix can be divided into four blocs S = S 00 S 01 S 10 S 11 and that the Schur complement S 00 -S 01 S 11 -1 S 10 is nothing but the matrix of the multiplication map by f 0 (x) in a canonical basis of the quotient ring R[x]/(f 1 , . . . , f n ). The point is that the eigenvalues of such a multiplication matrix are particularly interesting, they are the evaluation of f 0 at the common root of f 1 and f 2 . If f 0 is a linear form one can thus easily solve the polynomial system f 1 (x) = f 2 (x) = 0. Solving polynomial systems by hiding a variable. Another approach to solve a system of polynomial equations consists in hiding a variable (that is, in considering one of variables as a parameter), and in searching the value of this hidden variable for which the system has a solution. Typically, if we have n equations f 1 = 0, . . . , f n = 0 in n variables, we "hide" a variable, say x n , and apply one of resultant constructions described before to the over-determined system f 1 = 0, . . . , f n = 0 in the n -1 variables x 1 , . . . , x n-1 and a parameter

x n . This leads to a resultant matrix S(x n ) with polynomial entries in x n . It can be decomposed as

S(x n ) = S d x d n + S d-1 x d-1 n + • • • + S 0 ,
where S i has coefficients in R and the same size as S(x n ). We look for the values ζ n of x n for which the system has a solution ζ ′ = (ζ 1 , . . . , ζ n-1 ) in the corresponding variety X ′ (of dimension n -1) associated with the resultant formulation. This implies that

v(ζ ′ ) t S(ζ n ) = 0, (7)
where v(ζ ′ ) is the vector of monomials indexing the rows of S evaluated at ζ ′ . Conversely, for generic systems of the corresponding resultant formulation there is only one point ζ ′ above the value ζ n . Thus the vectors v satisfying S(ζ n ) t v = 0 are scalar multiples of v(ζ ′ ). From the entries of these vectors, we can deduce the other coordinates of the point ζ ′ . This will be assumed hereafter 6 . The relation [START_REF] Beier | Highlight-line algorithm for realtime surface quality assessment[END_REF] implies that v(ζ ′ ) is a genearlized eigenvector of S t (x n ). Computing such vectors can be transformed into the following linear generalized eigenproblem

           0 I • • • 0 . . . . . . . . . . . . 0 • • • 0 I S t 0 S t 1 . . . S t d-1      -ζ n       I 0 • • • 0 0 . . . . . . . . . . . . . . . I 0 0 • • • 0 -S t d             w = 0. ( 8
)
The set of eigenvalues of ( 8) contains the values of ζ n for which (7) has a solution. The corresponding eigenvectors w are decomposed as w = (w 0 , . . . , w d-1 ) so that the solution vector v(ζ ′ ) of ( 7) is

v(ζ ′ ) = w 0 + ζ n w 1 + • • • + ζ d-1 n w d-1 .

Normal forms

Gröbner basis is a powerful tool to handle a lot of computations on polynomial systems. However their construction is not numerically stable, they may introduce artificial discontinuities due to the choice of a monomial order. A recent generalization of this notion has been proposed in [134,138,139]. It is based on a new criterion which gives a necessary and sufficient condition for a projection onto a vector subspace of R to be a normal form modulo the ideal I. More precisely we have:

Theorem 3. Let B be a vector space in R = R[x 1 , . . . , x n ] connected to the constant polynomial 1 7 . If B + is the vector subspace generated by B ∪ x 1 B ∪ {1, . . . , n}.

. . . ∪ x n B, N : B + → B is a linear map such that N is the identity on B, we define for i = 1, . . . , n, the maps

M i : B → B b → M i (b) := N (x i b).
The two following properties are equivalent:

1. For all 1 ≤ i, j ≤ n, M i • M j = M j • M i . This leads to a completion-like algorithm which starts with the linear subspace K 0 generated by the polynomials f 1 , . . . , f m , that we want to solve and iterates the construction K i+1 = K + i ∩ L, where L is a fixed vector space. We stop when K i+1 = K i . See [134,138,186,139] for more details. This approach allows us to fix first the set of monomials on which we want to do linear operations and thus to treat more safely polynomials with approximate coefficients. It can be adapted very naturally to Laurent polynomials, which is not the case for Gröbner basis computation. Moreover it can be specialized very efficiently to systems of equations for which the basis of A is known a priori, such as in the case of a complete projective intersection [138]. Let us see how we can deduce the roots from this normal form computation. For this purpose, we will use the properties of the operators of multiplication by elements of A = R(f 1 , .., f m ). For any a ∈ A, we define We also consider its transpose operator

M t a : A → A Λ → M t a (Λ) = Λ • M a ,
where the dual space A is the set of R-linear forms from A to R. The matrix of M t a in the dual basis of a basis B of A is the transpose of the matrix of M a in B. The multiplication operators can be computed using a normal form algorithm, as described above.

Hereafter, x E = (x α ) α∈E denotes a monomial basis of A (for instance obtained by a Gröbner basis). Then any polynomial can be reduced modulo (f 1 , . . . , f m ) to a linear combination of monomials of x E .

The matrix approach to solve polynomial systems is based on the following fundamental theorem [START_REF] Auzinger | An elimination algorithm for the computation of all zeros of a system of multivariate polynomial equations[END_REF], [133]: 2. The common eigenvectors of (M t a ) a∈A are (up to a scalar) the evaluations 1 ζ1 , . . . , 1 ζ d .

Since x E = (x α ) α∈E is a basis of A, the coordinates of 1 ζi in the dual basis of x E are (ζ α i ) α∈E . Thus if x E contains 1, x 1 , . . . , x n (which is often the case), we can deduce directly all the coordinates of the roots. We have the following algorithm:

Algorithm 1 Solving in the case of simple roots. Let a ∈ A such that a(ζ i ) = a(ζ j ) for i = j (which is generically the case) and M a be the matrix of multiplication by a in the basis x E = (1, x 1 , . . . , x n , . . .) of A. .

The set of output points ζ contains the simple roots (i.e. roots with multiplicity 1) of f = 0, since for such a root the eigenspace associated to the eigenvalue a(ζ) is one-dimensional and contains 1 ζ . But as we will see in the next example, it can also yield in some cases the multiple roots.

In order to compute exactly the set of roots counted with their multiplicity, we use the following result. It is based on the fact that commuting matrices share common eigenspaces. [133,135,[START_REF] Corless | A reordered Schur factorization method for zero-dimensional polynomial systems with multiple roots[END_REF].

Theorem 5. There exists a basis of A such that for all a ∈ A, the matrix of M a in this basis is of the form 1. Compute a (Schur) decomposition P such that the matrices T i = PM xi P -1 , i = 1, . . . , n, are upper-triangular. 2. Compute and output the diagonal vectors t i = (t 1 i,i , . . . , t n i,i ) of triangular matrices T k = (t k i,j ) i,j .

output: Z(I) = {t i : i = 1, . . . , dim R (A)}.

The first step in this algorithm is performed by computing a Schur decomposition of M l (where l is a generic linear form) which yields a matrix P of bases change. Then we compute the triangular matrices T i = PM xi P -1 , i = 1, . . . , n, since they commute with M l . An implementation by Ph. Trébuchet of this algorithm is available in the synaps library (see solve(L,newmac<C>())).

Conclusion

Shape interrogation methods are still of increasing interest in geometric modeling as well as in computer graphics. Originating 20 years ago from CAD/CAM applications where "class A" surfaces are required and no surface imperfections are allowed, shape interrogation has become recently an important tool for various other types of surface representations such as triangulated or polygonal surfaces, subdivision surfaces, and algebraic surfaces. In this chapter, we presented the state-of-the-art of shape interrogation methods including methods for detecting surface imperfections, surface analysis tools and methods for visualizing intrinsic surface properties. Furthermore we focused on stable numerical and symbolic solving of algebraic systems of equations, a problem that arises in most shape interrogation methods. Nevertheless, many issues are still open promising intensive research in various areas of shape interrogation.

Let us focus on some of them now. Discrete geometry representations are frequently used in many applications, especially for shapes acquired from real-world objects. Typically, surfaces are approximated by polygonal meshes, and we showed how to estimate differential properties for piecewise linear surfaces. Various methods exist so far, and recent approaches prove approximation and convergence properties. The design of robust methods coming with certain guarantees is still an area of active research.

In the area of algebraic and numerical polynomial system solvers, that provide one of the basic tools for shape interrogation methods, improvements are indispensable. Many critical problems in Computer Aided Geometric Deisgn, such as shape interrogation, are reduced to finding the zero set of a system of polynomial equations. Several root-finding methods for polynomial systems exist, even if we mainly presented resultant-based methods and subdivision methods. A wide choice of techniques and algorithms to solve polynomial systems are thus now available, but as a main drawback, all of these methods have difficulties in handling roots with high multiplicities (or clusters of roots). They all have performance deterioration, lack of robustness in numerical computation and round-off errors during floating point arithmetic in such a situation. It is hence a crucial objective and an active research area to improve root-finding methods in this case which often occurs in practice.

In the area of symbolic curve and surface interrogation future/open problems can be addressed. Volumetric data sets are now used in many applications and serves as a prime candidate representation in medical applications. The extension of curve and surface interrogation methods to support volumetric representation, either as iso-surfaces in the volumes or direction analyze differential properties are highly desired. The degrees of many of these rational fields such as the Gaussian curvature, K, or the mean curvature square, M 2 , are high. Methods to robustly handle these fields, in a more stable way, could further improve the quality of the result.

Fig. 1 .

 1 Fig. 1. Frenet frame.

Fig. 2 .

 2 Fig. 2. Visualization of mean curvature M (a) and Gaussian curvature K (b) estimated on the Feline triangle mesh. Here, red, green and blue denote positive, zero and negative values, respectively, and lighting is enabled. (c) and (d) show the maximum curvature with lines of curvature on the Mannequin mesh.

Fig. 3 .

 3 Fig. 3. First order shape analysis by simulating light reflection.

Fig. 4 .

 4 Fig. 4. Pattern of computed reflection lines and isophotes on NURBS surfaces.

Fig. 5 .

 5 Fig. 5. Isophotes for curvature discontinuity.

Fig. 6 .

 6 Fig. 6. Bicubic surface patch with line of vanishing Gaussian curvature (left). Orthotomic analysis (right).

Fig. 7 .

 7 Fig. 7. Local shape of normal section curve is defined by curvature.

Fig. 8 .

 8 Fig. 8. Gaussian curvature determines local shape of surface. Left images: convex and concave regions (K > 0). Middle: saddle-shaped region (K < 0). Middle-right: a parabolic point (K = 0) Right: a region consisting of parabolic points.

Fig. 9 .

 9 Fig. 9. Gauss map and geometric meaning of Gaussian curvature.

Fig. 10 .

 10 Fig. 10. Typical singularities of the focal surface. From left to right: cuspidal edge (rib), swallowtail, pyramid, purse. At the swallowtail singularity the rib has a cusp. The pyramid and purse correspond to the umbilical point on the surface. The vertical lines at the bottom images are the surface normals at the corresponding umbilics.

Fig. 11 .

 11 Fig.11. Left: zoo of distance functions; thin lines are used to sketch typical profiles of the surface functions defined by the distance from a given point to the surface points. Center: the skeleton (blue), caustic (yellow), ridge (red) and an osculating sphere (brown) at a ridge point of the elliptic paraboloid. Right: schematic illustration of relationships between the cut locus, focal surface and ridges.

Fig. 12 .

 12 Fig. 12. Second order surface analysis with generalized focal surfaces.

Fig. 13 .

 13 Fig. 13. Color codings of Gaussian curvature.

Fig. 14 .

 14 Fig. 14. Gaussian curvature isoline. Left: parabolic lines. Right: isolines corresponding to different constant Gaussian curvature values.

•

  surfaces in implicit form and isosurfaces of 3D images [158, 129, 128, 184, 182, 10, 13]; • surfaces approximated by polygonal meshes [113, 9, 188, 82, 177, 26, 25, 145, 81]; • height data [65, 96, 95, 109];• surface given in parametric form[START_REF] Hosaka | Modeling of Curves and Surfaces in CAD/CAM[END_REF][START_REF] Hartmann | On the curvature of curves and surfaces defined by normalforms[END_REF].

Fig. 15 .

 15 Fig. 15. Various types of ridges detected on smooth surfaces. The images are taken from [13].

Fig. 16 .

 16 Fig. 16. The crest lines detected on various surfaces approximated by dense triangle meshes.

= 2 π ln κ 2 1 + κ 2 2

 22 • A ridge is tangent to the line of curvature of the same color. These surface landmarks corresponds to the swallowtail singularities of the focal surface.• A ridge crosses a ridge of other color. In[183] it was suggested to use these landmarks for 3D image registration. • A ridge crosses the parabolic line of the same color. The Gauss map has the the so-called pleat singularity at such a point[101].Koenderink [101] introduced two curvature-based measures of surface curvature: the curvedness C

Fig. 17 .

 17 Fig.17. Left: a free-form B-spline surface is presented, after being subdivided into convex (red), concave (green), and hyperbolic regions (yellow). The parabolic lines (white) separate the regions. Right: presents the function of K(u, v) (in yellow) and its zero set (the parabolic lines).

Fig. 18 .

 18 Fig.18. Left: a free-form B-spline surface is presented, after being subdivided into regions of different levels of κ 2 1 + κ 2 2 . Right: presents the rational surface κ 2 1 + κ 2 2 and its contouring (in white) at the different levels.

  [105], up to the Bezier clipping methods initiated by Nishita and al [144]. They combine a global control on the domain where the roots are searched with local and efficient refinements.

  the Bernstein basis on [r, s]. Hereafter, we are going to consider the sequence of values b together with the corresponding interval [r, s], as representing of our polynomial f .

  s[ is bounded by the number V (b) of sign changes of b = (b i ) i=0..n , and is equal modulo 2. As a consequence, if V (b) = 0 there is no root in ]r, s[ and if V (b) = 1, there is one root in ]r, s[. Another interesting property of this representation is the following (see e.g. [53], [164]): Theorem 2 (Convex hull). Let b = (b i ) i=0,...,d be the control coefficients of f (x) on the interval [r, s] and c = (c i ) i=0,...,d the corresponding control points. The graph {(x, f (x)); t ∈ [r, s]} is in the convex hull of the control points c.

  n and b = (b i1,...,in ) 0≤i1≤d1,...,0≤in≤dn are the control coefficients of f on D. The polynomial f is represented in this basis by the n th order tensor of control coefficients b = (b i1,...,in ) 0≤i≤d1,0≤j≤d2,0≤k≤d3 . The size of D, denoted by |D|, is |D| = max{|s i -r i |; i = 1, . . . , n}.

Fig. 21 .

 21 Fig. 21. Convex hull vs. extreme roots

2 .

 2 R = B ⊕ I, where I is the ideal generated by the kernel of N If this holds, the B-reduction along ker(N ) is canonical.

M

  a : A → A b → M a (b) := a b.

Theorem 4 .

 4 Assume that Z(I) = {ζ 1 , . . . , ζ d }. We have 1. Let a ∈ A. The eigenvalues of the operator M a (and its transpose M t a ) are a(ζ 1 ), . . . , a(ζ d ).

1 . 2 .

 12 Compute the eigenvectors Λ = (Λ 1 , Λ x1 , . . . , Λ xn , . . .) of M t a . For each eigenvector Λ with Λ 1 = 0, compute and output the point ζ = Λx 1 Λ1 , . . . , Λx n Λ1

Algorithm 2

 2 Solving by simultaneous triangulation. input: Matrices of multiplication M xi ,i = 1, . . . , n, in a basis of A.

  . Saddle points appear as "passes". The contour lines only cross in the exceptional case of a contour at the precise level of a saddle point.[141] describes systematically the distribution of other critical points on a surface. A disadvantage of contour lines is the fact that they are costly to compute. Several surface contouring methods exist, which are sometimes depending of the specific surface formulation [152, 169, 108]. Hartwig and Nowacki[START_REF] Hartwig | Isolinien und schnitte in coonschen flächen[END_REF] propose to subdivide the surface into sufficient small pieces which are then approximated by bilinear surfaces. Then the contour lines can easily be computed. Iso-contouring is the technique of extracting constant valued curves and surfaces from 2D and 3D scalar fields. Interactive display and quantitative interrogation helps understanding the overall structure of a scalar field and its evolution over time. Traditional iso-contouring techniques examine each cell of a mesh to test for intersection with the iso-contour of interest. For an overview see[168]. Extraction of isosurfaces from 3D scalar field is generally be done by the Marching Cubes algorithm and its variants[111, 143,[START_REF] Chan | A new tetrahedral tesselation scheme for isosurface generation[END_REF]

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

Notice however that this genericity condition can be relaxed by using duality, in order to compute the points ′ above ζn (when they form a zero-dimensional fiber) from the eigenspace of S(ζn).