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