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