use a specific traits when the plane is parallel to two axis

Polygon_mesh_processing/include/CGAL/Polygon_mesh_slicer_3.h

@@ -32,6 +32,7 @@
 #include <boost/foreach.hpp>
 #include <boost/graph/adjacency_list.hpp>
 #include <CGAL/internal/Polygon_mesh_slicer_3/Traversal_traits.h>
+#include <CGAL/internal/Polygon_mesh_slicer_3/Axis_parallel_plane_traits.h>
 
 #include <boost/variant.hpp>
 
@@ -50,7 +51,7 @@ namespace CGAL {
 /// \tparam AABBTree must be an instanciation of `CGAL::AABB_tree` able to handle
 ///         the edges of TriangleMesh, having its `edge_descriptor` as primitive id.
 /// Depends on \ref PkgAABB_treeSummary
-/// \todo use dedicated predicates if the plane is axis aligned
+/// \todo document usage of custom traits when plane is Plane(1,0,0,d), Plane(0,1,0,d), Plane(0,0,1,d)
 template<class TriangleMesh,
   class Traits,
   class VertexPointPmap = typename boost::property_map< TriangleMesh, vertex_point_t>::type,
@@ -69,8 +70,8 @@ class Polygon_mesh_slicer_3
 /// Geometric typedefs
   typedef typename Traits::Plane_3                                      Plane_3;
   typedef typename Traits::Segment_3                                  Segment_3;
-  typedef typename Traits::Intersect_3                              Intersect_3;
   typedef typename Traits::Point_3                                      Point_3;
+  typedef typename Traits::FT                                                FT;
 
 /// typedefs for internal graph to get connectivity of the polylines
   typedef boost::variant<vertex_descriptor, edge_descriptor>     AL_vertex_info;
@@ -83,13 +84,22 @@ class Polygon_mesh_slicer_3
   typedef std::pair<AL_vertex_descriptor, AL_vertex_descriptor>  AL_vertex_pair;
   typedef std::map<vertex_descriptor, AL_vertex_descriptor>        Vertices_map;
   typedef std::pair<const vertex_descriptor,AL_vertex_descriptor>   Vertex_pair;
-/// Traversal traits for the AABB-tree selecting edges and classifying them
-  typedef Polygon_mesh_slicer::Traversal_traits<  AL_graph,
-                                                  TriangleMesh,
-                                                  VertexPointPmap,
-                                                  typename AABBTree::AABB_traits,
-                                                  Traits >     Traversal_traits;
+/// Traversal traits
+  typedef Polygon_mesh_slicer::Traversal_traits<
+    AL_graph,
+    TriangleMesh,
+    VertexPointPmap,
+    typename AABBTree::AABB_traits,
+    Traits >                                           General_traversal_traits;
 
+  typedef Polygon_mesh_slicer::Traversal_traits<
+    AL_graph,
+    TriangleMesh,
+    VertexPointPmap,
+    typename AABBTree::AABB_traits,
+    Polygon_mesh_slicer::Axis_parallel_plane_traits<Traits>
+  >                                              Axis_parallel_traversal_traits;
+/// Auxiliary classes
   // compare the faces using the halfedge descriptors
   struct Compare_face{
     TriangleMesh& m_tmesh;
@@ -105,26 +115,26 @@ class Polygon_mesh_slicer_3
 
   typedef std::map< halfedge_descriptor, AL_vertex_pair, Compare_face > AL_edge_map;
 
-  template <class OutputIterator>
+  template <class OutputIterator, class Traits_>
   struct Polyline_visitor{
     AL_graph& al_graph;
     TriangleMesh& m_tmesh;
     const Plane_3& m_plane;
     VertexPointPmap m_vpmap;
-    const Traits& m_traits;
+    typename Traits_::Intersect_3 intersect_3;
     OutputIterator out;
 
     Polyline_visitor( TriangleMesh& tmesh,
                       AL_graph& al_graph,
                       const Plane_3& plane,
                       VertexPointPmap vpmap,
-                      const Traits& traits,
+                      const Traits_& traits,
                       OutputIterator out)
       : al_graph(al_graph)
       , m_tmesh(tmesh)
       , m_plane(plane)
       , m_vpmap(vpmap)
-      , m_traits(traits)
+      , intersect_3( traits.intersect_3_object() )
       , out(out)
     {}
 
@@ -147,9 +157,8 @@ class Polygon_mesh_slicer_3
           get(m_vpmap, source(ed, m_tmesh)),
           get(m_vpmap,target(ed, m_tmesh))
         );
-        Intersect_3 intersection = m_traits.intersect_3_object();
-        typename cpp11::result_of<Intersect_3(Plane_3, Segment_3)>::type
-          inter = intersection(m_plane, s);
+        typename cpp11::result_of<typename Traits_::Intersect_3(Plane_3, Segment_3)>::type
+          inter = intersect_3(m_plane, s);
         CGAL_assertion( inter );
         const Point_3* pt_ptr = boost::get<Point_3>(&(*inter));
         current_poly.push_back( *pt_ptr );
@@ -185,7 +194,7 @@ class Polygon_mesh_slicer_3
   {
     return face( opposite( halfedge(ed, m_tmesh), m_tmesh), m_tmesh);
   }
-
+/// Other private functions
   /// handle edge insertion in the adjacency_list graph
   /// we add an edge betweem two edge_descriptor if they
   /// share a common facet
@@ -232,6 +241,33 @@ class Polygon_mesh_slicer_3
       }
   }
 
+  std::pair<int, FT>
+  axis_parallel_plane_info(const Plane_3& plane) const
+  {
+    FT a = m_traits.compute_a_3_object()(plane);
+    FT b = m_traits.compute_b_3_object()(plane);
+    FT c = m_traits.compute_c_3_object()(plane);
+    FT d = m_traits.compute_d_3_object()(plane);
+
+    if (a==0)
+    {
+      if (b==0)
+      {
+        if (c==1 || c==-1)
+          return std::pair<int,FT>(2, -d*c); /// z=-d
+      }
+      else
+      {
+        if (c==0 && (b==1 || b==-1))
+          return std::pair<int,FT>(1, -d*b); /// y=-d
+      }
+    }
+    else
+      if (b==0 && c==0 && ( a==1 || a==-1))
+        return std::pair<int,FT>(0, -d*a); /// x=-d
+    return std::pair<int,FT>(-1, 0);
+  }
+
 public:
 
   /**
@@ -337,15 +373,35 @@ public:
     Vertices_map vertices;
 
     /// get all edges intersected by the plane and classify them
-    Traversal_traits ttraits(
-      all_coplanar_edges,
-      iedges,
-      vertices,
-      m_tmesh,
-      m_vpmap,
-      m_tree_ptr->traits(),
-      m_traits);
-    m_tree_ptr->traversal(plane, ttraits);
+
+    std::pair<int, FT> app_info = axis_parallel_plane_info(plane);
+    if (app_info.first==-1)
+    {
+      General_traversal_traits ttraits(
+        all_coplanar_edges,
+        iedges,
+        vertices,
+        m_tmesh,
+        m_vpmap,
+        m_tree_ptr->traits(),
+        m_traits);
+      m_tree_ptr->traversal(plane, ttraits);
+    }
+    else
+    {
+      Polygon_mesh_slicer::Axis_parallel_plane_traits<Traits>
+        traits(app_info.first, app_info.second, m_traits);
+
+      Axis_parallel_traversal_traits ttraits(
+        all_coplanar_edges,
+        iedges,
+        vertices,
+        m_tmesh,
+        m_vpmap,
+        m_tree_ptr->traits(),
+        traits);
+      m_tree_ptr->traversal(plane, ttraits);
+    }
 
     /// init output graph
     AL_graph al_graph;
@@ -435,9 +491,22 @@ public:
 
     /// now assemble the edges of al_graph to define polylines,
     /// putting them in the output iterator
-    Polyline_visitor<OutputIterator> visitor(m_tmesh, al_graph, plane, m_vpmap, m_traits, out);
-    split_graph_into_polylines(al_graph, visitor);
-    return visitor.out;
+    if (app_info.first==-1)
+    {
+      Polyline_visitor<OutputIterator, Traits> visitor(m_tmesh, al_graph, plane, m_vpmap, m_traits, out);
+      split_graph_into_polylines(al_graph, visitor);
+      return visitor.out;
+    }
+    else
+    {
+      typedef Polygon_mesh_slicer::Axis_parallel_plane_traits<Traits> App_traits;
+      App_traits app_traits(app_info.first, app_info.second, m_traits);
+
+      Polyline_visitor<OutputIterator, App_traits> visitor
+        (m_tmesh, al_graph, plane, m_vpmap, app_traits, out);
+      split_graph_into_polylines(al_graph, visitor);
+      return visitor.out;
+    }
   }
 
   ~Polygon_mesh_slicer_3()

Polygon_mesh_processing/include/CGAL/internal/Polygon_mesh_slicer_3/Axis_parallel_plane_traits.h

@@ -0,0 +1,160 @@
+// Copyright (c) 2015 GeometryFactory (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)     :  Sebastien Loriot
+
+#include <CGAL/enum.h>
+#include <CGAL/Bbox_3.h>
+#include <CGAL/array.h>
+
+#include <boost/optional.hpp>
+#include <boost/variant.hpp>
+
+
+#ifndef CGAL_INTERNAL_POLYGON_MESH_SLICER_3_AXIS_PARALLEL_PLANE_TRAITS_H
+#define CGAL_INTERNAL_POLYGON_MESH_SLICER_3_AXIS_PARALLEL_PLANE_TRAITS_H
+
+namespace CGAL{
+namespace Polygon_mesh_slicer{
+
+template <class Traits>
+class Axis_parallel_plane_traits
+{
+  typedef typename Traits::FT FT;
+
+  const Traits& m_traits;
+  const int m_cst_coord; // 0, 1 or 2 indicates which coordinates is constant
+  const FT m_value; // indicates the value of the constant coordinate
+
+public:
+  typedef typename Traits::Plane_3 Plane_3;
+
+  Axis_parallel_plane_traits(int cst_coord, FT value, const Traits& traits)
+    : m_traits(traits)
+    , m_cst_coord(cst_coord)
+    , m_value(value)
+  {}
+
+  struct Oriented_side_3
+  {
+    const int m_cst_coord;
+    const FT m_value;
+    const typename Traits::Construct_cartesian_const_iterator_3 m_coord_iterator;
+
+    typedef Oriented_side result_type;
+
+    Oriented_side_3(const Axis_parallel_plane_traits<Traits>& traits)
+      : m_cst_coord(traits.m_cst_coord)
+      , m_value(traits.m_value)
+      , m_coord_iterator(traits.m_traits.construct_cartesian_const_iterator_3_object())
+    {}
+
+    result_type
+    operator()(const typename Traits::Plane_3&, const typename Traits::Point_3& pt) const
+    {
+      if ( *( m_coord_iterator(pt)+m_cst_coord) > m_value ) return ON_POSITIVE_SIDE;
+      if ( *( m_coord_iterator(pt)+m_cst_coord) < m_value ) return ON_NEGATIVE_SIDE;
+      return ON_ORIENTED_BOUNDARY;
+    }
+  };
+
+  struct Do_intersect_3{
+    const int m_cst_coord;
+    const FT m_value;
+
+    typedef bool result_type;
+
+    Do_intersect_3(int cst_coord, FT value)
+      : m_cst_coord(cst_coord)
+      , m_value(value)
+    {}
+
+    result_type
+    operator()(const typename Traits::Plane_3&, const Bbox_3& bbox) const
+    {
+      return bbox.min(m_cst_coord) <= m_value &&
+             bbox.max(m_cst_coord) >= m_value;
+    }
+  };
+
+  struct Intersect_3{
+    const int m_cst_coord;
+    const FT m_value;
+    const typename Traits::Construct_cartesian_const_iterator_3 m_coord_iterator;
+    const typename Traits::Construct_point_3 m_point_3;
+    const typename Traits::Construct_source_3 m_source_3;
+    const typename Traits::Construct_target_3 m_target_3;
+
+    typedef boost::variant<typename Traits::Point_3, typename Traits::Segment_3> Variant_type;
+    typedef boost::optional< Variant_type > result_type;
+
+    Intersect_3(const Axis_parallel_plane_traits<Traits>& traits)
+      : m_cst_coord(traits.m_cst_coord)
+      , m_value(traits.m_value)
+      , m_coord_iterator(traits.m_traits.construct_cartesian_const_iterator_3_object())
+      , m_point_3(traits.m_traits.construct_point_3_object())
+      , m_source_3(traits.m_traits.construct_source_3_object())
+      , m_target_3(traits.m_traits.construct_target_3_object())
+    {}
+
+    result_type
+    operator()( const typename Traits::Plane_3&,
+                const typename Traits::Segment_3& s) const
+    {
+      const typename Traits::Point_3& src = m_source_3(s);
+      const typename Traits::Point_3& tgt = m_target_3(s);
+
+      cpp11::array<FT,3> src_coords = {{ *m_coord_iterator(src),
+                                         *(m_coord_iterator(src)+1),
+                                         *(m_coord_iterator(src)+2) }};
+      cpp11::array<FT,3> tgt_coords = {{ *m_coord_iterator(tgt),
+                                         *(m_coord_iterator(tgt)+1),
+                                         *(m_coord_iterator(tgt)+2) }};
+
+      FT alpha = ( m_value - src_coords[m_cst_coord] ) / ( tgt_coords[m_cst_coord] - src_coords[m_cst_coord] );
+      src_coords[m_cst_coord]=m_value;
+      for (int i=1;i<3;++i)
+      {
+        int index = (m_cst_coord+i)%3;
+        src_coords[index]+=(tgt_coords[index]-src_coords[index])*alpha;
+      }
+
+      return Variant_type( m_point_3(src_coords[0], src_coords[1], src_coords[2]) );
+    }
+  };
+
+  Oriented_side_3 oriented_side_3_object() const
+  {
+    return Oriented_side_3(*this);
+  }
+
+  Do_intersect_3 do_intersect_3_object() const
+  {
+    return Do_intersect_3(m_cst_coord, m_value);
+  }
+
+  Intersect_3 intersect_3_object() const
+  {
+    return Intersect_3(*this);
+  }
+};
+
+} } // end of namespace CGAL::Polygon_mesh_slicer
+
+
+#endif // CGAL_INTERNAL_POLYGON_MESH_SLICER_3_AXIS_PARALLEL_PLANE_TRAITS_H

Polygon_mesh_processing/include/CGAL/internal/Polygon_mesh_slicer_3/Traversal_traits.h

@@ -37,22 +37,25 @@ template <typename AL_graph,
           class Traits>
 class Traversal_traits
 {
-///typedefs
+/// typedefs
   typedef typename boost::graph_traits<TriangleMesh>::vertex_descriptor vertex_descriptor;
   typedef typename boost::graph_traits<TriangleMesh>::edge_descriptor     edge_descriptor;
   typedef typename AL_graph::vertex_descriptor                       AL_vertex_descriptor;
   typedef std::pair<const vertex_descriptor, AL_vertex_descriptor>            Vertex_pair;
   typedef std::map<vertex_descriptor, AL_vertex_descriptor>                  Vertices_map;
-///container filled by `intersection()`
+/// container filled by `intersection()`
   std::set<edge_descriptor>& m_all_coplanar_edges;
   Vertices_map& m_vertices;
   std::vector<edge_descriptor>& m_iedges;
-///data members
+/// data members
   TriangleMesh& m_tmesh;
   const VertexPointPmap& m_vpmap;
   const AABBTraits& m_aabb_traits;
   const Traits& m_traits;
   const typename AL_graph::vertex_descriptor null_vertex;
+/// predicates
+  typename Traits::Oriented_side_3 oriented_side_3;
+  typename Traits::Do_intersect_3 do_intersect_3;
 
 public:
 
@@ -71,17 +74,18 @@ public:
     , m_aabb_traits(aabb_traits)
     , m_traits(traits)
     , null_vertex( boost::graph_traits<AL_graph>::null_vertex() )
+    , oriented_side_3( m_traits.oriented_side_3_object() )
+    , do_intersect_3( m_traits.do_intersect_3_object() )
   {}
 
   bool go_further() const { return true; }
 
   void intersection(const typename Traits::Plane_3& plane, const typename AABBTraits::Primitive& primitive)
   {
-    typename Traits::Oriented_side_3 oriented_side = m_traits.oriented_side_3_object();
     typename boost::graph_traits<TriangleMesh>::edge_descriptor ed = primitive.id();
 
-    Oriented_side src = oriented_side(plane, get(m_vpmap, source(ed,m_tmesh)) );
-    Oriented_side tgt = oriented_side(plane, get(m_vpmap, target(ed,m_tmesh)) );
+    Oriented_side src = oriented_side_3(plane, get(m_vpmap, source(ed,m_tmesh)) );
+    Oriented_side tgt = oriented_side_3(plane, get(m_vpmap, target(ed,m_tmesh)) );
 
     if (src==ON_ORIENTED_BOUNDARY)
     {
@@ -99,10 +103,10 @@ public:
     }
   }
 
-  template<class Query,class Node>
-  bool do_intersect(const Query& query, const Node& node) const
+  template<class Node>
+  bool do_intersect(const typename Traits::Plane_3& plane, const Node& node) const
   {
-    return m_aabb_traits.do_intersect_object()(query, node.bbox());
+    return do_intersect_3(plane, node.bbox());
   }
 };