QGIS API Documentation 3.29.0-Master (da8bb1db43)
pointset.cpp
Go to the documentation of this file.
1/*
2 * libpal - Automated Placement of Labels Library
3 *
4 * Copyright (C) 2008 Maxence Laurent, MIS-TIC, HEIG-VD
5 * University of Applied Sciences, Western Switzerland
6 * http://www.hes-so.ch
7 *
8 * Contact:
9 * maxence.laurent <at> heig-vd <dot> ch
10 * or
11 * eric.taillard <at> heig-vd <dot> ch
12 *
13 * This file is part of libpal.
14 *
15 * libpal is free software: you can redistribute it and/or modify
16 * it under the terms of the GNU General Public License as published by
17 * the Free Software Foundation, either version 3 of the License, or
18 * (at your option) any later version.
19 *
20 * libpal is distributed in the hope that it will be useful,
21 * but WITHOUT ANY WARRANTY; without even the implied warranty of
22 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23 * GNU General Public License for more details.
24 *
25 * You should have received a copy of the GNU General Public License
26 * along with libpal. If not, see <http://www.gnu.org/licenses/>.
27 *
28 */
29
30#include "pointset.h"
31#include "util.h"
32#include "geomfunction.h"
33#include "qgsgeos.h"
34#include "qgsmessagelog.h"
35#include "qgsgeometryutils.h"
36#include <qglobal.h>
37
38using namespace pal;
39
40PointSet::PointSet()
41{
42 nbPoints = 0;
43 type = -1;
44}
45
46PointSet::PointSet( int nbPoints, double *x, double *y )
47 : nbPoints( nbPoints )
48 , type( GEOS_POLYGON )
49{
50 this->x.resize( nbPoints );
51 this->y.resize( nbPoints );
52 int i;
53
54 for ( i = 0; i < nbPoints; i++ )
55 {
56 this->x[i] = x[i];
57 this->y[i] = y[i];
58 }
59
60}
61
62PointSet::PointSet( double aX, double aY )
63 : type( GEOS_POINT )
64 , xmin( aX )
65 , xmax( aY )
66 , ymin( aX )
67 , ymax( aY )
68{
69 nbPoints = 1;
70 x.resize( 1 );
71 y.resize( 1 );
72 x[0] = aX;
73 y[0] = aY;
74}
75
77 : xmin( ps.xmin )
78 , xmax( ps.xmax )
79 , ymin( ps.ymin )
80 , ymax( ps.ymax )
81{
82 nbPoints = ps.nbPoints;
83 x = ps.x;
84 y = ps.y;
85
87
88 type = ps.type;
89
90 holeOf = ps.holeOf;
91
92 if ( ps.mGeos )
93 {
94 mGeos = GEOSGeom_clone_r( QgsGeos::getGEOSHandler(), ps.mGeos );
95 mOwnsGeom = true;
96 }
97}
98
100{
101 GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
102
103 bool needClose = false;
104 if ( type == GEOS_POLYGON && ( !qgsDoubleNear( x[0], x[ nbPoints - 1] ) || !qgsDoubleNear( y[0], y[ nbPoints - 1 ] ) ) )
105 {
106 needClose = true;
107 }
108
109 GEOSCoordSequence *coord = nullptr;
110#if GEOS_VERSION_MAJOR>3 || ( GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR>=10 )
111 if ( !needClose )
112 {
113 // use optimised method if we don't have to force close an open ring
114 coord = GEOSCoordSeq_copyFromArrays_r( geosctxt, x.data(), y.data(), nullptr, nullptr, nbPoints );
115 }
116#else
117 if ( false )
118 {
119 }
120#endif
121 else
122 {
123 coord = GEOSCoordSeq_create_r( geosctxt, nbPoints + ( needClose ? 1 : 0 ), 2 );
124 for ( int i = 0; i < nbPoints; ++i )
125 {
126 GEOSCoordSeq_setXY_r( geosctxt, coord, i, x[i], y[i] );
127 }
128
129 //close ring if needed
130 if ( needClose )
131 {
132 GEOSCoordSeq_setXY_r( geosctxt, coord, nbPoints, x[0], y[0] );
133 }
134 }
135
136 switch ( type )
137 {
138 case GEOS_POLYGON:
139 mGeos = GEOSGeom_createPolygon_r( geosctxt, GEOSGeom_createLinearRing_r( geosctxt, coord ), nullptr, 0 );
140 break;
141
142 case GEOS_LINESTRING:
143 mGeos = GEOSGeom_createLineString_r( geosctxt, coord );
144 break;
145
146 case GEOS_POINT:
147 mGeos = GEOSGeom_createPoint_r( geosctxt, coord );
148 break;
149 }
150
151 mOwnsGeom = true;
152}
153
154const GEOSPreparedGeometry *PointSet::preparedGeom() const
155{
156 if ( !mGeos )
158
159 if ( !mPreparedGeom )
160 {
161 mPreparedGeom = GEOSPrepare_r( QgsGeos::getGEOSHandler(), mGeos );
162 }
163 return mPreparedGeom;
164}
165
167{
168 GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
169 if ( mOwnsGeom ) // delete old geometry if we own it
170 GEOSGeom_destroy_r( geosctxt, mGeos );
171 mOwnsGeom = false;
172 mGeos = nullptr;
173
174 if ( mPreparedGeom )
175 {
176 GEOSPreparedGeom_destroy_r( geosctxt, mPreparedGeom );
177 mPreparedGeom = nullptr;
178 }
179
180 if ( mGeosPreparedBoundary )
181 {
182 GEOSPreparedGeom_destroy_r( geosctxt, mGeosPreparedBoundary );
183 mGeosPreparedBoundary = nullptr;
184 }
185
186 if ( mMultipartPreparedGeos )
187 {
188 GEOSPreparedGeom_destroy_r( geosctxt, mMultipartPreparedGeos );
189 mMultipartPreparedGeos = nullptr;
190 }
191 if ( mMultipartGeos )
192 {
193 GEOSGeom_destroy_r( geosctxt, mMultipartGeos );
194 mMultipartGeos = nullptr;
195 }
196
197 mLength = -1;
198 mArea = -1;
199}
200
202{
203 GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
204
205 if ( mGeos && mOwnsGeom )
206 {
207 GEOSGeom_destroy_r( geosctxt, mGeos );
208 mGeos = nullptr;
209 }
210 GEOSPreparedGeom_destroy_r( geosctxt, mPreparedGeom );
211
212 if ( mGeosPreparedBoundary )
213 {
214 GEOSPreparedGeom_destroy_r( geosctxt, mGeosPreparedBoundary );
215 mGeosPreparedBoundary = nullptr;
216 }
217
218 if ( mMultipartPreparedGeos )
219 {
220 GEOSPreparedGeom_destroy_r( geosctxt, mMultipartPreparedGeos );
221 mMultipartPreparedGeos = nullptr;
222 }
223 if ( mMultipartGeos )
224 {
225 GEOSGeom_destroy_r( geosctxt, mMultipartGeos );
226 mMultipartGeos = nullptr;
227 }
228
229 deleteCoords();
230}
231
233{
234 x.clear();
235 y.clear();
236}
237
238std::unique_ptr<PointSet> PointSet::extractShape( int nbPtSh, int imin, int imax, int fps, int fpe, double fptx, double fpty )
239{
240 int i, j;
241
242 std::unique_ptr<PointSet> newShape = std::make_unique< PointSet >();
243 newShape->type = GEOS_POLYGON;
244 newShape->nbPoints = nbPtSh;
245 newShape->x.resize( newShape->nbPoints );
246 newShape->y.resize( newShape->nbPoints );
247
248 // new shape # 1 from imin to imax
249 for ( j = 0, i = imin; i != ( imax + 1 ) % nbPoints; i = ( i + 1 ) % nbPoints, j++ )
250 {
251 newShape->x[j] = x[i];
252 newShape->y[j] = y[i];
253 }
254 // is the cutting point a new one ?
255 if ( fps != fpe )
256 {
257 // yes => so add it
258 newShape->x[j] = fptx;
259 newShape->y[j] = fpty;
260 }
261
262 return newShape;
263}
264
265std::unique_ptr<PointSet> PointSet::clone() const
266{
267 return std::unique_ptr< PointSet>( new PointSet( *this ) );
268}
269
270bool PointSet::containsPoint( double x, double y ) const
271{
272 GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
273 try
274 {
275 geos::unique_ptr point( GEOSGeom_createPointFromXY_r( geosctxt, x, y ) );
276 const bool result = ( GEOSPreparedContainsProperly_r( geosctxt, preparedGeom(), point.get() ) == 1 );
277
278 return result;
279 }
280 catch ( GEOSException &e )
281 {
282 qWarning( "GEOS exception: %s", e.what() );
283 QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
284 return false;
285 }
286
287}
288
289bool PointSet::containsLabelCandidate( double x, double y, double width, double height, double alpha ) const
290{
291 return GeomFunction::containsCandidate( preparedGeom(), x, y, width, height, alpha );
292}
293
294QLinkedList<PointSet *> PointSet::splitPolygons( PointSet *inputShape, double labelWidth, double labelHeight )
295{
296 int j;
297
298 double bestArea = 0;
299
300 double b;
301
302 int holeS = -1; // hole start and end points
303 int holeE = -1;
304
305 int retainedPt = -1;
306
307 const double labelArea = labelWidth * labelHeight;
308
309 QLinkedList<PointSet *> inputShapes;
310 inputShapes.push_back( inputShape );
311 QLinkedList<PointSet *> outputShapes;
312
313 while ( !inputShapes.isEmpty() )
314 {
315 PointSet *shape = inputShapes.takeFirst();
316
317 const std::vector< double > &x = shape->x;
318 const std::vector< double > &y = shape->y;
319 const int nbp = shape->nbPoints;
320 std::vector< int > pts( nbp );
321 for ( int i = 0; i < nbp; i++ )
322 {
323 pts[i] = i;
324 }
325
326 // compute convex hull
327 shape->convexHull = GeomFunction::convexHullId( pts, x, y );
328
329 bestArea = 0;
330 retainedPt = -1;
331
332 // lookup for a hole
333 for ( std::size_t ihs = 0; ihs < shape->convexHull.size(); ihs++ )
334 {
335 // ihs->ihn => cHull'seg
336 const std::size_t ihn = ( ihs + 1 ) % shape->convexHull.size();
337
338 const int ips = shape->convexHull[ihs];
339 const int ipn = ( ips + 1 ) % nbp;
340 if ( ipn != shape->convexHull[ihn] ) // next point on shape is not the next point on cHull => there is a hole here !
341 {
342 double bestcp = 0;
343 int pt = -1;
344 // lookup for the deepest point in the hole
345 for ( int i = ips; i != shape->convexHull[ihn]; i = ( i + 1 ) % nbp )
346 {
347 const double cp = std::fabs( GeomFunction::cross_product( x[shape->convexHull[ihs]], y[shape->convexHull[ihs]],
348 x[shape->convexHull[ihn]], y[shape->convexHull[ihn]],
349 x[i], y[i] ) );
350 if ( cp - bestcp > EPSILON )
351 {
352 bestcp = cp;
353 pt = i;
354 }
355 }
356
357 if ( pt != -1 )
358 {
359 // compute the ihs->ihn->pt triangle's area
360 const double base = GeomFunction::dist_euc2d( x[shape->convexHull[ihs]], y[shape->convexHull[ihs]],
361 x[shape->convexHull[ihn]], y[shape->convexHull[ihn]] );
362
363 b = GeomFunction::dist_euc2d( x[shape->convexHull[ihs]], y[shape->convexHull[ihs]],
364 x[pt], y[pt] );
365
366 const double c = GeomFunction::dist_euc2d( x[shape->convexHull[ihn]], y[shape->convexHull[ihn]],
367 x[pt], y[pt] );
368
369 const double s = ( base + b + c ) / 2; // s = half perimeter
370 double area = s * ( s - base ) * ( s - b ) * ( s - c );
371 if ( area < 0 )
372 area = -area;
373
374 // retain the biggest area
375 if ( area - bestArea > EPSILON )
376 {
377 bestArea = area;
378 retainedPt = pt;
379 holeS = ihs;
380 holeE = ihn;
381 }
382 }
383 }
384 }
385
386 // we have a hole, its area, and the deppest point in hole
387 // we're going to find the second point to cup the shape
388 // holeS = hole starting point
389 // holeE = hole ending point
390 // retainedPt = deppest point in hole
391 // bestArea = area of triangle HoleS->holeE->retainedPoint
392 bestArea = std::sqrt( bestArea );
393 double cx, cy, dx, dy, ex, ey, fx, fy, seg_length, ptx = 0, pty = 0, fptx = 0, fpty = 0;
394 int ps = -1, pe = -1, fps = -1, fpe = -1;
395 if ( retainedPt >= 0 && bestArea > labelArea ) // there is a hole so we'll cut the shape in two new shape (only if hole area is bigger than twice labelArea)
396 {
397 double c = std::numeric_limits<double>::max();
398
399 // iterate on all shape points except points which are in the hole
400 bool isValid;
401 int k, l;
402 for ( int i = ( shape->convexHull[holeE] + 1 ) % nbp; i != ( shape->convexHull[holeS] - 1 + nbp ) % nbp; i = j )
403 {
404 j = ( i + 1 ) % nbp; // i->j is shape segment not in hole
405
406 // compute distance between retainedPoint and segment
407 // whether perpendicular distance (if retaindPoint is fronting segment i->j)
408 // or distance between retainedPt and i or j (choose the nearest)
409 seg_length = GeomFunction::dist_euc2d( x[i], y[i], x[j], y[j] );
410 cx = ( x[i] + x[j] ) / 2.0;
411 cy = ( y[i] + y[j] ) / 2.0;
412 dx = cy - y[i];
413 dy = cx - x[i];
414
415 ex = cx - dx;
416 ey = cy + dy;
417 fx = cx + dx;
418 fy = cy - dy;
419
420 if ( seg_length < EPSILON || std::fabs( ( b = GeomFunction::cross_product( ex, ey, fx, fy, x[retainedPt], y[retainedPt] ) / ( seg_length ) ) ) > ( seg_length / 2 ) ) // retainedPt is not fronting i->j
421 {
422 if ( ( ex = GeomFunction::dist_euc2d_sq( x[i], y[i], x[retainedPt], y[retainedPt] ) ) < ( ey = GeomFunction::dist_euc2d_sq( x[j], y[j], x[retainedPt], y[retainedPt] ) ) )
423 {
424 b = ex;
425 ps = i;
426 pe = i;
427 }
428 else
429 {
430 b = ey;
431 ps = j;
432 pe = j;
433 }
434 }
435 else // point fronting i->j => compute pependicular distance => create a new point
436 {
437 b = GeomFunction::cross_product( x[i], y[i], x[j], y[j], x[retainedPt], y[retainedPt] ) / seg_length;
438 b *= b;
439 ps = i;
440 pe = j;
441
442 if ( !GeomFunction::computeLineIntersection( x[i], y[i], x[j], y[j], x[retainedPt], y[retainedPt], x[retainedPt] - dx, y[retainedPt] + dy, &ptx, &pty ) )
443 {
444 //error - it should intersect the line
445 }
446 }
447
448 isValid = true;
449 double pointX, pointY;
450 if ( ps == pe )
451 {
452 pointX = x[pe];
453 pointY = y[pe];
454 }
455 else
456 {
457 pointX = ptx;
458 pointY = pty;
459 }
460
461 for ( k = shape->convexHull[holeS]; k != shape->convexHull[holeE]; k = ( k + 1 ) % nbp )
462 {
463 l = ( k + 1 ) % nbp;
464 if ( GeomFunction::isSegIntersects( x[retainedPt], y[retainedPt], pointX, pointY, x[k], y[k], x[l], y[l] ) )
465 {
466 isValid = false;
467 break;
468 }
469 }
470
471
472 if ( isValid && b < c )
473 {
474 c = b;
475 fps = ps;
476 fpe = pe;
477 fptx = ptx;
478 fpty = pty;
479 }
480 } // for point which are not in hole
481
482 // we will cut the shapeu in two new shapes, one from [retainedPoint] to [newPoint] and one form [newPoint] to [retainedPoint]
483 const int imin = retainedPt;
484 int imax = ( ( ( fps < retainedPt && fpe < retainedPt ) || ( fps > retainedPt && fpe > retainedPt ) ) ? std::min( fps, fpe ) : std::max( fps, fpe ) );
485
486 int nbPtSh1, nbPtSh2; // how many points in new shapes ?
487 if ( imax > imin )
488 nbPtSh1 = imax - imin + 1 + ( fpe != fps );
489 else
490 nbPtSh1 = imax + nbp - imin + 1 + ( fpe != fps );
491
492 if ( ( imax == fps ? fpe : fps ) < imin )
493 nbPtSh2 = imin - ( imax == fps ? fpe : fps ) + 1 + ( fpe != fps );
494 else
495 nbPtSh2 = imin + nbp - ( imax == fps ? fpe : fps ) + 1 + ( fpe != fps );
496
497 if ( retainedPt == -1 || fps == -1 || fpe == -1 )
498 {
499 if ( shape->parent )
500 delete shape;
501 }
502 // check for useless splitting
503 else if ( imax == imin || nbPtSh1 <= 2 || nbPtSh2 <= 2 || nbPtSh1 == nbp || nbPtSh2 == nbp )
504 {
505 outputShapes.append( shape );
506 }
507 else
508 {
509
510 PointSet *newShape = shape->extractShape( nbPtSh1, imin, imax, fps, fpe, fptx, fpty ).release();
511
512 if ( shape->parent )
513 newShape->parent = shape->parent;
514 else
515 newShape->parent = shape;
516
517 inputShapes.append( newShape );
518
519 if ( imax == fps )
520 imax = fpe;
521 else
522 imax = fps;
523
524 newShape = shape->extractShape( nbPtSh2, imax, imin, fps, fpe, fptx, fpty ).release();
525
526 if ( shape->parent )
527 newShape->parent = shape->parent;
528 else
529 newShape->parent = shape;
530
531 inputShapes.append( newShape );
532
533 if ( shape->parent )
534 delete shape;
535 }
536 }
537 else
538 {
539 outputShapes.append( shape );
540 }
541 }
542 return outputShapes;
543}
544
545void PointSet::offsetCurveByDistance( double distance )
546{
547 if ( !mGeos )
549
550 if ( !mGeos || type != GEOS_LINESTRING )
551 return;
552
553 GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
554 geos::unique_ptr newGeos = nullptr;
555 try
556 {
557 newGeos.reset( GEOSOffsetCurve_r( geosctxt, mGeos, distance, 0, GEOSBUF_JOIN_MITRE, 2 ) );
558 if ( !newGeos )
559 return;
560
561 // happens sometime, if the offset curve self-intersects
562 if ( GEOSGeomTypeId_r( geosctxt, newGeos.get() ) == GEOS_MULTILINESTRING )
563 {
564 // we keep the longest part
565 const int nParts = GEOSGetNumGeometries_r( geosctxt, newGeos.get() );
566 double maximumLength = -1;
567 const GEOSGeometry *longestPart = nullptr;
568 for ( int i = 0; i < nParts; ++i )
569 {
570 const GEOSGeometry *part = GEOSGetGeometryN_r( geosctxt, newGeos.get(), i );
571 double partLength = -1;
572 if ( GEOSLength_r( geosctxt, part, &partLength ) == 1 )
573 {
574 if ( partLength > maximumLength )
575 {
576 maximumLength = partLength;
577 longestPart = part;
578 }
579 }
580 }
581
582 if ( !longestPart )
583 {
584 // something is really wrong!
585 return;
586 }
587
588 geos::unique_ptr longestPartClone( GEOSGeom_clone_r( geosctxt, longestPart ) );
589 newGeos = std::move( longestPartClone );
590 }
591
592#if GEOS_VERSION_MAJOR==3 && GEOS_VERSION_MINOR<11
593 if ( distance < 0 )
594 {
595 // geos < 3.11 reverses the direction of offset curves with negative distances -- we don't want that!
596 geos::unique_ptr reversed( GEOSReverse_r( geosctxt, newGeos.get() ) );
597 if ( !reversed )
598 return;
599
600 newGeos = std::move( reversed );
601 }
602#endif
603
604 const int newNbPoints = GEOSGeomGetNumPoints_r( geosctxt, newGeos.get() );
605 const GEOSCoordSequence *coordSeq = GEOSGeom_getCoordSeq_r( geosctxt, newGeos.get() );
606 std::vector< double > newX;
607 std::vector< double > newY;
608 newX.resize( newNbPoints );
609 newY.resize( newNbPoints );
610 for ( int i = 0; i < newNbPoints; i++ )
611 {
612 GEOSCoordSeq_getX_r( geosctxt, coordSeq, i, &newX[i] );
613 GEOSCoordSeq_getY_r( geosctxt, coordSeq, i, &newY[i] );
614 }
615 nbPoints = newNbPoints;
616 x = newX;
617 y = newY;
618 }
619 catch ( GEOSException &e )
620 {
621 qWarning( "GEOS exception: %s", e.what() );
622 QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
623 return;
624 }
625
627 mGeos = newGeos.release();
628 mOwnsGeom = true;
629}
630
631void PointSet::extendLineByDistance( double startDistance, double endDistance, double smoothDistance )
632{
633 if ( nbPoints < 2 )
634 return;
635
636 double x0 = x[0];
637 double y0 = y[0];
638 if ( startDistance > 0 )
639 {
640 // trace forward by smoothDistance
641 double x1 = x[1];
642 double y1 = y[1];
643
644 double distanceConsumed = 0;
645 double lastX = x0;
646 double lastY = y0;
647 for ( int i = 1; i < nbPoints; ++i )
648 {
649 const double thisX = x[i];
650 const double thisY = y[i];
651 const double thisSegmentLength = std::sqrt( ( thisX - lastX ) * ( thisX - lastX ) + ( thisY - lastY ) * ( thisY - lastY ) );
652 distanceConsumed += thisSegmentLength;
653 if ( distanceConsumed >= smoothDistance )
654 {
655 const double c = ( distanceConsumed - smoothDistance ) / thisSegmentLength;
656 x1 = lastX + c * ( thisX - lastX );
657 y1 = lastY + c * ( thisY - lastY );
658 break;
659 }
660 lastX = thisX;
661 lastY = thisY;
662 }
663
664 const double distance = std::sqrt( ( x1 - x0 ) * ( x1 - x0 ) + ( y1 - y0 ) * ( y1 - y0 ) );
665 const double extensionFactor = ( startDistance + distance ) / distance;
666 const QgsPointXY newStart = QgsGeometryUtils::interpolatePointOnLine( x1, y1, x0, y0, extensionFactor );
667 x0 = newStart.x();
668 y0 = newStart.y();
669 // defer actually changing the stored start until we've calculated the new end point
670 }
671
672 if ( endDistance > 0 )
673 {
674 const double xend0 = x[nbPoints - 1];
675 const double yend0 = y[nbPoints - 1];
676 double xend1 = x[nbPoints - 2];
677 double yend1 = y[nbPoints - 2];
678
679 // trace backward by smoothDistance
680 double distanceConsumed = 0;
681 double lastX = x0;
682 double lastY = y0;
683 for ( int i = nbPoints - 2; i >= 0; --i )
684 {
685 const double thisX = x[i];
686 const double thisY = y[i];
687 const double thisSegmentLength = std::sqrt( ( thisX - lastX ) * ( thisX - lastX ) + ( thisY - lastY ) * ( thisY - lastY ) );
688 distanceConsumed += thisSegmentLength;
689 if ( distanceConsumed >= smoothDistance )
690 {
691 const double c = ( distanceConsumed - smoothDistance ) / thisSegmentLength;
692 xend1 = lastX + c * ( thisX - lastX );
693 yend1 = lastY + c * ( thisY - lastY );
694 break;
695 }
696 lastX = thisX;
697 lastY = thisY;
698 }
699
700 const double distance = std::sqrt( ( xend1 - xend0 ) * ( xend1 - xend0 ) + ( yend1 - yend0 ) * ( yend1 - yend0 ) );
701 const double extensionFactor = ( endDistance + distance ) / distance;
702 const QgsPointXY newEnd = QgsGeometryUtils::interpolatePointOnLine( xend1, yend1, xend0, yend0, extensionFactor );
703 x.emplace_back( newEnd.x() );
704 y.emplace_back( newEnd.y() );
705 nbPoints++;
706 }
707
708 if ( startDistance > 0 )
709 {
710 x.insert( x.begin(), x0 );
711 y.insert( y.begin(), y0 );
712 nbPoints++;
713 }
714
716}
717
719{
720 ok = false;
721 double bbox[4]; // xmin, ymin, xmax, ymax
722
723 double alpha;
724 int alpha_d;
725
726 double alpha_seg;
727
728 double d1, d2;
729
730 double bb[16]; // {ax, ay, bx, by, cx, cy, dx, dy, ex, ey, fx, fy, gx, gy, hx, hy}}
731
732 double best_area = std::numeric_limits<double>::max();
733 double best_alpha = -1;
734 double best_bb[16];
735 double best_length = 0;
736 double best_width = 0;
737
738
739 bbox[0] = std::numeric_limits<double>::max();
740 bbox[1] = std::numeric_limits<double>::max();
741 bbox[2] = std::numeric_limits<double>::lowest();
742 bbox[3] = std::numeric_limits<double>::lowest();
743
744 for ( std::size_t i = 0; i < convexHull.size(); i++ )
745 {
746 if ( x[convexHull[i]] < bbox[0] )
747 bbox[0] = x[convexHull[i]];
748
749 if ( x[convexHull[i]] > bbox[2] )
750 bbox[2] = x[convexHull[i]];
751
752 if ( y[convexHull[i]] < bbox[1] )
753 bbox[1] = y[convexHull[i]];
754
755 if ( y[convexHull[i]] > bbox[3] )
756 bbox[3] = y[convexHull[i]];
757 }
758
760
761 const double dref = bbox[2] - bbox[0];
762 if ( qgsDoubleNear( dref, 0 ) )
763 {
764 ok = false;
765 return finalBb;
766 }
767
768 for ( alpha_d = 0; alpha_d < 90; alpha_d++ )
769 {
770 alpha = alpha_d * M_PI / 180.0;
771 d1 = std::cos( alpha ) * dref;
772 d2 = std::sin( alpha ) * dref;
773
774 bb[0] = bbox[0];
775 bb[1] = bbox[3]; // ax, ay
776
777 bb[4] = bbox[0];
778 bb[5] = bbox[1]; // cx, cy
779
780 bb[8] = bbox[2];
781 bb[9] = bbox[1]; // ex, ey
782
783 bb[12] = bbox[2];
784 bb[13] = bbox[3]; // gx, gy
785
786
787 bb[2] = bb[0] + d1;
788 bb[3] = bb[1] + d2; // bx, by
789 bb[6] = bb[4] - d2;
790 bb[7] = bb[5] + d1; // dx, dy
791 bb[10] = bb[8] - d1;
792 bb[11] = bb[9] - d2; // fx, fy
793 bb[14] = bb[12] + d2;
794 bb[15] = bb[13] - d1; // hx, hy
795
796 // adjust all points
797 for ( int i = 0; i < 16; i += 4 )
798 {
799
800 alpha_seg = ( ( i / 4 > 0 ? ( i / 4 ) - 1 : 3 ) ) * M_PI_2 + alpha;
801
802 double best_cp = std::numeric_limits<double>::max();
803
804 for ( std::size_t j = 0; j < convexHull.size(); j++ )
805 {
806 const double cp = GeomFunction::cross_product( bb[i + 2], bb[i + 3], bb[i], bb[i + 1], x[convexHull[j]], y[convexHull[j]] );
807 if ( cp < best_cp )
808 {
809 best_cp = cp;
810 }
811 }
812
813 const double distNearestPoint = best_cp / dref;
814
815 d1 = std::cos( alpha_seg ) * distNearestPoint;
816 d2 = std::sin( alpha_seg ) * distNearestPoint;
817
818 bb[i] += d1; // x
819 bb[i + 1] += d2; // y
820 bb[i + 2] += d1; // x
821 bb[i + 3] += d2; // y
822 }
823
824 // compute and compare AREA
825 const double width = GeomFunction::cross_product( bb[6], bb[7], bb[4], bb[5], bb[12], bb[13] ) / dref;
826 const double length = GeomFunction::cross_product( bb[2], bb[3], bb[0], bb[1], bb[8], bb[9] ) / dref;
827
828 double area = width * length;
829
830 if ( area < 0 )
831 area *= -1;
832
833
834 if ( best_area - area > EPSILON )
835 {
836 best_area = area;
837 best_length = length;
838 best_width = width;
839 best_alpha = alpha;
840 memcpy( best_bb, bb, sizeof( double ) * 16 );
841 }
842 }
843
844 // best bbox is defined
845 for ( int i = 0; i < 16; i = i + 4 )
846 {
847 GeomFunction::computeLineIntersection( best_bb[i], best_bb[i + 1], best_bb[i + 2], best_bb[i + 3],
848 best_bb[( i + 4 ) % 16], best_bb[( i + 5 ) % 16], best_bb[( i + 6 ) % 16], best_bb[( i + 7 ) % 16],
849 &finalBb.x[int ( i / 4 )], &finalBb.y[int ( i / 4 )] );
850 }
851
852 finalBb.alpha = best_alpha;
853 finalBb.width = best_width;
854 finalBb.length = best_length;
855
856 ok = true;
857 return finalBb;
858}
859
860double PointSet::minDistanceToPoint( double px, double py, double *rx, double *ry ) const
861{
862 if ( !mGeos )
864
865 if ( !mGeos )
866 return 0;
867
868 GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
869 try
870 {
871 geos::unique_ptr geosPt( GEOSGeom_createPointFromXY_r( geosctxt, px, py ) );
872 const int type = GEOSGeomTypeId_r( geosctxt, mGeos );
873 const GEOSGeometry *extRing = nullptr;
874 const GEOSPreparedGeometry *preparedExtRing = nullptr;
875
876 if ( type != GEOS_POLYGON )
877 {
878 extRing = mGeos;
879 preparedExtRing = preparedGeom();
880 }
881 else
882 {
883 //for polygons, we want distance to exterior ring (not an interior point)
884 extRing = GEOSGetExteriorRing_r( geosctxt, mGeos );
885 if ( ! mGeosPreparedBoundary )
886 {
887 mGeosPreparedBoundary = GEOSPrepare_r( geosctxt, extRing );
888 }
889 preparedExtRing = mGeosPreparedBoundary;
890 }
891
892 const geos::coord_sequence_unique_ptr nearestCoord( GEOSPreparedNearestPoints_r( geosctxt, preparedExtRing, geosPt.get() ) );
893 double nx;
894 double ny;
895 unsigned int nPoints = 0;
896 GEOSCoordSeq_getSize_r( geosctxt, nearestCoord.get(), &nPoints );
897 if ( nPoints == 0 )
898 return 0;
899
900 ( void )GEOSCoordSeq_getXY_r( geosctxt, nearestCoord.get(), 0, &nx, &ny );
901
902 if ( rx )
903 *rx = nx;
904 if ( ry )
905 *ry = ny;
906
907 return GeomFunction::dist_euc2d_sq( px, py, nx, ny );
908 }
909 catch ( GEOSException &e )
910 {
911 qWarning( "GEOS exception: %s", e.what() );
912 QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
913 return 0;
914 }
915}
916
917void PointSet::getCentroid( double &px, double &py, bool forceInside ) const
918{
919 if ( !mGeos )
921
922 if ( !mGeos )
923 return;
924
925 try
926 {
927 GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
928 geos::unique_ptr centroidGeom( GEOSGetCentroid_r( geosctxt, mGeos ) );
929 if ( centroidGeom )
930 {
931 const GEOSCoordSequence *coordSeq = GEOSGeom_getCoordSeq_r( geosctxt, centroidGeom.get() );
932 unsigned int nPoints = 0;
933 GEOSCoordSeq_getSize_r( geosctxt, coordSeq, &nPoints );
934 if ( nPoints == 0 )
935 return;
936 GEOSCoordSeq_getXY_r( geosctxt, coordSeq, 0, &px, &py );
937 }
938
939 // check if centroid inside in polygon
940 if ( forceInside && !containsPoint( px, py ) )
941 {
942 geos::unique_ptr pointGeom( GEOSPointOnSurface_r( geosctxt, mGeos ) );
943
944 if ( pointGeom )
945 {
946 const GEOSCoordSequence *coordSeq = GEOSGeom_getCoordSeq_r( geosctxt, pointGeom.get() );
947 unsigned int nPoints = 0;
948 GEOSCoordSeq_getSize_r( geosctxt, coordSeq, &nPoints );
949 if ( nPoints == 0 )
950 return;
951
952 GEOSCoordSeq_getXY_r( geosctxt, coordSeq, 0, &px, &py );
953 }
954 }
955 }
956 catch ( GEOSException &e )
957 {
958 qWarning( "GEOS exception: %s", e.what() );
959 QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
960 return;
961 }
962}
963
965{
966 const double x1 = ( xmin > other->xmin ? xmin : other->xmin );
967 const double x2 = ( xmax < other->xmax ? xmax : other->xmax );
968 if ( x1 > x2 )
969 return false;
970 const double y1 = ( ymin > other->ymin ? ymin : other->ymin );
971 const double y2 = ( ymax < other->ymax ? ymax : other->ymax );
972 return y1 <= y2;
973}
974
975void PointSet::getPointByDistance( double *d, double *ad, double dl, double *px, double *py ) const
976{
977 int i;
978 double dx, dy, di;
979 double distr;
980
981 i = 0;
982 if ( dl >= 0 )
983 {
984 while ( i < nbPoints && ad[i] <= dl ) i++;
985 i--;
986 }
987
988 if ( i < nbPoints - 1 )
989 {
990 if ( dl < 0 )
991 {
992 dx = x[nbPoints - 1] - x[0];
993 dy = y[nbPoints - 1] - y[0];
994 di = std::sqrt( dx * dx + dy * dy );
995 }
996 else
997 {
998 dx = x[i + 1] - x[i];
999 dy = y[i + 1] - y[i];
1000 di = d[i];
1001 }
1002
1003 distr = dl - ad[i];
1004 *px = x[i] + dx * distr / di;
1005 *py = y[i] + dy * distr / di;
1006 }
1007 else // just select last point...
1008 {
1009 *px = x[i];
1010 *py = y[i];
1011 }
1012}
1013
1015{
1016 const GEOSGeometry *thisGeos = geos();
1017 if ( !thisGeos )
1018 return nullptr;
1019
1020 try
1021 {
1022 geos::unique_ptr res( GEOSInterpolate_r( QgsGeos::getGEOSHandler(), thisGeos, distance ) );
1023 return res;
1024 }
1025 catch ( GEOSException &e )
1026 {
1027 qWarning( "GEOS exception: %s", e.what() );
1028 return nullptr;
1029 }
1030}
1031
1032double PointSet::lineLocatePoint( const GEOSGeometry *point ) const
1033{
1034 const GEOSGeometry *thisGeos = geos();
1035 if ( !thisGeos )
1036 return -1;
1037
1038 double distance = -1;
1039 try
1040 {
1041 distance = GEOSProject_r( QgsGeos::getGEOSHandler(), thisGeos, point );
1042 }
1043 catch ( GEOSException &e )
1044 {
1045 qWarning( "GEOS exception: %s", e.what() );
1046 return -1;
1047 }
1048
1049 return distance;
1050}
1051
1052const GEOSGeometry *PointSet::geos() const
1053{
1054 if ( !mGeos )
1056
1057 return mGeos;
1058}
1059
1060double PointSet::length() const
1061{
1062 if ( mLength >= 0 )
1063 return mLength;
1064
1065 if ( !mGeos )
1067
1068 if ( !mGeos )
1069 return -1;
1070
1071 GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
1072
1073 try
1074 {
1075 ( void )GEOSLength_r( geosctxt, mGeos, &mLength );
1076 return mLength;
1077 }
1078 catch ( GEOSException &e )
1079 {
1080 qWarning( "GEOS exception: %s", e.what() );
1081 QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
1082 return -1;
1083 }
1084}
1085
1086double PointSet::area() const
1087{
1088 if ( mArea >= 0 )
1089 return mArea;
1090
1091 if ( !mGeos )
1093
1094 if ( !mGeos )
1095 return -1;
1096
1097 GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
1098
1099 try
1100 {
1101 ( void )GEOSArea_r( geosctxt, mGeos, &mArea );
1102 mArea = std::fabs( mArea );
1103 return mArea;
1104 }
1105 catch ( GEOSException &e )
1106 {
1107 qWarning( "GEOS exception: %s", e.what() );
1108 QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
1109 return -1;
1110 }
1111}
1112
1114{
1115 return qgsDoubleNear( x[0], x[nbPoints - 1] ) && qgsDoubleNear( y[0], y[nbPoints - 1] );
1116}
1117
1118QString PointSet::toWkt() const
1119{
1120 if ( !mGeos )
1122
1123 GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
1124
1125 try
1126 {
1127 GEOSWKTWriter *writer = GEOSWKTWriter_create_r( geosctxt );
1128
1129 char *wkt = GEOSWKTWriter_write_r( geosctxt, writer, mGeos );
1130 const QString res( wkt );
1131
1132 GEOSFree_r( geosctxt, wkt );
1133
1134 GEOSWKTWriter_destroy_r( geosctxt, writer );
1135 writer = nullptr;
1136
1137 return res;
1138 }
1139 catch ( GEOSException &e )
1140 {
1141 qWarning( "GEOS exception: %s", e.what() );
1142 QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
1143 return QString();
1144 }
1145}
1146
1147std::tuple< std::vector< double >, double > PointSet::edgeDistances() const
1148{
1149 std::vector< double > distances( nbPoints );
1150 double totalDistance = 0;
1151 double oldX = -1.0, oldY = -1.0;
1152 for ( int i = 0; i < nbPoints; i++ )
1153 {
1154 if ( i == 0 )
1155 distances[i] = 0;
1156 else
1157 distances[i] = std::sqrt( std::pow( oldX - x[i], 2 ) + std::pow( oldY - y[i], 2 ) );
1158
1159 oldX = x[i];
1160 oldY = y[i];
1161 totalDistance += distances[i];
1162 }
1163 return std::make_tuple( std::move( distances ), totalDistance );
1164}
static QgsPointXY interpolatePointOnLine(double x1, double y1, double x2, double y2, double fraction) SIP_HOLDGIL
Interpolates the position of a point a fraction of the way along the line from (x1,...
static GEOSContextHandle_t getGEOSHandler()
Definition: qgsgeos.cpp:3446
static void logMessage(const QString &message, const QString &tag=QString(), Qgis::MessageLevel level=Qgis::MessageLevel::Warning, bool notifyUser=true)
Adds a message to the log instance (and creates it if necessary).
A class to represent a 2D point.
Definition: qgspointxy.h:59
double y
Definition: qgspointxy.h:63
Q_GADGET double x
Definition: qgspointxy.h:62
static double dist_euc2d(double x1, double y1, double x2, double y2)
Definition: geomfunction.h:67
static std::vector< int > convexHullId(std::vector< int > &id, const std::vector< double > &x, const std::vector< double > &y)
Compute the convex hull in O(n┬Ělog(n))
static bool computeLineIntersection(double x1, double y1, double x2, double y2, double x3, double y3, double x4, double y4, double *x, double *y)
Compute the point where two lines intersect.
static bool isSegIntersects(double x1, double y1, double x2, double y2, double x3, double y3, double x4, double y4)
Returns true if the two segments intersect.
static double cross_product(double x1, double y1, double x2, double y2, double x3, double y3)
Definition: geomfunction.h:62
static double dist_euc2d_sq(double x1, double y1, double x2, double y2)
Definition: geomfunction.h:72
static bool containsCandidate(const GEOSPreparedGeometry *geom, double x, double y, double width, double height, double alpha)
Returns true if a GEOS prepared geometry totally contains a label candidate.
The underlying raw pal geometry class.
Definition: pointset.h:77
geos::unique_ptr interpolatePoint(double distance) const
Returns a GEOS geometry representing the point interpolated on the shape by distance.
Definition: pointset.cpp:1014
std::unique_ptr< PointSet > clone() const
Returns a copy of the point set.
Definition: pointset.cpp:265
bool containsLabelCandidate(double x, double y, double width, double height, double alpha=0) const
Tests whether a possible label candidate will fit completely within the shape.
Definition: pointset.cpp:289
double lineLocatePoint(const GEOSGeometry *point) const
Returns the distance along the geometry closest to the specified GEOS point.
Definition: pointset.cpp:1032
OrientedConvexHullBoundingBox computeConvexHullOrientedBoundingBox(bool &ok)
Computes an oriented bounding box for the shape's convex hull.
Definition: pointset.cpp:718
double length() const
Returns length of line geometry.
Definition: pointset.cpp:1060
void deleteCoords()
Definition: pointset.cpp:232
void extendLineByDistance(double startDistance, double endDistance, double smoothDistance)
Extends linestrings by the specified amount at the start and end of the line, by extending the existi...
Definition: pointset.cpp:631
virtual ~PointSet()
Definition: pointset.cpp:201
double ymax
Definition: pointset.h:261
double ymin
Definition: pointset.h:260
double area() const
Returns area of polygon geometry.
Definition: pointset.cpp:1086
bool isClosed() const
Returns true if pointset is closed.
Definition: pointset.cpp:1113
PointSet * holeOf
Definition: pointset.h:241
void createGeosGeom() const
Definition: pointset.cpp:99
void getPointByDistance(double *d, double *ad, double dl, double *px, double *py) const
Gets a point a set distance along a line geometry.
Definition: pointset.cpp:975
std::vector< double > y
Definition: pointset.h:231
bool boundingBoxIntersects(const PointSet *other) const
Returns true if the bounding box of this pointset intersects the bounding box of another pointset.
Definition: pointset.cpp:964
void getCentroid(double &px, double &py, bool forceInside=false) const
Definition: pointset.cpp:917
std::vector< double > x
Definition: pointset.h:230
void offsetCurveByDistance(double distance)
Offsets linestrings by the specified distance.
Definition: pointset.cpp:545
std::vector< int > convexHull
Definition: pointset.h:237
const GEOSPreparedGeometry * preparedGeom() const
Definition: pointset.cpp:154
QString toWkt() const
Returns a WKT representation of the point set.
Definition: pointset.cpp:1118
GEOSGeometry * mGeos
Definition: pointset.h:234
double xmin
Definition: pointset.h:258
const GEOSGeometry * geos() const
Returns the point set's GEOS geometry.
Definition: pointset.cpp:1052
void invalidateGeos() const
Definition: pointset.cpp:166
double minDistanceToPoint(double px, double py, double *rx=nullptr, double *ry=nullptr) const
Returns the squared minimum distance between the point set geometry and the point (px,...
Definition: pointset.cpp:860
double xmax
Definition: pointset.h:259
std::unique_ptr< PointSet > extractShape(int nbPtSh, int imin, int imax, int fps, int fpe, double fptx, double fpty)
Does... something completely inscrutable.
Definition: pointset.cpp:238
bool containsPoint(double x, double y) const
Tests whether point set contains a specified point.
Definition: pointset.cpp:270
std::tuple< std::vector< double >, double > edgeDistances() const
Returns a vector of edge distances as well as its total length.
Definition: pointset.cpp:1147
PointSet * parent
Definition: pointset.h:242
double mLength
Definition: pointset.h:245
double mArea
Definition: pointset.h:244
bool mOwnsGeom
Definition: pointset.h:235
static QLinkedList< PointSet * > splitPolygons(PointSet *inputShape, double labelWidth, double labelHeight)
Split a polygon using some random logic into some other polygons.
Definition: pointset.cpp:294
std::unique_ptr< GEOSGeometry, GeosDeleter > unique_ptr
Scoped GEOS pointer.
Definition: qgsgeos.h:74
std::unique_ptr< GEOSCoordSequence, GeosDeleter > coord_sequence_unique_ptr
Scoped GEOS coordinate sequence pointer.
Definition: qgsgeos.h:89
As part of the API refactoring and improvements which landed in the Processing API was substantially reworked from the x version This was done in order to allow much of the underlying Processing framework to be ported into c
bool qgsDoubleNear(double a, double b, double epsilon=4 *std::numeric_limits< double >::epsilon())
Compare two doubles (but allow some difference)
Definition: qgis.h:3090
Represents the minimum area, oriented bounding box surrounding a convex hull.
Definition: pointset.h:60
#define EPSILON
Definition: util.h:78