Branch data Line data Source code
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 : 0 : void heapsort( std::vector< int > &sid, int *id, const std::vector< double > &x, std::size_t N )
41 : : {
42 : 0 : std::size_t n = N;
43 : 0 : std::size_t i = n / 2;
44 : : std::size_t parent;
45 : : std::size_t child;
46 : : int tx;
47 : 0 : for ( ;; )
48 : : {
49 : 0 : if ( i > 0 )
50 : : {
51 : 0 : i--;
52 : 0 : tx = sid[i];
53 : 0 : }
54 : : else
55 : : {
56 : 0 : n--;
57 : 0 : if ( n == 0 ) return;
58 : 0 : tx = sid[n];
59 : 0 : sid[n] = sid[0];
60 : : }
61 : 0 : parent = i;
62 : 0 : child = i * 2 + 1;
63 : 0 : while ( child < n )
64 : : {
65 : 0 : if ( child + 1 < n && x[id[sid[child + 1]]] > x[id[sid[child]]] )
66 : : {
67 : 0 : child++;
68 : 0 : }
69 : 0 : if ( x[id[sid[child]]] > x[id[tx]] )
70 : : {
71 : 0 : sid[parent] = sid[child];
72 : 0 : parent = child;
73 : 0 : child = parent * 2 + 1;
74 : 0 : }
75 : : else
76 : : {
77 : 0 : break;
78 : : }
79 : : }
80 : 0 : sid[parent] = tx;
81 : : }
82 : : }
83 : :
84 : :
85 : 0 : void heapsort2( int *x, double *heap, std::size_t N )
86 : : {
87 : 0 : std::size_t n = N;
88 : 0 : std::size_t i = n / 2;
89 : : std::size_t parent;
90 : : std::size_t child;
91 : : double t;
92 : : int tx;
93 : 0 : for ( ;; )
94 : : {
95 : 0 : if ( i > 0 )
96 : : {
97 : 0 : i--;
98 : 0 : t = heap[i];
99 : 0 : tx = x[i];
100 : 0 : }
101 : : else
102 : : {
103 : 0 : n--;
104 : 0 : if ( n == 0 ) return;
105 : 0 : t = heap[n];
106 : 0 : tx = x[n];
107 : 0 : heap[n] = heap[0];
108 : 0 : x[n] = x[0];
109 : : }
110 : 0 : parent = i;
111 : 0 : child = i * 2 + 1;
112 : 0 : while ( child < n )
113 : : {
114 : 0 : if ( child + 1 < n && heap[child + 1] > heap[child] )
115 : : {
116 : 0 : child++;
117 : 0 : }
118 : 0 : if ( heap[child] > t )
119 : : {
120 : 0 : heap[parent] = heap[child];
121 : 0 : x[parent] = x[child];
122 : 0 : parent = child;
123 : 0 : child = parent * 2 + 1;
124 : 0 : }
125 : : else
126 : : {
127 : 0 : break;
128 : : }
129 : : }
130 : 0 : heap[parent] = t;
131 : 0 : x[parent] = tx;
132 : : }
133 : : }
134 : :
135 : 0 : 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 : 0 : return ( cross_product( x1, y1, x2, y2, x3, y3 ) * cross_product( x1, y1, x2, y2, x4, y4 ) < 0
139 : 0 : && cross_product( x3, y3, x4, y4, x1, y1 ) * cross_product( x3, y3, x4, y4, x2, y2 ) < 0 );
140 : : }
141 : :
142 : 0 : 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 : 0 : a1 = y2 - y1;
151 : 0 : b1 = x1 - x2;
152 : 0 : c1 = x2 * y1 - x1 * y2;
153 : :
154 : 0 : a2 = y4 - y3;
155 : 0 : b2 = x3 - x4;
156 : 0 : c2 = x4 * y3 - x3 * y4;
157 : :
158 : 0 : denom = a1 * b2 - a2 * b1;
159 : 0 : if ( qgsDoubleNear( denom, 0.0 ) )
160 : : {
161 : 0 : return false;
162 : : }
163 : : else
164 : : {
165 : 0 : *x = ( b1 * c2 - b2 * c1 ) / denom;
166 : 0 : *y = ( a2 * c1 - a1 * c2 ) / denom;
167 : : }
168 : :
169 : 0 : return true;
170 : 0 : }
171 : :
172 : 0 : std::vector< int > GeomFunction::convexHullId( std::vector< int > &id, const std::vector< double > &x, const std::vector< double > &y )
173 : : {
174 : 0 : std::vector< int > convexHull( x.size() );
175 : 0 : for ( std::size_t i = 0; i < x.size(); i++ )
176 : : {
177 : 0 : convexHull[i] = static_cast< int >( i );
178 : 0 : }
179 : :
180 : 0 : if ( x.size() <= 3 )
181 : 0 : return convexHull;
182 : :
183 : 0 : std::vector< int > stack( x.size() );
184 : 0 : std::vector< double > tan( x.size() );
185 : :
186 : : // find the lowest y value
187 : 0 : heapsort( convexHull, id.data(), y, y.size() );
188 : :
189 : : // find the lowest x value from the lowest y
190 : 0 : std::size_t ref = 1;
191 : 0 : while ( ref < x.size() && qgsDoubleNear( y[id[convexHull[ref]]], y[id[convexHull[0]]] ) )
192 : 0 : ref++;
193 : :
194 : 0 : 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 : 0 : for ( std::size_t i = ref; i < x.size(); i++ )
198 : : {
199 : 0 : if ( qgsDoubleNear( y[id[convexHull[i]]], y[id[convexHull[0]]] ) )
200 : 0 : tan[i] = FLT_MAX;
201 : : else
202 : 0 : tan[i] = ( x[id[convexHull[0]]] - x[id[convexHull[i]]] ) / ( y[id[convexHull[i]]] - y[id[convexHull[0]]] );
203 : 0 : }
204 : :
205 : 0 : if ( ref < x.size() )
206 : 0 : heapsort2( convexHull.data() + ref, tan.data() + ref, x.size() - ref );
207 : :
208 : : // the second point is in too
209 : 0 : stack[0] = convexHull[0];
210 : 0 : if ( ref == 1 )
211 : : {
212 : 0 : stack[1] = convexHull[1];
213 : 0 : ref++;
214 : 0 : }
215 : : else
216 : 0 : stack[1] = convexHull[ref - 1];
217 : :
218 : 0 : std::size_t top = 1;
219 : 0 : std::size_t second = 0;
220 : :
221 : 0 : for ( std::size_t i = ref; i < x.size(); i++ )
222 : : {
223 : 0 : double result = cross_product( x[id[stack[second]]], y[id[stack[second]]],
224 : 0 : x[id[stack[top]]], y[id[stack[top]]], x[id[convexHull[i]]], y[id[convexHull[i]]] );
225 : : // Coolineaire !! garder le plus éloigné
226 : 0 : if ( qgsDoubleNear( result, 0.0 ) )
227 : : {
228 : 0 : if ( dist_euc2d_sq( x[id[stack[second]]], y[id[stack[second]]], x[id[convexHull[i]]], y[id[convexHull[i]]] )
229 : 0 : > dist_euc2d_sq( x[id[stack[second]]], y[id[stack[second]]], x[id[stack[top]]], y[id[stack[top]]] ) )
230 : : {
231 : 0 : stack[top] = convexHull[i];
232 : 0 : }
233 : 0 : }
234 : 0 : else if ( result > 0 ) //convexe
235 : : {
236 : 0 : second++;
237 : 0 : top++;
238 : 0 : stack[top] = convexHull[i];
239 : 0 : }
240 : : else
241 : : {
242 : 0 : while ( result < 0 && top > 1 )
243 : : {
244 : 0 : second--;
245 : 0 : top--;
246 : 0 : result = cross_product( x[id[stack[second]]],
247 : 0 : y[id[stack[second]]], x[id[stack[top]]],
248 : 0 : y[id[stack[top]]], x[id[convexHull[i]]], y[id[convexHull[i]]] );
249 : : }
250 : 0 : second++;
251 : 0 : top++;
252 : 0 : stack[top] = convexHull[i];
253 : : }
254 : 0 : }
255 : :
256 : 0 : for ( std::size_t i = 0; i <= top; i++ )
257 : : {
258 : 0 : convexHull[i] = stack[i];
259 : 0 : }
260 : :
261 : 0 : convexHull.resize( top + 1 );
262 : 0 : return convexHull;
263 : 0 : }
264 : :
265 : 0 : bool GeomFunction::reorderPolygon( std::vector<double> &x, std::vector<double> &y )
266 : : {
267 : 0 : std::vector< int > pts( x.size() );
268 : 0 : for ( std::size_t i = 0; i < x.size(); i++ )
269 : 0 : pts[i] = static_cast< int >( i );
270 : :
271 : 0 : std::vector< int > convexHull = convexHullId( pts, x, y );
272 : :
273 : 0 : int inc = 0;
274 : 0 : if ( pts[convexHull[0]] < pts[convexHull[1]] && pts[convexHull[1]] < pts[convexHull[2]] )
275 : 0 : inc = 1;
276 : 0 : else if ( pts[convexHull[0]] > pts[convexHull[1]] && pts[convexHull[1]] > pts[convexHull[2]] )
277 : 0 : inc = -1;
278 : 0 : else if ( pts[convexHull[0]] > pts[convexHull[1]] && pts[convexHull[1]] < pts[convexHull[2]] && pts[convexHull[2]] < pts[convexHull[0]] )
279 : 0 : inc = 1;
280 : 0 : else if ( pts[convexHull[0]] > pts[convexHull[1]] && pts[convexHull[1]] < pts[convexHull[2]] && pts[convexHull[2]] > pts[convexHull[0]] )
281 : 0 : inc = -1;
282 : 0 : else if ( pts[convexHull[0]] < pts[convexHull[1]] && pts[convexHull[1]] > pts[convexHull[2]] && pts[convexHull[2]] > pts[convexHull[0]] )
283 : 0 : inc = -1;
284 : 0 : else if ( pts[convexHull[0]] < pts[convexHull[1]] && pts[convexHull[1]] > pts[convexHull[2]] && pts[convexHull[2]] < pts[convexHull[0]] )
285 : 0 : inc = 1;
286 : : else
287 : : {
288 : : // wrong cHull
289 : 0 : return false;
290 : : }
291 : :
292 : 0 : if ( inc == -1 ) // re-order points
293 : : {
294 : 0 : for ( std::size_t i = 0, j = x.size() - 1; i <= j; i++, j-- )
295 : : {
296 : 0 : std::swap( x[i], x[j] );
297 : 0 : std::swap( y[i], y[j] );
298 : 0 : }
299 : 0 : }
300 : 0 : return true;
301 : 0 : }
302 : :
303 : 0 : bool GeomFunction::containsCandidate( const GEOSPreparedGeometry *geom, double x, double y, double width, double height, double alpha )
304 : : {
305 : 0 : if ( !geom )
306 : 0 : return false;
307 : :
308 : : try
309 : : {
310 : 0 : GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
311 : 0 : GEOSCoordSequence *coord = GEOSCoordSeq_create_r( geosctxt, 5, 2 );
312 : :
313 : : #if GEOS_VERSION_MAJOR>3 || GEOS_VERSION_MINOR>=8
314 : 0 : 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 : 0 : if ( !qgsDoubleNear( alpha, 0.0 ) )
320 : : {
321 : 0 : double beta = alpha + M_PI_2;
322 : 0 : double dx1 = std::cos( alpha ) * width;
323 : 0 : double dy1 = std::sin( alpha ) * width;
324 : 0 : double dx2 = std::cos( beta ) * height;
325 : 0 : double dy2 = std::sin( beta ) * height;
326 : : #if GEOS_VERSION_MAJOR>3 || GEOS_VERSION_MINOR>=8
327 : 0 : GEOSCoordSeq_setXY_r( geosctxt, coord, 1, x + dx1, y + dy1 );
328 : 0 : GEOSCoordSeq_setXY_r( geosctxt, coord, 2, x + dx1 + dx2, y + dy1 + dy2 );
329 : 0 : 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 : 0 : }
339 : : else
340 : : {
341 : : #if GEOS_VERSION_MAJOR>3 || GEOS_VERSION_MINOR>=8
342 : 0 : GEOSCoordSeq_setXY_r( geosctxt, coord, 1, x + width, y );
343 : 0 : GEOSCoordSeq_setXY_r( geosctxt, coord, 2, x + width, y + height );
344 : 0 : 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 : 0 : 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 : 0 : geos::unique_ptr bboxGeos( GEOSGeom_createLinearRing_r( geosctxt, coord ) );
363 : 0 : bool result = ( GEOSPreparedContainsProperly_r( geosctxt, geom, bboxGeos.get() ) == 1 );
364 : 0 : return result;
365 : 0 : }
366 : : catch ( GEOSException &e )
367 : : {
368 : 0 : qWarning( "GEOS exception: %s", e.what() );
369 : : Q_NOWARN_UNREACHABLE_PUSH
370 : 0 : QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
371 : 0 : return false;
372 : : Q_NOWARN_UNREACHABLE_POP
373 : 0 : }
374 : 0 : return false;
375 : 0 : }
376 : :
377 : 0 : 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 : 0 : double multiplier = 1;
382 : 0 : 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 : 0 : multiplier = 10000;
387 : 0 : x1 *= multiplier;
388 : 0 : y1 *= multiplier;
389 : 0 : x2 *= multiplier;
390 : 0 : y2 *= multiplier;
391 : 0 : cx *= multiplier;
392 : 0 : cy *= multiplier;
393 : 0 : radius *= multiplier;
394 : 0 : }
395 : :
396 : 0 : double dx = x2 - x1;
397 : 0 : double dy = y2 - y1;
398 : :
399 : 0 : double A = dx * dx + dy * dy;
400 : 0 : double B = 2 * ( dx * ( x1 - cx ) + dy * ( y1 - cy ) );
401 : 0 : double C = ( x1 - cx ) * ( x1 - cx ) + ( y1 - cy ) * ( y1 - cy ) - radius * radius;
402 : :
403 : 0 : double det = B * B - 4 * A * C;
404 : 0 : if ( A <= 0.000000000001 || det < 0 )
405 : : // Should never happen, No real solutions.
406 : 0 : return;
407 : :
408 : 0 : if ( qgsDoubleNear( det, 0.0 ) )
409 : : {
410 : : // Could potentially happen.... One solution.
411 : 0 : double t = -B / ( 2 * A );
412 : 0 : xRes = x1 + t * dx;
413 : 0 : yRes = y1 + t * dy;
414 : 0 : }
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 : 0 : double t = ( -B + std::sqrt( det ) ) / ( 2 * A );
421 : 0 : xRes = x1 + t * dx;
422 : 0 : yRes = y1 + t * dy;
423 : : }
424 : :
425 : 0 : if ( multiplier != 1 )
426 : : {
427 : 0 : xRes /= multiplier;
428 : 0 : yRes /= multiplier;
429 : 0 : }
430 : 0 : }
|
