1 // Copyright (c) 2018 GeometryFactory (France).
2 // All rights reserved.
3 //
4 // This file is part of CGAL (www.cgal.org).
5 //
6 // $URL: https://github.com/CGAL/cgal/blob/v5.3/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/smooth_mesh.h $
7 // $Id: smooth_mesh.h f55ef7d 2020-10-09T18:36:17+02:00 Mael Rouxel-Labbé
8 // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
9 //
10 //
11 // Author(s) : Mael Rouxel-Labbé
12 // Konstantinos Katrioplas (konst.katrioplas@gmail.com)
13
14 #ifndef CGAL_POLYGON_MESH_PROCESSING_SMOOTH_MESH_H
15 #define CGAL_POLYGON_MESH_PROCESSING_SMOOTH_MESH_H
16
17 #include <CGAL/license/Polygon_mesh_processing/meshing_hole_filling.h>
18
19 #include <CGAL/Polygon_mesh_processing/internal/Smoothing/mesh_smoothing_impl.h>
20 #include <CGAL/Polygon_mesh_processing/internal/Smoothing/smoothing_evaluation.h>
21 #include <CGAL/Polygon_mesh_processing/self_intersections.h>
22
23 #include <CGAL/boost/graph/Named_function_parameters.h>
24 #include <CGAL/boost/graph/named_params_helper.h>
25
26 #include <CGAL/property_map.h>
27
28 #ifdef DOXYGEN_RUNNING
29 #define CGAL_PMP_NP_TEMPLATE_PARAMETERS NamedParameters
30 #define CGAL_PMP_NP_CLASS NamedParameters
31 #endif
32
33 namespace CGAL {
34 namespace Polygon_mesh_processing {
35
36 /*!
37 * \ingroup PMP_meshing_grp
38 *
39 * \short smooths a triangulated region of a polygon mesh.
40 *
41 * This function attempts to make the triangle angle and area distributions as uniform as possible
42 * by moving (non-constrained) vertices.
43 *
44 * Angle-based smoothing does not change the combinatorial information of the mesh. Area-based smoothing
45 * might change the combinatorial information, unless specified otherwise. It is also possible
46 * to make the smoothing algorithm "safer" by rejecting moves that, when applied, would worsen the
47 * quality of the mesh, e.g. that would decrease the value of the smallest angle around a vertex or
48 * create self-intersections.
49 *
50 * Optionally, the points are reprojected after each iteration.
51 *
52 * @tparam TriangleMesh model of `MutableFaceGraph`.
53 * @tparam FaceRange range of `boost::graph_traits<TriangleMesh>::%face_descriptor`,
54 model of `Range`. Its iterator type is `ForwardIterator`.
55 * @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
56 *
57 * @param tmesh a polygon mesh with triangulated surface patches to be smoothed.
58 * @param faces the range of triangular faces defining one or several surface patches to be smoothed.
59 * @param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
60 *
61 * \cgalNamedParamsBegin
62 * \cgalParamNBegin{number_of_iterations}
63 * \cgalParamDescription{the number of iterations for the sequence of the smoothing iterations performed}
64 * \cgalParamType{unsigned int}
65 * \cgalParamDefault{`1`}
66 * \cgalParamNEnd
67 *
68 * \cgalParamNBegin{use_angle_smoothing}
69 * \cgalParamDescription{value to indicate whether angle-based smoothing should be used}
70 * \cgalParamType{Boolean}
71 * \cgalParamDefault{`true`}
72 * \cgalParamNEnd
73 *
74 * \cgalParamNBegin{use_area_smoothing}
75 * \cgalParamDescription{value to indicate whether area-based smoothing should be used}
76 * \cgalParamType{Boolean}
77 * \cgalParamDefault{`true`}
78 * \cgalParamNEnd
79 *
80 * \cgalParamNBegin{vertex_point_map}
81 * \cgalParamDescription{a property map associating points to the vertices of `tmesh`}
82 * \cgalParamType{a class model of `ReadWritePropertyMap` with `boost::graph_traits<TriangleMesh>::%vertex_descriptor`
83 * as key type and `%Point_3` as value type}
84 * \cgalParamDefault{`boost::get(CGAL::vertex_point, tmesh)`}
85 * \cgalParamExtra{If this parameter is omitted, an internal property map for `CGAL::vertex_point_t`
86 * must be available in `TriangleMesh`.}
87 * \cgalParamNEnd
88 *
89 * \cgalParamNBegin{geom_traits}
90 * \cgalParamDescription{an instance of a geometric traits class}
91 * \cgalParamType{a class model of `Kernel`}
92 * \cgalParamDefault{a \cgal Kernel deduced from the point type, using `CGAL::Kernel_traits`}
93 * \cgalParamExtra{The geometric traits class must be compatible with the vertex point type.}
94 * \cgalParamNEnd
95 *
96 * \cgalParamNBegin{use_safety_constraints}
97 * \cgalParamDescription{If `true`, vertex moves that would worsen the mesh are ignored.}
98 * \cgalParamType{Boolean}
99 * \cgalParamDefault{`false`}
100 * \cgalParamNEnd
101 *
102 * \cgalParamNBegin{use_Delaunay_flips}
103 * \cgalParamDescription{If `true`, area-based smoothing will be completed by a phase of
104 * Delaunay-based edge-flips to prevent the creation of elongated triangles.}
105 * \cgalParamType{Boolean}
106 * \cgalParamDefault{`true`}
107 * \cgalParamNEnd
108 *
109 * \cgalParamNBegin{do_project}
110 * \cgalParamDescription{If `true`, points are projected onto the initial surface after each iteration.}
111 * \cgalParamType{Boolean}
112 * \cgalParamDefault{`true`}
113 * \cgalParamNEnd
114 *
115 * \cgalParamNBegin{vertex_is_constrained_map}
116 * \cgalParamDescription{a property map containing the constrained-or-not status of each vertex of `tmesh`.}
117 * \cgalParamType{a class model of `ReadWritePropertyMap` with `boost::graph_traits<TriangleMesh>::%vertex_descriptor`
118 * as key type and `bool` as value type. It must be default constructible.}
119 * \cgalParamDefault{a default property map where no vertex is constrained}
120 * \cgalParamExtra{A constrained vertex cannot be modified at all during smoothing.}
121 * \cgalParamNEnd
122 *
123 * \cgalParamNBegin{edge_is_constrained_map}
124 * \cgalParamDescription{a property map containing the constrained-or-not status of each edge of `tmesh`.}
125 * \cgalParamType{a class model of `ReadWritePropertyMap` with `boost::graph_traits<TriangleMesh>::%edge_descriptor`
126 * as key type and `bool` as value type. It must be default constructible.}
127 * \cgalParamDefault{a default property map where no edge is constrained}
128 * \cgalParamExtra{A constrained edge cannot be modified at all during smoothing.}
129 * \cgalParamNEnd
130 * \cgalNamedParamsEnd
131 *
132 * @warning The third party library \link thirdpartyCeres Ceres \endlink is required
133 * to use area-based smoothing.
134 *
135 * @pre `tmesh` does not contain any degenerate faces
136 */
137 template<typename TriangleMesh, typename FaceRange, typename NamedParameters>
smooth_mesh(const FaceRange & faces,TriangleMesh & tmesh,const NamedParameters & np)138 void smooth_mesh(const FaceRange& faces,
139 TriangleMesh& tmesh,
140 const NamedParameters& np)
141 {
142 typedef typename boost::graph_traits<TriangleMesh>::vertex_descriptor vertex_descriptor;
143 typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor halfedge_descriptor;
144 typedef typename boost::graph_traits<TriangleMesh>::edge_descriptor edge_descriptor;
145 typedef typename boost::graph_traits<TriangleMesh>::face_descriptor face_descriptor;
146
147 typedef typename GetGeomTraits<TriangleMesh, NamedParameters>::type GeomTraits;
148 typedef typename GetVertexPointMap<TriangleMesh, NamedParameters>::type VertexPointMap;
149
150 // We need a default pmap that is not just 'constant_pmap(false)' because if an edge is constrained,
151 // its vertices are constrained.
152 typedef CGAL::dynamic_vertex_property_t<bool> Vertex_property_tag;
153 typedef typename boost::property_map<TriangleMesh, Vertex_property_tag>::type Default_VCMap;
154 typedef typename internal_np::Lookup_named_param_def<internal_np::vertex_is_constrained_t,
155 NamedParameters,
156 Default_VCMap
157 > ::type VCMap;
158
159 typedef typename internal_np::Lookup_named_param_def<internal_np::edge_is_constrained_t,
160 NamedParameters,
161 Static_boolean_property_map<edge_descriptor, false> // default
162 > ::type ECMap;
163
164 typedef internal::Area_smoother<TriangleMesh, VertexPointMap, GeomTraits> Area_optimizer;
165 typedef internal::Mesh_smoother<Area_optimizer, TriangleMesh,
166 VertexPointMap, VCMap, GeomTraits> Area_smoother;
167 typedef internal::Delaunay_edge_flipper<TriangleMesh, VertexPointMap,
168 ECMap, GeomTraits> Delaunay_flipper;
169
170 typedef internal::Angle_smoother<TriangleMesh, VertexPointMap, GeomTraits> Angle_optimizer;
171 typedef internal::Mesh_smoother<Angle_optimizer, TriangleMesh,
172 VertexPointMap, VCMap, GeomTraits> Angle_smoother;
173
174 typedef typename GeomTraits::Triangle_3 Triangle;
175 typedef std::vector<Triangle> Triangle_container;
176
177 typedef CGAL::AABB_triangle_primitive<GeomTraits,
178 typename Triangle_container::iterator> AABB_Primitive;
179 typedef CGAL::AABB_traits<GeomTraits, AABB_Primitive> AABB_Traits;
180 typedef CGAL::AABB_tree<AABB_Traits> Tree;
181
182 if(std::begin(faces) == std::end(faces))
183 return;
184
185 using parameters::choose_parameter;
186 using parameters::get_parameter;
187
188 // named parameters
189 GeomTraits gt = choose_parameter<GeomTraits>(get_parameter(np, internal_np::geom_traits));
190 VertexPointMap vpmap = choose_parameter(get_parameter(np, internal_np::vertex_point),
191 get_property_map(CGAL::vertex_point, tmesh));
192
193 const bool use_angle_smoothing = choose_parameter(get_parameter(np, internal_np::use_angle_smoothing), true);
194 bool use_area_smoothing = choose_parameter(get_parameter(np, internal_np::use_area_smoothing), true);
195
196 #ifndef CGAL_PMP_USE_CERES_SOLVER
197 std::cerr << "Area-based smoothing requires the Ceres Library, which is not available." << std::endl;
198 std::cerr << "No such smoothing will be performed!" << std::endl;
199 use_area_smoothing = false;
200 #endif
201
202 if(!use_angle_smoothing && !use_area_smoothing)
203 std::cerr << "Called PMP::smooth_mesh() without any smoothing method selected or available" << std::endl;
204
205 unsigned int nb_iterations = choose_parameter(get_parameter(np, internal_np::number_of_iterations), 1);
206 const bool do_project = choose_parameter(get_parameter(np, internal_np::do_project), true);
207 const bool use_safety_constraints = choose_parameter(get_parameter(np, internal_np::use_safety_constraints), true);
208 const bool use_Delaunay_flips = choose_parameter(get_parameter(np, internal_np::use_Delaunay_flips), true);
209
210 VCMap vcmap = choose_parameter(get_parameter(np, internal_np::vertex_is_constrained),
211 get(Vertex_property_tag(), tmesh));
212
213 // If it's the default vcmap, manually set everything to false because the dynamic pmap has no default initialization
214 if((std::is_same<VCMap, Default_VCMap>::value))
215 {
216 for(vertex_descriptor v : vertices(tmesh))
217 put(vcmap, v, false);
218 }
219
220 ECMap ecmap = choose_parameter(get_parameter(np, internal_np::edge_is_constrained),
221 Static_boolean_property_map<edge_descriptor, false>());
222
223 // a constrained edge has constrained extremities
224 for(face_descriptor f : faces)
225 {
226 if(f == boost::graph_traits<TriangleMesh>::null_face())
227 continue;
228
229 for(halfedge_descriptor h : CGAL::halfedges_around_face(halfedge(f, tmesh), tmesh))
230 {
231 if(get(ecmap, edge(h, tmesh)))
232 {
233 put(vcmap, source(h, tmesh), true);
234 put(vcmap, target(h, tmesh), true);
235 }
236 }
237 }
238
239 // Construct the AABB tree (if needed for reprojection)
240 std::vector<Triangle> input_triangles;
241
242 if(do_project)
243 {
244 input_triangles.reserve(faces.size());
245
246 for(face_descriptor f : faces)
247 {
248 halfedge_descriptor h = halfedge(f, tmesh);
249 if(is_border(h, tmesh)) // should not happen, but just in case
250 continue;
251
252 input_triangles.push_back(gt.construct_triangle_3_object()(get(vpmap, source(h, tmesh)),
253 get(vpmap, target(h, tmesh)),
254 get(vpmap, target(next(h, tmesh), tmesh))));
255 }
256 }
257
258 Tree aabb_tree(input_triangles.begin(), input_triangles.end());
259
260 // Setup the working ranges and check some preconditions
261 Angle_smoother angle_smoother(tmesh, vpmap, vcmap, gt);
262 Area_smoother area_smoother(tmesh, vpmap, vcmap, gt);
263 Delaunay_flipper delaunay_flipper(tmesh, vpmap, ecmap, gt);
264
265 if(use_angle_smoothing)
266 angle_smoother.init_smoothing(faces);
267
268 if(use_area_smoothing)
269 area_smoother.init_smoothing(faces);
270
271 for(unsigned int i=0; i<nb_iterations; ++i)
272 {
273 #ifdef CGAL_PMP_SMOOTHING_DEBUG
274 std::cout << "Iteration #" << i << std::endl;
275 #endif
276
277 if(use_area_smoothing)
278 {
279 #ifdef CGAL_PMP_SMOOTHING_DEBUG
280 std::cout << "Smooth areas..." << std::endl;
281 #endif
282
283 // First apply area smoothing...
284 area_smoother.optimize(use_safety_constraints /*check for bad faces*/,
285 false /*apply moves as soon as they're calculated*/,
286 false /*do not enforce a minimum angle improvement*/);
287 if(do_project)
288 {
289 if(use_safety_constraints && does_self_intersect(tmesh))
290 {
291 #ifdef CGAL_PMP_SMOOTHING_DEBUG
292 std::cerr << "Cannot re-project as there are self-intersections in the mesh!\n";
293 #endif
294 break;
295 }
296
297 area_smoother.project_to_surface(aabb_tree);
298 }
299
300 if(use_Delaunay_flips)
301 delaunay_flipper(faces);
302 }
303
304 // ... then angle smoothing
305 if(use_angle_smoothing)
306 {
307 #ifdef CGAL_PMP_SMOOTHING_DEBUG
308 std::cout << "Smooth angles..." << std::endl;
309 #endif
310
311 angle_smoother.optimize(use_safety_constraints /*check for bad faces*/,
312 true /*apply all moves at once*/,
313 use_safety_constraints /*check if the min angle is improved*/);
314
315 if(do_project)
316 {
317 if(use_safety_constraints && does_self_intersect(tmesh))
318 {
319 #ifdef CGAL_PMP_SMOOTHING_DEBUG
320 std::cerr << "Can't do re-projection, there are self-intersections in the mesh!\n";
321 #endif
322 break;
323 }
324
325 angle_smoother.project_to_surface(aabb_tree);
326 }
327 }
328 }
329 }
330
331 ///\cond SKIP_IN_MANUAL
332 template <typename FaceRange, typename TriangleMesh>
smooth_mesh(const FaceRange & face_range,TriangleMesh & tmesh)333 void smooth_mesh(const FaceRange& face_range, TriangleMesh& tmesh)
334 {
335 smooth_mesh(face_range, tmesh, parameters::all_default());
336 }
337
338 template <typename TriangleMesh, typename CGAL_PMP_NP_TEMPLATE_PARAMETERS>
smooth_mesh(TriangleMesh & tmesh,const CGAL_PMP_NP_CLASS & np)339 void smooth_mesh(TriangleMesh& tmesh, const CGAL_PMP_NP_CLASS& np)
340 {
341 smooth_mesh(faces(tmesh), tmesh, np);
342 }
343
344 template<typename TriangleMesh>
smooth_mesh(TriangleMesh & tmesh)345 void smooth_mesh(TriangleMesh& tmesh)
346 {
347 smooth_mesh(faces(tmesh), tmesh, parameters::all_default());
348 }
349
350 template<typename TriangleMesh, typename GeomTraits, typename Stream>
angles_evaluation(TriangleMesh & tmesh,GeomTraits traits,Stream & output)351 void angles_evaluation(TriangleMesh& tmesh, GeomTraits traits, Stream& output)
352 {
353 internal::Quality_evaluator<TriangleMesh, GeomTraits> evaluator(tmesh, traits);
354 evaluator.gather_angles();
355 evaluator.extract_angles(output);
356 }
357
358 template<typename TriangleMesh, typename GeomTraits, typename Stream>
areas_evaluation(TriangleMesh & tmesh,GeomTraits traits,Stream & output)359 void areas_evaluation(TriangleMesh& tmesh, GeomTraits traits, Stream& output)
360 {
361 internal::Quality_evaluator<TriangleMesh, GeomTraits> evaluator(tmesh, traits);
362 evaluator.measure_areas();
363 evaluator.extract_areas(output);
364 }
365
366 template<typename TriangleMesh, typename GeomTraits, typename Stream>
aspect_ratio_evaluation(TriangleMesh & tmesh,GeomTraits traits,Stream & output)367 void aspect_ratio_evaluation(TriangleMesh& tmesh, GeomTraits traits, Stream& output)
368 {
369 internal::Quality_evaluator<TriangleMesh, GeomTraits> evaluator(tmesh, traits);
370 evaluator.calc_aspect_ratios();
371 evaluator.extract_aspect_ratios(output);
372 }
373 ///\endcond SKIP_IN_MANUAL
374
375 } // namespace Polygon_mesh_processing
376 } // namespace CGAL
377
378 #endif // CGAL_POLYGON_MESH_PROCESSING_SMOOTH_MESH_H
379