Merge pull request #5349 from afabri/Polyline_simplification-simplify_overlaps-GF

2D Polyline simplification: Also simplify overlapping vertex sequences

Installation/CHANGES.md

@@ -32,6 +32,9 @@ A comprehensive list of the supported file formats is available in the Stream_su
 ### 2D Regularized Boolean Set Operations
 -   Added documentation for the free functions `oriented_side(const Point_2& p, ....)` that accept a point and a polygon, and improved the documentation of all other functions
 
+### [2D Polyline Simplification](https://doc.cgal.org/5.3/Manual/packages.html#PkgPolylineSimplification2)
+-   When polylines have common subsequences of vertices these subsequences may be simplifified simultaneously.
+
 [Release 5.2](https://github.com/CGAL/cgal/releases/tag/v5.2)
 -----------
 

Polyline_simplification_2/doc/Polyline_simplification_2/Polyline_simplification_2.txt

@@ -167,6 +167,16 @@ over all vertices of all polyline constraints.
 
 \cgalExample{Polyline_simplification_2/simplify.cpp}
 
+Note that when polylines share subsequences of polyline vertices they can get simplified simultaneously.
+
+\cgalFigureBegin{figure_overlapping_polylines, overlapping_polylines.png}
+Simplification of overlapping subsequence of polyline vertices
+\cgalFigureEnd
+
+\cgalExample{Polyline_simplification_2/simplify_overlapping_polylines.cpp}
+
+
+
 
 \subsection Subsection_PolylineSimplification_Keeping Keeping Points While Removing Vertices
 

Polyline_simplification_2/doc/Polyline_simplification_2/examples.txt

@@ -1,6 +1,7 @@
 /*!
 \example Polyline_simplification_2/simplify_polygon.cpp
 \example Polyline_simplification_2/simplify.cpp
+\example Polyline_simplification_2/simplify_overlapping_polylines.cpp
 \example Polyline_simplification_2/simplify_terrain.cpp
 \example Polyline_simplification_2/points_and_vertices.cpp
 */

Polyline_simplification_2/examples/Polyline_simplification_2/simplify_overlapping_polylines.cpp

@@ -0,0 +1,51 @@
+#include <iostream>
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Polygon_2.h>
+#include <CGAL/Constrained_Delaunay_triangulation_2.h>
+#include <CGAL/Constrained_triangulation_plus_2.h>
+#include <CGAL/Polyline_simplification_2/simplify.h>
+#include <CGAL/Polyline_simplification_2/Squared_distance_cost.h>
+
+namespace PS = CGAL::Polyline_simplification_2;
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+
+typedef PS::Vertex_base_2<K> Vb;
+typedef CGAL::Constrained_triangulation_face_base_2<K> Fb;
+typedef CGAL::Triangulation_data_structure_2<Vb, Fb> TDS;
+typedef CGAL::Constrained_Delaunay_triangulation_2<K, TDS, CGAL::Exact_predicates_tag> CDT;
+typedef CGAL::Constrained_triangulation_plus_2<CDT>     CT;
+
+typedef CT::Point                           Point;
+typedef CT::Constraint_iterator             Constraint_iterator;
+typedef CT::Vertices_in_constraint_iterator Vertices_in_constraint_iterator;
+typedef CT::Points_in_constraint_iterator   Points_in_constraint_iterator;
+typedef PS::Stop_below_count_ratio_threshold Stop;
+typedef PS::Squared_distance_cost Cost;
+
+int main()
+{
+  CT ct;
+
+  std::vector<Point> P = { Point(0,1), Point(1,1), Point(2,2), Point(3,1), Point(4,1),   Point(4,2), Point(4,1), Point(5,1) };
+  std::vector<Point> Q = { Point(5,0), Point(4,1), Point(3,1), Point(2,2), Point(1,1), Point(0,0) };
+  std::vector<Point> R = { Point(3,1), Point(4,1) };
+
+  ct.insert_constraint(P);
+  ct.insert_constraint(Q);
+  ct.insert_constraint(R);
+
+  PS::simplify(ct, Cost(), Stop(0.5));
+
+  for(Constraint_iterator cit = ct.constraints_begin();
+      cit != ct.constraints_end();
+      ++cit) {
+    std::cout << "simplified polyline" << std::endl;
+    for(Points_in_constraint_iterator vit =
+          ct.points_in_constraint_begin(*cit);
+        vit != ct.points_in_constraint_end(*cit);
+        ++vit)
+      std::cout << *vit << std::endl;
+  }
+  return 0;
+}

Polyline_simplification_2/include/CGAL/Polyline_simplification_2/Vertex_base_2.h

@@ -37,6 +37,9 @@ class Vertex_base_2
   bool m_removable;
   FT m_cost;
 
+public:
+  std::size_t ID;
+
 #ifndef DOXYGEN_RUNNING
 public:
   template < typename TDS2 >

Polyline_simplification_2/include/CGAL/Polyline_simplification_2/simplify.h

@@ -17,6 +17,7 @@
 #include <CGAL/disable_warnings.h>
 
 #include <list>
+#include <unordered_map>
 
 #include <CGAL/Polyline_simplification_2/Vertex_base_2.h>
 #include <CGAL/Polyline_simplification_2/Squared_distance_cost.h>
@@ -57,11 +58,15 @@ class Polyline_simplification_2
 public:
 
   typedef typename PCT::Point Point;
+  typedef typename PCT::Edge Edge;
   typedef typename PCT::Constraint_id Constraint_id;
+  typedef typename PCT::Constrained_edges_iterator Constrained_edges_iterator;
   typedef typename PCT::Constraint_iterator Constraint_iterator;
   typedef typename PCT::Vertices_in_constraint_iterator Vertices_in_constraint_iterator;
+  typedef typename PCT::Finite_vertices_iterator Finite_vertices_iterator;
   //typedef typename PCT::Points_in_constraint_iterator Points_in_constraint_iterator;
   typedef typename PCT::Vertex_handle Vertex_handle;
+    typedef typename PCT::Face_handle Face_handle;
   typedef typename PCT::Vertex_circulator Vertex_circulator;
 
   typedef typename PCT::Geom_traits::FT FT;
@@ -71,6 +76,7 @@ public:
   StopFunction stop;
   std::size_t pct_initial_number_of_vertices, number_of_unremovable_vertices;
 
+  std::unordered_map<Vertex_handle, std::list<Vertices_in_constraint_iterator> > vertex_to_iterator;
 
   struct Compare_cost
   {
@@ -79,6 +85,12 @@ public:
     {
       return (*x)->cost() < (*y)->cost();
     }
+
+    bool operator() (const Vertex_handle& x,const Vertex_handle& y) const
+    {
+      return x->cost() < y->cost();
+    }
+
   } ;
 
   struct Id_map : public boost::put_get_helper<std::size_t, Id_map>
@@ -86,12 +98,16 @@ public:
     typedef boost::readable_property_map_tag category;
     typedef std::size_t                      value_type;
     typedef value_type                       reference;
-    typedef Vertices_in_constraint_iterator  key_type;
+    typedef Vertex_handle  key_type;
 
-    reference operator[] ( key_type const& x ) const { return x.base()->id ; }
+
+    reference operator[] ( key_type const& x ) const
+    {
+        return x->ID;
+    }
   } ;
 
-  typedef CGAL::Modifiable_priority_queue<Vertices_in_constraint_iterator,Compare_cost,Id_map> MPQ ;
+  typedef CGAL::Modifiable_priority_queue<Vertex_handle,Compare_cost,Id_map> MPQ ;
 
   MPQ* mpq;
 
@@ -124,24 +140,64 @@ public:
     delete mpq;
   }
 
+  // endpoints of constraints are unremovable
+  // vertices which are not endpoint and have != 2 incident constrained edges are unremovable
   void initialize_unremovable()
   {
-    std::set<Vertex_handle> vertices;
     Constraint_iterator cit = pct.constraints_begin(), e = pct.constraints_end();
     for(; cit!=e; ++cit){
       Constraint_id cid = *cit;
-      Vertices_in_constraint_iterator it = pct.vertices_in_constraint_begin(cid);
+      Vertices_in_constraint_iterator it = pct.vertices_in_constraint_begin(cid),
+                                      ite = pct.vertices_in_constraint_end(cid);
       (*it)->set_removable(false);
-      for(; it != pct.vertices_in_constraint_end(cid); ++it){
-        if(vertices.find(*it) != vertices.end()){
+      ++it;
+      for(; it != ite; ++it){
+        if((boost::next(it) != ite) && (boost::prior(it)== boost::next(it))){
           (*it)->set_removable(false);
-        } else {
-          vertices.insert(*it);
         }
       }
       it = boost::prior(it);
       (*it)->set_removable(false);
     }
+
+    std::unordered_map<Vertex_handle, int> degrees;
+    for (Constrained_edges_iterator it = pct.constrained_edges_begin(); it != pct.constrained_edges_end(); ++it) {
+      Edge e = *it;
+      Face_handle fh = e.first;
+      int ei = e.second;
+      Vertex_handle vh = fh->vertex(pct.cw(ei));
+      ++degrees[vh];
+      vh = fh->vertex(pct.ccw(ei));
+      ++degrees[vh];
+    }
+
+    for(Finite_vertices_iterator it = pct.finite_vertices_begin(); it != pct.finite_vertices_end(); ++it){
+      if( it->is_removable() && (degrees[it] != 2) ){
+        it->set_removable(false);
+      }
+    }
+
+    cit = pct.constraints_begin(), e = pct.constraints_end();
+    for(; cit!=e; ++cit){
+      Constraint_id cid = *cit;
+      for(Vertices_in_constraint_iterator it = pct.vertices_in_constraint_begin(cid);
+          it != pct.vertices_in_constraint_end(cid);
+          ++it){
+        if((*it)->is_removable()){
+          typename std::unordered_map<Vertex_handle, std::list<Vertices_in_constraint_iterator> >::iterator lit;
+          lit = vertex_to_iterator.find(*it);
+
+          if(lit != vertex_to_iterator.end()){
+            std::list<Vertices_in_constraint_iterator>& ilist = lit->second;
+            if(std::find(ilist.begin(),ilist.end(),it) == ilist.end()){
+              ilist.push_back(it);
+            }
+          }else{
+            vertex_to_iterator[*it].push_back(it);
+          }
+        }
+      }
+    }
   }
 
   // For all polyline constraints we compute the cost of all unremovable and not removed vertices
@@ -156,7 +212,9 @@ public:
         boost::optional<FT> dist = cost(pct, it);
         if(dist){
           (*it)->set_cost(*dist);
-          (*mpq).push(it);
+          if(! (*mpq).contains(*it)){
+              (*mpq).push(*it);
+            }
           ++n;
         } else {
           // no need to set the costs as this vertex is not in the priority queue
@@ -230,7 +288,9 @@ public:
     for(Vertices_in_constraint_iterator it = pct.vertices_in_constraint_begin(cid);
         it != pct.vertices_in_constraint_end(cid);
         ++it){
-      it.base()->id = id++;
+        Vertex_handle vh = *it;
+        vh->ID = id++;
+        //vertex_index_map[vh] = id++;
     }
     return id;
   }
@@ -239,9 +299,9 @@ public:
   initialize_indices()
   {
     int id = 0;
-    Constraint_iterator b = pct.constraints_begin(), e = pct.constraints_end();
-    for(; b!=e; ++b){
-      id = initialize_indices(*b, id);
+
+    for(Finite_vertices_iterator it =  pct.finite_vertices_begin(); it != pct.finite_vertices_end(); ++it){
+      it->ID = id++;
     }
     return id;
   }
@@ -252,29 +312,31 @@ operator()()
   if((*mpq).empty()){
       return false;
     }
-  Vertices_in_constraint_iterator v = (*mpq).top();
+  Vertex_handle v = (*mpq).top();
   (*mpq).pop();
-  if(stop(pct, *v, (*v)->cost(), pct_initial_number_of_vertices, pct.number_of_vertices())){
+  if(stop(pct, v, v->cost(), pct_initial_number_of_vertices, pct.number_of_vertices())){
     return false;
   }
-  if(is_removable(v)){
-    Vertices_in_constraint_iterator u = boost::prior(v), w = boost::next(v);
-    pct.simplify(v);
+
+  Vertices_in_constraint_iterator vit = vertex_to_iterator[v].front();
+  if(is_removable(vit)){
+    Vertices_in_constraint_iterator u = boost::prior(vit), w = boost::next(vit);
+    pct.simplify(vit);
 
     if((*u)->is_removable()){
       boost::optional<FT> dist = cost(pct, u);
       if(! dist){
         // cost is undefined
-        if( mpq->contains(u) ){
-          mpq->erase(u);
+        if( mpq->contains(*u) ){
+          mpq->erase(*u);
         }
       } else {
         (*u)->set_cost(*dist);
-        if(mpq->contains(u)){
-          mpq->update(u, true);
+        if(mpq->contains(*u)){
+          mpq->update(*u, true);
         }
         else{
-          mpq->push(u);
+          mpq->push(*u);
         }
       }
     }
@@ -283,16 +345,16 @@ operator()()
       boost::optional<FT> dist = cost(pct, w);
       if(! dist){
         // cost is undefined
-        if( mpq->contains(w) ){
-          mpq->erase(w);
+        if( mpq->contains(*w) ){
+          mpq->erase(*w);
         }
       } else {
         (*w)->set_cost(*dist);
-        if(mpq->contains(w)){
-          mpq->update(w, true);
+        if(mpq->contains(*w)){
+          mpq->update(*w, true);
         }
         else{
-          mpq->push(w);
+          mpq->push(*w);
         }
 
       }

Triangulation_2/include/CGAL/Triangulation_2/internal/Polyline_constraint_hierarchy_2.h

@@ -46,14 +46,13 @@ private:
   class Node {
   public:
     explicit Node(Vertex_handle vh, bool input = false)
-      : vertex_(vh), id(-1), input(input)
+      : vertex_(vh), input(input)
     {}
     const Point& point() const { return vertex_->point(); }
     Vertex_handle vertex() const { return vertex_; }
   private:
     Vertex_handle vertex_;
   public:
-    int id;
     bool input;
   };
 
@@ -585,45 +584,58 @@ void Polyline_constraint_hierarchy_2<T,Compare,Point>::simplify(Vertex_it uc,
                                                        Vertex_it wc)
 
 {
+  // TODO: How do we (want to) deal with u == w ???
   Vertex_handle u = *uc, v = *vc, w = *wc;
   typename Sc_to_c_map::iterator uv_sc_iter = sc_to_c_map.find(make_edge(u, v));
-  CGAL_assertion_msg( uv_sc_iter != sc_to_c_map.end(), "not a subconstraint" );
+  typename Sc_to_c_map::iterator vw_sc_iter = sc_to_c_map.find(make_edge(v, w));
   Context_list*  uv_hcl = uv_sc_iter->second;
-  CGAL_assertion_msg((u == w) || (uv_hcl->size() == 1), "more than one constraint passing through the subconstraint" );
-
-  if(*(uv_hcl->front().current()) != u) {
-    std::swap(u,w);
-    uv_sc_iter = sc_to_c_map.find(make_edge(u, v));
-    CGAL_assertion_msg( uv_sc_iter != sc_to_c_map.end(), "not a subconstraint" );
-    uv_hcl = (*uv_sc_iter).second;
-    CGAL_assertion_msg((u == w) || (uv_hcl->size() == 1), "more than one constraint passing through the subconstraint" );
-  }
-  // now u,v, and w are ordered along the polyline constraint
+  Context_list*  vw_hcl = vw_sc_iter->second;
+  // AF:  what is input() about???
   if(vc.input()){
     uc.input() = true;
     wc.input() = true;
   }
-  typename Sc_to_c_map::iterator vw_sc_iter = sc_to_c_map.find(make_edge(v, w));
-  CGAL_assertion_msg( vw_sc_iter != sc_to_c_map.end(), "not a subconstraint" );
-  Context_list*  vw_hcl = vw_sc_iter->second;
-    CGAL_assertion_msg((u == w) || (vw_hcl->size() == 1), "more than one constraint passing through the subconstraint" );
 
-  Vertex_list* vertex_list = uv_hcl->front().id().vl_ptr();
-  CGAL_assertion_msg(vertex_list  == vw_hcl->front().id().vl_ptr(), "subconstraints from different polyline constraints" );
+  // Take contexts from the two context lists depending on the orientation of the constraints
+  // These are the contexts where current is either u or w
+  // remove from uv_hcl the contexts where current is not u
+  // remove from vw_hcl the contexts where current is not w
+  // splice into uv_hcl
+  typename Context_list::iterator it = uv_hcl->begin();
+  while(it != uv_hcl->end()){
+    if((*it->current()) != u){
+      it = uv_hcl->erase(it);
+    }else{
       // Remove the list item which points to v
-  vertex_list->skip(vc.base());
+      Vertex_list* vertex_list = it->id().vl_ptr();
+      Vertex_it vc_in_context = it->current();
+      vc_in_context = boost::next(vc_in_context);
+      vertex_list->skip(vc_in_context.base());
+      ++it;
+    }
+  }
+  it = vw_hcl->begin();
+  while(it != vw_hcl->end()){
+    if((*it->current()) != w){
+      it = vw_hcl->erase(it);
+    }else{
+      // Remove the list item which points to v
+      Vertex_list* vertex_list = it->id().vl_ptr();
+      Vertex_it vc_in_context = it->current();
+      vc_in_context = boost::next(vc_in_context);
+      vertex_list->skip(vc_in_context.base());
+      ++it;
+    }
+  }
+
+  uv_hcl->splice(uv_hcl->end(),*vw_hcl);
+  delete vw_hcl;
 
-  if(u != w){
-    // Remove the entries for [u,v] and [v,w]
   sc_to_c_map.erase(uv_sc_iter);
   sc_to_c_map.erase(vw_sc_iter);
-    delete vw_hcl;
+
   // reuse other context list
   sc_to_c_map[make_edge(u,w)] = uv_hcl;
-  }else{
-    sc_to_c_map.erase(uv_sc_iter);
-    delete vw_hcl;
-  }
 }