diff --git a/Installation/CHANGES.md b/Installation/CHANGES.md index d977d162f6b..bdae9dd56fb 100644 --- a/Installation/CHANGES.md +++ b/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. diff --git a/Kernel_23/include/CGAL/Kernel/function_objects.h b/Kernel_23/include/CGAL/Kernel/function_objects.h index 5a7b0774e6d..c49094aefb2 100644 --- a/Kernel_23/include/CGAL/Kernel/function_objects.h +++ b/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; - - if(cosine > 1.){ - cosine = 1.; - } - if(cosine < -1.){ - cosine = -1.; - } - + double cosine = std::clamp(dot / product, -1., 1.); return std::acos(cosine) * 180./CGAL_PI; } diff --git a/Lab/demo/Lab/Color_map.h b/Lab/demo/Lab/Color_map.h index a7a608112eb..909fc2f6a93 100644 --- a/Lab/demo/Lab/Color_map.h +++ b/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(1)) / nb_of_colors; + const qreal step = (static_cast(0.85)) / nb_of_colors; + qreal hue = base_color.hueF(); qreal h = (hue == -1) ? 0 : hue; for(std::size_t i=0; i #include #include +#include #include #include @@ -350,11 +351,11 @@ private: template void displayMapLegend(const std::vector& 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(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); - painter.drawText(text_rect, Qt::AlignCenter, QObject::tr("%1").arg(values[i], 0, 'f', 3, QLatin1Char(' '))); + 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', 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" || - property_name == "Interpolated Corrected Gaussian Curvature"); - dock_widget->expandingRadiusLabel->setVisible(is_curvature_property); - dock_widget->expandingRadiusSlider->setVisible(is_curvature_property); - dock_widget->expandingRadiusLabel->setEnabled(is_curvature_property); - dock_widget->expandingRadiusSlider->setEnabled(is_curvature_property); + const bool is_interpolated_curvature_property = (property_name == "Interpolated Corrected Mean Curvature" || + property_name == "Interpolated Corrected Gaussian Curvature"); + dock_widget->expandingRadiusLabel->setVisible(is_interpolated_curvature_property); + dock_widget->expandingRadiusSlider->setVisible(is_interpolated_curvature_property); + dock_widget->expandingRadiusLabel->setEnabled(is_interpolated_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("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 vdc; + std::tie(vdc, not_initialized) = sm->add_property_map(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(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; diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt b/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt index 6b74ce9696a..50f198e8317 100644 --- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt +++ b/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()` diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt b/Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt index 08d8d798f41..c4bf5a424fc 100644 --- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt +++ b/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 diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/curvature.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/curvature.h new file mode 100644 index 00000000000..f612535dfc5 --- /dev/null +++ b/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 + +#include +#include +#include +#include + +#include +#include + +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::%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 +#ifdef DOXYGEN_RUNNING +FT +#else +typename GetGeomTraits::type::FT +#endif +angle_sum(typename boost::graph_traits::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::type; + using FT = typename Geom_traits::FT; + + typename GetVertexPointMap::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(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 discrete Gaussian curvature. + * + * @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::%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 +#ifdef DOXYGEN_RUNNING +FT +#else +typename GetGeomTraits::type::FT +#endif +discrete_Gaussian_curvature(typename boost::graph_traits::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::type; + using FT = typename GeomTraits::FT; + using Vector_3 = typename GeomTraits::Vector_3; + + using VertexPointMap = typename GetVertexPointMap::const_type; + using halfedge_descriptor = typename boost::graph_traits::halfedge_descriptor; + + GeomTraits gt = choose_parameter(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 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 discrete Gaussian curvature. + * + * @tparam TriangleMesh a model of `FaceGraph` + * @tparam VertexCurvatureMap must be a model of `WritablePropertyMap` with key type + * `boost::graph_traits::%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::%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 +void discrete_Gaussian_curvatures(const TriangleMesh& tmesh, + VertexCurvatureMap vcm, + const CGAL_NP_CLASS& np = parameters::default_values()) +{ + using vertex_descriptor = typename boost::graph_traits::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 discrete mean curvature. + * + * @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::%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 +#ifdef DOXYGEN_RUNNING +FT +#else +typename GetGeomTraits::type::FT +#endif +discrete_mean_curvature(typename boost::graph_traits::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::type; + using FT = typename GeomTraits::FT; + using Vector_3 = typename GeomTraits::Vector_3; + + using VertexPointMap = typename GetVertexPointMap::const_type; + using Point_ref = typename boost::property_traits::reference; + + using vertex_descriptor = typename boost::graph_traits::vertex_descriptor; + using halfedge_descriptor = typename boost::graph_traits::halfedge_descriptor; + + GeomTraits gt = choose_parameter(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 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 discrete mean curvature. + * + * @tparam TriangleMesh a model of `FaceGraph` + * @tparam VertexCurvatureMap must be a model of `WritablePropertyMap` with key type + * `boost::graph_traits::%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::%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 +void discrete_mean_curvatures(const TriangleMesh& tmesh, + VertexCurvatureMap vcm, + const CGAL_NP_CLASS& np = parameters::default_values()) +{ + using vertex_descriptor = typename boost::graph_traits::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 diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt b/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt index 12cb26ca381..1e442ae768f 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt +++ b/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") diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_discrete_curvatures.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_discrete_curvatures.cpp new file mode 100644 index 00000000000..d37bf7471ad --- /dev/null +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_discrete_curvatures.cpp @@ -0,0 +1,146 @@ +#include +#include +#include + +#include + +#include + +#include +#include + +#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 SMesh; +typedef CGAL::Polyhedron_3 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 +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::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>::type + mean_curvature_map = get(CGAL::dynamic_vertex_property_t(), tmesh), + gaussian_curvature_map = get(CGAL::dynamic_vertex_property_t(), 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 +void test_angle_sums(const std::string mesh_path, + const std::vector& expected_values) +{ + typedef typename boost::graph_traits::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("meshes/sphere.off", Average_test_info(2, 4)); + test_curvatures("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("meshes/sphere966.off", Average_test_info(0.1, 0.01)); + test_curvatures("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("meshes/cylinder.off", Average_test_info(0.5, 0)); + // test_curvatures("meshes/cylinder.off", Average_test_info(0.5, 0)); + + test_angle_sums("meshes/quad.off", std::vector(4, 90)); + test_angle_sums("meshes/quad.off", std::vector(4, 90)); + + test_angle_sums("meshes/regular_tetrahedron.off", std::vector(4, 180)); + test_angle_sums("meshes/regular_tetrahedron.off", std::vector(4, 180)); + + test_angle_sums("meshes/cube_quad.off", std::vector(8, 270)); + test_angle_sums("meshes/cube_quad.off", std::vector(8, 270)); + + test_angle_sums("meshes/cube_poly.off", std::vector(8, 270)); + test_angle_sums("meshes/cube_poly.off", std::vector(8, 270)); + + std::cout << "Done." << std::endl; + return EXIT_SUCCESS; +} diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_interpolated_corrected_curvatures.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_interpolated_corrected_curvatures.cpp index 412443a42c7..a15258b0fc9 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_interpolated_corrected_curvatures.cpp +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_interpolated_corrected_curvatures.cpp @@ -9,8 +9,9 @@ #include +#include #include -#include +#include #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("meshes/sphere.off", Average_test_info(2, 4, 2), true); diff --git a/Weights/include/CGAL/Weights/cotangent_weights.h b/Weights/include/CGAL/Weights/cotangent_weights.h index 4f389e3b06b..34a6d79f178 100644 --- a/Weights/include/CGAL/Weights/cotangent_weights.h +++ b/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); diff --git a/Weights/include/CGAL/Weights/internal/utils.h b/Weights/include/CGAL/Weights/internal/utils.h index 850fccff1e8..80822b68c84 100644 --- a/Weights/include/CGAL/Weights/internal/utils.h +++ b/Weights/include/CGAL/Weights/internal/utils.h @@ -44,7 +44,7 @@ private: public: FT operator()(const FT value) const { - return static_cast(CGAL::sqrt(CGAL::to_double(CGAL::abs(value)))); + return CGAL::approximate_sqrt(CGAL::abs(value)); } };