QGIS API Documentation  3.0.2-Girona (307d082)
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 <qglobal.h>
37 
38 using namespace pal;
39 
41 {
42  nbPoints = cHullSize = 0;
43  x = nullptr;
44  y = nullptr;
45  cHull = nullptr;
46  type = -1;
47 }
48 
49 PointSet::PointSet( int nbPoints, double *x, double *y )
50  : cHullSize( 0 )
51 {
52  this->nbPoints = nbPoints;
53  this->x = new double[nbPoints];
54  this->y = new double[nbPoints];
55  int i;
56 
57  for ( i = 0; i < nbPoints; i++ )
58  {
59  this->x[i] = x[i];
60  this->y[i] = y[i];
61  }
62  type = GEOS_POLYGON;
63  cHull = nullptr;
64 }
65 
66 PointSet::PointSet( double aX, double aY )
67  : xmin( aX )
68  , xmax( aY )
69  , ymin( aX )
70  , ymax( aY )
71 {
72  nbPoints = cHullSize = 1;
73  x = new double[1];
74  y = new double[1];
75  x[0] = aX;
76  y[0] = aY;
77 
78  cHull = nullptr;
79  parent = nullptr;
80  holeOf = nullptr;
81 
82  type = GEOS_POINT;
83 }
84 
86  : xmin( ps.xmin )
87  , xmax( ps.xmax )
88  , ymin( ps.ymin )
89  , ymax( ps.ymax )
90 {
91  int i;
92 
93  nbPoints = ps.nbPoints;
94  x = new double[nbPoints];
95  y = new double[nbPoints];
96  memcpy( x, ps.x, sizeof( double )* nbPoints );
97  memcpy( y, ps.y, sizeof( double )* nbPoints );
98 
99  if ( ps.cHull )
100  {
101  cHullSize = ps.cHullSize;
102  cHull = new int[cHullSize];
103  for ( i = 0; i < cHullSize; i++ )
104  {
105  cHull[i] = ps.cHull[i];
106  }
107  }
108  else
109  {
110  cHull = nullptr;
111  cHullSize = 0;
112  }
113 
114  type = ps.type;
115 
116  holeOf = ps.holeOf;
117 
118  if ( ps.mGeos )
119  {
120  mGeos = GEOSGeom_clone_r( geosContext(), ps.mGeos );
121  mOwnsGeom = true;
122  }
123 }
124 
126 {
127  GEOSContextHandle_t geosctxt = geosContext();
128 
129  bool needClose = false;
130  if ( type == GEOS_POLYGON && ( !qgsDoubleNear( x[0], x[ nbPoints - 1] ) || !qgsDoubleNear( y[0], y[ nbPoints - 1 ] ) ) )
131  {
132  needClose = true;
133  }
134 
135  GEOSCoordSequence *coord = GEOSCoordSeq_create_r( geosctxt, nbPoints + ( needClose ? 1 : 0 ), 2 );
136  for ( int i = 0; i < nbPoints; ++i )
137  {
138  GEOSCoordSeq_setX_r( geosctxt, coord, i, x[i] );
139  GEOSCoordSeq_setY_r( geosctxt, coord, i, y[i] );
140  }
141 
142  //close ring if needed
143  if ( needClose )
144  {
145  GEOSCoordSeq_setX_r( geosctxt, coord, nbPoints, x[0] );
146  GEOSCoordSeq_setY_r( geosctxt, coord, nbPoints, y[0] );
147  }
148 
149  switch ( type )
150  {
151  case GEOS_POLYGON:
152  mGeos = GEOSGeom_createPolygon_r( geosctxt, GEOSGeom_createLinearRing_r( geosctxt, coord ), nullptr, 0 );
153  break;
154 
155  case GEOS_LINESTRING:
156  mGeos = GEOSGeom_createLineString_r( geosctxt, coord );
157  break;
158 
159  case GEOS_POINT:
160  mGeos = GEOSGeom_createPoint_r( geosctxt, coord );
161  break;
162  }
163 
164  mOwnsGeom = true;
165 }
166 
167 const GEOSPreparedGeometry *PointSet::preparedGeom() const
168 {
169  if ( !mGeos )
170  createGeosGeom();
171 
172  if ( !mPreparedGeom )
173  {
174  mPreparedGeom = GEOSPrepare_r( geosContext(), mGeos );
175  }
176  return mPreparedGeom;
177 }
178 
180 {
181  GEOSContextHandle_t geosctxt = geosContext();
182  if ( mOwnsGeom ) // delete old geometry if we own it
183  GEOSGeom_destroy_r( geosctxt, mGeos );
184  GEOSPreparedGeom_destroy_r( geosctxt, mPreparedGeom );
185  mOwnsGeom = false;
186  mGeos = nullptr;
187  mPreparedGeom = nullptr;
188 }
189 
191 {
192  GEOSContextHandle_t geosctxt = geosContext();
193 
194  if ( mGeos && mOwnsGeom )
195  {
196  GEOSGeom_destroy_r( geosctxt, mGeos );
197  mGeos = nullptr;
198  }
199  GEOSPreparedGeom_destroy_r( geosctxt, mPreparedGeom );
200 
201  deleteCoords();
202 
203  if ( cHull )
204  delete[] cHull;
205 }
206 
208 {
209  delete[] x;
210  delete[] y;
211  x = nullptr;
212  y = nullptr;
213 }
214 
215 PointSet *PointSet::extractShape( int nbPtSh, int imin, int imax, int fps, int fpe, double fptx, double fpty )
216 {
217 
218  int i, j;
219 
220  PointSet *newShape = new PointSet();
221  newShape->type = GEOS_POLYGON;
222  newShape->nbPoints = nbPtSh;
223  newShape->x = new double[newShape->nbPoints];
224  newShape->y = new double[newShape->nbPoints];
225 
226  // new shape # 1 from imin to imax
227  for ( j = 0, i = imin; i != ( imax + 1 ) % nbPoints; i = ( i + 1 ) % nbPoints, j++ )
228  {
229  newShape->x[j] = x[i];
230  newShape->y[j] = y[i];
231  }
232  // is the cutting point a new one ?
233  if ( fps != fpe )
234  {
235  // yes => so add it
236  newShape->x[j] = fptx;
237  newShape->y[j] = fpty;
238  }
239 
240  return newShape;
241 }
242 
243 bool PointSet::containsPoint( double x, double y ) const
244 {
245  GEOSContextHandle_t geosctxt = geosContext();
246  try
247  {
248  GEOSCoordSequence *seq = GEOSCoordSeq_create_r( geosctxt, 1, 2 );
249  GEOSCoordSeq_setX_r( geosctxt, seq, 0, x );
250  GEOSCoordSeq_setY_r( geosctxt, seq, 0, y );
251  geos::unique_ptr point( GEOSGeom_createPoint_r( geosctxt, seq ) );
252  bool result = ( GEOSPreparedContainsProperly_r( geosctxt, preparedGeom(), point.get() ) == 1 );
253 
254  return result;
255  }
256  catch ( GEOSException &e )
257  {
258  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
259  return false;
260  }
261 
262 }
263 
264 bool PointSet::containsLabelCandidate( double x, double y, double width, double height, double alpha ) const
265 {
266  return GeomFunction::containsCandidate( preparedGeom(), x, y, width, height, alpha );
267 }
268 
269 void PointSet::splitPolygons( QLinkedList<PointSet *> &shapes_toProcess,
270  QLinkedList<PointSet *> &shapes_final,
271  double xrm, double yrm )
272 {
273  int i, j;
274 
275  int nbp;
276  double *x = nullptr;
277  double *y = nullptr;
278 
279  int *pts = nullptr;
280 
281  int *cHull = nullptr;
282  int cHullSize;
283 
284  double cp;
285  double bestcp = 0;
286 
287  double bestArea = 0;
288  double area;
289 
290  double base;
291  double b, c;
292  double s;
293 
294  int ihs;
295  int ihn;
296 
297  int ips;
298  int ipn;
299 
300  int holeS = -1; // hole start and end points
301  int holeE = -1;
302 
303  int retainedPt = -1;
304  int pt = 0;
305 
306  double labelArea = xrm * yrm;
307 
308  PointSet *shape = nullptr;
309 
310  while ( !shapes_toProcess.isEmpty() )
311  {
312  shape = shapes_toProcess.takeFirst();
313 
314  x = shape->x;
315  y = shape->y;
316  nbp = shape->nbPoints;
317  pts = new int[nbp];
318 
319  for ( i = 0; i < nbp; i++ )
320  {
321  pts[i] = i;
322  }
323 
324  // conpute convex hull
325  shape->cHullSize = GeomFunction::convexHullId( pts, x, y, nbp, shape->cHull );
326 
327  cHull = shape->cHull;
328  cHullSize = shape->cHullSize;
329 
330  bestArea = 0;
331  retainedPt = -1;
332 
333  // lookup for a hole
334  for ( ihs = 0; ihs < cHullSize; ihs++ )
335  {
336  // ihs->ihn => cHull'seg
337  ihn = ( ihs + 1 ) % cHullSize;
338 
339  ips = cHull[ihs];
340  ipn = ( ips + 1 ) % nbp;
341  if ( ipn != cHull[ihn] ) // next point on shape is not the next point on cHull => there is a hole here !
342  {
343  bestcp = 0;
344  pt = -1;
345  // lookup for the deepest point in the hole
346  for ( i = ips; i != cHull[ihn]; i = ( i + 1 ) % nbp )
347  {
348  cp = std::fabs( GeomFunction::cross_product( x[cHull[ihs]], y[cHull[ihs]],
349  x[cHull[ihn]], y[cHull[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  base = GeomFunction::dist_euc2d( x[cHull[ihs]], y[cHull[ihs]],
362  x[cHull[ihn]], y[cHull[ihn]] );
363 
364  b = GeomFunction::dist_euc2d( x[cHull[ihs]], y[cHull[ihs]],
365  x[pt], y[pt] );
366 
367  c = GeomFunction::dist_euc2d( x[cHull[ihn]], y[cHull[ihn]],
368  x[pt], y[pt] );
369 
370  s = ( base + b + c ) / 2; // s = half perimeter
371  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  c = DBL_MAX;
399 
400  // iterate on all shape points except points which are in the hole
401  bool isValid;
402  int k, l;
403  for ( i = ( cHull[holeE] + 1 ) % nbp; i != ( cHull[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 = GeomFunction::dist_euc2d( 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 = 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] ) ) )
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 = cHull[holeS]; k != cHull[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  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 spliting
504  else if ( imax == imin || nbPtSh1 <= 2 || nbPtSh2 <= 2 || nbPtSh1 == nbp || nbPtSh2 == nbp )
505  {
506  shapes_final.append( shape );
507  }
508  else
509  {
510 
511  PointSet *newShape = shape->extractShape( nbPtSh1, imin, imax, fps, fpe, fptx, fpty );
512 
513  if ( shape->parent )
514  newShape->parent = shape->parent;
515  else
516  newShape->parent = shape;
517 
518  shapes_toProcess.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 );
526 
527  if ( shape->parent )
528  newShape->parent = shape->parent;
529  else
530  newShape->parent = shape;
531 
532  shapes_toProcess.append( newShape );
533 
534  if ( shape->parent )
535  delete shape;
536  }
537  }
538  else
539  {
540  shapes_final.append( shape );
541  }
542  delete[] pts;
543  }
544 }
545 
547 {
548  int i;
549  int j;
550 
551  double bbox[4]; // xmin, ymin, xmax, ymax
552 
553  double alpha;
554  int alpha_d;
555 
556  double alpha_seg;
557 
558  double dref;
559  double d1, d2;
560 
561  double bb[16]; // {ax, ay, bx, by, cx, cy, dx, dy, ex, ey, fx, fy, gx, gy, hx, hy}}
562 
563  double cp;
564  double best_cp;
565  double distNearestPoint;
566 
567  double area;
568  double width;
569  double length;
570 
571  double best_area = DBL_MAX;
572  double best_alpha = -1;
573  double best_bb[16];
574  double best_length = 0;
575  double best_width = 0;
576 
577 
578  bbox[0] = DBL_MAX;
579  bbox[1] = DBL_MAX;
580  bbox[2] = - DBL_MAX;
581  bbox[3] = - DBL_MAX;
582 
583  for ( i = 0; i < cHullSize; i++ )
584  {
585  if ( x[cHull[i]] < bbox[0] )
586  bbox[0] = x[cHull[i]];
587 
588  if ( x[cHull[i]] > bbox[2] )
589  bbox[2] = x[cHull[i]];
590 
591  if ( y[cHull[i]] < bbox[1] )
592  bbox[1] = y[cHull[i]];
593 
594  if ( y[cHull[i]] > bbox[3] )
595  bbox[3] = y[cHull[i]];
596  }
597 
598 
599  dref = bbox[2] - bbox[0];
600 
601  for ( alpha_d = 0; alpha_d < 90; alpha_d++ )
602  {
603  alpha = alpha_d * M_PI / 180.0;
604  d1 = std::cos( alpha ) * dref;
605  d2 = std::sin( alpha ) * dref;
606 
607  bb[0] = bbox[0];
608  bb[1] = bbox[3]; // ax, ay
609 
610  bb[4] = bbox[0];
611  bb[5] = bbox[1]; // cx, cy
612 
613  bb[8] = bbox[2];
614  bb[9] = bbox[1]; // ex, ey
615 
616  bb[12] = bbox[2];
617  bb[13] = bbox[3]; // gx, gy
618 
619 
620  bb[2] = bb[0] + d1;
621  bb[3] = bb[1] + d2; // bx, by
622  bb[6] = bb[4] - d2;
623  bb[7] = bb[5] + d1; // dx, dy
624  bb[10] = bb[8] - d1;
625  bb[11] = bb[9] - d2; // fx, fy
626  bb[14] = bb[12] + d2;
627  bb[15] = bb[13] - d1; // hx, hy
628 
629  // adjust all points
630  for ( i = 0; i < 16; i += 4 )
631  {
632 
633  alpha_seg = ( ( i / 4 > 0 ? ( i / 4 ) - 1 : 3 ) ) * M_PI_2 + alpha;
634 
635  best_cp = DBL_MAX;
636  for ( j = 0; j < nbPoints; j++ )
637  {
638  cp = GeomFunction::cross_product( bb[i + 2], bb[i + 3], bb[i], bb[i + 1], x[cHull[j]], y[cHull[j]] );
639  if ( cp < best_cp )
640  {
641  best_cp = cp;
642  }
643  }
644 
645  distNearestPoint = best_cp / dref;
646 
647  d1 = std::cos( alpha_seg ) * distNearestPoint;
648  d2 = std::sin( alpha_seg ) * distNearestPoint;
649 
650  bb[i] += d1; // x
651  bb[i + 1] += d2; // y
652  bb[i + 2] += d1; // x
653  bb[i + 3] += d2; // y
654  }
655 
656  // compute and compare AREA
657  width = GeomFunction::cross_product( bb[6], bb[7], bb[4], bb[5], bb[12], bb[13] ) / dref;
658  length = GeomFunction::cross_product( bb[2], bb[3], bb[0], bb[1], bb[8], bb[9] ) / dref;
659 
660  area = width * length;
661 
662  if ( area < 0 )
663  area *= -1;
664 
665 
666  if ( best_area - area > EPSILON )
667  {
668  best_area = area;
669  best_length = length;
670  best_width = width;
671  best_alpha = alpha;
672  memcpy( best_bb, bb, sizeof( double ) * 16 );
673  }
674  }
675 
676  // best bbox is defined
677 
678  CHullBox *finalBb = new CHullBox();
679 
680  for ( i = 0; i < 16; i = i + 4 )
681  {
682  GeomFunction::computeLineIntersection( best_bb[i], best_bb[i + 1], best_bb[i + 2], best_bb[i + 3],
683  best_bb[( i + 4 ) % 16], best_bb[( i + 5 ) % 16], best_bb[( i + 6 ) % 16], best_bb[( i + 7 ) % 16],
684  &finalBb->x[int ( i / 4 )], &finalBb->y[int ( i / 4 )] );
685  }
686 
687  finalBb->alpha = best_alpha;
688  finalBb->width = best_width;
689  finalBb->length = best_length;
690 
691  return finalBb;
692 }
693 
694 double PointSet::minDistanceToPoint( double px, double py, double *rx, double *ry ) const
695 {
696  if ( !mGeos )
697  createGeosGeom();
698 
699  if ( !mGeos )
700  return 0;
701 
702  GEOSContextHandle_t geosctxt = geosContext();
703  try
704  {
705  GEOSCoordSequence *coord = GEOSCoordSeq_create_r( geosctxt, 1, 2 );
706  GEOSCoordSeq_setX_r( geosctxt, coord, 0, px );
707  GEOSCoordSeq_setY_r( geosctxt, coord, 0, py );
708  geos::unique_ptr geosPt( GEOSGeom_createPoint_r( geosctxt, coord ) );
709 
710  int type = GEOSGeomTypeId_r( geosctxt, mGeos );
711  const GEOSGeometry *extRing = nullptr;
712  if ( type != GEOS_POLYGON )
713  {
714  extRing = mGeos;
715  }
716  else
717  {
718  //for polygons, we want distance to exterior ring (not an interior point)
719  extRing = GEOSGetExteriorRing_r( geosctxt, mGeos );
720  }
721  geos::coord_sequence_unique_ptr nearestCoord( GEOSNearestPoints_r( geosctxt, extRing, geosPt.get() ) );
722  double nx;
723  double ny;
724  ( void )GEOSCoordSeq_getX_r( geosctxt, nearestCoord.get(), 0, &nx );
725  ( void )GEOSCoordSeq_getY_r( geosctxt, nearestCoord.get(), 0, &ny );
726 
727  if ( rx )
728  *rx = nx;
729  if ( ry )
730  *ry = ny;
731 
732  return GeomFunction::dist_euc2d_sq( px, py, nx, ny );
733  }
734  catch ( GEOSException &e )
735  {
736  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
737  return 0;
738  }
739 }
740 
741 void PointSet::getCentroid( double &px, double &py, bool forceInside ) const
742 {
743  if ( !mGeos )
744  createGeosGeom();
745 
746  if ( !mGeos )
747  return;
748 
749  try
750  {
751  GEOSContextHandle_t geosctxt = geosContext();
752  geos::unique_ptr centroidGeom( GEOSGetCentroid_r( geosctxt, mGeos ) );
753  if ( centroidGeom )
754  {
755  const GEOSCoordSequence *coordSeq = GEOSGeom_getCoordSeq_r( geosctxt, centroidGeom.get() );
756  GEOSCoordSeq_getX_r( geosctxt, coordSeq, 0, &px );
757  GEOSCoordSeq_getY_r( geosctxt, coordSeq, 0, &py );
758  }
759 
760  // check if centroid inside in polygon
761  if ( forceInside && !containsPoint( px, py ) )
762  {
763  geos::unique_ptr pointGeom( GEOSPointOnSurface_r( geosctxt, mGeos ) );
764 
765  if ( pointGeom )
766  {
767  const GEOSCoordSequence *coordSeq = GEOSGeom_getCoordSeq_r( geosctxt, pointGeom.get() );
768  GEOSCoordSeq_getX_r( geosctxt, coordSeq, 0, &px );
769  GEOSCoordSeq_getY_r( geosctxt, coordSeq, 0, &py );
770  }
771  }
772  }
773  catch ( GEOSException &e )
774  {
775  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
776  return;
777  }
778 }
779 
780 void PointSet::getPointByDistance( double *d, double *ad, double dl, double *px, double *py )
781 {
782  int i;
783  double dx, dy, di;
784  double distr;
785 
786  i = 0;
787  if ( dl >= 0 )
788  {
789  while ( i < nbPoints && ad[i] <= dl ) i++;
790  i--;
791  }
792 
793  if ( i < nbPoints - 1 )
794  {
795  if ( dl < 0 )
796  {
797  dx = x[nbPoints - 1] - x[0];
798  dy = y[nbPoints - 1] - y[0];
799  di = std::sqrt( dx * dx + dy * dy );
800  }
801  else
802  {
803  dx = x[i + 1] - x[i];
804  dy = y[i + 1] - y[i];
805  di = d[i];
806  }
807 
808  distr = dl - ad[i];
809  *px = x[i] + dx * distr / di;
810  *py = y[i] + dy * distr / di;
811  }
812  else // just select last point...
813  {
814  *px = x[i];
815  *py = y[i];
816  }
817 }
818 
819 const GEOSGeometry *PointSet::geos() const
820 {
821  if ( !mGeos )
822  createGeosGeom();
823 
824  return mGeos;
825 }
826 
827 double PointSet::length() const
828 {
829  if ( !mGeos )
830  createGeosGeom();
831 
832  if ( !mGeos )
833  return -1;
834 
835  GEOSContextHandle_t geosctxt = geosContext();
836 
837  try
838  {
839  double len = 0;
840  ( void )GEOSLength_r( geosctxt, mGeos, &len );
841  return len;
842  }
843  catch ( GEOSException &e )
844  {
845  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
846  return -1;
847  }
848 }
849 
850 bool PointSet::isClosed() const
851 {
852  return qgsDoubleNear( x[0], x[nbPoints - 1] ) && qgsDoubleNear( y[0], y[nbPoints - 1] );
853 }
double length
Definition: pointset.h:60
void invalidateGeos()
Definition: pointset.cpp:179
bool containsPoint(double x, double y) const
Tests whether point set contains a specified point.
Definition: pointset.cpp:243
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:215
void createGeosGeom() const
Definition: pointset.cpp:125
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:850
struct pal::_cHullBox CHullBox
bool qgsDoubleNear(double a, double b, double epsilon=4 *DBL_EPSILON)
Compare two doubles (but allow some difference)
Definition: qgis.h:251
static void logMessage(const QString &message, const QString &tag=QString(), Qgis::MessageLevel level=Qgis::Warning)
add a message to the instance (and create it if necessary)
static double cross_product(double x1, double y1, double x2, double y2, double x3, double y3)
Definition: geomfunction.h:61
double width
Definition: pointset.h:59
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)
Get a point a set distance along a line geometry.
Definition: pointset.cpp:780
PointSet * parent
Definition: pointset.h:178
virtual ~PointSet()
Definition: pointset.cpp:190
double * x
Definition: pointset.h:169
double ymax
Definition: pointset.h:192
double xmin
Definition: pointset.h:189
PointSet * holeOf
Definition: pointset.h:177
std::unique_ptr< GEOSGeometry, GeosDeleter > unique_ptr
Scoped GEOS pointer.
Definition: qgsgeos.h:74
double ymin
Definition: pointset.h:191
int * cHull
Definition: pointset.h:172
static int convexHullId(int *id, const double *const x, const double *const y, int n, int *&cHull)
Compute the convex hull in O(n·log(n))
static double dist_euc2d(double x1, double y1, double x2, double y2)
Definition: geomfunction.h:66
GEOSContextHandle_t geosContext()
Get GEOS context handle to be used in all GEOS library calls with reentrant API.
Definition: pal.cpp:48
void deleteCoords()
Definition: pointset.cpp:207
double * y
Definition: pointset.h:170
CHullBox * compute_chull_bbox()
Definition: pointset.cpp:546
double x[4]
Definition: pointset.h:54
double length() const
Returns length of line geometry.
Definition: pointset.cpp:827
double y[4]
Definition: pointset.h:55
#define EPSILON
Definition: util.h:75
GEOSGeometry * mGeos
Definition: pointset.h:165
const GEOSPreparedGeometry * preparedGeom() const
Definition: pointset.cpp:167
double alpha
Definition: pointset.h:57
std::unique_ptr< GEOSCoordSequence, GeosDeleter > coord_sequence_unique_ptr
Scoped GEOS coordinate sequence pointer.
Definition: qgsgeos.h:89
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:741
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:269
const GEOSGeometry * geos() const
Returns the point set&#39;s GEOS geometry.
Definition: pointset.cpp:819
double xmax
Definition: pointset.h:190
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:694
bool mOwnsGeom
Definition: pointset.h:166
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:264