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