From b84bddbcbdaccaa8df65c1f0f4c9db958d5224cb Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Fri, 9 May 2003 23:02:21 +0000 Subject: [PATCH] Commit of changes for 0.9.3: split of Mesh_2.h into Mesh_2.h and Conform_2.h --- Packages/Mesh_2/changes.txt | 4 + .../demo/Mesh_2/Qt_layer_show_clusters.h | 159 ++ .../Mesh_2/demo/Mesh_2/Qt_layer_show_points.h | 114 +- Packages/Mesh_2/demo/Mesh_2/makefile | 12 +- Packages/Mesh_2/demo/Mesh_2/mesh_demo.C | 256 ++-- Packages/Mesh_2/examples/Mesh_2/Makefile | 18 +- Packages/Mesh_2/examples/Mesh_2/conform.C | 129 ++ Packages/Mesh_2/examples/Mesh_2/mesh.C | 13 +- Packages/Mesh_2/include/CGAL/Conform_2.h | 1140 ++++++++++++++ .../Constrained_Delaunay_triangulation_2.h | 664 -------- ...ltred_circulator.h => Filter_circulator.h} | 0 ...ltred_container.h => Filtered_container.h} | 6 +- Packages/Mesh_2/include/CGAL/Mesh_2.h | 1348 ++++------------- Packages/Mesh_2/include/CGAL/Read_write.h | 174 +++ .../test/Mesh_2/test_filtred_container.C | 4 +- 15 files changed, 2108 insertions(+), 1933 deletions(-) create mode 100644 Packages/Mesh_2/demo/Mesh_2/Qt_layer_show_clusters.h create mode 100644 Packages/Mesh_2/examples/Mesh_2/conform.C create mode 100644 Packages/Mesh_2/include/CGAL/Conform_2.h delete mode 100644 Packages/Mesh_2/include/CGAL/Constrained_Delaunay_triangulation_2.h rename Packages/Mesh_2/include/CGAL/{Filtred_circulator.h => Filter_circulator.h} (100%) rename Packages/Mesh_2/include/CGAL/{Filtred_container.h => Filtered_container.h} (87%) create mode 100644 Packages/Mesh_2/include/CGAL/Read_write.h diff --git a/Packages/Mesh_2/changes.txt b/Packages/Mesh_2/changes.txt index 4aef2787f2c..5d4893c7818 100644 --- a/Packages/Mesh_2/changes.txt +++ b/Packages/Mesh_2/changes.txt @@ -1,3 +1,7 @@ +version 0.9.3 (2003/05/10) +-------------------------- +o split of the code into Mesh_2 and Conform_2 + version 0.9.2 (2003/03/12) -------------------------- o added a bench sub-directory diff --git a/Packages/Mesh_2/demo/Mesh_2/Qt_layer_show_clusters.h b/Packages/Mesh_2/demo/Mesh_2/Qt_layer_show_clusters.h new file mode 100644 index 00000000000..65e64efe647 --- /dev/null +++ b/Packages/Mesh_2/demo/Mesh_2/Qt_layer_show_clusters.h @@ -0,0 +1,159 @@ +// ============================================================================ +// +// Copyright (c) 1997-2000 The CGAL Consortium +// +// This software and related documentation is part of an INTERNAL release +// of the Computational Geometry Algorithms Library (CGAL). It is not +// intended for general use. +// +// ---------------------------------------------------------------------------- +// +// file : +// package : +// author(s) : Laurent Rineau +// release : +// release_date : +// +// +// +// ============================================================================ + +#include +#include + +class Show_clusters_aux : public CGAL::Qt_widget_layer +{ + Q_OBJECT +private: + virtual voir reinit_clusters() {} + +public slots: + void reinitClusters() + { + reinit_clusters(); + } +}; + +template +class Show_clusters : public Show_clusters_aux +{ +public: + typedef typename Conform::Point Point; + typedef typename Conform::Triangulation DT; + typedef typename Conform::Segment Segment; + typedef typename Conform::Face_handle Face_handle; + typedef typename Conform::Vertex_handle Vertex_handle; + typedef typename Conform::Geom_traits::FT FT; + typedef typename Conform::Cluster_vertices_iterator CVIt; + typedef typename Conform::Vertices_in_cluster_iterator ViCIt + typedef std::list List_of_points; + typedef typename List_of_points::const_iterator Point_iterator; + + Show_clusters(Conform &conform, + Color c = CGAL::GREEN, + int pointsize = 3, + PointStyle pointstyle = CGAL::DISC, + Color lc = CGAL::RED, + int linewidth = 2) + : c(conform), first_time(true), dt(), _color(c), + size(pointsize), style(pointstyle), _line_color(lc), + width(linewidth) {} + + void reinit_clusters() + { + dt.clear(); + + for(CVIt it = c.clusters_vertices_begin(); + it != c.clusters_vertices_end(); + ++it) + dt.push_back( (*it)->point() ); + } + + void draw(){first_time = true;} + + void mouseMoveEvent(QMouseEvent *e) + { + if (dt.dimension()<1) return; + + FT x, y; + widget->x_real(e->x(), x); + widget->y_real(e->y(), y); + Point p(x, y); + + RasterOp old = widget->rasterOp(); //save the initial raster mode + widget->setRasterOp(XorROP); + widget->lock(); + + Vertex_handle v = dt.nearest_vertex(p); + *widget << _color << CGAL::PointSize (pointsize) + << CGAL::PointStyle (pointstyle); + + if(!first_time) + *widget << oldPoint; + *widget << v->point(); + + *widget << _line_color << CGAL::LineWidth(width); + if(!first_time) + for(Point_iterator pIt = oldPoints.begin(); + pIt != oldPoints.end(); + ++pIt) + *widget << Segment(oldPoint, *pIt); + oldPoints.clear(); + + typename Conform::Locate_type lt; + int i; + Face_handle fh = c.locate(v->point(), lt, t); + CGAL_assertion( lt == this->VERTEX ); + + Vertex_handle v2 = fh->vertex(i); + + int n = c.number_of_clusters_at_vertex(v2); + + for(int j = 0; j < n; ++j) + { + std::pair seq = c.vertices_in_cluster_sequence(v2, j); + for(ViCIt it = seq.first; + it != seq.second; + ++it) + { + oldPoints.push_back((*it)->point()); + *widget << Segment(*v2, (*it)->point()); + } + } + + widget->unlock(); + widget->setRasterOp(old); + oldPoint = v->point(); + first_time = false; + } + + void leaveEvent(QEvent *) + { + widget->lock(); + RasterOp old = widget->rasterOp(); //save the initial raster mode + widget->setRasterOp(XorROP); + + *widget << _color << CGAL::PointSize (pointsize) + << CGAL::PointStyle (pointstyle); + if(!first_time) *widget << oldPoint; + + widget->unlock(); + widget->setRasterOp(old); + first_time = true; + } + +private: + Conform& c; + DT dt; + Point oldPoint, newPoint; + List_of_points oldPoints; + bool first_time; + Color _color; + int size; + PointStyle style; + Color _line_color; + int width; +}; + +// moc_source_file: Show_clusters.h +#include "Show_clusters.moc" diff --git a/Packages/Mesh_2/demo/Mesh_2/Qt_layer_show_points.h b/Packages/Mesh_2/demo/Mesh_2/Qt_layer_show_points.h index afbda2aa4eb..d87251e1f84 100644 --- a/Packages/Mesh_2/demo/Mesh_2/Qt_layer_show_points.h +++ b/Packages/Mesh_2/demo/Mesh_2/Qt_layer_show_points.h @@ -1,61 +1,99 @@ -// ============================================================================ -// -// Copyright (c) 1997-2000 The CGAL Consortium -// -// This software and related documentation is part of an INTERNAL release -// of the Computational Geometry Algorithms Library (CGAL). It is not -// intended for general use. -// -// ---------------------------------------------------------------------------- -// -// file : include/CGAL/IO/Qt_layer_show_points.h -// package : Qt_widget -// author(s) : Radu Ursu -// release : -// release_date : -// -// coordinator : Laurent Rineau -// -// ============================================================================ - #ifndef CGAL_QT_LAYER_SHOW_POINTS_H #define CGAL_QT_LAYER_SHOW_POINTS_H #include +#include +#include namespace CGAL { -template +template > class Qt_layer_show_points : public Qt_widget_layer { public: -// typedef typename T::Point Point; -// typedef typename T::Segment Segment; -// typedef typename T::Vertex Vertex; - typedef typename T::Vertex_iterator Vertex_iterator; + typedef It (C::* iterator_function)() const; - Qt_layer_show_points(T *&t, Color c=CGAL::GREEN, int pointsize=3, + Qt_layer_show_points(C *&container, + iterator_function begin, + iterator_function end, + const Transform& t = Transform(), + Color c = CGAL::GREEN, + int pointsize = 3, PointStyle pointstyle = CGAL::DISC) - : tr(t), color(c), size(pointsize), style(pointstyle) {}; + : cont(container), _begin(begin), _end(end), _color(c), + size(pointsize), style(pointstyle), trans(t) {}; + + Qt_layer_show_points(C *&container, + iterator_function begin, + iterator_function end, + Color c = CGAL::GREEN, + int pointsize = 3, + PointStyle pointstyle = CGAL::DISC) + : cont(container), _begin(begin), _end(end), _color(c), + size(pointsize), style(pointstyle), trans(Transform()) {}; void draw() - { - Vertex_iterator it = tr->vertices_begin(), - beyond = tr->vertices_end(); - *widget << color << CGAL::PointSize (size) - << CGAL::PointStyle (style); - while(it != beyond) { - *widget << (*it).point(); - ++it; - } + { + widget->lock(); + QColor old_color = widget->color(); + int old_point_size = widget->pointSize(); + + *widget << _color << CGAL::PointSize (size) + << CGAL::PointStyle (style); + + for(It it = (cont->*_begin)(); + it!=(cont->*_end)(); + ++it) + *widget << trans(*it); + + widget->setPointSize(old_point_size); + widget->setColor(old_color); + widget->unlock(); }; + + QColor color() const { return _color; }; + + int pointSize() const { return size; }; + + void setColor(const QColor c) + { + _color = c; + widget->redraw(); + }; + + void setPointSize(const int s) + { + size = s; + widget->redraw(); + }; + private: - T *&tr; - Color color; + C *&cont; + iterator_function _begin; + iterator_function _end; + Color _color; int size; PointStyle style; + const Transform trans; };//end class +// template > +// // Qt_layer_show_points* +// void* +// make_show_points_layer(C *&container, +// It (C::*)() const begin, +// It (C::*)() const end, +// Transform& t = Transform(), +// Color c = CGAL::GREEN, +// int pointsize = 3, +// PointStyle pointstyle = CGAL::DISC) +// { +// return new CGAL::Qt_layer_show_points +// (container, begin, end, t, c, pointsize, pointstyle); +// }; + } // namespace CGAL #endif // CGAL_QT_LAYER_SHOW_POINTS_H diff --git a/Packages/Mesh_2/demo/Mesh_2/makefile b/Packages/Mesh_2/demo/Mesh_2/makefile index 6e24849998f..1ff71e5e095 100644 --- a/Packages/Mesh_2/demo/Mesh_2/makefile +++ b/Packages/Mesh_2/demo/Mesh_2/makefile @@ -13,8 +13,7 @@ include $(CGAL_MAKEFILE) # compiler flags #---------------------------------------------------------------------# -CXXFLAGS = -I../../../Qt_widget/include \ - -I../../include \ +CXXFLAGS = -I../../include \ $(CGAL_CXXFLAGS) \ $(LONG_NAME_PROBLEM_CXXFLAGS) \ $(DEBUG_OPT) @@ -35,6 +34,7 @@ LDFLAGS = \ #---------------------------------------------------------------------# all: \ + icons$(EXE_EXT) \ mesh_demo$(EXE_EXT) icons$(EXE_EXT): icons$(OBJ_EXT) @@ -43,10 +43,14 @@ icons$(EXE_EXT): icons$(OBJ_EXT) mesh_demo.moc: mesh_demo.C $(QT_MOC) -o mesh_demo.moc $< +Show_clusters.moc: Show_clusters.h + $(QT_MOC) -o Show_clusters.moc $< + mesh_demo$(OBJ_EXT): mesh_demo.moc \ + Show_clusters.moc \ ../../include/CGAL/Double_map.h \ - ../../include/CGAL/Filtred_circulator.h \ - ../../include/CGAL/Filtred_container.h \ + ../../include/CGAL/Filter_circulator.h \ + ../../include/CGAL/Filtered_container.h \ ../../include/CGAL/Mesh_2.h \ ../../include/CGAL/Mesh_default_traits_2.h \ ../../include/CGAL/Mesh_face_base_2.h \ diff --git a/Packages/Mesh_2/demo/Mesh_2/mesh_demo.C b/Packages/Mesh_2/demo/Mesh_2/mesh_demo.C index 4526acda1c3..5cfa0fea931 100644 --- a/Packages/Mesh_2/demo/Mesh_2/mesh_demo.C +++ b/Packages/Mesh_2/demo/Mesh_2/mesh_demo.C @@ -15,6 +15,10 @@ int main(int, char*) #include #include +#include +#include +#include +#include #include #include @@ -25,6 +29,8 @@ int main(int, char*) #include #include +#include + #include #include #include @@ -79,6 +85,7 @@ typedef K::Circle_2 Circle; typedef CGAL::Polygon_2 CGALPolygon; typedef CGAL::Mesh_2 Mesh; +typedef Mesh::Vertex Vertex; template class Show_marked_faces : public CGAL::Qt_widget_layer @@ -130,36 +137,13 @@ public: }; }; -namespace { - - class Seeds_wrapper_point : public Point { - public: - Seeds_wrapper_point(const Point& p): Point(p) {}; - - inline Point point() const { return *this; }; - }; - - class Seeds : public std::list { - public: - typedef std::list List; - typedef List::const_iterator const_iterator; - typedef const_iterator Vertex_iterator; - - inline const_iterator vertices_begin() { return begin(); } - inline const_iterator vertices_end() { return end(); } - }; - -}; - class MyWindow : public QMainWindow { Q_OBJECT public: MyWindow() : is_mesh_initialized(false), traits() { - triangulation = new Tr(traits); - mesh = 0; - seeds = new Seeds(); + mesh = new Mesh(traits); QFrame* mainframe = new QFrame(this, "mainframe"); QHBoxLayout *hbox = new QHBoxLayout(mainframe, 0, 0, @@ -198,22 +182,32 @@ public: // statusBar()->addWidget(aspect_ratio_label, 0, true); // LAYERS - show_points = new CGAL::Qt_layer_show_points(triangulation); - show_seeds = new CGAL::Qt_layer_show_points(seeds, - CGAL::BLUE, - 5, - CGAL::CROSS); + + show_points = + new Show_points_from_triangulation(mesh, + &Mesh::finite_vertices_begin, + &Mesh::finite_vertices_end, + std::mem_fun_ref(&Vertex::point)); + + show_seeds = new Show_seeds(mesh, + &Mesh::seeds_begin, + &Mesh::seeds_end, + CGAL::BLUE, + 5, + CGAL::CROSS); show_triangulation = - new CGAL::Qt_layer_show_triangulation(triangulation); + new CGAL::Qt_layer_show_triangulation(mesh); show_marked = - new Show_marked_faces(triangulation); + new Show_marked_faces(mesh); show_constraints = - new CGAL::Qt_layer_show_triangulation_constraints - (triangulation); + new CGAL::Qt_layer_show_triangulation_constraints + (mesh); show_circles = - new CGAL::Qt_layer_show_circles(triangulation); + new CGAL::Qt_layer_show_circles(mesh); show_mouse = new CGAL::Qt_layer_mouse_coordinates(*this); + show_clusters = new Show_Clusters(*mesh); + // layers order, first attached are "under" last attached widget->attach(show_marked); widget->attach(show_triangulation); @@ -450,7 +444,11 @@ public: CTRL+Key_L ); menuBar()->insertItem("&Advanced", this, SLOT(advanced())); - ; + + connect(this, SLOT(insertedInput()), + this, SIGNAL(set_not_initialized())); + connect(this, SLOT(initializedMesh()), + this, SIGNAL(set_initialized())); widget->set_window(-1.,1.,-1.,1.); widget->setMouseTracking(TRUE); @@ -460,11 +458,11 @@ public: void bounds(FT &xmin, FT &ymin, FT &xmax, FT &ymax) { - Mesh::Finite_vertices_iterator vi=triangulation->finite_vertices_begin(); + Mesh::Finite_vertices_iterator vi=mesh->finite_vertices_begin(); xmin=xmax=vi->point().x(); ymin=ymax=vi->point().y(); vi++; - while(vi != triangulation->finite_vertices_end()) + while(vi != mesh->finite_vertices_end()) { if(vi->point().x() < xmin) xmin=vi->point().x(); if(vi->point().x() > xmax) xmax=vi->point().x(); @@ -474,35 +472,21 @@ public: } } - void switch_to_mesh() - { - if(mesh == 0) - { - // std::cerr << "switch_to_mesh()" << std::endl; - // "construct" a mesh from triangulation without refining it - // (because of the third parameter dont_refine=true) - mesh = new Mesh(*triangulation, traits, true); - delete triangulation; - triangulation = mesh; - }; - } - - void switch_to_triangulation() - { - if(mesh != 0) - { - // std::cerr << "switch_to_triangulation()" << std::endl - // << mesh->number_of_vertices() << std::endl; - triangulation = new Tr(); // empty triangulation - triangulation->swap(*mesh); // swap *triangulation and *mesh - // *mesh is now empty - delete mesh; - mesh = 0; // in case delete does not the job - }; - } +signals: + void insertedInput(); + void initializedMesh(); public slots: + void set_initialized() + { + is_mesh_initialized = true; + } + + void set_not_initialized() + { + is_mesh_initialized = false; + } void get_cgal_object(CGAL::Object obj) { // TODO: when inputs arise, seeds should be cleared @@ -512,7 +496,6 @@ public slots: if(CGAL::assign(p,obj)) if(follow_mouse->is_active()) { - switch_to_mesh(); typedef Mesh::Face_handle Face_handle; std::list l; @@ -528,66 +511,45 @@ public slots: if(traits.is_bad_object().operator()(a, b, c)) l.push_back(fh); } - mesh->set_geom_traits(traits, l.begin(), l.end()); + mesh->set_geom_traits(traits); + mesh->set_bad_faces(l.begin(), l.end()); + mesh->calculate_bad_faces(); while( mesh->refine_step() ); widget->redraw(); } else if(get_seed->is_active()) { - seeds->push_back(p); - switch_to_mesh(); - mesh->mark_facets(seeds->begin(), seeds->end()); + Mesh::Seeds seeds; + std::copy(mesh->seeds_begin(), mesh->seeds_end(), + std::back_inserter(seeds)); + seeds.push_back(p); + mesh->set_seeds(seeds.begin(), seeds.end()); + mesh->mark_facets(); } else // get_point is active { - switch_to_triangulation(); - triangulation->insert(p); - switch_to_mesh(); - mesh->mark_facets(seeds->begin(), seeds->end()); - is_mesh_initialized = false; + mesh->insert(p); + mesh->mark_facets(); + emit( insertedInput() ); } else if (CGAL::assign(poly,obj)) { - switch_to_triangulation(); for(CGALPolygon::Edge_const_iterator it=poly.edges_begin(); it!=poly.edges_end(); it++) - triangulation->insert((*it).source(),(*it).target()); - switch_to_mesh(); - mesh->mark_facets(seeds->begin(), seeds->end()); - is_mesh_initialized = false; + mesh->insert((*it).source(),(*it).target()); + mesh->mark_facets(); + emit( insertedInput() ); } updatePointCouter(); widget->redraw(); } -// void redraw_win() -// { -// widget->lock(); -// widget->clear(); - -// for(Mesh::Edge_iterator it=mesh->edges_begin(); -// it!=mesh->edges_end(); -// it++) -// if(mesh->is_constrained(*it)) -// *widget << CGAL::RED << mesh->segment(*it); -// else -// *widget << CGAL::BLUE << mesh->segment(*it); -// *widget << CGAL::GREEN; -// for(Mesh::Vertex_iterator it=mesh->vertices_begin(); -// it!=mesh->vertices_end(); -// it++) -// *widget << it->point(); -// widget->unlock(); -// } - - //insert a bounding box around the mesh + //insert a bounding box around the mesh void insert_bounding_box() { - switch_to_triangulation(); - FT xmin, xmax, ymin, ymax; bounds(xmin, ymin, xmax, ymax); @@ -600,29 +562,28 @@ public slots: Point bb2(xcenter + 1.5*xspan, ycenter - 1.5*yspan); Point bb3(xcenter + 1.5*xspan, ycenter + 1.5*yspan); Point bb4(xcenter - 1.5*xspan, ycenter + 1.5*yspan); - triangulation->insert(bb1); - triangulation->insert(bb2); - triangulation->insert(bb3); - triangulation->insert(bb4); - triangulation->insert(bb1, bb2); - triangulation->insert(bb2, bb3); - triangulation->insert(bb3, bb4); - triangulation->insert(bb4, bb1); - switch_to_mesh(); - mesh->mark_facets(seeds->begin(), seeds->end()); + mesh->insert(bb1); + mesh->insert(bb2); + mesh->insert(bb3); + mesh->insert(bb4); + mesh->insert(bb1, bb2); + mesh->insert(bb2, bb3); + mesh->insert(bb3, bb4); + mesh->insert(bb4, bb1); + mesh->mark_facets(); + emit( insertedInput() ); widget->redraw(); } void updatePointCouter() { - nbre_de_points->setNum(triangulation->number_of_vertices()); + nbre_de_points->setNum(mesh->number_of_vertices()); } void refineMesh() { - switch_to_mesh(); saveTriangulationUrgently("last_input.edg"); - mesh->refine(seeds->begin(), seeds->end()); + mesh->refine(); is_mesh_initialized=true; updatePointCouter(); widget->redraw(); @@ -630,14 +591,13 @@ public slots: void conformMesh() { - switch_to_mesh(); if(!is_mesh_initialized) { saveTriangulationUrgently("last_input.edg"); - mesh->init(seeds->begin(), seeds->end()); + mesh->init(); is_mesh_initialized=true; } - mesh->conform(); + mesh->conform(CGAL::Gabriel_conform_inside_policy_2()); updatePointCouter(); widget->redraw(); } @@ -645,10 +605,10 @@ public slots: void refineMeshStep() { int counter = step_lenght; - switch_to_mesh(); if(!is_mesh_initialized) { - mesh->init(seeds->begin(), seeds->end()); + mesh->init(); + std::cerr << "mesh->init();"<< std::endl; is_mesh_initialized=true; saveTriangulationUrgently("last_input.edg"); } @@ -692,11 +652,8 @@ public slots: void clearMesh() { - switch_to_mesh(); mesh->clear(); - seeds->clear(); - is_mesh_initialized=false; - switch_to_triangulation(); + emit( insertedInput() ); updatePointCouter(); widget->clear_history(); widget->redraw(); @@ -704,14 +661,13 @@ public slots: void clearSeeds() { - seeds->clear(); - mesh->mark_facets(seeds->begin(), seeds->end()); + mesh->clear_seeds(); + mesh->mark_facets(); widget->redraw(); } void openTriangulation() { - switch_to_mesh(); QString s( QFileDialog::getOpenFileName( QString::null, my_filters, this ) ); if ( s.isEmpty() ) @@ -719,13 +675,8 @@ public slots: std::ifstream f(s); if(s.right(5) == ".poly") { - mesh->read_poly(f, true); - seeds->clear(); - for(Mesh::Seeds::iterator it = mesh->seeds.begin(); - it != mesh->seeds.end(); - ++it) - seeds->push_back(*it); - is_mesh_initialized=true; + read_poly(*mesh, f); + emit( insertedInput() ); } else if(s.right(5) == ".data") { @@ -810,13 +761,13 @@ public slots: mesh->insert_constraint(tl, bl); clearSeeds(); - is_mesh_initialized=false; + emit( insertedInput() ); } else { - mesh->read(f, true); + read(*mesh, f); clearSeeds(); - is_mesh_initialized=false; + emit( insertedInput() ); } // compute bounds @@ -838,23 +789,21 @@ public slots: void saveTriangulation() { - switch_to_mesh(); QString s( QFileDialog::getSaveFileName( "filename.edg", my_filters, this ) ); if ( s.isEmpty() ) return; std::ofstream of(s); if(s.right(5) == ".poly") - mesh->write_poly(of); + write_poly(*mesh, of); else - mesh->write(of); + write(*mesh, of); } void saveTriangulationUrgently(QString s=QString("dump.edg")) { - switch_to_mesh(); std::ofstream of(s); - mesh->write(of); + write(*mesh, of); } inline @@ -891,7 +840,6 @@ public slots: void setLocal() { - switch_to_mesh(); traits.set_local_size(!mesh->geom_traits().is_local_size()); pmCriteria->setItemChecked(menu_id, traits.is_local_size()); mesh->set_geom_traits(traits); @@ -913,11 +861,8 @@ private: static const QString my_filters; bool is_mesh_initialized; Meshtraits traits; - Tr* triangulation; Mesh* mesh; - Seeds* seeds; - QPopupMenu *pmCriteria; int menu_id; @@ -927,13 +872,22 @@ private: CGAL::Qt_widget_get_polygon* get_polygon; Follow_mouse* follow_mouse; - CGAL::Qt_layer_show_points* show_points; - CGAL::Qt_layer_show_points* show_seeds; - CGAL::Qt_layer_show_triangulation* show_triangulation; - CGAL::Qt_layer_show_triangulation_constraints* show_constraints; - CGAL::Qt_layer_show_circles* show_circles; + typedef CGAL::Qt_layer_show_points > + Show_points_from_triangulation; + + typedef CGAL::Qt_layer_show_points + Show_seeds; + + Show_points_from_triangulation* show_points; + Show_seeds* show_seeds; + CGAL::Qt_layer_show_triangulation* show_triangulation; + CGAL::Qt_layer_show_triangulation_constraints* show_constraints; + CGAL::Qt_layer_show_circles* show_circles; CGAL::Qt_layer_mouse_coordinates* show_mouse; - Show_marked_faces* show_marked; + Show_marked_faces* show_marked; + + Show_clusters* show_clusters; // QLabel* aspect_ratio_label; QLabel *nbre_de_points; diff --git a/Packages/Mesh_2/examples/Mesh_2/Makefile b/Packages/Mesh_2/examples/Mesh_2/Makefile index 4f0c49989f0..307e6659896 100644 --- a/Packages/Mesh_2/examples/Mesh_2/Makefile +++ b/Packages/Mesh_2/examples/Mesh_2/Makefile @@ -34,20 +34,30 @@ LDFLAGS = \ #---------------------------------------------------------------------# all: \ + conform$(EXE_EXT) \ mesh$(EXE_EXT) +conform$(OBJ_EXT): ../../include/CGAL/Conform_2.h \ + ../../include/CGAL/Mesh_default_traits_2.h \ + ../../include/CGAL/Double_map.h \ + ../../include/CGAL/Filter_circulator.h \ + ../../include/CGAL/Filtered_container.h + +conform$(EXE_EXT): conform$(OBJ_EXT) + $(CGAL_CXX) $(LIBPATH) $(EXE_OPT)conform conform$(OBJ_EXT) $(LDFLAGS) + mesh$(OBJ_EXT): ../../include/CGAL/Mesh_2.h \ ../../include/CGAL/Mesh_face_base_2.h \ ../../include/CGAL/Mesh_default_traits_2.h \ ../../include/CGAL/Double_map.h \ - ../../include/CGAL/Filtred_circulator.h \ - ../../include/CGAL/Filtred_container.h \ + ../../include/CGAL/Filter_circulator.h \ + ../../include/CGAL/Filtered_container.h \ mesh$(EXE_EXT): mesh$(OBJ_EXT) $(CGAL_CXX) $(LIBPATH) $(EXE_OPT)mesh mesh$(OBJ_EXT) $(LDFLAGS) clean: \ - mesh.clean + mesh.clean conform.clean #---------------------------------------------------------------------# # suffix rules @@ -56,3 +66,5 @@ clean: \ .C$(OBJ_EXT): $(CGAL_CXX) $(CXXFLAGS) $(OBJ_OPT) $< +%.h.gch: %.h + $(CGAL_CXX) $(CXXFLAGS) $(OBJ_OPT) $< diff --git a/Packages/Mesh_2/examples/Mesh_2/conform.C b/Packages/Mesh_2/examples/Mesh_2/conform.C new file mode 100644 index 00000000000..eb13d1b933d --- /dev/null +++ b/Packages/Mesh_2/examples/Mesh_2/conform.C @@ -0,0 +1,129 @@ +#include +#include +#include +#include +#include + +#include +#include + +#include + +#include +#include +#include + +typedef CGAL::Simple_cartesian K1; +typedef CGAL::Filtered_kernel K2; +struct K : public K2 {}; + +typedef CGAL::Triangulation_vertex_base_2 Vb; +typedef CGAL::Constrained_triangulation_face_base_2 Fb; +typedef CGAL::Triangulation_data_structure_2 Tds; +typedef CGAL::Mesh_default_traits_2 Meshtraits; +typedef CGAL::Constrained_Delaunay_triangulation_2 Tr; + +typedef CGAL::Conform_triangulation_2 Conform; + +typedef K::Point_2 Point; + +Conform conform; + +void usage(char** argv) +{ + std::cerr << "Usage: " << std::endl + << argv[0] << " [-Q] input.poly [output.poly]" << std::endl; +} + +int main(int argc, char** argv) +{ + int arg_count = 1; + bool terminal_output = true; + + if(argc < 2) + { + usage(argv); + return 1; + } + + while(argv[arg_count][0] == '-' && argv[arg_count] != "--") + { + if(std::string(argv[arg_count]) == "-Q") + terminal_output = false; + else + { + std::cerr << "Unknown option " << argv[arg_count] << std::endl; + usage(argv); + return 1; + } + ++arg_count; + } + if(argv[arg_count] == "--") + ++arg_count; + + if(argc < arg_count+1 || argc > arg_count+2) + { + usage(argv); + return 1; + }; + std::ifstream input(argv[arg_count]); + if(input) + { + read_poly(conform, input); + + conform.init(CGAL::Gabriel_conform_policy_2()); + + typedef Conform::cluster_vertices_iterator Vh_iterator; + for(Vh_iterator it = conform.clusters_vertices_begin(); + it != conform.clusters_vertices_end(); + ++it) + { + typedef std::pair It_pair; + + int n = conform.number_of_clusters_at_vertex(*it); + + std::cout << "Point(" << (*it)->point() << ")" << std::endl + << " " << n + << " cluster(s)." << std::endl; + + for(int i = 0; ipoint() << std::endl; + } + }; + + conform.gabriel_conform(); + } + else + { + std::cerr << "Bad file: " << argv[arg_count] << std::endl; + usage(argv); + return 1; + } + + if(argc==arg_count+1) + { + if(terminal_output) + write_poly(conform, std::cout); + } + else + { + std::ofstream output(argv[arg_count+1]); + write_poly(conform, output); + } + if(terminal_output) + std::cerr + << "Mesh points: " << conform.number_of_vertices() << std::endl + << "Mesh triangles: " << conform.number_of_faces () << std::endl; + + return 0; +}; diff --git a/Packages/Mesh_2/examples/Mesh_2/mesh.C b/Packages/Mesh_2/examples/Mesh_2/mesh.C index 4f563aee5a0..fbe2d059e2a 100644 --- a/Packages/Mesh_2/examples/Mesh_2/mesh.C +++ b/Packages/Mesh_2/examples/Mesh_2/mesh.C @@ -1,4 +1,4 @@ - #include +#include #include #include #include @@ -13,6 +13,8 @@ #include #include +#include + typedef CGAL::Simple_cartesian K1; typedef CGAL::Filtered_kernel K2; struct K : public K2 {}; @@ -88,7 +90,10 @@ int main(int argc, char** argv) }; std::ifstream input(argv[arg_count]); if(input) - mesh.read_poly(input); + { + read_poly(mesh, input); + mesh.refine(); + } else { std::cerr << "Bad file: " << argv[arg_count] << std::endl; @@ -99,12 +104,12 @@ int main(int argc, char** argv) if(argc==arg_count+1) { if(terminal_output) - mesh.write_poly(std::cout); + write_poly(mesh, std::cout); } else { std::ofstream output(argv[arg_count+1]); - mesh.write_poly(output); + write_poly(mesh, output); } if(terminal_output) std::cerr diff --git a/Packages/Mesh_2/include/CGAL/Conform_2.h b/Packages/Mesh_2/include/CGAL/Conform_2.h new file mode 100644 index 00000000000..4e0ed66a019 --- /dev/null +++ b/Packages/Mesh_2/include/CGAL/Conform_2.h @@ -0,0 +1,1140 @@ +#ifndef CGAL_CONFORM_2_H +#define CGAL_CONFORM_2_H +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#ifdef CGAL_USE_BOOST +#include +#endif + +CGAL_BEGIN_NAMESPACE + +class All_is_conform_policy_2 +{ +public: + template + struct Is_acceptable_in_a_cluster + { + bool operator()(const CTr&, + const typename CTr::Face_handle&) const + { + return true; + } + }; + + template + Is_acceptable_in_a_cluster is_acceptable_in_a_cluster_object() const + { + return Is_acceptable_in_a_cluster(); + } + + template + struct Is_locally_conform + { + bool operator()(const CTr&, + const typename CTr::Face_handle&, + const int) const + { + return true; + } + }; + + template + Is_locally_conform is_locally_conform_object() const + { + return Is_locally_conform(); + } +}; + +class Gabriel_conform_policy_2 : public All_is_conform_policy_2 +{ +public: + template + struct Is_locally_conform + { + bool operator()(const CTr& ct, + const typename CTr::Face_handle& fh, + const int i) const + { + typedef typename CTr::Geom_traits Geom_traits; + typedef typename Geom_traits::Angle_2 Angle_2; + typedef typename CTr::Point Point; + + const Angle_2 angle = ct.geom_traits().angle_2_object(); + + const Point& a = fh->vertex(ct.cw(i))->point(); + const Point& b = fh->vertex(ct.ccw(i))->point(); + const Point& c = fh->vertex(i)->point(); + const Point& d = fh->mirror_vertex(i)->point(); + + return( angle(a, c, b) != OBTUSE && + angle(a, d, b) != OBTUSE ); + } + }; + + template + Is_locally_conform is_locally_conform_object() const + { + return Is_locally_conform(); + } +}; + +class Delaunay_conform_policy_2 : public All_is_conform_policy_2 +{ +public: + template + struct Is_locally_conform + { + bool operator()(const CTr& ct, + const typename CTr::Face_handle& fh, + const int i) const + { + typedef typename CTr::Geom_traits Geom_traits; + typedef typename Geom_traits::Side_of_oriented_circle_2 + Side_of_oriented_circle_2; + typedef typename CTr::Point Point; + + const Side_of_oriented_circle_2 in_circle = + ct.geom_traits.side_of_oriented_circle_2_object(); + + const Point& a = fh->vertex(ct.cw(i))->point(); + const Point& b = fh->vertex(ct.ccw(i))->point(); + const Point& c = fh->vertex(i)->point(); + const Point& d = fh->mirror_vertex(i)->point(); + + return( in_circle(c, b, a, d) == ON_POSITIVE_SIDE ); + } + }; + + template + Is_locally_conform is_locally_conform_object() const + { + return Is_locally_conform(); + } +}; + +/** + - Tr is a Delaunay constrained triangulation (with intersections or + not). +*/ +template +class Conform_triangulation_2: public Tr +{ +public: + // -- public typedef -- + typedef Tr Triangulation; + typedef Conform_triangulation_2 Self; + + typedef typename Tr::Geom_traits Geom_traits; + typedef typename Tr::Triangulation_data_structure Tds; + + typedef typename Geom_traits::FT FT; + typedef FT Squared_length; + + // -- types inherited from the templated base class -- + // must be redefined for use in the implementation + typedef typename Tr::Vertex Vertex; + typedef typename Tr::Vertex_handle Vertex_handle; + typedef typename Tr::Vertex_circulator Vertex_circulator; + typedef typename Tr::Finite_vertices_iterator Finite_vertices_iterator; + + typedef typename Tr::Finite_edges_iterator Finite_edges_iterator; + + typedef typename Tr::Face_handle Face_handle; + typedef typename Tr::Face_circulator Face_circulator; + typedef typename Tr::Finite_faces_iterator Finite_faces_iterator; + typedef typename Tr::All_faces_iterator All_faces_iterator; + + typedef typename Tr::Locate_type Locate_type; + + typedef typename Tr::Point Point; + +protected: + // typedefs for private members types + typedef std::pair Constrained_edge; + + /** + Is_edge_constrained: function object class. + Take a Conform_triangulation_2 object m and a Vertex_handle v in + constructor. Is_edge_constrained(Vertex_handle v2) tells if + [v,v2] is a constrained edge in m. */ + class Is_edge_constrained { + Self* _m; + Vertex_handle _v; + public: + Is_edge_constrained(Self* m, Vertex_handle v) + : _m(m), _v(v) {} + + bool operator()(/*const*/ Vertex& v2) const + { + if(_m->is_infinite(v2.handle())) + return false; + Face_handle fh; + int i; + _m->is_edge(_v,v2.handle(),fh,i); + return fh->is_constrained(i); + } + }; + + class Is_really_a_constrained_edge { + const Self& _m; + public: + explicit Is_really_a_constrained_edge(const Self& m) : _m(m) {}; + bool operator()(const Constrained_edge& ce) const + { + Face_handle fh; + int i; + return _m.is_edge(ce.first, ce.second, fh,i) && + fh->is_constrained(i); + } + }; + + typedef Filtred_circulator + Constrained_vertex_circulator; + + typedef std::list List_of_constraints; + typedef CGAL::Filtered_container + Constrained_edges_queue; + + // Cluster register several informations about clusters. + // A cluster is a is a set of vertices v_i incident to one vertice + // v_0, so that angles between segments [v_0, v_i] is less than 60°. + struct Cluster { + bool reduced ; // Is the cluster reduced + + // smallest_angle gives the two vertices defining the + // smallest angle in the cluster + std::pair smallest_angle; + + FT rmin; // WARNING: rmin has no meaning if reduced=false!!! + Squared_length minimum_squared_length; + + // The following map tells what vertices are in the cluster and if + // the corresponding segment has been splitted once. + typedef std::map Vertices_map; + Vertices_map vertices; + + bool is_reduced() const { + return reduced; + } + + bool is_reduced(const Vertex_handle v) { + return vertices[v]; + } + }; + + typedef std::multimap Cluster_map; + +#ifdef CGAL_USE_BOOST +private: + template + struct Pair_get_first: public std::unary_function + { + typedef typename Pair::first_type result; + const result& operator()(const Pair& p) const + { + return p.first; + } + }; + +public: + typedef typename boost::projection_iterator_generator< + Pair_get_first, + typename Cluster_map::const_iterator>::type + Cluster_vertices_iterator; + +private: + typedef typename Cluster::Vertices_map Cluster_vertices_map; + +public: + typedef typename boost::projection_iterator_generator< + Pair_get_first, + typename Cluster_vertices_map::const_iterator>::type + Vertices_in_cluster_iterator; +#endif + +public: + // --- CONSTRUCTORS --- + + /** default constructor */ + explicit + Conform_triangulation_2(const Geom_traits& gt = Geom_traits()); + + // --- SURCHARGED INSERTION-DELETION FONCTIONS --- + // TODO! + + // --- ASSIGNEMENT --- + // TODO! + + // --- ACCESS FUNCTIONS --- + unsigned int number_of_constrained_edges() const; +#ifdef CGAL_USE_BOOST + Cluster_vertices_iterator clusters_vertices_begin() const + { + return Cluster_vertices_iterator(cluster_map.begin(), + Pair_get_first()); + } + + Cluster_vertices_iterator clusters_vertices_end() const + { + return Cluster_vertices_iterator(cluster_map.end(), + Pair_get_first()); + } + + unsigned int number_of_clusters_at_vertex(Vertex_handle vh) + { + typedef typename Cluster_map::iterator Iterator; + typedef std::pair Range; + Range range = cluster_map.equal_range(vh); + return std::distance(range.first, range.second); + } + + // returns the sequence of vertices bellonging to the n-th cluster of vh + std::pair + vertices_in_cluster_sequence(Vertex_handle vh, unsigned int n) + { + typedef typename Cluster_map::iterator Iterator; + typedef std::pair Range; + typedef typename Range::first_type Clusters_iterator; + typedef Pair_get_first + Get_first; + + Range range = cluster_map.equal_range(vh); + Iterator first = range.first; + std::advance(first, n); + const Cluster& c = first->second; + + return std::make_pair(Vertices_in_cluster_iterator(c.vertices.begin(), + Get_first()), + Vertices_in_cluster_iterator(c.vertices.end(), + Get_first())); + } + +#endif + + // --- HELPING FUNCTION --- + void clear(); + + // --- CONFORMING FUNCTIONS --- + + // Conform edges + void gabriel_conform(); + void delaunay_conform(); + + // -- other conform policies -- + template + void conform(const Conform_policy& conf_policy) + { + while( !edges_to_be_conformed.empty() ) + { + process_one_edge(conf_policy); + }; + } + + // --- STEP BY STEP FUNCTIONS --- + + /** + init(...): Initialize the data structures + (The call of this function is REQUIRED before any step by step + operation). + */ + template + void init(const Conform_policy&); + + /** Execute on step of the algorithm. + init() should have been called before. + */ + template + bool refine_step(const Conform_policy&); + + /** Tells if all constrained edges are conformed. */ + bool is_conformed() + // This function cannot be "const" because, as edges_to_be_conformed is + // filtred, its empty() method is not const. + { return edges_to_be_conformed.empty(); } + + // --- PROTECTED TYPES --- +protected: + typedef std::list List_of_face_handles; + + // -- traits type -- + typedef typename Geom_traits::Vector_2 Vector_2; + typedef typename Geom_traits::Construct_translated_point_2 + Construct_translated_point_2; + typedef typename Geom_traits::Compute_squared_distance_2 + Compute_squared_distance_2; + typedef typename Geom_traits::Angle_2 Angle_2; + typedef typename Geom_traits::Construct_vector_2 + Construct_vector_2; + typedef typename Geom_traits::Construct_scaled_vector_2 + Construct_scaled_vector_2; + typedef typename Geom_traits::Construct_midpoint_2 + Construct_midpoint_2; + typedef typename Geom_traits::Orientation_2 Orientation_2; + +private: + // PRIVATE MEMBER DATAS + + // edges_to_be_conformed: list of encroached constrained edges + // warning: some edges could be destroyed, use the same wrapper + // is_really_a_constrained_edge: tester to filter edges_to_be_conformed + const Is_really_a_constrained_edge is_really_a_constrained_edge; + Constrained_edges_queue edges_to_be_conformed; + + // Cluster_map: multimap Vertex_handle -> Cluster + // each vertex can have several clusters + Cluster_map cluster_map; + + // tell if init has been called since the last insertion of a + // constraint + bool initialized; + + // --- PRIVATE MEMBER FUNCTIONS --- +private: + + // -- auxiliary functions to handle clusters -- + + // for all vertices, call create_clusters_of_vertex + template + void create_clusters(const Conform_policy&); + + // compute clusters of the vertex v, using the auxiliary function + // construct_cluster + template + void create_clusters_of_vertex(const Vertex_handle v, + const Conform_policy&); + + // add the sequence [begin, end] to the cluster c and add it to the + // clusters of the vertex v + void construct_cluster(const Vertex_handle v, + Constrained_vertex_circulator begin, + const Constrained_vertex_circulator& end, + Cluster c = Cluster()); + + // see more functions about clusters in protected member functions + + + // -- functions that maintain the queue of encroached edges -- + + // scan all constrained edges and put them in the queue if they are + // encroached + template + void fill_edge_queue(const Conform_policy&); + + // update the queue with edges incident to vm + template + void update_edges_to_be_conformed(const Vertex_handle va, + const Vertex_handle vb, + const Vertex_handle vm, + const Conform_policy&); + + // -- inlined functions that compose the refinement process -- + + // take one edge in the queue and call refine_edge + template + void process_one_edge(const Conform_policy&); + + // handle the encroached edge, call cut_cluster_edge+update_cluster + // or insert_middle + template + void refine_edge(const Vertex_handle va, const Vertex_handle vb, + const Conform_policy&); + + // -- auxiliary functions that return a boolean -- + + // tell if [va,vb] is encroached, by looking for the two neighbors + // This function takes care of markers. +// bool is_encroached(const Vertex_handle va, +// const Vertex_handle vb) const; + +// // tell if [va,vb] is encroached by p +// bool is_encroached(const Vertex_handle va, +// const Vertex_handle vb, +// const Point& p) const; + + // tell if the angle is less than 60° + // Uses squared_cosine_of_angle_times_4 and used by + // create_clusters_of_vertex + bool is_small_angle(const Point& pleft, + const Point& pmiddle, + const Point& pright) const; + + // -- auxiliary functions that are called to split an edge or -- + // -- a face -- + + // cut [va,vb] knowing that it is in the cluster c + Vertex_handle cut_cluster_edge(const Vertex_handle va, + const Vertex_handle vb, + Cluster& c); + + // insert the midpoint of the edge (f,i) + Vertex_handle insert_middle(Face_handle f, const int i); + + + // -- functions that really insert points -- + + // virtual function that inserts the point p in the edge + // (fh,edge_index) + virtual + Vertex_handle insert_in_the_edge(Face_handle fh, + const int edge_index, + const Point& p); + + // -- helping computing functions -- + + // return the squared cosine of the angle + // times 4 + FT squared_cosine_of_angle_times_4(const Point& pleft, + const Point& pmiddle, + const Point& pright) const; + + // --- PROTECTED MEMBER FUNCTIONS --- +protected: + // -- functions to manage clusters from derived classes -- + + // update the cluster of [va,vb], putting vm instead of vb. If + // reduction=false, the edge [va,vm] is not set reduced. + void update_cluster(Cluster& c, const Vertex_handle va, + const Vertex_handle vb, const Vertex_handle vm, + bool reduction = true); + + // get_cluster returns the cluster of [va,vb] in c and return true + // if it is in a cluster. If erase=true, the cluster is remove from + // the cluster map + bool get_cluster(const Vertex_handle va, const Vertex_handle vb, + Cluster &c, bool erase = false); + + + void add_contrained_edge_to_be_conform(const Vertex_handle& va, + const Vertex_handle& vb) + { + edges_to_be_conformed.push_back(Constrained_edge(va,vb)); + } +}; // end of Conform_triangulation_2 + +// --- CONSTRUCTORS --- + +template +Conform_triangulation_2:: +Conform_triangulation_2(const Geom_traits& gt) + : Tr(gt), is_really_a_constrained_edge(*this), + edges_to_be_conformed(is_really_a_constrained_edge), + cluster_map(), + initialized(false) +{}; + +// --- ACCESS FUNCTIONS --- + +template +unsigned int Conform_triangulation_2:: +number_of_constrained_edges() const +{ + int nedges = 0; + for(Finite_edges_iterator eit = finite_edges_begin(); + eit != finite_edges_end(); + ++eit) + if((*eit).first->is_constrained((*eit).second)) + ++nedges; + return nedges; +} + + +// --- HELPING FUNCTIONS --- + +template +void Conform_triangulation_2:: +clear() +{ + cluster_map.clear(); + edges_to_be_conformed.clear(); + Triangulation::clear(); +} + +// --- CONFORMING FUNCTIONS --- + +template +inline +void Conform_triangulation_2:: +gabriel_conform() +{ + if(!initialized) init(Gabriel_conform_policy_2()); + while( !edges_to_be_conformed.empty() ) + { + process_one_edge(Gabriel_conform_policy_2()); + }; +} + +template +inline +void Conform_triangulation_2:: +delaunay_conform() +{ + if(!initialized) init(Delaunay_conform_policy_2()); + while( !edges_to_be_conformed.empty() ) + { + process_one_edge(Delaunay_conform_policy_2()); + }; +} + +// --- STEP BY STEP FUNCTIONS --- +template +template +inline +void Conform_triangulation_2:: +init(const Conform_policy& conf_policy) +{ + cluster_map.clear(); + edges_to_be_conformed.clear(); + create_clusters(conf_policy); + fill_edge_queue(conf_policy); + initialized = true; +} + +template +template +inline +bool Conform_triangulation_2:: +refine_step(const Conform_policy& conf_policy) +{ + if( !edges_to_be_conformed.empty() ) + process_one_edge(conf_policy); + else + return false; + return true; +} + + + +// --- PRIVATE MEMBER FUNCTIONS --- + +template +template +void Conform_triangulation_2:: +create_clusters(const Conform_policy& conf_policy) +{ + for(Finite_vertices_iterator vit = finite_vertices_begin(); + vit != finite_vertices_end(); + vit++) + create_clusters_of_vertex(vit, conf_policy); +} + +template +template +void Conform_triangulation_2:: +create_clusters_of_vertex(const Vertex_handle v, + const Conform_policy& conf_policy) +{ + Is_edge_constrained test(this, v); + Constrained_vertex_circulator begin(incident_vertices(v),test); + // This circulator represents all constrained edges around the + // vertex v. An edge [v,v'] is represented by the vertex v'. + + if(begin == 0) return; // if there is only one vertex + + Constrained_vertex_circulator + current(begin), next(begin), cluster_begin(begin); + ++next; // next is always just after current. + if(current == next) return; + + bool in_a_cluster = false; + do + { + typedef Conform_triangulation_2 Self; + typedef typename Conform_policy:: + template Is_acceptable_in_a_cluster Is_face_acceptable; + const Is_face_acceptable is_acceptable = + conf_policy.template is_acceptable_in_a_cluster_object(); + // is_acceptable tells if a face can be in a cluster, according + // to markers for example. + + Face_handle f; + int i; + is_edge(v, next, f, i); // put in f the face on the right side + // of (v,next) + if(is_acceptable(*this, f) && + is_small_angle(current->point(), v->point(), next->point())) + { + if(!in_a_cluster) + { + // at this point, current is the beginning of a cluster + in_a_cluster = true; + cluster_begin = current; + } + } + else + if(in_a_cluster) + { + // at this point, current is the end of a cluster and + // cluster_begin is its beginning + construct_cluster(v, cluster_begin, current); + in_a_cluster = false; + } + ++next; + ++current; + } while( current!=begin ); + if(in_a_cluster) + { + Cluster c; + if(get_cluster(v, begin, c, true)) + // get the cluster and erase it from the clusters map + construct_cluster(v, cluster_begin, begin, c); + else + construct_cluster(v, cluster_begin, current); + } +} + +template +void Conform_triangulation_2:: +construct_cluster(Vertex_handle v, + Constrained_vertex_circulator begin, + const Constrained_vertex_circulator& end, + Cluster c) +{ + Compute_squared_distance_2 squared_distance = + geom_traits().compute_squared_distance_2_object(); + + if(c.vertices.empty()) + { + c.reduced = false; + // c.rmin is not initialized because + // reduced=false! + c.minimum_squared_length = + squared_distance(v->point(), begin->point()); + Constrained_vertex_circulator second(begin); + ++second; + c.smallest_angle.first = begin; + c.smallest_angle.second = second; + } + + bool all_edges_in_cluster=false; // tell if all incident edges are + // in the cluster + if(begin==end) + all_edges_in_cluster=true; + + const Point& vp = v->point(); + + FT greatest_cosine = + squared_cosine_of_angle_times_4(c.smallest_angle.first->point(), + v->point(), + c.smallest_angle.second->point()); + + Constrained_vertex_circulator next(begin); + ++next; + do + { + c.vertices[begin] = false; + Squared_length l = squared_distance(vp, + begin->point()); + c.minimum_squared_length = + std::min(l,c.minimum_squared_length); + + if(all_edges_in_cluster || begin!=end) + { + FT cosine = + squared_cosine_of_angle_times_4(begin->point(), + v->point(), + next->point()); + if(cosine>greatest_cosine) + { + greatest_cosine = cosine; + c.smallest_angle.first = begin; + c.smallest_angle.second = next; + } + } + } + while(next++,begin++!=end); + cluster_map.insert(std::make_pair(v,c)); +} + +template +template +void Conform_triangulation_2:: +fill_edge_queue(const Conform_policy& conf_policy) +{ + typedef typename Conform_policy:: + template Is_locally_conform Is_locally_conform; + + Is_locally_conform is_locally_conform = + conf_policy.template is_locally_conform_object(); + + for(Finite_edges_iterator ei = finite_edges_begin(); + ei != finite_edges_end(); + ++ei) + { + if((*ei).first->is_constrained((*ei).second) && + !is_locally_conform(*this, (*ei).first, (*ei).second) ) + { + const Vertex_handle& va = (*ei).first->vertex(cw((*ei).second)); + const Vertex_handle& vb = (*ei).first->vertex(ccw((*ei).second)); + edges_to_be_conformed.push_back(std::make_pair(va, vb)); + } + } +} + +//update the encroached segments list +// TODO: perhaps we should remove destroyed edges too +// TODO: rewrite this function one day +template +template +void Conform_triangulation_2:: +update_edges_to_be_conformed(Vertex_handle va, + Vertex_handle vb, + Vertex_handle vm, + const Conform_policy& conf_policy) +{ + typedef typename Conform_policy:: + template Is_locally_conform Is_locally_conform; + + Is_locally_conform is_locally_conform = + conf_policy.template is_locally_conform_object(); + + Face_circulator fc = incident_faces(vm), fcbegin(fc); + + do { + for(int i = 0; i<3; i++) { + if( fc->is_constrained(i) && !is_infinite(fc,cw(i)) && + !is_infinite(fc,ccw(i)) && !is_locally_conform(*this, fc->handle(), i) ) + { + const Vertex_handle& v1 = fc->vertex(ccw(i)); + const Vertex_handle& v2 = fc->vertex(cw(i)); + edges_to_be_conformed.push_back(Constrained_edge(v1, v2)); + } + } + ++fc; + } while(fc != fcbegin); + + Face_handle fh; + int index; + is_edge(va, vm, fh, index); + if(!is_locally_conform(*this, fh, index)) { + edges_to_be_conformed.push_back(Constrained_edge(va, vm)); + } + is_edge(vb, vm, fh, index); + if(!is_locally_conform(*this, fh, index)) { + edges_to_be_conformed.push_back(Constrained_edge(vb, vm)); + } +} + +template +template +inline +void Conform_triangulation_2:: +process_one_edge(const Conform_policy& conf_policy) +{ + Constrained_edge ce = edges_to_be_conformed.front(); + edges_to_be_conformed.pop_front(); + + refine_edge(ce.first, ce.second, conf_policy); +} + +//this function split all the segments that are encroached +template +template +void Conform_triangulation_2:: +refine_edge(Vertex_handle va, Vertex_handle vb, + const Conform_policy& conf_policy) +{ + Face_handle f; + int i; + is_edge(va, vb, f, i); // get the edge (f,i) + CGAL_assertion(f->is_constrained(i)); + + Cluster c,c2; + Vertex_handle vm; + + if( get_cluster(va,vb,c,true) ) + if( get_cluster(vb,va,c2,true) ) + { // both ends are clusters + vm = insert_middle(f,i); + update_cluster(c,va,vb,vm,false); + update_cluster(c2,vb,va,vm,false); + } + else + // va only is a cluster + vm = cut_cluster_edge(va,vb,c); + else + if( get_cluster(vb,va,c,true) ) + // vb only is a cluster + vm = cut_cluster_edge(vb,va,c); + else + // no cluster + vm = insert_middle(f,i); + update_edges_to_be_conformed(va, vb, vm, conf_policy); +}; + +// template +// inline +// bool Conform_triangulation_2:: +// is_encroached(const Vertex_handle va, const Vertex_handle vb) const +// { +// Face_handle fh; +// int i; +// is_edge(va,vb,fh,i); + +// const Point& candidat_1 = fh->vertex(i)->point(); +// const Point& candidat_2 = fh->mirror_vertex(i)->point(); + +// return ( (/* fh->is_marked() && */ +// is_encroached(va, vb, candidat_1) ) || +// (/* fh->neighbor(i)->is_marked() && */ +// is_encroached(va, vb, candidat_2) ) +// ); +// } + +// -> traits? +// TODO, FIXME: not robust! +template +bool Conform_triangulation_2:: +is_small_angle(const Point& pleft, + const Point& pmiddle, + const Point& pright) const +{ + Angle_2 angle = geom_traits().angle_2_object(); + Orientation_2 orient = geom_traits().orientation_2_object(); + + if( angle(pleft, pmiddle, pright)==OBTUSE ) + return false; + if( orient(pmiddle,pleft,pright)==RIGHTTURN) + return false; + + FT cos_alpha = squared_cosine_of_angle_times_4(pleft, pmiddle, + pright); + + if(cos_alpha > 1) + { + return true; //the same cluster + } + else + { + return false; //another cluster + } +} + +template +typename Conform_triangulation_2::Vertex_handle +Conform_triangulation_2:: +cut_cluster_edge(Vertex_handle va, Vertex_handle vb, Cluster& c) +{ + Construct_vector_2 vector = + geom_traits().construct_vector_2_object(); + Construct_scaled_vector_2 scaled_vector = + geom_traits().construct_scaled_vector_2_object(); + Compute_squared_distance_2 squared_distance = + geom_traits().compute_squared_distance_2_object(); + Construct_midpoint_2 midpoint = + geom_traits().construct_midpoint_2_object(); + Construct_translated_point_2 translate = + geom_traits().construct_translated_point_2_object(); + + Vertex_handle vc; + + if(c.is_reduced(vb)) + { + Face_handle fh; + int i; + is_edge(va,vb,fh,i); + vc = insert_middle(fh,i); + } + else + { + const Point + & a = va->point(), + & b = vb->point(), + & m = midpoint(a, b); + + + + Vector_2 v = vector(a,m); + v = scaled_vector(v,CGAL_NTS sqrt(c.minimum_squared_length / + squared_distance(a,b))); + + Point i = translate(a,v), i2(i); + + do { + i = translate(a,v); + v = scaled_vector(v,FT(2)); + i2 = translate(a,v); + } while(squared_distance(a,i2) <= squared_distance(a,m)); + if( squared_distance(i,m) > squared_distance(m,i2) ) + i = i2; + //here i is the best point for splitting + Face_handle fh; + int index; + is_edge(va,vb,fh,index); + + vc = insert_in_the_edge(fh, index, i); + } + update_cluster(c, va, vb, vc); + return vc; +} + +template +typename Conform_triangulation_2::Vertex_handle +Conform_triangulation_2:: +insert_middle(Face_handle f, int i) +{ + Construct_midpoint_2 + midpoint = geom_traits().construct_midpoint_2_object(); + + const Vertex_handle + & va = f->vertex(cw(i)), + & vb = f->vertex(ccw(i)); + + const Point& mp = midpoint(va->point(), vb->point()); + + Vertex_handle vm = insert_in_the_edge(f, i, mp); + + return vm; +} + +template +inline +typename Conform_triangulation_2::Vertex_handle +Conform_triangulation_2:: +insert_in_the_edge(Face_handle fh, int edge_index, const Point& p) + // insert the point p in the edge (fh, edge_index). It updates seeds + // too. +{ + List_of_face_handles zone_of_p; + + // deconstrain the edge before finding the conflicts + fh->set_constraint(edge_index,false); + fh->neighbor(edge_index)->set_constraint(fh->mirror_index(edge_index), + false); + + get_conflicts_and_boundary(p, + std::back_inserter(zone_of_p), + Emptyset_iterator(), fh); + + // reconstrain the edge + fh->set_constraint(edge_index,true); + fh->neighbor(edge_index)->set_constraint(fh->mirror_index(edge_index),true); + + Vertex_handle vp = insert(p, Triangulation::EDGE, fh, edge_index); + // TODO, WARNING: this is not robust! + // We should deconstrained the constrained edge, insert the two + // subconstraints and re-constrain them + + return vp; +} + +// # used by: is_small_angle, create_clusters_of_vertex +// # compute 4 times the square of the cosine of the angle (ab,ac) +// # WARNING, TODO: this is not exact with doubles and can lead to crashes!! +template +typename Conform_triangulation_2::FT +Conform_triangulation_2:: +squared_cosine_of_angle_times_4(const Point& pb, const Point& pa, + const Point& pc) const +{ + Compute_squared_distance_2 squared_distance = + geom_traits().compute_squared_distance_2_object(); + + const FT + a = squared_distance(pb, pc), + b = squared_distance(pa, pb), + c = squared_distance(pa, pc); + + const FT num = a-(b+c); + + return (num*num)/(b*c); +} + + +// --- PROTECTED MEMBER FUNCTIONS --- + +template +void Conform_triangulation_2:: +update_cluster(Cluster& c, Vertex_handle va,Vertex_handle vb, + Vertex_handle vm, bool reduction) +{ + Compute_squared_distance_2 squared_distance = + geom_traits().compute_squared_distance_2_object(); + c.vertices.erase(vb); + c.vertices[vm] = reduction; + + if(vb==c.smallest_angle.first) + c.smallest_angle.first = vm; + if(vb==c.smallest_angle.second) + c.smallest_angle.second = vm; + + FT l = squared_distance(va->point(),vm->point()); + if(lfirst)) + ++it; // TODO: use std::find and an object class + if(it==c.vertices.end()) + c.reduced = true; + } + + if(c.is_reduced()) + c.rmin = squared_distance(c.smallest_angle.first->point(), + c.smallest_angle.second->point())/FT(4); + cluster_map.insert(std::make_pair(va,c)); +} + +// # used by refine_face and cut_cluster_edge +template +bool Conform_triangulation_2:: +get_cluster(Vertex_handle va, Vertex_handle vb, Cluster &c, bool erase) +{ + typedef typename Cluster_map::iterator Iterator; + typedef std::pair Range; + Range range = cluster_map.equal_range(va); + for(Iterator it = range.first; it != range.second; it++) + { + Cluster &cl = it->second; + if(cl.vertices.find(vb)!=cl.vertices.end()) { + c = it->second; + if(erase) + cluster_map.erase(it); + return true; + } + } + return false; +} + + +// --- GLOBAL FUNCTIONS --- +template +void +gabriel_conform(Tr& t) +{ + typedef Conform_triangulation_2 Conform; + + Conform conform; + conform.swap(t); + conform.gabriel_conform(); + t.swap(conform); +} + +template +void +delaunay_conform(Tr& t) +{ + typedef Conform_triangulation_2 Conform; + + Conform conform; + conform.swap(t); + conform.delaunay_conform(); + t.swap(conform); +} + +CGAL_END_NAMESPACE + + +#endif diff --git a/Packages/Mesh_2/include/CGAL/Constrained_Delaunay_triangulation_2.h b/Packages/Mesh_2/include/CGAL/Constrained_Delaunay_triangulation_2.h deleted file mode 100644 index 2002b177dfb..00000000000 --- a/Packages/Mesh_2/include/CGAL/Constrained_Delaunay_triangulation_2.h +++ /dev/null @@ -1,664 +0,0 @@ -// ====================================================================== -// -// Copyright (c) 1997 The CGAL Consortium -// -// This software and related documentation is part of an INTERNAL release -// of the Computational Geometry Algorithms Library (CGAL). It is not -// intended for general use. -// -// ---------------------------------------------------------------------- -// -// release : $CGAL_Revision: CGAL-2.5-I-22 $ -// release_date : $CGAL_Date: 2002/09/03 $ -// -// file : include/CGAL/Constrained_Delaunay_triangulation_2.h -// package : Triangulation_2 (7.46) -// maintainer : Mariette Yvinec -// source : $RCSfile$ -// revision : $Revision$ -// revision_date : $Date$ -// author(s) : Mariette Yvinec, Jean Daniel Boissonnat -// -// coordinator : Mariette Yvinec < Mariette Yvinec@sophia.inria.fr> -// -// ====================================================================== - -#ifndef CGAL_CONSTRAINED_DELAUNAY_TRIANGULATION_2_H -#define CGAL_CONSTRAINED_DELAUNAY_TRIANGULATION_2_H - - -#include -#include -#include - -CGAL_BEGIN_NAMESPACE - - -template , - Constrained_triangulation_face_base_2 >, - class Itag = No_intersection_tag > -class Constrained_Delaunay_triangulation_2 - : public Constrained_triangulation_2 -{ -public: - typedef Constrained_triangulation_2 Ctr; - typedef Constrained_Delaunay_triangulation_2 CDt; - typedef typename Ctr::Geom_traits Geom_traits; - typedef typename Ctr::Intersection_tag Intersection_tag; - - typedef typename Ctr::Constraint Constraint; - typedef typename Ctr::Vertex_handle Vertex_handle; - typedef typename Ctr::Face_handle Face_handle; - typedef typename Ctr::Edge Edge; - typedef typename Ctr::Finite_faces_iterator Finite_faces_iterator; - typedef typename Ctr::Face_circulator Face_circulator; - typedef typename Ctr::Locate_type Locate_type; - - typedef typename Ctr::List_edges List_edges; - typedef typename Ctr::List_faces List_faces; - typedef typename Ctr::List_vertices List_vertices; - typedef typename Ctr::List_constraints List_constraints; - typedef typename Ctr::Less_edge less_edge; - typedef typename Ctr::Edge_set Edge_set; - - typedef typename Geom_traits::Point_2 Point; - - - Constrained_Delaunay_triangulation_2(const Geom_traits& gt=Geom_traits()) - : Ctr(gt) { } - - Constrained_Delaunay_triangulation_2(const CDt& cdt) - : Ctr(cdt) {} - - Constrained_Delaunay_triangulation_2(List_constraints& lc, - const Geom_traits& gt=Geom_traits()) - : Ctr(gt) - { - typename List_constraints::iterator itc = lc.begin(); - for( ; itc != lc.end(); ++itc) { - insert((*itc).first, (*itc).second); - } - CGAL_triangulation_postcondition( is_valid() ); - } - - template - Constrained_Delaunay_triangulation_2(InputIterator it, - InputIterator last, - const Geom_traits& gt=Geom_traits() ) - : Ctr(gt) - { - for ( ; it != last; it++) { - insert((*it).first, (*it).second); - } - CGAL_triangulation_postcondition( is_valid() ); - } - - virtual ~Constrained_Delaunay_triangulation_2() {} - - - // FLIPS - bool is_flipable(Face_handle f, int i) const; - void flip(Face_handle& f, int i); - void flip_around(Vertex_handle va); - void flip_around(List_vertices & new_vertices); - void propagating_flip(Face_handle& f,int i); - void propagating_flip(List_edges & edges); - - // CONFLICTS - bool test_conflict(Face_handle fh, const Point& p) const; //deprecated - bool test_conflict(const Point& p, Face_handle fh) const; - void find_conflicts(const Point& p, std::list& le, //deprecated - Face_handle hint= Face_handle(NULL)) const; - // //template member functions, declared and defined at the end - // template - // bool get_conflicts_and_boundary(const Point &p, - // Out_it1 fit, - // Out_it2 eit, - // Face_handle start) const; - // template - // bool get_conflicts(const Point &p, - // Out_it1 fit, - // Face_handle start ) const; - // template - // bool get_boundary_of_conflicts(const Point &p, - // Out_it2 eit, - // Face_handle start ) const; - - - // INSERTION-REMOVAL - Vertex_handle insert(const Point & a, Face_handle start = Face_handle(NULL)); - Vertex_handle insert(const Point& p, - Locate_type lt, - Face_handle loc, int li ); - Vertex_handle push_back(const Point& a); -// template < class InputIterator > -// int insert(InputIterator first, InputIterator last); - - void remove(Vertex_handle v); - void remove_incident_constraints(Vertex_handle v); - void remove_constraint(Face_handle f, int i); - - //for backward compatibility - void insert(Point a, Point b) { insert_constraint(a, b);} - void insert(Vertex_handle va, Vertex_handle vb) {insert_constraint(va,vb);} - - // CHECK - bool is_valid(bool verbose = false, int level = 0) const; - -protected: - virtual Vertex_handle virtual_insert(const Point & a, - Face_handle start = Face_handle(NULL)); - virtual Vertex_handle virtual_insert(const Point& a, - Locate_type lt, - Face_handle loc, - int li ); -//Vertex_handle special_insert_in_edge(const Point & a, Face_handle f, int i); - void remove_2D(Vertex_handle v ); - virtual void triangulate_hole(List_faces& intersected_faces, - List_edges& conflict_boundary_ab, - List_edges& conflict_boundary_ba); - -public: - // MESHING - // suppressed meshing functions from here - - //template member functions -public: - template < class InputIterator > -#if defined(_MSC_VER) || defined(__SUNPRO_CC) - int insert(InputIterator first, InputIterator last, int i = 0) -#else - int insert(InputIterator first, InputIterator last) -#endif - { - int n = number_of_vertices(); - while(first != last){ - insert(*first); - ++first; - } - return number_of_vertices() - n; - } - -#warning "File patched by Laurent Rineau, for Mesh_2 robustness!" - - template - bool - get_conflicts_and_boundary(const Point &p, - Out_it1 fit, - Out_it2 eit, - Face_handle start = Face_handle(NULL)) const - { - CGAL_triangulation_precondition( dimension() == 2); - int li; - Locate_type lt; - Face_handle fh = locate(p,lt,li, start); - return get_conflicts_and_boundary(p, fit, eit, lt, fh, li); - } - - template - bool - get_conflicts_and_boundary(const Point &p, - Out_it1 fit, - Out_it2 eit, - Locate_type lt, - Face_handle fh, - int li) const - { - CGAL_triangulation_precondition( dimension() == 2); - switch(lt) { - case OUTSIDE_AFFINE_HULL: - case VERTEX: - return false; - case EDGE: - // std::cerr << "EDGE!" << std::endl; - if( fh->is_constrained(li) ) - { - Face_handle neighbor = fh->neighbor(li); - if( ! is_infinite(neighbor) ) - { - *fit++ = neighbor; - propagate_conflicts(p,neighbor,0,fit,eit); - propagate_conflicts(p,neighbor,1,fit,eit); - propagate_conflicts(p,neighbor,2,fit,eit); - } - } - case FACE: - case OUTSIDE_CONVEX_HULL: - // std::cerr << "OUTSIDE_CONVEX_HULL!" << std::endl; - *fit++ = fh; //put fh in Out_it1 - propagate_conflicts(p,fh,0,fit,eit); - propagate_conflicts(p,fh,1,fit,eit); - propagate_conflicts(p,fh,2,fit,eit); - return true; - } - CGAL_triangulation_assertion(false); - return false; - } - - template - bool - get_conflicts(const Point &p, - Out_it1 fit, - Face_handle start= Face_handle(NULL)) const - { - return get_conflicts_and_boundary(p, fit, Emptyset_iterator(), start); - } - - template - inline bool - get_boundary_of_conflicts(const Point &p, - Out_it2 eit, - Face_handle start= Face_handle(NULL)) const - { - return get_conflicts_and_boundary(p, Emptyset_iterator(), eit, start); - } - -private: - template - void propagate_conflicts (const Point &p, - Face_handle fh, - int i, - Out_it1 fit, - Out_it2 eit) const - { - Face_handle fn = fh->neighbor(i); - if ( fh->is_constrained(i) || ! test_conflict(p,fn)) { - *eit++ = Edge(fn, fn->index(fh)); - return; - } - *fit++ = fn; - int j = fn->index(fh); - propagate_conflicts(p,fn,ccw(j),fit,eit); - propagate_conflicts(p,fn,cw(j),fit,eit); - return; - } - - -}; - - -template < class Gt, class Tds, class Itag > -bool -Constrained_Delaunay_triangulation_2:: -is_flipable(Face_handle f, int i) const - // determines if edge (f,i) can be flipped -{ - Face_handle ni = f->neighbor(i); - if (is_infinite(f) || is_infinite(ni)) return false; - if (f->is_constrained(i)) return false; - return (side_of_oriented_circle(ni, f->vertex(i)->point()) - == ON_POSITIVE_SIDE); -} - -template < class Gt, class Tds, class Itag > -void -Constrained_Delaunay_triangulation_2:: -flip(Face_handle& f, int i) -{ - // The following precondition prevents the test suit - // of triangulation to work on constrained Delaunay triangulation - //CGAL_triangulation_precondition(is_flipable(f,i)); - Face_handle g = f->neighbor(i); - _tds.flip( &(*f), i); - int ig=g->index(f->vertex(i)); - // set constraints to new triangles - Face_handle nfi=f->neighbor(i); - Face_handle nfj=f->neighbor(cw(i)); - Face_handle ngi=g->neighbor(ig); - Face_handle ngj=g->neighbor(ccw(ig)); - f->set_constraint(ccw(i),false); - g->set_constraint(cw(ig),false); - if (nfi->is_constrained(f->mirror_index(i))) - f->set_constraint(i,true); - else f->set_constraint(i,false); - - if (nfj->is_constrained(f->mirror_index(cw(i)))) - f->set_constraint(cw(i),true); - else f->set_constraint(cw(i),false); - - if (ngi->is_constrained(g->mirror_index(ig))) - g->set_constraint(ig,true); - else g->set_constraint(ig,false); - - if (ngj->is_constrained(g->mirror_index(ccw(ig)))) - g->set_constraint(ccw(ig),true); - else g->set_constraint(ccw(ig),false); -} - -template < class Gt, class Tds, class Itag > -void -Constrained_Delaunay_triangulation_2:: -flip_around(Vertex_handle va) - // makes the triangles incident to vertex va Delaunay using flips -{ - if (dimension() <= 1) return; - Face_handle f=va->face(); - Face_handle next; - Face_handle start(f); - int i; - do { - i = f->index(va); // FRAGILE : DIM 1 - next = f->neighbor(ccw(i)); // turns ccw around a - propagating_flip(f,i); - f=next; - } while(next != start); -} - -template < class Gt, class Tds, class Itag > -void -Constrained_Delaunay_triangulation_2:: -flip_around(List_vertices& new_vertices) -{ - typename List_vertices::iterator itv=new_vertices.begin(); - for( ; itv != new_vertices.end(); itv++) { - flip_around(*itv); - } - return; -} - - -template < class Gt, class Tds, class Itag > -void -Constrained_Delaunay_triangulation_2:: -propagating_flip(Face_handle& f,int i) - // similar to the corresponding function in Delaunay_triangulation_2.h -{ - if (!is_flipable(f,i)) return; - Face_handle ni = f->neighbor(i); - flip(f, i); // flip for constrained triangulations - propagating_flip(f,i); - i = ni->index(f->vertex(i)); - propagating_flip(ni,i); -} - - -template < class Gt, class Tds, class Itag > -void -Constrained_Delaunay_triangulation_2:: -propagating_flip(List_edges & edges) - // makes the triangulation Delaunay by flipping - // List edges contains an initial list of edges to be flipped - // Precondition : the output triangulation is Delaunay if the list - // edges contains all edges of the input triangulation that need to be - // flipped (plus possibly others) -{ - int i, ii, indf, indn; - Face_handle ni, f,ff; - Edge ei,eni; - typename Ctr::Edge_set edge_set; - typename Ctr::Less_edge less_edge; - Edge e[4]; - typename List_edges::iterator itedge=edges.begin(); - - // initialization of the set of edges to be flip - while (itedge != edges.end()) { - f=(*itedge).first; - i=(*itedge).second; - if (is_flipable(f,i)) { - eni=Edge(f->neighbor(i),f->mirror_index(i)); - if (less_edge(*itedge,eni)) edge_set.insert(*itedge); - else edge_set.insert(eni); - } - ++itedge; - } - - // flip edges and updates the set of edges to be flipped - while (!(edge_set.empty())) { - f=(*(edge_set.begin())).first; - indf=(*(edge_set.begin())).second; - - // erase from edge_set the 4 edges of the wing of the edge to be - // flipped (edge_set.begin) , i.e. the edges of the faces f and - // f->neighbor(indf) that are distinct from the edge to be flipped - - ni = f->neighbor(indf); - indn=f->mirror_index(indf); - ei= Edge(f,indf); - edge_set.erase(ei); - e[0]= Edge(f,cw(indf)); - e[1]= Edge(f,ccw(indf)); - e[2]= Edge(ni,cw(indn)); - e[3]= Edge(ni,ccw(indn)); - - for(i=0;i<4;i++) { - ff=e[i].first; - ii=e[i].second; - if (is_flipable(ff,ii)) { - eni=Edge(ff->neighbor(ii),ff->mirror_index(ii)); - if (less_edge(e[i],eni)) { - edge_set.erase(e[i]);} - else { - edge_set.erase(eni);} - } - } - - // here is the flip - flip(f, indf); - - //insert in edge_set the 4 edges of the wing of the edge that - //have been flipped - e[0]= Edge(f,indf); - e[1]= Edge(f,cw(indf)); - e[2]= Edge(ni,indn); - e[3]= Edge(ni,cw(indn)); - - for(i=0;i<4;i++) { - ff=e[i].first; - ii=e[i].second; - if (is_flipable(ff,ii)) { - eni=Edge(ff->neighbor(ii),ff->mirror_index(ii)); - if (less_edge(e[i],eni)) { - edge_set.insert(e[i]);} - else { - edge_set.insert(eni);} - } - } - } -} - -template < class Gt, class Tds, class Itag > -inline bool -Constrained_Delaunay_triangulation_2:: -test_conflict(const Point& p, Face_handle fh) const - // true if point P lies inside the circle circumscribing face fh -{ - return ( side_of_oriented_circle(fh,p) == ON_POSITIVE_SIDE ); -} - -template < class Gt, class Tds, class Itag > -inline bool -Constrained_Delaunay_triangulation_2:: -test_conflict(Face_handle fh, const Point& p) const - // true if point P lies inside the circle circumscribing face fh -{ - return test_conflict(p,fh); -} - -template < class Gt, class Tds, class Itag > -void -Constrained_Delaunay_triangulation_2:: -find_conflicts(const Point& p, std::list& le, Face_handle hint) const -{ - // sets in le the counterclocwise list of the edges of the boundary of the - // union of the faces in conflict with p - // an edge is represented by the incident face that is not in conflict with p - get_boundary_of_conflicts(p, std::back_inserter(le), hint); -} - -template < class Gt, class Tds, class Itag > -inline -typename Constrained_Delaunay_triangulation_2::Vertex_handle -Constrained_Delaunay_triangulation_2:: -virtual_insert(const Point & a, Face_handle start) - // virtual version of the insertion -{ - return insert(a,start); -} - -template < class Gt, class Tds, class Itag > -inline -typename Constrained_Delaunay_triangulation_2::Vertex_handle -Constrained_Delaunay_triangulation_2:: -virtual_insert(const Point& a, - Locate_type lt, - Face_handle loc, - int li ) -// virtual version of insert -{ - return insert(a,lt,loc,li); -} - -template < class Gt, class Tds, class Itag > -inline -typename Constrained_Delaunay_triangulation_2::Vertex_handle -Constrained_Delaunay_triangulation_2:: -insert(const Point & a, Face_handle start) - // inserts a in the triangulation -// constrained edges are updated -// Delaunay property is restored -{ - Vertex_handle va= Ctr::insert(a, start); - flip_around(va); - return va; -} - -template < class Gt, class Tds, class Itag > -typename Constrained_Delaunay_triangulation_2::Vertex_handle -Constrained_Delaunay_triangulation_2:: -insert(const Point& a, Locate_type lt, Face_handle loc, int li) -// insert a point p, whose localisation is known (lt, f, i) -// constrained edges are updated -// Delaunay property is restored -{ - Vertex_handle va= Ctr::insert(a,lt,loc,li); - flip_around(va); - return va; -} - -template < class Gt, class Tds, class Itag > -inline -typename Constrained_Delaunay_triangulation_2::Vertex_handle -Constrained_Delaunay_triangulation_2:: -push_back(const Point &p) -{ - return insert(p); -} - -template < class Gt, class Tds, class Itag > -void -Constrained_Delaunay_triangulation_2:: -triangulate_hole(List_faces& intersected_faces, - List_edges& conflict_boundary_ab, - List_edges& conflict_boundary_ba) -{ - List_edges new_edges; - Ctr::triangulate_hole(intersected_faces, - conflict_boundary_ab, - conflict_boundary_ba, - new_edges); - propagating_flip(new_edges); - return; -} - - -template < class Gt, class Tds, class Itag > -inline void -Constrained_Delaunay_triangulation_2:: -remove(Vertex_handle v) - // remove a vertex and updates the constrained edges of the new faces - // precondition : there is no incident constraints -{ - CGAL_triangulation_precondition( v != NULL ); - CGAL_triangulation_precondition( ! is_infinite(v)); - CGAL_triangulation_precondition( ! are_there_incident_constraints(v)); - if (dimension() <= 1) Ctr::remove(v); - else remove_2D(v); - return; -} - -// template < class Gt, class Tds, class Itag > -// typename -// Constrained_Delaunay_triangulation_2::Vertex_handle -// Constrained_Delaunay_triangulation_2:: -// special_insert_in_edge(const Point & a, Face_handle f, int i) -// // insert point p in edge(f,i) -// // bypass the precondition for point a to be in edge(f,i) -// // update constrained status -// // this member fonction is not robust with exact predicates -// // and approximate construction. Should be removed -// { -// Vertex_handle vh=Ctr::special_insert_in_edge(a,f,i); -// flip_around(vh); -// return vh; -// } - -template < class Gt, class Tds, class Itag > -void -Constrained_Delaunay_triangulation_2:: -remove_2D(Vertex_handle v) -{ - if (test_dim_down(v)) { _tds.remove_dim_down(&(*v)); } - else { - std::list hole; - make_hole(v, hole); - std::list shell=hole; //because hole will be emptied by fill_hole - fill_hole_delaunay(hole); - update_constraints(shell); - delete_vertex(v); - } - return; -} - -template < class Gt, class Tds, class Itag > -void -Constrained_Delaunay_triangulation_2:: -remove_incident_constraints(Vertex_handle v) -{ - List_edges iconstraints; - if (are_there_incident_constraints(v, - std::back_inserter(iconstraints))) { - Ctr::remove_incident_constraints(v); - if (dimension()==2) propagating_flip(iconstraints); - } - return; -} - -template < class Gt, class Tds, class Itag > -void -Constrained_Delaunay_triangulation_2:: -remove_constraint(Face_handle f, int i) -{ - Ctr::remove_constraint(f,i); - if(dimension() == 2) { - List_edges le; - le.push_back(Edge(f,i)); - propagating_flip(le); - } - return; -} - - -template < class Gt, class Tds, class Itag > -bool -Constrained_Delaunay_triangulation_2:: -is_valid(bool verbose, int level) const -{ - bool result = Ctr::is_valid(verbose, level); - CGAL_triangulation_assertion( result ); - - Finite_faces_iterator fit= finite_faces_begin(); - for (; fit != finite_faces_end(); fit++) { - for(int i=0;i<3;i++) { - result = result && !is_flipable(fit,i); - CGAL_triangulation_assertion( result ); - } - } - return result; -} - - - -CGAL_END_NAMESPACE -#endif // CGAL_CONSTRAINED_DELAUNAY_TRIANGULATION_2_H diff --git a/Packages/Mesh_2/include/CGAL/Filtred_circulator.h b/Packages/Mesh_2/include/CGAL/Filter_circulator.h similarity index 100% rename from Packages/Mesh_2/include/CGAL/Filtred_circulator.h rename to Packages/Mesh_2/include/CGAL/Filter_circulator.h diff --git a/Packages/Mesh_2/include/CGAL/Filtred_container.h b/Packages/Mesh_2/include/CGAL/Filtered_container.h similarity index 87% rename from Packages/Mesh_2/include/CGAL/Filtred_container.h rename to Packages/Mesh_2/include/CGAL/Filtered_container.h index 38c0486e42a..5cc502de482 100644 --- a/Packages/Mesh_2/include/CGAL/Filtred_container.h +++ b/Packages/Mesh_2/include/CGAL/Filtered_container.h @@ -4,13 +4,13 @@ CGAL_BEGIN_NAMESPACE template -class Filtred_container +class Filtered_container { Cont cont; Pred test; public: - Filtred_container(Pred p=Pred()) : cont(), test(p) {}; - Filtred_container(Cont& c, Pred p=Pred()) : cont(c), test(p) {}; + Filtered_container(Pred p=Pred()) : cont(), test(p) {}; + Filtered_container(Cont& c, Pred p=Pred()) : cont(c), test(p) {}; typedef typename Cont::reference reference; typedef typename Cont::iterator iterator; diff --git a/Packages/Mesh_2/include/CGAL/Mesh_2.h b/Packages/Mesh_2/include/CGAL/Mesh_2.h index c8ab023921b..e3560c40b2c 100644 --- a/Packages/Mesh_2/include/CGAL/Mesh_2.h +++ b/Packages/Mesh_2/include/CGAL/Mesh_2.h @@ -1,383 +1,230 @@ #ifndef CGAL_MESH_2_H #define CGAL_MESH_2_H -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include +#include #include -#include -// to skip comments and EOF in the function read_poly - -#ifdef CGAL_MESH_2_USE_TIMERS -#include -#endif +#include CGAL_BEGIN_NAMESPACE -// Tr is a Delaunay constrained triangulation (with intersections or not) -template -class Mesh_2: public Tr +class Gabriel_conform_inside_policy_2 { public: - typedef Tr Triangulation; - typedef Mesh_2 Self; - - typedef typename Tr::Geom_traits Geom_traits; - typedef typename Tr::Triangulation_data_structure Tds; + template + struct Is_acceptable_in_a_cluster + { + bool operator()(const CTr&, + const typename CTr::Face_handle& fh) const + { + return fh->is_marked(); + } + }; + + template + Is_acceptable_in_a_cluster is_acceptable_in_a_cluster_object() const + { + return Is_acceptable_in_a_cluster(); + } + + template + struct Is_locally_conform + { + private: + typedef typename CTr::Geom_traits Geom_traits; + typedef typename Geom_traits::Angle_2 Angle_2; + public: + typedef typename CTr::Point Point; + + bool operator()(const CTr& ct, + const typename CTr::Face_handle& fh, + const int i) const + { + const Angle_2 angle = ct.geom_traits().angle_2_object(); + + const Point& a = fh->vertex(ct.cw(i))->point(); + const Point& b = fh->vertex(ct.ccw(i))->point(); + const Point& c = fh->vertex(i)->point(); + const Point& d = fh->mirror_vertex(i)->point(); + + return( (!fh->is_marked() || + angle(a, c, b) != OBTUSE) && + + (!fh->neighbor(i)->is_marked() || + angle(a, d, b) != OBTUSE) + ); + } + + bool operator()(const CTr& ct, + const typename CTr::Face_handle& fh, + const int i, + const Point& p) const + { + const Angle_2 angle = ct.geom_traits().angle_2_object(); + const Point& a = fh->vertex(ct.cw(i))->point(); + const Point& b = fh->vertex(ct.ccw(i))->point(); + + return( !fh->is_marked() || angle(a, b, p) != OBTUSE ); + } + + }; + + template + Is_locally_conform is_locally_conform_object() const + { + return Is_locally_conform(); + } +}; + +/** + Tr is a Delaunay constrained triangulation (with intersections or not) +*/ +template +class Mesh_2: public Conform_triangulation_2 +{ +public: + // --- public typedefs --- + typedef Conform_triangulation_2 Conform; + + typedef Conform Base; + typedef Mesh_2 Self; + + // -- types inherited from the templated base class -- + typedef typename Base::Geom_traits Geom_traits; typedef typename Geom_traits::FT FT; typedef FT Squared_length; - typedef typename Tr::Vertex Vertex; - typedef typename Tr::Edge Edge; - typedef typename Tr::Finite_edges_iterator Finite_edges_iterator; - typedef typename Tr::Edge_circulator Edge_circulator; - typedef typename Tr::Face_handle Face_handle; - typedef typename Tr::Finite_faces_iterator Finite_faces_iterator; - typedef typename Tr::All_faces_iterator All_faces_iterator; - typedef typename Tr::Face_circulator Face_circulator; - typedef typename Tr::Vertex_handle Vertex_handle; - typedef typename Tr::Vertex_circulator Vertex_circulator; - typedef typename Tr::Finite_vertices_iterator Finite_vertices_iterator; + typedef typename Base::Vertex_handle Vertex_handle; + typedef typename Base::Face_handle Face_handle; + + typedef typename Base::Face_circulator Face_circulator; + typedef typename Base::Finite_faces_iterator Finite_faces_iterator; + typedef typename Base::All_faces_iterator All_faces_iterator; + typedef typename Base::Point Point; - typedef typename Tr::Locate_type Locate_type; + // -- types needed to access member datas + typedef std::list Seeds; + typedef typename Seeds::iterator Seeds_iterator; + typedef typename Seeds::const_iterator Seeds_const_iterator; + + // --- CONSTRUCTORS --- - typedef typename Tr::Point Point; - - typedef typename Tr::List_constraints List_constraints; - -public: - // CONSTRUCTORS - - // default constructor, to be used with IO functions - Mesh_2(const Geom_traits& gt = Geom_traits()); - - // constructor from a Triangulation& t - explicit - Mesh_2(Triangulation& t, const Geom_traits& gt = Geom_traits(), - bool dont_refine = false); - - // the same with declaration of seeds. - // See the explanation of functions refine(...) for the - // signification of arguments Seed_it begin,end and bool mark. - // The parameter dont_refine should not be documented because it - // allow to construct a "mesh" that is not meshed. - template - explicit - Mesh_2(Triangulation& t, Seed_it begin, Seed_it end, - const Geom_traits& gt = Geom_traits(), bool mark = false, - bool dont_refine = false) - { - swap(t); - if(!dont_refine) - refine(begin, end, mark); - }; - - // ASSIGNEMENT - // TODO! - - // ACCESS FUNCTIONS - unsigned int number_of_constrained_edges() const; + explicit Mesh_2(const Geom_traits& gt = Geom_traits()): + Base(gt), initialized(false) {} + // --- ACCESS FUNCTION --- bool is_bad(const Face_handle fh) const; double squared_minimum_sine(const Face_handle fh) const; double squared_minimum_sine(const Vertex_handle& va, const Vertex_handle& vb, const Vertex_handle& vc) const; - // IO - // write and read the constrained edges in the format: - // number_of_edges - // segment1 - // segment2 - // ... - void write(std::ostream &f) const; - void read(std::istream &f, bool dont_refine = false); + Seeds_const_iterator seeds_begin() const; + Seeds_const_iterator seeds_end() const; - // write and read a mesh in the Triangle .poly format - // (see http://www-2.cs.cmu.edu/~quake/triangle.poly.html) - void write_poly(std::ostream &f) const; - void read_poly(std::istream &f, bool dont_refine = false); - - // HELPING FUNCTION + // --- HELPING FUNCTION --- void clear(); - // MESHING FUNCTIONS + // --- MARKING FUNCTIONS --- + /** The value type of InputIterator should be Point, and represents + seeds. Connected components of seeds are marked with the value of + "mark". Other components are marked with !mark. The connected + component of infinite faces is always marked with false. + */ + template + void set_seeds(InputIterator b, InputIterator e, + const bool mark = false, + const bool do_it_now = false) + { + seeds.clear(); + copy(b, e, std::back_inserter(seeds)); + _mark=mark; + if(do_it_now) mark_facets(); + } - // Perform meshing. All faces but those connected to the infinite - // faces are marked. + void clear_seeds() + { + seeds.clear(); + } + + /** Procedure that forces facets to be marked immediatelly. */ + void mark_facets(); + + // --- MESHING FUNCTIONS --- + + // Perform meshing. void refine(); - // Perform meshing. Seed_it is an iterator of points, representing - // seeds. Connected components of seeds are marked with the value of - // "mark". Other components are marked with !mark. The connected - // component of infinite faces is always marked with false. - template - void refine(Seed_it begin, Seed_it end, bool mark = false) - { - init(begin, end, mark); + // --- REMESHING FUNCTIONS --- -#ifdef CGAL_MESH_2_USE_TIMERS - timer.reset(); - timer.start(); -#endif - while(! (edges_to_be_conformed.empty() && bad_faces.empty()) ) - { - conform(); - if ( !bad_faces.empty() ) - process_one_face(); - } -#ifdef CGAL_MESH_2_USE_TIMERS - timer.stop(); - std::cout << "Refine mesh time: " << timer.time() << std::endl; -#endif - } - - // REMESHING FUNCTIONS - - // Set the geom_traits and recalculate the list of bad faces + // Set the geom_traits nut DO NOT recalculate the list of bad faces (must + // call set_bad_faces of calculate_bad_faces bellow) void set_geom_traits(const Geom_traits& gt); + // re-calculate the list of bad faces + void calculate_bad_faces(); + // Set the geom_traits and add the sequence [begin, end[ to the list // of bad faces. // Fh_it is a iterator of Face_Handle. // Use this overriden function if the list of bad faces can be // computed easily without testing all faces. template - void set_geom_traits(const Geom_traits& gt, - Fh_it begin, Fh_it end) + void set_bad_faces(Fh_it begin, Fh_it end) { - _gt = gt; + bad_faces.clear(); for(Fh_it pfit=begin; pfit!=end; ++pfit) push_in_bad_faces(*pfit); } - // STEP BY STEP MESHING + // --- STEP BY STEP FUNCTIONS --- - // two deprecated constructors - - Mesh_2(List_constraints& lc, const Geom_traits& gt = Geom_traits()); - - template - Mesh_2(InputIterator first, InputIterator last, - const Geom_traits& gt = Geom_traits()) - : Tr(gt), is_really_a_constrained_edge(*this), - edges_to_be_conformed(is_really_a_constrained_edge) - { - while(first != last){ - insert_constraint((*first).first, (*first).second); - ++first; - } - CGAL_triangulation_postcondition(is_valid()); - } - - // init(...): Initialize the data structures - // (The call of one of the following init function is REQUIRED - // before any step by step operation). - // See the explanations of functions refine(...) for the - // signification of arguments + /** + init(): Initialize the data structures + (The call of this function is REQUIRED before any step by step + operation). + */ void init(); - template - void init(Seed_it begin, Seed_it end, bool mark = false) - { - cluster_map.clear(); - edges_to_be_conformed.clear(); - bad_faces.clear(); - - mark_facets(begin, end, mark); - - create_clusters(); - fill_edge_queue(); - fill_facet_map(); - } - - // mark_facets(...): procedure called to mark facets, same arguments - // as init(...) above - template - void mark_facets(Seed_it begin, Seed_it end, bool mark = false) - { - if (dimension()<2) return; - if(begin!=end) - { - for(All_faces_iterator it=all_faces_begin(); - it!=all_faces_end(); - ++it) - it->set_marked(!mark); - - for(Seed_it sit=begin; sit!=end; ++sit) - { - Face_handle fh=locate(*sit); - if(fh!=NULL) - propagate_marks(fh, mark); - } - } - else - mark_convex_hull(); - propagate_marks(infinite_face(), false); - }; - - // Conform edges and doesn't refine faces - // Need init(...) see bellow - void conform(); - - // Execute on step of the algorithm. - // Need init(...) see bellow + /** Execute one step of the algorithm. + Needs init() see above */ bool refine_step(); private: - // PRIVATE TYPES - typedef std::pair - Constrained_edge; - typedef std::list List_of_edges; - typedef std::list List_of_face_handles; + // --- PRIVATE TYPES --- typedef CGAL::Triple Threevertices; - // traits type - typedef typename Geom_traits::Vector_2 Vector_2; - typedef typename Geom_traits::Construct_translated_point_2 - Construct_translated_point_2; - typedef typename Geom_traits::Compute_squared_distance_2 - Compute_squared_distance_2; - typedef typename Geom_traits::Angle_2 Angle_2; - typedef typename Geom_traits::Construct_vector_2 - Construct_vector_2; - typedef typename Geom_traits::Construct_scaled_vector_2 - Construct_scaled_vector_2; - typedef typename Geom_traits::Construct_midpoint_2 - Construct_midpoint_2; - typedef typename Geom_traits::Orientation_2 Orientation_2; + Vertex_handle> Threevertices; + + typedef std::list List_of_edges; + typedef std::list List_of_face_handles; + typedef typename Base::Cluster Cluster; + + // -- traits type -- + typedef typename Geom_traits::Is_bad Is_bad; typedef typename Geom_traits::Compute_squared_minimum_sine_2 Compute_squared_minimum_sine_2; - typedef typename Geom_traits::Is_bad Is_bad; - - // Cluster register several informations about clusters. - // A cluster is a is a set of vertices v_i incident to one vertice - // v_0, so that angles between segments [v_0, v_i] is less than 60°. - struct Cluster { - bool reduced ; // Is the cluster reduced - - // smallest_angle gives the two vertices defining the - // smallest angle in the cluster - std::pair smallest_angle; - - FT rmin; // WARNING: rmin has no meaning if reduced=false!!! - Squared_length minimum_squared_length; - - // The following map tells what vertices are in the cluster and if - // the corresponding segment has been splitted once. - typedef std::map Vertices_map; - Vertices_map vertices; - - bool is_reduced() const { - return reduced; - } - - bool is_reduced(const Vertex_handle v) { - return vertices[v]; - } - }; - - // Is_edge_constrained: function object class. - // Take a Mesh_2 object m and a Vertex_handle v in constructor. - // Is_edge_constrained(Vertex_handle v2) tells if [v,v2] is a - // constrained edge in m. - class Is_edge_constrained { - Self* _m; - Vertex_handle _v; - public: - Is_edge_constrained(Self* m, Vertex_handle v) - : _m(m), _v(v) {} - - Is_edge_constrained(const Is_edge_constrained& other) - : _m(other._m), _v(other._v) {} - - Is_edge_constrained& - operator=(const Is_edge_constrained& other) - { - _m = other._m; - _v = other._v; - return *this; - } - - bool operator()(/*const*/ Vertex& v2) const - { - if(_m->is_infinite(v2.handle())) - return false; - Face_handle fh; - int i; - _m->is_edge(_v,v2.handle(),fh,i); - return fh->is_constrained(i); - } - }; - - class Is_really_a_constrained_edge { - const Self& _m; - public: - explicit Is_really_a_constrained_edge(const Self& m) : _m(m) {}; - bool operator()(const Constrained_edge& ce) const - { - Face_handle fh; - int i; - return _m.is_edge(ce.first, ce.second, fh,i) && - fh->is_constrained(i); - } - }; - - // typedefs for private members types - typedef Filtred_circulator - Constrained_vertex_circulator; + typedef typename Geom_traits::Compute_squared_distance_2 + Compute_squared_distance_2; + // -- typedefs for private members types -- typedef CGAL::Double_map Bad_faces; - typedef std::list List_of_constraints; - typedef CGAL::Filtred_container - Constrained_edges_queue; - - typedef std::multimap Cluster_map; - private: - // PRIVATE MEMBER DATAS + // --- PRIVATE MEMBER DATAS --- // bad_faces: list of bad finite faces // warning: some faces could be recycled during insertion in the // triangulation, that's why we need to be able to remoce faces // from the map. Bad_faces bad_faces; + bool initialized; - // edges_to_be_conformed: list of encroached constrained edges - // warning: some edges could be destroyed, use the same wrapper - // is_really_a_constrained_edge: tester to filter edges_to_be_conformed - const Is_really_a_constrained_edge is_really_a_constrained_edge; - Constrained_edges_queue edges_to_be_conformed; - - // Cluster_map: multimap Vertex_handle -> Cluster - // each vertex can have several clusters - Cluster_map cluster_map; - -#ifdef CGAL_MESH_2_USE_TIMERS - // Timer to bench - Timer timer; -#endif - -public: - // PUBLIC DATA MEMBERS - typedef std::list Seeds; - Seeds seeds; // TODO, WARNING: have to think about this variable and - // about the default marker. + Seeds seeds; + bool _mark; private: - // PRIVATE MEMBER FUNCTIONS + // --- PRIVATE MEMBER FUNCTIONS --- // -- auxiliary functions to set markers -- @@ -388,43 +235,7 @@ private: // propagate the mark m recursivly void propagate_marks(Face_handle f, bool m); - // -- auxiliary functions to handle clusters -- - // for all vertices, call create_clusters_of_vertex - void create_clusters(); - - // compute clusters of the vertex v, using the auxiliary function - // construct_cluster - void create_clusters_of_vertex(Vertex_handle v); - - // add the sequence [begin, end] to the cluster c and add it to the - // clusters of the vertex v - void construct_cluster(Vertex_handle v, - Constrained_vertex_circulator begin, - const Constrained_vertex_circulator& end, - Cluster c = Cluster()); - - // update the cluster of [va,vb], putting vm instead of vb. If - // reduction=false, the edge [va,vm] is not set reduced. - void update_cluster(Cluster& c, Vertex_handle va, Vertex_handle vb, - Vertex_handle vm, bool reduction = true); - - // get_cluster returns the cluster of [va,vb] in c and return true - // if it is in a cluster. If erase=true, the cluster is remove from - // the cluster map - bool get_cluster(Vertex_handle va, Vertex_handle vb, Cluster &c, - bool erase = false); - - // -- functions that maintain the queue of encroached edges -- - - // scan all constrained edges and put them in the queue if they are - // encroached - void fill_edge_queue(); - - // update the queue with edges incident to vm - void update_edges_to_be_conformed(Vertex_handle va, - Vertex_handle vb, - Vertex_handle vm); // -- functions that maintain the map of bad faces @@ -443,109 +254,40 @@ private: // -- inlined functions that compose the refinement process -- - // take one edge in the queue and call refine_edge - void process_one_edge(); - // take one face in the queue and call refine_face void process_one_face(); - // handle the encroached edge, call cut_cluster_edge+update_cluster - // or insert_middle - void refine_edge(Vertex_handle va, Vertex_handle vb); - // handle one face; call split_face or put in the edges_to_be_conformed the // list of edges that would be encroached by the circum_center of f // This function uses Shewchuk's terminator criteria. void refine_face(Face_handle f); - // -- auxiliary functions that return a boolean -- - // tell if [va,vb] is encroached, by looking for the two neighbors - // This function takes of markers. - bool is_encroached(const Vertex_handle va, - const Vertex_handle vb) const; - - // tell if [va,vb] is encroached by p - bool is_encroached(const Vertex_handle va, - const Vertex_handle vb, - const Point& p) const; - - // tell if the angle is less than 60° - // Uses squared_cosine_of_angle_times_4 and used by - // create_clusters_of_vertex - bool is_small_angle(const Point& pleft, - const Point& pmiddle, - const Point& pright) const; - - // -- auxiliary functions that are called to split an edge or -- - // -- a face -- - - // cut [va,vb] knowing that it is in the cluster c - void cut_cluster_edge(Vertex_handle va, Vertex_handle vb, - Cluster& c); - - // insert the midpoint of the edge (f,i) - Vertex_handle insert_middle(Face_handle f, int i); // -- functions that really insert points -- // split the face f by inserting its circum center circum_center void split_face(const Face_handle& f, const Point& circum_center); - // insert the point p in the edge (fh,edge_index) - Vertex_handle insert_in_the_edge(Face_handle fh, int edge_index, + // overrideen functions that inserts the point p in the edge + // (fh,edge_index) + Vertex_handle insert_in_the_edge(Face_handle fh, const int edge_index, const Point& p); - - // -- helping computing functions -- - // return the squared cosine of the angle - // times 4 - FT squared_cosine_of_angle_times_4(const Point& pleft, - const Point& pmiddle, - const Point& pright) const; + // -- helping computing functions -- // return the squared length of the triangle corresponding to the // face f Squared_length shortest_edge_squared_length(Face_handle f); + // -- debugging functions -- + bool is_bad_faces_valid(); + }; // end of Mesh_2 -// CONSTRUCTORS - -template -Mesh_2:: -Mesh_2(const Geom_traits& gt) - : Tr(gt), is_really_a_constrained_edge(*this), - edges_to_be_conformed(is_really_a_constrained_edge) -{}; - -template -Mesh_2:: -Mesh_2(Triangulation& t, const Geom_traits& gt, bool dont_refine) - : Tr(gt), is_really_a_constrained_edge(*this), - edges_to_be_conformed(is_really_a_constrained_edge) -{ - swap(t); - if(!dont_refine) - refine(); -}; - -// ACCESS FUNCTIONS - -template -unsigned int Mesh_2:: -number_of_constrained_edges() const -{ - int nedges = 0; - for(Finite_edges_iterator eit = finite_edges_begin(); - eit != finite_edges_end(); - ++eit) - if((*eit).first->is_constrained((*eit).second)) - ++nedges; - return nedges; -} +// --- ACCESS FUNCTIONS --- // ????????????? // ->traits @@ -572,6 +314,7 @@ squared_minimum_sine(const Vertex_handle& va, const Vertex_handle& vb, { Compute_squared_minimum_sine_2 squared_sine = geom_traits().compute_squared_minimum_sine_2_object(); + return squared_sine(va->point(), vb->point(), vc->point()); } @@ -584,161 +327,68 @@ squared_minimum_sine(const Face_handle fh) const & va = fh->vertex(0), & vb = fh->vertex(1), & vc = fh->vertex(2); + return squared_minimum_sine(va, vb, vc); } -// IO -//the function that writes a file template -void Mesh_2:: -write(std::ostream &f) const +inline +typename Mesh_2::Seeds_const_iterator +Mesh_2:: +seeds_begin() const { - f << number_of_constrained_edges() << std::endl; - for(Finite_edges_iterator eit = finite_edges_begin(); - eit!=finite_edges_end(); - ++eit) - if((*eit).first->is_constrained((*eit).second)) - { - f << (*eit).first->vertex(cw((*eit).second))->point() << " " - << (*eit).first->vertex(ccw((*eit).second))->point() < -void Mesh_2:: -read(std::istream &f, bool dont_refine) +inline +typename Mesh_2::Seeds_const_iterator +Mesh_2:: +seeds_end() const { - int nedges = 0; - clear(); - f>>nedges; - for(int n = 0; n> p1 >> p2; - insert_constraint(p1, p2); - } - if(!dont_refine) - refine(); + return seeds.end(); } -//the function that write a Shewchuk Triangle .poly file -template -void Mesh_2:: -write_poly(std::ostream &f) const -{ - std::map index; - // write vertices - f << "# Shewchuk Triangle .poly file, produced by the CGAL::Mesh_2 package" - << std::endl - << "# Neither attributes nor boundary markers are used." << std::endl - << number_of_vertices() << " " << 2 << " " - << 0 << " " << 0 << std::endl; - f << std::endl; - - unsigned int vertices_counter = 0; - for(Finite_vertices_iterator vit = finite_vertices_begin(); - vit != finite_vertices_end(); - ++vit) - { - f << ++vertices_counter << " " << vit->point() << std::endl; - index[vit] = vertices_counter; - } - - f << std::endl; - - // write constrained edges - - f << number_of_constrained_edges() << " " << 0 << std::endl; - unsigned int edges_counter = 0; - for(Finite_edges_iterator eit = finite_edges_begin(); - eit != finite_edges_end(); - ++eit) - if((*eit).first->is_constrained((*eit).second)) - f << ++edges_counter << " " - << index[(*eit).first->vertex(cw((*eit).second))] << " " - << index[(*eit).first->vertex(ccw((*eit).second))] - << std::endl; - - f << std::endl; - - // write seeds, assuming that the seeds unmarks faces - unsigned int seeds_counter = 0; - f << seeds.size() << std::endl; - for(typename Seeds::const_iterator sit = seeds.begin(); - sit!=seeds.end(); ++sit) - f << ++seeds_counter << " " << *sit << std::endl; -} - -//the function that reads a Shewchuk Triangle .poly file -template -void Mesh_2:: -read_poly(std::istream &f, bool dont_refine) -{ - clear(); - - unsigned int number_of_points; - skip_comment_OFF(f); - f >> number_of_points; - skip_until_EOL(f); - skip_comment_OFF(f); - - // read vertices - std::vector vertices(number_of_points); - for(unsigned int i = 0; i < number_of_points; ++i) - { - unsigned int j; - Point p; - f >> j >> p; - skip_until_EOL(f); skip_comment_OFF(f); - vertices[--j] = insert(p); - } - - // read segments - unsigned int number_of_segments; - f >> number_of_segments; - skip_until_EOL(f); skip_comment_OFF(f); - for(unsigned int k = 0; k < number_of_segments; ++k) - { - unsigned int l, v1, v2; - f >> l >> v1 >> v2; - skip_until_EOL(f); skip_comment_OFF(f); - insert_constraint(vertices[--v1], vertices[--v2]); - } - - // read holes - unsigned int number_of_holes; - f >> number_of_holes; - for(unsigned int m = 0; m < number_of_holes; ++m) - { - unsigned int n; - Point p; - f >> n >> p; - skip_until_EOL(f); skip_comment_OFF(f); - seeds.push_back(p); - } - if(dont_refine) - init(seeds.begin(), seeds.end(), false); - else - refine(seeds.begin(), seeds.end(), false); -} - -// HELPING FUNCTIONS +// --- HELPING FUNCTIONS --- template void Mesh_2:: clear() { - cluster_map.clear(); - edges_to_be_conformed.clear(); bad_faces.clear(); seeds.clear(); - Triangulation::clear(); + Base::clear(); } -// MESHING FUNCTIONS +// --- MARKING FUNCTIONS --- +template +void Mesh_2:: +mark_facets() +{ + if (Base::dimension()<2) return; + if( seeds_begin() != seeds_end() ) + { + for(All_faces_iterator it=all_faces_begin(); + it!=all_faces_end(); + ++it) + it->set_marked(!_mark); + + for(Seeds_const_iterator sit=seeds_begin(); sit!=seeds_end(); ++sit) + { + Face_handle fh=locate(*sit); + if(fh!=NULL) + propagate_marks(fh, _mark); + } + } + else + mark_convex_hull(); + propagate_marks(infinite_face(), false); +}; + +// --- MESHING FUNCTIONS --- //the mesh refine function template @@ -746,55 +396,47 @@ inline void Mesh_2:: refine() { - std::list l; - refine(l.begin(), l.end()); -} + if(!initialized) init(); -// REMESHING FUNCTIONS - -// deprecated constructor -template -Mesh_2:: -Mesh_2(List_constraints& lc, const Geom_traits& gt) - : Tr(gt), is_really_a_constrained_edge(*this), - edges_to_be_conformed(is_really_a_constrained_edge) -{ - typename List_constraints::iterator lcit = lc.begin(); - for( ; lcit != lc.end(); ++lcit) + while(!Conform::is_conformed() || !bad_faces.empty() ) { - insert_constraint( (*lcit).first, (*lcit).second); + Conform::conform(Gabriel_conform_inside_policy_2()); + if ( !bad_faces.empty() ) + process_one_face(); } - CGAL_triangulation_postcondition(is_valid()); } +// --- REMESHING FUNCTIONS --- + template +inline void Mesh_2:: set_geom_traits(const Geom_traits& gt) { - _gt = gt; - bad_faces.clear(); + this->_gt = gt; +} + +template +inline +void Mesh_2:: +calculate_bad_faces() +{ fill_facet_map(); } -// STEP BY STEP MESHING +// --- STEP BY STEP FUNCTIONS --- + template inline void Mesh_2:: init() { - std::list l; - init(l.end(), l.end()); -} + bad_faces.clear(); + mark_facets(); // facets should be marked before the call to Base::init() -template -inline -void Mesh_2:: -conform() -{ - while( !edges_to_be_conformed.empty() ) - { - process_one_edge(); - }; + Base::init(Gabriel_conform_inside_policy_2()); // handles clusters and + // edges + fill_facet_map(); } template @@ -802,8 +444,8 @@ inline bool Mesh_2:: refine_step() { - if( !edges_to_be_conformed.empty() ) - process_one_edge(); + if( !Conform::is_conformed() ) + Conform::refine_step(Gabriel_conform_inside_policy_2()); else if ( !bad_faces.empty() ) process_one_face(); @@ -812,7 +454,7 @@ refine_step() return true; } -// PRIVATE MEMBER FUNCTIONS +// --- PRIVATE MEMBER FUNCTIONS --- template void Mesh_2:: @@ -830,8 +472,9 @@ void Mesh_2:: propagate_marks(const Face_handle fh, bool mark) { // std::queue only works with std::list on VC++6, and not with - // std::deque, wichi is the default - std::queue > face_queue; + // std::deque, which is the default + // But it should be fixed by VC++7 know. [Laurent Rineau 2003/03/24] + std::queue*/> face_queue; fh->set_marked(mark); face_queue.push(fh); while( !face_queue.empty() ) @@ -850,237 +493,6 @@ propagate_marks(const Face_handle fh, bool mark) } }; -template -void Mesh_2:: -create_clusters() -{ - for(Finite_vertices_iterator vit = finite_vertices_begin(); - vit != finite_vertices_end(); - vit++) - create_clusters_of_vertex(vit); -} - -template -void Mesh_2:: -create_clusters_of_vertex(Vertex_handle v) -{ - Is_edge_constrained test(this, v); - Constrained_vertex_circulator begin(incident_vertices(v),test); - // This circulator represents all constrained edges around the - // vertex v. An edge [v,v'] is represented by the vertex v'. - - if(begin == 0) return; // if there is only one vertex - - Constrained_vertex_circulator - current(begin), next(begin), cluster_begin(begin); - ++next; // next is always just after current. - if(current == next) return; - - bool in_a_cluster = false; - do - { - Face_handle f; - int i; - is_edge(v, next, f, i); // put in f the face on the right side - // of (v,next) - if(f->is_marked() && - is_small_angle(current->point(), v->point(), next->point())) - { - if(!in_a_cluster) - { - // at this point, current is the beginning of a cluster - in_a_cluster = true; - cluster_begin = current; - } - } - else - if(in_a_cluster) - { - // at this point, current is the end of a cluster and - // cluster_begin is its beginning - construct_cluster(v, cluster_begin, current); - in_a_cluster = false; - } - ++next; - ++current; - } while( current!=begin ); - if(in_a_cluster) - { - Cluster c; - if(get_cluster(v, begin, c, true)) - // get the cluster and erase it from the clusters map - construct_cluster(v, cluster_begin, begin, c); - else - construct_cluster(v, cluster_begin, current); - } -} - -template -void Mesh_2:: -construct_cluster(Vertex_handle v, - Constrained_vertex_circulator begin, - const Constrained_vertex_circulator& end, - Cluster c) -{ - Compute_squared_distance_2 squared_distance = - geom_traits().compute_squared_distance_2_object(); - - if(c.vertices.empty()) - { - c.reduced = false; - // c.rmin is not initialized because - // reduced=false! - c.minimum_squared_length = - squared_distance(v->point(), begin->point()); - Constrained_vertex_circulator second(begin); - ++second; - c.smallest_angle.first = begin; - c.smallest_angle.second = second; - } - - bool all_edges_in_cluster=false; // tell if all incident edges are - // in the cluster - if(begin==end) - all_edges_in_cluster=true; - - const Point& vp = v->point(); - - FT greatest_cosine = - squared_cosine_of_angle_times_4(c.smallest_angle.first->point(), - v->point(), - c.smallest_angle.second->point()); - - Constrained_vertex_circulator next(begin); - ++next; - do - { - c.vertices[begin] = false; - Squared_length l = squared_distance(vp, - begin->point()); - c.minimum_squared_length = - std::min(l,c.minimum_squared_length); - - if(all_edges_in_cluster || begin!=end) - { - FT cosine = - squared_cosine_of_angle_times_4(begin->point(), - v->point(), - next->point()); - if(cosine>greatest_cosine) - { - greatest_cosine = cosine; - c.smallest_angle.first = begin; - c.smallest_angle.second = next; - } - } - } - while(next++,begin++!=end); - cluster_map.insert(std::make_pair(v,c)); -} - -template -void Mesh_2:: -update_cluster(Cluster& c, Vertex_handle va,Vertex_handle vb, - Vertex_handle vm, bool reduction) -{ - Compute_squared_distance_2 squared_distance = - geom_traits().compute_squared_distance_2_object(); - c.vertices.erase(vb); - c.vertices[vm] = reduction; - - if(vb==c.smallest_angle.first) - c.smallest_angle.first = vm; - if(vb==c.smallest_angle.second) - c.smallest_angle.second = vm; - - FT l = squared_distance(va->point(),vm->point()); - if(lfirst)) - ++it; // TODO: use std::find and an object class - if(it==c.vertices.end()) - c.reduced = true; - } - - if(c.is_reduced()) - c.rmin = squared_distance(c.smallest_angle.first->point(), - c.smallest_angle.second->point())/FT(4); - cluster_map.insert(std::make_pair(va,c)); -} - -// # used by refine_face and cut_cluster_edge -template -bool Mesh_2:: -get_cluster(Vertex_handle va, Vertex_handle vb, Cluster &c, bool erase) -{ - typedef typename Cluster_map::iterator Iterator; - typedef std::pair Range; - Range range = cluster_map.equal_range(va); - for(Iterator it = range.first; it != range.second; it++) - { - Cluster &cl = it->second; - if(cl.vertices.find(vb)!=cl.vertices.end()) { - c = it->second; - if(erase) - cluster_map.erase(it); - return true; - } - } - return false; -} - -template -void Mesh_2:: -fill_edge_queue() -{ - for(Finite_edges_iterator ei = finite_edges_begin(); - ei != finite_edges_end(); - ++ei) - { - Vertex_handle va = (*ei).first->vertex(cw((*ei).second)); - Vertex_handle vb = (*ei).first->vertex(ccw((*ei).second)); - if((*ei).first->is_constrained((*ei).second) && - is_encroached(va, vb)) - { - edges_to_be_conformed.push_back(std::make_pair(va, vb)); - } - } -} - -//update the encroached segments list -// TODO: perhaps we should remove destroyed edges too -// TODO: rewrite this function one day -template -void Mesh_2:: -update_edges_to_be_conformed(Vertex_handle va, Vertex_handle vb, - Vertex_handle vm) -{ - Face_circulator fc = incident_faces(vm), fcbegin(fc); - - do { - for(int i = 0; i<3; i++) { - Vertex_handle v1 = fc->vertex(cw(i)); - Vertex_handle v2 = fc->vertex(ccw(i)); - if( fc->is_constrained(i) && !is_infinite(v1) && - !is_infinite(v2) && is_encroached(v1, v2) ) - edges_to_be_conformed.push_back(Constrained_edge(v1, v2)); - } - ++fc; - } while(fc != fcbegin); - - if(is_encroached(va, vm)) { - edges_to_be_conformed.push_back(Constrained_edge(va, vm)); - } - - if(is_encroached(vb, vm)) { - edges_to_be_conformed.push_back(Constrained_edge(vb, vm)); - } -} - template inline void Mesh_2:: @@ -1126,17 +538,6 @@ compute_new_bad_faces(Vertex_handle v) } while(fc!=fcbegin); } -template -inline -void Mesh_2:: -process_one_edge() -{ - Constrained_edge ce = edges_to_be_conformed.front(); - edges_to_be_conformed.pop_front(); - - refine_edge(ce.first, ce.second); -} - template inline void Mesh_2:: @@ -1147,42 +548,17 @@ process_one_face() refine_face(f); } -//this function split all the segments that are encroached -template -void Mesh_2:: -refine_edge(Vertex_handle va, Vertex_handle vb) -{ - Face_handle f; - int i; - is_edge(va, vb, f, i); // get the edge (f,i) - CGAL_assertion(f->is_constrained(i)); - - Cluster c,c2; - - if( get_cluster(va,vb,c,true) ) - if( get_cluster(vb,va,c2,true) ) - { // both ends are clusters - Vertex_handle vm = insert_middle(f,i); - update_cluster(c,va,vb,vm,false); - update_cluster(c2,vb,va,vm,false); - } - else - // va only is a cluster - cut_cluster_edge(va,vb,c); - else - if( get_cluster(vb,va,c,true) ) - // vb only is a cluster - cut_cluster_edge(vb,va,c); - else - // no cluster - insert_middle(f,i); -}; - - //split all the bad faces +//split all the bad faces template void Mesh_2:: refine_face(const Face_handle f) { + typedef typename Gabriel_conform_inside_policy_2:: + template Is_locally_conform Is_locally_conform; + + Is_locally_conform is_gabriel_conform_inside = + Gabriel_conform_inside_policy_2().is_locally_conform_object(); + const Point& pc = circumcenter(f); List_of_edges zone_of_pc_boundary; @@ -1208,7 +584,7 @@ refine_face(const Face_handle f) const Vertex_handle & va = fh->vertex(cw(i)), & vb = fh->vertex(ccw(i)); - if(fh->is_constrained(i) && is_encroached(va,vb,pc)) + if(fh->is_constrained(i) && !is_gabriel_conform_inside(*this, fh, i, pc)) { split_the_face = false; Cluster c,c2; @@ -1219,7 +595,7 @@ refine_face(const Face_handle f) (!is_cluster_at_va && !is_cluster_at_vb) ) { // two clusters or no cluster - edges_to_be_conformed.push_back(Constrained_edge(va,vb)); + add_contrained_edge_to_be_conform(va,vb); keep_the_face_bad = true; } else @@ -1238,12 +614,14 @@ refine_face(const Face_handle f) if( !c.is_reduced() || c.rmin >= shortest_edge_squared_length(f) ) { - edges_to_be_conformed.push_back(Constrained_edge(va,vb)); + add_contrained_edge_to_be_conform(va,vb); keep_the_face_bad = true; } } } - }; // after here edges_to_be_conformed contains edges encroached by pc + }; // after here edges encroached by pc are in the list of edges to + // be conformed. + const Vertex_handle & va = f->vertex(0), @@ -1252,10 +630,6 @@ refine_face(const Face_handle f) if(split_the_face) { -// int li; -// Locate_type lt; -// locate(pc,lt,li,f); -// if(lt!=OUTSIDE_CONVEX_HULL) CGAL_assertion(f->is_marked()); if(f->is_marked()) split_face(f, pc); @@ -1265,139 +639,6 @@ refine_face(const Face_handle f) push_in_bad_faces(va, vb, vc); } -template -bool Mesh_2:: -is_encroached(const Vertex_handle va, const Vertex_handle vb, - const Point& p) const -{ - Angle_2 angle = geom_traits().angle_2_object(); - - return angle(va->point(), p, vb->point())==OBTUSE; -} - -template -bool Mesh_2:: -is_encroached(const Vertex_handle va, const Vertex_handle vb) const -{ - Face_handle fh; - int i; - is_edge(va,vb,fh,i); - - const Point& candidat_1 = fh->vertex(i)->point(); - const Point& candidat_2 = fh->mirror_vertex(i)->point(); - - return ( ( fh->is_marked() && - is_encroached(va, vb, candidat_1) ) || - ( fh->neighbor(i)->is_marked() && - is_encroached(va, vb, candidat_2) ) - ); -} - -// -> traits? -template -bool Mesh_2:: -is_small_angle(const Point& pleft, - const Point& pmiddle, - const Point& pright) const -{ - Angle_2 angle = geom_traits().angle_2_object(); - Orientation_2 orient = geom_traits().orientation_2_object(); - - if( angle(pleft, pmiddle, pright)==OBTUSE ) - return false; - if( orient(pmiddle,pleft,pright)==RIGHTTURN) - return false; - - FT cos_alpha = squared_cosine_of_angle_times_4(pleft, pmiddle, - pright); - - if(cos_alpha > 1) - { - return true; //the same cluster - } - else - { - return false; //another cluster - } -} - -template -void Mesh_2:: -cut_cluster_edge(Vertex_handle va, Vertex_handle vb, Cluster& c) -{ - Construct_vector_2 vector = - geom_traits().construct_vector_2_object(); - Construct_scaled_vector_2 scaled_vector = - geom_traits().construct_scaled_vector_2_object(); - Compute_squared_distance_2 squared_distance = - geom_traits().compute_squared_distance_2_object(); - Construct_midpoint_2 midpoint = - geom_traits().construct_midpoint_2_object(); - Construct_translated_point_2 translate = - geom_traits().construct_translated_point_2_object(); - - Vertex_handle vc; - - if(c.is_reduced(vb)) - { - Face_handle fh; - int i; - is_edge(va,vb,fh,i); - vc = insert_middle(fh,i); - } - else - { - const Point - & a = va->point(), - & b = vb->point(), - & m = midpoint(a, b); - - - - Vector_2 v = vector(a,m); - v = scaled_vector(v,CGAL_NTS sqrt(c.minimum_squared_length / - squared_distance(a,b))); - - Point i = translate(a,v), i2(i); - - do { - i = translate(a,v); - v = scaled_vector(v,FT(2)); - i2 = translate(a,v); - } while(squared_distance(a,i2) <= squared_distance(a,m)); - if( squared_distance(i,m) > squared_distance(m,i2) ) - i = i2; - //here i is the best point for splitting - Face_handle fh; - int index; - is_edge(va,vb,fh,index); - - vc = insert_in_the_edge(fh, index, i); - } - update_edges_to_be_conformed(va, vb, vc); - compute_new_bad_faces(vc); - update_cluster(c, va, vb, vc); -} - -template -typename Mesh_2::Vertex_handle -Mesh_2:: -insert_middle(Face_handle f, int i) -{ - Construct_midpoint_2 - midpoint = geom_traits().construct_midpoint_2_object(); - - const Vertex_handle - & va = f->vertex(cw(i)), - & vb = f->vertex(ccw(i)); - - const Point& mp = midpoint(va->point(), vb->point()); - - Vertex_handle vm = insert_in_the_edge(f, i, mp); - - update_edges_to_be_conformed(va, vb, vm); - return vm; -} // # used by refine_face template inline @@ -1452,24 +693,23 @@ insert_in_the_edge(Face_handle fh, int edge_index, const Point& p) List_of_face_handles zone_of_p; // deconstrain the edge before finding the conflicts - // fh->set_constraint(edge_index,false); - // fh->neighbor(edge_index)->set_constraint(fh->mirror_index(edge_index),false); + fh->set_constraint(edge_index,false); + fh->neighbor(edge_index)->set_constraint(fh->mirror_index(edge_index),false); get_conflicts_and_boundary(p, std::back_inserter(zone_of_p), - Emptyset_iterator(), - EDGE, fh, edge_index); + Emptyset_iterator(), fh); // reconstrain the edge - // fh->set_constraint(edge_index,true); - // fh->neighbor(edge_index)->set_constraint(fh->mirror_index(edge_index),true); - + fh->set_constraint(edge_index,true); + fh->neighbor(edge_index)->set_constraint(fh->mirror_index(edge_index),true); + for(typename List_of_face_handles::iterator fh_it = zone_of_p.begin(); fh_it != zone_of_p.end(); ++fh_it) bad_faces.erase(*fh_it); - Vertex_handle vp = insert(p, EDGE, fh, edge_index); + Vertex_handle vp = insert(p, Base::EDGE, fh, edge_index); // TODO, WARNING: this is not robust! // We should deconstrained the constrained edge, insert the two // subconstraints and re-constrain them @@ -1477,11 +717,11 @@ insert_in_the_edge(Face_handle fh, int edge_index, const Point& p) // CGAL_assertion(is_bad_faces_valid()); is_edge(va, vp, fh, edge_index); - // set fh to the face at the right of [va,vb] + // set fh to the face at the right of [va,vp] Face_circulator fc = incident_faces(vp, fh), fcbegin(fc); // circulators are counter-clockwise, so we start at the right of - // [va,vb] + // [va,vp] do { if( !is_infinite(fc) ) fc->set_marked(mark_at_right); @@ -1493,51 +733,11 @@ insert_in_the_edge(Face_handle fh, int edge_index, const Point& p) fc->set_marked(mark_at_left); ++fc; } while ( fc != fcbegin ); - compute_new_bad_faces(vp); + return vp; } -// # used by: is_small_angle, create_clusters_of_vertex -// # compute 4 times the square of the cosine of the angle (ab,ac) -// # WARNING, TODO: this is not exact with doubles and can lead to crashes!! -template -typename Mesh_2::FT Mesh_2:: -squared_cosine_of_angle_times_4(const Point& pb, const Point& pa, - const Point& pc) const -{ - Compute_squared_distance_2 squared_distance = - geom_traits().compute_squared_distance_2_object(); - - const FT - a = squared_distance(pb, pc), - b = squared_distance(pa, pb), - c = squared_distance(pa, pc); - - const FT num = a-(b+c); - - return (num*num)/(b*c); -} - -// ->traits? -//the shortest edge that are in a triangle -// # used by: refine_face, squared_minimum_sine -template -typename Mesh_2::FT Mesh_2:: -shortest_edge_squared_length(Face_handle f) -{ - Compute_squared_distance_2 squared_distance = - geom_traits().compute_squared_distance_2_object(); - const Point - & pa = (f->vertex(0))->point(), - & pb = (f->vertex(1))->point(), - & pc = (f->vertex(2))->point(); - FT a, b, c; - a = squared_distance(pb, pc); - b = squared_distance(pc, pa); - c = squared_distance(pa, pb); - return (min(a, min(b, c))); -} template bool Mesh_2:: @@ -1580,7 +780,27 @@ is_bad_faces_valid() return result; } +// ->traits? +//the shortest edge that are in a triangle +// # used by: refine_face, squared_minimum_sine +template +typename Mesh_2::FT +Mesh_2:: +shortest_edge_squared_length(Face_handle f) +{ + Compute_squared_distance_2 squared_distance = + geom_traits().compute_squared_distance_2_object(); + const Point + & pa = (f->vertex(0))->point(), + & pb = (f->vertex(1))->point(), + & pc = (f->vertex(2))->point(); + FT a, b, c; + a = squared_distance(pb, pc); + b = squared_distance(pc, pa); + c = squared_distance(pa, pb); + return (min(a, min(b, c))); +} + CGAL_END_NAMESPACE - -#endif +#endif // CGAL_MESH_2_H diff --git a/Packages/Mesh_2/include/CGAL/Read_write.h b/Packages/Mesh_2/include/CGAL/Read_write.h new file mode 100644 index 00000000000..85a8f4c5027 --- /dev/null +++ b/Packages/Mesh_2/include/CGAL/Read_write.h @@ -0,0 +1,174 @@ +#include +#include + +#include +// to skip comments and EOF in the function read_poly + + + +// // IO + +// // write and read the constrained edges in the format: +// // number_of_edges +// // segment1 +// // segment2 +// // ... +// void write(std::ostream &f) const; +// void read(std::istream &f, bool dont_refine = false); + +// // write and read a mesh in the Triangle .poly format +// // (see http://www-2.cs.cmu.edu/~quake/triangle.poly.html) +// void write_poly(std::ostream &f) const; +// void read_poly(std::istream &f, bool dont_refine = false); + +CGAL_BEGIN_NAMESPACE + +// IO +//the function that writes a file +template +void +write(const Conform_triangulation_2& mesh, std::ostream &f) +{ + typedef typename Conform_triangulation_2::Finite_edges_iterator + Finite_edges_iterator; + + f << mesh.number_of_constrained_edges() << std::endl; + for(Finite_edges_iterator eit = mesh.finite_edges_begin(); + eit!=mesh.finite_edges_end(); + ++eit) + if((*eit).first->is_constrained((*eit).second)) + { + f << (*eit).first->vertex(mesh.cw((*eit).second))->point() << " " + << (*eit).first->vertex(mesh.ccw((*eit).second))->point() < +void +read(Conform_triangulation_2& mesh, std::istream &f) +{ + typedef typename Conform_triangulation_2::Point Point; + int nedges = 0; + mesh.clear(); + f>>nedges; + for(int n = 0; n> p1 >> p2; + mesh.insert_constraint(p1, p2); + } +} + +//the function that write a Shewchuk Triangle .poly file +template +void +write_poly(const Conform_triangulation_2& mesh, std::ostream &f) +{ + typedef Conform_triangulation_2 Triangulation; + typedef typename Triangulation::Vertex_handle Vertex_handle; + typedef typename Triangulation::Finite_vertices_iterator + Finite_vertices_iterator; + typedef typename Triangulation::Finite_edges_iterator + Finite_edges_iterator; + + + std::map index; + + // write vertices + f << "# Shewchuk Triangle .poly file, produced by the CGAL::Mesh_2 package" + << std::endl + << "# Neither attributes nor boundary markers are used." << std::endl + << mesh.number_of_vertices() << " " << 2 << " " + << 0 << " " << 0 << std::endl; + + f << std::endl; + + unsigned int vertices_counter = 0; + for(Finite_vertices_iterator vit = mesh.finite_vertices_begin(); + vit != mesh.finite_vertices_end(); + ++vit) + { + f << ++vertices_counter << " " << vit->point() << std::endl; + index[vit] = vertices_counter; + } + + f << std::endl; + + // write constrained edges + + f << mesh.number_of_constrained_edges() << " " << 0 << std::endl; + unsigned int edges_counter = 0; + for(Finite_edges_iterator eit = mesh.finite_edges_begin(); + eit != mesh.finite_edges_end(); + ++eit) + if((*eit).first->is_constrained((*eit).second)) + f << ++edges_counter << " " + << index[(*eit).first->vertex(mesh.cw((*eit).second))] << " " + << index[(*eit).first->vertex(mesh.ccw((*eit).second))] + << std::endl; + + f << std::endl; + +// // write seeds, assuming that the seeds unmarks faces +// unsigned int seeds_counter = 0; +// f << mesh.seeds.size() << std::endl; +// for(typename Seeds::const_iterator sit = seeds.begin(); +// sit!=seeds.end(); ++sit) +// f << ++seeds_counter << " " << *sit << std::endl; +} + +//the function that reads a Shewchuk Triangle .poly file +template +void +read_poly(Conform_triangulation_2& mesh, std::istream &f) +{ + typedef Conform_triangulation_2 Triangulation; + typedef typename Triangulation::Vertex_handle Vertex_handle; + typedef typename Triangulation::Point Point; + + mesh.clear(); + + unsigned int number_of_points; + skip_comment_OFF(f); + f >> number_of_points; + skip_until_EOL(f); + skip_comment_OFF(f); + + // read vertices + std::vector vertices(number_of_points); + for(unsigned int i = 0; i < number_of_points; ++i) + { + unsigned int j; + Point p; + f >> j >> p; + skip_until_EOL(f); skip_comment_OFF(f); + vertices[--j] = mesh.insert(p); + } + + // read segments + unsigned int number_of_segments; + f >> number_of_segments; + skip_until_EOL(f); skip_comment_OFF(f); + for(unsigned int k = 0; k < number_of_segments; ++k) + { + unsigned int l, v1, v2; + f >> l >> v1 >> v2; + skip_until_EOL(f); skip_comment_OFF(f); + mesh.insert_constraint(vertices[--v1], vertices[--v2]); + } + + // read holes + unsigned int number_of_holes; + f >> number_of_holes; + for(unsigned int m = 0; m < number_of_holes; ++m) + { + unsigned int n; + Point p; + f >> n >> p; + skip_until_EOL(f); skip_comment_OFF(f); + //seeds.push_back(p); + } +} + +CGAL_END_NAMESPACE + diff --git a/Packages/Mesh_2/test/Mesh_2/test_filtred_container.C b/Packages/Mesh_2/test/Mesh_2/test_filtred_container.C index 9264edd7a2b..d66af0cc7a1 100644 --- a/Packages/Mesh_2/test/Mesh_2/test_filtred_container.C +++ b/Packages/Mesh_2/test/Mesh_2/test_filtred_container.C @@ -1,5 +1,5 @@ #include -#include +#include #include #include @@ -13,7 +13,7 @@ public: } }; -typedef CGAL::Filtred_container, Is_odd> List; +typedef CGAL::Filtered_container, Is_odd> List; int main(int, char**) {