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