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 ac3e94db114..e532de3e728 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 @@ -973,6 +973,11 @@ quad meshes, and meshes with n-gon faces (for n-gons, the centroid must be insid The algorithms used prove to work well in general. Also, on meshes with noise on vertex positions, they give accurate results, under the condition that the correct vertex normals are provided. +It is worth noting that the Principal Curvatures and Directions can also be estimated using the +\ref PkgJetFitting3 package, which estimates the local differential quantities of a surface at a point +using a local polynomial fitting (fitting a d-jet). Unlike the Interpolated Corrected Curvatures, +the Jet Fitting method discards topological information and thus can be used on point clouds as well. + \subsection ICCBackground Brief Background Surface curvatures are quantities that describe the local geometry of a surface. They are important in many @@ -987,7 +992,7 @@ introduce a new way to compute curvatures on polygonal meshes. The main idea in based on decoupling the normal information from the position information, which is useful for dealing with digital surfaces, or meshes with noise on vertex positions. \cgalCite{cgal:lrtc-iccmps-20} introduces some extensions to this framework, as it uses linear interpolation on the corrected normal vector field -to derive new closed form equations for the corrected curvature measures. These interpolated +to derive new closed-form equations for the corrected curvature measures. These interpolated curvature measures are the first step for computing the curvatures. For a triangle \f$ \tau_{ijk} \f$, with vertices \a i, \a j, \a k: @@ -1011,8 +1016,13 @@ solver. The interpolated curvature measures are then computed for each vertex \f$ v \f$ as the sum of the curvature measures of the faces in a ball around \f$ v \f$ weighted by the inclusion ratio of the -triangle in the ball. If the ball radius is not specified, the sum is instead computed over the incident faces -of \f$ v \f$. +triangle in the ball. This ball radius is an optional (named) parameter of the function. There are 3 +cases for the ball radius passed value: +- A positive value is passed: it is naturally used as the radius of the ball. +- 0 is passed, a small epsilon (`average_edge_length * 1e-6`) is used +(to account for the convergence of curvatures at infinitely small balls). +- It is not specified (or negative), the sum is instead computed over the incident faces +of the vertex \f$ v \f$. To get the final curvature value for a vertex \f$ v \f$, the respective interpolated curvature measure is divided by the interpolated area measure. @@ -1030,7 +1040,7 @@ These computations are performed using (on all vertices of the mesh) `CGAL::Poly where function named parameters are used to select the curvatures (and possibly directions) to be computed. An overload function with the same name but taking a given vertex is also available in case the computation should be done only for that vertex. -\subsection ICCResults Results & Performance +\subsection ICCResults Results First, \cgalFigureRef{icc_measures} illustrates various curvature measures on a triangular mesh. @@ -1075,6 +1085,121 @@ When changing the integration ball radius, we obtain a scale space of curvature \cgalFigureCaptionEnd +\subsection ICCPerformance Performance + +The implemented algorithms exhibit a linear complexity in the number of faces of the mesh. It is worth noting that +we pre-computed the vertex normals and passed them as a named parameter to the function to better estimate the +performance of the curvature computation. For the data reported in the following table, we used a machine with an +Intel Core i7-8750H CPU @ 2.20GHz, 16GB of RAM, on Windows 11, 64 bits and compiled with Visual Studio 2019. + +\cgalFigureAnchor{icc_performance_table} +
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Ball
Radius		ComputationSpot
(6k faces)		Bunny
(144K faces)		Stanford Dragon
(871K faces)		Old Age or Winter
(6M faces)
vertex
1-ring faces
(default)		Mean Curvature < 0.001 s0.019 s0.11 s2.68 s
Gaussian Curvature < 0.001 s0.017 s0.10 s2.77 s
Principal Curvatures & Directions0.002 s0.044 s0.25 s3.98 s
All (optimized for shared computations)0.003 s0.049 s0.28 s4.52 s
r = 0.1
* avg_edge_length		Mean Curvature0.017 s0.401 s2.66 s22.29 s
Gaussian Curvature0.018 s0.406 s2.63 s21.61 s
Principal Curvatures & Directions0.019 s0.430 s2.85 s23.55 s
All (optimized for shared computations)0.017 s0.428 s2.89 s24.16 s
r = 0.5
* avg_edge_length		Mean Curvature0.024 s0.388 s3.18 s22.79 s
Gaussian Curvature0.024 s0.392 s3.21 s23.58 s
Principal Curvatures & Directions0.027 s0.428 s3.41 s24.44 s
All (optimized for shared computations)0.025 s0.417 s3.44 s23.93 s
+
+\cgalFigureCaptionBegin{icc_performance_table} +Performance of the curvature computation on various meshes (in seconds). The first 4 rows show the performance of the default value for the ball radius, +which is using the 1-ring of neighboring faces around each vertex, instead of actually approximating the inclusion ratio of the faces in a ball of certain radius. +The other rows show a ball radius of 0.1 (and 0.5) scaled by the average edge length of the mesh. It is clear that using the 1-ring of faces +is much faster, but it might not be as effective when used on a noisy input mesh. +\cgalFigureCaptionEnd + + \ref BGLPropertyMaps are used to record the computed curvatures as shown in examples. In the following examples, for each property map, we associate a curvature value to each vertex. @@ -1309,7 +1434,7 @@ available on 7th of October 2020. It only uses the high level algorithm of chec is covered by a set of prisms, where each prism is an offset for an input triangle. That is, the implementation in \cgal does not use indirect predicates. -The interpolated corrected curvatures were implemented during GSoC 2022. This was implemented by Hossam Saeed and under +The interpolated corrected curvatures were implemented during GSoC 2022. This was implemented by Hossam Saeed and under the supervision of David Coeurjolly, Jaques-Olivier Lachaud, and Sébastien Loriot. The implementation is based on \cgalCite{cgal:lrtc-iccmps-20}. DGtal's implementation was also used as a reference during the project. diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/interpolated_corrected_curvatures.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/interpolated_corrected_curvatures.h index a8d45d8b867..13cad81b410 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/interpolated_corrected_curvatures.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/interpolated_corrected_curvatures.h @@ -76,20 +76,6 @@ struct Principal_curvatures_and_directions { namespace internal { -template -typename GT::FT average_edge_length(const PolygonMesh& pmesh) { - const std::size_t n = edges(pmesh).size(); - if (n == 0) - return 0; - - typename GT::FT avg_edge_length = 0; - for (auto e : edges(pmesh)) - avg_edge_length += edge_length(e, pmesh); - - avg_edge_length /= static_cast(n); - return avg_edge_length; -} - template struct Vertex_measures { typename GT::FT area_measure = 0; @@ -650,7 +636,7 @@ void interpolated_corrected_curvatures_one_vertex( typename GT::FT radius = choose_parameter(get_parameter(np, internal_np::ball_radius), -1); // calculate avg_edge_length as it is used in case radius is 0, and in the principal curvature computation later - typename GT::FT avg_edge_length = average_edge_length(pmesh); + typename GT::FT avg_edge_length = average_edge_length(pmesh); // if the radius is 0, we use a small epsilon to expand the ball (scaled with the average edge length) if (is_zero(radius)) @@ -799,7 +785,7 @@ private: // if no radius is given, we pass -1 which will make the expansion be only on the incident faces instead of a ball const FT radius = choose_parameter(get_parameter(np, internal_np::ball_radius), -1); - avg_edge_length = average_edge_length(pmesh); + avg_edge_length = average_edge_length(pmesh); set_ball_radius(radius); // check which curvature maps are provided by the user (determines which curvatures are computed) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/measure.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/measure.h index 8f1c60347ee..520ebbec9be 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/measure.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/measure.h @@ -217,6 +217,25 @@ squared_edge_length(typename boost::graph_traits::edge_descriptor e return squared_edge_length(halfedge(e, pmesh), pmesh, np); } +template +typename GetGeomTraits::type::FT +average_edge_length(const PolygonMesh& pmesh, + const NamedParameters& np = parameters::default_values()) +{ + typedef typename GetGeomTraits::type GT; + + const std::size_t n = edges(pmesh).size(); + CGAL_assertion(n > 0); + + typename GT::FT avg_edge_length = 0; + for (auto e : edges(pmesh)) + avg_edge_length += edge_length(e, pmesh, np); + + avg_edge_length /= static_cast(n); + return avg_edge_length; +} + /** * \ingroup PMP_measure_grp * diff --git a/Polyhedron/demo/Polyhedron/Plugins/Display/Display_property_plugin.cpp b/Polyhedron/demo/Polyhedron/Plugins/Display/Display_property_plugin.cpp index cbff878abb7..68826d91862 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/Display/Display_property_plugin.cpp +++ b/Polyhedron/demo/Polyhedron/Plugins/Display/Display_property_plugin.cpp @@ -90,6 +90,7 @@ private: double bI = 0.; double expand_radius = 0.; + double expand_radius_updated = false; double maxEdgeLength = -1.; Color_ramp color_ramp; @@ -912,6 +913,8 @@ private Q_SLOTS: expand_radius = (pow(base, val) - 1) * outMax / (pow(base, sliderMax) - 1); dock_widget->expandingRadiusLabel->setText(tr("Expanding Radius: %1").arg(expand_radius)); + CGAL_assertion(expand_radius >= 0); + expand_radius_updated = true; } private: @@ -933,7 +936,7 @@ private: bool non_init; SMesh::Property_map mu_i_map; std::tie(mu_i_map, non_init) = smesh.add_property_map(tied_string, 0); - if(non_init) + if(non_init || expand_radius_updated) { if(vnm_exists) { @@ -967,6 +970,8 @@ private: .ball_radius(expand_radius)); } } + + expand_radius_updated = false; } displaySMProperty(tied_string, smesh);