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 *
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 )
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
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
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  {
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  {
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  {
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
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;
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  {
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;
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  {
1213  }
1214  else
1215  {
1216  std::unique_ptr< QgsPolygon > p( qgsgeometry_cast< QgsPolygon * >( ms->geometryN( i )->segmentize() ) );
1218  }
1219  }
1220  }
1221  else
1222  {
1223  if ( const QgsPolygon *poly = qgsgeometry_cast< const QgsPolygon * >( polygon.constGet() ) )
1224  {
1226  }
1227  else
1228  {
1229  std::unique_ptr< QgsPolygon > p( qgsgeometry_cast< QgsPolygon * >( polygon.constGet()->segmentize() ) );
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  {
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;
1390
1391  while ( i < numEdges - 2 )
1392  {
1393  unsigned int arcEdges = 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  {
1438  }
1439  else
1440  {
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 */
1454  {
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 );
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 );
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  {
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
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
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
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.