From e8cebcb6b11d7f6e71d43db424ad9cdf22dece5a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?=
 <mael.rouxel.labbe@geometryfactory.com>
Date: Tue, 18 Apr 2017 12:21:55 +0200
Subject: [PATCH] Fixed surface_neighbor_coordinates

Short version: the regular triangulation must be build using weighted points,
but the return type is (bare) points
---
 .../CGAL/surface_neighbor_coordinates_3.h     | 79 +++++++++++++++----
 1 file changed, 63 insertions(+), 16 deletions(-)

diff --git a/Interpolation/include/CGAL/surface_neighbor_coordinates_3.h b/Interpolation/include/CGAL/surface_neighbor_coordinates_3.h
index d2fd520ba87..64a0296f6e4 100644
--- a/Interpolation/include/CGAL/surface_neighbor_coordinates_3.h
+++ b/Interpolation/include/CGAL/surface_neighbor_coordinates_3.h
@@ -62,13 +62,36 @@ surface_neighbor_coordinates_3(InputIterator
 
   //build Voronoi intersection triangulation:
   I_triangulation it(traits);
-  it.insert(first,beyond);
 
-  return regular_neighbor_coordinates_2(it, p, out);
+  typename ITraits::Construct_weighted_point_2 p2wp =
+      it.geom_traits().construct_weighted_point_2_object();
+  typename ITraits::Construct_point_2 wp2p =
+      it.geom_traits().construct_point_2_object();
+
+  while(first != beyond){
+    it.insert(p2wp(*first++));
+  }
+
+  typedef std::vector<std::pair<typename ITraits::Weighted_point_2,
+                                typename ITraits::FT> > WPoint_coordinate_vector;
+  WPoint_coordinate_vector wpoint_coords;
+
+  CGAL::Triple<std::back_insert_iterator<WPoint_coordinate_vector>,
+               typename ITraits::FT, bool> res =
+    regular_neighbor_coordinates_2(it, p2wp(p), std::back_inserter(wpoint_coords));
+
+  // regular_neighbor_coordinates_2 returns weighted_point_2 (in fact, _3)
+  // but surface_neighbor_coordinates_3 returns point_3 so the weight must be
+  // dropped
+  typename WPoint_coordinate_vector::const_iterator wpit = wpoint_coords.begin(),
+                                                    wpend = wpoint_coords.end();
+  for(; wpit!=wpend; ++wpit)
+    *out++ = std::make_pair(wp2p(wpit->first), wpit->second);
+
+  return CGAL::make_triple(out, res.second, res.third);
 }
 
-
-//without Delaunay filtering but with certification:
+// Without Delaunay filtering but with certification:
 // a boolean is returned that indicates if a sufficiently large
 // neighborhood has been considered so that the
 // Voronoi cell of p is not affected by any point outside the smallest
@@ -164,7 +187,15 @@ surface_neighbor_coordinates_certified_3(InputIterator
 
   //build Voronoi intersection triangulation:
   I_triangulation it(traits);
-  it.insert(first,beyond);
+
+  typename ITraits::Construct_weighted_point_2 p2wp =
+      it.geom_traits().construct_weighted_point_2_object();
+  typename ITraits::Construct_point_2 wp2p =
+      it.geom_traits().construct_point_2_object();
+
+  while(first != beyond){
+    it.insert(p2wp(*first++));
+  }
 
   //collect the Voronoi vertices of the cell of p in order to
   //determine the furthest distance from p to the boundary of its cell
@@ -173,22 +204,38 @@ surface_neighbor_coordinates_certified_3(InputIterator
   //unfortunately, there is no function call without Face_handle
   // "start" because this would cause type conflicts because
   // of resembling function signatures (-> default constructor)
-  Triple< OutputIterator, typename ITraits::FT, bool >
-    res = regular_neighbor_coordinates_2
-    (it, p, out, std::back_inserter(vor_vertices),
-     typename I_triangulation::Face_handle());
+  typename ITraits::Weighted_point_2 wp =
+    traits.construct_weighted_point_2_object()(p);
+
+  typedef std::vector<std::pair<typename ITraits::Weighted_point_2,
+                                typename ITraits::FT> > WPoint_coordinate_vector;
+  WPoint_coordinate_vector wpoint_coords;
+
+  Triple< std::back_insert_iterator<WPoint_coordinate_vector>,
+          typename ITraits::FT, bool > res =
+      regular_neighbor_coordinates_2(it, wp, std::back_inserter(wpoint_coords),
+                                     std::back_inserter(vor_vertices),
+                                     typename I_triangulation::Face_handle());
+
+  // regular_neighbor_coordinates_2 returns weighted_point_2 (in fact, _3)
+  // but surface_neighbor_coordinates_3 returns point_3 so the weight must be
+  // dropped
+  typename WPoint_coordinate_vector::const_iterator wpit = wpoint_coords.begin(),
+                                                    wpend = wpoint_coords.end();
+  for(; wpit!=wpend; ++wpit)
+    *out++ = std::make_pair(wp2p(wpit->first), wpit->second);
 
-  typename ITraits::Point_2 furthest =
-    *std::max_element(vor_vertices.begin(),  vor_vertices.end(),
-		     closer_to_point<ITraits>(p,traits));
   // if the distance to the furthest sample point is smaller
   // than twice the distance to the furthest vertex, not all neighbors
   // might be found: return false
-  if(radius < 4* traits.compute_squared_distance_2_object()
-     (p, furthest))
-    return make_quadruple(res.first, res.second, res.third, false);
+  typename ITraits::Point_2 furthest =
+    *std::max_element(vor_vertices.begin(), vor_vertices.end(),
+                      closer_to_point<ITraits>(p,traits));
 
-  return make_quadruple(res.first, res.second,res.third, true);
+  if(radius < 4* traits.compute_squared_distance_2_object()(p, furthest))
+    return make_quadruple(out, res.second, res.third, false);
+
+  return make_quadruple(out, res.second,res.third, true);
 }
 
 // FIXME :