QGIS API Documentation  3.23.0-Master (eb871beae0)
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 
39 using namespace pal;
40 
42 {
43  nbPoints = 0;
44  type = -1;
45 }
46 
47 PointSet::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 
63 PointSet::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 
150 const GEOSPreparedGeometry *PointSet::preparedGeom() const
151 {
152  if ( !mGeos )
153  createGeosGeom();
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 
234 std::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 
261 std::unique_ptr<PointSet> PointSet::clone() const
262 {
263  return std::unique_ptr< PointSet>( new PointSet( *this ) );
264 }
265 
266 bool 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 
292 bool 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 
297 QLinkedList<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 
548 void PointSet::offsetCurveByDistance( double distance )
549 {
550  if ( !mGeos )
551  createGeosGeom();
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 
619  invalidateGeos();
620  mGeos = newGeos;
621  mOwnsGeom = true;
622 }
623 
624 void 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 
708  invalidateGeos();
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 
853 double PointSet::minDistanceToPoint( double px, double py, double *rx, double *ry ) const
854 {
855  if ( !mGeos )
856  createGeosGeom();
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 
932 void PointSet::getCentroid( double &px, double &py, bool forceInside ) const
933 {
934  if ( !mGeos )
935  createGeosGeom();
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 
989 bool PointSet::boundingBoxIntersects( const PointSet *other ) const
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 
1000 void PointSet::getPointByDistance( double *d, double *ad, double dl, double *px, double *py )
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 
1039 const GEOSGeometry *PointSet::geos() const
1040 {
1041  if ( !mGeos )
1042  createGeosGeom();
1043 
1044  return mGeos;
1045 }
1046 
1047 double PointSet::length() const
1048 {
1049  if ( mLength >= 0 )
1050  return mLength;
1051 
1052  if ( !mGeos )
1053  createGeosGeom();
1054 
1055  if ( !mGeos )
1056  return -1;
1057 
1058  GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
1059 
1060  try
1061  {
1062  ( void )GEOSLength_r( geosctxt, mGeos, &mLength );
1063  return mLength;
1064  }
1065  catch ( GEOSException &e )
1066  {
1067  qWarning( "GEOS exception: %s", e.what() );
1068  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
1069  return -1;
1070  }
1071 }
1072 
1073 double PointSet::area() const
1074 {
1075  if ( mArea >= 0 )
1076  return mArea;
1077 
1078  if ( !mGeos )
1079  createGeosGeom();
1080 
1081  if ( !mGeos )
1082  return -1;
1083 
1084  GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
1085 
1086  try
1087  {
1088  ( void )GEOSArea_r( geosctxt, mGeos, &mArea );
1089  mArea = std::fabs( mArea );
1090  return mArea;
1091  }
1092  catch ( GEOSException &e )
1093  {
1094  qWarning( "GEOS exception: %s", e.what() );
1095  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
1096  return -1;
1097  }
1098 }
1099 
1101 {
1102  return qgsDoubleNear( x[0], x[nbPoints - 1] ) && qgsDoubleNear( y[0], y[nbPoints - 1] );
1103 }
1104 
1105 QString PointSet::toWkt() const
1106 {
1107  if ( !mGeos )
1108  createGeosGeom();
1109 
1110  GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
1111 
1112  try
1113  {
1114  GEOSWKTWriter *writer = GEOSWKTWriter_create_r( geosctxt );
1115 
1116  char *wkt = GEOSWKTWriter_write_r( geosctxt, writer, mGeos );
1117  const QString res( wkt );
1118 
1119  GEOSFree_r( geosctxt, wkt );
1120 
1121  GEOSWKTWriter_destroy_r( geosctxt, writer );
1122  writer = nullptr;
1123 
1124  return res;
1125  }
1126  catch ( GEOSException &e )
1127  {
1128  qWarning( "GEOS exception: %s", e.what() );
1129  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
1130  return QString();
1131  }
1132 }
1133 
1134 std::tuple< std::vector< double >, double > PointSet::edgeDistances() const
1135 {
1136  std::vector< double > distances( nbPoints );
1137  double totalDistance = 0;
1138  double oldX = -1.0, oldY = -1.0;
1139  for ( int i = 0; i < nbPoints; i++ )
1140  {
1141  if ( i == 0 )
1142  distances[i] = 0;
1143  else
1144  distances[i] = std::sqrt( std::pow( oldX - x[i], 2 ) + std::pow( oldY - y[i], 2 ) );
1145 
1146  oldX = x[i];
1147  oldY = y[i];
1148  totalDistance += distances[i];
1149  }
1150  return std::make_tuple( std::move( distances ), totalDistance );
1151 }
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:3222
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:76
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
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:1047
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:250
double ymin
Definition: pointset.h:249
double area() const
Returns area of polygon geometry.
Definition: pointset.cpp:1073
bool isClosed() const
Returns true if pointset is closed.
Definition: pointset.cpp:1100
PointSet * holeOf
Definition: pointset.h:230
void createGeosGeom() const
Definition: pointset.cpp:100
std::vector< double > y
Definition: pointset.h:220
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:219
void offsetCurveByDistance(double distance)
Offsets linestrings by the specified distance.
Definition: pointset.cpp:548
std::vector< int > convexHull
Definition: pointset.h:226
const GEOSPreparedGeometry * preparedGeom() const
Definition: pointset.cpp:150
QString toWkt() const
Returns a WKT representation of the point set.
Definition: pointset.cpp:1105
GEOSGeometry * mGeos
Definition: pointset.h:223
double xmin
Definition: pointset.h:247
const GEOSGeometry * geos() const
Returns the point set's GEOS geometry.
Definition: pointset.cpp:1039
void invalidateGeos()
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:248
void getPointByDistance(double *d, double *ad, double dl, double *px, double *py)
Gets a point a set distance along a line geometry.
Definition: pointset.cpp:1000
std::unique_ptr< PointSet > extractShape(int nbPtSh, int imin, int imax, int fps, int fpe, double fptx, double fpty)
Does...
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:1134
PointSet * parent
Definition: pointset.h:231
double mLength
Definition: pointset.h:234
double mArea
Definition: pointset.h:233
bool mOwnsGeom
Definition: pointset.h:224
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:1504
Represents the minimum area, oriented bounding box surrounding a convex hull.
Definition: pointset.h:59
#define EPSILON
Definition: util.h:78