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