bug fix for natural_neighbor_coordinates_2 using a small feature in the conflict zone functions of Delaunay_triangulation_2

Interpolation/include/CGAL/natural_neighbor_coordinates_2.h

@@ -155,55 +155,7 @@ natural_neighbor_coordinates_vertex_2(const Dt& dt,
 
   std::list<Edge> hole;
 
-  dt.get_boundary_of_conflicts(p, std::back_inserter(hole), fh);
+  dt.get_boundary_of_conflicts(p, std::back_inserter(hole), fh, false);
-
-  // The symbolic perturbation makes that the conflict zone
-  // is too big for the interpolation. As a consequence 
-  // we would have points with lambda = 0 in the output
-  // Therefore we filter them out again 
-  typedef typename std::list<Edge>::iterator Edge_iterator;
-  Edge_iterator it = hole.begin(), nit = it;
-  ++nit;
-  
-  do {
-    Point_2 p0 = it->first->vertex(dt.cw(it->second))->point();
-    Point_2 p1 = it->first->vertex(dt.ccw(it->second))->point();
-    Point_2 p2 = nit->first->vertex(dt.ccw(nit->second))->point();
-    if(dt.geom_traits().side_of_oriented_circle_2_object()(p, p0, p1, p2) 
-       == ON_ORIENTED_BOUNDARY){
-      Face_handle nh = it->first->neighbor(it->second);
-      int ni = nh->index(it->first->vertex(dt.ccw(it->second)));
-      Edge_iterator eit = hole.insert(it, Edge(nh, ni));
-      hole.erase(it);
-      hole.erase(nit);
-      it = eit;
-      nit = it;
-    } else {
-      it = nit;
-    }
-    ++nit;
-  }while (nit != hole.end());
-
-  // now a similar test for the last and first edge
-  
-  nit = hole.begin();
-  do {
-    Point_2 p0 = it->first->vertex(dt.cw(it->second))->point();
-    Point_2 p1 = it->first->vertex(dt.ccw(it->second))->point();
-    Point_2 p2 = nit->first->vertex(dt.ccw(nit->second))->point();
-    if(dt.geom_traits().side_of_oriented_circle_2_object()(p, p0, p1, p2) 
-       == ON_ORIENTED_BOUNDARY){
-      Face_handle nh = it->first->neighbor(it->second);
-      int ni = nh->index(it->first->vertex(dt.ccw(it->second)));
-      Edge_iterator eit = hole.insert(it, Edge(nh, ni));
-      hole.erase(it);
-      hole.erase(nit);
-      it = eit;
-      nit = hole.begin();
-    } else {
-      ++it;
-    }
-  }while (it != hole.end());
 
   return natural_neighbor_coordinates_vertex_2
          (dt, p, out, hole.begin(), hole.end());

Triangulation_2/include/CGAL/Delaunay_triangulation_2.h

@@ -106,7 +106,7 @@ public:
   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 test_conflict(const Point  &p, Face_handle fh, bool strict = true) const;
   bool find_conflicts(const Point  &p,                //deprecated
 		      std::list<Face_handle>& conflicts,
 		      Face_handle start= Face_handle() ) const;
@@ -404,7 +404,8 @@ public:
   get_conflicts_and_boundary(const Point  &p, 
 			     OutputItFaces fit, 
 			     OutputItBoundaryEdges eit,
-			     Face_handle start = Face_handle()) const {
+			     Face_handle start = Face_handle(),
+                             bool strict = true) const {
     CGAL_triangulation_precondition( this->dimension() == 2);
     int li;
     Locate_type lt;
@@ -419,9 +420,9 @@ public:
       *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,0,pit, strict);
-      pit = propagate_conflicts(p,fh,1,pit);
+      pit = propagate_conflicts(p,fh,1,pit, strict);
-      pit = propagate_conflicts(p,fh,2,pit);
+      pit = propagate_conflicts(p,fh,2,pit, strict);
       return pit;    
     }
     CGAL_triangulation_assertion(false);
@@ -432,9 +433,10 @@ public:
   OutputItFaces
   get_conflicts (const Point  &p, 
 		 OutputItFaces fit, 
-		 Face_handle start= Face_handle()) const {
+		 Face_handle start= Face_handle(),
+                 bool strict = true) const {
     std::pair<OutputItFaces,Emptyset_iterator> pp = 
-      get_conflicts_and_boundary(p,fit,Emptyset_iterator(),start);
+      get_conflicts_and_boundary(p,fit,Emptyset_iterator(),start, strict);
     return pp.first;
   }
 
@@ -442,9 +444,10 @@ public:
   OutputItBoundaryEdges
   get_boundary_of_conflicts(const Point  &p, 
 			    OutputItBoundaryEdges eit, 
-			    Face_handle start= Face_handle()) const {
+			    Face_handle start= Face_handle(),
+                            bool strict = true) const {
     std::pair<Emptyset_iterator, OutputItBoundaryEdges> pp = 
-      get_conflicts_and_boundary(p,Emptyset_iterator(),eit,start);
+      get_conflicts_and_boundary(p,Emptyset_iterator(),eit,start,strict);
     return pp.second;
   }
 
@@ -460,15 +463,16 @@ private:
 		       const Face_handle fh, 
 		       const int i,
 		       std::pair<OutputItFaces,OutputItBoundaryEdges>
-		       pit)  const {
+		       pit,
+                       bool strict = true)  const {
     Face_handle fn = fh->neighbor(i);
-    if (! test_conflict(p,fn)) {
+    if (! test_conflict(p,fn,strict)) {
       *(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,ccw(j),pit,strict);
-      pit = propagate_conflicts(p,fn,cw(j), pit);
+      pit = propagate_conflicts(p,fn,cw(j), pit,strict);
     }
     return pit;
   }
@@ -478,7 +482,8 @@ private:
   non_recursive_propagate_conflicts ( const Point  &p,
                                       const Face_handle fh, 
                                       const int i,
-		                      std::pair<OutputItFaces,OutputItBoundaryEdges> pit)  const 
+		                      std::pair<OutputItFaces,OutputItBoundaryEdges> pit,
+                                      bool strict = true)  const 
   {
     std::stack<std::pair<Face_handle, int> > stack;
     stack.push( std::make_pair(fh,i) );
@@ -488,7 +493,7 @@ private:
       const int i=stack.top().second;
       stack.pop();
       Face_handle fn = fh->neighbor(i);
-      if (! test_conflict(p,fn)) {
+      if (! test_conflict(p,fn,strict)) {
         *(pit.second)++ = Edge(fn, fn->index(fh));
       } else {
         *(pit.first)++ = fn;
@@ -507,19 +512,20 @@ private:
 		       const int i,
 		       std::pair<OutputItFaces,OutputItBoundaryEdges>
 		       pit,
+                       bool strict = true,
                        int depth=0)  const 
   {
     if (depth == 100)
-      return non_recursive_propagate_conflicts(p, fh, i, pit);
+      return non_recursive_propagate_conflicts(p, fh, i, pit, strict);
 
     Face_handle fn = fh->neighbor(i);
-    if (! test_conflict(p,fn)) {
+    if (! test_conflict(p,fn,strict)) {
       *(pit.second)++ = Edge(fn, fn->index(fh));
     } else {
       *(pit.first)++ = fn;
       int j = fn->index(fh);
-      pit = propagate_conflicts(p,fn,ccw(j),pit,depth+1);
+      pit = propagate_conflicts(p,fn,ccw(j),pit, strict, depth+1);
-      pit = propagate_conflicts(p,fn,cw(j), pit,depth+1);
+      pit = propagate_conflicts(p,fn,cw(j), pit, strict, depth+1);
     }
     return pit;
   }
@@ -618,8 +624,12 @@ protected:
 template < class Gt, class Tds >
 inline bool
 Delaunay_triangulation_2<Gt,Tds>::
-test_conflict(const Point  &p, Face_handle fh) const
+test_conflict(const Point  &p, Face_handle fh, bool strict) const
 {
+  if(! strict){
+    Oriented_side os = side_of_oriented_circle(fh,p,false);
+    return os == ON_POSITIVE_SIDE;
+  }
   // 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