Poisson_implicit_surface_oracle_3.patch

--- Surface_mesher/include/CGAL/Surface_mesher/Implicit_surface_oracle_3.h	2024-06-06 16:15:39.923520749 +0200
+++ Poisson_surface_reconstruction_3/include/CGAL/Surface_mesher/Poisson_implicit_surface_oracle_3.h	2024-06-06 16:15:39.919520762 +0200
@@ -1,5 +1,5 @@
 // Copyright (c) 2003-2007  INRIA Sophia-Antipolis (France).
-// Copyright (c) 2008       GeometryFactory, Sophia Antipolis (France)
+// Copyright (c) 2008,2011  GeometryFactory Sarl (France)
 // All rights reserved.
 //
 // This file is part of CGAL (www.cgal.org).
@@ -11,20 +11,18 @@
 //
 // Author(s)     : Steve OUDOT, Laurent RINEAU
 
-#ifndef CGAL_SURFACE_MESHER_IMPLICIT_SURFACE_ORACLE_3_H
-#define CGAL_SURFACE_MESHER_IMPLICIT_SURFACE_ORACLE_3_H
+#ifndef CGAL_SURFACE_MESHER_POISSON_IMPLICIT_SURFACE_ORACLE_3_H
+#define CGAL_SURFACE_MESHER_POISSON_IMPLICIT_SURFACE_ORACLE_3_H
 
-#include <CGAL/license/Surface_mesher.h>
+#include <CGAL/license/Poisson_surface_reconstruction_3.h>
 
-#include <CGAL/disable_warnings.h>
 
 #include <CGAL/Surface_mesher/Null_oracle_visitor.h>
 #include <CGAL/point_generators_3.h>
 #include <CGAL/Surface_mesher/Sphere_oracle_3.h>
+#include <CGAL/Surface_mesher/Implicit_surface_oracle_3.h>
 #include <CGAL/Real_embeddable_traits.h>
-#include <CGAL/Surface_mesher/Profile_counter.h>
 #include <CGAL/squared_distance_3.h>
-#include <CGAL/Origin.h>
 
 #include <boost/mpl/has_xxx.hpp>
 
@@ -46,26 +44,8 @@
 // compute the value of the potential in any point of space
 namespace CGAL {
 
-#ifdef  BOOST_MPL_CFG_NO_HAS_XXX
-#  error "BOOST_MPL_HAS_XXX_TRAIT_DEF is needed!"
-#else
-   BOOST_MPL_HAS_XXX_TRAIT_NAMED_DEF(has_set_on_surface, Set_on_surface_tag, true)
-#endif
-
   namespace Surface_mesher {
 
-  namespace {
-
-  template <typename T>
-  struct Return_min : CGAL::cpp98::binary_function<T, T, T>
-  {
-    T operator()(const T& a, const T& b) const
-    {
-      return (std::min)(a, b);
-    }
-  };
-
-  } // end anonymous namespace
 
   template <
     class GT,
@@ -78,10 +58,10 @@
                                             typename GT::Point_3>,
     class Visitor = Null_oracle_visitor
     >
-  class Implicit_surface_oracle_3
+  class Poisson_implicit_surface_oracle_3
   {
     // private types
-    typedef Implicit_surface_oracle_3<GT,
+    typedef Poisson_implicit_surface_oracle_3<GT,
                                       Surface,
                                       Transform_functor_,
                                       Surface_identifiers_generator_,
@@ -124,17 +104,17 @@
   public:
 
     // Constructors
-//     Implicit_surface_oracle_3 (Visitor visitor_ = Visitor()) :
+//     Poisson_implicit_surface_oracle_3 (Visitor visitor_ = Visitor()) :
 //       visitor(visitor_),
 //       transform_functor(),
 //       surface_identifiers_generator()
 //     {
 // #ifdef CGAL_SURFACE_MESHER_DEBUG_CONSTRUCTORS
-//       std::cerr << "CONS: Implicit_surface_oracle_3\n";
+//       std::cerr << "CONS: Poisson_implicit_surface_oracle_3\n";
 // #endif
 //     }
 
-    Implicit_surface_oracle_3 (Transform_functor transform_functor
+    Poisson_implicit_surface_oracle_3 (Transform_functor transform_functor
                                  = Transform_functor(),
                                Surface_identifiers_generator
                                  surface_identifiers_generator
@@ -145,7 +125,7 @@
       surface_identifiers_generator(surface_identifiers_generator)
     {
 #ifdef CGAL_SURFACE_MESHER_DEBUG_CONSTRUCTORS
-      std::cerr << "CONS: Implicit_surface_oracle_3\n";
+      std::cerr << "CONS: Poisson_implicit_surface_oracle_3\n";
 #endif
     }
 
@@ -203,7 +183,7 @@
 
 #ifdef CGAL_SURFACE_MESHER_DEBUG_IMPLICIT_ORACLE
         std::cerr <<
-          boost::format("** Implicit_surface_oracle_3::operator()( (%1%), (%2%) )\n")
+          boost::format("** Poisson_implicit_surface_oracle_3::operator()( (%1%), (%2%) )\n")
           % a % b;
 #endif
         // if both extremities are on the same side of the surface, return
@@ -300,13 +280,19 @@
         typename GT::Construct_midpoint_3 midpoint =
           GT().construct_midpoint_3_object();
 
+        typedef typename Surface::Function::Cell_handle Cell_handle;
+
         // Cannot be const: those values are modified below.
-        result_type value_at_p1 = surf_equation(surface, p1);
-        result_type value_at_p2 = surf_equation(surface, p2);
+        FT value_at_p1, value_at_p2;
+        Cell_handle c1, c2;
+        bool c1_is_inf, c2_is_inf;
+
+        boost::tie(value_at_p1, c1, c1_is_inf) = surface.function().special_func(p1);
+        boost::tie(value_at_p2, c2, c2_is_inf) = surface.function().special_func(p2);
 
         // If both extremities are in the same volume component, returns
         // no intersection.
-        if(value_at_p1 == value_at_p2)
+        if(transform_functor(value_at_p1) == transform_functor(value_at_p2))
           return Object();
 
 #ifdef CGAL_SURFACE_MESHER_PROFILE
@@ -314,6 +300,15 @@
 #endif
         while(true)
         {
+          if(c1 == c2) {
+            if(c1_is_inf) {
+              return Object();
+            } else {
+              return make_object(Point(ORIGIN+ ((value_at_p2 * (p1 - ORIGIN)) - (value_at_p1 * (p2 - ORIGIN))) / (value_at_p2 - value_at_p1)));
+              }
+          }
+
+
 #ifdef CGAL_SURFACE_MESHER_PROFILE
           ++steps;
 #endif
@@ -323,7 +318,10 @@
 #endif
           // mid cannot be const, because it is modified below.
           Point mid = midpoint(p1, p2);
-          const result_type value_at_mid = surf_equation(surface, mid);
+          Cell_handle c_at_mid;
+          FT value_at_mid;
+          bool c_is_inf;
+          boost::tie(value_at_mid, c_at_mid, c_is_inf) = surface.function().special_func(mid);
 
           if ( squared_distance(p1, p2) < squared_distance_bound )
           // If the two points are close, then we must decide
@@ -334,8 +332,8 @@
             // the following function conditionnally call
             // mid.set_on_surface(...) if mid has such a function.
             set_on_surface(mid,
-                           surface_identifiers_generator(value_at_p1,
-                                                         value_at_p2));
+                           surface_identifiers_generator(transform_functor(value_at_p1),
+                                                         transform_functor(value_at_p2)));
 
             visitor.new_point(mid);
             CGAL_SURFACE_MESHER_HISTOGRAM_PROFILER("Implificit_surface_oracle::Intersect_3::operator(Segment_3) bissection steps", steps)
@@ -343,15 +341,19 @@
           }
 
           // Else we must go on
-          if ( value_at_p1 != value_at_mid )
+          if ( transform_functor(value_at_p1) != transform_functor(value_at_mid) )
           {
             p2 = mid;
             value_at_p2 = value_at_mid;
+            c2 = c_at_mid;
+            c2_is_inf = c_is_inf;
           }
           else
           {
             p1 = mid;
             value_at_p1 = value_at_mid;
+            c1 = c_at_mid;
+            c1_is_inf = c_is_inf;
           }
         }
       } // end intersect_clipped_segment
@@ -379,10 +381,11 @@
         const FT radius = CGAL::sqrt(squared_radius);
         typename Self::Intersect_3 intersect = oracle.intersect_3_object();
 
+        Random rng(0);
         typename CGAL::Random_points_on_sphere_3<Point,
-          Point_creator> random_point_on_sphere(CGAL::to_double(radius));
+          Point_creator> random_point_on_sphere(CGAL::to_double(radius), rng);
         typename CGAL::Random_points_in_sphere_3<Point,
-          Point_creator> random_point_in_sphere(CGAL::to_double(radius));
+          Point_creator> random_point_in_sphere(CGAL::to_double(radius), rng);
         typename GT::Construct_segment_3 segment_3 =
           GT().construct_segment_3_object();
         typename GT::Construct_vector_3 vector_3 =
@@ -445,12 +448,16 @@
       return transform_functor(surface(p)) != 0;
     }
 
-  };  // end Implicit_surface_oracle_3
+  };  // end Poisson_implicit_surface_oracle_3
+
+template <typename FT>
+FT approximate_sqrt(const FT x) {
+  return FT (CGAL_NTS sqrt(CGAL_NTS to_double(x)));
+}
 
   }  // namespace Surface_mesher
 
 } // namespace CGAL
 
-#include <CGAL/enable_warnings.h>
 
-#endif  // CGAL_SURFACE_MESHER_IMPLICIT_SURFACE_ORACLE_3_H
+#endif  // CGAL_SURFACE_MESHER_POISSON_IMPLICIT_SURFACE_ORACLE_3_H