Merge pull request #3910 from maxGimeno/Polgon_mesh_slicer-Fix_orientation-GF

PMP::Orient polylines of Polgon Mesh Slicer ccw.

Polygon_mesh_processing/include/CGAL/Polygon_mesh_slicer.h

@@ -40,6 +40,7 @@
 #include <boost/variant.hpp>
 
 #include <CGAL/boost/graph/split_graph_into_polylines.h>
+#include <CGAL/boost/graph/helpers.h>
 
 namespace CGAL {
 
@@ -151,6 +152,7 @@ class Polygon_mesh_slicer
     VertexPointMap m_vpmap;
     typename Traits_::Intersect_3 intersect_3;
     OutputIterator out;
+    std::pair<AL_vertex_descriptor, AL_vertex_descriptor> nodes_for_orient;
 
     Polyline_visitor( TriangleMesh& tmesh,
                       AL_graph& al_graph,
@@ -167,12 +169,102 @@ class Polygon_mesh_slicer
     {}
 
     std::vector< Point_3 > current_poly;
+
+    // returns true iff the polyline is not correctly oriented
+    // Using the first edge is oriented such that the normal induced by the face
+    // containing the edge, the oriented edge and the plane orthogonal vector defines
+    // a direct orthogonal basis
+    bool do_reverse_polyline()
+    {
+      if (current_poly.size() < 2)
+        return false;
+
+      AL_vertex_info v1 = al_graph[nodes_for_orient.first];
+      AL_vertex_info v2 = al_graph[nodes_for_orient.second];
+
+      if (const vertex_descriptor* vd1_ptr = boost::get<vertex_descriptor>(&v1) )
+      {
+        if (const vertex_descriptor* vd2_ptr = boost::get<vertex_descriptor>(&v2) )
+        {
+          CGAL_assertion( halfedge(*vd1_ptr, *vd2_ptr, m_tmesh).second );
+          halfedge_descriptor h_opp = halfedge(*vd1_ptr, *vd2_ptr, m_tmesh).first;
+          if ( !is_border(h_opp, m_tmesh) )
+          {
+            CGAL_assertion(source(h_opp, m_tmesh) == *vd1_ptr);
+            CGAL::Oriented_side os = m_plane.oriented_side( get(m_vpmap, target(next(h_opp, m_tmesh), m_tmesh) ) );
+            if (os != CGAL::ON_ORIENTED_BOUNDARY)
+              return m_plane.oriented_side( get(m_vpmap, target(next(h_opp, m_tmesh), m_tmesh)) ) == CGAL::ON_NEGATIVE_SIDE;
+          }
+          h_opp = opposite(h_opp, m_tmesh);
+          if ( !is_border(h_opp, m_tmesh) )
+          {
+            CGAL_assertion(source(h_opp, m_tmesh) == *vd2_ptr);
+            CGAL::Oriented_side os = m_plane.oriented_side( get(m_vpmap, target(next(h_opp, m_tmesh), m_tmesh) ) );
+            if (os != CGAL::ON_ORIENTED_BOUNDARY)
+              return m_plane.oriented_side( get(m_vpmap, target(next(h_opp, m_tmesh), m_tmesh)) ) == CGAL::ON_POSITIVE_SIDE;
+          }
+          return false; // since coplanar edges are filtered out, we should never end up here.
+        }
+        else
+        {
+          // e2 is intersected in its interior
+          edge_descriptor e2 = boost::get<edge_descriptor>(v2);
+          halfedge_descriptor h2 = halfedge(e2, m_tmesh);
+          if ( target(next(h2, m_tmesh), m_tmesh) != *vd1_ptr )
+            h2=opposite(h2, m_tmesh);
+          return m_plane.oriented_side( get(m_vpmap, source(h2, m_tmesh)) ) == CGAL::ON_POSITIVE_SIDE;
+        }
+      }
+      else
+      {
+        edge_descriptor e1 = boost::get<edge_descriptor>(v1);
+        halfedge_descriptor h1 = halfedge(e1, m_tmesh);
+        if (const vertex_descriptor* vd2_ptr = boost::get<vertex_descriptor>(&v2) )
+        {
+          // e1 is intersected in its interior
+          if ( target(next(h1, m_tmesh), m_tmesh) != *vd2_ptr )
+            h1=opposite(h1, m_tmesh);
+          CGAL_assertion( target(next(h1, m_tmesh), m_tmesh) == *vd2_ptr );
+        }
+        else
+        {
+          // intersection in the interior of both edges
+          edge_descriptor e2 = boost::get<edge_descriptor>(v2);
+          halfedge_descriptor h2 = halfedge(e2, m_tmesh);
+          if ( face(h1, m_tmesh) != face(h2,m_tmesh) )
+          {
+            halfedge_descriptor opp_h1 = opposite(h1, m_tmesh),
+                                opp_h2 = opposite(h2, m_tmesh);
+            if ( face(opp_h1, m_tmesh) == face(opp_h2, m_tmesh) )
+            {
+              h1=opp_h1;
+              h2=opp_h2;
+            }
+            else
+              if ( face(opp_h1, m_tmesh)==face(h2,m_tmesh) )
+                h1=opp_h1;
+              else
+                h2=opp_h2;
+          }
+          CGAL_assertion( face(h1, m_tmesh) == face(h2,m_tmesh) );
+        }
+        return m_plane.oriented_side( get(m_vpmap, source(h1, m_tmesh)) ) == CGAL::ON_NEGATIVE_SIDE;
+      }
+    }
+
     void start_new_polyline()
     {
       current_poly.clear();
     }
+
     void add_node(AL_vertex_descriptor node_id)
     {
+      if (current_poly.empty())
+        nodes_for_orient.first=node_id;
+      else
+        if (current_poly.size()==1)
+          nodes_for_orient.second=node_id;
+
       AL_vertex_info v = al_graph[node_id];
       if (const vertex_descriptor* vd_ptr = boost::get<vertex_descriptor>(&v) )
       {
@@ -196,6 +288,8 @@ class Polygon_mesh_slicer
     void end_polyline()
     {
       CGAL_assertion(!current_poly.empty());
+      if(do_reverse_polyline())
+        std::reverse(current_poly.begin(), current_poly.end());
       *out++=current_poly;
     }
   };
@@ -384,6 +478,18 @@ public:
 
   /**
    * Constructs the intersecting polylines of `plane` with the input triangulated surface mesh.
+   *
+   * If a polyline is closed, the first and last point of that polyline will be identical.
+   *
+   * Each resulting polyline `P` is oriented such that for two consecutive points `p` and `q` in `P`,
+   * the normal vector of the face(s) containing the segment `pq`, the vector `pq`, and the orthogonal vertor of `plane`
+   * is a direct orthogonal basis.
+   * The normal vector of each face is chosen  to point on the side of the face
+   * where its sequence of vertices is seen counterclockwise.
+   *
+   * Note that an edge shared by two faces included in `plane` will not be reported. For example,
+   * if `plane` passes though one face of a cube, only one closed polyline will be reported (the boundary of the face)
+   *
    * @tparam OutputIterator an output iterator accepting polylines.
    *              A polyline is provided as `std::vector<Traits::Point_3>`.
    *              A polyline is closed if its first and last point are identical.

Polygon_mesh_processing/test/Polygon_mesh_processing/polygon_mesh_slicer_test.cpp

@@ -10,29 +10,47 @@
 #include <CGAL/AABB_halfedge_graph_segment_primitive.h>
 
 #include <CGAL/Polygon_mesh_slicer.h>
+#include <CGAL/Polygon_mesh_processing/orientation.h>
 #include <CGAL/AABB_tree.h>
 #include <CGAL/AABB_traits.h>
-
+#include <CGAL/Polygon_2.h>
+#include <CGAL/use.h>
 
 #include <fstream>
 
 typedef CGAL::Exact_predicates_inexact_constructions_kernel Epic;
 typedef CGAL::Exact_predicates_inexact_constructions_kernel Epec;
 
+template<class K, class Polyline_type>
+bool is_ccw(int xi, int yi,
+    const Polyline_type& polyline)
+{
+  CGAL::Polygon_2<K> polygon;
+  if(polyline.front() == polyline.back())
+  {
+    BOOST_FOREACH(const typename K::Point_3& p, polyline)
+    {
+      polygon.push_back(typename K::Point_2(p[xi], p[yi]));
+    }
+    polygon.erase(polygon.vertices_end()-1);
+    return polygon.is_counterclockwise_oriented();
+  }
+  return true;
+}
 template <typename K>
 int test_slicer()
 {
 #ifdef USE_SURFACE_MESH
-typedef CGAL::Surface_mesh<K::Point_3> Mesh;
+  typedef CGAL::Surface_mesh<K::Point_3> Mesh;
 #else
-typedef CGAL::Polyhedron_3<K> Mesh;
+  typedef CGAL::Polyhedron_3<K> Mesh;
 #endif
 
-typedef CGAL::AABB_halfedge_graph_segment_primitive<Mesh> HGSP;
+  typedef CGAL::AABB_halfedge_graph_segment_primitive<Mesh> HGSP;
-typedef CGAL::AABB_traits<K, HGSP>    AABB_traits;
+  typedef CGAL::AABB_traits<K, HGSP>    AABB_traits;
-typedef CGAL::AABB_tree<AABB_traits>  AABB_tree;
+  typedef CGAL::AABB_tree<AABB_traits>  AABB_tree;
-typedef std::vector<typename K::Point_3> Polyline_type;
+  typedef std::vector<typename K::Point_3> Polyline_type;
-typedef std::list< Polyline_type > Polylines;
+  typedef std::vector< Polyline_type > Polylines;
 
   //API test
   {
@@ -95,6 +113,56 @@ typedef std::list< Polyline_type > Polylines;
   slicer(typename K::Plane_3(0,0,1,333), std::back_inserter(polylines));
   assert(polylines.empty());
 
+// Now test the orientation of polylines
+  polylines.clear();
+  slicer(typename K::Plane_3(0,1,0,0.5), std::back_inserter(polylines));
+  assert(polylines.size()==2); // two polylines
+  int closed_id = polylines.front().front()==polylines.front().back() ? 0 : 1;
+  CGAL_USE(closed_id);
+  CGAL_assertion( is_ccw<K>(0, 2 , polylines[closed_id]) );
+
+  polylines.clear();
+  slicer(typename K::Plane_3(0,-1,0,-0.5), std::back_inserter(polylines));
+  assert(polylines.size()==2); // two polylines
+  closed_id = polylines.front().front()==polylines.front().back() ? 0 : 1;
+  CGAL_assertion( !is_ccw<K>(0, 2, polylines[closed_id]) );
+
+  polylines.clear();
+  slicer(typename K::Plane_3(0,0,1,1), std::back_inserter(polylines));
+  assert(polylines.size()==1); // one polyline
+  CGAL_assertion( is_ccw<K>(0, 1 , polylines[0]) );
+
+  polylines.clear();
+  slicer(typename K::Plane_3(0,0,-1,-1), std::back_inserter(polylines));
+  assert(polylines.size()==1); // one polyline
+  CGAL_assertion( !is_ccw<K>(0, 1 , polylines[0]) );
+
+  // reverse face orientation (no need to rebuild the tree)
+  CGAL::Polygon_mesh_processing::reverse_face_orientations(m);
+  polylines.clear();
+  slicer(typename K::Plane_3(0,1,0,0.5), std::back_inserter(polylines));
+  assert(polylines.size()==2); // two polylines
+  closed_id = polylines.front().front()==polylines.front().back() ? 0 : 1;
+  CGAL_assertion( !is_ccw<K>(0, 2 , polylines[closed_id]) );
+
+  polylines.clear();
+  slicer(typename K::Plane_3(0,-1,0,-0.5), std::back_inserter(polylines));
+  assert(polylines.size()==2); // two polylines
+  closed_id = polylines.front().front()==polylines.front().back() ? 0 : 1;
+  CGAL_assertion( is_ccw<K>(0, 2, polylines[closed_id]) );
+
+  polylines.clear();
+  slicer(typename K::Plane_3(0,0,1,1), std::back_inserter(polylines));
+  assert(polylines.size()==1); // one polyline
+  CGAL_assertion( !is_ccw<K>(0, 1 , polylines[0]) );
+
+  polylines.clear();
+  slicer(typename K::Plane_3(0,0,-1,-1), std::back_inserter(polylines));
+  assert(polylines.size()==1); // one polyline
+  CGAL_assertion( is_ccw<K>(0, 1 , polylines[0]) );
+
+
+
   return 0;
 }