CGAL users discussion list

["Sebastien Loriot (GeometryFactory)" <>]

Thanks for the nice pictures.

With this patched file, you should be able to run
natural_neighbor_coordinates_2 on a object of type
Constrained_Delaunay_triangulation_2.

I am not sure of what you will get near constrained edges.

S.

mturner wrote:
> Hello,
>
> We are using the CGAL libraries to help analyze 
> surface motions of geological features on the Earth.  
> The geometry of the features define the boundary 
> and interior of specific surface regions.
>
> Because our geographic regions are small, 
> and to keep things simple for now,
> we consider the positions and motions only in 2D, 
> using coordinates of Latitude and Longitude.
>
> Each geologic feature is defined by a basic geometry (point, line, or
> polygon),
> and each geometry contributes to either the boundary of the region, or the
> interior.
>
> Following the example code in Chapter 35 "2D Triangulations", 
> we use the vertices from the feature geometry, to build a series of
> triangulations:
>
> - Delaunay Triangulation    ( CGAL::Delaunay_triangulation_2<...> ),
> - Constrained Triangulation (
> CGAL::Constrained_Delaunay_triangulation_2<...> ),
> - Mesh Triangulation        ( CGAL::Delaunay_mesher_2<...> ).
>
> We have posted images showing an example region,
> with the various feature geometry components,
> and the resulting triangulations here:
>
> http://web.gps.caltech.edu/~mturner/cgal_post/
>
> Figure 1 shows the basic feature data:
> the boundary of the region is defined by intersections of five lines
> (two red lines at the top, purple line on the right, red line at the bottom,
> green line on left).
> The Interior elements are defined by three points, three lines, and two
> polygons inside that region.
>
> Figure 2 shows the Delaunay Triangulation from that data in black,
> including the convex hull triangles formed where the boundary is concave.
>
> Figure 3 shows the Constrained Triangulation in gray.
>
> Figure 4 shows the Mesh Triangulation in light green, 
> with areas outside the total outer boundary un-meshed;
> and areas inside the interior polygons un-meshed.  
>
> The interior holes were generated by adding seed points and calling 
> the constrained mesher's set_seeds(...) function.
>
> Then the Mesh was generated with a call to:
> constrained_mesher->refine_mesh();
>
> This aspect of using CGAL to create triangulations and meshes 
> has been very successful for both test examples and real world data.
>
>
> Now we would like to use the final meshes to perform interpolation.
>
> We have measurements of surface velocity for the basic feature data,
> and we would like to use a refined CGAL mesh to interpolate 
> that velocity data for any given test point inside the region.
>
> Following the example code in Chapter 68 "2D and Surface Function
> Interpolation",
> we have successfully used the interpolation functions on the Delaunay
> Triangulation, 
> using the functions such as:
> CGAL::natural_neighbor_coordinates_2, 
> CGAL::sibson_gradient_fitting, 
> CGAL::Data_access, etc.
>
> Figure 5 shows the measured velocity at the feature vertices with white
> arrows, 
> and the interpolated velocity at the sample test points with black arrows.
>
> Unfortunately, interpolating with the Delaunay Triangulation 
> allows for interpolation outside the region (as shown by the left black
> arrow),
> and inside the holes of the mesh (as shown inside).
>
> We would like to use the Mesh Triangulation for interpolation, 
> not only because it is more refined,
> but also because the mesh is able to exclude, or 'mask out', 
> the interior polygon areas and external concave areas.
>
> However, I cannot find an example using a CGAL Mesh to interpolate.
>
> Is there a way to connect the CGAL interpolation code with the mesh code?
>
> It seems like there must be a way to utilize the refined mesh 
> to perform interpolation, but I could not find a way to make it work...
>
>
> Thanks for any help or guidance you can provide,
>
> Mark Turner
> California Institute of Technology
> Division of Geological and Planetary Sciences, Mail Code 100-23
> Pasadena, CA 91125
> Phone: 626-395-4543

// 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: svn+ssh://scm.gforge.inria.fr/svn/cgal/branches/CGAL-3.7-branch/Triangulation_2/include/CGAL/Constrained_Delaunay_triangulation_2.h $
// $Id: Constrained_Delaunay_triangulation_2.h 57014 2010-06-23 13:29:04Z afabri $
// 
//
// Author(s)     : Mariette Yvinec, Jean Daniel Boissonnat
 
#ifndef CGAL_CONSTRAINED_DELAUNAY_TRIANGULATION_2_H
#define CGAL_CONSTRAINED_DELAUNAY_TRIANGULATION_2_H

#include <CGAL/triangulation_assertions.h>
#include <CGAL/Constrained_triangulation_2.h>

namespace CGAL {


template <class Gt, 
          class Tds = Triangulation_data_structure_2 <
                      Triangulation_vertex_base_2<Gt>,
		      Constrained_triangulation_face_base_2<Gt> >,
	  class Itag = No_intersection_tag >		
class Constrained_Delaunay_triangulation_2
  : public  Constrained_triangulation_2<Gt, Tds, Itag> 
{
public:
  typedef Constrained_triangulation_2<Gt,Tds,Itag>             Ctr;
  typedef Constrained_Delaunay_triangulation_2<Gt,Tds,Itag>    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::size_type             size_type;
  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;

#ifndef CGAL_CFG_USING_BASE_MEMBER_BUG_2
  using Ctr::geom_traits;
  using Ctr::number_of_vertices;
  using Ctr::finite_faces_begin;
  using Ctr::finite_faces_end;
  using Ctr::dimension;
  using Ctr::cw;
  using Ctr::ccw;
  using Ctr::infinite_vertex;
  using Ctr::side_of_oriented_circle;
  using Ctr::is_infinite;
  using Ctr::collinear_between;
  using Ctr::are_there_incident_constraints;
  using Ctr::make_hole;
  using Ctr::insert_constraint;
  using Ctr::locate;
  using Ctr::test_dim_down;
  using Ctr::fill_hole_delaunay;
  using Ctr::update_constraints;
  using Ctr::delete_vertex;
  using Ctr::push_back;
#endif

  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<class InputIterator>
  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, bool perturb = true) 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<Edge>& le,  //deprecated
		      Face_handle hint= 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;
   

  // INSERTION-REMOVAL
  Vertex_handle insert(const Point & a, 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& a);
//   template < class InputIterator >
//   std::ptrdiff_t insert(InputIterator first, InputIterator last);

  void remove(Vertex_handle v);
  void remove_incident_constraints(Vertex_handle v);
  void remove_constrained_edge(Face_handle f, int i);
//  template <class OutputItFaces>
//  OutputItFaces
//  remove_constrained_edge(Face_handle f, int i, OutputItFaces out)
 
  //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);}
  void remove_constraint(Face_handle f, int i){remove_constrained_edge(f,i);}

  // 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());
  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)
  std::ptrdiff_t insert(InputIterator first, InputIterator last, int i = 0)
#else
    std::ptrdiff_t insert(InputIterator first, InputIterator last) 
#endif
    {
      size_type n = number_of_vertices();

      std::vector<Point> points (first, last);
      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 number_of_vertices() - n;
    }

  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( dimension() == 2);
    int li;
    Locate_type lt;
    Face_handle fh = locate(p,lt,li, start);
    switch(lt) {
    case Ctr::OUTSIDE_AFFINE_HULL:
    case Ctr::VERTEX:
      return std::make_pair(fit,eit);
    case Ctr::FACE:
    case Ctr::EDGE:
    case Ctr::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;
  }


public:
// made  public for the need of Mesh_2
// but not documented
 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 ( fh->is_constrained(i) || ! 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;
 }


public:
 template <class OutputItFaces>
 OutputItFaces
 propagating_flip(List_edges & edges, 
		  OutputItFaces out = Emptyset_iterator()) {
  // 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),this->mirror_index(f,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=this->mirror_index(f,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;
      eni=Edge(ff->neighbor(ii),this->mirror_index(ff,ii));
      if (less_edge(e[i],eni)) {edge_set.erase(e[i]);}
      else { edge_set.erase(eni);} 
    } 

    // here is the flip 
    *out++ = f;
    *out++ = f->neighbor(indf);
    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),this->mirror_index(ff,ii));
	if (less_edge(e[i],eni)) { 
	  edge_set.insert(e[i]);}
	else {
	  edge_set.insert(eni);} 
      }
    } 
  }
  return out;
 }

 template <class OutputItFaces>
 OutputItFaces
 remove_constrained_edge(Face_handle f, int i, OutputItFaces out) {
  Ctr::remove_constrained_edge(f,i);
  if(dimension() == 2) {
    List_edges le;
    le.push_back(Edge(f,i));
    propagating_flip(le,out);
  }
  return out;  
 }
 
  inline
  Point dual (Face_handle f) const
  {
    CGAL_triangulation_precondition(this->_tds.is_face(f));
    CGAL_triangulation_precondition (this->dimension()==2);
    return this->circumcenter(f);
  }
 
};


template < class Gt, class Tds, class Itag >
bool 
Constrained_Delaunay_triangulation_2<Gt,Tds,Itag>::
is_flipable(Face_handle f, int i, bool perturb) 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(), perturb) 
                                        == ON_POSITIVE_SIDE);
}

template < class Gt, class Tds, class Itag >
void 
Constrained_Delaunay_triangulation_2<Gt,Tds,Itag>::
flip (Face_handle& f, int i)
{
  Face_handle g = f->neighbor(i);
  int j = this->mirror_index(f,i);

  // save wings neighbors to be able to restore contraint status
  Face_handle f1 = f->neighbor(cw(i));
  int i1 = this->mirror_index(f,cw(i));
  Face_handle f2 = f->neighbor(ccw(i));
  int i2 = this->mirror_index(f,ccw(i));
  Face_handle f3 = g->neighbor(cw(j));
  int i3 = this->mirror_index(g,cw(j));
  Face_handle f4 = g->neighbor(ccw(j));
  int i4 = this->mirror_index(g,ccw(j));

  // The following precondition prevents the test suit 
  // of triangulation to work on constrained Delaunay triangulation
  //CGAL_triangulation_precondition(is_flipable(f,i));
  this->_tds.flip(f, i);
   
  // restore constraint status
  f->set_constraint(f->index(g), false);
  g->set_constraint(g->index(f), false);
  f1->neighbor(i1)->set_constraint(this->mirror_index(f1,i1),
				   f1->is_constrained(i1));
  f2->neighbor(i2)->set_constraint(this->mirror_index(f2,i2),
				   f2->is_constrained(i2));
  f3->neighbor(i3)->set_constraint(this->mirror_index(f3,i3),
				   f3->is_constrained(i3));
  f4->neighbor(i4)->set_constraint(this->mirror_index(f4,i4),
				   f4->is_constrained(i4));
  return;
}

template < class Gt, class Tds, class Itag >
void 
Constrained_Delaunay_triangulation_2<Gt,Tds,Itag>::
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<Gt,Tds,Itag>::
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<Gt,Tds,Itag>::
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<Gt,Tds,Itag>:: 
 propagating_flip(List_edges & edges) {
    propagating_flip(edges,Emptyset_iterator());
 }


template < class Gt, class Tds, class Itag >
inline bool
Constrained_Delaunay_triangulation_2<Gt,Tds,Itag>::
test_conflict(const Point& p, Face_handle fh) const
  // true if point P lies inside the circle circumscribing face fh
{
  // 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,true);
  if (os == ON_POSITIVE_SIDE) return true;
 
  if (os == ON_ORIENTED_BOUNDARY && is_infinite(fh)) {
    int i = fh->index(infinite_vertex());
    return collinear_between(fh->vertex(cw(i))->point(), p,
			     fh->vertex(ccw(i))->point() );
  }

  return false;
}

template < class Gt, class Tds, class Itag >
inline bool
Constrained_Delaunay_triangulation_2<Gt,Tds,Itag>::
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<Gt,Tds,Itag>::
find_conflicts(const Point& p, std::list<Edge>& 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<Gt,Tds,Itag>::Vertex_handle 
Constrained_Delaunay_triangulation_2<Gt,Tds,Itag>::
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<Gt,Tds,Itag>::Vertex_handle 
Constrained_Delaunay_triangulation_2<Gt,Tds,Itag>::
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<Gt,Tds,Itag>::Vertex_handle 
Constrained_Delaunay_triangulation_2<Gt,Tds,Itag>::
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<Gt,Tds,Itag>::Vertex_handle
Constrained_Delaunay_triangulation_2<Gt,Tds,Itag>::
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<Gt,Tds,Itag>::Vertex_handle
Constrained_Delaunay_triangulation_2<Gt,Tds,Itag>::
push_back(const Point &p)
{
  return insert(p);
}

template < class Gt, class Tds, class Itag >
void 
Constrained_Delaunay_triangulation_2<Gt,Tds,Itag>::
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<Gt,Tds,Itag>::
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 != Vertex_handle() );
  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<Gt,Tds,Itag>::Vertex_handle
// Constrained_Delaunay_triangulation_2<Gt,Tds,Itag>::
// 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<Gt,Tds,Itag>::
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);
    std::list<Edge> 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<Gt,Tds,Itag>::
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<Gt,Tds,Itag>::
remove_constrained_edge(Face_handle f, int i)
{
  remove_constrained_edge(f,i,Emptyset_iterator());
  return;  
}


template < class Gt, class Tds, class Itag >  
bool 
Constrained_Delaunay_triangulation_2<Gt,Tds,Itag>::
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, false);
	CGAL_triangulation_assertion( result );
      }
    }
    return result;
}



} //namespace CGAL
#endif // CGAL_CONSTRAINED_DELAUNAY_TRIANGULATION_2_H