songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity
cgal/Triangulation_2/include/CGAL/Delaunay_triangulation_2.h

729 lines
21 KiB
C++
Raw Blame History

// Copyright (c) 1997 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org); you may redistribute it under
// the terms of the Q Public License version 1.0.
// See the file LICENSE.QPL distributed with CGAL.
//
// Licensees holding a valid commercial license may use this file in
// accordance with the commercial license agreement provided with the software.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
//
// $URL$
// $Id$
//
//
// Author(s) : Mariette Yvinec
#ifndef CGAL_DELAUNAY_TRIANGULATION_2_H
#define CGAL_DELAUNAY_TRIANGULATION_2_H
#include <CGAL/Triangulation_2.h>
#include <CGAL/iterator.h>
CGAL_BEGIN_NAMESPACE
template < class Gt,
class Tds = Triangulation_data_structure_2 <
Triangulation_vertex_base_2<Gt> > >
class Delaunay_triangulation_2 : public Triangulation_2<Gt,Tds>
{
public:
typedef Gt Geom_traits;
typedef typename Geom_traits::Point_2 Point;
typedef typename Geom_traits::Segment_2 Segment;
typedef typename Geom_traits::Triangle_2 Triangle;
typedef typename Geom_traits::Orientation_2 Orientation_2;
typedef typename Geom_traits::Compare_x_2 Compare_x;
typedef typename Geom_traits::Compare_y_2 Compare_y;
typedef typename Geom_traits::Side_of_oriented_circle_2
Side_of_oriented_circle;
typedef Triangulation_2<Gt,Tds> Triangulation;
typedef typename Triangulation::Locate_type Locate_type;
typedef typename Triangulation::Face_handle Face_handle;
typedef typename Triangulation::Vertex_handle Vertex_handle;
typedef typename Triangulation::Edge Edge;
typedef typename Triangulation::Edge_circulator Edge_circulator;
typedef typename Triangulation::Face_circulator Face_circulator;
typedef typename Triangulation::Vertex_circulator Vertex_circulator;
typedef typename Triangulation::Finite_edges_iterator Finite_edges_iterator;
typedef typename Triangulation::Finite_faces_iterator Finite_faces_iterator;
typedef typename Triangulation::Finite_vertices_iterator
Finite_vertices_iterator;
#ifndef CGAL_CFG_USING_BASE_MEMBER_BUG_2
using Triangulation::side_of_oriented_circle;
using Triangulation::circumcenter;
#endif
Delaunay_triangulation_2(const Gt& gt = Gt())
: Triangulation_2<Gt,Tds>(gt) {}
Delaunay_triangulation_2(
const Delaunay_triangulation_2<Gt,Tds> &tr)
: Triangulation_2<Gt,Tds>(tr)
{ CGAL_triangulation_postcondition( is_valid() ); }
// CHECK -QUERY
bool is_valid(bool verbose = false, int level = 0) const;
Vertex_handle
nearest_vertex(const Point& p, Face_handle f= Face_handle()) const;
bool does_conflict(const Point &p, Face_handle fh) const;// deprecated
bool test_conflict(const Point &p, Face_handle fh) const;
bool find_conflicts(const Point &p, //deprecated
std::list<Face_handle>& conflicts,
Face_handle start= Face_handle() ) const;
// //template member functions, declared and defined at the end
// template <class OutputItFaces, class OutputItBoundaryEdges>
// std::pair<OutputItFaces,OutputItBoundaryEdges>
// get_conflicts_and_boundary(const Point &p,
// OutputItFaces fit,
// OutputItBoundaryEdges eit,
// Face_handle start) const;
// template <class OutputItFaces>
// OutputItFaces
// get_conflicts (const Point &p,
// OutputItFaces fit,
// Face_handle start ) const;
// template <class OutputItBoundaryEdges>
// OutputItBoundaryEdges
// get_boundary_of_conflicts(const Point &p,
// OutputItBoundaryEdges eit,
// Face_handle start ) const;
// DUAL
Point dual (Face_handle f) const;
Object dual(const Edge &e) const ;
Object dual(const Edge_circulator& ec) const;
Object dual(const Finite_edges_iterator& ei) const;
//INSERTION-REMOVAL
Vertex_handle insert(const Point &p,
Face_handle start = Face_handle() );
Vertex_handle insert(const Point& p,
Locate_type lt,
Face_handle loc, int li );
Vertex_handle push_back(const Point &p);
void remove(Vertex_handle v );
void restore_Delaunay(Vertex_handle v);
// MOVE
bool move(Vertex_handle v, const Point &p);
#ifndef CGAL_CFG_USING_BASE_MEMBER_BUG_2
using Triangulation::cw;
using Triangulation::ccw;
using Triangulation::geom_traits;
#endif
private:
void propagating_flip(Face_handle& f,int i);
void remove_2D(Vertex_handle v );
Vertex_handle nearest_vertex_2D(const Point& p, Face_handle f) const;
Vertex_handle nearest_vertex_1D(const Point& p) const;
void look_nearest_neighbor(const Point& p,
Face_handle f,
int i,
Vertex_handle& nn) const;
public:
template < class Stream>
Stream& draw_dual(Stream & ps)
{
Finite_edges_iterator eit= this->finite_edges_begin();
for (; eit != this->finite_edges_end(); ++eit) {
Object o = dual(eit);
typename Geom_traits::Line_2 l;
typename Geom_traits::Ray_2 r;
Segment s;
if (CGAL::assign(s,o)) ps << s;
if (CGAL::assign(r,o)) ps << r;
if (CGAL::assign(l,o)) ps << l;
}
return ps;
}
template < class InputIterator >
int
insert(InputIterator first, InputIterator last)
{
int n = this->number_of_vertices();
std::vector<Point> points (first, last);
std::random_shuffle (points.begin(), points.end());
spatial_sort (points.begin(), points.end(), geom_traits());
Face_handle f;
for (typename std::vector<Point>::const_iterator p = points.begin(), end = points.end();
p != end; ++p)
f = insert (*p, f)->face();
return this->number_of_vertices() - n;
}
template < class InputIterator >
bool move(InputIterator first, InputIterator last)
{
bool blocked = false;
std::map<Vertex_handle, int> hash;
std::list< std::pair<Vertex_handle, Point> > to_move(first, last);
while(!to_move.empty()) {
std::pair<Vertex_handle, Point> pp = to_move.front();
to_move.pop_front();
if(!move(pp.first, pp.second)) {
if(hash[pp.first] == 3) break;
else if(hash[pp.first] == 2) blocked = true;
hash[pp.first]++;
to_move.push_back(pp);
}
}
return !blocked;
}
template <class OutputItFaces, class OutputItBoundaryEdges>
std::pair<OutputItFaces,OutputItBoundaryEdges>
get_conflicts_and_boundary(const Point &p,
OutputItFaces fit,
OutputItBoundaryEdges eit,
Face_handle start = Face_handle()) const {
CGAL_triangulation_precondition( this->dimension() == 2);
int li;
Locate_type lt;
Face_handle fh = locate(p,lt,li, start);
switch(lt) {
case Triangulation::OUTSIDE_AFFINE_HULL:
case Triangulation::VERTEX:
return std::make_pair(fit,eit);
case Triangulation::FACE:
case Triangulation::EDGE:
case Triangulation::OUTSIDE_CONVEX_HULL:
*fit++ = fh; //put fh in OutputItFaces
std::pair<OutputItFaces,OutputItBoundaryEdges>
pit = std::make_pair(fit,eit);
pit = propagate_conflicts(p,fh,0,pit);
pit = propagate_conflicts(p,fh,1,pit);
pit = propagate_conflicts(p,fh,2,pit);
return pit;
}
CGAL_triangulation_assertion(false);
return std::make_pair(fit,eit);
}
template <class OutputItFaces>
OutputItFaces
get_conflicts (const Point &p,
OutputItFaces fit,
Face_handle start= Face_handle()) const {
std::pair<OutputItFaces,Emptyset_iterator> pp =
get_conflicts_and_boundary(p,fit,Emptyset_iterator(),start);
return pp.first;
}
template <class OutputItBoundaryEdges>
OutputItBoundaryEdges
get_boundary_of_conflicts(const Point &p,
OutputItBoundaryEdges eit,
Face_handle start= Face_handle()) const {
std::pair<Emptyset_iterator, OutputItBoundaryEdges> pp =
get_conflicts_and_boundary(p,Emptyset_iterator(),eit,start);
return pp.second;
}
private:
template <class OutputItFaces, class OutputItBoundaryEdges>
std::pair<OutputItFaces,OutputItBoundaryEdges>
propagate_conflicts (const Point &p,
Face_handle fh,
int i,
std::pair<OutputItFaces,OutputItBoundaryEdges>
pit) const {
Face_handle fn = fh->neighbor(i);
if (! test_conflict(p,fn)) {
*(pit.second)++ = Edge(fn, fn->index(fh));
} else {
*(pit.first)++ = fn;
int j = fn->index(fh);
pit = propagate_conflicts(p,fn,ccw(j),pit);
pit = propagate_conflicts(p,fn,cw(j), pit);
}
return pit;
}
protected:
void restore_edges(Vertex_handle v)
{
std::list<Edge> edges;
Face_circulator fc = this->incident_faces(v), done(fc);
int degree = 0;
do {
if((++degree) > 3) break;
} while(++fc != done);
fc = this->incident_faces(v);
done = fc;
if(degree == 3) {
do {
int i = fc->index(v);
edges.push_back(Edge(fc, i));
} while(++fc != done);
} else {
do {
int i = fc->index(v);
edges.push_back(Edge(fc, i));
edges.push_back(Edge(fc, this->cw(i)));
} while(++fc != done);
}
while(!edges.empty()) {
const Edge &e = edges.front();
Face_handle f = e.first;
int i = e.second;
edges.pop_front();
if(this->is_infinite(f->vertex(i))) continue;
Face_handle fi = f->neighbor(i);
int mi = this->_tds.mirror_index(f, i);
Vertex_handle vm = this->_tds.mirror_vertex(f, i);
if(this->is_infinite(vm)) continue;
if(this->side_of_oriented_circle(f, vm->point()) == ON_POSITIVE_SIDE) {
this->_tds.flip(f, i);
edges.push_back(Edge(f, i));
edges.push_back(Edge(f, this->cw(i)));
edges.push_back(Edge(fi, this->cw(mi)));
edges.push_back(Edge(fi, mi));
}
}
}
};
template < class Gt, class Tds >
inline bool
Delaunay_triangulation_2<Gt,Tds>::
test_conflict(const Point &p, Face_handle fh) const
{
// return true if P is inside the circumcircle of fh
// if fh is infinite, return true when p is in the positive
// halfspace or on the boundary and in the finite edge of fh
Oriented_side os = side_of_oriented_circle(fh,p);
if (os == ON_POSITIVE_SIDE) return true;
if (os == ON_ORIENTED_BOUNDARY && is_infinite(fh)) {
int i = fh->index(this->infinite_vertex());
return collinear_between(fh->vertex(cw(i))->point(), p,
fh->vertex(ccw(i))->point() );
}
return false;
}
template < class Gt, class Tds >
inline bool
Delaunay_triangulation_2<Gt,Tds>::
does_conflict(const Point &p, Face_handle fh) const
{
return test_conflict(p,fh);
}
template < class Gt, class Tds >
inline bool
Delaunay_triangulation_2<Gt,Tds>::
find_conflicts(const Point &p,
std::list<Face_handle>& conflicts,
Face_handle start ) const
{
get_conflicts(p, std::back_inserter(conflicts), start);
return (! conflicts.empty());
}
template < class Gt, class Tds >
bool
Delaunay_triangulation_2<Gt,Tds>::
is_valid(bool verbose, int level) const
{
bool result = Triangulation_2<Gt,Tds>::is_valid(verbose, level);
for( Finite_faces_iterator it = this->finite_faces_begin();
it != this->finite_faces_end() ; it++) {
for(int i=0; i<3; i++) {
if ( ! is_infinite( this->mirror_vertex(it,i))) {
result = result && ON_POSITIVE_SIDE !=
side_of_oriented_circle( it, this->mirror_vertex(it,i)->point());
}
CGAL_triangulation_assertion( result );
}
}
return result;
}
template < class Gt, class Tds >
typename Delaunay_triangulation_2<Gt,Tds>::Vertex_handle
Delaunay_triangulation_2<Gt,Tds>::
nearest_vertex(const Point &p, Face_handle f) const
{
switch (this->dimension()) {
case 0:
if (this->number_of_vertices() == 0) return Vertex_handle();
if (this->number_of_vertices() == 1) return this->finite_vertex();
//break;
case 1:
return nearest_vertex_1D(p);
//break;
case 2:
return nearest_vertex_2D(p,f);
//break;
}
return Vertex_handle();
}
template < class Gt, class Tds >
typename Delaunay_triangulation_2<Gt,Tds>::Vertex_handle
Delaunay_triangulation_2<Gt,Tds>::
nearest_vertex_2D(const Point& p, Face_handle f) const
{
CGAL_triangulation_precondition(this->dimension() == 2);
f = locate(p,f);
typename Geom_traits::Compare_distance_2
compare_distance = this->geom_traits().compare_distance_2_object();
Vertex_handle nn = !is_infinite(f->vertex(0)) ? f->vertex(0):f->vertex(1);
if ( !is_infinite(f->vertex(1)) && compare_distance(p,
f->vertex(1)->point(),
nn->point()) == SMALLER)
nn=f->vertex(1);
if ( !is_infinite(f->vertex(2)) && compare_distance(p,
f->vertex(2)->point(),
nn->point()) == SMALLER)
nn=f->vertex(2);
look_nearest_neighbor(p,f,0,nn);
look_nearest_neighbor(p,f,1,nn);
look_nearest_neighbor(p,f,2,nn);
return nn;
}
template < class Gt, class Tds >
typename Delaunay_triangulation_2<Gt,Tds>::Vertex_handle
Delaunay_triangulation_2<Gt,Tds>::
nearest_vertex_1D(const Point& p) const
{
typename Geom_traits::Compare_distance_2
compare_distance = this->geom_traits().compare_distance_2_object();
Vertex_handle nn;
Finite_vertices_iterator vit=this->finite_vertices_begin();
nn = vit;
for ( ; vit != this->finite_vertices_end(); ++vit){
if (compare_distance(p, vit->point(), nn->point()) == SMALLER)
nn = vit;
}
return nn;
}
template < class Gt, class Tds >
void
Delaunay_triangulation_2<Gt,Tds>::
look_nearest_neighbor(const Point& p,
Face_handle f,
int i,
Vertex_handle& nn) const
{
Face_handle ni=f->neighbor(i);
if ( ON_POSITIVE_SIDE != side_of_oriented_circle(ni,p) ) return;
typename Geom_traits::Compare_distance_2
compare_distance = this->geom_traits().compare_distance_2_object();
i = ni->index(f);
if ( !is_infinite(ni->vertex(i)) &&
compare_distance(p,
ni->vertex(i)->point(),
nn->point()) == SMALLER) nn=ni->vertex(i);
// recursive exploration of triangles whose circumcircle contains p
look_nearest_neighbor(p, ni, ccw(i), nn);
look_nearest_neighbor(p, ni, cw(i), nn);
}
//DUALITY
template<class Gt, class Tds>
inline
typename Delaunay_triangulation_2<Gt,Tds>::Point
Delaunay_triangulation_2<Gt,Tds>::
dual (Face_handle f) const
{
CGAL_triangulation_precondition (this->dimension()==2);
return circumcenter(f);
}
template < class Gt, class Tds >
Object
Delaunay_triangulation_2<Gt,Tds>::
dual(const Edge &e) const
{
typedef typename Geom_traits::Line_2 Line;
typedef typename Geom_traits::Ray_2 Ray;
CGAL_triangulation_precondition (!is_infinite(e));
if( this->dimension()== 1 ){
const Point& p = (e.first)->vertex(cw(e.second))->point();
const Point& q = (e.first)->vertex(ccw(e.second))->point();
Line l = this->geom_traits().construct_bisector_2_object()(p,q);
return make_object(l);
}
// dimension==2
if( (!is_infinite(e.first)) &&
(!is_infinite(e.first->neighbor(e.second))) ) {
Segment s = this->geom_traits().construct_segment_2_object()
(dual(e.first),dual(e.first->neighbor(e.second)));
return make_object(s);
}
// one of the adjacent faces is infinite
Face_handle f; int i;
if (is_infinite(e.first)) {
f=e.first->neighbor(e.second); i=f->index(e.first);
}
else {
f=e.first; i=e.second;
}
const Point& p = f->vertex(cw(i))->point();
const Point& q = f->vertex(ccw(i))->point();
Line l = this->geom_traits().construct_bisector_2_object()(p,q);
Ray r = this->geom_traits().construct_ray_2_object()(dual(f), l);
return make_object(r);
}
template < class Gt, class Tds >
inline Object
Delaunay_triangulation_2<Gt,Tds>::
dual(const Edge_circulator& ec) const
{
return dual(*ec);
}
template < class Gt, class Tds >
inline Object
Delaunay_triangulation_2<Gt,Tds>::
dual(const Finite_edges_iterator& ei) const
{
return dual(*ei);
}
template < class Gt, class Tds >
inline
typename Delaunay_triangulation_2<Gt,Tds>::Vertex_handle
Delaunay_triangulation_2<Gt,Tds>::
insert(const Point &p, Face_handle start)
{
Locate_type lt;
int li;
Face_handle loc = locate (p, lt, li, start);
return insert(p, lt, loc, li);
}
template < class Gt, class Tds >
inline
typename Delaunay_triangulation_2<Gt,Tds>::Vertex_handle
Delaunay_triangulation_2<Gt,Tds>::
push_back(const Point &p)
{
return insert(p);
}
template < class Gt, class Tds >
inline
typename Delaunay_triangulation_2<Gt,Tds>::Vertex_handle
Delaunay_triangulation_2<Gt,Tds>::
insert(const Point &p, Locate_type lt, Face_handle loc, int li)
{
Vertex_handle v = Triangulation_2<Gt,Tds>::insert(p,lt,loc,li);
restore_Delaunay(v);
return(v);
}
template < class Gt, class Tds >
void
Delaunay_triangulation_2<Gt,Tds>::
restore_Delaunay(Vertex_handle v)
{
if(this->dimension() <= 1) return;
Face_handle f=v->face();
Face_handle next;
int i;
Face_handle start(f);
do {
i = f->index(v);
next = f->neighbor(ccw(i)); // turn ccw around v
propagating_flip(f,i);
f=next;
} while(next != start);
return;
}
template < class Gt, class Tds >
void
Delaunay_triangulation_2<Gt,Tds>::
remove(Vertex_handle v)
{
CGAL_triangulation_precondition( v != Vertex_handle());
CGAL_triangulation_precondition( !is_infinite(v));
if ( this->dimension() <= 1) Triangulation::remove(v);
else remove_2D(v);
return;
}
template < class Gt, class Tds >
void
Delaunay_triangulation_2<Gt,Tds>::
propagating_flip(Face_handle& f,int i)
{
Face_handle n = f->neighbor(i);
if ( ON_POSITIVE_SIDE !=
side_of_oriented_circle(n, f->vertex(i)->point()) ) {
return;
}
flip(f, i);
propagating_flip(f,i);
i = n->index(f->vertex(i));
propagating_flip(n,i);
}
template <class Gt, class Tds >
void
Delaunay_triangulation_2<Gt,Tds>::
remove_2D(Vertex_handle v)
{
if (test_dim_down(v)) { this->_tds.remove_dim_down(v); }
else {
std::list<Edge> hole;
make_hole(v, hole);
fill_hole_delaunay(hole);
delete_vertex(v);
}
return;
}
template <class Gt, class Tds >
bool
Delaunay_triangulation_2<Gt,Tds>::
move(Vertex_handle v, const Point &p) {
CGAL_triangulation_precondition(!is_infinite(v));
const int dim = this->dimension();
if(dim == 2) {
Point ant = v->point();
v->set_point(p);
if(well_oriented(v)) {
restore_edges(v);
return true;
}
v->set_point(ant);
}
Locate_type lt;
int li;
Vertex_handle inserted;
Face_handle loc = locate(p, lt, li, v->face());
if(lt == Triangulation_2<Gt,Tds>::VERTEX) return false;
if(dim < 0) return true;
if(dim == 0) {
v->point() = p;
return true;
}
if((loc != NULL) && (dim == 1)) {
if(loc->has_vertex(v)) {
v->point() = p;
} else {
inserted = insert(p, lt, loc, li);
Face_handle f = v->face();
int i = f->index(v);
if (i==0) {f = f->neighbor(1);}
CGAL_triangulation_assertion(f->index(v) == 1);
Face_handle g= f->neighbor(0);
f->set_vertex(1, g->vertex(1));
f->set_neighbor(0,g->neighbor(0));
g->neighbor(0)->set_neighbor(1,f);
g->vertex(1)->set_face(f);
this->delete_face(g);
Face_handle f_ins = inserted->face();
i = f_ins->index(inserted);
if (i==0) {f_ins = f_ins->neighbor(1);}
CGAL_triangulation_assertion(f_ins->index(inserted) == 1);
Face_handle g_ins = f_ins->neighbor(0);
f_ins->set_vertex(1, v);
g_ins->set_vertex(0, v);
std::swap(*v, *inserted);
this->delete_vertex(inserted);
}
return true;
}
if((loc != NULL) && this->test_dim_down(v)) {
v->point() = p;
int i = loc->index(v);
Face_handle locl;
int i_locl;
if(is_infinite(loc)) {
int i_inf = loc->index(this->infinite_vertex());
locl = loc->neighbor(i_inf);
i_locl = locl->index(v);
} else { locl = loc; i_locl = i; }
if(this->orientation(p, locl->vertex(ccw(i_locl))->point(),
locl->vertex(cw(i_locl))->point()) == COLLINEAR) {
this->_tds.dim_2D_1D(loc, i);
}
return true;
}
inserted = insert(p, lt, loc, li);
std::list<Edge> hole;
make_hole(v, hole);
fill_hole_delaunay(hole);
// fixing pointer
Face_circulator fc = incident_faces(inserted), done(fc);
std::list<Face_handle> faces_pt;
do { faces_pt.push_back(fc); } while(++fc != done);
while(!faces_pt.empty()) {
Face_handle f = faces_pt.front();
faces_pt.pop_front();
int i = f->index(inserted);
f->set_vertex(i, v);
}
std::swap(*v, *inserted);
this->delete_vertex(inserted);
return true;
}
CGAL_END_NAMESPACE
#endif // CGAL_DELAUNAY_TRIANGULATION_2_H
Reference in New Issue View Git Blame Copy Permalink