QGIS API Documentation  2.11.0-Master
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups 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 #include "qgis.h"
34 
35 namespace pal
36 {
37 
38  void heapsort( int *sid, int *id, const double* const x, int N )
39  {
40  unsigned int n = N, i = n / 2, parent, child;
41  int tx;
42  for ( ;; )
43  {
44  if ( i > 0 )
45  {
46  i--;
47  tx = sid[i];
48  }
49  else
50  {
51  n--;
52  if ( n == 0 ) return;
53  tx = sid[n];
54  sid[n] = sid[0];
55  }
56  parent = i;
57  child = i * 2 + 1;
58  while ( child < n )
59  {
60  if ( child + 1 < n && x[id[sid[child + 1]]] > x[id[sid[child]]] )
61  {
62  child++;
63  }
64  if ( x[id[sid[child]]] > x[id[tx]] )
65  {
66  sid[parent] = sid[child];
67  parent = child;
68  child = parent * 2 + 1;
69  }
70  else
71  {
72  break;
73  }
74  }
75  sid[parent] = tx;
76  }
77  }
78 
79 
80  void heapsort2( int *x, double* heap, int N )
81  {
82  unsigned int n = N, i = n / 2, parent, child;
83  double t;
84  int tx;
85  for ( ;; )
86  {
87  if ( i > 0 )
88  {
89  i--;
90  t = heap[i];
91  tx = x[i];
92  }
93  else
94  {
95  n--;
96  if ( n == 0 ) return;
97  t = heap[n];
98  tx = x[n];
99  heap[n] = heap[0];
100  x[n] = x[0];
101  }
102  parent = i;
103  child = i * 2 + 1;
104  while ( child < n )
105  {
106  if ( child + 1 < n && heap[child + 1] > heap[child] )
107  {
108  child++;
109  }
110  if ( heap[child] > t )
111  {
112  heap[parent] = heap[child];
113  x[parent] = x[child];
114  parent = child;
115  child = parent * 2 + 1;
116  }
117  else
118  {
119  break;
120  }
121  }
122  heap[parent] = t;
123  x[parent] = tx;
124  }
125  }
126 
127 
128 
129 
130  /*
131  * \brief return true if the two seg intersect
132  */
133 
134  bool isSegIntersects( double x1, double y1, double x2, double y2, // 1st segment
135  double x3, double y3, double x4, double y4 ) // 2nd segment
136  {
137  /*
138  std::cout << "SegInrersect ? " << std::endl;
139  std::cout << " cp1 : " << cross_product (x1, y1, x2, y2, x3, y3) << std::endl;
140  std::cout << " cp2 : " << cross_product (x1, y1, x2, y2, x4, y4) << std::endl;
141  std::cout << " cp3 : " << cross_product (x3, y3, x4, y4, x1, y1) << std::endl;
142  std::cout << " cp4 : " << cross_product (x3, y3, x4, y4, x2, y2) << std::endl;
143  */
144  return ( cross_product( x1, y1, x2, y2, x3, y3 ) * cross_product( x1, y1, x2, y2, x4, y4 ) < 0
145  && cross_product( x3, y3, x4, y4, x1, y1 ) * cross_product( x3, y3, x4, y4, x2, y2 ) < 0 );
146  }
147 
148  /*
149  * \brief compute the point wherre two lines intersects
150  * \return true if the ok false if line are parallel
151  */
152  bool computeLineIntersection( double x1, double y1, double x2, double y2, // 1st line (segment)
153  double x3, double y3, double x4, double y4, // 2nd line segment
154  double *x, double *y )
155  {
156 
157  double a1, a2, b1, b2, c1, c2;
158  double denom;
159 
160  a1 = y2 - y1;
161  b1 = x1 - x2;
162  c1 = x2 * y1 - x1 * y2;
163 
164  a2 = y4 - y3;
165  b2 = x3 - x4;
166  c2 = x4 * y3 - x3 * y4;
167 
168  if (( denom = a1 * b2 - a2 * b1 ) == 0 )
169  {
170  return false;
171  }
172  else
173  {
174  *x = ( b1 * c2 - b2 * c1 ) / denom;
175  *y = ( a2 * c1 - a1 * c2 ) / denom;
176  }
177 
178  return true;
179 
180  }
181 
182 
183 
184  /*
185  * \brief Compute the convex hull in O(n·log(n))
186  * \param id set of point (i.e. point no 0 is (x,y) = x[id[0]],y[id[0]])
187  * \param x x coordinates
188  * \param y y coordinates
189  * \param n Size of subset (vector id)
190  * \param cHull returns the point id (id of id's vector...) whom are parts of the convex hull
191  * \return convexHull's size
192  */
193  int convexHullId( int *id, const double* const x, const double* const y, int n, int *&cHull )
194  {
195  int i;
196 
197  cHull = new int[n];
198  for ( i = 0; i < n; i++ )
199  {
200  cHull[i] = i;
201  }
202 
203 
204  if ( n <= 3 ) return n;
205 
206  int* stack = new int[n];
207  double* tan = new double [n];
208  int ref;
209 
210  int second, top;
211  double result;
212 
213  // find the lowest y value
214  heapsort( cHull, id, y, n );
215 
216  // find the lowest x value from the lowest y
217  ref = 1;
218  while ( ref < n && qgsDoubleNear( y[id[cHull[ref]]], y[id[cHull[0]]] ) ) ref++;
219 
220  heapsort( cHull, id, x, ref );
221 
222  // the first point is now for sure in the hull as well as the ref one
223  for ( i = ref; i < n; i++ )
224  {
225  if ( qgsDoubleNear( y[id[cHull[i]]], y[id[cHull[0]]] ) )
226  tan[i] = FLT_MAX;
227  else
228  tan[i] = ( x[id[cHull[0]]] - x[id[cHull[i]]] ) / ( y[id[cHull[i]]] - y[id[cHull[0]]] );
229  }
230 
231  if ( ref < n )
232  heapsort2( cHull + ref, tan + ref, n - ref );
233 
234  // the second point is in too
235  stack[0] = cHull[0];
236  if ( ref == 1 )
237  {
238  stack[1] = cHull[1];
239  ref++;
240  }
241  else
242  stack[1] = cHull[ref-1];
243 
244 
245  top = 1;
246  second = 0;
247 
248  for ( i = ref; i < n; i++ )
249  {
250  result = cross_product( x[id[stack[second]]], y[id[stack[second]]],
251  x[id[stack[top]]], y[id[stack[top]]], x[id[cHull[i]]], y[id[cHull[i]]] );
252  // Coolineaire !! garder le plus éloigné
253  if ( qgsDoubleNear( result, 0.0 ) )
254  {
255  if ( dist_euc2d_sq( x[id[stack[second]]], y[id[stack[second]]], x[id[cHull[i]]], y[id[cHull[i]]] )
256  > dist_euc2d_sq( x[id[stack[second]]], y[id[stack[second]]], x[id[stack[top]]], y[id[stack[top]]] ) )
257  {
258  stack[top] = cHull[i];
259  }
260  }
261  else if ( result > 0 ) //convexe
262  {
263  second++;
264  top++;
265  stack[top] = cHull[i];
266  }
267  else
268  {
269  while ( result < 0 && top > 1 )
270  {
271  second--;
272  top--;
273  result = cross_product( x[id[stack[second]]],
274  y[id[stack[second]]], x[id[stack[top]]],
275  y[id[stack[top]]], x[id[cHull[i]]], y[id[cHull[i]]] );
276  }
277  second++;
278  top++;
279  stack[top] = cHull[i];
280  }
281  }
282 
283  for ( i = 0; i <= top; i++ )
284  {
285  cHull[i] = stack[i];
286  }
287 
288  delete[] stack;
289  delete[] tan;
290 
291  return top + 1;
292  }
293 
294 // reorder points to have cross prod ((x,y)[i], (x,y)[i+1), point) > 0 when point is outside
295  int reorderPolygon( int nbPoints, double *x, double *y )
296  {
297  int inc = 0;
298  int *cHull;
299  int cHullSize;
300  int i;
301 
302  int *pts = new int[nbPoints];
303  for ( i = 0; i < nbPoints; i++ )
304  pts[i] = i;
305 
306 
307  cHullSize = convexHullId( pts, x, y, nbPoints, cHull );
308 
309  if ( pts[cHull[0]] < pts[cHull[1]] && pts[cHull[1]] < pts[cHull[2]] )
310  inc = 1;
311  else if ( pts[cHull[0]] > pts[cHull[1]] && pts[cHull[1]] > pts[cHull[2]] )
312  inc = -1;
313  else if ( pts[cHull[0]] > pts[cHull[1]] && pts[cHull[1]] < pts[cHull[2]] && pts[cHull[2]] < pts[cHull[0]] )
314  inc = 1;
315  else if ( pts[cHull[0]] > pts[cHull[1]] && pts[cHull[1]] < pts[cHull[2]] && pts[cHull[2]] > pts[cHull[0]] )
316  inc = -1;
317  else if ( pts[cHull[0]] < pts[cHull[1]] && pts[cHull[1]] > pts[cHull[2]] && pts[cHull[2]] > pts[cHull[0]] )
318  inc = -1;
319  else if ( pts[cHull[0]] < pts[cHull[1]] && pts[cHull[1]] > pts[cHull[2]] && pts[cHull[2]] < pts[cHull[0]] )
320  inc = 1;
321  else
322  {
323  std::cout << "Warning wrong cHull -> geometry: " << nbPoints << std::endl;
324  for ( i = 0; i < nbPoints; i++ )
325  {
326  std::cout << x[i] << ";" << y[i] << std::endl;
327  }
328  std::cout << "hull : " << cHullSize << std::endl;
329  for ( i = 0; i < cHullSize; i++ )
330  {
331  std::cout << pts[cHull[i]] << " ";
332  }
333  std::cout << std::endl;
334  delete[] cHull;
335  delete[] pts;
336  return -1;
337  }
338 
339  if ( inc == -1 ) // re-order points
340  {
341  double tmp;
342  int j;
343  for ( i = 0, j = nbPoints - 1; i <= j; i++, j-- )
344  {
345  tmp = x[i];
346  x[i] = x[j];
347  x[j] = tmp;
348 
349  tmp = y[i];
350  y[i] = y[j];
351  y[j] = tmp;
352  }
353  }
354 
355 
356  delete[] cHull;
357  delete[] pts;
358 
359  return 0;
360 
361  }
362 
363  void findLineCircleIntersection( double cx, double cy, double radius,
364  double x1, double y1, double x2, double y2,
365  double& xRes, double& yRes )
366  {
367  double dx = x2 - x1;
368  double dy = y2 - y1;
369 
370  double A = dx * dx + dy * dy;
371  double B = 2 * ( dx * ( x1 - cx ) + dy * ( y1 - cy ) );
372  double C = ( x1 - cx ) * ( x1 - cx ) + ( y1 - cy ) * ( y1 - cy ) - radius * radius;
373 
374  double det = B * B - 4 * A * C;
375  if ( A <= 0.0000001 || det < 0 )
376  // Should never happen, No real solutions.
377  return;
378 
379  if ( det == 0 )
380  {
381  // Could potentially happen.... One solution.
382  double t = -B / ( 2 * A );
383  xRes = x1 + t * dx;
384  yRes = y1 + t * dy;
385  }
386  else
387  {
388  // Two solutions.
389  // Always use the 1st one
390  // We only really have one solution here, as we know the line segment will start in the circle and end outside
391  double t = ( -B + sqrt( det ) ) / ( 2 * A );
392  xRes = x1 + t * dx;
393  yRes = y1 + t * dy;
394  }
395  }
396 
397 
398 } // 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 qgsDoubleNear(double a, double b, double epsilon=4 *DBL_EPSILON)
Definition: qgis.h:350
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
int convexHullId(int *id, const double *const x, const double *const y, int n, int *&cHull)
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)