QGIS API Documentation  3.23.0-Master (eb871beae0)
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 #include <vector>
37 
38 using namespace pal;
39 
40 void heapsort( std::vector< int > &sid, int *id, const std::vector< double > &x, std::size_t N )
41 {
42  std::size_t n = N;
43  std::size_t i = n / 2;
44  std::size_t parent;
45  std::size_t child;
46  int tx;
47  for ( ;; )
48  {
49  if ( i > 0 )
50  {
51  i--;
52  tx = sid[i];
53  }
54  else
55  {
56  n--;
57  if ( n == 0 ) return;
58  tx = sid[n];
59  sid[n] = sid[0];
60  }
61  parent = i;
62  child = i * 2 + 1;
63  while ( child < n )
64  {
65  if ( child + 1 < n && x[id[sid[child + 1]]] > x[id[sid[child]]] )
66  {
67  child++;
68  }
69  if ( x[id[sid[child]]] > x[id[tx]] )
70  {
71  sid[parent] = sid[child];
72  parent = child;
73  child = parent * 2 + 1;
74  }
75  else
76  {
77  break;
78  }
79  }
80  sid[parent] = tx;
81  }
82 }
83 
84 
85 void heapsort2( int *x, double *heap, std::size_t N )
86 {
87  std::size_t n = N;
88  std::size_t i = n / 2;
89  std::size_t parent;
90  std::size_t child;
91  double t;
92  int tx;
93  for ( ;; )
94  {
95  if ( i > 0 )
96  {
97  i--;
98  t = heap[i];
99  tx = x[i];
100  }
101  else
102  {
103  n--;
104  if ( n == 0 ) return;
105  t = heap[n];
106  tx = x[n];
107  heap[n] = heap[0];
108  x[n] = x[0];
109  }
110  parent = i;
111  child = i * 2 + 1;
112  while ( child < n )
113  {
114  if ( child + 1 < n && heap[child + 1] > heap[child] )
115  {
116  child++;
117  }
118  if ( heap[child] > t )
119  {
120  heap[parent] = heap[child];
121  x[parent] = x[child];
122  parent = child;
123  child = parent * 2 + 1;
124  }
125  else
126  {
127  break;
128  }
129  }
130  heap[parent] = t;
131  x[parent] = tx;
132  }
133 }
134 
135 bool GeomFunction::isSegIntersects( double x1, double y1, double x2, double y2, // 1st segment
136  double x3, double y3, double x4, double y4 ) // 2nd segment
137 {
138  return ( cross_product( x1, y1, x2, y2, x3, y3 ) * cross_product( x1, y1, x2, y2, x4, y4 ) < 0
139  && cross_product( x3, y3, x4, y4, x1, y1 ) * cross_product( x3, y3, x4, y4, x2, y2 ) < 0 );
140 }
141 
142 bool GeomFunction::computeLineIntersection( double x1, double y1, double x2, double y2, // 1st line (segment)
143  double x3, double y3, double x4, double y4, // 2nd line segment
144  double *x, double *y )
145 {
146 
147  double a1, a2, b1, b2, c1, c2;
148  double denom;
149 
150  a1 = y2 - y1;
151  b1 = x1 - x2;
152  c1 = x2 * y1 - x1 * y2;
153 
154  a2 = y4 - y3;
155  b2 = x3 - x4;
156  c2 = x4 * y3 - x3 * y4;
157 
158  denom = a1 * b2 - a2 * b1;
159  if ( qgsDoubleNear( denom, 0.0 ) )
160  {
161  return false;
162  }
163  else
164  {
165  *x = ( b1 * c2 - b2 * c1 ) / denom;
166  *y = ( a2 * c1 - a1 * c2 ) / denom;
167  }
168 
169  return true;
170 }
171 
172 std::vector< int > GeomFunction::convexHullId( std::vector< int > &id, const std::vector< double > &x, const std::vector< double > &y )
173 {
174  std::vector< int > convexHull( x.size() );
175  for ( std::size_t i = 0; i < x.size(); i++ )
176  {
177  convexHull[i] = static_cast< int >( i );
178  }
179 
180  if ( x.size() <= 3 )
181  return convexHull;
182 
183  std::vector< int > stack( x.size() );
184  std::vector< double > tan( x.size() );
185 
186  // find the lowest y value
187  heapsort( convexHull, id.data(), y, y.size() );
188 
189  // find the lowest x value from the lowest y
190  std::size_t ref = 1;
191  while ( ref < x.size() && qgsDoubleNear( y[id[convexHull[ref]]], y[id[convexHull[0]]] ) )
192  ref++;
193 
194  heapsort( convexHull, id.data(), x, ref );
195 
196  // the first point is now for sure in the hull as well as the ref one
197  for ( std::size_t i = ref; i < x.size(); i++ )
198  {
199  if ( qgsDoubleNear( y[id[convexHull[i]]], y[id[convexHull[0]]] ) )
200  tan[i] = FLT_MAX;
201  else
202  tan[i] = ( x[id[convexHull[0]]] - x[id[convexHull[i]]] ) / ( y[id[convexHull[i]]] - y[id[convexHull[0]]] );
203  }
204 
205  if ( ref < x.size() )
206  heapsort2( convexHull.data() + ref, tan.data() + ref, x.size() - ref );
207 
208  // the second point is in too
209  stack[0] = convexHull[0];
210  if ( ref == 1 )
211  {
212  stack[1] = convexHull[1];
213  ref++;
214  }
215  else
216  stack[1] = convexHull[ref - 1];
217 
218  std::size_t top = 1;
219  std::size_t second = 0;
220 
221  for ( std::size_t i = ref; i < x.size(); i++ )
222  {
223  double result = cross_product( x[id[stack[second]]], y[id[stack[second]]],
224  x[id[stack[top]]], y[id[stack[top]]], x[id[convexHull[i]]], y[id[convexHull[i]]] );
225  // Coolineaire !! garder le plus éloigné
226  if ( qgsDoubleNear( result, 0.0 ) )
227  {
228  if ( dist_euc2d_sq( x[id[stack[second]]], y[id[stack[second]]], x[id[convexHull[i]]], y[id[convexHull[i]]] )
229  > dist_euc2d_sq( x[id[stack[second]]], y[id[stack[second]]], x[id[stack[top]]], y[id[stack[top]]] ) )
230  {
231  stack[top] = convexHull[i];
232  }
233  }
234  else if ( result > 0 ) //convexe
235  {
236  second++;
237  top++;
238  stack[top] = convexHull[i];
239  }
240  else
241  {
242  while ( result < 0 && top > 1 )
243  {
244  second--;
245  top--;
246  result = cross_product( x[id[stack[second]]],
247  y[id[stack[second]]], x[id[stack[top]]],
248  y[id[stack[top]]], x[id[convexHull[i]]], y[id[convexHull[i]]] );
249  }
250  second++;
251  top++;
252  stack[top] = convexHull[i];
253  }
254  }
255 
256  for ( std::size_t i = 0; i <= top; i++ )
257  {
258  convexHull[i] = stack[i];
259  }
260 
261  convexHull.resize( top + 1 );
262  return convexHull;
263 }
264 
265 bool GeomFunction::reorderPolygon( std::vector<double> &x, std::vector<double> &y )
266 {
267  std::vector< int > pts( x.size() );
268  for ( std::size_t i = 0; i < x.size(); i++ )
269  pts[i] = static_cast< int >( i );
270 
271  std::vector< int > convexHull = convexHullId( pts, x, y );
272 
273  int inc = 0;
274  if ( pts[convexHull[0]] < pts[convexHull[1]] && pts[convexHull[1]] < pts[convexHull[2]] )
275  inc = 1;
276  else if ( pts[convexHull[0]] > pts[convexHull[1]] && pts[convexHull[1]] > pts[convexHull[2]] )
277  inc = -1;
278  else if ( pts[convexHull[0]] > pts[convexHull[1]] && pts[convexHull[1]] < pts[convexHull[2]] && pts[convexHull[2]] < pts[convexHull[0]] )
279  inc = 1;
280  else if ( pts[convexHull[0]] > pts[convexHull[1]] && pts[convexHull[1]] < pts[convexHull[2]] && pts[convexHull[2]] > pts[convexHull[0]] )
281  inc = -1;
282  else if ( pts[convexHull[0]] < pts[convexHull[1]] && pts[convexHull[1]] > pts[convexHull[2]] && pts[convexHull[2]] > pts[convexHull[0]] )
283  inc = -1;
284  else if ( pts[convexHull[0]] < pts[convexHull[1]] && pts[convexHull[1]] > pts[convexHull[2]] && pts[convexHull[2]] < pts[convexHull[0]] )
285  inc = 1;
286  else
287  {
288  // wrong cHull
289  return false;
290  }
291 
292  if ( inc == -1 ) // re-order points
293  {
294  for ( std::size_t i = 0, j = x.size() - 1; i <= j; i++, j-- )
295  {
296  std::swap( x[i], x[j] );
297  std::swap( y[i], y[j] );
298  }
299  }
300  return true;
301 }
302 
303 bool GeomFunction::containsCandidate( const GEOSPreparedGeometry *geom, double x, double y, double width, double height, double alpha )
304 {
305  if ( !geom )
306  return false;
307 
308  try
309  {
310  GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
311  GEOSCoordSequence *coord = GEOSCoordSeq_create_r( geosctxt, 5, 2 );
312 
313 #if GEOS_VERSION_MAJOR>3 || GEOS_VERSION_MINOR>=8
314  GEOSCoordSeq_setXY_r( geosctxt, coord, 0, x, y );
315 #else
316  GEOSCoordSeq_setX_r( geosctxt, coord, 0, x );
317  GEOSCoordSeq_setY_r( geosctxt, coord, 0, y );
318 #endif
319  if ( !qgsDoubleNear( alpha, 0.0 ) )
320  {
321  const double beta = alpha + M_PI_2;
322  const double dx1 = std::cos( alpha ) * width;
323  const double dy1 = std::sin( alpha ) * width;
324  const double dx2 = std::cos( beta ) * height;
325  const double dy2 = std::sin( beta ) * height;
326 #if GEOS_VERSION_MAJOR>3 || GEOS_VERSION_MINOR>=8
327  GEOSCoordSeq_setXY_r( geosctxt, coord, 1, x + dx1, y + dy1 );
328  GEOSCoordSeq_setXY_r( geosctxt, coord, 2, x + dx1 + dx2, y + dy1 + dy2 );
329  GEOSCoordSeq_setXY_r( geosctxt, coord, 3, x + dx2, y + dy2 );
330 #else
331  GEOSCoordSeq_setX_r( geosctxt, coord, 1, x + dx1 );
332  GEOSCoordSeq_setY_r( geosctxt, coord, 1, y + dy1 );
333  GEOSCoordSeq_setX_r( geosctxt, coord, 2, x + dx1 + dx2 );
334  GEOSCoordSeq_setY_r( geosctxt, coord, 2, y + dy1 + dy2 );
335  GEOSCoordSeq_setX_r( geosctxt, coord, 3, x + dx2 );
336  GEOSCoordSeq_setY_r( geosctxt, coord, 3, y + dy2 );
337 #endif
338  }
339  else
340  {
341 #if GEOS_VERSION_MAJOR>3 || GEOS_VERSION_MINOR>=8
342  GEOSCoordSeq_setXY_r( geosctxt, coord, 1, x + width, y );
343  GEOSCoordSeq_setXY_r( geosctxt, coord, 2, x + width, y + height );
344  GEOSCoordSeq_setXY_r( geosctxt, coord, 3, x, y + height );
345 #else
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 #endif
353  }
354  //close ring
355 #if GEOS_VERSION_MAJOR>3 || GEOS_VERSION_MINOR>=8
356  GEOSCoordSeq_setXY_r( geosctxt, coord, 4, x, y );
357 #else
358  GEOSCoordSeq_setX_r( geosctxt, coord, 4, x );
359  GEOSCoordSeq_setY_r( geosctxt, coord, 4, y );
360 #endif
361 
362  geos::unique_ptr bboxGeos( GEOSGeom_createLinearRing_r( geosctxt, coord ) );
363  const bool result = ( GEOSPreparedContainsProperly_r( geosctxt, geom, bboxGeos.get() ) == 1 );
364  return result;
365  }
366  catch ( GEOSException &e )
367  {
368  qWarning( "GEOS exception: %s", e.what() );
370  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
371  return false;
373  }
374  return false;
375 }
376 
377 void GeomFunction::findLineCircleIntersection( double cx, double cy, double radius,
378  double x1, double y1, double x2, double y2,
379  double &xRes, double &yRes )
380 {
381  double multiplier = 1;
382  if ( radius < 10 )
383  {
384  // these calculations get unstable for small coordinates differences, e.g. as a result of map labeling in a geographic
385  // CRS
386  multiplier = 10000;
387  x1 *= multiplier;
388  y1 *= multiplier;
389  x2 *= multiplier;
390  y2 *= multiplier;
391  cx *= multiplier;
392  cy *= multiplier;
393  radius *= multiplier;
394  }
395 
396  const double dx = x2 - x1;
397  const double dy = y2 - y1;
398 
399  const double A = dx * dx + dy * dy;
400  const double B = 2 * ( dx * ( x1 - cx ) + dy * ( y1 - cy ) );
401  const double C = ( x1 - cx ) * ( x1 - cx ) + ( y1 - cy ) * ( y1 - cy ) - radius * radius;
402 
403  const double det = B * B - 4 * A * C;
404  if ( A <= 0.000000000001 || det < 0 )
405  // Should never happen, No real solutions.
406  return;
407 
408  if ( qgsDoubleNear( det, 0.0 ) )
409  {
410  // Could potentially happen.... One solution.
411  const double t = -B / ( 2 * A );
412  xRes = x1 + t * dx;
413  yRes = y1 + t * dy;
414  }
415  else
416  {
417  // Two solutions.
418  // Always use the 1st one
419  // We only really have one solution here, as we know the line segment will start in the circle and end outside
420  const double t = ( -B + std::sqrt( det ) ) / ( 2 * A );
421  xRes = x1 + t * dx;
422  yRes = y1 + t * dy;
423  }
424 
425  if ( multiplier != 1 )
426  {
427  xRes /= multiplier;
428  yRes /= multiplier;
429  }
430 }
static GEOSContextHandle_t getGEOSHandler()
Definition: qgsgeos.cpp:3222
static void logMessage(const QString &message, const QString &tag=QString(), Qgis::MessageLevel level=Qgis::MessageLevel::Warning, bool notifyUser=true)
Adds a message to the log instance (and creates it if necessary).
static bool reorderPolygon(std::vector< double > &x, std::vector< double > &y)
Reorder points to have cross prod ((x,y)[i], (x,y)[i+1), point) > 0 when point is outside.
static std::vector< int > convexHullId(std::vector< int > &id, const std::vector< double > &x, const std::vector< double > &y)
Compute the convex hull in O(n·log(n))
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 double cross_product(double x1, double y1, double x2, double y2, double x3, double y3)
Definition: geomfunction.h:62
static double dist_euc2d_sq(double x1, double y1, double x2, double y2)
Definition: geomfunction.h:72
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.
static void findLineCircleIntersection(double cx, double cy, double radius, double x1, double y1, double x2, double y2, double &xRes, double &yRes)
void heapsort(std::vector< int > &sid, int *id, const std::vector< double > &x, std::size_t N)
void heapsort2(int *x, double *heap, std::size_t N)
std::unique_ptr< GEOSGeometry, GeosDeleter > unique_ptr
Scoped GEOS pointer.
Definition: qgsgeos.h:79
#define Q_NOWARN_UNREACHABLE_PUSH
Definition: qgis.h:2001
bool qgsDoubleNear(double a, double b, double epsilon=4 *std::numeric_limits< double >::epsilon())
Compare two doubles (but allow some difference)
Definition: qgis.h:1504
#define Q_NOWARN_UNREACHABLE_POP
Definition: qgis.h:2002