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