QGIS API Documentation 3.32.0-Lima (311a8cb8a6)
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
38using namespace pal;
39
40void 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
85void 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
135bool 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
142bool 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
172std::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
265bool 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
303bool 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 GEOSCoordSeq_setXY_r( geosctxt, coord, 0, x, y );
314 if ( !qgsDoubleNear( alpha, 0.0 ) )
315 {
316 const double beta = alpha + M_PI_2;
317 const double dx1 = std::cos( alpha ) * width;
318 const double dy1 = std::sin( alpha ) * width;
319 const double dx2 = std::cos( beta ) * height;
320 const double dy2 = std::sin( beta ) * height;
321 GEOSCoordSeq_setXY_r( geosctxt, coord, 1, x + dx1, y + dy1 );
322 GEOSCoordSeq_setXY_r( geosctxt, coord, 2, x + dx1 + dx2, y + dy1 + dy2 );
323 GEOSCoordSeq_setXY_r( geosctxt, coord, 3, x + dx2, y + dy2 );
324 }
325 else
326 {
327 GEOSCoordSeq_setXY_r( geosctxt, coord, 1, x + width, y );
328 GEOSCoordSeq_setXY_r( geosctxt, coord, 2, x + width, y + height );
329 GEOSCoordSeq_setXY_r( geosctxt, coord, 3, x, y + height );
330 }
331 //close ring
332 GEOSCoordSeq_setXY_r( geosctxt, coord, 4, x, y );
333
334 geos::unique_ptr bboxGeos( GEOSGeom_createLinearRing_r( geosctxt, coord ) );
335 const bool result = ( GEOSPreparedContainsProperly_r( geosctxt, geom, bboxGeos.get() ) == 1 );
336 return result;
337 }
338 catch ( GEOSException &e )
339 {
340 qWarning( "GEOS exception: %s", e.what() );
342 QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
343 return false;
345 }
346 return false;
347}
348
349void GeomFunction::findLineCircleIntersection( double cx, double cy, double radius,
350 double x1, double y1, double x2, double y2,
351 double &xRes, double &yRes )
352{
353 double multiplier = 1;
354 if ( radius < 10 )
355 {
356 // these calculations get unstable for small coordinates differences, e.g. as a result of map labeling in a geographic
357 // CRS
358 multiplier = 10000;
359 x1 *= multiplier;
360 y1 *= multiplier;
361 x2 *= multiplier;
362 y2 *= multiplier;
363 cx *= multiplier;
364 cy *= multiplier;
365 radius *= multiplier;
366 }
367
368 const double dx = x2 - x1;
369 const double dy = y2 - y1;
370
371 const double A = dx * dx + dy * dy;
372 const double B = 2 * ( dx * ( x1 - cx ) + dy * ( y1 - cy ) );
373 const double C = ( x1 - cx ) * ( x1 - cx ) + ( y1 - cy ) * ( y1 - cy ) - radius * radius;
374
375 const double det = B * B - 4 * A * C;
376 if ( A <= 0.000000000001 || det < 0 )
377 // Should never happen, No real solutions.
378 return;
379
380 if ( qgsDoubleNear( det, 0.0 ) )
381 {
382 // Could potentially happen.... One solution.
383 const double t = -B / ( 2 * A );
384 xRes = x1 + t * dx;
385 yRes = y1 + t * dy;
386 }
387 else
388 {
389 // Two solutions.
390 // Always use the 1st one
391 // We only really have one solution here, as we know the line segment will start in the circle and end outside
392 const double t = ( -B + std::sqrt( det ) ) / ( 2 * A );
393 xRes = x1 + t * dx;
394 yRes = y1 + t * dy;
395 }
396
397 if ( multiplier != 1 )
398 {
399 xRes /= multiplier;
400 yRes /= multiplier;
401 }
402}
static GEOSContextHandle_t getGEOSHandler()
Definition: qgsgeos.cpp:3451
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 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:74
#define Q_NOWARN_UNREACHABLE_PUSH
Definition: qgis.h:4573
bool qgsDoubleNear(double a, double b, double epsilon=4 *std::numeric_limits< double >::epsilon())
Compare two doubles (but allow some difference)
Definition: qgis.h:3988
#define Q_NOWARN_UNREACHABLE_POP
Definition: qgis.h:4574