diff --git a/Documentation/doc/biblio/cgal_manual.bib b/Documentation/doc/biblio/cgal_manual.bib index 26b897d23897..e16202c69ddf 100644 --- a/Documentation/doc/biblio/cgal_manual.bib +++ b/Documentation/doc/biblio/cgal_manual.bib @@ -3106,6 +3106,16 @@ @article{ecvp-bhhhk-14 bibsource = {dblp computer science bibliography, https://dblp.org/} } +@inproceedings {dunyach2013curvRemesh, + booktitle = {Eurographics 2013 - Short Papers}, + title = {{Adaptive Remeshing for Real-Time Mesh Deformation}}, + author = {Dunyach, Marion and Vanderhaeghe, David and Barthe, Loïc and Botsch, Mario}, + year = {2013}, + publisher = {The Eurographics Association}, + ISSN = {1017-4656}, + DOI = {10.2312/conf/EG2013/short/029-032} +} + @book{botsch2010PMP, title={Polygon mesh processing}, author={M. Botsch and L. Kobbelt and M. Pauly and P. Alliez and B. L{\'e}vy}, diff --git a/Documentation/doc/scripts/generate_how_to_cite.py b/Documentation/doc/scripts/generate_how_to_cite.py index a891edddebc6..4873b7065e78 100644 --- a/Documentation/doc/scripts/generate_how_to_cite.py +++ b/Documentation/doc/scripts/generate_how_to_cite.py @@ -157,7 +157,7 @@ def protect_upper_case(title): return title.replace("dD","{dD}").replace("2D","{2D}").replace("3D","{3D}").replace("CGAL","{CGAL}").replace("Qt","{Qt}").replace("Boost","{Boost}") def protect_accentuated_letters(authors): - res=authors.replace("é",r"{\'e}").replace("è",r"{\`e}").replace("É",r"{\'E}").replace("ä",r"{\"a}").replace("ö",r"{\"o}").replace("ñ",r"{\~n}").replace("ã",r"{\~a}").replace("ë",r"{\"e}").replace("ı",r"{\i}").replace("Ş",r"{\c{S}}").replace("ş",r"{\c{s}}").replace("%","") + res=authors.replace("é",r"{\'e}").replace("è",r"{\`e}").replace("É",r"{\'E}").replace("ä",r"{\"a}").replace("ö",r"{\"o}").replace("ñ",r"{\~n}").replace("ã",r"{\~a}").replace("ë",r"{\"e}").replace("ı",r"{\i}").replace("Ş",r"{\c{S}}").replace("ş",r"{\c{s}}").replace("%","").replace("đ",r"{\-d}") try: res.encode('ascii') except UnicodeEncodeError: diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/Concepts/PMPSizingField.h b/Polygon_mesh_processing/doc/Polygon_mesh_processing/Concepts/PMPSizingField.h new file mode 100644 index 000000000000..8641c30c6caf --- /dev/null +++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/Concepts/PMPSizingField.h @@ -0,0 +1,66 @@ +/// \ingroup PkgPolygonMeshProcessingConcepts +/// \cgalConcept +/// +/// The concept `PMPSizingField` defines the requirements for the sizing field +/// used in `CGAL::Polygon_mesh_processing::isotropic_remeshing()` to define +/// the target length for every individual edge during the remeshing process. +/// +/// \cgalHasModelsBegin +/// \cgalHasModels{CGAL::Polygon_mesh_processing::Uniform_sizing_field} +/// \cgalHasModels{CGAL::Polygon_mesh_processing::Adaptive_sizing_field} +/// \cgalHasModelsEnd +/// + + +class PMPSizingField{ +public: + +/// @name Types +/// @{ + +/// Vertex descriptor type +typedef boost::graph_traits::vertex_descriptor vertex_descriptor; + +/// Halfedge descriptor type +typedef boost::graph_traits::halfedge_descriptor halfedge_descriptor; + +/// 3D point type matching the value type of the vertex property map passed to `CGAL::Polygon_mesh_processing::isotropic_remeshing()` +typedef unspecified_type Point_3; + +/// Polygon mesh type matching the type passed to `CGAL::Polygon_mesh_processing::isotropic_remeshing()` +typedef unspecified_type PolygonMesh; + +/// Number type matching the `FT` type of the geometric traits passed to `CGAL::Polygon_mesh_processing::isotropic_remeshing()` +typedef unspecified_type FT; + +/// @} + +/// @name Functions +/// @{ + +/// a function that returns the sizing value at `v`. +FT at(const vertex_descriptor v) const; + +/// a function controlling edge split and edge collapse, +/// returning the ratio of the current edge length and the local target edge length between +/// the points of `va` and `vb` in case the current edge is too long, and `std::nullopt` otherwise. +std::optional is_too_long(const vertex_descriptor va, + const vertex_descriptor vb) const; + +/// a function controlling edge collapse by returning the ratio of the squared length of `h` and the +/// local target edge length if it is too short, and `std::nullopt` otherwise. +std::optional is_too_short(const halfedge_descriptor h, + const PolygonMesh& pmesh) const; + +/// a function returning the location of the split point of the edge of `h`. +Point_3 split_placement(const halfedge_descriptor h, + const PolygonMesh& pmesh) const; + +/// a function that updates the sizing field value at the vertex `v`. +void update(const vertex_descriptor v, + const PolygonMesh& pmesh); + +/// @} +}; + + diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/Doxyfile.in b/Polygon_mesh_processing/doc/Polygon_mesh_processing/Doxyfile.in index f6ed67912458..6369cf0c9e94 100644 --- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/Doxyfile.in +++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/Doxyfile.in @@ -37,3 +37,4 @@ HTML_EXTRA_FILES = ${CGAL_PACKAGE_DOC_DIR}/fig/selfintersections.jpg \ ${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 035781b35c7e..a8eae9eb2555 100644 --- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt +++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt @@ -79,7 +79,7 @@ \cgalPkgPicture{hole_filling_ico.png} \cgalPkgSummaryBegin -\cgalPkgAuthors{David Coeurjolly, Jaques-Olivier Lachaud, Konstantinos Katrioplas, Sébastien Loriot, Mael Rouxel-Labbé, Hossam Saeed, Jane Tournois, and Ilker %O. Yaz} +\cgalPkgAuthors{David Coeurjolly, Jaques-Olivier Lachaud, Konstantinos Katrioplas, Sébastien Loriot, Ivan Pađen, 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.} @@ -115,8 +115,6 @@ The page \ref bgl_namedparameters "Named Parameters" describes their usage. - `CGAL::Polygon_mesh_processing::split()` \cgalCRPSection{Meshing Functions} -- \link PMP_meshing_grp `CGAL::Polygon_mesh_processing::isotropic_remeshing()` \endlink -- \link PMP_meshing_grp `CGAL::Polygon_mesh_processing::split_long_edges()` \endlink - `CGAL::Polygon_mesh_processing::remesh_planar_patches()` - `CGAL::Polygon_mesh_processing::remesh_almost_planar_patches()` - `CGAL::Polygon_mesh_processing::refine()` @@ -134,6 +132,10 @@ The page \ref bgl_namedparameters "Named Parameters" describes their usage. - `CGAL::Polygon_mesh_processing::smooth_shape()` - `CGAL::Polygon_mesh_processing::random_perturbation()` +\cgalCRPSection{Sizing Fields} +- `CGAL::Polygon_mesh_processing::Uniform_sizing_field` +- `CGAL::Polygon_mesh_processing::Adaptive_sizing_field` + \cgalCRPSection{Orientation Functions} - `CGAL::Polygon_mesh_processing::orient_polygon_soup()` - `CGAL::Polygon_mesh_processing::orient()` @@ -214,6 +216,7 @@ The page \ref bgl_namedparameters "Named Parameters" describes their usage. \cgalCRPSection{Corrected Curvatures} - `CGAL::Polygon_mesh_processing::interpolated_corrected_curvatures()` +- `CGAL::Polygon_mesh_processing::Principal_curvatures_and_directions` \cgalCRPSection{Normal Computation Functions} - `CGAL::Polygon_mesh_processing::compute_face_normal()` 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 930a03523c90..26d5d9c64507 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 David Coeurjolly, Jaques-Olivier Lachaud, Konstantinos Katrioplas, Sébastien Loriot, Mael Rouxel-Labbé, Hossam Saeed, Jane Tournois, and Ilker %O. Yaz +\authors David Coeurjolly, Jaques-Olivier Lachaud, Konstantinos Katrioplas, Sébastien Loriot, Ivan Pađen, Mael Rouxel-Labbé, Hossam Saeed, Jane Tournois, and Ilker %O. Yaz \image html neptun_head.jpg \image latex neptun_head.jpg @@ -116,10 +116,18 @@ to the original surface to keep a good approximation of the input. A triangulated region of a polygon mesh can be remeshed using the function `CGAL::Polygon_mesh_processing::isotropic_remeshing()`, as illustrated -by \cgalFigureRef{iso_remeshing}. The algorithm has only two parameters : -the target edge length for the remeshed surface patch, and -the number of iterations of the abovementioned sequence of operations. The bigger -this number, the smoother and closer to target edge length the mesh will be. +by \cgalFigureRef{iso_remeshing}. The algorithm has two parameters: +the sizing field object for the remeshed surface patch, and +the number of iterations of the abovementioned sequence of operations. + +The sizing field establishes the local target edge length for the remeshed surface. Two sizing fields are +provided: a uniform and a curvature-adaptive sizing field. With `CGAL::Polygon_mesh_processing::Uniform_sizing_field`, +all triangle edges are targeted to have equal lengths. With `CGAL::Polygon_mesh_processing::Adaptive_sizing_field`, triangle edge lengths depend on the local curvature -- +shorter edges appear in regions with a higher curvature and vice versa. The outline of the adaptive sizing +field algorithm is available in \cgalCite{dunyach2013curvRemesh}. The distinction between uniform and adaptive +sizing fields is depicted in figure \cgalFigureRef{uniform_and_adaptive}. + +As the number of iterations increases, the mesh tends to be smoother and closer to the target edge length. An additional option has been added to \e protect (\e i.\e e. not modify) some given polylines. In some cases, those polylines are too long, and reaching the desired target edge length while protecting them is not @@ -134,6 +142,12 @@ Isotropic remeshing. (a) Triangulated input surface mesh. (d) Surface mesh with the selection uniformly remeshed. \cgalFigureEnd +\cgalFigureBegin{uniform_and_adaptive, uniform_and_adaptive.png} +Sizing fields in isotropic remeshing. +(a) Uniform sizing field. +(b) Curvature-based adaptive sizing field. +\cgalFigureEnd + \paragraph Delaunay-Based Surface Remeshing The mesh generation algorithm implemented in the \ref PkgMesh3 package can be used to remesh a given triangulated surface mesh. The algorithm, based on Delaunay refinement of a restricted Delaunay triangulation, @@ -1287,5 +1301,7 @@ supervision of David Coeurjolly, Jaques-Olivier Lachaud, and Sébastien Loriot. DGtal's implementation was also used as a reference during the project. +The curvature-based sizing field version of isotropic remeshing was added by Ivan Pađen during GSoC 2023, under the supervision of Sébastien Loriot and Jane Tournois. + */ } /* namespace CGAL */ diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/fig/uniform_and_adaptive.png b/Polygon_mesh_processing/doc/Polygon_mesh_processing/fig/uniform_and_adaptive.png new file mode 100644 index 000000000000..4aca8a398d26 Binary files /dev/null and b/Polygon_mesh_processing/doc/Polygon_mesh_processing/fig/uniform_and_adaptive.png differ diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt b/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt index b553b41f737b..590305c0c8ad 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt @@ -69,6 +69,8 @@ if(TARGET CGAL::Eigen3_support) target_link_libraries(hole_filling_example_LCC PUBLIC CGAL::Eigen3_support) create_single_source_cgal_program("mesh_smoothing_example.cpp") target_link_libraries(mesh_smoothing_example PUBLIC CGAL::Eigen3_support) + create_single_source_cgal_program("isotropic_remeshing_with_sizing_example.cpp") + target_link_libraries(isotropic_remeshing_with_sizing_example PUBLIC CGAL::Eigen3_support) create_single_source_cgal_program("delaunay_remeshing_example.cpp") target_link_libraries(delaunay_remeshing_example PUBLIC CGAL::Eigen3_support) create_single_source_cgal_program("remesh_almost_planar_patches.cpp") diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/isotropic_remeshing_with_sizing_example.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/isotropic_remeshing_with_sizing_example.cpp new file mode 100644 index 000000000000..b4e2475f39d2 --- /dev/null +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/isotropic_remeshing_with_sizing_example.cpp @@ -0,0 +1,45 @@ +#include +#include +#include +#include +#include + +#include + +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Surface_mesh Mesh; + +namespace PMP = CGAL::Polygon_mesh_processing; + +int main(int argc, char* argv[]) +{ + const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/nefertiti.off"); + + Mesh mesh; + if (!PMP::IO::read_polygon_mesh(filename, mesh) || !CGAL::is_triangle_mesh(mesh)) { + std::cerr << "Not a valid input file." << std::endl; + return 1; + } + + std::cout << "Start remeshing of " << filename + << " (" << num_faces(mesh) << " faces)..." << std::endl; + + const double tol = 0.001; + const std::pair edge_min_max{0.001, 0.5}; + PMP::Adaptive_sizing_field sizing_field(tol, edge_min_max, faces(mesh), mesh); + unsigned int nb_iter = 5; + + PMP::isotropic_remeshing( + faces(mesh), + sizing_field, + mesh, + CGAL::parameters::number_of_iterations(nb_iter) + .number_of_relaxation_steps(3) + ); + + CGAL::IO::write_polygon_mesh("out.off", mesh, CGAL::parameters::stream_precision(17)); + + std::cout << "Remeshing done." << std::endl; + + return 0; +} diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/Adaptive_sizing_field.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/Adaptive_sizing_field.h new file mode 100644 index 000000000000..a61edc699d06 --- /dev/null +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/Adaptive_sizing_field.h @@ -0,0 +1,280 @@ +// Copyright (c) 2023 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) : Ivan Paden + +#ifndef CGAL_PMP_REMESHING_ADAPTIVE_SIZING_FIELD_H +#define CGAL_PMP_REMESHING_ADAPTIVE_SIZING_FIELD_H + +#include + +#include +#include + +#include +#include + +#include + +namespace CGAL +{ +namespace Polygon_mesh_processing +{ +/*! +* \ingroup PMP_meshing_grp +* a sizing field describing variable target mesh edge lengths for +* `CGAL::Polygon_mesh_processing::isotropic_remeshing()`. +* This adaptive sizing field is a function of local discrete curvatures, +* computed using the +* `CGAL::Polygon_mesh_processing::interpolated_corrected_curvatures()` function. +* +* Edges too long with respect to the local target edge length are split in two, while +* edges that are too short are collapsed. +* +* \cgalModels{PMPSizingField} +* +* \sa `isotropic_remeshing()` +* \sa `Uniform_sizing_field` +* +* @tparam PolygonMesh model of `MutableFaceGraph` that +* has an internal property map for `CGAL::vertex_point_t`. +*/ +template ::const_type> +class Adaptive_sizing_field +#ifndef DOXYGEN_RUNNING +: public internal::Sizing_field_base +#endif +{ +private: + typedef internal::Sizing_field_base Base; + typedef typename CGAL::dynamic_vertex_property_t Vertex_property_tag; + typedef typename boost::property_map::type VertexSizingMap; + +public: + typedef typename Base::K K; + typedef typename Base::FT FT; + typedef typename Base::Point_3 Point_3; + typedef typename Base::face_descriptor face_descriptor; + typedef typename Base::halfedge_descriptor halfedge_descriptor; + typedef typename Base::vertex_descriptor vertex_descriptor; + + /// \name Creation + /// @{ + /*! + * Constructor + * + * @tparam FaceRange range of `boost::graph_traits::%face_descriptor`, + * model of `Range`. Its iterator type is `ForwardIterator`. + * @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters" + * + * @param tol the error tolerance, used together with curvature to derive target edge length. + * Lower tolerance values will result in shorter mesh edges. + * @param edge_len_min_max contains the bounds for minimum and maximum + * edge lengths + * @param face_range the range of triangular faces defining one or several surface patches + * to be remeshed. It should be the same as the range of faces passed to `isotropic_remeshing()`. + * @param pmesh a polygon mesh with triangulated surface patches to be remeshed. It should be the + * same mesh as the one passed to `isotropic_remeshing()`. + * @param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below + + * \cgalNamedParamsBegin + * \cgalParamNBegin{vertex_point_map} + * \cgalParamDescription{a property map associating points to the vertices of `pmesh`} + * \cgalParamType{a class model of `ReadWritePropertyMap` 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{ball_radius} + * \cgalParamDescription{`ball_radius` parameter passed to `interpolated_corrected_curvatures()`} + * \cgalParamNEnd + * \cgalNamedParamsEnd + */ + template + Adaptive_sizing_field(const FT tol + , const std::pair& edge_len_min_max + , const FaceRange& face_range + , PolygonMesh& pmesh + , const NamedParameters& np = parameters::default_values()) + : tol(tol) + , m_short(edge_len_min_max.first) + , m_long(edge_len_min_max.second) + , m_vpmap(parameters::choose_parameter( + parameters::get_parameter(np, internal_np::vertex_point), + get_property_map(vertex_point, pmesh))) + , m_vertex_sizing_map(get(Vertex_property_tag(), pmesh)) + { + if (face_range.size() == faces(pmesh).size()) + { + // calculate curvature from the whole mesh + calc_sizing_map(pmesh, np); + } + else + { + // expand face selection and calculate curvature from it + std::vector selection(face_range.begin(), face_range.end()); + auto is_selected = get(CGAL::dynamic_face_property_t(), pmesh); + for (face_descriptor f : faces(pmesh)) put(is_selected, f, false); + for (face_descriptor f : face_range) put(is_selected, f, true); + expand_face_selection(selection, pmesh, 1, + is_selected, std::back_inserter(selection)); + Face_filtered_graph ffg(pmesh, selection); + + calc_sizing_map(ffg, np); + } + } + + ///@} + +private: + template + void calc_sizing_map(FaceGraph& face_graph + , const NamedParameters& np) + { + typedef Principal_curvatures_and_directions Principal_curvatures; + typedef typename CGAL::dynamic_vertex_property_t Vertex_curvature_tag; + typedef typename boost::property_map::type Vertex_curvature_map; + + using parameters::choose_parameter; + using parameters::get_parameter; + typename Base::FT radius = choose_parameter(get_parameter(np, internal_np::ball_radius), -1); + +#ifdef CGAL_PMP_REMESHING_VERBOSE + int oversize = 0; + int undersize = 0; + int insize = 0; + std::cout << "Calculating sizing field..." << std::endl; +#endif + + Vertex_curvature_map vertex_curvature_map = get(Vertex_curvature_tag(), face_graph); + interpolated_corrected_curvatures(face_graph, + parameters::vertex_principal_curvatures_and_directions_map(vertex_curvature_map) + .ball_radius(radius)); + // calculate vertex sizing field L(x_i) from the curvature field + for(vertex_descriptor v : vertices(face_graph)) + { + auto vertex_curv = get(vertex_curvature_map, v); + const FT max_absolute_curv = (CGAL::max)(CGAL::abs(vertex_curv.max_curvature), + CGAL::abs(vertex_curv.min_curvature)); + const FT vertex_size_sq = 6 * tol / max_absolute_curv - 3 * CGAL::square(tol); + if (vertex_size_sq > CGAL::square(m_long)) + { + put(m_vertex_sizing_map, v, m_long); +#ifdef CGAL_PMP_REMESHING_VERBOSE + ++oversize; +#endif + } + else if (vertex_size_sq < CGAL::square(m_short)) + { + put(m_vertex_sizing_map, v, m_short); +#ifdef CGAL_PMP_REMESHING_VERBOSE + ++undersize; +#endif + } + else + { + put(m_vertex_sizing_map, v, CGAL::approximate_sqrt(vertex_size_sq)); +#ifdef CGAL_PMP_REMESHING_VERBOSE + ++insize; +#endif + } + } +#ifdef CGAL_PMP_REMESHING_VERBOSE + std::cout << " done (" << insize << " from curvature, " + << oversize << " set to max, " + << undersize << " set to min)" << std::endl; +#endif + } + + FT sqlength(const vertex_descriptor va, + const vertex_descriptor vb) const + { + return FT(CGAL::squared_distance(get(m_vpmap, va), get(m_vpmap, vb))); + } + + FT sqlength(const halfedge_descriptor& h, const PolygonMesh& pmesh) const + { + return sqlength(target(h, pmesh), source(h, pmesh)); + } + +public: + FT at(const vertex_descriptor v) const + { + CGAL_assertion(get(m_vertex_sizing_map, v)); + return get(m_vertex_sizing_map, v); + } + + std::optional is_too_long(const vertex_descriptor va, const vertex_descriptor vb) const + { + const FT sqlen = sqlength(va, vb); + FT sqtarg_len = CGAL::square(4./3. * (CGAL::min)(get(m_vertex_sizing_map, va), + get(m_vertex_sizing_map, vb))); + CGAL_assertion(get(m_vertex_sizing_map, va)); + CGAL_assertion(get(m_vertex_sizing_map, vb)); + if (sqlen > sqtarg_len) + return sqlen / sqtarg_len; + else + return std::nullopt; + } + + std::optional is_too_short(const halfedge_descriptor h, const PolygonMesh& pmesh) const + { + const FT sqlen = sqlength(h, pmesh); + FT sqtarg_len = CGAL::square(4./5. * (CGAL::min)(get(m_vertex_sizing_map, source(h, pmesh)), + get(m_vertex_sizing_map, target(h, pmesh)))); + CGAL_assertion(get(m_vertex_sizing_map, source(h, pmesh))); + CGAL_assertion(get(m_vertex_sizing_map, target(h, pmesh))); + if (sqlen < sqtarg_len) + return sqlen / sqtarg_len; + else + return std::nullopt; + } + + Point_3 split_placement(const halfedge_descriptor h, const PolygonMesh& pmesh) const + { + return midpoint(get(m_vpmap, target(h, pmesh)), + get(m_vpmap, source(h, pmesh))); + } + + void update(const vertex_descriptor v, const PolygonMesh& pmesh) + { + // calculating it as the average of two vertices on other ends + // of halfedges as updating is done during an edge split + FT vertex_size = 0; + CGAL_assertion(CGAL::halfedges_around_target(v, pmesh).size() == 2); + for (halfedge_descriptor ha: CGAL::halfedges_around_target(v, pmesh)) + { + vertex_size += get(m_vertex_sizing_map, source(ha, pmesh)); + } + vertex_size /= CGAL::halfedges_around_target(v, pmesh).size(); + + put(m_vertex_sizing_map, v, vertex_size); + } + +private: + const FT tol; + const FT m_short; + const FT m_long; + const VPMap m_vpmap; + VertexSizingMap m_vertex_sizing_map; +}; + +}//end namespace Polygon_mesh_processing +}//end namespace CGAL + +#endif //CGAL_PMP_REMESHING_ADAPTIVE_SIZING_FIELD_H diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/Uniform_sizing_field.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/Uniform_sizing_field.h new file mode 100644 index 000000000000..a4144ceee650 --- /dev/null +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/Uniform_sizing_field.h @@ -0,0 +1,149 @@ +// Copyright (c) 2020 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) : Jane Tournois + +#ifndef CGAL_PMP_REMESHING_UNIFORM_SIZING_FIELD_H +#define CGAL_PMP_REMESHING_UNIFORM_SIZING_FIELD_H + +#include + +#include + +#include + +namespace CGAL +{ +namespace Polygon_mesh_processing +{ +/*! +* \ingroup PMP_meshing_grp +* a sizing field describing a uniform target edge length for +* `CGAL::Polygon_mesh_processing::isotropic_remeshing()`. +* +* Edges longer than 4/3 of the target edge length will be split in half, while +* edges shorter than 4/5 of the target edge length will be collapsed. +* +* \cgalModels{PMPSizingField} +* +* \sa `isotropic_remeshing()` +* \sa `Adaptive_sizing_field` +* +* @tparam PolygonMesh model of `MutableFaceGraph` that +* has an internal property map for `CGAL::vertex_point_t`. +* @tparam VPMap property map associating points to the vertices of `pmesh`, +* model of `ReadWritePropertyMap` with `boost::graph_traits::%vertex_descriptor` +* as key type and `%Point_3` as value type. Default is `boost::get(CGAL::vertex_point, pmesh)`. +*/ +template ::const_type> +class Uniform_sizing_field +#ifndef DOXYGEN_RUNNING +: public internal::Sizing_field_base +#endif +{ +private: + typedef internal::Sizing_field_base Base; + +public: + typedef typename Base::FT FT; + typedef typename Base::Point_3 Point_3; + typedef typename Base::halfedge_descriptor halfedge_descriptor; + typedef typename Base::vertex_descriptor vertex_descriptor; + + /// \name Creation + /// @{ + + /*! + * Constructor. + * \param size the target edge length for isotropic remeshing. If set to 0, + * the criterion for edge length is ignored and edges are neither split nor collapsed. + * \param vpmap is the input vertex point map that associates points to the vertices of + * the input mesh. + */ + Uniform_sizing_field(const FT size, const VPMap& vpmap) + : m_size(size) + , m_sq_short( CGAL::square(4./5. * size)) + , m_sq_long( CGAL::square(4./3. * size)) + , m_vpmap(vpmap) + {} + + /*! + * Constructor using internal vertex point map of the input polygon mesh. + * + * @param size the target edge length for isotropic remeshing. If set to 0, + * the criterion for edge length is ignored and edges are neither split nor collapsed. + * @param pmesh a polygon mesh with triangulated surface patches to be remeshed. The default + * vertex point map of `pmesh` is used to construct the class. + */ + Uniform_sizing_field(const FT size, const PolygonMesh& pmesh) + : Uniform_sizing_field(size, get(CGAL::vertex_point, pmesh)) + {} + + /// @} + +private: + FT sqlength(const vertex_descriptor va, + const vertex_descriptor vb) const + { + return FT(squared_distance(get(m_vpmap, va), get(m_vpmap, vb))); + } + + FT sqlength(const halfedge_descriptor& h, const PolygonMesh& pmesh) const + { + return sqlength(target(h, pmesh), source(h, pmesh)); + } + +public: + FT at(const vertex_descriptor /* v */) const + { + return m_size; + } + + std::optional is_too_long(const vertex_descriptor va, const vertex_descriptor vb) const + { + const FT sqlen = sqlength(va, vb); + if (sqlen > m_sq_long) + //no need to return the ratio for the uniform field + return sqlen; + else + return std::nullopt; + } + + std::optional is_too_short(const halfedge_descriptor h, const PolygonMesh& pmesh) const + { + const FT sqlen = sqlength(h, pmesh); + if (sqlen < m_sq_short) + //no need to return the ratio for the uniform field + return sqlen; + else + return std::nullopt; + } + + Point_3 split_placement(const halfedge_descriptor h, const PolygonMesh& pmesh) const + { + return midpoint(get(m_vpmap, target(h, pmesh)), + get(m_vpmap, source(h, pmesh))); + } + + void update(const vertex_descriptor /* v */, const PolygonMesh& /* pmesh */) + {} + +private: + const FT m_size; + const FT m_sq_short; + const FT m_sq_long; + const VPMap m_vpmap; +}; + +}//end namespace Polygon_mesh_processing +}//end namespace CGAL + +#endif //CGAL_PMP_REMESHING_UNIFORM_SIZING_FIELD_H diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h index 458cdf008849..c1696428d756 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h @@ -43,7 +43,6 @@ #include #include #include -#include #include #include @@ -54,6 +53,7 @@ #include #include #include +#include #ifdef CGAL_PMP_REMESHING_DEBUG #include @@ -72,7 +72,11 @@ #endif namespace CGAL { + namespace Polygon_mesh_processing { + +template class Uniform_sizing_field; + namespace internal { enum Halfedge_status { @@ -226,15 +230,13 @@ namespace internal { template + typename FacePatchMap, + typename SizingFunction> bool constraints_are_short_enough(const PM& pmesh, EdgeConstraintMap ecmap, - VertexPointMap vpmap, const FacePatchMap& fpm, - const double& high) + const SizingFunction& sizing) { - double sqh = high*high; typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; typedef typename boost::graph_traits::edge_descriptor edge_descriptor; for(edge_descriptor e : edges(pmesh)) @@ -244,8 +246,7 @@ namespace internal { get(ecmap, e) || get(fpm, face(h,pmesh))!=get(fpm, face(opposite(h,pmesh),pmesh)) ) { - if (sqh < CGAL::squared_distance(get(vpmap, source(h, pmesh)), - get(vpmap, target(h, pmesh)))) + if (sizing.is_too_long(source(h, pmesh), target(h, pmesh))) { return false; } @@ -377,19 +378,17 @@ namespace internal { } } - // split edges of edge_range that have their length > high // Note: only used to split a range of edges provided as input - template + template void split_long_edges(const EdgeRange& edge_range, - const double& high) + SizingFunction& sizing) { #ifdef CGAL_PMP_REMESHING_VERBOSE - std::cout << "Split long edges (" << high << ")..."; + std::cout << "Split long edges..."; std::cout.flush(); #endif - double sq_high = high*high; //collect long edges typedef std::pair H_and_sql; @@ -400,9 +399,10 @@ namespace internal { ); for(edge_descriptor e : edge_range) { - double sqlen = sqlength(e); - if (sqlen > sq_high) - long_edges.emplace(halfedge(e, mesh_), sqlen); + const halfedge_descriptor he = halfedge(e, mesh_); + std::optional sqlen = sizing.is_too_long(source(he, mesh_), target(he, mesh_)); + if(sqlen != std::nullopt) + long_edges.emplace(he, sqlen.value()); } //split long edges @@ -414,7 +414,6 @@ namespace internal { //the edge with longest length auto eit = long_edges.begin(); halfedge_descriptor he = eit->first; - double sqlen = eit->second; long_edges.erase(eit); //split edge @@ -433,15 +432,19 @@ namespace internal { #ifdef CGAL_PMP_REMESHING_VERY_VERBOSE std::cout << " refinement point : " << refinement_point << std::endl; #endif + //update sizing field with the new point + sizing.update(vnew, mesh_); //check sub-edges - double sqlen_new = 0.25 * sqlen; - if (sqlen_new > sq_high) - { - //if it was more than twice the "long" threshold, insert them - long_edges.emplace(hnew, sqlen_new); - long_edges.emplace(next(hnew, mesh_), sqlen_new); - } + //if it was more than twice the "long" threshold, insert them + std::optional sqlen_new = sizing.is_too_long(source(hnew, mesh_), target(hnew, mesh_)); + if(sqlen_new != std::nullopt) + long_edges.emplace(hnew, sqlen_new.value()); + + const halfedge_descriptor hnext = next(hnew, mesh_); + sqlen_new = sizing.is_too_long(source(hnext, mesh_), target(hnext, mesh_)); + if (sqlen_new != std::nullopt) + long_edges.emplace(hnext, sqlen_new.value()); //insert new edges to keep triangular faces, and update long_edges if (!is_border(hnew, mesh_)) @@ -478,13 +481,12 @@ namespace internal { // "visits all edges of the mesh //if an edge is longer than the given threshold `high`, the edge //is split at its midpoint and the two adjacent triangles are bisected (2-4 split)" - void split_long_edges(const double& high) + template + void split_long_edges(SizingFunction& sizing) { #ifdef CGAL_PMP_REMESHING_VERBOSE - std::cout << "Split long edges (" << high << ")..." << std::endl; + std::cout << "Split long edges..." << std::endl; #endif - double sq_high = high*high; - //collect long edges typedef std::pair H_and_sql; std::multiset< H_and_sql, std::function > @@ -497,9 +499,10 @@ namespace internal { { if (!is_split_allowed(e)) continue; - double sqlen = sqlength(e); - if(sqlen > sq_high) - long_edges.emplace(halfedge(e, mesh_), sqlen); + const halfedge_descriptor he = halfedge(e, mesh_); + std::optional sqlen = sizing.is_too_long(source(he, mesh_), target(he, mesh_)); + if(sqlen != std::nullopt) + long_edges.emplace(halfedge(e, mesh_), sqlen.value()); } //split long edges @@ -511,7 +514,6 @@ namespace internal { //the edge with longest length auto eit = long_edges.begin(); halfedge_descriptor he = eit->first; - double sqlen = eit->second; long_edges.erase(eit); #ifdef CGAL_PMP_REMESHING_VERBOSE_PROGRESS @@ -528,7 +530,7 @@ namespace internal { Patch_id patch_id_opp = get_patch_id(face(opposite(he, mesh_), mesh_)); //split edge - Point refinement_point = this->midpoint(he); + Point refinement_point = sizing.split_placement(he, mesh_); halfedge_descriptor hnew = CGAL::Euler::split_edge(he, mesh_); CGAL_assertion(he == next(hnew, mesh_)); put(ecmap_, edge(hnew, mesh_), get(ecmap_, edge(he, mesh_)) ); @@ -547,14 +549,19 @@ namespace internal { halfedge_added(hnew, status(he)); halfedge_added(hnew_opp, status(opposite(he, mesh_))); + //update sizing field with the new point + sizing.update(vnew, mesh_); + //check sub-edges - double sqlen_new = 0.25 * sqlen; - if (sqlen_new > sq_high) - { - //if it was more than twice the "long" threshold, insert them - long_edges.emplace(hnew, sqlen_new); - long_edges.emplace(next(hnew, mesh_), sqlen_new); - } + //if it was more than twice the "long" threshold, insert them + std::optional sqlen_new = sizing.is_too_long(source(hnew, mesh_), target(hnew, mesh_)); + if(sqlen_new != std::nullopt) + long_edges.emplace(hnew, sqlen_new.value()); + + const halfedge_descriptor hnext = next(hnew, mesh_); + sqlen_new = sizing.is_too_long(source(hnext, mesh_), target(hnext, mesh_)); + if (sqlen_new != std::nullopt) + long_edges.emplace(hnext, sqlen_new.value()); //insert new edges to keep triangular faces, and update long_edges if (!is_on_border(hnew)) @@ -573,9 +580,9 @@ namespace internal { if (snew == PATCH) { - double sql = sqlength(hnew2); - if (sql > sq_high) - long_edges.emplace(hnew2, sql); + std::optional sql = sizing.is_too_long(source(hnew2, mesh_), target(hnew2, mesh_)); + if(sql != std::nullopt) + long_edges.emplace(hnew2, sql.value()); } } @@ -596,9 +603,9 @@ namespace internal { if (snew == PATCH) { - double sql = sqlength(hnew2); - if (sql > sq_high) - long_edges.emplace(hnew2, sql); + std::optional sql = sizing.is_too_long(source(hnew2, mesh_), target(hnew2, mesh_)); + if (sql != std::nullopt) + long_edges.emplace(hnew2, sql.value()); } } } @@ -620,8 +627,8 @@ namespace internal { // "collapses and thus removes all edges that are shorter than a // threshold `low`. [...] testing before each collapse whether the collapse // would produce an edge that is longer than `high`" - void collapse_short_edges(const double& low, - const double& high, + template + void collapse_short_edges(const SizingFunction& sizing, const bool collapse_constraints) { typedef boost::bimap< @@ -630,22 +637,21 @@ namespace internal { typedef typename Boost_bimap::value_type short_edge; #ifdef CGAL_PMP_REMESHING_VERBOSE - std::cout << "Collapse short edges (" << low << ", " << high << ")..." + std::cout << "Collapse short edges..." << std::endl; #endif #ifdef CGAL_PMP_REMESHING_VERBOSE_PROGRESS std::cout << "Fill bimap..."; std::cout.flush(); #endif - double sq_low = low*low; - double sq_high = high*high; Boost_bimap short_edges; for(edge_descriptor e : edges(mesh_)) { - double sqlen = sqlength(e); - if ((sqlen < sq_low) && is_collapse_allowed(e, collapse_constraints)) - short_edges.insert(short_edge(halfedge(e, mesh_), sqlen)); + std::optional sqlen = sizing.is_too_short(halfedge(e, mesh_), mesh_); + if(sqlen != std::nullopt + && is_collapse_allowed(e, collapse_constraints)) + short_edges.insert(short_edge(halfedge(e, mesh_), sqlen.value())); } #ifdef CGAL_PMP_REMESHING_VERBOSE_PROGRESS std::cout << "done." << std::endl; @@ -741,7 +747,8 @@ namespace internal { for(halfedge_descriptor ha : halfedges_around_target(va, mesh_)) { vertex_descriptor va_i = source(ha, mesh_); - if (sqlength(vb, va_i) > sq_high) + std::optional sqha = sizing.is_too_long(vb, va_i); + if (sqha != std::nullopt) { collapse_ok = false; break; @@ -796,7 +803,7 @@ namespace internal { //fix constrained case CGAL_assertion((is_constrained(vkept) || is_corner(vkept) || is_on_patch_border(vkept)) == (is_va_constrained || is_vb_constrained || is_va_on_constrained_polyline || is_vb_on_constrained_polyline)); - if (fix_degenerate_faces(vkept, short_edges, sq_low, collapse_constraints)) + if (fix_degenerate_faces(vkept, short_edges, sizing, collapse_constraints)) { #ifdef CGAL_PMP_REMESHING_DEBUG debug_status_map(); @@ -806,9 +813,10 @@ namespace internal { //insert new/remaining short edges for (halfedge_descriptor ht : halfedges_around_target(vkept, mesh_)) { - double sqlen = sqlength(ht); - if ((sqlen < sq_low) && is_collapse_allowed(edge(ht, mesh_), collapse_constraints)) - short_edges.insert(short_edge(ht, sqlen)); + std::optional sqlen = sizing.is_too_short(ht, mesh_); + if (sqlen != std::nullopt + && is_collapse_allowed(edge(ht, mesh_), collapse_constraints)) + short_edges.insert(short_edge(ht, sqlen.value())); } } }//end if(collapse_ok) @@ -1006,8 +1014,10 @@ namespace internal { // "applies an iterative smoothing filter to the mesh. // The vertex movement has to be constrained to the vertex tangent plane [...] // smoothing algorithm with uniform Laplacian weights" + template void tangential_relaxation_impl(const bool relax_constraints/*1d smoothing*/ - , const unsigned int nb_iterations) + , const unsigned int nb_iterations + , const SizingFunction& sizing) { #ifdef CGAL_PMP_REMESHING_VERBOSE std::cout << "Tangential relaxation (" << nb_iterations << " iter.)..."; @@ -1038,16 +1048,41 @@ namespace internal { auto constrained_vertices_pmap = boost::make_function_property_map(vertex_constraint); - tangential_relaxation( - vertices(mesh_), - mesh_, - CGAL::parameters::number_of_iterations(nb_iterations) - .vertex_point_map(vpmap_) - .geom_traits(gt_) - .edge_is_constrained_map(constrained_edges_pmap) - .vertex_is_constrained_map(constrained_vertices_pmap) - .relax_constraints(relax_constraints) - ); + if constexpr (std::is_same_v>) + { +#ifdef CGAL_PMP_REMESHING_VERBOSE + std::cout << " using tangential relaxation with weights equal to 1"; + std::cout << std::endl; +#endif + tangential_relaxation( + vertices(mesh_), + mesh_, + CGAL::parameters::number_of_iterations(nb_iterations) + .vertex_point_map(vpmap_) + .geom_traits(gt_) + .edge_is_constrained_map(constrained_edges_pmap) + .vertex_is_constrained_map(constrained_vertices_pmap) + .relax_constraints(relax_constraints) + ); + } + else + { +#ifdef CGAL_PMP_REMESHING_VERBOSE + std::cout << " using tangential relaxation weighted with the sizing field"; + std::cout << std::endl; +#endif + tangential_relaxation( + vertices(mesh_), + mesh_, + CGAL::parameters::number_of_iterations(nb_iterations) + .vertex_point_map(vpmap_) + .geom_traits(gt_) + .edge_is_constrained_map(constrained_edges_pmap) + .vertex_is_constrained_map(constrained_vertices_pmap) + .relax_constraints(relax_constraints) + .sizing_function(sizing) + ); + } CGAL_assertion(!input_mesh_is_valid_ || is_valid_polygon_mesh(mesh_)); @@ -1645,10 +1680,10 @@ namespace internal { // else keep current status for en and eno } - template + template bool fix_degenerate_faces(const vertex_descriptor& v, Bimap& short_edges, - const double& sq_low, + const SizingFunction& sizing, const bool collapse_constraints) { std::unordered_set degenerate_faces; @@ -1726,9 +1761,9 @@ namespace internal { //insert new edges in 'short_edges' if (is_collapse_allowed(edge(hf, mesh_), collapse_constraints)) { - double sqlen = sqlength(hf); - if (sqlen < sq_low) - short_edges.insert(typename Bimap::value_type(hf, sqlen)); + std::optional sqlen = sizing.is_too_short(hf, mesh_); + if (sqlen != std::nullopt) + short_edges.insert(typename Bimap::value_type(hf, sqlen.value())); } if(!is_border(hf, mesh_) && diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Sizing_field_base.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Sizing_field_base.h new file mode 100644 index 000000000000..f13135cd873e --- /dev/null +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Sizing_field_base.h @@ -0,0 +1,77 @@ +// Copyright (c) 2020 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) : Jane Tournois + +#ifndef CGAL_PMP_REMESHING_SIZING_FIELD_H +#define CGAL_PMP_REMESHING_SIZING_FIELD_H + +#include + +#include +#include + +#include + +namespace CGAL +{ +namespace Polygon_mesh_processing +{ +namespace internal +{ +/*! +* \ingroup PMP_meshing_grp +* pure virtual class serving as a base for sizing field classes used in isotropic +* remeshing. +* +* \cgalModels{PMPSizingField} +* +* \sa `isotropic_remeshing()` +* \sa `Uniform_sizing_field` +* \sa `Adaptive_sizing_field` +* +* @tparam PolygonMesh model of `MutableFaceGraph` that +* has an internal property map for `CGAL::vertex_point_t`. +* @tparam VPMap property map associating points to the vertices of `pmesh`, +* model of `ReadWritePropertyMap` with `boost::graph_traits::%vertex_descriptor` +* as key type and `%Point_3` as value type. Default is `boost::get(CGAL::vertex_point, pmesh)`. +*/ +template ::const_type> +class Sizing_field_base +{ +private: + typedef PolygonMesh PM; + typedef typename boost::property_traits::value_type Point; + +public: + typedef typename CGAL::Kernel_traits::Kernel K; + typedef typename boost::graph_traits::face_descriptor face_descriptor; + typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef Point Point_3; + typedef typename K::FT FT; + +public: + virtual FT at(const vertex_descriptor v) const = 0; + virtual std::optional is_too_long(const vertex_descriptor va, + const vertex_descriptor vb) const = 0; + virtual std::optional is_too_short(const halfedge_descriptor h, + const PolygonMesh& pmesh) const = 0; + virtual Point_3 split_placement(const halfedge_descriptor h, const PolygonMesh& pmesh) const = 0; + virtual void update(const vertex_descriptor v, const PolygonMesh& pmesh) = 0; + +}; + +}//end namespace internal +}//end namespace Polygon_mesh_processing +}//end namespace CGAL + +#endif //CGAL_PMP_REMESHING_SIZING_FIELD_H diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/interpolated_corrected_curvatures.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/interpolated_corrected_curvatures.h index 25de9ec9387f..a8d45d8b867f 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/interpolated_corrected_curvatures.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/interpolated_corrected_curvatures.h @@ -617,12 +617,11 @@ 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() - ) + typename NamedParameters> +void interpolated_corrected_curvatures_one_vertex( + const PolygonMesh& pmesh, + const typename boost::graph_traits::vertex_descriptor v, + const NamedParameters& np) { typedef typename GetGeomTraits::type GT; typedef typename GetVertexPointMap::const_type Vertex_position_map; @@ -716,7 +715,7 @@ template +template class Interpolated_corrected_curvatures_computer { typedef typename GetGeomTraits::type GT; @@ -780,7 +779,6 @@ class Interpolated_corrected_curvatures_computer 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) @@ -824,10 +822,8 @@ class Interpolated_corrected_curvatures_computer public: - Interpolated_corrected_curvatures_computer(const PolygonMesh& pmesh, - const NamedParameters& np = parameters::default_values() - ) : - pmesh(pmesh) + Interpolated_corrected_curvatures_computer(const PolygonMesh& pmesh,const NamedParameters& np) + : pmesh(pmesh) { set_named_params(np); diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/remesh.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/remesh.h index 014ae0af1410..7c5ff9ab2618 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/remesh.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/remesh.h @@ -18,6 +18,7 @@ #include #include +#include #include #include @@ -43,12 +44,15 @@ namespace Polygon_mesh_processing { * and `boost::graph_traits::%halfedge_descriptor` must be * models of `Hashable`. * @tparam FaceRange range of `boost::graph_traits::%face_descriptor`, - model of `Range`. Its iterator type is `ForwardIterator`. +* model of `Range`. Its iterator type is `ForwardIterator`. +* @tparam SizingFunction model of `PMPSizingField` * @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters" * * @param pmesh a polygon mesh with triangulated surface patches to be remeshed * @param faces the range of triangular faces defining one or several surface patches to be remeshed -* @param target_edge_length the edge length that is targeted in the remeshed patch. +* @param sizing field that determines a target length for individual edges. +* If a number convertible to a `double` is passed, `Uniform_sizing_field()` will be used, +* with the number as a target edge length. * If `0` is passed then only the edge-flip, tangential relaxation, and projection steps will be done. * @param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below * @@ -193,9 +197,11 @@ namespace Polygon_mesh_processing { */ template + , typename SizingFunction + , typename NamedParameters = parameters::Default_named_parameters + , typename = typename std::enable_if_t>> void isotropic_remeshing(const FaceRange& faces - , const double& target_edge_length + , SizingFunction& sizing , PolygonMesh& pmesh , const NamedParameters& np = parameters::default_values()) { @@ -261,9 +267,6 @@ void isotropic_remeshing(const FaceRange& faces #endif ) ) ); - double low = 4. / 5. * target_edge_length; - double high = 4. / 3. * target_edge_length; - #if !defined(CGAL_NO_PRECONDITIONS) if(protect) { @@ -271,7 +274,7 @@ void isotropic_remeshing(const FaceRange& faces msg.append(" true with constraints larger than 4/3 * target_edge_length."); msg.append(" Remeshing aborted."); CGAL_precondition_msg( - internal::constraints_are_short_enough(pmesh, ecmap, vpmap, fpmap, high), + internal::constraints_are_short_enough(pmesh, ecmap, fpmap, sizing), msg.c_str()); } #endif @@ -303,8 +306,7 @@ void isotropic_remeshing(const FaceRange& faces #ifdef CGAL_PMP_REMESHING_VERBOSE std::cout << std::endl; - std::cout << "Remeshing (size = " << target_edge_length; - std::cout << ", #iter = " << nb_iterations << ")..." << std::endl; + std::cout << "Remeshing (#iter = " << nb_iterations << ")..." << std::endl; t.reset(); t.start(); #endif @@ -313,16 +315,14 @@ void isotropic_remeshing(const FaceRange& faces #ifdef CGAL_PMP_REMESHING_VERBOSE std::cout << " * Iteration " << (i + 1) << " *" << std::endl; #endif - if (target_edge_length>0) - { - if(do_split) - remesher.split_long_edges(high); - if(do_collapse) - remesher.collapse_short_edges(low, high, collapse_constraints); - } + + if(do_split) + remesher.split_long_edges(sizing); + if(do_collapse) + remesher.collapse_short_edges(sizing, collapse_constraints); if(do_flip) remesher.flip_edges_for_valence_and_shape(); - remesher.tangential_relaxation_impl(smoothing_1d, nb_laplacian); + remesher.tangential_relaxation_impl(smoothing_1d, nb_laplacian, sizing); if ( choose_parameter(get_parameter(np, internal_np::do_project), true) ) remesher.project_to_surface(get_parameter(np, internal_np::projection_functor)); #ifdef CGAL_PMP_REMESHING_VERBOSE @@ -332,12 +332,36 @@ void isotropic_remeshing(const FaceRange& faces #ifdef CGAL_PMP_REMESHING_VERBOSE t.stop(); - std::cout << "Remeshing done (size = " << target_edge_length; - std::cout << ", #iter = " << nb_iterations; + std::cout << "Remeshing done (#iter = " << nb_iterations; std::cout << ", " << t.time() << " sec )." << std::endl; #endif } +// Overload when using target_edge_length for sizing +template>> +void isotropic_remeshing(const FaceRange& faces + , const SizingValue target_edge_length + , PolygonMesh& pmesh + , const NamedParameters& np = parameters::default_values()) +{ + using parameters::choose_parameter; + using parameters::get_parameter; + + typedef typename GetVertexPointMap::type VPMap; + VPMap vpmap = choose_parameter(get_parameter(np, internal_np::vertex_point), + get_property_map(vertex_point, pmesh)); + + Uniform_sizing_field sizing(target_edge_length, vpmap); + if (target_edge_length > 0) + isotropic_remeshing(faces, sizing, pmesh, np); + else + isotropic_remeshing(faces, sizing, pmesh, np.do_split(false).do_collapse(false)); +} + /*! * \ingroup PMP_meshing_grp * @brief splits the edges listed in `edges` into sub-edges @@ -352,12 +376,13 @@ void isotropic_remeshing(const FaceRange& faces * has an internal property map for `CGAL::vertex_point_t`. * @tparam EdgeRange range of `boost::graph_traits::%edge_descriptor`, * model of `Range`. Its iterator type is `InputIterator`. +* @tparam SizingFunction model of `PMPSizingField` * @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters" * * @param pmesh a polygon mesh * @param edges the range of edges to be split if they are longer than given threshold -* @param max_length the edge length above which an edge from `edges` is split -* into to sub-edges +* @param sizing the sizing function that is used to split edges from 'edges' list. If a number convertible to +* a `double` is passed, `Uniform_sizing_field()` will be used, with the number as target edge length. * @param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below * \cgalNamedParamsBegin @@ -402,9 +427,11 @@ void isotropic_remeshing(const FaceRange& faces */ template + , typename SizingFunction + , typename NamedParameters = parameters::Default_named_parameters + , typename = typename std::enable_if_t>> void split_long_edges(const EdgeRange& edges - , const double& max_length + , SizingFunction& sizing , PolygonMesh& pmesh , const NamedParameters& np = parameters::default_values()) { @@ -451,7 +478,23 @@ void split_long_edges(const EdgeRange& edges fimap, false/*need aabb_tree*/); - remesher.split_long_edges(edges, max_length); + remesher.split_long_edges(edges, sizing); +} + +// Convenience overload when using max_length for sizing +template>> +void split_long_edges(const EdgeRange& edges + , const SizingValue max_length + , PolygonMesh& pmesh + , const NamedParameters& np = parameters::default_values()) +{ + // construct the uniform field, 3/4 as m_sq_long is stored with 4/3 of length + Uniform_sizing_field sizing(3. / 4. * max_length, pmesh); + split_long_edges(edges, sizing, pmesh, np); } } //end namespace Polygon_mesh_processing diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/tangential_relaxation.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/tangential_relaxation.h index fe680d007183..6693524324a9 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/tangential_relaxation.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/tangential_relaxation.h @@ -1,4 +1,4 @@ -// Copyright (c) 2015-2021 GeometryFactory (France). +// Copyright (c) 2015-2023 GeometryFactory (France). // All rights reserved. // // This file is part of CGAL (www.cgal.org). @@ -8,7 +8,7 @@ // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial // // -// Author(s) : Jane Tournois +// Author(s) : Jane Tournois, Ivan Paden #ifndef CGAL_POLYGON_MESH_PROCESSING_TANGENTIAL_RELAXATION_H @@ -17,6 +17,7 @@ #include #include +#include #include #include @@ -115,14 +116,23 @@ struct Allow_all_moves{ * \cgalParamDefault{If not provided, all moves are allowed.} * \cgalParamNEnd * +* \cgalParamNBegin{sizing_function} +* \cgalParamDescription{An object containing sizing field for individual vertices. +* Used to derive smoothing weights.} +* \cgalParamType{A model of `PMPSizingField`.} +* \cgalParamDefault{If not provided, smoothing weights are the same for all vertices.} +* \cgalParamNEnd +* * \cgalNamedParamsEnd * * \todo check if it should really be a triangle mesh or if a polygon mesh is fine */ -template +template void tangential_relaxation(const VertexRange& vertices, TriangleMesh& tm, - const NamedParameters& np = parameters::default_values()) + const CGAL_NP_CLASS& np = parameters::default_values()) { typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; @@ -130,18 +140,19 @@ void tangential_relaxation(const VertexRange& vertices, using parameters::get_parameter; using parameters::choose_parameter; + using parameters::is_default_parameter; - typedef typename GetGeomTraits::type GT; + typedef typename GetGeomTraits::type GT; GT gt = choose_parameter(get_parameter(np, internal_np::geom_traits), GT()); - typedef typename GetVertexPointMap::type VPMap; + typedef typename GetVertexPointMap::type VPMap; VPMap vpm = choose_parameter(get_parameter(np, internal_np::vertex_point), get_property_map(vertex_point, tm)); typedef Static_boolean_property_map Default_ECM; typedef typename internal_np::Lookup_named_param_def < internal_np::edge_is_constrained_t, - NamedParameters, + CGAL_NP_CLASS, Static_boolean_property_map // default (no constraint) > ::type ECM; ECM ecm = choose_parameter(get_parameter(np, internal_np::edge_is_constrained), @@ -149,12 +160,20 @@ void tangential_relaxation(const VertexRange& vertices, typedef typename internal_np::Lookup_named_param_def < internal_np::vertex_is_constrained_t, - NamedParameters, + CGAL_NP_CLASS, Static_boolean_property_map // default (no constraint) > ::type VCM; VCM vcm = choose_parameter(get_parameter(np, internal_np::vertex_is_constrained), Static_boolean_property_map()); + typedef typename internal_np::Lookup_named_param_def < + internal_np::sizing_function_t, + CGAL_NP_CLASS, + Uniform_sizing_field + > ::type SizingFunction; + SizingFunction sizing = choose_parameter(get_parameter(np, internal_np::sizing_function), + Uniform_sizing_field(-1, vpm)); + const bool relax_constraints = choose_parameter(get_parameter(np, internal_np::relax_constraints), false); const unsigned int nb_iterations = choose_parameter(get_parameter(np, internal_np::number_of_iterations), 1); @@ -201,7 +220,7 @@ void tangential_relaxation(const VertexRange& vertices, typedef typename internal_np::Lookup_named_param_def < internal_np::allow_move_functor_t, - NamedParameters, + CGAL_NP_CLASS, internal::Allow_all_moves// default > ::type Shall_move; Shall_move shall_move = choose_parameter(get_parameter(np, internal_np::allow_move_functor), @@ -244,15 +263,41 @@ void tangential_relaxation(const VertexRange& vertices, { const Vector_3& vn = vnormals.at(v); Vector_3 move = CGAL::NULL_VECTOR; - unsigned int star_size = 0; - for(halfedge_descriptor h :interior_hedges) + if constexpr (std::is_same_v>) { - move = move + Vector_3(get(vpm, v), get(vpm, source(h, tm))); - ++star_size; + unsigned int star_size = 0; + for(halfedge_descriptor h :interior_hedges) + { + move = move + Vector_3(get(vpm, v), get(vpm, source(h, tm))); + ++star_size; + } + CGAL_assertion(star_size > 0); //isolated vertices have already been discarded + move = (1. / static_cast(star_size)) * move; + } + else + { + auto gt_centroid = gt.construct_centroid_3_object(); + auto gt_area = gt.compute_area_3_object(); + double weight = 0; + for(halfedge_descriptor h :interior_hedges) + { + // calculate weight + // need v, v1 and v2 + const vertex_descriptor v1 = target(next(h, tm), tm); + const vertex_descriptor v2 = source(h, tm); + + const double tri_area = gt_area(get(vpm, v), get(vpm, v1), get(vpm, v2)); + const double face_weight = tri_area + / (1. / 3. * (sizing.at(v) + + sizing.at(v1) + + sizing.at(v2))); + weight += face_weight; + + const Point_3 centroid = gt_centroid(get(vpm, v), get(vpm, v1), get(vpm, v2)); + move = move + Vector_3(get(vpm, v), centroid) * face_weight; + } + move = move / weight; //todo ip: what if weight ends up being close to 0? } - CGAL_assertion(star_size > 0); //isolated vertices have already been discarded - move = (1. / static_cast(star_size)) * move; - barycenters.emplace_back(v, vn, get(vpm, v) + move); } else @@ -281,8 +326,6 @@ void tangential_relaxation(const VertexRange& vertices, bary = squared_distance(p1, bary) 0.03 //5 attempts maximum - && ( !check_normals(vp.first) - || !shall_move(vp.first, initial_pos, get(vpm, vp.first)))) //if a face has been inverted + && ( !check_normals(vp.first) + || !shall_move(vp.first, initial_pos, get(vpm, vp.first)))) //if a face has been inverted { frac = 0.5 * frac; put(vpm, vp.first, initial_pos + frac * move);//shorten the move by 2 @@ -335,7 +378,8 @@ void tangential_relaxation(const VertexRange& vertices, */ template -void tangential_relaxation(TriangleMesh& tm, const CGAL_NP_CLASS& np = parameters::default_values()) +void tangential_relaxation(TriangleMesh& tm, + const CGAL_NP_CLASS& np = parameters::default_values()) { tangential_relaxation(vertices(tm), tm, np); } diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt b/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt index 8ba77a7ff15e..ec6cdacd0dcb 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt @@ -84,6 +84,8 @@ if(TARGET CGAL::Eigen3_support) 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) + create_single_source_cgal_program("remeshing_quality_test.cpp" ) + target_link_libraries(remeshing_quality_test PUBLIC CGAL::Eigen3_support) else() message(STATUS "NOTICE: Tests that use the Eigen library will not be compiled.") endif() diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/remeshing_quality_test.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/remeshing_quality_test.cpp new file mode 100644 index 000000000000..7ab20cda4246 --- /dev/null +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/remeshing_quality_test.cpp @@ -0,0 +1,110 @@ +#include +#include +#include +#include +#include + +#include +#include + +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef K::Point_3 Point_3; +typedef CGAL::Surface_mesh Mesh; + +namespace PMP = CGAL::Polygon_mesh_processing; + +int main(int argc, char* argv[]) +{ + const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/hand.off"); + + Mesh mesh; + if (!PMP::IO::read_polygon_mesh(filename, mesh) || !CGAL::is_triangle_mesh(mesh)) + { + std::cerr << "Not a valid input file." << std::endl; + return 1; + } + + std::cout << "Start remeshing of " << filename + << " (" << num_faces(mesh) << " faces)..." << std::endl; + + const double tol = 0.001; + const std::pair edge_min_max{0.001, 0.5}; + PMP::Adaptive_sizing_field sizing_field(tol, edge_min_max, faces(mesh), mesh); + const unsigned int nb_iter = 5; + + PMP::isotropic_remeshing( + faces(mesh), + sizing_field, + mesh, + CGAL::parameters::number_of_iterations(nb_iter).number_of_relaxation_steps(3) + ); + + /* + * More information on quality metrics can be found here: https://ieeexplore.ieee.org/document/9167456 + */ + std::cout << "Remeshing done, checking triangle quality...\n" << std::endl; + double qmin = (std::numeric_limits::max)(); // minimum triangle quality + double qavg = 0.; // average quality + double min_angle = (std::numeric_limits::max)(); // minimum angle + double avg_min_angle = 0.; // average minimum angle + for (auto face : mesh.faces()) + { + if (PMP::is_degenerate_triangle_face(face, mesh)) + { + std::cout << "Found degenerate triangle!" << std::endl; + continue; + } + + // Calculate Q(t) triangle quality indicator + std::vector pts; pts.reserve(3); + for (auto vx : mesh.vertices_around_face(mesh.halfedge(face))) + pts.push_back(mesh.point(vx)); + + double half_perim = 0.; // half-perimeter + double lmax = 0.; // longest edge + for (int i = 0; i < 3; ++i) + { + const double length = CGAL::sqrt(CGAL::squared_distance(pts[i], pts[(i + 1) % 2])); + + half_perim += length; + if (length > lmax) lmax = length; + } + half_perim /= 2.; + const double area = CGAL::sqrt(CGAL::squared_area(pts[0], pts[1], pts[2])); + + const double face_quality = 6. / CGAL::sqrt(3.) * area / half_perim / lmax; + + qavg += face_quality; + if (face_quality < qmin) qmin = face_quality; + + // Calculate minimum angle + const auto v0 = pts[1] - pts[0]; + const auto v1 = pts[2] - pts[0]; + const auto v2 = pts[2] - pts[1]; + + const double dotp0 = CGAL::scalar_product(v0,v1); + const double angle0 = acos(dotp0 / (sqrt(v0.squared_length()) * sqrt(v1.squared_length()))); + const double dotp1 = CGAL::scalar_product(-v0, v2); + const double angle1 = acos(dotp1 / (sqrt(v0.squared_length()) * sqrt(v2.squared_length()))); + const double angle2 = CGAL_PI - (angle0 + angle1); + + double curr_min_angle = angle1; + if (angle1 < curr_min_angle) curr_min_angle = angle1; + if (angle2 < curr_min_angle) curr_min_angle = angle2; + + avg_min_angle += curr_min_angle; + if (curr_min_angle < min_angle) min_angle = curr_min_angle; + } + qavg /= mesh.number_of_faces(); + avg_min_angle /= mesh.number_of_faces(); + + std::cout << "Mesh size: " << mesh.number_of_faces() << std::endl; + std::cout << "Average quality: " << qavg << std::endl; + std::cout << "Minimum quality: " << qmin << std::endl; + std::cout << "Average minimum angle: " << avg_min_angle * 180. / CGAL_PI + << " deg" << std::endl; + std::cout << "Minimum angle: " << min_angle * 180. / CGAL_PI + << " deg" << std::endl; + + return 0; +} diff --git a/Polyhedron/demo/Polyhedron/Plugins/Display/Display_property.ui b/Polyhedron/demo/Polyhedron/Plugins/Display/Display_property.ui index 656ec4827c55..0a824308300d 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/Display/Display_property.ui +++ b/Polyhedron/demo/Polyhedron/Plugins/Display/Display_property.ui @@ -3,7 +3,7 @@ DisplayPropertyWidget - false + true @@ -39,9 +39,18 @@ + + true + Property + + false + + + false + 6 @@ -82,13 +91,114 @@ + + + + true + + + 0 + + + 100 + + + true + + + Qt::Horizontal + + + false + + + false + + + QSlider::TicksAbove + + + + + + + true + + + Expanding Radius: 0 + + + + + + + + + + false + + + Extreme Values + + + + + + Min Value + + + + + + + Max Value + + + + + + + 2.00 + + + true + + + + + + + Zoom to max value + + + + + + + Zoom to min value + + + + + + + 0.00 + + + true + + + true + + + - - - + + + Qt::Vertical @@ -100,7 +210,33 @@ - + + + + true + + + + + 0 + 0 + 236 + 335 + + + + + + + RAMP DISPLAYING + + + + + + + + Color Visualization @@ -118,10 +254,17 @@ 0 - - + + - Random colors + + + + + + + + Max color @@ -135,20 +278,6 @@ - - - - First color - - - - - - - Max color - - - @@ -159,32 +288,42 @@ - + - - + + - + Min color - - + + - Min color + First color + + + + + + + Random colors - - + + + + + Qt::Vertical @@ -196,97 +335,8 @@ - - - - true - - - - - 0 - 0 - 236 - 397 - - - - - - - RAMP DISPLAYING - - - - - - - - - - - false - - - Extreme Values - - - - - - Min Value - - - - - - - Max Value - - - - - - - 2.00 - - - true - - - - - - - Zoom to max value - - - - - - - Zoom to min value - - - - - - - 0.00 - - - true - - - true - - - - - - diff --git a/Polyhedron/demo/Polyhedron/Plugins/Display/Display_property_plugin.cpp b/Polyhedron/demo/Polyhedron/Plugins/Display/Display_property_plugin.cpp index 687cb6a92dfd..fb06449e972d 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/Display/Display_property_plugin.cpp +++ b/Polyhedron/demo/Polyhedron/Plugins/Display/Display_property_plugin.cpp @@ -28,11 +28,11 @@ #include #include #include -#include #include #include #include #include +#include #include #include #include @@ -47,9 +47,9 @@ #define ARBITRARY_DBL_MIN 1.0E-17 #define ARBITRARY_DBL_MAX 1.0E+17 -using namespace CGAL::Three; - +namespace PMP = CGAL::Polygon_mesh_processing; +using namespace CGAL::Three; Viewer_interface* (&getActiveViewer)() = Three::activeViewer; @@ -89,6 +89,9 @@ class Display_property_plugin double gI = 1.; double bI = 0.; + double expand_radius = 0.; + double maxEdgeLength = -1.; + Color_ramp color_ramp; std::vector color_map; QPixmap legend; @@ -102,10 +105,10 @@ class Display_property_plugin MIN_VALUE, MAX_VALUE }; - - enum CurvatureType { + enum CurvatureType + { MEAN_CURVATURE, - GAUSSIAN_CURVATURE + GAUSSIAN_CURVATURE, }; public: @@ -205,6 +208,9 @@ class Display_property_plugin this, &Display_property_plugin::on_zoomToMinButton_pressed); connect(dock_widget->zoomToMaxButton, &QPushButton::pressed, this, &Display_property_plugin::on_zoomToMaxButton_pressed); + + connect(dock_widget->expandingRadiusSlider, &QSlider::valueChanged, + this, &Display_property_plugin::setExpandingRadius); } private Q_SLOTS: @@ -247,26 +253,32 @@ private Q_SLOTS: dock_widget->maxBox->setRange(0, 360); dock_widget->maxBox->setValue(0); } - else if (property_name == "Interpolated Corrected Gaussian Curvature" || - property_name == "Interpolated Corrected Mean Curvature") + else if(property_name == "Scaled Jacobian") { dock_widget->minBox->setRange(-1000, 1000); dock_widget->minBox->setValue(0); dock_widget->maxBox->setRange(-1000, 1000); dock_widget->maxBox->setValue(0); } - else if(property_name == "Scaled Jacobian") + else if(property_name == "Face Area") { dock_widget->minBox->setRange(-1000, 1000); dock_widget->minBox->setValue(0); dock_widget->maxBox->setRange(-1000, 1000); dock_widget->maxBox->setValue(0); } - else if(property_name == "Face Area") + else if (property_name == "Interpolated Corrected Mean Curvature") { - dock_widget->minBox->setRange(-1000, 1000); + dock_widget->minBox->setRange(-99999999, 99999999); dock_widget->minBox->setValue(0); - dock_widget->maxBox->setRange(-1000, 1000); + dock_widget->maxBox->setRange(-99999999, 99999999); + dock_widget->maxBox->setValue(0); + } + else if (property_name == "Interpolated Corrected Gaussian Curvature") + { + dock_widget->minBox->setRange(-99999999, 99999999); + dock_widget->minBox->setValue(0); + dock_widget->maxBox->setRange(-99999999, 99999999); dock_widget->maxBox->setValue(0); } else @@ -500,6 +512,15 @@ private Q_SLOTS: { dock_widget->setEnabled(true); disableExtremeValues(); // only available after coloring + + // Curvature property-specific slider + const std::string& property_name = dock_widget->propertyBox->currentText().toStdString(); + const bool is_curvature_property = (property_name == "Interpolated Corrected Mean Curvature" || + property_name == "Interpolated Corrected Gaussian Curvature"); + dock_widget->expandingRadiusLabel->setVisible(is_curvature_property); + dock_widget->expandingRadiusSlider->setVisible(is_curvature_property); + dock_widget->expandingRadiusLabel->setEnabled(is_curvature_property); + dock_widget->expandingRadiusSlider->setEnabled(is_curvature_property); } else // no or broken property { @@ -527,15 +548,10 @@ private Q_SLOTS: { CGAL_assertion(static_cast(dock_widget->propertyBox->count()) == property_simplex_types.size()); - const int property_index = dock_widget->propertyBox->currentIndex(); - // leave it flat if it was, otherwise set to flat+edges - if(sm_item->renderingMode() != Flat && sm_item->renderingMode() != FlatPlusEdges && property_simplex_types.at(property_index) == Property_simplex_type::FACE) + if(sm_item->renderingMode() != Flat && sm_item->renderingMode() != FlatPlusEdges) sm_item->setRenderingMode(FlatPlusEdges); - if(sm_item->renderingMode() != Gouraud && sm_item->renderingMode() != GouraudPlusEdges && property_simplex_types.at(property_index) == Property_simplex_type::VERTEX) - sm_item->setRenderingMode(GouraudPlusEdges); - const std::string& property_name = dock_widget->propertyBox->currentText().toStdString(); if(property_name == "Smallest Angle Per Face") { @@ -555,11 +571,13 @@ private Q_SLOTS: } else if(property_name == "Interpolated Corrected Mean Curvature") { - displayCurvature(sm_item, MEAN_CURVATURE); + displayInterpolatedCurvatureMeasure(sm_item, MEAN_CURVATURE); + sm_item->setRenderingMode(Gouraud); } else if(property_name == "Interpolated Corrected Gaussian Curvature") { - displayCurvature(sm_item, GAUSSIAN_CURVATURE); + displayInterpolatedCurvatureMeasure(sm_item, GAUSSIAN_CURVATURE); + sm_item->setRenderingMode(Gouraud); } else { @@ -663,8 +681,8 @@ private Q_SLOTS: removeDisplayPluginProperty(item, "f:display_plugin_largest_angle"); removeDisplayPluginProperty(item, "f:display_plugin_scaled_jacobian"); removeDisplayPluginProperty(item, "f:display_plugin_area"); - removeDisplayPluginProperty(item, "f:display_plugin_mean_curvature"); - removeDisplayPluginProperty(item, "f:display_plugin_gaussian_curvature"); + removeDisplayPluginProperty(item, "v:display_plugin_interpolated_corrected_mean_curvature"); + removeDisplayPluginProperty(item, "v:display_plugin_interpolated_corrected_Gaussian_curvature"); } void displayExtremumAnglePerFace(Scene_surface_mesh_item* sm_item, @@ -753,7 +771,7 @@ private Q_SLOTS: halfedge_descriptor local_border_h = opposite(halfedge(local_f, local_smesh), local_smesh); CGAL_assertion(is_border(local_border_h, local_smesh)); - CGAL::Polygon_mesh_processing::triangulate_faces(local_smesh); + PMP::triangulate_faces(local_smesh); double extremum_angle_in_face = ARBITRARY_DBL_MAX; halfedge_descriptor local_border_end_h = local_border_h; @@ -806,7 +824,7 @@ private Q_SLOTS: const SMesh& mesh) const { if(CGAL::is_triangle(halfedge(f, mesh), mesh)) - return CGAL::Polygon_mesh_processing::face_area(f, mesh); + return PMP::face_area(f, mesh); auto vpm = get(boost::vertex_point, mesh); @@ -822,8 +840,8 @@ private Q_SLOTS: } CGAL::Euler::add_face(local_vertices, local_smesh); - CGAL::Polygon_mesh_processing::triangulate_faces(local_smesh); - return CGAL::Polygon_mesh_processing::area(local_smesh); + PMP::triangulate_faces(local_smesh); + return PMP::area(local_smesh); } void displayArea(Scene_surface_mesh_item* sm_item) @@ -845,30 +863,113 @@ private Q_SLOTS: displaySMProperty("f:display_plugin_area", *sm); } - void displayCurvature(Scene_surface_mesh_item* sm_item, - const CurvatureType curvature_type) +private Q_SLOTS: + void setExpandingRadius() { - SMesh* sm = sm_item->face_graph(); - if(sm == nullptr) + double sliderMin = dock_widget->expandingRadiusSlider->minimum(); + double sliderMax = dock_widget->expandingRadiusSlider->maximum() - sliderMin; + double val = dock_widget->expandingRadiusSlider->value() - sliderMin; + sliderMin = 0; + + Scene_item* item = scene->item(scene->mainSelectionIndex()); + Scene_surface_mesh_item* sm_item = qobject_cast(item); + if(sm_item == nullptr) return; - bool not_initialized; - std::string tied_string = (curvature_type == MEAN_CURVATURE) ? - "v:display_plugin_mean_curvature" : "v:display_plugin_gaussian_curvature"; + SMesh& smesh = *(sm_item->face_graph()); - SMesh::Property_map vcurvature; - std::tie(vcurvature, not_initialized) = sm->add_property_map(tied_string, 0); + auto vpm = get(CGAL::vertex_point, smesh); - if (curvature_type == MEAN_CURVATURE) + // @todo use the upcoming PMP::longest_edge + if(maxEdgeLength < 0) { - CGAL::Polygon_mesh_processing::interpolated_corrected_curvatures(*sm, CGAL::parameters::vertex_mean_curvature(vcurvature)); + auto edge_range = CGAL::edges(smesh); + + if(num_edges(smesh) == 0) + { + expand_radius = 0; + dock_widget->expandingRadiusLabel->setText(tr("Expanding Radius: %1").arg(expand_radius)); + return; + } + + auto eit = 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); + }); + + CGAL_assertion(eit != edge_range.end()); + + maxEdgeLength = PMP::edge_length(*eit, smesh); } - else if (curvature_type == GAUSSIAN_CURVATURE) + + 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)); + } + +private: + 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:display_plugin_interpolated_corrected_mean_curvature" + : "v:display_plugin_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) { - CGAL::Polygon_mesh_processing::interpolated_corrected_curvatures(*sm, CGAL::parameters::vertex_Gaussian_curvature(vcurvature)); + if(vnm_exists) + { + if(mu_index == MEAN_CURVATURE) + { + PMP::interpolated_corrected_curvatures(smesh, + CGAL::parameters::vertex_mean_curvature_map(mu_i_map) + .ball_radius(expand_radius) + .vertex_normal_map(vnm)); + } + else + { + PMP::interpolated_corrected_curvatures(smesh, + CGAL::parameters::vertex_Gaussian_curvature_map(mu_i_map) + .ball_radius(expand_radius) + .vertex_normal_map(vnm)); + } + } + else + { + if(mu_index == MEAN_CURVATURE) + { + PMP::interpolated_corrected_curvatures(smesh, + CGAL::parameters::vertex_mean_curvature_map(mu_i_map) + .ball_radius(expand_radius)); + } + else + { + PMP::interpolated_corrected_curvatures(smesh, + CGAL::parameters::vertex_Gaussian_curvature_map(mu_i_map) + .ball_radius(expand_radius)); + } + } } - displaySMProperty(tied_string, *sm); + displaySMProperty(tied_string, smesh); } private: @@ -1027,9 +1128,9 @@ private Q_SLOTS: else if(property_name == "Face Area") zoomToSimplexWithPropertyExtremum(faces(mesh), mesh, "f:display_plugin_area", extremum); else if(property_name == "Interpolated Corrected Mean Curvature") - zoomToSimplexWithPropertyExtremum(vertices(mesh), mesh, "v:display_plugin_mean_curvature", extremum); + zoomToSimplexWithPropertyExtremum(vertices(mesh), mesh, "v:display_plugin_interpolated_corrected_mean_curvature", extremum); else if(property_name == "Interpolated Corrected Gaussian Curvature") - zoomToSimplexWithPropertyExtremum(vertices(mesh), mesh, "v:display_plugin_gaussian_curvature", extremum); + zoomToSimplexWithPropertyExtremum(vertices(mesh), mesh, "v:display_plugin_interpolated_corrected_Gaussian_curvature", extremum); else if(property_simplex_types.at(property_index) == Property_simplex_type::VERTEX) zoomToSimplexWithPropertyExtremum(vertices(mesh), mesh, property_name, extremum); else if(property_simplex_types.at(property_index) == Property_simplex_type::FACE) @@ -1312,7 +1413,7 @@ scaled_jacobian(const face_descriptor f, for(std::size_t i=0; i 0 0 - 376 - 369 + 381 + 620 @@ -59,18 +59,29 @@ Isotropic remeshing - - - + + + - Preserve duplicated edges + - - Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter + + + + + + Minimum edge length + + + + + + + - + 0 @@ -83,24 +94,106 @@ - - + + + + Error tolerance + + + + + + + true + - - + + - Target edge length + + + + 0.00 + + + + + + + -1 + + + QLayout::SetDefaultConstraint + + + 11 + + + 14 + + + + + Edge sizing + + + + + + + + 0 + 0 + + + + + 168 + 16777215 + + + + + Uniform + + + + + Adaptive + + + + + + + + + + Qt::Vertical + + + + 20 + 20 + + + + + + + + Allow 1D smoothing along borders Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter - + Protect borders/selected edges @@ -110,17 +203,24 @@ - - + + - Number of Smoothing iterations + + + + + + + + Preserve duplicated edges Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter - + @@ -133,17 +233,7 @@ - - - - Allow 1D smoothing along borders - - - Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter - - - - + Number of Main iterations @@ -156,24 +246,7 @@ - - - - - - - - - - - - - - true - - - - + Qt::Vertical @@ -184,12 +257,22 @@ 20 - 40 + 24 - + + + + Target edge length + + + Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter + + + + Qt::ImhNone @@ -205,6 +288,67 @@ + + + + Maximum edge length + + + + + + + 0.00 + + + + + + + 0.00 + + + + + + + Number of Smoothing iterations + + + Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter + + + + + + + Ball radius + + + Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter + + + + + + + Curvature smoothing + + + Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter + + + + + + + -1 + + + false + + + diff --git a/Polyhedron/demo/Polyhedron/Plugins/PMP/Isotropic_remeshing_plugin.cpp b/Polyhedron/demo/Polyhedron/Plugins/PMP/Isotropic_remeshing_plugin.cpp index 53693c316379..c87fb567631d 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/PMP/Isotropic_remeshing_plugin.cpp +++ b/Polyhedron/demo/Polyhedron/Plugins/PMP/Isotropic_remeshing_plugin.cpp @@ -14,6 +14,7 @@ #include #include +#include #include #include #include @@ -293,9 +294,14 @@ class Polyhedron_demo_isotropic_remeshing_plugin : PMP::triangulate_face(f_and_p.first, *f_and_p.second); } - void do_split_edges(Scene_polyhedron_selection_item* selection_item, + void do_split_edges(const int edge_sizing_type, + Scene_polyhedron_selection_item* selection_item, SMesh& pmesh, - double target_length) + const double target_length, + const double error_tol, + const double min_length, + const double max_length, + const double curv_ball_r) { std::vector p_edges; for(edge_descriptor e : edges(pmesh)) @@ -313,12 +319,33 @@ class Polyhedron_demo_isotropic_remeshing_plugin : } } if (!p_edges.empty()) - CGAL::Polygon_mesh_processing::split_long_edges( - p_edges - , target_length - , *selection_item->polyhedron() - , CGAL::parameters::geom_traits(EPICK()) + { + if (edge_sizing_type == 0) + { + CGAL::Polygon_mesh_processing::split_long_edges( + p_edges + , target_length + , *selection_item->polyhedron() + , CGAL::parameters::geom_traits(EPICK()) + .edge_is_constrained_map(selection_item->constrained_edges_pmap())); + } + else if (edge_sizing_type == 1) + { + std::pair edge_min_max{min_length, max_length}; + CGAL::Polygon_mesh_processing::Adaptive_sizing_field adaptive_sizing( + error_tol + , edge_min_max + , faces(*selection_item->polyhedron()) + , *selection_item->polyhedron() + , CGAL::parameters::ball_radius(curv_ball_r)); + CGAL::Polygon_mesh_processing::split_long_edges( + p_edges + , adaptive_sizing + , *selection_item->polyhedron() + , CGAL::parameters::geom_traits(EPICK()) .edge_is_constrained_map(selection_item->constrained_edges_pmap())); + } + } else std::cout << "No selected or boundary edges to be split" << std::endl; } @@ -358,13 +385,18 @@ public Q_SLOTS: return; } + int edge_sizing_type = ui.edgeSizing_type_combo_box->currentIndex(); bool edges_only = ui.splitEdgesOnly_checkbox->isChecked(); bool preserve_duplicates = ui.preserveDuplicates_checkbox->isChecked(); double target_length = ui.edgeLength_dspinbox->value(); + double error_tol = ui.errorTol_edit->value(); + double min_length = ui.minEdgeLength_edit->value(); + double max_length = ui.maxEdgeLength_edit->value(); unsigned int nb_iter = ui.nbIterations_spinbox->value(); unsigned int nb_smooth = ui.nbSmoothing_spinbox->value(); bool protect = ui.protect_checkbox->isChecked(); bool smooth_features = ui.smooth1D_checkbox->isChecked(); + double curv_ball_r = ui.curvSmoothBallR_edit->value(); // wait cursor QApplication::setOverrideCursor(Qt::WaitCursor); @@ -396,7 +428,8 @@ public Q_SLOTS: { if (edges_only) { - do_split_edges(selection_item, pmesh, target_length); + do_split_edges(edge_sizing_type, selection_item, pmesh, + target_length, error_tol, min_length, max_length, curv_ball_r); } else //not edges_only { @@ -404,9 +437,8 @@ public Q_SLOTS: !CGAL::Polygon_mesh_processing::internal::constraints_are_short_enough( *selection_item->polyhedron(), selection_item->constrained_edges_pmap(), - get(CGAL::vertex_point, *selection_item->polyhedron()), CGAL::Constant_property_map(1), - 4. / 3. * target_length)) + CGAL::Polygon_mesh_processing::Uniform_sizing_field( 4. / 3. * target_length, pmesh))) { QApplication::restoreOverrideCursor(); //If facets are selected, splitting edges will add facets that won't be selected, and it will mess up the rest. @@ -432,7 +464,8 @@ public Q_SLOTS: } else { - do_split_edges(selection_item, pmesh, target_length); + do_split_edges(edge_sizing_type, selection_item, pmesh, + target_length, error_tol, min_length, max_length, curv_ball_r); } } @@ -457,28 +490,62 @@ public Q_SLOTS: } } - if (fpmap_valid) - CGAL::Polygon_mesh_processing::isotropic_remeshing(faces(*selection_item->polyhedron()) - , target_length - , *selection_item->polyhedron() - , CGAL::parameters::number_of_iterations(nb_iter) - .protect_constraints(protect) - .edge_is_constrained_map(selection_item->constrained_edges_pmap()) - .relax_constraints(smooth_features) - .number_of_relaxation_steps(nb_smooth) - .vertex_is_constrained_map(selection_item->constrained_vertices_pmap()) - .face_patch_map(fpmap)); - else - CGAL::Polygon_mesh_processing::isotropic_remeshing(faces(*selection_item->polyhedron()) - , target_length - , *selection_item->polyhedron() - , CGAL::parameters::number_of_iterations(nb_iter) - .protect_constraints(protect) - .edge_is_constrained_map(selection_item->constrained_edges_pmap()) - .relax_constraints(smooth_features) - .number_of_relaxation_steps(nb_smooth) - .vertex_is_constrained_map(selection_item->constrained_vertices_pmap()) - ); + if (edge_sizing_type == 0) + { + if (fpmap_valid) + CGAL::Polygon_mesh_processing::isotropic_remeshing(faces(*selection_item->polyhedron()) + , target_length + , *selection_item->polyhedron() + , CGAL::parameters::number_of_iterations(nb_iter) + .protect_constraints(protect) + .edge_is_constrained_map(selection_item->constrained_edges_pmap()) + .relax_constraints(smooth_features) + .number_of_relaxation_steps(nb_smooth) + .vertex_is_constrained_map(selection_item->constrained_vertices_pmap()) + .face_patch_map(fpmap)); + else + CGAL::Polygon_mesh_processing::isotropic_remeshing(faces(*selection_item->polyhedron()) + , target_length + , *selection_item->polyhedron() + , CGAL::parameters::number_of_iterations(nb_iter) + .protect_constraints(protect) + .edge_is_constrained_map(selection_item->constrained_edges_pmap()) + .relax_constraints(smooth_features) + .number_of_relaxation_steps(nb_smooth) + .vertex_is_constrained_map(selection_item->constrained_vertices_pmap()) + ); + } + else if (edge_sizing_type == 1) + { + std::pair edge_min_max{min_length, max_length}; + PMP::Adaptive_sizing_field adaptive_sizing_field(error_tol + , edge_min_max + , faces(*selection_item->polyhedron()) + , *selection_item->polyhedron() + , CGAL::parameters::ball_radius(curv_ball_r)); + if (fpmap_valid) + CGAL::Polygon_mesh_processing::isotropic_remeshing(faces(*selection_item->polyhedron()) + , adaptive_sizing_field + , *selection_item->polyhedron() + , CGAL::parameters::number_of_iterations(nb_iter) + .protect_constraints(protect) + .edge_is_constrained_map(selection_item->constrained_edges_pmap()) + .relax_constraints(smooth_features) + .number_of_relaxation_steps(nb_smooth) + .vertex_is_constrained_map(selection_item->constrained_vertices_pmap()) + .face_patch_map(fpmap)); + else + CGAL::Polygon_mesh_processing::isotropic_remeshing(faces(*selection_item->polyhedron()) + , adaptive_sizing_field + , *selection_item->polyhedron() + , CGAL::parameters::number_of_iterations(nb_iter) + .protect_constraints(protect) + .edge_is_constrained_map(selection_item->constrained_edges_pmap()) + .relax_constraints(smooth_features) + .number_of_relaxation_steps(nb_smooth) + .vertex_is_constrained_map(selection_item->constrained_vertices_pmap()) + ); + } } else //selected_facets not empty { @@ -507,27 +574,62 @@ public Q_SLOTS: } } - if (fpmap_valid) - CGAL::Polygon_mesh_processing::isotropic_remeshing(selection_item->selected_facets - , target_length - , *selection_item->polyhedron() - , CGAL::parameters::number_of_iterations(nb_iter) - .protect_constraints(protect) - .edge_is_constrained_map(selection_item->constrained_edges_pmap()) - .relax_constraints(smooth_features) - .number_of_relaxation_steps(nb_smooth) - .vertex_is_constrained_map(selection_item->constrained_vertices_pmap()) - .face_patch_map(fpmap)); - else - CGAL::Polygon_mesh_processing::isotropic_remeshing(selection_item->selected_facets - , target_length - , *selection_item->polyhedron() - , CGAL::parameters::number_of_iterations(nb_iter) - .protect_constraints(protect) - .edge_is_constrained_map(selection_item->constrained_edges_pmap()) - .relax_constraints(smooth_features) - .number_of_relaxation_steps(nb_smooth) - .vertex_is_constrained_map(selection_item->constrained_vertices_pmap())); + if (edge_sizing_type == 0) + { + if (fpmap_valid) + CGAL::Polygon_mesh_processing::isotropic_remeshing(selection_item->selected_facets + , target_length + , *selection_item->polyhedron() + , CGAL::parameters::number_of_iterations(nb_iter) + .protect_constraints(protect) + .edge_is_constrained_map(selection_item->constrained_edges_pmap()) + .relax_constraints(smooth_features) + .number_of_relaxation_steps(nb_smooth) + .vertex_is_constrained_map(selection_item->constrained_vertices_pmap()) + .face_patch_map(fpmap)); + else + CGAL::Polygon_mesh_processing::isotropic_remeshing(selection_item->selected_facets + , target_length + , *selection_item->polyhedron() + , CGAL::parameters::number_of_iterations(nb_iter) + .protect_constraints(protect) + .edge_is_constrained_map(selection_item->constrained_edges_pmap()) + .relax_constraints(smooth_features) + .number_of_relaxation_steps(nb_smooth) + .vertex_is_constrained_map(selection_item->constrained_vertices_pmap()) + ); + } + else if (edge_sizing_type == 1) + { + std::pair edge_min_max{min_length, max_length}; + PMP::Adaptive_sizing_field adaptive_sizing_field(error_tol + , edge_min_max + , faces(*selection_item->polyhedron()) + , *selection_item->polyhedron() + , CGAL::parameters::ball_radius(curv_ball_r)); + if (fpmap_valid) + CGAL::Polygon_mesh_processing::isotropic_remeshing(selection_item->selected_facets + , adaptive_sizing_field + , *selection_item->polyhedron() + , CGAL::parameters::number_of_iterations(nb_iter) + .protect_constraints(protect) + .edge_is_constrained_map(selection_item->constrained_edges_pmap()) + .relax_constraints(smooth_features) + .number_of_relaxation_steps(nb_smooth) + .vertex_is_constrained_map(selection_item->constrained_vertices_pmap()) + .face_patch_map(fpmap)); + else + CGAL::Polygon_mesh_processing::isotropic_remeshing(selection_item->selected_facets + , adaptive_sizing_field + , *selection_item->polyhedron() + , CGAL::parameters::number_of_iterations(nb_iter) + .protect_constraints(protect) + .edge_is_constrained_map(selection_item->constrained_edges_pmap()) + .relax_constraints(smooth_features) + .number_of_relaxation_steps(nb_smooth) + .vertex_is_constrained_map(selection_item->constrained_vertices_pmap()) + ); + } } } @@ -564,6 +666,9 @@ public Q_SLOTS: if (!edges_to_split.empty()) { if (fpmap_valid) + { + if (edge_sizing_type == 0) + { CGAL::Polygon_mesh_processing::split_long_edges( edges_to_split , target_length @@ -571,13 +676,51 @@ public Q_SLOTS: , CGAL::parameters::geom_traits(EPICK()) . edge_is_constrained_map(eif) . face_patch_map(fpmap)); + } + else if (edge_sizing_type == 1) + { + std::pair edge_min_max{min_length, max_length}; + PMP::Adaptive_sizing_field adaptive_sizing_field(error_tol + , edge_min_max + , faces(pmesh) + , pmesh + , CGAL::parameters::ball_radius(curv_ball_r)); + CGAL::Polygon_mesh_processing::split_long_edges( + edges_to_split + , target_length + , pmesh + , CGAL::parameters::geom_traits(EPICK()) + . edge_is_constrained_map(eif) + . face_patch_map(fpmap)); + } + } else + { + if (edge_sizing_type == 0) + { CGAL::Polygon_mesh_processing::split_long_edges( edges_to_split , target_length , pmesh , CGAL::parameters::geom_traits(EPICK()) . edge_is_constrained_map(eif)); + } + else if (edge_sizing_type == 1) + { + std::pair edge_min_max{min_length, max_length}; + PMP::Adaptive_sizing_field adaptive_sizing_field(error_tol + , edge_min_max + , faces(pmesh) + , pmesh + , CGAL::parameters::ball_radius(curv_ball_r)); + CGAL::Polygon_mesh_processing::split_long_edges( + edges_to_split + , adaptive_sizing_field + , pmesh + , CGAL::parameters::geom_traits(EPICK()) + . edge_is_constrained_map(eif)); + } + } } else std::cout << "No border to be split" << std::endl; @@ -604,9 +747,8 @@ public Q_SLOTS: !CGAL::Polygon_mesh_processing::internal::constraints_are_short_enough( pmesh, ecm, - get(CGAL::vertex_point, pmesh), CGAL::Constant_property_map(1), - 4. / 3. * target_length)) + CGAL::Polygon_mesh_processing::Uniform_sizing_field(4. / 3. * target_length, pmesh))) { QApplication::restoreOverrideCursor(); QMessageBox::warning(mw, tr("Error"), @@ -634,10 +776,42 @@ public Q_SLOTS: } } - if (fpmap_valid) + if (edge_sizing_type == 0) + { + if (fpmap_valid) + CGAL::Polygon_mesh_processing::isotropic_remeshing( + faces(*poly_item->polyhedron()) + , target_length + , *poly_item->polyhedron() + , CGAL::parameters::number_of_iterations(nb_iter) + .protect_constraints(protect) + .number_of_relaxation_steps(nb_smooth) + .edge_is_constrained_map(ecm) + .relax_constraints(smooth_features) + .face_patch_map(fpmap)); + else + CGAL::Polygon_mesh_processing::isotropic_remeshing( + faces(*poly_item->polyhedron()) + , target_length + , *poly_item->polyhedron() + , CGAL::parameters::number_of_iterations(nb_iter) + .protect_constraints(protect) + .number_of_relaxation_steps(nb_smooth) + .edge_is_constrained_map(ecm) + .relax_constraints(smooth_features)); + } + else if (edge_sizing_type == 1) + { + std::pair edge_min_max{min_length, max_length}; + PMP::Adaptive_sizing_field adaptive_sizing_field(error_tol + , edge_min_max + , faces(*poly_item->polyhedron()) + , *poly_item->polyhedron() + , CGAL::parameters::ball_radius(curv_ball_r)); + if (fpmap_valid) CGAL::Polygon_mesh_processing::isotropic_remeshing( faces(*poly_item->polyhedron()) - , target_length + , adaptive_sizing_field , *poly_item->polyhedron() , CGAL::parameters::number_of_iterations(nb_iter) .protect_constraints(protect) @@ -645,16 +819,17 @@ public Q_SLOTS: .edge_is_constrained_map(ecm) .relax_constraints(smooth_features) .face_patch_map(fpmap)); - else - CGAL::Polygon_mesh_processing::isotropic_remeshing( - faces(*poly_item->polyhedron()) - , target_length - , *poly_item->polyhedron() - , CGAL::parameters::number_of_iterations(nb_iter) - .protect_constraints(protect) - .number_of_relaxation_steps(nb_smooth) - .edge_is_constrained_map(ecm) - .relax_constraints(smooth_features)); + else + CGAL::Polygon_mesh_processing::isotropic_remeshing( + faces(*poly_item->polyhedron()) + , adaptive_sizing_field + , *poly_item->polyhedron() + , CGAL::parameters::number_of_iterations(nb_iter) + .protect_constraints(protect) + .number_of_relaxation_steps(nb_smooth) + .edge_is_constrained_map(ecm) + .relax_constraints(smooth_features)); + } //recollect sharp edges for(edge_descriptor e : edges(pmesh)) @@ -687,7 +862,11 @@ public Q_SLOTS: { // Remeshing parameters bool edges_only = false, preserve_duplicates = false; + int edge_sizing_type = 0; double target_length = 0.; + double error_tol = 0.; + double min_length = 0.; + double max_length = 0.; unsigned int nb_iter = 1; bool protect = false; bool smooth_features = true; @@ -723,7 +902,11 @@ public Q_SLOTS: edges_only = ui.splitEdgesOnly_checkbox->isChecked(); preserve_duplicates = ui.preserveDuplicates_checkbox->isChecked(); + edge_sizing_type = ui.edgeSizing_type_combo_box->currentIndex(); target_length = ui.edgeLength_dspinbox->value(); + error_tol = ui.errorTol_edit->value(); + min_length = ui.minEdgeLength_edit->value(); + max_length = ui.maxEdgeLength_edit->value(); nb_iter = ui.nbIterations_spinbox->value(); protect = ui.protect_checkbox->isChecked(); smooth_features = ui.smooth1D_checkbox->isChecked(); @@ -779,14 +962,18 @@ public Q_SLOTS: tbb::parallel_for( tbb::blocked_range(0, selection.size()), Remesh_polyhedron_item_for_parallel_for( - selection, edges_to_protect, edges_only, target_length, nb_iter, protect, smooth_features)); + selection, edges_to_protect, edges_only + , edge_sizing_type, target_length, error_tol + , min_length , max_length, nb_iter + , protect, smooth_features) + ); total_time = time.elapsed(); #else - Remesh_polyhedron_item remesher(edges_only, - target_length, nb_iter, protect, smooth_features); + Remesh_polyhedron_item remesher(edges_only, edge_sizing_type, + target_length, error_tol, min_length, max_length, nb_iter, protect, smooth_features); for(Scene_facegraph_item* poly_item : selection) { QElapsedTimer time; @@ -823,11 +1010,16 @@ public Q_SLOTS: typedef boost::graph_traits::halfedge_descriptor halfedge_descriptor; typedef boost::graph_traits::face_descriptor face_descriptor; + int edge_sizing_type_; bool edges_only_; double target_length_; + double error_tol_; + double min_length_; + double max_length_; unsigned int nb_iter_; bool protect_; bool smooth_features_; + double curv_ball_r_; protected: void remesh(Scene_facegraph_item* poly_item, @@ -846,25 +1038,62 @@ public Q_SLOTS: for(halfedge_descriptor h : border) border_edges.push_back(edge(h, *poly_item->polyhedron())); + if (edge_sizing_type_ == 0) + { CGAL::Polygon_mesh_processing::split_long_edges( border_edges , target_length_ , *poly_item->polyhedron()); + } + else if (edge_sizing_type_ == 1) + { + std::pair edge_min_max{min_length_, max_length_}; + PMP::Adaptive_sizing_field adaptive_sizing_field(error_tol_ + , edge_min_max + , faces(*poly_item->polyhedron()) + , *poly_item->polyhedron() + , CGAL::parameters::ball_radius(curv_ball_r_)); + CGAL::Polygon_mesh_processing::split_long_edges( + border_edges + , target_length_ + , *poly_item->polyhedron()); + } } else { std::cout << "Isotropic remeshing of " << poly_item->name().toStdString() << " started..." << std::endl; Scene_polyhedron_selection_item::Is_constrained_map ecm(&edges_to_protect); - CGAL::Polygon_mesh_processing::isotropic_remeshing( - faces(*poly_item->polyhedron()) - , target_length_ - , *poly_item->polyhedron() - , CGAL::parameters::number_of_iterations(nb_iter_) - .protect_constraints(protect_) - .edge_is_constrained_map(ecm) - .face_patch_map(get(CGAL::face_patch_id_t(), *poly_item->polyhedron())) - .relax_constraints(smooth_features_)); + if (edge_sizing_type_ == 0) + { + CGAL::Polygon_mesh_processing::isotropic_remeshing( + faces(*poly_item->polyhedron()) + , target_length_ + , *poly_item->polyhedron() + , CGAL::parameters::number_of_iterations(nb_iter_) + .protect_constraints(protect_) + .edge_is_constrained_map(ecm) + .face_patch_map(get(CGAL::face_patch_id_t(), *poly_item->polyhedron())) + .relax_constraints(smooth_features_)); + } + else + { + std::pair edge_min_max{min_length_, max_length_}; + PMP::Adaptive_sizing_field adaptive_sizing_field(error_tol_ + , edge_min_max + , faces(*poly_item->polyhedron()) + , *poly_item->polyhedron() + , CGAL::parameters::ball_radius(curv_ball_r_)); + CGAL::Polygon_mesh_processing::isotropic_remeshing( + faces(*poly_item->polyhedron()) + , target_length_ + , *poly_item->polyhedron() + , CGAL::parameters::number_of_iterations(nb_iter_) + .protect_constraints(protect_) + .edge_is_constrained_map(ecm) + .face_patch_map(get(CGAL::face_patch_id_t(), *poly_item->polyhedron())) + .relax_constraints(smooth_features_)); + } std::cout << "Isotropic remeshing of " << poly_item->name().toStdString() << " done." << std::endl; } @@ -873,20 +1102,32 @@ public Q_SLOTS: public: Remesh_polyhedron_item( const bool edges_only, + const int edge_sizing_type, const double target_length, + const double error_tol, + const double min_length, + const double max_length, const unsigned int nb_iter, const bool protect, const bool smooth_features) - : edges_only_(edges_only) + : edge_sizing_type_(edge_sizing_type) + , edges_only_(edges_only) , target_length_(target_length) + , error_tol_(error_tol) + , min_length_(min_length) + , max_length_(max_length) , nb_iter_(nb_iter) , protect_(protect) , smooth_features_(smooth_features) {} Remesh_polyhedron_item(const Remesh_polyhedron_item& remesh) - : edges_only_(remesh.edges_only_) + : edge_sizing_type_(remesh.edge_sizing_type_) + , edges_only_(remesh.edges_only_) , target_length_(remesh.target_length_) + , error_tol_(remesh.error_tol_) + , min_length_(remesh.min_length_) + , max_length_(remesh.max_length_) , nb_iter_(remesh.nb_iter_) , protect_(remesh.protect_) , smooth_features_(remesh.smooth_features_) @@ -913,11 +1154,17 @@ public Q_SLOTS: const std::vector& selection, std::map& edges_to_protect, const bool edges_only, + const int edge_sizing_type, const double target_length, + const double error_tol, + const double min_length, + const double max_length, const unsigned int nb_iter, const bool protect, const bool smooth_features) - : RemeshFunctor(edges_only, target_length, nb_iter, protect, smooth_features) + : RemeshFunctor(edges_only, edge_sizing_type, target_length + , error_tol, min_length, max_length + , nb_iter, protect, smooth_features) , selection_(selection), edges_to_protect_(edges_to_protect) {} @@ -986,6 +1233,55 @@ public Q_SLOTS: } } + void on_edgeSizing_type_combo_box_changed(int index) + { + if (index == 0) + { + ui.edgeLength_label->show(); + ui.edgeLength_dspinbox->show(); + ui.errorTol_label->hide(); + ui.errorTol_edit->hide(); + ui.minEdgeLength_label->hide(); + ui.minEdgeLength_edit->hide(); + ui.maxEdgeLength_label->hide(); + ui.maxEdgeLength_edit->hide(); + ui.curvSmooth_checkbox->hide(); + ui.curvSmooth_label->hide(); + ui.curvSmoothBallR_edit->hide(); + ui.curvSmoothBallR_label->hide(); + } + else if (index == 1) + { + ui.edgeLength_label->hide(); + ui.edgeLength_dspinbox->hide(); + ui.errorTol_label->show(); + ui.errorTol_edit->show(); + ui.minEdgeLength_label->show(); + ui.minEdgeLength_edit->show(); + ui.maxEdgeLength_label->show(); + ui.maxEdgeLength_edit->show(); + ui.curvSmooth_checkbox->show(); + ui.curvSmooth_label->show(); + ui.curvSmoothBallR_edit->show(); + ui.curvSmoothBallR_label->show(); + } + } + + void update_after_curvSmooth_click() + { + if (ui.curvSmooth_checkbox->isChecked()) + { + ui.curvSmoothBallR_label->setEnabled(true); + ui.curvSmoothBallR_edit->setEnabled(true); + } + else + { + ui.curvSmoothBallR_label->setEnabled(false); + ui.curvSmoothBallR_edit->setValue(-1); + ui.curvSmoothBallR_edit->setEnabled(false); + } + } + public: void initialize_remeshing_dialog(QDialog* dialog, @@ -1004,6 +1300,9 @@ public Q_SLOTS: connect(ui.protect_checkbox, SIGNAL(clicked(bool)), this, SLOT(update_after_protect_checkbox_click())); connect(ui.splitEdgesOnly_checkbox, SIGNAL(clicked(bool)), this, SLOT(update_after_splitEdgesOnly_click())); + connect(ui.edgeSizing_type_combo_box, SIGNAL(currentIndexChanged(int)), + this, SLOT(on_edgeSizing_type_combo_box_changed(int))); + connect(ui.curvSmooth_checkbox, SIGNAL(clicked(bool)), this, SLOT(update_after_curvSmooth_click())); //Set default parameters Scene_interface::Bbox bbox = poly_item != nullptr ? poly_item->bbox() @@ -1025,19 +1324,34 @@ public Q_SLOTS: ui.edgeLength_dspinbox->setValue(0.05 * diago_length); - - std::ostringstream oss; - oss << "Diagonal length of the Bbox of the selection to remesh is "; - oss << diago_length << "." << std::endl; - oss << "Default is 5% of it" << std::endl; - ui.edgeLength_dspinbox->setToolTip(QString::fromStdString(oss.str())); + ui.errorTol_edit->setValue(0.001 * diago_length); + ui.minEdgeLength_edit->setValue(0.001 * diago_length); + ui.maxEdgeLength_edit->setValue(0.5 * diago_length); + + std::string diag_general_info = "Diagonal length of the Bbox of the selection to remesh is " + + std::to_string(diago_length) + ".\n"; + std::string specific_info; + specific_info = "Default is 5% of it\n"; + ui.edgeLength_dspinbox->setToolTip(QString::fromStdString(diag_general_info + specific_info)); + specific_info = "Default is 0.1% of it\n"; + ui.errorTol_edit->setToolTip(QString::fromStdString(diag_general_info + specific_info)); + specific_info = "Default is 0.1% of it\n"; + ui.minEdgeLength_edit->setToolTip(QString::fromStdString(diag_general_info + specific_info)); + specific_info = "Default is 50% of it\n"; + ui.maxEdgeLength_edit->setToolTip(QString::fromStdString(diag_general_info + specific_info)); ui.nbIterations_spinbox->setSingleStep(1); ui.nbIterations_spinbox->setRange(1/*min*/, 1000/*max*/); ui.nbIterations_spinbox->setValue(1); + ui.edgeSizing_type_combo_box->setCurrentIndex(0); + on_edgeSizing_type_combo_box_changed(0); ui.protect_checkbox->setChecked(false); ui.smooth1D_checkbox->setChecked(true); + ui.curvSmooth_checkbox->setChecked(false); + ui.curvSmoothBallR_label->setEnabled(false); + ui.curvSmoothBallR_edit->setEnabled(false); + ui.curvSmoothBallR_edit->setValue(-1); if (nullptr != selection_item) { @@ -1047,7 +1361,6 @@ public Q_SLOTS: } } - private: QAction* actionIsotropicRemeshing_; Ui::Isotropic_remeshing_dialog ui; 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 17f9700713be..e22b053fe642 100644 --- a/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h +++ b/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h @@ -161,6 +161,7 @@ CGAL_add_named_parameter(vertex_corner_map_t, vertex_corner_map, vertex_corner_m CGAL_add_named_parameter(patch_normal_map_t, patch_normal_map, patch_normal_map) CGAL_add_named_parameter(region_primitive_map_t, region_primitive_map, region_primitive_map) CGAL_add_named_parameter(postprocess_regions_t, postprocess_regions, postprocess_regions) +CGAL_add_named_parameter(sizing_function_t, sizing_function, sizing_function) // List of named parameters that we use in the package 'Surface Mesh Simplification' CGAL_add_named_parameter(get_cost_policy_t, get_cost_policy, get_cost)