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