QGIS API Documentation  2.14.0-Essen
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  GEOSContextHandle_t geosctxt = geosContext();
293  GEOSCoordSequence *coord = GEOSCoordSeq_create_r( geosctxt, 5, 2 );
294 
295  GEOSCoordSeq_setX_r( geosctxt, coord, 0, x );
296  GEOSCoordSeq_setY_r( geosctxt, coord, 0, y );
297  if ( !qgsDoubleNear( alpha, 0.0 ) )
298  {
299  double beta = alpha + ( M_PI / 2 );
300  double dx1 = cos( alpha ) * width;
301  double dy1 = sin( alpha ) * width;
302  double dx2 = cos( beta ) * height;
303  double dy2 = sin( beta ) * height;
304  GEOSCoordSeq_setX_r( geosctxt, coord, 1, x + dx1 );
305  GEOSCoordSeq_setY_r( geosctxt, coord, 1, y + dy1 );
306  GEOSCoordSeq_setX_r( geosctxt, coord, 2, x + dx1 + dx2 );
307  GEOSCoordSeq_setY_r( geosctxt, coord, 2, y + dy1 + dy2 );
308  GEOSCoordSeq_setX_r( geosctxt, coord, 3, x + dx2 );
309  GEOSCoordSeq_setY_r( geosctxt, coord, 3, y + dy2 );
310  }
311  else
312  {
313  GEOSCoordSeq_setX_r( geosctxt, coord, 1, x + width );
314  GEOSCoordSeq_setY_r( geosctxt, coord, 1, y );
315  GEOSCoordSeq_setX_r( geosctxt, coord, 2, x + width );
316  GEOSCoordSeq_setY_r( geosctxt, coord, 2, y + height );
317  GEOSCoordSeq_setX_r( geosctxt, coord, 3, x );
318  GEOSCoordSeq_setY_r( geosctxt, coord, 3, y + height );
319  }
320  //close ring
321  GEOSCoordSeq_setX_r( geosctxt, coord, 4, x );
322  GEOSCoordSeq_setY_r( geosctxt, coord, 4, y );
323 
324  try
325  {
326  GEOSGeometry* bboxGeos = GEOSGeom_createLinearRing_r( geosctxt, coord );
327  bool result = ( GEOSPreparedContainsProperly_r( geosctxt, preparedGeom(), bboxGeos ) == 1 );
328  GEOSGeom_destroy_r( geosctxt, bboxGeos );
329  return result;
330  }
331  catch ( GEOSException &e )
332  {
333  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
334  return false;
335  }
336 }
337 
339  QLinkedList<PointSet*> &shapes_final,
340  double xrm, double yrm )
341 {
342  int i, j;
343 
344  int nbp;
345  double *x;
346  double *y;
347 
348  int *pts;
349 
350  int *cHull;
351  int cHullSize;
352 
353  double cp;
354  double bestcp = 0;
355 
356  double bestArea = 0;
357  double area;
358 
359  double base;
360  double b, c;
361  double s;
362 
363  int ihs;
364  int ihn;
365 
366  int ips;
367  int ipn;
368 
369  int holeS = -1; // hole start and end points
370  int holeE = -1;
371 
372  int retainedPt = -1;
373  int pt = 0;
374 
375  double labelArea = xrm * yrm;
376 
377  PointSet *shape;
378 
379  while ( !shapes_toProcess.isEmpty() )
380  {
381  shape = shapes_toProcess.takeFirst();
382 
383  x = shape->x;
384  y = shape->y;
385  nbp = shape->nbPoints;
386  pts = new int[nbp];
387 
388  for ( i = 0; i < nbp; i++ )
389  {
390  pts[i] = i;
391  }
392 
393  // conpute convex hull
394  shape->cHullSize = GeomFunction::convexHullId( pts, x, y, nbp, shape->cHull );
395 
396  cHull = shape->cHull;
397  cHullSize = shape->cHullSize;
398 
399  bestArea = 0;
400  retainedPt = -1;
401 
402  // lookup for a hole
403  for ( ihs = 0; ihs < cHullSize; ihs++ )
404  {
405  // ihs->ihn => cHull'seg
406  ihn = ( ihs + 1 ) % cHullSize;
407 
408  ips = cHull[ihs];
409  ipn = ( ips + 1 ) % nbp;
410  if ( ipn != cHull[ihn] ) // next point on shape is not the next point on cHull => there is a hole here !
411  {
412  bestcp = 0;
413  pt = -1;
414  // lookup for the deepest point in the hole
415  for ( i = ips; i != cHull[ihn]; i = ( i + 1 ) % nbp )
416  {
417  cp = qAbs( GeomFunction::cross_product( x[cHull[ihs]], y[cHull[ihs]],
418  x[cHull[ihn]], y[cHull[ihn]],
419  x[i], y[i] ) );
420  if ( cp - bestcp > EPSILON )
421  {
422  bestcp = cp;
423  pt = i;
424  }
425  }
426 
427  if ( pt != -1 )
428  {
429  // compute the ihs->ihn->pt triangle's area
430  base = GeomFunction::dist_euc2d( x[cHull[ihs]], y[cHull[ihs]],
431  x[cHull[ihn]], y[cHull[ihn]] );
432 
433  b = GeomFunction::dist_euc2d( x[cHull[ihs]], y[cHull[ihs]],
434  x[pt], y[pt] );
435 
436  c = GeomFunction::dist_euc2d( x[cHull[ihn]], y[cHull[ihn]],
437  x[pt], y[pt] );
438 
439  s = ( base + b + c ) / 2; // s = half perimeter
440  area = s * ( s - base ) * ( s - b ) * ( s - c );
441  if ( area < 0 )
442  area = -area;
443 
444  // retain the biggest area
445  if ( area - bestArea > EPSILON )
446  {
447  bestArea = area;
448  retainedPt = pt;
449  holeS = ihs;
450  holeE = ihn;
451  }
452  }
453  }
454  }
455 
456  // we have a hole, its area, and the deppest point in hole
457  // we're going to find the second point to cup the shape
458  // holeS = hole starting point
459  // holeE = hole ending point
460  // retainedPt = deppest point in hole
461  // bestArea = area of triangle HoleS->holeE->retainedPoint
462  bestArea = sqrt( bestArea );
463  double cx, cy, dx, dy, ex, ey, fx, fy, seg_length, ptx = 0, pty = 0, fptx = 0, fpty = 0;
464  int ps = -1, pe = -1, fps = -1, fpe = -1;
465  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)
466  {
467  c = DBL_MAX;
468 
469  // iterate on all shape points except points which are in the hole
470  bool isValid;
471  int k, l;
472  for ( i = ( cHull[holeE] + 1 ) % nbp; i != ( cHull[holeS] - 1 + nbp ) % nbp; i = j )
473  {
474  j = ( i + 1 ) % nbp; // i->j is shape segment not in hole
475 
476  // compute distance between retainedPoint and segment
477  // whether perpendicular distance (if retaindPoint is fronting segment i->j)
478  // or distance between retainedPt and i or j (choose the nearest)
479  seg_length = GeomFunction::dist_euc2d( x[i], y[i], x[j], y[j] );
480  cx = ( x[i] + x[j] ) / 2.0;
481  cy = ( y[i] + y[j] ) / 2.0;
482  dx = cy - y[i];
483  dy = cx - x[i];
484 
485  ex = cx - dx;
486  ey = cy + dy;
487  fx = cx + dx;
488  fy = cy - dy;
489 
490  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
491  {
492  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] ) ) )
493  {
494  b = ex;
495  ps = i;
496  pe = i;
497  }
498  else
499  {
500  b = ey;
501  ps = j;
502  pe = j;
503  }
504  }
505  else // point fronting i->j => compute pependicular distance => create a new point
506  {
507  b = GeomFunction::cross_product( x[i], y[i], x[j], y[j], x[retainedPt], y[retainedPt] ) / seg_length;
508  b *= b;
509  ps = i;
510  pe = j;
511 
512  if ( !GeomFunction::computeLineIntersection( x[i], y[i], x[j], y[j], x[retainedPt], y[retainedPt], x[retainedPt] - dx, y[retainedPt] + dy, &ptx, &pty ) )
513  {
514  //error - it should intersect the line
515  }
516  }
517 
518  isValid = true;
519  double pointX, pointY;
520  if ( ps == pe )
521  {
522  pointX = x[pe];
523  pointY = y[pe];
524  }
525  else
526  {
527  pointX = ptx;
528  pointY = pty;
529  }
530 
531  for ( k = cHull[holeS]; k != cHull[holeE]; k = ( k + 1 ) % nbp )
532  {
533  l = ( k + 1 ) % nbp;
534  if ( GeomFunction::isSegIntersects( x[retainedPt], y[retainedPt], pointX, pointY, x[k], y[k], x[l], y[l] ) )
535  {
536  isValid = false;
537  break;
538  }
539  }
540 
541 
542  if ( isValid && b < c )
543  {
544  c = b;
545  fps = ps;
546  fpe = pe;
547  fptx = ptx;
548  fpty = pty;
549  }
550  } // for point which are not in hole
551 
552  // we will cut the shapeu in two new shapes, one from [retainedPoint] to [newPoint] and one form [newPoint] to [retainedPoint]
553  int imin = retainedPt;
554  int imax = ((( fps < retainedPt && fpe < retainedPt ) || ( fps > retainedPt && fpe > retainedPt ) ) ? qMin( fps, fpe ) : qMax( fps, fpe ) );
555 
556  int nbPtSh1, nbPtSh2; // how many points in new shapes ?
557  if ( imax > imin )
558  nbPtSh1 = imax - imin + 1 + ( fpe != fps );
559  else
560  nbPtSh1 = imax + nbp - imin + 1 + ( fpe != fps );
561 
562  if (( imax == fps ? fpe : fps ) < imin )
563  nbPtSh2 = imin - ( imax == fps ? fpe : fps ) + 1 + ( fpe != fps );
564  else
565  nbPtSh2 = imin + nbp - ( imax == fps ? fpe : fps ) + 1 + ( fpe != fps );
566 
567  if ( retainedPt == -1 || fps == -1 || fpe == -1 )
568  {
569  if ( shape->parent )
570  delete shape;
571  }
572  // check for useless spliting
573  else if ( imax == imin || nbPtSh1 <= 2 || nbPtSh2 <= 2 || nbPtSh1 == nbp || nbPtSh2 == nbp )
574  {
575  shapes_final.append( shape );
576  }
577  else
578  {
579 
580  PointSet *newShape = shape->extractShape( nbPtSh1, imin, imax, fps, fpe, fptx, fpty );
581 
582  if ( shape->parent )
583  newShape->parent = shape->parent;
584  else
585  newShape->parent = shape;
586 
587  shapes_toProcess.append( newShape );
588 
589  if ( imax == fps )
590  imax = fpe;
591  else
592  imax = fps;
593 
594  newShape = shape->extractShape( nbPtSh2, imax, imin, fps, fpe, fptx, fpty );
595 
596  if ( shape->parent )
597  newShape->parent = shape->parent;
598  else
599  newShape->parent = shape;
600 
601  shapes_toProcess.append( newShape );
602 
603  if ( shape->parent )
604  delete shape;
605  }
606  }
607  else
608  {
609  shapes_final.append( shape );
610  }
611  delete[] pts;
612  }
613 }
614 
616 {
617  int i;
618  int j;
619 
620  double bbox[4]; // xmin, ymin, xmax, ymax
621 
622  double alpha;
623  int alpha_d;
624 
625  double alpha_seg;
626 
627  double dref;
628  double d1, d2;
629 
630  double bb[16]; // {ax, ay, bx, by, cx, cy, dx, dy, ex, ey, fx, fy, gx, gy, hx, hy}}
631 
632  double cp;
633  double best_cp;
634  double distNearestPoint;
635 
636  double area;
637  double width;
638  double length;
639 
640  double best_area = DBL_MAX;
641  double best_alpha = -1;
642  double best_bb[16];
643  double best_length = 0;
644  double best_width = 0;
645 
646 
647  bbox[0] = DBL_MAX;
648  bbox[1] = DBL_MAX;
649  bbox[2] = - DBL_MAX;
650  bbox[3] = - DBL_MAX;
651 
652  for ( i = 0; i < cHullSize; i++ )
653  {
654  if ( x[cHull[i]] < bbox[0] )
655  bbox[0] = x[cHull[i]];
656 
657  if ( x[cHull[i]] > bbox[2] )
658  bbox[2] = x[cHull[i]];
659 
660  if ( y[cHull[i]] < bbox[1] )
661  bbox[1] = y[cHull[i]];
662 
663  if ( y[cHull[i]] > bbox[3] )
664  bbox[3] = y[cHull[i]];
665  }
666 
667 
668  dref = bbox[2] - bbox[0];
669 
670  for ( alpha_d = 0; alpha_d < 90; alpha_d++ )
671  {
672  alpha = alpha_d * M_PI / 180.0;
673  d1 = cos( alpha ) * dref;
674  d2 = sin( alpha ) * dref;
675 
676  bb[0] = bbox[0];
677  bb[1] = bbox[3]; // ax, ay
678 
679  bb[4] = bbox[0];
680  bb[5] = bbox[1]; // cx, cy
681 
682  bb[8] = bbox[2];
683  bb[9] = bbox[1]; // ex, ey
684 
685  bb[12] = bbox[2];
686  bb[13] = bbox[3]; // gx, gy
687 
688 
689  bb[2] = bb[0] + d1;
690  bb[3] = bb[1] + d2; // bx, by
691  bb[6] = bb[4] - d2;
692  bb[7] = bb[5] + d1; // dx, dy
693  bb[10] = bb[8] - d1;
694  bb[11] = bb[9] - d2; // fx, fy
695  bb[14] = bb[12] + d2;
696  bb[15] = bb[13] - d1; // hx, hy
697 
698  // adjust all points
699  for ( i = 0; i < 16; i += 4 )
700  {
701 
702  alpha_seg = (( i / 4 > 0 ? ( i / 4 ) - 1 : 3 ) ) * M_PI / 2 + alpha;
703 
704  best_cp = DBL_MAX;
705  for ( j = 0; j < nbPoints; j++ )
706  {
707  cp = GeomFunction::cross_product( bb[i+2], bb[i+3], bb[i], bb[i+1], x[cHull[j]], y[cHull[j]] );
708  if ( cp < best_cp )
709  {
710  best_cp = cp;
711  }
712  }
713 
714  distNearestPoint = best_cp / dref;
715 
716  d1 = cos( alpha_seg ) * distNearestPoint;
717  d2 = sin( alpha_seg ) * distNearestPoint;
718 
719  bb[i] += d1; // x
720  bb[i+1] += d2; // y
721  bb[i+2] += d1; // x
722  bb[i+3] += d2; // y
723  }
724 
725  // compute and compare AREA
726  width = GeomFunction::cross_product( bb[6], bb[7], bb[4], bb[5], bb[12], bb[13] ) / dref;
727  length = GeomFunction::cross_product( bb[2], bb[3], bb[0], bb[1], bb[8], bb[9] ) / dref;
728 
729  area = width * length;
730 
731  if ( area < 0 )
732  area *= -1;
733 
734 
735  if ( best_area - area > EPSILON )
736  {
737  best_area = area;
738  best_length = length;
739  best_width = width;
740  best_alpha = alpha;
741  memcpy( best_bb, bb, sizeof( double ) *16 );
742  }
743  }
744 
745  // best bbox is defined
746 
747  CHullBox * finalBb = new CHullBox();
748 
749  for ( i = 0; i < 16; i = i + 4 )
750  {
751  GeomFunction::computeLineIntersection( best_bb[i], best_bb[i+1], best_bb[i+2], best_bb[i+3],
752  best_bb[( i+4 ) %16], best_bb[( i+5 ) %16], best_bb[( i+6 ) %16], best_bb[( i+7 ) %16],
753  &finalBb->x[int ( i/4 )], &finalBb->y[int ( i/4 )] );
754  }
755 
756  finalBb->alpha = best_alpha;
757  finalBb->width = best_width;
758  finalBb->length = best_length;
759 
760  return finalBb;
761 }
762 
763 double PointSet::minDistanceToPoint( double px, double py, double *rx, double *ry ) const
764 {
765  if ( !mGeos )
766  createGeosGeom();
767 
768  if ( !mGeos )
769  return 0;
770 
771  GEOSContextHandle_t geosctxt = geosContext();
772  try
773  {
774  GEOSCoordSequence *coord = GEOSCoordSeq_create_r( geosctxt, 1, 2 );
775  GEOSCoordSeq_setX_r( geosctxt, coord, 0, px );
776  GEOSCoordSeq_setY_r( geosctxt, coord, 0, py );
777  GEOSGeometry* geosPt = GEOSGeom_createPoint_r( geosctxt, coord );
778 
779  int type = GEOSGeomTypeId_r( geosctxt, mGeos );
780  const GEOSGeometry* extRing = nullptr;
781  if ( type != GEOS_POLYGON )
782  {
783  extRing = mGeos;
784  }
785  else
786  {
787  //for polygons, we want distance to exterior ring (not an interior point)
788  extRing = GEOSGetExteriorRing_r( geosctxt, mGeos );
789  }
790  GEOSCoordSequence *nearestCoord = GEOSNearestPoints_r( geosctxt, extRing, geosPt );
791  double nx;
792  double ny;
793  ( void )GEOSCoordSeq_getX_r( geosctxt, nearestCoord, 0, &nx );
794  ( void )GEOSCoordSeq_getY_r( geosctxt, nearestCoord, 0, &ny );
795  GEOSCoordSeq_destroy_r( geosctxt, nearestCoord );
796  GEOSGeom_destroy_r( geosctxt, geosPt );
797 
798  if ( rx )
799  *rx = nx;
800  if ( ry )
801  *ry = ny;
802 
803  return GeomFunction::dist_euc2d_sq( px, py, nx, ny );
804  }
805  catch ( GEOSException &e )
806  {
807  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
808  return 0;
809  }
810 }
811 
812 void PointSet::getCentroid( double &px, double &py, bool forceInside ) const
813 {
814  if ( !mGeos )
815  createGeosGeom();
816 
817  if ( !mGeos )
818  return;
819 
820  try
821  {
822  GEOSContextHandle_t geosctxt = geosContext();
823  GEOSGeometry *centroidGeom = GEOSGetCentroid_r( geosctxt, mGeos );
824  if ( centroidGeom )
825  {
826  const GEOSCoordSequence *coordSeq = GEOSGeom_getCoordSeq_r( geosctxt, centroidGeom );
827  GEOSCoordSeq_getX_r( geosctxt, coordSeq, 0, &px );
828  GEOSCoordSeq_getY_r( geosctxt, coordSeq, 0, &py );
829  }
830 
831  // check if centroid inside in polygon
832  if ( forceInside && !containsPoint( px, py ) )
833  {
834  GEOSGeometry *pointGeom = GEOSPointOnSurface_r( geosctxt, mGeos );
835 
836  if ( pointGeom )
837  {
838  const GEOSCoordSequence *coordSeq = GEOSGeom_getCoordSeq_r( geosctxt, pointGeom );
839  GEOSCoordSeq_getX_r( geosctxt, coordSeq, 0, &px );
840  GEOSCoordSeq_getY_r( geosctxt, coordSeq, 0, &py );
841  GEOSGeom_destroy_r( geosctxt, pointGeom );
842  }
843  }
844 
845  GEOSGeom_destroy_r( geosctxt, centroidGeom );
846  }
847  catch ( GEOSException &e )
848  {
849  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
850  return;
851  }
852 }
853 
854 void PointSet::getPointByDistance( double *d, double *ad, double dl, double *px, double *py )
855 {
856  int i;
857  double dx, dy, di;
858  double distr;
859 
860  i = 0;
861  if ( dl >= 0 )
862  {
863  while ( i < nbPoints && ad[i] <= dl ) i++;
864  i--;
865  }
866 
867  if ( i < nbPoints - 1 )
868  {
869  if ( dl < 0 )
870  {
871  dx = x[nbPoints-1] - x[0];
872  dy = y[nbPoints-1] - y[0];
873  di = sqrt( dx * dx + dy * dy );
874  }
875  else
876  {
877  dx = x[i+1] - x[i];
878  dy = y[i+1] - y[i];
879  di = d[i];
880  }
881 
882  distr = dl - ad[i];
883  *px = x[i] + dx * distr / di;
884  *py = y[i] + dy * distr / di;
885  }
886  else // just select last point...
887  {
888  *px = x[i];
889  *py = y[i];
890  }
891 }
892 
893 const GEOSGeometry *PointSet::geos() const
894 {
895  if ( !mGeos )
896  createGeosGeom();
897 
898  return mGeos;
899 }
900 
901 double PointSet::length() const
902 {
903  if ( !mGeos )
904  createGeosGeom();
905 
906  if ( !mGeos )
907  return -1;
908 
909  GEOSContextHandle_t geosctxt = geosContext();
910 
911  try
912  {
913  double len = 0;
914  ( void )GEOSLength_r( geosctxt, mGeos, &len );
915  return len;
916  }
917  catch ( GEOSException &e )
918  {
919  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
920  return -1;
921  }
922 }
double length
Definition: pointset.h:56
void invalidateGeos()
Definition: pointset.cpp:204
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
const GEOSGeometry * geos() const
Returns the point set&#39;s GEOS geometry.
Definition: pointset.cpp:893
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:150
struct pal::_cHullBox CHullBox
void getCentroid(double &px, double &py, bool forceInside=false) const
Definition: pointset.cpp:812
QString tr(const char *sourceText, const char *disambiguation, int n)
bool qgsDoubleNear(double a, double b, double epsilon=4 *DBL_EPSILON)
Definition: qgis.h:285
static double cross_product(double x1, double y1, double x2, double y2, double x3, double y3)
Definition: geomfunction.h:54
double width
Definition: pointset.h:55
const GEOSPreparedGeometry * preparedGeom() const
Definition: pointset.cpp:192
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:854
PointSet * parent
Definition: pointset.h:161
virtual ~PointSet()
Definition: pointset.cpp:215
double * x
Definition: pointset.h:152
double ymax
Definition: pointset.h:175
double xmin
Definition: pointset.h:172
PointSet * holeOf
Definition: pointset.h:160
double ymin
Definition: pointset.h:174
int * cHull
Definition: pointset.h:155
#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)
double length() const
Returns length of line geometry.
Definition: pointset.cpp:901
static double dist_euc2d(double x1, double y1, double x2, double y2)
Definition: geomfunction.h:59
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
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
bool containsPoint(double x, double y) const
Tests whether point set contains a specified point.
Definition: pointset.cpp:268
double * y
Definition: pointset.h:153
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:338
CHullBox * compute_chull_bbox()
Definition: pointset.cpp:615
double x[4]
Definition: pointset.h:50
double y[4]
Definition: pointset.h:51
#define EPSILON
Definition: util.h:77
GEOSGeometry * mGeos
Definition: pointset.h:148
double alpha
Definition: pointset.h:53
static double dist_euc2d_sq(double x1, double y1, double x2, double y2)
Definition: geomfunction.h:64
double xmax
Definition: pointset.h:173
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:763
bool mOwnsGeom
Definition: pointset.h:149
void append(const T &value)