diff --git a/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/Disk_samplers.h b/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/Disk_samplers.h index 77d4402342a..27212b420e9 100644 --- a/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/Disk_samplers.h +++ b/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 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::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(i) / number_of_points); - double angle = atan(R / length_of_normal); - double weight = exp(-0.5 * (std::pow(angle / angle_st_dev, 2))); + double weight = exp(-0.5 * (std::pow(R / CGAL_ANGLE_ST_DEV_DIVIDER, 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 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::epsilon)(); - } - const int number_of_points_sqrt = static_cast(std::sqrt( static_cast(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(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 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::epsilon)(); - } - const int number_of_points_sqrt = static_cast(std::sqrt( static_cast(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); } } diff --git a/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/SDF_calculation.h b/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/SDF_calculation.h index 45b792910e0..a427b0ea9a6 100644 --- a/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/SDF_calculation.h +++ b/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 > 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 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(1.0 / CGAL::sqrt( - v1.squared_length()))); - v2 = scale_functor(v2, static_cast(1.0 / CGAL::sqrt( - v2.squared_length()))); + v1 = scale_functor(v1, static_cast(1.0 / CGAL::sqrt(v1.squared_length()))); + v2 = scale_functor(v2, static_cast(1.0 / CGAL::sqrt(v2.squared_length()))); std::vector > ray_distances; ray_distances.reserve(disk_samples.size()); - const typename GeomTraits::FT length_of_normal = - static_cast( 1.0 / tan(cone_angle / 2.0) ); - Vector scaled_normal = scale_functor(normal, length_of_normal); + const FT normal_multiplier( cos(cone_angle / 2.0) ); + const FT disk_multiplier ( sin(cone_angle / 2.0) ); + + 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(sample_it->get<0>())), - scale_functor(v2, static_cast(sample_it->get<1>())) ); + scale_functor(v1, static_cast(disk_multiplier * sample_it->get<0>())), + scale_functor(v2, static_cast(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);