QGIS API Documentation 3.28.0-Firenze (ed3ad0430f)
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 )
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 wkb << point.z();
1228 }
1229 if ( isMeasure )
1230 {
1231 wkb << point.m();
1232 }
1233 }
1234}
1235
1236QString QgsGeometryUtils::pointsToWKT( const QgsPointSequence &points, int precision, bool is3D, bool isMeasure )
1237{
1238 QString wkt = QStringLiteral( "(" );
1239 for ( const QgsPoint &p : points )
1240 {
1241 wkt += qgsDoubleToString( p.x(), precision );
1242 wkt += ' ' + qgsDoubleToString( p.y(), precision );
1243 if ( is3D )
1244 wkt += ' ' + qgsDoubleToString( p.z(), precision );
1245 if ( isMeasure )
1246 wkt += ' ' + qgsDoubleToString( p.m(), precision );
1247 wkt += QLatin1String( ", " );
1248 }
1249 if ( wkt.endsWith( QLatin1String( ", " ) ) )
1250 wkt.chop( 2 ); // Remove last ", "
1251 wkt += ')';
1252 return wkt;
1253}
1254
1255QDomElement QgsGeometryUtils::pointsToGML2( const QgsPointSequence &points, QDomDocument &doc, int precision, const QString &ns, QgsAbstractGeometry::AxisOrder axisOrder )
1256{
1257 QDomElement elemCoordinates = doc.createElementNS( ns, QStringLiteral( "coordinates" ) );
1258
1259 // coordinate separator
1260 const QString cs = QStringLiteral( "," );
1261 // tuple separator
1262 const QString ts = QStringLiteral( " " );
1263
1264 elemCoordinates.setAttribute( QStringLiteral( "cs" ), cs );
1265 elemCoordinates.setAttribute( QStringLiteral( "ts" ), ts );
1266
1267 QString strCoordinates;
1268
1269 for ( const QgsPoint &p : points )
1270 if ( axisOrder == QgsAbstractGeometry::AxisOrder::XY )
1271 strCoordinates += qgsDoubleToString( p.x(), precision ) + cs + qgsDoubleToString( p.y(), precision ) + ts;
1272 else
1273 strCoordinates += qgsDoubleToString( p.y(), precision ) + cs + qgsDoubleToString( p.x(), precision ) + ts;
1274
1275 if ( strCoordinates.endsWith( ts ) )
1276 strCoordinates.chop( 1 ); // Remove trailing space
1277
1278 elemCoordinates.appendChild( doc.createTextNode( strCoordinates ) );
1279 return elemCoordinates;
1280}
1281
1282QDomElement QgsGeometryUtils::pointsToGML3( const QgsPointSequence &points, QDomDocument &doc, int precision, const QString &ns, bool is3D, QgsAbstractGeometry::AxisOrder axisOrder )
1283{
1284 QDomElement elemPosList = doc.createElementNS( ns, QStringLiteral( "posList" ) );
1285 elemPosList.setAttribute( QStringLiteral( "srsDimension" ), is3D ? 3 : 2 );
1286
1287 QString strCoordinates;
1288 for ( const QgsPoint &p : points )
1289 {
1290 if ( axisOrder == QgsAbstractGeometry::AxisOrder::XY )
1291 strCoordinates += qgsDoubleToString( p.x(), precision ) + ' ' + qgsDoubleToString( p.y(), precision ) + ' ';
1292 else
1293 strCoordinates += qgsDoubleToString( p.y(), precision ) + ' ' + qgsDoubleToString( p.x(), precision ) + ' ';
1294 if ( is3D )
1295 strCoordinates += qgsDoubleToString( p.z(), precision ) + ' ';
1296 }
1297 if ( strCoordinates.endsWith( ' ' ) )
1298 strCoordinates.chop( 1 ); // Remove trailing space
1299
1300 elemPosList.appendChild( doc.createTextNode( strCoordinates ) );
1301 return elemPosList;
1302}
1303
1305{
1306 QString json = QStringLiteral( "[ " );
1307 for ( const QgsPoint &p : points )
1308 {
1309 json += '[' + qgsDoubleToString( p.x(), precision ) + QLatin1String( ", " ) + qgsDoubleToString( p.y(), precision ) + QLatin1String( "], " );
1310 }
1311 if ( json.endsWith( QLatin1String( ", " ) ) )
1312 {
1313 json.chop( 2 ); // Remove last ", "
1314 }
1315 json += ']';
1316 return json;
1317}
1318
1319
1321{
1322 json coordinates( json::array() );
1323 for ( const QgsPoint &p : points )
1324 {
1325 if ( p.is3D() )
1326 {
1327 coordinates.push_back( { qgsRound( p.x(), precision ), qgsRound( p.y(), precision ), qgsRound( p.z(), precision ) } );
1328 }
1329 else
1330 {
1331 coordinates.push_back( { qgsRound( p.x(), precision ), qgsRound( p.y(), precision ) } );
1332 }
1333 }
1334 return coordinates;
1335}
1336
1338{
1339 double clippedAngle = angle;
1340 if ( clippedAngle >= M_PI * 2 || clippedAngle <= -2 * M_PI )
1341 {
1342 clippedAngle = std::fmod( clippedAngle, 2 * M_PI );
1343 }
1344 if ( clippedAngle < 0.0 )
1345 {
1346 clippedAngle += 2 * M_PI;
1347 }
1348 return clippedAngle;
1349}
1350
1351QPair<QgsWkbTypes::Type, QString> QgsGeometryUtils::wktReadBlock( const QString &wkt )
1352{
1353 QString wktParsed = wkt;
1354 QString contents;
1355 if ( wkt.contains( QLatin1String( "EMPTY" ), Qt::CaseInsensitive ) )
1356 {
1357 const thread_local QRegularExpression sWktRegEx( QStringLiteral( "^\\s*(\\w+)\\s+(\\w+)\\s*$" ), QRegularExpression::DotMatchesEverythingOption );
1358 const QRegularExpressionMatch match = sWktRegEx.match( wkt );
1359 if ( match.hasMatch() )
1360 {
1361 wktParsed = match.captured( 1 );
1362 contents = match.captured( 2 ).toUpper();
1363 }
1364 }
1365 else
1366 {
1367 const int openedParenthesisCount = wktParsed.count( '(' );
1368 const int closedParenthesisCount = wktParsed.count( ')' );
1369 // closes missing parentheses
1370 for ( int i = 0 ; i < openedParenthesisCount - closedParenthesisCount; ++i )
1371 wktParsed.push_back( ')' );
1372 // removes extra parentheses
1373 wktParsed.truncate( wktParsed.size() - ( closedParenthesisCount - openedParenthesisCount ) );
1374
1375 const thread_local QRegularExpression cooRegEx( QStringLiteral( "^[^\\(]*\\((.*)\\)[^\\)]*$" ), QRegularExpression::DotMatchesEverythingOption );
1376 const QRegularExpressionMatch match = cooRegEx.match( wktParsed );
1377 contents = match.hasMatch() ? match.captured( 1 ) : QString();
1378 }
1379 const QgsWkbTypes::Type wkbType = QgsWkbTypes::parseType( wktParsed );
1380 return qMakePair( wkbType, contents );
1381}
1382
1383QStringList QgsGeometryUtils::wktGetChildBlocks( const QString &wkt, const QString &defaultType )
1384{
1385 int level = 0;
1386 QString block;
1387 block.reserve( wkt.size() );
1388 QStringList blocks;
1389
1390 const QChar *wktData = wkt.data();
1391 const int wktLength = wkt.length();
1392 for ( int i = 0, n = wktLength; i < n; ++i, ++wktData )
1393 {
1394 if ( ( wktData->isSpace() || *wktData == '\n' || *wktData == '\t' ) && level == 0 )
1395 continue;
1396
1397 if ( *wktData == ',' && level == 0 )
1398 {
1399 if ( !block.isEmpty() )
1400 {
1401 if ( block.startsWith( '(' ) && !defaultType.isEmpty() )
1402 block.prepend( defaultType + ' ' );
1403 blocks.append( block );
1404 }
1405 block.resize( 0 );
1406 continue;
1407 }
1408 if ( *wktData == '(' )
1409 ++level;
1410 else if ( *wktData == ')' )
1411 --level;
1412 block += *wktData;
1413 }
1414 if ( !block.isEmpty() )
1415 {
1416 if ( block.startsWith( '(' ) && !defaultType.isEmpty() )
1417 block.prepend( defaultType + ' ' );
1418 blocks.append( block );
1419 }
1420 return blocks;
1421}
1422
1423int QgsGeometryUtils::closestSideOfRectangle( double right, double bottom, double left, double top, double x, double y )
1424{
1425 // point outside rectangle
1426 if ( x <= left && y <= bottom )
1427 {
1428 const double dx = left - x;
1429 const double dy = bottom - y;
1430 if ( qgsDoubleNear( dx, dy ) )
1431 return 6;
1432 else if ( dx < dy )
1433 return 5;
1434 else
1435 return 7;
1436 }
1437 else if ( x >= right && y >= top )
1438 {
1439 const double dx = x - right;
1440 const double dy = y - top;
1441 if ( qgsDoubleNear( dx, dy ) )
1442 return 2;
1443 else if ( dx < dy )
1444 return 1;
1445 else
1446 return 3;
1447 }
1448 else if ( x >= right && y <= bottom )
1449 {
1450 const double dx = x - right;
1451 const double dy = bottom - y;
1452 if ( qgsDoubleNear( dx, dy ) )
1453 return 4;
1454 else if ( dx < dy )
1455 return 5;
1456 else
1457 return 3;
1458 }
1459 else if ( x <= left && y >= top )
1460 {
1461 const double dx = left - x;
1462 const double dy = y - top;
1463 if ( qgsDoubleNear( dx, dy ) )
1464 return 8;
1465 else if ( dx < dy )
1466 return 1;
1467 else
1468 return 7;
1469 }
1470 else if ( x <= left )
1471 return 7;
1472 else if ( x >= right )
1473 return 3;
1474 else if ( y <= bottom )
1475 return 5;
1476 else if ( y >= top )
1477 return 1;
1478
1479 // point is inside rectangle
1480 const double smallestX = std::min( right - x, x - left );
1481 const double smallestY = std::min( top - y, y - bottom );
1482 if ( smallestX < smallestY )
1483 {
1484 // closer to left/right side
1485 if ( right - x < x - left )
1486 return 3; // closest to right side
1487 else
1488 return 7;
1489 }
1490 else
1491 {
1492 // closer to top/bottom side
1493 if ( top - y < y - bottom )
1494 return 1; // closest to top side
1495 else
1496 return 5;
1497 }
1498}
1499
1501{
1503
1504
1505 const double x = ( pt1.x() + pt2.x() ) / 2.0;
1506 const double y = ( pt1.y() + pt2.y() ) / 2.0;
1507 double z = std::numeric_limits<double>::quiet_NaN();
1508 double m = std::numeric_limits<double>::quiet_NaN();
1509
1510 if ( pt1.is3D() || pt2.is3D() )
1511 {
1512 pType = QgsWkbTypes::addZ( pType );
1513 z = ( pt1.z() + pt2.z() ) / 2.0;
1514 }
1515
1516 if ( pt1.isMeasure() || pt2.isMeasure() )
1517 {
1518 pType = QgsWkbTypes::addM( pType );
1519 m = ( pt1.m() + pt2.m() ) / 2.0;
1520 }
1521
1522 return QgsPoint( pType, x, y, z, m );
1523}
1524
1525QgsPoint QgsGeometryUtils::interpolatePointOnLine( const QgsPoint &p1, const QgsPoint &p2, const double fraction )
1526{
1527 const double _fraction = 1 - fraction;
1528 return QgsPoint( p1.wkbType(),
1529 p1.x() * _fraction + p2.x() * fraction,
1530 p1.y() * _fraction + p2.y() * fraction,
1531 p1.is3D() ? p1.z() * _fraction + p2.z() * fraction : std::numeric_limits<double>::quiet_NaN(),
1532 p1.isMeasure() ? p1.m() * _fraction + p2.m() * fraction : std::numeric_limits<double>::quiet_NaN() );
1533}
1534
1535QgsPointXY QgsGeometryUtils::interpolatePointOnLine( const double x1, const double y1, const double x2, const double y2, const double fraction )
1536{
1537 const double deltaX = ( x2 - x1 ) * fraction;
1538 const double deltaY = ( y2 - y1 ) * fraction;
1539 return QgsPointXY( x1 + deltaX, y1 + deltaY );
1540}
1541
1542QgsPointXY QgsGeometryUtils::interpolatePointOnLineByValue( const double x1, const double y1, const double v1, const double x2, const double y2, const double v2, const double value )
1543{
1544 if ( qgsDoubleNear( v1, v2 ) )
1545 return QgsPointXY( x1, y1 );
1546
1547 const double fraction = ( value - v1 ) / ( v2 - v1 );
1548 return interpolatePointOnLine( x1, y1, x2, y2, fraction );
1549}
1550
1551double QgsGeometryUtils::gradient( const QgsPoint &pt1, const QgsPoint &pt2 )
1552{
1553 const double delta_x = pt2.x() - pt1.x();
1554 const double delta_y = pt2.y() - pt1.y();
1555 if ( qgsDoubleNear( delta_x, 0.0 ) )
1556 {
1557 return INFINITY;
1558 }
1559
1560 return delta_y / delta_x;
1561}
1562
1563void QgsGeometryUtils::coefficients( const QgsPoint &pt1, const QgsPoint &pt2, double &a, double &b, double &c )
1564{
1565 if ( qgsDoubleNear( pt1.x(), pt2.x() ) )
1566 {
1567 a = 1;
1568 b = 0;
1569 c = -pt1.x();
1570 }
1571 else if ( qgsDoubleNear( pt1.y(), pt2.y() ) )
1572 {
1573 a = 0;
1574 b = 1;
1575 c = -pt1.y();
1576 }
1577 else
1578 {
1579 a = pt1.y() - pt2.y();
1580 b = pt2.x() - pt1.x();
1581 c = pt1.x() * pt2.y() - pt1.y() * pt2.x();
1582 }
1583
1584}
1585
1587{
1588 QgsLineString line;
1589 QgsPoint p2;
1590
1591 if ( ( p == s1 ) || ( p == s2 ) )
1592 {
1593 return line;
1594 }
1595
1596 double a, b, c;
1597 coefficients( s1, s2, a, b, c );
1598
1599 if ( qgsDoubleNear( a, 0 ) )
1600 {
1601 p2 = QgsPoint( p.x(), s1.y() );
1602 }
1603 else if ( qgsDoubleNear( b, 0 ) )
1604 {
1605 p2 = QgsPoint( s1.x(), p.y() );
1606 }
1607 else
1608 {
1609 const double y = ( -c - a * p.x() ) / b;
1610 const double m = gradient( s1, s2 );
1611 const double d2 = 1 + m * m;
1612 const double H = p.y() - y;
1613 const double dx = m * H / d2;
1614 const double dy = m * dx;
1615 p2 = QgsPoint( p.x() + dx, y + dy );
1616 }
1617
1618 line.addVertex( p );
1619 line.addVertex( p2 );
1620
1621 return line;
1622}
1623
1624void QgsGeometryUtils::perpendicularCenterSegment( double pointx, double pointy, double segmentPoint1x, double segmentPoint1y, double segmentPoint2x, double segmentPoint2y, double &perpendicularSegmentPoint1x, double &perpendicularSegmentPoint1y, double &perpendicularSegmentPoint2x, double &perpendicularSegmentPoint2y, double desiredSegmentLength )
1625{
1626 QgsVector segmentVector = QgsVector( segmentPoint2x - segmentPoint1x, segmentPoint2y - segmentPoint1y );
1627 QgsVector perpendicularVector = segmentVector.perpVector();
1628 if ( desiredSegmentLength )
1629 {
1630 if ( desiredSegmentLength != 0 )
1631 {
1632 perpendicularVector = perpendicularVector.normalized() * ( desiredSegmentLength ) / 2;
1633 }
1634 }
1635 perpendicularSegmentPoint1x = pointx - perpendicularVector.x();
1636 perpendicularSegmentPoint1y = pointy - perpendicularVector.y();
1637 perpendicularSegmentPoint2x = pointx + perpendicularVector.x();
1638 perpendicularSegmentPoint2y = pointy + perpendicularVector.y();
1639}
1640
1641double QgsGeometryUtils::lineAngle( double x1, double y1, double x2, double y2 )
1642{
1643 const double at = std::atan2( y2 - y1, x2 - x1 );
1644 const double a = -at + M_PI_2;
1645 return normalizedAngle( a );
1646}
1647
1648double QgsGeometryUtils::angleBetweenThreePoints( double x1, double y1, double x2, double y2, double x3, double y3 )
1649{
1650 const double angle1 = std::atan2( y1 - y2, x1 - x2 );
1651 const double angle2 = std::atan2( y3 - y2, x3 - x2 );
1652 return normalizedAngle( angle1 - angle2 );
1653}
1654
1655double QgsGeometryUtils::linePerpendicularAngle( double x1, double y1, double x2, double y2 )
1656{
1657 double a = lineAngle( x1, y1, x2, y2 );
1658 a += M_PI_2;
1659 return normalizedAngle( a );
1660}
1661
1662double QgsGeometryUtils::averageAngle( double x1, double y1, double x2, double y2, double x3, double y3 )
1663{
1664 // calc average angle between the previous and next point
1665 const double a1 = lineAngle( x1, y1, x2, y2 );
1666 const double a2 = lineAngle( x2, y2, x3, y3 );
1667 return averageAngle( a1, a2 );
1668}
1669
1670double QgsGeometryUtils::averageAngle( double a1, double a2 )
1671{
1672 a1 = normalizedAngle( a1 );
1673 a2 = normalizedAngle( a2 );
1674 double clockwiseDiff = 0.0;
1675 if ( a2 >= a1 )
1676 {
1677 clockwiseDiff = a2 - a1;
1678 }
1679 else
1680 {
1681 clockwiseDiff = a2 + ( 2 * M_PI - a1 );
1682 }
1683 const double counterClockwiseDiff = 2 * M_PI - clockwiseDiff;
1684
1685 double resultAngle = 0;
1686 if ( clockwiseDiff <= counterClockwiseDiff )
1687 {
1688 resultAngle = a1 + clockwiseDiff / 2.0;
1689 }
1690 else
1691 {
1692 resultAngle = a1 - counterClockwiseDiff / 2.0;
1693 }
1694 return normalizedAngle( resultAngle );
1695}
1696
1698 const QgsVector3D &P2, const QgsVector3D &P22 )
1699{
1700 const QgsVector3D u1 = P12 - P1;
1701 const QgsVector3D u2 = P22 - P2;
1703 if ( u3.length() == 0 ) return 1;
1704 u3.normalize();
1705 const QgsVector3D dir = P1 - P2;
1706 return std::fabs( ( QgsVector3D::dotProduct( dir, u3 ) ) ); // u3 is already normalized
1707}
1708
1710 const QgsVector3D &P2, const QgsVector3D &P22,
1711 QgsVector3D &X1, double epsilon )
1712{
1713 const QgsVector3D d = P2 - P1;
1714 QgsVector3D u1 = P12 - P1;
1715 u1.normalize();
1716 QgsVector3D u2 = P22 - P2;
1717 u2.normalize();
1718 const QgsVector3D u3 = QgsVector3D::crossProduct( u1, u2 );
1719
1720 if ( std::fabs( u3.x() ) <= epsilon &&
1721 std::fabs( u3.y() ) <= epsilon &&
1722 std::fabs( u3.z() ) <= epsilon )
1723 {
1724 // The rays are almost parallel.
1725 return false;
1726 }
1727
1728 // X1 and X2 are the closest points on lines
1729 // we want to find X1 (lies on u1)
1730 // solving the linear equation in r1 and r2: Xi = Pi + ri*ui
1731 // we are only interested in X1 so we only solve for r1.
1732 float a1 = QgsVector3D::dotProduct( u1, u1 ), b1 = QgsVector3D::dotProduct( u1, u2 ), c1 = QgsVector3D::dotProduct( u1, d );
1733 float a2 = QgsVector3D::dotProduct( u1, u2 ), b2 = QgsVector3D::dotProduct( u2, u2 ), c2 = QgsVector3D::dotProduct( u2, d );
1734 if ( !( std::fabs( b1 ) > epsilon ) )
1735 {
1736 // Denominator is close to zero.
1737 return false;
1738 }
1739 if ( !( a2 != -1 && a2 != 1 ) )
1740 {
1741 // Lines are parallel
1742 return false;
1743 }
1744
1745 const double r1 = ( c2 - b2 * c1 / b1 ) / ( a2 - b2 * a1 / b1 );
1746 X1 = P1 + u1 * r1;
1747
1748 return true;
1749}
1750
1752 const QgsVector3D &Lb1, const QgsVector3D &Lb2,
1753 QgsVector3D &intersection )
1754{
1755
1756 // if all Vector are on the same plane (have the same Z), use the 2D intersection
1757 // else return a false result
1758 if ( qgsDoubleNear( La1.z(), La2.z() ) && qgsDoubleNear( La1.z(), Lb1.z() ) && qgsDoubleNear( La1.z(), Lb2.z() ) )
1759 {
1760 QgsPoint ptInter;
1761 bool isIntersection;
1762 segmentIntersection( QgsPoint( La1.x(), La1.y() ),
1763 QgsPoint( La2.x(), La2.y() ),
1764 QgsPoint( Lb1.x(), Lb1.y() ),
1765 QgsPoint( Lb2.x(), Lb2.y() ),
1766 ptInter,
1767 isIntersection,
1768 1e-8,
1769 true );
1770 intersection.set( ptInter.x(), ptInter.y(), La1.z() );
1771 return true;
1772 }
1773
1774 // first check if lines have an exact intersection point
1775 // do it by checking if the shortest distance is exactly 0
1776 const float distance = skewLinesDistance( La1, La2, Lb1, Lb2 );
1777 if ( qgsDoubleNear( distance, 0.0 ) )
1778 {
1779 // 3d lines have exact intersection point.
1780 const QgsVector3D C = La2;
1781 const QgsVector3D D = Lb2;
1782 const QgsVector3D e = La1 - La2;
1783 const QgsVector3D f = Lb1 - Lb2;
1784 const QgsVector3D g = D - C;
1785 if ( qgsDoubleNear( ( QgsVector3D::crossProduct( f, g ) ).length(), 0.0 ) || qgsDoubleNear( ( QgsVector3D::crossProduct( f, e ) ).length(), 0.0 ) )
1786 {
1787 // Lines have no intersection, are they parallel?
1788 return false;
1789 }
1790
1792 fgn.normalize();
1793
1795 fen.normalize();
1796
1797 int di = -1;
1798 if ( fgn == fen ) // same direction?
1799 di *= -1;
1800
1801 intersection = C + e * di * ( QgsVector3D::crossProduct( f, g ).length() / QgsVector3D::crossProduct( f, e ).length() );
1802 return true;
1803 }
1804
1805 // try to calculate the approximate intersection point
1806 QgsVector3D X1, X2;
1807 const bool firstIsDone = skewLinesProjection( La1, La2, Lb1, Lb2, X1 );
1808 const bool secondIsDone = skewLinesProjection( Lb1, Lb2, La1, La2, X2 );
1809
1810 if ( !firstIsDone || !secondIsDone )
1811 {
1812 // Could not obtain projection point.
1813 return false;
1814 }
1815
1816 intersection = ( X1 + X2 ) / 2.0;
1817 return true;
1818}
1819
1820double QgsGeometryUtils::triangleArea( double aX, double aY, double bX, double bY, double cX, double cY )
1821{
1822 return 0.5 * std::abs( ( aX - cX ) * ( bY - aY ) - ( aX - bX ) * ( cY - aY ) );
1823}
1824
1825void QgsGeometryUtils::weightedPointInTriangle( const double aX, const double aY, const double bX, const double bY, const double cX, const double cY,
1826 double weightB, double weightC, double &pointX, double &pointY )
1827{
1828 // if point will be outside of the triangle, invert weights
1829 if ( weightB + weightC > 1 )
1830 {
1831 weightB = 1 - weightB;
1832 weightC = 1 - weightC;
1833 }
1834
1835 const double rBx = weightB * ( bX - aX );
1836 const double rBy = weightB * ( bY - aY );
1837 const double rCx = weightC * ( cX - aX );
1838 const double rCy = weightC * ( cY - aY );
1839
1840 pointX = rBx + rCx + aX;
1841 pointY = rBy + rCy + aY;
1842}
1843
1845{
1846 bool rc = false;
1847
1848 for ( const QgsPoint &pt : points )
1849 {
1850 if ( pt.isMeasure() )
1851 {
1852 point.convertTo( QgsWkbTypes::addM( point.wkbType() ) );
1853 point.setM( pt.m() );
1854 rc = true;
1855 break;
1856 }
1857 }
1858
1859 return rc;
1860}
1861
1863{
1864 bool rc = false;
1865
1866 for ( const QgsPoint &pt : points )
1867 {
1868 if ( pt.is3D() )
1869 {
1870 point.convertTo( QgsWkbTypes::addZ( point.wkbType() ) );
1871 point.setZ( pt.z() );
1872 rc = true;
1873 break;
1874 }
1875 }
1876
1877 return rc;
1878}
1879
1880bool QgsGeometryUtils::angleBisector( double aX, double aY, double bX, double bY, double cX, double cY, double dX, double dY,
1881 double &pointX SIP_OUT, double &pointY SIP_OUT, double &angle SIP_OUT )
1882{
1883 const QgsPoint pA = QgsPoint( aX, aY );
1884 const QgsPoint pB = QgsPoint( bX, bY );
1885 const QgsPoint pC = QgsPoint( cX, cY );
1886 const QgsPoint pD = QgsPoint( dX, dY );
1887 angle = ( pA.azimuth( pB ) + pC.azimuth( pD ) ) / 2.0;
1888
1889 QgsPoint pOut;
1890 bool intersection = false;
1891 QgsGeometryUtils::segmentIntersection( pA, pB, pC, pD, pOut, intersection );
1892
1893 pointX = pOut.x();
1894 pointY = pOut.y();
1895
1896 return intersection;
1897}
1898
1899bool QgsGeometryUtils::bisector( double aX, double aY, double bX, double bY, double cX, double cY,
1900 double &pointX SIP_OUT, double &pointY SIP_OUT )
1901{
1902 const QgsPoint pA = QgsPoint( aX, aY );
1903 const QgsPoint pB = QgsPoint( bX, bY );
1904 const QgsPoint pC = QgsPoint( cX, cY );
1905 const double angle = ( pA.azimuth( pB ) + pA.azimuth( pC ) ) / 2.0;
1906
1907 QgsPoint pOut;
1908 bool intersection = false;
1909 QgsGeometryUtils::segmentIntersection( pB, pC, pA, pA.project( 1.0, angle ), pOut, intersection );
1910
1911 pointX = pOut.x();
1912 pointY = pOut.y();
1913
1914 return intersection;
1915}
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.
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 void pointsToWKB(QgsWkbPtr &wkb, const QgsPointSequence &points, bool is3D, bool isMeasure)
Returns a LinearRing { uint32 numPoints; Point points[numPoints]; }.
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 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:2466
double qgsRound(double number, int places)
Returns a double number, rounded (as close as possible) to the specified number of places.
Definition: qgis.h:2581
bool qgsDoubleNear(double a, double b, double epsilon=4 *std::numeric_limits< double >::epsilon())
Compare two doubles (but allow some difference)
Definition: qgis.h:2527
const double DEFAULT_SEGMENT_EPSILON
Default snapping tolerance for segments.
Definition: qgis.h:3020
#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