diff --git a/Packages/Stream_lines_2/include/CGAL/Euler_integrator_2.h b/Packages/Stream_lines_2/include/CGAL/Euler_integrator_2.h index a279df47df7..c972d03eea5 100644 --- a/Packages/Stream_lines_2/include/CGAL/Euler_integrator_2.h +++ b/Packages/Stream_lines_2/include/CGAL/Euler_integrator_2.h @@ -27,9 +27,10 @@ CGAL_BEGIN_NAMESPACE // The class Euler_integrator_2 is a model of the concept Integrator template -class Euler_integrator_2{ +class Euler_integrator_2 +{ public: - typedef Euler_integrator_2 self; + typedef Euler_integrator_2 self; typedef typename VectorField_2::FT FT; typedef typename VectorField_2::Point_2 Point_2; typedef typename VectorField_2::Vector_2 Vector_2; @@ -37,17 +38,32 @@ public: protected: FT default_integration_step; public: - Euler_integrator_2(); + Euler_integrator_2(); + Euler_integrator_2(const FT & integration_step); - inline Point_2 operator()(const Point_2 & p, const Vector_field_2 & vector_field_2, const bool & index) const; - inline Point_2 operator()(const Point_2 & p, const Vector_field_2 & vector_field_2, const FT & integration_step, const bool & index) const; - inline Point_2 operator()(const Point_2 & p, const Vector_field_2 & vector_field_2, const FT & integration_step, Vector_2 v, const bool & index) const; - inline FT get_default_integration_step(){return default_integration_step;} -// just for debugging -inline FT distance(const Point_2 & p, const Point_2 & q) -{ - return sqrt(((p.x() - q.x())*(p.x() - q.x()))+((p.y() - q.y())*(p.y() - q.y()))); -} + + inline Point_2 operator()(const Point_2 & p, const Vector_field_2 & + vector_field_2, const bool & index) const; + + inline Point_2 operator()(const Point_2 & p, const Vector_field_2 & + vector_field_2, const FT & + integration_step, const bool & index) + const; + + inline Point_2 operator()(const Point_2 & p, const Vector_field_2 & + vector_field_2, const FT & + integration_step, Vector_2 v, const bool & + index) const; + + inline FT get_default_integration_step() + { + return default_integration_step; + } + // just for debugging + inline FT distance(const Point_2 & p, const Point_2 & q) + { + return sqrt(((p.x() - q.x())*(p.x() - q.x()))+((p.y() - q.y())*(p.y() - q.y()))); + } }; template diff --git a/Packages/Stream_lines_2/include/CGAL/Regular_grid_2.h b/Packages/Stream_lines_2/include/CGAL/Regular_grid_2.h index 6f7b43839dd..39bdfff1210 100644 --- a/Packages/Stream_lines_2/include/CGAL/Regular_grid_2.h +++ b/Packages/Stream_lines_2/include/CGAL/Regular_grid_2.h @@ -31,9 +31,10 @@ CGAL_BEGIN_NAMESPACE // bilinear interpolation to extract the vector field values template -class Regular_grid_2{ +class Regular_grid_2 +{ public: - typedef Regular_grid_2 Vector_field_2; + typedef Regular_grid_2 Vector_field_2; typedef StreamLinesTraits_2 Geom_traits; typedef typename StreamLinesTraits_2::FT FT; typedef typename StreamLinesTraits_2::Point_2 Point_2; @@ -50,68 +51,107 @@ protected: bool is_in_samples(const int & i,const int & j) const; public: Regular_grid_2(const int & m, const int & n,const FT & x, const FT & y); -// Regular_grid_2(); - ~Regular_grid_2(){delete [] vector_field;} - void fill(std::ifstream & f, const FT & x, const FT & y); + // Regular_grid_2(); + ~Regular_grid_2() + { + delete [] vector_field; + } + void fill(std::ifstream & f, const FT & x, const FT & y); inline typename Geom_traits::Iso_rectangle_2 bbox() const; std::pair get_field(const Point_2 & p) const - { - CGAL_streamlines_precondition(is_in_domain(p)); - Vector_2 v = get_vector_field(p); - FT density = get_density_field(p); - return std::pair(v,density); - } + { + CGAL_streamlines_precondition(is_in_domain(p)); + Vector_2 v = get_vector_field(p); + FT density = get_density_field(p); + return std::pair(v,density); + } inline bool is_in_domain(const Point_2 & p) const; inline FT get_integration_step(const Point_2 &) const; - inline FT get_integration_step() const; + inline FT get_integration_step() const; inline void set_field(const int & i, const int & j, const Vector_2 & v); - inline std::pair get_dimension(){return std::pair(number_of_samples_x, number_of_samples_y);} - inline std::pair get_size(){return std::pair(domain_size_x, domain_size_y);}}; + inline std::pair get_dimension() + { + return std::pair(number_of_samples_x, number_of_samples_y); + } + inline std::pair get_size() + { + return std::pair(domain_size_x, + domain_size_y); + } +}; template inline typename Regular_grid_2::Geom_traits::Iso_rectangle_2 -Regular_grid_2::bbox() const{ - return typename Geom_traits::Iso_rectangle_2(0.0, 0.0, domain_size_x, domain_size_y);} +Regular_grid_2::bbox() const +{ + return typename Geom_traits::Iso_rectangle_2(0.0, 0.0, + domain_size_x, + domain_size_y); +} template inline int -Regular_grid_2::get_index(const int & i, const int & j) const{ - return 2*(number_of_samples_x*j + i);} +Regular_grid_2::get_index(const int & i, const + int & j) const +{ + return 2*(number_of_samples_x*j + i); +} template -Regular_grid_2::Regular_grid_2(const int & m, const int & n,const FT & x, const FT & y){ +Regular_grid_2::Regular_grid_2(const int & m, + const int & + n,const FT & x, + const FT & y) +{ number_of_samples_x = m; number_of_samples_y = n; domain_size_x = x; domain_size_y = y; - vector_field = new FT[number_of_samples_x*number_of_samples_y* 2];} + vector_field = new FT[number_of_samples_x*number_of_samples_y* 2]; +} template inline void -Regular_grid_2::set_field(const int & i, const int & j, const Vector_2 & v){ +Regular_grid_2::set_field(const int & i, const + int & j, const Vector_2 + & v) +{ CGAL_streamlines_precondition(is_in_samples(i,j)); int index = get_index(i,j); vector_field[index++] = v.x(); - vector_field[index] = v.y();} + vector_field[index] = v.y(); +} template inline bool -Regular_grid_2::is_in_domain(const Point_2 & p) const{ - return ((p.x()>=0.0) && (p.x()<=domain_size_x) && (p.y()>=0.0) && (p.y()<=domain_size_y));} +Regular_grid_2::is_in_domain(const Point_2 & p) + const +{ + return ((p.x()>=0.0) && (p.x()<=domain_size_x) && (p.y()>=0.0) && + (p.y()<=domain_size_y)); +} template bool -Regular_grid_2::is_in_samples(const int & i, const int & j) const{ - return ((i>=0) && (i<=number_of_samples_x-1) && (j>=0) && (j<=number_of_samples_y-1));} +Regular_grid_2::is_in_samples(const int & i, + const int & j) + const +{ + return ((i>=0) && (i<=number_of_samples_x-1) && (j>=0) && + (j<=number_of_samples_y-1)); +} template typename Regular_grid_2::Vector_2 -Regular_grid_2::get_vector_field(const Point_2 & p) const{ +Regular_grid_2::get_vector_field(const Point_2 & + p) const +{ FT fXv,fYv; FT x = (p.x() / domain_size_x) * (number_of_samples_x-1); FT y = (p.y() / domain_size_y) * (number_of_samples_y-1); @@ -146,24 +186,34 @@ Regular_grid_2::get_vector_field(const Point_2 & p) const{ fXv = fXv / normal; fYv = fYv / normal; Vector_2 v = Vector_2(fXv, fYv); - return v;} + return v; +} template typename Regular_grid_2::FT -Regular_grid_2::get_density_field(const Point_2 & p) const{ - return 1.0;} +Regular_grid_2::get_density_field(const Point_2 & + p) const +{ + return 1.0; +} template inline typename Regular_grid_2::FT -Regular_grid_2::get_integration_step(const Point_2 &) const{ - return 1.0;} +Regular_grid_2::get_integration_step(const + Point_2 &) + const +{ + return 1.0; +} template inline typename Regular_grid_2::FT -Regular_grid_2::get_integration_step() const{ - return 1.0;} +Regular_grid_2::get_integration_step() const +{ + return 1.0; +} CGAL_END_NAMESPACE diff --git a/Packages/Stream_lines_2/include/CGAL/Runge_kutta_integrator_2.h b/Packages/Stream_lines_2/include/CGAL/Runge_kutta_integrator_2.h index 5b94895e4e8..24dbc712c9b 100644 --- a/Packages/Stream_lines_2/include/CGAL/Runge_kutta_integrator_2.h +++ b/Packages/Stream_lines_2/include/CGAL/Runge_kutta_integrator_2.h @@ -28,36 +28,55 @@ CGAL_BEGIN_NAMESPACE // The class Runge_kutta_integrator_2 is a model of the concept Integrator template -class Runge_kutta_integrator_2{ +class Runge_kutta_integrator_2 +{ public: - typedef Runge_kutta_integrator_2 self; + typedef Runge_kutta_integrator_2 self; typedef typename VectorField_2::FT FT; typedef typename VectorField_2::Point_2 Point_2; typedef typename VectorField_2::Vector_2 Vector_2; typedef typename VectorField_2::Vector_field_2 Vector_field_2; protected: - typedef CGAL::Euler_integrator_2 Euler_integrator_2; + typedef CGAL::Euler_integrator_2 Euler_integrator_2; + Euler_integrator_2 * euler_integrator_2; + FT default_integration_step; + public: Runge_kutta_integrator_2(); + Runge_kutta_integrator_2(FT integration_step); + ~Runge_kutta_integrator_2(); + Point_2 operator()(const Point_2 & p, const Vector_field_2 & vector_field_2, const bool & index) const; + Point_2 operator()(const Point_2 & p, const Vector_field_2 & vector_field_2, const FT & integration_step, const bool & index) const; + Point_2 operator()(const Point_2 & p, const Vector_field_2 & vector_field_2, const FT & integration_step, Vector_2 v, const bool & index) const; - inline FT get_default_integration_step(){return default_integration_step;} - Euler_integrator_2 * get_euler_integrator_2(){return euler_integrator_2;} - inline FT distance(const Point_2 & p, const Point_2 & q){ - return sqrt(((p.x() - q.x())*(p.x() - q.x()))+((p.y() - q.y())*(p.y() - q.y())));} + inline FT get_default_integration_step() + { + return default_integration_step; + } + Euler_integrator_2 * get_euler_integrator_2() + { + return euler_integrator_2; + } + inline FT distance(const Point_2 & p, const Point_2 & q) + { + return sqrt(((p.x() - q.x())*(p.x() - q.x()))+((p.y() + - + q.y())*(p.y() - q.y()))); + } }; template Runge_kutta_integrator_2:: Runge_kutta_integrator_2() { - euler_integrator_2 = new Euler_integrator_2(); + euler_integrator_2 = new Euler_integrator_2(); default_integration_step = euler_integrator_2->get_default_integration_step(); } @@ -66,7 +85,7 @@ template Runge_kutta_integrator_2:: Runge_kutta_integrator_2(FT integration_step) { - euler_integrator_2 = new Euler_integrator_2(integration_step); + euler_integrator_2 = new Euler_integrator_2(integration_step); default_integration_step = euler_integrator_2->get_default_integration_step(); } @@ -75,7 +94,7 @@ template Runge_kutta_integrator_2:: ~Runge_kutta_integrator_2() { - delete euler_integrator_2; + delete euler_integrator_2; } template @@ -95,7 +114,9 @@ Point_2 Runge_kutta_integrator_2::operator() template inline typename Runge_kutta_integrator_2::Point_2 Runge_kutta_integrator_2::operator() -(const Point_2 & p, const Vector_field_2& vector_field_2, const FT & integration_step, const bool & index) const{ +(const Point_2 & p, const Vector_field_2& vector_field_2, const FT & + integration_step, const bool & index) const +{ Vector_2 v; v = vector_field_2.get_field(p).first; Runge_kutta_integrator_2 diff --git a/Packages/Stream_lines_2/include/CGAL/Stream_lines_2.h b/Packages/Stream_lines_2/include/CGAL/Stream_lines_2.h index 5a3a2cf0bda..bb552003c0c 100644 --- a/Packages/Stream_lines_2/include/CGAL/Stream_lines_2.h +++ b/Packages/Stream_lines_2/include/CGAL/Stream_lines_2.h @@ -39,208 +39,222 @@ CGAL_BEGIN_NAMESPACE template -class Stream_lines_2{ +class Stream_lines_2 +{ public: - typedef typename VectorField_2::Vector_field_2 Vector_field_2; - typedef typename VectorField_2::Geom_traits Geom_traits; - typedef typename VectorField_2::FT FT; - typedef typename VectorField_2::Point_2 Point_2; - typedef typename VectorField_2::Vector_2 Vector_2; + typedef typename VectorField_2::Vector_field_2 Vector_field_2; + typedef typename VectorField_2::Geom_traits Geom_traits; + typedef typename VectorField_2::FT FT; + typedef typename VectorField_2::Point_2 Point_2; + typedef typename VectorField_2::Vector_2 Vector_2; protected: -// typedef CGAL::Cartesian K1; -// typedef struct Kernel : public K1 {}; -// typedef CGAL::Filtered_kernel Kernel; - typedef Geom_traits Kernel; - typedef CGAL::Triangulation_vertex_base_2 Vb; - typedef CGAL::Triangulation_face_base_2 Fb; - typedef CGAL::Triangulation_data_structure_2 TDS; - typedef CGAL::Delaunay_triangulation_2 DT; - typedef typename DT::Vertex_handle Vertex_handle; - typedef typename DT::Face_handle Face_handle; - typedef typename DT::Face_circulator Face_circulator; - typedef typename DT::Edge_iterator Edge_iterator; - typedef std::pair Circle; - typedef - CGAL::Quadruple Pq_element; - Pq_element Biggest_circle; - FT distance(Point_2 p, Point_2 q){ - return sqrt(((p.x() - q.x())*(p.x() - q.x()))+((p.y() - q.y())*(p.y() - q.y())));} - int ir; - int il; - Pq_element Pq_element_max_r; - Pq_element Pq_previous_r,Pq_current_r,Pq_next_r; - Pq_element Pq_element_max_l; - Pq_element Pq_previous_l,Pq_current_l,Pq_next_l; + typedef Geom_traits Kernel; + typedef CGAL::Triangulation_vertex_base_2 Vb; + typedef CGAL::Triangulation_face_base_2 Fb; + typedef CGAL::Triangulation_data_structure_2 TDS; + typedef CGAL::Delaunay_triangulation_2 DT; + typedef typename DT::Vertex_handle Vertex_handle; + typedef typename DT::Face_handle Face_handle; + typedef typename DT::Face_circulator Face_circulator; + typedef typename DT::Edge_iterator Edge_iterator; + typedef std::pair Circle; + typedef + CGAL::Quadruple Pq_element; + Pq_element Biggest_circle; + FT distance(Point_2 p, Point_2 q) + { + return sqrt(((p.x() - q.x())*(p.x() - q.x()))+((p.y() + - + q.y())*(p.y() - q.y()))); + } + int ir; + int il; + Pq_element Pq_element_max_r; + Pq_element Pq_previous_r,Pq_current_r,Pq_next_r; + Pq_element Pq_element_max_l; + Pq_element Pq_previous_l,Pq_current_l,Pq_next_l; public: - DT m_DT; - typedef std::list Point_container_2; - typedef typename Point_container_2::iterator Point_iterator_2; - typedef std::list > Iterator_container_2; - Iterator_container_2 iterator_container; - typedef std::list Vertex_container_2; - typedef typename Iterator_container_2::iterator Stream_line_iterator_2; - typedef std::list Stream_line_container_2; + DT m_DT; + typedef std::list Point_container_2; + typedef typename Point_container_2::iterator Point_iterator_2; + typedef std::list > Iterator_container_2; + Iterator_container_2 iterator_container; + typedef std::list Vertex_container_2; + typedef typename Iterator_container_2::iterator Stream_line_iterator_2; + typedef std::list Stream_line_container_2; protected: - Stream_line_container_2 stl_container; - class C{ - public: - bool operator()(const Pq_element &a1, const Pq_element &a2){ - return a1.fourth.second < a2.fourth.second ;}}; - std::priority_queue, C> pq; - int iOrder_insertion; - FT fSepStl_seed; - FT separating_distance; - FT saturation_ratio; - Point_2 seed_point; - unsigned int _number_of_lines; + Stream_line_container_2 stl_container; + class C + { + public: + bool operator()(const Pq_element &a1, const Pq_element + &a2) + { + return a1.fourth.second < a2.fourth.second ; + } + }; + std::priority_queue, C> pq; + int iOrder_insertion; + FT fSepStl_seed; + FT separating_distance; + FT saturation_ratio; + Point_2 seed_point; + unsigned int _number_of_lines; protected: - void place_stream_lines(const Vector_field_2 & vector_field_2, const Integrator_2 & integrator, - const int & sampling_step, const bool & step_by_step = false); - bool get_next_seed_point(FT & distanceg, Point_2 & seed_point); - FT find_smallest_circle(const Vertex_handle & pVertex_handle); - Vertex_handle insert_point(const Point_2 & pPoint, FT& fDist,bool bDistanceCalculation); - Vertex_handle insert_point(const Point_2 & pPoint, const Face_handle & m_Face_handle, FT& fDist, - bool bDistanceCalculation); - void integrate_streamline(const Vector_field_2 & vector_field_2, const Integrator_2 & integrator, - Point_container_2& stl, Point_2& seed_point, Vertex_container_2& stl_vertices, const int & sampling_step); - void integrate_forward(const Vector_field_2 & vector_field_2, const Integrator_2 & integrator, - Point_container_2& stl,Point_2& seed_point, - Vertex_container_2& stl_vertices, const int & sampling_step); - void integrate_backward(const Vector_field_2 & vector_field_2, const Integrator_2 & integrator, Point_container_2& - stl, Vertex_container_2& stl_vertices, const int & sampling_step); - void insert_streamline(const Vector_field_2 & vector_field_2, Point_container_2 stl, + void place_stream_lines(const Vector_field_2 & vector_field_2, const Integrator_2 & integrator, + const int & sampling_step, const bool & step_by_step = false); + bool get_next_seed_point(FT & distanceg, Point_2 & seed_point); + FT find_smallest_circle(const Vertex_handle & pVertex_handle); + Vertex_handle insert_point(const Point_2 & pPoint, FT& fDist,bool bDistanceCalculation); + Vertex_handle insert_point(const Point_2 & pPoint, const Face_handle & m_Face_handle, FT& fDist, + bool bDistanceCalculation); + void integrate_streamline(const Vector_field_2 & vector_field_2, const Integrator_2 & integrator, + Point_container_2& stl, Point_2& seed_point, Vertex_container_2& stl_vertices, const int & sampling_step); + void integrate_forward(const Vector_field_2 & vector_field_2, const Integrator_2 & integrator, + Point_container_2& stl,Point_2& seed_point, + Vertex_container_2& stl_vertices, const int & sampling_step); + void integrate_backward(const Vector_field_2 & vector_field_2, const Integrator_2 & integrator, Point_container_2& + stl, Vertex_container_2& stl_vertices, const int & sampling_step); + void insert_streamline(const Vector_field_2 & vector_field_2, Point_container_2 stl, Vertex_container_2 stl_vertices); - void pq_elements(const Vector_field_2 & vector_field_2, Vertex_container_2 stl_vertices, int i, - const Vertex_handle & m_Vertex_handle, int before_end); - void make_iterator(); + void pq_elements(const Vector_field_2 & vector_field_2, Vertex_container_2 stl_vertices, int i, + const Vertex_handle & m_Vertex_handle, int before_end); + void make_iterator(); public: - Stream_lines_2(const Vector_field_2 & m_vector_field_2, const Integrator_2 & m_integrator, const FT - & m_separating_distance, const FT & m_saturation_ratio, const int & sampling_insertion = 0, const bool & - step_by_step = false); - bool continue_next(const Vector_field_2 & vector_field_2, const Integrator_2 & integrator, const int & sampling_step); - Stream_line_iterator_2 begin(); - Stream_line_iterator_2 end(); - // for visualizing streamlines - void print_stream_lines(std::ofstream & fw); - void print_stream_lines_eps(std::ofstream & fw); - std::list get_pq(); - unsigned int number_of_lines(){return _number_of_lines;} - std::list< std::pair > get_tr() + Stream_lines_2(const Vector_field_2 & m_vector_field_2, const Integrator_2 & m_integrator, const FT + & m_separating_distance, const FT & m_saturation_ratio, const int & sampling_insertion = 0, const bool & + step_by_step = false); + bool continue_next(const Vector_field_2 & vector_field_2, const Integrator_2 & integrator, const int & sampling_step); + Stream_line_iterator_2 begin(); + Stream_line_iterator_2 end(); + // for visualizing streamlines + void print_stream_lines(std::ofstream & fw); + void print_stream_lines_eps(std::ofstream & fw); + std::list get_pq(); + unsigned int number_of_lines() { return _number_of_lines; } + std::list< std::pair > get_tr() + { + std::list< std::pair > _list; + Edge_iterator eit = m_DT.edges_begin(); + for (;eit != m_DT.edges_end();eit++) { - std::list< std::pair > _list; - Edge_iterator eit = m_DT.edges_begin(); - for (;eit != m_DT.edges_end();eit++) - { - Point_2 p1 = (*eit).first->vertex(m_DT.ccw((*eit).second))->point(); - Point_2 p2 = (*eit).first->vertex(m_DT.cw((*eit).second))->point(); - _list.push_front(std::pair(p1, p2)); - } - return _list; - } - std::pair get_biggest_circle() - { - Pq_element m_Pq = Biggest_circle; - std::pair circle(m_Pq.fourth.first, m_Pq.fourth.second); - return circle; + Point_2 p1 = (*eit).first->vertex(m_DT.ccw((*eit).second))->point(); + Point_2 p2 = (*eit).first->vertex(m_DT.cw((*eit).second))->point(); + _list.push_front(std::pair(p1, p2)); } + return _list; + } + std::pair get_biggest_circle() + { + Pq_element m_Pq = Biggest_circle; + std::pair circle(m_Pq.fourth.first, m_Pq.fourth.second); + return circle; + } protected: - FT max_x; - FT min_x; - FT max_y; - FT min_y; + FT max_x; + FT min_x; + FT max_y; + FT min_y; protected: - Stream_line_iterator_2 begin_iterator; - Stream_line_iterator_2 end_iterator; + Stream_line_iterator_2 begin_iterator; + Stream_line_iterator_2 end_iterator; private: - int number_of_points; + int number_of_points; }; template Stream_lines_2::Stream_lines_2(const Vector_field_2 & vector_field_2, const Integrator_2 & m_integrator, const FT & m_separating_distance, const FT -& m_saturation_ratio, const int & sampling_step, const bool & step_by_step){ - separating_distance = m_separating_distance; - saturation_ratio = m_saturation_ratio; - ir = il = 0; // initialization - fSepStl_seed = separating_distance*saturation_ratio; - max_x = vector_field_2.bbox().xmax(); - min_x = vector_field_2.bbox().xmin(); - max_y = vector_field_2.bbox().ymax(); - min_y = vector_field_2.bbox().ymin(); - m_DT.clear(); - Point_2 pPoint; - pPoint = Point_2(min_x-separating_distance,min_y-separating_distance); - m_DT.insert(pPoint); - pPoint = Point_2(min_x-separating_distance,max_y+separating_distance); - m_DT.insert(pPoint); - pPoint = Point_2(max_x+separating_distance,min_y-separating_distance); - m_DT.insert(pPoint); - pPoint = Point_2(max_x+separating_distance,max_y+separating_distance); - m_DT.insert(pPoint); - for (int i=(int) (min_x-separating_distance);i void Stream_lines_2::place_stream_lines(const Vector_field_2 & vector_field_2, const Integrator_2 & integrator, const int & sampling_step, const bool & step_by_step) { - seed_point = Point_2((max_x+min_x)/2.0,(max_y+min_y)/2.0); -// the first chosen point can be not valid - FT xrange = max_x - min_x; - FT yrange = max_y - min_y; - while(!vector_field_2.is_in_domain(seed_point)) - { - std::cout << "searching valid seed point..\n"; - FT x = min_x + (FT) (((FT) rand() * xrange)/((FT) RAND_MAX)); - FT y = min_y + (FT) (((FT) rand() * yrange)/((FT) RAND_MAX)); - seed_point = Point_2(x, y); - } - std::cout << seed_point << " first seed point\n"; - std::cout << "creating the placement..\n"; - FT distance = (FT) max_x * (1.0/2.0); - bool b = (distance>fSepStl_seed); -// int i=0; - if (!step_by_step) - while(b) - { - Point_container_2 stl; - Vertex_container_2 stl_vertices; - integrate_streamline(vector_field_2, integrator, stl, seed_point, stl_vertices, sampling_step); - insert_streamline(vector_field_2, stl, stl_vertices); - _number_of_lines++; - b = get_next_seed_point(distance,seed_point); - } - else - { - Point_container_2 stl; - Vertex_container_2 stl_vertices; - integrate_streamline(vector_field_2, integrator, stl, seed_point, stl_vertices, sampling_step); - insert_streamline(vector_field_2, stl, stl_vertices); - _number_of_lines++; - b = get_next_seed_point(distance,seed_point); - } - make_iterator(); -} - -template -bool Stream_lines_2::continue_next(const Vector_field_2 & vector_field_2, const Integrator_2 & integrator, const int & sampling_step) -{ - FT distance; + seed_point = Point_2((max_x+min_x)/2.0,(max_y+min_y)/2.0); + // the first chosen point can be not valid + FT xrange = max_x - min_x; + FT yrange = max_y - min_y; + while(!vector_field_2.is_in_domain(seed_point)) + { + std::cout << "searching valid seed point..\n"; + FT x = min_x + (FT) (((FT) rand() * xrange)/((FT) RAND_MAX)); + FT y = min_y + (FT) (((FT) rand() * yrange)/((FT) RAND_MAX)); + seed_point = Point_2(x, y); + } + std::cout << seed_point << " first seed point\n"; + std::cout << "creating the placement..\n"; + FT distance = (FT) max_x * (1.0/2.0); + bool b = (distance>fSepStl_seed); + // int i=0; + if (!step_by_step) + while(b) + { Point_container_2 stl; Vertex_container_2 stl_vertices; integrate_streamline(vector_field_2, integrator, stl, seed_point, stl_vertices, sampling_step); insert_streamline(vector_field_2, stl, stl_vertices); _number_of_lines++; - make_iterator(); - return get_next_seed_point(distance,seed_point); + b = get_next_seed_point(distance,seed_point); + } + else + { + Point_container_2 stl; + Vertex_container_2 stl_vertices; + integrate_streamline(vector_field_2, integrator, stl, seed_point, stl_vertices, sampling_step); + insert_streamline(vector_field_2, stl, stl_vertices); + _number_of_lines++; + b = get_next_seed_point(distance,seed_point); + } + make_iterator(); +} + +template +bool Stream_lines_2::continue_next(const Vector_field_2 & vector_field_2, const Integrator_2 & integrator, const int & sampling_step) +{ + FT distance; + Point_container_2 stl; + Vertex_container_2 stl_vertices; + integrate_streamline(vector_field_2, integrator, stl, seed_point, stl_vertices, sampling_step); + insert_streamline(vector_field_2, stl, stl_vertices); + _number_of_lines++; + make_iterator(); + return get_next_seed_point(distance,seed_point); } // get the next seed point @@ -258,23 +272,23 @@ inline typename Stream_lines_2::FT Stream_lines_2::find_smallest_circle(const Vertex_handle & pVertex_handle) { - Face_circulator pFace_handle = m_DT.incident_faces(pVertex_handle); - Face_circulator pEnd = pFace_handle; - FT fMin = max_x; - CGAL_For_all(pFace_handle,pEnd) + Face_circulator pFace_handle = m_DT.incident_faces(pVertex_handle); + Face_circulator pEnd = pFace_handle; + FT fMin = max_x; + CGAL_For_all(pFace_handle,pEnd) + { + FT fDist = + CGAL::squared_radius( + pFace_handle->vertex(0)->point(), + pFace_handle->vertex(1)->point(), + pFace_handle->vertex(2)->point()) * 4.0; + fDist = sqrt(fDist); + if (fDist < fMin) { - FT fDist = - CGAL::squared_radius( - pFace_handle->vertex(0)->point(), - pFace_handle->vertex(1)->point(), - pFace_handle->vertex(2)->point()) * 4.0; - fDist = sqrt(fDist); - if (fDist < fMin) - { - fMin = fDist; - } + fMin = fDist; } - return fMin; + } + return fMin; } template @@ -282,12 +296,12 @@ inline typename Stream_lines_2::Vertex_handle Stream_lines_2::insert_point(const Point_2 & pPoint, FT& fDist,bool bDistanceCalculation) { - Vertex_handle pVertex_handle = m_DT.insert(pPoint); - if (bDistanceCalculation) - fDist = find_smallest_circle(pVertex_handle); - else - fDist = 0.0; - return (pVertex_handle); + Vertex_handle pVertex_handle = m_DT.insert(pPoint); + if (bDistanceCalculation) + fDist = find_smallest_circle(pVertex_handle); + else + fDist = 0.0; + return (pVertex_handle); } template @@ -295,418 +309,473 @@ inline typename Stream_lines_2::Vertex_handle Stream_lines_2::insert_point(const Point_2 & pPoint, const Face_handle & m_Face_handle, FT& fDist,bool bDistanceCalculation) { - Vertex_handle pVertex_handle = m_DT.insert(pPoint,m_Face_handle); - if (bDistanceCalculation) - fDist = find_smallest_circle(pVertex_handle); - else - fDist = 0.0; - return (pVertex_handle); + Vertex_handle pVertex_handle = m_DT.insert(pPoint,m_Face_handle); + if (bDistanceCalculation) + fDist = find_smallest_circle(pVertex_handle); + else + fDist = 0.0; + return (pVertex_handle); } template void Stream_lines_2::integrate_forward(const Vector_field_2 & vector_field_2, const Integrator_2 & integrator, Point_container_2& stl, Point_2& seed_point, Vertex_container_2& stl_vertices, const int & sampling_step) { - int sampling = 0; // sampling step; - int insertion = 0; // insertion order; - int insertion_step = 0; - Point_container_2 list_of_point; - Vertex_container_2 list_of_vertex; - number_of_points = 0; - Point_2 pPoint1; - bool bEnd = false; - FT dist; - Point_2 new_point = Point_2 (seed_point.x(), seed_point.y()); - Vertex_handle m_Vertex_handle = insert_point(new_point, dist, true); - stl_vertices.push_front(m_Vertex_handle); - stl.push_front(new_point); - number_of_points++; - Point_2 old_point; - insertion_step = (int) (((dist)-fSepStl_seed) / std::max((FT) sampling_step,vector_field_2.get_integration_step())); - if (insertion_step < 0) insertion_step = 0; - while (!bEnd) + int sampling = 0; // sampling step; + int insertion = 0; // insertion order; + int insertion_step = 0; + Point_container_2 list_of_point; + Vertex_container_2 list_of_vertex; + number_of_points = 0; + Point_2 pPoint1; + bool bEnd = false; + FT dist; + Point_2 new_point = Point_2 (seed_point.x(), seed_point.y()); + Vertex_handle m_Vertex_handle = insert_point(new_point, dist, true); + stl_vertices.push_front(m_Vertex_handle); + stl.push_front(new_point); + number_of_points++; + Point_2 old_point; + insertion_step = (int) (((dist)-fSepStl_seed) / std::max((FT) sampling_step,vector_field_2.get_integration_step())); + if (insertion_step < 0) insertion_step = 0; + while (!bEnd) + { + Point_2 ex_old_point = old_point; + old_point = new_point; + CGAL_streamlines_precondition(vector_field_2.is_in_domain(old_point)); + new_point = integrator(old_point,vector_field_2,true); + bEnd = !vector_field_2.is_in_domain(new_point); + bEnd = bEnd || (new_point == old_point);/* to review */ + if(number_of_points > 30) + bEnd = bEnd || ((distance(stl.front(), stl.back())) 30) - bEnd = bEnd || ((distance(stl.front(), stl.back()))face(), dist, true); - while ((dist <= separating_distance)&&(!list_of_point.empty())) - { - m_DT.remove(m_Vertex_handle); - for (int i=0;i<=sampling_step;i++){ - stl.pop_front(); - number_of_points--;} - new_point = list_of_point.front(); - list_of_point.pop_front(); - m_Vertex_handle = insert_point(new_point, stl_vertices.front()->face(), dist,true); - } - // adaptive insertion order coefficient - insertion_step = (int) (((dist)-fSepStl_seed) / - std::max((FT) sampling_step,vector_field_2.get_integration_step())); - if (insertion_step < 0) insertion_step = 0; - list_of_vertex.push_front(m_Vertex_handle); - (bEnd) = (((bEnd))||(distface(), dist, false); - list_of_vertex.push_front(m_Vertex_handle); - list_of_point.pop_front(); - } - while(!list_of_vertex.empty()) - { - stl_vertices.push_front(list_of_vertex.front()); - list_of_vertex.pop_front(); - } - insertion = 0; - } - sampling = 0; - } + stl.push_front(new_point); + number_of_points++; + insertion++; + list_of_point.push_front(new_point); } - else + else { - if (!list_of_point.empty()) - { - new_point = list_of_point.front(); - list_of_point.pop_front(); - Vertex_handle m_Vertex_handle = insert_point(new_point, dist, true); - while ((dist <= separating_distance)&&(!list_of_point.empty())) - { - m_DT.remove(m_Vertex_handle); - for (int i=0;i<=sampling_step;i++) - { - stl.pop_front(); - number_of_points--; - } - new_point = list_of_point.front(); - list_of_point.pop_front(); - m_Vertex_handle = insert_point(new_point, stl_vertices.front()->face(), dist, true); - } - insertion_step = (int) (((dist)-fSepStl_seed) / - std::max((FT) sampling_step,vector_field_2.get_integration_step())); - if (insertion_step < 0) insertion_step = 0; - (bEnd) = (((bEnd))||(distface(), dist, false); - list_of_vertex.push_front(m_Vertex_handle); - list_of_point.pop_front(); - } - while(!list_of_vertex.empty()) - { - stl_vertices.push_front(list_of_vertex.front()); - list_of_vertex.pop_front(); - } + stl.push_front(new_point); + number_of_points++; + list_of_point.push_front(new_point); + list_of_point.pop_front(); + m_Vertex_handle = insert_point(new_point, stl_vertices.front()->face(), dist, true); + while ((dist <= separating_distance)&&(!list_of_point.empty())) + { + m_DT.remove(m_Vertex_handle); + for (int i=0;i<=sampling_step;i++){ + stl.pop_front(); + number_of_points--;} + new_point = list_of_point.front(); + list_of_point.pop_front(); + m_Vertex_handle = insert_point(new_point, stl_vertices.front()->face(), dist,true); + } + // adaptive insertion order coefficient + insertion_step = (int) (((dist)-fSepStl_seed) / + std::max((FT) sampling_step,vector_field_2.get_integration_step())); + if (insertion_step < 0) insertion_step = 0; + list_of_vertex.push_front(m_Vertex_handle); + (bEnd) = (((bEnd))||(distface(), dist, false); + list_of_vertex.push_front(m_Vertex_handle); + list_of_point.pop_front(); + } + while(!list_of_vertex.empty()) + { + stl_vertices.push_front(list_of_vertex.front()); + list_of_vertex.pop_front(); + } + insertion = 0; } + sampling = 0; + } } + else + { + if (!list_of_point.empty()) + { + new_point = list_of_point.front(); + list_of_point.pop_front(); + Vertex_handle m_Vertex_handle = insert_point(new_point, dist, true); + while ((dist <= separating_distance)&&(!list_of_point.empty())) + { + m_DT.remove(m_Vertex_handle); + for (int i=0;i<=sampling_step;i++) + { + stl.pop_front(); + number_of_points--; + } + new_point = list_of_point.front(); + list_of_point.pop_front(); + m_Vertex_handle = insert_point(new_point, stl_vertices.front()->face(), dist, true); + } + insertion_step = (int) (((dist)-fSepStl_seed) / + std::max((FT) sampling_step,vector_field_2.get_integration_step())); + if (insertion_step < 0) insertion_step = 0; + (bEnd) = (((bEnd))||(distface(), dist, false); + list_of_vertex.push_front(m_Vertex_handle); + list_of_point.pop_front(); + } + while(!list_of_vertex.empty()) + { + stl_vertices.push_front(list_of_vertex.front()); + list_of_vertex.pop_front(); + } + } + } } template void Stream_lines_2::integrate_backward(const Vector_field_2 & vector_field_2, const Integrator_2 & integrator, Point_container_2& stl, Vertex_container_2& stl_vertices, const int & sampling_step) { - int sampling = 0; // sampling step; - int insertion = 0; // insertion order; - int insertion_step = 0; - Point_container_2 list_of_point; - Vertex_container_2 list_of_vertex; - Point_2 pPoint1; - bool bEnd = false; - FT dist; - Point_2 new_point = Point_2 (stl.back().x(),stl.back().y()); - Vertex_handle m_Vertex_handle = insert_point(new_point, stl_vertices.back()->face(), dist,true); - stl_vertices.push_back(m_Vertex_handle); - stl.push_back(new_point); - number_of_points++; - Point_2 old_point; - while (!bEnd) + int sampling = 0; // sampling step; + int insertion = 0; // insertion order; + int insertion_step = 0; + Point_container_2 list_of_point; + Vertex_container_2 list_of_vertex; + Point_2 pPoint1; + bool bEnd = false; + FT dist; + Point_2 new_point = Point_2 (stl.back().x(),stl.back().y()); + Vertex_handle m_Vertex_handle = insert_point(new_point, stl_vertices.back()->face(), dist,true); + stl_vertices.push_back(m_Vertex_handle); + stl.push_back(new_point); + number_of_points++; + Point_2 old_point; + while (!bEnd) + { + Point_2 ex_old_point = old_point; + old_point = new_point; + std::pair field_vector; + CGAL_streamlines_precondition(vector_field_2.is_in_domain(old_point)); + new_point = integrator(old_point,vector_field_2,false); + bEnd = !vector_field_2.is_in_domain(new_point); + FT dist_ = distance(ex_old_point,new_point); + bEnd = bEnd || dist_ <= vector_field_2.get_integration_step() || (new_point == old_point);/* to review */ + if(number_of_points > 30) + bEnd = bEnd || ((distance(stl.front(), stl.back())) 3000); + if (!bEnd) { - Point_2 ex_old_point = old_point; - old_point = new_point; - std::pair field_vector; - CGAL_streamlines_precondition(vector_field_2.is_in_domain(old_point)); - new_point = integrator(old_point,vector_field_2,false); - bEnd = !vector_field_2.is_in_domain(new_point); - FT dist_ = distance(ex_old_point,new_point); - bEnd = bEnd || dist_ <= vector_field_2.get_integration_step() || (new_point == old_point);/* to review */ - if(number_of_points > 30) - bEnd = bEnd || ((distance(stl.front(), stl.back())) 3000); - if (!bEnd) + if(sampling != sampling_step) + { + stl.push_back(new_point); + number_of_points ++; + sampling++; + } + else + { + if (insertion != insertion_step) { - if(sampling != sampling_step) - { - stl.push_back(new_point); - number_of_points ++; - sampling++; - } - else - { - if (insertion != insertion_step) - { - stl.push_back(new_point); - number_of_points++; - insertion++; - list_of_point.push_back(new_point); - } - else - { - stl.push_back(new_point); - number_of_points++; - list_of_point.push_back(new_point); - list_of_point.pop_back(); - m_Vertex_handle = insert_point(new_point, stl_vertices.back()->face(), dist,true); - while ((dist <= separating_distance)&&(!list_of_point.empty())) - { - m_DT.remove(m_Vertex_handle); - for (int i=0;i<=sampling_step;i++) - { - stl.pop_back(); - number_of_points--; - } - new_point = list_of_point.back(); - list_of_point.pop_back(); - m_Vertex_handle = insert_point(new_point, stl_vertices.back()->face(), dist,true); - } - // adaptive insertion order coefficient - insertion_step = (int) (((dist)-fSepStl_seed) / - std::max((FT) sampling_step,vector_field_2.get_integration_step())); - if (insertion_step < 0) insertion_step = 0; - list_of_vertex.push_back(m_Vertex_handle); - (bEnd) = (((bEnd))||(distface(), dist, false); - list_of_vertex.push_front(m_Vertex_handle); - list_of_point.pop_back(); - } - while(!list_of_vertex.empty()) - { - stl_vertices.push_back(list_of_vertex.back()); - list_of_vertex.pop_back(); - } - insertion = 0; - } - sampling = 0; - } + stl.push_back(new_point); + number_of_points++; + insertion++; + list_of_point.push_back(new_point); } - else + else { - if (!list_of_point.empty()) + stl.push_back(new_point); + number_of_points++; + list_of_point.push_back(new_point); + list_of_point.pop_back(); + m_Vertex_handle = insert_point(new_point, stl_vertices.back()->face(), dist,true); + while ((dist <= separating_distance)&&(!list_of_point.empty())) + { + m_DT.remove(m_Vertex_handle); + for (int i=0;i<=sampling_step;i++) { - new_point = list_of_point.back(); - list_of_point.pop_back(); - Vertex_handle m_Vertex_handle = insert_point(new_point, stl_vertices.back()->face(), dist, true); - while ((dist <= separating_distance)&&(!list_of_point.empty())) - { - m_DT.remove(m_Vertex_handle); - for (int i=0;i<=sampling_step;i++) - { - stl.pop_back(); - number_of_points--; - } - new_point = list_of_point.back(); - list_of_point.pop_back(); - m_Vertex_handle = insert_point(new_point, stl_vertices.back()->face(), dist, true); - } - // adaptive insertion order coefficient - insertion_step = (int) (((dist)-fSepStl_seed) / - std::max((FT) sampling_step,vector_field_2.get_integration_step())); - if (insertion_step < 0) insertion_step = 0; -// list_of_vertex.push_front(m_Vertex_handle); - (bEnd) = (((bEnd))||(distface(), dist, false); - list_of_vertex.push_back(m_Vertex_handle); - list_of_point.pop_back(); - } - while(!list_of_vertex.empty()) - { - stl_vertices.push_back(list_of_vertex.back()); - list_of_vertex.pop_back(); + stl.pop_back(); + number_of_points--; } + new_point = list_of_point.back(); + list_of_point.pop_back(); + m_Vertex_handle = insert_point(new_point, stl_vertices.back()->face(), dist,true); + } + // adaptive insertion order coefficient + insertion_step = (int) (((dist)-fSepStl_seed) / + std::max((FT) sampling_step,vector_field_2.get_integration_step())); + if (insertion_step < 0) insertion_step = 0; + list_of_vertex.push_back(m_Vertex_handle); + (bEnd) = (((bEnd))||(distface(), dist, false); + list_of_vertex.push_front(m_Vertex_handle); + list_of_point.pop_back(); + } + while(!list_of_vertex.empty()) + { + stl_vertices.push_back(list_of_vertex.back()); + list_of_vertex.pop_back(); + } + insertion = 0; } + sampling = 0; + } } + else + { + if (!list_of_point.empty()) + { + new_point = list_of_point.back(); + list_of_point.pop_back(); + Vertex_handle m_Vertex_handle = insert_point(new_point, stl_vertices.back()->face(), dist, true); + while ((dist <= separating_distance)&&(!list_of_point.empty())) + { + m_DT.remove(m_Vertex_handle); + for (int i=0;i<=sampling_step;i++) + { + stl.pop_back(); + number_of_points--; + } + new_point = list_of_point.back(); + list_of_point.pop_back(); + m_Vertex_handle = insert_point(new_point, stl_vertices.back()->face(), dist, true); + } + // adaptive insertion order coefficient + insertion_step = (int) (((dist)-fSepStl_seed) / + std::max((FT) sampling_step,vector_field_2.get_integration_step())); + if (insertion_step < 0) insertion_step = 0; + // list_of_vertex.push_front(m_Vertex_handle); + (bEnd) = (((bEnd))||(distface(), dist, false); + list_of_vertex.push_back(m_Vertex_handle); + list_of_point.pop_back(); + } + while(!list_of_vertex.empty()) + { + stl_vertices.push_back(list_of_vertex.back()); + list_of_vertex.pop_back(); + } + } + } } template inline void Stream_lines_2:: -insert_streamline(const Vector_field_2 & vector_field_2, Point_container_2 stl, Vertex_container_2 stl_vertices){ - stl_container.push_back(stl); - Vertex_handle m_Vertex_handle = NULL; - int i = 1; - unsigned int size_ = (int) (stl_vertices.size()); - ir = il = 0; - while (!stl_vertices.empty()){ - pq_elements(vector_field_2, stl_vertices, i, m_Vertex_handle, size_); - m_Vertex_handle = stl_vertices.front(); - stl_vertices.pop_front(); - i++;}} +insert_streamline(const Vector_field_2 & vector_field_2, + Point_container_2 stl, Vertex_container_2 + stl_vertices) +{ + stl_container.push_back(stl); + Vertex_handle m_Vertex_handle = NULL; + int i = 1; + unsigned int size_ = (int) (stl_vertices.size()); + ir = il = 0; + while (!stl_vertices.empty()) + { + pq_elements(vector_field_2, stl_vertices, i, m_Vertex_handle, size_); + m_Vertex_handle = stl_vertices.front(); + stl_vertices.pop_front(); + i++; + } +} template void Stream_lines_2:: pq_elements(const Vector_field_2 & vector_field_2, Vertex_container_2 stl_vertices, int i, - const Vertex_handle & m_Vertex_handle, int size_){ - if ((i!=0) && (i!=(size_))){// && (std::div(i,10).rem!=0)){ - Vertex_handle pVertex_handle = stl_vertices.front(); - Face_handle pFace_handle; - int iIndex; - if (m_DT.is_edge(pVertex_handle,m_Vertex_handle,pFace_handle,iIndex)){ - Point_2 p0 = pVertex_handle->point(); - Point_2 c = m_DT.circumcenter(pFace_handle); - FT fDist = distance(p0,c); - bool b = vector_field_2.is_in_domain(c) && (fDist >= fSepStl_seed); - if (b){ - Circle pCircle(c,fDist); - Pq_element m_Pq_element = Pq_element( - pFace_handle->vertex(0), - pFace_handle->vertex(1), - pFace_handle->vertex(2), - pCircle); - if (ir == 0){ - Pq_previous_r = m_Pq_element; - Pq_element_max_r = m_Pq_element; - ir++;} - else if (ir == 1){ - Pq_current_r = m_Pq_element;ir++;} - else{ - Pq_next_r = m_Pq_element; - if (Pq_element_max_r.fourth.second <= Pq_next_r.fourth.second) - Pq_element_max_r = Pq_next_r; - if ((Pq_current_r.fourth.second>=Pq_previous_r.fourth.second) - &&(Pq_current_r.fourth.second>=Pq_next_r.fourth.second)){ - pq.push(Pq_current_r);} - Pq_previous_r = Pq_current_r; - Pq_current_r = Pq_next_r; - ir++;}} - p0 = pFace_handle->neighbor(iIndex)->vertex(0)->point(); - c = m_DT.circumcenter(pFace_handle->neighbor(iIndex)); - fDist = distance(p0,c); - b = vector_field_2.is_in_domain(c) && (fDist >= fSepStl_seed); - if (b){ - Circle pCircle(c,fDist); - Pq_element m_Pq_element = Pq_element( - pFace_handle->neighbor(iIndex)->vertex(0), - pFace_handle->neighbor(iIndex)->vertex(1), - pFace_handle->neighbor(iIndex)->vertex(2), - pCircle); - if (il == 0){ - Pq_previous_l = m_Pq_element; - Pq_element_max_l = m_Pq_element; - il++;} - else if (il == 1){ - Pq_current_l = m_Pq_element;il++;} - else{ - Pq_next_l = m_Pq_element; - if (Pq_element_max_l.fourth.second <= Pq_next_l.fourth.second) - Pq_element_max_l = Pq_next_l; - if ((Pq_current_l.fourth.second>=Pq_previous_l.fourth.second) - &&(Pq_current_l.fourth.second>=Pq_next_l.fourth.second)){ - pq.push(Pq_current_l);} - Pq_previous_l = Pq_current_l; - Pq_current_l = Pq_next_l; - il++;}}} - if ((ir+il == (int) size_-2)&&(size_>2)){ - pq.push(Pq_element_max_l); - pq.push(Pq_element_max_r);}} - else{ - Vertex_handle pVertex_handle = stl_vertices.front(); - Face_circulator pFace_handle = m_DT.incident_faces(pVertex_handle); - Face_circulator pEnd = pFace_handle; - CGAL_For_all(pFace_handle,pEnd){ - Point_2 p0 = pFace_handle->vertex(0)->point(); - Point_2 c = m_DT.circumcenter(pFace_handle); - bool b = vector_field_2.is_in_domain(c); - if (b){ - FT fDist = distance(p0,c); - if (fDist >= fSepStl_seed){ - Circle pCircle(c,fDist); - Pq_element m_Pq_element = Pq_element( - pFace_handle->vertex(0), - pFace_handle->vertex(1), - pFace_handle->vertex(2), - pCircle); - pq.push(m_Pq_element);}}}}} + const Vertex_handle & m_Vertex_handle, int size_) +{ + if ((i!=0) && (i!=(size_))){// && (std::div(i,10).rem!=0)){ + Vertex_handle pVertex_handle = stl_vertices.front(); + Face_handle pFace_handle; + int iIndex; + if + (m_DT.is_edge(pVertex_handle,m_Vertex_handle,pFace_handle,iIndex)) + { + Point_2 p0 = pVertex_handle->point(); + Point_2 c = m_DT.circumcenter(pFace_handle); + FT fDist = distance(p0,c); + bool b = vector_field_2.is_in_domain(c) && (fDist >= fSepStl_seed); + if (b) + { + Circle pCircle(c,fDist); + Pq_element m_Pq_element = Pq_element( + pFace_handle->vertex(0), + pFace_handle->vertex(1), + pFace_handle->vertex(2), + pCircle); + if (ir == 0) + { + Pq_previous_r = m_Pq_element; + Pq_element_max_r = m_Pq_element; + ir++; + } + else if (ir == 1) + { + Pq_current_r = + m_Pq_element;ir++; + } + else + { + Pq_next_r = m_Pq_element; + if (Pq_element_max_r.fourth.second <= Pq_next_r.fourth.second) + Pq_element_max_r = Pq_next_r; + if ((Pq_current_r.fourth.second>=Pq_previous_r.fourth.second) + &&(Pq_current_r.fourth.second>=Pq_next_r.fourth.second)) + { + pq.push(Pq_current_r); + } + Pq_previous_r = Pq_current_r; + Pq_current_r = Pq_next_r; + ir++; + } + } + p0 = pFace_handle->neighbor(iIndex)->vertex(0)->point(); + c = m_DT.circumcenter(pFace_handle->neighbor(iIndex)); + fDist = distance(p0,c); + b = vector_field_2.is_in_domain(c) && (fDist >= fSepStl_seed); + if (b) + { + Circle pCircle(c,fDist); + Pq_element m_Pq_element = Pq_element( + pFace_handle->neighbor(iIndex)->vertex(0), + pFace_handle->neighbor(iIndex)->vertex(1), + pFace_handle->neighbor(iIndex)->vertex(2), + pCircle); + if (il == 0) + { + Pq_previous_l = m_Pq_element; + Pq_element_max_l = m_Pq_element; + il++; + } + else if (il == 1) + { + Pq_current_l = + m_Pq_element;il++; + } + else + { + Pq_next_l = m_Pq_element; + if (Pq_element_max_l.fourth.second <= Pq_next_l.fourth.second) + Pq_element_max_l = Pq_next_l; + if ((Pq_current_l.fourth.second>=Pq_previous_l.fourth.second) + &&(Pq_current_l.fourth.second>=Pq_next_l.fourth.second)) + { + pq.push(Pq_current_l); + } + Pq_previous_l = Pq_current_l; + Pq_current_l = Pq_next_l; + il++; + } + } + } + if ((ir+il == (int) size_-2)&&(size_>2)) + { + pq.push(Pq_element_max_l); + pq.push(Pq_element_max_r); + } + } + else{ + Vertex_handle pVertex_handle = stl_vertices.front(); + Face_circulator pFace_handle = m_DT.incident_faces(pVertex_handle); + Face_circulator pEnd = pFace_handle; + CGAL_For_all(pFace_handle,pEnd) + { + Point_2 p0 = pFace_handle->vertex(0)->point(); + Point_2 c = m_DT.circumcenter(pFace_handle); + bool b = vector_field_2.is_in_domain(c); + if (b){ + FT fDist = distance(p0,c); + if (fDist >= fSepStl_seed) + { + Circle pCircle(c,fDist); + Pq_element m_Pq_element = Pq_element( + pFace_handle->vertex(0), + pFace_handle->vertex(1), + pFace_handle->vertex(2), + pCircle); + pq.push(m_Pq_element); + } + } + } + } +} // get the next seed point template inline bool -Stream_lines_2::get_next_seed_point(FT & distance, Point_2 & seed_point){ - Vertex_handle v0, v1, v2; - Face_handle fr; - bool b0,b; - Pq_element m_Pq_element; - do{ - m_Pq_element = pq.top(); - v0 = m_Pq_element.first; - v1 = m_Pq_element.second; - v2 = m_Pq_element.third; - distance = m_Pq_element.fourth.second; - pq.pop(); - b0 = m_DT.is_face(v0,v1,v2,fr); - if (b0){ - seed_point = m_Pq_element.fourth.first;} - b = (!pq.empty()); - }while ((b)&&(!b0)); - Biggest_circle = m_Pq_element; - return b;} +Stream_lines_2::get_next_seed_point(FT & + distance, Point_2 & seed_point) +{ + Vertex_handle v0, v1, v2; + Face_handle fr; + bool b0,b; + Pq_element m_Pq_element; + do{ + m_Pq_element = pq.top(); + v0 = m_Pq_element.first; + v1 = m_Pq_element.second; + v2 = m_Pq_element.third; + distance = m_Pq_element.fourth.second; + pq.pop(); + b0 = m_DT.is_face(v0,v1,v2,fr); + if (b0){ + seed_point = m_Pq_element.fourth.first;} + b = (!pq.empty()); + }while ((b)&&(!b0)); + Biggest_circle = m_Pq_element; + return b; +} template typename Stream_lines_2::Stream_line_iterator_2 -Stream_lines_2::begin(){ - return begin_iterator;} +Stream_lines_2::begin() +{ + return begin_iterator; +} template typename Stream_lines_2::Stream_line_iterator_2 -Stream_lines_2::end(){ - return end_iterator;} +Stream_lines_2::end() +{ + return end_iterator; +} template inline -void Stream_lines_2::make_iterator(){ - for(typename Stream_line_container_2::iterator begin=stl_container.begin(); begin!=stl_container.end();begin++){ - std::pair - iterator_pair((*begin).begin(), (*begin).end()); - iterator_container.push_front(iterator_pair);} - begin_iterator = iterator_container.begin(); - end_iterator = iterator_container.end();} +void Stream_lines_2::make_iterator() +{ + for(typename Stream_line_container_2::iterator + begin=stl_container.begin(); + begin!=stl_container.end();begin++) + { + std::pair + iterator_pair((*begin).begin(), (*begin).end()); + iterator_container.push_front(iterator_pair); + } + begin_iterator = iterator_container.begin(); + end_iterator = iterator_container.end(); +} // output an eps file @@ -714,79 +783,88 @@ template void Stream_lines_2::print_stream_lines_eps(std::ofstream & fw) { - typename Stream_line_container_2::iterator begin_iterator; - Stream_line_container_2 stl_container_temp = stl_container; - Point_container_2 stl; - fw << "%!PS-Adobe-2.0 EPSF-2.0\n"; - fw << "%%BoundingBox: 0 0" << " " << max_x - min_x << " " << max_y - min_y << "\n"; - fw << "gsave\n"; - fw << "/L {moveto lineto stroke} bind def\n"; - fw << 0.5 << " setlinewidth\n"; - fw << 0.0 << " " << 0.0 << " " << 0.0 << " setrgbcolor\n"; + typename Stream_line_container_2::iterator begin_iterator; + Stream_line_container_2 stl_container_temp = stl_container; + Point_container_2 stl; + fw << "%!PS-Adobe-2.0 EPSF-2.0\n"; + fw << "%%BoundingBox: 0 0" << " " << max_x - min_x << " " << max_y - min_y << "\n"; + fw << "gsave\n"; + fw << "/L {moveto lineto stroke} bind def\n"; + fw << 0.5 << " setlinewidth\n"; + fw << 0.0 << " " << 0.0 << " " << 0.0 << " setrgbcolor\n"; - for(begin_iterator=stl_container_temp.begin();begin_iterator!=stl_container_temp.end();begin_iterator++){ - typename Point_container_2::iterator begin_point_iterator = (*begin_iterator).begin(); - typename Point_container_2::iterator end_point_iterator = (*begin_iterator).end(); + for(begin_iterator=stl_container_temp.begin();begin_iterator!=stl_container_temp.end();begin_iterator++) + { + typename Point_container_2::iterator begin_point_iterator = (*begin_iterator).begin(); + typename Point_container_2::iterator end_point_iterator = (*begin_iterator).end(); - FT i_prec = (*begin_point_iterator).x() - min_x; - FT j_prec = (*begin_point_iterator).y() - min_y; - begin_point_iterator++; + FT i_prec = (*begin_point_iterator).x() - min_x; + FT j_prec = (*begin_point_iterator).y() - min_y; + begin_point_iterator++; - FT i , j; - for(;begin_point_iterator!=end_point_iterator;begin_point_iterator++){ - Point_2 p = *begin_point_iterator; - i = p.x() - min_x; - j = p.y() - min_y; - fw << i_prec << " " << j_prec << " " << i << " " << j << " L\n"; - i_prec = i; - j_prec = j;}} - fw << "grestore\n"; - fw << "showpage\n"; - fw.close();} + FT i , j; + for(;begin_point_iterator!=end_point_iterator;begin_point_iterator++) + { + Point_2 p = *begin_point_iterator; + i = p.x() - min_x; + j = p.y() - min_y; + fw << i_prec << " " << j_prec << " " << i << " " << j << " L\n"; + i_prec = i; + j_prec = j; + } + } + fw << "grestore\n"; + fw << "showpage\n"; + fw.close(); +} // output an stl file template void Stream_lines_2::print_stream_lines(std::ofstream & fw) { - typename Stream_line_container_2::iterator begin_iterator; - Stream_line_container_2 stl_container_temp = stl_container; - Point_container_2 stl; - fw << max_x - min_x << " " << max_y - min_y << "\n"; - fw << stl_container.size() << "\n"; - for(begin_iterator=stl_container_temp.begin();begin_iterator!=stl_container_temp.end();begin_iterator++){ - fw << (*begin_iterator).size() << "\n"; - typename Point_container_2::iterator begin_point_iterator = (*begin_iterator).begin(); - typename Point_container_2::iterator end_point_iterator = (*begin_iterator).end(); - FT i , j; - for(;begin_point_iterator!=end_point_iterator;begin_point_iterator++){ - Point_2 p = *begin_point_iterator; - i = p.x() - min_x; - j = p.y() - min_y; - fw << i << " " << j << "\n";}} - fw.close();} + typename Stream_line_container_2::iterator begin_iterator; + Stream_line_container_2 stl_container_temp = stl_container; + Point_container_2 stl; + fw << max_x - min_x << " " << max_y - min_y << "\n"; + fw << stl_container.size() << "\n"; + for(begin_iterator=stl_container_temp.begin();begin_iterator!=stl_container_temp.end();begin_iterator++) + { + fw << (*begin_iterator).size() << "\n"; + typename Point_container_2::iterator begin_point_iterator = (*begin_iterator).begin(); + typename Point_container_2::iterator end_point_iterator = (*begin_iterator).end(); + FT i , j; + for(;begin_point_iterator!=end_point_iterator;begin_point_iterator++){ + Point_2 p = *begin_point_iterator; + i = p.x() - min_x; + j = p.y() - min_y; + fw << i << " " << j << "\n"; + } + } + fw.close(); +} template std::list::Point_2> Stream_lines_2::get_pq() { - std::list _list; - std::priority_queue, C> pq_temp; - pq_temp = pq; - while (!pq_temp.empty()) - { - Pq_element m_Pq_element = pq_temp.top(); - Vertex_handle v0 = m_Pq_element.first; - Vertex_handle v1 = m_Pq_element.second; - Vertex_handle v2 = m_Pq_element.third; - pq_temp.pop(); - Face_handle fr; - bool b0 = m_DT.is_face(v0,v1,v2,fr); - Point_2 sdPoint = m_Pq_element.fourth.first; - if (b0) - _list.push_front(sdPoint); - } - return _list; + std::list _list; + std::priority_queue, C> pq_temp; + pq_temp = pq; + while (!pq_temp.empty()) + { + Pq_element m_Pq_element = pq_temp.top(); + Vertex_handle v0 = m_Pq_element.first; + Vertex_handle v1 = m_Pq_element.second; + Vertex_handle v2 = m_Pq_element.third; + pq_temp.pop(); + Face_handle fr; + bool b0 = m_DT.is_face(v0,v1,v2,fr); + Point_2 sdPoint = m_Pq_element.fourth.first; + if (b0) + _list.push_front(sdPoint); + } + return _list; } CGAL_END_NAMESPACE diff --git a/Packages/Stream_lines_2/include/CGAL/Triangular_field_2.h b/Packages/Stream_lines_2/include/CGAL/Triangular_field_2.h index 4bc22c2b1d5..96ee8870335 100644 --- a/Packages/Stream_lines_2/include/CGAL/Triangular_field_2.h +++ b/Packages/Stream_lines_2/include/CGAL/Triangular_field_2.h @@ -20,19 +20,15 @@ #ifndef CGAL_TRIANGULAR_FIELD_2_H_ #define CGAL_TRIANGULAR_FIELD_2_H_ -#include -#include #include #include #include -// #include #include #include #include #include -#include #include #include @@ -41,143 +37,190 @@ CGAL_BEGIN_NAMESPACE template -class Triangular_field_2{ +class Triangular_field_2 +{ public: - typedef Triangular_field_2 Vector_field_2; - typedef StreamLinesTraits_2 Geom_traits; - typedef typename StreamLinesTraits_2::FT FT; - typedef typename StreamLinesTraits_2::Point_2 Point_2; - typedef typename StreamLinesTraits_2::Vector_2 Vector_2; + typedef Triangular_field_2 Vector_field_2; + typedef StreamLinesTraits_2 Geom_traits; + typedef typename StreamLinesTraits_2::FT FT; + typedef typename StreamLinesTraits_2::Point_2 Point_2; + typedef typename StreamLinesTraits_2::Vector_2 Vector_2; protected: - typedef CGAL::Triangulation_vertex_base_2 Vb; - typedef CGAL::Triangulation_face_base_with_info_2 Fb; - typedef CGAL::Triangulation_data_structure_2 TDS; - typedef CGAL::Delaunay_triangulation_2 D_Ttr; - typedef typename D_Ttr::Vertex_handle Vertex_handle; - typedef typename TDS::Face_handle Face_handle; + typedef CGAL::Triangulation_vertex_base_2 Vb; + typedef CGAL::Triangulation_face_base_with_info_2 Fb; + typedef CGAL::Triangulation_data_structure_2 TDS; + typedef CGAL::Delaunay_triangulation_2 D_Ttr; + typedef typename D_Ttr::Vertex_handle Vertex_handle; + typedef typename TDS::Face_handle Face_handle; public: - typedef typename D_Ttr::Vertex_iterator Vertex_iterator; - D_Ttr m_D_Ttr; + typedef typename D_Ttr::Vertex_iterator Vertex_iterator; + D_Ttr m_D_Ttr; protected: - Vector_2 get_vector_field(const Point_2 & p) const; - FT get_density_field(const Point_2 & p) const; - template - void fill(PointInputIterator pBegin, PointInputIterator pEnd, VectorInputIterator vBegin){ - std::cout << "reading file...\n"; - while(pBegin != pEnd){ - Point_2 p; - Vector_2 v; - p = (*pBegin); - v = (*vBegin); - Vertex_handle m_Vertex_handle = m_D_Ttr.insert(p); - field_map[p] = v; - if (m_D_Ttr.number_of_vertices() == 1){ - maxx = minx = p.x(); - maxy = miny = p.y();} - if(p.x()maxx) - maxx = p.x(); - if(p.y()>maxy) - maxy = p.y(); - pBegin++; - vBegin++;} - std::cout << "number of samples " << m_D_Ttr.number_of_vertices() << "\n";} -public: - template - Triangular_field_2(PointInputIterator pBegin, PointInputIterator pEnd, VectorInputIterator vBegin){ - fill(pBegin, pEnd, vBegin);} - inline typename Geom_traits::Iso_rectangle_2 bbox() const; + Vector_2 get_vector_field(const Point_2 & p) const; - std::pair get_field(const Point_2 & p) const - { - assert(is_in_domain(p)); - Vector_2 v = get_vector_field(p); - FT density = get_density_field(p); - std::pair field_value(v,density); - return field_value; + FT get_density_field(const Point_2 & p) const; + + template + void fill(PointInputIterator pBegin, PointInputIterator pEnd, VectorInputIterator vBegin) + { + std::cout << "reading file...\n"; + while(pBegin != pEnd) + { + Point_2 p; + Vector_2 v; + p = (*pBegin); + v = (*vBegin); + Vertex_handle m_Vertex_handle = m_D_Ttr.insert(p); + field_map[p] = v; + if (m_D_Ttr.number_of_vertices() == 1) + { + maxx = minx = p.x(); + maxy = miny = p.y(); + } + if(p.x()maxx) + { + maxx = p.x(); + } + if(p.y()>maxy) + { + maxy = p.y(); + } + ++pBegin; + ++vBegin; } + std::cout << "number of samples " << m_D_Ttr.number_of_vertices() << "\n"; + } +public: + template + Triangular_field_2(PointInputIterator pBegin, PointInputIterator + pEnd, VectorInputIterator vBegin) + { + fill(pBegin, pEnd, vBegin); + } + + inline typename Geom_traits::Iso_rectangle_2 bbox() const; + + std::pair get_field(const Point_2 & p) const + { + assert(is_in_domain(p)); + Vector_2 v = get_vector_field(p); + FT density = get_density_field(p); + return std::make_pair(v, density); + } - bool is_in_domain(const Point_2 & p) const; - FT get_integration_step(const Point_2 &) const; - FT get_integration_step() const; + bool is_in_domain(const Point_2 & p) const; + + FT get_integration_step(const Point_2 &) const; + + FT get_integration_step() const; protected: - FT minx; - FT miny; - FT maxx; - FT maxy; + FT minx; + FT miny; + FT maxx; + FT maxy; protected: - mutable std::map field_map; - FT distance(const Point_2 & p, const Point_2 & q){ - return sqrt(((p.x() - q.x()) * (p.x() - q.x())) + ((p.y() - q.y()) * (p.y() - q.y())));} + mutable std::map field_map; + + FT distance(const Point_2 & p, const Point_2 & q) + { + return sqrt( CGAL::squared_distance(p, q) ); + } }; template inline typename Triangular_field_2::Geom_traits::Iso_rectangle_2 -Triangular_field_2::bbox() const{ - return typename Geom_traits::Iso_rectangle_2(minx, miny, maxx, maxy);} +Triangular_field_2::bbox() const +{ + return typename Geom_traits::Iso_rectangle_2(minx, miny, maxx, + maxy); +} template bool -Triangular_field_2::is_in_domain(const Point_2 & p) const{ - Face_handle f = m_D_Ttr.locate(p); - return !m_D_Ttr.is_infinite(f);} +Triangular_field_2::is_in_domain(const Point_2 & + p) const +{ + Face_handle f = m_D_Ttr.locate(p); + return !m_D_Ttr.is_infinite(f); +} template typename Triangular_field_2::Vector_2 -Triangular_field_2::get_vector_field(const Point_2 & p) const{ - Face_handle m_Face_handle = m_D_Ttr.locate(p); - assert(is_in_domain(p)); - Point_2 p0 = m_Face_handle->vertex(0)->point(); - Point_2 p1 = m_Face_handle->vertex(1)->point(); - Point_2 p2 = m_Face_handle->vertex(2)->point(); - Vertex_handle v0 = m_Face_handle->vertex(0); - Vertex_handle v1 = m_Face_handle->vertex(1); - Vertex_handle v2 = m_Face_handle->vertex(2); - FT s0,s1,s2,s; - std::vector vec; - vec.push_back(p0); vec.push_back(p1); vec.push_back(p2); - s = polygon_area_2(vec.begin(), vec.end(), m_D_Ttr.geom_traits()); - vec.clear(); - vec.push_back(p); vec.push_back(p1); vec.push_back(p2); - s0 = polygon_area_2(vec.begin(), vec.end(), m_D_Ttr.geom_traits()); - vec.clear(); - vec.push_back(p0); vec.push_back(p); vec.push_back(p2); - s1 = polygon_area_2(vec.begin(), vec.end(), m_D_Ttr.geom_traits()); - vec.clear(); - vec.push_back(p0); vec.push_back(p1); vec.push_back(p); - s2 = polygon_area_2(vec.begin(), vec.end(), m_D_Ttr.geom_traits()); - vec.clear(); - s0 = s0 / s; s1 = s1 / s; s2 = s2 / s; - Vector_2 v_0 = field_map[p0]; - Vector_2 v_1 = field_map[p1]; - Vector_2 v_2 = field_map[p2]; - FT x = ((v_0.x()*s0)+(v_1.x()*s1)+(v_2.x()*s2)); - FT y = ((v_0.y()*s0)+(v_1.y()*s1)+(v_2.y()*s2)); - FT normal = sqrt((x)*(x) + (y)*(y)); - if (normal != 0){ - x = x / normal; - y = y / normal;} - Vector_2 v = Vector_2(x, y); -return v;} +Triangular_field_2::get_vector_field(const + Point_2 & p) + const +{ + Face_handle m_Face_handle = m_D_Ttr.locate(p); + assert(is_in_domain(p)); + Vertex_handle v0 = m_Face_handle->vertex(0); + Vertex_handle v1 = m_Face_handle->vertex(1); + Vertex_handle v2 = m_Face_handle->vertex(2); + const Point_2 & p0 = v0->point(); + const Point_2 & p1 = v1->point(); + const Point_2 & p2 = v2->point(); + FT s0,s1,s2,s; + std::vector vec; + vec.push_back(p0); vec.push_back(p1); vec.push_back(p2); + s = polygon_area_2(vec.begin(), vec.end(), m_D_Ttr.geom_traits()); + vec.clear(); + vec.push_back(p); vec.push_back(p1); vec.push_back(p2); + s0 = polygon_area_2(vec.begin(), vec.end(), m_D_Ttr.geom_traits()); + vec.clear(); + vec.push_back(p0); vec.push_back(p); vec.push_back(p2); + s1 = polygon_area_2(vec.begin(), vec.end(), m_D_Ttr.geom_traits()); + vec.clear(); + vec.push_back(p0); vec.push_back(p1); vec.push_back(p); + s2 = polygon_area_2(vec.begin(), vec.end(), m_D_Ttr.geom_traits()); + vec.clear(); + s0 = s0 / s; s1 = s1 / s; s2 = s2 / s; + Vector_2 v_0 = field_map[p0]; + Vector_2 v_1 = field_map[p1]; + Vector_2 v_2 = field_map[p2]; + FT x = ((v_0.x()*s0)+(v_1.x()*s1)+(v_2.x()*s2)); + FT y = ((v_0.y()*s0)+(v_1.y()*s1)+(v_2.y()*s2)); + FT normal = sqrt((x)*(x) + (y)*(y)); + if (normal != 0) + { + x = x / normal; + y = y / normal; + } + Vector_2 v = Vector_2(x, y); + return v; +} template typename Triangular_field_2::FT -Triangular_field_2::get_density_field(const Point_2 & p) const{ - return p.x();} +Triangular_field_2::get_density_field(const + Point_2 & + p) const +{ + return p.x(); +} template typename Triangular_field_2::FT -Triangular_field_2::get_integration_step(const Point_2 &) const{ - return 1.0;} +Triangular_field_2::get_integration_step(const + Point_2 + &) const +{ + return 1.0; +} template typename Triangular_field_2::FT -Triangular_field_2::get_integration_step() const{ - return 1.0;} +Triangular_field_2::get_integration_step() const +{ + return 1.0; +} CGAL_END_NAMESPACE