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