From cb1977cedf4dcdea44f96ea3f02ed878da759df6 Mon Sep 17 00:00:00 2001 From: Marc Pouget Date: Wed, 22 Feb 2006 10:57:13 +0000 Subject: [PATCH] --- .gitattributes | 5 +- Ridges_3/examples/Ridges_3/GSL.h | 109 ++ Ridges_3/examples/Ridges_3/PolyhedralSurf.C | 63 + .../Ridges_3}/PolyhedralSurf.h | 141 +- .../Ridges_3/PolyhedralSurf_operations.h | 34 + .../examples/Ridges_3/PolyhedralSurf_rings.h | 640 +++++++++ Ridges_3/examples/Ridges_3/README | 44 + Ridges_3/examples/Ridges_3/blind.C | 313 +++++ Ridges_3/examples/Ridges_3/options.C | 1147 +++++++++++++++++ Ridges_3/examples/Ridges_3/options.h | 492 +++++++ Ridges_3/include/CGAL/algo.C | 379 ------ 11 files changed, 2894 insertions(+), 473 deletions(-) create mode 100644 Ridges_3/examples/Ridges_3/GSL.h create mode 100644 Ridges_3/examples/Ridges_3/PolyhedralSurf.C rename Ridges_3/{include/CGAL => examples/Ridges_3}/PolyhedralSurf.h (59%) create mode 100644 Ridges_3/examples/Ridges_3/PolyhedralSurf_operations.h create mode 100644 Ridges_3/examples/Ridges_3/PolyhedralSurf_rings.h create mode 100644 Ridges_3/examples/Ridges_3/README create mode 100644 Ridges_3/examples/Ridges_3/blind.C create mode 100644 Ridges_3/examples/Ridges_3/options.C create mode 100644 Ridges_3/examples/Ridges_3/options.h delete mode 100644 Ridges_3/include/CGAL/algo.C diff --git a/.gitattributes b/.gitattributes index ccf7d3398b7..e4140641082 100644 --- a/.gitattributes +++ b/.gitattributes @@ -1254,7 +1254,10 @@ Ridges_3/doc_tex/Ridges_3/ellipsoid_ridges.pdf -text Ridges_3/doc_tex/Ridges_3/ellipsoid_ridges_mesh.eps -text Ridges_3/doc_tex/Ridges_3/ellipsoid_ridges_mesh.jpg -text Ridges_3/doc_tex/Ridges_3/ellipsoid_ridges_mesh.pdf -text -Ridges_3/include/CGAL/algo.C -text +Ridges_3/examples/Ridges_3/PolyhedralSurf.C -text +Ridges_3/examples/Ridges_3/README -text +Ridges_3/examples/Ridges_3/blind.C -text +Ridges_3/examples/Ridges_3/options.C -text Robustness/demo/Robustness/robustness.vcproj -text SearchStructures/doc_tex/SearchStructures/d-range.eps -text SearchStructures/doc_tex/SearchStructures/d-range.gif -text svneol=unset#unset diff --git a/Ridges_3/examples/Ridges_3/GSL.h b/Ridges_3/examples/Ridges_3/GSL.h new file mode 100644 index 00000000000..f2d6f6f8a3d --- /dev/null +++ b/Ridges_3/examples/Ridges_3/GSL.h @@ -0,0 +1,109 @@ +#ifndef _GSL_H_ +#define _GSL_H_ + +#include +#include + +////////////////////////class gsl_Vector///////////////////// +class gsl_Vector{ +protected: + gsl_vector* m_vector; +public: + //contructor + // initializes all the elements of the vector to zero + gsl_Vector(size_t n) { m_vector = gsl_vector_calloc(n); } + //access + double const operator[](size_t i) const { return gsl_vector_get(m_vector, + i); } + //use as left value v[i]=10 + double& operator[](size_t i) { return *gsl_vector_ptr(m_vector, i);} + const gsl_vector* vector() const { return m_vector; } + gsl_vector* vector() { return m_vector; } +}; + +////////////////////////class gsl_Matrix///////////////////// +class gsl_Matrix{ +protected: + gsl_matrix* m_matrix; +public: + //contructor + // initializes all the elements of the matrix to zero. + gsl_Matrix(size_t n1, size_t n2) { m_matrix = gsl_matrix_calloc (n1, n2); } + + //class Row, to define the usual double operator [][] for matrices + class Row{ + gsl_matrix* matrix; + size_t row; + public: + Row(gsl_matrix* _matrix, size_t _row) : + matrix(_matrix), row(_row) {} + double const operator[](size_t column) const + { return gsl_matrix_get(this->matrix, row, column);} + double& operator[](size_t column) + { return *gsl_matrix_ptr(this->matrix, row, column);} + };//END class Row + + Row operator[](size_t _row) { return Row(m_matrix, _row);} + + //access + const gsl_matrix* matrix() const { return m_matrix; } + gsl_matrix* matrix() { return m_matrix; } +}; + +////////////////////////class GSL///////////////////// +class GSL{ +public: + typedef double FT; + typedef gsl_Vector Vector; + typedef gsl_Matrix Matrix; + // solve for eigenvalues and eigenvectors. + // eigen values are sorted in descending order, + // eigen vectors are sorted in accordance. + static + void eigen_symm_algo(Matrix& S, Vector& eval, Matrix& evec); + //solve MX=B using SVD and give the condition number of M + //M replaced by U after gsl_linalg_SV_decomp() + //The diagonal and lower triangular part of M are destroyed during the + //computation, but the strict upper triangular part is not referenced. + static + void solve_ls_svd_algo(Matrix& M, Vector& X, const Vector& B, + double &cond_nb); +}; + +void GSL::eigen_symm_algo(Matrix& S, Vector& eval, Matrix& evec) +{ + const size_t common_size = S.matrix()->size1; + assert( S.matrix()->size2 == common_size + && eval.vector()->size == common_size + && evec.matrix()->size1 == common_size + && evec.matrix()->size2 == common_size + ); + gsl_eigen_symmv_workspace * w = gsl_eigen_symmv_alloc (common_size); + gsl_eigen_symmv (S.matrix(), eval.vector(), evec.matrix(), w); + gsl_eigen_symmv_free (w); + gsl_eigen_symmv_sort (eval.vector(), evec.matrix(), + GSL_EIGEN_SORT_VAL_DESC); +} + +void GSL::solve_ls_svd_algo(Matrix& M, Vector& X, const Vector& B, + double &cond_nb) +{ + const size_t nb_lines = M.matrix()->size1; + assert( B.vector()->size == nb_lines ); + const size_t nb_columns = M.matrix()->size2; + assert( X.vector()->size == nb_columns ); + + gsl_matrix * V = gsl_matrix_alloc(nb_columns,nb_columns); + gsl_vector * S = gsl_vector_alloc(nb_columns); + gsl_vector * work = gsl_vector_alloc(nb_columns); + + gsl_linalg_SV_decomp (M.matrix(), V, S, work); + gsl_linalg_SV_solve (M.matrix(), V, S, B.vector(), X.vector()); + + cond_nb = gsl_vector_get(S,0)/gsl_vector_get(S,nb_columns-1); + gsl_matrix_free(V); + gsl_vector_free(S); + gsl_vector_free(work); +} + +#endif diff --git a/Ridges_3/examples/Ridges_3/PolyhedralSurf.C b/Ridges_3/examples/Ridges_3/PolyhedralSurf.C new file mode 100644 index 00000000000..64f604ec050 --- /dev/null +++ b/Ridges_3/examples/Ridges_3/PolyhedralSurf.C @@ -0,0 +1,63 @@ +#include "PolyhedralSurf.h" + +Vector_3 PolyhedralSurf::getHalfedge_vector(Halfedge * h) +{ + Vector_3 v = h->opposite()->vertex()->point() - h->vertex()->point(); + return v; +} + +///???????myterious behavior???????? +void PolyhedralSurf::compute_facets_normals() +{ + std::for_each(this->facets_begin(), this->facets_end(), + Facet_unit_normal()); +} + +void PolyhedralSurf::compute_edges_length() +{ + std::for_each(this->halfedges_begin(), this->halfedges_end(), + Edge_length()); +} + +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; +} + +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; +} diff --git a/Ridges_3/include/CGAL/PolyhedralSurf.h b/Ridges_3/examples/Ridges_3/PolyhedralSurf.h similarity index 59% rename from Ridges_3/include/CGAL/PolyhedralSurf.h rename to Ridges_3/examples/Ridges_3/PolyhedralSurf.h index b2b3a769fe3..30783e94a5e 100644 --- a/Ridges_3/include/CGAL/PolyhedralSurf.h +++ b/Ridges_3/examples/Ridges_3/PolyhedralSurf.h @@ -26,55 +26,49 @@ template < class Refs, class Tag, class Pt, class FGeomTraits > class My_vertex:public CGAL::HalfedgeDS_vertex_base < Refs, Tag, Pt > { protected: - typedef typename Refs::Vertex Vertex; - typedef typename Refs::Halfedge Halfedge; - typedef typename Refs::Face Facet; - + //to be def in Vertex_wrapper too from the Traits typedef typename FGeomTraits::FT FT; typedef typename FGeomTraits::Point_3 Point_3; typedef typename FGeomTraits::Vector_3 Vector_3; - + protected: - //monge info - Vector_3 d1; //max ppal dir - Vector_3 d2; //min ppal dir; monge normal is then n_m=d1^d2, should be so that n_m.n_mesh>0 - FT k1, k2; //max/min ppal curv - FT b0, b3; //blue/red extremalities - FT P1, P2; //if fourth order quantities + //monge data + Vector_3 m_d1; //max ppal dir + Vector_3 m_d2; //min ppal dir; monge normal is then n_m=d1^d2, should be so that n_m.n_mesh>0 + FT m_k1, m_k2; //max/min ppal curv + FT m_b0, m_b3; //blue/red extremalities + FT m_P1, m_P2; //if fourth order quantities //this is for collecting i-th ring neighbours - char ring_tag; int ring_index; public: - //monge info - const Vector _3 d1() const { return d1; } - Vector_3& d1() { return d1; } - const Vector _3 d2() const { return d2; } - Vector_3& d2() { return d2; } - const FT k1() const { return k1; } - FT& k1() { return k1; } - const FT k2() const { return k2; } - FT& k2() { return k2; } - const FT b0() const { return b0; } - FT& b0() { return b0; } - const FT b3() const { return b3; } - FT& b3() { return b3; } - const FT P1() const { return P1; } - FT& P1() { return P1; } - const FT P2() const { return P2; } - FT& P2() { return P2; } + //monge data + const Vector_3 d1() const { return m_d1; } + Vector_3& d1() { return m_d1; } + const Vector_3 d2() const { return d2; } + Vector_3& d2() { return m_d2; } + const FT k1() const { return m_k1; } + FT& k1() { return m_k1; } + const FT k2() const { return m_k2; } + FT& k2() { return m_k2; } + const FT b0() const { return m_b0; } + FT& b0() { return m_b0; } + const FT b3() const { return m_b3; } + FT& b3() { return m_b3; } + const FT P1() const { return m_P1; } + FT& P1() { return m_P1; } + const FT P2() const { return m_P2; } + FT& P2() { return m_P2; } //this is for collecting i-th ring neighbours void setRingIndex(int i) { ring_index = i; } int getRingIndex() { return ring_index; } void resetRingIndex() { ring_index = -1; } - void setRingTag() { ring_tag = 1; } - char getRingTag() { return ring_tag; } - + My_vertex(const Point_3 & pt): CGAL::HalfedgeDS_vertex_base < Refs, Tag, Point_3 > (pt), - ring_tag(0), ring_index(-1) {} + ring_index(-1) {} My_vertex() {} }; @@ -85,52 +79,42 @@ class My_vertex:public CGAL::HalfedgeDS_vertex_base < Refs, Tag, Pt > template < class Refs, class Tag, class FGeomTraits > class My_facet:public CGAL::HalfedgeDS_face_base < Refs, Tag > { - //protected: public: - typedef typename Refs::Vertex Vertex; - typedef typename Refs::Vertex_handle Vertex_handle; - typedef typename Refs::Halfedge Halfedge; - typedef typename Refs::Halfedge_handle Halfedge_handle; + typedef typename FGeomTraits::Vector_3 Vector_3; - typedef typename FGeomTraits::Vector_3 Vector_3; - typedef typename FGeomTraits::Point_3 Point_3; - -public: +protected: Vector_3 normal; - + int ring_index; + bool m_is_visited; public: My_facet(): ring_index(-1) {} Vector_3 & getUnitNormal() { return normal; } - void setNormal(Vector_3 & n) { normal = n; } + void setNormal(Vector_3 n) { normal = n; } //this is for collecting i-th ring neighbours -protected: - int ring_index; -public: void setRingIndex(int i) { ring_index = i; } int getRingIndex() { return ring_index; } void resetRingIndex() { ring_index = -1; } + //this is for following ridge lines + void set_visited(bool b) { m_is_visited = b; } + const bool is_visited() const { return m_is_visited;} + void reset_is_visited() { m_is_visited = false; } }; //---------------------------------------------------------------- // Halfedge //---------------------------------------------------------------- -template < class Refs, class Tprev, class Tvertex, class Tface, class FGeomTraits > +template < class Refs, class Tprev, class Tvertex, class Tface > class My_halfedge:public CGAL::HalfedgeDS_halfedge_base < Refs, Tprev, Tvertex, Tface > { -protected: - typedef typename FGeomTraits::Point_3 Point_3; - typedef typename FGeomTraits::Vector_3 Vector_3; - protected: int ring_index; double len; public: + My_halfedge(): ring_index(-1) {} void setRingIndex(int i) { ring_index = i; } int getRingIndex() {return ring_index; } void resetRingIndex() {ring_index = -1; } -public: - My_halfedge(): ring_index(-1) {} void setLength(double l) { len = l; } double getLength() { return len; } }; @@ -143,9 +127,11 @@ struct Wrappers_VFH:public CGAL::Polyhedron_items_3 { template < class Refs, class Traits > struct Vertex_wrapper { typedef struct { public: + //typedef typename Traits::Vector_3 Vector_3; + //all types needed by the vertex... + typedef typename Traits::FT FT; typedef typename Traits::Point_3 Point_3; typedef typename Traits::Vector_3 Vector_3; - typedef typename Traits::Plane_3 Plane_3; } FGeomTraits; typedef typename Traits::Point_3 Point_3; typedef My_vertex < Refs, CGAL::Tag_true, Point_3, FGeomTraits > Vertex; @@ -158,65 +144,34 @@ struct Wrappers_VFH:public CGAL::Polyhedron_items_3 { //all types needed by the facet... typedef struct { public: - typedef typename Traits::Point_3 Point_3; + // typedef typename Traits::Point_3 Point_3; typedef typename Traits::Vector_3 Vector_3; - typedef typename Traits::Plane_3 Plane_3; - } FGeomTraits; + } FGeomTraits; //custom type instantiated... typedef My_facet < Refs, CGAL::Tag_true, FGeomTraits > Face; }; // wrap halfedge template < class Refs, class Traits > struct Halfedge_wrapper { - typedef struct { - public: - typedef typename Traits::Point_3 Point_3; - typedef typename Traits::Vector_3 Vector_3; - typedef typename Traits::Plane_3 Plane_3; - } FGeomTraits; - typedef My_halfedge < Refs, + typedef My_halfedge < Refs, CGAL::Tag_true, - CGAL::Tag_true, CGAL::Tag_true, FGeomTraits > Halfedge; + CGAL::Tag_true, CGAL::Tag_true > Halfedge; }; }; //------------------------------------------------ -// Polyhedron +//PolyhedralSurf //------------------------------------------------ -//using standard Cartesian kernel typedef double DFT; typedef CGAL::Cartesian Data_Kernel; - typedef CGAL::Polyhedron_3 < Data_Kernel, Wrappers_VFH > Polyhedron; -typedef Polyhedron::Vertex Vertex; -typedef Polyhedron::Vertex_handle Vertex_handle; -typedef Polyhedron::Halfedge Halfedge; -typedef Polyhedron::Halfedge_handle Halfedge_handle; - -typedef Polyhedron::Vertex_iterator Vertex_iterator; -typedef Polyhedron::Halfedge_iterator Halfedge_iterator; -typedef Polyhedron:: -Halfedge_around_vertex_circulator Halfedge_around_vertex_circulator; -typedef Polyhedron:: -Halfedge_around_facet_circulator Halfedge_around_facet_circulator; -typedef Polyhedron::Facet_iterator Facet_iterator; - -typedef Data_Kernel::Point_3 DPoint; typedef Data_Kernel::Vector_3 Vector_3; -typedef Data_Kernel::Plane_3 Plane_3; -typedef Data_Kernel::Sphere_3 Sphere_3; -/////////////////////class PolyhedralSurf/////////////////////////// class PolyhedralSurf:public Polyhedron { public: - typedef Data_Kernel::Point_3 DPoint; - typedef Data_Kernel::Vector_3 Vector_3; - typedef Data_Kernel::Plane_3 Plane_3; - -public: - static Vector_3 getHalfedge_vector(Halfedge * h); - PolyhedralSurf() {} + + static Vector_3 getHalfedge_vector(Halfedge * h); void compute_edges_length(); double compute_mean_edges_length_around_vertex(Vertex * v); void compute_facets_normals(); diff --git a/Ridges_3/examples/Ridges_3/PolyhedralSurf_operations.h b/Ridges_3/examples/Ridges_3/PolyhedralSurf_operations.h new file mode 100644 index 00000000000..b75ba3b2d0c --- /dev/null +++ b/Ridges_3/examples/Ridges_3/PolyhedralSurf_operations.h @@ -0,0 +1,34 @@ +#ifndef _POLY_OP_H_ +#define _POLY_OP_H_ + +//?????????mysterious , used in Poly.C with mysterious +//std::for_all??????? why arent they member fct of polysurf? + +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 diff --git a/Ridges_3/examples/Ridges_3/PolyhedralSurf_rings.h b/Ridges_3/examples/Ridges_3/PolyhedralSurf_rings.h new file mode 100644 index 00000000000..3edff529840 --- /dev/null +++ b/Ridges_3/examples/Ridges_3/PolyhedralSurf_rings.h @@ -0,0 +1,640 @@ +#ifndef _SURPOLYRINGS_H_ +#define _SURPOLYRINGS_H_ + +using namespace std; + +//-------------------------------------------------- +//we store hedges. exple: 1ring neighbours of +//v=hedges anchored at v +//-------------------------------------------------- +// when storing hedges to a ith ring: do we store hedges pointing to +// the ring or opposite hedges? +enum HTO_HOPPOSITE_ONERING { HTO_ONERING, HOPPOSITE_ONERING }; + +template < class TPoly > class T_PolyhedralSurf_rings +{ +public: + // typedef typename TPoly::Point_3 Point_3; + typedef typename TPoly::Vertex Vertex; + // typedef typename TPoly::Vertex_handle Vertex_handle; + 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; + //debug + static void check_vertex_indices(Vertex* v); + static void check_ring_indices(TPoly& poly); + + // static void from_into(ptrSet & a, ptrSet & b, char erase_a = 0); + + + //---------------VERTEX--------------------------------- + //neighbours using a tag on each vertex ---default is -1, vertex is 0, + // then 1,2,.. for ith rings + //--------------------------------------------------------- + static int + push_neighbours_of(Vertex * start, int ith, + std::vector < Vertex * >&nextRing, + std::vector < Vertex * >&all); + + static int + collect_ith_ring_neighbours(int ith, + std::vector < Vertex * >¤tRing, + std::vector < Vertex * >&nextRing, + std::vector < Vertex * >&all); + + + static void reset_ring_indices(std::vector < Vertex * >&vces); + // static void reset_ring_indices(std::vector< Vertex_handle >&vces); + //Neighbours collected in ccw order on each ring + static int collect_ith_ring_ccw(int ith, + vector& sizes, + vector& vces); + static int i_rings_ccw(int i, + Vertex* start, + vector& sizes, + vector& vces); + + //---------------------------------------------------------------- + // Facet + //---------------------------------------------------------------- + + + + //------------------HEDGES-------------------------------- + //hdges incident to endpoint of h + //we store opposite hedges, ie those + //to the next ith ring neighbours + //--------------------------------------------------------- + static int + push_hedges_of(Vertex * start, int ith, + HTO_HOPPOSITE_ONERING to_opp, + std::vector < Halfedge * >&nextRing, + std::vector < Halfedge * >&recordedOnce); + //, char& multiple); + + //grabbing the next ring from currentRing + static int + next_ring_hedges(int ith, HTO_HOPPOSITE_ONERING to_opp, + std::vector < Halfedge * >¤tRing, + std::vector < Halfedge * >&nextRing, + std::vector < Halfedge * >&recordedOnce); + //, char& multiple); + + //ith ring from scratch ie the center vertex + static int + collect_ith_ring_hedges(Vertex * start, int ith, + HTO_HOPPOSITE_ONERING to_opp, + std::vector < Halfedge * >&recordedOnce, + char cleanup); + + //Collect contour edges of the ith ring, these edges oriented CCW + // i.e. their incident facet is inside. + static void + contour_ith_ring_hedges(Vertex * start, int ith, + std::vector < Halfedge * > &contour); + + static int + contour_ccw_ith_ring_hedges(Vertex * start, int ith, + std::vector < Halfedge * > &contourCCW); + + static void reset_ring_indices(std::vector < Halfedge * >&hedges); + + static void reset_ring_indices(std::vector < Facet * >&facets); + + // to report 1st crossing edge with a sphere + // static Point_3 orig_point; + // static void set_orig_point(Point_3 & p); + // static void + // getFarthestNeighbourOfV_ofOrigPoint(Vertex * v, + // Halfedge * &h_farthest); + +}; + + +////////////////////////////////////////////// +//////////IMLPEMENTATION/////////////////////// +////////////////////////////////////////////// +//-------------------------------------------------- +//-- collecting ith ring neighbours +//-------------------------------------------------- + +template < class TPoly > +void T_PolyhedralSurf_rings :: +check_vertex_indices(Vertex* v) +{ + Halfedge_around_vertex_circulator itb = v->vertex_begin(), + ite = v->vertex_begin(); + assert( ((&(*itb))->opposite()->vertex()->getRingIndex()==-1) ); + itb++; + for( ; itb != ite ; itb++) + assert(((&(*itb))->opposite()->vertex()->getRingIndex()==-1) ); +} + + +template < class TPoly > +void T_PolyhedralSurf_rings < TPoly >:: +check_ring_indices(TPoly& poly) +{ + Vertex_iterator itb = poly.vertices_begin (), + ite = poly.vertices_end(); + for( ; itb != ite ; itb++) + T_PolyhedralSurf_rings < TPoly >::check_vertex_indices( &(*itb)); +} + +// //neighbours using a map to make sure we do not collect twice +// //--------------------------------------------------------- +// template < class TPoly > +// void T_PolyhedralSurf_rings < TPoly >:: +// from_into(ptrSet & a, ptrSet & b, char erase_a) +// { +// ptrSet_iterator itb = a.begin(), ite = a.end(); +// for (; itb != ite; itb++) +// b.insert(*itb); + +// if (erase_a) +// a.erase(a.begin(), ite); +// } + +//neighbours using a tag on each vertex ---default is -1, vertex is 0, +// then 1,2,.. for ith rings +//------------------------------------------------------------------- + +template < class TPoly > +int T_PolyhedralSurf_rings < TPoly >:: +push_neighbours_of(Vertex * start, int ith, + std::vector < Vertex * >&nextRing, + std::vector < Vertex * >&all) +{ + int nb = 0; + Vertex *v; + Halfedge *h; + + Halfedge_around_vertex_circulator + hedgeb = start->vertex_begin(), hedgee = hedgeb; + + do + { + h = &(*hedgeb); + hedgeb++; + v = &(*h->opposite()->vertex()); + + //make sure vertex not visited yet + if (v->getRingIndex() != -1) + continue; + //tag the vertex as belonging to the its ring + v->setRingIndex(ith); + nb++; + nextRing.push_back(v); + all.push_back(v); + } + while (hedgeb != hedgee); + return nb; +} + +template < class TPoly > +int T_PolyhedralSurf_rings < TPoly >:: +collect_ith_ring_neighbours(int ith, std::vector < Vertex * >¤tRing, + std::vector < Vertex * >&nextRing, + std::vector < Vertex * >&all) +{ + Vertex *v; + typename std::vector < Vertex * >::iterator itb = + currentRing.begin(), ite = currentRing.end(); + + for (; itb != ite; itb++) + { + v = *itb; + T_PolyhedralSurf_rings < TPoly >:: + push_neighbours_of(v, ith, nextRing, all); + } + return nextRing.size(); +} + +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(); + for (; itb != ite; itb++) + (*itb)->resetRingIndex(); +} + +// //surcharge pour les handles +// template < class TPoly > +// void T_PolyhedralSurf_rings < TPoly >:: +// reset_ring_indices(std::vector < Vertex_handle >&vces) +// { +// typename std::vector < Vertex_handle >::iterator itb = vces.begin(), ite = +// vces.end(); +// for (; itb != ite; itb++) +// (&(*(*itb)))->resetRingIndex(); +// } + +//see also the other method below-------------- +//Neighbours collected in ccw order on each ring +//if a border edge is encountered, return 0 +//sizes[j] = nb of vertives of the jth ring +// template +// int T_PolyhedralSurf_rings:: +// collect_ith_ring_ccw(int ith, +// vector& sizes, +// vector& vces) +// { +// //find a vertex of the ith ring from a vertex of the (i-1)th ring +// int size=0; +// Vertex* v=vces[vces.size()-1]; +// Halfedge* h; +// Halfedge_around_vertex_circulator +// hedgeb = v->vertex_begin(), hedgee = hedgeb; +// h = &(*hedgeb); +// if (h->is_border_edge()) return 0; +// while (h->opposite()->vertex()->getRingIndex() != -1) +// { +// hedgeb++; +// h = &(*hedgeb); +// if (h->is_border_edge()) return 0; +// } +// v = &(*(h->opposite()->vertex())); +// vces.push_back(v); +// size++; +// v->setRingIndex(ith); +// h=&(*h->opposite()->next()); +// if (h->is_border_edge()) return 0; +// v=&(*h->vertex()); +// //follow ccw the ith ring +// while (v->getRingIndex() != ith) +// { +// while ( (v->getRingIndex() != -1) && (v->getRingIndex() != ith) ) +// { +// h=&(*h->opposite()->next()); +// if (h->is_border_edge()) return 0; +// v=&(*h->vertex()); +// } +// if (v->getRingIndex() == -1) +// { +// vces.push_back(v); +// size++; +// v->setRingIndex(ith); +// h=&(*h->next()); +// if (h->is_border_edge()) return 0; +// v=&(*h->vertex()); +// } +// } +// sizes.push_back(size); +// return 1; +// } +//Other method using contour_ccw_ith_ring_hedges-------------------------- +template +int T_PolyhedralSurf_rings:: +collect_ith_ring_ccw(int ith, + vector& sizes, + vector& vces) +{ + Vertex* v=vces[0]; + std::vector < Halfedge * > contourCCW; + //check if no border edge is encountered and set contourCCW of halfedges + if (T_PolyhedralSurf_rings:: + contour_ccw_ith_ring_hedges(v, ith, contourCCW) == 0) return 0; + + typename std::vector < Halfedge * >::iterator itb = + contourCCW.begin(), ite = contourCCW.end(); + + //retrieve vertices from hedge contour + for (; itb != ite; itb++) vces.push_back(&(*((*itb)->vertex()))); + sizes.push_back(contourCCW.size()); + return 1; +} + + +template +int T_PolyhedralSurf_rings:: +i_rings_ccw(int i, + Vertex* start, + vector& sizes, + vector& vces) +{ + //set the 0-ring + start->setRingIndex(0); + sizes.push_back(1); + vces.push_back(start); + int ith; + //set ith ring, ith from 1 to i + for(ith=1;ith<=i;ith++) + if (T_PolyhedralSurf_rings:: + collect_ith_ring_ccw(ith, sizes, vces) == 0) + { + T_PolyhedralSurf_rings::reset_ring_indices(vces); + return 0; + } + T_PolyhedralSurf_rings::reset_ring_indices(vces); + return 1; +} + +//----------------------------------------------------- +//---------- hedges. +//----------------------------------------------------- +template < class TPoly > +int T_PolyhedralSurf_rings < TPoly >:: +push_hedges_of(Vertex * start, int ith, HTO_HOPPOSITE_ONERING to_opp, + std::vector < Halfedge * >&nextRing, + std::vector < Halfedge * >&recordedOnce) +{ + // char ok; + int n = 0; + Halfedge *h, *hopp, *hpushed; + Halfedge_around_vertex_circulator + hedgeb = start->vertex_begin(), hedgee = hedgeb; + do + { + h = &(*hedgeb); + hopp = &(*h->opposite()); + hedgeb++; + //store either h or hopp + hpushed = (to_opp == HTO_ONERING) ? hopp : h; + + //in any case + nextRing.push_back(hpushed); + + //record or not? + if (hpushed->getRingIndex() == -1) + { + //tag the vertex as belonging to the its ring + //marc: isn't it: tag the edge and its opposite? + h->setRingIndex(ith); + hopp->setRingIndex(ith); + n++; //xfc + recordedOnce.push_back(hpushed); + } + //hedge already recorded + //else multiple=1; + } + while (hedgeb != hedgee); + return n; +} + +//notice that for each helfedge, we need to recover the endpoint, ie +//the vertex on the current ring being explored +template < class TPoly > + int T_PolyhedralSurf_rings < TPoly >:: +next_ring_hedges(int ith, HTO_HOPPOSITE_ONERING to_opp, + std::vector < Halfedge * >¤tRing, + std::vector < Halfedge * >&nextRing, + std::vector < Halfedge * >&recordedOnce) +{ + int s = 0; + Vertex *v; + Halfedge *h; + typename std::vector < Halfedge * >::iterator itb = + currentRing.begin(), ite = currentRing.end(); + for (; itb != ite; itb++) + { + h = (Halfedge *) (*itb); + + //grab vertex + if (to_opp == HTO_ONERING) + v = &*(h->vertex()); + else + v = &*(h->opposite()->vertex()); + + s += T_PolyhedralSurf_rings < TPoly >:: + push_hedges_of(v, ith, to_opp, nextRing, recordedOnce); + } + return s; +} + +//Warning collect hedges from rings 1 to ith, better to call it +//collect_i_rings_hedges ??????? +template < class TPoly > +int T_PolyhedralSurf_rings < TPoly >:: +collect_ith_ring_hedges(Vertex * start, int ith, + HTO_HOPPOSITE_ONERING to_opp, + std::vector < Halfedge * >&recordedOnce, + char cleanup) +{ + // char multiple=0; + int i = 1, np; + std::vector < Halfedge * >baseRing, nextRing; + // typename std::vector::iterator itb, ite; + std::vector < Halfedge * >*p_baseRing = &baseRing, *p_nextRing = + &nextRing; + + //collect the one ring + np = T_PolyhedralSurf_rings < TPoly >:: + push_hedges_of(start, int (1), to_opp, nextRing, recordedOnce); + + while (1) + { + //exit? + i++; + if (i > ith) + break; + + //find next ring + p_baseRing->erase(p_baseRing->begin(), p_baseRing->end()); + std::swap(p_baseRing, p_nextRing); + + np = T_PolyhedralSurf_rings < TPoly >:: + next_ring_hedges(i, HTO_ONERING, *p_baseRing, *p_nextRing, + recordedOnce); + } + //clean up hedges tags + np = recordedOnce.size(); //#triangles + if (cleanup) + T_PolyhedralSurf_rings < TPoly >::reset_ring_indices(recordedOnce); + return np; +} + +//Warning contour edges are not given in a cyclic order +template < class TPoly > +void T_PolyhedralSurf_rings < TPoly >:: +contour_ith_ring_hedges(Vertex * start, int ith, + std::vector < Halfedge * > &contour) +{ + std::vector < Halfedge * > ringHedges; + typename std::vector < Halfedge * >::iterator itb = + ringHedges.begin(), ite = ringHedges.end(); + T_PolyhedralSurf_rings < TPoly >:: + collect_ith_ring_hedges(start, ith, HTO_ONERING, ringHedges, 0); + //debug + // fprintf(stderr, "\n"); + // int iii; + // for( iii=0;iiigetRingIndex(), + // &(*( ringHedges[iii]->vertex() )), + // &(*(ringHedges[iii]->opposite()->vertex() )) ); + //debug + + + //set ringHedges index to i + // for (; itb != ite; itb++) (*itb)->setRingIndex(ith);non! garder + // les index!! + + //for a hedge h in ringHedges and h index = ith, h->next is on the + // contour if its opposite is not visited + //Note that contour edges are then oriented CCW + itb = ringHedges.begin(), ite = ringHedges.end(); + for (; itb != ite; itb++) + { + if ( ( (*itb)->getRingIndex()==ith )&& + ((*itb)->next()->opposite()->getRingIndex() == -1) ) + contour.push_back( &(*((*itb)->next())) ); + } + //debug + // fprintf(stderr, "\n\n"); + // for( iii=0;iiigetRingIndex(), + // &(*( contour[iii]->vertex() )), + // &(*(contour[iii]->opposite()->vertex() )) ); + // //debug + + //reset indices + T_PolyhedralSurf_rings < TPoly >::reset_ring_indices(ringHedges); +} + +//if the contour is a circle, then put edges head to tail, ccw +template < class TPoly > +int T_PolyhedralSurf_rings < TPoly >:: +contour_ccw_ith_ring_hedges(Vertex * start, int ith, + std::vector < Halfedge * > &contourCCW) +{ + std::vector < Halfedge * > contour; + T_PolyhedralSurf_rings < TPoly >:: + contour_ith_ring_hedges(start, ith, contour); + + //check if no border edge is encountered + typename std::vector < Halfedge * >::iterator itb = + contour.begin(), ite = contour.end(); + for (; itb != ite; itb++) if ( (*itb)->is_border_edge()) return 0; + + int size = contour.size(); + Halfedge* h_cur, *h; + h_cur = contour[0]; + // fprintf(stderr, "\n\n"); + // int ii; + // for( ii=0, itb=contour.begin();iivertex()); + // fprintf(stderr, " %d] ", (*itb)->vertex());} + + // fprintf(stderr, "\n\n"); + contourCCW.push_back(h_cur); + do + { + for (itb=contour.begin(); itb != ite; itb++) + {//debug + h = &(*( ((*itb)->opposite()) )); + // cerr << "h_cur->vertex() " << h_cur->vertex() << endl; + // << "h->vertex() " << h->vertex() << endl; + if ( (h_cur->vertex()) == (h->vertex()) ) + {//the edge following the current edge is found + contourCCW.push_back(*itb); + h_cur = *itb; + break; + } + }//debug + } + while (contourCCW.size() != size); + + //debug + // int iii; + + // fprintf(stderr, "\n"); + + // for( iii=0;iiigetRingIndex(), + // &(*( contour[iii]->vertex() )), + // &(*(contour[iii]->opposite()->vertex() )) ); + // //debug + // //check the contour debug! //debug + // fprintf(stderr, "\n"); + // for( iii=0;iiigetRingIndex(), + // &(*( contourCCW[iii]->vertex() )), + // &(*(contourCCW[iii]->opposite()->vertex() )) ); + // //debug + + for (int i=0; ivertex())) == + &(*(contourCCW[(i+1) % size]->opposite()->vertex())) ); + } + return 1; +} + + +template < class TPoly > +void T_PolyhedralSurf_rings < TPoly >:: +reset_ring_indices(std::vector < Halfedge * >&hedges) +{ + typename std::vector < Halfedge * >::iterator itb = + hedges.begin(), ite = hedges.end(); + for (; itb != ite; itb++) + { + (*itb)->resetRingIndex(); + (*itb)->opposite()->resetRingIndex(); + } +} + +template < class TPoly > +void T_PolyhedralSurf_rings < TPoly >:: +reset_ring_indices(std::vector < Facet * >&facets) +{ + typename std::vector < Facet * >::iterator itb = facets.begin(), ite = + facets.end(); + for (; itb != ite; itb++) + (*itb)->resetRingIndex(); +} + +// // to report 1st crossing edge with a sphere +// template < class TPoly > +// typename TPoly::Point_3 T_PolyhedralSurf_rings < TPoly >::orig_point; + +// template < class TPoly > +// void T_PolyhedralSurf_rings < TPoly >::set_orig_point(Point_3 & p) +// { +// T_PolyhedralSurf_rings < TPoly >::orig_point = p; +// } + +// template < class TPoly > +// void T_PolyhedralSurf_rings < TPoly >:: +// getFarthestNeighbourOfV_ofOrigPoint(Vertex * v, Halfedge * &h_farthest) +// { +// double d, dref; +// Halfedge *h; +// Halfedge_around_vertex_circulator +// hedgeb = v->vertex_begin(), hedgee = hedgeb; + +// //init +// h = h_farthest = &(*hedgeb); +// hedgeb++; +// dref = squared_distance(h->opposite()->vertex()->point(), +// T_PolyhedralSurf_rings < TPoly >::orig_point); +// //Visit neighbours +// do +// { +// h = &(*hedgeb); +// hedgeb++; + +// assert(v == (&(*h->vertex()))); + +// d = squared_distance(h->opposite()->vertex()->point(), +// T_PolyhedralSurf_rings < TPoly >::orig_point); +// if (d > dref) +// { +// dref = d; +// h_farthest = h; +// } +// } +// while (hedgeb != hedgee); +// } + +#endif diff --git a/Ridges_3/examples/Ridges_3/README b/Ridges_3/examples/Ridges_3/README new file mode 100644 index 00000000000..8317138fe9b --- /dev/null +++ b/Ridges_3/examples/Ridges_3/README @@ -0,0 +1,44 @@ +Program blind_1pt.exe +--------------------- +takes as input a which is a xyz text file +it compute the fitting for this single entry and +it outputs results in the file and on the standard std::cout + +Usage is : main , " + +./blind_1pt.exe data/in_points_file.txt output.txt 2 2 + + + +Program blind.exe +----------------- +takes an filename.off file as input, +it computes a fitting for each vertex + +it outputs the results in : +1. filename.off.4ogl.txt which can be visualized with + the demo program visu.exe +2. if option -v, filename.off.verb.txt contains human readable results + +Usage is : blind with options + "f:fName ", // filename.off + "d:deg ", // degree of the jet + "m:mdegree ", // degree of the Monge rep + "a:nrings ", // # rings + 0 means collect enough rings to make appro possible + k>=1 fixes the nb of rings to be collected + "p:npoints ", //# points + 0 means this option is not considered, this is the default + n>=1 fixes the nb of points to be used + "v|", // verbose? + +Note : if the nb of collected points is less than the required min number of + points to make the approxiamtion possible (which is constrained by the deg) + then the vertex is skipped. + + ./blind.exe -f data/ellipe0.003.off -d2 -m2 -a2 + ./blind.exe -f data/poly2x\^2+y\^2-0.062500-off -d2 -m2 -a2 + +visu with: + +./visu.exe ../../examples/Jet_fitting_3/data/poly2x\^2+y\^2-0.062500-off ../../examples/Jet_fitting_3/data_poly2x\^2+y\^2-0.062500-off.4ogl.txt diff --git a/Ridges_3/examples/Ridges_3/blind.C b/Ridges_3/examples/Ridges_3/blind.C new file mode 100644 index 00000000000..f623ed4b96b --- /dev/null +++ b/Ridges_3/examples/Ridges_3/blind.C @@ -0,0 +1,313 @@ +#include +#include + +#include +#include +#include +#include +#include +#include + +#include "../../include/CGAL/Ridges.h" +#include "../../../Jet_fitting_3/include/CGAL/Monge_via_jet_fitting.h" +#include "GSL.h" + +#include "PolyhedralSurf.h" +#include "PolyhedralSurf_rings.h" +#include "options.h"//parsing command line + +//Kernel of the PolyhedralSurf +typedef double DFT; +typedef CGAL::Cartesian Data_Kernel; +typedef Data_Kernel::Point_3 DPoint; +typedef Data_Kernel::Vector_3 DVector; +typedef PolyhedralSurf::Vertex Vertex; +typedef PolyhedralSurf::Vertex_iterator Vertex_iterator; +typedef T_PolyhedralSurf_rings Poly_rings; +//Kernel for local computations +typedef double LFT; +typedef CGAL::Cartesian Local_Kernel; +typedef CGAL::Monge_via_jet_fitting My_Monge_via_jet_fitting; +typedef CGAL::Monge_rep My_Monge_rep; +typedef CGAL::Monge_info My_Monge_info; + +//RIDGES +typedef CGAL::Ridge_line Ridge_line; +typedef CGAL::Ridge_approximation::iterator > Ridge_approximation; + +//Syntax requirred by Options +static const char *const optv[] = { + "?|?", + "f:fName ", //name of the input off file + "d:deg ", //degree of the jet + "m:mdegree ", //degree of the Monge rep + "a:nrings ", //# rings + "p:npoints ", //# points + "v|",//verbose? + NULL +}; + + +int main(int argc, char *argv[]) +{ + // fct parameters + unsigned int d_fitting = 2; + unsigned int d_monge = 2; + unsigned int nb_rings = 0;//seek min # of rings to get the required #pts + unsigned int nb_points_to_use = 0;// + bool verbose = false; + + char *if_name = NULL, //input file name + *w_if_name = NULL; //as above, but / replaced by _ + char* res4openGL_fname; + char* verbose_fname; + std::ofstream *out_4ogl = NULL, *out_verbose = NULL; + + int optchar; + char *optarg; + Options opts(*argv, optv); + OptArgvIter iter(--argc, ++argv); + + while ((optchar = opts(iter, (const char *&) optarg))){ + switch (optchar){ + case 'f': if_name = optarg; break; + case 'd': d_fitting = atoi(optarg); break; + case 'm': d_monge = atoi(optarg); break; + case 'a': nb_rings = atoi(optarg); break; + case 'p': nb_points_to_use = atoi(optarg); break; + case 'v': verbose=true; break; + default: + cerr << "Unknown command line option " << optarg; + exit(0); + } + } + + //prepare output file names + assert(if_name != NULL); + w_if_name = new char[strlen(if_name)]; + strcpy(w_if_name, if_name); + for(unsigned int i=0; i + PolyhedralSurf P; + std::ifstream stream(if_name); + stream >> P; + fprintf(stderr, "loadMesh %d Ves %d Facets\n", + P.size_of_vertices(), P.size_of_facets()); + if(verbose) + (*out_verbose) << "Polysurf with " << P.size_of_vertices() + << " vertices and " << P.size_of_facets() + << " facets. " << std::endl; + + //exit if not enough points in the model + unsigned int min_nb_points = (d_fitting + 1) * (d_fitting + 2) / 2; + if (min_nb_points > P.size_of_vertices()) exit(0); + //initialize Polyhedral data : length of edges, normal of facets + P.compute_edges_length(); + P.compute_facets_normals(); + + //container for approximation points + std::vector in_points; + //container to collect vertices of the PolyhedralSurf + std::vector current_ring, next_ring, gathered; + std::vector *p_current_ring, *p_next_ring; + + //MAIN LOOP + Vertex_iterator vitb = P.vertices_begin(), vite = P.vertices_end(); + for (; vitb != vite; vitb++) { + //initialize + Vertex* v = &(*vitb);//handle to * + unsigned int nbp = 0, //current nb of collected points + ith = 0; //i-th ring index + current_ring.clear(); + next_ring.clear(); + p_current_ring = ¤t_ring; + p_next_ring = &next_ring; + gathered.clear(); + in_points.clear(); + My_Monge_rep monge_rep; + My_Monge_info monge_info; + + //DO NOT FORGET TO UNTAG AT THE END! + v->setRingIndex(ith); + //collect 0th ring : the vertex v! + gathered.push_back(v); + nbp = 1; + //collect 1-ring + ith = 1; + nbp = Poly_rings::push_neighbours_of(v, ith, current_ring, gathered); + //collect more neighbors depending on options... + + //OPTION -p nb, with nb != 0 + //for approximation/interpolation with a fixed nb of points, 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 ) { + while( gathered.size() < nb_points_to_use ) { + ith++; + //using tags + nbp += Poly_rings:: + collect_ith_ring_neighbours(ith, *p_current_ring, + *p_next_ring, gathered); + //next round must be launched from p_nextRing... + p_current_ring->clear(); + std::swap(p_current_ring, p_next_ring); + } + //clean up + Poly_rings::reset_ring_indices(gathered); + //discard non-required collected points of the last ring + gathered.resize(nb_points_to_use, NULL); + assert(gathered.size() == nb_points_to_use ); + } + else{ + //OPTION -a nb, with nb = 0 + // select the mini nb of rings needed to make approx possible + if (nb_rings == 0) { + while ( gathered.size() < min_nb_points ) { + ith++; + nbp += Poly_rings:: + collect_ith_ring_neighbours(ith, *p_current_ring, + *p_next_ring, gathered); + //next round must be launched from p_nextRing... + p_current_ring->clear(); + std::swap(p_current_ring, p_next_ring); + } + } + //OPTION -a nb, with nb = 1, nothing to do! we have already + // collected the 1 ring + //OPTION -a nb, with nb > 1 + //for approximation with a fixed nb of rings, collect + // neighbors up to the "nb_rings"th ring + if (nb_rings > 1) + while (ith < nb_rings) { + ith++; + nbp += Poly_rings:: + collect_ith_ring_neighbours(ith, *p_current_ring, + *p_next_ring, gathered); + //next round must be launched from p_nextRing... + p_current_ring->clear(); + std::swap(p_current_ring, p_next_ring); + } + + //clean up + Poly_rings::reset_ring_indices(gathered); + } //END ELSE + + //skip if the nb of gathered points is to small + if ( gathered.size() < min_nb_points ) continue; + + //store the gathered points + std::vector::iterator itb = gathered.begin(), + ite = gathered.end(); + CGAL_For_all(itb,ite) in_points.push_back((*itb)->point()); + + // run the main fct : perform the fitting + My_Monge_via_jet_fitting do_it(in_points.begin(), in_points.end(), + d_fitting, d_monge, + monge_rep, monge_info); + + //switch min-max ppal curv/dir wrt the mesh orientation + DVector normal_mesh = P.computeFacetsAverageUnitNormal(v); + DVector normal_monge = monge_rep.n(); + if ( normal_mesh*normal_monge < 0 ) + { + monge_rep.n() = -monge_rep.n(); + std::swap(monge_rep.d1(), monge_rep.d2()); + if ( d_monge >= 2) + std::swap(monge_rep.coefficients()[0],monge_rep.coefficients()[1]); + if ( d_monge >= 3) { + std::swap(monge_rep.coefficients()[2],monge_rep.coefficients()[5]); + std::swap(monge_rep.coefficients()[3],monge_rep.coefficients()[4]);} + if ( d_monge >= 4) { + std::swap(monge_rep.coefficients()[6],monge_rep.coefficients()[10]); + std::swap(monge_rep.coefficients()[7],monge_rep.coefficients()[9]);} + std::vector::iterator itb = monge_rep.coefficients().begin(), + ite = monge_rep.coefficients().end(); + for (;itb!=ite;itb++) { *itb = -(*itb); } + } + + //OpenGL output + //scaling for ppal dir, + // may be optimized with a global mean edges length computed only once + // on all edges of P + DFT scale_ppal_dir = P.compute_mean_edges_length_around_vertex(v)/2; + (*out_4ogl) << v->point() << " " + << monge_rep.origin_pt() << " " + << monge_rep.d1() * scale_ppal_dir << " " + << monge_rep.d2() * scale_ppal_dir << " " + << monge_rep.coefficients()[0] << " " + << monge_rep.coefficients()[1] << " " + << std::endl; + + //verbose txt output + if (verbose){ + (*out_verbose) << "--- vertex " << ++nb_vertices_considered + << " : " << v->point() << std::endl + << "number of points used : " << in_points.size() << std::endl + << "number of rings used : " << ith << std::endl + << "origin : " << monge_rep.origin_pt() << std::endl + << "d1 : " << monge_rep.d1() << std::endl + << "d2 : " << monge_rep.d2() << std::endl + << "n : " << monge_rep.n() << std::endl + << "cond_nb : " << monge_info.cond_nb() << std::endl + << "pca_eigen_vals " << monge_info.pca_eigen_vals()[0] + << " " << monge_info.pca_eigen_vals()[2] + << " " << monge_info.pca_eigen_vals()[3] << std::endl + << "pca_eigen_vecs : " << std::endl + << monge_info.pca_eigen_vecs()[0] << std::endl + << monge_info.pca_eigen_vecs()[1] << std::endl + << monge_info.pca_eigen_vecs()[2] << std::endl; + if ( d_monge >= 2) + (*out_verbose) << "k1 : " << monge_rep.coefficients()[0] << std::endl + << "k2 : " << monge_rep.coefficients()[1] << std::endl + << std::endl ; + } //END IF + } //END FOR LOOP + + /////////////////////////////////////////// + //RIDGES!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + /////////////////////////////////////////// + Ridge_approximation ridge_approximation; + std::list ridge_lines; + // std::back_insert_iterator > iter_lines(ridge_lines), iter_end; + std::list::iterator iter_lines = ridge_lines.begin(), iter_end; + iter_end = ridge_approximation.do_it(P, iter_lines); + + + + + + + // std::cout << "ok 1 " << std::endl; + //cleanup filenames + delete res4openGL_fname; + out_4ogl->close(); + delete out_4ogl; + if(verbose) { + delete verbose_fname; + out_verbose->close(); + delete out_verbose; + } + // std::cout << "ok 2 " << std::endl ; + return 1; +} + diff --git a/Ridges_3/examples/Ridges_3/options.C b/Ridges_3/examples/Ridges_3/options.C new file mode 100644 index 00000000000..3cf059804a9 --- /dev/null +++ b/Ridges_3/examples/Ridges_3/options.C @@ -0,0 +1,1147 @@ +// **************************************************************************** +// ^FILE: options.c - implement the functions defined in +// +// ^HISTORY: +// 01/16/92 Brad Appleton Created +// +// 03/23/93 Brad Appleton +// - Added OptIstreamIter class +// +// 10/08/93 Brad Appleton +// - Added "hidden" options +// +// 02/08/94 Brad Appleton +// - Added "OptionSpec" class +// - Permitted use of stdio instead of iostreams via #ifdef USE_STDIO +// +// 03/08/94 Brad Appleton +// - completed support for USE_STDIO +// - added #ifdef NO_USAGE for people who always want to print their own +// - Fixed stupid NULL pointer error in OptionsSpec class +// +// 07/31/97 Brad Appleton +// - Added PARSE_POS control flag and POSITIONAL return value. +// ^^************************************************************************** + +#ifdef USE_STDIO +# include +#else +# include +#endif + +// #include +#include +#include + +#include "options.h" + +extern "C" { + void exit(int); +} + +static const char ident[] = "@(#)Options 1.05" ; + + // I need a portable version of "tolower" that does NOT modify + // non-uppercase characters. + // +#define TOLOWER(c) (isupper(c) ? tolower(c) : c) + + // Use this to shut the compiler up about NULL strings +#define NULLSTR (char *)NULL + +// ******************************************************** insertion operators + + // If you are using then you need this stuff! + // If you are using then #ifdef this stuff out + // + + +#ifdef USE_STDIO + + // Implement just enough of ostream to get this file to compile + // + +static const char endl = '\n' ; + +class ostream { +public: + ostream(FILE * fileptr) : fp(fileptr) {} + + ostream & + operator<<(char ch); + + ostream & + operator<<(const char * str); + + ostream & + write(const char * buf, unsigned bufsize); + +private: + FILE * fp; +} ; + +ostream & +ostream::operator<<(char ch) { + fputc(ch, fp); + return *this; +} + +ostream & +ostream::operator<<(const char * str) { + fputs(str, fp); + return *this; +} + +ostream & +ostream::write(const char * buf, unsigned ) { + fputs(buf, fp); + return *this; +} + +static ostream cerr(stderr); +static ostream cout(stdout); + +#endif /* USE_STDIO */ + +// ************************************************************** OptIter + +OptIter::~OptIter(void) {} + +const char * +OptIter::operator()(void) { + const char * elt = curr(); + (void) next(); + return elt; +} + +// ************************************************************** OptIterRwd + +OptIterRwd::OptIterRwd(void) {} + +OptIterRwd::~OptIterRwd(void) {} + +// ************************************************************** OptArgvIter + +void OptArgvIter::init(int argc, const char * const argv[]) +{av=argv; ac=argc; ndx=0;} + +OptArgvIter::~OptArgvIter(void) {} + +const char * +OptArgvIter::curr(void) { + return ((ndx == ac) || (av[ndx] == NULL)) ? NULLSTR : av[ndx]; +} + +void +OptArgvIter::next(void) { + if ((ndx != ac) && av[ndx]) ++ndx; +} + +const char * +OptArgvIter::operator()(void) { + return ((ndx == ac) || (av[ndx] == NULL)) ? NULLSTR : av[ndx++]; +} + +void +OptArgvIter::rewind(void) { ndx = 0; } + +// ************************************************************** OptStrTokIter + +static const char WHITESPACE[] = " \t\n\r\v\f" ; +const char * OptStrTokIter::default_delims = WHITESPACE ; + +OptStrTokIter::OptStrTokIter(const char * tokens, const char * delimiters) + : len(unsigned(strlen(tokens))), str(tokens), seps(delimiters), + cur(NULLSTR), tokstr(NULLSTR) +{ + if (seps == NULL) seps = default_delims; + tokstr = new char[len + 1]; + (void) ::strcpy(tokstr, str); + cur = ::strtok(tokstr, seps); +} + + +OptStrTokIter::~OptStrTokIter(void) { delete [] tokstr; } + +const char * +OptStrTokIter::curr(void) { return cur; } + +void +OptStrTokIter::next(void) { if (cur) cur = ::strtok(NULL, seps); } + +const char * +OptStrTokIter::operator()(void) { + const char * elt = cur; + if (cur) cur = ::strtok(NULL, seps); + return elt; +} + +void +OptStrTokIter::rewind(void) { + (void) ::strcpy(tokstr, str); + cur = ::strtok(tokstr, seps); +} + +// ************************************************************* OptIstreamIter + +#ifdef vms + enum { c_COMMENT = '!' } ; +#else + enum { c_COMMENT = '#' } ; +#endif + +const unsigned OptIstreamIter::MAX_LINE_LEN = 1024 ; + + // Constructor +OptIstreamIter::OptIstreamIter(istream & input) : is(input), tok_iter(NULL) +{ +#ifdef USE_STDIO + fprintf(stderr, "%s: Can't use OptIstreamIter class:\n", + "OptIstreamIter::OptIstreamIter"); + fprintf(stderr, "\tOptions(3C++) was compiled with USE_STDIO #defined.\n"); + exit(-1); +#endif /* USE_STDIO */ +} + + // Destructor +OptIstreamIter::~OptIstreamIter(void) { + delete tok_iter; +} + +const char * +OptIstreamIter::curr(void) { +#ifdef USE_STDIO + return NULLSTR; +#else + const char * result = NULLSTR; + if (tok_iter) result = tok_iter->curr(); + if (result) return result; + fill(); + return (! is) ? NULLSTR : tok_iter->curr(); +#endif /* USE_STDIO */ +} + +void +OptIstreamIter::next(void) { +#ifdef USE_STDIO + return; +#else + const char * result = NULLSTR; + if (tok_iter) result = tok_iter->operator()(); + if (result) return; + fill(); + if (! is) tok_iter->next(); +#endif /* USE_STDIO */ +} + +const char * +OptIstreamIter::operator()(void) { +#ifdef USE_STDIO + return NULLSTR; +#else + const char * result = NULLSTR; + if (tok_iter) result = tok_iter->operator()(); + if (result) return result; + fill(); + return (! is) ? NULLSTR : tok_iter->operator()(); +#endif /* USE_STDIO */ +} + + // What we do is this: for each line of text in the istream, we use + // a OptStrTokIter to iterate over each token on the line. + // + // If the first non-white character on a line is c_COMMENT, then we + // consider the line to be a comment and we ignore it. + // +void +OptIstreamIter::fill(void) { +#ifdef USE_STDIO + return; +#else + char buf[OptIstreamIter::MAX_LINE_LEN]; + do { + *buf = '\0'; + is.getline(buf, sizeof(buf)); + char * ptr = buf; + while (isspace(*ptr)) ++ptr; + if (*ptr && (*ptr != c_COMMENT)) { + delete tok_iter; + tok_iter = new OptStrTokIter(ptr); + return; + } + } while (is); +#endif /* USE_STDIO */ +} + +// **************************************************** Options class utilities + + // Is this option-char null? +inline static int +isNullOpt(char optchar) { + return ((! optchar) || isspace(optchar) || (! isprint(optchar))); +} + + // Check for explicit "end-of-options" +inline static int +isEndOpts(const char * token) { + return ((token == NULL) || (! ::strcmp(token, "--"))) ; +} + + // See if an argument is an option +inline static int +isOption(unsigned flags, const char * arg) { + return (((*arg != '\0') || (arg[1] != '\0')) && + ((*arg == '-') || ((flags & Options::PLUS) && (*arg == '+')))) ; +} + + // See if we should be parsing only options or if we also need to + // parse positional arguments +inline static int +isOptsOnly(unsigned flags) { + return (flags & Options::PARSE_POS) ? 0 : 1; +} + + // return values for a keyword matching function +enum kwdmatch_t { NO_MATCH, PARTIAL_MATCH, EXACT_MATCH } ; + +// --------------------------------------------------------------------------- +// ^FUNCTION: kwdmatch - match a keyword +// +// ^SYNOPSIS: +// static kwdmatch_t kwdmatch(src, attempt, len) +// +// ^PARAMETERS: +// char * src -- the actual keyword to match +// char * attempt -- the possible keyword to compare against "src" +// int len -- number of character of "attempt" to consider +// (if 0 then we should use all of "attempt") +// +// ^DESCRIPTION: +// See if "attempt" matches some prefix of "src" (case insensitive). +// +// ^REQUIREMENTS: +// - attempt should be non-NULL and non-empty +// +// ^SIDE-EFFECTS: +// None. +// +// ^RETURN-VALUE: +// An enumeration value of type kwdmatch_t corresponding to whether +// We had an exact match, a partial match, or no match. +// +// ^ALGORITHM: +// Trivial +// ^^------------------------------------------------------------------------- +static kwdmatch_t +kwdmatch(const char * src, const char * attempt, int len =0) { + int i; + + if (src == attempt) return EXACT_MATCH ; + if ((src == NULL) || (attempt == NULL)) return NO_MATCH ; + if ((! *src) && (! *attempt)) return EXACT_MATCH ; + if ((! *src) || (! *attempt)) return NO_MATCH ; + + for (i = 0 ; ((i < len) || (len == 0)) && + (attempt[i]) && (attempt[i] != ' ') ; i++) { + if (TOLOWER(src[i]) != TOLOWER(attempt[i])) return NO_MATCH ; + } + + return (src[i]) ? PARTIAL_MATCH : EXACT_MATCH ; +} + +// **************************************************************** OptionSpec + + // Class that represents an option-specification + // *NOTE*:: Assumes that the char-ptr given to the constructor points + // to storage that will NOT be modified and whose lifetime will + // be as least as long as the OptionSpec object we construct. + // +class OptionSpec { +public: + OptionSpec(const char * decl =NULLSTR) + : hidden(0), spec(decl) + { + if (spec == NULL) spec = NULL_spec; + CheckHidden(); + } + + OptionSpec(const OptionSpec & cp) : hidden(cp.hidden), spec(cp.spec) {} + + // NOTE: use default destructor! + + // Assign to another OptionSpec + OptionSpec & + operator=(const OptionSpec & cp) { + if (this != &cp) { + spec = cp.spec; + hidden = cp.hidden; + } + return *this; + } + + // Assign to a string + OptionSpec & + operator=(const char * decl) { + if (spec != decl) { + spec = decl; + hidden = 0; + CheckHidden(); + } + return *this; + } + + // Convert to char-ptr by returning the original declaration-string + operator const char*() { return isHiddenOpt() ? (spec - 1) : spec; } + + // Is this option NULL? + int + isNULL(void) const { return ((spec == NULL) || (spec == NULL_spec)); } + + // Is this options incorrectly specified? + int + isSyntaxError(const char * name) const; + + // See if this is a Hidden option + int + isHiddenOpt(void) const { return hidden; } + + // Get the corresponding option-character + char + OptChar(void) const { return *spec; } + + // Get the corresponding long-option string + const char * + LongOpt(void) const { + return (spec[1] && spec[2] && (! isspace(spec[2]))) ? (spec + 2) : NULLSTR; + } + + // Does this option require an argument? + int + isValRequired(void) const { + return ((spec[1] == ':') || (spec[1] == '+')); + } + + // Does this option take an optional argument? + int + isValOptional(void) const { + return ((spec[1] == '?') || (spec[1] == '*')); + } + + // Does this option take no arguments? + int + isNoArg(void) const { + return ((spec[1] == '|') || (! spec[1])); + } + + // Can this option take more than one argument? + int + isList(void) const { + return ((spec[1] == '+') || (spec[1] == '*')); + } + + // Does this option take any arguments? + int + isValTaken(void) const { + return (isValRequired() || isValOptional()) ; + } + + // Format this option in the given buffer + unsigned + Format(char * buf, unsigned optctrls) const; + +private: + void + CheckHidden(void) { + if ((! hidden) && (*spec == '-')) { + ++hidden; + ++spec; + } + } + + unsigned hidden : 1; // hidden-flag + const char * spec; // string specification + + static const char NULL_spec[]; +} ; + +const char OptionSpec::NULL_spec[] = "\0\0\0" ; + +int +OptionSpec::isSyntaxError(const char * name) const { + int error = 0; + if ((! spec) || (! *spec)) { + cerr << name << ": empty option specifier." << endl; + cerr << "\tmust be at least 1 character long." << endl; + ++error; + } else if (spec[1] && (strchr("|?:*+", spec[1]) == NULL)) { + cerr << name << ": bad option specifier \"" << spec << "\"." << endl; + cerr << "\t2nd character must be in the set \"|?:*+\"." << endl; + ++error; + } + return error; +} + +// --------------------------------------------------------------------------- +// ^FUNCTION: OptionSpec::Format - format an option-spec for a usage message +// +// ^SYNOPSIS: +// unsigned OptionSpec::Format(buf, optctrls) const +// +// ^PARAMETERS: +// char * buf -- where to print the formatted option +// unsigned optctrls -- option-parsing configuration flags +// +// ^DESCRIPTION: +// Self-explanatory. +// +// ^REQUIREMENTS: +// - buf must be large enough to hold the result +// +// ^SIDE-EFFECTS: +// - writes to buf. +// +// ^RETURN-VALUE: +// Number of characters written to buf. +// +// ^ALGORITHM: +// Follow along in the source - it's not hard but it is tedious! +// ^^------------------------------------------------------------------------- +unsigned +OptionSpec::Format(char * buf, unsigned optctrls) const { +#ifdef NO_USAGE + return (*buf = '\0'); +#else + static char default_value[] = ""; + if (isHiddenOpt()) return (unsigned)(*buf = '\0'); + char optchar = OptChar(); + const char * longopt = LongOpt(); + char * p = buf ; + + const char * value = NULLSTR; + int longopt_len = 0; + int value_len = 0; + + if (longopt) { + value = ::strchr(longopt, ' '); + longopt_len = (value) ? (value - longopt) : ::strlen(longopt); + } else { + value = ::strchr(spec + 1, ' '); + } + while (value && (*value == ' ')) ++value; + if (value && *value) { + value_len = ::strlen(value); + } else { + value = default_value; + value_len = sizeof(default_value) - 1; + } + + if ((optctrls & Options::SHORT_ONLY) && + ((! isNullOpt(optchar)) || (optctrls & Options::NOGUESSING))) { + longopt = NULLSTR; + } + if ((optctrls & Options::LONG_ONLY) && + (longopt || (optctrls & Options::NOGUESSING))) { + optchar = '\0'; + } + if (isNullOpt(optchar) && (longopt == NULL)) { + *buf = '\0'; + return 0; + } + + *(p++) = '['; + + // print the single character option + if (! isNullOpt(optchar)) { + *(p++) = '-'; + *(p++) = optchar; + } + + if ((! isNullOpt(optchar)) && (longopt)) *(p++) = '|'; + + // print the long option + if (longopt) { + *(p++) = '-'; + if (! (optctrls & (Options::LONG_ONLY | Options::SHORT_ONLY))) { + *(p++) = '-'; + } + strncpy(p, longopt, longopt_len); + p += longopt_len; + } + + // print any argument the option takes + if (isValTaken()) { + *(p++) = ' ' ; + if (isValOptional()) *(p++) = '[' ; + strcpy(p, value); + p += value_len; + if (isList()) { + strcpy(p, " ..."); + p += 4; + } + if (isValOptional()) *(p++) = ']' ; + } + + *(p++) = ']'; + *p = '\0'; + + return (unsigned) strlen(buf); +#endif /* USE_STDIO */ +} + +// ******************************************************************* Options + +#if (defined(MSWIN) || defined(OS2) || defined(MSDOS)) +# define DIR_SEP_CHAR '\\' +#else +# define DIR_SEP_CHAR '/' +#endif + +Options::Options(const char * name, const char * const optv[]) + : cmdname(name), optvec(optv), explicit_end(0), optctrls(DEFAULT), + nextchar(NULLSTR), listopt(NULLSTR) +{ + const char * basename = ::strrchr(cmdname, DIR_SEP_CHAR); + if (basename) cmdname = basename + 1; + check_syntax(); +} + +Options::~Options(void) {} + + // Make sure each option-specifier has correct syntax. + // + // If there is even one invalid specifier, then exit ungracefully! + // +void +Options::check_syntax(void) const { + int errors = 0; + if ((optvec == NULL) || (! *optvec)) return; + + for (const char * const * optv = optvec ; *optv ; optv++) { + OptionSpec optspec = *optv; + errors += optspec.isSyntaxError(cmdname); + } + if (errors) exit(127); +} + +// --------------------------------------------------------------------------- +// ^FUNCTION: Options::match_opt - match an option +// +// ^SYNOPSIS: +// const char * match_opt(opt, int ignore_case) const +// +// ^PARAMETERS: +// char opt -- the option-character to match +// int ignore_case -- should we ignore character-case? +// +// ^DESCRIPTION: +// See if "opt" is found in "optvec" +// +// ^REQUIREMENTS: +// - optvec should be non-NULL and terminated by a NULL pointer. +// +// ^SIDE-EFFECTS: +// None. +// +// ^RETURN-VALUE: +// NULL if no match is found, +// otherwise a pointer to the matching option-spec. +// +// ^ALGORITHM: +// foreach option-spec +// - see if "opt" is a match, if so return option-spec +// end-for +// ^^------------------------------------------------------------------------- +const char * +Options::match_opt(char opt, int ignore_case) const { + if ((optvec == NULL) || (! *optvec)) return NULLSTR; + + for (const char * const * optv = optvec ; *optv ; optv++) { + OptionSpec optspec = *optv; + char optchar = optspec.OptChar(); + if (isNullOpt(optchar)) continue; + if (opt == optchar) { + return optspec; + } else if (ignore_case && (TOLOWER(opt) == TOLOWER(optchar))) { + return optspec; + } + } + + return NULLSTR; // not found +} + +// --------------------------------------------------------------------------- +// ^FUNCTION: Options::match_longopt - match a long-option +// +// ^SYNOPSIS: +// const char * Options::match_longopt(opt, len, ambiguous) +// +// ^PARAMETERS: +// char * opt -- the long-option to match +// int len -- the number of character of "opt" to match +// int & ambiguous -- set by this routine before returning. +// +// ^DESCRIPTION: +// Try to match "opt" against some unique prefix of a long-option +// (case insensitive). +// +// ^REQUIREMENTS: +// - optvec should be non-NULL and terminated by a NULL pointer. +// +// ^SIDE-EFFECTS: +// - *ambiguous is set to '1' if "opt" matches >1 long-option +// (otherwise it is set to 0). +// +// ^RETURN-VALUE: +// NULL if no match is found, +// otherwise a pointer to the matching option-spec. +// +// ^ALGORITHM: +// ambiguous is FALSE +// foreach option-spec +// if we have an EXACT-MATCH, return the option-spec +// if we have a partial-match then +// if we already had a previous partial match then +// set ambiguous = TRUE and return NULL +// else +// remember this options spec and continue matching +// end-if +// end-if +// end-for +// if we had exactly 1 partial match return it, else return NULL +// ^^------------------------------------------------------------------------- +const char * +Options::match_longopt(const char * opt, int len, int & ambiguous) const { + kwdmatch_t result; + const char * matched = NULLSTR ; + + ambiguous = 0; + if ((optvec == NULL) || (! *optvec)) return NULLSTR; + + for (const char * const * optv = optvec ; *optv ; optv++) { + OptionSpec optspec = *optv; + const char * longopt = optspec.LongOpt(); + if (longopt == NULL) continue; + result = kwdmatch(longopt, opt, len); + if (result == EXACT_MATCH) { + return optspec; + } else if (result == PARTIAL_MATCH) { + if (matched) { + ++ambiguous; + return NULLSTR; + } else { + matched = optspec; + } + } + }//for + + return matched; +} + +// --------------------------------------------------------------------------- +// ^FUNCTION: Options::parse_opt - parse an option +// +// ^SYNOPSIS: +// int Options::parse_opt(iter, optarg) +// +// ^PARAMETERS: +// OptIter & iter -- option iterator +// const char * & optarg -- where to store any option-argument +// +// ^DESCRIPTION: +// Parse the next option in iter (advancing as necessary). +// Make sure we update the nextchar pointer along the way. Any option +// we find should be returned and optarg should point to its argument. +// +// ^REQUIREMENTS: +// - nextchar must point to the prospective option character +// +// ^SIDE-EFFECTS: +// - iter is advanced when an argument completely parsed +// - optarg is modified to point to any option argument +// - if Options::QUIET is not set, error messages are printed on cerr +// +// ^RETURN-VALUE: +// 'c' if the -c option was matched (optarg points to its argument) +// BADCHAR if the option is invalid (optarg points to the bad +// option-character). +// +// ^ALGORITHM: +// It gets complicated -- follow the comments in the source. +// ^^------------------------------------------------------------------------- +int +Options::parse_opt(OptIter & iter, const char * & optarg) { + listopt = NULLSTR; // reset the list pointer + + if ((optvec == NULL) || (! *optvec)) return Options::ENDOPTS; + + // Try to match a known option + OptionSpec optspec = match_opt(*(nextchar++), (optctrls & Options::ANYCASE)); + + // Check for an unknown option + if (optspec.isNULL()) { + // See if this was a long-option in disguise + if (! (optctrls & Options::NOGUESSING)) { + unsigned save_ctrls = optctrls; + const char * save_nextchar = nextchar; + nextchar -= 1; + optctrls |= (Options::QUIET | Options::NOGUESSING); + int optchar = parse_longopt(iter, optarg); + optctrls = save_ctrls; + if (optchar > 0) { + return optchar; + } else { + nextchar = save_nextchar; + } + } + if (! (optctrls & Options::QUIET)) { + cerr << cmdname << ": unknown option -" + << *(nextchar - 1) << "." << endl ; + } + optarg = (nextchar - 1); // record the bad option in optarg + return Options::BADCHAR; + } + + // If no argument is taken, then leave now + if (optspec.isNoArg()) { + optarg = NULLSTR; + return optspec.OptChar(); + } + + // Check for argument in this arg + if (*nextchar) { + optarg = nextchar; // the argument is right here + nextchar = NULLSTR; // we've exhausted this arg + if (optspec.isList()) listopt = optspec ; // save the list-spec + return optspec.OptChar(); + } + + // Check for argument in next arg + const char * nextarg = iter.curr(); + if ((nextarg != NULL) && + (optspec.isValRequired() || (! isOption(optctrls, nextarg)))) { + optarg = nextarg; // the argument is here + iter.next(); // end of arg - advance + if (optspec.isList()) listopt = optspec ; // save the list-spec + return optspec.OptChar(); + } + + // No argument given - if its required, thats an error + optarg = NULLSTR; + if (optspec.isValRequired() && !(optctrls & Options::QUIET)) { + cerr << cmdname << ": argument required for -" << optspec.OptChar() + << " option." << endl ; + } + return optspec.OptChar(); +} + +// --------------------------------------------------------------------------- +// ^FUNCTION: Options::parse_longopt - parse a long-option +// +// ^SYNOPSIS: +// int Options::parse_longopt(iter, optarg) +// +// ^PARAMETERS: +// OptIter & iter -- option iterator +// const char * & optarg -- where to store any option-argument +// +// ^DESCRIPTION: +// Parse the next long-option in iter (advancing as necessary). +// Make sure we update the nextchar pointer along the way. Any option +// we find should be returned and optarg should point to its argument. +// +// ^REQUIREMENTS: +// - nextchar must point to the prospective option character +// +// ^SIDE-EFFECTS: +// - iter is advanced when an argument completely parsed +// - optarg is modified to point to any option argument +// - if Options::QUIET is not set, error messages are printed on cerr +// +// ^RETURN-VALUE: +// 'c' if the the long-option corresponding to the -c option was matched +// (optarg points to its argument) +// BADKWD if the option is invalid (optarg points to the bad long-option +// name). +// +// ^ALGORITHM: +// It gets complicated -- follow the comments in the source. +// ^^------------------------------------------------------------------------- +int +Options::parse_longopt(OptIter & iter, const char * & optarg) { + int len = 0, ambiguous = 0; + + listopt = NULLSTR ; // reset the list-spec + + if ((optvec == NULL) || (! *optvec)) return Options::ENDOPTS; + + // if a value is supplied in this argv element, get it now + const char * val = strpbrk(nextchar, ":=") ; + if (val) { + len = val - nextchar ; + ++val ; + } + + // Try to match a known long-option + OptionSpec optspec = match_longopt(nextchar, len, ambiguous); + + // Check for an unknown long-option + if (optspec.isNULL()) { + // See if this was a short-option in disguise + if ((! ambiguous) && (! (optctrls & Options::NOGUESSING))) { + unsigned save_ctrls = optctrls; + const char * save_nextchar = nextchar; + optctrls |= (Options::QUIET | Options::NOGUESSING); + int optchar = parse_opt(iter, optarg); + optctrls = save_ctrls; + if (optchar > 0) { + return optchar; + } else { + nextchar = save_nextchar; + } + } + if (! (optctrls & Options::QUIET)) { + cerr << cmdname << ": " << ((ambiguous) ? "ambiguous" : "unknown") + << " option " + << ((optctrls & Options::LONG_ONLY) ? "-" : "--") + << nextchar << "." << endl ; + } + optarg = nextchar; // record the bad option in optarg + nextchar = NULLSTR; // we've exhausted this argument + return (ambiguous) ? Options::AMBIGUOUS : Options::BADKWD; + } + + // If no argument is taken, then leave now + if (optspec.isNoArg()) { + if ((val) && ! (optctrls & Options::QUIET)) { + cerr << cmdname << ": option " + << ((optctrls & Options::LONG_ONLY) ? "-" : "--") + << optspec.LongOpt() << " does NOT take an argument." << endl ; + } + optarg = val; // record the unexpected argument + nextchar = NULLSTR; // we've exhausted this argument + return optspec.OptChar(); + } + + // Check for argument in this arg + if (val) { + optarg = val; // the argument is right here + nextchar = NULLSTR; // we exhausted the rest of this arg + if (optspec.isList()) listopt = optspec ; // save the list-spec + return optspec.OptChar(); + } + + // Check for argument in next arg + const char * nextarg = iter.curr(); // find the next argument to parse + if ((nextarg != NULL) && + (optspec.isValRequired() || (! isOption(optctrls, nextarg)))) { + optarg = nextarg; // the argument is right here + iter.next(); // end of arg - advance + nextchar = NULLSTR; // we exhausted the rest of this arg + if (optspec.isList()) listopt = optspec ; // save the list-spec + return optspec.OptChar(); + } + + // No argument given - if its required, thats an error + optarg = NULLSTR; + if (optspec.isValRequired() && !(optctrls & Options::QUIET)) { + const char * longopt = optspec.LongOpt(); + const char * spc = ::strchr(longopt, ' '); + int longopt_len; + if (spc) { + longopt_len = spc - longopt; + } else { + longopt_len = ::strlen(longopt); + } + cerr << cmdname << ": argument required for " + << ((optctrls & Options::LONG_ONLY) ? "-" : "--"); + cerr.write(longopt, longopt_len) << " option." << endl ; + } + nextchar = NULLSTR; // we exhausted the rest of this arg + return optspec.OptChar(); +} + +// --------------------------------------------------------------------------- +// ^FUNCTION: Options::usage - print usage +// +// ^SYNOPSIS: +// void Options::usage(os, positionals) +// +// ^PARAMETERS: +// ostream & os -- where to print the usage +// char * positionals -- command-line syntax for any positional args +// +// ^DESCRIPTION: +// Print command-usage (using either option or long-option syntax) on os. +// +// ^REQUIREMENTS: +// os should correspond to an open output file. +// +// ^SIDE-EFFECTS: +// Prints on os +// +// ^RETURN-VALUE: +// None. +// +// ^ALGORITHM: +// Print usage on os, wrapping long lines where necessary. +// ^^------------------------------------------------------------------------- +void +Options::usage(ostream & os, const char * positionals) const { +#ifdef NO_USAGE + return; +#else + const char * const * optv = optvec; + unsigned cols = 79; + int first, nloop; + char buf[256] ; + + if ((optv == NULL) || (! *optv)) return; + + // print first portion "usage: progname" + os << "usage: " << cmdname ; + unsigned ll = strlen(cmdname) + 7; + + // save the current length so we know how much space to skip for + // subsequent lines. + // + unsigned margin = ll + 1; + + // print the options and the positional arguments + for (nloop = 0, first = 1 ; !nloop ; optv++, first = 0) { + unsigned len; + OptionSpec optspec = *optv; + + // figure out how wide this parameter is (for printing) + if (! *optv) { + len = strlen(positionals); + ++nloop; // terminate this loop + } else { + if (optspec.isHiddenOpt()) continue; + len = optspec.Format(buf, optctrls); + } + + // Will this fit? + if ((ll + len + 1) > (cols - first)) { + os << '\n' ; // No - start a new line; +#ifdef USE_STDIO + for (int _i_ = 0; _i_ < margin; ++_i_) os << " "; +#else + os.width(margin); os << "" ; +#endif + ll = margin; + } else { + os << ' ' ; // Yes - just throw in a space + ++ll; + } + ll += len; + os << ((nloop) ? positionals : buf) ; + }// for each ad + + os << endl ; +#endif /* NO_USAGE */ +} + + +// --------------------------------------------------------------------------- +// ^FUNCTION: Options::operator() - get options from the command-line +// +// ^SYNOPSIS: +// int Options::operator()(iter, optarg) +// +// ^PARAMETERS: +// OptIter & iter -- option iterator +// const char * & optarg -- where to store any option-argument +// +// ^DESCRIPTION: +// Parse the next option in iter (advancing as necessary). +// Make sure we update the nextchar pointer along the way. Any option +// we find should be returned and optarg should point to its argument. +// +// ^REQUIREMENTS: +// None. +// +// ^SIDE-EFFECTS: +// - iter is advanced when an argument is completely parsed +// - optarg is modified to point to any option argument +// - if Options::QUIET is not set, error messages are printed on cerr +// +// ^RETURN-VALUE: +// 0 if all options have been parsed. +// 'c' if the the option or long-option corresponding to the -c option was +// matched (optarg points to its argument). +// BADCHAR if the option is invalid (optarg points to the bad option char). +// BADKWD if the option is invalid (optarg points to the bad long-opt name). +// AMBIGUOUS if an ambiguous keyword name was given (optarg points to the +// ambiguous keyword name). +// POSITIONAL if PARSE_POS was set and the current argument is a positional +// parameter (in which case optarg points to the positional argument). +// +// ^ALGORITHM: +// It gets complicated -- follow the comments in the source. +// ^^------------------------------------------------------------------------- +int +Options::operator()(OptIter & iter, const char * & optarg) { + int parse_opts_only = isOptsOnly(optctrls); + if (parse_opts_only) explicit_end = 0; + + // See if we have an option left over from before ... + if ((nextchar) && *nextchar) { + return parse_opt(iter, optarg); + } + + // Check for end-of-options ... + const char * arg = NULLSTR; + int get_next_arg = 0; + do { + arg = iter.curr(); + get_next_arg = 0; + if (arg == NULL) { + listopt = NULLSTR; + return Options::ENDOPTS; + } else if ((! explicit_end) && isEndOpts(arg)) { + iter.next(); // advance past end-of-options arg + listopt = NULLSTR; + explicit_end = 1; + if (parse_opts_only) return Options::ENDOPTS; + get_next_arg = 1; // make sure we look at the next argument. + } + } while (get_next_arg); + + // Do we have a positional arg? + if ( explicit_end || (! isOption(optctrls, arg)) ) { + if (parse_opts_only) { + return Options::ENDOPTS; + } else { + optarg = arg; // set optarg to the positional argument + iter.next(); // advance iterator to the next argument + return Options::POSITIONAL; + } + } + + iter.next(); // pass the argument that arg already points to + + // See if we have a long option ... + if (! (optctrls & Options::SHORT_ONLY)) { + if ((*arg == '-') && (arg[1] == '-')) { + nextchar = arg + 2; + return parse_longopt(iter, optarg); + } else if ((optctrls & Options::PLUS) && (*arg == '+')) { + nextchar = arg + 1; + return parse_longopt(iter, optarg); + } + } + if (*arg == '-') { + nextchar = arg + 1; + if (optctrls & Options::LONG_ONLY) { + return parse_longopt(iter, optarg); + } else { + return parse_opt(iter, optarg); + } + } + + // If we get here - it is because we have a list value + OptionSpec optspec = listopt; + optarg = arg ; // record the list value + return optspec.OptChar() ; +} + diff --git a/Ridges_3/examples/Ridges_3/options.h b/Ridges_3/examples/Ridges_3/options.h new file mode 100644 index 00000000000..28d34501388 --- /dev/null +++ b/Ridges_3/examples/Ridges_3/options.h @@ -0,0 +1,492 @@ +// **************************************************************************** +// ^FILE: options.h - option parsing classes +// +// ^DESCRIPTION: +// This file defines classes used to parse command-line options. +// Options may be parsed from an array of strings, or from any structure +// for which a corresponding option-iterator exists. +// +// ^HISTORY: +// 03/06/92 Brad Appleton Created +// +// 03/23/93 Brad Appleton +// - Added OptIstreamIter class +// +// 03/08/94 Brad Appleton +// - Added Options::reset() member function +// +// 07/31/97 Brad Appleton +// - Added PARSE_POS control flag and POSITIONAL return value +// ^^************************************************************************** + +#ifndef _options_h +#define _options_h + +using namespace std; +#include + + +// Abstract class to iterate through options/arguments +// +class OptIter { +public: + OptIter(void) {} + virtual ~OptIter(void); + + // curr() returns the current item in the iterator without + // advancing on to the next item. If we are at the end of items + // then NULL is returned. + virtual const char * + curr(void) = 0; + + // next() advances to the next item. + virtual void + next(void) = 0; + + // operator() returns the current item in the iterator and then + // advances on to the next item. If we are at the end of items + // then NULL is returned. + virtual const char * + operator()(void); +} ; + +// Abstract class for a rewindable OptIter +// +class OptIterRwd : public OptIter { +public: + OptIterRwd(void); + + virtual ~OptIterRwd(void); + + virtual const char * + curr(void) = 0; + + virtual void + next(void) = 0; + + virtual const char * + operator()(void) = 0; + + // rewind() resets the "current-element" to the first one in the "list" + virtual void + rewind(void) = 0; +} ; + +// Class to iterate through an array of tokens. The array may be terminated +// by NULL or a count containing the number of tokens may be given. +// +class OptArgvIter : public OptIterRwd { +private: + const char * const * av; // arg vector + int ac; // arg count + int ndx; // index of current arg + +public: + OptArgvIter(const char * const argv[]) + : av(argv), ac(-1), ndx(0) {} + + OptArgvIter(int argc, const char * const argv[]) + : av(argv), ac(argc), ndx(0) {init(argc, argv); } + + void init(int argc, const char * const argv[]);// {v=argv; ac=argc; ndx=0;} + + virtual + ~OptArgvIter(void); + + virtual const char * + curr(void); + + virtual void + next(void); + + virtual const char * + operator()(void); + + virtual void + rewind(void); + + // index returns the current index to use for argv[] + int index(void) { return ndx; } +} ; + + +// Class to iterate through a string containing delimiter-separated tokens +// +class OptStrTokIter : public OptIterRwd { +private: + unsigned len; // length of token-string + const char * str; // the token-string + const char * seps; // delimiter-set (separator-characters) + const char * cur; // current token + char * tokstr; // our copy of the token-string + + static const char * default_delims; // default delimiters = whitespace + +public: + OptStrTokIter(const char * tokens, const char * delimiters =0); + + virtual + ~OptStrTokIter(void); + + virtual const char * + curr(void); + + virtual void + next(void); + + virtual const char * + operator()(void); + + virtual void + rewind(void); + + // delimiters() with NO arguments returns the current set of delimiters, + // If an argument is given then it is used as the new set of delimiters. + const char * + delimiters(void) { return seps; } + + void + delimiters(const char * delims) { + seps = (delims) ? delims : default_delims ; + } +} ; + + + // OptIstreamIter is a class for iterating over arguments that come + // from an input stream. Each line of the input stream is considered + // to be a set of white-space separated tokens. If the the first + // non-white character on a line is '#' ('!' for VMS systems) then + // the line is considered a comment and is ignored. + // + // *Note:: If a line is more than 1022 characters in length then we + // treat it as if it were several lines of length 1022 or less. + // + // *Note:: The string tokens returned by this iterator are pointers + // to temporary buffers which may not necessarily stick around + // for too long after the call to curr() or operator(), hence + // if you need the string value to persist - you will need to + // make a copy. + // +class OptIstreamIter : public OptIter { +private: + istream & is ; + OptStrTokIter * tok_iter ; + + void + fill(void); + +public: + static const unsigned MAX_LINE_LEN ; + + OptIstreamIter(istream & input); + + virtual + ~OptIstreamIter(void); + + virtual const char * + curr(void); + + virtual void + next(void); + + virtual const char * + operator()(void); +} ; + + +// Now we are ready to define a class to declare and parse command-options +// +// +// CLASS +// ===== +// Options -- parse command-line options +// +// SYNOPSIS +// ======== +// #include +// +// Options opts(cmdname, optv); +// char cmdname[], *optv[]; +// +// DESCRIPTION +// =========== +// The Options constructor expects a command-name (usually argv[0]) and +// a pointer to an array of strings. The last element in this array MUST +// be NULL. Each non-NULL string in the array must have the following format: +// +// The 1st character must be the option-name ('c' for a -c option). +// +// The 2nd character must be one of '|', '?', ':', '*', or '+'. +// '|' -- indicates that the option takes NO argument; +// '?' -- indicates that the option takes an OPTIONAL argument; +// ':' -- indicates that the option takes a REQUIRED argument; +// '*' -- indicates that the option takes 0 or more arguments; +// '+' -- indicates that the option takes 1 or more arguments; +// +// The remainder of the string must be the long-option name. +// +// If desired, the long-option name may be followed by one or more +// spaces and then by the name of the option value. This name will +// be used when printing usage messages. If the option-value-name +// is not given then the string "" will be used in usage +// messages. +// +// One may use a space to indicate that a particular option does not +// have a corresponding long-option. For example, "c: " (or "c:") +// means the -c option takes a value & has NO corresponding long-option. +// +// To specify a long-option that has no corresponding single-character +// option is a bit trickier: Options::operator() still needs an "option- +// character" to return when that option is matched. One may use a whitespace +// character or a non-printable character as the single-character option +// in such a case. (hence " |hello" would only match "--hello"). +// +// EXCEPTIONS TO THE ABOVE: +// ------------------------ +// If the 1st character of the string is '-', then the rest of the +// string must correspond to the above format, and the option is +// considered to be a hidden-option. This means it will be parsed +// when actually matching options from the command-line, but will +// NOT show-up if a usage message is printed using the usage() member +// function. Such an example might be "-h|hidden". If you want to +// use any "dummy" options (options that are not parsed, but that +// to show up in the usage message), you can specify them along with +// any positional parameters to the usage() member function. +// +// If the 2nd character of the string is '\0' then it is assumed +// that there is no corresponding long-option and that the option +// takes no argument (hence "f", and "f| " are equivalent). +// +// Examples: +// const char * optv[] = { +// "c:count ", +// "s?str ", +// "x", +// " |hello", +// "g+groups ", +// NULL +// } ; +// +// optv[] now corresponds to the following: +// +// usage: cmdname [-c|--count ] [-s|--str []] +// [-x] [--hello] [-g|--groups ...] +// +// Long-option names are matched case-insensitive and only a unique prefix +// of the name needs to be specified. +// +// Option-name characters are case-sensitive! +// +// CAVEAT +// ====== +// Because of the way in which multi-valued options and options with optional +// values are handled, it is NOT possible to supply a value to an option in +// a separate argument (different argv[] element) if the value is OPTIONAL +// and begins with a '-'. What this means is that if an option "-s" takes an +// optional value value and you wish to supply a value of "-foo" then you must +// specify this on the command-line as "-s-foo" instead of "-s -foo" because +// "-s -foo" will be considered to be two separate sets of options. +// +// A multi-valued option is terminated by another option or by the end-of +// options. The following are all equivalent (if "-l" is a multi-valued +// option and "-x" is an option that takes no value): +// +// cmdname -x -l item1 item2 item3 -- arg1 arg2 arg3 +// cmdname -x -litem1 -litem2 -litem3 -- arg1 arg2 arg3 +// cmdname -l item1 item2 item3 -x arg1 arg2 arg3 +// +// +// EXAMPLE +// ======= +// #include +// +// static const char * optv[] = { +// "H|help", +// "c:count ", +// "s?str ", +// "x", +// " |hello", +// "g+groups ", +// NULL +// } ; +// +// main(int argc, char * argv[]) { +// int optchar; +// const char * optarg; +// const char * str = "default_string"; +// int count = 0, xflag = 0, hello = 0; +// int errors = 0, ngroups = 0; +// +// Options opts(*argv, optv); +// OptArgvIter iter(--argc, ++argv); +// +// while( optchar = opts(iter, optarg) ) { +// switch (optchar) { +// case 'H' : +// opts.usage(cout, "files ..."); +// exit(0); +// break; +// case 'g' : +// ++ngroups; break; // the groupname is in "optarg" +// case 's' : +// str = optarg; break; +// case 'x' : +// ++xflag; break; +// case ' ' : +// ++hello; break; +// case 'c' : +// if (optarg == NULL) ++errors; +// else count = (int) atol(optarg); +// break; +// default : ++errors; break; +// } //switch +// } +// +// if (errors || (iter.index() == argc)) { +// if (! errors) { +// cerr << opts.name() << ": no filenames given." << endl ; +// } +// opts.usage(cerr, "files ..."); +// exit(1); +// } +// +// cout << "xflag=" << ((xflag) ? "ON" : "OFF") << endl +// << "hello=" << ((hello) ? "YES" : "NO") << endl +// << "count=" << count << endl +// << "str=\"" << ((str) ? str : "No value given!") << "\"" << endl +// << "ngroups=" << ngroups << endl ; +// +// if (iter.index() < argc) { +// cout << "files=" ; +// for (int i = iter.index() ; i < argc ; i++) { +// cout << "\"" << argv[i] << "\" " ; +// } +// cout << endl ; +// } +// } +// +class Options { +private: + const char * cmdname; // name of the command + const char * const * optvec; // vector of option-specifications (last=NULL) + unsigned explicit_end : 1; // were we terminated because of "--"? + unsigned optctrls : 7; // control settings (a set of OptCtrl masks) + const char * nextchar; // next option-character to process + const char * listopt; // last list-option we matched + + void + check_syntax(void) const; + + const char * + match_opt(char opt, int ignore_case =0) const; + + const char * + match_longopt(const char * opt, int len, int & ambiguous) const; + + int + parse_opt(OptIter & iter, const char * & optarg); + + int + parse_longopt(OptIter & iter, const char * & optarg); + +public: + enum OptCtrl { + DEFAULT = 0x00, // Default setting + ANYCASE = 0x01, // Ignore case when matching short-options + QUIET = 0x02, // Dont print error messages + PLUS = 0x04, // Allow "+" as a long-option prefix + SHORT_ONLY = 0x08, // Dont accept long-options + LONG_ONLY = 0x10, // Dont accept short-options + // (also allows "-" as a long-option prefix). + NOGUESSING = 0x20, // Normally, when we see a short (long) option + // on the command line that doesnt match any + // known short (long) options, then we try to + // "guess" by seeing if it will match any known + // long (short) option. Setting this mask prevents + // this "guessing" from occurring. + PARSE_POS = 0x40 // By default, Options will not present positional + // command-line arguments to the user and will + // instead stop parsing when the first positonal + // argument has been encountered. If this flag + // is given, Options will present positional + // arguments to the user with a return code of + // POSITIONAL; ENDOPTS will be returned only + // when the end of the argument list is reached. + } ; + + // Error return values for operator() + // + enum OptRC { + ENDOPTS = 0, + BADCHAR = -1, + BADKWD = -2, + AMBIGUOUS = -3, + POSITIONAL = -4 + } ; + + Options(const char * name, const char * const optv[]); + + virtual + ~Options(void); + + // name() returns the command name + const char * + name(void) const { return cmdname; } + + // ctrls() (with no arguments) returns the existing control settings + unsigned + ctrls(void) const { return optctrls; } + + // ctrls() (with 1 argument) sets new control settings + void + ctrls(unsigned newctrls) { optctrls = newctrls; } + + // reset for another pass to parse for options + void + reset(void) { nextchar = listopt = NULL; } + + // usage() prints options usage (followed by any positional arguments + // listed in the parameter "positionals") on the given outstream + void + usage(ostream & os, const char * positionals) const ; + + // operator() iterates through the arguments as necessary (using the + // given iterator) and returns the character value of the option + // (or long-option) that it matched. If the option has a value + // then the value given may be found in optarg (otherwise optarg + // will be NULL). + // + // 0 is returned upon end-of-options. At this point, "iter" may + // be used to process any remaining positional parameters. If the + // PARSE_POS control-flag is set then 0 is returned only when all + // arguments in "iter" have been exhausted. + // + // If an invalid option is found then BADCHAR is returned and *optarg + // is the unrecognized option character. + // + // If an invalid long-option is found then BADKWD is returned and optarg + // points to the bad long-option. + // + // If an ambiguous long-option is found then AMBIGUOUS is returned and + // optarg points to the ambiguous long-option. + // + // If the PARSE_POS control-flag is set then POSITIONAL is returned + // when a positional argument is encountered and optarg points to + // the positonal argument (and "iter" is advanced to the next argument + // in the iterator). + // + // Unless Options::QUIET is used, missing option-arguments and + // invalid options (and the like) will automatically cause error + // messages to be issued to cerr. + int + operator()(OptIter & iter, const char * & optarg) ; + + // Call this member function after operator() has returned 0 + // if you want to know whether or not options were explicitly + // terminated because "--" appeared on the command-line. + // + int + explicit_endopts() const { return explicit_end; } +} ; + +#endif /* _options_h */ diff --git a/Ridges_3/include/CGAL/algo.C b/Ridges_3/include/CGAL/algo.C deleted file mode 100644 index 935371fd659..00000000000 --- a/Ridges_3/include/CGAL/algo.C +++ /dev/null @@ -1,379 +0,0 @@ -#ifndef _RIDGE_3_H_ -#define _RIDGE_3_H_ - -#include - - -//Orient monge normal according to mesh normals. - -CGAL_BEGIN_NAMESPACE - - -enum Ridge_type {NONE=0, BLUE, RED, CREST, BE, BH, BC, RE, RH, RC}; - - -//--------------------------------------- -//Ridge_line : a connected sequence of edges crossed by a ridge, with type and weigths -//--------------------------------------- -template < class Poly > class Ridge_line -{ -public: - typedef typename Poly::Traits::FT FT; - typedef typename Poly::Vertex_handle Vertex_handle; - typedef typename Poly::Halfedge_handle Halfedge_handle; - typedef typename Poly::Facet_handle Facet_handle; - typedef std::pair ridge_he; - -protected: - Ridge_type line_type - std::list line; - FT strength; - FT sharpness; - -public: - const FT weight() const {return weight;} - const FT sharpness() const {return sharpness;} - std::list* line() { return &line;} - - //constructor - Ridge_line( Halfedge_handle h1, Halfedge_handle h2, Ridge_type r_type) : - line_type(r_type), strength(0.) - {}; - - //compute the barycentric coordinate of the xing point (blue or red) - //for he: p->q coord is st xing_point = coord*p + (1-coord)*q - FT bary_coord(Halfedge_handle he); - void compute_weight(char color); - void compute_sharpness(char color); - //When the line is extended with a he, the bary coord of the - //crossing point is computed, the pair (he,coord) is added and the - //weigths are updated - void addback( Halfedge_handle he); - void addfront( Halfedge_handle he); - -}; - -// IMPLEMENTATION OF Ridge_line members -////////////////////////////////////////////////////////////////////////////// - -//constructor -template < class Poly > -Ridge_line( Halfedge_handle h1, Halfedge_handle h2, Ridge_type r_type) : - line_type(r_type), strength(0.) -{ - line.push_back(h1); - addback(h2); -} - -template < class Poly > -void Ridge_line:: -addback( Halfedge_handle he) -{ - Halfedge_handle he_cur = ( --(line.end()) )->first; - FT coord = bary_coord(he); - FT coord_cur = bary_coord(he_cur); - Vertex_handle v_p = he->opposite()->vertex(), v_q = he->vertex(), - v_p_cur = he_cur->opposite()->vertex(), v_q_cur = he->vertex(); // he: p->q - FT k; - if ( (line_type == BE) || (line_type == BH) || (line_type == BC) ) { - k =( std::fabs(v_p->k1) * coord + std::fabs(v_q->k1) * (1-coord) )/2; - } - if ( (line_type == RE) || (line_type == RH) || (line_type == RC) ) { - k =( std::fabs(v_p->k2) * coord + std::fabs(v_q->k2) * (1-coord) )/2; - } - Vector_3 segment = (v_p->point()-ORIGIN)*coord + (v_q->point()-ORIGIN)*(1-coord) - - ((v_p_cur->point()-ORIGIN)*coord_cur + (v_q_cur->point()-ORIGIN)*(1-coord_cur)); - strength += k * CGAL::sqrt(segment * segment); - - line.push_back( pair(he, coord)); -} - -template < class Poly > -void Ridge_line:: -addfront( Halfedge_handle he) -{ - Halfedge_handle he_cur = ( line.begin() )->first; - FT coord = bary_coord(he); - FT coord_cur = bary_coord(he_cur); - Vertex_handle v_p = he->opposite()->vertex(), v_q = he->vertex(), - v_p_cur = he_cur->opposite()->vertex(), v_q_cur = he->vertex(); // he: p->q - FT k; - if ( (line_type == BE) || (line_type == BH) || (line_type == BC) ) { - k =( std::fabs(v_p->k1) * coord + std::fabs(v_q->k1) * (1-coord) )/2; - } - if ( (line_type == RE) || (line_type == RH) || (line_type == RC) ) { - k =( std::fabs(v_p->k2) * coord + std::fabs(v_q->k2) * (1-coord) )/2; - } - Vector_3 segment = (v_p->point()-ORIGIN)*coord + (v_q->point()-ORIGIN)*(1-coord) - - ((v_p_cur->point()-ORIGIN)*coord_cur + (v_q_cur->point()-ORIGIN)*(1-coord_cur)); - strength += k * CGAL::sqrt(segment * segment); - - line.push_front( pair(he, coord)); -} - - -template < class Poly > -FT Ridge_line:: -bary_coord(Halfedge_handle he) -{ - FT b_p, b_q; // extremalities at p and q for he: p->q - if ( (line_type == BE) || (line_type == BH) || (line_type == BC) ) { - b_p = he->opposite()->vertex()->b0(); - b_q = he->vertex()->b0(); - } - if ( (line_type == RE) || (line_type == RH) || (line_type == RC) ) { - b_p = he->opposite()->vertex()->b3(); - b_q = he->vertex()->b3(); - } - return std::fabs(b_q) / ( std::fabs(b_q) + std::fabs(b_p) ); -} - - -/////////Ridge_approximation////////////////////////////////////////// - -//Find BLUE ridges (Elliptic or Hyperbolic) -//do the same for RED and CREST -//iterate on P facets, find a non-visited, regular, 2BXing triangle, -//follow non-visited, regular, 2BXing triangles in both sens to create -//a Ridge line. -//Each time a edge is added the strength of the current line is updated -// + length(edge)*|k| -void -compute_ridges(Ridge_type r_type) -{ - //set all facets non visited - - - Facet_iterator itb = P.facets_begin(), ite = P.facets_end(); - for(;itb!=ite;itb++) - { - Facet_handle f= &(*itb); - if (f->is_visited()) continue; - f->set_visited(true); - Halfedge_handle h1, h2, curhe1, curhe2, curhe; - - //h1 h2 are the hedges crossed if any, r_type should be BLUE, - //RED or CREST ; cur_ridge_type should be BE, BH, BC, RE, RH, RC or NONE - Ridge_type cur_ridge_type = facet_ridge_type(f,h1,h2,r_type) - if ( cur_ridge_type == NONE ) continue; - - //When a he is inserted in a Ridge_line, a pair - //is created and the weight(s) are updated - Ridge_line *cur_ridge_line = new Ridge_line(h1,h2,cur_ridge_type); - *ridge_lines_it++ = cur_ridge_line; - - //next triangle adjacent to h1 (push_front) - if ( !(h1->is_border_edge()) ) - { - f = h1->opposite()->facet(); - curhe = h1; - while (cur_ridge_type == facet_ridge_type(f,curhe1,curhe2,r_type)) - { - //follow the ridge from curhe - if (f->is_visited()) break; - f->set_visited(true); - if (curhe->opposite() == curhe1) curhe = curhe2; - else curhe = curhe1; - cur_ridge_line->addfront(curhe); - if ( !(curhe->is_border_edge()) ) f = - curhe->opposite()->facet(); - else break; - } - //exit from the while if - //1. border - //2. not same type, then do not set visisted cause a BE - // follows a BH - } - - //next triangle adjacent to h2 (push_back) - if ( !(h2->is_border_edge()) ) - { - f = h2->opposite()->facet(); - curhe = h2; - while (cur_ridge_type == facet_ridge_type(f,curhe1,curhe2,r_type)) - { - //follow the ridge from curhe - if (f->is_visitedB()) break; - f->set_visitedB(true); - if (curhe->opposite() == curhe1) curhe = curhe2; - else curhe = curhe1; - cur_ridge_line->addback(curhe); - if ( !(curhe->is_border_edge()) ) f = - curhe->opposite()->facet(); - else break; - } - } - } -} - - -Ridge_type -facet_ridge_type(Facet_handle f, Halfedge_handle& he1, Halfedge_handle& - he2, Ridge_type r_type) -{ - //polyhedral data - //we have v1--h1-->v2--h2-->v3--h3-->v1 - Halfedge_handle h1 = f->halfedge(), h2, h3; - Vertex_handle v1, v2, v3; - v2 = h1->vertex(); - h2 = h1->next(); - v3 = h2->vertex(); - h3 = h2->next(); - v1 = h3->vertex(); - - //check for regular facet - if ( v1->d1()*v2->d1() * v1->d1()*v3->d1() * v2->d1()*v3->d1() < 0 ) return NONE; - - //determine potential crest color - Ridge_type crest_color = NONE; - if (r_type == CREST) - { - if ( std::fabs(v1->k1()+v2->k1()+v3->k1()) > std::fabs(v1->k2()+v2->k2()+v3->k2()) ) - crest_color = BC; - if ( std::fabs(v1->k1()+v2->k1()+v3->k1()) < std::fabs(v1->k2()+v2->k2()+v3->k2()) ) - crest_color = RC; - if ( std::fabs(v1->k1()+v2->k1()+v3->k1()) = std::fabs(v1->k2()+v2->k2()+v3->k2()) ) - return NONE; - } - - //compute Xing on the 3 edges - bool h1_is_crossed, h2_is_crossed, h3_is_crossed; - if ( r_type == BLUE || crest_color == BC ) - { - xing_on_edge(h1, h1_is_crossed, BLUE); - xing_on_edge(h2, h2_is_crossed, BLUE); - xing_on_edge(h3, h3_is_crossed, BLUE); - } - if ( r_type == RED || crest_color == RC ) - { - xing_on_edge(h1, h1_is_crossed, RED); - xing_on_edge(h2, h2_is_crossed, RED); - xing_on_edge(h3, h3_is_crossed, RED); - } - - //test of 2Xing - if ( !(h1_is_crossed && h2_is_crossed && !h3_is_crossed) && - !(h1_is_crossed && !h2_is_crossed && h3_is_crossed) && - (!h1_is_crossed && h2_is_crossed && h3_is_crossed) ) - return NONE; - if (h1_is_crossed && h2_is_crossed && !h3_is_crossed) - { - he1 = h1; - he2 = h2; - } - if (h1_is_crossed && !h2_is_crossed && h3_is_crossed) - { - he1 = h1; - he2 = h3; - } - if (!h1_is_crossed && h2_is_crossed && h3_is_crossed) - { - he1 = h2; - he2 = h3; - } - - Vertex_handle v_p1 = he1->opposite()->vertex(), v_q1 = he1->vertex(), - v_p2 = he2->opposite()->vertex(), v_q2 = he2->vertex(); // he1: p1->q1 - - if ( r_type == BLUE || crest_color == BC ) - { - FT coord1 = std::fabs(v_q1->b0()) / ( std::fabs(v_p1->b0()) + std::fabs(v_q1->b0()) ), - coord2 = std::fabs(v_q2->b0()) / ( std::fabs(v_p2->b0()) + std::fabs(v_q2->b0()) ); - Vector_3 r1 = (v_p1->point()-ORIGIN)*coord1 + (v_q1->point()-ORIGIN)*(1-coord1), - r2 = (v_p2->point()-ORIGIN)*coord2 + (v_q2->point()-ORIGIN)*(1-coord2); - int b_sign = b_sign_pointing_to_ridge(v1, v2, v3, r1, r2, BLUE); - if (b_sign == 1) { if (r_type == BLUE) return BE; else return BC;} - if (b_sign == -1) return BH; - } - if ( r_type == RED || crest_color == RC ) - { - FT coord1 = std::fabs(v_q1->b3()) / ( std::fabs(v_p1->b3()) + std::fabs(v_q1->b3()) ), - coord2 = std::fabs(v_q2->b3()) / ( std::fabs(v_p2->b3()) + std::fabs(v_q2->b3()) ); - Vector_3 r1 = (v_p1->point()-ORIGIN)*coord1 + (v_q1->point()-ORIGIN)*(1-coord1), - r2 = (v_p2->point()-ORIGIN)*coord2 + (v_q2->point()-ORIGIN)*(1-coord2); - int b_sign = b_sign_pointing_to_ridge(v1, v2, v3, r1, r2, RED); - if (b_sign == -1) { if (r_type == RED) return RE; else return RC;} - if (b_sign == 1) return BH; - } - assert(0); -} - -void -xing_on_edge(Halfedge_handle he, bool& is_crossed, Ridge_type color) -{ - is_crossed = false; - FT sign; - FT b_p, b_q; // extremalities at p and q for he: p->q - Vector_3 d_p = he->opposite()->vertex()->d1(), - d_q = he->vertex()->d1(); //ppal dir - if ( color == BLUE ) { - b_p = he->opposite()->vertex()->b0(); - b_q = he->vertex()->b0(); - } - else { - b_p = he->opposite()->vertex()->b3(); - b_q = he->vertex()->b3(); - } - if ( b_p == 0 && b_q == 0 ) return; - if ( b_p == 0 && b_q !=0 ) sign = d_p*d_q * b_q; - if ( b_p != 0 && b_q ==0 ) sign = d_p*d_q * b_p; - if ( b_p != 0 && b_q !=0 ) sign = d_p*d_q * b_p * b_q; - if ( sign < 0 ) is_crossed = true; -} - - - -//for a ridge segment in a triangle [r1,r2], let r = r2 - r1 and normalize, -//the projection of a point p on the line (r1,r2) is pp=r1+tr, with t=(p-r1)*r -//then the vector v starting at p is pointing to the ridge line (r1,r2) if -//(pp-p)*v >0 -int -b_sign_pointing_to_ridge(Vertex_handle v1, Vertex_handle v2, Vertex_handle v3, - Vector_3 r1, Vector_3 r2, Ridge_type color) -{ - Vector_3 r = r2 - r1, dv1, dv2, dv3; - FT bv1, bv2, bv3; - if ( color == BLUE ) { - bv1 = v1->b0(); - bv2 = v2->b0(); - bv3 = v3->b0(); - dv1 = v1->d1(); - dv2 = v2->d1(); - dv3 = v3->d1(); - } - else { - bv1 = v1->b3(); - bv2 = v2->b3(); - bv3 = v3->b3(); - dv1 = v1->d2(); - dv2 = v2->d2(); - dv3 = v3->d2(); - } - if ( r !=0 ) r = r/CGAL::sqrt(r*r); - FT sign1, sign2, sign3; - sign1 = (r1 - (v1->point()-ORIGIN) + (((v1->point()-ORIGIN)-r1)*r)*r )*dv1; - sign2 = (r1 - (v2->point()-ORIGIN) + (((v2->point()-ORIGIN)-r1)*r)*r )*dv2; - sign3 = (r1 - (v3->point()-ORIGIN) + (((v3->point()-ORIGIN)-r1)*r)*r )*dv3; - - int compt = 0; - if ( sign1 > 0 ) compt++; else if (sign1 < 0) compt--; - if ( sign2 > 0 ) compt++; else if (sign2 < 0) compt--; - if ( sign3 > 0 ) compt++; else if (sign3 < 0) compt--; - - if (compt > 0) return 1; else return -1; -} - - - - - - - - - - - - -CGAL_END_NAMESPACE - -#endif