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