QGIS API Documentation  3.6.0-Noosa (5873452)
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
bool containsPoint(double x, double y) const
Tests whether point set contains a specified point.
Definition: pointset.cpp:242
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:265
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:849
struct pal::_cHullBox CHullBox
static GEOSContextHandle_t getGEOSHandler()
Definition: qgsgeos.cpp:2777
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
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
static double dist_euc2d(double x1, double y1, double x2, double y2)
Definition: geomfunction.h:66
void deleteCoords()
Definition: pointset.cpp:206
double * y
Definition: pointset.h:170
CHullBox * compute_chull_bbox()
Definition: pointset.cpp:545
double x[4]
Definition: pointset.h:54
double length() const
Returns length of line geometry.
Definition: pointset.cpp:826
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: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:740
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
const GEOSGeometry * geos() const
Returns the point set&#39;s GEOS geometry.
Definition: pointset.cpp:818
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))
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