diff --git a/Documentation/doc/biblio/geom.bib b/Documentation/doc/biblio/geom.bib index 969a1a79b5b..f072116a2e2 100644 --- a/Documentation/doc/biblio/geom.bib +++ b/Documentation/doc/biblio/geom.bib @@ -152043,6 +152043,7 @@ pages = {179--189} Pages = {215--224}, Year = {2012}, Url = {https://monge.univ-mlv.fr/~colinde/pub/09edgewidth.pdf} +} @inproceedings{tang2009interactive, title={Interactive Hausdorff distance computation for general polygonal models}, @@ -152054,3 +152055,25 @@ pages = {179--189} year={2009}, organization={ACM} } + +@article{lachaud2020, + author = {Jacques-Olivier Lachaud and Pascal Romon and Boris Thibert and David Coeurjolly}, + journal = {Computer Graphics Forum (Proceedings of Symposium on Geometry Processing)}, + number = {5}, + title = {Interpolated corrected curvature measures for polygonal surfaces}, + volume = {39}, + month = jul, + year = {2020} +} + +@article{lachaud2022 + author = {Jacques-Olivier Lachaud and Pascal Romon and Boris Thibert}, + journal = {Discrete & Computational Geometry}, + title = {Corrected Curvature Measures}, + volume = {68}, + pages = {477-524}, + month = jul, + year = {2022}, + url = {https://doi.org/10.1007/s00454-022-00399-4} +} + diff --git a/Installation/include/CGAL/license/Polygon_mesh_processing/interpolated_corrected_curvatures.h b/Installation/include/CGAL/license/Polygon_mesh_processing/interpolated_corrected_curvatures.h new file mode 100644 index 00000000000..e484daaf5ea --- /dev/null +++ b/Installation/include/CGAL/license/Polygon_mesh_processing/interpolated_corrected_curvatures.h @@ -0,0 +1,54 @@ +// Copyright (c) 2016 GeometryFactory SARL (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org) +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial +// +// Author(s) : Andreas Fabri +// +// Warning: this file is generated, see include/CGAL/licence/README.md + +#ifndef CGAL_POLYGON_MESH_PROCESSING_INTERPOLATED_CORRECTED_CURVATURES_H +#define CGAL_POLYGON_MESH_PROCESSING_INTERPOLATED_CORRECTED_CURVATURES_H + +#include +#include + +#ifdef CGAL_POLYGON_MESH_PROCESSING_INTERPOLATED_CORRECTED_CURVATURES_COMMERCIAL_LICENSE + +# if CGAL_POLYGON_MESH_PROCESSING_INTERPOLATED_CORRECTED_CURVATURES_COMMERCIAL_LICENSE < CGAL_RELEASE_DATE + +# if defined(CGAL_LICENSE_WARNING) + + CGAL_pragma_warning("Your commercial license for CGAL does not cover " + "this release of the Polygon Mesh Processing - Interpolated Corrected Curvatures package.") +# endif + +# ifdef CGAL_LICENSE_ERROR +# error "Your commercial license for CGAL does not cover this release \ + of the Polygon Mesh Processing - Interpolated Corrected Curvatures package. \ + You get this error, as you defined CGAL_LICENSE_ERROR." +# endif // CGAL_LICENSE_ERROR + +# endif // CGAL_POLYGON_MESH_PROCESSING_INTERPOLATED_CORRECTED_CURVATURES_COMMERCIAL_LICENSE < CGAL_RELEASE_DATE + +#else // no CGAL_POLYGON_MESH_PROCESSING_INTERPOLATED_CORRECTED_CURVATURES_COMMERCIAL_LICENSE + +# if defined(CGAL_LICENSE_WARNING) + CGAL_pragma_warning("\nThe macro CGAL_POLYGON_MESH_PROCESSING_INTERPOLATED_CORRECTED_CURVATURES_COMMERCIAL_LICENSE is not defined." + "\nYou use the CGAL Polygon Mesh Processing - Interpolated Corrected Curvatures package under " + "the terms of the GPLv3+.") +# endif // CGAL_LICENSE_WARNING + +# ifdef CGAL_LICENSE_ERROR +# error "The macro CGAL_POLYGON_MESH_PROCESSING_INTERPOLATED_CORRECTED_CURVATURES_COMMERCIAL_LICENSE is not defined.\ + You use the CGAL Polygon Mesh Processing - Interpolated Corrected Curvatures package under the terms of \ + the GPLv3+. You get this error, as you defined CGAL_LICENSE_ERROR." +# endif // CGAL_LICENSE_ERROR + +#endif // no CGAL_POLYGON_MESH_PROCESSING_INTERPOLATED_CORRECTED_CURVATURES_COMMERCIAL_LICENSE + +#endif // CGAL_POLYGON_MESH_PROCESSING_INTERPOLATED_CORRECTED_CURVATURES_H diff --git a/Installation/include/CGAL/license/gpl_package_list.txt b/Installation/include/CGAL/license/gpl_package_list.txt index cd274cf329b..41612d5b9fc 100644 --- a/Installation/include/CGAL/license/gpl_package_list.txt +++ b/Installation/include/CGAL/license/gpl_package_list.txt @@ -51,6 +51,7 @@ Polygon_mesh_processing/connected_components Polygon Mesh Processing - Connected Polygon_mesh_processing/corefinement Polygon Mesh Processing - Corefinement Polygon_mesh_processing/core Polygon Mesh Processing - Core Polygon_mesh_processing/distance Polygon Mesh Processing - Distance +Polygon_mesh_processing/interpolated_corrected_curvatures Polygon Mesh Processing - Interpolated Corrected Curvatures Polygon_mesh_processing/measure Polygon Mesh Processing - Geometric Measure Polygon_mesh_processing/meshing_hole_filling Polygon Mesh Processing - Meshing and Hole Filling Polygon_mesh_processing/orientation Polygon Mesh Processing - Orientation diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/Doxyfile.in b/Polygon_mesh_processing/doc/Polygon_mesh_processing/Doxyfile.in index 937b43f1de6..c27c024648f 100644 --- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/Doxyfile.in +++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/Doxyfile.in @@ -23,6 +23,8 @@ EXCLUDE_SYMBOLS += experimental HTML_EXTRA_FILES = ${CGAL_PACKAGE_DOC_DIR}/fig/selfintersections.jpg \ ${CGAL_PACKAGE_DOC_DIR}/fig/mesh_smoothing.png \ ${CGAL_PACKAGE_DOC_DIR}/fig/shape_smoothing.png \ + ${CGAL_PACKAGE_DOC_DIR}/fig/icc_diff_radius.png \ ${CGAL_PACKAGE_DOC_DIR}/fig/decimate_cheese.png \ ${CGAL_PACKAGE_DOC_DIR}/fig/decimate_colors.png \ ${CGAL_PACKAGE_DOC_DIR}/fig/decimate_rg_joint.png + diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt b/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt index b17c41399d0..83929a69981 100644 --- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt +++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt @@ -16,6 +16,10 @@ /// Functions to triangulate faces, and to refine and fair regions of a polygon mesh. /// \ingroup PkgPolygonMeshProcessingRef +/// \defgroup PMP_corrected_curvatures_grp Corrected Curvature Computation +/// Functions to compute the corrected curvatures of a polygon mesh. +/// \ingroup PkgPolygonMeshProcessingRef + /// \defgroup PMP_normal_grp Normal Computation /// Functions to compute unit normals for individual/all vertices or faces. /// \ingroup PkgPolygonMeshProcessingRef @@ -75,7 +79,7 @@ \cgalPkgPicture{hole_filling_ico.png} \cgalPkgSummaryBegin -\cgalPkgAuthors{Sébastien Loriot, Mael Rouxel-Labbé, Jane Tournois, and Ilker %O. Yaz} +\cgalPkgAuthors{David Coeurjolly, Jaques-Olivier Lachaud, Konstantinos Katriopla, Sébastien Loriot, Mael Rouxel-Labbé, Hossam Saeed, Jane Tournois, and Ilker %O. Yaz} \cgalPkgDesc{This package provides a collection of methods and classes for polygon mesh processing, ranging from basic operations on simplices, to complex geometry processing algorithms such as Boolean operations, remeshing, repairing, collision and intersection detection, and more.} @@ -208,6 +212,17 @@ The page \ref bgl_namedparameters "Named Parameters" describes their usage. - \link PMP_locate_grp Nearest Face Location Queries \endlink - \link PMP_locate_grp Random Location Generation \endlink +\cgalCRPSection{Corrected Curvatures} +- `CGAL::Polygon_mesh_processing::interpolated_corrected_mean_curvature()` +- `CGAL::Polygon_mesh_processing::interpolated_corrected_Gaussian_curvature()` +- `CGAL::Polygon_mesh_processing::interpolated_corrected_principal_curvatures_and_directions()` +- `CGAL::Polygon_mesh_processing::interpolated_corrected_curvatures()` +- `CGAL::Polygon_mesh_processing::interpolated_corrected_mean_curvature_one_vertex()` +- `CGAL::Polygon_mesh_processing::interpolated_corrected_Gaussian_curvature_one_vertex()` +- `CGAL::Polygon_mesh_processing::interpolated_corrected_principal_curvatures_and_directions_one_vertex()` +- `CGAL::Polygon_mesh_processing::interpolated_corrected_curvatures_one_vertex()` +- `CGAL::Polygon_mesh_processing::Principal_curvatures_and_directions` + \cgalCRPSection{Normal Computation Functions} - `CGAL::Polygon_mesh_processing::compute_face_normal()` - `CGAL::Polygon_mesh_processing::compute_face_normals()` 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 f414898362f..67f2405290a 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 @@ -4,7 +4,7 @@ namespace CGAL { \anchor Chapter_PolygonMeshProcessing \cgalAutoToc -\authors Sébastien Loriot, Mael Rouxel-Labbé, Jane Tournois, Ilker %O. Yaz +\authors David Coeurjolly, Jaques-Olivier Lachaud, Konstantinos Katriopla, Sébastien Loriot, Mael Rouxel-Labbé, Hossam Saeed, Jane Tournois, and Ilker %O. Yaz \image html neptun_head.jpg \image latex neptun_head.jpg @@ -49,6 +49,7 @@ mesh, which includes point location and self intersection tests. - \ref PMPCombinatorialRepair : repair of polygon meshes and polygon soups. - \ref PMPGeometricRepair : repair of the geometry of polygon meshes. - \ref PMPNormalComp : normal computation at vertices and on faces of a polygon mesh. +- \ref PMPICC : computing curvatures (mean, gaussian, principal) on a polygon mesh. - \ref PMPSlicer : functor able to compute the intersections of a polygon mesh with arbitrary planes (slicer). - \ref PMPConnectedComponents : methods to deal with connected components of a polygon mesh (extraction, marks, removal, ...). @@ -935,6 +936,128 @@ not provide storage for the normals. \cgalExample{Polygon_mesh_processing/compute_normals_example_Polyhedron.cpp} +**************************************** +\section PMPICC Computing Curvatures + +This package provides methods to compute curvatures on polygonal meshes based on Interpolated +Corrected Curvatures on Polyhedral Surfaces \cgalCite{lachaud2020}. This includes mean curvature, +Gaussian curvature, principal curvatures and directions. These can be computed on triangle meshes, +quad meshes, and meshes with n-gon faces (for n-gons, the centroid must be inside the n-gon face). +The algorithms used prove to work well in general. Also, on meshes with noise on vertex positions, +they give accurate results, on the condition that the correct vertex normals are provided. + +\subsection ICCBackground Brief Background + +Surface curvatures are quantities that describe the local geometry of a surface. They are important in many +geometry processing applications. As surfaces are 2-dimensional objects (embedded in 3D), they can bend +in 2 independent directions. These directions are called principal directions, and the amount of bending +in each direction is called the principal curvature: \f$ k_1 \f$ and \f$ k_2 \f$ (denoting max and min +curvatures). Curvature is usually expressed as scalar quantities like the mean curvature \f$ H \f$ and +the Gaussian curvature \f$ K \f$ which are defined in terms of the principal curvatures. + +The algorithms are based on the two papers \cgalCite{lachaud2022} and \cgalCite{lachaud2020}. They +introduce a new way to compute curvatures on polygonal meshes. The main idea in \cgalCite{lachaud2022} is +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{lachaud2020} 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 +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: + +\f[ + \begin{align*} + \mu^{(0)}(\tau_{ijk}) = &\frac{1}{2} \langle \bar{\mathbf{u}} \mid (\mathbf{x}_j - \mathbf{x}_i) \times (\mathbf{x}_k - \mathbf{x}_i) \rangle, \\ + \mu^{(1)}(\tau_{ijk}) = &\frac{1}{2} \langle \bar{\mathbf{u}} \mid (\mathbf{u}_k - \mathbf{u}_j) \times \mathbf{x}_i + (\mathbf{u}_i - \mathbf{u}_k) \times \mathbf{x}_j + (\mathbf{u}_j - \mathbf{u}_i) \times \mathbf{x}_k \rangle, \\ + \mu^{(2)}(\tau_{ijk}) = &\frac{1}{2} \langle \mathbf{u}_i \mid \mathbf{u}_j \times \mathbf{u}_k \rangle, \\ + \mu^{\mathbf{X},\mathbf{Y}}(\tau_{ijk}) = & \frac{1}{2} \big\langle \bar{\mathbf{u}} \big| \langle \mathbf{Y} | \mathbf{u}_k -\mathbf{u}_i \rangle \mathbf{X} \times (\mathbf{x}_j - \mathbf{x}_i) \big\rangle + -\frac{1}{2} \big\langle \bar{\mathbf{u}} \big| \langle \mathbf{Y} | \mathbf{u}_j -\mathbf{u}_i \rangle \mathbf{X} \times (\mathbf{x}_k - \mathbf{x}_i) \big\rangle, + \end{align*} +\f] +where \f$ \langle \cdot \mid \cdot \rangle \f$ denotes the usual scalar product, +\f$ \bar{\mathbf{u}}=\frac{1}{3}( \mathbf{u}_i + \mathbf{u}_j + \mathbf{u}_k )\f$. + +The first measure \f$ \mu^{(0)} \f$ is the area measure of the triangle, and the measures \f$ \mu^{(1)} \f$ and +\f$ \mu^{(2)} \f$ are the mean and Gaussian corrected curvature measures of the triangle. The last measure +\f$ \mu^{\mathbf{X},\mathbf{Y}} \f$ is the anisotropic corrected curvature measure of the triangle. The +anisotropic measure is later used to compute the principal curvatures and directions through an eigenvalue +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 over the incident faces +of \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. + +\f[ +\mu^{(k)}( B ) = \sum_{\tau : \text{triangle} } \mu^{(k)}( \tau ) \frac{\mathrm{Area}( \tau \cap B )}{\mathrm{Area}(\tau)}. +\f] + +\subsection ICCAPI API + +The implementation is generic in terms of mesh data structure. It can be used on `Surface_mesh`, +`Polyhedron_3` and other polygonal mesh structures based on the concept `FaceGraph`. + +These computations are performed using (on all vertices of the mesh): +- `CGAL::Polygon_mesh_processing::interpolated_corrected_mean_curvature()` +- `CGAL::Polygon_mesh_processing::interpolated_corrected_Gaussian_curvature()` +- `CGAL::Polygon_mesh_processing::interpolated_corrected_principal_curvatures_and_directions()` +- `CGAL::Polygon_mesh_processing::interpolated_corrected_curvatures()` + +Where it is recommended to use the last function for computing multiple curvatures (for example: mean and +Gaussian) as the implementation performs the shared computations only once, making it more efficient. + +Similarly, we can use the following functions to compute curvatures on a specific vertex: +- `CGAL::Polygon_mesh_processing::interpolated_corrected_mean_curvature_one_vertex()` +- `CGAL::Polygon_mesh_processing::interpolated_corrected_Gaussian_curvature_one_vertex()` +- `CGAL::Polygon_mesh_processing::interpolated_corrected_principal_curvatures_and_directions_one_vertex()` +- `CGAL::Polygon_mesh_processing::interpolated_corrected_curvatures_one_vertex()` + +\subsection ICCResults Results & Performance + +**To be updated** + +\cgalFigureRef{icc_diff_radius} shows how the mean curvature changes depending on +the named parameter `ball_radius`, which can be set to a value > 0 to get a smoother +distribution of values and "diffuses" the extreme values of curvatures across the mesh. + +\cgalFigureAnchor{icc_diff_radius} +
+ +
+\cgalFigureCaptionBegin{icc_diff_radius} +The mean curvature on a mesh with different values for the ball radius +parameter: (a) R = 0, (b) R = 0.025, (c) R = 0.05, (d) R = 0.16. Note that the max +edge length is 0.031 and the size of the bounding box of the mesh is 1 x .7 x .8. +\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. + +\subsection ICCExampleSM Interpolated Corrected Curvatures on a Surface Mesh Example + +The following example illustrates how to +compute the curvatures on vertices +and store them in the property maps provided by the class `Surface_mesh`. + +\cgalExample{Polygon_mesh_processing/interpolated_corrected_curvatures_SM.cpp} + +\subsection ICCExamplePH Interpolated Corrected Curvatures on a Polyhedron Example + +The following example illustrates how to +compute the curvatures on vertices +and store them in dynamic property maps as the class `Polyhedron_3` does +not provide storage for the curvatures. + +\cgalExample{Polygon_mesh_processing/interpolated_corrected_curvatures_PH.cpp} + +\subsection ICCExampleSV Interpolated Corrected Curvatures on a Vertex Example + +The following example illustrates how to +compute the curvatures on a specific vertex + +\cgalExample{Polygon_mesh_processing/interpolated_corrected_curvatures_vertex.cpp} **************************************** \section PMPSlicer Slicer @@ -1143,6 +1266,10 @@ 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 +supervision of David Coeurjolly, Jaques-Olivier Lachaud, and Sébastien Loriot. The implementation is based on \cgalCite{lachaud2020}. +DGtal's implementation was also +used as a reference during the project. */ } /* namespace CGAL */ diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/examples.txt b/Polygon_mesh_processing/doc/Polygon_mesh_processing/examples.txt index 146e1261ff4..6fdeda48ca4 100644 --- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/examples.txt +++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/examples.txt @@ -19,6 +19,9 @@ \example Polygon_mesh_processing/refine_fair_example.cpp \example Polygon_mesh_processing/mesh_slicer_example.cpp \example Polygon_mesh_processing/isotropic_remeshing_example.cpp +\example Polygon_mesh_processing/interpolated_corrected_curvatures_SM.cpp +\example Polygon_mesh_processing/interpolated_corrected_curvatures_PH.cpp +\example Polygon_mesh_processing/interpolated_corrected_curvatures_vertex.cpp \example Polygon_mesh_processing/delaunay_remeshing_example.cpp \example Polygon_mesh_processing/compute_normals_example_Polyhedron.cpp \example Polygon_mesh_processing/hausdorff_distance_remeshing_example.cpp diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt b/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt index cf1424b05c9..28d0550a31c 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt @@ -74,6 +74,12 @@ if(TARGET CGAL::Eigen3_support) target_link_libraries(delaunay_remeshing_example PUBLIC CGAL::Eigen3_support) create_single_source_cgal_program("remesh_almost_planar_patches.cpp") target_link_libraries(remesh_almost_planar_patches PUBLIC CGAL::Eigen3_support) + create_single_source_cgal_program("interpolated_corrected_curvatures_SM.cpp") + target_link_libraries(interpolated_corrected_curvatures_SM PUBLIC CGAL::Eigen3_support) + create_single_source_cgal_program("interpolated_corrected_curvatures_PH.cpp") + target_link_libraries(interpolated_corrected_curvatures_PH PUBLIC CGAL::Eigen3_support) + create_single_source_cgal_program("interpolated_corrected_curvatures_vertex.cpp") + target_link_libraries(interpolated_corrected_curvatures_vertex PUBLIC CGAL::Eigen3_support) else() message(STATUS "NOTICE: Examples that use Eigen will not be compiled.") endif() diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/interpolated_corrected_curvatures_PH.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/interpolated_corrected_curvatures_PH.cpp new file mode 100644 index 00000000000..c93de26f26c --- /dev/null +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/interpolated_corrected_curvatures_PH.cpp @@ -0,0 +1,72 @@ +#include +#include +#include +#include +#include + +#include + +#include +#include + +namespace PMP = CGAL::Polygon_mesh_processing; + +typedef CGAL::Exact_predicates_inexact_constructions_kernel Epic_kernel; +typedef CGAL::Polyhedron_3 Polyhedron; +typedef boost::graph_traits::vertex_descriptor vertex_descriptor; + +int main(int argc, char* argv[]) +{ + Polyhedron polyhedron; + const std::string filename = (argc > 1) ? + argv[1] : + CGAL::data_file_path("meshes/sphere.off"); + + if (!CGAL::IO::read_polygon_mesh(filename, polyhedron)) + { + std::cerr << "Invalid input file." << std::endl; + return EXIT_FAILURE; + } + + boost::property_map>::type + mean_curvature_map = get(CGAL::dynamic_vertex_property_t(), polyhedron), + Gaussian_curvature_map = get(CGAL::dynamic_vertex_property_t(), polyhedron); + + boost::property_map>>::type + principal_curvatures_and_directions_map = + get(CGAL::dynamic_vertex_property_t>(), polyhedron); + + PMP::interpolated_corrected_mean_curvature(polyhedron, mean_curvature_map); + + PMP::interpolated_corrected_Gaussian_curvature(polyhedron, Gaussian_curvature_map); + + PMP::interpolated_corrected_principal_curvatures_and_directions(polyhedron, principal_curvatures_and_directions_map); + + // uncomment this to compute a curvature while specifying named parameters + // Example: an expansion ball radius of 0.5 and a vertex normals map (does not have to depend on positions) + + /*std::unordered_map vnm; + + PMP::interpolated_corrected_mean_curvature( + polyhedron, + mean_curvature_map, + CGAL::parameters::ball_radius(0.5).vertex_normal_map(boost::make_assoc_property_map(vnm)) + );*/ + + // This function can be used to compute multiple curvature types by specifiying them as named parameters + // This is more efficient than computing each one separately (shared computations). + PMP::interpolated_corrected_curvatures( + polyhedron, + CGAL::parameters::vertex_mean_curvature_map(mean_curvature_map) + .vertex_principal_curvatures_and_directions_map(principal_curvatures_and_directions_map)); + + int i = 0; + for (vertex_descriptor v : vertices(polyhedron)) + { + auto PC = get(principal_curvatures_and_directions_map, v); + std::cout << i << ": HC = " << get(mean_curvature_map, v) + << ", GC = " << get(Gaussian_curvature_map, v) << "\n" + << ", PC = [ " << PC.min_curvature << " , " << PC.max_curvature << " ]\n"; + i++; + } +} diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/interpolated_corrected_curvatures_SM.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/interpolated_corrected_curvatures_SM.cpp new file mode 100644 index 00000000000..dae2ec0716e --- /dev/null +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/interpolated_corrected_curvatures_SM.cpp @@ -0,0 +1,94 @@ +#include +#include +#include +#include +#include + +#include + +#include + +namespace PMP = CGAL::Polygon_mesh_processing; + +typedef CGAL::Exact_predicates_inexact_constructions_kernel Epic_kernel; +typedef CGAL::Surface_mesh Surface_Mesh; +typedef boost::graph_traits::vertex_descriptor vertex_descriptor; + +int main(int argc, char* argv[]) +{ + Surface_Mesh smesh; + const std::string filename = (argc > 1) ? + argv[1] : + CGAL::data_file_path("meshes/sphere.off"); + + if (!CGAL::IO::read_polygon_mesh(filename, smesh)) + { + std::cerr << "Invalid input file." << std::endl; + return EXIT_FAILURE; + } + + // creating and tying surface mesh property maps for curvatures (with defaults = 0) + bool created = false; + Surface_Mesh::Property_map + mean_curvature_map, Gaussian_curvature_map; + + boost::tie(mean_curvature_map, created) = + smesh.add_property_map("v:mean_curvature_map", 0); + assert(created); + + boost::tie(Gaussian_curvature_map, created) = + smesh.add_property_map("v:Gaussian_curvature_map", 0); + assert(created); + + // we use a tuple of 2 scalar values and 2 vectors for principal curvatures and directions + Surface_Mesh::Property_map> + principal_curvatures_and_directions_map; + + boost::tie(principal_curvatures_and_directions_map, created) = + smesh.add_property_map> + ("v:principal_curvatures_and_directions_map", { 0, 0, + Epic_kernel::Vector_3(0,0,0), + Epic_kernel::Vector_3(0,0,0) }); + assert(created); + + // user can call these fucntions to compute a specfic curvature type on each vertex. + // (Note: if no ball radius is specified, the measure expansion of each vertex happens by + // summing measures on faces adjacent to each vertex.) + PMP::interpolated_corrected_mean_curvature(smesh, mean_curvature_map); + + PMP::interpolated_corrected_Gaussian_curvature(smesh, Gaussian_curvature_map); + + PMP::interpolated_corrected_principal_curvatures_and_directions(smesh, principal_curvatures_and_directions_map); + + // uncomment this to compute a curvature while specifying named parameters + // Example: an expansion ball radius of 0.5 and a vertex normals map (does not have to depend on positions) + + /*Surface_Mesh::Property_map vnm; + boost::tie(vnm, created) = smesh.add_property_map( + "v:vnm", Epic_kernel::Vector_3(0, 0, 0) + ); + + assert(created); + + PMP::interpolated_corrected_mean_curvature( + smesh, + mean_curvature_map, + CGAL::parameters::ball_radius(0.5).vertex_normal_map(vnm) + );*/ + + // This function can be used to compute multiple curvature types by specifiying them as named parameters + // This is more efficient than computing each one separately (shared computations). + PMP::interpolated_corrected_curvatures( + smesh, + CGAL::parameters::vertex_mean_curvature_map(mean_curvature_map) + .vertex_principal_curvatures_and_directions_map(principal_curvatures_and_directions_map) + ); + + for (vertex_descriptor v : vertices(smesh)) + { + auto PC = principal_curvatures_and_directions_map[v]; + std::cout << v.idx() << ": HC = " << mean_curvature_map[v] + << ", GC = " << Gaussian_curvature_map[v] << "\n" + << ", PC = [ " << PC.min_curvature << " , " << PC.max_curvature << " ]\n"; + } +} diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/interpolated_corrected_curvatures_vertex.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/interpolated_corrected_curvatures_vertex.cpp new file mode 100644 index 00000000000..c6bef5faa74 --- /dev/null +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/interpolated_corrected_curvatures_vertex.cpp @@ -0,0 +1,58 @@ +#include +#include +#include +#include + +#include + +#include + +namespace PMP = CGAL::Polygon_mesh_processing; + +typedef CGAL::Exact_predicates_inexact_constructions_kernel Epic_kernel; +typedef CGAL::Surface_mesh Surface_Mesh; +typedef boost::graph_traits::vertex_descriptor vertex_descriptor; + +int main(int argc, char* argv[]) +{ + // instantiating and reading mesh + Surface_Mesh smesh; + const std::string filename = (argc > 1) ? + argv[1] : + CGAL::data_file_path("meshes/sphere.off"); + + if (!CGAL::IO::read_polygon_mesh(filename, smesh)) + { + std::cerr << "Invalid input file." << std::endl; + return EXIT_FAILURE; + } + + // loop over vertices and use vertex_descriptor to compute a curvature on one vertex + for (vertex_descriptor v : vertices(smesh)) + { + double h = PMP::interpolated_corrected_mean_curvature_one_vertex(smesh, v); + double g = PMP::interpolated_corrected_Gaussian_curvature_one_vertex(smesh, v); + PMP::Principal_curvatures_and_directions p = + PMP::interpolated_corrected_principal_curvatures_and_directions_one_vertex(smesh, v); + + // we can also specify a ball radius for expansion and a user defined vertex normals map using + // named parameters. Refer to interpolated_corrected_curvatures_SM.cpp to see example usage. + + // Can also use interpolated_corrected_curvatures_one_vertex() to compute multiple curvatures + // on the vertex at the same time. This is more efficient than computing each one separately. + // The following commented lines show this (all mentioned named parameters work on it as well) + // we specify which curvatures we want to compute by passing pointers as named parameters + // as shown. These pointers are used for storing the result as well. in this example we + // selected mean and Gaussian curvatures + // PMP::interpolated_corrected_curvatures_one_vertex( + // smesh, + // v, + // CGAL::parameters::vertex_mean_curvature(&h) + // .vertex_Gaussian_curvature(&g) + // ); + + std::cout << v.idx() << ": HC = " << h + << ", GC = " << g << "\n" + << ", PC = [ " << p.min_curvature << " , " << p.max_curvature << " ]\n"; + } +} 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 new file mode 100644 index 00000000000..4b940f157e8 --- /dev/null +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/interpolated_corrected_curvatures.h @@ -0,0 +1,1659 @@ +// Copyright (c) 2022 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) : Hossam Saeed +// + +#ifndef CGAL_POLYGON_MESH_PROCESSING_INTERPOLATED_CORRECTED_CURVATURES_H +#define CGAL_POLYGON_MESH_PROCESSING_INTERPOLATED_CORRECTED_CURVATURES_H + +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +namespace CGAL { + +namespace Polygon_mesh_processing { + +/** + * \ingroup PMP_corrected_curvatures_grp + * + * \brief a struct for storing principal curvatures and directions. + * + * @tparam GT is the geometric traits class, model of `Kernel`. +*/ +template +struct Principal_curvatures_and_directions { + + /// min curvature magnitude + typename GT::FT min_curvature; + + /// max curvature magnitude + typename GT::FT max_curvature; + + /// min curvature direction vector + typename GT::Vector_3 min_direction; + + /// max curvature direction vector + typename GT::Vector_3 max_direction; + + Principal_curvatures_and_directions() { + min_curvature = 0; + max_curvature = 0; + min_direction = typename GT::Vector_3(0, 0, 0); + max_direction = typename GT::Vector_3(0, 0, 0); + } + + Principal_curvatures_and_directions( + typename GT::FT min_curvature, + typename GT::FT max_curvature, + typename GT::Vector_3 min_direction, + typename GT::Vector_3 max_direction + ) { + this->min_curvature = min_curvature; + this->max_curvature = max_curvature; + this->min_direction = min_direction; + this->max_direction = max_direction; + } +}; + +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; + typename GT::FT mean_curvature_measure = 0; + typename GT::FT gaussian_curvature_measure = 0; + std::array anisotropic_measure = { 0, 0, 0, + 0, 0, 0, + 0, 0, 0 }; +}; + +template +typename GT::FT interpolated_corrected_area_measure_face(const std::vector& u, + const std::vector& x) +{ + const std::size_t n = x.size(); + CGAL_precondition(u.size() == n); + CGAL_precondition(n >= 3); + + typename GT::Construct_cross_product_vector_3 cross_product; + + // Triangle: use triangle formula + if (n == 3) + { + const typename GT::Vector_3 um = (u[0] + u[1] + u[2]) / 3.0; + return 0.5 * um * cross_product(x[1] - x[0], x[2] - x[0]); + } + // Quad: use the bilinear interpolation formula + else if (n == 4) + { + // for the formulas below, values of verices 2 & 3 are swapped (compared to paper) to correct order. + // the indices in paper vs in here are: 00 = 0, 10 = 1, 11 = 2, 01 = 3 + return (1.0 / 36.0) * ( + (4 * u[0] + 2 * u[1] + 2 * u[3] + u[2]) * cross_product(x[1] - x[0], x[3] - x[0]) + + (2 * u[0] + 4 * u[1] + u[3] + 2 * u[2]) * cross_product(x[1] - x[0], x[2] - x[1]) + + (2 * u[0] + u[1] + 4 * u[3] + 2 * u[2]) * cross_product(x[2] - x[3], x[3] - x[0]) + + (u[0] + 2 * u[1] + 2 * u[3] + 4 * u[2]) * cross_product(x[2] - x[3], x[2] - x[1]) + ); + } + // N-gon: split into n triangles by polygon center and use triangle formula for each + else + { + typename GT::FT mu0 = 0; + + // getting center of points + typename GT::Vector_3 xc = + std::accumulate(x.begin(), x.end(), typename GT::Vector_3(0, 0, 0)); + xc /= static_cast(n); + + // getting unit average normal of points + typename GT::Vector_3 uc = + std::accumulate(u.begin(), u.end(), typename GT::Vector_3(0, 0, 0)); + uc /= sqrt(uc * uc); + + // summing each triangle's measure after triangulation by barycenter split. + for (std::size_t i = 0; i < n; i++) + { + mu0 += interpolated_corrected_area_measure_face( + std::vector {u[i], u[(i + 1) % n], uc}, + std::vector {x[i], x[(i + 1) % n], xc} + ); + } + return mu0; + } +} + +template +typename GT::FT interpolated_corrected_mean_curvature_measure_face(const std::vector& u, + const std::vector& x) +{ + const std::size_t n = x.size(); + CGAL_precondition(u.size() == n); + CGAL_precondition(n >= 3); + + typename GT::Construct_cross_product_vector_3 cross_product; + + // Triangle: use triangle formula + if (n == 3) + { + const typename GT::Vector_3 um = (u[0] + u[1] + u[2]) / 3.0; + + return 0.5 * um * (cross_product(u[2] - u[1], x[0]) + + cross_product(u[0] - u[2], x[1]) + + cross_product(u[1] - u[0], x[2])); + } + // Quad: use the bilinear interpolation formula + else if (n == 4) + { + // for the formulas below, values of verices 2 & 3 are swapped (compared to paper) to correct order. + // the indices in paper vs in here are: 00 = 0, 10 = 1, 11 = 2, 01 = 3 + + const typename GT::Vector_3 u02 = u[2] - u[0]; + const typename GT::Vector_3 u13 = u[3] - u[1]; + const typename GT::Vector_3 x0_cross = cross_product(u13, x[0]); + const typename GT::Vector_3 x1_cross = -cross_product(u02, x[1]); + const typename GT::Vector_3 x3_cross = cross_product(u02, x[3]); + const typename GT::Vector_3 x2_cross = -cross_product(u13, x[2]); + + return (1.0 / 12.0) * ( + u[0] * (2 * x0_cross - cross_product((u[3] + u[2]), x[1]) + cross_product((u[1] + u[2]), x[3]) + x2_cross) + + u[1] * (cross_product((u[3] + u[2]), x[0]) + 2 * x1_cross + x3_cross - cross_product((u[0] + u[3]), x[2])) + + u[3] * (-cross_product((u[1] + u[2]), x[0]) + x1_cross + 2 * x3_cross + cross_product((u[0] + u[1]), x[2])) + + u[2] * (x0_cross + cross_product((u[0] + u[3]), x[1]) - cross_product((u[0] + u[1]), x[3]) + 2 * x2_cross) + ); + } + // N-gon: split into n triangles by polygon center and use triangle formula for each + else + { + typename GT::FT mu1 = 0; + + // getting center of points + typename GT::Vector_3 xc = + std::accumulate(x.begin(), x.end(), typename GT::Vector_3(0, 0, 0)); + xc /= static_cast(n); + + // getting unit average normal of points + typename GT::Vector_3 uc = + std::accumulate(u.begin(), u.end(), typename GT::Vector_3(0, 0, 0)); + uc /= sqrt(uc * uc); + + // summing each triangle's measure after triangulation by barycenter split. + for (std::size_t i = 0; i < n; i++) + { + mu1 += interpolated_corrected_mean_curvature_measure_face( + std::vector {u[i], u[(i + 1) % n], uc}, + std::vector {x[i], x[(i + 1) % n], xc} + ); + } + return mu1; + } +} + +template +typename GT::FT interpolated_corrected_Gaussian_curvature_measure_face(const std::vector& u) +{ + const std::size_t n = u.size(); + CGAL_precondition(n >= 3); + + typename GT::Construct_cross_product_vector_3 cross_product; + + // Triangle: use triangle formula + if (n == 3) + { + return 0.5 * u[0] * cross_product(u[1], u[2]); + } + // Quad: use the bilinear interpolation formula + else if (n == 4) + { + // for the formulas below, values of verices 2 & 3 are swapped (compared to paper) to correct order. + // the indices in paper vs in here are: 00 = 0, 10 = 1, 11 = 2, 01 = 3 + return (1.0 / 36.0) * ( + (4 * u[0] + 2 * u[1] + 2 * u[3] + u[2]) * cross_product(u[1] - u[0], u[3] - u[0]) + + (2 * u[0] + 4 * u[1] + u[3] + 2 * u[2]) * cross_product(u[1] - u[0], u[2] - u[1]) + + (2 * u[0] + u[1] + 4 * u[3] + 2 * u[2]) * cross_product(u[2] - u[3], u[3] - u[0]) + + (u[0] + 2 * u[1] + 2 * u[3] + 4 * u[2]) * cross_product(u[2] - u[3], u[2] - u[1]) + ); + } + // N-gon: split into n triangles by polygon center and use triangle formula for each + else + { + typename GT::FT mu2 = 0; + + // getting unit average normal of points + typename GT::Vector_3 uc = + std::accumulate(u.begin(), u.end(), typename GT::Vector_3(0, 0, 0)); + uc /= sqrt(uc * uc); + + // summing each triangle's measure after triangulation by barycenter split. + for (std::size_t i = 0; i < n; i++) + { + mu2 += interpolated_corrected_Gaussian_curvature_measure_face( + std::vector {u[i], u[(i + 1) % n], uc} + ); + } + return mu2; + } +} + +template +std::array interpolated_corrected_anisotropic_measure_face(const std::vector& u, + const std::vector& x) +{ + const std::size_t n = x.size(); + CGAL_precondition(u.size() == n); + CGAL_precondition(n >= 3); + + typename GT::Construct_cross_product_vector_3 cross_product; + std::array muXY{ 0 }; + + // Triangle: use triangle formula + if (n == 3) + { + const typename GT::Vector_3 u01 = u[1] - u[0]; + const typename GT::Vector_3 u02 = u[2] - u[0]; + const typename GT::Vector_3 x01 = x[1] - x[0]; + const typename GT::Vector_3 x02 = x[2] - x[0]; + const typename GT::Vector_3 um = (u[0] + u[1] + u[2]) / 3.0; + + for (unsigned int ix = 0; ix < 3; ix++) + { + typename GT::Vector_3 X; + if (ix == 0) + X = typename GT::Vector_3(1, 0, 0); + if (ix == 1) + X = typename GT::Vector_3(0, 1, 0); + if (ix == 2) + X = typename GT::Vector_3(0, 0, 1); + + for (unsigned int iy = 0; iy < 3; iy++) + muXY[ix * 3 + iy] = 0.5 * um * (cross_product(u02[iy] * X, x01) - cross_product(u01[iy] * X, x02)); + } + } + // Quad: use the bilinear interpolation formula + else if (n == 4) + { + // for the formulas below, values of verices 2 & 3 are swapped (compared to paper) to correct order. + // the indices in paper vs in here are: 00 = 0, 10 = 1, 11 = 2, 01 = 3 + for (unsigned int ix = 0; ix < 3; ix++) + { + typename GT::Vector_3 X; + if (ix == 0) + X = typename GT::Vector_3(1, 0, 0); + if (ix == 1) + X = typename GT::Vector_3(0, 1, 0); + if (ix == 2) + X = typename GT::Vector_3(0, 0, 1); + + const typename GT::Vector_3 u0xX = cross_product(u[0], X); + const typename GT::Vector_3 u1xX = cross_product(u[1], X); + const typename GT::Vector_3 u2xX = cross_product(u[2], X); + const typename GT::Vector_3 u3xX = cross_product(u[3], X); + + for (unsigned int iy = 0; iy < 3; iy++) + muXY[ix * 3 + iy] = (1.0 / 72.0) * ( + + u[0][iy] * (u0xX * (-x[0] - 11 * x[1] + 13 * x[3] - x[2]) + + u1xX * (-5 * x[0] - 7 * x[1] + 11 * x[3] + x[2]) + + u3xX * (x[0] - 7 * x[1] + 11 * x[3] - 5 * x[2]) + + u2xX * (-x[0] - 5 * x[1] + 7 * x[3] - x[2]) + ) + + u[1][iy] * (u0xX * (13 * x[0] - x[1] - 7 * x[3] - 5 * x[2]) + + u1xX * (17 * x[0] - 5 * x[1] - 5 * x[3] - 7 * x[2]) + + u3xX * (5 * x[0] + x[1] + x[3] - 7 * x[2]) + + u2xX * (7 * x[0] - x[1] + 5 * x[3] - 11 * x[2]) + ) + + u[2][iy] * (u0xX * (-11 * x[0] + 5 * x[1] - x[3] + 7 * x[2]) + + u1xX * (-7 * x[0] + x[1] + x[3] + 5 * x[2]) + + u3xX * (-7 * x[0] - 5 * x[1] - 5 * x[3] + 17 * x[2]) + + u2xX * (-5 * x[0] - 7 * x[1] - x[3] + 13 * x[2]) + ) + + u[3][iy] * (u0xX * (-x[0] + 7 * x[1] - 5 * x[3] - x[2]) + + u1xX * (-5 * x[0] + 11 * x[1] - 7 * x[3] + x[2]) + + u3xX * (x[0] + 11 * x[1] - 7 * x[3] - 5 * x[2]) + + u2xX * (-x[0] + 13 * x[1] - 11 * x[3] - x[2]) + ) + + ); + } + } + // N-gon: split into n triangles by polygon center and use triangle formula for each + else + { + // getting center of points + typename GT::Vector_3 xc = + std::accumulate(x.begin(), x.end(), typename GT::Vector_3(0, 0, 0)); + xc /= static_cast(n); + + // getting unit average normal of points + typename GT::Vector_3 uc = + std::accumulate(u.begin(), u.end(), typename GT::Vector_3(0, 0, 0)); + uc /= sqrt(uc * uc); + + // summing each triangle's measure after triangulation by barycenter split. + for (std::size_t i = 0; i < n; i++) + { + std::array muXY_curr_triangle = + interpolated_corrected_anisotropic_measure_face( + std::vector {u[i], u[(i + 1) % n], uc}, + std::vector {x[i], x[(i + 1) % n], xc} + ); + + for (std::size_t ix = 0; ix < 3; ix++) + for (std::size_t iy = 0; iy < 3; iy++) + muXY[ix * 3 + iy] += muXY_curr_triangle[ix * 3 + iy]; + } + } + return muXY; +} + +//template +//typename GT::FT triangle_in_ball_ratio(const typename GT::Vector_3 x1, +// const typename GT::Vector_3 x2, +// const typename GT::Vector_3 x3, +// const typename GT::FT r, +// const typename GT::Vector_3 c, +// const std::size_t res = 3) +//{ +// const typename GT::FT R = r * r; +// const typename GT::FT acc = 1.0 / res; +// std::size_t samples_in = 0; +// for (typename GT::FT alpha = acc / 3; alpha < 1; alpha += acc) +// for (typename GT::FT beta = acc / 3; beta < 1 - alpha; beta += acc) +// { +// if ((alpha * x1 + beta * x2 + (1 - alpha - beta) * x3 - c).squared_length() < R) +// samples_in++; +// } +// return samples_in / (typename GT::FT)(res * (res + 1) / 2); +//} + +template +typename GT::FT face_in_ball_ratio(const std::vector& x, + const typename GT::FT r, + const typename GT::Vector_3 c) +{ + const std::size_t n = x.size(); + + // getting center of points + typename GT::Vector_3 xm = + std::accumulate(x.begin(), x.end(), typename GT::Vector_3(0, 0, 0)); + xm /= static_cast(n); + + // computing squared distance of furthest and closest point to ball center + typename GT::FT d_min = (xm - c).squared_length(); + typename GT::FT d_max = d_min; + + for (const typename GT::Vector_3 xi : x) + { + const typename GT::FT d_sq = (xi - c).squared_length(); + d_max = (std::max)(d_sq, d_max); + d_min = (std::min)(d_sq, d_min); + } + + // if the furthest point is inside ball, return 1 + if (d_max <= r * r) return 1.0; + // if the closest point is outside ball, return 0 + else if (r * r <= d_min) return 0.0; + + // else, approximate inclusion ratio of the triangle: + d_max = sqrt(d_max); + d_min = sqrt(d_min); + return (r - d_min) / (d_max - d_min); +} + +template +Principal_curvatures_and_directions principal_curvatures_and_directions_from_anisotropic_measures( + const std::array anisotropic_measure, + const typename GT::FT v_mu0, + const typename GT::Vector_3 u_GT, + const typename GT::FT avg_edge_length +) +{ + // putting anisotropic measure in matrix form + Eigen::Matrix v_muXY = Eigen::Matrix::Zero(); + + for (std::size_t ix = 0; ix < 3; ix++) + for (std::size_t iy = 0; iy < 3; iy++) + v_muXY(ix, iy) = anisotropic_measure[ix * 3 + iy]; + + // constant factor K to force the principal direction eigenvectors to be tangential to the surface + Eigen::Matrix u(u_GT.x(), u_GT.y(), u_GT.z()); + const typename GT::FT K = 1000 * avg_edge_length * v_mu0; + + // symmetrizing and adding the constant term + v_muXY = 0.5 * (v_muXY + v_muXY.transpose()) + K * u * u.transpose(); + + + // computing eigenvalues and eigenvectors + Eigen::SelfAdjointEigenSolver> eigensolver; + + eigensolver.computeDirect(v_muXY); + + if (eigensolver.info() != Eigen::Success) + return Principal_curvatures_and_directions(); + + const Eigen::Matrix eig_vals = eigensolver.eigenvalues(); + const Eigen::Matrix eig_vecs = eigensolver.eigenvectors(); + + const typename GT::Vector_3 min_eig_vec(eig_vecs(0, 1), eig_vecs(1, 1), eig_vecs(2, 1)); + const typename GT::Vector_3 max_eig_vec(eig_vecs(0, 0), eig_vecs(1, 0), eig_vecs(2, 0)); + + // returning principal curvatures and directions (with the correct sign) + return Principal_curvatures_and_directions( + (!is_zero(v_mu0)) ? -eig_vals[1] / v_mu0 : 0.0, + (!is_zero(v_mu0)) ? -eig_vals[0] / v_mu0 : 0.0, + min_eig_vec, + max_eig_vec + ); +} + +// measures are computed for faces only if they are adjacent to the vertex +template +Vertex_measures interpolated_corrected_measures_one_vertex_no_radius( + const PolygonMesh& pmesh, + const typename boost::graph_traits::vertex_descriptor v, + const bool is_mean_curvature_selected, + const bool is_Gaussian_curvature_selected, + const bool is_principal_curvatures_and_directions_selected, + const VPM vpm, + const VNM vnm +) +{ + typedef typename boost::graph_traits::face_descriptor Face_descriptor; + typedef typename boost::graph_traits::vertex_descriptor Vertex_descriptor; + typedef typename GT::Point_3 Point_3; + typedef typename GT::Vector_3 Vector_3; + typedef typename GT::FT FT; + + Vertex_measures vertex_measures; + + std::vector x; + std::vector u; + + // compute for each face around the vertex (except the null (boundary) face) + for (Face_descriptor f : faces_around_target(halfedge(v, pmesh), pmesh)) { + if (f != boost::graph_traits::null_face()) + { + // looping over vertices in face to get point coordinates and normal vectors + for (Vertex_descriptor vi : vertices_around_face(halfedge(f, pmesh), pmesh)) + { + Point_3 pi = get(vpm, vi); + Vector_3 ui = get(vnm, vi); + x.push_back(Vector_3(pi.x(), pi.y(), pi.z())); + u.push_back(ui); + } + + // compute measures for selected curvatures (area is always computed) + vertex_measures.area_measure += interpolated_corrected_area_measure_face(u, x); + + if (is_mean_curvature_selected) + vertex_measures.mean_curvature_measure += interpolated_corrected_mean_curvature_measure_face(u, x); + + if (is_Gaussian_curvature_selected) + vertex_measures.gaussian_curvature_measure += interpolated_corrected_Gaussian_curvature_measure_face(u); + + if (is_principal_curvatures_and_directions_selected) + { + const std::array face_anisotropic_measure = interpolated_corrected_anisotropic_measure_face(u, x); + for (std::size_t i = 0; i < 3 * 3; i++) + vertex_measures.anisotropic_measure[i] += face_anisotropic_measure[i]; + } + } + + x.clear(); + u.clear(); + } + return vertex_measures; +} + +// measures are computed for faces only if they are in the ball of radius 'radius' centered at the vertex +template +Vertex_measures interpolated_corrected_measures_one_vertex( + const PolygonMesh& pmesh, + const typename boost::graph_traits::vertex_descriptor v, + const typename GT::FT radius, + const bool is_mean_curvature_selected, + const bool is_Gaussian_curvature_selected, + const bool is_principal_curvatures_and_directions_selected, + const VPM vpm, + const VNM vnm +) +{ + typedef typename boost::graph_traits::face_descriptor Face_descriptor; + typedef typename boost::graph_traits::vertex_descriptor Vertex_descriptor; + typedef typename GT::Point_3 Point_3; + typedef typename GT::Vector_3 Vector_3; + typedef typename GT::FT FT; + + // the ball expansion is done using a BFS traversal from the vertex + std::queue bfs_queue; + std::unordered_set bfs_visited; + + Vertex_measures vertex_measures; + + Point_3 vp = get(vpm, v); + Vector_3 c = Vector_3(vp.x(), vp.y(), vp.z()); + + for (Face_descriptor f : faces_around_target(halfedge(v, pmesh), pmesh)) { + if (f != boost::graph_traits::null_face()) + { + bfs_queue.push(f); + bfs_visited.insert(f); + } + } + std::vector x; + std::vector u; + while (!bfs_queue.empty()) { + Face_descriptor fi = bfs_queue.front(); + bfs_queue.pop(); + + // looping over vertices in face to get point coordinates and normal vectors + for (Vertex_descriptor vi : vertices_around_face(halfedge(fi, pmesh), pmesh)) + { + Point_3 pi = get(vpm, vi); + Vector_3 ui = get(vnm, vi); + x.push_back(Vector_3(pi.x(), pi.y(), pi.z())); + u.push_back(ui); + } + + // approximate inclusion ratio of the face in the ball + const FT f_ratio = face_in_ball_ratio(x, radius, c); + + // if it is not 0 (not completely outside), compute measures for selected curvatures (area is always computed) + // and add neighboring faces to the bfs queue + if (!is_zero(f_ratio)) + { + vertex_measures.area_measure += f_ratio * interpolated_corrected_area_measure_face(u, x); + + if (is_mean_curvature_selected) + vertex_measures.mean_curvature_measure += f_ratio * interpolated_corrected_mean_curvature_measure_face(u, x); + + if (is_Gaussian_curvature_selected) + vertex_measures.gaussian_curvature_measure += f_ratio * interpolated_corrected_Gaussian_curvature_measure_face(u); + + if (is_principal_curvatures_and_directions_selected) + { + const std::array face_anisotropic_measure = interpolated_corrected_anisotropic_measure_face(u, x); + for (std::size_t i = 0; i < 3 * 3; i++) + vertex_measures.anisotropic_measure[i] += f_ratio * face_anisotropic_measure[i]; + } + + for (Face_descriptor fj : faces_around_face(halfedge(fi, pmesh), pmesh)) + { + if (bfs_visited.find(fj) == bfs_visited.end() && fj != boost::graph_traits::null_face()) + { + bfs_queue.push(fj); + bfs_visited.insert(fj); + } + } + } + + x.clear(); + u.clear(); + } + return vertex_measures; +} + +template +void set_value(const T& value, std::reference_wrapper rw) +{ + rw.get()=value; +} + +template +void set_value(const T&, internal_np::Param_not_found) +{} + +// computes selected curvatures for one specific vertex +template + void interpolated_corrected_curvatures_one_vertex( + const PolygonMesh& pmesh, + const typename boost::graph_traits::vertex_descriptor v, + const NamedParameters& np = parameters::default_values() + ) +{ + typedef typename GetGeomTraits::type GT; + typedef typename GetVertexPointMap::const_type Vertex_position_map; + + typedef dynamic_vertex_property_t Vector_map_tag; + typedef typename boost::property_map::const_type Default_vector_map; + typedef typename internal_np::Lookup_named_param_def::type Vertex_normal_map; + + using parameters::choose_parameter; + using parameters::get_parameter; + using parameters::is_default_parameter; + + Vertex_position_map vpm = choose_parameter(get_parameter(np, CGAL::vertex_point), + get_const_property_map(CGAL::vertex_point, pmesh)); + + Vertex_normal_map vnm = choose_parameter(get_parameter(np, internal_np::vertex_normal_map), + get(Vector_map_tag(), pmesh)); + + // if the normal map is not provided, compute it + if (is_default_parameter::value) + compute_vertex_normals(pmesh, vnm, np); + + // if no radius is given, we pass -1 which will make the expansion be only on the incident faces instead of a ball + 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); + + // if the radius is 0, we use a small epsilon to expand the ball (scaled with the average edge length) + if (is_zero(radius)) + radius = avg_edge_length * 1e-6; + + // determine which curvatures are selected + const bool is_mean_curvature_selected = !is_default_parameter::value; + const bool is_Gaussian_curvature_selected = !is_default_parameter::value; + const bool is_principal_curvatures_and_directions_selected = !is_default_parameter::value; + + Vertex_measures vertex_measures; + + // if the radius is negative, we do not expand the ball (only the incident faces) + if (is_negative(radius)) + { + vertex_measures = interpolated_corrected_measures_one_vertex_no_radius( + pmesh, + v, + is_mean_curvature_selected, + is_Gaussian_curvature_selected, + is_principal_curvatures_and_directions_selected, + vpm, + vnm + ); + } + else + { + vertex_measures = interpolated_corrected_measures_one_vertex( + pmesh, + v, + radius, + is_mean_curvature_selected, + is_Gaussian_curvature_selected, + is_principal_curvatures_and_directions_selected, + vpm, + vnm + ); + } + + // compute the selected curvatures from expanded measures + if (is_mean_curvature_selected) { + set_value(!is_zero(vertex_measures.area_measure) ? 0.5 * vertex_measures.mean_curvature_measure / vertex_measures.area_measure : 0, + get_parameter(np, internal_np::vertex_mean_curvature)); + } + + if (is_Gaussian_curvature_selected) { + set_value(!is_zero(vertex_measures.area_measure) ? vertex_measures.gaussian_curvature_measure / vertex_measures.area_measure : 0, + get_parameter(np, internal_np::vertex_Gaussian_curvature)); + } + + if (is_principal_curvatures_and_directions_selected) { + // compute the principal curvatures and directions from the anisotropic measures + const typename GT::Vector_3 v_normal = get(vnm, v); + const Principal_curvatures_and_directions principal_curvatures_and_directions = principal_curvatures_and_directions_from_anisotropic_measures( + vertex_measures.anisotropic_measure, + vertex_measures.area_measure, + v_normal, + avg_edge_length + ); + set_value(principal_curvatures_and_directions, get_parameter(np, internal_np::vertex_principal_curvatures_and_directions)); + } +} + + +template +class Interpolated_corrected_curvatures_computer +{ + typedef typename GetGeomTraits::type GT; + + typedef typename GT::FT FT; + typedef typename GT::Point_3 Point_3; + typedef typename GT::Vector_3 Vector_3; + + typedef typename boost::graph_traits::halfedge_descriptor Halfedge_descriptor; + typedef typename boost::graph_traits::edge_descriptor Edge_descriptor; + typedef typename boost::graph_traits::face_descriptor Face_descriptor; + typedef typename boost::graph_traits::vertex_descriptor Vertex_descriptor; + + typedef typename GetVertexPointMap::const_type Vertex_position_map; + + typedef dynamic_vertex_property_t Vector_map_tag; + typedef typename boost::property_map::const_type Default_vector_map; + typedef typename internal_np::Lookup_named_param_def::type Vertex_normal_map; + + typedef Constant_property_map Default_scalar_map; + typedef typename internal_np::Lookup_named_param_def::type Vertex_mean_curvature_map; + + typedef typename internal_np::Lookup_named_param_def::type Vertex_Gaussian_curvature_map; + + typedef Constant_property_map> Default_principal_map; + typedef typename internal_np::Lookup_named_param_def::type Vertex_principal_curvatures_and_directions_map; + + typedef typename boost::property_map>::const_type Face_scalar_measure_map; + typedef typename boost::property_map>>::const_type Face_anisotropic_measure_map; + +private: + const PolygonMesh& pmesh; + Vertex_position_map vpm; + Vertex_normal_map vnm; + FT ball_radius; + FT avg_edge_length; + + bool is_mean_curvature_selected; + bool is_Gaussian_curvature_selected; + bool is_principal_curvatures_and_directions_selected; + + Vertex_mean_curvature_map mean_curvature_map; + Vertex_Gaussian_curvature_map gaussian_curvature_map; + Vertex_principal_curvatures_and_directions_map principal_curvatures_and_directions_map; + + Face_scalar_measure_map mu0_map, mu1_map, mu2_map; + Face_anisotropic_measure_map muXY_map; + + void set_property_maps() { + mu0_map = get(CGAL::dynamic_face_property_t(), pmesh); + mu1_map = get(CGAL::dynamic_face_property_t(), pmesh); + mu2_map = get(CGAL::dynamic_face_property_t(), pmesh); + muXY_map = get(CGAL::dynamic_face_property_t>(), pmesh); + + } + + void set_named_params(const NamedParameters& np) + { + using parameters::choose_parameter; + using parameters::get_parameter; + using parameters::is_default_parameter; + + vpm = choose_parameter(get_parameter(np, CGAL::vertex_point), + get_const_property_map(CGAL::vertex_point, pmesh)); + + vnm = choose_parameter(get_parameter(np, internal_np::vertex_normal_map), + get(Vector_map_tag(), pmesh)); + + // if no normal map is given, compute normals + if (is_default_parameter::value) + compute_vertex_normals(pmesh, vnm, np); + + // 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); + set_ball_radius(radius); + + // check which curvature maps are provided by the user (determines which curvatures are computed) + is_mean_curvature_selected = !is_default_parameter::value; + is_Gaussian_curvature_selected = !is_default_parameter::value; + is_principal_curvatures_and_directions_selected = !is_default_parameter::value; + + mean_curvature_map = choose_parameter(get_parameter(np, internal_np::vertex_mean_curvature_map), Default_scalar_map()); + gaussian_curvature_map = choose_parameter(get_parameter(np, internal_np::vertex_Gaussian_curvature_map), Default_scalar_map()); + principal_curvatures_and_directions_map = choose_parameter(get_parameter(np, internal_np::vertex_principal_curvatures_and_directions_map), Default_principal_map()); + } + + void set_ball_radius(const FT radius) { + // if given radius is 0, we use a small epsilon to expand the ball (scaled by the average edge length) + if (is_zero(radius)) + ball_radius = avg_edge_length * 1e-6; + else + ball_radius = radius; + } + +public: + + Interpolated_corrected_curvatures_computer(const PolygonMesh& pmesh, + const NamedParameters& np = parameters::default_values() + ) : + pmesh(pmesh) + { + set_named_params(np); + + if (is_mean_curvature_selected || is_Gaussian_curvature_selected || is_principal_curvatures_and_directions_selected) + { + set_property_maps(); + + compute_selected_curvatures(); + } + } + +private: + + // Computes the (selected) interpolated corrected measures for all faces + // and stores them in the property maps + void interpolated_corrected_selected_measures_all_faces() + { + std::vector x; + std::vector u; + + for (Face_descriptor f : faces(pmesh)) + { + for (Vertex_descriptor v : vertices_around_face(halfedge(f, pmesh), pmesh)) + { + Point_3 p = get(vpm, v); + x.push_back(Vector_3(p.x(), p.y(), p.z())); + u.push_back(get(vnm, v)); + } + put(mu0_map, f, interpolated_corrected_area_measure_face(u, x)); + + if (is_mean_curvature_selected) + put(mu1_map, f, interpolated_corrected_mean_curvature_measure_face(u, x)); + + if (is_Gaussian_curvature_selected) + put(mu2_map, f, interpolated_corrected_Gaussian_curvature_measure_face(u)); + + if (is_principal_curvatures_and_directions_selected) + put(muXY_map, f, interpolated_corrected_anisotropic_measure_face(u, x)); + + x.clear(); + u.clear(); + } + } + + // expand the measures of the faces incident to v + Vertex_measures expand_interpolated_corrected_measure_vertex_no_radius(Vertex_descriptor v) + { + Vertex_measures vertex_measures; + + // add the measures of the faces incident to v (excluding the null (boundary) face) + for (Face_descriptor f : faces_around_target(halfedge(v, pmesh), pmesh)) { + if (f == boost::graph_traits::null_face()) + continue; + + // only add the measures for the selected curvatures (area measure is always added) + vertex_measures.area_measure += get(mu0_map, f); + + if (is_mean_curvature_selected) + vertex_measures.mean_curvature_measure += get(mu1_map, f); + + if (is_Gaussian_curvature_selected) + vertex_measures.gaussian_curvature_measure += get(mu2_map, f); + + if (is_principal_curvatures_and_directions_selected) + { + const std::array face_anisotropic_measure = get(muXY_map, f); + for (std::size_t i = 0; i < 3 * 3; i++) + vertex_measures.anisotropic_measure[i] += face_anisotropic_measure[i]; + } + } + + return vertex_measures; + } + + // expand the measures of the faces inside the ball of radius r around v + Vertex_measures expand_interpolated_corrected_measure_vertex(Vertex_descriptor v) + { + // the ball expansion is done using a BFS traversal from the vertex + std::queue bfs_queue; + std::unordered_set bfs_visited; + + Point_3 vp = get(vpm, v); + Vector_3 c = Vector_3(vp.x(), vp.y(), vp.z()); + + Vertex_measures vertex_measures; + + // add the measures of the faces incident to v (excluding the null (boundary) face) + for (Face_descriptor f : faces_around_target(halfedge(v, pmesh), pmesh)) { + if (f != boost::graph_traits::null_face()) + { + bfs_queue.push(f); + bfs_visited.insert(f); + } + } + while (!bfs_queue.empty()) { + Face_descriptor fi = bfs_queue.front(); + bfs_queue.pop(); + + // looping over vertices in face to get point coordinates + std::vector x; + for (Vertex_descriptor vi : vertices_around_face(halfedge(fi, pmesh), pmesh)) + { + Point_3 pi = get(vpm, vi); + x.push_back(Vector_3(pi.x(), pi.y(), pi.z())); + } + + // compute the inclusion ratio of the face in the ball + const FT f_ratio = face_in_ball_ratio(x, ball_radius, c); + + // if the face is inside the ball, add the measures + // only add the measures for the selected curvatures (area measure is always added) + if (!is_zero(f_ratio)) + { + vertex_measures.area_measure += f_ratio * get(mu0_map, fi); + + if (is_mean_curvature_selected) + vertex_measures.mean_curvature_measure += f_ratio * get(mu1_map, fi); + + if (is_Gaussian_curvature_selected) + vertex_measures.gaussian_curvature_measure += f_ratio * get(mu2_map, fi); + + if (is_principal_curvatures_and_directions_selected) + { + const std::array face_anisotropic_measure = get(muXY_map, fi); + for (std::size_t i = 0; i < 3 * 3; i++) + vertex_measures.anisotropic_measure[i] += f_ratio * face_anisotropic_measure[i]; + } + + for (Face_descriptor fj : faces_around_face(halfedge(fi, pmesh), pmesh)) + { + if (bfs_visited.find(fj) == bfs_visited.end() && fj != boost::graph_traits::null_face()) + { + bfs_queue.push(fj); + bfs_visited.insert(fj); + } + } + } + } + return vertex_measures; + } + + void compute_selected_curvatures() { + interpolated_corrected_selected_measures_all_faces(); + + for (Vertex_descriptor v : vertices(pmesh)) + { + // expand the computed measures (on faces) to the vertices + Vertex_measures vertex_measures = (is_negative(ball_radius)) ? + expand_interpolated_corrected_measure_vertex_no_radius(v) : + expand_interpolated_corrected_measure_vertex(v); + + // compute the selected curvatures from the expanded measures and store them in the property maps + // if the area measure is zero, the curvature is set to zero + if (is_mean_curvature_selected) { + !is_zero(vertex_measures.area_measure) ? + put(mean_curvature_map, v, 0.5 * vertex_measures.mean_curvature_measure / vertex_measures.area_measure) : + put(mean_curvature_map, v, 0); + } + + if (is_Gaussian_curvature_selected) { + !is_zero(vertex_measures.area_measure) ? + put(gaussian_curvature_map, v, vertex_measures.gaussian_curvature_measure / vertex_measures.area_measure) : + put(gaussian_curvature_map, v, 0); + } + + if (is_principal_curvatures_and_directions_selected) { + // compute the principal curvatures and directions from the anisotropic measure + const Vector_3 v_normal = get(vnm, v); + const Principal_curvatures_and_directions principal_curvatures_and_directions = principal_curvatures_and_directions_from_anisotropic_measures( + vertex_measures.anisotropic_measure, + vertex_measures.area_measure, + v_normal, + avg_edge_length + ); + put(principal_curvatures_and_directions_map, v, principal_curvatures_and_directions); + } + } + } +}; + +} // namespace internal + +/** +* \ingroup PMP_corrected_curvatures_grp +* +* computes the interpolated corrected curvatures across the mesh `pmesh`. +* By providing mean, Gaussian and/or principal curvature and direction property maps as named parameters, the user +* can choose which quantites to compute. +* +* @tparam PolygonMesh a model of `FaceListGraph`. +* @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters". +* +* @param pmesh the polygon mesh. +* @param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below. +* `GT` stands for the type of the object provided to the named parameter `geom_traits()`. +* +* \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 `GT::Point_3` as value type.} +* \cgalParamDefault{`boost::get(CGAL::vertex_point, pmesh)`.} +* \cgalParamExtra{If this parameter is omitted, an internal property map for +* `CGAL::vertex_point_t` must be available in `PolygonMesh`.} +* \cgalParamNEnd +* +* \cgalParamNBegin{vertex_normal_map} +* \cgalParamDescription{a property map associating normal vectors to the vertices of `pmesh`.} +* \cgalParamType{a class model of `ReadablePropertyMap` with +* `boost::graph_traits::%vertex_descriptor` +* as key type and `GT::Vector_3` as value type.} +* \cgalParamDefault{`get(dynamic_vertex_property_t(), pmesh)`.} +* \cgalParamExtra{If this parameter is omitted, vertex normals will be +* computed using `compute_vertex_normals()`.} +* \cgalParamNEnd +* +* \cgalParamNBegin{geom_traits} +* \cgalParamDescription{an instance of a geometric traits class.} +* \cgalParamType{a class 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 +* +* \cgalParamNBegin{vertex_mean_curvature_map} +* \cgalParamDescription{a property map associating mean curvatures to the vertices of `pmesh`.} +* \cgalParamType{a class model of `WritablePropertyMap` with +* `boost::graph_traits::%vertex_descriptor` +* as key type and `GT::FT` as value type} +* \cgalParamExtra{If this parameter is omitted, mean curvatures will not be computed.} +* \cgalParamNEnd +* +* \cgalParamNBegin{vertex_Gaussian_curvature_map} +* \cgalParamDescription{a property map associating Gaussian curvatures to the vertices of `pmesh`.} +* \cgalParamType{a class model of `WritablePropertyMap` with +* `boost::graph_traits::%vertex_descriptor` +* as key type and `GT::FT` as value type.} +* \cgalParamExtra{If this parameter is omitted, Gaussian curvatures will not be computed.} +* \cgalParamNEnd +* +* \cgalParamNBegin{vertex_principal_curvatures_and_directions_map} +* \cgalParamDescription{a property map associating principal curvatures and directions to the vertices of `pmesh`.} +* \cgalParamType{a class model of `WritablePropertyMap` with +* `boost::graph_traits::%vertex_descriptor` +* as key type and `Principal_curvatures_and_directions` as value type.} +* \cgalParamExtra{If this parameter is omitted, principal curvatures and directions will not be computed.} +* \cgalParamNEnd +* +* \cgalParamNBegin{ball_radius} +* \cgalParamDescription{a strictly positive scalar value specifying the radius used for expanding curvature measures +* by summing measures of faces inside a ball of this radius centered at the +* vertex expanded from. The summed face measures are weighted by their +* inclusion ratio inside this ball.} +* \cgalParamType{`GT::FT`} +* \cgalParamDefault{`-1`} +* \cgalParamExtra{If this parameter is omitted (`-1`), the expansion is then just a sum of +* measures on faces around the vertex.} +* \cgalParamNEnd +* +* \cgalNamedParamsEnd +* +* @see `interpolated_corrected_mean_curvature()` +* @see `interpolated_corrected_Gaussian_curvature()` +* @see `interpolated_corrected_principal_curvatures_and_directions()` +*/ +template +void interpolated_corrected_curvatures(const PolygonMesh& pmesh, + const NamedParameters& np = parameters::default_values()) +{ + internal::Interpolated_corrected_curvatures_computer(pmesh, np); +} + +/** +* \ingroup PMP_corrected_curvatures_grp +* +* computes the interpolated corrected mean curvature across the mesh `pmesh`. +* +* @tparam PolygonMesh a model of `FaceListGraph`. +* @tparam VertexCurvatureMap a model of `WritablePropertyMap` with +* `boost::graph_traits::%vertex_descriptor` as key type and `GT::FT` as value type. +* @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters". +* +* @param pmesh the polygon mesh. +* @param vcm the vertex property map in which the computed mean curvatures are stored. +* @param np an optional sequence of \ref bgl_namedparameters "Named Parameters". +* `GT` stands for the type of the object provided to the named parameter `geom_traits()`. +* +* \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 `GT::Point_3` as value type.} +* \cgalParamDefault{`boost::get(CGAL::vertex_point, pmesh)`.} +* \cgalParamExtra{If this parameter is omitted, an internal property map for +* `CGAL::vertex_point_t` must be available in `PolygonMesh`.} +* \cgalParamNEnd +* +* \cgalParamNBegin{vertex_normal_map} +* \cgalParamDescription{a property map associating normal vectors to the vertices of `pmesh`.} +* \cgalParamType{a class model of `ReadablePropertyMap` with +* `boost::graph_traits::%vertex_descriptor` +* as key type and `GT::Vector_3` as value type.} +* \cgalParamDefault{`get(dynamic_vertex_property_t(), pmesh)`.} +* \cgalParamExtra{If this parameter is omitted, vertex normals will be +* computed using `compute_vertex_normals()`.} +* \cgalParamNEnd +* +* \cgalParamNBegin{geom_traits} +* \cgalParamDescription{an instance of a geometric traits class.} +* \cgalParamType{a class 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 +* +* \cgalParamNBegin{ball_radius} +* \cgalParamDescription{a scalar value specifying the radius used for expanding curvature measures +* by summing measures of faces inside a ball of this radius centered at the +* vertex expanded from. The summed face measures are weighted by their +* inclusion ratio inside this ball.} +* \cgalParamType{`GT::FT`} +* \cgalParamDefault{`-1`} +* \cgalParamExtra{If this parameter is omitted (`-1`), the expansion will just by a weightless sum of +* measures on faces around the vertex.} +* \cgalParamNEnd +* +* \cgalNamedParamsEnd +* +* @see `interpolated_corrected_Gaussian_curvature()` +* @see `interpolated_corrected_principal_curvatures_and_directions()` +* @see `interpolated_corrected_curvatures()` +*/ + +template +void interpolated_corrected_mean_curvature(const PolygonMesh& pmesh, + VertexCurvatureMap vcm, + const NamedParameters& np = parameters::default_values()) +{ + interpolated_corrected_curvatures(pmesh, np.vertex_mean_curvature_map(vcm)); +} + +/** +* \ingroup PMP_corrected_curvatures_grp +* +* computes the interpolated corrected Gaussian curvature across the mesh `pmesh`. +* +* @tparam PolygonMesh a model of `FaceListGraph`. +* @tparam VertexCurvatureMap a model of `WritablePropertyMap` with +* `boost::graph_traits::%vertex_descriptor` as key type and `GT::FT` as value type. +* @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters". +* +* @param pmesh the polygon mesh. +* @param vcm the vertex property map in which the computed Gaussian curvatures are stored. +* @param np an optional sequence of \ref bgl_namedparameters "Named Parameters". +* `GT` stands for the type of the object provided to the named parameter `geom_traits()`. +* +* \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)`.} +* \cgalParamExtra{If this parameter is omitted, an internal property map for +* `CGAL::vertex_point_t` must be available in `PolygonMesh`.} +* \cgalParamNEnd +* +* \cgalParamNBegin{vertex_normal_map} +* \cgalParamDescription{a property map associating normal vectors to the vertices of `pmesh`.} +* \cgalParamType{a class model of `GT::ReadablePropertyMap` with +* `boost::graph_traits::%vertex_descriptor` +* as key type and `%Vector_3` as value type.} +* \cgalParamDefault{`get(dynamic_vertex_property_t(), pmesh)`.} +* \cgalParamExtra{If this parameter is omitted, vertex normals will be +* computed using `compute_vertex_normals()`.} +* \cgalParamNEnd +* +* \cgalParamNBegin{geom_traits} +* \cgalParamDescription{an instance of a geometric traits class.} +* \cgalParamType{a class 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 +* +* \cgalParamNBegin{ball_radius} +* \cgalParamDescription{a scalar value specifying the radius used for expanding curvature measures +* by summing measures of faces inside a ball of this radius centered at the +* vertex expanded from. The summed face measures are weighted by their +* inclusion ratio inside this ball.} +* \cgalParamType{`GT::FT`} +* \cgalParamDefault{`-1`} +* \cgalParamExtra{If this parameter is omitted (`-1`), the expansion will just by a weightless sum of +* measures on faces around the vertex.} +* \cgalParamNEnd +* +* \cgalNamedParamsEnd +* +* @see `interpolated_corrected_mean_curvature()` +* @see `interpolated_corrected_principal_curvatures_and_directions()` +* @see `interpolated_corrected_curvatures()` +*/ +template +void interpolated_corrected_Gaussian_curvature(const PolygonMesh& pmesh, + VertexCurvatureMap vcm, + const NamedParameters& np = parameters::default_values()) +{ + interpolated_corrected_curvatures(pmesh, np.vertex_Gaussian_curvature_map(vcm)); +} + +/** +* \ingroup PMP_corrected_curvatures_grp +* +* computes the interpolated corrected principal curvatures and directions across the mesh `pmesh`. +* +* @tparam PolygonMesh a model of `FaceListGraph`. +* @tparam VertexCurvatureMap a model of `WritablePropertyMap` with +* `boost::graph_traits::%vertex_descriptor` as key type and +* `Principal_curvatures_and_directions` as value type. +* @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters". +* +* @param pmesh the polygon mesh. +* @param vcm the vertex property map in which the computed principal curvatures and directions are stored. +* @param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below +* `GT` stands for the type of the object provided to the named parameter `geom_traits()`. +* +* \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 `GT::Point_3` as value type.} +* \cgalParamDefault{`boost::get(CGAL::vertex_point, pmesh)`.} +* \cgalParamExtra{If this parameter is omitted, an internal property map for +* `CGAL::vertex_point_t` must be available in `PolygonMesh`.} +* \cgalParamNEnd +* +* \cgalParamNBegin{vertex_normal_map} +* \cgalParamDescription{a property map associating normal vectors to the vertices of `pmesh`.} +* \cgalParamType{a class model of `ReadablePropertyMap` with +* `boost::graph_traits::%vertex_descriptor` +* as key type and `GT::Vector_3` as value type.} +* \cgalParamDefault{`get(dynamic_vertex_property_t(), pmesh)`.} +* \cgalParamExtra{If this parameter is omitted, vertex normals will be +* computed using `compute_vertex_normals()`.} +* \cgalParamNEnd +* +* \cgalParamNBegin{geom_traits} +* \cgalParamDescription{an instance of a geometric traits class.} +* \cgalParamType{a class 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 +* +* \cgalParamNBegin{ball_radius} +* \cgalParamDescription{a scalar value specifying the radius used for expanding curvature measures +* by summing measures of faces inside a ball of this radius centered at the +* vertex expanded from. The summed face measures are weighted by their +* inclusion ratio inside this ball.} +* \cgalParamType{`GT::FT`} +* \cgalParamDefault{`-1`} +* \cgalParamExtra{If this parameter is omitted (`-1`), the expansion will just by a weightless sum of +* measures on faces around the vertex.} +* \cgalParamNEnd +* +* \cgalNamedParamsEnd +* +* @see `interpolated_corrected_mean_curvature()` +* @see `interpolated_corrected_Gaussian_curvature()` +* @see `interpolated_corrected_curvatures()` +*/ +template +void interpolated_corrected_principal_curvatures_and_directions(const PolygonMesh& pmesh, + VertexCurvatureMap vcm, + const NamedParameters& np = parameters::default_values()) +{ + interpolated_corrected_curvatures(pmesh, np.vertex_principal_curvatures_and_directions_map(vcm)); +} + + +/** +* \ingroup PMP_corrected_curvatures_grp +* computes the interpolated corrected curvatures at a vertex `v`. +* By providing mean, Gaussian and/or principal curvature and direction property maps as named parameters, the user +* can choose which quantites to compute. +* +* @tparam PolygonMesh a model of `FaceListGraph`. +* @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters". +* +* @param pmesh the polygon mesh. +* @param v the vertex of `pmesh` to compute the curvatures at. +* @param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below +* `GT` stands for the type of the object provided to the named parameter `geom_traits()`. +* +* \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 `GT::Point_3` as value type.} +* \cgalParamDefault{`boost::get(CGAL::vertex_point, pmesh)`.} +* \cgalParamExtra{If this parameter is omitted, an internal property map for +* `CGAL::vertex_point_t` must be available in `PolygonMesh`.} +* \cgalParamNEnd +* +* \cgalParamNBegin{vertex_normal_map} +* \cgalParamDescription{a property map associating normal vectors to the vertices of `pmesh`.} +* \cgalParamType{a class model of `ReadablePropertyMap` with +* `boost::graph_traits::%vertex_descriptor` +* as key type and `GT::Vector_3` as value type.} +* \cgalParamDefault{`get(dynamic_vertex_property_t(), pmesh)`.} +* \cgalParamExtra{If this parameter is omitted, vertex normals will be +* computed using `compute_vertex_normals()`.} +* \cgalParamNEnd +* +* \cgalParamNBegin{geom_traits} +* \cgalParamDescription{an instance of a geometric traits class.} +* \cgalParamType{a class 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 +* +* \cgalParamNBegin{vertex_mean_curvature} +* \cgalParamDescription{a reference to a scalar value to store the mean curvature at the vertex `v`.} +* \cgalParamType{`std::reference_wrapper`.} +* \cgalParamExtra{If this parameter is omitted, mean curvature will not be computed.} +* \cgalParamNEnd +* +* \cgalParamNBegin{vertex_Gaussian_curvature} +* \cgalParamDescription{a reference to a scalar value to store the Gaussian curvature at the vertex `v`.} +* \cgalParamType{`std::reference_wrapper`.} +* \cgalParamExtra{If this parameter is omitted, Gaussian curvature will not be computed.} +* \cgalParamNEnd +* +* \cgalParamNBegin{vertex_principal_curvatures_and_directions} +* \cgalParamDescription{a reference to an`Principal_curvatures_and_directions` object to store the principal curvatures and directions at the vertex `v`.} +* \cgalParamType{`std::reference_wrapper>`.} +* \cgalParamExtra{If this parameter is omitted, principal curvatures and directions will not be computed.} +* \cgalParamNEnd +* +* \cgalParamNBegin{ball_radius} +* \cgalParamDescription{a scalar value specifying the radius used for expanding curvature measures +* by summing measures of faces inside a ball of this radius centered at the +* vertex expanded from. The summed face measures are weighted by their +* inclusion ratio inside this ball.} +* \cgalParamType{`GT::FT`} +* \cgalParamDefault{`-1`} +* \cgalParamExtra{If this parameter is omitted (`-1`), the epansion is then just a sum of +* measures on faces around the vertex.} +* \cgalParamNEnd +* +* \cgalNamedParamsEnd +* +* @see `CGAL::Polygon_mesh_processing::interpolated_corrected_curvatures()` +* @see `CGAL::Polygon_mesh_processing::interpolated_corrected_mean_curvature_one_vertex()` +* @see `CGAL::Polygon_mesh_processing::interpolated_corrected_Gaussian_curvature_one_vertex()` +* @see `CGAL::Polygon_mesh_processing::interpolated_corrected_principal_curvatures_and_directions_one_vertex()` +*/ +template +void interpolated_corrected_curvatures_one_vertex(const PolygonMesh& pmesh, + typename boost::graph_traits::vertex_descriptor v, + const NamedParameters& np = parameters::default_values()) +{ + internal::interpolated_corrected_curvatures_one_vertex(pmesh, v, np); +} + +/** +* \ingroup PMP_corrected_curvatures_grp +* computes the interpolated corrected mean curvature at vertex `v` of mesh `pmesh`. +* +* @tparam PolygonMesh a model of `FaceListGraph`. +* @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters". +* +* @param pmesh the polygon mesh. +* @param v the vertex of `pmesh` to compute the mean curvature at. +* @param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below +* `GT` stands for the type of the object provided to the named parameter `geom_traits()`. +* +* \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 `GT::Point_3` as value type.} +* \cgalParamDefault{`boost::get(CGAL::vertex_point, pmesh)`.} +* \cgalParamExtra{If this parameter is omitted, an internal property map for +* `CGAL::vertex_point_t` must be available in `PolygonMesh`.} +* \cgalParamNEnd +* +* \cgalParamNBegin{vertex_normal_map} +* \cgalParamDescription{a property map associating normal vectors to the vertices of `pmesh`.} +* \cgalParamType{a class model of `ReadablePropertyMap` with +* `boost::graph_traits::%vertex_descriptor` +* as key type and `GT::Vector_3` as value type.} +* \cgalParamDefault{`get(dynamic_vertex_property_t(), pmesh)`.} +* \cgalParamExtra{If this parameter is omitted, vertex normals will be +* computed using `compute_vertex_normals()`.} +* \cgalParamNEnd +* +* \cgalParamNBegin{geom_traits} +* \cgalParamDescription{an instance of a geometric traits class.} +* \cgalParamType{a class 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 +* +* \cgalParamNBegin{ball_radius} +* \cgalParamDescription{a scalar value specifying the radius used for expanding curvature measures +* by summing measures of faces inside a ball of this radius centered at the +* vertex expanded from. The summed face measures are weighted by their +* inclusion ratio inside this ball.} +* \cgalParamType{`GT::FT`} +* \cgalParamDefault{`-1`} +* \cgalParamExtra{If this parameter is omitted (`-1`), the expansion will just by a weightless sum of +* measures on faces around the vertex.} +* \cgalParamNEnd +* +* \cgalNamedParamsEnd +* +* @return the interpolated corrected mean curvature at the vertex `v`. +* +* @see `interpolated_corrected_mean_curvature()` +* @see `interpolated_corrected_Gaussian_curvature_one_vertex()` +* @see `interpolated_corrected_principal_curvatures_and_directions_one_vertex()` +* @see `interpolated_corrected_curvatures_one_vertex()` +*/ + +template +#ifdef DOXYGEN_RUNNING +typename GT::FT +#else +typename GetGeomTraits::type::FT +#endif +interpolated_corrected_mean_curvature_one_vertex(const PolygonMesh& pmesh, + typename boost::graph_traits::vertex_descriptor v, + const NamedParameters& np = parameters::default_values()) +{ + typename GetGeomTraits::type::FT mean_curvature; + interpolated_corrected_curvatures_one_vertex(pmesh, v, np.vertex_mean_curvature(std::ref(mean_curvature))); + return mean_curvature; +} + +/** +* \ingroup PMP_corrected_curvatures_grp +* computes the interpolated corrected Gaussian curvature at vertex `v` of mesh `pmesh`. +* +* @tparam PolygonMesh a model of `FaceListGraph`. +* @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters". +* +* @param pmesh the polygon mesh. +* @param v the vertex of `pmesh` to compute the Gaussian curvature at. +* @param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below +* `GT` stands for the type of the object provided to the named parameter `geom_traits()`. +* +* \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 `GT::Point_3` as value type.} +* \cgalParamDefault{`boost::get(CGAL::vertex_point, pmesh)`.} +* \cgalParamExtra{If this parameter is omitted, an internal property map for +* `CGAL::vertex_point_t` must be available in `PolygonMesh`.} +* \cgalParamNEnd +* +* \cgalParamNBegin{vertex_normal_map} +* \cgalParamDescription{a property map associating normal vectors to the vertices of `pmesh`.} +* \cgalParamType{a class model of `ReadablePropertyMap` with +* `boost::graph_traits::%vertex_descriptor` +* as key type and `GT::Vector_3` as value type.} +* \cgalParamDefault{`get(dynamic_vertex_property_t(), pmesh)`.} +* \cgalParamExtra{If this parameter is omitted, vertex normals will be +* computed using `compute_vertex_normals()`.} +* \cgalParamNEnd +* +* \cgalParamNBegin{geom_traits} +* \cgalParamDescription{an instance of a geometric traits class.} +* \cgalParamType{a class 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 +* +* \cgalParamNBegin{ball_radius} +* \cgalParamDescription{a scalar value specifying the radius used for expanding curvature measures +* by summing measures of faces inside a ball of this radius centered at the +* vertex expanded from. The summed face measures are weighted by their +* inclusion ratio inside this ball.} +* \cgalParamType{`GT::FT`} +* \cgalParamDefault{`-1`} +* \cgalParamExtra{If this parameter is omitted (`-1`), the expansion will just by a weightless sum of +* measures on faces around the vertex.} +* \cgalParamNEnd +* +* \cgalNamedParamsEnd +* +* @return the interpolated corrected Gaussian curvature at the vertex `v`. +* +* @see `interpolated_corrected_Gaussian_curvature()` +* @see `interpolated_corrected_mean_curvature_one_vertex()` +* @see `interpolated_corrected_principal_curvatures_and_directions_one_vertex()` +* @see `interpolated_corrected_curvatures_one_vertex()` +*/ + +template +#ifdef DOXYGEN_RUNNING +typename GT::FT +#else +typename GetGeomTraits::type::FT +#endif +interpolated_corrected_Gaussian_curvature_one_vertex(const PolygonMesh& pmesh, + typename boost::graph_traits::vertex_descriptor v, + const NamedParameters& np = parameters::default_values()) +{ + typename GetGeomTraits::type::FT gc; + interpolated_corrected_curvatures_one_vertex(pmesh, v, np.vertex_Gaussian_curvature(std::ref(gc))); + return gc; +} + +/** +* \ingroup PMP_corrected_curvatures_grp +* computes the interpolated corrected principal curvatures and directions at vertex `v` of mesh `pmesh`. +* +* @tparam PolygonMesh a model of `FaceListGraph`. +* @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters". +* +* @param pmesh the polygon mesh. +* @param v the vertex of `pmesh` to compute the principal curvatures and directions at. +* @param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below +* `GT` stands for the type of the object provided to the named parameter `geom_traits()`. +* +* \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)`.} +* \cgalParamExtra{If this parameter is omitted, an internal property map for +* `CGAL::vertex_point_t` must be available in `PolygonMesh`.} +* \cgalParamNEnd +* +* \cgalParamNBegin{vertex_normal_map} +* \cgalParamDescription{a property map associating normal vectors to the vertices of `pmesh`.} +* \cgalParamType{a class model of `ReadablePropertyMap` with +* `boost::graph_traits::%vertex_descriptor` +* as key type and `%Vector_3` as value type.} +* \cgalParamDefault{`get(dynamic_vertex_property_t(), pmesh)`.} +* \cgalParamExtra{If this parameter is omitted, vertex normals will be +* computed using `compute_vertex_normals()`.} +* \cgalParamNEnd +* +* \cgalParamNBegin{geom_traits} +* \cgalParamDescription{an instance of a geometric traits class.} +* \cgalParamType{a class 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 +* +* \cgalParamNBegin{ball_radius} +* \cgalParamDescription{a scalar value specifying the radius used for expanding curvature measures +* by summing measures of faces inside a ball of this radius centered at the +* vertex expanded from. The summed face measures are weighted by their +* inclusion ratio inside this ball.} +* \cgalParamType{`GT::FT`} +* \cgalParamDefault{`-1`} +* \cgalParamExtra{If this parameter is omitted (`-1`), the expansion will just by a weightless sum of +* measures on faces around the vertex.} +* \cgalParamNEnd +* \cgalNamedParamsEnd +* +* @return the interpolated corrected principal curvatures and directions at the vertex `v`. +* +* @see `interpolated_corrected_principal_curvatures_and_directions()` +* @see `interpolated_corrected_mean_curvature_one_vertex()` +* @see `interpolated_corrected_Gaussian_curvature_one_vertex()` +* @see `interpolated_corrected_curvatures_one_vertex()` +*/ + +template +#ifdef DOXYGEN_RUNNING +Principal_curvatures_and_directions +#else +Principal_curvatures_and_directions::type> +#endif +interpolated_corrected_principal_curvatures_and_directions_one_vertex(const PolygonMesh& pmesh, + typename boost::graph_traits::vertex_descriptor v, + const NamedParameters& np = parameters::default_values()) +{ + using GT=typename GetGeomTraits::type; + Principal_curvatures_and_directions pcd; + interpolated_corrected_curvatures_one_vertex(pmesh, v, np.vertex_principal_curvatures_and_directions(std::ref(pcd))); + return pcd; +} + +} // namespace Polygon_mesh_processing +} // namespace CGAL + +#endif // CGAL_POLYGON_MESH_PROCESSING_INTERPOLATED_CORRECTED_CURVATURES_H diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt b/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt index eb9b98cbe59..8ba77a7ff15 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt @@ -80,6 +80,8 @@ if(TARGET CGAL::Eigen3_support) target_link_libraries(test_shape_smoothing PUBLIC CGAL::Eigen3_support) create_single_source_cgal_program("delaunay_remeshing_test.cpp") target_link_libraries(delaunay_remeshing_test PUBLIC CGAL::Eigen3_support) + create_single_source_cgal_program("test_interpolated_corrected_curvatures.cpp") + target_link_libraries(test_interpolated_corrected_curvatures PUBLIC CGAL::Eigen3_support) create_single_source_cgal_program("test_decimation_of_planar_patches.cpp") target_link_libraries(test_decimation_of_planar_patches PUBLIC CGAL::Eigen3_support) else() 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 new file mode 100644 index 00000000000..2f1af8adb10 --- /dev/null +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_interpolated_corrected_curvatures.cpp @@ -0,0 +1,255 @@ +#include +#include +#include +#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 Epic_kernel; +typedef CGAL::Surface_mesh SMesh; +typedef CGAL::Polyhedron_3 Polyhedron; + +struct Average_test_info { + Epic_kernel::FT expansion_radius = -1; + Epic_kernel::FT mean_curvature_avg; + Epic_kernel::FT gaussian_curvature_avg; + Epic_kernel::FT principal_curvature_avg; + Epic_kernel::FT tolerance = 0.9; + + Average_test_info( + Epic_kernel::FT mean_curvature_avg, + Epic_kernel::FT gaussian_curvature_avg, + Epic_kernel::FT principal_curvature_avg, + Epic_kernel::FT expansion_radius = -1, + Epic_kernel::FT tolerance = 0.9 + ) : + expansion_radius(expansion_radius), + mean_curvature_avg(mean_curvature_avg), + gaussian_curvature_avg(gaussian_curvature_avg), + principal_curvature_avg(principal_curvature_avg), + tolerance(tolerance) + { + } + +}; + +bool passes_comparison(Epic_kernel::FT result, Epic_kernel::FT expected, Epic_kernel::FT tolerance) { + 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_average_curvatures(std::string mesh_path, + Average_test_info test_info, + bool compare_single_vertex = false, + int scale_factor_exponent = 0 +) { + 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; + } + + // The following part is used to scale the given mesh and expected curvatures by a constant factor + // this is used to test the stability of the implementation across different scales + if (scale_factor_exponent) { + Epic_kernel::FT factor = pow(10, scale_factor_exponent); + + test_info.expansion_radius *= factor; + test_info.mean_curvature_avg /= factor; + test_info.gaussian_curvature_avg /= factor * factor; + test_info.principal_curvature_avg /= factor; + + auto vpm = get(CGAL::vertex_point, pmesh); + + for (vertex_descriptor vi : vertices(pmesh)) { + Epic_kernel::Point_3 pi = get(vpm, vi); + Epic_kernel::Point_3 pi_new(pi.x() * factor, pi.y() * factor, pi.z() * factor); + put(vpm, vi, pi_new); + } + } + + + typename boost::property_map>::type + mean_curvature_map = get(CGAL::dynamic_vertex_property_t(), pmesh), + gaussian_curvature_map = get(CGAL::dynamic_vertex_property_t(), pmesh); + typename boost::property_map + >>::type + principal_curvatures_and_directions_map = + get(CGAL::dynamic_vertex_property_t>(), pmesh); + // test_info.expansion_radius -> test if no radius is provided by user. + if (test_info.expansion_radius < 0) { + PMP::interpolated_corrected_mean_curvature(pmesh, mean_curvature_map); + PMP::interpolated_corrected_Gaussian_curvature(pmesh, gaussian_curvature_map); + PMP::interpolated_corrected_principal_curvatures_and_directions(pmesh, principal_curvatures_and_directions_map); + } + else { + PMP::interpolated_corrected_mean_curvature( + pmesh, + mean_curvature_map, + CGAL::parameters::ball_radius(test_info.expansion_radius) + ); + + PMP::interpolated_corrected_Gaussian_curvature( + pmesh, + gaussian_curvature_map, + CGAL::parameters::ball_radius(test_info.expansion_radius) + ); + + PMP::interpolated_corrected_principal_curvatures_and_directions( + pmesh, + principal_curvatures_and_directions_map, + CGAL::parameters::ball_radius(test_info.expansion_radius) + ); + } + + Epic_kernel::FT mean_curvature_avg = 0, gaussian_curvature_avg = 0, principal_curvature_avg = 0; + + for (vertex_descriptor v : vertices(pmesh)) { + mean_curvature_avg += get(mean_curvature_map, v); + gaussian_curvature_avg += get(gaussian_curvature_map, v); + principal_curvature_avg += get(principal_curvatures_and_directions_map, v).min_curvature + + get(principal_curvatures_and_directions_map, v).max_curvature; + } + + mean_curvature_avg /= vertices(pmesh).size(); + gaussian_curvature_avg /= vertices(pmesh).size(); + principal_curvature_avg /= vertices(pmesh).size() * 2; + + // are average curvatures equal to expected? + assert(passes_comparison(mean_curvature_avg, test_info.mean_curvature_avg, test_info.tolerance)); + assert(passes_comparison(gaussian_curvature_avg, test_info.gaussian_curvature_avg, test_info.tolerance)); + assert(passes_comparison(principal_curvature_avg, test_info.principal_curvature_avg, test_info.tolerance)); + + // computing curvatures together from interpolated_corrected_curvatures() + PMP::interpolated_corrected_curvatures( + pmesh, + CGAL::parameters::ball_radius(test_info.expansion_radius) + .vertex_mean_curvature_map(mean_curvature_map) + .vertex_Gaussian_curvature_map(gaussian_curvature_map) + .vertex_principal_curvatures_and_directions_map(principal_curvatures_and_directions_map) + ); + + Epic_kernel::FT new_mean_curvature_avg = 0, new_Gaussian_curvature_avg = 0, new_principal_curvature_avg = 0; + + for (vertex_descriptor v : vertices(pmesh)) { + new_mean_curvature_avg += get(mean_curvature_map, v); + new_Gaussian_curvature_avg += get(gaussian_curvature_map, v); + new_principal_curvature_avg += get(principal_curvatures_and_directions_map, v).min_curvature + + get(principal_curvatures_and_directions_map, v).max_curvature; + } + + new_mean_curvature_avg /= vertices(pmesh).size(); + new_Gaussian_curvature_avg /= vertices(pmesh).size(); + new_principal_curvature_avg /= vertices(pmesh).size() * 2; + + // are average curvatures computed from interpolated_corrected_curvatures() + // equal to average curvatures each computed on its own? + assert(passes_comparison(mean_curvature_avg, new_mean_curvature_avg, 0.99)); + assert(passes_comparison(gaussian_curvature_avg, new_Gaussian_curvature_avg, 0.99)); + assert(passes_comparison(principal_curvature_avg, new_principal_curvature_avg, 0.99)); + + if (compare_single_vertex) { + // computing curvatures together from interpolated_corrected_curvatures() + + Epic_kernel::FT single_vertex_mean_curvature_avg = 0, + single_vertex_Gaussian_curvature_avg = 0, + single_vertex_principal_curvature_avg = 0; + + Epic_kernel::FT h, g; + PMP::Principal_curvatures_and_directions p; + + for (vertex_descriptor v : vertices(pmesh)) { + PMP::interpolated_corrected_curvatures_one_vertex( + pmesh, + v, + CGAL::parameters::vertex_Gaussian_curvature(std::ref(g)) + .vertex_mean_curvature(std::ref(h)) + .vertex_principal_curvatures_and_directions(std::ref(p)) + .ball_radius(test_info.expansion_radius) + ); + + single_vertex_mean_curvature_avg += h; + single_vertex_Gaussian_curvature_avg += g; + single_vertex_principal_curvature_avg += p.min_curvature + p.max_curvature; + } + + single_vertex_mean_curvature_avg /= vertices(pmesh).size(); + single_vertex_Gaussian_curvature_avg /= vertices(pmesh).size(); + single_vertex_principal_curvature_avg /= vertices(pmesh).size() * 2; + + assert(passes_comparison(mean_curvature_avg, single_vertex_mean_curvature_avg, 0.99)); + assert(passes_comparison(gaussian_curvature_avg, single_vertex_Gaussian_curvature_avg, 0.99)); + assert(passes_comparison(principal_curvature_avg, single_vertex_principal_curvature_avg, 0.99)); + } + +} + +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 + // 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); + test_average_curvatures("meshes/sphere.off", Average_test_info(2, 4, 2), true); + + // Same mesh but with specified expansion radii of 0 and 0.25 (half radius of sphere) + test_average_curvatures("meshes/sphere.off", Average_test_info(2, 4, 2, 0), true); + test_average_curvatures("meshes/sphere.off", Average_test_info(2, 4, 2, 0.25), true); + + // testing on a simple sphere(r = 10), on both Polyhedron & SurfaceMesh: + // Expected: Mean Curvature = 0.1, Gaussian Curvature = 0.01, Principal Curvatures = 0.1 & 0.1 so 0.1 on avg. + test_average_curvatures("meshes/sphere966.off", Average_test_info(0.1, 0.01, 0.1)); + test_average_curvatures("meshes/sphere966.off", Average_test_info(0.1, 0.01, 0.1)); + + // Same mesh but with specified expansion radii of 0 and 5 (half radius of sphere) + test_average_curvatures("meshes/sphere966.off", Average_test_info(0.1, 0.01, 0.1, 0)); + test_average_curvatures("meshes/sphere966.off", Average_test_info(0.1, 0.01, 0.1, 5)); + + // testing on a simple half cylinder(r = 1), on both Polyhedron & SurfaceMesh: + // Expected: Mean Curvature = 0.5, Gaussian Curvature = 0, Principal Curvatures = 0 & 1 so 0.5 on avg. + test_average_curvatures("meshes/cylinder.off", Average_test_info(0.5, 0, 0.5)); + test_average_curvatures("meshes/cylinder.off", Average_test_info(0.5, 0, 0.5)); + + // Same mesh but with specified expansion radii of 0 and 0.5 (half radius of cylinder) + test_average_curvatures("meshes/cylinder.off", Average_test_info(0.5, 0, 0.5, 0)); + test_average_curvatures("meshes/cylinder.off", Average_test_info(0.5, 0, 0.5, 0.5)); + + // Same tests as last one, but with a scaling on the mesh with different values to check for scale stability + test_average_curvatures("meshes/cylinder.off", Average_test_info(0.5, 0, 0.5, 0), false, -6); + test_average_curvatures("meshes/cylinder.off", Average_test_info(0.5, 0, 0.5, 0.5), false, -6); + + test_average_curvatures("meshes/cylinder.off", Average_test_info(0.5, 0, 0.5, 0), false, -3); + test_average_curvatures("meshes/cylinder.off", Average_test_info(0.5, 0, 0.5, 0.5), false, -3); + + test_average_curvatures("meshes/cylinder.off", Average_test_info(0.5, 0, 0.5, 0), false, -1); + test_average_curvatures("meshes/cylinder.off", Average_test_info(0.5, 0, 0.5, 0.5), false, -1); + + test_average_curvatures("meshes/cylinder.off", Average_test_info(0.5, 0, 0.5, 0), false, 1); + test_average_curvatures("meshes/cylinder.off", Average_test_info(0.5, 0, 0.5, 0.5), false, 1); + + test_average_curvatures("meshes/cylinder.off", Average_test_info(0.5, 0, 0.5, 0), false, 3); + test_average_curvatures("meshes/cylinder.off", Average_test_info(0.5, 0, 0.5, 0.5), false, 3); + + test_average_curvatures("meshes/cylinder.off", Average_test_info(0.5, 0, 0.5, 0), false, 6); + test_average_curvatures("meshes/cylinder.off", Average_test_info(0.5, 0, 0.5, 0.5), false, 6); +} diff --git a/Polyhedron/demo/Polyhedron/Plugins/Display/Display_property.ui b/Polyhedron/demo/Polyhedron/Plugins/Display/Display_property.ui index 7a5a19621f4..70c22f8348c 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/Display/Display_property.ui +++ b/Polyhedron/demo/Polyhedron/Plugins/Display/Display_property.ui @@ -84,18 +84,18 @@ Ramp Colors - + - Color Min... + Min - Color Max... + Max @@ -135,7 +135,7 @@ Zoom - + @@ -158,7 +158,7 @@ Ramp Extrema - + @@ -186,6 +186,32 @@ + + + + 0 + + + 100 + + + true + + + Qt::Horizontal + + + QSlider::TicksAbove + + + + + + + Expanding Radius: 0 + + + diff --git a/Polyhedron/demo/Polyhedron/Plugins/Display/Display_property_plugin.cpp b/Polyhedron/demo/Polyhedron/Plugins/Display/Display_property_plugin.cpp index af1eef9e115..050487c9fec 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/Display/Display_property_plugin.cpp +++ b/Polyhedron/demo/Polyhedron/Plugins/Display/Display_property_plugin.cpp @@ -8,6 +8,7 @@ #include #include #include +#include #include #include #include @@ -15,6 +16,7 @@ #include #include #include +#include #include "Scene_points_with_normal_item.h" @@ -31,6 +33,7 @@ #include #include #include +#include #define ARBITRARY_DBL_MIN 1.0E-30 #define ARBITRARY_DBL_MAX 1.0E+30 @@ -41,6 +44,8 @@ typedef CGAL::Three::Triangle_container Tri; typedef CGAL::Three::Viewer_interface VI; +namespace PMP = CGAL::Polygon_mesh_processing; + class Scene_heat_item : public CGAL::Three::Scene_item_rendering_helper { @@ -333,7 +338,10 @@ class DisplayPropertyPlugin : typedef CGAL::Heat_method_3::Surface_mesh_geodesic_distances_3 Heat_method_idt; typedef CGAL::dynamic_vertex_property_t Vertex_source_tag; typedef boost::property_map::type Vertex_source_map; - + enum CurvatureType { + MEAN_CURVATURE, + GAUSSIAN_CURVATURE, +}; public: bool applicable(QAction* action) const Q_DECL_OVERRIDE @@ -482,8 +490,10 @@ public: connect(scene_obj, SIGNAL(itemIndexSelected(int)), this,SLOT(detectScalarProperties(int))); - on_propertyBox_currentIndexChanged(0); + connect(dock_widget->expandingRadiusSlider, SIGNAL(valueChanged(int)), + this, SLOT(setExpandingRadius(int))); + on_propertyBox_currentIndexChanged(0); } private: @@ -529,6 +539,8 @@ private Q_SLOTS: dock_widget->propertyBox->addItem("Scaled Jacobian"); dock_widget->propertyBox->addItem("Heat Intensity"); dock_widget->propertyBox->addItem("Heat Intensity (Intrinsic Delaunay)"); + dock_widget->propertyBox->addItem("Interpolated Corrected Mean Curvature"); + dock_widget->propertyBox->addItem("Interpolated Corrected Gaussian Curvature"); detectSMScalarProperties(sm_item->face_graph()); } @@ -603,6 +615,14 @@ private Q_SLOTS: return; sm_item->setRenderingMode(Gouraud); break; + case 4: // Interpolated Corrected Mean Curvature + displayInterpolatedCurvatureMeasure(sm_item, MEAN_CURVATURE); + sm_item->setRenderingMode(Gouraud); + break; + case 5: // Interpolated Corrected Gaussian Curvature + displayInterpolatedCurvatureMeasure(sm_item, GAUSSIAN_CURVATURE); + sm_item->setRenderingMode(Gouraud); + break; default: if(dock_widget->propertyBox->currentText().contains("v:")) { @@ -637,6 +657,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:interpolated_corrected_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:interpolated_corrected_Gaussian_curvature"); + if (does_exist) + sm_item->face_graph()->remove_property_map(pmap); }); QApplication::restoreOverrideCursor(); sm_item->invalidateOpenGLBuffers(); @@ -668,13 +696,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(mean_curvature_min.count(sm_item) > 0); + dock_widget->zoomToMaxButton->setEnabled(mean_curvature_max.count(sm_item) > 0); + break; + case 5: + dock_widget->zoomToMinButton->setEnabled(gaussian_curvature_min.count(sm_item) > 0); + dock_widget->zoomToMaxButton->setEnabled(gaussian_curvature_max.count(sm_item) > 0); + break; default: break; } @@ -685,6 +721,10 @@ private Q_SLOTS: { Scene_surface_mesh_item* item = qobject_cast(sender()); + + maxEdgeLength = -1; + setExpandingRadius(dock_widget->expandingRadiusSlider->value()); + if(!item) return; SMesh& smesh = *item->face_graph(); @@ -701,11 +741,22 @@ private Q_SLOTS: { smesh.remove_property_map(angles); } + SMesh::Property_map mean_curvature; + std::tie(mean_curvature, found) = smesh.property_map("v:interpolated_corrected_mean_curvature"); + if (found) + { + smesh.remove_property_map(mean_curvature); + } + SMesh::Property_map gaussian_curvature; + std::tie(gaussian_curvature, found) = smesh.property_map("v:interpolated_corrected_Gaussian_curvature"); + if (found) + { + smesh.remove_property_map(gaussian_curvature); + } } void displayScaledJacobian(Scene_surface_mesh_item* item) { - SMesh& smesh = *item->face_graph(); //compute and store the jacobian per face bool non_init; @@ -742,16 +793,118 @@ private Q_SLOTS: treat_sm_property("f:jacobian", item->face_graph()); } - bool resetScaledJacobian(Scene_surface_mesh_item* item) + void setExpandingRadius(int val_int) { - SMesh& smesh = *item->face_graph(); - if(!smesh.property_map("f:jacobian").second) + double sliderMin = dock_widget->expandingRadiusSlider->minimum(); + double sliderMax = dock_widget->expandingRadiusSlider->maximum() - sliderMin; + double val = val_int - sliderMin; + sliderMin = 0; + + SMesh& smesh = *(qobject_cast(scene->item(scene->mainSelectionIndex())))->face_graph(); + + auto vpm = get(CGAL::vertex_point, smesh); + + if (maxEdgeLength < 0) { - return false; + auto edge_range = CGAL::edges(smesh); + + if (edge_range.begin() == edge_range.end()) + { + expand_radius = 0; + dock_widget->expandingRadiusLabel->setText(tr("Expanding Radius : %1").arg(expand_radius)); + return; + } + + auto edge_reference = std::max_element(edge_range.begin(), edge_range.end(), [&, vpm, smesh](auto l, auto r) { + auto res = EPICK().compare_squared_distance_3_object()( + get(vpm, source((l), smesh)), + get(vpm, target((l), smesh)), + get(vpm, source((r), smesh)), + get(vpm, target((r), smesh))); + return res == CGAL::SMALLER; + }); + + // if edge_reference is not derefrenceble + if (edge_reference == edge_range.end()) + { + expand_radius = 0; + dock_widget->expandingRadiusLabel->setText(tr("Expanding Radius : %1").arg(expand_radius)); + return; + } + + maxEdgeLength = sqrt( + (get(vpm, source((*edge_reference), smesh)) - get(vpm, target((*edge_reference), smesh))) + .squared_length() + ); + } - dock_widget->minBox->setValue(jacobian_min[item].first-0.01); - dock_widget->maxBox->setValue(jacobian_max[item].first); - return true; + + double outMax = 5 * maxEdgeLength, base = 1.2; + + expand_radius = (pow(base, val) - 1) * outMax / (pow(base, sliderMax) - 1); + dock_widget->expandingRadiusLabel->setText(tr("Expanding Radius : %1").arg(expand_radius)); + + } + + void displayInterpolatedCurvatureMeasure(Scene_surface_mesh_item* item, CurvatureType mu_index) + { + if (mu_index != MEAN_CURVATURE && mu_index != GAUSSIAN_CURVATURE) return; + + std::string tied_string = (mu_index == MEAN_CURVATURE)? + "v:interpolated_corrected_mean_curvature": "v:interpolated_corrected_Gaussian_curvature"; + SMesh& smesh = *item->face_graph(); + + const auto vnm = smesh.property_map("v:normal_before_perturbation").first; + const bool vnm_exists = smesh.property_map("v:normal_before_perturbation").second; + + //compute once and store the value per vertex + 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 (vnm_exists) { + if (mu_index == MEAN_CURVATURE) + PMP::interpolated_corrected_mean_curvature(smesh, mu_i_map, CGAL::parameters::ball_radius(expand_radius).vertex_normal_map(vnm)); + else + PMP::interpolated_corrected_Gaussian_curvature(smesh, mu_i_map, CGAL::parameters::ball_radius(expand_radius).vertex_normal_map(vnm)); + } + else { + if (mu_index == MEAN_CURVATURE) + PMP::interpolated_corrected_mean_curvature(smesh, mu_i_map, CGAL::parameters::ball_radius(expand_radius)); + else + PMP::interpolated_corrected_Gaussian_curvature(smesh, mu_i_map, CGAL::parameters::ball_radius(expand_radius)); + } + 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 (mu_i_map[v] > res_max) + { + res_max = mu_i_map[v]; + max_index = v; + } + if (mu_i_map[v] < res_min) + { + res_min = mu_i_map[v]; + min_index = v; + } + } + if (mu_index == MEAN_CURVATURE){ + mean_curvature_max[item] = std::make_pair(res_max, max_index); + mean_curvature_min[item] = std::make_pair(res_min, min_index); + } + else { + gaussian_curvature_max[item] = std::make_pair(res_max, max_index); + gaussian_curvature_min[item] = std::make_pair(res_min, min_index); + } + + connect(item, &Scene_surface_mesh_item::itemChanged, + this, &DisplayPropertyPlugin::resetProperty); + } + treat_sm_property(tied_string, item->face_graph()); } @@ -1054,6 +1207,8 @@ private Q_SLOTS: break; } case 1: + case 4: + case 5: dock_widget->groupBox-> setEnabled(true); dock_widget->groupBox_3->setEnabled(true); @@ -1113,6 +1268,26 @@ private Q_SLOTS: dummy_fd, dummy_p); } + break; + case 4: + { + ::zoomToId(*item->face_graph(), + QString("v%1").arg(mean_curvature_min[item].second), + getActiveViewer(), + dummy_fd, + dummy_p); + } + break; + case 5: + { + ::zoomToId(*item->face_graph(), + QString("v%1").arg(gaussian_curvature_min[item].second), + getActiveViewer(), + dummy_fd, + dummy_p); + } + break; + break; break; default: break; @@ -1146,7 +1321,25 @@ private Q_SLOTS: dummy_fd, dummy_p); } - break; + break; + case 4: + { + ::zoomToId(*item->face_graph(), + QString("v%1").arg(mean_curvature_max[item].second), + getActiveViewer(), + dummy_fd, + dummy_p); + } + break; + case 5: + { + ::zoomToId(*item->face_graph(), + QString("v%1").arg(gaussian_curvature_max[item].second), + getActiveViewer(), + dummy_fd, + dummy_p); + } + break; default: break; } @@ -1315,10 +1508,10 @@ private: void displayMapLegend(const std::vector& values) { // Create a legend_ and display it - const std::size_t size = (std::min)(color_map.size(), (std::size_t)256); + const std::size_t size = (std::min)(color_map.size(), (std::size_t)2048); const int text_height = 20; const int height = text_height*static_cast(size) + text_height; - const int width = 140; + const int width = 170; const int cell_width = width/3; const int top_margin = 15; const int left_margin = 5; @@ -1344,8 +1537,8 @@ 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, tr("%1").arg(values[i], 0, 'f', 3, QLatin1Char(' '))); + 100, text_height); + painter.drawText(text_rect, Qt::AlignCenter, tr("%1").arg(values[i], 0, 'f', 7, QLatin1Char(' '))); } if(color_map.size() > size){ QRect text_rect(left_margin + cell_width+10, 0, @@ -1364,7 +1557,7 @@ private: { // Create a legend_ and display it const int height = 256; - const int width = 140; + const int width = 170; const int cell_width = width/3; const int top_margin = 5; const int left_margin = 5; @@ -1408,11 +1601,11 @@ private: painter.setPen(Qt::blue); QRect min_text_rect(left_margin + cell_width+10,drawing_height - top_margin, 100, text_height); - painter.drawText(min_text_rect, Qt::AlignCenter, tr("%1").arg(min_value, 0, 'f', 1)); + painter.drawText(min_text_rect, Qt::AlignCenter, tr("%1").arg(min_value, 0, 'f', 7)); QRect max_text_rect(left_margin + cell_width+10, drawing_height - top_margin - 200, 100, text_height); - painter.drawText(max_text_rect, Qt::AlignCenter, tr("%1").arg(max_value, 0, 'f', 1)); + painter.drawText(max_text_rect, Qt::AlignCenter, tr("%1").arg(max_value, 0, 'f', 7)); dock_widget->legendLabel->setPixmap(legend_); } @@ -1436,9 +1629,17 @@ private: std::unordered_map > angles_min; std::unordered_map > angles_max; + + std::unordered_map > mean_curvature_min; + std::unordered_map > mean_curvature_max; + + std::unordered_map > gaussian_curvature_min; + std::unordered_map > gaussian_curvature_max; + std::unordered_map is_source; - + double expand_radius = 0; + double maxEdgeLength = -1; double minBox; double maxBox; QPixmap legend_; @@ -1718,7 +1919,7 @@ private: } - EPICK::Vector_3 unit_center_normal = CGAL::Polygon_mesh_processing::compute_face_normal(f, mesh); + EPICK::Vector_3 unit_center_normal = PMP::compute_face_normal(f, mesh); unit_center_normal *= 1.0/CGAL::approximate_sqrt(unit_center_normal.squared_length()); for(std::size_t i = 0; i < corner_areas.size(); ++i) diff --git a/Polyhedron/demo/Polyhedron/Plugins/PMP/CMakeLists.txt b/Polyhedron/demo/Polyhedron/Plugins/PMP/CMakeLists.txt index 152d84952ae..6d5b6cf16d4 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/PMP/CMakeLists.txt +++ b/Polyhedron/demo/Polyhedron/Plugins/PMP/CMakeLists.txt @@ -9,6 +9,20 @@ else() message(STATUS "NOTICE: Eigen 3.1 (or greater) was not found. Jet fitting plugin will not be available.") endif() +if(TARGET CGAL::Eigen3_support) + + polyhedron_demo_plugin(interpolated_corrected_principal_curvatures_plugin Interpolated_corrected_principal_curvatures_plugin) + target_link_libraries( + interpolated_corrected_principal_curvatures_plugin PUBLIC scene_surface_mesh_item scene_polylines_item + CGAL::Eigen3_support) + +else() + message( + STATUS + "NOTICE: Eigen 3.1 (or greater) was not found. Interpolated corrected principal curvatures plugin will not be available." + ) +endif() + polyhedron_demo_plugin(extrude_plugin Extrude_plugin KEYWORDS PMP) target_link_libraries(extrude_plugin PUBLIC scene_surface_mesh_item scene_selection_item) diff --git a/Polyhedron/demo/Polyhedron/Plugins/PMP/Interpolated_corrected_principal_curvatures_plugin.cpp b/Polyhedron/demo/Polyhedron/Plugins/PMP/Interpolated_corrected_principal_curvatures_plugin.cpp new file mode 100644 index 00000000000..8bb6a229967 --- /dev/null +++ b/Polyhedron/demo/Polyhedron/Plugins/PMP/Interpolated_corrected_principal_curvatures_plugin.cpp @@ -0,0 +1,180 @@ +#include +#include +#include +#include "Scene_surface_mesh_item.h" +#include "Scene_polylines_item.h" + +#include + +#include "Scene.h" +#include +#include + +#include +#include + +using namespace CGAL::Three; +class Polyhedron_demo_interpolated_corrected_principal_curvatures_and_directions_plugin : + public QObject, + public Polyhedron_demo_plugin_interface +{ + Q_OBJECT + Q_INTERFACES(CGAL::Three::Polyhedron_demo_plugin_interface) + Q_PLUGIN_METADATA(IID "com.geometryfactory.PolyhedronDemo.PluginInterface/1.0") + +public: + + QList actions() const { + return _actions; + } + void init(QMainWindow* mw, + Scene_interface* scene_interface, + Messages_interface*) + { + scene = scene_interface; + QAction *actionEstimateCurvature = new QAction(tr("Interpolated Corrected Principal Curvatures"), mw); + connect(actionEstimateCurvature, SIGNAL(triggered()), this, SLOT(on_actionEstimateCurvature_triggered())); + _actions <(scene->item(scene->mainSelectionIndex())); + } + +public Q_SLOTS: + void on_actionEstimateCurvature_triggered(); +private : + Scene_interface *scene; + QList _actions; +}; // end Polyhedron_demo_interpolated_corrected_principal_curvatures_and_directions_plugin + + +void compute(SMesh* sMesh, + Scene_polylines_item* max_curv, + Scene_polylines_item* min_curv, + Scene_polylines_item* max_negative_curv, + Scene_polylines_item* min_negative_curv) +{ + namespace PMP = CGAL::Polygon_mesh_processing; + typedef CGAL::Exact_predicates_inexact_constructions_kernel Epic_kernel; + typedef Epic_kernel::Point_3 Point; + typedef Epic_kernel::Point_3 Point; + typedef Epic_kernel::Vector_3 Vector; + typedef boost::graph_traits::vertex_descriptor Vertex_descriptor; + + typename boost::property_map::type vpmap = get(CGAL::vertex_point, *sMesh); + + bool created = false; + SMesh::Property_map> principal_curvatures_and_directions_map; + + boost::tie(principal_curvatures_and_directions_map, created) = sMesh->add_property_map> + ("v:principal_curvatures_and_directions_map", { 0, 0, + Vector(0,0,0), + Vector(0,0,0) }); + assert(created); + + + PMP::interpolated_corrected_principal_curvatures_and_directions( + *sMesh, + principal_curvatures_and_directions_map, + CGAL::parameters::ball_radius(0) + ); + + double max_curvature_magnitude_on_mesh = 0; + for (Vertex_descriptor v : vertices(*sMesh)) + { + const PMP::Principal_curvatures_and_directions pc = principal_curvatures_and_directions_map[v]; + max_curvature_magnitude_on_mesh = std::max(max_curvature_magnitude_on_mesh, std::max(abs(pc.min_curvature), abs(pc.max_curvature))); + } + + for(Vertex_descriptor v : vertices(*sMesh)) + { + std::vector points; + + // pick central point + const Point& central_point = get(vpmap,v); + points.push_back(central_point); + + // compute min edge len around central vertex + // to scale the ribbons used to display the directions + + const std::size_t n = CGAL::edges(*sMesh).size(); + + double avg_edge_length = 0; + if (n > 0) { + for (auto e : CGAL::edges(*sMesh)) + avg_edge_length += PMP::edge_length(e, *sMesh); + avg_edge_length /= n; + } + + const PMP::Principal_curvatures_and_directions pc = principal_curvatures_and_directions_map[v]; + + Vector umin = (pc.min_curvature / max_curvature_magnitude_on_mesh) * pc.min_direction * avg_edge_length; + Vector umax = (pc.max_curvature / max_curvature_magnitude_on_mesh) * pc.max_direction * avg_edge_length; + + Scene_polylines_item::Polyline max_segment(2), min_segment(2); + + const double du = 0.4; + + min_segment[0] = central_point + du * umin; + min_segment[1] = central_point - du * umin; + max_segment[0] = central_point + du * umax; + max_segment[1] = central_point - du * umax; + + (pc.min_curvature > 0 ? min_curv : min_negative_curv)->polylines.push_back(min_segment); + (pc.max_curvature > 0 ? max_curv : max_negative_curv)->polylines.push_back(max_segment); + } +} + +void Polyhedron_demo_interpolated_corrected_principal_curvatures_and_directions_plugin::on_actionEstimateCurvature_triggered() +{ + // get active polyhedron + const CGAL::Three::Scene_interface::Item_id index = scene->mainSelectionIndex(); + QString name = scene->item(index)->name(); + Scene_surface_mesh_item* sm_item = + qobject_cast(scene->item(index)); + if(! sm_item){ + return; + } + // wait cursor + QApplication::setOverrideCursor(Qt::WaitCursor); + + // types + Scene_polylines_item* max_curv = new Scene_polylines_item; + max_curv->setColor(Qt::red); + max_curv->setWidth(3); + max_curv->setName(tr("%1 (max curvatures)").arg(name)); + + Scene_polylines_item* min_curv = new Scene_polylines_item; + min_curv->setColor(QColor(255,210,0)); + min_curv->setWidth(4); + min_curv->setName(tr("%1 (min curvatures)").arg(name)); + + Scene_polylines_item* max_negative_curv = new Scene_polylines_item; + max_negative_curv->setColor(Qt::cyan); + max_negative_curv->setWidth(4); + max_negative_curv->setName(tr("%1 (max negative curvatures)").arg(name)); + + Scene_polylines_item* min_negative_curv = new Scene_polylines_item; + min_negative_curv->setColor(Qt::blue); + min_negative_curv->setWidth(3); + min_negative_curv->setName(tr("%1 (min negative curvatures)").arg(name)); + + SMesh* pMesh = sm_item->polyhedron(); + compute(pMesh, max_curv, min_curv, max_negative_curv, min_negative_curv); + + scene->addItem(max_curv); + scene->addItem(min_curv); + max_curv->invalidateOpenGLBuffers(); + min_curv->invalidateOpenGLBuffers(); + scene->addItem(max_negative_curv); + scene->addItem(min_negative_curv); + max_negative_curv->invalidateOpenGLBuffers(); + min_negative_curv->invalidateOpenGLBuffers(); + + // default cursor + QApplication::restoreOverrideCursor(); +} + +#include "Interpolated_corrected_principal_curvatures_plugin.moc" diff --git a/Polyhedron/demo/Polyhedron/Plugins/PMP/Random_perturbation_dialog.ui b/Polyhedron/demo/Polyhedron/Plugins/PMP/Random_perturbation_dialog.ui index 1d368c3dd59..3427d5b12e6 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/PMP/Random_perturbation_dialog.ui +++ b/Polyhedron/demo/Polyhedron/Plugins/PMP/Random_perturbation_dialog.ui @@ -135,6 +135,26 @@ + + + keep vertex normals unperturbed + + + Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter + + + seed_spinbox + + + + + + + + + + + Random seed @@ -147,7 +167,7 @@ - + @@ -163,6 +183,8 @@ seed_spinbox seed_label + keep_normal_label + keep_normal_checkbox deterministic_label deterministic_checkbox diff --git a/Polyhedron/demo/Polyhedron/Plugins/PMP/Random_perturbation_plugin.cpp b/Polyhedron/demo/Polyhedron/Plugins/PMP/Random_perturbation_plugin.cpp index 53499b82336..004768739cd 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/PMP/Random_perturbation_plugin.cpp +++ b/Polyhedron/demo/Polyhedron/Plugins/PMP/Random_perturbation_plugin.cpp @@ -10,6 +10,7 @@ #include #include +#include #include #include @@ -102,6 +103,12 @@ public Q_SLOTS: if (poly_item) { SMesh& pmesh = *poly_item->face_graph(); + if(ui.keep_normal_checkbox->isChecked()) + { + SMesh::Property_map vnormals = + pmesh.add_property_map("v:normal_before_perturbation").first; + CGAL::Polygon_mesh_processing::compute_vertex_normals(pmesh,vnormals); + } if(ui.deterministic_checkbox->isChecked()) { unsigned int seed = static_cast(ui.seed_spinbox->value()); diff --git a/Polyhedron/demo/Polyhedron/cgal_test_with_cmake b/Polyhedron/demo/Polyhedron/cgal_test_with_cmake index 6fea7949176..59eb5c1d021 100755 --- a/Polyhedron/demo/Polyhedron/cgal_test_with_cmake +++ b/Polyhedron/demo/Polyhedron/cgal_test_with_cmake @@ -149,6 +149,7 @@ hole_filling_plugin \ hole_filling_sm_plugin \ hole_filling_polyline_plugin \ inside_out_plugin \ +interpolated_corrected_principal_curvatures_plugin\ surface_intersection_plugin \ surface_intersection_sm_plugin \ io_image_plugin \ diff --git a/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h b/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h index ea2325c0921..71428a87a83 100644 --- a/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h +++ b/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h @@ -93,6 +93,13 @@ CGAL_add_named_parameter(number_of_points_per_edge_t, number_of_points_per_edge, CGAL_add_named_parameter(number_of_points_on_edges_t, number_of_points_on_edges, number_of_points_on_edges) CGAL_add_named_parameter(nb_points_per_area_unit_t, nb_points_per_area_unit, number_of_points_per_area_unit) CGAL_add_named_parameter(nb_points_per_distance_unit_t, nb_points_per_distance_unit, number_of_points_per_distance_unit) +CGAL_add_named_parameter(vertex_mean_curvature_map_t, vertex_mean_curvature_map, vertex_mean_curvature_map) +CGAL_add_named_parameter(vertex_Gaussian_curvature_map_t, vertex_Gaussian_curvature_map, vertex_Gaussian_curvature_map) +CGAL_add_named_parameter(vertex_principal_curvatures_and_directions_map_t, vertex_principal_curvatures_and_directions_map, vertex_principal_curvatures_and_directions_map) +CGAL_add_named_parameter(vertex_mean_curvature_t, vertex_mean_curvature, vertex_mean_curvature) +CGAL_add_named_parameter(vertex_Gaussian_curvature_t, vertex_Gaussian_curvature, vertex_Gaussian_curvature) +CGAL_add_named_parameter(vertex_principal_curvatures_and_directions_t, vertex_principal_curvatures_and_directions, vertex_principal_curvatures_and_directions) +CGAL_add_named_parameter(ball_radius_t, ball_radius, ball_radius) CGAL_add_named_parameter(outward_orientation_t, outward_orientation, outward_orientation) CGAL_add_named_parameter(overlap_test_t, overlap_test, do_overlap_test_of_bounded_sides) CGAL_add_named_parameter(preserve_genus_t, preserve_genus, preserve_genus)