QGIS API Documentation  3.10.0-A Coruña (6c816b4204)
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 = cHullSize = 0;
44  cHull = nullptr;
45  type = -1;
46 }
47 
48 PointSet::PointSet( int nbPoints, double *x, double *y )
49  : nbPoints( nbPoints )
50  , type( GEOS_POLYGON )
51 {
52  this->x.resize( nbPoints );
53  this->y.resize( nbPoints );
54  int i;
55 
56  for ( i = 0; i < nbPoints; i++ )
57  {
58  this->x[i] = x[i];
59  this->y[i] = y[i];
60  }
61 
62 }
63 
64 PointSet::PointSet( double aX, double aY )
65  : type( GEOS_POINT )
66  , xmin( aX )
67  , xmax( aY )
68  , ymin( aX )
69  , ymax( aY )
70 
71 {
72  nbPoints = cHullSize = 1;
73  x.resize( 1 );
74  y.resize( 1 );
75  x[0] = aX;
76  y[0] = aY;
77 }
78 
80  : xmin( ps.xmin )
81  , xmax( ps.xmax )
82  , ymin( ps.ymin )
83  , ymax( ps.ymax )
84 {
85  int i;
86 
87  nbPoints = ps.nbPoints;
88  x = ps.x;
89  y = ps.y;
90 
91  if ( ps.cHull )
92  {
93  cHullSize = ps.cHullSize;
94  cHull = new int[cHullSize];
95  for ( i = 0; i < cHullSize; i++ )
96  {
97  cHull[i] = ps.cHull[i];
98  }
99  }
100  else
101  {
102  cHull = nullptr;
103  cHullSize = 0;
104  }
105 
106  type = ps.type;
107 
108  holeOf = ps.holeOf;
109 
110  if ( ps.mGeos )
111  {
112  mGeos = GEOSGeom_clone_r( QgsGeos::getGEOSHandler(), ps.mGeos );
113  mOwnsGeom = true;
114  }
115 }
116 
118 {
119  GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
120 
121  bool needClose = false;
122  if ( type == GEOS_POLYGON && ( !qgsDoubleNear( x[0], x[ nbPoints - 1] ) || !qgsDoubleNear( y[0], y[ nbPoints - 1 ] ) ) )
123  {
124  needClose = true;
125  }
126 
127  GEOSCoordSequence *coord = GEOSCoordSeq_create_r( geosctxt, nbPoints + ( needClose ? 1 : 0 ), 2 );
128  for ( int i = 0; i < nbPoints; ++i )
129  {
130  GEOSCoordSeq_setX_r( geosctxt, coord, i, x[i] );
131  GEOSCoordSeq_setY_r( geosctxt, coord, i, y[i] );
132  }
133 
134  //close ring if needed
135  if ( needClose )
136  {
137  GEOSCoordSeq_setX_r( geosctxt, coord, nbPoints, x[0] );
138  GEOSCoordSeq_setY_r( geosctxt, coord, nbPoints, y[0] );
139  }
140 
141  switch ( type )
142  {
143  case GEOS_POLYGON:
144  mGeos = GEOSGeom_createPolygon_r( geosctxt, GEOSGeom_createLinearRing_r( geosctxt, coord ), nullptr, 0 );
145  break;
146 
147  case GEOS_LINESTRING:
148  mGeos = GEOSGeom_createLineString_r( geosctxt, coord );
149  break;
150 
151  case GEOS_POINT:
152  mGeos = GEOSGeom_createPoint_r( geosctxt, coord );
153  break;
154  }
155 
156  mOwnsGeom = true;
157 }
158 
159 const GEOSPreparedGeometry *PointSet::preparedGeom() const
160 {
161  if ( !mGeos )
162  createGeosGeom();
163 
164  if ( !mPreparedGeom )
165  {
166  mPreparedGeom = GEOSPrepare_r( QgsGeos::getGEOSHandler(), mGeos );
167  }
168  return mPreparedGeom;
169 }
170 
172 {
173  GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
174  if ( mOwnsGeom ) // delete old geometry if we own it
175  GEOSGeom_destroy_r( geosctxt, mGeos );
176  GEOSPreparedGeom_destroy_r( geosctxt, mPreparedGeom );
177  mOwnsGeom = false;
178  mGeos = nullptr;
179  mPreparedGeom = nullptr;
180 }
181 
183 {
184  GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
185 
186  if ( mGeos && mOwnsGeom )
187  {
188  GEOSGeom_destroy_r( geosctxt, mGeos );
189  mGeos = nullptr;
190  }
191  GEOSPreparedGeom_destroy_r( geosctxt, mPreparedGeom );
192 
193  deleteCoords();
194 
195  delete[] cHull;
196 }
197 
199 {
200  x.clear();
201  y.clear();
202 }
203 
204 PointSet *PointSet::extractShape( int nbPtSh, int imin, int imax, int fps, int fpe, double fptx, double fpty )
205 {
206 
207  int i, j;
208 
209  PointSet *newShape = new PointSet();
210  newShape->type = GEOS_POLYGON;
211  newShape->nbPoints = nbPtSh;
212  newShape->x.resize( newShape->nbPoints );
213  newShape->y.resize( newShape->nbPoints );
214 
215  // new shape # 1 from imin to imax
216  for ( j = 0, i = imin; i != ( imax + 1 ) % nbPoints; i = ( i + 1 ) % nbPoints, j++ )
217  {
218  newShape->x[j] = x[i];
219  newShape->y[j] = y[i];
220  }
221  // is the cutting point a new one ?
222  if ( fps != fpe )
223  {
224  // yes => so add it
225  newShape->x[j] = fptx;
226  newShape->y[j] = fpty;
227  }
228 
229  return newShape;
230 }
231 
232 std::unique_ptr<PointSet> PointSet::clone() const
233 {
234  return std::unique_ptr< PointSet>( new PointSet( *this ) );
235 }
236 
237 bool PointSet::containsPoint( double x, double y ) const
238 {
239  GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
240  try
241  {
242  GEOSCoordSequence *seq = GEOSCoordSeq_create_r( geosctxt, 1, 2 );
243  GEOSCoordSeq_setX_r( geosctxt, seq, 0, x );
244  GEOSCoordSeq_setY_r( geosctxt, seq, 0, y );
245  geos::unique_ptr point( GEOSGeom_createPoint_r( geosctxt, seq ) );
246  bool result = ( GEOSPreparedContainsProperly_r( geosctxt, preparedGeom(), point.get() ) == 1 );
247 
248  return result;
249  }
250  catch ( GEOSException &e )
251  {
252  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
253  return false;
254  }
255 
256 }
257 
258 bool PointSet::containsLabelCandidate( double x, double y, double width, double height, double alpha ) const
259 {
260  return GeomFunction::containsCandidate( preparedGeom(), x, y, width, height, alpha );
261 }
262 
263 void PointSet::splitPolygons( QLinkedList<PointSet *> &shapes_toProcess,
264  QLinkedList<PointSet *> &shapes_final,
265  double xrm, double yrm )
266 {
267  int i, j;
268 
269  int nbp;
270  std::vector< double > x;
271  std::vector< double > y;
272 
273  int *pts = nullptr;
274 
275  int *cHull = nullptr;
276  int cHullSize;
277 
278  double cp;
279  double bestcp = 0;
280 
281  double bestArea = 0;
282  double area;
283 
284  double base;
285  double b, c;
286  double s;
287 
288  int ihs;
289  int ihn;
290 
291  int ips;
292  int ipn;
293 
294  int holeS = -1; // hole start and end points
295  int holeE = -1;
296 
297  int retainedPt = -1;
298  int pt = 0;
299 
300  double labelArea = xrm * yrm;
301 
302  PointSet *shape = nullptr;
303 
304  while ( !shapes_toProcess.isEmpty() )
305  {
306  shape = shapes_toProcess.takeFirst();
307 
308  x = shape->x;
309  y = shape->y;
310  nbp = shape->nbPoints;
311  pts = new int[nbp];
312 
313  for ( i = 0; i < nbp; i++ )
314  {
315  pts[i] = i;
316  }
317 
318  // conpute convex hull
319  shape->cHullSize = GeomFunction::convexHullId( pts, x, y, nbp, shape->cHull );
320 
321  cHull = shape->cHull;
322  cHullSize = shape->cHullSize;
323 
324  bestArea = 0;
325  retainedPt = -1;
326 
327  // lookup for a hole
328  for ( ihs = 0; ihs < cHullSize; ihs++ )
329  {
330  // ihs->ihn => cHull'seg
331  ihn = ( ihs + 1 ) % cHullSize;
332 
333  ips = cHull[ihs];
334  ipn = ( ips + 1 ) % nbp;
335  if ( ipn != cHull[ihn] ) // next point on shape is not the next point on cHull => there is a hole here !
336  {
337  bestcp = 0;
338  pt = -1;
339  // lookup for the deepest point in the hole
340  for ( i = ips; i != cHull[ihn]; i = ( i + 1 ) % nbp )
341  {
342  cp = std::fabs( GeomFunction::cross_product( x[cHull[ihs]], y[cHull[ihs]],
343  x[cHull[ihn]], y[cHull[ihn]],
344  x[i], y[i] ) );
345  if ( cp - bestcp > EPSILON )
346  {
347  bestcp = cp;
348  pt = i;
349  }
350  }
351 
352  if ( pt != -1 )
353  {
354  // compute the ihs->ihn->pt triangle's area
355  base = GeomFunction::dist_euc2d( x[cHull[ihs]], y[cHull[ihs]],
356  x[cHull[ihn]], y[cHull[ihn]] );
357 
358  b = GeomFunction::dist_euc2d( x[cHull[ihs]], y[cHull[ihs]],
359  x[pt], y[pt] );
360 
361  c = GeomFunction::dist_euc2d( x[cHull[ihn]], y[cHull[ihn]],
362  x[pt], y[pt] );
363 
364  s = ( base + b + c ) / 2; // s = half perimeter
365  area = s * ( s - base ) * ( s - b ) * ( s - c );
366  if ( area < 0 )
367  area = -area;
368 
369  // retain the biggest area
370  if ( area - bestArea > EPSILON )
371  {
372  bestArea = area;
373  retainedPt = pt;
374  holeS = ihs;
375  holeE = ihn;
376  }
377  }
378  }
379  }
380 
381  // we have a hole, its area, and the deppest point in hole
382  // we're going to find the second point to cup the shape
383  // holeS = hole starting point
384  // holeE = hole ending point
385  // retainedPt = deppest point in hole
386  // bestArea = area of triangle HoleS->holeE->retainedPoint
387  bestArea = std::sqrt( bestArea );
388  double cx, cy, dx, dy, ex, ey, fx, fy, seg_length, ptx = 0, pty = 0, fptx = 0, fpty = 0;
389  int ps = -1, pe = -1, fps = -1, fpe = -1;
390  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)
391  {
392  c = std::numeric_limits<double>::max();
393 
394  // iterate on all shape points except points which are in the hole
395  bool isValid;
396  int k, l;
397  for ( i = ( cHull[holeE] + 1 ) % nbp; i != ( cHull[holeS] - 1 + nbp ) % nbp; i = j )
398  {
399  j = ( i + 1 ) % nbp; // i->j is shape segment not in hole
400 
401  // compute distance between retainedPoint and segment
402  // whether perpendicular distance (if retaindPoint is fronting segment i->j)
403  // or distance between retainedPt and i or j (choose the nearest)
404  seg_length = GeomFunction::dist_euc2d( x[i], y[i], x[j], y[j] );
405  cx = ( x[i] + x[j] ) / 2.0;
406  cy = ( y[i] + y[j] ) / 2.0;
407  dx = cy - y[i];
408  dy = cx - x[i];
409 
410  ex = cx - dx;
411  ey = cy + dy;
412  fx = cx + dx;
413  fy = cy - dy;
414 
415  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
416  {
417  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] ) ) )
418  {
419  b = ex;
420  ps = i;
421  pe = i;
422  }
423  else
424  {
425  b = ey;
426  ps = j;
427  pe = j;
428  }
429  }
430  else // point fronting i->j => compute pependicular distance => create a new point
431  {
432  b = GeomFunction::cross_product( x[i], y[i], x[j], y[j], x[retainedPt], y[retainedPt] ) / seg_length;
433  b *= b;
434  ps = i;
435  pe = j;
436 
437  if ( !GeomFunction::computeLineIntersection( x[i], y[i], x[j], y[j], x[retainedPt], y[retainedPt], x[retainedPt] - dx, y[retainedPt] + dy, &ptx, &pty ) )
438  {
439  //error - it should intersect the line
440  }
441  }
442 
443  isValid = true;
444  double pointX, pointY;
445  if ( ps == pe )
446  {
447  pointX = x[pe];
448  pointY = y[pe];
449  }
450  else
451  {
452  pointX = ptx;
453  pointY = pty;
454  }
455 
456  for ( k = cHull[holeS]; k != cHull[holeE]; k = ( k + 1 ) % nbp )
457  {
458  l = ( k + 1 ) % nbp;
459  if ( GeomFunction::isSegIntersects( x[retainedPt], y[retainedPt], pointX, pointY, x[k], y[k], x[l], y[l] ) )
460  {
461  isValid = false;
462  break;
463  }
464  }
465 
466 
467  if ( isValid && b < c )
468  {
469  c = b;
470  fps = ps;
471  fpe = pe;
472  fptx = ptx;
473  fpty = pty;
474  }
475  } // for point which are not in hole
476 
477  // we will cut the shapeu in two new shapes, one from [retainedPoint] to [newPoint] and one form [newPoint] to [retainedPoint]
478  int imin = retainedPt;
479  int imax = ( ( ( fps < retainedPt && fpe < retainedPt ) || ( fps > retainedPt && fpe > retainedPt ) ) ? std::min( fps, fpe ) : std::max( fps, fpe ) );
480 
481  int nbPtSh1, nbPtSh2; // how many points in new shapes ?
482  if ( imax > imin )
483  nbPtSh1 = imax - imin + 1 + ( fpe != fps );
484  else
485  nbPtSh1 = imax + nbp - imin + 1 + ( fpe != fps );
486 
487  if ( ( imax == fps ? fpe : fps ) < imin )
488  nbPtSh2 = imin - ( imax == fps ? fpe : fps ) + 1 + ( fpe != fps );
489  else
490  nbPtSh2 = imin + nbp - ( imax == fps ? fpe : fps ) + 1 + ( fpe != fps );
491 
492  if ( retainedPt == -1 || fps == -1 || fpe == -1 )
493  {
494  if ( shape->parent )
495  delete shape;
496  }
497  // check for useless splitting
498  else if ( imax == imin || nbPtSh1 <= 2 || nbPtSh2 <= 2 || nbPtSh1 == nbp || nbPtSh2 == nbp )
499  {
500  shapes_final.append( shape );
501  }
502  else
503  {
504 
505  PointSet *newShape = shape->extractShape( nbPtSh1, imin, imax, fps, fpe, fptx, fpty );
506 
507  if ( shape->parent )
508  newShape->parent = shape->parent;
509  else
510  newShape->parent = shape;
511 
512  shapes_toProcess.append( newShape );
513 
514  if ( imax == fps )
515  imax = fpe;
516  else
517  imax = fps;
518 
519  newShape = shape->extractShape( nbPtSh2, imax, imin, fps, fpe, fptx, fpty );
520 
521  if ( shape->parent )
522  newShape->parent = shape->parent;
523  else
524  newShape->parent = shape;
525 
526  shapes_toProcess.append( newShape );
527 
528  if ( shape->parent )
529  delete shape;
530  }
531  }
532  else
533  {
534  shapes_final.append( shape );
535  }
536  delete[] pts;
537  }
538 }
539 
540 void PointSet::extendLineByDistance( double startDistance, double endDistance, double smoothDistance )
541 {
542  if ( nbPoints < 2 )
543  return;
544 
545  double x0 = x[0];
546  double y0 = y[0];
547  if ( startDistance > 0 )
548  {
549  // trace forward by smoothDistance
550  double x1 = x[1];
551  double y1 = y[1];
552 
553  double distanceConsumed = 0;
554  double lastX = x0;
555  double lastY = y0;
556  for ( int i = 1; i < nbPoints; ++i )
557  {
558  double thisX = x[i];
559  double thisY = y[i];
560  const double thisSegmentLength = std::sqrt( ( thisX - lastX ) * ( thisX - lastX ) + ( thisY - lastY ) * ( thisY - lastY ) );
561  distanceConsumed += thisSegmentLength;
562  if ( distanceConsumed >= smoothDistance )
563  {
564  double c = ( distanceConsumed - smoothDistance ) / thisSegmentLength;
565  x1 = lastX + c * ( thisX - lastX );
566  y1 = lastY + c * ( thisY - lastY );
567  break;
568  }
569  lastX = thisX;
570  lastY = thisY;
571  }
572 
573  const double distance = std::sqrt( ( x1 - x0 ) * ( x1 - x0 ) + ( y1 - y0 ) * ( y1 - y0 ) );
574  const double extensionFactor = ( startDistance + distance ) / distance;
575  const QgsPointXY newStart = QgsGeometryUtils::interpolatePointOnLine( x1, y1, x0, y0, extensionFactor );
576  x0 = newStart.x();
577  y0 = newStart.y();
578  // defer actually changing the stored start until we've calculated the new end point
579  }
580 
581  if ( endDistance > 0 )
582  {
583  double xend0 = x[nbPoints - 1];
584  double yend0 = y[nbPoints - 1];
585  double xend1 = x[nbPoints - 2];
586  double yend1 = y[nbPoints - 2];
587 
588  // trace backward by smoothDistance
589  double distanceConsumed = 0;
590  double lastX = x0;
591  double lastY = y0;
592  for ( int i = nbPoints - 2; i >= 0; --i )
593  {
594  double thisX = x[i];
595  double thisY = y[i];
596  const double thisSegmentLength = std::sqrt( ( thisX - lastX ) * ( thisX - lastX ) + ( thisY - lastY ) * ( thisY - lastY ) );
597  distanceConsumed += thisSegmentLength;
598  if ( distanceConsumed >= smoothDistance )
599  {
600  double c = ( distanceConsumed - smoothDistance ) / thisSegmentLength;
601  xend1 = lastX + c * ( thisX - lastX );
602  yend1 = lastY + c * ( thisY - lastY );
603  break;
604  }
605  lastX = thisX;
606  lastY = thisY;
607  }
608 
609  const double distance = std::sqrt( ( xend1 - xend0 ) * ( xend1 - xend0 ) + ( yend1 - yend0 ) * ( yend1 - yend0 ) );
610  const double extensionFactor = ( endDistance + distance ) / distance;
611  const QgsPointXY newEnd = QgsGeometryUtils::interpolatePointOnLine( xend1, yend1, xend0, yend0, extensionFactor );
612  x.emplace_back( newEnd.x() );
613  y.emplace_back( newEnd.y() );
614  nbPoints++;
615  }
616 
617  if ( startDistance > 0 )
618  {
619  x.insert( x.begin(), x0 );
620  y.insert( y.begin(), y0 );
621  nbPoints++;
622  }
623 
624  invalidateGeos();
625 }
626 
628 {
629  int i;
630  int j;
631 
632  double bbox[4]; // xmin, ymin, xmax, ymax
633 
634  double alpha;
635  int alpha_d;
636 
637  double alpha_seg;
638 
639  double dref;
640  double d1, d2;
641 
642  double bb[16]; // {ax, ay, bx, by, cx, cy, dx, dy, ex, ey, fx, fy, gx, gy, hx, hy}}
643 
644  double cp;
645  double best_cp;
646  double distNearestPoint;
647 
648  double area;
649  double width;
650  double length;
651 
652  double best_area = std::numeric_limits<double>::max();
653  double best_alpha = -1;
654  double best_bb[16];
655  double best_length = 0;
656  double best_width = 0;
657 
658 
659  bbox[0] = std::numeric_limits<double>::max();
660  bbox[1] = std::numeric_limits<double>::max();
661  bbox[2] = std::numeric_limits<double>::lowest();
662  bbox[3] = std::numeric_limits<double>::lowest();
663 
664  for ( i = 0; i < cHullSize; i++ )
665  {
666  if ( x[cHull[i]] < bbox[0] )
667  bbox[0] = x[cHull[i]];
668 
669  if ( x[cHull[i]] > bbox[2] )
670  bbox[2] = x[cHull[i]];
671 
672  if ( y[cHull[i]] < bbox[1] )
673  bbox[1] = y[cHull[i]];
674 
675  if ( y[cHull[i]] > bbox[3] )
676  bbox[3] = y[cHull[i]];
677  }
678 
679 
680  dref = bbox[2] - bbox[0];
681 
682  for ( alpha_d = 0; alpha_d < 90; alpha_d++ )
683  {
684  alpha = alpha_d * M_PI / 180.0;
685  d1 = std::cos( alpha ) * dref;
686  d2 = std::sin( alpha ) * dref;
687 
688  bb[0] = bbox[0];
689  bb[1] = bbox[3]; // ax, ay
690 
691  bb[4] = bbox[0];
692  bb[5] = bbox[1]; // cx, cy
693 
694  bb[8] = bbox[2];
695  bb[9] = bbox[1]; // ex, ey
696 
697  bb[12] = bbox[2];
698  bb[13] = bbox[3]; // gx, gy
699 
700 
701  bb[2] = bb[0] + d1;
702  bb[3] = bb[1] + d2; // bx, by
703  bb[6] = bb[4] - d2;
704  bb[7] = bb[5] + d1; // dx, dy
705  bb[10] = bb[8] - d1;
706  bb[11] = bb[9] - d2; // fx, fy
707  bb[14] = bb[12] + d2;
708  bb[15] = bb[13] - d1; // hx, hy
709 
710  // adjust all points
711  for ( i = 0; i < 16; i += 4 )
712  {
713 
714  alpha_seg = ( ( i / 4 > 0 ? ( i / 4 ) - 1 : 3 ) ) * M_PI_2 + alpha;
715 
716  best_cp = std::numeric_limits<double>::max();
717  for ( j = 0; j < nbPoints; j++ )
718  {
719  cp = GeomFunction::cross_product( bb[i + 2], bb[i + 3], bb[i], bb[i + 1], x[cHull[j]], y[cHull[j]] );
720  if ( cp < best_cp )
721  {
722  best_cp = cp;
723  }
724  }
725 
726  distNearestPoint = best_cp / dref;
727 
728  d1 = std::cos( alpha_seg ) * distNearestPoint;
729  d2 = std::sin( alpha_seg ) * distNearestPoint;
730 
731  bb[i] += d1; // x
732  bb[i + 1] += d2; // y
733  bb[i + 2] += d1; // x
734  bb[i + 3] += d2; // y
735  }
736 
737  // compute and compare AREA
738  width = GeomFunction::cross_product( bb[6], bb[7], bb[4], bb[5], bb[12], bb[13] ) / dref;
739  length = GeomFunction::cross_product( bb[2], bb[3], bb[0], bb[1], bb[8], bb[9] ) / dref;
740 
741  area = width * length;
742 
743  if ( area < 0 )
744  area *= -1;
745 
746 
747  if ( best_area - area > EPSILON )
748  {
749  best_area = area;
750  best_length = length;
751  best_width = width;
752  best_alpha = alpha;
753  memcpy( best_bb, bb, sizeof( double ) * 16 );
754  }
755  }
756 
757  // best bbox is defined
758 
759  CHullBox *finalBb = new CHullBox();
760 
761  for ( i = 0; i < 16; i = i + 4 )
762  {
763  GeomFunction::computeLineIntersection( best_bb[i], best_bb[i + 1], best_bb[i + 2], best_bb[i + 3],
764  best_bb[( i + 4 ) % 16], best_bb[( i + 5 ) % 16], best_bb[( i + 6 ) % 16], best_bb[( i + 7 ) % 16],
765  &finalBb->x[int ( i / 4 )], &finalBb->y[int ( i / 4 )] );
766  }
767 
768  finalBb->alpha = best_alpha;
769  finalBb->width = best_width;
770  finalBb->length = best_length;
771 
772  return finalBb;
773 }
774 
775 double PointSet::minDistanceToPoint( double px, double py, double *rx, double *ry ) const
776 {
777  if ( !mGeos )
778  createGeosGeom();
779 
780  if ( !mGeos )
781  return 0;
782 
783  GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
784  try
785  {
786  GEOSCoordSequence *coord = GEOSCoordSeq_create_r( geosctxt, 1, 2 );
787  GEOSCoordSeq_setX_r( geosctxt, coord, 0, px );
788  GEOSCoordSeq_setY_r( geosctxt, coord, 0, py );
789  geos::unique_ptr geosPt( GEOSGeom_createPoint_r( geosctxt, coord ) );
790 
791  int type = GEOSGeomTypeId_r( geosctxt, mGeos );
792  const GEOSGeometry *extRing = nullptr;
793  if ( type != GEOS_POLYGON )
794  {
795  extRing = mGeos;
796  }
797  else
798  {
799  //for polygons, we want distance to exterior ring (not an interior point)
800  extRing = GEOSGetExteriorRing_r( geosctxt, mGeos );
801  }
802  geos::coord_sequence_unique_ptr nearestCoord( GEOSNearestPoints_r( geosctxt, extRing, geosPt.get() ) );
803  double nx;
804  double ny;
805  ( void )GEOSCoordSeq_getX_r( geosctxt, nearestCoord.get(), 0, &nx );
806  ( void )GEOSCoordSeq_getY_r( geosctxt, nearestCoord.get(), 0, &ny );
807 
808  if ( rx )
809  *rx = nx;
810  if ( ry )
811  *ry = ny;
812 
813  return GeomFunction::dist_euc2d_sq( px, py, nx, ny );
814  }
815  catch ( GEOSException &e )
816  {
817  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
818  return 0;
819  }
820 }
821 
822 void PointSet::getCentroid( double &px, double &py, bool forceInside ) const
823 {
824  if ( !mGeos )
825  createGeosGeom();
826 
827  if ( !mGeos )
828  return;
829 
830  try
831  {
832  GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
833  geos::unique_ptr centroidGeom( GEOSGetCentroid_r( geosctxt, mGeos ) );
834  if ( centroidGeom )
835  {
836  const GEOSCoordSequence *coordSeq = GEOSGeom_getCoordSeq_r( geosctxt, centroidGeom.get() );
837  GEOSCoordSeq_getX_r( geosctxt, coordSeq, 0, &px );
838  GEOSCoordSeq_getY_r( geosctxt, coordSeq, 0, &py );
839  }
840 
841  // check if centroid inside in polygon
842  if ( forceInside && !containsPoint( px, py ) )
843  {
844  geos::unique_ptr pointGeom( GEOSPointOnSurface_r( geosctxt, mGeos ) );
845 
846  if ( pointGeom )
847  {
848  const GEOSCoordSequence *coordSeq = GEOSGeom_getCoordSeq_r( geosctxt, pointGeom.get() );
849  GEOSCoordSeq_getX_r( geosctxt, coordSeq, 0, &px );
850  GEOSCoordSeq_getY_r( geosctxt, coordSeq, 0, &py );
851  }
852  }
853  }
854  catch ( GEOSException &e )
855  {
856  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
857  return;
858  }
859 }
860 
861 bool PointSet::boundingBoxIntersects( const PointSet *other ) const
862 {
863  double x1 = ( xmin > other->xmin ? xmin : other->xmin );
864  double x2 = ( xmax < other->xmax ? xmax : other->xmax );
865  if ( x1 > x2 )
866  return false;
867  double y1 = ( ymin > other->ymin ? ymin : other->ymin );
868  double y2 = ( ymax < other->ymax ? ymax : other->ymax );
869  return y1 <= y2;
870 }
871 
872 void PointSet::getPointByDistance( double *d, double *ad, double dl, double *px, double *py )
873 {
874  int i;
875  double dx, dy, di;
876  double distr;
877 
878  i = 0;
879  if ( dl >= 0 )
880  {
881  while ( i < nbPoints && ad[i] <= dl ) i++;
882  i--;
883  }
884 
885  if ( i < nbPoints - 1 )
886  {
887  if ( dl < 0 )
888  {
889  dx = x[nbPoints - 1] - x[0];
890  dy = y[nbPoints - 1] - y[0];
891  di = std::sqrt( dx * dx + dy * dy );
892  }
893  else
894  {
895  dx = x[i + 1] - x[i];
896  dy = y[i + 1] - y[i];
897  di = d[i];
898  }
899 
900  distr = dl - ad[i];
901  *px = x[i] + dx * distr / di;
902  *py = y[i] + dy * distr / di;
903  }
904  else // just select last point...
905  {
906  *px = x[i];
907  *py = y[i];
908  }
909 }
910 
911 const GEOSGeometry *PointSet::geos() const
912 {
913  if ( !mGeos )
914  createGeosGeom();
915 
916  return mGeos;
917 }
918 
919 double PointSet::length() const
920 {
921  if ( !mGeos )
922  createGeosGeom();
923 
924  if ( !mGeos )
925  return -1;
926 
927  GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
928 
929  try
930  {
931  double len = 0;
932  ( void )GEOSLength_r( geosctxt, mGeos, &len );
933  return len;
934  }
935  catch ( GEOSException &e )
936  {
937  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
938  return -1;
939  }
940 }
941 
942 bool PointSet::isClosed() const
943 {
944  return qgsDoubleNear( x[0], x[nbPoints - 1] ) && qgsDoubleNear( y[0], y[nbPoints - 1] );
945 }
double length
Definition: pointset.h:62
void invalidateGeos()
Definition: pointset.cpp:171
bool containsPoint(double x, double y) const
Tests whether point set contains a specified point.
Definition: pointset.cpp:237
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.
PointSet * extractShape(int nbPtSh, int imin, int imax, int fps, int fpe, double fptx, double fpty)
Definition: pointset.cpp:204
double y
Definition: qgspointxy.h:48
A class to represent a 2D point.
Definition: qgspointxy.h:43
bool qgsDoubleNear(double a, double b, double epsilon=4 *std::numeric_limits< double >::epsilon())
Compare two doubles (but allow some difference)
Definition: qgis.h:280
std::vector< double > x
Definition: pointset.h:188
void createGeosGeom() const
Definition: pointset.cpp:117
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.
bool isClosed() const
Returns true if pointset is closed.
Definition: pointset.cpp:942
struct pal::_cHullBox CHullBox
static int convexHullId(int *id, const std::vector< double > &x, const std::vector< double > &y, int n, int *&cHull)
Compute the convex hull in O(n·log(n))
static GEOSContextHandle_t getGEOSHandler()
Definition: qgsgeos.cpp:2825
static double cross_product(double x1, double y1, double x2, double y2, double x3, double y3)
Definition: geomfunction.h:61
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:861
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
double width
Definition: pointset.h:61
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.
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:872
PointSet * parent
Definition: pointset.h:201
virtual ~PointSet()
Definition: pointset.cpp:182
double ymax
Definition: pointset.h:215
double xmin
Definition: pointset.h:212
PointSet * holeOf
Definition: pointset.h:200
std::unique_ptr< GEOSGeometry, GeosDeleter > unique_ptr
Scoped GEOS pointer.
Definition: qgsgeos.h:79
static void logMessage(const QString &message, const QString &tag=QString(), Qgis::MessageLevel level=Qgis::Warning, bool notifyUser=true)
Adds a message to the log instance (and creates it if necessary).
double ymin
Definition: pointset.h:214
int * cHull
Definition: pointset.h:195
static double dist_euc2d(double x1, double y1, double x2, double y2)
Definition: geomfunction.h:66
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, y1) to (x2, y2).
void deleteCoords()
Definition: pointset.cpp:198
double x
Definition: qgspointxy.h:47
CHullBox * compute_chull_bbox()
Definition: pointset.cpp:627
double x[4]
Definition: pointset.h:56
std::unique_ptr< PointSet > clone() const
Returns a copy of the point set.
Definition: pointset.cpp:232
double length() const
Returns length of line geometry.
Definition: pointset.cpp:919
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:540
double y[4]
Definition: pointset.h:57
#define EPSILON
Definition: util.h:75
GEOSGeometry * mGeos
Definition: pointset.h:192
const GEOSPreparedGeometry * preparedGeom() const
Definition: pointset.cpp:159
double alpha
Definition: pointset.h:59
std::unique_ptr< GEOSCoordSequence, GeosDeleter > coord_sequence_unique_ptr
Scoped GEOS coordinate sequence pointer.
Definition: qgsgeos.h:94
static double dist_euc2d_sq(double x1, double y1, double x2, double y2)
Definition: geomfunction.h:71
void getCentroid(double &px, double &py, bool forceInside=false) const
Definition: pointset.cpp:822
static void splitPolygons(QLinkedList< PointSet *> &shapes_toProcess, QLinkedList< PointSet *> &shapes_final, double xrm, double yrm)
Split a concave shape into several convex shapes.
Definition: pointset.cpp:263
const GEOSGeometry * geos() const
Returns the point set&#39;s GEOS geometry.
Definition: pointset.cpp:911
double xmax
Definition: pointset.h:213
std::vector< double > y
Definition: pointset.h:189
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:775
bool mOwnsGeom
Definition: pointset.h:193
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:258