QGIS API Documentation  3.6.0-Noosa (5873452)
qgsinternalgeometryengine.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgsinternalgeometryengine.cpp - QgsInternalGeometryEngine
3 
4  ---------------------
5  begin : 13.1.2016
6  copyright : (C) 2016 by Matthias Kuhn
7  email : [email protected]
8  ***************************************************************************
9  * *
10  * This program is free software; you can redistribute it and/or modify *
11  * it under the terms of the GNU General Public License as published by *
12  * the Free Software Foundation; either version 2 of the License, or *
13  * (at your option) any later version. *
14  * *
15  ***************************************************************************/
16 
18 
19 #include "qgslinestring.h"
20 #include "qgsmultipolygon.h"
21 #include "qgspolygon.h"
22 #include "qgsmulticurve.h"
23 #include "qgsgeometry.h"
24 #include "qgsgeometryutils.h"
25 #include "qgslinesegment.h"
26 #include "qgscircle.h"
27 #include "qgslogger.h"
28 #include <QTransform>
29 #include <functional>
30 #include <memory>
31 #include <queue>
32 
34  : mGeometry( geometry.constGet() )
35 {
36 
37 }
38 
39 /***************************************************************************
40  * This class is considered CRITICAL and any change MUST be accompanied with
41  * full unit tests.
42  * See details in QEP #17
43  ****************************************************************************/
44 
46 {
47  QVector<QgsLineString *> linesToProcess;
48 
49  const QgsMultiCurve *multiCurve = qgsgeometry_cast< const QgsMultiCurve * >( mGeometry );
50  if ( multiCurve )
51  {
52  linesToProcess.reserve( multiCurve->partCount() );
53  for ( int i = 0; i < multiCurve->partCount(); ++i )
54  {
55  linesToProcess << static_cast<QgsLineString *>( multiCurve->geometryN( i )->clone() );
56  }
57  }
58 
59  const QgsCurve *curve = qgsgeometry_cast< const QgsCurve * >( mGeometry );
60  if ( curve )
61  {
62  linesToProcess << static_cast<QgsLineString *>( curve->segmentize() );
63  }
64 
65  std::unique_ptr<QgsMultiPolygon> multipolygon( linesToProcess.size() > 1 ? new QgsMultiPolygon() : nullptr );
66  QgsPolygon *polygon = nullptr;
67 
68  if ( !linesToProcess.empty() )
69  {
70  std::unique_ptr< QgsLineString > secondline;
71  for ( QgsLineString *line : qgis::as_const( linesToProcess ) )
72  {
73  QTransform transform = QTransform::fromTranslate( x, y );
74 
75  secondline.reset( line->reversed() );
76  secondline->transform( transform );
77 
78  line->append( secondline.get() );
79  line->addVertex( line->pointN( 0 ) );
80 
81  polygon = new QgsPolygon();
82  polygon->setExteriorRing( line );
83 
84  if ( multipolygon )
85  multipolygon->addGeometry( polygon );
86  }
87 
88  if ( multipolygon )
89  return QgsGeometry( multipolygon.release() );
90  else
91  return QgsGeometry( polygon );
92  }
93 
94  return QgsGeometry();
95 }
96 
97 
98 
99 // polylabel implementation
100 // ported from the original JavaScript implementation developed by Vladimir Agafonkin
101 // originally licensed under the ISC License
102 
104 class Cell
105 {
106  public:
107  Cell( double x, double y, double h, const QgsPolygon *polygon )
108  : x( x )
109  , y( y )
110  , h( h )
111  , d( polygon->pointDistanceToBoundary( x, y ) )
112  , max( d + h * M_SQRT2 )
113  {}
114 
116  double x;
118  double y;
120  double h;
122  double d;
124  double max;
125 };
126 
127 struct GreaterThanByMax
128 {
129  bool operator()( const Cell *lhs, const Cell *rhs )
130  {
131  return rhs->max > lhs->max;
132  }
133 };
134 
135 Cell *getCentroidCell( const QgsPolygon *polygon )
136 {
137  double area = 0;
138  double x = 0;
139  double y = 0;
140 
141  const QgsLineString *exterior = static_cast< const QgsLineString *>( polygon->exteriorRing() );
142  int len = exterior->numPoints() - 1; //assume closed
143  for ( int i = 0, j = len - 1; i < len; j = i++ )
144  {
145  double aX = exterior->xAt( i );
146  double aY = exterior->yAt( i );
147  double bX = exterior->xAt( j );
148  double bY = exterior->yAt( j );
149  double f = aX * bY - bX * aY;
150  x += ( aX + bX ) * f;
151  y += ( aY + bY ) * f;
152  area += f * 3;
153  }
154  if ( area == 0 )
155  return new Cell( exterior->xAt( 0 ), exterior->yAt( 0 ), 0, polygon );
156  else
157  return new Cell( x / area, y / area, 0.0, polygon );
158 }
159 
160 QgsPoint surfacePoleOfInaccessibility( const QgsSurface *surface, double precision, double &distanceFromBoundary )
161 {
162  std::unique_ptr< QgsPolygon > segmentizedPoly;
163  const QgsPolygon *polygon = qgsgeometry_cast< const QgsPolygon * >( surface );
164  if ( !polygon )
165  {
166  segmentizedPoly.reset( static_cast< QgsPolygon *>( surface->segmentize() ) );
167  polygon = segmentizedPoly.get();
168  }
169 
170  // start with the bounding box
171  QgsRectangle bounds = polygon->boundingBox();
172 
173  // initial parameters
174  double cellSize = std::min( bounds.width(), bounds.height() );
175 
176  if ( qgsDoubleNear( cellSize, 0.0 ) )
177  return QgsPoint( bounds.xMinimum(), bounds.yMinimum() );
178 
179  double h = cellSize / 2.0;
180  std::priority_queue< Cell *, std::vector<Cell *>, GreaterThanByMax > cellQueue;
181 
182  // cover polygon with initial cells
183  for ( double x = bounds.xMinimum(); x < bounds.xMaximum(); x += cellSize )
184  {
185  for ( double y = bounds.yMinimum(); y < bounds.yMaximum(); y += cellSize )
186  {
187  cellQueue.push( new Cell( x + h, y + h, h, polygon ) );
188  }
189  }
190 
191  // take centroid as the first best guess
192  std::unique_ptr< Cell > bestCell( getCentroidCell( polygon ) );
193 
194  // special case for rectangular polygons
195  std::unique_ptr< Cell > bboxCell( new Cell( bounds.xMinimum() + bounds.width() / 2.0,
196  bounds.yMinimum() + bounds.height() / 2.0,
197  0, polygon ) );
198  if ( bboxCell->d > bestCell->d )
199  {
200  bestCell = std::move( bboxCell );
201  }
202 
203  while ( !cellQueue.empty() )
204  {
205  // pick the most promising cell from the queue
206  std::unique_ptr< Cell > cell( cellQueue.top() );
207  cellQueue.pop();
208  Cell *currentCell = cell.get();
209 
210  // update the best cell if we found a better one
211  if ( currentCell->d > bestCell->d )
212  {
213  bestCell = std::move( cell );
214  }
215 
216  // do not drill down further if there's no chance of a better solution
217  if ( currentCell->max - bestCell->d <= precision )
218  continue;
219 
220  // split the cell into four cells
221  h = currentCell->h / 2.0;
222  cellQueue.push( new Cell( currentCell->x - h, currentCell->y - h, h, polygon ) );
223  cellQueue.push( new Cell( currentCell->x + h, currentCell->y - h, h, polygon ) );
224  cellQueue.push( new Cell( currentCell->x - h, currentCell->y + h, h, polygon ) );
225  cellQueue.push( new Cell( currentCell->x + h, currentCell->y + h, h, polygon ) );
226  }
227 
228  distanceFromBoundary = bestCell->d;
229 
230  return QgsPoint( bestCell->x, bestCell->y );
231 }
232 
234 
235 QgsGeometry QgsInternalGeometryEngine::poleOfInaccessibility( double precision, double *distanceFromBoundary ) const
236 {
237  if ( distanceFromBoundary )
238  *distanceFromBoundary = std::numeric_limits<double>::max();
239 
240  if ( !mGeometry || mGeometry->isEmpty() )
241  return QgsGeometry();
242 
243  if ( precision <= 0 )
244  return QgsGeometry();
245 
246  if ( const QgsGeometryCollection *gc = qgsgeometry_cast< const QgsGeometryCollection *>( mGeometry ) )
247  {
248  int numGeom = gc->numGeometries();
249  double maxDist = 0;
250  QgsPoint bestPoint;
251  bool found = false;
252  for ( int i = 0; i < numGeom; ++i )
253  {
254  const QgsSurface *surface = qgsgeometry_cast< const QgsSurface * >( gc->geometryN( i ) );
255  if ( !surface )
256  continue;
257 
258  found = true;
259  double dist = std::numeric_limits<double>::max();
260  QgsPoint p = surfacePoleOfInaccessibility( surface, precision, dist );
261  if ( dist > maxDist )
262  {
263  maxDist = dist;
264  bestPoint = p;
265  }
266  }
267 
268  if ( !found )
269  return QgsGeometry();
270 
271  if ( distanceFromBoundary )
272  *distanceFromBoundary = maxDist;
273  return QgsGeometry( new QgsPoint( bestPoint ) );
274  }
275  else
276  {
277  const QgsSurface *surface = qgsgeometry_cast< const QgsSurface * >( mGeometry );
278  if ( !surface )
279  return QgsGeometry();
280 
281  double dist = std::numeric_limits<double>::max();
282  QgsPoint p = surfacePoleOfInaccessibility( surface, precision, dist );
283  if ( distanceFromBoundary )
284  *distanceFromBoundary = dist;
285  return QgsGeometry( new QgsPoint( p ) );
286  }
287 }
288 
289 
290 // helpers for orthogonalize
291 // adapted from original code in potlatch/id osm editor
292 
293 bool dotProductWithinAngleTolerance( double dotProduct, double lowerThreshold, double upperThreshold )
294 {
295  return lowerThreshold > std::fabs( dotProduct ) || std::fabs( dotProduct ) > upperThreshold;
296 }
297 
298 double normalizedDotProduct( const QgsPoint &a, const QgsPoint &b, const QgsPoint &c )
299 {
300  QgsVector p = a - b;
301  QgsVector q = c - b;
302 
303  if ( p.length() > 0 )
304  p = p.normalized();
305  if ( q.length() > 0 )
306  q = q.normalized();
307 
308  return p * q;
309 }
310 
311 double squareness( QgsLineString *ring, double lowerThreshold, double upperThreshold )
312 {
313  double sum = 0.0;
314 
315  bool isClosed = ring->isClosed();
316  int numPoints = ring->numPoints();
317  QgsPoint a;
318  QgsPoint b;
319  QgsPoint c;
320 
321  for ( int i = 0; i < numPoints; ++i )
322  {
323  if ( !isClosed && i == numPoints - 1 )
324  break;
325  else if ( !isClosed && i == 0 )
326  {
327  b = ring->pointN( 0 );
328  c = ring->pointN( 1 );
329  }
330  else
331  {
332  if ( i == 0 )
333  {
334  a = ring->pointN( numPoints - 1 );
335  b = ring->pointN( 0 );
336  }
337  if ( i == numPoints - 1 )
338  c = ring->pointN( 0 );
339  else
340  c = ring->pointN( i + 1 );
341 
342  double dotProduct = normalizedDotProduct( a, b, c );
343  if ( !dotProductWithinAngleTolerance( dotProduct, lowerThreshold, upperThreshold ) )
344  continue;
345 
346  sum += 2.0 * std::min( std::fabs( dotProduct - 1.0 ), std::min( std::fabs( dotProduct ), std::fabs( dotProduct + 1 ) ) );
347  }
348  a = b;
349  b = c;
350  }
351 
352  return sum;
353 }
354 
355 QgsVector calcMotion( const QgsPoint &a, const QgsPoint &b, const QgsPoint &c,
356  double lowerThreshold, double upperThreshold )
357 {
358  QgsVector p = a - b;
359  QgsVector q = c - b;
360 
361  if ( qgsDoubleNear( p.length(), 0.0 ) || qgsDoubleNear( q.length(), 0.0 ) )
362  return QgsVector( 0, 0 );
363 
364  // 2.0 is a magic number from the original JOSM source code
365  double scale = 2.0 * std::min( p.length(), q.length() );
366 
367  p = p.normalized();
368  q = q.normalized();
369  double dotProduct = p * q;
370 
371  if ( !dotProductWithinAngleTolerance( dotProduct, lowerThreshold, upperThreshold ) )
372  {
373  return QgsVector( 0, 0 );
374  }
375 
376  // wonderful nasty hack which has survived through JOSM -> id -> QGIS
377  // to deal with almost-straight segments (angle is closer to 180 than to 90/270).
378  if ( dotProduct < -M_SQRT1_2 )
379  dotProduct += 1.0;
380 
381  QgsVector new_v = p + q;
382  // 0.1 magic number from JOSM implementation - think this is to limit each iterative step
383  return new_v.normalized() * ( 0.1 * dotProduct * scale );
384 }
385 
386 QgsLineString *doOrthogonalize( QgsLineString *ring, int iterations, double tolerance, double lowerThreshold, double upperThreshold )
387 {
388  double minScore = std::numeric_limits<double>::max();
389 
390  bool isClosed = ring->isClosed();
391  int numPoints = ring->numPoints();
392 
393  std::unique_ptr< QgsLineString > best( ring->clone() );
394 
395  QVector< QgsVector > /* yep */ motions;
396  motions.reserve( numPoints );
397 
398  for ( int it = 0; it < iterations; ++it )
399  {
400  motions.resize( 0 ); // avoid re-allocations
401 
402  // first loop through an calculate all motions
403  QgsPoint a;
404  QgsPoint b;
405  QgsPoint c;
406  for ( int i = 0; i < numPoints; ++i )
407  {
408  if ( isClosed && i == numPoints - 1 )
409  motions << motions.at( 0 ); //closed ring, so last point follows first point motion
410  else if ( !isClosed && ( i == 0 || i == numPoints - 1 ) )
411  {
412  b = ring->pointN( 0 );
413  c = ring->pointN( 1 );
414  motions << QgsVector( 0, 0 ); //non closed line, leave first/last vertex untouched
415  }
416  else
417  {
418  if ( i == 0 )
419  {
420  a = ring->pointN( numPoints - 1 );
421  b = ring->pointN( 0 );
422  }
423  if ( i == numPoints - 1 )
424  c = ring->pointN( 0 );
425  else
426  c = ring->pointN( i + 1 );
427 
428  motions << calcMotion( a, b, c, lowerThreshold, upperThreshold );
429  }
430  a = b;
431  b = c;
432  }
433 
434  // then apply them
435  for ( int i = 0; i < numPoints; ++i )
436  {
437  ring->setXAt( i, ring->xAt( i ) + motions.at( i ).x() );
438  ring->setYAt( i, ring->yAt( i ) + motions.at( i ).y() );
439  }
440 
441  double newScore = squareness( ring, lowerThreshold, upperThreshold );
442  if ( newScore < minScore )
443  {
444  best.reset( ring->clone() );
445  minScore = newScore;
446  }
447 
448  if ( minScore < tolerance )
449  break;
450  }
451 
452  delete ring;
453 
454  return best.release();
455 }
456 
457 
458 QgsAbstractGeometry *orthogonalizeGeom( const QgsAbstractGeometry *geom, int maxIterations, double tolerance, double lowerThreshold, double upperThreshold )
459 {
460  std::unique_ptr< QgsAbstractGeometry > segmentizedCopy;
461  if ( QgsWkbTypes::isCurvedType( geom->wkbType() ) )
462  {
463  segmentizedCopy.reset( geom->segmentize() );
464  geom = segmentizedCopy.get();
465  }
466 
468  {
469  return doOrthogonalize( static_cast< QgsLineString * >( geom->clone() ),
470  maxIterations, tolerance, lowerThreshold, upperThreshold );
471  }
472  else
473  {
474  // polygon
475  const QgsPolygon *polygon = static_cast< const QgsPolygon * >( geom );
476  QgsPolygon *result = new QgsPolygon();
477 
478  result->setExteriorRing( doOrthogonalize( static_cast< QgsLineString * >( polygon->exteriorRing()->clone() ),
479  maxIterations, tolerance, lowerThreshold, upperThreshold ) );
480  for ( int i = 0; i < polygon->numInteriorRings(); ++i )
481  {
482  result->addInteriorRing( doOrthogonalize( static_cast< QgsLineString * >( polygon->interiorRing( i )->clone() ),
483  maxIterations, tolerance, lowerThreshold, upperThreshold ) );
484  }
485 
486  return result;
487  }
488 }
489 
490 QgsGeometry QgsInternalGeometryEngine::orthogonalize( double tolerance, int maxIterations, double angleThreshold ) const
491 {
492  if ( !mGeometry || ( QgsWkbTypes::geometryType( mGeometry->wkbType() ) != QgsWkbTypes::LineGeometry
493  && QgsWkbTypes::geometryType( mGeometry->wkbType() ) != QgsWkbTypes::PolygonGeometry ) )
494  {
495  return QgsGeometry();
496  }
497 
498  double lowerThreshold = std::cos( ( 90 - angleThreshold ) * M_PI / 180.00 );
499  double upperThreshold = std::cos( angleThreshold * M_PI / 180.0 );
500 
501  if ( const QgsGeometryCollection *gc = qgsgeometry_cast< const QgsGeometryCollection *>( mGeometry ) )
502  {
503  int numGeom = gc->numGeometries();
504  QVector< QgsAbstractGeometry * > geometryList;
505  geometryList.reserve( numGeom );
506  for ( int i = 0; i < numGeom; ++i )
507  {
508  geometryList << orthogonalizeGeom( gc->geometryN( i ), maxIterations, tolerance, lowerThreshold, upperThreshold );
509  }
510 
511  QgsGeometry first = QgsGeometry( geometryList.takeAt( 0 ) );
512  for ( QgsAbstractGeometry *g : qgis::as_const( geometryList ) )
513  {
514  first.addPart( g );
515  }
516  return first;
517  }
518  else
519  {
520  return QgsGeometry( orthogonalizeGeom( mGeometry, maxIterations, tolerance, lowerThreshold, upperThreshold ) );
521  }
522 }
523 
524 // if extraNodesPerSegment < 0, then use distance based mode
525 QgsLineString *doDensify( const QgsLineString *ring, int extraNodesPerSegment = -1, double distance = 1 )
526 {
527  QVector< double > outX;
528  QVector< double > outY;
529  QVector< double > outZ;
530  QVector< double > outM;
531  double multiplier = 1.0 / double( extraNodesPerSegment + 1 );
532 
533  int nPoints = ring->numPoints();
534  outX.reserve( ( extraNodesPerSegment + 1 ) * nPoints );
535  outY.reserve( ( extraNodesPerSegment + 1 ) * nPoints );
536  bool withZ = ring->is3D();
537  if ( withZ )
538  outZ.reserve( ( extraNodesPerSegment + 1 ) * nPoints );
539  bool withM = ring->isMeasure();
540  if ( withM )
541  outM.reserve( ( extraNodesPerSegment + 1 ) * nPoints );
542  double x1 = 0;
543  double x2 = 0;
544  double y1 = 0;
545  double y2 = 0;
546  double z1 = 0;
547  double z2 = 0;
548  double m1 = 0;
549  double m2 = 0;
550  double xOut = 0;
551  double yOut = 0;
552  double zOut = 0;
553  double mOut = 0;
554  int extraNodesThisSegment = extraNodesPerSegment;
555  for ( int i = 0; i < nPoints - 1; ++i )
556  {
557  x1 = ring->xAt( i );
558  x2 = ring->xAt( i + 1 );
559  y1 = ring->yAt( i );
560  y2 = ring->yAt( i + 1 );
561  if ( withZ )
562  {
563  z1 = ring->zAt( i );
564  z2 = ring->zAt( i + 1 );
565  }
566  if ( withM )
567  {
568  m1 = ring->mAt( i );
569  m2 = ring->mAt( i + 1 );
570  }
571 
572  outX << x1;
573  outY << y1;
574  if ( withZ )
575  outZ << z1;
576  if ( withM )
577  outM << m1;
578 
579  if ( extraNodesPerSegment < 0 )
580  {
581  // distance mode
582  extraNodesThisSegment = std::floor( std::sqrt( ( x2 - x1 ) * ( x2 - x1 ) + ( y2 - y1 ) * ( y2 - y1 ) ) / distance );
583  if ( extraNodesThisSegment >= 1 )
584  multiplier = 1.0 / ( extraNodesThisSegment + 1 );
585  }
586 
587  for ( int j = 0; j < extraNodesThisSegment; ++j )
588  {
589  double delta = multiplier * ( j + 1 );
590  xOut = x1 + delta * ( x2 - x1 );
591  yOut = y1 + delta * ( y2 - y1 );
592  if ( withZ )
593  zOut = z1 + delta * ( z2 - z1 );
594  if ( withM )
595  mOut = m1 + delta * ( m2 - m1 );
596 
597  outX << xOut;
598  outY << yOut;
599  if ( withZ )
600  outZ << zOut;
601  if ( withM )
602  outM << mOut;
603  }
604  }
605  outX << ring->xAt( nPoints - 1 );
606  outY << ring->yAt( nPoints - 1 );
607  if ( withZ )
608  outZ << ring->zAt( nPoints - 1 );
609  if ( withM )
610  outM << ring->mAt( nPoints - 1 );
611 
612  QgsLineString *result = new QgsLineString( outX, outY, outZ, outM );
613  return result;
614 }
615 
616 QgsAbstractGeometry *densifyGeometry( const QgsAbstractGeometry *geom, int extraNodesPerSegment = 1, double distance = 1 )
617 {
618  std::unique_ptr< QgsAbstractGeometry > segmentizedCopy;
619  if ( QgsWkbTypes::isCurvedType( geom->wkbType() ) )
620  {
621  segmentizedCopy.reset( geom->segmentize() );
622  geom = segmentizedCopy.get();
623  }
624 
626  {
627  return doDensify( static_cast< const QgsLineString * >( geom ), extraNodesPerSegment, distance );
628  }
629  else
630  {
631  // polygon
632  const QgsPolygon *polygon = static_cast< const QgsPolygon * >( geom );
633  QgsPolygon *result = new QgsPolygon();
634 
635  result->setExteriorRing( doDensify( static_cast< const QgsLineString * >( polygon->exteriorRing() ),
636  extraNodesPerSegment, distance ) );
637  for ( int i = 0; i < polygon->numInteriorRings(); ++i )
638  {
639  result->addInteriorRing( doDensify( static_cast< const QgsLineString * >( polygon->interiorRing( i ) ),
640  extraNodesPerSegment, distance ) );
641  }
642 
643  return result;
644  }
645 }
646 
648 {
649  if ( !mGeometry )
650  {
651  return QgsGeometry();
652  }
653 
654  if ( QgsWkbTypes::geometryType( mGeometry->wkbType() ) == QgsWkbTypes::PointGeometry )
655  {
656  return QgsGeometry( mGeometry->clone() ); // point geometry, nothing to do
657  }
658 
659  if ( const QgsGeometryCollection *gc = qgsgeometry_cast< const QgsGeometryCollection *>( mGeometry ) )
660  {
661  int numGeom = gc->numGeometries();
662  QVector< QgsAbstractGeometry * > geometryList;
663  geometryList.reserve( numGeom );
664  for ( int i = 0; i < numGeom; ++i )
665  {
666  geometryList << densifyGeometry( gc->geometryN( i ), extraNodesPerSegment );
667  }
668 
669  QgsGeometry first = QgsGeometry( geometryList.takeAt( 0 ) );
670  for ( QgsAbstractGeometry *g : qgis::as_const( geometryList ) )
671  {
672  first.addPart( g );
673  }
674  return first;
675  }
676  else
677  {
678  return QgsGeometry( densifyGeometry( mGeometry, extraNodesPerSegment ) );
679  }
680 }
681 
683 {
684  if ( !mGeometry )
685  {
686  return QgsGeometry();
687  }
688 
689  if ( QgsWkbTypes::geometryType( mGeometry->wkbType() ) == QgsWkbTypes::PointGeometry )
690  {
691  return QgsGeometry( mGeometry->clone() ); // point geometry, nothing to do
692  }
693 
694  if ( const QgsGeometryCollection *gc = qgsgeometry_cast< const QgsGeometryCollection *>( mGeometry ) )
695  {
696  int numGeom = gc->numGeometries();
697  QVector< QgsAbstractGeometry * > geometryList;
698  geometryList.reserve( numGeom );
699  for ( int i = 0; i < numGeom; ++i )
700  {
701  geometryList << densifyGeometry( gc->geometryN( i ), -1, distance );
702  }
703 
704  QgsGeometry first = QgsGeometry( geometryList.takeAt( 0 ) );
705  for ( QgsAbstractGeometry *g : qgis::as_const( geometryList ) )
706  {
707  first.addPart( g );
708  }
709  return first;
710  }
711  else
712  {
713  return QgsGeometry( densifyGeometry( mGeometry, -1, distance ) );
714  }
715 }
716 
718 //
719 // QgsLineSegmentDistanceComparer
720 //
721 
722 // adapted for QGIS geometry classes from original work at https://github.com/trylock/visibility by trylock
723 bool QgsLineSegmentDistanceComparer::operator()( QgsLineSegment2D ab, QgsLineSegment2D cd ) const
724 {
725  Q_ASSERT_X( ab.pointLeftOfLine( mOrigin ) != 0,
726  "line_segment_dist_comparer",
727  "AB must not be collinear with the origin." );
728  Q_ASSERT_X( cd.pointLeftOfLine( mOrigin ) != 0,
729  "line_segment_dist_comparer",
730  "CD must not be collinear with the origin." );
731 
732  // flip the segments so that if there are common endpoints,
733  // they will be the segment's start points
734  if ( ab.end() == cd.start() || ab.end() == cd.end() )
735  ab.reverse();
736  if ( ab.start() == cd.end() )
737  cd.reverse();
738 
739  // cases with common endpoints
740  if ( ab.start() == cd.start() )
741  {
742  const int oad = QgsGeometryUtils::leftOfLine( cd.endX(), cd.endY(), mOrigin.x(), mOrigin.y(), ab.startX(), ab.startY() );
743  const int oab = ab.pointLeftOfLine( mOrigin );
744  if ( ab.end() == cd.end() || oad != oab )
745  return false;
746  else
747  return ab.pointLeftOfLine( cd.end() ) != oab;
748  }
749  else
750  {
751  // cases without common endpoints
752  const int cda = cd.pointLeftOfLine( ab.start() );
753  const int cdb = cd.pointLeftOfLine( ab.end() );
754  if ( cdb == 0 && cda == 0 )
755  {
756  return mOrigin.sqrDist( ab.start() ) < mOrigin.sqrDist( cd.start() );
757  }
758  else if ( cda == cdb || cda == 0 || cdb == 0 )
759  {
760  const int cdo = cd.pointLeftOfLine( mOrigin );
761  return cdo == cda || cdo == cdb;
762  }
763  else
764  {
765  const int abo = ab.pointLeftOfLine( mOrigin );
766  return abo != ab.pointLeftOfLine( cd.start() );
767  }
768  }
769 }
770 
771 //
772 // QgsClockwiseAngleComparer
773 //
774 
775 bool QgsClockwiseAngleComparer::operator()( const QgsPointXY &a, const QgsPointXY &b ) const
776 {
777  const bool aIsLeft = a.x() < mVertex.x();
778  const bool bIsLeft = b.x() < mVertex.x();
779  if ( aIsLeft != bIsLeft )
780  return bIsLeft;
781 
782  if ( qgsDoubleNear( a.x(), mVertex.x() ) && qgsDoubleNear( b.x(), mVertex.x() ) )
783  {
784  if ( a.y() >= mVertex.y() || b.y() >= mVertex.y() )
785  {
786  return b.y() < a.y();
787  }
788  else
789  {
790  return a.y() < b.y();
791  }
792  }
793  else
794  {
795  const QgsVector oa = a - mVertex;
796  const QgsVector ob = b - mVertex;
797  const double det = oa.crossProduct( ob );
798  if ( qgsDoubleNear( det, 0.0 ) )
799  {
800  return oa.lengthSquared() < ob.lengthSquared();
801  }
802  else
803  {
804  return det < 0;
805  }
806  }
807 }
808 
810 
811 //
812 // QgsRay2D
813 //
814 
815 bool QgsRay2D::intersects( const QgsLineSegment2D &segment, QgsPointXY &intersectPoint ) const
816 {
817  const QgsVector ao = origin - segment.start();
818  const QgsVector ab = segment.end() - segment.start();
819  const double det = ab.crossProduct( direction );
820  if ( qgsDoubleNear( det, 0.0 ) )
821  {
822  const int abo = segment.pointLeftOfLine( origin );
823  if ( abo != 0 )
824  {
825  return false;
826  }
827  else
828  {
829  const double distA = ao * direction;
830  const double distB = ( origin - segment.end() ) * direction;
831 
832  if ( distA > 0 && distB > 0 )
833  {
834  return false;
835  }
836  else
837  {
838  if ( ( distA > 0 ) != ( distB > 0 ) )
839  intersectPoint = origin;
840  else if ( distA > distB ) // at this point, both distances are negative
841  intersectPoint = segment.start(); // hence the nearest point is A
842  else
843  intersectPoint = segment.end();
844  return true;
845  }
846  }
847  }
848  else
849  {
850  const double u = ao.crossProduct( direction ) / det;
851  if ( u < 0.0 || 1.0 < u )
852  {
853  return false;
854  }
855  else
856  {
857  const double t = -ab.crossProduct( ao ) / det;
858  intersectPoint = origin + direction * t;
859  return qgsDoubleNear( t, 0.0 ) || t > 0;
860  }
861  }
862 }
863 
864 QVector<QgsPointXY> generateSegmentCurve( const QgsPoint &center1, const double radius1, const QgsPoint &center2, const double radius2 )
865 {
866  // ensure that first circle is smaller than second
867  if ( radius1 > radius2 )
868  return generateSegmentCurve( center2, radius2, center1, radius1 );
869 
870  QgsPointXY t1;
871  QgsPointXY t2;
872  QgsPointXY t3;
873  QgsPointXY t4;
874  QVector<QgsPointXY> points;
875  if ( QgsGeometryUtils::circleCircleOuterTangents( center1, radius1, center2, radius2, t1, t2, t3, t4 ) )
876  {
877  points << t1
878  << t2
879  << t4
880  << t3;
881  }
882  return points;
883 }
884 
885 // partially ported from JTS VariableWidthBuffer,
886 // https://github.com/topobyte/jts/blob/master/jts-lab/src/main/java/com/vividsolutions/jts/operation/buffer/VariableWidthBuffer.java
887 
888 QgsGeometry QgsInternalGeometryEngine::variableWidthBuffer( int segments, const std::function< std::unique_ptr< double[] >( const QgsLineString *line ) > &widthFunction ) const
889 {
890  if ( !mGeometry )
891  {
892  return QgsGeometry();
893  }
894 
895  std::vector< std::unique_ptr<QgsLineString > > linesToProcess;
896 
897  const QgsMultiCurve *multiCurve = qgsgeometry_cast< const QgsMultiCurve * >( mGeometry );
898  if ( multiCurve )
899  {
900  for ( int i = 0; i < multiCurve->partCount(); ++i )
901  {
902  if ( static_cast< const QgsCurve * >( multiCurve->geometryN( i ) )->nCoordinates() == 0 )
903  continue; // skip 0 length lines
904 
905  linesToProcess.emplace_back( static_cast<QgsLineString *>( multiCurve->geometryN( i )->clone() ) );
906  }
907  }
908 
909  const QgsCurve *curve = qgsgeometry_cast< const QgsCurve * >( mGeometry );
910  if ( curve )
911  {
912  if ( curve->nCoordinates() > 0 )
913  linesToProcess.emplace_back( static_cast<QgsLineString *>( curve->segmentize() ) );
914  }
915 
916  if ( linesToProcess.empty() )
917  {
918  QgsGeometry g;
919  g.mLastError = QStringLiteral( "Input geometry was not a curve type geometry" );
920  return g;
921  }
922 
923  QVector<QgsGeometry> bufferedLines;
924  bufferedLines.reserve( linesToProcess.size() );
925 
926  for ( std::unique_ptr< QgsLineString > &line : linesToProcess )
927  {
928  QVector<QgsGeometry> parts;
929  QgsPoint prevPoint;
930  double prevRadius = 0;
931  QgsGeometry prevCircle;
932 
933  std::unique_ptr< double[] > widths = widthFunction( line.get() ) ;
934  for ( int i = 0; i < line->nCoordinates(); ++i )
935  {
936  QgsPoint thisPoint = line->pointN( i );
937  QgsGeometry thisCircle;
938  double thisRadius = widths[ i ] / 2.0;
939  if ( thisRadius > 0 )
940  {
941  QgsGeometry p = QgsGeometry( thisPoint.clone() );
942 
943  QgsCircle circ( thisPoint, thisRadius );
944  thisCircle = QgsGeometry( circ.toPolygon( segments * 4 ) );
945  parts << thisCircle;
946  }
947  else
948  {
949  thisCircle = QgsGeometry( thisPoint.clone() );
950  }
951 
952  if ( i > 0 )
953  {
954  if ( prevRadius > 0 || thisRadius > 0 )
955  {
956  QVector< QgsPointXY > points = generateSegmentCurve( prevPoint, prevRadius, thisPoint, thisRadius );
957  if ( !points.empty() )
958  {
959  // snap points to circle vertices
960 
961  int atVertex = 0;
962  int beforeVertex = 0;
963  int afterVertex = 0;
964  double sqrDist = 0;
965  double sqrDistPrev = 0;
966  for ( int j = 0; j < points.count(); ++j )
967  {
968  QgsPointXY pA = prevCircle.closestVertex( points.at( j ), atVertex, beforeVertex, afterVertex, sqrDistPrev );
969  QgsPointXY pB = thisCircle.closestVertex( points.at( j ), atVertex, beforeVertex, afterVertex, sqrDist );
970  points[j] = sqrDistPrev < sqrDist ? pA : pB;
971  }
972  // close ring
973  points.append( points.at( 0 ) );
974 
975  std::unique_ptr< QgsPolygon > poly = qgis::make_unique< QgsPolygon >();
976  poly->setExteriorRing( new QgsLineString( points ) );
977  if ( poly->area() > 0 )
978  parts << QgsGeometry( std::move( poly ) );
979  }
980  }
981  }
982  prevPoint = thisPoint;
983  prevRadius = thisRadius;
984  prevCircle = thisCircle;
985  }
986 
987  bufferedLines << QgsGeometry::unaryUnion( parts );
988  }
989 
990  return QgsGeometry::collectGeometry( bufferedLines );
991 }
992 
993 QgsGeometry QgsInternalGeometryEngine::taperedBuffer( double start, double end, int segments ) const
994 {
995  start = std::fabs( start );
996  end = std::fabs( end );
997 
998  auto interpolateWidths = [ start, end ]( const QgsLineString * line )->std::unique_ptr< double [] >
999  {
1000  // ported from JTS VariableWidthBuffer,
1001  // https://github.com/topobyte/jts/blob/master/jts-lab/src/main/java/com/vividsolutions/jts/operation/buffer/VariableWidthBuffer.java
1002  std::unique_ptr< double [] > widths( new double[ line->nCoordinates() ] );
1003  widths[0] = start;
1004  widths[line->nCoordinates() - 1] = end;
1005 
1006  double lineLength = line->length();
1007  double currentLength = 0;
1008  QgsPoint prevPoint = line->pointN( 0 );
1009  for ( int i = 1; i < line->nCoordinates() - 1; ++i )
1010  {
1011  QgsPoint point = line->pointN( i );
1012  double segmentLength = point.distance( prevPoint );
1013  currentLength += segmentLength;
1014  double lengthFraction = lineLength > 0 ? currentLength / lineLength : 1;
1015  double delta = lengthFraction * ( end - start );
1016  widths[i] = start + delta;
1017  prevPoint = point;
1018  }
1019  return widths;
1020  };
1021 
1022  return variableWidthBuffer( segments, interpolateWidths );
1023 }
1024 
1026 {
1027  auto widthByM = []( const QgsLineString * line )->std::unique_ptr< double [] >
1028  {
1029  std::unique_ptr< double [] > widths( new double[ line->nCoordinates() ] );
1030  for ( int i = 0; i < line->nCoordinates(); ++i )
1031  {
1032  widths[ i ] = line->mAt( i );
1033  }
1034  return widths;
1035  };
1036 
1037  return variableWidthBuffer( segments, widthByM );
1038 }
bool isMeasure() const
Returns true if the geometry contains m values.
QgsGeometry taperedBuffer(double startWidth, double endWidth, int segments) const
Calculates a tapered width buffer for a (multi)curve geometry.
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
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)
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 ...
double distance(double x, double y) const
Returns the distance between this point and a specified x, y coordinate.
Definition: qgspoint.h:276
double y
Definition: qgspointxy.h:48
A class to represent a 2D point.
Definition: qgspointxy.h:43
double endY() const
Returns the segment&#39;s end y-coordinate.
bool qgsDoubleNear(double a, double b, double epsilon=4 *std::numeric_limits< double >::epsilon())
Compare two doubles (but allow some difference)
Definition: qgis.h:265
const QgsCurve * interiorRing(int i) const
Retrieves an interior ring from the curve polygon.
bool intersects(const QgsLineSegment2D &segment, QgsPointXY &intersectPoint) const
Finds the closest intersection point of the ray and a line segment.
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:106
double crossProduct(QgsVector v) const
Returns the 2D cross product of this vector and another vector v.
Definition: qgsvector.cpp:102
void setYAt(int index, double y)
Sets the y-coordinate of the specified node in the line string.
double squareness(QgsLineString *ring, double lowerThreshold, double upperThreshold)
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:236
virtual QgsAbstractGeometry * clone() const =0
Clones the geometry by performing a deep copy.
static int circleCircleOuterTangents(const QgsPointXY &center1, double radius1, const QgsPointXY &center2, double radius2, QgsPointXY &line1P1, QgsPointXY &line1P2, QgsPointXY &line2P1, QgsPointXY &line2P2)
Calculates the outer tangent points for two circles, centered at center1 and center2 and with radii o...
void setXAt(int index, double x)
Sets the x-coordinate of the specified node in the line string.
void reverse()
Reverses the line segment, so that the start and end points are flipped.
int numPoints() const override
Returns the number of points in the curve.
QgsPointXY closestVertex(const QgsPointXY &point, int &atVertex, int &beforeVertex, int &afterVertex, double &sqrDist) const
Returns the vertex closest to the given point, the corresponding vertex index, squared distance snap ...
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
QgsPointXY end() const
Returns the segment&#39;s end point.
void addInteriorRing(QgsCurve *ring) override
Adds an interior ring to the geometry (takes ownership)
Definition: qgspolygon.cpp:148
int numInteriorRings() const
Returns the number of interior rings contained with the curve polygon.
double startX() const
Returns the segment&#39;s start x-coordinate.
static int leftOfLine(const double x, const double y, const double x1, const double y1, const double x2, const double y2)
Returns a value < 0 if the point (x, y) is left of the line from (x1, y1) -> ( x2, y2).
Geometry collection.
double width() const
Returns the width of the rectangle.
Definition: qgsrectangle.h:202
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.
static GeometryType geometryType(Type type)
Returns the geometry type for a WKB type, e.g., both MultiPolygon and CurvePolygon would have a Polyg...
Definition: qgswkbtypes.h:665
T qgsgeometry_cast(const QgsAbstractGeometry *geom)
Multi curve geometry collection.
Definition: qgsmulticurve.h:29
QgsRectangle boundingBox() const override
Returns the minimal bounding box for the geometry.
Definition: qgssurface.h:45
int pointLeftOfLine(const QgsPointXY &point) const
Tests if a point is to the left of the line segment.
QgsPoint * clone() const override
Clones the geometry by performing a deep copy.
Definition: qgspoint.cpp:94
QgsGeometry densifyByDistance(double distance) const
Densifies the geometry by adding regularly placed extra nodes inside each segment so that the maximum...
Abstract base class for curved geometry type.
Definition: qgscurve.h:35
Abstract base class for all geometries.
QgsWkbTypes::Type wkbType() const
Returns the WKB type of the geometry.
Point geometry type, with support for z-dimension and m-values.
Definition: qgspoint.h:37
double length() const
Returns the length of the vector.
Definition: qgsvector.cpp:71
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
virtual bool isClosed() const
Returns true if the curve is closed.
Definition: qgscurve.cpp:39
double yMinimum() const
Returns the y minimum value (bottom side of rectangle).
Definition: qgsrectangle.h:177
double xMaximum() const
Returns the x maximum value (right side of rectangle).
Definition: qgsrectangle.h:162
QgsCurve * segmentize(double tolerance=M_PI_2/90, SegmentationToleranceType toleranceType=MaximumAngle) const override
Returns a geometry without curves.
Definition: qgscurve.cpp:146
bool dotProductWithinAngleTolerance(double dotProduct, double lowerThreshold, double upperThreshold)
int partCount() const override
Returns count of parts contained in the geometry.
double endX() const
Returns the segment&#39;s end x-coordinate.
void setExteriorRing(QgsCurve *ring) override
Sets the exterior ring of the polygon.
Definition: qgspolygon.cpp:179
Multi polygon geometry collection.
QgsCurve * clone() const override=0
Clones the geometry by performing a deep copy.
static bool isCurvedType(Type type)
Returns true if the WKB type is a curved type or can contain curved geometries.
Definition: qgswkbtypes.h:608
Line string geometry type, with support for z-dimension and m-values.
Definition: qgslinestring.h:43
const QgsAbstractGeometry * geometryN(int n) const
Returns a const reference to a geometry from within the collection.
QgsPoint pointN(int i) const
Returns the specified point from inside the line string.
double xMinimum() const
Returns the x minimum value (left side of rectangle).
Definition: qgsrectangle.h:167
double xAt(int index) const override
Returns the x-coordinate of the specified node in the line string.
double yMaximum() const
Returns the y maximum value (top side of rectangle).
Definition: qgsrectangle.h:172
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:31
const QgsCurve * exteriorRing() const
Returns the curve polygon&#39;s exterior ring.
QgsGeometry densifyByCount(int extraNodesPerSegment) const
Densifies the geometry by adding the specified number of extra nodes within each segment of the geome...
double startY() const
Returns the segment&#39;s start y-coordinate.
int nCoordinates() const override
Returns the number of nodes contained in the geometry.
double lengthSquared() const
Returns the length of the vector.
Definition: qgsvector.h:103
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.
bool is3D() const
Returns true if the geometry is 3D and contains a z-value.
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...
QgsPointXY start() const
Returns the segment&#39;s start point.
QgsVector normalized() const
Returns the vector&#39;s normalized (or "unit") vector (ie same angle but length of 1.0).
Definition: qgsvector.cpp:114
double height() const
Returns the height of the rectangle.
Definition: qgsrectangle.h:209
QgsLineString * doOrthogonalize(QgsLineString *ring, int iterations, double tolerance, double lowerThreshold, double upperThreshold)