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