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