QGIS API Documentation  3.25.0-Master (dec16ba68b)
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
7  email : [email protected]
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 
51 enum class Direction
52 {
53  Up,
54  Right,
55  Down,
56  Left,
57  None
58 };
59 
65 Direction 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 
92 std::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 
130 bool 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 
144 bool 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 
154 bool 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 
161 bool 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 
249 class 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 
272 struct GreaterThanByMax
273 {
274  bool operator()( const Cell *lhs, const Cell *rhs ) const
275  {
276  return rhs->max > lhs->max;
277  }
278 };
279 
280 Cell *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 
305 QgsPoint 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 
380 QgsGeometry 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 
440 bool dotProductWithinAngleTolerance( double dotProduct, double lowerThreshold, double upperThreshold )
441 {
442  return lowerThreshold > std::fabs( dotProduct ) || std::fabs( dotProduct ) > upperThreshold;
443 }
444 
445 double 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 
458 double 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 
502 QgsVector 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  // 0.1 magic number from JOSM implementation - think this is to limit each iterative step
530  return new_v.normalized() * ( 0.1 * dotProduct * scale );
531 }
532 
533 QgsLineString *doOrthogonalize( QgsLineString *ring, int iterations, double tolerance, double lowerThreshold, double upperThreshold )
534 {
535  double minScore = std::numeric_limits<double>::max();
536 
537  bool isClosed = ring->isClosed();
538  int numPoints = ring->numPoints();
539 
540  std::unique_ptr< QgsLineString > best( ring->clone() );
541 
542  QVector< QgsVector > /* yep */ motions;
543  motions.reserve( numPoints );
544 
545  for ( int it = 0; it < iterations; ++it )
546  {
547  motions.resize( 0 ); // avoid re-allocations
548 
549  // first loop through an calculate all motions
550  QgsPoint a;
551  QgsPoint b;
552  QgsPoint c;
553  for ( int i = 0; i < numPoints; ++i )
554  {
555  if ( isClosed && i == numPoints - 1 )
556  motions << motions.at( 0 ); //closed ring, so last point follows first point motion
557  else if ( !isClosed && ( i == 0 || i == numPoints - 1 ) )
558  {
559  b = ring->pointN( 0 );
560  c = ring->pointN( 1 );
561  motions << QgsVector( 0, 0 ); //non closed line, leave first/last vertex untouched
562  }
563  else
564  {
565  if ( i == 0 )
566  {
567  a = ring->pointN( numPoints - 1 );
568  b = ring->pointN( 0 );
569  }
570  if ( i == numPoints - 1 )
571  c = ring->pointN( 0 );
572  else
573  c = ring->pointN( i + 1 );
574 
575  motions << calcMotion( a, b, c, lowerThreshold, upperThreshold );
576  }
577  a = b;
578  b = c;
579  }
580 
581  // then apply them
582  for ( int i = 0; i < numPoints; ++i )
583  {
584  ring->setXAt( i, ring->xAt( i ) + motions.at( i ).x() );
585  ring->setYAt( i, ring->yAt( i ) + motions.at( i ).y() );
586  }
587 
588  double newScore = squareness( ring, lowerThreshold, upperThreshold );
589  if ( newScore < minScore )
590  {
591  best.reset( ring->clone() );
592  minScore = newScore;
593  }
594 
595  if ( minScore < tolerance )
596  break;
597  }
598 
599  delete ring;
600 
601  return best.release();
602 }
603 
604 
605 QgsAbstractGeometry *orthogonalizeGeom( const QgsAbstractGeometry *geom, int maxIterations, double tolerance, double lowerThreshold, double upperThreshold )
606 {
607  std::unique_ptr< QgsAbstractGeometry > segmentizedCopy;
608  if ( QgsWkbTypes::isCurvedType( geom->wkbType() ) )
609  {
610  segmentizedCopy.reset( geom->segmentize() );
611  geom = segmentizedCopy.get();
612  }
613 
615  {
616  return doOrthogonalize( static_cast< QgsLineString * >( geom->clone() ),
617  maxIterations, tolerance, lowerThreshold, upperThreshold );
618  }
619  else
620  {
621  // polygon
622  const QgsPolygon *polygon = static_cast< const QgsPolygon * >( geom );
623  QgsPolygon *result = new QgsPolygon();
624 
625  result->setExteriorRing( doOrthogonalize( static_cast< QgsLineString * >( polygon->exteriorRing()->clone() ),
626  maxIterations, tolerance, lowerThreshold, upperThreshold ) );
627  for ( int i = 0; i < polygon->numInteriorRings(); ++i )
628  {
629  result->addInteriorRing( doOrthogonalize( static_cast< QgsLineString * >( polygon->interiorRing( i )->clone() ),
630  maxIterations, tolerance, lowerThreshold, upperThreshold ) );
631  }
632 
633  return result;
634  }
635 }
636 
637 QgsGeometry QgsInternalGeometryEngine::orthogonalize( double tolerance, int maxIterations, double angleThreshold ) const
638 {
639  mLastError.clear();
640  if ( !mGeometry || ( QgsWkbTypes::geometryType( mGeometry->wkbType() ) != QgsWkbTypes::LineGeometry
642  {
643  return QgsGeometry();
644  }
645 
646  double lowerThreshold = std::cos( ( 90 - angleThreshold ) * M_PI / 180.00 );
647  double upperThreshold = std::cos( angleThreshold * M_PI / 180.0 );
648 
649  if ( const QgsGeometryCollection *gc = qgsgeometry_cast< const QgsGeometryCollection *>( mGeometry ) )
650  {
651  int numGeom = gc->numGeometries();
652  QVector< QgsAbstractGeometry * > geometryList;
653  geometryList.reserve( numGeom );
654  for ( int i = 0; i < numGeom; ++i )
655  {
656  geometryList << orthogonalizeGeom( gc->geometryN( i ), maxIterations, tolerance, lowerThreshold, upperThreshold );
657  }
658 
659  QgsGeometry first = QgsGeometry( geometryList.takeAt( 0 ) );
660  for ( QgsAbstractGeometry *g : std::as_const( geometryList ) )
661  {
662  first.addPart( g );
663  }
664  return first;
665  }
666  else
667  {
668  return QgsGeometry( orthogonalizeGeom( mGeometry, maxIterations, tolerance, lowerThreshold, upperThreshold ) );
669  }
670 }
671 
672 // if extraNodesPerSegment < 0, then use distance based mode
673 QgsLineString *doDensify( const QgsLineString *ring, int extraNodesPerSegment = -1, double distance = 1 )
674 {
675  QVector< double > outX;
676  QVector< double > outY;
677  QVector< double > outZ;
678  QVector< double > outM;
679  double multiplier = 1.0 / double( extraNodesPerSegment + 1 );
680 
681  int nPoints = ring->numPoints();
682  outX.reserve( ( extraNodesPerSegment + 1 ) * nPoints );
683  outY.reserve( ( extraNodesPerSegment + 1 ) * nPoints );
684  bool withZ = ring->is3D();
685  if ( withZ )
686  outZ.reserve( ( extraNodesPerSegment + 1 ) * nPoints );
687  bool withM = ring->isMeasure();
688  if ( withM )
689  outM.reserve( ( extraNodesPerSegment + 1 ) * nPoints );
690  double x1 = 0;
691  double x2 = 0;
692  double y1 = 0;
693  double y2 = 0;
694  double z1 = 0;
695  double z2 = 0;
696  double m1 = 0;
697  double m2 = 0;
698  double xOut = 0;
699  double yOut = 0;
700  double zOut = 0;
701  double mOut = 0;
702  int extraNodesThisSegment = extraNodesPerSegment;
703  for ( int i = 0; i < nPoints - 1; ++i )
704  {
705  x1 = ring->xAt( i );
706  x2 = ring->xAt( i + 1 );
707  y1 = ring->yAt( i );
708  y2 = ring->yAt( i + 1 );
709  if ( withZ )
710  {
711  z1 = ring->zAt( i );
712  z2 = ring->zAt( i + 1 );
713  }
714  if ( withM )
715  {
716  m1 = ring->mAt( i );
717  m2 = ring->mAt( i + 1 );
718  }
719 
720  outX << x1;
721  outY << y1;
722  if ( withZ )
723  outZ << z1;
724  if ( withM )
725  outM << m1;
726 
727  if ( extraNodesPerSegment < 0 )
728  {
729  // distance mode
730  extraNodesThisSegment = std::floor( std::sqrt( ( x2 - x1 ) * ( x2 - x1 ) + ( y2 - y1 ) * ( y2 - y1 ) ) / distance );
731  if ( extraNodesThisSegment >= 1 )
732  multiplier = 1.0 / ( extraNodesThisSegment + 1 );
733  }
734 
735  for ( int j = 0; j < extraNodesThisSegment; ++j )
736  {
737  double delta = multiplier * ( j + 1 );
738  xOut = x1 + delta * ( x2 - x1 );
739  yOut = y1 + delta * ( y2 - y1 );
740  if ( withZ )
741  zOut = z1 + delta * ( z2 - z1 );
742  if ( withM )
743  mOut = m1 + delta * ( m2 - m1 );
744 
745  outX << xOut;
746  outY << yOut;
747  if ( withZ )
748  outZ << zOut;
749  if ( withM )
750  outM << mOut;
751  }
752  }
753  outX << ring->xAt( nPoints - 1 );
754  outY << ring->yAt( nPoints - 1 );
755  if ( withZ )
756  outZ << ring->zAt( nPoints - 1 );
757  if ( withM )
758  outM << ring->mAt( nPoints - 1 );
759 
760  QgsLineString *result = new QgsLineString( outX, outY, outZ, outM );
761  return result;
762 }
763 
764 QgsAbstractGeometry *densifyGeometry( const QgsAbstractGeometry *geom, int extraNodesPerSegment = 1, double distance = 1 )
765 {
766  std::unique_ptr< QgsAbstractGeometry > segmentizedCopy;
767  if ( QgsWkbTypes::isCurvedType( geom->wkbType() ) )
768  {
769  segmentizedCopy.reset( geom->segmentize() );
770  geom = segmentizedCopy.get();
771  }
772 
774  {
775  return doDensify( static_cast< const QgsLineString * >( geom ), extraNodesPerSegment, distance );
776  }
777  else
778  {
779  // polygon
780  const QgsPolygon *polygon = static_cast< const QgsPolygon * >( geom );
781  QgsPolygon *result = new QgsPolygon();
782 
783  result->setExteriorRing( doDensify( static_cast< const QgsLineString * >( polygon->exteriorRing() ),
784  extraNodesPerSegment, distance ) );
785  for ( int i = 0; i < polygon->numInteriorRings(); ++i )
786  {
787  result->addInteriorRing( doDensify( static_cast< const QgsLineString * >( polygon->interiorRing( i ) ),
788  extraNodesPerSegment, distance ) );
789  }
790 
791  return result;
792  }
793 }
794 
796 {
797  mLastError.clear();
798  if ( !mGeometry )
799  {
800  return QgsGeometry();
801  }
802 
804  {
805  return QgsGeometry( mGeometry->clone() ); // point geometry, nothing to do
806  }
807 
808  if ( const QgsGeometryCollection *gc = qgsgeometry_cast< const QgsGeometryCollection *>( mGeometry ) )
809  {
810  int numGeom = gc->numGeometries();
811  QVector< QgsAbstractGeometry * > geometryList;
812  geometryList.reserve( numGeom );
813  for ( int i = 0; i < numGeom; ++i )
814  {
815  geometryList << densifyGeometry( gc->geometryN( i ), extraNodesPerSegment );
816  }
817 
818  QgsGeometry first = QgsGeometry( geometryList.takeAt( 0 ) );
819  for ( QgsAbstractGeometry *g : std::as_const( geometryList ) )
820  {
821  first.addPart( g );
822  }
823  return first;
824  }
825  else
826  {
827  return QgsGeometry( densifyGeometry( mGeometry, extraNodesPerSegment ) );
828  }
829 }
830 
832 {
833  mLastError.clear();
834  if ( !mGeometry )
835  {
836  return QgsGeometry();
837  }
838 
840  {
841  return QgsGeometry( mGeometry->clone() ); // point geometry, nothing to do
842  }
843 
844  if ( const QgsGeometryCollection *gc = qgsgeometry_cast< const QgsGeometryCollection *>( mGeometry ) )
845  {
846  int numGeom = gc->numGeometries();
847  QVector< QgsAbstractGeometry * > geometryList;
848  geometryList.reserve( numGeom );
849  for ( int i = 0; i < numGeom; ++i )
850  {
851  geometryList << densifyGeometry( gc->geometryN( i ), -1, distance );
852  }
853 
854  QgsGeometry first = QgsGeometry( geometryList.takeAt( 0 ) );
855  for ( QgsAbstractGeometry *g : std::as_const( geometryList ) )
856  {
857  first.addPart( g );
858  }
859  return first;
860  }
861  else
862  {
863  return QgsGeometry( densifyGeometry( mGeometry, -1, distance ) );
864  }
865 }
866 
868 //
869 // QgsLineSegmentDistanceComparer
870 //
871 
872 // adapted for QGIS geometry classes from original work at https://github.com/trylock/visibility by trylock
873 bool QgsLineSegmentDistanceComparer::operator()( QgsLineSegment2D ab, QgsLineSegment2D cd ) const
874 {
875  Q_ASSERT_X( ab.pointLeftOfLine( mOrigin ) != 0,
876  "line_segment_dist_comparer",
877  "AB must not be collinear with the origin." );
878  Q_ASSERT_X( cd.pointLeftOfLine( mOrigin ) != 0,
879  "line_segment_dist_comparer",
880  "CD must not be collinear with the origin." );
881 
882  // flip the segments so that if there are common endpoints,
883  // they will be the segment's start points
884  if ( ab.end() == cd.start() || ab.end() == cd.end() )
885  ab.reverse();
886  if ( ab.start() == cd.end() )
887  cd.reverse();
888 
889  // cases with common endpoints
890  if ( ab.start() == cd.start() )
891  {
892  const int oad = QgsGeometryUtils::leftOfLine( cd.endX(), cd.endY(), mOrigin.x(), mOrigin.y(), ab.startX(), ab.startY() );
893  const int oab = ab.pointLeftOfLine( mOrigin );
894  if ( ab.end() == cd.end() || oad != oab )
895  return false;
896  else
897  return ab.pointLeftOfLine( cd.end() ) != oab;
898  }
899  else
900  {
901  // cases without common endpoints
902  const int cda = cd.pointLeftOfLine( ab.start() );
903  const int cdb = cd.pointLeftOfLine( ab.end() );
904  if ( cdb == 0 && cda == 0 )
905  {
906  return mOrigin.sqrDist( ab.start() ) < mOrigin.sqrDist( cd.start() );
907  }
908  else if ( cda == cdb || cda == 0 || cdb == 0 )
909  {
910  const int cdo = cd.pointLeftOfLine( mOrigin );
911  return cdo == cda || cdo == cdb;
912  }
913  else
914  {
915  const int abo = ab.pointLeftOfLine( mOrigin );
916  return abo != ab.pointLeftOfLine( cd.start() );
917  }
918  }
919 }
920 
921 //
922 // QgsClockwiseAngleComparer
923 //
924 
925 bool QgsClockwiseAngleComparer::operator()( const QgsPointXY &a, const QgsPointXY &b ) const
926 {
927  const bool aIsLeft = a.x() < mVertex.x();
928  const bool bIsLeft = b.x() < mVertex.x();
929  if ( aIsLeft != bIsLeft )
930  return bIsLeft;
931 
932  if ( qgsDoubleNear( a.x(), mVertex.x() ) && qgsDoubleNear( b.x(), mVertex.x() ) )
933  {
934  if ( a.y() >= mVertex.y() || b.y() >= mVertex.y() )
935  {
936  return b.y() < a.y();
937  }
938  else
939  {
940  return a.y() < b.y();
941  }
942  }
943  else
944  {
945  const QgsVector oa = a - mVertex;
946  const QgsVector ob = b - mVertex;
947  const double det = oa.crossProduct( ob );
948  if ( qgsDoubleNear( det, 0.0 ) )
949  {
950  return oa.lengthSquared() < ob.lengthSquared();
951  }
952  else
953  {
954  return det < 0;
955  }
956  }
957 }
958 
960 
961 //
962 // QgsRay2D
963 //
964 
965 bool QgsRay2D::intersects( const QgsLineSegment2D &segment, QgsPointXY &intersectPoint ) const
966 {
967  const QgsVector ao = origin - segment.start();
968  const QgsVector ab = segment.end() - segment.start();
969  const double det = ab.crossProduct( direction );
970  if ( qgsDoubleNear( det, 0.0 ) )
971  {
972  const int abo = segment.pointLeftOfLine( origin );
973  if ( abo != 0 )
974  {
975  return false;
976  }
977  else
978  {
979  const double distA = ao * direction;
980  const double distB = ( origin - segment.end() ) * direction;
981 
982  if ( distA > 0 && distB > 0 )
983  {
984  return false;
985  }
986  else
987  {
988  if ( ( distA > 0 ) != ( distB > 0 ) )
989  intersectPoint = origin;
990  else if ( distA > distB ) // at this point, both distances are negative
991  intersectPoint = segment.start(); // hence the nearest point is A
992  else
993  intersectPoint = segment.end();
994  return true;
995  }
996  }
997  }
998  else
999  {
1000  const double u = ao.crossProduct( direction ) / det;
1001  if ( u < 0.0 || 1.0 < u )
1002  {
1003  return false;
1004  }
1005  else
1006  {
1007  const double t = -ab.crossProduct( ao ) / det;
1008  intersectPoint = origin + direction * t;
1009  return qgsDoubleNear( t, 0.0 ) || t > 0;
1010  }
1011  }
1012 }
1013 
1014 QVector<QgsPointXY> generateSegmentCurve( const QgsPoint &center1, const double radius1, const QgsPoint &center2, const double radius2 )
1015 {
1016  // ensure that first circle is smaller than second
1017  if ( radius1 > radius2 )
1018  return generateSegmentCurve( center2, radius2, center1, radius1 );
1019 
1020  QgsPointXY t1;
1021  QgsPointXY t2;
1022  QgsPointXY t3;
1023  QgsPointXY t4;
1024  QVector<QgsPointXY> points;
1025  if ( QgsGeometryUtils::circleCircleOuterTangents( center1, radius1, center2, radius2, t1, t2, t3, t4 ) )
1026  {
1027  points << t1
1028  << t2
1029  << t4
1030  << t3;
1031  }
1032  return points;
1033 }
1034 
1035 // partially ported from JTS VariableWidthBuffer,
1036 // https://github.com/topobyte/jts/blob/master/jts-lab/src/main/java/com/vividsolutions/jts/operation/buffer/VariableWidthBuffer.java
1037 
1038 QgsGeometry QgsInternalGeometryEngine::variableWidthBuffer( int segments, const std::function< std::unique_ptr< double[] >( const QgsLineString *line ) > &widthFunction ) const
1039 {
1040  mLastError.clear();
1041  if ( !mGeometry )
1042  {
1043  return QgsGeometry();
1044  }
1045 
1046  std::vector< std::unique_ptr<QgsLineString > > linesToProcess;
1047 
1048  const QgsMultiCurve *multiCurve = qgsgeometry_cast< const QgsMultiCurve * >( mGeometry );
1049  if ( multiCurve )
1050  {
1051  for ( int i = 0; i < multiCurve->partCount(); ++i )
1052  {
1053  if ( static_cast< const QgsCurve * >( multiCurve->geometryN( i ) )->nCoordinates() == 0 )
1054  continue; // skip 0 length lines
1055 
1056  linesToProcess.emplace_back( static_cast<QgsLineString *>( multiCurve->geometryN( i )->clone() ) );
1057  }
1058  }
1059 
1060  const QgsCurve *curve = qgsgeometry_cast< const QgsCurve * >( mGeometry );
1061  if ( curve )
1062  {
1063  if ( curve->nCoordinates() > 0 )
1064  linesToProcess.emplace_back( static_cast<QgsLineString *>( curve->segmentize() ) );
1065  }
1066 
1067  if ( linesToProcess.empty() )
1068  {
1069  QgsGeometry g;
1070  g.mLastError = QStringLiteral( "Input geometry was not a curve type geometry" );
1071  return g;
1072  }
1073 
1074  QVector<QgsGeometry> bufferedLines;
1075  bufferedLines.reserve( linesToProcess.size() );
1076 
1077  for ( std::unique_ptr< QgsLineString > &line : linesToProcess )
1078  {
1079  QVector<QgsGeometry> parts;
1080  QgsPoint prevPoint;
1081  double prevRadius = 0;
1082  QgsGeometry prevCircle;
1083 
1084  std::unique_ptr< double[] > widths = widthFunction( line.get() ) ;
1085  for ( int i = 0; i < line->nCoordinates(); ++i )
1086  {
1087  QgsPoint thisPoint = line->pointN( i );
1088  QgsGeometry thisCircle;
1089  double thisRadius = widths[ i ] / 2.0;
1090  if ( thisRadius > 0 )
1091  {
1092  QgsGeometry p = QgsGeometry( thisPoint.clone() );
1093 
1094  QgsCircle circ( thisPoint, thisRadius );
1095  thisCircle = QgsGeometry( circ.toPolygon( segments * 4 ) );
1096  parts << thisCircle;
1097  }
1098  else
1099  {
1100  thisCircle = QgsGeometry( thisPoint.clone() );
1101  }
1102 
1103  if ( i > 0 )
1104  {
1105  if ( prevRadius > 0 || thisRadius > 0 )
1106  {
1107  QVector< QgsPointXY > points = generateSegmentCurve( prevPoint, prevRadius, thisPoint, thisRadius );
1108  if ( !points.empty() )
1109  {
1110  // snap points to circle vertices
1111 
1112  int atVertex = 0;
1113  int beforeVertex = 0;
1114  int afterVertex = 0;
1115  double sqrDist = 0;
1116  double sqrDistPrev = 0;
1117  for ( int j = 0; j < points.count(); ++j )
1118  {
1119  QgsPointXY pA = prevCircle.closestVertex( points.at( j ), atVertex, beforeVertex, afterVertex, sqrDistPrev );
1120  QgsPointXY pB = thisCircle.closestVertex( points.at( j ), atVertex, beforeVertex, afterVertex, sqrDist );
1121  points[j] = sqrDistPrev < sqrDist ? pA : pB;
1122  }
1123  // close ring
1124  points.append( points.at( 0 ) );
1125 
1126  std::unique_ptr< QgsPolygon > poly = std::make_unique< QgsPolygon >();
1127  poly->setExteriorRing( new QgsLineString( points ) );
1128  if ( poly->area() > 0 )
1129  parts << QgsGeometry( std::move( poly ) );
1130  }
1131  }
1132  }
1133  prevPoint = thisPoint;
1134  prevRadius = thisRadius;
1135  prevCircle = thisCircle;
1136  }
1137 
1138  bufferedLines << QgsGeometry::unaryUnion( parts );
1139  }
1140 
1141  return QgsGeometry::collectGeometry( bufferedLines );
1142 }
1143 
1144 QgsGeometry QgsInternalGeometryEngine::taperedBuffer( double start, double end, int segments ) const
1145 {
1146  mLastError.clear();
1147  start = std::fabs( start );
1148  end = std::fabs( end );
1149 
1150  auto interpolateWidths = [ start, end ]( const QgsLineString * line )->std::unique_ptr< double [] >
1151  {
1152  // ported from JTS VariableWidthBuffer,
1153  // https://github.com/topobyte/jts/blob/master/jts-lab/src/main/java/com/vividsolutions/jts/operation/buffer/VariableWidthBuffer.java
1154  std::unique_ptr< double [] > widths( new double[ line->nCoordinates() ] );
1155  widths[0] = start;
1156  widths[line->nCoordinates() - 1] = end;
1157 
1158  double lineLength = line->length();
1159  double currentLength = 0;
1160  QgsPoint prevPoint = line->pointN( 0 );
1161  for ( int i = 1; i < line->nCoordinates() - 1; ++i )
1162  {
1163  QgsPoint point = line->pointN( i );
1164  double segmentLength = point.distance( prevPoint );
1165  currentLength += segmentLength;
1166  double lengthFraction = lineLength > 0 ? currentLength / lineLength : 1;
1167  double delta = lengthFraction * ( end - start );
1168  widths[i] = start + delta;
1169  prevPoint = point;
1170  }
1171  return widths;
1172  };
1173 
1174  return variableWidthBuffer( segments, interpolateWidths );
1175 }
1176 
1178 {
1179  mLastError.clear();
1180  auto widthByM = []( const QgsLineString * line )->std::unique_ptr< double [] >
1181  {
1182  std::unique_ptr< double [] > widths( new double[ line->nCoordinates() ] );
1183  for ( int i = 0; i < line->nCoordinates(); ++i )
1184  {
1185  widths[ i ] = line->mAt( i );
1186  }
1187  return widths;
1188  };
1189 
1190  return variableWidthBuffer( segments, widthByM );
1191 }
1192 
1193 QVector<QgsPointXY> QgsInternalGeometryEngine::randomPointsInPolygon( const QgsGeometry &polygon, int count,
1194  const std::function< bool( const QgsPointXY & ) > &acceptPoint, unsigned long seed, QgsFeedback *feedback, int maxTriesPerPoint )
1195 {
1196  if ( polygon.type() != QgsWkbTypes::PolygonGeometry || count == 0 )
1197  return QVector< QgsPointXY >();
1198 
1199  // step 1 - tessellate the polygon to triangles
1200  QgsRectangle bounds = polygon.boundingBox();
1201  QgsTessellator t( bounds, false, false, false, true );
1202 
1203  if ( polygon.isMultipart() )
1204  {
1205  const QgsMultiSurface *ms = qgsgeometry_cast< const QgsMultiSurface * >( polygon.constGet() );
1206  for ( int i = 0; i < ms->numGeometries(); ++i )
1207  {
1208  if ( feedback && feedback->isCanceled() )
1209  return QVector< QgsPointXY >();
1210 
1211  if ( QgsPolygon *poly = qgsgeometry_cast< QgsPolygon * >( ms->geometryN( i ) ) )
1212  {
1213  t.addPolygon( *poly, 0 );
1214  }
1215  else
1216  {
1217  std::unique_ptr< QgsPolygon > p( qgsgeometry_cast< QgsPolygon * >( ms->geometryN( i )->segmentize() ) );
1218  t.addPolygon( *p, 0 );
1219  }
1220  }
1221  }
1222  else
1223  {
1224  if ( const QgsPolygon *poly = qgsgeometry_cast< const QgsPolygon * >( polygon.constGet() ) )
1225  {
1226  t.addPolygon( *poly, 0 );
1227  }
1228  else
1229  {
1230  std::unique_ptr< QgsPolygon > p( qgsgeometry_cast< QgsPolygon * >( polygon.constGet()->segmentize() ) );
1231  t.addPolygon( *p, 0 );
1232  }
1233  }
1234 
1235  if ( feedback && feedback->isCanceled() )
1236  return QVector< QgsPointXY >();
1237 
1238  const QVector<float> triangleData = t.data();
1239  if ( triangleData.empty() )
1240  return QVector< QgsPointXY >(); //hm
1241 
1242  // calculate running sum of triangle areas
1243  std::vector< double > cumulativeAreas;
1244  cumulativeAreas.reserve( t.dataVerticesCount() / 3 );
1245  double totalArea = 0;
1246  for ( auto it = triangleData.constBegin(); it != triangleData.constEnd(); )
1247  {
1248  if ( feedback && feedback->isCanceled() )
1249  return QVector< QgsPointXY >();
1250 
1251  const float aX = *it++;
1252  ( void )it++; // z
1253  const float aY = -( *it++ );
1254  const float bX = *it++;
1255  ( void )it++; // z
1256  const float bY = -( *it++ );
1257  const float cX = *it++;
1258  ( void )it++; // z
1259  const float cY = -( *it++ );
1260 
1261  const double area = QgsGeometryUtils::triangleArea( aX, aY, bX, bY, cX, cY );
1262  totalArea += area;
1263  cumulativeAreas.emplace_back( totalArea );
1264  }
1265 
1266  std::random_device rd;
1267  std::mt19937 mt( seed == 0 ? rd() : seed );
1268  std::uniform_real_distribution<> uniformDist( 0, 1 );
1269 
1270  // selects a random triangle, weighted by triangle area
1271  auto selectRandomTriangle = [&cumulativeAreas, totalArea]( double random )->int
1272  {
1273  int triangle = 0;
1274  const double target = random * totalArea;
1275  for ( auto it = cumulativeAreas.begin(); it != cumulativeAreas.end(); ++it, triangle++ )
1276  {
1277  if ( *it > target )
1278  return triangle;
1279  }
1280  Q_ASSERT_X( false, "QgsInternalGeometryEngine::randomPointsInPolygon", "Invalid random triangle index" );
1281  return 0; // no warnings
1282  };
1283 
1284 
1285  QVector<QgsPointXY> result;
1286  result.reserve( count );
1287  int tries = 0;
1288  for ( int i = 0; i < count; )
1289  {
1290  if ( feedback && feedback->isCanceled() )
1291  return QVector< QgsPointXY >();
1292 
1293  const double triangleIndexRnd = uniformDist( mt );
1294  // pick random triangle, weighted by triangle area
1295  const int triangleIndex = selectRandomTriangle( triangleIndexRnd );
1296 
1297  // generate a random point inside this triangle
1298  const double weightB = uniformDist( mt );
1299  const double weightC = uniformDist( mt );
1300  double x;
1301  double y;
1302 
1303  // get triangle
1304  const double aX = triangleData.at( triangleIndex * 9 ) + bounds.xMinimum();
1305  const double aY = -triangleData.at( triangleIndex * 9 + 2 ) + bounds.yMinimum();
1306  const double bX = triangleData.at( triangleIndex * 9 + 3 ) + bounds.xMinimum();
1307  const double bY = -triangleData.at( triangleIndex * 9 + 5 ) + bounds.yMinimum();
1308  const double cX = triangleData.at( triangleIndex * 9 + 6 ) + bounds.xMinimum();
1309  const double cY = -triangleData.at( triangleIndex * 9 + 8 ) + bounds.yMinimum();
1310  QgsGeometryUtils::weightedPointInTriangle( aX, aY, bX, bY, cX, cY, weightB, weightC, x, y );
1311 
1312  QgsPointXY candidate( x, y );
1313  if ( acceptPoint( candidate ) )
1314  {
1315  result << QgsPointXY( x, y );
1316  i++;
1317  tries = 0;
1318  }
1319  else if ( maxTriesPerPoint != 0 )
1320  {
1321  tries++;
1322  // Skip this point if maximum tries is reached
1323  if ( tries == maxTriesPerPoint )
1324  {
1325  tries = 0;
1326  i++;
1327  }
1328  }
1329  }
1330  return result;
1331 }
1332 
1333 // ported from PostGIS' lwgeom pta_unstroke
1334 
1335 std::unique_ptr< QgsCompoundCurve > lineToCurve( const QgsLineString *lineString, double distanceTolerance,
1336  double pointSpacingAngleTolerance )
1337 {
1338  std::unique_ptr< QgsCompoundCurve > out = std::make_unique< QgsCompoundCurve >();
1339 
1340  /* Minimum number of edges, per quadrant, required to define an arc */
1341  const unsigned int minQuadEdges = 2;
1342 
1343  /* Die on null input */
1344  if ( !lineString )
1345  return nullptr;
1346 
1347  /* Null on empty input? */
1348  if ( lineString->nCoordinates() == 0 )
1349  return nullptr;
1350 
1351  /* We can't desegmentize anything shorter than four points */
1352  if ( lineString->nCoordinates() < 4 )
1353  {
1354  out->addCurve( lineString->clone() );
1355  return out;
1356  }
1357 
1358  /* Allocate our result array of vertices that are part of arcs */
1359  int numEdges = lineString->nCoordinates() - 1;
1360  QVector< int > edgesInArcs( numEdges + 1, 0 );
1361 
1362  auto arcAngle = []( const QgsPoint & a, const QgsPoint & b, const QgsPoint & c )->double
1363  {
1364  double abX = b.x() - a.x();
1365  double abY = b.y() - a.y();
1366 
1367  double cbX = b.x() - c.x();
1368  double cbY = b.y() - c.y();
1369 
1370  double dot = ( abX * cbX + abY * cbY ); /* dot product */
1371  double cross = ( abX * cbY - abY * cbX ); /* cross product */
1372 
1373  double alpha = std::atan2( cross, dot );
1374 
1375  return alpha;
1376  };
1377 
1378  /* We make a candidate arc of the first two edges, */
1379  /* And then see if the next edge follows it */
1380  int i = 0;
1381  int j = 0;
1382  int k = 0;
1383  int currentArc = 1;
1384  QgsPoint a1;
1385  QgsPoint a2;
1386  QgsPoint a3;
1387  QgsPoint b;
1388  double centerX = 0.0;
1389  double centerY = 0.0;
1390  double radius = 0;
1391 
1392  while ( i < numEdges - 2 )
1393  {
1394  unsigned int arcEdges = 0;
1395  double numQuadrants = 0;
1396  double angle;
1397 
1398  bool foundArc = false;
1399  /* Make candidate arc */
1400  a1 = lineString->pointN( i );
1401  a2 = lineString->pointN( i + 1 );
1402  a3 = lineString->pointN( i + 2 );
1403  QgsPoint first = a1;
1404 
1405  for ( j = i + 3; j < numEdges + 1; j++ )
1406  {
1407  b = lineString->pointN( j );
1408 
1409  /* Does this point fall on our candidate arc? */
1410  if ( QgsGeometryUtils::pointContinuesArc( a1, a2, a3, b, distanceTolerance, pointSpacingAngleTolerance ) )
1411  {
1412  /* Yes. Mark this edge and the two preceding it as arc components */
1413  foundArc = true;
1414  for ( k = j - 1; k > j - 4; k-- )
1415  edgesInArcs[k] = currentArc;
1416  }
1417  else
1418  {
1419  /* No. So we're done with this candidate arc */
1420  currentArc++;
1421  break;
1422  }
1423 
1424  a1 = a2;
1425  a2 = a3;
1426  a3 = b;
1427  }
1428  /* Jump past all the edges that were added to the arc */
1429  if ( foundArc )
1430  {
1431  /* Check if an arc was composed by enough edges to be
1432  * really considered an arc
1433  * See http://trac.osgeo.org/postgis/ticket/2420
1434  */
1435  arcEdges = j - 1 - i;
1436  if ( first.x() == b.x() && first.y() == b.y() )
1437  {
1438  numQuadrants = 4;
1439  }
1440  else
1441  {
1442  QgsGeometryUtils::circleCenterRadius( first, b, a1, radius, centerX, centerY );
1443 
1444  angle = arcAngle( first, QgsPoint( centerX, centerY ), b );
1445  int p2Side = QgsGeometryUtils::leftOfLine( b.x(), b.y(), first.x(), first.y(), a1.x(), a1.y() );
1446  if ( p2Side >= 0 )
1447  angle = -angle;
1448 
1449  if ( angle < 0 )
1450  angle = 2 * M_PI + angle;
1451  numQuadrants = ( 4 * angle ) / ( 2 * M_PI );
1452  }
1453  /* a1 is first point, b is last point */
1454  if ( arcEdges < minQuadEdges * numQuadrants )
1455  {
1456  // LWDEBUGF( 4, "Not enough edges for a %g quadrants arc, %g needed", num_quadrants, min_quad_edges * num_quadrants );
1457  for ( k = j - 1; k >= i; k-- )
1458  edgesInArcs[k] = 0;
1459  }
1460 
1461  i = j - 1;
1462  }
1463  else
1464  {
1465  /* Mark this edge as a linear edge */
1466  edgesInArcs[i] = 0;
1467  i = i + 1;
1468  }
1469  }
1470 
1471  int start = 0;
1472  int end = 0;
1473  /* non-zero if edge is part of an arc */
1474  int edgeType = edgesInArcs[0];
1475 
1476  auto addPointsToCurve = [ lineString, &out ]( int start, int end, int type )
1477  {
1478  if ( type == 0 )
1479  {
1480  // straight segment
1481  QVector< QgsPoint > points;
1482  for ( int j = start; j < end + 2; ++ j )
1483  {
1484  points.append( lineString->pointN( j ) );
1485  }
1486  std::unique_ptr< QgsCurve > straightSegment = std::make_unique< QgsLineString >( points );
1487  out->addCurve( straightSegment.release() );
1488  }
1489  else
1490  {
1491  // curved segment
1492  QVector< QgsPoint > points;
1493  points.append( lineString->pointN( start ) );
1494  points.append( lineString->pointN( ( start + end + 1 ) / 2 ) );
1495  points.append( lineString->pointN( end + 1 ) );
1496  std::unique_ptr< QgsCircularString > curvedSegment = std::make_unique< QgsCircularString >();
1497  curvedSegment->setPoints( points );
1498  out->addCurve( curvedSegment.release() );
1499  }
1500  };
1501 
1502  for ( int i = 1; i < numEdges; i++ )
1503  {
1504  if ( edgeType != edgesInArcs[i] )
1505  {
1506  end = i - 1;
1507  addPointsToCurve( start, end, edgeType );
1508  start = i;
1509  edgeType = edgesInArcs[i];
1510  }
1511  }
1512 
1513  /* Roll out last item */
1514  end = numEdges - 1;
1515  addPointsToCurve( start, end, edgeType );
1516 
1517  return out;
1518 }
1519 
1520 std::unique_ptr< QgsAbstractGeometry > convertGeometryToCurves( const QgsAbstractGeometry *geom, double distanceTolerance, double angleTolerance )
1521 {
1523  {
1524  return lineToCurve( static_cast< const QgsLineString * >( geom ), distanceTolerance, angleTolerance );
1525  }
1526  else
1527  {
1528  // polygon
1529  const QgsPolygon *polygon = static_cast< const QgsPolygon * >( geom );
1530  std::unique_ptr< QgsCurvePolygon > result = std::make_unique< QgsCurvePolygon>();
1531 
1532  result->setExteriorRing( lineToCurve( static_cast< const QgsLineString * >( polygon->exteriorRing() ),
1533  distanceTolerance, angleTolerance ).release() );
1534  for ( int i = 0; i < polygon->numInteriorRings(); ++i )
1535  {
1536  result->addInteriorRing( lineToCurve( static_cast< const QgsLineString * >( polygon->interiorRing( i ) ),
1537  distanceTolerance, angleTolerance ).release() );
1538  }
1539 
1540  return result;
1541  }
1542 }
1543 
1544 QgsGeometry QgsInternalGeometryEngine::convertToCurves( double distanceTolerance, double angleTolerance ) const
1545 {
1546  mLastError.clear();
1547  if ( !mGeometry )
1548  {
1549  return QgsGeometry();
1550  }
1551 
1553  {
1554  return QgsGeometry( mGeometry->clone() ); // point geometry, nothing to do
1555  }
1556 
1557  if ( const QgsGeometryCollection *gc = qgsgeometry_cast< const QgsGeometryCollection *>( mGeometry ) )
1558  {
1559  int numGeom = gc->numGeometries();
1560  QVector< QgsAbstractGeometry * > geometryList;
1561  geometryList.reserve( numGeom );
1562  for ( int i = 0; i < numGeom; ++i )
1563  {
1564  geometryList << convertGeometryToCurves( gc->geometryN( i ), distanceTolerance, angleTolerance ).release();
1565  }
1566 
1567  QgsGeometry first = QgsGeometry( geometryList.takeAt( 0 ) );
1568  for ( QgsAbstractGeometry *g : std::as_const( geometryList ) )
1569  {
1570  first.addPart( g );
1571  }
1572  return first;
1573  }
1574  else
1575  {
1576  return QgsGeometry( convertGeometryToCurves( mGeometry, distanceTolerance, angleTolerance ) );
1577  }
1578 }
1579 
1580 QgsGeometry QgsInternalGeometryEngine::orientedMinimumBoundingBox( double &area, double &angle, double &width, double &height ) const
1581 {
1582  mLastError.clear();
1583 
1584  QgsRectangle minRect;
1585  area = std::numeric_limits<double>::max();
1586  angle = 0;
1587  width = std::numeric_limits<double>::max();
1588  height = std::numeric_limits<double>::max();
1589 
1590  if ( !mGeometry || mGeometry->nCoordinates() < 2 )
1591  return QgsGeometry();
1592 
1593  std::unique_ptr< QgsGeometryEngine >engine( QgsGeometry::createGeometryEngine( mGeometry ) );
1594  std::unique_ptr< QgsAbstractGeometry > hull( engine->convexHull( &mLastError ) );
1595  if ( !hull )
1596  return QgsGeometry();
1597 
1598  QgsVertexId vertexId;
1599  QgsPoint pt0;
1600  QgsPoint pt1;
1601  QgsPoint pt2;
1602  // get first point
1603  hull->nextVertex( vertexId, pt0 );
1604  pt1 = pt0;
1605  double totalRotation = 0;
1606  while ( hull->nextVertex( vertexId, pt2 ) )
1607  {
1608  double currentAngle = QgsGeometryUtils::lineAngle( pt1.x(), pt1.y(), pt2.x(), pt2.y() );
1609  double rotateAngle = 180.0 / M_PI * currentAngle;
1610  totalRotation += rotateAngle;
1611 
1612  QTransform t = QTransform::fromTranslate( pt0.x(), pt0.y() );
1613  t.rotate( rotateAngle );
1614  t.translate( -pt0.x(), -pt0.y() );
1615 
1616  hull->transform( t );
1617 
1618  QgsRectangle bounds = hull->boundingBox();
1619  double currentArea = bounds.width() * bounds.height();
1620  if ( currentArea < area )
1621  {
1622  minRect = bounds;
1623  area = currentArea;
1624  angle = totalRotation;
1625  width = bounds.width();
1626  height = bounds.height();
1627  }
1628 
1629  pt1 = hull->vertexAt( vertexId );
1630  }
1631 
1632  QgsGeometry minBounds = QgsGeometry::fromRect( minRect );
1633  minBounds.rotate( angle, QgsPointXY( pt0.x(), pt0.y() ) );
1634 
1635  if ( width > height )
1636  {
1637  width = minRect.height();
1638  height = minRect.width();
1639  angle = angle + 90.0;
1640  }
1641 
1642  // constrain angle to 0 - 180
1643  if ( angle > 180.0 )
1644  angle = std::fmod( angle, 180.0 );
1645 
1646  return minBounds;
1647 }
1648 
1649 std::unique_ptr< QgsLineString > triangularWavesAlongLine( const QgsLineString *line, double wavelength, const double amplitude, const bool strictWavelength )
1650 {
1651  const int totalPoints = line->numPoints();
1652  if ( totalPoints < 2 )
1653  return nullptr;
1654 
1655  const double *x = line->xData();
1656  const double *y = line->yData();
1657 
1658  double prevX = *x++;
1659  double prevY = *y++;
1660 
1661  QVector< double > outX;
1662  QVector< double > outY;
1663  const double totalLength = line->length();
1664 
1665  const int maxRepetitions = std::ceil( totalLength / wavelength );
1666  if ( !strictWavelength )
1667  wavelength = totalLength / maxRepetitions;
1668 
1669  const int estimatedPoints = maxRepetitions * 2 + 2;
1670  outX.reserve( estimatedPoints );
1671  outY.reserve( estimatedPoints );
1672  outX.append( prevX );
1673  outY.append( prevY );
1674 
1675  double distanceToNextPointFromStartOfSegment = wavelength / 4;
1676  int side = -1;
1677  for ( int i = 1; i < totalPoints; ++i )
1678  {
1679  double thisX = *x++;
1680  double thisY = *y++;
1681 
1682  const double segmentAngleRadians = QgsGeometryUtils::lineAngle( prevX, prevY, thisX, thisY );
1683  const double segmentLength = std::sqrt( ( thisX - prevX ) * ( thisX - prevX ) + ( thisY - prevY ) * ( thisY - prevY ) );
1684  while ( distanceToNextPointFromStartOfSegment < segmentLength || qgsDoubleNear( distanceToNextPointFromStartOfSegment, segmentLength ) )
1685  {
1686  // point falls on this segment - truncate to segment length if qgsDoubleNear test was actually > segment length
1687  const double distanceToPoint = std::min( distanceToNextPointFromStartOfSegment, segmentLength );
1688  double pX, pY;
1689  QgsGeometryUtils::pointOnLineWithDistance( prevX, prevY, thisX, thisY, distanceToPoint, pX, pY );
1690 
1691  // project point on line out by amplitude
1692  const double outPointX = pX + side * amplitude * std::sin( segmentAngleRadians + M_PI_2 );
1693  const double outPointY = pY + side * amplitude * std::cos( segmentAngleRadians + M_PI_2 );
1694 
1695  outX.append( outPointX );
1696  outY.append( outPointY );
1697 
1698  distanceToNextPointFromStartOfSegment += wavelength / 2;
1699  side = -side;
1700  }
1701 
1702  prevX = thisX;
1703  prevY = thisY;
1704  distanceToNextPointFromStartOfSegment -= segmentLength;
1705  }
1706 
1707  outX.append( prevX );
1708  outY.append( prevY );
1709 
1710  return std::make_unique< QgsLineString >( outX, outY );
1711 }
1712 
1713 std::unique_ptr< QgsLineString > triangularWavesRandomizedAlongLine( const QgsLineString *line,
1714  const double minimumWavelength, const double maximumWavelength,
1715  const double minimumAmplitude, const double maximumAmplitude,
1716  std::uniform_real_distribution<> &uniformDist, std::mt19937 &mt )
1717 {
1718  const int totalPoints = line->numPoints();
1719  if ( totalPoints < 2 )
1720  return nullptr;
1721 
1722  const double *x = line->xData();
1723  const double *y = line->yData();
1724 
1725  double prevX = *x++;
1726  double prevY = *y++;
1727 
1728  QVector< double > outX;
1729  QVector< double > outY;
1730  const double totalLength = line->length();
1731 
1732  const int maxRepetitions = std::ceil( totalLength / minimumWavelength );
1733 
1734  const int estimatedPoints = maxRepetitions * 2 + 2;
1735  outX.reserve( estimatedPoints );
1736  outY.reserve( estimatedPoints );
1737  outX.append( prevX );
1738  outY.append( prevY );
1739 
1740  double wavelength = uniformDist( mt ) * ( maximumWavelength - minimumWavelength ) + minimumWavelength;
1741  double distanceToNextPointFromStartOfSegment = wavelength / 4;
1742  double amplitude = uniformDist( mt ) * ( maximumAmplitude - minimumAmplitude ) + minimumAmplitude;
1743 
1744  int side = -1;
1745  for ( int i = 1; i < totalPoints; ++i )
1746  {
1747  double thisX = *x++;
1748  double thisY = *y++;
1749 
1750  const double segmentAngleRadians = QgsGeometryUtils::lineAngle( prevX, prevY, thisX, thisY );
1751  const double segmentLength = std::sqrt( ( thisX - prevX ) * ( thisX - prevX ) + ( thisY - prevY ) * ( thisY - prevY ) );
1752  while ( distanceToNextPointFromStartOfSegment < segmentLength || qgsDoubleNear( distanceToNextPointFromStartOfSegment, segmentLength ) )
1753  {
1754  // point falls on this segment - truncate to segment length if qgsDoubleNear test was actually > segment length
1755  const double distanceToPoint = std::min( distanceToNextPointFromStartOfSegment, segmentLength );
1756  double pX, pY;
1757  QgsGeometryUtils::pointOnLineWithDistance( prevX, prevY, thisX, thisY, distanceToPoint, pX, pY );
1758 
1759  // project point on line out by amplitude
1760  const double outPointX = pX + side * amplitude * std::sin( segmentAngleRadians + M_PI_2 );
1761  const double outPointY = pY + side * amplitude * std::cos( segmentAngleRadians + M_PI_2 );
1762  outX.append( outPointX );
1763  outY.append( outPointY );
1764 
1765  wavelength = uniformDist( mt ) * ( maximumWavelength - minimumWavelength ) + minimumWavelength;
1766  amplitude = uniformDist( mt ) * ( maximumAmplitude - minimumAmplitude ) + minimumAmplitude;
1767 
1768  distanceToNextPointFromStartOfSegment += wavelength / 2;
1769  side = -side;
1770  }
1771 
1772  prevX = thisX;
1773  prevY = thisY;
1774  distanceToNextPointFromStartOfSegment -= segmentLength;
1775  }
1776 
1777  outX.append( prevX );
1778  outY.append( prevY );
1779 
1780  return std::make_unique< QgsLineString >( outX, outY );
1781 }
1782 
1783 std::unique_ptr< QgsAbstractGeometry > triangularWavesPrivate( const QgsAbstractGeometry *geom, double wavelength, double amplitude, bool strictWavelength )
1784 {
1785  std::unique_ptr< QgsAbstractGeometry > segmentizedCopy;
1786  if ( QgsWkbTypes::isCurvedType( geom->wkbType() ) )
1787  {
1788  segmentizedCopy.reset( geom->segmentize() );
1789  geom = segmentizedCopy.get();
1790  }
1791 
1793  {
1794  return triangularWavesAlongLine( static_cast< const QgsLineString * >( geom ), wavelength, amplitude, strictWavelength );
1795  }
1796  else
1797  {
1798  // polygon
1799  const QgsPolygon *polygon = static_cast< const QgsPolygon * >( geom );
1800  std::unique_ptr< QgsPolygon > result = std::make_unique< QgsPolygon >();
1801 
1802  result->setExteriorRing( triangularWavesAlongLine( static_cast< const QgsLineString * >( polygon->exteriorRing() ),
1803  wavelength, amplitude, strictWavelength ).release() );
1804  for ( int i = 0; i < polygon->numInteriorRings(); ++i )
1805  {
1806  result->addInteriorRing( triangularWavesAlongLine( static_cast< const QgsLineString * >( polygon->interiorRing( i ) ),
1807  wavelength, amplitude, strictWavelength ).release() );
1808  }
1809 
1810  return result;
1811  }
1812 }
1813 
1814 std::unique_ptr< QgsAbstractGeometry > triangularWavesRandomizedPrivate( const QgsAbstractGeometry *geom, double minimumWavelength, double maximumWavelength, double minimumAmplitude, double maximumAmplitude, std::uniform_real_distribution<> &uniformDist, std::mt19937 &mt )
1815 {
1816  std::unique_ptr< QgsAbstractGeometry > segmentizedCopy;
1817  if ( QgsWkbTypes::isCurvedType( geom->wkbType() ) )
1818  {
1819  segmentizedCopy.reset( geom->segmentize() );
1820  geom = segmentizedCopy.get();
1821  }
1822 
1824  {
1825  return triangularWavesRandomizedAlongLine( static_cast< const QgsLineString * >( geom ), minimumWavelength, maximumWavelength, minimumAmplitude, maximumAmplitude, uniformDist, mt );
1826  }
1827  else
1828  {
1829  // polygon
1830  const QgsPolygon *polygon = static_cast< const QgsPolygon * >( geom );
1831  std::unique_ptr< QgsPolygon > result = std::make_unique< QgsPolygon >();
1832 
1833  result->setExteriorRing( triangularWavesRandomizedAlongLine( static_cast< const QgsLineString * >( polygon->exteriorRing() ),
1834  minimumWavelength, maximumWavelength, minimumAmplitude, maximumAmplitude, uniformDist, mt ).release() );
1835  for ( int i = 0; i < polygon->numInteriorRings(); ++i )
1836  {
1837  result->addInteriorRing( triangularWavesRandomizedAlongLine( static_cast< const QgsLineString * >( polygon->interiorRing( i ) ),
1838  minimumWavelength, maximumWavelength, minimumAmplitude, maximumAmplitude, uniformDist, mt ).release() );
1839  }
1840 
1841  return result;
1842  }
1843 }
1844 
1845 QgsGeometry QgsInternalGeometryEngine::triangularWaves( double wavelength, double amplitude, bool strictWavelength ) const
1846 {
1847  if ( wavelength < 0 || qgsDoubleNear( wavelength, 0 ) )
1848  return QgsGeometry();
1849 
1850  mLastError.clear();
1851  if ( !mGeometry )
1852  {
1853  return QgsGeometry();
1854  }
1855 
1857  {
1858  return QgsGeometry( mGeometry->clone() ); // point geometry, nothing to do
1859  }
1860 
1861  if ( const QgsGeometryCollection *gc = qgsgeometry_cast< const QgsGeometryCollection *>( mGeometry ) )
1862  {
1863  int numGeom = gc->numGeometries();
1864  QVector< QgsAbstractGeometry * > geometryList;
1865  geometryList.reserve( numGeom );
1866  for ( int i = 0; i < numGeom; ++i )
1867  {
1868  geometryList << triangularWavesPrivate( gc->geometryN( i ), wavelength, amplitude, strictWavelength ).release();
1869  }
1870 
1871  QgsGeometry first = QgsGeometry( geometryList.takeAt( 0 ) );
1872  for ( QgsAbstractGeometry *g : std::as_const( geometryList ) )
1873  {
1874  first.addPart( g );
1875  }
1876  return first;
1877  }
1878  else
1879  {
1880  return QgsGeometry( triangularWavesPrivate( mGeometry, wavelength, amplitude, strictWavelength ) );
1881  }
1882 }
1883 
1884 QgsGeometry QgsInternalGeometryEngine::triangularWavesRandomized( double minimumWavelength, double maximumWavelength, double minimumAmplitude, double maximumAmplitude, unsigned long seed ) const
1885 {
1886  if ( minimumWavelength < 0 || qgsDoubleNear( minimumWavelength, 0 ) || maximumWavelength < 0 || qgsDoubleNear( maximumWavelength, 0 ) || maximumWavelength < minimumWavelength )
1887  return QgsGeometry();
1888 
1889  mLastError.clear();
1890  if ( !mGeometry )
1891  {
1892  return QgsGeometry();
1893  }
1894 
1895  std::random_device rd;
1896  std::mt19937 mt( seed == 0 ? rd() : seed );
1897  std::uniform_real_distribution<> uniformDist( 0, 1 );
1898 
1900  {
1901  return QgsGeometry( mGeometry->clone() ); // point geometry, nothing to do
1902  }
1903 
1904  if ( const QgsGeometryCollection *gc = qgsgeometry_cast< const QgsGeometryCollection *>( mGeometry ) )
1905  {
1906  int numGeom = gc->numGeometries();
1907  QVector< QgsAbstractGeometry * > geometryList;
1908  geometryList.reserve( numGeom );
1909  for ( int i = 0; i < numGeom; ++i )
1910  {
1911  geometryList << triangularWavesRandomizedPrivate( gc->geometryN( i ), minimumWavelength, maximumWavelength, minimumAmplitude, maximumAmplitude, uniformDist, mt ).release();
1912  }
1913 
1914  QgsGeometry first = QgsGeometry( geometryList.takeAt( 0 ) );
1915  for ( QgsAbstractGeometry *g : std::as_const( geometryList ) )
1916  {
1917  first.addPart( g );
1918  }
1919  return first;
1920  }
1921  else
1922  {
1923  return QgsGeometry( triangularWavesRandomizedPrivate( mGeometry, minimumWavelength, maximumWavelength, minimumAmplitude, maximumAmplitude, uniformDist, mt ) );
1924  }
1925 }
1926 
1927 std::unique_ptr< QgsLineString > squareWavesAlongLine( const QgsLineString *line, double wavelength, const double amplitude, const bool strictWavelength )
1928 {
1929  const int totalPoints = line->numPoints();
1930  if ( totalPoints < 2 )
1931  return nullptr;
1932 
1933  const double *x = line->xData();
1934  const double *y = line->yData();
1935 
1936  double prevX = *x++;
1937  double prevY = *y++;
1938 
1939  QVector< double > outX;
1940  QVector< double > outY;
1941  const double totalLength = line->length();
1942 
1943  const int maxRepetitions = std::ceil( totalLength / wavelength );
1944  if ( !strictWavelength )
1945  wavelength = totalLength / maxRepetitions;
1946 
1947  const int estimatedPoints = maxRepetitions * 4 + 2;
1948  outX.reserve( estimatedPoints );
1949  outY.reserve( estimatedPoints );
1950  outX.append( prevX );
1951  outY.append( prevY );
1952 
1953  const double startSegmentAngleRadians = QgsGeometryUtils::lineAngle( prevX, prevY, *x, *y );
1954  outX.append( prevX - amplitude * std::sin( startSegmentAngleRadians + M_PI_2 ) );
1955  outY.append( prevY - amplitude * std::cos( startSegmentAngleRadians + M_PI_2 ) );
1956 
1957  double distanceToNextPointFromStartOfSegment = wavelength / 2;
1958 
1959  int side = -1;
1960  for ( int i = 1; i < totalPoints; ++i )
1961  {
1962  double thisX = *x++;
1963  double thisY = *y++;
1964 
1965  const double segmentAngleRadians = QgsGeometryUtils::lineAngle( prevX, prevY, thisX, thisY );
1966  const double segmentLength = std::sqrt( ( thisX - prevX ) * ( thisX - prevX ) + ( thisY - prevY ) * ( thisY - prevY ) );
1967  while ( distanceToNextPointFromStartOfSegment < segmentLength || qgsDoubleNear( distanceToNextPointFromStartOfSegment, segmentLength ) )
1968  {
1969  // point falls on this segment - truncate to segment length if qgsDoubleNear test was actually > segment length
1970  const double distanceToPoint = std::min( distanceToNextPointFromStartOfSegment, segmentLength );
1971  double pX, pY;
1972  QgsGeometryUtils::pointOnLineWithDistance( prevX, prevY, thisX, thisY, distanceToPoint, pX, pY );
1973 
1974  // project point on line out by amplitude
1975  const double sinAngle = std::sin( segmentAngleRadians + M_PI_2 );
1976  const double cosAngle = std::cos( segmentAngleRadians + M_PI_2 );
1977  outX.append( pX + side * amplitude * sinAngle );
1978  outY.append( pY + side * amplitude * cosAngle );
1979  outX.append( pX - side * amplitude * sinAngle );
1980  outY.append( pY - side * amplitude * cosAngle );
1981 
1982  distanceToNextPointFromStartOfSegment += wavelength / 2;
1983  side = -side;
1984  }
1985 
1986  prevX = thisX;
1987  prevY = thisY;
1988  distanceToNextPointFromStartOfSegment -= segmentLength;
1989  }
1990 
1991  // replace last point, which will be projected out to amplitude of wave, with the actual end point of the line
1992  outX.pop_back();
1993  outY.pop_back();
1994  outX.append( prevX );
1995  outY.append( prevY );
1996 
1997  return std::make_unique< QgsLineString >( outX, outY );
1998 }
1999 
2000 std::unique_ptr< QgsLineString > squareWavesRandomizedAlongLine( const QgsLineString *line,
2001  const double minimumWavelength, const double maximumWavelength,
2002  const double minimumAmplitude, const double maximumAmplitude,
2003  std::uniform_real_distribution<> &uniformDist, std::mt19937 &mt )
2004 {
2005  const int totalPoints = line->numPoints();
2006  if ( totalPoints < 2 )
2007  return nullptr;
2008 
2009  const double *x = line->xData();
2010  const double *y = line->yData();
2011 
2012  double prevX = *x++;
2013  double prevY = *y++;
2014 
2015  QVector< double > outX;
2016  QVector< double > outY;
2017  const double totalLength = line->length();
2018 
2019  const int maxRepetitions = std::ceil( totalLength / minimumWavelength );
2020 
2021  const int estimatedPoints = maxRepetitions * 4 + 2;
2022  outX.reserve( estimatedPoints );
2023  outY.reserve( estimatedPoints );
2024  outX.append( prevX );
2025  outY.append( prevY );
2026 
2027  double amplitude = uniformDist( mt ) * ( maximumAmplitude - minimumAmplitude ) + minimumAmplitude;
2028 
2029  double segmentAngleRadians = QgsGeometryUtils::lineAngle( prevX, prevY, *x, *y );
2030  outX.append( prevX - amplitude * std::sin( segmentAngleRadians + M_PI_2 ) );
2031  outY.append( prevY - amplitude * std::cos( segmentAngleRadians + M_PI_2 ) );
2032 
2033  double wavelength = uniformDist( mt ) * ( maximumWavelength - minimumWavelength ) + minimumWavelength;
2034  double distanceToNextPointFromStartOfSegment = wavelength / 2;
2035 
2036  int side = -1;
2037  for ( int i = 1; i < totalPoints; ++i )
2038  {
2039  double thisX = *x++;
2040  double thisY = *y++;
2041 
2042  segmentAngleRadians = QgsGeometryUtils::lineAngle( prevX, prevY, thisX, thisY );
2043  const double segmentLength = std::sqrt( ( thisX - prevX ) * ( thisX - prevX ) + ( thisY - prevY ) * ( thisY - prevY ) );
2044  while ( distanceToNextPointFromStartOfSegment < segmentLength || qgsDoubleNear( distanceToNextPointFromStartOfSegment, segmentLength ) )
2045  {
2046  // point falls on this segment - truncate to segment length if qgsDoubleNear test was actually > segment length
2047  const double distanceToPoint = std::min( distanceToNextPointFromStartOfSegment, segmentLength );
2048  double pX, pY;
2049  QgsGeometryUtils::pointOnLineWithDistance( prevX, prevY, thisX, thisY, distanceToPoint, pX, pY );
2050 
2051  // project point on line out by amplitude
2052  const double sinAngle = std::sin( segmentAngleRadians + M_PI_2 );
2053  const double cosAngle = std::cos( segmentAngleRadians + M_PI_2 );
2054  outX.append( pX + side * amplitude * sinAngle );
2055  outY.append( pY + side * amplitude * cosAngle );
2056 
2057  amplitude = uniformDist( mt ) * ( maximumAmplitude - minimumAmplitude ) + minimumAmplitude;
2058  outX.append( pX - side * amplitude * sinAngle );
2059  outY.append( pY - side * amplitude * cosAngle );
2060 
2061  wavelength = uniformDist( mt ) * ( maximumWavelength - minimumWavelength ) + minimumWavelength;
2062  distanceToNextPointFromStartOfSegment += wavelength / 2;
2063  side = -side;
2064  }
2065 
2066  prevX = thisX;
2067  prevY = thisY;
2068  distanceToNextPointFromStartOfSegment -= segmentLength;
2069  }
2070 
2071  outX.append( prevX + side * amplitude * std::sin( segmentAngleRadians + M_PI_2 ) );
2072  outY.append( prevY + side * amplitude * std::cos( segmentAngleRadians + M_PI_2 ) );
2073  outX.append( prevX );
2074  outY.append( prevY );
2075 
2076  return std::make_unique< QgsLineString >( outX, outY );
2077 }
2078 
2079 std::unique_ptr< QgsAbstractGeometry > squareWavesPrivate( const QgsAbstractGeometry *geom, double wavelength, double amplitude, bool strictWavelength )
2080 {
2081  std::unique_ptr< QgsAbstractGeometry > segmentizedCopy;
2082  if ( QgsWkbTypes::isCurvedType( geom->wkbType() ) )
2083  {
2084  segmentizedCopy.reset( geom->segmentize() );
2085  geom = segmentizedCopy.get();
2086  }
2087 
2089  {
2090  return squareWavesAlongLine( static_cast< const QgsLineString * >( geom ), wavelength, amplitude, strictWavelength );
2091  }
2092  else
2093  {
2094  // polygon
2095  const QgsPolygon *polygon = static_cast< const QgsPolygon * >( geom );
2096  std::unique_ptr< QgsPolygon > result = std::make_unique< QgsPolygon >();
2097 
2098  result->setExteriorRing( squareWavesAlongLine( static_cast< const QgsLineString * >( polygon->exteriorRing() ),
2099  wavelength, amplitude, strictWavelength ).release() );
2100  for ( int i = 0; i < polygon->numInteriorRings(); ++i )
2101  {
2102  result->addInteriorRing( squareWavesAlongLine( static_cast< const QgsLineString * >( polygon->interiorRing( i ) ),
2103  wavelength, amplitude, strictWavelength ).release() );
2104  }
2105 
2106  return result;
2107  }
2108 }
2109 
2110 std::unique_ptr< QgsAbstractGeometry > squareWavesRandomizedPrivate( const QgsAbstractGeometry *geom, double minimumWavelength, double maximumWavelength, double minimumAmplitude, double maximumAmplitude, std::uniform_real_distribution<> &uniformDist, std::mt19937 &mt )
2111 {
2112  std::unique_ptr< QgsAbstractGeometry > segmentizedCopy;
2113  if ( QgsWkbTypes::isCurvedType( geom->wkbType() ) )
2114  {
2115  segmentizedCopy.reset( geom->segmentize() );
2116  geom = segmentizedCopy.get();
2117  }
2118 
2120  {
2121  return squareWavesRandomizedAlongLine( static_cast< const QgsLineString * >( geom ), minimumWavelength, maximumWavelength, minimumAmplitude, maximumAmplitude, uniformDist, mt );
2122  }
2123  else
2124  {
2125  // polygon
2126  const QgsPolygon *polygon = static_cast< const QgsPolygon * >( geom );
2127  std::unique_ptr< QgsPolygon > result = std::make_unique< QgsPolygon >();
2128 
2129  result->setExteriorRing( squareWavesRandomizedAlongLine( static_cast< const QgsLineString * >( polygon->exteriorRing() ),
2130  minimumWavelength, maximumWavelength, minimumAmplitude, maximumAmplitude, uniformDist, mt ).release() );
2131  for ( int i = 0; i < polygon->numInteriorRings(); ++i )
2132  {
2133  result->addInteriorRing( squareWavesRandomizedAlongLine( static_cast< const QgsLineString * >( polygon->interiorRing( i ) ),
2134  minimumWavelength, maximumWavelength, minimumAmplitude, maximumAmplitude, uniformDist, mt ).release() );
2135  }
2136 
2137  return result;
2138  }
2139 }
2140 
2141 QgsGeometry QgsInternalGeometryEngine::squareWaves( double wavelength, double amplitude, bool strictWavelength ) const
2142 {
2143  if ( wavelength < 0 || qgsDoubleNear( wavelength, 0 ) )
2144  return QgsGeometry();
2145 
2146  mLastError.clear();
2147  if ( !mGeometry )
2148  {
2149  return QgsGeometry();
2150  }
2151 
2153  {
2154  return QgsGeometry( mGeometry->clone() ); // point geometry, nothing to do
2155  }
2156 
2157  if ( const QgsGeometryCollection *gc = qgsgeometry_cast< const QgsGeometryCollection *>( mGeometry ) )
2158  {
2159  int numGeom = gc->numGeometries();
2160  QVector< QgsAbstractGeometry * > geometryList;
2161  geometryList.reserve( numGeom );
2162  for ( int i = 0; i < numGeom; ++i )
2163  {
2164  geometryList << squareWavesPrivate( gc->geometryN( i ), wavelength, amplitude, strictWavelength ).release();
2165  }
2166 
2167  QgsGeometry first = QgsGeometry( geometryList.takeAt( 0 ) );
2168  for ( QgsAbstractGeometry *g : std::as_const( geometryList ) )
2169  {
2170  first.addPart( g );
2171  }
2172  return first;
2173  }
2174  else
2175  {
2176  return QgsGeometry( squareWavesPrivate( mGeometry, wavelength, amplitude, strictWavelength ) );
2177  }
2178 }
2179 
2180 QgsGeometry QgsInternalGeometryEngine::squareWavesRandomized( double minimumWavelength, double maximumWavelength, double minimumAmplitude, double maximumAmplitude, unsigned long seed ) const
2181 {
2182  if ( minimumWavelength < 0 || qgsDoubleNear( minimumWavelength, 0 ) || maximumWavelength < 0 || qgsDoubleNear( maximumWavelength, 0 ) || maximumWavelength < minimumWavelength )
2183  return QgsGeometry();
2184 
2185  mLastError.clear();
2186  if ( !mGeometry )
2187  {
2188  return QgsGeometry();
2189  }
2190 
2191  std::random_device rd;
2192  std::mt19937 mt( seed == 0 ? rd() : seed );
2193  std::uniform_real_distribution<> uniformDist( 0, 1 );
2194 
2196  {
2197  return QgsGeometry( mGeometry->clone() ); // point geometry, nothing to do
2198  }
2199 
2200  if ( const QgsGeometryCollection *gc = qgsgeometry_cast< const QgsGeometryCollection *>( mGeometry ) )
2201  {
2202  int numGeom = gc->numGeometries();
2203  QVector< QgsAbstractGeometry * > geometryList;
2204  geometryList.reserve( numGeom );
2205  for ( int i = 0; i < numGeom; ++i )
2206  {
2207  geometryList << squareWavesRandomizedPrivate( gc->geometryN( i ), minimumWavelength, maximumWavelength, minimumAmplitude, maximumAmplitude, uniformDist, mt ).release();
2208  }
2209 
2210  QgsGeometry first = QgsGeometry( geometryList.takeAt( 0 ) );
2211  for ( QgsAbstractGeometry *g : std::as_const( geometryList ) )
2212  {
2213  first.addPart( g );
2214  }
2215  return first;
2216  }
2217  else
2218  {
2219  return QgsGeometry( squareWavesRandomizedPrivate( mGeometry, minimumWavelength, maximumWavelength, minimumAmplitude, maximumAmplitude, uniformDist, mt ) );
2220  }
2221 }
2222 
2223 std::unique_ptr< QgsLineString > roundWavesAlongLine( const QgsLineString *line, double wavelength, const double amplitude, const bool strictWavelength )
2224 {
2225  const int totalPoints = line->numPoints();
2226  if ( totalPoints < 2 )
2227  return nullptr;
2228 
2229  const double *x = line->xData();
2230  const double *y = line->yData();
2231 
2232  double prevX = *x++;
2233  double prevY = *y++;
2234 
2235  const double totalLength = line->length();
2236 
2237  const int maxRepetitions = std::ceil( totalLength / wavelength );
2238  if ( !strictWavelength )
2239  wavelength = totalLength / maxRepetitions;
2240 
2241  const int segments = 10;
2242 
2243  int side = -1;
2244 
2245  double xOutBuffer[4] { prevX, prevX, prevX, prevX };
2246  double yOutBuffer[4] { prevY, prevY, prevY, prevY };
2247  bool isFirstPart = true;
2248 
2249  double distanceToNextPointFromStartOfSegment = wavelength / 8;
2250  int bufferIndex = 1;
2251  std::unique_ptr< QgsLineString > out = std::make_unique< QgsLineString >();
2252 
2253  double segmentAngleRadians = 0;
2254  double remainingDistance = totalLength;
2255  double totalCoveredDistance = 0;
2256 
2257  for ( int i = 1; i < totalPoints; ++i )
2258  {
2259  double thisX = *x++;
2260  double thisY = *y++;
2261 
2262  segmentAngleRadians = QgsGeometryUtils::lineAngle( prevX, prevY, thisX, thisY );
2263  const double segmentLength = std::sqrt( ( thisX - prevX ) * ( thisX - prevX ) + ( thisY - prevY ) * ( thisY - prevY ) );
2264  while ( distanceToNextPointFromStartOfSegment < segmentLength || qgsDoubleNear( distanceToNextPointFromStartOfSegment, segmentLength ) )
2265  {
2266  // point falls on this segment - truncate to segment length if qgsDoubleNear test was actually > segment length
2267  const double distanceToPoint = std::min( distanceToNextPointFromStartOfSegment, segmentLength );
2268  double pX, pY;
2269  QgsGeometryUtils::pointOnLineWithDistance( prevX, prevY, thisX, thisY, distanceToPoint, pX, pY );
2270  remainingDistance = totalLength - totalCoveredDistance - distanceToPoint;
2271 
2272  const double sinAngle = std::sin( segmentAngleRadians + M_PI_2 );
2273  const double cosAngle = std::cos( segmentAngleRadians + M_PI_2 );
2274 
2275  if ( bufferIndex == 0 )
2276  {
2277  xOutBuffer[0] = pX + side * amplitude * sinAngle;
2278  yOutBuffer[0] = pY + side * amplitude * cosAngle;
2279  bufferIndex = 1;
2280  distanceToNextPointFromStartOfSegment += wavelength / 4;
2281  }
2282  else if ( bufferIndex == 1 && isFirstPart )
2283  {
2284  xOutBuffer[1] = ( xOutBuffer[0] + pX - side * amplitude * sinAngle ) * 0.5;
2285  yOutBuffer[1] = ( yOutBuffer[0] + pY - side * amplitude * cosAngle ) * 0.5;
2286  xOutBuffer[2] = pX - side * amplitude * sinAngle;
2287  yOutBuffer[2] = pY - side * amplitude * cosAngle;
2288  bufferIndex = 2;
2289  distanceToNextPointFromStartOfSegment += wavelength / 8;
2290  }
2291  else if ( bufferIndex == 1 )
2292  {
2293  xOutBuffer[1] = pX + side * amplitude * sinAngle;
2294  yOutBuffer[1] = pY + side * amplitude * cosAngle;
2295  xOutBuffer[2] = pX - side * amplitude * sinAngle;
2296  yOutBuffer[2] = pY - side * amplitude * cosAngle;
2297  bufferIndex = 2;
2298  distanceToNextPointFromStartOfSegment += wavelength / 4;
2299  }
2300  else if ( bufferIndex == 2 )
2301  {
2302  xOutBuffer[3] = pX - side * amplitude * sinAngle;
2303  yOutBuffer[3] = pY - side * amplitude * cosAngle;
2304 
2305  if ( isFirstPart )
2306  {
2307  xOutBuffer[2] = ( xOutBuffer[2] + xOutBuffer[3] ) * 0.5;
2308  yOutBuffer[2] = ( yOutBuffer[2] + yOutBuffer[3] ) * 0.5;
2309  isFirstPart = false;
2310  }
2311 
2312  // flush curve
2313  std::unique_ptr< QgsLineString > bezier( QgsLineString::fromBezierCurve( QgsPoint( xOutBuffer[0], yOutBuffer[0] ),
2314  QgsPoint( xOutBuffer[1], yOutBuffer[1] ),
2315  QgsPoint( xOutBuffer[2], yOutBuffer[2] ),
2316  QgsPoint( xOutBuffer[3], yOutBuffer[3] ),
2317  segments ) );
2318  out->append( bezier.get() );
2319 
2320  // shuffle buffer alone
2321  xOutBuffer[0] = xOutBuffer[3];
2322  yOutBuffer[0] = yOutBuffer[3];
2323  bufferIndex = 1;
2324  side = -side;
2325  distanceToNextPointFromStartOfSegment += wavelength / 4;
2326  }
2327  }
2328  totalCoveredDistance += segmentLength;
2329  prevX = thisX;
2330  prevY = thisY;
2331  distanceToNextPointFromStartOfSegment -= segmentLength;
2332  }
2333 
2334  const double midX = prevX - remainingDistance / 2 * std::sin( segmentAngleRadians );
2335  const double midY = prevY - remainingDistance / 2 * std::cos( segmentAngleRadians );
2336  const double pX = midX + side * amplitude * std::sin( segmentAngleRadians + M_PI_2 );
2337  const double pY = midY + side * amplitude * std::cos( segmentAngleRadians + M_PI_2 );
2338 
2339  std::unique_ptr< QgsLineString > bezier( QgsLineString::fromBezierCurve( out->endPoint(),
2340  QgsPoint( ( out->endPoint().x() + pX ) * 0.5, ( out->endPoint().y() + pY ) * 0.5 ),
2341  QgsPoint( ( prevX + pX ) * 0.5, ( prevY + pY ) * 0.5 ),
2342  QgsPoint( prevX, prevY ),
2343  segments / 2 ) );
2344  out->append( bezier.get() );
2345 
2346  return out;
2347 }
2348 
2349 std::unique_ptr< QgsLineString > roundWavesRandomizedAlongLine( const QgsLineString *line,
2350  const double minimumWavelength, const double maximumWavelength,
2351  const double minimumAmplitude, const double maximumAmplitude,
2352  std::uniform_real_distribution<> &uniformDist, std::mt19937 &mt )
2353 {
2354  const int totalPoints = line->numPoints();
2355  if ( totalPoints < 2 )
2356  return nullptr;
2357 
2358  const double *x = line->xData();
2359  const double *y = line->yData();
2360 
2361  double prevX = *x++;
2362  double prevY = *y++;
2363 
2364  const double totalLength = line->length();
2365 
2366  const int segments = 10;
2367 
2368  int side = -1;
2369 
2370  double xOutBuffer[4] { prevX, prevX, prevX, prevX };
2371  double yOutBuffer[4] { prevY, prevY, prevY, prevY };
2372  bool isFirstPart = true;
2373 
2374  double amplitude = uniformDist( mt ) * ( maximumAmplitude - minimumAmplitude ) + minimumAmplitude;
2375  double wavelength = uniformDist( mt ) * ( maximumWavelength - minimumWavelength ) + minimumWavelength;
2376 
2377  double distanceToNextPointFromStartOfSegment = wavelength / 8;
2378  int bufferIndex = 1;
2379  std::unique_ptr< QgsLineString > out = std::make_unique< QgsLineString >();
2380 
2381  double segmentAngleRadians = 0;
2382 
2383  double remainingDistance = totalLength;
2384  double totalCoveredDistance = 0;
2385 
2386  for ( int i = 1; i < totalPoints; ++i )
2387  {
2388  double thisX = *x++;
2389  double thisY = *y++;
2390 
2391  segmentAngleRadians = QgsGeometryUtils::lineAngle( prevX, prevY, thisX, thisY );
2392  const double segmentLength = std::sqrt( ( thisX - prevX ) * ( thisX - prevX ) + ( thisY - prevY ) * ( thisY - prevY ) );
2393  while ( distanceToNextPointFromStartOfSegment < segmentLength || qgsDoubleNear( distanceToNextPointFromStartOfSegment, segmentLength ) )
2394  {
2395  // point falls on this segment - truncate to segment length if qgsDoubleNear test was actually > segment length
2396  const double distanceToPoint = std::min( distanceToNextPointFromStartOfSegment, segmentLength );
2397  double pX, pY;
2398  QgsGeometryUtils::pointOnLineWithDistance( prevX, prevY, thisX, thisY, distanceToPoint, pX, pY );
2399  remainingDistance = totalLength - totalCoveredDistance - distanceToPoint;
2400 
2401  const double sinAngle = std::sin( segmentAngleRadians + M_PI_2 );
2402  const double cosAngle = std::cos( segmentAngleRadians + M_PI_2 );
2403 
2404  if ( bufferIndex == 0 )
2405  {
2406  xOutBuffer[0] = pX + side * amplitude * sinAngle;
2407  yOutBuffer[0] = pY + side * amplitude * cosAngle;
2408  bufferIndex = 1;
2409  distanceToNextPointFromStartOfSegment += wavelength / 4;
2410  }
2411  else if ( bufferIndex == 1 && isFirstPart )
2412  {
2413  xOutBuffer[1] = ( xOutBuffer[0] + pX - side * amplitude * sinAngle ) * 0.5;
2414  yOutBuffer[1] = ( yOutBuffer[0] + pY - side * amplitude * cosAngle ) * 0.5;
2415  xOutBuffer[2] = pX - side * amplitude * sinAngle;
2416  yOutBuffer[2] = pY - side * amplitude * cosAngle;
2417  bufferIndex = 2;
2418  distanceToNextPointFromStartOfSegment += wavelength / 8;
2419  }
2420  else if ( bufferIndex == 1 )
2421  {
2422  xOutBuffer[1] = pX + side * amplitude * sinAngle;
2423  yOutBuffer[1] = pY + side * amplitude * cosAngle;
2424  amplitude = uniformDist( mt ) * ( maximumAmplitude - minimumAmplitude ) + minimumAmplitude;
2425  xOutBuffer[2] = pX - side * amplitude * sinAngle;
2426  bufferIndex = 2;
2427  yOutBuffer[2] = pY - side * amplitude * cosAngle;
2428  distanceToNextPointFromStartOfSegment += wavelength / 4;
2429  }
2430  else if ( bufferIndex == 2 )
2431  {
2432  xOutBuffer[3] = pX - side * amplitude * sinAngle;
2433  yOutBuffer[3] = pY - side * amplitude * cosAngle;
2434 
2435  if ( isFirstPart )
2436  {
2437  xOutBuffer[2] = ( xOutBuffer[2] + xOutBuffer[3] ) * 0.5;
2438  yOutBuffer[2] = ( yOutBuffer[2] + yOutBuffer[3] ) * 0.5;
2439  isFirstPart = false;
2440  }
2441 
2442  // flush curve
2443  std::unique_ptr< QgsLineString > bezier( QgsLineString::fromBezierCurve( QgsPoint( xOutBuffer[0], yOutBuffer[0] ),
2444  QgsPoint( xOutBuffer[1], yOutBuffer[1] ),
2445  QgsPoint( xOutBuffer[2], yOutBuffer[2] ),
2446  QgsPoint( xOutBuffer[3], yOutBuffer[3] ),
2447  segments ) );
2448  out->append( bezier.get() );
2449 
2450  // shuffle buffer alone
2451  xOutBuffer[0] = xOutBuffer[3];
2452  yOutBuffer[0] = yOutBuffer[3];
2453  bufferIndex = 1;
2454  side = -side;
2455 
2456  wavelength = uniformDist( mt ) * ( maximumWavelength - minimumWavelength ) + minimumWavelength;
2457 
2458  distanceToNextPointFromStartOfSegment += wavelength / 4;
2459  }
2460  }
2461  totalCoveredDistance += segmentLength;
2462 
2463  prevX = thisX;
2464  prevY = thisY;
2465  distanceToNextPointFromStartOfSegment -= segmentLength;
2466  }
2467 
2468  const double midX = prevX - remainingDistance / 2 * std::sin( segmentAngleRadians );
2469  const double midY = prevY - remainingDistance / 2 * std::cos( segmentAngleRadians );
2470  const double pX = midX + side * amplitude * std::sin( segmentAngleRadians + M_PI_2 );
2471  const double pY = midY + side * amplitude * std::cos( segmentAngleRadians + M_PI_2 );
2472 
2473  if ( out->isEmpty() )
2474  {
2475  std::unique_ptr< QgsLineString > bezier( QgsLineString::fromBezierCurve( line->startPoint(),
2476  QgsPoint( ( line->startPoint().x() + pX ) * 0.5, ( line->startPoint().y() + pY ) * 0.5 ),
2477  QgsPoint( ( prevX + pX ) * 0.5, ( prevY + pY ) * 0.5 ),
2478  QgsPoint( prevX, prevY ),
2479  segments ) );
2480  out->append( bezier.get() );
2481  }
2482  else
2483  {
2484  std::unique_ptr< QgsLineString > bezier( QgsLineString::fromBezierCurve( out->endPoint(),
2485  QgsPoint( ( out->endPoint().x() + pX ) * 0.5, ( out->endPoint().y() + pY ) * 0.5 ),
2486  QgsPoint( ( prevX + pX ) * 0.5, ( prevY + pY ) * 0.5 ),
2487  QgsPoint( prevX, prevY ),
2488  segments ) );
2489  out->append( bezier.get() );
2490  }
2491 
2492  return out;
2493 }
2494 
2495 std::unique_ptr< QgsAbstractGeometry > roundWavesPrivate( const QgsAbstractGeometry *geom, double wavelength, double amplitude, bool strictWavelength )
2496 {
2497  std::unique_ptr< QgsAbstractGeometry > segmentizedCopy;
2498  if ( QgsWkbTypes::isCurvedType( geom->wkbType() ) )
2499  {
2500  segmentizedCopy.reset( geom->segmentize() );
2501  geom = segmentizedCopy.get();
2502  }
2503 
2505  {
2506  return roundWavesAlongLine( static_cast< const QgsLineString * >( geom ), wavelength, amplitude, strictWavelength );
2507  }
2508  else
2509  {
2510  // polygon
2511  const QgsPolygon *polygon = static_cast< const QgsPolygon * >( geom );
2512  std::unique_ptr< QgsPolygon > result = std::make_unique< QgsPolygon >();
2513 
2514  result->setExteriorRing( roundWavesAlongLine( static_cast< const QgsLineString * >( polygon->exteriorRing() ),
2515  wavelength, amplitude, strictWavelength ).release() );
2516  for ( int i = 0; i < polygon->numInteriorRings(); ++i )
2517  {
2518  result->addInteriorRing( roundWavesAlongLine( static_cast< const QgsLineString * >( polygon->interiorRing( i ) ),
2519  wavelength, amplitude, strictWavelength ).release() );
2520  }
2521 
2522  return result;
2523  }
2524 }
2525 
2526 std::unique_ptr< QgsAbstractGeometry > roundWavesRandomizedPrivate( const QgsAbstractGeometry *geom, double minimumWavelength, double maximumWavelength, double minimumAmplitude, double maximumAmplitude, std::uniform_real_distribution<> &uniformDist, std::mt19937 &mt )
2527 {
2528  std::unique_ptr< QgsAbstractGeometry > segmentizedCopy;
2529  if ( QgsWkbTypes::isCurvedType( geom->wkbType() ) )
2530  {
2531  segmentizedCopy.reset( geom->segmentize() );
2532  geom = segmentizedCopy.get();
2533  }
2534 
2536  {
2537  return roundWavesRandomizedAlongLine( static_cast< const QgsLineString * >( geom ), minimumWavelength, maximumWavelength, minimumAmplitude, maximumAmplitude, uniformDist, mt );
2538  }
2539  else
2540  {
2541  // polygon
2542  const QgsPolygon *polygon = static_cast< const QgsPolygon * >( geom );
2543  std::unique_ptr< QgsPolygon > result = std::make_unique< QgsPolygon >();
2544 
2545  result->setExteriorRing( roundWavesRandomizedAlongLine( static_cast< const QgsLineString * >( polygon->exteriorRing() ),
2546  minimumWavelength, maximumWavelength, minimumAmplitude, maximumAmplitude, uniformDist, mt ).release() );
2547  for ( int i = 0; i < polygon->numInteriorRings(); ++i )
2548  {
2549  result->addInteriorRing( roundWavesRandomizedAlongLine( static_cast< const QgsLineString * >( polygon->interiorRing( i ) ),
2550  minimumWavelength, maximumWavelength, minimumAmplitude, maximumAmplitude, uniformDist, mt ).release() );
2551  }
2552 
2553  return result;
2554  }
2555 }
2556 
2557 QgsGeometry QgsInternalGeometryEngine::roundWaves( double wavelength, double amplitude, bool strictWavelength ) const
2558 {
2559  if ( wavelength < 0 || qgsDoubleNear( wavelength, 0 ) )
2560  return QgsGeometry();
2561 
2562  mLastError.clear();
2563  if ( !mGeometry )
2564  {
2565  return QgsGeometry();
2566  }
2567 
2569  {
2570  return QgsGeometry( mGeometry->clone() ); // point geometry, nothing to do
2571  }
2572 
2573  if ( const QgsGeometryCollection *gc = qgsgeometry_cast< const QgsGeometryCollection *>( mGeometry ) )
2574  {
2575  int numGeom = gc->numGeometries();
2576  QVector< QgsAbstractGeometry * > geometryList;
2577  geometryList.reserve( numGeom );
2578  for ( int i = 0; i < numGeom; ++i )
2579  {
2580  geometryList << roundWavesPrivate( gc->geometryN( i ), wavelength, amplitude, strictWavelength ).release();
2581  }
2582 
2583  QgsGeometry first = QgsGeometry( geometryList.takeAt( 0 ) );
2584  for ( QgsAbstractGeometry *g : std::as_const( geometryList ) )
2585  {
2586  first.addPart( g );
2587  }
2588  return first;
2589  }
2590  else
2591  {
2592  return QgsGeometry( roundWavesPrivate( mGeometry, wavelength, amplitude, strictWavelength ) );
2593  }
2594 }
2595 
2596 QgsGeometry QgsInternalGeometryEngine::roundWavesRandomized( double minimumWavelength, double maximumWavelength, double minimumAmplitude, double maximumAmplitude, unsigned long seed ) const
2597 {
2598  if ( minimumWavelength < 0 || qgsDoubleNear( minimumWavelength, 0 ) || maximumWavelength < 0 || qgsDoubleNear( maximumWavelength, 0 ) || maximumWavelength < minimumWavelength )
2599  return QgsGeometry();
2600 
2601  mLastError.clear();
2602  if ( !mGeometry )
2603  {
2604  return QgsGeometry();
2605  }
2606 
2607  std::random_device rd;
2608  std::mt19937 mt( seed == 0 ? rd() : seed );
2609  std::uniform_real_distribution<> uniformDist( 0, 1 );
2610 
2612  {
2613  return QgsGeometry( mGeometry->clone() ); // point geometry, nothing to do
2614  }
2615 
2616  if ( const QgsGeometryCollection *gc = qgsgeometry_cast< const QgsGeometryCollection *>( mGeometry ) )
2617  {
2618  int numGeom = gc->numGeometries();
2619  QVector< QgsAbstractGeometry * > geometryList;
2620  geometryList.reserve( numGeom );
2621  for ( int i = 0; i < numGeom; ++i )
2622  {
2623  geometryList << roundWavesRandomizedPrivate( gc->geometryN( i ), minimumWavelength, maximumWavelength, minimumAmplitude, maximumAmplitude, uniformDist, mt ).release();
2624  }
2625 
2626  QgsGeometry first = QgsGeometry( geometryList.takeAt( 0 ) );
2627  for ( QgsAbstractGeometry *g : std::as_const( geometryList ) )
2628  {
2629  first.addPart( g );
2630  }
2631  return first;
2632  }
2633  else
2634  {
2635  return QgsGeometry( roundWavesRandomizedPrivate( mGeometry, minimumWavelength, maximumWavelength, minimumAmplitude, maximumAmplitude, uniformDist, mt ) );
2636  }
2637 }
2638 
2639 std::unique_ptr< QgsMultiLineString > dashPatternAlongLine( const QgsLineString *line,
2640  const QVector< double> &pattern,
2644  double patternOffset )
2645 {
2646  const int totalPoints = line->numPoints();
2647  if ( totalPoints < 2 )
2648  return nullptr;
2649 
2650  const int patternSize = pattern.size();
2651 
2652  const double *x = line->xData();
2653  const double *y = line->yData();
2654 
2655  double prevX = *x++;
2656  double prevY = *y++;
2657 
2658  std::unique_ptr< QgsMultiLineString > result = std::make_unique< QgsMultiLineString >();
2659 
2660  QVector< double > outX;
2661  QVector< double > outY;
2662  const double totalLength = line->length();
2663 
2664  double patternLength = 0;
2665  double patternDashLength = 0;
2666  double patternGapLength = 0;
2667  for ( int i = 0; i < pattern.size(); ++i )
2668  {
2669  patternLength += pattern.at( i );
2670  if ( i % 2 == 0 )
2671  patternDashLength += pattern.at( i );
2672  else
2673  patternGapLength += pattern.at( i );
2674  }
2675 
2676  double firstPatternLength = 0;
2677  double firstPatternDashLength = 0;
2678  double firstPatternGapLength = 0;
2679  switch ( startRule )
2680  {
2683  firstPatternLength = patternLength;
2684  firstPatternDashLength = patternDashLength;
2685  firstPatternGapLength = patternGapLength;
2686  break;
2688  firstPatternLength = patternLength - pattern.at( 0 ) * 0.5; // remove half of first dash
2689  firstPatternDashLength = patternDashLength - pattern.at( 0 ) * 0.5;
2690  firstPatternGapLength = patternGapLength;
2691  break;
2693  firstPatternLength = pattern.at( patternSize - 1 ); // start at LAST dash
2694  firstPatternDashLength = 0;
2695  firstPatternGapLength = pattern.at( patternSize - 1 );
2696  break;
2698  firstPatternLength = pattern.at( patternSize - 1 ) * 0.5; // start at half of last dash
2699  firstPatternDashLength = 0;
2700  firstPatternGapLength = pattern.at( patternSize - 1 ) * 0.5;
2701  break;
2702  }
2703 
2704  const bool isSmallEnoughForSinglePattern = ( totalLength - firstPatternLength ) < patternLength * 0.5;
2705 
2706  double lastPatternLength = isSmallEnoughForSinglePattern ? firstPatternLength : patternLength;
2707  double lastPatternDashLength = isSmallEnoughForSinglePattern ? firstPatternDashLength : patternDashLength;
2708  double lastPatternGapLength = isSmallEnoughForSinglePattern ? firstPatternGapLength : patternGapLength;
2709  switch ( endRule )
2710  {
2712  lastPatternLength = 0;
2713  lastPatternDashLength = 0;
2714  lastPatternGapLength = 0;
2715  break;
2717  lastPatternLength -= pattern.at( patternSize - 1 ); // remove last gap
2718  lastPatternGapLength -= pattern.at( patternSize - 1 );
2719  break;
2721  lastPatternLength -= pattern.at( patternSize - 1 ) + pattern.at( patternSize - 2 ) * 0.5; // remove last gap, half of last dash
2722  lastPatternDashLength -= pattern.at( patternSize - 2 ) * 0.5;
2723  lastPatternGapLength -= pattern.at( patternSize - 1 );
2724  break;
2726  lastPatternGapLength = patternGapLength;
2727  break;
2729  lastPatternLength -= pattern.at( patternSize - 1 ) * 0.5; // remove half of last gap
2730  lastPatternGapLength -= pattern.at( patternSize - 1 ) * 0.5;
2731  break;
2732  }
2733 
2734  const double remainingLengthForCompletePatterns = totalLength - ( !isSmallEnoughForSinglePattern ? firstPatternLength : 0 ) - lastPatternLength;
2735  const int middlePatternRepetitions = std::max( static_cast< int >( std::round( remainingLengthForCompletePatterns / patternLength ) ), 0 );
2736 
2737  const double totalUnscaledLengthOfPatterns = ( !isSmallEnoughForSinglePattern ? firstPatternLength : 0 ) + middlePatternRepetitions * patternLength + lastPatternLength;
2738  const double totalUnscaledLengthOfDashes = ( !isSmallEnoughForSinglePattern ? firstPatternDashLength : 0 ) + middlePatternRepetitions * patternDashLength + lastPatternDashLength;
2739  const double totalUnscaledLengthOfGaps = ( !isSmallEnoughForSinglePattern ? firstPatternGapLength : 0 ) + middlePatternRepetitions * patternGapLength + lastPatternGapLength;
2740 
2741  double dashLengthScalingFactor = 1;
2742  double gapLengthScalingFactor = 1;
2743  if ( endRule != Qgis::DashPatternLineEndingRule::NoRule )
2744  {
2745  // calculate scaling factors
2746  const double lengthToShrinkBy = totalUnscaledLengthOfPatterns - totalLength;
2747 
2748  switch ( adjustment )
2749  {
2751  dashLengthScalingFactor = totalLength / totalUnscaledLengthOfPatterns;
2752  gapLengthScalingFactor = dashLengthScalingFactor;
2753  break;
2755  dashLengthScalingFactor = ( totalUnscaledLengthOfDashes - lengthToShrinkBy ) / totalUnscaledLengthOfDashes;
2756  break;
2758  gapLengthScalingFactor = ( totalUnscaledLengthOfGaps - lengthToShrinkBy ) / totalUnscaledLengthOfGaps;
2759  break;
2760  }
2761  }
2762 
2763  dashLengthScalingFactor = std::max( dashLengthScalingFactor, 0.0 );
2764  gapLengthScalingFactor = std::max( gapLengthScalingFactor, 0.0 );
2765 
2766  const int maxPatterns = middlePatternRepetitions + 2;
2767  result->reserve( maxPatterns );
2768 
2769  int patternIndex = 0;
2770  double distanceToNextPointFromStartOfSegment = pattern.at( 0 ) * dashLengthScalingFactor;
2771  bool isDash = true;
2772  switch ( startRule )
2773  {
2776  break;
2778  distanceToNextPointFromStartOfSegment *= 0.5; // skip to half way through first dash
2779  break;
2781  patternIndex = patternSize - 1; // start at last gap
2782  isDash = false;
2783  distanceToNextPointFromStartOfSegment = pattern.at( patternSize - 1 ) * gapLengthScalingFactor;
2784  break;
2786  patternIndex = patternSize - 1; // skip straight to half way through last gap
2787  isDash = false;
2788  distanceToNextPointFromStartOfSegment = 0.5 * pattern.at( patternSize - 1 ) * gapLengthScalingFactor;
2789  break;
2790  }
2791 
2792  const double adjustedOffset = fmod( patternOffset, patternLength );
2793  const double scaledOffset = ( adjustedOffset < 0 ? ( adjustedOffset + patternLength ) : adjustedOffset ) * ( gapLengthScalingFactor + dashLengthScalingFactor ) / 2;
2794  if ( !qgsDoubleNear( scaledOffset, 0 ) )
2795  {
2796  // shuffle pattern along by offset
2797  double remainingOffset = scaledOffset;
2798  while ( remainingOffset > 0 )
2799  {
2800  if ( distanceToNextPointFromStartOfSegment > remainingOffset )
2801  {
2802  distanceToNextPointFromStartOfSegment -= remainingOffset;
2803  break;
2804  }
2805 
2806  remainingOffset -= distanceToNextPointFromStartOfSegment;
2807  isDash = !isDash;
2808  patternIndex++;
2809  if ( patternIndex == patternSize )
2810  patternIndex = 0;
2811 
2812  distanceToNextPointFromStartOfSegment = pattern.at( patternIndex ) * ( isDash ? dashLengthScalingFactor : gapLengthScalingFactor );
2813  }
2814  }
2815 
2816  if ( isDash )
2817  {
2818  outX.append( prevX );
2819  outY.append( prevY );
2820  }
2821 
2822  for ( int i = 1; i < totalPoints; ++i )
2823  {
2824  double thisX = *x++;
2825  double thisY = *y++;
2826 
2827  const double segmentLength = std::sqrt( ( thisX - prevX ) * ( thisX - prevX ) + ( thisY - prevY ) * ( thisY - prevY ) );
2828  while ( distanceToNextPointFromStartOfSegment < segmentLength || qgsDoubleNear( distanceToNextPointFromStartOfSegment, segmentLength ) )
2829  {
2830  // point falls on this segment - truncate to segment length if qgsDoubleNear test was actually > segment length
2831  const double distanceToPoint = std::min( distanceToNextPointFromStartOfSegment, segmentLength );
2832  double pX, pY;
2833  QgsGeometryUtils::pointOnLineWithDistance( prevX, prevY, thisX, thisY, distanceToPoint, pX, pY );
2834 
2835  outX.append( pX );
2836  outY.append( pY );
2837  if ( isDash )
2838  {
2839  result->addGeometry( new QgsLineString( outX, outY ) );
2840  outX.resize( 0 );
2841  outY.resize( 0 );
2842  }
2843 
2844  isDash = !isDash;
2845  patternIndex++;
2846  if ( patternIndex >= patternSize )
2847  patternIndex = 0;
2848 
2849  distanceToNextPointFromStartOfSegment += pattern.at( patternIndex ) * ( isDash ? dashLengthScalingFactor : gapLengthScalingFactor );
2850  }
2851 
2852  if ( isDash )
2853  {
2854  outX.append( thisX );
2855  outY.append( thisY );
2856  }
2857 
2858  prevX = thisX;
2859  prevY = thisY;
2860  distanceToNextPointFromStartOfSegment -= segmentLength;
2861  }
2862 
2863  if ( isDash )
2864  {
2865  outX.append( prevX );
2866  outY.append( prevY );
2867  result->addGeometry( new QgsLineString( outX, outY ) );
2868  }
2869 
2870  return result;
2871 }
2872 
2873 std::unique_ptr< QgsAbstractGeometry > applyDashPatternPrivate( const QgsAbstractGeometry *geom,
2874  const QVector<double> &pattern,
2878  double patternOffset )
2879 {
2880  std::unique_ptr< QgsAbstractGeometry > segmentizedCopy;
2881  if ( QgsWkbTypes::isCurvedType( geom->wkbType() ) )
2882  {
2883  segmentizedCopy.reset( geom->segmentize() );
2884  geom = segmentizedCopy.get();
2885  }
2886 
2888  {
2889  return dashPatternAlongLine( static_cast< const QgsLineString * >( geom ), pattern, startRule, endRule, adjustment, patternOffset );
2890  }
2891  else
2892  {
2893  // polygon
2894  const QgsPolygon *polygon = static_cast< const QgsPolygon * >( geom );
2895  std::unique_ptr< QgsMultiLineString > result = std::make_unique< QgsMultiLineString >();
2896 
2897  std::unique_ptr< QgsMultiLineString > exteriorParts = dashPatternAlongLine( static_cast< const QgsLineString * >( polygon->exteriorRing() ), pattern, startRule, endRule, adjustment, patternOffset );
2898  for ( int i = 0; i < exteriorParts->numGeometries(); ++i )
2899  result->addGeometry( exteriorParts->geometryN( i )->clone() );
2900 
2901  for ( int i = 0; i < polygon->numInteriorRings(); ++i )
2902  {
2903  std::unique_ptr< QgsMultiLineString > ringParts = dashPatternAlongLine( static_cast< const QgsLineString * >( polygon->interiorRing( i ) ), pattern, startRule, endRule, adjustment, patternOffset );
2904  for ( int j = 0; j < ringParts->numGeometries(); ++j )
2905  result->addGeometry( ringParts->geometryN( j )->clone() );
2906  }
2907 
2908  return result;
2909  }
2910 }
2911 
2913 {
2914  if ( pattern.size() < 2 )
2915  return QgsGeometry( mGeometry->clone() );
2916 
2917  mLastError.clear();
2918  if ( !mGeometry || mGeometry->isEmpty() )
2919  {
2920  return QgsGeometry();
2921  }
2922 
2924  {
2925  return QgsGeometry( mGeometry->clone() ); // point geometry, nothing to do
2926  }
2927 
2928  if ( const QgsGeometryCollection *gc = qgsgeometry_cast< const QgsGeometryCollection *>( mGeometry ) )
2929  {
2930  int numGeom = gc->numGeometries();
2931  QVector< QgsAbstractGeometry * > geometryList;
2932  geometryList.reserve( numGeom );
2933  for ( int i = 0; i < numGeom; ++i )
2934  {
2935  geometryList << applyDashPatternPrivate( gc->geometryN( i ), pattern, startRule, endRule, adjustment, patternOffset ).release();
2936  }
2937 
2938  QgsGeometry first = QgsGeometry( geometryList.takeAt( 0 ) );
2939  for ( QgsAbstractGeometry *g : std::as_const( geometryList ) )
2940  {
2941  if ( const QgsGeometryCollection *collection = qgsgeometry_cast< QgsGeometryCollection * >( g ) )
2942  {
2943  for ( int j = 0; j < collection->numGeometries(); ++j )
2944  {
2945  first.addPart( collection->geometryN( j )->clone() );
2946  }
2947  delete collection;
2948  }
2949  else
2950  {
2951  first.addPart( g );
2952  }
2953  }
2954  return first;
2955  }
2956  else
2957  {
2958  return QgsGeometry( applyDashPatternPrivate( mGeometry, pattern, startRule, endRule, adjustment, patternOffset ) );
2959  }
2960 }
DashPatternSizeAdjustment
Dash pattern size adjustment options.
Definition: qgis.h:1382
@ 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:1367
@ 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 QgsAbstractGeometry * clone() const =0
Clones the geometry by performing a deep copy.
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.
Circle geometry type.
Definition: qgscircle.h:44
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.
const QgsAbstractGeometry * geometryN(int n) const
Returns a const reference to a geometry from within the collection.
int partCount() const override
Returns count of parts contained in the geometry.
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:125
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.
static QgsGeometry unaryUnion(const QVector< QgsGeometry > &geometries)
Compute the unary union on a list of geometries.
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:128
static QgsGeometryEngine * createGeometryEngine(const QgsAbstractGeometry *geometry)
Creates and returns a new geometry engine representing the specified geometry.
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:44
QgsPoint startPoint() const override SIP_HOLDGIL
Returns the starting point of the curve.
bool isClosed() const override SIP_HOLDGIL
Returns true if the curve is closed.
const double * yData() const
Returns a const pointer to the y 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.
const double * xData() const
Returns a const pointer to the x vertex data.
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
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:79
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:1990
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)
QgsAbstractGeometry * densifyGeometry(const QgsAbstractGeometry *geom, int extraNodesPerSegment=1, double distance=1)
std::unique_ptr< QgsLineString > squareWavesAlongLine(const QgsLineString *line, double wavelength, const double amplitude, const bool strictWavelength)
std::unique_ptr< QgsCompoundCurve > lineToCurve(const QgsLineString *lineString, double distanceTolerance, double pointSpacingAngleTolerance)
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)
Direction getEdgeDirection(const QgsPoint &p1, const QgsPoint &p2, double maxDev)
Determines the direction of an edge from p1 to p2.
QVector< QgsPointXY > generateSegmentCurve(const QgsPoint &center1, const double radius1, const QgsPoint &center2, const double radius2)
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)
bool dotProductWithinAngleTolerance(double dotProduct, double lowerThreshold, double upperThreshold)
std::pair< bool, std::array< Direction, 4 > > getEdgeDirections(const QgsPolygon *g, double maxDev)
Checks whether the polygon consists of four nearly axis-parallel sides.
double squareness(QgsLineString *ring, double lowerThreshold, double upperThreshold)
QgsAbstractGeometry * orthogonalizeGeom(const QgsAbstractGeometry *geom, int maxIterations, double tolerance, 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.
std::unique_ptr< QgsLineString > roundWavesAlongLine(const QgsLineString *line, double wavelength, const double amplitude, const bool strictWavelength)
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)
QgsLineString * doOrthogonalize(QgsLineString *ring, int iterations, double tolerance, double lowerThreshold, double upperThreshold)
std::unique_ptr< QgsLineString > triangularWavesAlongLine(const QgsLineString *line, double wavelength, const double amplitude, const bool strictWavelength)
QgsLineString * doDensify(const QgsLineString *ring, int extraNodesPerSegment=-1, double distance=1)
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