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 "pointset.h"
31 : : #include "util.h"
32 : : #include "pal.h"
33 : : #include "geomfunction.h"
34 : : #include "qgsgeos.h"
35 : : #include "qgsmessagelog.h"
36 : : #include "qgsgeometryutils.h"
37 : : #include <qglobal.h>
38 : :
39 : : using namespace pal;
40 : :
41 : 0 : PointSet::PointSet()
42 : 0 : {
43 : 0 : nbPoints = 0;
44 : 0 : type = -1;
45 : 0 : }
46 : :
47 : 0 : PointSet::PointSet( int nbPoints, double *x, double *y )
48 : 0 : : nbPoints( nbPoints )
49 : 0 : , type( GEOS_POLYGON )
50 : 0 : {
51 : 0 : this->x.resize( nbPoints );
52 : 0 : this->y.resize( nbPoints );
53 : : int i;
54 : :
55 : 0 : for ( i = 0; i < nbPoints; i++ )
56 : : {
57 : 0 : this->x[i] = x[i];
58 : 0 : this->y[i] = y[i];
59 : 0 : }
60 : :
61 : 0 : }
62 : :
63 : 0 : PointSet::PointSet( double aX, double aY )
64 : 0 : : type( GEOS_POINT )
65 : 0 : , xmin( aX )
66 : 0 : , xmax( aY )
67 : 0 : , ymin( aX )
68 : 0 : , ymax( aY )
69 : 0 : {
70 : 0 : nbPoints = 1;
71 : 0 : x.resize( 1 );
72 : 0 : y.resize( 1 );
73 : 0 : x[0] = aX;
74 : 0 : y[0] = aY;
75 : 0 : }
76 : :
77 : 0 : PointSet::PointSet( const PointSet &ps )
78 : 0 : : xmin( ps.xmin )
79 : 0 : , xmax( ps.xmax )
80 : 0 : , ymin( ps.ymin )
81 : 0 : , ymax( ps.ymax )
82 : 0 : {
83 : 0 : nbPoints = ps.nbPoints;
84 : 0 : x = ps.x;
85 : 0 : y = ps.y;
86 : :
87 : 0 : convexHull = ps.convexHull;
88 : :
89 : 0 : type = ps.type;
90 : :
91 : 0 : holeOf = ps.holeOf;
92 : :
93 : 0 : if ( ps.mGeos )
94 : : {
95 : 0 : mGeos = GEOSGeom_clone_r( QgsGeos::getGEOSHandler(), ps.mGeos );
96 : 0 : mOwnsGeom = true;
97 : 0 : }
98 : 0 : }
99 : :
100 : 0 : void PointSet::createGeosGeom() const
101 : : {
102 : 0 : GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
103 : :
104 : 0 : bool needClose = false;
105 : 0 : if ( type == GEOS_POLYGON && ( !qgsDoubleNear( x[0], x[ nbPoints - 1] ) || !qgsDoubleNear( y[0], y[ nbPoints - 1 ] ) ) )
106 : : {
107 : 0 : needClose = true;
108 : 0 : }
109 : :
110 : 0 : GEOSCoordSequence *coord = GEOSCoordSeq_create_r( geosctxt, nbPoints + ( needClose ? 1 : 0 ), 2 );
111 : 0 : for ( int i = 0; i < nbPoints; ++i )
112 : : {
113 : : #if GEOS_VERSION_MAJOR>3 || GEOS_VERSION_MINOR>=8
114 : 0 : GEOSCoordSeq_setXY_r( geosctxt, coord, i, x[i], y[i] );
115 : : #else
116 : : GEOSCoordSeq_setX_r( geosctxt, coord, i, x[i] );
117 : : GEOSCoordSeq_setY_r( geosctxt, coord, i, y[i] );
118 : : #endif
119 : 0 : }
120 : :
121 : : //close ring if needed
122 : 0 : if ( needClose )
123 : : {
124 : : #if GEOS_VERSION_MAJOR>3 || GEOS_VERSION_MINOR>=8
125 : 0 : GEOSCoordSeq_setXY_r( geosctxt, coord, nbPoints, x[0], y[0] );
126 : : #else
127 : : GEOSCoordSeq_setX_r( geosctxt, coord, nbPoints, x[0] );
128 : : GEOSCoordSeq_setY_r( geosctxt, coord, nbPoints, y[0] );
129 : : #endif
130 : 0 : }
131 : :
132 : 0 : switch ( type )
133 : : {
134 : : case GEOS_POLYGON:
135 : 0 : mGeos = GEOSGeom_createPolygon_r( geosctxt, GEOSGeom_createLinearRing_r( geosctxt, coord ), nullptr, 0 );
136 : 0 : break;
137 : :
138 : : case GEOS_LINESTRING:
139 : 0 : mGeos = GEOSGeom_createLineString_r( geosctxt, coord );
140 : 0 : break;
141 : :
142 : : case GEOS_POINT:
143 : 0 : mGeos = GEOSGeom_createPoint_r( geosctxt, coord );
144 : 0 : break;
145 : : }
146 : :
147 : 0 : mOwnsGeom = true;
148 : 0 : }
149 : :
150 : 0 : const GEOSPreparedGeometry *PointSet::preparedGeom() const
151 : : {
152 : 0 : if ( !mGeos )
153 : 0 : createGeosGeom();
154 : :
155 : 0 : if ( !mPreparedGeom )
156 : : {
157 : 0 : mPreparedGeom = GEOSPrepare_r( QgsGeos::getGEOSHandler(), mGeos );
158 : 0 : }
159 : 0 : return mPreparedGeom;
160 : : }
161 : :
162 : 0 : void PointSet::invalidateGeos()
163 : : {
164 : 0 : GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
165 : 0 : if ( mOwnsGeom ) // delete old geometry if we own it
166 : 0 : GEOSGeom_destroy_r( geosctxt, mGeos );
167 : 0 : GEOSPreparedGeom_destroy_r( geosctxt, mPreparedGeom );
168 : 0 : mOwnsGeom = false;
169 : 0 : mGeos = nullptr;
170 : 0 : if ( mGeosPreparedBoundary )
171 : : {
172 : 0 : GEOSPreparedGeom_destroy_r( geosctxt, mGeosPreparedBoundary );
173 : 0 : mGeosPreparedBoundary = nullptr;
174 : 0 : }
175 : :
176 : 0 : if ( mMultipartPreparedGeos )
177 : : {
178 : 0 : GEOSPreparedGeom_destroy_r( geosctxt, mMultipartPreparedGeos );
179 : 0 : mMultipartPreparedGeos = nullptr;
180 : 0 : }
181 : 0 : if ( mMultipartGeos )
182 : : {
183 : 0 : GEOSGeom_destroy_r( geosctxt, mMultipartGeos );
184 : 0 : mMultipartGeos = nullptr;
185 : 0 : }
186 : :
187 : 0 : mPreparedGeom = nullptr;
188 : 0 : mLength = -1;
189 : 0 : mArea = -1;
190 : 0 : }
191 : :
192 : 0 : PointSet::~PointSet()
193 : 0 : {
194 : 0 : GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
195 : :
196 : 0 : if ( mGeos && mOwnsGeom )
197 : : {
198 : 0 : GEOSGeom_destroy_r( geosctxt, mGeos );
199 : 0 : mGeos = nullptr;
200 : 0 : }
201 : 0 : GEOSPreparedGeom_destroy_r( geosctxt, mPreparedGeom );
202 : :
203 : 0 : if ( mGeosPreparedBoundary )
204 : : {
205 : 0 : GEOSPreparedGeom_destroy_r( geosctxt, mGeosPreparedBoundary );
206 : 0 : mGeosPreparedBoundary = nullptr;
207 : 0 : }
208 : :
209 : 0 : if ( mMultipartPreparedGeos )
210 : : {
211 : 0 : GEOSPreparedGeom_destroy_r( geosctxt, mMultipartPreparedGeos );
212 : 0 : mMultipartPreparedGeos = nullptr;
213 : 0 : }
214 : 0 : if ( mMultipartGeos )
215 : : {
216 : 0 : GEOSGeom_destroy_r( geosctxt, mMultipartGeos );
217 : 0 : mMultipartGeos = nullptr;
218 : 0 : }
219 : :
220 : 0 : deleteCoords();
221 : 0 : }
222 : :
223 : 0 : void PointSet::deleteCoords()
224 : 0 : {
225 : 0 : x.clear();
226 : 0 : y.clear();
227 : 0 : }
228 : :
229 : 0 : std::unique_ptr<PointSet> PointSet::extractShape( int nbPtSh, int imin, int imax, int fps, int fpe, double fptx, double fpty )
230 : : {
231 : : int i, j;
232 : :
233 : 0 : std::unique_ptr<PointSet> newShape = std::make_unique< PointSet >();
234 : 0 : newShape->type = GEOS_POLYGON;
235 : 0 : newShape->nbPoints = nbPtSh;
236 : 0 : newShape->x.resize( newShape->nbPoints );
237 : 0 : newShape->y.resize( newShape->nbPoints );
238 : 0 :
239 : 0 : // new shape # 1 from imin to imax
240 : 0 : for ( j = 0, i = imin; i != ( imax + 1 ) % nbPoints; i = ( i + 1 ) % nbPoints, j++ )
241 : : {
242 : 0 : newShape->x[j] = x[i];
243 : 0 : newShape->y[j] = y[i];
244 : 0 : }
245 : 0 : // is the cutting point a new one ?
246 : 0 : if ( fps != fpe )
247 : 0 : {
248 : 0 : // yes => so add it
249 : 0 : newShape->x[j] = fptx;
250 : 0 : newShape->y[j] = fpty;
251 : 0 : }
252 : :
253 : 0 : return newShape;
254 : 0 : }
255 : :
256 : 0 : std::unique_ptr<PointSet> PointSet::clone() const
257 : : {
258 : 0 : return std::unique_ptr< PointSet>( new PointSet( *this ) );
259 : 0 : }
260 : :
261 : 0 : bool PointSet::containsPoint( double x, double y ) const
262 : : {
263 : 0 : GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
264 : : try
265 : : {
266 : : #if GEOS_VERSION_MAJOR>3 || GEOS_VERSION_MINOR>=8
267 : 0 : geos::unique_ptr point( GEOSGeom_createPointFromXY_r( geosctxt, x, y ) );
268 : : #else
269 : : GEOSCoordSequence *seq = GEOSCoordSeq_create_r( geosctxt, 1, 2 );
270 : : GEOSCoordSeq_setX_r( geosctxt, seq, 0, x );
271 : : GEOSCoordSeq_setY_r( geosctxt, seq, 0, y );
272 : : geos::unique_ptr point( GEOSGeom_createPoint_r( geosctxt, seq ) );
273 : : #endif
274 : 0 : bool result = ( GEOSPreparedContainsProperly_r( geosctxt, preparedGeom(), point.get() ) == 1 );
275 : :
276 : 0 : return result;
277 : 0 : }
278 : : catch ( GEOSException &e )
279 : : {
280 : 0 : qWarning( "GEOS exception: %s", e.what() );
281 : 0 : QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
282 : 0 : return false;
283 : 0 : }
284 : :
285 : 0 : }
286 : :
287 : 0 : bool PointSet::containsLabelCandidate( double x, double y, double width, double height, double alpha ) const
288 : : {
289 : 0 : return GeomFunction::containsCandidate( preparedGeom(), x, y, width, height, alpha );
290 : : }
291 : :
292 : 0 : QLinkedList<PointSet *> PointSet::splitPolygons( PointSet *inputShape, double labelWidth, double labelHeight )
293 : : {
294 : : int j;
295 : :
296 : 0 : double bestArea = 0;
297 : :
298 : : double b;
299 : :
300 : 0 : int holeS = -1; // hole start and end points
301 : 0 : int holeE = -1;
302 : :
303 : 0 : int retainedPt = -1;
304 : :
305 : 0 : const double labelArea = labelWidth * labelHeight;
306 : :
307 : 0 : QLinkedList<PointSet *> inputShapes;
308 : 0 : inputShapes.push_back( inputShape );
309 : 0 : QLinkedList<PointSet *> outputShapes;
310 : :
311 : 0 : while ( !inputShapes.isEmpty() )
312 : : {
313 : 0 : PointSet *shape = inputShapes.takeFirst();
314 : :
315 : 0 : const std::vector< double > &x = shape->x;
316 : 0 : const std::vector< double > &y = shape->y;
317 : 0 : const int nbp = shape->nbPoints;
318 : 0 : std::vector< int > pts( nbp );
319 : 0 : for ( int i = 0; i < nbp; i++ )
320 : : {
321 : 0 : pts[i] = i;
322 : 0 : }
323 : :
324 : : // compute convex hull
325 : 0 : shape->convexHull = GeomFunction::convexHullId( pts, x, y );
326 : :
327 : 0 : bestArea = 0;
328 : 0 : retainedPt = -1;
329 : :
330 : : // lookup for a hole
331 : 0 : for ( std::size_t ihs = 0; ihs < shape->convexHull.size(); ihs++ )
332 : : {
333 : : // ihs->ihn => cHull'seg
334 : 0 : std::size_t ihn = ( ihs + 1 ) % shape->convexHull.size();
335 : :
336 : 0 : int ips = shape->convexHull[ihs];
337 : 0 : int ipn = ( ips + 1 ) % nbp;
338 : 0 : if ( ipn != shape->convexHull[ihn] ) // next point on shape is not the next point on cHull => there is a hole here !
339 : : {
340 : 0 : double bestcp = 0;
341 : 0 : int pt = -1;
342 : : // lookup for the deepest point in the hole
343 : 0 : for ( int i = ips; i != shape->convexHull[ihn]; i = ( i + 1 ) % nbp )
344 : : {
345 : 0 : double cp = std::fabs( GeomFunction::cross_product( x[shape->convexHull[ihs]], y[shape->convexHull[ihs]],
346 : 0 : x[shape->convexHull[ihn]], y[shape->convexHull[ihn]],
347 : 0 : x[i], y[i] ) );
348 : 0 : if ( cp - bestcp > EPSILON )
349 : : {
350 : 0 : bestcp = cp;
351 : 0 : pt = i;
352 : 0 : }
353 : 0 : }
354 : :
355 : 0 : if ( pt != -1 )
356 : : {
357 : : // compute the ihs->ihn->pt triangle's area
358 : 0 : const double base = GeomFunction::dist_euc2d( x[shape->convexHull[ihs]], y[shape->convexHull[ihs]],
359 : 0 : x[shape->convexHull[ihn]], y[shape->convexHull[ihn]] );
360 : :
361 : 0 : b = GeomFunction::dist_euc2d( x[shape->convexHull[ihs]], y[shape->convexHull[ihs]],
362 : 0 : x[pt], y[pt] );
363 : :
364 : 0 : const double c = GeomFunction::dist_euc2d( x[shape->convexHull[ihn]], y[shape->convexHull[ihn]],
365 : 0 : x[pt], y[pt] );
366 : :
367 : 0 : const double s = ( base + b + c ) / 2; // s = half perimeter
368 : 0 : double area = s * ( s - base ) * ( s - b ) * ( s - c );
369 : 0 : if ( area < 0 )
370 : 0 : area = -area;
371 : :
372 : : // retain the biggest area
373 : 0 : if ( area - bestArea > EPSILON )
374 : : {
375 : 0 : bestArea = area;
376 : 0 : retainedPt = pt;
377 : 0 : holeS = ihs;
378 : 0 : holeE = ihn;
379 : 0 : }
380 : 0 : }
381 : 0 : }
382 : 0 : }
383 : :
384 : : // we have a hole, its area, and the deppest point in hole
385 : : // we're going to find the second point to cup the shape
386 : : // holeS = hole starting point
387 : : // holeE = hole ending point
388 : : // retainedPt = deppest point in hole
389 : : // bestArea = area of triangle HoleS->holeE->retainedPoint
390 : 0 : bestArea = std::sqrt( bestArea );
391 : 0 : double cx, cy, dx, dy, ex, ey, fx, fy, seg_length, ptx = 0, pty = 0, fptx = 0, fpty = 0;
392 : 0 : int ps = -1, pe = -1, fps = -1, fpe = -1;
393 : 0 : if ( retainedPt >= 0 && bestArea > labelArea ) // there is a hole so we'll cut the shape in two new shape (only if hole area is bigger than twice labelArea)
394 : : {
395 : 0 : double c = std::numeric_limits<double>::max();
396 : :
397 : : // iterate on all shape points except points which are in the hole
398 : : bool isValid;
399 : : int k, l;
400 : 0 : for ( int i = ( shape->convexHull[holeE] + 1 ) % nbp; i != ( shape->convexHull[holeS] - 1 + nbp ) % nbp; i = j )
401 : : {
402 : 0 : j = ( i + 1 ) % nbp; // i->j is shape segment not in hole
403 : :
404 : : // compute distance between retainedPoint and segment
405 : : // whether perpendicular distance (if retaindPoint is fronting segment i->j)
406 : : // or distance between retainedPt and i or j (choose the nearest)
407 : 0 : seg_length = GeomFunction::dist_euc2d( x[i], y[i], x[j], y[j] );
408 : 0 : cx = ( x[i] + x[j] ) / 2.0;
409 : 0 : cy = ( y[i] + y[j] ) / 2.0;
410 : 0 : dx = cy - y[i];
411 : 0 : dy = cx - x[i];
412 : :
413 : 0 : ex = cx - dx;
414 : 0 : ey = cy + dy;
415 : 0 : fx = cx + dx;
416 : 0 : fy = cy - dy;
417 : :
418 : 0 : if ( seg_length < EPSILON || std::fabs( ( b = GeomFunction::cross_product( ex, ey, fx, fy, x[retainedPt], y[retainedPt] ) / ( seg_length ) ) ) > ( seg_length / 2 ) ) // retainedPt is not fronting i->j
419 : : {
420 : 0 : if ( ( ex = GeomFunction::dist_euc2d_sq( x[i], y[i], x[retainedPt], y[retainedPt] ) ) < ( ey = GeomFunction::dist_euc2d_sq( x[j], y[j], x[retainedPt], y[retainedPt] ) ) )
421 : : {
422 : 0 : b = ex;
423 : 0 : ps = i;
424 : 0 : pe = i;
425 : 0 : }
426 : : else
427 : : {
428 : 0 : b = ey;
429 : 0 : ps = j;
430 : 0 : pe = j;
431 : : }
432 : 0 : }
433 : : else // point fronting i->j => compute pependicular distance => create a new point
434 : : {
435 : 0 : b = GeomFunction::cross_product( x[i], y[i], x[j], y[j], x[retainedPt], y[retainedPt] ) / seg_length;
436 : 0 : b *= b;
437 : 0 : ps = i;
438 : 0 : pe = j;
439 : :
440 : 0 : if ( !GeomFunction::computeLineIntersection( x[i], y[i], x[j], y[j], x[retainedPt], y[retainedPt], x[retainedPt] - dx, y[retainedPt] + dy, &ptx, &pty ) )
441 : : {
442 : : //error - it should intersect the line
443 : 0 : }
444 : : }
445 : :
446 : 0 : isValid = true;
447 : : double pointX, pointY;
448 : 0 : if ( ps == pe )
449 : : {
450 : 0 : pointX = x[pe];
451 : 0 : pointY = y[pe];
452 : 0 : }
453 : : else
454 : : {
455 : 0 : pointX = ptx;
456 : 0 : pointY = pty;
457 : : }
458 : :
459 : 0 : for ( k = shape->convexHull[holeS]; k != shape->convexHull[holeE]; k = ( k + 1 ) % nbp )
460 : : {
461 : 0 : l = ( k + 1 ) % nbp;
462 : 0 : if ( GeomFunction::isSegIntersects( x[retainedPt], y[retainedPt], pointX, pointY, x[k], y[k], x[l], y[l] ) )
463 : : {
464 : 0 : isValid = false;
465 : 0 : break;
466 : : }
467 : 0 : }
468 : :
469 : :
470 : 0 : if ( isValid && b < c )
471 : : {
472 : 0 : c = b;
473 : 0 : fps = ps;
474 : 0 : fpe = pe;
475 : 0 : fptx = ptx;
476 : 0 : fpty = pty;
477 : 0 : }
478 : 0 : } // for point which are not in hole
479 : :
480 : : // we will cut the shapeu in two new shapes, one from [retainedPoint] to [newPoint] and one form [newPoint] to [retainedPoint]
481 : 0 : int imin = retainedPt;
482 : 0 : int imax = ( ( ( fps < retainedPt && fpe < retainedPt ) || ( fps > retainedPt && fpe > retainedPt ) ) ? std::min( fps, fpe ) : std::max( fps, fpe ) );
483 : :
484 : : int nbPtSh1, nbPtSh2; // how many points in new shapes ?
485 : 0 : if ( imax > imin )
486 : 0 : nbPtSh1 = imax - imin + 1 + ( fpe != fps );
487 : : else
488 : 0 : nbPtSh1 = imax + nbp - imin + 1 + ( fpe != fps );
489 : :
490 : 0 : if ( ( imax == fps ? fpe : fps ) < imin )
491 : 0 : nbPtSh2 = imin - ( imax == fps ? fpe : fps ) + 1 + ( fpe != fps );
492 : : else
493 : 0 : nbPtSh2 = imin + nbp - ( imax == fps ? fpe : fps ) + 1 + ( fpe != fps );
494 : :
495 : 0 : if ( retainedPt == -1 || fps == -1 || fpe == -1 )
496 : : {
497 : 0 : if ( shape->parent )
498 : 0 : delete shape;
499 : 0 : }
500 : : // check for useless splitting
501 : 0 : else if ( imax == imin || nbPtSh1 <= 2 || nbPtSh2 <= 2 || nbPtSh1 == nbp || nbPtSh2 == nbp )
502 : : {
503 : 0 : outputShapes.append( shape );
504 : 0 : }
505 : : else
506 : : {
507 : :
508 : 0 : PointSet *newShape = shape->extractShape( nbPtSh1, imin, imax, fps, fpe, fptx, fpty ).release();
509 : :
510 : 0 : if ( shape->parent )
511 : 0 : newShape->parent = shape->parent;
512 : : else
513 : 0 : newShape->parent = shape;
514 : :
515 : 0 : inputShapes.append( newShape );
516 : :
517 : 0 : if ( imax == fps )
518 : 0 : imax = fpe;
519 : : else
520 : 0 : imax = fps;
521 : :
522 : 0 : newShape = shape->extractShape( nbPtSh2, imax, imin, fps, fpe, fptx, fpty ).release();
523 : :
524 : 0 : if ( shape->parent )
525 : 0 : newShape->parent = shape->parent;
526 : : else
527 : 0 : newShape->parent = shape;
528 : :
529 : 0 : inputShapes.append( newShape );
530 : :
531 : 0 : if ( shape->parent )
532 : 0 : delete shape;
533 : : }
534 : 0 : }
535 : : else
536 : : {
537 : 0 : outputShapes.append( shape );
538 : : }
539 : 0 : }
540 : 0 : return outputShapes;
541 : 0 : }
542 : :
543 : 0 : void PointSet::extendLineByDistance( double startDistance, double endDistance, double smoothDistance )
544 : : {
545 : 0 : if ( nbPoints < 2 )
546 : 0 : return;
547 : :
548 : 0 : double x0 = x[0];
549 : 0 : double y0 = y[0];
550 : 0 : if ( startDistance > 0 )
551 : : {
552 : : // trace forward by smoothDistance
553 : 0 : double x1 = x[1];
554 : 0 : double y1 = y[1];
555 : :
556 : 0 : double distanceConsumed = 0;
557 : 0 : double lastX = x0;
558 : 0 : double lastY = y0;
559 : 0 : for ( int i = 1; i < nbPoints; ++i )
560 : : {
561 : 0 : double thisX = x[i];
562 : 0 : double thisY = y[i];
563 : 0 : const double thisSegmentLength = std::sqrt( ( thisX - lastX ) * ( thisX - lastX ) + ( thisY - lastY ) * ( thisY - lastY ) );
564 : 0 : distanceConsumed += thisSegmentLength;
565 : 0 : if ( distanceConsumed >= smoothDistance )
566 : : {
567 : 0 : double c = ( distanceConsumed - smoothDistance ) / thisSegmentLength;
568 : 0 : x1 = lastX + c * ( thisX - lastX );
569 : 0 : y1 = lastY + c * ( thisY - lastY );
570 : 0 : break;
571 : : }
572 : 0 : lastX = thisX;
573 : 0 : lastY = thisY;
574 : 0 : }
575 : :
576 : 0 : const double distance = std::sqrt( ( x1 - x0 ) * ( x1 - x0 ) + ( y1 - y0 ) * ( y1 - y0 ) );
577 : 0 : const double extensionFactor = ( startDistance + distance ) / distance;
578 : 0 : const QgsPointXY newStart = QgsGeometryUtils::interpolatePointOnLine( x1, y1, x0, y0, extensionFactor );
579 : 0 : x0 = newStart.x();
580 : 0 : y0 = newStart.y();
581 : : // defer actually changing the stored start until we've calculated the new end point
582 : 0 : }
583 : :
584 : 0 : if ( endDistance > 0 )
585 : : {
586 : 0 : double xend0 = x[nbPoints - 1];
587 : 0 : double yend0 = y[nbPoints - 1];
588 : 0 : double xend1 = x[nbPoints - 2];
589 : 0 : double yend1 = y[nbPoints - 2];
590 : :
591 : : // trace backward by smoothDistance
592 : 0 : double distanceConsumed = 0;
593 : 0 : double lastX = x0;
594 : 0 : double lastY = y0;
595 : 0 : for ( int i = nbPoints - 2; i >= 0; --i )
596 : : {
597 : 0 : double thisX = x[i];
598 : 0 : double thisY = y[i];
599 : 0 : const double thisSegmentLength = std::sqrt( ( thisX - lastX ) * ( thisX - lastX ) + ( thisY - lastY ) * ( thisY - lastY ) );
600 : 0 : distanceConsumed += thisSegmentLength;
601 : 0 : if ( distanceConsumed >= smoothDistance )
602 : : {
603 : 0 : double c = ( distanceConsumed - smoothDistance ) / thisSegmentLength;
604 : 0 : xend1 = lastX + c * ( thisX - lastX );
605 : 0 : yend1 = lastY + c * ( thisY - lastY );
606 : 0 : break;
607 : : }
608 : 0 : lastX = thisX;
609 : 0 : lastY = thisY;
610 : 0 : }
611 : :
612 : 0 : const double distance = std::sqrt( ( xend1 - xend0 ) * ( xend1 - xend0 ) + ( yend1 - yend0 ) * ( yend1 - yend0 ) );
613 : 0 : const double extensionFactor = ( endDistance + distance ) / distance;
614 : 0 : const QgsPointXY newEnd = QgsGeometryUtils::interpolatePointOnLine( xend1, yend1, xend0, yend0, extensionFactor );
615 : 0 : x.emplace_back( newEnd.x() );
616 : 0 : y.emplace_back( newEnd.y() );
617 : 0 : nbPoints++;
618 : 0 : }
619 : :
620 : 0 : if ( startDistance > 0 )
621 : : {
622 : 0 : x.insert( x.begin(), x0 );
623 : 0 : y.insert( y.begin(), y0 );
624 : 0 : nbPoints++;
625 : 0 : }
626 : :
627 : 0 : invalidateGeos();
628 : 0 : }
629 : :
630 : 0 : OrientedConvexHullBoundingBox PointSet::computeConvexHullOrientedBoundingBox( bool &ok )
631 : : {
632 : 0 : ok = false;
633 : : double bbox[4]; // xmin, ymin, xmax, ymax
634 : :
635 : : double alpha;
636 : : int alpha_d;
637 : :
638 : : double alpha_seg;
639 : :
640 : : double d1, d2;
641 : :
642 : : double bb[16]; // {ax, ay, bx, by, cx, cy, dx, dy, ex, ey, fx, fy, gx, gy, hx, hy}}
643 : :
644 : 0 : double best_area = std::numeric_limits<double>::max();
645 : 0 : double best_alpha = -1;
646 : : double best_bb[16];
647 : 0 : double best_length = 0;
648 : 0 : double best_width = 0;
649 : :
650 : :
651 : 0 : bbox[0] = std::numeric_limits<double>::max();
652 : 0 : bbox[1] = std::numeric_limits<double>::max();
653 : 0 : bbox[2] = std::numeric_limits<double>::lowest();
654 : 0 : bbox[3] = std::numeric_limits<double>::lowest();
655 : :
656 : 0 : for ( std::size_t i = 0; i < convexHull.size(); i++ )
657 : : {
658 : 0 : if ( x[convexHull[i]] < bbox[0] )
659 : 0 : bbox[0] = x[convexHull[i]];
660 : :
661 : 0 : if ( x[convexHull[i]] > bbox[2] )
662 : 0 : bbox[2] = x[convexHull[i]];
663 : :
664 : 0 : if ( y[convexHull[i]] < bbox[1] )
665 : 0 : bbox[1] = y[convexHull[i]];
666 : :
667 : 0 : if ( y[convexHull[i]] > bbox[3] )
668 : 0 : bbox[3] = y[convexHull[i]];
669 : 0 : }
670 : :
671 : : OrientedConvexHullBoundingBox finalBb;
672 : :
673 : 0 : const double dref = bbox[2] - bbox[0];
674 : 0 : if ( qgsDoubleNear( dref, 0 ) )
675 : : {
676 : 0 : ok = false;
677 : 0 : return finalBb;
678 : : }
679 : :
680 : 0 : for ( alpha_d = 0; alpha_d < 90; alpha_d++ )
681 : : {
682 : 0 : alpha = alpha_d * M_PI / 180.0;
683 : 0 : d1 = std::cos( alpha ) * dref;
684 : 0 : d2 = std::sin( alpha ) * dref;
685 : :
686 : 0 : bb[0] = bbox[0];
687 : 0 : bb[1] = bbox[3]; // ax, ay
688 : :
689 : 0 : bb[4] = bbox[0];
690 : 0 : bb[5] = bbox[1]; // cx, cy
691 : :
692 : 0 : bb[8] = bbox[2];
693 : 0 : bb[9] = bbox[1]; // ex, ey
694 : :
695 : 0 : bb[12] = bbox[2];
696 : 0 : bb[13] = bbox[3]; // gx, gy
697 : :
698 : :
699 : 0 : bb[2] = bb[0] + d1;
700 : 0 : bb[3] = bb[1] + d2; // bx, by
701 : 0 : bb[6] = bb[4] - d2;
702 : 0 : bb[7] = bb[5] + d1; // dx, dy
703 : 0 : bb[10] = bb[8] - d1;
704 : 0 : bb[11] = bb[9] - d2; // fx, fy
705 : 0 : bb[14] = bb[12] + d2;
706 : 0 : bb[15] = bb[13] - d1; // hx, hy
707 : :
708 : : // adjust all points
709 : 0 : for ( int i = 0; i < 16; i += 4 )
710 : : {
711 : :
712 : 0 : alpha_seg = ( ( i / 4 > 0 ? ( i / 4 ) - 1 : 3 ) ) * M_PI_2 + alpha;
713 : :
714 : 0 : double best_cp = std::numeric_limits<double>::max();
715 : :
716 : 0 : for ( std::size_t j = 0; j < convexHull.size(); j++ )
717 : : {
718 : 0 : const double cp = GeomFunction::cross_product( bb[i + 2], bb[i + 3], bb[i], bb[i + 1], x[convexHull[j]], y[convexHull[j]] );
719 : 0 : if ( cp < best_cp )
720 : : {
721 : 0 : best_cp = cp;
722 : 0 : }
723 : 0 : }
724 : :
725 : 0 : const double distNearestPoint = best_cp / dref;
726 : :
727 : 0 : d1 = std::cos( alpha_seg ) * distNearestPoint;
728 : 0 : d2 = std::sin( alpha_seg ) * distNearestPoint;
729 : :
730 : 0 : bb[i] += d1; // x
731 : 0 : bb[i + 1] += d2; // y
732 : 0 : bb[i + 2] += d1; // x
733 : 0 : bb[i + 3] += d2; // y
734 : 0 : }
735 : :
736 : : // compute and compare AREA
737 : 0 : const double width = GeomFunction::cross_product( bb[6], bb[7], bb[4], bb[5], bb[12], bb[13] ) / dref;
738 : 0 : const double length = GeomFunction::cross_product( bb[2], bb[3], bb[0], bb[1], bb[8], bb[9] ) / dref;
739 : :
740 : 0 : double area = width * length;
741 : :
742 : 0 : if ( area < 0 )
743 : 0 : area *= -1;
744 : :
745 : :
746 : 0 : if ( best_area - area > EPSILON )
747 : : {
748 : 0 : best_area = area;
749 : 0 : best_length = length;
750 : 0 : best_width = width;
751 : 0 : best_alpha = alpha;
752 : 0 : memcpy( best_bb, bb, sizeof( double ) * 16 );
753 : 0 : }
754 : 0 : }
755 : :
756 : : // best bbox is defined
757 : 0 : for ( int i = 0; i < 16; i = i + 4 )
758 : : {
759 : 0 : GeomFunction::computeLineIntersection( best_bb[i], best_bb[i + 1], best_bb[i + 2], best_bb[i + 3],
760 : 0 : best_bb[( i + 4 ) % 16], best_bb[( i + 5 ) % 16], best_bb[( i + 6 ) % 16], best_bb[( i + 7 ) % 16],
761 : 0 : &finalBb.x[int ( i / 4 )], &finalBb.y[int ( i / 4 )] );
762 : 0 : }
763 : :
764 : 0 : finalBb.alpha = best_alpha;
765 : 0 : finalBb.width = best_width;
766 : 0 : finalBb.length = best_length;
767 : :
768 : 0 : ok = true;
769 : 0 : return finalBb;
770 : 0 : }
771 : :
772 : 0 : double PointSet::minDistanceToPoint( double px, double py, double *rx, double *ry ) const
773 : : {
774 : 0 : if ( !mGeos )
775 : 0 : createGeosGeom();
776 : :
777 : 0 : if ( !mGeos )
778 : 0 : return 0;
779 : :
780 : 0 : GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
781 : : try
782 : : {
783 : : #if GEOS_VERSION_MAJOR>3 || GEOS_VERSION_MINOR>=8
784 : 0 : geos::unique_ptr geosPt( GEOSGeom_createPointFromXY_r( geosctxt, px, py ) );
785 : : #else
786 : : GEOSCoordSequence *coord = GEOSCoordSeq_create_r( geosctxt, 1, 2 );
787 : : GEOSCoordSeq_setX_r( geosctxt, coord, 0, px );
788 : : GEOSCoordSeq_setY_r( geosctxt, coord, 0, py );
789 : : geos::unique_ptr geosPt( GEOSGeom_createPoint_r( geosctxt, coord ) );
790 : : #endif
791 : 0 : int type = GEOSGeomTypeId_r( geosctxt, mGeos );
792 : 0 : const GEOSGeometry *extRing = nullptr;
793 : : #if GEOS_VERSION_MAJOR>3 || GEOS_VERSION_MINOR>=9
794 : 0 : const GEOSPreparedGeometry *preparedExtRing = nullptr;
795 : : #endif
796 : :
797 : 0 : if ( type != GEOS_POLYGON )
798 : : {
799 : 0 : extRing = mGeos;
800 : : #if GEOS_VERSION_MAJOR>3 || GEOS_VERSION_MINOR>=9
801 : 0 : preparedExtRing = preparedGeom();
802 : : #endif
803 : 0 : }
804 : : else
805 : : {
806 : : //for polygons, we want distance to exterior ring (not an interior point)
807 : 0 : extRing = GEOSGetExteriorRing_r( geosctxt, mGeos );
808 : : #if GEOS_VERSION_MAJOR>3 || ( GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR>=9 )
809 : 0 : if ( ! mGeosPreparedBoundary )
810 : : {
811 : 0 : mGeosPreparedBoundary = GEOSPrepare_r( geosctxt, extRing );
812 : 0 : }
813 : 0 : preparedExtRing = mGeosPreparedBoundary;
814 : : #endif
815 : : }
816 : :
817 : : #if GEOS_VERSION_MAJOR>3 || ( GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR>=9 )
818 : 0 : geos::coord_sequence_unique_ptr nearestCoord( GEOSPreparedNearestPoints_r( geosctxt, preparedExtRing, geosPt.get() ) );
819 : : #else
820 : : geos::coord_sequence_unique_ptr nearestCoord( GEOSNearestPoints_r( geosctxt, extRing, geosPt.get() ) );
821 : : #endif
822 : : double nx;
823 : : double ny;
824 : : #if GEOS_VERSION_MAJOR>3 || GEOS_VERSION_MINOR>=8
825 : 0 : unsigned int nPoints = 0;
826 : 0 : GEOSCoordSeq_getSize_r( geosctxt, nearestCoord.get(), &nPoints );
827 : 0 : if ( nPoints == 0 )
828 : 0 : return 0;
829 : :
830 : 0 : ( void )GEOSCoordSeq_getXY_r( geosctxt, nearestCoord.get(), 0, &nx, &ny );
831 : : #else
832 : : ( void )GEOSCoordSeq_getX_r( geosctxt, nearestCoord.get(), 0, &nx );
833 : : ( void )GEOSCoordSeq_getY_r( geosctxt, nearestCoord.get(), 0, &ny );
834 : : #endif
835 : :
836 : 0 : if ( rx )
837 : 0 : *rx = nx;
838 : 0 : if ( ry )
839 : 0 : *ry = ny;
840 : :
841 : 0 : return GeomFunction::dist_euc2d_sq( px, py, nx, ny );
842 : 0 : }
843 : : catch ( GEOSException &e )
844 : : {
845 : 0 : qWarning( "GEOS exception: %s", e.what() );
846 : 0 : QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
847 : 0 : return 0;
848 : 0 : }
849 : 0 : }
850 : :
851 : 0 : void PointSet::getCentroid( double &px, double &py, bool forceInside ) const
852 : : {
853 : 0 : if ( !mGeos )
854 : 0 : createGeosGeom();
855 : :
856 : 0 : if ( !mGeos )
857 : 0 : return;
858 : :
859 : : try
860 : : {
861 : 0 : GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
862 : 0 : geos::unique_ptr centroidGeom( GEOSGetCentroid_r( geosctxt, mGeos ) );
863 : 0 : if ( centroidGeom )
864 : : {
865 : 0 : const GEOSCoordSequence *coordSeq = GEOSGeom_getCoordSeq_r( geosctxt, centroidGeom.get() );
866 : : #if GEOS_VERSION_MAJOR>3 || GEOS_VERSION_MINOR>=8
867 : 0 : unsigned int nPoints = 0;
868 : 0 : GEOSCoordSeq_getSize_r( geosctxt, coordSeq, &nPoints );
869 : 0 : if ( nPoints == 0 )
870 : 0 : return;
871 : 0 : GEOSCoordSeq_getXY_r( geosctxt, coordSeq, 0, &px, &py );
872 : : #else
873 : : GEOSCoordSeq_getX_r( geosctxt, coordSeq, 0, &px );
874 : : GEOSCoordSeq_getY_r( geosctxt, coordSeq, 0, &py );
875 : : #endif
876 : 0 : }
877 : :
878 : : // check if centroid inside in polygon
879 : 0 : if ( forceInside && !containsPoint( px, py ) )
880 : : {
881 : 0 : geos::unique_ptr pointGeom( GEOSPointOnSurface_r( geosctxt, mGeos ) );
882 : :
883 : 0 : if ( pointGeom )
884 : : {
885 : 0 : const GEOSCoordSequence *coordSeq = GEOSGeom_getCoordSeq_r( geosctxt, pointGeom.get() );
886 : : #if GEOS_VERSION_MAJOR>3 || GEOS_VERSION_MINOR>=8
887 : 0 : unsigned int nPoints = 0;
888 : 0 : GEOSCoordSeq_getSize_r( geosctxt, coordSeq, &nPoints );
889 : 0 : if ( nPoints == 0 )
890 : 0 : return;
891 : :
892 : 0 : GEOSCoordSeq_getXY_r( geosctxt, coordSeq, 0, &px, &py );
893 : : #else
894 : : GEOSCoordSeq_getX_r( geosctxt, coordSeq, 0, &px );
895 : : GEOSCoordSeq_getY_r( geosctxt, coordSeq, 0, &py );
896 : : #endif
897 : 0 : }
898 : 0 : }
899 : 0 : }
900 : : catch ( GEOSException &e )
901 : : {
902 : 0 : qWarning( "GEOS exception: %s", e.what() );
903 : 0 : QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
904 : : return;
905 : 0 : }
906 : 0 : }
907 : :
908 : 0 : bool PointSet::boundingBoxIntersects( const PointSet *other ) const
909 : : {
910 : 0 : double x1 = ( xmin > other->xmin ? xmin : other->xmin );
911 : 0 : double x2 = ( xmax < other->xmax ? xmax : other->xmax );
912 : 0 : if ( x1 > x2 )
913 : 0 : return false;
914 : 0 : double y1 = ( ymin > other->ymin ? ymin : other->ymin );
915 : 0 : double y2 = ( ymax < other->ymax ? ymax : other->ymax );
916 : 0 : return y1 <= y2;
917 : 0 : }
918 : :
919 : 0 : void PointSet::getPointByDistance( double *d, double *ad, double dl, double *px, double *py )
920 : : {
921 : : int i;
922 : : double dx, dy, di;
923 : : double distr;
924 : :
925 : 0 : i = 0;
926 : 0 : if ( dl >= 0 )
927 : : {
928 : 0 : while ( i < nbPoints && ad[i] <= dl ) i++;
929 : 0 : i--;
930 : 0 : }
931 : :
932 : 0 : if ( i < nbPoints - 1 )
933 : : {
934 : 0 : if ( dl < 0 )
935 : : {
936 : 0 : dx = x[nbPoints - 1] - x[0];
937 : 0 : dy = y[nbPoints - 1] - y[0];
938 : 0 : di = std::sqrt( dx * dx + dy * dy );
939 : 0 : }
940 : : else
941 : : {
942 : 0 : dx = x[i + 1] - x[i];
943 : 0 : dy = y[i + 1] - y[i];
944 : 0 : di = d[i];
945 : : }
946 : :
947 : 0 : distr = dl - ad[i];
948 : 0 : *px = x[i] + dx * distr / di;
949 : 0 : *py = y[i] + dy * distr / di;
950 : 0 : }
951 : : else // just select last point...
952 : : {
953 : 0 : *px = x[i];
954 : 0 : *py = y[i];
955 : : }
956 : 0 : }
957 : :
958 : 0 : const GEOSGeometry *PointSet::geos() const
959 : : {
960 : 0 : if ( !mGeos )
961 : 0 : createGeosGeom();
962 : :
963 : 0 : return mGeos;
964 : : }
965 : :
966 : 0 : double PointSet::length() const
967 : : {
968 : 0 : if ( mLength >= 0 )
969 : 0 : return mLength;
970 : :
971 : 0 : if ( !mGeos )
972 : 0 : createGeosGeom();
973 : :
974 : 0 : if ( !mGeos )
975 : 0 : return -1;
976 : :
977 : 0 : GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
978 : :
979 : : try
980 : : {
981 : 0 : ( void )GEOSLength_r( geosctxt, mGeos, &mLength );
982 : 0 : return mLength;
983 : 0 : }
984 : : catch ( GEOSException &e )
985 : : {
986 : 0 : qWarning( "GEOS exception: %s", e.what() );
987 : 0 : QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
988 : 0 : return -1;
989 : 0 : }
990 : 0 : }
991 : :
992 : 0 : double PointSet::area() const
993 : : {
994 : 0 : if ( mArea >= 0 )
995 : 0 : return mArea;
996 : :
997 : 0 : if ( !mGeos )
998 : 0 : createGeosGeom();
999 : :
1000 : 0 : if ( !mGeos )
1001 : 0 : return -1;
1002 : :
1003 : 0 : GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
1004 : :
1005 : : try
1006 : : {
1007 : 0 : ( void )GEOSArea_r( geosctxt, mGeos, &mArea );
1008 : 0 : mArea = std::fabs( mArea );
1009 : 0 : return mArea;
1010 : 0 : }
1011 : : catch ( GEOSException &e )
1012 : : {
1013 : 0 : qWarning( "GEOS exception: %s", e.what() );
1014 : 0 : QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
1015 : 0 : return -1;
1016 : 0 : }
1017 : 0 : }
1018 : :
1019 : 0 : bool PointSet::isClosed() const
1020 : : {
1021 : 0 : return qgsDoubleNear( x[0], x[nbPoints - 1] ) && qgsDoubleNear( y[0], y[nbPoints - 1] );
1022 : : }
1023 : :
1024 : 0 : QString PointSet::toWkt() const
1025 : : {
1026 : 0 : if ( !mGeos )
1027 : 0 : createGeosGeom();
1028 : :
1029 : 0 : GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
1030 : :
1031 : : try
1032 : : {
1033 : 0 : GEOSWKTWriter *writer = GEOSWKTWriter_create_r( geosctxt );
1034 : :
1035 : 0 : char *wkt = GEOSWKTWriter_write_r( geosctxt, writer, mGeos );
1036 : 0 : const QString res( wkt );
1037 : :
1038 : 0 : GEOSFree_r( geosctxt, wkt );
1039 : :
1040 : 0 : GEOSWKTWriter_destroy_r( geosctxt, writer );
1041 : 0 : writer = nullptr;
1042 : :
1043 : 0 : return res;
1044 : 0 : }
1045 : : catch ( GEOSException &e )
1046 : : {
1047 : 0 : qWarning( "GEOS exception: %s", e.what() );
1048 : 0 : QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
1049 : 0 : return QString();
1050 : 0 : }
1051 : 0 : }
|
