Chapter 41
Approximation of Ridges and Umbilics on Triangulated Surface Meshes

Table of Contents

This chapter describes the CGAL package for the approximating the ridges and umbilics of a smooth surface discretized by a triangle mesh. Given a smooth surface, a ridge is a curve along which one of the principal curvatures has an extremum along its curvature line. An umbilic is a point at which both principal curvatures are equal. Ridges define a singular curve - i.e. a self-intersecting curve, and umbilics are are special points on this curve. Ridges are curves of extremal curvature and therefore encode important informations used in segmentation, registration, matching and surface analysis. Based on the results of the article [CP05c], we propose algorithms to identify and extract different parts of this singular ridge curve as well as umbilics on a surface given as a triangulated surface mesh. Differential quantities associated to the mesh vertices are assumed to be given for these algorithms; such quantities may be computed by the package Estimation of Local Differential Properties of Sampled Surfaces via Polynomial Fitting.

Note that this package needs the third party libraries Lapack and Blas for linear algebra operations.

41.0.1   Overview

Section 41.1 presents the basics of the theory of ridges and umbilics on smooth surfaces. Sections 41.2 and 41.3 present algorithms for approximating the ridges and umbilics (of a smooth surface) from a triangle mesh. Section 41.4 gives the package specifications, while example calls to functions of the package are provided in Section 41.5.

41.1   Ridges and Umbilics of a Smooth Surface

For a detailed introduction to ridges and related topics, the reader may consult [HGY$$⁺99, Por01], as well as the following survey article [CP05b]. In the sequel, we just introduce the basic notions so as to explain our algorithms. Consider a smooth embedded surface, and denote $$k₁ and $$k₂ the principal curvatures, with $$k₁ k₂. Umbilics are the points where $$k₁=k₂. For any non umbilical point, the corresponding principal directions of curvature are well defined, and we denote them $$d₁ and $$d₂. In local coordinates, we denote $$ , the inner product induced by the ambient Euclidean space, and $$dk₁, $$dk₂ the gradients of the principal curvatures. Ridges, illustrated on figure 41.2 for the standard ellipsoid, are defined by:

Definition.   A non umbilical point is called

• a max ridge point if the extremality coefficient $$b₀= dk₁,d₁ vanishes, i.e. $$b₀=0.

• a min ridge point if the extremality coefficient $$b₃= dk₂,d₂ vanishes, i.e. $$b₃=0 ¹.

The previous characterization of ridges involves third-order differential properties. Using fourth-order differential quantities, a ridge point can further be qualified as elliptic if it corresponds to a maximum of $$k₁ or a minimum of $$k₂, or hyperbolic otherwise. Hence we end with four types of ridges, namely : MAX_ELLIPTIC_RIDGE, MAX_HYPERBOLIC_RIDGE, MIN_ELLIPTIC_RIDGE, MIN_HYPERBOLIC_RIDGE, see figure 41.2. Also of interest are the crest lines, a crest line being an elliptic ridge which is a maximum of $$max$$(|k₁|,|k₂|). Crest lines form a subset of elliptic ridges, and can be seen as the visually most salient curves on a surface. Hence we provide two additional ridge types : MAX_CREST_RIDGE and MIN_CREST_RIDGE, see figure 41.

At any point of the surface which is not an umbilic, principal directions $$d₁, d₂ are well defined, and these (non oriented) directions together with the normal vector $$n define two direct orthonormal frames. If $$v₁ is a unit vector of direction $$d₁ then there exists a unique unit vector $$v₂ so that $$(v₁,v₂,n) is direct, that is has the same orientation as the canonical basis of the ambient $$3d space (and the other possible frame is $$(-v₁,-v₂,n)). In the coordinate systems $$(v₁,v₂,n), the surface has the following canonical form, known as the Monge form :

The Taylor expansion of $$k₁ (resp. $$k₂) along the max (resp. min) curvature line going through the origin and parameterized by $$x (resp. $$y) are:

$$ k₁(x) = k₁ + b₀x + (P₁)/(2(k₁-k₂))x² +h.o.t,       P₁= 3b₁²+(k₁-k₂)(c₀-3k₁³).

$$ k₂(y) = k₂ + b₃y + (P₂)/(2(k₂-k₁))y² +h.o.t,       P₂= 3b₂²+(k₂-k₁)(c₄-3k₂³).

Notice also that switching from one to the other of the two afore-mentioned coordinate systems reverts the sign of all the odd coefficients on the Monge form of the surface.

Hence ridge types are characterized by

• max ridge if $$b₀=0
• MAX_ELLIPTIC_RIDGE if $$b₀=0 and $$P₁<0
• MAX_HYPERBOLIC_RIDGE if $$b₀=0 and $$P₁>0
• min ridge if $$b₃=0
• MIN_ELLIPTIC_RIDGE if $$b₃=0 and $$P₂<0
• MIN_HYPERBOLIC_RIDGE if $$b₃=0 and $$P₂>0
• MAX_CREST_RIDGE if $$b₀=0 and $$P₁<0 and $$|k₁|>|k₂|
• MIN_CREST_RIDGE if $$b₃=0 and $$P₂<0 and $$|k₂|>|k₁|

As illustrated on Fig. 41.3 and 41.4, the patterns made by curvature lines around an umbilic can be characterized using the notion of index associated to the principal directions - see also [CP05b]. As depicted on Fig. 41.3, consider a small circuit $$C around the umbilic, and a point $$p C. Starting from an initial orientation $$u of a tangent vector to the curvature line through point $$p, propagate by continuity this orientation around the circuit. The index is defined by the angle swept by $$u around this revolution, normalized by $$2. On our example, the index is thus 1/2.

If the index of the principal direction field is $$1/2 then it is called a ELLIPTIC_UMBILIC, if it is $$-1/2 it is called a HYPERBOLIC_UMBILIC. Otherwise the umbilic is qualified NON_GENERIC_UMBILIC.

41.2   Approximating Ridges on Triangulated Surface Meshes

Our method aims at reporting ridges as polygonal lines living on the mesh. It assumes differential quantities are available for each vertex of the mesh (principal curvatures and directions together with third order quantities $$b₀, b₃ and optionally fourth order quantities $$P₁, P₂). These differential quantities may be computed for the smooth surface the mesh is inscribed in (analytically or using approximation methods), or may be estimated for a mesh given without reference to a underlying smooth surface. Although the ridge approximation algorithm is the same in both cases, one cannot ambition to ask for the same certificates. This distinction calls for the notion of compliant mesh.

Compliant meshes. Ridges of a smooth surface are points with prescribed differential properties, and reporting them from a mesh inscribed in the surface requires delicate hypothesis on the geometry of that mesh so as to get a certified result. In this paragraph, we assume the mesh provided complies with a number of hypothesis, which guarantee the topology of the ridges reported matches that of the ridges on the smooth surface. To summarize things, a compliant mesh is a mesh dense enough so that (i) umbilics are properly isolated (ii) ridges running next to one another are also properly separated. See [CP05c] for a detailed discussion of compliant meshes.

As 0-level set of the extremality coefficients $$b₀ and $$b₃, ridges are extracted by a marching triangles algorithm ²

As the signs of these extremality coefficients depend on the orientation of the principal directions, we expect both extremalities and vectors orienting the principal direction to be given at each point vertex of the mesh. Except in the neighborhood of umbilics, if the mesh is dense enough, a coherent orientation of principal directions at both endpoints of an edge is chosen such than the angle between the two vectors is acute. This rule, the acute rule, is precisely analyzed in [CP05c]. Moreover, we only seek ridges in triangles such than one can find an orientation of its three vertices such that the three edges are coherently oriented by the acute rule. Such triangles are termed regular. This said, two remarks are in order.

- Regular triangles and ridge segments. A regular triangle has 0 or 2 edges crossed by a max (resp. min) ridge, which is tantamount to a sign change of $$b₀ (resp. $$b₃) along the corresponding edges. In the later case, we say that the triangle contains a ridge segment. Two methods are provided to compute its type, be it elliptic or hyperbolic. First, if fourth order differential quantities are provided, one can use the $$P₁ ($$P₂) values of Eq. (41.2) ((41.2)) for a max (min) ridge. The value of $$P_i for a ridge segment is defined as the mean value of the $$P_i values of the two crossing points on edges (while the value at a crossing point on an edge is the $$b_i-weighted mean value of the values at endpoints). Alternatively, if third order differential quantities only are available, one may use the geometric method developed in [CP05c].

Using the notion of ridge segment, a Ridge_line is defined as a maximal connected sequence of ridge segments of the same type and connected together. Notice that the topology of a Ridge_line is either that of an interval or a circle.

- Non-regular triangles. In the neighborhood of umbilics, triangle are less likely to be regular and the detection of ridges cannot be relevant by this method. This is why we propose another method to detect umbilics independently.

Non compliant meshes : filtering ridges on strength and sharpness. For real world applications dealing with coarse meshes, or meshes featuring degenerate regions or sharp features, or meshes conveying some amount of noise, the compliance hypothesis [CP05c] cannot be met. In that case, it still makes sense to seek loci of points featuring extremality of estimated principal curvatures, but the ridges reported may require filtering. For example, if the principal curvatures are constant - which is the case on a plane or a cylinder, then all points are ridge points. In this context, an appealing notion is that of sharp ridge or prominent ridge. Since ridges are witnessed by zero crossings of $$b₀ and $$b₃, one can expect erroneous detections as long as these coefficients remain small. In order to select the most prominent ridge points, we focus on points where the variation of the curvature is fast along the curvature line. One can observe that, at a ridge point, according to equation 41.2, the second derivative of $$k₁ along its curvature line satisfies $$k₁^''(0) = P₁/(k₁-k₂). Using this observation, one can define the sharpness of a ridge as the integral of the absolute value of $$P₁/(k₁-k₂) along the ridge. As the second derivative of the curvature is homogeneous to the inverse of the cube of a length, the sharpness is homogeneous to the inverse of the square of a length. Multiplying the sharpness by the square of the model size gives a threshold and an associated sharpness-filter which are scale independent. Another filtering is also available with the strength which is the integral of the curvature along the ridge line [OBS04].

41.3   Approximating Umbilics on Triangulated Surface Meshes

The method aims at identifying some vertices of a mesh as umbilics. It assumes principal curvatures and directions are given at each vertex on the mesh.

Algorithm. Assume each vertex $$v of the mesh comes with a patch (a topological disk) around it. Checking whether vertex $$v is an umbilic is a two stages process, which are respectively concerned with the variation of the function $$k₁-k₂ over the patch, and the index of the vertex computed from the boundary of the patch. More precisely, vertex $$v is declared to be an umbilic if the following two conditions are met:

• the function $$k₁-k₂ has its minimum at $$v amongst all the vertices of the patch;
• the deviation $$ of any principal direction along the patch boundary, traversed counter-clock-wise (CCW), has prescribed properties:
  • $$ ]/2,3/2[, then the umbilic is called elliptic,
  • $$ ]-3/2,-/2[, then the umbilic is called a hyperbolic,
  • otherwise the umbilic is called non-generic.

Finding patches around vertices. Given a vertex $$v and a parameter $$t, we aim at defining a collection of triangles around $$v so that (i) this collection defines a topological disk on the triangulation $$T and (ii) its size depends on $$t. First we collect the 1-ring triangles. We define the size $$s of this 1-ring patch as the (Euclidean) distance from $$v to its farthest 1-ring vertex neighbor. We then collect recursively adjacent triangles so that the patch remains a topological disk and such that these triangles are at distance less than $$s × t. Parameter $$t is the only parameter of the algorithm.

Umbilical patches versus ridges. On a generic surface, generic umbilics are traversed by one or three ridges. For compliant meshes, an umbilic can thus be connected to the ridge points located on the boundary of its patch. This functionality is not provided, and the interested reader is referred to [CP05c] for more details.

41.4   Software Design

All classes of this package are templated by the parameter TriangulatedSurfaceMesh, which defines the type of the mesh to which the approximation algorithms operate.

The differential quantities are provided at vertices of this mesh via property maps, a concept commonly used in the Boost library. Scalar data (curvatures and their derivatives) are provided via Vertex2FTPropertyMap concepts, while $$3d vector data (principal directions of curvatures) are provided via Vertex2VectorPropertyMap concepts. The rationale for introducing these concepts is that properties are used independently from the way they are stored. This enables the user to store them internally in extended vertices or externally with maps. We provide a class Vertex2Data_Property_Map_with_std_map to adapt std::map with a boost::associative_property_map to model these concepts.

Output of ridges or umbilics are provided via output iterator.

Approximation of ridges and umbilics are performed by two independent classes, which we now further describe.

41.4.1   Ridge Approximation

The main class is Ridge_approximation<TriangulatedSurfaceMesh,Vertex2FTPropertyMap,Vertex2VectorPropertyMap>. Its construction requires the mesh and the property maps defining the differential quantities for principal curvatures $$k₁ and $$k₂, the third order extremalities $$b₀ and $$b₃, the principal directions of curvature $$d₁ and $$d₂, and the fourth order quantities $$P₁ and $$P₂ if the tagging of ridges as elliptic or hyperbolic is to be done using the polynomials $$P₁ and $$P₂.

Three functions (provided as members and also as global functions) enable the computation of different types of ridges :

• compute_max_ridges (resp. compute_min_ridges) outputs ridges of types MAX_ELLIPTIC_RIDGE and MAX_HYPERBOLIC_RIDGE (resp. MIN_ELLIPTIC_RIDGE and MIN_HYPERBOLIC_RIDGE).
• compute_crest_ridges outputs ridges of types MAX_CREST_RIDGE and MIN_CREST_RIDGE.

These functions allows the user to specify how the elliptic/hyperbolic tagging is carried out. Notice the rationale for the choice of these three functions is simple: each computation needs a single pass over the triangles of the mesh. This should be clear for the min and max ridges. For crests, just notice max and min crests cannot intersect over a triangle.

The ridge lines are stored in Ridge_line objects and output through an iterator. Each ridge line is represented as a list of halfedges of the mesh it crosses with a scalar defining the barycentric coordinate of the crossing point with respect to the half-egde endpoints. Each ridge line comes with its type Ridge_type, its strength and sharpness.

If one chooses to use only third order quantities, the quantities $$P_i do not have to be defined. Then the sharpness will not be defined.

41.4.2   Umbilic Approximation

Umbilics are stored in Umbilic objects, they come with their type : ELLIPTIC_UMBILIC, HYPERBOLIC_UMBILIC or NON_GENERIC_UMBILIC; the vertex of the mesh they are associated to and the list of half-edges representing the contour of the neighborhood.

41.4.3   Models for the Property Map Concepts

41.5   Illustrations - Examples

Example program. ³ The following program computes ridges and umbilics from an off file. It uses the Jet fitting package to estimate the differential quantities. The default output file gives rough data for visualization purpose, a verbose output file may also be asked for. Parameters are

• $$d, the degree of the jet for the Monge_via_jet_fitting class, $$d must be greater or equal to 3 which is the default value;
• $$m, the degree of the Monge representation for the Monge_via_jet_fitting class, $$m must be 3 (the default value) or 4 and smaller than $$d;
• $$a, the number of rings of neighbors collected for the Monge_via_jet_fitting class, in addition the number of vertices collected must be greater than $$Nd:=(d+1)(d+2)/2 to make the approximation possible. $$a may be an integer greater than 1, the value 0 (which is the default) means that the minimum number of rings is collected to make the approximation possible. (Alternatively option $$p allows the specification of a constant number of neighbors);
• $$t, the CGAL::Ridge_order for the distinction between elliptic and hyperbolic ridges, $$t is 3 (default) or 4;
• $$u, the parameter for umbilic patch size, $$u 1 (default is 1).

 
#include <CGAL/Ridges.h> 
#include <CGAL/Umbilics.h>

//this is an enriched Polyhedron with facets' normal
#include "PolyhedralSurf.h"

typedef PolyhedralSurf::Traits          Kernel;
typedef Kernel::FT                      FT;
typedef Kernel::Point_3                 Point_3;
typedef Kernel::Vector_3                Vector_3;

typedef PolyhedralSurf::Vertex_const_handle   Vertex_const_handle;
typedef PolyhedralSurf::Vertex_const_iterator Vertex_const_iterator;
    
typedef CGAL::Vertex2Data_Property_Map_with_std_map<PolyhedralSurf> Vertex2Data_Property_Map_with_std_map;
typedef Vertex2Data_Property_Map_with_std_map::Vertex2FT_map Vertex2FT_map;
typedef Vertex2Data_Property_Map_with_std_map::Vertex2Vector_map Vertex2Vector_map;
typedef Vertex2Data_Property_Map_with_std_map::Vertex2FT_property_map Vertex2FT_property_map;
typedef Vertex2Data_Property_Map_with_std_map::Vertex2Vector_property_map Vertex2Vector_property_map;

//RIDGES
typedef CGAL::Ridge_line<PolyhedralSurf> Ridge_line;
typedef CGAL::Ridge_approximation < PolyhedralSurf,
				    back_insert_iterator< std::vector<Ridge_line*> >,
				    Vertex2FT_property_map,
				    Vertex2Vector_property_map > Ridge_approximation;
//UMBILICS
typedef CGAL::Umbilic<PolyhedralSurf> Umbilic;
typedef CGAL::Umbilic_approximation < PolyhedralSurf,
				      back_insert_iterator< std::vector<Umbilic*> >, 
				      Vertex2FT_property_map, 
				      Vertex2Vector_property_map > Umbilic_approximation;

//create property maps
Vertex2FT_map vertex2k1_map, vertex2k2_map, 
  vertex2b0_map, vertex2b3_map, 
  vertex2P1_map, vertex2P2_map;
Vertex2Vector_map vertex2d1_map, vertex2d2_map;

Vertex2FT_property_map vertex2k1_pm(vertex2k1_map), vertex2k2_pm(vertex2k2_map), 
  vertex2b0_pm(vertex2b0_map), vertex2b3_pm(vertex2b3_map), 
  vertex2P1_pm(vertex2P1_map), vertex2P2_pm(vertex2P2_map);
Vertex2Vector_property_map vertex2d1_pm(vertex2d1_map), vertex2d2_pm(vertex2d2_map);

int main(int argc, char *argv[])
{  
  //compute differential quantities with the jet fitting package
	...
  //initialize the property maps
	...
   //---------------------------------------------------------------------------
  //Ridges
  //--------------------------------------------------------------------------
  Ridge_approximation ridge_approximation(P, 
					  vertex2k1_pm, vertex2k2_pm,
					  vertex2b0_pm, vertex2b3_pm,
					  vertex2P1_pm, vertex2P2_pm,
					  vertex2d1_pm, vertex2d2_pm);
  std::vector<Ridge_line*> ridge_lines;
  back_insert_iterator<std::vector<Ridge_line*> > ii(ridge_lines);
  
  //Find MAX_RIDGE, MIN_RIDGE, CREST or all ridges
  ridge_approximation.compute_max_ridges(ii, tag_order);  
  ridge_approximation.compute_min_ridges(ii, tag_order);  
  ridge_approximation.compute_crest_ridges(ii, tag_order);  

  //---------------------------------------------------------------------------
  // UMBILICS
  //--------------------------------------------------------------------------
  Umbilic_approximation umbilic_approximation(P, 
					      vertex2k1_pm, vertex2k2_pm,
					      vertex2d1_pm, vertex2d2_pm);
  std::vector<Umbilic*> umbilics;
  back_insert_iterator<std::vector<Umbilic*> > umb_it(umbilics);
  umbilic_approximation.compute(umb_it, umb_size);
 }

Example: Fig. 41.5. For that figure, the data have been computed as follows

./Compute_Ridges_Umbilics -f data/ellipsoid_u_0.02.off -d 4 -m 4 -a 3 -t 3

nb of umbilics 4
Umbilic at location (-0.80899 0.00426003 0.293896) of type elliptic
Umbilic at location (-0.811197 0.0122098 -0.292259) of type elliptic
Umbilic at location (0.808372 -0.00551307 -0.29431) of type elliptic
Umbilic at location (0.81413 0.0018689 0.290339) of type elliptic

Example: Fig. 41.6. This Fig. illustrates the filtering of crest ridges on a mechanical model, and has been computed as follows:

./Compute_Ridges_Umbilics -f data/mecanic.off -d 4 -m 4 -a 4 -t 4

Footnotes