add function to simplify a 3D polyline (closed or open) with approximation guarantees

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/simplify_polyline.h

@@ -0,0 +1,195 @@
+// Copyright (c) 2020 GeometryFactory (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
+//
+//
+// Author(s)     : Sébastien Loriot
+
+#ifndef CGAL_POLYGON_MESH_PROCESSING_SIMPLIFY_POLYLINE_H
+#define CGAL_POLYGON_MESH_PROCESSING_SIMPLIFY_POLYLINE_H
+
+#include <CGAL/license/Polygon_mesh_processing/repair.h>
+
+#include <CGAL/boost/graph/Named_function_parameters.h>
+#include <CGAL/boost/graph/named_params_helper.h>
+
+#include <type_traits>
+
+
+namespace CGAL {
+namespace Polygon_mesh_processing {
+namespace experimental {
+
+enum Polyline_simplification_algorithms { DOUGLAS_PEUCKER, ITERATIVE };
+
+template <typename PointRangeIn, typename PointRangeOut,
+          typename NamedParametersIn, typename NamedParametersOut>
+void simplify_polyline(const PointRangeIn& input,
+                             PointRangeOut& output,
+                       const double max_squared_frechet_distance,
+                       const NamedParametersIn& np_in,
+                       const NamedParametersOut& np_out)
+{
+  using parameters::choose_parameter;
+  using parameters::get_parameter;
+
+  static_assert(std::is_same<
+    typename std::iterator_traits<typename PointRangeIn::const_iterator>::value_type,
+    typename std::iterator_traits<typename PointRangeOut::const_iterator>::value_type >::value );
+
+  typedef typename GetPointMap<PointRangeIn, NamedParametersIn>::type Point_map_in;
+  typedef typename GetPointMap<PointRangeOut, NamedParametersOut>::type Point_map_out;
+  typedef typename Point_set_processing_3::GetK<PointRangeIn, NamedParametersIn>::Kernel Kernel;
+
+  Point_map_in in_pm = choose_parameter<Point_map_in>(get_parameter(np_in, internal_np::point_map));
+  Point_map_out out_pm = choose_parameter<Point_map_out>(get_parameter(np_out, internal_np::point_map));
+
+  const Polyline_simplification_algorithms algorithm =
+    choose_parameter(get_parameter(np_in, internal_np::algorithm), DOUGLAS_PEUCKER);
+
+  switch(algorithm)
+  {
+    case ITERATIVE:
+    {
+      // skip points in the input range that do not contains any information
+      if (input.size()<=2)
+      {
+        output.reserve(input.size());
+        for (const auto& p : input)
+        {
+          output.push_back(p);
+          put(out_pm, output.back(), get(in_pm, p));
+        }
+        return;
+      }
+
+      auto is_valid_approx = [&input, &in_pm, max_squared_frechet_distance](
+        std::size_t b, std::size_t e,
+        const typename Kernel::Line_3& line)
+      {
+        typename Kernel::Compare_squared_distance_3 compare_squared_distance;
+        for (std::size_t i=b+1; i<e; ++i)
+        {
+          if (compare_squared_distance(get(in_pm, input[i]), line, max_squared_frechet_distance) == LARGER)
+            return false;
+        }
+        return true;
+      };
+
+      std::size_t bi=0;
+      while(bi!=input.size())
+      {
+        std::size_t ei=bi+2;
+        output.push_back(input[bi]);
+        put(out_pm, output.back(), get(in_pm, input[bi]));
+
+        while(ei<input.size())
+        {
+          typename Kernel::Line_3 sl(get(in_pm,input[bi]), get(in_pm,input[ei]));
+          if (is_valid_approx(bi,ei,sl))
+            ++ei; // we skip ei-1
+          else
+          {
+            bi=ei-1; // ei-1 shall not be skipt
+            break;
+          }
+        }
+        if(ei>=input.size()) break;
+      }
+      output.push_back(input.back());
+      put(out_pm, output.back(), get(in_pm, input.back()));
+      return;
+    }
+    case DOUGLAS_PEUCKER:
+    {
+      const bool is_closed = input.front()==input.back();
+
+      std::size_t nb_points = is_closed ? input.size()-1 : input.size();
+
+      if (nb_points<=2)
+      {
+        output.reserve(input.size());
+        for (const auto& p : input)
+        {
+          output.push_back(p);
+          put(out_pm, output.back(), get(in_pm, p));
+        }
+        return;
+      }
+
+      std::vector< std::pair<std::size_t, std::size_t> > ranges;
+      ranges.push_back(std::make_pair(0, nb_points-1));
+
+      std::vector<bool> kept(input.size(), false);
+      if (is_closed) kept[nb_points]=true;
+      while( !ranges.empty() )
+      {
+        std::size_t rb, re;
+        std::tie(rb, re) = ranges.back();
+        ranges.pop_back();
+        kept[rb]=true;
+        kept[re]=true;
+        if (rb+1==re) continue;
+
+        typename Kernel::Line_3 line(get(in_pm, input[rb]), get(in_pm, input[re]));
+        double max_d = max_squared_frechet_distance;
+        std::size_t max_i = 0;
+        for (std::size_t i=rb; i<re; ++i)
+        {
+          double d = squared_distance(line, get(in_pm, input[i]));
+          if (d > max_d)
+          {
+            max_d = d;
+            max_i = i;
+          }
+        }
+        if (max_i != 0)
+        {
+          ranges.push_back( std::make_pair(max_i, re) );
+          ranges.push_back( std::make_pair(rb, max_i) );
+        }
+      }
+      std::size_t nb_kept=0;
+      for (std::size_t i=0; i<input.size(); ++i)
+        if (kept[i]) ++nb_kept;
+
+      output.reserve(nb_kept);
+      for (std::size_t i=0; i<input.size(); ++i)
+        if (kept[i])
+        {
+          output.push_back(input[i]);
+          put(out_pm, output.back(), get(in_pm, input[i]));
+        }
+      //TODO if is_closed-==true, shall we add en extra step to see if we can remove output.front() and output[output.size()-2] (inital endpoints)
+    }
+  }
+}
+
+
+template <typename PointRangeIn, typename PointRangeOut>
+void simplify_polyline(const PointRangeIn& input,
+                             PointRangeOut& output,
+                       const double max_squared_frechet_distance)
+{
+  simplify_polyline(input, output, max_squared_frechet_distance,
+                    parameters::all_default(), parameters::all_default());
+}
+
+template <typename PointRangeIn, typename PointRangeOut, typename NamedParametersIn>
+void simplify_polyline(const PointRangeIn& input,
+                             PointRangeOut& output,
+                       const double max_squared_frechet_distance,
+                       const NamedParametersIn& np_in)
+{
+  simplify_polyline(input, output, max_squared_frechet_distance, np_in, parameters::all_default());
+}
+
+} } } // end of CGAL::Polygon_mesh_processing::experimental namespace
+
+
+#endif // CGAL_POLYGON_MESH_PROCESSING_SIMPLIFY_POLYLINE_H

Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt

@@ -103,6 +103,7 @@ endif()
   create_single_source_cgal_program("test_pmp_manifoldness.cpp")
   create_single_source_cgal_program("test_mesh_smoothing.cpp")
   create_single_source_cgal_program("test_remove_caps_needles.cpp")
+  create_single_source_cgal_program("test_simplify_polylines_pmp.cpp")
 #  create_single_source_cgal_program("test_pmp_repair_self_intersections.cpp")
 
   if( TARGET CGAL::TBB_support )

Polygon_mesh_processing/test/Polygon_mesh_processing/test_simplify_polylines_pmp.cmd

@@ -0,0 +1 @@
+data/U.polylines.txt 0.000001

Polygon_mesh_processing/test/Polygon_mesh_processing/test_simplify_polylines_pmp.cpp

@@ -0,0 +1,60 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Polygon_mesh_processing/internal/simplify_polyline.h>
+
+#include <vector>
+#include <fstream>
+#include <iostream>
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+typedef K::Point_3 Point_3;
+
+namespace PMP_exp=CGAL::Polygon_mesh_processing::experimental;
+
+int main(int argc, char **argv)
+{
+  if (argc!=3)
+  {
+    std::cerr << "Usage: " << argv[0] << " in.polylines.txt max_squared_frechet_distance\n";
+    return 1;
+  }
+
+  std::ifstream in(argv[1]);
+  const double max_squared_frechet_distance = atof(argv[2]);
+
+  //TODO read more than 1 poly or use CGAL I/O function
+  int n;
+  in >> n;
+
+  std::vector<Point_3> polyline(n);
+  for(int i=0; i<n; ++i)
+    in >> polyline[i];
+
+  // first run
+  std::vector<Point_3> simplified_polyline;
+  PMP_exp::simplify_polyline(polyline, simplified_polyline, max_squared_frechet_distance, 
+                             CGAL::parameters::algorithm(PMP_exp::ITERATIVE));
+
+  std::cout << "Polyline of " << n << " points simplified into a polyline of " << simplified_polyline.size() << " points\n";
+
+  std::ofstream out("out_iterative.polylines.txt");
+  out << simplified_polyline.size();
+  for (auto p : simplified_polyline)
+    out << " " << p;
+  out << "\n";
+  out.close();
+
+  // second run
+  simplified_polyline.clear();
+  PMP_exp::simplify_polyline(polyline, simplified_polyline, max_squared_frechet_distance, 
+                             CGAL::parameters::algorithm(PMP_exp::DOUGLAS_PEUCKER));
+
+  std::cout << "Polyline of " << n << " points simplified into a polyline of " << simplified_polyline.size() << " points\n";
+
+  out.open("out_dp.polylines.txt");
+  out << simplified_polyline.size();
+  for (auto p : simplified_polyline)
+    out << " " << p;
+  out << "\n";
+
+  return 0;
+}