1 /***************************************************************************
2 qgsinternalgeometryengine.cpp - QgsInternalGeometryEngine
3
4 ---------------------
5 begin : 13.1.2016
6 copyright : (C) 2016 by Matthias Kuhn
7 email : matthias@opengis.ch
8 ***************************************************************************
9 * *
10 * This program is free software; you can redistribute it and/or modify *
11 * it under the terms of the GNU General Public License as published by *
12 * the Free Software Foundation; either version 2 of the License, or *
13 * (at your option) any later version. *
14 * *
15 ***************************************************************************/
16
17 #include "qgsinternalgeometryengine.h"
18
19 #include "qgslinestring.h"
20 #include "qgsmultipolygon.h"
21 #include "qgspolygon.h"
22 #include "qgsmulticurve.h"
23 #include "qgscircularstring.h"
24 #include "qgsgeometry.h"
25 #include "qgsgeometryutils.h"
26 #include "qgslinesegment.h"
27 #include "qgscircle.h"
28 #include "qgslogger.h"
29 #include "qgstessellator.h"
30 #include "qgsfeedback.h"
31 #include "qgsgeometryengine.h"
32 #include <QTransform>
33 #include <functional>
34 #include <memory>
35 #include <queue>
36 #include <random>
37
QgsInternalGeometryEngine(const QgsGeometry & geometry)38 QgsInternalGeometryEngine::QgsInternalGeometryEngine( const QgsGeometry &geometry )
39 : mGeometry( geometry.constGet() )
40 {
41
42 }
43
lastError() const44 QString QgsInternalGeometryEngine::lastError() const
45 {
46 return mLastError;
47 }
48
49
50 enum class Direction
51 {
52 Up,
53 Right,
54 Down,
55 Left,
56 None
57 };
58
59 /**
60 * Determines the direction of an edge from p1 to p2. maxDev is the tangent of
61 * the maximum allowed edge deviation angle. If the edge deviates more than
62 * the allowed angle, Direction::None will be returned.
63 */
getEdgeDirection(const QgsPoint & p1,const QgsPoint & p2,double maxDev)64 Direction getEdgeDirection( const QgsPoint &p1, const QgsPoint &p2, double maxDev )
65 {
66 double dx = p2.x() - p1.x();
67 double dy = p2.y() - p1.y();
68 if ( ( dx == 0.0 ) && ( dy == 0.0 ) )
69 return Direction::None;
70 if ( fabs( dx ) >= fabs( dy ) )
71 {
72 double dev = fabs( dy ) / fabs( dx );
73 if ( dev > maxDev )
74 return Direction::None;
75 return dx > 0.0 ? Direction::Right : Direction::Left;
76 }
77 else
78 {
79 double dev = fabs( dx ) / fabs( dy );
80 if ( dev > maxDev )
81 return Direction::None;
82 return dy > 0.0 ? Direction::Up : Direction::Down;
83 }
84 }
85
86 /**
87 * Checks whether the polygon consists of four nearly axis-parallel sides. All
88 * consecutive edges having the same direction are considered to belong to the
89 * same side.
90 */
getEdgeDirections(const QgsPolygon * g,double maxDev)91 std::pair<bool, std::array<Direction, 4>> getEdgeDirections( const QgsPolygon *g, double maxDev )
92 {
93 std::pair<bool, std::array<Direction, 4>> ret = { false, { Direction::None, Direction::None, Direction::None, Direction::None } };
94 // The polygon might start in the middle of a side. Hence, we need a fifth
95 // direction to record the beginning of the side when we went around the
96 // polygon.
97 std::array<Direction, 5> dirs;
98
99 int idx = 0;
100 QgsAbstractGeometry::vertex_iterator previous = g->vertices_begin();
101 QgsAbstractGeometry::vertex_iterator current = previous;
102 ++current;
103 QgsAbstractGeometry::vertex_iterator end = g->vertices_end();
104 while ( current != end )
105 {
106 Direction dir = getEdgeDirection( *previous, *current, maxDev );
107 if ( dir == Direction::None )
108 return ret;
109 if ( idx == 0 )
110 {
111 dirs[0] = dir;
112 ++idx;
113 }
114 else if ( dir != dirs[idx - 1] )
115 {
116 if ( idx == 5 )
117 return ret;
118 dirs[idx] = dir;
119 ++idx;
120 }
121 previous = current;
122 ++current;
123 }
124 ret.first = ( idx == 5 ) ? ( dirs[0] == dirs[4] ) : ( idx == 4 );
125 std::copy( dirs.begin(), dirs.begin() + 4, ret.second.begin() );
126 return ret;
127 }
128
matchesOrientation(std::array<Direction,4> dirs,std::array<Direction,4> oriented)129 bool matchesOrientation( std::array<Direction, 4> dirs, std::array<Direction, 4> oriented )
130 {
131 int idx = std::find( oriented.begin(), oriented.end(), dirs[0] ) - oriented.begin();
132 for ( int i = 1; i < 4; ++i )
133 {
134 if ( dirs[i] != oriented[( idx + i ) % 4] )
135 return false;
136 }
137 return true;
138 }
139
140 /**
141 * Checks whether the 4 directions in dirs make up a clockwise rectangle.
142 */
isClockwise(std::array<Direction,4> dirs)143 bool isClockwise( std::array<Direction, 4> dirs )
144 {
145 const std::array<Direction, 4> cwdirs = { Direction::Up, Direction::Right, Direction::Down, Direction::Left };
146 return matchesOrientation( dirs, cwdirs );
147 }
148
149 /**
150 * Checks whether the 4 directions in dirs make up a counter-clockwise
151 * rectangle.
152 */
isCounterClockwise(std::array<Direction,4> dirs)153 bool isCounterClockwise( std::array<Direction, 4> dirs )
154 {
155 const std::array<Direction, 4> ccwdirs = { Direction::Right, Direction::Up, Direction::Left, Direction::Down };
156 return matchesOrientation( dirs, ccwdirs );
157 }
158
159
isAxisParallelRectangle(double maximumDeviation,bool simpleRectanglesOnly) const160 bool QgsInternalGeometryEngine::isAxisParallelRectangle( double maximumDeviation, bool simpleRectanglesOnly ) const
161 {
162 if ( QgsWkbTypes::flatType( mGeometry->wkbType() ) != QgsWkbTypes::Polygon )
163 return false;
164
165 const QgsPolygon *polygon = qgsgeometry_cast< const QgsPolygon * >( mGeometry );
166 if ( !polygon->exteriorRing() || polygon->numInteriorRings() > 0 )
167 return false;
168
169 const int vertexCount = polygon->exteriorRing()->numPoints();
170 if ( vertexCount < 4 )
171 return false;
172 else if ( simpleRectanglesOnly && ( vertexCount != 5 || !polygon->exteriorRing()->isClosed() ) )
173 return false;
174
175 bool found4Dirs;
176 std::array<Direction, 4> dirs;
177 std::tie( found4Dirs, dirs ) = getEdgeDirections( polygon, std::tan( maximumDeviation * M_PI / 180 ) );
178
179 return found4Dirs && ( isCounterClockwise( dirs ) || isClockwise( dirs ) );
180 }
181
182 /***************************************************************************
183 * This class is considered CRITICAL and any change MUST be accompanied with
184 * full unit tests.
185 * See details in QEP #17
186 ****************************************************************************/
187
extrude(double x,double y) const188 QgsGeometry QgsInternalGeometryEngine::extrude( double x, double y ) const
189 {
190 mLastError.clear();
191 QVector<QgsLineString *> linesToProcess;
192
193 const QgsMultiCurve *multiCurve = qgsgeometry_cast< const QgsMultiCurve * >( mGeometry );
194 if ( multiCurve )
195 {
196 linesToProcess.reserve( multiCurve->partCount() );
197 for ( int i = 0; i < multiCurve->partCount(); ++i )
198 {
199 linesToProcess << static_cast<QgsLineString *>( multiCurve->geometryN( i )->clone() );
200 }
201 }
202
203 const QgsCurve *curve = qgsgeometry_cast< const QgsCurve * >( mGeometry );
204 if ( curve )
205 {
206 linesToProcess << static_cast<QgsLineString *>( curve->segmentize() );
207 }
208
209 std::unique_ptr<QgsMultiPolygon> multipolygon( linesToProcess.size() > 1 ? new QgsMultiPolygon() : nullptr );
210 QgsPolygon *polygon = nullptr;
211
212 if ( !linesToProcess.empty() )
213 {
214 std::unique_ptr< QgsLineString > secondline;
215 for ( QgsLineString *line : std::as_const( linesToProcess ) )
216 {
217 QTransform transform = QTransform::fromTranslate( x, y );
218
219 secondline.reset( line->reversed() );
220 secondline->transform( transform );
221
222 line->append( secondline.get() );
223 line->addVertex( line->pointN( 0 ) );
224
225 polygon = new QgsPolygon();
226 polygon->setExteriorRing( line );
227
228 if ( multipolygon )
229 multipolygon->addGeometry( polygon );
230 }
231
232 if ( multipolygon )
233 return QgsGeometry( multipolygon.release() );
234 else
235 return QgsGeometry( polygon );
236 }
237
238 return QgsGeometry();
239 }
240
241
242
243 // polylabel implementation
244 // ported from the original JavaScript implementation developed by Vladimir Agafonkin
245 // originally licensed under the ISC License
246
247 /// @cond PRIVATE
248 class Cell
249 {
250 public:
Cell(double x,double y,double h,const QgsPolygon * polygon)251 Cell( double x, double y, double h, const QgsPolygon *polygon )
252 : x( x )
253 , y( y )
254 , h( h )
255 , d( polygon->pointDistanceToBoundary( x, y ) )
256 , max( d + h * M_SQRT2 )
257 {}
258
259 //! Cell center x
260 double x;
261 //! Cell center y
262 double y;
263 //! Half the cell size
264 double h;
265 //! Distance from cell center to polygon
266 double d;
267 //! Maximum distance to polygon within a cell
268 double max;
269 };
270
271 struct GreaterThanByMax
272 {
operator ()GreaterThanByMax273 bool operator()( const Cell *lhs, const Cell *rhs )
274 {
275 return rhs->max > lhs->max;
276 }
277 };
278
getCentroidCell(const QgsPolygon * polygon)279 Cell *getCentroidCell( const QgsPolygon *polygon )
280 {
281 double area = 0;
282 double x = 0;
283 double y = 0;
284
285 const QgsLineString *exterior = static_cast< const QgsLineString *>( polygon->exteriorRing() );
286 int len = exterior->numPoints() - 1; //assume closed
287 for ( int i = 0, j = len - 1; i < len; j = i++ )
288 {
289 double aX = exterior->xAt( i );
290 double aY = exterior->yAt( i );
291 double bX = exterior->xAt( j );
292 double bY = exterior->yAt( j );
293 double f = aX * bY - bX * aY;
294 x += ( aX + bX ) * f;
295 y += ( aY + bY ) * f;
296 area += f * 3;
297 }
298 if ( area == 0 )
299 return new Cell( exterior->xAt( 0 ), exterior->yAt( 0 ), 0, polygon );
300 else
301 return new Cell( x / area, y / area, 0.0, polygon );
302 }
303
surfacePoleOfInaccessibility(const QgsSurface * surface,double precision,double & distanceFromBoundary)304 QgsPoint surfacePoleOfInaccessibility( const QgsSurface *surface, double precision, double &distanceFromBoundary )
305 {
306 std::unique_ptr< QgsPolygon > segmentizedPoly;
307 const QgsPolygon *polygon = qgsgeometry_cast< const QgsPolygon * >( surface );
308 if ( !polygon )
309 {
310 segmentizedPoly.reset( static_cast< QgsPolygon *>( surface->segmentize() ) );
311 polygon = segmentizedPoly.get();
312 }
313
314 // start with the bounding box
315 QgsRectangle bounds = polygon->boundingBox();
316
317 // initial parameters
318 double cellSize = std::min( bounds.width(), bounds.height() );
319
320 if ( qgsDoubleNear( cellSize, 0.0 ) )
321 return QgsPoint( bounds.xMinimum(), bounds.yMinimum() );
322
323 double h = cellSize / 2.0;
324 std::priority_queue< Cell *, std::vector<Cell *>, GreaterThanByMax > cellQueue;
325
326 // cover polygon with initial cells
327 for ( double x = bounds.xMinimum(); x < bounds.xMaximum(); x += cellSize )
328 {
329 for ( double y = bounds.yMinimum(); y < bounds.yMaximum(); y += cellSize )
330 {
331 cellQueue.push( new Cell( x + h, y + h, h, polygon ) );
332 }
333 }
334
335 // take centroid as the first best guess
336 std::unique_ptr< Cell > bestCell( getCentroidCell( polygon ) );
337
338 // special case for rectangular polygons
339 std::unique_ptr< Cell > bboxCell( new Cell( bounds.xMinimum() + bounds.width() / 2.0,
340 bounds.yMinimum() + bounds.height() / 2.0,
341 0, polygon ) );
342 if ( bboxCell->d > bestCell->d )
343 {
344 bestCell = std::move( bboxCell );
345 }
346
347 while ( !cellQueue.empty() )
348 {
349 // pick the most promising cell from the queue
350 std::unique_ptr< Cell > cell( cellQueue.top() );
351 cellQueue.pop();
352 Cell *currentCell = cell.get();
353
354 // update the best cell if we found a better one
355 if ( currentCell->d > bestCell->d )
356 {
357 bestCell = std::move( cell );
358 }
359
360 // do not drill down further if there's no chance of a better solution
361 if ( currentCell->max - bestCell->d <= precision )
362 continue;
363
364 // split the cell into four cells
365 h = currentCell->h / 2.0;
366 cellQueue.push( new Cell( currentCell->x - h, currentCell->y - h, h, polygon ) );
367 cellQueue.push( new Cell( currentCell->x + h, currentCell->y - h, h, polygon ) );
368 cellQueue.push( new Cell( currentCell->x - h, currentCell->y + h, h, polygon ) );
369 cellQueue.push( new Cell( currentCell->x + h, currentCell->y + h, h, polygon ) );
370 }
371
372 distanceFromBoundary = bestCell->d;
373
374 return QgsPoint( bestCell->x, bestCell->y );
375 }
376
377 ///@endcond
378
poleOfInaccessibility(double precision,double * distanceFromBoundary) const379 QgsGeometry QgsInternalGeometryEngine::poleOfInaccessibility( double precision, double *distanceFromBoundary ) const
380 {
381 mLastError.clear();
382 if ( distanceFromBoundary )
383 *distanceFromBoundary = std::numeric_limits<double>::max();
384
385 if ( !mGeometry || mGeometry->isEmpty() )
386 return QgsGeometry();
387
388 if ( precision <= 0 )
389 return QgsGeometry();
390
391 if ( const QgsGeometryCollection *gc = qgsgeometry_cast< const QgsGeometryCollection *>( mGeometry ) )
392 {
393 int numGeom = gc->numGeometries();
394 double maxDist = 0;
395 QgsPoint bestPoint;
396 for ( int i = 0; i < numGeom; ++i )
397 {
398 const QgsSurface *surface = qgsgeometry_cast< const QgsSurface * >( gc->geometryN( i ) );
399 if ( !surface )
400 continue;
401
402 double dist = std::numeric_limits<double>::max();
403 QgsPoint p = surfacePoleOfInaccessibility( surface, precision, dist );
404 if ( dist > maxDist )
405 {
406 maxDist = dist;
407 bestPoint = p;
408 }
409 }
410
411 if ( bestPoint.isEmpty() )
412 return QgsGeometry();
413
414 if ( distanceFromBoundary )
415 *distanceFromBoundary = maxDist;
416 return QgsGeometry( new QgsPoint( bestPoint ) );
417 }
418 else
419 {
420 const QgsSurface *surface = qgsgeometry_cast< const QgsSurface * >( mGeometry );
421 if ( !surface )
422 return QgsGeometry();
423
424 double dist = std::numeric_limits<double>::max();
425 QgsPoint p = surfacePoleOfInaccessibility( surface, precision, dist );
426 if ( p.isEmpty() )
427 return QgsGeometry();
428
429 if ( distanceFromBoundary )
430 *distanceFromBoundary = dist;
431 return QgsGeometry( new QgsPoint( p ) );
432 }
433 }
434
435
436 // helpers for orthogonalize
437 // adapted from original code in potlatch/id osm editor
438
dotProductWithinAngleTolerance(double dotProduct,double lowerThreshold,double upperThreshold)439 bool dotProductWithinAngleTolerance( double dotProduct, double lowerThreshold, double upperThreshold )
440 {
441 return lowerThreshold > std::fabs( dotProduct ) || std::fabs( dotProduct ) > upperThreshold;
442 }
443
normalizedDotProduct(const QgsPoint & a,const QgsPoint & b,const QgsPoint & c)444 double normalizedDotProduct( const QgsPoint &a, const QgsPoint &b, const QgsPoint &c )
445 {
446 QgsVector p = a - b;
447 QgsVector q = c - b;
448
449 if ( p.length() > 0 )
450 p = p.normalized();
451 if ( q.length() > 0 )
452 q = q.normalized();
453
454 return p * q;
455 }
456
squareness(QgsLineString * ring,double lowerThreshold,double upperThreshold)457 double squareness( QgsLineString *ring, double lowerThreshold, double upperThreshold )
458 {
459 double sum = 0.0;
460
461 bool isClosed = ring->isClosed();
462 int numPoints = ring->numPoints();
463 QgsPoint a;
464 QgsPoint b;
465 QgsPoint c;
466
467 for ( int i = 0; i < numPoints; ++i )
468 {
469 if ( !isClosed && i == numPoints - 1 )
470 break;
471 else if ( !isClosed && i == 0 )
472 {
473 b = ring->pointN( 0 );
474 c = ring->pointN( 1 );
475 }
476 else
477 {
478 if ( i == 0 )
479 {
480 a = ring->pointN( numPoints - 1 );
481 b = ring->pointN( 0 );
482 }
483 if ( i == numPoints - 1 )
484 c = ring->pointN( 0 );
485 else
486 c = ring->pointN( i + 1 );
487
488 double dotProduct = normalizedDotProduct( a, b, c );
489 if ( !dotProductWithinAngleTolerance( dotProduct, lowerThreshold, upperThreshold ) )
490 continue;
491
492 sum += 2.0 * std::min( std::fabs( dotProduct - 1.0 ), std::min( std::fabs( dotProduct ), std::fabs( dotProduct + 1 ) ) );
493 }
494 a = b;
495 b = c;
496 }
497
498 return sum;
499 }
500
calcMotion(const QgsPoint & a,const QgsPoint & b,const QgsPoint & c,double lowerThreshold,double upperThreshold)501 QgsVector calcMotion( const QgsPoint &a, const QgsPoint &b, const QgsPoint &c,
502 double lowerThreshold, double upperThreshold )
503 {
504 QgsVector p = a - b;
505 QgsVector q = c - b;
506
507 if ( qgsDoubleNear( p.length(), 0.0 ) || qgsDoubleNear( q.length(), 0.0 ) )
508 return QgsVector( 0, 0 );
509
510 // 2.0 is a magic number from the original JOSM source code
511 double scale = 2.0 * std::min( p.length(), q.length() );
512
513 p = p.normalized();
514 q = q.normalized();
515 double dotProduct = p * q;
516
517 if ( !dotProductWithinAngleTolerance( dotProduct, lowerThreshold, upperThreshold ) )
518 {
519 return QgsVector( 0, 0 );
520 }
521
522 // wonderful nasty hack which has survived through JOSM -> id -> QGIS
523 // to deal with almost-straight segments (angle is closer to 180 than to 90/270).
524 if ( dotProduct < -M_SQRT1_2 )
525 dotProduct += 1.0;
526
527 QgsVector new_v = p + q;
528 // 0.1 magic number from JOSM implementation - think this is to limit each iterative step
529 return new_v.normalized() * ( 0.1 * dotProduct * scale );
530 }
531
doOrthogonalize(QgsLineString * ring,int iterations,double tolerance,double lowerThreshold,double upperThreshold)532 QgsLineString *doOrthogonalize( QgsLineString *ring, int iterations, double tolerance, double lowerThreshold, double upperThreshold )
533 {
534 double minScore = std::numeric_limits<double>::max();
535
536 bool isClosed = ring->isClosed();
537 int numPoints = ring->numPoints();
538
539 std::unique_ptr< QgsLineString > best( ring->clone() );
540
541 QVector< QgsVector > /* yep */ motions;
542 motions.reserve( numPoints );
543
544 for ( int it = 0; it < iterations; ++it )
545 {
546 motions.resize( 0 ); // avoid re-allocations
547
548 // first loop through an calculate all motions
549 QgsPoint a;
550 QgsPoint b;
551 QgsPoint c;
552 for ( int i = 0; i < numPoints; ++i )
553 {
554 if ( isClosed && i == numPoints - 1 )
555 motions << motions.at( 0 ); //closed ring, so last point follows first point motion
556 else if ( !isClosed && ( i == 0 || i == numPoints - 1 ) )
557 {
558 b = ring->pointN( 0 );
559 c = ring->pointN( 1 );
560 motions << QgsVector( 0, 0 ); //non closed line, leave first/last vertex untouched
561 }
562 else
563 {
564 if ( i == 0 )
565 {
566 a = ring->pointN( numPoints - 1 );
567 b = ring->pointN( 0 );
568 }
569 if ( i == numPoints - 1 )
570 c = ring->pointN( 0 );
571 else
572 c = ring->pointN( i + 1 );
573
574 motions << calcMotion( a, b, c, lowerThreshold, upperThreshold );
575 }
576 a = b;
577 b = c;
578 }
579
580 // then apply them
581 for ( int i = 0; i < numPoints; ++i )
582 {
583 ring->setXAt( i, ring->xAt( i ) + motions.at( i ).x() );
584 ring->setYAt( i, ring->yAt( i ) + motions.at( i ).y() );
585 }
586
587 double newScore = squareness( ring, lowerThreshold, upperThreshold );
588 if ( newScore < minScore )
589 {
590 best.reset( ring->clone() );
591 minScore = newScore;
592 }
593
594 if ( minScore < tolerance )
595 break;
596 }
597
598 delete ring;
599
600 return best.release();
601 }
602
603
orthogonalizeGeom(const QgsAbstractGeometry * geom,int maxIterations,double tolerance,double lowerThreshold,double upperThreshold)604 QgsAbstractGeometry *orthogonalizeGeom( const QgsAbstractGeometry *geom, int maxIterations, double tolerance, double lowerThreshold, double upperThreshold )
605 {
606 std::unique_ptr< QgsAbstractGeometry > segmentizedCopy;
607 if ( QgsWkbTypes::isCurvedType( geom->wkbType() ) )
608 {
609 segmentizedCopy.reset( geom->segmentize() );
610 geom = segmentizedCopy.get();
611 }
612
613 if ( QgsWkbTypes::geometryType( geom->wkbType() ) == QgsWkbTypes::LineGeometry )
614 {
615 return doOrthogonalize( static_cast< QgsLineString * >( geom->clone() ),
616 maxIterations, tolerance, lowerThreshold, upperThreshold );
617 }
618 else
619 {
620 // polygon
621 const QgsPolygon *polygon = static_cast< const QgsPolygon * >( geom );
622 QgsPolygon *result = new QgsPolygon();
623
624 result->setExteriorRing( doOrthogonalize( static_cast< QgsLineString * >( polygon->exteriorRing()->clone() ),
625 maxIterations, tolerance, lowerThreshold, upperThreshold ) );
626 for ( int i = 0; i < polygon->numInteriorRings(); ++i )
627 {
628 result->addInteriorRing( doOrthogonalize( static_cast< QgsLineString * >( polygon->interiorRing( i )->clone() ),
629 maxIterations, tolerance, lowerThreshold, upperThreshold ) );
630 }
631
632 return result;
633 }
634 }
635
orthogonalize(double tolerance,int maxIterations,double angleThreshold) const636 QgsGeometry QgsInternalGeometryEngine::orthogonalize( double tolerance, int maxIterations, double angleThreshold ) const
637 {
638 mLastError.clear();
639 if ( !mGeometry || ( QgsWkbTypes::geometryType( mGeometry->wkbType() ) != QgsWkbTypes::LineGeometry
640 && QgsWkbTypes::geometryType( mGeometry->wkbType() ) != QgsWkbTypes::PolygonGeometry ) )
641 {
642 return QgsGeometry();
643 }
644
645 double lowerThreshold = std::cos( ( 90 - angleThreshold ) * M_PI / 180.00 );
646 double upperThreshold = std::cos( angleThreshold * M_PI / 180.0 );
647
648 if ( const QgsGeometryCollection *gc = qgsgeometry_cast< const QgsGeometryCollection *>( mGeometry ) )
649 {
650 int numGeom = gc->numGeometries();
651 QVector< QgsAbstractGeometry * > geometryList;
652 geometryList.reserve( numGeom );
653 for ( int i = 0; i < numGeom; ++i )
654 {
655 geometryList << orthogonalizeGeom( gc->geometryN( i ), maxIterations, tolerance, lowerThreshold, upperThreshold );
656 }
657
658 QgsGeometry first = QgsGeometry( geometryList.takeAt( 0 ) );
659 for ( QgsAbstractGeometry *g : std::as_const( geometryList ) )
660 {
661 first.addPart( g );
662 }
663 return first;
664 }
665 else
666 {
667 return QgsGeometry( orthogonalizeGeom( mGeometry, maxIterations, tolerance, lowerThreshold, upperThreshold ) );
668 }
669 }
670
671 // if extraNodesPerSegment < 0, then use distance based mode
doDensify(const QgsLineString * ring,int extraNodesPerSegment=-1,double distance=1)672 QgsLineString *doDensify( const QgsLineString *ring, int extraNodesPerSegment = -1, double distance = 1 )
673 {
674 QVector< double > outX;
675 QVector< double > outY;
676 QVector< double > outZ;
677 QVector< double > outM;
678 double multiplier = 1.0 / double( extraNodesPerSegment + 1 );
679
680 int nPoints = ring->numPoints();
681 outX.reserve( ( extraNodesPerSegment + 1 ) * nPoints );
682 outY.reserve( ( extraNodesPerSegment + 1 ) * nPoints );
683 bool withZ = ring->is3D();
684 if ( withZ )
685 outZ.reserve( ( extraNodesPerSegment + 1 ) * nPoints );
686 bool withM = ring->isMeasure();
687 if ( withM )
688 outM.reserve( ( extraNodesPerSegment + 1 ) * nPoints );
689 double x1 = 0;
690 double x2 = 0;
691 double y1 = 0;
692 double y2 = 0;
693 double z1 = 0;
694 double z2 = 0;
695 double m1 = 0;
696 double m2 = 0;
697 double xOut = 0;
698 double yOut = 0;
699 double zOut = 0;
700 double mOut = 0;
701 int extraNodesThisSegment = extraNodesPerSegment;
702 for ( int i = 0; i < nPoints - 1; ++i )
703 {
704 x1 = ring->xAt( i );
705 x2 = ring->xAt( i + 1 );
706 y1 = ring->yAt( i );
707 y2 = ring->yAt( i + 1 );
708 if ( withZ )
709 {
710 z1 = ring->zAt( i );
711 z2 = ring->zAt( i + 1 );
712 }
713 if ( withM )
714 {
715 m1 = ring->mAt( i );
716 m2 = ring->mAt( i + 1 );
717 }
718
719 outX << x1;
720 outY << y1;
721 if ( withZ )
722 outZ << z1;
723 if ( withM )
724 outM << m1;
725
726 if ( extraNodesPerSegment < 0 )
727 {
728 // distance mode
729 extraNodesThisSegment = std::floor( std::sqrt( ( x2 - x1 ) * ( x2 - x1 ) + ( y2 - y1 ) * ( y2 - y1 ) ) / distance );
730 if ( extraNodesThisSegment >= 1 )
731 multiplier = 1.0 / ( extraNodesThisSegment + 1 );
732 }
733
734 for ( int j = 0; j < extraNodesThisSegment; ++j )
735 {
736 double delta = multiplier * ( j + 1 );
737 xOut = x1 + delta * ( x2 - x1 );
738 yOut = y1 + delta * ( y2 - y1 );
739 if ( withZ )
740 zOut = z1 + delta * ( z2 - z1 );
741 if ( withM )
742 mOut = m1 + delta * ( m2 - m1 );
743
744 outX << xOut;
745 outY << yOut;
746 if ( withZ )
747 outZ << zOut;
748 if ( withM )
749 outM << mOut;
750 }
751 }
752 outX << ring->xAt( nPoints - 1 );
753 outY << ring->yAt( nPoints - 1 );
754 if ( withZ )
755 outZ << ring->zAt( nPoints - 1 );
756 if ( withM )
757 outM << ring->mAt( nPoints - 1 );
758
759 QgsLineString *result = new QgsLineString( outX, outY, outZ, outM );
760 return result;
761 }
762
densifyGeometry(const QgsAbstractGeometry * geom,int extraNodesPerSegment=1,double distance=1)763 QgsAbstractGeometry *densifyGeometry( const QgsAbstractGeometry *geom, int extraNodesPerSegment = 1, double distance = 1 )
764 {
765 std::unique_ptr< QgsAbstractGeometry > segmentizedCopy;
766 if ( QgsWkbTypes::isCurvedType( geom->wkbType() ) )
767 {
768 segmentizedCopy.reset( geom->segmentize() );
769 geom = segmentizedCopy.get();
770 }
771
772 if ( QgsWkbTypes::geometryType( geom->wkbType() ) == QgsWkbTypes::LineGeometry )
773 {
774 return doDensify( static_cast< const QgsLineString * >( geom ), extraNodesPerSegment, distance );
775 }
776 else
777 {
778 // polygon
779 const QgsPolygon *polygon = static_cast< const QgsPolygon * >( geom );
780 QgsPolygon *result = new QgsPolygon();
781
782 result->setExteriorRing( doDensify( static_cast< const QgsLineString * >( polygon->exteriorRing() ),
783 extraNodesPerSegment, distance ) );
784 for ( int i = 0; i < polygon->numInteriorRings(); ++i )
785 {
786 result->addInteriorRing( doDensify( static_cast< const QgsLineString * >( polygon->interiorRing( i ) ),
787 extraNodesPerSegment, distance ) );
788 }
789
790 return result;
791 }
792 }
793
densifyByCount(int extraNodesPerSegment) const794 QgsGeometry QgsInternalGeometryEngine::densifyByCount( int extraNodesPerSegment ) const
795 {
796 mLastError.clear();
797 if ( !mGeometry )
798 {
799 return QgsGeometry();
800 }
801
802 if ( QgsWkbTypes::geometryType( mGeometry->wkbType() ) == QgsWkbTypes::PointGeometry )
803 {
804 return QgsGeometry( mGeometry->clone() ); // point geometry, nothing to do
805 }
806
807 if ( const QgsGeometryCollection *gc = qgsgeometry_cast< const QgsGeometryCollection *>( mGeometry ) )
808 {
809 int numGeom = gc->numGeometries();
810 QVector< QgsAbstractGeometry * > geometryList;
811 geometryList.reserve( numGeom );
812 for ( int i = 0; i < numGeom; ++i )
813 {
814 geometryList << densifyGeometry( gc->geometryN( i ), extraNodesPerSegment );
815 }
816
817 QgsGeometry first = QgsGeometry( geometryList.takeAt( 0 ) );
818 for ( QgsAbstractGeometry *g : std::as_const( geometryList ) )
819 {
820 first.addPart( g );
821 }
822 return first;
823 }
824 else
825 {
826 return QgsGeometry( densifyGeometry( mGeometry, extraNodesPerSegment ) );
827 }
828 }
829
densifyByDistance(double distance) const830 QgsGeometry QgsInternalGeometryEngine::densifyByDistance( double distance ) const
831 {
832 mLastError.clear();
833 if ( !mGeometry )
834 {
835 return QgsGeometry();
836 }
837
838 if ( QgsWkbTypes::geometryType( mGeometry->wkbType() ) == QgsWkbTypes::PointGeometry )
839 {
840 return QgsGeometry( mGeometry->clone() ); // point geometry, nothing to do
841 }
842
843 if ( const QgsGeometryCollection *gc = qgsgeometry_cast< const QgsGeometryCollection *>( mGeometry ) )
844 {
845 int numGeom = gc->numGeometries();
846 QVector< QgsAbstractGeometry * > geometryList;
847 geometryList.reserve( numGeom );
848 for ( int i = 0; i < numGeom; ++i )
849 {
850 geometryList << densifyGeometry( gc->geometryN( i ), -1, distance );
851 }
852
853 QgsGeometry first = QgsGeometry( geometryList.takeAt( 0 ) );
854 for ( QgsAbstractGeometry *g : std::as_const( geometryList ) )
855 {
856 first.addPart( g );
857 }
858 return first;
859 }
860 else
861 {
862 return QgsGeometry( densifyGeometry( mGeometry, -1, distance ) );
863 }
864 }
865
866 ///@cond PRIVATE
867 //
868 // QgsLineSegmentDistanceComparer
869 //
870
871 // adapted for QGIS geometry classes from original work at https://github.com/trylock/visibility by trylock
operator ()(QgsLineSegment2D ab,QgsLineSegment2D cd) const872 bool QgsLineSegmentDistanceComparer::operator()( QgsLineSegment2D ab, QgsLineSegment2D cd ) const
873 {
874 Q_ASSERT_X( ab.pointLeftOfLine( mOrigin ) != 0,
875 "line_segment_dist_comparer",
876 "AB must not be collinear with the origin." );
877 Q_ASSERT_X( cd.pointLeftOfLine( mOrigin ) != 0,
878 "line_segment_dist_comparer",
879 "CD must not be collinear with the origin." );
880
881 // flip the segments so that if there are common endpoints,
882 // they will be the segment's start points
883 if ( ab.end() == cd.start() || ab.end() == cd.end() )
884 ab.reverse();
885 if ( ab.start() == cd.end() )
886 cd.reverse();
887
888 // cases with common endpoints
889 if ( ab.start() == cd.start() )
890 {
891 const int oad = QgsGeometryUtils::leftOfLine( cd.endX(), cd.endY(), mOrigin.x(), mOrigin.y(), ab.startX(), ab.startY() );
892 const int oab = ab.pointLeftOfLine( mOrigin );
893 if ( ab.end() == cd.end() || oad != oab )
894 return false;
895 else
896 return ab.pointLeftOfLine( cd.end() ) != oab;
897 }
898 else
899 {
900 // cases without common endpoints
901 const int cda = cd.pointLeftOfLine( ab.start() );
902 const int cdb = cd.pointLeftOfLine( ab.end() );
903 if ( cdb == 0 && cda == 0 )
904 {
905 return mOrigin.sqrDist( ab.start() ) < mOrigin.sqrDist( cd.start() );
906 }
907 else if ( cda == cdb || cda == 0 || cdb == 0 )
908 {
909 const int cdo = cd.pointLeftOfLine( mOrigin );
910 return cdo == cda || cdo == cdb;
911 }
912 else
913 {
914 const int abo = ab.pointLeftOfLine( mOrigin );
915 return abo != ab.pointLeftOfLine( cd.start() );
916 }
917 }
918 }
919
920 //
921 // QgsClockwiseAngleComparer
922 //
923
operator ()(const QgsPointXY & a,const QgsPointXY & b) const924 bool QgsClockwiseAngleComparer::operator()( const QgsPointXY &a, const QgsPointXY &b ) const
925 {
926 const bool aIsLeft = a.x() < mVertex.x();
927 const bool bIsLeft = b.x() < mVertex.x();
928 if ( aIsLeft != bIsLeft )
929 return bIsLeft;
930
931 if ( qgsDoubleNear( a.x(), mVertex.x() ) && qgsDoubleNear( b.x(), mVertex.x() ) )
932 {
933 if ( a.y() >= mVertex.y() || b.y() >= mVertex.y() )
934 {
935 return b.y() < a.y();
936 }
937 else
938 {
939 return a.y() < b.y();
940 }
941 }
942 else
943 {
944 const QgsVector oa = a - mVertex;
945 const QgsVector ob = b - mVertex;
946 const double det = oa.crossProduct( ob );
947 if ( qgsDoubleNear( det, 0.0 ) )
948 {
949 return oa.lengthSquared() < ob.lengthSquared();
950 }
951 else
952 {
953 return det < 0;
954 }
955 }
956 }
957
958 ///@endcond PRIVATE
959
960 //
961 // QgsRay2D
962 //
963
intersects(const QgsLineSegment2D & segment,QgsPointXY & intersectPoint) const964 bool QgsRay2D::intersects( const QgsLineSegment2D &segment, QgsPointXY &intersectPoint ) const
965 {
966 const QgsVector ao = origin - segment.start();
967 const QgsVector ab = segment.end() - segment.start();
968 const double det = ab.crossProduct( direction );
969 if ( qgsDoubleNear( det, 0.0 ) )
970 {
971 const int abo = segment.pointLeftOfLine( origin );
972 if ( abo != 0 )
973 {
974 return false;
975 }
976 else
977 {
978 const double distA = ao * direction;
979 const double distB = ( origin - segment.end() ) * direction;
980
981 if ( distA > 0 && distB > 0 )
982 {
983 return false;
984 }
985 else
986 {
987 if ( ( distA > 0 ) != ( distB > 0 ) )
988 intersectPoint = origin;
989 else if ( distA > distB ) // at this point, both distances are negative
990 intersectPoint = segment.start(); // hence the nearest point is A
991 else
992 intersectPoint = segment.end();
993 return true;
994 }
995 }
996 }
997 else
998 {
999 const double u = ao.crossProduct( direction ) / det;
1000 if ( u < 0.0 || 1.0 < u )
1001 {
1002 return false;
1003 }
1004 else
1005 {
1006 const double t = -ab.crossProduct( ao ) / det;
1007 intersectPoint = origin + direction * t;
1008 return qgsDoubleNear( t, 0.0 ) || t > 0;
1009 }
1010 }
1011 }
1012
generateSegmentCurve(const QgsPoint & center1,const double radius1,const QgsPoint & center2,const double radius2)1013 QVector<QgsPointXY> generateSegmentCurve( const QgsPoint ¢er1, const double radius1, const QgsPoint ¢er2, const double radius2 )
1014 {
1015 // ensure that first circle is smaller than second
1016 if ( radius1 > radius2 )
1017 return generateSegmentCurve( center2, radius2, center1, radius1 );
1018
1019 QgsPointXY t1;
1020 QgsPointXY t2;
1021 QgsPointXY t3;
1022 QgsPointXY t4;
1023 QVector<QgsPointXY> points;
1024 if ( QgsGeometryUtils::circleCircleOuterTangents( center1, radius1, center2, radius2, t1, t2, t3, t4 ) )
1025 {
1026 points << t1
1027 << t2
1028 << t4
1029 << t3;
1030 }
1031 return points;
1032 }
1033
1034 // partially ported from JTS VariableWidthBuffer,
1035 // https://github.com/topobyte/jts/blob/master/jts-lab/src/main/java/com/vividsolutions/jts/operation/buffer/VariableWidthBuffer.java
1036
variableWidthBuffer(int segments,const std::function<std::unique_ptr<double[]> (const QgsLineString * line)> & widthFunction) const1037 QgsGeometry QgsInternalGeometryEngine::variableWidthBuffer( int segments, const std::function< std::unique_ptr< double[] >( const QgsLineString *line ) > &widthFunction ) const
1038 {
1039 mLastError.clear();
1040 if ( !mGeometry )
1041 {
1042 return QgsGeometry();
1043 }
1044
1045 std::vector< std::unique_ptr<QgsLineString > > linesToProcess;
1046
1047 const QgsMultiCurve *multiCurve = qgsgeometry_cast< const QgsMultiCurve * >( mGeometry );
1048 if ( multiCurve )
1049 {
1050 for ( int i = 0; i < multiCurve->partCount(); ++i )
1051 {
1052 if ( static_cast< const QgsCurve * >( multiCurve->geometryN( i ) )->nCoordinates() == 0 )
1053 continue; // skip 0 length lines
1054
1055 linesToProcess.emplace_back( static_cast<QgsLineString *>( multiCurve->geometryN( i )->clone() ) );
1056 }
1057 }
1058
1059 const QgsCurve *curve = qgsgeometry_cast< const QgsCurve * >( mGeometry );
1060 if ( curve )
1061 {
1062 if ( curve->nCoordinates() > 0 )
1063 linesToProcess.emplace_back( static_cast<QgsLineString *>( curve->segmentize() ) );
1064 }
1065
1066 if ( linesToProcess.empty() )
1067 {
1068 QgsGeometry g;
1069 g.mLastError = QStringLiteral( "Input geometry was not a curve type geometry" );
1070 return g;
1071 }
1072
1073 QVector<QgsGeometry> bufferedLines;
1074 bufferedLines.reserve( linesToProcess.size() );
1075
1076 for ( std::unique_ptr< QgsLineString > &line : linesToProcess )
1077 {
1078 QVector<QgsGeometry> parts;
1079 QgsPoint prevPoint;
1080 double prevRadius = 0;
1081 QgsGeometry prevCircle;
1082
1083 std::unique_ptr< double[] > widths = widthFunction( line.get() ) ;
1084 for ( int i = 0; i < line->nCoordinates(); ++i )
1085 {
1086 QgsPoint thisPoint = line->pointN( i );
1087 QgsGeometry thisCircle;
1088 double thisRadius = widths[ i ] / 2.0;
1089 if ( thisRadius > 0 )
1090 {
1091 QgsGeometry p = QgsGeometry( thisPoint.clone() );
1092
1093 QgsCircle circ( thisPoint, thisRadius );
1094 thisCircle = QgsGeometry( circ.toPolygon( segments * 4 ) );
1095 parts << thisCircle;
1096 }
1097 else
1098 {
1099 thisCircle = QgsGeometry( thisPoint.clone() );
1100 }
1101
1102 if ( i > 0 )
1103 {
1104 if ( prevRadius > 0 || thisRadius > 0 )
1105 {
1106 QVector< QgsPointXY > points = generateSegmentCurve( prevPoint, prevRadius, thisPoint, thisRadius );
1107 if ( !points.empty() )
1108 {
1109 // snap points to circle vertices
1110
1111 int atVertex = 0;
1112 int beforeVertex = 0;
1113 int afterVertex = 0;
1114 double sqrDist = 0;
1115 double sqrDistPrev = 0;
1116 for ( int j = 0; j < points.count(); ++j )
1117 {
1118 QgsPointXY pA = prevCircle.closestVertex( points.at( j ), atVertex, beforeVertex, afterVertex, sqrDistPrev );
1119 QgsPointXY pB = thisCircle.closestVertex( points.at( j ), atVertex, beforeVertex, afterVertex, sqrDist );
1120 points[j] = sqrDistPrev < sqrDist ? pA : pB;
1121 }
1122 // close ring
1123 points.append( points.at( 0 ) );
1124
1125 std::unique_ptr< QgsPolygon > poly = std::make_unique< QgsPolygon >();
1126 poly->setExteriorRing( new QgsLineString( points ) );
1127 if ( poly->area() > 0 )
1128 parts << QgsGeometry( std::move( poly ) );
1129 }
1130 }
1131 }
1132 prevPoint = thisPoint;
1133 prevRadius = thisRadius;
1134 prevCircle = thisCircle;
1135 }
1136
1137 bufferedLines << QgsGeometry::unaryUnion( parts );
1138 }
1139
1140 return QgsGeometry::collectGeometry( bufferedLines );
1141 }
1142
taperedBuffer(double start,double end,int segments) const1143 QgsGeometry QgsInternalGeometryEngine::taperedBuffer( double start, double end, int segments ) const
1144 {
1145 mLastError.clear();
1146 start = std::fabs( start );
1147 end = std::fabs( end );
1148
1149 auto interpolateWidths = [ start, end ]( const QgsLineString * line )->std::unique_ptr< double [] >
1150 {
1151 // ported from JTS VariableWidthBuffer,
1152 // https://github.com/topobyte/jts/blob/master/jts-lab/src/main/java/com/vividsolutions/jts/operation/buffer/VariableWidthBuffer.java
1153 std::unique_ptr< double [] > widths( new double[ line->nCoordinates() ] );
1154 widths[0] = start;
1155 widths[line->nCoordinates() - 1] = end;
1156
1157 double lineLength = line->length();
1158 double currentLength = 0;
1159 QgsPoint prevPoint = line->pointN( 0 );
1160 for ( int i = 1; i < line->nCoordinates() - 1; ++i )
1161 {
1162 QgsPoint point = line->pointN( i );
1163 double segmentLength = point.distance( prevPoint );
1164 currentLength += segmentLength;
1165 double lengthFraction = lineLength > 0 ? currentLength / lineLength : 1;
1166 double delta = lengthFraction * ( end - start );
1167 widths[i] = start + delta;
1168 prevPoint = point;
1169 }
1170 return widths;
1171 };
1172
1173 return variableWidthBuffer( segments, interpolateWidths );
1174 }
1175
variableWidthBufferByM(int segments) const1176 QgsGeometry QgsInternalGeometryEngine::variableWidthBufferByM( int segments ) const
1177 {
1178 mLastError.clear();
1179 auto widthByM = []( const QgsLineString * line )->std::unique_ptr< double [] >
1180 {
1181 std::unique_ptr< double [] > widths( new double[ line->nCoordinates() ] );
1182 for ( int i = 0; i < line->nCoordinates(); ++i )
1183 {
1184 widths[ i ] = line->mAt( i );
1185 }
1186 return widths;
1187 };
1188
1189 return variableWidthBuffer( segments, widthByM );
1190 }
1191
randomPointsInPolygon(const QgsGeometry & polygon,int count,const std::function<bool (const QgsPointXY &)> & acceptPoint,unsigned long seed,QgsFeedback * feedback,int maxTriesPerPoint)1192 QVector<QgsPointXY> QgsInternalGeometryEngine::randomPointsInPolygon( const QgsGeometry &polygon, int count,
1193 const std::function< bool( const QgsPointXY & ) > &acceptPoint, unsigned long seed, QgsFeedback *feedback, int maxTriesPerPoint )
1194 {
1195 if ( polygon.type() != QgsWkbTypes::PolygonGeometry || count == 0 )
1196 return QVector< QgsPointXY >();
1197
1198 // step 1 - tessellate the polygon to triangles
1199 QgsRectangle bounds = polygon.boundingBox();
1200 QgsTessellator t( bounds, false, false, false, true );
1201
1202 if ( polygon.isMultipart() )
1203 {
1204 const QgsMultiSurface *ms = qgsgeometry_cast< const QgsMultiSurface * >( polygon.constGet() );
1205 for ( int i = 0; i < ms->numGeometries(); ++i )
1206 {
1207 if ( feedback && feedback->isCanceled() )
1208 return QVector< QgsPointXY >();
1209
1210 if ( QgsPolygon *poly = qgsgeometry_cast< QgsPolygon * >( ms->geometryN( i ) ) )
1211 {
1212 t.addPolygon( *poly, 0 );
1213 }
1214 else
1215 {
1216 std::unique_ptr< QgsPolygon > p( qgsgeometry_cast< QgsPolygon * >( ms->geometryN( i )->segmentize() ) );
1217 t.addPolygon( *p, 0 );
1218 }
1219 }
1220 }
1221 else
1222 {
1223 if ( const QgsPolygon *poly = qgsgeometry_cast< const QgsPolygon * >( polygon.constGet() ) )
1224 {
1225 t.addPolygon( *poly, 0 );
1226 }
1227 else
1228 {
1229 std::unique_ptr< QgsPolygon > p( qgsgeometry_cast< QgsPolygon * >( polygon.constGet()->segmentize() ) );
1230 t.addPolygon( *p, 0 );
1231 }
1232 }
1233
1234 if ( feedback && feedback->isCanceled() )
1235 return QVector< QgsPointXY >();
1236
1237 const QVector<float> triangleData = t.data();
1238 if ( triangleData.empty() )
1239 return QVector< QgsPointXY >(); //hm
1240
1241 // calculate running sum of triangle areas
1242 std::vector< double > cumulativeAreas;
1243 cumulativeAreas.reserve( t.dataVerticesCount() / 3 );
1244 double totalArea = 0;
1245 for ( auto it = triangleData.constBegin(); it != triangleData.constEnd(); )
1246 {
1247 if ( feedback && feedback->isCanceled() )
1248 return QVector< QgsPointXY >();
1249
1250 const float aX = *it++;
1251 ( void )it++; // z
1252 const float aY = -( *it++ );
1253 const float bX = *it++;
1254 ( void )it++; // z
1255 const float bY = -( *it++ );
1256 const float cX = *it++;
1257 ( void )it++; // z
1258 const float cY = -( *it++ );
1259
1260 const double area = QgsGeometryUtils::triangleArea( aX, aY, bX, bY, cX, cY );
1261 totalArea += area;
1262 cumulativeAreas.emplace_back( totalArea );
1263 }
1264
1265 std::random_device rd;
1266 std::mt19937 mt( seed == 0 ? rd() : seed );
1267 std::uniform_real_distribution<> uniformDist( 0, 1 );
1268
1269 // selects a random triangle, weighted by triangle area
1270 auto selectRandomTriangle = [&cumulativeAreas, totalArea]( double random )->int
1271 {
1272 int triangle = 0;
1273 const double target = random * totalArea;
1274 for ( auto it = cumulativeAreas.begin(); it != cumulativeAreas.end(); ++it, triangle++ )
1275 {
1276 if ( *it > target )
1277 return triangle;
1278 }
1279 Q_ASSERT_X( false, "QgsInternalGeometryEngine::randomPointsInPolygon", "Invalid random triangle index" );
1280 return 0; // no warnings
1281 };
1282
1283
1284 QVector<QgsPointXY> result;
1285 result.reserve( count );
1286 int tries = 0;
1287 for ( int i = 0; i < count; )
1288 {
1289 if ( feedback && feedback->isCanceled() )
1290 return QVector< QgsPointXY >();
1291
1292 const double triangleIndexRnd = uniformDist( mt );
1293 // pick random triangle, weighted by triangle area
1294 const int triangleIndex = selectRandomTriangle( triangleIndexRnd );
1295
1296 // generate a random point inside this triangle
1297 const double weightB = uniformDist( mt );
1298 const double weightC = uniformDist( mt );
1299 double x;
1300 double y;
1301
1302 // get triangle
1303 const double aX = triangleData.at( triangleIndex * 9 ) + bounds.xMinimum();
1304 const double aY = -triangleData.at( triangleIndex * 9 + 2 ) + bounds.yMinimum();
1305 const double bX = triangleData.at( triangleIndex * 9 + 3 ) + bounds.xMinimum();
1306 const double bY = -triangleData.at( triangleIndex * 9 + 5 ) + bounds.yMinimum();
1307 const double cX = triangleData.at( triangleIndex * 9 + 6 ) + bounds.xMinimum();
1308 const double cY = -triangleData.at( triangleIndex * 9 + 8 ) + bounds.yMinimum();
1309 QgsGeometryUtils::weightedPointInTriangle( aX, aY, bX, bY, cX, cY, weightB, weightC, x, y );
1310
1311 QgsPointXY candidate( x, y );
1312 if ( acceptPoint( candidate ) )
1313 {
1314 result << QgsPointXY( x, y );
1315 i++;
1316 tries = 0;
1317 }
1318 else if ( maxTriesPerPoint != 0 )
1319 {
1320 tries++;
1321 // Skip this point if maximum tries is reached
1322 if ( tries == maxTriesPerPoint )
1323 {
1324 tries = 0;
1325 i++;
1326 }
1327 }
1328 }
1329 return result;
1330 }
1331
1332 // ported from PostGIS' lwgeom pta_unstroke
1333
lineToCurve(const QgsLineString * lineString,double distanceTolerance,double pointSpacingAngleTolerance)1334 std::unique_ptr< QgsCompoundCurve > lineToCurve( const QgsLineString *lineString, double distanceTolerance,
1335 double pointSpacingAngleTolerance )
1336 {
1337 std::unique_ptr< QgsCompoundCurve > out = std::make_unique< QgsCompoundCurve >();
1338
1339 /* Minimum number of edges, per quadrant, required to define an arc */
1340 const unsigned int minQuadEdges = 2;
1341
1342 /* Die on null input */
1343 if ( !lineString )
1344 return nullptr;
1345
1346 /* Null on empty input? */
1347 if ( lineString->nCoordinates() == 0 )
1348 return nullptr;
1349
1350 /* We can't desegmentize anything shorter than four points */
1351 if ( lineString->nCoordinates() < 4 )
1352 {
1353 out->addCurve( lineString->clone() );
1354 return out;
1355 }
1356
1357 /* Allocate our result array of vertices that are part of arcs */
1358 int numEdges = lineString->nCoordinates() - 1;
1359 QVector< int > edgesInArcs( numEdges + 1, 0 );
1360
1361 auto arcAngle = []( const QgsPoint & a, const QgsPoint & b, const QgsPoint & c )->double
1362 {
1363 double abX = b.x() - a.x();
1364 double abY = b.y() - a.y();
1365
1366 double cbX = b.x() - c.x();
1367 double cbY = b.y() - c.y();
1368
1369 double dot = ( abX * cbX + abY * cbY ); /* dot product */
1370 double cross = ( abX * cbY - abY * cbX ); /* cross product */
1371
1372 double alpha = std::atan2( cross, dot );
1373
1374 return alpha;
1375 };
1376
1377 /* We make a candidate arc of the first two edges, */
1378 /* And then see if the next edge follows it */
1379 int i = 0;
1380 int j = 0;
1381 int k = 0;
1382 int currentArc = 1;
1383 QgsPoint a1;
1384 QgsPoint a2;
1385 QgsPoint a3;
1386 QgsPoint b;
1387 double centerX = 0.0;
1388 double centerY = 0.0;
1389 double radius = 0;
1390
1391 while ( i < numEdges - 2 )
1392 {
1393 unsigned int arcEdges = 0;
1394 double numQuadrants = 0;
1395 double angle;
1396
1397 bool foundArc = false;
1398 /* Make candidate arc */
1399 a1 = lineString->pointN( i );
1400 a2 = lineString->pointN( i + 1 );
1401 a3 = lineString->pointN( i + 2 );
1402 QgsPoint first = a1;
1403
1404 for ( j = i + 3; j < numEdges + 1; j++ )
1405 {
1406 b = lineString->pointN( j );
1407
1408 /* Does this point fall on our candidate arc? */
1409 if ( QgsGeometryUtils::pointContinuesArc( a1, a2, a3, b, distanceTolerance, pointSpacingAngleTolerance ) )
1410 {
1411 /* Yes. Mark this edge and the two preceding it as arc components */
1412 foundArc = true;
1413 for ( k = j - 1; k > j - 4; k-- )
1414 edgesInArcs[k] = currentArc;
1415 }
1416 else
1417 {
1418 /* No. So we're done with this candidate arc */
1419 currentArc++;
1420 break;
1421 }
1422
1423 a1 = a2;
1424 a2 = a3;
1425 a3 = b;
1426 }
1427 /* Jump past all the edges that were added to the arc */
1428 if ( foundArc )
1429 {
1430 /* Check if an arc was composed by enough edges to be
1431 * really considered an arc
1432 * See http://trac.osgeo.org/postgis/ticket/2420
1433 */
1434 arcEdges = j - 1 - i;
1435 if ( first.x() == b.x() && first.y() == b.y() )
1436 {
1437 numQuadrants = 4;
1438 }
1439 else
1440 {
1441 QgsGeometryUtils::circleCenterRadius( first, b, a1, radius, centerX, centerY );
1442
1443 angle = arcAngle( first, QgsPoint( centerX, centerY ), b );
1444 int p2Side = QgsGeometryUtils::leftOfLine( b.x(), b.y(), first.x(), first.y(), a1.x(), a1.y() );
1445 if ( p2Side >= 0 )
1446 angle = -angle;
1447
1448 if ( angle < 0 )
1449 angle = 2 * M_PI + angle;
1450 numQuadrants = ( 4 * angle ) / ( 2 * M_PI );
1451 }
1452 /* a1 is first point, b is last point */
1453 if ( arcEdges < minQuadEdges * numQuadrants )
1454 {
1455 // LWDEBUGF( 4, "Not enough edges for a %g quadrants arc, %g needed", num_quadrants, min_quad_edges * num_quadrants );
1456 for ( k = j - 1; k >= i; k-- )
1457 edgesInArcs[k] = 0;
1458 }
1459
1460 i = j - 1;
1461 }
1462 else
1463 {
1464 /* Mark this edge as a linear edge */
1465 edgesInArcs[i] = 0;
1466 i = i + 1;
1467 }
1468 }
1469
1470 int start = 0;
1471 int end = 0;
1472 /* non-zero if edge is part of an arc */
1473 int edgeType = edgesInArcs[0];
1474
1475 auto addPointsToCurve = [ lineString, &out ]( int start, int end, int type )
1476 {
1477 if ( type == 0 )
1478 {
1479 // straight segment
1480 QVector< QgsPoint > points;
1481 for ( int j = start; j < end + 2; ++ j )
1482 {
1483 points.append( lineString->pointN( j ) );
1484 }
1485 std::unique_ptr< QgsCurve > straightSegment = std::make_unique< QgsLineString >( points );
1486 out->addCurve( straightSegment.release() );
1487 }
1488 else
1489 {
1490 // curved segment
1491 QVector< QgsPoint > points;
1492 points.append( lineString->pointN( start ) );
1493 points.append( lineString->pointN( ( start + end + 1 ) / 2 ) );
1494 points.append( lineString->pointN( end + 1 ) );
1495 std::unique_ptr< QgsCircularString > curvedSegment = std::make_unique< QgsCircularString >();
1496 curvedSegment->setPoints( points );
1497 out->addCurve( curvedSegment.release() );
1498 }
1499 };
1500
1501 for ( int i = 1; i < numEdges; i++ )
1502 {
1503 if ( edgeType != edgesInArcs[i] )
1504 {
1505 end = i - 1;
1506 addPointsToCurve( start, end, edgeType );
1507 start = i;
1508 edgeType = edgesInArcs[i];
1509 }
1510 }
1511
1512 /* Roll out last item */
1513 end = numEdges - 1;
1514 addPointsToCurve( start, end, edgeType );
1515
1516 return out;
1517 }
1518
convertGeometryToCurves(const QgsAbstractGeometry * geom,double distanceTolerance,double angleTolerance)1519 std::unique_ptr< QgsAbstractGeometry > convertGeometryToCurves( const QgsAbstractGeometry *geom, double distanceTolerance, double angleTolerance )
1520 {
1521 if ( QgsWkbTypes::geometryType( geom->wkbType() ) == QgsWkbTypes::LineGeometry )
1522 {
1523 return lineToCurve( static_cast< const QgsLineString * >( geom ), distanceTolerance, angleTolerance );
1524 }
1525 else
1526 {
1527 // polygon
1528 const QgsPolygon *polygon = static_cast< const QgsPolygon * >( geom );
1529 std::unique_ptr< QgsCurvePolygon > result = std::make_unique< QgsCurvePolygon>();
1530
1531 result->setExteriorRing( lineToCurve( static_cast< const QgsLineString * >( polygon->exteriorRing() ),
1532 distanceTolerance, angleTolerance ).release() );
1533 for ( int i = 0; i < polygon->numInteriorRings(); ++i )
1534 {
1535 result->addInteriorRing( lineToCurve( static_cast< const QgsLineString * >( polygon->interiorRing( i ) ),
1536 distanceTolerance, angleTolerance ).release() );
1537 }
1538
1539 return result;
1540 }
1541 }
1542
convertToCurves(double distanceTolerance,double angleTolerance) const1543 QgsGeometry QgsInternalGeometryEngine::convertToCurves( double distanceTolerance, double angleTolerance ) const
1544 {
1545 mLastError.clear();
1546 if ( !mGeometry )
1547 {
1548 return QgsGeometry();
1549 }
1550
1551 if ( QgsWkbTypes::geometryType( mGeometry->wkbType() ) == QgsWkbTypes::PointGeometry )
1552 {
1553 return QgsGeometry( mGeometry->clone() ); // point geometry, nothing to do
1554 }
1555
1556 if ( const QgsGeometryCollection *gc = qgsgeometry_cast< const QgsGeometryCollection *>( mGeometry ) )
1557 {
1558 int numGeom = gc->numGeometries();
1559 QVector< QgsAbstractGeometry * > geometryList;
1560 geometryList.reserve( numGeom );
1561 for ( int i = 0; i < numGeom; ++i )
1562 {
1563 geometryList << convertGeometryToCurves( gc->geometryN( i ), distanceTolerance, angleTolerance ).release();
1564 }
1565
1566 QgsGeometry first = QgsGeometry( geometryList.takeAt( 0 ) );
1567 for ( QgsAbstractGeometry *g : std::as_const( geometryList ) )
1568 {
1569 first.addPart( g );
1570 }
1571 return first;
1572 }
1573 else
1574 {
1575 return QgsGeometry( convertGeometryToCurves( mGeometry, distanceTolerance, angleTolerance ) );
1576 }
1577 }
1578
orientedMinimumBoundingBox(double & area,double & angle,double & width,double & height) const1579 QgsGeometry QgsInternalGeometryEngine::orientedMinimumBoundingBox( double &area, double &angle, double &width, double &height ) const
1580 {
1581 mLastError.clear();
1582
1583 QgsRectangle minRect;
1584 area = std::numeric_limits<double>::max();
1585 angle = 0;
1586 width = std::numeric_limits<double>::max();
1587 height = std::numeric_limits<double>::max();
1588
1589 if ( !mGeometry || mGeometry->nCoordinates() < 2 )
1590 return QgsGeometry();
1591
1592 std::unique_ptr< QgsGeometryEngine >engine( QgsGeometry::createGeometryEngine( mGeometry ) );
1593 QString error;
1594 std::unique_ptr< QgsAbstractGeometry > hull( engine->convexHull( &mLastError ) );
1595 if ( !hull )
1596 return QgsGeometry();
1597
1598 QgsVertexId vertexId;
1599 QgsPoint pt0;
1600 QgsPoint pt1;
1601 QgsPoint pt2;
1602 // get first point
1603 hull->nextVertex( vertexId, pt0 );
1604 pt1 = pt0;
1605 double totalRotation = 0;
1606 while ( hull->nextVertex( vertexId, pt2 ) )
1607 {
1608 double currentAngle = QgsGeometryUtils::lineAngle( pt1.x(), pt1.y(), pt2.x(), pt2.y() );
1609 double rotateAngle = 180.0 / M_PI * currentAngle;
1610 totalRotation += rotateAngle;
1611
1612 QTransform t = QTransform::fromTranslate( pt0.x(), pt0.y() );
1613 t.rotate( rotateAngle );
1614 t.translate( -pt0.x(), -pt0.y() );
1615
1616 hull->transform( t );
1617
1618 QgsRectangle bounds = hull->boundingBox();
1619 double currentArea = bounds.width() * bounds.height();
1620 if ( currentArea < area )
1621 {
1622 minRect = bounds;
1623 area = currentArea;
1624 angle = totalRotation;
1625 width = bounds.width();
1626 height = bounds.height();
1627 }
1628
1629 pt1 = hull->vertexAt( vertexId );
1630 }
1631
1632 if ( width > height )
1633 {
1634 width = minRect.height();
1635 height = minRect.width();
1636 angle = angle + 90.0;
1637 }
1638
1639 QgsGeometry minBounds = QgsGeometry::fromRect( minRect );
1640 minBounds.rotate( angle, QgsPointXY( pt0.x(), pt0.y() ) );
1641
1642 // constrain angle to 0 - 180
1643 if ( angle > 180.0 )
1644 angle = std::fmod( angle, 180.0 );
1645
1646 return minBounds;
1647 }
1648