QGIS API Documentation 3.37.0-Master (fdefdf9c27f)
qgsinternalgeometryengine.cpp
Go to the documentation of this file.
1/***************************************************************************
2 qgsinternalgeometryengine.cpp - QgsInternalGeometryEngine
3
4 ---------------------
5 begin : 13.1.2016
6 copyright : (C) 2016 by Matthias Kuhn
8 ***************************************************************************
9 * *
10 * This program is free software; you can redistribute it and/or modify *
12 * the Free Software Foundation; either version 2 of the License, or *
13 * (at your option) any later version. *
14 * *
15 ***************************************************************************/
16
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 "qgstessellator.h"
29#include "qgsfeedback.h"
30#include "qgsgeometryengine.h"
31#include "qgsmultilinestring.h"
32#include "qgsgeos.h"
33#include <QTransform>
34#include <functional>
35#include <memory>
36#include <queue>
37#include <random>
38#include <geos_c.h>
39
41 : mGeometry( geometry.constGet() )
42{
43
44}
45
47{
48 return mLastError;
49}
50
51
52enum class Direction
53{
54 Up,
55 Right,
56 Down,
57 Left,
58 None
59};
60
66Direction getEdgeDirection( const QgsPoint &p1, const QgsPoint &p2, double maxDev )
67{
68 double dx = p2.x() - p1.x();
69 double dy = p2.y() - p1.y();
70 if ( ( dx == 0.0 ) && ( dy == 0.0 ) )
71 return Direction::None;
72 if ( fabs( dx ) >= fabs( dy ) )
73 {
74 double dev = fabs( dy ) / fabs( dx );
75 if ( dev > maxDev )
76 return Direction::None;
77 return dx > 0.0 ? Direction::Right : Direction::Left;
78 }
79 else
80 {
81 double dev = fabs( dx ) / fabs( dy );
82 if ( dev > maxDev )
83 return Direction::None;
84 return dy > 0.0 ? Direction::Up : Direction::Down;
85 }
86}
87
93std::pair<bool, std::array<Direction, 4>> getEdgeDirections( const QgsPolygon *g, double maxDev )
94{
95 std::pair<bool, std::array<Direction, 4>> ret = { false, { Direction::None, Direction::None, Direction::None, Direction::None } };
96 // The polygon might start in the middle of a side. Hence, we need a fifth
97 // direction to record the beginning of the side when we went around the
98 // polygon.
99 std::array<Direction, 5> dirs;
100
101 int idx = 0;
103 QgsAbstractGeometry::vertex_iterator current = previous;
104 ++current;
106 while ( current != end )
107 {
108 Direction dir = getEdgeDirection( *previous, *current, maxDev );
109 if ( dir == Direction::None )
110 return ret;
111 if ( idx == 0 )
112 {
113 dirs[0] = dir;
114 ++idx;
115 }
116 else if ( dir != dirs[idx - 1] )
117 {
118 if ( idx == 5 )
119 return ret;
120 dirs[idx] = dir;
121 ++idx;
122 }
123 previous = current;
124 ++current;
125 }
126 ret.first = ( idx == 5 ) ? ( dirs[0] == dirs[4] ) : ( idx == 4 );
127 std::copy( dirs.begin(), dirs.begin() + 4, ret.second.begin() );
128 return ret;
129}
130
131bool matchesOrientation( std::array<Direction, 4> dirs, std::array<Direction, 4> oriented )
132{
133 int idx = std::find( oriented.begin(), oriented.end(), dirs[0] ) - oriented.begin();
134 for ( int i = 1; i < 4; ++i )
135 {
136 if ( dirs[i] != oriented[( idx + i ) % 4] )
137 return false;
138 }
139 return true;
140}
141
145bool isClockwise( std::array<Direction, 4> dirs )
146{
147 const std::array<Direction, 4> cwdirs = { Direction::Up, Direction::Right, Direction::Down, Direction::Left };
148 return matchesOrientation( dirs, cwdirs );
149}
150
155bool isCounterClockwise( std::array<Direction, 4> dirs )
156{
157 const std::array<Direction, 4> ccwdirs = { Direction::Right, Direction::Up, Direction::Left, Direction::Down };
158 return matchesOrientation( dirs, ccwdirs );
159}
160
161
162bool QgsInternalGeometryEngine::isAxisParallelRectangle( double maximumDeviation, bool simpleRectanglesOnly ) const
163{
165 return false;
166
167 const QgsPolygon *polygon = qgsgeometry_cast< const QgsPolygon * >( mGeometry );
168 if ( !polygon->exteriorRing() || polygon->numInteriorRings() > 0 )
169 return false;
170
171 const int vertexCount = polygon->exteriorRing()->numPoints();
172 if ( vertexCount < 4 )
173 return false;
174 else if ( simpleRectanglesOnly && ( vertexCount != 5 || !polygon->exteriorRing()->isClosed() ) )
175 return false;
176
177 bool found4Dirs;
178 std::array<Direction, 4> dirs;
179 std::tie( found4Dirs, dirs ) = getEdgeDirections( polygon, std::tan( maximumDeviation * M_PI / 180 ) );
180
181 return found4Dirs && ( isCounterClockwise( dirs ) || isClockwise( dirs ) );
182}
183
184/***************************************************************************
185 * This class is considered CRITICAL and any change MUST be accompanied with
186 * full unit tests.
187 * See details in QEP #17
188 ****************************************************************************/
189
191{
192 mLastError.clear();
193 QVector<QgsLineString *> linesToProcess;
194
195 const QgsMultiCurve *multiCurve = qgsgeometry_cast< const QgsMultiCurve * >( mGeometry );
196 if ( multiCurve )
197 {
198 linesToProcess.reserve( multiCurve->partCount() );
199 for ( int i = 0; i < multiCurve->partCount(); ++i )
200 {
201 linesToProcess << static_cast<QgsLineString *>( multiCurve->geometryN( i )->clone() );
202 }
203 }
204
205 const QgsCurve *curve = qgsgeometry_cast< const QgsCurve * >( mGeometry );
206 if ( curve )
207 {
208 linesToProcess << static_cast<QgsLineString *>( curve->segmentize() );
209 }
210
211 std::unique_ptr<QgsMultiPolygon> multipolygon( linesToProcess.size() > 1 ? new QgsMultiPolygon() : nullptr );
212 QgsPolygon *polygon = nullptr;
213
214 if ( !linesToProcess.empty() )
215 {
216 std::unique_ptr< QgsLineString > secondline;
217 for ( QgsLineString *line : std::as_const( linesToProcess ) )
218 {
219 QTransform transform = QTransform::fromTranslate( x, y );
220
221 secondline.reset( line->reversed() );
222 secondline->transform( transform );
223
224 line->append( secondline.get() );
225 line->addVertex( line->pointN( 0 ) );
226
227 polygon = new QgsPolygon();
228 polygon->setExteriorRing( line );
229
230 if ( multipolygon )
232 }
233
234 if ( multipolygon )
235 return QgsGeometry( multipolygon.release() );
236 else
237 return QgsGeometry( polygon );
238 }
239
240 return QgsGeometry();
241}
242
243
244
245// polylabel implementation
246// ported from the original JavaScript implementation developed by Vladimir Agafonkin
248
250class Cell
251{
252 public:
253 Cell( double x, double y, double h, const QgsPolygon *polygon )
254 : x( x )
255 , y( y )
256 , h( h )
257 , d( polygon->pointDistanceToBoundary( x, y ) )
258 , max( d + h * M_SQRT2 )
259 {}
260
262 double x;
264 double y;
266 double h;
268 double d;
270 double max;
271};
272
273struct GreaterThanByMax
274{
275 bool operator()( const Cell *lhs, const Cell *rhs ) const
276 {
277 return rhs->max > lhs->max;
278 }
279};
280
281Cell *getCentroidCell( const QgsPolygon *polygon )
282{
283 double area = 0;
284 double x = 0;
285 double y = 0;
286
287 const QgsLineString *exterior = static_cast< const QgsLineString *>( polygon->exteriorRing() );
288 int len = exterior->numPoints() - 1; //assume closed
289 for ( int i = 0, j = len - 1; i < len; j = i++ )
290 {
291 double aX = exterior->xAt( i );
292 double aY = exterior->yAt( i );
293 double bX = exterior->xAt( j );
294 double bY = exterior->yAt( j );
295 double f = aX * bY - bX * aY;
296 x += ( aX + bX ) * f;
297 y += ( aY + bY ) * f;
298 area += f * 3;
299 }
300 if ( area == 0 )
301 return new Cell( exterior->xAt( 0 ), exterior->yAt( 0 ), 0, polygon );
302 else
303 return new Cell( x / area, y / area, 0.0, polygon );
304}
305
306QgsPoint surfacePoleOfInaccessibility( const QgsSurface *surface, double precision, double &distanceFromBoundary )
307{
308 std::unique_ptr< QgsPolygon > segmentizedPoly;
309 const QgsPolygon *polygon = qgsgeometry_cast< const QgsPolygon * >( surface );
310 if ( !polygon )
311 {
312 segmentizedPoly.reset( static_cast< QgsPolygon *>( surface->segmentize() ) );
313 polygon = segmentizedPoly.get();
314 }
315
317 QgsRectangle bounds = polygon->boundingBox();
318
319 // initial parameters
320 double cellSize = std::min( bounds.width(), bounds.height() );
321
322 if ( qgsDoubleNear( cellSize, 0.0 ) )
323 return QgsPoint( bounds.xMinimum(), bounds.yMinimum() );
324
325 double h = cellSize / 2.0;
326 std::priority_queue< Cell *, std::vector<Cell *>, GreaterThanByMax > cellQueue;
327
328 // cover polygon with initial cells
329 for ( double x = bounds.xMinimum(); x < bounds.xMaximum(); x += cellSize )
330 {
331 for ( double y = bounds.yMinimum(); y < bounds.yMaximum(); y += cellSize )
332 {
333 cellQueue.push( new Cell( x + h, y + h, h, polygon ) );
334 }
335 }
336
337 // take centroid as the first best guess
338 std::unique_ptr< Cell > bestCell( getCentroidCell( polygon ) );
339
340 // special case for rectangular polygons
341 std::unique_ptr< Cell > bboxCell( new Cell( bounds.xMinimum() + bounds.width() / 2.0,
342 bounds.yMinimum() + bounds.height() / 2.0,
343 0, polygon ) );
344 if ( bboxCell->d > bestCell->d )
345 {
346 bestCell = std::move( bboxCell );
347 }
348
349 while ( !cellQueue.empty() )
350 {
351 // pick the most promising cell from the queue
352 std::unique_ptr< Cell > cell( cellQueue.top() );
353 cellQueue.pop();
354 Cell *currentCell = cell.get();
355
356 // update the best cell if we found a better one
357 if ( currentCell->d > bestCell->d )
358 {
359 bestCell = std::move( cell );
360 }
361
362 // do not drill down further if there's no chance of a better solution
363 if ( currentCell->max - bestCell->d <= precision )
364 continue;
365
366 // split the cell into four cells
367 h = currentCell->h / 2.0;
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 cellQueue.push( new Cell( currentCell->x - h, currentCell->y + h, h, polygon ) );
371 cellQueue.push( new Cell( currentCell->x + h, currentCell->y + h, h, polygon ) );
372 }
373
374 distanceFromBoundary = bestCell->d;
375
376 return QgsPoint( bestCell->x, bestCell->y );
377}
378
380
381QgsGeometry QgsInternalGeometryEngine::poleOfInaccessibility( double precision, double *distanceFromBoundary ) const
382{
383 mLastError.clear();
384 if ( distanceFromBoundary )
385 *distanceFromBoundary = std::numeric_limits<double>::max();
386
387 if ( !mGeometry || mGeometry->isEmpty() )
388 return QgsGeometry();
389
390 if ( precision <= 0 )
391 return QgsGeometry();
392
393 if ( const QgsGeometryCollection *gc = qgsgeometry_cast< const QgsGeometryCollection *>( mGeometry ) )
394 {
395 int numGeom = gc->numGeometries();
396 double maxDist = 0;
397 QgsPoint bestPoint;
398 for ( int i = 0; i < numGeom; ++i )
399 {
400 const QgsSurface *surface = qgsgeometry_cast< const QgsSurface * >( gc->geometryN( i ) );
401 if ( !surface )
402 continue;
403
404 double dist = std::numeric_limits<double>::max();
405 QgsPoint p = surfacePoleOfInaccessibility( surface, precision, dist );
406 if ( dist > maxDist )
407 {
408 maxDist = dist;
409 bestPoint = p;
410 }
411 }
412
413 if ( bestPoint.isEmpty() )
414 return QgsGeometry();
415
416 if ( distanceFromBoundary )
417 *distanceFromBoundary = maxDist;
418 return QgsGeometry( new QgsPoint( bestPoint ) );
419 }
420 else
421 {
422 const QgsSurface *surface = qgsgeometry_cast< const QgsSurface * >( mGeometry );
423 if ( !surface )
424 return QgsGeometry();
425
426 double dist = std::numeric_limits<double>::max();
427 QgsPoint p = surfacePoleOfInaccessibility( surface, precision, dist );
428 if ( p.isEmpty() )
429 return QgsGeometry();
430
431 if ( distanceFromBoundary )
432 *distanceFromBoundary = dist;
433 return QgsGeometry( new QgsPoint( p ) );
434 }
435}
436
437
438// helpers for orthogonalize
439// adapted from original code in potlatch/id osm editor
440
441bool dotProductWithinAngleTolerance( double dotProduct, double lowerThreshold, double upperThreshold )
442{
443 return lowerThreshold > std::fabs( dotProduct ) || std::fabs( dotProduct ) > upperThreshold;
444}
445
446double normalizedDotProduct( const QgsPoint &a, const QgsPoint &b, const QgsPoint &c )
447{
448 QgsVector p = a - b;
449 QgsVector q = c - b;
450
451 if ( p.length() > 0 )
452 p = p.normalized();
453 if ( q.length() > 0 )
454 q = q.normalized();
455
456 return p * q;
457}
458
459double squareness( QgsLineString *ring, double lowerThreshold, double upperThreshold )
460{
461 double sum = 0.0;
462
463 bool isClosed = ring->isClosed();
464 int numPoints = ring->numPoints();
465 QgsPoint a;
466 QgsPoint b;
467 QgsPoint c;
468
469 for ( int i = 0; i < numPoints; ++i )
470 {
471 if ( !isClosed && i == numPoints - 1 )
472 break;
473 else if ( !isClosed && i == 0 )
474 {
475 b = ring->pointN( 0 );
476 c = ring->pointN( 1 );
477 }
478 else
479 {
480 if ( i == 0 )
481 {
482 a = ring->pointN( numPoints - 1 );
483 b = ring->pointN( 0 );
484 }
485 if ( i == numPoints - 1 )
486 c = ring->pointN( 0 );
487 else
488 c = ring->pointN( i + 1 );
489
490 double dotProduct = normalizedDotProduct( a, b, c );
491 if ( !dotProductWithinAngleTolerance( dotProduct, lowerThreshold, upperThreshold ) )
492 continue;
493
494 sum += 2.0 * std::min( std::fabs( dotProduct - 1.0 ), std::min( std::fabs( dotProduct ), std::fabs( dotProduct + 1 ) ) );
495 }
496 a = b;
497 b = c;
498 }
499
500 return sum;
501}
502
503QgsVector calcMotion( const QgsPoint &a, const QgsPoint &b, const QgsPoint &c,
504 double lowerThreshold, double upperThreshold )
505{
506 QgsVector p = a - b;
507 QgsVector q = c - b;
508
509 if ( qgsDoubleNear( p.length(), 0.0 ) || qgsDoubleNear( q.length(), 0.0 ) )
510 return QgsVector( 0, 0 );
511
512 // 2.0 is a magic number from the original JOSM source code
513 double scale = 2.0 * std::min( p.length(), q.length() );
514
515 p = p.normalized();
516 q = q.normalized();
517 double dotProduct = p * q;
518
519 if ( !dotProductWithinAngleTolerance( dotProduct, lowerThreshold, upperThreshold ) )
520 {
521 return QgsVector( 0, 0 );
522 }
523
524 // wonderful nasty hack which has survived through JOSM -> id -> QGIS
525 // to deal with almost-straight segments (angle is closer to 180 than to 90/270).
526 if ( dotProduct < -M_SQRT1_2 )
527 dotProduct += 1.0;
528
529 QgsVector new_v = p + q;
530 if ( qgsDoubleNear( new_v.length(), 0.0 ) )
531 {
532 return QgsVector( 0, 0 );
533 }
534 // 0.1 magic number from JOSM implementation - think this is to limit each iterative step
535 return new_v.normalized() * ( 0.1 * dotProduct * scale );
536}
537
538QgsLineString *doOrthogonalize( QgsLineString *ring, int iterations, double tolerance, double lowerThreshold, double upperThreshold )
539{
540 double minScore = std::numeric_limits<double>::max();
541
542 bool isClosed = ring->isClosed();
543 int numPoints = ring->numPoints();
544
545 std::unique_ptr< QgsLineString > best( ring->clone() );
546
547 QVector< QgsVector > /* yep */ motions;
548 motions.reserve( numPoints );
549
550 for ( int it = 0; it < iterations; ++it )
551 {
552 motions.resize( 0 ); // avoid re-allocations
553
554 // first loop through an calculate all motions
555 QgsPoint a;
556 QgsPoint b;
557 QgsPoint c;
558 for ( int i = 0; i < numPoints; ++i )
559 {
560 if ( isClosed && i == numPoints - 1 )
561 motions << motions.at( 0 ); //closed ring, so last point follows first point motion
562 else if ( !isClosed && ( i == 0 || i == numPoints - 1 ) )
563 {
564 b = ring->pointN( 0 );
565 c = ring->pointN( 1 );
566 motions << QgsVector( 0, 0 ); //non closed line, leave first/last vertex untouched
567 }
568 else
569 {
570 if ( i == 0 )
571 {
572 a = ring->pointN( numPoints - 1 );
573 b = ring->pointN( 0 );
574 }
575 if ( i == numPoints - 1 )
576 c = ring->pointN( 0 );
577 else
578 c = ring->pointN( i + 1 );
579
580 motions << calcMotion( a, b, c, lowerThreshold, upperThreshold );
581 }
582 a = b;
583 b = c;
584 }
585
586 // then apply them
587 for ( int i = 0; i < numPoints; ++i )
588 {
589 ring->setXAt( i, ring->xAt( i ) + motions.at( i ).x() );
590 ring->setYAt( i, ring->yAt( i ) + motions.at( i ).y() );
591 }
592
593 double newScore = squareness( ring, lowerThreshold, upperThreshold );
594 if ( newScore < minScore )
595 {
596 best.reset( ring->clone() );
597 minScore = newScore;
598 }
599
600 if ( minScore < tolerance )
601 break;
602 }
603
604 delete ring;
605
606 return best.release();
607}
608
609
610QgsAbstractGeometry *orthogonalizeGeom( const QgsAbstractGeometry *geom, int maxIterations, double tolerance, double lowerThreshold, double upperThreshold )
611{
612 std::unique_ptr< QgsAbstractGeometry > segmentizedCopy;
613 if ( QgsWkbTypes::isCurvedType( geom->wkbType() ) )
614 {
615 segmentizedCopy.reset( geom->segmentize() );
616 geom = segmentizedCopy.get();
617 }
618
620 {
621 return doOrthogonalize( static_cast< QgsLineString * >( geom->clone() ),
622 maxIterations, tolerance, lowerThreshold, upperThreshold );
623 }
624 else
625 {
626 // polygon
627 const QgsPolygon *polygon = static_cast< const QgsPolygon * >( geom );
628 QgsPolygon *result = new QgsPolygon();
629
630 result->setExteriorRing( doOrthogonalize( static_cast< QgsLineString * >( polygon->exteriorRing()->clone() ),
631 maxIterations, tolerance, lowerThreshold, upperThreshold ) );
632 for ( int i = 0; i < polygon->numInteriorRings(); ++i )
633 {
634 result->addInteriorRing( doOrthogonalize( static_cast< QgsLineString * >( polygon->interiorRing( i )->clone() ),
635 maxIterations, tolerance, lowerThreshold, upperThreshold ) );
636 }
637
638 return result;
639 }
640}
641
642QgsGeometry QgsInternalGeometryEngine::orthogonalize( double tolerance, int maxIterations, double angleThreshold ) const
643{
644 mLastError.clear();
645 if ( !mGeometry || ( QgsWkbTypes::geometryType( mGeometry->wkbType() ) != Qgis::GeometryType::Line
647 {
648 return QgsGeometry();
649 }
650
651 double lowerThreshold = std::cos( ( 90 - angleThreshold ) * M_PI / 180.00 );
652 double upperThreshold = std::cos( angleThreshold * M_PI / 180.0 );
653
654 if ( const QgsGeometryCollection *gc = qgsgeometry_cast< const QgsGeometryCollection *>( mGeometry ) )
655 {
656 int numGeom = gc->numGeometries();
657 QVector< QgsAbstractGeometry * > geometryList;
658 geometryList.reserve( numGeom );
659 for ( int i = 0; i < numGeom; ++i )
660 {
661 geometryList << orthogonalizeGeom( gc->geometryN( i ), maxIterations, tolerance, lowerThreshold, upperThreshold );
662 }
663
664 QgsGeometry first = QgsGeometry( geometryList.takeAt( 0 ) );
665 for ( QgsAbstractGeometry *g : std::as_const( geometryList ) )
666 {
668 }
669 return first;
670 }
671 else
672 {
673 return QgsGeometry( orthogonalizeGeom( mGeometry, maxIterations, tolerance, lowerThreshold, upperThreshold ) );
674 }
675}
676
677// if extraNodesPerSegment < 0, then use distance based mode
678QgsLineString *doDensify( const QgsLineString *ring, int extraNodesPerSegment = -1, double distance = 1 )
679{
680 QVector< double > outX;
681 QVector< double > outY;
682 QVector< double > outZ;
683 QVector< double > outM;
684 double multiplier = 1.0 / double( extraNodesPerSegment + 1 );
685
686 int nPoints = ring->numPoints();
687 outX.reserve( ( extraNodesPerSegment + 1 ) * nPoints );
688 outY.reserve( ( extraNodesPerSegment + 1 ) * nPoints );
689 bool withZ = ring->is3D();
690 if ( withZ )
691 outZ.reserve( ( extraNodesPerSegment + 1 ) * nPoints );
692 bool withM = ring->isMeasure();
693 if ( withM )
694 outM.reserve( ( extraNodesPerSegment + 1 ) * nPoints );
695 double x1 = 0;
696 double x2 = 0;
697 double y1 = 0;
698 double y2 = 0;
699 double z1 = 0;
700 double z2 = 0;
701 double m1 = 0;
702 double m2 = 0;
703 double xOut = 0;
704 double yOut = 0;
705 double zOut = 0;
706 double mOut = 0;
707 int extraNodesThisSegment = extraNodesPerSegment;
708 for ( int i = 0; i < nPoints - 1; ++i )
709 {
710 x1 = ring->xAt( i );
711 x2 = ring->xAt( i + 1 );
712 y1 = ring->yAt( i );
713 y2 = ring->yAt( i + 1 );
714 if ( withZ )
715 {
716 z1 = ring->zAt( i );
717 z2 = ring->zAt( i + 1 );
718 }
719 if ( withM )
720 {
721 m1 = ring->mAt( i );
722 m2 = ring->mAt( i + 1 );
723 }
724
725 outX << x1;
726 outY << y1;
727 if ( withZ )
728 outZ << z1;
729 if ( withM )
730 outM << m1;
731
732 if ( extraNodesPerSegment < 0 )
733 {
734 // distance mode
735 extraNodesThisSegment = std::floor( QgsGeometryUtilsBase::distance2D( x2, y2, x1, y1 ) / distance );
736 if ( extraNodesThisSegment >= 1 )
737 multiplier = 1.0 / ( extraNodesThisSegment + 1 );
738 }
739
740 for ( int j = 0; j < extraNodesThisSegment; ++j )
741 {
742 double delta = multiplier * ( j + 1 );
743 xOut = x1 + delta * ( x2 - x1 );
744 yOut = y1 + delta * ( y2 - y1 );
745 if ( withZ )
746 zOut = z1 + delta * ( z2 - z1 );
747 if ( withM )
748 mOut = m1 + delta * ( m2 - m1 );
749
750 outX << xOut;
751 outY << yOut;
752 if ( withZ )
753 outZ << zOut;
754 if ( withM )
755 outM << mOut;
756 }
757 }
758 outX << ring->xAt( nPoints - 1 );
759 outY << ring->yAt( nPoints - 1 );
760 if ( withZ )
761 outZ << ring->zAt( nPoints - 1 );
762 if ( withM )
763 outM << ring->mAt( nPoints - 1 );
764
765 QgsLineString *result = new QgsLineString( outX, outY, outZ, outM );
766 return result;
767}
768
769QgsAbstractGeometry *densifyGeometry( const QgsAbstractGeometry *geom, int extraNodesPerSegment = 1, double distance = 1 )
770{
771 std::unique_ptr< QgsAbstractGeometry > segmentizedCopy;
772 if ( QgsWkbTypes::isCurvedType( geom->wkbType() ) )
773 {
774 segmentizedCopy.reset( geom->segmentize() );
775 geom = segmentizedCopy.get();
776 }
777
779 {
780 return doDensify( static_cast< const QgsLineString * >( geom ), extraNodesPerSegment, distance );
781 }
782 else
783 {
784 // polygon
785 const QgsPolygon *polygon = static_cast< const QgsPolygon * >( geom );
786 QgsPolygon *result = new QgsPolygon();
787
788 result->setExteriorRing( doDensify( static_cast< const QgsLineString * >( polygon->exteriorRing() ),
789 extraNodesPerSegment, distance ) );
790 for ( int i = 0; i < polygon->numInteriorRings(); ++i )
791 {
792 result->addInteriorRing( doDensify( static_cast< const QgsLineString * >( polygon->interiorRing( i ) ),
793 extraNodesPerSegment, distance ) );
794 }
795
796 return result;
797 }
798}
799
801{
802 mLastError.clear();
803 if ( !mGeometry )
804 {
805 return QgsGeometry();
806 }
807
809 {
810 return QgsGeometry( mGeometry->clone() ); // point geometry, nothing to do
811 }
812
813 if ( const QgsGeometryCollection *gc = qgsgeometry_cast< const QgsGeometryCollection *>( mGeometry ) )
814 {
815 int numGeom = gc->numGeometries();
816 QVector< QgsAbstractGeometry * > geometryList;
817 geometryList.reserve( numGeom );
818 for ( int i = 0; i < numGeom; ++i )
819 {
820 geometryList << densifyGeometry( gc->geometryN( i ), extraNodesPerSegment );
821 }
822
823 QgsGeometry first = QgsGeometry( geometryList.takeAt( 0 ) );
824 for ( QgsAbstractGeometry *g : std::as_const( geometryList ) )
825 {
827 }
828 return first;
829 }
830 else
831 {
832 return QgsGeometry( densifyGeometry( mGeometry, extraNodesPerSegment ) );
833 }
834}
835
837{
838 mLastError.clear();
839 if ( !mGeometry )
840 {
841 return QgsGeometry();
842 }
843
845 {
846 return QgsGeometry( mGeometry->clone() ); // point geometry, nothing to do
847 }
848
849 if ( const QgsGeometryCollection *gc = qgsgeometry_cast< const QgsGeometryCollection *>( mGeometry ) )
850 {
851 int numGeom = gc->numGeometries();
852 QVector< QgsAbstractGeometry * > geometryList;
853 geometryList.reserve( numGeom );
854 for ( int i = 0; i < numGeom; ++i )
855 {
856 geometryList << densifyGeometry( gc->geometryN( i ), -1, distance );
857 }
858
859 QgsGeometry first = QgsGeometry( geometryList.takeAt( 0 ) );
860 for ( QgsAbstractGeometry *g : std::as_const( geometryList ) )
861 {
863 }
864 return first;
865 }
866 else
867 {
868 return QgsGeometry( densifyGeometry( mGeometry, -1, distance ) );
869 }
870}
871
873//
874// QgsLineSegmentDistanceComparer
875//
876
877// adapted for QGIS geometry classes from original work at https://github.com/trylock/visibility by trylock
878bool QgsLineSegmentDistanceComparer::operator()( QgsLineSegment2D ab, QgsLineSegment2D cd ) const
879{
880 Q_ASSERT_X( ab.pointLeftOfLine( mOrigin ) != 0,
881 "line_segment_dist_comparer",
882 "AB must not be collinear with the origin." );
883 Q_ASSERT_X( cd.pointLeftOfLine( mOrigin ) != 0,
884 "line_segment_dist_comparer",
885 "CD must not be collinear with the origin." );
886
887 // flip the segments so that if there are common endpoints,
888 // they will be the segment's start points
889 // cppcheck-suppress mismatchingContainerExpression
890 if ( ab.end() == cd.start() || ab.end() == cd.end() )
891 ab.reverse();
892 // cppcheck-suppress mismatchingContainerExpression
893 if ( ab.start() == cd.end() )
894 cd.reverse();
895
896 // cases with common endpoints
897 if ( ab.start() == cd.start() )
898 {
899 const int oad = QgsGeometryUtilsBase::leftOfLine( cd.endX(), cd.endY(), mOrigin.x(), mOrigin.y(), ab.startX(), ab.startY() );
900 const int oab = ab.pointLeftOfLine( mOrigin );
901 // cppcheck-suppress mismatchingContainerExpression
902 if ( ab.end() == cd.end() || oad != oab )
903 return false;
904 else
905 return ab.pointLeftOfLine( cd.end() ) != oab;
906 }
907 else
908 {
909 // cases without common endpoints
910 const int cda = cd.pointLeftOfLine( ab.start() );
911 const int cdb = cd.pointLeftOfLine( ab.end() );
912 if ( cdb == 0 && cda == 0 )
913 {
914 return mOrigin.sqrDist( ab.start() ) < mOrigin.sqrDist( cd.start() );
915 }
916 else if ( cda == cdb || cda == 0 || cdb == 0 )
917 {
918 const int cdo = cd.pointLeftOfLine( mOrigin );
919 return cdo == cda || cdo == cdb;
920 }
921 else
922 {
923 const int abo = ab.pointLeftOfLine( mOrigin );
924 return abo != ab.pointLeftOfLine( cd.start() );
925 }
926 }
927}
928
929//
930// QgsClockwiseAngleComparer
931//
932
933bool QgsClockwiseAngleComparer::operator()( const QgsPointXY &a, const QgsPointXY &b ) const
934{
935 const bool aIsLeft = a.x() < mVertex.x();
936 const bool bIsLeft = b.x() < mVertex.x();
937 if ( aIsLeft != bIsLeft )
938 return bIsLeft;
939
940 if ( qgsDoubleNear( a.x(), mVertex.x() ) && qgsDoubleNear( b.x(), mVertex.x() ) )
941 {
942 if ( a.y() >= mVertex.y() || b.y() >= mVertex.y() )
943 {
944 return b.y() < a.y();
945 }
946 else
947 {
948 return a.y() < b.y();
949 }
950 }
951 else
952 {
953 const QgsVector oa = a - mVertex;
954 const QgsVector ob = b - mVertex;
955 const double det = oa.crossProduct( ob );
956 if ( qgsDoubleNear( det, 0.0 ) )
957 {
958 return oa.lengthSquared() < ob.lengthSquared();
959 }
960 else
961 {
962 return det < 0;
963 }
964 }
965}
966
968
969//
970// QgsRay2D
971//
972
973bool QgsRay2D::intersects( const QgsLineSegment2D &segment, QgsPointXY &intersectPoint ) const
974{
975 const QgsVector ao = origin - segment.start();
976 const QgsVector ab = segment.end() - segment.start();
977 const double det = ab.crossProduct( direction );
978 if ( qgsDoubleNear( det, 0.0 ) )
979 {
980 const int abo = segment.pointLeftOfLine( origin );
981 if ( abo != 0 )
982 {
983 return false;
984 }
985 else
986 {
987 const double distA = ao * direction;
988 const double distB = ( origin - segment.end() ) * direction;
989
990 if ( distA > 0 && distB > 0 )
991 {
992 return false;
993 }
994 else
995 {
996 if ( ( distA > 0 ) != ( distB > 0 ) )
997 intersectPoint = origin;
998 else if ( distA > distB ) // at this point, both distances are negative
999 intersectPoint = segment.start(); // hence the nearest point is A
1000 else
1001 intersectPoint = segment.end();
1002 return true;
1003 }
1004 }
1005 }
1006 else
1007 {
1008 const double u = ao.crossProduct( direction ) / det;
1009 if ( u < 0.0 || 1.0 < u )
1010 {
1011 return false;
1012 }
1013 else
1014 {
1015 const double t = -ab.crossProduct( ao ) / det;
1016 intersectPoint = origin + direction * t;
1017 return qgsDoubleNear( t, 0.0 ) || t > 0;
1018 }
1019 }
1020}
1021
1022QVector<QgsPointXY> generateSegmentCurve( const QgsPoint &center1, const double radius1, const QgsPoint &center2, const double radius2 )
1023{
1024 // ensure that first circle is smaller than second
1027
1028 QgsPointXY t1;
1029 QgsPointXY t2;
1030 QgsPointXY t3;
1031 QgsPointXY t4;
1032 QVector<QgsPointXY> points;
1033 if ( QgsGeometryUtils::circleCircleOuterTangents( center1, radius1, center2, radius2, t1, t2, t3, t4 ) )
1034 {
1035 points << t1
1036 << t2
1037 << t4
1038 << t3;
1039 }
1040 return points;
1041}
1042
1043// partially ported from JTS VariableWidthBuffer,
1044// https://github.com/topobyte/jts/blob/master/jts-lab/src/main/java/com/vividsolutions/jts/operation/buffer/VariableWidthBuffer.java
1045
1046QgsGeometry QgsInternalGeometryEngine::variableWidthBuffer( int segments, const std::function< std::unique_ptr< double[] >( const QgsLineString *line ) > &widthFunction ) const
1047{
1048 mLastError.clear();
1049 if ( !mGeometry )
1050 {
1051 return QgsGeometry();
1052 }
1053
1054 std::vector< std::unique_ptr<QgsLineString > > linesToProcess;
1055
1056 const QgsMultiCurve *multiCurve = qgsgeometry_cast< const QgsMultiCurve * >( mGeometry );
1057 if ( multiCurve )
1058 {
1059 for ( int i = 0; i < multiCurve->partCount(); ++i )
1060 {
1061 if ( static_cast< const QgsCurve * >( multiCurve->geometryN( i ) )->nCoordinates() == 0 )
1062 continue; // skip 0 length lines
1063
1064 linesToProcess.emplace_back( static_cast<QgsLineString *>( multiCurve->geometryN( i )->clone() ) );
1065 }
1066 }
1067
1068 const QgsCurve *curve = qgsgeometry_cast< const QgsCurve * >( mGeometry );
1069 if ( curve )
1070 {
1071 if ( curve->nCoordinates() > 0 )
1072 linesToProcess.emplace_back( static_cast<QgsLineString *>( curve->segmentize() ) );
1073 }
1074
1075 if ( linesToProcess.empty() )
1076 {
1077 QgsGeometry g;
1078 g.mLastError = QStringLiteral( "Input geometry was not a curve type geometry" );
1079 return g;
1080 }
1081
1082 QVector<QgsGeometry> bufferedLines;
1083 bufferedLines.reserve( linesToProcess.size() );
1084
1085 for ( std::unique_ptr< QgsLineString > &line : linesToProcess )
1086 {
1087 QVector<QgsGeometry> parts;
1088 QgsPoint prevPoint;
1090 QgsGeometry prevCircle;
1091
1092 std::unique_ptr< double[] > widths = widthFunction( line.get() ) ;
1093 for ( int i = 0; i < line->nCoordinates(); ++i )
1094 {
1095 QgsPoint thisPoint = line->pointN( i );
1096 QgsGeometry thisCircle;
1097 double thisRadius = widths[ i ] / 2.0;
1098 if ( thisRadius > 0 )
1099 {
1100 QgsGeometry p = QgsGeometry( thisPoint.clone() );
1101
1102 QgsCircle circ( thisPoint, thisRadius );
1103 thisCircle = QgsGeometry( circ.toPolygon( segments * 4 ) );
1104 parts << thisCircle;
1105 }
1106 else
1107 {
1108 thisCircle = QgsGeometry( thisPoint.clone() );
1109 }
1110
1111 if ( i > 0 )
1112 {
1114 {
1115 QVector< QgsPointXY > points = generateSegmentCurve( prevPoint, prevRadius, thisPoint, thisRadius );
1116 if ( !points.empty() )
1117 {
1118 // snap points to circle vertices
1119
1120 int atVertex = 0;
1121 int beforeVertex = 0;
1122 int afterVertex = 0;
1123 double sqrDist = 0;
1124 double sqrDistPrev = 0;
1125 for ( int j = 0; j < points.count(); ++j )
1126 {
1127 QgsPointXY pA = prevCircle.closestVertex( points.at( j ), atVertex, beforeVertex, afterVertex, sqrDistPrev );
1128 QgsPointXY pB = thisCircle.closestVertex( points.at( j ), atVertex, beforeVertex, afterVertex, sqrDist );
1129 points[j] = sqrDistPrev < sqrDist ? pA : pB;
1130 }
1131 // close ring
1132 points.append( points.at( 0 ) );
1133
1134 std::unique_ptr< QgsPolygon > poly = std::make_unique< QgsPolygon >();
1135 poly->setExteriorRing( new QgsLineString( points ) );
1136 if ( poly->area() > 0 )
1137 parts << QgsGeometry( std::move( poly ) );
1138 }
1139 }
1140 }
1141 prevPoint = thisPoint;
1143 prevCircle = thisCircle;
1144 }
1145
1146 bufferedLines << QgsGeometry::unaryUnion( parts );
1147 }
1148
1149 QgsGeometry res = QgsGeometry::collectGeometry( bufferedLines );
1150 // happens on some GEOS versions...
1152 return res;
1153}
1154
1155QgsGeometry QgsInternalGeometryEngine::taperedBuffer( double start, double end, int segments ) const
1156{
1157 mLastError.clear();
1158 start = std::fabs( start );
1159 end = std::fabs( end );
1160
1161 auto interpolateWidths = [ start, end ]( const QgsLineString * line )->std::unique_ptr< double [] >
1162 {
1163 // ported from JTS VariableWidthBuffer,
1164 // https://github.com/topobyte/jts/blob/master/jts-lab/src/main/java/com/vividsolutions/jts/operation/buffer/VariableWidthBuffer.java
1165 std::unique_ptr< double [] > widths( new double[ line->nCoordinates() ] );
1166 widths[0] = start;
1167 widths[line->nCoordinates() - 1] = end;
1168
1169 double lineLength = line->length();
1170 double currentLength = 0;
1171 QgsPoint prevPoint = line->pointN( 0 );
1172 for ( int i = 1; i < line->nCoordinates() - 1; ++i )
1173 {
1174 QgsPoint point = line->pointN( i );
1175 double segmentLength = point.distance( prevPoint );
1176 currentLength += segmentLength;
1177 double lengthFraction = lineLength > 0 ? currentLength / lineLength : 1;
1178 double delta = lengthFraction * ( end - start );
1179 widths[i] = start + delta;
1180 prevPoint = point;
1181 }
1182 return widths;
1183 };
1184
1185 return variableWidthBuffer( segments, interpolateWidths );
1186}
1187
1189{
1190 mLastError.clear();
1191 auto widthByM = []( const QgsLineString * line )->std::unique_ptr< double [] >
1192 {
1193 std::unique_ptr< double [] > widths( new double[ line->nCoordinates() ] );
1194 for ( int i = 0; i < line->nCoordinates(); ++i )
1195 {
1196 widths[ i ] = line->mAt( i );
1197 }
1198 return widths;
1199 };
1200
1201 return variableWidthBuffer( segments, widthByM );
1202}
1203
1204QVector<QgsPointXY> randomPointsInPolygonPoly2TriBackend( const QgsAbstractGeometry *geometry, int count,
1205 const std::function< bool( const QgsPointXY & ) > &acceptPoint, unsigned long seed, QgsFeedback *feedback, int maxTriesPerPoint, QString &error )
1206{
1207 // step 1 - tessellate the polygon to triangles
1208 QgsRectangle bounds = geometry->boundingBox();
1209 QgsTessellator t( bounds, false, false, false, true );
1210
1211 if ( const QgsMultiSurface *ms = qgsgeometry_cast< const QgsMultiSurface * >( geometry ) )
1212 {
1213 for ( int i = 0; i < ms->numGeometries(); ++i )
1214 {
1215 if ( feedback && feedback->isCanceled() )
1216 return QVector< QgsPointXY >();
1217
1218 if ( QgsPolygon *poly = qgsgeometry_cast< QgsPolygon * >( ms->geometryN( i ) ) )
1219 {
1221 }
1222 else
1223 {
1224 std::unique_ptr< QgsPolygon > p( qgsgeometry_cast< QgsPolygon * >( ms->geometryN( i )->segmentize() ) );
1226 }
1227 }
1228 }
1229 else
1230 {
1231 if ( const QgsPolygon *poly = qgsgeometry_cast< const QgsPolygon * >( geometry ) )
1232 {
1234 }
1235 else
1236 {
1237 std::unique_ptr< QgsPolygon > p( qgsgeometry_cast< QgsPolygon * >( geometry->segmentize() ) );
1239 }
1240 }
1241
1242 if ( feedback && feedback->isCanceled() )
1243 return QVector< QgsPointXY >();
1244
1245 if ( !t.error().isEmpty() )
1246 {
1247 error = t.error();
1248 return QVector< QgsPointXY >();
1249 }
1250
1251 const QVector<float> triangleData = t.data();
1252 if ( triangleData.empty() )
1253 return QVector< QgsPointXY >(); //hm
1254
1255 // calculate running sum of triangle areas
1256 std::vector< double > cumulativeAreas;
1257 cumulativeAreas.reserve( t.dataVerticesCount() / 3 );
1258 double totalArea = 0;
1259 for ( auto it = triangleData.constBegin(); it != triangleData.constEnd(); )
1260 {
1261 if ( feedback && feedback->isCanceled() )
1262 return QVector< QgsPointXY >();
1263
1264 const float aX = *it++;
1265 ( void )it++; // z
1266 const float aY = -( *it++ );
1267 const float bX = *it++;
1268 ( void )it++; // z
1269 const float bY = -( *it++ );
1270 const float cX = *it++;
1271 ( void )it++; // z
1272 const float cY = -( *it++ );
1273
1274 const double area = QgsGeometryUtilsBase::triangleArea( aX, aY, bX, bY, cX, cY );
1275 totalArea += area;
1276 cumulativeAreas.emplace_back( totalArea );
1277 }
1278
1279 std::random_device rd;
1280 std::mt19937 mt( seed == 0 ? rd() : seed );
1281 std::uniform_real_distribution<> uniformDist( 0, 1 );
1282
1283 // selects a random triangle, weighted by triangle area
1284 auto selectRandomTriangle = [&cumulativeAreas, totalArea]( double random )->int
1285 {
1286 int triangle = 0;
1287 const double target = random * totalArea;
1288 for ( auto it = cumulativeAreas.begin(); it != cumulativeAreas.end(); ++it, triangle++ )
1289 {
1290 if ( *it > target )
1291 return triangle;
1292 }
1293 Q_ASSERT_X( false, "QgsInternalGeometryEngine::randomPointsInPolygon", "Invalid random triangle index" );
1294 return 0; // no warnings
1295 };
1296
1297
1298 QVector<QgsPointXY> result;
1299 result.reserve( count );
1300 int tries = 0;
1301 for ( int i = 0; i < count; )
1302 {
1303 if ( feedback && feedback->isCanceled() )
1304 return QVector< QgsPointXY >();
1305
1306 const double triangleIndexRnd = uniformDist( mt );
1307 // pick random triangle, weighted by triangle area
1308 const int triangleIndex = selectRandomTriangle( triangleIndexRnd );
1309
1310 // generate a random point inside this triangle
1311 const double weightB = uniformDist( mt );
1312 const double weightC = uniformDist( mt );
1313 double x;
1314 double y;
1315
1316 // get triangle
1317 const double aX = triangleData.at( triangleIndex * 9 ) + bounds.xMinimum();
1318 const double aY = -triangleData.at( triangleIndex * 9 + 2 ) + bounds.yMinimum();
1319 const double bX = triangleData.at( triangleIndex * 9 + 3 ) + bounds.xMinimum();
1320 const double bY = -triangleData.at( triangleIndex * 9 + 5 ) + bounds.yMinimum();
1321 const double cX = triangleData.at( triangleIndex * 9 + 6 ) + bounds.xMinimum();
1322 const double cY = -triangleData.at( triangleIndex * 9 + 8 ) + bounds.yMinimum();
1323 QgsGeometryUtilsBase::weightedPointInTriangle( aX, aY, bX, bY, cX, cY, weightB, weightC, x, y );
1324
1325 QgsPointXY candidate( x, y );
1326 if ( acceptPoint( candidate ) )
1327 {
1328 result << QgsPointXY( x, y );
1329 i++;
1330 tries = 0;
1331 }
1332 else if ( maxTriesPerPoint != 0 )
1333 {
1334 tries++;
1335 // Skip this point if maximum tries is reached
1336 if ( tries == maxTriesPerPoint )
1337 {
1338 tries = 0;
1339 i++;
1340 }
1341 }
1342 }
1343 return result;
1344}
1345
1346QVector<QgsPointXY> randomPointsInPolygonGeosBackend( const QgsAbstractGeometry *geometry, int count,
1347 const std::function< bool( const QgsPointXY & ) > &acceptPoint, unsigned long seed, QgsFeedback *feedback, int maxTriesPerPoint, QString &error )
1348{
1349 // step 1 - tessellate the polygon to triangles
1350 QgsGeos geos( geometry );
1351 std::unique_ptr<QgsAbstractGeometry> triangulation = geos.constrainedDelaunayTriangulation( &error );
1352 if ( !triangulation || triangulation->isEmpty( ) )
1353 return {};
1354
1355 if ( feedback && feedback->isCanceled() )
1356 return {};
1357
1358 const QgsMultiPolygon *mp = qgsgeometry_cast< const QgsMultiPolygon * >( triangulation.get() );
1359 if ( !mp )
1360 return {};
1361
1362 // calculate running sum of triangle areas
1363 std::vector< double > cumulativeAreas;
1364 cumulativeAreas.reserve( mp->numGeometries() );
1365 double totalArea = 0;
1366 std::vector< double > vertices( static_cast< std::size_t >( mp->numGeometries() ) * 6 );
1367 double *vertexData = vertices.data();
1368 for ( auto it = mp->const_parts_begin(); it != mp->const_parts_end(); ++it )
1369 {
1370 if ( feedback && feedback->isCanceled() )
1371 return {};
1372
1373 const QgsPolygon *part = qgsgeometry_cast< const QgsPolygon * >( *it );
1374 if ( !part )
1375 return {};
1376
1377 const QgsLineString *exterior = qgsgeometry_cast< const QgsLineString * >( part->exteriorRing() );
1378 if ( !exterior )
1379 return {};
1380
1381 const double aX = exterior->xAt( 0 );
1382 const double aY = exterior->yAt( 0 );
1383 const double bX = exterior->xAt( 1 );
1384 const double bY = exterior->yAt( 1 );
1385 const double cX = exterior->xAt( 2 );
1386 const double cY = exterior->yAt( 2 );
1387
1388 const double area = QgsGeometryUtilsBase::triangleArea( aX, aY, bX, bY, cX, cY );
1389 *vertexData++ = aX;
1390 *vertexData++ = aY;
1391 *vertexData++ = bX;
1392 *vertexData++ = bY;
1393 *vertexData++ = cX;
1394 *vertexData++ = cY;
1395 totalArea += area;
1396 cumulativeAreas.emplace_back( totalArea );
1397 }
1398
1399 std::random_device rd;
1400 std::mt19937 mt( seed == 0 ? rd() : seed );
1401 std::uniform_real_distribution<> uniformDist( 0, 1 );
1402
1403 // selects a random triangle, weighted by triangle area
1404 auto selectRandomTriangle = [&cumulativeAreas, totalArea]( double random )->int
1405 {
1406 int triangle = 0;
1407 const double target = random * totalArea;
1408 for ( auto it = cumulativeAreas.begin(); it != cumulativeAreas.end(); ++it, triangle++ )
1409 {
1410 if ( *it > target )
1411 return triangle;
1412 }
1413 Q_ASSERT_X( false, "QgsInternalGeometryEngine::randomPointsInPolygon", "Invalid random triangle index" );
1414 return 0; // no warnings
1415 };
1416
1417
1418 QVector<QgsPointXY> result;
1419 result.reserve( count );
1420 int tries = 0;
1421 vertexData = vertices.data();
1422 for ( int i = 0; i < count; )
1423 {
1424 if ( feedback && feedback->isCanceled() )
1425 return QVector< QgsPointXY >();
1426
1427 const double triangleIndexRnd = uniformDist( mt );
1428 // pick random triangle, weighted by triangle area
1429 const std::size_t triangleIndex = selectRandomTriangle( triangleIndexRnd );
1430
1431 // generate a random point inside this triangle
1432 const double weightB = uniformDist( mt );
1433 const double weightC = uniformDist( mt );
1434 double x;
1435 double y;
1436
1437 // get triangle data
1438 const double aX = vertexData[ triangleIndex * 6 ];
1439 const double aY = vertexData[ triangleIndex * 6 + 1 ];
1440 const double bX = vertexData[ triangleIndex * 6 + 2 ];
1441 const double bY = vertexData[ triangleIndex * 6 + 3 ];
1442 const double cX = vertexData[ triangleIndex * 6 + 4 ];
1443 const double cY = vertexData[ triangleIndex * 6 + 5 ];
1444
1445 QgsGeometryUtilsBase::weightedPointInTriangle( aX, aY, bX, bY, cX, cY, weightB, weightC, x, y );
1446
1447 QgsPointXY candidate( x, y );
1448 if ( acceptPoint( candidate ) )
1449 {
1450 result << QgsPointXY( x, y );
1451 i++;
1452 tries = 0;
1453 }
1454 else if ( maxTriesPerPoint != 0 )
1455 {
1456 tries++;
1457 // Skip this point if maximum tries is reached
1458 if ( tries == maxTriesPerPoint )
1459 {
1460 tries = 0;
1461 i++;
1462 }
1463 }
1464 }
1465 return result;
1466}
1467
1469 const std::function< bool( const QgsPointXY & ) > &acceptPoint, unsigned long seed, QgsFeedback *feedback, int maxTriesPerPoint )
1470{
1471 if ( QgsWkbTypes::geometryType( mGeometry->wkbType() ) != Qgis::GeometryType::Polygon || count == 0 )
1472 return QVector< QgsPointXY >();
1473
1474 // prefer more stable GEOS implementation if available
1475#if (GEOS_VERSION_MAJOR==3 && GEOS_VERSION_MINOR>=11) || GEOS_VERSION_MAJOR>3
1476 return randomPointsInPolygonGeosBackend( mGeometry, count, acceptPoint, seed, feedback, maxTriesPerPoint, mLastError );
1477#else
1478 return randomPointsInPolygonPoly2TriBackend( mGeometry, count, acceptPoint, seed, feedback, maxTriesPerPoint, mLastError );
1479#endif
1480}
1481
1482// ported from PostGIS' lwgeom pta_unstroke
1483
1484std::unique_ptr< QgsCurve > lineToCurve( const QgsCurve *curve, double distanceTolerance,
1485 double pointSpacingAngleTolerance )
1486{
1487 const Qgis::WkbType flatType = QgsWkbTypes::flatType( curve->wkbType() );
1488 if ( flatType == Qgis::WkbType::CircularString )
1489 {
1490 // already curved, so just return a copy
1491 std::unique_ptr< QgsCircularString > out;
1492 out.reset( qgsgeometry_cast< QgsCircularString * >( curve )->clone() );
1493 return out;
1494 }
1495 else if ( flatType == Qgis::WkbType::CompoundCurve )
1496 {
1497 std::unique_ptr< QgsCompoundCurve > out = std::make_unique< QgsCompoundCurve >();
1498 const QgsCompoundCurve *in = qgsgeometry_cast< const QgsCompoundCurve * >( curve );
1499 for ( int i = 0; i < in->nCurves(); i ++ )
1500 {
1501 std::unique_ptr< QgsCurve > processed = lineToCurve( in->curveAt( i ), distanceTolerance, pointSpacingAngleTolerance );
1502 if ( processed )
1503 {
1504 if ( const QgsCompoundCurve *processedCompoundCurve = qgsgeometry_cast< const QgsCompoundCurve *>( processed.get() ) )
1505 {
1506 for ( int i = 0; i < processedCompoundCurve->nCurves(); ++i )
1507 {
1508 out->addCurve( processedCompoundCurve->curveAt( i )->clone() );
1509 }
1510 }
1511 else
1512 {
1514 }
1515 }
1516 }
1517 return out;
1518 }
1519 else if ( flatType == Qgis::WkbType::LineString )
1520 {
1521 const QgsLineString *lineString = qgsgeometry_cast< QgsLineString * >( curve );
1522
1523 std::unique_ptr< QgsCompoundCurve > out = std::make_unique< QgsCompoundCurve >();
1524
1525 /* Minimum number of edges, per quadrant, required to define an arc */
1526 const unsigned int minQuadEdges = 2;
1527
1528 /* Die on null input */
1529 if ( !lineString )
1530 return nullptr;
1531
1532 /* Null on empty input? */
1533 if ( lineString->nCoordinates() == 0 )
1534 return nullptr;
1535
1536 /* We can't desegmentize anything shorter than four points */
1537 if ( lineString->nCoordinates() < 4 )
1538 {
1540 return out;
1541 }
1542
1543 /* Allocate our result array of vertices that are part of arcs */
1544 int numEdges = lineString->nCoordinates() - 1;
1545 QVector< int > edgesInArcs( numEdges + 1, 0 );
1546
1547 auto arcAngle = []( const QgsPoint & a, const QgsPoint & b, const QgsPoint & c )->double
1548 {
1549 double abX = b.x() - a.x();
1550 double abY = b.y() - a.y();
1551
1552 double cbX = b.x() - c.x();
1553 double cbY = b.y() - c.y();
1554
1555 double dot = ( abX * cbX + abY * cbY ); /* dot product */
1556 double cross = ( abX * cbY - abY * cbX ); /* cross product */
1557
1558 double alpha = std::atan2( cross, dot );
1559
1560 return alpha;
1561 };
1562
1563 /* We make a candidate arc of the first two edges, */
1564 /* And then see if the next edge follows it */
1565 int i = 0;
1566 int j = 0;
1567 int k = 0;
1568 int currentArc = 1;
1569 QgsPoint a1;
1570 QgsPoint a2;
1571 QgsPoint a3;
1572 QgsPoint b;
1573 double centerX = 0.0;
1574 double centerY = 0.0;
1576
1577 while ( i < numEdges - 2 )
1578 {
1579 unsigned int arcEdges = 0;
1581 double angle;
1582
1583 bool foundArc = false;
1584 /* Make candidate arc */
1585 a1 = lineString->pointN( i );
1586 a2 = lineString->pointN( i + 1 );
1587 a3 = lineString->pointN( i + 2 );
1588 QgsPoint first = a1;
1589
1590 for ( j = i + 3; j < numEdges + 1; j++ )
1591 {
1592 b = lineString->pointN( j );
1593
1594 /* Does this point fall on our candidate arc? */
1595 if ( QgsGeometryUtils::pointContinuesArc( a1, a2, a3, b, distanceTolerance, pointSpacingAngleTolerance ) )
1596 {
1597 /* Yes. Mark this edge and the two preceding it as arc components */
1598 foundArc = true;
1599 for ( k = j - 1; k > j - 4; k-- )
1600 edgesInArcs[k] = currentArc;
1601 }
1602 else
1603 {
1604 /* No. So we're done with this candidate arc */
1605 currentArc++;
1606 break;
1607 }
1608
1609 a1 = a2;
1610 a2 = a3;
1611 a3 = b;
1612 }
1613 /* Jump past all the edges that were added to the arc */
1614 if ( foundArc )
1615 {
1616 /* Check if an arc was composed by enough edges to be
1617 * really considered an arc
1618 * See http://trac.osgeo.org/postgis/ticket/2420
1619 */
1620 arcEdges = j - 1 - i;
1621 if ( first.x() == b.x() && first.y() == b.y() )
1622 {
1624 }
1625 else
1626 {
1628
1629 angle = arcAngle( first, QgsPoint( centerX, centerY ), b );
1630 int p2Side = QgsGeometryUtilsBase::leftOfLine( b.x(), b.y(), first.x(), first.y(), a1.x(), a1.y() );
1631 if ( p2Side >= 0 )
1632 angle = -angle;
1633
1634 if ( angle < 0 )
1635 angle = 2 * M_PI + angle;
1636 numQuadrants = ( 4 * angle ) / ( 2 * M_PI );
1637 }
1638 /* a1 is first point, b is last point */
1640 {
1642 for ( k = j - 1; k >= i; k-- )
1643 edgesInArcs[k] = 0;
1644 }
1645
1646 i = j - 1;
1647 }
1648 else
1649 {
1650 /* Mark this edge as a linear edge */
1651 edgesInArcs[i] = 0;
1652 i = i + 1;
1653 }
1654 }
1655
1656 int start = 0;
1657 int end = 0;
1658 /* non-zero if edge is part of an arc */
1659 int edgeType = edgesInArcs[0];
1660
1661 auto addPointsToCurve = [ lineString, &out ]( int start, int end, int type )
1662 {
1663 if ( type == 0 )
1664 {
1665 // straight segment
1666 QVector< QgsPoint > points;
1667 for ( int j = start; j < end + 2; ++ j )
1668 {
1669 points.append( lineString->pointN( j ) );
1670 }
1671 std::unique_ptr< QgsCurve > straightSegment = std::make_unique< QgsLineString >( points );
1673 }
1674 else
1675 {
1676 // curved segment
1677 QVector< QgsPoint > points;
1678 points.append( lineString->pointN( start ) );
1679 points.append( lineString->pointN( ( start + end + 1 ) / 2 ) );
1680 points.append( lineString->pointN( end + 1 ) );
1681 std::unique_ptr< QgsCircularString > curvedSegment = std::make_unique< QgsCircularString >();
1682 curvedSegment->setPoints( points );
1684 }
1685 };
1686
1687 for ( int i = 1; i < numEdges; i++ )
1688 {
1689 if ( edgeType != edgesInArcs[i] )
1690 {
1691 end = i - 1;
1692 addPointsToCurve( start, end, edgeType );
1693 start = i;
1694 edgeType = edgesInArcs[i];
1695 }
1696 }
1697
1698 /* Roll out last item */
1699 end = numEdges - 1;
1700 addPointsToCurve( start, end, edgeType );
1701
1702 // return a simplified type if it doesn't need to be a compound curve resu
1703 if ( QgsWkbTypes::flatType( out->simplifiedTypeRef()->wkbType() ) == Qgis::WkbType::CircularString )
1704 {
1705 std::unique_ptr< QgsCircularString > res;
1706 res.reset( qgsgeometry_cast< const QgsCircularString * >( out->simplifiedTypeRef() )->clone() );
1707 return res;
1708 }
1709
1710 return out;
1711 }
1712
1713 return nullptr;
1714}
1715
1716std::unique_ptr< QgsAbstractGeometry > convertGeometryToCurves( const QgsAbstractGeometry *geom, double distanceTolerance, double angleTolerance )
1717{
1719 {
1720 return lineToCurve( qgsgeometry_cast< const QgsCurve * >( geom ), distanceTolerance, angleTolerance );
1721 }
1722 else
1723 {
1724 // polygon
1725 const QgsCurvePolygon *polygon = qgsgeometry_cast< const QgsCurvePolygon * >( geom );
1726 std::unique_ptr< QgsCurvePolygon > result = std::make_unique< QgsCurvePolygon>();
1727
1728 result->setExteriorRing( lineToCurve( polygon->exteriorRing(),
1729 distanceTolerance, angleTolerance ).release() );
1730 for ( int i = 0; i < polygon->numInteriorRings(); ++i )
1731 {
1732 result->addInteriorRing( lineToCurve( polygon->interiorRing( i ),
1733 distanceTolerance, angleTolerance ).release() );
1734 }
1735
1736 return result;
1737 }
1738}
1739
1740QgsGeometry QgsInternalGeometryEngine::convertToCurves( double distanceTolerance, double angleTolerance ) const
1741{
1742 mLastError.clear();
1743 if ( !mGeometry )
1744 {
1745 return QgsGeometry();
1746 }
1747
1749 {
1750 return QgsGeometry( mGeometry->clone() ); // point geometry, nothing to do
1751 }
1752
1753 if ( const QgsGeometryCollection *gc = qgsgeometry_cast< const QgsGeometryCollection *>( mGeometry ) )
1754 {
1755 int numGeom = gc->numGeometries();
1756 QVector< QgsAbstractGeometry * > geometryList;
1757 geometryList.reserve( numGeom );
1758 for ( int i = 0; i < numGeom; ++i )
1759 {
1760 geometryList << convertGeometryToCurves( gc->geometryN( i ), distanceTolerance, angleTolerance ).release();
1761 }
1762
1763 QgsGeometry first = QgsGeometry( geometryList.takeAt( 0 ) );
1764 for ( QgsAbstractGeometry *g : std::as_const( geometryList ) )
1765 {
1767 }
1768 return first;
1769 }
1770 else
1771 {
1772 return QgsGeometry( convertGeometryToCurves( mGeometry, distanceTolerance, angleTolerance ) );
1773 }
1774}
1775
1776QgsGeometry QgsInternalGeometryEngine::orientedMinimumBoundingBox( double &area, double &angle, double &width, double &height ) const
1777{
1778 mLastError.clear();
1779
1780 QgsRectangle minRect;
1781 area = std::numeric_limits<double>::max();
1782 angle = 0;
1783 width = std::numeric_limits<double>::max();
1784 height = std::numeric_limits<double>::max();
1785
1786 if ( !mGeometry || mGeometry->nCoordinates() < 2 )
1787 return QgsGeometry();
1788
1789 std::unique_ptr< QgsGeometryEngine >engine( QgsGeometry::createGeometryEngine( mGeometry ) );
1790 std::unique_ptr< QgsAbstractGeometry > hull( engine->convexHull( &mLastError ) );
1791 if ( !hull )
1792 return QgsGeometry();
1793
1794 QgsVertexId vertexId;
1795 QgsPoint pt0;
1796 QgsPoint pt1;
1797 QgsPoint pt2;
1798 // get first point
1799 hull->nextVertex( vertexId, pt0 );
1800 pt1 = pt0;
1801 double totalRotation = 0;
1802 while ( hull->nextVertex( vertexId, pt2 ) )
1803 {
1804 double currentAngle = QgsGeometryUtilsBase::lineAngle( pt1.x(), pt1.y(), pt2.x(), pt2.y() );
1805 double rotateAngle = 180.0 / M_PI * currentAngle;
1806 totalRotation += rotateAngle;
1807
1808 QTransform t = QTransform::fromTranslate( pt0.x(), pt0.y() );
1809 t.rotate( rotateAngle );
1810 t.translate( -pt0.x(), -pt0.y() );
1811
1812 hull->transform( t );
1813
1814 QgsRectangle bounds = hull->boundingBox();
1815 double currentArea = bounds.width() * bounds.height();
1816 if ( currentArea < area )
1817 {
1818 minRect = bounds;
1819 area = currentArea;
1820 angle = totalRotation;
1821 width = bounds.width();
1822 height = bounds.height();
1823 }
1824
1825 pt1 = hull->vertexAt( vertexId );
1826 }
1827
1828 QgsGeometry minBounds = QgsGeometry::fromRect( minRect );
1829 minBounds.rotate( angle, QgsPointXY( pt0.x(), pt0.y() ) );
1830
1831 if ( width > height )
1832 {
1833 width = minRect.height();
1834 height = minRect.width();
1835 angle = angle + 90.0;
1836 }
1837
1838 // constrain angle to 0 - 180
1839 if ( angle > 180.0 )
1840 angle = std::fmod( angle, 180.0 );
1841
1842 return minBounds;
1843}
1844
1845std::unique_ptr< QgsLineString > triangularWavesAlongLine( const QgsLineString *line, double wavelength, const double amplitude, const bool strictWavelength )
1846{
1847 const int totalPoints = line->numPoints();
1848 if ( totalPoints < 2 )
1849 return nullptr;
1850
1851 const double *x = line->xData();
1852 const double *y = line->yData();
1853
1854 double prevX = *x++;
1855 double prevY = *y++;
1856
1857 QVector< double > outX;
1858 QVector< double > outY;
1859 const double totalLength = line->length();
1860
1861 const int maxRepetitions = std::ceil( totalLength / wavelength );
1862 if ( !strictWavelength )
1863 wavelength = totalLength / maxRepetitions;
1864
1865 const int estimatedPoints = maxRepetitions * 2 + 2;
1866 outX.reserve( estimatedPoints );
1867 outY.reserve( estimatedPoints );
1868 outX.append( prevX );
1869 outY.append( prevY );
1870
1871 double distanceToNextPointFromStartOfSegment = wavelength / 4;
1872 int side = -1;
1873 for ( int i = 1; i < totalPoints; ++i )
1874 {
1875 double thisX = *x++;
1876 double thisY = *y++;
1877
1878 const double segmentAngleRadians = QgsGeometryUtilsBase::lineAngle( prevX, prevY, thisX, thisY );
1879 const double segmentLength = QgsGeometryUtilsBase::distance2D( thisX, thisY, prevX, prevY );
1880 while ( distanceToNextPointFromStartOfSegment < segmentLength || qgsDoubleNear( distanceToNextPointFromStartOfSegment, segmentLength ) )
1881 {
1882 // point falls on this segment - truncate to segment length if qgsDoubleNear test was actually > segment length
1883 const double distanceToPoint = std::min( distanceToNextPointFromStartOfSegment, segmentLength );
1884 double pX, pY;
1885 QgsGeometryUtilsBase::pointOnLineWithDistance( prevX, prevY, thisX, thisY, distanceToPoint, pX, pY );
1886
1887 // project point on line out by amplitude
1888 const double outPointX = pX + side * amplitude * std::sin( segmentAngleRadians + M_PI_2 );
1889 const double outPointY = pY + side * amplitude * std::cos( segmentAngleRadians + M_PI_2 );
1890
1891 outX.append( outPointX );
1892 outY.append( outPointY );
1893
1894 distanceToNextPointFromStartOfSegment += wavelength / 2;
1895 side = -side;
1896 }
1897
1898 prevX = thisX;
1899 prevY = thisY;
1900 distanceToNextPointFromStartOfSegment -= segmentLength;
1901 }
1902
1903 outX.append( prevX );
1904 outY.append( prevY );
1905
1906 return std::make_unique< QgsLineString >( outX, outY );
1907}
1908
1909std::unique_ptr< QgsLineString > triangularWavesRandomizedAlongLine( const QgsLineString *line,
1910 const double minimumWavelength, const double maximumWavelength,
1911 const double minimumAmplitude, const double maximumAmplitude,
1912 std::uniform_real_distribution<> &uniformDist, std::mt19937 &mt )
1913{
1914 const int totalPoints = line->numPoints();
1915 if ( totalPoints < 2 )
1916 return nullptr;
1917
1918 const double *x = line->xData();
1919 const double *y = line->yData();
1920
1921 double prevX = *x++;
1922 double prevY = *y++;
1923
1924 QVector< double > outX;
1925 QVector< double > outY;
1926 const double totalLength = line->length();
1927
1928 const int maxRepetitions = std::ceil( totalLength / minimumWavelength );
1929
1930 const int estimatedPoints = maxRepetitions * 2 + 2;
1931 outX.reserve( estimatedPoints );
1932 outY.reserve( estimatedPoints );
1933 outX.append( prevX );
1934 outY.append( prevY );
1935
1936 double wavelength = uniformDist( mt ) * ( maximumWavelength - minimumWavelength ) + minimumWavelength;
1937 double distanceToNextPointFromStartOfSegment = wavelength / 4;
1938 double amplitude = uniformDist( mt ) * ( maximumAmplitude - minimumAmplitude ) + minimumAmplitude;
1939
1940 int side = -1;
1941 for ( int i = 1; i < totalPoints; ++i )
1942 {
1943 double thisX = *x++;
1944 double thisY = *y++;
1945
1946 const double segmentAngleRadians = QgsGeometryUtilsBase::lineAngle( prevX, prevY, thisX, thisY );
1947 const double segmentLength = QgsGeometryUtilsBase::distance2D( thisX, thisY, prevX, prevY );
1948 while ( distanceToNextPointFromStartOfSegment < segmentLength || qgsDoubleNear( distanceToNextPointFromStartOfSegment, segmentLength ) )
1949 {
1950 // point falls on this segment - truncate to segment length if qgsDoubleNear test was actually > segment length
1951 const double distanceToPoint = std::min( distanceToNextPointFromStartOfSegment, segmentLength );
1952 double pX, pY;
1953 QgsGeometryUtilsBase::pointOnLineWithDistance( prevX, prevY, thisX, thisY, distanceToPoint, pX, pY );
1954
1955 // project point on line out by amplitude
1956 const double outPointX = pX + side * amplitude * std::sin( segmentAngleRadians + M_PI_2 );
1957 const double outPointY = pY + side * amplitude * std::cos( segmentAngleRadians + M_PI_2 );
1958 outX.append( outPointX );
1959 outY.append( outPointY );
1960
1961 wavelength = uniformDist( mt ) * ( maximumWavelength - minimumWavelength ) + minimumWavelength;
1962 amplitude = uniformDist( mt ) * ( maximumAmplitude - minimumAmplitude ) + minimumAmplitude;
1963
1964 distanceToNextPointFromStartOfSegment += wavelength / 2;
1965 side = -side;
1966 }
1967
1968 prevX = thisX;
1969 prevY = thisY;
1970 distanceToNextPointFromStartOfSegment -= segmentLength;
1971 }
1972
1973 outX.append( prevX );
1974 outY.append( prevY );
1975
1976 return std::make_unique< QgsLineString >( outX, outY );
1977}
1978
1979std::unique_ptr< QgsAbstractGeometry > triangularWavesPrivate( const QgsAbstractGeometry *geom, double wavelength, double amplitude, bool strictWavelength )
1980{
1981 std::unique_ptr< QgsAbstractGeometry > segmentizedCopy;
1982 if ( QgsWkbTypes::isCurvedType( geom->wkbType() ) )
1983 {
1984 segmentizedCopy.reset( geom->segmentize() );
1985 geom = segmentizedCopy.get();
1986 }
1987
1989 {
1990 return triangularWavesAlongLine( static_cast< const QgsLineString * >( geom ), wavelength, amplitude, strictWavelength );
1991 }
1992 else
1993 {
1994 // polygon
1995 const QgsPolygon *polygon = static_cast< const QgsPolygon * >( geom );
1996 std::unique_ptr< QgsPolygon > result = std::make_unique< QgsPolygon >();
1997
1998 result->setExteriorRing( triangularWavesAlongLine( static_cast< const QgsLineString * >( polygon->exteriorRing() ),
1999 wavelength, amplitude, strictWavelength ).release() );
2000 for ( int i = 0; i < polygon->numInteriorRings(); ++i )
2001 {
2002 result->addInteriorRing( triangularWavesAlongLine( static_cast< const QgsLineString * >( polygon->interiorRing( i ) ),
2003 wavelength, amplitude, strictWavelength ).release() );
2004 }
2005
2006 return result;
2007 }
2008}
2009
2010std::unique_ptr< QgsAbstractGeometry > triangularWavesRandomizedPrivate( const QgsAbstractGeometry *geom, double minimumWavelength, double maximumWavelength, double minimumAmplitude, double maximumAmplitude, std::uniform_real_distribution<> &uniformDist, std::mt19937 &mt )
2011{
2012 std::unique_ptr< QgsAbstractGeometry > segmentizedCopy;
2013 if ( QgsWkbTypes::isCurvedType( geom->wkbType() ) )
2014 {
2015 segmentizedCopy.reset( geom->segmentize() );
2016 geom = segmentizedCopy.get();
2017 }
2018
2020 {
2021 return triangularWavesRandomizedAlongLine( static_cast< const QgsLineString * >( geom ), minimumWavelength, maximumWavelength, minimumAmplitude, maximumAmplitude, uniformDist, mt );
2022 }
2023 else
2024 {
2025 // polygon
2026 const QgsPolygon *polygon = static_cast< const QgsPolygon * >( geom );
2027 std::unique_ptr< QgsPolygon > result = std::make_unique< QgsPolygon >();
2028
2029 result->setExteriorRing( triangularWavesRandomizedAlongLine( static_cast< const QgsLineString * >( polygon->exteriorRing() ),
2030 minimumWavelength, maximumWavelength, minimumAmplitude, maximumAmplitude, uniformDist, mt ).release() );
2031 for ( int i = 0; i < polygon->numInteriorRings(); ++i )
2032 {
2033 result->addInteriorRing( triangularWavesRandomizedAlongLine( static_cast< const QgsLineString * >( polygon->interiorRing( i ) ),
2034 minimumWavelength, maximumWavelength, minimumAmplitude, maximumAmplitude, uniformDist, mt ).release() );
2035 }
2036
2037 return result;
2038 }
2039}
2040
2041QgsGeometry QgsInternalGeometryEngine::triangularWaves( double wavelength, double amplitude, bool strictWavelength ) const
2042{
2043 if ( wavelength < 0 || qgsDoubleNear( wavelength, 0 ) )
2044 return QgsGeometry();
2045
2046 mLastError.clear();
2047 if ( !mGeometry )
2048 {
2049 return QgsGeometry();
2050 }
2051
2053 {
2054 return QgsGeometry( mGeometry->clone() ); // point geometry, nothing to do
2055 }
2056
2057 if ( const QgsGeometryCollection *gc = qgsgeometry_cast< const QgsGeometryCollection *>( mGeometry ) )
2058 {
2059 int numGeom = gc->numGeometries();
2060 QVector< QgsAbstractGeometry * > geometryList;
2061 geometryList.reserve( numGeom );
2062 for ( int i = 0; i < numGeom; ++i )
2063 {
2064 geometryList << triangularWavesPrivate( gc->geometryN( i ), wavelength, amplitude, strictWavelength ).release();
2065 }
2066
2067 QgsGeometry first = QgsGeometry( geometryList.takeAt( 0 ) );
2068 for ( QgsAbstractGeometry *g : std::as_const( geometryList ) )
2069 {
2071 }
2072 return first;
2073 }
2074 else
2075 {
2076 return QgsGeometry( triangularWavesPrivate( mGeometry, wavelength, amplitude, strictWavelength ) );
2077 }
2078}
2079
2080QgsGeometry QgsInternalGeometryEngine::triangularWavesRandomized( double minimumWavelength, double maximumWavelength, double minimumAmplitude, double maximumAmplitude, unsigned long seed ) const
2081{
2082 if ( minimumWavelength < 0 || qgsDoubleNear( minimumWavelength, 0 ) || maximumWavelength < 0 || qgsDoubleNear( maximumWavelength, 0 ) || maximumWavelength < minimumWavelength )
2083 return QgsGeometry();
2084
2085 mLastError.clear();
2086 if ( !mGeometry )
2087 {
2088 return QgsGeometry();
2089 }
2090
2091 std::random_device rd;
2092 std::mt19937 mt( seed == 0 ? rd() : seed );
2093 std::uniform_real_distribution<> uniformDist( 0, 1 );
2094
2096 {
2097 return QgsGeometry( mGeometry->clone() ); // point geometry, nothing to do
2098 }
2099
2100 if ( const QgsGeometryCollection *gc = qgsgeometry_cast< const QgsGeometryCollection *>( mGeometry ) )
2101 {
2102 int numGeom = gc->numGeometries();
2103 QVector< QgsAbstractGeometry * > geometryList;
2104 geometryList.reserve( numGeom );
2105 for ( int i = 0; i < numGeom; ++i )
2106 {
2107 geometryList << triangularWavesRandomizedPrivate( gc->geometryN( i ), minimumWavelength, maximumWavelength, minimumAmplitude, maximumAmplitude, uniformDist, mt ).release();
2108 }
2109
2110 QgsGeometry first = QgsGeometry( geometryList.takeAt( 0 ) );
2111 for ( QgsAbstractGeometry *g : std::as_const( geometryList ) )
2112 {
2114 }
2115 return first;
2116 }
2117 else
2118 {
2119 return QgsGeometry( triangularWavesRandomizedPrivate( mGeometry, minimumWavelength, maximumWavelength, minimumAmplitude, maximumAmplitude, uniformDist, mt ) );
2120 }
2121}
2122
2123std::unique_ptr< QgsLineString > squareWavesAlongLine( const QgsLineString *line, double wavelength, const double amplitude, const bool strictWavelength )
2124{
2125 const int totalPoints = line->numPoints();
2126 if ( totalPoints < 2 )
2127 return nullptr;
2128
2129 const double *x = line->xData();
2130 const double *y = line->yData();
2131
2132 double prevX = *x++;
2133 double prevY = *y++;
2134
2135 QVector< double > outX;
2136 QVector< double > outY;
2137 const double totalLength = line->length();
2138
2139 const int maxRepetitions = std::ceil( totalLength / wavelength );
2140 if ( !strictWavelength )
2141 wavelength = totalLength / maxRepetitions;
2142
2143 const int estimatedPoints = maxRepetitions * 4 + 2;
2144 outX.reserve( estimatedPoints );
2145 outY.reserve( estimatedPoints );
2146 outX.append( prevX );
2147 outY.append( prevY );
2148
2149 const double startSegmentAngleRadians = QgsGeometryUtilsBase::lineAngle( prevX, prevY, *x, *y );
2150 outX.append( prevX - amplitude * std::sin( startSegmentAngleRadians + M_PI_2 ) );
2151 outY.append( prevY - amplitude * std::cos( startSegmentAngleRadians + M_PI_2 ) );
2152
2153 double distanceToNextPointFromStartOfSegment = wavelength / 2;
2154
2155 int side = -1;
2156 for ( int i = 1; i < totalPoints; ++i )
2157 {
2158 double thisX = *x++;
2159 double thisY = *y++;
2160
2161 const double segmentAngleRadians = QgsGeometryUtilsBase::lineAngle( prevX, prevY, thisX, thisY );
2162 const double segmentLength = QgsGeometryUtilsBase::distance2D( thisX, thisY, prevX, prevY );
2163 while ( distanceToNextPointFromStartOfSegment < segmentLength || qgsDoubleNear( distanceToNextPointFromStartOfSegment, segmentLength ) )
2164 {
2165 // point falls on this segment - truncate to segment length if qgsDoubleNear test was actually > segment length
2166 const double distanceToPoint = std::min( distanceToNextPointFromStartOfSegment, segmentLength );
2167 double pX, pY;
2168 QgsGeometryUtilsBase::pointOnLineWithDistance( prevX, prevY, thisX, thisY, distanceToPoint, pX, pY );
2169
2170 // project point on line out by amplitude
2171 const double sinAngle = std::sin( segmentAngleRadians + M_PI_2 );
2172 const double cosAngle = std::cos( segmentAngleRadians + M_PI_2 );
2173 outX.append( pX + side * amplitude * sinAngle );
2174 outY.append( pY + side * amplitude * cosAngle );
2175 outX.append( pX - side * amplitude * sinAngle );
2176 outY.append( pY - side * amplitude * cosAngle );
2177
2178 distanceToNextPointFromStartOfSegment += wavelength / 2;
2179 side = -side;
2180 }
2181
2182 prevX = thisX;
2183 prevY = thisY;
2184 distanceToNextPointFromStartOfSegment -= segmentLength;
2185 }
2186
2187 // replace last point, which will be projected out to amplitude of wave, with the actual end point of the line
2188 outX.pop_back();
2189 outY.pop_back();
2190 outX.append( prevX );
2191 outY.append( prevY );
2192
2193 return std::make_unique< QgsLineString >( outX, outY );
2194}
2195
2196std::unique_ptr< QgsLineString > squareWavesRandomizedAlongLine( const QgsLineString *line,
2197 const double minimumWavelength, const double maximumWavelength,
2198 const double minimumAmplitude, const double maximumAmplitude,
2199 std::uniform_real_distribution<> &uniformDist, std::mt19937 &mt )
2200{
2201 const int totalPoints = line->numPoints();
2202 if ( totalPoints < 2 )
2203 return nullptr;
2204
2205 const double *x = line->xData();
2206 const double *y = line->yData();
2207
2208 double prevX = *x++;
2209 double prevY = *y++;
2210
2211 QVector< double > outX;
2212 QVector< double > outY;
2213 const double totalLength = line->length();
2214
2215 const int maxRepetitions = std::ceil( totalLength / minimumWavelength );
2216
2217 const int estimatedPoints = maxRepetitions * 4 + 2;
2218 outX.reserve( estimatedPoints );
2219 outY.reserve( estimatedPoints );
2220 outX.append( prevX );
2221 outY.append( prevY );
2222
2223 double amplitude = uniformDist( mt ) * ( maximumAmplitude - minimumAmplitude ) + minimumAmplitude;
2224
2225 double segmentAngleRadians = QgsGeometryUtilsBase::lineAngle( prevX, prevY, *x, *y );
2226 outX.append( prevX - amplitude * std::sin( segmentAngleRadians + M_PI_2 ) );
2227 outY.append( prevY - amplitude * std::cos( segmentAngleRadians + M_PI_2 ) );
2228
2229 double wavelength = uniformDist( mt ) * ( maximumWavelength - minimumWavelength ) + minimumWavelength;
2230 double distanceToNextPointFromStartOfSegment = wavelength / 2;
2231
2232 int side = -1;
2233 for ( int i = 1; i < totalPoints; ++i )
2234 {
2235 double thisX = *x++;
2236 double thisY = *y++;
2237
2238 segmentAngleRadians = QgsGeometryUtilsBase::lineAngle( prevX, prevY, thisX, thisY );
2239 const double segmentLength = QgsGeometryUtilsBase::distance2D( thisX, thisY, prevX, prevY );
2240 while ( distanceToNextPointFromStartOfSegment < segmentLength || qgsDoubleNear( distanceToNextPointFromStartOfSegment, segmentLength ) )
2241 {
2242 // point falls on this segment - truncate to segment length if qgsDoubleNear test was actually > segment length
2243 const double distanceToPoint = std::min( distanceToNextPointFromStartOfSegment, segmentLength );
2244 double pX, pY;
2245 QgsGeometryUtilsBase::pointOnLineWithDistance( prevX, prevY, thisX, thisY, distanceToPoint, pX, pY );
2246
2247 // project point on line out by amplitude
2248 const double sinAngle = std::sin( segmentAngleRadians + M_PI_2 );
2249 const double cosAngle = std::cos( segmentAngleRadians + M_PI_2 );
2250 outX.append( pX + side * amplitude * sinAngle );
2251 outY.append( pY + side * amplitude * cosAngle );
2252
2253 amplitude = uniformDist( mt ) * ( maximumAmplitude - minimumAmplitude ) + minimumAmplitude;
2254 outX.append( pX - side * amplitude * sinAngle );
2255 outY.append( pY - side * amplitude * cosAngle );
2256
2257 wavelength = uniformDist( mt ) * ( maximumWavelength - minimumWavelength ) + minimumWavelength;
2258 distanceToNextPointFromStartOfSegment += wavelength / 2;
2259 side = -side;
2260 }
2261
2262 prevX = thisX;
2263 prevY = thisY;
2264 distanceToNextPointFromStartOfSegment -= segmentLength;
2265 }
2266
2267 outX.append( prevX + side * amplitude * std::sin( segmentAngleRadians + M_PI_2 ) );
2268 outY.append( prevY + side * amplitude * std::cos( segmentAngleRadians + M_PI_2 ) );
2269 outX.append( prevX );
2270 outY.append( prevY );
2271
2272 return std::make_unique< QgsLineString >( outX, outY );
2273}
2274
2275std::unique_ptr< QgsAbstractGeometry > squareWavesPrivate( const QgsAbstractGeometry *geom, double wavelength, double amplitude, bool strictWavelength )
2276{
2277 std::unique_ptr< QgsAbstractGeometry > segmentizedCopy;
2278 if ( QgsWkbTypes::isCurvedType( geom->wkbType() ) )
2279 {
2280 segmentizedCopy.reset( geom->segmentize() );
2281 geom = segmentizedCopy.get();
2282 }
2283
2285 {
2286 return squareWavesAlongLine( static_cast< const QgsLineString * >( geom ), wavelength, amplitude, strictWavelength );
2287 }
2288 else
2289 {
2290 // polygon
2291 const QgsPolygon *polygon = static_cast< const QgsPolygon * >( geom );
2292 std::unique_ptr< QgsPolygon > result = std::make_unique< QgsPolygon >();
2293
2294 result->setExteriorRing( squareWavesAlongLine( static_cast< const QgsLineString * >( polygon->exteriorRing() ),
2295 wavelength, amplitude, strictWavelength ).release() );
2296 for ( int i = 0; i < polygon->numInteriorRings(); ++i )
2297 {
2298 result->addInteriorRing( squareWavesAlongLine( static_cast< const QgsLineString * >( polygon->interiorRing( i ) ),
2299 wavelength, amplitude, strictWavelength ).release() );
2300 }
2301
2302 return result;
2303 }
2304}
2305
2306std::unique_ptr< QgsAbstractGeometry > squareWavesRandomizedPrivate( const QgsAbstractGeometry *geom, double minimumWavelength, double maximumWavelength, double minimumAmplitude, double maximumAmplitude, std::uniform_real_distribution<> &uniformDist, std::mt19937 &mt )
2307{
2308 std::unique_ptr< QgsAbstractGeometry > segmentizedCopy;
2309 if ( QgsWkbTypes::isCurvedType( geom->wkbType() ) )
2310 {
2311 segmentizedCopy.reset( geom->segmentize() );
2312 geom = segmentizedCopy.get();
2313 }
2314
2316 {
2317 return squareWavesRandomizedAlongLine( static_cast< const QgsLineString * >( geom ), minimumWavelength, maximumWavelength, minimumAmplitude, maximumAmplitude, uniformDist, mt );
2318 }
2319 else
2320 {
2321 // polygon
2322 const QgsPolygon *polygon = static_cast< const QgsPolygon * >( geom );
2323 std::unique_ptr< QgsPolygon > result = std::make_unique< QgsPolygon >();
2324
2325 result->setExteriorRing( squareWavesRandomizedAlongLine( static_cast< const QgsLineString * >( polygon->exteriorRing() ),
2326 minimumWavelength, maximumWavelength, minimumAmplitude, maximumAmplitude, uniformDist, mt ).release() );
2327 for ( int i = 0; i < polygon->numInteriorRings(); ++i )
2328 {
2329 result->addInteriorRing( squareWavesRandomizedAlongLine( static_cast< const QgsLineString * >( polygon->interiorRing( i ) ),
2330 minimumWavelength, maximumWavelength, minimumAmplitude, maximumAmplitude, uniformDist, mt ).release() );
2331 }
2332
2333 return result;
2334 }
2335}
2336
2337QgsGeometry QgsInternalGeometryEngine::squareWaves( double wavelength, double amplitude, bool strictWavelength ) const
2338{
2339 if ( wavelength < 0 || qgsDoubleNear( wavelength, 0 ) )
2340 return QgsGeometry();
2341
2342 mLastError.clear();
2343 if ( !mGeometry )
2344 {
2345 return QgsGeometry();
2346 }
2347
2349 {
2350 return QgsGeometry( mGeometry->clone() ); // point geometry, nothing to do
2351 }
2352
2353 if ( const QgsGeometryCollection *gc = qgsgeometry_cast< const QgsGeometryCollection *>( mGeometry ) )
2354 {
2355 int numGeom = gc->numGeometries();
2356 QVector< QgsAbstractGeometry * > geometryList;
2357 geometryList.reserve( numGeom );
2358 for ( int i = 0; i < numGeom; ++i )
2359 {
2360 geometryList << squareWavesPrivate( gc->geometryN( i ), wavelength, amplitude, strictWavelength ).release();
2361 }
2362
2363 QgsGeometry first = QgsGeometry( geometryList.takeAt( 0 ) );
2364 for ( QgsAbstractGeometry *g : std::as_const( geometryList ) )
2365 {
2367 }
2368 return first;
2369 }
2370 else
2371 {
2372 return QgsGeometry( squareWavesPrivate( mGeometry, wavelength, amplitude, strictWavelength ) );
2373 }
2374}
2375
2376QgsGeometry QgsInternalGeometryEngine::squareWavesRandomized( double minimumWavelength, double maximumWavelength, double minimumAmplitude, double maximumAmplitude, unsigned long seed ) const
2377{
2378 if ( minimumWavelength < 0 || qgsDoubleNear( minimumWavelength, 0 ) || maximumWavelength < 0 || qgsDoubleNear( maximumWavelength, 0 ) || maximumWavelength < minimumWavelength )
2379 return QgsGeometry();
2380
2381 mLastError.clear();
2382 if ( !mGeometry )
2383 {
2384 return QgsGeometry();
2385 }
2386
2387 std::random_device rd;
2388 std::mt19937 mt( seed == 0 ? rd() : seed );
2389 std::uniform_real_distribution<> uniformDist( 0, 1 );
2390
2392 {
2393 return QgsGeometry( mGeometry->clone() ); // point geometry, nothing to do
2394 }
2395
2396 if ( const QgsGeometryCollection *gc = qgsgeometry_cast< const QgsGeometryCollection *>( mGeometry ) )
2397 {
2398 int numGeom = gc->numGeometries();
2399 QVector< QgsAbstractGeometry * > geometryList;
2400 geometryList.reserve( numGeom );
2401 for ( int i = 0; i < numGeom; ++i )
2402 {
2403 geometryList << squareWavesRandomizedPrivate( gc->geometryN( i ), minimumWavelength, maximumWavelength, minimumAmplitude, maximumAmplitude, uniformDist, mt ).release();
2404 }
2405
2406 QgsGeometry first = QgsGeometry( geometryList.takeAt( 0 ) );
2407 for ( QgsAbstractGeometry *g : std::as_const( geometryList ) )
2408 {
2410 }
2411 return first;
2412 }
2413 else
2414 {
2415 return QgsGeometry( squareWavesRandomizedPrivate( mGeometry, minimumWavelength, maximumWavelength, minimumAmplitude, maximumAmplitude, uniformDist, mt ) );
2416 }
2417}
2418
2419std::unique_ptr< QgsLineString > roundWavesAlongLine( const QgsLineString *line, double wavelength, const double amplitude, const bool strictWavelength )
2420{
2421 const int totalPoints = line->numPoints();
2422 if ( totalPoints < 2 )
2423 return nullptr;
2424
2425 const double *x = line->xData();
2426 const double *y = line->yData();
2427
2428 double prevX = *x++;
2429 double prevY = *y++;
2430
2431 const double totalLength = line->length();
2432
2433 const int maxRepetitions = std::ceil( totalLength / wavelength );
2434 if ( !strictWavelength )
2435 wavelength = totalLength / maxRepetitions;
2436
2437 const int segments = 10;
2438
2439 int side = -1;
2440
2441 double xOutBuffer[4] { prevX, prevX, prevX, prevX };
2442 double yOutBuffer[4] { prevY, prevY, prevY, prevY };
2443 bool isFirstPart = true;
2444
2445 double distanceToNextPointFromStartOfSegment = wavelength / 8;
2446 int bufferIndex = 1;
2447 std::unique_ptr< QgsLineString > out = std::make_unique< QgsLineString >();
2448
2450 double remainingDistance = totalLength;
2451 double totalCoveredDistance = 0;
2452
2453 for ( int i = 1; i < totalPoints; ++i )
2454 {
2455 double thisX = *x++;
2456 double thisY = *y++;
2457
2458 segmentAngleRadians = QgsGeometryUtilsBase::lineAngle( prevX, prevY, thisX, thisY );
2459 const double segmentLength = QgsGeometryUtilsBase::distance2D( thisX, thisY, prevX, prevY );
2460 while ( distanceToNextPointFromStartOfSegment < segmentLength || qgsDoubleNear( distanceToNextPointFromStartOfSegment, segmentLength ) )
2461 {
2462 // point falls on this segment - truncate to segment length if qgsDoubleNear test was actually > segment length
2463 const double distanceToPoint = std::min( distanceToNextPointFromStartOfSegment, segmentLength );
2464 double pX, pY;
2465 QgsGeometryUtilsBase::pointOnLineWithDistance( prevX, prevY, thisX, thisY, distanceToPoint, pX, pY );
2466 remainingDistance = totalLength - totalCoveredDistance - distanceToPoint;
2467
2468 const double sinAngle = std::sin( segmentAngleRadians + M_PI_2 );
2469 const double cosAngle = std::cos( segmentAngleRadians + M_PI_2 );
2470
2471 if ( bufferIndex == 0 )
2472 {
2473 xOutBuffer[0] = pX + side * amplitude * sinAngle;
2474 yOutBuffer[0] = pY + side * amplitude * cosAngle;
2475 bufferIndex = 1;
2476 distanceToNextPointFromStartOfSegment += wavelength / 4;
2477 }
2478 else if ( bufferIndex == 1 && isFirstPart )
2479 {
2480 xOutBuffer[1] = ( xOutBuffer[0] + pX - side * amplitude * sinAngle ) * 0.5;
2481 yOutBuffer[1] = ( yOutBuffer[0] + pY - side * amplitude * cosAngle ) * 0.5;
2482 xOutBuffer[2] = pX - side * amplitude * sinAngle;
2483 yOutBuffer[2] = pY - side * amplitude * cosAngle;
2484 bufferIndex = 2;
2485 distanceToNextPointFromStartOfSegment += wavelength / 8;
2486 }
2487 else if ( bufferIndex == 1 )
2488 {
2489 xOutBuffer[1] = pX + side * amplitude * sinAngle;
2490 yOutBuffer[1] = pY + side * amplitude * cosAngle;
2491 xOutBuffer[2] = pX - side * amplitude * sinAngle;
2492 yOutBuffer[2] = pY - side * amplitude * cosAngle;
2493 bufferIndex = 2;
2494 distanceToNextPointFromStartOfSegment += wavelength / 4;
2495 }
2496 else if ( bufferIndex == 2 )
2497 {
2498 xOutBuffer[3] = pX - side * amplitude * sinAngle;
2499 yOutBuffer[3] = pY - side * amplitude * cosAngle;
2500
2501 if ( isFirstPart )
2502 {
2503 xOutBuffer[2] = ( xOutBuffer[2] + xOutBuffer[3] ) * 0.5;
2504 yOutBuffer[2] = ( yOutBuffer[2] + yOutBuffer[3] ) * 0.5;
2505 isFirstPart = false;
2506 }
2507
2508 // flush curve
2509 std::unique_ptr< QgsLineString > bezier( QgsLineString::fromBezierCurve( QgsPoint( xOutBuffer[0], yOutBuffer[0] ),
2510 QgsPoint( xOutBuffer[1], yOutBuffer[1] ),
2511 QgsPoint( xOutBuffer[2], yOutBuffer[2] ),
2512 QgsPoint( xOutBuffer[3], yOutBuffer[3] ),
2513 segments ) );
2514 out->append( bezier.get() );
2515
2516 // shuffle buffer alone
2517 xOutBuffer[0] = xOutBuffer[3];
2518 yOutBuffer[0] = yOutBuffer[3];
2519 bufferIndex = 1;
2520 side = -side;
2521 distanceToNextPointFromStartOfSegment += wavelength / 4;
2522 }
2523 }
2524 totalCoveredDistance += segmentLength;
2525 prevX = thisX;
2526 prevY = thisY;
2527 distanceToNextPointFromStartOfSegment -= segmentLength;
2528 }
2529
2530 const double midX = prevX - remainingDistance / 2 * std::sin( segmentAngleRadians );
2531 const double midY = prevY - remainingDistance / 2 * std::cos( segmentAngleRadians );
2532 const double pX = midX + side * amplitude * std::sin( segmentAngleRadians + M_PI_2 );
2533 const double pY = midY + side * amplitude * std::cos( segmentAngleRadians + M_PI_2 );
2534
2535 std::unique_ptr< QgsLineString > bezier( QgsLineString::fromBezierCurve( out->endPoint(),
2536 QgsPoint( ( out->endPoint().x() + pX ) * 0.5, ( out->endPoint().y() + pY ) * 0.5 ),
2537 QgsPoint( ( prevX + pX ) * 0.5, ( prevY + pY ) * 0.5 ),
2538 QgsPoint( prevX, prevY ),
2539 segments / 2 ) );
2540 out->append( bezier.get() );
2541
2542 return out;
2543}
2544
2545std::unique_ptr< QgsLineString > roundWavesRandomizedAlongLine( const QgsLineString *line,
2546 const double minimumWavelength, const double maximumWavelength,
2547 const double minimumAmplitude, const double maximumAmplitude,
2548 std::uniform_real_distribution<> &uniformDist, std::mt19937 &mt )
2549{
2550 const int totalPoints = line->numPoints();
2551 if ( totalPoints < 2 )
2552 return nullptr;
2553
2554 const double *x = line->xData();
2555 const double *y = line->yData();
2556
2557 double prevX = *x++;
2558 double prevY = *y++;
2559
2560 const double totalLength = line->length();
2561
2562 const int segments = 10;
2563
2564 int side = -1;
2565
2566 double xOutBuffer[4] { prevX, prevX, prevX, prevX };
2567 double yOutBuffer[4] { prevY, prevY, prevY, prevY };
2568 bool isFirstPart = true;
2569
2570 double amplitude = uniformDist( mt ) * ( maximumAmplitude - minimumAmplitude ) + minimumAmplitude;
2571 double wavelength = uniformDist( mt ) * ( maximumWavelength - minimumWavelength ) + minimumWavelength;
2572
2573 double distanceToNextPointFromStartOfSegment = wavelength / 8;
2574 int bufferIndex = 1;
2575 std::unique_ptr< QgsLineString > out = std::make_unique< QgsLineString >();
2576
2578
2579 double remainingDistance = totalLength;
2580 double totalCoveredDistance = 0;
2581
2582 for ( int i = 1; i < totalPoints; ++i )
2583 {
2584 double thisX = *x++;
2585 double thisY = *y++;
2586
2587 segmentAngleRadians = QgsGeometryUtilsBase::lineAngle( prevX, prevY, thisX, thisY );
2588 const double segmentLength = QgsGeometryUtilsBase::distance2D( thisX, thisY, prevX, prevY );
2589 while ( distanceToNextPointFromStartOfSegment < segmentLength || qgsDoubleNear( distanceToNextPointFromStartOfSegment, segmentLength ) )
2590 {
2591 // point falls on this segment - truncate to segment length if qgsDoubleNear test was actually > segment length
2592 const double distanceToPoint = std::min( distanceToNextPointFromStartOfSegment, segmentLength );
2593 double pX, pY;
2594 QgsGeometryUtilsBase::pointOnLineWithDistance( prevX, prevY, thisX, thisY, distanceToPoint, pX, pY );
2595 remainingDistance = totalLength - totalCoveredDistance - distanceToPoint;
2596
2597 const double sinAngle = std::sin( segmentAngleRadians + M_PI_2 );
2598 const double cosAngle = std::cos( segmentAngleRadians + M_PI_2 );
2599
2600 if ( bufferIndex == 0 )
2601 {
2602 xOutBuffer[0] = pX + side * amplitude * sinAngle;
2603 yOutBuffer[0] = pY + side * amplitude * cosAngle;
2604 bufferIndex = 1;
2605 distanceToNextPointFromStartOfSegment += wavelength / 4;
2606 }
2607 else if ( bufferIndex == 1 && isFirstPart )
2608 {
2609 xOutBuffer[1] = ( xOutBuffer[0] + pX - side * amplitude * sinAngle ) * 0.5;
2610 yOutBuffer[1] = ( yOutBuffer[0] + pY - side * amplitude * cosAngle ) * 0.5;
2611 xOutBuffer[2] = pX - side * amplitude * sinAngle;
2612 yOutBuffer[2] = pY - side * amplitude * cosAngle;
2613 bufferIndex = 2;
2614 distanceToNextPointFromStartOfSegment += wavelength / 8;
2615 }
2616 else if ( bufferIndex == 1 )
2617 {
2618 xOutBuffer[1] = pX + side * amplitude * sinAngle;
2619 yOutBuffer[1] = pY + side * amplitude * cosAngle;
2620 amplitude = uniformDist( mt ) * ( maximumAmplitude - minimumAmplitude ) + minimumAmplitude;
2621 xOutBuffer[2] = pX - side * amplitude * sinAngle;
2622 bufferIndex = 2;
2623 yOutBuffer[2] = pY - side * amplitude * cosAngle;
2624 distanceToNextPointFromStartOfSegment += wavelength / 4;
2625 }
2626 else if ( bufferIndex == 2 )
2627 {
2628 xOutBuffer[3] = pX - side * amplitude * sinAngle;
2629 yOutBuffer[3] = pY - side * amplitude * cosAngle;
2630
2631 if ( isFirstPart )
2632 {
2633 xOutBuffer[2] = ( xOutBuffer[2] + xOutBuffer[3] ) * 0.5;
2634 yOutBuffer[2] = ( yOutBuffer[2] + yOutBuffer[3] ) * 0.5;
2635 isFirstPart = false;
2636 }
2637
2638 // flush curve
2639 std::unique_ptr< QgsLineString > bezier( QgsLineString::fromBezierCurve( QgsPoint( xOutBuffer[0], yOutBuffer[0] ),
2640 QgsPoint( xOutBuffer[1], yOutBuffer[1] ),
2641 QgsPoint( xOutBuffer[2], yOutBuffer[2] ),
2642 QgsPoint( xOutBuffer[3], yOutBuffer[3] ),
2643 segments ) );
2644 out->append( bezier.get() );
2645
2646 // shuffle buffer alone
2647 xOutBuffer[0] = xOutBuffer[3];
2648 yOutBuffer[0] = yOutBuffer[3];
2649 bufferIndex = 1;
2650 side = -side;
2651
2652 wavelength = uniformDist( mt ) * ( maximumWavelength - minimumWavelength ) + minimumWavelength;
2653
2654 distanceToNextPointFromStartOfSegment += wavelength / 4;
2655 }
2656 }
2657 totalCoveredDistance += segmentLength;
2658
2659 prevX = thisX;
2660 prevY = thisY;
2661 distanceToNextPointFromStartOfSegment -= segmentLength;
2662 }
2663
2664 const double midX = prevX - remainingDistance / 2 * std::sin( segmentAngleRadians );
2665 const double midY = prevY - remainingDistance / 2 * std::cos( segmentAngleRadians );
2666 const double pX = midX + side * amplitude * std::sin( segmentAngleRadians + M_PI_2 );
2667 const double pY = midY + side * amplitude * std::cos( segmentAngleRadians + M_PI_2 );
2668
2669 if ( out->isEmpty() )
2670 {
2671 std::unique_ptr< QgsLineString > bezier( QgsLineString::fromBezierCurve( line->startPoint(),
2672 QgsPoint( ( line->startPoint().x() + pX ) * 0.5, ( line->startPoint().y() + pY ) * 0.5 ),
2673 QgsPoint( ( prevX + pX ) * 0.5, ( prevY + pY ) * 0.5 ),
2674 QgsPoint( prevX, prevY ),
2675 segments ) );
2676 out->append( bezier.get() );
2677 }
2678 else
2679 {
2680 std::unique_ptr< QgsLineString > bezier( QgsLineString::fromBezierCurve( out->endPoint(),
2681 QgsPoint( ( out->endPoint().x() + pX ) * 0.5, ( out->endPoint().y() + pY ) * 0.5 ),
2682 QgsPoint( ( prevX + pX ) * 0.5, ( prevY + pY ) * 0.5 ),
2683 QgsPoint( prevX, prevY ),
2684 segments ) );
2685 out->append( bezier.get() );
2686 }
2687
2688 return out;
2689}
2690
2691std::unique_ptr< QgsAbstractGeometry > roundWavesPrivate( const QgsAbstractGeometry *geom, double wavelength, double amplitude, bool strictWavelength )
2692{
2693 std::unique_ptr< QgsAbstractGeometry > segmentizedCopy;
2694 if ( QgsWkbTypes::isCurvedType( geom->wkbType() ) )
2695 {
2696 segmentizedCopy.reset( geom->segmentize() );
2697 geom = segmentizedCopy.get();
2698 }
2699
2701 {
2702 return roundWavesAlongLine( static_cast< const QgsLineString * >( geom ), wavelength, amplitude, strictWavelength );
2703 }
2704 else
2705 {
2706 // polygon
2707 const QgsPolygon *polygon = static_cast< const QgsPolygon * >( geom );
2708 std::unique_ptr< QgsPolygon > result = std::make_unique< QgsPolygon >();
2709
2710 result->setExteriorRing( roundWavesAlongLine( static_cast< const QgsLineString * >( polygon->exteriorRing() ),
2711 wavelength, amplitude, strictWavelength ).release() );
2712 for ( int i = 0; i < polygon->numInteriorRings(); ++i )
2713 {
2714 result->addInteriorRing( roundWavesAlongLine( static_cast< const QgsLineString * >( polygon->interiorRing( i ) ),
2715 wavelength, amplitude, strictWavelength ).release() );
2716 }
2717
2718 return result;
2719 }
2720}
2721
2722std::unique_ptr< QgsAbstractGeometry > roundWavesRandomizedPrivate( const QgsAbstractGeometry *geom, double minimumWavelength, double maximumWavelength, double minimumAmplitude, double maximumAmplitude, std::uniform_real_distribution<> &uniformDist, std::mt19937 &mt )
2723{
2724 std::unique_ptr< QgsAbstractGeometry > segmentizedCopy;
2725 if ( QgsWkbTypes::isCurvedType( geom->wkbType() ) )
2726 {
2727 segmentizedCopy.reset( geom->segmentize() );
2728 geom = segmentizedCopy.get();
2729 }
2730
2732 {
2733 return roundWavesRandomizedAlongLine( static_cast< const QgsLineString * >( geom ), minimumWavelength, maximumWavelength, minimumAmplitude, maximumAmplitude, uniformDist, mt );
2734 }
2735 else
2736 {
2737 // polygon
2738 const QgsPolygon *polygon = static_cast< const QgsPolygon * >( geom );
2739 std::unique_ptr< QgsPolygon > result = std::make_unique< QgsPolygon >();
2740
2741 result->setExteriorRing( roundWavesRandomizedAlongLine( static_cast< const QgsLineString * >( polygon->exteriorRing() ),
2742 minimumWavelength, maximumWavelength, minimumAmplitude, maximumAmplitude, uniformDist, mt ).release() );
2743 for ( int i = 0; i < polygon->numInteriorRings(); ++i )
2744 {
2745 result->addInteriorRing( roundWavesRandomizedAlongLine( static_cast< const QgsLineString * >( polygon->interiorRing( i ) ),
2746 minimumWavelength, maximumWavelength, minimumAmplitude, maximumAmplitude, uniformDist, mt ).release() );
2747 }
2748
2749 return result;
2750 }
2751}
2752
2753QgsGeometry QgsInternalGeometryEngine::roundWaves( double wavelength, double amplitude, bool strictWavelength ) const
2754{
2755 if ( wavelength < 0 || qgsDoubleNear( wavelength, 0 ) )
2756 return QgsGeometry();
2757
2758 mLastError.clear();
2759 if ( !mGeometry )
2760 {
2761 return QgsGeometry();
2762 }
2763
2765 {
2766 return QgsGeometry( mGeometry->clone() ); // point geometry, nothing to do
2767 }
2768
2769 if ( const QgsGeometryCollection *gc = qgsgeometry_cast< const QgsGeometryCollection *>( mGeometry ) )
2770 {
2771 int numGeom = gc->numGeometries();
2772 QVector< QgsAbstractGeometry * > geometryList;
2773 geometryList.reserve( numGeom );
2774 for ( int i = 0; i < numGeom; ++i )
2775 {
2776 geometryList << roundWavesPrivate( gc->geometryN( i ), wavelength, amplitude, strictWavelength ).release();
2777 }
2778
2779 QgsGeometry first = QgsGeometry( geometryList.takeAt( 0 ) );
2780 for ( QgsAbstractGeometry *g : std::as_const( geometryList ) )
2781 {
2783 }
2784 return first;
2785 }
2786 else
2787 {
2788 return QgsGeometry( roundWavesPrivate( mGeometry, wavelength, amplitude, strictWavelength ) );
2789 }
2790}
2791
2792QgsGeometry QgsInternalGeometryEngine::roundWavesRandomized( double minimumWavelength, double maximumWavelength, double minimumAmplitude, double maximumAmplitude, unsigned long seed ) const
2793{
2794 if ( minimumWavelength < 0 || qgsDoubleNear( minimumWavelength, 0 ) || maximumWavelength < 0 || qgsDoubleNear( maximumWavelength, 0 ) || maximumWavelength < minimumWavelength )
2795 return QgsGeometry();
2796
2797 mLastError.clear();
2798 if ( !mGeometry )
2799 {
2800 return QgsGeometry();
2801 }
2802
2803 std::random_device rd;
2804 std::mt19937 mt( seed == 0 ? rd() : seed );
2805 std::uniform_real_distribution<> uniformDist( 0, 1 );
2806
2808 {
2809 return QgsGeometry( mGeometry->clone() ); // point geometry, nothing to do
2810 }
2811
2812 if ( const QgsGeometryCollection *gc = qgsgeometry_cast< const QgsGeometryCollection *>( mGeometry ) )
2813 {
2814 int numGeom = gc->numGeometries();
2815 QVector< QgsAbstractGeometry * > geometryList;
2816 geometryList.reserve( numGeom );
2817 for ( int i = 0; i < numGeom; ++i )
2818 {
2819 geometryList << roundWavesRandomizedPrivate( gc->geometryN( i ), minimumWavelength, maximumWavelength, minimumAmplitude, maximumAmplitude, uniformDist, mt ).release();
2820 }
2821
2822 QgsGeometry first = QgsGeometry( geometryList.takeAt( 0 ) );
2823 for ( QgsAbstractGeometry *g : std::as_const( geometryList ) )
2824 {
2826 }
2827 return first;
2828 }
2829 else
2830 {
2831 return QgsGeometry( roundWavesRandomizedPrivate( mGeometry, minimumWavelength, maximumWavelength, minimumAmplitude, maximumAmplitude, uniformDist, mt ) );
2832 }
2833}
2834
2835std::unique_ptr< QgsMultiLineString > dashPatternAlongLine( const QgsLineString *line,
2836 const QVector< double> &pattern,
2840 double patternOffset )
2841{
2842 const int totalPoints = line->numPoints();
2843 if ( totalPoints < 2 )
2844 return nullptr;
2845
2846 const int patternSize = pattern.size();
2847
2848 const double *x = line->xData();
2849 const double *y = line->yData();
2850
2851 double prevX = *x++;
2852 double prevY = *y++;
2853
2854 std::unique_ptr< QgsMultiLineString > result = std::make_unique< QgsMultiLineString >();
2855
2856 QVector< double > outX;
2857 QVector< double > outY;
2858 const double totalLength = line->length();
2859
2860 double patternLength = 0;
2861 double patternDashLength = 0;
2862 double patternGapLength = 0;
2863 for ( int i = 0; i < pattern.size(); ++i )
2864 {
2865 patternLength += pattern.at( i );
2866 if ( i % 2 == 0 )
2867 patternDashLength += pattern.at( i );
2868 else
2869 patternGapLength += pattern.at( i );
2870 }
2871
2872 double firstPatternLength = 0;
2873 double firstPatternDashLength = 0;
2874 double firstPatternGapLength = 0;
2875 switch ( startRule )
2876 {
2879 firstPatternLength = patternLength;
2880 firstPatternDashLength = patternDashLength;
2881 firstPatternGapLength = patternGapLength;
2882 break;
2884 firstPatternLength = patternLength - pattern.at( 0 ) * 0.5; // remove half of first dash
2885 firstPatternDashLength = patternDashLength - pattern.at( 0 ) * 0.5;
2886 firstPatternGapLength = patternGapLength;
2887 break;
2889 firstPatternLength = pattern.at( patternSize - 1 ); // start at LAST dash
2890 firstPatternDashLength = 0;
2891 firstPatternGapLength = pattern.at( patternSize - 1 );
2892 break;
2894 firstPatternLength = pattern.at( patternSize - 1 ) * 0.5; // start at half of last dash
2895 firstPatternDashLength = 0;
2896 firstPatternGapLength = pattern.at( patternSize - 1 ) * 0.5;
2897 break;
2898 }
2899
2900 const bool isSmallEnoughForSinglePattern = ( totalLength - firstPatternLength ) < patternLength * 0.5;
2901
2902 double lastPatternLength = isSmallEnoughForSinglePattern ? firstPatternLength : patternLength;
2903 double lastPatternDashLength = isSmallEnoughForSinglePattern ? firstPatternDashLength : patternDashLength;
2904 double lastPatternGapLength = isSmallEnoughForSinglePattern ? firstPatternGapLength : patternGapLength;
2905 switch ( endRule )
2906 {
2908 lastPatternLength = 0;
2909 lastPatternDashLength = 0;
2910 lastPatternGapLength = 0;
2911 break;
2913 lastPatternLength -= pattern.at( patternSize - 1 ); // remove last gap
2914 lastPatternGapLength -= pattern.at( patternSize - 1 );
2915 break;
2917 lastPatternLength -= pattern.at( patternSize - 1 ) + pattern.at( patternSize - 2 ) * 0.5; // remove last gap, half of last dash
2918 lastPatternDashLength -= pattern.at( patternSize - 2 ) * 0.5;
2919 lastPatternGapLength -= pattern.at( patternSize - 1 );
2920 break;
2922 lastPatternGapLength = patternGapLength;
2923 break;
2925 lastPatternLength -= pattern.at( patternSize - 1 ) * 0.5; // remove half of last gap
2926 lastPatternGapLength -= pattern.at( patternSize - 1 ) * 0.5;
2927 break;
2928 }
2929
2930 const double remainingLengthForCompletePatterns = totalLength - ( !isSmallEnoughForSinglePattern ? firstPatternLength : 0 ) - lastPatternLength;
2931 const int middlePatternRepetitions = std::max( static_cast< int >( std::round( remainingLengthForCompletePatterns / patternLength ) ), 0 );
2932
2933 const double totalUnscaledLengthOfPatterns = ( !isSmallEnoughForSinglePattern ? firstPatternLength : 0 ) + middlePatternRepetitions * patternLength + lastPatternLength;
2934 const double totalUnscaledLengthOfDashes = ( !isSmallEnoughForSinglePattern ? firstPatternDashLength : 0 ) + middlePatternRepetitions * patternDashLength + lastPatternDashLength;
2935 const double totalUnscaledLengthOfGaps = ( !isSmallEnoughForSinglePattern ? firstPatternGapLength : 0 ) + middlePatternRepetitions * patternGapLength + lastPatternGapLength;
2936
2937 double dashLengthScalingFactor = 1;
2938 double gapLengthScalingFactor = 1;
2940 {
2941 // calculate scaling factors
2942 const double lengthToShrinkBy = totalUnscaledLengthOfPatterns - totalLength;
2943
2945 {
2947 dashLengthScalingFactor = totalLength / totalUnscaledLengthOfPatterns;
2948 gapLengthScalingFactor = dashLengthScalingFactor;
2949 break;
2951 dashLengthScalingFactor = ( totalUnscaledLengthOfDashes - lengthToShrinkBy ) / totalUnscaledLengthOfDashes;
2952 break;
2954 gapLengthScalingFactor = ( totalUnscaledLengthOfGaps - lengthToShrinkBy ) / totalUnscaledLengthOfGaps;
2955 break;
2956 }
2957 }
2958
2959 dashLengthScalingFactor = std::max( dashLengthScalingFactor, 0.0 );
2960 gapLengthScalingFactor = std::max( gapLengthScalingFactor, 0.0 );
2961
2962 const int maxPatterns = middlePatternRepetitions + 2;
2963 result->reserve( maxPatterns );
2964
2965 int patternIndex = 0;
2966 double distanceToNextPointFromStartOfSegment = pattern.at( 0 ) * dashLengthScalingFactor;
2967 bool isDash = true;
2968 switch ( startRule )
2969 {
2972 break;
2974 distanceToNextPointFromStartOfSegment *= 0.5; // skip to half way through first dash
2975 break;
2977 patternIndex = patternSize - 1; // start at last gap
2978 isDash = false;
2979 distanceToNextPointFromStartOfSegment = pattern.at( patternSize - 1 ) * gapLengthScalingFactor;
2980 break;
2982 patternIndex = patternSize - 1; // skip straight to half way through last gap
2983 isDash = false;
2984 distanceToNextPointFromStartOfSegment = 0.5 * pattern.at( patternSize - 1 ) * gapLengthScalingFactor;
2985 break;
2986 }
2987
2988 const double adjustedOffset = fmod( patternOffset, patternLength );
2989 const double scaledOffset = ( adjustedOffset < 0 ? ( adjustedOffset + patternLength ) : adjustedOffset ) * ( gapLengthScalingFactor + dashLengthScalingFactor ) / 2;
2990 if ( !qgsDoubleNear( scaledOffset, 0 ) )
2991 {
2992 // shuffle pattern along by offset
2993 double remainingOffset = scaledOffset;
2994 while ( remainingOffset > 0 )
2995 {
2996 if ( distanceToNextPointFromStartOfSegment > remainingOffset )
2997 {
2998 distanceToNextPointFromStartOfSegment -= remainingOffset;
2999 break;
3000 }
3001
3002 remainingOffset -= distanceToNextPointFromStartOfSegment;
3003 isDash = !isDash;
3004 patternIndex++;
3005 if ( patternIndex == patternSize )
3006 patternIndex = 0;
3007
3008 distanceToNextPointFromStartOfSegment = pattern.at( patternIndex ) * ( isDash ? dashLengthScalingFactor : gapLengthScalingFactor );
3009 }
3010 }
3011
3012 if ( isDash )
3013 {
3014 outX.append( prevX );
3015 outY.append( prevY );
3016 }
3017
3018 for ( int i = 1; i < totalPoints; ++i )
3019 {
3020 double thisX = *x++;
3021 double thisY = *y++;
3022
3023 const double segmentLength = QgsGeometryUtilsBase::distance2D( thisX, thisY, prevX, prevY );
3024 while ( distanceToNextPointFromStartOfSegment < segmentLength || qgsDoubleNear( distanceToNextPointFromStartOfSegment, segmentLength ) )
3025 {
3026 // point falls on this segment - truncate to segment length if qgsDoubleNear test was actually > segment length
3027 const double distanceToPoint = std::min( distanceToNextPointFromStartOfSegment, segmentLength );
3028 double pX, pY;
3029 QgsGeometryUtilsBase::pointOnLineWithDistance( prevX, prevY, thisX, thisY, distanceToPoint, pX, pY );
3030
3031 outX.append( pX );
3032 outY.append( pY );
3033 if ( isDash )
3034 {
3035 result->addGeometry( new QgsLineString( outX, outY ) );
3036 outX.resize( 0 );
3037 outY.resize( 0 );
3038 }
3039
3040 isDash = !isDash;
3041 patternIndex++;
3042 if ( patternIndex >= patternSize )
3043 patternIndex = 0;
3044
3045 distanceToNextPointFromStartOfSegment += pattern.at( patternIndex ) * ( isDash ? dashLengthScalingFactor : gapLengthScalingFactor );
3046 }
3047
3048 if ( isDash )
3049 {
3050 outX.append( thisX );
3051 outY.append( thisY );
3052 }
3053
3054 prevX = thisX;
3055 prevY = thisY;
3056 distanceToNextPointFromStartOfSegment -= segmentLength;
3057 }
3058
3059 if ( isDash )
3060 {
3061 outX.append( prevX );
3062 outY.append( prevY );
3063 result->addGeometry( new QgsLineString( outX, outY ) );
3064 }
3065
3066 return result;
3067}
3068
3069std::unique_ptr< QgsAbstractGeometry > applyDashPatternPrivate( const QgsAbstractGeometry *geom,
3070 const QVector<double> &pattern,
3074 double patternOffset )
3075{
3076 std::unique_ptr< QgsAbstractGeometry > segmentizedCopy;
3077 if ( QgsWkbTypes::isCurvedType( geom->wkbType() ) )
3078 {
3079 segmentizedCopy.reset( geom->segmentize() );
3080 geom = segmentizedCopy.get();
3081 }
3082
3084 {
3085 return dashPatternAlongLine( static_cast< const QgsLineString * >( geom ), pattern, startRule, endRule, adjustment, patternOffset );
3086 }
3087 else
3088 {
3089 // polygon
3090 const QgsPolygon *polygon = static_cast< const QgsPolygon * >( geom );
3091 std::unique_ptr< QgsMultiLineString > result = std::make_unique< QgsMultiLineString >();
3092
3093 std::unique_ptr< QgsMultiLineString > exteriorParts = dashPatternAlongLine( static_cast< const QgsLineString * >( polygon->exteriorRing() ), pattern, startRule, endRule, adjustment, patternOffset );
3094 for ( int i = 0; i < exteriorParts->numGeometries(); ++i )
3095 result->addGeometry( exteriorParts->geometryN( i )->clone() );
3096
3097 for ( int i = 0; i < polygon->numInteriorRings(); ++i )
3098 {
3099 std::unique_ptr< QgsMultiLineString > ringParts = dashPatternAlongLine( static_cast< const QgsLineString * >( polygon->interiorRing( i ) ), pattern, startRule, endRule, adjustment, patternOffset );
3100 for ( int j = 0; j < ringParts->numGeometries(); ++j )
3101 result->addGeometry( ringParts->geometryN( j )->clone() );
3102 }
3103
3104 return result;
3105 }
3106}
3107
3109{
3110 if ( pattern.size() < 2 )
3111 return QgsGeometry( mGeometry->clone() );
3112
3113 mLastError.clear();
3114 if ( !mGeometry || mGeometry->isEmpty() )
3115 {
3116 return QgsGeometry();
3117 }
3118
3120 {
3121 return QgsGeometry( mGeometry->clone() ); // point geometry, nothing to do
3122 }
3123
3124 if ( const QgsGeometryCollection *gc = qgsgeometry_cast< const QgsGeometryCollection *>( mGeometry ) )
3125 {
3126 int numGeom = gc->numGeometries();
3127 QVector< QgsAbstractGeometry * > geometryList;
3128 geometryList.reserve( numGeom );
3129 for ( int i = 0; i < numGeom; ++i )
3130 {
3131 geometryList << applyDashPatternPrivate( gc->geometryN( i ), pattern, startRule, endRule, adjustment, patternOffset ).release();
3132 }
3133
3134 QgsGeometry first = QgsGeometry( geometryList.takeAt( 0 ) );
3135 for ( QgsAbstractGeometry *g : std::as_const( geometryList ) )
3136 {
3137 if ( const QgsGeometryCollection *collection = qgsgeometry_cast< QgsGeometryCollection * >( g ) )
3138 {
3139 for ( int j = 0; j < collection->numGeometries(); ++j )
3140 {
3141 first.addPart( collection->geometryN( j )->clone() );
3142 }
3143 delete collection;
3144 }
3145 else
3146 {
3148 }
3149 }
3150 return first;
3151 }
3152 else
3153 {
3154 return QgsGeometry( applyDashPatternPrivate( mGeometry, pattern, startRule, endRule, adjustment, patternOffset ) );
3155 }
3156}
Definition: qgis.h:2675
@ ScaleDashOnly
@ ScaleBothDashAndGap
Both the dash and gap lengths are adjusted equally.
@ ScaleGapOnly
@ Polygon
Polygons.
DashPatternLineEndingRule
Dash pattern line ending rules.
Definition: qgis.h:2660
@ HalfDash
Start or finish the pattern with a half length dash.
@ HalfGap
Start or finish the pattern with a half length gap.
@ FullGap
Start or finish the pattern with a full gap.
@ FullDash
Start or finish the pattern with a full dash.
WkbType
The WKB type describes the number of dimensions a geometry has.
Definition: qgis.h:182
@ CompoundCurve
CompoundCurve.
@ LineString
LineString.
@ Polygon
Polygon.
@ CircularString
CircularString.
The vertex_iterator class provides STL-style iterator for vertices.
Abstract base class for all geometries.
vertex_iterator vertices_end() const
Returns STL-style iterator pointing to the imaginary vertex after the last vertex of the geometry.
virtual QgsAbstractGeometry * segmentize(double tolerance=M_PI/180., SegmentationToleranceType toleranceType=MaximumAngle) const
Returns a version of the geometry without curves.
bool isMeasure() const
Returns true if the geometry contains m values.
virtual QgsRectangle boundingBox() const
Returns the minimal bounding box for the geometry.
bool is3D() const
Returns true if the geometry is 3D and contains a z-value.
virtual int nCoordinates() const
Returns the number of nodes contained in the geometry.
Qgis::WkbType wkbType() const
Returns the WKB type of the geometry.
const_part_iterator const_parts_end() const
Returns STL-style iterator pointing to the imaginary const part after the last part of the geometry.
virtual bool isEmpty() const
Returns true if the geometry is empty.
vertex_iterator vertices_begin() const
Returns STL-style iterator pointing to the first vertex of the geometry.
const_part_iterator const_parts_begin() const
Returns STL-style iterator pointing to the const first part of the geometry.
virtual QgsAbstractGeometry * clone() const =0
Clones the geometry by performing a deep copy.
Circle geometry type.
Definition: qgscircle.h:43
Compound curve geometry type.
int nCurves() const
Returns the number of curves in the geometry.
const QgsCurve * curveAt(int i) const
Returns the curve at the specified index.
Curve polygon geometry type.
int numInteriorRings() const
Returns the number of interior rings contained with the curve polygon.
const QgsCurve * exteriorRing() const
Returns the curve polygon's exterior ring.
const QgsCurve * interiorRing(int i) const
Retrieves an interior ring from the curve polygon.
Abstract base class for curved geometry type.
Definition: qgscurve.h:35
virtual int numPoints() const =0
Returns the number of points in the curve.
QgsCurve * segmentize(double tolerance=M_PI_2/90, SegmentationToleranceType toleranceType=MaximumAngle) const override
Returns a geometry without curves.
Definition: qgscurve.cpp:175
virtual bool isClosed() const
Returns true if the curve is closed.
Definition: qgscurve.cpp:53
QgsCurve * clone() const override=0
Clones the geometry by performing a deep copy.
virtual QgsPolygon * toPolygon(unsigned int segments=36) const
Returns a segmented polygon.
Definition: qgsellipse.cpp:224
Base class for feedback objects to be used for cancellation of something running in a worker thread.
Definition: qgsfeedback.h:44
bool isCanceled() const
Tells whether the operation has been canceled already.
Definition: qgsfeedback.h:53
Geometry collection.
int partCount() const override
Returns count of parts contained in the geometry.
int numGeometries() const
Returns the number of geometries within the collection.
const QgsAbstractGeometry * geometryN(int n) const
Returns a const reference to a geometry from within the collection.
static void pointOnLineWithDistance(double x1, double y1, double x2, double y2, double distance, double &x, double &y, double *z1=nullptr, double *z2=nullptr, double *z=nullptr, double *m1=nullptr, double *m2=nullptr, double *m=nullptr)
Calculates the point a specified distance from (x1, y1) toward a second point (x2,...
static double distance2D(double x1, double y1, double x2, double y2)
Returns the 2D distance between (x1, y1) and (x2, y2).
static double lineAngle(double x1, double y1, double x2, double y2)
Calculates the direction of line joining two points in radians, clockwise from the north direction.
static void weightedPointInTriangle(double aX, double aY, double bX, double bY, double cX, double cY, double weightB, double weightC, double &pointX, double &pointY)
Returns a weighted point inside the triangle denoted by the points (aX, aY), (bX, bY) and (cX,...
static double triangleArea(double aX, double aY, double bX, double bY, double cX, double cY)
Returns the area of the triangle denoted by the points (aX, aY), (bX, bY) and (cX,...
static int leftOfLine(const double x, const double y, const double x1, const double y1, const double x2, const double y2)
Returns a value < 0 if the point (x, y) is left of the line from (x1, y1) -> (x2, y2).
static void circleCenterRadius(const QgsPoint &pt1, const QgsPoint &pt2, const QgsPoint &pt3, double &radius, double &centerX, double &centerY)
Returns radius and center of the circle through pt1, pt2, pt3.
static bool pointContinuesArc(const QgsPoint &a1, const QgsPoint &a2, const QgsPoint &a3, const QgsPoint &b, double distanceTolerance, double pointSpacingAngleTolerance)
Returns true if point b is on the arc formed by points a1, a2, and a3, but not within that arc portio...
static int circleCircleOuterTangents(const QgsPointXY &center1, double radius1, const QgsPointXY &center2, double radius2, QgsPointXY &line1P1, QgsPointXY &line1P2, QgsPointXY &line2P1, QgsPointXY &line2P2)
Calculates the outer tangent points for two circles, centered at center1 and center2 and with radii o...
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:162
static QgsGeometry fromRect(const QgsRectangle &rect)
Creates a new geometry from a QgsRectangle.
static QgsGeometry collectGeometry(const QVector< QgsGeometry > &geometries)
Creates a new multipart geometry from a list of QgsGeometry objects.
QgsPointXY closestVertex(const QgsPointXY &point, int &closestVertexIndex, int &previousVertexIndex, int &nextVertexIndex, double &sqrDist) const
Returns the vertex closest to the given point, the corresponding vertex index, squared distance snap ...
Qgis::GeometryOperationResult addPart(const QVector< QgsPointXY > &points, Qgis::GeometryType geomType=Qgis::GeometryType::Unknown)
Adds a new part to a the geometry.
bool removeDuplicateNodes(double epsilon=4 *std::numeric_limits< double >::epsilon(), bool useZValues=false)
Removes duplicate nodes from the geometry, wherever removing the nodes does not result in a degenerat...
static QgsGeometry unaryUnion(const QVector< QgsGeometry > &geometries, const QgsGeometryParameters &parameters=QgsGeometryParameters())
Compute the unary union on a list of geometries.
Qgis::GeometryOperationResult rotate(double rotation, const QgsPointXY &center)
Rotate this geometry around the Z axis.
static QgsGeometryEngine * createGeometryEngine(const QgsAbstractGeometry *geometry, double precision=0.0)
Creates and returns a new geometry engine representing the specified geometry using precision on a gr...
Does vector analysis using the geos library and handles import, export, exception handling*.
Definition: qgsgeos.h:98
QgsGeometry triangularWavesRandomized(double minimumWavelength, double maximumWavelength, double minimumAmplitude, double maximumAmplitude, unsigned long seed=0) const
Constructs randomized triangular waves along the boundary of the geometry, with the specified wavelen...
QgsGeometry triangularWaves(double wavelength, double amplitude, bool strictWavelength=false) const
Constructs triangular waves along the boundary of the geometry, with the specified wavelength and amp...
QgsInternalGeometryEngine(const QgsGeometry &geometry)
The caller is responsible that the geometry is available and unchanged for the whole lifetime of this...
QgsGeometry roundWaves(double wavelength, double amplitude, bool strictWavelength=false) const
Constructs rounded (sine-like) waves along the boundary of the geometry, with the specified wavelengt...
QgsGeometry poleOfInaccessibility(double precision, double *distanceFromBoundary=nullptr) const
Calculates the approximate pole of inaccessibility for a surface, which is the most distant internal ...
QgsGeometry squareWaves(double wavelength, double amplitude, bool strictWavelength=false) const
Constructs square waves along the boundary of the geometry, with the specified wavelength and amplitu...
QgsGeometry variableWidthBufferByM(int segments) const
Calculates a variable width buffer using the m-values from a (multi)line geometry.
QgsGeometry extrude(double x, double y) const
Will extrude a line or (segmentized) curve by a given offset and return a polygon representation of i...
QgsGeometry roundWavesRandomized(double minimumWavelength, double maximumWavelength, double minimumAmplitude, double maximumAmplitude, unsigned long seed=0) const
Constructs randomized rounded (sine-like) waves along the boundary of the geometry,...
QgsGeometry orthogonalize(double tolerance=1.0E-8, int maxIterations=1000, double angleThreshold=15.0) const
Attempts to orthogonalize a line or polygon geometry by shifting vertices to make the geometries angl...
QgsGeometry variableWidthBuffer(int segments, const std::function< std::unique_ptr< double[] >(const QgsLineString *line) > &widthFunction) const
Calculates a variable width buffer for a (multi)curve geometry.
QString lastError() const
Returns an error string referring to the last error encountered.
QgsGeometry orientedMinimumBoundingBox(double &area, double &angle, double &width, double &height) const
Returns the oriented minimum bounding box for the geometry, which is the smallest (by area) rotated r...
QgsGeometry densifyByDistance(double distance) const
Densifies the geometry by adding regularly placed extra nodes inside each segment so that the maximum...
QgsGeometry taperedBuffer(double startWidth, double endWidth, int segments) const
Calculates a tapered width buffer for a (multi)curve geometry.
QVector< QgsPointXY > randomPointsInPolygon(int count, const std::function< bool(const QgsPointXY &) > &acceptPoint, unsigned long seed=0, QgsFeedback *feedback=nullptr, int maxTriesPerPoint=0)
Returns a list of count random points generated inside a polygon geometry (if acceptPoint is specifie...
QgsGeometry densifyByCount(int extraNodesPerSegment) const
Densifies the geometry by adding the specified number of extra nodes within each segment of the geome...
Applies a dash pattern to a geometry, returning a MultiLineString geometry which is the input geometr...
QgsGeometry squareWavesRandomized(double minimumWavelength, double maximumWavelength, double minimumAmplitude, double maximumAmplitude, unsigned long seed=0) const
Constructs randomized square waves along the boundary of the geometry, with the specified wavelength ...
QgsGeometry convertToCurves(double distanceTolerance, double angleTolerance) const
Attempts to convert a non-curved geometry into a curved geometry type (e.g.
bool isAxisParallelRectangle(double maximumDeviation, bool simpleRectanglesOnly=false) const
Returns true if the geometry is a polygon that is almost an axis-parallel rectangle.
Represents a single 2D line segment, consisting of a 2D start and end vertex only.
int pointLeftOfLine(const QgsPointXY &point) const
Tests if a point is to the left of the line segment.
double endY() const
Returns the segment's end y-coordinate.
QgsPointXY end() const
Returns the segment's end point.
double endX() const
Returns the segment's end x-coordinate.
QgsPointXY start() const
Returns the segment's start point.
double startX() const
Returns the segment's start x-coordinate.
void reverse()
Reverses the line segment, so that the start and end points are flipped.
double startY() const
Returns the segment's start y-coordinate.
Line string geometry type, with support for z-dimension and m-values.
Definition: qgslinestring.h:45
bool isClosed() const override
Returns true if the curve is closed.
const double * yData() const
Returns a const pointer to the y vertex data.
const double * xData() const
Returns a const pointer to the x vertex data.
double length() const override
Returns the planar, 2-dimensional length of the geometry.
QgsPoint startPoint() const override
Returns the starting point of the curve.
static QgsLineString * fromBezierCurve(const QgsPoint &start, const QgsPoint &controlPoint1, const QgsPoint &controlPoint2, const QgsPoint &end, int segments=30)
Returns a new linestring created by segmentizing the bezier curve between start and end,...
int numPoints() const override
Returns the number of points in the curve.
QgsPoint pointN(int i) const
Returns the specified point from inside the line string.
int nCoordinates() const override
Returns the number of nodes contained in the geometry.
double yAt(int index) const override
Returns the y-coordinate of the specified node in the line string.
void setYAt(int index, double y)
Sets the y-coordinate of the specified node in the line string.
double mAt(int index) const override
Returns the m value of the specified node in the line string.
void setXAt(int index, double x)
Sets the x-coordinate of the specified node in the line string.
QgsLineString * clone() const override
Clones the geometry by performing a deep copy.
double zAt(int index) const override
Returns the z-coordinate of the specified node in the line string.
double xAt(int index) const override
Returns the x-coordinate of the specified node in the line string.
Multi curve geometry collection.
Definition: qgsmulticurve.h:29
Multi polygon geometry collection.
Multi surface geometry collection.
A class to represent a 2D point.
Definition: qgspointxy.h:60
double y
Definition: qgspointxy.h:64
Definition: qgspointxy.h:63
Point geometry type, with support for z-dimension and m-values.
Definition: qgspoint.h:49
Definition: qgspoint.h:52
QgsPoint * clone() const override
Clones the geometry by performing a deep copy.
Definition: qgspoint.cpp:105
QgsPoint vertexAt(QgsVertexId) const override
Returns the point corresponding to a specified vertex id.
Definition: qgspoint.cpp:522
bool isEmpty() const override
Returns true if the geometry is empty.
Definition: qgspoint.cpp:733
double distance(double x, double y) const
Returns the Cartesian 2D distance between this point and a specified x, y coordinate.
Definition: qgspoint.h:393
double y
Definition: qgspoint.h:53
bool nextVertex(QgsVertexId &id, QgsPoint &vertex) const override
Returns next vertex id and coordinates.
Definition: qgspoint.cpp:473
Polygon geometry type.
Definition: qgspolygon.h:33
void setExteriorRing(QgsCurve *ring) override
Sets the exterior ring of the polygon.
Definition: qgspolygon.cpp:263
Adds an interior ring to the geometry (takes ownership)
Definition: qgspolygon.cpp:232
bool intersects(const QgsLineSegment2D &segment, QgsPointXY &intersectPoint) const
Finds the closest intersection point of the ray and a line segment.
A rectangle specified with double values.
Definition: qgsrectangle.h:42
double xMinimum() const
Returns the x minimum value (left side of rectangle).
Definition: qgsrectangle.h:201
double yMinimum() const
Returns the y minimum value (bottom side of rectangle).
Definition: qgsrectangle.h:211
double width() const
Returns the width of the rectangle.
Definition: qgsrectangle.h:236
double xMaximum() const
Returns the x maximum value (right side of rectangle).
Definition: qgsrectangle.h:196
double yMaximum() const
Returns the y maximum value (top side of rectangle).
Definition: qgsrectangle.h:206
double height() const
Returns the height of the rectangle.
Definition: qgsrectangle.h:243
Surface geometry type.
Definition: qgssurface.h:34
Class that takes care of tessellation of polygons into triangles.
QVector< float > data() const
Returns array of triangle vertex data.
void addPolygon(const QgsPolygon &polygon, float extrusionHeight)
Tessellates a triangle and adds its vertex entries to the output data array.
QString error() const
Returns a descriptive error string if the tessellation failed.
int dataVerticesCount() const
Returns the number of vertices stored in the output data array.
A class to represent a vector.
Definition: qgsvector.h:30
double lengthSquared() const
Returns the length of the vector.
Definition: qgsvector.h:134
double crossProduct(QgsVector v) const
Returns the 2D cross product of this vector and another vector v.
Definition: qgsvector.h:188
QgsVector normalized() const
Returns the vector's normalized (or "unit") vector (ie same angle but length of 1....
Definition: qgsvector.cpp:28
double length() const
Returns the length of the vector.
Definition: qgsvector.h:124
static Qgis::GeometryType geometryType(Qgis::WkbType type)
Returns the geometry type for a WKB type, e.g., both MultiPolygon and CurvePolygon would have a Polyg...
Definition: qgswkbtypes.h:862
static bool isCurvedType(Qgis::WkbType type)
Returns true if the WKB type is a curved type or can contain curved geometries.
Definition: qgswkbtypes.h:806
static Qgis::WkbType flatType(Qgis::WkbType type)
Returns the flat type for a WKB type.
Definition: qgswkbtypes.h:628
double ANALYSIS_EXPORT angle(QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, QgsPoint *p4)
Calculates the angle between two segments (in 2 dimension, z-values are ignored)
Definition: MathUtils.cpp:716
Contains geos related utilities and functions.
Definition: qgsgeos.h:36
std::unique_ptr< GEOSGeometry, GeosDeleter > unique_ptr
Scoped GEOS pointer.
Definition: qgsgeos.h:73
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
bool qgsDoubleNear(double a, double b, double epsilon=4 *std::numeric_limits< double >::epsilon())
Compare two doubles (but allow some difference)
Definition: qgis.h:5207
QVector< QgsPointXY > randomPointsInPolygonGeosBackend(const QgsAbstractGeometry *geometry, int count, const std::function< bool(const QgsPointXY &) > &acceptPoint, unsigned long seed, QgsFeedback *feedback, int maxTriesPerPoint, QString &error)
std::unique_ptr< QgsCurve > lineToCurve(const QgsCurve *curve, double distanceTolerance, double pointSpacingAngleTolerance)
bool matchesOrientation(std::array< Direction, 4 > dirs, std::array< Direction, 4 > oriented)
double normalizedDotProduct(const QgsPoint &a, const QgsPoint &b, const QgsPoint &c)
bool isClockwise(std::array< Direction, 4 > dirs)
Checks whether the 4 directions in dirs make up a clockwise rectangle.
std::unique_ptr< QgsLineString > squareWavesRandomizedAlongLine(const QgsLineString *line, const double minimumWavelength, const double maximumWavelength, const double minimumAmplitude, const double maximumAmplitude, std::uniform_real_distribution<> &uniformDist, std::mt19937 &mt)
std::unique_ptr< QgsLineString > squareWavesAlongLine(const QgsLineString *line, double wavelength, const double amplitude, const bool strictWavelength)
QgsAbstractGeometry * orthogonalizeGeom(const QgsAbstractGeometry *geom, int maxIterations, double tolerance, double lowerThreshold, double upperThreshold)
QVector< QgsPointXY > randomPointsInPolygonPoly2TriBackend(const QgsAbstractGeometry *geometry, int count, const std::function< bool(const QgsPointXY &) > &acceptPoint, unsigned long seed, QgsFeedback *feedback, int maxTriesPerPoint, QString &error)
std::unique_ptr< QgsLineString > triangularWavesRandomizedAlongLine(const QgsLineString *line, const double minimumWavelength, const double maximumWavelength, const double minimumAmplitude, const double maximumAmplitude, std::uniform_real_distribution<> &uniformDist, std::mt19937 &mt)
QgsLineString * doDensify(const QgsLineString *ring, int extraNodesPerSegment=-1, double distance=1)
Direction getEdgeDirection(const QgsPoint &p1, const QgsPoint &p2, double maxDev)
Determines the direction of an edge from p1 to p2.
std::unique_ptr< QgsAbstractGeometry > applyDashPatternPrivate(const QgsAbstractGeometry *geom, const QVector< double > &pattern, Qgis::DashPatternLineEndingRule startRule, Qgis::DashPatternLineEndingRule endRule, Qgis::DashPatternSizeAdjustment adjustment, double patternOffset)
std::unique_ptr< QgsLineString > roundWavesRandomizedAlongLine(const QgsLineString *line, const double minimumWavelength, const double maximumWavelength, const double minimumAmplitude, const double maximumAmplitude, std::uniform_real_distribution<> &uniformDist, std::mt19937 &mt)
QgsVector calcMotion(const QgsPoint &a, const QgsPoint &b, const QgsPoint &c, double lowerThreshold, double upperThreshold)
QVector< QgsPointXY > generateSegmentCurve(const QgsPoint &center1, const double radius1, const QgsPoint &center2, const double radius2)
bool dotProductWithinAngleTolerance(double dotProduct, double lowerThreshold, double upperThreshold)
QgsLineString * doOrthogonalize(QgsLineString *ring, int iterations, double tolerance, double lowerThreshold, double upperThreshold)
double squareness(QgsLineString *ring, double lowerThreshold, double upperThreshold)
std::unique_ptr< QgsMultiLineString > dashPatternAlongLine(const QgsLineString *line, const QVector< double > &pattern, Qgis::DashPatternLineEndingRule startRule, Qgis::DashPatternLineEndingRule endRule, Qgis::DashPatternSizeAdjustment adjustment, double patternOffset)
std::unique_ptr< QgsAbstractGeometry > convertGeometryToCurves(const QgsAbstractGeometry *geom, double distanceTolerance, double angleTolerance)
std::unique_ptr< QgsAbstractGeometry > roundWavesPrivate(const QgsAbstractGeometry *geom, double wavelength, double amplitude, bool strictWavelength)
bool isCounterClockwise(std::array< Direction, 4 > dirs)
Checks whether the 4 directions in dirs make up a counter-clockwise rectangle.
QgsAbstractGeometry * densifyGeometry(const QgsAbstractGeometry *geom, int extraNodesPerSegment=1, double distance=1)
std::unique_ptr< QgsLineString > roundWavesAlongLine(const QgsLineString *line, double wavelength, const double amplitude, const bool strictWavelength)
std::pair< bool, std::array< Direction, 4 > > getEdgeDirections(const QgsPolygon *g, double maxDev)
Checks whether the polygon consists of four nearly axis-parallel sides.
std::unique_ptr< QgsAbstractGeometry > squareWavesRandomizedPrivate(const QgsAbstractGeometry *geom, double minimumWavelength, double maximumWavelength, double minimumAmplitude, double maximumAmplitude, std::uniform_real_distribution<> &uniformDist, std::mt19937 &mt)
std::unique_ptr< QgsAbstractGeometry > triangularWavesRandomizedPrivate(const QgsAbstractGeometry *geom, double minimumWavelength, double maximumWavelength, double minimumAmplitude, double maximumAmplitude, std::uniform_real_distribution<> &uniformDist, std::mt19937 &mt)
std::unique_ptr< QgsAbstractGeometry > triangularWavesPrivate(const QgsAbstractGeometry *geom, double wavelength, double amplitude, bool strictWavelength)
std::unique_ptr< QgsAbstractGeometry > squareWavesPrivate(const QgsAbstractGeometry *geom, double wavelength, double amplitude, bool strictWavelength)
std::unique_ptr< QgsLineString > triangularWavesAlongLine(const QgsLineString *line, double wavelength, const double amplitude, const bool strictWavelength)
std::unique_ptr< QgsAbstractGeometry > roundWavesRandomizedPrivate(const QgsAbstractGeometry *geom, double minimumWavelength, double maximumWavelength, double minimumAmplitude, double maximumAmplitude, std::uniform_real_distribution<> &uniformDist, std::mt19937 &mt)
QLineF segment(int index, QRectF rect, double radius)
int precision
Utility class for identifying a unique vertex within a geometry.
Definition: qgsvertexid.h:30