QGIS API Documentation  3.4.15-Madeira (e83d02e274)
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( QgsGeos::getGEOSHandler(), ps.mGeos );
121  mOwnsGeom = true;
122  }
123 }
124 
126 {
127  GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
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( QgsGeos::getGEOSHandler(), mGeos );
175  }
176  return mPreparedGeom;
177 }
178 
180 {
181  GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
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 = QgsGeos::getGEOSHandler();
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  delete[] cHull;
204 }
205 
207 {
208  delete[] x;
209  delete[] y;
210  x = nullptr;
211  y = nullptr;
212 }
213 
214 PointSet *PointSet::extractShape( int nbPtSh, int imin, int imax, int fps, int fpe, double fptx, double fpty )
215 {
216 
217  int i, j;
218 
219  PointSet *newShape = new PointSet();
220  newShape->type = GEOS_POLYGON;
221  newShape->nbPoints = nbPtSh;
222  newShape->x = new double[newShape->nbPoints];
223  newShape->y = new double[newShape->nbPoints];
224 
225  // new shape # 1 from imin to imax
226  for ( j = 0, i = imin; i != ( imax + 1 ) % nbPoints; i = ( i + 1 ) % nbPoints, j++ )
227  {
228  newShape->x[j] = x[i];
229  newShape->y[j] = y[i];
230  }
231  // is the cutting point a new one ?
232  if ( fps != fpe )
233  {
234  // yes => so add it
235  newShape->x[j] = fptx;
236  newShape->y[j] = fpty;
237  }
238 
239  return newShape;
240 }
241 
242 bool PointSet::containsPoint( double x, double y ) const
243 {
244  GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
245  try
246  {
247  GEOSCoordSequence *seq = GEOSCoordSeq_create_r( geosctxt, 1, 2 );
248  GEOSCoordSeq_setX_r( geosctxt, seq, 0, x );
249  GEOSCoordSeq_setY_r( geosctxt, seq, 0, y );
250  geos::unique_ptr point( GEOSGeom_createPoint_r( geosctxt, seq ) );
251  bool result = ( GEOSPreparedContainsProperly_r( geosctxt, preparedGeom(), point.get() ) == 1 );
252 
253  return result;
254  }
255  catch ( GEOSException &e )
256  {
257  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
258  return false;
259  }
260 
261 }
262 
263 bool PointSet::containsLabelCandidate( double x, double y, double width, double height, double alpha ) const
264 {
265  return GeomFunction::containsCandidate( preparedGeom(), x, y, width, height, alpha );
266 }
267 
268 void PointSet::splitPolygons( QLinkedList<PointSet *> &shapes_toProcess,
269  QLinkedList<PointSet *> &shapes_final,
270  double xrm, double yrm )
271 {
272  int i, j;
273 
274  int nbp;
275  double *x = nullptr;
276  double *y = nullptr;
277 
278  int *pts = nullptr;
279 
280  int *cHull = nullptr;
281  int cHullSize;
282 
283  double cp;
284  double bestcp = 0;
285 
286  double bestArea = 0;
287  double area;
288 
289  double base;
290  double b, c;
291  double s;
292 
293  int ihs;
294  int ihn;
295 
296  int ips;
297  int ipn;
298 
299  int holeS = -1; // hole start and end points
300  int holeE = -1;
301 
302  int retainedPt = -1;
303  int pt = 0;
304 
305  double labelArea = xrm * yrm;
306 
307  PointSet *shape = nullptr;
308 
309  while ( !shapes_toProcess.isEmpty() )
310  {
311  shape = shapes_toProcess.takeFirst();
312 
313  x = shape->x;
314  y = shape->y;
315  nbp = shape->nbPoints;
316  pts = new int[nbp];
317 
318  for ( i = 0; i < nbp; i++ )
319  {
320  pts[i] = i;
321  }
322 
323  // conpute convex hull
324  shape->cHullSize = GeomFunction::convexHullId( pts, x, y, nbp, shape->cHull );
325 
326  cHull = shape->cHull;
327  cHullSize = shape->cHullSize;
328 
329  bestArea = 0;
330  retainedPt = -1;
331 
332  // lookup for a hole
333  for ( ihs = 0; ihs < cHullSize; ihs++ )
334  {
335  // ihs->ihn => cHull'seg
336  ihn = ( ihs + 1 ) % cHullSize;
337 
338  ips = cHull[ihs];
339  ipn = ( ips + 1 ) % nbp;
340  if ( ipn != cHull[ihn] ) // next point on shape is not the next point on cHull => there is a hole here !
341  {
342  bestcp = 0;
343  pt = -1;
344  // lookup for the deepest point in the hole
345  for ( i = ips; i != cHull[ihn]; i = ( i + 1 ) % nbp )
346  {
347  cp = std::fabs( GeomFunction::cross_product( x[cHull[ihs]], y[cHull[ihs]],
348  x[cHull[ihn]], y[cHull[ihn]],
349  x[i], y[i] ) );
350  if ( cp - bestcp > EPSILON )
351  {
352  bestcp = cp;
353  pt = i;
354  }
355  }
356 
357  if ( pt != -1 )
358  {
359  // compute the ihs->ihn->pt triangle's area
360  base = GeomFunction::dist_euc2d( x[cHull[ihs]], y[cHull[ihs]],
361  x[cHull[ihn]], y[cHull[ihn]] );
362 
363  b = GeomFunction::dist_euc2d( x[cHull[ihs]], y[cHull[ihs]],
364  x[pt], y[pt] );
365 
366  c = GeomFunction::dist_euc2d( x[cHull[ihn]], y[cHull[ihn]],
367  x[pt], y[pt] );
368 
369  s = ( base + b + c ) / 2; // s = half perimeter
370  area = s * ( s - base ) * ( s - b ) * ( s - c );
371  if ( area < 0 )
372  area = -area;
373 
374  // retain the biggest area
375  if ( area - bestArea > EPSILON )
376  {
377  bestArea = area;
378  retainedPt = pt;
379  holeS = ihs;
380  holeE = ihn;
381  }
382  }
383  }
384  }
385 
386  // we have a hole, its area, and the deppest point in hole
387  // we're going to find the second point to cup the shape
388  // holeS = hole starting point
389  // holeE = hole ending point
390  // retainedPt = deppest point in hole
391  // bestArea = area of triangle HoleS->holeE->retainedPoint
392  bestArea = std::sqrt( bestArea );
393  double cx, cy, dx, dy, ex, ey, fx, fy, seg_length, ptx = 0, pty = 0, fptx = 0, fpty = 0;
394  int ps = -1, pe = -1, fps = -1, fpe = -1;
395  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)
396  {
397  c = std::numeric_limits<double>::max();
398 
399  // iterate on all shape points except points which are in the hole
400  bool isValid;
401  int k, l;
402  for ( i = ( cHull[holeE] + 1 ) % nbp; i != ( cHull[holeS] - 1 + nbp ) % nbp; i = j )
403  {
404  j = ( i + 1 ) % nbp; // i->j is shape segment not in hole
405 
406  // compute distance between retainedPoint and segment
407  // whether perpendicular distance (if retaindPoint is fronting segment i->j)
408  // or distance between retainedPt and i or j (choose the nearest)
409  seg_length = GeomFunction::dist_euc2d( x[i], y[i], x[j], y[j] );
410  cx = ( x[i] + x[j] ) / 2.0;
411  cy = ( y[i] + y[j] ) / 2.0;
412  dx = cy - y[i];
413  dy = cx - x[i];
414 
415  ex = cx - dx;
416  ey = cy + dy;
417  fx = cx + dx;
418  fy = cy - dy;
419 
420  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
421  {
422  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] ) ) )
423  {
424  b = ex;
425  ps = i;
426  pe = i;
427  }
428  else
429  {
430  b = ey;
431  ps = j;
432  pe = j;
433  }
434  }
435  else // point fronting i->j => compute pependicular distance => create a new point
436  {
437  b = GeomFunction::cross_product( x[i], y[i], x[j], y[j], x[retainedPt], y[retainedPt] ) / seg_length;
438  b *= b;
439  ps = i;
440  pe = j;
441 
442  if ( !GeomFunction::computeLineIntersection( x[i], y[i], x[j], y[j], x[retainedPt], y[retainedPt], x[retainedPt] - dx, y[retainedPt] + dy, &ptx, &pty ) )
443  {
444  //error - it should intersect the line
445  }
446  }
447 
448  isValid = true;
449  double pointX, pointY;
450  if ( ps == pe )
451  {
452  pointX = x[pe];
453  pointY = y[pe];
454  }
455  else
456  {
457  pointX = ptx;
458  pointY = pty;
459  }
460 
461  for ( k = cHull[holeS]; k != cHull[holeE]; k = ( k + 1 ) % nbp )
462  {
463  l = ( k + 1 ) % nbp;
464  if ( GeomFunction::isSegIntersects( x[retainedPt], y[retainedPt], pointX, pointY, x[k], y[k], x[l], y[l] ) )
465  {
466  isValid = false;
467  break;
468  }
469  }
470 
471 
472  if ( isValid && b < c )
473  {
474  c = b;
475  fps = ps;
476  fpe = pe;
477  fptx = ptx;
478  fpty = pty;
479  }
480  } // for point which are not in hole
481 
482  // we will cut the shapeu in two new shapes, one from [retainedPoint] to [newPoint] and one form [newPoint] to [retainedPoint]
483  int imin = retainedPt;
484  int imax = ( ( ( fps < retainedPt && fpe < retainedPt ) || ( fps > retainedPt && fpe > retainedPt ) ) ? std::min( fps, fpe ) : std::max( fps, fpe ) );
485 
486  int nbPtSh1, nbPtSh2; // how many points in new shapes ?
487  if ( imax > imin )
488  nbPtSh1 = imax - imin + 1 + ( fpe != fps );
489  else
490  nbPtSh1 = imax + nbp - imin + 1 + ( fpe != fps );
491 
492  if ( ( imax == fps ? fpe : fps ) < imin )
493  nbPtSh2 = imin - ( imax == fps ? fpe : fps ) + 1 + ( fpe != fps );
494  else
495  nbPtSh2 = imin + nbp - ( imax == fps ? fpe : fps ) + 1 + ( fpe != fps );
496 
497  if ( retainedPt == -1 || fps == -1 || fpe == -1 )
498  {
499  if ( shape->parent )
500  delete shape;
501  }
502  // check for useless splitting
503  else if ( imax == imin || nbPtSh1 <= 2 || nbPtSh2 <= 2 || nbPtSh1 == nbp || nbPtSh2 == nbp )
504  {
505  shapes_final.append( shape );
506  }
507  else
508  {
509 
510  PointSet *newShape = shape->extractShape( nbPtSh1, imin, imax, fps, fpe, fptx, fpty );
511 
512  if ( shape->parent )
513  newShape->parent = shape->parent;
514  else
515  newShape->parent = shape;
516 
517  shapes_toProcess.append( newShape );
518 
519  if ( imax == fps )
520  imax = fpe;
521  else
522  imax = fps;
523 
524  newShape = shape->extractShape( nbPtSh2, imax, imin, fps, fpe, fptx, fpty );
525 
526  if ( shape->parent )
527  newShape->parent = shape->parent;
528  else
529  newShape->parent = shape;
530 
531  shapes_toProcess.append( newShape );
532 
533  if ( shape->parent )
534  delete shape;
535  }
536  }
537  else
538  {
539  shapes_final.append( shape );
540  }
541  delete[] pts;
542  }
543 }
544 
546 {
547  int i;
548  int j;
549 
550  double bbox[4]; // xmin, ymin, xmax, ymax
551 
552  double alpha;
553  int alpha_d;
554 
555  double alpha_seg;
556 
557  double dref;
558  double d1, d2;
559 
560  double bb[16]; // {ax, ay, bx, by, cx, cy, dx, dy, ex, ey, fx, fy, gx, gy, hx, hy}}
561 
562  double cp;
563  double best_cp;
564  double distNearestPoint;
565 
566  double area;
567  double width;
568  double length;
569 
570  double best_area = std::numeric_limits<double>::max();
571  double best_alpha = -1;
572  double best_bb[16];
573  double best_length = 0;
574  double best_width = 0;
575 
576 
577  bbox[0] = std::numeric_limits<double>::max();
578  bbox[1] = std::numeric_limits<double>::max();
579  bbox[2] = std::numeric_limits<double>::lowest();
580  bbox[3] = std::numeric_limits<double>::lowest();
581 
582  for ( i = 0; i < cHullSize; i++ )
583  {
584  if ( x[cHull[i]] < bbox[0] )
585  bbox[0] = x[cHull[i]];
586 
587  if ( x[cHull[i]] > bbox[2] )
588  bbox[2] = x[cHull[i]];
589 
590  if ( y[cHull[i]] < bbox[1] )
591  bbox[1] = y[cHull[i]];
592 
593  if ( y[cHull[i]] > bbox[3] )
594  bbox[3] = y[cHull[i]];
595  }
596 
597 
598  dref = bbox[2] - bbox[0];
599 
600  for ( alpha_d = 0; alpha_d < 90; alpha_d++ )
601  {
602  alpha = alpha_d * M_PI / 180.0;
603  d1 = std::cos( alpha ) * dref;
604  d2 = std::sin( alpha ) * dref;
605 
606  bb[0] = bbox[0];
607  bb[1] = bbox[3]; // ax, ay
608 
609  bb[4] = bbox[0];
610  bb[5] = bbox[1]; // cx, cy
611 
612  bb[8] = bbox[2];
613  bb[9] = bbox[1]; // ex, ey
614 
615  bb[12] = bbox[2];
616  bb[13] = bbox[3]; // gx, gy
617 
618 
619  bb[2] = bb[0] + d1;
620  bb[3] = bb[1] + d2; // bx, by
621  bb[6] = bb[4] - d2;
622  bb[7] = bb[5] + d1; // dx, dy
623  bb[10] = bb[8] - d1;
624  bb[11] = bb[9] - d2; // fx, fy
625  bb[14] = bb[12] + d2;
626  bb[15] = bb[13] - d1; // hx, hy
627 
628  // adjust all points
629  for ( i = 0; i < 16; i += 4 )
630  {
631 
632  alpha_seg = ( ( i / 4 > 0 ? ( i / 4 ) - 1 : 3 ) ) * M_PI_2 + alpha;
633 
634  best_cp = std::numeric_limits<double>::max();
635  for ( j = 0; j < nbPoints; j++ )
636  {
637  cp = GeomFunction::cross_product( bb[i + 2], bb[i + 3], bb[i], bb[i + 1], x[cHull[j]], y[cHull[j]] );
638  if ( cp < best_cp )
639  {
640  best_cp = cp;
641  }
642  }
643 
644  distNearestPoint = best_cp / dref;
645 
646  d1 = std::cos( alpha_seg ) * distNearestPoint;
647  d2 = std::sin( alpha_seg ) * distNearestPoint;
648 
649  bb[i] += d1; // x
650  bb[i + 1] += d2; // y
651  bb[i + 2] += d1; // x
652  bb[i + 3] += d2; // y
653  }
654 
655  // compute and compare AREA
656  width = GeomFunction::cross_product( bb[6], bb[7], bb[4], bb[5], bb[12], bb[13] ) / dref;
657  length = GeomFunction::cross_product( bb[2], bb[3], bb[0], bb[1], bb[8], bb[9] ) / dref;
658 
659  area = width * length;
660 
661  if ( area < 0 )
662  area *= -1;
663 
664 
665  if ( best_area - area > EPSILON )
666  {
667  best_area = area;
668  best_length = length;
669  best_width = width;
670  best_alpha = alpha;
671  memcpy( best_bb, bb, sizeof( double ) * 16 );
672  }
673  }
674 
675  // best bbox is defined
676 
677  CHullBox *finalBb = new CHullBox();
678 
679  for ( i = 0; i < 16; i = i + 4 )
680  {
681  GeomFunction::computeLineIntersection( best_bb[i], best_bb[i + 1], best_bb[i + 2], best_bb[i + 3],
682  best_bb[( i + 4 ) % 16], best_bb[( i + 5 ) % 16], best_bb[( i + 6 ) % 16], best_bb[( i + 7 ) % 16],
683  &finalBb->x[int ( i / 4 )], &finalBb->y[int ( i / 4 )] );
684  }
685 
686  finalBb->alpha = best_alpha;
687  finalBb->width = best_width;
688  finalBb->length = best_length;
689 
690  return finalBb;
691 }
692 
693 double PointSet::minDistanceToPoint( double px, double py, double *rx, double *ry ) const
694 {
695  if ( !mGeos )
696  createGeosGeom();
697 
698  if ( !mGeos )
699  return 0;
700 
701  GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
702  try
703  {
704  GEOSCoordSequence *coord = GEOSCoordSeq_create_r( geosctxt, 1, 2 );
705  GEOSCoordSeq_setX_r( geosctxt, coord, 0, px );
706  GEOSCoordSeq_setY_r( geosctxt, coord, 0, py );
707  geos::unique_ptr geosPt( GEOSGeom_createPoint_r( geosctxt, coord ) );
708 
709  int type = GEOSGeomTypeId_r( geosctxt, mGeos );
710  const GEOSGeometry *extRing = nullptr;
711  if ( type != GEOS_POLYGON )
712  {
713  extRing = mGeos;
714  }
715  else
716  {
717  //for polygons, we want distance to exterior ring (not an interior point)
718  extRing = GEOSGetExteriorRing_r( geosctxt, mGeos );
719  }
720  geos::coord_sequence_unique_ptr nearestCoord( GEOSNearestPoints_r( geosctxt, extRing, geosPt.get() ) );
721  double nx;
722  double ny;
723  ( void )GEOSCoordSeq_getX_r( geosctxt, nearestCoord.get(), 0, &nx );
724  ( void )GEOSCoordSeq_getY_r( geosctxt, nearestCoord.get(), 0, &ny );
725 
726  if ( rx )
727  *rx = nx;
728  if ( ry )
729  *ry = ny;
730 
731  return GeomFunction::dist_euc2d_sq( px, py, nx, ny );
732  }
733  catch ( GEOSException &e )
734  {
735  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
736  return 0;
737  }
738 }
739 
740 void PointSet::getCentroid( double &px, double &py, bool forceInside ) const
741 {
742  if ( !mGeos )
743  createGeosGeom();
744 
745  if ( !mGeos )
746  return;
747 
748  try
749  {
750  GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
751  geos::unique_ptr centroidGeom( GEOSGetCentroid_r( geosctxt, mGeos ) );
752  if ( centroidGeom )
753  {
754  const GEOSCoordSequence *coordSeq = GEOSGeom_getCoordSeq_r( geosctxt, centroidGeom.get() );
755  GEOSCoordSeq_getX_r( geosctxt, coordSeq, 0, &px );
756  GEOSCoordSeq_getY_r( geosctxt, coordSeq, 0, &py );
757  }
758 
759  // check if centroid inside in polygon
760  if ( forceInside && !containsPoint( px, py ) )
761  {
762  geos::unique_ptr pointGeom( GEOSPointOnSurface_r( geosctxt, mGeos ) );
763 
764  if ( pointGeom )
765  {
766  const GEOSCoordSequence *coordSeq = GEOSGeom_getCoordSeq_r( geosctxt, pointGeom.get() );
767  GEOSCoordSeq_getX_r( geosctxt, coordSeq, 0, &px );
768  GEOSCoordSeq_getY_r( geosctxt, coordSeq, 0, &py );
769  }
770  }
771  }
772  catch ( GEOSException &e )
773  {
774  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
775  return;
776  }
777 }
778 
779 void PointSet::getPointByDistance( double *d, double *ad, double dl, double *px, double *py )
780 {
781  int i;
782  double dx, dy, di;
783  double distr;
784 
785  i = 0;
786  if ( dl >= 0 )
787  {
788  while ( i < nbPoints && ad[i] <= dl ) i++;
789  i--;
790  }
791 
792  if ( i < nbPoints - 1 )
793  {
794  if ( dl < 0 )
795  {
796  dx = x[nbPoints - 1] - x[0];
797  dy = y[nbPoints - 1] - y[0];
798  di = std::sqrt( dx * dx + dy * dy );
799  }
800  else
801  {
802  dx = x[i + 1] - x[i];
803  dy = y[i + 1] - y[i];
804  di = d[i];
805  }
806 
807  distr = dl - ad[i];
808  *px = x[i] + dx * distr / di;
809  *py = y[i] + dy * distr / di;
810  }
811  else // just select last point...
812  {
813  *px = x[i];
814  *py = y[i];
815  }
816 }
817 
818 const GEOSGeometry *PointSet::geos() const
819 {
820  if ( !mGeos )
821  createGeosGeom();
822 
823  return mGeos;
824 }
825 
826 double PointSet::length() const
827 {
828  if ( !mGeos )
829  createGeosGeom();
830 
831  if ( !mGeos )
832  return -1;
833 
834  GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
835 
836  try
837  {
838  double len = 0;
839  ( void )GEOSLength_r( geosctxt, mGeos, &len );
840  return len;
841  }
842  catch ( GEOSException &e )
843  {
844  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
845  return -1;
846  }
847 }
848 
849 bool PointSet::isClosed() const
850 {
851  return qgsDoubleNear( x[0], x[nbPoints - 1] ) && qgsDoubleNear( y[0], y[nbPoints - 1] );
852 }
double length
Definition: pointset.h:60
void invalidateGeos()
Definition: pointset.cpp:179
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:214
bool qgsDoubleNear(double a, double b, double epsilon=4 *std::numeric_limits< double >::epsilon())
Compare two doubles (but allow some difference)
Definition: qgis.h:278
const GEOSGeometry * geos() const
Returns the point set&#39;s GEOS geometry.
Definition: pointset.cpp:818
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.
void createGeosGeom() const
Definition: pointset.cpp:125
struct pal::_cHullBox CHullBox
void getCentroid(double &px, double &py, bool forceInside=false) const
Definition: pointset.cpp:740
static GEOSContextHandle_t getGEOSHandler()
Definition: qgsgeos.cpp:2821
static double cross_product(double x1, double y1, double x2, double y2, double x3, double y3)
Definition: geomfunction.h:61
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:59
bool isClosed() const
Returns true if pointset is closed.
Definition: pointset.cpp:849
const GEOSPreparedGeometry * preparedGeom() const
Definition: pointset.cpp:167
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:779
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: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:191
int * cHull
Definition: pointset.h:172
double length() const
Returns length of line geometry.
Definition: pointset.cpp:826
static double dist_euc2d(double x1, double y1, double x2, double y2)
Definition: geomfunction.h:66
void deleteCoords()
Definition: pointset.cpp:206
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:263
bool containsPoint(double x, double y) const
Tests whether point set contains a specified point.
Definition: pointset.cpp:242
double * y
Definition: pointset.h:170
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:268
CHullBox * compute_chull_bbox()
Definition: pointset.cpp:545
double x[4]
Definition: pointset.h:54
double y[4]
Definition: pointset.h:55
#define EPSILON
Definition: util.h:75
GEOSGeometry * mGeos
Definition: pointset.h:165
double alpha
Definition: pointset.h:57
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
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:693
bool mOwnsGeom
Definition: pointset.h:166
static int convexHullId(int *id, const double *x, const double *y, int n, int *&cHull)
Compute the convex hull in O(n┬Ělog(n))