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