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:346
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:333
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:327
void OptimizeDelaunay(DART_LIST_TYPE &aElist)
Optimizes the edges in the given sequence according to the Delaunay criterion, i.e....
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
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 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:340
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)
size_t i
Definition: json11.cpp:597
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