From ef8cd2930d28644f1f2e599e1dd544ccda9325fc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Thu, 28 Oct 2021 14:18:27 +0200 Subject: [PATCH 01/27] add functions to compute discrete curvature --- .../internal/curvature.h | 184 ++++++++++++++++++ 1 file changed, 184 insertions(+) create mode 100644 Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h new file mode 100644 index 00000000000..aace605e6ca --- /dev/null +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h @@ -0,0 +1,184 @@ +// 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) : Mael Rouxel-Labbé + +#ifndef CGAL_PMP_INTERNAL_CURVATURE_H +#define CGAL_PMP_INTERNAL_CURVATURE_H + +#include + +#include +#include +#include + +namespace CGAL { +namespace Polygon_mesh_processing { + +template +typename GetGeomTraits::type::FT +vertex_discrete_gaussian_curvature(typename boost::graph_traits::vertex_descriptor vd, + const TriangleMesh& tm, + const NamedParameters& np) +{ + typedef typename GetGeomTraits::type GeomTraits; + typedef typename GetVertexPointMap::const_type VertexPointMap; + typedef typename GeomTraits::FT FT; + typedef typename GeomTraits::Vector_3 Vector_3; + typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + + using parameters::choose_parameter; + using parameters::get_parameter; + + 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, tm)); + + typename GeomTraits::Construct_vector_3 vector = gt.construct_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(); + typename GeomTraits::Construct_cross_product_vector_3 cross_product = gt.construct_cross_product_vector_3_object(); + typename GeomTraits::Compute_approximate_angle_3 approximate_angle = gt.compute_approximate_angle_3_object(); + + const FT two_pi = 2 * CGAL_PI; + + FT ki = 0; + for(halfedge_descriptor h : CGAL::halfedges_around_target(vd, tm)) + { + if(is_border(h, tm)) + continue; + + const Vector_3 v0 = vector(get(vpm, vd), get(vpm, target(next(h, tm), tm))); // p1p2 + const Vector_3 v1 = vector(get(vpm, vd), get(vpm, source(h, tm))); // p1p0 + + const FT dot = scalar_product(v0, v1); + const Vector_3 cross = cross_product(v0, v1); + const FT sqcn = squared_length(cross); + + if(dot != FT(0) && sqcn != FT(0)) + { + const FT loc_ki = std::atan2(CGAL::approximate_sqrt(sqcn), dot); + + FT ia = approximate_angle(get(vpm, source(h, tm)), + get(vpm, target(h, tm)), + get(vpm, target(next(h, tm), tm))); + ia *= CGAL_PI / FT(180.); + + ki += loc_ki; + } + } + + return two_pi - ki; +} + +template +void discrete_gaussian_curvature(const TriangleMesh& tm, + VertexCurvatureMap vcm, + const NamedParameters& np) +{ + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + + for(vertex_descriptor v : vertices(tm)) + put(vcm, v, vertex_discrete_gaussian_curvature(v, tm, np)); +} + + +template +typename GetGeomTraits::type::FT +vertex_discrete_mean_curvature(typename boost::graph_traits::vertex_descriptor vd, + const TriangleMesh& tm, + const NamedParameters& np) +{ + typedef typename GetGeomTraits::type GeomTraits; + typedef typename GetVertexPointMap::const_type VertexPointMap; + typedef typename GeomTraits::FT FT; + typedef typename GeomTraits::Point_3 Point_3; + typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + + using parameters::choose_parameter; + using parameters::get_parameter; + + 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, tm)); + + typename GeomTraits::Compute_squared_distance_3 squared_distance = gt.compute_squared_distance_3_object(); + typename GeomTraits::Compute_approximate_dihedral_angle_3 approximate_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(vd, tm)) + { + const Point_3& p = get(vpm, source(h, tm)); + const Point_3& q = get(vpm, target(h, tm)); + const Point_3& r = get(vpm, target(next(h, tm), tm)); + const Point_3& s = get(vpm, target(next(opposite(h, tm), tm), tm)); + const FT l = squared_distance(p,q); + + FT phi = CGAL_PI * approximate_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; +} + +template +void discrete_mean_curvature(const TriangleMesh& tm, + VertexCurvatureMap vcm, + const NamedParameters& np) +{ + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + + for(vertex_descriptor v : vertices(tm)) + put(vcm, v, vertex_discrete_mean_curvature(v, tm, np)); +} + +/// convenience overloads + +template +auto +vertex_discrete_gaussian_curvature(typename boost::graph_traits::vertex_descriptor vd, + const TriangleMesh& tm) +{ + return vertex_discrete_gaussian_curvature(vd, tm, parameters::all_default()); +} + +template +void discrete_gaussian_curvature(const TriangleMesh& tm, + VertexCurvatureMap vcm) +{ + discrete_gaussian_curvature(tm, vcm, parameters::all_default()); +} + +template +auto +vertex_discrete_mean_curvature(typename boost::graph_traits::vertex_descriptor vd, + const TriangleMesh& tm) +{ + return vertex_discrete_mean_curvature(vd, tm, parameters::all_default()); +} + +template +void discrete_mean_curvature(const TriangleMesh& tm, + VertexCurvatureMap vcm) +{ + discrete_mean_curvature(tm, vcm, parameters::all_default()); +} + +} } // CGAL::Polygon_mesh_processing::experimental + +#endif //CGAL_PMP_INTERNAL_CURVATURE_H From da3179abe0b94160d7a4478ba4365c8f4e5d4b0b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Thu, 28 Oct 2021 15:56:52 +0200 Subject: [PATCH 02/27] remove debug variable --- .../CGAL/Polygon_mesh_processing/internal/curvature.h | 7 ------- 1 file changed, 7 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h index aace605e6ca..6b2ffbb01cc 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h @@ -45,7 +45,6 @@ vertex_discrete_gaussian_curvature(typename boost::graph_traits::v 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(); typename GeomTraits::Construct_cross_product_vector_3 cross_product = gt.construct_cross_product_vector_3_object(); - typename GeomTraits::Compute_approximate_angle_3 approximate_angle = gt.compute_approximate_angle_3_object(); const FT two_pi = 2 * CGAL_PI; @@ -65,12 +64,6 @@ vertex_discrete_gaussian_curvature(typename boost::graph_traits::v if(dot != FT(0) && sqcn != FT(0)) { const FT loc_ki = std::atan2(CGAL::approximate_sqrt(sqcn), dot); - - FT ia = approximate_angle(get(vpm, source(h, tm)), - get(vpm, target(h, tm)), - get(vpm, target(next(h, tm), tm))); - ia *= CGAL_PI / FT(180.); - ki += loc_ki; } } From 2b7246a14d02a67177c3b43a6557b6fc3104b2b1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Thu, 28 Oct 2021 16:58:19 +0200 Subject: [PATCH 03/27] handle border cases --- .../Polygon_mesh_processing/internal/curvature.h | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h index 6b2ffbb01cc..7b976fee3a6 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h @@ -60,11 +60,18 @@ vertex_discrete_gaussian_curvature(typename boost::graph_traits::v const FT dot = scalar_product(v0, v1); const Vector_3 cross = cross_product(v0, v1); const FT sqcn = squared_length(cross); - - if(dot != FT(0) && sqcn != FT(0)) + if(dot == FT(0)) + ki += CGAL_PI/FT(2); + else { - const FT loc_ki = std::atan2(CGAL::approximate_sqrt(sqcn), dot); - ki += loc_ki; + if (sqcn == FT(0)) + { + if (dot < 0) + ki += CGAL_PI; + } + else{ + ki += std::atan2(CGAL::approximate_sqrt(sqcn), dot); + } } } From 8a52d7c65cba07edad57c616ecc87e5b34ca8a1e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Fri, 29 Oct 2021 11:44:18 +0200 Subject: [PATCH 04/27] remove unused function --- .../Plugins/Display/Display_property_plugin.cpp | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/Polyhedron/demo/Polyhedron/Plugins/Display/Display_property_plugin.cpp b/Polyhedron/demo/Polyhedron/Plugins/Display/Display_property_plugin.cpp index c787c901e5b..a767ebf31c3 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/Display/Display_property_plugin.cpp +++ b/Polyhedron/demo/Polyhedron/Plugins/Display/Display_property_plugin.cpp @@ -859,18 +859,6 @@ private Q_SLOTS: treat_sm_property("f:angle", item->face_graph()); } - bool resetAngles(Scene_surface_mesh_item* item) - { - SMesh& smesh = *item->face_graph(); - if(!smesh.property_map("f:angle").second) - { - return false; - } - dock_widget->minBox->setValue(angles_min[item].first); - dock_widget->maxBox->setValue(angles_max[item].first); - return true; - } - // AF: This function gets called when we click on the button "Colorize" bool displayHeatIntensity(Scene_surface_mesh_item* item, bool iDT = false) { From 80ff6bfb67bf73fe3f85f6291a2099826f47613c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Fri, 29 Oct 2021 11:52:25 +0200 Subject: [PATCH 05/27] add display for discrete mean and gaussian curvature --- .../Display/Display_property_plugin.cpp | 190 ++++++++++++++++-- 1 file changed, 177 insertions(+), 13 deletions(-) diff --git a/Polyhedron/demo/Polyhedron/Plugins/Display/Display_property_plugin.cpp b/Polyhedron/demo/Polyhedron/Plugins/Display/Display_property_plugin.cpp index a767ebf31c3..1b00f6c074a 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/Display/Display_property_plugin.cpp +++ b/Polyhedron/demo/Polyhedron/Plugins/Display/Display_property_plugin.cpp @@ -15,6 +15,7 @@ #include #include #include +#include #include "Scene_points_with_normal_item.h" @@ -526,10 +527,12 @@ private Q_SLOTS: if(sm_item) { - dock_widget->propertyBox->addItem("Smallest Angle Per Face"); - dock_widget->propertyBox->addItem("Scaled Jacobian"); - dock_widget->propertyBox->addItem("Heat Intensity"); - dock_widget->propertyBox->addItem("Heat Intensity (Intrinsic Delaunay)"); + dock_widget->propertyBox->addItem("Smallest Angle Per Face"); + dock_widget->propertyBox->addItem("Scaled Jacobian"); + dock_widget->propertyBox->addItem("Heat Intensity"); + dock_widget->propertyBox->addItem("Heat Intensity (Intrinsic Delaunay)"); + dock_widget->propertyBox->addItem("Discrete Mean Curvature"); + dock_widget->propertyBox->addItem("Discrete Gaussian Curvature"); detectSMScalarProperties(sm_item->face_graph()); } @@ -604,6 +607,14 @@ private Q_SLOTS: return; sm_item->setRenderingMode(Gouraud); break; + case 4: // Discrete Mean Curvature + displayDiscreteMeanCurvature(sm_item); + sm_item->setRenderingMode(Gouraud); + break; + case 5: // Discrete Gaussian Curvature + displayDiscreteGaussianCurvature(sm_item); + sm_item->setRenderingMode(Gouraud); + break; default: if(dock_widget->propertyBox->currentText().contains("v:")) { @@ -638,6 +649,14 @@ private Q_SLOTS: sm_item->face_graph()->property_map("f:angle"); if(does_exist) sm_item->face_graph()->remove_property_map(pmap); + std::tie(pmap, does_exist) = + sm_item->face_graph()->property_map("v:discrete_mean_curvature"); + if(does_exist) + sm_item->face_graph()->remove_property_map(pmap); + std::tie(pmap, does_exist) = + sm_item->face_graph()->property_map("v:discrete_gaussian_curvature"); + if(does_exist) + sm_item->face_graph()->remove_property_map(pmap); }); QApplication::restoreOverrideCursor(); sm_item->invalidateOpenGLBuffers(); @@ -669,13 +688,21 @@ private Q_SLOTS: switch(dock_widget->propertyBox->currentIndex()) { case 0: - dock_widget->zoomToMinButton->setEnabled(angles_max.count(sm_item)>0 ); + dock_widget->zoomToMinButton->setEnabled(angles_min.count(sm_item)>0 ); dock_widget->zoomToMaxButton->setEnabled(angles_max.count(sm_item)>0 ); break; case 1: - dock_widget->zoomToMinButton->setEnabled(jacobian_max.count(sm_item)>0); + dock_widget->zoomToMinButton->setEnabled(jacobian_min.count(sm_item)>0); dock_widget->zoomToMaxButton->setEnabled(jacobian_max.count(sm_item)>0); break; + case 4: + dock_widget->zoomToMinButton->setEnabled(dmc_min.count(sm_item)>0); + dock_widget->zoomToMaxButton->setEnabled(dmc_max.count(sm_item)>0); + break; + case 5: + dock_widget->zoomToMinButton->setEnabled(dgc_min.count(sm_item)>0); + dock_widget->zoomToMaxButton->setEnabled(dgc_max.count(sm_item)>0); + break; default: break; } @@ -702,11 +729,22 @@ private Q_SLOTS: { smesh.remove_property_map(angles); } + SMesh::Property_map dmc; + std::tie(dmc, found) = smesh.property_map("v:discrete_mean_curvature"); + if(found) + { + smesh.remove_property_map(dmc); + } + SMesh::Property_map dgc; + std::tie(dgc, found) = smesh.property_map("v:discrete_gaussian_curvature"); + if(found) + { + smesh.remove_property_map(dgc); + } } void displayScaledJacobian(Scene_surface_mesh_item* item) { - SMesh& smesh = *item->face_graph(); //compute and store the jacobian per face bool non_init; @@ -743,18 +781,77 @@ private Q_SLOTS: treat_sm_property("f:jacobian", item->face_graph()); } - bool resetScaledJacobian(Scene_surface_mesh_item* item) + void displayDiscreteMeanCurvature(Scene_surface_mesh_item* item) { SMesh& smesh = *item->face_graph(); - if(!smesh.property_map("f:jacobian").second) + //compute once and store the value per vertex + bool non_init; + SMesh::Property_map mcurvature; + std::tie(mcurvature, non_init) = smesh.add_property_map("v:discrete_mean_curvature", 0); + if(non_init) { - return false; + CGAL::Polygon_mesh_processing::discrete_mean_curvature(smesh, mcurvature); + + double res_min = ARBITRARY_DBL_MAX, + res_max = -ARBITRARY_DBL_MAX; + SMesh::Vertex_index min_index, max_index; + for (SMesh::Vertex_index v : vertices(smesh)) + { + if (mcurvature[v] > res_max) + { + res_max = mcurvature[v]; + max_index = v; + } + if (mcurvature[v] < res_min) + { + res_min = mcurvature[v]; + min_index = v; + } + } + dmc_max[item]=std::make_pair(res_max, max_index); + dmc_min[item]=std::make_pair(res_min, min_index); + + connect(item, &Scene_surface_mesh_item::itemChanged, + this, &DisplayPropertyPlugin::resetProperty); } - dock_widget->minBox->setValue(jacobian_min[item].first-0.01); - dock_widget->maxBox->setValue(jacobian_max[item].first); - return true; + treat_sm_property("v:discrete_mean_curvature", item->face_graph()); } + void displayDiscreteGaussianCurvature(Scene_surface_mesh_item* item) + { + SMesh& smesh = *item->face_graph(); + //compute once and store the value per vertex + bool non_init; + SMesh::Property_map mcurvature; + std::tie(mcurvature, non_init) = smesh.add_property_map("v:discrete_gaussian_curvature", 0); + if(non_init) + { + CGAL::Polygon_mesh_processing::discrete_gaussian_curvature(smesh, mcurvature); + + double res_min = ARBITRARY_DBL_MAX, + res_max = -ARBITRARY_DBL_MAX; + SMesh::Vertex_index min_index, max_index; + for (SMesh::Vertex_index v : vertices(smesh)) + { + if (mcurvature[v] > res_max) + { + res_max = mcurvature[v]; + max_index = v; + } + if (mcurvature[v] < res_min) + { + res_min = mcurvature[v]; + min_index = v; + } + } + dgc_max[item]=std::make_pair(res_max, max_index); + dgc_min[item]=std::make_pair(res_min, min_index); + + connect(item, &Scene_surface_mesh_item::itemChanged, + this, &DisplayPropertyPlugin::resetProperty); + } + treat_sm_property("v:discrete_gaussian_curvature", item->face_graph()); + } void displayAngles(Scene_surface_mesh_item* item) { @@ -1058,6 +1155,30 @@ private Q_SLOTS: case 3: dock_widget->sourcePointsButton->setEnabled(true); CGAL_FALLTHROUGH; + case 4: + dock_widget->groupBox-> setEnabled(true); + dock_widget->groupBox_3->setEnabled(true); + + dock_widget->minBox->setMinimum(-1000); + dock_widget->minBox->setMaximum(1000); + dock_widget->minBox->setValue(0); + + dock_widget->maxBox->setMinimum(-1000); + dock_widget->maxBox->setMaximum(1000); + dock_widget->maxBox->setValue(2); + break; + case 5: + dock_widget->groupBox-> setEnabled(true); + dock_widget->groupBox_3->setEnabled(true); + + dock_widget->minBox->setMinimum(-1000); + dock_widget->minBox->setMaximum(1000); + dock_widget->minBox->setValue(0); + + dock_widget->maxBox->setMinimum(-1000); + dock_widget->maxBox->setMaximum(1000); + dock_widget->maxBox->setValue(2); + break; default: dock_widget->maxBox->setMinimum(-99999999); dock_widget->maxBox->setMaximum(99999999); @@ -1103,6 +1224,24 @@ private Q_SLOTS: dummy_p); } break; + case 4: + { + ::zoomToId(*item->face_graph(), + QString("v%1").arg(dmc_min[item].second), + getActiveViewer(), + dummy_fd, + dummy_p); + } + break; + case 5: + { + ::zoomToId(*item->face_graph(), + QString("v%1").arg(dgc_min[item].second), + getActiveViewer(), + dummy_fd, + dummy_p); + } + break; default: break; } @@ -1136,6 +1275,24 @@ private Q_SLOTS: dummy_p); } break; + case 4: + { + ::zoomToId(*item->face_graph(), + QString("v%1").arg(dmc_max[item].second), + getActiveViewer(), + dummy_fd, + dummy_p); + } + break; + case 5: + { + ::zoomToId(*item->face_graph(), + QString("v%1").arg(dgc_max[item].second), + getActiveViewer(), + dummy_fd, + dummy_p); + } + break; default: break; } @@ -1425,6 +1582,13 @@ private: boost::unordered_map > angles_min; boost::unordered_map > angles_max; + + boost::unordered_map > dmc_min; + boost::unordered_map > dmc_max; + + boost::unordered_map > dgc_min; + boost::unordered_map > dgc_max; + boost::unordered_map is_source; From 0fe940555dbd08d4ddb332c0612c2e2cded6216e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Sat, 7 Dec 2024 23:41:38 +0100 Subject: [PATCH 06/27] Fix namespace --- .../CGAL/Polygon_mesh_processing/internal/curvature.h | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h index 7b976fee3a6..aad4304c3fc 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h @@ -21,6 +21,7 @@ namespace CGAL { namespace Polygon_mesh_processing { +namespace exprimental { template typename GetGeomTraits::type::FT @@ -179,6 +180,8 @@ void discrete_mean_curvature(const TriangleMesh& tm, discrete_mean_curvature(tm, vcm, parameters::all_default()); } -} } // CGAL::Polygon_mesh_processing::experimental +} // namespace experimental +} // namespace Polygon_mesh_processing +} // namespace CGAL #endif //CGAL_PMP_INTERNAL_CURVATURE_H From 693fd251747f050ce328fd2010be52ef7760b885 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Sat, 7 Dec 2024 23:42:31 +0100 Subject: [PATCH 07/27] Use a capital 'G' --- .../Polygon_mesh_processing/internal/curvature.h | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h index aad4304c3fc..67f0a79b559 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h @@ -25,7 +25,7 @@ namespace exprimental { template typename GetGeomTraits::type::FT -vertex_discrete_gaussian_curvature(typename boost::graph_traits::vertex_descriptor vd, +vertex_discrete_Gaussian_curvature(typename boost::graph_traits::vertex_descriptor vd, const TriangleMesh& tm, const NamedParameters& np) { @@ -80,14 +80,14 @@ vertex_discrete_gaussian_curvature(typename boost::graph_traits::v } template -void discrete_gaussian_curvature(const TriangleMesh& tm, +void discrete_Gaussian_curvature(const TriangleMesh& tm, VertexCurvatureMap vcm, const NamedParameters& np) { typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; for(vertex_descriptor v : vertices(tm)) - put(vcm, v, vertex_discrete_gaussian_curvature(v, tm, np)); + put(vcm, v, vertex_discrete_Gaussian_curvature(v, tm, np)); } @@ -152,17 +152,17 @@ void discrete_mean_curvature(const TriangleMesh& tm, template auto -vertex_discrete_gaussian_curvature(typename boost::graph_traits::vertex_descriptor vd, +vertex_discrete_Gaussian_curvature(typename boost::graph_traits::vertex_descriptor vd, const TriangleMesh& tm) { - return vertex_discrete_gaussian_curvature(vd, tm, parameters::all_default()); + return vertex_discrete_Gaussian_curvature(vd, tm, parameters::all_default()); } template -void discrete_gaussian_curvature(const TriangleMesh& tm, +void discrete_Gaussian_curvature(const TriangleMesh& tm, VertexCurvatureMap vcm) { - discrete_gaussian_curvature(tm, vcm, parameters::all_default()); + discrete_Gaussian_curvature(tm, vcm, parameters::all_default()); } template From f2299f844f641b44c810d740c9530f74a413540d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Sat, 7 Dec 2024 23:42:49 +0100 Subject: [PATCH 08/27] Use plural for the all-vertices function --- .../CGAL/Polygon_mesh_processing/internal/curvature.h | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h index 67f0a79b559..1cae4dd1ab9 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h @@ -80,9 +80,9 @@ vertex_discrete_Gaussian_curvature(typename boost::graph_traits::v } template -void discrete_Gaussian_curvature(const TriangleMesh& tm, - VertexCurvatureMap vcm, - const NamedParameters& np) +void discrete_Gaussian_curvatures(const TriangleMesh& tm, + VertexCurvatureMap vcm, + const NamedParameters& np) { typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; @@ -159,10 +159,10 @@ vertex_discrete_Gaussian_curvature(typename boost::graph_traits::v } template -void discrete_Gaussian_curvature(const TriangleMesh& tm, +void discrete_Gaussian_curvatures(const TriangleMesh& tm, VertexCurvatureMap vcm) { - discrete_Gaussian_curvature(tm, vcm, parameters::all_default()); + discrete_Gaussian_curvatures(tm, vcm, parameters::all_default()); } template From a0e319e67b1d41e8c7784c9ecfe9aed43927df63 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Sat, 7 Dec 2024 23:43:07 +0100 Subject: [PATCH 09/27] Re-attach convenience overloads to their respective main function --- .../internal/curvature.h | 31 +++++++++---------- 1 file changed, 14 insertions(+), 17 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h index 1cae4dd1ab9..68ccebaa75e 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h @@ -79,6 +79,14 @@ vertex_discrete_Gaussian_curvature(typename boost::graph_traits::v return two_pi - ki; } +template +auto +vertex_discrete_Gaussian_curvature(typename boost::graph_traits::vertex_descriptor vd, + const TriangleMesh& tm) +{ + return vertex_discrete_Gaussian_curvature(vd, tm, parameters::all_default()); +} + template void discrete_Gaussian_curvatures(const TriangleMesh& tm, VertexCurvatureMap vcm, @@ -90,6 +98,12 @@ void discrete_Gaussian_curvatures(const TriangleMesh& tm, put(vcm, v, vertex_discrete_Gaussian_curvature(v, tm, np)); } +template +void discrete_Gaussian_curvatures(const TriangleMesh& tm, + VertexCurvatureMap vcm) +{ + discrete_Gaussian_curvatures(tm, vcm, parameters::all_default()); +} template typename GetGeomTraits::type::FT @@ -148,23 +162,6 @@ void discrete_mean_curvature(const TriangleMesh& tm, put(vcm, v, vertex_discrete_mean_curvature(v, tm, np)); } -/// convenience overloads - -template -auto -vertex_discrete_Gaussian_curvature(typename boost::graph_traits::vertex_descriptor vd, - const TriangleMesh& tm) -{ - return vertex_discrete_Gaussian_curvature(vd, tm, parameters::all_default()); -} - -template -void discrete_Gaussian_curvatures(const TriangleMesh& tm, - VertexCurvatureMap vcm) -{ - discrete_Gaussian_curvatures(tm, vcm, parameters::all_default()); -} - template auto vertex_discrete_mean_curvature(typename boost::graph_traits::vertex_descriptor vd, From 25cd25c477cc234337b56ba2ede0894a90a01300 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Sat, 7 Dec 2024 23:44:03 +0100 Subject: [PATCH 10/27] Tiny indentation fix --- .../internal/curvature.h | 41 ++++++++++--------- 1 file changed, 22 insertions(+), 19 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h index 68ccebaa75e..9eca4794b80 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h @@ -25,9 +25,9 @@ namespace exprimental { template typename GetGeomTraits::type::FT -vertex_discrete_Gaussian_curvature(typename boost::graph_traits::vertex_descriptor vd, - const TriangleMesh& tm, - const NamedParameters& np) +vertex_discrete_Gaussian_curvature(typename boost::graph_traits::vertex_descriptor v, + const TriangleMesh& tm, + const NamedParameters& np) { typedef typename GetGeomTraits::type GeomTraits; typedef typename GetVertexPointMap::const_type VertexPointMap; @@ -50,27 +50,30 @@ vertex_discrete_Gaussian_curvature(typename boost::graph_traits::v const FT two_pi = 2 * CGAL_PI; FT ki = 0; - for(halfedge_descriptor h : CGAL::halfedges_around_target(vd, tm)) + for(halfedge_descriptor h : CGAL::halfedges_around_target(v, tm)) { if(is_border(h, tm)) continue; - const Vector_3 v0 = vector(get(vpm, vd), get(vpm, target(next(h, tm), tm))); // p1p2 - const Vector_3 v1 = vector(get(vpm, vd), get(vpm, source(h, tm))); // p1p0 + const Vector_3 v0 = vector(get(vpm, v), get(vpm, target(next(h, tm), tm))); // p1p2 + const Vector_3 v1 = vector(get(vpm, v), get(vpm, source(h, tm))); // p1p0 const FT dot = scalar_product(v0, v1); const Vector_3 cross = cross_product(v0, v1); const FT sqcn = squared_length(cross); if(dot == FT(0)) + { ki += CGAL_PI/FT(2); + } else { - if (sqcn == FT(0)) + if(sqcn == FT(0)) { - if (dot < 0) + if(dot < 0) ki += CGAL_PI; } - else{ + else + { ki += std::atan2(CGAL::approximate_sqrt(sqcn), dot); } } @@ -81,10 +84,10 @@ vertex_discrete_Gaussian_curvature(typename boost::graph_traits::v template auto -vertex_discrete_Gaussian_curvature(typename boost::graph_traits::vertex_descriptor vd, +vertex_discrete_Gaussian_curvature(typename boost::graph_traits::vertex_descriptor v, const TriangleMesh& tm) { - return vertex_discrete_Gaussian_curvature(vd, tm, parameters::all_default()); + return vertex_discrete_Gaussian_curvature(v, tm, parameters::all_default()); } template @@ -107,9 +110,9 @@ void discrete_Gaussian_curvatures(const TriangleMesh& tm, template typename GetGeomTraits::type::FT -vertex_discrete_mean_curvature(typename boost::graph_traits::vertex_descriptor vd, - const TriangleMesh& tm, - const NamedParameters& np) +vertex_discrete_mean_curvature(typename boost::graph_traits::vertex_descriptor v, + const TriangleMesh& tm, + const NamedParameters& np) { typedef typename GetGeomTraits::type GeomTraits; typedef typename GetVertexPointMap::const_type VertexPointMap; @@ -130,7 +133,7 @@ vertex_discrete_mean_curvature(typename boost::graph_traits::verte const FT two_pi = 2 * CGAL_PI; FT hi = 0; - for(halfedge_descriptor h : CGAL::halfedges_around_target(vd, tm)) + for(halfedge_descriptor h : CGAL::halfedges_around_target(v, tm)) { const Point_3& p = get(vpm, source(h, tm)); const Point_3& q = get(vpm, target(h, tm)); @@ -153,8 +156,8 @@ vertex_discrete_mean_curvature(typename boost::graph_traits::verte template void discrete_mean_curvature(const TriangleMesh& tm, - VertexCurvatureMap vcm, - const NamedParameters& np) + VertexCurvatureMap vcm, + const NamedParameters& np) { typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; @@ -165,14 +168,14 @@ void discrete_mean_curvature(const TriangleMesh& tm, template auto vertex_discrete_mean_curvature(typename boost::graph_traits::vertex_descriptor vd, - const TriangleMesh& tm) + const TriangleMesh& tm) { return vertex_discrete_mean_curvature(vd, tm, parameters::all_default()); } template void discrete_mean_curvature(const TriangleMesh& tm, - VertexCurvatureMap vcm) + VertexCurvatureMap vcm) { discrete_mean_curvature(tm, vcm, parameters::all_default()); } From e8602e1f4eaf77444a20c42274574bfd78c3bee1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Sat, 7 Dec 2024 23:50:02 +0100 Subject: [PATCH 11/27] Use modern typedefs --- .../internal/curvature.h | 30 ++++++++++--------- 1 file changed, 16 insertions(+), 14 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h index 9eca4794b80..992db0a56a3 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h @@ -29,15 +29,16 @@ vertex_discrete_Gaussian_curvature(typename boost::graph_traits::v const TriangleMesh& tm, const NamedParameters& np) { - typedef typename GetGeomTraits::type GeomTraits; - typedef typename GetVertexPointMap::const_type VertexPointMap; - typedef typename GeomTraits::FT FT; - typedef typename GeomTraits::Vector_3 Vector_3; - typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; - using parameters::choose_parameter; using parameters::get_parameter; + using GeomTraits = typename GetGeomTraits::type; + using tyFT = 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, tm)); @@ -95,7 +96,7 @@ void discrete_Gaussian_curvatures(const TriangleMesh& tm, VertexCurvatureMap vcm, const NamedParameters& np) { - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + using vertex_descriptor = typename boost::graph_traits::vertex_descriptor; for(vertex_descriptor v : vertices(tm)) put(vcm, v, vertex_discrete_Gaussian_curvature(v, tm, np)); @@ -114,15 +115,16 @@ vertex_discrete_mean_curvature(typename boost::graph_traits::verte const TriangleMesh& tm, const NamedParameters& np) { - typedef typename GetGeomTraits::type GeomTraits; - typedef typename GetVertexPointMap::const_type VertexPointMap; - typedef typename GeomTraits::FT FT; - typedef typename GeomTraits::Point_3 Point_3; - typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; - using parameters::choose_parameter; using parameters::get_parameter; + using GeomTraits = typename GetGeomTraits::type; + using FT = typename GeomTraits::FT; + using Point_3 = typename GeomTraits::Point_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, tm)); @@ -159,7 +161,7 @@ void discrete_mean_curvature(const TriangleMesh& tm, VertexCurvatureMap vcm, const NamedParameters& np) { - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + using vertex_descriptor = typename boost::graph_traits::vertex_descriptor; for(vertex_descriptor v : vertices(tm)) put(vcm, v, vertex_discrete_mean_curvature(v, tm, np)); From 2838ad24d9c537edb9abd64c5fc067882a516f10 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Tue, 10 Dec 2024 17:09:31 +0100 Subject: [PATCH 12/27] Delay to_double call --- Kernel_23/include/CGAL/Kernel/function_objects.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Kernel_23/include/CGAL/Kernel/function_objects.h b/Kernel_23/include/CGAL/Kernel/function_objects.h index 5a7b0774e6d..71d0f4d5130 100644 --- a/Kernel_23/include/CGAL/Kernel/function_objects.h +++ b/Kernel_23/include/CGAL/Kernel/function_objects.h @@ -965,7 +965,7 @@ 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; From 20115d3dcb70e4b1c1677837c1b8e3a087d4b98c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Tue, 10 Dec 2024 17:09:54 +0100 Subject: [PATCH 13/27] Use std::clamp --- Kernel_23/include/CGAL/Kernel/function_objects.h | 10 +--------- 1 file changed, 1 insertion(+), 9 deletions(-) diff --git a/Kernel_23/include/CGAL/Kernel/function_objects.h b/Kernel_23/include/CGAL/Kernel/function_objects.h index 71d0f4d5130..c49094aefb2 100644 --- a/Kernel_23/include/CGAL/Kernel/function_objects.h +++ b/Kernel_23/include/CGAL/Kernel/function_objects.h @@ -972,15 +972,7 @@ namespace CommonKernelFunctors { // 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; } From f7c36ea81d5c3b889df889ce798ef1cb9f363127 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Tue, 10 Dec 2024 17:10:11 +0100 Subject: [PATCH 14/27] Make the bottom and top of the ramp not have same colors --- Lab/demo/Lab/Color_map.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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 Date: Tue, 10 Dec 2024 21:17:36 +0100 Subject: [PATCH 15/27] Improve discrete curvature implementation --- .../internal/curvature.h | 268 +++++++++++++----- .../include/CGAL/Weights/cotangent_weights.h | 1 - 2 files changed, 196 insertions(+), 73 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h index 992db0a56a3..94860c496d3 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h @@ -15,42 +15,129 @@ #include -#include -#include #include +#include +#include +#include + +#include +#include namespace CGAL { namespace Polygon_mesh_processing { -namespace exprimental { -template -typename GetGeomTraits::type::FT -vertex_discrete_Gaussian_curvature(typename boost::graph_traits::vertex_descriptor v, - const TriangleMesh& tm, - const NamedParameters& np) +// Discrete Gaussian Curvature + +/** + * \ingroup PMP_vertex_angle_grp + * + * computes the sum of the angles around a vertex. + * + * @tparam PolygonMesh a model of `HalfedgeGraph` + * @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 GeomTraits = typename GetGeomTraits::type; - using tyFT = typename GeomTraits::FT; + 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::Construct_vector_3 approx_angle = gt.compute_approximate_angle_3_object(); + + FT angle_sum = 0; + for(auto h : halfedges_around_source(v, pmesh)) + { + angle_sum += approx_angle(get(vpm, target(h, pmesh)), + get(vpm, source(h, pmesh)), + get(vpm, source(prev(h,pmesh), pmesh))); + } + + return angle_sum; +} + +template +#ifdef DOXYGEN_RUNNING +FT +#else +typename GetGeomTraits::type::FT +#endif +discrete_Gaussian_curvature(typename boost::graph_traits::vertex_descriptor v, + const TriangleMesh& tm, + 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 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, tm)); - typename GeomTraits::Construct_vector_3 vector = gt.construct_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(); - typename GeomTraits::Construct_cross_product_vector_3 cross_product = gt.construct_cross_product_vector_3_object(); + 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(); - const FT two_pi = 2 * CGAL_PI; + FT angle_sum = 0; - FT ki = 0; for(halfedge_descriptor h : CGAL::halfedges_around_target(v, tm)) { if(is_border(h, tm)) @@ -62,75 +149,81 @@ vertex_discrete_Gaussian_curvature(typename boost::graph_traits::v const FT dot = scalar_product(v0, v1); const Vector_3 cross = cross_product(v0, v1); const FT sqcn = squared_length(cross); - if(dot == FT(0)) + if(is_zero(dot)) { - ki += CGAL_PI/FT(2); + angle_sum += CGAL_PI/FT(2); } else { - if(sqcn == FT(0)) + if(is_zero(sqcn)) // collinear { if(dot < 0) - ki += CGAL_PI; + angle_sum += CGAL_PI; + // else + // angle_sum += 0; } else { - ki += std::atan2(CGAL::approximate_sqrt(sqcn), dot); + angle_sum += std::atan2(CGAL::approximate_sqrt(sqcn), dot); } } } - return two_pi - ki; + Weights::Secure_cotangent_weight_with_voronoi_area wc(tm, vpm, gt); + + const FT gaussian_curvature = (2 * CGAL_PI - angle_sum) / wc.voronoi(v); + + return gaussian_curvature; } -template -auto -vertex_discrete_Gaussian_curvature(typename boost::graph_traits::vertex_descriptor v, - const TriangleMesh& tm) -{ - return vertex_discrete_Gaussian_curvature(v, tm, parameters::all_default()); -} - -template +template void discrete_Gaussian_curvatures(const TriangleMesh& tm, VertexCurvatureMap vcm, - const NamedParameters& np) + const CGAL_NP_CLASS& np = parameters::default_values()) { using vertex_descriptor = typename boost::graph_traits::vertex_descriptor; for(vertex_descriptor v : vertices(tm)) - put(vcm, v, vertex_discrete_Gaussian_curvature(v, tm, np)); + put(vcm, v, discrete_Gaussian_curvature(v, tm, np)); } -template -void discrete_Gaussian_curvatures(const TriangleMesh& tm, - VertexCurvatureMap vcm) -{ - discrete_Gaussian_curvatures(tm, vcm, parameters::all_default()); -} +// Discrete Mean Curvature -template -typename GetGeomTraits::type::FT -vertex_discrete_mean_curvature(typename boost::graph_traits::vertex_descriptor v, - const TriangleMesh& tm, - const NamedParameters& np) +template +#ifdef DOXYGEN_RUNNING +FT +#else +typename GetGeomTraits::type::FT +#endif +discrete_mean_curvature(typename boost::graph_traits::vertex_descriptor v, + const TriangleMesh& tm, + const CGAL_NP_CLASS& np = parameters::default_values()) { using parameters::choose_parameter; using parameters::get_parameter; - using GeomTraits = typename GetGeomTraits::type; + using GeomTraits = typename GetGeomTraits::type; using FT = typename GeomTraits::FT; - using Point_3 = typename GeomTraits::Point_3; + using Vector_3 = typename GeomTraits::Vector_3; - using VertexPointMap = typename GetVertexPointMap::const_type; + 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, tm)); - typename GeomTraits::Compute_squared_distance_3 squared_distance = gt.compute_squared_distance_3_object(); - typename GeomTraits::Compute_approximate_dihedral_angle_3 approximate_dihedral_angle= gt.compute_approximate_dihedral_angle_3_object(); +#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; @@ -143,7 +236,7 @@ vertex_discrete_mean_curvature(typename boost::graph_traits::verte const Point_3& s = get(vpm, target(next(opposite(h, tm), tm), tm)); const FT l = squared_distance(p,q); - FT phi = CGAL_PI * approximate_dihedral_angle(p, q, r, s) / FT(180); + FT phi = CGAL_PI * dihedral_angle(p, q, r, s) / FT(180); if(phi < 0) phi += two_pi; @@ -154,35 +247,66 @@ vertex_discrete_mean_curvature(typename boost::graph_traits::verte } 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(tm, vpm, gt); + + Vector_3 kh = vector(CGAL::NULL_VECTOR); + for(halfedge_descriptor h : CGAL::halfedges_around_target(v, tm)) + { + const vertex_descriptor v1 = source(h, tm); + + const Point_ref p0 = get(vpm, v); + const Point_ref p1 = get(vpm, v1); + + FT local_c = 0; + if(!is_border(h, tm)) + { + const vertex_descriptor v2 = target(next(h, tm), tm); + const Point_ref p2 = get(vpm, v2); + local_c += Weights::cotangent_3_clamped(p0, p2, p1, gt); + } + + if(!is_border(opposite(h, tm), tm)) + { + const vertex_descriptor v3 = target(next(opposite(h, tm), tm), tm); + 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 } -template -void discrete_mean_curvature(const TriangleMesh& tm, - VertexCurvatureMap vcm, - const NamedParameters& np) +template +void discrete_mean_curvatures(const TriangleMesh& tm, + VertexCurvatureMap vcm, + const CGAL_NP_CLASS& np = parameters::default_values()) { using vertex_descriptor = typename boost::graph_traits::vertex_descriptor; for(vertex_descriptor v : vertices(tm)) - put(vcm, v, vertex_discrete_mean_curvature(v, tm, np)); + put(vcm, v, discrete_mean_curvature(v, tm, np)); } -template -auto -vertex_discrete_mean_curvature(typename boost::graph_traits::vertex_descriptor vd, - const TriangleMesh& tm) -{ - return vertex_discrete_mean_curvature(vd, tm, parameters::all_default()); -} - -template -void discrete_mean_curvature(const TriangleMesh& tm, - VertexCurvatureMap vcm) -{ - discrete_mean_curvature(tm, vcm, parameters::all_default()); -} - -} // namespace experimental } // namespace Polygon_mesh_processing } // namespace CGAL diff --git a/Weights/include/CGAL/Weights/cotangent_weights.h b/Weights/include/CGAL/Weights/cotangent_weights.h index 4f389e3b06b..0d971f0fd17 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(); From 1189356739f455768591dd8ffb95ac220c181159 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Tue, 10 Dec 2024 21:17:57 +0100 Subject: [PATCH 16/27] Add discrete curvatures to display plugin --- .../Display/Display_property_plugin.cpp | 74 ++++++++++++++++--- 1 file changed, 63 insertions(+), 11 deletions(-) diff --git a/Lab/demo/Lab/Plugins/Display/Display_property_plugin.cpp b/Lab/demo/Lab/Plugins/Display/Display_property_plugin.cpp index f0eb9086215..4156305559e 100644 --- a/Lab/demo/Lab/Plugins/Display/Display_property_plugin.cpp +++ b/Lab/demo/Lab/Plugins/Display/Display_property_plugin.cpp @@ -22,6 +22,7 @@ #include #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; From 5ece39385f2ba40b66965a952781d62d2b7031d4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Tue, 10 Dec 2024 21:18:11 +0100 Subject: [PATCH 17/27] Fix placement of assertion --- Weights/include/CGAL/Weights/cotangent_weights.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Weights/include/CGAL/Weights/cotangent_weights.h b/Weights/include/CGAL/Weights/cotangent_weights.h index 0d971f0fd17..34a6d79f178 100644 --- a/Weights/include/CGAL/Weights/cotangent_weights.h +++ b/Weights/include/CGAL/Weights/cotangent_weights.h @@ -353,11 +353,12 @@ public: 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); From 38142206b69a71f3cec646cd7869ab347f56050f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Tue, 10 Dec 2024 21:18:28 +0100 Subject: [PATCH 18/27] Do not use an inexact SQRT if an exact one exists --- Weights/include/CGAL/Weights/internal/utils.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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)); } }; From d7791980d17ff1a69a148f7adb630cbccd00d43f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Fri, 10 Jan 2025 11:40:10 +0100 Subject: [PATCH 19/27] Move separator --- .../include/CGAL/Polygon_mesh_processing/internal/curvature.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h index 94860c496d3..622e97e32ca 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h @@ -26,8 +26,6 @@ namespace CGAL { namespace Polygon_mesh_processing { -// Discrete Gaussian Curvature - /** * \ingroup PMP_vertex_angle_grp * @@ -102,6 +100,8 @@ angle_sum(typename boost::graph_traits::vertex_descriptor v, return angle_sum; } +// Discrete Gaussian Curvature + template #ifdef DOXYGEN_RUNNING From 6d68861f1e445fedf31413dffe0b3d319646ea59 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Fri, 10 Jan 2025 11:40:50 +0100 Subject: [PATCH 20/27] Move PMP/curvature.h outside of 'internal' --- .../CGAL/Polygon_mesh_processing/{internal => }/curvature.h | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/{internal => }/curvature.h (100%) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/curvature.h similarity index 100% rename from Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/curvature.h rename to Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/curvature.h From 26005e59f3f1cb5f0997e003dbae94da0e821eee Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Fri, 10 Jan 2025 13:10:16 +0100 Subject: [PATCH 21/27] Fix compilation --- .../include/CGAL/Polygon_mesh_processing/curvature.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/curvature.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/curvature.h index 622e97e32ca..819d3c8a638 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/curvature.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/curvature.h @@ -87,7 +87,7 @@ angle_sum(typename boost::graph_traits::vertex_descriptor v, CGAL_precondition(is_valid_vertex_descriptor(v, pmesh)); - typename Geom_traits::Construct_vector_3 approx_angle = gt.compute_approximate_angle_3_object(); + 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)) From db297d7eff6c059538f36344efd209a02f36a8c6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Fri, 10 Jan 2025 13:10:33 +0100 Subject: [PATCH 22/27] Handle boundary vertices in angle_sum --- .../include/CGAL/Polygon_mesh_processing/curvature.h | 3 +++ 1 file changed, 3 insertions(+) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/curvature.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/curvature.h index 819d3c8a638..288ea9ab3f7 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/curvature.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/curvature.h @@ -92,6 +92,9 @@ angle_sum(typename boost::graph_traits::vertex_descriptor v, 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))); From c418aea7c25111436dedb873b1e085d34db4ff83 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Fri, 10 Jan 2025 13:11:21 +0100 Subject: [PATCH 23/27] Document PMP/curvature.h --- .../PackageDescription.txt | 8 +- .../Polygon_mesh_processing.txt | 11 + .../CGAL/Polygon_mesh_processing/curvature.h | 225 +++++++++++++++--- 3 files changed, 209 insertions(+), 35 deletions(-) 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 index 288ea9ab3f7..cdcf687d678 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/curvature.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/curvature.h @@ -8,10 +8,11 @@ // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial // // -// Author(s) : Mael Rouxel-Labbé +// Author(s) : Andreas Fabri, +// Mael Rouxel-Labbé -#ifndef CGAL_PMP_INTERNAL_CURVATURE_H -#define CGAL_PMP_INTERNAL_CURVATURE_H +#ifndef CGAL_PMP_CURVATURE_H +#define CGAL_PMP_CURVATURE_H #include @@ -27,10 +28,12 @@ namespace CGAL { namespace Polygon_mesh_processing { /** - * \ingroup PMP_vertex_angle_grp + * \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 `HalfedgeGraph` * @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters" * @@ -60,7 +63,6 @@ namespace Polygon_mesh_processing { * or the geometric traits class deduced from the point property map of `pmesh`. * * \warning This function involves trigonometry. - * */ template @@ -105,6 +107,45 @@ angle_sum(typename boost::graph_traits::vertex_descriptor v, // 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 `HalfedgeGraph` + * @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 @@ -113,7 +154,7 @@ FT typename GetGeomTraits::type::FT #endif discrete_Gaussian_curvature(typename boost::graph_traits::vertex_descriptor v, - const TriangleMesh& tm, + const TriangleMesh& tmesh, const CGAL_NP_CLASS& np = parameters::default_values()) { using parameters::choose_parameter; @@ -128,7 +169,7 @@ discrete_Gaussian_curvature(typename boost::graph_traits::vertex_d 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, tm)); + get_const_property_map(vertex_point, tmesh)); typename GeomTraits::Construct_vector_3 vector = gt.construct_vector_3_object(); @@ -141,13 +182,13 @@ discrete_Gaussian_curvature(typename boost::graph_traits::vertex_d FT angle_sum = 0; - for(halfedge_descriptor h : CGAL::halfedges_around_target(v, tm)) + for(halfedge_descriptor h : CGAL::halfedges_around_target(v, tmesh)) { - if(is_border(h, tm)) + if(is_border(h, tmesh)) continue; - const Vector_3 v0 = vector(get(vpm, v), get(vpm, target(next(h, tm), tm))); // p1p2 - const Vector_3 v1 = vector(get(vpm, v), get(vpm, source(h, tm))); // p1p0 + 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); @@ -172,28 +213,108 @@ discrete_Gaussian_curvature(typename boost::graph_traits::vertex_d } } - Weights::Secure_cotangent_weight_with_voronoi_area wc(tm, vpm, gt); + 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 `HalfedgeGraph` + * @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& tm, +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(tm)) - put(vcm, v, discrete_Gaussian_curvature(v, tm, np)); + 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 `HalfedgeGraph` + * @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 @@ -202,7 +323,7 @@ FT typename GetGeomTraits::type::FT #endif discrete_mean_curvature(typename boost::graph_traits::vertex_descriptor v, - const TriangleMesh& tm, + const TriangleMesh& tmesh, const CGAL_NP_CLASS& np = parameters::default_values()) { using parameters::choose_parameter; @@ -220,7 +341,7 @@ discrete_mean_curvature(typename boost::graph_traits::vertex_descr 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, tm)); + get_const_property_map(vertex_point, tmesh)); #if 0 typename GeomTraits::Compute_squared_distance_3 squared_distance = @@ -231,12 +352,12 @@ discrete_mean_curvature(typename boost::graph_traits::vertex_descr const FT two_pi = 2 * CGAL_PI; FT hi = 0; - for(halfedge_descriptor h : CGAL::halfedges_around_target(v, tm)) + for(halfedge_descriptor h : CGAL::halfedges_around_target(v, tmesh)) { - const Point_3& p = get(vpm, source(h, tm)); - const Point_3& q = get(vpm, target(h, tm)); - const Point_3& r = get(vpm, target(next(h, tm), tm)); - const Point_3& s = get(vpm, target(next(opposite(h, tm), tm), tm)); + 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); @@ -260,27 +381,27 @@ discrete_mean_curvature(typename boost::graph_traits::vertex_descr typename GeomTraits::Compute_squared_length_3 squared_length = gt.compute_squared_length_3_object(); - Weights::Secure_cotangent_weight_with_voronoi_area wc(tm, vpm, gt); + 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, tm)) + for(halfedge_descriptor h : CGAL::halfedges_around_target(v, tmesh)) { - const vertex_descriptor v1 = source(h, tm); + 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, tm)) + if(!is_border(h, tmesh)) { - const vertex_descriptor v2 = target(next(h, tm), tm); + 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, tm), tm)) + if(!is_border(opposite(h, tmesh), tmesh)) { - const vertex_descriptor v3 = target(next(opposite(h, tm), tm), tm); + 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); } @@ -297,20 +418,58 @@ discrete_mean_curvature(typename boost::graph_traits::vertex_descr #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 `HalfedgeGraph` + * @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& tm, +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(tm)) - put(vcm, v, discrete_mean_curvature(v, tm, np)); + for(vertex_descriptor v : vertices(tmesh)) + put(vcm, v, discrete_mean_curvature(v, tmesh, np)); } } // namespace Polygon_mesh_processing } // namespace CGAL -#endif //CGAL_PMP_INTERNAL_CURVATURE_H +#endif //CGAL_PMP_CURVATURE_H From 8e3befc09a65135fb21b74b622028c5807f0d29f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Fri, 10 Jan 2025 13:11:37 +0100 Subject: [PATCH 24/27] Add tests for discrete curvatures --- .../Polygon_mesh_processing/CMakeLists.txt | 1 + .../test_discrete_curvatures.cpp | 146 ++++++++++++++++++ ...test_interpolated_corrected_curvatures.cpp | 5 +- 3 files changed, 150 insertions(+), 2 deletions(-) create mode 100644 Polygon_mesh_processing/test/Polygon_mesh_processing/test_discrete_curvatures.cpp diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt b/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt index 0dd21ab22e1..e50408d7169 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); From 76e47d40ed03f7eda2edd6cee81db7146847cd4c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Fri, 10 Jan 2025 13:51:57 +0100 Subject: [PATCH 25/27] Fix header path --- Lab/demo/Lab/Plugins/Display/Display_property_plugin.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Lab/demo/Lab/Plugins/Display/Display_property_plugin.cpp b/Lab/demo/Lab/Plugins/Display/Display_property_plugin.cpp index 4156305559e..ae298d99940 100644 --- a/Lab/demo/Lab/Plugins/Display/Display_property_plugin.cpp +++ b/Lab/demo/Lab/Plugins/Display/Display_property_plugin.cpp @@ -22,7 +22,7 @@ #include #include #include -#include +#include #include #include From 07d7615d91e6c6ceec012e78ad4e9cd413d04de3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Thu, 6 Feb 2025 10:16:10 +0100 Subject: [PATCH 26/27] Require a FaceGraph for boundary checks --- .../include/CGAL/Polygon_mesh_processing/curvature.h | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/curvature.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/curvature.h index cdcf687d678..f612535dfc5 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/curvature.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/curvature.h @@ -34,7 +34,7 @@ namespace Polygon_mesh_processing { * * The angle sum is given in degrees. * - * @tparam PolygonMesh a model of `HalfedgeGraph` + * @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 @@ -114,7 +114,7 @@ angle_sum(typename boost::graph_traits::vertex_descriptor v, * * We refer to Meyer et al. \cgalCite{cgal:mdsb-ddgot-02} for the definition of discrete Gaussian curvature. * - * @tparam TriangleMesh a model of `HalfedgeGraph` + * @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 @@ -227,7 +227,7 @@ discrete_Gaussian_curvature(typename boost::graph_traits::vertex_d * * We refer to Meyer et al. \cgalCite{cgal:mdsb-ddgot-02} for the definition of discrete Gaussian curvature. * - * @tparam TriangleMesh a model of `HalfedgeGraph` + * @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, @@ -284,7 +284,7 @@ void discrete_Gaussian_curvatures(const TriangleMesh& tmesh, * * We refer to Meyer et al. \cgalCite{cgal:mdsb-ddgot-02} for the definition of discrete mean curvature. * - * @tparam TriangleMesh a model of `HalfedgeGraph` + * @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 @@ -425,7 +425,7 @@ discrete_mean_curvature(typename boost::graph_traits::vertex_descr * * We refer to Meyer et al. \cgalCite{cgal:mdsb-ddgot-02} for the definition of discrete mean curvature. * - * @tparam TriangleMesh a model of `HalfedgeGraph` + * @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, From 4cc3eef67b90a46edcaac3d0654002d6c24449c2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Thu, 6 Feb 2025 20:17:19 +0100 Subject: [PATCH 27/27] Update CHANGES.md --- Installation/CHANGES.md | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/Installation/CHANGES.md b/Installation/CHANGES.md index 7a159e087bf..fbc287d2fd1 100644 --- a/Installation/CHANGES.md +++ b/Installation/CHANGES.md @@ -3,6 +3,10 @@ ## [Release 6.1](https://github.com/CGAL/cgal/releases/tag/v6.1) +### [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.