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