KiCad PCB EDA Suite
ttl.h
Go to the documentation of this file.
1 /*
2  * Copyright (C) 1998, 2000-2007, 2010, 2011, 2012, 2013 SINTEF ICT,
3  * Applied Mathematics, Norway.
4  *
5  * Contact information: E-mail: tor.dokken@sintef.no
6  * SINTEF ICT, Department of Applied Mathematics,
7  * P.O. Box 124 Blindern,
8  * 0314 Oslo, Norway.
9  *
10  * This file is part of TTL.
11  *
12  * TTL is free software: you can redistribute it and/or modify
13  * it under the terms of the GNU Affero General Public License as
14  * published by the Free Software Foundation, either version 3 of the
15  * License, or (at your option) any later version.
16  *
17  * TTL is distributed in the hope that it will be useful,
18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20  * GNU Affero General Public License for more details.
21  *
22  * You should have received a copy of the GNU Affero General Public
23  * License along with TTL. If not, see
24  * <http://www.gnu.org/licenses/>.
25  *
26  * In accordance with Section 7(b) of the GNU Affero General Public
27  * License, a covered work must retain the producer line in every data
28  * file that is created or manipulated using TTL.
29  *
30  * Other Usage
31  * You can be released from the requirements of the license by purchasing
32  * a commercial license. Buying such a license is mandatory as soon as you
33  * develop commercial activities involving the TTL library without
34  * disclosing the source code of your own applications.
35  *
36  * This file may be used in accordance with the terms contained in a
37  * written agreement between you and SINTEF ICT.
38  */
39 
40 #ifndef _TTL_H_
41 #define _TTL_H_
42 
43 #include <list>
44 #include <iterator>
45 
46 // Debugging
47 #ifdef DEBUG_TTL
48 static void errorAndExit( char* aMessage )
49 {
50  cout << "\n!!! ERROR: " << aMessage << " !!!\n" << endl;
51  exit(-1);
52 }
53 #endif
54 
55 // Next on TOPOLOGY:
56 // - get triangle strips
57 // - weighted graph, algorithms using a weight (real) for each edge,
58 // e.g. an "abstract length". Use for minimum spanning tree
59 // or some arithmetics on weights?
60 // - Circulators as defined in CGAL with more STL compliant code
61 
62 // - analyze in detail locateFace: e.g. detect 0-orbit in case of infinite loop
63 // around a node etc.
64 
113 namespace ttl
114 {
116 {
117 #ifndef DOXYGEN_SHOULD_SKIP_THIS
118 
119 public:
121  m_triangulation( aTriang )
122  {
123  }
124 
125  // Delaunay Triangulation
126  template <class TRAITS_TYPE, class DART_TYPE, class POINT_TYPE>
127  bool InsertNode( DART_TYPE& aDart, POINT_TYPE& aPoint );
128 
129  template <class TRAITS_TYPE, class DART_TYPE>
130  void RemoveRectangularBoundary( DART_TYPE& aDart );
131 
132  template <class TRAITS_TYPE, class DART_TYPE>
133  void RemoveNode( DART_TYPE& aDart );
134 
135  template <class TRAITS_TYPE, class DART_TYPE>
136  void RemoveBoundaryNode( DART_TYPE& aDart );
137 
138  template <class TRAITS_TYPE, class DART_TYPE>
139  void RemoveInteriorNode( DART_TYPE& aDart );
140 
141  // Topological and Geometric Queries
142  // ---------------------------------
143  template <class TRAITS_TYPE, class POINT_TYPE, class DART_TYPE>
144  static bool LocateFaceSimplest( const POINT_TYPE& aPoint, DART_TYPE& aDart );
145 
146  template <class TRAITS_TYPE, class POINT_TYPE, class DART_TYPE>
147  static bool LocateTriangle( const POINT_TYPE& aPoint, DART_TYPE& aDart );
148 
149  template <class TRAITS_TYPE, class POINT_TYPE, class DART_TYPE>
150  static bool InTriangle( const POINT_TYPE& aPoint, const DART_TYPE& aDart );
151 
152  template <class DART_TYPE, class DART_LIST_TYPE>
153  static void GetBoundary( const DART_TYPE& aDart, DART_LIST_TYPE& aBoundary );
154 
155  template <class DART_TYPE>
156  static bool IsBoundaryEdge( const DART_TYPE& aDart );
157 
158  template <class DART_TYPE>
159  static bool IsBoundaryFace( const DART_TYPE& aDart );
160 
161  template <class DART_TYPE>
162  static bool IsBoundaryNode( const DART_TYPE& aDart );
163 
164  template <class DART_TYPE>
165  static int GetDegreeOfNode( const DART_TYPE& aDart );
166 
167  template <class DART_TYPE, class DART_LIST_TYPE>
168  static void Get0OrbitInterior( const DART_TYPE& aDart, DART_LIST_TYPE& aOrbit );
169 
170  template <class DART_TYPE, class DART_LIST_TYPE>
171  static void Get0OrbitBoundary( const DART_TYPE& aDart, DART_LIST_TYPE& aOrbit );
172 
173  template <class DART_TYPE>
174  static bool Same0Orbit( const DART_TYPE& aD1, const DART_TYPE& aD2 );
175 
176  template <class DART_TYPE>
177  static bool Same1Orbit( const DART_TYPE& aD1, const DART_TYPE& aD2 );
178 
179  template <class DART_TYPE>
180  static bool Same2Orbit( const DART_TYPE& aD1, const DART_TYPE& aD2 );
181 
182  template <class TRAITS_TYPE, class DART_TYPE>
183  static bool SwappableEdge( const DART_TYPE& aDart, bool aAllowDegeneracy = false );
184 
185  template <class DART_TYPE>
186  static void PositionAtNextBoundaryEdge( DART_TYPE& aDart );
187 
188  template <class TRAITS_TYPE, class DART_TYPE>
189  static bool ConvexBoundary( const DART_TYPE& aDart );
190 
191  // Utilities for Delaunay Triangulation
192  // ------------------------------------
193  template <class TRAITS_TYPE, class DART_TYPE, class DART_LIST_TYPE>
194  void OptimizeDelaunay( DART_LIST_TYPE& aElist );
195 
196  template <class TRAITS_TYPE, class DART_TYPE, class DART_LIST_TYPE>
197  void OptimizeDelaunay( DART_LIST_TYPE& aElist, const typename DART_LIST_TYPE::iterator aEnd );
198 
199  template <class TRAITS_TYPE, class DART_TYPE>
200  bool SwapTestDelaunay( const DART_TYPE& aDart, bool aCyclingCheck = false ) const;
201 
202  template <class TRAITS_TYPE, class DART_TYPE>
203  void RecSwapDelaunay( DART_TYPE& aDiagonal );
204 
205  template <class TRAITS_TYPE, class DART_TYPE, class LIST_TYPE>
206  void SwapEdgesAwayFromInteriorNode( DART_TYPE& aDart, LIST_TYPE& aSwappedEdges );
207 
208  template <class TRAITS_TYPE, class DART_TYPE, class LIST_TYPE>
209  void SwapEdgesAwayFromBoundaryNode( DART_TYPE& aDart, LIST_TYPE& aSwappedEdges );
210 
211  template <class TRAITS_TYPE, class DART_TYPE, class DART_LIST_TYPE>
212  void SwapEdgeInList( const typename DART_LIST_TYPE::iterator& aIt, DART_LIST_TYPE& aElist );
213 
214  // Constrained Triangulation
215  // -------------------------
216  template <class TRAITS_TYPE, class DART_TYPE>
217  static DART_TYPE InsertConstraint( DART_TYPE& aDStart, DART_TYPE& aDEnd, bool aOptimizeDelaunay );
218 
219 private:
221 
222  template <class TRAITS_TYPE, class FORWARD_ITERATOR, class DART_TYPE>
223  void insertNodes( FORWARD_ITERATOR aFirst, FORWARD_ITERATOR aLast, DART_TYPE& aDart );
224 
225  template <class TOPOLOGY_ELEMENT_TYPE, class DART_TYPE>
226  static bool isMemberOfFace( const TOPOLOGY_ELEMENT_TYPE& aTopologyElement, const DART_TYPE& aDart );
227 
228  template <class TRAITS_TYPE, class NODE_TYPE, class DART_TYPE>
229  static bool locateFaceWithNode( const NODE_TYPE& aNode, DART_TYPE& aDartIter );
230 
231  template <class DART_TYPE>
232  static void getAdjacentTriangles( const DART_TYPE& aDart, DART_TYPE& aT1, DART_TYPE& aT2,
233  DART_TYPE& aT3 );
234 
235  template <class DART_TYPE>
236  static void getNeighborNodes( const DART_TYPE& aDart, std::list<DART_TYPE>& aNodeList,
237  bool& aBoundary );
238 
239  template <class TRAITS_TYPE, class DART_TYPE>
240  static bool degenerateTriangle( const DART_TYPE& aDart );
241 };
242 
243 #endif // DOXYGEN_SHOULD_SKIP_THIS
244 
245 
288 template <class TRAITS_TYPE, class DART_TYPE, class POINT_TYPE>
289 bool TRIANGULATION_HELPER::InsertNode( DART_TYPE& aDart, POINT_TYPE& aPoint )
290 {
291  bool found = LocateTriangle<TRAITS_TYPE>( aPoint, aDart );
292 
293  if( !found )
294  {
295 #ifdef DEBUG_TTL
296  cout << "ERROR: Triangulation::insertNode: NO triangle found. /n";
297 #endif
298  return false;
299  }
300 
301  // ??? can we hide the dart? this is not possible if one triangle only
302  m_triangulation.splitTriangle( aDart, aPoint );
303 
304  DART_TYPE d1 = aDart;
305  d1.Alpha2().Alpha1().Alpha2().Alpha0().Alpha1();
306 
307  DART_TYPE d2 = aDart;
308  d2.Alpha1().Alpha0().Alpha1();
309 
310  // Preserve a dart as output incident to the node and CCW
311  DART_TYPE d3 = aDart;
312  d3.Alpha2();
313  aDart = d3; // and see below
314  //DART_TYPE dsav = d3;
315  d3.Alpha0().Alpha1();
316 
317  //if (!TRAITS_TYPE::fixedEdge(d1) && !IsBoundaryEdge(d1)) {
318  if( !IsBoundaryEdge( d1 ) )
319  {
320  d1.Alpha2();
321  RecSwapDelaunay<TRAITS_TYPE>( d1 );
322  }
323 
324  //if (!TRAITS_TYPE::fixedEdge(d2) && !IsBoundaryEdge(d2)) {
325  if( !IsBoundaryEdge( d2 ) )
326  {
327  d2.Alpha2();
328  RecSwapDelaunay<TRAITS_TYPE>( d2 );
329  }
330 
331  // Preserve the incoming dart as output incident to the node and CCW
332  //d = dsav.Alpha2();
333  aDart.Alpha2();
334  //if (!TRAITS_TYPE::fixedEdge(d3) && !IsBoundaryEdge(d3)) {
335  if( !IsBoundaryEdge( d3 ) )
336  {
337  d3.Alpha2();
338  RecSwapDelaunay<TRAITS_TYPE>( d3 );
339  }
340 
341  return true;
342 }
343 
344 //------------------------------------------------------------------------------------------------
345 // Private/Hidden function (might change later)
346 template <class TRAITS_TYPE, class FORWARD_ITERATOR, class DART_TYPE>
347 void TRIANGULATION_HELPER::insertNodes( FORWARD_ITERATOR aFirst, FORWARD_ITERATOR aLast,
348  DART_TYPE& aDart )
349 {
350 
351  // Assumes that the dereferenced point objects are pointers.
352  // References to the point objects are then passed to TTL.
353 
354  FORWARD_ITERATOR it;
355  for( it = aFirst; it != aLast; ++it )
356  {
357  InsertNode<TRAITS_TYPE>( aDart, **it );
358  }
359 }
360 
361 
378 template <class TRAITS_TYPE, class DART_TYPE>
380 {
381  DART_TYPE d_next = aDart;
382  DART_TYPE d_iter;
383 
384  for( int i = 0; i < 4; i++ )
385  {
386  d_iter = d_next;
387  d_next.Alpha0();
388  PositionAtNextBoundaryEdge( d_next );
389  RemoveBoundaryNode<TRAITS_TYPE>( d_iter );
390  }
391 
392  aDart = d_next; // Return a dart at the new boundary
393 }
394 
406 template <class TRAITS_TYPE, class DART_TYPE>
407 void TRIANGULATION_HELPER::RemoveNode( DART_TYPE& aDart )
408 {
409 
410  if( isBoundaryNode( aDart ) )
411  RemoveBoundaryNode<TRAITS_TYPE>( aDart );
412  else
413  RemoveInteriorNode<TRAITS_TYPE>( aDart );
414 }
415 
426 template <class TRAITS_TYPE, class DART_TYPE>
428 {
429 
430  // ... and update Delaunay
431  // - CCW dart must be given (for remove)
432  // - No dart is delivered back now (but this is possible if
433  // we assume that there is not only one triangle left in the m_triangulation.
434 
435  // Position at boundary edge and CCW
436  if( !IsBoundaryEdge( aDart ) )
437  {
438  aDart.Alpha1(); // ensures that next function delivers back a CCW dart (if the given dart is CCW)
440  }
441 
442  std::list<DART_TYPE> swapped_edges;
443  SwapEdgesAwayFromBoundaryNode<TRAITS_TYPE>( aDart, swapped_edges );
444 
445  // Remove boundary triangles and remove the new boundary from the list
446  // of swapped edges, see below.
447  DART_TYPE d_iter = aDart;
448  DART_TYPE dnext = aDart;
449  bool bend = false;
450  while( bend == false )
451  {
452  dnext.Alpha1().Alpha2();
453  if( IsBoundaryEdge( dnext ) )
454  bend = true; // Stop when boundary
455 
456  // Generic: Also remove the new boundary from the list of swapped edges
457  DART_TYPE n_bedge = d_iter;
458  n_bedge.Alpha1().Alpha0().Alpha1().Alpha2(); // new boundary edge
459 
460  // ??? can we avoid find if we do this in swap away?
461  typename std::list<DART_TYPE>::iterator it;
462  it = find( swapped_edges.begin(), swapped_edges.end(), n_bedge );
463 
464  if( it != swapped_edges.end() )
465  swapped_edges.erase( it );
466 
467  // Remove the boundary triangle
469  d_iter = dnext;
470  }
471 
472  // Optimize Delaunay
473  typedef std::list<DART_TYPE> DART_LIST_TYPE;
474  OptimizeDelaunay<TRAITS_TYPE, DART_TYPE, DART_LIST_TYPE>( swapped_edges );
475 }
476 
477 
492 template <class TRAITS_TYPE, class DART_TYPE>
494 {
495  // ... and update to Delaunay.
496  // Must allow degeneracy temporarily, see comments in swap edges away
497  // Assumes:
498  // - revese_splitTriangle does not affect darts
499  // outside the resulting triangle.
500 
501  // 1) Swaps edges away from the node until degree=3 (generic)
502  // 2) Removes the remaining 3 triangles and creates a new to fill the hole
503  // unsplitTriangle which is required
504  // 3) Runs LOP on the platelet to obtain a Delaunay m_triangulation
505  // (No dart is delivered as output)
506 
507  // Assumes dart is counterclockwise
508 
509  std::list<DART_TYPE> swapped_edges;
510  SwapEdgesAwayFromInteriorNode<TRAITS_TYPE>( aDart, swapped_edges );
511 
512  // The reverse operation of split triangle:
513  // Make one triangle of the three triangles at the node associated with dart
514  // TRAITS_TYPE::
516 
517  // ???? Not generic yet if we are very strict:
518  // When calling unsplit triangle, darts at the three opposite sides may
519  // change!
520  // Should we hide them longer away??? This is possible since they cannot
521  // be boundary edges.
522  // ----> Or should we just require that they are not changed???
523 
524  // Make the swapped-away edges Delaunay.
525  // Note the theoretical result: if there are no edges in the list,
526  // the triangulation is Delaunay already
527 
528  OptimizeDelaunay<TRAITS_TYPE, DART_TYPE>( swapped_edges );
529 }
530 
532 
535 //------------------------------------------------------------------------------------------------
536 // Private/Hidden function (might change later)
537 template <class TOPOLOGY_ELEMENT_TYPE, class DART_TYPE>
538 bool TRIANGULATION_HELPER::isMemberOfFace( const TOPOLOGY_ELEMENT_TYPE& aTopologyElement,
539  const DART_TYPE& aDart )
540 {
541  // Check if the given topology element (node, edge or face) is a member of the face
542  // Assumes:
543  // - DART_TYPE::isMember(TOPOLOGY_ELEMENT_TYPE)
544 
545  DART_TYPE dart_iter = aDart;
546 
547  do
548  {
549  if( dart_iter.isMember( aTopologyElement ) )
550  return true;
551  dart_iter.Alpha0().Alpha1();
552  }
553  while( dart_iter != aDart );
554 
555  return false;
556 }
557 
558 //------------------------------------------------------------------------------------------------
559 // Private/Hidden function (might change later)
560 template <class TRAITS_TYPE, class NODE_TYPE, class DART_TYPE>
561 bool TRIANGULATION_HELPER::locateFaceWithNode( const NODE_TYPE& aNode, DART_TYPE& aDartIter )
562 {
563  // Locate a face in the topology structure with the given node as a member
564  // Assumes:
565  // - TRAITS_TYPE::Orient2D(DART_TYPE, DART_TYPE, NODE_TYPE)
566  // - DART_TYPE::isMember(NODE_TYPE)
567  // - Note that if false is returned, the node might still be in the
568  // topology structure. Application programmer
569  // should check all if by hypothesis the node is in the topology structure;
570  // see doc. on LocateTriangle.
571 
572  bool status = LocateFaceSimplest<TRAITS_TYPE>( aNode, aDartIter );
573 
574  if( status == false )
575  return status;
576 
577  // True was returned from LocateFaceSimplest, but if the located triangle is
578  // degenerate and the node is on the extension of the edges,
579  // the node might still be inside. Check if node is a member and return false
580  // if not. (Still the node might be in the topology structure, see doc. above
581  // and in locateTriangle(const POINT_TYPE& point, DART_TYPE& dart_iter)
582 
583  return isMemberOfFace( aNode, aDartIter );
584 }
585 
610 template <class TRAITS_TYPE, class POINT_TYPE, class DART_TYPE>
611 bool TRIANGULATION_HELPER::LocateFaceSimplest( const POINT_TYPE& aPoint, DART_TYPE& aDart )
612 {
613  // Not degenerate triangles if point is on the extension of the edges
614  // But inTriangle may be called in case of true (may update to inFace2)
615  // Convex boundary
616  // no holes
617  // convex faces (works for general convex faces)
618  // Not specialized for triangles, but ok?
619  //
620  // TRAITS_TYPE::orint2d(POINT_TYPE) is the half open half-plane defined
621  // by the dart:
622  // n1 = dart.node()
623  // n2 = dart.Alpha0().node
624  // Only the following gives true:
625  // ((n2->x()-n1->x())*(point.y()-n1->y()) >= (point.x()-n1->x())*(n2->y()-n1->y()))
626 
627  DART_TYPE dart_start;
628  dart_start = aDart;
629  DART_TYPE dart_prev;
630 
631  DART_TYPE d0;
632  for( ;; )
633  {
634  d0 = aDart;
635  d0.Alpha0();
636 
637  if( TRAITS_TYPE::Orient2D( aDart, d0, aPoint ) >= 0 )
638  {
639  aDart.Alpha0().Alpha1();
640  if( aDart == dart_start )
641  return true; // left to all edges in face
642  }
643  else
644  {
645  dart_prev = aDart;
646  aDart.Alpha2();
647 
648  if( aDart == dart_prev )
649  return false; // iteration to outside boundary
650 
651  dart_start = aDart;
652  dart_start.Alpha0();
653 
654  aDart.Alpha1(); // avoid twice on same edge and ccw in next
655  }
656  }
657 }
658 
659 
682 template <class TRAITS_TYPE, class POINT_TYPE, class DART_TYPE>
683 bool TRIANGULATION_HELPER::LocateTriangle( const POINT_TYPE& aPoint, DART_TYPE& aDart )
684 {
685  // The purpose is to have a fast and stable procedure that
686  // i) avoids concluding that a point is inside a triangle if it is not inside
687  // ii) avoids infinite loops
688 
689  // Thus, if false is returned, the point might still be inside a triangle in
690  // the triangulation. But this will probably only occur in the following cases:
691  // i) There are holes in the triangulation which causes the procedure to stop.
692  // ii) The boundary of the m_triangulation is not convex.
693  // ii) There might be degenerate triangles interior to the triangulation, or on the
694  // the boundary, which in some cases might cause the procedure to stop there due
695  // to the logic of the algorithm.
696 
697  // It is the application programmer's responsibility to check further if false is
698  // returned. For example, if by hypothesis the point is inside a triangle
699  // in the triangulation and and false is returned, then all triangles in the
700  // triangulation should be checked by the application. This can be done using
701  // the function:
702  // bool inTriangle(const POINT_TYPE& point, const DART_TYPE& dart).
703 
704  // Assumes:
705  // - CrossProduct2D, ScalarProduct2D etc., see functions called
706 
707  bool status = LocateFaceSimplest<TRAITS_TYPE>( aPoint, aDart );
708 
709  if( status == false )
710  return status;
711 
712  // There may be degeneracy, i.e., the point might be outside the triangle
713  // on the extension of the edges of a degenerate triangle.
714 
715  // The next call returns true if inside a non-degenerate or a degenerate triangle,
716  // but false if the point coincides with the "supernode" in the case where all
717  // edges are degenerate.
718  return InTriangle<TRAITS_TYPE>( aPoint, aDart );
719 }
720 
735 template <class TRAITS_TYPE, class POINT_TYPE, class DART_TYPE>
736 bool TRIANGULATION_HELPER::InTriangle( const POINT_TYPE& aPoint, const DART_TYPE& aDart )
737 {
738 
739  // SHOULD WE INCLUDE A STRATEGY WITH EDGE X e_1 ETC? TO GUARANTEE THAT
740  // ONLY ON ONE EDGE? BUT THIS DOES NOT SOLVE PROBLEMS WITH
741  // notInE1 && notInE1.neghbour ?
742 
743  // Returns true if inside (but not necessarily strictly inside)
744  // Works for degenerate triangles, but not when all edges are degenerate,
745  // and the aPoint coincides with all nodes;
746  // then false is always returned.
747 
748  typedef typename TRAITS_TYPE::REAL_TYPE REAL_TYPE;
749 
750  DART_TYPE dart_iter = aDart;
751 
752  REAL_TYPE cr1 = TRAITS_TYPE::CrossProduct2D( dart_iter, aPoint );
753  if( cr1 < 0 )
754  return false;
755 
756  dart_iter.Alpha0().Alpha1();
757  REAL_TYPE cr2 = TRAITS_TYPE::CrossProduct2D( dart_iter, aPoint );
758 
759  if( cr2 < 0 )
760  return false;
761 
762  dart_iter.Alpha0().Alpha1();
763  REAL_TYPE cr3 = TRAITS_TYPE::CrossProduct2D( dart_iter, aPoint );
764  if( cr3 < 0 )
765  return false;
766 
767  // All cross products are >= 0
768  // Check for degeneracy
769  if( cr1 != 0 || cr2 != 0 || cr3 != 0 )
770  return true; // inside non-degenerate face
771 
772  // All cross-products are zero, i.e. degenerate triangle, check if inside
773  // Strategy: d.ScalarProduct2D >= 0 && alpha0(d).d.ScalarProduct2D >= 0 for one of
774  // the edges. But if all edges are degenerate and the aPoint is on (all) the nodes,
775  // then "false is returned".
776 
777  DART_TYPE dart_tmp = dart_iter;
778  REAL_TYPE sc1 = TRAITS_TYPE::ScalarProduct2D( dart_tmp, aPoint );
779  REAL_TYPE sc2 = TRAITS_TYPE::ScalarProduct2D( dart_tmp.Alpha0(), aPoint );
780 
781  if( sc1 >= 0 && sc2 >= 0 )
782  {
783  // test for degenerate edge
784  if( sc1 != 0 || sc2 != 0 )
785  return true; // interior to this edge or on a node (but see comment above)
786  }
787 
788  dart_tmp = dart_iter.Alpha0().Alpha1();
789  sc1 = TRAITS_TYPE::ScalarProduct2D( dart_tmp, aPoint );
790  sc2 = TRAITS_TYPE::ScalarProduct2D( dart_tmp.Alpha0(), aPoint );
791 
792  if( sc1 >= 0 && sc2 >= 0 )
793  {
794  // test for degenerate edge
795  if( sc1 != 0 || sc2 != 0 )
796  return true; // interior to this edge or on a node (but see comment above)
797  }
798 
799  dart_tmp = dart_iter.Alpha1();
800  sc1 = TRAITS_TYPE::ScalarProduct2D( dart_tmp, aPoint );
801  sc2 = TRAITS_TYPE::ScalarProduct2D( dart_tmp.Alpha0(), aPoint );
802 
803  if( sc1 >= 0 && sc2 >= 0 )
804  {
805  // test for degenerate edge
806  if( sc1 != 0 || sc2 != 0 )
807  return true; // interior to this edge or on a node (but see comment above)
808  }
809 
810  // Not on any of the edges of the degenerate triangle.
811  // The only possibility for the aPoint to be "inside" is that all edges are degenerate
812  // and the point coincide with all nodes. So false is returned in this case.
813 
814  return false;
815 }
816 
817 
818  //------------------------------------------------------------------------------------------------
819 // Private/Hidden function (might change later)
820 template <class DART_TYPE>
821 void TRIANGULATION_HELPER::getAdjacentTriangles( const DART_TYPE& aDart, DART_TYPE& aT1,
822  DART_TYPE& aT2, DART_TYPE& aT3 )
823 {
824 
825  DART_TYPE dart_iter = aDart;
826 
827  // add first
828  if( dart_iter.Alpha2() != aDart )
829  {
830  aT1 = dart_iter;
831  dart_iter = aDart;
832  }
833 
834  // add second
835  dart_iter.Alpha0();
836  dart_iter.Alpha1();
837  DART_TYPE dart_prev = dart_iter;
838 
839  if( ( dart_iter.Alpha2() ) != dart_prev )
840  {
841  aT2 = dart_iter;
842  dart_iter = dart_prev;
843  }
844 
845  // add third
846  dart_iter.Alpha0();
847  dart_iter.Alpha1();
848  dart_prev = dart_iter;
849 
850  if( ( dart_iter.Alpha2() ) != dart_prev )
851  aT3 = dart_iter;
852 }
853 
854 //------------------------------------------------------------------------------------------------
870 template <class DART_TYPE, class DART_LIST_TYPE>
871 void TRIANGULATION_HELPER::GetBoundary( const DART_TYPE& aDart, DART_LIST_TYPE& aBoundary )
872 {
873  // assumes the given dart is at the boundary (by edge)
874 
875  DART_TYPE dart_iter( aDart );
876  aBoundary.push_back( dart_iter ); // Given dart as first element
877  dart_iter.Alpha0();
878  PositionAtNextBoundaryEdge( dart_iter );
879 
880  while( dart_iter != aDart )
881  {
882  aBoundary.push_back( dart_iter );
883  dart_iter.Alpha0();
884  PositionAtNextBoundaryEdge( dart_iter );
885  }
886 }
887 
900 template <class DART_TYPE>
901 bool TRIANGULATION_HELPER::IsBoundaryEdge( const DART_TYPE& aDart )
902 {
903  DART_TYPE dart_iter = aDart;
904 
905  if( dart_iter.Alpha2() == aDart )
906  return true;
907  else
908  return false;
909 }
910 
914 template <class DART_TYPE>
915 bool TRIANGULATION_HELPER::IsBoundaryFace( const DART_TYPE& aDart )
916 {
917  // Strategy: boundary if alpha2(d)=d
918 
919  DART_TYPE dart_iter( aDart );
920  DART_TYPE dart_prev;
921 
922  do
923  {
924  dart_prev = dart_iter;
925 
926  if( dart_iter.Alpha2() == dart_prev )
927  return true;
928  else
929  dart_iter = dart_prev; // back again
930 
931  dart_iter.Alpha0();
932  dart_iter.Alpha1();
933 
934  } while( dart_iter != aDart );
935 
936  return false;
937 }
938 
942 template <class DART_TYPE>
943 bool TRIANGULATION_HELPER::IsBoundaryNode( const DART_TYPE& aDart )
944 {
945  // Strategy: boundary if alpha2(d)=d
946 
947  DART_TYPE dart_iter( aDart );
948  DART_TYPE dart_prev;
949 
950  // If input dart is reached again, then internal node
951  // If alpha2(d)=d, then boundary
952 
953  do
954  {
955  dart_iter.Alpha1();
956  dart_prev = dart_iter;
957  dart_iter.Alpha2();
958 
959  if( dart_iter == dart_prev )
960  return true;
961 
962  } while( dart_iter != aDart );
963 
964  return false;
965 }
966 
974 template <class DART_TYPE>
975 int TRIANGULATION_HELPER::GetDegreeOfNode( const DART_TYPE& aDart )
976 {
977  DART_TYPE dart_iter( aDart );
978  DART_TYPE dart_prev;
979 
980  // If input dart is reached again, then interior node
981  // If alpha2(d)=d, then boundary
982 
983  int degree = 0;
984  bool boundaryVisited = false;
985  do
986  {
987  dart_iter.Alpha1();
988  degree++;
989  dart_prev = dart_iter;
990 
991  dart_iter.Alpha2();
992 
993  if( dart_iter == dart_prev )
994  {
995  if( !boundaryVisited )
996  {
997  boundaryVisited = true;
998  // boundary is reached first time, count in the reversed direction
999  degree++; // count the start since it is not done above
1000  dart_iter = aDart;
1001  dart_iter.Alpha2();
1002  } else
1003  return degree;
1004  }
1005 
1006  } while( dart_iter != aDart );
1007 
1008  return degree;
1009 }
1010 
1011 // Modification of GetDegreeOfNode:
1012 // Strategy, reverse the list and start in the other direction if the boundary
1013 // is reached. NB. copying of darts but ok., or we could have collected pointers,
1014 // but the memory management.
1015 
1016 // NOTE: not symmetry if we choose to collect opposite edges
1017 // now we collect darts with radiating edges
1018 
1019 // Remember that we must also copy the node, but ok with push_back
1020 // The size of the list will be the degree of the node
1021 
1022 // No CW/CCW since topology only
1023 
1024 // Each dart consists of an incident edge and an adjacent node.
1025 // But note that this is only how we interpret the dart in this implementation.
1026 // Given this list, how can we find the opposite edges:
1027 // We can perform alpha1 on each, but for boundary nodes we will get one edge twice.
1028 // But this is will always be the last dart!
1029 // The darts in the list are in sequence and starts with the alpha0(dart)
1030 // alpha0, alpha1 and alpha2
1031 
1032 // Private/Hidden function
1033 template <class DART_TYPE>
1034 void TRIANGULATION_HELPER::getNeighborNodes( const DART_TYPE& aDart,
1035  std::list<DART_TYPE>& aNodeList, bool& aBoundary )
1036 {
1037  DART_TYPE dart_iter( aDart );
1038  dart_iter.Alpha0(); // position the dart at an opposite node
1039 
1040  DART_TYPE dart_prev = dart_iter;
1041  bool start_at_boundary = false;
1042  dart_iter.Alpha2();
1043 
1044  if( dart_iter == dart_prev )
1045  start_at_boundary = true;
1046  else
1047  dart_iter = dart_prev; // back again
1048 
1049  DART_TYPE dart_start = dart_iter;
1050 
1051  do
1052  {
1053  aNodeList.push_back( dart_iter );
1054  dart_iter.Alpha1();
1055  dart_iter.Alpha0();
1056  dart_iter.Alpha1();
1057  dart_prev = dart_iter;
1058  dart_iter.Alpha2();
1059 
1060  if( dart_iter == dart_prev )
1061  {
1062  // boundary reached
1063  aBoundary = true;
1064 
1065  if( start_at_boundary == true )
1066  {
1067  // add the dart which now is positioned at the opposite boundary
1068  aNodeList.push_back( dart_iter );
1069  return;
1070  }
1071  else
1072  {
1073  // call the function again such that we start at the boundary
1074  // first clear the list and reposition to the initial node
1075  dart_iter.Alpha0();
1076  aNodeList.clear();
1077  getNeighborNodes( dart_iter, aNodeList, aBoundary );
1078 
1079  return; // after one recursive step
1080  }
1081  }
1082  }
1083  while( dart_iter != dart_start );
1084 
1085  aBoundary = false;
1086 }
1087 
1104 template <class DART_TYPE, class DART_LIST_TYPE>
1105 void TRIANGULATION_HELPER::Get0OrbitInterior( const DART_TYPE& aDart, DART_LIST_TYPE& aOrbit )
1106 {
1107  DART_TYPE d_iter = aDart;
1108  aOrbit.push_back( d_iter );
1109  d_iter.Alpha1().Alpha2();
1110 
1111  while( d_iter != aDart )
1112  {
1113  aOrbit.push_back( d_iter );
1114  d_iter.Alpha1().Alpha2();
1115  }
1116 }
1117 
1137 template <class DART_TYPE, class DART_LIST_TYPE>
1138 void TRIANGULATION_HELPER::Get0OrbitBoundary( const DART_TYPE& aDart, DART_LIST_TYPE& aOrbit )
1139 {
1140  DART_TYPE dart_prev;
1141  DART_TYPE d_iter = aDart;
1142 
1143  do
1144  {
1145  aOrbit.push_back( d_iter );
1146  d_iter.Alpha1();
1147  dart_prev = d_iter;
1148  d_iter.Alpha2();
1149  }
1150  while( d_iter != dart_prev );
1151 
1152  aOrbit.push_back( d_iter ); // the last one with opposite orientation
1153 }
1154 
1165 template <class DART_TYPE>
1166 bool TRIANGULATION_HELPER::Same0Orbit( const DART_TYPE& aD1, const DART_TYPE& aD2 )
1167 {
1168  // Two copies of the same dart
1169  DART_TYPE d_iter = aD2;
1170  DART_TYPE d_end = aD2;
1171 
1172  if( isBoundaryNode( d_iter ) )
1173  {
1174  // position at both boundary edges
1175  PositionAtNextBoundaryEdge( d_iter );
1176  d_end.Alpha1();
1177  PositionAtNextBoundaryEdge( d_end );
1178  }
1179 
1180  for( ;; )
1181  {
1182  if( d_iter == aD1 )
1183  return true;
1184 
1185  d_iter.Alpha1();
1186 
1187  if( d_iter == aD1 )
1188  return true;
1189 
1190  d_iter.Alpha2();
1191 
1192  if( d_iter == d_end )
1193  break;
1194  }
1195 
1196  return false;
1197 }
1198 
1203 template <class DART_TYPE>
1204 bool TRIANGULATION_HELPER::Same1Orbit( const DART_TYPE& aD1, const DART_TYPE& aD2 )
1205 {
1206  DART_TYPE d_iter = aD2;
1207 
1208  // (Also works at the boundary)
1209  return ( d_iter == aD1 || d_iter.Alpha0() == aD1 ||
1210  d_iter.Alpha2() == aD1 || d_iter.Alpha0() == aD1 );
1211 }
1212 
1213 //------------------------------------------------------------------------------------------------
1218 template <class DART_TYPE>
1219 bool TRIANGULATION_HELPER::Same2Orbit( const DART_TYPE& aD1, const DART_TYPE& aD2 )
1220 {
1221  DART_TYPE d_iter = aD2;
1222 
1223  return ( d_iter == aD1 || d_iter.Alpha0() == aD1 || d_iter.Alpha1() == aD1 ||
1224  d_iter.Alpha0() == aD1 || d_iter.Alpha1() == aD1 || d_iter.Alpha0() == aD1 );
1225 }
1226 
1227 // Private/Hidden function
1228 template <class TRAITS_TYPE, class DART_TYPE>
1229 bool TRIANGULATION_HELPER::degenerateTriangle( const DART_TYPE& aDart )
1230 {
1231  // Check if triangle is degenerate
1232  // Assumes CCW dart
1233 
1234  DART_TYPE d1 = aDart;
1235  DART_TYPE d2 = d1;
1236  d2.Alpha1();
1237 
1238  return ( TRAITS_TYPE::CrossProduct2D( d1, d2 ) == 0 );
1239 }
1240 
1252 template <class TRAITS_TYPE, class DART_TYPE>
1253 bool TRIANGULATION_HELPER::SwappableEdge( const DART_TYPE& aDart, bool aAllowDegeneracy )
1254 {
1255  // How "safe" is it?
1256 
1257  if( IsBoundaryEdge( aDart ) )
1258  return false;
1259 
1260  // "angles" are at the diagonal
1261  DART_TYPE d1 = aDart;
1262  d1.Alpha2().Alpha1();
1263  DART_TYPE d2 = aDart;
1264  d2.Alpha1();
1265 
1266  if( aAllowDegeneracy )
1267  {
1268  if( TRAITS_TYPE::CrossProduct2D( d1, d2 ) < 0.0 )
1269  return false;
1270  }
1271  else
1272  {
1273  if( TRAITS_TYPE::CrossProduct2D( d1, d2 ) <= 0.0 )
1274  return false;
1275  }
1276 
1277  // Opposite side (still angle at the diagonal)
1278  d1 = aDart;
1279  d1.Alpha0();
1280  d2 = d1;
1281  d1.Alpha1();
1282  d2.Alpha2().Alpha1();
1283 
1284  if( aAllowDegeneracy )
1285  {
1286  if( TRAITS_TYPE::CrossProduct2D( d1, d2 ) < 0.0 )
1287  return false;
1288  }
1289  else
1290  {
1291  if( TRAITS_TYPE::CrossProduct2D( d1, d2 ) <= 0.0 )
1292  return false;
1293  }
1294 
1295  return true;
1296 }
1297 
1309 template <class DART_TYPE>
1311 {
1312  DART_TYPE dart_prev;
1313 
1314  // If alpha2(d)=d, then boundary
1315 
1316  //old convention: dart.Alpha0();
1317  do
1318  {
1319  aDart.Alpha1();
1320  dart_prev = aDart;
1321  aDart.Alpha2();
1322  }
1323  while( aDart != dart_prev );
1324 }
1325 
1334 template <class TRAITS_TYPE, class DART_TYPE>
1335 bool TRIANGULATION_HELPER::ConvexBoundary( const DART_TYPE& aDart )
1336 {
1337  std::list<DART_TYPE> blist;
1338  getBoundary( aDart, blist );
1339 
1340  int no;
1341  no = (int) blist.size();
1342  typename std::list<DART_TYPE>::const_iterator bit = blist.begin();
1343  DART_TYPE d1 = *bit;
1344  ++bit;
1345  DART_TYPE d2;
1346  bool convex = true;
1347 
1348  for( ; bit != blist.end(); ++bit )
1349  {
1350  d2 = *bit;
1351  double crossProd = TRAITS_TYPE::CrossProduct2D( d1, d2 );
1352 
1353  if( crossProd < 0.0 )
1354  {
1355  //cout << "!!! Boundary is NOT convex: crossProd = " << crossProd << endl;
1356  convex = false;
1357  return convex;
1358  }
1359 
1360  d1 = d2;
1361  }
1362 
1363  // Check the last angle
1364  d2 = *blist.begin();
1365  double crossProd = TRAITS_TYPE::CrossProduct2D( d1, d2 );
1366 
1367  if( crossProd < 0.0 )
1368  {
1369  //cout << "!!! Boundary is NOT convex: crossProd = " << crossProd << endl;
1370  convex = false;
1371  }
1372 
1373  //if (convex)
1374  // cout << "\n---> Boundary is convex\n" << endl;
1375  //cout << endl;
1376  return convex;
1377 }
1378 
1380 
1383 //------------------------------------------------------------------------------------------------
1401 template <class TRAITS_TYPE, class DART_TYPE, class DART_LIST_TYPE>
1402 void TRIANGULATION_HELPER::OptimizeDelaunay( DART_LIST_TYPE& aElist )
1403 {
1404  OptimizeDelaunay<TRAITS_TYPE, DART_TYPE, DART_LIST_TYPE>( aElist, aElist.end() );
1405 }
1406 
1407 //------------------------------------------------------------------------------------------------
1408 template <class TRAITS_TYPE, class DART_TYPE, class DART_LIST_TYPE>
1409 void TRIANGULATION_HELPER::OptimizeDelaunay( DART_LIST_TYPE& aElist,
1410  const typename DART_LIST_TYPE::iterator aEnd )
1411 {
1412  // CCW darts
1413  // Optimize here means Delaunay, but could be any criterion by
1414  // requiring a "should swap" in the traits class, or give
1415  // a function object?
1416  // Assumes that elist has only one dart for each arc.
1417  // Darts outside the quadrilateral are preserved
1418 
1419  // For some data structures it is possible to preserve
1420  // all darts when swapping. Thus a preserve_darts_when swapping
1421  // ccould be given to indicate this and we would gain performance by avoiding
1422  // find in list.
1423 
1424  // Requires that swap retuns a dart in the "same position when rotated CCW"
1425  // (A vector instead of a list may be better.)
1426 
1427  // First check that elist is not empty
1428  if( aElist.empty() )
1429  return;
1430 
1431  // Avoid cycling by more extensive circumcircle test
1432  bool cycling_check = true;
1433  bool optimal = false;
1434  typename DART_LIST_TYPE::iterator it;
1435 
1436  typename DART_LIST_TYPE::iterator end_opt = aEnd;
1437 
1438  // Hmm... The following code is trying to derefence an iterator that may
1439  // be invalid. This may lead to debug error on Windows, so we comment out
1440  // this code. Checking elist.empty() above will prevent some
1441  // problems...
1442  //
1443  // last_opt is passed the end of the "active list"
1444  //typename DART_LIST_TYPE::iterator end_opt;
1445  //if (*end != NULL)
1446  // end_opt = end;
1447  //else
1448  // end_opt = elist.end();
1449 
1450  while( !optimal )
1451  {
1452  optimal = true;
1453  for( it = aElist.begin(); it != end_opt; ++it )
1454  {
1455  if( SwapTestDelaunay<TRAITS_TYPE>( *it, cycling_check ) )
1456  {
1457  // Preserve darts. Potential darts in the list are:
1458  // - The current dart
1459  // - the four CCW darts on the boundary of the quadrilateral
1460  // (the current arc has only one dart)
1461 
1462  SwapEdgeInList<TRAITS_TYPE, DART_TYPE>( it, aElist );
1463 
1464  optimal = false;
1465  } // end if should swap
1466  } // end for
1467  } // end pass
1468 }
1469 
1483 template <class TRAITS_TYPE, class DART_TYPE>
1484 #if ((_MSC_VER > 0) && (_MSC_VER < 1300))//#ifdef _MSC_VER
1485 bool TRIANGULATION_HELPER::SwapTestDelaunay(const DART_TYPE& aDart, bool aCyclingCheck = false) const
1486 {
1487 #else
1488 bool TRIANGULATION_HELPER::SwapTestDelaunay( const DART_TYPE& aDart, bool aCyclingCheck ) const
1489 {
1490 #endif
1491  // The general strategy is taken from Cline & Renka. They claim that
1492  // their algorithm insure numerical stability, but experiments show
1493  // that this is not correct for neutral, or almost neutral cases.
1494  // I have extended this strategy (without using tolerances) to avoid
1495  // cycling and infinit loops when used in connection with LOP algorithms;
1496  // see the comments below.
1497 
1498  typedef typename TRAITS_TYPE::REAL_TYPE REAL_TYPE;
1499 
1500  if( IsBoundaryEdge( aDart ) )
1501  return false;
1502 
1503  DART_TYPE v11 = aDart;
1504  v11.Alpha1().Alpha0();
1505  DART_TYPE v12 = v11;
1506  v12.Alpha1();
1507 
1508  DART_TYPE v22 = aDart;
1509  v22.Alpha2().Alpha1().Alpha0();
1510  DART_TYPE v21 = v22;
1511  v21.Alpha1();
1512 
1513  REAL_TYPE cos1 = TRAITS_TYPE::ScalarProduct2D( v11, v12 );
1514  REAL_TYPE cos2 = TRAITS_TYPE::ScalarProduct2D( v21, v22 );
1515 
1516  // "Angles" are opposite to the diagonal.
1517  // The diagonals should be swapped iff (t1+t2) .gt. 180
1518  // degrees. The following two tests insure numerical
1519  // stability according to Cline & Renka. But experiments show
1520  // that cycling may still happen; see the aditional test below.
1521  if( cos1 >= 0 && cos2 >= 0 ) // both angles are grater or equual 90
1522  return false;
1523 
1524  if( cos1 < 0 && cos2 < 0 ) // both angles are less than 90
1525  return true;
1526 
1527  REAL_TYPE sin1 = TRAITS_TYPE::CrossProduct2D( v11, v12 );
1528  REAL_TYPE sin2 = TRAITS_TYPE::CrossProduct2D( v21, v22 );
1529  REAL_TYPE sin12 = sin1 * cos2 + cos1 * sin2;
1530 
1531  if( sin12 >= 0 ) // equality represents a neutral case
1532  return false;
1533 
1534  if( aCyclingCheck )
1535  {
1536  // situation so far is sin12 < 0. Test if this also
1537  // happens for the swapped edge.
1538 
1539  // The numerical calculations so far indicate that the edge is
1540  // not Delaunay and should not be swapped. But experiments show that
1541  // in neutral cases, or almost neutral cases, it may happen that
1542  // the swapped edge may again be found to be not Delaunay and thus
1543  // be swapped if we return true here. This may lead to cycling and
1544  // an infinte loop when used, e.g., in connection with OptimizeDelaunay.
1545  //
1546  // In an attempt to avoid this we test if the swapped edge will
1547  // also be found to be not Delaunay by repeating the last test above
1548  // for the swapped edge.
1549  // We now rely on the general requirement for TRAITS_TYPE::swapEdge which
1550  // should deliver CCW dart back in "the same position"; see the general
1551  // description. This will insure numerical stability as the next calculation
1552  // is the same as if this function was called again with the swapped edge.
1553  // Cycling is thus impossible provided that the initial tests above does
1554  // not result in ambiguity (and they should probably not do so).
1555 
1556  v11.Alpha0();
1557  v12.Alpha0();
1558  v21.Alpha0();
1559  v22.Alpha0();
1560  // as if the edge was swapped/rotated CCW
1561  cos1 = TRAITS_TYPE::ScalarProduct2D( v22, v11 );
1562  cos2 = TRAITS_TYPE::ScalarProduct2D( v12, v21 );
1563  sin1 = TRAITS_TYPE::CrossProduct2D( v22, v11 );
1564  sin2 = TRAITS_TYPE::CrossProduct2D( v12, v21 );
1565  sin12 = sin1 * cos2 + cos1 * sin2;
1566 
1567  if( sin12 < 0 )
1568  {
1569  // A neutral case, but the tests above lead to swapping
1570  return false;
1571  }
1572  }
1573 
1574  return true;
1575 }
1576 
1577 //-----------------------------------------------------------------------
1578 //
1579 // x
1580 //" / \ "
1581 // / | \ Darts:
1582 //oe2 / | \ oe2 = oppEdge2
1583 // x....|....x
1584 // \ d| d/ d = diagonal (input and output)
1585 // \ | /
1586 // oe1 \ / oe1 = oppEdge1
1587 // x
1588 //
1589 //-----------------------------------------------------------------------
1603 template <class TRAITS_TYPE, class DART_TYPE>
1604 void TRIANGULATION_HELPER::RecSwapDelaunay( DART_TYPE& aDiagonal )
1605 {
1606  if( !SwapTestDelaunay<TRAITS_TYPE>( aDiagonal ) )
1607  // ??? swapTestDelaunay also checks if boundary, so this can be optimized
1608  return;
1609 
1610  // Get the other "edges" of the current triangle; see illustration above.
1611  DART_TYPE oppEdge1 = aDiagonal;
1612  oppEdge1.Alpha1();
1613  bool b1;
1614 
1615  if( IsBoundaryEdge( oppEdge1 ) )
1616  b1 = true;
1617  else
1618  {
1619  b1 = false;
1620  oppEdge1.Alpha2();
1621  }
1622 
1623  DART_TYPE oppEdge2 = aDiagonal;
1624  oppEdge2.Alpha0().Alpha1().Alpha0();
1625  bool b2;
1626 
1627  if( IsBoundaryEdge( oppEdge2 ) )
1628  b2 = true;
1629  else
1630  {
1631  b2 = false;
1632  oppEdge2.Alpha2();
1633  }
1634 
1635  // Swap the given diagonal
1636  m_triangulation.swapEdge( aDiagonal );
1637 
1638  if( !b1 )
1639  RecSwapDelaunay<TRAITS_TYPE>( oppEdge1 );
1640 
1641  if( !b2 )
1642  RecSwapDelaunay<TRAITS_TYPE>( oppEdge2 );
1643 }
1644 
1670 template <class TRAITS_TYPE, class DART_TYPE, class LIST_TYPE>
1672  LIST_TYPE& aSwappedEdges )
1673 {
1674 
1675  // Same iteration as in fixEdgesAtCorner, but not boundary
1676  DART_TYPE dnext = aDart;
1677 
1678  // Allow degeneracy, otherwise we might end up with degree=4.
1679  // For example, the reverse operation of inserting a point on an
1680  // existing edge gives a situation where all edges are non-swappable.
1681  // Ideally, degeneracy in this case should be along the actual node,
1682  // but there is no strategy for this now.
1683  // ??? An alternative here is to wait with degeneracy till we get an
1684  // infinite loop with degree > 3.
1685  bool allowDegeneracy = true;
1686 
1687  int degree = getDegreeOfNode( aDart );
1688  DART_TYPE d_iter;
1689 
1690  while( degree > 3 )
1691  {
1692  d_iter = dnext;
1693  dnext.Alpha1().Alpha2();
1694 
1695  if( SwappableEdge<TRAITS_TYPE>( d_iter, allowDegeneracy ) )
1696  {
1697  m_triangulation.swapEdge( d_iter ); // swap the edge away
1698  // Collect swapped edges in the list
1699  // "Hide" the dart on the other side of the edge to avoid it being changed for
1700  // other swaps
1701  DART_TYPE swapped_edge = d_iter; // it was delivered back
1702  swapped_edge.Alpha2().Alpha0(); // CCW (if not at boundary)
1703  aSwappedEdges.push_back( swapped_edge );
1704 
1705  degree--;
1706  }
1707  }
1708 
1709  // Output, incident to the node
1710  aDart = dnext;
1711 }
1712 
1732 template <class TRAITS_TYPE, class DART_TYPE, class LIST_TYPE>
1734  LIST_TYPE& aSwappedEdges )
1735 {
1736  // All darts that are swappable.
1737  // To treat collinear nodes at an existing boundary, we must allow degeneracy
1738  // when swapping to the boundary.
1739  // dart is CCW and at the boundary.
1740  // The 0-orbit runs CCW
1741  // Deliver the dart back in the "same position".
1742  // Assume for the swap in the traits class:
1743  // - A dart on the swapped edge is delivered back in a position as
1744  // seen if it was glued to the edge when swapping (rotating) the edge CCW
1745 
1746  //int degree = getDegreeOfNode(dart);
1747 
1748  passes:
1749  // Swap swappable edges that radiate from the node away
1750  DART_TYPE d_iter = aDart; // ???? can simply use dart
1751  d_iter.Alpha1().Alpha2(); // first not at boundary
1752  DART_TYPE d_next = d_iter;
1753  bool bend = false;
1754  bool swapped_next_to_boundary = false;
1755  bool swapped_in_pass = false;
1756 
1757  bool allowDegeneracy; // = true;
1758  DART_TYPE tmp1, tmp2;
1759 
1760  while( !bend )
1761  {
1762  d_next.Alpha1().Alpha2();
1763 
1764  if( IsBoundaryEdge( d_next ) )
1765  bend = true; // then it is CW since alpha2
1766 
1767  // To allow removing among collinear nodes at the boundary,
1768  // degenerate triangles must be allowed
1769  // (they will be removed when used in connection with RemoveBoundaryNode)
1770  tmp1 = d_iter;
1771  tmp1.Alpha1();
1772  tmp2 = d_iter;
1773  tmp2.Alpha2().Alpha1(); // don't bother with boundary (checked later)
1774 
1775  if( IsBoundaryEdge( tmp1 ) && IsBoundaryEdge( tmp2 ) )
1776  allowDegeneracy = true;
1777  else
1778  allowDegeneracy = false;
1779 
1780  if( SwappableEdge<TRAITS_TYPE>( d_iter, allowDegeneracy ) )
1781  {
1782  m_triangulation.swapEdge( d_iter );
1783 
1784  // Collect swapped edges in the list
1785  // "Hide" the dart on the other side of the edge to avoid it being changed for
1786  // other swapps
1787  DART_TYPE swapped_edge = d_iter; // it was delivered back
1788  swapped_edge.Alpha2().Alpha0(); // CCW
1789  aSwappedEdges.push_back( swapped_edge );
1790 
1791  //degree--; // if degree is 2, or bend=true, we are done
1792  swapped_in_pass = true;
1793  if( bend )
1794  swapped_next_to_boundary = true;
1795  }
1796 
1797  if( !bend )
1798  d_iter = d_next;
1799  }
1800 
1801  // Deliver a dart as output in the same position as the incoming dart
1802  if( swapped_next_to_boundary )
1803  {
1804  // Assume that "swapping is CCW and dart is preserved in the same position
1805  d_iter.Alpha1().Alpha0().Alpha1(); // CW and see below
1806  }
1807  else
1808  {
1809  d_iter.Alpha1(); // CW and see below
1810  }
1811  PositionAtNextBoundaryEdge( d_iter ); // CCW
1812 
1813  aDart = d_iter; // for next pass or output
1814 
1815  // If a dart was swapped in this iteration we must run it more
1816  if( swapped_in_pass )
1817  goto passes;
1818 }
1819 
1826 template <class TRAITS_TYPE, class DART_TYPE, class DART_LIST_TYPE>
1827 void TRIANGULATION_HELPER::SwapEdgeInList( const typename DART_LIST_TYPE::iterator& aIt,
1828  DART_LIST_TYPE& aElist )
1829 {
1830 
1831  typename DART_LIST_TYPE::iterator it1, it2, it3, it4;
1832  DART_TYPE dart( *aIt );
1833 
1834  //typename TRAITS_TYPE::DART_TYPE d1 = dart; d1.Alpha2().Alpha1();
1835  //typename TRAITS_TYPE::DART_TYPE d2 = d1; d2.Alpha0().Alpha1();
1836  //typename TRAITS_TYPE::DART_TYPE d3 = dart; d3.Alpha0().Alpha1();
1837  //typename TRAITS_TYPE::DART_TYPE d4 = d3; d4.Alpha0().Alpha1();
1838  DART_TYPE d1 = dart;
1839  d1.Alpha2().Alpha1();
1840  DART_TYPE d2 = d1;
1841  d2.Alpha0().Alpha1();
1842  DART_TYPE d3 = dart;
1843  d3.Alpha0().Alpha1();
1844  DART_TYPE d4 = d3;
1845  d4.Alpha0().Alpha1();
1846 
1847  // Find pinters to the darts that may change.
1848  // ??? Note, this is not very efficient since we must use find, which is O(N),
1849  // four times.
1850  // - Solution?: replace elist with a vector of pair (dart,number)
1851  // and avoid find?
1852  // - make a function for swapping generically?
1853  // - sould we use another container type or,
1854  // - erase them and reinsert?
1855  // - or use two lists?
1856  it1 = find( aElist.begin(), aElist.end(), d1 );
1857  it2 = find( aElist.begin(), aElist.end(), d2 );
1858  it3 = find( aElist.begin(), aElist.end(), d3 );
1859  it4 = find( aElist.begin(), aElist.end(), d4 );
1860 
1861  m_triangulation.swapEdge( dart );
1862  // Update the current dart which may have changed
1863  *aIt = dart;
1864 
1865  // Update darts that may have changed again (if they were present)
1866  // Note that dart is delivered back after swapping
1867  if( it1 != aElist.end() )
1868  {
1869  d1 = dart;
1870  d1.Alpha1().Alpha0();
1871  *it1 = d1;
1872  }
1873 
1874  if( it2 != aElist.end() )
1875  {
1876  d2 = dart;
1877  d2.Alpha2().Alpha1();
1878  *it2 = d2;
1879  }
1880 
1881  if( it3 != aElist.end() )
1882  {
1883  d3 = dart;
1884  d3.Alpha2().Alpha1().Alpha0().Alpha1();
1885  *it3 = d3;
1886  }
1887 
1888  if( it4 != aElist.end() )
1889  {
1890  d4 = dart;
1891  d4.Alpha0().Alpha1();
1892  *it4 = d4;
1893  }
1894 }
1895 
1897 
1898 }
1899 // End of ttl namespace scope (but other files may also contain functions for ttl)
1900 
1901 #endif // _TTL_H_
static bool LocateFaceSimplest(const POINT_TYPE &aPoint, DART_TYPE &aDart)
Locates the face containing a given point.
Definition: ttl.h:611
void removeBoundaryTriangle(DART &aDart)
Removes a triangle with an edge at the boundary of the triangulation in the actual data structure...
Definition: hetriang.cpp:355
bool SwapTestDelaunay(const DART_TYPE &aDart, bool aCyclingCheck=false) const
Checks if the edge associated with dart should be swapped according to the Delaunay criterion...
Definition: ttl.h:1488
static bool Same0Orbit(const DART_TYPE &aD1, const DART_TYPE &aD2)
Checks if the two darts belong to the same 0-orbit, i.e., if they share a node.
Definition: ttl.h:1166
void splitTriangle(DART &aDart, const NODE_PTR &aPoint)
Splits the triangle associated with dart in the actual data structure into three new triangles joinin...
Definition: hetriang.cpp:342
static bool ConvexBoundary(const DART_TYPE &aDart)
Checks if the boundary of a triangulation is convex.
Definition: ttl.h:1335
void swapEdge(DART &aDart)
Swaps the edge associated with dart in the actual data structure.
Definition: hetriang.cpp:336
void OptimizeDelaunay(DART_LIST_TYPE &aElist)
Optimizes the edges in the given sequence according to the Delaunay criterion, i.e., such that the edge will fullfill the circumcircle criterion (or equivalently the MaxMin angle criterion) with respect to the quadrilaterals where they are diagonals.
Definition: ttl.h:1402
void RecSwapDelaunay(DART_TYPE &aDiagonal)
Recursively swaps edges in the triangulation according to the Delaunay criterion. ...
Definition: ttl.h:1604
static bool degenerateTriangle(const DART_TYPE &aDart)
Definition: ttl.h:1229
static bool IsBoundaryNode(const DART_TYPE &aDart)
Checks if the node associated with dart is at the boundary of the m_triangulation.
Definition: ttl.h:943
static bool InTriangle(const POINT_TYPE &aPoint, const DART_TYPE &aDart)
Checks if point is inside the triangle associated with dart.
Definition: ttl.h:736
static bool Same2Orbit(const DART_TYPE &aD1, const DART_TYPE &aD2)
Checks if the two darts belong to the same 2-orbit, i.e., if they lie in the same triangle...
Definition: ttl.h:1219
void RemoveNode(DART_TYPE &aDart)
Removes the node associated with dart and updates the triangulation to be Delaunay.
Definition: ttl.h:407
void SwapEdgesAwayFromBoundaryNode(DART_TYPE &aDart, LIST_TYPE &aSwappedEdges)
Swaps edges away from the (boundary) node associated with dart in such a way that when removing the e...
Definition: ttl.h:1733
void SwapEdgesAwayFromInteriorNode(DART_TYPE &aDart, LIST_TYPE &aSwappedEdges)
Swaps edges away from the (interior) node associated with dart such that that exactly three edges rem...
Definition: ttl.h:1671
void SwapEdgeInList(const typename DART_LIST_TYPE::iterator &aIt, DART_LIST_TYPE &aElist)
Swap the the edge associated with iterator it and update affected darts in elist accordingly.
Definition: ttl.h:1827
static void GetBoundary(const DART_TYPE &aDart, DART_LIST_TYPE &aBoundary)
Gets the boundary as sequence of darts, where the edges associated with the darts are boundary edges...
Definition: ttl.h:871
static void Get0OrbitBoundary(const DART_TYPE &aDart, DART_LIST_TYPE &aOrbit)
Gets the 0-orbit around a node at the boundary.
Definition: ttl.h:1138
static bool IsBoundaryEdge(const DART_TYPE &aDart)
Checks if the edge associated with dart is at the boundary of the m_triangulation.
Definition: ttl.h:901
static bool IsBoundaryFace(const DART_TYPE &aDart)
Checks if the face associated with dart is at the boundary of the m_triangulation.
Definition: ttl.h:915
void insertNodes(FORWARD_ITERATOR aFirst, FORWARD_ITERATOR aLast, DART_TYPE &aDart)
Definition: ttl.h:347
Main interface to TTL.
Definition: hetriang.h:59
hed::TRIANGULATION & m_triangulation
Definition: ttl.h:220
void reverseSplitTriangle(DART &aDart)
The reverse operation of TTLtraits::splitTriangle.
Definition: hetriang.cpp:349
static bool SwappableEdge(const DART_TYPE &aDart, bool aAllowDegeneracy=false)
Checks if the edge associated with dart is swappable, i.e., if the edge is a diagonal in a strictly c...
Definition: ttl.h:1253
static void getAdjacentTriangles(const DART_TYPE &aDart, DART_TYPE &aT1, DART_TYPE &aT2, DART_TYPE &aT3)
Definition: ttl.h:821
static bool Same1Orbit(const DART_TYPE &aD1, const DART_TYPE &aD2)
Checks if the two darts belong to the same 1-orbit, i.e., if they share an edge.
Definition: ttl.h:1204
static bool LocateTriangle(const POINT_TYPE &aPoint, DART_TYPE &aDart)
Locates the triangle containing a given point.
Definition: ttl.h:683
REAL_TYPE ScalarProduct2D(REAL_TYPE aDX1, REAL_TYPE aDY1, REAL_TYPE aDX2, REAL_TYPE aDY2)
Scalar product between two 2D vectors.
Definition: ttl_util.h:89
void RemoveInteriorNode(DART_TYPE &aDart)
Removes the interior node associated with dart and updates the triangulation to be Delaunay...
Definition: ttl.h:493
void RemoveBoundaryNode(DART_TYPE &aDart)
Removes the boundary node associated with dart and updates the triangulation to be Delaunay...
Definition: ttl.h:427
static DART_TYPE InsertConstraint(DART_TYPE &aDStart, DART_TYPE &aDEnd, bool aOptimizeDelaunay)
static long bit[8][9]
Definition: solve.cpp:826
static void PositionAtNextBoundaryEdge(DART_TYPE &aDart)
Given a dart, CCW or CW, positioned in a 0-orbit at the boundary of a tessellation.
Definition: ttl.h:1310
TRIANGULATION_HELPER(hed::TRIANGULATION &aTriang)
Definition: ttl.h:120
static void getNeighborNodes(const DART_TYPE &aDart, std::list< DART_TYPE > &aNodeList, bool &aBoundary)
Definition: ttl.h:1034
static void Get0OrbitInterior(const DART_TYPE &aDart, DART_LIST_TYPE &aOrbit)
Gets the 0-orbit around an interior node.
Definition: ttl.h:1105
bool InsertNode(DART_TYPE &aDart, POINT_TYPE &aPoint)
Inserts a new node in an existing Delaunay triangulation and swaps edges to obtain a new Delaunay tri...
Definition: ttl.h:289
void RemoveRectangularBoundary(DART_TYPE &aDart)
Removes the rectangular boundary of a triangulation as a final step of an incremental Delaunay triang...
Definition: ttl.h:379
static bool isMemberOfFace(const TOPOLOGY_ELEMENT_TYPE &aTopologyElement, const DART_TYPE &aDart)
Definition: ttl.h:538
REAL_TYPE CrossProduct2D(REAL_TYPE aDX1, REAL_TYPE aDY1, REAL_TYPE aDX2, REAL_TYPE aDY2)
Cross product between two 2D vectors.
Definition: ttl_util.h:102
static bool locateFaceWithNode(const NODE_TYPE &aNode, DART_TYPE &aDartIter)
Definition: ttl.h:561
static int GetDegreeOfNode(const DART_TYPE &aDart)
Returns the degree of the node associated with dart.
Definition: ttl.h:975
Triangulation class for the half-edge data structure with adaption to TTL.
Definition: hetriang.h:384