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