local feature size for point set processing

Point_set_processing_3/doc/Point_set_processing_3/PackageDescription.txt

@@ -78,6 +78,7 @@ format.
 - `CGAL::vcm_estimate_normals()`
 - `CGAL::vcm_is_on_feature_edge()`
 - `CGAL::structure_point_set()`
+- `CGAL::estimate_local_feature_size()`
 
 \cgalCRPSection{I/O (All Formats)}
 

Point_set_processing_3/doc/Point_set_processing_3/Point_set_processing_3.txt

@@ -926,6 +926,18 @@ The following example reads a point set from a file, estimates the
 points that are on sharp edges:
 \cgalExample{Point_set_processing_3/edges_example.cpp}
 
+\section Point_set_processing_3LFSEstimation Local Feature Size Estimation
+
+The function `estimate_local_feature_size()` can estimate the local feature size for each point by giving a raw point set. If the point has no normals, the function can also estimate normals using jet-fitting.
+
+The function will estimate the curvature and shape diameter. Then, the local feature size will be the minimum of the curvature radius and half of the shape diameter.
+
+\subsection Point_set_processing_3Example_LFS Example
+
+The following example reads a point set from a file, estimates the
+local feature size for each point:
+\cgalExample{Point_set_processing_3/lfs_example_tuple.cpp}
+
 
 \section Point_set_processing_3Structuring Structuring
 

Point_set_processing_3/doc/Point_set_processing_3/examples.txt

@@ -21,6 +21,8 @@
 \example Point_set_processing_3/bilateral_smooth_point_set_example.cpp
 \example Point_set_processing_3/edge_aware_upsample_point_set_example.cpp
 \example Point_set_processing_3/edges_example.cpp
+\example Point_set_processing_3/lfs_example_tuple.cpp
+\example Point_set_processing_3/lfs_example_pointset.cpp
 \example Point_set_processing_3/structuring_example.cpp
 \example Point_set_processing_3/callback_example.cpp
 */

Point_set_processing_3/examples/Point_set_processing_3/CMakeLists.txt

@@ -72,6 +72,8 @@ if(TARGET CGAL::Eigen3_support)
   # Executables that require Eigen
   foreach(
     target
+    lfs_example_tuple
+    lfs_example_pointset
     jet_smoothing_example
     normal_estimation
     clustering_example

Point_set_processing_3/examples/Point_set_processing_3/lfs_example_pointset.cpp

@@ -0,0 +1,61 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/estimate_lfs.h>
+#include <CGAL/IO/read_points.h>
+#include <CGAL/Point_set_3.h>
+#include <CGAL/Point_set_3/IO.h>
+
+#include <vector>
+
+
+
+// types
+typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
+typedef Kernel::FT FT;
+typedef Kernel::Point_3 Point;
+typedef Kernel::Vector_3 Vector;
+
+typedef CGAL::Point_set_3<Point, Vector> Point_set;
+typedef Point_set::Property_map<Point> Point_map;
+typedef Point_set::Property_map<FT> FT_map;
+
+// types for K nearest neighbors search structure
+typedef CGAL::Point_set_processing_3::internal::Neighbor_query<Kernel, Point_set&, Point_map> Neighbor_query;
+
+// Concurrency
+typedef CGAL::Parallel_if_available_tag Concurrency_tag;
+
+int main(void)
+{
+  // read xyz
+  const std::string filename = CGAL::data_file_path("points_3/kitten.xyz");
+
+  Point_set point_set;
+  point_set.add_normal_map();
+
+  FT_map lfs_map;
+  boost::tie (lfs_map, boost::tuples::ignore) = point_set.add_property_map<FT> ("LFS", 0.);
+
+  if (!CGAL::IO::read_point_set(filename, point_set))
+  {
+    std::cerr << "Error: cannot read file " << filename<< std::endl;
+    return EXIT_FAILURE;
+  }
+
+
+  unsigned int jet_k = 24;
+  std::size_t N_rays = 60;
+  FT apex_angle = 30;
+  CGAL::estimate_local_feature_size<Concurrency_tag>(point_set, jet_k, N_rays, apex_angle, lfs_map,
+    CGAL::parameters::point_map(point_set.point_map())
+                                                .normal_map(point_set.normal_map()));
+  // print
+  for (Point_set::iterator it = point_set.begin(); it != point_set.end(); ++ it)
+   {
+      // Point p = point_set.point(*it);
+      // Vector n = point_set.normal(*it);
+      FT lfs = lfs_map[*it];
+      std::cerr << lfs << std::endl;
+   }
+
+  return EXIT_SUCCESS;
+}

Point_set_processing_3/examples/Point_set_processing_3/lfs_example_tuple.cpp

@@ -0,0 +1,52 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/estimate_lfs.h>
+#include <CGAL/IO/read_points.h>
+
+#include <vector>
+#include <utility>
+#include <tuple>
+
+// types
+typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
+typedef Kernel::FT FT;
+typedef Kernel::Point_3 Point;
+typedef Kernel::Vector_3 Vector;
+
+typedef std::tuple<Point, Vector, FT> Point_with_normal_and_lfs;
+
+// Concurrency
+typedef CGAL::Parallel_if_available_tag Concurrency_tag;
+
+int main(void)
+{
+  // read xyz
+  const std::string filename = CGAL::data_file_path("points_3/kitten.xyz");
+
+  std::vector<Point_with_normal_and_lfs> points;
+  if(!CGAL::IO::read_points(filename,
+                            std::back_inserter(points),
+                            CGAL::parameters::point_map(CGAL::Nth_of_tuple_property_map<0, Point_with_normal_and_lfs>())
+                            .normal_map(CGAL::Nth_of_tuple_property_map<1, Point_with_normal_and_lfs>())))
+  {
+    std::cerr << "Error: cannot read file " << filename<< std::endl;
+    return EXIT_FAILURE;
+  }
+
+  unsigned int jet_k = 24;
+  std::size_t N_rays = 60;
+  FT apex_angle = 30;
+  auto lfs_map = CGAL::Nth_of_tuple_property_map<2, Point_with_normal_and_lfs>();
+  CGAL::estimate_local_feature_size<Concurrency_tag>(points,
+                                                     jet_k,
+                                                     N_rays,
+                                                     apex_angle,
+                                                     lfs_map,
+                                                     CGAL::parameters::point_map(CGAL::Nth_of_tuple_property_map<0, Point_with_normal_and_lfs>())
+                                                                      .normal_map(CGAL::Nth_of_tuple_property_map<1, Point_with_normal_and_lfs>()));
+
+  for (const auto &pts : points){
+    std::cerr << std::get<2>(pts) << std::endl;
+  }
+
+  return EXIT_SUCCESS;
+}