QGIS API Documentation
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 using namespace pal;
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 bool GeomFunction::isSegIntersects( double x1, double y1, double x2, double y2, // 1st segment
127  double x3, double y3, double x4, double y4 ) // 2nd segment
128 {
129  return ( cross_product( x1, y1, x2, y2, x3, y3 ) * cross_product( x1, y1, x2, y2, x4, y4 ) < 0
130  && cross_product( x3, y3, x4, y4, x1, y1 ) * cross_product( x3, y3, x4, y4, x2, y2 ) < 0 );
131 }
132 
133 bool GeomFunction::computeLineIntersection( double x1, double y1, double x2, double y2, // 1st line (segment)
134  double x3, double y3, double x4, double y4, // 2nd line segment
135  double *x, double *y )
136 {
137 
138  double a1, a2, b1, b2, c1, c2;
139  double denom;
140 
141  a1 = y2 - y1;
142  b1 = x1 - x2;
143  c1 = x2 * y1 - x1 * y2;
144 
145  a2 = y4 - y3;
146  b2 = x3 - x4;
147  c2 = x4 * y3 - x3 * y4;
148 
149  denom = a1 * b2 - a2 * b1;
150  if ( qgsDoubleNear( denom, 0.0 ) )
151  {
152  return false;
153  }
154  else
155  {
156  *x = ( b1 * c2 - b2 * c1 ) / denom;
157  *y = ( a2 * c1 - a1 * c2 ) / denom;
158  }
159 
160  return true;
161 }
162 
163 int GeomFunction::convexHullId( int *id, const double* const x, const double* const y, int n, int *&cHull )
164 {
165  int i;
166 
167  cHull = new int[n];
168  for ( i = 0; i < n; i++ )
169  {
170  cHull[i] = i;
171  }
172 
173 
174  if ( n <= 3 ) return n;
175 
176  int* stack = new int[n];
177  double* tan = new double [n];
178  int ref;
179 
180  int second, top;
181  double result;
182 
183  // find the lowest y value
184  heapsort( cHull, id, y, n );
185 
186  // find the lowest x value from the lowest y
187  ref = 1;
188  while ( ref < n && qgsDoubleNear( y[id[cHull[ref]]], y[id[cHull[0]]] ) ) ref++;
189 
190  heapsort( cHull, id, x, ref );
191 
192  // the first point is now for sure in the hull as well as the ref one
193  for ( i = ref; i < n; i++ )
194  {
195  if ( qgsDoubleNear( y[id[cHull[i]]], y[id[cHull[0]]] ) )
196  tan[i] = FLT_MAX;
197  else
198  tan[i] = ( x[id[cHull[0]]] - x[id[cHull[i]]] ) / ( y[id[cHull[i]]] - y[id[cHull[0]]] );
199  }
200 
201  if ( ref < n )
202  heapsort2( cHull + ref, tan + ref, n - ref );
203 
204  // the second point is in too
205  stack[0] = cHull[0];
206  if ( ref == 1 )
207  {
208  stack[1] = cHull[1];
209  ref++;
210  }
211  else
212  stack[1] = cHull[ref-1];
213 
214 
215  top = 1;
216  second = 0;
217 
218  for ( i = ref; i < n; i++ )
219  {
220  result = cross_product( x[id[stack[second]]], y[id[stack[second]]],
221  x[id[stack[top]]], y[id[stack[top]]], x[id[cHull[i]]], y[id[cHull[i]]] );
222  // Coolineaire !! garder le plus éloigné
223  if ( qgsDoubleNear( result, 0.0 ) )
224  {
225  if ( dist_euc2d_sq( x[id[stack[second]]], y[id[stack[second]]], x[id[cHull[i]]], y[id[cHull[i]]] )
226  > dist_euc2d_sq( x[id[stack[second]]], y[id[stack[second]]], x[id[stack[top]]], y[id[stack[top]]] ) )
227  {
228  stack[top] = cHull[i];
229  }
230  }
231  else if ( result > 0 ) //convexe
232  {
233  second++;
234  top++;
235  stack[top] = cHull[i];
236  }
237  else
238  {
239  while ( result < 0 && top > 1 )
240  {
241  second--;
242  top--;
243  result = cross_product( x[id[stack[second]]],
244  y[id[stack[second]]], x[id[stack[top]]],
245  y[id[stack[top]]], x[id[cHull[i]]], y[id[cHull[i]]] );
246  }
247  second++;
248  top++;
249  stack[top] = cHull[i];
250  }
251  }
252 
253  for ( i = 0; i <= top; i++ )
254  {
255  cHull[i] = stack[i];
256  }
257 
258  delete[] stack;
259  delete[] tan;
260 
261  return top + 1;
262 }
263 
264 int GeomFunction::reorderPolygon( int nbPoints, double *x, double *y )
265 {
266  int inc = 0;
267  int *cHull;
268  int i;
269 
270  int *pts = new int[nbPoints];
271  for ( i = 0; i < nbPoints; i++ )
272  pts[i] = i;
273 
274  ( void )convexHullId( pts, x, y, nbPoints, cHull );
275 
276  if ( pts[cHull[0]] < pts[cHull[1]] && pts[cHull[1]] < pts[cHull[2]] )
277  inc = 1;
278  else if ( pts[cHull[0]] > pts[cHull[1]] && pts[cHull[1]] > pts[cHull[2]] )
279  inc = -1;
280  else if ( pts[cHull[0]] > pts[cHull[1]] && pts[cHull[1]] < pts[cHull[2]] && pts[cHull[2]] < pts[cHull[0]] )
281  inc = 1;
282  else if ( pts[cHull[0]] > pts[cHull[1]] && pts[cHull[1]] < pts[cHull[2]] && pts[cHull[2]] > pts[cHull[0]] )
283  inc = -1;
284  else if ( pts[cHull[0]] < pts[cHull[1]] && pts[cHull[1]] > pts[cHull[2]] && pts[cHull[2]] > pts[cHull[0]] )
285  inc = -1;
286  else if ( pts[cHull[0]] < pts[cHull[1]] && pts[cHull[1]] > pts[cHull[2]] && pts[cHull[2]] < pts[cHull[0]] )
287  inc = 1;
288  else
289  {
290  // wrong cHull
291  delete[] cHull;
292  delete[] pts;
293  return -1;
294  }
295 
296  if ( inc == -1 ) // re-order points
297  {
298  double tmp;
299  int j;
300  for ( i = 0, j = nbPoints - 1; i <= j; i++, j-- )
301  {
302  tmp = x[i];
303  x[i] = x[j];
304  x[j] = tmp;
305 
306  tmp = y[i];
307  y[i] = y[j];
308  y[j] = tmp;
309  }
310  }
311 
312  delete[] cHull;
313  delete[] pts;
314 
315  return 0;
316 }
317 
318 void GeomFunction::findLineCircleIntersection( double cx, double cy, double radius,
319  double x1, double y1, double x2, double y2,
320  double& xRes, double& yRes )
321 {
322  double dx = x2 - x1;
323  double dy = y2 - y1;
324 
325  double A = dx * dx + dy * dy;
326  double B = 2 * ( dx * ( x1 - cx ) + dy * ( y1 - cy ) );
327  double C = ( x1 - cx ) * ( x1 - cx ) + ( y1 - cy ) * ( y1 - cy ) - radius * radius;
328 
329  double det = B * B - 4 * A * C;
330  if ( A <= 0.0000001 || det < 0 )
331  // Should never happen, No real solutions.
332  return;
333 
334  if ( qgsDoubleNear( det, 0.0 ) )
335  {
336  // Could potentially happen.... One solution.
337  double t = -B / ( 2 * A );
338  xRes = x1 + t * dx;
339  yRes = y1 + t * dy;
340  }
341  else
342  {
343  // Two solutions.
344  // Always use the 1st one
345  // We only really have one solution here, as we know the line segment will start in the circle and end outside
346  double t = ( -B + sqrt( det ) ) / ( 2 * A );
347  xRes = x1 + t * dx;
348  yRes = y1 + t * dy;
349  }
350 }
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.
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.
static void findLineCircleIntersection(double cx, double cy, double radius, double x1, double y1, double x2, double y2, double &xRes, double &yRes)
bool qgsDoubleNear(double a, double b, double epsilon=4 *DBL_EPSILON)
Compare two doubles (but allow some difference)
Definition: qgis.h:348
static double cross_product(double x1, double y1, double x2, double y2, double x3, double y3)
Definition: geomfunction.h:54
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))
void heapsort(int *sid, int *id, const double *const x, int N)
void heapsort2(int *x, double *heap, int N)
static int reorderPolygon(int nbPoints, double *x, double *y)
Reorder points to have cross prod ((x,y)[i], (x,y)[i+1), point) > 0 when point is outside...
static double dist_euc2d_sq(double x1, double y1, double x2, double y2)
Definition: geomfunction.h:64