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