Merge pull request #6088 from sloriot/PMP-add_discrete_curvature

Add functions to compute discrete curvatures

Installation/CHANGES.md

@@ -6,6 +6,11 @@
 ### General Changes
 - The minimal supported version of Boost is now 1.74.0.
 
+### [Polygon Mesh Processing](https://doc.cgal.org/6.1/Manual/packages.html#PkgPolygonMeshProcessing)
+- Added the function `CGAL::Polygon_mesh_processing::discrete_mean_curvature` and `CGAL::Polygon_mesh_processing::discrete_Guassian_curvature` to evaluate the discrete curvature at a vertex of a mesh.
+- Added the function `CGAL::Polygon_mesh_processing::angle_sum` to compute the sum of the angles around a vertex.
+
+
 ### [Algebraic Kernel](https://doc.cgal.org/6.1/Manual/packages.html#PkgAlgebraicKernelD)
 
 -   **Breaking change**: Classes based on the RS Library are no longer provided.

Kernel_23/include/CGAL/Kernel/function_objects.h

@@ -965,22 +965,14 @@ namespace CommonKernelFunctors {
      typename K::Compute_scalar_product_3 scalar_product =
        k.compute_scalar_product_3_object();
 
-     double product = CGAL::sqrt(to_double(scalar_product(u,u)) * to_double(scalar_product(v,v)));
+     double product = to_double(approximate_sqrt(scalar_product(u,u) * scalar_product(v,v)));
 
      if(product == 0)
        return 0;
 
      // cosine
      double dot = to_double(scalar_product(u,v));
-     double cosine = dot / product;
+     double cosine = std::clamp(dot / product, -1., 1.);
-
-     if(cosine > 1.){
-       cosine = 1.;
-     }
-     if(cosine < -1.){
-       cosine = -1.;
-     }
-
      return std::acos(cosine) * 180./CGAL_PI;
    }
 

Lab/demo/Lab/Color_map.h

@@ -23,9 +23,9 @@ compute_color_map(QColor base_color,
                   std::size_t nb_of_colors,
                   Output_color_iterator out)
 {
-  qreal hue = base_color.hueF();
+  const qreal step = (static_cast<qreal>(0.85)) / nb_of_colors;
-  const qreal step = (static_cast<qreal>(1)) / nb_of_colors;
 
+  qreal hue = base_color.hueF();
   qreal h = (hue == -1) ? 0 : hue;
   for(std::size_t i=0; i<nb_of_colors; ++i)
   {

Lab/demo/Lab/Plugins/Display/Display_property_plugin.cpp

@@ -22,6 +22,7 @@
 #include <CGAL/Polygon_mesh_processing/compute_normal.h>
 #include <CGAL/Polygon_mesh_processing/measure.h>
 #include <CGAL/Polygon_mesh_processing/triangulate_faces.h>
+#include <CGAL/Polygon_mesh_processing/curvature.h>
 #include <CGAL/Polygon_mesh_processing/interpolated_corrected_curvatures.h>
 
 #include <QAbstractItemView>
@@ -350,11 +351,11 @@ private:
   template<typename ValueType>
   void displayMapLegend(const std::vector<ValueType>& values)
   {
-    const std::size_t size = (std::min)(color_map.size(), std::size_t(1024));
+    const std::size_t size = (std::min)(color_map.size(), std::size_t(4096));
 
     const int text_height = 20;
     const int height = text_height * static_cast<int>(size) + text_height;
-    const int width = 140;
+    const int width = 200;
     const int cell_width = width / 3;
     const int top_margin = 15;
     const int left_margin = 5;
@@ -381,13 +382,13 @@ private:
                        tick_height,
                        color);
 
-      QRect text_rect(left_margin + cell_width + 10, drawing_height - top_margin - j, 50, text_height);
+      QRect text_rect(left_margin + cell_width + 10, drawing_height - top_margin - j, 100, text_height);
-      painter.drawText(text_rect, Qt::AlignCenter, QObject::tr("%1").arg(values[i], 0, 'f', 3, QLatin1Char(' ')));
+      painter.drawText(text_rect, Qt::AlignCenter, QObject::tr("%1").arg(values[i], 0, 'f', 6, QLatin1Char(' ')));
     }
 
     if(color_map.size() > size)
     {
-      QRect text_rect(left_margin + cell_width + 10, 0, 50, text_height);
+      QRect text_rect(left_margin + cell_width + 10, 0, 100, text_height);
       painter.drawText(text_rect, Qt::AlignCenter, QObject::tr("[...]"));
     }
 
@@ -463,6 +464,8 @@ private:
                                           "Largest Angle Per Face",
                                           "Scaled Jacobian",
                                           "Face Area",
+                                          "Discrete Mean Curvature",
+                                          "Discrete Gaussian Curvature",
                                           "Interpolated Corrected Mean Curvature",
                                           "Interpolated Corrected Gaussian Curvature"});
       property_simplex_types = { Property_simplex_type::FACE,
@@ -470,6 +473,8 @@ private:
                                  Property_simplex_type::FACE,
                                  Property_simplex_type::FACE,
                                  Property_simplex_type::VERTEX,
+                                 Property_simplex_type::VERTEX,
+                                 Property_simplex_type::VERTEX,
                                  Property_simplex_type::VERTEX };
       detectSMScalarProperties(*(sm_item->face_graph()));
     }
@@ -516,12 +521,12 @@ private Q_SLOTS:
 
       // Curvature property-specific slider
       const std::string& property_name = dock_widget->propertyBox->currentText().toStdString();
-      const bool is_curvature_property = (property_name == "Interpolated Corrected Mean Curvature" ||
+      const bool is_interpolated_curvature_property = (property_name == "Interpolated Corrected Mean Curvature" ||
-                                          property_name == "Interpolated Corrected Gaussian Curvature");
+                                                       property_name == "Interpolated Corrected Gaussian Curvature");
-      dock_widget->expandingRadiusLabel->setVisible(is_curvature_property);
+      dock_widget->expandingRadiusLabel->setVisible(is_interpolated_curvature_property);
-      dock_widget->expandingRadiusSlider->setVisible(is_curvature_property);
+      dock_widget->expandingRadiusSlider->setVisible(is_interpolated_curvature_property);
-      dock_widget->expandingRadiusLabel->setEnabled(is_curvature_property);
+      dock_widget->expandingRadiusLabel->setEnabled(is_interpolated_curvature_property);
-      dock_widget->expandingRadiusSlider->setEnabled(is_curvature_property);
+      dock_widget->expandingRadiusSlider->setEnabled(is_interpolated_curvature_property);
     }
     else // no or broken property
     {
@@ -570,6 +575,16 @@ private:
     {
       displayArea(sm_item);
     }
+    else if(property_name == "Discrete Mean Curvature")
+    {
+      displayDiscreteCurvatureMeasure(sm_item, MEAN_CURVATURE);
+      sm_item->setRenderingMode(Gouraud);
+    }
+    else if(property_name == "Discrete Gaussian Curvature")
+    {
+      displayDiscreteCurvatureMeasure(sm_item, GAUSSIAN_CURVATURE);
+      sm_item->setRenderingMode(Gouraud);
+    }
     else if(property_name == "Interpolated Corrected Mean Curvature")
     {
       displayInterpolatedCurvatureMeasure(sm_item, MEAN_CURVATURE);
@@ -682,6 +697,8 @@ private:
     removeDisplayPluginProperty(item, "f:display_plugin_largest_angle");
     removeDisplayPluginProperty(item, "f:display_plugin_scaled_jacobian");
     removeDisplayPluginProperty(item, "f:display_plugin_area");
+    removeDisplayPluginProperty(item, "v:display_plugin_discrete_mean_curvature");
+    removeDisplayPluginProperty(item, "v:display_plugin_discrete_Gaussian_curvature");
     removeDisplayPluginProperty(item, "v:display_plugin_interpolated_corrected_mean_curvature");
     removeDisplayPluginProperty(item, "v:display_plugin_interpolated_corrected_Gaussian_curvature");
   }
@@ -864,6 +881,35 @@ private:
     displaySMProperty<face_descriptor>("f:display_plugin_area", *sm);
   }
 
+private:
+  void displayDiscreteCurvatureMeasure(Scene_surface_mesh_item* sm_item,
+                                       CurvatureType mu_index)
+  {
+    SMesh* sm = sm_item->face_graph();
+    if(sm == nullptr)
+      return;
+
+    if(mu_index != MEAN_CURVATURE && mu_index != GAUSSIAN_CURVATURE)
+      return;
+
+    std::string vdc_name = (mu_index == MEAN_CURVATURE) ? "v:display_plugin_discrete_mean_curvature"
+                                                        : "v:display_plugin_discrete_Gaussian_curvature";
+
+    bool not_initialized;
+    SMesh::Property_map<vertex_descriptor, double> vdc;
+    std::tie(vdc, not_initialized) = sm->add_property_map<vertex_descriptor, double>(vdc_name, 0);
+
+    if(not_initialized)
+    {
+      if(mu_index == MEAN_CURVATURE)
+        PMP::discrete_mean_curvatures(*sm, vdc);
+      else
+        PMP::discrete_Gaussian_curvatures(*sm, vdc);
+    }
+
+    displaySMProperty<vertex_descriptor>(vdc_name, *sm);
+  }
+
 private Q_SLOTS:
   void setExpandingRadius()
   {
@@ -1131,6 +1177,10 @@ private:
         zoomToSimplexWithPropertyExtremum(faces(mesh), mesh, "f:display_plugin_scaled_jacobian", extremum);
       else if(property_name == "Face Area")
         zoomToSimplexWithPropertyExtremum(faces(mesh), mesh, "f:display_plugin_area", extremum);
+      else if(property_name == "Discrete Mean Curvature")
+        zoomToSimplexWithPropertyExtremum(vertices(mesh), mesh, "v:display_plugin_discrete_mean_curvature", extremum);
+      else if(property_name == "Discrete Gaussian Curvature")
+        zoomToSimplexWithPropertyExtremum(vertices(mesh), mesh, "v:display_plugin_discrete_Gaussian_curvature", extremum);
       else if(property_name == "Interpolated Corrected Mean Curvature")
         zoomToSimplexWithPropertyExtremum(vertices(mesh), mesh, "v:display_plugin_interpolated_corrected_mean_curvature", extremum);
       else if(property_name == "Interpolated Corrected Gaussian Curvature")
@@ -1470,6 +1520,8 @@ isSMPropertyScalar(const std::string& name,
      name == "f:display_plugin_largest_angle" ||
      name == "f:display_plugin_scaled_jacobian" ||
      name == "f:display_plugin_area" ||
+     name == "v:display_plugin_discrete_mean_curvature" ||
+     name == "v:display_plugin_discrete_Gaussian_curvature" ||
      name == "v:display_plugin_interpolated_corrected_mean_curvature" ||
      name == "v:display_plugin_interpolated_corrected_Gaussian_curvature")
     return false;

Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt

@@ -25,7 +25,7 @@
 /// \ingroup PkgPolygonMeshProcessingRef
 
 /// \defgroup PMP_measure_grp Geometric Measure Functions
-/// Functions to compute lengths of edges and borders, areas of faces and patches, as well as volumes of closed meshes.
+/// Functions to compute discrete curvatures, lengths of edges and borders, areas of faces and patches, volumes of closed meshes.
 /// \ingroup PkgPolygonMeshProcessingRef
 
 /// \defgroup PMP_orientation_grp Orientation Functions
@@ -239,6 +239,11 @@ The page \ref bgl_namedparameters "Named Parameters" describes their usage.
 - `CGAL::Polygon_mesh_processing::sample_triangle_mesh()`
 
 \cgalCRPSection{Geometric Measure Functions}
+- \link PMP_measure_grp `CGAL::Polygon_mesh_processing::angle_sum()` \endlink
+- \link PMP_measure_grp `CGAL::Polygon_mesh_processing::discrete_Gaussian_curvatures()` \endlink
+- \link PMP_measure_grp `CGAL::Polygon_mesh_processing::discrete_Gaussian_curvature()` \endlink
+- \link PMP_measure_grp `CGAL::Polygon_mesh_processing::discrete_mean_curvature()` \endlink
+- \link PMP_measure_grp `CGAL::Polygon_mesh_processing::discrete_mean_curvatures()` \endlink
 - \link PMP_measure_grp `CGAL::Polygon_mesh_processing::edge_length()` \endlink
 - \link PMP_measure_grp `CGAL::Polygon_mesh_processing::squared_edge_length()` \endlink
 - \link PMP_measure_grp `CGAL::Polygon_mesh_processing::face_area()` \endlink
@@ -248,7 +253,6 @@ The page \ref bgl_namedparameters "Named Parameters" describes their usage.
 - \link PMP_measure_grp `CGAL::Polygon_mesh_processing::face_border_length()` \endlink
 - \link PMP_measure_grp `CGAL::Polygon_mesh_processing::longest_border()` \endlink
 - \link PMP_measure_grp `CGAL::Polygon_mesh_processing::centroid()` \endlink
-- \link PMP_measure_grp `CGAL::Polygon_mesh_processing::match_faces()` \endlink
 
 \cgalCRPSection{Feature Detection Functions}
 - `CGAL::Polygon_mesh_processing::sharp_edges_segmentation()`

Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt

@@ -1225,6 +1225,17 @@ compute the curvatures on a specific vertex.
 
 \cgalExample{Polygon_mesh_processing/interpolated_corrected_curvatures_vertex.cpp}
 
+\subsection DCurvartures Discrete Curvatures
+
+The package also provides methods to compute the standard, non-interpolated discrete mean and Gaussian
+curvatures on triangle meshes, based on the work of Meyer et al. \cgalCite{cgal:mdsb-ddgot-02}.
+These curvatures are computed at each vertex of the mesh, and are based on the angles of the incident
+triangles. The functions are:
+- `CGAL::Polygon_mesh_processing::discrete_mean_curvature()`
+- `CGAL::Polygon_mesh_processing::discrete_mean_curvatures()`
+- `CGAL::Polygon_mesh_processing::discrete_Gaussian_curvature()`
+- `CGAL::Polygon_mesh_processing::discrete_Gaussian_curvatures()`
+
 ****************************************
 \section PMPSlicer Slicer
 

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/curvature.h

@@ -0,0 +1,475 @@
+// Copyright (c) 2021 GeometryFactory (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
+//
+//
+// Author(s)     : Andreas Fabri,
+//                 Mael Rouxel-Labbé
+
+#ifndef CGAL_PMP_CURVATURE_H
+#define CGAL_PMP_CURVATURE_H
+
+#include <CGAL/license/Polygon_mesh_processing/measure.h>
+
+#include <CGAL/boost/graph/named_params_helper.h>
+#include <CGAL/Named_function_parameters.h>
+#include <CGAL/Polygon_mesh_processing/measure.h>
+#include <CGAL/Weights/cotangent_weights.h>
+
+#include <cmath>
+#include <algorithm>
+
+namespace CGAL {
+namespace Polygon_mesh_processing {
+
+/**
+  * \ingroup PMP_measure_grp
+  *
+  * computes the sum of the angles around a vertex.
+  *
+  * The angle sum is given in degrees.
+  *
+  * @tparam PolygonMesh a model of `FaceGraph`
+  * @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
+  *
+  * @param v the vertex whose sum of angles is computed
+  * @param pmesh the polygon mesh to which `v` belongs
+  * @param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
+  *
+  * \cgalNamedParamsBegin
+  *   \cgalParamNBegin{vertex_point_map}
+  *     \cgalParamDescription{a property map associating points to the vertices of `pmesh`}
+  *     \cgalParamType{a class model of `ReadablePropertyMap` with `boost::graph_traits<PolygonMesh>::%vertex_descriptor`
+  *                    as key type and `%Point_3` as value type}
+  *     \cgalParamDefault{`boost::get(CGAL::vertex_point, pmesh)`}
+  *   \cgalParamNEnd
+  *
+ *   \cgalParamNBegin{geom_traits}
+ *     \cgalParamDescription{an instance of a geometric traits class}
+ *     \cgalParamType{The traits class must provide the nested functor `Compute_approximate_angle_3`,
+ *                    model of `Kernel::ComputeApproximateAngle_3`.}
+ *     \cgalParamDefault{a \cgal kernel deduced from the point type, using `CGAL::Kernel_traits`}
+ *     \cgalParamExtra{The geometric traits class must be compatible with the vertex point type.}
+ *   \cgalParamNEnd
+  * \cgalNamedParamsEnd
+  *
+  * @return the sum of angles around `v`. The return type `FT` is a number type either deduced
+  * from the `geom_traits` \ref bgl_namedparameters "Named Parameters" if provided,
+  * or the geometric traits class deduced from the point property map of `pmesh`.
+  *
+  * \warning This function involves trigonometry.
+  */
+template<typename PolygonMesh,
+         typename CGAL_NP_TEMPLATE_PARAMETERS>
+#ifdef DOXYGEN_RUNNING
+FT
+#else
+typename GetGeomTraits<PolygonMesh, CGAL_NP_CLASS>::type::FT
+#endif
+angle_sum(typename boost::graph_traits<PolygonMesh>::vertex_descriptor v,
+          const PolygonMesh& pmesh,
+          const CGAL_NP_CLASS& np = parameters::default_values())
+{
+  using parameters::choose_parameter;
+  using parameters::get_parameter;
+
+  using Geom_traits = typename GetGeomTraits<PolygonMesh, CGAL_NP_CLASS>::type;
+  using FT = typename Geom_traits::FT;
+
+  typename GetVertexPointMap<PolygonMesh, CGAL_NP_CLASS>::const_type
+      vpm = choose_parameter(get_parameter(np, internal_np::vertex_point),
+                             get_const_property_map(CGAL::vertex_point, pmesh));
+
+  Geom_traits gt = choose_parameter<Geom_traits>(get_parameter(np, internal_np::geom_traits));
+
+  CGAL_precondition(is_valid_vertex_descriptor(v, pmesh));
+
+  typename Geom_traits::Compute_approximate_angle_3 approx_angle = gt.compute_approximate_angle_3_object();
+
+  FT angle_sum = 0;
+  for(auto h : halfedges_around_source(v, pmesh))
+  {
+    if(is_border(h, pmesh))
+      continue;
+
+    angle_sum += approx_angle(get(vpm, target(h, pmesh)),
+                              get(vpm, source(h, pmesh)),
+                              get(vpm, source(prev(h,pmesh), pmesh)));
+  }
+
+  return angle_sum;
+}
+
+// Discrete Gaussian Curvature
+
+/**
+  * \ingroup PMP_measure_grp
+  *
+  * computes the discrete Gaussian curvature at a vertex.
+  *
+  * We refer to Meyer et al. \cgalCite{cgal:mdsb-ddgot-02} for the definition of <i>discrete Gaussian curvature</i>.
+  *
+  * @tparam TriangleMesh a model of `FaceGraph`
+  * @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
+  *
+  * @param v the vertex whose discrete Gaussian curvature is being computed
+  * @param tmesh the triangle mesh to which `v` belongs
+  * @param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
+  *
+  * \cgalNamedParamsBegin
+  *   \cgalParamNBegin{vertex_point_map}
+  *     \cgalParamDescription{a property map associating points to the vertices of `tmesh`}
+  *     \cgalParamType{a class model of `ReadablePropertyMap` with `boost::graph_traits<TriangleMesh>::%vertex_descriptor`
+  *                    as key type and `%Point_3` as value type}
+  *     \cgalParamDefault{`boost::get(CGAL::vertex_point, tmesh)`}
+  *   \cgalParamNEnd
+  *
+ *   \cgalParamNBegin{geom_traits}
+ *     \cgalParamDescription{an instance of a geometric traits class}
+ *     \cgalParamType{The traits class must be a model of `Kernel`.}
+ *     \cgalParamDefault{a \cgal kernel deduced from the point type, using `CGAL::Kernel_traits`}
+ *     \cgalParamExtra{The geometric traits class must be compatible with the vertex point type.}
+ *   \cgalParamNEnd
+  * \cgalNamedParamsEnd
+  *
+  * @return the discrete Gaussian curvature at `v`. The return type `FT` is a number type either deduced
+  * from the `geom_traits` \ref bgl_namedparameters "Named Parameters" if provided,
+  * or the geometric traits class deduced from the point property map of `tmesh`.
+  *
+  * \warning This function involves trigonometry.
+  * \warning The current formulation is not well defined for border vertices.
+  *
+  * \pre `tmesh` is a triangle mesh
+  */
+template <typename TriangleMesh,
+          typename CGAL_NP_TEMPLATE_PARAMETERS>
+#ifdef DOXYGEN_RUNNING
+FT
+#else
+typename GetGeomTraits<TriangleMesh, CGAL_NP_CLASS>::type::FT
+#endif
+discrete_Gaussian_curvature(typename boost::graph_traits<TriangleMesh>::vertex_descriptor v,
+                            const TriangleMesh& tmesh,
+                            const CGAL_NP_CLASS& np = parameters::default_values())
+{
+  using parameters::choose_parameter;
+  using parameters::get_parameter;
+
+  using GeomTraits = typename GetGeomTraits<TriangleMesh, CGAL_NP_CLASS>::type;
+  using FT = typename GeomTraits::FT;
+  using Vector_3 = typename GeomTraits::Vector_3;
+
+  using VertexPointMap = typename GetVertexPointMap<TriangleMesh, CGAL_NP_CLASS>::const_type;
+  using halfedge_descriptor = typename boost::graph_traits<TriangleMesh>::halfedge_descriptor;
+
+  GeomTraits gt = choose_parameter<GeomTraits>(get_parameter(np, internal_np::geom_traits));
+  VertexPointMap vpm = choose_parameter(get_parameter(np, internal_np::vertex_point),
+                                        get_const_property_map(vertex_point, tmesh));
+
+  typename GeomTraits::Construct_vector_3 vector =
+    gt.construct_vector_3_object();
+  typename GeomTraits::Construct_cross_product_vector_3 cross_product =
+    gt.construct_cross_product_vector_3_object();
+  typename GeomTraits::Compute_scalar_product_3 scalar_product =
+    gt.compute_scalar_product_3_object();
+  typename GeomTraits::Compute_squared_length_3 squared_length =
+    gt.compute_squared_length_3_object();
+
+  FT angle_sum = 0;
+
+  for(halfedge_descriptor h : CGAL::halfedges_around_target(v, tmesh))
+  {
+    if(is_border(h, tmesh))
+      continue;
+
+    const Vector_3 v0 = vector(get(vpm, v), get(vpm, target(next(h, tmesh), tmesh))); // p1p2
+    const Vector_3 v1 = vector(get(vpm, v), get(vpm, source(h, tmesh))); // p1p0
+
+    const FT dot = scalar_product(v0, v1);
+    const Vector_3 cross = cross_product(v0, v1);
+    const FT sqcn = squared_length(cross);
+    if(is_zero(dot))
+    {
+      angle_sum += CGAL_PI/FT(2);
+    }
+    else
+    {
+      if(is_zero(sqcn)) // collinear
+      {
+        if(dot < 0)
+          angle_sum += CGAL_PI;
+        // else
+        //   angle_sum += 0;
+      }
+      else
+      {
+        angle_sum += std::atan2(CGAL::approximate_sqrt(sqcn), dot);
+      }
+    }
+  }
+
+  Weights::Secure_cotangent_weight_with_voronoi_area<TriangleMesh, VertexPointMap, GeomTraits> wc(tmesh, vpm, gt);
+
+  const FT gaussian_curvature = (2 * CGAL_PI - angle_sum) / wc.voronoi(v);
+
+  return gaussian_curvature;
+}
+
+/**
+  * \ingroup PMP_measure_grp
+  *
+  * computes the discrete Gaussian curvatures at the vertices of a mesh.
+  *
+  * We refer to Meyer et al. \cgalCite{cgal:mdsb-ddgot-02} for the definition of <i>discrete Gaussian curvature</i>.
+  *
+  * @tparam TriangleMesh a model of `FaceGraph`
+  * @tparam VertexCurvatureMap must be a model of `WritablePropertyMap` with key type
+  *                            `boost::graph_traits<TriangleMesh>::%vertex_descriptor` and value type `FT`,
+  *                            which is either `geom_traits::FT` if this named parameter is provided,
+  *                            or `kernel::FT` with the kernel deduced from from the point property map of `tmesh`.
+  * @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
+  *
+  * @param tmesh the triangle mesh to which `v` belongs
+  * @param vcm the property map that contains the computed discrete curvatures
+  * @param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
+  *
+  * \cgalNamedParamsBegin
+  *   \cgalParamNBegin{vertex_point_map}
+  *     \cgalParamDescription{a property map associating points to the vertices of `tmesh`}
+  *     \cgalParamType{a class model of `ReadablePropertyMap` with `boost::graph_traits<TriangleMesh>::%vertex_descriptor`
+  *                    as key type and `%Point_3` as value type}
+  *     \cgalParamDefault{`boost::get(CGAL::vertex_point, tmesh)`}
+  *   \cgalParamNEnd
+  *
+ *   \cgalParamNBegin{geom_traits}
+ *     \cgalParamDescription{an instance of a geometric traits class}
+ *     \cgalParamType{The traits class must be a model of `Kernel`.}
+ *     \cgalParamDefault{a \cgal kernel deduced from the point type, using `CGAL::Kernel_traits`}
+ *     \cgalParamExtra{The geometric traits class must be compatible with the vertex point type.}
+ *   \cgalParamNEnd
+  * \cgalNamedParamsEnd
+  *
+  * \warning This function involves trigonometry.
+  * \warning The current formulation is not well defined for border vertices.
+  *
+  * \pre `tmesh` is a triangle mesh
+  */
+template <typename TriangleMesh,
+          typename VertexCurvatureMap,
+          typename CGAL_NP_TEMPLATE_PARAMETERS>
+void discrete_Gaussian_curvatures(const TriangleMesh& tmesh,
+                                  VertexCurvatureMap vcm,
+                                  const CGAL_NP_CLASS& np = parameters::default_values())
+{
+  using vertex_descriptor = typename boost::graph_traits<TriangleMesh>::vertex_descriptor;
+
+  for(vertex_descriptor v : vertices(tmesh))
+  {
+    put(vcm, v, discrete_Gaussian_curvature(v, tmesh, np));
+    // std::cout << "curvature: " << get(vcm, v) << std::endl;
+  }
+}
+
+// Discrete Mean Curvature
+
+/**
+  * \ingroup PMP_measure_grp
+  *
+  * computes the discrete mean curvature at a vertex.
+  *
+  * We refer to Meyer et al. \cgalCite{cgal:mdsb-ddgot-02} for the definition of <i>discrete mean curvature</i>.
+  *
+  * @tparam TriangleMesh a model of `FaceGraph`
+  * @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
+  *
+  * @param v the vertex whose discrete mean curvature is being computed
+  * @param tmesh the triangle mesh to which `v` belongs
+  * @param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
+  *
+  * \cgalNamedParamsBegin
+  *   \cgalParamNBegin{vertex_point_map}
+  *     \cgalParamDescription{a property map associating points to the vertices of `tmesh`}
+  *     \cgalParamType{a class model of `ReadablePropertyMap` with `boost::graph_traits<TriangleMesh>::%vertex_descriptor`
+  *                    as key type and `%Point_3` as value type}
+  *     \cgalParamDefault{`boost::get(CGAL::vertex_point, tmesh)`}
+  *   \cgalParamNEnd
+  *
+ *   \cgalParamNBegin{geom_traits}
+ *     \cgalParamDescription{an instance of a geometric traits class}
+ *     \cgalParamType{The traits class must be a model of `Kernel`.}
+ *     \cgalParamDefault{a \cgal kernel deduced from the point type, using `CGAL::Kernel_traits`}
+ *     \cgalParamExtra{The geometric traits class must be compatible with the vertex point type.}
+ *   \cgalParamNEnd
+  * \cgalNamedParamsEnd
+  *
+  * @return the discrete mean curvature at `v`. The return type `FT` is a number type either deduced
+  * from the `geom_traits` \ref bgl_namedparameters "Named Parameters" if provided,
+  * or the geometric traits class deduced from the point property map of `tmesh`.
+  *
+  * \warning The current formulation is not well defined for border vertices.
+  *
+  * \pre `tmesh` is a triangle mesh
+  */
+template <typename TriangleMesh,
+          typename CGAL_NP_TEMPLATE_PARAMETERS>
+#ifdef DOXYGEN_RUNNING
+FT
+#else
+typename GetGeomTraits<TriangleMesh, CGAL_NP_CLASS>::type::FT
+#endif
+discrete_mean_curvature(typename boost::graph_traits<TriangleMesh>::vertex_descriptor v,
+                        const TriangleMesh& tmesh,
+                        const CGAL_NP_CLASS& np = parameters::default_values())
+{
+  using parameters::choose_parameter;
+  using parameters::get_parameter;
+
+  using GeomTraits = typename GetGeomTraits<TriangleMesh, CGAL_NP_CLASS>::type;
+  using FT = typename GeomTraits::FT;
+  using Vector_3 = typename GeomTraits::Vector_3;
+
+  using VertexPointMap = typename GetVertexPointMap<TriangleMesh, CGAL_NP_CLASS>::const_type;
+  using Point_ref = typename boost::property_traits<VertexPointMap>::reference;
+
+  using vertex_descriptor = typename boost::graph_traits<TriangleMesh>::vertex_descriptor;
+  using halfedge_descriptor = typename boost::graph_traits<TriangleMesh>::halfedge_descriptor;
+
+  GeomTraits gt = choose_parameter<GeomTraits>(get_parameter(np, internal_np::geom_traits));
+  VertexPointMap vpm = choose_parameter(get_parameter(np, internal_np::vertex_point),
+                                        get_const_property_map(vertex_point, tmesh));
+
+#if 0
+  typename GeomTraits::Compute_squared_distance_3 squared_distance =
+    gt.compute_squared_distance_3_object();
+  typename GeomTraits::Compute_approximate_dihedral_angle_3 dihedral_angle =
+    gt.compute_approximate_dihedral_angle_3_object();
+
+  const FT two_pi = 2 * CGAL_PI;
+
+  FT hi = 0;
+  for(halfedge_descriptor h : CGAL::halfedges_around_target(v, tmesh))
+  {
+    const Point_3& p = get(vpm, source(h, tmesh));
+    const Point_3& q = get(vpm, target(h, tmesh));
+    const Point_3& r = get(vpm, target(next(h, tmesh), tmesh));
+    const Point_3& s = get(vpm, target(next(opposite(h, tmesh), tmesh), tmesh));
+    const FT l = squared_distance(p,q);
+
+    FT phi = CGAL_PI * dihedral_angle(p, q, r, s) / FT(180);
+
+    if(phi < 0)
+      phi += two_pi;
+    if(phi > two_pi)
+      phi = two_pi;
+
+    hi += FT(0.5) * l * (CGAL_PI - phi);
+  }
+
+  return FT(0.5) * hi;
+#else
+  typename GeomTraits::Construct_vector_3 vector =
+    gt.construct_vector_3_object();
+  typename GeomTraits::Construct_sum_of_vectors_3 vector_sum =
+    gt.construct_sum_of_vectors_3_object();
+  typename GeomTraits::Construct_scaled_vector_3 scaled_vector =
+    gt.construct_scaled_vector_3_object();
+  typename GeomTraits::Compute_squared_length_3 squared_length =
+    gt.compute_squared_length_3_object();
+
+  Weights::Secure_cotangent_weight_with_voronoi_area<TriangleMesh, VertexPointMap, GeomTraits> wc(tmesh, vpm, gt);
+
+  Vector_3 kh = vector(CGAL::NULL_VECTOR);
+  for(halfedge_descriptor h : CGAL::halfedges_around_target(v, tmesh))
+  {
+    const vertex_descriptor v1 = source(h, tmesh);
+
+    const Point_ref p0 = get(vpm, v);
+    const Point_ref p1 = get(vpm, v1);
+
+    FT local_c = 0;
+    if(!is_border(h, tmesh))
+    {
+      const vertex_descriptor v2 = target(next(h, tmesh), tmesh);
+      const Point_ref p2 = get(vpm, v2);
+      local_c += Weights::cotangent_3_clamped(p0, p2, p1, gt);
+    }
+
+    if(!is_border(opposite(h, tmesh), tmesh))
+    {
+      const vertex_descriptor v3 = target(next(opposite(h, tmesh), tmesh), tmesh);
+      const Point_ref p3 = get(vpm, v3);
+      local_c += Weights::cotangent_3_clamped(p1, p3, p0, gt);
+    }
+
+    kh = vector_sum(kh, scaled_vector(vector(p0, p1), local_c));
+  }
+
+  const FT khn = CGAL::approximate_sqrt(squared_length(kh));
+  const FT va = wc.voronoi(v);
+  CGAL_assertion(!is_zero(va));
+
+  const FT mean_curvature = khn / (FT(4) * va);
+  return mean_curvature;
+#endif
+}
+
+/**
+  * \ingroup PMP_measure_grp
+  *
+  * computes the discrete mean curvatures at the vertices of a mesh.
+  *
+  * We refer to Meyer et al. \cgalCite{cgal:mdsb-ddgot-02} for the definition of <i>discrete mean curvature</i>.
+  *
+  * @tparam TriangleMesh a model of `FaceGraph`
+  * @tparam VertexCurvatureMap must be a model of `WritablePropertyMap` with key type
+  *                            `boost::graph_traits<TriangleMesh>::%vertex_descriptor` and value type `FT`,
+  *                            which is either `geom_traits::FT` if this named parameter is provided,
+  *                            or `kernel::FT` with the kernel deduced from from the point property map of `tmesh`.
+  * @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
+  *
+  * @param tmesh the triangle mesh to which `v` belongs
+  * @param vcm the property map that contains the computed discrete curvatures
+  * @param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
+  *
+  * \cgalNamedParamsBegin
+  *   \cgalParamNBegin{vertex_point_map}
+  *     \cgalParamDescription{a property map associating points to the vertices of `tmesh`}
+  *     \cgalParamType{a class model of `ReadablePropertyMap` with `boost::graph_traits<TriangleMesh>::%vertex_descriptor`
+  *                    as key type and `%Point_3` as value type}
+  *     \cgalParamDefault{`boost::get(CGAL::vertex_point, tmesh)`}
+  *   \cgalParamNEnd
+  *
+ *   \cgalParamNBegin{geom_traits}
+ *     \cgalParamDescription{an instance of a geometric traits class}
+ *     \cgalParamType{The traits class must be a model of `Kernel`.}
+ *     \cgalParamDefault{a \cgal kernel deduced from the point type, using `CGAL::Kernel_traits`}
+ *     \cgalParamExtra{The geometric traits class must be compatible with the vertex point type.}
+ *   \cgalParamNEnd
+  * \cgalNamedParamsEnd
+  *
+  * \warning The current formulation is not well defined for border vertices.
+  *
+  * \pre `tmesh` is a triangle mesh
+  */
+template <typename TriangleMesh,
+          typename VertexCurvatureMap,
+          typename CGAL_NP_TEMPLATE_PARAMETERS>
+void discrete_mean_curvatures(const TriangleMesh& tmesh,
+                              VertexCurvatureMap vcm,
+                              const CGAL_NP_CLASS& np = parameters::default_values())
+{
+  using vertex_descriptor = typename boost::graph_traits<TriangleMesh>::vertex_descriptor;
+
+  for(vertex_descriptor v : vertices(tmesh))
+    put(vcm, v, discrete_mean_curvature(v, tmesh, np));
+}
+
+} // namespace Polygon_mesh_processing
+} // namespace CGAL
+
+#endif //CGAL_PMP_CURVATURE_H

Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt

@@ -28,6 +28,7 @@ create_single_source_cgal_program("test_stitching.cpp")
 create_single_source_cgal_program("remeshing_test.cpp")
 create_single_source_cgal_program("remeshing_with_isolated_constraints_test.cpp" )
 create_single_source_cgal_program("measures_test.cpp")
+create_single_source_cgal_program("test_discrete_curvatures.cpp")
 create_single_source_cgal_program("triangulate_faces_test.cpp")
 create_single_source_cgal_program("triangulate_faces_hole_filling_dt3_test.cpp")
 create_single_source_cgal_program("triangulate_faces_hole_filling_all_search_test.cpp")

Polygon_mesh_processing/test/Polygon_mesh_processing/test_discrete_curvatures.cpp

@@ -0,0 +1,146 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/Polyhedron_3.h>
+
+#include <CGAL/Polygon_mesh_processing/curvature.h>
+
+#include <boost/graph/graph_traits.hpp>
+
+#include <iostream>
+#include <string>
+
+#define ABS_ERROR 1e-6
+
+namespace PMP = CGAL::Polygon_mesh_processing;
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+typedef K::FT FT;
+typedef CGAL::Surface_mesh<K::Point_3> SMesh;
+typedef CGAL::Polyhedron_3<K> Polyhedron;
+
+struct Average_test_info
+{
+  FT mean_curvature_avg;
+  FT gaussian_curvature_avg;
+  FT tolerance = 0.9;
+
+  Average_test_info(FT mean_curvature_avg,
+                    FT gaussian_curvature_avg)
+    : mean_curvature_avg(mean_curvature_avg),
+      gaussian_curvature_avg(gaussian_curvature_avg)
+  { }
+};
+
+bool passes_comparison(FT result, FT expected, FT tolerance)
+{
+  std::cout << "result: " << result << std::endl;
+  std::cout << "expected: " << expected << std::endl;
+
+  if(abs(expected) < ABS_ERROR && abs(result) < ABS_ERROR)
+    return true; // expected 0, got 0
+  else if (abs(expected) < ABS_ERROR)
+    return false; // expected 0, got non-0
+
+  return (std::min)(result, expected) / (std::max)(result, expected) > tolerance;
+}
+
+template <typename TriangleMesh>
+void test_curvatures(std::string mesh_path,
+                     Average_test_info test_info)
+{
+  std::cout << "test discrete curvatures of " << mesh_path << std::endl;
+  std::cout << "mesh type: " << typeid(mesh_path).name() << std::endl;
+
+  typedef typename boost::graph_traits<TriangleMesh>::vertex_descriptor vertex_descriptor;
+
+  TriangleMesh tmesh;
+  const std::string filename = CGAL::data_file_path(mesh_path);
+
+  if(!CGAL::IO::read_polygon_mesh(filename, tmesh) || faces(tmesh).size() == 0)
+  {
+    std::cerr << "Invalid input file." << std::endl;
+    std::exit(1);
+  }
+
+  typename boost::property_map<TriangleMesh, CGAL::dynamic_vertex_property_t<FT>>::type
+    mean_curvature_map = get(CGAL::dynamic_vertex_property_t<FT>(), tmesh),
+    gaussian_curvature_map = get(CGAL::dynamic_vertex_property_t<FT>(), tmesh);
+
+  PMP::discrete_mean_curvatures(tmesh, mean_curvature_map);
+  PMP::discrete_Gaussian_curvatures(tmesh, gaussian_curvature_map);
+
+  FT mean_curvature_avg = 0, gaussian_curvature_avg = 0;
+  for(vertex_descriptor v : vertices(tmesh))
+  {
+    mean_curvature_avg += get(mean_curvature_map, v);
+    gaussian_curvature_avg += get(gaussian_curvature_map, v);
+  }
+
+  mean_curvature_avg /= vertices(tmesh).size();
+  gaussian_curvature_avg /= vertices(tmesh).size();
+
+  std::cout << "checking mean curvature..." << std::endl;
+  assert(passes_comparison(mean_curvature_avg, test_info.mean_curvature_avg, test_info.tolerance));
+
+  std::cout << "checking Gaussian curvature..." << std::endl;
+  assert(passes_comparison(gaussian_curvature_avg, test_info.gaussian_curvature_avg, test_info.tolerance));
+}
+
+template <typename PolygonMesh>
+void test_angle_sums(const std::string mesh_path,
+                     const std::vector<FT>& expected_values)
+{
+  typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
+
+  PolygonMesh pmesh;
+  const std::string filename = CGAL::data_file_path(mesh_path);
+
+  if(!CGAL::IO::read_polygon_mesh(filename, pmesh) || faces(pmesh).size() == 0)
+  {
+    std::cerr << "Invalid input file." << std::endl;
+    std::exit(1);
+  }
+
+  std::size_t pos = 0;
+  for(vertex_descriptor v : vertices(pmesh))
+  {
+    FT angle_sum = PMP::angle_sum(v, pmesh,
+                                  CGAL::parameters::geom_traits(K())
+                                                   .vertex_point_map(get(CGAL::vertex_point, pmesh)));
+    assert(passes_comparison(angle_sum, expected_values[pos++], 0.9));
+  }
+}
+
+int main(int, char**)
+{
+  // testing on a simple sphere(r = 0.5), on both Polyhedron & SurfaceMesh:
+  // Expected: Mean Curvature = 2, Gaussian Curvature = 4
+  test_curvatures<Polyhedron>("meshes/sphere.off", Average_test_info(2, 4));
+  test_curvatures<SMesh>("meshes/sphere.off", Average_test_info(2, 4));
+
+  // testing on a simple sphere(r = 10), on both Polyhedron & SurfaceMesh:
+  // Expected: Mean Curvature = 0.1, Gaussian Curvature = 0.01
+  test_curvatures<Polyhedron>("meshes/sphere966.off", Average_test_info(0.1, 0.01));
+  test_curvatures<SMesh>("meshes/sphere966.off", Average_test_info(0.1, 0.01));
+
+  // testing on a simple half cylinder(r = 1), on both Polyhedron & SurfaceMesh:
+  // Expected: Mean Curvature = 0.5, Gaussian Curvature = 0
+  // To be tested once the discrete curvatures are well defined for boundary vertices
+  // test_curvatures<Polyhedron>("meshes/cylinder.off", Average_test_info(0.5, 0));
+  // test_curvatures<SMesh>("meshes/cylinder.off", Average_test_info(0.5, 0));
+
+  test_angle_sums<Polyhedron>("meshes/quad.off", std::vector<FT>(4, 90));
+  test_angle_sums<SMesh>("meshes/quad.off", std::vector<FT>(4, 90));
+
+  test_angle_sums<Polyhedron>("meshes/regular_tetrahedron.off", std::vector<FT>(4, 180));
+  test_angle_sums<SMesh>("meshes/regular_tetrahedron.off", std::vector<FT>(4, 180));
+
+  test_angle_sums<Polyhedron>("meshes/cube_quad.off", std::vector<FT>(8, 270));
+  test_angle_sums<SMesh>("meshes/cube_quad.off", std::vector<FT>(8, 270));
+
+  test_angle_sums<Polyhedron>("meshes/cube_poly.off", std::vector<FT>(8, 270));
+  test_angle_sums<SMesh>("meshes/cube_poly.off", std::vector<FT>(8, 270));
+
+  std::cout << "Done." << std::endl;
+  return EXIT_SUCCESS;
+}

Polygon_mesh_processing/test/Polygon_mesh_processing/test_interpolated_corrected_curvatures.cpp

@@ -9,8 +9,9 @@
 
 #include <boost/graph/graph_traits.hpp>
 
+#include <functional>
 #include <iostream>
-#include <unordered_map>
+#include <string>
 
 #define ABS_ERROR 1e-6
 
@@ -181,7 +182,7 @@ void test_average_curvatures(std::string mesh_path,
 int main()
 {
   // testing on a simple sphere(r = 0.5), on both Polyhedron & SurfaceMesh:
-  // For this mesh, ina addition to the whole mesh functions, we also compare against the single vertex
+  // For this mesh, in addition to the whole mesh functions, we also compare against the single vertex
   // curvature functions to make sure the produce the same results
   // Expected: Mean Curvature = 2, Gaussian Curvature = 4, Principal Curvatures = 2 & 2 so 2 on avg.
   test_average_curvatures<Polyhedron>("meshes/sphere.off", Average_test_info(2, 4, 2), true);

Weights/include/CGAL/Weights/cotangent_weights.h

@@ -344,7 +344,6 @@ public:
     return cotangent_weight_calculator(he);
   }
 
-private:
   FT voronoi(const vertex_descriptor v0) const
   {
     auto squared_length_3 = m_traits.compute_squared_length_3_object();
@@ -354,11 +353,12 @@ private:
     for (const halfedge_descriptor he : halfedges_around_target(halfedge(v0, m_pmesh), m_pmesh))
     {
       CGAL_assertion(v0 == target(he, m_pmesh));
-      CGAL_assertion(CGAL::is_triangle(he, m_pmesh));
 
       if (is_border(he, m_pmesh))
         continue;
 
+      CGAL_assertion(CGAL::is_triangle(he, m_pmesh));
+
       const vertex_descriptor v1 = source(he, m_pmesh);
       const vertex_descriptor v2 = target(next(he, m_pmesh), m_pmesh);
 

Weights/include/CGAL/Weights/internal/utils.h

@@ -44,7 +44,7 @@ private:
 public:
   FT operator()(const FT value) const
   {
-    return static_cast<FT>(CGAL::sqrt(CGAL::to_double(CGAL::abs(value))));
+    return CGAL::approximate_sqrt(CGAL::abs(value));
   }
 };