Branch data Line data Source code
1 : : /***************************************************************************
2 : : qgscircularstring.cpp
3 : : -----------------------
4 : : begin : September 2014
5 : : copyright : (C) 2014 by Marco Hugentobler
6 : : email : marco at sourcepole dot ch
7 : : ***************************************************************************/
8 : :
9 : : /***************************************************************************
10 : : * *
11 : : * This program is free software; you can redistribute it and/or modify *
12 : : * it under the terms of the GNU General Public License as published by *
13 : : * the Free Software Foundation; either version 2 of the License, or *
14 : : * (at your option) any later version. *
15 : : * *
16 : : ***************************************************************************/
17 : :
18 : : #include "qgscircularstring.h"
19 : : #include "qgsapplication.h"
20 : : #include "qgscoordinatetransform.h"
21 : : #include "qgsgeometryutils.h"
22 : : #include "qgslinestring.h"
23 : : #include "qgsmaptopixel.h"
24 : : #include "qgspoint.h"
25 : : #include "qgswkbptr.h"
26 : : #include "qgslogger.h"
27 : : #include "qgsgeometrytransformer.h"
28 : : #include "qgsfeedback.h"
29 : :
30 : : #include <QJsonObject>
31 : : #include <QPainter>
32 : : #include <QPainterPath>
33 : : #include <memory>
34 : : #include <nlohmann/json.hpp>
35 : :
36 : 1316 : QgsCircularString::QgsCircularString()
37 : 2632 : {
38 : 1316 : mWkbType = QgsWkbTypes::CircularString;
39 : 1316 : }
40 : :
41 : 9 : QgsCircularString::QgsCircularString( const QgsPoint &p1, const QgsPoint &p2, const QgsPoint &p3 )
42 : 18 : {
43 : : //get wkb type from first point
44 : 9 : bool hasZ = p1.is3D();
45 : 9 : bool hasM = p1.isMeasure();
46 : 9 : mWkbType = QgsWkbTypes::CircularString;
47 : :
48 : 9 : mX.resize( 3 );
49 : 9 : mX[ 0 ] = p1.x();
50 : 9 : mX[ 1 ] = p2.x();
51 : 9 : mX[ 2 ] = p3.x();
52 : 9 : mY.resize( 3 );
53 : 9 : mY[ 0 ] = p1.y();
54 : 9 : mY[ 1 ] = p2.y();
55 : 9 : mY[ 2 ] = p3.y();
56 : 9 : if ( hasZ )
57 : : {
58 : 4 : mWkbType = QgsWkbTypes::addZ( mWkbType );
59 : 4 : mZ.resize( 3 );
60 : 4 : mZ[ 0 ] = p1.z();
61 : 4 : mZ[ 1 ] = p2.z();
62 : 4 : mZ[ 2 ] = p3.z();
63 : 4 : }
64 : 9 : if ( hasM )
65 : : {
66 : 4 : mWkbType = QgsWkbTypes::addM( mWkbType );
67 : 4 : mM.resize( 3 );
68 : 4 : mM[ 0 ] = p1.m();
69 : 4 : mM[ 1 ] = p2.m();
70 : 4 : mM[ 2 ] = p3.m();
71 : 4 : }
72 : 9 : }
73 : :
74 : 5 : QgsCircularString QgsCircularString::fromTwoPointsAndCenter( const QgsPoint &p1, const QgsPoint &p2, const QgsPoint ¢er, const bool useShortestArc )
75 : : {
76 : 5 : const QgsPoint midPoint = QgsGeometryUtils::segmentMidPointFromCenter( p1, p2, center, useShortestArc );
77 : 5 : return QgsCircularString( p1, midPoint, p2 );
78 : 5 : }
79 : :
80 : 117 : bool QgsCircularString::equals( const QgsCurve &other ) const
81 : : {
82 : 117 : const QgsCircularString *otherLine = dynamic_cast< const QgsCircularString * >( &other );
83 : 117 : if ( !otherLine )
84 : 1 : return false;
85 : :
86 : 116 : if ( mWkbType != otherLine->mWkbType )
87 : 2 : return false;
88 : :
89 : 114 : if ( mX.count() != otherLine->mX.count() )
90 : 2 : return false;
91 : :
92 : 484 : for ( int i = 0; i < mX.count(); ++i )
93 : : {
94 : 382 : if ( !qgsDoubleNear( mX.at( i ), otherLine->mX.at( i ) )
95 : 382 : || !qgsDoubleNear( mY.at( i ), otherLine->mY.at( i ) ) )
96 : 6 : return false;
97 : :
98 : 376 : if ( is3D() && !qgsDoubleNear( mZ.at( i ), otherLine->mZ.at( i ) ) )
99 : 2 : return false;
100 : :
101 : 374 : if ( isMeasure() && !qgsDoubleNear( mM.at( i ), otherLine->mM.at( i ) ) )
102 : 2 : return false;
103 : 372 : }
104 : :
105 : 102 : return true;
106 : 117 : }
107 : :
108 : 8 : QgsCircularString *QgsCircularString::createEmptyWithSameType() const
109 : : {
110 : 8 : auto result = std::make_unique< QgsCircularString >();
111 : 8 : result->mWkbType = mWkbType;
112 : 8 : return result.release();
113 : 8 : }
114 : :
115 : 503 : QString QgsCircularString::geometryType() const
116 : : {
117 : 1006 : return QStringLiteral( "CircularString" );
118 : : }
119 : :
120 : 2 : int QgsCircularString::dimension() const
121 : : {
122 : 2 : return 1;
123 : : }
124 : :
125 : 309 : QgsCircularString *QgsCircularString::clone() const
126 : : {
127 : 309 : return new QgsCircularString( *this );
128 : 0 : }
129 : :
130 : 473 : void QgsCircularString::clear()
131 : : {
132 : 473 : mWkbType = QgsWkbTypes::CircularString;
133 : 473 : mX.clear();
134 : 473 : mY.clear();
135 : 473 : mZ.clear();
136 : 473 : mM.clear();
137 : 473 : clearCache();
138 : 473 : }
139 : :
140 : 27 : QgsRectangle QgsCircularString::calculateBoundingBox() const
141 : : {
142 : 27 : QgsRectangle bbox;
143 : 27 : int nPoints = numPoints();
144 : 53 : for ( int i = 0; i < ( nPoints - 2 ) ; i += 2 )
145 : : {
146 : 26 : if ( i == 0 )
147 : : {
148 : 17 : bbox = segmentBoundingBox( QgsPoint( mX[i], mY[i] ), QgsPoint( mX[i + 1], mY[i + 1] ), QgsPoint( mX[i + 2], mY[i + 2] ) );
149 : 17 : }
150 : : else
151 : : {
152 : 9 : QgsRectangle segmentBox = segmentBoundingBox( QgsPoint( mX[i], mY[i] ), QgsPoint( mX[i + 1], mY[i + 1] ), QgsPoint( mX[i + 2], mY[i + 2] ) );
153 : 9 : bbox.combineExtentWith( segmentBox );
154 : : }
155 : 26 : }
156 : :
157 : 27 : if ( nPoints > 0 && nPoints % 2 == 0 )
158 : : {
159 : 8 : if ( nPoints == 2 )
160 : : {
161 : 8 : bbox.combineExtentWith( mX[ 0 ], mY[ 0 ] );
162 : 8 : }
163 : 8 : bbox.combineExtentWith( mX[ nPoints - 1 ], mY[ nPoints - 1 ] );
164 : 8 : }
165 : 27 : return bbox;
166 : 0 : }
167 : :
168 : 26 : QgsRectangle QgsCircularString::segmentBoundingBox( const QgsPoint &pt1, const QgsPoint &pt2, const QgsPoint &pt3 )
169 : : {
170 : : double centerX, centerY, radius;
171 : 26 : QgsGeometryUtils::circleCenterRadius( pt1, pt2, pt3, radius, centerX, centerY );
172 : :
173 : 26 : double p1Angle = QgsGeometryUtils::ccwAngle( pt1.y() - centerY, pt1.x() - centerX );
174 : 26 : double p2Angle = QgsGeometryUtils::ccwAngle( pt2.y() - centerY, pt2.x() - centerX );
175 : 26 : double p3Angle = QgsGeometryUtils::ccwAngle( pt3.y() - centerY, pt3.x() - centerX );
176 : :
177 : : //start point, end point and compass points in between can be on bounding box
178 : 26 : QgsRectangle bbox( pt1.x(), pt1.y(), pt1.x(), pt1.y() );
179 : 26 : bbox.combineExtentWith( pt3.x(), pt3.y() );
180 : :
181 : 26 : QgsPointSequence compassPoints = compassPointsOnSegment( p1Angle, p2Angle, p3Angle, centerX, centerY, radius );
182 : 26 : QgsPointSequence::const_iterator cpIt = compassPoints.constBegin();
183 : 93 : for ( ; cpIt != compassPoints.constEnd(); ++cpIt )
184 : : {
185 : 67 : bbox.combineExtentWith( cpIt->x(), cpIt->y() );
186 : 67 : }
187 : : return bbox;
188 : 26 : }
189 : :
190 : 26 : QgsPointSequence QgsCircularString::compassPointsOnSegment( double p1Angle, double p2Angle, double p3Angle, double centerX, double centerY, double radius )
191 : : {
192 : 26 : QgsPointSequence pointList;
193 : :
194 : 26 : QgsPoint nPoint( centerX, centerY + radius );
195 : 26 : QgsPoint ePoint( centerX + radius, centerY );
196 : 26 : QgsPoint sPoint( centerX, centerY - radius );
197 : 26 : QgsPoint wPoint( centerX - radius, centerY );
198 : :
199 : 26 : if ( p3Angle >= p1Angle )
200 : : {
201 : 14 : if ( p2Angle > p1Angle && p2Angle < p3Angle )
202 : : {
203 : 2 : if ( p1Angle <= 90 && p3Angle >= 90 )
204 : : {
205 : 0 : pointList.append( nPoint );
206 : 0 : }
207 : 2 : if ( p1Angle <= 180 && p3Angle >= 180 )
208 : : {
209 : 2 : pointList.append( wPoint );
210 : 2 : }
211 : 2 : if ( p1Angle <= 270 && p3Angle >= 270 )
212 : : {
213 : 0 : pointList.append( sPoint );
214 : 0 : }
215 : 2 : }
216 : : else
217 : : {
218 : 12 : pointList.append( ePoint );
219 : 12 : if ( p1Angle >= 90 || p3Angle <= 90 )
220 : : {
221 : 10 : pointList.append( nPoint );
222 : 10 : }
223 : 12 : if ( p1Angle >= 180 || p3Angle <= 180 )
224 : : {
225 : 6 : pointList.append( wPoint );
226 : 6 : }
227 : 12 : if ( p1Angle >= 270 || p3Angle <= 270 )
228 : : {
229 : 8 : pointList.append( sPoint );
230 : 8 : }
231 : : }
232 : 14 : }
233 : : else
234 : : {
235 : 12 : if ( p2Angle < p1Angle && p2Angle > p3Angle )
236 : : {
237 : 8 : if ( p1Angle >= 270 && p3Angle <= 270 )
238 : : {
239 : 5 : pointList.append( sPoint );
240 : 5 : }
241 : 8 : if ( p1Angle >= 180 && p3Angle <= 180 )
242 : : {
243 : 8 : pointList.append( wPoint );
244 : 8 : }
245 : 8 : if ( p1Angle >= 90 && p3Angle <= 90 )
246 : : {
247 : 6 : pointList.append( nPoint );
248 : 6 : }
249 : 8 : }
250 : : else
251 : : {
252 : 4 : pointList.append( ePoint );
253 : 4 : if ( p1Angle <= 270 || p3Angle >= 270 )
254 : : {
255 : 2 : pointList.append( sPoint );
256 : 2 : }
257 : 4 : if ( p1Angle <= 180 || p3Angle >= 180 )
258 : : {
259 : 2 : pointList.append( wPoint );
260 : 2 : }
261 : 4 : if ( p1Angle <= 90 || p3Angle >= 90 )
262 : : {
263 : 2 : pointList.append( nPoint );
264 : 2 : }
265 : : }
266 : : }
267 : 26 : return pointList;
268 : 26 : }
269 : :
270 : 33 : bool QgsCircularString::fromWkb( QgsConstWkbPtr &wkbPtr )
271 : : {
272 : 33 : if ( !wkbPtr )
273 : 1 : return false;
274 : :
275 : 32 : QgsWkbTypes::Type type = wkbPtr.readHeader();
276 : 32 : if ( QgsWkbTypes::flatType( type ) != QgsWkbTypes::CircularString )
277 : : {
278 : 1 : return false;
279 : : }
280 : 31 : clearCache();
281 : 31 : mWkbType = type;
282 : :
283 : : //type
284 : 31 : bool hasZ = is3D();
285 : 31 : bool hasM = isMeasure();
286 : 31 : int nVertices = 0;
287 : 31 : wkbPtr >> nVertices;
288 : 31 : mX.resize( nVertices );
289 : 31 : mY.resize( nVertices );
290 : 31 : hasZ ? mZ.resize( nVertices ) : mZ.clear();
291 : 31 : hasM ? mM.resize( nVertices ) : mM.clear();
292 : 148 : for ( int i = 0; i < nVertices; ++i )
293 : : {
294 : 117 : wkbPtr >> mX[i];
295 : 117 : wkbPtr >> mY[i];
296 : 117 : if ( hasZ )
297 : : {
298 : 67 : wkbPtr >> mZ[i];
299 : 67 : }
300 : 117 : if ( hasM )
301 : : {
302 : 52 : wkbPtr >> mM[i];
303 : 52 : }
304 : 117 : }
305 : :
306 : 31 : return true;
307 : 33 : }
308 : :
309 : 451 : bool QgsCircularString::fromWkt( const QString &wkt )
310 : : {
311 : 451 : clear();
312 : :
313 : 451 : QPair<QgsWkbTypes::Type, QString> parts = QgsGeometryUtils::wktReadBlock( wkt );
314 : :
315 : 451 : if ( QgsWkbTypes::flatType( parts.first ) != QgsWkbTypes::CircularString )
316 : 1 : return false;
317 : 450 : mWkbType = parts.first;
318 : :
319 : 450 : parts.second = parts.second.remove( '(' ).remove( ')' );
320 : 450 : QString secondWithoutParentheses = parts.second;
321 : 450 : secondWithoutParentheses = secondWithoutParentheses.simplified().remove( ' ' );
322 : 500 : if ( ( parts.second.compare( QLatin1String( "EMPTY" ), Qt::CaseInsensitive ) == 0 ) ||
323 : 50 : secondWithoutParentheses.isEmpty() )
324 : 400 : return true;
325 : :
326 : 50 : QgsPointSequence points = QgsGeometryUtils::pointsFromWKT( parts.second, is3D(), isMeasure() );
327 : 50 : if ( points.isEmpty() )
328 : 2 : return false;
329 : :
330 : 48 : setPoints( points );
331 : 48 : return true;
332 : 451 : }
333 : :
334 : 76 : int QgsCircularString::wkbSize( QgsAbstractGeometry::WkbFlags ) const
335 : : {
336 : 76 : int binarySize = sizeof( char ) + sizeof( quint32 ) + sizeof( quint32 );
337 : 76 : binarySize += numPoints() * ( 2 + is3D() + isMeasure() ) * sizeof( double );
338 : 76 : return binarySize;
339 : : }
340 : :
341 : 31 : QByteArray QgsCircularString::asWkb( WkbFlags flags ) const
342 : : {
343 : 31 : QByteArray wkbArray;
344 : 31 : wkbArray.resize( QgsCircularString::wkbSize( flags ) );
345 : 31 : QgsWkbPtr wkb( wkbArray );
346 : 31 : wkb << static_cast<char>( QgsApplication::endian() );
347 : 31 : wkb << static_cast<quint32>( wkbType() );
348 : 31 : QgsPointSequence pts;
349 : 31 : points( pts );
350 : 31 : QgsGeometryUtils::pointsToWKB( wkb, pts, is3D(), isMeasure() );
351 : 31 : return wkbArray;
352 : 31 : }
353 : :
354 : 498 : QString QgsCircularString::asWkt( int precision ) const
355 : : {
356 : 498 : QString wkt = wktTypeStr() + ' ';
357 : :
358 : 498 : if ( isEmpty() )
359 : 400 : wkt += QLatin1String( "EMPTY" );
360 : : else
361 : : {
362 : 98 : QgsPointSequence pts;
363 : 98 : points( pts );
364 : 98 : wkt += QgsGeometryUtils::pointsToWKT( pts, precision, is3D(), isMeasure() );
365 : 98 : }
366 : 498 : return wkt;
367 : 498 : }
368 : :
369 : 4 : QDomElement QgsCircularString::asGml2( QDomDocument &doc, int precision, const QString &ns, const AxisOrder axisOrder ) const
370 : : {
371 : : // GML2 does not support curves
372 : 4 : std::unique_ptr< QgsLineString > line( curveToLine() );
373 : 4 : QDomElement gml = line->asGml2( doc, precision, ns, axisOrder );
374 : 4 : return gml;
375 : 4 : }
376 : :
377 : 15 : QDomElement QgsCircularString::asGml3( QDomDocument &doc, int precision, const QString &ns, const QgsAbstractGeometry::AxisOrder axisOrder ) const
378 : : {
379 : 15 : QgsPointSequence pts;
380 : 15 : points( pts );
381 : :
382 : 30 : QDomElement elemCurve = doc.createElementNS( ns, QStringLiteral( "Curve" ) );
383 : :
384 : 15 : if ( isEmpty() )
385 : 1 : return elemCurve;
386 : :
387 : 28 : QDomElement elemSegments = doc.createElementNS( ns, QStringLiteral( "segments" ) );
388 : 28 : QDomElement elemArcString = doc.createElementNS( ns, QStringLiteral( "ArcString" ) );
389 : 14 : elemArcString.appendChild( QgsGeometryUtils::pointsToGML3( pts, doc, precision, ns, is3D(), axisOrder ) );
390 : 14 : elemSegments.appendChild( elemArcString );
391 : 14 : elemCurve.appendChild( elemSegments );
392 : 14 : return elemCurve;
393 : 15 : }
394 : :
395 : :
396 : 3 : json QgsCircularString::asJsonObject( int precision ) const
397 : : {
398 : : // GeoJSON does not support curves
399 : 3 : std::unique_ptr< QgsLineString > line( curveToLine() );
400 : 3 : return line->asJsonObject( precision );
401 : 3 : }
402 : :
403 : 612 : bool QgsCircularString::isEmpty() const
404 : : {
405 : 612 : return mX.isEmpty();
406 : : }
407 : :
408 : 5 : bool QgsCircularString::isValid( QString &error, int flags ) const
409 : : {
410 : 5 : if ( !isEmpty() && ( numPoints() < 3 ) )
411 : : {
412 : 1 : error = QObject::tr( "CircularString has less than 3 points and is not empty." );
413 : 1 : return false;
414 : : }
415 : 4 : return QgsCurve::isValid( error, flags );
416 : 5 : }
417 : :
418 : : //curve interface
419 : 28 : double QgsCircularString::length() const
420 : : {
421 : 28 : int nPoints = numPoints();
422 : 28 : double length = 0;
423 : 59 : for ( int i = 0; i < ( nPoints - 2 ) ; i += 2 )
424 : : {
425 : 31 : length += QgsGeometryUtils::circleLength( mX[i], mY[i], mX[i + 1], mY[i + 1], mX[i + 2], mY[i + 2] );
426 : 31 : }
427 : 28 : return length;
428 : : }
429 : :
430 : 50 : QgsPoint QgsCircularString::startPoint() const
431 : : {
432 : 50 : if ( numPoints() < 1 )
433 : : {
434 : 1 : return QgsPoint();
435 : : }
436 : 49 : return pointN( 0 );
437 : 50 : }
438 : :
439 : 48 : QgsPoint QgsCircularString::endPoint() const
440 : : {
441 : 48 : if ( numPoints() < 1 )
442 : : {
443 : 1 : return QgsPoint();
444 : : }
445 : 47 : return pointN( numPoints() - 1 );
446 : 48 : }
447 : :
448 : 69 : QgsLineString *QgsCircularString::curveToLine( double tolerance, SegmentationToleranceType toleranceType ) const
449 : : {
450 : 69 : QgsLineString *line = new QgsLineString();
451 : 69 : QgsPointSequence points;
452 : 69 : int nPoints = numPoints();
453 : :
454 : 162 : for ( int i = 0; i < ( nPoints - 2 ) ; i += 2 )
455 : : {
456 : 93 : QgsGeometryUtils::segmentizeArc( pointN( i ), pointN( i + 1 ), pointN( i + 2 ), points, tolerance, toleranceType, is3D(), isMeasure() );
457 : 93 : }
458 : :
459 : 69 : line->setPoints( points );
460 : 69 : return line;
461 : 69 : }
462 : :
463 : 6 : QgsCircularString *QgsCircularString::snappedToGrid( double hSpacing, double vSpacing, double dSpacing, double mSpacing ) const
464 : : {
465 : : // prepare result
466 : 6 : std::unique_ptr<QgsCircularString> result { createEmptyWithSameType() };
467 : :
468 : 12 : bool res = snapToGridPrivate( hSpacing, vSpacing, dSpacing, mSpacing, mX, mY, mZ, mM,
469 : 6 : result->mX, result->mY, result->mZ, result->mM );
470 : 6 : if ( res )
471 : 6 : return result.release();
472 : : else
473 : 0 : return nullptr;
474 : 6 : }
475 : :
476 : 21 : bool QgsCircularString::removeDuplicateNodes( double epsilon, bool useZValues )
477 : : {
478 : 21 : if ( mX.count() <= 3 )
479 : 9 : return false; // don't create degenerate lines
480 : 12 : bool result = false;
481 : 12 : double prevX = mX.at( 0 );
482 : 12 : double prevY = mY.at( 0 );
483 : 12 : bool hasZ = is3D();
484 : 12 : bool useZ = hasZ && useZValues;
485 : 12 : double prevZ = useZ ? mZ.at( 0 ) : 0;
486 : 12 : int i = 1;
487 : 12 : int remaining = mX.count();
488 : : // we have to consider points in pairs, since a segment can validly have the same start and
489 : : // end if it has a different curve point
490 : 32 : while ( i + 1 < remaining )
491 : : {
492 : 20 : double currentCurveX = mX.at( i );
493 : 20 : double currentCurveY = mY.at( i );
494 : 20 : double currentX = mX.at( i + 1 );
495 : 20 : double currentY = mY.at( i + 1 );
496 : 20 : double currentZ = useZ ? mZ.at( i + 1 ) : 0;
497 : 21 : if ( qgsDoubleNear( currentCurveX, prevX, epsilon ) &&
498 : 12 : qgsDoubleNear( currentCurveY, prevY, epsilon ) &&
499 : 4 : qgsDoubleNear( currentX, prevX, epsilon ) &&
500 : 4 : qgsDoubleNear( currentY, prevY, epsilon ) &&
501 : 4 : ( !useZ || qgsDoubleNear( currentZ, prevZ, epsilon ) ) )
502 : : {
503 : 3 : result = true;
504 : : // remove point
505 : 3 : mX.removeAt( i );
506 : 3 : mX.removeAt( i );
507 : 3 : mY.removeAt( i );
508 : 3 : mY.removeAt( i );
509 : 3 : if ( hasZ )
510 : : {
511 : 1 : mZ.removeAt( i );
512 : 1 : mZ.removeAt( i );
513 : 1 : }
514 : 3 : remaining -= 2;
515 : 3 : }
516 : : else
517 : : {
518 : 17 : prevX = currentX;
519 : 17 : prevY = currentY;
520 : 17 : prevZ = currentZ;
521 : 17 : i += 2;
522 : : }
523 : : }
524 : 12 : return result;
525 : 21 : }
526 : :
527 : 1801 : int QgsCircularString::numPoints() const
528 : : {
529 : 1801 : return std::min( mX.size(), mY.size() );
530 : : }
531 : :
532 : 1635 : QgsPoint QgsCircularString::pointN( int i ) const
533 : : {
534 : 1635 : if ( i < 0 || std::min( mX.size(), mY.size() ) <= i )
535 : : {
536 : 2 : return QgsPoint();
537 : : }
538 : :
539 : 1633 : double x = mX.at( i );
540 : 1633 : double y = mY.at( i );
541 : 1633 : double z = 0;
542 : 1633 : double m = 0;
543 : :
544 : 1633 : if ( is3D() )
545 : : {
546 : 669 : z = mZ.at( i );
547 : 669 : }
548 : 1633 : if ( isMeasure() )
549 : : {
550 : 635 : m = mM.at( i );
551 : 635 : }
552 : :
553 : 1633 : QgsWkbTypes::Type t = QgsWkbTypes::Point;
554 : 1633 : if ( is3D() && isMeasure() )
555 : : {
556 : 522 : t = QgsWkbTypes::PointZM;
557 : 522 : }
558 : 1111 : else if ( is3D() )
559 : : {
560 : 147 : t = QgsWkbTypes::PointZ;
561 : 147 : }
562 : 964 : else if ( isMeasure() )
563 : : {
564 : 113 : t = QgsWkbTypes::PointM;
565 : 113 : }
566 : 1633 : return QgsPoint( t, x, y, z, m );
567 : 1635 : }
568 : :
569 : 64 : double QgsCircularString::xAt( int index ) const
570 : : {
571 : 64 : if ( index >= 0 && index < mX.size() )
572 : 60 : return mX.at( index );
573 : : else
574 : 4 : return 0.0;
575 : 64 : }
576 : :
577 : 64 : double QgsCircularString::yAt( int index ) const
578 : : {
579 : 64 : if ( index >= 0 && index < mY.size() )
580 : 60 : return mY.at( index );
581 : : else
582 : 4 : return 0.0;
583 : 64 : }
584 : :
585 : 6 : bool QgsCircularString::transform( QgsAbstractGeometryTransformer *transformer, QgsFeedback *feedback )
586 : : {
587 : 6 : if ( !transformer )
588 : 0 : return false;
589 : :
590 : 6 : bool hasZ = is3D();
591 : 6 : bool hasM = isMeasure();
592 : 6 : int size = mX.size();
593 : :
594 : 6 : double *srcX = mX.data();
595 : 6 : double *srcY = mY.data();
596 : 6 : double *srcM = hasM ? mM.data() : nullptr;
597 : 6 : double *srcZ = hasZ ? mZ.data() : nullptr;
598 : :
599 : 6 : bool res = true;
600 : 19 : for ( int i = 0; i < size; ++i )
601 : : {
602 : 15 : double x = *srcX;
603 : 15 : double y = *srcY;
604 : 15 : double z = hasZ ? *srcZ : std::numeric_limits<double>::quiet_NaN();
605 : 15 : double m = hasM ? *srcM : std::numeric_limits<double>::quiet_NaN();
606 : 15 : if ( !transformer->transformPoint( x, y, z, m ) )
607 : : {
608 : 2 : res = false;
609 : 2 : break;
610 : : }
611 : :
612 : 13 : *srcX++ = x;
613 : 13 : *srcY++ = y;
614 : 13 : if ( hasM )
615 : 13 : *srcM++ = m;
616 : 13 : if ( hasZ )
617 : 13 : *srcZ++ = z;
618 : :
619 : 13 : if ( feedback && feedback->isCanceled() )
620 : : {
621 : 0 : res = false;
622 : 0 : break;
623 : : }
624 : 13 : }
625 : 6 : clearCache();
626 : 6 : return res;
627 : 6 : }
628 : :
629 : 4 : void QgsCircularString::filterVertices( const std::function<bool ( const QgsPoint & )> &filter )
630 : : {
631 : 4 : bool hasZ = is3D();
632 : 4 : bool hasM = isMeasure();
633 : 4 : int size = mX.size();
634 : :
635 : 4 : double *srcX = mX.data(); // clazy:exclude=detaching-member
636 : 4 : double *srcY = mY.data(); // clazy:exclude=detaching-member
637 : 4 : double *srcM = hasM ? mM.data() : nullptr; // clazy:exclude=detaching-member
638 : 4 : double *srcZ = hasZ ? mZ.data() : nullptr; // clazy:exclude=detaching-member
639 : :
640 : 4 : double *destX = srcX;
641 : 4 : double *destY = srcY;
642 : 4 : double *destM = srcM;
643 : 4 : double *destZ = srcZ;
644 : :
645 : 4 : int filteredPoints = 0;
646 : 15 : for ( int i = 0; i < size; ++i )
647 : : {
648 : 11 : double x = *srcX++;
649 : 11 : double y = *srcY++;
650 : 11 : double z = hasZ ? *srcZ++ : std::numeric_limits<double>::quiet_NaN();
651 : 11 : double m = hasM ? *srcM++ : std::numeric_limits<double>::quiet_NaN();
652 : :
653 : 11 : if ( filter( QgsPoint( x, y, z, m ) ) )
654 : : {
655 : 8 : filteredPoints++;
656 : 8 : *destX++ = x;
657 : 8 : *destY++ = y;
658 : 8 : if ( hasM )
659 : 8 : *destM++ = m;
660 : 8 : if ( hasZ )
661 : 8 : *destZ++ = z;
662 : 8 : }
663 : 11 : }
664 : :
665 : 4 : mX.resize( filteredPoints );
666 : 4 : mY.resize( filteredPoints );
667 : 4 : if ( hasZ )
668 : 3 : mZ.resize( filteredPoints );
669 : 4 : if ( hasM )
670 : 3 : mM.resize( filteredPoints );
671 : :
672 : 4 : clearCache();
673 : 4 : }
674 : :
675 : 4 : void QgsCircularString::transformVertices( const std::function<QgsPoint( const QgsPoint & )> &transform )
676 : : {
677 : 4 : bool hasZ = is3D();
678 : 4 : bool hasM = isMeasure();
679 : 4 : int size = mX.size();
680 : :
681 : 4 : double *srcX = mX.data();
682 : 4 : double *srcY = mY.data();
683 : 4 : double *srcM = hasM ? mM.data() : nullptr;
684 : 4 : double *srcZ = hasZ ? mZ.data() : nullptr;
685 : :
686 : 17 : for ( int i = 0; i < size; ++i )
687 : : {
688 : 13 : double x = *srcX;
689 : 13 : double y = *srcY;
690 : 13 : double z = hasZ ? *srcZ : std::numeric_limits<double>::quiet_NaN();
691 : 13 : double m = hasM ? *srcM : std::numeric_limits<double>::quiet_NaN();
692 : 13 : QgsPoint res = transform( QgsPoint( x, y, z, m ) );
693 : 13 : *srcX++ = res.x();
694 : 13 : *srcY++ = res.y();
695 : 13 : if ( hasM )
696 : 13 : *srcM++ = res.m();
697 : 13 : if ( hasZ )
698 : 13 : *srcZ++ = res.z();
699 : 13 : }
700 : 4 : clearCache();
701 : 4 : }
702 : :
703 : 244 : void QgsCircularString::points( QgsPointSequence &pts ) const
704 : : {
705 : 244 : pts.clear();
706 : 244 : int nPts = numPoints();
707 : 993 : for ( int i = 0; i < nPts; ++i )
708 : : {
709 : 749 : pts.push_back( pointN( i ) );
710 : 749 : }
711 : 244 : }
712 : :
713 : 466 : void QgsCircularString::setPoints( const QgsPointSequence &points )
714 : : {
715 : 466 : clearCache();
716 : :
717 : 466 : if ( points.empty() )
718 : : {
719 : 2 : mWkbType = QgsWkbTypes::CircularString;
720 : 2 : mX.clear();
721 : 2 : mY.clear();
722 : 2 : mZ.clear();
723 : 2 : mM.clear();
724 : 2 : return;
725 : : }
726 : :
727 : : //get wkb type from first point
728 : 464 : const QgsPoint &firstPt = points.at( 0 );
729 : 464 : bool hasZ = firstPt.is3D();
730 : 464 : bool hasM = firstPt.isMeasure();
731 : :
732 : 464 : setZMTypeFromSubGeometry( &firstPt, QgsWkbTypes::CircularString );
733 : :
734 : 464 : mX.resize( points.size() );
735 : 464 : mY.resize( points.size() );
736 : 464 : if ( hasZ )
737 : : {
738 : 188 : mZ.resize( points.size() );
739 : 188 : }
740 : : else
741 : : {
742 : 276 : mZ.clear();
743 : : }
744 : 464 : if ( hasM )
745 : : {
746 : 168 : mM.resize( points.size() );
747 : 168 : }
748 : : else
749 : : {
750 : 296 : mM.clear();
751 : : }
752 : :
753 : 1980 : for ( int i = 0; i < points.size(); ++i )
754 : : {
755 : 1516 : mX[i] = points[i].x();
756 : 1516 : mY[i] = points[i].y();
757 : 1516 : if ( hasZ )
758 : : {
759 : 632 : double z = points.at( i ).z();
760 : 632 : mZ[i] = std::isnan( z ) ? 0 : z;
761 : 632 : }
762 : 1516 : if ( hasM )
763 : : {
764 : 552 : double m = points.at( i ).m();
765 : 552 : mM[i] = std::isnan( m ) ? 0 : m;
766 : 552 : }
767 : 1516 : }
768 : 466 : }
769 : :
770 : 0 : void QgsCircularString::draw( QPainter &p ) const
771 : : {
772 : 0 : QPainterPath path;
773 : 0 : addToPainterPath( path );
774 : 0 : p.drawPath( path );
775 : 0 : }
776 : :
777 : 6 : void QgsCircularString::transform( const QgsCoordinateTransform &ct, QgsCoordinateTransform::TransformDirection d, bool transformZ )
778 : : {
779 : 6 : clearCache();
780 : :
781 : 6 : double *zArray = mZ.data();
782 : :
783 : 6 : bool hasZ = is3D();
784 : 6 : int nPoints = numPoints();
785 : 6 : bool useDummyZ = !hasZ || !transformZ;
786 : 6 : if ( useDummyZ )
787 : : {
788 : 6 : zArray = new double[nPoints];
789 : 18 : for ( int i = 0; i < nPoints; ++i )
790 : : {
791 : 12 : zArray[i] = 0;
792 : 12 : }
793 : 6 : }
794 : 6 : ct.transformCoords( nPoints, mX.data(), mY.data(), zArray, d );
795 : 6 : if ( useDummyZ )
796 : : {
797 : 6 : delete[] zArray;
798 : 6 : }
799 : 6 : }
800 : :
801 : 3 : void QgsCircularString::transform( const QTransform &t, double zTranslate, double zScale, double mTranslate, double mScale )
802 : : {
803 : 3 : clearCache();
804 : :
805 : 3 : int nPoints = numPoints();
806 : 3 : bool hasZ = is3D();
807 : 3 : bool hasM = isMeasure();
808 : 9 : for ( int i = 0; i < nPoints; ++i )
809 : : {
810 : : qreal x, y;
811 : 6 : t.map( mX.at( i ), mY.at( i ), &x, &y );
812 : 6 : mX[i] = x;
813 : 6 : mY[i] = y;
814 : 6 : if ( hasZ )
815 : : {
816 : 6 : mZ[i] = mZ.at( i ) * zScale + zTranslate;
817 : 6 : }
818 : 6 : if ( hasM )
819 : : {
820 : 6 : mM[i] = mM.at( i ) * mScale + mTranslate;
821 : 6 : }
822 : 6 : }
823 : 3 : }
824 : :
825 : 6 : void QgsCircularString::addToPainterPath( QPainterPath &path ) const
826 : : {
827 : 6 : int nPoints = numPoints();
828 : 6 : if ( nPoints < 1 )
829 : : {
830 : 1 : return;
831 : : }
832 : :
833 : 5 : if ( path.isEmpty() || path.currentPosition() != QPointF( mX[0], mY[0] ) )
834 : : {
835 : 5 : path.moveTo( QPointF( mX[0], mY[0] ) );
836 : 5 : }
837 : :
838 : 8 : for ( int i = 0; i < ( nPoints - 2 ) ; i += 2 )
839 : : {
840 : 3 : QgsPointSequence pt;
841 : 3 : QgsGeometryUtils::segmentizeArc( QgsPoint( mX[i], mY[i] ), QgsPoint( mX[i + 1], mY[i + 1] ), QgsPoint( mX[i + 2], mY[i + 2] ), pt );
842 : 543 : for ( int j = 1; j < pt.size(); ++j )
843 : : {
844 : 540 : path.lineTo( pt.at( j ).x(), pt.at( j ).y() );
845 : 540 : }
846 : : #if 0
847 : : //arcTo( path, QPointF( mX[i], mY[i] ), QPointF( mX[i + 1], mY[i + 1] ), QPointF( mX[i + 2], mY[i + 2] ) );
848 : : #endif
849 : 3 : }
850 : :
851 : : //if number of points is even, connect to last point with straight line (even though the circular string is not valid)
852 : 5 : if ( nPoints % 2 == 0 )
853 : : {
854 : 2 : path.lineTo( mX[ nPoints - 1 ], mY[ nPoints - 1 ] );
855 : 2 : }
856 : 6 : }
857 : :
858 : : #if 0
859 : : void QgsCircularString::arcTo( QPainterPath &path, QPointF pt1, QPointF pt2, QPointF pt3 )
860 : : {
861 : : double centerX, centerY, radius;
862 : : QgsGeometryUtils::circleCenterRadius( QgsPoint( pt1.x(), pt1.y() ), QgsPoint( pt2.x(), pt2.y() ), QgsPoint( pt3.x(), pt3.y() ),
863 : : radius, centerX, centerY );
864 : :
865 : : double p1Angle = QgsGeometryUtils::ccwAngle( pt1.y() - centerY, pt1.x() - centerX );
866 : : double sweepAngle = QgsGeometryUtils::sweepAngle( centerX, centerY, pt1.x(), pt1.y(), pt2.x(), pt2.y(), pt3.x(), pt3.y() );
867 : :
868 : : double diameter = 2 * radius;
869 : : path.arcTo( centerX - radius, centerY - radius, diameter, diameter, p1Angle, sweepAngle );
870 : : }
871 : : #endif
872 : :
873 : 0 : void QgsCircularString::drawAsPolygon( QPainter &p ) const
874 : : {
875 : 0 : draw( p );
876 : 0 : }
877 : :
878 : 20 : bool QgsCircularString::insertVertex( QgsVertexId position, const QgsPoint &vertex )
879 : : {
880 : 20 : if ( position.vertex >= mX.size() || position.vertex < 1 )
881 : : {
882 : 6 : return false;
883 : : }
884 : :
885 : 14 : mX.insert( position.vertex, vertex.x() );
886 : 14 : mY.insert( position.vertex, vertex.y() );
887 : 14 : if ( is3D() )
888 : : {
889 : 4 : mZ.insert( position.vertex, vertex.z() );
890 : 4 : }
891 : 14 : if ( isMeasure() )
892 : : {
893 : 4 : mM.insert( position.vertex, vertex.m() );
894 : 4 : }
895 : :
896 : 14 : bool vertexNrEven = ( position.vertex % 2 == 0 );
897 : 14 : if ( vertexNrEven )
898 : : {
899 : 2 : insertVertexBetween( position.vertex - 2, position.vertex - 1, position.vertex );
900 : 2 : }
901 : : else
902 : : {
903 : 12 : insertVertexBetween( position.vertex, position.vertex + 1, position.vertex - 1 );
904 : : }
905 : 14 : clearCache(); //set bounding box invalid
906 : 14 : return true;
907 : 20 : }
908 : :
909 : 38 : bool QgsCircularString::moveVertex( QgsVertexId position, const QgsPoint &newPos )
910 : : {
911 : 38 : if ( position.vertex < 0 || position.vertex >= mX.size() )
912 : : {
913 : 5 : return false;
914 : : }
915 : :
916 : 33 : mX[position.vertex] = newPos.x();
917 : 33 : mY[position.vertex] = newPos.y();
918 : 33 : if ( is3D() && newPos.is3D() )
919 : : {
920 : 8 : mZ[position.vertex] = newPos.z();
921 : 8 : }
922 : 33 : if ( isMeasure() && newPos.isMeasure() )
923 : : {
924 : 6 : mM[position.vertex] = newPos.m();
925 : 6 : }
926 : 33 : clearCache(); //set bounding box invalid
927 : 33 : return true;
928 : 38 : }
929 : :
930 : 12 : bool QgsCircularString::deleteVertex( QgsVertexId position )
931 : : {
932 : 12 : int nVertices = this->numPoints();
933 : 12 : if ( nVertices < 4 ) //circular string must have at least 3 vertices
934 : : {
935 : 6 : clear();
936 : 6 : return true;
937 : : }
938 : 6 : if ( position.vertex < 0 || position.vertex > ( nVertices - 1 ) )
939 : : {
940 : 2 : return false;
941 : : }
942 : :
943 : 4 : if ( position.vertex < ( nVertices - 2 ) )
944 : : {
945 : : //remove this and the following vertex
946 : 4 : deleteVertex( position.vertex + 1 );
947 : 4 : deleteVertex( position.vertex );
948 : 4 : }
949 : : else //remove this and the preceding vertex
950 : : {
951 : 0 : deleteVertex( position.vertex );
952 : 0 : deleteVertex( position.vertex - 1 );
953 : : }
954 : :
955 : 4 : clearCache(); //set bounding box invalid
956 : 4 : return true;
957 : 12 : }
958 : :
959 : 8 : void QgsCircularString::deleteVertex( int i )
960 : : {
961 : 8 : mX.remove( i );
962 : 8 : mY.remove( i );
963 : 8 : if ( is3D() )
964 : : {
965 : 4 : mZ.remove( i );
966 : 4 : }
967 : 8 : if ( isMeasure() )
968 : : {
969 : 4 : mM.remove( i );
970 : 4 : }
971 : 8 : clearCache();
972 : 8 : }
973 : :
974 : 46 : double QgsCircularString::closestSegment( const QgsPoint &pt, QgsPoint &segmentPt, QgsVertexId &vertexAfter, int *leftOf, double epsilon ) const
975 : : {
976 : 46 : double minDist = std::numeric_limits<double>::max();
977 : 46 : QgsPoint minDistSegmentPoint;
978 : 46 : QgsVertexId minDistVertexAfter;
979 : 46 : int minDistLeftOf = 0;
980 : :
981 : 46 : double currentDist = 0.0;
982 : :
983 : 46 : int nPoints = numPoints();
984 : 89 : for ( int i = 0; i < ( nPoints - 2 ) ; i += 2 )
985 : : {
986 : 43 : currentDist = closestPointOnArc( mX[i], mY[i], mX[i + 1], mY[i + 1], mX[i + 2], mY[i + 2], pt, segmentPt, vertexAfter, leftOf, epsilon );
987 : 43 : if ( currentDist < minDist )
988 : : {
989 : 43 : minDist = currentDist;
990 : 43 : minDistSegmentPoint = segmentPt;
991 : 43 : minDistVertexAfter.vertex = vertexAfter.vertex + i;
992 : 43 : if ( leftOf )
993 : : {
994 : 43 : minDistLeftOf = *leftOf;
995 : 43 : }
996 : 43 : }
997 : 43 : }
998 : :
999 : 46 : if ( minDist == std::numeric_limits<double>::max() )
1000 : 3 : return -1; // error: no segments
1001 : :
1002 : 43 : segmentPt = minDistSegmentPoint;
1003 : 43 : vertexAfter = minDistVertexAfter;
1004 : 43 : vertexAfter.part = 0;
1005 : 43 : vertexAfter.ring = 0;
1006 : 43 : if ( leftOf )
1007 : : {
1008 : 43 : *leftOf = qgsDoubleNear( minDist, 0.0 ) ? 0 : minDistLeftOf;
1009 : 43 : }
1010 : 43 : return minDist;
1011 : 46 : }
1012 : :
1013 : 228 : bool QgsCircularString::pointAt( int node, QgsPoint &point, QgsVertexId::VertexType &type ) const
1014 : : {
1015 : 228 : if ( node < 0 || node >= numPoints() )
1016 : : {
1017 : 16 : return false;
1018 : : }
1019 : 212 : point = pointN( node );
1020 : 212 : type = ( node % 2 == 0 ) ? QgsVertexId::SegmentVertex : QgsVertexId::CurveVertex;
1021 : 212 : return true;
1022 : 228 : }
1023 : :
1024 : 20 : void QgsCircularString::sumUpArea( double &sum ) const
1025 : : {
1026 : 20 : int maxIndex = numPoints() - 2;
1027 : :
1028 : 40 : for ( int i = 0; i < maxIndex; i += 2 )
1029 : : {
1030 : 20 : QgsPoint p1( mX[i], mY[i] );
1031 : 20 : QgsPoint p2( mX[i + 1], mY[i + 1] );
1032 : 20 : QgsPoint p3( mX[i + 2], mY[i + 2] );
1033 : :
1034 : : //segment is a full circle, p2 is the center point
1035 : 20 : if ( p1 == p3 )
1036 : : {
1037 : 2 : double r2 = QgsGeometryUtils::sqrDistance2D( p1, p2 ) / 4.0;
1038 : 2 : sum += M_PI * r2;
1039 : 2 : continue;
1040 : : }
1041 : :
1042 : 18 : sum += 0.5 * ( mX[i] * mY[i + 2] - mY[i] * mX[i + 2] );
1043 : :
1044 : : //calculate area between circle and chord, then sum / subtract from total area
1045 : 18 : double midPointX = ( p1.x() + p3.x() ) / 2.0;
1046 : 18 : double midPointY = ( p1.y() + p3.y() ) / 2.0;
1047 : :
1048 : : double radius, centerX, centerY;
1049 : 18 : QgsGeometryUtils::circleCenterRadius( p1, p2, p3, radius, centerX, centerY );
1050 : :
1051 : 18 : double d = std::sqrt( QgsGeometryUtils::sqrDistance2D( QgsPoint( centerX, centerY ), QgsPoint( midPointX, midPointY ) ) );
1052 : 18 : double r2 = radius * radius;
1053 : :
1054 : 18 : if ( d > radius )
1055 : : {
1056 : : //d cannot be greater than radius, something must be wrong...
1057 : 0 : continue;
1058 : : }
1059 : :
1060 : 18 : bool circlePointLeftOfLine = QgsGeometryUtils::leftOfLine( p2.x(), p2.y(), p1.x(), p1.y(), p3.x(), p3.y() ) < 0;
1061 : 18 : bool centerPointLeftOfLine = QgsGeometryUtils::leftOfLine( centerX, centerY, p1.x(), p1.y(), p3.x(), p3.y() ) < 0;
1062 : :
1063 : 18 : double cov = 0.5 - d * std::sqrt( r2 - d * d ) / ( M_PI * r2 ) - M_1_PI * std::asin( d / radius );
1064 : 18 : double circleChordArea = 0;
1065 : 18 : if ( circlePointLeftOfLine == centerPointLeftOfLine )
1066 : : {
1067 : 7 : circleChordArea = M_PI * r2 * ( 1 - cov );
1068 : 7 : }
1069 : : else
1070 : : {
1071 : 11 : circleChordArea = M_PI * r2 * cov;
1072 : : }
1073 : :
1074 : 18 : if ( !circlePointLeftOfLine )
1075 : : {
1076 : 11 : sum += circleChordArea;
1077 : 11 : }
1078 : : else
1079 : : {
1080 : 7 : sum -= circleChordArea;
1081 : : }
1082 : 20 : }
1083 : 20 : }
1084 : :
1085 : 17 : bool QgsCircularString::hasCurvedSegments() const
1086 : : {
1087 : 17 : return true;
1088 : : }
1089 : :
1090 : 43 : double QgsCircularString::closestPointOnArc( double x1, double y1, double x2, double y2, double x3, double y3,
1091 : : const QgsPoint &pt, QgsPoint &segmentPt, QgsVertexId &vertexAfter, int *leftOf, double epsilon )
1092 : : {
1093 : : double radius, centerX, centerY;
1094 : 43 : QgsPoint pt1( x1, y1 );
1095 : 43 : QgsPoint pt2( x2, y2 );
1096 : 43 : QgsPoint pt3( x3, y3 );
1097 : :
1098 : 43 : QgsGeometryUtils::circleCenterRadius( pt1, pt2, pt3, radius, centerX, centerY );
1099 : 43 : double angle = QgsGeometryUtils::ccwAngle( pt.y() - centerY, pt.x() - centerX );
1100 : 43 : double angle1 = QgsGeometryUtils::ccwAngle( pt1.y() - centerY, pt1.x() - centerX );
1101 : 43 : double angle2 = QgsGeometryUtils::ccwAngle( pt2.y() - centerY, pt2.x() - centerX );
1102 : 43 : double angle3 = QgsGeometryUtils::ccwAngle( pt3.y() - centerY, pt3.x() - centerX );
1103 : :
1104 : 43 : bool clockwise = QgsGeometryUtils::circleClockwise( angle1, angle2, angle3 );
1105 : :
1106 : 43 : if ( QgsGeometryUtils::angleOnCircle( angle, angle1, angle2, angle3 ) )
1107 : : {
1108 : : //get point on line center -> pt with distance radius
1109 : 29 : segmentPt = QgsGeometryUtils::pointOnLineWithDistance( QgsPoint( centerX, centerY ), pt, radius );
1110 : :
1111 : : //vertexAfter
1112 : 29 : vertexAfter.vertex = QgsGeometryUtils::circleAngleBetween( angle, angle1, angle2, clockwise ) ? 1 : 2;
1113 : 29 : }
1114 : : else
1115 : : {
1116 : 14 : double distPtPt1 = QgsGeometryUtils::sqrDistance2D( pt, pt1 );
1117 : 14 : double distPtPt3 = QgsGeometryUtils::sqrDistance2D( pt, pt3 );
1118 : 14 : segmentPt = ( distPtPt1 <= distPtPt3 ) ? pt1 : pt3;
1119 : 14 : vertexAfter.vertex = ( distPtPt1 <= distPtPt3 ) ? 1 : 2;
1120 : : }
1121 : :
1122 : 43 : double sqrDistance = QgsGeometryUtils::sqrDistance2D( segmentPt, pt );
1123 : : //prevent rounding errors if the point is directly on the segment
1124 : 43 : if ( qgsDoubleNear( sqrDistance, 0.0, epsilon ) )
1125 : : {
1126 : 3 : segmentPt.setX( pt.x() );
1127 : 3 : segmentPt.setY( pt.y() );
1128 : 3 : sqrDistance = 0.0;
1129 : 3 : }
1130 : :
1131 : 43 : if ( leftOf )
1132 : : {
1133 : 43 : double sqrDistancePointToCenter = ( pt.x() - centerX ) * ( pt.x() - centerX ) + ( pt.y() - centerY ) * ( pt.y() - centerY );
1134 : 43 : *leftOf = clockwise ? ( sqrDistancePointToCenter > radius * radius ? -1 : 1 )
1135 : 32 : : ( sqrDistancePointToCenter < radius * radius ? -1 : 1 );
1136 : 43 : }
1137 : :
1138 : 43 : return sqrDistance;
1139 : 43 : }
1140 : :
1141 : 14 : void QgsCircularString::insertVertexBetween( int after, int before, int pointOnCircle )
1142 : : {
1143 : 14 : double xAfter = mX.at( after );
1144 : 14 : double yAfter = mY.at( after );
1145 : 14 : double xBefore = mX.at( before );
1146 : 14 : double yBefore = mY.at( before );
1147 : 14 : double xOnCircle = mX.at( pointOnCircle );
1148 : 14 : double yOnCircle = mY.at( pointOnCircle );
1149 : :
1150 : : double radius, centerX, centerY;
1151 : 14 : QgsGeometryUtils::circleCenterRadius( QgsPoint( xAfter, yAfter ), QgsPoint( xBefore, yBefore ), QgsPoint( xOnCircle, yOnCircle ), radius, centerX, centerY );
1152 : :
1153 : 14 : double x = ( xAfter + xBefore ) / 2.0;
1154 : 14 : double y = ( yAfter + yBefore ) / 2.0;
1155 : :
1156 : 14 : QgsPoint newVertex = QgsGeometryUtils::pointOnLineWithDistance( QgsPoint( centerX, centerY ), QgsPoint( x, y ), radius );
1157 : 14 : mX.insert( before, newVertex.x() );
1158 : 14 : mY.insert( before, newVertex.y() );
1159 : :
1160 : 14 : if ( is3D() )
1161 : : {
1162 : 4 : mZ.insert( before, ( mZ[after] + mZ[before] ) / 2.0 );
1163 : 4 : }
1164 : 14 : if ( isMeasure() )
1165 : : {
1166 : 4 : mM.insert( before, ( mM[after] + mM[before] ) / 2.0 );
1167 : 4 : }
1168 : 14 : clearCache();
1169 : 14 : }
1170 : :
1171 : 87 : double QgsCircularString::vertexAngle( QgsVertexId vId ) const
1172 : : {
1173 : 87 : if ( numPoints() < 3 )
1174 : : {
1175 : : //undefined
1176 : 8 : return 0.0;
1177 : : }
1178 : :
1179 : 79 : int before = vId.vertex - 1;
1180 : 79 : int vertex = vId.vertex;
1181 : 79 : int after = vId.vertex + 1;
1182 : :
1183 : 79 : if ( vId.vertex % 2 != 0 ) // a curve vertex
1184 : : {
1185 : 28 : if ( vId.vertex >= 1 && vId.vertex < numPoints() - 1 )
1186 : : {
1187 : 56 : return QgsGeometryUtils::circleTangentDirection( QgsPoint( mX[vertex], mY[vertex] ), QgsPoint( mX[before], mY[before] ),
1188 : 28 : QgsPoint( mX[vertex], mY[vertex] ), QgsPoint( mX[after], mY[after] ) );
1189 : : }
1190 : 0 : }
1191 : : else //a point vertex
1192 : : {
1193 : 51 : if ( vId.vertex == 0 )
1194 : : {
1195 : 40 : return QgsGeometryUtils::circleTangentDirection( QgsPoint( mX[0], mY[0] ), QgsPoint( mX[0], mY[0] ),
1196 : 20 : QgsPoint( mX[1], mY[1] ), QgsPoint( mX[2], mY[2] ) );
1197 : : }
1198 : 31 : if ( vId.vertex >= numPoints() - 1 )
1199 : : {
1200 : 23 : int a = numPoints() - 3;
1201 : 23 : int b = numPoints() - 2;
1202 : 23 : int c = numPoints() - 1;
1203 : 46 : return QgsGeometryUtils::circleTangentDirection( QgsPoint( mX[c], mY[c] ), QgsPoint( mX[a], mY[a] ),
1204 : 23 : QgsPoint( mX[b], mY[b] ), QgsPoint( mX[c], mY[c] ) );
1205 : : }
1206 : : else
1207 : : {
1208 : 8 : if ( vId.vertex + 2 > numPoints() - 1 )
1209 : : {
1210 : 0 : return 0.0;
1211 : : }
1212 : :
1213 : 8 : int vertex1 = vId.vertex - 2;
1214 : 8 : int vertex2 = vId.vertex - 1;
1215 : 8 : int vertex3 = vId.vertex;
1216 : 16 : double angle1 = QgsGeometryUtils::circleTangentDirection( QgsPoint( mX[vertex3], mY[vertex3] ),
1217 : 8 : QgsPoint( mX[vertex1], mY[vertex1] ), QgsPoint( mX[vertex2], mY[vertex2] ), QgsPoint( mX[vertex3], mY[vertex3] ) );
1218 : 8 : int vertex4 = vId.vertex + 1;
1219 : 8 : int vertex5 = vId.vertex + 2;
1220 : 16 : double angle2 = QgsGeometryUtils::circleTangentDirection( QgsPoint( mX[vertex3], mY[vertex3] ),
1221 : 8 : QgsPoint( mX[vertex3], mY[vertex3] ), QgsPoint( mX[vertex4], mY[vertex4] ), QgsPoint( mX[vertex5], mY[vertex5] ) );
1222 : 8 : return QgsGeometryUtils::averageAngle( angle1, angle2 );
1223 : : }
1224 : : }
1225 : 0 : return 0.0;
1226 : 87 : }
1227 : :
1228 : 56 : double QgsCircularString::segmentLength( QgsVertexId startVertex ) const
1229 : : {
1230 : 56 : if ( startVertex.vertex < 0 || startVertex.vertex >= mX.count() - 2 )
1231 : 26 : return 0.0;
1232 : :
1233 : 30 : if ( startVertex.vertex % 2 == 1 )
1234 : 3 : return 0.0; // curve point?
1235 : :
1236 : 27 : double x1 = mX.at( startVertex.vertex );
1237 : 27 : double y1 = mY.at( startVertex.vertex );
1238 : 27 : double x2 = mX.at( startVertex.vertex + 1 );
1239 : 27 : double y2 = mY.at( startVertex.vertex + 1 );
1240 : 27 : double x3 = mX.at( startVertex.vertex + 2 );
1241 : 27 : double y3 = mY.at( startVertex.vertex + 2 );
1242 : 27 : return QgsGeometryUtils::circleLength( x1, y1, x2, y2, x3, y3 );
1243 : 56 : }
1244 : :
1245 : 5 : QgsCircularString *QgsCircularString::reversed() const
1246 : : {
1247 : 5 : QgsCircularString *copy = clone();
1248 : 5 : std::reverse( copy->mX.begin(), copy->mX.end() );
1249 : 5 : std::reverse( copy->mY.begin(), copy->mY.end() );
1250 : 5 : if ( is3D() )
1251 : : {
1252 : 4 : std::reverse( copy->mZ.begin(), copy->mZ.end() );
1253 : 4 : }
1254 : 5 : if ( isMeasure() )
1255 : : {
1256 : 4 : std::reverse( copy->mM.begin(), copy->mM.end() );
1257 : 4 : }
1258 : 5 : return copy;
1259 : : }
1260 : :
1261 : 17 : QgsPoint *QgsCircularString::interpolatePoint( const double distance ) const
1262 : : {
1263 : 17 : if ( distance < 0 )
1264 : 1 : return nullptr;
1265 : :
1266 : 16 : double distanceTraversed = 0;
1267 : 16 : const int totalPoints = numPoints();
1268 : 16 : if ( totalPoints == 0 )
1269 : 1 : return nullptr;
1270 : :
1271 : 15 : QgsWkbTypes::Type pointType = QgsWkbTypes::Point;
1272 : 15 : if ( is3D() )
1273 : 13 : pointType = QgsWkbTypes::PointZ;
1274 : 15 : if ( isMeasure() )
1275 : 13 : pointType = QgsWkbTypes::addM( pointType );
1276 : :
1277 : 15 : const double *x = mX.constData();
1278 : 15 : const double *y = mY.constData();
1279 : 15 : const double *z = is3D() ? mZ.constData() : nullptr;
1280 : 15 : const double *m = isMeasure() ? mM.constData() : nullptr;
1281 : :
1282 : 15 : double prevX = *x++;
1283 : 15 : double prevY = *y++;
1284 : 15 : double prevZ = z ? *z++ : 0.0;
1285 : 15 : double prevM = m ? *m++ : 0.0;
1286 : :
1287 : 15 : if ( qgsDoubleNear( distance, 0.0 ) )
1288 : : {
1289 : 1 : return new QgsPoint( pointType, prevX, prevY, prevZ, prevM );
1290 : : }
1291 : :
1292 : 19 : for ( int i = 0; i < ( totalPoints - 2 ) ; i += 2 )
1293 : : {
1294 : 17 : double x1 = prevX;
1295 : 17 : double y1 = prevY;
1296 : 17 : double z1 = prevZ;
1297 : 17 : double m1 = prevM;
1298 : :
1299 : 17 : double x2 = *x++;
1300 : 17 : double y2 = *y++;
1301 : 17 : double z2 = z ? *z++ : 0.0;
1302 : 17 : double m2 = m ? *m++ : 0.0;
1303 : :
1304 : 17 : double x3 = *x++;
1305 : 17 : double y3 = *y++;
1306 : 17 : double z3 = z ? *z++ : 0.0;
1307 : 17 : double m3 = m ? *m++ : 0.0;
1308 : :
1309 : 17 : const double segmentLength = QgsGeometryUtils::circleLength( x1, y1, x2, y2, x3, y3 );
1310 : 17 : if ( distance < distanceTraversed + segmentLength || qgsDoubleNear( distance, distanceTraversed + segmentLength ) )
1311 : : {
1312 : : // point falls on this segment - truncate to segment length if qgsDoubleNear test was actually > segment length
1313 : 12 : const double distanceToPoint = std::min( distance - distanceTraversed, segmentLength );
1314 : 24 : return new QgsPoint( QgsGeometryUtils::interpolatePointOnArc( QgsPoint( pointType, x1, y1, z1, m1 ),
1315 : 12 : QgsPoint( pointType, x2, y2, z2, m2 ),
1316 : 12 : QgsPoint( pointType, x3, y3, z3, m3 ), distanceToPoint ) );
1317 : : }
1318 : :
1319 : 5 : distanceTraversed += segmentLength;
1320 : :
1321 : 5 : prevX = x3;
1322 : 5 : prevY = y3;
1323 : 5 : prevZ = z3;
1324 : 5 : prevM = m3;
1325 : 5 : }
1326 : :
1327 : 2 : return nullptr;
1328 : 17 : }
1329 : :
1330 : 26 : QgsCircularString *QgsCircularString::curveSubstring( double startDistance, double endDistance ) const
1331 : : {
1332 : 26 : if ( startDistance < 0 && endDistance < 0 )
1333 : 1 : return createEmptyWithSameType();
1334 : :
1335 : 25 : endDistance = std::max( startDistance, endDistance );
1336 : :
1337 : 25 : const int totalPoints = numPoints();
1338 : 25 : if ( totalPoints == 0 )
1339 : 1 : return clone();
1340 : :
1341 : 24 : QVector< QgsPoint > substringPoints;
1342 : 24 : substringPoints.reserve( totalPoints );
1343 : :
1344 : 24 : QgsWkbTypes::Type pointType = QgsWkbTypes::Point;
1345 : 24 : if ( is3D() )
1346 : 20 : pointType = QgsWkbTypes::PointZ;
1347 : 24 : if ( isMeasure() )
1348 : 20 : pointType = QgsWkbTypes::addM( pointType );
1349 : :
1350 : 24 : const double *x = mX.constData();
1351 : 24 : const double *y = mY.constData();
1352 : 24 : const double *z = is3D() ? mZ.constData() : nullptr;
1353 : 24 : const double *m = isMeasure() ? mM.constData() : nullptr;
1354 : :
1355 : 24 : double distanceTraversed = 0;
1356 : 24 : double prevX = *x++;
1357 : 24 : double prevY = *y++;
1358 : 24 : double prevZ = z ? *z++ : 0.0;
1359 : 24 : double prevM = m ? *m++ : 0.0;
1360 : 24 : bool foundStart = false;
1361 : :
1362 : 24 : if ( startDistance < 0 )
1363 : 8 : startDistance = 0;
1364 : :
1365 : 44 : for ( int i = 0; i < ( totalPoints - 2 ) ; i += 2 )
1366 : : {
1367 : 33 : double x1 = prevX;
1368 : 33 : double y1 = prevY;
1369 : 33 : double z1 = prevZ;
1370 : 33 : double m1 = prevM;
1371 : :
1372 : 33 : double x2 = *x++;
1373 : 33 : double y2 = *y++;
1374 : 33 : double z2 = z ? *z++ : 0.0;
1375 : 33 : double m2 = m ? *m++ : 0.0;
1376 : :
1377 : 33 : double x3 = *x++;
1378 : 33 : double y3 = *y++;
1379 : 33 : double z3 = z ? *z++ : 0.0;
1380 : 33 : double m3 = m ? *m++ : 0.0;
1381 : :
1382 : 33 : bool addedSegmentEnd = false;
1383 : 33 : const double segmentLength = QgsGeometryUtils::circleLength( x1, y1, x2, y2, x3, y3 );
1384 : 33 : if ( distanceTraversed <= startDistance && startDistance < distanceTraversed + segmentLength )
1385 : : {
1386 : : // start point falls on this segment
1387 : 23 : const double distanceToStart = startDistance - distanceTraversed;
1388 : 46 : const QgsPoint startPoint = QgsGeometryUtils::interpolatePointOnArc( QgsPoint( pointType, x1, y1, z1, m1 ),
1389 : 23 : QgsPoint( pointType, x2, y2, z2, m2 ),
1390 : 23 : QgsPoint( pointType, x3, y3, z3, m3 ), distanceToStart );
1391 : :
1392 : : // does end point also fall on this segment?
1393 : 23 : const bool endPointOnSegment = distanceTraversed + segmentLength > endDistance;
1394 : 23 : if ( endPointOnSegment )
1395 : : {
1396 : 10 : const double distanceToEnd = endDistance - distanceTraversed;
1397 : 10 : const double midPointDistance = ( distanceToEnd - distanceToStart ) * 0.5 + distanceToStart;
1398 : 20 : substringPoints << startPoint
1399 : 20 : << QgsGeometryUtils::interpolatePointOnArc( QgsPoint( pointType, x1, y1, z1, m1 ),
1400 : 10 : QgsPoint( pointType, x2, y2, z2, m2 ),
1401 : 10 : QgsPoint( pointType, x3, y3, z3, m3 ), midPointDistance )
1402 : 20 : << QgsGeometryUtils::interpolatePointOnArc( QgsPoint( pointType, x1, y1, z1, m1 ),
1403 : 10 : QgsPoint( pointType, x2, y2, z2, m2 ),
1404 : 10 : QgsPoint( pointType, x3, y3, z3, m3 ), distanceToEnd );
1405 : 10 : addedSegmentEnd = true;
1406 : 10 : }
1407 : : else
1408 : : {
1409 : 13 : const double midPointDistance = ( segmentLength - distanceToStart ) * 0.5 + distanceToStart;
1410 : 26 : substringPoints << startPoint
1411 : 26 : << QgsGeometryUtils::interpolatePointOnArc( QgsPoint( pointType, x1, y1, z1, m1 ),
1412 : 13 : QgsPoint( pointType, x2, y2, z2, m2 ),
1413 : 13 : QgsPoint( pointType, x3, y3, z3, m3 ), midPointDistance )
1414 : 13 : << QgsPoint( pointType, x3, y3, z3, m3 );
1415 : 13 : addedSegmentEnd = true;
1416 : : }
1417 : 23 : foundStart = true;
1418 : 23 : }
1419 : 33 : if ( !addedSegmentEnd && foundStart && ( distanceTraversed + segmentLength > endDistance ) )
1420 : : {
1421 : : // end point falls on this segment
1422 : 2 : const double distanceToEnd = endDistance - distanceTraversed;
1423 : : // add mid point, at half way along this arc, then add the interpolated end point
1424 : 6 : substringPoints << QgsGeometryUtils::interpolatePointOnArc( QgsPoint( pointType, x1, y1, z1, m1 ),
1425 : 2 : QgsPoint( pointType, x2, y2, z2, m2 ),
1426 : 2 : QgsPoint( pointType, x3, y3, z3, m3 ), distanceToEnd / 2.0 )
1427 : :
1428 : 4 : << QgsGeometryUtils::interpolatePointOnArc( QgsPoint( pointType, x1, y1, z1, m1 ),
1429 : 2 : QgsPoint( pointType, x2, y2, z2, m2 ),
1430 : 2 : QgsPoint( pointType, x3, y3, z3, m3 ), distanceToEnd );
1431 : 2 : }
1432 : 31 : else if ( !addedSegmentEnd && foundStart )
1433 : : {
1434 : 8 : substringPoints << QgsPoint( pointType, x2, y2, z2, m2 )
1435 : 4 : << QgsPoint( pointType, x3, y3, z3, m3 );
1436 : 4 : }
1437 : :
1438 : 33 : prevX = x3;
1439 : 33 : prevY = y3;
1440 : 33 : prevZ = z3;
1441 : 33 : prevM = m3;
1442 : 33 : distanceTraversed += segmentLength;
1443 : 33 : if ( distanceTraversed >= endDistance )
1444 : 13 : break;
1445 : 20 : }
1446 : :
1447 : : // start point is the last node
1448 : 24 : if ( !foundStart && qgsDoubleNear( distanceTraversed, startDistance ) )
1449 : : {
1450 : 0 : substringPoints << QgsPoint( pointType, prevX, prevY, prevZ, prevM )
1451 : 0 : << QgsPoint( pointType, prevX, prevY, prevZ, prevM )
1452 : 0 : << QgsPoint( pointType, prevX, prevY, prevZ, prevM );
1453 : 0 : }
1454 : :
1455 : 24 : std::unique_ptr< QgsCircularString > result = std::make_unique< QgsCircularString >();
1456 : 24 : result->setPoints( substringPoints );
1457 : 24 : return result.release();
1458 : 26 : }
1459 : :
1460 : 28 : bool QgsCircularString::addZValue( double zValue )
1461 : : {
1462 : 28 : if ( QgsWkbTypes::hasZ( mWkbType ) )
1463 : 1 : return false;
1464 : :
1465 : 27 : clearCache();
1466 : 27 : mWkbType = QgsWkbTypes::addZ( mWkbType );
1467 : :
1468 : 27 : int nPoints = numPoints();
1469 : 27 : mZ.clear();
1470 : 27 : mZ.reserve( nPoints );
1471 : 92 : for ( int i = 0; i < nPoints; ++i )
1472 : : {
1473 : 65 : mZ << zValue;
1474 : 65 : }
1475 : 27 : return true;
1476 : 28 : }
1477 : :
1478 : 28 : bool QgsCircularString::addMValue( double mValue )
1479 : : {
1480 : 28 : if ( QgsWkbTypes::hasM( mWkbType ) )
1481 : 1 : return false;
1482 : :
1483 : 27 : clearCache();
1484 : 27 : mWkbType = QgsWkbTypes::addM( mWkbType );
1485 : :
1486 : 27 : int nPoints = numPoints();
1487 : 27 : mM.clear();
1488 : 27 : mM.reserve( nPoints );
1489 : 96 : for ( int i = 0; i < nPoints; ++i )
1490 : : {
1491 : 69 : mM << mValue;
1492 : 69 : }
1493 : 27 : return true;
1494 : 28 : }
1495 : :
1496 : 41 : bool QgsCircularString::dropZValue()
1497 : : {
1498 : 41 : if ( !QgsWkbTypes::hasZ( mWkbType ) )
1499 : 19 : return false;
1500 : :
1501 : 22 : clearCache();
1502 : :
1503 : 22 : mWkbType = QgsWkbTypes::dropZ( mWkbType );
1504 : 22 : mZ.clear();
1505 : 22 : return true;
1506 : 41 : }
1507 : :
1508 : 45 : bool QgsCircularString::dropMValue()
1509 : : {
1510 : 45 : if ( !QgsWkbTypes::hasM( mWkbType ) )
1511 : 25 : return false;
1512 : :
1513 : 20 : clearCache();
1514 : :
1515 : 20 : mWkbType = QgsWkbTypes::dropM( mWkbType );
1516 : 20 : mM.clear();
1517 : 20 : return true;
1518 : 45 : }
1519 : :
1520 : 4 : void QgsCircularString::swapXy()
1521 : : {
1522 : 4 : std::swap( mX, mY );
1523 : 4 : clearCache();
1524 : 4 : }
|
