30 const double x2 = point2.
x() - point1.
x();
31 const double y2 = point2.
y() - point1.
y();
32 const double x3 = point3.
x() - point1.
x();
33 const double y3 = point3.
y() - point1.
y();
35 const double denom = x2 * y3 - y2 * x3;
41 frac = ( x2 * ( x2 - x3 ) + y2 * ( y2 - y3 ) ) / denom;
42 const double cx = ( x3 + frac * y3 ) / 2;
43 const double cy = ( y3 - frac * x3 ) / 2;
44 const double squaredRadius = ( cx * cx + cy * cy );
45 const QgsPoint center( cx + point1.
x(), cy + point1.
y() );
52 qDeleteAll( mPointVector );
53 qDeleteAll( mHalfEdge );
60 for (
int i = 0; i < mHalfEdge.count(); i++ )
62 const int a = mHalfEdge[mHalfEdge[i]->getDual()]->getDual();
63 const int b = mHalfEdge[mHalfEdge[mHalfEdge[i]->getNext()]->getNext()]->getNext();
66 QgsDebugError( QStringLiteral(
"warning, first test failed" ) );
70 QgsDebugError( QStringLiteral(
"warning, second test failed" ) );
79 int currentpoint = -10;
86 if ( actpoint != -100 )
92 if ( actpoint == -100 )
97 for ( ; i < points.size(); ++i )
99 currentpoint =
addPoint( points.at( i ) );
100 if ( currentpoint != -100 && actpoint != -100 && currentpoint != actpoint )
102 insertForcedSegment( actpoint, currentpoint, lineType );
104 actpoint = currentpoint;
111 if ( mPointVector.isEmpty() )
120 mXMin = std::min( p.
x(), mXMin );
121 mXMax = std::max( p.
x(), mXMax );
122 mYMin = std::min( p.
y(), mYMin );
123 mYMax = std::max( p.
y(), mYMax );
127 mPointVector.append(
new QgsPoint( p ) );
130 if ( mDimension == -1 )
132 const unsigned int zedge = insertEdge( -10, -10, -1,
false,
false );
133 const unsigned int fedge = insertEdge(
static_cast<int>( zedge ),
static_cast<int>( zedge ), 0,
false,
false );
134 ( mHalfEdge.at( zedge ) )->setDual(
static_cast<int>( fedge ) );
135 ( mHalfEdge.at( zedge ) )->setNext(
static_cast<int>( fedge ) );
138 else if ( mDimension == 0 )
141 if ( p.
x() == mPointVector[0]->x() && p.
y() == mPointVector[0]->y() )
148 const unsigned int edgeFromPoint0ToPoint1 = insertEdge( -10, -10, 1,
false,
false );
149 const unsigned int edgeFromPoint1ToPoint0 = insertEdge( edgeFromPoint0ToPoint1, -10, 0,
false,
false );
150 const unsigned int edgeFromVirtualToPoint1Side1 = insertEdge( -10, -10, 1,
false,
false );
151 const unsigned int edgeFromPoint1ToVirtualSide1 = insertEdge( edgeFromVirtualToPoint1Side1, 1, -1,
false,
false );
152 const unsigned int edgeFromVirtualToPoint1Side2 = insertEdge( -10, edgeFromPoint1ToPoint0, 1,
false,
false );
153 const unsigned int edgeFromPoint1ToVirtualSide2 = insertEdge( edgeFromVirtualToPoint1Side2, edgeFromVirtualToPoint1Side1, -1,
false,
false );
154 const unsigned int edgeFromVirtualToPoint0Side2 = insertEdge( -10, -10, 0,
false,
false );
155 const unsigned int edgeFromPoint0ToVirtualSide2 = insertEdge( edgeFromVirtualToPoint0Side2, edgeFromVirtualToPoint1Side2, -1,
false,
false );
156 mHalfEdge.at( edgeFromPoint1ToPoint0 )->setNext( edgeFromPoint0ToVirtualSide2 );
157 mHalfEdge.at( edgeFromPoint0ToPoint1 )->setDual( edgeFromPoint1ToPoint0 );
158 mHalfEdge.at( edgeFromPoint0ToPoint1 )->setNext( edgeFromPoint1ToVirtualSide1 );
159 mHalfEdge.at( edgeFromVirtualToPoint1Side1 )->setDual( edgeFromPoint1ToVirtualSide1 );
160 mHalfEdge.at( edgeFromVirtualToPoint1Side1 )->setNext( edgeFromPoint1ToVirtualSide2 );
161 mHalfEdge.at( 0 )->setNext(
static_cast<int>( edgeFromVirtualToPoint0Side2 ) );
162 mHalfEdge.at( 1 )->setNext(
static_cast<int>( edgeFromPoint0ToPoint1 ) );
163 mHalfEdge.at( edgeFromVirtualToPoint1Side2 )->setDual( edgeFromPoint1ToVirtualSide2 );
164 mHalfEdge.at( edgeFromVirtualToPoint0Side2 )->setDual( edgeFromPoint0ToVirtualSide2 );
165 mHalfEdge.at( edgeFromVirtualToPoint0Side2 )->setNext( 0 );
167 mEdgeOutside = edgeFromPoint0ToPoint1;
170 else if ( mDimension == 1 )
172 if ( mEdgeOutside < 0 || mHalfEdge[mEdgeOutside]->getPoint() < 0 || mHalfEdge[mHalfEdge[mEdgeOutside]->getDual()]->getPoint() < 0 )
173 mEdgeOutside = firstEdgeOutSide();
174 if ( mEdgeOutside < 0 || mHalfEdge[mEdgeOutside]->getPoint() < 0 || mHalfEdge[mHalfEdge[mEdgeOutside]->getDual()]->getPoint() < 0 )
177 const double leftOfNumber =
MathUtils::leftOf( p, mPointVector[mHalfEdge[mHalfEdge[mEdgeOutside]->getDual()]->getPoint()], mPointVector[mHalfEdge[mEdgeOutside]->getPoint()] );
183 int closestEdge = -1;
184 double distance = std::numeric_limits<double>::max();
185 const int firstEdge = mEdgeOutside;
188 const int point1 = mHalfEdge[mEdgeOutside]->getPoint();
189 const int point2 = mHalfEdge[mHalfEdge[mEdgeOutside]->getDual()]->getPoint();
190 const double distance1 = p.
distance( *mPointVector[point1] );
196 const double distance2 = p.
distance( *mPointVector[point2] );
203 const double edgeLength = mPointVector[point1]->distance( *mPointVector[point2] );
205 if ( distance1 < edgeLength && distance2 < edgeLength )
208 const int newPoint = mPointVector.count() - 1;
211 const int edgeFromNewPointToPoint1 = mEdgeOutside;
212 const int edgeFromNewPointToPoint2 = mHalfEdge[mEdgeOutside]->getDual();
214 const int edgeFromPoint1ToVirtualSide2 = mHalfEdge[edgeFromNewPointToPoint1]->getNext();
215 const int edgeFromVirtualToPoint1Side1 = mHalfEdge[mHalfEdge[edgeFromNewPointToPoint2]->getNext()]->getNext();
216 const int edgeFromPoint2ToVirtualSide1 = mHalfEdge[edgeFromNewPointToPoint2]->getNext();
217 const int edgeFromVirtualToPoint2Side2 = mHalfEdge[mHalfEdge[edgeFromNewPointToPoint1]->getNext()]->getNext();
219 const int edgeFromVirtualToNewPointSide1 = insertEdge( -10, edgeFromNewPointToPoint2, newPoint,
false,
false );
220 const int edgeFromNewPointToVirtualSide1 = insertEdge( edgeFromVirtualToNewPointSide1, edgeFromVirtualToPoint1Side1, -1,
false,
false );
221 const int edgeFromVirtualToNewPointSide2 = insertEdge( -10, edgeFromNewPointToPoint1, newPoint,
false,
false );
222 const int edgeFromNewPointToVirtualSide2 = insertEdge( edgeFromVirtualToNewPointSide2, edgeFromVirtualToPoint2Side2, -1,
false,
false );
223 const int edgeFromPoint1ToNewPoint = insertEdge( edgeFromNewPointToPoint1, edgeFromNewPointToVirtualSide1, newPoint,
false,
false );
224 const int edgeFromPoint2ToNewPoint = insertEdge( edgeFromNewPointToPoint2, edgeFromNewPointToVirtualSide2, newPoint,
false,
false );
225 mHalfEdge.at( edgeFromVirtualToNewPointSide1 )->setDual( edgeFromNewPointToVirtualSide1 );
226 mHalfEdge.at( edgeFromVirtualToNewPointSide2 )->setDual( edgeFromNewPointToVirtualSide2 );
228 mHalfEdge.at( edgeFromPoint1ToVirtualSide2 )->setNext( edgeFromVirtualToNewPointSide2 );
229 mHalfEdge.at( edgeFromVirtualToPoint1Side1 )->setNext( edgeFromPoint1ToNewPoint );
230 mHalfEdge.at( edgeFromPoint2ToVirtualSide1 )->setNext( edgeFromVirtualToNewPointSide1 );
231 mHalfEdge.at( edgeFromVirtualToPoint2Side2 )->setNext( edgeFromPoint2ToNewPoint );
232 mHalfEdge.at( edgeFromNewPointToPoint1 )->setDual( edgeFromPoint1ToNewPoint );
233 mHalfEdge.at( edgeFromNewPointToPoint2 )->setDual( edgeFromPoint2ToNewPoint );
238 if ( distance1 < distance )
240 closestEdge = mEdgeOutside;
241 distance = distance1;
243 else if ( distance2 < distance )
245 closestEdge = mHalfEdge[mEdgeOutside]->getDual();
246 distance = distance2;
249 mEdgeOutside = mHalfEdge[mHalfEdge[mHalfEdge[mEdgeOutside]->getNext()]->getDual()]->getNext();
251 while ( mEdgeOutside != firstEdge && mHalfEdge[mEdgeOutside]->getPoint() != -1 && mHalfEdge[mHalfEdge[mEdgeOutside]->getDual()]->getPoint() != -1 );
253 if ( closestEdge < 0 )
257 const int extremePoint = mHalfEdge[closestEdge]->getPoint();
258 const int newPoint = mPointVector.count() - 1;
260 const int edgeFromExtremeToOpposite = mHalfEdge[closestEdge]->getDual();
262 const int edgeFromVirtualToExtremeSide1 = mHalfEdge[mHalfEdge[closestEdge]->getNext()]->getDual();
263 const int edgeFromVirtualToExtremeSide2 = mHalfEdge[mHalfEdge[mHalfEdge[closestEdge]->getDual()]->getNext()]->getNext();
264 const int edgeFromExtremeToVirtualSide2 = mHalfEdge[edgeFromVirtualToExtremeSide2]->getDual();
266 const int edgeFromExtremeToNewPoint = insertEdge( -10, -10, newPoint,
false,
false );
267 const int edgeFromNewPointToExtreme = insertEdge( edgeFromExtremeToNewPoint, edgeFromExtremeToVirtualSide2, extremePoint,
false,
false );
268 const int edgeFromNewPointToVirtualSide1 = insertEdge( -10, edgeFromVirtualToExtremeSide1, -1,
false,
false );
269 const int edgeFromVirtualToNewPointSide1 = insertEdge( edgeFromNewPointToVirtualSide1, -10, newPoint,
false,
false );
270 const int edgeFromNewPointToVirtualSide2 = insertEdge( -10, edgeFromVirtualToNewPointSide1, -1,
false,
false );
271 const int edgeFromVirtualToNewPointSide2 = insertEdge( edgeFromNewPointToVirtualSide2, edgeFromNewPointToExtreme, newPoint,
false,
false );
272 mHalfEdge.at( edgeFromExtremeToNewPoint )->setDual( edgeFromNewPointToExtreme );
273 mHalfEdge.at( edgeFromExtremeToNewPoint )->setNext( edgeFromNewPointToVirtualSide1 );
274 mHalfEdge.at( edgeFromNewPointToVirtualSide1 )->setDual( edgeFromVirtualToNewPointSide1 );
275 mHalfEdge.at( edgeFromNewPointToVirtualSide2 )->setDual( edgeFromVirtualToNewPointSide2 );
276 mHalfEdge.at( edgeFromVirtualToNewPointSide1 )->setNext( edgeFromNewPointToVirtualSide2 );
278 mHalfEdge.at( edgeFromVirtualToExtremeSide1 )->setNext( edgeFromExtremeToNewPoint );
279 mHalfEdge.at( edgeFromVirtualToExtremeSide2 )->setNext( edgeFromExtremeToOpposite );
280 mHalfEdge.at( edgeFromExtremeToVirtualSide2 )->setNext( edgeFromVirtualToNewPointSide2 );
287 mEdgeOutside = mHalfEdge[mEdgeOutside]->getDual();
290 const int newPoint = mPointVector.count() - 1;
293 int cwEdge = mEdgeOutside;
294 while ( mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[cwEdge]->getNext()]->getDual()]->getNext()]->getPoint() != -1 )
296 mHalfEdge[mHalfEdge[ cwEdge ]->getNext()]->setPoint( newPoint );
297 cwEdge = mHalfEdge[mHalfEdge[mHalfEdge[cwEdge]->getNext()]->getDual()]->getNext();
300 const int edgeFromLastCwPointToVirtualPoint = mHalfEdge[mHalfEdge[mHalfEdge[cwEdge]->getNext()]->getDual()]->getNext();
301 const int edgeFromLastCwPointToNewPointPoint = mHalfEdge[ cwEdge ]->getNext();
302 const int edgeFromNewPointPointToLastCwPoint = mHalfEdge[ edgeFromLastCwPointToNewPointPoint ]->getDual();
304 const int edgeFromNewPointtoVirtualPoint = insertEdge( -10, -10, -1,
false,
false );
305 const int edgeFromVirtualPointToNewPoint = insertEdge( edgeFromNewPointtoVirtualPoint, edgeFromNewPointPointToLastCwPoint, newPoint,
false,
false );
306 mHalfEdge.at( edgeFromLastCwPointToNewPointPoint )->setPoint( newPoint );
307 mHalfEdge.at( edgeFromNewPointtoVirtualPoint )->setDual( edgeFromVirtualPointToNewPoint );
308 mHalfEdge.at( edgeFromLastCwPointToVirtualPoint )->setNext( edgeFromVirtualPointToNewPoint );
311 int ccwEdge = mEdgeOutside;
312 while ( mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[ccwEdge]->getDual()]->getNext()]->getDual()]->getNext()]->getPoint() != -1 )
314 mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[ ccwEdge ]->getNext()]->getNext()]->getDual()]->setPoint( newPoint );
315 ccwEdge = mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[ccwEdge]->getDual()]->getNext()]->getDual()]->getNext()]->getDual();
318 const int edgeToLastCcwPoint = mHalfEdge[mHalfEdge[mHalfEdge[ccwEdge]->getDual()]->getNext()]->getDual();
319 const int edgeFromLastCcwPointToNewPoint = mHalfEdge[edgeToLastCcwPoint]->getNext();
320 mHalfEdge.at( edgeFromNewPointtoVirtualPoint )->setNext( edgeToLastCcwPoint );
321 mHalfEdge.at( edgeFromLastCcwPointToNewPoint )->setNext( edgeFromNewPointtoVirtualPoint );
322 mHalfEdge.at( edgeFromLastCcwPointToNewPoint )->setPoint( newPoint );
326 const int number = baseEdgeOfTriangle( p );
331 unsigned int cwEdge = mEdgeOutside;
332 unsigned int ccwEdge = mEdgeOutside;
335 mHalfEdge[mHalfEdge[mEdgeOutside]->getNext()]->setPoint( mPointVector.count() - 1 );
338 while (
MathUtils::leftOf( *mPointVector[ mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[cwEdge]->getNext()]->getDual()]->getNext()]->getPoint()],
339 &p, mPointVector[ mHalfEdge[cwEdge]->getPoint()] ) < ( -
leftOfTresh ) )
342 mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[cwEdge]->getNext()]->getDual()]->getNext()]->getNext()]->setPoint( mPointVector.count() - 1 );
344 cwEdge = (
unsigned int )mHalfEdge[mHalfEdge[mHalfEdge[cwEdge]->getNext()]->getDual()]->getNext();
348 const unsigned int edge1 = insertEdge( mHalfEdge[cwEdge]->getNext(), -10, mHalfEdge[cwEdge]->getPoint(),
false,
false );
349 const unsigned int edge2 = insertEdge( mHalfEdge[mHalfEdge[cwEdge]->getNext()]->getDual(), -10, -1,
false,
false );
350 const unsigned int edge3 = insertEdge( -10, edge1, mPointVector.count() - 1,
false,
false );
353 mHalfEdge[mHalfEdge[mHalfEdge[cwEdge]->getNext()]->getDual()]->setDual( edge2 );
354 mHalfEdge[mHalfEdge[cwEdge]->getNext()]->setDual( edge1 );
355 mHalfEdge[edge1]->setNext( edge2 );
356 mHalfEdge[edge2]->setNext( edge3 );
359 while (
MathUtils::leftOf( *mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[ccwEdge]->getNext()]->getNext()]->getPoint()], mPointVector[mPointVector.count() - 1], mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[ccwEdge]->getNext()]->getNext()]->getDual()]->getNext()]->getPoint()] ) < ( -
leftOfTresh ) )
362 mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[ccwEdge]->getNext()]->getNext()]->getDual()]->setPoint( mPointVector.count() - 1 );
364 ccwEdge = mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[ccwEdge]->getNext()]->getNext()]->getDual()]->getNext()]->getNext();
368 const unsigned int edge4 = insertEdge( mHalfEdge[mHalfEdge[ccwEdge]->getNext()]->getNext(), -10, mPointVector.count() - 1,
false,
false );
369 const unsigned int edge5 = insertEdge( edge3, -10, -1,
false,
false );
370 const unsigned int edge6 = insertEdge( mHalfEdge[mHalfEdge[mHalfEdge[ccwEdge]->getNext()]->getNext()]->getDual(), edge4, mHalfEdge[mHalfEdge[ccwEdge]->getDual()]->getPoint(),
false,
false );
375 mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[ccwEdge]->getNext()]->getNext()]->getDual()]->setDual( edge6 );
376 mHalfEdge[mHalfEdge[mHalfEdge[ccwEdge]->getNext()]->getNext()]->setDual( edge4 );
377 mHalfEdge[edge4]->setNext( edge5 );
378 mHalfEdge[edge5]->setNext( edge6 );
379 mHalfEdge[edge3]->setDual( edge5 );
382 unsigned int index = ccwEdge;
387 index = mHalfEdge[mHalfEdge[mHalfEdge[index]->getNext()]->getDual()]->getNext();
388 checkSwapRecursively( toswap, 0 );
389 if ( toswap == cwEdge )
395 else if ( number >= 0 )
397 const int nextnumber = mHalfEdge[number]->getNext();
398 const int nextnextnumber = mHalfEdge[mHalfEdge[number]->getNext()]->getNext();
401 const unsigned int edge1 = insertEdge( -10, nextnumber, mHalfEdge[number]->getPoint(),
false,
false );
402 const unsigned int edge2 = insertEdge(
static_cast<int>( edge1 ), -10, mPointVector.count() - 1,
false,
false );
403 const unsigned int edge3 = insertEdge( -10, nextnextnumber, mHalfEdge[nextnumber]->getPoint(),
false,
false );
404 const unsigned int edge4 = insertEdge(
static_cast<int>( edge3 ),
static_cast<int>( edge1 ), mPointVector.count() - 1,
false,
false );
405 const unsigned int edge5 = insertEdge( -10, number, mHalfEdge[nextnextnumber]->getPoint(),
false,
false );
406 const unsigned int edge6 = insertEdge(
static_cast<int>( edge5 ),
static_cast<int>( edge3 ), mPointVector.count() - 1,
false,
false );
409 mHalfEdge.at( edge1 )->setDual(
static_cast<int>( edge2 ) );
410 mHalfEdge.at( edge2 )->setNext(
static_cast<int>( edge5 ) );
411 mHalfEdge.at( edge3 )->setDual(
static_cast<int>( edge4 ) );
412 mHalfEdge.at( edge5 )->setDual(
static_cast<int>( edge6 ) );
413 mHalfEdge.at( number )->setNext(
static_cast<int>( edge2 ) );
414 mHalfEdge.at( nextnumber )->setNext(
static_cast<int>( edge4 ) );
415 mHalfEdge.at( nextnextnumber )->setNext(
static_cast<int>( edge6 ) );
418 checkSwapRecursively( number, 0 );
419 checkSwapRecursively( nextnumber, 0 );
420 checkSwapRecursively( nextnextnumber, 0 );
423 else if ( number == -20 )
428 const int point1 = mHalfEdge[mEdgeWithPoint]->getPoint();
429 const int point2 = mHalfEdge[mHalfEdge[mEdgeWithPoint]->getDual()]->getPoint();
430 const double distance1 = p.
distance( *mPointVector[point1] );
436 const double distance2 = p.
distance( *mPointVector[point2] );
443 const int edgea = mEdgeWithPoint;
444 const int edgeb = mHalfEdge[mEdgeWithPoint]->getDual();
445 const int edgec = mHalfEdge[edgea]->getNext();
446 const int edged = mHalfEdge[edgec]->getNext();
447 const int edgee = mHalfEdge[edgeb]->getNext();
448 const int edgef = mHalfEdge[edgee]->getNext();
451 const int nedge1 = insertEdge( -10, mHalfEdge[edgea]->getNext(), mHalfEdge[edgea]->getPoint(),
false,
false );
452 const int nedge2 = insertEdge( nedge1, -10, mPointVector.count() - 1,
false,
false );
453 const int nedge3 = insertEdge( -10, edged, mHalfEdge[edgec]->getPoint(),
false,
false );
454 const int nedge4 = insertEdge( nedge3, nedge1, mPointVector.count() - 1,
false,
false );
455 const int nedge5 = insertEdge( -10, edgef, mHalfEdge[edgee]->getPoint(),
false,
false );
456 const int nedge6 = insertEdge( nedge5, edgeb, mPointVector.count() - 1,
false,
false );
459 mHalfEdge[nedge1]->setDual( nedge2 );
460 mHalfEdge[nedge2]->setNext( nedge5 );
461 mHalfEdge[nedge3]->setDual( nedge4 );
462 mHalfEdge[nedge5]->setDual( nedge6 );
463 mHalfEdge[edgea]->setPoint( mPointVector.count() - 1 );
464 mHalfEdge[edgea]->setNext( nedge3 );
465 mHalfEdge[edgec]->setNext( nedge4 );
466 mHalfEdge[edgee]->setNext( nedge6 );
467 mHalfEdge[edgef]->setNext( nedge2 );
470 checkSwapRecursively( edgec, 0 );
471 checkSwapRecursively( edged, 0 );
472 checkSwapRecursively( edgee, 0 );
473 checkSwapRecursively( edgef, 0 );
476 else if ( number == -100 || number == -5 )
482 else if ( number == -25 )
485 QgsPoint *newPoint = mPointVector[mPointVector.count() - 1];
486 QgsPoint *existingPoint = mPointVector[mTwiceInsPoint];
487 existingPoint->
setZ( std::max( newPoint->
z(), existingPoint->
z() ) );
490 return mTwiceInsPoint;
494 return ( mPointVector.count() - 1 );
497int QgsDualEdgeTriangulation::baseEdgeOfPoint(
int point )
499 unsigned int actedge = mEdgeInside;
501 if ( mPointVector.count() < 4 ||
point == -1 || mDimension == 1 )
503 int fromVirtualPoint = -1;
505 for (
int i = 0; i < mHalfEdge.count(); i++ )
507 if ( mHalfEdge[i]->getPoint() ==
point )
509 if ( mHalfEdge[mHalfEdge[i]->getDual()]->getPoint() != -1 )
512 fromVirtualPoint = i;
515 return fromVirtualPoint;
523 if ( control > 1000000 )
529 for (
int i = 0; i < mHalfEdge.count(); i++ )
531 if ( mHalfEdge[i]->getPoint() ==
point && mHalfEdge[mHalfEdge[i]->getNext()]->getPoint() != -1 )
538 const int fromPoint = mHalfEdge[mHalfEdge[actedge]->getDual()]->getPoint();
539 const int toPoint = mHalfEdge[actedge]->getPoint();
541 if ( fromPoint == -1 || toPoint == -1 )
543 for (
int i = 0; i < mHalfEdge.count(); i++ )
545 if ( mHalfEdge[i]->getPoint() ==
point && mHalfEdge[mHalfEdge[i]->getNext()]->getPoint() != -1 )
553 const double leftOfNumber =
MathUtils::leftOf( *mPointVector[
point], mPointVector[mHalfEdge[mHalfEdge[actedge]->getDual()]->getPoint()], mPointVector[mHalfEdge[actedge]->getPoint()] );
556 if ( mHalfEdge[actedge]->getPoint() ==
point && mHalfEdge[mHalfEdge[actedge]->getNext()]->getPoint() != -1 )
558 mEdgeInside = actedge;
561 else if ( leftOfNumber <= 0.0 )
563 actedge = mHalfEdge[actedge]->getNext();
567 actedge = mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[actedge]->getDual()]->getNext()]->getNext()]->getDual();
572int QgsDualEdgeTriangulation::baseEdgeOfTriangle(
const QgsPoint &point )
574 unsigned int actEdge = mEdgeInside;
575 if ( mHalfEdge.at( actEdge )->getPoint() < 0 )
576 actEdge = mHalfEdge.at( mHalfEdge.at( mHalfEdge.at( actEdge )->getDual() )->getNext() )->getDual();
577 if ( mHalfEdge.at( mHalfEdge.at( actEdge )->getDual() )->getPoint() < 0 )
578 actEdge = mHalfEdge.at( mHalfEdge.at( actEdge )->getNext() )->getDual();
584 int firstEndPoint = 0, secEndPoint = 0, thEndPoint = 0, fouEndPoint = 0;
588 if ( runs > MAX_BASE_ITERATIONS )
594 const double leftOfValue =
MathUtils::leftOf(
point, mPointVector.at( mHalfEdge.at( mHalfEdge.at( actEdge )->getDual() )->getPoint() ), mPointVector.at( mHalfEdge.at( actEdge )->getPoint() ) );
609 firstEndPoint = mHalfEdge.at( mHalfEdge.at( actEdge )->getDual() )->getPoint();
610 secEndPoint = mHalfEdge.at( actEdge )->getPoint();
612 else if ( nulls == 1 )
615 thEndPoint = mHalfEdge.at( mHalfEdge.at( actEdge )->getDual() )->getPoint();
616 fouEndPoint = mHalfEdge.at( actEdge )->getPoint();
619 mEdgeWithPoint = actEdge;
628 actEdge = mHalfEdge.at( actEdge )->getDual();
633 actEdge = mHalfEdge.at( actEdge )->getNext();
634 if ( mHalfEdge.at( actEdge )->getPoint() == -1 )
640 mEdgeOutside = (
unsigned int )mHalfEdge.at( mHalfEdge.at( actEdge )->getNext() )->getNext();
641 mEdgeInside = mHalfEdge.at( mHalfEdge.at( mEdgeOutside )->getDual() )->getNext();
647 if ( numInstabs > 0 )
650 mUnstableEdge = actEdge;
657 if ( firstEndPoint == thEndPoint || firstEndPoint == fouEndPoint )
660 mTwiceInsPoint = firstEndPoint;
663 else if ( secEndPoint == thEndPoint || secEndPoint == fouEndPoint )
666 mTwiceInsPoint = secEndPoint;
678 mEdgeInside = actEdge;
681 nr1 = mHalfEdge.at( actEdge )->getPoint();
682 nr2 = mHalfEdge.at( mHalfEdge.at( actEdge )->getNext() )->getPoint();
683 nr3 = mHalfEdge.at( mHalfEdge.at( mHalfEdge.at( actEdge )->getNext() )->getNext() )->getPoint();
684 const double x1 = mPointVector.at( nr1 )->x();
685 const double y1 = mPointVector.at( nr1 )->y();
686 const double x2 = mPointVector.at( nr2 )->x();
687 const double y2 = mPointVector.at( nr2 )->y();
688 const double x3 = mPointVector.at( nr3 )->x();
689 const double y3 = mPointVector.at( nr3 )->y();
692 if ( x1 < x2 && x1 < x3 )
696 else if ( x2 < x1 && x2 < x3 )
698 return mHalfEdge.at( actEdge )->getNext();
700 else if ( x3 < x1 && x3 < x2 )
702 return mHalfEdge.at( mHalfEdge.at( actEdge )->getNext() )->getNext();
713 return mHalfEdge.at( actEdge )->getNext();
720 return mHalfEdge.at( actEdge )->getNext();
724 return mHalfEdge.at( mHalfEdge.at( actEdge )->getNext() )->getNext();
735 return mHalfEdge.at( mHalfEdge.at( actEdge )->getNext() )->getNext();
743 if ( mTriangleInterpolator )
745 return mTriangleInterpolator->
calcNormVec( x, y, result );
756 if ( mTriangleInterpolator )
758 return mTriangleInterpolator->
calcPoint( x, y, result );
767bool QgsDualEdgeTriangulation::checkSwapRecursively(
unsigned int edge,
unsigned int recursiveDeep )
769 if ( swapPossible( edge ) )
771 QgsPoint *pta = mPointVector.at( mHalfEdge.at( edge )->getPoint() );
772 QgsPoint *ptb = mPointVector.at( mHalfEdge.at( mHalfEdge.at( edge )->getNext() )->getPoint() );
773 QgsPoint *ptc = mPointVector.at( mHalfEdge.at( mHalfEdge.at( mHalfEdge.at( edge )->getNext( ) )->getNext() )->getPoint() );
774 QgsPoint *ptd = mPointVector.at( mHalfEdge.at( mHalfEdge.at( mHalfEdge.at( edge )->getDual() )->getNext() )->getPoint() );
775 if ( inCircle( *ptd, *pta, *ptb, *ptc ) )
777 doSwapRecursively( edge, recursiveDeep );
784bool QgsDualEdgeTriangulation::isEdgeNeedSwap(
unsigned int edge )
const
786 if ( swapPossible( edge ) )
788 QgsPoint *pta = mPointVector[mHalfEdge[edge]->getPoint()];
789 QgsPoint *ptb = mPointVector[mHalfEdge[mHalfEdge[edge]->getNext()]->getPoint()];
790 QgsPoint *ptc = mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[edge]->getNext()]->getNext()]->getPoint()];
791 QgsPoint *ptd = mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[edge]->getDual()]->getNext()]->getPoint()];
792 return inCircle( *ptd, *pta, *ptb, *ptc );
799void QgsDualEdgeTriangulation::doOnlySwap(
unsigned int edge )
801 const unsigned int edge1 = edge;
802 const unsigned int edge2 = mHalfEdge[edge]->getDual();
803 const unsigned int edge3 = mHalfEdge[edge]->getNext();
804 const unsigned int edge4 = mHalfEdge[mHalfEdge[edge]->getNext()]->getNext();
805 const unsigned int edge5 = mHalfEdge[mHalfEdge[edge]->getDual()]->getNext();
806 const unsigned int edge6 = mHalfEdge[mHalfEdge[mHalfEdge[edge]->getDual()]->getNext()]->getNext();
807 mHalfEdge[edge1]->setNext( edge4 );
808 mHalfEdge[edge2]->setNext( edge6 );
809 mHalfEdge[edge3]->setNext( edge2 );
810 mHalfEdge[edge4]->setNext( edge5 );
811 mHalfEdge[edge5]->setNext( edge1 );
812 mHalfEdge[edge6]->setNext( edge3 );
813 mHalfEdge[edge1]->setPoint( mHalfEdge[edge3]->getPoint() );
814 mHalfEdge[edge2]->setPoint( mHalfEdge[edge5]->getPoint() );
817void QgsDualEdgeTriangulation::doSwapRecursively(
unsigned int edge,
unsigned int recursiveDeep )
819 const unsigned int edge1 = edge;
820 const unsigned int edge2 = mHalfEdge.at( edge )->getDual();
821 const unsigned int edge3 = mHalfEdge.at( edge )->getNext();
822 const unsigned int edge4 = mHalfEdge.at( mHalfEdge.at( edge )->getNext() )->getNext();
823 const unsigned int edge5 = mHalfEdge.at( mHalfEdge.at( edge )->getDual() )->getNext();
824 const unsigned int edge6 = mHalfEdge.at( mHalfEdge.at( mHalfEdge.at( edge )->getDual() )->getNext() )->getNext();
825 mHalfEdge.at( edge1 )->setNext( edge4 );
826 mHalfEdge.at( edge2 )->setNext( edge6 );
827 mHalfEdge.at( edge3 )->setNext( edge2 );
828 mHalfEdge.at( edge4 )->setNext( edge5 );
829 mHalfEdge.at( edge5 )->setNext( edge1 );
830 mHalfEdge.at( edge6 )->setNext( edge3 );
831 mHalfEdge.at( edge1 )->setPoint( mHalfEdge.at( edge3 )->getPoint() );
832 mHalfEdge.at( edge2 )->setPoint( mHalfEdge.at( edge5 )->getPoint() );
835 if ( recursiveDeep < 100 )
837 checkSwapRecursively( edge3, recursiveDeep );
838 checkSwapRecursively( edge6, recursiveDeep );
839 checkSwapRecursively( edge4, recursiveDeep );
840 checkSwapRecursively( edge5, recursiveDeep );
844 QStack<int> edgesToSwap;
845 edgesToSwap.push( edge3 );
846 edgesToSwap.push( edge6 );
847 edgesToSwap.push( edge4 );
848 edgesToSwap.push( edge5 );
850 while ( !edgesToSwap.isEmpty() && loopCount < 10000 )
853 const unsigned int e1 = edgesToSwap.pop();
854 if ( isEdgeNeedSwap( e1 ) )
856 const unsigned int e2 = mHalfEdge.at( e1 )->getDual();
857 const unsigned int e3 = mHalfEdge.at( e1 )->getNext();
858 const unsigned int e4 = mHalfEdge.at( mHalfEdge.at( e1 )->getNext() )->getNext();
859 const unsigned int e5 = mHalfEdge.at( mHalfEdge.at( e1 )->getDual() )->getNext();
860 const unsigned int e6 = mHalfEdge.at( mHalfEdge.at( mHalfEdge.at( e1 )->getDual() )->getNext() )->getNext();
861 mHalfEdge.at( e1 )->setNext( e4 );
862 mHalfEdge.at( e2 )->setNext( e6 );
863 mHalfEdge.at( e3 )->setNext( e2 );
864 mHalfEdge.at( e4 )->setNext( e5 );
865 mHalfEdge.at( e5 )->setNext( e1 );
866 mHalfEdge.at( e6 )->setNext( e3 );
867 mHalfEdge.at( e1 )->setPoint( mHalfEdge.at( e3 )->getPoint() );
868 mHalfEdge.at( e2 )->setPoint( mHalfEdge.at( e5 )->getPoint() );
870 edgesToSwap.push( e3 );
871 edgesToSwap.push( e6 );
872 edgesToSwap.push( e4 );
873 edgesToSwap.push( e5 );
881void DualEdgeTriangulation::draw( QPainter *p,
double xlowleft,
double ylowleft,
double xupright,
double yupright,
double width,
double height )
const
884 if ( mPointVector.isEmpty() )
889 p->setPen( mEdgeColor );
891 bool *control =
new bool[mHalfEdge.count()];
892 bool *control2 =
new bool[mHalfEdge.count()];
894 for (
unsigned int i = 0; i <= mHalfEdge.count() - 1; i++ )
900 if ( ( ( xupright - xlowleft ) / width ) > ( ( yupright - ylowleft ) / height ) )
902 double lowerborder = -( height * ( xupright - xlowleft ) / width - yupright );
903 for (
unsigned int i = 0; i < mHalfEdge.count() - 1; i++ )
905 if ( mHalfEdge[i]->getPoint() == -1 || mHalfEdge[mHalfEdge[i]->getDual()]->getPoint() == -1 )
912 if ( mHalfEdge[i]->getPoint() != -1 && mHalfEdge[mHalfEdge[i]->getNext()]->getPoint() != -1 && mHalfEdge[mHalfEdge[mHalfEdge[i]->getNext()]->getNext()]->getPoint() != -1 )
914 p1 = mPointVector[mHalfEdge[i]->getPoint()]->getZ();
915 p2 = mPointVector[mHalfEdge[mHalfEdge[i]->getNext()]->getPoint()]->getZ();
916 p3 = mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[i]->getNext()]->getNext()]->getPoint()]->getZ();
917 if ( p1 == p2 && p2 == p3 && halfEdgeBBoxTest( i, xlowleft, lowerborder, xupright, yupright ) && halfEdgeBBoxTest( mHalfEdge[i]->getNext(), xlowleft, lowerborder, xupright, yupright ) && halfEdgeBBoxTest( mHalfEdge[mHalfEdge[i]->getNext()]->getNext(), xlowleft, lowerborder, xupright, yupright ) )
920 pa.setPoint( 0, ( mPointVector[mHalfEdge[i]->getPoint()]->getX() - xlowleft ) / ( xupright - xlowleft )*width, ( yupright - mPointVector[mHalfEdge[i]->getPoint()]->getY() ) / ( xupright - xlowleft )*width );
921 pa.setPoint( 1, ( mPointVector[mHalfEdge[mHalfEdge[i]->getNext()]->getPoint()]->getX() - xlowleft ) / ( xupright - xlowleft )*width, ( yupright - mPointVector[mHalfEdge[mHalfEdge[i]->getNext()]->getPoint()]->getY() ) / ( xupright - xlowleft )*width );
922 pa.setPoint( 2, ( mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[i]->getNext()]->getNext()]->getPoint()]->getX() - xlowleft ) / ( xupright - xlowleft )*width, ( yupright - mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[i]->getNext()]->getNext()]->getPoint()]->getY() ) / ( xupright - xlowleft )*width );
923 QColor
c( 255, 0, 0 );
925 p->drawPolygon( pa );
930 control2[mHalfEdge[i]->getNext()] =
true;
931 control2[mHalfEdge[mHalfEdge[i]->getNext()]->getNext()] =
true;
938 if ( halfEdgeBBoxTest( i, xlowleft, lowerborder, xupright, yupright ) )
940 if ( mHalfEdge[i]->getBreak() )
942 p->setPen( mBreakEdgeColor );
944 else if ( mHalfEdge[i]->getForced() )
946 p->setPen( mForcedEdgeColor );
950 p->drawLine( ( mPointVector[mHalfEdge[i]->getPoint()]->getX() - xlowleft ) / ( xupright - xlowleft )*width, ( yupright - mPointVector[mHalfEdge[i]->getPoint()]->getY() ) / ( xupright - xlowleft )*width, ( mPointVector[mHalfEdge[mHalfEdge[i]->getDual()]->getPoint()]->getX() - xlowleft ) / ( xupright - xlowleft )*width, ( yupright - mPointVector[mHalfEdge[mHalfEdge[i]->getDual()]->getPoint()]->getY() ) / ( xupright - xlowleft )*width );
952 if ( mHalfEdge[i]->getForced() )
954 p->setPen( mEdgeColor );
960 control[mHalfEdge[i]->getDual()] =
true;
965 double rightborder = width * ( yupright - ylowleft ) / height + xlowleft;
966 for (
unsigned int i = 0; i < mHalfEdge.count() - 1; i++ )
968 if ( mHalfEdge[i]->getPoint() == -1 || mHalfEdge[mHalfEdge[i]->getDual()]->getPoint() == -1 )
975 if ( mHalfEdge[i]->getPoint() != -1 && mHalfEdge[mHalfEdge[i]->getNext()]->getPoint() != -1 && mHalfEdge[mHalfEdge[mHalfEdge[i]->getNext()]->getNext()]->getPoint() != -1 )
977 p1 = mPointVector[mHalfEdge[i]->getPoint()]->getZ();
978 p2 = mPointVector[mHalfEdge[mHalfEdge[i]->getNext()]->getPoint()]->getZ();
979 p3 = mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[i]->getNext()]->getNext()]->getPoint()]->getZ();
980 if ( p1 == p2 && p2 == p3 && halfEdgeBBoxTest( i, xlowleft, ylowleft, rightborder, yupright ) && halfEdgeBBoxTest( mHalfEdge[i]->getNext(), xlowleft, ylowleft, rightborder, yupright ) && halfEdgeBBoxTest( mHalfEdge[mHalfEdge[i]->getNext()]->getNext(), xlowleft, ylowleft, rightborder, yupright ) )
983 pa.setPoint( 0, ( mPointVector[mHalfEdge[i]->getPoint()]->getX() - xlowleft ) / ( yupright - ylowleft )*height, ( yupright - mPointVector[mHalfEdge[i]->getPoint()]->getY() ) / ( yupright - ylowleft )*height );
984 pa.setPoint( 1, ( mPointVector[mHalfEdge[mHalfEdge[i]->getNext()]->getPoint()]->getX() - xlowleft ) / ( yupright - ylowleft )*height, ( yupright - mPointVector[mHalfEdge[mHalfEdge[i]->getNext()]->getPoint()]->getY() ) / ( yupright - ylowleft )*height );
985 pa.setPoint( 2, ( mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[i]->getNext()]->getNext()]->getPoint()]->getX() - xlowleft ) / ( yupright - ylowleft )*height, ( yupright - mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[i]->getNext()]->getNext()]->getPoint()]->getY() ) / ( yupright - ylowleft )*height );
986 QColor
c( 255, 0, 0 );
988 p->drawPolygon( pa );
993 control2[mHalfEdge[i]->getNext()] =
true;
994 control2[mHalfEdge[mHalfEdge[i]->getNext()]->getNext()] =
true;
1002 if ( halfEdgeBBoxTest( i, xlowleft, ylowleft, rightborder, yupright ) )
1004 if ( mHalfEdge[i]->getBreak() )
1006 p->setPen( mBreakEdgeColor );
1008 else if ( mHalfEdge[i]->getForced() )
1010 p->setPen( mForcedEdgeColor );
1013 p->drawLine( ( mPointVector[mHalfEdge[i]->getPoint()]->getX() - xlowleft ) / ( yupright - ylowleft )*height, ( yupright - mPointVector[mHalfEdge[i]->getPoint()]->getY() ) / ( yupright - ylowleft )*height, ( mPointVector[mHalfEdge[mHalfEdge[i]->getDual()]->getPoint()]->getX() - xlowleft ) / ( yupright - ylowleft )*height, ( yupright - mPointVector[mHalfEdge[mHalfEdge[i]->getDual()]->getPoint()]->getY() ) / ( yupright - ylowleft )*height );
1015 if ( mHalfEdge[i]->getForced() )
1017 p->setPen( mEdgeColor );
1022 control[mHalfEdge[i]->getDual()] =
true;
1035 const int firstedge = baseEdgeOfPoint( p2 );
1039 int nextnextedge = firstedge;
1043 edge = mHalfEdge[nextnextedge]->getDual();
1044 if ( mHalfEdge[edge]->getPoint() == p1 )
1046 theedge = nextnextedge;
1049 nextedge = mHalfEdge[edge]->getNext();
1050 nextnextedge = mHalfEdge[nextedge]->getNext();
1052 while ( nextnextedge != firstedge );
1054 if ( theedge == -10 )
1061 return mHalfEdge[mHalfEdge[mHalfEdge[theedge]->getDual()]->getNext()]->getPoint();
1067 const int firstedge = baseEdgeOfPoint( pointno );
1070 if ( firstedge == -1 )
1075 int actedge = firstedge;
1076 int edge, nextedge, nextnextedge;
1079 edge = mHalfEdge[actedge]->getDual();
1080 vlist.append( mHalfEdge[edge]->getPoint() );
1081 nextedge = mHalfEdge[edge]->getNext();
1082 vlist.append( mHalfEdge[nextedge]->getPoint() );
1083 nextnextedge = mHalfEdge[nextedge]->getNext();
1084 vlist.append( mHalfEdge[nextnextedge]->getPoint() );
1085 if ( mHalfEdge[nextnextedge]->getBreak() )
1087 vlist.append( -10 );
1091 vlist.append( -20 );
1093 actedge = nextnextedge;
1095 while ( nextnextedge != firstedge );
1103 if ( mPointVector.size() < 3 )
1109 const int edge = baseEdgeOfTriangle(
point );
1115 else if ( edge >= 0 )
1117 const int ptnr1 = mHalfEdge[edge]->getPoint();
1118 const int ptnr2 = mHalfEdge[mHalfEdge[edge]->getNext()]->getPoint();
1119 const int ptnr3 = mHalfEdge[mHalfEdge[mHalfEdge[edge]->getNext()]->getNext()]->getPoint();
1120 p1.
setX( mPointVector[ptnr1]->x() );
1121 p1.
setY( mPointVector[ptnr1]->y() );
1122 p1.
setZ( mPointVector[ptnr1]->z() );
1123 p2.
setX( mPointVector[ptnr2]->x() );
1124 p2.
setY( mPointVector[ptnr2]->y() );
1125 p2.
setZ( mPointVector[ptnr2]->z() );
1126 p3.
setX( mPointVector[ptnr3]->x() );
1127 p3.
setY( mPointVector[ptnr3]->y() );
1128 p3.
setZ( mPointVector[ptnr3]->z() );
1134 else if ( edge == -20 )
1136 const int ptnr1 = mHalfEdge[mEdgeWithPoint]->getPoint();
1137 const int ptnr2 = mHalfEdge[mHalfEdge[mEdgeWithPoint]->getNext()]->getPoint();
1138 const int ptnr3 = mHalfEdge[mHalfEdge[mHalfEdge[mEdgeWithPoint]->getNext()]->getNext()]->getPoint();
1139 if ( ptnr1 == -1 || ptnr2 == -1 || ptnr3 == -1 )
1143 p1.
setX( mPointVector[ptnr1]->x() );
1144 p1.
setY( mPointVector[ptnr1]->y() );
1145 p1.
setZ( mPointVector[ptnr1]->z() );
1146 p2.
setX( mPointVector[ptnr2]->x() );
1147 p2.
setY( mPointVector[ptnr2]->y() );
1148 p2.
setZ( mPointVector[ptnr2]->z() );
1149 p3.
setX( mPointVector[ptnr3]->x() );
1150 p3.
setY( mPointVector[ptnr3]->y() );
1151 p3.
setZ( mPointVector[ptnr3]->z() );
1157 else if ( edge == -25 )
1159 const int edge1 = baseEdgeOfPoint( mTwiceInsPoint );
1160 const int edge2 = mHalfEdge[edge1]->getNext();
1161 const int edge3 = mHalfEdge[edge2]->getNext();
1162 const int ptnr1 = mHalfEdge[edge1]->getPoint();
1163 const int ptnr2 = mHalfEdge[edge2]->getPoint();
1164 const int ptnr3 = mHalfEdge[edge3]->getPoint();
1165 p1.
setX( mPointVector[ptnr1]->x() );
1166 p1.
setY( mPointVector[ptnr1]->y() );
1167 p1.
setZ( mPointVector[ptnr1]->z() );
1168 p2.
setX( mPointVector[ptnr2]->x() );
1169 p2.
setY( mPointVector[ptnr2]->y() );
1170 p2.
setZ( mPointVector[ptnr2]->z() );
1171 p3.
setX( mPointVector[ptnr3]->x() );
1172 p3.
setY( mPointVector[ptnr3]->y() );
1173 p3.
setZ( mPointVector[ptnr3]->z() );
1179 else if ( edge == -5 )
1181 const int ptnr1 = mHalfEdge[mUnstableEdge]->getPoint();
1182 const int ptnr2 = mHalfEdge[mHalfEdge[mUnstableEdge]->getNext()]->getPoint();
1183 const int ptnr3 = mHalfEdge[mHalfEdge[mHalfEdge[mUnstableEdge]->getNext()]->getNext()]->getPoint();
1184 if ( ptnr1 == -1 || ptnr2 == -1 || ptnr3 == -1 )
1188 p1.
setX( mPointVector[ptnr1]->x() );
1189 p1.
setY( mPointVector[ptnr1]->y() );
1190 p1.
setZ( mPointVector[ptnr1]->z() );
1191 p2.
setX( mPointVector[ptnr2]->x() );
1192 p2.
setY( mPointVector[ptnr2]->y() );
1193 p2.
setZ( mPointVector[ptnr2]->z() );
1194 p3.
setX( mPointVector[ptnr3]->x() );
1195 p3.
setY( mPointVector[ptnr3]->y() );
1196 p3.
setZ( mPointVector[ptnr3]->z() );
1211 if ( mPointVector.size() < 3 )
1217 const int edge = baseEdgeOfTriangle(
point );
1222 else if ( edge >= 0 )
1224 const int ptnr1 = mHalfEdge[edge]->getPoint();
1225 const int ptnr2 = mHalfEdge[mHalfEdge[edge]->getNext()]->getPoint();
1226 const int ptnr3 = mHalfEdge[mHalfEdge[mHalfEdge[edge]->getNext()]->getNext()]->getPoint();
1227 p1.
setX( mPointVector[ptnr1]->x() );
1228 p1.
setY( mPointVector[ptnr1]->y() );
1229 p1.
setZ( mPointVector[ptnr1]->z() );
1230 p2.
setX( mPointVector[ptnr2]->x() );
1231 p2.
setY( mPointVector[ptnr2]->y() );
1232 p2.
setZ( mPointVector[ptnr2]->z() );
1233 p3.
setX( mPointVector[ptnr3]->x() );
1234 p3.
setY( mPointVector[ptnr3]->y() );
1235 p3.
setZ( mPointVector[ptnr3]->z() );
1238 else if ( edge == -20 )
1240 const int ptnr1 = mHalfEdge[mEdgeWithPoint]->getPoint();
1241 const int ptnr2 = mHalfEdge[mHalfEdge[mEdgeWithPoint]->getNext()]->getPoint();
1242 const int ptnr3 = mHalfEdge[mHalfEdge[mHalfEdge[mEdgeWithPoint]->getNext()]->getNext()]->getPoint();
1243 if ( ptnr1 == -1 || ptnr2 == -1 || ptnr3 == -1 )
1247 p1.
setX( mPointVector[ptnr1]->x() );
1248 p1.
setY( mPointVector[ptnr1]->y() );
1249 p1.
setZ( mPointVector[ptnr1]->z() );
1250 p2.
setX( mPointVector[ptnr2]->x() );
1251 p2.
setY( mPointVector[ptnr2]->y() );
1252 p2.
setZ( mPointVector[ptnr2]->z() );
1253 p3.
setX( mPointVector[ptnr3]->x() );
1254 p3.
setY( mPointVector[ptnr3]->y() );
1255 p3.
setZ( mPointVector[ptnr3]->z() );
1258 else if ( edge == -25 )
1260 const int edge1 = baseEdgeOfPoint( mTwiceInsPoint );
1261 const int edge2 = mHalfEdge[edge1]->getNext();
1262 const int edge3 = mHalfEdge[edge2]->getNext();
1263 const int ptnr1 = mHalfEdge[edge1]->getPoint();
1264 const int ptnr2 = mHalfEdge[edge2]->getPoint();
1265 const int ptnr3 = mHalfEdge[edge3]->getPoint();
1266 if ( ptnr1 == -1 || ptnr2 == -1 || ptnr3 == -1 )
1270 p1.
setX( mPointVector[ptnr1]->x() );
1271 p1.
setY( mPointVector[ptnr1]->y() );
1272 p1.
setZ( mPointVector[ptnr1]->z() );
1273 p2.
setX( mPointVector[ptnr2]->x() );
1274 p2.
setY( mPointVector[ptnr2]->y() );
1275 p2.
setZ( mPointVector[ptnr2]->z() );
1276 p3.
setX( mPointVector[ptnr3]->x() );
1277 p3.
setY( mPointVector[ptnr3]->y() );
1278 p3.
setZ( mPointVector[ptnr3]->z() );
1281 else if ( edge == -5 )
1283 const int ptnr1 = mHalfEdge[mUnstableEdge]->getPoint();
1284 const int ptnr2 = mHalfEdge[mHalfEdge[mUnstableEdge]->getNext()]->getPoint();
1285 const int ptnr3 = mHalfEdge[mHalfEdge[mHalfEdge[mUnstableEdge]->getNext()]->getNext()]->getPoint();
1286 if ( ptnr1 == -1 || ptnr2 == -1 || ptnr3 == -1 )
1290 p1.
setX( mPointVector[ptnr1]->x() );
1291 p1.
setY( mPointVector[ptnr1]->y() );
1292 p1.
setZ( mPointVector[ptnr1]->z() );
1293 p2.
setX( mPointVector[ptnr2]->x() );
1294 p2.
setY( mPointVector[ptnr2]->y() );
1295 p2.
setZ( mPointVector[ptnr2]->z() );
1296 p3.
setX( mPointVector[ptnr3]->x() );
1297 p3.
setY( mPointVector[ptnr3]->y() );
1298 p3.
setZ( mPointVector[ptnr3]->z() );
1307unsigned int QgsDualEdgeTriangulation::insertEdge(
int dual,
int next,
int point,
bool mbreak,
bool forced )
1310 mHalfEdge.append( edge );
1311 return mHalfEdge.count() - 1;
1315static bool altitudeTriangleIsSmall(
const QgsPoint &pointBase1,
const QgsPoint &pointBase2,
const QgsPoint &pt3,
double tolerance )
1318 const double x1 = pointBase1.
x();
1319 const double y1 = pointBase1.
y();
1320 const double x2 = pointBase2.
x();
1321 const double y2 = pointBase2.
y();
1322 const double X = pt3.
x();
1323 const double Y = pt3.
y();
1332 t = ( X * ny - Y * nx - x1 * ny + y1 * nx ) / ( ( x2 - x1 ) * ny - ( y2 - y1 ) * nx );
1333 projectedPoint.
setX( x1 + t * ( x2 - x1 ) );
1334 projectedPoint.
setY( y1 + t * ( y2 - y1 ) );
1336 return pt3.
distance( projectedPoint ) < tolerance;
1346 QgsPoint *point1 = mPointVector.at( p1 );
1347 QgsPoint *point2 = mPointVector.at( p2 );
1350 QList<int> crossedEdges;
1353 int pointingEdge = 0;
1355 pointingEdge = baseEdgeOfPoint( p1 );
1357 if ( pointingEdge == -1 )
1363 int actEdge = mHalfEdge[pointingEdge]->getDual();
1364 const int firstActEdge = actEdge;
1371 if ( control > 100 )
1377 if ( mHalfEdge[actEdge]->getPoint() == -1 )
1379 actEdge = mHalfEdge[mHalfEdge[mHalfEdge[actEdge]->getNext()]->getNext()]->getDual();
1384 if ( mHalfEdge[actEdge]->getPoint() == p2 )
1386 mHalfEdge[actEdge]->setForced(
true );
1388 mHalfEdge[mHalfEdge[actEdge]->getDual()]->setForced(
true );
1394 if ( ( point2->
y() - point1->
y() ) / ( mPointVector[mHalfEdge[actEdge]->getPoint()]->y() - point1->
y() ) == ( point2->
x() - point1->
x() ) / ( mPointVector[mHalfEdge[actEdge]->getPoint()]->x() - point1->
x() )
1395 && ( ( point2->
y() - point1->
y() ) >= 0 ) == ( ( mPointVector[mHalfEdge[actEdge]->getPoint()]->y() - point1->
y() ) > 0 )
1396 && ( ( point2->
x() - point1->
x() ) >= 0 ) == ( ( mPointVector[mHalfEdge[actEdge]->getPoint()]->x() - point1->
x() ) > 0 ) )
1399 mHalfEdge[actEdge]->setForced(
true );
1401 mHalfEdge[mHalfEdge[actEdge]->getDual()]->setForced(
true );
1403 const int a = insertForcedSegment( mHalfEdge[actEdge]->getPoint(), p2, segmentType );
1408 const int oppositeEdge = mHalfEdge[actEdge]->getNext();
1410 if ( mHalfEdge[oppositeEdge]->getPoint() == -1 || mHalfEdge[mHalfEdge[oppositeEdge]->getDual()]->getPoint() == -1 )
1412 actEdge = mHalfEdge[mHalfEdge[oppositeEdge]->getNext()]->getDual();
1416 QgsPoint *oppositePoint1 = mPointVector[mHalfEdge[oppositeEdge]->getPoint()];
1417 QgsPoint *oppositePoint2 = mPointVector[mHalfEdge[mHalfEdge[oppositeEdge]->getDual()]->getPoint()];
1419 if ( altitudeTriangleIsSmall( *oppositePoint1, *oppositePoint2, *point1, oppositePoint1->
distance( *oppositePoint2 ) / 500 ) )
1427 mPointVector[mHalfEdge[oppositeEdge]->getPoint()],
1428 mPointVector[mHalfEdge[mHalfEdge[oppositeEdge]->getDual()]->getPoint()] ) )
1434 p3 = mHalfEdge[mHalfEdge[actEdge]->getNext()]->getPoint();
1435 p4 = mHalfEdge[mHalfEdge[mHalfEdge[actEdge]->getNext()]->getDual()]->getPoint();
1437 const double dista = crosspoint.distance( *mPointVector[p3] );
1438 const double distb = crosspoint.distance( *mPointVector[p4] );
1439 if ( dista <= distb )
1441 insertForcedSegment( p1, p3, segmentType );
1442 const int e = insertForcedSegment( p3, p2, segmentType );
1445 else if ( distb <= dista )
1447 insertForcedSegment( p1, p4, segmentType );
1448 const int e = insertForcedSegment( p4, p2, segmentType );
1457 p3 = mHalfEdge[mHalfEdge[actEdge]->getNext()]->getPoint();
1458 p4 = mHalfEdge[mHalfEdge[mHalfEdge[actEdge]->getNext()]->getDual()]->getPoint();
1460 const double distpart = crosspoint.distance( *mPointVector[p4] );
1461 const double disttot = mPointVector[p3]->distance( *mPointVector[p4] );
1462 const float frac = distpart / disttot;
1464 if ( frac == 0 || frac == 1 )
1469 mHalfEdge[actEdge]->setForced(
true );
1471 mHalfEdge[mHalfEdge[actEdge]->getDual()]->setForced(
true );
1473 const int a = insertForcedSegment( p4, p2, segmentType );
1476 else if ( frac == 1 )
1479 mHalfEdge[actEdge]->setForced(
true );
1481 mHalfEdge[mHalfEdge[actEdge]->getDual()]->setForced(
true );
1485 const int a = insertForcedSegment( p3, p2, segmentType );
1497 const int newpoint = splitHalfEdge( mHalfEdge[actEdge]->getNext(), frac );
1498 insertForcedSegment( p1, newpoint, segmentType );
1499 const int e = insertForcedSegment( newpoint, p2, segmentType );
1505 crossedEdges.append( oppositeEdge );
1508 actEdge = mHalfEdge[mHalfEdge[oppositeEdge]->getNext()]->getDual();
1510 while ( actEdge != firstActEdge );
1512 if ( crossedEdges.isEmpty() )
1514 int lastEdgeOppositePointIndex = mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getPoint();
1517 while ( lastEdgeOppositePointIndex != p2 )
1519 QgsPoint *lastEdgePoint1 = mPointVector[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getPoint()];
1520 QgsPoint *lastEdgePoint2 = mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getNext()]->getPoint()];
1521 QgsPoint *lastEdgeOppositePoint = mPointVector[lastEdgeOppositePointIndex];
1530 p3 = mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getPoint();
1531 p4 = mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getPoint();
1533 const double dista = crosspoint.distance( *mPointVector[p3] );
1534 const double distb = crosspoint.distance( *mPointVector[p4] );
1535 if ( dista <= distb )
1537 insertForcedSegment( p1, p3, segmentType );
1538 const int e = insertForcedSegment( p3, p2, segmentType );
1541 else if ( distb <= dista )
1543 insertForcedSegment( p1, p4, segmentType );
1544 const int e = insertForcedSegment( p4, p2, segmentType );
1548 else if ( mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getForced() && mForcedCrossBehavior ==
QgsTriangulation::InsertVertex )
1552 p3 = mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getPoint();
1553 p4 = mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getPoint();
1555 const double distpart = crosspoint.distance( *mPointVector[p3] );
1556 const double disttot = mPointVector[p3]->distance( *mPointVector[p4] );
1557 const float frac = distpart / disttot;
1558 if ( frac == 0 || frac == 1 )
1562 const int newpoint = splitHalfEdge( mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext(), frac );
1563 insertForcedSegment( p1, newpoint, segmentType );
1564 const int e = insertForcedSegment( newpoint, p2, segmentType );
1568 crossedEdges.append( mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext() );
1573 if ( mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getNext()]->getForced() && mForcedCrossBehavior ==
QgsTriangulation::SnappingTypeVertex )
1577 p3 = mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getPoint();
1578 p4 = mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getNext()]->getPoint();
1580 const double dista = crosspoint.distance( *mPointVector[p3] );
1581 const double distb = crosspoint.distance( *mPointVector[p4] );
1582 if ( dista <= distb )
1584 insertForcedSegment( p1, p3, segmentType );
1585 const int e = insertForcedSegment( p3, p2, segmentType );
1588 else if ( distb <= dista )
1590 insertForcedSegment( p1, p4, segmentType );
1591 const int e = insertForcedSegment( p4, p2, segmentType );
1595 else if ( mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getNext()]->getForced() && mForcedCrossBehavior ==
QgsTriangulation::InsertVertex )
1599 p3 = mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getPoint();
1600 p4 = mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getNext()]->getPoint();
1602 const double distpart = crosspoint.distance( *mPointVector[p3] );
1603 const double disttot = mPointVector[p3]->distance( *mPointVector[p4] );
1604 const float frac = distpart / disttot;
1605 if ( frac == 0 || frac == 1 )
1609 const int newpoint = splitHalfEdge( mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getNext(), frac );
1610 insertForcedSegment( p1, newpoint, segmentType );
1611 const int e = insertForcedSegment( newpoint, p2, segmentType );
1615 crossedEdges.append( mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getNext() );
1622 lastEdgeOppositePointIndex = mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getPoint();
1626 QgsPoint *lastEdgePoint1 = mPointVector[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getPoint()];
1627 QgsPoint *lastEdgePoint2 = mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getNext()]->getPoint()];
1628 QgsPoint *lastEdgeOppositePoint = mPointVector[lastEdgeOppositePointIndex];
1629 if ( altitudeTriangleIsSmall( *lastEdgePoint1, *lastEdgePoint2, *lastEdgeOppositePoint, lastEdgePoint1->
distance( *lastEdgePoint2 ) / 500 ) )
1633 QList<int>::const_iterator iter;
1634 for ( iter = crossedEdges.constBegin(); iter != crossedEdges.constEnd(); ++iter )
1636 mHalfEdge[( *( iter ) )]->setForced(
false );
1637 mHalfEdge[( *( iter ) )]->setBreak(
false );
1638 mHalfEdge[mHalfEdge[( *( iter ) )]->getDual()]->setForced(
false );
1639 mHalfEdge[mHalfEdge[( *( iter ) )]->getDual()]->setBreak(
false );
1644 QList<int> freelist = crossedEdges;
1647 QList<int> leftPolygon;
1648 QList<int> rightPolygon;
1651 const int firstedge = freelist.first();
1652 mHalfEdge[firstedge]->setForced(
true );
1654 leftPolygon.append( firstedge );
1655 const int dualfirstedge = mHalfEdge[freelist.first()]->getDual();
1656 mHalfEdge[dualfirstedge]->setForced(
true );
1658 rightPolygon.append( dualfirstedge );
1659 freelist.pop_front();
1663 QList<int>::const_iterator leftiter;
1664 leftiter = crossedEdges.constEnd();
1668 const int newpoint = mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[( *leftiter )]->getDual()]->getNext()]->getNext()]->getPoint();
1669 if ( newpoint != actpointl )
1672 actpointl = newpoint;
1673 const int theedge = mHalfEdge[mHalfEdge[mHalfEdge[( *leftiter )]->getDual()]->getNext()]->getNext();
1674 leftPolygon.append( theedge );
1676 if ( leftiter == crossedEdges.constBegin() )
1682 leftPolygon.append( mHalfEdge[crossedEdges.first()]->getNext() );
1685 QList<int>::const_iterator rightiter;
1687 for ( rightiter = crossedEdges.constBegin(); rightiter != crossedEdges.constEnd(); ++rightiter )
1689 const int newpoint = mHalfEdge[mHalfEdge[mHalfEdge[( *rightiter )]->getNext()]->getNext()]->getPoint();
1690 if ( newpoint != actpointr )
1693 actpointr = newpoint;
1694 const int theedge = mHalfEdge[mHalfEdge[( *rightiter )]->getNext()]->getNext();
1695 rightPolygon.append( theedge );
1701 rightPolygon.append( mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext() );
1702 mHalfEdge[rightPolygon.last()]->setNext( dualfirstedge );
1705 int actedgel = leftPolygon[1];
1706 leftiter = leftPolygon.constBegin();
1708 for ( ; leftiter != leftPolygon.constEnd(); ++leftiter )
1710 mHalfEdge[actedgel]->setNext( ( *leftiter ) );
1711 actedgel = ( *leftiter );
1715 int actedger = rightPolygon[1];
1716 rightiter = rightPolygon.constBegin();
1718 for ( ; rightiter != rightPolygon.constEnd(); ++rightiter )
1720 mHalfEdge[actedger]->setNext( ( *rightiter ) );
1721 actedger = ( *( rightiter ) );
1726 mHalfEdge[leftPolygon.first()]->setNext( ( *( ++( leftiter = leftPolygon.constBegin() ) ) ) );
1727 mHalfEdge[leftPolygon.first()]->setPoint( p2 );
1728 mHalfEdge[leftPolygon.last()]->setNext( firstedge );
1729 mHalfEdge[rightPolygon.first()]->setNext( ( *( ++( rightiter = rightPolygon.constBegin() ) ) ) );
1730 mHalfEdge[rightPolygon.first()]->setPoint( p1 );
1731 mHalfEdge[rightPolygon.last()]->setNext( dualfirstedge );
1733 triangulatePolygon( &leftPolygon, &freelist, firstedge );
1734 triangulatePolygon( &rightPolygon, &freelist, dualfirstedge );
1737 for ( iter = crossedEdges.constBegin(); iter != crossedEdges.constEnd(); ++iter )
1739 checkSwapRecursively( ( *( iter ) ), 0 );
1742 return leftPolygon.first();
1747 mForcedCrossBehavior = b;
1752 mTriangleInterpolator = interpolator;
1758 const double minangle = 0;
1762 bool swapped =
false;
1763 bool *control =
new bool[mHalfEdge.count()];
1765 for (
int i = 0; i <= mHalfEdge.count() - 1; i++ )
1771 for (
int i = 0; i <= mHalfEdge.count() - 1; i++ )
1780 e2 = mHalfEdge[e1]->getNext();
1781 e3 = mHalfEdge[e2]->getNext();
1784 p1 = mHalfEdge[e1]->getPoint();
1785 p2 = mHalfEdge[e2]->getPoint();
1786 p3 = mHalfEdge[e3]->getPoint();
1789 if ( p1 == -1 || p2 == -1 || p3 == -1 )
1797 double el1, el2, el3;
1798 el1 = mPointVector[p1]->z();
1799 el2 = mPointVector[p2]->z();
1800 el3 = mPointVector[p3]->z();
1802 if ( el1 == el2 && el2 == el3 )
1805 if ( swapPossible( ( uint )e1 ) && mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[e1]->getDual()]->getNext()]->getPoint()]->z() != el1 && swapMinAngle( e1 ) > minangle )
1807 doOnlySwap( ( uint )e1 );
1810 else if ( swapPossible( ( uint )e2 ) && mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[e2]->getDual()]->getNext()]->getPoint()]->z() != el2 && swapMinAngle( e2 ) > minangle )
1812 doOnlySwap( ( uint )e2 );
1815 else if ( swapPossible( ( uint )e3 ) && mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[e3]->getDual()]->getNext()]->getPoint()]->z() != el3 && swapMinAngle( e3 ) > minangle )
1817 doOnlySwap( ( uint )e3 );
1847 const double mintol = 17;
1850 std::map<int, double> edge_angle;
1851 std::multimap<double, int> angle_edge;
1852 QSet<int> dontexamine;
1861 const int nhalfedges = mHalfEdge.count();
1863 for (
int i = 0; i < nhalfedges - 1; i++ )
1865 const int next = mHalfEdge[i]->getNext();
1866 const int nextnext = mHalfEdge[next]->getNext();
1868 if ( mHalfEdge[next]->getPoint() != -1 && ( mHalfEdge[i]->getForced() || mHalfEdge[mHalfEdge[mHalfEdge[i]->getDual()]->getNext()]->getPoint() == -1 ) )
1870 if ( !( ( mHalfEdge[next]->getForced() || edgeOnConvexHull( next ) ) || ( mHalfEdge[nextnext]->getForced() || edgeOnConvexHull( nextnext ) ) ) )
1873 while (
MathUtils::inDiametral( mPointVector[mHalfEdge[mHalfEdge[i]->getDual()]->getPoint()], mPointVector[mHalfEdge[i]->getPoint()], mPointVector[mHalfEdge[next]->getPoint()] ) )
1876 const int pointno = splitHalfEdge( i, 0.5 );
1888 for (
int i = 0; i < mHalfEdge.count() - 1; i++ )
1890 p1 = mHalfEdge[mHalfEdge[i]->getDual()]->getPoint();
1891 p2 = mHalfEdge[i]->getPoint();
1892 p3 = mHalfEdge[mHalfEdge[i]->getNext()]->getPoint();
1894 if ( p1 == -1 || p2 == -1 || p3 == -1 )
1898 angle =
MathUtils::angle( mPointVector[p1], mPointVector[p2], mPointVector[p3], mPointVector[p2] );
1900 bool twoforcededges;
1903 twoforcededges = ( mHalfEdge[i]->getForced() || edgeOnConvexHull( i ) ) && ( mHalfEdge[mHalfEdge[i]->getNext()]->getForced() || edgeOnConvexHull( mHalfEdge[i]->getNext() ) );
1905 if ( angle < mintol && !twoforcededges )
1907 edge_angle.insert( std::make_pair( i, angle ) );
1908 angle_edge.insert( std::make_pair( angle, i ) );
1913 for ( std::multimap<double, int>::const_iterator it = angle_edge.begin(); it != angle_edge.end(); ++it )
1919 double minangle = 0;
1922 int minedgenextnext;
1926 while ( !edge_angle.empty() )
1928 minangle = angle_edge.begin()->first;
1930 minedge = angle_edge.begin()->second;
1931 minedgenext = mHalfEdge[minedge]->getNext();
1932 minedgenextnext = mHalfEdge[minedgenext]->getNext();
1935 if ( !
MathUtils::circumcenter( mPointVector[mHalfEdge[minedge]->getPoint()], mPointVector[mHalfEdge[minedgenext]->getPoint()], mPointVector[mHalfEdge[minedgenextnext]->getPoint()], &circumcenter ) )
1937 QgsDebugError( QStringLiteral(
"warning, calculation of circumcenter failed" ) );
1939 dontexamine.insert( minedge );
1940 edge_angle.erase( minedge );
1941 std::multimap<double, int>::iterator minedgeiter = angle_edge.find( minangle );
1942 while ( minedgeiter->second != minedge )
1946 angle_edge.erase( minedgeiter );
1950 if ( !
pointInside( circumcenter.x(), circumcenter.y() ) )
1953 QgsDebugError( QStringLiteral(
"put circumcenter %1//%2 on dontexamine list because it is outside the convex hull" )
1954 .arg( circumcenter.x() ).arg( circumcenter.y() ) );
1955 dontexamine.insert( minedge );
1956 edge_angle.erase( minedge );
1957 std::multimap<double, int>::iterator minedgeiter = angle_edge.find( minangle );
1958 while ( minedgeiter->second != minedge )
1962 angle_edge.erase( minedgeiter );
1968 bool encroached =
false;
1971 int numhalfedges = mHalfEdge.count();
1972 for (
int i = 0; i < numhalfedges; i++ )
1974 if ( mHalfEdge[i]->getForced() || edgeOnConvexHull( i ) )
1976 if (
MathUtils::inDiametral( mPointVector[mHalfEdge[i]->getPoint()], mPointVector[mHalfEdge[mHalfEdge[i]->getDual()]->getPoint()], &circumcenter ) )
1981 int pointno = splitHalfEdge( i, 0.5 );
1984 int pointingedge = baseEdgeOfPoint( pointno );
1986 int actedge = pointingedge;
1989 double angle1, angle2, angle3;
1993 ed1 = mHalfEdge[actedge]->getDual();
1994 pt1 = mHalfEdge[ed1]->getPoint();
1995 ed2 = mHalfEdge[ed1]->getNext();
1996 pt2 = mHalfEdge[ed2]->getPoint();
1997 ed3 = mHalfEdge[ed2]->getNext();
1998 pt3 = mHalfEdge[ed3]->getPoint();
2001 if ( pt1 == -1 || pt2 == -1 || pt3 == -1 )
2006 angle1 =
MathUtils::angle( mPointVector[pt3], mPointVector[pt1], mPointVector[pt2], mPointVector[pt1] );
2007 angle2 =
MathUtils::angle( mPointVector[pt1], mPointVector[pt2], mPointVector[pt3], mPointVector[pt2] );
2008 angle3 =
MathUtils::angle( mPointVector[pt2], mPointVector[pt3], mPointVector[pt1], mPointVector[pt3] );
2011 bool twoforcededges1, twoforcededges2, twoforcededges3;
2013 if ( ( mHalfEdge[ed1]->getForced() || edgeOnConvexHull( ed1 ) ) && ( mHalfEdge[ed2]->getForced() || edgeOnConvexHull( ed2 ) ) )
2015 twoforcededges1 =
true;
2019 twoforcededges1 =
false;
2022 if ( ( mHalfEdge[ed2]->getForced() || edgeOnConvexHull( ed2 ) ) && ( mHalfEdge[ed3]->getForced() || edgeOnConvexHull( ed3 ) ) )
2024 twoforcededges2 =
true;
2028 twoforcededges2 =
false;
2031 if ( ( mHalfEdge[ed3]->getForced() || edgeOnConvexHull( ed3 ) ) && ( mHalfEdge[ed1]->getForced() || edgeOnConvexHull( ed1 ) ) )
2033 twoforcededges3 =
true;
2037 twoforcededges3 =
false;
2041 QSet<int>::iterator ed1iter = dontexamine.find( ed1 );
2042 if ( ed1iter != dontexamine.end() )
2045 dontexamine.erase( ed1iter );
2050 std::map<int, double>::iterator tempit1;
2051 tempit1 = edge_angle.find( ed1 );
2052 if ( tempit1 != edge_angle.end() )
2055 double angle = tempit1->second;
2056 edge_angle.erase( ed1 );
2057 std::multimap<double, int>::iterator tempit2 = angle_edge.lower_bound( angle );
2058 while ( tempit2->second != ed1 )
2062 angle_edge.erase( tempit2 );
2066 if ( angle1 < mintol && !twoforcededges1 )
2068 edge_angle.insert( std::make_pair( ed1, angle1 ) );
2069 angle_edge.insert( std::make_pair( angle1, ed1 ) );
2073 QSet<int>::iterator ed2iter = dontexamine.find( ed2 );
2074 if ( ed2iter != dontexamine.end() )
2077 dontexamine.erase( ed2iter );
2082 std::map<int, double>::iterator tempit1;
2083 tempit1 = edge_angle.find( ed2 );
2084 if ( tempit1 != edge_angle.end() )
2087 double angle = tempit1->second;
2088 edge_angle.erase( ed2 );
2089 std::multimap<double, int>::iterator tempit2 = angle_edge.lower_bound( angle );
2090 while ( tempit2->second != ed2 )
2094 angle_edge.erase( tempit2 );
2098 if ( angle2 < mintol && !twoforcededges2 )
2100 edge_angle.insert( std::make_pair( ed2, angle2 ) );
2101 angle_edge.insert( std::make_pair( angle2, ed2 ) );
2105 QSet<int>::iterator ed3iter = dontexamine.find( ed3 );
2106 if ( ed3iter != dontexamine.end() )
2109 dontexamine.erase( ed3iter );
2114 std::map<int, double>::iterator tempit1;
2115 tempit1 = edge_angle.find( ed3 );
2116 if ( tempit1 != edge_angle.end() )
2119 double angle = tempit1->second;
2120 edge_angle.erase( ed3 );
2121 std::multimap<double, int>::iterator tempit2 = angle_edge.lower_bound( angle );
2122 while ( tempit2->second != ed3 )
2126 angle_edge.erase( tempit2 );
2130 if ( angle3 < mintol && !twoforcededges3 )
2132 edge_angle.insert( std::make_pair( ed3, angle3 ) );
2133 angle_edge.insert( std::make_pair( angle3, ed3 ) );
2138 while ( actedge != pointingedge );
2146 QSet<int> influenceedges;
2147 int baseedge = baseEdgeOfTriangle( circumcenter );
2148 if ( baseedge == -5 )
2151 edge_angle.erase( minedge );
2152 std::multimap<double, int>::iterator minedgeiter = angle_edge.find( minangle );
2153 while ( minedgeiter->second != minedge )
2157 angle_edge.erase( minedgeiter );
2160 else if ( baseedge == -25 )
2163 edge_angle.erase( minedge );
2164 std::multimap<double, int>::iterator minedgeiter = angle_edge.find( minangle );
2165 while ( minedgeiter->second != minedge )
2169 angle_edge.erase( minedgeiter );
2172 else if ( baseedge == -20 )
2174 baseedge = mEdgeWithPoint;
2177 evaluateInfluenceRegion( &circumcenter, baseedge, influenceedges );
2178 evaluateInfluenceRegion( &circumcenter, mHalfEdge[baseedge]->getNext(), influenceedges );
2179 evaluateInfluenceRegion( &circumcenter, mHalfEdge[mHalfEdge[baseedge]->getNext()]->getNext(), influenceedges );
2181 for ( QSet<int>::iterator it = influenceedges.begin(); it != influenceedges.end(); ++it )
2183 if ( ( mHalfEdge[*it]->getForced() || edgeOnConvexHull( *it ) ) &&
MathUtils::inDiametral( mPointVector[mHalfEdge[*it]->getPoint()], mPointVector[mHalfEdge[mHalfEdge[*it]->getDual()]->getPoint()], &circumcenter ) )
2187 const int pointno = splitHalfEdge( *it, 0.5 );
2191 const int pointingedge = baseEdgeOfPoint( pointno );
2193 int actedge = pointingedge;
2196 double angle1, angle2, angle3;
2200 ed1 = mHalfEdge[actedge]->getDual();
2201 pt1 = mHalfEdge[ed1]->getPoint();
2202 ed2 = mHalfEdge[ed1]->getNext();
2203 pt2 = mHalfEdge[ed2]->getPoint();
2204 ed3 = mHalfEdge[ed2]->getNext();
2205 pt3 = mHalfEdge[ed3]->getPoint();
2208 if ( pt1 == -1 || pt2 == -1 || pt3 == -1 )
2213 angle1 =
MathUtils::angle( mPointVector[pt3], mPointVector[pt1], mPointVector[pt2], mPointVector[pt1] );
2214 angle2 =
MathUtils::angle( mPointVector[pt1], mPointVector[pt2], mPointVector[pt3], mPointVector[pt2] );
2215 angle3 =
MathUtils::angle( mPointVector[pt2], mPointVector[pt3], mPointVector[pt1], mPointVector[pt3] );
2218 bool twoforcededges1, twoforcededges2, twoforcededges3;
2222 twoforcededges1 = ( mHalfEdge[ed1]->getForced() || edgeOnConvexHull( ed1 ) ) && ( mHalfEdge[ed2]->getForced() || edgeOnConvexHull( ed2 ) );
2224 twoforcededges2 = ( mHalfEdge[ed2]->getForced() || edgeOnConvexHull( ed2 ) ) && ( mHalfEdge[ed3]->getForced() || edgeOnConvexHull( ed3 ) );
2226 twoforcededges3 = ( mHalfEdge[ed3]->getForced() || edgeOnConvexHull( ed3 ) ) && ( mHalfEdge[ed1]->getForced() || edgeOnConvexHull( ed1 ) );
2230 const QSet<int>::iterator ed1iter = dontexamine.find( ed1 );
2231 if ( ed1iter != dontexamine.end() )
2234 dontexamine.erase( ed1iter );
2239 std::map<int, double>::iterator tempit1;
2240 tempit1 = edge_angle.find( ed1 );
2241 if ( tempit1 != edge_angle.end() )
2244 const double angle = tempit1->second;
2245 edge_angle.erase( ed1 );
2246 std::multimap<double, int>::iterator tempit2 = angle_edge.lower_bound( angle );
2247 while ( tempit2->second != ed1 )
2251 angle_edge.erase( tempit2 );
2255 if ( angle1 < mintol && !twoforcededges1 )
2257 edge_angle.insert( std::make_pair( ed1, angle1 ) );
2258 angle_edge.insert( std::make_pair( angle1, ed1 ) );
2262 const QSet<int>::iterator ed2iter = dontexamine.find( ed2 );
2263 if ( ed2iter != dontexamine.end() )
2266 dontexamine.erase( ed2iter );
2271 std::map<int, double>::iterator tempit1;
2272 tempit1 = edge_angle.find( ed2 );
2273 if ( tempit1 != edge_angle.end() )
2276 const double angle = tempit1->second;
2277 edge_angle.erase( ed2 );
2278 std::multimap<double, int>::iterator tempit2 = angle_edge.lower_bound( angle );
2279 while ( tempit2->second != ed2 )
2283 angle_edge.erase( tempit2 );
2287 if ( angle2 < mintol && !twoforcededges2 )
2289 edge_angle.insert( std::make_pair( ed2, angle2 ) );
2290 angle_edge.insert( std::make_pair( angle2, ed2 ) );
2294 const QSet<int>::iterator ed3iter = dontexamine.find( ed3 );
2295 if ( ed3iter != dontexamine.end() )
2298 dontexamine.erase( ed3iter );
2303 std::map<int, double>::iterator tempit1;
2304 tempit1 = edge_angle.find( ed3 );
2305 if ( tempit1 != edge_angle.end() )
2308 const double angle = tempit1->second;
2309 edge_angle.erase( ed3 );
2310 std::multimap<double, int>::iterator tempit2 = angle_edge.lower_bound( angle );
2311 while ( tempit2->second != ed3 )
2315 angle_edge.erase( tempit2 );
2319 if ( angle3 < mintol && !twoforcededges3 )
2321 edge_angle.insert( std::make_pair( ed3, angle3 ) );
2322 angle_edge.insert( std::make_pair( angle3, ed3 ) );
2327 while ( actedge != pointingedge );
2340 calcPoint( circumcenter.x(), circumcenter.y(), p );
2343 if ( pointno == -100 || pointno == mTwiceInsPoint )
2345 if ( pointno == -100 )
2347 QgsDebugError( QStringLiteral(
"put circumcenter %1//%2 on dontexamine list because of numerical instabilities" )
2348 .arg( circumcenter.x() ).arg( circumcenter.y() ) );
2350 else if ( pointno == mTwiceInsPoint )
2352 QgsDebugError( QStringLiteral(
"put circumcenter %1//%2 on dontexamine list because it is already inserted" )
2353 .arg( circumcenter.x() ).arg( circumcenter.y() ) );
2356 for (
int i = 0; i < mPointVector.count(); i++ )
2358 if ( mPointVector[i]->x() == circumcenter.x() && mPointVector[i]->y() == circumcenter.y() )
2365 QgsDebugError( QStringLiteral(
"point is not present in the triangulation" ) );
2369 dontexamine.insert( minedge );
2370 edge_angle.erase( minedge );
2371 std::multimap<double, int>::iterator minedgeiter = angle_edge.lower_bound( minangle );
2372 while ( minedgeiter->second != minedge )
2376 angle_edge.erase( minedgeiter );
2385 const int pointingedge = baseEdgeOfPoint( pointno );
2387 int actedge = pointingedge;
2390 double angle1, angle2, angle3;
2394 ed1 = mHalfEdge[actedge]->getDual();
2395 pt1 = mHalfEdge[ed1]->getPoint();
2396 ed2 = mHalfEdge[ed1]->getNext();
2397 pt2 = mHalfEdge[ed2]->getPoint();
2398 ed3 = mHalfEdge[ed2]->getNext();
2399 pt3 = mHalfEdge[ed3]->getPoint();
2402 if ( pt1 == -1 || pt2 == -1 || pt3 == -1 )
2407 angle1 =
MathUtils::angle( mPointVector[pt3], mPointVector[pt1], mPointVector[pt2], mPointVector[pt1] );
2408 angle2 =
MathUtils::angle( mPointVector[pt1], mPointVector[pt2], mPointVector[pt3], mPointVector[pt2] );
2409 angle3 =
MathUtils::angle( mPointVector[pt2], mPointVector[pt3], mPointVector[pt1], mPointVector[pt3] );
2412 bool twoforcededges1, twoforcededges2, twoforcededges3;
2414 twoforcededges1 = ( mHalfEdge[ed1]->getForced() || edgeOnConvexHull( ed1 ) ) && ( mHalfEdge[ed2]->getForced() || edgeOnConvexHull( ed2 ) );
2416 twoforcededges2 = ( mHalfEdge[ed2]->getForced() || edgeOnConvexHull( ed2 ) ) && ( mHalfEdge[ed3]->getForced() || edgeOnConvexHull( ed3 ) );
2418 twoforcededges3 = ( mHalfEdge[ed3]->getForced() || edgeOnConvexHull( ed3 ) ) && ( mHalfEdge[ed1]->getForced() || edgeOnConvexHull( ed1 ) );
2422 const QSet<int>::iterator ed1iter = dontexamine.find( ed1 );
2423 if ( ed1iter != dontexamine.end() )
2426 dontexamine.erase( ed1iter );
2431 std::map<int, double>::iterator tempit1;
2432 tempit1 = edge_angle.find( ed1 );
2433 if ( tempit1 != edge_angle.end() )
2436 const double angle = tempit1->second;
2437 edge_angle.erase( ed1 );
2438 std::multimap<double, int>::iterator tempit2 = angle_edge.lower_bound( angle );
2439 while ( tempit2->second != ed1 )
2443 angle_edge.erase( tempit2 );
2447 if ( angle1 < mintol && !twoforcededges1 )
2449 edge_angle.insert( std::make_pair( ed1, angle1 ) );
2450 angle_edge.insert( std::make_pair( angle1, ed1 ) );
2455 const QSet<int>::iterator ed2iter = dontexamine.find( ed2 );
2456 if ( ed2iter != dontexamine.end() )
2459 dontexamine.erase( ed2iter );
2464 std::map<int, double>::iterator tempit1;
2465 tempit1 = edge_angle.find( ed2 );
2466 if ( tempit1 != edge_angle.end() )
2469 const double angle = tempit1->second;
2470 edge_angle.erase( ed2 );
2471 std::multimap<double, int>::iterator tempit2 = angle_edge.lower_bound( angle );
2472 while ( tempit2->second != ed2 )
2476 angle_edge.erase( tempit2 );
2480 if ( angle2 < mintol && !twoforcededges2 )
2482 edge_angle.insert( std::make_pair( ed2, angle2 ) );
2483 angle_edge.insert( std::make_pair( angle2, ed2 ) );
2487 const QSet<int>::iterator ed3iter = dontexamine.find( ed3 );
2488 if ( ed3iter != dontexamine.end() )
2491 dontexamine.erase( ed3iter );
2496 std::map<int, double>::iterator tempit1;
2497 tempit1 = edge_angle.find( ed3 );
2498 if ( tempit1 != edge_angle.end() )
2501 const double angle = tempit1->second;
2502 edge_angle.erase( ed3 );
2503 std::multimap<double, int>::iterator tempit2 = angle_edge.lower_bound( angle );
2504 while ( tempit2->second != ed3 )
2508 angle_edge.erase( tempit2 );
2512 if ( angle3 < mintol && !twoforcededges3 )
2514 edge_angle.insert( std::make_pair( ed3, angle3 ) );
2515 angle_edge.insert( std::make_pair( angle3, ed3 ) );
2520 while ( actedge != pointingedge );
2526 for ( QSet<int>::iterator it = dontexamine.begin(); it != dontexamine.end(); ++it )
2528 QgsDebugMsgLevel( QStringLiteral(
"edge nr. %1 is in dontexamine" ).arg( *it ), 2 );
2534bool QgsDualEdgeTriangulation::swapPossible(
unsigned int edge )
const
2537 if ( mHalfEdge[edge]->getForced() )
2543 if ( mHalfEdge[edge]->getPoint() == -1 || mHalfEdge[mHalfEdge[edge]->getNext()]->getPoint() == -1 || mHalfEdge[mHalfEdge[mHalfEdge[edge]->getDual()]->getNext()]->getPoint() == -1 || mHalfEdge[mHalfEdge[edge]->getDual()]->getPoint() == -1 )
2548 QgsPoint *pta = mPointVector[mHalfEdge[edge]->getPoint()];
2549 QgsPoint *ptb = mPointVector[mHalfEdge[mHalfEdge[edge]->getNext()]->getPoint()];
2550 QgsPoint *ptc = mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[edge]->getNext()]->getNext()]->getPoint()];
2551 QgsPoint *ptd = mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[edge]->getDual()]->getNext()]->getPoint()];
2571void QgsDualEdgeTriangulation::triangulatePolygon( QList<int> *poly, QList<int> *free,
int mainedge )
2575 if ( poly->count() == 3 )
2581 QList<int>::const_iterator iterator = ++( poly->constBegin() );
2582 double distance =
MathUtils::distPointFromLine( mPointVector[mHalfEdge[( *iterator )]->getPoint()], mPointVector[mHalfEdge[mHalfEdge[mainedge]->getDual()]->getPoint()], mPointVector[mHalfEdge[mainedge]->getPoint()] );
2583 int distedge = ( *iterator );
2584 int nextdistedge = mHalfEdge[( *iterator )]->getNext();
2587 while ( iterator != --( poly->constEnd() ) )
2589 if (
MathUtils::distPointFromLine( mPointVector[mHalfEdge[( *iterator )]->getPoint()], mPointVector[mHalfEdge[mHalfEdge[mainedge]->getDual()]->getPoint()], mPointVector[mHalfEdge[mainedge]->getPoint()] ) < distance )
2591 distedge = ( *iterator );
2592 nextdistedge = mHalfEdge[( *iterator )]->getNext();
2593 distance =
MathUtils::distPointFromLine( mPointVector[mHalfEdge[( *iterator )]->getPoint()], mPointVector[mHalfEdge[mHalfEdge[mainedge]->getDual()]->getPoint()], mPointVector[mHalfEdge[mainedge]->getPoint()] );
2598 if ( nextdistedge == ( *( --poly->end() ) ) )
2600 const int inserta = free->first();
2601 const int insertb = mHalfEdge[inserta]->getDual();
2604 mHalfEdge[inserta]->setNext( ( poly->at( 1 ) ) );
2605 mHalfEdge[inserta]->setPoint( mHalfEdge[mainedge]->getPoint() );
2606 mHalfEdge[insertb]->setNext( nextdistedge );
2607 mHalfEdge[insertb]->setPoint( mHalfEdge[distedge]->getPoint() );
2608 mHalfEdge[distedge]->setNext( inserta );
2609 mHalfEdge[mainedge]->setNext( insertb );
2612 for ( iterator = ( ++( poly->constBegin() ) ); ( *iterator ) != nextdistedge; ++iterator )
2614 polya.append( ( *iterator ) );
2616 polya.prepend( inserta );
2620 for ( iterator = polya.begin(); iterator != polya.end(); ++iterator )
2626 triangulatePolygon( &polya, free, inserta );
2629 else if ( distedge == ( *( ++poly->begin() ) ) )
2631 const int inserta = free->first();
2632 const int insertb = mHalfEdge[inserta]->getDual();
2635 mHalfEdge[inserta]->setNext( ( poly->at( 2 ) ) );
2636 mHalfEdge[inserta]->setPoint( mHalfEdge[distedge]->getPoint() );
2637 mHalfEdge[insertb]->setNext( mainedge );
2638 mHalfEdge[insertb]->setPoint( mHalfEdge[mHalfEdge[mainedge]->getDual()]->getPoint() );
2639 mHalfEdge[distedge]->setNext( insertb );
2640 mHalfEdge[( *( --poly->end() ) )]->setNext( inserta );
2643 iterator = poly->constBegin();
2645 while ( iterator != poly->constEnd() )
2647 polya.append( ( *iterator ) );
2650 polya.prepend( inserta );
2652 triangulatePolygon( &polya, free, inserta );
2657 const int inserta = free->first();
2658 const int insertb = mHalfEdge[inserta]->getDual();
2661 const int insertc = free->first();
2662 const int insertd = mHalfEdge[insertc]->getDual();
2665 mHalfEdge[inserta]->setNext( ( poly->at( 1 ) ) );
2666 mHalfEdge[inserta]->setPoint( mHalfEdge[mainedge]->getPoint() );
2667 mHalfEdge[insertb]->setNext( insertd );
2668 mHalfEdge[insertb]->setPoint( mHalfEdge[distedge]->getPoint() );
2669 mHalfEdge[insertc]->setNext( nextdistedge );
2670 mHalfEdge[insertc]->setPoint( mHalfEdge[distedge]->getPoint() );
2671 mHalfEdge[insertd]->setNext( mainedge );
2672 mHalfEdge[insertd]->setPoint( mHalfEdge[mHalfEdge[mainedge]->getDual()]->getPoint() );
2674 mHalfEdge[distedge]->setNext( inserta );
2675 mHalfEdge[mainedge]->setNext( insertb );
2676 mHalfEdge[( *( --poly->end() ) )]->setNext( insertc );
2682 for ( iterator = ++( poly->constBegin() ); ( *iterator ) != nextdistedge; ++iterator )
2684 polya.append( ( *iterator ) );
2686 polya.prepend( inserta );
2689 while ( iterator != poly->constEnd() )
2691 polyb.append( ( *iterator ) );
2694 polyb.prepend( insertc );
2696 triangulatePolygon( &polya, free, inserta );
2697 triangulatePolygon( &polyb, free, insertc );
2711 unsigned int actedge = mEdgeInside;
2719 if ( runs > MAX_BASE_ITERATIONS )
2721 QgsDebugError( QStringLiteral(
"warning, instability detected: Point coordinates: %1//%2" ).arg( x ).arg( y ) );
2725 if (
MathUtils::leftOf(
point, mPointVector[mHalfEdge[mHalfEdge[actedge]->getDual()]->getPoint()], mPointVector[mHalfEdge[actedge]->getPoint()] ) < ( -
leftOfTresh ) )
2734 else if ( fabs(
MathUtils::leftOf(
point, mPointVector[mHalfEdge[mHalfEdge[actedge]->getDual()]->getPoint()], mPointVector[mHalfEdge[actedge]->getPoint()] ) ) <=
leftOfTresh )
2737 mEdgeWithPoint = actedge;
2747 actedge = mHalfEdge[actedge]->getDual();
2753 actedge = mHalfEdge[actedge]->getNext();
2754 if ( mHalfEdge[actedge]->getPoint() == -1 )
2760 mEdgeOutside = (
unsigned int )mHalfEdge[mHalfEdge[actedge]->getNext()]->getNext();
2770 if ( numinstabs > 0 )
2772 QgsDebugError( QStringLiteral(
"numerical instabilities" ) );
2780 mEdgeInside = actedge;
2785bool DualEdgeTriangulation::readFromTAFF( QString filename )
2787 QApplication::setOverrideCursor( QCursor( Qt::WaitCursor ) );
2789 QFile file( filename );
2790 file.open( IO_Raw | IO_ReadOnly );
2791 QBuffer buffer( file.readAll() );
2794 QTextStream textstream( &buffer );
2795 buffer.open( IO_ReadOnly );
2798 int numberofhalfedges;
2802 while ( buff.mid( 0, 4 ) !=
"TRIA" )
2804 buff = textstream.readLine();
2806 while ( buff.mid( 0, 4 ) !=
"NEDG" )
2808 buff = textstream.readLine();
2810 numberofhalfedges = buff.section(
' ', 1, 1 ).toInt();
2811 mHalfEdge.resize( numberofhalfedges );
2814 while ( buff.mid( 0, 4 ) !=
"DATA" )
2819 int nr1, nr2, dual1, dual2, point1, point2, next1, next2, fo1, fo2, b1, b2;
2820 bool forced1, forced2, break1, break2;
2824 QProgressBar *edgebar =
new QProgressBar();
2825 edgebar->setCaption( tr(
"Reading edges…" ) );
2826 edgebar->setTotalSteps( numberofhalfedges / 2 );
2827 edgebar->setMinimumWidth( 400 );
2828 edgebar->move( 500, 500 );
2831 for (
int i = 0; i < numberofhalfedges / 2; i++ )
2833 if ( i % 1000 == 0 )
2835 edgebar->setProgress( i );
2839 textstream >> point1;
2840 textstream >> next1;
2864 textstream >> point2;
2865 textstream >> next2;
2903 mHalfEdge.insert( nr1, hf1 );
2904 mHalfEdge.insert( nr2, hf2 );
2908 edgebar->setProgress( numberofhalfedges / 2 );
2912 for (
int i = 0; i < numberofhalfedges; i++ )
2915 a = mHalfEdge[i]->getPoint();
2916 b = mHalfEdge[mHalfEdge[i]->getDual()]->getPoint();
2917 c = mHalfEdge[mHalfEdge[i]->getNext()]->getPoint();
2918 d = mHalfEdge[mHalfEdge[mHalfEdge[i]->getDual()]->getNext()]->getPoint();
2919 if ( a != -1 && b != -1 &&
c != -1 && d != -1 )
2927 while ( buff.mid( 0, 4 ) !=
"POIN" )
2929 buff = textstream.readLine();
2932 while ( buff.mid( 0, 4 ) !=
"NPTS" )
2934 buff = textstream.readLine();
2937 numberofpoints = buff.section(
' ', 1, 1 ).toInt();
2938 mPointVector.resize( numberofpoints );
2940 while ( buff.mid( 0, 4 ) !=
"DATA" )
2945 QProgressBar *pointbar =
new QProgressBar();
2946 pointbar->setCaption( tr(
"Reading points…" ) );
2947 pointbar->setTotalSteps( numberofpoints );
2948 pointbar->setMinimumWidth( 400 );
2949 pointbar->move( 500, 500 );
2954 for (
int i = 0; i < numberofpoints; i++ )
2956 if ( i % 1000 == 0 )
2958 pointbar->setProgress( i );
2968 mPointVector.insert( i, p );
2984 else if ( x > xMax )
2993 else if ( y > yMax )
3000 pointbar->setProgress( numberofpoints );
3002 QApplication::restoreOverrideCursor();
3006bool DualEdgeTriangulation::saveToTAFF( QString filename )
const
3008 QFile outputfile( filename );
3009 if ( !outputfile.open( IO_WriteOnly ) )
3011 QMessageBox::warning( 0, tr(
"Warning" ), tr(
"File could not be written." ), QMessageBox::Ok, QMessageBox::NoButton, QMessageBox::NoButton );
3015 QTextStream outstream( &outputfile );
3016 outstream.precision( 9 );
3019 outstream <<
"TRIA" << std::endl << std::flush;
3020 outstream <<
"NEDG " << mHalfEdge.count() << std::endl << std::flush;
3021 outstream <<
"PANO 1" << std::endl << std::flush;
3022 outstream <<
"DATA ";
3024 bool *cont =
new bool[mHalfEdge.count()];
3025 for (
unsigned int i = 0; i <= mHalfEdge.count() - 1; i++ )
3030 for (
unsigned int i = 0; i < mHalfEdge.count(); i++ )
3037 int dual = mHalfEdge[i]->getDual();
3038 outstream << i <<
" " << mHalfEdge[i]->getPoint() <<
" " << mHalfEdge[i]->getNext() <<
" " << mHalfEdge[i]->getForced() <<
" " << mHalfEdge[i]->getBreak() <<
" ";
3039 outstream << dual <<
" " << mHalfEdge[dual]->getPoint() <<
" " << mHalfEdge[dual]->getNext() <<
" " << mHalfEdge[dual]->getForced() <<
" " << mHalfEdge[dual]->getBreak() <<
" ";
3043 outstream << std::endl << std::flush;
3044 outstream << std::endl << std::flush;
3049 outstream <<
"POIN" << std::endl << std::flush;
3050 outstream <<
"NPTS " << getNumberOfPoints() << std::endl << std::flush;
3051 outstream <<
"PATT 3" << std::endl << std::flush;
3052 outstream <<
"DATA ";
3054 for (
int i = 0; i < getNumberOfPoints(); i++ )
3057 outstream << p->getX() <<
" " << p->getY() <<
" " << p->getZ() <<
" ";
3059 outstream << std::endl << std::flush;
3060 outstream << std::endl << std::flush;
3069 const int edge1 = baseEdgeOfTriangle( p );
3076 edge2 = mHalfEdge[edge1]->getNext();
3077 edge3 = mHalfEdge[edge2]->getNext();
3078 point1 =
point( mHalfEdge[edge1]->getPoint() );
3079 point2 =
point( mHalfEdge[edge2]->getPoint() );
3080 point3 =
point( mHalfEdge[edge3]->getPoint() );
3081 if ( point1 && point2 && point3 )
3084 double dist1, dist2, dist3;
3088 if ( dist1 <= dist2 && dist1 <= dist3 )
3091 if ( swapPossible( edge1 ) )
3093 doOnlySwap( edge1 );
3096 else if ( dist2 <= dist1 && dist2 <= dist3 )
3099 if ( swapPossible( edge2 ) )
3101 doOnlySwap( edge2 );
3104 else if ( dist3 <= dist1 && dist3 <= dist2 )
3107 if ( swapPossible( edge3 ) )
3109 doOnlySwap( edge3 );
3131 const int edge1 = baseEdgeOfTriangle( p );
3134 const int edge2 = mHalfEdge[edge1]->getNext();
3135 const int edge3 = mHalfEdge[edge2]->getNext();
3139 if ( point1 && point2 && point3 )
3145 if ( dist1 <= dist2 && dist1 <= dist3 )
3147 p1 = mHalfEdge[edge1]->getPoint();
3148 p2 = mHalfEdge[mHalfEdge[edge1]->getNext()]->getPoint();
3149 p3 = mHalfEdge[mHalfEdge[edge1]->getDual()]->getPoint();
3150 p4 = mHalfEdge[mHalfEdge[mHalfEdge[edge1]->getDual()]->getNext()]->getPoint();
3152 else if ( dist2 <= dist1 && dist2 <= dist3 )
3154 p1 = mHalfEdge[edge2]->getPoint();
3155 p2 = mHalfEdge[mHalfEdge[edge2]->getNext()]->getPoint();
3156 p3 = mHalfEdge[mHalfEdge[edge2]->getDual()]->getPoint();
3157 p4 = mHalfEdge[mHalfEdge[mHalfEdge[edge2]->getDual()]->getNext()]->getPoint();
3161 p1 = mHalfEdge[edge3]->getPoint();
3162 p2 = mHalfEdge[mHalfEdge[edge3]->getNext()]->getPoint();
3163 p3 = mHalfEdge[mHalfEdge[edge3]->getDual()]->getPoint();
3164 p4 = mHalfEdge[mHalfEdge[mHalfEdge[edge3]->getDual()]->getNext()]->getPoint();
3192 bool *alreadyVisitedEdges =
new bool[mHalfEdge.size()];
3193 if ( !alreadyVisitedEdges )
3199 for (
int i = 0; i < mHalfEdge.size(); ++i )
3201 alreadyVisitedEdges[i] =
false;
3204 for (
int i = 0; i < mHalfEdge.size(); ++i )
3209 HalfEdge *currentEdge = mHalfEdge[i];
3210 if ( currentEdge->
getPoint() != -1 && mHalfEdge[currentEdge->
getDual()]->getPoint() != -1 && !alreadyVisitedEdges[currentEdge->
getDual()] )
3216 QgsPoint *p2 = mPointVector[mHalfEdge[currentEdge->
getDual()]->getPoint()];
3224 QString attributeString;
3229 attributeString = QStringLiteral(
"break line" );
3233 attributeString = QStringLiteral(
"structure line" );
3240 alreadyVisitedEdges[i] =
true;
3243 delete [] alreadyVisitedEdges;
3253 QVector< bool> edgeToTreat( mHalfEdge.count(),
true );
3254 QHash<HalfEdge *, int > edgesHash;
3255 for (
int i = 0; i < mHalfEdge.count(); ++i )
3257 edgesHash.insert( mHalfEdge[i], i );
3266 const int edgeCount = edgeToTreat.count();
3267 for (
int i = 0 ; i < edgeCount; ++i )
3269 bool containVirtualPoint =
false;
3270 if ( edgeToTreat[i] )
3272 HalfEdge *currentEdge = mHalfEdge[i];
3277 edgeToTreat[edgesHash.value( currentEdge )] =
false;
3278 face.append( currentEdge->
getPoint() );
3279 containVirtualPoint |= currentEdge->
getPoint() == -1;
3280 currentEdge = mHalfEdge.at( currentEdge->
getNext() );
3282 while ( currentEdge != firstEdge && !containVirtualPoint && ( !feedback || !feedback->
isCanceled() ) );
3283 if ( !containVirtualPoint )
3284 mesh.
faces.append( face );
3288 feedback->
setProgress( ( 100 * i ) / edgeCount ) ;
3295double QgsDualEdgeTriangulation::swapMinAngle(
int edge )
const
3298 QgsPoint *p2 =
point( mHalfEdge[mHalfEdge[edge]->getNext()]->getPoint() );
3299 QgsPoint *p3 =
point( mHalfEdge[mHalfEdge[edge]->getDual()]->getPoint() );
3300 QgsPoint *p4 =
point( mHalfEdge[mHalfEdge[mHalfEdge[edge]->getDual()]->getNext()]->getPoint() );
3307 if ( angle2 < minangle )
3312 if ( angle3 < minangle )
3317 if ( angle4 < minangle )
3322 if ( angle5 < minangle )
3327 if ( angle6 < minangle )
3335int QgsDualEdgeTriangulation::splitHalfEdge(
int edge,
float position )
3338 if ( position < 0 || position > 1 )
3340 QgsDebugError( QStringLiteral(
"warning, position is not between 0 and 1" ) );
3344 QgsPoint *p =
new QgsPoint( mPointVector[mHalfEdge[edge]->getPoint()]->x()*position + mPointVector[mHalfEdge[mHalfEdge[edge]->getDual()]->getPoint()]->x() * ( 1 - position ), mPointVector[mHalfEdge[edge]->getPoint()]->y()*position + mPointVector[mHalfEdge[mHalfEdge[edge]->getDual()]->getPoint()]->y() * ( 1 - position ), 0 );
3349 p->
setZ( zvaluepoint.z() );
3352 if ( mPointVector.count() >= mPointVector.size() )
3354 mPointVector.resize( mPointVector.count() + 1 );
3356 QgsDebugMsgLevel( QStringLiteral(
"inserting point nr. %1, %2//%3//%4" ).arg( mPointVector.count() ).arg( p->
x() ).arg( p->
y() ).arg( p->
z() ), 2 );
3357 mPointVector.insert( mPointVector.count(), p );
3360 const int dualedge = mHalfEdge[edge]->getDual();
3361 const int edge1 = insertEdge( -10, -10, mPointVector.count() - 1,
false,
false );
3362 const int edge2 = insertEdge( edge1, mHalfEdge[mHalfEdge[edge]->getNext()]->getNext(), mHalfEdge[mHalfEdge[edge]->getNext()]->getPoint(),
false,
false );
3363 const int edge3 = insertEdge( -10, mHalfEdge[mHalfEdge[dualedge]->getNext()]->getNext(), mHalfEdge[mHalfEdge[dualedge]->getNext()]->getPoint(),
false,
false );
3364 const int edge4 = insertEdge( edge3, dualedge, mPointVector.count() - 1,
false,
false );
3365 const int edge5 = insertEdge( -10, mHalfEdge[edge]->getNext(), mHalfEdge[edge]->getPoint(), mHalfEdge[edge]->getBreak(), mHalfEdge[edge]->getForced() );
3366 const int edge6 = insertEdge( edge5, edge3, mPointVector.count() - 1, mHalfEdge[dualedge]->getBreak(), mHalfEdge[dualedge]->getForced() );
3367 mHalfEdge[edge1]->setDual( edge2 );
3368 mHalfEdge[edge1]->setNext( edge5 );
3369 mHalfEdge[edge3]->setDual( edge4 );
3370 mHalfEdge[edge5]->setDual( edge6 );
3373 mHalfEdge[mHalfEdge[edge]->getNext()]->setNext( edge1 );
3374 mHalfEdge[mHalfEdge[dualedge]->getNext()]->setNext( edge4 );
3375 mHalfEdge[edge]->setNext( edge2 );
3376 mHalfEdge[edge]->setPoint( mPointVector.count() - 1 );
3377 mHalfEdge[mHalfEdge[edge3]->getNext()]->setNext( edge6 );
3380 checkSwapRecursively( mHalfEdge[edge5]->getNext(), 0 );
3381 checkSwapRecursively( mHalfEdge[edge2]->getNext(), 0 );
3382 checkSwapRecursively( mHalfEdge[dualedge]->getNext(), 0 );
3383 checkSwapRecursively( mHalfEdge[edge3]->getNext(), 0 );
3387 return mPointVector.count() - 1;
3390bool QgsDualEdgeTriangulation::edgeOnConvexHull(
int edge )
3392 return ( mHalfEdge[mHalfEdge[edge]->getNext()]->getPoint() == -1 || mHalfEdge[mHalfEdge[mHalfEdge[edge]->getDual()]->getNext()]->getPoint() == -1 );
3395void QgsDualEdgeTriangulation::evaluateInfluenceRegion(
QgsPoint *point,
int edge, QSet<int> &set )
3397 if ( set.find( edge ) == set.end() )
3406 if ( !mHalfEdge[edge]->getForced() && !edgeOnConvexHull( edge ) )
3409 if (
inCircle( *
point, *mPointVector[mHalfEdge[mHalfEdge[edge]->getDual()]->getPoint()], *mPointVector[mHalfEdge[edge]->getPoint()], *mPointVector[mHalfEdge[mHalfEdge[edge]->getNext()]->getPoint()] ) )
3411 evaluateInfluenceRegion(
point, mHalfEdge[mHalfEdge[edge]->getDual()]->getNext(), set );
3412 evaluateInfluenceRegion(
point, mHalfEdge[mHalfEdge[mHalfEdge[edge]->getDual()]->getNext()]->getNext(), set );
3417int QgsDualEdgeTriangulation::firstEdgeOutSide()
3420 while ( ( mHalfEdge[mHalfEdge[edge]->getNext()]->getPoint() != -1 ||
3421 mHalfEdge[edge]->getPoint() == -1 || mHalfEdge[mHalfEdge[edge]->getDual()]->getPoint() == -1 ) &&
3422 edge < mHalfEdge.count() )
3427 if ( edge >= mHalfEdge.count() )
3433void QgsDualEdgeTriangulation::removeLastPoint()
3435 if ( mPointVector.isEmpty() )
3437 QgsPoint *p = mPointVector.takeLast();
bool getForced() const
Returns, whether the HalfEdge belongs to a constrained edge or not.
int getNext() const
Returns the number of the next HalfEdge.
void setNext(int n)
Sets the number of the next HalfEdge.
void setPoint(int p)
Sets the number of point at which this HalfEdge points.
int getPoint() const
Returns the number of the point at which this HalfEdge points.
int getDual() const
Returns the number of the dual HalfEdge.
void setDual(int d)
Sets the number of the dual HalfEdge.
void setForced(bool f)
Sets the forced flag.
bool getBreak() const
Returns, whether the HalfEdge belongs to a break line or not.
void setBreak(bool b)
Sets the break flag.
bool saveTriangulation(QgsFeatureSink *sink, QgsFeedback *feedback=nullptr) const override
Saves the triangulation features to a feature sink.
QgsPoint * point(int i) const override
Draws the points, edges and the forced lines.
QList< int > surroundingTriangles(int pointno) override
Returns a value list with the information of the triangles surrounding (counterclockwise) a point.
void ruppertRefinement() override
Adds points to make the triangles better shaped (algorithm of ruppert)
bool pointInside(double x, double y) override
Returns true, if the point with coordinates x and y is inside the convex hull and false otherwise.
void setTriangleInterpolator(TriangleInterpolator *interpolator) override
Sets an interpolator object.
void performConsistencyTest() override
Performs a consistency check, remove this later.
int oppositePoint(int p1, int p2) override
Returns the number of the point opposite to the triangle points p1, p2 (which have to be on a halfedg...
int addPoint(const QgsPoint &p) override
Adds a point to the triangulation.
virtual QgsMesh triangulationToMesh(QgsFeedback *feedback=nullptr) const override
Returns a QgsMesh corresponding to the triangulation.
void setForcedCrossBehavior(QgsTriangulation::ForcedCrossBehavior b) override
Sets the behavior of the triangulation in case of crossing forced lines.
bool calcNormal(double x, double y, QgsPoint &result) override
Calculates the normal at a point on the surface.
void addLine(const QVector< QgsPoint > &points, QgsInterpolator::SourceType lineType) override
Adds a line (e.g.
bool calcPoint(double x, double y, QgsPoint &result) override
Calculates x-, y and z-value of the point on the surface and assigns it to 'result'.
void eliminateHorizontalTriangles() override
Eliminates the horizontal triangles by swapping or by insertion of new points.
bool triangleVertices(double x, double y, QgsPoint &p1, int &n1, QgsPoint &p2, int &n2, QgsPoint &p3, int &n3) override
Finds out in which triangle the point with coordinates x and y is and assigns the numbers of the vert...
~QgsDualEdgeTriangulation() override
bool swapEdge(double x, double y) override
Swaps the edge which is closest to the point with x and y coordinates (if this is possible)
QList< int > pointsAroundEdge(double x, double y) override
Returns a value list with the numbers of the four points, which would be affected by an edge swap....
An interface for objects which accept features via addFeature(s) methods.
virtual bool addFeature(QgsFeature &feature, QgsFeatureSink::Flags flags=QgsFeatureSink::Flags())
Adds a single feature to the sink.
@ FastInsert
Use faster inserts, at the cost of updating the passed features to reflect changes made at the provid...
The feature class encapsulates a single feature including its unique ID, geometry and a list of field...
Q_INVOKABLE bool setAttribute(int field, const QVariant &attr)
Sets an attribute's value by field index.
void initAttributes(int fieldCount)
Initialize this feature with the given number of fields.
void setGeometry(const QgsGeometry &geometry)
Set the feature's geometry.
Base class for feedback objects to be used for cancellation of something running in a worker thread.
bool isCanceled() const
Tells whether the operation has been canceled already.
void setProgress(double progress)
Sets the current progress for the feedback object.
static QgsGeometry fromPolylineXY(const QgsPolylineXY &polyline)
Creates a new LineString geometry from a list of QgsPointXY points.
SourceType
Describes the type of input data.
@ SourceBreakLines
Break lines.
A class to represent a 2D point.
Point geometry type, with support for z-dimension and m-values.
void setY(double y)
Sets the point's y-coordinate.
void setX(double x)
Sets the point's x-coordinate.
double distance(double x, double y) const
Returns the Cartesian 2D distance between this point and a specified x, y coordinate.
void setZ(double z)
Sets the point's z-coordinate.
double distanceSquared(double x, double y) const
Returns the Cartesian 2D squared distance between this point a specified x, y coordinate.
ForcedCrossBehavior
Enumeration describing the behavior, if two forced lines cross.
@ SnappingTypeVertex
The second inserted forced line is snapped to the closest vertice of the first inserted forced line.
This is an interface for interpolator classes for triangulations.
virtual bool calcNormVec(double x, double y, QgsPoint &result)=0
Calculates the normal vector and assigns it to vec.
virtual bool calcPoint(double x, double y, QgsPoint &result)=0
Performs a linear interpolation in a triangle and assigns the x-,y- and z-coordinates to point.
bool ANALYSIS_EXPORT inDiametral(QgsPoint *p1, QgsPoint *p2, QgsPoint *point)
Tests, whether 'point' is inside the diametral circle through 'p1' and 'p2'.
double ANALYSIS_EXPORT angle(QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, QgsPoint *p4)
Calculates the angle between two segments (in 2 dimension, z-values are ignored)
bool ANALYSIS_EXPORT lineIntersection(QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, QgsPoint *p4)
Returns true, if line1 (p1 to p2) and line2 (p3 to p4) intersect. If the lines have an endpoint in co...
bool ANALYSIS_EXPORT circumcenter(QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, QgsPoint *result)
Calculates the center of the circle passing through p1, p2 and p3. Returns true in case of success an...
double ANALYSIS_EXPORT leftOf(const QgsPoint &thepoint, const QgsPoint *p1, const QgsPoint *p2)
Returns whether 'thepoint' is left or right of the line from 'p1' to 'p2'. Negative values mean left ...
double ANALYSIS_EXPORT distPointFromLine(QgsPoint *thepoint, QgsPoint *p1, QgsPoint *p2)
Calculates the (2 dimensional) distance from 'thepoint' to the line defined by p1 and p2.
bool ANALYSIS_EXPORT inCircle(QgsPoint *testp, QgsPoint *p1, QgsPoint *p2, QgsPoint *p3)
Tests, whether 'testp' is inside the circle through 'p1', 'p2' and 'p3'.
As part of the API refactoring and improvements which landed in the Processing API was substantially reworked from the x version This was done in order to allow much of the underlying Processing framework to be ported into c
QVector< QgsPointXY > QgsPolylineXY
Polyline as represented as a vector of two-dimensional points.
#define QgsDebugMsgLevel(str, level)
#define QgsDebugError(str)
QVector< int > QgsMeshFace
List of vertex indexes.
Mesh - vertices, edges and faces.
QVector< QgsMeshVertex > vertices
QVector< QgsMeshFace > faces