QGIS API Documentation  3.21.0-Master (56b4176581)
qgsdualedgetriangulation.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  QgsDualEdgeTriangulation.cpp
3  -------------------------
4  copyright : (C) 2004 by Marco Hugentobler
5  email : [email protected]
6  ***************************************************************************/
7 
8 /***************************************************************************
9  * *
10  * This program is free software; you can redistribute it and/or modify *
11  * it under the terms of the GNU General Public License as published by *
12  * the Free Software Foundation; either version 2 of the License, or *
13  * (at your option) any later version. *
14  * *
15  ***************************************************************************/
16 
17 #include <QStack>
19 #include <map>
20 #include "MathUtils.h"
21 #include "qgsgeometry.h"
22 #include "qgslogger.h"
23 #include "qgsvectorfilewriter.h"
24 #include "qgsinterpolator.h"
25 
26 double leftOfTresh = 0;
27 
28 static bool inCircle( const QgsPoint &testedPoint, const QgsPoint &point1, const QgsPoint &point2, const QgsPoint &point3 )
29 {
30  const double x2 = point2.x() - point1.x();
31  const double y2 = point2.y() - point1.y();
32  const double x3 = point3.x() - point1.x();
33  const double y3 = point3.y() - point1.y();
34 
35  const double denom = x2 * y3 - y2 * x3;
36  double frac;
37 
38  if ( denom == 0 )
39  return false;
40 
41  frac = ( x2 * ( x2 - x3 ) + y2 * ( y2 - y3 ) ) / denom;
42  const double cx = ( x3 + frac * y3 ) / 2;
43  const double cy = ( y3 - frac * x3 ) / 2;
44  const double squaredRadius = ( cx * cx + cy * cy );
45  const QgsPoint center( cx + point1.x(), cy + point1.y() );
46 
47  return center.distanceSquared( testedPoint ) < squaredRadius ;
48 }
49 
51 {
52  //remove all the points
53  if ( !mPointVector.isEmpty() )
54  {
55  for ( int i = 0; i < mPointVector.count(); i++ )
56  {
57  delete mPointVector[i];
58  }
59  }
60 
61  //remove all the HalfEdge
62  if ( !mHalfEdge.isEmpty() )
63  {
64  for ( int i = 0; i < mHalfEdge.count(); i++ )
65  {
66  delete mHalfEdge[i];
67  }
68  }
69 }
70 
72 {
73  QgsDebugMsg( QStringLiteral( "performing consistency test" ) );
74 
75  for ( int i = 0; i < mHalfEdge.count(); i++ )
76  {
77  const int a = mHalfEdge[mHalfEdge[i]->getDual()]->getDual();
78  const int b = mHalfEdge[mHalfEdge[mHalfEdge[i]->getNext()]->getNext()]->getNext();
79  if ( i != a )
80  {
81  QgsDebugMsg( QStringLiteral( "warning, first test failed" ) );
82  }
83  if ( i != b )
84  {
85  QgsDebugMsg( QStringLiteral( "warning, second test failed" ) );
86  }
87  }
88  QgsDebugMsg( QStringLiteral( "consistency test finished" ) );
89 }
90 
91 void QgsDualEdgeTriangulation::addLine( const QVector<QgsPoint> &points, QgsInterpolator::SourceType lineType )
92 {
93  int actpoint = -10;//number of the last point, which has been inserted from the line
94  int currentpoint = -10;//number of the point, which is currently inserted from the line
95 
96  int i = 0;
97  for ( const QgsPoint &point : points )
98  {
99  actpoint = addPoint( point );
100  i++;
101  if ( actpoint != -100 )
102  {
103  break;
104  }
105  }
106 
107  if ( actpoint == -100 )//no point of the line could be inserted
108  {
109  return;
110  }
111 
112  for ( ; i < points.size(); ++i )
113  {
114  currentpoint = addPoint( points.at( i ) );
115  if ( currentpoint != -100 && actpoint != -100 && currentpoint != actpoint )//-100 is the return value if the point could not be not inserted
116  {
117  insertForcedSegment( actpoint, currentpoint, lineType );
118  }
119  actpoint = currentpoint;
120  }
121 }
122 
124 {
125  //first update the bounding box
126  if ( mPointVector.isEmpty() )//update bounding box when the first point is inserted
127  {
128  mXMin = p.x();
129  mYMin = p.y();
130  mXMax = p.x();
131  mYMax = p.y();
132  }
133  else //update bounding box else
134  {
135  mXMin = std::min( p.x(), mXMin );
136  mXMax = std::max( p.x(), mXMax );
137  mYMin = std::min( p.y(), mYMin );
138  mYMax = std::max( p.y(), mYMax );
139  }
140 
141  //then update mPointVector
142  mPointVector.append( new QgsPoint( p ) );
143 
144  //then update the HalfEdgeStructure
145  if ( mDimension == -1 )//insert the first point into the triangulation
146  {
147  const unsigned int zedge /* 0 */ = insertEdge( -10, -10, -1, false, false ); //edge pointing from p to the virtual point
148  const unsigned int fedge /* 1 */ = insertEdge( static_cast<int>( zedge ), static_cast<int>( zedge ), 0, false, false ); //edge pointing from the virtual point to p
149  ( mHalfEdge.at( zedge ) )->setDual( static_cast<int>( fedge ) );
150  ( mHalfEdge.at( zedge ) )->setNext( static_cast<int>( fedge ) );
151  mDimension = 0;
152  }
153  else if ( mDimension == 0 )//insert the second point into the triangulation
154  {
155  //test, if it is the same point as the first point
156  if ( p.x() == mPointVector[0]->x() && p.y() == mPointVector[0]->y() )
157  {
158  //second point is the same as the first point
159  removeLastPoint();
160  return 0;
161  }
162 
163  const unsigned int edgeFromPoint0ToPoint1 /* 2 */ = insertEdge( -10, -10, 1, false, false );//edge pointing from point 0 to point 1
164  const unsigned int edgeFromPoint1ToPoint0 /* 3 */ = insertEdge( edgeFromPoint0ToPoint1, -10, 0, false, false ); //edge pointing from point 1 to point 0
165  const unsigned int edgeFromVirtualToPoint1Side1 /* 4 */ = insertEdge( -10, -10, 1, false, false ); //edge pointing from the virtual point to point 1
166  const unsigned int edgeFromPoint1ToVirtualSide1 /* 5 */ = insertEdge( edgeFromVirtualToPoint1Side1, 1, -1, false, false ); //edge pointing from point 1 to the virtual point
167  const unsigned int edgeFromVirtualToPoint1Side2 /* 6 */ = insertEdge( -10, edgeFromPoint1ToPoint0, 1, false, false );
168  const unsigned int edgeFromPoint1ToVirtualSide2 /* 7 */ = insertEdge( edgeFromVirtualToPoint1Side2, edgeFromVirtualToPoint1Side1, -1, false, false );
169  const unsigned int edgeFromVirtualToPoint0Side2 /* 8 */ = insertEdge( -10, -10, 0, false, false );
170  const unsigned int edgeFromPoint0ToVirtualSide2 /* 9 */ = insertEdge( edgeFromVirtualToPoint0Side2, edgeFromVirtualToPoint1Side2, -1, false, false );
171  mHalfEdge.at( edgeFromPoint1ToPoint0 )->setNext( edgeFromPoint0ToVirtualSide2 );
172  mHalfEdge.at( edgeFromPoint0ToPoint1 )->setDual( edgeFromPoint1ToPoint0 );
173  mHalfEdge.at( edgeFromPoint0ToPoint1 )->setNext( edgeFromPoint1ToVirtualSide1 );
174  mHalfEdge.at( edgeFromVirtualToPoint1Side1 )->setDual( edgeFromPoint1ToVirtualSide1 );
175  mHalfEdge.at( edgeFromVirtualToPoint1Side1 )->setNext( edgeFromPoint1ToVirtualSide2 );
176  mHalfEdge.at( 0 )->setNext( static_cast<int>( edgeFromVirtualToPoint0Side2 ) );
177  mHalfEdge.at( 1 )->setNext( static_cast<int>( edgeFromPoint0ToPoint1 ) );
178  mHalfEdge.at( edgeFromVirtualToPoint1Side2 )->setDual( edgeFromPoint1ToVirtualSide2 );
179  mHalfEdge.at( edgeFromVirtualToPoint0Side2 )->setDual( edgeFromPoint0ToVirtualSide2 );
180  mHalfEdge.at( edgeFromVirtualToPoint0Side2 )->setNext( 0 );
181  mEdgeInside = 3;
182  mEdgeOutside = edgeFromPoint0ToPoint1;
183  mDimension = 1;
184  }
185  else if ( mDimension == 1 )
186  {
187  if ( mEdgeOutside < 0 || mHalfEdge[mEdgeOutside]->getPoint() < 0 || mHalfEdge[mHalfEdge[mEdgeOutside]->getDual()]->getPoint() < 0 )
188  mEdgeOutside = firstEdgeOutSide();
189  if ( mEdgeOutside < 0 || mHalfEdge[mEdgeOutside]->getPoint() < 0 || mHalfEdge[mHalfEdge[mEdgeOutside]->getDual()]->getPoint() < 0 )
190  return -100;
191 
192  const double leftOfNumber = MathUtils::leftOf( p, mPointVector[mHalfEdge[mHalfEdge[mEdgeOutside]->getDual()]->getPoint()], mPointVector[mHalfEdge[mEdgeOutside]->getPoint()] );
193  if ( fabs( leftOfNumber ) <= leftOfTresh )
194  {
195  // new point colinear with existing points
196  mDimension = 1;
197  // First find the edge which has the point in or the closest edge if the new point is outside
198  int closestEdge = -1;
199  double distance = std::numeric_limits<double>::max();
200  const int firstEdge = mEdgeOutside;
201  do
202  {
203  const int point1 = mHalfEdge[mEdgeOutside]->getPoint();
204  const int point2 = mHalfEdge[mHalfEdge[mEdgeOutside]->getDual()]->getPoint();
205  const double distance1 = p.distance( *mPointVector[point1] );
206  if ( distance1 <= leftOfTresh ) // point1 == new point
207  {
208  removeLastPoint();
209  return point1;
210  }
211  const double distance2 = p.distance( *mPointVector[point2] );
212  if ( distance2 <= leftOfTresh ) // point2 == new point
213  {
214  removeLastPoint();
215  return point2;
216  }
217 
218  const double edgeLength = mPointVector[point1]->distance( *mPointVector[point2] );
219 
220  if ( distance1 < edgeLength && distance2 < edgeLength )
221  {
222  //new point include in mEdgeOutside
223  const int newPoint = mPointVector.count() - 1;
224 
225  //edges that do not change
226  const int edgeFromNewPointToPoint1 = mEdgeOutside;
227  const int edgeFromNewPointToPoint2 = mHalfEdge[mEdgeOutside]->getDual();
228  //edges to modify
229  const int edgeFromPoint1ToVirtualSide2 = mHalfEdge[edgeFromNewPointToPoint1]->getNext();
230  const int edgeFromVirtualToPoint1Side1 = mHalfEdge[mHalfEdge[edgeFromNewPointToPoint2]->getNext()]->getNext();
231  const int edgeFromPoint2ToVirtualSide1 = mHalfEdge[edgeFromNewPointToPoint2]->getNext();
232  const int edgeFromVirtualToPoint2Side2 = mHalfEdge[mHalfEdge[edgeFromNewPointToPoint1]->getNext()]->getNext();
233  //insert new edges
234  const int edgeFromVirtualToNewPointSide1 = insertEdge( -10, edgeFromNewPointToPoint2, newPoint, false, false );
235  const int edgeFromNewPointToVirtualSide1 = insertEdge( edgeFromVirtualToNewPointSide1, edgeFromVirtualToPoint1Side1, -1, false, false );
236  const int edgeFromVirtualToNewPointSide2 = insertEdge( -10, edgeFromNewPointToPoint1, newPoint, false, false );
237  const int edgeFromNewPointToVirtualSide2 = insertEdge( edgeFromVirtualToNewPointSide2, edgeFromVirtualToPoint2Side2, -1, false, false );
238  const int edgeFromPoint1ToNewPoint = insertEdge( edgeFromNewPointToPoint1, edgeFromNewPointToVirtualSide1, newPoint, false, false );
239  const int edgeFromPoint2ToNewPoint = insertEdge( edgeFromNewPointToPoint2, edgeFromNewPointToVirtualSide2, newPoint, false, false );
240  mHalfEdge.at( edgeFromVirtualToNewPointSide1 )->setDual( edgeFromNewPointToVirtualSide1 );
241  mHalfEdge.at( edgeFromVirtualToNewPointSide2 )->setDual( edgeFromNewPointToVirtualSide2 );
242  //modify existing edges
243  mHalfEdge.at( edgeFromPoint1ToVirtualSide2 )->setNext( edgeFromVirtualToNewPointSide2 );
244  mHalfEdge.at( edgeFromVirtualToPoint1Side1 )->setNext( edgeFromPoint1ToNewPoint );
245  mHalfEdge.at( edgeFromPoint2ToVirtualSide1 )->setNext( edgeFromVirtualToNewPointSide1 );
246  mHalfEdge.at( edgeFromVirtualToPoint2Side2 )->setNext( edgeFromPoint2ToNewPoint );
247  mHalfEdge.at( edgeFromNewPointToPoint1 )->setDual( edgeFromPoint1ToNewPoint );
248  mHalfEdge.at( edgeFromNewPointToPoint2 )->setDual( edgeFromPoint2ToNewPoint );
249  return newPoint;
250  }
251  else
252  {
253  if ( distance1 < distance )
254  {
255  closestEdge = mEdgeOutside;
256  distance = distance1;
257  }
258  else if ( distance2 < distance )
259  {
260  closestEdge = mHalfEdge[mEdgeOutside]->getDual();
261  distance = distance2;
262  }
263  }
264  mEdgeOutside = mHalfEdge[mHalfEdge[mHalfEdge[mEdgeOutside]->getNext()]->getDual()]->getNext();
265  }
266  while ( mEdgeOutside != firstEdge && mHalfEdge[mEdgeOutside]->getPoint() != -1 && mHalfEdge[mHalfEdge[mEdgeOutside]->getDual()]->getPoint() != -1 );
267 
268  if ( closestEdge < 0 )
269  return -100; //something gets wrong
270 
271  //add the new colinear point linking it to the extremity of closest edge
272  const int extremPoint = mHalfEdge[closestEdge]->getPoint();
273  const int newPoint = mPointVector.count() - 1;
274  //edges that do not change
275  const int edgeFromExtremeToOpposite = mHalfEdge[closestEdge]->getDual();
276  //edges to modify
277  const int edgeFromVirtualToExtremeSide1 = mHalfEdge[mHalfEdge[closestEdge]->getNext()]->getDual();
278  const int edgeFromVirtualToExtremeSide2 = mHalfEdge[mHalfEdge[mHalfEdge[closestEdge]->getDual()]->getNext()]->getNext();
279  const int edgeFromExtremeToVirtualSide2 = mHalfEdge[edgeFromVirtualToExtremeSide2]->getDual();
280  //insert new edge
281  const int edgeFromExtremeToNewPoint = insertEdge( -10, -10, newPoint, false, false );
282  const int edgeFromNewPointToExtrem = insertEdge( edgeFromExtremeToNewPoint, edgeFromExtremeToVirtualSide2, extremPoint, false, false );
283  const int edgeFromNewPointToVirtualSide1 = insertEdge( -10, edgeFromVirtualToExtremeSide1, -1, false, false );
284  const int edgeFromVirtualToNewPointSide1 = insertEdge( edgeFromNewPointToVirtualSide1, -10, newPoint, false, false );
285  const int edgeFromNewPointToVirtualSide2 = insertEdge( -10, edgeFromVirtualToNewPointSide1, -1, false, false );
286  const int edgeFromVirtualToNewPointSide2 = insertEdge( edgeFromNewPointToVirtualSide2, edgeFromNewPointToExtrem, newPoint, false, false );
287  mHalfEdge.at( edgeFromExtremeToNewPoint )->setDual( edgeFromNewPointToExtrem );
288  mHalfEdge.at( edgeFromExtremeToNewPoint )->setNext( edgeFromNewPointToVirtualSide1 );
289  mHalfEdge.at( edgeFromNewPointToVirtualSide1 )->setDual( edgeFromVirtualToNewPointSide1 );
290  mHalfEdge.at( edgeFromNewPointToVirtualSide2 )->setDual( edgeFromVirtualToNewPointSide2 );
291  mHalfEdge.at( edgeFromVirtualToNewPointSide1 )->setNext( edgeFromNewPointToVirtualSide2 );
292  //modify existing edges
293  mHalfEdge.at( edgeFromVirtualToExtremeSide1 )->setNext( edgeFromExtremeToNewPoint );
294  mHalfEdge.at( edgeFromVirtualToExtremeSide2 )->setNext( edgeFromExtremeToOpposite );
295  mHalfEdge.at( edgeFromExtremeToVirtualSide2 )->setNext( edgeFromVirtualToNewPointSide2 );
296 
297  return newPoint;
298  }
299  else if ( leftOfNumber >= leftOfTresh )
300  {
301  // new point on the right of mEdgeOutside
302  mEdgeOutside = mHalfEdge[mEdgeOutside]->getDual();
303  }
304  mDimension = 2;
305  const int newPoint = mPointVector.count() - 1;
306  //buil the 2D dimension triangulation
307  //First clock wise
308  int cwEdge = mEdgeOutside;
309  while ( mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[cwEdge]->getNext()]->getDual()]->getNext()]->getPoint() != -1 )
310  {
311  mHalfEdge[mHalfEdge[ cwEdge ]->getNext()]->setPoint( newPoint );
312  cwEdge = mHalfEdge[mHalfEdge[mHalfEdge[cwEdge]->getNext()]->getDual()]->getNext();
313  }
314 
315  const int edgeFromLastCwPointToVirtualPoint = mHalfEdge[mHalfEdge[mHalfEdge[cwEdge]->getNext()]->getDual()]->getNext();
316  const int edgeFromLastCwPointToNewPointPoint = mHalfEdge[ cwEdge ]->getNext();
317  const int edgeFromNewPointPointToLastCwPoint = mHalfEdge[ edgeFromLastCwPointToNewPointPoint ]->getDual();
318  //connect the new point
319  const int edgeFromNewPointtoVirtualPoint = insertEdge( -10, -10, -1, false, false );
320  const int edgeFromVirtualPointToNewPoint = insertEdge( edgeFromNewPointtoVirtualPoint, edgeFromNewPointPointToLastCwPoint, newPoint, false, false );
321  mHalfEdge.at( edgeFromLastCwPointToNewPointPoint )->setPoint( newPoint );
322  mHalfEdge.at( edgeFromNewPointtoVirtualPoint )->setDual( edgeFromVirtualPointToNewPoint );
323  mHalfEdge.at( edgeFromLastCwPointToVirtualPoint )->setNext( edgeFromVirtualPointToNewPoint );
324 
325  //First counter clock wise
326  int ccwEdge = mEdgeOutside;
327  while ( mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[ccwEdge]->getDual()]->getNext()]->getDual()]->getNext()]->getPoint() != -1 )
328  {
329  mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[ ccwEdge ]->getNext()]->getNext()]->getDual()]->setPoint( newPoint );
330  ccwEdge = mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[ccwEdge]->getDual()]->getNext()]->getDual()]->getNext()]->getDual();
331  }
332 
333  const int edgeToLastCcwPoint = mHalfEdge[mHalfEdge[mHalfEdge[ccwEdge]->getDual()]->getNext()]->getDual();
334  const int edgeFromLastCcwPointToNewPoint = mHalfEdge[edgeToLastCcwPoint]->getNext();
335  mHalfEdge.at( edgeFromNewPointtoVirtualPoint )->setNext( edgeToLastCcwPoint );
336  mHalfEdge.at( edgeFromLastCcwPointToNewPoint )->setNext( edgeFromNewPointtoVirtualPoint );
337  mHalfEdge.at( edgeFromLastCcwPointToNewPoint )->setPoint( newPoint );
338  }
339  else
340  {
341  const int number = baseEdgeOfTriangle( p );
342 
343  //point is outside the convex hull----------------------------------------------------
344  if ( number == -10 )
345  {
346  unsigned int cwEdge = mEdgeOutside;//the last visible edge clockwise from mEdgeOutside
347  unsigned int ccwEdge = mEdgeOutside;//the last visible edge counterclockwise from mEdgeOutside
348 
349  //mEdgeOutside is in each case visible
350  mHalfEdge[mHalfEdge[mEdgeOutside]->getNext()]->setPoint( mPointVector.count() - 1 );
351 
352  //find cwEdge and replace the virtual point with the new point when necessary (equivalent to while the hull is not convex going clock wise)
353  while ( MathUtils::leftOf( *mPointVector[ mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[cwEdge]->getNext()]->getDual()]->getNext()]->getPoint()],
354  &p, mPointVector[ mHalfEdge[cwEdge]->getPoint()] ) < ( -leftOfTresh ) )
355  {
356  //set the point number of the necessary edge to the actual point instead of the virtual point
357  mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[cwEdge]->getNext()]->getDual()]->getNext()]->getNext()]->setPoint( mPointVector.count() - 1 );
358  //advance cwedge one edge further clockwise
359  cwEdge = ( unsigned int )mHalfEdge[mHalfEdge[mHalfEdge[cwEdge]->getNext()]->getDual()]->getNext();
360  }
361 
362  //build the necessary connections with the virtual point
363  const unsigned int edge1 = insertEdge( mHalfEdge[cwEdge]->getNext(), -10, mHalfEdge[cwEdge]->getPoint(), false, false );//edge pointing from the new point to the last visible point clockwise
364  const unsigned int edge2 = insertEdge( mHalfEdge[mHalfEdge[cwEdge]->getNext()]->getDual(), -10, -1, false, false );//edge pointing from the last visible point to the virtual point
365  const unsigned int edge3 = insertEdge( -10, edge1, mPointVector.count() - 1, false, false );//edge pointing from the virtual point to new point
366 
367  //adjust the other pointers
368  mHalfEdge[mHalfEdge[mHalfEdge[cwEdge]->getNext()]->getDual()]->setDual( edge2 );
369  mHalfEdge[mHalfEdge[cwEdge]->getNext()]->setDual( edge1 );
370  mHalfEdge[edge1]->setNext( edge2 );
371  mHalfEdge[edge2]->setNext( edge3 );
372 
373  //find ccwedge and replace the virtual point with the new point when necessary
374  while ( MathUtils::leftOf( *mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[ccwEdge]->getNext()]->getNext()]->getPoint()], mPointVector[mPointVector.count() - 1], mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[ccwEdge]->getNext()]->getNext()]->getDual()]->getNext()]->getPoint()] ) < ( -leftOfTresh ) )
375  {
376  //set the point number of the necessary edge to the actual point instead of the virtual point
377  mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[ccwEdge]->getNext()]->getNext()]->getDual()]->setPoint( mPointVector.count() - 1 );
378  //advance ccwedge one edge further counterclockwise
379  ccwEdge = mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[ccwEdge]->getNext()]->getNext()]->getDual()]->getNext()]->getNext();
380  }
381 
382  //build the necessary connections with the virtual point
383  const unsigned int edge4 = insertEdge( mHalfEdge[mHalfEdge[ccwEdge]->getNext()]->getNext(), -10, mPointVector.count() - 1, false, false );//points from the last visible point counterclockwise to the new point
384  const unsigned int edge5 = insertEdge( edge3, -10, -1, false, false );//points from the new point to the virtual point
385  const unsigned int edge6 = insertEdge( mHalfEdge[mHalfEdge[mHalfEdge[ccwEdge]->getNext()]->getNext()]->getDual(), edge4, mHalfEdge[mHalfEdge[ccwEdge]->getDual()]->getPoint(), false, false );//points from the virtual point to the last visible point counterclockwise
386 
387 
388 
389  //adjust the other pointers
390  mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[ccwEdge]->getNext()]->getNext()]->getDual()]->setDual( edge6 );
391  mHalfEdge[mHalfEdge[mHalfEdge[ccwEdge]->getNext()]->getNext()]->setDual( edge4 );
392  mHalfEdge[edge4]->setNext( edge5 );
393  mHalfEdge[edge5]->setNext( edge6 );
394  mHalfEdge[edge3]->setDual( edge5 );
395 
396  //now test the HalfEdge at the former convex hull for swappint
397  unsigned int index = ccwEdge;
398  unsigned int toswap;
399  while ( true )
400  {
401  toswap = index;
402  index = mHalfEdge[mHalfEdge[mHalfEdge[index]->getNext()]->getDual()]->getNext();
403  checkSwapRecursively( toswap, 0 );
404  if ( toswap == cwEdge )
405  {
406  break;
407  }
408  }
409  }
410  else if ( number >= 0 ) //point is inside the convex hull------------------------------------------------------
411  {
412  const int nextnumber = mHalfEdge[number]->getNext();
413  const int nextnextnumber = mHalfEdge[mHalfEdge[number]->getNext()]->getNext();
414 
415  //insert 6 new HalfEdges for the connections to the vertices of the triangle
416  const unsigned int edge1 = insertEdge( -10, nextnumber, mHalfEdge[number]->getPoint(), false, false );
417  const unsigned int edge2 = insertEdge( static_cast<int>( edge1 ), -10, mPointVector.count() - 1, false, false );
418  const unsigned int edge3 = insertEdge( -10, nextnextnumber, mHalfEdge[nextnumber]->getPoint(), false, false );
419  const unsigned int edge4 = insertEdge( static_cast<int>( edge3 ), static_cast<int>( edge1 ), mPointVector.count() - 1, false, false );
420  const unsigned int edge5 = insertEdge( -10, number, mHalfEdge[nextnextnumber]->getPoint(), false, false );
421  const unsigned int edge6 = insertEdge( static_cast<int>( edge5 ), static_cast<int>( edge3 ), mPointVector.count() - 1, false, false );
422 
423 
424  mHalfEdge.at( edge1 )->setDual( static_cast<int>( edge2 ) );
425  mHalfEdge.at( edge2 )->setNext( static_cast<int>( edge5 ) );
426  mHalfEdge.at( edge3 )->setDual( static_cast<int>( edge4 ) );
427  mHalfEdge.at( edge5 )->setDual( static_cast<int>( edge6 ) );
428  mHalfEdge.at( number )->setNext( static_cast<int>( edge2 ) );
429  mHalfEdge.at( nextnumber )->setNext( static_cast<int>( edge4 ) );
430  mHalfEdge.at( nextnextnumber )->setNext( static_cast<int>( edge6 ) );
431 
432  //check, if there are swaps necessary
433  checkSwapRecursively( number, 0 );
434  checkSwapRecursively( nextnumber, 0 );
435  checkSwapRecursively( nextnextnumber, 0 );
436  }
437  //the point is exactly on an existing edge (the number of the edge is stored in the variable 'mEdgeWithPoint'---------------
438  else if ( number == -20 )
439  {
440  //point exactly on edge;
441 
442  //check if new point is the same than one extremity
443  const int point1 = mHalfEdge[mEdgeWithPoint]->getPoint();
444  const int point2 = mHalfEdge[mHalfEdge[mEdgeWithPoint]->getDual()]->getPoint();
445  const double distance1 = p.distance( *mPointVector[point1] );
446  if ( distance1 <= leftOfTresh ) // point1 == new point
447  {
448  removeLastPoint();
449  return point1;
450  }
451  const double distance2 = p.distance( *mPointVector[point2] );
452  if ( distance2 <= leftOfTresh ) // point2 == new point
453  {
454  removeLastPoint();
455  return point2;
456  }
457 
458  const int edgea = mEdgeWithPoint;
459  const int edgeb = mHalfEdge[mEdgeWithPoint]->getDual();
460  const int edgec = mHalfEdge[edgea]->getNext();
461  const int edged = mHalfEdge[edgec]->getNext();
462  const int edgee = mHalfEdge[edgeb]->getNext();
463  const int edgef = mHalfEdge[edgee]->getNext();
464 
465  //insert the six new edges
466  const int nedge1 = insertEdge( -10, mHalfEdge[edgea]->getNext(), mHalfEdge[edgea]->getPoint(), false, false );
467  const int nedge2 = insertEdge( nedge1, -10, mPointVector.count() - 1, false, false );
468  const int nedge3 = insertEdge( -10, edged, mHalfEdge[edgec]->getPoint(), false, false );
469  const int nedge4 = insertEdge( nedge3, nedge1, mPointVector.count() - 1, false, false );
470  const int nedge5 = insertEdge( -10, edgef, mHalfEdge[edgee]->getPoint(), false, false );
471  const int nedge6 = insertEdge( nedge5, edgeb, mPointVector.count() - 1, false, false );
472 
473  //adjust the triangular structure
474  mHalfEdge[nedge1]->setDual( nedge2 );
475  mHalfEdge[nedge2]->setNext( nedge5 );
476  mHalfEdge[nedge3]->setDual( nedge4 );
477  mHalfEdge[nedge5]->setDual( nedge6 );
478  mHalfEdge[edgea]->setPoint( mPointVector.count() - 1 );
479  mHalfEdge[edgea]->setNext( nedge3 );
480  mHalfEdge[edgec]->setNext( nedge4 );
481  mHalfEdge[edgee]->setNext( nedge6 );
482  mHalfEdge[edgef]->setNext( nedge2 );
483 
484  //swap edges if necessary
485  checkSwapRecursively( edgec, 0 );
486  checkSwapRecursively( edged, 0 );
487  checkSwapRecursively( edgee, 0 );
488  checkSwapRecursively( edgef, 0 );
489  }
490 
491  else if ( number == -100 || number == -5 )//this means unknown problems or a numerical error occurred in 'baseEdgeOfTriangle'
492  {
493  //QgsDebugMsg( "point has not been inserted because of unknown problems" );
494  removeLastPoint();
495  return -100;
496  }
497  else if ( number == -25 )//this means that the point has already been inserted in the triangulation
498  {
499  //Take the higher z-Value in case of two equal points
500  QgsPoint *newPoint = mPointVector[mPointVector.count() - 1];
501  QgsPoint *existingPoint = mPointVector[mTwiceInsPoint];
502  existingPoint->setZ( std::max( newPoint->z(), existingPoint->z() ) );
503 
504  removeLastPoint();
505  return mTwiceInsPoint;
506  }
507  }
508 
509  return ( mPointVector.count() - 1 );
510 }
511 
512 int QgsDualEdgeTriangulation::baseEdgeOfPoint( int point )
513 {
514  unsigned int actedge = mEdgeInside;//starting edge
515 
516  if ( mPointVector.count() < 4 || point == -1 || mDimension == 1 ) //at the beginning, mEdgeInside is not defined yet
517  {
518  int fromVirtualPoint = -1;
519  //first find pointingedge(an edge pointing to p1, priority to edge that no come from virtual point)
520  for ( int i = 0; i < mHalfEdge.count(); i++ )
521  {
522  if ( mHalfEdge[i]->getPoint() == point )//we found one
523  {
524  if ( mHalfEdge[mHalfEdge[i]->getDual()]->getPoint() != -1 )
525  return i;
526  else
527  fromVirtualPoint = i;
528  }
529  }
530  return fromVirtualPoint;
531  }
532 
533  int control = 0;
534 
535  while ( true )//otherwise, start the search
536  {
537  control += 1;
538  if ( control > 1000000 )
539  {
540  //QgsDebugMsg( QStringLiteral( "warning, endless loop" ) );
541 
542  //use the secure and slow method
543  //qWarning( "******************warning, using the slow method in baseEdgeOfPoint****************************************" );
544  for ( int i = 0; i < mHalfEdge.count(); i++ )
545  {
546  if ( mHalfEdge[i]->getPoint() == point && mHalfEdge[mHalfEdge[i]->getNext()]->getPoint() != -1 )//we found it
547  {
548  return i;
549  }
550  }
551  }
552 
553  const int fromPoint = mHalfEdge[mHalfEdge[actedge]->getDual()]->getPoint();
554  const int toPoint = mHalfEdge[actedge]->getPoint();
555 
556  if ( fromPoint == -1 || toPoint == -1 )//this would cause a crash. Therefore we use the slow method in this case
557  {
558  for ( int i = 0; i < mHalfEdge.count(); i++ )
559  {
560  if ( mHalfEdge[i]->getPoint() == point && mHalfEdge[mHalfEdge[i]->getNext()]->getPoint() != -1 )//we found it
561  {
562  mEdgeInside = i;
563  return i;
564  }
565  }
566  }
567 
568  const double leftOfNumber = MathUtils::leftOf( *mPointVector[point], mPointVector[mHalfEdge[mHalfEdge[actedge]->getDual()]->getPoint()], mPointVector[mHalfEdge[actedge]->getPoint()] );
569 
570 
571  if ( mHalfEdge[actedge]->getPoint() == point && mHalfEdge[mHalfEdge[actedge]->getNext()]->getPoint() != -1 )//we found the edge
572  {
573  mEdgeInside = actedge;
574  return actedge;
575  }
576  else if ( leftOfNumber <= 0.0 )
577  {
578  actedge = mHalfEdge[actedge]->getNext();
579  }
580  else
581  {
582  actedge = mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[actedge]->getDual()]->getNext()]->getNext()]->getDual();
583  }
584  }
585 }
586 
587 int QgsDualEdgeTriangulation::baseEdgeOfTriangle( const QgsPoint &point )
588 {
589  unsigned int actEdge = mEdgeInside;//start with an edge which does not point to the virtual point
590  if ( mHalfEdge.at( actEdge )->getPoint() < 0 )
591  actEdge = mHalfEdge.at( mHalfEdge.at( mHalfEdge.at( actEdge )->getDual() )->getNext() )->getDual(); //get an real inside edge
592  if ( mHalfEdge.at( mHalfEdge.at( actEdge )->getDual() )->getPoint() < 0 )
593  actEdge = mHalfEdge.at( mHalfEdge.at( actEdge )->getNext() )->getDual();
594 
595  int counter = 0;//number of consecutive successful left-of-tests
596  int nulls = 0;//number of left-of-tests, which returned 0. 1 means, that the point is on a line, 2 means that it is on an existing point
597  int numInstabs = 0;//number of suspect left-of-tests due to 'leftOfTresh'
598  int runs = 0;//counter for the number of iterations in the loop to prevent an endless loop
599  int firstEndPoint = 0, secEndPoint = 0, thEndPoint = 0, fouEndPoint = 0; //four numbers of endpoints in cases when two left-of-test are 0
600 
601  while ( true )
602  {
603  if ( runs > MAX_BASE_ITERATIONS )//prevents endless loops
604  {
605  //QgsDebugMsg( "warning, probable endless loop detected" );
606  return -100;
607  }
608 
609  const double leftOfValue = MathUtils::leftOf( point, mPointVector.at( mHalfEdge.at( mHalfEdge.at( actEdge )->getDual() )->getPoint() ), mPointVector.at( mHalfEdge.at( actEdge )->getPoint() ) );
610 
611  if ( leftOfValue < ( -leftOfTresh ) )//point is on the left side
612  {
613  counter += 1;
614  if ( counter == 3 )//three successful passes means that we have found the triangle
615  {
616  break;
617  }
618  }
619  else if ( fabs( leftOfValue ) <= leftOfTresh ) //point is exactly in the line of the edge
620  {
621  if ( nulls == 0 )
622  {
623  //store the numbers of the two endpoints of the line
624  firstEndPoint = mHalfEdge.at( mHalfEdge.at( actEdge )->getDual() )->getPoint();
625  secEndPoint = mHalfEdge.at( actEdge )->getPoint();
626  }
627  else if ( nulls == 1 )
628  {
629  //store the numbers of the two endpoints of the line
630  thEndPoint = mHalfEdge.at( mHalfEdge.at( actEdge )->getDual() )->getPoint();
631  fouEndPoint = mHalfEdge.at( actEdge )->getPoint();
632  }
633  counter += 1;
634  mEdgeWithPoint = actEdge;
635  nulls += 1;
636  if ( counter == 3 )//three successful passes means that we have found the triangle
637  {
638  break;
639  }
640  }
641  else//point is on the right side
642  {
643  actEdge = mHalfEdge.at( actEdge )->getDual();
644  counter = 1;
645  nulls = 0;
646  numInstabs = 0;
647  }
648  actEdge = mHalfEdge.at( actEdge )->getNext();
649  if ( mHalfEdge.at( actEdge )->getPoint() == -1 )//the half edge points to the virtual point
650  {
651  if ( nulls == 1 )//point is exactly on the convex hull
652  {
653  return -20;
654  }
655  mEdgeOutside = ( unsigned int )mHalfEdge.at( mHalfEdge.at( actEdge )->getNext() )->getNext();
656  mEdgeInside = mHalfEdge.at( mHalfEdge.at( mEdgeOutside )->getDual() )->getNext();
657  return -10;//the point is outside the convex hull
658  }
659  runs++;
660  }
661 
662  if ( numInstabs > 0 )//we hit an existing point or a numerical instability occurred
663  {
664  // QgsDebugMsg("numerical instability occurred");
665  mUnstableEdge = actEdge;
666  return -5;
667  }
668 
669  if ( nulls == 2 )
670  {
671  //find out the number of the point, which has already been inserted
672  if ( firstEndPoint == thEndPoint || firstEndPoint == fouEndPoint )
673  {
674  //firstendp is the number of the point which has been inserted twice
675  mTwiceInsPoint = firstEndPoint;
676  // QgsDebugMsg(QString("point nr %1 already inserted").arg(firstendp));
677  }
678  else if ( secEndPoint == thEndPoint || secEndPoint == fouEndPoint )
679  {
680  //secendp is the number of the point which has been inserted twice
681  mTwiceInsPoint = secEndPoint;
682  // QgsDebugMsg(QString("point nr %1 already inserted").arg(secendp));
683  }
684 
685  return -25;//return the code for a point that is already contained in the triangulation
686  }
687 
688  if ( nulls == 1 )//point is on an existing edge
689  {
690  return -20;
691  }
692 
693  mEdgeInside = actEdge;
694 
695  int nr1, nr2, nr3;
696  nr1 = mHalfEdge.at( actEdge )->getPoint();
697  nr2 = mHalfEdge.at( mHalfEdge.at( actEdge )->getNext() )->getPoint();
698  nr3 = mHalfEdge.at( mHalfEdge.at( mHalfEdge.at( actEdge )->getNext() )->getNext() )->getPoint();
699  const double x1 = mPointVector.at( nr1 )->x();
700  const double y1 = mPointVector.at( nr1 )->y();
701  const double x2 = mPointVector.at( nr2 )->x();
702  const double y2 = mPointVector.at( nr2 )->y();
703  const double x3 = mPointVector.at( nr3 )->x();
704  const double y3 = mPointVector.at( nr3 )->y();
705 
706  //now make sure that always the same edge is returned
707  if ( x1 < x2 && x1 < x3 )//return the edge which points to the point with the lowest x-coordinate
708  {
709  return actEdge;
710  }
711  else if ( x2 < x1 && x2 < x3 )
712  {
713  return mHalfEdge.at( actEdge )->getNext();
714  }
715  else if ( x3 < x1 && x3 < x2 )
716  {
717  return mHalfEdge.at( mHalfEdge.at( actEdge )->getNext() )->getNext();
718  }
719  //in case two x-coordinates are the same, the edge pointing to the point with the lower y-coordinate is returned
720  else if ( x1 == x2 )
721  {
722  if ( y1 < y2 )
723  {
724  return actEdge;
725  }
726  else if ( y2 < y1 )
727  {
728  return mHalfEdge.at( actEdge )->getNext();
729  }
730  }
731  else if ( x2 == x3 )
732  {
733  if ( y2 < y3 )
734  {
735  return mHalfEdge.at( actEdge )->getNext();
736  }
737  else if ( y3 < y2 )
738  {
739  return mHalfEdge.at( mHalfEdge.at( actEdge )->getNext() )->getNext();
740  }
741  }
742  else if ( x1 == x3 )
743  {
744  if ( y1 < y3 )
745  {
746  return actEdge;
747  }
748  else if ( y3 < y1 )
749  {
750  return mHalfEdge.at( mHalfEdge.at( actEdge )->getNext() )->getNext();
751  }
752  }
753  return -100;//this means a bug happened
754 }
755 
756 bool QgsDualEdgeTriangulation::calcNormal( double x, double y, QgsPoint &result )
757 {
758  if ( mTriangleInterpolator )
759  {
760  return mTriangleInterpolator->calcNormVec( x, y, result );
761  }
762  else
763  {
764  QgsDebugMsg( QStringLiteral( "warning, null pointer" ) );
765  return false;
766  }
767 }
768 
769 bool QgsDualEdgeTriangulation::calcPoint( double x, double y, QgsPoint &result )
770 {
771  if ( mTriangleInterpolator )
772  {
773  return mTriangleInterpolator->calcPoint( x, y, result );
774  }
775  else
776  {
777  QgsDebugMsg( QStringLiteral( "warning, null pointer" ) );
778  return false;
779  }
780 }
781 
782 bool QgsDualEdgeTriangulation::checkSwapRecursively( unsigned int edge, unsigned int recursiveDeep )
783 {
784  if ( swapPossible( edge ) )
785  {
786  QgsPoint *pta = mPointVector.at( mHalfEdge.at( edge )->getPoint() );
787  QgsPoint *ptb = mPointVector.at( mHalfEdge.at( mHalfEdge.at( edge )->getNext() )->getPoint() );
788  QgsPoint *ptc = mPointVector.at( mHalfEdge.at( mHalfEdge.at( mHalfEdge.at( edge )->getNext( ) )->getNext() )->getPoint() );
789  QgsPoint *ptd = mPointVector.at( mHalfEdge.at( mHalfEdge.at( mHalfEdge.at( edge )->getDual() )->getNext() )->getPoint() );
790  if ( inCircle( *ptd, *pta, *ptb, *ptc ) /*&& recursiveDeep < 100*/ ) //empty circle criterion violated
791  {
792  doSwapRecursively( edge, recursiveDeep );//swap the edge (recursive)
793  return true;
794  }
795  }
796  return false;
797 }
798 
799 bool QgsDualEdgeTriangulation::isEdgeNeedSwap( unsigned int edge ) const
800 {
801  if ( swapPossible( edge ) )
802  {
803  QgsPoint *pta = mPointVector[mHalfEdge[edge]->getPoint()];
804  QgsPoint *ptb = mPointVector[mHalfEdge[mHalfEdge[edge]->getNext()]->getPoint()];
805  QgsPoint *ptc = mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[edge]->getNext()]->getNext()]->getPoint()];
806  QgsPoint *ptd = mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[edge]->getDual()]->getNext()]->getPoint()];
807  return inCircle( *ptd, *pta, *ptb, *ptc );
808  }
809 
810  return false;
811 
812 }
813 
814 void QgsDualEdgeTriangulation::doOnlySwap( unsigned int edge )
815 {
816  const unsigned int edge1 = edge;
817  const unsigned int edge2 = mHalfEdge[edge]->getDual();
818  const unsigned int edge3 = mHalfEdge[edge]->getNext();
819  const unsigned int edge4 = mHalfEdge[mHalfEdge[edge]->getNext()]->getNext();
820  const unsigned int edge5 = mHalfEdge[mHalfEdge[edge]->getDual()]->getNext();
821  const unsigned int edge6 = mHalfEdge[mHalfEdge[mHalfEdge[edge]->getDual()]->getNext()]->getNext();
822  mHalfEdge[edge1]->setNext( edge4 );//set the necessary nexts
823  mHalfEdge[edge2]->setNext( edge6 );
824  mHalfEdge[edge3]->setNext( edge2 );
825  mHalfEdge[edge4]->setNext( edge5 );
826  mHalfEdge[edge5]->setNext( edge1 );
827  mHalfEdge[edge6]->setNext( edge3 );
828  mHalfEdge[edge1]->setPoint( mHalfEdge[edge3]->getPoint() );//change the points to which edge1 and edge2 point
829  mHalfEdge[edge2]->setPoint( mHalfEdge[edge5]->getPoint() );
830 }
831 
832 void QgsDualEdgeTriangulation::doSwapRecursively( unsigned int edge, unsigned int recursiveDeep )
833 {
834  const unsigned int edge1 = edge;
835  const unsigned int edge2 = mHalfEdge.at( edge )->getDual();
836  const unsigned int edge3 = mHalfEdge.at( edge )->getNext();
837  const unsigned int edge4 = mHalfEdge.at( mHalfEdge.at( edge )->getNext() )->getNext();
838  const unsigned int edge5 = mHalfEdge.at( mHalfEdge.at( edge )->getDual() )->getNext();
839  const unsigned int edge6 = mHalfEdge.at( mHalfEdge.at( mHalfEdge.at( edge )->getDual() )->getNext() )->getNext();
840  mHalfEdge.at( edge1 )->setNext( edge4 );//set the necessary nexts
841  mHalfEdge.at( edge2 )->setNext( edge6 );
842  mHalfEdge.at( edge3 )->setNext( edge2 );
843  mHalfEdge.at( edge4 )->setNext( edge5 );
844  mHalfEdge.at( edge5 )->setNext( edge1 );
845  mHalfEdge.at( edge6 )->setNext( edge3 );
846  mHalfEdge.at( edge1 )->setPoint( mHalfEdge.at( edge3 )->getPoint() );//change the points to which edge1 and edge2 point
847  mHalfEdge.at( edge2 )->setPoint( mHalfEdge.at( edge5 )->getPoint() );
848  recursiveDeep++;
849 
850  if ( recursiveDeep < 100 )
851  {
852  checkSwapRecursively( edge3, recursiveDeep );
853  checkSwapRecursively( edge6, recursiveDeep );
854  checkSwapRecursively( edge4, recursiveDeep );
855  checkSwapRecursively( edge5, recursiveDeep );
856  }
857  else
858  {
859  QStack<int> edgesToSwap;
860  edgesToSwap.push( edge3 );
861  edgesToSwap.push( edge6 );
862  edgesToSwap.push( edge4 );
863  edgesToSwap.push( edge5 );
864  int loopCount = 0;
865  while ( !edgesToSwap.isEmpty() && loopCount < 10000 )
866  {
867  loopCount++;
868  const unsigned int e1 = edgesToSwap.pop();
869  if ( isEdgeNeedSwap( e1 ) )
870  {
871  const unsigned int e2 = mHalfEdge.at( e1 )->getDual();
872  const unsigned int e3 = mHalfEdge.at( e1 )->getNext();
873  const unsigned int e4 = mHalfEdge.at( mHalfEdge.at( e1 )->getNext() )->getNext();
874  const unsigned int e5 = mHalfEdge.at( mHalfEdge.at( e1 )->getDual() )->getNext();
875  const unsigned int e6 = mHalfEdge.at( mHalfEdge.at( mHalfEdge.at( e1 )->getDual() )->getNext() )->getNext();
876  mHalfEdge.at( e1 )->setNext( e4 );//set the necessary nexts
877  mHalfEdge.at( e2 )->setNext( e6 );
878  mHalfEdge.at( e3 )->setNext( e2 );
879  mHalfEdge.at( e4 )->setNext( e5 );
880  mHalfEdge.at( e5 )->setNext( e1 );
881  mHalfEdge.at( e6 )->setNext( e3 );
882  mHalfEdge.at( e1 )->setPoint( mHalfEdge.at( e3 )->getPoint() );//change the points to which edge1 and edge2 point
883  mHalfEdge.at( e2 )->setPoint( mHalfEdge.at( e5 )->getPoint() );
884 
885  edgesToSwap.push( e3 );
886  edgesToSwap.push( e6 );
887  edgesToSwap.push( e4 );
888  edgesToSwap.push( e5 );
889  }
890  }
891  }
892 
893 }
894 
895 #if 0
896 void DualEdgeTriangulation::draw( QPainter *p, double xlowleft, double ylowleft, double xupright, double yupright, double width, double height ) const
897 {
898  //if mPointVector is empty, there is nothing to do
899  if ( mPointVector.isEmpty() )
900  {
901  return;
902  }
903 
904  p->setPen( mEdgeColor );
905 
906  bool *control = new bool[mHalfEdge.count()];//controllarray that no edge is painted twice
907  bool *control2 = new bool[mHalfEdge.count()];//controllarray for the flat triangles
908 
909  for ( unsigned int i = 0; i <= mHalfEdge.count() - 1; i++ )
910  {
911  control[i] = false;
912  control2[i] = false;
913  }
914 
915  if ( ( ( xupright - xlowleft ) / width ) > ( ( yupright - ylowleft ) / height ) )
916  {
917  double lowerborder = -( height * ( xupright - xlowleft ) / width - yupright );//real world coordinates of the lower widget border. This is useful to know because of the HalfEdge bounding box test
918  for ( unsigned int i = 0; i < mHalfEdge.count() - 1; i++ )
919  {
920  if ( mHalfEdge[i]->getPoint() == -1 || mHalfEdge[mHalfEdge[i]->getDual()]->getPoint() == -1 )
921  {continue;}
922 
923  //check, if the edge belongs to a flat triangle, remove this later
924  if ( !control2[i] )
925  {
926  double p1, p2, p3;
927  if ( mHalfEdge[i]->getPoint() != -1 && mHalfEdge[mHalfEdge[i]->getNext()]->getPoint() != -1 && mHalfEdge[mHalfEdge[mHalfEdge[i]->getNext()]->getNext()]->getPoint() != -1 )
928  {
929  p1 = mPointVector[mHalfEdge[i]->getPoint()]->getZ();
930  p2 = mPointVector[mHalfEdge[mHalfEdge[i]->getNext()]->getPoint()]->getZ();
931  p3 = mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[i]->getNext()]->getNext()]->getPoint()]->getZ();
932  if ( p1 == p2 && p2 == p3 && halfEdgeBBoxTest( i, xlowleft, lowerborder, xupright, yupright ) && halfEdgeBBoxTest( mHalfEdge[i]->getNext(), xlowleft, lowerborder, xupright, yupright ) && halfEdgeBBoxTest( mHalfEdge[mHalfEdge[i]->getNext()]->getNext(), xlowleft, lowerborder, xupright, yupright ) )//draw the triangle
933  {
934  QPointArray pa( 3 );
935  pa.setPoint( 0, ( mPointVector[mHalfEdge[i]->getPoint()]->getX() - xlowleft ) / ( xupright - xlowleft )*width, ( yupright - mPointVector[mHalfEdge[i]->getPoint()]->getY() ) / ( xupright - xlowleft )*width );
936  pa.setPoint( 1, ( mPointVector[mHalfEdge[mHalfEdge[i]->getNext()]->getPoint()]->getX() - xlowleft ) / ( xupright - xlowleft )*width, ( yupright - mPointVector[mHalfEdge[mHalfEdge[i]->getNext()]->getPoint()]->getY() ) / ( xupright - xlowleft )*width );
937  pa.setPoint( 2, ( mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[i]->getNext()]->getNext()]->getPoint()]->getX() - xlowleft ) / ( xupright - xlowleft )*width, ( yupright - mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[i]->getNext()]->getNext()]->getPoint()]->getY() ) / ( xupright - xlowleft )*width );
938  QColor c( 255, 0, 0 );
939  p->setBrush( c );
940  p->drawPolygon( pa );
941  }
942  }
943 
944  control2[i] = true;
945  control2[mHalfEdge[i]->getNext()] = true;
946  control2[mHalfEdge[mHalfEdge[i]->getNext()]->getNext()] = true;
947  }//end of the section, which has to be removed later
948 
949  if ( control[i] )//check, if edge has already been drawn
950  {continue;}
951 
952  //draw the edge;
953  if ( halfEdgeBBoxTest( i, xlowleft, lowerborder, xupright, yupright ) )//only draw the halfedge if its bounding box intersects the painted area
954  {
955  if ( mHalfEdge[i]->getBreak() )//change the color it the edge is a breakline
956  {
957  p->setPen( mBreakEdgeColor );
958  }
959  else if ( mHalfEdge[i]->getForced() )//change the color if the edge is forced
960  {
961  p->setPen( mForcedEdgeColor );
962  }
963 
964 
965  p->drawLine( ( mPointVector[mHalfEdge[i]->getPoint()]->getX() - xlowleft ) / ( xupright - xlowleft )*width, ( yupright - mPointVector[mHalfEdge[i]->getPoint()]->getY() ) / ( xupright - xlowleft )*width, ( mPointVector[mHalfEdge[mHalfEdge[i]->getDual()]->getPoint()]->getX() - xlowleft ) / ( xupright - xlowleft )*width, ( yupright - mPointVector[mHalfEdge[mHalfEdge[i]->getDual()]->getPoint()]->getY() ) / ( xupright - xlowleft )*width );
966 
967  if ( mHalfEdge[i]->getForced() )
968  {
969  p->setPen( mEdgeColor );
970  }
971 
972 
973  }
974  control[i] = true;
975  control[mHalfEdge[i]->getDual()] = true;
976  }
977  }
978  else
979  {
980  double rightborder = width * ( yupright - ylowleft ) / height + xlowleft;//real world coordinates of the right widget border. This is useful to know because of the HalfEdge bounding box test
981  for ( unsigned int i = 0; i < mHalfEdge.count() - 1; i++ )
982  {
983  if ( mHalfEdge[i]->getPoint() == -1 || mHalfEdge[mHalfEdge[i]->getDual()]->getPoint() == -1 )
984  {continue;}
985 
986  //check, if the edge belongs to a flat triangle, remove this section later
987  if ( !control2[i] )
988  {
989  double p1, p2, p3;
990  if ( mHalfEdge[i]->getPoint() != -1 && mHalfEdge[mHalfEdge[i]->getNext()]->getPoint() != -1 && mHalfEdge[mHalfEdge[mHalfEdge[i]->getNext()]->getNext()]->getPoint() != -1 )
991  {
992  p1 = mPointVector[mHalfEdge[i]->getPoint()]->getZ();
993  p2 = mPointVector[mHalfEdge[mHalfEdge[i]->getNext()]->getPoint()]->getZ();
994  p3 = mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[i]->getNext()]->getNext()]->getPoint()]->getZ();
995  if ( p1 == p2 && p2 == p3 && halfEdgeBBoxTest( i, xlowleft, ylowleft, rightborder, yupright ) && halfEdgeBBoxTest( mHalfEdge[i]->getNext(), xlowleft, ylowleft, rightborder, yupright ) && halfEdgeBBoxTest( mHalfEdge[mHalfEdge[i]->getNext()]->getNext(), xlowleft, ylowleft, rightborder, yupright ) )//draw the triangle
996  {
997  QPointArray pa( 3 );
998  pa.setPoint( 0, ( mPointVector[mHalfEdge[i]->getPoint()]->getX() - xlowleft ) / ( yupright - ylowleft )*height, ( yupright - mPointVector[mHalfEdge[i]->getPoint()]->getY() ) / ( yupright - ylowleft )*height );
999  pa.setPoint( 1, ( mPointVector[mHalfEdge[mHalfEdge[i]->getNext()]->getPoint()]->getX() - xlowleft ) / ( yupright - ylowleft )*height, ( yupright - mPointVector[mHalfEdge[mHalfEdge[i]->getNext()]->getPoint()]->getY() ) / ( yupright - ylowleft )*height );
1000  pa.setPoint( 2, ( mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[i]->getNext()]->getNext()]->getPoint()]->getX() - xlowleft ) / ( yupright - ylowleft )*height, ( yupright - mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[i]->getNext()]->getNext()]->getPoint()]->getY() ) / ( yupright - ylowleft )*height );
1001  QColor c( 255, 0, 0 );
1002  p->setBrush( c );
1003  p->drawPolygon( pa );
1004  }
1005  }
1006 
1007  control2[i] = true;
1008  control2[mHalfEdge[i]->getNext()] = true;
1009  control2[mHalfEdge[mHalfEdge[i]->getNext()]->getNext()] = true;
1010  }//end of the section, which has to be removed later
1011 
1012 
1013  if ( control[i] )//check, if edge has already been drawn
1014  {continue;}
1015 
1016  //draw the edge
1017  if ( halfEdgeBBoxTest( i, xlowleft, ylowleft, rightborder, yupright ) )//only draw the edge if its bounding box intersects with the painted area
1018  {
1019  if ( mHalfEdge[i]->getBreak() )//change the color if the edge is a breakline
1020  {
1021  p->setPen( mBreakEdgeColor );
1022  }
1023  else if ( mHalfEdge[i]->getForced() )//change the color if the edge is forced
1024  {
1025  p->setPen( mForcedEdgeColor );
1026  }
1027 
1028  p->drawLine( ( mPointVector[mHalfEdge[i]->getPoint()]->getX() - xlowleft ) / ( yupright - ylowleft )*height, ( yupright - mPointVector[mHalfEdge[i]->getPoint()]->getY() ) / ( yupright - ylowleft )*height, ( mPointVector[mHalfEdge[mHalfEdge[i]->getDual()]->getPoint()]->getX() - xlowleft ) / ( yupright - ylowleft )*height, ( yupright - mPointVector[mHalfEdge[mHalfEdge[i]->getDual()]->getPoint()]->getY() ) / ( yupright - ylowleft )*height );
1029 
1030  if ( mHalfEdge[i]->getForced() )
1031  {
1032  p->setPen( mEdgeColor );
1033  }
1034 
1035  }
1036  control[i] = true;
1037  control[mHalfEdge[i]->getDual()] = true;
1038  }
1039  }
1040 
1041  delete[] control;
1042  delete[] control2;
1043 }
1044 #endif
1045 
1047 {
1048 
1049  //first find a half edge which points to p2
1050  const int firstedge = baseEdgeOfPoint( p2 );
1051 
1052  //then find the edge which comes from p1 or print an error message if there isn't any
1053  int theedge = -10;
1054  int nextnextedge = firstedge;
1055  int edge, nextedge;
1056  do
1057  {
1058  edge = mHalfEdge[nextnextedge]->getDual();
1059  if ( mHalfEdge[edge]->getPoint() == p1 )
1060  {
1061  theedge = nextnextedge;
1062  break;
1063  }//we found the edge
1064  nextedge = mHalfEdge[edge]->getNext();
1065  nextnextedge = mHalfEdge[nextedge]->getNext();
1066  }
1067  while ( nextnextedge != firstedge );
1068 
1069  if ( theedge == -10 )//there is no edge between p1 and p2
1070  {
1071  //QgsDebugMsg( QStringLiteral( "warning, error: the points are: %1 and %2" ).arg( p1 ).arg( p2 ) );
1072  return -10;
1073  }
1074 
1075  //finally find the opposite point
1076  return mHalfEdge[mHalfEdge[mHalfEdge[theedge]->getDual()]->getNext()]->getPoint();
1077 
1078 }
1079 
1081 {
1082  const int firstedge = baseEdgeOfPoint( pointno );
1083 
1084  QList<int> vlist;
1085  if ( firstedge == -1 )//an error occurred
1086  {
1087  return vlist;
1088  }
1089 
1090  int actedge = firstedge;
1091  int edge, nextedge, nextnextedge;
1092  do
1093  {
1094  edge = mHalfEdge[actedge]->getDual();
1095  vlist.append( mHalfEdge[edge]->getPoint() );//add the number of the endpoint of the first edge to the value list
1096  nextedge = mHalfEdge[edge]->getNext();
1097  vlist.append( mHalfEdge[nextedge]->getPoint() );//add the number of the endpoint of the second edge to the value list
1098  nextnextedge = mHalfEdge[nextedge]->getNext();
1099  vlist.append( mHalfEdge[nextnextedge]->getPoint() );//add the number of endpoint of the third edge to the value list
1100  if ( mHalfEdge[nextnextedge]->getBreak() )//add, whether the third edge is a breakline or not
1101  {
1102  vlist.append( -10 );
1103  }
1104  else
1105  {
1106  vlist.append( -20 );
1107  }
1108  actedge = nextnextedge;
1109  }
1110  while ( nextnextedge != firstedge );
1111 
1112  return vlist;
1113 
1114 }
1115 
1116 bool QgsDualEdgeTriangulation::triangleVertices( double x, double y, QgsPoint &p1, int &n1, QgsPoint &p2, int &n2, QgsPoint &p3, int &n3 )
1117 {
1118  if ( mPointVector.size() < 3 )
1119  {
1120  return false;
1121  }
1122 
1123  const QgsPoint point( x, y, 0 );
1124  const int edge = baseEdgeOfTriangle( point );
1125  if ( edge == -10 )//the point is outside the convex hull
1126  {
1127  return false;
1128  }
1129 
1130  else if ( edge >= 0 )//the point is inside the convex hull
1131  {
1132  const int ptnr1 = mHalfEdge[edge]->getPoint();
1133  const int ptnr2 = mHalfEdge[mHalfEdge[edge]->getNext()]->getPoint();
1134  const int ptnr3 = mHalfEdge[mHalfEdge[mHalfEdge[edge]->getNext()]->getNext()]->getPoint();
1135  p1.setX( mPointVector[ptnr1]->x() );
1136  p1.setY( mPointVector[ptnr1]->y() );
1137  p1.setZ( mPointVector[ptnr1]->z() );
1138  p2.setX( mPointVector[ptnr2]->x() );
1139  p2.setY( mPointVector[ptnr2]->y() );
1140  p2.setZ( mPointVector[ptnr2]->z() );
1141  p3.setX( mPointVector[ptnr3]->x() );
1142  p3.setY( mPointVector[ptnr3]->y() );
1143  p3.setZ( mPointVector[ptnr3]->z() );
1144  n1 = ptnr1;
1145  n2 = ptnr2;
1146  n3 = ptnr3;
1147  return true;
1148  }
1149  else if ( edge == -20 )//the point is exactly on an edge
1150  {
1151  const int ptnr1 = mHalfEdge[mEdgeWithPoint]->getPoint();
1152  const int ptnr2 = mHalfEdge[mHalfEdge[mEdgeWithPoint]->getNext()]->getPoint();
1153  const int ptnr3 = mHalfEdge[mHalfEdge[mHalfEdge[mEdgeWithPoint]->getNext()]->getNext()]->getPoint();
1154  if ( ptnr1 == -1 || ptnr2 == -1 || ptnr3 == -1 )
1155  {
1156  return false;
1157  }
1158  p1.setX( mPointVector[ptnr1]->x() );
1159  p1.setY( mPointVector[ptnr1]->y() );
1160  p1.setZ( mPointVector[ptnr1]->z() );
1161  p2.setX( mPointVector[ptnr2]->x() );
1162  p2.setY( mPointVector[ptnr2]->y() );
1163  p2.setZ( mPointVector[ptnr2]->z() );
1164  p3.setX( mPointVector[ptnr3]->x() );
1165  p3.setY( mPointVector[ptnr3]->y() );
1166  p3.setZ( mPointVector[ptnr3]->z() );
1167  n1 = ptnr1;
1168  n2 = ptnr2;
1169  n3 = ptnr3;
1170  return true;
1171  }
1172  else if ( edge == -25 )//x and y are the coordinates of an existing point
1173  {
1174  const int edge1 = baseEdgeOfPoint( mTwiceInsPoint );
1175  const int edge2 = mHalfEdge[edge1]->getNext();
1176  const int edge3 = mHalfEdge[edge2]->getNext();
1177  const int ptnr1 = mHalfEdge[edge1]->getPoint();
1178  const int ptnr2 = mHalfEdge[edge2]->getPoint();
1179  const int ptnr3 = mHalfEdge[edge3]->getPoint();
1180  p1.setX( mPointVector[ptnr1]->x() );
1181  p1.setY( mPointVector[ptnr1]->y() );
1182  p1.setZ( mPointVector[ptnr1]->z() );
1183  p2.setX( mPointVector[ptnr2]->x() );
1184  p2.setY( mPointVector[ptnr2]->y() );
1185  p2.setZ( mPointVector[ptnr2]->z() );
1186  p3.setX( mPointVector[ptnr3]->x() );
1187  p3.setY( mPointVector[ptnr3]->y() );
1188  p3.setZ( mPointVector[ptnr3]->z() );
1189  n1 = ptnr1;
1190  n2 = ptnr2;
1191  n3 = ptnr3;
1192  return true;
1193  }
1194  else if ( edge == -5 )//numerical problems in 'baseEdgeOfTriangle'
1195  {
1196  const int ptnr1 = mHalfEdge[mUnstableEdge]->getPoint();
1197  const int ptnr2 = mHalfEdge[mHalfEdge[mUnstableEdge]->getNext()]->getPoint();
1198  const int ptnr3 = mHalfEdge[mHalfEdge[mHalfEdge[mUnstableEdge]->getNext()]->getNext()]->getPoint();
1199  if ( ptnr1 == -1 || ptnr2 == -1 || ptnr3 == -1 )
1200  {
1201  return false;
1202  }
1203  p1.setX( mPointVector[ptnr1]->x() );
1204  p1.setY( mPointVector[ptnr1]->y() );
1205  p1.setZ( mPointVector[ptnr1]->z() );
1206  p2.setX( mPointVector[ptnr2]->x() );
1207  p2.setY( mPointVector[ptnr2]->y() );
1208  p2.setZ( mPointVector[ptnr2]->z() );
1209  p3.setX( mPointVector[ptnr3]->x() );
1210  p3.setY( mPointVector[ptnr3]->y() );
1211  p3.setZ( mPointVector[ptnr3]->z() );
1212  n1 = ptnr1;
1213  n2 = ptnr2;
1214  n3 = ptnr3;
1215  return true;
1216  }
1217  else//problems
1218  {
1219  //QgsDebugMsg( QStringLiteral( "problem: the edge is: %1" ).arg( edge ) );
1220  return false;
1221  }
1222 }
1223 
1225 {
1226  if ( mPointVector.size() < 3 )
1227  {
1228  return false;
1229  }
1230 
1231  const QgsPoint point( x, y, 0 );
1232  const int edge = baseEdgeOfTriangle( point );
1233  if ( edge == -10 )//the point is outside the convex hull
1234  {
1235  return false;
1236  }
1237  else if ( edge >= 0 )//the point is inside the convex hull
1238  {
1239  const int ptnr1 = mHalfEdge[edge]->getPoint();
1240  const int ptnr2 = mHalfEdge[mHalfEdge[edge]->getNext()]->getPoint();
1241  const int ptnr3 = mHalfEdge[mHalfEdge[mHalfEdge[edge]->getNext()]->getNext()]->getPoint();
1242  p1.setX( mPointVector[ptnr1]->x() );
1243  p1.setY( mPointVector[ptnr1]->y() );
1244  p1.setZ( mPointVector[ptnr1]->z() );
1245  p2.setX( mPointVector[ptnr2]->x() );
1246  p2.setY( mPointVector[ptnr2]->y() );
1247  p2.setZ( mPointVector[ptnr2]->z() );
1248  p3.setX( mPointVector[ptnr3]->x() );
1249  p3.setY( mPointVector[ptnr3]->y() );
1250  p3.setZ( mPointVector[ptnr3]->z() );
1251  return true;
1252  }
1253  else if ( edge == -20 )//the point is exactly on an edge
1254  {
1255  const int ptnr1 = mHalfEdge[mEdgeWithPoint]->getPoint();
1256  const int ptnr2 = mHalfEdge[mHalfEdge[mEdgeWithPoint]->getNext()]->getPoint();
1257  const int ptnr3 = mHalfEdge[mHalfEdge[mHalfEdge[mEdgeWithPoint]->getNext()]->getNext()]->getPoint();
1258  if ( ptnr1 == -1 || ptnr2 == -1 || ptnr3 == -1 )
1259  {
1260  return false;
1261  }
1262  p1.setX( mPointVector[ptnr1]->x() );
1263  p1.setY( mPointVector[ptnr1]->y() );
1264  p1.setZ( mPointVector[ptnr1]->z() );
1265  p2.setX( mPointVector[ptnr2]->x() );
1266  p2.setY( mPointVector[ptnr2]->y() );
1267  p2.setZ( mPointVector[ptnr2]->z() );
1268  p3.setX( mPointVector[ptnr3]->x() );
1269  p3.setY( mPointVector[ptnr3]->y() );
1270  p3.setZ( mPointVector[ptnr3]->z() );
1271  return true;
1272  }
1273  else if ( edge == -25 )//x and y are the coordinates of an existing point
1274  {
1275  const int edge1 = baseEdgeOfPoint( mTwiceInsPoint );
1276  const int edge2 = mHalfEdge[edge1]->getNext();
1277  const int edge3 = mHalfEdge[edge2]->getNext();
1278  const int ptnr1 = mHalfEdge[edge1]->getPoint();
1279  const int ptnr2 = mHalfEdge[edge2]->getPoint();
1280  const int ptnr3 = mHalfEdge[edge3]->getPoint();
1281  if ( ptnr1 == -1 || ptnr2 == -1 || ptnr3 == -1 )
1282  {
1283  return false;
1284  }
1285  p1.setX( mPointVector[ptnr1]->x() );
1286  p1.setY( mPointVector[ptnr1]->y() );
1287  p1.setZ( mPointVector[ptnr1]->z() );
1288  p2.setX( mPointVector[ptnr2]->x() );
1289  p2.setY( mPointVector[ptnr2]->y() );
1290  p2.setZ( mPointVector[ptnr2]->z() );
1291  p3.setX( mPointVector[ptnr3]->x() );
1292  p3.setY( mPointVector[ptnr3]->y() );
1293  p3.setZ( mPointVector[ptnr3]->z() );
1294  return true;
1295  }
1296  else if ( edge == -5 )//numerical problems in 'baseEdgeOfTriangle'
1297  {
1298  const int ptnr1 = mHalfEdge[mUnstableEdge]->getPoint();
1299  const int ptnr2 = mHalfEdge[mHalfEdge[mUnstableEdge]->getNext()]->getPoint();
1300  const int ptnr3 = mHalfEdge[mHalfEdge[mHalfEdge[mUnstableEdge]->getNext()]->getNext()]->getPoint();
1301  if ( ptnr1 == -1 || ptnr2 == -1 || ptnr3 == -1 )
1302  {
1303  return false;
1304  }
1305  p1.setX( mPointVector[ptnr1]->x() );
1306  p1.setY( mPointVector[ptnr1]->y() );
1307  p1.setZ( mPointVector[ptnr1]->z() );
1308  p2.setX( mPointVector[ptnr2]->x() );
1309  p2.setY( mPointVector[ptnr2]->y() );
1310  p2.setZ( mPointVector[ptnr2]->z() );
1311  p3.setX( mPointVector[ptnr3]->x() );
1312  p3.setY( mPointVector[ptnr3]->y() );
1313  p3.setZ( mPointVector[ptnr3]->z() );
1314  return true;
1315  }
1316  else//problems
1317  {
1318  return false;
1319  }
1320 }
1321 
1322 unsigned int QgsDualEdgeTriangulation::insertEdge( int dual, int next, int point, bool mbreak, bool forced )
1323 {
1324  HalfEdge *edge = new HalfEdge( dual, next, point, mbreak, forced );
1325  mHalfEdge.append( edge );
1326  return mHalfEdge.count() - 1;
1327 
1328 }
1329 
1330 static bool altitudeTriangleIsSmall( const QgsPoint &pointBase1, const QgsPoint &pointBase2, const QgsPoint &pt3, double tolerance )
1331 {
1332  // Compare the altitude of the triangle defined by base points and a third point with tolerance. Return true if the altitude < tolerance
1333  const double x1 = pointBase1.x();
1334  const double y1 = pointBase1.y();
1335  const double x2 = pointBase2.x();
1336  const double y2 = pointBase2.y();
1337  const double X = pt3.x();
1338  const double Y = pt3.y();
1339  QgsPoint projectedPoint;
1340 
1341  double nx, ny; //normal vector
1342 
1343  nx = y2 - y1;
1344  ny = -( x2 - x1 );
1345 
1346  double t;
1347  t = ( X * ny - Y * nx - x1 * ny + y1 * nx ) / ( ( x2 - x1 ) * ny - ( y2 - y1 ) * nx );
1348  projectedPoint.setX( x1 + t * ( x2 - x1 ) );
1349  projectedPoint.setY( y1 + t * ( y2 - y1 ) );
1350 
1351  return pt3.distance( projectedPoint ) < tolerance;
1352 }
1353 
1354 int QgsDualEdgeTriangulation::insertForcedSegment( int p1, int p2, QgsInterpolator::SourceType segmentType )
1355 {
1356  if ( p1 == p2 )
1357  {
1358  return 0;
1359  }
1360 
1361  QgsPoint *point1 = mPointVector.at( p1 );
1362  QgsPoint *point2 = mPointVector.at( p2 );
1363 
1364  //list with (half of) the crossed edges
1365  QList<int> crossedEdges;
1366 
1367  //an edge pointing to p1
1368  int pointingEdge = 0;
1369 
1370  pointingEdge = baseEdgeOfPoint( p1 );
1371 
1372  if ( pointingEdge == -1 )
1373  {
1374  return -100;//return an error code
1375  }
1376 
1377  //go around p1 and find out, if the segment already exists and if not, which is the first cutted edge
1378  int actEdge = mHalfEdge[pointingEdge]->getDual();
1379  const int firstActEdge = actEdge;
1380  //number to prevent endless loops
1381  int control = 0;
1382 
1383  do //if it's an endless loop, something went wrong
1384  {
1385  control += 1;
1386  if ( control > 100 )
1387  {
1388  //QgsDebugMsg( QStringLiteral( "warning, endless loop" ) );
1389  return -100;//return an error code
1390  }
1391 
1392  if ( mHalfEdge[actEdge]->getPoint() == -1 )//actEdge points to the virtual point
1393  {
1394  actEdge = mHalfEdge[mHalfEdge[mHalfEdge[actEdge]->getNext()]->getNext()]->getDual();
1395  continue;
1396  }
1397 
1398  //test, if actEdge is already the forced edge
1399  if ( mHalfEdge[actEdge]->getPoint() == p2 )
1400  {
1401  mHalfEdge[actEdge]->setForced( true );
1402  mHalfEdge[actEdge]->setBreak( segmentType == QgsInterpolator::SourceBreakLines );
1403  mHalfEdge[mHalfEdge[actEdge]->getDual()]->setForced( true );
1404  mHalfEdge[mHalfEdge[actEdge]->getDual()]->setBreak( segmentType == QgsInterpolator::SourceBreakLines );
1405  return actEdge;
1406  }
1407 
1408  //test, if the forced segment is a multiple of actEdge and if the direction is the same
1409  if ( /*lines are parallel*/( point2->y() - point1->y() ) / ( mPointVector[mHalfEdge[actEdge]->getPoint()]->y() - point1->y() ) == ( point2->x() - point1->x() ) / ( mPointVector[mHalfEdge[actEdge]->getPoint()]->x() - point1->x() )
1410  && ( ( point2->y() - point1->y() ) >= 0 ) == ( ( mPointVector[mHalfEdge[actEdge]->getPoint()]->y() - point1->y() ) > 0 )
1411  && ( ( point2->x() - point1->x() ) >= 0 ) == ( ( mPointVector[mHalfEdge[actEdge]->getPoint()]->x() - point1->x() ) > 0 ) )
1412  {
1413  //mark actedge and Dual(actedge) as forced, reset p1 and start the method from the beginning
1414  mHalfEdge[actEdge]->setForced( true );
1415  mHalfEdge[actEdge]->setBreak( segmentType == QgsInterpolator::SourceBreakLines );
1416  mHalfEdge[mHalfEdge[actEdge]->getDual()]->setForced( true );
1417  mHalfEdge[mHalfEdge[actEdge]->getDual()]->setBreak( segmentType == QgsInterpolator::SourceBreakLines );
1418  const int a = insertForcedSegment( mHalfEdge[actEdge]->getPoint(), p2, segmentType );
1419  return a;
1420  }
1421 
1422  //test, if the forced segment intersects Next(actEdge)
1423  const int oppositeEdge = mHalfEdge[actEdge]->getNext();
1424 
1425  if ( mHalfEdge[oppositeEdge]->getPoint() == -1 || mHalfEdge[mHalfEdge[oppositeEdge]->getDual()]->getPoint() == -1 ) //intersection with line to the virtual point makes no sense
1426  {
1427  actEdge = mHalfEdge[mHalfEdge[oppositeEdge]->getNext()]->getDual(); //continue with the next edge around p1
1428  continue;
1429  }
1430 
1431  QgsPoint *oppositePoint1 = mPointVector[mHalfEdge[oppositeEdge]->getPoint()];
1432  QgsPoint *oppositePoint2 = mPointVector[mHalfEdge[mHalfEdge[oppositeEdge]->getDual()]->getPoint()];
1433 
1434  if ( altitudeTriangleIsSmall( *oppositePoint1, *oppositePoint2, *point1, oppositePoint1->distance( *oppositePoint2 ) / 500 ) )
1435  {
1436  // to much risks to do something, go away
1437  return -100;
1438  }
1439 
1440  if ( MathUtils::lineIntersection( point1,
1441  point2,
1442  mPointVector[mHalfEdge[oppositeEdge]->getPoint()],
1443  mPointVector[mHalfEdge[mHalfEdge[oppositeEdge]->getDual()]->getPoint()] ) )
1444  {
1445  if ( mHalfEdge[mHalfEdge[actEdge]->getNext()]->getForced() && mForcedCrossBehavior == QgsTriangulation::SnappingTypeVertex )//if the crossed edge is a forced edge, we have to snap the forced line to the next node
1446  {
1447  QgsPoint crosspoint( 0, 0, 0 );
1448  int p3, p4;
1449  p3 = mHalfEdge[mHalfEdge[actEdge]->getNext()]->getPoint();
1450  p4 = mHalfEdge[mHalfEdge[mHalfEdge[actEdge]->getNext()]->getDual()]->getPoint();
1451  MathUtils::lineIntersection( point1, point2, mPointVector[p3], mPointVector[p4], &crosspoint );
1452  const double dista = std::sqrt( ( crosspoint.x() - mPointVector[p3]->x() ) * ( crosspoint.x() - mPointVector[p3]->x() ) + ( crosspoint.y() - mPointVector[p3]->y() ) * ( crosspoint.y() - mPointVector[p3]->y() ) );
1453  const double distb = std::sqrt( ( crosspoint.x() - mPointVector[p4]->x() ) * ( crosspoint.x() - mPointVector[p4]->x() ) + ( crosspoint.y() - mPointVector[p4]->y() ) * ( crosspoint.y() - mPointVector[p4]->y() ) );
1454  if ( dista <= distb )
1455  {
1456  insertForcedSegment( p1, p3, segmentType );
1457  const int e = insertForcedSegment( p3, p2, segmentType );
1458  return e;
1459  }
1460  else if ( distb <= dista )
1461  {
1462  insertForcedSegment( p1, p4, segmentType );
1463  const int e = insertForcedSegment( p4, p2, segmentType );
1464  return e;
1465  }
1466  }
1467 
1468  if ( mHalfEdge[mHalfEdge[actEdge]->getNext()]->getForced() && mForcedCrossBehavior == QgsTriangulation::InsertVertex )//if the crossed edge is a forced edge, we have to insert a new vertice on this edge
1469  {
1470  QgsPoint crosspoint( 0, 0, 0 );
1471  int p3, p4;
1472  p3 = mHalfEdge[mHalfEdge[actEdge]->getNext()]->getPoint();
1473  p4 = mHalfEdge[mHalfEdge[mHalfEdge[actEdge]->getNext()]->getDual()]->getPoint();
1474  MathUtils::lineIntersection( mPointVector[p1], mPointVector[p2], mPointVector[p3], mPointVector[p4], &crosspoint );
1475  const double distpart = std::sqrt( ( crosspoint.x() - mPointVector[p4]->x() ) * ( crosspoint.x() - mPointVector[p4]->x() ) + ( crosspoint.y() - mPointVector[p4]->y() ) * ( crosspoint.y() - mPointVector[p4]->y() ) );
1476  const double disttot = std::sqrt( ( mPointVector[p3]->x() - mPointVector[p4]->x() ) * ( mPointVector[p3]->x() - mPointVector[p4]->x() ) + ( mPointVector[p3]->y() - mPointVector[p4]->y() ) * ( mPointVector[p3]->y() - mPointVector[p4]->y() ) );
1477  const float frac = distpart / disttot;
1478 
1479  if ( frac == 0 || frac == 1 )//just in case MathUtils::lineIntersection does not work as expected
1480  {
1481  if ( frac == 0 )
1482  {
1483  //mark actEdge and Dual(actEdge) as forced, reset p1 and start the method from the beginning
1484  mHalfEdge[actEdge]->setForced( true );
1485  mHalfEdge[actEdge]->setBreak( segmentType == QgsInterpolator::SourceBreakLines );
1486  mHalfEdge[mHalfEdge[actEdge]->getDual()]->setForced( true );
1487  mHalfEdge[mHalfEdge[actEdge]->getDual()]->setBreak( segmentType == QgsInterpolator::SourceBreakLines );
1488  const int a = insertForcedSegment( p4, p2, segmentType );
1489  return a;
1490  }
1491  else if ( frac == 1 )
1492  {
1493  //mark actEdge and Dual(actEdge) as forced, reset p1 and start the method from the beginning
1494  mHalfEdge[actEdge]->setForced( true );
1495  mHalfEdge[actEdge]->setBreak( segmentType == QgsInterpolator::SourceBreakLines );
1496  mHalfEdge[mHalfEdge[actEdge]->getDual()]->setForced( true );
1497  mHalfEdge[mHalfEdge[actEdge]->getDual()]->setBreak( segmentType == QgsInterpolator::SourceBreakLines );
1498  if ( p3 != p2 )
1499  {
1500  const int a = insertForcedSegment( p3, p2, segmentType );
1501  return a;
1502  }
1503  else
1504  {
1505  return actEdge;
1506  }
1507  }
1508 
1509  }
1510  else
1511  {
1512  const int newpoint = splitHalfEdge( mHalfEdge[actEdge]->getNext(), frac );
1513  insertForcedSegment( p1, newpoint, segmentType );
1514  const int e = insertForcedSegment( newpoint, p2, segmentType );
1515  return e;
1516  }
1517  }
1518 
1519  //add the first HalfEdge to the list of crossed edges
1520  crossedEdges.append( oppositeEdge );
1521  break;
1522  }
1523  actEdge = mHalfEdge[mHalfEdge[oppositeEdge]->getNext()]->getDual(); //continue with the next edge around p1
1524  }
1525  while ( actEdge != firstActEdge );
1526 
1527  if ( crossedEdges.isEmpty() ) //nothing found due to rounding error, better to go away!
1528  return -100;
1529  int lastEdgeOppositePointIndex = mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getPoint();
1530 
1531  //we found the first edge, terminated the method or called the method with other points. Lets search for all the other crossed edges
1532  while ( lastEdgeOppositePointIndex != p2 )
1533  {
1534  QgsPoint *lastEdgePoint1 = mPointVector[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getPoint()];
1535  QgsPoint *lastEdgePoint2 = mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getNext()]->getPoint()];
1536  QgsPoint *lastEdgeOppositePoint = mPointVector[lastEdgeOppositePointIndex];
1537 
1538  if ( MathUtils::lineIntersection( lastEdgePoint1, lastEdgeOppositePoint,
1539  point1, point2 ) )
1540  {
1541  if ( mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getForced() && mForcedCrossBehavior == QgsTriangulation::SnappingTypeVertex )//if the crossed edge is a forced edge and mForcedCrossBehavior is SnappingType_VERTICE, we have to snap the forced line to the next node
1542  {
1543  QgsPoint crosspoint( 0, 0, 0 );
1544  int p3, p4;
1545  p3 = mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getPoint();
1546  p4 = mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getPoint();
1547  MathUtils::lineIntersection( point1, point2, mPointVector[p3], mPointVector[p4], &crosspoint );
1548  const double dista = std::sqrt( ( crosspoint.x() - mPointVector[p3]->x() ) * ( crosspoint.x() - mPointVector[p3]->x() ) + ( crosspoint.y() - mPointVector[p3]->y() ) * ( crosspoint.y() - mPointVector[p3]->y() ) );
1549  const double distb = std::sqrt( ( crosspoint.x() - mPointVector[p4]->x() ) * ( crosspoint.x() - mPointVector[p4]->x() ) + ( crosspoint.y() - mPointVector[p4]->y() ) * ( crosspoint.y() - mPointVector[p4]->y() ) );
1550  if ( dista <= distb )
1551  {
1552  insertForcedSegment( p1, p3, segmentType );
1553  const int e = insertForcedSegment( p3, p2, segmentType );
1554  return e;
1555  }
1556  else if ( distb <= dista )
1557  {
1558  insertForcedSegment( p1, p4, segmentType );
1559  const int e = insertForcedSegment( p4, p2, segmentType );
1560  return e;
1561  }
1562  }
1563  else if ( mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getForced() && mForcedCrossBehavior == QgsTriangulation::InsertVertex )//if the crossed edge is a forced edge, we have to insert a new vertice on this edge
1564  {
1565  QgsPoint crosspoint( 0, 0, 0 );
1566  int p3, p4;
1567  p3 = mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getPoint();
1568  p4 = mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getPoint();
1569  MathUtils::lineIntersection( point1, point2, mPointVector[p3], mPointVector[p4], &crosspoint );
1570  const double distpart = std::sqrt( ( crosspoint.x() - mPointVector[p3]->x() ) * ( crosspoint.x() - mPointVector[p3]->x() ) + ( crosspoint.y() - mPointVector[p3]->y() ) * ( crosspoint.y() - mPointVector[p3]->y() ) );
1571  const double disttot = std::sqrt( ( mPointVector[p3]->x() - mPointVector[p4]->x() ) * ( mPointVector[p3]->x() - mPointVector[p4]->x() ) + ( mPointVector[p3]->y() - mPointVector[p4]->y() ) * ( mPointVector[p3]->y() - mPointVector[p4]->y() ) );
1572  const float frac = distpart / disttot;
1573  if ( frac == 0 || frac == 1 )
1574  {
1575  break;//seems that a roundoff error occurred. We found the endpoint
1576  }
1577  const int newpoint = splitHalfEdge( mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext(), frac );
1578  insertForcedSegment( p1, newpoint, segmentType );
1579  const int e = insertForcedSegment( newpoint, p2, segmentType );
1580  return e;
1581  }
1582 
1583  crossedEdges.append( mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext() );
1584  }
1585  else if ( MathUtils::lineIntersection( lastEdgePoint2, lastEdgeOppositePoint,
1586  point1, point2 ) )
1587  {
1588  if ( mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getNext()]->getForced() && mForcedCrossBehavior == QgsTriangulation::SnappingTypeVertex )//if the crossed edge is a forced edge and mForcedCrossBehavior is SnappingType_VERTICE, we have to snap the forced line to the next node
1589  {
1590  QgsPoint crosspoint( 0, 0, 0 );
1591  int p3, p4;
1592  p3 = mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getPoint();
1593  p4 = mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getNext()]->getPoint();
1594  MathUtils::lineIntersection( mPointVector[p1], mPointVector[p2], mPointVector[p3], mPointVector[p4], &crosspoint );
1595  const double dista = std::sqrt( ( crosspoint.x() - mPointVector[p3]->x() ) * ( crosspoint.x() - mPointVector[p3]->x() ) + ( crosspoint.y() - mPointVector[p3]->y() ) * ( crosspoint.y() - mPointVector[p3]->y() ) );
1596  const double distb = std::sqrt( ( crosspoint.x() - mPointVector[p4]->x() ) * ( crosspoint.x() - mPointVector[p4]->x() ) + ( crosspoint.y() - mPointVector[p4]->y() ) * ( crosspoint.y() - mPointVector[p4]->y() ) );
1597  if ( dista <= distb )
1598  {
1599  insertForcedSegment( p1, p3, segmentType );
1600  const int e = insertForcedSegment( p3, p2, segmentType );
1601  return e;
1602  }
1603  else if ( distb <= dista )
1604  {
1605  insertForcedSegment( p1, p4, segmentType );
1606  const int e = insertForcedSegment( p4, p2, segmentType );
1607  return e;
1608  }
1609  }
1610  else if ( mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getNext()]->getForced() && mForcedCrossBehavior == QgsTriangulation::InsertVertex )//if the crossed edge is a forced edge, we have to insert a new vertice on this edge
1611  {
1612  QgsPoint crosspoint( 0, 0, 0 );
1613  int p3, p4;
1614  p3 = mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getPoint();
1615  p4 = mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getNext()]->getPoint();
1616  MathUtils::lineIntersection( point1, point2, mPointVector[p3], mPointVector[p4], &crosspoint );
1617  const double distpart = std::sqrt( ( crosspoint.x() - mPointVector[p3]->x() ) * ( crosspoint.x() - mPointVector[p3]->x() ) + ( crosspoint.y() - mPointVector[p3]->y() ) * ( crosspoint.y() - mPointVector[p3]->y() ) );
1618  const double disttot = std::sqrt( ( mPointVector[p3]->x() - mPointVector[p4]->x() ) * ( mPointVector[p3]->x() - mPointVector[p4]->x() ) + ( mPointVector[p3]->y() - mPointVector[p4]->y() ) * ( mPointVector[p3]->y() - mPointVector[p4]->y() ) );
1619  const float frac = distpart / disttot;
1620  if ( frac == 0 || frac == 1 )
1621  {
1622  break;//seems that a roundoff error occurred. We found the endpoint
1623  }
1624  const int newpoint = splitHalfEdge( mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getNext(), frac );
1625  insertForcedSegment( p1, newpoint, segmentType );
1626  const int e = insertForcedSegment( newpoint, p2, segmentType );
1627  return e;
1628  }
1629 
1630  crossedEdges.append( mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getNext() );
1631  }
1632  else
1633  {
1634  //no intersection found, surely due to rounding error or something else wrong, better to give up!
1635  return -100;
1636  }
1637  lastEdgeOppositePointIndex = mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getPoint();
1638  }
1639 
1640  // Last check before construct the breakline
1641  QgsPoint *lastEdgePoint1 = mPointVector[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getPoint()];
1642  QgsPoint *lastEdgePoint2 = mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getNext()]->getPoint()];
1643  QgsPoint *lastEdgeOppositePoint = mPointVector[lastEdgeOppositePointIndex];
1644  if ( altitudeTriangleIsSmall( *lastEdgePoint1, *lastEdgePoint2, *lastEdgeOppositePoint, lastEdgePoint1->distance( *lastEdgePoint2 ) / 500 ) )
1645  return -100;
1646 
1647  //set the flags 'forced' and 'break' to false for every edge and dualedge of 'crossEdges'
1648  QList<int>::const_iterator iter;
1649  for ( iter = crossedEdges.constBegin(); iter != crossedEdges.constEnd(); ++iter )
1650  {
1651  mHalfEdge[( *( iter ) )]->setForced( false );
1652  mHalfEdge[( *( iter ) )]->setBreak( false );
1653  mHalfEdge[mHalfEdge[( *( iter ) )]->getDual()]->setForced( false );
1654  mHalfEdge[mHalfEdge[( *( iter ) )]->getDual()]->setBreak( false );
1655  }
1656 
1657  //crossed edges is filled, now the two polygons to be retriangulated can be build
1658 
1659  QList<int> freelist = crossedEdges;//copy the list with the crossed edges to remove the edges already reused
1660 
1661  //create the left polygon as a list of the numbers of the halfedges
1662  QList<int> leftPolygon;
1663  QList<int> rightPolygon;
1664 
1665  //insert the forced edge and enter the corresponding halfedges as the first edges in the left and right polygons. The nexts and points are set later because of the algorithm to build two polygons from 'crossedEdges'
1666  const int firstedge = freelist.first();//edge pointing from p1 to p2
1667  mHalfEdge[firstedge]->setForced( true );
1668  mHalfEdge[firstedge]->setBreak( segmentType == QgsInterpolator::SourceBreakLines );
1669  leftPolygon.append( firstedge );
1670  const int dualfirstedge = mHalfEdge[freelist.first()]->getDual();//edge pointing from p2 to p1
1671  mHalfEdge[dualfirstedge]->setForced( true );
1672  mHalfEdge[dualfirstedge]->setBreak( segmentType == QgsInterpolator::SourceBreakLines );
1673  rightPolygon.append( dualfirstedge );
1674  freelist.pop_front();//delete the first entry from the freelist
1675 
1676  //finish the polygon on the left side
1677  int actpointl = p2;
1678  QList<int>::const_iterator leftiter; //todo: is there a better way to set an iterator to the last list element?
1679  leftiter = crossedEdges.constEnd();
1680  --leftiter;
1681  while ( true )
1682  {
1683  const int newpoint = mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[( *leftiter )]->getDual()]->getNext()]->getNext()]->getPoint();
1684  if ( newpoint != actpointl )
1685  {
1686  //insert the edge into the leftPolygon
1687  actpointl = newpoint;
1688  const int theedge = mHalfEdge[mHalfEdge[mHalfEdge[( *leftiter )]->getDual()]->getNext()]->getNext();
1689  leftPolygon.append( theedge );
1690  }
1691  if ( leftiter == crossedEdges.constBegin() )
1692  {break;}
1693  --leftiter;
1694  }
1695 
1696  //insert the last element into leftPolygon
1697  leftPolygon.append( mHalfEdge[crossedEdges.first()]->getNext() );
1698 
1699  //finish the polygon on the right side
1700  QList<int>::const_iterator rightiter;
1701  int actpointr = p1;
1702  for ( rightiter = crossedEdges.constBegin(); rightiter != crossedEdges.constEnd(); ++rightiter )
1703  {
1704  const int newpoint = mHalfEdge[mHalfEdge[mHalfEdge[( *rightiter )]->getNext()]->getNext()]->getPoint();
1705  if ( newpoint != actpointr )
1706  {
1707  //insert the edge into the right polygon
1708  actpointr = newpoint;
1709  const int theedge = mHalfEdge[mHalfEdge[( *rightiter )]->getNext()]->getNext();
1710  rightPolygon.append( theedge );
1711  }
1712  }
1713 
1714 
1715  //insert the last element into rightPolygon
1716  rightPolygon.append( mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext() );
1717  mHalfEdge[rightPolygon.last()]->setNext( dualfirstedge );//set 'Next' of the last edge to dualfirstedge
1718 
1719  //set the necessary nexts of leftPolygon(except the first)
1720  int actedgel = leftPolygon[1];
1721  leftiter = leftPolygon.constBegin();
1722  leftiter += 2;
1723  for ( ; leftiter != leftPolygon.constEnd(); ++leftiter )
1724  {
1725  mHalfEdge[actedgel]->setNext( ( *leftiter ) );
1726  actedgel = ( *leftiter );
1727  }
1728 
1729  //set all the necessary nexts of rightPolygon
1730  int actedger = rightPolygon[1];
1731  rightiter = rightPolygon.constBegin();
1732  rightiter += 2;
1733  for ( ; rightiter != rightPolygon.constEnd(); ++rightiter )
1734  {
1735  mHalfEdge[actedger]->setNext( ( *rightiter ) );
1736  actedger = ( *( rightiter ) );
1737  }
1738 
1739 
1740  //setNext and setPoint for the forced edge because this would disturb the building of 'leftpoly' and 'rightpoly' otherwise
1741  mHalfEdge[leftPolygon.first()]->setNext( ( *( ++( leftiter = leftPolygon.constBegin() ) ) ) );
1742  mHalfEdge[leftPolygon.first()]->setPoint( p2 );
1743  mHalfEdge[leftPolygon.last()]->setNext( firstedge );
1744  mHalfEdge[rightPolygon.first()]->setNext( ( *( ++( rightiter = rightPolygon.constBegin() ) ) ) );
1745  mHalfEdge[rightPolygon.first()]->setPoint( p1 );
1746  mHalfEdge[rightPolygon.last()]->setNext( dualfirstedge );
1747 
1748  triangulatePolygon( &leftPolygon, &freelist, firstedge );
1749  triangulatePolygon( &rightPolygon, &freelist, dualfirstedge );
1750 
1751  //optimisation of the new edges
1752  for ( iter = crossedEdges.constBegin(); iter != crossedEdges.constEnd(); ++iter )
1753  {
1754  checkSwapRecursively( ( *( iter ) ), 0 );
1755  }
1756 
1757  return leftPolygon.first();
1758 }
1759 
1761 {
1762  mForcedCrossBehavior = b;
1763 }
1764 
1766 {
1767  mTriangleInterpolator = interpolator;
1768 }
1769 
1771 {
1772  //QgsDebugMsg( QStringLiteral( "am in eliminateHorizontalTriangles" ) );
1773  const double minangle = 0;//minimum angle for swapped triangles. If triangles generated by a swap would have a minimum angle (in degrees) below that value, the swap will not be done.
1774 
1775  while ( true )
1776  {
1777  bool swapped = false;//flag which allows exiting the loop
1778  bool *control = new bool[mHalfEdge.count()];//controlarray
1779 
1780  for ( int i = 0; i <= mHalfEdge.count() - 1; i++ )
1781  {
1782  control[i] = false;
1783  }
1784 
1785 
1786  for ( int i = 0; i <= mHalfEdge.count() - 1; i++ )
1787  {
1788  if ( control[i] )//edge has already been examined
1789  {
1790  continue;
1791  }
1792 
1793  int e1, e2, e3;//numbers of the three edges
1794  e1 = i;
1795  e2 = mHalfEdge[e1]->getNext();
1796  e3 = mHalfEdge[e2]->getNext();
1797 
1798  int p1, p2, p3;//numbers of the three points
1799  p1 = mHalfEdge[e1]->getPoint();
1800  p2 = mHalfEdge[e2]->getPoint();
1801  p3 = mHalfEdge[e3]->getPoint();
1802 
1803  //skip the iteration, if one point is the virtual point
1804  if ( p1 == -1 || p2 == -1 || p3 == -1 )
1805  {
1806  control[e1] = true;
1807  control[e2] = true;
1808  control[e3] = true;
1809  continue;
1810  }
1811 
1812  double el1, el2, el3;//elevations of the points
1813  el1 = mPointVector[p1]->z();
1814  el2 = mPointVector[p2]->z();
1815  el3 = mPointVector[p3]->z();
1816 
1817  if ( el1 == el2 && el2 == el3 )//we found a horizontal triangle
1818  {
1819  //swap edges if it is possible, if it would remove the horizontal triangle and if the minimum angle generated by the swap is high enough
1820  if ( swapPossible( ( uint )e1 ) && mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[e1]->getDual()]->getNext()]->getPoint()]->z() != el1 && swapMinAngle( e1 ) > minangle )
1821  {
1822  doOnlySwap( ( uint )e1 );
1823  swapped = true;
1824  }
1825  else if ( swapPossible( ( uint )e2 ) && mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[e2]->getDual()]->getNext()]->getPoint()]->z() != el2 && swapMinAngle( e2 ) > minangle )
1826  {
1827  doOnlySwap( ( uint )e2 );
1828  swapped = true;
1829  }
1830  else if ( swapPossible( ( uint )e3 ) && mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[e3]->getDual()]->getNext()]->getPoint()]->z() != el3 && swapMinAngle( e3 ) > minangle )
1831  {
1832  doOnlySwap( ( uint )e3 );
1833  swapped = true;
1834  }
1835  control[e1] = true;
1836  control[e2] = true;
1837  control[e3] = true;
1838  continue;
1839  }
1840  else//no horizontal triangle, go to the next one
1841  {
1842  control[e1] = true;
1843  control[e2] = true;
1844  control[e3] = true;
1845  continue;
1846  }
1847  }
1848  if ( !swapped )
1849  {
1850  delete[] control;
1851  break;
1852  }
1853  delete[] control;
1854  }
1855 
1856  //QgsDebugMsg( QStringLiteral( "end of method" ) );
1857 }
1858 
1860 {
1861  //minimum angle
1862  const double mintol = 17;//refinement stops after the minimum angle reached this tolerance
1863 
1864  //data structures
1865  std::map<int, double> edge_angle;//search tree with the edge number as key
1866  std::multimap<double, int> angle_edge;//multimap (map with not unique keys) with angle as key
1867  QSet<int> dontexamine;//search tree containing the edges which do not have to be examined (because of numerical problems)
1868 
1869 
1870  //first, go through all the forced edges and subdivide if they are encroached by a point
1871  bool stop = false;//flag to ensure that the for-loop is repeated until no half edge is split any more
1872 
1873  while ( !stop )
1874  {
1875  stop = true;
1876  const int nhalfedges = mHalfEdge.count();
1877 
1878  for ( int i = 0; i < nhalfedges - 1; i++ )
1879  {
1880  const int next = mHalfEdge[i]->getNext();
1881  const int nextnext = mHalfEdge[next]->getNext();
1882 
1883  if ( mHalfEdge[next]->getPoint() != -1 && ( mHalfEdge[i]->getForced() || mHalfEdge[mHalfEdge[mHalfEdge[i]->getDual()]->getNext()]->getPoint() == -1 ) )//check for encroached points on forced segments and segments on the inner side of the convex hull, but don't consider edges on the outer side of the convex hull
1884  {
1885  if ( !( ( mHalfEdge[next]->getForced() || edgeOnConvexHull( next ) ) || ( mHalfEdge[nextnext]->getForced() || edgeOnConvexHull( nextnext ) ) ) ) //don't consider triangles where all three edges are forced edges or hull edges
1886  {
1887  //test for encroachment
1888  while ( MathUtils::inDiametral( mPointVector[mHalfEdge[mHalfEdge[i]->getDual()]->getPoint()], mPointVector[mHalfEdge[i]->getPoint()], mPointVector[mHalfEdge[next]->getPoint()] ) )
1889  {
1890  //split segment
1891  const int pointno = splitHalfEdge( i, 0.5 );
1892  Q_UNUSED( pointno )
1893  stop = false;
1894  }
1895  }
1896  }
1897  }
1898  }
1899 
1900  //examine the triangulation for angles below the minimum and insert the edges into angle_edge and edge_angle, except the small angle is between forced segments or convex hull edges
1901  double angle;//angle between edge i and the consecutive edge
1902  int p1, p2, p3;//numbers of the triangle points
1903  for ( int i = 0; i < mHalfEdge.count() - 1; i++ )
1904  {
1905  p1 = mHalfEdge[mHalfEdge[i]->getDual()]->getPoint();
1906  p2 = mHalfEdge[i]->getPoint();
1907  p3 = mHalfEdge[mHalfEdge[i]->getNext()]->getPoint();
1908 
1909  if ( p1 == -1 || p2 == -1 || p3 == -1 )//don't consider triangles with the virtual point
1910  {
1911  continue;
1912  }
1913  angle = MathUtils::angle( mPointVector[p1], mPointVector[p2], mPointVector[p3], mPointVector[p2] );
1914 
1915  bool twoforcededges;//flag to decide, if edges should be added to the maps. Do not add them if true
1916 
1917 
1918  twoforcededges = ( mHalfEdge[i]->getForced() || edgeOnConvexHull( i ) ) && ( mHalfEdge[mHalfEdge[i]->getNext()]->getForced() || edgeOnConvexHull( mHalfEdge[i]->getNext() ) );
1919 
1920  if ( angle < mintol && !twoforcededges )
1921  {
1922  edge_angle.insert( std::make_pair( i, angle ) );
1923  angle_edge.insert( std::make_pair( angle, i ) );
1924  }
1925  }
1926 
1927  //debugging: print out all the angles below the minimum for a test
1928  for ( std::multimap<double, int>::const_iterator it = angle_edge.begin(); it != angle_edge.end(); ++it )
1929  {
1930  QgsDebugMsg( QStringLiteral( "angle: %1" ).arg( it->first ) );
1931  }
1932 
1933 
1934  double minangle = 0;//actual minimum angle
1935  int minedge;//first edge adjacent to the minimum angle
1936  int minedgenext;
1937  int minedgenextnext;
1938 
1939  QgsPoint circumcenter( 0, 0, 0 );
1940 
1941  while ( !edge_angle.empty() )
1942  {
1943  minangle = angle_edge.begin()->first;
1944  QgsDebugMsg( QStringLiteral( "minangle: %1" ).arg( minangle ) );
1945  minedge = angle_edge.begin()->second;
1946  minedgenext = mHalfEdge[minedge]->getNext();
1947  minedgenextnext = mHalfEdge[minedgenext]->getNext();
1948 
1949  //calculate the circumcenter
1950  if ( !MathUtils::circumcenter( mPointVector[mHalfEdge[minedge]->getPoint()], mPointVector[mHalfEdge[minedgenext]->getPoint()], mPointVector[mHalfEdge[minedgenextnext]->getPoint()], &circumcenter ) )
1951  {
1952  QgsDebugMsg( QStringLiteral( "warning, calculation of circumcenter failed" ) );
1953  //put all three edges to dontexamine and remove them from the other maps
1954  dontexamine.insert( minedge );
1955  edge_angle.erase( minedge );
1956  std::multimap<double, int>::iterator minedgeiter = angle_edge.find( minangle );
1957  while ( minedgeiter->second != minedge )
1958  {
1959  ++minedgeiter;
1960  }
1961  angle_edge.erase( minedgeiter );
1962  continue;
1963  }
1964 
1965  if ( !pointInside( circumcenter.x(), circumcenter.y() ) )
1966  {
1967  //put all three edges to dontexamine and remove them from the other maps
1968  QgsDebugMsg( QStringLiteral( "put circumcenter %1//%2 on dontexamine list because it is outside the convex hull" )
1969  .arg( circumcenter.x() ).arg( circumcenter.y() ) );
1970  dontexamine.insert( minedge );
1971  edge_angle.erase( minedge );
1972  std::multimap<double, int>::iterator minedgeiter = angle_edge.find( minangle );
1973  while ( minedgeiter->second != minedge )
1974  {
1975  ++minedgeiter;
1976  }
1977  angle_edge.erase( minedgeiter );
1978  continue;
1979  }
1980 
1981  /*******find out, if any forced segment or convex hull segment is in the influence region of the circumcenter. In this case, split the segment********************/
1982 
1983  bool encroached = false;
1984 
1985 #if 0 //slow version
1986  int numhalfedges = mHalfEdge.count();//begin slow version
1987  for ( int i = 0; i < numhalfedges; i++ )
1988  {
1989  if ( mHalfEdge[i]->getForced() || edgeOnConvexHull( i ) )
1990  {
1991  if ( MathUtils::inDiametral( mPointVector[mHalfEdge[i]->getPoint()], mPointVector[mHalfEdge[mHalfEdge[i]->getDual()]->getPoint()], &circumcenter ) )
1992  {
1993  encroached = true;
1994  //split segment
1995  QgsDebugMsg( QStringLiteral( "segment split" ) );
1996  int pointno = splitHalfEdge( i, 0.5 );
1997 
1998  //update dontexmine, angle_edge, edge_angle
1999  int pointingedge = baseEdgeOfPoint( pointno );
2000 
2001  int actedge = pointingedge;
2002  int ed1, ed2, ed3;//numbers of the three edges
2003  int pt1, pt2, pt3;//numbers of the three points
2004  double angle1, angle2, angle3;
2005 
2006  do
2007  {
2008  ed1 = mHalfEdge[actedge]->getDual();
2009  pt1 = mHalfEdge[ed1]->getPoint();
2010  ed2 = mHalfEdge[ed1]->getNext();
2011  pt2 = mHalfEdge[ed2]->getPoint();
2012  ed3 = mHalfEdge[ed2]->getNext();
2013  pt3 = mHalfEdge[ed3]->getPoint();
2014  actedge = ed3;
2015 
2016  if ( pt1 == -1 || pt2 == -1 || pt3 == -1 )//don't consider triangles with the virtual point
2017  {
2018  continue;
2019  }
2020 
2021  angle1 = MathUtils::angle( mPointVector[pt3], mPointVector[pt1], mPointVector[pt2], mPointVector[pt1] );
2022  angle2 = MathUtils::angle( mPointVector[pt1], mPointVector[pt2], mPointVector[pt3], mPointVector[pt2] );
2023  angle3 = MathUtils::angle( mPointVector[pt2], mPointVector[pt3], mPointVector[pt1], mPointVector[pt3] );
2024 
2025  //don't put the edges on the maps if two segments are forced or on a hull
2026  bool twoforcededges1, twoforcededges2, twoforcededges3;//flag to indicate, if angle1, angle2 and angle3 are between forced edges or hull edges
2027 
2028  if ( ( mHalfEdge[ed1]->getForced() || edgeOnConvexHull( ed1 ) ) && ( mHalfEdge[ed2]->getForced() || edgeOnConvexHull( ed2 ) ) )
2029  {
2030  twoforcededges1 = true;
2031  }
2032  else
2033  {
2034  twoforcededges1 = false;
2035  }
2036 
2037  if ( ( mHalfEdge[ed2]->getForced() || edgeOnConvexHull( ed2 ) ) && ( mHalfEdge[ed3]->getForced() || edgeOnConvexHull( ed3 ) ) )
2038  {
2039  twoforcededges2 = true;
2040  }
2041  else
2042  {
2043  twoforcededges2 = false;
2044  }
2045 
2046  if ( ( mHalfEdge[ed3]->getForced() || edgeOnConvexHull( ed3 ) ) && ( mHalfEdge[ed1]->getForced() || edgeOnConvexHull( ed1 ) ) )
2047  {
2048  twoforcededges3 = true;
2049  }
2050  else
2051  {
2052  twoforcededges3 = false;
2053  }
2054 
2055  //update the settings related to ed1
2056  QSet<int>::iterator ed1iter = dontexamine.find( ed1 );
2057  if ( ed1iter != dontexamine.end() )
2058  {
2059  //edge number is on dontexamine list
2060  dontexamine.erase( ed1iter );
2061  }
2062  else
2063  {
2064  //test, if it is on edge_angle and angle_edge and erase them if yes
2065  std::map<int, double>::iterator tempit1;
2066  tempit1 = edge_angle.find( ed1 );
2067  if ( tempit1 != edge_angle.end() )
2068  {
2069  //erase the entries
2070  double angle = tempit1->second;
2071  edge_angle.erase( ed1 );
2072  std::multimap<double, int>::iterator tempit2 = angle_edge.lower_bound( angle );
2073  while ( tempit2->second != ed1 )
2074  {
2075  ++tempit2;
2076  }
2077  angle_edge.erase( tempit2 );
2078  }
2079  }
2080 
2081  if ( angle1 < mintol && !twoforcededges1 )
2082  {
2083  edge_angle.insert( std::make_pair( ed1, angle1 ) );
2084  angle_edge.insert( std::make_pair( angle1, ed1 ) );
2085  }
2086 
2087  //update the settings related to ed2
2088  QSet<int>::iterator ed2iter = dontexamine.find( ed2 );
2089  if ( ed2iter != dontexamine.end() )
2090  {
2091  //edge number is on dontexamine list
2092  dontexamine.erase( ed2iter );
2093  }
2094  else
2095  {
2096  //test, if it is on edge_angle and angle_edge and erase them if yes
2097  std::map<int, double>::iterator tempit1;
2098  tempit1 = edge_angle.find( ed2 );
2099  if ( tempit1 != edge_angle.end() )
2100  {
2101  //erase the entries
2102  double angle = tempit1->second;
2103  edge_angle.erase( ed2 );
2104  std::multimap<double, int>::iterator tempit2 = angle_edge.lower_bound( angle );
2105  while ( tempit2->second != ed2 )
2106  {
2107  ++tempit2;
2108  }
2109  angle_edge.erase( tempit2 );
2110  }
2111  }
2112 
2113  if ( angle2 < mintol && !twoforcededges2 )
2114  {
2115  edge_angle.insert( std::make_pair( ed2, angle2 ) );
2116  angle_edge.insert( std::make_pair( angle2, ed2 ) );
2117  }
2118 
2119  //update the settings related to ed3
2120  QSet<int>::iterator ed3iter = dontexamine.find( ed3 );
2121  if ( ed3iter != dontexamine.end() )
2122  {
2123  //edge number is on dontexamine list
2124  dontexamine.erase( ed3iter );
2125  }
2126  else
2127  {
2128  //test, if it is on edge_angle and angle_edge and erase them if yes
2129  std::map<int, double>::iterator tempit1;
2130  tempit1 = edge_angle.find( ed3 );
2131  if ( tempit1 != edge_angle.end() )
2132  {
2133  //erase the entries
2134  double angle = tempit1->second;
2135  edge_angle.erase( ed3 );
2136  std::multimap<double, int>::iterator tempit2 = angle_edge.lower_bound( angle );
2137  while ( tempit2->second != ed3 )
2138  {
2139  ++tempit2;
2140  }
2141  angle_edge.erase( tempit2 );
2142  }
2143  }
2144 
2145  if ( angle3 < mintol && !twoforcededges3 )
2146  {
2147  edge_angle.insert( std::make_pair( ed3, angle3 ) );
2148  angle_edge.insert( std::make_pair( angle3, ed3 ) );
2149  }
2150 
2151 
2152  }
2153  while ( actedge != pointingedge );
2154  }
2155  }
2156  }
2157 #endif //end slow version
2158 
2159 
2160  //fast version. Maybe this does not work
2161  QSet<int> influenceedges;//begin fast method
2162  int baseedge = baseEdgeOfTriangle( circumcenter );
2163  if ( baseedge == -5 )//a numerical instability occurred or the circumcenter already exists in the triangulation
2164  {
2165  //delete minedge from edge_angle and minangle from angle_edge
2166  edge_angle.erase( minedge );
2167  std::multimap<double, int>::iterator minedgeiter = angle_edge.find( minangle );
2168  while ( minedgeiter->second != minedge )
2169  {
2170  ++minedgeiter;
2171  }
2172  angle_edge.erase( minedgeiter );
2173  continue;
2174  }
2175  else if ( baseedge == -25 )//circumcenter already exists in the triangulation
2176  {
2177  //delete minedge from edge_angle and minangle from angle_edge
2178  edge_angle.erase( minedge );
2179  std::multimap<double, int>::iterator minedgeiter = angle_edge.find( minangle );
2180  while ( minedgeiter->second != minedge )
2181  {
2182  ++minedgeiter;
2183  }
2184  angle_edge.erase( minedgeiter );
2185  continue;
2186  }
2187  else if ( baseedge == -20 )
2188  {
2189  baseedge = mEdgeWithPoint;
2190  }
2191 
2192  evaluateInfluenceRegion( &circumcenter, baseedge, influenceedges );
2193  evaluateInfluenceRegion( &circumcenter, mHalfEdge[baseedge]->getNext(), influenceedges );
2194  evaluateInfluenceRegion( &circumcenter, mHalfEdge[mHalfEdge[baseedge]->getNext()]->getNext(), influenceedges );
2195 
2196  for ( QSet<int>::iterator it = influenceedges.begin(); it != influenceedges.end(); ++it )
2197  {
2198  if ( ( mHalfEdge[*it]->getForced() || edgeOnConvexHull( *it ) ) && MathUtils::inDiametral( mPointVector[mHalfEdge[*it]->getPoint()], mPointVector[mHalfEdge[mHalfEdge[*it]->getDual()]->getPoint()], &circumcenter ) )
2199  {
2200  //split segment
2201  QgsDebugMsg( QStringLiteral( "segment split" ) );
2202  const int pointno = splitHalfEdge( *it, 0.5 );
2203  encroached = true;
2204 
2205  //update dontexmine, angle_edge, edge_angle
2206  const int pointingedge = baseEdgeOfPoint( pointno );
2207 
2208  int actedge = pointingedge;
2209  int ed1, ed2, ed3;//numbers of the three edges
2210  int pt1, pt2, pt3;//numbers of the three points
2211  double angle1, angle2, angle3;
2212 
2213  do
2214  {
2215  ed1 = mHalfEdge[actedge]->getDual();
2216  pt1 = mHalfEdge[ed1]->getPoint();
2217  ed2 = mHalfEdge[ed1]->getNext();
2218  pt2 = mHalfEdge[ed2]->getPoint();
2219  ed3 = mHalfEdge[ed2]->getNext();
2220  pt3 = mHalfEdge[ed3]->getPoint();
2221  actedge = ed3;
2222 
2223  if ( pt1 == -1 || pt2 == -1 || pt3 == -1 )//don't consider triangles with the virtual point
2224  {
2225  continue;
2226  }
2227 
2228  angle1 = MathUtils::angle( mPointVector[pt3], mPointVector[pt1], mPointVector[pt2], mPointVector[pt1] );
2229  angle2 = MathUtils::angle( mPointVector[pt1], mPointVector[pt2], mPointVector[pt3], mPointVector[pt2] );
2230  angle3 = MathUtils::angle( mPointVector[pt2], mPointVector[pt3], mPointVector[pt1], mPointVector[pt3] );
2231 
2232  //don't put the edges on the maps if two segments are forced or on a hull
2233  bool twoforcededges1, twoforcededges2, twoforcededges3;//flag to decide, if edges should be added to the maps. Do not add them if true
2234 
2235 
2236 
2237  twoforcededges1 = ( mHalfEdge[ed1]->getForced() || edgeOnConvexHull( ed1 ) ) && ( mHalfEdge[ed2]->getForced() || edgeOnConvexHull( ed2 ) );
2238 
2239  twoforcededges2 = ( mHalfEdge[ed2]->getForced() || edgeOnConvexHull( ed2 ) ) && ( mHalfEdge[ed3]->getForced() || edgeOnConvexHull( ed3 ) );
2240 
2241  twoforcededges3 = ( mHalfEdge[ed3]->getForced() || edgeOnConvexHull( ed3 ) ) && ( mHalfEdge[ed1]->getForced() || edgeOnConvexHull( ed1 ) );
2242 
2243 
2244  //update the settings related to ed1
2245  const QSet<int>::iterator ed1iter = dontexamine.find( ed1 );
2246  if ( ed1iter != dontexamine.end() )
2247  {
2248  //edge number is on dontexamine list
2249  dontexamine.erase( ed1iter );
2250  }
2251  else
2252  {
2253  //test, if it is on edge_angle and angle_edge and erase them if yes
2254  std::map<int, double>::iterator tempit1;
2255  tempit1 = edge_angle.find( ed1 );
2256  if ( tempit1 != edge_angle.end() )
2257  {
2258  //erase the entries
2259  const double angle = tempit1->second;
2260  edge_angle.erase( ed1 );
2261  std::multimap<double, int>::iterator tempit2 = angle_edge.lower_bound( angle );
2262  while ( tempit2->second != ed1 )
2263  {
2264  ++tempit2;
2265  }
2266  angle_edge.erase( tempit2 );
2267  }
2268  }
2269 
2270  if ( angle1 < mintol && !twoforcededges1 )
2271  {
2272  edge_angle.insert( std::make_pair( ed1, angle1 ) );
2273  angle_edge.insert( std::make_pair( angle1, ed1 ) );
2274  }
2275 
2276  //update the settings related to ed2
2277  const QSet<int>::iterator ed2iter = dontexamine.find( ed2 );
2278  if ( ed2iter != dontexamine.end() )
2279  {
2280  //edge number is on dontexamine list
2281  dontexamine.erase( ed2iter );
2282  }
2283  else
2284  {
2285  //test, if it is on edge_angle and angle_edge and erase them if yes
2286  std::map<int, double>::iterator tempit1;
2287  tempit1 = edge_angle.find( ed2 );
2288  if ( tempit1 != edge_angle.end() )
2289  {
2290  //erase the entries
2291  const double angle = tempit1->second;
2292  edge_angle.erase( ed2 );
2293  std::multimap<double, int>::iterator tempit2 = angle_edge.lower_bound( angle );
2294  while ( tempit2->second != ed2 )
2295  {
2296  ++tempit2;
2297  }
2298  angle_edge.erase( tempit2 );
2299  }
2300  }
2301 
2302  if ( angle2 < mintol && !twoforcededges2 )
2303  {
2304  edge_angle.insert( std::make_pair( ed2, angle2 ) );
2305  angle_edge.insert( std::make_pair( angle2, ed2 ) );
2306  }
2307 
2308  //update the settings related to ed3
2309  const QSet<int>::iterator ed3iter = dontexamine.find( ed3 );
2310  if ( ed3iter != dontexamine.end() )
2311  {
2312  //edge number is on dontexamine list
2313  dontexamine.erase( ed3iter );
2314  }
2315  else
2316  {
2317  //test, if it is on edge_angle and angle_edge and erase them if yes
2318  std::map<int, double>::iterator tempit1;
2319  tempit1 = edge_angle.find( ed3 );
2320  if ( tempit1 != edge_angle.end() )
2321  {
2322  //erase the entries
2323  const double angle = tempit1->second;
2324  edge_angle.erase( ed3 );
2325  std::multimap<double, int>::iterator tempit2 = angle_edge.lower_bound( angle );
2326  while ( tempit2->second != ed3 )
2327  {
2328  ++tempit2;
2329  }
2330  angle_edge.erase( tempit2 );
2331  }
2332  }
2333 
2334  if ( angle3 < mintol && !twoforcededges3 )
2335  {
2336  edge_angle.insert( std::make_pair( ed3, angle3 ) );
2337  angle_edge.insert( std::make_pair( angle3, ed3 ) );
2338  }
2339 
2340 
2341  }
2342  while ( actedge != pointingedge );
2343  }
2344  } //end fast method
2345 
2346 
2347  if ( encroached )
2348  {
2349  continue;
2350  }
2351 
2352  /*******otherwise, try to add the circumcenter to the triangulation************************************************************************************************/
2353 
2354  QgsPoint p( 0, 0, 0 );
2355  calcPoint( circumcenter.x(), circumcenter.y(), p );
2356  const int pointno = addPoint( p );
2357 
2358  if ( pointno == -100 || pointno == mTwiceInsPoint )
2359  {
2360  if ( pointno == -100 )
2361  {
2362  QgsDebugMsg( QStringLiteral( "put circumcenter %1//%2 on dontexamine list because of numerical instabilities" )
2363  .arg( circumcenter.x() ).arg( circumcenter.y() ) );
2364  }
2365  else if ( pointno == mTwiceInsPoint )
2366  {
2367  QgsDebugMsg( QStringLiteral( "put circumcenter %1//%2 on dontexamine list because it is already inserted" )
2368  .arg( circumcenter.x() ).arg( circumcenter.y() ) );
2369  //test, if the point is present in the triangulation
2370  bool flag = false;
2371  for ( int i = 0; i < mPointVector.count(); i++ )
2372  {
2373  if ( mPointVector[i]->x() == circumcenter.x() && mPointVector[i]->y() == circumcenter.y() )
2374  {
2375  flag = true;
2376  }
2377  }
2378  if ( !flag )
2379  {
2380  QgsDebugMsg( QStringLiteral( "point is not present in the triangulation" ) );
2381  }
2382  }
2383  //put all three edges to dontexamine and remove them from the other maps
2384  dontexamine.insert( minedge );
2385  edge_angle.erase( minedge );
2386  std::multimap<double, int>::iterator minedgeiter = angle_edge.lower_bound( minangle );
2387  while ( minedgeiter->second != minedge )
2388  {
2389  ++minedgeiter;
2390  }
2391  angle_edge.erase( minedgeiter );
2392  continue;
2393  }
2394  else//insertion successful
2395  {
2396  QgsDebugMsg( QStringLiteral( "circumcenter added" ) );
2397 
2398  //update the maps
2399  //go around the inserted point and make changes for every half edge
2400  const int pointingedge = baseEdgeOfPoint( pointno );
2401 
2402  int actedge = pointingedge;
2403  int ed1, ed2, ed3;//numbers of the three edges
2404  int pt1, pt2, pt3;//numbers of the three points
2405  double angle1, angle2, angle3;
2406 
2407  do
2408  {
2409  ed1 = mHalfEdge[actedge]->getDual();
2410  pt1 = mHalfEdge[ed1]->getPoint();
2411  ed2 = mHalfEdge[ed1]->getNext();
2412  pt2 = mHalfEdge[ed2]->getPoint();
2413  ed3 = mHalfEdge[ed2]->getNext();
2414  pt3 = mHalfEdge[ed3]->getPoint();
2415  actedge = ed3;
2416 
2417  if ( pt1 == -1 || pt2 == -1 || pt3 == -1 )//don't consider triangles with the virtual point
2418  {
2419  continue;
2420  }
2421 
2422  angle1 = MathUtils::angle( mPointVector[pt3], mPointVector[pt1], mPointVector[pt2], mPointVector[pt1] );
2423  angle2 = MathUtils::angle( mPointVector[pt1], mPointVector[pt2], mPointVector[pt3], mPointVector[pt2] );
2424  angle3 = MathUtils::angle( mPointVector[pt2], mPointVector[pt3], mPointVector[pt1], mPointVector[pt3] );
2425 
2426  //todo: put all three edges on the dontexamine list if two edges are forced or convex hull edges
2427  bool twoforcededges1, twoforcededges2, twoforcededges3;
2428 
2429  twoforcededges1 = ( mHalfEdge[ed1]->getForced() || edgeOnConvexHull( ed1 ) ) && ( mHalfEdge[ed2]->getForced() || edgeOnConvexHull( ed2 ) );
2430 
2431  twoforcededges2 = ( mHalfEdge[ed2]->getForced() || edgeOnConvexHull( ed2 ) ) && ( mHalfEdge[ed3]->getForced() || edgeOnConvexHull( ed3 ) );
2432 
2433  twoforcededges3 = ( mHalfEdge[ed3]->getForced() || edgeOnConvexHull( ed3 ) ) && ( mHalfEdge[ed1]->getForced() || edgeOnConvexHull( ed1 ) );
2434 
2435 
2436  //update the settings related to ed1
2437  const QSet<int>::iterator ed1iter = dontexamine.find( ed1 );
2438  if ( ed1iter != dontexamine.end() )
2439  {
2440  //edge number is on dontexamine list
2441  dontexamine.erase( ed1iter );
2442  }
2443  else
2444  {
2445  //test, if it is on edge_angle and angle_edge and erase them if yes
2446  std::map<int, double>::iterator tempit1;
2447  tempit1 = edge_angle.find( ed1 );
2448  if ( tempit1 != edge_angle.end() )
2449  {
2450  //erase the entries
2451  const double angle = tempit1->second;
2452  edge_angle.erase( ed1 );
2453  std::multimap<double, int>::iterator tempit2 = angle_edge.lower_bound( angle );
2454  while ( tempit2->second != ed1 )
2455  {
2456  ++tempit2;
2457  }
2458  angle_edge.erase( tempit2 );
2459  }
2460  }
2461 
2462  if ( angle1 < mintol && !twoforcededges1 )
2463  {
2464  edge_angle.insert( std::make_pair( ed1, angle1 ) );
2465  angle_edge.insert( std::make_pair( angle1, ed1 ) );
2466  }
2467 
2468 
2469  //update the settings related to ed2
2470  const QSet<int>::iterator ed2iter = dontexamine.find( ed2 );
2471  if ( ed2iter != dontexamine.end() )
2472  {
2473  //edge number is on dontexamine list
2474  dontexamine.erase( ed2iter );
2475  }
2476  else
2477  {
2478  //test, if it is on edge_angle and angle_edge and erase them if yes
2479  std::map<int, double>::iterator tempit1;
2480  tempit1 = edge_angle.find( ed2 );
2481  if ( tempit1 != edge_angle.end() )
2482  {
2483  //erase the entries
2484  const double angle = tempit1->second;
2485  edge_angle.erase( ed2 );
2486  std::multimap<double, int>::iterator tempit2 = angle_edge.lower_bound( angle );
2487  while ( tempit2->second != ed2 )
2488  {
2489  ++tempit2;
2490  }
2491  angle_edge.erase( tempit2 );
2492  }
2493  }
2494 
2495  if ( angle2 < mintol && !twoforcededges2 )
2496  {
2497  edge_angle.insert( std::make_pair( ed2, angle2 ) );
2498  angle_edge.insert( std::make_pair( angle2, ed2 ) );
2499  }
2500 
2501  //update the settings related to ed3
2502  const QSet<int>::iterator ed3iter = dontexamine.find( ed3 );
2503  if ( ed3iter != dontexamine.end() )
2504  {
2505  //edge number is on dontexamine list
2506  dontexamine.erase( ed3iter );
2507  }
2508  else
2509  {
2510  //test, if it is on edge_angle and angle_edge and erase them if yes
2511  std::map<int, double>::iterator tempit1;
2512  tempit1 = edge_angle.find( ed3 );
2513  if ( tempit1 != edge_angle.end() )
2514  {
2515  //erase the entries
2516  const double angle = tempit1->second;
2517  edge_angle.erase( ed3 );
2518  std::multimap<double, int>::iterator tempit2 = angle_edge.lower_bound( angle );
2519  while ( tempit2->second != ed3 )
2520  {
2521  ++tempit2;
2522  }
2523  angle_edge.erase( tempit2 );
2524  }
2525  }
2526 
2527  if ( angle3 < mintol && !twoforcededges3 )
2528  {
2529  edge_angle.insert( std::make_pair( ed3, angle3 ) );
2530  angle_edge.insert( std::make_pair( angle3, ed3 ) );
2531  }
2532 
2533 
2534  }
2535  while ( actedge != pointingedge );
2536  }
2537  }
2538 
2539 #if 0
2540  //debugging: print out all edge of dontexamine
2541  for ( QSet<int>::iterator it = dontexamine.begin(); it != dontexamine.end(); ++it )
2542  {
2543  QgsDebugMsg( QStringLiteral( "edge nr. %1 is in dontexamine" ).arg( *it ) );
2544  }
2545 #endif
2546 }
2547 
2548 
2549 bool QgsDualEdgeTriangulation::swapPossible( unsigned int edge ) const
2550 {
2551  //test, if edge belongs to a forced edge
2552  if ( mHalfEdge[edge]->getForced() )
2553  {
2554  return false;
2555  }
2556 
2557  //test, if the edge is on the convex hull or is connected to the virtual point
2558  if ( mHalfEdge[edge]->getPoint() == -1 || mHalfEdge[mHalfEdge[edge]->getNext()]->getPoint() == -1 || mHalfEdge[mHalfEdge[mHalfEdge[edge]->getDual()]->getNext()]->getPoint() == -1 || mHalfEdge[mHalfEdge[edge]->getDual()]->getPoint() == -1 )
2559  {
2560  return false;
2561  }
2562  //then, test, if the edge is in the middle of a not convex quad
2563  QgsPoint *pta = mPointVector[mHalfEdge[edge]->getPoint()];
2564  QgsPoint *ptb = mPointVector[mHalfEdge[mHalfEdge[edge]->getNext()]->getPoint()];
2565  QgsPoint *ptc = mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[edge]->getNext()]->getNext()]->getPoint()];
2566  QgsPoint *ptd = mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[edge]->getDual()]->getNext()]->getPoint()];
2567  if ( MathUtils::leftOf( *ptc, pta, ptb ) > leftOfTresh )
2568  {
2569  return false;
2570  }
2571  else if ( MathUtils::leftOf( *ptd, ptb, ptc ) > leftOfTresh )
2572  {
2573  return false;
2574  }
2575  else if ( MathUtils::leftOf( *pta, ptc, ptd ) > leftOfTresh )
2576  {
2577  return false;
2578  }
2579  else if ( MathUtils::leftOf( *ptb, ptd, pta ) > leftOfTresh )
2580  {
2581  return false;
2582  }
2583  return true;
2584 }
2585 
2586 void QgsDualEdgeTriangulation::triangulatePolygon( QList<int> *poly, QList<int> *free, int mainedge )
2587 {
2588  if ( poly && free )
2589  {
2590  if ( poly->count() == 3 )//polygon is already a triangle
2591  {
2592  return;
2593  }
2594 
2595  //search for the edge pointing on the closest point(distedge) and for the next(nextdistedge)
2596  QList<int>::const_iterator iterator = ++( poly->constBegin() );//go to the second edge
2597  double distance = MathUtils::distPointFromLine( mPointVector[mHalfEdge[( *iterator )]->getPoint()], mPointVector[mHalfEdge[mHalfEdge[mainedge]->getDual()]->getPoint()], mPointVector[mHalfEdge[mainedge]->getPoint()] );
2598  int distedge = ( *iterator );
2599  int nextdistedge = mHalfEdge[( *iterator )]->getNext();
2600  ++iterator;
2601 
2602  while ( iterator != --( poly->constEnd() ) )
2603  {
2604  if ( MathUtils::distPointFromLine( mPointVector[mHalfEdge[( *iterator )]->getPoint()], mPointVector[mHalfEdge[mHalfEdge[mainedge]->getDual()]->getPoint()], mPointVector[mHalfEdge[mainedge]->getPoint()] ) < distance )
2605  {
2606  distedge = ( *iterator );
2607  nextdistedge = mHalfEdge[( *iterator )]->getNext();
2608  distance = MathUtils::distPointFromLine( mPointVector[mHalfEdge[( *iterator )]->getPoint()], mPointVector[mHalfEdge[mHalfEdge[mainedge]->getDual()]->getPoint()], mPointVector[mHalfEdge[mainedge]->getPoint()] );
2609  }
2610  ++iterator;
2611  }
2612 
2613  if ( nextdistedge == ( *( --poly->end() ) ) )//the nearest point is connected to the endpoint of mainedge
2614  {
2615  const int inserta = free->first();//take an edge from the freelist
2616  const int insertb = mHalfEdge[inserta]->getDual();
2617  free->pop_front();
2618 
2619  mHalfEdge[inserta]->setNext( ( poly->at( 1 ) ) );
2620  mHalfEdge[inserta]->setPoint( mHalfEdge[mainedge]->getPoint() );
2621  mHalfEdge[insertb]->setNext( nextdistedge );
2622  mHalfEdge[insertb]->setPoint( mHalfEdge[distedge]->getPoint() );
2623  mHalfEdge[distedge]->setNext( inserta );
2624  mHalfEdge[mainedge]->setNext( insertb );
2625 
2626  QList<int> polya;
2627  for ( iterator = ( ++( poly->constBegin() ) ); ( *iterator ) != nextdistedge; ++iterator )
2628  {
2629  polya.append( ( *iterator ) );
2630  }
2631  polya.prepend( inserta );
2632 
2633 #if 0
2634  //print out all the elements of polya for a test
2635  for ( iterator = polya.begin(); iterator != polya.end(); ++iterator )
2636  {
2637  QgsDebugMsg( *iterator );
2638  }
2639 #endif
2640 
2641  triangulatePolygon( &polya, free, inserta );
2642  }
2643 
2644  else if ( distedge == ( *( ++poly->begin() ) ) )//the nearest point is connected to the beginpoint of mainedge
2645  {
2646  const int inserta = free->first();//take an edge from the freelist
2647  const int insertb = mHalfEdge[inserta]->getDual();
2648  free->pop_front();
2649 
2650  mHalfEdge[inserta]->setNext( ( poly->at( 2 ) ) );
2651  mHalfEdge[inserta]->setPoint( mHalfEdge[distedge]->getPoint() );
2652  mHalfEdge[insertb]->setNext( mainedge );
2653  mHalfEdge[insertb]->setPoint( mHalfEdge[mHalfEdge[mainedge]->getDual()]->getPoint() );
2654  mHalfEdge[distedge]->setNext( insertb );
2655  mHalfEdge[( *( --poly->end() ) )]->setNext( inserta );
2656 
2657  QList<int> polya;
2658  iterator = poly->constBegin();
2659  iterator += 2;
2660  while ( iterator != poly->constEnd() )
2661  {
2662  polya.append( ( *iterator ) );
2663  ++iterator;
2664  }
2665  polya.prepend( inserta );
2666 
2667  triangulatePolygon( &polya, free, inserta );
2668  }
2669 
2670  else//the nearest point is not connected to an endpoint of mainedge
2671  {
2672  const int inserta = free->first();//take an edge from the freelist
2673  const int insertb = mHalfEdge[inserta]->getDual();
2674  free->pop_front();
2675 
2676  const int insertc = free->first();
2677  const int insertd = mHalfEdge[insertc]->getDual();
2678  free->pop_front();
2679 
2680  mHalfEdge[inserta]->setNext( ( poly->at( 1 ) ) );
2681  mHalfEdge[inserta]->setPoint( mHalfEdge[mainedge]->getPoint() );
2682  mHalfEdge[insertb]->setNext( insertd );
2683  mHalfEdge[insertb]->setPoint( mHalfEdge[distedge]->getPoint() );
2684  mHalfEdge[insertc]->setNext( nextdistedge );
2685  mHalfEdge[insertc]->setPoint( mHalfEdge[distedge]->getPoint() );
2686  mHalfEdge[insertd]->setNext( mainedge );
2687  mHalfEdge[insertd]->setPoint( mHalfEdge[mHalfEdge[mainedge]->getDual()]->getPoint() );
2688 
2689  mHalfEdge[distedge]->setNext( inserta );
2690  mHalfEdge[mainedge]->setNext( insertb );
2691  mHalfEdge[( *( --poly->end() ) )]->setNext( insertc );
2692 
2693  //build two new polygons for recursive triangulation
2694  QList<int> polya;
2695  QList<int> polyb;
2696 
2697  for ( iterator = ++( poly->constBegin() ); ( *iterator ) != nextdistedge; ++iterator )
2698  {
2699  polya.append( ( *iterator ) );
2700  }
2701  polya.prepend( inserta );
2702 
2703 
2704  while ( iterator != poly->constEnd() )
2705  {
2706  polyb.append( ( *iterator ) );
2707  ++iterator;
2708  }
2709  polyb.prepend( insertc );
2710 
2711  triangulatePolygon( &polya, free, inserta );
2712  triangulatePolygon( &polyb, free, insertc );
2713  }
2714  }
2715 
2716  else
2717  {
2718  QgsDebugMsg( QStringLiteral( "warning, null pointer" ) );
2719  }
2720 
2721 }
2722 
2723 bool QgsDualEdgeTriangulation::pointInside( double x, double y )
2724 {
2725  const QgsPoint point( x, y, 0 );
2726  unsigned int actedge = mEdgeInside;//start with an edge which does not point to the virtual point
2727  int counter = 0;//number of consecutive successful left-of-tests
2728  int nulls = 0;//number of left-of-tests, which returned 0. 1 means, that the point is on a line, 2 means that it is on an existing point
2729  int numinstabs = 0;//number of suspect left-of-tests due to 'leftOfTresh'
2730  int runs = 0;//counter for the number of iterations in the loop to prevent an endless loop
2731 
2732  while ( true )
2733  {
2734  if ( runs > MAX_BASE_ITERATIONS )//prevents endless loops
2735  {
2736  QgsDebugMsg( QStringLiteral( "warning, instability detected: Point coordinates: %1//%2" ).arg( x ).arg( y ) );
2737  return false;
2738  }
2739 
2740  if ( MathUtils::leftOf( point, mPointVector[mHalfEdge[mHalfEdge[actedge]->getDual()]->getPoint()], mPointVector[mHalfEdge[actedge]->getPoint()] ) < ( -leftOfTresh ) )//point is on the left side
2741  {
2742  counter += 1;
2743  if ( counter == 3 )//three successful passes means that we have found the triangle
2744  {
2745  break;
2746  }
2747  }
2748 
2749  else if ( fabs( MathUtils::leftOf( point, mPointVector[mHalfEdge[mHalfEdge[actedge]->getDual()]->getPoint()], mPointVector[mHalfEdge[actedge]->getPoint()] ) ) <= leftOfTresh ) //point is exactly in the line of the edge
2750  {
2751  counter += 1;
2752  mEdgeWithPoint = actedge;
2753  nulls += 1;
2754  if ( counter == 3 )//three successful passes means that we have found the triangle
2755  {
2756  break;
2757  }
2758  }
2759 
2760  else//point is on the right side
2761  {
2762  actedge = mHalfEdge[actedge]->getDual();
2763  counter = 1;
2764  nulls = 0;
2765  numinstabs = 0;
2766  }
2767 
2768  actedge = mHalfEdge[actedge]->getNext();
2769  if ( mHalfEdge[actedge]->getPoint() == -1 )//the half edge points to the virtual point
2770  {
2771  if ( nulls == 1 )//point is exactly on the convex hull
2772  {
2773  return true;
2774  }
2775  mEdgeOutside = ( unsigned int )mHalfEdge[mHalfEdge[actedge]->getNext()]->getNext();
2776  return false;//the point is outside the convex hull
2777  }
2778  runs++;
2779  }
2780 
2781  if ( nulls == 2 )//we hit an existing point
2782  {
2783  return true;
2784  }
2785  if ( numinstabs > 0 )//a numerical instability occurred
2786  {
2787  QgsDebugMsg( QStringLiteral( "numerical instabilities" ) );
2788  return true;
2789  }
2790 
2791  if ( nulls == 1 )//point is on an existing edge
2792  {
2793  return true;
2794  }
2795  mEdgeInside = actedge;
2796  return true;
2797 }
2798 
2799 #if 0
2800 bool DualEdgeTriangulation::readFromTAFF( QString filename )
2801 {
2802  QApplication::setOverrideCursor( QCursor( Qt::WaitCursor ) );//change the cursor
2803 
2804  QFile file( filename );//the file to be read
2805  file.open( IO_Raw | IO_ReadOnly );
2806  QBuffer buffer( file.readAll() );//the buffer to copy the file to
2807  file.close();
2808 
2809  QTextStream textstream( &buffer );
2810  buffer.open( IO_ReadOnly );
2811 
2812  QString buff;
2813  int numberofhalfedges;
2814  int numberofpoints;
2815 
2816  //edge section
2817  while ( buff.mid( 0, 4 ) != "TRIA" )//search for the TRIA-section
2818  {
2819  buff = textstream.readLine();
2820  }
2821  while ( buff.mid( 0, 4 ) != "NEDG" )
2822  {
2823  buff = textstream.readLine();
2824  }
2825  numberofhalfedges = buff.section( ' ', 1, 1 ).toInt();
2826  mHalfEdge.resize( numberofhalfedges );
2827 
2828 
2829  while ( buff.mid( 0, 4 ) != "DATA" )//search for the data section
2830  {
2831  textstream >> buff;
2832  }
2833 
2834  int nr1, nr2, dual1, dual2, point1, point2, next1, next2, fo1, fo2, b1, b2;
2835  bool forced1, forced2, break1, break2;
2836 
2837 
2838  //import all the DualEdges
2839  QProgressBar *edgebar = new QProgressBar();//use a progress dialog so that it is not boring for the user
2840  edgebar->setCaption( tr( "Reading edges…" ) );
2841  edgebar->setTotalSteps( numberofhalfedges / 2 );
2842  edgebar->setMinimumWidth( 400 );
2843  edgebar->move( 500, 500 );
2844  edgebar->show();
2845 
2846  for ( int i = 0; i < numberofhalfedges / 2; i++ )
2847  {
2848  if ( i % 1000 == 0 )
2849  {
2850  edgebar->setProgress( i );
2851  }
2852 
2853  textstream >> nr1;
2854  textstream >> point1;
2855  textstream >> next1;
2856  textstream >> fo1;
2857 
2858  if ( fo1 == 0 )
2859  {
2860  forced1 = false;
2861  }
2862  else
2863  {
2864  forced1 = true;
2865  }
2866 
2867  textstream >> b1;
2868 
2869  if ( b1 == 0 )
2870  {
2871  break1 = false;
2872  }
2873  else
2874  {
2875  break1 = true;
2876  }
2877 
2878  textstream >> nr2;
2879  textstream >> point2;
2880  textstream >> next2;
2881  textstream >> fo2;
2882 
2883  if ( fo2 == 0 )
2884  {
2885  forced2 = false;
2886  }
2887  else
2888  {
2889  forced2 = true;
2890  }
2891 
2892  textstream >> b2;
2893 
2894  if ( b2 == 0 )
2895  {
2896  break2 = false;
2897  }
2898  else
2899  {
2900  break2 = true;
2901  }
2902 
2903  HalfEdge *hf1 = new HalfEdge();
2904  hf1->setDual( nr2 );
2905  hf1->setNext( next1 );
2906  hf1->setPoint( point1 );
2907  hf1->setBreak( break1 );
2908  hf1->setForced( forced1 );
2909 
2910  HalfEdge *hf2 = new HalfEdge();
2911  hf2->setDual( nr1 );
2912  hf2->setNext( next2 );
2913  hf2->setPoint( point2 );
2914  hf2->setBreak( break2 );
2915  hf2->setForced( forced2 );
2916 
2917  // QgsDebugMsg( QStringLiteral( "inserting half edge pair %1" ).arg( i ) );
2918  mHalfEdge.insert( nr1, hf1 );
2919  mHalfEdge.insert( nr2, hf2 );
2920 
2921  }
2922 
2923  edgebar->setProgress( numberofhalfedges / 2 );
2924  delete edgebar;
2925 
2926  //set mEdgeInside to a reasonable value
2927  for ( int i = 0; i < numberofhalfedges; i++ )
2928  {
2929  int a, b, c, d;
2930  a = mHalfEdge[i]->getPoint();
2931  b = mHalfEdge[mHalfEdge[i]->getDual()]->getPoint();
2932  c = mHalfEdge[mHalfEdge[i]->getNext()]->getPoint();
2933  d = mHalfEdge[mHalfEdge[mHalfEdge[i]->getDual()]->getNext()]->getPoint();
2934  if ( a != -1 && b != -1 && c != -1 && d != -1 )
2935  {
2936  mEdgeInside = i;
2937  break;
2938  }
2939  }
2940 
2941  //point section
2942  while ( buff.mid( 0, 4 ) != "POIN" )
2943  {
2944  buff = textstream.readLine();
2945  QgsDebugMsg( buff );
2946  }
2947  while ( buff.mid( 0, 4 ) != "NPTS" )
2948  {
2949  buff = textstream.readLine();
2950  QgsDebugMsg( buff );
2951  }
2952  numberofpoints = buff.section( ' ', 1, 1 ).toInt();
2953  mPointVector.resize( numberofpoints );
2954 
2955  while ( buff.mid( 0, 4 ) != "DATA" )
2956  {
2957  textstream >> buff;
2958  }
2959 
2960  QProgressBar *pointbar = new QProgressBar();
2961  pointbar->setCaption( tr( "Reading points…" ) );
2962  pointbar->setTotalSteps( numberofpoints );
2963  pointbar->setMinimumWidth( 400 );
2964  pointbar->move( 500, 500 );
2965  pointbar->show();
2966 
2967 
2968  double x, y, z;
2969  for ( int i = 0; i < numberofpoints; i++ )
2970  {
2971  if ( i % 1000 == 0 )
2972  {
2973  pointbar->setProgress( i );
2974  }
2975 
2976  textstream >> x;
2977  textstream >> y;
2978  textstream >> z;
2979 
2980  QgsPoint *p = new QgsPoint( x, y, z );
2981 
2982  // QgsDebugMsg( QStringLiteral( "inserting point %1" ).arg( i ) );
2983  mPointVector.insert( i, p );
2984 
2985  if ( i == 0 )
2986  {
2987  xMin = x;
2988  xMax = x;
2989  yMin = y;
2990  yMax = y;
2991  }
2992  else
2993  {
2994  //update the bounding box
2995  if ( x < xMin )
2996  {
2997  xMin = x;
2998  }
2999  else if ( x > xMax )
3000  {
3001  xMax = x;
3002  }
3003 
3004  if ( y < yMin )
3005  {
3006  yMin = y;
3007  }
3008  else if ( y > yMax )
3009  {
3010  yMax = y;
3011  }
3012  }
3013  }
3014 
3015  pointbar->setProgress( numberofpoints );
3016  delete pointbar;
3017  QApplication::restoreOverrideCursor();
3018 
3019 }
3020 
3021 bool DualEdgeTriangulation::saveToTAFF( QString filename ) const
3022 {
3023  QFile outputfile( filename );
3024  if ( !outputfile.open( IO_WriteOnly ) )
3025  {
3026  QMessageBox::warning( 0, tr( "Warning" ), tr( "File could not be written." ), QMessageBox::Ok, QMessageBox::NoButton, QMessageBox::NoButton );
3027  return false;
3028  }
3029 
3030  QTextStream outstream( &outputfile );
3031  outstream.precision( 9 );
3032 
3033  //export the edges. Attention, dual edges must be adjacent in the TAFF-file
3034  outstream << "TRIA" << std::endl << std::flush;
3035  outstream << "NEDG " << mHalfEdge.count() << std::endl << std::flush;
3036  outstream << "PANO 1" << std::endl << std::flush;
3037  outstream << "DATA ";
3038 
3039  bool *cont = new bool[mHalfEdge.count()];
3040  for ( unsigned int i = 0; i <= mHalfEdge.count() - 1; i++ )
3041  {
3042  cont[i] = false;
3043  }
3044 
3045  for ( unsigned int i = 0; i < mHalfEdge.count(); i++ )
3046  {
3047  if ( cont[i] )
3048  {
3049  continue;
3050  }
3051 
3052  int dual = mHalfEdge[i]->getDual();
3053  outstream << i << " " << mHalfEdge[i]->getPoint() << " " << mHalfEdge[i]->getNext() << " " << mHalfEdge[i]->getForced() << " " << mHalfEdge[i]->getBreak() << " ";
3054  outstream << dual << " " << mHalfEdge[dual]->getPoint() << " " << mHalfEdge[dual]->getNext() << " " << mHalfEdge[dual]->getForced() << " " << mHalfEdge[dual]->getBreak() << " ";
3055  cont[i] = true;
3056  cont[dual] = true;
3057  }
3058  outstream << std::endl << std::flush;
3059  outstream << std::endl << std::flush;
3060 
3061  delete[] cont;
3062 
3063  //export the points to the file
3064  outstream << "POIN" << std::endl << std::flush;
3065  outstream << "NPTS " << getNumberOfPoints() << std::endl << std::flush;
3066  outstream << "PATT 3" << std::endl << std::flush;
3067  outstream << "DATA ";
3068 
3069  for ( int i = 0; i < getNumberOfPoints(); i++ )
3070  {
3071  QgsPoint *p = mPointVector[i];
3072  outstream << p->getX() << " " << p->getY() << " " << p->getZ() << " ";
3073  }
3074  outstream << std::endl << std::flush;
3075  outstream << std::endl << std::flush;
3076 
3077  return true;
3078 }
3079 #endif //0
3080 
3081 bool QgsDualEdgeTriangulation::swapEdge( double x, double y )
3082 {
3083  QgsPoint p( x, y, 0 );
3084  const int edge1 = baseEdgeOfTriangle( p );
3085  if ( edge1 >= 0 )
3086  {
3087  int edge2, edge3;
3088  QgsPoint *point1 = nullptr;
3089  QgsPoint *point2 = nullptr;
3090  QgsPoint *point3 = nullptr;
3091  edge2 = mHalfEdge[edge1]->getNext();
3092  edge3 = mHalfEdge[edge2]->getNext();
3093  point1 = point( mHalfEdge[edge1]->getPoint() );
3094  point2 = point( mHalfEdge[edge2]->getPoint() );
3095  point3 = point( mHalfEdge[edge3]->getPoint() );
3096  if ( point1 && point2 && point3 )
3097  {
3098  //find out the closest edge to the point and swap this edge
3099  double dist1, dist2, dist3;
3100  dist1 = MathUtils::distPointFromLine( &p, point3, point1 );
3101  dist2 = MathUtils::distPointFromLine( &p, point1, point2 );
3102  dist3 = MathUtils::distPointFromLine( &p, point2, point3 );
3103  if ( dist1 <= dist2 && dist1 <= dist3 )
3104  {
3105  //qWarning("edge "+QString::number(edge1)+" is closest");
3106  if ( swapPossible( edge1 ) )
3107  {
3108  doOnlySwap( edge1 );
3109  }
3110  }
3111  else if ( dist2 <= dist1 && dist2 <= dist3 )
3112  {
3113  //qWarning("edge "+QString::number(edge2)+" is closest");
3114  if ( swapPossible( edge2 ) )
3115  {
3116  doOnlySwap( edge2 );
3117  }
3118  }
3119  else if ( dist3 <= dist1 && dist3 <= dist2 )
3120  {
3121  //qWarning("edge "+QString::number(edge3)+" is closest");
3122  if ( swapPossible( edge3 ) )
3123  {
3124  doOnlySwap( edge3 );
3125  }
3126  }
3127  return true;
3128  }
3129  else
3130  {
3131  // QgsDebugMsg("warning: null pointer");
3132  return false;
3133  }
3134  }
3135  else
3136  {
3137  // QgsDebugMsg("Edge number negative");
3138  return false;
3139  }
3140 }
3141 
3142 QList<int> QgsDualEdgeTriangulation::pointsAroundEdge( double x, double y )
3143 {
3144  QgsPoint p( x, y, 0 );
3145  QList<int> list;
3146  const int edge1 = baseEdgeOfTriangle( p );
3147  if ( edge1 >= 0 )
3148  {
3149  const int edge2 = mHalfEdge[edge1]->getNext();
3150  const int edge3 = mHalfEdge[edge2]->getNext();
3151  QgsPoint *point1 = point( mHalfEdge[edge1]->getPoint() );
3152  QgsPoint *point2 = point( mHalfEdge[edge2]->getPoint() );
3153  QgsPoint *point3 = point( mHalfEdge[edge3]->getPoint() );
3154  if ( point1 && point2 && point3 )
3155  {
3156  int p1, p2, p3, p4;
3157  const double dist1 = MathUtils::distPointFromLine( &p, point3, point1 );
3158  const double dist2 = MathUtils::distPointFromLine( &p, point1, point2 );
3159  const double dist3 = MathUtils::distPointFromLine( &p, point2, point3 );
3160  if ( dist1 <= dist2 && dist1 <= dist3 )
3161  {
3162  p1 = mHalfEdge[edge1]->getPoint();
3163  p2 = mHalfEdge[mHalfEdge[edge1]->getNext()]->getPoint();
3164  p3 = mHalfEdge[mHalfEdge[edge1]->getDual()]->getPoint();
3165  p4 = mHalfEdge[mHalfEdge[mHalfEdge[edge1]->getDual()]->getNext()]->getPoint();
3166  }
3167  else if ( dist2 <= dist1 && dist2 <= dist3 )
3168  {
3169  p1 = mHalfEdge[edge2]->getPoint();
3170  p2 = mHalfEdge[mHalfEdge[edge2]->getNext()]->getPoint();
3171  p3 = mHalfEdge[mHalfEdge[edge2]->getDual()]->getPoint();
3172  p4 = mHalfEdge[mHalfEdge[mHalfEdge[edge2]->getDual()]->getNext()]->getPoint();
3173  }
3174  else /* if ( dist3 <= dist1 && dist3 <= dist2 ) */
3175  {
3176  p1 = mHalfEdge[edge3]->getPoint();
3177  p2 = mHalfEdge[mHalfEdge[edge3]->getNext()]->getPoint();
3178  p3 = mHalfEdge[mHalfEdge[edge3]->getDual()]->getPoint();
3179  p4 = mHalfEdge[mHalfEdge[mHalfEdge[edge3]->getDual()]->getNext()]->getPoint();
3180  }
3181 
3182 
3183  list.append( p1 );
3184  list.append( p2 );
3185  list.append( p3 );
3186  list.append( p4 );
3187  }
3188  else
3189  {
3190  QgsDebugMsg( QStringLiteral( "warning: null pointer" ) );
3191  }
3192  }
3193  else
3194  {
3195  QgsDebugMsg( QStringLiteral( "Edge number negative" ) );
3196  }
3197  return list;
3198 }
3199 
3201 {
3202  if ( !sink )
3203  {
3204  return false;
3205  }
3206 
3207  bool *alreadyVisitedEdges = new bool[mHalfEdge.size()];
3208  if ( !alreadyVisitedEdges )
3209  {
3210  QgsDebugMsg( QStringLiteral( "out of memory" ) );
3211  return false;
3212  }
3213 
3214  for ( int i = 0; i < mHalfEdge.size(); ++i )
3215  {
3216  alreadyVisitedEdges[i] = false;
3217  }
3218 
3219  for ( int i = 0; i < mHalfEdge.size(); ++i )
3220  {
3221  if ( feedback && feedback->isCanceled() )
3222  break;
3223 
3224  HalfEdge *currentEdge = mHalfEdge[i];
3225  if ( currentEdge->getPoint() != -1 && mHalfEdge[currentEdge->getDual()]->getPoint() != -1 && !alreadyVisitedEdges[currentEdge->getDual()] )
3226  {
3227  QgsFeature edgeLineFeature;
3228 
3229  //geometry
3230  QgsPoint *p1 = mPointVector[currentEdge->getPoint()];
3231  QgsPoint *p2 = mPointVector[mHalfEdge[currentEdge->getDual()]->getPoint()];
3232  QgsPolylineXY lineGeom;
3233  lineGeom.push_back( QgsPointXY( p1->x(), p1->y() ) );
3234  lineGeom.push_back( QgsPointXY( p2->x(), p2->y() ) );
3235  edgeLineFeature.setGeometry( QgsGeometry::fromPolylineXY( lineGeom ) );
3236  edgeLineFeature.initAttributes( 1 );
3237 
3238  //attributes
3239  QString attributeString;
3240  if ( currentEdge->getForced() )
3241  {
3242  if ( currentEdge->getBreak() )
3243  {
3244  attributeString = QStringLiteral( "break line" );
3245  }
3246  else
3247  {
3248  attributeString = QStringLiteral( "structure line" );
3249  }
3250  }
3251  edgeLineFeature.setAttribute( 0, attributeString );
3252 
3253  sink->addFeature( edgeLineFeature, QgsFeatureSink::FastInsert );
3254  }
3255  alreadyVisitedEdges[i] = true;
3256  }
3257 
3258  delete [] alreadyVisitedEdges;
3259 
3260  return !feedback || !feedback->isCanceled();
3261 }
3262 
3264 {
3265  const QVector<bool> alreadyVisitedEdges( mHalfEdge.count(), false );
3266 
3267  if ( feedback )
3268  feedback->setProgress( 0 );
3269 
3270  QVector< bool> edgeToTreat( mHalfEdge.count(), true );
3271  QHash<HalfEdge *, int > edgesHash;
3272  for ( int i = 0; i < mHalfEdge.count(); ++i )
3273  {
3274  edgesHash.insert( mHalfEdge[i], i );
3275  }
3276 
3277  QgsMesh mesh;
3278  for ( const QgsPoint *point : mPointVector )
3279  {
3280  mesh.vertices.append( *point );
3281  }
3282 
3283  const int edgeCount = edgeToTreat.count();
3284  for ( int i = 0 ; i < edgeCount; ++i )
3285  {
3286  bool containVirtualPoint = false;
3287  if ( edgeToTreat[i] )
3288  {
3289  HalfEdge *currentEdge = mHalfEdge[i];
3290  HalfEdge *firstEdge = currentEdge;
3291  QgsMeshFace face;
3292  do
3293  {
3294  edgeToTreat[edgesHash.value( currentEdge )] = false;
3295  face.append( currentEdge->getPoint() );
3296  containVirtualPoint |= currentEdge->getPoint() == -1;
3297  currentEdge = mHalfEdge.at( currentEdge->getNext() );
3298  }
3299  while ( currentEdge != firstEdge && !containVirtualPoint && ( !feedback || !feedback->isCanceled() ) );
3300  if ( !containVirtualPoint )
3301  mesh.faces.append( face );
3302  }
3303  if ( feedback )
3304  {
3305  feedback->setProgress( ( 100 * i ) / edgeCount ) ;
3306  }
3307  }
3308 
3309  return mesh;
3310 }
3311 
3312 double QgsDualEdgeTriangulation::swapMinAngle( int edge ) const
3313 {
3314  QgsPoint *p1 = point( mHalfEdge[edge]->getPoint() );
3315  QgsPoint *p2 = point( mHalfEdge[mHalfEdge[edge]->getNext()]->getPoint() );
3316  QgsPoint *p3 = point( mHalfEdge[mHalfEdge[edge]->getDual()]->getPoint() );
3317  QgsPoint *p4 = point( mHalfEdge[mHalfEdge[mHalfEdge[edge]->getDual()]->getNext()]->getPoint() );
3318 
3319  //search for the minimum angle (it is important, which directions the lines have!)
3320  double minangle;
3321  const double angle1 = MathUtils::angle( p1, p2, p4, p2 );
3322  minangle = angle1;
3323  const double angle2 = MathUtils::angle( p3, p2, p4, p2 );
3324  if ( angle2 < minangle )
3325  {
3326  minangle = angle2;
3327  }
3328  const double angle3 = MathUtils::angle( p2, p3, p4, p3 );
3329  if ( angle3 < minangle )
3330  {
3331  minangle = angle3;
3332  }
3333  const double angle4 = MathUtils::angle( p3, p4, p2, p4 );
3334  if ( angle4 < minangle )
3335  {
3336  minangle = angle4;
3337  }
3338  const double angle5 = MathUtils::angle( p2, p4, p1, p4 );
3339  if ( angle5 < minangle )
3340  {
3341  minangle = angle5;
3342  }
3343  const double angle6 = MathUtils::angle( p4, p1, p2, p1 );
3344  if ( angle6 < minangle )
3345  {
3346  minangle = angle6;
3347  }
3348 
3349  return minangle;
3350 }
3351 
3352 int QgsDualEdgeTriangulation::splitHalfEdge( int edge, float position )
3353 {
3354  //just a short test if position is between 0 and 1
3355  if ( position < 0 || position > 1 )
3356  {
3357  QgsDebugMsg( QStringLiteral( "warning, position is not between 0 and 1" ) );
3358  }
3359 
3360  //create the new point on the heap
3361  QgsPoint *p = new QgsPoint( mPointVector[mHalfEdge[edge]->getPoint()]->x()*position + mPointVector[mHalfEdge[mHalfEdge[edge]->getDual()]->getPoint()]->x() * ( 1 - position ), mPointVector[mHalfEdge[edge]->getPoint()]->y()*position + mPointVector[mHalfEdge[mHalfEdge[edge]->getDual()]->getPoint()]->y() * ( 1 - position ), 0 );
3362 
3363  //calculate the z-value of the point to insert
3364  QgsPoint zvaluepoint( 0, 0, 0 );
3365  calcPoint( p->x(), p->y(), zvaluepoint );
3366  p->setZ( zvaluepoint.z() );
3367 
3368  //insert p into mPointVector
3369  if ( mPointVector.count() >= mPointVector.size() )
3370  {
3371  mPointVector.resize( mPointVector.count() + 1 );
3372  }
3373  QgsDebugMsg( QStringLiteral( "inserting point nr. %1, %2//%3//%4" ).arg( mPointVector.count() ).arg( p->x() ).arg( p->y() ).arg( p->z() ) );
3374  mPointVector.insert( mPointVector.count(), p );
3375 
3376  //insert the six new halfedges
3377  const int dualedge = mHalfEdge[edge]->getDual();
3378  const int edge1 = insertEdge( -10, -10, mPointVector.count() - 1, false, false );
3379  const int edge2 = insertEdge( edge1, mHalfEdge[mHalfEdge[edge]->getNext()]->getNext(), mHalfEdge[mHalfEdge[edge]->getNext()]->getPoint(), false, false );
3380  const int edge3 = insertEdge( -10, mHalfEdge[mHalfEdge[dualedge]->getNext()]->getNext(), mHalfEdge[mHalfEdge[dualedge]->getNext()]->getPoint(), false, false );
3381  const int edge4 = insertEdge( edge3, dualedge, mPointVector.count() - 1, false, false );
3382  const int edge5 = insertEdge( -10, mHalfEdge[edge]->getNext(), mHalfEdge[edge]->getPoint(), mHalfEdge[edge]->getBreak(), mHalfEdge[edge]->getForced() );
3383  const int edge6 = insertEdge( edge5, edge3, mPointVector.count() - 1, mHalfEdge[dualedge]->getBreak(), mHalfEdge[dualedge]->getForced() );
3384  mHalfEdge[edge1]->setDual( edge2 );
3385  mHalfEdge[edge1]->setNext( edge5 );
3386  mHalfEdge[edge3]->setDual( edge4 );
3387  mHalfEdge[edge5]->setDual( edge6 );
3388 
3389  //adjust the already existing halfedges
3390  mHalfEdge[mHalfEdge[edge]->getNext()]->setNext( edge1 );
3391  mHalfEdge[mHalfEdge[dualedge]->getNext()]->setNext( edge4 );
3392  mHalfEdge[edge]->setNext( edge2 );
3393  mHalfEdge[edge]->setPoint( mPointVector.count() - 1 );
3394  mHalfEdge[mHalfEdge[edge3]->getNext()]->setNext( edge6 );
3395 
3396  //test four times recursively for swapping
3397  checkSwapRecursively( mHalfEdge[edge5]->getNext(), 0 );
3398  checkSwapRecursively( mHalfEdge[edge2]->getNext(), 0 );
3399  checkSwapRecursively( mHalfEdge[dualedge]->getNext(), 0 );
3400  checkSwapRecursively( mHalfEdge[edge3]->getNext(), 0 );
3401 
3402  addPoint( QgsPoint( p->x(), p->y(), 0 ) );//dirty hack to enforce update of decorators
3403 
3404  return mPointVector.count() - 1;
3405 }
3406 
3407 bool QgsDualEdgeTriangulation::edgeOnConvexHull( int edge )
3408 {
3409  return ( mHalfEdge[mHalfEdge[edge]->getNext()]->getPoint() == -1 || mHalfEdge[mHalfEdge[mHalfEdge[edge]->getDual()]->getNext()]->getPoint() == -1 );
3410 }
3411 
3412 void QgsDualEdgeTriangulation::evaluateInfluenceRegion( QgsPoint *point, int edge, QSet<int> &set )
3413 {
3414  if ( set.find( edge ) == set.end() )
3415  {
3416  set.insert( edge );
3417  }
3418  else//prevent endless loops
3419  {
3420  return;
3421  }
3422 
3423  if ( !mHalfEdge[edge]->getForced() && !edgeOnConvexHull( edge ) )
3424  {
3425  //test, if point is in the circle through both endpoints of edge and the endpoint of edge->dual->next->point
3426  if ( inCircle( *point, *mPointVector[mHalfEdge[mHalfEdge[edge]->getDual()]->getPoint()], *mPointVector[mHalfEdge[edge]->getPoint()], *mPointVector[mHalfEdge[mHalfEdge[edge]->getNext()]->getPoint()] ) )
3427  {
3428  evaluateInfluenceRegion( point, mHalfEdge[mHalfEdge[edge]->getDual()]->getNext(), set );
3429  evaluateInfluenceRegion( point, mHalfEdge[mHalfEdge[mHalfEdge[edge]->getDual()]->getNext()]->getNext(), set );
3430  }
3431  }
3432 }
3433 
3434 int QgsDualEdgeTriangulation::firstEdgeOutSide()
3435 {
3436  int edge = 0;
3437  while ( ( mHalfEdge[mHalfEdge[edge]->getNext()]->getPoint() != -1 ||
3438  mHalfEdge[edge]->getPoint() == -1 || mHalfEdge[mHalfEdge[edge]->getDual()]->getPoint() == -1 ) &&
3439  edge < mHalfEdge.count() )
3440  {
3441  edge++;
3442  }
3443 
3444  if ( edge >= mHalfEdge.count() )
3445  return -1;
3446  else
3447  return edge;
3448 }
3449 
3450 void QgsDualEdgeTriangulation::removeLastPoint()
3451 {
3452  if ( mPointVector.isEmpty() )
3453  return;
3454  QgsPoint *p = mPointVector.takeLast();
3455  delete p;
3456 }
HalfEdge.
Definition: HalfEdge.h:31
bool getForced() const
Returns, whether the HalfEdge belongs to a constrained edge or not.
Definition: HalfEdge.h:98
int getNext() const
Returns the number of the next HalfEdge.
Definition: HalfEdge.h:83
void setNext(int n)
Sets the number of the next HalfEdge.
Definition: HalfEdge.h:108
void setPoint(int p)
Sets the number of point at which this HalfEdge points.
Definition: HalfEdge.h:113
int getPoint() const
Returns the number of the point at which this HalfEdge points.
Definition: HalfEdge.h:88
int getDual() const
Returns the number of the dual HalfEdge.
Definition: HalfEdge.h:78
void setDual(int d)
Sets the number of the dual HalfEdge.
Definition: HalfEdge.h:103
void setForced(bool f)
Sets the forced flag.
Definition: HalfEdge.h:123
bool getBreak() const
Returns, whether the HalfEdge belongs to a break line or not.
Definition: HalfEdge.h:93
void setBreak(bool b)
Sets the break flag.
Definition: HalfEdge.h:118
bool saveTriangulation(QgsFeatureSink *sink, QgsFeedback *feedback=nullptr) const override
Saves the triangulation features to a feature sink.
QgsPoint * point(int i) const override
Draws the points, edges and the forced lines.
QList< int > surroundingTriangles(int pointno) override
Returns a value list with the information of the triangles surrounding (counterclockwise) a point.
void ruppertRefinement() override
Adds points to make the triangles better shaped (algorithm of ruppert)
bool pointInside(double x, double y) override
Returns true, if the point with coordinates x and y is inside the convex hull and false otherwise.
void setTriangleInterpolator(TriangleInterpolator *interpolator) override
Sets an interpolator object.
void performConsistencyTest() override
Performs a consistency check, remove this later.
int oppositePoint(int p1, int p2) override
Returns the number of the point opposite to the triangle points p1, p2 (which have to be on a halfedg...
int addPoint(const QgsPoint &p) override
Adds a point to the triangulation.
virtual QgsMesh triangulationToMesh(QgsFeedback *feedback=nullptr) const override
Returns a QgsMesh corresponding to the triangulation.
void setForcedCrossBehavior(QgsTriangulation::ForcedCrossBehavior b) override
Sets the behavior of the triangulation in case of crossing forced lines.
bool calcNormal(double x, double y, QgsPoint &result) override
Calculates the normal at a point on the surface.
void addLine(const QVector< QgsPoint > &points, QgsInterpolator::SourceType lineType) override
Adds a line (e.g.
bool calcPoint(double x, double y, QgsPoint &result) override
Calculates x-, y and z-value of the point on the surface and assigns it to 'result'.
void eliminateHorizontalTriangles() override
Eliminates the horizontal triangles by swapping or by insertion of new points.
bool triangleVertices(double x, double y, QgsPoint &p1, int &n1, QgsPoint &p2, int &n2, QgsPoint &p3, int &n3) override
Finds out in which triangle the point with coordinates x and y is and assigns the numbers of the vert...
bool swapEdge(double x, double y) override
Swaps the edge which is closest to the point with x and y coordinates (if this is possible)
QList< int > pointsAroundEdge(double x, double y) override
Returns a value list with the numbers of the four points, which would be affected by an edge swap....
An interface for objects which accept features via addFeature(s) methods.
virtual bool addFeature(QgsFeature &feature, QgsFeatureSink::Flags flags=QgsFeatureSink::Flags())
Adds a single feature to the sink.
@ FastInsert
Use faster inserts, at the cost of updating the passed features to reflect changes made at the provid...
The feature class encapsulates a single feature including its unique ID, geometry and a list of field...
Definition: qgsfeature.h:56
bool setAttribute(int field, const QVariant &attr)
Sets an attribute's value by field index.
Definition: qgsfeature.cpp:237
void initAttributes(int fieldCount)
Initialize this feature with the given number of fields.
Definition: qgsfeature.cpp:210
void setGeometry(const QgsGeometry &geometry)
Set the feature's geometry.
Definition: qgsfeature.cpp:145
Base class for feedback objects to be used for cancellation of something running in a worker thread.
Definition: qgsfeedback.h:45
bool isCanceled() const SIP_HOLDGIL
Tells whether the operation has been canceled already.
Definition: qgsfeedback.h:54
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition: qgsfeedback.h:63
static QgsGeometry fromPolylineXY(const QgsPolylineXY &polyline)
Creates a new LineString geometry from a list of QgsPointXY points.
SourceType
Describes the type of input data.
@ SourceBreakLines
Break lines.
A class to represent a 2D point.
Definition: qgspointxy.h:59
Point geometry type, with support for z-dimension and m-values.
Definition: qgspoint.h:49
double distance(double x, double y) const SIP_HOLDGIL
Returns the Cartesian 2D distance between this point and a specified x, y coordinate.
Definition: qgspoint.h:343
void setX(double x) SIP_HOLDGIL
Sets the point's x-coordinate.
Definition: qgspoint.h:280
Q_GADGET double x
Definition: qgspoint.h:52
void setY(double y) SIP_HOLDGIL
Sets the point's y-coordinate.
Definition: qgspoint.h:291
void setZ(double z) SIP_HOLDGIL
Sets the point's z-coordinate.
Definition: qgspoint.h:304
double z
Definition: qgspoint.h:54
double distanceSquared(double x, double y) const SIP_HOLDGIL
Returns the Cartesian 2D squared distance between this point a specified x, y coordinate.
Definition: qgspoint.h:367
double y
Definition: qgspoint.h:53
ForcedCrossBehavior
Enumeration describing the behavior, if two forced lines cross.
@ SnappingTypeVertex
The second inserted forced line is snapped to the closest vertice of the first inserted forced line.
This is an interface for interpolator classes for triangulations.
virtual bool calcNormVec(double x, double y, QgsPoint &result)=0
Calculates the normal vector and assigns it to vec.
virtual bool calcPoint(double x, double y, QgsPoint &result)=0
Performs a linear interpolation in a triangle and assigns the x-,y- and z-coordinates to point.
bool ANALYSIS_EXPORT inDiametral(QgsPoint *p1, QgsPoint *p2, QgsPoint *point)
Tests, whether 'point' is inside the diametral circle through 'p1' and 'p2'.
Definition: MathUtils.cpp:266
double ANALYSIS_EXPORT angle(QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, QgsPoint *p4)
Calculates the angle between two segments (in 2 dimension, z-values are ignored)
Definition: MathUtils.cpp:786
bool ANALYSIS_EXPORT lineIntersection(QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, QgsPoint *p4)
Returns true, if line1 (p1 to p2) and line2 (p3 to p4) intersect. If the lines have an endpoint in co...
Definition: MathUtils.cpp:311
bool ANALYSIS_EXPORT circumcenter(QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, QgsPoint *result)
Calculates the center of the circle passing through p1, p2 and p3. Returns true in case of success an...
Definition: MathUtils.cpp:116
double ANALYSIS_EXPORT leftOf(const QgsPoint &thepoint, const QgsPoint *p1, const QgsPoint *p2)
Returns whether 'thepoint' is left or right of the line from 'p1' to 'p2'. Negative values mean left ...
Definition: MathUtils.cpp:292
double ANALYSIS_EXPORT distPointFromLine(QgsPoint *thepoint, QgsPoint *p1, QgsPoint *p2)
Calculates the (2 dimensional) distance from 'thepoint' to the line defined by p1 and p2.
Definition: MathUtils.cpp:201
bool ANALYSIS_EXPORT inCircle(QgsPoint *testp, QgsPoint *p1, QgsPoint *p2, QgsPoint *p3)
Tests, whether 'testp' is inside the circle through 'p1', 'p2' and 'p3'.
Definition: MathUtils.cpp:226
As part of the API refactoring and improvements which landed in the Processing API was substantially reworked from the x version This was done in order to allow much of the underlying Processing framework to be ported into c
double leftOfTresh
QVector< QgsPointXY > QgsPolylineXY
Polyline as represented as a vector of two-dimensional points.
Definition: qgsgeometry.h:51
#define QgsDebugMsg(str)
Definition: qgslogger.h:38
QVector< int > QgsMeshFace
List of vertex indexes.
Mesh - vertices, edges and faces.
QVector< QgsMeshVertex > vertices
QVector< QgsMeshFace > faces