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