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