Developer Documentation
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
HoleFillerT.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 #define HOLEFILLER_CC
52 #include "HoleFillerT.hh"
53 //=============================================================================
54 #include <QInputDialog>
55 
56 template< class MeshT >
58  : mesh_( _mesh )
59 {
60  mesh_.request_vertex_status();
61  mesh_.request_edge_status();
62 
63  if (! mesh_.get_property_handle(scale_,"scale") )
64  mesh_.add_property( scale_ , "scale" );
65 }
66 
67 
68 //=============================================================================
69 
70 
71 
72 template< class MeshT >
74 {
75  mesh_.release_vertex_status();
76  mesh_.release_edge_status();
77 
78  if ( mesh_.get_property_handle(scale_,"scale") )
79  mesh_.remove_property( scale_ );
80 }
81 
82 
83 //=============================================================================
84 //
85 // Identify and fill all holes of the mesh.
86 //
87 //=============================================================================
88 
89 
90 template< class MeshT >
91 void
93 {
94  // Collect all boundary edges
95 
96  std::vector< EH > bdry_edge;
97 
98  for ( EI ei = mesh_.edges_begin();
99  ei != mesh_.edges_end(); ++ei )
100  if ( mesh_.is_boundary( *ei ) )
101  bdry_edge.push_back( *ei );
102 
103  // Fill holes
104 
105  int cnt = 0;
106  for ( typename std::vector< EH >::iterator i = bdry_edge.begin();
107  i != bdry_edge.end(); ++i )
108  if ( mesh_.is_boundary( *i ) )
109  {
110  ++cnt;
111  std::cerr << "Filling hole " << cnt << "\n";
112  fill_hole( *i, _stages );
113  }
114 
115 
116  // Smooth fillings
117 
118  if ( _stages <= 2 )
119  return;
120 
121  std::cerr << "Stage 3 : Smoothing the hole fillings ... ";
122 
124  smoother.initialize( OpenMesh::Smoother::SmootherT< Mesh >::
125  Tangential_and_Normal,
127 
128  smoother.smooth( 500 );
129  std::cerr << "ok\n";
130 }
131 
132 
133 //=============================================================================
134 //
135 // Fill a hole which is identified by one of its boundary edges.
136 //
137 //=============================================================================
138 
139 
140 template< class MeshT >
141 void
142 HoleFiller< MeshT >::fill_hole( EH _eh, int _stages )
143 {
144  std::cerr << " Stage 1 : Computing a minimal triangulation ... ";
145 
146  //remember last vertex for selection of new ones
147  typename MeshT::VertexHandle old_last_handle = *(--mesh_.vertices_end());
148 
149  // No boundary edge, no hole
150 
151  if ( ! mesh_.is_boundary( _eh ) )
152  return;
153 
154 
155  // Get boundary halfedge
156 
157  HH hh = mesh_.halfedge_handle( _eh, 0 );
158  if ( ! mesh_.is_boundary( hh ) )
159  hh = mesh_.opposite_halfedge_handle( hh );
160 
161 
162  // Collect boundary vertices
163 
164  boundary_vertex_.clear();
165  opposite_vertex_.clear();
166 
167  HH ch = hh;
168 
169  do {
170  boundary_vertex_.push_back( mesh_.from_vertex_handle( ch ) );
171  opposite_vertex_.push_back( mesh_.to_vertex_handle
172  (mesh_.next_halfedge_handle( mesh_.opposite_halfedge_handle( ch ) ) ) );
173 
174  //check number of outgoing boundary HEH's at Vertex
175  int c = 0;
176  VH vh = mesh_.to_vertex_handle(ch);
177 
178  for ( typename MeshT::VertexOHalfedgeIter voh_it(mesh_,vh); voh_it.is_valid(); ++voh_it)
179  if ( mesh_.is_boundary( *voh_it ) )
180  c++;
181 
182  if ( c >= 2){
183  HH op = mesh_.opposite_halfedge_handle( ch );
184  typename MeshT::VertexOHalfedgeIter voh_it(mesh_,op);
185 
186  ch = *(++voh_it);
187  }else
188  ch = mesh_.next_halfedge_handle( ch );
189 
190  } while ( ch != hh );
191 
192 
193  int nv = boundary_vertex_.size();
194 
195  // Compute an initial triangulation
196 
197  w_.clear();
198  w_.resize( nv, std::vector<Weight>( nv, Weight() ) );
199 
200  l_.clear();
201  l_.resize( nv, std::vector<int>( nv, 0 ) );
202 
203 
204  for ( int i = 0; i < nv - 1; ++i )
205  w_[i][i+1] = Weight( 0, 0 );
206 
207  for ( int j = 2; j < nv; ++j )
208  {
209  #pragma omp parallel for shared(j, nv)
210  for(int i = 0; i < nv-j; ++i)
211  {
212  Weight valmin;
213  int argmin = -1;
214  for ( int m = i + 1; m < i + j; ++m )
215  {
216  Weight newval = w_[i][m] + w_[m][i+j] + weight( i, m, i+j );
217  if ( newval < valmin )
218  {
219  valmin = newval;
220  argmin = m;
221  }
222  }
223 
224  w_[i][i+j] = valmin;
225  l_[i][i+j] = argmin;
226  }
227  }
228 
229 
230  // Actually fill the hole. We collect all triangles and edges of
231  // this filling for further processing.
232 
233  hole_edge_.clear();
234  hole_triangle_.clear();
235  if ( fill( 0, nv - 1 ) ){
236 
237  std::cerr << "ok\n";
238 
239  if ( _stages <= 1 )
240  return;
241 
242  std::cerr << " Stage 2 : Fairing the filling ... ";
243 
244  std::vector< FH > handles = hole_triangle_;
245 
246  fairing(handles);
247 
248  //select all new vertices
249  typename MeshT::VertexIter old_end = ++typename MeshT::VertexIter(mesh_,old_last_handle);
250  typename MeshT::VertexIter v_end = mesh_.vertices_end();
251 
252  for(; old_end != v_end; ++old_end)
253  if ( !mesh_.status(*old_end).deleted() )
254  mesh_.status(*old_end).set_selected( true );
255 
256  std::cerr << "ok\n";
257  }else
258  std::cerr << "Could not create triangulation" << std::endl;
259 }
260 
261 
263 template< class MeshT >
264 void
265 HoleFiller< MeshT >::fairing( std::vector< FH >& _faceHandles ){
266 
267  //generate vector of all edges
268  hole_edge_.clear();
269 
270  hole_triangle_ = _faceHandles;
271 
276 
277  if (! mesh_.get_property_handle(edgeProp,"edgeProp") )
278  mesh_.add_property( edgeProp, "edgeProp" );
279  if (! mesh_.get_property_handle(vertexProp,"vertexProp") )
280  mesh_.add_property( vertexProp, "vertexProp" );
281  if (! mesh_.get_property_handle(faceProp,"faceProp") )
282  mesh_.add_property( faceProp, "faceProp" );
283  if (! mesh_.get_property_handle(orderProp,"orderProp") )
284  mesh_.add_property( orderProp, "orderProp" );
285 
286  EI eIt;
287  EI eEnd = mesh_.edges_end();
288  VI vIt;
289  VI vEnd = mesh_.vertices_end();
290  FI fIt;
291  FI fEnd = mesh_.faces_end();
292 
293  //init properties by setting all of them to false
294  for (fIt = mesh_.faces_begin(); fIt != fEnd; ++fIt){
295  mesh_.property( orderProp, *fIt ) = false;
296  mesh_.property( faceProp, *fIt ) = false;
297  }
298 
299  for (eIt = mesh_.edges_begin(); eIt != eEnd; ++eIt)
300  mesh_.property( edgeProp, *eIt ) = false;
301 
302  for (vIt = mesh_.vertices_begin(); vIt != vEnd; ++vIt){
303  mesh_.property( vertexProp, *vIt ) = false;
304  }
305 
306  //set face property
307  for (uint i = 0; i < hole_triangle_.size(); i++){
308  mesh_.property( faceProp, hole_triangle_[i] ) = true;
309  }
310 
311  //set properties
312  for (unsigned int i = 0; i < hole_triangle_.size(); i++){
313  for (FEI fei = mesh_.fe_iter( hole_triangle_[i] ); fei.is_valid(); ++fei){
314  mesh_.status( *fei ).set_locked(true);
315  //set edge property for all edges inside the hole (eg not on the hole boundary)
316  if (mesh_.property( faceProp, mesh_.face_handle(mesh_.halfedge_handle(*fei, 0) ) ) &&
317  mesh_.property( faceProp, mesh_.face_handle(mesh_.halfedge_handle(*fei, 1) ) ) ){
318 
319  mesh_.property( edgeProp, *fei ) = true;
320  hole_edge_.push_back( *fei );
321  mesh_.status( *fei ).set_locked(false);
322  }
323  }
324 
326  for (FVI fvi = mesh_.fv_iter( hole_triangle_[i] ); fvi.is_valid(); ++fvi){
327  //set vertex property for all vertices of the hole
328  for ( VFI vfi = mesh_.vf_iter( *fvi ); vfi.is_valid(); ++vfi )
329  mesh_.property( vertexProp, *fvi ) = true;
330  }
331  }
332 
333  //calculate scaling weights for vertices
334  for (vIt = mesh_.vertices_begin(); vIt != vEnd; ++vIt)
335  if (mesh_.property(vertexProp, *vIt)){
336 
337  Scalar cnt = 0;
338  Scalar scale = 0;
339 
340  for ( VOHI voh_it = mesh_.voh_iter( *vIt ); voh_it.is_valid(); ++voh_it )
341  {
342  HH h2 = mesh_.opposite_halfedge_handle( *voh_it );
343 
344  if (mesh_.face_handle(*voh_it).is_valid() &&
345  mesh_.face_handle(h2).is_valid() &&
346  mesh_.property(faceProp, mesh_.face_handle( *voh_it ) ) &&
347  mesh_.property(faceProp, mesh_.face_handle(h2) ))
348  continue;
349 
350  cnt += 1.0f;
351  Point p0 = mesh_.point( *vIt );
352  Point p1 = mesh_.point( mesh_.to_vertex_handle( *voh_it ) );
353  scale += ( p1 - p0 ).norm();
354 
355  }
356 
357  scale /= cnt;
358 
359  mesh_.property( scale_, *vIt ) = scale;
360  }
361 
362  mesh_.remove_property(edgeProp);
363  mesh_.remove_property(vertexProp);
364  mesh_.remove_property(faceProp);
365  mesh_.remove_property(orderProp);
366 
367  // Do the patch fairing
368 
369  bool did_refine = true;
370 
371  for ( int k = 0; k < 40 && did_refine; ++k )
372  {
373  uint end = hole_triangle_.size();
374 
375  did_refine = false;
376  for ( unsigned int j = 0; j < end; ++j )
377  did_refine |= refine( hole_triangle_[j] );
378 
379  for ( int i = 0; i < 10; ++i )
380  for ( unsigned int j = 0; j < hole_edge_.size(); ++j )
381  relax_edge( hole_edge_[j] );
382  }
383 
384  // unlock everything
385  for ( EI ei = mesh_.edges_begin(); ei != mesh_.edges_end(); ++ei )
386  mesh_.status( *ei ).set_locked( false );
387 }
388 
389 //=============================================================================
390 //
391 // Refine a face: Apply a 1-to-3 split if the edge lengths of the
392 // face do not match the interpolated edge lengths of the hole
393 // boundary.
394 //
395 //=============================================================================
396 
397 
398 template< class MeshT >
399 bool
401 {
402 
403  // Collect the three edges of the face into e0, e1, e2
404 
405  FEI fei = mesh_.fe_iter( _fh );
406  EH e0 = *fei; ++fei;
407  EH e1 = *fei; ++fei;
408  EH e2 = *fei; ++fei;
409 
410 
411  // Collect the vertices, vertex positions and scale factors of the face.
412 
413  FVI fvi = mesh_.fv_iter( _fh );
414 
415  VH v0 = *fvi; ++fvi;
416  VH v1 = *fvi; ++fvi;
417  VH v2 = *fvi; ++fvi;
418 
419  Point p0 = mesh_.point( v0 );
420  Point p1 = mesh_.point( v1 );
421  Point p2 = mesh_.point( v2 );
422 
423  Scalar scale0 = mesh_.property( scale_, v0 );
424  Scalar scale1 = mesh_.property( scale_, v1 );
425  Scalar scale2 = mesh_.property( scale_, v2 );
426 
427  // Interpolate the scale factor.
428 
429  Scalar scale = ( scale0 + scale1 + scale2 ) / 3.0f;
430  Point center = typename MeshT::Scalar(1.0/3.0) * ( p0 + p1 + p2 );
431 
432  Scalar d0 = 1.0f * ( p0 - center ).norm();
433  Scalar d1 = 1.0f * ( p1 - center ).norm();
434  Scalar d2 = 1.0f * ( p2 - center ).norm();
435 
436 
437  //dont split triangles which tend to degenerate
438  if ( (d0 + d1 + d2) / 3.0f < scale) return false;
439 
440 
441  // If the edge lengths differ too much from the scale, perform a
442  // triangle split.
443 
444  if ( ( d0 > scale && d0 > scale0 ) ||
445  ( d1 > scale && d1 > scale1 ) ||
446  ( d2 > scale && d2 > scale2 ) )
447  {
448  // Split the face ...
449 
450  VH ch = mesh_.add_vertex( center );
451  mesh_.split( _fh, ch );
452 
453  // ... put the new triangles into the global triangle list ...
454 
455  for ( VFI vfi = mesh_.vf_iter( ch ); vfi.is_valid(); ++vfi )
456  if ( *vfi != _fh )
457  hole_triangle_.push_back( *vfi );
458 
459  // ... put the new edges into the global edge list ...
460 
461  for ( VEI vei = mesh_.ve_iter( ch ); vei.is_valid(); ++vei )
462  hole_edge_.push_back( *vei );
463 
464  // ... and set the appropriate scale factor for the new vertex.
465 
466  mesh_.property( scale_, ch ) = scale;
467 
468  relax_edge( e0 );
469  relax_edge( e1 );
470  relax_edge( e2 );
471 
472  return true;
473  }
474 
475  return false;
476 }
477 
478 
479 //=============================================================================
480 //
481 // Relax an edge: Flip it if one of its opposing vertices lies in
482 // the circumsphere of the other three vertices.
483 //
484 //=============================================================================
485 
486 
487 template< class MeshT >
488 bool
490 {
491  if ( mesh_.status( _eh ).locked() )
492  return false;
493 
494  // Abbreviations for the two halfedges of _eh
495 
496  HH h0 = mesh_.halfedge_handle( _eh, 0 );
497  HH h1 = mesh_.halfedge_handle( _eh, 1 );
498 
499  // Get the two end-vertices u and v of the edge
500 
501  Point u( mesh_.point( mesh_.to_vertex_handle( h0 ) ) );
502  Point v( mesh_.point( mesh_.to_vertex_handle( h1 ) ) );
503 
504  // Get the two opposing vertices a and b
505 
506  Point a( mesh_.point
507  ( mesh_.to_vertex_handle
508  ( mesh_.next_halfedge_handle( h0 ) ) ) );
509 
510  Point b( mesh_.point
511  ( mesh_.to_vertex_handle
512  ( mesh_.next_halfedge_handle( h1 ) ) ) );
513 
514  // If the circumsphere criterion is fullfilled AND if the flip is
515  // topologically admissible, we do it.
516 
517  if ( in_circumsphere( a, u, v, b ) || in_circumsphere( b, u, v, a ) ){
518  if ( mesh_.is_flip_ok( _eh ) )
519  {
520 
521  mesh_.flip( _eh );
522 
523  return true;
524  }else
525  mesh_.status(_eh).set_selected( true );
526 }
527  return false;
528 }
529 
530 
531 //=============================================================================
532 //
533 // Test whether a point _x lies in the circumsphere of _a,_b,_c.
534 //
535 //=============================================================================
536 
537 
538 template< class MeshT >
539 bool
540 HoleFiller< MeshT >::in_circumsphere( const Point & _x,
541  const Point & _a,
542  const Point & _b,
543  const Point & _c ) const
544 {
545  Point ab = _b - _a;
546  Point ac = _c - _a;
547 
548  Scalar a00 = -2.0f * ( ab | _a );
549  Scalar a01 = -2.0f * ( ab | _b );
550  Scalar a02 = -2.0f * ( ab | _c );
551  Scalar b0 = _a.sqrnorm() - _b.sqrnorm();
552 
553  Scalar a10 = -2.0f * ( ac | _a );
554  Scalar a11 = -2.0f * ( ac | _b );
555  Scalar a12 = -2.0f * ( ac | _c );
556  Scalar b1 = _a.sqrnorm() - _c.sqrnorm();
557 
558  typename MeshT::Scalar alpha = -(-a11*a02+a01*a12-a12*b0+b1*a02+a11*b0-a01*b1)
559  / (-a11*a00+a11*a02-a10*a02+a00*a12+a01*a10-a01*a12);
560  typename MeshT::Scalar beta = (a10*b0-a10*a02-a12*b0+a00*a12+b1*a02-a00*b1)
561  / (-a11*a00+a11*a02-a10*a02+a00*a12+a01*a10-a01*a12);
562  typename MeshT::Scalar gamma = (-a11*a00-a10*b0+a00*b1+a11*b0+a01*a10-a01*b1)
563  / (-a11*a00+a11*a02-a10*a02+a00*a12+a01*a10-a01*a12);
564 
565  Point center = alpha * _a + beta * _b + gamma * _c;
566 
567  return ( _x - center ).sqrnorm() < ( _a - center ).sqrnorm();
568 }
569 
570 
571 //=============================================================================
572 //
573 // Create the triangulation
574 //
575 // Recursively creates the triangulation for polygon (_i,...,_j).
576 //
577 //=============================================================================
578 
579 
580 template< class MeshT >
581 bool
582 HoleFiller< MeshT >::fill( int _i, int _j )
583 {
584  // If the two vertices _i and _j are adjacent, there is nothing to do.
585 
586  if ( _i + 1 == _j )
587  return true;
588 
589 
590  // Create and store the middle triangle, store its edges.
591 
592  FH fh = mesh_.add_face( boundary_vertex_[_i],
593  boundary_vertex_[ l_[_i][_j] ],
594  boundary_vertex_[_j] );
595  hole_triangle_.push_back( fh );
596 
597  if (!fh.is_valid())
598  return false;
599 
600  hole_edge_.push_back( mesh_.edge_handle
601  ( mesh_.find_halfedge( boundary_vertex_[_i],
602  boundary_vertex_[ l_[_i][_j] ] ) ) );
603  hole_edge_.push_back( mesh_.edge_handle
604  ( mesh_.find_halfedge( boundary_vertex_[ l_[_i][_j] ],
605  boundary_vertex_[_j] ) ) );
606 
607 
608  // Recursively create the left and right side of the
609  // triangulation.
610 
611  if (!fill( _i, l_[_i][_j] ) || !fill( l_[_i][_j], _j ))
612  return false;
613  else
614  return true;
615 }
616 
617 
618 
619 //=============================================================================
620 //
621 // Compute the weight of the triangle (_i,_j,_k).
622 //
623 //=============================================================================
624 
625 
626 template< class MeshT >
628 HoleFiller< MeshT >::weight( int _i, int _j, int _k )
629 {
630  // Return an infinite weight if the insertion of this triangle
631  // would create complex edges.
632 
633  if ( exists_edge( boundary_vertex_[_i], boundary_vertex_[_j] ) ||
634  exists_edge( boundary_vertex_[_j], boundary_vertex_[_k] ) ||
635  exists_edge( boundary_vertex_[_k], boundary_vertex_[_i] ) )
636  return Weight();
637 
638 
639  // Return an infinite weight, if one of the neighboring patches
640  // could not be created.
641 
642 
643  if ( l_[_i][_j] == -1 ) return Weight();
644  if ( l_[_j][_k] == -1 ) return Weight();
645 
646 
647  // Compute the maxmimum dihedral angles to the adjacent triangles
648  // (if they exist)
649 
650  Scalar angle = 0.0f;
651 
652  if ( _i + 1 == _j )
653  angle = std::max( angle, dihedral_angle( boundary_vertex_[_i],
654  boundary_vertex_[_j],
655  boundary_vertex_[_k],
656  opposite_vertex_[_i] ) );
657  else
658  angle = std::max( angle, dihedral_angle( boundary_vertex_[_i],
659  boundary_vertex_[_j],
660  boundary_vertex_[_k],
661  boundary_vertex_[l_[_i][_j]] ) );
662 
663  if ( _j + 1 == _k )
664  angle = std::max( angle, dihedral_angle( boundary_vertex_[_j],
665  boundary_vertex_[_k],
666  boundary_vertex_[_i],
667  opposite_vertex_[_j] ) );
668  else
669  angle = std::max( angle, dihedral_angle( boundary_vertex_[_j],
670  boundary_vertex_[_k],
671  boundary_vertex_[_i],
672  boundary_vertex_[l_[_j][_k]] ) );
673 
674  if ( _i == 0 && _k == (int) boundary_vertex_.size() - 1 )
675  angle = std::max( angle, dihedral_angle( boundary_vertex_[_k],
676  boundary_vertex_[_i],
677  boundary_vertex_[_j],
678  opposite_vertex_[_k] ) );
679 
680 
681  return Weight( angle, area( boundary_vertex_[_i],
682  boundary_vertex_[_j],
683  boundary_vertex_[_k] ) );
684 }
685 
686 
687 
688 //=============================================================================
689 //
690 // Does an edge from vertex _u to _v exist?
691 //
692 //=============================================================================
693 
694 
695 template< class MeshT >
696 bool
697 HoleFiller< MeshT >::exists_edge( VH _u, VH _w )
698 {
699  for ( VOHI vohi = mesh_.voh_iter( _u ); vohi.is_valid(); ++vohi )
700  if ( ! mesh_.is_boundary( mesh_.edge_handle( *vohi ) ) )
701  if ( mesh_.to_vertex_handle( *vohi ) == _w )
702  return true;
703 
704  return false;
705 }
706 
707 
708 
709 //=============================================================================
710 //
711 // Compute the area of the triangle (_a,_b,_c).
712 //
713 //=============================================================================
714 
715 
716 template< class MeshT >
717 typename MeshT::Scalar
718 HoleFiller< MeshT >::area( VH _a, VH _b, VH _c )
719 {
720  Point a( mesh_.point( _a ) );
721  Point b( mesh_.point( _b ) );
722  Point c( mesh_.point( _c ) );
723 
724  Point n( ( b - a ) % ( c - b ) );
725 
726  return 0.5 * n.norm();
727 }
728 
729 
730 //=============================================================================
731 //
732 // Compute a dihedral angle
733 //
734 // Computes the dihedral angle (in degrees) between triangle
735 // (_u,_v,_a) and triangle (_v,_u,_b), no matter whether these
736 // triangles actually exist in the current mesh or not).
737 //
738 //=============================================================================
739 
740 
741 template< class MeshT >
742 typename MeshT::Scalar
743 HoleFiller< MeshT >::dihedral_angle( VH _u, VH _v, VH _a, VH _b )
744 {
745  Point u( mesh_.point( _u ) );
746  Point v( mesh_.point( _v ) );
747  Point a( mesh_.point( _a ) );
748  Point b( mesh_.point( _b ) );
749 
750  Point n0( ( v - u ) % ( a - v ) );
751  Point n1( ( u - v ) % ( b - u ) );
752 
753  n0.normalize();
754  n1.normalize();
755 
756  return acos( n0 | n1 ) * 180.0 / M_PI;
757 }
758 
759 
761 template< class MeshT >
762 void
763 HoleFiller< MeshT >::removeDegeneratedFaces( std::vector< FH >& _faceHandles ){
764 
765  for (int i = _faceHandles.size()-1; i >= 0 ; i--){
766 
767  if ( mesh_.status( _faceHandles[i] ).deleted() ){
768  // face might be deleted because of a previous edge collapse
769  // erase the face from the vector
770  _faceHandles.erase( _faceHandles.begin() + i );
771  continue;
772  }
773 
774  //get the vertices (works only on triMeshes)
775  FVI fvi = mesh_.fv_iter( _faceHandles[i] );
776  Point v0 = mesh_.point( *fvi);
777  ++fvi;
778  Point v1 = mesh_.point( *fvi );
779  ++fvi;
780  Point v2 = mesh_.point( *fvi );
781 
782  //check if its degenerated
783  Point v0v1 = v1 - v0;
784  Point v0v2 = v2 - v0;
785  Point n = v0v1 % v0v2; // not normalized !
786  double d = n.sqrnorm();
787 
788  if (d < FLT_MIN && d > -FLT_MIN) {
789  // degenerated face found
790  FHI hIt = mesh_.fh_iter( _faceHandles[i] );
791 
792  //try to collapse one of the edges
793  while (hIt.is_valid()){
794  if ( mesh_.is_collapse_ok( *hIt ) ){
795  // collapse the edge to remove the triangle
796  mesh_.collapse( *hIt );
797  // and erase the corresponding face from the vector
798  _faceHandles.erase( _faceHandles.begin() + i );
799  break;
800  } else {
801  ++hIt;
802  }
803  }
804  }
805  }
806 }
807 
808 //=============================================================================
809 
T angle(T _cos_angle, T _sin_angle)
Definition: MathDefs.hh:145