Changes to guarantee there will be no border case problems with cone angle in [0, PI].

Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/Disk_samplers.h

@@ -58,32 +58,21 @@ public:
   /**
    * Samples points from unit-disk.
    * @param number_of_points number of points to be picked
-   * @param cone_angle opening angle of cone (might be necessary for weighting)
    * @param[out] out_it sampled points from disk, each point is tuple of:
    *   - coordinate-x
    *   - coordinate-y
-   *   - weight (proportional to angle between cone-normal)
+   *   - weight (proportional to distance from disk center)
    */
   template<class OutputIterator>
   void operator()(int number_of_points,
-                  double cone_angle,
                   OutputIterator out_it) const {
-    if(cone_angle <= 0.0) {
-      CGAL_warning(false
-                   && "Warning: cone angle smaller than or equal to zero will be assigned to epsilon!");
-      cone_angle = (std::numeric_limits<double>::epsilon)();
-    }
-
     const double golden_ratio = 3.0 - std::sqrt(5.0);
 
     if(uniform) {
-      const double length_of_normal = 1.0 / tan(cone_angle / 2.0);
-      const double angle_st_dev = cone_angle / CGAL_ANGLE_ST_DEV_DIVIDER;
       for(int i = 0; i < number_of_points; ++i) {
         double Q = i * golden_ratio * CGAL_PI;
         double R = std::sqrt(static_cast<double>(i) / number_of_points);
-        double angle = atan(R / length_of_normal);
+        double weight =  exp(-0.5 * (std::pow(R / CGAL_ANGLE_ST_DEV_DIVIDER, 2)));
-        double weight =  exp(-0.5 * (std::pow(angle / angle_st_dev, 2)));
         *out_it++ = Tuple(R * cos(Q), R * sin(Q), weight);
       }
     } else {
@@ -137,29 +126,19 @@ public:
   /**
    * Samples points from unit-disk.
    * @param number_of_points number of points to be picked
-   * @param cone_angle opening angle of cone (might be necessary for weighting)
    * @param[out] out_it sampled points from disk, each point is tuple of:
    *   - coordinate-x
    *   - coordinate-y
-   *   - weight (proportional to angle between cone-normal)
+   *   - weight (proportional to distance from disk center)
    *
    * Note: returned samples size = \f$ \lfloor \sqrt {number\_of\_points} \rfloor ^ 2 \f$
    */
   template<class OutputIterator>
   void operator()(int number_of_points,
-                  double cone_angle,
                   OutputIterator out_it) const {
-    if(cone_angle <= 0.0) {
-      CGAL_warning(false
-                   && "Warning: cone angle smaller than or equal to zero will be assigned to epsilon!");
-      cone_angle = (std::numeric_limits<double>::epsilon)();
-    }
-
     const int number_of_points_sqrt = static_cast<int>(std::sqrt(
                                         static_cast<double>(number_of_points)));
-    const double length_of_normal = 1.0 / tan(cone_angle / 2.0);
     // use cone_angle / 3 as one standard deviation while weighting.
-    const double angle_st_dev = cone_angle / CGAL_ANGLE_ST_DEV_DIVIDER;
 
     for(int i = 1; i <= number_of_points_sqrt; ++i)
       for(int j = 1; j <= number_of_points_sqrt; ++j) {
@@ -167,8 +146,8 @@ public:
         double w2 = static_cast<double>(j) / number_of_points_sqrt;
         double R = w1;
         double Q = 2 * w2 * CGAL_PI;
-        double angle = atan(R / length_of_normal);
+
-        double weight = exp(-0.5 * (pow(angle / angle_st_dev, 2)));
+        double weight = exp(-0.5 * (pow(R / CGAL_ANGLE_ST_DEV_DIVIDER, 2)));
         *out_it++ = Tuple(R * cos(Q), R * sin(Q), weight);
       }
   }
@@ -219,21 +198,11 @@ public:
    */
   template<class OutputIterator>
   void operator()(int number_of_points,
-                  double cone_angle,
                   OutputIterator out_it) const {
-    if(cone_angle <= 0.0) {
-      CGAL_warning(false
-                   && "Warning: cone angle smaller than or equal to zero will be assigned to epsilon!");
-      cone_angle = (std::numeric_limits<double>::epsilon)();
-    }
-
     const int number_of_points_sqrt = static_cast<int>(std::sqrt(
                                         static_cast<double>(number_of_points)));
-    const double length_of_normal = 1.0 / tan(cone_angle / 2.0);
     const double fraction = (number_of_points_sqrt == 1) ? 0.0
                             : 2.0 / (number_of_points_sqrt -1);
-    // use cone_angle / 3 as one standard deviation while weighting.
-    const double angle_st_dev = cone_angle / CGAL_ANGLE_ST_DEV_DIVIDER;
 
     for(int i = 0; i < number_of_points_sqrt; ++i)
       for(int j = 0; j < number_of_points_sqrt; ++j) {
@@ -261,8 +230,8 @@ public:
           }
         }
         Q *= (CGAL_PI / 4.0);
-        double angle = atan(R / length_of_normal);
+
-        double weight = exp(-0.5 * (pow(angle / angle_st_dev, 2)));
+        double weight = exp(-0.5 * (pow(R / CGAL_ANGLE_ST_DEV_DIVIDER, 2)));
         *out_it++ = Tuple(R * cos(Q), R * sin(Q), weight);
       }
   }

Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/SDF_calculation.h

@@ -60,6 +60,7 @@ private:
   typedef typename GeomTraits::Ray_3      Ray;
   typedef typename GeomTraits::Plane_3    Plane;
   typedef typename GeomTraits::Segment_3  Segment;
+  typedef typename GeomTraits::FT         FT;
 
   typedef typename Polyhedron::Facet  Facet;
 
@@ -80,6 +81,7 @@ private:
 
   typedef Vogel_disk_sampling<boost::tuple<double, double, double> >
   Default_sampler;
+
 // member variables
 private:
   GeomTraits traits;
@@ -168,7 +170,7 @@ public:
     FacetValueMap sdf_values,
     DiskSampling disk_sampler) const {
     Disk_samples_list disk_samples;
-    disk_sampler(number_of_rays, cone_angle, std::back_inserter(disk_samples));
+    disk_sampler(number_of_rays, std::back_inserter(disk_samples));
 
     for( ; facet_begin != facet_end; ++facet_begin) {
       boost::optional<double> sdf_value = calculate_sdf_value_of_facet(facet_begin,
@@ -232,7 +234,7 @@ public:
     bool accept_if_acute,
     DiskSampling disk_sampler) const {
     Disk_samples_list disk_samples;
-    disk_sampler(number_of_rays, cone_angle, std::back_inserter(disk_samples));
+    disk_sampler(number_of_rays, std::back_inserter(disk_samples));
 
     return calculate_sdf_value_of_point(center, normal, skip, visitor, cone_angle,
                                         accept_if_acute, disk_samples);
@@ -250,19 +252,23 @@ public:
     double cone_angle,
     bool accept_if_acute,
     const Disk_samples_list& disk_samples) const {
+    if(cone_angle < 0.0 || cone_angle > CGAL_PI) {
+      CGAL_warning(false && "Cone angle is clamped between [0, CGAL_PI].");
+      cone_angle = (std::min)(CGAL_PI, (std::max)(0.0, cone_angle));
+    }
+
     Plane plane(center, normal);
     Vector v1 = plane.base1(), v2 = plane.base2();
-    v1 = scale_functor(v1, static_cast<typename GeomTraits::FT>(1.0 / CGAL::sqrt(
+    v1 = scale_functor(v1, static_cast<FT>(1.0 / CGAL::sqrt(v1.squared_length())));
-                         v1.squared_length())));
+    v2 = scale_functor(v2, static_cast<FT>(1.0 / CGAL::sqrt(v2.squared_length())));
-    v2 = scale_functor(v2, static_cast<typename GeomTraits::FT>(1.0 / CGAL::sqrt(
-                         v2.squared_length())));
 
     std::vector<std::pair<double, double> > ray_distances;
     ray_distances.reserve(disk_samples.size());
 
-    const typename GeomTraits::FT length_of_normal =
+    const FT normal_multiplier( cos(cone_angle / 2.0) );
-      static_cast<typename GeomTraits::FT>( 1.0 / tan(cone_angle / 2.0) );
+    const FT disk_multiplier  ( sin(cone_angle / 2.0) );
-    Vector scaled_normal = scale_functor(normal, length_of_normal);
+
+    Vector scaled_normal = scale_functor(normal, normal_multiplier);
 
     for(Disk_samples_list::const_iterator sample_it = disk_samples.begin();
         sample_it != disk_samples.end(); ++sample_it) {
@@ -271,11 +277,17 @@ public:
       Primitive_id closest_id;
 
       Vector disk_vector = sum_functor(
-                             scale_functor(v1, static_cast<typename GeomTraits::FT>(sample_it->get<0>())),
+                             scale_functor(v1, static_cast<FT>(disk_multiplier * sample_it->get<0>())),
-                             scale_functor(v2, static_cast<typename GeomTraits::FT>(sample_it->get<1>())) );
+                             scale_functor(v2, static_cast<FT>(disk_multiplier * sample_it->get<1>())) );
       Vector ray_direction = sum_functor(scaled_normal, disk_vector);
 
       Ray ray(center, ray_direction);
+
+      if(traits.is_degenerate_3_object()(ray)) {
+        CGAL_warning(false &&
+                     "A degenerate ray is constructed. Most probable reason is using CGAL_PI as cone_angle parameter and also picking center of disk as a sample.");
+      }
+
       boost::tie(is_intersected, intersection_is_acute, min_distance, closest_id)
         = cast_and_return_minimum(ray, skip, accept_if_acute);