QGIS API Documentation 3.28.0-Firenze (ed3ad0430f)
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 ( distance < 0 )
593 {
594 // geos reverses the direction of offset curves with negative distances -- we don't want that!
595 geos::unique_ptr reversed( GEOSReverse_r( geosctxt, newGeos.get() ) );
596 if ( !reversed )
597 return;
598
599 newGeos = std::move( reversed );
600 }
601
602 const int newNbPoints = GEOSGeomGetNumPoints_r( geosctxt, newGeos.get() );
603 const GEOSCoordSequence *coordSeq = GEOSGeom_getCoordSeq_r( geosctxt, newGeos.get() );
604 std::vector< double > newX;
605 std::vector< double > newY;
606 newX.resize( newNbPoints );
607 newY.resize( newNbPoints );
608 for ( int i = 0; i < newNbPoints; i++ )
609 {
610 GEOSCoordSeq_getX_r( geosctxt, coordSeq, i, &newX[i] );
611 GEOSCoordSeq_getY_r( geosctxt, coordSeq, i, &newY[i] );
612 }
613 nbPoints = newNbPoints;
614 x = newX;
615 y = newY;
616 }
617 catch ( GEOSException &e )
618 {
619 qWarning( "GEOS exception: %s", e.what() );
620 QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
621 return;
622 }
623
625 mGeos = newGeos.release();
626 mOwnsGeom = true;
627}
628
629void PointSet::extendLineByDistance( double startDistance, double endDistance, double smoothDistance )
630{
631 if ( nbPoints < 2 )
632 return;
633
634 double x0 = x[0];
635 double y0 = y[0];
636 if ( startDistance > 0 )
637 {
638 // trace forward by smoothDistance
639 double x1 = x[1];
640 double y1 = y[1];
641
642 double distanceConsumed = 0;
643 double lastX = x0;
644 double lastY = y0;
645 for ( int i = 1; i < nbPoints; ++i )
646 {
647 const double thisX = x[i];
648 const double thisY = y[i];
649 const double thisSegmentLength = std::sqrt( ( thisX - lastX ) * ( thisX - lastX ) + ( thisY - lastY ) * ( thisY - lastY ) );
650 distanceConsumed += thisSegmentLength;
651 if ( distanceConsumed >= smoothDistance )
652 {
653 const double c = ( distanceConsumed - smoothDistance ) / thisSegmentLength;
654 x1 = lastX + c * ( thisX - lastX );
655 y1 = lastY + c * ( thisY - lastY );
656 break;
657 }
658 lastX = thisX;
659 lastY = thisY;
660 }
661
662 const double distance = std::sqrt( ( x1 - x0 ) * ( x1 - x0 ) + ( y1 - y0 ) * ( y1 - y0 ) );
663 const double extensionFactor = ( startDistance + distance ) / distance;
664 const QgsPointXY newStart = QgsGeometryUtils::interpolatePointOnLine( x1, y1, x0, y0, extensionFactor );
665 x0 = newStart.x();
666 y0 = newStart.y();
667 // defer actually changing the stored start until we've calculated the new end point
668 }
669
670 if ( endDistance > 0 )
671 {
672 const double xend0 = x[nbPoints - 1];
673 const double yend0 = y[nbPoints - 1];
674 double xend1 = x[nbPoints - 2];
675 double yend1 = y[nbPoints - 2];
676
677 // trace backward by smoothDistance
678 double distanceConsumed = 0;
679 double lastX = x0;
680 double lastY = y0;
681 for ( int i = nbPoints - 2; i >= 0; --i )
682 {
683 const double thisX = x[i];
684 const double thisY = y[i];
685 const double thisSegmentLength = std::sqrt( ( thisX - lastX ) * ( thisX - lastX ) + ( thisY - lastY ) * ( thisY - lastY ) );
686 distanceConsumed += thisSegmentLength;
687 if ( distanceConsumed >= smoothDistance )
688 {
689 const double c = ( distanceConsumed - smoothDistance ) / thisSegmentLength;
690 xend1 = lastX + c * ( thisX - lastX );
691 yend1 = lastY + c * ( thisY - lastY );
692 break;
693 }
694 lastX = thisX;
695 lastY = thisY;
696 }
697
698 const double distance = std::sqrt( ( xend1 - xend0 ) * ( xend1 - xend0 ) + ( yend1 - yend0 ) * ( yend1 - yend0 ) );
699 const double extensionFactor = ( endDistance + distance ) / distance;
700 const QgsPointXY newEnd = QgsGeometryUtils::interpolatePointOnLine( xend1, yend1, xend0, yend0, extensionFactor );
701 x.emplace_back( newEnd.x() );
702 y.emplace_back( newEnd.y() );
703 nbPoints++;
704 }
705
706 if ( startDistance > 0 )
707 {
708 x.insert( x.begin(), x0 );
709 y.insert( y.begin(), y0 );
710 nbPoints++;
711 }
712
714}
715
717{
718 ok = false;
719 double bbox[4]; // xmin, ymin, xmax, ymax
720
721 double alpha;
722 int alpha_d;
723
724 double alpha_seg;
725
726 double d1, d2;
727
728 double bb[16]; // {ax, ay, bx, by, cx, cy, dx, dy, ex, ey, fx, fy, gx, gy, hx, hy}}
729
730 double best_area = std::numeric_limits<double>::max();
731 double best_alpha = -1;
732 double best_bb[16];
733 double best_length = 0;
734 double best_width = 0;
735
736
737 bbox[0] = std::numeric_limits<double>::max();
738 bbox[1] = std::numeric_limits<double>::max();
739 bbox[2] = std::numeric_limits<double>::lowest();
740 bbox[3] = std::numeric_limits<double>::lowest();
741
742 for ( std::size_t i = 0; i < convexHull.size(); i++ )
743 {
744 if ( x[convexHull[i]] < bbox[0] )
745 bbox[0] = x[convexHull[i]];
746
747 if ( x[convexHull[i]] > bbox[2] )
748 bbox[2] = x[convexHull[i]];
749
750 if ( y[convexHull[i]] < bbox[1] )
751 bbox[1] = y[convexHull[i]];
752
753 if ( y[convexHull[i]] > bbox[3] )
754 bbox[3] = y[convexHull[i]];
755 }
756
758
759 const double dref = bbox[2] - bbox[0];
760 if ( qgsDoubleNear( dref, 0 ) )
761 {
762 ok = false;
763 return finalBb;
764 }
765
766 for ( alpha_d = 0; alpha_d < 90; alpha_d++ )
767 {
768 alpha = alpha_d * M_PI / 180.0;
769 d1 = std::cos( alpha ) * dref;
770 d2 = std::sin( alpha ) * dref;
771
772 bb[0] = bbox[0];
773 bb[1] = bbox[3]; // ax, ay
774
775 bb[4] = bbox[0];
776 bb[5] = bbox[1]; // cx, cy
777
778 bb[8] = bbox[2];
779 bb[9] = bbox[1]; // ex, ey
780
781 bb[12] = bbox[2];
782 bb[13] = bbox[3]; // gx, gy
783
784
785 bb[2] = bb[0] + d1;
786 bb[3] = bb[1] + d2; // bx, by
787 bb[6] = bb[4] - d2;
788 bb[7] = bb[5] + d1; // dx, dy
789 bb[10] = bb[8] - d1;
790 bb[11] = bb[9] - d2; // fx, fy
791 bb[14] = bb[12] + d2;
792 bb[15] = bb[13] - d1; // hx, hy
793
794 // adjust all points
795 for ( int i = 0; i < 16; i += 4 )
796 {
797
798 alpha_seg = ( ( i / 4 > 0 ? ( i / 4 ) - 1 : 3 ) ) * M_PI_2 + alpha;
799
800 double best_cp = std::numeric_limits<double>::max();
801
802 for ( std::size_t j = 0; j < convexHull.size(); j++ )
803 {
804 const double cp = GeomFunction::cross_product( bb[i + 2], bb[i + 3], bb[i], bb[i + 1], x[convexHull[j]], y[convexHull[j]] );
805 if ( cp < best_cp )
806 {
807 best_cp = cp;
808 }
809 }
810
811 const double distNearestPoint = best_cp / dref;
812
813 d1 = std::cos( alpha_seg ) * distNearestPoint;
814 d2 = std::sin( alpha_seg ) * distNearestPoint;
815
816 bb[i] += d1; // x
817 bb[i + 1] += d2; // y
818 bb[i + 2] += d1; // x
819 bb[i + 3] += d2; // y
820 }
821
822 // compute and compare AREA
823 const double width = GeomFunction::cross_product( bb[6], bb[7], bb[4], bb[5], bb[12], bb[13] ) / dref;
824 const double length = GeomFunction::cross_product( bb[2], bb[3], bb[0], bb[1], bb[8], bb[9] ) / dref;
825
826 double area = width * length;
827
828 if ( area < 0 )
829 area *= -1;
830
831
832 if ( best_area - area > EPSILON )
833 {
834 best_area = area;
835 best_length = length;
836 best_width = width;
837 best_alpha = alpha;
838 memcpy( best_bb, bb, sizeof( double ) * 16 );
839 }
840 }
841
842 // best bbox is defined
843 for ( int i = 0; i < 16; i = i + 4 )
844 {
845 GeomFunction::computeLineIntersection( best_bb[i], best_bb[i + 1], best_bb[i + 2], best_bb[i + 3],
846 best_bb[( i + 4 ) % 16], best_bb[( i + 5 ) % 16], best_bb[( i + 6 ) % 16], best_bb[( i + 7 ) % 16],
847 &finalBb.x[int ( i / 4 )], &finalBb.y[int ( i / 4 )] );
848 }
849
850 finalBb.alpha = best_alpha;
851 finalBb.width = best_width;
852 finalBb.length = best_length;
853
854 ok = true;
855 return finalBb;
856}
857
858double PointSet::minDistanceToPoint( double px, double py, double *rx, double *ry ) const
859{
860 if ( !mGeos )
862
863 if ( !mGeos )
864 return 0;
865
866 GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
867 try
868 {
869 geos::unique_ptr geosPt( GEOSGeom_createPointFromXY_r( geosctxt, px, py ) );
870 const int type = GEOSGeomTypeId_r( geosctxt, mGeos );
871 const GEOSGeometry *extRing = nullptr;
872 const GEOSPreparedGeometry *preparedExtRing = nullptr;
873
874 if ( type != GEOS_POLYGON )
875 {
876 extRing = mGeos;
877 preparedExtRing = preparedGeom();
878 }
879 else
880 {
881 //for polygons, we want distance to exterior ring (not an interior point)
882 extRing = GEOSGetExteriorRing_r( geosctxt, mGeos );
883 if ( ! mGeosPreparedBoundary )
884 {
885 mGeosPreparedBoundary = GEOSPrepare_r( geosctxt, extRing );
886 }
887 preparedExtRing = mGeosPreparedBoundary;
888 }
889
890 const geos::coord_sequence_unique_ptr nearestCoord( GEOSPreparedNearestPoints_r( geosctxt, preparedExtRing, geosPt.get() ) );
891 double nx;
892 double ny;
893 unsigned int nPoints = 0;
894 GEOSCoordSeq_getSize_r( geosctxt, nearestCoord.get(), &nPoints );
895 if ( nPoints == 0 )
896 return 0;
897
898 ( void )GEOSCoordSeq_getXY_r( geosctxt, nearestCoord.get(), 0, &nx, &ny );
899
900 if ( rx )
901 *rx = nx;
902 if ( ry )
903 *ry = ny;
904
905 return GeomFunction::dist_euc2d_sq( px, py, nx, ny );
906 }
907 catch ( GEOSException &e )
908 {
909 qWarning( "GEOS exception: %s", e.what() );
910 QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
911 return 0;
912 }
913}
914
915void PointSet::getCentroid( double &px, double &py, bool forceInside ) const
916{
917 if ( !mGeos )
919
920 if ( !mGeos )
921 return;
922
923 try
924 {
925 GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
926 geos::unique_ptr centroidGeom( GEOSGetCentroid_r( geosctxt, mGeos ) );
927 if ( centroidGeom )
928 {
929 const GEOSCoordSequence *coordSeq = GEOSGeom_getCoordSeq_r( geosctxt, centroidGeom.get() );
930 unsigned int nPoints = 0;
931 GEOSCoordSeq_getSize_r( geosctxt, coordSeq, &nPoints );
932 if ( nPoints == 0 )
933 return;
934 GEOSCoordSeq_getXY_r( geosctxt, coordSeq, 0, &px, &py );
935 }
936
937 // check if centroid inside in polygon
938 if ( forceInside && !containsPoint( px, py ) )
939 {
940 geos::unique_ptr pointGeom( GEOSPointOnSurface_r( geosctxt, mGeos ) );
941
942 if ( pointGeom )
943 {
944 const GEOSCoordSequence *coordSeq = GEOSGeom_getCoordSeq_r( geosctxt, pointGeom.get() );
945 unsigned int nPoints = 0;
946 GEOSCoordSeq_getSize_r( geosctxt, coordSeq, &nPoints );
947 if ( nPoints == 0 )
948 return;
949
950 GEOSCoordSeq_getXY_r( geosctxt, coordSeq, 0, &px, &py );
951 }
952 }
953 }
954 catch ( GEOSException &e )
955 {
956 qWarning( "GEOS exception: %s", e.what() );
957 QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
958 return;
959 }
960}
961
963{
964 const double x1 = ( xmin > other->xmin ? xmin : other->xmin );
965 const double x2 = ( xmax < other->xmax ? xmax : other->xmax );
966 if ( x1 > x2 )
967 return false;
968 const double y1 = ( ymin > other->ymin ? ymin : other->ymin );
969 const double y2 = ( ymax < other->ymax ? ymax : other->ymax );
970 return y1 <= y2;
971}
972
973void PointSet::getPointByDistance( double *d, double *ad, double dl, double *px, double *py ) const
974{
975 int i;
976 double dx, dy, di;
977 double distr;
978
979 i = 0;
980 if ( dl >= 0 )
981 {
982 while ( i < nbPoints && ad[i] <= dl ) i++;
983 i--;
984 }
985
986 if ( i < nbPoints - 1 )
987 {
988 if ( dl < 0 )
989 {
990 dx = x[nbPoints - 1] - x[0];
991 dy = y[nbPoints - 1] - y[0];
992 di = std::sqrt( dx * dx + dy * dy );
993 }
994 else
995 {
996 dx = x[i + 1] - x[i];
997 dy = y[i + 1] - y[i];
998 di = d[i];
999 }
1000
1001 distr = dl - ad[i];
1002 *px = x[i] + dx * distr / di;
1003 *py = y[i] + dy * distr / di;
1004 }
1005 else // just select last point...
1006 {
1007 *px = x[i];
1008 *py = y[i];
1009 }
1010}
1011
1013{
1014 const GEOSGeometry *thisGeos = geos();
1015 if ( !thisGeos )
1016 return nullptr;
1017
1018 try
1019 {
1020 geos::unique_ptr res( GEOSInterpolate_r( QgsGeos::getGEOSHandler(), thisGeos, distance ) );
1021 return res;
1022 }
1023 catch ( GEOSException &e )
1024 {
1025 qWarning( "GEOS exception: %s", e.what() );
1026 return nullptr;
1027 }
1028}
1029
1030double PointSet::lineLocatePoint( const GEOSGeometry *point ) const
1031{
1032 const GEOSGeometry *thisGeos = geos();
1033 if ( !thisGeos )
1034 return -1;
1035
1036 double distance = -1;
1037 try
1038 {
1039 distance = GEOSProject_r( QgsGeos::getGEOSHandler(), thisGeos, point );
1040 }
1041 catch ( GEOSException &e )
1042 {
1043 qWarning( "GEOS exception: %s", e.what() );
1044 return -1;
1045 }
1046
1047 return distance;
1048}
1049
1050const GEOSGeometry *PointSet::geos() const
1051{
1052 if ( !mGeos )
1054
1055 return mGeos;
1056}
1057
1058double PointSet::length() const
1059{
1060 if ( mLength >= 0 )
1061 return mLength;
1062
1063 if ( !mGeos )
1065
1066 if ( !mGeos )
1067 return -1;
1068
1069 GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
1070
1071 try
1072 {
1073 ( void )GEOSLength_r( geosctxt, mGeos, &mLength );
1074 return mLength;
1075 }
1076 catch ( GEOSException &e )
1077 {
1078 qWarning( "GEOS exception: %s", e.what() );
1079 QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
1080 return -1;
1081 }
1082}
1083
1084double PointSet::area() const
1085{
1086 if ( mArea >= 0 )
1087 return mArea;
1088
1089 if ( !mGeos )
1091
1092 if ( !mGeos )
1093 return -1;
1094
1095 GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
1096
1097 try
1098 {
1099 ( void )GEOSArea_r( geosctxt, mGeos, &mArea );
1100 mArea = std::fabs( mArea );
1101 return mArea;
1102 }
1103 catch ( GEOSException &e )
1104 {
1105 qWarning( "GEOS exception: %s", e.what() );
1106 QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
1107 return -1;
1108 }
1109}
1110
1112{
1113 return qgsDoubleNear( x[0], x[nbPoints - 1] ) && qgsDoubleNear( y[0], y[nbPoints - 1] );
1114}
1115
1116QString PointSet::toWkt() const
1117{
1118 if ( !mGeos )
1120
1121 GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
1122
1123 try
1124 {
1125 GEOSWKTWriter *writer = GEOSWKTWriter_create_r( geosctxt );
1126
1127 char *wkt = GEOSWKTWriter_write_r( geosctxt, writer, mGeos );
1128 const QString res( wkt );
1129
1130 GEOSFree_r( geosctxt, wkt );
1131
1132 GEOSWKTWriter_destroy_r( geosctxt, writer );
1133 writer = nullptr;
1134
1135 return res;
1136 }
1137 catch ( GEOSException &e )
1138 {
1139 qWarning( "GEOS exception: %s", e.what() );
1140 QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
1141 return QString();
1142 }
1143}
1144
1145std::tuple< std::vector< double >, double > PointSet::edgeDistances() const
1146{
1147 std::vector< double > distances( nbPoints );
1148 double totalDistance = 0;
1149 double oldX = -1.0, oldY = -1.0;
1150 for ( int i = 0; i < nbPoints; i++ )
1151 {
1152 if ( i == 0 )
1153 distances[i] = 0;
1154 else
1155 distances[i] = std::sqrt( std::pow( oldX - x[i], 2 ) + std::pow( oldY - y[i], 2 ) );
1156
1157 oldX = x[i];
1158 oldY = y[i];
1159 totalDistance += distances[i];
1160 }
1161 return std::make_tuple( std::move( distances ), totalDistance );
1162}
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:1012
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:1030
OrientedConvexHullBoundingBox computeConvexHullOrientedBoundingBox(bool &ok)
Computes an oriented bounding box for the shape's convex hull.
Definition: pointset.cpp:716
double length() const
Returns length of line geometry.
Definition: pointset.cpp:1058
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:629
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:1084
bool isClosed() const
Returns true if pointset is closed.
Definition: pointset.cpp:1111
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:973
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:962
void getCentroid(double &px, double &py, bool forceInside=false) const
Definition: pointset.cpp:915
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:1116
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:1050
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:858
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:1145
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:2527
Represents the minimum area, oriented bounding box surrounding a convex hull.
Definition: pointset.h:60
#define EPSILON
Definition: util.h:78