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