QGIS API Documentation  2.17.0-Master (3a3b9ab7)
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 #include "pal.h"
35 #include "qgsmessagelog.h"
36 
37 using namespace pal;
38 
39 void heapsort( int *sid, int *id, const double* const x, int N )
40 {
41  unsigned int n = N, i = n / 2, parent, child;
42  int tx;
43  for ( ;; )
44  {
45  if ( i > 0 )
46  {
47  i--;
48  tx = sid[i];
49  }
50  else
51  {
52  n--;
53  if ( n == 0 ) return;
54  tx = sid[n];
55  sid[n] = sid[0];
56  }
57  parent = i;
58  child = i * 2 + 1;
59  while ( child < n )
60  {
61  if ( child + 1 < n && x[id[sid[child + 1]]] > x[id[sid[child]]] )
62  {
63  child++;
64  }
65  if ( x[id[sid[child]]] > x[id[tx]] )
66  {
67  sid[parent] = sid[child];
68  parent = child;
69  child = parent * 2 + 1;
70  }
71  else
72  {
73  break;
74  }
75  }
76  sid[parent] = tx;
77  }
78 }
79 
80 
81 void heapsort2( int *x, double* heap, int N )
82 {
83  unsigned int n = N, i = n / 2, parent, child;
84  double t;
85  int tx;
86  for ( ;; )
87  {
88  if ( i > 0 )
89  {
90  i--;
91  t = heap[i];
92  tx = x[i];
93  }
94  else
95  {
96  n--;
97  if ( n == 0 ) return;
98  t = heap[n];
99  tx = x[n];
100  heap[n] = heap[0];
101  x[n] = x[0];
102  }
103  parent = i;
104  child = i * 2 + 1;
105  while ( child < n )
106  {
107  if ( child + 1 < n && heap[child + 1] > heap[child] )
108  {
109  child++;
110  }
111  if ( heap[child] > t )
112  {
113  heap[parent] = heap[child];
114  x[parent] = x[child];
115  parent = child;
116  child = parent * 2 + 1;
117  }
118  else
119  {
120  break;
121  }
122  }
123  heap[parent] = t;
124  x[parent] = tx;
125  }
126 }
127 
128 bool GeomFunction::isSegIntersects( double x1, double y1, double x2, double y2, // 1st segment
129  double x3, double y3, double x4, double y4 ) // 2nd segment
130 {
131  return ( cross_product( x1, y1, x2, y2, x3, y3 ) * cross_product( x1, y1, x2, y2, x4, y4 ) < 0
132  && cross_product( x3, y3, x4, y4, x1, y1 ) * cross_product( x3, y3, x4, y4, x2, y2 ) < 0 );
133 }
134 
135 bool GeomFunction::computeLineIntersection( double x1, double y1, double x2, double y2, // 1st line (segment)
136  double x3, double y3, double x4, double y4, // 2nd line segment
137  double *x, double *y )
138 {
139 
140  double a1, a2, b1, b2, c1, c2;
141  double denom;
142 
143  a1 = y2 - y1;
144  b1 = x1 - x2;
145  c1 = x2 * y1 - x1 * y2;
146 
147  a2 = y4 - y3;
148  b2 = x3 - x4;
149  c2 = x4 * y3 - x3 * y4;
150 
151  denom = a1 * b2 - a2 * b1;
152  if ( qgsDoubleNear( denom, 0.0 ) )
153  {
154  return false;
155  }
156  else
157  {
158  *x = ( b1 * c2 - b2 * c1 ) / denom;
159  *y = ( a2 * c1 - a1 * c2 ) / denom;
160  }
161 
162  return true;
163 }
164 
165 int GeomFunction::convexHullId( int *id, const double* const x, const double* const y, int n, int *&cHull )
166 {
167  int i;
168 
169  cHull = new int[n];
170  for ( i = 0; i < n; i++ )
171  {
172  cHull[i] = i;
173  }
174 
175 
176  if ( n <= 3 ) return n;
177 
178  int* stack = new int[n];
179  double* tan = new double [n];
180  int ref;
181 
182  int second, top;
183  double result;
184 
185  // find the lowest y value
186  heapsort( cHull, id, y, n );
187 
188  // find the lowest x value from the lowest y
189  ref = 1;
190  while ( ref < n && qgsDoubleNear( y[id[cHull[ref]]], y[id[cHull[0]]] ) ) ref++;
191 
192  heapsort( cHull, id, x, ref );
193 
194  // the first point is now for sure in the hull as well as the ref one
195  for ( i = ref; i < n; i++ )
196  {
197  if ( qgsDoubleNear( y[id[cHull[i]]], y[id[cHull[0]]] ) )
198  tan[i] = FLT_MAX;
199  else
200  tan[i] = ( x[id[cHull[0]]] - x[id[cHull[i]]] ) / ( y[id[cHull[i]]] - y[id[cHull[0]]] );
201  }
202 
203  if ( ref < n )
204  heapsort2( cHull + ref, tan + ref, n - ref );
205 
206  // the second point is in too
207  stack[0] = cHull[0];
208  if ( ref == 1 )
209  {
210  stack[1] = cHull[1];
211  ref++;
212  }
213  else
214  stack[1] = cHull[ref-1];
215 
216 
217  top = 1;
218  second = 0;
219 
220  for ( i = ref; i < n; i++ )
221  {
222  result = cross_product( x[id[stack[second]]], y[id[stack[second]]],
223  x[id[stack[top]]], y[id[stack[top]]], x[id[cHull[i]]], y[id[cHull[i]]] );
224  // Coolineaire !! garder le plus éloigné
225  if ( qgsDoubleNear( result, 0.0 ) )
226  {
227  if ( dist_euc2d_sq( x[id[stack[second]]], y[id[stack[second]]], x[id[cHull[i]]], y[id[cHull[i]]] )
228  > dist_euc2d_sq( x[id[stack[second]]], y[id[stack[second]]], x[id[stack[top]]], y[id[stack[top]]] ) )
229  {
230  stack[top] = cHull[i];
231  }
232  }
233  else if ( result > 0 ) //convexe
234  {
235  second++;
236  top++;
237  stack[top] = cHull[i];
238  }
239  else
240  {
241  while ( result < 0 && top > 1 )
242  {
243  second--;
244  top--;
245  result = cross_product( x[id[stack[second]]],
246  y[id[stack[second]]], x[id[stack[top]]],
247  y[id[stack[top]]], x[id[cHull[i]]], y[id[cHull[i]]] );
248  }
249  second++;
250  top++;
251  stack[top] = cHull[i];
252  }
253  }
254 
255  for ( i = 0; i <= top; i++ )
256  {
257  cHull[i] = stack[i];
258  }
259 
260  delete[] stack;
261  delete[] tan;
262 
263  return top + 1;
264 }
265 
266 int GeomFunction::reorderPolygon( int nbPoints, double *x, double *y )
267 {
268  int inc = 0;
269  int *cHull;
270  int i;
271 
272  int *pts = new int[nbPoints];
273  for ( i = 0; i < nbPoints; i++ )
274  pts[i] = i;
275 
276  ( void )convexHullId( pts, x, y, nbPoints, cHull );
277 
278  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]] )
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 if ( pts[cHull[0]] < pts[cHull[1]] && pts[cHull[1]] > pts[cHull[2]] && pts[cHull[2]] < pts[cHull[0]] )
289  inc = 1;
290  else
291  {
292  // wrong cHull
293  delete[] cHull;
294  delete[] pts;
295  return -1;
296  }
297 
298  if ( inc == -1 ) // re-order points
299  {
300  double tmp;
301  int j;
302  for ( i = 0, j = nbPoints - 1; i <= j; i++, j-- )
303  {
304  tmp = x[i];
305  x[i] = x[j];
306  x[j] = tmp;
307 
308  tmp = y[i];
309  y[i] = y[j];
310  y[j] = tmp;
311  }
312  }
313 
314  delete[] cHull;
315  delete[] pts;
316 
317  return 0;
318 }
319 
320 bool GeomFunction::containsCandidate( const GEOSPreparedGeometry *geom, double x, double y, double width, double height, double alpha )
321 {
322  if ( !geom )
323  return false;
324 
325  GEOSContextHandle_t geosctxt = geosContext();
326  GEOSCoordSequence *coord = GEOSCoordSeq_create_r( geosctxt, 5, 2 );
327 
328  GEOSCoordSeq_setX_r( geosctxt, coord, 0, x );
329  GEOSCoordSeq_setY_r( geosctxt, coord, 0, y );
330  if ( !qgsDoubleNear( alpha, 0.0 ) )
331  {
332  double beta = alpha + ( M_PI / 2 );
333  double dx1 = cos( alpha ) * width;
334  double dy1 = sin( alpha ) * width;
335  double dx2 = cos( beta ) * height;
336  double dy2 = sin( beta ) * height;
337  GEOSCoordSeq_setX_r( geosctxt, coord, 1, x + dx1 );
338  GEOSCoordSeq_setY_r( geosctxt, coord, 1, y + dy1 );
339  GEOSCoordSeq_setX_r( geosctxt, coord, 2, x + dx1 + dx2 );
340  GEOSCoordSeq_setY_r( geosctxt, coord, 2, y + dy1 + dy2 );
341  GEOSCoordSeq_setX_r( geosctxt, coord, 3, x + dx2 );
342  GEOSCoordSeq_setY_r( geosctxt, coord, 3, y + dy2 );
343  }
344  else
345  {
346  GEOSCoordSeq_setX_r( geosctxt, coord, 1, x + width );
347  GEOSCoordSeq_setY_r( geosctxt, coord, 1, y );
348  GEOSCoordSeq_setX_r( geosctxt, coord, 2, x + width );
349  GEOSCoordSeq_setY_r( geosctxt, coord, 2, y + height );
350  GEOSCoordSeq_setX_r( geosctxt, coord, 3, x );
351  GEOSCoordSeq_setY_r( geosctxt, coord, 3, y + height );
352  }
353  //close ring
354  GEOSCoordSeq_setX_r( geosctxt, coord, 4, x );
355  GEOSCoordSeq_setY_r( geosctxt, coord, 4, y );
356 
357  try
358  {
359  GEOSGeometry* bboxGeos = GEOSGeom_createLinearRing_r( geosctxt, coord );
360  bool result = ( GEOSPreparedContainsProperly_r( geosctxt, geom, bboxGeos ) == 1 );
361  GEOSGeom_destroy_r( geosctxt, bboxGeos );
362  return result;
363  }
364  catch ( GEOSException &e )
365  {
366  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
367  return false;
368  }
369 }
370 
371 void GeomFunction::findLineCircleIntersection( double cx, double cy, double radius,
372  double x1, double y1, double x2, double y2,
373  double& xRes, double& yRes )
374 {
375  double dx = x2 - x1;
376  double dy = y2 - y1;
377 
378  double A = dx * dx + dy * dy;
379  double B = 2 * ( dx * ( x1 - cx ) + dy * ( y1 - cy ) );
380  double C = ( x1 - cx ) * ( x1 - cx ) + ( y1 - cy ) * ( y1 - cy ) - radius * radius;
381 
382  double det = B * B - 4 * A * C;
383  if ( A <= 0.0000001 || det < 0 )
384  // Should never happen, No real solutions.
385  return;
386 
387  if ( qgsDoubleNear( det, 0.0 ) )
388  {
389  // Could potentially happen.... One solution.
390  double t = -B / ( 2 * A );
391  xRes = x1 + t * dx;
392  yRes = y1 + t * dy;
393  }
394  else
395  {
396  // Two solutions.
397  // Always use the 1st one
398  // We only really have one solution here, as we know the line segment will start in the circle and end outside
399  double t = ( -B + sqrt( det ) ) / ( 2 * A );
400  xRes = x1 + t * dx;
401  yRes = y1 + t * dy;
402  }
403 }
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)
QString tr(const char *sourceText, const char *disambiguation, int n)
bool qgsDoubleNear(double a, double b, double epsilon=4 *DBL_EPSILON)
Compare two doubles (but allow some difference)
Definition: qgis.h:353
static double cross_product(double x1, double y1, double x2, double y2, double x3, double y3)
Definition: geomfunction.h:56
static bool containsCandidate(const GEOSPreparedGeometry *geom, double x, double y, double width, double height, double alpha)
Returns true if a GEOS prepared geometry totally contains a label candidate.
#define M_PI
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))
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)
GEOSContextHandle_t geosContext()
Get GEOS context handle to be used in all GEOS library calls with reentrant API.
Definition: pal.cpp:48
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:66