From 48185ac1533f327ffbd87ac8d3ccfc51610d763b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?=
 <mael.rouxel.labbe@geometryfactory.com>
Date: Mon, 14 Jan 2019 11:20:16 +0100
Subject: [PATCH] Filter out irrelevant (<= 0) natural and regular neighbor
 coordinates

- coord == 0 is possible in theory, but not interesting to output
- coord < 0 is not possible in theory, but since we do construction and use
  signed (polygon_area_2()) area computations, we could in theory get
  negative coordinates in output
---
 .../CGAL/natural_neighbor_coordinates_2.h     | 34 +++++++++++++------
 .../CGAL/regular_neighbor_coordinates_2.h     | 14 +++++---
 2 files changed, 33 insertions(+), 15 deletions(-)

diff --git a/Interpolation/include/CGAL/natural_neighbor_coordinates_2.h b/Interpolation/include/CGAL/natural_neighbor_coordinates_2.h
index 0a18eb7d1a6..d2df19d72f0 100644
--- a/Interpolation/include/CGAL/natural_neighbor_coordinates_2.h
+++ b/Interpolation/include/CGAL/natural_neighbor_coordinates_2.h
@@ -99,21 +99,28 @@ natural_neighbors_2(const Dt& dt,
     Point_2 p1(v1->point()), p2(v2->point());
 
     Coord_type coef1(0);
-    Coord_type coef2(0);
     Equal_x_2 equal_x_2;
     if(!equal_x_2(p1,p2))
-    {
       coef1 = (p.x() - p2.x()) / (p1.x() - p2.x());
-      coef2 = 1 - coef1;
-      *out++ = std::make_pair(v1,coef1);
-      *out++ = std::make_pair(v2,coef2);
-    } else {
+    else
       coef1 = (p.y() - p2.y()) / (p1.y() - p2.y());
-      coef2 = 1 - coef1;
-      *out++ = std::make_pair(v1,coef1);
-      *out++ = std::make_pair(v2,coef2);
+
+    if(coef1 == 0)
+    {
+      *out++ = std::make_pair(v2, Coord_type(1));
+      return make_triple(out, Coord_type(1), true);
     }
 
+    Coord_type coef2 = 1 - coef1;
+    if(coef2 == 0)
+    {
+      *out++ = std::make_pair(v1, Coord_type(1));
+      return make_triple(out, Coord_type(1), true);
+    }
+
+    *out++ = std::make_pair(v1,coef1);
+    *out++ = std::make_pair(v2,coef2);
+
     return make_triple(out, coef1+coef2, true);
   }
 
@@ -142,6 +149,7 @@ natural_neighbors_2(const Dt& dt,
                     EdgeIterator hole_begin, EdgeIterator hole_end)
 {
   CGAL_precondition(dt.dimension() == 2);
+
   typedef typename Dt::Geom_traits       Traits;
   typedef typename Traits::FT            Coord_type;
   typedef typename Traits::Point_2       Point_2;
@@ -187,8 +195,12 @@ natural_neighbors_2(const Dt& dt,
                                                                 p);
 
     area += polygon_area_2(vor.begin(), vor.end(), dt.geom_traits());
-    *out++ = std::make_pair(current,area);
-    area_sum += area;
+
+    if(area > 0)
+    {
+      *out++ = std::make_pair(current,area);
+      area_sum += area;
+    }
 
     //update prev and hit:
     prev = current;
diff --git a/Interpolation/include/CGAL/regular_neighbor_coordinates_2.h b/Interpolation/include/CGAL/regular_neighbor_coordinates_2.h
index bbb03227e5c..5a83dd079c1 100644
--- a/Interpolation/include/CGAL/regular_neighbor_coordinates_2.h
+++ b/Interpolation/include/CGAL/regular_neighbor_coordinates_2.h
@@ -131,9 +131,12 @@ regular_neighbor_coordinates_vertex_2(const Rt& rt,
     *vor_vertices++ = vor[2];
 
     area += polygon_area_2(vor.begin(), vor.end(), rt.geom_traits());
-    *out++= std::make_pair(current, area);
 
-    area_sum += area;
+    if(area > 0)
+    {
+      *out++= std::make_pair(current, area);
+      area_sum += area;
+    }
 
     //update prev and hit:
     prev = current;
@@ -163,8 +166,11 @@ regular_neighbor_coordinates_vertex_2(const Rt& rt,
       ++fc;
     }
 
-    *out++ = std::make_pair((*hidden_vertices_begin), area);
-    area_sum += area;
+    if(area > 0)
+    {
+      *out++ = std::make_pair((*hidden_vertices_begin), area);
+      area_sum += area;
+    }
   }
 
   return make_triple(out, area_sum, true);