From cd0ccf07c4d85e9f5c8b03a1330504cc123f51b8 Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Thu, 4 Aug 2016 15:27:14 +0200
Subject: [PATCH 01/26] New algorithm: automatic scale selection for
 reconstruction

---
 .../Point_set_processing_3/CMakeLists.txt     |   1 +
 .../scale_estimation_example.cpp              |  39 ++++
 .../include/CGAL/estimate_scale.h             | 204 ++++++++++++++++++
 3 files changed, 244 insertions(+)
 create mode 100644 Point_set_processing_3/examples/Point_set_processing_3/scale_estimation_example.cpp
 create mode 100644 Point_set_processing_3/include/CGAL/estimate_scale.h

diff --git a/Point_set_processing_3/examples/Point_set_processing_3/CMakeLists.txt b/Point_set_processing_3/examples/Point_set_processing_3/CMakeLists.txt
index 3a432c08a1c..e9e37fd0808 100644
--- a/Point_set_processing_3/examples/Point_set_processing_3/CMakeLists.txt
+++ b/Point_set_processing_3/examples/Point_set_processing_3/CMakeLists.txt
@@ -53,6 +53,7 @@ if ( CGAL_FOUND )
   create_single_source_cgal_program( "read_write_xyz_point_set_example.cpp" )
   create_single_source_cgal_program( "read_ply_points_with_colors_example.cpp" )
   create_single_source_cgal_program( "remove_outliers_example.cpp" )
+  create_single_source_cgal_program( "scale_estimation_example.cpp" )
   create_single_source_cgal_program( "wlop_simplify_and_regularize_point_set_example.cpp" )
   create_single_source_cgal_program( "edge_aware_upsample_point_set_example.cpp" )
   
diff --git a/Point_set_processing_3/examples/Point_set_processing_3/scale_estimation_example.cpp b/Point_set_processing_3/examples/Point_set_processing_3/scale_estimation_example.cpp
new file mode 100644
index 00000000000..a01575817f7
--- /dev/null
+++ b/Point_set_processing_3/examples/Point_set_processing_3/scale_estimation_example.cpp
@@ -0,0 +1,39 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/estimate_scale.h>
+#include <CGAL/IO/read_xyz_points.h>
+
+#include <vector>
+#include <fstream>
+#include <boost/tuple/tuple.hpp>
+
+// Types
+typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
+typedef Kernel::FT FT;
+typedef Kernel::Point_3 Point;
+
+int main (int argc, char** argv)
+{
+  const char* fname = (argc>1)?argv[1]:"data/sphere_20k.xyz";
+  // Reads a .xyz point set file in points.
+  // As the point is the second element of the tuple (that is with index 1)
+  // we use a property map that accesses the 1st element of the tuple.
+    
+  std::vector<Point> points;
+  std::ifstream stream(fname);
+  if (!stream ||
+      !CGAL::read_xyz_points(
+                             stream, std::back_inserter(points)))
+    {
+      std::cerr << "Error: cannot read file " << fname << std::endl;
+      return EXIT_FAILURE;
+    }
+
+  // estimate global scale
+  std::size_t scale = CGAL::estimate_global_k_neighbor_scale (points.begin(), points.end(),
+                                                              CGAL::Identity_property_map<Point>(),
+                                                              Kernel());
+  std::cout << "Global K scale: " << scale << std::endl;
+
+  return EXIT_SUCCESS;
+}
+
diff --git a/Point_set_processing_3/include/CGAL/estimate_scale.h b/Point_set_processing_3/include/CGAL/estimate_scale.h
new file mode 100644
index 00000000000..971ae84dcff
--- /dev/null
+++ b/Point_set_processing_3/include/CGAL/estimate_scale.h
@@ -0,0 +1,204 @@
+#ifndef CGAL_ESTIMATE_SCALE_H
+#define CGAL_ESTIMATE_SCALE_H
+
+#include <CGAL/Search_traits_3.h>
+#include <CGAL/squared_distance_3.h>
+#include <CGAL/Orthogonal_k_neighbor_search.h>
+#include <CGAL/property_map.h>
+#include <CGAL/point_set_processing_assertions.h>
+#include <CGAL/assertions.h>
+#include <CGAL/hierarchy_simplify_point_set.h>
+
+#include <iterator>
+#include <list>
+
+#ifdef CGAL_LINKED_WITH_TBB
+#include <tbb/parallel_for.h>
+#include <tbb/blocked_range.h>
+#include <tbb/scalable_allocator.h>  
+#endif // CGAL_LINKED_WITH_TBB
+
+namespace CGAL {
+
+
+// ----------------------------------------------------------------------------
+// Private section
+// ----------------------------------------------------------------------------
+/// \cond SKIP_IN_MANUAL
+namespace internal {
+
+template <class Kernel>
+class Quick_multiscale_approximate_knn_distance
+{
+  typedef typename Kernel::FT FT;
+  typedef Search_traits_3<Kernel> Tree_traits;
+  typedef Orthogonal_k_neighbor_search<Tree_traits> Neighbor_search;
+  typedef typename Neighbor_search::Tree Tree;
+  typedef typename Neighbor_search::iterator Iterator;
+
+  std::size_t m_cluster_size;
+  std::vector<Tree*> m_trees;
+
+public:
+
+  template <typename InputIterator, typename PointPMap>
+  Quick_multiscale_approximate_knn_distance (InputIterator first,
+                                             InputIterator beyond,
+                                             PointPMap point_pmap,
+                                             std::size_t cluster_size = 10)
+    : m_cluster_size (cluster_size)
+  {
+    m_trees.push_back (new Tree ());
+    std::size_t nb_pts = 0;
+    for (InputIterator it = first; it != beyond; ++ it)
+      {
+        m_trees[0]->insert (get(point_pmap, *it));
+        ++ nb_pts;
+      }
+    
+    std::size_t nb_trees = 0;
+    while (nb_pts > m_cluster_size)
+      {
+        nb_trees ++;
+        nb_pts /= m_cluster_size;
+      }
+
+    m_trees.reserve (nb_trees);
+
+    InputIterator first_unused = beyond;
+
+    for (std::size_t i = 1; i < nb_trees; ++ i)
+      {
+        m_trees.push_back (new Tree());
+        first_unused
+          = CGAL::hierarchy_simplify_point_set (first, first_unused, point_pmap, cluster_size, 1./3.);
+
+        for (InputIterator it = first; it != first_unused; ++ it)
+          m_trees.back()->insert (get(point_pmap, *it));
+      }
+  }
+
+  ~Quick_multiscale_approximate_knn_distance()
+  {
+    for (std::size_t i = 0; i < m_trees.size(); ++ i)
+      delete m_trees[i];
+  }
+
+  template <typename InputIterator, typename PointPMap>
+  std::size_t compute_scale (InputIterator query, PointPMap point_pmap)
+  {
+    std::size_t out = 0;
+
+    std::size_t weight = 1;
+    FT dist_min = std::numeric_limits<FT>::max();
+    FT sum_sq_distances = 0.;
+    std::size_t nb = 0;
+    
+    for (std::size_t t = 0; t < m_trees.size(); ++ t)
+      {
+        Neighbor_search search (*(m_trees[t]), get(point_pmap, *query), m_cluster_size);
+        Iterator it = search.begin();
+        
+        if (t != 0) // Skip first point except on first scale
+          ++ it;
+
+        for (; it != search.end(); ++ it)
+          {
+            sum_sq_distances += weight * it->second;
+            nb += weight;
+            
+            if (nb < 6) // do not consider values under 6
+              continue;
+
+            FT dist = std::sqrt (sum_sq_distances / nb)
+              / std::pow (nb, 0.375); // nb^(5/12)
+
+            if (dist < dist_min)
+              {
+                dist_min = dist;
+                out = nb;
+              }
+          }
+        weight *= m_cluster_size;
+      }
+    return out;
+  }
+
+};
+
+} /* namespace internal */
+/// \endcond
+
+
+
+// ----------------------------------------------------------------------------
+// Public section
+// ----------------------------------------------------------------------------
+
+template <typename SamplesInputIterator,
+          typename SamplesPointPMap,
+          typename QueriesInputIterator,
+          typename QueriesPointPMap,
+          typename OutputIterator,
+          typename Kernel
+>
+void
+estimate_local_k_neighbor_scales(
+  SamplesInputIterator first,
+  SamplesInputIterator beyond,
+  SamplesPointPMap samples_pmap,
+  QueriesInputIterator first_query,
+  QueriesInputIterator beyond_query,
+  QueriesPointPMap queries_pmap,
+  OutputIterator output,
+  const Kernel& /*kernel*/) ///< geometric traits.
+{
+  // Build multi-scale KD-tree
+  internal::Quick_multiscale_approximate_knn_distance<Kernel> kdtree (first, beyond, samples_pmap);
+
+  // Compute local scales everywhere
+  for (QueriesInputIterator it = first_query; it != beyond_query; ++ it)
+    *(output ++) = kdtree.compute_scale (it, queries_pmap);
+
+}
+
+template <typename InputIterator,
+          typename PointPMap,
+          typename Kernel
+>
+std::size_t
+estimate_global_k_neighbor_scale(
+  InputIterator first,  ///< iterator over the first input point.
+  InputIterator beyond, ///< past-the-end iterator over the input points.
+  PointPMap point_pmap, ///< property map: value_type of InputIterator -> Point_3
+  const Kernel& kernel) ///< geometric traits.
+{
+  std::vector<std::size_t> scales;
+  estimate_local_k_neighbor_scales (first, beyond, point_pmap,
+                                    first, beyond, point_pmap,
+                                    std::back_inserter (scales),
+                                    kernel);
+  std::sort (scales.begin(), scales.end());
+  return scales[scales.size() / 2];
+}
+
+template <typename InputIterator,
+          typename PointPMap,
+          typename Kernel
+>
+typename Kernel::FT
+estimate_global_range_scale(
+  InputIterator first,  ///< iterator over the first input point.
+  InputIterator beyond, ///< past-the-end iterator over the input points.
+  PointPMap point_pmap, ///< property map: value_type of InputIterator -> Point_3
+  const Kernel& /*kernel*/) ///< geometric traits.
+{
+
+  return 0.;
+}
+
+
+
+} //namespace CGAL
+
+#endif // CGAL_ESTIMATE_SCALE_3_H

From 11d59640a91d5df08db65515055ed81b50f50173 Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Thu, 4 Aug 2016 15:34:56 +0200
Subject: [PATCH 02/26] Add possibility to look for a range scale (instead of a
 K scale)

---
 .../scale_estimation_example.cpp              | 13 ++--
 .../include/CGAL/estimate_scale.h             | 69 ++++++++++++++++---
 2 files changed, 70 insertions(+), 12 deletions(-)

diff --git a/Point_set_processing_3/examples/Point_set_processing_3/scale_estimation_example.cpp b/Point_set_processing_3/examples/Point_set_processing_3/scale_estimation_example.cpp
index a01575817f7..9f79d46f340 100644
--- a/Point_set_processing_3/examples/Point_set_processing_3/scale_estimation_example.cpp
+++ b/Point_set_processing_3/examples/Point_set_processing_3/scale_estimation_example.cpp
@@ -29,10 +29,15 @@ int main (int argc, char** argv)
     }
 
   // estimate global scale
-  std::size_t scale = CGAL::estimate_global_k_neighbor_scale (points.begin(), points.end(),
-                                                              CGAL::Identity_property_map<Point>(),
-                                                              Kernel());
-  std::cout << "Global K scale: " << scale << std::endl;
+  std::size_t k_scale = CGAL::estimate_global_k_neighbor_scale (points.begin(), points.end(),
+                                                                CGAL::Identity_property_map<Point>(),
+                                                                Kernel());
+  std::cout << "Global K scale: " << k_scale << std::endl;
+
+  FT range_scale = CGAL::estimate_global_range_scale (points.begin(), points.end(),
+                                                      CGAL::Identity_property_map<Point>(),
+                                                      Kernel());
+  std::cout << "Global range scale: " << range_scale << std::endl;
 
   return EXIT_SUCCESS;
 }
diff --git a/Point_set_processing_3/include/CGAL/estimate_scale.h b/Point_set_processing_3/include/CGAL/estimate_scale.h
index 971ae84dcff..e38f82cb014 100644
--- a/Point_set_processing_3/include/CGAL/estimate_scale.h
+++ b/Point_set_processing_3/include/CGAL/estimate_scale.h
@@ -85,9 +85,29 @@ public:
   }
 
   template <typename InputIterator, typename PointPMap>
-  std::size_t compute_scale (InputIterator query, PointPMap point_pmap)
+  std::size_t compute_k_scale (InputIterator query, PointPMap point_pmap)
   {
-    std::size_t out = 0;
+    std::size_t out;
+    FT dummy;
+    compute_scale (query, point_pmap, out, dummy);
+    return out;
+  }
+
+  template <typename InputIterator, typename PointPMap>
+  FT compute_range_scale (InputIterator query, PointPMap point_pmap)
+  {
+    std::size_t dummy;
+    FT out;
+    compute_scale (query, point_pmap, dummy, out);
+    return out;
+  }
+
+  template <typename InputIterator, typename PointPMap>
+  void compute_scale (InputIterator query, PointPMap point_pmap,
+                      std::size_t& k, FT& d)
+  {
+    k = 0;
+    d = 0.;
 
     std::size_t weight = 1;
     FT dist_min = std::numeric_limits<FT>::max();
@@ -116,12 +136,12 @@ public:
             if (dist < dist_min)
               {
                 dist_min = dist;
-                out = nb;
+                k = nb;
+                d = it->second;
               }
           }
         weight *= m_cluster_size;
       }
-    return out;
   }
 
 };
@@ -158,7 +178,7 @@ estimate_local_k_neighbor_scales(
 
   // Compute local scales everywhere
   for (QueriesInputIterator it = first_query; it != beyond_query; ++ it)
-    *(output ++) = kdtree.compute_scale (it, queries_pmap);
+    *(output ++) = kdtree.compute_k_scale (it, queries_pmap);
 
 }
 
@@ -182,6 +202,34 @@ estimate_global_k_neighbor_scale(
   return scales[scales.size() / 2];
 }
 
+  
+template <typename SamplesInputIterator,
+          typename SamplesPointPMap,
+          typename QueriesInputIterator,
+          typename QueriesPointPMap,
+          typename OutputIterator,
+          typename Kernel
+>
+void
+estimate_local_range_scales(
+  SamplesInputIterator first,
+  SamplesInputIterator beyond,
+  SamplesPointPMap samples_pmap,
+  QueriesInputIterator first_query,
+  QueriesInputIterator beyond_query,
+  QueriesPointPMap queries_pmap,
+  OutputIterator output,
+  const Kernel& /*kernel*/) ///< geometric traits.
+{
+  // Build multi-scale KD-tree
+  internal::Quick_multiscale_approximate_knn_distance<Kernel> kdtree (first, beyond, samples_pmap);
+
+  // Compute local scales everywhere
+  for (QueriesInputIterator it = first_query; it != beyond_query; ++ it)
+    *(output ++) = kdtree.compute_range_scale (it, queries_pmap);
+
+}
+
 template <typename InputIterator,
           typename PointPMap,
           typename Kernel
@@ -191,10 +239,15 @@ estimate_global_range_scale(
   InputIterator first,  ///< iterator over the first input point.
   InputIterator beyond, ///< past-the-end iterator over the input points.
   PointPMap point_pmap, ///< property map: value_type of InputIterator -> Point_3
-  const Kernel& /*kernel*/) ///< geometric traits.
+  const Kernel& kernel) ///< geometric traits.
 {
-
-  return 0.;
+  std::vector<typename Kernel::FT> scales;
+  estimate_local_range_scales (first, beyond, point_pmap,
+                               first, beyond, point_pmap,
+                               std::back_inserter (scales),
+                               kernel);
+  std::sort (scales.begin(), scales.end());
+  return std::sqrt (scales[scales.size() / 2]);
 }
 
 

From 5b7d0a99cbb3b2e120f243d68d79e23d717c2988 Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Fri, 5 Aug 2016 08:30:59 +0200
Subject: [PATCH 03/26] Add 2D case

---
 .../scale_estimation_example.cpp              |  77 +++++---
 .../include/CGAL/estimate_scale.h             | 172 +++++++++++++++++-
 2 files changed, 220 insertions(+), 29 deletions(-)

diff --git a/Point_set_processing_3/examples/Point_set_processing_3/scale_estimation_example.cpp b/Point_set_processing_3/examples/Point_set_processing_3/scale_estimation_example.cpp
index 9f79d46f340..d10028701ba 100644
--- a/Point_set_processing_3/examples/Point_set_processing_3/scale_estimation_example.cpp
+++ b/Point_set_processing_3/examples/Point_set_processing_3/scale_estimation_example.cpp
@@ -9,36 +9,65 @@
 // Types
 typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
 typedef Kernel::FT FT;
-typedef Kernel::Point_3 Point;
+typedef Kernel::Point_3 Point_3;
+typedef Kernel::Point_2 Point_2;
 
 int main (int argc, char** argv)
 {
-  const char* fname = (argc>1)?argv[1]:"data/sphere_20k.xyz";
-  // Reads a .xyz point set file in points.
-  // As the point is the second element of the tuple (that is with index 1)
-  // we use a property map that accesses the 1st element of the tuple.
-    
-  std::vector<Point> points;
-  std::ifstream stream(fname);
-  if (!stream ||
-      !CGAL::read_xyz_points(
-                             stream, std::back_inserter(points)))
-    {
-      std::cerr << "Error: cannot read file " << fname << std::endl;
-      return EXIT_FAILURE;
-    }
 
-  // estimate global scale
-  std::size_t k_scale = CGAL::estimate_global_k_neighbor_scale (points.begin(), points.end(),
-                                                                CGAL::Identity_property_map<Point>(),
-                                                                Kernel());
-  std::cout << "Global K scale: " << k_scale << std::endl;
+  // 3D CASE
+  {
+    const char* fname = (argc>1)?argv[1]:"data/sphere_20k.xyz";
+    // Reads a .xyz point set file in points.
+    // As the point is the second element of the tuple (that is with index 1)
+    // we use a property map that accesses the 1st element of the tuple.
+  
+    std::vector<Point_3> points;
+    std::ifstream stream(fname);
+    if (!stream ||
+        !CGAL::read_xyz_points(
+                               stream, std::back_inserter(points)))
+      {
+        std::cerr << "Error: cannot read file " << fname << std::endl;
+        return EXIT_FAILURE;
+      }
 
-  FT range_scale = CGAL::estimate_global_range_scale (points.begin(), points.end(),
-                                                      CGAL::Identity_property_map<Point>(),
-                                                      Kernel());
-  std::cout << "Global range scale: " << range_scale << std::endl;
+    // estimate global scale
+    std::size_t k_scale = CGAL::estimate_global_k_neighbor_scale (points.begin(), points.end(),
+                                                                  CGAL::Identity_property_map<Point_3>(),
+                                                                  Kernel());
+    std::cout << "Global K scale: " << k_scale << std::endl;
 
+    FT range_scale = CGAL::estimate_global_range_scale (points.begin(), points.end(),
+                                                        CGAL::Identity_property_map<Point_3>(),
+                                                        Kernel());
+    std::cout << "Global range scale: " << range_scale << std::endl;
+  }
+
+  // 2D CASE
+  {
+    std::vector<Point_2> points;
+    for (std::size_t i = 0; i < 5000; ++ i)
+      {
+        double angle = CGAL_PI * 2. * (rand() / (double)RAND_MAX);
+        double radius = 1. + 0.1 * (rand() / (double)RAND_MAX);
+        points.push_back (Point_2 (radius * std::cos (angle),
+                                   radius * std::sin (angle)));
+      }
+
+        
+    // estimate global scale
+    std::size_t k_scale = CGAL::estimate_global_k_neighbor_scale (points.begin(), points.end(),
+                                                                  CGAL::Identity_property_map<Point_2>(),
+                                                                  Kernel());
+    std::cout << "Global K scale: " << k_scale << std::endl;
+
+    FT range_scale = CGAL::estimate_global_range_scale (points.begin(), points.end(),
+                                                        CGAL::Identity_property_map<Point_2>(),
+                                                        Kernel());
+    std::cout << "Global range scale: " << range_scale << std::endl;
+
+  }
   return EXIT_SUCCESS;
 }
 
diff --git a/Point_set_processing_3/include/CGAL/estimate_scale.h b/Point_set_processing_3/include/CGAL/estimate_scale.h
index e38f82cb014..747de5b5522 100644
--- a/Point_set_processing_3/include/CGAL/estimate_scale.h
+++ b/Point_set_processing_3/include/CGAL/estimate_scale.h
@@ -8,6 +8,8 @@
 #include <CGAL/point_set_processing_assertions.h>
 #include <CGAL/assertions.h>
 #include <CGAL/hierarchy_simplify_point_set.h>
+#include <CGAL/random_simplify_point_set.h>
+#include <CGAL/Point_set_2.h>
 
 #include <iterator>
 #include <list>
@@ -27,8 +29,15 @@ namespace CGAL {
 /// \cond SKIP_IN_MANUAL
 namespace internal {
 
-template <class Kernel>
+template <class Kernel, class PointType>
 class Quick_multiscale_approximate_knn_distance
+{
+
+};
+
+  
+template <class Kernel>
+class Quick_multiscale_approximate_knn_distance<Kernel, typename Kernel::Point_3>
 {
   typedef typename Kernel::FT FT;
   typedef Search_traits_3<Kernel> Tree_traits;
@@ -131,7 +140,7 @@ public:
               continue;
 
             FT dist = std::sqrt (sum_sq_distances / nb)
-              / std::pow (nb, 0.375); // nb^(5/12)
+              / std::pow (nb, 0.41666); // nb^(5/12)
 
             if (dist < dist_min)
               {
@@ -146,6 +155,156 @@ public:
 
 };
 
+  
+template <class Kernel>
+class Quick_multiscale_approximate_knn_distance<Kernel, typename Kernel::Point_2>
+{
+  typedef typename Kernel::FT FT;
+  typedef CGAL::Point_set_2<Kernel> Point_set;
+  typedef typename Point_set::Vertex_handle Vertex_handle;
+
+  template <typename ValueType, typename PointPMap>
+  struct Pmap_unary_function : public std::unary_function<ValueType, typename Kernel::Point_2>
+  {
+    PointPMap point_pmap;
+    Pmap_unary_function (PointPMap point_pmap) : point_pmap (point_pmap) { }
+    const typename Kernel::Point_2& operator() (const ValueType& v) const { return get(point_pmap, v); }
+  };
+
+  template <typename PointPMap>
+  struct Pmap_to_3d
+  {
+    PointPMap point_pmap;
+    typedef typename Kernel::Point_3 value_type;
+    typedef const value_type& reference;
+    typedef typename boost::property_traits<PointPMap>::key_type key_type;
+    typedef boost::lvalue_property_map_tag category;
+    Pmap_to_3d () { }
+    Pmap_to_3d (PointPMap point_pmap)
+      : point_pmap (point_pmap) { }
+    friend inline value_type get (const Pmap_to_3d& ppmap, key_type i) 
+    {
+      typename Kernel::Point_2 p2 = get(ppmap.point_pmap, i);
+      return value_type (p2.x(), p2.y(), 0.);
+    }
+
+  };
+
+
+  std::size_t m_cluster_size;
+  std::vector<Point_set*> m_point_sets;
+  
+public:
+
+  template <typename InputIterator, typename PointPMap>
+  Quick_multiscale_approximate_knn_distance (InputIterator first,
+                                             InputIterator beyond,
+                                             PointPMap point_pmap,
+                                             std::size_t cluster_size = 10)
+    : m_cluster_size (cluster_size)
+  {
+    typedef Pmap_unary_function<typename std::iterator_traits<InputIterator>::value_type,
+                                PointPMap> Unary_f;
+
+    m_point_sets.push_back (new Point_set (boost::make_transform_iterator (first, Unary_f(point_pmap)),
+                                           boost::make_transform_iterator (beyond, Unary_f(point_pmap))));
+
+    std::size_t nb_pts = m_point_sets[0]->number_of_vertices();
+    std::size_t nb_trees = 0;
+    while (nb_pts > m_cluster_size)
+      {
+        nb_trees ++;
+        nb_pts /= m_cluster_size;
+      }
+
+    m_point_sets.reserve (nb_trees);
+
+    InputIterator first_unused = beyond;
+    nb_pts = m_point_sets[0]->number_of_vertices();
+
+    for (std::size_t i = 1; i < nb_trees; ++ i)
+      {
+        first_unused
+          = CGAL::hierarchy_simplify_point_set (first, beyond, Pmap_to_3d<PointPMap> (point_pmap),
+                                                m_cluster_size, 1./3.);
+
+        m_point_sets.push_back (new Point_set (boost::make_transform_iterator (first, Unary_f(point_pmap)),
+                                               boost::make_transform_iterator (first_unused, Unary_f(point_pmap))));
+
+        m_cluster_size *= cluster_size;
+      }
+    
+    m_cluster_size = cluster_size;
+  }
+
+  ~Quick_multiscale_approximate_knn_distance()
+  {
+    for (std::size_t i = 0; i < m_point_sets.size(); ++ i)
+      delete m_point_sets[i];
+  }
+
+  template <typename InputIterator, typename PointPMap>
+  std::size_t compute_k_scale (InputIterator query, PointPMap point_pmap)
+  {
+    std::size_t out;
+    FT dummy;
+    compute_scale (query, point_pmap, out, dummy);
+    return out;
+  }
+
+  template <typename InputIterator, typename PointPMap>
+  FT compute_range_scale (InputIterator query, PointPMap point_pmap)
+  {
+    std::size_t dummy;
+    FT out;
+    compute_scale (query, point_pmap, dummy, out);
+    return out;
+  }
+
+  template <typename InputIterator, typename PointPMap>
+  void compute_scale (InputIterator query, PointPMap point_pmap,
+                      std::size_t& k, FT& d)
+  {
+    k = 0;
+    d = 0.;
+
+    std::size_t weight = 1;
+    FT dist_min = std::numeric_limits<FT>::max();
+    FT sum_sq_distances = 0.;
+    std::size_t nb = 0;
+
+    const typename Kernel::Point_2& pquery = get(point_pmap, *query);
+    for (std::size_t t = 0; t < m_point_sets.size(); ++ t)
+      {
+        std::vector<Vertex_handle> neighbors;
+        m_point_sets[t]->nearest_neighbors (pquery, m_cluster_size,
+                                            std::back_inserter (neighbors));
+
+        for (std::size_t n = 1; n < neighbors.size(); ++ n)
+          {
+            FT sq_dist = CGAL::squared_distance (pquery, neighbors[n]->point());
+            sum_sq_distances += weight * sq_dist;
+            nb += weight;
+            
+            if (nb < 6) // do not consider values under 6
+              continue;
+
+            FT dist = std::sqrt (sum_sq_distances / nb)
+              / std::pow (nb, 0.75); // nb^(3/4)
+
+            if (dist < dist_min)
+              {
+                dist_min = dist;
+                k = nb;
+                d = sq_dist;
+              }
+          }
+        weight *= m_cluster_size;
+      }
+  }
+
+};
+
 } /* namespace internal */
 /// \endcond
 
@@ -173,13 +332,14 @@ estimate_local_k_neighbor_scales(
   OutputIterator output,
   const Kernel& /*kernel*/) ///< geometric traits.
 {
+  typedef typename boost::property_traits<SamplesPointPMap>::value_type Point_d;
+
   // Build multi-scale KD-tree
-  internal::Quick_multiscale_approximate_knn_distance<Kernel> kdtree (first, beyond, samples_pmap);
+  internal::Quick_multiscale_approximate_knn_distance<Kernel, Point_d> kdtree (first, beyond, samples_pmap);
 
   // Compute local scales everywhere
   for (QueriesInputIterator it = first_query; it != beyond_query; ++ it)
     *(output ++) = kdtree.compute_k_scale (it, queries_pmap);
-
 }
 
 template <typename InputIterator,
@@ -221,8 +381,10 @@ estimate_local_range_scales(
   OutputIterator output,
   const Kernel& /*kernel*/) ///< geometric traits.
 {
+  typedef typename boost::property_traits<SamplesPointPMap>::value_type Point_d;
+
   // Build multi-scale KD-tree
-  internal::Quick_multiscale_approximate_knn_distance<Kernel> kdtree (first, beyond, samples_pmap);
+  internal::Quick_multiscale_approximate_knn_distance<Kernel, Point_d> kdtree (first, beyond, samples_pmap);
 
   // Compute local scales everywhere
   for (QueriesInputIterator it = first_query; it != beyond_query; ++ it)

From cbad1f375cccf9acfe2e0f5bd4cc38d325e3d30b Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Fri, 5 Aug 2016 15:41:06 +0200
Subject: [PATCH 04/26] Increase precision by keeping track of clusters'
 weights

---
 .../include/CGAL/estimate_scale.h             | 119 +++++++++++-------
 1 file changed, 72 insertions(+), 47 deletions(-)

diff --git a/Point_set_processing_3/include/CGAL/estimate_scale.h b/Point_set_processing_3/include/CGAL/estimate_scale.h
index 747de5b5522..49b808e3064 100644
--- a/Point_set_processing_3/include/CGAL/estimate_scale.h
+++ b/Point_set_processing_3/include/CGAL/estimate_scale.h
@@ -11,14 +11,11 @@
 #include <CGAL/random_simplify_point_set.h>
 #include <CGAL/Point_set_2.h>
 
+#include <fstream>
+
 #include <iterator>
 #include <list>
 
-#ifdef CGAL_LINKED_WITH_TBB
-#include <tbb/parallel_for.h>
-#include <tbb/blocked_range.h>
-#include <tbb/scalable_allocator.h>  
-#endif // CGAL_LINKED_WITH_TBB
 
 namespace CGAL {
 
@@ -47,6 +44,7 @@ class Quick_multiscale_approximate_knn_distance<Kernel, typename Kernel::Point_3
 
   std::size_t m_cluster_size;
   std::vector<Tree*> m_trees;
+  std::vector<FT> m_weights;
 
 public:
 
@@ -58,6 +56,7 @@ public:
     : m_cluster_size (cluster_size)
   {
     m_trees.push_back (new Tree ());
+    m_weights.push_back (1.);
     std::size_t nb_pts = 0;
     for (InputIterator it = first; it != beyond; ++ it)
       {
@@ -73,9 +72,11 @@ public:
       }
 
     m_trees.reserve (nb_trees);
+    m_weights.reserve (nb_trees);
 
     InputIterator first_unused = beyond;
-
+    std::cerr << nb_trees << std::endl;
+    nb_pts = m_trees[0]->size();
     for (std::size_t i = 1; i < nb_trees; ++ i)
       {
         m_trees.push_back (new Tree());
@@ -84,6 +85,8 @@ public:
 
         for (InputIterator it = first; it != first_unused; ++ it)
           m_trees.back()->insert (get(point_pmap, *it));
+
+        m_weights.push_back (m_trees[0]->size() / (FT)(m_trees.back()->size()));
       }
   }
 
@@ -113,19 +116,22 @@ public:
 
   template <typename InputIterator, typename PointPMap>
   void compute_scale (InputIterator query, PointPMap point_pmap,
-                      std::size_t& k, FT& d)
+                      std::size_t& k, FT& d,
+                      std::ostream& log)
   {
     k = 0;
     d = 0.;
 
-    std::size_t weight = 1;
     FT dist_min = std::numeric_limits<FT>::max();
     FT sum_sq_distances = 0.;
     std::size_t nb = 0;
     
     for (std::size_t t = 0; t < m_trees.size(); ++ t)
       {
-        Neighbor_search search (*(m_trees[t]), get(point_pmap, *query), m_cluster_size);
+        Neighbor_search search (*(m_trees[t]), get(point_pmap, *query),
+                                (t == (m_trees.size() - 1)
+                                 ? m_trees[t]->size()
+                                 : m_weights[t+1] / m_weights[t]));
         Iterator it = search.begin();
         
         if (t != 0) // Skip first point except on first scale
@@ -133,15 +139,16 @@ public:
 
         for (; it != search.end(); ++ it)
           {
-            sum_sq_distances += weight * it->second;
-            nb += weight;
-            
+            sum_sq_distances += m_weights[t] * it->second;
+            nb += m_weights[t];
+
             if (nb < 6) // do not consider values under 6
               continue;
 
+            
             FT dist = std::sqrt (sum_sq_distances / nb)
               / std::pow (nb, 0.41666); // nb^(5/12)
-
+            log << nb << " " << dist << std::endl;
             if (dist < dist_min)
               {
                 dist_min = dist;
@@ -149,7 +156,6 @@ public:
                 d = it->second;
               }
           }
-        weight *= m_cluster_size;
       }
   }
 
@@ -190,9 +196,21 @@ class Quick_multiscale_approximate_knn_distance<Kernel, typename Kernel::Point_2
 
   };
 
+  struct Sort_by_distance_to_point
+  {
+    const typename Kernel::Point_2& ref;
+    Sort_by_distance_to_point (const typename Kernel::Point_2& ref) : ref (ref) { }
+    bool operator() (const Vertex_handle& a, const Vertex_handle& b)
+    {
+      return (CGAL::squared_distance (a->point(), ref)
+              < CGAL::squared_distance (b->point(), ref));
+    }
+  };
+
 
   std::size_t m_cluster_size;
   std::vector<Point_set*> m_point_sets;
+  std::vector<FT> m_weights;
   
 public:
 
@@ -200,7 +218,7 @@ public:
   Quick_multiscale_approximate_knn_distance (InputIterator first,
                                              InputIterator beyond,
                                              PointPMap point_pmap,
-                                             std::size_t cluster_size = 10)
+                                             std::size_t cluster_size = 20)
     : m_cluster_size (cluster_size)
   {
     typedef Pmap_unary_function<typename std::iterator_traits<InputIterator>::value_type,
@@ -208,7 +226,8 @@ public:
 
     m_point_sets.push_back (new Point_set (boost::make_transform_iterator (first, Unary_f(point_pmap)),
                                            boost::make_transform_iterator (beyond, Unary_f(point_pmap))));
-
+    m_weights.push_back (1.);
+    
     std::size_t nb_pts = m_point_sets[0]->number_of_vertices();
     std::size_t nb_trees = 0;
     while (nb_pts > m_cluster_size)
@@ -218,22 +237,23 @@ public:
       }
 
     m_point_sets.reserve (nb_trees);
+    m_weights.reserve (nb_trees);
 
     InputIterator first_unused = beyond;
     nb_pts = m_point_sets[0]->number_of_vertices();
-
+    std::cerr << nb_trees << std::endl;
     for (std::size_t i = 1; i < nb_trees; ++ i)
       {
         first_unused
-          = CGAL::hierarchy_simplify_point_set (first, beyond, Pmap_to_3d<PointPMap> (point_pmap),
+          = CGAL::hierarchy_simplify_point_set (first, first_unused, Pmap_to_3d<PointPMap> (point_pmap),
                                                 m_cluster_size, 1./3.);
 
         m_point_sets.push_back (new Point_set (boost::make_transform_iterator (first, Unary_f(point_pmap)),
                                                boost::make_transform_iterator (first_unused, Unary_f(point_pmap))));
 
-        m_cluster_size *= cluster_size;
+        m_weights.push_back (nb_pts / (FT)(m_point_sets.back()->number_of_vertices()));
       }
-    
+    std::cerr << std::endl;
     m_cluster_size = cluster_size;
   }
 
@@ -268,25 +288,31 @@ public:
     k = 0;
     d = 0.;
 
-    std::size_t weight = 1;
     FT dist_min = std::numeric_limits<FT>::max();
     FT sum_sq_distances = 0.;
-    std::size_t nb = 0;
+    FT nb = 0;
 
     const typename Kernel::Point_2& pquery = get(point_pmap, *query);
     for (std::size_t t = 0; t < m_point_sets.size(); ++ t)
       {
         std::vector<Vertex_handle> neighbors;
-        m_point_sets[t]->nearest_neighbors (pquery, m_cluster_size,
-                                            std::back_inserter (neighbors));
+        if (t == m_point_sets.size() - 1)
+          m_point_sets[t]->nearest_neighbors (pquery, m_point_sets[t]->number_of_vertices(),
+                                              std::back_inserter (neighbors));
+        else
+          m_point_sets[t]->nearest_neighbors (pquery, m_weights[t+1] / m_weights[t],
+                                              std::back_inserter (neighbors));
 
-        for (std::size_t n = 1; n < neighbors.size(); ++ n)
+        std::sort (neighbors.begin(), neighbors.end(),
+                   Sort_by_distance_to_point (pquery));
+        for (std::size_t n = (t == 0 ? 0 : 1); n < neighbors.size(); ++ n)
           {
             FT sq_dist = CGAL::squared_distance (pquery, neighbors[n]->point());
-            sum_sq_distances += weight * sq_dist;
-            nb += weight;
-            
-            if (nb < 6) // do not consider values under 6
+
+            sum_sq_distances += m_weights[t] * sq_dist;
+            nb += m_weights[t];
+
+            if (nb < 6.) // do not consider values under 6
               continue;
 
             FT dist = std::sqrt (sum_sq_distances / nb)
@@ -295,11 +321,10 @@ public:
             if (dist < dist_min)
               {
                 dist_min = dist;
-                k = nb;
+                k = (std::size_t)nb;
                 d = sq_dist;
               }
           }
-        weight *= m_cluster_size;
       }
   }
 
@@ -323,13 +348,13 @@ template <typename SamplesInputIterator,
 >
 void
 estimate_local_k_neighbor_scales(
-  SamplesInputIterator first,
-  SamplesInputIterator beyond,
-  SamplesPointPMap samples_pmap,
-  QueriesInputIterator first_query,
-  QueriesInputIterator beyond_query,
-  QueriesPointPMap queries_pmap,
-  OutputIterator output,
+  SamplesInputIterator first, ///< iterator over the first input sample.
+  SamplesInputIterator beyond, ///< past-the-end iterator over the input samples.
+  SamplesPointPMap samples_pmap, ///< property map: value_type of InputIterator -> Point_3 or Point_2
+  QueriesInputIterator first_query, ///< iterator over the first point where scale must be estimated
+  QueriesInputIterator beyond_query, ///< past-the-end iterator over the points where scale must be estimated
+  QueriesPointPMap queries_pmap, ///< property map: value_type of InputIterator -> Point_3 or Point_2
+  OutputIterator output, ///< output iterator to store the computed scales
   const Kernel& /*kernel*/) ///< geometric traits.
 {
   typedef typename boost::property_traits<SamplesPointPMap>::value_type Point_d;
@@ -350,7 +375,7 @@ std::size_t
 estimate_global_k_neighbor_scale(
   InputIterator first,  ///< iterator over the first input point.
   InputIterator beyond, ///< past-the-end iterator over the input points.
-  PointPMap point_pmap, ///< property map: value_type of InputIterator -> Point_3
+  PointPMap point_pmap, ///< property map: value_type of InputIterator -> Point_3 or Point_2
   const Kernel& kernel) ///< geometric traits.
 {
   std::vector<std::size_t> scales;
@@ -372,13 +397,13 @@ template <typename SamplesInputIterator,
 >
 void
 estimate_local_range_scales(
-  SamplesInputIterator first,
-  SamplesInputIterator beyond,
-  SamplesPointPMap samples_pmap,
-  QueriesInputIterator first_query,
-  QueriesInputIterator beyond_query,
-  QueriesPointPMap queries_pmap,
-  OutputIterator output,
+  SamplesInputIterator first, ///< iterator over the first input sample.
+  SamplesInputIterator beyond, ///< past-the-end iterator over the input samples.
+  SamplesPointPMap samples_pmap, ///< property map: value_type of InputIterator -> Point_3 or Point_2
+  QueriesInputIterator first_query, ///< iterator over the first point where scale must be estimated
+  QueriesInputIterator beyond_query, ///< past-the-end iterator over the points where scale must be estimated
+  QueriesPointPMap queries_pmap, ///< property map: value_type of InputIterator -> Point_3 or Point_2
+  OutputIterator output, ///< output iterator to store the computed scales
   const Kernel& /*kernel*/) ///< geometric traits.
 {
   typedef typename boost::property_traits<SamplesPointPMap>::value_type Point_d;
@@ -400,7 +425,7 @@ typename Kernel::FT
 estimate_global_range_scale(
   InputIterator first,  ///< iterator over the first input point.
   InputIterator beyond, ///< past-the-end iterator over the input points.
-  PointPMap point_pmap, ///< property map: value_type of InputIterator -> Point_3
+  PointPMap point_pmap, ///< property map: value_type of InputIterator -> Point_3 or Point_3
   const Kernel& kernel) ///< geometric traits.
 {
   std::vector<typename Kernel::FT> scales;

From e99d74aefc72750424aa0ca54f21d54f120fb781 Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Mon, 8 Aug 2016 10:30:33 +0200
Subject: [PATCH 05/26] Minor corrections

---
 Point_set_processing_3/include/CGAL/estimate_scale.h | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/Point_set_processing_3/include/CGAL/estimate_scale.h b/Point_set_processing_3/include/CGAL/estimate_scale.h
index 49b808e3064..2b9d06f9b85 100644
--- a/Point_set_processing_3/include/CGAL/estimate_scale.h
+++ b/Point_set_processing_3/include/CGAL/estimate_scale.h
@@ -81,7 +81,8 @@ public:
       {
         m_trees.push_back (new Tree());
         first_unused
-          = CGAL::hierarchy_simplify_point_set (first, first_unused, point_pmap, cluster_size, 1./3.);
+          = CGAL::hierarchy_simplify_point_set (first, first_unused, point_pmap,
+                                                m_cluster_size, 1./3.);
 
         for (InputIterator it = first; it != first_unused; ++ it)
           m_trees.back()->insert (get(point_pmap, *it));
@@ -124,7 +125,7 @@ public:
 
     FT dist_min = std::numeric_limits<FT>::max();
     FT sum_sq_distances = 0.;
-    std::size_t nb = 0;
+    FT nb = 0.;
     
     for (std::size_t t = 0; t < m_trees.size(); ++ t)
       {
@@ -142,7 +143,7 @@ public:
             sum_sq_distances += m_weights[t] * it->second;
             nb += m_weights[t];
 
-            if (nb < 6) // do not consider values under 6
+            if (nb < 6.) // do not consider values under 6
               continue;
 
             
@@ -290,7 +291,7 @@ public:
 
     FT dist_min = std::numeric_limits<FT>::max();
     FT sum_sq_distances = 0.;
-    FT nb = 0;
+    FT nb = 0.;
 
     const typename Kernel::Point_2& pquery = get(point_pmap, *query);
     for (std::size_t t = 0; t < m_point_sets.size(); ++ t)

From c1fc95dc9f4b9d900aca2dce05612c1d32d4c82d Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Mon, 8 Aug 2016 16:52:21 +0200
Subject: [PATCH 06/26] Code optimization

---
 .../include/CGAL/estimate_scale.h             | 56 ++++++++++++-------
 1 file changed, 37 insertions(+), 19 deletions(-)

diff --git a/Point_set_processing_3/include/CGAL/estimate_scale.h b/Point_set_processing_3/include/CGAL/estimate_scale.h
index 2b9d06f9b85..2500f300faf 100644
--- a/Point_set_processing_3/include/CGAL/estimate_scale.h
+++ b/Point_set_processing_3/include/CGAL/estimate_scale.h
@@ -42,9 +42,18 @@ class Quick_multiscale_approximate_knn_distance<Kernel, typename Kernel::Point_3
   typedef typename Neighbor_search::Tree Tree;
   typedef typename Neighbor_search::iterator Iterator;
 
+  template <typename ValueType, typename PointPMap>
+  struct Pmap_unary_function : public std::unary_function<ValueType, typename Kernel::Point_3>
+  {
+    PointPMap point_pmap;
+    Pmap_unary_function (PointPMap point_pmap) : point_pmap (point_pmap) { }
+    const typename Kernel::Point_3& operator() (const ValueType& v) const { return get(point_pmap, v); }
+  };
+  
   std::size_t m_cluster_size;
   std::vector<Tree*> m_trees;
   std::vector<FT> m_weights;
+  std::vector<FT> m_precomputed_factor;
 
 public:
 
@@ -55,14 +64,13 @@ public:
                                              std::size_t cluster_size = 10)
     : m_cluster_size (cluster_size)
   {
-    m_trees.push_back (new Tree ());
+    typedef Pmap_unary_function<typename std::iterator_traits<InputIterator>::value_type,
+                                PointPMap> Unary_f;
+
+    m_trees.push_back (new Tree (boost::make_transform_iterator (first, Unary_f(point_pmap)),
+                                 boost::make_transform_iterator (beyond, Unary_f(point_pmap))));
     m_weights.push_back (1.);
-    std::size_t nb_pts = 0;
-    for (InputIterator it = first; it != beyond; ++ it)
-      {
-        m_trees[0]->insert (get(point_pmap, *it));
-        ++ nb_pts;
-      }
+    std::size_t nb_pts = m_trees[0]->size();
     
     std::size_t nb_trees = 0;
     while (nb_pts > m_cluster_size)
@@ -75,17 +83,16 @@ public:
     m_weights.reserve (nb_trees);
 
     InputIterator first_unused = beyond;
-    std::cerr << nb_trees << std::endl;
+
     nb_pts = m_trees[0]->size();
     for (std::size_t i = 1; i < nb_trees; ++ i)
       {
-        m_trees.push_back (new Tree());
         first_unused
           = CGAL::hierarchy_simplify_point_set (first, first_unused, point_pmap,
                                                 m_cluster_size, 1./3.);
 
-        for (InputIterator it = first; it != first_unused; ++ it)
-          m_trees.back()->insert (get(point_pmap, *it));
+        m_trees.push_back (new Tree(boost::make_transform_iterator (first, Unary_f(point_pmap)),
+                                    boost::make_transform_iterator (first_unused, Unary_f(point_pmap))));
 
         m_weights.push_back (m_trees[0]->size() / (FT)(m_trees.back()->size()));
       }
@@ -115,10 +122,17 @@ public:
     return out;
   }
 
+  FT precomputed_factor (std::size_t index, FT nb)
+  {
+    if (m_precomputed_factor.size() == index)
+      m_precomputed_factor.push_back (0.91666666 * std::log (nb));
+    return m_precomputed_factor[index];
+  }
+  
+  
   template <typename InputIterator, typename PointPMap>
   void compute_scale (InputIterator query, PointPMap point_pmap,
-                      std::size_t& k, FT& d,
-                      std::ostream& log)
+                      std::size_t& k, FT& d)
   {
     k = 0;
     d = 0.;
@@ -126,7 +140,7 @@ public:
     FT dist_min = std::numeric_limits<FT>::max();
     FT sum_sq_distances = 0.;
     FT nb = 0.;
-    
+    std::size_t index = 0;
     for (std::size_t t = 0; t < m_trees.size(); ++ t)
       {
         Neighbor_search search (*(m_trees[t]), get(point_pmap, *query),
@@ -146,10 +160,14 @@ public:
             if (nb < 6.) // do not consider values under 6
               continue;
 
+            // FT dist = std::sqrt (sum_sq_distances / nb)
+            //   / std::pow (nb, 0.41666); // nb^(5/12)
+            
+            // sqrt(sum_sq_distances / nb) / nb^(5/12)
+            // Computed in log space for time optimization
+            FT dist = 0.5 * std::log (sum_sq_distances) - precomputed_factor (index, nb);
+            ++ index;
             
-            FT dist = std::sqrt (sum_sq_distances / nb)
-              / std::pow (nb, 0.41666); // nb^(5/12)
-            log << nb << " " << dist << std::endl;
             if (dist < dist_min)
               {
                 dist_min = dist;
@@ -242,7 +260,7 @@ public:
 
     InputIterator first_unused = beyond;
     nb_pts = m_point_sets[0]->number_of_vertices();
-    std::cerr << nb_trees << std::endl;
+
     for (std::size_t i = 1; i < nb_trees; ++ i)
       {
         first_unused
@@ -254,7 +272,7 @@ public:
 
         m_weights.push_back (nb_pts / (FT)(m_point_sets.back()->number_of_vertices()));
       }
-    std::cerr << std::endl;
+
     m_cluster_size = cluster_size;
   }
 

From e9464de8f28dc3c1024b137589a5407ce65ef587 Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Mon, 8 Aug 2016 17:29:17 +0200
Subject: [PATCH 07/26] A bit more optimizations

---
 .../include/CGAL/estimate_scale.h             | 75 +++++++++++++------
 1 file changed, 54 insertions(+), 21 deletions(-)

diff --git a/Point_set_processing_3/include/CGAL/estimate_scale.h b/Point_set_processing_3/include/CGAL/estimate_scale.h
index 2500f300faf..de7fea399de 100644
--- a/Point_set_processing_3/include/CGAL/estimate_scale.h
+++ b/Point_set_processing_3/include/CGAL/estimate_scale.h
@@ -122,11 +122,22 @@ public:
     return out;
   }
 
-  FT precomputed_factor (std::size_t index, FT nb)
+  void precompute_factors ()
   {
-    if (m_precomputed_factor.size() == index)
-      m_precomputed_factor.push_back (0.91666666 * std::log (nb));
-    return m_precomputed_factor[index];
+    FT nb = 0.;
+    for (std::size_t t = 0; t < m_trees.size(); ++ t)
+      {
+        std::size_t size = (t == (m_trees.size() - 1)
+                            ? m_trees[t]->size()
+                            : m_weights[t+1] / m_weights[t]);
+        for (std::size_t i = (t == 0 ? 0 : 1); i < size; ++ i)
+          {
+            nb += m_weights[t];
+            if (nb < 6.) // do not consider values under 6
+              continue;
+            m_precomputed_factor.push_back (0.91666666 * std::log (nb));
+          }
+      }
   }
   
   
@@ -134,6 +145,9 @@ public:
   void compute_scale (InputIterator query, PointPMap point_pmap,
                       std::size_t& k, FT& d)
   {
+    if (m_precomputed_factor.empty())
+      precompute_factors();
+    
     k = 0;
     d = 0.;
 
@@ -159,19 +173,15 @@ public:
 
             if (nb < 6.) // do not consider values under 6
               continue;
-
-            // FT dist = std::sqrt (sum_sq_distances / nb)
-            //   / std::pow (nb, 0.41666); // nb^(5/12)
             
             // sqrt(sum_sq_distances / nb) / nb^(5/12)
-            // Computed in log space for time optimization
-            FT dist = 0.5 * std::log (sum_sq_distances) - precomputed_factor (index, nb);
-            ++ index;
+            // Computed in log space with precomputed factor for time optimization
+            FT dist = 0.5 * std::log (sum_sq_distances) - m_precomputed_factor[index ++];
             
             if (dist < dist_min)
               {
                 dist_min = dist;
-                k = nb;
+                k = (std::size_t)nb;
                 d = it->second;
               }
           }
@@ -230,6 +240,7 @@ class Quick_multiscale_approximate_knn_distance<Kernel, typename Kernel::Point_2
   std::size_t m_cluster_size;
   std::vector<Point_set*> m_point_sets;
   std::vector<FT> m_weights;
+  std::vector<FT> m_precomputed_factor;
   
 public:
 
@@ -300,27 +311,48 @@ public:
     return out;
   }
 
+  void precompute_factors ()
+  {
+    FT nb = 0.;
+    for (std::size_t t = 0; t < m_point_sets.size(); ++ t)
+      {
+        std::size_t size = (t == m_point_sets.size() - 1
+                            ? m_point_sets[t]->number_of_vertices()
+                            : m_weights[t+1] / m_weights[t]);
+        for (std::size_t i = (t == 0 ? 0 : 1); i < size; ++ i)
+          {
+            nb += m_weights[t];
+            if (nb < 6.) // do not consider values under 6
+              continue;
+            m_precomputed_factor.push_back (1.25 * std::log (nb));
+          }
+      }
+  }
+  
   template <typename InputIterator, typename PointPMap>
   void compute_scale (InputIterator query, PointPMap point_pmap,
                       std::size_t& k, FT& d)
   {
+    if (m_precomputed_factor.empty())
+      precompute_factors();
+
     k = 0;
     d = 0.;
 
     FT dist_min = std::numeric_limits<FT>::max();
     FT sum_sq_distances = 0.;
     FT nb = 0.;
-
+    std::size_t index = 0;
+    
     const typename Kernel::Point_2& pquery = get(point_pmap, *query);
     for (std::size_t t = 0; t < m_point_sets.size(); ++ t)
       {
+        std::size_t size = ((t == m_point_sets.size() - 1)
+                            ? m_point_sets[t]->number_of_vertices()
+                            : m_weights[t+1] / m_weights[t]);
         std::vector<Vertex_handle> neighbors;
-        if (t == m_point_sets.size() - 1)
-          m_point_sets[t]->nearest_neighbors (pquery, m_point_sets[t]->number_of_vertices(),
-                                              std::back_inserter (neighbors));
-        else
-          m_point_sets[t]->nearest_neighbors (pquery, m_weights[t+1] / m_weights[t],
-                                              std::back_inserter (neighbors));
+        neighbors.reserve (size);
+        m_point_sets[t]->nearest_neighbors (pquery, size, std::back_inserter (neighbors));
 
         std::sort (neighbors.begin(), neighbors.end(),
                    Sort_by_distance_to_point (pquery));
@@ -334,9 +366,10 @@ public:
             if (nb < 6.) // do not consider values under 6
               continue;
 
-            FT dist = std::sqrt (sum_sq_distances / nb)
-              / std::pow (nb, 0.75); // nb^(3/4)
-
+            // sqrt(sum_sq_distances / nb) / nb^(3/4)
+            // Computed in log space with precomputed factor for time optimization
+            FT dist = 0.5 * std::log (sum_sq_distances) - m_precomputed_factor[index ++];
+            
             if (dist < dist_min)
               {
                 dist_min = dist;

From 0de675ff3b0197b970c9e6af43fc182d2c946701 Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Tue, 9 Aug 2016 11:37:45 +0200
Subject: [PATCH 08/26] Add doxygen doc and useful overloads

---
 .../include/CGAL/estimate_scale.h             | 245 +++++++++++++++++-
 1 file changed, 243 insertions(+), 2 deletions(-)

diff --git a/Point_set_processing_3/include/CGAL/estimate_scale.h b/Point_set_processing_3/include/CGAL/estimate_scale.h
index de7fea399de..79e76204cb9 100644
--- a/Point_set_processing_3/include/CGAL/estimate_scale.h
+++ b/Point_set_processing_3/include/CGAL/estimate_scale.h
@@ -61,7 +61,7 @@ public:
   Quick_multiscale_approximate_knn_distance (InputIterator first,
                                              InputIterator beyond,
                                              PointPMap point_pmap,
-                                             std::size_t cluster_size = 10)
+                                             std::size_t cluster_size = 25)
     : m_cluster_size (cluster_size)
   {
     typedef Pmap_unary_function<typename std::iterator_traits<InputIterator>::value_type,
@@ -248,7 +248,7 @@ public:
   Quick_multiscale_approximate_knn_distance (InputIterator first,
                                              InputIterator beyond,
                                              PointPMap point_pmap,
-                                             std::size_t cluster_size = 20)
+                                             std::size_t cluster_size = 25)
     : m_cluster_size (cluster_size)
   {
     typedef Pmap_unary_function<typename std::iterator_traits<InputIterator>::value_type,
@@ -391,6 +391,32 @@ public:
 // Public section
 // ----------------------------------------------------------------------------
 
+/// \ingroup PkgPointSetProcessing
+
+/// Estimates the local scale in a K nearest neighbors sense on a set
+/// of user-defined locations. The computed scales correspond to the
+/// smallest scales such that the K subsets of points have the
+/// appearance of a surface in 3D or the appearance of a curve in 2D.
+///
+///
+/// @tparam SamplesInputIterator iterator over input sample points.
+/// @tparam SamplesPointPMap is a model of `ReadablePropertyMap` with
+///        value type `Point_3<Kernel>` or `Point_2<Kernel>`.  It can
+///        be omitted if the value type of `SamplesInputIterator` is
+///        convertible to `Point_3<Kernel>` or to `Point_2<Kernel>`.
+/// @tparam QueriesInputIterator iterator over points where scale
+///        should be computed.
+/// @tparam QueriesInputIterator is a model of `ReadablePropertyMap`
+///        with value type `Point_3<Kernel>` or `Point_2<Kernel>`.  It
+///        can be omitted if the value type of `QueriesInputIterator` is
+///        convertible to `Point_3<Kernel>` or to `Point_2<Kernel>`.
+/// @tparam Kernel Geometric traits class.  It can be omitted and
+///        deduced automatically from the value type of `SamplesPointPMap`.
+///
+/// @note This function accepts both 2D and 3D points, but sample
+///      points and query must have the same dimension.
+
+// This variant requires all parameters.
 template <typename SamplesInputIterator,
           typename SamplesPointPMap,
           typename QueriesInputIterator,
@@ -419,6 +445,26 @@ estimate_local_k_neighbor_scales(
     *(output ++) = kdtree.compute_k_scale (it, queries_pmap);
 }
 
+  
+/// \ingroup PkgPointSetProcessing
+
+/// Estimates the gloabl scale in a K nearest neighbors sense. The
+/// computed scale corresponds to the smallest scale such that the K
+/// subsets of points have the appearance of a surface in 3D or the
+/// appearance of a curve in 2D.
+///
+///
+/// @tparam InputIterator iterator over input points.
+/// @tparam PointPMap is a model of `ReadablePropertyMap` with
+///        value type `Point_3<Kernel>` or `Point_2<Kernel>`.  It can
+///        be omitted if the value type of `InputIterator` is
+///        convertible to `Point_3<Kernel>` or to `Point_2<Kernel>`.
+/// @tparam Kernel Geometric traits class.  It can be omitted and
+///        deduced automatically from the value type of `PointPMap`.
+///
+/// @note This function accepts both 2D and 3D points.
+
+// This variant requires all parameters.
 template <typename InputIterator,
           typename PointPMap,
           typename Kernel
@@ -440,6 +486,33 @@ estimate_global_k_neighbor_scale(
 }
 
   
+/// \ingroup PkgPointSetProcessing
+
+/// Estimates the local scale in a range sense on a set of
+/// user-defined locations. The computed scales correspond to the
+/// smallest scales such that the subsets of points included in the
+/// sphere range have the appearance of a surface in 3D or the
+/// appearance of a curve in 2D.
+///
+///
+/// @tparam SamplesInputIterator iterator over input sample points.
+/// @tparam SamplesPointPMap is a model of `ReadablePropertyMap` with
+///        value type `Point_3<Kernel>` or `Point_2<Kernel>`.  It can
+///        be omitted if the value type of `SamplesInputIterator` is
+///        convertible to `Point_3<Kernel>` or to `Point_2<Kernel>`.
+/// @tparam QueriesInputIterator iterator over points where scale
+///        should be computed.
+/// @tparam QueriesInputIterator is a model of `ReadablePropertyMap`
+///        with value type `Point_3<Kernel>` or `Point_2<Kernel>`.  It
+///        can be omitted if the value type of `QueriesInputIterator` is
+///        convertible to `Point_3<Kernel>` or to `Point_2<Kernel>`.
+/// @tparam Kernel Geometric traits class.  It can be omitted and
+///        deduced automatically from the value type of `SamplesPointPMap`.
+///
+/// @note This function accepts both 2D and 3D points, but sample
+///      points and query must have the same dimension.
+
+// This variant requires all parameters.
 template <typename SamplesInputIterator,
           typename SamplesPointPMap,
           typename QueriesInputIterator,
@@ -469,6 +542,26 @@ estimate_local_range_scales(
 
 }
 
+  
+/// \ingroup PkgPointSetProcessing
+
+/// Estimates the gloabl scale in a range sense. The computed scale
+/// corresponds to the smallest scale such that the subsets of points
+/// inside the sphere range have the appearance of a surface in 3D or
+/// the appearance of a curve in 2D.
+///
+///
+/// @tparam InputIterator iterator over input points.
+/// @tparam PointPMap is a model of `ReadablePropertyMap` with
+///        value type `Point_3<Kernel>` or `Point_2<Kernel>`.  It can
+///        be omitted if the value type of `InputIterator` is
+///        convertible to `Point_3<Kernel>` or to `Point_2<Kernel>`.
+/// @tparam Kernel Geometric traits class.  It can be omitted and
+///        deduced automatically from the value type of `PointPMap`.
+///
+/// @note This function accepts both 2D and 3D points.
+
+// This variant requires all parameters.
 template <typename InputIterator,
           typename PointPMap,
           typename Kernel
@@ -490,6 +583,154 @@ estimate_global_range_scale(
 }
 
 
+// ----------------------------------------------------------------------------
+// Useful overloads
+// ----------------------------------------------------------------------------
+/// \cond SKIP_IN_MANUAL
+
+template <typename SamplesInputIterator,
+          typename SamplesPointPMap,
+          typename QueriesInputIterator,
+          typename QueriesPointPMap,
+          typename OutputIterator
+>
+void
+estimate_local_k_neighbor_scales(
+  SamplesInputIterator first, ///< iterator over the first input sample.
+  SamplesInputIterator beyond, ///< past-the-end iterator over the input samples.
+  SamplesPointPMap samples_pmap, ///< property map: value_type of InputIterator -> Point_3 or Point_2
+  QueriesInputIterator first_query, ///< iterator over the first point where scale must be estimated
+  QueriesInputIterator beyond_query, ///< past-the-end iterator over the points where scale must be estimated
+  QueriesPointPMap queries_pmap, ///< property map: value_type of InputIterator -> Point_3 or Point_2
+  OutputIterator output) ///< output iterator to store the computed scales
+{
+  typedef typename boost::property_traits<SamplesPointPMap>::value_type Point;
+  typedef typename Kernel_traits<Point>::Kernel Kernel;
+  estimate_local_k_neighbor_scales (first, beyond, samples_pmap, first_query, beyond_query, queries_pmap,
+                                    output, Kernel());
+}
+
+template <typename SamplesInputIterator,
+          typename SamplesPointPMap,
+          typename QueriesInputIterator,
+          typename QueriesPointPMap,
+          typename OutputIterator
+>
+void
+estimate_local_k_neighbor_scales(
+  SamplesInputIterator first, ///< iterator over the first input sample.
+  SamplesInputIterator beyond, ///< past-the-end iterator over the input samples.
+  QueriesInputIterator first_query, ///< iterator over the first point where scale must be estimated
+  QueriesInputIterator beyond_query, ///< past-the-end iterator over the points where scale must be estimated
+  OutputIterator output) ///< output iterator to store the computed scales
+{
+  estimate_local_k_neighbor_scales
+    (first, beyond,
+     make_identity_property_map (typename std::iterator_traits<SamplesInputIterator>::value_type()),
+     first_query, beyond_query,
+     make_identity_property_map (typename std::iterator_traits<QueriesInputIterator>::value_type()),
+     output);
+ }
+
+
+template <typename InputIterator,
+          typename PointPMap
+>
+std::size_t
+estimate_global_k_neighbor_scale(
+  InputIterator first,  ///< iterator over the first input point.
+  InputIterator beyond, ///< past-the-end iterator over the input points.
+  PointPMap point_pmap) ///< property map: value_type of InputIterator -> Point_3 or Point_2
+{
+  typedef typename boost::property_traits<PointPMap>::value_type Point;
+  typedef typename Kernel_traits<Point>::Kernel Kernel;
+  return estimate_global_k_neighbor_scale (first, beyond, point_pmap, Kernel());
+}
+
+template <typename InputIterator
+>
+std::size_t
+estimate_global_k_neighbor_scale(
+  InputIterator first,  ///< iterator over the first input point.
+  InputIterator beyond) ///< past-the-end iterator over the input points.
+{
+  return estimate_global_k_neighbor_scale
+    (first, beyond, make_identity_property_map (typename std::iterator_traits<InputIterator>::value_type()));
+}
+
+
+template <typename SamplesInputIterator,
+          typename SamplesPointPMap,
+          typename QueriesInputIterator,
+          typename QueriesPointPMap,
+          typename OutputIterator
+>
+void
+estimate_local_range_scales(
+  SamplesInputIterator first, ///< iterator over the first input sample.
+  SamplesInputIterator beyond, ///< past-the-end iterator over the input samples.
+  SamplesPointPMap samples_pmap, ///< property map: value_type of InputIterator -> Point_3 or Point_2
+  QueriesInputIterator first_query, ///< iterator over the first point where scale must be estimated
+  QueriesInputIterator beyond_query, ///< past-the-end iterator over the points where scale must be estimated
+  QueriesPointPMap queries_pmap, ///< property map: value_type of InputIterator -> Point_3 or Point_2
+  OutputIterator output) ///< output iterator to store the computed scales
+{
+  typedef typename boost::property_traits<SamplesPointPMap>::value_type Point;
+  typedef typename Kernel_traits<Point>::Kernel Kernel;
+  estimate_local_range_scales(first, beyond, samples_pmap, first_query, beyond_query, queries_pmap,
+                              output, Kernel());
+}
+
+
+template <typename SamplesInputIterator,
+          typename QueriesInputIterator,
+          typename OutputIterator
+>
+void
+estimate_local_range_scales(
+  SamplesInputIterator first, ///< iterator over the first input sample.
+  SamplesInputIterator beyond, ///< past-the-end iterator over the input samples.
+  QueriesInputIterator first_query, ///< iterator over the first point where scale must be estimated
+  QueriesInputIterator beyond_query, ///< past-the-end iterator over the points where scale must be estimated
+  OutputIterator output) ///< output iterator to store the computed scales
+{
+  estimate_local_range_scales
+    (first, beyond,
+     make_identity_property_map (typename std::iterator_traits<SamplesInputIterator>::value_type()),
+     first_query, beyond_query,
+     make_identity_property_map (typename std::iterator_traits<QueriesInputIterator>::value_type()),
+     output);
+}
+
+
+
+template <typename InputIterator,
+          typename PointPMap
+>
+typename Kernel_traits<typename boost::property_traits<PointPMap>::value_type>::Kernel::FT
+estimate_global_range_scale(
+  InputIterator first,  ///< iterator over the first input point.
+  InputIterator beyond, ///< past-the-end iterator over the input points.
+  PointPMap point_pmap) ///< property map: value_type of InputIterator -> Point_3 or Point_3
+{
+  typedef typename boost::property_traits<PointPMap>::value_type Point;
+  typedef typename Kernel_traits<Point>::Kernel Kernel;
+  return estimate_global_range_scale (first, beyond, point_pmap, Kernel());
+}
+
+
+
+template <typename InputIterator>
+typename Kernel_traits<typename std::iterator_traits<InputIterator>::value_type>::Kernel::FT
+estimate_global_range_scale(
+  InputIterator first,  ///< iterator over the first input point.
+  InputIterator beyond) ///< past-the-end iterator over the input points.
+{
+  return estimate_global_range_scale
+    (first, beyond, make_identity_property_map (typename std::iterator_traits<InputIterator>::value_type()));
+                                      
+}
+/// \endcond  
 
 } //namespace CGAL
 

From 195e4083dc90c8f76ef2c51162a4d436961cdedb Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Wed, 10 Aug 2016 16:10:22 +0200
Subject: [PATCH 09/26] Update changes.html

---
 Installation/changes.html | 11 +++++++++++
 1 file changed, 11 insertions(+)

diff --git a/Installation/changes.html b/Installation/changes.html
index bae95e7af8f..45fb4affeaf 100644
--- a/Installation/changes.html
+++ b/Installation/changes.html
@@ -127,6 +127,17 @@ and <code>src/</code> directories).
 </table>
 
 <HR>
+<h2 id="release4.10">Release 4.10 </h2>
+<div>
+  <p>Release date: </p>
+  <h3>Point Set Processing</h3>
+  <ul>
+    <li>New functions for automatic scale estimations: either a global
+    scale or multiple local scales can be estimated for both 2D and 3D
+    point sets based on the assumption that they sample a curve in 2D
+      or a surface in 3D.</li>
+  </ul>
+</div>
 <h2 id="release4.9">Release 4.9 </h2>
 <div>
   <p>Release date: Sept 2016 </p>

From 0d9eda15aaa62bafb3a570fe58d5b31820217232 Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Wed, 10 Aug 2016 16:10:52 +0200
Subject: [PATCH 10/26] User manual on scale estimation

---
 .../PackageDescription.txt                    |  4 ++
 .../Point_set_processing_3.txt                | 37 ++++++++++++++++++-
 .../doc/Point_set_processing_3/examples.txt   |  1 +
 3 files changed, 41 insertions(+), 1 deletion(-)

diff --git a/Point_set_processing_3/doc/Point_set_processing_3/PackageDescription.txt b/Point_set_processing_3/doc/Point_set_processing_3/PackageDescription.txt
index 0daeb79d209..3af1268711a 100644
--- a/Point_set_processing_3/doc/Point_set_processing_3/PackageDescription.txt
+++ b/Point_set_processing_3/doc/Point_set_processing_3/PackageDescription.txt
@@ -36,6 +36,10 @@
 - `CGAL::read_ply_custom_points()`
 - `CGAL::read_xyz_points()`
 - `CGAL::compute_average_spacing()`
+- `CGAL::estimate_global_k_neighbor_scale()`
+- `CGAL::estimate_global_range_scale()`
+- `CGAL::estimate_local_k_neighbor_scales()`
+- `CGAL::estimate_local_range_scales()`
 - `CGAL::remove_outliers()`
 - `CGAL::grid_simplify_point_set()`
 - `CGAL::random_simplify_point_set()`
diff --git a/Point_set_processing_3/doc/Point_set_processing_3/Point_set_processing_3.txt b/Point_set_processing_3/doc/Point_set_processing_3/Point_set_processing_3.txt
index e4ba85448c4..86b9c362c7c 100644
--- a/Point_set_processing_3/doc/Point_set_processing_3/Point_set_processing_3.txt
+++ b/Point_set_processing_3/doc/Point_set_processing_3/Point_set_processing_3.txt
@@ -124,7 +124,7 @@ these attributes in user-defined containers.
 
 \cgalExample{Point_set_processing_3/read_ply_points_with_colors_example.cpp}
 
-\section Point_set_processing_3Analysis Analysis
+\section Point_set_processing_3Spacing Average Spacing
 
 Function `compute_average_spacing()` computes the average
 spacing of all input points to their `k` nearest neighbor points,
@@ -147,6 +147,41 @@ found in other \cgal components:
 - `bounding_box()`
 - `Min_sphere_of_spheres_d<Traits>`
 
+\section Point_set_processing_3Scale Automatic Scale Estimation
+
+Point sets are often used to sample objects with a higher dimension,
+typically a curve in 2D or a surface in 3D. In such cases, finding the
+scale of the objet is crucial, that is to say finding the minimal
+number of points (or the minimal local range) such that the subset of
+points has the appearance of a curve in 2D or a surface in 3D.
+
+\cgal provides 2 functions that automatically estimate the scale of a
+2D point set sampling a curve or a 3D point set sampling a surface:
+
+- `estimate_global_k_neighbor_scale()`
+- `estimate_global_range_scale()`
+
+In some specific cases, the scale of a point set might not be
+homogeneous (for example if the point set contains variable
+noise). \cgal also provides 2 functions that automatically estimate
+the scales of a point set at a set of user-defined locations:
+
+- `estimate_local_k_neighbor_scales()`
+- `estimate_local_range_scales()`
+
+Notice that the 4 functions presented here work both with 2D points
+and 3D points and that they shouldn't be used if the point sets do not
+sample a curve in 2D or a surface in 3D.
+
+\subsection Point_set_processing_3Example_scale_estimation Example
+
+The following example reads a point set in the `xyz` format and
+estimates the global scale (both in K neighbor sense and in a range
+sense). If no `xyz` point set can be read, it tries to read a 2D point
+set and perform the same scales estimation.
+
+\cgalExample{Point_set_processing_3/scale_estimation_example.cpp}
+
 
 
 \section Point_set_processing_3OutlierRemoval Outlier Removal
diff --git a/Point_set_processing_3/doc/Point_set_processing_3/examples.txt b/Point_set_processing_3/doc/Point_set_processing_3/examples.txt
index c3996a09dfe..bfc21ed1d29 100644
--- a/Point_set_processing_3/doc/Point_set_processing_3/examples.txt
+++ b/Point_set_processing_3/doc/Point_set_processing_3/examples.txt
@@ -2,6 +2,7 @@
 \example Point_set_processing_3/read_write_xyz_point_set_example.cpp
 \example Point_set_processing_3/read_ply_points_with_colors_example.cpp
 \example Point_set_processing_3/average_spacing_example.cpp
+\example Point_set_processing_3/scale_estimation_example.cpp
 \example Point_set_processing_3/remove_outliers_example.cpp
 \example Point_set_processing_3/grid_simplification_example.cpp
 \example Point_set_processing_3/grid_simplify_indices.cpp

From ba1812642b6706be6954dfdc75a4d9fdfa359b9d Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Wed, 10 Aug 2016 16:12:36 +0200
Subject: [PATCH 11/26] Update example

---
 .../scale_estimation_example.cpp              | 109 ++++++++++--------
 1 file changed, 63 insertions(+), 46 deletions(-)

diff --git a/Point_set_processing_3/examples/Point_set_processing_3/scale_estimation_example.cpp b/Point_set_processing_3/examples/Point_set_processing_3/scale_estimation_example.cpp
index d10028701ba..7a0c5fed97e 100644
--- a/Point_set_processing_3/examples/Point_set_processing_3/scale_estimation_example.cpp
+++ b/Point_set_processing_3/examples/Point_set_processing_3/scale_estimation_example.cpp
@@ -1,6 +1,8 @@
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 #include <CGAL/estimate_scale.h>
 #include <CGAL/IO/read_xyz_points.h>
+#include <CGAL/Timer.h>
+#include <CGAL/Memory_sizer.h>
 
 #include <vector>
 #include <fstream>
@@ -15,59 +17,74 @@ typedef Kernel::Point_2 Point_2;
 int main (int argc, char** argv)
 {
 
-  // 3D CASE
-  {
-    const char* fname = (argc>1)?argv[1]:"data/sphere_20k.xyz";
-    // Reads a .xyz point set file in points.
-    // As the point is the second element of the tuple (that is with index 1)
-    // we use a property map that accesses the 1st element of the tuple.
+  const char* fname = (argc>1)?argv[1]:"data/sphere_20k.xyz";
   
-    std::vector<Point_3> points;
-    std::ifstream stream(fname);
-    if (!stream ||
-        !CGAL::read_xyz_points(
-                               stream, std::back_inserter(points)))
-      {
-        std::cerr << "Error: cannot read file " << fname << std::endl;
-        return EXIT_FAILURE;
-      }
+  CGAL::Timer task_timer; 
+  
+  std::ifstream stream(fname);
+  if (stream)
+    {
+      std::vector<Point_3> points;
+      std::cerr << "TRYING 3D CASE" << std::endl;
+      if (CGAL::read_xyz_points(stream, std::back_inserter(points)))
+        {
+          task_timer.start();
+          
+          // estimate global scale
+          std::size_t k_scale = CGAL::estimate_global_k_neighbor_scale (points.begin(), points.end());
+          FT range_scale = CGAL::estimate_global_range_scale (points.begin(), points.end());
 
-    // estimate global scale
-    std::size_t k_scale = CGAL::estimate_global_k_neighbor_scale (points.begin(), points.end(),
-                                                                  CGAL::Identity_property_map<Point_3>(),
-                                                                  Kernel());
-    std::cout << "Global K scale: " << k_scale << std::endl;
+          std::size_t memory = CGAL::Memory_sizer().virtual_size();
+          double time = task_timer.time();
 
-    FT range_scale = CGAL::estimate_global_range_scale (points.begin(), points.end(),
-                                                        CGAL::Identity_property_map<Point_3>(),
-                                                        Kernel());
-    std::cout << "Global range scale: " << range_scale << std::endl;
-  }
+          std::cout << "Scales computed in " << time << " second(s) using "
+                    << (memory>>20) << " MiB of memory:" << std::endl;
+          std::cout << " * Global K scale: " << k_scale << std::endl;
+          std::cout << " * Global range scale: " << range_scale << std::endl;
+        }
+      else
+        {
+          stream.seekg(0);
+          std::cerr << "TRYING 2D CASE" << std::endl;
+          std::vector<Point_2> points_2;
 
-  // 2D CASE
-  {
-    std::vector<Point_2> points;
-    for (std::size_t i = 0; i < 5000; ++ i)
-      {
-        double angle = CGAL_PI * 2. * (rand() / (double)RAND_MAX);
-        double radius = 1. + 0.1 * (rand() / (double)RAND_MAX);
-        points.push_back (Point_2 (radius * std::cos (angle),
-                                   radius * std::sin (angle)));
-      }
+          // Read ASCII 2D point set file
+          std::string str;
+          while (getline (stream, str))
+            {
+              std::istringstream iss (str);
+              double x = 0., y = 0.;
+              iss >> x >> y;
+              points_2.push_back (Point_2 (x, y));
+            }
 
-        
-    // estimate global scale
-    std::size_t k_scale = CGAL::estimate_global_k_neighbor_scale (points.begin(), points.end(),
-                                                                  CGAL::Identity_property_map<Point_2>(),
-                                                                  Kernel());
-    std::cout << "Global K scale: " << k_scale << std::endl;
+          if (points_2.empty())
+            {
+              std::cerr << "Error: cannot read file " << fname << std::endl;
+              return EXIT_FAILURE;
+            }
+          
+          task_timer.start();
 
-    FT range_scale = CGAL::estimate_global_range_scale (points.begin(), points.end(),
-                                                        CGAL::Identity_property_map<Point_2>(),
-                                                        Kernel());
-    std::cout << "Global range scale: " << range_scale << std::endl;
+          // estimate global scale
+          std::size_t k_scale = CGAL::estimate_global_k_neighbor_scale (points_2.begin(), points_2.end());
+          FT range_scale = CGAL::estimate_global_range_scale (points_2.begin(), points_2.end());
+
+          std::size_t memory = CGAL::Memory_sizer().virtual_size();
+          double time = task_timer.time();
+
+          std::cout << "Scales computed in " << time << " second(s) using "
+                    << (memory>>20) << " MiB of memory:" << std::endl;
+          std::cout << " * Global K scale: " << k_scale << std::endl;
+          std::cout << " * Global range scale: " << range_scale << std::endl;
+        }
+    }
+  else
+    {
+      std::cerr << "Error: cannot read file " << fname << std::endl;
+      return EXIT_FAILURE;
+    }
 
-  }
   return EXIT_SUCCESS;
 }
 

From 08e4f7e6c83032a745bdf0f33a349ab2c7f2582a Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Wed, 10 Aug 2016 16:13:10 +0200
Subject: [PATCH 12/26] Minor corrections in doxygen doc

---
 Point_set_processing_3/include/CGAL/estimate_scale.h | 10 ++++++----
 1 file changed, 6 insertions(+), 4 deletions(-)

diff --git a/Point_set_processing_3/include/CGAL/estimate_scale.h b/Point_set_processing_3/include/CGAL/estimate_scale.h
index 79e76204cb9..f8f74a24c63 100644
--- a/Point_set_processing_3/include/CGAL/estimate_scale.h
+++ b/Point_set_processing_3/include/CGAL/estimate_scale.h
@@ -448,7 +448,7 @@ estimate_local_k_neighbor_scales(
   
 /// \ingroup PkgPointSetProcessing
 
-/// Estimates the gloabl scale in a K nearest neighbors sense. The
+/// Estimates the global scale in a K nearest neighbors sense. The
 /// computed scale corresponds to the smallest scale such that the K
 /// subsets of points have the appearance of a surface in 3D or the
 /// appearance of a curve in 2D.
@@ -463,7 +463,8 @@ estimate_local_k_neighbor_scales(
 ///        deduced automatically from the value type of `PointPMap`.
 ///
 /// @note This function accepts both 2D and 3D points.
-
+///
+/// @return The estimated scale in the K nearest neighbors sense.
 // This variant requires all parameters.
 template <typename InputIterator,
           typename PointPMap,
@@ -545,7 +546,7 @@ estimate_local_range_scales(
   
 /// \ingroup PkgPointSetProcessing
 
-/// Estimates the gloabl scale in a range sense. The computed scale
+/// Estimates the global scale in a range sense. The computed scale
 /// corresponds to the smallest scale such that the subsets of points
 /// inside the sphere range have the appearance of a surface in 3D or
 /// the appearance of a curve in 2D.
@@ -560,7 +561,8 @@ estimate_local_range_scales(
 ///        deduced automatically from the value type of `PointPMap`.
 ///
 /// @note This function accepts both 2D and 3D points.
-
+///
+/// @return The estimated scale in the range sense.
 // This variant requires all parameters.
 template <typename InputIterator,
           typename PointPMap,

From fee10522bb07603180ac8a9575e3de93bb4c57b5 Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Wed, 10 Aug 2016 17:29:53 +0200
Subject: [PATCH 13/26] Add biblio reference

---
 Documentation/doc/biblio/cgal_manual.bib              | 11 +++++++++++
 .../Point_set_processing_3/Point_set_processing_3.txt |  3 ++-
 2 files changed, 13 insertions(+), 1 deletion(-)

diff --git a/Documentation/doc/biblio/cgal_manual.bib b/Documentation/doc/biblio/cgal_manual.bib
index 7e461f528be..5ef4ed06080 100644
--- a/Documentation/doc/biblio/cgal_manual.bib
+++ b/Documentation/doc/biblio/cgal_manual.bib
@@ -715,6 +715,17 @@ note         = {\url{ttp://hal.inria.fr/inria-00090522}}
   ,pages        = "325--338"
 }
 
+@article{cgal:gcsa-nasr-13,
+   journal = {Computer Graphics Forum},
+   title = {{Noise-Adaptive Shape Reconstruction from Raw Point Sets}},
+   author = {Simon Giraudot and David Cohen-Steiner and Pierre Alliez },
+   pages = {229-238},
+   volume= {32},
+   number= {5},
+   year = {2013},
+   DOI = {10.1111/cgf.12189},
+}
+
 @manual{ cgal:g-gmpal-96
   ,author       = "T. Granlund"
   ,title        = "{GNU MP}, The {GNU} Multiple Precision Arithmetic Library,
diff --git a/Point_set_processing_3/doc/Point_set_processing_3/Point_set_processing_3.txt b/Point_set_processing_3/doc/Point_set_processing_3/Point_set_processing_3.txt
index 86b9c362c7c..528f2f9f7aa 100644
--- a/Point_set_processing_3/doc/Point_set_processing_3/Point_set_processing_3.txt
+++ b/Point_set_processing_3/doc/Point_set_processing_3/Point_set_processing_3.txt
@@ -153,7 +153,8 @@ Point sets are often used to sample objects with a higher dimension,
 typically a curve in 2D or a surface in 3D. In such cases, finding the
 scale of the objet is crucial, that is to say finding the minimal
 number of points (or the minimal local range) such that the subset of
-points has the appearance of a curve in 2D or a surface in 3D.
+points has the appearance of a curve in 2D or a surface in 3D
+\cgalCite{cgal:gcsa-nasr-13}.
 
 \cgal provides 2 functions that automatically estimate the scale of a
 2D point set sampling a curve or a 3D point set sampling a surface:

From 7c633e8acc0069ec4c38b9deacd154abcacfdc36 Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Thu, 18 Aug 2016 14:52:38 +0200
Subject: [PATCH 14/26] Remove useless and ambiguous template typenames

---
 Point_set_processing_3/include/CGAL/estimate_scale.h | 2 --
 1 file changed, 2 deletions(-)

diff --git a/Point_set_processing_3/include/CGAL/estimate_scale.h b/Point_set_processing_3/include/CGAL/estimate_scale.h
index f8f74a24c63..690e21c2ec9 100644
--- a/Point_set_processing_3/include/CGAL/estimate_scale.h
+++ b/Point_set_processing_3/include/CGAL/estimate_scale.h
@@ -613,9 +613,7 @@ estimate_local_k_neighbor_scales(
 }
 
 template <typename SamplesInputIterator,
-          typename SamplesPointPMap,
           typename QueriesInputIterator,
-          typename QueriesPointPMap,
           typename OutputIterator
 >
 void

From 307a64d45a62d8ead98961daf4f7935885cb3de6 Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Thu, 18 Aug 2016 14:53:45 +0200
Subject: [PATCH 15/26] Separate 2D and 3D examples and improve these examples

---
 .../Point_set_processing_3/CMakeLists.txt     |   1 +
 .../scale_estimation_2d_example.cpp           |  52 ++++++++
 .../scale_estimation_example.cpp              | 111 ++++++++----------
 3 files changed, 102 insertions(+), 62 deletions(-)
 create mode 100644 Point_set_processing_3/examples/Point_set_processing_3/scale_estimation_2d_example.cpp

diff --git a/Point_set_processing_3/examples/Point_set_processing_3/CMakeLists.txt b/Point_set_processing_3/examples/Point_set_processing_3/CMakeLists.txt
index e9e37fd0808..60d29d9a76d 100644
--- a/Point_set_processing_3/examples/Point_set_processing_3/CMakeLists.txt
+++ b/Point_set_processing_3/examples/Point_set_processing_3/CMakeLists.txt
@@ -54,6 +54,7 @@ if ( CGAL_FOUND )
   create_single_source_cgal_program( "read_ply_points_with_colors_example.cpp" )
   create_single_source_cgal_program( "remove_outliers_example.cpp" )
   create_single_source_cgal_program( "scale_estimation_example.cpp" )
+  create_single_source_cgal_program( "scale_estimation_2d_example.cpp" )
   create_single_source_cgal_program( "wlop_simplify_and_regularize_point_set_example.cpp" )
   create_single_source_cgal_program( "edge_aware_upsample_point_set_example.cpp" )
   
diff --git a/Point_set_processing_3/examples/Point_set_processing_3/scale_estimation_2d_example.cpp b/Point_set_processing_3/examples/Point_set_processing_3/scale_estimation_2d_example.cpp
new file mode 100644
index 00000000000..9f1599eb2a7
--- /dev/null
+++ b/Point_set_processing_3/examples/Point_set_processing_3/scale_estimation_2d_example.cpp
@@ -0,0 +1,52 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/estimate_scale.h>
+#include <CGAL/IO/read_xyz_points.h>
+#include <CGAL/Random.h>
+
+#include <vector>
+#include <fstream>
+#include <boost/tuple/tuple.hpp>
+
+// Types
+typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
+typedef Kernel::FT FT;
+typedef Kernel::Point_2 Point_2;
+
+int main (int, char**)
+{
+  std::vector<Point_2> samples;
+  samples.reserve (100000);
+
+  // Generate circle with gradually variable noise
+  //  - noise-free for points with x close to (-1)
+  //  - noisy for points with x close to (+1)
+  for (std::size_t i = 0; i < 100000; ++ i)
+    {
+      FT theta = CGAL::get_default_random().get_double(0, 2. * CGAL_PI);
+      FT noise = 0.5 * (std::cos(theta) + 1.) * CGAL::get_default_random().get_double(0., 0.2);
+      int mult = (CGAL::get_default_random().get_bool() ? 1 : -1);
+      samples.push_back (Point_2 (std::cos(theta) * (1. + mult * noise * noise),
+                                  std::sin(theta) * (1. + mult * noise * noise)));
+    }
+
+  // Search for local scales on 3 different locations
+  std::vector<Point_2> queries;
+  queries.reserve (3);
+  queries.push_back (Point_2 (-1., 0.));
+  queries.push_back (Point_2 (0., 1.));
+  queries.push_back (Point_2 (1., 0.));
+
+  std::vector<std::size_t> k_scales;
+  CGAL::estimate_local_k_neighbor_scales (samples.begin(), samples.end(),
+                                          queries.begin(), queries.end(),
+                                          std::back_inserter (k_scales));
+
+  // Display results
+  std::cerr << "K-Scales found:" << std::endl
+            << " - On noise-free region: " << k_scales[0] << std::endl // Should be small
+            << " - On moderately noisy region: " << k_scales[1] << std::endl // Should be higher
+            << " - On very noisy region: " << k_scales[2] << std::endl; // Should be even higher
+
+  return EXIT_SUCCESS;
+}
+
diff --git a/Point_set_processing_3/examples/Point_set_processing_3/scale_estimation_example.cpp b/Point_set_processing_3/examples/Point_set_processing_3/scale_estimation_example.cpp
index 7a0c5fed97e..438b452112b 100644
--- a/Point_set_processing_3/examples/Point_set_processing_3/scale_estimation_example.cpp
+++ b/Point_set_processing_3/examples/Point_set_processing_3/scale_estimation_example.cpp
@@ -1,6 +1,10 @@
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
-#include <CGAL/estimate_scale.h>
 #include <CGAL/IO/read_xyz_points.h>
+
+#include <CGAL/estimate_scale.h>
+#include <CGAL/jet_smooth_point_set.h>
+#include <CGAL/grid_simplify_point_set.h>
+
 #include <CGAL/Timer.h>
 #include <CGAL/Memory_sizer.h>
 
@@ -8,82 +12,65 @@
 #include <fstream>
 #include <boost/tuple/tuple.hpp>
 
+// Concurrency
+#ifdef CGAL_LINKED_WITH_TBB
+typedef CGAL::Parallel_tag Concurrency_tag;
+#else
+typedef CGAL::Sequential_tag Concurrency_tag;
+#endif
+
 // Types
 typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
 typedef Kernel::FT FT;
 typedef Kernel::Point_3 Point_3;
-typedef Kernel::Point_2 Point_2;
 
 int main (int argc, char** argv)
 {
 
   const char* fname = (argc>1)?argv[1]:"data/sphere_20k.xyz";
   
-  CGAL::Timer task_timer; 
+  CGAL::Timer task_timer;
   
+  std::vector<Point_3> points;  
   std::ifstream stream(fname);
-  if (stream)
+
+  // read input
+  if (!(stream
+        && CGAL::read_xyz_points(stream, std::back_inserter(points))))
     {
-      std::vector<Point_3> points;
-      std::cerr << "TRYING 3D CASE" << std::endl;
-      if (CGAL::read_xyz_points(stream, std::back_inserter(points)))
-        {
-          task_timer.start();
-          
-          // estimate global scale
-          std::size_t k_scale = CGAL::estimate_global_k_neighbor_scale (points.begin(), points.end());
-          FT range_scale = CGAL::estimate_global_range_scale (points.begin(), points.end());
-
-          std::size_t memory = CGAL::Memory_sizer().virtual_size();
-          double time = task_timer.time();
-
-          std::cout << "Scales computed in " << time << " second(s) using "
-                    << (memory>>20) << " MiB of memory:" << std::endl;
-          std::cout << " * Global K scale: " << k_scale << std::endl;
-          std::cout << " * Global range scale: " << range_scale << std::endl;
-        }
-      else
-        {
-          stream.seekg(0);
-          std::cerr << "TRYING 2D CASE" << std::endl;
-          std::vector<Point_2> points_2;
-
-          // Read ASCII 2D point set file
-          std::string str;
-          while (getline (stream, str))
-            {
-              std::istringstream iss (str);
-              double x = 0., y = 0.;
-              iss >> x >> y;
-              points_2.push_back (Point_2 (x, y));
-            }
-
-          if (points_2.empty())
-            {
-              std::cerr << "Error: cannot read file " << fname << std::endl;
-              return EXIT_FAILURE;
-            }
-          
-          task_timer.start();
-
-          // estimate global scale
-          std::size_t k_scale = CGAL::estimate_global_k_neighbor_scale (points_2.begin(), points_2.end());
-          FT range_scale = CGAL::estimate_global_range_scale (points_2.begin(), points_2.end());
-
-          std::size_t memory = CGAL::Memory_sizer().virtual_size();
-          double time = task_timer.time();
-
-          std::cout << "Scales computed in " << time << " second(s) using "
-                    << (memory>>20) << " MiB of memory:" << std::endl;
-          std::cout << " * Global K scale: " << k_scale << std::endl;
-          std::cout << " * Global range scale: " << range_scale << std::endl;
-        }
-    }
-  else
-    {
-      std::cerr << "Error: cannot read file " << fname << std::endl;
+      std::cerr << "Error: can't read input file" << std::endl;
       return EXIT_FAILURE;
     }
+          
+  // estimate k scale
+  task_timer.start();
+  std::size_t k_scale = CGAL::estimate_global_k_neighbor_scale (points.begin(), points.end());
+  task_timer.stop();
+  
+  // Example: use estimated k as scale for jet smoothing
+  CGAL::jet_smooth_point_set<Concurrency_tag>
+    (points.begin(), points.end(),
+     k_scale);
+
+  // estimate range scale
+  task_timer.start();
+  FT range_scale = CGAL::estimate_global_range_scale (points.begin(), points.end());
+  task_timer.stop();
+
+  // Example: use estimated range for grid simplification
+  points.erase (CGAL::grid_simplify_point_set (points.begin(), points.end(), range_scale),
+                points.end());
+
+  
+  // print some informations on runtime
+  std::size_t memory = CGAL::Memory_sizer().virtual_size();
+  double time = task_timer.time();
+
+  std::cout << "Scales computed in " << time << " second(s) using "
+            << (memory>>20) << " MiB of memory:" << std::endl;
+  std::cout << " * Global K scale: " << k_scale << std::endl;
+  std::cout << " * Global range scale: " << range_scale << std::endl;
+
 
   return EXIT_SUCCESS;
 }

From 1a56459f993dedf9be3bd166b2ef2da8a0b8f58f Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Fri, 19 Aug 2016 14:12:36 +0200
Subject: [PATCH 16/26] Update doc

---
 .../Point_set_processing_3.txt                | 30 ++++++++++++++-----
 .../doc/Point_set_processing_3/examples.txt   |  1 +
 2 files changed, 23 insertions(+), 8 deletions(-)

diff --git a/Point_set_processing_3/doc/Point_set_processing_3/Point_set_processing_3.txt b/Point_set_processing_3/doc/Point_set_processing_3/Point_set_processing_3.txt
index 528f2f9f7aa..239cadab87a 100644
--- a/Point_set_processing_3/doc/Point_set_processing_3/Point_set_processing_3.txt
+++ b/Point_set_processing_3/doc/Point_set_processing_3/Point_set_processing_3.txt
@@ -162,6 +162,12 @@ points has the appearance of a curve in 2D or a surface in 3D
 - `estimate_global_k_neighbor_scale()`
 - `estimate_global_range_scale()`
 
+Functions such as `jet_estimate_normals()` or
+`remove_outliers()` require a K neighbor scale while others such as
+`grid_simplify_point_set()` require a range
+scale. `vcm_estimate_normals()` is an example of a function that
+accepts both.
+
 In some specific cases, the scale of a point set might not be
 homogeneous (for example if the point set contains variable
 noise). \cgal also provides 2 functions that automatically estimate
@@ -170,19 +176,27 @@ the scales of a point set at a set of user-defined locations:
 - `estimate_local_k_neighbor_scales()`
 - `estimate_local_range_scales()`
 
-Notice that the 4 functions presented here work both with 2D points
-and 3D points and that they shouldn't be used if the point sets do not
-sample a curve in 2D or a surface in 3D.
+The 4 functions presented here work both with 2D points and 3D points
+and they shouldn't be used if the point sets do not sample a curve in
+2D or a surface in 3D.
 
-\subsection Point_set_processing_3Example_scale_estimation Example
+\subsection Point_set_processing_3Examples_scale_estimation Examples
 
-The following example reads a point set in the `xyz` format and
-estimates the global scale (both in K neighbor sense and in a range
-sense). If no `xyz` point set can be read, it tries to read a 2D point
-set and perform the same scales estimation.
+The following example reads a 3D point set in the `xyz` format and:
+
+- Estimates the K neighbor global scale
+- Uses it to smooth the point set
+- Estimates the range global scale
+- Uses it to simplify the point set
 
 \cgalExample{Point_set_processing_3/scale_estimation_example.cpp}
 
+This second example generates a 2D point set sampling a circle with
+variable noise. It then estimates the scale at 3 different locations
+of the domain.
+
+\cgalExample{Point_set_processing_3/scale_estimation_2d_example.cpp}
+
 
 
 \section Point_set_processing_3OutlierRemoval Outlier Removal
diff --git a/Point_set_processing_3/doc/Point_set_processing_3/examples.txt b/Point_set_processing_3/doc/Point_set_processing_3/examples.txt
index bfc21ed1d29..bd9a495d190 100644
--- a/Point_set_processing_3/doc/Point_set_processing_3/examples.txt
+++ b/Point_set_processing_3/doc/Point_set_processing_3/examples.txt
@@ -3,6 +3,7 @@
 \example Point_set_processing_3/read_ply_points_with_colors_example.cpp
 \example Point_set_processing_3/average_spacing_example.cpp
 \example Point_set_processing_3/scale_estimation_example.cpp
+\example Point_set_processing_3/scale_estimation_2d_example.cpp
 \example Point_set_processing_3/remove_outliers_example.cpp
 \example Point_set_processing_3/grid_simplification_example.cpp
 \example Point_set_processing_3/grid_simplify_indices.cpp

From 00381feae40e5f7ed5ac944af436b3789bdbc453 Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Thu, 25 Aug 2016 07:53:00 +0200
Subject: [PATCH 17/26] Return output iterators

---
 .../include/CGAL/estimate_scale.h             | 27 ++++++++++---------
 1 file changed, 15 insertions(+), 12 deletions(-)

diff --git a/Point_set_processing_3/include/CGAL/estimate_scale.h b/Point_set_processing_3/include/CGAL/estimate_scale.h
index 690e21c2ec9..fe79bf1b3e3 100644
--- a/Point_set_processing_3/include/CGAL/estimate_scale.h
+++ b/Point_set_processing_3/include/CGAL/estimate_scale.h
@@ -424,7 +424,7 @@ template <typename SamplesInputIterator,
           typename OutputIterator,
           typename Kernel
 >
-void
+OutputIterator
 estimate_local_k_neighbor_scales(
   SamplesInputIterator first, ///< iterator over the first input sample.
   SamplesInputIterator beyond, ///< past-the-end iterator over the input samples.
@@ -443,6 +443,8 @@ estimate_local_k_neighbor_scales(
   // Compute local scales everywhere
   for (QueriesInputIterator it = first_query; it != beyond_query; ++ it)
     *(output ++) = kdtree.compute_k_scale (it, queries_pmap);
+
+  return output;
 }
 
   
@@ -521,7 +523,7 @@ template <typename SamplesInputIterator,
           typename OutputIterator,
           typename Kernel
 >
-void
+OutputIterator
 estimate_local_range_scales(
   SamplesInputIterator first, ///< iterator over the first input sample.
   SamplesInputIterator beyond, ///< past-the-end iterator over the input samples.
@@ -541,6 +543,7 @@ estimate_local_range_scales(
   for (QueriesInputIterator it = first_query; it != beyond_query; ++ it)
     *(output ++) = kdtree.compute_range_scale (it, queries_pmap);
 
+  return output;
 }
 
   
@@ -596,7 +599,7 @@ template <typename SamplesInputIterator,
           typename QueriesPointPMap,
           typename OutputIterator
 >
-void
+OutputIterator
 estimate_local_k_neighbor_scales(
   SamplesInputIterator first, ///< iterator over the first input sample.
   SamplesInputIterator beyond, ///< past-the-end iterator over the input samples.
@@ -608,15 +611,15 @@ estimate_local_k_neighbor_scales(
 {
   typedef typename boost::property_traits<SamplesPointPMap>::value_type Point;
   typedef typename Kernel_traits<Point>::Kernel Kernel;
-  estimate_local_k_neighbor_scales (first, beyond, samples_pmap, first_query, beyond_query, queries_pmap,
-                                    output, Kernel());
+  return estimate_local_k_neighbor_scales (first, beyond, samples_pmap, first_query, beyond_query,
+                                           queries_pmap, output, Kernel());
 }
 
 template <typename SamplesInputIterator,
           typename QueriesInputIterator,
           typename OutputIterator
 >
-void
+OutputIterator
 estimate_local_k_neighbor_scales(
   SamplesInputIterator first, ///< iterator over the first input sample.
   SamplesInputIterator beyond, ///< past-the-end iterator over the input samples.
@@ -624,7 +627,7 @@ estimate_local_k_neighbor_scales(
   QueriesInputIterator beyond_query, ///< past-the-end iterator over the points where scale must be estimated
   OutputIterator output) ///< output iterator to store the computed scales
 {
-  estimate_local_k_neighbor_scales
+  return estimate_local_k_neighbor_scales
     (first, beyond,
      make_identity_property_map (typename std::iterator_traits<SamplesInputIterator>::value_type()),
      first_query, beyond_query,
@@ -665,7 +668,7 @@ template <typename SamplesInputIterator,
           typename QueriesPointPMap,
           typename OutputIterator
 >
-void
+OutputIterator
 estimate_local_range_scales(
   SamplesInputIterator first, ///< iterator over the first input sample.
   SamplesInputIterator beyond, ///< past-the-end iterator over the input samples.
@@ -677,8 +680,8 @@ estimate_local_range_scales(
 {
   typedef typename boost::property_traits<SamplesPointPMap>::value_type Point;
   typedef typename Kernel_traits<Point>::Kernel Kernel;
-  estimate_local_range_scales(first, beyond, samples_pmap, first_query, beyond_query, queries_pmap,
-                              output, Kernel());
+  return estimate_local_range_scales(first, beyond, samples_pmap, first_query, beyond_query,
+                                     queries_pmap, output, Kernel());
 }
 
 
@@ -686,7 +689,7 @@ template <typename SamplesInputIterator,
           typename QueriesInputIterator,
           typename OutputIterator
 >
-void
+OutputIterator
 estimate_local_range_scales(
   SamplesInputIterator first, ///< iterator over the first input sample.
   SamplesInputIterator beyond, ///< past-the-end iterator over the input samples.
@@ -694,7 +697,7 @@ estimate_local_range_scales(
   QueriesInputIterator beyond_query, ///< past-the-end iterator over the points where scale must be estimated
   OutputIterator output) ///< output iterator to store the computed scales
 {
-  estimate_local_range_scales
+  return estimate_local_range_scales
     (first, beyond,
      make_identity_property_map (typename std::iterator_traits<SamplesInputIterator>::value_type()),
      first_query, beyond_query,

From d94a7a63167f447e0996cfa45e8083763fe2e81a Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Thu, 25 Aug 2016 08:06:24 +0200
Subject: [PATCH 18/26] Unify: location -> query

---
 .../doc/Point_set_processing_3/Point_set_processing_3.txt   | 6 +++---
 Point_set_processing_3/include/CGAL/estimate_scale.h        | 6 +++---
 2 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/Point_set_processing_3/doc/Point_set_processing_3/Point_set_processing_3.txt b/Point_set_processing_3/doc/Point_set_processing_3/Point_set_processing_3.txt
index 239cadab87a..8fe27ca17c7 100644
--- a/Point_set_processing_3/doc/Point_set_processing_3/Point_set_processing_3.txt
+++ b/Point_set_processing_3/doc/Point_set_processing_3/Point_set_processing_3.txt
@@ -171,7 +171,7 @@ accepts both.
 In some specific cases, the scale of a point set might not be
 homogeneous (for example if the point set contains variable
 noise). \cgal also provides 2 functions that automatically estimate
-the scales of a point set at a set of user-defined locations:
+the scales of a point set at a set of user-defined query points:
 
 - `estimate_local_k_neighbor_scales()`
 - `estimate_local_range_scales()`
@@ -192,8 +192,8 @@ The following example reads a 3D point set in the `xyz` format and:
 \cgalExample{Point_set_processing_3/scale_estimation_example.cpp}
 
 This second example generates a 2D point set sampling a circle with
-variable noise. It then estimates the scale at 3 different locations
-of the domain.
+variable noise. It then estimates the scale at 3 different query
+points in the domain.
 
 \cgalExample{Point_set_processing_3/scale_estimation_2d_example.cpp}
 
diff --git a/Point_set_processing_3/include/CGAL/estimate_scale.h b/Point_set_processing_3/include/CGAL/estimate_scale.h
index fe79bf1b3e3..b1e04d4533c 100644
--- a/Point_set_processing_3/include/CGAL/estimate_scale.h
+++ b/Point_set_processing_3/include/CGAL/estimate_scale.h
@@ -394,8 +394,8 @@ public:
 /// \ingroup PkgPointSetProcessing
 
 /// Estimates the local scale in a K nearest neighbors sense on a set
-/// of user-defined locations. The computed scales correspond to the
-/// smallest scales such that the K subsets of points have the
+/// of user-defined query points. The computed scales correspond to
+/// the smallest scales such that the K subsets of points have the
 /// appearance of a surface in 3D or the appearance of a curve in 2D.
 ///
 ///
@@ -492,7 +492,7 @@ estimate_global_k_neighbor_scale(
 /// \ingroup PkgPointSetProcessing
 
 /// Estimates the local scale in a range sense on a set of
-/// user-defined locations. The computed scales correspond to the
+/// user-defined query points. The computed scales correspond to the
 /// smallest scales such that the subsets of points included in the
 /// sphere range have the appearance of a surface in 3D or the
 /// appearance of a curve in 2D.

From cb7914f1ec912e7a2a6d8d84d142b5ce85691d3c Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Thu, 25 Aug 2016 08:22:58 +0200
Subject: [PATCH 19/26] Document output iterators

---
 Point_set_processing_3/include/CGAL/estimate_scale.h | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/Point_set_processing_3/include/CGAL/estimate_scale.h b/Point_set_processing_3/include/CGAL/estimate_scale.h
index b1e04d4533c..a008293ef0b 100644
--- a/Point_set_processing_3/include/CGAL/estimate_scale.h
+++ b/Point_set_processing_3/include/CGAL/estimate_scale.h
@@ -410,6 +410,8 @@ public:
 ///        with value type `Point_3<Kernel>` or `Point_2<Kernel>`.  It
 ///        can be omitted if the value type of `QueriesInputIterator` is
 ///        convertible to `Point_3<Kernel>` or to `Point_2<Kernel>`.
+/// @tparam OutputIterator is used to store the computed scales. Its
+///        value type is `std::size_t`.
 /// @tparam Kernel Geometric traits class.  It can be omitted and
 ///        deduced automatically from the value type of `SamplesPointPMap`.
 ///
@@ -509,6 +511,8 @@ estimate_global_k_neighbor_scale(
 ///        with value type `Point_3<Kernel>` or `Point_2<Kernel>`.  It
 ///        can be omitted if the value type of `QueriesInputIterator` is
 ///        convertible to `Point_3<Kernel>` or to `Point_2<Kernel>`.
+/// @tparam OutputIterator is used to store the computed scales. Its
+///        value type is `Kernel::FT`.
 /// @tparam Kernel Geometric traits class.  It can be omitted and
 ///        deduced automatically from the value type of `SamplesPointPMap`.
 ///

From a3e4a8b927fc7c883d2de8837af979f74ead8ac6 Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Thu, 25 Aug 2016 09:05:02 +0200
Subject: [PATCH 20/26] Terminology correction (value type)

---
 Point_set_processing_3/include/CGAL/estimate_scale.h | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/Point_set_processing_3/include/CGAL/estimate_scale.h b/Point_set_processing_3/include/CGAL/estimate_scale.h
index a008293ef0b..509ac036c17 100644
--- a/Point_set_processing_3/include/CGAL/estimate_scale.h
+++ b/Point_set_processing_3/include/CGAL/estimate_scale.h
@@ -410,8 +410,8 @@ public:
 ///        with value type `Point_3<Kernel>` or `Point_2<Kernel>`.  It
 ///        can be omitted if the value type of `QueriesInputIterator` is
 ///        convertible to `Point_3<Kernel>` or to `Point_2<Kernel>`.
-/// @tparam OutputIterator is used to store the computed scales. Its
-///        value type is `std::size_t`.
+/// @tparam OutputIterator is used to store the computed scales. It accepts
+///        values of type `std::size_t`.
 /// @tparam Kernel Geometric traits class.  It can be omitted and
 ///        deduced automatically from the value type of `SamplesPointPMap`.
 ///
@@ -511,8 +511,8 @@ estimate_global_k_neighbor_scale(
 ///        with value type `Point_3<Kernel>` or `Point_2<Kernel>`.  It
 ///        can be omitted if the value type of `QueriesInputIterator` is
 ///        convertible to `Point_3<Kernel>` or to `Point_2<Kernel>`.
-/// @tparam OutputIterator is used to store the computed scales. Its
-///        value type is `Kernel::FT`.
+/// @tparam OutputIterator is used to store the computed scales. It accepts
+///        values of type `Kernel::FT`.
 /// @tparam Kernel Geometric traits class.  It can be omitted and
 ///        deduced automatically from the value type of `SamplesPointPMap`.
 ///

From ea24e3ccdeb89bca72b29772dce0d5a1408cf38e Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Thu, 25 Aug 2016 09:31:20 +0200
Subject: [PATCH 21/26] Add reference to user manual in ref manual

---
 Point_set_processing_3/include/CGAL/estimate_scale.h | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/Point_set_processing_3/include/CGAL/estimate_scale.h b/Point_set_processing_3/include/CGAL/estimate_scale.h
index 509ac036c17..996363bd203 100644
--- a/Point_set_processing_3/include/CGAL/estimate_scale.h
+++ b/Point_set_processing_3/include/CGAL/estimate_scale.h
@@ -396,7 +396,8 @@ public:
 /// Estimates the local scale in a K nearest neighbors sense on a set
 /// of user-defined query points. The computed scales correspond to
 /// the smallest scales such that the K subsets of points have the
-/// appearance of a surface in 3D or the appearance of a curve in 2D.
+/// appearance of a surface in 3D or the appearance of a curve in 2D
+/// (see \ref Point_set_processing_3Scale).
 ///
 ///
 /// @tparam SamplesInputIterator iterator over input sample points.
@@ -455,7 +456,7 @@ estimate_local_k_neighbor_scales(
 /// Estimates the global scale in a K nearest neighbors sense. The
 /// computed scale corresponds to the smallest scale such that the K
 /// subsets of points have the appearance of a surface in 3D or the
-/// appearance of a curve in 2D.
+/// appearance of a curve in 2D (see \ref Point_set_processing_3Scale).
 ///
 ///
 /// @tparam InputIterator iterator over input points.
@@ -497,7 +498,7 @@ estimate_global_k_neighbor_scale(
 /// user-defined query points. The computed scales correspond to the
 /// smallest scales such that the subsets of points included in the
 /// sphere range have the appearance of a surface in 3D or the
-/// appearance of a curve in 2D.
+/// appearance of a curve in 2D (see \ref Point_set_processing_3Scale).
 ///
 ///
 /// @tparam SamplesInputIterator iterator over input sample points.
@@ -556,7 +557,7 @@ estimate_local_range_scales(
 /// Estimates the global scale in a range sense. The computed scale
 /// corresponds to the smallest scale such that the subsets of points
 /// inside the sphere range have the appearance of a surface in 3D or
-/// the appearance of a curve in 2D.
+/// the appearance of a curve in 2D (see \ref Point_set_processing_3Scale).
 ///
 ///
 /// @tparam InputIterator iterator over input points.

From 8e67eb7dd2d0f46e9a017ef2f78e8f579d7d545d Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Fri, 14 Oct 2016 09:12:21 +0200
Subject: [PATCH 22/26] Separate examples in 2 sections

---
 .../doc/Point_set_processing_3/Point_set_processing_3.txt     | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/Point_set_processing_3/doc/Point_set_processing_3/Point_set_processing_3.txt b/Point_set_processing_3/doc/Point_set_processing_3/Point_set_processing_3.txt
index 8fe27ca17c7..b6ed3e5d1db 100644
--- a/Point_set_processing_3/doc/Point_set_processing_3/Point_set_processing_3.txt
+++ b/Point_set_processing_3/doc/Point_set_processing_3/Point_set_processing_3.txt
@@ -180,7 +180,7 @@ The 4 functions presented here work both with 2D points and 3D points
 and they shouldn't be used if the point sets do not sample a curve in
 2D or a surface in 3D.
 
-\subsection Point_set_processing_3Examples_scale_estimation Examples
+\subsection Point_set_processing_3Example_scale_estimation_global Global Scale Example 
 
 The following example reads a 3D point set in the `xyz` format and:
 
@@ -191,6 +191,8 @@ The following example reads a 3D point set in the `xyz` format and:
 
 \cgalExample{Point_set_processing_3/scale_estimation_example.cpp}
 
+\subsection Point_set_processing_3Example_scale_estimation_local Local Scales Example
+
 This second example generates a 2D point set sampling a circle with
 variable noise. It then estimates the scale at 3 different query
 points in the domain.

From 2a37790e475d089154fa08c65eb4997c76a72a17 Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Wed, 21 Dec 2016 13:50:18 +0100
Subject: [PATCH 23/26] Add license header

---
 .../include/CGAL/estimate_scale.h             | 20 +++++++++++++++++++
 1 file changed, 20 insertions(+)

diff --git a/Point_set_processing_3/include/CGAL/estimate_scale.h b/Point_set_processing_3/include/CGAL/estimate_scale.h
index 996363bd203..48a62882228 100644
--- a/Point_set_processing_3/include/CGAL/estimate_scale.h
+++ b/Point_set_processing_3/include/CGAL/estimate_scale.h
@@ -1,3 +1,23 @@
+// Copyright (c) 2013 INRIA Sophia-Antipolis (France).
+// Copyright (c) 2016 GeometryFactory Sarl (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+// You can redistribute it and/or modify it under the terms of the GNU
+// General Public License as published by the Free Software Foundation,
+// either version 3 of the License, or (at your option) any later version.
+//
+// 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) : Simon Giraudot
+
 #ifndef CGAL_ESTIMATE_SCALE_H
 #define CGAL_ESTIMATE_SCALE_H
 

From d423e5e6f8b1651d308d2dc6a02a9eedcd6eff9e Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Thu, 29 Dec 2016 09:28:36 +0100
Subject: [PATCH 24/26] Fix warning with explicit casts

---
 .../scale_estimation_example.cpp                     |  2 +-
 Point_set_processing_3/include/CGAL/estimate_scale.h | 12 ++++++------
 2 files changed, 7 insertions(+), 7 deletions(-)

diff --git a/Point_set_processing_3/examples/Point_set_processing_3/scale_estimation_example.cpp b/Point_set_processing_3/examples/Point_set_processing_3/scale_estimation_example.cpp
index 438b452112b..bac9900fecf 100644
--- a/Point_set_processing_3/examples/Point_set_processing_3/scale_estimation_example.cpp
+++ b/Point_set_processing_3/examples/Point_set_processing_3/scale_estimation_example.cpp
@@ -50,7 +50,7 @@ int main (int argc, char** argv)
   // Example: use estimated k as scale for jet smoothing
   CGAL::jet_smooth_point_set<Concurrency_tag>
     (points.begin(), points.end(),
-     k_scale);
+     static_cast<unsigned int>(k_scale));
 
   // estimate range scale
   task_timer.start();
diff --git a/Point_set_processing_3/include/CGAL/estimate_scale.h b/Point_set_processing_3/include/CGAL/estimate_scale.h
index 48a62882228..8b55821ed49 100644
--- a/Point_set_processing_3/include/CGAL/estimate_scale.h
+++ b/Point_set_processing_3/include/CGAL/estimate_scale.h
@@ -109,7 +109,7 @@ public:
       {
         first_unused
           = CGAL::hierarchy_simplify_point_set (first, first_unused, point_pmap,
-                                                m_cluster_size, 1./3.);
+                                                static_cast<unsigned int>(m_cluster_size), 1./3.);
 
         m_trees.push_back (new Tree(boost::make_transform_iterator (first, Unary_f(point_pmap)),
                                     boost::make_transform_iterator (first_unused, Unary_f(point_pmap))));
@@ -149,7 +149,7 @@ public:
       {
         std::size_t size = (t == (m_trees.size() - 1)
                             ? m_trees[t]->size()
-                            : m_weights[t+1] / m_weights[t]);
+                            : static_cast<std::size_t>(m_weights[t+1] / m_weights[t]));
         for (std::size_t i = (t == 0 ? 0 : 1); i < size; ++ i)
           {
             nb += m_weights[t];
@@ -180,7 +180,7 @@ public:
         Neighbor_search search (*(m_trees[t]), get(point_pmap, *query),
                                 (t == (m_trees.size() - 1)
                                  ? m_trees[t]->size()
-                                 : m_weights[t+1] / m_weights[t]));
+                                 : static_cast<unsigned int>(m_weights[t+1] / m_weights[t])));
         Iterator it = search.begin();
         
         if (t != 0) // Skip first point except on first scale
@@ -296,7 +296,7 @@ public:
       {
         first_unused
           = CGAL::hierarchy_simplify_point_set (first, first_unused, Pmap_to_3d<PointPMap> (point_pmap),
-                                                m_cluster_size, 1./3.);
+                                                static_cast<unsigned int>(m_cluster_size), 1./3.);
 
         m_point_sets.push_back (new Point_set (boost::make_transform_iterator (first, Unary_f(point_pmap)),
                                                boost::make_transform_iterator (first_unused, Unary_f(point_pmap))));
@@ -338,7 +338,7 @@ public:
       {
         std::size_t size = (t == m_point_sets.size() - 1
                             ? m_point_sets[t]->number_of_vertices()
-                            : m_weights[t+1] / m_weights[t]);
+                            : static_cast<std::size_t>(m_weights[t+1] / m_weights[t]));
         for (std::size_t i = (t == 0 ? 0 : 1); i < size; ++ i)
           {
             nb += m_weights[t];
@@ -369,7 +369,7 @@ public:
       {
         std::size_t size = ((t == m_point_sets.size() - 1)
                             ? m_point_sets[t]->number_of_vertices()
-                            : m_weights[t+1] / m_weights[t]);
+                            : static_cast<std::size_t>(m_weights[t+1] / m_weights[t]));
         std::vector<Vertex_handle> neighbors;
         neighbors.reserve (size);
         m_point_sets[t]->nearest_neighbors (pquery, size, std::back_inserter (neighbors));

From 5edb6e1459556ebb56e445f2cf04ef30476ca7a4 Mon Sep 17 00:00:00 2001
From: Simon Giraudot <simon.giraudot@geometryfactory.com>
Date: Fri, 30 Dec 2016 07:38:44 +0100
Subject: [PATCH 25/26] Warning fix: explicit conversion from size_t to
 unsigned int

---
 Point_set_processing_3/include/CGAL/estimate_scale.h | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/Point_set_processing_3/include/CGAL/estimate_scale.h b/Point_set_processing_3/include/CGAL/estimate_scale.h
index 8b55821ed49..8e01da0a342 100644
--- a/Point_set_processing_3/include/CGAL/estimate_scale.h
+++ b/Point_set_processing_3/include/CGAL/estimate_scale.h
@@ -178,9 +178,9 @@ public:
     for (std::size_t t = 0; t < m_trees.size(); ++ t)
       {
         Neighbor_search search (*(m_trees[t]), get(point_pmap, *query),
-                                (t == (m_trees.size() - 1)
-                                 ? m_trees[t]->size()
-                                 : static_cast<unsigned int>(m_weights[t+1] / m_weights[t])));
+                                static_cast<unsigned int>((t == (m_trees.size() - 1)
+                                                           ? m_trees[t]->size()
+                                                           : m_weights[t+1] / m_weights[t])));
         Iterator it = search.begin();
         
         if (t != 0) // Skip first point except on first scale

From ab784644703259a63e3b3326bce6e53909e5281f Mon Sep 17 00:00:00 2001
From: Andreas Fabri <andreas.fabri@geometryfactory.com>
Date: Fri, 30 Dec 2016 10:59:07 +0100
Subject: [PATCH 26/26] Fix VC++ min/max problem

---
 Point_set_processing_3/include/CGAL/estimate_scale.h | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/Point_set_processing_3/include/CGAL/estimate_scale.h b/Point_set_processing_3/include/CGAL/estimate_scale.h
index 8e01da0a342..1ac7897e0f5 100644
--- a/Point_set_processing_3/include/CGAL/estimate_scale.h
+++ b/Point_set_processing_3/include/CGAL/estimate_scale.h
@@ -171,7 +171,7 @@ public:
     k = 0;
     d = 0.;
 
-    FT dist_min = std::numeric_limits<FT>::max();
+    FT dist_min = (std::numeric_limits<FT>::max)();
     FT sum_sq_distances = 0.;
     FT nb = 0.;
     std::size_t index = 0;
@@ -359,7 +359,7 @@ public:
     k = 0;
     d = 0.;
 
-    FT dist_min = std::numeric_limits<FT>::max();
+    FT dist_min = (std::numeric_limits<FT>::max)();
     FT sum_sq_distances = 0.;
     FT nb = 0.;
     std::size_t index = 0;