Merge branch 'PMP-smoothing-kkatrio-old' into PMP-smoothing-kkatrio

This commit is contained in:
Mael Rouxel-Labbé 2019-06-19 11:34:48 +02:00
commit 77711e4e71
31 changed files with 3356 additions and 266 deletions

View File

@ -82,6 +82,11 @@ CGAL_add_named_parameter(use_compact_clipper_t, use_compact_clipper, use_compact
CGAL_add_named_parameter(output_iterator_t, output_iterator, output_iterator)
CGAL_add_named_parameter(erase_all_duplicates_t, erase_all_duplicates, erase_all_duplicates)
CGAL_add_named_parameter(require_same_orientation_t, require_same_orientation, require_same_orientation)
CGAL_add_named_parameter(use_safety_constraints_t, use_safety_constraints, use_safety_constraints)
CGAL_add_named_parameter(use_angle_smoothing_t, use_angle_smoothing, use_angle_smoothing)
CGAL_add_named_parameter(use_area_smoothing_t, use_area_smoothing, use_area_smoothing)
CGAL_add_named_parameter(use_Delaunay_flips_t, use_Delaunay_flips, use_Delaunay_flips)
CGAL_add_named_parameter(do_project_t, do_project, do_project)
// 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)
@ -90,7 +95,6 @@ CGAL_add_named_parameter(get_placement_policy_t, get_placement_policy, get_place
//to be documented
CGAL_add_named_parameter(face_normal_t, face_normal, face_normal_map)
CGAL_add_named_parameter(random_seed_t, random_seed, random_seed)
CGAL_add_named_parameter(do_project_t, do_project, do_project)
CGAL_add_named_parameter(tolerance_map_t, tolerance_map, tolerance_map)
//internal

View File

@ -88,6 +88,11 @@ void test(const NamedParameters& np)
assert(get_param(np, CGAL::internal_np::erase_all_duplicates).v == 48);
assert(get_param(np, CGAL::internal_np::require_same_orientation).v == 49);
assert(get_param(np, CGAL::internal_np::use_bool_op_to_clip_surface).v == 50);
assert(get_param(np, CGAL::internal_np::use_angle_smoothing).v == 51);
assert(get_param(np, CGAL::internal_np::use_area_smoothing).v == 52);
assert(get_param(np, CGAL::internal_np::use_safety_constraints).v == 53);
assert(get_param(np, CGAL::internal_np::use_bool_op_to_clip_surface).v == 54);
assert(get_param(np, CGAL::internal_np::use_safety_constraints).v == 56);
// Named parameters that we use in the package 'Surface Mesh Simplification'
assert(get_param(np, CGAL::internal_np::get_cost_policy).v == 34);
@ -169,6 +174,11 @@ void test(const NamedParameters& np)
check_same_type<48>(get_param(np, CGAL::internal_np::erase_all_duplicates));
check_same_type<49>(get_param(np, CGAL::internal_np::require_same_orientation));
check_same_type<50>(get_param(np, CGAL::internal_np::use_bool_op_to_clip_surface));
check_same_type<51>(get_param(np, CGAL::internal_np::use_angle_smoothing));
check_same_type<52>(get_param(np, CGAL::internal_np::use_area_smoothing));
check_same_type<53>(get_param(np, CGAL::internal_np::use_Delaunay_flips));
check_same_type<54>(get_param(np, CGAL::internal_np::use_safety_constraints));
check_same_type<56>(get_param(np, CGAL::internal_np::use_safety_constraints));
// Named parameters that we use in the package 'Surface Mesh Simplification'
check_same_type<34>(get_param(np, CGAL::internal_np::get_cost_policy));
@ -243,16 +253,20 @@ int main()
.weight_calculator(A<39>(39))
.preserve_genus(A<40>(40))
.verbosity_level(A<41>(41))
.use_binary_mode(A<51>(51))
.projection_functor(A<42>(42))
.throw_on_self_intersection(A<43>(43))
.clip_volume(A<44>(44))
.use_compact_clipper(A<45>(45))
.use_bool_op_to_clip_surface(A<50>(50))
.apply_per_connected_component(A<46>(46))
.output_iterator(A<47>(47))
.erase_all_duplicates(A<48>(48))
.require_same_orientation(A<49>(49))
.use_bool_op_to_clip_surface(A<50>(50))
.use_binary_mode(A<51>(51))
.use_angle_smoothing_t(A<52>(52))
.use_area_smoothing_t(A<53>(53))
.use_Delaunay_flips_t(A<54>(54))
.use_safety_constraints(A<56>(56))
);
return EXIT_SUCCESS;

View File

@ -606,6 +606,16 @@ and, by extension, of surface meshes (see Section \ref BGLPartitioning of the pa
More information is available on the METIS library
at <A HREF="http://glaros.dtc.umn.edu/gkhome/metis/metis/overview">`http://glaros.dtc.umn.edu/gkhome/metis/metis/overview`</A>.
\subsection thirdpartyCeres Ceres Solver
\sc{Ceres} is an open source C++ library for modeling and solving large, complicated optimization problems.
In \cgal, \sc{Ceres} is used by the \ref PkgPolygonMeshProcessingRef package for mesh smoothing, which
requires solving complex non-linear least squares problems.
Visit the official website of the library at <A HREF="http://ceres-solver.org/index.html">`ceres-solver.org`</A>
for more information.
\section secbuilding Building CGAL
The results of a successful configuration are build files that control the build step.

View File

@ -3013,6 +3013,36 @@ pages = "207--221"
HAL_VERSION = {v1}
}
@article{cgal:sg-hqct-04,
title={High quality compatible triangulations},
author={V. Surazhsky and C. Gotsman},
journal={Engineering with Computers},
volume={20},
number={2},
pages={147--156},
year={2004},
publisher={springer}
}
@inproceedings{cgal:dmsb-ifamdcf-99
,author = {M. Desbrun and M. Meyer and P. Schr{\"o}der and A. H. Barr}
,title = {Implicit Fairing of Arbitrary Meshes using
Diffusion and Curvature Flow}
,booktitle = "Computer Graphics (Proc. SIGGRAPH '99)"
,year = 1999
,pages = {317--324}
}
@inproceedings{kazhdan2012can,
title={Can Mean-Curvature Flow be Modified to be Non-singular?},
author={Kazhdan, Michael and Solomon, Jake and Ben-Chen, Mirela},
booktitle={Computer Graphics Forum},
volume={31},
number={5},
pages={1745--1754},
year={2012},
organization={Wiley Online Library}
}
% ----------------------------------------------------------------------------
% END OF BIBFILE
% ----------------------------------------------------------------------------

View File

@ -6,6 +6,12 @@ Release 5.0
Release date: September 2019
### 2D and 3D Linear Geometry Kernel
- Added `ComputeApproximateAngle_3` in the 2D/3D Kernel concept to compute
the approximate dihedral angle between 2 vectors. Corresponding functors
in the model (`Compute_approximate_angle_3`) and free function (`approximate_angle`)
are also added.
### 2D Triangulations
- **Breaking change**: Removed the functions `CGAL::Constrained_triangulation_plus_2::
@ -17,12 +23,12 @@ Release date: September 2019
- **Breaking change**: The constructor and the `insert()` function of `CGAL::Triangulation_2` which takes
a range of points as argument no longer performs a `spatial_sort()` of the points.
- Add constructor and `insert()` function to `CGAL::Triangulation_2` that takes a range of points with info.
### 3D Triangulations
- **Breaking change**: The constructor and the `insert()` function of `CGAL::Triangulation_3` which takes
a range of points as argument no longer performs a `spatial_sort()` of the points.
- Add constructor and `insert()` function to `CGAL::Triangulation_3` that takes a range of points with info.
### Surface Mesh
- New functions to read and write using the PLY format,
`CGAL::read_ply()` and `CGAL::write_ply()`, allowing to save and
@ -41,10 +47,14 @@ Release date: September 2019
### Polygon Mesh Processing
- Added the function `CGAL::Polygon_mesh_processing::centroid()`, which computes
the centroid of a closed triangle mesh.
- Added the functions `CGAL::Polygon_mesh_processing::stitch_boundary_cycle()` and
- Added the functions `CGAL::Polygon_mesh_processing::stitch_boundary_cycle()` and
`CGAL::Polygon_mesh_processing::stitch_boundary_cycles()`, which can be used
to try and merge together geometrically compatible but combinatorially different halfedges
that belong to the same boundary cycle.
- Added the mesh smoothing function `smooth_mesh()`, which can be used to
improve the quality of triangle elements based on various geometric characteristics.
- Added the shape smoothing function `smooth_shape()`, which can be used to
smooth the surface of a triangle mesh, using the mean curvature flow to perform noise removal.
### IO Streams
- **Breaking change:** The API of `CGAL::Color` has been cleaned up.
@ -58,12 +68,6 @@ Release date: September 2019
### 3D Boolean Operations on Nef Polyhedra
- Added a function to convert a Nef_polyhedron_3 to a polygon soup: `CGAL::convert_nef_to_polygon_soup()`
### 2D and 3D Linear Geometry Kernel
- Add `ComputeApproximateAngle_3` in the 2D/3D Kernel concept to compute
the approximate dihedral angle between 2 vectors. Corresponding functors
in the model (`Compute_approximate_angle_3`) and free function (`approximate_angle`)
are also added.
### IO Streams
- Added new functions to support some parts of the WKT file format:
- `CGAL::read_point_WKT()`

View File

@ -136,12 +136,15 @@ and C<sup>2</sup> continuity.\n
\cgalNPEnd
\cgalNPBegin{sparse_linear_solver} \anchor PMP_sparse_linear_solver
is the solver used in `fair()`.\n
is the solver used.\n
<b>Type:</b> a class model of `SparseLinearAlgebraWithFactorTraits_d`.\n
<b>Default:</b> if \ref thirdpartyEigen "Eigen" 3.2 (or greater) is available and
`CGAL_EIGEN3_ENABLED` is defined, then the following overload of `Eigen_solver_traits`
is provided as default value:\n
in `fair()`:\n
\code CGAL::Eigen_solver_traits<Eigen::SparseLU<CGAL::Eigen_sparse_matrix<double>::EigenType, Eigen::COLAMDOrdering<int> > > \endcode
in `smooth_shape()`:\n
\code CGAL::Eigen_solver_traits<Eigen::BiCGSTAB<CGAL::Eigen_sparse_matrix<double>::EigenType, Eigen::IncompleteLUT<double> > > \endcode
\cgalNPEnd
\cgalNPBegin{number_of_iterations} \anchor PMP_number_of_iterations
@ -390,6 +393,47 @@ if orientation should matter when determining whether two faces are duplicates.
<b>Default:</b> `false`
\cgalNPEnd
\cgalNPBegin{use_angle_smoothing_t} \anchor PMP_use_angle_smoothing
Parameter used in the function `smooth_mesh()` to indicate if angle-based smoothing should be used.
When this type of smoothing is used, the algorithm attempts to equalize angles incident to each vertex.
\n
<b>Type:</b> `bool` \n
<b>Default:</b> `true`
\cgalNPEnd
\cgalNPBegin{use_area_smoothing_t} \anchor PMP_use_area_smoothing
Parameter used in the function `smooth_mesh()` to indicate if area-based smoothing should be used.
When this type of smoothing is used, the algorithm attempts to equalize the areas of the triangles
incident to each vertex. Since this can create elongated triangles, a second phase uses Delaunay-based
flips to recover good shapes, unless specified otherwise (see below).
\n
<b>Type:</b> `bool` \n
<b>Default:</b> `true`
\cgalNPEnd
\cgalNPBegin{use_Delaunay_flips_t} \anchor PMP_use_Delaunay_flips
Parameter used in the function `smooth_mesh()` to indicate if Delaunay-based flips should be used
after area-based smoothing has been performed. A user wishing to preserve combinatorial information
can set this parameter to `false`, but the mesh might have elongated elements.
\n
<b>Type:</b> `bool` \n
<b>Default:</b> `true`
\cgalNPEnd
\cgalNPBegin{use_safety_constraints} \anchor PMP_use_safety_constraints
Parameter used in the function `smooth_mesh()` to indicate if some sanity checks should be used to decide
if the move of a vertex should be applied or rejected. These sanity checks consists of checking that
no face incident to the vertex becomes inverted and that the minimum angle of the incident faces
is not decreased by the move.
\n
<b>Type:</b> `bool` \n
<b>Default:</b> `true`
\cgalNPEnd
\cgalNPTableEnd
*/

View File

@ -96,6 +96,8 @@ and provides a list of the parameters that are used in this package.
\cgalCRPSection{Meshing Functions}
- `CGAL::Polygon_mesh_processing::fair()`
- `CGAL::Polygon_mesh_processing::refine()`
- `CGAL::Polygon_mesh_processing::smooth_mesh()`
- `CGAL::Polygon_mesh_processing::smooth_shape()`
- `CGAL::Polygon_mesh_processing::triangulate_face()`
- `CGAL::Polygon_mesh_processing::triangulate_faces()`
- \link PMP_meshing_grp `CGAL::Polygon_mesh_processing::isotropic_remeshing()` \endlink

View File

@ -39,7 +39,8 @@ and provides a list of the parameters that are used in this package.
\subsection PMPOutline Outline
The algorithms described in this manual are organized in sections:
- \ref PMPMeshing : meshing algorithms, including triangulation of non-triangulated
meshes, refinement, optimization by fairing, and isotropic remeshing of triangulated surface meshes.
meshes, refinement, optimization by fairing, isotropic remeshing of triangulated surface meshes
and smoothing algorithms.
- \ref Coref_section : methods to corefine triangle meshes and to compute
boolean operations out of corefined closed triangle meshes.
- \ref PMPHoleFilling : available hole filling algorithms, which can possibly be combined with refinement and fairing.
@ -122,6 +123,58 @@ Isotropic remeshing. (a) Triangulated input surface mesh.
(d) Surface mesh with the selection uniformly remeshed.
\cgalFigureEnd
\subsubsection Smoothing
Smoothing of a triangulated mesh region can be achieved with algorithms that aim at either mesh smoothing or shape smoothing.
While mesh smoothing is achieved by improving the quality of triangles based on criteria such as angle and area,
shape smoothing is designed to be \e intrinsic, depending as little as possible on the discretization
and smoothing the shape alone without optimizing the shape of the triangles.
- Mesh smoothing: `CGAL::Polygon_mesh_processing::smooth_mesh()` moves vertices to optimize geometry around each vertex:
it can try to equalize the angles between incident edges, or (and) move vertices so that areas of adjacent triangles tend to equalize.
Border vertices are considered constrained and do not move at any step of the procedure. No vertices are inserted at any time.
Angle and area smoothing algorithms are based on Surazhsky and Gotsman \cgalCite{cgal:sg-hqct-04}.
Since area smoothing considers only areas as a smoothing criterion, it may result in long and skinny
triangles. However, area smoothing is guaranteed to improve the spatial distribution of the vertices
over the area that is being smoothed. A simple example can be found in \ref Polygon_mesh_processing/mesh_smoothing_example.cpp.
\cgalFigureAnchor{PMPFigMeshSmoothing}
<center>
<img src="mesh_smoothing.png" style="max-width:70%;"/>
</center>
\cgalFigureCaptionBegin{PMPFigMeshSmoothing}
Mesh smoothing of the closed surface <i>blobby</i>, containing self-intersections (circled in red). For each smoothing combination,
10 iterations were applied. From left to right:
(a) Input mesh;
(b) Smoothing based on areas without using Delaunay flips;
(c) Smoothing based on areas with Delaunay flips;
(d) Smoothing based on angles;
(e) Smoothing based on angles and areas, with Delaunay flips.
\cgalFigureCaptionEnd
\cgalFigureBegin{Fig_smooth_stats, smooth_statistics.png}
Statistics for the various combinations of mesh smoothing.
\cgalFigureEnd
- Shape smoothing: `CGAL::Polygon_mesh_processing::smooth_shape()`
incrementally moves vertices towards a weighted barycenter of their neighbors along the mean curvature flow.
The curvature flow algorithm for shape smoothing is based on Desbrun et al. \cgalCite{cgal:dmsb-ifamdcf-99} and on Kazhdan et al. \cgalCite{kazhdan2012can}.
The algorithm uses the mean curvature flow to calculate the translation of vertices along the surface normal with a speed equal
to the mean curvature of the area that is being smoothed. This means that vertices on sharp corners slide faster.
If the region around a vertex is flat, this vertex does not move (zero curvature). To avoid the formation
of undesirable neck pinches (cylindrical surface areas that form singularities) the algorithm slows
down the evolution in cylindrical regions. The smoothed shape converges to a sphere while staying
conformally equivalent to its original shape.
A simple example can be found in \ref Polygon_mesh_processing/shape_smoothing_example.cpp.
\cgalFigureAnchor{PMPFigShapeSmoothing}
<center>
<img src="shape_smoothing.png" style="max-width:70%;"/>
</center>
\cgalFigureCaptionBegin{PMPFigShapeSmoothing}
Shape smoothing of the devil model, using the mean curvature flow with a time step equal to 0.05 and
constraining border vertices (located at the neck, where the mesh is open).
\cgalFigureCaptionEnd
\subsection MeshingExamples Meshing Examples

View File

@ -23,4 +23,6 @@
\example Polygon_mesh_processing/detect_features_example.cpp
\example Polygon_mesh_processing/manifoldness_repair_example.cpp
\example Polygon_mesh_processing/repair_polygon_soup_example.cpp
\example Polygon_mesh_processing/mesh_smoothing_example.cpp
\example Polygon_mesh_processing/shape_smoothing_example.cpp
*/

Binary file not shown.

After

Width:  |  Height:  |  Size: 813 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 1.5 MiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 148 KiB

View File

@ -68,6 +68,9 @@ if (EIGEN3_FOUND)
create_single_source_cgal_program( "hole_filling_example.cpp" )
create_single_source_cgal_program( "hole_filling_example_SM.cpp" )
create_single_source_cgal_program( "refine_fair_example.cpp")
create_single_source_cgal_program( "shape_smoothing_example.cpp")
create_single_source_cgal_program( "mesh_smoothing_example.cpp")
endif(EIGEN3_FOUND)
create_single_source_cgal_program( "self_intersections_example.cpp" )
@ -99,6 +102,33 @@ create_single_source_cgal_program( "detect_features_example.cpp" )
create_single_source_cgal_program( "manifoldness_repair_example.cpp" )
create_single_source_cgal_program( "repair_polygon_soup_example.cpp" )
set(SuiteSparse_USE_LAPACK_BLAS ON)
find_package(SuiteSparse QUIET NO_MODULE) # 1st: Try to locate the *config.cmake file.
if(NOT SuiteSparse_FOUND)
set(SuiteSparse_VERBOSE ON)
find_package(SuiteSparse QUIET) # 2nd: Use FindSuiteSparse.cmake module
if(SuiteSparse_FOUND)
include_directories(${SuiteSparse_INCLUDE_DIRS})
endif(SuiteSparse_FOUND)
else(NOT SuiteSparse_FOUND)
include(${USE_SuiteSparse})
endif(NOT SuiteSparse_FOUND)
if(SuiteSparse_FOUND)
if (SuiteSparse_UMFPACK_FOUND)
message(STATUS "SuiteSparse_LIBS: ${SuiteSparse_LIBRARIES}")
message(STATUS "Orbifold Tutte Embeddings will use UmfPackLU")
add_definitions(-DEIGEN_DONT_ALIGN_STATICALLY)
add_definitions(-DCGAL_SMP_USE_SPARSESUITE_SOLVERS)
else()
message(STATUS "NOTICE: The example `orbifold.cpp` will be compiled without the Sparsesuite library and UmfPack. Try setting SuiteSparse_UMF_INCLUDE_DIR and at least one of SuiteSparse_UMFPACK_LIBRARY_RELEASE and SuiteSparse_UMFPACK_LIBRARY_DEBUG to you UMFPACK installation.")
endif()
else(SuiteSparse_FOUND)
message(STATUS "NOTICE: The example `orbifold.cpp` will be compiled without the Sparsesuite library.")
endif(SuiteSparse_FOUND)
target_link_libraries( mesh_smoothing_example PRIVATE ceres glog ${SuiteSparse_LIBRARIES})
if(OpenMesh_FOUND)
create_single_source_cgal_program( "compute_normals_example_OM.cpp" )
target_link_libraries( compute_normals_example_OM PRIVATE ${OPENMESH_LIBRARIES} )

View File

@ -0,0 +1,49 @@
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/smooth_mesh.h>
#include <iostream>
#include <fstream>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Surface_mesh<K::Point_3> Mesh;
namespace PMP = CGAL::Polygon_mesh_processing;
int main(int argc, char** argv)
{
const char* filename = argc > 1 ? argv[1] : "data/mech-holes-shark.off";
std::ifstream input(filename);
Mesh mesh;
if (!input || !(input >> mesh) || mesh.is_empty())
{
std::cerr << "Not a valid .off file." << std::endl;
return EXIT_FAILURE;
}
std::set<Mesh::Vertex_index> constrained_vertices;
for(Mesh::Vertex_index v : vertices(mesh))
{
if(is_border(v, mesh))
constrained_vertices.insert(v);
}
std::cout << "Constraining: " << constrained_vertices.size() << " border vertices" << std::endl;
const unsigned int nb_iterations = 5;
CGAL::Boolean_property_map<std::set<Mesh::Vertex_index> > vcmap(constrained_vertices);
std::cout << "Smoothing... (" << nb_iterations << " iterations)" << std::endl;
PMP::smooth_mesh(mesh, PMP::parameters::number_of_iterations(nb_iterations)
.use_safety_constraints(false) // authorize all moves
.vertex_is_constrained_map(vcmap));
std::ofstream output("mesh_smoothed.off");
output << mesh;
std::cout << "Done!" << std::endl;
return EXIT_SUCCESS;
}

View File

@ -0,0 +1,35 @@
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/smooth_shape.h>
#include <iostream>
#include <fstream>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Surface_mesh<K::Point_3> Mesh;
namespace PMP = CGAL::Polygon_mesh_processing;
int main(int argc, char* argv[])
{
const char* filename = argc > 1 ? argv[1] : "data/pig.off";
std::ifstream input(filename);
Mesh mesh;
if (!input || !(input >> mesh) || mesh.is_empty())
{
std::cerr << "Not a valid .off file." << std::endl;
return EXIT_FAILURE;
}
const double time = 0.1;
const unsigned int nb_iterations = 2;
PMP::smooth_along_curvature_flow(mesh, time, PMP::parameters::number_of_iterations(nb_iterations));
std::ofstream output("mesh_shape_smoothed.off");
output << mesh;
return EXIT_SUCCESS;
}

View File

@ -271,9 +271,9 @@ struct Cotangent_value_minimum_zero_impl : CotangentValue
typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
template <class VertexPointMap>
double operator()(vertex_descriptor v0, vertex_descriptor v1, vertex_descriptor v2, const VertexPointMap& ppmap)
double operator()(vertex_descriptor v0, vertex_descriptor v1, vertex_descriptor v2, const VertexPointMap ppmap)
{
double value = CotangentValue::operator()(v0, v1, v2,ppmap);
double value = CotangentValue::operator()(v0, v1, v2, ppmap);
return (std::max)(0.0, value);
}
};
@ -284,18 +284,11 @@ template<class PolygonMesh
= Cotangent_value_Meyer<PolygonMesh, VertexPointMap> >
class Voronoi_area : CotangentValue
{
//Voronoi_area()
//{}
public:
Voronoi_area(PolygonMesh& pmesh_, VertexPointMap vpmap_)
: CotangentValue(pmesh_, vpmap_)
{}
//Voronoi_area(PolygonMesh& pmesh_)
// : CotangentValue(pmesh_, get(CGAL::vertex_point, pmesh_))
//{}
PolygonMesh& pmesh()
{
return CotangentValue::pmesh();
@ -391,8 +384,10 @@ public:
double operator()(vertex_descriptor v0, vertex_descriptor v1, vertex_descriptor v2)
{
return CotangentValue::operator()(v0, v1, v2)
/ CGAL::sqrt(squared_area(v0->point(), v1->point(), v2->point()));
return CotangentValue::operator()(v0, v1, v2) /
CGAL::sqrt(squared_area(get(this->ppmap(), v0),
get(this->ppmap(), v1),
get(this->ppmap(), v2)));
}
};
/////////////////////////////////////////////////////////////////////////////////////////
@ -683,7 +678,7 @@ public:
{
vertex_descriptor v0 = target(he, pmesh());
vertex_descriptor v1 = source(he, pmesh());
Vector vec = v0->point() - v1->point();
Vector vec = get(vpmap_, v0) - get(vpmap_, v1);
double norm = CGAL::sqrt( vec.squared_length() );
// Only one triangle for border edges
@ -714,9 +709,9 @@ private:
// Returns the tangent value of half angle v0_v1_v2/2
double half_tan_value(vertex_descriptor v0, vertex_descriptor v1, vertex_descriptor v2)
{
Vector vec0 = v1->point() - v2->point();
Vector vec1 = v2->point() - v0->point();
Vector vec2 = v0->point() - v1->point();
Vector vec0 = get(vpmap_, v1) - get(vpmap_, v2);
Vector vec1 = get(vpmap_, v2) - get(vpmap_, v0);
Vector vec2 = get(vpmap_, v0) - get(vpmap_, v1);
double e0_square = vec0.squared_length();
double e1_square = vec1.squared_length();
double e2_square = vec2.squared_length();
@ -732,8 +727,8 @@ private:
// My deviation built on Meyer_02
double half_tan_value_2(vertex_descriptor v0, vertex_descriptor v1, vertex_descriptor v2)
{
Vector a = v0->point() - v1->point();
Vector b = v2->point() - v1->point();
Vector a = get(vpmap_, v0) - get(vpmap_, v1);
Vector b = get(vpmap_, v2) - get(vpmap_, v1);
double dot_ab = a[0]*b[0] + a[1]*b[1] + a[2]*b[2];
double dot_aa = a.squared_length();
double dot_bb = b.squared_length();

View File

@ -0,0 +1,428 @@
// Copyright (c) 2018 GeometryFactory (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
// You can redistribute it and/or modify it under the terms of the GNU
// General Public License as published by the Free Software Foundation,
// either version 3 of the License, or (at your option) any later version.
//
// Licensees holding a valid commercial license may use this file in
// accordance with the commercial license agreement provided with the software.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0+
//
//
// Author(s) : Mael Rouxel-Labbé
// Konstantinos Katrioplas (konst.katrioplas@gmail.com)
#ifndef CGAL_POLYGON_MESH_PROCESSING_INTERNAL_CURVATURE_FLOW_IMPL_H
#define CGAL_POLYGON_MESH_PROCESSING_INTERNAL_CURVATURE_FLOW_IMPL_H
#include <CGAL/license/Polygon_mesh_processing/meshing_hole_filling.h>
#include <CGAL/Polygon_mesh_processing/measure.h>
#include <CGAL/Polygon_mesh_processing/Weights.h>
#include <CGAL/Polygon_mesh_processing/internal/named_function_params.h>
#include <CGAL/Polygon_mesh_processing/internal/named_params_helper.h>
#include <CGAL/Polygon_mesh_processing/internal/Smoothing/smoothing_helpers.h>
#include <CGAL/Dynamic_property_map.h>
#include <CGAL/utility.h>
#if defined(CGAL_EIGEN3_ENABLED)
#include <CGAL/Eigen_solver_traits.h>
#endif
#include <boost/graph/graph_traits.hpp>
#include <boost/property_map/property_map.hpp>
#include <algorithm>
#include <iostream>
#include <utility>
#include <vector>
namespace CGAL {
namespace Polygon_mesh_processing {
namespace internal {
// Empirically, _Meyer seems to produce the best results from the various weights available in Weights.h
template<typename TriangleMesh,
typename VertexPointMap,
typename CotangentValue = CGAL::internal::Cotangent_value_Meyer<TriangleMesh, VertexPointMap> >
struct Edge_cotangent_weight
: public CotangentValue
{
typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor halfedge_descriptor;
typedef typename boost::graph_traits<TriangleMesh>::vertex_descriptor vertex_descriptor;
Edge_cotangent_weight(TriangleMesh& pmesh_, VertexPointMap vpmap_) : CotangentValue(pmesh_, vpmap_) {}
TriangleMesh& pmesh() { return CotangentValue::pmesh(); }
double operator()(halfedge_descriptor he)
{
if(is_border_edge(he, pmesh()))
{
halfedge_descriptor h1 = next(he, pmesh());
vertex_descriptor vs = source(he, pmesh());
vertex_descriptor vt = target(he, pmesh());
vertex_descriptor v1 = target(h1, pmesh());
return CotangentValue::operator()(vs, v1, vt);
}
else
{
halfedge_descriptor h1 = next(he, pmesh());
halfedge_descriptor h2 = prev(opposite(he, pmesh()), pmesh());
vertex_descriptor vs = source(he, pmesh());
vertex_descriptor vt = target(he, pmesh());
vertex_descriptor v1 = target(h1, pmesh());
vertex_descriptor v2 = source(h2, pmesh());
return CotangentValue::operator()(vs, v1, vt) + CotangentValue::operator()(vs, v2, vt);
}
}
};
template<typename TriangleMesh,
typename VertexPointMap,
typename VertexConstraintMap,
typename SparseLinearSolver,
typename GeomTraits>
class Shape_smoother
{
typedef typename GeomTraits::FT FT;
typedef typename GeomTraits::Point_3 Point;
typedef typename GeomTraits::Vector_3 Vector;
typedef typename boost::property_traits<VertexPointMap>::reference Point_ref;
typedef CGAL::Triple<std::size_t, std::size_t, double> Triplet;
typedef typename boost::graph_traits<TriangleMesh>::vertex_descriptor vertex_descriptor;
typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor halfedge_descriptor;
typedef typename boost::graph_traits<TriangleMesh>::edge_descriptor edge_descriptor;
typedef typename boost::graph_traits<TriangleMesh>::face_descriptor face_descriptor;
typedef CGAL::dynamic_vertex_property_t<std::size_t> Vertex_local_index;
typedef typename boost::property_map<TriangleMesh, Vertex_local_index>::type IndexMap;
// linear system
typedef typename SparseLinearSolver::Matrix Eigen_matrix;
typedef typename SparseLinearSolver::Vector Eigen_vector;
public:
Shape_smoother(TriangleMesh& mesh,
VertexPointMap& vpmap,
VertexConstraintMap& vcmap,
const GeomTraits& traits)
:
mesh_(mesh),
vpmap_(vpmap),
vcmap_(vcmap),
vimap_(get(Vertex_local_index(), mesh_)),
scale_volume_after_smoothing(true),
traits_(traits),
weight_calculator_(mesh, vpmap)
{ }
template<typename FaceRange>
void init_smoothing(const FaceRange& face_range)
{
set_face_range(face_range);
std::size_t id = 0;
for(vertex_descriptor v : vertices(mesh_))
put(vimap_, v, id++);
// vertices that are not in the range or are constrained still need value '1' in D because the RHS is D * X^n
diagonal_.assign(vertices(mesh_).size(), 1.);
constrained_flags_.assign(vertices(mesh_).size(), false);
for(vertex_descriptor v : vertices(mesh_))
{
if(is_constrained(v))
{
constrained_flags_[get(vimap_, v)] = true;
// scaling things cannot preserve the position of more than a single constrained point
if(anchor_point == boost::none)
anchor_point = get(vpmap_, v);
else
scale_volume_after_smoothing = false;
}
}
if(!CGAL::is_closed(mesh_))
scale_volume_after_smoothing = false;
}
void setup_system(Eigen_matrix& A,
Eigen_vector& bx, Eigen_vector& by, Eigen_vector& bz,
std::vector<Triplet>& stiffness_elements,
const double& time)
{
compute_coefficient_matrix(A, stiffness_elements, time);
compute_rhs(bx, by, bz);
}
bool solve_system(const Eigen_matrix& A,
Eigen_vector& Xx, Eigen_vector& Xy, Eigen_vector& Xz,
const Eigen_vector& bx, const Eigen_vector& by, const Eigen_vector& bz,
SparseLinearSolver& solver)
{
FT D;
// calls compute once to factorize with the preconditioner
if(!solver.factor(A, D))
{
#ifdef CGAL_PMP_SMOOTHING_VERBOSE
std::cerr << "Could not factorize linear system with preconditioner." << std::endl;
#endif
return false;
}
if(!solver.linear_solver(bx, Xx) ||
!solver.linear_solver(by, Xy) ||
!solver.linear_solver(bz, Xz))
{
#ifdef CGAL_PMP_SMOOTHING_VERBOSE
std::cerr << "Could not solve linear system." << std::endl;
#endif
return false;
}
return true;
}
void calculate_stiffness_matrix_elements(std::vector<Triplet>& stiffness_elements)
{
CGAL_assertion(stiffness_elements.empty());
stiffness_elements.reserve(8 * vrange_.size());
std::unordered_map<std::size_t, double> diag_coeff;
for(face_descriptor f : frange_)
{
for(halfedge_descriptor hi : halfedges_around_face(halfedge(f, mesh_), mesh_))
{
// Get a single canonical non-border halfedge per edge
if(is_border(hi, mesh_))
continue;
const halfedge_descriptor hi_opp = opposite(hi, mesh_);
if(!is_border(hi_opp, mesh_) && hi < hi_opp)
continue;
const vertex_descriptor v_source = source(hi, mesh_);
const vertex_descriptor v_target = target(hi, mesh_);
const bool is_source_constrained = is_constrained(v_source);
const bool is_target_constrained = is_constrained(v_target);
if(is_source_constrained && is_target_constrained)
continue;
const FT Lij = weight_calculator_(hi);
const std::size_t i_source = get(vimap_, v_source);
const std::size_t i_target = get(vimap_, v_target);
// note that these constraints create asymmetry in the matrix
if(!is_source_constrained)
{
stiffness_elements.push_back(Triplet(i_source, i_target, Lij));
diag_coeff.insert(std::make_pair(i_source, 0)).first->second -= Lij;
}
if(!is_target_constrained)
{
stiffness_elements.push_back(Triplet(i_target, i_source, Lij));
diag_coeff.insert(std::make_pair(i_target, 0)).first->second -= Lij;
}
}
}
typename std::unordered_map<std::size_t, double>::iterator it = diag_coeff.begin(),
end = diag_coeff.end();
for(; it!=end; ++it)
stiffness_elements.push_back(Triplet(it->first, it->first, it->second));
}
void update_mesh_no_scaling(const Eigen_vector& Xx, const Eigen_vector& Xy, const Eigen_vector& Xz)
{
for(vertex_descriptor v : vrange_)
{
std::size_t index = get(vimap_, v);
const FT x_new = Xx[index];
const FT y_new = Xy[index];
const FT z_new = Xz[index];
Point new_pos(x_new, y_new, z_new);
put(vpmap_, v, new_pos);
}
}
void update_mesh(const Eigen_vector& Xx, const Eigen_vector& Xy, const Eigen_vector& Xz)
{
namespace PMP = CGAL::Polygon_mesh_processing;
if(!scale_volume_after_smoothing)
return update_mesh_no_scaling(Xx, Xy, Xz);
const FT old_vol = volume(mesh_, parameters::vertex_point_map(vpmap_).geom_traits(traits_));
// If no vertex is constrained, then the smoothed mesh will simply share the same centroid as the input mesh
Point pre_smooth_anchor_point;
if(anchor_point != boost::none)
pre_smooth_anchor_point = *anchor_point;
else
pre_smooth_anchor_point = PMP::centroid(mesh_, parameters::vertex_point_map(vpmap_).geom_traits(traits_));
for(vertex_descriptor v : vrange_)
{
std::size_t index = get(vimap_, v);
const FT x_new = Xx[index];
const FT y_new = Xy[index];
const FT z_new = Xz[index];
Point new_pos(x_new, y_new, z_new);
put(vpmap_, v, new_pos);
}
Point post_smooth_anchor_point;
if(anchor_point != boost::none)
post_smooth_anchor_point = *anchor_point;
else
post_smooth_anchor_point = PMP::centroid(mesh_, parameters::vertex_point_map(vpmap_).geom_traits(traits_));
const FT new_vol = volume(mesh_, parameters::vertex_point_map(vpmap_));
CGAL_assertion(new_vol != 0);
const FT inflating_factor = std::cbrt(CGAL::abs(old_vol / new_vol));
for(vertex_descriptor v : vertices(mesh_))
{
Vector d = traits_.construct_vector_3_object()(post_smooth_anchor_point, get(vpmap_, v));
Point new_pos = traits_.construct_translated_point_3_object()(pre_smooth_anchor_point,
traits_.construct_scaled_vector_3_object()(d, inflating_factor));
put(vpmap_, v, new_pos);
}
}
private:
bool is_constrained(const vertex_descriptor& v)
{
return get(vcmap_, v);
}
template<typename FaceRange>
void set_face_range(const FaceRange& face_range)
{
frange_.assign(face_range.begin(), face_range.end());
vrange_.reserve(3 * face_range.size());
for(face_descriptor f : face_range)
{
for(vertex_descriptor v : vertices_around_face(halfedge(f, mesh_), mesh_))
vrange_.push_back(v);
}
// get rid of duplicate vertices
std::sort(vrange_.begin(), vrange_.end());
vrange_.erase(std::unique(vrange_.begin(), vrange_.end()), vrange_.end());
}
void compute_coefficient_matrix(Eigen_matrix& A,
std::vector<Triplet>& stiffness_elements,
const double& time)
{
fill_mass_matrix();
// fill A = Mass - time * Laplacian
for(const Triplet& t : stiffness_elements)
{
std::size_t i = t.get<0>(), j = t.get<1>();
if(i != j)
A.set_coef(i, j, - time * t.get<2>(), true);
else if(!constrained_flags_[i]) // && i==j
A.set_coef(i, i, diagonal_[t.get<0>()] - time * t.get<2>(), true);
}
for(vertex_descriptor v : vrange_)
{
std::size_t index = get(vimap_, v);
if(constrained_flags_[index])
A.set_coef(index, index, 1., true);
}
// we do not call A.assemble_matrix here
// Eigen's compute during factorization does the building correctly,
// and without assemble_matrix the reference A can be used in the next iterations.
}
void fill_mass_matrix()
{
for(vertex_descriptor v : vrange_)
{
std::size_t index = get(vimap_, v);
if(!is_constrained(v))
diagonal_[index] = 0.;
}
for(face_descriptor f : frange_)
{
const double area = face_area(f, mesh_, parameters::vertex_point_map(vpmap_).geom_traits(traits_));
for(vertex_descriptor v : vertices_around_face(halfedge(f, mesh_), mesh_))
{
if(!is_constrained(v))
diagonal_[get(vimap_, v)] += area / 6.;
}
}
}
void compute_rhs(Eigen_vector& bx, Eigen_vector& by, Eigen_vector& bz)
{
for(vertex_descriptor vi : vertices(mesh_))
{
std::size_t index = get(vimap_, vi);
Point_ref p = get(vpmap_, vi);
bx.set(index, diagonal_[index] * p.x());
by.set(index, diagonal_[index] * p.y());
bz.set(index, diagonal_[index] * p.z());
}
}
private:
std::vector<vertex_descriptor> vrange_;
std::vector<face_descriptor> frange_;
TriangleMesh& mesh_;
VertexPointMap vpmap_;
VertexConstraintMap vcmap_;
IndexMap vimap_;
// Smoothing has a tendency to reduce volumes, so we can scale things back up based on the change
// of volume. We need an anchor point to scale up, either a constrained point or the centroid
// of the initial mesh if no vertex is constrained. If there is more than a constrained vertex,
// then no scaling can be done without violating the constraint.
bool scale_volume_after_smoothing;
boost::optional<Point> anchor_point;
// linear system data
std::vector<double> diagonal_; // index of vector -> index of vimap_
std::vector<bool> constrained_flags_;
const GeomTraits& traits_;
Edge_cotangent_weight<TriangleMesh, VertexPointMap> weight_calculator_;
};
} // internal
} // PMP
} // CGAL
#endif // CGAL_POLYGON_MESH_PROCESSING_INTERNAL_CURVATURE_FLOW_NEW_IMPL_H

View File

@ -0,0 +1,709 @@
// Copyright (c) 2018 GeometryFactory (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
// You can redistribute it and/or modify it under the terms of the GNU
// General Public License as published by the Free Software Foundation,
// either version 3 of the License, or (at your option) any later version.
//
// Licensees holding a valid commercial license may use this file in
// accordance with the commercial license agreement provided with the software.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0+
//
//
// Author(s) : Mael Rouxel-Labbé
// Konstantinos Katrioplas (konst.katrioplas@gmail.com)
#ifndef CGAL_POLYGON_MESH_PROCESSING_INTERNAL_MESH_SMOOTHING_IMPL_H
#define CGAL_POLYGON_MESH_PROCESSING_INTERNAL_MESH_SMOOTHING_IMPL_H
#include <CGAL/license/Polygon_mesh_processing/meshing_hole_filling.h>
#include <CGAL/Polygon_mesh_processing/compute_normal.h>
#include <CGAL/Polygon_mesh_processing/internal/Smoothing/smoothing_helpers.h>
#include <CGAL/Polygon_mesh_processing/repair.h>
#include <CGAL/AABB_tree.h>
#include <CGAL/AABB_traits.h>
#include <CGAL/AABB_triangle_primitive.h>
#include <CGAL/boost/graph/Euler_operations.h>
#include <CGAL/Dynamic_property_map.h>
#include <CGAL/Kernel/global_functions_3.h>
#include <CGAL/iterator.h>
#include <CGAL/number_type_config.h>
#include <CGAL/Origin.h>
#include <CGAL/property_map.h>
#include <CGAL/utils.h>
#include <boost/graph/graph_traits.hpp>
#include "ceres/ceres.h"
#include <algorithm>
#include <cmath>
#include <iterator>
#include <map>
#include <utility>
#include <vector>
namespace CGAL {
namespace Polygon_mesh_processing {
namespace internal {
template <typename V, typename GT>
double get_radian_angle(const V& v1, const V& v2, const GT& gt)
{
return gt.compute_approximate_angle_3_object()(v1, v2) * CGAL_PI / 180.;
}
// super naive for now. Not sure it even makes sense to do something like that for surfaces
template<typename TriangleMesh,
typename VertexPointMap,
typename ECMap,
typename GeomTraits>
class Delaunay_edge_flipper
{
typedef typename boost::graph_traits<TriangleMesh>::vertex_descriptor vertex_descriptor;
typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor halfedge_descriptor;
typedef typename boost::graph_traits<TriangleMesh>::edge_descriptor edge_descriptor;
typedef typename boost::graph_traits<TriangleMesh>::face_descriptor face_descriptor;
typedef typename boost::property_traits<VertexPointMap>::reference Point_ref;
typedef typename GeomTraits::Vector_3 Vector;
public:
Delaunay_edge_flipper(TriangleMesh& mesh,
const VertexPointMap vpmap,
const ECMap ecmap,
const GeomTraits& traits)
: mesh_(mesh), vpmap_(vpmap), ecmap_(ecmap), traits_(traits)
{ }
bool should_be_flipped(const edge_descriptor e) const
{
if(is_border(e, mesh_) || get(ecmap_, e))
return false;
const halfedge_descriptor h = halfedge(e, mesh_);
const halfedge_descriptor opp_h = opposite(h, mesh_);
vertex_descriptor v0 = source(h, mesh_);
vertex_descriptor v1 = target(h, mesh_);
vertex_descriptor v2 = target(next(h, mesh_), mesh_);
vertex_descriptor v3 = target(next(opp_h, mesh_), mesh_);
const Point_ref p0 = get(vpmap_, v0);
const Point_ref p1 = get(vpmap_, v1);
const Point_ref p2 = get(vpmap_, v2);
const Point_ref p3 = get(vpmap_, v3);
double alpha = get_radian_angle(Vector(p0 - p2), Vector(p1 - p2), traits_);
double beta = get_radian_angle(Vector(p1 - p3), Vector(p0 - p3), traits_);
// not local Delaunay if the sum of the angles is greater than pi
if(alpha + beta <= CGAL_PI)
return false;
// Don't want to flip if the other diagonal already exists
// @todo remeshing can be used to still flip those
std::pair<edge_descriptor, bool> other_hd_already_exists = edge(v2, v3, mesh_);
if(other_hd_already_exists.second)
return false;
return true;
}
template <typename Marked_edges_map, typename EdgeRange>
void add_to_stack_if_unmarked(const edge_descriptor e,
const Marked_edges_map marks,
EdgeRange& edge_range)
{
if(!get(marks, e))
{
put(marks, e, true);
edge_range.push_back(e);
}
}
public:
template <typename FaceRange>
void operator()(const FaceRange& face_range)
{
// edges to consider
std::vector<edge_descriptor> edge_range;
edge_range.reserve(3 * face_range.size());
for(face_descriptor f : face_range)
{
for(halfedge_descriptor h : halfedges_around_face(halfedge(f, mesh_), mesh_))
edge_range.push_back(edge(h, mesh_));
}
// keep unique elements
std::sort(edge_range.begin(), edge_range.end());
edge_range.erase(std::unique(edge_range.begin(), edge_range.end()), edge_range.end());
// Mark edges that are in the stack
typedef CGAL::dynamic_edge_property_t<bool> Edge_property_tag;
typedef typename boost::property_map<TriangleMesh,
Edge_property_tag>::type Marked_edges_map;
Marked_edges_map marks = get(Edge_property_tag(), mesh_);
// dynamic pmaps do not have default values...
for(edge_descriptor e : edges(mesh_))
put(marks, e, false);
for(edge_descriptor e : edge_range)
put(marks, e, true);
int flipped_n = 0;
while(!edge_range.empty())
{
edge_descriptor e = edge_range.back();
edge_range.pop_back();
put(marks, e, false);
if(should_be_flipped(e))
{
++flipped_n;
halfedge_descriptor h = halfedge(e, mesh_);
Euler::flip_edge(h, mesh_);
add_to_stack_if_unmarked(edge(next(h, mesh_), mesh_), marks, edge_range);
add_to_stack_if_unmarked(edge(prev(h, mesh_), mesh_), marks, edge_range);
add_to_stack_if_unmarked(edge(next(opposite(h, mesh_), mesh_), mesh_), marks, edge_range);
add_to_stack_if_unmarked(edge(prev(opposite(h, mesh_), mesh_), mesh_), marks, edge_range);
}
}
}
private:
TriangleMesh& mesh_;
const VertexPointMap vpmap_;
const ECMap ecmap_;
const GeomTraits& traits_;
};
template<typename TriangleMesh, typename VertexPointMap, typename GeomTraits>
class Angle_smoother
{
typedef typename boost::graph_traits<TriangleMesh>::vertex_descriptor vertex_descriptor;
typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor halfedge_descriptor;
typedef typename boost::property_traits<VertexPointMap>::reference Point_ref;
typedef typename GeomTraits::Vector_3 Vector;
typedef std::pair<halfedge_descriptor, halfedge_descriptor> He_pair;
public:
Angle_smoother(const TriangleMesh& mesh,
const VertexPointMap vpmap,
const GeomTraits& traits)
: mesh_(mesh), vpmap_(vpmap), traits_(traits)
{ }
private:
Vector rotate_edge(const halfedge_descriptor main_he,
const He_pair& incident_pair) const
{
// get common vertex around which the edge is rotated
Point_ref pt = get(vpmap_, target(main_he, mesh_));
Point_ref left_pt = get(vpmap_, source(incident_pair.first, mesh_));
Point_ref right_pt = get(vpmap_, target(incident_pair.second, mesh_));
CGAL_assertion(target(incident_pair.first, mesh_) == source(incident_pair.second, mesh_));
Vector edge1(pt, left_pt);
Vector edge2(pt, right_pt);
// find bisector
internal::normalize(edge1, traits_);
internal::normalize(edge2, traits_);
Vector bisector = traits_.construct_sum_of_vectors_3_object()(edge1, edge2);
internal::normalize(bisector, traits_);
return bisector;
}
public:
// If it's ever allowed to move vertices on the border, the min angle computations will be missing
// some values (angles incident to the border)
Vector operator()(const vertex_descriptor v) const
{
Vector move = CGAL::NULL_VECTOR;
double weights_sum = 0.;
for(halfedge_descriptor main_he : halfedges_around_source(v, mesh_))
{
He_pair incident_pair = std::make_pair(prev(opposite(main_he, mesh_), mesh_),
next(main_he, mesh_));
// avoid zero angles
Point_ref ps = get(vpmap_, source(main_he, mesh_));
Point_ref pt = get(vpmap_, target(main_he, mesh_));
Point_ref left_pt = get(vpmap_, source(incident_pair.first, mesh_));
Point_ref right_pt = get(vpmap_, target(incident_pair.second, mesh_));
CGAL_assertion(target(incident_pair.first, mesh_) == source(incident_pair.second, mesh_));
Vector left_v(pt, left_pt);
Vector right_v(pt, right_pt);
// rotate
double angle = get_radian_angle(right_v, left_v, traits_);
CGAL_warning(angle != 0.); // no degenerate faces is a precondition
if(angle == 0.)
continue;
Vector bisector = rotate_edge(main_he, incident_pair);
double scaling_factor = CGAL::approximate_sqrt(
traits_.compute_squared_distance_3_object()(get(vpmap_, source(main_he, mesh_)),
get(vpmap_, target(main_he, mesh_))));
bisector = traits_.construct_scaled_vector_3_object()(bisector, scaling_factor);
Vector ps_psi(ps, traits_.construct_translated_point_3_object()(pt, bisector));
// small angles carry more weight
double weight = 1. / (angle*angle);
weights_sum += weight;
move += weight * ps_psi;
}
if(weights_sum != 0.)
move /= weights_sum;
return move;
}
private:
const TriangleMesh& mesh_;
const VertexPointMap vpmap_;
const GeomTraits& traits_;
};
template<typename TriangleMesh, typename VertexPointMap, typename GeomTraits>
class Area_smoother
{
typedef typename boost::graph_traits<TriangleMesh>::vertex_descriptor vertex_descriptor;
typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor halfedge_descriptor;
typedef typename boost::property_traits<VertexPointMap>::value_type Point;
typedef typename boost::property_traits<VertexPointMap>::reference Point_ref;
typedef typename GeomTraits::Vector_3 Vector;
public:
Area_smoother(const TriangleMesh& mesh,
const VertexPointMap vpmap,
const GeomTraits& traits)
: mesh_(mesh), vpmap_(vpmap), traits_(traits)
{ }
private:
double element_area(const vertex_descriptor v1,
const vertex_descriptor v2,
const vertex_descriptor v3) const
{
return CGAL::to_double(CGAL::approximate_sqrt(traits_.compute_squared_area_3_object()(get(vpmap_, v1),
get(vpmap_, v2),
get(vpmap_, v3))));
}
double element_area(const Point& P,
const vertex_descriptor v2,
const vertex_descriptor v3) const
{
return CGAL::to_double(CGAL::approximate_sqrt(traits_.compute_squared_area_3_object()(P,
get(vpmap_, v2),
get(vpmap_, v3))));
}
double compute_average_area_around(const vertex_descriptor v) const
{
double sum_areas = 0.;
unsigned int number_of_edges = 0;
for(halfedge_descriptor h : halfedges_around_source(v, mesh_))
{
// opposite vertices
vertex_descriptor vi = source(next(h, mesh_), mesh_);
vertex_descriptor vj = target(next(h, mesh_), mesh_);
double S = element_area(v, vi, vj);
sum_areas += S;
++number_of_edges;
}
return sum_areas / number_of_edges;
}
struct Face_energy
{
Face_energy(const Point& pi, const Point& pj, const double s_av)
:
qx(pi.x()), qy(pi.y()), qz(pi.z()),
rx(pj.x()), ry(pj.y()), rz(pj.z()),
s_av(s_av)
{ }
// next two functions are just for convencience, the only thing ceres cares about is the operator()
template <typename T>
double area(const T x, const T y, const T z) const
{
return CGAL::approximate_sqrt(CGAL::squared_area(Point(x, y, z),
Point(qx, qy, qz),
Point(rx, ry, rz)));
}
template <typename T>
double evaluate(const T x, const T y, const T z) const { return area(x, y, z) - s_av; }
template <typename T>
bool operator()(const T* const x, const T* const y, const T* const z,
T* residual) const
{
#define CGAL_CERES_USE_NUMERIC_DIFFERENCIATION // @todo test this
#ifdef CGAL_CERES_USE_NUMERIC_DIFFERENCIATION
residual[0] = evaluate(x[0], y[0], z[0]);
#else
// Computations must be explicit so automatic differenciation can be used
T dqx = qx - x[0];
T dqy = qy - y[0];
T dqz = qz - z[0];
T drx = rx - x[0];
T dry = ry - y[0];
T drz = rz - z[0];
T vx = dqy*drz - dqz*dry;
T vy = dqz*drx - dqx*drz;
T vz = dqx*dry - dqy*drx;
T squared_area = 0.25 * (vx*vx + vy*vy + vz*vz);
T area = sqrt(squared_area);
residual[0] = area - s_av;
#endif
return true;
}
private:
const double qx, qy, qz;
const double rx, ry, rz;
const double s_av;
};
public:
Vector operator()(const vertex_descriptor v) const
{
const Point_ref vp = get(vpmap_, v);
const double S_av = compute_average_area_around(v);
const double initial_x = vp.x();
const double initial_y = vp.y();
const double initial_z = vp.z();
double x = initial_x, y = initial_y, z = initial_z;
ceres::Problem problem;
for(halfedge_descriptor h : halfedges_around_source(v, mesh_))
{
CGAL_assertion(!is_border(h, mesh_));
vertex_descriptor vi = source(next(h, mesh_), mesh_);
vertex_descriptor vj = target(next(h, mesh_), mesh_);
const Point_ref vip = get(vpmap_, vi);
const Point_ref vjp = get(vpmap_, vj);
#ifdef CGAL_CERES_USE_NUMERIC_DIFFERENCIATION
ceres::CostFunction* cost_function =
new ceres::NumericDiffCostFunction<Face_energy, ceres::CENTRAL, 1, 1, 1, 1>(new Face_energy(vip, vjp, S_av));
#else
ceres::CostFunction* cost_function =
new ceres::AutoDiffCostFunction<Face_energy, 1, 1, 1, 1>(new Face_energy(vip, vjp, S_av));
#endif
problem.AddResidualBlock(cost_function, NULL, &x, &y, &z);
}
ceres::Solver::Options options;
// options.minimizer_progress_to_stdout = true;
ceres::Solver::Summary summary;
ceres::Solve(options, &problem, &summary);
// std::cout << summary.BriefReport() << "\n";
// std::cout << "x : " << initial_x << " -> " << x << "\n";
// std::cout << "y : " << initial_y << " -> " << y << "\n";
// std::cout << "z : " << initial_z << " -> " << z << "\n";
return Vector(x - initial_x, y - initial_y, z - initial_z);
}
private:
const TriangleMesh& mesh_;
const VertexPointMap vpmap_;
const GeomTraits& traits_;
};
template<typename Optimizer, typename TriangleMesh,
typename VertexPointMap, typename VertexConstraintMap,
typename GeomTraits>
class Mesh_smoother
{
typedef typename boost::graph_traits<TriangleMesh>::vertex_descriptor vertex_descriptor;
typedef typename boost::graph_traits<TriangleMesh>::edge_descriptor edge_descriptor;
typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor halfedge_descriptor;
typedef typename boost::graph_traits<TriangleMesh>::face_descriptor face_descriptor;
typedef typename boost::property_traits<VertexPointMap>::value_type Point;
typedef typename boost::property_traits<VertexPointMap>::reference Point_ref;
typedef typename GeomTraits::FT FT;
typedef typename GeomTraits::Vector_3 Vector;
public:
Mesh_smoother(TriangleMesh& pmesh,
VertexPointMap& vpmap,
VertexConstraintMap& vcmap,
const GeomTraits& traits)
:
mesh_(pmesh), vpmap_(vpmap), vcmap_(vcmap), traits_(traits)
{}
public:
template<typename FaceRange>
void set_vertex_range(const FaceRange& face_range)
{
vrange_.reserve(3 * face_range.size());
for(face_descriptor f : face_range)
{
for(vertex_descriptor v : vertices_around_face(halfedge(f, mesh_), mesh_))
vrange_.push_back(v);
}
// get rid of duplicate vertices
std::sort(vrange_.begin(), vrange_.end());
vrange_.erase(std::unique(vrange_.begin(), vrange_.end()), vrange_.end());
}
template<typename FaceRange>
void init_smoothing(const FaceRange& face_range)
{
CGAL_precondition(CGAL::is_triangle_mesh(mesh_));
CGAL_precondition_code(std::set<typename boost::graph_traits<TriangleMesh>::face_descriptor> degen_faces;)
CGAL_precondition_code(CGAL::Polygon_mesh_processing::degenerate_faces(
mesh_, std::inserter(degen_faces, degen_faces.begin()),
CGAL::parameters::vertex_point_map(vpmap_).geom_traits(traits_));)
CGAL_precondition(degen_faces.empty());
set_vertex_range(face_range);
}
// generic optimizer, the move is computed by 'Optimizer'
std::size_t optimize(const bool use_sanity_checks = true,
const bool apply_moves_in_single_batch = false,
const bool enforce_no_min_angle_regression = false)
{
typedef CGAL::dynamic_vertex_property_t<Point> Vertex_property_tag;
typedef typename boost::property_map<TriangleMesh,
Vertex_property_tag>::type Position_map;
Position_map new_positions = get(Vertex_property_tag(), mesh_);
Optimizer compute_move(mesh_, vpmap_, traits_);
#ifdef CGAL_PMP_SMOOTHING_DEBUG
double total_displacement = 0;
#endif
std::size_t moved_points = 0;
for(vertex_descriptor v : vrange_)
{
if(is_border(v, mesh_) || is_constrained(v))
continue;
// compute normal to v
Vector vn = compute_vertex_normal(v, mesh_, CGAL::parameters::vertex_point_map(vpmap_)
.geom_traits(traits_));
// calculate movement
const Point_ref pos = get(vpmap_, v);
Vector move = compute_move(v);
// Gram Schmidt so that the new location is on the tangent plane of v (i.e. do mv -= (mv*n)*n)
const FT sp = traits_.compute_scalar_product_3_object()(vn, move);
move = traits_.construct_sum_of_vectors_3_object()(move,
traits_.construct_scaled_vector_3_object()(vn, - sp));
Point new_pos = pos + move;
if((!use_sanity_checks || !does_move_create_bad_faces(v, new_pos)) &&
(!enforce_no_min_angle_regression || does_improve_min_angle_in_star(v, new_pos)))
{
#ifdef CGAL_PMP_SMOOTHING_DEBUG
std::cout << "moving " << get(vpmap_, v) << " to " << new_pos << std::endl;
total_displacement += CGAL::approximate_sqrt(traits_.compute_squared_length_3_object()(move));
#endif
if(apply_moves_in_single_batch)
put(new_positions, v, new_pos);
else
put(vpmap_, v, new_pos);
++moved_points;
}
else // some sanity check failed
{
#ifdef CGAL_PMP_SMOOTHING_DEBUG
std::cout << "move is rejected!" << std::endl;
#endif
if(apply_moves_in_single_batch)
put(new_positions, v, pos);
}
}
// update locations
if(apply_moves_in_single_batch)
{
for(vertex_descriptor v : vrange_)
{
if(is_border(v, mesh_) || is_constrained(v))
continue;
put(vpmap_, v, get(new_positions, v));
}
}
#ifdef CGAL_PMP_SMOOTHING_DEBUG
std::cout << "moved: " << moved_points << " points based on angle." << std::endl;
std::cout << "total displacement: " << total_displacement << std::endl;
std::cout << "not improved min angle: " << vrange_.size() - moved_points << " times." << std::endl;
#endif
return moved_points;
}
template <typename AABBTree>
void project_to_surface(const AABBTree& tree)
{
#ifdef CGAL_PMP_SMOOTHING_DEBUG
std::cout << "Projecting back to the surface" << std::endl;
#endif
for(vertex_descriptor v : vrange_)
{
if(is_border(v, mesh_) || is_constrained(v))
continue;
Point_ref p_query = get(vpmap_, v);
const Point projected = tree.closest_point(p_query);
put(vpmap_, v, projected);
}
}
private:
bool is_constrained(const vertex_descriptor v)
{
return get(vcmap_, v);
}
// check for degenerate or inversed faces
bool does_move_create_bad_faces(const vertex_descriptor v,
const Point& new_pos) const
{
// check for null faces and face inversions
for(halfedge_descriptor main_he : halfedges_around_source(v, mesh_))
{
const halfedge_descriptor prev_he = prev(main_he, mesh_);
const Point_ref lpt = get(vpmap_, target(main_he, mesh_));
const Point_ref rpt = get(vpmap_, source(prev_he, mesh_));
if(traits_.collinear_3_object()(lpt, rpt, new_pos))
return true;
const Point_ref old_pos = get(vpmap_, v);
Vector ov_1 = traits_.construct_vector_3_object()(old_pos, lpt);
Vector ov_2 = traits_.construct_vector_3_object()(old_pos, rpt);
Vector old_n = traits_.construct_cross_product_vector_3_object()(ov_1, ov_2);
Vector nv_1 = traits_.construct_vector_3_object()(new_pos, lpt);
Vector nv_2 = traits_.construct_vector_3_object()(new_pos, rpt);
Vector new_n = traits_.construct_cross_product_vector_3_object()(nv_1, nv_2);
if(!is_positive(traits_.compute_scalar_product_3_object()(old_n, new_n)))
{
#ifdef CGAL_PMP_SMOOTHING_DEBUG
std::cout << "Moving vertex would result in the inversion of a face normal!" << std::endl;
#endif
return true;
}
}
return false;
}
bool does_improve_min_angle_in_star(const vertex_descriptor v,
const Point& new_pos) const
{
// check if the minimum angle of the star has not deteriorated
double old_min_angle = CGAL_PI;
for(halfedge_descriptor main_he : halfedges_around_source(v, mesh_))
{
const Point_ref old_pos = get(vpmap_, v);
const halfedge_descriptor prev_he = prev(main_he, mesh_);
const Point_ref lpt = get(vpmap_, target(main_he, mesh_));
const Point_ref rpt = get(vpmap_, source(prev_he, mesh_));
old_min_angle = (std::min)(old_min_angle,
(std::min)(get_radian_angle(Vector(old_pos, lpt),
Vector(old_pos, rpt), traits_),
(std::min)(get_radian_angle(Vector(lpt, rpt),
Vector(lpt, old_pos), traits_),
get_radian_angle(Vector(rpt, old_pos),
Vector(rpt, lpt), traits_))));
}
for(halfedge_descriptor main_he : halfedges_around_source(v, mesh_))
{
const halfedge_descriptor prev_he = prev(main_he, mesh_);
const Point_ref lpt = get(vpmap_, target(main_he, mesh_));
const Point_ref rpt = get(vpmap_, source(prev_he, mesh_));
if(get_radian_angle(Vector(new_pos, lpt), Vector(new_pos, rpt), traits_) < old_min_angle ||
get_radian_angle(Vector(lpt, rpt), Vector(lpt, new_pos), traits_) < old_min_angle ||
get_radian_angle(Vector(rpt, new_pos), Vector(rpt, lpt), traits_) < old_min_angle)
{
#ifdef CGAL_PMP_SMOOTHING_DEBUG
const Point_ref old_pos = get(vpmap_, v);
std::cout << "deterioration of min angle in the star!" << std::endl;
std::cout << "old/new positions: " << old_pos << " " << new_pos << std::endl;;
std::cout << "old min angle: " << old_min_angle << std::endl;
std::cout << "new angles: " << std::endl;
std::cout << get_radian_angle(Vector(new_pos, lpt), Vector(new_pos, rpt), traits_) << " ";
std::cout << get_radian_angle(Vector(lpt, rpt), Vector(lpt, new_pos), traits_) << " ";
std::cout << get_radian_angle(Vector(rpt, new_pos), Vector(rpt, lpt), traits_) << std::endl;
#endif
return false;
}
}
return true;
}
private:
TriangleMesh& mesh_;
VertexPointMap& vpmap_;
VertexConstraintMap vcmap_;
GeomTraits traits_;
std::vector<vertex_descriptor> vrange_;
};
} // namespace internal
} // namespace Polygon_mesh_processing
} // namespace CGAL
#endif // CGAL_POLYGON_MESH_PROCESSING_INTERNAL_MESH_SMOOTHING_IMPL_H

View File

@ -0,0 +1,136 @@
// Copyright (c) 2018 GeometryFactory (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
// You can redistribute it and/or modify it under the terms of the GNU
// General Public License as published by the Free Software Foundation,
// either version 3 of the License, or (at your option) any later version.
//
// Licensees holding a valid commercial license may use this file in
// accordance with the commercial license agreement provided with the software.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0+
//
//
// Author(s) : Konstantinos Katrioplas (konst.katrioplas@gmail.com)
#ifndef CGAL_POLYGON_MESH_PROCESSING_INTERNAL_SMOOTHING_EVALUATION_H
#define CGAL_POLYGON_MESH_PROCESSING_INTERNAL_SMOOTHING_EVALUATION_H
#include <CGAL/license/Polygon_mesh_processing/meshing_hole_filling.h>
#include <CGAL/Polygon_mesh_processing/measure.h>
#include <boost/graph/graph_traits.hpp>
#include <boost/property_map/property_map.hpp>
#include <set>
#include <vector>
namespace CGAL {
namespace Polygon_mesh_processing {
namespace internal {
template<typename PolygonMesh, typename GeomTraits>
class Quality_evaluator
{
typedef typename GeomTraits::Point_3 Point;
typedef typename GeomTraits::Vector_3 Vector;
typedef typename GeomTraits::Line_3 Line;
typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
typedef typename boost::graph_traits<PolygonMesh>::face_descriptor face_descriptor;
typedef typename boost::property_map<PolygonMesh, CGAL::vertex_point_t>::type VertexPointMap;
typedef typename boost::property_traits<VertexPointMap>::reference Point_ref;
public:
Quality_evaluator(PolygonMesh& pmesh,
const GeomTraits& traits)
: mesh_(pmesh), traits_(traits)
{
std::size_t number_of_triangles = faces(mesh_).size();
angles_.reserve(number_of_triangles * 3);
areas_.reserve(number_of_triangles);
aspect_ratios_.reserve(number_of_triangles);
vpmap_ = get(CGAL::vertex_point, mesh_);
}
void gather_angles()
{
for(halfedge_descriptor hi : halfedges(mesh_))
{
const Point_ref a = get(vpmap_, source(hi, mesh_));
const Point_ref b = get(vpmap_, target(hi, mesh_));
const Point_ref c = get(vpmap_, target(next(hi, mesh_), mesh_));
angles_.push_back(traits_.compute_approximate_angle_3_object()(a, b, c));
}
#ifdef CGAL_PMP_SMOOTHING_VERBOSE
std::cout << "angles_ size = " << angles_.size() << std::endl;
#endif
}
template <typename Stream>
void extract_angles(Stream& output)
{
for(unsigned int i=0; i!=angles_.size(); ++i)
output << angles_[i] << std::endl;
output.close();
}
void measure_areas()
{
for(face_descriptor f : faces(mesh_))
areas_.push_back(face_area(f, mesh_));
#ifdef CGAL_PMP_SMOOTHING_VERBOSE
std::cout << "areas_ size = " << areas_.size() << std::endl;
#endif
}
template <typename Stream>
void extract_areas(Stream& output)
{
for(unsigned int i=0; i!=areas_.size(); ++i)
output << areas_[i] << std::endl;
output.close();
}
void calc_aspect_ratios()
{
for(face_descriptor f : faces(mesh_))
aspect_ratios_.push_back(CGAL::Polygon_mesh_processing::face_aspect_ratio(f, mesh_));
#ifdef CGAL_PMP_SMOOTHING_VERBOSE
std::cout << "aspect_ratios_ size = " << aspect_ratios_.size() << std::endl;
#endif
}
template <typename Stream>
void extract_aspect_ratios(Stream& output)
{
for(unsigned int i=0; i!=aspect_ratios_.size(); ++i)
output << aspect_ratios_[i] << std::endl;
output.close();
}
private:
PolygonMesh& mesh_;
const GeomTraits& traits_;
VertexPointMap vpmap_;
std::vector<double> angles_;
std::vector<double> areas_;
std::vector<double> aspect_ratios_;
};
} // namespace internal
} // namespace Polygon_mesh_processing
} // namespace CGAL
#endif // CGAL_POLYGON_MESH_PROCESSING_INTERNAL_SMOOTHING_EVALUATION_H

View File

@ -0,0 +1,41 @@
// Copyright (c) 2018 GeometryFactory (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
// You can redistribute it and/or modify it under the terms of the GNU
// General Public License as published by the Free Software Foundation,
// either version 3 of the License, or (at your option) any later version.
//
// Licensees holding a valid commercial license may use this file in
// accordance with the commercial license agreement provided with the software.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0+
//
//
// Author(s) : Konstantinos Katrioplas (konst.katrioplas@gmail.com)
#ifndef CGAL_POLYGON_MESH_PROCESSING_INTERNAL_SMOOTHING_HELPERS_H
#define CGAL_POLYGON_MESH_PROCESSING_INTERNAL_SMOOTHING_HELPERS_H
#include <CGAL/license/Polygon_mesh_processing/meshing_hole_filling.h>
#include <CGAL/Polygon_mesh_processing/Weights.h>
#include <CGAL/Polygon_mesh_processing/internal/named_function_params.h>
#include <CGAL/Polygon_mesh_processing/internal/named_params_helper.h>
namespace CGAL {
namespace Polygon_mesh_processing {
namespace internal {
} // internal
} // Polygon_mesh_processing
} // CGAL
#endif //CGAL_POLYGON_MESH_PROCESSING_INTERNAL_SMOOTHING_HELPERS_H

View File

@ -26,12 +26,10 @@
#include <CGAL/disable_warnings.h>
#include <CGAL/assertions.h>
#include <CGAL/boost/graph/iterator.h>
#include <CGAL/boost/graph/helpers.h>
#include <CGAL/boost/graph/properties.h>
#include <boost/unordered_set.hpp>
#include <boost/graph/graph_traits.hpp>
#include <CGAL/Polygon_mesh_processing/internal/named_function_params.h>
#include <CGAL/Polygon_mesh_processing/internal/named_params_helper.h>
#include <CGAL/squared_distance_3.h>
@ -39,6 +37,9 @@
#include <CGAL/Lazy.h> // needed for CGAL::exact(FT)/CGAL::exact(Lazy_exact_nt<T>)
#include <boost/unordered_set.hpp>
#include <boost/graph/graph_traits.hpp>
#include <utility>
#ifdef DOXYGEN_RUNNING
@ -58,7 +59,7 @@ public:
struct type{};
};
/**
/**
* \ingroup measure_grp
* computes the length of an edge of a given polygon mesh.
* The edge is given by one of its halfedges, or the edge itself.
@ -88,58 +89,57 @@ public:
*
* @sa `face_border_length()`
*/
template<typename PolygonMesh,
typename NamedParameters>
template<typename PolygonMesh,
typename NamedParameters>
#ifdef DOXYGEN_RUNNING
FT
FT
#else
typename GetGeomTraits<PolygonMesh, NamedParameters>::type::FT
typename GetGeomTraits<PolygonMesh, NamedParameters>::type::FT
#endif
edge_length(typename boost::graph_traits<PolygonMesh>::halfedge_descriptor h
, const PolygonMesh& pmesh
, const NamedParameters& np)
{
using boost::choose_param;
using boost::get_param;
edge_length(typename boost::graph_traits<PolygonMesh>::halfedge_descriptor h,
const PolygonMesh& pmesh,
const NamedParameters& np)
{
using boost::choose_param;
using boost::get_param;
typename GetVertexPointMap<PolygonMesh, NamedParameters>::const_type
vpm = choose_param(get_param(np, internal_np::vertex_point),
get_const_property_map(CGAL::vertex_point, pmesh));
typename GetVertexPointMap<PolygonMesh, NamedParameters>::const_type
vpm = choose_param(get_param(np, internal_np::vertex_point),
get_const_property_map(CGAL::vertex_point, pmesh));
return CGAL::approximate_sqrt(CGAL::squared_distance(get(vpm, source(h, pmesh)),
get(vpm, target(h, pmesh))));
}
return CGAL::approximate_sqrt(CGAL::squared_distance(get(vpm, source(h, pmesh)),
get(vpm, target(h, pmesh))));
}
template<typename PolygonMesh>
typename CGAL::Kernel_traits<typename property_map_value<PolygonMesh,
CGAL::vertex_point_t>::type>::Kernel::FT
edge_length(typename boost::graph_traits<PolygonMesh>::halfedge_descriptor h
, const PolygonMesh& pmesh)
{
return edge_length(h, pmesh,
CGAL::Polygon_mesh_processing::parameters::all_default());
}
// edge overloads
template<typename PolygonMesh,
typename NamedParameters>
typename GetGeomTraits<PolygonMesh, NamedParameters>::type::FT
edge_length(typename boost::graph_traits<PolygonMesh>::edge_descriptor e
, const PolygonMesh& pmesh
, const NamedParameters& np)
{
return edge_length(halfedge(e,pmesh), pmesh, np);
}
template<typename PolygonMesh>
typename CGAL::Kernel_traits<typename property_map_value<PolygonMesh,
CGAL::vertex_point_t>::type>::Kernel::FT
edge_length(typename boost::graph_traits<PolygonMesh>::halfedge_descriptor h,
const PolygonMesh& pmesh)
{
return edge_length(h, pmesh, CGAL::Polygon_mesh_processing::parameters::all_default());
}
// edge overloads
template<typename PolygonMesh,
typename NamedParameters>
typename GetGeomTraits<PolygonMesh, NamedParameters>::type::FT
edge_length(typename boost::graph_traits<PolygonMesh>::edge_descriptor e,
const PolygonMesh& pmesh,
const NamedParameters& np)
{
return edge_length(halfedge(e, pmesh), pmesh, np);
}
template<typename PolygonMesh>
typename CGAL::Kernel_traits<typename property_map_value<PolygonMesh,
CGAL::vertex_point_t>::type>::Kernel::FT
edge_length(typename boost::graph_traits<PolygonMesh>::edge_descriptor e
, const PolygonMesh& pmesh)
{
return edge_length(halfedge(e,pmesh), pmesh);
}
template<typename PolygonMesh>
typename CGAL::Kernel_traits<typename property_map_value<PolygonMesh,
CGAL::vertex_point_t>::type>::Kernel::FT
edge_length(typename boost::graph_traits<PolygonMesh>::edge_descriptor e,
const PolygonMesh& pmesh)
{
return edge_length(halfedge(e, pmesh), pmesh);
}
/**
/**
* \ingroup measure_grp
* computes the length of the border polyline
* that contains a given halfedge.
@ -170,40 +170,38 @@ public:
*
* @sa `edge_length()`
*/
template<typename PolygonMesh,
typename NamedParameters>
template<typename PolygonMesh,
typename NamedParameters>
#ifdef DOXYGEN_RUNNING
FT
FT
#else
typename GetGeomTraits<PolygonMesh, NamedParameters>::type::FT
typename GetGeomTraits<PolygonMesh, NamedParameters>::type::FT
#endif
face_border_length(
typename boost::graph_traits<PolygonMesh>::halfedge_descriptor h
, const PolygonMesh& pmesh
, const NamedParameters& np)
face_border_length(typename boost::graph_traits<PolygonMesh>::halfedge_descriptor h,
const PolygonMesh& pmesh,
const NamedParameters& np)
{
typename GetGeomTraits<PolygonMesh, NamedParameters>::type::FT result = 0;
for(typename boost::graph_traits<PolygonMesh>::halfedge_descriptor haf : halfedges_around_face(h, pmesh))
{
typename GetGeomTraits<PolygonMesh, NamedParameters>::type::FT result = 0.;
for(typename boost::graph_traits<PolygonMesh>::halfedge_descriptor haf :
halfedges_around_face(h, pmesh))
{
result += edge_length(haf, pmesh, np);
exact(result);
}
return result;
result += edge_length(haf, pmesh, np);
exact(result);
}
template<typename PolygonMesh>
typename CGAL::Kernel_traits<typename property_map_value<PolygonMesh,
CGAL::vertex_point_t>::type>::Kernel::FT
face_border_length(
typename boost::graph_traits<PolygonMesh>::halfedge_descriptor h
, const PolygonMesh& pmesh)
{
return face_border_length(h, pmesh,
CGAL::Polygon_mesh_processing::parameters::all_default());
}
return result;
}
/**
template<typename PolygonMesh>
typename CGAL::Kernel_traits<typename property_map_value<PolygonMesh,
CGAL::vertex_point_t>::type>::Kernel::FT
face_border_length(typename boost::graph_traits<PolygonMesh>::halfedge_descriptor h,
const PolygonMesh& pmesh)
{
return face_border_length(h, pmesh, CGAL::Polygon_mesh_processing::parameters::all_default());
}
/**
* \ingroup measure_grp
* finds the longest border of a given triangulated surface and returns
* a halfedge that is part of this border and the length of this border.
@ -232,53 +230,58 @@ public:
* of `pmesh`
*
*/
template<typename PolygonMesh,
typename NamedParameters>
template<typename PolygonMesh,
typename NamedParameters>
#ifdef DOXYGEN_RUNNING
std::pair<halfedge_descriptor, FT>
std::pair<halfedge_descriptor, FT>
#else
std::pair<typename boost::graph_traits<PolygonMesh>::halfedge_descriptor,
typename GetGeomTraits<PolygonMesh, NamedParameters>::type::FT>
std::pair<typename boost::graph_traits<PolygonMesh>::halfedge_descriptor,
typename GetGeomTraits<PolygonMesh, NamedParameters>::type::FT>
#endif
longest_border(const PolygonMesh& pmesh
, const NamedParameters& np)
{
typedef typename CGAL::Kernel_traits<typename property_map_value<PolygonMesh,
CGAL::vertex_point_t>::type>::Kernel::FT FT;
typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
longest_border(const PolygonMesh& pmesh,
const NamedParameters& np)
{
typedef typename CGAL::Kernel_traits<
typename property_map_value<PolygonMesh, CGAL::vertex_point_t>::type>::Kernel::FT FT;
typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
boost::unordered_set<halfedge_descriptor> visited;
halfedge_descriptor result_halfedge = boost::graph_traits<PolygonMesh>::null_halfedge();
FT result_len = 0;
for(halfedge_descriptor h : halfedges(pmesh)){
if(visited.find(h)== visited.end()){
if(is_border(h,pmesh)){
FT len = 0;
for(halfedge_descriptor haf : halfedges_around_face(h, pmesh)){
len += edge_length(haf, pmesh, np);
visited.insert(haf);
}
if(result_len < len){
result_len = len;
result_halfedge = h;
}
boost::unordered_set<halfedge_descriptor> visited;
halfedge_descriptor result_halfedge = boost::graph_traits<PolygonMesh>::null_halfedge();
FT result_len = 0;
for(halfedge_descriptor h : halfedges(pmesh))
{
if(visited.find(h)== visited.end())
{
if(is_border(h, pmesh))
{
FT len = 0;
for(halfedge_descriptor haf : halfedges_around_face(h, pmesh))
{
len += edge_length(haf, pmesh, np);
visited.insert(haf);
}
if(result_len < len)
{
result_len = len;
result_halfedge = h;
}
}
}
return std::make_pair(result_halfedge, result_len);
}
template<typename PolygonMesh>
std::pair<typename boost::graph_traits<PolygonMesh>::halfedge_descriptor,
typename CGAL::Kernel_traits<typename property_map_value<PolygonMesh,
CGAL::vertex_point_t>::type>::Kernel::FT>
longest_border(const PolygonMesh& pmesh)
{
return longest_border(pmesh,
CGAL::Polygon_mesh_processing::parameters::all_default());
}
return std::make_pair(result_halfedge, result_len);
}
/**
template<typename PolygonMesh>
std::pair<typename boost::graph_traits<PolygonMesh>::halfedge_descriptor,
typename CGAL::Kernel_traits<typename property_map_value<PolygonMesh,
CGAL::vertex_point_t>::type>::Kernel::FT>
longest_border(const PolygonMesh& pmesh)
{
return longest_border(pmesh, CGAL::Polygon_mesh_processing::parameters::all_default());
}
/**
* \ingroup measure_grp
* computes the area of a face of a given
* triangulated surface mesh.
@ -307,49 +310,48 @@ public:
*
* @sa `area()`
*/
template<typename TriangleMesh,
typename CGAL_PMP_NP_TEMPLATE_PARAMETERS>
template<typename TriangleMesh,
typename CGAL_PMP_NP_TEMPLATE_PARAMETERS>
#ifdef DOXYGEN_RUNNING
FT
FT
#else
typename GetGeomTraits<TriangleMesh, CGAL_PMP_NP_CLASS>::type::FT
typename GetGeomTraits<TriangleMesh, CGAL_PMP_NP_CLASS>::type::FT
#endif
face_area(typename boost::graph_traits<TriangleMesh>::face_descriptor f
, const TriangleMesh& tmesh
, const CGAL_PMP_NP_CLASS& np)
{
using boost::choose_param;
using boost::get_param;
face_area(typename boost::graph_traits<TriangleMesh>::face_descriptor f,
const TriangleMesh& tmesh,
const CGAL_PMP_NP_CLASS& np)
{
using boost::choose_param;
using boost::get_param;
CGAL_precondition(boost::graph_traits<TriangleMesh>::null_face() != f);
typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor halfedge_descriptor;
typename GetVertexPointMap<TriangleMesh, CGAL_PMP_NP_CLASS>::const_type
vpm = choose_param(get_param(np, internal_np::vertex_point),
get_const_property_map(CGAL::vertex_point, tmesh));
CGAL_precondition(boost::graph_traits<TriangleMesh>::null_face() != f);
typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor halfedge_descriptor;
halfedge_descriptor hd = halfedge(f, tmesh);
halfedge_descriptor nhd = next(hd, tmesh);
typename GetVertexPointMap<TriangleMesh, CGAL_PMP_NP_CLASS>::const_type
vpm = choose_param(get_param(np, internal_np::vertex_point),
get_const_property_map(CGAL::vertex_point, tmesh));
typename GetGeomTraits<TriangleMesh, CGAL_PMP_NP_CLASS>::type traits;
halfedge_descriptor hd = halfedge(f, tmesh);
halfedge_descriptor nhd = next(hd, tmesh);
return approximate_sqrt(
traits.compute_squared_area_3_object()(get(vpm, source(hd, tmesh)),
get(vpm, target(hd, tmesh)),
get(vpm, target(nhd, tmesh))));
}
typename GetGeomTraits<TriangleMesh, CGAL_PMP_NP_CLASS>::type traits;
template<typename TriangleMesh>
typename CGAL::Kernel_traits<typename property_map_value<TriangleMesh,
CGAL::vertex_point_t>::type>::Kernel::FT
face_area(typename boost::graph_traits<TriangleMesh>::face_descriptor f
, const TriangleMesh& tmesh)
{
return face_area(f, tmesh,
CGAL::Polygon_mesh_processing::parameters::all_default());
}
return approximate_sqrt(traits.compute_squared_area_3_object()(get(vpm, source(hd, tmesh)),
get(vpm, target(hd, tmesh)),
get(vpm, target(nhd, tmesh))));
}
/**
template<typename TriangleMesh>
typename CGAL::Kernel_traits<typename property_map_value<TriangleMesh,
CGAL::vertex_point_t>::type>::Kernel::FT
face_area(typename boost::graph_traits<TriangleMesh>::face_descriptor f,
const TriangleMesh& tmesh)
{
return face_area(f, tmesh, CGAL::Polygon_mesh_processing::parameters::all_default());
}
/**
* \ingroup measure_grp
* computes the area of a range of faces of a given
* triangulated surface mesh.
@ -383,37 +385,38 @@ public:
*
* @sa `face_area()`
*/
template<typename FaceRange,
typename TriangleMesh,
typename CGAL_PMP_NP_TEMPLATE_PARAMETERS>
template<typename FaceRange,
typename TriangleMesh,
typename CGAL_PMP_NP_TEMPLATE_PARAMETERS>
#ifdef DOXYGEN_RUNNING
FT
FT
#else
typename GetGeomTraits<TriangleMesh, CGAL_PMP_NP_CLASS>::type::FT
typename GetGeomTraits<TriangleMesh, CGAL_PMP_NP_CLASS>::type::FT
#endif
area(FaceRange face_range
, const TriangleMesh& tmesh
, const CGAL_PMP_NP_CLASS& np)
area(FaceRange face_range,
const TriangleMesh& tmesh,
const CGAL_PMP_NP_CLASS& np)
{
typedef typename boost::graph_traits<TriangleMesh>::face_descriptor face_descriptor;
typename GetGeomTraits<TriangleMesh, CGAL_PMP_NP_CLASS>::type::FT result = 0;
for(face_descriptor f : face_range)
{
typedef typename boost::graph_traits<TriangleMesh>::face_descriptor face_descriptor;
typename GetGeomTraits<TriangleMesh, CGAL_PMP_NP_CLASS>::type::FT result = 0.;
for(face_descriptor f : face_range)
{
result += face_area(f, tmesh, np);
exact(result);
}
return result;
result += face_area(f, tmesh, np);
exact(result);
}
template<typename FaceRange, typename TriangleMesh>
typename GetGeomTraits<TriangleMesh>::type::FT
area(FaceRange face_range, const TriangleMesh& tmesh)
{
return area(face_range, tmesh,
CGAL::Polygon_mesh_processing::parameters::all_default());
}
return result;
}
/**
template<typename FaceRange, typename TriangleMesh>
typename GetGeomTraits<TriangleMesh>::type::FT
area(FaceRange face_range, const TriangleMesh& tmesh)
{
return area(face_range, tmesh, CGAL::Polygon_mesh_processing::parameters::all_default());
}
/**
* \ingroup measure_grp
* computes the surface area of a triangulated surface mesh.
*
@ -442,28 +445,27 @@ public:
*
* @sa `face_area()`
*/
template<typename TriangleMesh
, typename CGAL_PMP_NP_TEMPLATE_PARAMETERS>
template<typename TriangleMesh,
typename CGAL_PMP_NP_TEMPLATE_PARAMETERS>
#ifdef DOXYGEN_RUNNING
FT
FT
#else
typename GetGeomTraits<TriangleMesh, CGAL_PMP_NP_CLASS>::type::FT
typename GetGeomTraits<TriangleMesh, CGAL_PMP_NP_CLASS>::type::FT
#endif
area(const TriangleMesh& tmesh, const CGAL_PMP_NP_CLASS& np)
{
return area(faces(tmesh), tmesh, np);
}
area(const TriangleMesh& tmesh, const CGAL_PMP_NP_CLASS& np)
{
return area(faces(tmesh), tmesh, np);
}
template<typename TriangleMesh>
typename CGAL::Kernel_traits<typename property_map_value<TriangleMesh,
CGAL::vertex_point_t>::type>::Kernel::FT
area(const TriangleMesh& tmesh)
{
return area(faces(tmesh), tmesh
, CGAL::Polygon_mesh_processing::parameters::all_default());
}
template<typename TriangleMesh>
typename CGAL::Kernel_traits<typename property_map_value<TriangleMesh,
CGAL::vertex_point_t>::type>::Kernel::FT
area(const TriangleMesh& tmesh)
{
return area(faces(tmesh), tmesh, CGAL::Polygon_mesh_processing::parameters::all_default());
}
/**
/**
* \ingroup measure_grp
* computes the volume of the domain bounded by
* a closed triangulated surface mesh.
@ -489,14 +491,15 @@ public:
* or the geometric traits class deduced from the point property map
* of `tmesh`.
*/
template<typename TriangleMesh
, typename CGAL_PMP_NP_TEMPLATE_PARAMETERS>
template<typename TriangleMesh,
typename CGAL_PMP_NP_TEMPLATE_PARAMETERS>
#ifdef DOXYGEN_RUNNING
FT
FT
#else
typename GetGeomTraits<TriangleMesh, CGAL_PMP_NP_CLASS>::type::FT
typename GetGeomTraits<TriangleMesh, CGAL_PMP_NP_CLASS>::type::FT
#endif
volume(const TriangleMesh& tmesh, const CGAL_PMP_NP_CLASS& np)
volume(const TriangleMesh& tmesh,
const CGAL_PMP_NP_CLASS& np)
{
CGAL_assertion(is_triangle_mesh(tmesh));
CGAL_assertion(is_closed(tmesh));
@ -505,37 +508,148 @@ volume(const TriangleMesh& tmesh, const CGAL_PMP_NP_CLASS& np)
using boost::get_param;
typename GetVertexPointMap<TriangleMesh, CGAL_PMP_NP_CLASS>::const_type
vpm = choose_param(get_param(np, internal_np::vertex_point),
get_const_property_map(CGAL::vertex_point, tmesh));
typename GetGeomTraits<TriangleMesh, CGAL_PMP_NP_CLASS>::type::Point_3
origin(0, 0, 0);
vpm = choose_param(get_param(np, internal_np::vertex_point),
get_const_property_map(CGAL::vertex_point, tmesh));
typename GetGeomTraits<TriangleMesh, CGAL_PMP_NP_CLASS>::type::Point_3 origin(0, 0, 0);
typedef typename boost::graph_traits<TriangleMesh>::face_descriptor face_descriptor;
typename GetGeomTraits<TriangleMesh, CGAL_PMP_NP_CLASS>::type::FT volume = 0.;
typename GetGeomTraits<TriangleMesh, CGAL_PMP_NP_CLASS>::type::FT volume = 0;
typename CGAL::Kernel_traits<typename property_map_value<TriangleMesh,
CGAL::vertex_point_t>::type>::Kernel::Compute_volume_3 cv3;
CGAL::vertex_point_t>::type>::Kernel::Compute_volume_3 cv3;
for(face_descriptor f : faces(tmesh))
{
volume += cv3(origin,
get(vpm, target(halfedge(f, tmesh), tmesh)),
get(vpm, target(next(halfedge(f, tmesh), tmesh), tmesh)),
get(vpm, target(prev(halfedge(f, tmesh), tmesh), tmesh)));
get(vpm, target(halfedge(f, tmesh), tmesh)),
get(vpm, target(next(halfedge(f, tmesh), tmesh), tmesh)),
get(vpm, target(prev(halfedge(f, tmesh), tmesh), tmesh)));
exact(volume);
}
return volume;
}
template<typename TriangleMesh>
typename CGAL::Kernel_traits<typename property_map_value<TriangleMesh,
CGAL::vertex_point_t>::type>::Kernel::FT
volume(const TriangleMesh& tmesh)
template<typename TriangleMesh>
typename CGAL::Kernel_traits<typename property_map_value<TriangleMesh,
CGAL::vertex_point_t>::type>::Kernel::FT
volume(const TriangleMesh& tmesh)
{
return volume(tmesh, CGAL::Polygon_mesh_processing::parameters::all_default());
}
/**
* \ingroup measure_grp
* computes the aspect ratio of a face of a given triangulated surface mesh.
*
* @tparam TriangleMesh a model of `HalfedgeGraph`
* @tparam NamedParameters a sequence of \ref pmp_namedparameters "Named Parameters"
*
* @param f the face of which the aspect ratio is computed
* @param tmesh the triangulated surface mesh to which `f` belongs
* @param np optional sequence of \ref pmp_namedparameters "Named Parameters" among the ones listed below
*
* \cgalNamedParamsBegin
* \cgalParamBegin{vertex_point_map} the property map with the points associated to the vertices of `pmesh`.
* If this parameter is omitted, an internal property map for
* `CGAL::vertex_point_t` must be available in `TriangleMesh`\cgalParamEnd
* \cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `Kernel`\cgalParamEnd
* \cgalNamedParamsEnd
*
* @pre `f != boost::graph_traits<TriangleMesh>::%null_face()`
*
* @return the aspect ratio of `f`. The return type `FT` is a number type. It is
* either deduced from the `geom_traits` \ref pmp_namedparameters "Named Parameters" if provided,
* or the geometric traits class deduced from the point property map of `tmesh`.
*
* @warning
*
*/
template<typename TriangleMesh,
typename CGAL_PMP_NP_TEMPLATE_PARAMETERS>
#ifdef DOXYGEN_RUNNING
FT
#else
typename GetGeomTraits<TriangleMesh, CGAL_PMP_NP_CLASS>::type::FT
#endif
face_aspect_ratio(typename boost::graph_traits<TriangleMesh>::face_descriptor f,
const TriangleMesh& tmesh,
const CGAL_PMP_NP_CLASS& np)
{
CGAL_precondition(is_triangle(f, tmesh));
typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor halfedge_descriptor;
typedef typename GetGeomTraits<TriangleMesh, CGAL_PMP_NP_CLASS>::type Geom_traits;
typedef typename Geom_traits::FT FT;
using boost::choose_param;
using boost::get_param;
typename GetVertexPointMap<TriangleMesh, CGAL_PMP_NP_CLASS>::const_type
vpm = choose_param(get_param(np, internal_np::vertex_point),
get_const_property_map(CGAL::vertex_point, tmesh));
halfedge_descriptor h = halfedge(f, tmesh);
Geom_traits gt = choose_param(get_param(np, internal_np::geom_traits), Geom_traits());
#if 0
const FT sq_triangle_area = gt.compute_squared_area_3_object()(get(vpm, source(h, tmesh)),
get(vpm, target(h, tmesh)),
get(vpm, target(next(h, tmesh), tmesh)));
const FT sq_d12 = gt.compute_squared_distance_2_object()(get(vpm, source(h, tmesh)),
get(vpm, target(h, tmesh)));
const FT sq_d13 = gt.compute_squared_distance_2_object()(get(vpm, source(h, tmesh)),
get(vpm, target(next(h, tmesh), tmesh)));
const FT sq_d23 = gt.compute_squared_distance_2_object()(get(vpm, target(h, tmesh)),
get(vpm, target(next(h, tmesh), tmesh)));
const FT min_sq_d123 = (std::min)(sq_d12, (std::min)(sq_d13, sq_d23));
const FT aspect_ratio = 4*sq_triangle_area*min_sq_d123 / (sq_d12*sq_d13*sq_d23);
#else // below requires SQRT
typedef typename Geom_traits::Line_3 Line_3;
FT sq_max_edge_length = gt.compute_squared_distance_3_object()(get(vpm, source(h, tmesh)),
get(vpm, target(h, tmesh)));
FT sq_min_alt = gt.compute_squared_distance_3_object()(get(vpm, target(next(h, tmesh), tmesh)),
Line_3(get(vpm, source(h, tmesh)),
get(vpm, target(h, tmesh))));
h = next(h, tmesh);
for(int i=1; i<3; ++i)
{
return volume(tmesh,
CGAL::Polygon_mesh_processing::parameters::all_default());
FT sq_edge_length = gt.compute_squared_distance_3_object()(get(vpm, source(h, tmesh)),
get(vpm, target(h, tmesh)));
FT sq_alt = gt.compute_squared_distance_3_object()(get(vpm, target(next(h, tmesh), tmesh)),
Line_3(get(vpm, source(h, tmesh)),
get(vpm, target(h, tmesh))));
if(sq_alt < sq_min_alt)
sq_min_alt = sq_alt;
if(sq_edge_length > sq_max_edge_length)
sq_max_edge_length = sq_edge_length;
h = next(h, tmesh);
}
CGAL_assertion(sq_min_alt > 0);
const FT aspect_ratio = CGAL::approximate_sqrt(sq_max_edge_length / sq_min_alt);
#endif
return aspect_ratio;
}
template<typename TriangleMesh>
typename CGAL::Kernel_traits<typename property_map_value<TriangleMesh,
CGAL::vertex_point_t>::type>::Kernel::FT
face_aspect_ratio(typename boost::graph_traits<TriangleMesh>::face_descriptor f,
const TriangleMesh& tmesh)
{
return face_aspect_ratio(f, tmesh, CGAL::Polygon_mesh_processing::parameters::all_default());
}
/**
* \ingroup measure_grp
* computes the centroid of a volume bounded by
@ -558,17 +672,16 @@ volume(const TriangleMesh& tmesh, const CGAL_PMP_NP_CLASS& np)
*
* @return the centroid of the domain bounded by `tmesh`.
*/
template<typename TriangleMesh
, typename CGAL_PMP_NP_TEMPLATE_PARAMETERS>
template<typename TriangleMesh, typename CGAL_PMP_NP_TEMPLATE_PARAMETERS>
#ifdef DOXYGEN_RUNNING
Point_3
Point_3
#else
typename GetGeomTraits<TriangleMesh, CGAL_PMP_NP_CLASS>::type::Point_3
typename GetGeomTraits<TriangleMesh, CGAL_PMP_NP_CLASS>::type::Point_3
#endif
centroid(const TriangleMesh& tmesh, const CGAL_PMP_NP_CLASS& np)
{
// See: http://www2.imperial.ac.uk/~rn/centroid.pdf
CGAL_assertion(is_triangle_mesh(tmesh));
CGAL_assertion(is_closed(tmesh));
@ -578,18 +691,20 @@ centroid(const TriangleMesh& tmesh, const CGAL_PMP_NP_CLASS& np)
typedef typename GetVertexPointMap<TriangleMesh, CGAL_PMP_NP_CLASS>::const_type Vpm;
Vpm vpm = choose_param(get_param(np, internal_np::vertex_point),
get_const_property_map(CGAL::vertex_point, tmesh));
typedef typename GetGeomTraits<TriangleMesh, CGAL_PMP_NP_CLASS>::type Kernel;
typedef typename Kernel::Point_3 Point_3;
typedef typename Kernel::Vector_3 Vector_3;
typedef typename Kernel::Construct_translated_point_3 Construct_translated_point_3;
typedef typename Kernel::Construct_vector_3 Construct_vector_3;
typedef typename Kernel::Construct_normal_3 Construct_normal_3;
typedef typename Kernel::Compute_scalar_product_3 Scalar_product;
typedef typename Kernel::Construct_scaled_vector_3 Scale;
typedef typename Kernel::Construct_sum_of_vectors_3 Sum;
typedef typename boost::graph_traits<TriangleMesh>::face_descriptor face_descriptor;
typedef typename Kernel::FT FT;
FT volume = 0.;
typedef typename Kernel::Point_3 Point_3;
typedef typename Kernel::Vector_3 Vector_3;
typedef typename Kernel::Construct_translated_point_3 Construct_translated_point_3;
typedef typename Kernel::Construct_vector_3 Construct_vector_3;
typedef typename Kernel::Construct_normal_3 Construct_normal_3;
typedef typename Kernel::Compute_scalar_product_3 Scalar_product;
typedef typename Kernel::Construct_scaled_vector_3 Scale;
typedef typename Kernel::Construct_sum_of_vectors_3 Sum;
typedef typename boost::graph_traits<TriangleMesh>::face_descriptor face_descriptor;
typedef typename Kernel::FT FT;
FT volume = 0;
Vector_3 centroid(NULL_VECTOR);
@ -600,13 +715,14 @@ typedef typename Kernel::FT FT;
Scale scale;
Sum sum;
for(face_descriptor fd : faces(tmesh)){
for(face_descriptor fd : faces(tmesh))
{
const Point_3& p = get(vpm, target(halfedge(fd, tmesh), tmesh));
const Point_3& q = get(vpm, target(next(halfedge(fd, tmesh), tmesh), tmesh));
const Point_3& r = get(vpm, target(prev(halfedge(fd, tmesh), tmesh), tmesh));
Vector_3 vp = vector(ORIGIN, p),
vq = vector(ORIGIN, q),
vr = vector(ORIGIN, r);
vq = vector(ORIGIN, q),
vr = vector(ORIGIN, r);
Vector_3 n = normal(p, q, r);
volume += (scalar_product(n,vp))/FT(6);
n = scale(n, FT(1)/FT(24));
@ -617,26 +733,25 @@ typedef typename Kernel::FT FT;
v3 = sum(v3, Vector_3(square(v2.x()), square(v2.y()), square(v2.z())));
v2 = sum(vp, vr);
v3 = sum(v3, Vector_3(square(v2.x()), square(v2.y()), square(v2.z())));
centroid = sum(centroid, Vector_3(n.x() * v3.x(), n.y() * v3.y(), n.z() * v3.z()));
centroid = sum(centroid, Vector_3(n.x() * v3.x(), n.y() * v3.y(), n.z() * v3.z()));
}
centroid = scale(centroid, FT(1)/(FT(2)*volume));
return point(ORIGIN, centroid);
}
template<typename TriangleMesh>
typename CGAL::Kernel_traits<typename property_map_value<TriangleMesh,
CGAL::vertex_point_t>::type>::Kernel::Point_3
CGAL::vertex_point_t>::type>::Kernel::Point_3
centroid(const TriangleMesh& tmesh)
{
return centroid(tmesh,
CGAL::Polygon_mesh_processing::parameters::all_default());
}
}
return centroid(tmesh, CGAL::Polygon_mesh_processing::parameters::all_default());
}
} // namespace Polygon_mesh_processing
} // namespace CGAL
#include <CGAL/enable_warnings.h>
#endif // CGAL_POLYGON_MESH_PROCESSING_MEASURE_H

View File

@ -0,0 +1,300 @@
// Copyright (c) 2018 GeometryFactory (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
// You can redistribute it and/or modify it under the terms of the GNU
// General Public License as published by the Free Software Foundation,
// either version 3 of the License, or (at your option) any later version.
//
// Licensees holding a valid commercial license may use this file in
// accordance with the commercial license agreement provided with the software.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0+
//
//
// Author(s) : Mael Rouxel-Labbé
// Konstantinos Katrioplas (konst.katrioplas@gmail.com)
#ifndef CGAL_POLYGON_MESH_PROCESSING_SMOOTH_MESH_H
#define CGAL_POLYGON_MESH_PROCESSING_SMOOTH_MESH_H
#include <CGAL/license/Polygon_mesh_processing/meshing_hole_filling.h>
#include <CGAL/Polygon_mesh_processing/internal/Smoothing/mesh_smoothing_impl.h>
#include <CGAL/Polygon_mesh_processing/internal/Smoothing/smoothing_evaluation.h>
#include <CGAL/Polygon_mesh_processing/self_intersections.h>
#include <CGAL/boost/graph/named_function_params.h>
#include <CGAL/boost/graph/named_params_helper.h>
#include <CGAL/property_map.h>
namespace CGAL {
namespace Polygon_mesh_processing {
/*!
* \ingroup PMP_meshing_grp
*
* \short smoothes a triangulated region of a polygon mesh.
*
* This function attempts to make the triangle angle and area distributions as uniform as possible
* by moving (non-constrained) vertices.
*
* Angle-based smoothing does not change the combinatorial information of the mesh. Area-based smoothing
* might change the combinatorial information, unless specified otherwise. It is also possible
* to make the smoothing algorithm "safer" by rejecting moves that, when applied, would worsen the
* quality of the mesh, e.g. that would decrease the value of the smallest angle around a vertex or
* create self-intersections.
*
* Optionally, the points are reprojected after each iteration.
*
* @tparam TriangleMesh model of `MutableFaceGraph`.
* @tparam FaceRange range of `boost::graph_traits<TriangleMesh>::%face_descriptor`,
model of `Range`. Its iterator type is `ForwardIterator`.
* @tparam NamedParameters a sequence of \ref pmp_namedparameters "Named Parameters"
*
* @param tmesh a polygon mesh with triangulated surface patches to be smoothed.
* @param faces the range of triangular faces defining one or several surface patches to be smoothed.
* @param np optional sequence of \ref pmp_namedparameters "Named Parameters" among the ones listed below.
*
* \cgalNamedParamsBegin
* \cgalParamBegin{use_angle_smoothing} Boolean value to indicate whether angle-based smoothing should be used.
* %Default is `true`.
* \cgalParamEnd
* \cgalParamBegin{use_area_smoothing} Boolean value to indicate whether area-based smoothing should be used.
* %Default is `true`.
* \cgalParamEnd
* \cgalParamBegin{number_of_iterations} the number of iterations for the
* sequence of the smoothing iterations performed (default is 1).
* \cgalParamEnd
* \cgalParamBegin{use_safety_constraints} if `true`, vertex moves that would worsen the mesh
* are ignored. %Default is `false`.
* \cgalParamEnd
* \cgalParamBegin{use_Delaunay_flips} if `true` (default value), area-based smoothing will be completed
* by a phase of Delaunay-based edge-flips to prevent the creation of elongated triangles.
* \cgalParamEnd
* \cgalParamBegin{do_project} if `true` (default value), points are projected onto the initial surface
* after each iteration.
* \cgalParamEnd
* \cgalParamBegin{vertex_is_constrained_map} a property map containing the
* constrained-or-not status of each vertex of `tmesh`. A constrained vertex
* cannot be modified at all during smoothing.
* \cgalParamEnd
* \cgalParamBegin{edge_is_constrained_map} a property map, model of `ReadWritePropertyMap`, containing the
* constrained-or-not status of each edge of `tmesh`. A constrained edge cannot be flipped.
* \cgalParamEnd
* \cgalParamBegin{vertex_point_map} the property map, model of `ReadWritePropertyMap`, with the points
* associated to the vertices of `tmesh`.
* \cgalParamEnd
* \cgalParamBegin{geom_traits} a geometric traits class instance, model of `Kernel`.
* Exact constructions kernels are not supported by this function.
* \cgalParamEnd
* \cgalNamedParamsEnd
*
* @warning The third party libraries \link thirdpartyCeres Ceres \endlink (and \link installation_eigen Eigen \endlink)
* are required to use the area-based smoothing.
*
* @pre `tmesh` does not contain any degenerate faces
*/
template<typename TriangleMesh, typename FaceRange, typename NamedParameters>
void smooth_mesh(const FaceRange& faces,
TriangleMesh& tmesh,
const NamedParameters& np)
{
typedef typename boost::graph_traits<TriangleMesh>::vertex_descriptor vertex_descriptor;
typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor halfedge_descriptor;
typedef typename boost::graph_traits<TriangleMesh>::edge_descriptor edge_descriptor;
typedef typename boost::graph_traits<TriangleMesh>::face_descriptor face_descriptor;
typedef typename GetGeomTraits<TriangleMesh, NamedParameters>::type GeomTraits;
typedef typename GetVertexPointMap<TriangleMesh, NamedParameters>::type VertexPointMap;
typedef typename boost::lookup_named_param_def<internal_np::vertex_is_constrained_t,
NamedParameters,
Constant_property_map<vertex_descriptor, bool> // default
> ::type VCMap;
typedef typename boost::lookup_named_param_def<internal_np::edge_is_constrained_t,
NamedParameters,
Constant_property_map<edge_descriptor, bool> // default
> ::type ECMap;
typedef internal::Area_smoother<TriangleMesh, VertexPointMap, GeomTraits> Area_optimizer;
typedef internal::Mesh_smoother<Area_optimizer, TriangleMesh,
VertexPointMap, VCMap, GeomTraits> Area_smoother;
typedef internal::Delaunay_edge_flipper<TriangleMesh, VertexPointMap,
ECMap, GeomTraits> Delaunay_flipper;
typedef internal::Angle_smoother<TriangleMesh, VertexPointMap, GeomTraits> Angle_optimizer;
typedef internal::Mesh_smoother<Angle_optimizer, TriangleMesh,
VertexPointMap, VCMap, GeomTraits> Angle_smoother;
typedef typename GeomTraits::Triangle_3 Triangle;
typedef std::vector<Triangle> Triangle_container;
typedef CGAL::AABB_triangle_primitive<GeomTraits,
typename Triangle_container::iterator> AABB_Primitive;
typedef CGAL::AABB_traits<GeomTraits, AABB_Primitive> AABB_Traits;
typedef CGAL::AABB_tree<AABB_Traits> Tree;
if(std::begin(faces) == std::end(faces))
return;
using boost::choose_param;
using boost::get_param;
// named parameters
GeomTraits gt = choose_param(get_param(np, internal_np::geom_traits),
GeomTraits());
VertexPointMap vpmap = choose_param(get_param(np, internal_np::vertex_point),
get_property_map(CGAL::vertex_point, tmesh));
VCMap vcmap = choose_param(get_param(np, internal_np::vertex_is_constrained),
Constant_property_map<vertex_descriptor, bool>());
const bool use_angle_smoothing = choose_param(get_param(np, internal_np::use_angle_smoothing), true);
const bool use_area_smoothing = choose_param(get_param(np, internal_np::use_area_smoothing), true);
if(!use_angle_smoothing && !use_area_smoothing)
std::cerr << "Warning: called PMP::smooth_mesh() but no smoothing method is being used" << std::endl;
std::size_t nb_iterations = choose_param(get_param(np, internal_np::number_of_iterations), 1);
const bool do_project = choose_param(get_param(np, internal_np::do_project), true);
const bool use_safety_constraints = choose_param(get_param(np, internal_np::use_safety_constraints), true);
const bool use_Delaunay_flips = choose_param(get_param(np, internal_np::use_Delaunay_flips), true);
// Construct the AABB tree (if needed for reprojection)
std::vector<Triangle> input_triangles;
if(do_project)
{
input_triangles.reserve(faces.size());
for(face_descriptor f : faces)
{
halfedge_descriptor h = halfedge(f, tmesh);
input_triangles.push_back(gt.construct_triangle_3_object()(get(vpmap, source(h, tmesh)),
get(vpmap, target(h, tmesh)),
get(vpmap, target(next(h, tmesh), tmesh))));
}
}
Tree aabb_tree(input_triangles.begin(), input_triangles.end());
aabb_tree.accelerate_distance_queries();
// Setup the working ranges and check some preconditions
Area_smoother area_smoother(tmesh, vpmap, vcmap, gt);
Angle_smoother angle_smoother(tmesh, vpmap, vcmap, gt);
if(use_area_smoothing)
area_smoother.init_smoothing(faces);
if(use_angle_smoothing)
angle_smoother.init_smoothing(faces);
for(std::size_t i=0; i<nb_iterations; ++i)
{
if(use_area_smoothing)
{
// First apply area smoothing...
area_smoother.optimize(use_safety_constraints /*check for bad faces*/,
false /*apply moves as soon as they're calculated*/,
false /*do not enforce a minimum angle improvement*/);
if(do_project)
{
if(use_safety_constraints && does_self_intersect(tmesh))
{
#ifdef CGAL_PMP_SMOOTHING_VERBOSE
std::cerr << "Cannot re-project as there are self-intersections in the mesh!\n";
#endif
break;
}
area_smoother.project_to_surface(aabb_tree);
}
if(use_Delaunay_flips)
{
ECMap ecmap = choose_param(get_param(np, internal_np::edge_is_constrained),
Constant_property_map<edge_descriptor, bool>(false));
Delaunay_flipper delaunay_flipper(tmesh, vpmap, ecmap, gt);
delaunay_flipper(faces);
}
}
// ... then angle smoothing
if(use_angle_smoothing)
{
angle_smoother.optimize(use_safety_constraints /*check for bad faces*/,
true /*apply all moves at once*/,
use_safety_constraints /*check if the min angle is improved*/);
if(do_project)
{
if(use_safety_constraints && does_self_intersect(tmesh))
{
#ifdef CGAL_PMP_SMOOTHING_VERBOSE
std::cerr << "Can't do re-projection, there are self-intersections in the mesh!\n";
#endif
break;
}
angle_smoother.project_to_surface(aabb_tree);
}
}
}
}
template <typename FaceRange, typename TriangleMesh>
void smooth_mesh(const FaceRange& face_range, TriangleMesh& tmesh)
{
smooth_mesh(face_range, tmesh, parameters::all_default());
}
template <typename TriangleMesh, typename CGAL_PMP_NP_TEMPLATE_PARAMETERS>
void smooth_mesh(TriangleMesh& tmesh, const CGAL_PMP_NP_CLASS& np)
{
smooth_mesh(faces(tmesh), tmesh, np);
}
template<typename TriangleMesh>
void smooth_mesh(TriangleMesh& tmesh)
{
smooth_mesh(faces(tmesh), tmesh, parameters::all_default());
}
///\cond SKIP_IN_MANUAL
template<typename TriangleMesh, typename GeomTraits, typename Stream>
void angles_evaluation(TriangleMesh& tmesh, GeomTraits traits, Stream& output)
{
internal::Quality_evaluator<TriangleMesh, GeomTraits> evaluator(tmesh, traits);
evaluator.gather_angles();
evaluator.extract_angles(output);
}
template<typename TriangleMesh, typename GeomTraits, typename Stream>
void areas_evaluation(TriangleMesh& tmesh, GeomTraits traits, Stream& output)
{
internal::Quality_evaluator<TriangleMesh, GeomTraits> evaluator(tmesh, traits);
evaluator.measure_areas();
evaluator.extract_areas(output);
}
template<typename TriangleMesh, typename GeomTraits, typename Stream>
void aspect_ratio_evaluation(TriangleMesh& tmesh, GeomTraits traits, Stream& output)
{
internal::Quality_evaluator<TriangleMesh, GeomTraits> evaluator(tmesh, traits);
evaluator.calc_aspect_ratios();
evaluator.extract_aspect_ratios(output);
}
///\endcond SKIP_IN_MANUAL
} // namespace Polygon_mesh_processing
} // namespace CGAL
#endif // CGAL_POLYGON_MESH_PROCESSING_SMOOTH_MESH_H

View File

@ -0,0 +1,211 @@
// Copyright (c) 2018-2019 GeometryFactory (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
// You can redistribute it and/or modify it under the terms of the GNU
// General Public License as published by the Free Software Foundation,
// either version 3 of the License, or (at your option) any later version.
//
// Licensees holding a valid commercial license may use this file in
// accordance with the commercial license agreement provided with the software.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0+
//
//
// Author(s) : Mael Rouxel-Labbé
// Konstantinos Katrioplas (konst.katrioplas@gmail.com)
#ifndef CGAL_POLYGON_MESH_PROCESSING_SMOOTH_SHAPE_H
#define CGAL_POLYGON_MESH_PROCESSING_SMOOTH_SHAPE_H
#include <CGAL/license/Polygon_mesh_processing/meshing_hole_filling.h>
#if defined(CGAL_EIGEN3_ENABLED)
#include <CGAL/Eigen_solver_traits.h>
#endif
#include <CGAL/Polygon_mesh_processing/internal/Smoothing/curvature_flow_impl.h>
#include <CGAL/utility.h>
#include <CGAL/property_map.h>
#ifdef CGAL_PMP_SMOOTHING_OUTPUT_INTERMEDIARY_STEPS
#include <fstream>
#include <sstream>
#endif
namespace CGAL {
namespace Polygon_mesh_processing {
/*!
* \ingroup PMP_meshing_grp
* smooths the overall shape of the mesh by using the mean curvature flow.
* The effect depends on the curvature of each area and on a time step which
* represents the amount by which vertices are allowed to move.
* The result conformally maps the initial surface to a sphere.
*
* @tparam TriangleMesh model of `MutableFaceGraph`.
* @tparam FaceRange range of `boost::graph_traits<TriangleMesh>::%face_descriptor`,
* model of `Range`. Its iterator type is `ForwardIterator`.
* @tparam NamedParameters a sequence of \ref pmp_namedparameters "Named Parameters"
*
* @param tmesh a polygon mesh with triangulated surface patches to be smoothed.
* @param faces the range of triangular faces defining one or several surface patches to be smoothed.
* @param time a time step that corresponds to the speed by which the surface is smoothed.
* A larger time step results in faster convergence but details may be distorted to have a larger extent
* compared to more iterations with a smaller step. Typical values scale in the interval (1e-6, 1].
* @param np optional sequence of \ref pmp_namedparameters "Named Parameters" among the ones listed below.
*
* \cgalNamedParamsBegin
* \cgalParamBegin{vertex_point_map} the property map with the points associated
* to the vertices of `tmesh`. Instance of a class model of `ReadWritePropertyMap`.
* \cgalParamEnd
* \cgalParamBegin{geom_traits} a geometric traits class instance, model of `Kernel`.
* \cgalParamEnd
* \cgalParamBegin{vertex_is_constrained_map} a property map containing the
* constrained-or-not status of each vertex of `tmesh`. A constrained vertex
* cannot be modified at all during smoothing.
* \cgalParamEnd
* \cgalParamBegin{number_of_iterations} the number of iterations for the
* sequence of the smoothing iterations performed. Each iteration is performed
* with the given time step.
* \cgalParamEnd
* \cgalParamBegin{sparse_linear_solver} an instance of the sparse linear solver used for smoothing \cgalParamEnd
* \cgalParamEnd
* \cgalNamedParamsEnd
*
* @warning This function involves linear algebra, that is computed using a non-exact floating-point arithmetic.
*/
template<typename TriangleMesh, typename FaceRange, typename NamedParameters>
void smooth_shape(const FaceRange& faces,
TriangleMesh& tmesh,
const double time,
const NamedParameters& np)
{
if(std::begin(faces) == std::end(faces))
return;
typedef typename boost::graph_traits<TriangleMesh>::vertex_descriptor vertex_descriptor;
typedef typename GetGeomTraits<TriangleMesh, NamedParameters>::type GeomTraits;
typedef typename GetVertexPointMap<TriangleMesh, NamedParameters>::type VertexPointMap;
typedef typename boost::lookup_named_param_def<
internal_np::vertex_is_constrained_t,
NamedParameters,
Constant_property_map<vertex_descriptor, bool> >::type VCMap;
using boost::choose_param;
using boost::get_param;
GeomTraits gt = choose_param(get_param(np, internal_np::geom_traits), GeomTraits());
VertexPointMap vpmap = choose_param(get_param(np, internal_np::vertex_point),
get_property_map(CGAL::vertex_point, tmesh));
VCMap vcmap = choose_param(get_param(np, internal_np::vertex_is_constrained),
Constant_property_map<vertex_descriptor, bool>(false));
const std::size_t nb_iterations = choose_param(get_param(np, internal_np::number_of_iterations), 1);
#if defined(CGAL_EIGEN3_ENABLED)
#if EIGEN_VERSION_AT_LEAST(3,2,0)
typedef typename Eigen::SparseMatrix<double> Eigen_sparse_matrix;
typedef typename Eigen::BiCGSTAB<Eigen_sparse_matrix, Eigen::IncompleteLUT<double> > Eigen_solver;
typedef CGAL::Eigen_solver_traits<Eigen_solver> Default_solver;
#else
typedef bool Default_solver;//compilation should crash
//if no solver is provided and Eigen version < 3.2
#endif
#else
typedef bool Default_solver;//compilation should crash
// if no solver is provided and Eigen version < 3.2
#endif
#if defined(CGAL_EIGEN3_ENABLED)
CGAL_static_assertion_msg(
(!boost::is_same<typename GetSolver<NamedParameters, Default_solver>::type, bool>::value) || EIGEN_VERSION_AT_LEAST(3, 2, 0),
"Eigen3 version 3.2 or later is required.");
#else
CGAL_static_assertion_msg(
(!boost::is_same<typename GetSolver<NamedParameters, Default_solver>::type, bool>::value),
"Eigen3 version 3.2 or later is required.");
#endif
typedef typename GetSolver<NamedParameters, Default_solver>::type Sparse_solver;
typedef typename Sparse_solver::Matrix Eigen_matrix;
typedef typename Sparse_solver::Vector Eigen_vector;
Sparse_solver solver = choose_param(get_param(np, internal_np::sparse_linear_solver), Default_solver());
std::size_t n = vertices(tmesh).size();
Eigen_matrix A(n, n);
Eigen_vector bx(n), by(n), bz(n), Xx(n), Xy(n), Xz(n);
std::vector<CGAL::Triple<std::size_t, std::size_t, double> > stiffness;
internal::Shape_smoother<TriangleMesh, VertexPointMap, VCMap, Sparse_solver, GeomTraits> smoother(tmesh, vpmap, vcmap, gt);
smoother.init_smoothing(faces);
// For robustness reasons, the laplacian coefficients are computed only once (only the mass
// matrix is updated at every iteration). See Kazdhan et al. "Can Mean-Curvature Flow Be Made Non-Singular?".
smoother.calculate_stiffness_matrix_elements(stiffness);
for(std::size_t iter=0; iter<nb_iterations; ++iter)
{
#ifdef CGAL_PMP_SMOOTHING_VERBOSE
std::cout << "iteration #" << iter << std::endl;
#endif
smoother.setup_system(A, bx, by, bz, stiffness, time);
if(smoother.solve_system(A, Xx, Xy, Xz, bx, by, bz, solver))
{
smoother.update_mesh(Xx, Xy, Xz);
}
else
{
#ifdef CGAL_PMP_SMOOTHING_VERBOSE
std::cerr << "Failed to solve system!" << std::endl;
#endif
break;
}
#ifdef CGAL_PMP_SMOOTHING_OUTPUT_INTERMEDIARY_STEPS
std::stringstream oss;
oss << "smoothed_mesh_step_" << iter << ".off" << std::ends;
std::ofstream out(oss.str().c_str());
out.precision(17);
out << tmesh;
out.close();
#endif
}
}
template<typename TriangleMesh, typename FaceRange>
void smooth_shape(const FaceRange& faces,
TriangleMesh& tmesh,
const double time)
{
smooth_shape(faces, tmesh, time, parameters::all_default());
}
template <typename TriangleMesh, typename CGAL_PMP_NP_TEMPLATE_PARAMETERS>
void smooth_shape(TriangleMesh& tmesh,
const double time,
const CGAL_PMP_NP_CLASS& np)
{
smooth_shape(faces(tmesh), tmesh, time, np);
}
template<typename TriangleMesh>
void smooth_shape(TriangleMesh& tmesh,
const double time)
{
smooth_shape(faces(tmesh), tmesh, time, parameters::all_default());
}
} // Polygon_mesh_processing
} // CGAL
#endif // CGAL_POLYGON_MESH_PROCESSING_SMOOTH_SHAPE_H

View File

@ -52,6 +52,8 @@
#include <CGAL/Polygon_mesh_processing/shape_predicates.h>
#include <CGAL/Polygon_mesh_processing/repair_polygon_soup.h>
#include <CGAL/Polygon_mesh_processing/merge_border_vertices.h>
#include <CGAL/Polygon_mesh_processing/smooth_mesh.h>
#include <CGAL/Polygon_mesh_processing/smooth_shape.h>
// the named parameter header being not documented the doc is put here for now
#ifdef DOXYGEN_RUNNING

View File

@ -50,6 +50,7 @@ if (EIGEN3_FOUND)
create_single_source_cgal_program("fairing_test.cpp")
create_single_source_cgal_program("triangulate_hole_Polyhedron_3_no_delaunay_test.cpp" )
create_single_source_cgal_program("triangulate_hole_Polyhedron_3_test.cpp")
create_single_source_cgal_program("test_shape_smoothing.cpp")
endif(EIGEN3_FOUND)
find_package( OpenMesh QUIET )
@ -92,6 +93,7 @@ endif()
create_single_source_cgal_program("triangulate_hole_polyline_test.cpp")
create_single_source_cgal_program("surface_intersection_sm_poly.cpp" )
create_single_source_cgal_program("test_orient_cc.cpp")
create_single_source_cgal_program("test_mesh_smoothing.cpp")
create_single_source_cgal_program("test_pmp_transform.cpp")
create_single_source_cgal_program("extrude_test.cpp")
create_single_source_cgal_program("test_merging_border_vertices.cpp")
@ -103,6 +105,34 @@ endif()
create_single_source_cgal_program("remove_degeneracies_test.cpp")
create_single_source_cgal_program("test_pmp_manifoldness.cpp")
set(SuiteSparse_USE_LAPACK_BLAS ON)
find_package(SuiteSparse QUIET NO_MODULE) # 1st: Try to locate the *config.cmake file.
if(NOT SuiteSparse_FOUND)
set(SuiteSparse_VERBOSE ON)
find_package(SuiteSparse QUIET) # 2nd: Use FindSuiteSparse.cmake module
if(SuiteSparse_FOUND)
include_directories(${SuiteSparse_INCLUDE_DIRS})
endif(SuiteSparse_FOUND)
else(NOT SuiteSparse_FOUND)
include(${USE_SuiteSparse})
endif(NOT SuiteSparse_FOUND)
if(SuiteSparse_FOUND)
if (SuiteSparse_UMFPACK_FOUND)
message(STATUS "SuiteSparse_LIBS: ${SuiteSparse_LIBRARIES}")
message(STATUS "Orbifold Tutte Embeddings will use UmfPackLU")
add_definitions(-DEIGEN_DONT_ALIGN_STATICALLY)
add_definitions(-DCGAL_SMP_USE_SPARSESUITE_SOLVERS)
else()
message(STATUS "NOTICE: The example `orbifold.cpp` will be compiled without the Sparsesuite library and UmfPack. Try setting SuiteSparse_UMF_INCLUDE_DIR and at least one of SuiteSparse_UMFPACK_LIBRARY_RELEASE and SuiteSparse_UMFPACK_LIBRARY_DEBUG to you UMFPACK installation.")
endif()
else(SuiteSparse_FOUND)
message(STATUS "NOTICE: The example `orbifold.cpp` will be compiled without the Sparsesuite library.")
endif(SuiteSparse_FOUND)
target_link_libraries( test_mesh_smoothing PRIVATE ceres glog ${SuiteSparse_LIBRARIES})
if( TBB_FOUND )
CGAL_target_use_TBB(test_pmp_distance)
else()

View File

@ -0,0 +1,138 @@
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/Polygon_mesh_processing/smooth_mesh.h>
#include <CGAL/property_map.h>
#include <boost/graph/graph_traits.hpp>
#include <fstream>
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point;
typedef CGAL::Surface_mesh<Point> SurfaceMesh;
typedef CGAL::Polyhedron_3<Kernel> Polyhedron;
namespace PMP = CGAL::Polygon_mesh_processing;
template <typename Mesh>
void read_mesh(const char* filename, Mesh& mesh)
{
std::ifstream input(filename);
if (!input || !(input >> mesh))
{
std::cerr << "Error: can not read file.";
std::exit(1);
}
}
template <typename Mesh>
void test_smoothing(const char* filename)
{
Mesh mesh;
read_mesh(filename, mesh);
PMP::smooth_mesh(mesh);
PMP::smooth_mesh(mesh, CGAL::parameters::number_of_iterations(10));
}
template <typename Mesh>
void test_angle_smoothing(const char* filename)
{
Mesh mesh;
read_mesh(filename, mesh);
PMP::smooth_mesh(mesh);
PMP::smooth_mesh(mesh, CGAL::parameters::number_of_iterations(10)
.use_area_smoothing(false));
}
template <typename Mesh>
void test_area_smoothing(const char* filename)
{
Mesh mesh;
read_mesh(filename, mesh);
PMP::smooth_mesh(mesh);
PMP::smooth_mesh(mesh, CGAL::parameters::number_of_iterations(10)
.use_angle_smoothing(false));
}
template <typename Mesh>
void test_angle_smoothing_without_projection(const char* filename)
{
Mesh mesh;
read_mesh(filename, mesh);
PMP::smooth_mesh(mesh, CGAL::parameters::do_project(false)
.use_area_smoothing(false));
}
template <typename Mesh>
void test_area_smoothing_without_projection(const char* filename)
{
Mesh mesh;
read_mesh(filename, mesh);
PMP::smooth_mesh(mesh, CGAL::parameters::do_project(false)
.use_angle_smoothing(false));
}
template <typename Mesh>
void test_constrained_vertices(const char* filename)
{
Mesh mesh;
read_mesh(filename, mesh);
typedef typename boost::graph_traits<Mesh>::vertex_descriptor vertex_descriptor;
typename boost::property_map<Mesh, CGAL::vertex_point_t>::type vpmap = get(CGAL::vertex_point, mesh);
std::set<vertex_descriptor> selected_vertices;
std::map<vertex_descriptor, Point> initial_positions;
for(vertex_descriptor v : vertices(mesh))
{
if(!is_border(v, mesh))
{
selected_vertices.insert(v);
initial_positions[v] = get(vpmap, v);
}
}
CGAL::Boolean_property_map<std::set<vertex_descriptor> > vcmap(selected_vertices);
PMP::smooth_mesh(mesh, CGAL::parameters::vertex_is_constrained_map(vcmap));
for(vertex_descriptor v : vertices(mesh))
{
if(!is_border(v, mesh))
assert(initial_positions.at(v) == get(vpmap, v));
}
}
int main(int /*argc*/, char** /*argv*/)
{
const char* filename_elephant = "data/elephant.off";
const char* filename_mannequin = "data/mannequin-devil.off";
// test with Surface_mesh
test_smoothing<SurfaceMesh>(filename_elephant);
test_angle_smoothing<SurfaceMesh>(filename_elephant);
test_area_smoothing<SurfaceMesh>(filename_mannequin);
test_angle_smoothing_without_projection<SurfaceMesh>(filename_elephant);
test_area_smoothing_without_projection<SurfaceMesh>(filename_mannequin);
test_constrained_vertices<SurfaceMesh>(filename_elephant);
// test with Polyhedron
test_smoothing<Polyhedron>(filename_elephant);
test_angle_smoothing<Polyhedron>(filename_elephant);
test_area_smoothing<Polyhedron>(filename_mannequin);
test_angle_smoothing_without_projection<Polyhedron>(filename_elephant);
test_area_smoothing_without_projection<Polyhedron>(filename_mannequin);
test_constrained_vertices<Polyhedron>(filename_mannequin);
return EXIT_SUCCESS;
}

View File

@ -0,0 +1,205 @@
#define CGAL_PMP_SMOOTHING_VERBOSE
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/Polyhedron_items_with_id_3.h>
#include <CGAL/Polygon_mesh_processing/smooth_shape.h>
#include <CGAL/utility.h>
#include <fstream>
#include <iostream>
#include <set>
#include <vector>
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point;
typedef CGAL::Surface_mesh<Point> SurfaceMesh;
typedef CGAL::Polyhedron_3<Kernel, CGAL::Polyhedron_items_with_id_3> Mesh_with_id;
namespace PMP = CGAL::Polygon_mesh_processing;
bool equal_doubles(double d1, double d2, double e)
{
return (d1 > d2 - e) && (d1 < d2 + e);
}
template <typename Mesh>
void test_implicit_constrained_devil(Mesh mesh)
{
#ifdef CGAL_PMP_SMOOTHING_VERBOSE
std::cout << "-- test_implicit_constrained_devil --" << std::endl;
#endif
typedef typename boost::graph_traits<Mesh>::vertex_descriptor vertex_descriptor;
typename boost::property_map<Mesh, CGAL::vertex_point_t>::type vpmap = get(CGAL::vertex_point, mesh);
// z max is 20 in the devil
std::set<vertex_descriptor> selected_vertices;
for(vertex_descriptor v : vertices(mesh))
{
const double z = get(vpmap, v).z();
if(is_border(v, mesh) || z > 19.0)
selected_vertices.insert(v);
}
CGAL::Boolean_property_map<std::set<vertex_descriptor> > vcmap(selected_vertices);
std::vector<Point> fixed_points(selected_vertices.size());
std::size_t i = 0;
for(vertex_descriptor v : selected_vertices)
fixed_points[i++] = get(vpmap, v);
const double time_step = 1.0;
PMP::smooth_shape(mesh, time_step, CGAL::parameters::vertex_is_constrained_map(vcmap)
.number_of_iterations(2));
i = 0;
for(vertex_descriptor v : selected_vertices)
{
const Point p = get(vpmap, v);
CGAL_assertion(equal_doubles(p.x(), fixed_points[i].x(), 1e-10));
CGAL_assertion(equal_doubles(p.y(), fixed_points[i].y(), 1e-10));
CGAL_assertion(equal_doubles(p.z(), fixed_points[i].z(), 1e-10));
++i;
}
#ifdef CGAL_PMP_SMOOTHING_VERBOSE
std::ofstream out("output_implicit_constrained_devil.off");
out << mesh;
out.close();
#endif
}
template <typename Mesh>
void test_implicit_constrained_elephant(Mesh mesh)
{
#ifdef CGAL_PMP_SMOOTHING_VERBOSE
std::cout << "-- test_implicit_constrained_elephant --" << std::endl;
#endif
typedef typename boost::graph_traits<Mesh>::vertex_descriptor vertex_descriptor;
typename boost::property_map<Mesh, CGAL::vertex_point_t>::type vpmap = get(CGAL::vertex_point, mesh);
std::set<vertex_descriptor> selected_vertices;
for(vertex_descriptor v : vertices(mesh))
{
const double y = get(vpmap, v).z();
if(y > 0.14)
selected_vertices.insert(v);
}
CGAL::Boolean_property_map<std::set<vertex_descriptor> > vcmap(selected_vertices);
Point fixed_point = get(vpmap, *selected_vertices.begin());
const double time_step = 1.0;
PMP::smooth_shape(mesh, time_step,
CGAL::parameters::vertex_is_constrained_map(vcmap)
.number_of_iterations(1));
CGAL_assertion(equal_doubles(get(vpmap, *selected_vertices.begin()).x(), fixed_point.x(), 1e-14));
CGAL_assertion(equal_doubles(get(vpmap, *selected_vertices.begin()).y(), fixed_point.y(), 1e-14));
CGAL_assertion(equal_doubles(get(vpmap, *selected_vertices.begin()).z(), fixed_point.z(), 1e-14));
#ifdef CGAL_PMP_SMOOTHING_VERBOSE
std::ofstream out("output_implicit_constrained_elephant.off");
out << mesh;
out.close();
#endif
}
template <typename Mesh>
void test_curvature_flow_time_step(Mesh mesh)
{
#ifdef CGAL_PMP_SMOOTHING_VERBOSE
std::cout << "-- test_curvature_flow_time_step --" << std::endl;
#endif
const double time_step = 1e-15;
PMP::smooth_shape(mesh, time_step);
#ifdef CGAL_PMP_SMOOTHING_VERBOSE
std::ofstream out("output_devil_time_step.off");
out << mesh;
out.close();
#endif
}
template <typename Mesh>
void test_curvature_flow(Mesh mesh)
{
#ifdef CGAL_PMP_SMOOTHING_VERBOSE
std::cout << "-- test_curvature_flow --" << std::endl;
#endif
const double time_step = 1.0;
PMP::smooth_shape(mesh, time_step);
#ifdef CGAL_PMP_SMOOTHING_VERBOSE
std::ofstream out("output_precision_elephant.off");
out << mesh;
out.close();
#endif
}
int main(int, char**)
{
const char* filename_devil = "data/mannequin-devil.off";
const char* filename_elephant = "data/elephant.off";
std::ifstream input1(filename_devil);
SurfaceMesh mesh_devil;
if(!input1 || !(input1 >> mesh_devil))
{
std::cerr << "Error: can not read file.";
return 1;
}
input1.close();
std::ifstream input2(filename_elephant);
SurfaceMesh mesh_elephant;
if(!input2 || !(input2 >> mesh_elephant))
{
std::cerr << "Error: can not read file.";
return 1;
}
input2.close();
test_curvature_flow_time_step<SurfaceMesh>(mesh_devil);
test_curvature_flow<SurfaceMesh>(mesh_elephant);
test_implicit_constrained_elephant<SurfaceMesh>(mesh_elephant);
test_implicit_constrained_devil<SurfaceMesh>(mesh_devil);
input1.open(filename_devil);
Mesh_with_id pl_mesh_devil;
if(!input1 || !(input1 >> pl_mesh_devil)){
std::cerr << "Error: can not read file.";
return EXIT_FAILURE;
}
input1.close();
// Polyhedron
input2.open(filename_elephant);
Mesh_with_id pl_mesh_elephant;
if(!input2 || !(input2 >> pl_mesh_elephant))
{
std::cerr << "Error: can not read file.";
return EXIT_FAILURE;
}
input2.close();
set_halfedgeds_items_id(pl_mesh_devil);
set_halfedgeds_items_id(pl_mesh_elephant);
test_curvature_flow_time_step<Mesh_with_id>(pl_mesh_devil);
test_curvature_flow<Mesh_with_id>(pl_mesh_elephant);
test_implicit_constrained_elephant<Mesh_with_id>(pl_mesh_elephant);
test_implicit_constrained_devil<Mesh_with_id>(pl_mesh_devil);
return EXIT_SUCCESS;
}

View File

@ -115,3 +115,32 @@ target_link_libraries(engrave_text_plugin PUBLIC scene_surface_mesh_item scene_s
polyhedron_demo_plugin(extrude_plugin Extrude_plugin)
target_link_libraries(extrude_plugin PUBLIC scene_surface_mesh_item scene_selection_item)
set(SuiteSparse_USE_LAPACK_BLAS ON)
find_package(SuiteSparse QUIET NO_MODULE) # 1st: Try to locate the *config.cmake file.
if(NOT SuiteSparse_FOUND)
set(SuiteSparse_VERBOSE ON)
find_package(SuiteSparse QUIET) # 2nd: Use FindSuiteSparse.cmake module
if(SuiteSparse_FOUND)
include_directories(${SuiteSparse_INCLUDE_DIRS})
endif(SuiteSparse_FOUND)
else(NOT SuiteSparse_FOUND)
include(${USE_SuiteSparse})
endif(NOT SuiteSparse_FOUND)
if(SuiteSparse_FOUND)
if (SuiteSparse_UMFPACK_FOUND)
message(STATUS "SuiteSparse_LIBS: ${SuiteSparse_LIBRARIES}")
message(STATUS "Orbifold Tutte Embeddings will use UmfPackLU")
add_definitions(-DEIGEN_DONT_ALIGN_STATICALLY)
add_definitions(-DCGAL_SMP_USE_SPARSESUITE_SOLVERS)
else()
message(STATUS "NOTICE: The example `orbifold.cpp` will be compiled without the Sparsesuite library and UmfPack. Try setting SuiteSparse_UMF_INCLUDE_DIR and at least one of SuiteSparse_UMFPACK_LIBRARY_RELEASE and SuiteSparse_UMFPACK_LIBRARY_DEBUG to you UMFPACK installation.")
endif()
else(SuiteSparse_FOUND)
message(STATUS "NOTICE: The example `orbifold.cpp` will be compiled without the Sparsesuite library.")
endif(SuiteSparse_FOUND)
qt5_wrap_ui( smoothingUI_FILES Smoothing_plugin.ui)
polyhedron_demo_plugin(smoothing_plugin Smoothing_plugin ${smoothingUI_FILES})
target_link_libraries(smoothing_plugin PUBLIC ceres glog scene_surface_mesh_item scene_selection_item ${SuiteSparse_LIBRARIES})

View File

@ -0,0 +1,295 @@
#include <QtCore/qglobal.h>
#include <CGAL/Three/Polyhedron_demo_plugin_interface.h>
#include <CGAL/Three/Polyhedron_demo_plugin_helper.h>
#include <CGAL/iterator.h>
#include <CGAL/utility.h>
#include <CGAL/property_map.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/smooth_mesh.h>
#include <CGAL/Polygon_mesh_processing/smooth_shape.h>
#include "Scene.h"
#include "Scene_surface_mesh_item.h"
#include "Scene_polyhedron_selection_item.h"
#include "ui_Smoothing_plugin.h"
#include <QTime>
#include <QAction>
#include <QMainWindow>
#include <QApplication>
#include <QDockWidget>
#include <QString>
#include <QInputDialog>
#include <QtPlugin>
#include <QMessageBox>
using namespace CGAL::Polygon_mesh_processing;
using namespace CGAL::Three;
typedef Scene_surface_mesh_item Scene_face_graph_item;
typedef Scene_face_graph_item::Face_graph Face_graph;
typedef boost::graph_traits<Face_graph>::vertex_descriptor vertex_descriptor;
typedef boost::graph_traits<Face_graph>::edge_descriptor edge_descriptor;
typedef boost::graph_traits<Face_graph>::face_descriptor face_descriptor;
typedef CGAL::dynamic_vertex_property_t<bool> Vertex_bool_property;
typedef typename boost::property_map<Face_graph, Vertex_bool_property>::type VCMap;
class Polyhedron_demo_smothing_plugin
: public QObject, public Polyhedron_demo_plugin_helper
{
Q_OBJECT
Q_INTERFACES(CGAL::Three::Polyhedron_demo_plugin_interface)
Q_PLUGIN_METADATA(IID "com.geometryfactory.PolyhedronDemo.PluginInterface/1.0")
public:
void init(QMainWindow* mainWindow, Scene_interface* scene_interface, Messages_interface*)
{
scene = scene_interface;
mw = mainWindow;
actionSmoothing_ = new QAction(tr("Mesh and Shape Smoothing"), mw);
actionSmoothing_->setProperty("subMenuName", "Polygon Mesh Processing");
connect(actionSmoothing_, SIGNAL(triggered()), this, SLOT(smoothing_action()));
dock_widget = new QDockWidget("Smoothing", mw);
dock_widget->setVisible(false);
ui_widget.setupUi(dock_widget);
addDockWidget(dock_widget);
connect(ui_widget.area_smoothing_checkBox, SIGNAL(toggled(bool)),
ui_widget.flip_checkBox, SLOT(setEnabled(bool)));
connect(ui_widget.mesh_smoothing_button, SIGNAL(clicked()), this, SLOT(on_mesh_smoothing_clicked()));
connect(ui_widget.shape_smoothing_button, SIGNAL(clicked()), this, SLOT(on_shape_smoothing_clicked()));
}
QList<QAction*> actions() const
{
return QList<QAction*>() << actionSmoothing_;
}
bool applicable(QAction*) const
{
const Scene_interface::Item_id index = scene->mainSelectionIndex();
if(qobject_cast<Scene_face_graph_item*>(scene->item(index)))
return true;
else if(qobject_cast<Scene_polyhedron_selection_item*>(scene->item(index)))
return true;
else
return false;
}
virtual void closure()
{
dock_widget->hide();
}
void init_ui()
{
ui_widget.time_step_spinBox->setValue(0.00001);
ui_widget.time_step_spinBox->setSingleStep(0.00001);
ui_widget.time_step_spinBox->setMinimum(1e-6);
ui_widget.smooth_iter_spinBox->setValue(1);
ui_widget.projection_checkBox->setChecked(true);
ui_widget.area_smoothing_checkBox->setChecked(false);
ui_widget.flip_checkBox->setDisabled(true);
}
void mark_border_vertices(const VCMap vcmap, const Face_graph& pmesh) const
{
for(vertex_descriptor v : vertices(pmesh))
put(vcmap, v, false);
for(halfedge_descriptor h : halfedges(pmesh))
{
if(CGAL::is_border(h, pmesh))
put(vcmap, target(h, pmesh), true);
}
}
void mark_selected_vertices(const VCMap vcmap,
const Face_graph& pmesh,
Scene_polyhedron_selection_item* selection_item) const
{
for(vertex_descriptor v : selection_item->selected_vertices)
put(vcmap, v, true);
for(edge_descriptor e : selection_item->selected_edges)
{
put(vcmap, source(e, pmesh), true);
put(vcmap, target(e, pmesh), true);
}
}
public Q_SLOTS:
void smoothing_action()
{
dock_widget->show();
dock_widget->raise();
const Scene_interface::Item_id index = scene->mainSelectionIndex();
Scene_face_graph_item* poly_item = qobject_cast<Scene_face_graph_item*>(scene->item(index));
Scene_polyhedron_selection_item* selection_item =
qobject_cast<Scene_polyhedron_selection_item*>(scene->item(index));
if(poly_item || selection_item)
{
init_ui();
}
}
void on_mesh_smoothing_clicked()
{
const Scene_interface::Item_id index = scene->mainSelectionIndex();
Scene_face_graph_item* poly_item = qobject_cast<Scene_face_graph_item*>(scene->item(index));
Scene_polyhedron_selection_item* selection_item = qobject_cast<Scene_polyhedron_selection_item*>(scene->item(index));
if(!poly_item && !selection_item)
return;
Face_graph& pmesh = (poly_item != nullptr) ? * poly_item->polyhedron() : * selection_item->polyhedron();
const unsigned int nb_iter = ui_widget.smooth_iter_spinBox->value();
const bool projection = ui_widget.projection_checkBox->isChecked();
const bool use_safety_measures = ui_widget.sanity_checkBox->isChecked();
const bool constrain_border_vertices = ui_widget.border_button->isChecked() && !CGAL::is_closed(pmesh);
const bool use_angle_smoothing = ui_widget.angle_smoothing_checkBox->isChecked();
const bool use_area_smoothing = ui_widget.area_smoothing_checkBox->isChecked();
const bool use_Delaunay_flips = ui_widget.flip_checkBox->isChecked();
QApplication::setOverrideCursor(Qt::WaitCursor);
VCMap vcmap = get(Vertex_bool_property(), pmesh);
for(vertex_descriptor v : vertices(pmesh))
put(vcmap, v, false);
if(constrain_border_vertices)
mark_border_vertices(vcmap, pmesh);
if(poly_item)
{
smooth_mesh(pmesh, parameters::do_project(projection)
.number_of_iterations(nb_iter)
.vertex_is_constrained_map(vcmap)
.use_safety_constraints(use_safety_measures)
.use_angle_smoothing(use_angle_smoothing)
.use_area_smoothing(use_area_smoothing)
.use_Delaunay_flips(use_Delaunay_flips));
poly_item->invalidateOpenGLBuffers();
poly_item->itemChanged();
}
else if(selection_item)
{
mark_selected_vertices(vcmap, pmesh, selection_item);
// No faces selected --> use all faces
if(std::begin(selection_item->selected_facets) == std::end(selection_item->selected_facets))
{
smooth_mesh(pmesh, parameters::do_project(projection)
.number_of_iterations(nb_iter)
.vertex_is_constrained_map(vcmap)
.edge_is_constrained_map(selection_item->constrained_edges_pmap())
.use_safety_constraints(use_safety_measures)
.use_angle_smoothing(use_angle_smoothing)
.use_area_smoothing(use_area_smoothing)
.use_Delaunay_flips(use_Delaunay_flips));
}
else // some faces exist in the selection
{
smooth_mesh(selection_item->selected_facets, pmesh, parameters::do_project(projection)
.number_of_iterations(nb_iter)
.vertex_is_constrained_map(vcmap)
.edge_is_constrained_map(selection_item->constrained_edges_pmap())
.use_safety_constraints(use_safety_measures)
.use_angle_smoothing(use_angle_smoothing)
.use_area_smoothing(use_area_smoothing)
.use_Delaunay_flips(use_Delaunay_flips));
}
selection_item->poly_item_changed();
selection_item->changed_with_poly_item();
}
QApplication::restoreOverrideCursor();
}
void on_shape_smoothing_clicked()
{
const Scene_interface::Item_id index = scene->mainSelectionIndex();
Scene_face_graph_item* poly_item = qobject_cast<Scene_face_graph_item*>(scene->item(index));
Scene_polyhedron_selection_item* selection_item = qobject_cast<Scene_polyhedron_selection_item*>(scene->item(index));
if(!poly_item && !selection_item)
return;
Face_graph& pmesh = (poly_item != nullptr) ? * poly_item->polyhedron() : * selection_item->polyhedron();
const double time_step = ui_widget.time_step_spinBox->value();
const unsigned int nb_iter = ui_widget.smooth_iter_spinBox->value();
const bool constrain_border_vertices = ui_widget.border_button->isChecked() && !CGAL::is_closed(pmesh);
QApplication::setOverrideCursor(Qt::WaitCursor);
VCMap vcmap = get(Vertex_bool_property(), pmesh);
for(vertex_descriptor v : vertices(pmesh))
put(vcmap, v, false);
if(constrain_border_vertices)
mark_border_vertices(vcmap, pmesh);
if(poly_item)
{
smooth_shape(pmesh, time_step, parameters::number_of_iterations(nb_iter)
.vertex_is_constrained_map(vcmap));
poly_item->invalidateOpenGLBuffers();
poly_item->itemChanged();
}
else if(selection_item)
{
mark_selected_vertices(vcmap, pmesh, selection_item);
if(std::begin(selection_item->selected_facets) == std::end(selection_item->selected_facets))
{
smooth_shape(pmesh, time_step, parameters::number_of_iterations(nb_iter)
.vertex_is_constrained_map(vcmap));
}
else
{
smooth_shape(selection_item->selected_facets, pmesh, time_step,
parameters::number_of_iterations(nb_iter)
.vertex_is_constrained_map(vcmap));
}
selection_item->poly_item_changed();
selection_item->changed_with_poly_item();
}
else
{
std::cerr << "Something's gone wrong.\n";
CGAL_assertion(false);
}
QApplication::restoreOverrideCursor();
}
private:
QAction* actionSmoothing_;
QDockWidget* dock_widget;
Ui::Smoothing ui_widget;
};
#include "Smoothing_plugin.moc"

View File

@ -0,0 +1,178 @@
<?xml version="1.0" encoding="UTF-8"?>
<ui version="4.0">
<class>Smoothing</class>
<widget class="QDockWidget" name="Smoothing">
<property name="geometry">
<rect>
<x>0</x>
<y>0</y>
<width>518</width>
<height>478</height>
</rect>
</property>
<property name="windowTitle">
<string>Smoothing</string>
</property>
<widget class="QWidget" name="dockWidgetContents">
<layout class="QGridLayout" name="gridLayout_3">
<item row="0" column="1">
<widget class="QSpinBox" name="smooth_iter_spinBox">
<property name="minimum">
<number>1</number>
</property>
<property name="maximum">
<number>1000</number>
</property>
</widget>
</item>
<item row="9" column="0" colspan="2">
<widget class="QGroupBox" name="Shape_groupBox">
<property name="title">
<string>Shape Smoothing</string>
</property>
<layout class="QGridLayout" name="gridLayout">
<item row="1" column="2">
<widget class="QDoubleSpinBox" name="time_step_spinBox">
<property name="decimals">
<number>6</number>
</property>
<property name="minimum">
<double>0.000001000000000</double>
</property>
</widget>
</item>
<item row="6" column="2">
<widget class="QPushButton" name="shape_smoothing_button">
<property name="text">
<string>Smooth Shape</string>
</property>
</widget>
</item>
<item row="1" column="0">
<widget class="QLabel" name="time_step_label">
<property name="toolTip">
<string>The time step should not be too large or the smoothing will be unsatisfactory</string>
</property>
<property name="text">
<string>Time Step:</string>
</property>
</widget>
</item>
</layout>
</widget>
</item>
<item row="2" column="0">
<widget class="QCheckBox" name="border_button">
<property name="toolTip">
<string>If you are using a selection item, this will combine your selected vertices and the border vertices</string>
</property>
<property name="text">
<string>Constrain Border Vertices</string>
</property>
</widget>
</item>
<item row="4" column="0" colspan="2">
<spacer name="verticalSpacer">
<property name="orientation">
<enum>Qt::Vertical</enum>
</property>
<property name="sizeHint" stdset="0">
<size>
<width>20</width>
<height>40</height>
</size>
</property>
</spacer>
</item>
<item row="0" column="0">
<widget class="QLabel" name="label">
<property name="text">
<string>Number of iterations</string>
</property>
</widget>
</item>
<item row="5" column="0" colspan="2">
<widget class="QGroupBox" name="Mesh_groupBox">
<property name="title">
<string>Mesh Smoothing</string>
</property>
<layout class="QGridLayout" name="gridLayout_2">
<item row="4" column="1">
<widget class="QPushButton" name="mesh_smoothing_button">
<property name="enabled">
<bool>true</bool>
</property>
<property name="text">
<string>Smooth Mesh</string>
</property>
<property name="checkable">
<bool>false</bool>
</property>
</widget>
</item>
<item row="2" column="0">
<widget class="QCheckBox" name="sanity_checkBox">
<property name="toolTip">
<string>By default, a serie of checks are performed when a displacement is computed (e.g. is the minimul angle around the vertex not reduced?). However, sometimes local bad moves are necessary to achieve a better result overall.</string>
</property>
<property name="text">
<string>Do not apply worsening moves</string>
</property>
</widget>
</item>
<item row="5" column="1">
<spacer name="verticalSpacer_2">
<property name="orientation">
<enum>Qt::Vertical</enum>
</property>
<property name="sizeHint" stdset="0">
<size>
<width>20</width>
<height>40</height>
</size>
</property>
</spacer>
</item>
<item row="0" column="0">
<widget class="QCheckBox" name="area_smoothing_checkBox">
<property name="text">
<string>Use Area Smoothing</string>
</property>
</widget>
</item>
<item row="1" column="0">
<widget class="QCheckBox" name="angle_smoothing_checkBox">
<property name="text">
<string>Use Angle Smoothing</string>
</property>
</widget>
</item>
<item row="3" column="0">
<widget class="QCheckBox" name="projection_checkBox">
<property name="toolTip">
<string>Reprojection on the initial surface will not happen if there are self-intersections in the mesh, except if worsening moves are allowed (the option just above).</string>
</property>
<property name="text">
<string>Re-project Vertices</string>
</property>
</widget>
</item>
<item row="0" column="1">
<widget class="QCheckBox" name="flip_checkBox">
<property name="text">
<string>Allow Delaunay flips</string>
</property>
<property name="checkable">
<bool>true</bool>
</property>
</widget>
</item>
</layout>
</widget>
</item>
</layout>
</widget>
</widget>
<resources/>
<connections/>
</ui>

View File

@ -214,6 +214,7 @@ triangulate_facets_plugin \
trivial_plugin \
vtk_plugin \
xyz_plugin \
smoothing_plugin \
all
do
if ${MAKE_CMD} -f Makefile help | grep "$target" > /dev/null; then