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