added test and revised doc

.gitattributes

@@ -1696,6 +1696,8 @@ Ridges_3/doc_tex/Ridges_3/david_eye_crest.eps -text svneol=unset#application/pos
 Ridges_3/doc_tex/Ridges_3/david_eye_crest.png -text svneol=unset#image/png
 Ridges_3/doc_tex/Ridges_3/ellipsoid_ridges.eps -text svneol=unset#application/postscript
 Ridges_3/doc_tex/Ridges_3/ellipsoid_ridges.png -text svneol=unset#image/png
+Ridges_3/doc_tex/Ridges_3/index_umbilic.eps -text svneol=unset#application/postscript
+Ridges_3/doc_tex/Ridges_3/index_umbilic.png -text svneol=unset#image/png
 Ridges_3/doc_tex/Ridges_3/lemon.eps -text svneol=unset#application/postscript
 Ridges_3/doc_tex/Ridges_3/lemon.png -text svneol=unset#image/png
 Ridges_3/doc_tex/Ridges_3/mecanic-sub1_crest-jpg.eps -text svneol=unset#application/postscript

Ridges_3/doc_tex/Ridges_3/Ridges_3_user.tex

@@ -1,9 +1,4 @@
 %note on pictures
-
-% ellipsoid_ridges 
-% ./blind -f data/ellipsoid_u_0.02.off -d4 -m4 -a3 -t3 
-% ../../demo/Ridges_3/introspect-qt data/ellipsoid_u_0.02.off data/data_ellipsoid_u_0.02.offRIDGES-d4-m4-t3-a3-p0.4ogl.txt 0 0
-
 %david pict for old algo, don't know parameters, todo... and get better pict
 
 %\newcommand{\makeremark}[2]{{
@@ -64,12 +59,10 @@ functions of the package are provided in section \ref{examples}.
 \label{smooth}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-\FC{I did not find cp-ssulc-05 in cgal\_manual.bib???????}
-
 For a detailed introduction to ridges and related topics, the reader
 may consult 
 \cite{cgal:hgygm-ttdpf-99,cgal:p-gd-01}, as well as
-the following survey article \cite{cgal:cp-dtges-05}.
+the following survey article \cite{cgal:cp-ssulc-05}.
 %%
 In the sequel, we just introduce the basic notions so as to explain
 our algorithms.  Consider a smooth embedded surface, and denote $k_1$
@@ -189,21 +182,18 @@ Hence ridge types are characterized by
 \ccc{RED_CREST_RIDGE} if  $b_3=0$ and $P_2<0$ and $|k_2|>|k_1|$
 \end{itemize}
 
-\FC{picture for the index}
-\FC{I removed  $k_1=k_2$ as umbilics have been defined several times before}
 
-
-As illustrated on Fig. \ref{umbilics} and \ref{}, the patterns made by
-curvature lines around an umbilic can be characterized using the
-notion of {\em index} associated to the principal directions ---see
-also \cite{cgal:cp-dtges-05}.
+As illustrated on Fig. \ref{umbilics} and \ref{index_umbilic}, the
+patterns made by curvature lines around an umbilic can be
+characterized using the notion of {\em index} associated to the
+principal directions ---see also \cite{cgal:cp-ssulc-05}.
 %%
-As depicted on Fig. \ref{}, consider a small circuit $C$ around the
+As depicted on Fig. \ref{index_umbilic}, consider a small circuit $C$ around the
 umbilic, and a point $p \in C$. Starting from an initial orientation
 $u$ of a tangent vector to the curvature line through point $p$,
 propagate {\em by continuity} this orientation around the circuit.  The
 index is defined by the angle swept by $u$ around this revolution,
-normalized by $2\pi$. On our example, the index is thus ???????.
+normalized by $2\pi$. 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
@@ -211,7 +201,20 @@ then it is called a
  Otherwise the umbilic is qualified
 \ccc{UMBILIC_NON_GENERIC}.
 
-\FC{find the right figure: asymptotes are missing}
+
+\begin{figure}[!ht]
+\begin{ccTexOnly}
+\centerline{
+\includegraphics[width=.5\linewidth]{Ridges_3/index_umbilic}}
+\end{ccTexOnly}
+\caption{Index $1/2$ umbilic.}
+\label{index_umbilic}
+\begin{ccHtmlOnly}
+<CENTER> <img border=0 src="./index_umbilic.png" width=200>
+</CENTER>
+\end{ccHtmlOnly}
+\end{figure}
+
 
 \begin{figure}[!ht]
 \begin{ccTexOnly}
@@ -234,8 +237,6 @@ then it is called a
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \label{ridge-mesh}
 
-\FC{marching line -> marching triangles???}
-
 \paragraph{Compliant meshes.}
 Ridges of a smooth surface are points with prescribed differential
 properties, and reporting them from a mesh requires delicate
@@ -253,25 +254,20 @@ that a one-manifold is reported.}.
 %%
 %%on the surface.  
 %%
-As the signs of these extremality coefficients depends on the
+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 {\em acute rule}, is
-precisely analyzed in \cite{cgal:tr-cp-tdare-05}.
+precisely analyzed in \cite{cgal:cp-tdare-05}.
 %%
 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
 {\em regular}. This said, two remarks are in order.
 
-\FC{above: why do we expect the extremalities to be given since we use
-the acute rule???????}
-
-\FC{I introduced the def of a ridge segment type below}
-
 \noindent {\em ---Regular triangles and ridge segments.}
 A regular triangle has 0 or 2 edges crossed by a ridge of a given
 color (meaning that $b_0$ or $b_3$ change sign on the edge). In the
@@ -279,9 +275,12 @@ later case, we say that the triangle contains a ridge segment.
 %%
 Two methods are provided to compute its type, as elliptic or
 hyperbolic. First, if fourth order differential quantities are
-provided, one can use the $P_1$ ($P_2$) polynomial of Eq.
-(\ref{{eq:taylor_along_line}}) ((\ref{eq:taylor_along_red_line})) for
-a blue (red) ridge.
+provided, one can use the $P_1$ ($P_2$) values of Eq.
+(\ref{eq:taylor_along_line}) ((\ref{eq:taylor_along_red_line})) for
+a blue (red) 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
@@ -290,35 +289,19 @@ available, one may use the geometric method developed in
 
 Using the notion of ridge segment, a \ccc{Ridge_line} is defined as a
 maximal sequence of ridge segments of the same type and connected
-together.  Notice the the topology of a \ccc{Ridge_line} is either
+together.  Notice that the topology of a \ccc{Ridge_line} is either
 that of an interval or a circle.
 
-%We provide two methods to distinguish elliptic and hyperbolic
-%ridges. The first one uses the signs of fourth order quantities
-%$P_i$. The other uses only third order quantities, details can be
-%found in \cite{cgal:cp-tdare-05}.
-
-
-\FC{above: if fourth order taging used, type = majority of the tags of
-the triangle vertices???????}
-
 \noindent {\em ---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.
 
-
-%We provide two methods to distinguish elliptic and hyperbolic
-%ridges. The first one uses the signs of fourth order quantities
-%$P_i$. The other uses only third order quantities, details can be
-%found in \cite{cgal:cp-tdare-05}.
-
-
 \paragraph{Non compliant meshes.}
 %%
 For real world applications dealing with coarse meshes, or meshes
 featuring degenerate regions or sharp features, meshes conveying some
-amount of noise, the hypothesis \cite{cgal:tr-cp-tdare-05} cannot be
+amount of noise, the hypothesis \cite{cgal:cp-tdare-05} 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
@@ -341,7 +324,7 @@ the model size gives a threshold and an associated sharpness-filter
 which are scale independent. Another filtering is also available with
 the {\em strength } which is the integral of the curvature along the
 ridge line
-\cite{cgal:yo-ab-hps-04}.
+\cite{cgal:ybs-rvlmi-04}.
 
 
 
@@ -370,7 +353,7 @@ declared to be an umbilic if the following two conditions are met:
 the function $k_1-k_2$ has its minimum at $v$ amongst all the
 vertices of the patch;
 \item
-the deviation $\delta$ of any principal direction along the the patch
+the deviation $\delta$ of any principal direction along the patch
 boundary, traversed counter-clock-wise (CCW), has prescribed
 properties:
 %%
@@ -384,7 +367,6 @@ $\delta \geq 3\pi/2$ or $\delta \leq -3\pi/2$ then the umbilic is called non-gen
 \end{itemize}
 \end{itemize}
 
-\FC{of any  ppal direction I guess...}
 
 \paragraph{Finding patches around vertices.}
 Given a vertex $v$, we aim at defining a collection of triangles
@@ -452,7 +434,7 @@ Two functions enable to either compute all ridges
 These functions allows the user to specify how the elliptic/hyperbolic
 tagging is carried out.  For the second function, one may further
 specify a ridge intersection type, that is an entry in the following
-enumeration \ccc{ BLUE_RIDGE, RED_RIDGE, CREST}. 
+enumeration \ccc{ BLUE_RIDGE, RED_RIDGE, CREST_RIDGE}. 
 %%
 Notice the rationale for allowing to specify an entry in this enum is
 simple: all ridges of such a type can be computed by a single pass
@@ -461,10 +443,6 @@ blue ridges. For crests, just notice red and blue crests cannot
 intersect over a triangle.
 \medskip
 
-\FC{in the interrogation type, why do we have crest .. while it is not
-in the enum gathering all types????????}
-
-
 The ridge lines are stored in
 \ccc{Ridge_line} objects and output through an iterator. 
 Each ridge line is represented as a list of half-edges of the mesh it
@@ -504,23 +482,31 @@ using \ccc{std::maps} and \ccc{boost::associative_property_map}.
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \label{examples}
 
-\FC{provide ranges and default values}
-
 The example program computes ridges and umbilics from an off file. It
 uses the Jet fitting package to estimate the differential quantities.
 The output file contains data to be visualized with the demo program introspect-qt.
 Parameters are 
 \begin{itemize}
 \item
-d, the degree of the jet for the \ccc{Monge_via_jet_fitting} class;
+d, the degree of the jet for the \ccc{Monge_via_jet_fitting} class, d
+must be greater or equal to 3 which is the default value;
 \item
-m, the degree of the Monge representation for the \ccc{Monge_via_jet_fitting} class;
+m, the degree of the Monge representation for the
+\ccc{Monge_via_jet_fitting} class, m must be 3 (the default value) or
+4 and smaller than d;
 \item
-a, the number of rings of neighbors collected for the \ccc{Monge_via_jet_fitting} class;
+a, the number of rings of neighbors collected for the
+\ccc{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. (Alternativelly option p
+enable to ask for a constant number of neighbors);
 \item
-t, the \ccc{Tag_order} for the distinction between elliptic and hyperbolic ridges;
+t, the \ccc{Tag_order} for the distinction between elliptic and
+hyperbolic ridges, t is 3 (default) or 4;
 \item
-u, the parameter for umbilic patch size.
+u, the parameter for umbilic patch size, $u \geq 1$ (default is 1).
 \end{itemize}
 
 \begin{ccExampleCode} 

Ridges_3/doc_tex/Ridges_3_ref/TriangularPolyhedralSurface.tex

@@ -18,8 +18,8 @@
 
 \ccDefinition
   
-The concept \ccRefName\ describes a polyhedral surface, based on an
-edge data structure, whose faces are all triangular.
+The concept \ccRefName\ describes a polyhedral surface, based on a
+halfedge data structure, whose faces are all triangular.
 
 
 %\ccTypes

Ridges_3/examples/Ridges_3/README

@@ -35,6 +35,7 @@ EXAMPLES
 ./blind -f data/ellipsoid_u_0.02.off -d4 -m4 -a3 -t3 -v
 
 visu with:
+
  ../../demo/Ridges_3/introspect-qt data/ellipsoid_u_0.02.off data/data_ellipsoid_u_0.02.offRIDGES-d4-m4-t4-a3-p0.4ogl.txt 0 0
 
 

Ridges_3/examples/Ridges_3/blind.cpp

@@ -309,6 +309,7 @@ int main(int argc, char *argv[])
     for (iter_lines = ridge_lines.begin();iter_lines!=iter_end;iter_lines++) 
       out_verb << **iter_lines; 
 
+  std::cout << "Compute umbilics..." << std::endl;
   //---------------------------------------------------------------------------
   // UMBILICS
   //--------------------------------------------------------------------------

Ridges_3/include/CGAL/PolyhedralSurf_neighbors.h

@@ -188,7 +188,23 @@ compute_neighbors(Vertex_handle v,
   FT d_current = firstGate.d();
   // main loop
   while ( !GatePQ.empty() && d_current <= d_max ) {
-    Gate gate = GatePQ.top();
+ 
+  //debug check if the contour is closed  	 
+        typename std::list<Halfedge_handle>::iterator itbc = contour.begin(), 	 
+                                           itec = contour.end(), 	 
+        h_cur, h_next; 	 
+      for (; itbc != itec; itbc++) 	 
+        { 	 
+          h_cur = itbc; 	 
+          if ( h_cur != (--contour.end()) ) {h_next = ++h_cur; h_cur--;} 	 
+          else h_next = contour.begin(); 	 
+          assert( (*h_cur)->vertex() == (*h_next)->opposite()->vertex() ); 	 
+        //cout << endl << &**itbc ; 	 
+      } 	 
+  //debug  	 
+    //cout << endl; cout << endl;
+
+   Gate gate = GatePQ.top();
     GatePQ.pop();
     d_current = gate.d();
     Halfedge_handle he = gate.he(), he1, he2;

Ridges_3/test/Ridges_3/PolyhedralSurf.cpp

@@ -0,0 +1,64 @@
+#include "PolyhedralSurf.h"
+
+// Vector_3 PolyhedralSurf::getHalfedge_vector(Halfedge * h)
+// {
+//   Vector_3 v = h->opposite()->vertex()->point() - h->vertex()->point();
+//   return v;
+// }
+
+// double PolyhedralSurf::
+// compute_mean_edges_length_around_vertex(Vertex* v)
+// {
+//   Halfedge_around_vertex_circulator
+//     hedgeb = v->vertex_begin(), hedgee = hedgeb;
+//   int  count_he = 0;
+//   double sum = 0.;
+//   CGAL_For_all(hedgeb, hedgee)
+//     {
+//       sum += hedgeb->getLength();
+//       count_he++;
+//     }
+//   return sum/count_he;
+// }
+
+// void PolyhedralSurf::compute_edges_length()
+// {
+//   std::for_each(this->halfedges_begin(), this->halfedges_end(),
+// 		Edge_length());
+// }
+
+
+void PolyhedralSurf::compute_facets_normals()
+{
+  std::for_each(this->facets_begin(), this->facets_end(),
+		Facet_unit_normal()); 
+}
+
+Vector_3 PolyhedralSurf::computeFacetsAverageUnitNormal(Vertex * v)
+{
+  Halfedge *h;
+  Facet *f;
+  Vector_3 sum(0., 0., 0.), n;
+
+  Halfedge_around_vertex_circulator
+    hedgeb = v->vertex_begin(), hedgee = hedgeb;
+
+  do
+    {
+      h = &(*hedgeb);
+      if (h->is_border_edge())
+	{
+	  hedgeb++;
+	  continue;
+	}
+
+      f = &(*h->facet());
+      n = f->getUnitNormal();
+      sum = (sum + n);
+      hedgeb++;
+    }
+  while (hedgeb != hedgee);
+  sum = sum / std::sqrt(sum * sum);
+  return sum;
+}
+

Ridges_3/test/Ridges_3/PolyhedralSurf.h

@@ -0,0 +1,82 @@
+#ifndef _POLYHEDRALSURF_H_
+#define _POLYHEDRALSURF_H_
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <algorithm>
+#include <vector>
+#include <list>
+
+#include <CGAL/Cartesian.h>
+#include <CGAL/Polyhedron_3.h>
+#include <CGAL/IO/Polyhedron_iostream.h>
+
+#include "PolyhedralSurf_operations.h" 
+
+//----------------------------------------------------------------
+// A redefined items class for the Polyhedron_3 with 
+// a refined facet with a normal vector
+//---------------------------------------------------------------
+
+//----------------------------------------------------------------
+// Facet with normal and possibly more types. types are recovered
+//from the FGeomTraits template arg
+//----------------------------------------------------------------
+template < class Refs, class Tag, class FGeomTraits >
+class My_facet:public CGAL::HalfedgeDS_face_base < Refs, Tag >
+{
+public:
+ typedef typename FGeomTraits::Vector_3 Vector_3;
+
+protected:
+  Vector_3 normal;
+public:
+  My_facet() {}
+  Vector_3 & getUnitNormal() { return normal; }
+  void setNormal(Vector_3 n) { normal = n; }
+};
+
+//------------------------------------------------
+// Wrappers [Vertex, Face, Halfedge]
+//------------------------------------------------
+struct Wrappers_VFH:public CGAL::Polyhedron_items_3 {
+  // wrap face
+  //NOTE: [HDS, Face] renamed [Polyhedron, Facet]
+  template < class Refs, class Traits > struct Face_wrapper {
+    //typedef typename Traits::Vector_3 Vector_3;
+    //all types needed by the facet...
+    typedef struct {
+    public:
+       typedef typename Traits::Vector_3 Vector_3;
+     } FGeomTraits;
+    //custom type instantiated...
+    typedef My_facet < Refs, CGAL::Tag_true, FGeomTraits > Face;
+  };
+};
+
+//------------------------------------------------
+//PolyhedralSurf with facet normal operations
+//------------------------------------------------
+typedef double                FT;
+typedef CGAL::Cartesian<FT>  Kernel;
+typedef CGAL::Polyhedron_3 < Kernel, Wrappers_VFH > Polyhedron;
+typedef Kernel::Vector_3 Vector_3;
+
+class PolyhedralSurf:public Polyhedron {
+public:
+  PolyhedralSurf() {}
+  
+  //static Vector_3 getHalfedge_vector(Halfedge * h);
+  //double compute_mean_edges_length_around_vertex(Vertex * v);
+  //void compute_edges_length();
+
+  void compute_facets_normals();
+  Vector_3 computeFacetsAverageUnitNormal(Vertex * v);
+};
+
+#endif
+
+
+
+
+

Ridges_3/test/Ridges_3/PolyhedralSurf_operations.h

@@ -0,0 +1,31 @@
+#ifndef _POLY_OP_H_
+#define _POLY_OP_H_
+
+/* struct Edge_length { */
+/*   template < class HalfEdge >  */
+/*   void operator() (HalfEdge & h)  */
+/*     { */
+/*       double d = */
+/* 	CGAL::squared_distance(h.prev()->vertex()->point(), */
+/* 			 h.vertex()->point()); */
+/*       h.setLength(CGAL::sqrt(d)); */
+/*   } */
+/* }; */
+
+//the facet stores the normal
+struct Facet_unit_normal {
+  template < class Facet > 
+  void operator() (Facet & f) 
+    {
+      typename Facet::Halfedge_handle h = f.halfedge();
+      typename Facet::Vector_3 normal =
+	CGAL::cross_product(h->vertex()->point() - 
+			    h->opposite()->vertex()->point(),
+			    h->next()->vertex()->point() -
+			    h->opposite()->vertex()->point());
+      f.setNormal( normal / CGAL::sqrt(normal * normal));
+    }
+};
+
+
+#endif

Ridges_3/test/Ridges_3/PolyhedralSurf_rings.h

@@ -0,0 +1,162 @@
+#ifndef _PSURF_RINGS_H_
+#define _PSURF_RINGS_H_
+
+using namespace std;
+
+//---------------------------------------------------------------------------
+//T_PolyhedralSurf_rings
+//---------------------------------------------------------------------------
+template < class TPoly > class T_PolyhedralSurf_rings
+{
+protected:
+  //Polyhedron
+  typedef typename TPoly::Vertex Vertex;
+  typedef typename TPoly::Halfedge Halfedge;
+  typedef typename TPoly::Facet Facet;
+  typedef typename TPoly::Halfedge_around_vertex_circulator  Halfedge_around_vertex_circulator;
+  typedef typename TPoly::Vertex_iterator Vertex_iterator;
+  typedef std::map<Vertex*, int> Vertex2int_map_type;
+  Vertex2int_map_type ring_index_map;
+  
+  //vertex indices are initialised to -1
+  void reset_ring_indices(std::vector < Vertex * >&vces);
+
+  //i >= 1; from a start vertex on the current i-1 ring, push non-visited neighbors
+  //of start in the nextRing and set indices to i. Also add these vertices in all.
+  void push_neighbours_of(Vertex * start, int ith,
+			  std::vector < Vertex * >&nextRing,
+			  std::vector < Vertex * >&all);
+
+  //i >= 1, from a currentRing i-1, collect all neighbors, set indices
+  //to i and store them in nextRing and all.
+  void collect_ith_ring(int ith,
+			std::vector < Vertex * >&currentRing,
+			std::vector < Vertex * >&nextRing,
+			std::vector < Vertex * >&all);
+ 
+ public:  
+  T_PolyhedralSurf_rings(TPoly& P);
+  
+  //collect i>=1 rings : all neighbours up to the ith ring,
+  void collect_i_rings(Vertex* v, 
+		       int ring_i, 
+		       std::vector < Vertex * >& all);
+
+  //collect enough rings (at least 1), to get at least min_nb of neighbors
+  void collect_enough_rings(Vertex* v,
+			    unsigned int min_nb,
+			    std::vector < Vertex * >& all);
+};
+
+////IMPLEMENTATION/////////////////////////////////////////////////////////////////////
+
+template < class TPoly >
+T_PolyhedralSurf_rings <TPoly>::
+T_PolyhedralSurf_rings(TPoly& P)
+{
+  //init the ring_index_map
+  Vertex_iterator itb = P.vertices_begin(), ite = P.vertices_end();
+  for(;itb!=ite;itb++) ring_index_map[&*itb] = -1; 
+}
+
+template < class TPoly >
+void T_PolyhedralSurf_rings <TPoly>::
+push_neighbours_of(Vertex * start, int ith,
+		   std::vector < Vertex * >&nextRing,
+		   std::vector < Vertex * >&all)
+{
+  Vertex *v;  
+  Halfedge_around_vertex_circulator
+    hedgeb = start->vertex_begin(), hedgee = hedgeb;
+ 
+ CGAL_For_all(hedgeb, hedgee)
+  {
+    v = &*(hedgeb->opposite()->vertex());
+    if (ring_index_map[v] != -1)  continue;//if visited: next
+    
+    ring_index_map[v] = ith;
+    nextRing.push_back(v);
+    all.push_back(v);
+  }
+}
+
+template <class TPoly>
+void T_PolyhedralSurf_rings <TPoly>::
+collect_ith_ring(int ith, std::vector < Vertex * >&currentRing,
+		 std::vector < Vertex * >&nextRing,
+		 std::vector < Vertex * >&all)
+{
+  typename std::vector < Vertex * >::iterator 
+    itb =    currentRing.begin(), ite = currentRing.end();
+  
+  CGAL_For_all(itb, ite) push_neighbours_of(*itb, ith, nextRing, all);
+}
+
+template <class TPoly>
+  void T_PolyhedralSurf_rings <TPoly>::
+reset_ring_indices(std::vector < Vertex * >&vces)
+{
+  typename std::vector < Vertex * >::iterator 
+    itb = vces.begin(), ite = vces.end();
+  CGAL_For_all(itb, ite)  ring_index_map[*itb] = -1;
+}
+ 
+template <class TPoly>
+  void T_PolyhedralSurf_rings <TPoly>::
+collect_i_rings(Vertex* v, 
+		int ring_i, 
+		std::vector < Vertex * >& all)
+{
+  std::vector<Vertex*> current_ring, next_ring; 
+  std::vector<Vertex*> *p_current_ring, *p_next_ring;
+  assert(ring_i >= 1);
+  //initialize
+  p_current_ring = &current_ring;
+  p_next_ring = &next_ring;
+  ring_index_map[v] = 0;
+  current_ring.push_back(v);
+  all.push_back(v);
+
+  for (int i=1; i<=ring_i; i++)
+    {
+      collect_ith_ring(i, *p_current_ring, *p_next_ring, all);
+      //next round must be launched from p_nextRing...
+      p_current_ring->clear();
+      std::swap(p_current_ring, p_next_ring);
+    }
+  //clean up
+  reset_ring_indices(all);
+}
+
+template <class TPoly>
+  void T_PolyhedralSurf_rings <TPoly>::
+collect_enough_rings(Vertex* v,
+		     unsigned int min_nb,
+		     std::vector < Vertex * >& all)
+{
+  std::vector<Vertex*> current_ring, next_ring; 
+  std::vector<Vertex*> *p_current_ring, *p_next_ring;
+
+  //initialize
+  p_current_ring = &current_ring;
+  p_next_ring = &next_ring;
+  ring_index_map[v] = 0;
+  current_ring.push_back(v);
+  all.push_back(v);
+
+  int i = 1;
+
+  while ( (all.size() < min_nb) &&  (p_current_ring->size() != 0) ) 
+    {
+      collect_ith_ring(i, *p_current_ring, *p_next_ring, all);
+      //next round must be launched from p_nextRing...
+      p_current_ring->clear();
+      std::swap(p_current_ring, p_next_ring);
+      i++;
+    }
+  //clean up
+  reset_ring_indices(all);
+}
+
+
+#endif

Ridges_3/test/Ridges_3/makefile

@@ -0,0 +1,94 @@
+# Created by the script cgal_create_makefile
+# This is the makefile for compiling a CGAL application.
+
+#---------------------------------------------------------------------#
+#                    include platform specific settings
+#---------------------------------------------------------------------#
+# Choose the right include file from the <cgalroot>/make directory.
+include $(CGAL_MAKEFILE)
+
+PROFOPT=  -g #-O3 # 
+CFLAGS=${PROFOPT} 
+#---------------------------------------------------------------------#
+#                    compiler flags
+#---------------------------------------------------------------------#
+LAPACK_INCDIRS = -I$(LAPACK_DIR)/SRC -I$(LAPACK_DIR)
+
+#note line -I../../../Jet_fitting_3/include needed since 
+#jet_fitting_3 not in the release yet
+CXXFLAGS = \
+           -I../../include \
+           -I../../../Jet_fitting_3/include \
+           $(CGAL_CXXFLAGS) \
+           $(LONG_NAME_PROBLEM_CXXFLAGS) \
+	   $(LAPACK_INCDIRS) \
+	   $(CFLAGS)
+
+#---------------------------------------------------------------------#
+#                    linker flags
+#---------------------------------------------------------------------#
+F2CDIR  = $(LAPACK_DIR)/F2CLIBS
+#LAPACK_LDLIBS  = -L$(LAPACK_DIR) -L$(F2CDIR) \
+#	-llapack_LINUX \
+#          -lblas_LINUX -ltmglib_LINUX \
+#          -lF77 -lI77 -lm -lc
+
+#LAPACK_LDLIBS  = 	-l$(LAPACK_DIR)/lapack_LINUX.a \
+#          -l$(LAPACK_DIR)/blas_LINUX.a -l$(LAPACK_DIR)/tmglib_LINUX.a \
+#          $(F2CDIR)/libF77.a $(F2CDIR)/libI77.a -lm -lc
+
+#this is the only way to get through the linking... for me, Marc
+F2CDIR  = $(LAPACK_DIR)/F2CLIBS
+LAPACK_LDLIBS  = $(LAPACK_DIR)/lapack_LINUX.a \
+           $(LAPACK_DIR)/blas_LINUX.a $(LAPACK_DIR)/tmglib_LINUX.a \
+           $(F2CDIR)/libF77.a $(F2CDIR)/libI77.a -lm -lc
+
+LIBPATH = \
+          $(CGAL_LIBPATH)
+
+LDFLAGS = \
+	${LAPACK_LDLIBS} \
+	$(LONG_NAME_PROBLEM_LDFLAGS) \
+        $(CGAL_LDFLAGS) 
+
+#---------------------------------------------------------------------#
+#                    target entries
+#---------------------------------------------------------------------#
+
+all:            \
+                ridge_test 
+
+RIDGE_TEST_OBJS=PolyhedralSurf.o ridge_test.o
+ridge_test: $(RIDGE_TEST_OBJS)
+	$(CGAL_CXX) $(LIBPATH) -o ridge_test $(RIDGE_TEST_OBJS) \
+	$(LDFLAGS) 
+
+clean: \
+        ridge_test.clean PolyhedralSurf.clean 
+
+depend: 
+	makedepend *.[Ch]
+
+#---------------------------------------------------------------------#
+#                    suffix rules
+#---------------------------------------------------------------------#
+
+.cpp$(OBJ_EXT):
+	$(CGAL_CXX) $(CXXFLAGS) $(OBJ_OPT) $<
+
+# DO NOT DELETE
+
+PolyhedralSurf.o: /usr/include/stdlib.h /usr/include/features.h
+PolyhedralSurf.o: /usr/include/sys/cdefs.h /usr/include/gnu/stubs.h
+PolyhedralSurf.o: /usr/include/sys/types.h /usr/include/bits/types.h
+PolyhedralSurf.o: /usr/include/bits/wordsize.h /usr/include/bits/typesizes.h
+PolyhedralSurf.o: /usr/include/time.h /usr/include/endian.h
+PolyhedralSurf.o: /usr/include/bits/endian.h /usr/include/sys/select.h
+PolyhedralSurf.o: /usr/include/bits/select.h /usr/include/bits/sigset.h
+PolyhedralSurf.o: /usr/include/bits/time.h /usr/include/sys/sysmacros.h
+PolyhedralSurf.o: /usr/include/bits/pthreadtypes.h /usr/include/bits/sched.h
+PolyhedralSurf.o: /usr/include/alloca.h /usr/include/stdio.h
+PolyhedralSurf.o: /usr/include/libio.h /usr/include/_G_config.h
+PolyhedralSurf.o: /usr/include/wchar.h /usr/include/bits/wchar.h
+PolyhedralSurf.o: /usr/include/gconv.h /usr/include/bits/stdio_lim.h
+PolyhedralSurf.o: /usr/include/bits/sys_errlist.h PolyhedralSurf_operations.h

Ridges_3/test/Ridges_3/ridge_test.cpp

@@ -0,0 +1,251 @@
+#include <CGAL/basic.h>
+#include <CGAL/Cartesian.h>
+
+#include <cstdio>
+#include <iostream>
+#include <fstream>
+#include <cstdlib>
+#include <vector>
+#include <list>
+
+#include <CGAL/Ridges.h> 
+#include <CGAL/Umbilic.h>
+
+// #include <CGAL/Monge_via_jet_fitting.h> //does not work since not in the release yet
+// #include <CGAL/Lapack/Linear_algebra_lapack.h>
+#include "../../../Jet_fitting_3/include/CGAL/Monge_via_jet_fitting.h" 
+#include "../../../Jet_fitting_3/include/CGAL/Lapack/Linear_algebra_lapack.h" 
+ 
+//this is an enriched Polyhedron with facets' normal
+#include "PolyhedralSurf.h"
+#include "PolyhedralSurf_rings.h"
+
+typedef PolyhedralSurf::Traits          Kernel;
+typedef Kernel::FT                      FT;
+typedef Kernel::Point_3                 Point_3;
+typedef Kernel::Vector_3                Vector_3;
+
+typedef PolyhedralSurf::Vertex          Vertex;
+typedef PolyhedralSurf::Vertex_handle   Vertex_handle;
+typedef PolyhedralSurf::Vertex_iterator Vertex_iterator;
+
+typedef T_PolyhedralSurf_rings<PolyhedralSurf> Poly_rings;
+typedef CGAL::Monge_via_jet_fitting<Kernel>    Monge_via_jet_fitting;
+typedef Monge_via_jet_fitting::Monge_form      Monge_form;
+typedef Monge_via_jet_fitting::Monge_form_condition_numbers Monge_form_condition_numbers;
+      
+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, to be moved in main?
+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);
+
+  // default fct parameter values and global variables
+  unsigned int d_fitting = 4;
+  unsigned int d_monge = 4;
+  unsigned int nb_rings = 0;//seek min # of rings to get the required #pts
+  unsigned int nb_points_to_use = 0;//
+  Ridge_approximation::Tag_order tag_order = Ridge_approximation::Tag_3;
+  double umb_size = 1;
+  bool verbose = false;
+  unsigned int min_nb_points = (d_fitting + 1) * (d_fitting + 2) / 2;
+
+/* gather points around the vertex v using rings on the
+   polyhedralsurf. the collection of points resorts to 3 alternatives:
+   1. the exact number of points to be used
+   2. the exact number of rings to be used
+   3. nothing is specified
+*/
+void gather_fitting_points(Vertex* v, 
+			   std::vector<Point_3> &in_points,
+			   Poly_rings& poly_rings)
+{
+  //container to collect vertices of v on the PolyhedralSurf
+  std::vector<Vertex*> gathered; 
+  //initialize
+  in_points.clear();  
+  
+  //OPTION -p nb_points_to_use, with nb_points_to_use != 0. Collect
+  //enough rings and discard some points of the last collected ring to
+  //get the exact "nb_points_to_use" 
+  if ( nb_points_to_use != 0 ) {
+    poly_rings.collect_enough_rings(v, nb_points_to_use, gathered);//, vpm);
+    if ( gathered.size() > nb_points_to_use ) gathered.resize(nb_points_to_use);
+  }
+  else { // nb_points_to_use=0, this is the default and the option -p is not considered;
+    // then option -a nb_rings is checked. If nb_rings=0, collect
+    // enough rings to get the min_nb_points required for the fitting
+    // else collect the nb_rings required
+    if ( nb_rings == 0 ) 
+      poly_rings.collect_enough_rings(v, min_nb_points, gathered);//, vpm);
+    else poly_rings.collect_i_rings(v, nb_rings, gathered);//, vpm);
+  }
+     
+  //store the gathered points
+  std::vector<Vertex*>::iterator 
+    itb = gathered.begin(), ite = gathered.end();
+  CGAL_For_all(itb,ite) in_points.push_back((*itb)->point());
+}
+
+/* Use the jet_fitting package and the class Poly_rings to compute
+   diff quantities.
+*/
+void compute_differential_quantities(PolyhedralSurf& P, Poly_rings& poly_rings)
+{
+  //container for approximation points
+  std::vector<Point_3> in_points;
+ 
+  //MAIN LOOP
+  Vertex_iterator vitb = P.vertices_begin(), vite = P.vertices_end();
+  for (; vitb != vite; vitb++) {
+    //initialize
+    Vertex* v = &(*vitb);
+    in_points.clear();  
+    Monge_form monge_form;
+    Monge_form_condition_numbers monge_form_condition_numbers;
+      
+    //gather points around the vertex using rings
+    gather_fitting_points(v, in_points, poly_rings);
+
+    //exit if the nb of points is too small 
+    if ( in_points.size() < min_nb_points )
+      {std::cerr << "Too few points to perform the fitting" << std::endl; exit(1);}
+
+    //For Ridges we need at least 3rd order info
+    assert( d_monge >= 3);
+    // run the main fct : perform the fitting
+    Monge_via_jet_fitting do_it(in_points.begin(), in_points.end(),
+				d_fitting, d_monge, 
+				monge_form, monge_form_condition_numbers);
+    
+    //switch min-max ppal curv/dir wrt the mesh orientation
+    const Vector_3 normal_mesh = P.computeFacetsAverageUnitNormal(v);
+    monge_form.comply_wrt_given_normal(normal_mesh);
+       
+    //Store monge data needed for ridge computations in property maps
+    vertex2d1_pm[v] = monge_form.d1();
+    vertex2d2_pm[v] = monge_form.d2();
+    vertex2k1_pm[v] = monge_form.coefficients()[0];
+    vertex2k2_pm[v] = monge_form.coefficients()[1];
+    vertex2b0_pm[v] = monge_form.coefficients()[2];
+    vertex2b3_pm[v] = monge_form.coefficients()[5];
+    if ( d_monge >= 4) {
+      //= 3*b1^2+(k1-k2)(c0-3k1^3)
+      vertex2P1_pm[v] =
+	3*monge_form.coefficients()[3]*monge_form.coefficients()[3]
+	+(monge_form.coefficients()[0]-monge_form.coefficients()[1])
+	*(monge_form.coefficients()[6]
+	  -3*monge_form.coefficients()[0]*monge_form.coefficients()[0]
+	  *monge_form.coefficients()[0]); 
+      //= 3*b2^2+(k2-k1)(c4-3k2^3)
+      vertex2P2_pm[v] = 
+	3*monge_form.coefficients()[4]*monge_form.coefficients()[4]
+	+(-monge_form.coefficients()[0]+monge_form.coefficients()[1])
+	*(monge_form.coefficients()[10]
+	  -3*monge_form.coefficients()[1]*monge_form.coefficients()[1]
+	  *monge_form.coefficients()[1]); 
+    }
+  } //END FOR LOOP
+}
+
+
+int main()
+{  
+  //load the model from <mesh.off>
+  PolyhedralSurf P;
+  std::ifstream stream("data/ellipsoid.off");
+  stream >> P;
+  fprintf(stderr, "loadMesh %d Ves %d Facets\n",
+	  P.size_of_vertices(), P.size_of_facets());
+  
+  //exit if not enough points in the model
+  if (min_nb_points > P.size_of_vertices())  
+    {std::cerr << "not enough points in the model" << std::endl;   return 1;}
+
+  //initialize Polyhedral data : normal of facets
+  P.compute_facets_normals();
+  
+  //create a Poly_rings object
+  Poly_rings poly_rings(P);
+
+  std::cout << "Compute differential quantities via jet fitting..." << std::endl;
+  //initialize the diff quantities property maps
+  compute_differential_quantities(P, poly_rings);
+  
+  std::cout << "Compute ridges... tag_3" << std::endl;
+  //---------------------------------------------------------------------------
+  //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 BLUE_RIDGE, RED_RIDGE, CREST or all ridges
+  ridge_approximation.compute_ridges(CGAL::BLUE_RIDGE, ii, tag_order);  
+  ridge_approximation.compute_ridges(CGAL::RED_RIDGE, ii, tag_order);  
+  ridge_approximation.compute_ridges(CGAL::CREST_RIDGE, ii, tag_order);  
+  ridge_approximation.compute_all_ridges(ii, tag_order);  
+ 
+  std::cout << "Compute ridges... tag_4" << std::endl;
+  tag_order = Ridge_approximation::Tag_4;
+   //Find BLUE_RIDGE, RED_RIDGE, CREST or all ridges
+  ridge_approximation.compute_ridges(CGAL::BLUE_RIDGE, ii, tag_order);  
+  ridge_approximation.compute_ridges(CGAL::RED_RIDGE, ii, tag_order);  
+  ridge_approximation.compute_ridges(CGAL::CREST_RIDGE, ii, tag_order);  
+  ridge_approximation.compute_all_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);
+  std::cout << "compute umbilics u=1" << std::endl;
+  umbilic_approximation.compute(umb_it, umb_size);
+  umb_size=2;
+  std::cout << "compute umbilics u=2" << std::endl;
+  umb_size=5;
+  std::cout << "compute umbilics u=2" << std::endl;
+ 
+
+ std::vector<Umbilic*>::iterator iter_umb = umbilics.begin(), 
+   iter_umb_end = umbilics.end();
+  // output
+  std::cout << "nb of umbilics " << umbilics.size() << std::endl;
+  for (;iter_umb!=iter_umb_end;iter_umb++) std::cout << **iter_umb;
+
+  return 1;
+}
+