Developer Documentation
BaseRemesherT_impl.hh
1 /*===========================================================================*\
2 * *
3 * OpenFlipper *
4  * Copyright (c) 2001-2015, RWTH-Aachen University *
5  * Department of Computer Graphics and Multimedia *
6  * All rights reserved. *
7  * www.openflipper.org *
8  * *
9  *---------------------------------------------------------------------------*
10  * This file is part of OpenFlipper. *
11  *---------------------------------------------------------------------------*
12  * *
13  * Redistribution and use in source and binary forms, with or without *
14  * modification, are permitted provided that the following conditions *
15  * are met: *
16  * *
17  * 1. Redistributions of source code must retain the above copyright notice, *
18  * this list of conditions and the following disclaimer. *
19  * *
20  * 2. Redistributions in binary form must reproduce the above copyright *
21  * notice, this list of conditions and the following disclaimer in the *
22  * documentation and/or other materials provided with the distribution. *
23  * *
24  * 3. Neither the name of the copyright holder nor the names of its *
25  * contributors may be used to endorse or promote products derived from *
26  * this software without specific prior written permission. *
27  * *
28  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS *
29  * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED *
30  * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A *
31  * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER *
32  * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, *
33  * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, *
34  * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR *
35  * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF *
36  * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING *
37  * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS *
38  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *
39 * *
40 \*===========================================================================*/
41 
42 //=============================================================================
43 //
44 // CLASS BaseRemesherT - IMPLEMENTATION
45 //
46 //=============================================================================
47 
48 #define BASE_REMESHERT_C
49 
50 
51 //== INCLUDES =================================================================
52 
53 
54 #include "BaseRemesherT.hh"
55 
56 // Remshing
57 #include "DiffGeoT.hh"
58 
59 // OpenMesh
60 #include <OpenMesh/Core/Utils/Property.hh>
62 
63 // Selection Stuff
65 
66 
67 //== NAMESPACES ===============================================================
68 
69 namespace Remeshing {
70 
71 
72 //== IMPLEMENTATION ==========================================================
73 
74 
75 template <class Mesh>
76 BaseRemesherT<Mesh>::
77 BaseRemesherT(Mesh& _mesh, ProgressEmitter* _progress)
78  : mesh_(_mesh), refmesh_(0), bsp_(0), nothing_selected_(true),progress_(_progress) {
79 
80 }
81 
82 
83 //-----------------------------------------------------------------------------
84 
85 
86 template <class Mesh>
87 BaseRemesherT<Mesh>::
88 ~BaseRemesherT() {
89 
90  delete bsp_;
91  delete refmesh_;
92 }
93 
94 
95 //-----------------------------------------------------------------------------
96 
97 
98 template <class Mesh>
99 void
100 BaseRemesherT<Mesh>::
101 init_reference()
102 {
103  typename Mesh::VIter v_it, v_end;
104  typename Mesh::FIter f_it, f_end;
105  typename Mesh::CFVIter fv_it;
106  typename Mesh::VHandle v0, v1, v2;
107 
108 
109  // clear old stuff first
110  if (bsp_) delete_reference();
111 
112  // tag non-locked vertices + one ring
113  for (v_it=mesh_.vertices_begin(), v_end=mesh_.vertices_end();
114  v_it!=v_end; ++v_it)
115  mesh_.status(*v_it).set_tagged(!mesh_.status(*v_it).locked());
116 
118 
119  // create reference mesh
120  refmesh_ = new Mesh();
121 
123  mesh_.add_property(vhandles);
124 
125  for (v_it=mesh_.vertices_begin(), v_end=mesh_.vertices_end();
126  v_it!=v_end; ++v_it)
127  if (mesh_.status(*v_it).tagged())
128  mesh_.property(vhandles, *v_it) = refmesh_->add_vertex(mesh_.point(*v_it));
129 
130  for (f_it=mesh_.faces_begin(), f_end=mesh_.faces_end();
131  f_it!=f_end; ++f_it)
132  {
133  fv_it = mesh_.cfv_iter(*f_it);
134  v0 = *fv_it;
135  v1 = *(++fv_it);
136  v2 = *(++fv_it);
137 
138  if (mesh_.status(v0).tagged() &&
139  mesh_.status(v1).tagged() &&
140  mesh_.status(v2).tagged())
141  {
142  refmesh_->add_face( mesh_.property(vhandles, v0),
143  mesh_.property(vhandles, v1),
144  mesh_.property(vhandles, v2) );
145  }
146  }
147 
149  mesh_.remove_property(vhandles);
150 
151  // need vertex normals
152  refmesh_->request_face_normals();
153  refmesh_->request_vertex_normals();
154  refmesh_->update_normals(/*Mesh::ANGLE_WEIGHTED*/);
155 
156 
157  // build BSP
158  bsp_ = new BSP(*refmesh_);
159  bsp_->reserve(refmesh_->n_faces());
160  for (f_it=refmesh_->faces_begin(), f_end=refmesh_->faces_end();
161  f_it!=f_end; ++f_it)
162  bsp_->push_back(*f_it);
163  bsp_->build(10, 100);
164 }
165 
166 
167 //-----------------------------------------------------------------------------
168 
169 
170 template <class Mesh>
171 void
172 BaseRemesherT<Mesh>::
173 delete_reference()
174 {
175  delete bsp_; bsp_ = 0;
176  delete refmesh_; refmesh_ = 0;
177 }
178 
179 
180 //-----------------------------------------------------------------------------
181 
182 
183 template <class Mesh>
184 void
185 BaseRemesherT<Mesh>::
186 project_to_reference(typename Mesh::VertexHandle _vh) const
187 {
188  typename Mesh::Point p = mesh_.point(_vh);
189  typename Mesh::FaceHandle fh = bsp_->nearest(p).handle;
190 
191  if ( ! fh.is_valid() )
192  return;
193 
194  typename Mesh::CFVIter fv_it = refmesh_->cfv_iter(fh);
195 
196 
197  const typename Mesh::Point& p0 = refmesh_->point(*fv_it);
198  typename Mesh::Normal n0 = refmesh_->normal(*fv_it);
199  const typename Mesh::Point& p1 = refmesh_->point(*(++fv_it));
200  typename Mesh::Normal n1 = refmesh_->normal(*fv_it);
201  const typename Mesh::Point& p2 = refmesh_->point(*(++fv_it));
202  typename Mesh::Normal n2 = refmesh_->normal(*fv_it);
203 
204 
205  // project
206  ACG::Geometry::distPointTriangle(p, p0, p1, p2, p);
207  mesh_.set_point(_vh, p);
208 
209 
210  // get barycentric coordinates
211  if (!ACG::Geometry::baryCoord(p, p0, p1, p2, p))
212  p[0] = p[1] = p[2] = 1.0/3.0;
213 
214 
215  // interpolate normal
216  typename Mesh::Normal n;
217  n = (n0 *= p[0]);
218  n += (n1 *= p[1]);
219  n += (n2 *= p[2]);
220  n.normalize();
221  mesh_.set_normal(_vh, n);
222 }
223 
224 
225 //-----------------------------------------------------------------------------
226 
227 
228 template <class Mesh>
229 void
230 BaseRemesherT<Mesh>::
231 remesh(unsigned int _iters,
232  unsigned int _area_iters,
233  bool _use_projection,
234  Selection _selection) {
235 
236  try
237  {
238  if (_selection == VERTEX_SELECTION)
240  else if (_selection == FACE_SELECTION)
242  remeshh(_iters, _area_iters, _use_projection);
243  }
244  catch (std::bad_alloc&)
245  {
246  mesh_.clear();
247  omerr() << "Remeshig: Out of memory\n";
248  }
249 
250  cleanup();
251 }
252 
253 //-----------------------------------------------------------------------------
254 
255 template <class Mesh>
256 void
259 {
260  typename Mesh::EIter e_it, e_end;
261  typename Mesh::VIter v_it, v_end;
262  typename Mesh::FIter f_it, f_end;
263  typename Mesh::CFVIter fv_it;
264  typename Mesh::CVOHIter vh_it;
265  typename Mesh::VHandle v0, v1;
266  typename Mesh::FHandle f0, f1;
267 
268 
269  // need vertex and edge status
270  mesh_.request_vertex_status();
271  mesh_.request_edge_status();
272  mesh_.request_face_status();
273 
274 
275 
276  // if nothing selected -> select all
277  nothing_selected_ = true;
278  for (v_it=mesh_.vertices_begin(), v_end=mesh_.vertices_end();
279  v_it!=v_end; ++v_it)
280  if (mesh_.status(*v_it).selected())
281  { nothing_selected_ = false; break; }
282 
283  if (nothing_selected_)
284  for (v_it=mesh_.vertices_begin(), v_end=mesh_.vertices_end();
285  v_it!=v_end; ++v_it)
286  mesh_.status(*v_it).set_selected(true);
287 
288 
289 
290  // lock un-selected vertices & edges
291  for (v_it=mesh_.vertices_begin(), v_end=mesh_.vertices_end();
292  v_it!=v_end; ++v_it)
293  mesh_.status(*v_it).set_locked(!mesh_.status(*v_it).selected());
294 
295  for (e_it=mesh_.edges_begin(), e_end=mesh_.edges_end();
296  e_it!=e_end; ++e_it)
297  {
298  v0 = mesh_.to_vertex_handle(mesh_.halfedge_handle(*e_it, 0));
299  v1 = mesh_.to_vertex_handle(mesh_.halfedge_handle(*e_it, 1));
300  mesh_.status(*e_it).set_locked(mesh_.status(v0).locked() ||
301  mesh_.status(v1).locked());
302  }
303 
304 
305 
306  // handle feature corners:
307  // lock corner vertices (>2 feature edges) and endpoints
308  // of feature lines (1 feature edge)
309  for (v_it=mesh_.vertices_begin(), v_end=mesh_.vertices_end();
310  v_it!=v_end; ++v_it)
311  {
312  if (mesh_.status(*v_it).feature())
313  {
314  int c=0;
315  for (vh_it=mesh_.cvoh_iter(*v_it); vh_it.is_valid(); ++vh_it)
316  if (mesh_.status(mesh_.edge_handle(*vh_it)).feature())
317  ++c;
318  if (c!=2) mesh_.status(*v_it).set_locked(true);
319  }
320  }
321 
322 
323 
324  // build reference mesh
325  init_reference();
326 // if (emit_progress_) Progress().step(5);
327 
328 
329  // add properties
330  mesh_.add_property(valences_);
331  mesh_.add_property(update_);
332  mesh_.add_property(area_);
333 }
334 
335 //-----------------------------------------------------------------------------
336 
337 
338 template <class Mesh>
339 void
342 {
343  typename Mesh::EIter e_it, e_end;
344  typename Mesh::VIter v_it, v_end;
345  typename Mesh::FIter f_it, f_end;
346  typename Mesh::CFVIter fv_it;
347  typename Mesh::CVOHIter vh_it;
348  typename Mesh::VHandle v0, v1;
349  typename Mesh::FHandle f0, f1;
350 
351 
352  // need vertex and edge status
353  mesh_.request_vertex_status();
354  mesh_.request_edge_status();
355  mesh_.request_face_status();
356 
357  // if nothing selected -> select all
358  nothing_selected_ = true;
359  for (f_it = mesh_.faces_begin(), f_end = mesh_.faces_end(); f_it != f_end;
360  ++f_it)
361  {
362  if (mesh_.status(*f_it).selected())
363  {
364  nothing_selected_ = false;
365  break;
366  }
367  }
368 
369  if (nothing_selected_)
370  MeshSelection::selectAllFaces(&mesh_);
371 
372 
373 
374  // lock un-selected vertices & edges
375  for (v_it = mesh_.vertices_begin(), v_end = mesh_.vertices_end();
376  v_it != v_end; ++v_it)
377  {
378  bool all_faces_selected = true;
379 
380  for (typename Mesh::ConstVertexFaceIter vf_it = mesh_.cvf_iter(*v_it);
381  vf_it.is_valid(); ++vf_it)
382  {
383  if (!mesh_.status(*vf_it).selected())
384  {
385  all_faces_selected = false;
386  break;
387  }
388  }
389  mesh_.status(*v_it).set_locked(!all_faces_selected);
390  }
391 
392  for (e_it = mesh_.edges_begin(), e_end = mesh_.edges_end();
393  e_it != e_end; ++e_it)
394  {
395  if (mesh_.is_boundary(*e_it))
396  {
397  mesh_.status(*e_it).set_locked(true);
398  }
399  else
400  {
401  f0 = mesh_.face_handle(mesh_.halfedge_handle(*e_it, 0));
402  f1 = mesh_.face_handle(mesh_.halfedge_handle(*e_it, 1));
403 
404  mesh_.status(*e_it).set_locked(!(mesh_.status(f0).selected() && mesh_.status(f1).selected()));
405  }
406  }
407 
408  MeshSelection::clearFaceSelection(&mesh_);
409 
410  // handle feature corners:
411  // lock corner vertices (>2 feature edges) and endpoints
412  // of feature lines (1 feature edge)
413  for (v_it=mesh_.vertices_begin(), v_end=mesh_.vertices_end();
414  v_it!=v_end; ++v_it)
415  {
416  if (mesh_.status(*v_it).feature())
417  {
418  int c=0;
419  for (vh_it=mesh_.cvoh_iter(*v_it); vh_it.is_valid(); ++vh_it)
420  if (mesh_.status(mesh_.edge_handle(*vh_it)).feature())
421  ++c;
422  if (c!=2) mesh_.status(*v_it).set_locked(true);
423  }
424  }
425 
426 
427 
428  // build reference mesh
429  init_reference();
430 // if (emit_progress_) Progress().step(5);
431 
432 
433  // add properties
434  mesh_.add_property(valences_);
435  mesh_.add_property(update_);
436  mesh_.add_property(area_);
437 }
438 
439 
440 //-----------------------------------------------------------------------------
441 
442 
443 template <class Mesh>
444 void
446 remeshh(unsigned int _iters,
447  unsigned int _area_iters, bool _use_projection) {
448 
449  double progress_step = 100.0 / (((double)_iters/4.0 + (double)_area_iters));
450  double total_progress = 0.0;
451 
452  for (unsigned int i = 0; i < _iters; ++i) {
453 
454  split_long_edges();
455  if(progress_ != NULL) {
456  total_progress += progress_step;
457  progress_->sendProgressSignal(total_progress);
458  }
459 
460  collapse_short_edges();
461  if(progress_ != NULL) {
462  total_progress += progress_step;
463  progress_->sendProgressSignal(total_progress);
464  }
465 
466  flip_edges();
467  if(progress_ != NULL) {
468  total_progress += progress_step;
469  progress_->sendProgressSignal(total_progress);
470  }
471 
472  tangential_smoothing(_use_projection);
473  if(progress_ != NULL) {
474  total_progress += progress_step;
475  progress_->sendProgressSignal(total_progress);
476  }
477  }
478 
479  if (_area_iters) {
480  balanace_area(_area_iters, _use_projection);
481  if(progress_ != NULL) {
482  total_progress += progress_step;
483  progress_->sendProgressSignal(total_progress);
484  }
485  }
486 
487  // feature edges block flips -> might lead to caps
488  remove_caps();
489 }
490 
491 
492 //-----------------------------------------------------------------------------
493 
494 
495 template <class Mesh>
496 void
498 cleanup()
499 {
500  typename Mesh::EIter e_it, e_end;
501  typename Mesh::VIter v_it, v_end;
502  typename Mesh::FIter f_it, f_end;
503  typename Mesh::CFVIter fv_it;
504  typename Mesh::CVOHIter vh_it;
505  typename Mesh::VHandle v0, v1;
506  typename Mesh::FHandle f0, f1;
507 
508 
509  // remove properties
510  mesh_.remove_property(valences_);
511  mesh_.remove_property(update_);
512  mesh_.remove_property(area_);
513 
514 
515  // remove reference again
516  delete_reference();
517 
518 
519  // unlock all vertices & edges
520  for (v_it=mesh_.vertices_begin(), v_end=mesh_.vertices_end();
521  v_it!=v_end; ++v_it)
522  {
523  mesh_.status(*v_it).set_locked(false);
524  }
525  for (e_it=mesh_.edges_begin(), e_end=mesh_.edges_end();
526  e_it!=e_end; ++e_it)
527  {
528  mesh_.status(*e_it).set_locked(false);
529  }
530 
531 
532  // de-select if nothing was selected before
533  if (nothing_selected_)
534  for (v_it=mesh_.vertices_begin(), v_end=mesh_.vertices_end();
535  v_it!=v_end; ++v_it)
536  mesh_.status(*v_it).set_selected(false);
537 
538 
539 
540  // free vertex and edge status
541  mesh_.release_vertex_status();
542  mesh_.release_edge_status();
543  mesh_.release_face_status();
544 }
545 
546 
547 //-----------------------------------------------------------------------------
548 
549 
550 template <class Mesh>
551 void
554 {
555  typename Mesh::VIter v_it, v_end;
556  typename Mesh::EIter e_it, e_end;
557  typename Mesh::VHandle v0, v1, vh;
558  typename Mesh::EHandle eh, e0, e1;
559  typename Mesh::FHandle f0, f1, f2, f3;
560 
561  bool ok, is_feature, is_boundary;
562  int i;
563 
564  // un-tagg
565  for (v_it = mesh_.vertices_begin(), v_end = mesh_.vertices_end(); v_it != v_end; ++v_it)
566  mesh_.status(*v_it).set_tagged(false);
567 
568  // handle Nastran PIDs during refinement
570  mesh_.get_property_handle(pids, "Nastran PIDs");
571 
572  // split long edges
573  for (ok = false, i = 0; !ok && i < 100; ++i) {
574  ok = true;
575 
576  for (e_it = mesh_.edges_begin(), e_end = mesh_.edges_end(); e_it != e_end; ++e_it) {
577  v0 = mesh_.to_vertex_handle(mesh_.halfedge_handle(*e_it, 0));
578  v1 = mesh_.to_vertex_handle(mesh_.halfedge_handle(*e_it, 1));
579 
580  if (!mesh_.status(*e_it).locked() && is_too_long(v0, v1)) {
581  const typename Mesh::Point& p0 = mesh_.point(v0);
582  const typename Mesh::Point& p1 = mesh_.point(v1);
583 
584  is_feature = mesh_.status(*e_it).feature();
585  is_boundary = mesh_.is_boundary(*e_it);
586 
587  vh = mesh_.add_vertex((p0 + p1) * 0.5);
588  mesh_.split(*e_it, vh);
589 
590  mesh_.status(vh).set_selected(true);
591 
592  if (is_feature) {
593  eh = (is_boundary ? mesh_.edge_handle(mesh_.n_edges() - 2) : mesh_.edge_handle(mesh_.n_edges() - 3));
594 
595  mesh_.status(eh).set_feature(true);
596  } else {
597  project_to_reference(vh);
598  }
599 
600  if (pids.is_valid()) {
601  e0 = *e_it;
602  e1 = (is_boundary ? mesh_.edge_handle(mesh_.n_edges() - 2) : mesh_.edge_handle(mesh_.n_edges() - 3));
603 
604  f0 = mesh_.face_handle(mesh_.halfedge_handle(e0, 0));
605  f1 = mesh_.face_handle(mesh_.halfedge_handle(e0, 1));
606  f2 = mesh_.face_handle(mesh_.halfedge_handle(e1, 0));
607  f3 = mesh_.face_handle(mesh_.halfedge_handle(e1, 1));
608 
609  if (f0.is_valid() && f3.is_valid()) mesh_.property(pids, f3) = mesh_.property(pids, f0);
610 
611  if (f1.is_valid() && f2.is_valid()) mesh_.property(pids, f2) = mesh_.property(pids, f1);
612  }
613 
614  ok = false;
615  }
616  }
617  }
618 
619  if (i == 100) omlog() << "split break\n";
620 }
621 
622 
623 //-----------------------------------------------------------------------------
624 
625 
626 template <class Mesh>
627 void
630 {
631  typename Mesh::EIter e_it, e_end;
632  typename Mesh::CVVIter vv_it;
633  typename Mesh::VHandle v0, v1;
634  typename Mesh::HHandle h0, h1, h01, h10;
635  typename Mesh::Point p;
636  bool ok, skip, b0, b1, l0, l1, f0, f1;
637  int i;
638  bool hcol01, hcol10;
639 
640  for (ok = false, i = 0; !ok && i < 100; ++i) {
641  ok = true;
642 
643  for (e_it = mesh_.edges_begin(), e_end = mesh_.edges_end(); e_it != e_end; ++e_it) {
644  if (!mesh_.status(*e_it).deleted() && !mesh_.status(*e_it).locked()) {
645  h10 = mesh_.halfedge_handle(*e_it, 0);
646  h01 = mesh_.halfedge_handle(*e_it, 1);
647  v0 = mesh_.to_vertex_handle(h10);
648  v1 = mesh_.to_vertex_handle(h01);
649 
650  if (is_too_short(v0, v1)) {
651  // get status
652  b0 = mesh_.is_boundary(v0);
653  b1 = mesh_.is_boundary(v1);
654  l0 = mesh_.status(v0).locked();
655  l1 = mesh_.status(v1).locked();
656  f0 = mesh_.status(v0).feature();
657  f1 = mesh_.status(v1).feature();
658  hcol01 = hcol10 = true;
659 
660  // boundary rules
661  if (b0 && b1) {
662  if (!mesh_.is_boundary(*e_it))
663  continue;
664  } else if (b0)
665  hcol01 = false;
666  else if (b1)
667  hcol10 = false;
668 
669  // locked rules
670  if (l0 && l1)
671  continue;
672  else if (l0)
673  hcol01 = false;
674  else if (l1)
675  hcol10 = false;
676 
677  // feature vertex rules
678  if (f0 && f1)
679  continue;
680  else if (f0)
681  hcol01 = false;
682  else if (f1)
683  hcol10 = false;
684 
685  // feature edge rules
686  // if vertex is incident to feature edge but collapse edge is not feature dont collapse
687 
688  // hcol01
689  {
690  int feature_valence_0 = 0;
691  for (auto eh : mesh_.ve_range(v0))
692  if (mesh_.status(eh).feature())
693  ++feature_valence_0;
694  if (feature_valence_0 == 2)
695  {
696  if (!mesh_.status(*e_it).feature())
697  hcol01 = false;
698  }
699  else if (feature_valence_0 != 0)
700  hcol01 = false;
701  }
702  // hcol10
703  {
704  int feature_valence_1 = 0;
705  for (auto eh : mesh_.ve_range(v1))
706  if (mesh_.status(eh).feature())
707  ++feature_valence_1;
708  if (feature_valence_1 == 2)
709  {
710  if (!mesh_.status(*e_it).feature())
711  hcol10 = false;
712  }
713  else if (feature_valence_1 != 0)
714  hcol10 = false;
715  }
716 
717  // the other two edges removed by collapse must not be features
718  h0 = mesh_.prev_halfedge_handle(h01);
719  h1 = mesh_.next_halfedge_handle(h10);
720  if (mesh_.status(mesh_.edge_handle(h0)).feature() || mesh_.status(mesh_.edge_handle(h1)).feature())
721  hcol01 = false;
722  h0 = mesh_.prev_halfedge_handle(h10);
723  h1 = mesh_.next_halfedge_handle(h01);
724  if (mesh_.status(mesh_.edge_handle(h0)).feature() || mesh_.status(mesh_.edge_handle(h1)).feature())
725  hcol10 = false;
726 
727  // topological rules
728  if (hcol01)
729  hcol01 = mesh_.is_collapse_ok(h01);
730  if (hcol10)
731  hcol10 = mesh_.is_collapse_ok(h10);
732 
733  // both collapses possible: collapse into vertex w/ higher valence
734  if (hcol01 && hcol10) {
735  if (mesh_.valence(v0) < mesh_.valence(v1))
736  hcol10 = false;
737  else
738  hcol01 = false;
739  }
740 
741  // check if collapse flips triangles
742  if (collapse_flips_normal(h01))
743  hcol01 = false;
744  if (collapse_flips_normal(h10))
745  hcol10 = false;
746 
747  // try v1 -> v0
748  if (hcol10) {
749  // don't create too long edges
750  skip = false;
751  for (vv_it = mesh_.cvv_iter(v1); vv_it.is_valid() && !skip; ++vv_it)
752  if (is_too_long(v0, *vv_it))
753  skip = true;
754 
755  if (!skip) {
756  mesh_.collapse(h10);
757  ok = false;
758  }
759  }
760 
761  // try v0 -> v1
762  else if (hcol01) {
763  // don't create too long edges
764  skip = false;
765  for (vv_it = mesh_.cvv_iter(v0); vv_it.is_valid() && !skip; ++vv_it)
766  if (is_too_long(v1, *vv_it))
767  skip = true;
768 
769  if (!skip) {
770  mesh_.collapse(h01);
771  ok = false;
772  }
773  }
774  }
775  }
776  }
777  }
778 
779  mesh_.garbage_collection();
780 
781  if (i == 100)
782  omlog() << "collapse break\n";
783 }
784 
785 
786 //-----------------------------------------------------------------------------
787 
788 
789 template <class Mesh>
790 void
792 flip_edges()
793 {
794  typename Mesh::EIter e_it, e_end;
795  typename Mesh::VIter v_it, v_end;
796  typename Mesh::VHandle v0, v1, v2, v3, vh;
797  typename Mesh::HHandle hh;
798  typename Mesh::Point p;
799  typename Mesh::FHandle fh;
800 
801  int val0, val1, val2, val3;
802  int val_opt0, val_opt1, val_opt2, val_opt3;
803  int ve0, ve1, ve2, ve3, ve_before, ve_after;
804  bool ok;
805  int i;
806 
807  // compute vertex valences
808  for (v_it = mesh_.vertices_begin(), v_end = mesh_.vertices_end(); v_it != v_end; ++v_it)
809  mesh_.property(valences_, *v_it) = mesh_.valence(*v_it);
810 
811  // flip all edges
812  for (ok = false, i = 0; !ok && i < 100; ++i) {
813  ok = true;
814 
815  for (e_it = mesh_.edges_begin(), e_end = mesh_.edges_end(); e_it != e_end; ++e_it) {
816  if (!mesh_.status(*e_it).locked() && !mesh_.status(*e_it).feature()) {
817  hh = mesh_.halfedge_handle(*e_it, 0);
818  v0 = mesh_.to_vertex_handle(hh);
819  v2 = mesh_.to_vertex_handle(mesh_.next_halfedge_handle(hh));
820  if ( !mesh_.next_halfedge_handle(hh).is_valid() ) {
821  std::cerr << "Error v2" << std::endl;
822  continue;
823  }
824  hh = mesh_.halfedge_handle(*e_it, 1);
825  v1 = mesh_.to_vertex_handle(hh);
826  v3 = mesh_.to_vertex_handle(mesh_.next_halfedge_handle(hh));
827  if ( !mesh_.next_halfedge_handle(hh).is_valid() ) {
828  std::cerr << "Error v3" << std::endl;
829  continue;
830  }
831 
832  if ( !v2.is_valid())
833  continue;
834 
835  if (!mesh_.status(v0).locked() && !mesh_.status(v1).locked() && !mesh_.status(v2).locked()
836  && !mesh_.status(v3).locked()) {
837  val0 = mesh_.property(valences_, v0);
838  val1 = mesh_.property(valences_, v1);
839  val2 = mesh_.property(valences_, v2);
840  val3 = mesh_.property(valences_, v3);
841 
842  val_opt0 = (mesh_.is_boundary(v0) ? 4 : 6);
843  val_opt1 = (mesh_.is_boundary(v1) ? 4 : 6);
844  val_opt2 = (mesh_.is_boundary(v2) ? 4 : 6);
845  val_opt3 = (mesh_.is_boundary(v3) ? 4 : 6);
846 
847  ve0 = (val0 - val_opt0);
848  ve0 *= ve0;
849  ve1 = (val1 - val_opt1);
850  ve1 *= ve1;
851  ve2 = (val2 - val_opt2);
852  ve2 *= ve2;
853  ve3 = (val3 - val_opt3);
854  ve3 *= ve3;
855 
856  ve_before = ve0 + ve1 + ve2 + ve3;
857 
858  --val0;
859  --val1;
860  ++val2;
861  ++val3;
862 
863  ve0 = (val0 - val_opt0);
864  ve0 *= ve0;
865  ve1 = (val1 - val_opt1);
866  ve1 *= ve1;
867  ve2 = (val2 - val_opt2);
868  ve2 *= ve2;
869  ve3 = (val3 - val_opt3);
870  ve3 *= ve3;
871 
872  ve_after = ve0 + ve1 + ve2 + ve3;
873 
874  if (ve_before > ve_after && mesh_.is_flip_ok(*e_it) && !edge_flip_flips_normal(*e_it)) {
875  mesh_.flip(*e_it);
876  --mesh_.property(valences_, v0);
877  --mesh_.property(valences_, v1);
878  ++mesh_.property(valences_, v2);
879  ++mesh_.property(valences_, v3);
880  ok = false;
881  }
882  }
883  }
884  }
885  }
886 
887  if (i == 100)
888  omlog() << "flip break\n";
889 }
890 
891 
892 //-----------------------------------------------------------------------------
893 
894 
895 template <class Mesh>
896 void
898 tangential_smoothing(bool _use_projection)
899 {
900  typename Mesh::VIter v_it, v_end(mesh_.vertices_end());
901  typename Mesh::CVVIter vv_it;
902  typename Mesh::Point u, n;
903 
904 
905  // tag active vertices
906  for (v_it=mesh_.vertices_begin(); v_it!=v_end; ++v_it)
907  mesh_.status(*v_it).set_tagged( !mesh_.status(*v_it).locked() &&
908  !mesh_.status(*v_it).feature() &&
909  !mesh_.is_boundary(*v_it) );
910 
911  mesh_.update_face_normals();
912 
913  double damping = 1.0;
914 
915  // smooth
916  for (int iters=0; iters<10; ++iters)
917  {
918  for (v_it=mesh_.vertices_begin(); v_it!=v_end; ++v_it)
919  {
920  if (mesh_.status(*v_it).tagged())
921  {
922  u.vectorize(0.0);
923  int valence = 0;
924  int feature_valence = 0;
925 
926  for (auto heh : mesh_.voh_range(*v_it))
927  {
928  u += mesh_.point(mesh_.to_vertex_handle(heh));
929  if (mesh_.status(mesh_.edge_handle(heh)).feature())
930  ++feature_valence;
931  ++valence;
932  }
933 
934  if (valence)
935  {
936  u *= (1.0/valence);
937  u -= mesh_.point(*v_it);
938  n = mesh_.normal(*v_it);
939  n *= (u | n);
940  u -= n;
941  }
942 
943  if (feature_valence == 2)
944  {
945  // project update onto feature directions
946  typename Mesh::Point feature_dir(0.0,0.0,0.0);
947  for (auto heh : mesh_.voh_range(*v_it))
948  {
949  if (mesh_.status(mesh_.edge_handle(heh)).feature())
950  {
951  auto dir = mesh_.point(mesh_.to_vertex_handle(heh)) - mesh_.point(mesh_.from_vertex_handle(heh));
952  if ((dir | feature_dir) > 0)
953  feature_dir += dir;
954  else
955  feature_dir -= dir;
956  }
957  }
958  if (feature_dir.sqrnorm() > 0)
959  feature_dir.normalize();
960  u = (feature_dir | u) * feature_dir;
961  }
962  else if (feature_valence != 0)
963  {
964  // more than two or exactly one incident feature edge means vertex should be preserved
965  u.vectorize(0.0);
966  }
967 
968  mesh_.property(update_, *v_it) = u;
969  }
970  }
971 
972  for (v_it=mesh_.vertices_begin(); v_it!=v_end; ++v_it)
973  if (mesh_.status(*v_it).tagged())
974  mesh_.point(*v_it) += damping*mesh_.property(update_, *v_it);
975 
976  // check if normals changed
977  bool normals_changed = false;
978  for (auto fh : mesh_.faces())
979  {
980  if ((mesh_.calc_face_normal(fh) | mesh_.normal(fh)) < 0)
981  {
982  normals_changed = true;
983  break;
984  }
985  }
986 
987  if (normals_changed)
988  {
989  // revert update and try again with smaller step
990  for (v_it=mesh_.vertices_begin(); v_it!=v_end; ++v_it)
991  if (mesh_.status(*v_it).tagged())
992  mesh_.point(*v_it) -= damping*mesh_.property(update_, *v_it);
993  damping *= 0.5;
994  }
995  }
996 
997 
998  // reset tagged status
999  for (v_it=mesh_.vertices_begin(); v_it!=v_end; ++v_it)
1000  mesh_.status(*v_it).set_tagged(false);
1001 
1002 
1003  // project
1004  if (_use_projection)
1005  for (v_it=mesh_.vertices_begin(); v_it!=v_end; ++v_it)
1006  if (!mesh_.status(*v_it).locked() && !mesh_.status(*v_it).feature())
1007  project_to_reference(*v_it);
1008 }
1009 
1010 
1011 //-----------------------------------------------------------------------------
1012 
1013 
1014 template <class Mesh>
1015 void
1017 balanace_area(unsigned int _iters, bool _use_projection)
1018 {
1019  typename Mesh::VIter v_it, v_end(mesh_.vertices_end());
1020  typename Mesh::CVVIter vv_it;
1021  typename Mesh::Scalar w, ww;
1022  typename Mesh::Point u, n;
1023 
1024 
1025  DiffGeoT<Mesh> diffgeo(mesh_);
1026 
1027 
1028  // tag active vertices
1029  for (v_it=mesh_.vertices_begin(); v_it!=v_end; ++v_it)
1030  {
1031  bool active = ( !mesh_.status(*v_it).locked() &&
1032  !mesh_.status(*v_it).feature() &&
1033  !mesh_.is_boundary(*v_it) );
1034 
1035  // don't move neighbors of boundary vertices
1036  for (vv_it=mesh_.cvv_iter(*v_it); active && vv_it.is_valid(); ++vv_it)
1037  if (mesh_.is_boundary(*vv_it))
1038  active = false;
1039 
1040  mesh_.status(*v_it).set_tagged( active );
1041  }
1042 
1043 
1044  // tag2 vertices for which area has to be computed
1045  for (v_it=mesh_.vertices_begin(); v_it!=v_end; ++v_it)
1046  {
1047  mesh_.status(*v_it).set_tagged2(false);
1048  for (vv_it=mesh_.cvv_iter(*v_it); vv_it.is_valid(); ++vv_it)
1049  {
1050  if (mesh_.status(*vv_it).tagged())
1051  {
1052  mesh_.status(*v_it).set_tagged2(true);
1053  break;
1054  }
1055  }
1056  }
1057 
1058 
1059 
1060  for (unsigned int bla=0; bla<_iters; ++bla)
1061  {
1062  // compute vertex areas
1063  for (v_it=mesh_.vertices_begin(); v_it!=v_end; ++v_it)
1064  if (mesh_.status(*v_it).tagged2())
1065  mesh_.property(area_, *v_it) = pow(diffgeo.compute_area(*v_it), 2);
1066 
1067 
1068 
1069  // smooth
1070  for (int iters=0; iters<10; ++iters)
1071  {
1072  for (v_it=mesh_.vertices_begin(); v_it!=v_end; ++v_it)
1073  {
1074  if (mesh_.status(*v_it).tagged())
1075  {
1076  u.vectorize(0.0);
1077  ww = 0;
1078 
1079  for (vv_it=mesh_.cvv_iter(*v_it); vv_it.is_valid(); ++vv_it)
1080  {
1081  w = mesh_.property(area_, *vv_it);
1082  u += mesh_.point(*vv_it) * w;
1083  ww += w;
1084  }
1085 
1086  if (ww)
1087  {
1088  u *= (1.0/ww);
1089  u -= mesh_.point(*v_it);
1090  n = mesh_.normal(*v_it);
1091  n *= (u | n);
1092  u -= n;
1093 
1094  u *= 0.1;
1095  }
1096 
1097  mesh_.property(update_, *v_it) = u;
1098  }
1099  }
1100 
1101  for (v_it=mesh_.vertices_begin(); v_it!=v_end; ++v_it)
1102  if (!mesh_.status(*v_it).locked() &&
1103  !mesh_.status(*v_it).feature() &&
1104  !mesh_.is_boundary(*v_it))
1105  mesh_.point(*v_it) += mesh_.property(update_, *v_it);
1106  }
1107 
1108 
1109  // project
1110  if (_use_projection)
1111  for (v_it=mesh_.vertices_begin(); v_it!=v_end; ++v_it)
1112  if (!mesh_.status(*v_it).locked() && !mesh_.status(*v_it).feature())
1113  project_to_reference(*v_it);
1114 
1115 
1116 // if (emit_progress_)
1117 // if (!Progress().step(3))
1118 // break;
1119  }
1120 }
1121 
1122 
1123 //-----------------------------------------------------------------------------
1124 
1125 
1126 template <class Mesh>
1127 void
1129 remove_caps()
1130 {
1131  typename Mesh::EdgeIter e_it, e_end(mesh_.edges_end());
1132  typename Mesh::HalfedgeHandle h;
1133  typename Mesh::Scalar a0, a1, amin, aa(cos(170.0 * M_PI / 180.0));
1134  typename Mesh::VHandle v, vb, vd;
1135  typename Mesh::FHandle fb, fd;
1136  typename Mesh::Point a, b, c, d;
1137 
1138 
1139  // handle Nastran PIDs for edge flips
1141  mesh_.get_property_handle(pids, "Nastran PIDs");
1142 
1143 
1144  for (e_it=mesh_.edges_begin(); e_it!=e_end; ++e_it)
1145  {
1146  if (!mesh_.status(*e_it).locked() &&
1147  mesh_.is_flip_ok(*e_it))
1148  {
1149  h = mesh_.halfedge_handle(*e_it, 0);
1150  a = mesh_.point(mesh_.to_vertex_handle(h));
1151 
1152  h = mesh_.next_halfedge_handle(h);
1153  b = mesh_.point(vb=mesh_.to_vertex_handle(h));
1154 
1155  h = mesh_.halfedge_handle(*e_it, 1);
1156  c = mesh_.point(mesh_.to_vertex_handle(h));
1157 
1158  h = mesh_.next_halfedge_handle(h);
1159  d = mesh_.point(vd=mesh_.to_vertex_handle(h));
1160 
1161  a0 = ( (a-b).normalize() | (c-b).normalize() );
1162  a1 = ( (a-d).normalize() | (c-d).normalize() );
1163 
1164  if (a0 < a1) { amin = a0; v = vb; }
1165  else { amin = a1; v = vd; }
1166 
1167  // is it a cap?
1168  if (amin < aa)
1169  {
1170  // feature edge and feature vertex -> seems to be intended
1171  if (mesh_.status(*e_it).feature() && mesh_.status(v).feature())
1172  continue;
1173 
1174  // handle PIDs: flip = split + collapse
1175  if (pids.is_valid())
1176  {
1177  fb = mesh_.face_handle(mesh_.halfedge_handle(*e_it, 0));
1178  fd = mesh_.face_handle(mesh_.halfedge_handle(*e_it, 1));
1179 
1180  if (v == vb)
1181  mesh_.property(pids, fb) = mesh_.property(pids, fd);
1182  else
1183  mesh_.property(pids, fd) = mesh_.property(pids, fb);
1184  }
1185 
1186  // project v onto feature edge
1187  if (mesh_.status(*e_it).feature())
1188  mesh_.set_point(v, (a+c)*0.5);
1189 
1190  // flip
1191  mesh_.flip(*e_it);
1192  }
1193  }
1194  }
1195 }
1196 
1197 
1198 template<class Mesh>
1199 bool
1201 edge_flip_flips_normal(EdgeHandle _eh)
1202 {
1203  if (mesh_.is_boundary(_eh))
1204  return true;
1205 
1206  auto heh = mesh_.halfedge_handle(_eh, 0);
1207 
1208  // get the four points of the two incident faces in ccw order
1209  auto p0 = mesh_.point(mesh_.to_vertex_handle(heh));
1210  auto p1 = mesh_.point(mesh_.to_vertex_handle(mesh_.next_halfedge_handle(heh)));
1211  auto p2 = mesh_.point(mesh_.from_vertex_handle(heh));
1212  auto p3 = mesh_.point(mesh_.to_vertex_handle(mesh_.next_halfedge_handle(mesh_.opposite_halfedge_handle(heh))));
1213 
1214  // compute normals before flip
1215  auto n0_before = (p1-p0) % (p2-p0);
1216  auto n1_before = (p2-p0) % (p3-p0);
1217 
1218  // compute normals after flip
1219  auto n0_after = (p1-p0) % (p3-p0);
1220  auto n1_after = (p3-p2) % (p1-p2);
1221 
1222  // compare all pairs of before and after normals
1223  if ((n0_before | n0_after) < 0 || (n1_before | n1_after) < 0 ||
1224  (n0_before | n1_after) < 0 || (n1_before | n0_after) < 0)
1225  return true;
1226 
1227  return false;
1228 }
1229 
1230 
1231 template<class Mesh>
1232 bool
1234 collapse_flips_normal(HalfedgeHandle _heh)
1235 {
1236  if (!mesh_.is_collapse_ok(_heh))
1237  return true;
1238 
1239  // get faces that are removed by the collapse
1240  auto fh0 = mesh_.face_handle(_heh);
1241  auto fh1 = mesh_.face_handle(mesh_.opposite_halfedge_handle(_heh));
1242 
1243  auto collapsing_vertex = mesh_.from_vertex_handle(_heh);
1244  auto point_before = mesh_.point(collapsing_vertex);
1245 
1246  // compute normals before collapse
1247  std::vector<decltype (mesh_.calc_face_normal(OpenMesh::FaceHandle(0)))> normals_before;
1248  for (auto fh : mesh_.vf_range(collapsing_vertex))
1249  normals_before.push_back(mesh_.calc_face_normal(fh));
1250 
1251  // move point
1252  mesh_.point(collapsing_vertex) = mesh_.point(mesh_.to_vertex_handle(_heh));
1253 
1254  bool collapse_ok = true;
1255 
1256  int i = 0;
1257  for (auto fh : mesh_.vf_range(collapsing_vertex))
1258  {
1259  if (fh != fh0 && fh != fh1) // we don't care about faces that are removec by the collapse
1260  {
1261  auto normal_after = mesh_.calc_face_normal(fh);
1262  if ((normal_after | normals_before[i]) <= 0)
1263  collapse_ok = false;
1264  }
1265  ++i;
1266  }
1267 
1268  // move point back
1269  mesh_.point(collapsing_vertex) = point_before;
1270 
1271  return !collapse_ok;
1272 }
1273 
1274 
1275 //=============================================================================
1276 } // namespace Remeshing
1277 //=============================================================================
Handle for a face entity.
Definition: Handles.hh:141
Kernel::ConstVertexFaceIter ConstVertexFaceIter
Circulator.
Definition: PolyMeshT.hh:176
Kernel::Normal Normal
Normal type.
Definition: PolyMeshT.hh:114
Kernel::Point Point
Coordinate type.
Definition: PolyMeshT.hh:112
void clearVertexSelection(MeshT *_mesh)
Set all vertices to unselected.
void build(unsigned int _max_handles, unsigned int _max_depth)
bool is_valid() const
The handle is valid iff the index is not negative.
Definition: Handles.hh:72
void push_back(Handle _h)
Add a handle to the BSP.
void prepare_face_selection()
prepare for remeshing only vertices which are fully surrounded by selected faces (if no face was sele...
bool baryCoord(const VectorT< Scalar, 3 > &_p, const VectorT< Scalar, 3 > &_u, const VectorT< Scalar, 3 > &_v, const VectorT< Scalar, 3 > &_w, VectorT< Scalar, 3 > &_result)
Definition: Algorithms.cc:83
Normal calc_face_normal(FaceHandle _fh) const
Functions for selection on a mesh.
SmartVertexHandle add_vertex(const Point &_p)
Definition: PolyMeshT.hh:245
Vec::value_type distPointTriangle(const Vec &_p, const Vec &_v0, const Vec &_v1, const Vec &_v2, Vec &_nearestPoint)
distance from point _p to triangle (_v0, _v1, _v2)
Definition: Algorithms.hh:326
SmartVertexHandle split(EdgeHandle _eh, const Point &_p)
Edge split (= 2-to-4 split)
Definition: TriMeshT.hh:275
void reserve(size_t _n)
Reserve memory for _n entries.
void update_normals()
Compute normals for all primitives.
NearestNeighbor nearest(const Point &_p) const
Return handle of the nearest neighbor face.
Kernel::Scalar Scalar
Scalar type.
Definition: PolyMeshT.hh:110
void prepare_vertex_selection()
prepare for remeshing only selected vertices (if no vertex was selected, remesh whole mesh) ...
void update_face_normals()
Update normal vectors for all faces.
void growVertexSelection(MeshT *_mesh)
Grow vertex selection.
Kernel::VertexHandle VertexHandle
Handle for referencing the corresponding item.
Definition: PolyMeshT.hh:136