Add a geom-traits concept for distance functions

Update the code and the doc accordingly

AABB_tree/doc/AABB_tree/Concepts/AABBGeomTraits.h

@@ -65,7 +65,6 @@ A functor object to construct the sphere centered at one point and passing throu
 typedef unspecified_type Construct_sphere_3; 
 
 /*!
-\todo This is not correct! that is not used!
 A functor object to compute the point on a geometric primitive which is closest from a query. Provides the operator: 
 `Point_3 operator()(const Type_2& type_2, const Point_3& p);` where `Type_2` is any type among `Segment_3` and `Triangle_3`. The operator returns the point on `type_2` which is closest to `p`. 
 */ 
@@ -76,7 +75,7 @@ A functor object to compare the distance of two points wrt a third one.
 Provides the operator:
 `CGAL::Comparision_result operator()(const Point_3& p1, const Point_3& p2, const Point_3& p3)`. The operator compare the distance between `p1 and `p2`, and between `p2` and `p3`.
 */
-typedef unspecified_type Compare_distance_3
+typedef unspecified_type Compare_distance_3;
 
 /*!
 A functor object to detect if a point lies inside a sphere or not. 

Poisson_surface_reconstruction_3/examples/Poisson_surface_reconstruction_3/poisson_reconstruction_example.cpp

@@ -109,7 +109,7 @@ int main(void)
 
     // computes the approximation error of the reconstruction
     double max_dist =
-      CGAL::Polygon_mesh_processing::max_distance_to_point_set(output_mesh,
+      CGAL::Polygon_mesh_processing::approximate_max_distance_to_point_set(output_mesh,
                                                                points,
                                                                4000);
     std::cout << "Max distance to point_set: " << max_dist << std::endl;

Polygon_mesh_processing/doc/Polygon_mesh_processing/Concepts/PMPDistanceTraits.h

@@ -0,0 +1,73 @@
+/// \ingroup PkgPolygonMeshProcessingConcepts
+/// \cgalConcept
+///
+/// The concept `PMPDistanceTraits` is a refinement of the
+/// concepts `AABBGeomTraits` and `SpatialSortingTraits_3`. In addition to the types required by
+/// those concepts,  it also requires types and functors needed by the functions `approximate_max_distance_to_point_set()`,
+/// `sample_triangle_mesh()`, `approximate_Hausdorff_distance()` and `max_distance_to_triangle_mesh()`
+///
+/// \cgalRefines `AABBGeomTraits`
+/// \cgalRefines `SpatialSortingTraits_3`
+/// \cgalHasModel Any 3D Kernel is a model of this concept.
+
+class PMPDistanceTraits{
+public:
+    /*!
+     * A number type model of `Field` and `RealEmbeddable`
+     */
+    typedef unspecified_type FT;
+    ///
+    /*! 3D point type
+     * It must be default constructible, and can be constructed from 3 objects of type `FT`.
+     * `bool operator<(Point_3, Point_3)` to lexicographically compare two points must be available.
+     * Access to Cartesian coordinates must be possible using `Point_3::x()`, `Point_3::y(), Point_3::z()` and
+     * `FT operator[](int i)` with  `0 <= i < 3`.
+     *
+     * There must be a specialization of `CGAL::Kernel_traits` such that
+     * `CGAL::Kernel_traits<Point_3>::Kernel` is the model implementing this concept.
+     */
+    typedef unspecified_type Point_3;
+
+    /// 3D vector type
+    typedef unspecified_type Vector_3;
+
+    /*!
+     * 3D triangle type
+     * It must be constructible from three points and access to its vertices
+     * must be possible with operator `const Point_3& operator[](int i) const`
+     * with  `0 <= i < 3`.
+     * In addition, it must provide `BBox_3 bbox() const` that returns a
+     * bounding box of the triangle.
+     */
+    typedef unspecified_type Triangle_3;
+
+    /// @name Functors
+    /// @{
+
+    /// Functor for computing squared area of a triangle.
+    /// It provides `FT operator()(const Point_3&, const Point_3&, const Point_3&) const`
+    /// and `FT operator()(const Triangle_3&) const` and has `FT` as result_type.
+    typedef unspecified_type Compute_squared_area_3;
+    /// Functor for constructing translated points.
+    /// It provides `Point_3 operator()(const Point_3 &, const Vector_3 &)`
+    typedef unspecified_type Construct_translated_point_3;
+    /// Functor for constructing vectors.
+    /// It provides `Vector_3 operator()(const Point_3 &, const Point_3 &)`
+    typedef unspecified_type Construct_vector_3;
+    /// Functor for constructing scaled vectors.
+    /// It provides `Vector_3 operator()(const Vector_3 &, const FT &)`
+    typedef unspecified_type Construct_scaled_vector_3;
+    /// Functor for constructing the circumcenter of three points.
+    /// It provides `Point_3 operator()(const Point_3&, const Point_3&, const Point_3&)`
+    typedef unspecified_type Construct_circumcenter_3;
+    /// @}
+
+    /// @name Functions
+    /// @{
+    Compute_squared_area_3 compute_squared_area_3_object();
+    Construct_translated_point_3 construct_translated_point_3_object();
+    Construct_vector_3 construct_vector_3_object();
+    Construct_scaled_vector_3 construct_scaled_vector_3_object();
+    Construct_circumcenter_3 construct_circumcenter_3_object();
+    /// @}
+};

Polygon_mesh_processing/doc/Polygon_mesh_processing/Concepts/PMPSelfIntersectionTraits.h

@@ -2,7 +2,7 @@
 /// \cgalConcept
 ///
 /// Geometric traits concept for the functions `CGAL::self_intersections()` and `CGAL::does_self_intersect()`.
-class SelfIntersectionTraits{
+class PMPSelfIntersectionTraits{
 public:
   /// @name Geometric Types
   /// @{

Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt

@@ -133,9 +133,9 @@ and provides a list of the parameters that are used in this package.
 - \link measure_grp `CGAL::Polygon_mesh_processing::face_border_length()` \endlink
 
 ## Distance Functions ##
-- `CGAL::Polygon_mesh_processing::approximated_Hausdorff_distance()`
-- `CGAL::Polygon_mesh_processing::approximated_symmetric_Hausdorff_distance()`
-- `CGAL::Polygon_mesh_processing::max_distance_to_point_set()`
+- `CGAL::Polygon_mesh_processing::approximate_Hausdorff_distance()`
+- `CGAL::Polygon_mesh_processing::approximate_symmetric_Hausdorff_distance()`
+- `CGAL::Polygon_mesh_processing::approximate_max_distance_to_point_set()`
 - `CGAL::Polygon_mesh_processing::max_distance_to_triangle_mesh()`
 - `CGAL::Polygon_mesh_processing::sample_triangle_mesh()`
 

Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt

@@ -486,9 +486,9 @@ the propagation of a connected component index to cross it.
 This package provides methods to compute approximated Hausdorff distances.
 These distances can be computed between two meshes, a mesh and a point set, or a point set and a mesh.
 These computations are performed with :
-- `CGAL::Polygon_mesh_processing::approximated_Hausdorff_distance()`
-- `CGAL::Polygon_mesh_processing::approximated_symmetric_Hausdorff_distance()`
-- `CGAL::Polygon_mesh_processing::max_distance_to_point_set()`
+- `CGAL::Polygon_mesh_processing::approximate_Hausdorff_distance()`
+- `CGAL::Polygon_mesh_processing::approximate_symmetric_Hausdorff_distance()`
+- `CGAL::Polygon_mesh_processing::approximate_max_distance_to_point_set()`
 - `CGAL::Polygon_mesh_processing::max_distance_to_triangle_mesh()`
 - `CGAL::Polygon_mesh_processing::sample_triangle_mesh()`
 

Polygon_mesh_processing/doc/Polygon_mesh_processing/dependencies

@@ -11,3 +11,4 @@ Surface_mesh
 Surface_mesh_deformation
 AABB_tree
 Triangulation_2
+Spatial_sorting

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h

@@ -55,8 +55,9 @@ triangle_grid_sampling( const typename Kernel::Point_3& p0,
                         double distance,
                         OutputIterator out)
 {
-  const double d_p0p1 = CGAL::sqrt( CGAL::squared_distance(p0, p1) );
-  const double d_p0p2 = CGAL::sqrt( CGAL::squared_distance(p0, p2) );
+    typename Kernel::Compute_squared_distance_3 squared_distance;
+  const double d_p0p1 = to_double(approximate_sqrt( squared_distance(p0, p1) ));
+  const double d_p0p2 = to_double(approximate_sqrt( squared_distance(p0, p2) ));
 
   const double n = (std::max)(std::ceil( d_p0p1 / distance ),
                               std::ceil( d_p0p2 / distance ));
@@ -102,13 +103,14 @@ sample_triangles(const TriangleRange& triangles, double distance, OutputIterator
       const Point_3& p1=t[(i+1)%3];
       if ( sampled_edges.insert(CGAL::make_sorted_pair(p0, p1)).second )
       {
-        const double d_p0p1 = CGAL::sqrt( CGAL::squared_distance(p0, p1) );
+          typename Kernel::Compute_squared_distance_3 squared_distance;
+        const double d_p0p1 = to_double(approximate_sqrt( squared_distance(p0, p1) ));
 
         const double nb_pts = std::ceil( d_p0p1 / distance );
-        const Vector_3 step_vec = (p1 - p0) / nb_pts;
+        const Vector_3 step_vec =  typename Kernel::Construct_scaled_vector_3()(typename Kernel::Construct_vector_3()(p1, p0), typename Kernel::FT(1)/typename Kernel::FT(nb_pts));
         for (double i=1; i<nb_pts; ++i)
         {
-          *out++=p0 + step_vec * i;
+          *out++=typename Kernel::Construct_translated_point_3()(p0, typename Kernel::Construct_scaled_vector_3()(step_vec , typename Kernel::FT(i)));
         }
       }
     }
@@ -142,8 +144,8 @@ struct Distance_computation{
           cpp11::atomic<double>* d)
     : tree(tree)
     , sample_points(sample_points)
-    , distance(d)
     , initial_hint(p)
+    , distance(d)
   {}
 
   void
@@ -154,7 +156,8 @@ struct Distance_computation{
     for( std::size_t i = range.begin(); i != range.end(); ++i)
     {
       hint = tree.closest_point(sample_points[i], hint);
-      double d = CGAL::approximate_sqrt( squared_distance(hint,sample_points[i]) );
+      typename Kernel_traits<Point_3>::Kernel::Compute_squared_distance_3 squared_distance;
+      double d = to_double(CGAL::approximate_sqrt( squared_distance(hint,sample_points[i]) ));
       if (d>hdist) hdist=d;
     }
 
@@ -198,7 +201,7 @@ enum Sampling_method{
  *    \cgalParamBegin{vertex_point_map}
  *    the property map with the points associated to the vertices of `tm`. If this parameter is omitted,
  *    an internal property map for `CGAL::vertex_point_t` should be available for `TriangleMesh` \cgalParamEnd
- *    \cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `Kernel`\cgalParamEnd
+ *    \cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `PMPDistanceTraits`\cgalParamEnd
  * \cgalNamedParamsEnd
  */
 template<class OutputIterator, class TriangleMesh, class NamedParameters>
@@ -217,56 +220,59 @@ sample_triangle_mesh(const TriangleMesh& tm,
 
     Vpm pmap = choose_param(get_param(np, vertex_point),
                            get_const_property_map(vertex_point, tm));
+    typedef Creator_uniform_3<typename Geom_traits::FT,
+                              typename Geom_traits::Point_3> Creator;
   switch(method)
   {
-  case RANDOM_UNIFORM:
-  {
-    std::size_t nb_points = std::ceil(
-       parameter * area(tm, parameters::geom_traits(Geom_traits())));
-    Random_points_in_triangle_mesh_3<TriangleMesh, Vpm>
-        g(tm, pmap);
-    CGAL::cpp11::copy_n(g, nb_points, out);
-    return out;
-  }
-  case GRID:
-  {
-    std::vector<typename Geom_traits::Triangle_3> triangles;
-    BOOST_FOREACH(typename boost::graph_traits<TriangleMesh>::face_descriptor f, faces(tm))
+    case RANDOM_UNIFORM:
     {
-      //create the triangles and store them
-      typename Geom_traits::Point_3 points[3];
-      typename boost::graph_traits<TriangleMesh>::halfedge_descriptor hd(halfedge(f,tm));
-      for(int i=0; i<3; ++i)
-      {
-        points[i] = get(pmap, target(hd, tm));
-        hd = next(hd, tm);
-      }
-      triangles.push_back(typename Geom_traits::Triangle_3(points[0], points[1], points[2]));
-    }
-    sample_triangles<Geom_traits>(triangles, parameter, out);
-    return out;
-  }
-  case MONTE_CARLO:
-  {
-    BOOST_FOREACH(typename boost::graph_traits<TriangleMesh>::face_descriptor f, faces(tm))
-    {
-      std::size_t nb_points =  (std::max)(
-                  (int)std::ceil(parameter * face_area(f,tm,parameters::geom_traits(Geom_traits()))),
-                                       1);
-      //create the triangles and store them
-      typename Geom_traits::Point_3 points[3];
-      typename boost::graph_traits<TriangleMesh>::halfedge_descriptor hd(halfedge(f,tm));
-      for(int i=0; i<3; ++i)
-      {
-        points[i] = get(pmap, target(hd, tm));
-        hd = next(hd, tm);
-      }
-      Random_points_in_triangle_3<typename Geom_traits::Point_3> g(points[0], points[1], points[2]);
+      std::size_t nb_points = std::ceil( to_double(
+         area(tm, parameters::geom_traits(Geom_traits()))*parameter) );
+      Random_points_in_triangle_mesh_3<TriangleMesh, Vpm, Creator>
+          g(tm, pmap);
       CGAL::cpp11::copy_n(g, nb_points, out);
+      return out;
+    }
+    case GRID:
+    {
+      std::vector<typename Geom_traits::Triangle_3> triangles;
+      BOOST_FOREACH(typename boost::graph_traits<TriangleMesh>::face_descriptor f, faces(tm))
+      {
+        //create the triangles and store them
+        typename Geom_traits::Point_3 points[3];
+        typename boost::graph_traits<TriangleMesh>::halfedge_descriptor hd(halfedge(f,tm));
+        for(int i=0; i<3; ++i)
+        {
+          points[i] = get(pmap, target(hd, tm));
+          hd = next(hd, tm);
+        }
+        triangles.push_back(typename Geom_traits::Triangle_3(points[0], points[1], points[2]));
+      }
+      sample_triangles<Geom_traits>(triangles, parameter, out);
+      return out;
+    }
+    case MONTE_CARLO:
+    {
+      BOOST_FOREACH(typename boost::graph_traits<TriangleMesh>::face_descriptor f, faces(tm))
+      {
+        std::size_t nb_points( (std::max)(
+                    std::ceil(to_double(face_area(f,tm,parameters::geom_traits(Geom_traits())))*parameter),
+                                         1.) );
+        //create the triangles and store them
+        typename Geom_traits::Point_3 points[3];
+        typename boost::graph_traits<TriangleMesh>::halfedge_descriptor hd(halfedge(f,tm));
+        for(int i=0; i<3; ++i)
+        {
+          points[i] = get(pmap, target(hd, tm));
+          hd = next(hd, tm);
+        }
+        Random_points_in_triangle_3<typename Geom_traits::Point_3, Creator>
+          g(points[0], points[1], points[2]);
+        CGAL::cpp11::copy_n(g, nb_points, out);
+      }
     }
-    return out;
-  }
   }
+  return out;
 }
 
 template<class OutputIterator, class TriangleMesh>
@@ -328,6 +334,7 @@ double approximate_Hausdorff_distance(
     BOOST_FOREACH(const typename Kernel::Point_3& pt, sample_points)
     {
       hint = tree.closest_point(pt, hint);
+      typename Kernel::Compute_squared_distance_3 squared_distance;
       typename Kernel::FT dist = squared_distance(hint,pt);
       double d = to_double(CGAL::approximate_sqrt(dist));
       if(d>hdist)
@@ -390,7 +397,7 @@ double approximate_Hausdorff_distance(
  *    \cgalParamBegin{vertex_point_map}
  *    the property map with the points associated to the vertices of `tm1` (`tm2`). If this parameter is omitted,
  *    an internal property map for `CGAL::vertex_point_t` should be available for `TriangleMesh` \cgalParamEnd
- *    \cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `Kernel`\cgalParamEnd
+ *    \cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `PMPDistanceTraits`\cgalParamEnd
  * \cgalNamedParamsEnd
  */
 template< class Concurrency_tag,
@@ -458,7 +465,7 @@ double approximate_symmetric_Hausdorff_distance(
  *    the property map with the points associated to the vertices of `tm`. If this parameter is omitted,
  *    an internal property map for `CGAL::vertex_point_t` should be available for the
       vertices of `tm` \cgalParamEnd
- *    \cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `Kernel`\cgalParamEnd
+ *    \cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `PMPDistanceTraits`\cgalParamEnd
  * \cgalNamedParamsEnd
  */
 template< class Concurrency_tag,
@@ -494,13 +501,13 @@ double max_distance_to_triangle_mesh(const PointRange& points,
  *    the property map with the points associated to the vertices of `tm`. If this parameter is omitted,
  *    an internal property map for `CGAL::vertex_point_t` should be available for the
       vertices of `tm` \cgalParamEnd
- *    \cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `Kernel`. \cgalParamEnd
+ *    \cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `PMPDistanceTraits`. \cgalParamEnd
  * \cgalNamedParamsEnd
  */
 template< class TriangleMesh,
           class PointRange,
           class NamedParameters>
-double max_distance_to_point_set(const TriangleMesh& tm,
+double approximate_max_distance_to_point_set(const TriangleMesh& tm,
                                  const PointRange& points,
                                  const double precision,
                                  const NamedParameters& np)
@@ -525,7 +532,7 @@ double max_distance_to_point_set(const TriangleMesh& tm,
     }
     ref.add(points[0], points[1], points[2], tree);
   }
-  return ref.refine(precision, tree);
+  return to_double(ref.refine(precision, tree));
 }
 
 // convenience functions with default parameters
@@ -544,11 +551,11 @@ double max_distance_to_triangle_mesh(const PointRange& points,
 
 template< class TriangleMesh,
           class PointRange>
-double max_distance_to_point_set(const TriangleMesh& tm,
+double approximate_max_distance_to_point_set(const TriangleMesh& tm,
                                  const PointRange& points,
                                  const double precision)
 {
-  return max_distance_to_point_set(tm, points, precision, parameters::all_default());
+  return approximate_max_distance_to_point_set(tm, points, precision, parameters::all_default());
 }
 
 template< class Concurrency_tag,

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/mesh_to_point_set_hausdorff_distance.h

@@ -25,23 +25,21 @@
 #include <CGAL/Search_traits_3.h>
 #include <CGAL/squared_distance_3.h>
 #include <queue>
+#include <iostream>
+
 namespace CGAL{
 namespace internal{
 template <typename Kernel>
 class CPointH
 {
-private:
-
-  typedef typename Kernel::Point_3 Point;
-  typedef typename Kernel::FT FT;
-  Point m_point;
-  FT m_hausdorff;
 
 public:
+    typedef typename Kernel::Point_3 Point;
+    typedef typename Kernel::FT FT;
 
   CPointH ()
   {
-    m_point = Point (0., 0., 0.);
+    m_point = Point (FT(0.), FT(0.), FT(0.));
     m_hausdorff = 0.;
   }
   CPointH (const Point& p, FT h = 0.)
@@ -67,6 +65,10 @@ public:
   }
 
 
+private:
+
+    Point m_point;
+    FT m_hausdorff;
 };
 
 template <typename Kernel>
@@ -75,9 +77,9 @@ class CRefTriangle
 
 private:
 
-  typedef typename Kernel::Point_3 Point;
   typedef typename Kernel::FT FT;
   typedef CPointH<Kernel> PointH;
+  typedef typename PointH::Point Point;
   PointH m_point[3];
   int m_edge;
   FT m_upper_bound;
@@ -93,9 +95,10 @@ public:
     m_point[2] = c;
 
     m_edge = -1;
-    for (unsigned int i = 0; i < 3; ++ i)
+    //should disappear when Simon reworks the code
+    /*for (unsigned int i = 0; i < 3; ++ i)
     {
-      if (CGAL::angle (m_point[(i+1)%3](), m_point[i](), m_point[(i+2)%3]())
+      if (CGAL::angle(m_point[(i+1)%3](), m_point[i](), m_point[(i+2)%3]())
           == CGAL::OBTUSE)
       {
         m_edge = i;
@@ -103,14 +106,14 @@ public:
       }
     }
 
+*/
     if (m_edge == -1)
-      m_bisector = CGAL::circumcenter (a(), b(), c());
+        m_bisector = typename Kernel::Construct_circumcenter_3() (a(), b(), c());
     else
     {
       m_bisector = PointH::mid_point (m_point[(m_edge+1)%3],
                                       m_point[(m_edge+2)%3]);
     }
-
     m_lower_bound = 0.;
     m_upper_bound = 0.;
     for (unsigned int i = 0; i < 3; ++ i)
@@ -118,8 +121,9 @@ public:
       if (m_point[i].hausdorff () > m_lower_bound)
         m_lower_bound = m_point[i].hausdorff ();
 
+      typename Kernel::Compute_squared_distance_3 squared_distance;
       FT up = m_point[i].hausdorff ()
-        + std::sqrt (CGAL::squared_distance (m_point[i](), m_bisector));
+        + CGAL::approximate_sqrt (squared_distance (m_point[i](), m_bisector));
 
       if (up > m_upper_bound)
         m_upper_bound = up;
@@ -160,6 +164,7 @@ public:
   const PointH* points () const { return m_point; }
   int edge () const { return m_edge; }
 
+  #ifdef CGAL_MTPS_HD_DEBUG
   void print () const
   {
     std::cerr << "[Refinement triangle]" << std::endl
@@ -170,6 +175,7 @@ public:
               << " * " << m_point[2]() << std::endl
               << " -> " << m_bisector << std::endl;
   }
+  #endif
 };
 }//internal
 
@@ -240,9 +246,9 @@ public:
   bool add (const Point& a, const Point& b, const Point& c, const Tree& tree)
   {
 
-    RefTriangle r (PointH (a, std::sqrt (Knn(tree, a, 1).begin()->second)),
-                   PointH (b, std::sqrt (Knn(tree, b, 1).begin()->second)),
-                   PointH (c, std::sqrt (Knn(tree, c, 1).begin()->second)));
+    RefTriangle r (PointH (a, CGAL::approximate_sqrt (Knn(tree, a, 1).begin()->second)),
+                   PointH (b, CGAL::approximate_sqrt (Knn(tree, b, 1).begin()->second)),
+                   PointH (c, CGAL::approximate_sqrt (Knn(tree, c, 1).begin()->second)));
 
     if (r.lower_bound () > m_lower_bound)
       {
@@ -280,7 +286,9 @@ public:
     {
       if (m_queue.size () > 100000)
       {
+        #ifdef CGAL_MTPS_HD_DEBUG
         m_queue.top ().print ();
+        #endif
         ++ nb_clean;
         if (nb_clean > 5)
           return m_upper_bound;
@@ -295,9 +303,9 @@ public:
 
       if (current.lower_bound () > m_lower_bound)
         m_lower_bound = current.lower_bound ();
+      typename Kernel::Compute_squared_area_3 squared_area;
 
-
-      if(CGAL::squared_area(current.points()[0](),
+      if(squared_area(current.points()[0](),
                             current.points()[1](),
                             current.points()[2]()
                             ) < 1e-20)
@@ -308,7 +316,7 @@ public:
       const Point& bisector = current.bisector ();
       m_point = bisector;
       //squared distance between bisector and its closst point in the mesh
-      FT hausdorff = std::sqrt (Knn(tree, bisector, 1).begin()->second);
+      FT hausdorff = CGAL::approximate_sqrt (Knn(tree, bisector, 1).begin()->second);
       if (hausdorff > m_lower_bound)
         m_lower_bound = hausdorff;
 

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/self_intersections.h

@@ -241,7 +241,7 @@ self_intersections( const FaceRange& face_range,
  *    \cgalParamBegin{vertex_point_map} the property map with the points associated to the vertices of `pmesh`.
  *   If this parameter is omitted, an internal property map for
  *   `CGAL::vertex_point_t` should be available in `TriangleMesh`\cgalParamEnd
- *    \cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `SelfIntersectionTraits` \cgalParamEnd
+ *    \cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `PMPSelfIntersectionTraits` \cgalParamEnd
  * \cgalNamedParamsEnd
  *
  * @return `out`
@@ -300,7 +300,7 @@ self_intersections(const TriangleMesh& tmesh, OutputIterator out)
  *    \cgalParamBegin{vertex_point_map} the property map with the points associated to the vertices of `pmesh`.
  *   If this parameter is omitted, an internal property map for
  *   `CGAL::vertex_point_t` should be available in `TriangleMesh`\cgalParamEnd
- *    \cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `SelfIntersectionTraits` \cgalParamEnd
+ *    \cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `PMPSelfIntersectionTraits` \cgalParamEnd
  * \cgalNamedParamsEnd
 
  */
@@ -386,7 +386,7 @@ OutputIterator self_intersections(const FaceRange& face_range,
  *    \cgalParamBegin{vertex_point_map} the property map with the points associated to the vertices of `tmesh`.
  *   If this parameter is omitted, an internal property map for
  *   `CGAL::vertex_point_t` should be available in `TriangleMesh`\cgalParamEnd
- *    \cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `SelfIntersectionTraits` \cgalParamEnd
+ *    \cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `PMPSelfIntersectionTraits` \cgalParamEnd
  * \cgalNamedParamsEnd
  *
  * @return true if `tmesh` self-intersects

Polygon_mesh_processing/test/Polygon_mesh_processing/test_pmp_distance.cmd

@@ -0,0 +1 @@
+../data/elephant.off ../data/blobby_3cc.off

Polygon_mesh_processing/test/Polygon_mesh_processing/test_pmp_distance.cpp

@@ -6,49 +6,259 @@
 
 #include <CGAL/boost/graph/property_maps.h>
 
-#include <CGAL/Polygon_mesh_processing/distance.h>
 
 #include <fstream>
 #include <ostream>
 
-// typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
-typedef CGAL::Simple_cartesian<double> K;
+struct Custom_traits_Hausdorff
+{
+// New requirements {
+  struct FT
+  {
+    FT(){}
+    FT(double){}
+    FT operator/(FT)const{return FT();}
+    FT operator-(const FT)const{return FT();}
+    FT operator+(FT)const{return FT();}
+    FT operator*(FT)const{return FT();}
+    bool operator<=(FT)const{return false;}
+    bool operator>=(FT)const{return false;}
+    bool operator!=(FT)const{return false;}
+    bool operator<(FT)const{return false;}
+    bool operator>(FT)const{return false;}
+    FT& operator+=(FT){ return *this; }
+  };
+
+  struct Point_3
+  {
+    Point_3(){}
+    Point_3(FT, FT, FT){}
+    FT operator[](int)const{return FT();}
+    bool operator<(Point_3)const{return false;}
+    double x()const{return 0;}
+    double y()const{return 0;}
+    double z()const{return 0;}
+  };
+
+  struct Vector_3{};
+
+  struct Triangle_3
+  {
+    Triangle_3(){}
+    Triangle_3(const Point_3&,  const Point_3&, const Point_3&){}
+    Point_3 operator[](int)const{return Point_3();}
+    CGAL::Bbox_3 bbox(){return CGAL::Bbox_3();}
+  };
+
+  struct Compute_squared_area_3
+  {
+    typedef FT result_type;
+    FT operator()(Point_3,Point_3,Point_3)const{return FT();}
+    FT operator()(Triangle_3)const{return FT();}
+  };
+
+  struct Construct_translated_point_3
+  {
+    Point_3 operator() (const Point_3 &, const Vector_3 &){return Point_3();}
+  };
+
+  struct Construct_vector_3{
+    Vector_3 	operator() (const Point_3 &, const Point_3 &){return Vector_3();}
+  };
+
+  struct Construct_scaled_vector_3
+  {
+    Vector_3 	operator() (const Vector_3 &, const FT &)
+    {return Vector_3();}
+
+  };
+
+  struct Construct_circumcenter_3
+  {
+    Point_3 operator()(const Point_3&, const Point_3&, const Point_3&){return Point_3();}
+  };
+
+  Compute_squared_area_3 compute_squared_area_3_object(){return Compute_squared_area_3();}
+// } end of new requirements
+
+// requirements from AABBGeomTraits {
+  struct Sphere_3
+  {};
+
+  struct Equal_3
+  {
+    bool operator()(Point_3,Point_3) const {return false;}
+  };
+
+  struct Do_intersect_3
+  {
+    CGAL::Comparison_result operator()(const Sphere_3& , const CGAL::Bbox_3&){ return CGAL::ZERO;}
+  };
+
+  struct Intersect_3
+  {
+    struct result_type{};
+  };
+
+  struct Construct_sphere_3
+  {
+    Sphere_3 operator()(const Point_3& , FT ){return Sphere_3();}
+  };
+
+  struct Construct_projected_point_3
+  {
+    const Point_3 operator()(Triangle_3,  Point_3)const{return Point_3();}
+  };
+
+  struct Compare_distance_3
+  {
+    CGAL::Comparison_result operator()(Point_3, Point_3, Point_3)
+    {return CGAL::ZERO;}
+  };
+
+  struct Has_on_bounded_side_3
+  {
+    //documented as Comparision_result
+    CGAL::Comparison_result operator()(const Point_3&, const Point_3&, const Point_3&)
+    {return CGAL::ZERO;}
+    bool operator()(const Sphere_3&, const Point_3&)
+    {return false;}
+  };
+
+  struct Compute_squared_radius_3
+  {
+    FT operator()(const Sphere_3&){return FT();}
+  };
+
+  struct Compute_squared_distance_3
+  {
+    FT operator()(const Point_3& , const Point_3& ){return FT();}
+  };
+
+  Compare_distance_3 compare_distance_3_object(){return Compare_distance_3();}
+  Construct_sphere_3 construct_sphere_3_object(){return Construct_sphere_3();}
+  Construct_projected_point_3 construct_projected_point_3_object(){return Construct_projected_point_3();}
+  Compute_squared_distance_3 compute_squared_distance_3_object(){return Compute_squared_distance_3();}
+  Do_intersect_3 do_intersect_3_object(){return Do_intersect_3();}
+  Equal_3 equal_3_object(){return Equal_3();}
+// } end of requirments from AABBGeomTraits
+
+
+// requirements from SearchGeomTraits_3 {
+  struct Iso_cuboid_3{};
+
+  struct Construct_iso_cuboid_3{};
+
+  struct Construct_min_vertex_3{};
+
+  struct Construct_max_vertex_3{};
+
+  struct Construct_center_3{};
+
+
+  typedef const FT* Cartesian_const_iterator_3;
+
+
+  struct Construct_cartesian_const_iterator_3{
+    Construct_cartesian_const_iterator_3(){}
+    Construct_cartesian_const_iterator_3(const Point_3&){}
+    const FT* operator()(const Point_3&) const
+    { return 0; }
+
+    const FT* operator()(const Point_3&, int)  const
+    { return 0; }
+    typedef const FT* result_type;
+  };
+// } end of requirements from SearchGeomTraits_3
+
+// requirements from SpatialSortingTraits_3 {
+  struct Less_x_3{
+    bool operator()(Point_3, Point_3){return false;}
+  };
+  struct Less_y_3{
+    bool operator()(Point_3, Point_3){return false;}
+  };
+  struct Less_z_3{
+    bool operator()(Point_3, Point_3){return false;}
+  };
+  Less_x_3 less_x_3_object()const{return Less_x_3();}
+  Less_y_3 less_y_3_object()const{return Less_y_3();}
+  Less_z_3 less_z_3_object()const{return Less_z_3();}
+// } end of requirements from SpatialSortingTraits_3
+};
+
+namespace CGAL{
+template<>struct Kernel_traits<Custom_traits_Hausdorff::Point_3>
+{
+  typedef Custom_traits_Hausdorff Kernel;
+};
+
+template<>
+struct Algebraic_structure_traits< Custom_traits_Hausdorff::FT >
+: public  Algebraic_structure_traits_base< Custom_traits_Hausdorff::FT, Field_tag >
+{
+};
+
+
+template<>
+struct Real_embeddable_traits< Custom_traits_Hausdorff::FT >
+  : public INTERN_RET::Real_embeddable_traits_base< Custom_traits_Hausdorff::FT , CGAL::Tag_true>
+{
+  class To_double
+    : public std::unary_function< Custom_traits_Hausdorff::FT, double > {
+    public:
+      double operator()( const Custom_traits_Hausdorff::FT&  ) const { return 0; }
+  };
+};
+
+}//end CGAL
+
+
+void exact(Custom_traits_Hausdorff::FT){}
+
+//}missing
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+
+#include <CGAL/Polygon_mesh_processing/distance.h>
+
 typedef CGAL::Surface_mesh<K::Point_3> Mesh;
 
-bool general_tests(const Mesh& m1,
-                   const Mesh& m2 )
+template <class GeomTraits, class TriangleMesh>
+void general_tests(const TriangleMesh& m1,
+                   const TriangleMesh& m2 )
 {
 
-    std::cout << "Symetric distance between meshes (sequential) "
+  if (m1.number_of_vertices()==0 || m2.number_of_vertices()==0) return;
+
+  std::vector<typename GeomTraits::Point_3> points;
+    points.push_back(typename GeomTraits::Point_3(0,0,0));
+    points.push_back(typename GeomTraits::Point_3(0,1,0));
+    points.push_back(typename GeomTraits::Point_3(1,0,0));
+    points.push_back(typename GeomTraits::Point_3(0,0,1));
+    points.push_back(typename GeomTraits::Point_3(0,2,0));
+
+  std::cout << "Symmetric distance between meshes (sequential) "
               << CGAL::Polygon_mesh_processing::approximate_symmetric_Hausdorff_distance<CGAL::Sequential_tag>(m1,m2,4000)
               << "\n";
 
-    std::vector<K::Point_3> points;
+  std::cout << "Max distance to point set "
+            << CGAL::Polygon_mesh_processing::approximate_max_distance_to_point_set(m1,points,4000)
+            << "\n";
 
-    points.push_back(K::Point_3(0,0,0));
-    points.push_back(K::Point_3(0,1,0));
-    points.push_back(K::Point_3(1,0,0));
-    points.push_back(K::Point_3(0,0,1));
-    points.push_back(K::Point_3(0,2,0));
-
-    Mesh m;
-    CGAL::make_tetrahedron(points[0],
-                           points[1],
-                           points[2],
-                           points[3],
-                           m);
-
-    std::cout << "Max distance to point set "
-              << CGAL::Polygon_mesh_processing::max_distance_to_point_set(m,points,4000)
-              << "\n";
-
-    std::cout << "Max distance to triangle mesh (sequential) "
-              << CGAL::Polygon_mesh_processing::max_distance_to_triangle_mesh<CGAL::Sequential_tag>(points,m)
-              << "\n";
+  std::cout << "Max distance to triangle mesh (sequential) "
+            << CGAL::Polygon_mesh_processing::max_distance_to_triangle_mesh<CGAL::Sequential_tag>(points,m1)
+            << "\n";
 
 }
 
-int main(int argc, char** argv)
+void test_concept()
+{
+  typedef Custom_traits_Hausdorff CK;
+  CGAL::Surface_mesh<CK::Point_3> m1, m2;
+  general_tests<CK>(m1,m2);
+}
+
+int main(int, char** argv)
 {
   Mesh m1,m2;
 
@@ -64,9 +274,9 @@ int main(int argc, char** argv)
 
   CGAL::Timer time;
   time.start();
-  std::cout << "Distance between meshes (parallel) "
-            << CGAL::Polygon_mesh_processing::approximate_Hausdorff_distance<CGAL::Parallel_tag>(m1,m2,40000)
-            << "\n";
+   std::cout << "Distance between meshes (parallel) "
+             << CGAL::Polygon_mesh_processing::approximate_Hausdorff_distance<CGAL::Parallel_tag>(m1,m2,40000)
+             << "\n";
   time.stop();
   std::cout << "done in " << time.time() << "s.\n";
 
@@ -78,7 +288,11 @@ int main(int argc, char** argv)
   time.stop();
   std::cout << "done in " << time.time() << "s.\n";
 
-  general_tests(m1,m2);
+  general_tests<K>(m1,m2);
+
+  test_concept();
+
+  return 0;
 }
 
 

Spatial_searching/include/CGAL/Orthogonal_k_neighbor_search.h

@@ -145,7 +145,7 @@ private:
       FT diff2 = val - node->upper_low_value();
       if ( (diff1 + diff2 >= FT(0.0)) ) 
       {
-          new_off = 2*val < node->upper_low_value()+node->upper_high_value() ?
+          new_off = node->upper_low_value()+node->upper_high_value() > val*2?
                     val - node->upper_high_value():
                     val - node->upper_low_value();
           bestChild = node->lower();
@@ -153,7 +153,7 @@ private:
       }
       else // compute new distance
       {
-          new_off = 2*val < node->lower_low_value()+node->lower_high_value() ?
+          new_off = node->lower_low_value()+node->lower_high_value() > val*2 ?
                     val - node->lower_high_value():
                     val - node->lower_low_value();
           bestChild = node->upper();