QGIS API Documentation  2.11.0-Master
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
geomfunction.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 "geomfunction.h"
31 #include "feature.h"
32 #include "util.h"
33 
34 namespace pal
35 {
36 
37  void heapsort( int *sid, int *id, const double* const x, int N )
38  {
39  unsigned int n = N, i = n / 2, parent, child;
40  int tx;
41  for ( ;; )
42  {
43  if ( i > 0 )
44  {
45  i--;
46  tx = sid[i];
47  }
48  else
49  {
50  n--;
51  if ( n == 0 ) return;
52  tx = sid[n];
53  sid[n] = sid[0];
54  }
55  parent = i;
56  child = i * 2 + 1;
57  while ( child < n )
58  {
59  if ( child + 1 < n && x[id[sid[child + 1]]] > x[id[sid[child]]] )
60  {
61  child++;
62  }
63  if ( x[id[sid[child]]] > x[id[tx]] )
64  {
65  sid[parent] = sid[child];
66  parent = child;
67  child = parent * 2 + 1;
68  }
69  else
70  {
71  break;
72  }
73  }
74  sid[parent] = tx;
75  }
76  }
77 
78 
79  void heapsort2( int *x, double* heap, int N )
80  {
81  unsigned int n = N, i = n / 2, parent, child;
82  double t;
83  int tx;
84  for ( ;; )
85  {
86  if ( i > 0 )
87  {
88  i--;
89  t = heap[i];
90  tx = x[i];
91  }
92  else
93  {
94  n--;
95  if ( n == 0 ) return;
96  t = heap[n];
97  tx = x[n];
98  heap[n] = heap[0];
99  x[n] = x[0];
100  }
101  parent = i;
102  child = i * 2 + 1;
103  while ( child < n )
104  {
105  if ( child + 1 < n && heap[child + 1] > heap[child] )
106  {
107  child++;
108  }
109  if ( heap[child] > t )
110  {
111  heap[parent] = heap[child];
112  x[parent] = x[child];
113  parent = child;
114  child = parent * 2 + 1;
115  }
116  else
117  {
118  break;
119  }
120  }
121  heap[parent] = t;
122  x[parent] = tx;
123  }
124  }
125 
126 
127 
128 
129  /*
130  * \brief return true if the two seg intersect
131  */
132 
133  bool isSegIntersects( double x1, double y1, double x2, double y2, // 1st segment
134  double x3, double y3, double x4, double y4 ) // 2nd segment
135  {
136  /*
137  std::cout << "SegInrersect ? " << std::endl;
138  std::cout << " cp1 : " << cross_product (x1, y1, x2, y2, x3, y3) << std::endl;
139  std::cout << " cp2 : " << cross_product (x1, y1, x2, y2, x4, y4) << std::endl;
140  std::cout << " cp3 : " << cross_product (x3, y3, x4, y4, x1, y1) << std::endl;
141  std::cout << " cp4 : " << cross_product (x3, y3, x4, y4, x2, y2) << std::endl;
142  */
143  return ( cross_product( x1, y1, x2, y2, x3, y3 ) * cross_product( x1, y1, x2, y2, x4, y4 ) < 0
144  && cross_product( x3, y3, x4, y4, x1, y1 ) * cross_product( x3, y3, x4, y4, x2, y2 ) < 0 );
145  }
146 
147 
148 
149 
150  /*
151  */
152 
153  bool computeSegIntersectionExt( double x1, double y1, double x2, double y2, double xs1, double ys1, // 1st (segment)
154  double x3, double y3, double x4, double y4, double xs2, double ys2, // 2nd segment
155  double *x, double *y )
156  {
157  double cp1, cp2, cp3, cp4;
158  cp1 = cross_product( x1, y1, x2, y2, x3, y3 );
159  cp2 = cross_product( x1, y1, x2, y2, x4, y4 );
160  cp3 = cross_product( x3, y3, x4, y4, x1, y1 );
161  cp4 = cross_product( x3, y3, x4, y4, x2, y2 );
162 
163 
164  if ( cp1 == 0 && cp2 == 0 && cp3 == 0 && cp4 == 0 )
165  {
166 #ifdef _DEBUG_FULL_
167  std::cout << "coolineaire..." << std::endl;
168 #endif
169  return false;
170  }
171 
172  // 1 ter
173  if ( cp1 == 0 && cp3 == 0 )
174  {
175 #ifdef _DEBUG_FULL_
176  std::cout << "cp1 = cp3 = 0 => ignoring..." << std::endl;
177 #endif
178  return false;
179  }
180 
181  // 1 bis
182  if ( cp1 == 0 && cp4 == 0 )
183  {
184 #ifdef _DEBUG_FULL_
185  std::cout << "cp1 = cp4 = 0 => ignoring..." << std::endl;
186 #endif
187  return false;
188  }
189 
190  // 1 bis
191  if ( cp2 == 0 && cp3 == 0 )
192  {
193 #ifdef _DEBUG_FULL_
194  std::cout << "cp2 = cp3 = 0 => ignoring..." << std::endl;
195 #endif
196  return false;
197  }
198 
199  // 2bis and 3bis
200  if ( cp1 == 0 || cp3 == 0 )
201  {
202 #ifdef _DEBUG_FULL_
203  std::cout << "skip..." << std::endl;
204 #endif
205  return false;
206  }
207 
208  // case 3
209  if ( cp4 == 0 && cp1 * cp1 < 0 )
210  {
211  if ( cross_product( x3, y3, x4, y4, xs1, ys1 ) * cp3 < 0 )
212  {
213  *x = x2;
214  *y = y2;
215  return true;
216  }
217  else
218  return false;
219  }
220 
221  // case 2
222  if ( cp2 == 0 && cp3 * cp4 < 0 )
223  {
224  if ( cross_product( x1, y1, x2, y2, xs2, ys2 ) * cp1 < 0 )
225  {
226  *x = x4;
227  *y = y4;
228  return true;
229  }
230  else
231  return false;
232  }
233 
234  // case 1
235  if ( cp2 == 0 && cp4 == 0 )
236  {
237  double distance[4];
238  double cx, cy;
239  double dx, dy;
240  double nx[4], ny[4];
241  double toDist;
242  double ratio;
243  int i;
244 
245  cx = x2;
246  cy = y2;
247 
248  nx[0] = x1;
249  ny[0] = y1;
250 
251  nx[1] = xs1;
252  ny[1] = ys1;
253 
254  nx[2] = x3;
255  ny[2] = y3;
256 
257  nx[3] = xs2;
258  ny[3] = ys2;
259 
260  distance[0] = dist_euc2d( cx, cy, x1, y1 ); // i
261  toDist = distance[0];
262 
263  distance[1] = dist_euc2d( cx, cy, xs1, ys1 );// j2
264  toDist = max( toDist, distance[1] );
265 
266  distance[2] = dist_euc2d( cx, cy, x3, y3 );// k
267  toDist = max( toDist, distance[2] );
268 
269  distance[3] = dist_euc2d( cx, cy, xs2, ys2 ); // l2
270  toDist = max( toDist, distance[3] );
271 
272  for ( i = 0; i < 4; i++ )
273  {
274  dx = nx[i] - cx;
275  dy = ny[i] - cy;
276 
277  ratio = toDist / distance[i];
278 
279  nx[i] = cx + dx * ratio;
280  ny[i] = cy + dy * ratio;
281  }
282 
283  bool return_val = computeSegIntersection( nx[0], ny[0], nx[1], ny[1], nx[2], ny[2], nx[3], ny[3], x, y );
284 
285  return return_val;
286  }
287 
288  if ( cp1 * cp2 <= 0
289  && cp3 *cp4 <= 0 )
290  {
291  return computeLineIntersection( x1, y1, x2, y2, x3, y3, x4, y4, x, y );
292  }
293 
294  return false;
295  }
296 
297 
298 
299 
300  /*
301  * \brief Intersection bw a line and a segment
302  * \return true if the point exist false otherwise
303  */
304  bool computeLineSegIntersection( double x1, double y1, double x2, double y2, // 1st line
305  double x3, double y3, double x4, double y4, // 2nd segment
306  double *x, double *y )
307  {
308  double cp1, cp2;
309  cp1 = cross_product( x1, y1, x2, y2, x3, y3 );
310  cp2 = cross_product( x1, y1, x2, y2, x4, y4 );
311 
312  if ( cp1 * cp2 <= 0 )
313  return computeLineIntersection( x1, y1, x2, y2, x3, y3, x4, y4, x, y );
314 
315  return false;
316  }
317 
318 
319 
320  /*
321  * \brief compute the point wherre two segment intersects
322  * \return true if the point exist false otherwise
323  */
324 
325  bool computeSegIntersection( double x1, double y1, double x2, double y2, // 1st (segment)
326  double x3, double y3, double x4, double y4, // 2nd segment
327  double *x, double *y )
328  {
329  double cp1, cp2, cp3, cp4;
330  cp1 = cross_product( x1, y1, x2, y2, x3, y3 );
331  cp2 = cross_product( x1, y1, x2, y2, x4, y4 );
332  cp3 = cross_product( x3, y3, x4, y4, x1, y1 );
333  cp4 = cross_product( x3, y3, x4, y4, x2, y2 );
334 
335  if ( cp1 * cp2 <= 0
336  && cp3 *cp4 <= 0 )
337  return computeLineIntersection( x1, y1, x2, y2, x3, y3, x4, y4, x, y );
338 
339  return false;
340  }
341 
342  /*
343  * \brief compute the point wherre two lines intersects
344  * \return true if the ok false if line are parallel
345  */
346  bool computeLineIntersection( double x1, double y1, double x2, double y2, // 1st line (segment)
347  double x3, double y3, double x4, double y4, // 2nd line segment
348  double *x, double *y )
349  {
350 
351  double a1, a2, b1, b2, c1, c2;
352  double denom;
353 
354  a1 = y2 - y1;
355  b1 = x1 - x2;
356  c1 = x2 * y1 - x1 * y2;
357 
358  a2 = y4 - y3;
359  b2 = x3 - x4;
360  c2 = x4 * y3 - x3 * y4;
361 
362  if (( denom = a1 * b2 - a2 * b1 ) == 0 )
363  {
364  return false;
365  }
366  else
367  {
368  *x = ( b1 * c2 - b2 * c1 ) / denom;
369  *y = ( a2 * c1 - a1 * c2 ) / denom;
370  }
371 
372  return true;
373 
374  }
375 
376 
377 
378  /*
379  * \brief Compute the convex hull in O(n·log(n))
380  * \param id set of point (i.e. point no 0 is (x,y) = x[id[0]],y[id[0]])
381  * \param x x coordinates
382  * \param y y coordinates
383  * \param n Size of subset (vector id)
384  * \param cHull returns the point id (id of id's vector...) whom are parts of the convex hull
385  * \return convexHull's size
386  */
387  int convexHullId( int *id, const double* const x, const double* const y, int n, int *&cHull )
388  {
389  int i;
390 
391  cHull = new int[n];
392  for ( i = 0; i < n; i++ )
393  {
394  cHull[i] = i;
395  }
396 
397 
398  if ( n <= 3 ) return n;
399 
400  int* stack = new int[n];
401  double* tan = new double [n];
402  int ref;
403 
404  int second, top;
405  double result;
406 
407  // find the lowest y value
408  heapsort( cHull, id, y, n );
409 
410  // find the lowest x value from the lowest y
411  ref = 1;
412  while ( ref < n && vabs( y[id[cHull[ref]]] - y[id[cHull[0]]] ) < EPSILON ) ref++;
413 
414  heapsort( cHull, id, x, ref );
415 
416  // the first point is now for sure in the hull as well as the ref one
417  for ( i = ref; i < n; i++ )
418  {
419  if ( vabs( y[id[cHull[i]]] - y[id[cHull[0]]] ) < EPSILON )
420  tan[i] = FLT_MAX;
421  else
422  tan[i] = ( x[id[cHull[0]]] - x[id[cHull[i]]] ) / ( y[id[cHull[i]]] - y[id[cHull[0]]] );
423  }
424 
425  if ( ref < n )
426  heapsort2( cHull + ref, tan + ref, n - ref );
427 
428  // the second point is in too
429  stack[0] = cHull[0];
430  if ( ref == 1 )
431  {
432  stack[1] = cHull[1];
433  ref++;
434  }
435  else
436  stack[1] = cHull[ref-1];
437 
438 
439  top = 1;
440  second = 0;
441 
442  for ( i = ref; i < n; i++ )
443  {
444  result = cross_product( x[id[stack[second]]], y[id[stack[second]]],
445  x[id[stack[top]]], y[id[stack[top]]], x[id[cHull[i]]], y[id[cHull[i]]] );
446  // Coolineaire !! garder le plus éloigné
447  if ( vabs( result ) < EPSILON )
448  {
449  if ( dist_euc2d_sq( x[id[stack[second]]], y[id[stack[second]]], x[id[cHull[i]]], y[id[cHull[i]]] )
450  > dist_euc2d_sq( x[id[stack[second]]], y[id[stack[second]]], x[id[stack[top]]], y[id[stack[top]]] ) )
451  {
452  stack[top] = cHull[i];
453  }
454  }
455  else if ( result > 0 ) //convexe
456  {
457  second++;
458  top++;
459  stack[top] = cHull[i];
460  }
461  else
462  {
463  while ( result < 0 && top > 1 )
464  {
465  second--;
466  top--;
467  result = cross_product( x[id[stack[second]]],
468  y[id[stack[second]]], x[id[stack[top]]],
469  y[id[stack[top]]], x[id[cHull[i]]], y[id[cHull[i]]] );
470  }
471  second++;
472  top++;
473  stack[top] = cHull[i];
474  }
475  }
476 
477  for ( i = 0; i <= top; i++ )
478  {
479  cHull[i] = stack[i];
480  }
481 
482  delete[] stack;
483  delete[] tan;
484 
485  return top + 1;
486  }
487 
488 // reorder points to have cross prod ((x,y)[i], (x,y)[i+1), point) > 0 when point is outside
489  int reorderPolygon( int nbPoints, double *x, double *y )
490  {
491  int inc = 0;
492  int *cHull;
493  int cHullSize;
494  int i;
495 
496  int *pts = new int[nbPoints];
497  for ( i = 0; i < nbPoints; i++ )
498  pts[i] = i;
499 
500 
501  cHullSize = convexHullId( pts, x, y, nbPoints, cHull );
502 
503  if ( pts[cHull[0]] < pts[cHull[1]] && pts[cHull[1]] < pts[cHull[2]] )
504  inc = 1;
505  else if ( pts[cHull[0]] > pts[cHull[1]] && pts[cHull[1]] > pts[cHull[2]] )
506  inc = -1;
507  else if ( pts[cHull[0]] > pts[cHull[1]] && pts[cHull[1]] < pts[cHull[2]] && pts[cHull[2]] < pts[cHull[0]] )
508  inc = 1;
509  else if ( pts[cHull[0]] > pts[cHull[1]] && pts[cHull[1]] < pts[cHull[2]] && pts[cHull[2]] > pts[cHull[0]] )
510  inc = -1;
511  else if ( pts[cHull[0]] < pts[cHull[1]] && pts[cHull[1]] > pts[cHull[2]] && pts[cHull[2]] > pts[cHull[0]] )
512  inc = -1;
513  else if ( pts[cHull[0]] < pts[cHull[1]] && pts[cHull[1]] > pts[cHull[2]] && pts[cHull[2]] < pts[cHull[0]] )
514  inc = 1;
515  else
516  {
517  std::cout << "Warning wrong cHull -> geometry: " << nbPoints << std::endl;
518  for ( i = 0; i < nbPoints; i++ )
519  {
520  std::cout << x[i] << ";" << y[i] << std::endl;
521  }
522  std::cout << "hull : " << cHullSize << std::endl;
523  for ( i = 0; i < cHullSize; i++ )
524  {
525  std::cout << pts[cHull[i]] << " ";
526  }
527  std::cout << std::endl;
528  delete[] cHull;
529  delete[] pts;
530  return -1;
531  }
532 
533  if ( inc == -1 ) // re-order points
534  {
535  double tmp;
536  int j;
537  for ( i = 0, j = nbPoints - 1; i <= j; i++, j-- )
538  {
539  tmp = x[i];
540  x[i] = x[j];
541  x[j] = tmp;
542 
543  tmp = y[i];
544  y[i] = y[j];
545  y[j] = tmp;
546  }
547  }
548 
549 
550  delete[] cHull;
551  delete[] pts;
552 
553  return 0;
554 
555  }
556 
557 
558  bool isPointInPolygon( int npol, double *xp, double *yp, double x, double y )
559  {
560  // code from Randolph Franklin (found at http://local.wasp.uwa.edu.au/~pbourke/geometry/insidepoly/)
561  int i, j;
562  bool c = false;
563 
564  for ( i = 0, j = npol - 1; i < npol; j = i++ )
565  {
566  if (((( yp[i] <= y ) && ( y < yp[j] ) ) ||
567  (( yp[j] <= y ) && ( y < yp[i] ) ) )
568  && ( x < ( xp[j] - xp[i] ) *( y - yp[i] ) / ( yp[j] - yp[i] ) + xp[i] ) )
569  {
570  c = !c;
571  }
572  }
573  return c;
574  }
575 
576  void findLineCircleIntersection( double cx, double cy, double radius,
577  double x1, double y1, double x2, double y2,
578  double& xRes, double& yRes )
579  {
580  double dx = x2 - x1;
581  double dy = y2 - y1;
582 
583  double A = dx * dx + dy * dy;
584  double B = 2 * ( dx * ( x1 - cx ) + dy * ( y1 - cy ) );
585  double C = ( x1 - cx ) * ( x1 - cx ) + ( y1 - cy ) * ( y1 - cy ) - radius * radius;
586 
587  double det = B * B - 4 * A * C;
588  if ( A <= 0.0000001 || det < 0 )
589  // Should never happen, No real solutions.
590  return;
591 
592  if ( det == 0 )
593  {
594  // Could potentially happen.... One solution.
595  double t = -B / ( 2 * A );
596  xRes = x1 + t * dx;
597  yRes = y1 + t * dy;
598  }
599  else
600  {
601  // Two solutions.
602  // Always use the 1st one
603  // We only really have one solution here, as we know the line segment will start in the circle and end outside
604  double t = ( -B + sqrt( det ) ) / ( 2 * A );
605  xRes = x1 + t * dx;
606  yRes = y1 + t * dy;
607  }
608  }
609 
610 
611 } // end namespace
int reorderPolygon(int nbPoints, double *x, double *y)
bool computeLineIntersection(double x1, double y1, double x2, double y2, double x3, double y3, double x4, double y4, double *x, double *y)
void heapsort(int *sid, int *id, const double *const x, int N)
double dist_euc2d_sq(double x1, double y1, double x2, double y2)
Definition: geomfunction.h:57
bool computeSegIntersection(double x1, double y1, double x2, double y2, double x3, double y3, double x4, double y4, double *x, double *y)
bool isPointInPolygon(int npol, double *xp, double *yp, double x, double y)
bool computeSegIntersectionExt(double x1, double y1, double x2, double y2, double xs1, double ys1, double x3, double y3, double x4, double y4, double xs2, double ys2, double *x, double *y)
double dist_euc2d(double x1, double y1, double x2, double y2)
Definition: geomfunction.h:52
void heapsort2(int *x, double *heap, int N)
double cross_product(double x1, double y1, double x2, double y2, double x3, double y3)
Definition: geomfunction.h:47
bool computeLineSegIntersection(double x1, double y1, double x2, double y2, double x3, double y3, double x4, double y4, double *x, double *y)
int convexHullId(int *id, const double *const x, const double *const y, int n, int *&cHull)
#define EPSILON
Definition: util.h:79
bool isSegIntersects(double x1, double y1, double x2, double y2, double x3, double y3, double x4, double y4)
void findLineCircleIntersection(double cx, double cy, double radius, double x1, double y1, double x2, double y2, double &xRes, double &yRes)
int max(int a, int b)
Definition: util.h:82
double vabs(double x)
Definition: util.h:94