1 // Copyright (c) 2020 GeometryFactory (France) and Telecom Paris (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/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/smooth_vertices.h $ 7 // $Id: smooth_vertices.h 086299c 2021-01-08T10:39:24+01:00 Dmitry Anisimov 8 // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial 9 // 10 // 11 // Author(s) : Jane Tournois, Noura Faraj, Jean-Marc Thiery, Tamy Boubekeur 12 13 #ifndef CGAL_INTERNAL_SMOOTH_VERTICES_H 14 #define CGAL_INTERNAL_SMOOTH_VERTICES_H 15 16 #include <CGAL/license/Tetrahedral_remeshing.h> 17 18 #include <CGAL/Vector_3.h> 19 20 #include <CGAL/Tetrahedral_remeshing/internal/tetrahedral_remeshing_helpers.h> 21 #include <CGAL/Tetrahedral_remeshing/internal/FMLS.h> 22 23 #include <boost/unordered_map.hpp> 24 #include <boost/optional.hpp> 25 #include <boost/container/small_vector.hpp> 26 27 #include <vector> 28 #include <cmath> 29 #include <list> 30 31 namespace CGAL 32 { 33 namespace Tetrahedral_remeshing 34 { 35 namespace internal 36 { 37 template<typename C3t3> 38 class Tetrahedral_remeshing_smoother 39 { 40 typedef typename C3t3::Triangulation Tr; 41 typedef typename C3t3::Surface_patch_index Surface_patch_index; 42 typedef typename Tr::Vertex_handle Vertex_handle; 43 typedef typename Tr::Edge Edge; 44 typedef typename Tr::Facet Facet; 45 46 typedef typename Tr::Geom_traits Gt; 47 typedef typename Gt::Vector_3 Vector_3; 48 typedef typename Gt::Point_3 Point_3; 49 50 private: 51 typedef CGAL::Tetrahedral_remeshing::internal::FMLS<Gt> FMLS; 52 std::vector<FMLS> subdomain_FMLS; 53 boost::unordered_map<Surface_patch_index, std::size_t> subdomain_FMLS_indices; 54 bool m_smooth_constrained_edges; 55 56 public: 57 template<typename CellSelector> init(const C3t3 & c3t3,const CellSelector & cell_selector,const bool smooth_constrained_edges)58 void init(const C3t3& c3t3, 59 const CellSelector& cell_selector, 60 const bool smooth_constrained_edges) 61 { 62 //collect a map of vertices surface indices 63 boost::unordered_map<Vertex_handle, std::vector<Surface_patch_index> > vertices_surface_indices; 64 collect_vertices_surface_indices(c3t3, vertices_surface_indices); 65 66 //collect a map of normals at surface vertices 67 boost::unordered_map<Vertex_handle, 68 boost::unordered_map<Surface_patch_index, Vector_3> > vertices_normals; 69 compute_vertices_normals(c3t3, vertices_normals, cell_selector); 70 71 // Build MLS Surfaces 72 createMLSSurfaces(subdomain_FMLS, 73 subdomain_FMLS_indices, 74 vertices_normals, 75 vertices_surface_indices, 76 c3t3); 77 78 m_smooth_constrained_edges = smooth_constrained_edges; 79 } 80 81 private: 82 project_on_tangent_plane(const Vector_3 & gi,const Vector_3 & pi,const Vector_3 & normal)83 Vector_3 project_on_tangent_plane(const Vector_3& gi, 84 const Vector_3& pi, 85 const Vector_3& normal) 86 { 87 Vector_3 diff = pi - gi; 88 return gi + (normal * diff) * normal; 89 } 90 91 template<typename CellSelector> 92 boost::optional<Facet> find_adjacent_facet_on_surface(const Facet & f,const Edge & edge,const C3t3 & c3t3,const CellSelector & cell_selector)93 find_adjacent_facet_on_surface(const Facet& f, 94 const Edge& edge, 95 const C3t3& c3t3, 96 const CellSelector& cell_selector) 97 { 98 CGAL_assertion(is_boundary(c3t3, f, cell_selector)); 99 100 typedef typename Tr::Facet_circulator Facet_circulator; 101 102 if (c3t3.is_in_complex(edge)) 103 return {}; //do not "cross" complex edges 104 //they are likely to be sharp and not to follow the > 0 dot product criterion 105 106 const Surface_patch_index& patch = c3t3.surface_patch_index(f); 107 const Facet& mf = c3t3.triangulation().mirror_facet(f); 108 109 Facet_circulator fcirc = c3t3.triangulation().incident_facets(edge); 110 Facet_circulator fend = fcirc; 111 do 112 { 113 const Facet fi = *fcirc; 114 if (f != fi 115 && mf != fi 116 && is_boundary(c3t3, fi, cell_selector) 117 && patch == c3t3.surface_patch_index(fi)) 118 { 119 return canonical_facet(fi); //"canonical" is important 120 } 121 } while (++fcirc != fend); 122 123 return {}; 124 } 125 126 template<typename Gt> compute_normal(const Facet & f,const Vector_3 & reference_normal,const Gt & gt)127 Vector_3 compute_normal(const Facet& f, 128 const Vector_3& reference_normal, 129 const Gt& gt) 130 { 131 typename Gt::Construct_opposite_vector_3 132 opp = gt.construct_opposite_vector_3_object(); 133 typename Gt::Compute_scalar_product_3 134 scalar_product = gt.compute_scalar_product_3_object(); 135 136 Vector_3 n = CGAL::Tetrahedral_remeshing::normal(f, gt); 137 if (scalar_product(n, reference_normal) < 0.) 138 n = opp(n); 139 140 return n; 141 } 142 143 template<typename VertexNormalsMap, typename CellSelector> compute_vertices_normals(const C3t3 & c3t3,VertexNormalsMap & normals_map,const CellSelector & cell_selector)144 void compute_vertices_normals(const C3t3& c3t3, 145 VertexNormalsMap& normals_map, 146 const CellSelector& cell_selector) 147 { 148 typename Tr::Geom_traits gt = c3t3.triangulation().geom_traits(); 149 typename Tr::Geom_traits::Construct_opposite_vector_3 150 opp = gt.construct_opposite_vector_3_object(); 151 152 const Tr& tr = c3t3.triangulation(); 153 154 //collect all facet normals 155 boost::unordered_map<Facet, Vector_3> fnormals; 156 for (const Facet& f : tr.finite_facets()) 157 { 158 if (is_boundary(c3t3, f, cell_selector)) 159 { 160 const Facet cf = canonical_facet(f); 161 fnormals[cf] = CGAL::NULL_VECTOR; 162 } 163 } 164 165 for (const auto& fn : fnormals) 166 { 167 if(fn.second != CGAL::NULL_VECTOR) 168 continue; 169 170 const Facet& f = fn.first; 171 const Facet& mf = tr.mirror_facet(f); 172 CGAL_assertion(is_boundary(c3t3, f, cell_selector)); 173 174 Vector_3 start_ref = CGAL::Tetrahedral_remeshing::normal(f, tr.geom_traits()); 175 if (c3t3.triangulation().is_infinite(mf.first) 176 || c3t3.subdomain_index(mf.first) < c3t3.subdomain_index(f.first)) 177 start_ref = opp(start_ref); 178 fnormals[f] = start_ref; 179 180 std::list<Facet> facets; 181 facets.push_back(f); 182 while (!facets.empty()) 183 { 184 const Facet ff = facets.front(); 185 facets.pop_front(); 186 187 const typename C3t3::Cell_handle ch = f.first; 188 const std::array<std::array<int, 2>, 3> edges 189 = {{ {{(ff.second + 1) % 4, (ff.second + 2) % 4}}, //edge 1-2 190 {{(ff.second + 2) % 4, (ff.second + 3) % 4}}, //edge 2-3 191 {{(ff.second + 3) % 4, (ff.second + 1) % 4}} //edge 3-1 192 }}; //vertex indices in cells 193 194 const Vector_3& ref = fnormals[f]; 195 for (const std::array<int, 2>& ei : edges) 196 { 197 Edge edge(ch, ei[0], ei[1]); 198 if (boost::optional<Facet> neighbor 199 = find_adjacent_facet_on_surface(f, edge, c3t3, cell_selector)) 200 { 201 const Facet neigh = *neighbor; //already a canonical_facet 202 if (fnormals[neigh] == CGAL::NULL_VECTOR) //check it's not already computed 203 { 204 fnormals[neigh] = compute_normal(neigh, ref, gt); 205 facets.push_back(neigh); 206 } 207 } 208 } 209 } 210 } 211 212 #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG 213 std::ofstream osf("dump_facet_normals.polylines.txt"); 214 #endif 215 for (const auto& fn : fnormals) 216 { 217 const Facet& f = fn.first; 218 const Vector_3& n = fn.second; 219 220 #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG 221 typename Tr::Geom_traits::Point_3 fc 222 = CGAL::centroid(point(f.first->vertex(indices(f.second, 0))->point()), 223 point(f.first->vertex(indices(f.second, 1))->point()), 224 point(f.first->vertex(indices(f.second, 2))->point())); 225 osf << "2 " << fc << " " << (fc + n) << std::endl; 226 #endif 227 const Surface_patch_index& surf_i = c3t3.surface_patch_index(f); 228 229 for (int i = 0; i < 3; ++i) 230 { 231 const Vertex_handle vi = f.first->vertex(indices(f.second, i)); 232 typename VertexNormalsMap::iterator patch_vector_it = normals_map.find(vi); 233 234 if (patch_vector_it == normals_map.end() 235 || patch_vector_it->second.find(surf_i) == patch_vector_it->second.end()) 236 { 237 normals_map[vi][surf_i] = n; 238 } 239 else 240 { 241 normals_map[vi][surf_i] += n; 242 } 243 } 244 } 245 246 #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG 247 osf.close(); 248 std::ofstream os("dump_normals.polylines.txt"); 249 boost::unordered_map<Surface_patch_index, 250 std::vector<typename Tr::Geom_traits::Segment_3 > > ons_map; 251 #endif 252 253 //normalize the computed normals 254 for (typename VertexNormalsMap::iterator vnm_it = normals_map.begin(); 255 vnm_it != normals_map.end(); ++vnm_it) 256 { 257 //value type is map<Surface_patch_index, Vector_3> 258 for (typename VertexNormalsMap::mapped_type::iterator it = vnm_it->second.begin(); 259 it != vnm_it->second.end(); ++it) 260 { 261 Vector_3& n = it->second; 262 263 #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG 264 auto p = point(vnm_it->first->point()); 265 os << "2 " << p << " " << (p + n) << std::endl; 266 #endif 267 268 CGAL::Tetrahedral_remeshing::normalize(n, c3t3.triangulation().geom_traits()); 269 270 #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG 271 const Surface_patch_index si = it->first; 272 if (ons_map.find(si) == ons_map.end()) 273 ons_map[si] = std::vector<typename Tr::Geom_traits::Segment_3>(); 274 ons_map[si].push_back(typename Tr::Geom_traits::Segment_3(p, p + n)); 275 #endif 276 } 277 } 278 279 #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG 280 os.close(); 281 for (auto& kv : ons_map) 282 { 283 std::ostringstream oss; 284 oss << "dump_normals_normalized_" << kv.first << ".polylines.txt"; 285 std::ofstream ons(oss.str()); 286 for (auto s : kv.second) 287 ons << "2 " << s.source() << " " << s.target() << std::endl; 288 ons.close(); 289 } 290 #endif 291 } 292 project(const Surface_patch_index & si,const Vector_3 & gi)293 boost::optional<Vector_3> project(const Surface_patch_index& si, 294 const Vector_3& gi) 295 { 296 CGAL_assertion(subdomain_FMLS_indices.find(si) != subdomain_FMLS_indices.end()); 297 CGAL_assertion(!std::isnan(gi.x()) && !std::isnan(gi.y()) && !std::isnan(gi.z())); 298 299 Vector_3 point(gi.x(), gi.y(), gi.z()); 300 Vector_3 res_normal = CGAL::NULL_VECTOR; 301 Vector_3 result(point); 302 303 const FMLS& fmls = subdomain_FMLS[subdomain_FMLS_indices.at(si)]; 304 305 int it_nb = 0; 306 const int max_it_nb = 5; 307 const double epsilon = fmls.getPNScale() / 1000.; 308 const double sq_eps = CGAL::square(epsilon); 309 310 do 311 { 312 point = result; 313 314 fmls.fastProjectionCPU(point, result, res_normal); 315 316 if (std::isnan(result[0]) || std::isnan(result[1]) || std::isnan(result[2])) { 317 std::cout << "MLS error detected si " //<< si 318 << "\t(size : " << fmls.getPNSize() << ")" 319 << "\t(point = " << point << " )" << std::endl; 320 return {}; 321 } 322 } while ((result - point).squared_length() > sq_eps && ++it_nb < max_it_nb); 323 324 return Vector_3(result[0], result[1], result[2]); 325 } 326 327 template<typename CellRange, typename Tr> check_inversion_and_move(const typename Tr::Vertex_handle v,const typename Tr::Point & final_pos,const CellRange & inc_cells,const Tr &,double & total_move)328 bool check_inversion_and_move(const typename Tr::Vertex_handle v, 329 const typename Tr::Point& final_pos, 330 const CellRange& inc_cells, 331 const Tr& /* tr */, 332 #ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE 333 double& total_move) 334 #else 335 double&) 336 #endif 337 { 338 const typename Tr::Point backup = v->point(); //backup v's position 339 const typename Tr::Geom_traits::Point_3 pv = point(backup); 340 341 bool valid_orientation = false; 342 double frac = 1.0; 343 typename Tr::Geom_traits::Vector_3 move(pv, point(final_pos)); 344 345 do 346 { 347 v->set_point(typename Tr::Point(pv + frac * move)); 348 349 bool valid_try = true; 350 for (const typename Tr::Cell_handle& ci : inc_cells) 351 { 352 if (CGAL::POSITIVE != CGAL::orientation(point(ci->vertex(0)->point()), 353 point(ci->vertex(1)->point()), 354 point(ci->vertex(2)->point()), 355 point(ci->vertex(3)->point()))) 356 { 357 frac = 0.9 * frac; 358 valid_try = false; 359 break; 360 } 361 } 362 valid_orientation = valid_try; 363 } 364 while(!valid_orientation && frac > 0.1); 365 366 if (!valid_orientation) //move failed 367 v->set_point(backup); 368 369 #ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE 370 else 371 total_move += CGAL::approximate_sqrt(CGAL::squared_distance(pv, point(v->point()))); 372 #endif 373 374 return valid_orientation; 375 } 376 collect_vertices_surface_indices(const C3t3 & c3t3,boost::unordered_map<Vertex_handle,std::vector<Surface_patch_index>> & vertices_surface_indices)377 void collect_vertices_surface_indices( 378 const C3t3& c3t3, 379 boost::unordered_map<Vertex_handle, 380 std::vector<Surface_patch_index> >& vertices_surface_indices) 381 { 382 for (typename C3t3::Facets_in_complex_iterator 383 fit = c3t3.facets_in_complex_begin(); 384 fit != c3t3.facets_in_complex_end(); ++fit) 385 { 386 const Surface_patch_index& surface_index = c3t3.surface_patch_index(*fit); 387 388 for (int i = 0; i < 3; i++) 389 { 390 const Vertex_handle vi = fit->first->vertex(indices(fit->second, i)); 391 392 std::vector<Surface_patch_index>& v_surface_indices = vertices_surface_indices[vi]; 393 if (std::find(v_surface_indices.begin(), v_surface_indices.end(), surface_index) == v_surface_indices.end()) 394 v_surface_indices.push_back(surface_index); 395 } 396 } 397 } 398 399 public: 400 template<typename C3T3, typename CellSelector> smooth_vertices(C3T3 & c3t3,const bool protect_boundaries,const CellSelector & cell_selector)401 void smooth_vertices(C3T3& c3t3, 402 const bool protect_boundaries, 403 const CellSelector& cell_selector) 404 { 405 typedef typename C3T3::Cell_handle Cell_handle; 406 typedef typename Gt::FT FT; 407 408 #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG 409 std::ofstream os_surf("smooth_surfaces.polylines.txt"); 410 std::ofstream os_surf0("smooth_surfaces0.polylines.txt"); 411 std::ofstream os_vol("smooth_volume.polylines.txt"); 412 #endif 413 414 #ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE 415 std::cout << "Smooth vertices..."; 416 std::cout.flush(); 417 #endif 418 std::size_t nb_done_3d = 0; 419 std::size_t nb_done_2d = 0; 420 std::size_t nb_done_1d = 0; 421 FT total_move = 0.; 422 423 Tr& tr = c3t3.triangulation(); 424 425 //collect a map of vertices surface indices 426 boost::unordered_map<Vertex_handle, std::vector<Surface_patch_index> > vertices_surface_indices; 427 if(m_smooth_constrained_edges) 428 collect_vertices_surface_indices(c3t3, vertices_surface_indices); 429 430 //collect a map of normals at surface vertices 431 boost::unordered_map<Vertex_handle, 432 boost::unordered_map<Surface_patch_index, Vector_3> > vertices_normals; 433 compute_vertices_normals(c3t3, vertices_normals, cell_selector); 434 435 //smooth() 436 const std::size_t nbv = tr.number_of_vertices(); 437 boost::unordered_map<Vertex_handle, std::size_t> vertex_id; 438 std::vector<Vector_3> smoothed_positions(nbv, CGAL::NULL_VECTOR); 439 std::vector<int> neighbors(nbv, -1); 440 std::vector<bool> free_vertex(nbv, false);//are vertices free to move? indices are in `vertex_id` 441 442 //collect ids 443 std::size_t id = 0; 444 for (const Vertex_handle v : tr.finite_vertex_handles()) 445 { 446 vertex_id[v] = id++; 447 } 448 449 //collect incident cells 450 std::vector<boost::container::small_vector<Cell_handle, 40> > 451 inc_cells(nbv, boost::container::small_vector<Cell_handle, 40>()); 452 for (const Cell_handle c : tr.finite_cell_handles()) 453 { 454 const bool cell_is_selected = cell_selector(c); 455 456 for (int i = 0; i < 4; ++i) 457 { 458 const std::size_t idi = vertex_id[c->vertex(i)]; 459 inc_cells[idi].push_back(c); 460 if(cell_is_selected) 461 free_vertex[idi] = true; 462 } 463 } 464 465 if (!protect_boundaries && m_smooth_constrained_edges) 466 { 467 /////////////// EDGES IN COMPLEX ////////////////// 468 //collect neighbors 469 for (const Edge& e : tr.finite_edges()) 470 { 471 if (c3t3.is_in_complex(e)) 472 { 473 const Vertex_handle vh0 = e.first->vertex(e.second); 474 const Vertex_handle vh1 = e.first->vertex(e.third); 475 476 const std::size_t& i0 = vertex_id.at(vh0); 477 const std::size_t& i1 = vertex_id.at(vh1); 478 479 const bool on_feature_v0 = is_on_feature(vh0); 480 const bool on_feature_v1 = is_on_feature(vh1); 481 482 if (!c3t3.is_in_complex(vh0)) 483 neighbors[i0] = (std::max)(0, neighbors[i0]); 484 if (!c3t3.is_in_complex(vh1)) 485 neighbors[i1] = (std::max)(0, neighbors[i1]); 486 487 if (!c3t3.is_in_complex(vh0) && on_feature_v1) 488 { 489 const Point_3& p1 = point(vh1->point()); 490 smoothed_positions[i0] = smoothed_positions[i0] + Vector_3(p1.x(), p1.y(), p1.z()); 491 neighbors[i0]++; 492 } 493 if (!c3t3.is_in_complex(vh1) && on_feature_v0) 494 { 495 const Point_3& p0 = point(vh0->point()); 496 smoothed_positions[i1] = smoothed_positions[i1] + Vector_3(p0.x(), p0.y(), p0.z()); 497 neighbors[i1]++; 498 } 499 } 500 } 501 502 // Smooth 503 for (Vertex_handle v : tr.finite_vertex_handles()) 504 { 505 const std::size_t& vid = vertex_id.at(v); 506 if (!free_vertex[vid]) 507 continue; 508 509 if (neighbors[vid] > 1) 510 { 511 Vector_3 smoothed_position = smoothed_positions[vid] / neighbors[vid]; 512 Vector_3 final_position = CGAL::NULL_VECTOR; 513 514 std::size_t count = 0; 515 const Vector_3 current_pos(CGAL::ORIGIN, point(v->point())); 516 517 const std::vector<Surface_patch_index>& v_surface_indices = vertices_surface_indices[v]; 518 for (const Surface_patch_index& si : v_surface_indices) 519 { 520 Vector_3 normal_projection 521 = project_on_tangent_plane(smoothed_position, current_pos, vertices_normals[v][si]); 522 523 //Check if the mls surface exists to avoid degenerated cases 524 if (boost::optional<Vector_3> mls_projection = project(si, normal_projection)) { 525 final_position = final_position + *mls_projection; 526 } 527 else { 528 final_position = final_position + normal_projection; 529 } 530 count++; 531 } 532 533 if (count > 0) 534 final_position = final_position / static_cast<FT>(count); 535 else 536 final_position = smoothed_position; 537 538 #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG 539 os_surf << "2 " << current_pos << " " << final_position << std::endl; 540 #endif 541 // move vertex 542 const typename Tr::Point new_pos(final_position.x(), final_position.y(), final_position.z()); 543 if(check_inversion_and_move(v, new_pos, inc_cells[vid], tr, total_move)) 544 nb_done_1d++; 545 } 546 else if (neighbors[vid] > 0) 547 { 548 Vector_3 final_position = CGAL::NULL_VECTOR; 549 550 int count = 0; 551 const Vector_3 current_pos(CGAL::ORIGIN, point(v->point())); 552 553 const std::vector<Surface_patch_index>& v_surface_indices = vertices_surface_indices[v]; 554 for (const Surface_patch_index& si : v_surface_indices) 555 { 556 //Check if the mls surface exists to avoid degenerated cases 557 558 if (boost::optional<Vector_3> mls_projection = project(si, current_pos)) { 559 final_position = final_position + *mls_projection; 560 } 561 else { 562 final_position = final_position + current_pos; 563 } 564 count++; 565 } 566 567 if (count > 0) 568 final_position = final_position / static_cast<FT>(count); 569 else 570 final_position = current_pos; 571 572 #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG 573 os_surf << "2 " << current_pos << " " << final_position << std::endl; 574 #endif 575 // move vertex 576 const typename Tr::Point new_pos(final_position.x(), final_position.y(), final_position.z()); 577 if(check_inversion_and_move(v, new_pos, inc_cells[vid], tr, total_move)) 578 nb_done_1d++; 579 } 580 } 581 } 582 583 smoothed_positions.assign(nbv, CGAL::NULL_VECTOR); 584 neighbors.assign(nbv, -1); 585 586 /////////////// EDGES ON SURFACE, BUT NOT IN COMPLEX ////////////////// 587 if (!protect_boundaries) 588 { 589 for (const Edge& e : tr.finite_edges()) 590 { 591 if (is_boundary(c3t3, e, cell_selector) && !c3t3.is_in_complex(e)) 592 { 593 const Vertex_handle vh0 = e.first->vertex(e.second); 594 const Vertex_handle vh1 = e.first->vertex(e.third); 595 596 const std::size_t& i0 = vertex_id.at(vh0); 597 const std::size_t& i1 = vertex_id.at(vh1); 598 599 const bool on_feature_v0 = is_on_feature(vh0); 600 const bool on_feature_v1 = is_on_feature(vh1); 601 602 if (!on_feature_v0) 603 neighbors[i0] = (std::max)(0, neighbors[i0]); 604 if (!on_feature_v1) 605 neighbors[i1] = (std::max)(0, neighbors[i1]); 606 607 if (!on_feature_v0) 608 { 609 const Point_3& p1 = point(vh1->point()); 610 smoothed_positions[i0] = smoothed_positions[i0] + Vector_3(p1.x(), p1.y(), p1.z()); 611 neighbors[i0]++; 612 } 613 if (!on_feature_v1) 614 { 615 const Point_3& p0 = point(vh0->point()); 616 smoothed_positions[i1] = smoothed_positions[i1] + Vector_3(p0.x(), p0.y(), p0.z()); 617 neighbors[i1]++; 618 } 619 } 620 } 621 622 for (Vertex_handle v : tr.finite_vertex_handles()) 623 { 624 const std::size_t& vid = vertex_id.at(v); 625 if (!free_vertex[vid] || v->in_dimension() != 2) 626 continue; 627 628 if (neighbors[vid] > 1) 629 { 630 Vector_3 smoothed_position = smoothed_positions[vid] / static_cast<FT>(neighbors[vid]); 631 const Vector_3 current_pos(CGAL::ORIGIN, point(v->point())); 632 Vector_3 final_position = CGAL::NULL_VECTOR; 633 634 const Surface_patch_index si = surface_patch_index(v, c3t3); 635 CGAL_assertion(si != Surface_patch_index()); 636 637 Vector_3 normal_projection = project_on_tangent_plane(smoothed_position, 638 current_pos, 639 vertices_normals[v][si]); 640 641 if (boost::optional<Vector_3> mls_projection = project(si, normal_projection)) 642 final_position = final_position + *mls_projection; 643 else 644 final_position = smoothed_position; 645 646 #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG 647 os_surf << "2 " << current_pos << " " << final_position << std::endl; 648 #endif 649 const typename Tr::Point new_pos(final_position.x(), final_position.y(), final_position.z()); 650 if(check_inversion_and_move(v, new_pos, inc_cells[vid], tr, total_move)) 651 nb_done_2d++; 652 } 653 else if (neighbors[vid] > 0) 654 { 655 const Surface_patch_index si = surface_patch_index(v, c3t3); 656 CGAL_assertion(si != Surface_patch_index()); 657 658 const Vector_3 current_pos(CGAL::ORIGIN, point(v->point())); 659 660 if (boost::optional<Vector_3> mls_projection = project(si, current_pos)) 661 { 662 const typename Tr::Point new_pos(CGAL::ORIGIN + *mls_projection); 663 if(check_inversion_and_move(v, new_pos, inc_cells[vid], tr, total_move)) 664 nb_done_2d++; 665 666 #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG 667 os_surf0 << "2 " << current_pos << " " << new_pos << std::endl; 668 #endif 669 } 670 } 671 } 672 } 673 CGAL_assertion(CGAL::Tetrahedral_remeshing::debug::are_cell_orientations_valid(tr)); 674 //// end if(!protect_boundaries) 675 676 smoothed_positions.assign(nbv, CGAL::NULL_VECTOR); 677 neighbors.assign(nbv, 0/*for dim 3 vertices, start counting directly from 0*/); 678 679 ////////////// INTERNAL VERTICES /////////////////////// 680 for (const Edge& e : tr.finite_edges()) 681 { 682 if (!is_outside(e, c3t3, cell_selector)) 683 { 684 const Vertex_handle vh0 = e.first->vertex(e.second); 685 const Vertex_handle vh1 = e.first->vertex(e.third); 686 687 const std::size_t& i0 = vertex_id.at(vh0); 688 const std::size_t& i1 = vertex_id.at(vh1); 689 690 if (c3t3.in_dimension(vh0) == 3) 691 { 692 const Point_3& p1 = point(vh1->point()); 693 smoothed_positions[i0] = smoothed_positions[i0] + Vector_3(CGAL::ORIGIN, p1); 694 neighbors[i0]++; 695 } 696 if (c3t3.in_dimension(vh1) == 3) 697 { 698 const Point_3& p0 = point(vh0->point()); 699 smoothed_positions[i1] = smoothed_positions[i1] + Vector_3(CGAL::ORIGIN, p0); 700 neighbors[i1]++; 701 } 702 } 703 } 704 705 for (Vertex_handle v : tr.finite_vertex_handles()) 706 { 707 const std::size_t& vid = vertex_id.at(v); 708 if (!free_vertex[vid]) 709 continue; 710 711 if (c3t3.in_dimension(v) == 3 && neighbors[vid] > 1) 712 { 713 #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG 714 os_vol << "2 " << point(v->point()); 715 #endif 716 const Vector_3 p = smoothed_positions[vid] / static_cast<FT>(neighbors[vid]); 717 typename Tr::Point new_pos(p.x(), p.y(), p.z()); 718 if(check_inversion_and_move(v, new_pos, inc_cells[vid], tr, total_move)) 719 nb_done_3d++; 720 721 #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG 722 os_vol << " " << point(v->point()) << std::endl; 723 #endif 724 } 725 } 726 CGAL_assertion(CGAL::Tetrahedral_remeshing::debug::are_cell_orientations_valid(tr)); 727 728 #ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE 729 std::size_t nb_done = nb_done_3d + nb_done_2d + nb_done_1d; 730 std::cout << " done (" 731 << nb_done_3d << "/" << nb_done_2d << "/" << nb_done_1d << " vertices smoothed," 732 << " average move = " << (total_move / nb_done) 733 << ")." << std::endl; 734 #endif 735 #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG 736 CGAL::Tetrahedral_remeshing::debug::dump_vertices_by_dimension( 737 c3t3.triangulation(), "c3t3_vertices_after_smoothing"); 738 os_surf.close(); 739 os_vol.close(); 740 #endif 741 } 742 743 };//end class Tetrahedral_remeshing_smoother 744 }//namespace internal 745 }//namespace Tetrahedral_adaptive_remeshing 746 }//namespace CGAL 747 748 #endif //CGAL_INTERNAL_SMOOTH_VERTICES_H 749