QGIS API Documentation 3.34.0-Prizren (ffbdd678812)
Loading...
Searching...
No Matches
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:3465
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)
static double dist_euc2d_sq(double x1, double y1, double x2, double y2)
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:4917
bool qgsDoubleNear(double a, double b, double epsilon=4 *std::numeric_limits< double >::epsilon())
Compare two doubles (but allow some difference)
Definition qgis.h:4332
#define Q_NOWARN_UNREACHABLE_POP
Definition qgis.h:4918