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