diff --git a/Ridges_3/include/CGAL/PolyhedralSurf_neighbors.h b/Ridges_3/include/CGAL/PolyhedralSurf_neighbors.h index f044755ab84..2d98405f278 100644 --- a/Ridges_3/include/CGAL/PolyhedralSurf_neighbors.h +++ b/Ridges_3/include/CGAL/PolyhedralSurf_neighbors.h @@ -6,9 +6,11 @@ CGAL_BEGIN_NAMESPACE -//Element of the priority queue. A gate is a halfedge and a number -//giving the max distance from v to the vertices of the triangle -//incident to the halfedge. +//--------------------------------------------------------------------------- +//T_Gate : element of the priority queue. A gate is a halfedge and a +//number giving the max distance from v to the vertices of the +//triangle incident to the halfedge. +//--------------------------------------------------------------------------- template < class TPoly > class T_Gate { public: @@ -44,8 +46,10 @@ public: Halfedge_handle he() {return m_he;} }; +//--------------------------------------------------------------------------- // functor for priority queue // order so than the top element is the smallest in the queue +//--------------------------------------------------------------------------- template struct compare_gates { @@ -56,8 +60,11 @@ struct compare_gates } }; -//MAIN class for computation, it uses the class Gate and the functor -//compare_gates for the definition of a priority queue +//--------------------------------------------------------------------------- +//T_PolyhedralSurf_neighbors : MAIN class for computation, it uses the +//class Gate and the functor compare_gates for the definition of a +//priority queue +//--------------------------------------------------------------------------- template < class TPoly > class T_PolyhedralSurf_neighbors { public: @@ -68,31 +75,53 @@ public: typedef typename TPoly::Halfedge_handle Halfedge_handle; typedef typename TPoly::Halfedge_around_vertex_circulator Halfedge_around_vertex_circulator; + typedef typename TPoly::Vertex_iterator Vertex_iterator; typedef T_Gate Gate; public: - // vertex_neigh stores the vertex v and its 1Ring neighbors - // contour stores halfedges, oriented CW, following the 1Ring disk border - // OneRingSize is the max distance from v to its OneRing neighbors - static void compute_one_ring(Vertex_handle v, + T_PolyhedralSurf_neighbors(TPoly& P); + // vertex_neigh stores the vertex v and its 1Ring neighbors contour + // stores halfedges, oriented CW, following the 1Ring disk border + // OneRingSize is the max distance from v to its OneRing + // neighbors. (the tag is_visited is not mofified) + void compute_one_ring(Vertex_handle v, std::vector &vertex_neigh, std::list &contour, FT &OneRingSize); // call compute_one_ring and expand the contour (circle of halfedges // CW), vertex_neigh are vertices on and inside the contour (there - // ring index is set to 1, but reset at the end), size is such that - // gates with distance less than size*OneRingSize are processed - static void compute_neighbors(Vertex_handle v, + // tag is_visited is set to true, but reset to false at the end), + // size is such that gates with distance less than size*OneRingSize + // are processed + void compute_neighbors(Vertex_handle v, std::vector &vertex_neigh, std::list &contour, FT size); - //vertex indices are initialised to -1 - static void reset_ring_indices(std::vector &vces); + //vertex tags is_visited are set to false + void reset_is_visited_map(std::vector &vces); + protected: + //tag to visit vertices + struct Vertex_cmp{//comparison is wrt vertex addresses + bool operator()(Vertex_handle a, Vertex_handle b) const{ + return &*a < &*b; + } + }; + typedef std::map Vertex2bool_map_type; + Vertex2bool_map_type is_visited_map; }; //////////////IMPLEMENTATION////////////////////////// ////////////////////////////////////////////////////// +template < class TPoly > +T_PolyhedralSurf_neighbors < TPoly >:: +T_PolyhedralSurf_neighbors(TPoly& P) +{ + //init the is_visited_map + Vertex_iterator itb = P.vertices_begin(), ite = P.vertices_end(); + for(;itb!=ite;itb++) is_visited_map[itb] = false; +} + template < class TPoly > void T_PolyhedralSurf_neighbors < TPoly >:: compute_one_ring(Vertex_handle v, @@ -151,13 +180,11 @@ compute_neighbors(Vertex_handle v, FT OneRingSize; compute_one_ring(v, vertex_neigh, contour, OneRingSize); FT d_max = OneRingSize*size; - std::priority_queue< Gate, - std::vector< Gate >, - compare_gates< Gate > > GatePQ; + std::priority_queue< Gate, std::vector< Gate >, compare_gates< Gate > > GatePQ; // tag neighbors typename std::vector::iterator itbv = vertex_neigh.begin(), itev = vertex_neigh.end(); - for (; itbv != itev; itbv++) (*itbv)->setRingIndex(1); + for (; itbv != itev; itbv++) is_visited_map.find(*itbv)->second = true; // init GatePQ typename std::list::iterator itb = contour.begin(), @@ -243,10 +270,10 @@ compute_neighbors(Vertex_handle v, continue; } v1 = he->next()->vertex(); - if ( v1->getRingIndex() == -1 ) + if ( !is_visited_map.find(v1)->second ) { // case 1 //vertex - v1->setRingIndex(1); + is_visited_map.find(v1)->second = true; vertex_neigh.push_back(v1); //contour he1 = he->prev()->opposite(); @@ -275,16 +302,16 @@ compute_neighbors(Vertex_handle v, /* } */ /* //debug */ - reset_ring_indices(vertex_neigh); + reset_is_visited_map(vertex_neigh); } template < class TPoly > void T_PolyhedralSurf_neighbors < TPoly >:: -reset_ring_indices(std::vector &vces) +reset_is_visited_map(std::vector &vces) { typename std::vector::iterator itb = vces.begin(), ite = vces.end(); - for (;itb != ite; itb++) (*itb)->resetRingIndex(); + for (;itb != ite; itb++) is_visited_map[*itb] = false; } CGAL_END_NAMESPACE diff --git a/Ridges_3/include/CGAL/Ridges.h b/Ridges_3/include/CGAL/Ridges.h index b1f7f225e26..10403858f92 100644 --- a/Ridges_3/include/CGAL/Ridges.h +++ b/Ridges_3/include/CGAL/Ridges.h @@ -1,138 +1,78 @@ #ifndef _RIDGE_3_H_ #define _RIDGE_3_H_ -#include #include #include -#include "PolyhedralSurf_neighbors.h" -#include "Umbilic.h" +#include +#include + +#include +#include +#include //note : one has to orient monge normals according to mesh normals to //define min/max curv CGAL_BEGIN_NAMESPACE - -enum Ridge_type {NONE=0, BLUE_RIDGE, RED_RIDGE, CREST, BE, BH, BC, RE, RH, RC}; + +enum Ridge_type {NONE=0, BLUE_RIDGE, RED_RIDGE, CREST, + BLUE_ELLIPTIC_RIDGE, BLUE_HYPERBOLIC_RIDGE, BLUE_CREST, + RED_ELLIPTIC_RIDGE, RED_HYPERBOLIC_RIDGE, RED_CREST}; //--------------------------------------------------------------------------- -//Ridge_line : a connected sequence of edges crossed by a ridge, with -//type and weigths +//Ridge_line : a connected sequence of edges of a Poly crossed by a +//ridge (with a barycentric coordinate to compute the crossing point), +//with type and weigths. //-------------------------------------------------------------------------- template < class Poly > class Ridge_line { public: - typedef typename Poly::Traits::FT FT; - typedef typename Poly::Traits::Vector_3 Vector_3; - typedef typename Poly::Traits::Point_3 Point_3; - typedef typename Poly::Vertex_handle Vertex_handle; - typedef typename Poly::Halfedge_handle Halfedge_handle; - typedef std::pair< Halfedge_handle, FT> ridge_he; + typedef typename Poly::Traits::FT FT; + typedef typename Poly::Traits::Vector_3 Vector_3; + typedef typename Poly::Traits::Point_3 Point_3; + typedef typename Poly::Halfedge_handle Halfedge_handle; + typedef std::pair< Halfedge_handle, FT> ridge_he; protected: - Ridge_type m_line_type;//one of BE, BH, BC, RE, RH or RC + //one of BLUE_ELLIPTIC_RIDGE, BLUE_HYPERBOLIC_RIDGE, BLUE_CREST, + //RED_ELLIPTIC_RIDGE, RED_HYPERBOLIC_RIDGE or RED_CREST + Ridge_type m_line_type; std::list m_line; FT m_strength; FT m_sharpness; public: const Ridge_type line_type() const {return m_line_type;} + Ridge_type& line_type() {return m_line_type;} + const FT strength() const {return m_strength;} + FT& strength() {return m_strength;} + const FT sharpness() const {return m_sharpness;} - std::list* line() { return &m_line;} + FT& sharpness() {return m_sharpness;} + const std::list* line() const { return &m_line;} + std::list* line() { return &m_line;} //constructor - //a ridge line begins with a segment in a triangle - Ridge_line( Halfedge_handle h1, Halfedge_handle h2, Ridge_type r_type); - - //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); - - //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); + Ridge_line(); - void dump_4ogl(std::ostream& out_stream); + /* The output is : + line_type, strength, sharpness, list of points of the polyline. + */ + void dump_4ogl(std::ostream& out_stream) ; + void dump_verbose(std::ostream& out_stream) const ; }; +//-------------------------------------------------------------------------- // IMPLEMENTATION OF Ridge_line members -////////////////////////////////////////////////////////////////////////////// +//-------------------------------------------------------------------------- //constructor - template < class Poly > - Ridge_line:: - Ridge_line( Halfedge_handle h1, Halfedge_handle h2, Ridge_type r_type) - : m_line_type(r_type), m_strength(0.), m_sharpness(0.) - { - m_line.push_back(ridge_he(h1, bary_coord(h1))); - addback(h2); - } - template < class Poly > -void Ridge_line:: -addback( Halfedge_handle he) -{ - Halfedge_handle he_cur = ( --(m_line.end()) )->first; - FT coord_cur = ( --(m_line.end()) )->second;//bary_coord(he_cur); - FT coord = bary_coord(he); - 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;//abs value of the ppal curvature at the Xing point on he. - if ( (m_line_type == BE) || (m_line_type == BH) || (m_line_type == BC) ) { - k = CGAL::abs(v_p->k1()) * coord + CGAL::abs(v_q->k1()) * (1-coord) ; - } - if ( (m_line_type == RE) || (m_line_type == RH) || (m_line_type == RC) ) { - k = CGAL::abs(v_p->k2()) * coord + CGAL::abs(v_q->k2()) * (1-coord) ; - } - 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)); - m_strength += k * CGAL::sqrt(segment * segment); - //TODO update sharpness - m_line.push_back( ridge_he(he, coord)); -} - -template < class Poly > -void Ridge_line:: -addfront( Halfedge_handle he) -{ - Halfedge_handle he_cur = ( m_line.begin() )->first; - FT coord_cur = ( m_line.begin() )->second; - FT coord = bary_coord(he); - 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 ( (m_line_type == BE) || (m_line_type == BH) || (m_line_type == BC) ) { - k = CGAL::abs(v_p->k1()) * coord + CGAL::abs(v_q->k1()) * (1-coord) ; - } - if ( (m_line_type == RE) || (m_line_type == RH) || (m_line_type == RC) ) { - k = CGAL::abs(v_p->k2()) * coord + CGAL::abs(v_q->k2()) * (1-coord) ; - } - 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)); - m_strength += k * CGAL::sqrt(segment * segment); - //TODO update sharpness - m_line.push_front( ridge_he(he, coord)); -} - - -template < class Poly > -typename Poly::Traits::FT Ridge_line:: -bary_coord( Halfedge_handle he) -{ - FT b_p, b_q; // extremalities at p and q for he: p->q - if ( (m_line_type == BE) || (m_line_type == BH) || (m_line_type == BC) ) { - b_p = he->opposite()->vertex()->b0(); - b_q = he->vertex()->b0(); - } - if ( (m_line_type == RE) || (m_line_type == RH) || (m_line_type == RC) ) { - b_p = he->opposite()->vertex()->b3(); - b_q = he->vertex()->b3(); - } - return CGAL::abs(b_q) / ( CGAL::abs(b_q) + CGAL::abs(b_p) ); -} +Ridge_line:: +Ridge_line() : m_strength(0.), m_sharpness(0.) {} + template < class Poly > void Ridge_line:: @@ -156,29 +96,80 @@ dump_4ogl(std::ostream& out_stream) out_stream << std::endl; } +//verbose output +template < class Poly > +void Ridge_line:: +dump_verbose(std::ostream& out_stream) const +{ + out_stream << "Line type is : " << line_type() << std::endl + << "Strength is : " << strength() << std::endl + << "Sharpness is : " << sharpness() << std::endl + << "Polyline point coordinates are : " << std::endl; + + typename std::list::const_iterator + iter = line()->begin(), + ite = line()->end(); + for (;iter!=ite;iter++){ + //he: p->q, r is the crossing point + Point_3 p = iter->first->opposite()->vertex()->point(), + q = iter->first->vertex()->point(); + Vector_3 r = (p-CGAL::ORIGIN)*iter->second + + (q-CGAL::ORIGIN)*(1-iter->second); + out_stream << r << std::endl; + } +} + +template +std::ostream& +operator<<(std::ostream& out_stream, const Ridge_line& ridge_line) +{ + ridge_line.dump_verbose(out_stream); + return out_stream; +} + +//--------------------------------------------------------------------------- +//Differential_quantities +//-------------------------------------------------------------------------- +template < class Vertex2FTPropertyMap, class Vertex2VectorPropertyMap> + class Differential_quantities { + public : + Vertex2FTPropertyMap vertex2k1_pm, vertex2k2_pm, + vertex2b0_pm, vertex2b3_pm, + vertex2P1_pm, vertex2P2_pm; + Vertex2VectorPropertyMap vertex2d1_pm, vertex2d2_pm; +}; + + + //--------------------------------------------------------------------------- //Ridge_approximation //-------------------------------------------------------------------------- -template < class Poly, class OutputIt > -class Ridge_approximation +template < class Poly, class OutputIt, class Vertex2FTPropertyMap, class Vertex2VectorPropertyMap > + class Ridge_approximation { -public: - typedef typename Poly::Traits::FT FT; - typedef typename Poly::Traits::Vector_3 Vector_3; - typedef typename Poly::Vertex_handle Vertex_handle; - typedef typename Poly::Halfedge_handle Halfedge_handle; - typedef typename Poly::Facet_handle Facet_handle; - typedef typename Poly::Facet_iterator Facet_iterator; - typedef Ridge_line Ridge_line; - // typedef T_PolyhedralSurf_neighbors Poly_neighbors;//for umbilics + public: + typedef typename Poly::Traits::FT FT; + typedef typename Poly::Traits::Vector_3 Vector_3; + typedef typename Poly::Vertex_handle Vertex_handle; + typedef typename Poly::Halfedge_handle Halfedge_handle; + typedef typename Poly::Facet_handle Facet_handle; + typedef typename Poly::Facet_iterator Facet_iterator; + + typedef std::pair< Halfedge_handle, FT> ridge_he; + typedef Ridge_line Ridge_line; + //are ridges tagged as elliptic or hyperbolic using 3rd or 4th order //differential quantitities? + //tag_3 is not working or badly, to be removed or improved? enum Tag_order {Tag_3 = 3, Tag_4 = 4}; - -public: - Ridge_approximation(){}; - + + public: + Ridge_approximation(Poly &P, + Vertex2FTPropertyMap vertex2k1_pm, Vertex2FTPropertyMap vertex2k2_pm, + Vertex2FTPropertyMap vertex2b0_pm, Vertex2FTPropertyMap vertex2b3_pm, + Vertex2FTPropertyMap vertex2P1_pm, Vertex2FTPropertyMap vertex2P2_pm, + Vertex2VectorPropertyMap vertex2d1_pm, Vertex2VectorPropertyMap vertex2d2_pm); OutputIt compute_all_ridges(Poly &P, OutputIt it, Tag_order ord = Tag_3); //Find BLUE_RIDGE, RED_RIDGE or CREST ridges @@ -194,9 +185,26 @@ public: // void compute_umbilics(Poly &P );//container, class for umbilics? protected: + Poly* P; + FT model_size;//radius of the smallest enclosing sphere of the Poly + //used to make the sharpness scale independant and iso indep + //tag to visit faces + struct Facet_cmp{ //comparison is wrt facet addresses + bool operator()(Facet_handle a, Facet_handle b) const{ + return &*a < &*b; + } + }; + typedef std::map Facet2bool_map_type; + Facet2bool_map_type is_visited_map; + + //Property maps + Vertex2FTPropertyMap k1, k2, b0, b3, P1, P2; + Vertex2VectorPropertyMap d1, d2; + //is a facet crossed by a BLUE, RED or CREST ridge? if so, return - //the crossed edges and more precise type from BE, BH, BC, RE, RH, - //RC or NONE + //the crossed edges and more precise type from BLUE_ELLIPTIC_RIDGE, + //BLUE_HYPERBOLIC_RIDGE, BLUE_CREST, RED_ELLIPTIC_RIDGE, + //RED_HYPERBOLIC_RIDGE, RED_CREST or NONE Ridge_type facet_ridge_type(Facet_handle f, Halfedge_handle& he1, Halfedge_handle& he2, @@ -223,20 +231,62 @@ protected: //pointing to the ridge line (r1,r2) if (pp-p)*v >0. Return the sign //of b, for a ppal direction pointing to the ridge segment, //appearing at least at two vertices of the facet. - // for color = BLUE_RIDGE, sign = 1 if BE, -1 if BH - // for color = RED_RIDGE, sign = -1 if RE, 1 if RH + // + // for color = BLUE_RIDGE, sign = 1 if BLUE_ELLIPTIC_RIDGE, -1 if + // BLUE_HYPERBOLIC_RIDGE + // + // for color = RED_RIDGE, sign = -1 if RED_ELLIPTIC_RIDGE, 1 if + // RED_HYPERBOLIC_RIDGE int b_sign_pointing_to_ridge(Vertex_handle v1, Vertex_handle v2, Vertex_handle v3, Vector_3 r1, Vector_3 r2, Ridge_type color); + + //a ridge line begins with a segment in a triangle + void init_ridge_line(Ridge_line* ridge_line, + Halfedge_handle h1, Halfedge_handle h2, + Ridge_type r_type); + + //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(Ridge_line* ridge_line, Halfedge_handle he, Ridge_type r_type); + void addfront(Ridge_line* ridge_line, Halfedge_handle he, Ridge_type r_type); + + //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, Ridge_type r_type); }; // IMPLEMENTATION OF Ridge_approximation members ///////////////////////////////////////////////////////////////////////////// -template < class Poly, class OutputIt > -OutputIt Ridge_approximation:: +//contructor +template < class Poly, class OutputIt, class Vertex2FTPropertyMap, class Vertex2VectorPropertyMap > + Ridge_approximation:: +Ridge_approximation(Poly &P, + Vertex2FTPropertyMap vertex2k1_pm, Vertex2FTPropertyMap vertex2k2_pm, + Vertex2FTPropertyMap vertex2b0_pm, Vertex2FTPropertyMap vertex2b3_pm, + Vertex2FTPropertyMap vertex2P1_pm, Vertex2FTPropertyMap vertex2P2_pm, + Vertex2VectorPropertyMap vertex2d1_pm, Vertex2VectorPropertyMap vertex2d2_pm) +: P(&P), k1(vertex2k1_pm), k2(vertex2k2_pm), b0(vertex2b0_pm), b3(vertex2b3_pm), + P1(vertex2P1_pm), P2(vertex2P2_pm), d1(vertex2d1_pm), d2(vertex2d2_pm) +{ + CGAL::Min_sphere_d > + min_sphere(P.points_begin(), P.points_end()); + model_size = min_sphere.squared_radius(); + //maybe better to use CGAL::Min_sphere_of_spheres_d ?? but need to create spheres? + + //init the is_visited_map + Facet_iterator itb = P.facets_begin(), ite = P.facets_end(); + for(;itb!=ite;itb++) is_visited_map[itb] = false; +} + + + +template < class Poly, class OutputIt, class Vertex2FTPropertyMap, class Vertex2VectorPropertyMap > +OutputIt Ridge_approximation:: compute_all_ridges(Poly &P, OutputIt it, Tag_order ord) { compute_ridges(P, BLUE_RIDGE, it, ord); @@ -245,37 +295,35 @@ compute_all_ridges(Poly &P, OutputIt it, Tag_order ord) return it; } -template < class Poly, class OutputIt > -void Ridge_approximation:: +template < class Poly, class OutputIt, class Vertex2FTPropertyMap, class Vertex2VectorPropertyMap > + void Ridge_approximation:: compute_ridges(Poly &P, Ridge_type r_type, OutputIt ridge_lines_it, Tag_order ord) { //set all facets non visited Facet_iterator itb = P.facets_begin(), ite = P.facets_end(); - for(;itb!=ite;itb++) itb->reset_is_visited(); + for(;itb!=ite;itb++) is_visited_map[itb] = false; itb = P.facets_begin(); for(;itb!=ite;itb++) { - Facet_handle f= &(*itb); - if (f->is_visited()) continue; - f->set_visited(true); + Facet_handle f = itb; + if (is_visited_map.find(f)->second) continue; + is_visited_map.find(f)->second = true; Halfedge_handle h1, h2, curhe1, curhe2, curhe; //h1 h2 are the hedges crossed if any, r_type should be - //BLUE_RIDGE, RED_RIDGE or CREST ; cur_ridge_type should be BE, - //BH, BC, RE, RH, RC or NONE + //BLUE_RIDGE, RED_RIDGE or CREST ; cur_ridge_type should be + //BLUE_ELLIPTIC_RIDGE, BLUE_HYPERBOLIC_RIDGE, BLUE_CREST, + //RED_ELLIPTIC_RIDGE, RED_HYPERBOLIC_RIDGE, RED_CREST or NONE Ridge_type cur_ridge_type = facet_ridge_type(f,h1,h2,r_type, ord); if ( cur_ridge_type == NONE ) continue; //a ridge_line is begining and stored - // Ridge_line cur_ridge_line(h1,h2,cur_ridge_type); - Ridge_line* cur_ridge_line = new Ridge_line(h1,h2,cur_ridge_type); + Ridge_line* cur_ridge_line = new Ridge_line(); + init_ridge_line(cur_ridge_line, h1, h2, cur_ridge_type); *ridge_lines_it++ = cur_ridge_line; - //debug - // cur_ridge_line->dump_4ogl(std::cout); - // std::cout << "??????????????????????????" << endl; - + //next triangle adjacent to h1 (push_front) if ( !(h1->is_border_edge()) ) { @@ -285,19 +333,19 @@ compute_ridges(Poly &P, Ridge_type r_type, r_type, ord)) { //follow the ridge from curhe - if (f->is_visited()) break; - f->set_visited(true); + if (is_visited_map.find(f)->second) break; + is_visited_map.find(f)->second = true; if (curhe->opposite() == curhe1) curhe = curhe2; else curhe = curhe1;//curhe stays at the ridge extremity - cur_ridge_line->addfront(curhe); + addfront(cur_ridge_line, curhe, cur_ridge_type); if ( !(curhe->is_border_edge()) ) f = curhe->opposite()->facet(); else break; } //exit from the while if //1. border or already visited (this is a ridge loop) - //2. not same type, then do not set visisted cause a BE - // follows a BH + //2. not same type, then do not set visisted cause a BLUE_ELLIPTIC_RIDGE + // follows a BLUE_HYPERBOLIC_RIDGE } //next triangle adjacent to h2 (push_back) @@ -309,11 +357,11 @@ compute_ridges(Poly &P, Ridge_type r_type, facet_ridge_type(f,curhe1,curhe2,r_type, ord)) { //follow the ridge from curhe - if (f->is_visited()) break; - f->set_visited(true); + if (is_visited_map.find(f)->second) break; + is_visited_map.find(f)->second = true; if (curhe->opposite() == curhe1) curhe = curhe2; else curhe = curhe1; - cur_ridge_line->addback(curhe); + addback(cur_ridge_line, curhe, cur_ridge_type); if ( !(curhe->is_border_edge()) ) f = curhe->opposite()->facet(); else break; @@ -322,8 +370,8 @@ compute_ridges(Poly &P, Ridge_type r_type, } } -template < class Poly, class OutputIt > -Ridge_type Ridge_approximation:: +template < class Poly, class OutputIt, class Vertex2FTPropertyMap, class Vertex2VectorPropertyMap > +Ridge_type Ridge_approximation:: facet_ridge_type(Facet_handle f, Halfedge_handle& he1, Halfedge_handle& he2, Ridge_type r_type, Tag_order ord) { @@ -339,32 +387,32 @@ facet_ridge_type(Facet_handle f, Halfedge_handle& he1, Halfedge_handle& //check for regular facet //i.e. if there is a coherent orientation of ppal dir at the facet vertices - if ( v1->d1()*v2->d1() * v1->d1()*v3->d1() * v2->d1()*v3->d1() < 0 ) + if ( d1[v1]*d1[v2] * d1[v1]*d1[v3] * d1[v2]*d1[v3] < 0 ) return NONE; //determine potential crest color - //BC if |sum(k1)|>|sum(k2)| sum over facet vertices vi - //RC if |sum(k1)|<|sum(k2)| + //BLUE_CREST if |sum(k1)|>|sum(k2)| sum over facet vertices vi + //RED_CREST if |sum(k1)|<|sum(k2)| Ridge_type crest_color = NONE; if (r_type == CREST) { - if ( CGAL::abs(v1->k1()+v2->k1()+v3->k1()) > CGAL::abs(v1->k2()+v2->k2()+v3->k2()) ) - crest_color = BC; - if ( CGAL::abs(v1->k1()+v2->k1()+v3->k1()) < CGAL::abs(v1->k2()+v2->k2()+v3->k2()) ) - crest_color = RC; - if ( CGAL::abs(v1->k1()+v2->k1()+v3->k1()) == CGAL::abs(v1->k2()+v2->k2()+v3->k2()) ) + if ( CGAL::abs(k1[v1]+k1[v2]+k1[v3]) > CGAL::abs(k2[v1]+k2[v2]+k2[v3]) ) + crest_color = BLUE_CREST; + if ( CGAL::abs(k1[v1]+k1[v2]+k1[v3]) < CGAL::abs(k2[v1]+k2[v2]+k2[v3]) ) + crest_color = RED_CREST; + if ( CGAL::abs(k1[v1]+k1[v2]+k1[v3]) == CGAL::abs(k2[v1]+k2[v2]+k2[v3]) ) return NONE; } //compute Xing on the 3 edges bool h1_is_crossed = false, h2_is_crossed = false, h3_is_crossed = false; - if ( r_type == BLUE_RIDGE || crest_color == BC ) + if ( r_type == BLUE_RIDGE || crest_color == BLUE_CREST ) { xing_on_edge(h1, h1_is_crossed, BLUE_RIDGE); xing_on_edge(h2, h2_is_crossed, BLUE_RIDGE); xing_on_edge(h3, h3_is_crossed, BLUE_RIDGE); } - if ( r_type == RED_RIDGE || crest_color == RC ) + if ( r_type == RED_RIDGE || crest_color == RED_CREST ) { xing_on_edge(h1, h1_is_crossed, RED_RIDGE); xing_on_edge(h2, h2_is_crossed, RED_RIDGE); @@ -398,70 +446,72 @@ facet_ridge_type(Facet_handle f, Halfedge_handle& he1, Halfedge_handle& 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_RIDGE || crest_color == BC ) { - FT coord1 = CGAL::abs(v_q1->b0()) / ( CGAL::abs(v_p1->b0()) + - CGAL::abs(v_q1->b0()) ), - coord2 = CGAL::abs(v_q2->b0()) / ( CGAL::abs(v_p2->b0()) + - CGAL::abs(v_q2->b0()) ); + if ( r_type == BLUE_RIDGE || crest_color == BLUE_CREST ) { + FT coord1 = CGAL::abs(b0[v_q1]) / ( CGAL::abs(b0[v_p1]) + + CGAL::abs(b0[v_q1]) ), + coord2 = CGAL::abs(b0[v_q2]) / ( CGAL::abs(b0[v_p2]) + + CGAL::abs(b0[v_q2]) ); if ( ord == Tag_3 ) { 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_RIDGE); - if (r_type == CREST) {if (b_sign == 1) return BC; else return NONE;} - if (b_sign == 1) return BE; else return BH; + if (r_type == CREST) {if (b_sign == 1) return BLUE_CREST; + else return NONE;} + if (b_sign == 1) return BLUE_ELLIPTIC_RIDGE; + else return BLUE_HYPERBOLIC_RIDGE; } else {//ord == Tag_4, check the sign of the meanvalue of the signs // of P1 at the two crossing points - FT sign_P1 = v_p1->P1()*coord1 + v_q1->P1()*(1-coord1) - + v_p2->P1()*coord2 + v_q2->P1()*(1-coord2); - if (r_type == CREST) {if ( sign_P1 < 0 ) return BC; else return NONE;} - if ( sign_P1 < 0 ) return BE; else return BH; + FT sign_P1 = P1[v_p1]*coord1 + P1[v_q1]*(1-coord1) + + P1[v_p2]*coord2 + P1[v_q2]*(1-coord2); + if (r_type == CREST) {if ( sign_P1 < 0 ) return BLUE_CREST; else return NONE;} + if ( sign_P1 < 0 ) return BLUE_ELLIPTIC_RIDGE; else return BLUE_HYPERBOLIC_RIDGE; } } - if ( r_type == RED_RIDGE || crest_color == RC ) { - FT coord1 = CGAL::abs(v_q1->b3()) / ( CGAL::abs(v_p1->b3()) + - CGAL::abs(v_q1->b3()) ), - coord2 = CGAL::abs(v_q2->b3()) / ( CGAL::abs(v_p2->b3()) + - CGAL::abs(v_q2->b3()) ); + if ( r_type == RED_RIDGE || crest_color == RED_CREST ) { + FT coord1 = CGAL::abs(b3[v_q1]) / ( CGAL::abs(b3[v_p1]) + + CGAL::abs(b3[v_q1]) ), + coord2 = CGAL::abs(b3[v_q2]) / ( CGAL::abs(b3[v_p2]) + + CGAL::abs(b3[v_q2]) ); if ( ord == Tag_3 ) { 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_RIDGE); - if (r_type == CREST) {if (b_sign == -1) return RC; else return NONE;} - if (b_sign == -1) return RE; else return RH; + if (r_type == CREST) {if (b_sign == -1) return RED_CREST; else return NONE;} + if (b_sign == -1) return RED_ELLIPTIC_RIDGE; else return RED_HYPERBOLIC_RIDGE; } else {//ord == Tag_4, check the sign of the meanvalue of the signs // of P2 at the two crossing points - FT sign_P2 = v_p1->P2()*coord1 + v_q1->P2()*(1-coord1) - + v_p2->P2()*coord2 + v_q2->P2()*(1-coord2); - if (r_type == CREST) {if ( sign_P2 < 0 ) return RC; else return NONE;} - if ( sign_P2 < 0 ) return RE; else return RH; + FT sign_P2 = P2[v_p1]*coord1 + P2[v_q1]*(1-coord1) + + P2[v_p2]*coord2 + P2[v_q2]*(1-coord2); + if (r_type == CREST) {if ( sign_P2 < 0 ) return RED_CREST; else return NONE;} + if ( sign_P2 < 0 ) return RED_ELLIPTIC_RIDGE; else return RED_HYPERBOLIC_RIDGE; } } assert(0);//should return before! } -template < class Poly, class OutputIt > -void Ridge_approximation:: +template < class Poly, class OutputIt, class Vertex2FTPropertyMap, class Vertex2VectorPropertyMap > +void Ridge_approximation:: 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 + Vector_3 d_p = d1[he->opposite()->vertex()], + d_q = d1[he->vertex()]; //ppal dir if ( color == BLUE_RIDGE ) { - b_p = he->opposite()->vertex()->b0(); - b_q = he->vertex()->b0(); + b_p = b0[he->opposite()->vertex()]; + b_q = b0[he->vertex()]; } else { - b_p = he->opposite()->vertex()->b3(); - b_q = he->vertex()->b3(); + b_p = b3[he->opposite()->vertex()]; + b_q = b3[he->vertex()]; } if ( b_p == 0 && b_q == 0 ) return; if ( b_p == 0 && b_q !=0 ) sign = d_p*d_q * b_q; @@ -471,28 +521,28 @@ xing_on_edge(Halfedge_handle he, bool& is_crossed, Ridge_type color) } -template < class Poly, class OutputIt > -int Ridge_approximation:: -b_sign_pointing_to_ridge(Vertex_handle v1, Vertex_handle v2, Vertex_handle v3, - Vector_3 r1, Vector_3 r2, Ridge_type color) +template < class Poly, class OutputIt, class Vertex2FTPropertyMap, class Vertex2VectorPropertyMap > + int Ridge_approximation:: + 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_RIDGE ) { - bv1 = v1->b0(); - bv2 = v2->b0(); - bv3 = v3->b0(); - dv1 = v1->d1(); - dv2 = v2->d1(); - dv3 = v3->d1(); + bv1 = b0[v1]; + bv2 = b0[v2]; + bv3 = b0[v3]; + dv1 = d1[v1]; + dv2 = d1[v2]; + dv3 = d1[v3]; } else { - bv1 = v1->b3(); - bv2 = v2->b3(); - bv3 = v3->b3(); - dv1 = v1->d2(); - dv2 = v2->d2(); - dv3 = v3->d2(); + bv1 = b3[v1]; + bv2 = b3[v2]; + bv3 = b3[v3]; + dv1 = d2[v1]; + dv2 = d2[v2]; + dv3 = d2[v3]; } if ( r != CGAL::NULL_VECTOR ) r = r/CGAL::sqrt(r*r); FT sign1, sign2, sign3; @@ -508,83 +558,112 @@ b_sign_pointing_to_ridge(Vertex_handle v1, Vertex_handle v2, Vertex_handle v3, if (compt > 0) return 1; else return -1; } -//--------------------------------------------------------------------------- -//Umbilic_approximation -//-------------------------------------------------------------------------- -template < class Poly, class OutputIt > -class Umbilic_approximation +template < class Poly, class OutputIt, class Vertex2FTPropertyMap, class Vertex2VectorPropertyMap > + void Ridge_approximation:: +init_ridge_line(Ridge_line* ridge_line, + Halfedge_handle h1, Halfedge_handle h2, + Ridge_type r_type) { -public: - typedef typename Poly::Traits::FT FT; - typedef typename Poly::Traits::Vector_3 Vector_3; - typedef typename Poly::Vertex_handle Vertex_handle; - typedef typename Poly::Halfedge_handle Halfedge_handle; - typedef typename Poly::Facet_handle Facet_handle; - typedef typename Poly::Facet_iterator Facet_iterator; - typedef typename Poly::Vertex_iterator Vertex_iterator; - typedef T_PolyhedralSurf_neighbors Poly_neighbors; - typedef Umbilic Umbilic; - CGAL::Abs cgal_abs; - static FT neigh_size;//the size of neighbourhood for umbilic - // computation is (neigh_size * OneRingSize) - - Umbilic_approximation(){}; - OutputIt compute(Poly &P, OutputIt it, FT size); - -}; - -template < class Poly, class OutputIt > - OutputIt Umbilic_approximation< Poly, OutputIt >:: -compute(Poly &P, OutputIt umbilics_it, FT size) -{ - std::vector vces; - std::list contour; - double umbilicEstimatorVertex, umbilicEstimatorNeigh; - - bool is_umbilic = true; - - //MAIN loop on P vertices - Vertex_iterator itb = P.vertices_begin(), ite = P.vertices_end(); - for (;itb != ite; itb++) { - Vertex_handle vh = itb; - umbilicEstimatorVertex = cgal_abs(vh->k1()-vh->k2()); - //reset vector, list and bool - vces.clear(); - contour.clear(); - is_umbilic = true; - Poly_neighbors::compute_neighbors(vh, vces, contour, size); - - - // OPTIONAL: avoid umbilics whose contours touch the border - typename std::list::iterator itb_cont = contour.begin(), - ite_cont = contour.end(); - for (; itb_cont != ite_cont; itb_cont++) - if ( (*itb_cont)->is_border() ) {is_umbilic = false; continue;} - if (is_umbilic == false) continue; - - //is v an umbilic? - //a priori is_umbilic = true, and it switches to false as soon as a - // neigh vertex has a lower umbilicEstimator value - typename std::vector::iterator itbv = vces.begin(), - itev = vces.end(); - assert(*itbv == vh); - itbv++; - for (; itbv != itev; itbv++) - { umbilicEstimatorNeigh = cgal_abs( (*itbv)->k1() - (*itbv)->k2() ); - if ( umbilicEstimatorNeigh < umbilicEstimatorVertex ) - {is_umbilic = false; break;} - } - if (is_umbilic == false) continue; - - //v is an umbilic, compute the index - Umbilic* cur_umbilic = new Umbilic(vh, contour); - cur_umbilic->compute_type(); - *umbilics_it++ = cur_umbilic; - } - return umbilics_it; + ridge_line->line_type() = r_type; + ridge_line->line()->push_back(ridge_he(h1, bary_coord(h1,r_type))); + addback(ridge_line, h2, r_type); } +template < class Poly, class OutputIt, class Vertex2FTPropertyMap, class Vertex2VectorPropertyMap > + void Ridge_approximation:: +addback(Ridge_line* ridge_line, Halfedge_handle he, Ridge_type r_type) +{ + Halfedge_handle he_cur = ( --(ridge_line->line()->end()) )->first; + FT coord_cur = ( --(ridge_line->line()->end()) )->second;//bary_coord(he_cur); + FT coord = bary_coord(he,r_type); + 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 + 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)); + + FT k1x, k2x; //abs value of the ppal curvatures at the Xing point on he. + FT k_second; // abs value of the second derivative of the curvature + // along the line of curvature + k1x = CGAL::abs(k1[v_p]) * coord + CGAL::abs(k1[v_q]) * (1-coord) ; + k2x = CGAL::abs(k2[v_p]) * coord + CGAL::abs(k2[v_q]) * (1-coord) ; + + if ( (ridge_line->line_type() == BLUE_ELLIPTIC_RIDGE) + || (ridge_line->line_type() == BLUE_HYPERBOLIC_RIDGE) + || (ridge_line->line_type() == BLUE_CREST) ) { + ridge_line->strength() += k1x * CGAL::sqrt(segment * segment); + k_second =CGAL::abs(( CGAL::abs(P1[v_p]) * coord + CGAL::abs(P1[v_q]) * (1-coord) )/(k1x-k2x)); + ridge_line->sharpness() += k_second * CGAL::sqrt(segment * segment) * model_size; + } + if ( (ridge_line->line_type() == RED_ELLIPTIC_RIDGE) + || (ridge_line->line_type() == RED_HYPERBOLIC_RIDGE) + || (ridge_line->line_type() == RED_CREST) ) { + ridge_line->strength() += k2x * CGAL::sqrt(segment * segment); + k_second =CGAL::abs(( CGAL::abs(P2[v_p]) * coord + CGAL::abs(P2[v_q]) * (1-coord) )/(k1x-k2x)); + ridge_line->sharpness() += k_second * CGAL::sqrt(segment * segment) * model_size; + } + ridge_line->line()->push_back( ridge_he(he, coord)); +} + + + +template < class Poly, class OutputIt, class Vertex2FTPropertyMap, class Vertex2VectorPropertyMap > + void Ridge_approximation:: +addfront(Ridge_line* ridge_line, Halfedge_handle he, Ridge_type r_type) +{ + Halfedge_handle he_cur = ( ridge_line->line()->begin() )->first; + FT coord_cur = ( ridge_line->line()->begin() )->second; + FT coord = bary_coord(he,r_type); + 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 + 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)); + + FT k1x, k2x; //abs value of the ppal curvatures at the Xing point on he. + FT k_second; // abs value of the second derivative of the curvature + // along the line of curvature + k1x = CGAL::abs(k1[v_p]) * coord + CGAL::abs(k1[v_q]) * (1-coord) ; + k2x = CGAL::abs(k2[v_p]) * coord + CGAL::abs(k2[v_q]) * (1-coord) ; + + if ( (ridge_line->line_type() == BLUE_ELLIPTIC_RIDGE) + || (ridge_line->line_type() == BLUE_HYPERBOLIC_RIDGE) + || (ridge_line->line_type() == BLUE_CREST) ) { + ridge_line->strength() += k1x * CGAL::sqrt(segment * segment); + k_second =CGAL::abs(( CGAL::abs(P1[v_p]) * coord + CGAL::abs(P1[v_q]) * (1-coord) )/(k1x-k2x)); + ridge_line->sharpness() += k_second * CGAL::sqrt(segment * segment) * model_size; + } + if ( (ridge_line->line_type() == RED_ELLIPTIC_RIDGE) + || (ridge_line->line_type() == RED_HYPERBOLIC_RIDGE) + || (ridge_line->line_type() == RED_CREST) ) { + ridge_line->strength() += k2x * CGAL::sqrt(segment * segment); + k_second =CGAL::abs(( CGAL::abs(P2[v_p]) * coord + CGAL::abs(P2[v_q]) * (1-coord) )/(k1x-k2x)); + ridge_line->sharpness() += k_second * CGAL::sqrt(segment * segment) * model_size; + } + ridge_line->line()->push_front( ridge_he(he, coord)); +} + + +template < class Poly, class OutputIt, class Vertex2FTPropertyMap, class Vertex2VectorPropertyMap > +typename Poly::Traits::FT Ridge_approximation:: +bary_coord(Halfedge_handle he, Ridge_type r_type) +{ + FT b_p, b_q; // extremalities at p and q for he: p->q + if ( (r_type == BLUE_ELLIPTIC_RIDGE) + || (r_type == BLUE_HYPERBOLIC_RIDGE) + || (r_type == BLUE_CREST) ) { + b_p = b0[he->opposite()->vertex()]; + b_q = b0[he->vertex()]; + } + if ( (r_type == RED_ELLIPTIC_RIDGE) + || (r_type == RED_HYPERBOLIC_RIDGE) + || (r_type == RED_CREST) ) { + b_p = b3[he->opposite()->vertex()]; + b_q = b3[he->vertex()]; + } + return CGAL::abs(b_q) / ( CGAL::abs(b_q) + CGAL::abs(b_p) ); +} + + CGAL_END_NAMESPACE #endif diff --git a/Ridges_3/include/CGAL/Umbilic.h b/Ridges_3/include/CGAL/Umbilic.h index 550e8bce859..2268786d43c 100644 --- a/Ridges_3/include/CGAL/Umbilic.h +++ b/Ridges_3/include/CGAL/Umbilic.h @@ -3,11 +3,15 @@ #include #include +#include CGAL_BEGIN_NAMESPACE enum Umbilic_type { NON_GENERIC = 0, WEDGE, TRISECTOR}; +//------------------------------------------------------------------- +// Umbilic : stores umbilic data +//------------------------------------------------------------------ template < class Poly > class Umbilic { @@ -22,46 +26,136 @@ public: std::list contour; //contructor Umbilic(Vertex_handle v_init, - std::list contour_init); + std::list contour_init) + : v(v_init), contour(contour_init) {}; +}; + + +//--------------------------------------------------------------------------- +//Umbilic_approximation +//-------------------------------------------------------------------------- +template < class Poly, class OutputIt, class Vertex2FTPropertyMap, class Vertex2VectorPropertyMap > + class Umbilic_approximation +{ +public: + typedef typename Poly::Traits::FT FT; + typedef typename Poly::Traits::Vector_3 Vector_3; + typedef typename Poly::Vertex_handle Vertex_handle; + typedef typename Poly::Halfedge_handle Halfedge_handle; + typedef typename Poly::Facet_handle Facet_handle; + typedef typename Poly::Facet_iterator Facet_iterator; + typedef typename Poly::Vertex_iterator Vertex_iterator; + typedef Umbilic Umbilic; + static FT neigh_size;//the size of neighbourhood for umbilic + // computation is (neigh_size * OneRingSize) + + Umbilic_approximation(Poly &P, + Vertex2FTPropertyMap vertex2k1_pm, Vertex2FTPropertyMap vertex2k2_pm, + Vertex2VectorPropertyMap vertex2d1_pm, Vertex2VectorPropertyMap vertex2d2_pm); + + OutputIt compute(Poly &P, OutputIt it, FT size); + + protected: + typedef T_PolyhedralSurf_neighbors Poly_neighbors; + Poly_neighbors* poly_neighbors; + + CGAL::Abs cgal_abs; + + //Property maps + Vertex2FTPropertyMap k1, k2; + Vertex2VectorPropertyMap d1, d2; + // index: following CW the contour, we choose an orientation for the // max dir of an arbitrary starting point, the max dir field is // oriented on the next point so that the scalar product of the // consecutive vectors is positive. Adding all the angles between // consecutive vectors around the contour gives -/+180 for a // wedge/trisector - void compute_type(); + void compute_type(Umbilic& umb); }; -//------------------------------------------------------------------- -// Implementation -//------------------------------------------------------------------ +template < class Poly, class OutputIt, class Vertex2FTPropertyMap, class Vertex2VectorPropertyMap > + Umbilic_approximation< Poly, OutputIt, Vertex2FTPropertyMap, Vertex2VectorPropertyMap >:: +Umbilic_approximation(Poly &P, + Vertex2FTPropertyMap vertex2k1_pm, Vertex2FTPropertyMap vertex2k2_pm, + Vertex2VectorPropertyMap vertex2d1_pm, Vertex2VectorPropertyMap vertex2d2_pm) + : k1(vertex2k1_pm), k2(vertex2k2_pm), + d1(vertex2d1_pm), d2(vertex2d2_pm) +{ + poly_neighbors = new Poly_neighbors(P); +} -template < class Poly > -Umbilic:: -Umbilic(Vertex_handle v_init, - std::list contour_init) - : v(v_init), contour(contour_init) -{} +template < class Poly, class OutputIt, class Vertex2FTPropertyMap, class Vertex2VectorPropertyMap > + OutputIt Umbilic_approximation< Poly, OutputIt, Vertex2FTPropertyMap, Vertex2VectorPropertyMap >:: +compute(Poly &P, OutputIt umbilics_it, FT size) +{ + std::vector vces; + std::list contour; + double umbilicEstimatorVertex, umbilicEstimatorNeigh; + + bool is_umbilic = true; -template < class Poly > -void Umbilic:: -compute_type() + //MAIN loop on P vertices + Vertex_iterator itb = P.vertices_begin(), ite = P.vertices_end(); + for (;itb != ite; itb++) { + Vertex_handle vh = itb; + umbilicEstimatorVertex = cgal_abs(k1[vh]-k2[vh]); + //reset vector, list and bool + vces.clear(); + contour.clear(); + is_umbilic = true; + poly_neighbors->compute_neighbors(vh, vces, contour, size); + + + // OPTIONAL: avoid umbilics whose contours touch the border + typename std::list::iterator itb_cont = contour.begin(), + ite_cont = contour.end(); + for (; itb_cont != ite_cont; itb_cont++) + if ( (*itb_cont)->is_border() ) {is_umbilic = false; continue;} + if (is_umbilic == false) continue; + + //is v an umbilic? + //a priori is_umbilic = true, and it switches to false as soon as a + // neigh vertex has a lower umbilicEstimator value + typename std::vector::iterator itbv = vces.begin(), + itev = vces.end(); + assert(*itbv == vh); + itbv++; + for (; itbv != itev; itbv++) + { umbilicEstimatorNeigh = cgal_abs( k1[*itbv] - k2[*itbv] ); + if ( umbilicEstimatorNeigh < umbilicEstimatorVertex ) + {is_umbilic = false; break;} + } + if (is_umbilic == false) continue; + + //v is an umbilic, compute the index + Umbilic* cur_umbilic = new Umbilic(vh, contour); + compute_type(*cur_umbilic); + *umbilics_it++ = cur_umbilic; + } + return umbilics_it; +} + +template < class Poly, class OutputIt, class Vertex2FTPropertyMap, class Vertex2VectorPropertyMap > + void Umbilic_approximation< Poly, OutputIt, Vertex2FTPropertyMap, Vertex2VectorPropertyMap >:: +compute_type(Umbilic& umb) { Vector_3 dir, dirnext, normal; - double cosinus, angle=0, pi=3.141592653589793, angleSum=0; + double cosinus, angle=0, angleSum=0; + const double pi=3.141592653589793; Vertex_handle v; - typename std::list::iterator itb = contour.begin(), - itlast = --contour.end(); + typename std::list::iterator itb = umb.contour.begin(), + itlast = --umb.contour.end(); v = (*itb)->vertex(); - dir = v->d1(); - normal = CGAL::cross_product(v->d1(), v->d2()); + dir = d1[v]; + normal = CGAL::cross_product(d1[v], d2[v]); //sum angles along the contour do{ itb++; v=(*itb)->vertex(); - dirnext = v->d1(); + dirnext = d1[v]; cosinus = dir*dirnext; if (cosinus < 0) {dirnext = dirnext*(-1); cosinus *= -1;} if (cosinus>1) cosinus = 1; @@ -70,13 +164,13 @@ compute_type() else angle = -acos(cosinus); angleSum += angle; dir = dirnext; - normal = CGAL::cross_product(v->d1(), v->d2()); + normal = CGAL::cross_product(d1[v], d2[v]); } while (itb != (itlast)); //angle (v_last, v_0) - v=(*contour.begin())->vertex(); - dirnext = v->d1(); + v=(*umb.contour.begin())->vertex(); + dirnext = d1[v]; cosinus = dir*dirnext; if (cosinus < 0) {dirnext = dirnext*(-1); cosinus *= -1;} if (cosinus>1) cosinus = 1; @@ -84,9 +178,9 @@ compute_type() else angle = -acos(cosinus); angleSum += angle; - if ((angleSum > (pi/2)) && (angleSum < (3*pi/2))) umb_type = TRISECTOR ; - else if ((angleSum < (-pi/2)) && (angleSum > (-3*pi/2))) umb_type = WEDGE; - else umb_type = NON_GENERIC; + if ((angleSum > (pi/2)) && (angleSum < (3*pi/2))) umb.umb_type = TRISECTOR ; + else if ((angleSum < (-pi/2)) && (angleSum > (-3*pi/2))) umb.umb_type = WEDGE; + else umb.umb_type = NON_GENERIC; } CGAL_END_NAMESPACE