Developer Documentation
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
BaseRemesherT.cc
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 * $Revision$ *
45 * $LastChangedBy$ *
46 * $Date$ *
47 * *
48 \*===========================================================================*/
49 
50 //=============================================================================
51 //
52 // CLASS BaseRemesherT - IMPLEMENTATION
53 //
54 //=============================================================================
55 
56 #define BASE_REMESHERT_C
57 
58 
59 //== INCLUDES =================================================================
60 
61 
62 #include "BaseRemesherT.hh"
63 
64 // Remshing
65 #include "DiffGeoT.hh"
66 
67 // OpenMesh
68 #include <OpenMesh/Core/Utils/Property.hh>
70 
71 // Selection Stuff
73 
74 
75 //== NAMESPACES ===============================================================
76 
77 namespace Remeshing {
78 
79 
80 //== IMPLEMENTATION ==========================================================
81 
82 
83 template <class Mesh>
84 BaseRemesherT<Mesh>::
85 BaseRemesherT(Mesh& _mesh, ProgressEmitter* _progress)
86  : mesh_(_mesh), refmesh_(0), bsp_(0), progress_(_progress) {
87 
88 }
89 
90 
91 //-----------------------------------------------------------------------------
92 
93 
94 template <class Mesh>
95 BaseRemesherT<Mesh>::
96 ~BaseRemesherT() {
97 
98  delete bsp_;
99  delete refmesh_;
100 }
101 
102 
103 //-----------------------------------------------------------------------------
104 
105 
106 template <class Mesh>
107 void
108 BaseRemesherT<Mesh>::
109 init_reference()
110 {
111  typename Mesh::VIter v_it, v_end;
112  typename Mesh::FIter f_it, f_end;
113  typename Mesh::CFVIter fv_it;
114  typename Mesh::VHandle v0, v1, v2;
115 
116 
117  // clear old stuff first
118  if (bsp_) delete_reference();
119 
120  // tag non-locked vertices + one ring
121  for (v_it=mesh_.vertices_begin(), v_end=mesh_.vertices_end();
122  v_it!=v_end; ++v_it)
123  mesh_.status(*v_it).set_tagged(!mesh_.status(*v_it).locked());
124 
125  MeshSelection::growVertexSelection(&mesh_);
126 
127  // create reference mesh
128  refmesh_ = new Mesh();
129 
131  mesh_.add_property(vhandles);
132 
133  for (v_it=mesh_.vertices_begin(), v_end=mesh_.vertices_end();
134  v_it!=v_end; ++v_it)
135  if (mesh_.status(*v_it).tagged())
136  mesh_.property(vhandles, *v_it) = refmesh_->add_vertex(mesh_.point(*v_it));
137 
138  for (f_it=mesh_.faces_begin(), f_end=mesh_.faces_end();
139  f_it!=f_end; ++f_it)
140  {
141  fv_it = mesh_.cfv_iter(*f_it);
142  v0 = *fv_it;
143  v1 = *(++fv_it);
144  v2 = *(++fv_it);
145 
146  if (mesh_.status(v0).tagged() &&
147  mesh_.status(v1).tagged() &&
148  mesh_.status(v2).tagged())
149  {
150  refmesh_->add_face( mesh_.property(vhandles, v0),
151  mesh_.property(vhandles, v1),
152  mesh_.property(vhandles, v2) );
153  }
154  }
155 
156  MeshSelection::clearVertexSelection(&mesh_);
157  mesh_.remove_property(vhandles);
158 
159  // need vertex normals
160  refmesh_->request_face_normals();
161  refmesh_->request_vertex_normals();
162  refmesh_->update_normals(/*Mesh::ANGLE_WEIGHTED*/);
163 
164 
165  // build BSP
166  bsp_ = new BSP(*refmesh_);
167  bsp_->reserve(refmesh_->n_faces());
168  for (f_it=refmesh_->faces_begin(), f_end=refmesh_->faces_end();
169  f_it!=f_end; ++f_it)
170  bsp_->push_back(*f_it);
171  bsp_->build(10, 100);
172 }
173 
174 
175 //-----------------------------------------------------------------------------
176 
177 
178 template <class Mesh>
179 void
180 BaseRemesherT<Mesh>::
181 delete_reference()
182 {
183  delete bsp_; bsp_ = 0;
184  delete refmesh_; refmesh_ = 0;
185 }
186 
187 
188 //-----------------------------------------------------------------------------
189 
190 
191 template <class Mesh>
192 void
193 BaseRemesherT<Mesh>::
194 project_to_reference(typename Mesh::VertexHandle _vh) const
195 {
196  typename Mesh::Point p = mesh_.point(_vh);
197  typename Mesh::FaceHandle fh = bsp_->nearest(p).handle;
198 
199  if ( ! fh.is_valid() )
200  return;
201 
202  typename Mesh::CFVIter fv_it = refmesh_->cfv_iter(fh);
203 
204 
205  const typename Mesh::Point& p0 = refmesh_->point(*fv_it);
206  typename Mesh::Normal n0 = refmesh_->normal(*fv_it);
207  const typename Mesh::Point& p1 = refmesh_->point(*(++fv_it));
208  typename Mesh::Normal n1 = refmesh_->normal(*fv_it);
209  const typename Mesh::Point& p2 = refmesh_->point(*(++fv_it));
210  typename Mesh::Normal n2 = refmesh_->normal(*fv_it);
211 
212 
213  // project
214  ACG::Geometry::distPointTriangle(p, p0, p1, p2, p);
215  mesh_.set_point(_vh, p);
216 
217 
218  // get barycentric coordinates
219  if (!ACG::Geometry::baryCoord(p, p0, p1, p2, p))
220  p[0] = p[1] = p[2] = 1.0/3.0;
221 
222 
223  // interpolate normal
224  typename Mesh::Normal n;
225  n = (n0 *= p[0]);
226  n += (n1 *= p[1]);
227  n += (n2 *= p[2]);
228  n.normalize();
229  mesh_.set_normal(_vh, n);
230 }
231 
232 
233 //-----------------------------------------------------------------------------
234 
235 
236 template <class Mesh>
237 void
238 BaseRemesherT<Mesh>::
239 remesh(unsigned int _iters,
240  unsigned int _area_iters,
241  bool _use_projection,
242  Selection _selection) {
243 
244  try
245  {
246  if (_selection == VERTEX_SELECTION)
247  prepare_vertex_selection();
248  else if (_selection == FACE_SELECTION)
249  prepare_face_selection();
250  remeshh(_iters, _area_iters, _use_projection);
251  }
252  catch (std::bad_alloc&)
253  {
254  mesh_.clear();
255  omerr() << "Remeshig: Out of memory\n";
256  }
257 
258  cleanup();
259 }
260 
261 //-----------------------------------------------------------------------------
262 
263 template <class Mesh>
264 void
267 {
268  typename Mesh::EIter e_it, e_end;
269  typename Mesh::VIter v_it, v_end;
270  typename Mesh::FIter f_it, f_end;
271  typename Mesh::CFVIter fv_it;
272  typename Mesh::CVOHIter vh_it;
273  typename Mesh::VHandle v0, v1;
274  typename Mesh::FHandle f0, f1;
275 
276 
277  // need vertex and edge status
278  mesh_.request_vertex_status();
279  mesh_.request_edge_status();
280  mesh_.request_face_status();
281 
282 
283 
284  // if nothing selected -> select all
285  nothing_selected_ = true;
286  for (v_it=mesh_.vertices_begin(), v_end=mesh_.vertices_end();
287  v_it!=v_end; ++v_it)
288  if (mesh_.status(*v_it).selected())
289  { nothing_selected_ = false; break; }
290 
291  if (nothing_selected_)
292  for (v_it=mesh_.vertices_begin(), v_end=mesh_.vertices_end();
293  v_it!=v_end; ++v_it)
294  mesh_.status(*v_it).set_selected(true);
295 
296 
297 
298  // lock un-selected vertices & edges
299  for (v_it=mesh_.vertices_begin(), v_end=mesh_.vertices_end();
300  v_it!=v_end; ++v_it)
301  mesh_.status(*v_it).set_locked(!mesh_.status(*v_it).selected());
302 
303  for (e_it=mesh_.edges_begin(), e_end=mesh_.edges_end();
304  e_it!=e_end; ++e_it)
305  {
306  v0 = mesh_.to_vertex_handle(mesh_.halfedge_handle(*e_it, 0));
307  v1 = mesh_.to_vertex_handle(mesh_.halfedge_handle(*e_it, 1));
308  mesh_.status(*e_it).set_locked(mesh_.status(v0).locked() ||
309  mesh_.status(v1).locked());
310  }
311 
312 
313 
314  // handle feature corners:
315  // lock corner vertices (>2 feature edges) and endpoints
316  // of feature lines (1 feature edge)
317  for (v_it=mesh_.vertices_begin(), v_end=mesh_.vertices_end();
318  v_it!=v_end; ++v_it)
319  {
320  if (mesh_.status(*v_it).feature())
321  {
322  int c=0;
323  for (vh_it=mesh_.cvoh_iter(*v_it); vh_it.is_valid(); ++vh_it)
324  if (mesh_.status(mesh_.edge_handle(*vh_it)).feature())
325  ++c;
326  if (c!=2) mesh_.status(*v_it).set_locked(true);
327  }
328  }
329 
330 
331 
332  // build reference mesh
333  init_reference();
334 // if (emit_progress_) Progress().step(5);
335 
336 
337  // add properties
338  mesh_.add_property(valences_);
339  mesh_.add_property(update_);
340  mesh_.add_property(area_);
341 }
342 
343 //-----------------------------------------------------------------------------
344 
345 
346 template <class Mesh>
347 void
350 {
351  typename Mesh::EIter e_it, e_end;
352  typename Mesh::VIter v_it, v_end;
353  typename Mesh::FIter f_it, f_end;
354  typename Mesh::CFVIter fv_it;
355  typename Mesh::CVOHIter vh_it;
356  typename Mesh::VHandle v0, v1;
357  typename Mesh::FHandle f0, f1;
358 
359 
360  // need vertex and edge status
361  mesh_.request_vertex_status();
362  mesh_.request_edge_status();
363  mesh_.request_face_status();
364 
365  // if nothing selected -> select all
366  nothing_selected_ = true;
367  for (f_it = mesh_.faces_begin(), f_end = mesh_.faces_end(); f_it != f_end;
368  ++f_it)
369  {
370  if (mesh_.status(*f_it).selected())
371  {
372  nothing_selected_ = false;
373  break;
374  }
375  }
376 
377  if (nothing_selected_)
378  MeshSelection::selectAllFaces(&mesh_);
379 
380 
381 
382  // lock un-selected vertices & edges
383  for (v_it = mesh_.vertices_begin(), v_end = mesh_.vertices_end();
384  v_it != v_end; ++v_it)
385  {
386  bool all_faces_selected = true;
387 
388  for (typename Mesh::ConstVertexFaceIter vf_it = mesh_.cvf_iter(*v_it);
389  vf_it.is_valid(); ++vf_it)
390  {
391  if (!mesh_.status(*vf_it).selected())
392  {
393  all_faces_selected = false;
394  break;
395  }
396  }
397  mesh_.status(*v_it).set_locked(!all_faces_selected);
398  }
399 
400  for (e_it = mesh_.edges_begin(), e_end = mesh_.edges_end();
401  e_it != e_end; ++e_it)
402  {
403  if (mesh_.is_boundary(*e_it))
404  {
405  mesh_.status(*e_it).set_locked(true);
406  }
407  else
408  {
409  f0 = mesh_.face_handle(mesh_.halfedge_handle(*e_it, 0));
410  f1 = mesh_.face_handle(mesh_.halfedge_handle(*e_it, 1));
411 
412  mesh_.status(*e_it).set_locked(!(mesh_.status(f0).selected() && mesh_.status(f1).selected()));
413  }
414  }
415 
416  MeshSelection::clearFaceSelection(&mesh_);
417 
418  // handle feature corners:
419  // lock corner vertices (>2 feature edges) and endpoints
420  // of feature lines (1 feature edge)
421  for (v_it=mesh_.vertices_begin(), v_end=mesh_.vertices_end();
422  v_it!=v_end; ++v_it)
423  {
424  if (mesh_.status(*v_it).feature())
425  {
426  int c=0;
427  for (vh_it=mesh_.cvoh_iter(*v_it); vh_it.is_valid(); ++vh_it)
428  if (mesh_.status(mesh_.edge_handle(*vh_it)).feature())
429  ++c;
430  if (c!=2) mesh_.status(*v_it).set_locked(true);
431  }
432  }
433 
434 
435 
436  // build reference mesh
437  init_reference();
438 // if (emit_progress_) Progress().step(5);
439 
440 
441  // add properties
442  mesh_.add_property(valences_);
443  mesh_.add_property(update_);
444  mesh_.add_property(area_);
445 }
446 
447 
448 //-----------------------------------------------------------------------------
449 
450 
451 template <class Mesh>
452 void
454 remeshh(unsigned int _iters,
455  unsigned int _area_iters, bool _use_projection) {
456 
457  double progress_step = 100.0 / (((double)_iters/4.0 + (double)_area_iters));
458  double total_progress = 0.0;
459 
460  for (unsigned int i = 0; i < _iters; ++i) {
461 
462  split_long_edges();
463  if(progress_ != NULL) {
464  total_progress += progress_step;
465  progress_->sendProgressSignal(total_progress);
466  }
467 
468  collapse_short_edges();
469  if(progress_ != NULL) {
470  total_progress += progress_step;
471  progress_->sendProgressSignal(total_progress);
472  }
473 
474  flip_edges();
475  if(progress_ != NULL) {
476  total_progress += progress_step;
477  progress_->sendProgressSignal(total_progress);
478  }
479 
480  tangential_smoothing(_use_projection);
481  if(progress_ != NULL) {
482  total_progress += progress_step;
483  progress_->sendProgressSignal(total_progress);
484  }
485  }
486 
487  if (_area_iters) {
488  balanace_area(_area_iters, _use_projection);
489  if(progress_ != NULL) {
490  total_progress += progress_step;
491  progress_->sendProgressSignal(total_progress);
492  }
493  }
494 
495  // feature edges block flips -> might lead to caps
496  remove_caps();
497 }
498 
499 
500 //-----------------------------------------------------------------------------
501 
502 
503 template <class Mesh>
504 void
505 BaseRemesherT<Mesh>::
506 cleanup()
507 {
508  typename Mesh::EIter e_it, e_end;
509  typename Mesh::VIter v_it, v_end;
510  typename Mesh::FIter f_it, f_end;
511  typename Mesh::CFVIter fv_it;
512  typename Mesh::CVOHIter vh_it;
513  typename Mesh::VHandle v0, v1;
514  typename Mesh::FHandle f0, f1;
515 
516 
517  // remove properties
518  mesh_.remove_property(valences_);
519  mesh_.remove_property(update_);
520  mesh_.remove_property(area_);
521 
522 
523  // remove reference again
524  delete_reference();
525 
526 
527  // unlock all vertices & edges
528  for (v_it=mesh_.vertices_begin(), v_end=mesh_.vertices_end();
529  v_it!=v_end; ++v_it)
530  {
531  mesh_.status(*v_it).set_locked(false);
532  }
533  for (e_it=mesh_.edges_begin(), e_end=mesh_.edges_end();
534  e_it!=e_end; ++e_it)
535  {
536  mesh_.status(*e_it).set_locked(false);
537  }
538 
539 
540  // de-select if nothing was selected before
541  if (nothing_selected_)
542  for (v_it=mesh_.vertices_begin(), v_end=mesh_.vertices_end();
543  v_it!=v_end; ++v_it)
544  mesh_.status(*v_it).set_selected(false);
545 
546 
547 
548  // free vertex and edge status
549  mesh_.release_vertex_status();
550  mesh_.release_edge_status();
551  mesh_.release_face_status();
552 }
553 
554 
555 //-----------------------------------------------------------------------------
556 
557 
558 template <class Mesh>
559 void
560 BaseRemesherT<Mesh>::
561 split_long_edges()
562 {
563  typename Mesh::VIter v_it, v_end;
564  typename Mesh::EIter e_it, e_end;
565  typename Mesh::VHandle v0, v1, vh;
566  typename Mesh::EHandle eh, e0, e1;
567  typename Mesh::FHandle f0, f1, f2, f3;
568 
569  bool ok, is_feature, is_boundary;
570  int i;
571 
572  // un-tagg
573  for (v_it = mesh_.vertices_begin(), v_end = mesh_.vertices_end(); v_it != v_end; ++v_it)
574  mesh_.status(*v_it).set_tagged(false);
575 
576  // handle Nastran PIDs during refinement
578  mesh_.get_property_handle(pids, "Nastran PIDs");
579 
580  // split long edges
581  for (ok = false, i = 0; !ok && i < 100; ++i) {
582  ok = true;
583 
584  for (e_it = mesh_.edges_begin(), e_end = mesh_.edges_end(); e_it != e_end; ++e_it) {
585  v0 = mesh_.to_vertex_handle(mesh_.halfedge_handle(*e_it, 0));
586  v1 = mesh_.to_vertex_handle(mesh_.halfedge_handle(*e_it, 1));
587 
588  if (!mesh_.status(*e_it).locked() && is_too_long(v0, v1)) {
589  const typename Mesh::Point& p0 = mesh_.point(v0);
590  const typename Mesh::Point& p1 = mesh_.point(v1);
591 
592  is_feature = mesh_.status(*e_it).feature();
593  is_boundary = mesh_.is_boundary(*e_it);
594 
595  vh = mesh_.add_vertex((p0 + p1) * 0.5);
596  mesh_.split(*e_it, vh);
597 
598  mesh_.status(vh).set_selected(true);
599 
600  if (is_feature) {
601  eh = (is_boundary ? mesh_.edge_handle(mesh_.n_edges() - 2) : mesh_.edge_handle(mesh_.n_edges() - 3));
602 
603  mesh_.status(eh).set_feature(true);
604  mesh_.status(vh).set_feature(true);
605  } else {
606  project_to_reference(vh);
607  }
608 
609  if (pids.is_valid()) {
610  e0 = *e_it;
611  e1 = (is_boundary ? mesh_.edge_handle(mesh_.n_edges() - 2) : mesh_.edge_handle(mesh_.n_edges() - 3));
612 
613  f0 = mesh_.face_handle(mesh_.halfedge_handle(e0, 0));
614  f1 = mesh_.face_handle(mesh_.halfedge_handle(e0, 1));
615  f2 = mesh_.face_handle(mesh_.halfedge_handle(e1, 0));
616  f3 = mesh_.face_handle(mesh_.halfedge_handle(e1, 1));
617 
618  if (f0.is_valid() && f3.is_valid()) mesh_.property(pids, f3) = mesh_.property(pids, f0);
619 
620  if (f1.is_valid() && f2.is_valid()) mesh_.property(pids, f2) = mesh_.property(pids, f1);
621  }
622 
623  ok = false;
624  }
625  }
626  }
627 
628  if (i == 100) omlog() << "split break\n";
629 }
630 
631 
632 //-----------------------------------------------------------------------------
633 
634 
635 template <class Mesh>
636 void
637 BaseRemesherT<Mesh>::
638 collapse_short_edges()
639 {
640  typename Mesh::EIter e_it, e_end;
641  typename Mesh::CVVIter vv_it;
642  typename Mesh::VHandle v0, v1;
643  typename Mesh::HHandle h0, h1, h01, h10;
644  typename Mesh::Point p;
645  bool ok, skip, b0, b1, l0, l1, f0, f1;
646  int i;
647  bool hcol01, hcol10;
648 
649  for (ok = false, i = 0; !ok && i < 100; ++i) {
650  ok = true;
651 
652  for (e_it = mesh_.edges_begin(), e_end = mesh_.edges_end(); e_it != e_end; ++e_it) {
653  if (!mesh_.status(*e_it).deleted() && !mesh_.status(*e_it).locked()) {
654  h10 = mesh_.halfedge_handle(*e_it, 0);
655  h01 = mesh_.halfedge_handle(*e_it, 1);
656  v0 = mesh_.to_vertex_handle(h10);
657  v1 = mesh_.to_vertex_handle(h01);
658 
659  if (is_too_short(v0, v1)) {
660  // get status
661  b0 = mesh_.is_boundary(v0);
662  b1 = mesh_.is_boundary(v1);
663  l0 = mesh_.status(v0).locked();
664  l1 = mesh_.status(v1).locked();
665  f0 = mesh_.status(v0).feature();
666  f1 = mesh_.status(v1).feature();
667  hcol01 = hcol10 = true;
668 
669  if (mesh_.status(*e_it).feature() && !f0 && !f1)
670  std::cerr << "Bad luck" << std::endl;
671 
672  // boundary rules
673  if (b0 && b1) {
674  if (!mesh_.is_boundary(*e_it))
675  continue;
676  } else if (b0)
677  hcol01 = false;
678  else if (b1)
679  hcol10 = false;
680 
681  // locked rules
682  if (l0 && l1)
683  continue;
684  else if (l0)
685  hcol01 = false;
686  else if (l1)
687  hcol10 = false;
688 
689  // feature rules
690 
691  // the other two edges removed by collapse must not be features
692  h0 = mesh_.prev_halfedge_handle(h01);
693  h1 = mesh_.next_halfedge_handle(h10);
694  if (mesh_.status(mesh_.edge_handle(h0)).feature() || mesh_.status(mesh_.edge_handle(h1)).feature())
695  hcol01 = false;
696  h0 = mesh_.prev_halfedge_handle(h10);
697  h1 = mesh_.next_halfedge_handle(h01);
698  if (mesh_.status(mesh_.edge_handle(h0)).feature() || mesh_.status(mesh_.edge_handle(h1)).feature())
699  hcol10 = false;
700 
701  if (f0 && f1) {
702  // edge must be feature
703  if (!mesh_.status(*e_it).feature())
704  continue;
705 
706  // the other two edges removed by collapse must not be features
707  h0 = mesh_.prev_halfedge_handle(h01);
708  h1 = mesh_.next_halfedge_handle(h10);
709  if (mesh_.status(mesh_.edge_handle(h0)).feature() || mesh_.status(mesh_.edge_handle(h1)).feature())
710  hcol01 = false;
711  h0 = mesh_.prev_halfedge_handle(h10);
712  h1 = mesh_.next_halfedge_handle(h01);
713  if (mesh_.status(mesh_.edge_handle(h0)).feature() || mesh_.status(mesh_.edge_handle(h1)).feature())
714  hcol10 = false;
715  } else if (f0)
716  hcol01 = false;
717  else if (f1)
718  hcol10 = false;
719 
720  // topological rules
721  if (hcol01)
722  hcol01 = mesh_.is_collapse_ok(h01);
723  if (hcol10)
724  hcol10 = mesh_.is_collapse_ok(h10);
725 
726  // both collapses possible: collapse into vertex w/ higher valence
727  if (hcol01 && hcol10) {
728  if (mesh_.valence(v0) < mesh_.valence(v1))
729  hcol10 = false;
730  else
731  hcol01 = false;
732  }
733 
734  // try v1 -> v0
735  if (hcol10) {
736  // don't create too long edges
737  skip = false;
738  for (vv_it = mesh_.cvv_iter(v1); vv_it.is_valid() && !skip; ++vv_it)
739  if (is_too_long(v0, *vv_it))
740  skip = true;
741 
742  if (!skip) {
743  mesh_.collapse(h10);
744  ok = false;
745  }
746  }
747 
748  // try v0 -> v1
749  else if (hcol01) {
750  // don't create too long edges
751  skip = false;
752  for (vv_it = mesh_.cvv_iter(v0); vv_it.is_valid() && !skip; ++vv_it)
753  if (is_too_long(v1, *vv_it))
754  skip = true;
755 
756  if (!skip) {
757  mesh_.collapse(h01);
758  ok = false;
759  }
760  }
761  }
762  }
763  }
764  }
765 
766  mesh_.garbage_collection();
767 
768  if (i == 100)
769  omlog() << "collapse break\n";
770 }
771 
772 
773 //-----------------------------------------------------------------------------
774 
775 
776 template <class Mesh>
777 void
778 BaseRemesherT<Mesh>::
779 flip_edges()
780 {
781  typename Mesh::EIter e_it, e_end;
782  typename Mesh::VIter v_it, v_end;
783  typename Mesh::VHandle v0, v1, v2, v3, vh;
784  typename Mesh::HHandle hh;
785  typename Mesh::Point p;
786  typename Mesh::FHandle fh;
787 
788  int val0, val1, val2, val3;
789  int val_opt0, val_opt1, val_opt2, val_opt3;
790  int ve0, ve1, ve2, ve3, ve_before, ve_after;
791  bool ok;
792  int i;
793 
794  // compute vertex valences
795  for (v_it = mesh_.vertices_begin(), v_end = mesh_.vertices_end(); v_it != v_end; ++v_it)
796  mesh_.property(valences_, *v_it) = mesh_.valence(*v_it);
797 
798  // flip all edges
799  for (ok = false, i = 0; !ok && i < 100; ++i) {
800  ok = true;
801 
802  for (e_it = mesh_.edges_begin(), e_end = mesh_.edges_end(); e_it != e_end; ++e_it) {
803  if (!mesh_.status(*e_it).locked() && !mesh_.status(*e_it).feature()) {
804  hh = mesh_.halfedge_handle(*e_it, 0);
805  v0 = mesh_.to_vertex_handle(hh);
806  v2 = mesh_.to_vertex_handle(mesh_.next_halfedge_handle(hh));
807  if ( !mesh_.next_halfedge_handle(hh).is_valid() ) {
808  std::cerr << "Error v2" << std::endl;
809  continue;
810  }
811  hh = mesh_.halfedge_handle(*e_it, 1);
812  v1 = mesh_.to_vertex_handle(hh);
813  v3 = mesh_.to_vertex_handle(mesh_.next_halfedge_handle(hh));
814  if ( !mesh_.next_halfedge_handle(hh).is_valid() ) {
815  std::cerr << "Error v3" << std::endl;
816  continue;
817  }
818 
819  if ( !v2.is_valid())
820  continue;
821 
822  if (!mesh_.status(v0).locked() && !mesh_.status(v1).locked() && !mesh_.status(v2).locked()
823  && !mesh_.status(v3).locked()) {
824  val0 = mesh_.property(valences_, v0);
825  val1 = mesh_.property(valences_, v1);
826  val2 = mesh_.property(valences_, v2);
827  val3 = mesh_.property(valences_, v3);
828 
829  val_opt0 = (mesh_.is_boundary(v0) ? 4 : 6);
830  val_opt1 = (mesh_.is_boundary(v1) ? 4 : 6);
831  val_opt2 = (mesh_.is_boundary(v2) ? 4 : 6);
832  val_opt3 = (mesh_.is_boundary(v3) ? 4 : 6);
833 
834  ve0 = (val0 - val_opt0);
835  ve0 *= ve0;
836  ve1 = (val1 - val_opt1);
837  ve1 *= ve1;
838  ve2 = (val2 - val_opt2);
839  ve2 *= ve2;
840  ve3 = (val3 - val_opt3);
841  ve3 *= ve3;
842 
843  ve_before = ve0 + ve1 + ve2 + ve3;
844 
845  --val0;
846  --val1;
847  ++val2;
848  ++val3;
849 
850  ve0 = (val0 - val_opt0);
851  ve0 *= ve0;
852  ve1 = (val1 - val_opt1);
853  ve1 *= ve1;
854  ve2 = (val2 - val_opt2);
855  ve2 *= ve2;
856  ve3 = (val3 - val_opt3);
857  ve3 *= ve3;
858 
859  ve_after = ve0 + ve1 + ve2 + ve3;
860 
861  if (ve_before > ve_after && mesh_.is_flip_ok(*e_it)) {
862  mesh_.flip(*e_it);
863  --mesh_.property(valences_, v0);
864  --mesh_.property(valences_, v1);
865  ++mesh_.property(valences_, v2);
866  ++mesh_.property(valences_, v3);
867  ok = false;
868  }
869  }
870  }
871  }
872  }
873 
874  if (i == 100)
875  omlog() << "flip break\n";
876 }
877 
878 
879 //-----------------------------------------------------------------------------
880 
881 
882 template <class Mesh>
883 void
884 BaseRemesherT<Mesh>::
885 tangential_smoothing(bool _use_projection)
886 {
887  typename Mesh::VIter v_it, v_end(mesh_.vertices_end());
888  typename Mesh::CVVIter vv_it;
889  typename Mesh::Scalar valence;
890  typename Mesh::Point u, n;
891 
892 
893  // tag active vertices
894  for (v_it=mesh_.vertices_begin(); v_it!=v_end; ++v_it)
895  mesh_.status(*v_it).set_tagged( !mesh_.status(*v_it).locked() &&
896  !mesh_.status(*v_it).feature() &&
897  !mesh_.is_boundary(*v_it) );
898 
899 
900  // smooth
901  for (int iters=0; iters<10; ++iters)
902  {
903  for (v_it=mesh_.vertices_begin(); v_it!=v_end; ++v_it)
904  {
905  if (mesh_.status(*v_it).tagged())
906  {
907  u.vectorize(0.0);
908  valence = 0;
909 
910  for (vv_it=mesh_.cvv_iter(*v_it); vv_it.is_valid(); ++vv_it)
911  {
912  u += mesh_.point(*vv_it);
913  ++valence;
914  }
915 
916  if (valence)
917  {
918  u *= (1.0/valence);
919  u -= mesh_.point(*v_it);
920  n = mesh_.normal(*v_it);
921  n *= (u | n);
922  u -= n;
923  }
924 
925  mesh_.property(update_, *v_it) = u;
926  }
927  }
928 
929  for (v_it=mesh_.vertices_begin(); v_it!=v_end; ++v_it)
930  if (mesh_.status(*v_it).tagged())
931  mesh_.point(*v_it) += mesh_.property(update_, *v_it);
932  }
933 
934 
935  // reset tagged status
936  for (v_it=mesh_.vertices_begin(); v_it!=v_end; ++v_it)
937  mesh_.status(*v_it).set_tagged(false);
938 
939 
940  // project
941  if (_use_projection)
942  for (v_it=mesh_.vertices_begin(); v_it!=v_end; ++v_it)
943  if (!mesh_.status(*v_it).locked() && !mesh_.status(*v_it).feature())
944  project_to_reference(*v_it);
945 }
946 
947 
948 //-----------------------------------------------------------------------------
949 
950 
951 template <class Mesh>
952 void
953 BaseRemesherT<Mesh>::
954 balanace_area(unsigned int _iters, bool _use_projection)
955 {
956  typename Mesh::VIter v_it, v_end(mesh_.vertices_end());
957  typename Mesh::CVVIter vv_it;
958  typename Mesh::Scalar w, ww;
959  typename Mesh::Point u, n;
960 
961 
962  DiffGeoT<Mesh> diffgeo(mesh_);
963 
964 
965  // tag active vertices
966  for (v_it=mesh_.vertices_begin(); v_it!=v_end; ++v_it)
967  {
968  bool active = ( !mesh_.status(*v_it).locked() &&
969  !mesh_.status(*v_it).feature() &&
970  !mesh_.is_boundary(*v_it) );
971 
972  // don't move neighbors of boundary vertices
973  for (vv_it=mesh_.cvv_iter(*v_it); active && vv_it.is_valid(); ++vv_it)
974  if (mesh_.is_boundary(*vv_it))
975  active = false;
976 
977  mesh_.status(*v_it).set_tagged( active );
978  }
979 
980 
981  // tag2 vertices for which area has to be computed
982  for (v_it=mesh_.vertices_begin(); v_it!=v_end; ++v_it)
983  {
984  mesh_.status(*v_it).set_tagged2(false);
985  for (vv_it=mesh_.cvv_iter(*v_it); vv_it.is_valid(); ++vv_it)
986  {
987  if (mesh_.status(*vv_it).tagged())
988  {
989  mesh_.status(*v_it).set_tagged2(true);
990  break;
991  }
992  }
993  }
994 
995 
996 
997  for (unsigned int bla=0; bla<_iters; ++bla)
998  {
999  // compute vertex areas
1000  for (v_it=mesh_.vertices_begin(); v_it!=v_end; ++v_it)
1001  if (mesh_.status(*v_it).tagged2())
1002  mesh_.property(area_, *v_it) = pow(diffgeo.compute_area(*v_it), 2);
1003 
1004 
1005 
1006  // smooth
1007  for (int iters=0; iters<10; ++iters)
1008  {
1009  for (v_it=mesh_.vertices_begin(); v_it!=v_end; ++v_it)
1010  {
1011  if (mesh_.status(*v_it).tagged())
1012  {
1013  u.vectorize(0.0);
1014  ww = 0;
1015 
1016  for (vv_it=mesh_.cvv_iter(*v_it); vv_it.is_valid(); ++vv_it)
1017  {
1018  w = mesh_.property(area_, *vv_it);
1019  u += mesh_.point(*vv_it) * w;
1020  ww += w;
1021  }
1022 
1023  if (ww)
1024  {
1025  u *= (1.0/ww);
1026  u -= mesh_.point(*v_it);
1027  n = mesh_.normal(*v_it);
1028  n *= (u | n);
1029  u -= n;
1030 
1031  u *= 0.1;
1032  }
1033 
1034  mesh_.property(update_, *v_it) = u;
1035  }
1036  }
1037 
1038  for (v_it=mesh_.vertices_begin(); v_it!=v_end; ++v_it)
1039  if (!mesh_.status(*v_it).locked() &&
1040  !mesh_.status(*v_it).feature() &&
1041  !mesh_.is_boundary(*v_it))
1042  mesh_.point(*v_it) += mesh_.property(update_, *v_it);
1043  }
1044 
1045 
1046  // project
1047  if (_use_projection)
1048  for (v_it=mesh_.vertices_begin(); v_it!=v_end; ++v_it)
1049  if (!mesh_.status(*v_it).locked() && !mesh_.status(*v_it).feature())
1050  project_to_reference(*v_it);
1051 
1052 
1053 // if (emit_progress_)
1054 // if (!Progress().step(3))
1055 // break;
1056  }
1057 }
1058 
1059 
1060 //-----------------------------------------------------------------------------
1061 
1062 
1063 template <class Mesh>
1064 void
1065 BaseRemesherT<Mesh>::
1066 remove_caps()
1067 {
1068  typename Mesh::EdgeIter e_it, e_end(mesh_.edges_end());
1069  typename Mesh::HalfedgeHandle h;
1070  typename Mesh::Scalar a0, a1, amin, aa(cos(170.0 * M_PI / 180.0));
1071  typename Mesh::VHandle v, vb, vd;
1072  typename Mesh::FHandle fb, fd;
1073  typename Mesh::Point a, b, c, d;
1074 
1075 
1076  // handle Nastran PIDs for edge flips
1078  mesh_.get_property_handle(pids, "Nastran PIDs");
1079 
1080 
1081  for (e_it=mesh_.edges_begin(); e_it!=e_end; ++e_it)
1082  {
1083  if (!mesh_.status(*e_it).locked() &&
1084  mesh_.is_flip_ok(*e_it))
1085  {
1086  h = mesh_.halfedge_handle(*e_it, 0);
1087  a = mesh_.point(mesh_.to_vertex_handle(h));
1088 
1089  h = mesh_.next_halfedge_handle(h);
1090  b = mesh_.point(vb=mesh_.to_vertex_handle(h));
1091 
1092  h = mesh_.halfedge_handle(*e_it, 1);
1093  c = mesh_.point(mesh_.to_vertex_handle(h));
1094 
1095  h = mesh_.next_halfedge_handle(h);
1096  d = mesh_.point(vd=mesh_.to_vertex_handle(h));
1097 
1098  a0 = ( (a-b).normalize() | (c-b).normalize() );
1099  a1 = ( (a-d).normalize() | (c-d).normalize() );
1100 
1101  if (a0 < a1) { amin = a0; v = vb; }
1102  else { amin = a1; v = vd; }
1103 
1104  // is it a cap?
1105  if (amin < aa)
1106  {
1107  // feature edge and feature vertex -> seems to be intended
1108  if (mesh_.status(*e_it).feature() && mesh_.status(v).feature())
1109  continue;
1110 
1111  // handle PIDs: flip = split + collapse
1112  if (pids.is_valid())
1113  {
1114  fb = mesh_.face_handle(mesh_.halfedge_handle(*e_it, 0));
1115  fd = mesh_.face_handle(mesh_.halfedge_handle(*e_it, 1));
1116 
1117  if (v == vb)
1118  mesh_.property(pids, fb) = mesh_.property(pids, fd);
1119  else
1120  mesh_.property(pids, fd) = mesh_.property(pids, fb);
1121  }
1122 
1123  // project v onto feature edge
1124  if (mesh_.status(*e_it).feature())
1125  mesh_.set_point(v, (a+c)*0.5);
1126 
1127  // flip
1128  mesh_.flip(*e_it);
1129  }
1130  }
1131  }
1132 }
1133 
1134 
1135 //=============================================================================
1136 } // namespace Remeshing
1137 //=============================================================================
void prepare_face_selection()
prepare for remeshing only vertices which are fully surrounded by selected faces (if no face was sele...
void prepare_vertex_selection()
prepare for remeshing only selected vertices (if no vertex was selected, remesh whole mesh) ...
bool is_valid() const
The handle is valid iff the index is not equal to -1.
Definition: Handles.hh:77
Kernel::Point Point
Coordinate type.
Definition: PolyMeshT.hh:115
VertexHandle add_vertex(const Point &_p)
Alias for new_vertex(const Point&).
Definition: PolyMeshT.hh:236
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:93
Kernel::Scalar Scalar
Scalar type.
Definition: PolyMeshT.hh:113
Kernel::ConstVertexFaceIter ConstVertexFaceIter
Circulator.
Definition: PolyMeshT.hh:179
Kernel::VertexHandle VertexHandle
Handle for referencing the corresponding item.
Definition: PolyMeshT.hh:139
VertexHandle split(EdgeHandle _eh, const Point &_p)
Edge split (= 2-to-4 split)
Definition: TriMeshT.hh:266
void update_normals()
Compute normals for all primitives.
Definition: PolyMeshT.cc:241
Functions for selection on a mesh.
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:333
Kernel::Normal Normal
Normal type.
Definition: PolyMeshT.hh:117