QGIS API Documentation 3.34.0-Prizren (ffbdd678812)
Loading...
Searching...
No Matches
qgsgeometryutils.cpp
Go to the documentation of this file.
1/***************************************************************************
2 qgsgeometryutils.cpp
3 -------------------------------------------------------------------
4Date : 21 Nov 2014
5Copyright : (C) 2014 by Marco Hugentobler
6email : marco.hugentobler at sourcepole dot com
7 ***************************************************************************
8 * *
9 * This program is free software; you can redistribute it and/or modify *
10 * it under the terms of the GNU General Public License as published by *
11 * the Free Software Foundation; either version 2 of the License, or *
12 * (at your option) any later version. *
13 * *
14 ***************************************************************************/
15
16#include "qgsgeometryutils.h"
17
18#include "qgscurve.h"
19#include "qgscurvepolygon.h"
21#include "qgslinestring.h"
22#include "qgswkbptr.h"
23
24#include <memory>
25#include <QStringList>
26#include <QVector>
27#include <QRegularExpression>
28#include <nlohmann/json.hpp>
29
30QVector<QgsLineString *> QgsGeometryUtils::extractLineStrings( const QgsAbstractGeometry *geom )
31{
32 QVector< QgsLineString * > linestrings;
33 if ( !geom )
34 return linestrings;
35
36 QVector< const QgsAbstractGeometry * > geometries;
37 geometries << geom;
38 while ( ! geometries.isEmpty() )
39 {
40 const QgsAbstractGeometry *g = geometries.takeFirst();
41 if ( const QgsCurve *curve = qgsgeometry_cast< const QgsCurve * >( g ) )
42 {
43 linestrings << static_cast< QgsLineString * >( curve->segmentize() );
44 }
45 else if ( const QgsGeometryCollection *collection = qgsgeometry_cast< const QgsGeometryCollection * >( g ) )
46 {
47 for ( int i = 0; i < collection->numGeometries(); ++i )
48 {
49 geometries.append( collection->geometryN( i ) );
50 }
51 }
52 else if ( const QgsCurvePolygon *curvePolygon = qgsgeometry_cast< const QgsCurvePolygon * >( g ) )
53 {
54 if ( curvePolygon->exteriorRing() )
55 linestrings << static_cast< QgsLineString * >( curvePolygon->exteriorRing()->segmentize() );
56
57 for ( int i = 0; i < curvePolygon->numInteriorRings(); ++i )
58 {
59 linestrings << static_cast< QgsLineString * >( curvePolygon->interiorRing( i )->segmentize() );
60 }
61 }
62 }
63 return linestrings;
64}
65
67{
68 double minDist = std::numeric_limits<double>::max();
69 double currentDist = 0;
70 QgsPoint minDistPoint;
71 id = QgsVertexId(); // set as invalid
72
73 if ( geom.isEmpty() || pt.isEmpty() )
74 return minDistPoint;
75
76 QgsVertexId vertexId;
77 QgsPoint vertex;
78 while ( geom.nextVertex( vertexId, vertex ) )
79 {
80 currentDist = QgsGeometryUtils::sqrDistance2D( pt, vertex );
81 // The <= is on purpose: for geometries with closing vertices, this ensures
82 // that the closing vertex is returned. For the vertex tool, the rubberband
83 // of the closing vertex is above the opening vertex, hence with the <=
84 // situations where the covered opening vertex rubberband is selected are
85 // avoided.
86 if ( currentDist <= minDist )
87 {
88 minDist = currentDist;
89 minDistPoint = vertex;
90 id.part = vertexId.part;
91 id.ring = vertexId.ring;
92 id.vertex = vertexId.vertex;
93 id.type = vertexId.type;
94 }
95 }
96
97 return minDistPoint;
98}
99
101{
103 QgsVertexId vertexAfter;
104 geometry.closestSegment( point, closestPoint, vertexAfter, nullptr, DEFAULT_SEGMENT_EPSILON );
105 if ( vertexAfter.isValid() )
106 {
107 const QgsPoint pointAfter = geometry.vertexAt( vertexAfter );
108 if ( vertexAfter.vertex > 0 )
109 {
110 QgsVertexId vertexBefore = vertexAfter;
111 vertexBefore.vertex--;
112 const QgsPoint pointBefore = geometry.vertexAt( vertexBefore );
113 const double length = pointBefore.distance( pointAfter );
114 const double distance = pointBefore.distance( closestPoint );
115
116 if ( qgsDoubleNear( distance, 0.0 ) )
117 closestPoint = pointBefore;
118 else if ( qgsDoubleNear( distance, length ) )
119 closestPoint = pointAfter;
120 else
121 {
122 if ( QgsWkbTypes::hasZ( geometry.wkbType() ) && length )
123 closestPoint.addZValue( pointBefore.z() + ( pointAfter.z() - pointBefore.z() ) * distance / length );
124 if ( QgsWkbTypes::hasM( geometry.wkbType() ) )
125 closestPoint.addMValue( pointBefore.m() + ( pointAfter.m() - pointBefore.m() ) * distance / length );
126 }
127 }
128 }
129
130 return closestPoint;
131}
132
134{
135 double currentDist = 0;
136 QgsVertexId vertexId;
137 QgsPoint vertex;
138 while ( geom.nextVertex( vertexId, vertex ) )
139 {
140 if ( vertexId == id )
141 {
142 //found target vertex
143 return currentDist;
144 }
145 currentDist += geom.segmentLength( vertexId );
146 }
147
148 //could not find target vertex
149 return -1;
150}
151
152bool QgsGeometryUtils::verticesAtDistance( const QgsAbstractGeometry &geometry, double distance, QgsVertexId &previousVertex, QgsVertexId &nextVertex )
153{
154 double currentDist = 0;
155 previousVertex = QgsVertexId();
156 nextVertex = QgsVertexId();
157
158 if ( geometry.isEmpty() )
159 {
160 return false;
161 }
162
163 QgsPoint point;
164 QgsPoint previousPoint;
165
166 if ( qgsDoubleNear( distance, 0.0 ) )
167 {
168 geometry.nextVertex( previousVertex, point );
169 nextVertex = previousVertex;
170 return true;
171 }
172
173 bool first = true;
174 while ( currentDist < distance && geometry.nextVertex( nextVertex, point ) )
175 {
176 if ( !first && nextVertex.part == previousVertex.part && nextVertex.ring == previousVertex.ring )
177 {
178 currentDist += std::sqrt( QgsGeometryUtils::sqrDistance2D( previousPoint, point ) );
179 }
180
181 if ( qgsDoubleNear( currentDist, distance ) )
182 {
183 // exact hit!
184 previousVertex = nextVertex;
185 return true;
186 }
187
188 if ( currentDist > distance )
189 {
190 return true;
191 }
192
193 previousVertex = nextVertex;
194 previousPoint = point;
195 first = false;
196 }
197
198 //could not find target distance
199 return false;
200}
201
202double QgsGeometryUtils::sqrDistance2D( const QgsPoint &pt1, const QgsPoint &pt2 )
203{
204 return ( pt1.x() - pt2.x() ) * ( pt1.x() - pt2.x() ) + ( pt1.y() - pt2.y() ) * ( pt1.y() - pt2.y() );
205}
206
207double QgsGeometryUtils::sqrDistToLine( double ptX, double ptY, double x1, double y1, double x2, double y2, double &minDistX, double &minDistY, double epsilon )
208{
209 minDistX = x1;
210 minDistY = y1;
211
212 double dx = x2 - x1;
213 double dy = y2 - y1;
214
215 if ( !qgsDoubleNear( dx, 0.0 ) || !qgsDoubleNear( dy, 0.0 ) )
216 {
217 const double t = ( ( ptX - x1 ) * dx + ( ptY - y1 ) * dy ) / ( dx * dx + dy * dy );
218 if ( t > 1 )
219 {
220 minDistX = x2;
221 minDistY = y2;
222 }
223 else if ( t > 0 )
224 {
225 minDistX += dx * t;
226 minDistY += dy * t;
227 }
228 }
229
230 dx = ptX - minDistX;
231 dy = ptY - minDistY;
232
233 const double dist = dx * dx + dy * dy;
234
235 //prevent rounding errors if the point is directly on the segment
236 if ( qgsDoubleNear( dist, 0.0, epsilon ) )
237 {
238 minDistX = ptX;
239 minDistY = ptY;
240 return 0.0;
241 }
242
243 return dist;
244}
245
246double QgsGeometryUtils::distToInfiniteLine( const QgsPoint &point, const QgsPoint &linePoint1, const QgsPoint &linePoint2, double epsilon )
247{
248 const double area = std::abs(
249 ( linePoint1.x() - linePoint2.x() ) * ( point.y() - linePoint2.y() ) -
250 ( linePoint1.y() - linePoint2.y() ) * ( point.x() - linePoint2.x() )
251 );
252
253 const double length = std::sqrt(
254 std::pow( linePoint1.x() - linePoint2.x(), 2 ) +
255 std::pow( linePoint1.y() - linePoint2.y(), 2 )
256 );
257
258 const double distance = area / length;
259 return qgsDoubleNear( distance, 0.0, epsilon ) ? 0.0 : distance;
260}
261
262bool QgsGeometryUtils::lineIntersection( const QgsPoint &p1, QgsVector v1, const QgsPoint &p2, QgsVector v2, QgsPoint &intersection )
263{
264 const double d = v1.y() * v2.x() - v1.x() * v2.y();
265
266 if ( qgsDoubleNear( d, 0 ) )
267 return false;
268
269 const double dx = p2.x() - p1.x();
270 const double dy = p2.y() - p1.y();
271 const double k = ( dy * v2.x() - dx * v2.y() ) / d;
272
273 intersection = QgsPoint( p1.x() + v1.x() * k, p1.y() + v1.y() * k );
274
275 // z and m support for intersection point
277
278 return true;
279}
280
281bool QgsGeometryUtils::segmentIntersection( const QgsPoint &p1, const QgsPoint &p2, const QgsPoint &q1, const QgsPoint &q2, QgsPoint &intersectionPoint, bool &isIntersection, const double tolerance, bool acceptImproperIntersection )
282{
283 isIntersection = false;
284
285 QgsVector v( p2.x() - p1.x(), p2.y() - p1.y() );
286 QgsVector w( q2.x() - q1.x(), q2.y() - q1.y() );
287 const double vl = v.length();
288 const double wl = w.length();
289
290 if ( qgsDoubleNear( vl, 0.0, tolerance ) || qgsDoubleNear( wl, 0.0, tolerance ) )
291 {
292 return false;
293 }
294 v = v / vl;
295 w = w / wl;
296
297 if ( !QgsGeometryUtils::lineIntersection( p1, v, q1, w, intersectionPoint ) )
298 {
299 return false;
300 }
301
302 isIntersection = true;
303 if ( acceptImproperIntersection )
304 {
305 if ( ( p1 == q1 ) || ( p1 == q2 ) )
306 {
307 intersectionPoint = p1;
308 return true;
309 }
310 else if ( ( p2 == q1 ) || ( p2 == q2 ) )
311 {
312 intersectionPoint = p2;
313 return true;
314 }
315
316 double x, y;
317 if (
318 // intersectionPoint = p1
319 qgsDoubleNear( QgsGeometryUtils::sqrDistToLine( p1.x(), p1.y(), q1.x(), q1.y(), q2.x(), q2.y(), x, y, tolerance ), 0.0, tolerance ) ||
320 // intersectionPoint = p2
321 qgsDoubleNear( QgsGeometryUtils::sqrDistToLine( p2.x(), p2.y(), q1.x(), q1.y(), q2.x(), q2.y(), x, y, tolerance ), 0.0, tolerance ) ||
322 // intersectionPoint = q1
323 qgsDoubleNear( QgsGeometryUtils::sqrDistToLine( q1.x(), q1.y(), p1.x(), p1.y(), p2.x(), p2.y(), x, y, tolerance ), 0.0, tolerance ) ||
324 // intersectionPoint = q2
325 qgsDoubleNear( QgsGeometryUtils::sqrDistToLine( q2.x(), q2.y(), p1.x(), p1.y(), p2.x(), p2.y(), x, y, tolerance ), 0.0, tolerance )
326 )
327 {
328 return true;
329 }
330 }
331
332 const double lambdav = QgsVector( intersectionPoint.x() - p1.x(), intersectionPoint.y() - p1.y() ) * v;
333 if ( lambdav < 0. + tolerance || lambdav > vl - tolerance )
334 return false;
335
336 const double lambdaw = QgsVector( intersectionPoint.x() - q1.x(), intersectionPoint.y() - q1.y() ) * w;
337 return !( lambdaw < 0. + tolerance || lambdaw >= wl - tolerance );
338
339}
340
341bool QgsGeometryUtils::lineCircleIntersection( const QgsPointXY &center, const double radius,
342 const QgsPointXY &linePoint1, const QgsPointXY &linePoint2,
343 QgsPointXY &intersection )
344{
345 // formula taken from http://mathworld.wolfram.com/Circle-LineIntersection.html
346
347 const double x1 = linePoint1.x() - center.x();
348 const double y1 = linePoint1.y() - center.y();
349 const double x2 = linePoint2.x() - center.x();
350 const double y2 = linePoint2.y() - center.y();
351 const double dx = x2 - x1;
352 const double dy = y2 - y1;
353
354 const double dr2 = std::pow( dx, 2 ) + std::pow( dy, 2 );
355 const double d = x1 * y2 - x2 * y1;
356
357 const double disc = std::pow( radius, 2 ) * dr2 - std::pow( d, 2 );
358
359 if ( disc < 0 )
360 {
361 //no intersection or tangent
362 return false;
363 }
364 else
365 {
366 // two solutions
367 const int sgnDy = dy < 0 ? -1 : 1;
368
369 const double sqrDisc = std::sqrt( disc );
370
371 const double ax = center.x() + ( d * dy + sgnDy * dx * sqrDisc ) / dr2;
372 const double ay = center.y() + ( -d * dx + std::fabs( dy ) * sqrDisc ) / dr2;
373 const QgsPointXY p1( ax, ay );
374
375 const double bx = center.x() + ( d * dy - sgnDy * dx * sqrDisc ) / dr2;
376 const double by = center.y() + ( -d * dx - std::fabs( dy ) * sqrDisc ) / dr2;
377 const QgsPointXY p2( bx, by );
378
379 // snap to nearest intersection
380
381 if ( intersection.sqrDist( p1 ) < intersection.sqrDist( p2 ) )
382 {
383 intersection.set( p1.x(), p1.y() );
384 }
385 else
386 {
387 intersection.set( p2.x(), p2.y() );
388 }
389 return true;
390 }
391}
392
393// based on public domain work by 3/26/2005 Tim Voght
394// see http://paulbourke.net/geometry/circlesphere/tvoght.c
395int QgsGeometryUtils::circleCircleIntersections( const QgsPointXY &center1, const double r1, const QgsPointXY &center2, const double r2, QgsPointXY &intersection1, QgsPointXY &intersection2 )
396{
397 // determine the straight-line distance between the centers
398 const double d = center1.distance( center2 );
399
400 // check if the circles intersect at only 1 point, either "externally" or "internally"
401 const bool singleSolutionExt = qgsDoubleNear( d, r1 + r2 );
402 const bool singleSolutionInt = qgsDoubleNear( d, std::fabs( r1 - r2 ) );
403
404 // check for solvability
405 if ( !singleSolutionExt && d > ( r1 + r2 ) )
406 {
407 // no solution. circles do not intersect.
408 return 0;
409 }
410 else if ( !singleSolutionInt && d < std::fabs( r1 - r2 ) )
411 {
412 // no solution. one circle is contained in the other
413 return 0;
414 }
415 else if ( qgsDoubleNear( d, 0 ) && ( qgsDoubleNear( r1, r2 ) ) )
416 {
417 // no solutions, the circles coincide
418 return 0;
419 }
420
421 /* 'point 2' is the point where the line through the circle
422 * intersection points crosses the line between the circle
423 * centers.
424 */
425
426 /* Determine the distance 'a' from point 0 to point 2.
427 * In the general case, a = ( ( r1 * r1 ) - ( r2 * r2 ) + ( d * d ) ) / ( 2.0 * d ).
428 * If d = r1 + r2 or d = r1 - r2 (i.e. r1 > r2), then a = r1; if d = r2 - r1 (i.e. r2 > r1), then a = -r1.
429 */
430 const double a = singleSolutionExt ? r1 : ( singleSolutionInt ? ( r1 > r2 ? r1 : -r1 ) : ( ( r1 * r1 ) - ( r2 * r2 ) + ( d * d ) ) / ( 2.0 * d ) );
431
432 /* dx and dy are the vertical and horizontal distances between
433 * the circle centers.
434 */
435 const double dx = center2.x() - center1.x();
436 const double dy = center2.y() - center1.y();
437
438 // Determine the coordinates of point 2.
439 const double x2 = center1.x() + ( dx * a / d );
440 const double y2 = center1.y() + ( dy * a / d );
441
442 // only 1 solution
443 if ( singleSolutionExt || singleSolutionInt )
444 {
445 intersection1 = QgsPointXY( x2, y2 );
446 intersection2 = QgsPointXY( x2, y2 );
447
448 return 1;
449 }
450
451 // 2 solutions
452
453 /* Determine the distance from point 2 to either of the
454 * intersection points.
455 */
456 const double h = std::sqrt( ( r1 * r1 ) - ( a * a ) );
457
458 /* Now determine the offsets of the intersection points from
459 * point 2.
460 */
461 const double rx = dy * ( h / d );
462 const double ry = dx * ( h / d );
463
464 // determine the absolute intersection points
465 intersection1 = QgsPointXY( x2 + rx, y2 - ry );
466 intersection2 = QgsPointXY( x2 - rx, y2 + ry );
467
468 return 2;
469}
470
471// Using https://stackoverflow.com/a/1351794/1861260
472// and inspired by http://csharphelper.com/blog/2014/11/find-the-tangent-lines-between-a-point-and-a-circle-in-c/
473bool QgsGeometryUtils::tangentPointAndCircle( const QgsPointXY &center, double radius, const QgsPointXY &p, QgsPointXY &pt1, QgsPointXY &pt2 )
474{
475 // distance from point to center of circle
476 const double dx = center.x() - p.x();
477 const double dy = center.y() - p.y();
478 const double distanceSquared = dx * dx + dy * dy;
479 const double radiusSquared = radius * radius;
480 if ( distanceSquared < radiusSquared )
481 {
482 // point is inside circle!
483 return false;
484 }
485
486 // distance from point to tangent point, using pythagoras
487 const double distanceToTangent = std::sqrt( distanceSquared - radiusSquared );
488
489 // tangent points are those where the original circle intersects a circle centered
490 // on p with radius distanceToTangent
491 circleCircleIntersections( center, radius, p, distanceToTangent, pt1, pt2 );
492
493 return true;
494}
495
496// inspired by http://csharphelper.com/blog/2014/12/find-the-tangent-lines-between-two-circles-in-c/
497int QgsGeometryUtils::circleCircleOuterTangents( const QgsPointXY &center1, double radius1, const QgsPointXY &center2, double radius2, QgsPointXY &line1P1, QgsPointXY &line1P2, QgsPointXY &line2P1, QgsPointXY &line2P2 )
498{
499 if ( radius1 > radius2 )
500 return circleCircleOuterTangents( center2, radius2, center1, radius1, line1P1, line1P2, line2P1, line2P2 );
501
502 const double radius2a = radius2 - radius1;
503 if ( !tangentPointAndCircle( center2, radius2a, center1, line1P2, line2P2 ) )
504 {
505 // there are no tangents
506 return 0;
507 }
508
509 // get the vector perpendicular to the
510 // first tangent with length radius1
511 QgsVector v1( -( line1P2.y() - center1.y() ), line1P2.x() - center1.x() );
512 const double v1Length = v1.length();
513 v1 = v1 * ( radius1 / v1Length );
514
515 // offset the tangent vector's points
516 line1P1 = center1 + v1;
517 line1P2 = line1P2 + v1;
518
519 // get the vector perpendicular to the
520 // second tangent with length radius1
521 QgsVector v2( line2P2.y() - center1.y(), -( line2P2.x() - center1.x() ) );
522 const double v2Length = v2.length();
523 v2 = v2 * ( radius1 / v2Length );
524
525 // offset the tangent vector's points
526 line2P1 = center1 + v2;
527 line2P2 = line2P2 + v2;
528
529 return 2;
530}
531
532// inspired by http://csharphelper.com/blog/2014/12/find-the-tangent-lines-between-two-circles-in-c/
533int QgsGeometryUtils::circleCircleInnerTangents( const QgsPointXY &center1, double radius1, const QgsPointXY &center2, double radius2, QgsPointXY &line1P1, QgsPointXY &line1P2, QgsPointXY &line2P1, QgsPointXY &line2P2 )
534{
535 if ( radius1 > radius2 )
536 return circleCircleInnerTangents( center2, radius2, center1, radius1, line1P1, line1P2, line2P1, line2P2 );
537
538 // determine the straight-line distance between the centers
539 const double d = center1.distance( center2 );
540 const double radius1a = radius1 + radius2;
541
542 // check for solvability
543 if ( d <= radius1a || qgsDoubleNear( d, radius1a ) )
544 {
545 // no solution. circles intersect or touch.
546 return 0;
547 }
548
549 if ( !tangentPointAndCircle( center1, radius1a, center2, line1P2, line2P2 ) )
550 {
551 // there are no tangents
552 return 0;
553 }
554
555 // get the vector perpendicular to the
556 // first tangent with length radius2
557 QgsVector v1( ( line1P2.y() - center2.y() ), -( line1P2.x() - center2.x() ) );
558 const double v1Length = v1.length();
559 v1 = v1 * ( radius2 / v1Length );
560
561 // offset the tangent vector's points
562 line1P1 = center2 + v1;
563 line1P2 = line1P2 + v1;
564
565 // get the vector perpendicular to the
566 // second tangent with length radius2
567 QgsVector v2( -( line2P2.y() - center2.y() ), line2P2.x() - center2.x() );
568 const double v2Length = v2.length();
569 v2 = v2 * ( radius2 / v2Length );
570
571 // offset the tangent vector's points in opposite direction
572 line2P1 = center2 + v2;
573 line2P2 = line2P2 + v2;
574
575 return 2;
576}
577
578QVector<QgsGeometryUtils::SelfIntersection> QgsGeometryUtils::selfIntersections( const QgsAbstractGeometry *geom, int part, int ring, double tolerance )
579{
580 QVector<SelfIntersection> intersections;
581
582 const int n = geom->vertexCount( part, ring );
583 const bool isClosed = geom->vertexAt( QgsVertexId( part, ring, 0 ) ) == geom->vertexAt( QgsVertexId( part, ring, n - 1 ) );
584
585 // Check every pair of segments for intersections
586 for ( int i = 0, j = 1; j < n; i = j++ )
587 {
588 const QgsPoint pi = geom->vertexAt( QgsVertexId( part, ring, i ) );
589 const QgsPoint pj = geom->vertexAt( QgsVertexId( part, ring, j ) );
590 if ( QgsGeometryUtils::sqrDistance2D( pi, pj ) < tolerance * tolerance ) continue;
591
592 // Don't test neighboring edges
593 const int start = j + 1;
594 const int end = i == 0 && isClosed ? n - 1 : n;
595 for ( int k = start, l = start + 1; l < end; k = l++ )
596 {
597 const QgsPoint pk = geom->vertexAt( QgsVertexId( part, ring, k ) );
598 const QgsPoint pl = geom->vertexAt( QgsVertexId( part, ring, l ) );
599
600 QgsPoint inter;
601 bool intersection = false;
602 if ( !QgsGeometryUtils::segmentIntersection( pi, pj, pk, pl, inter, intersection, tolerance ) ) continue;
603
605 s.segment1 = i;
606 s.segment2 = k;
607 if ( s.segment1 > s.segment2 )
608 {
609 std::swap( s.segment1, s.segment2 );
610 }
611 s.point = inter;
612 intersections.append( s );
613 }
614 }
615 return intersections;
616}
617
618int QgsGeometryUtils::leftOfLine( const QgsPoint &point, const QgsPoint &p1, const QgsPoint &p2 )
619{
620 return leftOfLine( point.x(), point.y(), p1.x(), p1.y(), p2.x(), p2.y() );
621}
622
623int QgsGeometryUtils::leftOfLine( const double x, const double y, const double x1, const double y1, const double x2, const double y2 )
624{
625 const double f1 = x - x1;
626 const double f2 = y2 - y1;
627 const double f3 = y - y1;
628 const double f4 = x2 - x1;
629 const double test = ( f1 * f2 - f3 * f4 );
630 // return -1, 0, or 1
631 return qgsDoubleNear( test, 0.0 ) ? 0 : ( test < 0 ? -1 : 1 );
632}
633
634QgsPoint QgsGeometryUtils::pointOnLineWithDistance( const QgsPoint &startPoint, const QgsPoint &directionPoint, double distance )
635{
636 double x, y;
637 pointOnLineWithDistance( startPoint.x(), startPoint.y(), directionPoint.x(), directionPoint.y(), distance, x, y );
638 return QgsPoint( x, y );
639}
640
641void QgsGeometryUtils::pointOnLineWithDistance( double x1, double y1, double x2, double y2, double distance, double &x, double &y, double *z1, double *z2, double *z, double *m1, double *m2, double *m )
642{
643 const double dx = x2 - x1;
644 const double dy = y2 - y1;
645 const double length = std::sqrt( dx * dx + dy * dy );
646
647 if ( qgsDoubleNear( length, 0.0 ) )
648 {
649 x = x1;
650 y = y1;
651 if ( z && z1 )
652 *z = *z1;
653 if ( m && m1 )
654 *m = *m1;
655 }
656 else
657 {
658 const double scaleFactor = distance / length;
659 x = x1 + dx * scaleFactor;
660 y = y1 + dy * scaleFactor;
661 if ( z && z1 && z2 )
662 *z = *z1 + ( *z2 - *z1 ) * scaleFactor;
663 if ( m && m1 && m2 )
664 *m = *m1 + ( *m2 - *m1 ) * scaleFactor;
665 }
666}
667
668void QgsGeometryUtils::perpendicularOffsetPointAlongSegment( double x1, double y1, double x2, double y2, double proportion, double offset, double *x, double *y )
669{
670 // calculate point along segment
671 const double mX = x1 + ( x2 - x1 ) * proportion;
672 const double mY = y1 + ( y2 - y1 ) * proportion;
673 const double pX = x1 - x2;
674 const double pY = y1 - y2;
675 double normalX = -pY;
676 double normalY = pX; //#spellok
677 const double normalLength = sqrt( ( normalX * normalX ) + ( normalY * normalY ) ); //#spellok
678 normalX /= normalLength;
679 normalY /= normalLength; //#spellok
680
681 *x = mX + offset * normalX;
682 *y = mY + offset * normalY; //#spellok
683}
684
685QgsPoint QgsGeometryUtils::interpolatePointOnArc( const QgsPoint &pt1, const QgsPoint &pt2, const QgsPoint &pt3, double distance )
686{
687 double centerX, centerY, radius;
688 circleCenterRadius( pt1, pt2, pt3, radius, centerX, centerY );
689
690 const double theta = distance / radius; // angle subtended
691 const double anglePt1 = std::atan2( pt1.y() - centerY, pt1.x() - centerX );
692 const double anglePt2 = std::atan2( pt2.y() - centerY, pt2.x() - centerX );
693 const double anglePt3 = std::atan2( pt3.y() - centerY, pt3.x() - centerX );
694 const bool isClockwise = circleClockwise( anglePt1, anglePt2, anglePt3 );
695 const double angleDest = anglePt1 + ( isClockwise ? -theta : theta );
696
697 const double x = centerX + radius * ( std::cos( angleDest ) );
698 const double y = centerY + radius * ( std::sin( angleDest ) );
699
700 const double z = pt1.is3D() ?
701 interpolateArcValue( angleDest, anglePt1, anglePt2, anglePt3, pt1.z(), pt2.z(), pt3.z() )
702 : 0;
703 const double m = pt1.isMeasure() ?
704 interpolateArcValue( angleDest, anglePt1, anglePt2, anglePt3, pt1.m(), pt2.m(), pt3.m() )
705 : 0;
706
707 return QgsPoint( pt1.wkbType(), x, y, z, m );
708}
709
710double QgsGeometryUtils::ccwAngle( double dy, double dx )
711{
712 const double angle = std::atan2( dy, dx ) * 180 / M_PI;
713 if ( angle < 0 )
714 {
715 return 360 + angle;
716 }
717 else if ( angle > 360 )
718 {
719 return 360 - angle;
720 }
721 return angle;
722}
723
724void QgsGeometryUtils::circleCenterRadius( const QgsPoint &pt1, const QgsPoint &pt2, const QgsPoint &pt3, double &radius, double &centerX, double &centerY )
725{
726 double dx21, dy21, dx31, dy31, h21, h31, d;
727
728 //closed circle
729 if ( qgsDoubleNear( pt1.x(), pt3.x() ) && qgsDoubleNear( pt1.y(), pt3.y() ) )
730 {
731 centerX = ( pt1.x() + pt2.x() ) / 2.0;
732 centerY = ( pt1.y() + pt2.y() ) / 2.0;
733 radius = std::sqrt( std::pow( centerX - pt1.x(), 2.0 ) + std::pow( centerY - pt1.y(), 2.0 ) );
734 return;
735 }
736
737 // Using Cartesian circumcenter eguations from page https://en.wikipedia.org/wiki/Circumscribed_circle
738 dx21 = pt2.x() - pt1.x();
739 dy21 = pt2.y() - pt1.y();
740 dx31 = pt3.x() - pt1.x();
741 dy31 = pt3.y() - pt1.y();
742
743 h21 = std::pow( dx21, 2.0 ) + std::pow( dy21, 2.0 );
744 h31 = std::pow( dx31, 2.0 ) + std::pow( dy31, 2.0 );
745
746 // 2*Cross product, d<0 means clockwise and d>0 counterclockwise sweeping angle
747 d = 2 * ( dx21 * dy31 - dx31 * dy21 );
748
749 // Check colinearity, Cross product = 0
750 if ( qgsDoubleNear( std::fabs( d ), 0.0, 0.00000000001 ) )
751 {
752 radius = -1.0;
753 return;
754 }
755
756 // Calculate centroid coordinates and radius
757 centerX = pt1.x() + ( h21 * dy31 - h31 * dy21 ) / d;
758 centerY = pt1.y() - ( h21 * dx31 - h31 * dx21 ) / d;
759 radius = std::sqrt( std::pow( centerX - pt1.x(), 2.0 ) + std::pow( centerY - pt1.y(), 2.0 ) );
760}
761
762bool QgsGeometryUtils::circleClockwise( double angle1, double angle2, double angle3 )
763{
764 if ( angle3 >= angle1 )
765 {
766 return !( angle2 > angle1 && angle2 < angle3 );
767 }
768 else
769 {
770 return !( angle2 > angle1 || angle2 < angle3 );
771 }
772}
773
774bool QgsGeometryUtils::circleAngleBetween( double angle, double angle1, double angle2, bool clockwise )
775{
776 if ( clockwise )
777 {
778 if ( angle2 < angle1 )
779 {
780 return ( angle <= angle1 && angle >= angle2 );
781 }
782 else
783 {
784 return ( angle <= angle1 || angle >= angle2 );
785 }
786 }
787 else
788 {
789 if ( angle2 > angle1 )
790 {
791 return ( angle >= angle1 && angle <= angle2 );
792 }
793 else
794 {
795 return ( angle >= angle1 || angle <= angle2 );
796 }
797 }
798}
799
800bool QgsGeometryUtils::angleOnCircle( double angle, double angle1, double angle2, double angle3 )
801{
802 const bool clockwise = circleClockwise( angle1, angle2, angle3 );
803 return circleAngleBetween( angle, angle1, angle3, clockwise );
804}
805
806double QgsGeometryUtils::circleLength( double x1, double y1, double x2, double y2, double x3, double y3 )
807{
808 double centerX, centerY, radius;
809 circleCenterRadius( QgsPoint( x1, y1 ), QgsPoint( x2, y2 ), QgsPoint( x3, y3 ), radius, centerX, centerY );
810 double length = M_PI / 180.0 * radius * sweepAngle( centerX, centerY, x1, y1, x2, y2, x3, y3 );
811 if ( length < 0 )
812 {
813 length = -length;
814 }
815 return length;
816}
817
818double QgsGeometryUtils::sweepAngle( double centerX, double centerY, double x1, double y1, double x2, double y2, double x3, double y3 )
819{
820 const double p1Angle = QgsGeometryUtils::ccwAngle( y1 - centerY, x1 - centerX );
821 const double p2Angle = QgsGeometryUtils::ccwAngle( y2 - centerY, x2 - centerX );
822 const double p3Angle = QgsGeometryUtils::ccwAngle( y3 - centerY, x3 - centerX );
823
824 if ( p3Angle >= p1Angle )
825 {
826 if ( p2Angle > p1Angle && p2Angle < p3Angle )
827 {
828 return ( p3Angle - p1Angle );
829 }
830 else
831 {
832 return ( - ( p1Angle + ( 360 - p3Angle ) ) );
833 }
834 }
835 else
836 {
837 if ( p2Angle < p1Angle && p2Angle > p3Angle )
838 {
839 return ( -( p1Angle - p3Angle ) );
840 }
841 else
842 {
843 return ( p3Angle + ( 360 - p1Angle ) );
844 }
845 }
846}
847
848bool QgsGeometryUtils::segmentMidPoint( const QgsPoint &p1, const QgsPoint &p2, QgsPoint &result, double radius, const QgsPoint &mousePos )
849{
850 const QgsPoint midPoint( ( p1.x() + p2.x() ) / 2.0, ( p1.y() + p2.y() ) / 2.0 );
851 const double midDist = std::sqrt( sqrDistance2D( p1, midPoint ) );
852 if ( radius < midDist )
853 {
854 return false;
855 }
856 const double centerMidDist = std::sqrt( radius * radius - midDist * midDist );
857 const double dist = radius - centerMidDist;
858
859 const double midDx = midPoint.x() - p1.x();
860 const double midDy = midPoint.y() - p1.y();
861
862 //get the four possible midpoints
863 QVector<QgsPoint> possibleMidPoints;
864 possibleMidPoints.append( pointOnLineWithDistance( midPoint, QgsPoint( midPoint.x() - midDy, midPoint.y() + midDx ), dist ) );
865 possibleMidPoints.append( pointOnLineWithDistance( midPoint, QgsPoint( midPoint.x() - midDy, midPoint.y() + midDx ), 2 * radius - dist ) );
866 possibleMidPoints.append( pointOnLineWithDistance( midPoint, QgsPoint( midPoint.x() + midDy, midPoint.y() - midDx ), dist ) );
867 possibleMidPoints.append( pointOnLineWithDistance( midPoint, QgsPoint( midPoint.x() + midDy, midPoint.y() - midDx ), 2 * radius - dist ) );
868
869 //take the closest one
870 double minDist = std::numeric_limits<double>::max();
871 int minDistIndex = -1;
872 for ( int i = 0; i < possibleMidPoints.size(); ++i )
873 {
874 const double currentDist = sqrDistance2D( mousePos, possibleMidPoints.at( i ) );
875 if ( currentDist < minDist )
876 {
877 minDistIndex = i;
878 minDist = currentDist;
879 }
880 }
881
882 if ( minDistIndex == -1 )
883 {
884 return false;
885 }
886
887 result = possibleMidPoints.at( minDistIndex );
888
889 // add z and m support if necessary
891
892 return true;
893}
894
895QgsPoint QgsGeometryUtils::segmentMidPointFromCenter( const QgsPoint &p1, const QgsPoint &p2, const QgsPoint &center, const bool useShortestArc )
896{
897 double midPointAngle = averageAngle( lineAngle( center.x(), center.y(), p1.x(), p1.y() ),
898 lineAngle( center.x(), center.y(), p2.x(), p2.y() ) );
899 if ( !useShortestArc )
900 midPointAngle += M_PI;
901 return center.project( center.distance( p1 ), midPointAngle * 180 / M_PI );
902}
903
904double QgsGeometryUtils::circleTangentDirection( const QgsPoint &tangentPoint, const QgsPoint &cp1,
905 const QgsPoint &cp2, const QgsPoint &cp3 )
906{
907 //calculate circle midpoint
908 double mX, mY, radius;
909 circleCenterRadius( cp1, cp2, cp3, radius, mX, mY );
910
911 const double p1Angle = QgsGeometryUtils::ccwAngle( cp1.y() - mY, cp1.x() - mX );
912 const double p2Angle = QgsGeometryUtils::ccwAngle( cp2.y() - mY, cp2.x() - mX );
913 const double p3Angle = QgsGeometryUtils::ccwAngle( cp3.y() - mY, cp3.x() - mX );
914 double angle = 0;
915 if ( circleClockwise( p1Angle, p2Angle, p3Angle ) )
916 {
917 angle = lineAngle( tangentPoint.x(), tangentPoint.y(), mX, mY ) - M_PI_2;
918 }
919 else
920 {
921 angle = lineAngle( mX, mY, tangentPoint.x(), tangentPoint.y() ) - M_PI_2;
922 }
923 if ( angle < 0 )
924 angle += 2 * M_PI;
925 return angle;
926}
927
928// Ported from PostGIS' pt_continues_arc
929bool QgsGeometryUtils::pointContinuesArc( const QgsPoint &a1, const QgsPoint &a2, const QgsPoint &a3, const QgsPoint &b, double distanceTolerance, double pointSpacingAngleTolerance )
930{
931 double centerX = 0;
932 double centerY = 0;
933 double radius = 0;
934 circleCenterRadius( a1, a2, a3, radius, centerX, centerY );
935
936 // Co-linear a1/a2/a3
937 if ( radius < 0.0 )
938 return false;
939
940 // distance of candidate point to center of arc a1-a2-a3
941 const double bDistance = std::sqrt( ( b.x() - centerX ) * ( b.x() - centerX ) +
942 ( b.y() - centerY ) * ( b.y() - centerY ) );
943
944 double diff = std::fabs( radius - bDistance );
945
946 auto arcAngle = []( const QgsPoint & a, const QgsPoint & b, const QgsPoint & c )->double
947 {
948 const double abX = b.x() - a.x();
949 const double abY = b.y() - a.y();
950
951 const double cbX = b.x() - c.x();
952 const double cbY = b.y() - c.y();
953
954 const double dot = ( abX * cbX + abY * cbY ); /* dot product */
955 const double cross = ( abX * cbY - abY * cbX ); /* cross product */
956
957 const double alpha = std::atan2( cross, dot );
958
959 return alpha;
960 };
961
962 // Is the point b on the circle?
963 if ( diff < distanceTolerance )
964 {
965 const double angle1 = arcAngle( a1, a2, a3 );
966 const double angle2 = arcAngle( a2, a3, b );
967
968 // Is the sweep angle similar to the previous one?
969 // We only consider a segment replaceable by an arc if the points within
970 // it are regularly spaced
971 diff = std::fabs( angle1 - angle2 );
972 if ( diff > pointSpacingAngleTolerance )
973 {
974 return false;
975 }
976
977 const int a2Side = leftOfLine( a2.x(), a2.y(), a1.x(), a1.y(), a3.x(), a3.y() );
978 const int bSide = leftOfLine( b.x(), b.y(), a1.x(), a1.y(), a3.x(), a3.y() );
979
980 // Is the point b on the same side of a1/a3 as the mid-point a2 is?
981 // If not, it's in the unbounded part of the circle, so it continues the arc, return true.
982 if ( bSide != a2Side )
983 return true;
984 }
985 return false;
986}
987
988void QgsGeometryUtils::segmentizeArc( const QgsPoint &p1, const QgsPoint &p2, const QgsPoint &p3, QgsPointSequence &points, double tolerance, QgsAbstractGeometry::SegmentationToleranceType toleranceType, bool hasZ, bool hasM )
989{
990 bool reversed = false;
991 const int segSide = segmentSide( p1, p3, p2 );
992
993 QgsPoint circlePoint1;
994 const QgsPoint &circlePoint2 = p2;
995 QgsPoint circlePoint3;
996
997 if ( segSide == -1 )
998 {
999 // Reverse !
1000 circlePoint1 = p3;
1001 circlePoint3 = p1;
1002 reversed = true;
1003 }
1004 else
1005 {
1006 circlePoint1 = p1;
1007 circlePoint3 = p3;
1008 }
1009
1010 //adapted code from PostGIS
1011 double radius = 0;
1012 double centerX = 0;
1013 double centerY = 0;
1014 circleCenterRadius( circlePoint1, circlePoint2, circlePoint3, radius, centerX, centerY );
1015
1016 if ( circlePoint1 != circlePoint3 && ( radius < 0 || qgsDoubleNear( segSide, 0.0 ) ) ) //points are colinear
1017 {
1018 points.append( p1 );
1019 points.append( p2 );
1020 points.append( p3 );
1021 return;
1022 }
1023
1024 double increment = tolerance; //one segment per degree
1025 if ( toleranceType == QgsAbstractGeometry::MaximumDifference )
1026 {
1027 // Ensure tolerance is not higher than twice the radius
1028 tolerance = std::min( tolerance, radius * 2 );
1029 const double halfAngle = std::acos( -tolerance / radius + 1 );
1030 increment = 2 * halfAngle;
1031 }
1032
1033 //angles of pt1, pt2, pt3
1034 const double a1 = std::atan2( circlePoint1.y() - centerY, circlePoint1.x() - centerX );
1035 double a2 = std::atan2( circlePoint2.y() - centerY, circlePoint2.x() - centerX );
1036 double a3 = std::atan2( circlePoint3.y() - centerY, circlePoint3.x() - centerX );
1037
1038 // Make segmentation symmetric
1039 const bool symmetric = true;
1040 if ( symmetric )
1041 {
1042 double angle = a3 - a1;
1043 // angle == 0 when full circle
1044 if ( angle <= 0 ) angle += M_PI * 2;
1045
1046 /* Number of segments in output */
1047 const int segs = ceil( angle / increment );
1048 /* Tweak increment to be regular for all the arc */
1049 increment = angle / segs;
1050 }
1051
1052 /* Adjust a3 up so we can increment from a1 to a3 cleanly */
1053 // a3 == a1 when full circle
1054 if ( a3 <= a1 )
1055 a3 += 2.0 * M_PI;
1056 if ( a2 < a1 )
1057 a2 += 2.0 * M_PI;
1058
1059 double x, y;
1060 double z = 0;
1061 double m = 0;
1062
1063 QVector<QgsPoint> stringPoints;
1064 stringPoints.insert( 0, circlePoint1 );
1065 if ( circlePoint2 != circlePoint3 && circlePoint1 != circlePoint2 ) //draw straight line segment if two points have the same position
1066 {
1067 Qgis::WkbType pointWkbType = Qgis::WkbType::Point;
1068 if ( hasZ )
1069 pointWkbType = QgsWkbTypes::addZ( pointWkbType );
1070 if ( hasM )
1071 pointWkbType = QgsWkbTypes::addM( pointWkbType );
1072
1073 // As we're adding the last point in any case, we'll avoid
1074 // including a point which is at less than 1% increment distance
1075 // from it (may happen to find them due to numbers approximation).
1076 // NOTE that this effectively allows in output some segments which
1077 // are more distant than requested. This is at most 1% off
1078 // from requested MaxAngle and less for MaxError.
1079 const double tolError = increment / 100;
1080 const double stopAngle = a3 - tolError;
1081 for ( double angle = a1 + increment; angle < stopAngle; angle += increment )
1082 {
1083 x = centerX + radius * std::cos( angle );
1084 y = centerY + radius * std::sin( angle );
1085
1086 if ( hasZ )
1087 {
1088 z = interpolateArcValue( angle, a1, a2, a3, circlePoint1.z(), circlePoint2.z(), circlePoint3.z() );
1089 }
1090 if ( hasM )
1091 {
1092 m = interpolateArcValue( angle, a1, a2, a3, circlePoint1.m(), circlePoint2.m(), circlePoint3.m() );
1093 }
1094
1095 stringPoints.insert( stringPoints.size(), QgsPoint( pointWkbType, x, y, z, m ) );
1096 }
1097 }
1098 stringPoints.insert( stringPoints.size(), circlePoint3 );
1099
1100 // TODO: check if or implement QgsPointSequence directly taking an iterator to append
1101 if ( reversed )
1102 {
1103 std::reverse( stringPoints.begin(), stringPoints.end() );
1104 }
1105 if ( ! points.empty() && stringPoints.front() == points.back() ) stringPoints.pop_front();
1106 points.append( stringPoints );
1107}
1108
1109int QgsGeometryUtils::segmentSide( const QgsPoint &pt1, const QgsPoint &pt3, const QgsPoint &pt2 )
1110{
1111 const double side = ( ( pt2.x() - pt1.x() ) * ( pt3.y() - pt1.y() ) - ( pt3.x() - pt1.x() ) * ( pt2.y() - pt1.y() ) );
1112 if ( side == 0.0 )
1113 {
1114 return 0;
1115 }
1116 else
1117 {
1118 if ( side < 0 )
1119 {
1120 return -1;
1121 }
1122 if ( side > 0 )
1123 {
1124 return 1;
1125 }
1126 return 0;
1127 }
1128}
1129
1130double QgsGeometryUtils::interpolateArcValue( double angle, double a1, double a2, double a3, double zm1, double zm2, double zm3 )
1131{
1132 /* Counter-clockwise sweep */
1133 if ( a1 < a2 )
1134 {
1135 if ( angle <= a2 )
1136 return zm1 + ( zm2 - zm1 ) * ( angle - a1 ) / ( a2 - a1 );
1137 else
1138 return zm2 + ( zm3 - zm2 ) * ( angle - a2 ) / ( a3 - a2 );
1139 }
1140 /* Clockwise sweep */
1141 else
1142 {
1143 if ( angle >= a2 )
1144 return zm1 + ( zm2 - zm1 ) * ( a1 - angle ) / ( a1 - a2 );
1145 else
1146 return zm2 + ( zm3 - zm2 ) * ( a2 - angle ) / ( a2 - a3 );
1147 }
1148}
1149
1150QgsPointSequence QgsGeometryUtils::pointsFromWKT( const QString &wktCoordinateList, bool is3D, bool isMeasure )
1151{
1152 const int dim = 2 + is3D + isMeasure;
1153 QgsPointSequence points;
1154
1155#if QT_VERSION < QT_VERSION_CHECK(5, 15, 0)
1156 const QStringList coordList = wktCoordinateList.split( ',', QString::SkipEmptyParts );
1157#else
1158 const QStringList coordList = wktCoordinateList.split( ',', Qt::SkipEmptyParts );
1159#endif
1160
1161 //first scan through for extra unexpected dimensions
1162 bool foundZ = false;
1163 bool foundM = false;
1164 const thread_local QRegularExpression rx( QStringLiteral( "\\s" ) );
1165 const thread_local QRegularExpression rxIsNumber( QStringLiteral( "^[+-]?(\\d\\.?\\d*[Ee][+\\-]?\\d+|(\\d+\\.\\d*|\\d*\\.\\d+)|\\d+)$" ) );
1166 for ( const QString &pointCoordinates : coordList )
1167 {
1168#if QT_VERSION < QT_VERSION_CHECK(5, 15, 0)
1169 QStringList coordinates = pointCoordinates.split( rx, QString::SkipEmptyParts );
1170#else
1171 const QStringList coordinates = pointCoordinates.split( rx, Qt::SkipEmptyParts );
1172#endif
1173
1174 // exit with an empty set if one list contains invalid value.
1175 if ( coordinates.filter( rxIsNumber ).size() != coordinates.size() )
1176 return points;
1177
1178 if ( coordinates.size() == 3 && !foundZ && !foundM && !is3D && !isMeasure )
1179 {
1180 // 3 dimensional coordinates, but not specifically marked as such. We allow this
1181 // anyway and upgrade geometry to have Z dimension
1182 foundZ = true;
1183 }
1184 else if ( coordinates.size() >= 4 && ( !( is3D || foundZ ) || !( isMeasure || foundM ) ) )
1185 {
1186 // 4 (or more) dimensional coordinates, but not specifically marked as such. We allow this
1187 // anyway and upgrade geometry to have Z&M dimensions
1188 foundZ = true;
1189 foundM = true;
1190 }
1191 }
1192
1193 for ( const QString &pointCoordinates : coordList )
1194 {
1195#if QT_VERSION < QT_VERSION_CHECK(5, 15, 0)
1196 QStringList coordinates = pointCoordinates.split( rx, QString::SkipEmptyParts );
1197#else
1198 QStringList coordinates = pointCoordinates.split( rx, Qt::SkipEmptyParts );
1199#endif
1200 if ( coordinates.size() < dim )
1201 continue;
1202
1203 int idx = 0;
1204 const double x = coordinates[idx++].toDouble();
1205 const double y = coordinates[idx++].toDouble();
1206
1207 double z = 0;
1208 if ( ( is3D || foundZ ) && coordinates.length() > idx )
1209 z = coordinates[idx++].toDouble();
1210
1211 double m = 0;
1212 if ( ( isMeasure || foundM ) && coordinates.length() > idx )
1213 m = coordinates[idx++].toDouble();
1214
1216 if ( is3D || foundZ )
1217 {
1218 if ( isMeasure || foundM )
1220 else
1222 }
1223 else
1224 {
1225 if ( isMeasure || foundM )
1227 else
1229 }
1230
1231 points.append( QgsPoint( t, x, y, z, m ) );
1232 }
1233
1234 return points;
1235}
1236
1237void QgsGeometryUtils::pointsToWKB( QgsWkbPtr &wkb, const QgsPointSequence &points, bool is3D, bool isMeasure, QgsAbstractGeometry::WkbFlags flags )
1238{
1239 wkb << static_cast<quint32>( points.size() );
1240 for ( const QgsPoint &point : points )
1241 {
1242 wkb << point.x() << point.y();
1243 if ( is3D )
1244 {
1245 double z = point.z();
1247 && std::isnan( z ) )
1248 z = -std::numeric_limits<double>::max();
1249
1250 wkb << z;
1251 }
1252 if ( isMeasure )
1253 {
1254 double m = point.m();
1256 && std::isnan( m ) )
1257 m = -std::numeric_limits<double>::max();
1258
1259 wkb << m;
1260 }
1261 }
1262}
1263
1264QString QgsGeometryUtils::pointsToWKT( const QgsPointSequence &points, int precision, bool is3D, bool isMeasure )
1265{
1266 QString wkt = QStringLiteral( "(" );
1267 for ( const QgsPoint &p : points )
1268 {
1269 wkt += qgsDoubleToString( p.x(), precision );
1270 wkt += ' ' + qgsDoubleToString( p.y(), precision );
1271 if ( is3D )
1272 wkt += ' ' + qgsDoubleToString( p.z(), precision );
1273 if ( isMeasure )
1274 wkt += ' ' + qgsDoubleToString( p.m(), precision );
1275 wkt += QLatin1String( ", " );
1276 }
1277 if ( wkt.endsWith( QLatin1String( ", " ) ) )
1278 wkt.chop( 2 ); // Remove last ", "
1279 wkt += ')';
1280 return wkt;
1281}
1282
1283QDomElement QgsGeometryUtils::pointsToGML2( const QgsPointSequence &points, QDomDocument &doc, int precision, const QString &ns, QgsAbstractGeometry::AxisOrder axisOrder )
1284{
1285 QDomElement elemCoordinates = doc.createElementNS( ns, QStringLiteral( "coordinates" ) );
1286
1287 // coordinate separator
1288 const QString cs = QStringLiteral( "," );
1289 // tuple separator
1290 const QString ts = QStringLiteral( " " );
1291
1292 elemCoordinates.setAttribute( QStringLiteral( "cs" ), cs );
1293 elemCoordinates.setAttribute( QStringLiteral( "ts" ), ts );
1294
1295 QString strCoordinates;
1296
1297 for ( const QgsPoint &p : points )
1298 if ( axisOrder == QgsAbstractGeometry::AxisOrder::XY )
1299 strCoordinates += qgsDoubleToString( p.x(), precision ) + cs + qgsDoubleToString( p.y(), precision ) + ts;
1300 else
1301 strCoordinates += qgsDoubleToString( p.y(), precision ) + cs + qgsDoubleToString( p.x(), precision ) + ts;
1302
1303 if ( strCoordinates.endsWith( ts ) )
1304 strCoordinates.chop( 1 ); // Remove trailing space
1305
1306 elemCoordinates.appendChild( doc.createTextNode( strCoordinates ) );
1307 return elemCoordinates;
1308}
1309
1310QDomElement QgsGeometryUtils::pointsToGML3( const QgsPointSequence &points, QDomDocument &doc, int precision, const QString &ns, bool is3D, QgsAbstractGeometry::AxisOrder axisOrder )
1311{
1312 QDomElement elemPosList = doc.createElementNS( ns, QStringLiteral( "posList" ) );
1313 elemPosList.setAttribute( QStringLiteral( "srsDimension" ), is3D ? 3 : 2 );
1314
1315 QString strCoordinates;
1316 for ( const QgsPoint &p : points )
1317 {
1318 if ( axisOrder == QgsAbstractGeometry::AxisOrder::XY )
1319 strCoordinates += qgsDoubleToString( p.x(), precision ) + ' ' + qgsDoubleToString( p.y(), precision ) + ' ';
1320 else
1321 strCoordinates += qgsDoubleToString( p.y(), precision ) + ' ' + qgsDoubleToString( p.x(), precision ) + ' ';
1322 if ( is3D )
1323 strCoordinates += qgsDoubleToString( p.z(), precision ) + ' ';
1324 }
1325 if ( strCoordinates.endsWith( ' ' ) )
1326 strCoordinates.chop( 1 ); // Remove trailing space
1327
1328 elemPosList.appendChild( doc.createTextNode( strCoordinates ) );
1329 return elemPosList;
1330}
1331
1333{
1334 QString json = QStringLiteral( "[ " );
1335 for ( const QgsPoint &p : points )
1336 {
1337 json += '[' + qgsDoubleToString( p.x(), precision ) + QLatin1String( ", " ) + qgsDoubleToString( p.y(), precision ) + QLatin1String( "], " );
1338 }
1339 if ( json.endsWith( QLatin1String( ", " ) ) )
1340 {
1341 json.chop( 2 ); // Remove last ", "
1342 }
1343 json += ']';
1344 return json;
1345}
1346
1347
1349{
1350 json coordinates( json::array() );
1351 for ( const QgsPoint &p : points )
1352 {
1353 if ( p.is3D() )
1354 {
1355 coordinates.push_back( { qgsRound( p.x(), precision ), qgsRound( p.y(), precision ), qgsRound( p.z(), precision ) } );
1356 }
1357 else
1358 {
1359 coordinates.push_back( { qgsRound( p.x(), precision ), qgsRound( p.y(), precision ) } );
1360 }
1361 }
1362 return coordinates;
1363}
1364
1366{
1367 double clippedAngle = angle;
1368 if ( clippedAngle >= M_PI * 2 || clippedAngle <= -2 * M_PI )
1369 {
1370 clippedAngle = std::fmod( clippedAngle, 2 * M_PI );
1371 }
1372 if ( clippedAngle < 0.0 )
1373 {
1374 clippedAngle += 2 * M_PI;
1375 }
1376 return clippedAngle;
1377}
1378
1379QPair<Qgis::WkbType, QString> QgsGeometryUtils::wktReadBlock( const QString &wkt )
1380{
1381 QString wktParsed = wkt;
1382 QString contents;
1383 const QLatin1String empty { "EMPTY" };
1384 if ( wkt.contains( empty, Qt::CaseInsensitive ) )
1385 {
1386 const thread_local QRegularExpression whiteSpaces( "\\s" );
1387 wktParsed.remove( whiteSpaces );
1388 const int index = wktParsed.indexOf( empty, 0, Qt::CaseInsensitive );
1389
1390 if ( index == wktParsed.length() - empty.size() )
1391 {
1392 // "EMPTY" found at the end of the QString
1393 // Extract the part of the QString to the left of "EMPTY"
1394 wktParsed = wktParsed.left( index );
1395 contents = empty;
1396 }
1397 }
1398 else
1399 {
1400 const int openedParenthesisCount = wktParsed.count( '(' );
1401 const int closedParenthesisCount = wktParsed.count( ')' );
1402 // closes missing parentheses
1403 for ( int i = 0 ; i < openedParenthesisCount - closedParenthesisCount; ++i )
1404 wktParsed.push_back( ')' );
1405 // removes extra parentheses
1406 wktParsed.truncate( wktParsed.size() - ( closedParenthesisCount - openedParenthesisCount ) );
1407
1408 const thread_local QRegularExpression cooRegEx( QStringLiteral( "^[^\\(]*\\((.*)\\)[^\\)]*$" ), QRegularExpression::DotMatchesEverythingOption );
1409 const QRegularExpressionMatch match = cooRegEx.match( wktParsed );
1410 contents = match.hasMatch() ? match.captured( 1 ) : QString();
1411 }
1412 const Qgis::WkbType wkbType = QgsWkbTypes::parseType( wktParsed );
1413 return qMakePair( wkbType, contents );
1414}
1415
1416QStringList QgsGeometryUtils::wktGetChildBlocks( const QString &wkt, const QString &defaultType )
1417{
1418 int level = 0;
1419 QString block;
1420 block.reserve( wkt.size() );
1421 QStringList blocks;
1422
1423 const QChar *wktData = wkt.data();
1424 const int wktLength = wkt.length();
1425 for ( int i = 0, n = wktLength; i < n; ++i, ++wktData )
1426 {
1427 if ( ( wktData->isSpace() || *wktData == '\n' || *wktData == '\t' ) && level == 0 )
1428 continue;
1429
1430 if ( *wktData == ',' && level == 0 )
1431 {
1432 if ( !block.isEmpty() )
1433 {
1434 if ( block.startsWith( '(' ) && !defaultType.isEmpty() )
1435 block.prepend( defaultType + ' ' );
1436 blocks.append( block );
1437 }
1438 block.resize( 0 );
1439 continue;
1440 }
1441 if ( *wktData == '(' )
1442 ++level;
1443 else if ( *wktData == ')' )
1444 --level;
1445 block += *wktData;
1446 }
1447 if ( !block.isEmpty() )
1448 {
1449 if ( block.startsWith( '(' ) && !defaultType.isEmpty() )
1450 block.prepend( defaultType + ' ' );
1451 blocks.append( block );
1452 }
1453 return blocks;
1454}
1455
1456int QgsGeometryUtils::closestSideOfRectangle( double right, double bottom, double left, double top, double x, double y )
1457{
1458 // point outside rectangle
1459 if ( x <= left && y <= bottom )
1460 {
1461 const double dx = left - x;
1462 const double dy = bottom - y;
1463 if ( qgsDoubleNear( dx, dy ) )
1464 return 6;
1465 else if ( dx < dy )
1466 return 5;
1467 else
1468 return 7;
1469 }
1470 else if ( x >= right && y >= top )
1471 {
1472 const double dx = x - right;
1473 const double dy = y - top;
1474 if ( qgsDoubleNear( dx, dy ) )
1475 return 2;
1476 else if ( dx < dy )
1477 return 1;
1478 else
1479 return 3;
1480 }
1481 else if ( x >= right && y <= bottom )
1482 {
1483 const double dx = x - right;
1484 const double dy = bottom - y;
1485 if ( qgsDoubleNear( dx, dy ) )
1486 return 4;
1487 else if ( dx < dy )
1488 return 5;
1489 else
1490 return 3;
1491 }
1492 else if ( x <= left && y >= top )
1493 {
1494 const double dx = left - x;
1495 const double dy = y - top;
1496 if ( qgsDoubleNear( dx, dy ) )
1497 return 8;
1498 else if ( dx < dy )
1499 return 1;
1500 else
1501 return 7;
1502 }
1503 else if ( x <= left )
1504 return 7;
1505 else if ( x >= right )
1506 return 3;
1507 else if ( y <= bottom )
1508 return 5;
1509 else if ( y >= top )
1510 return 1;
1511
1512 // point is inside rectangle
1513 const double smallestX = std::min( right - x, x - left );
1514 const double smallestY = std::min( top - y, y - bottom );
1515 if ( smallestX < smallestY )
1516 {
1517 // closer to left/right side
1518 if ( right - x < x - left )
1519 return 3; // closest to right side
1520 else
1521 return 7;
1522 }
1523 else
1524 {
1525 // closer to top/bottom side
1526 if ( top - y < y - bottom )
1527 return 1; // closest to top side
1528 else
1529 return 5;
1530 }
1531}
1532
1534{
1536
1537
1538 const double x = ( pt1.x() + pt2.x() ) / 2.0;
1539 const double y = ( pt1.y() + pt2.y() ) / 2.0;
1540 double z = std::numeric_limits<double>::quiet_NaN();
1541 double m = std::numeric_limits<double>::quiet_NaN();
1542
1543 if ( pt1.is3D() || pt2.is3D() )
1544 {
1545 pType = QgsWkbTypes::addZ( pType );
1546 z = ( pt1.z() + pt2.z() ) / 2.0;
1547 }
1548
1549 if ( pt1.isMeasure() || pt2.isMeasure() )
1550 {
1551 pType = QgsWkbTypes::addM( pType );
1552 m = ( pt1.m() + pt2.m() ) / 2.0;
1553 }
1554
1555 return QgsPoint( pType, x, y, z, m );
1556}
1557
1558QgsPoint QgsGeometryUtils::interpolatePointOnLine( const QgsPoint &p1, const QgsPoint &p2, const double fraction )
1559{
1560 const double _fraction = 1 - fraction;
1561 return QgsPoint( p1.wkbType(),
1562 p1.x() * _fraction + p2.x() * fraction,
1563 p1.y() * _fraction + p2.y() * fraction,
1564 p1.is3D() ? p1.z() * _fraction + p2.z() * fraction : std::numeric_limits<double>::quiet_NaN(),
1565 p1.isMeasure() ? p1.m() * _fraction + p2.m() * fraction : std::numeric_limits<double>::quiet_NaN() );
1566}
1567
1568QgsPointXY QgsGeometryUtils::interpolatePointOnLine( const double x1, const double y1, const double x2, const double y2, const double fraction )
1569{
1570 const double deltaX = ( x2 - x1 ) * fraction;
1571 const double deltaY = ( y2 - y1 ) * fraction;
1572 return QgsPointXY( x1 + deltaX, y1 + deltaY );
1573}
1574
1575QgsPointXY QgsGeometryUtils::interpolatePointOnLineByValue( const double x1, const double y1, const double v1, const double x2, const double y2, const double v2, const double value )
1576{
1577 if ( qgsDoubleNear( v1, v2 ) )
1578 return QgsPointXY( x1, y1 );
1579
1580 const double fraction = ( value - v1 ) / ( v2 - v1 );
1581 return interpolatePointOnLine( x1, y1, x2, y2, fraction );
1582}
1583
1584double QgsGeometryUtils::gradient( const QgsPoint &pt1, const QgsPoint &pt2 )
1585{
1586 const double delta_x = pt2.x() - pt1.x();
1587 const double delta_y = pt2.y() - pt1.y();
1588 if ( qgsDoubleNear( delta_x, 0.0 ) )
1589 {
1590 return INFINITY;
1591 }
1592
1593 return delta_y / delta_x;
1594}
1595
1596void QgsGeometryUtils::coefficients( const QgsPoint &pt1, const QgsPoint &pt2, double &a, double &b, double &c )
1597{
1598 if ( qgsDoubleNear( pt1.x(), pt2.x() ) )
1599 {
1600 a = 1;
1601 b = 0;
1602 c = -pt1.x();
1603 }
1604 else if ( qgsDoubleNear( pt1.y(), pt2.y() ) )
1605 {
1606 a = 0;
1607 b = 1;
1608 c = -pt1.y();
1609 }
1610 else
1611 {
1612 a = pt1.y() - pt2.y();
1613 b = pt2.x() - pt1.x();
1614 c = pt1.x() * pt2.y() - pt1.y() * pt2.x();
1615 }
1616
1617}
1618
1620{
1621 QgsLineString line;
1622 QgsPoint p2;
1623
1624 if ( ( p == s1 ) || ( p == s2 ) )
1625 {
1626 return line;
1627 }
1628
1629 double a, b, c;
1630 coefficients( s1, s2, a, b, c );
1631
1632 if ( qgsDoubleNear( a, 0 ) )
1633 {
1634 p2 = QgsPoint( p.x(), s1.y() );
1635 }
1636 else if ( qgsDoubleNear( b, 0 ) )
1637 {
1638 p2 = QgsPoint( s1.x(), p.y() );
1639 }
1640 else
1641 {
1642 const double y = ( -c - a * p.x() ) / b;
1643 const double m = gradient( s1, s2 );
1644 const double d2 = 1 + m * m;
1645 const double H = p.y() - y;
1646 const double dx = m * H / d2;
1647 const double dy = m * dx;
1648 p2 = QgsPoint( p.x() + dx, y + dy );
1649 }
1650
1651 line.addVertex( p );
1652 line.addVertex( p2 );
1653
1654 return line;
1655}
1656
1657void QgsGeometryUtils::perpendicularCenterSegment( double pointx, double pointy, double segmentPoint1x, double segmentPoint1y, double segmentPoint2x, double segmentPoint2y, double &perpendicularSegmentPoint1x, double &perpendicularSegmentPoint1y, double &perpendicularSegmentPoint2x, double &perpendicularSegmentPoint2y, double desiredSegmentLength )
1658{
1659 QgsVector segmentVector = QgsVector( segmentPoint2x - segmentPoint1x, segmentPoint2y - segmentPoint1y );
1660 QgsVector perpendicularVector = segmentVector.perpVector();
1661 if ( desiredSegmentLength != 0 )
1662 {
1663 perpendicularVector = perpendicularVector.normalized() * ( desiredSegmentLength ) / 2;
1664 }
1665
1666 perpendicularSegmentPoint1x = pointx - perpendicularVector.x();
1667 perpendicularSegmentPoint1y = pointy - perpendicularVector.y();
1668 perpendicularSegmentPoint2x = pointx + perpendicularVector.x();
1669 perpendicularSegmentPoint2y = pointy + perpendicularVector.y();
1670}
1671
1672double QgsGeometryUtils::lineAngle( double x1, double y1, double x2, double y2 )
1673{
1674 const double at = std::atan2( y2 - y1, x2 - x1 );
1675 const double a = -at + M_PI_2;
1676 return normalizedAngle( a );
1677}
1678
1679double QgsGeometryUtils::angleBetweenThreePoints( double x1, double y1, double x2, double y2, double x3, double y3 )
1680{
1681 const double angle1 = std::atan2( y1 - y2, x1 - x2 );
1682 const double angle2 = std::atan2( y3 - y2, x3 - x2 );
1683 return normalizedAngle( angle1 - angle2 );
1684}
1685
1686double QgsGeometryUtils::linePerpendicularAngle( double x1, double y1, double x2, double y2 )
1687{
1688 double a = lineAngle( x1, y1, x2, y2 );
1689 a += M_PI_2;
1690 return normalizedAngle( a );
1691}
1692
1693double QgsGeometryUtils::averageAngle( double x1, double y1, double x2, double y2, double x3, double y3 )
1694{
1695 // calc average angle between the previous and next point
1696 const double a1 = lineAngle( x1, y1, x2, y2 );
1697 const double a2 = lineAngle( x2, y2, x3, y3 );
1698 return averageAngle( a1, a2 );
1699}
1700
1701double QgsGeometryUtils::averageAngle( double a1, double a2 )
1702{
1703 a1 = normalizedAngle( a1 );
1704 a2 = normalizedAngle( a2 );
1705 double clockwiseDiff = 0.0;
1706 if ( a2 >= a1 )
1707 {
1708 clockwiseDiff = a2 - a1;
1709 }
1710 else
1711 {
1712 clockwiseDiff = a2 + ( 2 * M_PI - a1 );
1713 }
1714 const double counterClockwiseDiff = 2 * M_PI - clockwiseDiff;
1715
1716 double resultAngle = 0;
1717 if ( clockwiseDiff <= counterClockwiseDiff )
1718 {
1719 resultAngle = a1 + clockwiseDiff / 2.0;
1720 }
1721 else
1722 {
1723 resultAngle = a1 - counterClockwiseDiff / 2.0;
1724 }
1725 return normalizedAngle( resultAngle );
1726}
1727
1729 const QgsVector3D &P2, const QgsVector3D &P22 )
1730{
1731 const QgsVector3D u1 = P12 - P1;
1732 const QgsVector3D u2 = P22 - P2;
1734 if ( u3.length() == 0 ) return 1;
1735 u3.normalize();
1736 const QgsVector3D dir = P1 - P2;
1737 return std::fabs( ( QgsVector3D::dotProduct( dir, u3 ) ) ); // u3 is already normalized
1738}
1739
1741 const QgsVector3D &P2, const QgsVector3D &P22,
1742 QgsVector3D &X1, double epsilon )
1743{
1744 const QgsVector3D d = P2 - P1;
1745 QgsVector3D u1 = P12 - P1;
1746 u1.normalize();
1747 QgsVector3D u2 = P22 - P2;
1748 u2.normalize();
1749 const QgsVector3D u3 = QgsVector3D::crossProduct( u1, u2 );
1750
1751 if ( std::fabs( u3.x() ) <= epsilon &&
1752 std::fabs( u3.y() ) <= epsilon &&
1753 std::fabs( u3.z() ) <= epsilon )
1754 {
1755 // The rays are almost parallel.
1756 return false;
1757 }
1758
1759 // X1 and X2 are the closest points on lines
1760 // we want to find X1 (lies on u1)
1761 // solving the linear equation in r1 and r2: Xi = Pi + ri*ui
1762 // we are only interested in X1 so we only solve for r1.
1763 float a1 = QgsVector3D::dotProduct( u1, u1 ), b1 = QgsVector3D::dotProduct( u1, u2 ), c1 = QgsVector3D::dotProduct( u1, d );
1764 float a2 = QgsVector3D::dotProduct( u1, u2 ), b2 = QgsVector3D::dotProduct( u2, u2 ), c2 = QgsVector3D::dotProduct( u2, d );
1765 if ( !( std::fabs( b1 ) > epsilon ) )
1766 {
1767 // Denominator is close to zero.
1768 return false;
1769 }
1770 if ( !( a2 != -1 && a2 != 1 ) )
1771 {
1772 // Lines are parallel
1773 return false;
1774 }
1775
1776 const double r1 = ( c2 - b2 * c1 / b1 ) / ( a2 - b2 * a1 / b1 );
1777 X1 = P1 + u1 * r1;
1778
1779 return true;
1780}
1781
1783 const QgsVector3D &Lb1, const QgsVector3D &Lb2,
1784 QgsVector3D &intersection )
1785{
1786
1787 // if all Vector are on the same plane (have the same Z), use the 2D intersection
1788 // else return a false result
1789 if ( qgsDoubleNear( La1.z(), La2.z() ) && qgsDoubleNear( La1.z(), Lb1.z() ) && qgsDoubleNear( La1.z(), Lb2.z() ) )
1790 {
1791 QgsPoint ptInter;
1792 bool isIntersection;
1793 segmentIntersection( QgsPoint( La1.x(), La1.y() ),
1794 QgsPoint( La2.x(), La2.y() ),
1795 QgsPoint( Lb1.x(), Lb1.y() ),
1796 QgsPoint( Lb2.x(), Lb2.y() ),
1797 ptInter,
1798 isIntersection,
1799 1e-8,
1800 true );
1801 intersection.set( ptInter.x(), ptInter.y(), La1.z() );
1802 return true;
1803 }
1804
1805 // first check if lines have an exact intersection point
1806 // do it by checking if the shortest distance is exactly 0
1807 const float distance = skewLinesDistance( La1, La2, Lb1, Lb2 );
1808 if ( qgsDoubleNear( distance, 0.0 ) )
1809 {
1810 // 3d lines have exact intersection point.
1811 const QgsVector3D C = La2;
1812 const QgsVector3D D = Lb2;
1813 const QgsVector3D e = La1 - La2;
1814 const QgsVector3D f = Lb1 - Lb2;
1815 const QgsVector3D g = D - C;
1816 if ( qgsDoubleNear( ( QgsVector3D::crossProduct( f, g ) ).length(), 0.0 ) || qgsDoubleNear( ( QgsVector3D::crossProduct( f, e ) ).length(), 0.0 ) )
1817 {
1818 // Lines have no intersection, are they parallel?
1819 return false;
1820 }
1821
1823 fgn.normalize();
1824
1826 fen.normalize();
1827
1828 int di = -1;
1829 if ( fgn == fen ) // same direction?
1830 di *= -1;
1831
1832 intersection = C + e * di * ( QgsVector3D::crossProduct( f, g ).length() / QgsVector3D::crossProduct( f, e ).length() );
1833 return true;
1834 }
1835
1836 // try to calculate the approximate intersection point
1837 QgsVector3D X1, X2;
1838 const bool firstIsDone = skewLinesProjection( La1, La2, Lb1, Lb2, X1 );
1839 const bool secondIsDone = skewLinesProjection( Lb1, Lb2, La1, La2, X2 );
1840
1841 if ( !firstIsDone || !secondIsDone )
1842 {
1843 // Could not obtain projection point.
1844 return false;
1845 }
1846
1847 intersection = ( X1 + X2 ) / 2.0;
1848 return true;
1849}
1850
1851double QgsGeometryUtils::triangleArea( double aX, double aY, double bX, double bY, double cX, double cY )
1852{
1853 return 0.5 * std::abs( ( aX - cX ) * ( bY - aY ) - ( aX - bX ) * ( cY - aY ) );
1854}
1855
1856double QgsGeometryUtils::pointFractionAlongLine( double x1, double y1, double x2, double y2, double px, double py )
1857{
1858 const double dxp = px - x1;
1859 const double dyp = py - y1;
1860
1861 const double dxl = x2 - x1;
1862 const double dyl = y2 - y1;
1863
1864 return std::sqrt( ( dxp * dxp ) + ( dyp * dyp ) ) / std::sqrt( ( dxl * dxl ) + ( dyl * dyl ) );
1865}
1866
1867void QgsGeometryUtils::weightedPointInTriangle( const double aX, const double aY, const double bX, const double bY, const double cX, const double cY,
1868 double weightB, double weightC, double &pointX, double &pointY )
1869{
1870 // if point will be outside of the triangle, invert weights
1871 if ( weightB + weightC > 1 )
1872 {
1873 weightB = 1 - weightB;
1874 weightC = 1 - weightC;
1875 }
1876
1877 const double rBx = weightB * ( bX - aX );
1878 const double rBy = weightB * ( bY - aY );
1879 const double rCx = weightC * ( cX - aX );
1880 const double rCy = weightC * ( cY - aY );
1881
1882 pointX = rBx + rCx + aX;
1883 pointY = rBy + rCy + aY;
1884}
1885
1886bool QgsGeometryUtils::pointsAreCollinear( double x1, double y1, double x2, double y2, double x3, double y3, double epsilon )
1887{
1888 return qgsDoubleNear( x1 * ( y2 - y3 ) + x2 * ( y3 - y1 ) + x3 * ( y1 - y2 ), 0, epsilon );
1889};
1890
1892{
1893 bool rc = false;
1894
1895 for ( const QgsPoint &pt : points )
1896 {
1897 if ( pt.isMeasure() )
1898 {
1899 point.convertTo( QgsWkbTypes::addM( point.wkbType() ) );
1900 point.setM( pt.m() );
1901 rc = true;
1902 break;
1903 }
1904 }
1905
1906 return rc;
1907}
1908
1910{
1911 bool rc = false;
1912
1913 for ( const QgsPoint &pt : points )
1914 {
1915 if ( pt.is3D() )
1916 {
1917 point.convertTo( QgsWkbTypes::addZ( point.wkbType() ) );
1918 point.setZ( pt.z() );
1919 rc = true;
1920 break;
1921 }
1922 }
1923
1924 return rc;
1925}
1926
1927bool QgsGeometryUtils::angleBisector( double aX, double aY, double bX, double bY, double cX, double cY, double dX, double dY,
1928 double &pointX SIP_OUT, double &pointY SIP_OUT, double &angle SIP_OUT )
1929{
1930 const QgsPoint pA = QgsPoint( aX, aY );
1931 const QgsPoint pB = QgsPoint( bX, bY );
1932 const QgsPoint pC = QgsPoint( cX, cY );
1933 const QgsPoint pD = QgsPoint( dX, dY );
1934 angle = ( pA.azimuth( pB ) + pC.azimuth( pD ) ) / 2.0;
1935
1936 QgsPoint pOut;
1937 bool intersection = false;
1938 QgsGeometryUtils::segmentIntersection( pA, pB, pC, pD, pOut, intersection );
1939
1940 pointX = pOut.x();
1941 pointY = pOut.y();
1942
1943 return intersection;
1944}
1945
1946bool QgsGeometryUtils::bisector( double aX, double aY, double bX, double bY, double cX, double cY,
1947 double &pointX SIP_OUT, double &pointY SIP_OUT )
1948{
1949 const QgsPoint pA = QgsPoint( aX, aY );
1950 const QgsPoint pB = QgsPoint( bX, bY );
1951 const QgsPoint pC = QgsPoint( cX, cY );
1952 const double angle = ( pA.azimuth( pB ) + pA.azimuth( pC ) ) / 2.0;
1953
1954 QgsPoint pOut;
1955 bool intersection = false;
1956 QgsGeometryUtils::segmentIntersection( pB, pC, pA, pA.project( 1.0, angle ), pOut, intersection );
1957
1958 pointX = pOut.x();
1959 pointY = pOut.y();
1960
1961 return intersection;
1962}
WkbType
The WKB type describes the number of dimensions a geometry has.
Definition qgis.h:182
@ PointM
PointM.
@ PointZ
PointZ.
@ PointZM
PointZM.
Abstract base class for all geometries.
SegmentationToleranceType
Segmentation tolerance as maximum angle or maximum difference between approximation and circle.
@ MaximumDifference
Maximum distance between an arbitrary point on the original curve and closest point on its approximat...
virtual int vertexCount(int part=0, int ring=0) const =0
Returns the number of vertices of which this geometry is built.
bool isMeasure() const
Returns true if the geometry contains m values.
bool is3D() const
Returns true if the geometry is 3D and contains a z-value.
AxisOrder
Axis order for GML generation.
@ XY
X comes before Y (or lon before lat)
virtual QgsPoint vertexAt(QgsVertexId id) const =0
Returns the point corresponding to a specified vertex id.
Qgis::WkbType wkbType() const
Returns the WKB type of the geometry.
virtual bool isEmpty() const
Returns true if the geometry is empty.
@ FlagExportNanAsDoubleMin
Use -DOUBLE_MAX to represent NaN (since QGIS 3.30)
virtual double segmentLength(QgsVertexId startVertex) const =0
Returns the length of the segment of the geometry which begins at startVertex.
virtual bool nextVertex(QgsVertexId &id, QgsPoint &vertex) const =0
Returns next vertex id and coordinates.
virtual double closestSegment(const QgsPoint &pt, QgsPoint &segmentPt, QgsVertexId &vertexAfter, int *leftOf=nullptr, double epsilon=4 *std::numeric_limits< double >::epsilon()) const =0
Searches for the closest segment of the geometry to a given point.
Curve polygon geometry type.
Abstract base class for curved geometry type.
Definition qgscurve.h:36
static int circleCircleIntersections(const QgsPointXY &center1, double radius1, const QgsPointXY &center2, double radius2, QgsPointXY &intersection1, QgsPointXY &intersection2)
Calculates the intersections points between the circle with center center1 and radius radius1 and the...
static bool linesIntersection3D(const QgsVector3D &La1, const QgsVector3D &La2, const QgsVector3D &Lb1, const QgsVector3D &Lb2, QgsVector3D &intersection)
An algorithm to calculate an (approximate) intersection of two lines in 3D.
static QString pointsToJSON(const QgsPointSequence &points, int precision)
Returns a geoJSON coordinates string.
static double sqrDistance2D(const QgsPoint &pt1, const QgsPoint &pt2)
Returns the squared 2D distance between two points.
static QgsPointXY interpolatePointOnLine(double x1, double y1, double x2, double y2, double fraction)
Interpolates the position of a point a fraction of the way along the line from (x1,...
static double circleTangentDirection(const QgsPoint &tangentPoint, const QgsPoint &cp1, const QgsPoint &cp2, const QgsPoint &cp3)
Calculates the direction angle of a circle tangent (clockwise from north in radians)
static double normalizedAngle(double angle)
Ensures that an angle is in the range 0 <= angle < 2 pi.
static QgsPoint pointOnLineWithDistance(const QgsPoint &startPoint, const QgsPoint &directionPoint, double distance)
Returns a point a specified distance toward a second point.
static double gradient(const QgsPoint &pt1, const QgsPoint &pt2)
Returns the gradient of a line defined by points pt1 and pt2.
static json pointsToJson(const QgsPointSequence &points, int precision)
Returns coordinates as json object.
static double interpolateArcValue(double angle, double a1, double a2, double a3, double zm1, double zm2, double zm3)
Interpolate a value at given angle on circular arc given values (zm1, zm2, zm3) at three different an...
static QVector< QgsLineString * > extractLineStrings(const QgsAbstractGeometry *geom)
Returns list of linestrings extracted from the passed geometry.
static bool lineIntersection(const QgsPoint &p1, QgsVector v1, const QgsPoint &p2, QgsVector v2, QgsPoint &intersection)
Computes the intersection between two lines.
static bool angleOnCircle(double angle, double angle1, double angle2, double angle3)
Returns true if an angle is between angle1 and angle3 on a circle described by angle1,...
static QVector< SelfIntersection > selfIntersections(const QgsAbstractGeometry *geom, int part, int ring, double tolerance)
Find self intersections in a polyline.
static bool transferFirstZValueToPoint(const QgsPointSequence &points, QgsPoint &point)
A Z dimension is added to point if one of the point in the list points is in 3D.
static bool segmentMidPoint(const QgsPoint &p1, const QgsPoint &p2, QgsPoint &result, double radius, const QgsPoint &mousePos)
Calculates midpoint on circle passing through p1 and p2, closest to the given coordinate mousePos.
static QgsPointXY interpolatePointOnLineByValue(double x1, double y1, double v1, double x2, double y2, double v2, double value)
Interpolates the position of a point along the line from (x1, y1) to (x2, y2).
static bool segmentIntersection(const QgsPoint &p1, const QgsPoint &p2, const QgsPoint &q1, const QgsPoint &q2, QgsPoint &intersectionPoint, bool &isIntersection, double tolerance=1e-8, bool acceptImproperIntersection=false)
Compute the intersection between two segments.
static QgsPoint closestPoint(const QgsAbstractGeometry &geometry, const QgsPoint &point)
Returns the nearest point on a segment of a geometry for the specified point.
static bool circleClockwise(double angle1, double angle2, double angle3)
Returns true if the circle defined by three angles is ordered clockwise.
static double pointFractionAlongLine(double x1, double y1, double x2, double y2, double px, double py)
Given the line (x1, y1) to (x2, y2) and a point (px, py) returns the fraction of the line length at w...
static double triangleArea(double aX, double aY, double bX, double bY, double cX, double cY)
Returns the area of the triangle denoted by the points (aX, aY), (bX, bY) and (cX,...
static double sweepAngle(double centerX, double centerY, double x1, double y1, double x2, double y2, double x3, double y3)
Calculates angle of a circular string part defined by pt1, pt2, pt3.
static bool verticesAtDistance(const QgsAbstractGeometry &geometry, double distance, QgsVertexId &previousVertex, QgsVertexId &nextVertex)
Retrieves the vertices which are before and after the interpolated point at a specified distance alon...
static QStringList wktGetChildBlocks(const QString &wkt, const QString &defaultType=QString())
Parses a WKT string and returns of list of blocks contained in the WKT.
static double circleLength(double x1, double y1, double x2, double y2, double x3, double y3)
Length of a circular string segment defined by pt1, pt2, pt3.
static QgsLineString perpendicularSegment(const QgsPoint &p, const QgsPoint &s1, const QgsPoint &s2)
Create a perpendicular line segment from p to segment [s1, s2].
static double sqrDistToLine(double ptX, double ptY, double x1, double y1, double x2, double y2, double &minDistX, double &minDistY, double epsilon)
Returns the squared distance between a point and a line.
static int segmentSide(const QgsPoint &pt1, const QgsPoint &pt3, const QgsPoint &pt2)
For line defined by points pt1 and pt3, find out on which side of the line is point pt3.
static void pointsToWKB(QgsWkbPtr &wkb, const QgsPointSequence &points, bool is3D, bool isMeasure, QgsAbstractGeometry::WkbFlags flags)
Returns a LinearRing { uint32 numPoints; Point points[numPoints]; }.
static double distanceToVertex(const QgsAbstractGeometry &geom, QgsVertexId id)
Returns the distance along a geometry from its first vertex to the specified vertex.
static void weightedPointInTriangle(double aX, double aY, double bX, double bY, double cX, double cY, double weightB, double weightC, double &pointX, double &pointY)
Returns a weighted point inside the triangle denoted by the points (aX, aY), (bX, bY) and (cX,...
static QPair< Qgis::WkbType, QString > wktReadBlock(const QString &wkt)
Parses a WKT block of the format "TYPE( contents )" and returns a pair of geometry type to contents (...
static bool lineCircleIntersection(const QgsPointXY &center, double radius, const QgsPointXY &linePoint1, const QgsPointXY &linePoint2, QgsPointXY &intersection)
Compute the intersection of a line and a circle.
static double distToInfiniteLine(const QgsPoint &point, const QgsPoint &linePoint1, const QgsPoint &linePoint2, double epsilon=1e-7)
Returns the distance between a point and an infinite line.
static double lineAngle(double x1, double y1, double x2, double y2)
Calculates the direction of line joining two points in radians, clockwise from the north direction.
static void perpendicularOffsetPointAlongSegment(double x1, double y1, double x2, double y2, double proportion, double offset, double *x, double *y)
Calculates a point a certain proportion of the way along the segment from (x1, y1) to (x2,...
static void coefficients(const QgsPoint &pt1, const QgsPoint &pt2, double &a, double &b, double &c)
Returns the coefficients (a, b, c for equation "ax + by + c = 0") of a line defined by points pt1 and...
static double linePerpendicularAngle(double x1, double y1, double x2, double y2)
Calculates the perpendicular angle to a line joining two points.
static bool transferFirstMValueToPoint(const QgsPointSequence &points, QgsPoint &point)
A M dimension is added to point if one of the points in the list points contains an M value.
static QgsPoint midpoint(const QgsPoint &pt1, const QgsPoint &pt2)
Returns a middle point between points pt1 and pt2.
static QgsPoint closestVertex(const QgsAbstractGeometry &geom, const QgsPoint &pt, QgsVertexId &id)
Returns the closest vertex to a geometry for a specified point.
static bool pointContinuesArc(const QgsPoint &a1, const QgsPoint &a2, const QgsPoint &a3, const QgsPoint &b, double distanceTolerance, double pointSpacingAngleTolerance)
Returns true if point b is on the arc formed by points a1, a2, and a3, but not within that arc portio...
static QgsPoint interpolatePointOnArc(const QgsPoint &pt1, const QgsPoint &pt2, const QgsPoint &pt3, double distance)
Interpolates a point on an arc defined by three points, pt1, pt2 and pt3.
static bool circleAngleBetween(double angle, double angle1, double angle2, bool clockwise)
Returns true if, in a circle, angle is between angle1 and angle2.
static double angleBetweenThreePoints(double x1, double y1, double x2, double y2, double x3, double y3)
Calculates the angle between the lines AB and BC, where AB and BC described by points a,...
static void perpendicularCenterSegment(double centerPointX, double centerPointY, double segmentPoint1x, double segmentPoint1y, double segmentPoint2x, double segmentPoint2y, double &perpendicularSegmentPoint1x, double &perpendicularSegmentPoint1y, double &perpendicularSegmentPoint2x, double &perpendicularSegmentPoint2y, double segmentLength=0)
Create a perpendicular line segment to a given segment [segmentPoint1,segmentPoint2] with its center ...
static bool tangentPointAndCircle(const QgsPointXY &center, double radius, const QgsPointXY &p, QgsPointXY &pt1, QgsPointXY &pt2)
Calculates the tangent points between the circle with the specified center and radius and the point p...
static bool bisector(double aX, double aY, double bX, double bY, double cX, double cY, double &pointX, double &pointY)
Returns the point (pointX, pointY) forming the bisector from point (aX, aY) to the segment (bX,...
static double ccwAngle(double dy, double dx)
Returns the counter clockwise angle between a line with components dx, dy and the line with dx > 0 an...
static int circleCircleOuterTangents(const QgsPointXY &center1, double radius1, const QgsPointXY &center2, double radius2, QgsPointXY &line1P1, QgsPointXY &line1P2, QgsPointXY &line2P1, QgsPointXY &line2P2)
Calculates the outer tangent points for two circles, centered at center1 and center2 and with radii o...
static int closestSideOfRectangle(double right, double bottom, double left, double top, double x, double y)
Returns a number representing the closest side of a rectangle defined by /a right,...
static double skewLinesDistance(const QgsVector3D &P1, const QgsVector3D &P12, const QgsVector3D &P2, const QgsVector3D &P22)
An algorithm to calculate the shortest distance between two skew lines.
static QgsPointSequence pointsFromWKT(const QString &wktCoordinateList, bool is3D, bool isMeasure)
Returns a list of points contained in a WKT string.
static void segmentizeArc(const QgsPoint &p1, const QgsPoint &p2, const QgsPoint &p3, QgsPointSequence &points, double tolerance=M_PI_2/90, QgsAbstractGeometry::SegmentationToleranceType toleranceType=QgsAbstractGeometry::MaximumAngle, bool hasZ=false, bool hasM=false)
Convert circular arc defined by p1, p2, p3 (p1/p3 being start resp.
static QDomElement pointsToGML2(const QgsPointSequence &points, QDomDocument &doc, int precision, const QString &ns, QgsAbstractGeometry::AxisOrder axisOrder=QgsAbstractGeometry::AxisOrder::XY)
Returns a gml::coordinates DOM element.
static bool transferFirstZOrMValueToPoint(Iterator verticesBegin, Iterator verticesEnd, QgsPoint &point)
A Z or M dimension is added to point if one of the points in the list points contains Z or M value.
static QDomElement pointsToGML3(const QgsPointSequence &points, QDomDocument &doc, int precision, const QString &ns, bool is3D, QgsAbstractGeometry::AxisOrder axisOrder=QgsAbstractGeometry::AxisOrder::XY)
Returns a gml::posList DOM element.
static int circleCircleInnerTangents(const QgsPointXY &center1, double radius1, const QgsPointXY &center2, double radius2, QgsPointXY &line1P1, QgsPointXY &line1P2, QgsPointXY &line2P1, QgsPointXY &line2P2)
Calculates the inner tangent points for two circles, centered at center1 and center2 and with radii o...
static bool angleBisector(double aX, double aY, double bX, double bY, double cX, double cY, double dX, double dY, double &pointX, double &pointY, double &angle)
Returns the point (pointX, pointY) forming the bisector from segment (aX aY) (bX bY) and segment (bX,...
static bool skewLinesProjection(const QgsVector3D &P1, const QgsVector3D &P12, const QgsVector3D &P2, const QgsVector3D &P22, QgsVector3D &X1, double epsilon=0.0001)
A method to project one skew line onto another.
static bool pointsAreCollinear(double x1, double y1, double x2, double y2, double x3, double y3, double epsilon)
Given the points (x1, y1), (x2, y2) and (x3, y3) returns true if these points can be considered colli...
static QString pointsToWKT(const QgsPointSequence &points, int precision, bool is3D, bool isMeasure)
Returns a WKT coordinate list.
static QgsPoint segmentMidPointFromCenter(const QgsPoint &p1, const QgsPoint &p2, const QgsPoint &center, bool useShortestArc=true)
Calculates the midpoint on the circle passing through p1 and p2, with the specified center coordinate...
static int leftOfLine(const double x, const double y, const double x1, const double y1, const double x2, const double y2)
Returns a value < 0 if the point (x, y) is left of the line from (x1, y1) -> (x2, y2).
static double averageAngle(double x1, double y1, double x2, double y2, double x3, double y3)
Calculates the average angle (in radians) between the two linear segments from (x1,...
static void circleCenterRadius(const QgsPoint &pt1, const QgsPoint &pt2, const QgsPoint &pt3, double &radius, double &centerX, double &centerY)
Returns radius and center of the circle through pt1, pt2, pt3.
Line string geometry type, with support for z-dimension and m-values.
void addVertex(const QgsPoint &pt)
Adds a new vertex to the end of the line string.
A class to represent a 2D point.
Definition qgspointxy.h:59
double sqrDist(double x, double y) const
Returns the squared distance between this point a specified x, y coordinate.
Definition qgspointxy.h:190
double distance(double x, double y) const
Returns the distance between this point and a specified x, y coordinate.
Definition qgspointxy.h:211
void set(double x, double y)
Sets the x and y value of the point.
Definition qgspointxy.h:139
double y
Definition qgspointxy.h:63
double x
Definition qgspointxy.h:62
Point geometry type, with support for z-dimension and m-values.
Definition qgspoint.h:49
bool addMValue(double mValue=0) override
Adds a measure to the geometry, initialized to a preset value.
Definition qgspoint.cpp:568
double azimuth(const QgsPoint &other) const
Calculates Cartesian azimuth between this point and other one (clockwise in degree,...
Definition qgspoint.cpp:722
bool addZValue(double zValue=0) override
Adds a z-dimension to the geometry, initialized to a preset value.
Definition qgspoint.cpp:557
double z
Definition qgspoint.h:54
double x
Definition qgspoint.h:52
void setM(double m)
Sets the point's m-value.
Definition qgspoint.h:319
bool convertTo(Qgis::WkbType type) override
Converts the geometry to a specified type.
Definition qgspoint.cpp:626
bool isEmpty() const override
Returns true if the geometry is empty.
Definition qgspoint.cpp:773
double distance(double x, double y) const
Returns the Cartesian 2D distance between this point and a specified x, y coordinate.
Definition qgspoint.h:343
void setZ(double z)
Sets the point's z-coordinate.
Definition qgspoint.h:304
double m
Definition qgspoint.h:55
QgsPoint project(double distance, double azimuth, double inclination=90.0) const
Returns a new point which corresponds to this point projected by a specified distance with specified ...
Definition qgspoint.cpp:741
double y
Definition qgspoint.h:53
Class for storage of 3D vectors similar to QVector3D, with the difference that it uses double precisi...
Definition qgsvector3d.h:32
double y() const
Returns Y coordinate.
Definition qgsvector3d.h:51
double z() const
Returns Z coordinate.
Definition qgsvector3d.h:53
static double dotProduct(const QgsVector3D &v1, const QgsVector3D &v2)
Returns the dot product of two vectors.
double x() const
Returns X coordinate.
Definition qgsvector3d.h:49
void normalize()
Normalizes the current vector in place.
static QgsVector3D crossProduct(const QgsVector3D &v1, const QgsVector3D &v2)
Returns the cross product of two vectors.
void set(double x, double y, double z)
Sets vector coordinates.
Definition qgsvector3d.h:74
double length() const
Returns the length of the vector.
A class to represent a vector.
Definition qgsvector.h:30
double y() const
Returns the vector's y-component.
Definition qgsvector.h:156
QgsVector normalized() const
Returns the vector's normalized (or "unit") vector (ie same angle but length of 1....
Definition qgsvector.cpp:28
QgsVector perpVector() const
Returns the perpendicular vector to this vector (rotated 90 degrees counter-clockwise)
Definition qgsvector.h:164
double x() const
Returns the vector's x-component.
Definition qgsvector.h:147
double length() const
Returns the length of the vector.
Definition qgsvector.h:128
WKB pointer handler.
Definition qgswkbptr.h:44
static Qgis::WkbType parseType(const QString &wktStr)
Attempts to extract the WKB type from a WKT string.
static Qgis::WkbType addM(Qgis::WkbType type)
Adds the m dimension to a WKB type and returns the new type.
static Qgis::WkbType addZ(Qgis::WkbType type)
Adds the z dimension to a WKB type and returns the new type.
static bool hasZ(Qgis::WkbType type)
Tests whether a WKB type contains the z-dimension.
static bool hasM(Qgis::WkbType type)
Tests whether a WKB type contains m values.
As part of the API refactoring and improvements which landed in the Processing API was substantially reworked from the x version This was done in order to allow much of the underlying Processing framework to be ported into c
QString qgsDoubleToString(double a, int precision=17)
Returns a string representation of a double.
Definition qgis.h:4271
double qgsRound(double number, int places)
Returns a double number, rounded (as close as possible) to the specified number of places.
Definition qgis.h:4386
bool qgsDoubleNear(double a, double b, double epsilon=4 *std::numeric_limits< double >::epsilon())
Compare two doubles (but allow some difference)
Definition qgis.h:4332
const double DEFAULT_SEGMENT_EPSILON
Default snapping tolerance for segments.
Definition qgis.h:4875
#define SIP_OUT
Definition qgis_sip.h:58
QVector< QgsPoint > QgsPointSequence
bool isClockwise(std::array< Direction, 4 > dirs)
Checks whether the 4 directions in dirs make up a clockwise rectangle.
int precision
Utility class for identifying a unique vertex within a geometry.
Definition qgsvertexid.h:31
int vertex
Vertex number.
Definition qgsvertexid.h:95
bool isValid() const
Returns true if the vertex id is valid.
Definition qgsvertexid.h:46
int part
Part number.
Definition qgsvertexid.h:89
Qgis::VertexType type
Vertex type.
Definition qgsvertexid.h:98
int ring
Ring number.
Definition qgsvertexid.h:92