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