QGIS API Documentation  3.27.0-Master (597e8eebd4)
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< QgsCurve > lineToCurve( const QgsCurve *curve, double distanceTolerance,
1336  double pointSpacingAngleTolerance )
1337 {
1338  const QgsWkbTypes::Type flatType = QgsWkbTypes::flatType( curve->wkbType() );
1339  if ( flatType == QgsWkbTypes::CircularString )
1340  {
1341  // already curved, so just return a copy
1342  std::unique_ptr< QgsCircularString > out;
1343  out.reset( qgsgeometry_cast< QgsCircularString * >( curve )->clone() );
1344  return out;
1345  }
1346  else if ( flatType == QgsWkbTypes::CompoundCurve )
1347  {
1348  std::unique_ptr< QgsCompoundCurve > out = std::make_unique< QgsCompoundCurve >();
1349  const QgsCompoundCurve *in = qgsgeometry_cast< const QgsCompoundCurve * >( curve );
1350  for ( int i = 0; i < in->nCurves(); i ++ )
1351  {
1352  std::unique_ptr< QgsCurve > processed = lineToCurve( in->curveAt( i ), distanceTolerance, pointSpacingAngleTolerance );
1353  if ( processed )
1354  {
1355  // check the simplified type ref of the geometry to see if it consists of just a single linestring component.
1356  // If so, we extract this line string along, so that we don't create geometry components which consits of a
1357  // compound curve with only linestring parts.
1358  if ( QgsWkbTypes::flatType( processed->simplifiedTypeRef()->wkbType() ) == QgsWkbTypes::LineString )
1359  {
1360  out->addCurve( qgsgeometry_cast< const QgsLineString * >( processed->simplifiedTypeRef() )->clone() );
1361  }
1362  else
1363  {
1364  out->addCurve( processed.release() );
1365  }
1366  }
1367  }
1368  return out;
1369  }
1370  else if ( flatType == QgsWkbTypes::LineString )
1371  {
1372  const QgsLineString *lineString = qgsgeometry_cast< QgsLineString * >( curve );
1373 
1374  std::unique_ptr< QgsCompoundCurve > out = std::make_unique< QgsCompoundCurve >();
1375 
1376  /* Minimum number of edges, per quadrant, required to define an arc */
1377  const unsigned int minQuadEdges = 2;
1378 
1379  /* Die on null input */
1380  if ( !lineString )
1381  return nullptr;
1382 
1383  /* Null on empty input? */
1384  if ( lineString->nCoordinates() == 0 )
1385  return nullptr;
1386 
1387  /* We can't desegmentize anything shorter than four points */
1388  if ( lineString->nCoordinates() < 4 )
1389  {
1390  out->addCurve( lineString->clone() );
1391  return out;
1392  }
1393 
1394  /* Allocate our result array of vertices that are part of arcs */
1395  int numEdges = lineString->nCoordinates() - 1;
1396  QVector< int > edgesInArcs( numEdges + 1, 0 );
1397 
1398  auto arcAngle = []( const QgsPoint & a, const QgsPoint & b, const QgsPoint & c )->double
1399  {
1400  double abX = b.x() - a.x();
1401  double abY = b.y() - a.y();
1402 
1403  double cbX = b.x() - c.x();
1404  double cbY = b.y() - c.y();
1405 
1406  double dot = ( abX * cbX + abY * cbY ); /* dot product */
1407  double cross = ( abX * cbY - abY * cbX ); /* cross product */
1408 
1409  double alpha = std::atan2( cross, dot );
1410 
1411  return alpha;
1412  };
1413 
1414  /* We make a candidate arc of the first two edges, */
1415  /* And then see if the next edge follows it */
1416  int i = 0;
1417  int j = 0;
1418  int k = 0;
1419  int currentArc = 1;
1420  QgsPoint a1;
1421  QgsPoint a2;
1422  QgsPoint a3;
1423  QgsPoint b;
1424  double centerX = 0.0;
1425  double centerY = 0.0;
1426  double radius = 0;
1427 
1428  while ( i < numEdges - 2 )
1429  {
1430  unsigned int arcEdges = 0;
1431  double numQuadrants = 0;
1432  double angle;
1433 
1434  bool foundArc = false;
1435  /* Make candidate arc */
1436  a1 = lineString->pointN( i );
1437  a2 = lineString->pointN( i + 1 );
1438  a3 = lineString->pointN( i + 2 );
1439  QgsPoint first = a1;
1440 
1441  for ( j = i + 3; j < numEdges + 1; j++ )
1442  {
1443  b = lineString->pointN( j );
1444 
1445  /* Does this point fall on our candidate arc? */
1446  if ( QgsGeometryUtils::pointContinuesArc( a1, a2, a3, b, distanceTolerance, pointSpacingAngleTolerance ) )
1447  {
1448  /* Yes. Mark this edge and the two preceding it as arc components */
1449  foundArc = true;
1450  for ( k = j - 1; k > j - 4; k-- )
1451  edgesInArcs[k] = currentArc;
1452  }
1453  else
1454  {
1455  /* No. So we're done with this candidate arc */
1456  currentArc++;
1457  break;
1458  }
1459 
1460  a1 = a2;
1461  a2 = a3;
1462  a3 = b;
1463  }
1464  /* Jump past all the edges that were added to the arc */
1465  if ( foundArc )
1466  {
1467  /* Check if an arc was composed by enough edges to be
1468  * really considered an arc
1469  * See http://trac.osgeo.org/postgis/ticket/2420
1470  */
1471  arcEdges = j - 1 - i;
1472  if ( first.x() == b.x() && first.y() == b.y() )
1473  {
1474  numQuadrants = 4;
1475  }
1476  else
1477  {
1478  QgsGeometryUtils::circleCenterRadius( first, b, a1, radius, centerX, centerY );
1479 
1480  angle = arcAngle( first, QgsPoint( centerX, centerY ), b );
1481  int p2Side = QgsGeometryUtils::leftOfLine( b.x(), b.y(), first.x(), first.y(), a1.x(), a1.y() );
1482  if ( p2Side >= 0 )
1483  angle = -angle;
1484 
1485  if ( angle < 0 )
1486  angle = 2 * M_PI + angle;
1487  numQuadrants = ( 4 * angle ) / ( 2 * M_PI );
1488  }
1489  /* a1 is first point, b is last point */
1490  if ( arcEdges < minQuadEdges * numQuadrants )
1491  {
1492  // LWDEBUGF( 4, "Not enough edges for a %g quadrants arc, %g needed", num_quadrants, min_quad_edges * num_quadrants );
1493  for ( k = j - 1; k >= i; k-- )
1494  edgesInArcs[k] = 0;
1495  }
1496 
1497  i = j - 1;
1498  }
1499  else
1500  {
1501  /* Mark this edge as a linear edge */
1502  edgesInArcs[i] = 0;
1503  i = i + 1;
1504  }
1505  }
1506 
1507  int start = 0;
1508  int end = 0;
1509  /* non-zero if edge is part of an arc */
1510  int edgeType = edgesInArcs[0];
1511 
1512  auto addPointsToCurve = [ lineString, &out ]( int start, int end, int type )
1513  {
1514  if ( type == 0 )
1515  {
1516  // straight segment
1517  QVector< QgsPoint > points;
1518  for ( int j = start; j < end + 2; ++ j )
1519  {
1520  points.append( lineString->pointN( j ) );
1521  }
1522  std::unique_ptr< QgsCurve > straightSegment = std::make_unique< QgsLineString >( points );
1523  out->addCurve( straightSegment.release() );
1524  }
1525  else
1526  {
1527  // curved segment
1528  QVector< QgsPoint > points;
1529  points.append( lineString->pointN( start ) );
1530  points.append( lineString->pointN( ( start + end + 1 ) / 2 ) );
1531  points.append( lineString->pointN( end + 1 ) );
1532  std::unique_ptr< QgsCircularString > curvedSegment = std::make_unique< QgsCircularString >();
1533  curvedSegment->setPoints( points );
1534  out->addCurve( curvedSegment.release() );
1535  }
1536  };
1537 
1538  for ( int i = 1; i < numEdges; i++ )
1539  {
1540  if ( edgeType != edgesInArcs[i] )
1541  {
1542  end = i - 1;
1543  addPointsToCurve( start, end, edgeType );
1544  start = i;
1545  edgeType = edgesInArcs[i];
1546  }
1547  }
1548 
1549  /* Roll out last item */
1550  end = numEdges - 1;
1551  addPointsToCurve( start, end, edgeType );
1552 
1553  // return a simplified type if it doesn't need to be a compound curve resu
1554  if ( QgsWkbTypes::flatType( out->simplifiedTypeRef()->wkbType() ) == QgsWkbTypes::CircularString )
1555  {
1556  std::unique_ptr< QgsCircularString > res;
1557  res.reset( qgsgeometry_cast< const QgsCircularString * >( out->simplifiedTypeRef() )->clone() );
1558  return res;
1559  }
1560 
1561  return out;
1562  }
1563 
1564  return nullptr;
1565 }
1566 
1567 std::unique_ptr< QgsAbstractGeometry > convertGeometryToCurves( const QgsAbstractGeometry *geom, double distanceTolerance, double angleTolerance )
1568 {
1570  {
1571  return lineToCurve( qgsgeometry_cast< const QgsCurve * >( geom ), distanceTolerance, angleTolerance );
1572  }
1573  else
1574  {
1575  // polygon
1576  const QgsCurvePolygon *polygon = qgsgeometry_cast< const QgsCurvePolygon * >( geom );
1577  std::unique_ptr< QgsCurvePolygon > result = std::make_unique< QgsCurvePolygon>();
1578 
1579  result->setExteriorRing( lineToCurve( polygon->exteriorRing(),
1580  distanceTolerance, angleTolerance ).release() );
1581  for ( int i = 0; i < polygon->numInteriorRings(); ++i )
1582  {
1583  result->addInteriorRing( lineToCurve( polygon->interiorRing( i ),
1584  distanceTolerance, angleTolerance ).release() );
1585  }
1586 
1587  return result;
1588  }
1589 }
1590 
1591 QgsGeometry QgsInternalGeometryEngine::convertToCurves( double distanceTolerance, double angleTolerance ) const
1592 {
1593  mLastError.clear();
1594  if ( !mGeometry )
1595  {
1596  return QgsGeometry();
1597  }
1598 
1600  {
1601  return QgsGeometry( mGeometry->clone() ); // point geometry, nothing to do
1602  }
1603 
1604  if ( const QgsGeometryCollection *gc = qgsgeometry_cast< const QgsGeometryCollection *>( mGeometry ) )
1605  {
1606  int numGeom = gc->numGeometries();
1607  QVector< QgsAbstractGeometry * > geometryList;
1608  geometryList.reserve( numGeom );
1609  for ( int i = 0; i < numGeom; ++i )
1610  {
1611  geometryList << convertGeometryToCurves( gc->geometryN( i ), distanceTolerance, angleTolerance ).release();
1612  }
1613 
1614  QgsGeometry first = QgsGeometry( geometryList.takeAt( 0 ) );
1615  for ( QgsAbstractGeometry *g : std::as_const( geometryList ) )
1616  {
1617  first.addPart( g );
1618  }
1619  return first;
1620  }
1621  else
1622  {
1623  return QgsGeometry( convertGeometryToCurves( mGeometry, distanceTolerance, angleTolerance ) );
1624  }
1625 }
1626 
1627 QgsGeometry QgsInternalGeometryEngine::orientedMinimumBoundingBox( double &area, double &angle, double &width, double &height ) const
1628 {
1629  mLastError.clear();
1630 
1631  QgsRectangle minRect;
1632  area = std::numeric_limits<double>::max();
1633  angle = 0;
1634  width = std::numeric_limits<double>::max();
1635  height = std::numeric_limits<double>::max();
1636 
1637  if ( !mGeometry || mGeometry->nCoordinates() < 2 )
1638  return QgsGeometry();
1639 
1640  std::unique_ptr< QgsGeometryEngine >engine( QgsGeometry::createGeometryEngine( mGeometry ) );
1641  std::unique_ptr< QgsAbstractGeometry > hull( engine->convexHull( &mLastError ) );
1642  if ( !hull )
1643  return QgsGeometry();
1644 
1645  QgsVertexId vertexId;
1646  QgsPoint pt0;
1647  QgsPoint pt1;
1648  QgsPoint pt2;
1649  // get first point
1650  hull->nextVertex( vertexId, pt0 );
1651  pt1 = pt0;
1652  double totalRotation = 0;
1653  while ( hull->nextVertex( vertexId, pt2 ) )
1654  {
1655  double currentAngle = QgsGeometryUtils::lineAngle( pt1.x(), pt1.y(), pt2.x(), pt2.y() );
1656  double rotateAngle = 180.0 / M_PI * currentAngle;
1657  totalRotation += rotateAngle;
1658 
1659  QTransform t = QTransform::fromTranslate( pt0.x(), pt0.y() );
1660  t.rotate( rotateAngle );
1661  t.translate( -pt0.x(), -pt0.y() );
1662 
1663  hull->transform( t );
1664 
1665  QgsRectangle bounds = hull->boundingBox();
1666  double currentArea = bounds.width() * bounds.height();
1667  if ( currentArea < area )
1668  {
1669  minRect = bounds;
1670  area = currentArea;
1671  angle = totalRotation;
1672  width = bounds.width();
1673  height = bounds.height();
1674  }
1675 
1676  pt1 = hull->vertexAt( vertexId );
1677  }
1678 
1679  QgsGeometry minBounds = QgsGeometry::fromRect( minRect );
1680  minBounds.rotate( angle, QgsPointXY( pt0.x(), pt0.y() ) );
1681 
1682  if ( width > height )
1683  {
1684  width = minRect.height();
1685  height = minRect.width();
1686  angle = angle + 90.0;
1687  }
1688 
1689  // constrain angle to 0 - 180
1690  if ( angle > 180.0 )
1691  angle = std::fmod( angle, 180.0 );
1692 
1693  return minBounds;
1694 }
1695 
1696 std::unique_ptr< QgsLineString > triangularWavesAlongLine( const QgsLineString *line, double wavelength, const double amplitude, const bool strictWavelength )
1697 {
1698  const int totalPoints = line->numPoints();
1699  if ( totalPoints < 2 )
1700  return nullptr;
1701 
1702  const double *x = line->xData();
1703  const double *y = line->yData();
1704 
1705  double prevX = *x++;
1706  double prevY = *y++;
1707 
1708  QVector< double > outX;
1709  QVector< double > outY;
1710  const double totalLength = line->length();
1711 
1712  const int maxRepetitions = std::ceil( totalLength / wavelength );
1713  if ( !strictWavelength )
1714  wavelength = totalLength / maxRepetitions;
1715 
1716  const int estimatedPoints = maxRepetitions * 2 + 2;
1717  outX.reserve( estimatedPoints );
1718  outY.reserve( estimatedPoints );
1719  outX.append( prevX );
1720  outY.append( prevY );
1721 
1722  double distanceToNextPointFromStartOfSegment = wavelength / 4;
1723  int side = -1;
1724  for ( int i = 1; i < totalPoints; ++i )
1725  {
1726  double thisX = *x++;
1727  double thisY = *y++;
1728 
1729  const double segmentAngleRadians = QgsGeometryUtils::lineAngle( prevX, prevY, thisX, thisY );
1730  const double segmentLength = std::sqrt( ( thisX - prevX ) * ( thisX - prevX ) + ( thisY - prevY ) * ( thisY - prevY ) );
1731  while ( distanceToNextPointFromStartOfSegment < segmentLength || qgsDoubleNear( distanceToNextPointFromStartOfSegment, segmentLength ) )
1732  {
1733  // point falls on this segment - truncate to segment length if qgsDoubleNear test was actually > segment length
1734  const double distanceToPoint = std::min( distanceToNextPointFromStartOfSegment, segmentLength );
1735  double pX, pY;
1736  QgsGeometryUtils::pointOnLineWithDistance( prevX, prevY, thisX, thisY, distanceToPoint, pX, pY );
1737 
1738  // project point on line out by amplitude
1739  const double outPointX = pX + side * amplitude * std::sin( segmentAngleRadians + M_PI_2 );
1740  const double outPointY = pY + side * amplitude * std::cos( segmentAngleRadians + M_PI_2 );
1741 
1742  outX.append( outPointX );
1743  outY.append( outPointY );
1744 
1745  distanceToNextPointFromStartOfSegment += wavelength / 2;
1746  side = -side;
1747  }
1748 
1749  prevX = thisX;
1750  prevY = thisY;
1751  distanceToNextPointFromStartOfSegment -= segmentLength;
1752  }
1753 
1754  outX.append( prevX );
1755  outY.append( prevY );
1756 
1757  return std::make_unique< QgsLineString >( outX, outY );
1758 }
1759 
1760 std::unique_ptr< QgsLineString > triangularWavesRandomizedAlongLine( const QgsLineString *line,
1761  const double minimumWavelength, const double maximumWavelength,
1762  const double minimumAmplitude, const double maximumAmplitude,
1763  std::uniform_real_distribution<> &uniformDist, std::mt19937 &mt )
1764 {
1765  const int totalPoints = line->numPoints();
1766  if ( totalPoints < 2 )
1767  return nullptr;
1768 
1769  const double *x = line->xData();
1770  const double *y = line->yData();
1771 
1772  double prevX = *x++;
1773  double prevY = *y++;
1774 
1775  QVector< double > outX;
1776  QVector< double > outY;
1777  const double totalLength = line->length();
1778 
1779  const int maxRepetitions = std::ceil( totalLength / minimumWavelength );
1780 
1781  const int estimatedPoints = maxRepetitions * 2 + 2;
1782  outX.reserve( estimatedPoints );
1783  outY.reserve( estimatedPoints );
1784  outX.append( prevX );
1785  outY.append( prevY );
1786 
1787  double wavelength = uniformDist( mt ) * ( maximumWavelength - minimumWavelength ) + minimumWavelength;
1788  double distanceToNextPointFromStartOfSegment = wavelength / 4;
1789  double amplitude = uniformDist( mt ) * ( maximumAmplitude - minimumAmplitude ) + minimumAmplitude;
1790 
1791  int side = -1;
1792  for ( int i = 1; i < totalPoints; ++i )
1793  {
1794  double thisX = *x++;
1795  double thisY = *y++;
1796 
1797  const double segmentAngleRadians = QgsGeometryUtils::lineAngle( prevX, prevY, thisX, thisY );
1798  const double segmentLength = std::sqrt( ( thisX - prevX ) * ( thisX - prevX ) + ( thisY - prevY ) * ( thisY - prevY ) );
1799  while ( distanceToNextPointFromStartOfSegment < segmentLength || qgsDoubleNear( distanceToNextPointFromStartOfSegment, segmentLength ) )
1800  {
1801  // point falls on this segment - truncate to segment length if qgsDoubleNear test was actually > segment length
1802  const double distanceToPoint = std::min( distanceToNextPointFromStartOfSegment, segmentLength );
1803  double pX, pY;
1804  QgsGeometryUtils::pointOnLineWithDistance( prevX, prevY, thisX, thisY, distanceToPoint, pX, pY );
1805 
1806  // project point on line out by amplitude
1807  const double outPointX = pX + side * amplitude * std::sin( segmentAngleRadians + M_PI_2 );
1808  const double outPointY = pY + side * amplitude * std::cos( segmentAngleRadians + M_PI_2 );
1809  outX.append( outPointX );
1810  outY.append( outPointY );
1811 
1812  wavelength = uniformDist( mt ) * ( maximumWavelength - minimumWavelength ) + minimumWavelength;
1813  amplitude = uniformDist( mt ) * ( maximumAmplitude - minimumAmplitude ) + minimumAmplitude;
1814 
1815  distanceToNextPointFromStartOfSegment += wavelength / 2;
1816  side = -side;
1817  }
1818 
1819  prevX = thisX;
1820  prevY = thisY;
1821  distanceToNextPointFromStartOfSegment -= segmentLength;
1822  }
1823 
1824  outX.append( prevX );
1825  outY.append( prevY );
1826 
1827  return std::make_unique< QgsLineString >( outX, outY );
1828 }
1829 
1830 std::unique_ptr< QgsAbstractGeometry > triangularWavesPrivate( const QgsAbstractGeometry *geom, double wavelength, double amplitude, bool strictWavelength )
1831 {
1832  std::unique_ptr< QgsAbstractGeometry > segmentizedCopy;
1833  if ( QgsWkbTypes::isCurvedType( geom->wkbType() ) )
1834  {
1835  segmentizedCopy.reset( geom->segmentize() );
1836  geom = segmentizedCopy.get();
1837  }
1838 
1840  {
1841  return triangularWavesAlongLine( static_cast< const QgsLineString * >( geom ), wavelength, amplitude, strictWavelength );
1842  }
1843  else
1844  {
1845  // polygon
1846  const QgsPolygon *polygon = static_cast< const QgsPolygon * >( geom );
1847  std::unique_ptr< QgsPolygon > result = std::make_unique< QgsPolygon >();
1848 
1849  result->setExteriorRing( triangularWavesAlongLine( static_cast< const QgsLineString * >( polygon->exteriorRing() ),
1850  wavelength, amplitude, strictWavelength ).release() );
1851  for ( int i = 0; i < polygon->numInteriorRings(); ++i )
1852  {
1853  result->addInteriorRing( triangularWavesAlongLine( static_cast< const QgsLineString * >( polygon->interiorRing( i ) ),
1854  wavelength, amplitude, strictWavelength ).release() );
1855  }
1856 
1857  return result;
1858  }
1859 }
1860 
1861 std::unique_ptr< QgsAbstractGeometry > triangularWavesRandomizedPrivate( const QgsAbstractGeometry *geom, double minimumWavelength, double maximumWavelength, double minimumAmplitude, double maximumAmplitude, std::uniform_real_distribution<> &uniformDist, std::mt19937 &mt )
1862 {
1863  std::unique_ptr< QgsAbstractGeometry > segmentizedCopy;
1864  if ( QgsWkbTypes::isCurvedType( geom->wkbType() ) )
1865  {
1866  segmentizedCopy.reset( geom->segmentize() );
1867  geom = segmentizedCopy.get();
1868  }
1869 
1871  {
1872  return triangularWavesRandomizedAlongLine( static_cast< const QgsLineString * >( geom ), minimumWavelength, maximumWavelength, minimumAmplitude, maximumAmplitude, uniformDist, mt );
1873  }
1874  else
1875  {
1876  // polygon
1877  const QgsPolygon *polygon = static_cast< const QgsPolygon * >( geom );
1878  std::unique_ptr< QgsPolygon > result = std::make_unique< QgsPolygon >();
1879 
1880  result->setExteriorRing( triangularWavesRandomizedAlongLine( static_cast< const QgsLineString * >( polygon->exteriorRing() ),
1881  minimumWavelength, maximumWavelength, minimumAmplitude, maximumAmplitude, uniformDist, mt ).release() );
1882  for ( int i = 0; i < polygon->numInteriorRings(); ++i )
1883  {
1884  result->addInteriorRing( triangularWavesRandomizedAlongLine( static_cast< const QgsLineString * >( polygon->interiorRing( i ) ),
1885  minimumWavelength, maximumWavelength, minimumAmplitude, maximumAmplitude, uniformDist, mt ).release() );
1886  }
1887 
1888  return result;
1889  }
1890 }
1891 
1892 QgsGeometry QgsInternalGeometryEngine::triangularWaves( double wavelength, double amplitude, bool strictWavelength ) const
1893 {
1894  if ( wavelength < 0 || qgsDoubleNear( wavelength, 0 ) )
1895  return QgsGeometry();
1896 
1897  mLastError.clear();
1898  if ( !mGeometry )
1899  {
1900  return QgsGeometry();
1901  }
1902 
1904  {
1905  return QgsGeometry( mGeometry->clone() ); // point geometry, nothing to do
1906  }
1907 
1908  if ( const QgsGeometryCollection *gc = qgsgeometry_cast< const QgsGeometryCollection *>( mGeometry ) )
1909  {
1910  int numGeom = gc->numGeometries();
1911  QVector< QgsAbstractGeometry * > geometryList;
1912  geometryList.reserve( numGeom );
1913  for ( int i = 0; i < numGeom; ++i )
1914  {
1915  geometryList << triangularWavesPrivate( gc->geometryN( i ), wavelength, amplitude, strictWavelength ).release();
1916  }
1917 
1918  QgsGeometry first = QgsGeometry( geometryList.takeAt( 0 ) );
1919  for ( QgsAbstractGeometry *g : std::as_const( geometryList ) )
1920  {
1921  first.addPart( g );
1922  }
1923  return first;
1924  }
1925  else
1926  {
1927  return QgsGeometry( triangularWavesPrivate( mGeometry, wavelength, amplitude, strictWavelength ) );
1928  }
1929 }
1930 
1931 QgsGeometry QgsInternalGeometryEngine::triangularWavesRandomized( double minimumWavelength, double maximumWavelength, double minimumAmplitude, double maximumAmplitude, unsigned long seed ) const
1932 {
1933  if ( minimumWavelength < 0 || qgsDoubleNear( minimumWavelength, 0 ) || maximumWavelength < 0 || qgsDoubleNear( maximumWavelength, 0 ) || maximumWavelength < minimumWavelength )
1934  return QgsGeometry();
1935 
1936  mLastError.clear();
1937  if ( !mGeometry )
1938  {
1939  return QgsGeometry();
1940  }
1941 
1942  std::random_device rd;
1943  std::mt19937 mt( seed == 0 ? rd() : seed );
1944  std::uniform_real_distribution<> uniformDist( 0, 1 );
1945 
1947  {
1948  return QgsGeometry( mGeometry->clone() ); // point geometry, nothing to do
1949  }
1950 
1951  if ( const QgsGeometryCollection *gc = qgsgeometry_cast< const QgsGeometryCollection *>( mGeometry ) )
1952  {
1953  int numGeom = gc->numGeometries();
1954  QVector< QgsAbstractGeometry * > geometryList;
1955  geometryList.reserve( numGeom );
1956  for ( int i = 0; i < numGeom; ++i )
1957  {
1958  geometryList << triangularWavesRandomizedPrivate( gc->geometryN( i ), minimumWavelength, maximumWavelength, minimumAmplitude, maximumAmplitude, uniformDist, mt ).release();
1959  }
1960 
1961  QgsGeometry first = QgsGeometry( geometryList.takeAt( 0 ) );
1962  for ( QgsAbstractGeometry *g : std::as_const( geometryList ) )
1963  {
1964  first.addPart( g );
1965  }
1966  return first;
1967  }
1968  else
1969  {
1970  return QgsGeometry( triangularWavesRandomizedPrivate( mGeometry, minimumWavelength, maximumWavelength, minimumAmplitude, maximumAmplitude, uniformDist, mt ) );
1971  }
1972 }
1973 
1974 std::unique_ptr< QgsLineString > squareWavesAlongLine( const QgsLineString *line, double wavelength, const double amplitude, const bool strictWavelength )
1975 {
1976  const int totalPoints = line->numPoints();
1977  if ( totalPoints < 2 )
1978  return nullptr;
1979 
1980  const double *x = line->xData();
1981  const double *y = line->yData();
1982 
1983  double prevX = *x++;
1984  double prevY = *y++;
1985 
1986  QVector< double > outX;
1987  QVector< double > outY;
1988  const double totalLength = line->length();
1989 
1990  const int maxRepetitions = std::ceil( totalLength / wavelength );
1991  if ( !strictWavelength )
1992  wavelength = totalLength / maxRepetitions;
1993 
1994  const int estimatedPoints = maxRepetitions * 4 + 2;
1995  outX.reserve( estimatedPoints );
1996  outY.reserve( estimatedPoints );
1997  outX.append( prevX );
1998  outY.append( prevY );
1999 
2000  const double startSegmentAngleRadians = QgsGeometryUtils::lineAngle( prevX, prevY, *x, *y );
2001  outX.append( prevX - amplitude * std::sin( startSegmentAngleRadians + M_PI_2 ) );
2002  outY.append( prevY - amplitude * std::cos( startSegmentAngleRadians + M_PI_2 ) );
2003 
2004  double distanceToNextPointFromStartOfSegment = wavelength / 2;
2005 
2006  int side = -1;
2007  for ( int i = 1; i < totalPoints; ++i )
2008  {
2009  double thisX = *x++;
2010  double thisY = *y++;
2011 
2012  const double segmentAngleRadians = QgsGeometryUtils::lineAngle( prevX, prevY, thisX, thisY );
2013  const double segmentLength = std::sqrt( ( thisX - prevX ) * ( thisX - prevX ) + ( thisY - prevY ) * ( thisY - prevY ) );
2014  while ( distanceToNextPointFromStartOfSegment < segmentLength || qgsDoubleNear( distanceToNextPointFromStartOfSegment, segmentLength ) )
2015  {
2016  // point falls on this segment - truncate to segment length if qgsDoubleNear test was actually > segment length
2017  const double distanceToPoint = std::min( distanceToNextPointFromStartOfSegment, segmentLength );
2018  double pX, pY;
2019  QgsGeometryUtils::pointOnLineWithDistance( prevX, prevY, thisX, thisY, distanceToPoint, pX, pY );
2020 
2021  // project point on line out by amplitude
2022  const double sinAngle = std::sin( segmentAngleRadians + M_PI_2 );
2023  const double cosAngle = std::cos( segmentAngleRadians + M_PI_2 );
2024  outX.append( pX + side * amplitude * sinAngle );
2025  outY.append( pY + side * amplitude * cosAngle );
2026  outX.append( pX - side * amplitude * sinAngle );
2027  outY.append( pY - side * amplitude * cosAngle );
2028 
2029  distanceToNextPointFromStartOfSegment += wavelength / 2;
2030  side = -side;
2031  }
2032 
2033  prevX = thisX;
2034  prevY = thisY;
2035  distanceToNextPointFromStartOfSegment -= segmentLength;
2036  }
2037 
2038  // replace last point, which will be projected out to amplitude of wave, with the actual end point of the line
2039  outX.pop_back();
2040  outY.pop_back();
2041  outX.append( prevX );
2042  outY.append( prevY );
2043 
2044  return std::make_unique< QgsLineString >( outX, outY );
2045 }
2046 
2047 std::unique_ptr< QgsLineString > squareWavesRandomizedAlongLine( const QgsLineString *line,
2048  const double minimumWavelength, const double maximumWavelength,
2049  const double minimumAmplitude, const double maximumAmplitude,
2050  std::uniform_real_distribution<> &uniformDist, std::mt19937 &mt )
2051 {
2052  const int totalPoints = line->numPoints();
2053  if ( totalPoints < 2 )
2054  return nullptr;
2055 
2056  const double *x = line->xData();
2057  const double *y = line->yData();
2058 
2059  double prevX = *x++;
2060  double prevY = *y++;
2061 
2062  QVector< double > outX;
2063  QVector< double > outY;
2064  const double totalLength = line->length();
2065 
2066  const int maxRepetitions = std::ceil( totalLength / minimumWavelength );
2067 
2068  const int estimatedPoints = maxRepetitions * 4 + 2;
2069  outX.reserve( estimatedPoints );
2070  outY.reserve( estimatedPoints );
2071  outX.append( prevX );
2072  outY.append( prevY );
2073 
2074  double amplitude = uniformDist( mt ) * ( maximumAmplitude - minimumAmplitude ) + minimumAmplitude;
2075 
2076  double segmentAngleRadians = QgsGeometryUtils::lineAngle( prevX, prevY, *x, *y );
2077  outX.append( prevX - amplitude * std::sin( segmentAngleRadians + M_PI_2 ) );
2078  outY.append( prevY - amplitude * std::cos( segmentAngleRadians + M_PI_2 ) );
2079 
2080  double wavelength = uniformDist( mt ) * ( maximumWavelength - minimumWavelength ) + minimumWavelength;
2081  double distanceToNextPointFromStartOfSegment = wavelength / 2;
2082 
2083  int side = -1;
2084  for ( int i = 1; i < totalPoints; ++i )
2085  {
2086  double thisX = *x++;
2087  double thisY = *y++;
2088 
2089  segmentAngleRadians = QgsGeometryUtils::lineAngle( prevX, prevY, thisX, thisY );
2090  const double segmentLength = std::sqrt( ( thisX - prevX ) * ( thisX - prevX ) + ( thisY - prevY ) * ( thisY - prevY ) );
2091  while ( distanceToNextPointFromStartOfSegment < segmentLength || qgsDoubleNear( distanceToNextPointFromStartOfSegment, segmentLength ) )
2092  {
2093  // point falls on this segment - truncate to segment length if qgsDoubleNear test was actually > segment length
2094  const double distanceToPoint = std::min( distanceToNextPointFromStartOfSegment, segmentLength );
2095  double pX, pY;
2096  QgsGeometryUtils::pointOnLineWithDistance( prevX, prevY, thisX, thisY, distanceToPoint, pX, pY );
2097 
2098  // project point on line out by amplitude
2099  const double sinAngle = std::sin( segmentAngleRadians + M_PI_2 );
2100  const double cosAngle = std::cos( segmentAngleRadians + M_PI_2 );
2101  outX.append( pX + side * amplitude * sinAngle );
2102  outY.append( pY + side * amplitude * cosAngle );
2103 
2104  amplitude = uniformDist( mt ) * ( maximumAmplitude - minimumAmplitude ) + minimumAmplitude;
2105  outX.append( pX - side * amplitude * sinAngle );
2106  outY.append( pY - side * amplitude * cosAngle );
2107 
2108  wavelength = uniformDist( mt ) * ( maximumWavelength - minimumWavelength ) + minimumWavelength;
2109  distanceToNextPointFromStartOfSegment += wavelength / 2;
2110  side = -side;
2111  }
2112 
2113  prevX = thisX;
2114  prevY = thisY;
2115  distanceToNextPointFromStartOfSegment -= segmentLength;
2116  }
2117 
2118  outX.append( prevX + side * amplitude * std::sin( segmentAngleRadians + M_PI_2 ) );
2119  outY.append( prevY + side * amplitude * std::cos( segmentAngleRadians + M_PI_2 ) );
2120  outX.append( prevX );
2121  outY.append( prevY );
2122 
2123  return std::make_unique< QgsLineString >( outX, outY );
2124 }
2125 
2126 std::unique_ptr< QgsAbstractGeometry > squareWavesPrivate( const QgsAbstractGeometry *geom, double wavelength, double amplitude, bool strictWavelength )
2127 {
2128  std::unique_ptr< QgsAbstractGeometry > segmentizedCopy;
2129  if ( QgsWkbTypes::isCurvedType( geom->wkbType() ) )
2130  {
2131  segmentizedCopy.reset( geom->segmentize() );
2132  geom = segmentizedCopy.get();
2133  }
2134 
2136  {
2137  return squareWavesAlongLine( static_cast< const QgsLineString * >( geom ), wavelength, amplitude, strictWavelength );
2138  }
2139  else
2140  {
2141  // polygon
2142  const QgsPolygon *polygon = static_cast< const QgsPolygon * >( geom );
2143  std::unique_ptr< QgsPolygon > result = std::make_unique< QgsPolygon >();
2144 
2145  result->setExteriorRing( squareWavesAlongLine( static_cast< const QgsLineString * >( polygon->exteriorRing() ),
2146  wavelength, amplitude, strictWavelength ).release() );
2147  for ( int i = 0; i < polygon->numInteriorRings(); ++i )
2148  {
2149  result->addInteriorRing( squareWavesAlongLine( static_cast< const QgsLineString * >( polygon->interiorRing( i ) ),
2150  wavelength, amplitude, strictWavelength ).release() );
2151  }
2152 
2153  return result;
2154  }
2155 }
2156 
2157 std::unique_ptr< QgsAbstractGeometry > squareWavesRandomizedPrivate( const QgsAbstractGeometry *geom, double minimumWavelength, double maximumWavelength, double minimumAmplitude, double maximumAmplitude, std::uniform_real_distribution<> &uniformDist, std::mt19937 &mt )
2158 {
2159  std::unique_ptr< QgsAbstractGeometry > segmentizedCopy;
2160  if ( QgsWkbTypes::isCurvedType( geom->wkbType() ) )
2161  {
2162  segmentizedCopy.reset( geom->segmentize() );
2163  geom = segmentizedCopy.get();
2164  }
2165 
2167  {
2168  return squareWavesRandomizedAlongLine( static_cast< const QgsLineString * >( geom ), minimumWavelength, maximumWavelength, minimumAmplitude, maximumAmplitude, uniformDist, mt );
2169  }
2170  else
2171  {
2172  // polygon
2173  const QgsPolygon *polygon = static_cast< const QgsPolygon * >( geom );
2174  std::unique_ptr< QgsPolygon > result = std::make_unique< QgsPolygon >();
2175 
2176  result->setExteriorRing( squareWavesRandomizedAlongLine( static_cast< const QgsLineString * >( polygon->exteriorRing() ),
2177  minimumWavelength, maximumWavelength, minimumAmplitude, maximumAmplitude, uniformDist, mt ).release() );
2178  for ( int i = 0; i < polygon->numInteriorRings(); ++i )
2179  {
2180  result->addInteriorRing( squareWavesRandomizedAlongLine( static_cast< const QgsLineString * >( polygon->interiorRing( i ) ),
2181  minimumWavelength, maximumWavelength, minimumAmplitude, maximumAmplitude, uniformDist, mt ).release() );
2182  }
2183 
2184  return result;
2185  }
2186 }
2187 
2188 QgsGeometry QgsInternalGeometryEngine::squareWaves( double wavelength, double amplitude, bool strictWavelength ) const
2189 {
2190  if ( wavelength < 0 || qgsDoubleNear( wavelength, 0 ) )
2191  return QgsGeometry();
2192 
2193  mLastError.clear();
2194  if ( !mGeometry )
2195  {
2196  return QgsGeometry();
2197  }
2198 
2200  {
2201  return QgsGeometry( mGeometry->clone() ); // point geometry, nothing to do
2202  }
2203 
2204  if ( const QgsGeometryCollection *gc = qgsgeometry_cast< const QgsGeometryCollection *>( mGeometry ) )
2205  {
2206  int numGeom = gc->numGeometries();
2207  QVector< QgsAbstractGeometry * > geometryList;
2208  geometryList.reserve( numGeom );
2209  for ( int i = 0; i < numGeom; ++i )
2210  {
2211  geometryList << squareWavesPrivate( gc->geometryN( i ), wavelength, amplitude, strictWavelength ).release();
2212  }
2213 
2214  QgsGeometry first = QgsGeometry( geometryList.takeAt( 0 ) );
2215  for ( QgsAbstractGeometry *g : std::as_const( geometryList ) )
2216  {
2217  first.addPart( g );
2218  }
2219  return first;
2220  }
2221  else
2222  {
2223  return QgsGeometry( squareWavesPrivate( mGeometry, wavelength, amplitude, strictWavelength ) );
2224  }
2225 }
2226 
2227 QgsGeometry QgsInternalGeometryEngine::squareWavesRandomized( double minimumWavelength, double maximumWavelength, double minimumAmplitude, double maximumAmplitude, unsigned long seed ) const
2228 {
2229  if ( minimumWavelength < 0 || qgsDoubleNear( minimumWavelength, 0 ) || maximumWavelength < 0 || qgsDoubleNear( maximumWavelength, 0 ) || maximumWavelength < minimumWavelength )
2230  return QgsGeometry();
2231 
2232  mLastError.clear();
2233  if ( !mGeometry )
2234  {
2235  return QgsGeometry();
2236  }
2237 
2238  std::random_device rd;
2239  std::mt19937 mt( seed == 0 ? rd() : seed );
2240  std::uniform_real_distribution<> uniformDist( 0, 1 );
2241 
2243  {
2244  return QgsGeometry( mGeometry->clone() ); // point geometry, nothing to do
2245  }
2246 
2247  if ( const QgsGeometryCollection *gc = qgsgeometry_cast< const QgsGeometryCollection *>( mGeometry ) )
2248  {
2249  int numGeom = gc->numGeometries();
2250  QVector< QgsAbstractGeometry * > geometryList;
2251  geometryList.reserve( numGeom );
2252  for ( int i = 0; i < numGeom; ++i )
2253  {
2254  geometryList << squareWavesRandomizedPrivate( gc->geometryN( i ), minimumWavelength, maximumWavelength, minimumAmplitude, maximumAmplitude, uniformDist, mt ).release();
2255  }
2256 
2257  QgsGeometry first = QgsGeometry( geometryList.takeAt( 0 ) );
2258  for ( QgsAbstractGeometry *g : std::as_const( geometryList ) )
2259  {
2260  first.addPart( g );
2261  }
2262  return first;
2263  }
2264  else
2265  {
2266  return QgsGeometry( squareWavesRandomizedPrivate( mGeometry, minimumWavelength, maximumWavelength, minimumAmplitude, maximumAmplitude, uniformDist, mt ) );
2267  }
2268 }
2269 
2270 std::unique_ptr< QgsLineString > roundWavesAlongLine( const QgsLineString *line, double wavelength, const double amplitude, const bool strictWavelength )
2271 {
2272  const int totalPoints = line->numPoints();
2273  if ( totalPoints < 2 )
2274  return nullptr;
2275 
2276  const double *x = line->xData();
2277  const double *y = line->yData();
2278 
2279  double prevX = *x++;
2280  double prevY = *y++;
2281 
2282  const double totalLength = line->length();
2283 
2284  const int maxRepetitions = std::ceil( totalLength / wavelength );
2285  if ( !strictWavelength )
2286  wavelength = totalLength / maxRepetitions;
2287 
2288  const int segments = 10;
2289 
2290  int side = -1;
2291 
2292  double xOutBuffer[4] { prevX, prevX, prevX, prevX };
2293  double yOutBuffer[4] { prevY, prevY, prevY, prevY };
2294  bool isFirstPart = true;
2295 
2296  double distanceToNextPointFromStartOfSegment = wavelength / 8;
2297  int bufferIndex = 1;
2298  std::unique_ptr< QgsLineString > out = std::make_unique< QgsLineString >();
2299 
2300  double segmentAngleRadians = 0;
2301  double remainingDistance = totalLength;
2302  double totalCoveredDistance = 0;
2303 
2304  for ( int i = 1; i < totalPoints; ++i )
2305  {
2306  double thisX = *x++;
2307  double thisY = *y++;
2308 
2309  segmentAngleRadians = QgsGeometryUtils::lineAngle( prevX, prevY, thisX, thisY );
2310  const double segmentLength = std::sqrt( ( thisX - prevX ) * ( thisX - prevX ) + ( thisY - prevY ) * ( thisY - prevY ) );
2311  while ( distanceToNextPointFromStartOfSegment < segmentLength || qgsDoubleNear( distanceToNextPointFromStartOfSegment, segmentLength ) )
2312  {
2313  // point falls on this segment - truncate to segment length if qgsDoubleNear test was actually > segment length
2314  const double distanceToPoint = std::min( distanceToNextPointFromStartOfSegment, segmentLength );
2315  double pX, pY;
2316  QgsGeometryUtils::pointOnLineWithDistance( prevX, prevY, thisX, thisY, distanceToPoint, pX, pY );
2317  remainingDistance = totalLength - totalCoveredDistance - distanceToPoint;
2318 
2319  const double sinAngle = std::sin( segmentAngleRadians + M_PI_2 );
2320  const double cosAngle = std::cos( segmentAngleRadians + M_PI_2 );
2321 
2322  if ( bufferIndex == 0 )
2323  {
2324  xOutBuffer[0] = pX + side * amplitude * sinAngle;
2325  yOutBuffer[0] = pY + side * amplitude * cosAngle;
2326  bufferIndex = 1;
2327  distanceToNextPointFromStartOfSegment += wavelength / 4;
2328  }
2329  else if ( bufferIndex == 1 && isFirstPart )
2330  {
2331  xOutBuffer[1] = ( xOutBuffer[0] + pX - side * amplitude * sinAngle ) * 0.5;
2332  yOutBuffer[1] = ( yOutBuffer[0] + pY - side * amplitude * cosAngle ) * 0.5;
2333  xOutBuffer[2] = pX - side * amplitude * sinAngle;
2334  yOutBuffer[2] = pY - side * amplitude * cosAngle;
2335  bufferIndex = 2;
2336  distanceToNextPointFromStartOfSegment += wavelength / 8;
2337  }
2338  else if ( bufferIndex == 1 )
2339  {
2340  xOutBuffer[1] = pX + side * amplitude * sinAngle;
2341  yOutBuffer[1] = pY + side * amplitude * cosAngle;
2342  xOutBuffer[2] = pX - side * amplitude * sinAngle;
2343  yOutBuffer[2] = pY - side * amplitude * cosAngle;
2344  bufferIndex = 2;
2345  distanceToNextPointFromStartOfSegment += wavelength / 4;
2346  }
2347  else if ( bufferIndex == 2 )
2348  {
2349  xOutBuffer[3] = pX - side * amplitude * sinAngle;
2350  yOutBuffer[3] = pY - side * amplitude * cosAngle;
2351 
2352  if ( isFirstPart )
2353  {
2354  xOutBuffer[2] = ( xOutBuffer[2] + xOutBuffer[3] ) * 0.5;
2355  yOutBuffer[2] = ( yOutBuffer[2] + yOutBuffer[3] ) * 0.5;
2356  isFirstPart = false;
2357  }
2358 
2359  // flush curve
2360  std::unique_ptr< QgsLineString > bezier( QgsLineString::fromBezierCurve( QgsPoint( xOutBuffer[0], yOutBuffer[0] ),
2361  QgsPoint( xOutBuffer[1], yOutBuffer[1] ),
2362  QgsPoint( xOutBuffer[2], yOutBuffer[2] ),
2363  QgsPoint( xOutBuffer[3], yOutBuffer[3] ),
2364  segments ) );
2365  out->append( bezier.get() );
2366 
2367  // shuffle buffer alone
2368  xOutBuffer[0] = xOutBuffer[3];
2369  yOutBuffer[0] = yOutBuffer[3];
2370  bufferIndex = 1;
2371  side = -side;
2372  distanceToNextPointFromStartOfSegment += wavelength / 4;
2373  }
2374  }
2375  totalCoveredDistance += segmentLength;
2376  prevX = thisX;
2377  prevY = thisY;
2378  distanceToNextPointFromStartOfSegment -= segmentLength;
2379  }
2380 
2381  const double midX = prevX - remainingDistance / 2 * std::sin( segmentAngleRadians );
2382  const double midY = prevY - remainingDistance / 2 * std::cos( segmentAngleRadians );
2383  const double pX = midX + side * amplitude * std::sin( segmentAngleRadians + M_PI_2 );
2384  const double pY = midY + side * amplitude * std::cos( segmentAngleRadians + M_PI_2 );
2385 
2386  std::unique_ptr< QgsLineString > bezier( QgsLineString::fromBezierCurve( out->endPoint(),
2387  QgsPoint( ( out->endPoint().x() + pX ) * 0.5, ( out->endPoint().y() + pY ) * 0.5 ),
2388  QgsPoint( ( prevX + pX ) * 0.5, ( prevY + pY ) * 0.5 ),
2389  QgsPoint( prevX, prevY ),
2390  segments / 2 ) );
2391  out->append( bezier.get() );
2392 
2393  return out;
2394 }
2395 
2396 std::unique_ptr< QgsLineString > roundWavesRandomizedAlongLine( const QgsLineString *line,
2397  const double minimumWavelength, const double maximumWavelength,
2398  const double minimumAmplitude, const double maximumAmplitude,
2399  std::uniform_real_distribution<> &uniformDist, std::mt19937 &mt )
2400 {
2401  const int totalPoints = line->numPoints();
2402  if ( totalPoints < 2 )
2403  return nullptr;
2404 
2405  const double *x = line->xData();
2406  const double *y = line->yData();
2407 
2408  double prevX = *x++;
2409  double prevY = *y++;
2410 
2411  const double totalLength = line->length();
2412 
2413  const int segments = 10;
2414 
2415  int side = -1;
2416 
2417  double xOutBuffer[4] { prevX, prevX, prevX, prevX };
2418  double yOutBuffer[4] { prevY, prevY, prevY, prevY };
2419  bool isFirstPart = true;
2420 
2421  double amplitude = uniformDist( mt ) * ( maximumAmplitude - minimumAmplitude ) + minimumAmplitude;
2422  double wavelength = uniformDist( mt ) * ( maximumWavelength - minimumWavelength ) + minimumWavelength;
2423 
2424  double distanceToNextPointFromStartOfSegment = wavelength / 8;
2425  int bufferIndex = 1;
2426  std::unique_ptr< QgsLineString > out = std::make_unique< QgsLineString >();
2427 
2428  double segmentAngleRadians = 0;
2429 
2430  double remainingDistance = totalLength;
2431  double totalCoveredDistance = 0;
2432 
2433  for ( int i = 1; i < totalPoints; ++i )
2434  {
2435  double thisX = *x++;
2436  double thisY = *y++;
2437 
2438  segmentAngleRadians = QgsGeometryUtils::lineAngle( prevX, prevY, thisX, thisY );
2439  const double segmentLength = std::sqrt( ( thisX - prevX ) * ( thisX - prevX ) + ( thisY - prevY ) * ( thisY - prevY ) );
2440  while ( distanceToNextPointFromStartOfSegment < segmentLength || qgsDoubleNear( distanceToNextPointFromStartOfSegment, segmentLength ) )
2441  {
2442  // point falls on this segment - truncate to segment length if qgsDoubleNear test was actually > segment length
2443  const double distanceToPoint = std::min( distanceToNextPointFromStartOfSegment, segmentLength );
2444  double pX, pY;
2445  QgsGeometryUtils::pointOnLineWithDistance( prevX, prevY, thisX, thisY, distanceToPoint, pX, pY );
2446  remainingDistance = totalLength - totalCoveredDistance - distanceToPoint;
2447 
2448  const double sinAngle = std::sin( segmentAngleRadians + M_PI_2 );
2449  const double cosAngle = std::cos( segmentAngleRadians + M_PI_2 );
2450 
2451  if ( bufferIndex == 0 )
2452  {
2453  xOutBuffer[0] = pX + side * amplitude * sinAngle;
2454  yOutBuffer[0] = pY + side * amplitude * cosAngle;
2455  bufferIndex = 1;
2456  distanceToNextPointFromStartOfSegment += wavelength / 4;
2457  }
2458  else if ( bufferIndex == 1 && isFirstPart )
2459  {
2460  xOutBuffer[1] = ( xOutBuffer[0] + pX - side * amplitude * sinAngle ) * 0.5;
2461  yOutBuffer[1] = ( yOutBuffer[0] + pY - side * amplitude * cosAngle ) * 0.5;
2462  xOutBuffer[2] = pX - side * amplitude * sinAngle;
2463  yOutBuffer[2] = pY - side * amplitude * cosAngle;
2464  bufferIndex = 2;
2465  distanceToNextPointFromStartOfSegment += wavelength / 8;
2466  }
2467  else if ( bufferIndex == 1 )
2468  {
2469  xOutBuffer[1] = pX + side * amplitude * sinAngle;
2470  yOutBuffer[1] = pY + side * amplitude * cosAngle;
2471  amplitude = uniformDist( mt ) * ( maximumAmplitude - minimumAmplitude ) + minimumAmplitude;
2472  xOutBuffer[2] = pX - side * amplitude * sinAngle;
2473  bufferIndex = 2;
2474  yOutBuffer[2] = pY - side * amplitude * cosAngle;
2475  distanceToNextPointFromStartOfSegment += wavelength / 4;
2476  }
2477  else if ( bufferIndex == 2 )
2478  {
2479  xOutBuffer[3] = pX - side * amplitude * sinAngle;
2480  yOutBuffer[3] = pY - side * amplitude * cosAngle;
2481 
2482  if ( isFirstPart )
2483  {
2484  xOutBuffer[2] = ( xOutBuffer[2] + xOutBuffer[3] ) * 0.5;
2485  yOutBuffer[2] = ( yOutBuffer[2] + yOutBuffer[3] ) * 0.5;
2486  isFirstPart = false;
2487  }
2488 
2489  // flush curve
2490  std::unique_ptr< QgsLineString > bezier( QgsLineString::fromBezierCurve( QgsPoint( xOutBuffer[0], yOutBuffer[0] ),
2491  QgsPoint( xOutBuffer[1], yOutBuffer[1] ),
2492  QgsPoint( xOutBuffer[2], yOutBuffer[2] ),
2493  QgsPoint( xOutBuffer[3], yOutBuffer[3] ),
2494  segments ) );
2495  out->append( bezier.get() );
2496 
2497  // shuffle buffer alone
2498  xOutBuffer[0] = xOutBuffer[3];
2499  yOutBuffer[0] = yOutBuffer[3];
2500  bufferIndex = 1;
2501  side = -side;
2502 
2503  wavelength = uniformDist( mt ) * ( maximumWavelength - minimumWavelength ) + minimumWavelength;
2504 
2505  distanceToNextPointFromStartOfSegment += wavelength / 4;
2506  }
2507  }
2508  totalCoveredDistance += segmentLength;
2509 
2510  prevX = thisX;
2511  prevY = thisY;
2512  distanceToNextPointFromStartOfSegment -= segmentLength;
2513  }
2514 
2515  const double midX = prevX - remainingDistance / 2 * std::sin( segmentAngleRadians );
2516  const double midY = prevY - remainingDistance / 2 * std::cos( segmentAngleRadians );
2517  const double pX = midX + side * amplitude * std::sin( segmentAngleRadians + M_PI_2 );
2518  const double pY = midY + side * amplitude * std::cos( segmentAngleRadians + M_PI_2 );
2519 
2520  if ( out->isEmpty() )
2521  {
2522  std::unique_ptr< QgsLineString > bezier( QgsLineString::fromBezierCurve( line->startPoint(),
2523  QgsPoint( ( line->startPoint().x() + pX ) * 0.5, ( line->startPoint().y() + pY ) * 0.5 ),
2524  QgsPoint( ( prevX + pX ) * 0.5, ( prevY + pY ) * 0.5 ),
2525  QgsPoint( prevX, prevY ),
2526  segments ) );
2527  out->append( bezier.get() );
2528  }
2529  else
2530  {
2531  std::unique_ptr< QgsLineString > bezier( QgsLineString::fromBezierCurve( out->endPoint(),
2532  QgsPoint( ( out->endPoint().x() + pX ) * 0.5, ( out->endPoint().y() + pY ) * 0.5 ),
2533  QgsPoint( ( prevX + pX ) * 0.5, ( prevY + pY ) * 0.5 ),
2534  QgsPoint( prevX, prevY ),
2535  segments ) );
2536  out->append( bezier.get() );
2537  }
2538 
2539  return out;
2540 }
2541 
2542 std::unique_ptr< QgsAbstractGeometry > roundWavesPrivate( const QgsAbstractGeometry *geom, double wavelength, double amplitude, bool strictWavelength )
2543 {
2544  std::unique_ptr< QgsAbstractGeometry > segmentizedCopy;
2545  if ( QgsWkbTypes::isCurvedType( geom->wkbType() ) )
2546  {
2547  segmentizedCopy.reset( geom->segmentize() );
2548  geom = segmentizedCopy.get();
2549  }
2550 
2552  {
2553  return roundWavesAlongLine( static_cast< const QgsLineString * >( geom ), wavelength, amplitude, strictWavelength );
2554  }
2555  else
2556  {
2557  // polygon
2558  const QgsPolygon *polygon = static_cast< const QgsPolygon * >( geom );
2559  std::unique_ptr< QgsPolygon > result = std::make_unique< QgsPolygon >();
2560 
2561  result->setExteriorRing( roundWavesAlongLine( static_cast< const QgsLineString * >( polygon->exteriorRing() ),
2562  wavelength, amplitude, strictWavelength ).release() );
2563  for ( int i = 0; i < polygon->numInteriorRings(); ++i )
2564  {
2565  result->addInteriorRing( roundWavesAlongLine( static_cast< const QgsLineString * >( polygon->interiorRing( i ) ),
2566  wavelength, amplitude, strictWavelength ).release() );
2567  }
2568 
2569  return result;
2570  }
2571 }
2572 
2573 std::unique_ptr< QgsAbstractGeometry > roundWavesRandomizedPrivate( const QgsAbstractGeometry *geom, double minimumWavelength, double maximumWavelength, double minimumAmplitude, double maximumAmplitude, std::uniform_real_distribution<> &uniformDist, std::mt19937 &mt )
2574 {
2575  std::unique_ptr< QgsAbstractGeometry > segmentizedCopy;
2576  if ( QgsWkbTypes::isCurvedType( geom->wkbType() ) )
2577  {
2578  segmentizedCopy.reset( geom->segmentize() );
2579  geom = segmentizedCopy.get();
2580  }
2581 
2583  {
2584  return roundWavesRandomizedAlongLine( static_cast< const QgsLineString * >( geom ), minimumWavelength, maximumWavelength, minimumAmplitude, maximumAmplitude, uniformDist, mt );
2585  }
2586  else
2587  {
2588  // polygon
2589  const QgsPolygon *polygon = static_cast< const QgsPolygon * >( geom );
2590  std::unique_ptr< QgsPolygon > result = std::make_unique< QgsPolygon >();
2591 
2592  result->setExteriorRing( roundWavesRandomizedAlongLine( static_cast< const QgsLineString * >( polygon->exteriorRing() ),
2593  minimumWavelength, maximumWavelength, minimumAmplitude, maximumAmplitude, uniformDist, mt ).release() );
2594  for ( int i = 0; i < polygon->numInteriorRings(); ++i )
2595  {
2596  result->addInteriorRing( roundWavesRandomizedAlongLine( static_cast< const QgsLineString * >( polygon->interiorRing( i ) ),
2597  minimumWavelength, maximumWavelength, minimumAmplitude, maximumAmplitude, uniformDist, mt ).release() );
2598  }
2599 
2600  return result;
2601  }
2602 }
2603 
2604 QgsGeometry QgsInternalGeometryEngine::roundWaves( double wavelength, double amplitude, bool strictWavelength ) const
2605 {
2606  if ( wavelength < 0 || qgsDoubleNear( wavelength, 0 ) )
2607  return QgsGeometry();
2608 
2609  mLastError.clear();
2610  if ( !mGeometry )
2611  {
2612  return QgsGeometry();
2613  }
2614 
2616  {
2617  return QgsGeometry( mGeometry->clone() ); // point geometry, nothing to do
2618  }
2619 
2620  if ( const QgsGeometryCollection *gc = qgsgeometry_cast< const QgsGeometryCollection *>( mGeometry ) )
2621  {
2622  int numGeom = gc->numGeometries();
2623  QVector< QgsAbstractGeometry * > geometryList;
2624  geometryList.reserve( numGeom );
2625  for ( int i = 0; i < numGeom; ++i )
2626  {
2627  geometryList << roundWavesPrivate( gc->geometryN( i ), wavelength, amplitude, strictWavelength ).release();
2628  }
2629 
2630  QgsGeometry first = QgsGeometry( geometryList.takeAt( 0 ) );
2631  for ( QgsAbstractGeometry *g : std::as_const( geometryList ) )
2632  {
2633  first.addPart( g );
2634  }
2635  return first;
2636  }
2637  else
2638  {
2639  return QgsGeometry( roundWavesPrivate( mGeometry, wavelength, amplitude, strictWavelength ) );
2640  }
2641 }
2642 
2643 QgsGeometry QgsInternalGeometryEngine::roundWavesRandomized( double minimumWavelength, double maximumWavelength, double minimumAmplitude, double maximumAmplitude, unsigned long seed ) const
2644 {
2645  if ( minimumWavelength < 0 || qgsDoubleNear( minimumWavelength, 0 ) || maximumWavelength < 0 || qgsDoubleNear( maximumWavelength, 0 ) || maximumWavelength < minimumWavelength )
2646  return QgsGeometry();
2647 
2648  mLastError.clear();
2649  if ( !mGeometry )
2650  {
2651  return QgsGeometry();
2652  }
2653 
2654  std::random_device rd;
2655  std::mt19937 mt( seed == 0 ? rd() : seed );
2656  std::uniform_real_distribution<> uniformDist( 0, 1 );
2657 
2659  {
2660  return QgsGeometry( mGeometry->clone() ); // point geometry, nothing to do
2661  }
2662 
2663  if ( const QgsGeometryCollection *gc = qgsgeometry_cast< const QgsGeometryCollection *>( mGeometry ) )
2664  {
2665  int numGeom = gc->numGeometries();
2666  QVector< QgsAbstractGeometry * > geometryList;
2667  geometryList.reserve( numGeom );
2668  for ( int i = 0; i < numGeom; ++i )
2669  {
2670  geometryList << roundWavesRandomizedPrivate( gc->geometryN( i ), minimumWavelength, maximumWavelength, minimumAmplitude, maximumAmplitude, uniformDist, mt ).release();
2671  }
2672 
2673  QgsGeometry first = QgsGeometry( geometryList.takeAt( 0 ) );
2674  for ( QgsAbstractGeometry *g : std::as_const( geometryList ) )
2675  {
2676  first.addPart( g );
2677  }
2678  return first;
2679  }
2680  else
2681  {
2682  return QgsGeometry( roundWavesRandomizedPrivate( mGeometry, minimumWavelength, maximumWavelength, minimumAmplitude, maximumAmplitude, uniformDist, mt ) );
2683  }
2684 }
2685 
2686 std::unique_ptr< QgsMultiLineString > dashPatternAlongLine( const QgsLineString *line,
2687  const QVector< double> &pattern,
2691  double patternOffset )
2692 {
2693  const int totalPoints = line->numPoints();
2694  if ( totalPoints < 2 )
2695  return nullptr;
2696 
2697  const int patternSize = pattern.size();
2698 
2699  const double *x = line->xData();
2700  const double *y = line->yData();
2701 
2702  double prevX = *x++;
2703  double prevY = *y++;
2704 
2705  std::unique_ptr< QgsMultiLineString > result = std::make_unique< QgsMultiLineString >();
2706 
2707  QVector< double > outX;
2708  QVector< double > outY;
2709  const double totalLength = line->length();
2710 
2711  double patternLength = 0;
2712  double patternDashLength = 0;
2713  double patternGapLength = 0;
2714  for ( int i = 0; i < pattern.size(); ++i )
2715  {
2716  patternLength += pattern.at( i );
2717  if ( i % 2 == 0 )
2718  patternDashLength += pattern.at( i );
2719  else
2720  patternGapLength += pattern.at( i );
2721  }
2722 
2723  double firstPatternLength = 0;
2724  double firstPatternDashLength = 0;
2725  double firstPatternGapLength = 0;
2726  switch ( startRule )
2727  {
2730  firstPatternLength = patternLength;
2731  firstPatternDashLength = patternDashLength;
2732  firstPatternGapLength = patternGapLength;
2733  break;
2735  firstPatternLength = patternLength - pattern.at( 0 ) * 0.5; // remove half of first dash
2736  firstPatternDashLength = patternDashLength - pattern.at( 0 ) * 0.5;
2737  firstPatternGapLength = patternGapLength;
2738  break;
2740  firstPatternLength = pattern.at( patternSize - 1 ); // start at LAST dash
2741  firstPatternDashLength = 0;
2742  firstPatternGapLength = pattern.at( patternSize - 1 );
2743  break;
2745  firstPatternLength = pattern.at( patternSize - 1 ) * 0.5; // start at half of last dash
2746  firstPatternDashLength = 0;
2747  firstPatternGapLength = pattern.at( patternSize - 1 ) * 0.5;
2748  break;
2749  }
2750 
2751  const bool isSmallEnoughForSinglePattern = ( totalLength - firstPatternLength ) < patternLength * 0.5;
2752 
2753  double lastPatternLength = isSmallEnoughForSinglePattern ? firstPatternLength : patternLength;
2754  double lastPatternDashLength = isSmallEnoughForSinglePattern ? firstPatternDashLength : patternDashLength;
2755  double lastPatternGapLength = isSmallEnoughForSinglePattern ? firstPatternGapLength : patternGapLength;
2756  switch ( endRule )
2757  {
2759  lastPatternLength = 0;
2760  lastPatternDashLength = 0;
2761  lastPatternGapLength = 0;
2762  break;
2764  lastPatternLength -= pattern.at( patternSize - 1 ); // remove last gap
2765  lastPatternGapLength -= pattern.at( patternSize - 1 );
2766  break;
2768  lastPatternLength -= pattern.at( patternSize - 1 ) + pattern.at( patternSize - 2 ) * 0.5; // remove last gap, half of last dash
2769  lastPatternDashLength -= pattern.at( patternSize - 2 ) * 0.5;
2770  lastPatternGapLength -= pattern.at( patternSize - 1 );
2771  break;
2773  lastPatternGapLength = patternGapLength;
2774  break;
2776  lastPatternLength -= pattern.at( patternSize - 1 ) * 0.5; // remove half of last gap
2777  lastPatternGapLength -= pattern.at( patternSize - 1 ) * 0.5;
2778  break;
2779  }
2780 
2781  const double remainingLengthForCompletePatterns = totalLength - ( !isSmallEnoughForSinglePattern ? firstPatternLength : 0 ) - lastPatternLength;
2782  const int middlePatternRepetitions = std::max( static_cast< int >( std::round( remainingLengthForCompletePatterns / patternLength ) ), 0 );
2783 
2784  const double totalUnscaledLengthOfPatterns = ( !isSmallEnoughForSinglePattern ? firstPatternLength : 0 ) + middlePatternRepetitions * patternLength + lastPatternLength;
2785  const double totalUnscaledLengthOfDashes = ( !isSmallEnoughForSinglePattern ? firstPatternDashLength : 0 ) + middlePatternRepetitions * patternDashLength + lastPatternDashLength;
2786  const double totalUnscaledLengthOfGaps = ( !isSmallEnoughForSinglePattern ? firstPatternGapLength : 0 ) + middlePatternRepetitions * patternGapLength + lastPatternGapLength;
2787 
2788  double dashLengthScalingFactor = 1;
2789  double gapLengthScalingFactor = 1;
2790  if ( endRule != Qgis::DashPatternLineEndingRule::NoRule )
2791  {
2792  // calculate scaling factors
2793  const double lengthToShrinkBy = totalUnscaledLengthOfPatterns - totalLength;
2794 
2795  switch ( adjustment )
2796  {
2798  dashLengthScalingFactor = totalLength / totalUnscaledLengthOfPatterns;
2799  gapLengthScalingFactor = dashLengthScalingFactor;
2800  break;
2802  dashLengthScalingFactor = ( totalUnscaledLengthOfDashes - lengthToShrinkBy ) / totalUnscaledLengthOfDashes;
2803  break;
2805  gapLengthScalingFactor = ( totalUnscaledLengthOfGaps - lengthToShrinkBy ) / totalUnscaledLengthOfGaps;
2806  break;
2807  }
2808  }
2809 
2810  dashLengthScalingFactor = std::max( dashLengthScalingFactor, 0.0 );
2811  gapLengthScalingFactor = std::max( gapLengthScalingFactor, 0.0 );
2812 
2813  const int maxPatterns = middlePatternRepetitions + 2;
2814  result->reserve( maxPatterns );
2815 
2816  int patternIndex = 0;
2817  double distanceToNextPointFromStartOfSegment = pattern.at( 0 ) * dashLengthScalingFactor;
2818  bool isDash = true;
2819  switch ( startRule )
2820  {
2823  break;
2825  distanceToNextPointFromStartOfSegment *= 0.5; // skip to half way through first dash
2826  break;
2828  patternIndex = patternSize - 1; // start at last gap
2829  isDash = false;
2830  distanceToNextPointFromStartOfSegment = pattern.at( patternSize - 1 ) * gapLengthScalingFactor;
2831  break;
2833  patternIndex = patternSize - 1; // skip straight to half way through last gap
2834  isDash = false;
2835  distanceToNextPointFromStartOfSegment = 0.5 * pattern.at( patternSize - 1 ) * gapLengthScalingFactor;
2836  break;
2837  }
2838 
2839  const double adjustedOffset = fmod( patternOffset, patternLength );
2840  const double scaledOffset = ( adjustedOffset < 0 ? ( adjustedOffset + patternLength ) : adjustedOffset ) * ( gapLengthScalingFactor + dashLengthScalingFactor ) / 2;
2841  if ( !qgsDoubleNear( scaledOffset, 0 ) )
2842  {
2843  // shuffle pattern along by offset
2844  double remainingOffset = scaledOffset;
2845  while ( remainingOffset > 0 )
2846  {
2847  if ( distanceToNextPointFromStartOfSegment > remainingOffset )
2848  {
2849  distanceToNextPointFromStartOfSegment -= remainingOffset;
2850  break;
2851  }
2852 
2853  remainingOffset -= distanceToNextPointFromStartOfSegment;
2854  isDash = !isDash;
2855  patternIndex++;
2856  if ( patternIndex == patternSize )
2857  patternIndex = 0;
2858 
2859  distanceToNextPointFromStartOfSegment = pattern.at( patternIndex ) * ( isDash ? dashLengthScalingFactor : gapLengthScalingFactor );
2860  }
2861  }
2862 
2863  if ( isDash )
2864  {
2865  outX.append( prevX );
2866  outY.append( prevY );
2867  }
2868 
2869  for ( int i = 1; i < totalPoints; ++i )
2870  {
2871  double thisX = *x++;
2872  double thisY = *y++;
2873 
2874  const double segmentLength = std::sqrt( ( thisX - prevX ) * ( thisX - prevX ) + ( thisY - prevY ) * ( thisY - prevY ) );
2875  while ( distanceToNextPointFromStartOfSegment < segmentLength || qgsDoubleNear( distanceToNextPointFromStartOfSegment, segmentLength ) )
2876  {
2877  // point falls on this segment - truncate to segment length if qgsDoubleNear test was actually > segment length
2878  const double distanceToPoint = std::min( distanceToNextPointFromStartOfSegment, segmentLength );
2879  double pX, pY;
2880  QgsGeometryUtils::pointOnLineWithDistance( prevX, prevY, thisX, thisY, distanceToPoint, pX, pY );
2881 
2882  outX.append( pX );
2883  outY.append( pY );
2884  if ( isDash )
2885  {
2886  result->addGeometry( new QgsLineString( outX, outY ) );
2887  outX.resize( 0 );
2888  outY.resize( 0 );
2889  }
2890 
2891  isDash = !isDash;
2892  patternIndex++;
2893  if ( patternIndex >= patternSize )
2894  patternIndex = 0;
2895 
2896  distanceToNextPointFromStartOfSegment += pattern.at( patternIndex ) * ( isDash ? dashLengthScalingFactor : gapLengthScalingFactor );
2897  }
2898 
2899  if ( isDash )
2900  {
2901  outX.append( thisX );
2902  outY.append( thisY );
2903  }
2904 
2905  prevX = thisX;
2906  prevY = thisY;
2907  distanceToNextPointFromStartOfSegment -= segmentLength;
2908  }
2909 
2910  if ( isDash )
2911  {
2912  outX.append( prevX );
2913  outY.append( prevY );
2914  result->addGeometry( new QgsLineString( outX, outY ) );
2915  }
2916 
2917  return result;
2918 }
2919 
2920 std::unique_ptr< QgsAbstractGeometry > applyDashPatternPrivate( const QgsAbstractGeometry *geom,
2921  const QVector<double> &pattern,
2925  double patternOffset )
2926 {
2927  std::unique_ptr< QgsAbstractGeometry > segmentizedCopy;
2928  if ( QgsWkbTypes::isCurvedType( geom->wkbType() ) )
2929  {
2930  segmentizedCopy.reset( geom->segmentize() );
2931  geom = segmentizedCopy.get();
2932  }
2933 
2935  {
2936  return dashPatternAlongLine( static_cast< const QgsLineString * >( geom ), pattern, startRule, endRule, adjustment, patternOffset );
2937  }
2938  else
2939  {
2940  // polygon
2941  const QgsPolygon *polygon = static_cast< const QgsPolygon * >( geom );
2942  std::unique_ptr< QgsMultiLineString > result = std::make_unique< QgsMultiLineString >();
2943 
2944  std::unique_ptr< QgsMultiLineString > exteriorParts = dashPatternAlongLine( static_cast< const QgsLineString * >( polygon->exteriorRing() ), pattern, startRule, endRule, adjustment, patternOffset );
2945  for ( int i = 0; i < exteriorParts->numGeometries(); ++i )
2946  result->addGeometry( exteriorParts->geometryN( i )->clone() );
2947 
2948  for ( int i = 0; i < polygon->numInteriorRings(); ++i )
2949  {
2950  std::unique_ptr< QgsMultiLineString > ringParts = dashPatternAlongLine( static_cast< const QgsLineString * >( polygon->interiorRing( i ) ), pattern, startRule, endRule, adjustment, patternOffset );
2951  for ( int j = 0; j < ringParts->numGeometries(); ++j )
2952  result->addGeometry( ringParts->geometryN( j )->clone() );
2953  }
2954 
2955  return result;
2956  }
2957 }
2958 
2960 {
2961  if ( pattern.size() < 2 )
2962  return QgsGeometry( mGeometry->clone() );
2963 
2964  mLastError.clear();
2965  if ( !mGeometry || mGeometry->isEmpty() )
2966  {
2967  return QgsGeometry();
2968  }
2969 
2971  {
2972  return QgsGeometry( mGeometry->clone() ); // point geometry, nothing to do
2973  }
2974 
2975  if ( const QgsGeometryCollection *gc = qgsgeometry_cast< const QgsGeometryCollection *>( mGeometry ) )
2976  {
2977  int numGeom = gc->numGeometries();
2978  QVector< QgsAbstractGeometry * > geometryList;
2979  geometryList.reserve( numGeom );
2980  for ( int i = 0; i < numGeom; ++i )
2981  {
2982  geometryList << applyDashPatternPrivate( gc->geometryN( i ), pattern, startRule, endRule, adjustment, patternOffset ).release();
2983  }
2984 
2985  QgsGeometry first = QgsGeometry( geometryList.takeAt( 0 ) );
2986  for ( QgsAbstractGeometry *g : std::as_const( geometryList ) )
2987  {
2988  if ( const QgsGeometryCollection *collection = qgsgeometry_cast< QgsGeometryCollection * >( g ) )
2989  {
2990  for ( int j = 0; j < collection->numGeometries(); ++j )
2991  {
2992  first.addPart( collection->geometryN( j )->clone() );
2993  }
2994  delete collection;
2995  }
2996  else
2997  {
2998  first.addPart( g );
2999  }
3000  }
3001  return first;
3002  }
3003  else
3004  {
3005  return QgsGeometry( applyDashPatternPrivate( mGeometry, pattern, startRule, endRule, adjustment, patternOffset ) );
3006  }
3007 }
DashPatternSizeAdjustment
Dash pattern size adjustment options.
Definition: qgis.h:1565
@ 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:1550
@ 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
Compound curve geometry type.
const QgsCurve * curveAt(int i) const SIP_HOLDGIL
Returns the curve at the specified index.
int nCurves() const SIP_HOLDGIL
Returns the number of curves in the geometry.
Curve polygon geometry type.
const QgsCurve * interiorRing(int i) const SIP_HOLDGIL
Retrieves an interior ring from the curve polygon.
const QgsCurve * exteriorRing() const SIP_HOLDGIL
Returns the curve polygon's exterior ring.
int numInteriorRings() const SIP_HOLDGIL
Returns the number of interior rings contained with the curve polygon.
Abstract base class for curved geometry type.
Definition: qgscurve.h:36
virtual int numPoints() const =0
Returns the number of points in the curve.
QgsCurve * segmentize(double tolerance=M_PI_2/90, SegmentationToleranceType toleranceType=MaximumAngle) const override
Returns a geometry without curves.
Definition: qgscurve.cpp:175
virtual bool isClosed() const SIP_HOLDGIL
Returns true if the curve is closed.
Definition: qgscurve.cpp:53
QgsCurve * clone() const override=0
Clones the geometry by performing a deep copy.
virtual QgsPolygon * toPolygon(unsigned int segments=36) const
Returns a segmented polygon.
Definition: qgsellipse.cpp:224
Base class for feedback objects to be used for cancellation of something running in a worker thread.
Definition: qgsfeedback.h:45
bool isCanceled() const SIP_HOLDGIL
Tells whether the operation has been canceled already.
Definition: qgsfeedback.h:54
Geometry collection.
int numGeometries() const SIP_HOLDGIL
Returns the number of geometries within the collection.
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:45
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
Type
The WKB type describes the number of dimensions a geometry has.
Definition: qgswkbtypes.h:70
static bool isCurvedType(Type type) SIP_HOLDGIL
Returns true if the WKB type is a curved type or can contain curved geometries.
Definition: qgswkbtypes.h:911
static Type flatType(Type type) SIP_HOLDGIL
Returns the flat type for a WKB type.
Definition: qgswkbtypes.h:732
double ANALYSIS_EXPORT angle(QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, QgsPoint *p4)
Calculates the angle between two segments (in 2 dimension, z-values are ignored)
Definition: MathUtils.cpp:786
std::unique_ptr< GEOSGeometry, GeosDeleter > unique_ptr
Scoped GEOS pointer.
Definition: qgsgeos.h: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:2260
std::unique_ptr< QgsCurve > lineToCurve(const QgsCurve *curve, double distanceTolerance, double pointSpacingAngleTolerance)
bool matchesOrientation(std::array< Direction, 4 > dirs, std::array< Direction, 4 > oriented)
double normalizedDotProduct(const QgsPoint &a, const QgsPoint &b, const QgsPoint &c)
bool isClockwise(std::array< Direction, 4 > dirs)
Checks whether the 4 directions in dirs make up a clockwise rectangle.
std::unique_ptr< QgsLineString > squareWavesRandomizedAlongLine(const QgsLineString *line, const double minimumWavelength, const double maximumWavelength, const double minimumAmplitude, const double maximumAmplitude, std::uniform_real_distribution<> &uniformDist, std::mt19937 &mt)
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< 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