single vertex computation

implemented single vertex curvature computation
compared against the whole mesh computation on both cases of no radius and with radius on some vertices
still need to add tests, documentation and an example file

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/Curvatures/interpolated_corrected_curvature_measures.h

@@ -79,8 +79,8 @@ struct Principal_curvatures_and_directions {
 
 namespace internal {
 
-  template<typename PolygonMesh, typename GT>
+template<typename PolygonMesh, typename GT>
-  typename GT::FT average_edge_length(const PolygonMesh& pmesh) {
+typename GT::FT average_edge_length(const PolygonMesh& pmesh) {
   const std::size_t n = edges(pmesh).size();
   if (n == 0)
     return 0;
@@ -91,28 +91,28 @@ namespace internal {
 
   avg_edge_length /= n;
   return avg_edge_length;
-  }
+}
 
-  enum Curvature_measure_index {
+enum Curvature_measure_index {
   MU0_AREA_MEASURE, ///< corrected area density
   MU1_MEAN_CURVATURE_MEASURE, ///< corrected mean curvature density
   MU2_GAUSSIAN_CURVATURE_MEASURE ///< corrected gaussian curvature density
-  };
+};
 
-  template<typename GT>
+template<typename GT>
-  struct Vertex_measures {
+struct Vertex_measures {
   typename GT::FT area_measure = 0;
   typename GT::FT mean_curvature_measure = 0;
   typename GT::FT gaussian_curvature_measure = 0;
   std::array<typename GT::FT, 3 * 3> anisotropic_measure = { 0, 0, 0,
                                                              0, 0, 0,
                                                              0, 0, 0 };
-  };
+};
 
-  template<typename GT>
+template<typename GT>
-  typename GT::FT interpolated_corrected_area_measure_face(const std::vector<typename GT::Vector_3>& u,
+typename GT::FT interpolated_corrected_area_measure_face(const std::vector<typename GT::Vector_3>& u,
   const std::vector<typename GT::Vector_3>& x)
-  {
+{
   const std::size_t n = x.size();
   CGAL_precondition(u.size() == n);
   CGAL_precondition(n >= 3);
@@ -162,12 +162,12 @@ namespace internal {
     }
     return mu0;
   }
-  }
+}
 
-  template<typename GT>
+template<typename GT>
-  typename GT::FT interpolated_corrected_mean_curvature_measure_face(const std::vector<typename GT::Vector_3>& u,
+typename GT::FT interpolated_corrected_mean_curvature_measure_face(const std::vector<typename GT::Vector_3>& u,
   const std::vector<typename GT::Vector_3>& x)
-  {
+{
   const std::size_t n = x.size();
   CGAL_precondition(u.size() == n);
   CGAL_precondition(n >= 3);
@@ -228,12 +228,12 @@ namespace internal {
     }
     return mu1;
   }
-  }
+}
 
-  template<typename GT>
+template<typename GT>
-  typename GT::FT interpolated_corrected_gaussian_curvature_measure_face(const std::vector<typename GT::Vector_3>& u,
+typename GT::FT interpolated_corrected_gaussian_curvature_measure_face(const std::vector<typename GT::Vector_3>& u,
   const std::vector<typename GT::Vector_3>& x = {})
-  {
+{
   const std::size_t n = u.size();
   CGAL_precondition(n >= 3);
 
@@ -275,12 +275,12 @@ namespace internal {
     }
     return mu2;
   }
-  }
+}
 
-  template<typename GT>
+template<typename GT>
-  std::array<typename GT::FT, 3 * 3> interpolated_corrected_anisotropic_measure_face(const std::vector<typename GT::Vector_3>& u,
+std::array<typename GT::FT, 3 * 3> interpolated_corrected_anisotropic_measure_face(const std::vector<typename GT::Vector_3>& u,
   const std::vector<typename GT::Vector_3>& x)
-  {
+{
   const std::size_t n = x.size();
   CGAL_precondition(u.size() == n);
   CGAL_precondition(n >= 3);
@@ -386,33 +386,33 @@ namespace internal {
     }
   }
   return muXY;
-  }
+}
 
-  //template<typename GT>
+//template<typename GT>
-  //typename GT::FT triangle_in_ball_ratio(const typename GT::Vector_3 x1,
+//typename GT::FT triangle_in_ball_ratio(const typename GT::Vector_3 x1,
-  //                                       const typename GT::Vector_3 x2,
+//                                       const typename GT::Vector_3 x2,
-  //                                       const typename GT::Vector_3 x3,
+//                                       const typename GT::Vector_3 x3,
-  //                                       const typename GT::FT r,
+//                                       const typename GT::FT r,
-  //                                       const typename GT::Vector_3 c,
+//                                       const typename GT::Vector_3 c,
-  //                                       const std::size_t res = 3)
+//                                       const std::size_t res = 3)
-  //{
+//{
-  //    const typename GT::FT R = r * r;
+//    const typename GT::FT R = r * r;
-  //    const typename GT::FT acc = 1.0 / res;
+//    const typename GT::FT acc = 1.0 / res;
-  //    std::size_t samples_in = 0;
+//    std::size_t samples_in = 0;
-  //    for (GT::FT alpha = acc / 3; alpha < 1; alpha += acc)
+//    for (GT::FT alpha = acc / 3; alpha < 1; alpha += acc)
-  //        for (GT::FT beta = acc / 3; beta < 1 - alpha; beta += acc)
+//        for (GT::FT beta = acc / 3; beta < 1 - alpha; beta += acc)
-  //        {
+//        {
-  //            if ((alpha * x1 + beta * x2 + (1 - alpha - beta) * x3 - c).squared_length() < R)
+//            if ((alpha * x1 + beta * x2 + (1 - alpha - beta) * x3 - c).squared_length() < R)
-  //                samples_in++;
+//                samples_in++;
-  //        }
+//        }
-  //    return samples_in / (typename GT::FT)(res * (res + 1) / 2);
+//    return samples_in / (typename GT::FT)(res * (res + 1) / 2);
-  //}
+//}
 
-  template<typename GT>
+template<typename GT>
-  typename GT::FT face_in_ball_ratio(const std::vector<typename GT::Vector_3>& x,
+typename GT::FT face_in_ball_ratio(const std::vector<typename GT::Vector_3>& x,
   const typename GT::FT r,
   const typename GT::Vector_3 c)
-  {
+{
   const std::size_t n = x.size();
 
   // getting center of points
@@ -437,15 +437,15 @@ namespace internal {
   d_min = sqrt(d_min);
 
   return (r - d_min) / (d_max - d_min);
-  }
+}
 
-  template<typename GT>
+template<typename GT>
-  Principal_curvatures_and_directions<GT> principal_curvatures_and_directions_from_anisotropic_measures(
+Principal_curvatures_and_directions<GT> principal_curvatures_and_directions_from_anisotropic_measures(
   const std::array<typename GT::FT, 3 * 3> anisotropic_measure,
   const typename GT::FT v_mu0,
   const typename GT::Vector_3 u_GT
-  )
+)
-  {
+{
   Eigen::Matrix<typename GT::FT, 3, 3> v_muXY = Eigen::Matrix<typename GT::FT, 3, 3>::Zero();
 
   for (std::size_t ix = 0; ix < 3; ix++)
@@ -476,11 +476,247 @@ namespace internal {
     min_eig_vec,
     max_eig_vec
     );
+}
+
+template<typename GT, typename PolygonMesh, typename VPM, typename VNM>
+typename Vertex_measures<GT> interpolated_corrected_measures_one_vertex_no_radius(
+  const PolygonMesh pmesh,
+  const typename boost::graph_traits<PolygonMesh>::vertex_descriptor v,
+  const bool is_mean_curvature_selected,
+  const bool is_gaussian_curvature_selected,
+  const bool is_principal_curvatures_and_directions_selected,
+  const VPM vpm,
+  const VNM vnm
+)
+{
+  typedef typename boost::graph_traits<PolygonMesh>::face_descriptor Face_descriptor;
+  typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor Vertex_descriptor;
+  typedef typename GT::Point_3 Point_3;
+  typedef typename GT::Vector_3 Vector_3;
+  typedef typename GT::FT FT;
+
+  std::queue<Face_descriptor> bfs_queue;
+  std::unordered_set<Face_descriptor> bfs_visited;
+
+  typename Vertex_measures<GT> vertex_measures;
+
+  std::vector<Vector_3> x;
+  std::vector<Vector_3> u;
+
+  for (Face_descriptor f : faces_around_target(halfedge(v, pmesh), pmesh)) {
+    if (f != boost::graph_traits<PolygonMesh>::null_face())
+    {
+      for (Vertex_descriptor vi : vertices_around_face(halfedge(f, pmesh), pmesh))
+      {
+        Point_3 pi = get(vpm, vi);
+        Vector_3 ui = get(vnm, vi);
+        x.push_back(Vector_3(pi.x(), pi.y(), pi.z()));
+        u.push_back(ui);
       }
 
-  template<typename PolygonMesh, typename NamedParameters = parameters::Default_named_parameters>
+      vertex_measures.area_measure += interpolated_corrected_area_measure_face<GT>(u, x);
-  class Interpolated_corrected_curvatures_computer
+
+      if (is_mean_curvature_selected)
+        vertex_measures.mean_curvature_measure += interpolated_corrected_mean_curvature_measure_face<GT>(u, x);
+
+      if (is_gaussian_curvature_selected)
+        vertex_measures.gaussian_curvature_measure += interpolated_corrected_gaussian_curvature_measure_face<GT>(u, x);
+
+      if (is_principal_curvatures_and_directions_selected)
       {
+        const std::array<FT, 3 * 3> face_anisotropic_measure = interpolated_corrected_anisotropic_measure_face<GT>(u, x);
+        for (std::size_t i = 0; i < 3 * 3; i++)
+          vertex_measures.anisotropic_measure[i] += face_anisotropic_measure[i];
+      }
+    }
+
+    x.clear();
+    u.clear();
+  }
+  return vertex_measures;
+}
+
+
+template<typename GT, typename PolygonMesh, typename VPM, typename VNM>
+typename Vertex_measures<GT> interpolated_corrected_measures_one_vertex(
+  const PolygonMesh pmesh,
+  const typename boost::graph_traits<PolygonMesh>::vertex_descriptor v,
+  const typename GT::FT radius,
+  const bool is_mean_curvature_selected,
+  const bool is_gaussian_curvature_selected,
+  const bool is_principal_curvatures_and_directions_selected,
+  const VPM vpm,
+  const VNM vnm
+)
+{
+  typedef typename boost::graph_traits<PolygonMesh>::face_descriptor Face_descriptor;
+  typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor Vertex_descriptor;
+  typedef typename GT::Point_3 Point_3;
+  typedef typename GT::Vector_3 Vector_3;
+  typedef typename GT::FT FT;
+
+  std::queue<Face_descriptor> bfs_queue;
+  std::unordered_set<Face_descriptor> bfs_visited;
+
+  typename Vertex_measures<GT> vertex_measures;
+
+  typename Point_3 vp = get(vpm, v);
+  typename Vector_3 c = Vector_3(vp.x(), vp.y(), vp.z());
+
+  for (Face_descriptor f : faces_around_target(halfedge(v, pmesh), pmesh)) {
+    if (f != boost::graph_traits<PolygonMesh>::null_face())
+    {
+      bfs_queue.push(f);
+      bfs_visited.insert(f);
+    }
+  }
+  std::vector<Vector_3> x;
+  std::vector<Vector_3> u;
+  while (!bfs_queue.empty()) {
+    Face_descriptor fi = bfs_queue.front();
+    bfs_queue.pop();
+
+    // looping over vertices in face to get point coordinates
+
+    for (Vertex_descriptor vi : vertices_around_face(halfedge(fi, pmesh), pmesh))
+    {
+      Point_3 pi = get(vpm, vi);
+      Vector_3 ui = get(vnm, vi);
+      x.push_back(Vector_3(pi.x(), pi.y(), pi.z()));
+      u.push_back(ui);
+    }
+
+    const FT f_ratio = face_in_ball_ratio<GT>(x, radius, c);
+
+    if (f_ratio != 0.0)
+    {
+      vertex_measures.area_measure += f_ratio * interpolated_corrected_area_measure_face<GT>(u, x);
+
+      if (is_mean_curvature_selected)
+        vertex_measures.mean_curvature_measure += f_ratio * interpolated_corrected_mean_curvature_measure_face<GT>(u, x);
+
+      if (is_gaussian_curvature_selected)
+        vertex_measures.gaussian_curvature_measure += f_ratio * interpolated_corrected_gaussian_curvature_measure_face<GT>(u, x);
+
+      if (is_principal_curvatures_and_directions_selected)
+      {
+        const std::array<FT, 3 * 3> face_anisotropic_measure = interpolated_corrected_anisotropic_measure_face<GT>(u, x);
+        for (std::size_t i = 0; i < 3 * 3; i++)
+          vertex_measures.anisotropic_measure[i] += f_ratio * face_anisotropic_measure[i];
+      }
+
+      for (Face_descriptor fj : faces_around_face(halfedge(fi, pmesh), pmesh))
+      {
+        if (bfs_visited.find(fj) == bfs_visited.end() && fj != boost::graph_traits<PolygonMesh>::null_face())
+        {
+          bfs_queue.push(fj);
+          bfs_visited.insert(fj);
+        }
+      }
+    }
+
+    x.clear();
+    u.clear();
+  }
+  return vertex_measures;
+}
+
+template<typename PolygonMesh,
+  typename NamedParameters = parameters::Default_named_parameters>
+  void interpolated_corrected_curvatures_one_vertex(
+    const PolygonMesh pmesh,
+    const typename boost::graph_traits<PolygonMesh>::vertex_descriptor v,
+    NamedParameters& np = parameters::default_values()
+  )
+{
+  typedef typename GetGeomTraits<PolygonMesh, NamedParameters>::type GT;
+  typedef typename GetVertexPointMap<PolygonMesh, NamedParameters>::const_type Vertex_position_map;
+
+  typedef dynamic_vertex_property_t<GT::Vector_3> Vector_map_tag;
+  typedef typename boost::property_map<PolygonMesh, Vector_map_tag>::const_type Default_vector_map;
+  typedef typename internal_np::Lookup_named_param_def<internal_np::vertex_normal_map_t,
+    NamedParameters,
+    Default_vector_map>::type Vertex_normal_map;
+
+  using parameters::choose_parameter;
+  using parameters::get_parameter;
+  using parameters::is_default_parameter;
+
+  Vertex_position_map vpm = choose_parameter(get_parameter(np, CGAL::vertex_point),
+    get_const_property_map(CGAL::vertex_point, pmesh));
+
+  Vertex_normal_map vnm = choose_parameter(get_parameter(np, internal_np::vertex_normal_map),
+    get(Vector_map_tag(), pmesh));
+
+  if (is_default_parameter<NamedParameters, internal_np::vertex_normal_map_t>::value)
+    compute_vertex_normals(pmesh, vnm, np);
+
+  typename GT::FT radius = choose_parameter(get_parameter(np, internal_np::ball_radius), -1);
+  if (radius == 0)
+    radius = average_edge_length<PolygonMesh, GT>(pmesh) * EXPANDING_RADIUS_EPSILON;
+  else
+    radius = radius;
+
+  typename GT::FT* vertex_mean_curvature = choose_parameter(get_parameter(np, internal_np::vertex_mean_curvature), nullptr);
+  typename GT::FT* vertex_gaussian_curvature = choose_parameter(get_parameter(np, internal_np::vertex_gaussian_curvature), nullptr);
+  typename Principal_curvatures_and_directions<GT>* vertex_principal_curvatures_and_directions = choose_parameter(get_parameter(np, internal_np::vertex_principal_curvatures_and_directions), nullptr);
+
+  const bool is_mean_curvature_selected = (vertex_mean_curvature != nullptr);
+  const bool is_gaussian_curvature_selected = (vertex_gaussian_curvature != nullptr);
+  const bool is_principal_curvatures_and_directions_selected = (vertex_principal_curvatures_and_directions != nullptr);
+
+  std::cout << is_mean_curvature_selected << is_gaussian_curvature_selected << is_principal_curvatures_and_directions_selected << std::endl;
+
+  Vertex_measures<GT> vertex_measures;
+
+  if (radius < 0)
+    vertex_measures = interpolated_corrected_measures_one_vertex_no_radius<GT>(
+      pmesh,
+      v,
+      is_mean_curvature_selected,
+      is_gaussian_curvature_selected,
+      is_principal_curvatures_and_directions_selected,
+      vpm,
+      vnm
+      );
+  else
+    vertex_measures = interpolated_corrected_measures_one_vertex<GT>(
+    pmesh,
+    v,
+    radius,
+    is_mean_curvature_selected,
+    is_gaussian_curvature_selected,
+    is_principal_curvatures_and_directions_selected,
+    vpm,
+    vnm
+    );
+
+
+  if (is_mean_curvature_selected) {
+    *vertex_mean_curvature = vertex_measures.area_measure != 0 ?
+      0.5 * vertex_measures.mean_curvature_measure / vertex_measures.area_measure : 0;
+  }
+
+  if (is_gaussian_curvature_selected) {
+    *vertex_gaussian_curvature = vertex_measures.area_measure != 0 ?
+      vertex_measures.gaussian_curvature_measure / vertex_measures.area_measure : 0;
+  }
+
+  if (is_principal_curvatures_and_directions_selected) {
+    const GT::Vector_3  v_normal = get(vnm, v);
+    const Principal_curvatures_and_directions<GT> principal_curvatures_and_directions = principal_curvatures_and_directions_from_anisotropic_measures<GT>(
+      vertex_measures.anisotropic_measure,
+      vertex_measures.area_measure,
+      v_normal
+      );
+    *vertex_principal_curvatures_and_directions = principal_curvatures_and_directions;
+  }
+}
+
+
+template<typename PolygonMesh, typename NamedParameters = parameters::Default_named_parameters>
+class Interpolated_corrected_curvatures_computer
+{
   typedef typename GetGeomTraits<PolygonMesh, NamedParameters>::type GT;
 
   typedef typename GT::FT FT;
@@ -519,7 +755,7 @@ namespace internal {
   typedef typename boost::property_map<PolygonMesh,
     CGAL::dynamic_face_property_t<std::array<FT, 3 * 3>>>::const_type Face_anisotropic_measure_map;
 
-  private:
+private:
   const PolygonMesh& pmesh;
   Vertex_position_map vpm;
   Vertex_normal_map vnm;
@@ -579,7 +815,7 @@ namespace internal {
       ball_radius = radius;
   }
 
-  public:
+public:
 
   Interpolated_corrected_curvatures_computer(const PolygonMesh& pmesh,
     const NamedParameters& np = parameters::default_values()
@@ -596,7 +832,7 @@ namespace internal {
     }
   }
 
-  private:
+private:
 
   void interpolated_corrected_selected_measures_all_faces()
   {
@@ -747,7 +983,7 @@ namespace internal {
       }
     }
   }
-  };
+};
 
 } // namespace internal
 
@@ -1030,6 +1266,50 @@ template<typename PolygonMesh,
   internal::Interpolated_corrected_curvatures_computer<PolygonMesh, NamedParameters>(pmesh, np);
 }
 
+template<typename GT, typename PolygonMesh,
+  typename NamedParameters = parameters::Default_named_parameters>
+  typename GT::FT interpolated_corrected_mean_curvature_at_vertex(const PolygonMesh& pmesh,
+    typename boost::graph_traits<PolygonMesh>::vertex_descriptor v,
+    const NamedParameters& np = parameters::default_values())
+{
+  // use interpolated_corrected_curvatures_at_vertex to compute mean curvature
+  typename GT::FT* mean_curvature = new typename GT::FT();
+  interpolated_corrected_curvatures_at_vertex(pmesh, v, np.vertex_mean_curvature(mean_curvature));
+  return *mean_curvature;
+}
+
+template<typename GT, typename PolygonMesh,
+  typename NamedParameters = parameters::Default_named_parameters>
+  typename GT::FT interpolated_corrected_gaussian_curvature_at_vertex(const PolygonMesh& pmesh,
+    typename boost::graph_traits<PolygonMesh>::vertex_descriptor v,
+    const NamedParameters& np = parameters::default_values())
+{
+  // use interpolated_corrected_curvatures_at_vertex to compute gaussian curvature
+  typename GT::FT* gc = new typename GT::FT();
+  interpolated_corrected_curvatures_at_vertex(pmesh, v, np.vertex_gaussian_curvature(gc));
+  return *gc;
+}
+
+template<typename GT, typename PolygonMesh,
+  typename NamedParameters = parameters::Default_named_parameters>
+  Principal_curvatures_and_directions<GT> interpolated_corrected_principal_curvatures_and_directions_at_vertex(const PolygonMesh& pmesh,
+    typename boost::graph_traits<PolygonMesh>::vertex_descriptor v,
+    const NamedParameters& np = parameters::default_values())
+{
+  // use interpolated_corrected_curvatures_at_vertex to compute principal curvatures
+  Principal_curvatures_and_directions<GT>* pcd = new Principal_curvatures_and_directions<GT>();
+  interpolated_corrected_curvatures_at_vertex(pmesh, v, np.vertex_principal_curvatures_and_directions(pcd));
+  return *pcd;
+}
+
+template<typename PolygonMesh,
+  typename NamedParameters = parameters::Default_named_parameters>
+  void interpolated_corrected_curvatures_at_vertex(const PolygonMesh& pmesh,
+    typename boost::graph_traits<PolygonMesh>::vertex_descriptor v,
+    const NamedParameters& np = parameters::default_values())
+{
+  internal::interpolated_corrected_curvatures_one_vertex(pmesh, v, np);
+}
 
 } // namespace Polygon_mesh_processing
 } // namespace CGAL

STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h

@@ -90,6 +90,9 @@ CGAL_add_named_parameter(nb_points_per_distance_unit_t, nb_points_per_distance_u
 CGAL_add_named_parameter(vertex_mean_curvature_map_t, vertex_mean_curvature_map, vertex_mean_curvature_map)
 CGAL_add_named_parameter(vertex_gaussian_curvature_map_t, vertex_gaussian_curvature_map, vertex_gaussian_curvature_map)
 CGAL_add_named_parameter(vertex_principal_curvatures_and_directions_map_t, vertex_principal_curvatures_and_directions_map, vertex_principal_curvatures_and_directions_map)
+CGAL_add_named_parameter(vertex_mean_curvature_t, vertex_mean_curvature, vertex_mean_curvature)
+CGAL_add_named_parameter(vertex_gaussian_curvature_t, vertex_gaussian_curvature, vertex_gaussian_curvature)
+CGAL_add_named_parameter(vertex_principal_curvatures_and_directions_t, vertex_principal_curvatures_and_directions, vertex_principal_curvatures_and_directions)
 CGAL_add_named_parameter(ball_radius_t, ball_radius, ball_radius)
 CGAL_add_named_parameter(outward_orientation_t, outward_orientation, outward_orientation)
 CGAL_add_named_parameter(overlap_test_t, overlap_test, do_overlap_test_of_bounded_sides)