Add a robust unfolding traits class

This is an option for getting around the lack of robustness when
unfolding adjacent faces, without the full overhead of an exact precision
kernel.

Author: Hippocrates

Polyhedron_shortest_path/doc/Polyhedron_shortest_path/Polyhedron_shortest_path.txt

@@ -89,6 +89,8 @@ In order to use the algorithm, one creates an instance of the `Polyhedron_shorte
 \subsection Polyhedron_shortest_pathUnfoldingMethods Face Unfolding/Layout Methods
   In order to effectively compute the distance along the surface, it is necessary to unfold sequences of faces, edge-to-edge, out into a common plane.  The functor `Project_triangle_3_to_triangle_2` provides an initial layout of the first face in a sequence, by rotating a given face into the XY plane.  `Flatten_triangle_3_along_segment_2` unfolds a triangle into the plane, using a specified segment as a base.  Since this results in a chain of constructed triangles in the plane, exact representation types such as CORE::Expr or leda::real will process extremely slow (effectively adding an \f$O(n)\f$ factor to any computation) and therefore their use with this algorithm is not recommended.
 
+  If more robustness is needed, consider using the `Polyhedron_shortest_path_default_traits_with_robust_unfolding` traits class; it implements the unfolding methods using an exact kernel in isolation, avoiding intermediate roundoff error within each operation, while still keeping the overall algorithm from becoming too slow (however it will still take much longer than with an inexact kernel).
+  
 \subsection Polyhedron_shortest_pathBarycentric Barycentric Coordinates and Face Locations
 
   The source locations of the algorithm are specified either using vertices (specifically `boost::GraphTraits<FaceGraph>::vertex_descriptor>`, or face/location pairs, where the location is represented using a <em>Barycentric Coordinate</em> in the face. Given a point \f$p\f$ that lies inside and on the plane of triangle \f$(A,B,C)\f$, its barycentric coordinate is a real-valued triple \f$(b_0,b_1,b_2)\f$ such that \f$p = b_0\cdot~A + b_1\cdot~B + b_2\cdot~C\f$. 

Polyhedron_shortest_path/include/CGAL/Polyhedron_shortest_path/Polyhedron_shortest_path_traits.h

@@ -27,7 +27,7 @@ model as required by the Polyhedron_shortest_path algorithm
 
 \tparam F The faceGraph type the algorithm is to act on
 
-\cgalModels `FaceGraphShortestPathTraits`
+\cgalModels `PolyhedronShortestPathTraits`
 */
 template <
   class K, 
@@ -155,6 +155,53 @@ std::ostream& operator<<(std::ostream& os, typename Polyhedron_shortest_path_def
   return os << b[0] << " " << b[1] << " " << b[2];
 }
 
+
+
+/*!
+\ingroup PkgPolyhedronShortestPathTraitsClasses
+
+\brief Provides an implementation of the FaceGraphShortestPathTraits 
+model as required by the Polyhedron_shortest_path algorithm
+
+\tparam K The kernel type whose geometric primitives to use
+
+\tparam F The faceGraph type the algorithm is to act on
+
+\cgalModels `PolyhedronShortestPathTraits`
+*/
+template <
+  class K, 
+  class F>
+class Polyhedron_shortest_path_default_traits_with_robust_unfolding : public Polyhedron_shortest_path_default_traits<K,F>
+{
+public:
+  typedef K Kernel;
+  typedef typename PolyhedronShortestPath::Robust_project_triangle_3_to_triangle_2<K> Project_triangle_3_to_triangle_2;
+  typedef typename PolyhedronShortestPath::Robust_flatten_triangle_3_along_segment_2<K> Flatten_triangle_3_along_segment_2;
+  
+private:
+
+  Project_triangle_3_to_triangle_2 m_robust_project_triangle_3_to_triangle_2_object;
+  Flatten_triangle_3_along_segment_2 m_robust_flatten_triangle_3_along_segment_2;
+  
+public:
+
+  Polyhedron_shortest_path_default_traits_with_robust_unfolding()
+  {
+  }
+  
+  Polyhedron_shortest_path_default_traits_with_robust_unfolding(const Kernel& kernel)
+    : Polyhedron_shortest_path_default_traits<K,F>(kernel)
+    , m_robust_project_triangle_3_to_triangle_2_object(kernel)
+    , m_robust_flatten_triangle_3_along_segment_2(kernel)
+  {
+  }
+  
+  Project_triangle_3_to_triangle_2 project_triangle_3_to_triangle_2_object() const { return m_robust_project_triangle_3_to_triangle_2_object; }
+  Flatten_triangle_3_along_segment_2 flatten_triangle_3_along_segment_2_object() const { return m_robust_flatten_triangle_3_along_segment_2; }
+};
+
 } // namespace CGAL
 
+
 #endif // CGAL_POLYHEDRON_SHORTEST_PATH_TRAITS_H

Polyhedron_shortest_path/include/CGAL/Polyhedron_shortest_path/function_objects.h

@@ -13,6 +13,8 @@
 
 #include <CGAL/Polyhedron_shortest_path/internal/misc_functions.h>
 
+#include <CGAL/Exact_predicates_exact_constructions_kernel_with_sqrt.h>
+
 #ifndef CGAL_POLYHEDRON_SHORTEST_PATH_INTERNAL_FUNCTION_OBJECTS_H
 #define CGAL_POLYHEDRON_SHORTEST_PATH_INTERNAL_FUNCTION_OBJECTS_H
 
@@ -194,7 +196,7 @@ class Project_triangle_3_to_triangle_2
 private:
   Parametric_distance_along_segment_3<K> m_parametric_distance_along_segment_3;
   Compute_squared_distance_3 m_compute_squared_distance_3;
-  Construct_line_3 m_constuct_line_3;
+  Construct_line_3 m_construct_line_3;
   Construct_projected_point_3 m_construct_projected_point_3;
   mutable Construct_vertex_3 m_construct_vertex_3;
   Construct_vector_3 m_construct_vector_3;
@@ -208,7 +210,7 @@ class Project_triangle_3_to_triangle_2
 
   Project_triangle_3_to_triangle_2(const K& kernel)
     : m_compute_squared_distance_3(kernel.compute_squared_distance_3_object())
-    , m_constuct_line_3(kernel.constuct_line_3_object())
+    , m_construct_line_3(kernel.construct_line_3_object())
     , m_construct_vertex_3(kernel.construct_vertex_3_object())
     , m_construct_projected_point_3(kernel.construct_projected_point_3_object())
     , m_construct_vector_3(kernel.construct_vector_3_object())
@@ -222,7 +224,7 @@ class Project_triangle_3_to_triangle_2
     Vector_3 v01(m_construct_vector_3(m_construct_vertex_3(t3, 0), m_construct_vertex_3(t3, 1)));
     Vector_3 v02(m_construct_vector_3(m_construct_vertex_3(t3, 0), m_construct_vertex_3(t3, 2)));
 
-    Line_3 baseSegment(m_constuct_line_3(m_construct_vertex_3(t3, 0), m_construct_vertex_3(t3, 1)));
+    Line_3 baseSegment(m_construct_line_3(m_construct_vertex_3(t3, 0), m_construct_vertex_3(t3, 1)));
 
     //FT scalePoint = (v01 * v02) / (v01 * v01);
     Point_3 projectedLocation3d(m_construct_projected_point_3(baseSegment, m_construct_vertex_3(t3, 2)));
@@ -238,6 +240,36 @@ class Project_triangle_3_to_triangle_2
   }
 };
 
+template<class K>
+class Robust_project_triangle_3_to_triangle_2
+{
+public:
+  typedef typename K::Triangle_3 Triangle_3;
+  typedef typename K::Triangle_2 Triangle_2;
+  typedef Exact_predicates_exact_constructions_kernel_with_sqrt EKSQRT;
+  typedef Project_triangle_3_to_triangle_2<EKSQRT> Exact_project_triangle_3_to_triangle_2;
+  typedef Cartesian_converter<K, EKSQRT>  To_exact;
+  typedef Cartesian_converter<EKSQRT, K>  Back_from_exact;
+
+public:
+  Robust_project_triangle_3_to_triangle_2()
+  {
+  }
+  
+  Robust_project_triangle_3_to_triangle_2(const K& kernel)
+  {
+  }
+
+  Triangle_2 operator() (const Triangle_3& t3) const
+  {
+    Exact_project_triangle_3_to_triangle_2 ept3t2;
+    To_exact to_exact;
+    Back_from_exact back_from_exact;
+    
+    return back_from_exact(ept3t2(to_exact(t3)));
+  }
+};
+
 template<class K>
 class Flatten_triangle_3_along_segment_2
 {
@@ -273,7 +305,7 @@ class Flatten_triangle_3_along_segment_2
 
   Parametric_distance_along_segment_3<K> m_parametric_distance_along_segment_3;
   Compute_squared_distance_3 m_compute_squared_distance_3;
-  Construct_projected_point_3 m_constrct_projected_point_3;
+  Construct_projected_point_3 m_construct_projected_point_3;
   Construct_perpendicular_vector_2 m_construct_perpendicular_vector_2;
   Construct_sum_of_vectors_2 m_construct_sum_of_vectors_2;
   Construct_scaled_vector_2 m_construct_scaled_vector_2;
@@ -294,7 +326,7 @@ class Flatten_triangle_3_along_segment_2
 
   Flatten_triangle_3_along_segment_2(const K& kernel)
     : m_compute_squared_distance_3(kernel.compute_squared_distance_3_object())
-    , m_constrct_projected_point_3(kernel.constrct_projected_point_3_object())
+    , m_construct_projected_point_3(kernel.construct_projected_point_3_object())
     , m_construct_perpendicular_vector_2(kernel.construct_perpendicular_vector_2_object())
     , m_construct_sum_of_vectors_2(kernel.construct_sum_of_vectors_2_object())
     , m_construct_scaled_vector_2(kernel.construct_scaled_vector_2_object())
@@ -312,7 +344,7 @@ class Flatten_triangle_3_along_segment_2
 
   result_type operator() (const Triangle_3& t3, size_t edgeIndex, const Segment_2& segment) const
   {
-    Point_3 projectedLocation3d(m_constrct_projected_point_3(m_construct_line_3(m_construct_vertex_3(t3, edgeIndex), m_construct_vertex_3(t3, edgeIndex + 1)), m_construct_vertex_3(t3, edgeIndex + 2)));
+    Point_3 projectedLocation3d(m_construct_projected_point_3(m_construct_line_3(m_construct_vertex_3(t3, edgeIndex), m_construct_vertex_3(t3, edgeIndex + 1)), m_construct_vertex_3(t3, edgeIndex + 2)));
     FT scalePoint = m_parametric_distance_along_segment_3(m_construct_segment_3(m_construct_vertex_3(t3, edgeIndex), m_construct_vertex_3(t3, edgeIndex + 1)), projectedLocation3d);
     FT triangleHeight = CGAL::internal::select_sqrt(m_compute_squared_distance_3(projectedLocation3d, m_construct_vertex_3(t3, edgeIndex + 2)));
 
@@ -331,6 +363,38 @@ class Flatten_triangle_3_along_segment_2
   }
 };
 
+template<class K>
+class Robust_flatten_triangle_3_along_segment_2
+{
+public:
+  typedef typename K::Triangle_3 Triangle_3;
+  typedef typename K::Segment_2 Segment_2;
+  typedef typename K::Triangle_2 Triangle_2;
+  
+  typedef Exact_predicates_exact_constructions_kernel_with_sqrt EKSQRT;
+  typedef Flatten_triangle_3_along_segment_2<EKSQRT> Exact_flatten_triangle_3_along_segment_2;
+  typedef Cartesian_converter<K, EKSQRT>  To_exact;
+  typedef Cartesian_converter<EKSQRT, K>  Back_from_exact;
+
+public:
+  Robust_flatten_triangle_3_along_segment_2()
+  {
+  }
+  
+  Robust_flatten_triangle_3_along_segment_2(const K& kernel)
+  {
+  }
+
+  Triangle_2 operator() (const Triangle_3& t3, size_t edgeIndex, const Segment_2& segment) const
+  {
+    Exact_flatten_triangle_3_along_segment_2 eft3as2;
+    To_exact to_exact;
+    Back_from_exact back_from_exact;
+    
+    return back_from_exact(eft3as2(to_exact(t3), edgeIndex, to_exact(segment)));
+  }
+};
+
 template <class K>
 class Compare_relative_intersection_along_segment_2
 {
@@ -473,7 +537,7 @@ class Is_saddle_vertex
   {
   }
 
-  Is_saddle_vertex(const Kernel& kernel, const Project_triangle_3_to_triangle_2& pt3tt2, const Flatten_triangle_3_along_segment_2& ft3as2, const Orientation_2& o2)
+  Is_saddle_vertex(const Kernel& kernel, const Project_triangle_3_to_triangle_2& pt3tt2, const Flatten_triangle_3_along_segment_2& ft3as2)
     : m_orientation_2(kernel.orientation_2_object())
     , m_construct_triangle_3(kernel.construct_triangle_3_object())
     , m_construct_vertex_2(kernel.construct_vertex_2_object())