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