diff --git a/Documentation/doc/biblio/cgal_manual.bib b/Documentation/doc/biblio/cgal_manual.bib index d1600747b484..1330f605d52f 100644 --- a/Documentation/doc/biblio/cgal_manual.bib +++ b/Documentation/doc/biblio/cgal_manual.bib @@ -2386,6 +2386,46 @@ @article{ cgal:v-tm-95 ,update = "98.01 kettner" } +@article{cgal:vc-acvdupmc-04, + author = {Valette, Sébastien and Chassery, Jean-Marc}, + doi = {10.1111/j.1467-8659.2004.00769.x}, + hal_id = {hal-00534535}, + hal_version = {v1}, + journal = {Computer Graphics Forum}, + note = {Eurographics 2004 proceedings}, + number = {3}, + pages = {381-389}, + pdf = {https://hal.archives-ouvertes.fr/hal-00534535/file/valette.pdf}, + publisher = {Wiley}, + title = {Approximated Centroidal Voronoi Diagrams for Uniform Polygonal Mesh Coarsening}, + url = {https://hal.archives-ouvertes.fr/hal-00534535}, + volume = {23}, + year = {2004} +} + +@inproceedings{cgal:audette2011approach, + title={Approach-guided controlled resolution brain meshing for fe-based interactive neurosurgery simulation}, + author={Audette, M and Rivi{\`e}re, D and Ewend, M and Enquobahrie, A and Valette, S}, + booktitle={Workshop on Mesh Processing in Medical Image Analysis, in conjunction with MICCAI 2011}, + pages={176--186}, + year={2011} +} + +@article{cgal:vcp-grtmmdvd-08, + author = {Valette, Sébastien and Chassery, Jean-Marc and Prost, Rémy}, + doi = {10.1109/TVCG.2007.70430}, + hal_id = {hal-00537025}, + hal_version = {v1}, + journal = {IEEE Transactions on Visualization and Computer Graphics}, + number = {2}, + pages = {369--381}, + pdf = {https://hal.archives-ouvertes.fr/hal-00537025/file/VCP08TVCG.pdf}, + publisher = {Institute of Electrical and Electronics Engineers}, + title = {Generic remeshing of 3D triangular meshes with metric-dependent discrete Voronoi Diagrams}, + url = {https://hal.archives-ouvertes.fr/hal-00537025}, + volume = {14}, + year = {2008} +} @techreport{cgal:vla-lod-15, title={LOD Generation for Urban Scenes}, @@ -3111,14 +3151,12 @@ @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} +@inproceedings{dunyach2013curvRemesh, + title={Adaptive remeshing for real-time mesh deformation}, + author={Dunyach, Marion and Vanderhaeghe, David and Barthe, Loïc and Botsch, Mario}, + booktitle={Eurographics 2013}, + year={2013}, + organization={The Eurographics Association} } @book{botsch2010PMP, diff --git a/Installation/include/CGAL/license/Polygon_mesh_processing/acvd.h b/Installation/include/CGAL/license/Polygon_mesh_processing/acvd.h new file mode 100644 index 000000000000..eb65e39c3b80 --- /dev/null +++ b/Installation/include/CGAL/license/Polygon_mesh_processing/acvd.h @@ -0,0 +1,54 @@ +// Copyright (c) 2016 GeometryFactory SARL (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org) +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial +// +// Author(s) : Andreas Fabri +// +// Warning: this file is generated, see include/CGAL/license/README.md + +#ifndef CGAL_LICENSE_POLYGON_MESH_PROCESSING_ACVD_H +#define CGAL_LICENSE_POLYGON_MESH_PROCESSING_ACVD_H + +#include +#include + +#ifdef CGAL_POLYGON_MESH_PROCESSING_ACVD_COMMERCIAL_LICENSE + +# if CGAL_POLYGON_MESH_PROCESSING_ACVD_COMMERCIAL_LICENSE < CGAL_RELEASE_DATE + +# if defined(CGAL_LICENSE_WARNING) + + CGAL_pragma_warning("Your commercial license for CGAL does not cover " + "this release of the Polygon Mesh Processing - ACVD remeshing package.") +# endif + +# ifdef CGAL_LICENSE_ERROR +# error "Your commercial license for CGAL does not cover this release \ + of the Polygon Mesh Processing - ACVD remeshing package. \ + You get this error, as you defined CGAL_LICENSE_ERROR." +# endif // CGAL_LICENSE_ERROR + +# endif // CGAL_POLYGON_MESH_PROCESSING_ACVD_COMMERCIAL_LICENSE < CGAL_RELEASE_DATE + +#else // no CGAL_POLYGON_MESH_PROCESSING_ACVD_COMMERCIAL_LICENSE + +# if defined(CGAL_LICENSE_WARNING) + CGAL_pragma_warning("\nThe macro CGAL_POLYGON_MESH_PROCESSING_ACVD_COMMERCIAL_LICENSE is not defined." + "\nYou use the CGAL Polygon Mesh Processing - ACVD remeshing package under " + "the terms of the GPLv3+.") +# endif // CGAL_LICENSE_WARNING + +# ifdef CGAL_LICENSE_ERROR +# error "The macro CGAL_POLYGON_MESH_PROCESSING_ACVD_COMMERCIAL_LICENSE is not defined.\ + You use the CGAL Polygon Mesh Processing - ACVD remeshing package under the terms of \ + the GPLv3+. You get this error, as you defined CGAL_LICENSE_ERROR." +# endif // CGAL_LICENSE_ERROR + +#endif // no CGAL_POLYGON_MESH_PROCESSING_ACVD_COMMERCIAL_LICENSE + +#endif // CGAL_LICENSE_POLYGON_MESH_PROCESSING_ACVD_H diff --git a/Installation/include/CGAL/license/gpl_package_list.txt b/Installation/include/CGAL/license/gpl_package_list.txt index 75ff9d5c9ed9..91104d936f5e 100644 --- a/Installation/include/CGAL/license/gpl_package_list.txt +++ b/Installation/include/CGAL/license/gpl_package_list.txt @@ -67,6 +67,7 @@ Polygon_mesh_processing/miscellaneous Polygon Mesh Processing - Miscellaneous Polygon_mesh_processing/detect_features Polygon Mesh Processing - Feature Detection Polygon_mesh_processing/collision_detection Polygon Mesh Processing - Collision Detection Polygon_mesh_processing/autorefinement Polygon Mesh Processing - Autorefinement +Polygon_mesh_processing/acvd Polygon Mesh Processing - ACVD remeshing Polyhedron 3D Polyhedral Surface Polyline_simplification_2 2D Polyline Simplification Polytope_distance_d Optimal Distances diff --git a/Lab/demo/Lab/Plugins/PMP/ACVD_remeshing_dialog.ui b/Lab/demo/Lab/Plugins/PMP/ACVD_remeshing_dialog.ui new file mode 100644 index 000000000000..8c119c36ba6f --- /dev/null +++ b/Lab/demo/Lab/Plugins/PMP/ACVD_remeshing_dialog.ui @@ -0,0 +1,128 @@ + + + ACVDDialog + + + + 0 + 0 + 400 + 241 + + + + ACVD Remeshing + + + + + + Target number of Vertices (Clusters): + + + + + + + 100 + + + + + + + Cluster Metric: + + + + + + + Isotropic Clustering + + + true + + + buttonGroup + + + + + + + QEM Clustering + + + buttonGroup + + + + + + + Anisotropic Clustering + + + buttonGroup + + + + + + + Qt::Horizontal + + + QDialogButtonBox::Cancel|QDialogButtonBox::Ok + + + + + + + + DoubleEdit + QLineEdit +
CGAL_double_edit.h
+ + + + + + buttonBox + accepted() + ACVDDialog + accept() + + + 248 + 254 + + + 157 + 274 + + + + + buttonBox + rejected() + ACVDDialog + reject() + + + 316 + 260 + + + 286 + 274 + + + + + + + + diff --git a/Lab/demo/Lab/Plugins/PMP/ACVD_remeshing_plugin.cpp b/Lab/demo/Lab/Plugins/PMP/ACVD_remeshing_plugin.cpp new file mode 100644 index 000000000000..719586c710f2 --- /dev/null +++ b/Lab/demo/Lab/Plugins/PMP/ACVD_remeshing_plugin.cpp @@ -0,0 +1,124 @@ +#include + +#include "Messages_interface.h" +#include +#include +#include "Scene_polyhedron_selection_item.h" +#include "ui_ACVD_remeshing_dialog.h" + +#include "SMesh_type.h" +typedef Scene_surface_mesh_item Scene_facegraph_item; + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +typedef Scene_facegraph_item::Face_graph FaceGraph; + +using namespace CGAL::Three; +namespace PMP = CGAL::Polygon_mesh_processing; + +class CGAL_Lab_acvd_remeshing_plugin : + public QObject, + public CGAL_Lab_plugin_helper +{ + Q_OBJECT + Q_INTERFACES(CGAL::Three::CGAL_Lab_plugin_interface) + Q_PLUGIN_METADATA(IID "com.geometryfactory.PolyhedronDemo.PluginInterface/1.0" FILE "acvd_remeshing_plugin.json") +public: + bool applicable(QAction*) const { + return qobject_cast(scene->item(scene->mainSelectionIndex())); + } + void print_message(QString message) { CGAL::Three::Three::information(message);} + QList actions() const { return QList() << actionACVD; } + + + void init(QMainWindow* mainWindow, CGAL::Three::Scene_interface* scene_interface, Messages_interface*) { + mw = mainWindow; + scene = scene_interface; + actionACVD = new QAction(tr( + "ACVD Remeshing" + ), mw); + //actionACVD->setProperty("subMenuName", ""); + if (actionACVD) + connect(actionACVD, SIGNAL(triggered()), this, SLOT(acvd_action())); + } + +public Q_SLOTS: + void acvd_action() { + + // Create dialog box + QDialog dialog(mw); + ui.setupUi(&dialog); + + ui.nb_clusters_spin_box->setMinimum(10); + ui.nb_clusters_spin_box->setMaximum(1000000); + // get number of vertices of the mesh + Scene_facegraph_item* item = getSelectedItem(); + if(!item) { return; } + // set default number of clusters to 5% of the number of vertices + ui.nb_clusters_spin_box->setValue((std::max)(item->face_graph()->number_of_faces() / 20, (unsigned int)10)); + + int i = dialog.exec(); + if (i == QDialog::Rejected) + { + std::cout << "Remeshing aborted" << std::endl; + return; + } + + + Face_graph* graph = item->face_graph(); + if(!graph) return; + QElapsedTimer time; + time.start(); + CGAL::Three::Three::information("ACVD Remeshing..."); + + + QApplication::setOverrideCursor(Qt::WaitCursor); + FaceGraph remeshed = *graph; + + // TODO add post-processing + + if (ui.IsotropicClustering->isChecked()) + PMP::approximated_centroidal_Voronoi_diagram_remeshing(remeshed, ui.nb_clusters_spin_box->value()); + else if (ui.QEMClustering->isChecked()) + PMP::approximated_centroidal_Voronoi_diagram_remeshing(remeshed, ui.nb_clusters_spin_box->value(), CGAL::parameters::use_qem_based_energy(true)); + else + PMP::approximated_centroidal_Voronoi_diagram_remeshing(remeshed, ui.nb_clusters_spin_box->value(), CGAL::parameters::gradation_factor(0.8)); + + // update the scene + CGAL::Three::Three::information(QString("ok (%1 ms)").arg(time.elapsed())); + QApplication::restoreOverrideCursor(); + item->invalidateOpenGLBuffers(); + + // add the simplified mesh to the scene + Scene_facegraph_item* remeshed_item = new Scene_facegraph_item(remeshed); + remeshed_item->setName(QString(item->name() + "_simplified")); + scene->addItem(remeshed_item); + remeshed_item->invalidateOpenGLBuffers(); + scene->itemChanged(remeshed_item); + item->setVisible(false); + } + +private: + QAction* actionACVD; + Ui::ACVDDialog ui; +}; // end CGAL_Lab_acvd_remeshing_plugin + + +#include "ACVD_remeshing_plugin.moc" diff --git a/Lab/demo/Lab/Plugins/PMP/CMakeLists.txt b/Lab/demo/Lab/Plugins/PMP/CMakeLists.txt index f7efd5451286..d88e9f232f4c 100644 --- a/Lab/demo/Lab/Plugins/PMP/CMakeLists.txt +++ b/Lab/demo/Lab/Plugins/PMP/CMakeLists.txt @@ -29,6 +29,11 @@ target_link_libraries(extrude_plugin PRIVATE scene_surface_mesh_item if(TARGET CGAL::Eigen3_support) if("${EIGEN3_VERSION}" VERSION_GREATER "3.1.90") + + qt6_wrap_ui( acvd_remeshingUI_FILES ACVD_remeshing_dialog.ui) + cgal_lab_plugin(acvd_remeshing_plugin ACVD_remeshing_plugin ${acvd_remeshingUI_FILES} KEYWORDS PMP) + target_link_libraries(acvd_remeshing_plugin PUBLIC scene_surface_mesh_item scene_points_with_normal_item scene_polygon_soup_item CGAL::Eigen3_support) + qt6_wrap_ui( hole_fillingUI_FILES Hole_filling_widget.ui) cgal_lab_plugin(hole_filling_plugin Hole_filling_plugin ${hole_fillingUI_FILES} KEYWORDS PMP) target_link_libraries(hole_filling_plugin PRIVATE scene_surface_mesh_item scene_polylines_item scene_selection_item CGAL::Eigen3_support) diff --git a/Lab/demo/Lab/Plugins/Subdivision_methods/Subdivision_methods_plugin.cpp b/Lab/demo/Lab/Plugins/Subdivision_methods/Subdivision_methods_plugin.cpp index 8ec139d05589..538fca0de764 100644 --- a/Lab/demo/Lab/Plugins/Subdivision_methods/Subdivision_methods_plugin.cpp +++ b/Lab/demo/Lab/Plugins/Subdivision_methods/Subdivision_methods_plugin.cpp @@ -37,14 +37,17 @@ class CGAL_Lab_subdivision_methods_plugin : messages = m; scene = scene_interface; QAction *actionLoop = new QAction("Loop", mw); + QAction *actionUpsample = new QAction("Upsample", mw); QAction *actionCatmullClark = new QAction("Catmull Clark", mw); QAction *actionSqrt3 = new QAction("Sqrt3", mw); QAction *actionDooSabin = new QAction("Doo Sabin", mw); actionLoop->setObjectName("actionLoop"); + actionUpsample->setObjectName("actionUpsample"); actionCatmullClark->setObjectName("actionCatmullClark"); actionSqrt3->setObjectName("actionSqrt3"); actionDooSabin->setObjectName("actionDooSabin"); _actions << actionLoop + << actionUpsample << actionCatmullClark << actionSqrt3 << actionDooSabin; @@ -60,6 +63,7 @@ class CGAL_Lab_subdivision_methods_plugin : public Q_SLOTS: void on_actionLoop_triggered(); + void on_actionUpsample_triggered(); void on_actionCatmullClark_triggered(); void on_actionSqrt3_triggered(); void on_actionDooSabin_triggered(); @@ -69,6 +73,8 @@ private : template void apply_loop(FaceGraphItem* item, int nb_steps); template + void apply_upsample(FaceGraphItem* item, int nb_steps); + template void apply_catmullclark(FaceGraphItem* item, int nb_steps); template void apply_sqrt3(FaceGraphItem* item, int nb_steps); @@ -92,6 +98,22 @@ void CGAL_Lab_subdivision_methods_plugin::apply_loop(FaceGraphItem* item, int nb scene->itemChanged(item); } +template +void CGAL_Lab_subdivision_methods_plugin::apply_upsample(FaceGraphItem* item, int nb_steps) +{ + typename FaceGraphItem::Face_graph* graph = item->face_graph(); + if(!graph) return; + QElapsedTimer time; + time.start(); + CGAL::Three::Three::information("Upsample subdivision..."); + QApplication::setOverrideCursor(Qt::WaitCursor); + CGAL::Subdivision_method_3::Upsample_subdivision(*graph, params::number_of_iterations(nb_steps)); + CGAL::Three::Three::information(QString("ok (%1 ms)").arg(time.elapsed())); + QApplication::restoreOverrideCursor(); + item->invalidateOpenGLBuffers(); + scene->itemChanged(item); +} + void CGAL_Lab_subdivision_methods_plugin::on_actionLoop_triggered() { CGAL::Three::Scene_interface::Item_id index = scene->mainSelectionIndex(); @@ -111,6 +133,25 @@ void CGAL_Lab_subdivision_methods_plugin::on_actionLoop_triggered() apply_loop(sm_item, nb_steps); } +void CGAL_Lab_subdivision_methods_plugin::on_actionUpsample_triggered() +{ + CGAL::Three::Scene_interface::Item_id index = scene->mainSelectionIndex(); + Scene_surface_mesh_item* sm_item = qobject_cast(scene->item(index)); + if(!sm_item) + return; + + bool ok = true; + int nb_steps = QInputDialog::getInt(mw, + QString("Number of Iterations"), + QString("Choose number of iterations"), + 1 /* value */, 1 /* min */, (std::numeric_limits::max)() /* max */, 1 /*step*/, + &ok); + if(!ok) + return; + + apply_upsample(sm_item, nb_steps); +} + template void CGAL_Lab_subdivision_methods_plugin::apply_catmullclark(FaceGraphItem* item, int nb_steps) { diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt b/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt index 50f198e8317f..bf12fdb0cd10 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, Ivan Pađen, 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, Sébastien Valette, 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.} @@ -127,6 +127,7 @@ The page \ref bgl_namedparameters "Named Parameters" describes their usage. - \link PMP_meshing_grp `CGAL::Polygon_mesh_processing::isotropic_remeshing()` \endlink - \link PMP_meshing_grp `CGAL::Polygon_mesh_processing::surface_Delaunay_remeshing()` \endlink - \link PMP_meshing_grp `CGAL::Polygon_mesh_processing::split_long_edges()` \endlink +- `CGAL::Polygon_mesh_processing::approximated_centroidal_Voronoi_diagram_remeshing()` - `CGAL::Polygon_mesh_processing::extrude_mesh()` - `CGAL::Polygon_mesh_processing::smooth_mesh()` (deprecated) - `CGAL::Polygon_mesh_processing::angle_and_area_smoothing()` @@ -134,6 +135,7 @@ 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` 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 c4bf5a424fc8..158c2bcc40d2 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, Ivan Pađen, 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, Sébastien Valette, and Ilker %O. Yaz \image html neptun_head.jpg \image latex neptun_head.jpg @@ -107,7 +107,7 @@ maximizes the minimum angle of all the angles of the triangles in the triangulat \subsubsection Remeshing -\paragraph Local Isotropic Incremental Remeshing +\paragraph isorem Local Isotropic Incremental Remeshing The incremental triangle-based isotropic remeshing algorithm introduced by Botsch et al \cgalCite{botsch2004remeshing}, \cgalCite{botsch2010PMP} is implemented in this package. This algorithm incrementally performs simple operations such as edge splits, edge collapses, @@ -148,7 +148,7 @@ Sizing fields in isotropic remeshing. (b) Curvature-based adaptive sizing field. \cgalFigureEnd -\paragraph Delaunay-Based Surface Remeshing +\paragraph mesh3rem 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, generates a triangle surface mesh that provably has the required properties on simplices size, @@ -167,7 +167,32 @@ An example of how to remesh a given triangulated surface mesh with the Delaunay algorithm while preserving the detected sharp edges is given in \ref Polygon_mesh_processing/delaunay_remeshing_example.cpp -\paragraph Decimate Remeshing of (Almost) Planar Patches +\paragraph acvdrem Approximated Discrete Centroidal Voronoi Diagram (ACVD) Remeshing +This method is a vertex clustering based remeshing algorithm. It is based on +\cgalCite{cgal:vc-acvdupmc-04} and extended in \cgalCite{cgal:audette2011approach} and \cgalCite{cgal:vcp-grtmmdvd-08}. + +The function `CGAL::Polygon_mesh_processing::approximated_centroidal_Voronoi_diagram_remeshing()` takes as input a triangle mesh and +the expected number of output vertices, and it updates the mesh with the remeshed version. +The number of vertices in the output mesh might be larger than the input parameters if the input is not closed +or if the budget of points provided is too low to generate a manifold output mesh. +Note that providing an initial low number of vertices will affect the uniformity of the output triangles +in case some extra points are added to make the output manifold. The algorithm is similar to Lloyd's algorithm +(or k-means) where some random input vertices are picked to initialize clusters of vertices, and then clusters +are grown to minimize a particular energy, until convergence. Upon convergence, output vertices are computed +from the vertices of each cluster. +Note that there is no guarantee that the output mesh will have the same topology as the input. For example, if +a small handle is entirely contained in a cluster then this handle will not be part of the output. +If the input mesh contains some sharp features or corners, it is possible +to use the quadric error metrics to either move output vertices (fast but can produce bad looking triangles), +of to use the quadric error metrics directly in the energy formulation of each cluster (slower but produces +better quality triangles). It is also possible to use an adaptive remeshing based on curvature. + +\cgalFigureBegin{acvd_remeshing, acvd_remeshing.png} +ACVD remeshing. +From left to right: Input; with default parameters; using adaptative option; using QEM based energy. +\cgalFigureEnd + +\paragraph decimate Remeshing of (Almost) Planar Patches When many triangles are used to describe a planar region of a model, one might wish to simplify the mesh in this region to use few elements, or even a single large polygonal face when the region makes up a simply connected patch. This can be achieved using the function `CGAL::Polygon_mesh_processing::remesh_planar_patches()`. @@ -1450,5 +1475,11 @@ 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. +The implementation of the ACVD Remeshing method was initiated by Hossam Saeed during GSoC 2023, under +supervision of Sébastien Valette and Sébastien Loriot, who later finalized the work. +The implementation is based on \cgalCite{cgal:vcp-grtmmdvd-08}. and +preceeding work. ACVD's implementation was also +used as a reference during the project. + */ } /* namespace CGAL */ diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/examples.txt b/Polygon_mesh_processing/doc/Polygon_mesh_processing/examples.txt index dba765d7f6a1..f2005fa8b4fd 100644 --- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/examples.txt +++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/examples.txt @@ -50,6 +50,7 @@ \example Polygon_mesh_processing/cc_compatible_orientations.cpp \example Polygon_mesh_processing/remesh_planar_patches.cpp \example Polygon_mesh_processing/remesh_almost_planar_patches.cpp +\example Polygon_mesh_processing/acvd_remeshing_example.cpp \example Polygon_mesh_processing/sample_example.cpp \example Polygon_mesh_processing/soup_autorefinement.cpp */ diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/fig/acvd_remeshing.png b/Polygon_mesh_processing/doc/Polygon_mesh_processing/fig/acvd_remeshing.png new file mode 100644 index 000000000000..2bed8627c842 Binary files /dev/null and b/Polygon_mesh_processing/doc/Polygon_mesh_processing/fig/acvd_remeshing.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 cbeb71cb6b9a..f5035cd158c8 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt @@ -62,6 +62,8 @@ create_single_source_cgal_program("soup_autorefinement.cpp") find_package(Eigen3 3.2.0 QUIET) #(requires 3.2.0 or greater) include(CGAL_Eigen3_support) if(TARGET CGAL::Eigen3_support) + create_single_source_cgal_program("acvd_remeshing_example.cpp") + target_link_libraries(acvd_remeshing_example PUBLIC CGAL::Eigen3_support) create_single_source_cgal_program("hole_filling_example.cpp") target_link_libraries(hole_filling_example PRIVATE CGAL::Eigen3_support) create_single_source_cgal_program("hole_filling_example_SM.cpp") diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/acvd_remeshing_example.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/acvd_remeshing_example.cpp new file mode 100644 index 000000000000..ca7a965d9945 --- /dev/null +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/acvd_remeshing_example.cpp @@ -0,0 +1,88 @@ +#include +#include +#include +#include +#include + +#include + +#include +#include + +namespace PMP = CGAL::Polygon_mesh_processing; + +typedef CGAL::Exact_predicates_inexact_constructions_kernel Epic_kernel; +typedef CGAL::Surface_mesh Mesh; +typedef boost::graph_traits::vertex_descriptor vertex_descriptor; + +namespace params = CGAL::parameters; + +int main(int argc, char* argv[]) +{ + Mesh smesh; + //CGAL::get_default_random() = CGAL::Random( 1739197120 ); //connexity constraint issue + infinite loop with cheese + //CGAL::get_default_random() = CGAL::Random( 1739199762 ); //one small edge + //CGAL::get_default_random() = CGAL::Random( 1739264620 ); //one very small edge + //CGAL::get_default_random() = CGAL::Random( 1739293586 ); // seed elephant non-manifold with 300 clusters + regression for fandisk 3000 with qem + + std::cout << "Seed : " << CGAL::get_default_random().get_seed() << std::endl; + const std::string filename = (argc > 1) ? + argv[1] : + CGAL::data_file_path("meshes/fandisk.off"); + + const std::string stem = std::filesystem::path(filename).stem().string(); + const std::string extension = std::filesystem::path(filename).extension().string(); + + const int nb_clusters = (argc > 2) ? atoi(argv[2]) : 3000; + const std::string nbc = std::to_string(nb_clusters); + + if (!CGAL::IO::read_polygon_mesh(filename, smesh)) + { + std::cerr << "Invalid input file." << std::endl; + return EXIT_FAILURE; + } + + ///// Uniform Isotropic ACVD +#if 0 + std::cout << "Uniform Isotropic ACVD ...." << std::endl; + Mesh acvd_mesh = smesh; + PMP::approximated_centroidal_Voronoi_diagram_remeshing(acvd_mesh, nb_clusters); + CGAL::IO::write_polygon_mesh(stem+"_acvd_"+nbc+extension, acvd_mesh); + + std::cout << "Completed" << std::endl; +#endif + + //// With Post-Processing QEM Optimization +#if 0 + std::cout << "Uniform Isotropic ACVD with QEM optimization ...." << std::endl; + + Mesh acvd_mesh_qem_pp = smesh; + PMP::approximated_centroidal_Voronoi_diagram_remeshing(acvd_mesh_qem_pp, nb_clusters, params::use_postprocessing_qem(true)); + CGAL::IO::write_polygon_mesh( stem +"_acvd_qem-pp_"+nbc+extension, acvd_mesh_qem_pp); + + std::cout << "Completed" << std::endl; +#endif + +#if 1 + // With QEM Energy Minimization + std::cout << "Uniform QEM ACVD ...." << std::endl; + auto acvd_mesh_qem = smesh; + PMP::approximated_centroidal_Voronoi_diagram_remeshing(acvd_mesh_qem, nb_clusters, params::use_qem_based_energy(true)); + CGAL::IO::write_polygon_mesh(stem +"_acvd_qem_"+ std::to_string(nb_clusters) + extension, acvd_mesh_qem); +#endif + +#if 0 + /// Adaptive Isotropic ACVD + std::cout << "Adaptive Isotropic ACVD ...." << std::endl; + const double gradation_factor = 2; + Mesh adaptive_acvd_mesh = smesh; + PMP::approximated_centroidal_Voronoi_diagram_remeshing(adaptive_acvd_mesh, nb_clusters, params::gradation_factor(gradation_factor)); + + CGAL::IO::write_OFF(stem +"_acvd_adaptative_"+ std::to_string(nb_clusters) + extension, acvd_mesh_qem, adaptive_acvd_mesh); +#endif + + std::cout << "Completed" << std::endl; + + return 0; +} + diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/approximated_centroidal_Voronoi_diagram_remeshing.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/approximated_centroidal_Voronoi_diagram_remeshing.h new file mode 100644 index 000000000000..e5541214dd5a --- /dev/null +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/approximated_centroidal_Voronoi_diagram_remeshing.h @@ -0,0 +1,1340 @@ +// 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) : Hossam Saeed, Sébastien Valette, and Sébastien Loriot +// + +#ifndef CGAL_PMP_ACVD_REMESHING_H +#define CGAL_PMP_ACVD_REMESHING_H + +#include + +#include +#include +#include +#include + +#include +#include +#ifdef CGAL_DEBUG_ACVD +#include +#endif +#include +#include +#include + +#ifndef CGAL_ACVD_DOES_NOT_USE_INTERPOLATED_CORRECTED_CURVATURES +#include +#endif +#include +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include + + +#define CGAL_TO_QEM_MODIFICATION_THRESHOLD 1e-3 +#define CGAL_WEIGHT_CLAMP_RATIO_THRESHOLD 10000 + +namespace CGAL { + +namespace Polygon_mesh_processing { + +namespace internal { +#ifdef CGAL_DEBUG_ACVD +template +void dump_mesh_with_cluster_colors(TriangleMesh pmesh, ClusterMap cluster_map, std::string fname) +{ + std::vector palette {{ CGAL::IO::red(), + CGAL::IO::green(), + CGAL::IO::blue(), + CGAL::IO::purple(), + CGAL::IO::orange(), + CGAL::IO::deep_blue(), + CGAL::IO::yellow(), + CGAL::IO::violet(), + CGAL::IO::gray(), + CGAL::IO::white() }}; + + auto vcm = pmesh.template add_property_map("f:color").first; + + for (auto v : vertices(pmesh)) + { + int cluster_id=get(cluster_map, v); + if (cluster_id!=-1) + put(vcm, v, palette[ cluster_id % palette.size() ]); + else + put(vcm, v, CGAL::IO::black()); + } + + std::ofstream out(fname); + CGAL::IO::write_PLY(out, pmesh, CGAL::parameters::use_binary_mode(false) + .stream_precision(17) + .vertex_color_map(vcm)); + +} +#endif + +template +void compute_qem_face(const typename GT::Vector_3& p1, const typename GT::Vector_3& p2, const typename GT::Vector_3& p3, Eigen::Matrix& quadric) +{ + typename GT::Construct_cross_product_vector_3 cross_product; + + typename GT::Vector_3 crossX1X2 = cross_product(p1, p2); + typename GT::Vector_3 crossX2X3 = cross_product(p2, p3); + typename GT::Vector_3 crossX3X1 = cross_product(p3, p1); + typename GT::FT determinantABC = CGAL::determinant(p1, p2, p3); + + typename GT::FT n[4] = { + crossX1X2.x() + crossX2X3.x() + crossX3X1.x(), + crossX1X2.y() + crossX2X3.y() + crossX3X1.y(), + crossX1X2.z() + crossX2X3.z() + crossX3X1.z(), + -determinantABC + }; + + for (int i = 0; i < 4; ++i) + for (int j = 0; j < 4; ++j) + quadric(i, j) = n[i] * n[j]; +} + +template +void compute_qem_vertex(std::vector> cluster_tris, Eigen::Matrix& quadric) +{ + quadric.setZero(); + + for (int i = 0; i < cluster_tris.size(); ++i) + { + Eigen::Matrix q; + compute_qem_face(cluster_tris[i][0], cluster_tris[i][1], cluster_tris[i][2], q); + quadric = quadric + q; + } +} + +template +typename GT::Vector_3 compute_displacement(const Eigen::Matrix quadric, const typename GT::Point_3& p, int& rank_deficiency) +{ + using Matrix3d = Eigen::Matrix; + + int max_nb_of_singular_values_used = 3; + Matrix3d A; + A(0, 0) = quadric(0, 0); + A(0, 1) = A(1, 0) = quadric(0, 1); + A(0, 2) = A(2, 0) = quadric(0, 2); + A(1, 1) = quadric(1, 1); + A(1, 2) = A(2, 1) = quadric(1, 2); + A(2, 2) = quadric(2, 2); + + Matrix3d U, V; + typename GT::Vector_3 w; + Eigen::JacobiSVD svd(A, Eigen::ComputeFullU | Eigen::ComputeFullV); + U = svd.matrixU(); + V = svd.matrixV(); + auto w_temp = svd.singularValues(); + w = typename GT::Vector_3(w_temp(0), w_temp(1), w_temp(2)); + + // compute all eigen values absolute values + typename GT::FT AbsolutesEigenValues[3]; + typename GT::FT maxW = -1.0; + for (int j = 0; j < 3; ++j) + { + typename GT::FT AbsoluteEigenValue = fabs(w[j]); + AbsolutesEigenValues[j] = AbsoluteEigenValue; + if (AbsoluteEigenValue > maxW) + maxW = AbsoluteEigenValue; + } + typename GT::FT invmaxW = 1.0 / maxW; + + Matrix3d tempMatrix, tempMatrix2; + + for (int i = 0; i < 3; ++i) + { + typename GT::FT LocalMaxW = -1; + int IndexMax = -1; + + // find the remaining eigenvalue with highest absolute value + for (int j = 0; j < 3; ++j) + { + if (LocalMaxW < AbsolutesEigenValues[j]) + { + LocalMaxW = AbsolutesEigenValues[j]; + IndexMax = j; + } + } + + if ((AbsolutesEigenValues[IndexMax] * invmaxW > 1e-3) + && (max_nb_of_singular_values_used > 0)) + { + // If this is true, then w[i] != 0, so this division is ok. + double Inv = 1.0 / w[IndexMax]; + tempMatrix(IndexMax, 0) = U(0, IndexMax) * Inv; + tempMatrix(IndexMax, 1) = U(1, IndexMax) * Inv; + tempMatrix(IndexMax, 2) = U(2, IndexMax) * Inv; + } + else + { + tempMatrix(IndexMax, 0) = 0; + tempMatrix(IndexMax, 1) = 0; + tempMatrix(IndexMax, 2) = 0; + ++rank_deficiency; + } + + // set the eigenvalu to -2 to remove it from subsequent tests + AbsolutesEigenValues[IndexMax] = -2; + --max_nb_of_singular_values_used; + } + + tempMatrix2 = V * tempMatrix; + Eigen::Matrix p_temp, b; + b << -quadric(0, 3), -quadric(1, 3), -quadric(2, 3); + p_temp << p.x(), p.y(), p.z(); + Eigen::Matrix displacement_temp = tempMatrix2 * (b - A * p_temp); + typename GT::Vector_3 displacement = { displacement_temp(0), displacement_temp(1), displacement_temp(2) }; + return displacement; +} + +template +void compute_representative_point(const Eigen::Matrix& quadric, typename GT::Point_3& p, int& rank_deficiency) +{ + // average point + typename GT::Vector_3 displacement = compute_displacement(quadric, p, rank_deficiency); + p = p + displacement; +} + + +template +struct QEMClusterData { + typename GT::Vector_3 site_sum; + typename GT::FT weight_sum; + Eigen::Matrix quadric_sum; + typename GT::Vector_3 representative_point_; + typename GT::FT energy; + bool modified = true; + int last_modification_iteration = 0; + size_t nb_vertices = 0; + + QEMClusterData() : + site_sum(0, 0, 0), + weight_sum(0), + representative_point_(0, 0, 0), + energy(0) + { + quadric_sum.setZero(); + } + + void add_vertex(const typename GT::Vector_3 vertex_position, const typename GT::FT weight, const Eigen::Matrix& quadric) + { + this->site_sum += vertex_position * weight; + this->weight_sum += weight; + this->quadric_sum += quadric; + ++this->nb_vertices; + this->modified=true; + } + + void add_vertex(const typename GT::Point_3 vertex_position, const typename GT::FT weight, const Eigen::Matrix& quadric) + { + add_vertex( vertex_position - ORIGIN, weight, quadric ); + } + + void remove_vertex(const typename GT::Vector_3 vertex_position, const typename GT::FT weight, const Eigen::Matrix& quadric) + { + this->site_sum -= vertex_position * weight; + this->weight_sum -= weight; + this->quadric_sum -= quadric; + this->nb_vertices--; + this->modified=true; + } + + void compute_representative(bool qem) + { + if (this->weight_sum > 0) + { + if (qem) + { + int rank_deficiency = 0; + typename GT::Point_3 point = CGAL::ORIGIN + (this->site_sum / this->weight_sum); + compute_representative_point(this->quadric_sum, point, rank_deficiency); + this->representative_point_ = { point.x(), point.y(), point.z() }; + } + else + this->representative_point_ = (this->site_sum) / this->weight_sum; + } + } + + typename GT::FT compute_energy(bool qem) + { + if (modified) + { + if (qem) + { + compute_representative(qem); + auto dot_product = GT().compute_scalar_product_3_object(); + + this->energy = (this->representative_point_).squared_length() * this->weight_sum + - 2 * dot_product(this->representative_point_, this->site_sum); + } + else + this->energy = - (this->site_sum).squared_length() / this->weight_sum; + } + return this->energy; + } + + typename GT::Point_3 representative_point(bool qem) + { + if (modified) compute_representative(qem); + return ORIGIN+representative_point_; + } +}; + +#ifndef CGAL_ACVD_DOES_NOT_USE_INTERPOLATED_CORRECTED_CURVATURES +template +void upsample_subdivision_property(TriangleMesh& pmesh, VPCDM vpcd_map, const NamedParameters& np) { + using GT = typename GetGeomTraits::type; + using Vertex_descriptor = typename boost::graph_traits::vertex_descriptor; + using Halfedge_descriptor = typename boost::graph_traits::halfedge_descriptor; + + using parameters::choose_parameter; + using parameters::get_parameter; + using parameters::is_default_parameter; + + using VPM = typename CGAL::GetVertexPointMap::type; + VPM vpm = choose_parameter(get_parameter(np, internal_np::vertex_point), + get_property_map(CGAL::vertex_point, pmesh)); + + // unordered_set of old vertices + std::unordered_set old_vertices; + + unsigned int step = choose_parameter(get_parameter(np, internal_np::number_of_iterations), 1); + Upsample_mask_3 mask(&pmesh, vpm); + + for (unsigned int i = 0; i < step; ++i){ + for (Vertex_descriptor vd : vertices(pmesh)) + old_vertices.insert(vd); + + Subdivision_method_3::internal::PTQ_1step(pmesh, vpm, mask); + // interpolate curvature values + for (Vertex_descriptor vd : vertices(pmesh)) + { + if (old_vertices.find(vd) == old_vertices.end()) + { + Principal_curvatures_and_directions pcd; + pcd.min_curvature = 0; + pcd.max_curvature = 0; + pcd.min_direction = typename GT::Vector_3(0, 0, 0); + pcd.max_direction = typename GT::Vector_3(0, 0, 0); + for (Halfedge_descriptor hd : halfedges_around_target(vd, pmesh)) + { + Vertex_descriptor v1 = source(hd, pmesh); + if (old_vertices.find(v1) != old_vertices.end()) + { + Principal_curvatures_and_directions pcd1 = get(vpcd_map, v1); + pcd.min_curvature += pcd1.min_curvature; + pcd.max_curvature += pcd1.max_curvature; + pcd.min_direction += pcd1.min_direction; + pcd.max_direction += pcd1.max_direction; + } + } + pcd.min_curvature = pcd.min_curvature / 2; + pcd.max_curvature = pcd.max_curvature / 2; + pcd.min_direction = pcd.min_direction / sqrt(pcd.min_direction.squared_length()); + pcd.max_direction = pcd.max_direction / sqrt(pcd.max_direction.squared_length()); + put(vpcd_map, vd, pcd); + } + } + } +} +#endif + +template +std::pair< + std::vector::type::Point_3>, + std::vector> +> acvd_impl( + TriangleMesh& pmesh, + int nb_clusters, + const NamedParameters& np = parameters::default_values() + ) +{ + using GT = typename GetGeomTraits::type; + using Matrix4x4 = typename Eigen::Matrix; + using Vertex_position_map = typename GetVertexPointMap::const_type; + using Halfedge_descriptor = typename boost::graph_traits::halfedge_descriptor; + using Edge_descriptor = typename boost::graph_traits::edge_descriptor; + using Vertex_descriptor = typename boost::graph_traits::vertex_descriptor; + using Face_descriptor = typename boost::graph_traits::face_descriptor; + using VertexClusterMap = typename boost::property_map >::type; + using VertexVisitedMap = typename boost::property_map >::type; + using VertexWeightMap = typename boost::property_map >::type; + using VertexQuadricMap = typename boost::property_map >::type; +#ifndef CGAL_ACVD_DOES_NOT_USE_INTERPOLATED_CORRECTED_CURVATURES + using Default_principal_map = typename boost::property_map> >::type; + using Vertex_principal_curvatures_and_directions_map = + typename internal_np::Lookup_named_param_def::type ; +#endif + using parameters::choose_parameter; + using parameters::get_parameter; + using parameters::is_default_parameter; + + Vertex_position_map vpm = choose_parameter(get_parameter(np, CGAL::vertex_point), + get_property_map(CGAL::vertex_point, pmesh)); + + // get curvature related parameters +#ifdef CGAL_ACVD_DOES_NOT_USE_INTERPOLATED_CORRECTED_CURVATURES + static_assert(is_default_parameter::value, "gradation_factor option is disabled"); +#endif + const typename GT::FT gradation_factor = choose_parameter(get_parameter(np, internal_np::gradation_factor), 0); + + // using QEM ? + bool use_postprocessing_qem = choose_parameter(get_parameter(np, internal_np::use_postprocessing_qem), false); + const bool use_qem_based_energy = choose_parameter(get_parameter(np, internal_np::use_qem_based_energy), false); + if (use_qem_based_energy) use_postprocessing_qem = false; + +#ifndef CGAL_ACVD_DOES_NOT_USE_INTERPOLATED_CORRECTED_CURVATURES + Vertex_principal_curvatures_and_directions_map vpcd_map; + + // if adaptive clustering + if (gradation_factor > 0) + { + if constexpr (is_default_parameter::value) + { + vpcd_map = get(CGAL::dynamic_vertex_property_t>(), pmesh); + interpolated_corrected_curvatures(pmesh, parameters::vertex_principal_curvatures_and_directions_map(vpcd_map)); + } + else + vpcd_map = get_parameter(np, internal_np::vertex_principal_curvatures_and_directions_map); + } +#endif + + CGAL_precondition(CGAL::is_triangle_mesh(pmesh)); + + const double vertex_count_ratio = choose_parameter(get_parameter(np, internal_np::vertex_count_ratio), 0.1); + // TODO compute less qem things if not used + // TODO use symmetric matrices? + + //TriangleMesh pmesh = pmesh_org; + std::size_t nb_vertices = vertices(pmesh).size(); + + + // For remeshing, we might need to subdivide the mesh before clustering + // This should always hold: nb_clusters <= nb_vertices * vertex_count_ratio + // So do the following until the condition is met + if (gradation_factor == 0 && nb_clusters > nb_vertices * vertex_count_ratio) + { + std::vector lengths; + lengths.reserve(num_edges(pmesh)); + typename GT::FT cum=0; + int nbe=0; + for (Edge_descriptor e : edges(pmesh)) + { + lengths.push_back(std::sqrt(squared_distance(get(vpm, source(e, pmesh)), get(vpm, target(e, pmesh))))); + cum += lengths.back(); + ++nbe; + } + typename GT::FT threshold = 3. * cum / nbe; + + std::vector> to_split; + int ei=-1; + for (Edge_descriptor e : edges(pmesh)) + { + typename GT::FT l = lengths[++ei]; + if (l >= threshold) + { + double nb_subsegments = std::ceil(l / threshold); + to_split.emplace_back(e, nb_subsegments); + } + } + + while(!to_split.empty()) + { + std::vector> to_split_new; + for (auto [e, nb_subsegments] : to_split) + { + Halfedge_descriptor h = halfedge(e, pmesh); + typename GT::Point_3 s = get(vpm, source(h,pmesh)); + typename GT::Point_3 t = get(vpm, target(h,pmesh)); + + for (double k=1; k= threshold) + { + double nb_subsegments = std::ceil(l / threshold); + to_split_new.emplace_back(edge(hf,pmesh), nb_subsegments); + } + } + } + } + } + to_split.swap(to_split_new); + } + nb_vertices = vertices(pmesh).size(); + } + + while (nb_clusters > nb_vertices * vertex_count_ratio) + { + typename GT::FT curr_factor = nb_clusters / (nb_vertices * vertex_count_ratio); + int subdivide_steps = (std::max)((int)ceil(log(curr_factor) / log(4)), 0); + + if (subdivide_steps > 0) + { + if (gradation_factor == 0) // no adaptive clustering + Subdivision_method_3::Upsample_subdivision( + pmesh, + CGAL::parameters::number_of_iterations(subdivide_steps) + ); +#ifndef CGAL_ACVD_DOES_NOT_USE_INTERPOLATED_CORRECTED_CURVATURES + else // adaptive clustering + upsample_subdivision_property( + pmesh, vpcd_map, + CGAL::parameters::number_of_iterations(subdivide_steps).vertex_principal_curvatures_and_directions_map(vpcd_map) + ); +#endif + vpm = get_property_map(CGAL::vertex_point, pmesh); + nb_vertices = num_vertices(pmesh); + } + } + + // creating needed property maps + VertexClusterMap vertex_cluster_pmap = get(CGAL::dynamic_vertex_property_t(), pmesh, -1); + VertexWeightMap vertex_weight_pmap = get(CGAL::dynamic_vertex_property_t(), pmesh, typename GT::FT(0)); + Matrix4x4 zero_mat; zero_mat.setZero(); + VertexQuadricMap vertex_quadric_pmap = get(CGAL::dynamic_vertex_property_t(), pmesh, zero_mat); + + std::vector> clusters(nb_clusters); + std::queue clusters_edges_active; + std::queue clusters_edges_new; + + // compute vertex weights (dual area), and quadrics + typename GT::FT weight_avg = 0; + for (Face_descriptor fd : faces(pmesh)) + { + typename GT::FT weight = abs(CGAL::Polygon_mesh_processing::face_area(fd, pmesh)) / 3; + + // get points of the face + Halfedge_descriptor hd = halfedge(fd, pmesh); + typename GT::Point_3 pi = get(vpm, source(hd, pmesh)); + typename GT::Vector_3 vp1(pi.x(), pi.y(), pi.z()); + hd = next(hd, pmesh); + typename GT::Point_3 pj = get(vpm, source(hd, pmesh)); + typename GT::Vector_3 vp2(pj.x(), pj.y(), pj.z()); + hd = next(hd, pmesh); + typename GT::Point_3 pk = get(vpm, source(hd, pmesh)); + typename GT::Vector_3 vp3(pk.x(), pk.y(), pk.z()); + + // compute quadric for the face + Matrix4x4 face_quadric; + if (use_qem_based_energy) + compute_qem_face(vp1, vp2, vp3, face_quadric); + + for (Vertex_descriptor vd : vertices_around_face(halfedge(fd, pmesh), pmesh)) + { + typename GT::FT vertex_weight = get(vertex_weight_pmap, vd); + + if (gradation_factor == 0) // no adaptive clustering + vertex_weight += weight; +#ifndef CGAL_ACVD_DOES_NOT_USE_INTERPOLATED_CORRECTED_CURVATURES + else // adaptive clustering + { + typename GT::FT k1 = get(vpcd_map, vd).min_curvature; + typename GT::FT k2 = get(vpcd_map, vd).max_curvature; + typename GT::FT k_sq = (k1 * k1 + k2 * k2); + vertex_weight += weight * pow(k_sq, gradation_factor / 2.0); // /2.0 because k_sq is squared + } +#endif + weight_avg += vertex_weight; + put(vertex_weight_pmap, vd, vertex_weight); + + if (use_qem_based_energy) + { + Matrix4x4 vertex_quadric = get(vertex_quadric_pmap, vd); + vertex_quadric += face_quadric; + put(vertex_quadric_pmap, vd, vertex_quadric); + } + } + } + weight_avg /= nb_vertices; + + // clamp the weights up and below by a ratio (like 10,000) * avg_weights + for (Vertex_descriptor vd : vertices(pmesh)) + { + typename GT::FT vertex_weight = get(vertex_weight_pmap, vd); + if (vertex_weight > CGAL_WEIGHT_CLAMP_RATIO_THRESHOLD * weight_avg) + put(vertex_weight_pmap, vd, CGAL_WEIGHT_CLAMP_RATIO_THRESHOLD * weight_avg); + else if (vertex_weight < 1.0 / CGAL_WEIGHT_CLAMP_RATIO_THRESHOLD * weight_avg) + put(vertex_weight_pmap, vd, 1.0 / CGAL_WEIGHT_CLAMP_RATIO_THRESHOLD * weight_avg); + } + + Random rnd = is_default_parameter::value + ? get_default_random() + : CGAL::Random(choose_parameter(get_parameter(np, internal_np::random_seed),0)); + + // randomly initialize clusters + //TODO: std::lower_bound with vertex_weight_pmap for better sampling + for (int ci = 0; ci < nb_clusters; ++ci) + { + int vi; + Vertex_descriptor vd; + do { + vi = rnd.get_int(0, num_vertices(pmesh)); + vd = *(vertices(pmesh).begin() + vi); // TODO: bad with Polyhedron + } while (get(vertex_cluster_pmap, vd) != -1); + + put(vertex_cluster_pmap, vd, ci); + typename GT::Point_3 vp = get(vpm, vd); + typename GT::Vector_3 vpv(vp.x(), vp.y(), vp.z()); + clusters[ci].add_vertex(vpv, get(vertex_weight_pmap, vd), get(vertex_quadric_pmap, vd)); + + for (Halfedge_descriptor hd : halfedges_around_source(vd, pmesh)) + clusters_edges_active.push(hd); + } + + // the energy minimization loop (clustering loop) + int nb_modifications = 0; + int nb_disconnected = 0; + + // Turned on once nb_modifications < nb_vertices * CGAL_TO_QEM_MODIFICATION_THRESHOLD + bool qem_energy_minimization = false; + int nb_loops = 0; + + QEMClusterData cluster1_v1_to_c2, cluster2_v1_to_c2, cluster1_v2_to_c1, cluster2_v2_to_c1; + std::vector frozen_clusters(nb_clusters, false); + do + { + do + { + // reset cluster data + for (int ci=0; ci(); + for ( Vertex_descriptor v : vertices( pmesh ) ) { + typename GT::FT v_weight = get(vertex_weight_pmap, v); + Matrix4x4 v_qem = get (vertex_quadric_pmap, v); + int cluster_id = get(vertex_cluster_pmap, v ); + if ( cluster_id != -1 && !frozen_clusters[cluster_id]) + clusters[ cluster_id ].add_vertex( get(vpm, v), v_weight, v_qem); + } + + CGAL_assertion_code(for (int ci=0; ci edge_seen(num_edges(pmesh), -1); + do + { + ++nb_iterations; + nb_modifications = 0; + + while (!clusters_edges_active.empty()) { + Halfedge_descriptor hi = clusters_edges_active.front(); + clusters_edges_active.pop(); + + if (edge_seen[edge(hi,pmesh)]==nb_iterations) continue; + edge_seen[edge(hi,pmesh)]=nb_iterations; + + Vertex_descriptor v1 = source(hi, pmesh); + Vertex_descriptor v2 = target(hi, pmesh); + + // add all halfedges around v1 except hi to the queue + auto push_vertex_edge_ring_to_queue = [&](Vertex_descriptor v) + { + for (Halfedge_descriptor hd : halfedges_around_source(v, pmesh)) + //TODO: if (hd != hi && hd != opposite(hi, pmesh)) + clusters_edges_new.push(hd); + }; + + + + + int c1 = get(vertex_cluster_pmap, v1); + int c2 = get(vertex_cluster_pmap, v2); + + if ( (c1!=-1 && frozen_clusters[c1]) || (c2!=-1 && frozen_clusters[c2]) ) + { + clusters_edges_new.push(hi); + continue; + } + + if (c1==c2) continue; + + if (c1 == -1) + { + // expand cluster c2 (add v1 to c2) + put(vertex_cluster_pmap, v1, c2); + typename GT::Point_3 vp1 = get(vpm, v1); + typename GT::Vector_3 vpv(vp1.x(), vp1.y(), vp1.z()); + clusters[c2].add_vertex(vpv, get(vertex_weight_pmap, v1), get(vertex_quadric_pmap, v1)); + clusters[c2].last_modification_iteration = nb_iterations; + push_vertex_edge_ring_to_queue(v1); + ++nb_modifications; + } + else if (c2 == -1) + { + // expand cluster c1 (add v2 to c1) + put(vertex_cluster_pmap, v2, c1); + typename GT::Point_3 vp2 = get(vpm, v2); + typename GT::Vector_3 vpv(vp2.x(), vp2.y(), vp2.z()); + clusters[c1].add_vertex(vpv, get(vertex_weight_pmap, v2), get(vertex_quadric_pmap, v2)); + clusters[c1].last_modification_iteration = nb_iterations; + push_vertex_edge_ring_to_queue(v2); + ++nb_modifications; + } + else + { + if ( ( clusters[ c1 ].last_modification_iteration < nb_iterations - 1 ) && + ( clusters[ c2 ].last_modification_iteration < nb_iterations - 1 ) ) + { + clusters_edges_new.push(hi); + continue; + } + + // topological test to avoid creating disconnected clusters + auto is_topologically_valid_merge = [&](Halfedge_descriptor hv, int cluster_id) + { + CGAL_assertion(get(vertex_cluster_pmap,target(hv, pmesh))==cluster_id); + Halfedge_descriptor h=hv; + bool in_cluster=false; + int nb_cc_cluster=0; + do{ + h=next(h, pmesh); + + int ci = get(vertex_cluster_pmap, target(h,pmesh)); + if (in_cluster) + { + if (ci!=cluster_id) in_cluster=false; + } + else + { + if (ci==cluster_id) + { + in_cluster=true; + if (++nb_cc_cluster>1) return false; + } + } + h=opposite(h, pmesh); + } + while(h!=hv); + + return true; + }; + + // compare the energy of the 3 cases + typename GT::Point_3 vp1 = get(vpm, v1); + typename GT::Vector_3 vpv1(vp1.x(), vp1.y(), vp1.z()); + typename GT::Point_3 vp2 = get(vpm, v2); + typename GT::Vector_3 vpv2(vp2.x(), vp2.y(), vp2.z()); + typename GT::FT v1_weight = get(vertex_weight_pmap, v1); + typename GT::FT v2_weight = get(vertex_weight_pmap, v2); + Matrix4x4 v1_qem = get (vertex_quadric_pmap, v1); + Matrix4x4 v2_qem = get (vertex_quadric_pmap, v2); + + cluster1_v2_to_c1 = clusters[c1]; + cluster2_v2_to_c1 = clusters[c2]; + cluster1_v1_to_c2 = clusters[c1]; + cluster2_v1_to_c2 = clusters[c2]; + cluster1_v2_to_c1.last_modification_iteration = nb_iterations; + cluster2_v2_to_c1.last_modification_iteration = nb_iterations; + cluster1_v1_to_c2.last_modification_iteration = nb_iterations; + cluster2_v1_to_c2.last_modification_iteration = nb_iterations; + + typename GT::FT e_no_change = clusters[c1].compute_energy(qem_energy_minimization) + + clusters[c2].compute_energy(qem_energy_minimization); + typename GT::FT e_v1_to_c2 = (std::numeric_limits< double >::max)(); + typename GT::FT e_v2_to_c1 = (std::numeric_limits< double >::max)(); + + if ( ( clusters[ c1 ].nb_vertices > 1 ) && ( !qem_energy_minimization || nb_loops<2 || is_topologically_valid_merge(opposite(hi, pmesh), c1) ) ) + { + cluster1_v1_to_c2.remove_vertex(vpv1, v1_weight, v1_qem); + cluster2_v1_to_c2.add_vertex(vpv1, v1_weight, v1_qem); + e_v1_to_c2 = cluster1_v1_to_c2.compute_energy(qem_energy_minimization) + + cluster2_v1_to_c2.compute_energy(qem_energy_minimization); + } + + if ( ( clusters[ c2 ].nb_vertices > 1 ) && ( !qem_energy_minimization || nb_loops<2 || is_topologically_valid_merge(hi, c2) ) ) + { + cluster1_v2_to_c1.add_vertex(vpv2, v2_weight, v2_qem); + cluster2_v2_to_c1.remove_vertex(vpv2, v2_weight, v2_qem); + e_v2_to_c1 = cluster1_v2_to_c1.compute_energy(qem_energy_minimization) + + cluster2_v2_to_c1.compute_energy(qem_energy_minimization); + } + + + if (e_v2_to_c1 < e_no_change && e_v2_to_c1 < e_v1_to_c2) + { + // move v2 to c1 + put(vertex_cluster_pmap, v2, c1); + clusters[c1] = cluster1_v2_to_c1; + clusters[c2] = cluster2_v2_to_c1; + push_vertex_edge_ring_to_queue(v2); + ++nb_modifications; + } + else if (e_v1_to_c2 < e_no_change) // > 2 as 1 vertex was added to c1 + { + // move v1 to c2 + put(vertex_cluster_pmap, v1, c2); + clusters[c1] = cluster1_v1_to_c2; + clusters[c2] = cluster2_v1_to_c2; + push_vertex_edge_ring_to_queue(v1); + ++nb_modifications; + } + else + { + // no change + clusters_edges_new.push(hi); + } + } + } +#ifdef CGAL_DEBUG_ACVD + std::cout << "# Modifications: " << nb_modifications << "\n"; +#endif + + clusters_edges_active.swap(clusters_edges_new); + if (use_qem_based_energy && nb_modifications < nb_vertices * CGAL_TO_QEM_MODIFICATION_THRESHOLD && nb_loops<2) + { + qem_energy_minimization = true; + break; + } + + } while (nb_modifications > 0); + + // Disconnected clusters handling + // the goal is to delete clusters with multiple connected components and only keep the largest connected component of each cluster + // For each cluster, do a BFS from a vertex in the cluster + + std::vector>> cluster_components(nb_clusters); + std::queue vertex_queue; + + // loop over vertices to compute cluster components + VertexVisitedMap vertex_visited_pmap = get(CGAL::dynamic_vertex_property_t(), pmesh, false); + for (Vertex_descriptor vd : vertices(pmesh)) + { + if (get(vertex_visited_pmap, vd)) continue; + int c = get(vertex_cluster_pmap, vd); + if (c != -1) + { + cluster_components[c].emplace_back(); + + vertex_queue.push(vd); + put(vertex_visited_pmap, vd, true); + while (!vertex_queue.empty()) + { + Vertex_descriptor v = vertex_queue.front(); + vertex_queue.pop(); + cluster_components[c].back().push_back(v); + + for (Halfedge_descriptor hd : halfedges_around_source(v, pmesh)) + { + Vertex_descriptor v2 = target(hd, pmesh); + int c2 = get(vertex_cluster_pmap, v2); + if (c2 == c && !get(vertex_visited_pmap, v2)) + { + vertex_queue.push(v2); + put(vertex_visited_pmap, v2, true); + } + } + } + } + } + + // loop over clusters to delete disconnected components except the largest one + for (int c = 0; c < nb_clusters; ++c) + { + if (cluster_components[c].size() <= 1) continue; // only one component, no need to do anything + ++nb_disconnected; + std::size_t max_component_size = 0; + std::size_t max_component_index = -1; + for (std::size_t component_i = 0; component_i < cluster_components[c].size(); ++component_i) + { + if (cluster_components[c][component_i].size() > max_component_size) + { + max_component_size = cluster_components[c][component_i].size(); + max_component_index = component_i; + } + } + // set cluster to -1 for all components except the largest one + for (std::size_t component_i = 0; component_i < cluster_components[c].size(); ++component_i) + { + if (component_i != max_component_index) + { + for (Vertex_descriptor vd : cluster_components[c][component_i]) + { + put(vertex_cluster_pmap, vd, -1); + // remove vd from cluster c + typename GT::Point_3 vp = get(vpm, vd); + typename GT::Vector_3 vpv(vp.x(), vp.y(), vp.z()); + clusters[c].remove_vertex(vpv, get(vertex_weight_pmap, vd), get(vertex_quadric_pmap, vd)); + // add all halfedges around v except hi to the queue + for (Halfedge_descriptor hd : halfedges_around_source(vd, pmesh)) + { + // add hd to the queue if its target is not in the same cluster + Vertex_descriptor v2 = target(hd, pmesh); + int c2 = get(vertex_cluster_pmap, v2); + if (c2 != c) + clusters_edges_active.push(hd); + } + } + } + } + } + +#ifdef CGAL_DEBUG_ACVD + std::cout << "# nb_disconnected: " << nb_disconnected << "\n"; +#endif + ++nb_loops; + + } while (nb_disconnected > 0 || nb_loops < 3 ); + + /// Construct new Mesh + std::vector valid_cluster_map(nb_clusters, -1); + std::vector points; + + std::vector> polygons; + for (int i = 0; i < nb_clusters; ++i) + { + // if (clusters[i].weight_sum > 0) + { + valid_cluster_map[i] = points.size(); + points.push_back(clusters[i].representative_point(qem_energy_minimization)); + } + } + + if (use_postprocessing_qem){ + // create a point for each cluster + std::vector> cluster_quadrics(clusters.size()); + + // initialize quadrics + for (int i = 0; i < nb_clusters; ++i) + cluster_quadrics[i].setZero(); + + // for storing the vertex_descriptor of each face + std::vector face_vertices; + + for (Face_descriptor fd : faces(pmesh)) { + // get Vs for fd + // compute qem from Vs->"vector_3"s + // add to the 3 indices of the cluster + Halfedge_descriptor hd = halfedge(fd, pmesh); + do { + Vertex_descriptor vd = target(hd, pmesh); + face_vertices.push_back(vd); + hd = next(hd, pmesh); + } while (hd != halfedge(fd, pmesh)); + + auto p_i = get(vpm, face_vertices[0]); + typename GT::Vector_3 vec_1 = typename GT::Vector_3(p_i.x(), p_i.y(), p_i.z()); + p_i = get(vpm, face_vertices[1]); + typename GT::Vector_3 vec_2 = typename GT::Vector_3(p_i.x(), p_i.y(), p_i.z()); + p_i = get(vpm, face_vertices[2]); + typename GT::Vector_3 vec_3 = typename GT::Vector_3(p_i.x(), p_i.y(), p_i.z()); + + int c_1 = get(vertex_cluster_pmap, face_vertices[0]); + int c_2 = get(vertex_cluster_pmap, face_vertices[1]); + int c_3 = get(vertex_cluster_pmap, face_vertices[2]); + + if (c_1 != -1 && c_2 != -1 && c_3 != -1) + { + Eigen::Matrix q; + compute_qem_face(vec_1, vec_2, vec_3, q); + cluster_quadrics[c_1] += q; + cluster_quadrics[c_2] += q; + cluster_quadrics[c_3] += q; + } + + face_vertices.clear(); + } + + int valid_index = 0; + + for (int i = 0; i < nb_clusters; ++i) + { + if (clusters[i].weight_sum > 0) + { + int rank_deficiency = 0; + compute_representative_point(cluster_quadrics[i], points[valid_index], rank_deficiency); + ++valid_index; + } + } + } + + // extract boundary cycles + std::vector border_hedges; + extract_boundary_cycles(pmesh, std::back_inserter(border_hedges)); + + // loop over boundary loops + for (Halfedge_descriptor hd : border_hedges) + { + Halfedge_descriptor hd1 = hd; + + int cb_first = -1; + + do + { + // 1- get the target and source vertices vt, vs + // 2- if the target and source vertices are in different clusters, create a new vertex vb between them vb = (vt + vs) / 2 + // 3- make a new face with the new vertex vb and the centers of the clusters of vt and vs + // 4- also make a new face with vb, the next vb, and the center of the cluster of vt + + Vertex_descriptor vt = target(hd1, pmesh); + Vertex_descriptor vs = source(hd1, pmesh); + + int ct = get(vertex_cluster_pmap, vt); + int cs = get(vertex_cluster_pmap, vs); + + if (ct != cs) + { + typename GT::Point_3 vt_p = get(vpm, vt); + typename GT::Point_3 vs_p = get(vpm, vs); + typename GT::Vector_3 vt_v(vt_p.x(), vt_p.y(), vt_p.z()); + typename GT::Vector_3 vs_v(vs_p.x(), vs_p.y(), vs_p.z()); + + typename GT::Vector_3 vb_v = (vt_v + vs_v) / 2; + typename GT::Point_3 vb_p(vb_v.x(), vb_v.y(), vb_v.z()); + + points.push_back(vb_p); + + std::size_t cb = points.size() - 1; + + if (cb_first == -1) + cb_first = cb; + + std::size_t ct_mapped = valid_cluster_map[ct], cs_mapped = valid_cluster_map[cs]; + + if (ct_mapped != std::size_t(-1) && cs_mapped != std::size_t(-1)) + { + std::array polygon = {ct_mapped, cb, cs_mapped}; + polygons.push_back(polygon); + + // after the loop, the last cb+1 should be modified to the first cb + polygon = {cb, ct_mapped, cb + 1}; + polygons.push_back(polygon); + } + } + hd1 = next(hd1, pmesh); + } while (hd1 != hd); + polygons[polygons.size() - 1][2] = cb_first; + } + + // create a triangle for each face with all vertices in 3 different clusters + for (Face_descriptor fd : faces(pmesh)) + { + Halfedge_descriptor hd1 = halfedge(fd, pmesh); + Vertex_descriptor v1 = source(hd1, pmesh); + Halfedge_descriptor hd2 = next(hd1, pmesh); + Vertex_descriptor v2 = source(hd2, pmesh); + Halfedge_descriptor hd3 = next(hd2, pmesh); + Vertex_descriptor v3 = source(hd3, pmesh); + + int c1 = get(vertex_cluster_pmap, v1); + int c2 = get(vertex_cluster_pmap, v2); + int c3 = get(vertex_cluster_pmap, v3); + + if (c1 != c2 && c1 != c3 && c2 != c3) + { + std::size_t c1_mapped = valid_cluster_map[c1], c2_mapped = valid_cluster_map[c2], c3_mapped = valid_cluster_map[c3]; + if (c1_mapped != std::size_t(-1) && c2_mapped != std::size_t(-1) && c3_mapped != std::size_t(-1)) + { + std::array polygon = {c1_mapped, c2_mapped, c3_mapped}; + polygons.push_back(polygon); + } + } + } + +#ifdef CGAL_DEBUG_ACVD + static int kkk=-1; + CGAL::IO::write_polygon_soup("/tmp/soup_"+std::to_string(++kkk)+".off", points, polygons); + dump_mesh_with_cluster_colors(pmesh, vertex_cluster_pmap, "/tmp/cluster_"+std::to_string(kkk)+".ply"); +#endif + + std::vector > edge_map(points.size()); + for (const std::array & p : polygons) + { + for (int i=0; i<3; ++i) + { + std::pair edge = make_sorted_pair(p[i], p[(i+1)%3]); + edge_map[edge.first].emplace(edge.second,0).first->second+=1; + } + } + + std::vector nm_clusters; + for (std::size_t i=0; i2) + { +#ifdef CGAL_DEBUG_ACVD + std::cout << "non-manifold edge : " << i << " " << j << "\n"; +#endif + nm_clusters.push_back(i); + nm_clusters.push_back(j); + } + } + + + // detect isolated graph edges + for (Edge_descriptor ed : edges(pmesh)) + { + int c1 = get(vertex_cluster_pmap, source(ed, pmesh)); + int c2 = get(vertex_cluster_pmap, target(ed, pmesh)); + if (c1==-1 || c2 ==-1) throw std::runtime_error("non assigned vertex"); + + if (c1==c2) continue; + if (c2 > > link_edges(points.size()); + for (const std::array & p : polygons) + { + for (int i=0; i<3; ++i) + //TODO: skip if in nm_clusters + link_edges[p[i]].emplace_back(p[(i+1)%3], p[(i+2)%3]); + } + + using Graph = boost::adjacency_list ; + for (std::size_t i=0; i< points.size(); ++i) + { + if (std::binary_search(nm_clusters.begin(), nm_clusters.end(), i)) continue; + std::vector descriptors(points.size(), Graph::null_vertex()); + + + Graph graph; + for (const auto& p : link_edges[i]) + { + if (descriptors[p.first] == Graph::null_vertex()) descriptors[p.first] = add_vertex(graph); + if (descriptors[p.second] == Graph::null_vertex()) descriptors[p.second] = add_vertex(graph); + add_edge(descriptors[p.first], descriptors[p.second], graph); + } + + std::map the_map; + + if (boost::connected_components(graph, boost::make_assoc_property_map(the_map)) > 1) + { +#ifdef CGAL_DEBUG_ACVD + std::cout << "non-manifold vertex " << i << "\n"; +#endif + nm_clusters.push_back(i); + } + } + + if (nm_clusters.empty()) + { + //dump_mesh_with_cluster_colors(pmesh, vertex_cluster_pmap, "/tmp/debug.ply"); + //CGAL::IO::write_polygon_soup("/tmp/soup.off", points, polygons); + + // orient_polygon_soup(points, polygons); + + return std::make_pair(points, polygons); + } + + + std::vector one_ring; + for (std::size_t nmi : nm_clusters) + { + for ( auto [n, s] : edge_map[nmi] ) + one_ring.push_back(n); + } + //~ nm_clusters.insert(nm_clusters.end(), one_ring.begin(), one_ring.end()); + //~ std::sort(nm_clusters.begin(), nm_clusters.end()); + //~ nm_clusters.erase(std::unique(nm_clusters.begin(), nm_clusters.end()), nm_clusters.end()); + + + std::set nm_clusters_picked; // TODO: use cluster data instead? + + for (Vertex_descriptor v : vertices(pmesh)) + { + int c = get(vertex_cluster_pmap, v); + if (std::binary_search(nm_clusters.begin(), nm_clusters.end(), c)) + { + if (!nm_clusters_picked.insert(c).second) continue; + + if (clusters[c].nb_vertices==1) continue; + + put(vertex_cluster_pmap, v, clusters.size()); + CGAL_assertion(get(vertex_cluster_pmap, v) == clusters.size()); + clusters.emplace_back(); + } + + if (nm_clusters_picked.size()==nm_clusters.size()) break; + } + + frozen_clusters = std::vector(nb_clusters, true); + for (std::size_t nmi : nm_clusters) + frozen_clusters[nmi]=false; + for (std::size_t nmi : one_ring) + frozen_clusters[nmi]=false; + + nb_clusters=clusters.size(); + frozen_clusters.resize(nb_clusters, false); + nb_loops=0; + qem_energy_minimization=false; + + } while(true); +} + + +} // namespace internal + +/** +* \ingroup PkgPolygonMeshProcessingRef +* +* performs Approximated Centroidal Voronoi Diagram (ACVD) remeshing on a triangle mesh. The remeshing is either uniform or adaptative. +* Note that there is no guarantee that the output mesh will have the same topology as the input. +* +* @tparam TriangleMesh a model of `FaceListGraph` and `MutableFaceGraph` +* @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters". +* +* @param tmesh triangle mesh to be remeshed +* @param nb_vertices lower bound on the number of target vertices in the output mesh. +* The requested number of vertices in the output will be respected except if the input mesh is not closed +* (extra vertices will be used on the boundary), or if the number of points is too low +* and no manifold mesh can be produced with that budget of points (extra points are added to get a manifold output). +* @param np optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below. +* `GT` stands for the type of the object provided to the named parameter `geom_traits()`. +* +* \cgalNamedParamsBegin +* +* \cgalParamNBegin{vertex_principal_curvatures_and_directions_map} +* \cgalParamDescription{a property map associating principal curvatures and directions to the vertices of `tmesh`, used for adaptive clustering.} +* \cgalParamType{a class model of `ReadWritePropertyMap` with +* `boost::graph_traits::%vertex_descriptor` +* as key type and `Principal_curvatures_and_directions` as value type.} +* \cgalParamExtra{If this parameter is omitted, but `gradation_factor` is provided, an internal property map +* will be created and curvature values will be computed using the function `interpolated_corrected_curvatures()`.} +* \cgalParamNEnd +* +* \cgalParamNBegin{gradation_factor} +* \cgalParamDescription{a factor used to gradate the weights of the vertices based on their curvature values. +* The larger the value is, the more the curvature will impact the distribution of output vertices +* on the surface. The original paper recommends the value `1.5`.} +* \cgalParamType{`GT::FT`} +* \cgalParamDefault{0} +* \cgalParamExtra{If this parameter is omitted, no adaptive clustering will be performed.} +* \cgalParamNEnd +* +* \cgalParamNBegin{use_postprocessing_qem} +* \cgalParamDescription{indicates if a projection step using quadric error metric should be applied to cluster centers at the end of the minimization, +* in order for example to recover sharp features. +* This is a fast method but can result in badly shaped triangles and even self-intersections.} +* \cgalParamType{Boolean} +* \cgalParamDefault{false} +* \cgalParamNEnd +* +* \cgalParamNBegin{use_qem_based_energy} +* \cgalParamDescription{indicates if quadric error metric should be applied during the minimization algorithm in order for example to recover sharp features. +* This is slower than using `use_postprocessing_qem(true)`, but it is more accurate.} +* \cgalParamType{Boolean} +* \cgalParamDefault{false} +* \cgalParamExtra{If this parameter is `true` then `use_postprocessing_qem` will be automatically set to `false`.} +* \cgalParamNEnd +* +* \cgalParamNBegin{vertex_count_ratio} +* \cgalParamDescription{a ratio used to control the subdivision of the input mesh in case it does not have enough vertices compared to `nb_vertices`. +* More precisely, the number of vertices of the input mesh should at least be the ratio times `nb_vertices`. +* If not, the mesh will first be subdivided until the aforementioned criterium is met.} +* \cgalParamType{`GT::FT`} +* \cgalParamDefault{0.1} +* \cgalParamExtra{A value between 0.1 and 0.01 is recommended, the smaller the better the approximation will be, but it will increase the runtime.} +* \cgalParamNEnd +* +* \cgalParamNBegin{vertex_point_map} +* \cgalParamDescription{a property map associating points to the vertices of `tmesh`.} +* \cgalParamType{a class model of `ReadWritePropertyMap` with +* `boost::graph_traits::%vertex_descriptor` +* as key type and `GT::Point_3` as value type.} +* \cgalParamDefault{`boost::get(CGAL::vertex_point, tmesh)`.} +* \cgalParamExtra{If this parameter is omitted, an internal property map for +* `CGAL::vertex_point_t` must be available in `TriangleMesh`.} +* \cgalParamNEnd +* +* \cgalParamNBegin{random_seed} +* \cgalParamDescription{a value to seed the random number generator} +* \cgalParamType{unsigned int} +* \cgalParamDefault{a value generated with `std::time()`} +* \cgalParamNEnd +* +* \cgalParamNBegin{geom_traits} +* \cgalParamDescription{an instance of a geometric traits class.} +* \cgalParamType{a class model of `Kernel`, with `GT::FT` being either `float` or `double`.} +* \cgalParamDefault{a \cgal Kernel deduced from the point type, using `CGAL::Kernel_traits`.} +* \cgalParamExtra{The geometric traits class must be compatible with the vertex point type.} +* \cgalParamNEnd +* +* \cgalNamedParamsEnd +* +* @return `true` if `nb_vertices` was sufficiently large for remeshing the input, and `false` if more points were used +* +*/ +template +bool approximated_centroidal_Voronoi_diagram_remeshing( + TriangleMesh& tmesh, + std::size_t nb_vertices, + const NamedParameters& np = parameters::default_values() + ) +{ + auto ps = internal::acvd_impl(tmesh, nb_vertices, np); + CGAL_assertion(is_polygon_soup_a_polygon_mesh(ps.second)); + + auto vpm = parameters::choose_parameter( + parameters::get_parameter(np, CGAL::vertex_point), + get_property_map(CGAL::vertex_point, tmesh)); + + remove_all_elements(tmesh); + polygon_soup_to_polygon_mesh(ps.first, ps.second, tmesh, parameters::vertex_point_map(vpm)); + return ps.first.size()==nb_vertices; +} + +} // namespace Polygon_mesh_processing + +} // namespace CGAL + +#undef CGAL_WEIGHT_CLAMP_RATIO_THRESHOLD +#undef CGAL_TO_QEM_MODIFICATION_THRESHOLD + +#endif // CGAL_PMP_ACVD_REMESHING_H diff --git a/Polygon_mesh_processing/package_info/Polygon_mesh_processing/dependencies b/Polygon_mesh_processing/package_info/Polygon_mesh_processing/dependencies index f44718e3c9bc..31ee94b381a7 100644 --- a/Polygon_mesh_processing/package_info/Polygon_mesh_processing/dependencies +++ b/Polygon_mesh_processing/package_info/Polygon_mesh_processing/dependencies @@ -31,6 +31,7 @@ Solver_interface Spatial_searching Spatial_sorting Stream_support +Subdivision_method_3 TDS_2 TDS_3 Triangulation_2 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 8563e9f9a4ee..5feb680e3a22 100644 --- a/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h +++ b/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h @@ -101,6 +101,10 @@ CGAL_add_named_parameter(vertex_mean_curvature_t, vertex_mean_curvature, vertex_ CGAL_add_named_parameter(vertex_Gaussian_curvature_t, vertex_Gaussian_curvature, vertex_Gaussian_curvature) CGAL_add_named_parameter(vertex_principal_curvatures_and_directions_t, vertex_principal_curvatures_and_directions, vertex_principal_curvatures_and_directions) CGAL_add_named_parameter(ball_radius_t, ball_radius, ball_radius) +CGAL_add_named_parameter(gradation_factor_t, gradation_factor, gradation_factor) +CGAL_add_named_parameter(use_postprocessing_qem_t, use_postprocessing_qem, use_postprocessing_qem) +CGAL_add_named_parameter(use_qem_based_energy_t, use_qem_based_energy, use_qem_based_energy) +CGAL_add_named_parameter(vertex_count_ratio_t, vertex_count_ratio, vertex_count_ratio) CGAL_add_named_parameter(outward_orientation_t, outward_orientation, outward_orientation) CGAL_add_named_parameter(overlap_test_t, overlap_test, do_overlap_test_of_bounded_sides) CGAL_add_named_parameter(preserve_genus_t, preserve_genus, preserve_genus) @@ -165,6 +169,7 @@ CGAL_add_named_parameter(postprocess_regions_t, postprocess_regions, postprocess CGAL_add_named_parameter(sizing_function_t, sizing_function, sizing_function) CGAL_add_named_parameter(bbox_scaling_t, bbox_scaling, bbox_scaling) + // 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) CGAL_add_named_parameter(get_placement_policy_t, get_placement_policy, get_placement) diff --git a/Subdivision_method_3/include/CGAL/Subdivision_method_3/subdivision_masks_3.h b/Subdivision_method_3/include/CGAL/Subdivision_method_3/subdivision_masks_3.h index af9edb119681..502677fc0a45 100644 --- a/Subdivision_method_3/include/CGAL/Subdivision_method_3/subdivision_masks_3.h +++ b/Subdivision_method_3/include/CGAL/Subdivision_method_3/subdivision_masks_3.h @@ -373,6 +373,99 @@ class Loop_mask_3 : public PQQ_stencil_3 { /// @} }; +// ====================================================================== +/*! +\ingroup PkgSurfaceSubdivisionMethod3Ref + +TODO: need to find the correct name for this scheme, if it exists. And also complete the documentation. + +The geometry mask of Upsample subdivision. This does not change the geometry of the mesh. +In other words, the original vertices are not moved, and the new vertices are +the average of the incident vertices. This is useful for increasing the resolution +of a mesh without changing its geometry. + +A stencil determines a source neighborhood +whose points contribute to the position of a refined point. +The geometry mask of a stencil specifies +the computation on the nodes of the stencil. +`Upsample_mask_3` implements the geometry masks of +Upsample subdivision on a triangulated model of `MutableFaceGraph`, +such as `Polyhedron_3` and `Surface_mesh`. + +\tparam PolygonMesh must be a model of the concept `MutableFaceGraph`. Additionally all faces must be triangles. +\tparam VertexPointMap must be a model of `WritablePropertyMap` with value type `Point_3` + +\cgalModels{PTQMask_3} + +\sa `CGAL::Subdivision_method_3` + +*/ +template ::type > +class Upsample_mask_3 : public PQQ_stencil_3 { + typedef PQQ_stencil_3 Base; +public: + typedef PolygonMesh Mesh; + +#ifndef DOXYGEN_RUNNING + typedef typename Base::vertex_descriptor vertex_descriptor; + typedef typename Base::halfedge_descriptor halfedge_descriptor; + + typedef typename Base::Kernel Kernel; + + typedef typename Base::FT FT; + typedef typename Base::Point Point; + typedef typename boost::property_traits::reference Point_ref; +#endif + + typedef Halfedge_around_target_circulator Halfedge_around_vertex_circulator; + +public: +/// \name Creation +/// @{ + + /// Constructor. + /// The default vertex point property map, `get(vertex_point, pmesh)`, is used. + Upsample_mask_3(Mesh* pmesh) + : Base(pmesh, get(vertex_point, pmesh)) + { } + + /// Constructor with a custom vertex point property map. + Upsample_mask_3(Mesh* pmesh, VertexPointMap vpmap) + : Base(pmesh, vpmap) + { } + +/// @} + +/// \name Stencil functions +/// @{ + + /// computes the Upsample edge-point `pt` of the edge `edge`, simply the midpoint of the edge. + void edge_node(halfedge_descriptor edge, Point& pt) { + Point_ref p1 = get(this->vpmap,target(edge, *(this->pmesh))); + Point_ref p2 = get(this->vpmap,target(opposite(edge, *(this->pmesh)), *(this->pmesh))); + + pt = Point(0.5 * (p1[0]+p2[0]), + 0.5 * (p1[1]+p2[1]), + 0.5 * (p1[2]+p2[2]) ); + } + + /// returns Upsample vertex-point `pt` of the vertex `vertex`, simply the original vertex. + void vertex_node(vertex_descriptor vertex, Point& pt) { + pt = get(this->vpmap,vertex); + } + + /// computes the Upsample edge-point `ept` and the Upsample vertex-point `vpt` of the border edge `edge`. + void border_node(halfedge_descriptor edge, Point& ept, Point& vpt) { + edge_node(edge, ept); + + Halfedge_around_vertex_circulator vcir(edge, *(this->pmesh)); + vpt = get(this->vpmap,target(*vcir, *(this->pmesh))); + } + +/// @} +}; //========================================================================== /// The stencil of the Dual-Quadrilateral-Quadrisection diff --git a/Subdivision_method_3/include/CGAL/Subdivision_method_3/subdivision_methods_3.h b/Subdivision_method_3/include/CGAL/Subdivision_method_3/subdivision_methods_3.h index 4aebef94d494..d035bf3cf0c5 100644 --- a/Subdivision_method_3/include/CGAL/Subdivision_method_3/subdivision_methods_3.h +++ b/Subdivision_method_3/include/CGAL/Subdivision_method_3/subdivision_methods_3.h @@ -199,6 +199,56 @@ void Loop_subdivision(PolygonMesh& pmesh, const NamedParameters& np = parameters } // ----------------------------------------------------------------------------- +/*! + * + * applies Upsample subdivision several times on the control mesh `pmesh`. + * The geometry of the refined mesh is computed by the geometry policy mask `Upsample_mask_3`. + * Which is similar to Loop subdivision but does not change the shape of the mesh (only the connectivity). + * The new vertices are trivially computed as the average of the incident vertices (midpoint of edge). + * This function overwrites the control mesh `pmesh` with the subdivided mesh. + + * @tparam PolygonMesh a model of `MutableFaceGraph` + * @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters" + * + * @param pmesh a polygon mesh + * @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` + * should be available for the vertices of `pmesh`.} + * \cgalParamNEnd + * + * \cgalParamNBegin{number_of_iterations} + * \cgalParamDescription{the number of subdivision steps} + * \cgalParamType{unsigned int} + * \cgalParamDefault{`1`} + * \cgalParamNEnd + * \cgalNamedParamsEnd + * + * \pre `pmesh` must be a triangle mesh. + **/ +template +void Upsample_subdivision(PolygonMesh& pmesh, const NamedParameters& np = parameters::default_values()) { + using parameters::choose_parameter; + using parameters::get_parameter; + + typedef typename CGAL::GetVertexPointMap::type Vpm; + Vpm vpm = choose_parameter(get_parameter(np, internal_np::vertex_point), + get_property_map(CGAL::vertex_point, pmesh)); + + unsigned int step = choose_parameter(get_parameter(np, internal_np::number_of_iterations), 1); + Upsample_mask_3 mask(&pmesh, vpm); + + for(unsigned int i = 0; i < step; i++) + internal::PTQ_1step(pmesh, vpm, mask); +} +// ----------------------------------------------------------------------------- + #ifndef DOXYGEN_RUNNING // backward compatibility #ifndef CGAL_NO_DEPRECATED_CODE