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