QGIS API Documentation  3.17.0-Master (a035f434f4)
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  QgsPoint pointAfter = geometry.vertexAt( vertexAfter );
109  if ( vertexAfter.vertex > 0 )
110  {
111  QgsVertexId vertexBefore = vertexAfter;
112  vertexBefore.vertex--;
113  QgsPoint pointBefore = geometry.vertexAt( vertexBefore );
114  double length = pointBefore.distance( pointAfter );
115  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  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  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  double d = v1.y() * v2.x() - v1.x() * v2.y();
245 
246  if ( qgsDoubleNear( d, 0 ) )
247  return false;
248 
249  double dx = p2.x() - p1.x();
250  double dy = p2.y() - p1.y();
251  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 support for intersection point
256  QgsGeometryUtils::setZValueFromPoints( QgsPointSequence() << p1 << p2, intersection );
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  double vl = v.length();
268  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  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  double lambdaw = QgsVector( intersectionPoint.x() - q1.x(), intersectionPoint.y() - q1.y() ) * w;
317  return !( lambdaw < 0. + tolerance || lambdaw >= wl - tolerance );
318 
319 }
320 
321 bool QgsGeometryUtils::lineCircleIntersection( const QgsPointXY &center, const double radius,
322  const QgsPointXY &linePoint1, const QgsPointXY &linePoint2,
323  QgsPointXY &intersection )
324 {
325  // formula taken from http://mathworld.wolfram.com/Circle-LineIntersection.html
326 
327  const double x1 = linePoint1.x() - center.x();
328  const double y1 = linePoint1.y() - center.y();
329  const double x2 = linePoint2.x() - center.x();
330  const double y2 = linePoint2.y() - center.y();
331  const double dx = x2 - x1;
332  const double dy = y2 - y1;
333 
334  const double dr2 = std::pow( dx, 2 ) + std::pow( dy, 2 );
335  const double d = x1 * y2 - x2 * y1;
336 
337  const double disc = std::pow( radius, 2 ) * dr2 - std::pow( d, 2 );
338 
339  if ( disc < 0 )
340  {
341  //no intersection or tangent
342  return false;
343  }
344  else
345  {
346  // two solutions
347  const int sgnDy = dy < 0 ? -1 : 1;
348 
349  const double sqrDisc = std::sqrt( disc );
350 
351  const double ax = center.x() + ( d * dy + sgnDy * dx * sqrDisc ) / dr2;
352  const double ay = center.y() + ( -d * dx + std::fabs( dy ) * sqrDisc ) / dr2;
353  const QgsPointXY p1( ax, ay );
354 
355  const double bx = center.x() + ( d * dy - sgnDy * dx * sqrDisc ) / dr2;
356  const double by = center.y() + ( -d * dx - std::fabs( dy ) * sqrDisc ) / dr2;
357  const QgsPointXY p2( bx, by );
358 
359  // snap to nearest intersection
360 
361  if ( intersection.sqrDist( p1 ) < intersection.sqrDist( p2 ) )
362  {
363  intersection.set( p1.x(), p1.y() );
364  }
365  else
366  {
367  intersection.set( p2.x(), p2.y() );
368  }
369  return true;
370  }
371 }
372 
373 // based on public domain work by 3/26/2005 Tim Voght
374 // see http://paulbourke.net/geometry/circlesphere/tvoght.c
375 int QgsGeometryUtils::circleCircleIntersections( QgsPointXY center1, const double r1, QgsPointXY center2, const double r2, QgsPointXY &intersection1, QgsPointXY &intersection2 )
376 {
377  // determine the straight-line distance between the centers
378  const double d = center1.distance( center2 );
379 
380  // check for solvability
381  if ( d > ( r1 + r2 ) )
382  {
383  // no solution. circles do not intersect.
384  return 0;
385  }
386  else if ( d < std::fabs( r1 - r2 ) )
387  {
388  // no solution. one circle is contained in the other
389  return 0;
390  }
391  else if ( qgsDoubleNear( d, 0 ) && ( qgsDoubleNear( r1, r2 ) ) )
392  {
393  // no solutions, the circles coincide
394  return 0;
395  }
396 
397  /* 'point 2' is the point where the line through the circle
398  * intersection points crosses the line between the circle
399  * centers.
400  */
401 
402  // Determine the distance from point 0 to point 2.
403  const double a = ( ( r1 * r1 ) - ( r2 * r2 ) + ( d * d ) ) / ( 2.0 * d ) ;
404 
405  /* dx and dy are the vertical and horizontal distances between
406  * the circle centers.
407  */
408  const double dx = center2.x() - center1.x();
409  const double dy = center2.y() - center1.y();
410 
411  // Determine the coordinates of point 2.
412  const double x2 = center1.x() + ( dx * a / d );
413  const double y2 = center1.y() + ( dy * a / d );
414 
415  /* Determine the distance from point 2 to either of the
416  * intersection points.
417  */
418  const double h = std::sqrt( ( r1 * r1 ) - ( a * a ) );
419 
420  /* Now determine the offsets of the intersection points from
421  * point 2.
422  */
423  const double rx = dy * ( h / d );
424  const double ry = dx * ( h / d );
425 
426  // determine the absolute intersection points
427  intersection1 = QgsPointXY( x2 + rx, y2 - ry );
428  intersection2 = QgsPointXY( x2 - rx, y2 + ry );
429 
430  // see if we have 1 or 2 solutions
431  if ( qgsDoubleNear( d, r1 + r2 ) )
432  return 1;
433 
434  return 2;
435 }
436 
437 // Using https://stackoverflow.com/a/1351794/1861260
438 // and inspired by http://csharphelper.com/blog/2014/11/find-the-tangent-lines-between-a-point-and-a-circle-in-c/
439 bool QgsGeometryUtils::tangentPointAndCircle( const QgsPointXY &center, double radius, const QgsPointXY &p, QgsPointXY &pt1, QgsPointXY &pt2 )
440 {
441  // distance from point to center of circle
442  const double dx = center.x() - p.x();
443  const double dy = center.y() - p.y();
444  const double distanceSquared = dx * dx + dy * dy;
445  const double radiusSquared = radius * radius;
446  if ( distanceSquared < radiusSquared )
447  {
448  // point is inside circle!
449  return false;
450  }
451 
452  // distance from point to tangent point, using pythagoras
453  const double distanceToTangent = std::sqrt( distanceSquared - radiusSquared );
454 
455  // tangent points are those where the original circle intersects a circle centered
456  // on p with radius distanceToTangent
457  circleCircleIntersections( center, radius, p, distanceToTangent, pt1, pt2 );
458 
459  return true;
460 }
461 
462 // inspired by http://csharphelper.com/blog/2014/12/find-the-tangent-lines-between-two-circles-in-c/
463 int QgsGeometryUtils::circleCircleOuterTangents( const QgsPointXY &center1, double radius1, const QgsPointXY &center2, double radius2, QgsPointXY &line1P1, QgsPointXY &line1P2, QgsPointXY &line2P1, QgsPointXY &line2P2 )
464 {
465  if ( radius1 > radius2 )
466  return circleCircleOuterTangents( center2, radius2, center1, radius1, line1P1, line1P2, line2P1, line2P2 );
467 
468  const double radius2a = radius2 - radius1;
469  if ( !tangentPointAndCircle( center2, radius2a, center1, line1P2, line2P2 ) )
470  {
471  // there are no tangents
472  return 0;
473  }
474 
475  // get the vector perpendicular to the
476  // first tangent with length radius1
477  QgsVector v1( -( line1P2.y() - center1.y() ), line1P2.x() - center1.x() );
478  const double v1Length = v1.length();
479  v1 = v1 * ( radius1 / v1Length );
480 
481  // offset the tangent vector's points
482  line1P1 = center1 + v1;
483  line1P2 = line1P2 + v1;
484 
485  // get the vector perpendicular to the
486  // second tangent with length radius1
487  QgsVector v2( line2P2.y() - center1.y(), -( line2P2.x() - center1.x() ) );
488  const double v2Length = v2.length();
489  v2 = v2 * ( radius1 / v2Length );
490 
491  // offset the tangent vector's points
492  line2P1 = center1 + v2;
493  line2P2 = line2P2 + v2;
494 
495  return 2;
496 }
497 
498 // inspired by http://csharphelper.com/blog/2014/12/find-the-tangent-lines-between-two-circles-in-c/
499 int QgsGeometryUtils::circleCircleInnerTangents( const QgsPointXY &center1, double radius1, const QgsPointXY &center2, double radius2, QgsPointXY &line1P1, QgsPointXY &line1P2, QgsPointXY &line2P1, QgsPointXY &line2P2 )
500 {
501  if ( radius1 > radius2 )
502  return circleCircleInnerTangents( center2, radius2, center1, radius1, line1P1, line1P2, line2P1, line2P2 );
503 
504  // determine the straight-line distance between the centers
505  const double d = center1.distance( center2 );
506  const double radius1a = radius1 + radius2;
507 
508  // check for solvability
509  if ( d <= radius1a || qgsDoubleNear( d, radius1a ) )
510  {
511  // no solution. circles intersect or touch.
512  return 0;
513  }
514 
515  if ( !tangentPointAndCircle( center1, radius1a, center2, line1P2, line2P2 ) )
516  {
517  // there are no tangents
518  return 0;
519  }
520 
521  // get the vector perpendicular to the
522  // first tangent with length radius2
523  QgsVector v1( ( line1P2.y() - center2.y() ), -( line1P2.x() - center2.x() ) );
524  const double v1Length = v1.length();
525  v1 = v1 * ( radius2 / v1Length );
526 
527  // offset the tangent vector's points
528  line1P1 = center2 + v1;
529  line1P2 = line1P2 + v1;
530 
531  // get the vector perpendicular to the
532  // second tangent with length radius2
533  QgsVector v2( -( line2P2.y() - center2.y() ), line2P2.x() - center2.x() );
534  const double v2Length = v2.length();
535  v2 = v2 * ( radius2 / v2Length );
536 
537  // offset the tangent vector's points in opposite direction
538  line2P1 = center2 + v2;
539  line2P2 = line2P2 + v2;
540 
541  return 2;
542 }
543 
544 QVector<QgsGeometryUtils::SelfIntersection> QgsGeometryUtils::selfIntersections( const QgsAbstractGeometry *geom, int part, int ring, double tolerance )
545 {
546  QVector<SelfIntersection> intersections;
547 
548  int n = geom->vertexCount( part, ring );
549  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  QgsPoint pi = geom->vertexAt( QgsVertexId( part, ring, i ) );
555  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  int start = j + 1;
560  int end = i == 0 && isClosed ? n - 1 : n;
561  for ( int k = start, l = start + 1; l < end; k = l++ )
562  {
563  QgsPoint pk = geom->vertexAt( QgsVertexId( part, ring, k ) );
564  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  double f1 = x - x1;
592  double f2 = y2 - y1;
593  double f3 = y - y1;
594  double f4 = x2 - x1;
595  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 QgsPoint QgsGeometryUtils::interpolatePointOnArc( const QgsPoint &pt1, const QgsPoint &pt2, const QgsPoint &pt3, double distance )
635 {
636  double centerX, centerY, radius;
637  circleCenterRadius( pt1, pt2, pt3, radius, centerX, centerY );
638 
639  const double theta = distance / radius; // angle subtended
640  const double anglePt1 = std::atan2( pt1.y() - centerY, pt1.x() - centerX );
641  const double anglePt2 = std::atan2( pt2.y() - centerY, pt2.x() - centerX );
642  const double anglePt3 = std::atan2( pt3.y() - centerY, pt3.x() - centerX );
643  const bool isClockwise = circleClockwise( anglePt1, anglePt2, anglePt3 );
644  const double angleDest = anglePt1 + ( isClockwise ? -theta : theta );
645 
646  const double x = centerX + radius * ( std::cos( angleDest ) );
647  const double y = centerY + radius * ( std::sin( angleDest ) );
648 
649  const double z = pt1.is3D() ?
650  interpolateArcValue( angleDest, anglePt1, anglePt2, anglePt3, pt1.z(), pt2.z(), pt3.z() )
651  : 0;
652  const double m = pt1.isMeasure() ?
653  interpolateArcValue( angleDest, anglePt1, anglePt2, anglePt3, pt1.m(), pt2.m(), pt3.m() )
654  : 0;
655 
656  return QgsPoint( pt1.wkbType(), x, y, z, m );
657 }
658 
659 double QgsGeometryUtils::ccwAngle( double dy, double dx )
660 {
661  double angle = std::atan2( dy, dx ) * 180 / M_PI;
662  if ( angle < 0 )
663  {
664  return 360 + angle;
665  }
666  else if ( angle > 360 )
667  {
668  return 360 - angle;
669  }
670  return angle;
671 }
672 
673 void QgsGeometryUtils::circleCenterRadius( const QgsPoint &pt1, const QgsPoint &pt2, const QgsPoint &pt3, double &radius, double &centerX, double &centerY )
674 {
675  double dx21, dy21, dx31, dy31, h21, h31, d;
676 
677  //closed circle
678  if ( qgsDoubleNear( pt1.x(), pt3.x() ) && qgsDoubleNear( pt1.y(), pt3.y() ) )
679  {
680  centerX = ( pt1.x() + pt2.x() ) / 2.0;
681  centerY = ( pt1.y() + pt2.y() ) / 2.0;
682  radius = std::sqrt( std::pow( centerX - pt1.x(), 2.0 ) + std::pow( centerY - pt1.y(), 2.0 ) );
683  return;
684  }
685 
686  // Using Cartesian circumcenter eguations from page https://en.wikipedia.org/wiki/Circumscribed_circle
687  dx21 = pt2.x() - pt1.x();
688  dy21 = pt2.y() - pt1.y();
689  dx31 = pt3.x() - pt1.x();
690  dy31 = pt3.y() - pt1.y();
691 
692  h21 = std::pow( dx21, 2.0 ) + std::pow( dy21, 2.0 );
693  h31 = std::pow( dx31, 2.0 ) + std::pow( dy31, 2.0 );
694 
695  // 2*Cross product, d<0 means clockwise and d>0 counterclockwise sweeping angle
696  d = 2 * ( dx21 * dy31 - dx31 * dy21 );
697 
698  // Check colinearity, Cross product = 0
699  if ( qgsDoubleNear( std::fabs( d ), 0.0, 0.00000000001 ) )
700  {
701  radius = -1.0;
702  return;
703  }
704 
705  // Calculate centroid coordinates and radius
706  centerX = pt1.x() + ( h21 * dy31 - h31 * dy21 ) / d;
707  centerY = pt1.y() - ( h21 * dx31 - h31 * dx21 ) / d;
708  radius = std::sqrt( std::pow( centerX - pt1.x(), 2.0 ) + std::pow( centerY - pt1.y(), 2.0 ) );
709 }
710 
711 bool QgsGeometryUtils::circleClockwise( double angle1, double angle2, double angle3 )
712 {
713  if ( angle3 >= angle1 )
714  {
715  return !( angle2 > angle1 && angle2 < angle3 );
716  }
717  else
718  {
719  return !( angle2 > angle1 || angle2 < angle3 );
720  }
721 }
722 
723 bool QgsGeometryUtils::circleAngleBetween( double angle, double angle1, double angle2, bool clockwise )
724 {
725  if ( clockwise )
726  {
727  if ( angle2 < angle1 )
728  {
729  return ( angle <= angle1 && angle >= angle2 );
730  }
731  else
732  {
733  return ( angle <= angle1 || angle >= angle2 );
734  }
735  }
736  else
737  {
738  if ( angle2 > angle1 )
739  {
740  return ( angle >= angle1 && angle <= angle2 );
741  }
742  else
743  {
744  return ( angle >= angle1 || angle <= angle2 );
745  }
746  }
747 }
748 
749 bool QgsGeometryUtils::angleOnCircle( double angle, double angle1, double angle2, double angle3 )
750 {
751  bool clockwise = circleClockwise( angle1, angle2, angle3 );
752  return circleAngleBetween( angle, angle1, angle3, clockwise );
753 }
754 
755 double QgsGeometryUtils::circleLength( double x1, double y1, double x2, double y2, double x3, double y3 )
756 {
757  double centerX, centerY, radius;
758  circleCenterRadius( QgsPoint( x1, y1 ), QgsPoint( x2, y2 ), QgsPoint( x3, y3 ), radius, centerX, centerY );
759  double length = M_PI / 180.0 * radius * sweepAngle( centerX, centerY, x1, y1, x2, y2, x3, y3 );
760  if ( length < 0 )
761  {
762  length = -length;
763  }
764  return length;
765 }
766 
767 double QgsGeometryUtils::sweepAngle( double centerX, double centerY, double x1, double y1, double x2, double y2, double x3, double y3 )
768 {
769  double p1Angle = QgsGeometryUtils::ccwAngle( y1 - centerY, x1 - centerX );
770  double p2Angle = QgsGeometryUtils::ccwAngle( y2 - centerY, x2 - centerX );
771  double p3Angle = QgsGeometryUtils::ccwAngle( y3 - centerY, x3 - centerX );
772 
773  if ( p3Angle >= p1Angle )
774  {
775  if ( p2Angle > p1Angle && p2Angle < p3Angle )
776  {
777  return ( p3Angle - p1Angle );
778  }
779  else
780  {
781  return ( - ( p1Angle + ( 360 - p3Angle ) ) );
782  }
783  }
784  else
785  {
786  if ( p2Angle < p1Angle && p2Angle > p3Angle )
787  {
788  return ( -( p1Angle - p3Angle ) );
789  }
790  else
791  {
792  return ( p3Angle + ( 360 - p1Angle ) );
793  }
794  }
795 }
796 
797 bool QgsGeometryUtils::segmentMidPoint( const QgsPoint &p1, const QgsPoint &p2, QgsPoint &result, double radius, const QgsPoint &mousePos )
798 {
799  QgsPoint midPoint( ( p1.x() + p2.x() ) / 2.0, ( p1.y() + p2.y() ) / 2.0 );
800  double midDist = std::sqrt( sqrDistance2D( p1, midPoint ) );
801  if ( radius < midDist )
802  {
803  return false;
804  }
805  double centerMidDist = std::sqrt( radius * radius - midDist * midDist );
806  double dist = radius - centerMidDist;
807 
808  double midDx = midPoint.x() - p1.x();
809  double midDy = midPoint.y() - p1.y();
810 
811  //get the four possible midpoints
812  QVector<QgsPoint> possibleMidPoints;
813  possibleMidPoints.append( pointOnLineWithDistance( midPoint, QgsPoint( midPoint.x() - midDy, midPoint.y() + midDx ), dist ) );
814  possibleMidPoints.append( pointOnLineWithDistance( midPoint, QgsPoint( midPoint.x() - midDy, midPoint.y() + midDx ), 2 * radius - dist ) );
815  possibleMidPoints.append( pointOnLineWithDistance( midPoint, QgsPoint( midPoint.x() + midDy, midPoint.y() - midDx ), dist ) );
816  possibleMidPoints.append( pointOnLineWithDistance( midPoint, QgsPoint( midPoint.x() + midDy, midPoint.y() - midDx ), 2 * radius - dist ) );
817 
818  //take the closest one
819  double minDist = std::numeric_limits<double>::max();
820  int minDistIndex = -1;
821  for ( int i = 0; i < possibleMidPoints.size(); ++i )
822  {
823  double currentDist = sqrDistance2D( mousePos, possibleMidPoints.at( i ) );
824  if ( currentDist < minDist )
825  {
826  minDistIndex = i;
827  minDist = currentDist;
828  }
829  }
830 
831  if ( minDistIndex == -1 )
832  {
833  return false;
834  }
835 
836  result = possibleMidPoints.at( minDistIndex );
837 
838  // add z support if necessary
840 
841  return true;
842 }
843 
844 QgsPoint QgsGeometryUtils::segmentMidPointFromCenter( const QgsPoint &p1, const QgsPoint &p2, const QgsPoint &center, const bool useShortestArc )
845 {
846  double midPointAngle = averageAngle( lineAngle( center.x(), center.y(), p1.x(), p1.y() ),
847  lineAngle( center.x(), center.y(), p2.x(), p2.y() ) );
848  if ( !useShortestArc )
849  midPointAngle += M_PI;
850  return center.project( center.distance( p1 ), midPointAngle * 180 / M_PI );
851 }
852 
853 double QgsGeometryUtils::circleTangentDirection( const QgsPoint &tangentPoint, const QgsPoint &cp1,
854  const QgsPoint &cp2, const QgsPoint &cp3 )
855 {
856  //calculate circle midpoint
857  double mX, mY, radius;
858  circleCenterRadius( cp1, cp2, cp3, radius, mX, mY );
859 
860  double p1Angle = QgsGeometryUtils::ccwAngle( cp1.y() - mY, cp1.x() - mX );
861  double p2Angle = QgsGeometryUtils::ccwAngle( cp2.y() - mY, cp2.x() - mX );
862  double p3Angle = QgsGeometryUtils::ccwAngle( cp3.y() - mY, cp3.x() - mX );
863  double angle = 0;
864  if ( circleClockwise( p1Angle, p2Angle, p3Angle ) )
865  {
866  angle = lineAngle( tangentPoint.x(), tangentPoint.y(), mX, mY ) - M_PI_2;
867  }
868  else
869  {
870  angle = lineAngle( mX, mY, tangentPoint.x(), tangentPoint.y() ) - M_PI_2;
871  }
872  if ( angle < 0 )
873  angle += 2 * M_PI;
874  return angle;
875 }
876 
877 // Ported from PostGIS' pt_continues_arc
878 bool QgsGeometryUtils::pointContinuesArc( const QgsPoint &a1, const QgsPoint &a2, const QgsPoint &a3, const QgsPoint &b, double distanceTolerance, double pointSpacingAngleTolerance )
879 {
880  double centerX = 0;
881  double centerY = 0;
882  double radius = 0;
883  circleCenterRadius( a1, a2, a3, radius, centerX, centerY );
884 
885  // Co-linear a1/a2/a3
886  if ( radius < 0.0 )
887  return false;
888 
889  // distance of candidate point to center of arc a1-a2-a3
890  double bDistance = std::sqrt( ( b.x() - centerX ) * ( b.x() - centerX ) +
891  ( b.y() - centerY ) * ( b.y() - centerY ) );
892 
893  double diff = std::fabs( radius - bDistance );
894 
895  auto arcAngle = []( const QgsPoint & a, const QgsPoint & b, const QgsPoint & c )->double
896  {
897  double abX = b.x() - a.x();
898  double abY = b.y() - a.y();
899 
900  double cbX = b.x() - c.x();
901  double cbY = b.y() - c.y();
902 
903  double dot = ( abX * cbX + abY * cbY ); /* dot product */
904  double cross = ( abX * cbY - abY * cbX ); /* cross product */
905 
906  double alpha = std::atan2( cross, dot );
907 
908  return alpha;
909  };
910 
911  // Is the point b on the circle?
912  if ( diff < distanceTolerance )
913  {
914  double angle1 = arcAngle( a1, a2, a3 );
915  double angle2 = arcAngle( a2, a3, b );
916 
917  // Is the sweep angle similar to the previous one?
918  // We only consider a segment replaceable by an arc if the points within
919  // it are regularly spaced
920  diff = std::fabs( angle1 - angle2 );
921  if ( diff > pointSpacingAngleTolerance )
922  {
923  return false;
924  }
925 
926  int a2Side = leftOfLine( a2.x(), a2.y(), a1.x(), a1.y(), a3.x(), a3.y() );
927  int bSide = leftOfLine( b.x(), b.y(), a1.x(), a1.y(), a3.x(), a3.y() );
928 
929  // Is the point b on the same side of a1/a3 as the mid-point a2 is?
930  // If not, it's in the unbounded part of the circle, so it continues the arc, return true.
931  if ( bSide != a2Side )
932  return true;
933  }
934  return false;
935 }
936 
937 void QgsGeometryUtils::segmentizeArc( const QgsPoint &p1, const QgsPoint &p2, const QgsPoint &p3, QgsPointSequence &points, double tolerance, QgsAbstractGeometry::SegmentationToleranceType toleranceType, bool hasZ, bool hasM )
938 {
939  bool reversed = false;
940  int segSide = segmentSide( p1, p3, p2 );
941 
942  QgsPoint circlePoint1;
943  const QgsPoint circlePoint2 = p2;
944  QgsPoint circlePoint3;
945 
946  if ( segSide == -1 )
947  {
948  // Reverse !
949  circlePoint1 = p3;
950  circlePoint3 = p1;
951  reversed = true;
952  }
953  else
954  {
955  circlePoint1 = p1;
956  circlePoint3 = p3;
957  }
958 
959  //adapted code from PostGIS
960  double radius = 0;
961  double centerX = 0;
962  double centerY = 0;
963  circleCenterRadius( circlePoint1, circlePoint2, circlePoint3, radius, centerX, centerY );
964 
965  if ( circlePoint1 != circlePoint3 && ( radius < 0 || qgsDoubleNear( segSide, 0.0 ) ) ) //points are colinear
966  {
967  points.append( p1 );
968  points.append( p2 );
969  points.append( p3 );
970  return;
971  }
972 
973  double increment = tolerance; //one segment per degree
974  if ( toleranceType == QgsAbstractGeometry::MaximumDifference )
975  {
976  // Ensure tolerance is not higher than twice the radius
977  tolerance = std::min( tolerance, radius * 2 );
978  double halfAngle = std::acos( -tolerance / radius + 1 );
979  increment = 2 * halfAngle;
980  }
981 
982  //angles of pt1, pt2, pt3
983  double a1 = std::atan2( circlePoint1.y() - centerY, circlePoint1.x() - centerX );
984  double a2 = std::atan2( circlePoint2.y() - centerY, circlePoint2.x() - centerX );
985  double a3 = std::atan2( circlePoint3.y() - centerY, circlePoint3.x() - centerX );
986 
987  // Make segmentation symmetric
988  const bool symmetric = true;
989  if ( symmetric )
990  {
991  double angle = a3 - a1;
992  // angle == 0 when full circle
993  if ( angle <= 0 ) angle += M_PI * 2;
994 
995  /* Number of segments in output */
996  int segs = ceil( angle / increment );
997  /* Tweak increment to be regular for all the arc */
998  increment = angle / segs;
999  }
1000 
1001  /* Adjust a3 up so we can increment from a1 to a3 cleanly */
1002  // a3 == a1 when full circle
1003  if ( a3 <= a1 )
1004  a3 += 2.0 * M_PI;
1005  if ( a2 < a1 )
1006  a2 += 2.0 * M_PI;
1007 
1008  double x, y;
1009  double z = 0;
1010  double m = 0;
1011 
1012  QVector<QgsPoint> stringPoints;
1013  stringPoints.insert( 0, circlePoint1 );
1014  if ( circlePoint2 != circlePoint3 && circlePoint1 != circlePoint2 ) //draw straight line segment if two points have the same position
1015  {
1016  QgsWkbTypes::Type pointWkbType = QgsWkbTypes::Point;
1017  if ( hasZ )
1018  pointWkbType = QgsWkbTypes::addZ( pointWkbType );
1019  if ( hasM )
1020  pointWkbType = QgsWkbTypes::addM( pointWkbType );
1021 
1022  // As we're adding the last point in any case, we'll avoid
1023  // including a point which is at less than 1% increment distance
1024  // from it (may happen to find them due to numbers approximation).
1025  // NOTE that this effectively allows in output some segments which
1026  // are more distant than requested. This is at most 1% off
1027  // from requested MaxAngle and less for MaxError.
1028  double tolError = increment / 100;
1029  double stopAngle = a3 - tolError;
1030  for ( double angle = a1 + increment; angle < stopAngle; angle += increment )
1031  {
1032  x = centerX + radius * std::cos( angle );
1033  y = centerY + radius * std::sin( angle );
1034 
1035  if ( hasZ )
1036  {
1037  z = interpolateArcValue( angle, a1, a2, a3, circlePoint1.z(), circlePoint2.z(), circlePoint3.z() );
1038  }
1039  if ( hasM )
1040  {
1041  m = interpolateArcValue( angle, a1, a2, a3, circlePoint1.m(), circlePoint2.m(), circlePoint3.m() );
1042  }
1043 
1044  stringPoints.insert( stringPoints.size(), QgsPoint( pointWkbType, x, y, z, m ) );
1045  }
1046  }
1047  stringPoints.insert( stringPoints.size(), circlePoint3 );
1048 
1049  // TODO: check if or implement QgsPointSequence directly taking an iterator to append
1050  if ( reversed )
1051  {
1052  std::reverse( stringPoints.begin(), stringPoints.end() );
1053  }
1054  if ( ! points.empty() && stringPoints.front() == points.back() ) stringPoints.pop_front();
1055  points.append( stringPoints );
1056 }
1057 
1058 int QgsGeometryUtils::segmentSide( const QgsPoint &pt1, const QgsPoint &pt3, const QgsPoint &pt2 )
1059 {
1060  double side = ( ( pt2.x() - pt1.x() ) * ( pt3.y() - pt1.y() ) - ( pt3.x() - pt1.x() ) * ( pt2.y() - pt1.y() ) );
1061  if ( side == 0.0 )
1062  {
1063  return 0;
1064  }
1065  else
1066  {
1067  if ( side < 0 )
1068  {
1069  return -1;
1070  }
1071  if ( side > 0 )
1072  {
1073  return 1;
1074  }
1075  return 0;
1076  }
1077 }
1078 
1079 double QgsGeometryUtils::interpolateArcValue( double angle, double a1, double a2, double a3, double zm1, double zm2, double zm3 )
1080 {
1081  /* Counter-clockwise sweep */
1082  if ( a1 < a2 )
1083  {
1084  if ( angle <= a2 )
1085  return zm1 + ( zm2 - zm1 ) * ( angle - a1 ) / ( a2 - a1 );
1086  else
1087  return zm2 + ( zm3 - zm2 ) * ( angle - a2 ) / ( a3 - a2 );
1088  }
1089  /* Clockwise sweep */
1090  else
1091  {
1092  if ( angle >= a2 )
1093  return zm1 + ( zm2 - zm1 ) * ( a1 - angle ) / ( a1 - a2 );
1094  else
1095  return zm2 + ( zm3 - zm2 ) * ( a2 - angle ) / ( a2 - a3 );
1096  }
1097 }
1098 
1099 QgsPointSequence QgsGeometryUtils::pointsFromWKT( const QString &wktCoordinateList, bool is3D, bool isMeasure )
1100 {
1101  int dim = 2 + is3D + isMeasure;
1102  QgsPointSequence points;
1103 
1104  const QStringList coordList = wktCoordinateList.split( ',', QString::SkipEmptyParts );
1105 
1106  //first scan through for extra unexpected dimensions
1107  bool foundZ = false;
1108  bool foundM = false;
1109  QRegularExpression rx( QStringLiteral( "\\s" ) );
1110  QRegularExpression rxIsNumber( QStringLiteral( "^[+-]?(\\d\\.?\\d*[Ee][+\\-]?\\d+|(\\d+\\.\\d*|\\d*\\.\\d+)|\\d+)$" ) );
1111  for ( const QString &pointCoordinates : coordList )
1112  {
1113  QStringList coordinates = pointCoordinates.split( rx, QString::SkipEmptyParts );
1114 
1115  // exit with an empty set if one list contains invalid value.
1116  if ( coordinates.filter( rxIsNumber ).size() != coordinates.size() )
1117  return points;
1118 
1119  if ( coordinates.size() == 3 && !foundZ && !foundM && !is3D && !isMeasure )
1120  {
1121  // 3 dimensional coordinates, but not specifically marked as such. We allow this
1122  // anyway and upgrade geometry to have Z dimension
1123  foundZ = true;
1124  }
1125  else if ( coordinates.size() >= 4 && ( !( is3D || foundZ ) || !( isMeasure || foundM ) ) )
1126  {
1127  // 4 (or more) dimensional coordinates, but not specifically marked as such. We allow this
1128  // anyway and upgrade geometry to have Z&M dimensions
1129  foundZ = true;
1130  foundM = true;
1131  }
1132  }
1133 
1134  for ( const QString &pointCoordinates : coordList )
1135  {
1136  QStringList coordinates = pointCoordinates.split( rx, QString::SkipEmptyParts );
1137  if ( coordinates.size() < dim )
1138  continue;
1139 
1140  int idx = 0;
1141  double x = coordinates[idx++].toDouble();
1142  double y = coordinates[idx++].toDouble();
1143 
1144  double z = 0;
1145  if ( ( is3D || foundZ ) && coordinates.length() > idx )
1146  z = coordinates[idx++].toDouble();
1147 
1148  double m = 0;
1149  if ( ( isMeasure || foundM ) && coordinates.length() > idx )
1150  m = coordinates[idx++].toDouble();
1151 
1153  if ( is3D || foundZ )
1154  {
1155  if ( isMeasure || foundM )
1157  else
1158  t = QgsWkbTypes::PointZ;
1159  }
1160  else
1161  {
1162  if ( isMeasure || foundM )
1163  t = QgsWkbTypes::PointM;
1164  else
1165  t = QgsWkbTypes::Point;
1166  }
1167 
1168  points.append( QgsPoint( t, x, y, z, m ) );
1169  }
1170 
1171  return points;
1172 }
1173 
1175 {
1176  wkb << static_cast<quint32>( points.size() );
1177  for ( const QgsPoint &point : points )
1178  {
1179  wkb << point.x() << point.y();
1180  if ( is3D )
1181  {
1182  wkb << point.z();
1183  }
1184  if ( isMeasure )
1185  {
1186  wkb << point.m();
1187  }
1188  }
1189 }
1190 
1191 QString QgsGeometryUtils::pointsToWKT( const QgsPointSequence &points, int precision, bool is3D, bool isMeasure )
1192 {
1193  QString wkt = QStringLiteral( "(" );
1194  for ( const QgsPoint &p : points )
1195  {
1196  wkt += qgsDoubleToString( p.x(), precision );
1197  wkt += ' ' + qgsDoubleToString( p.y(), precision );
1198  if ( is3D )
1199  wkt += ' ' + qgsDoubleToString( p.z(), precision );
1200  if ( isMeasure )
1201  wkt += ' ' + qgsDoubleToString( p.m(), precision );
1202  wkt += QLatin1String( ", " );
1203  }
1204  if ( wkt.endsWith( QLatin1String( ", " ) ) )
1205  wkt.chop( 2 ); // Remove last ", "
1206  wkt += ')';
1207  return wkt;
1208 }
1209 
1210 QDomElement QgsGeometryUtils::pointsToGML2( const QgsPointSequence &points, QDomDocument &doc, int precision, const QString &ns, QgsAbstractGeometry::AxisOrder axisOrder )
1211 {
1212  QDomElement elemCoordinates = doc.createElementNS( ns, QStringLiteral( "coordinates" ) );
1213 
1214  // coordinate separator
1215  QString cs = QStringLiteral( "," );
1216  // tuple separator
1217  QString ts = QStringLiteral( " " );
1218 
1219  elemCoordinates.setAttribute( QStringLiteral( "cs" ), cs );
1220  elemCoordinates.setAttribute( QStringLiteral( "ts" ), ts );
1221 
1222  QString strCoordinates;
1223 
1224  for ( const QgsPoint &p : points )
1225  if ( axisOrder == QgsAbstractGeometry::AxisOrder::XY )
1226  strCoordinates += qgsDoubleToString( p.x(), precision ) + cs + qgsDoubleToString( p.y(), precision ) + ts;
1227  else
1228  strCoordinates += qgsDoubleToString( p.y(), precision ) + cs + qgsDoubleToString( p.x(), precision ) + ts;
1229 
1230  if ( strCoordinates.endsWith( ts ) )
1231  strCoordinates.chop( 1 ); // Remove trailing space
1232 
1233  elemCoordinates.appendChild( doc.createTextNode( strCoordinates ) );
1234  return elemCoordinates;
1235 }
1236 
1237 QDomElement QgsGeometryUtils::pointsToGML3( const QgsPointSequence &points, QDomDocument &doc, int precision, const QString &ns, bool is3D, QgsAbstractGeometry::AxisOrder axisOrder )
1238 {
1239  QDomElement elemPosList = doc.createElementNS( ns, QStringLiteral( "posList" ) );
1240  elemPosList.setAttribute( QStringLiteral( "srsDimension" ), is3D ? 3 : 2 );
1241 
1242  QString strCoordinates;
1243  for ( const QgsPoint &p : points )
1244  {
1245  if ( axisOrder == QgsAbstractGeometry::AxisOrder::XY )
1246  strCoordinates += qgsDoubleToString( p.x(), precision ) + ' ' + qgsDoubleToString( p.y(), precision ) + ' ';
1247  else
1248  strCoordinates += qgsDoubleToString( p.y(), precision ) + ' ' + qgsDoubleToString( p.x(), precision ) + ' ';
1249  if ( is3D )
1250  strCoordinates += qgsDoubleToString( p.z(), precision ) + ' ';
1251  }
1252  if ( strCoordinates.endsWith( ' ' ) )
1253  strCoordinates.chop( 1 ); // Remove trailing space
1254 
1255  elemPosList.appendChild( doc.createTextNode( strCoordinates ) );
1256  return elemPosList;
1257 }
1258 
1260 {
1261  QString json = QStringLiteral( "[ " );
1262  for ( const QgsPoint &p : points )
1263  {
1264  json += '[' + qgsDoubleToString( p.x(), precision ) + QLatin1String( ", " ) + qgsDoubleToString( p.y(), precision ) + QLatin1String( "], " );
1265  }
1266  if ( json.endsWith( QLatin1String( ", " ) ) )
1267  {
1268  json.chop( 2 ); // Remove last ", "
1269  }
1270  json += ']';
1271  return json;
1272 }
1273 
1274 
1276 {
1277  json coordinates( json::array() );
1278  for ( const QgsPoint &p : points )
1279  {
1280  if ( p.is3D() )
1281  {
1282  coordinates.push_back( { qgsRound( p.x(), precision ), qgsRound( p.y(), precision ), qgsRound( p.z(), precision ) } );
1283  }
1284  else
1285  {
1286  coordinates.push_back( { qgsRound( p.x(), precision ), qgsRound( p.y(), precision ) } );
1287  }
1288  }
1289  return coordinates;
1290 }
1291 
1293 {
1294  double clippedAngle = angle;
1295  if ( clippedAngle >= M_PI * 2 || clippedAngle <= -2 * M_PI )
1296  {
1297  clippedAngle = std::fmod( clippedAngle, 2 * M_PI );
1298  }
1299  if ( clippedAngle < 0.0 )
1300  {
1301  clippedAngle += 2 * M_PI;
1302  }
1303  return clippedAngle;
1304 }
1305 
1306 QPair<QgsWkbTypes::Type, QString> QgsGeometryUtils::wktReadBlock( const QString &wkt )
1307 {
1308  QString wktParsed = wkt;
1309  QString contents;
1310  if ( wkt.contains( QString( "EMPTY" ), Qt::CaseInsensitive ) )
1311  {
1312  QRegularExpression wktRegEx( QStringLiteral( "^\\s*(\\w+)\\s+(\\w+)\\s*$" ) );
1313  wktRegEx.setPatternOptions( QRegularExpression::DotMatchesEverythingOption );
1314  QRegularExpressionMatch match = wktRegEx.match( wkt );
1315  if ( match.hasMatch() )
1316  {
1317  wktParsed = match.captured( 1 );
1318  contents = match.captured( 2 ).toUpper();
1319  }
1320  }
1321  else
1322  {
1323  const int openedParenthesisCount = wktParsed.count( '(' );
1324  const int closedParenthesisCount = wktParsed.count( ')' );
1325  // closes missing parentheses
1326  for ( int i = 0 ; i < openedParenthesisCount - closedParenthesisCount; ++i )
1327  wktParsed.push_back( ')' );
1328  // removes extra parentheses
1329  wktParsed.truncate( wktParsed.size() - ( closedParenthesisCount - openedParenthesisCount ) );
1330 
1331  QRegularExpression cooRegEx( QStringLiteral( "^[^\\(]*\\((.*)\\)[^\\)]*$" ) );
1332  cooRegEx.setPatternOptions( QRegularExpression::DotMatchesEverythingOption );
1333  QRegularExpressionMatch match = cooRegEx.match( wktParsed );
1334  contents = match.hasMatch() ? match.captured( 1 ) : QString();
1335  }
1337  return qMakePair( wkbType, contents );
1338 }
1339 
1340 QStringList QgsGeometryUtils::wktGetChildBlocks( const QString &wkt, const QString &defaultType )
1341 {
1342  int level = 0;
1343  QString block;
1344  QStringList blocks;
1345  for ( int i = 0, n = wkt.length(); i < n; ++i )
1346  {
1347  if ( ( wkt[i].isSpace() || wkt[i] == '\n' || wkt[i] == '\t' ) && level == 0 )
1348  continue;
1349 
1350  if ( wkt[i] == ',' && level == 0 )
1351  {
1352  if ( !block.isEmpty() )
1353  {
1354  if ( block.startsWith( '(' ) && !defaultType.isEmpty() )
1355  block.prepend( defaultType + ' ' );
1356  blocks.append( block );
1357  }
1358  block.clear();
1359  continue;
1360  }
1361  if ( wkt[i] == '(' )
1362  ++level;
1363  else if ( wkt[i] == ')' )
1364  --level;
1365  block += wkt[i];
1366  }
1367  if ( !block.isEmpty() )
1368  {
1369  if ( block.startsWith( '(' ) && !defaultType.isEmpty() )
1370  block.prepend( defaultType + ' ' );
1371  blocks.append( block );
1372  }
1373  return blocks;
1374 }
1375 
1377 {
1379 
1380 
1381  double x = ( pt1.x() + pt2.x() ) / 2.0;
1382  double y = ( pt1.y() + pt2.y() ) / 2.0;
1383  double z = std::numeric_limits<double>::quiet_NaN();
1384  double m = std::numeric_limits<double>::quiet_NaN();
1385 
1386  if ( pt1.is3D() || pt2.is3D() )
1387  {
1388  pType = QgsWkbTypes::addZ( pType );
1389  z = ( pt1.z() + pt2.z() ) / 2.0;
1390  }
1391 
1392  if ( pt1.isMeasure() || pt2.isMeasure() )
1393  {
1394  pType = QgsWkbTypes::addM( pType );
1395  m = ( pt1.m() + pt2.m() ) / 2.0;
1396  }
1397 
1398  return QgsPoint( pType, x, y, z, m );
1399 }
1400 
1401 QgsPoint QgsGeometryUtils::interpolatePointOnLine( const QgsPoint &p1, const QgsPoint &p2, const double fraction )
1402 {
1403  const double _fraction = 1 - fraction;
1404  return QgsPoint( p1.wkbType(),
1405  p1.x() * _fraction + p2.x() * fraction,
1406  p1.y() * _fraction + p2.y() * fraction,
1407  p1.is3D() ? p1.z() * _fraction + p2.z() * fraction : std::numeric_limits<double>::quiet_NaN(),
1408  p1.isMeasure() ? p1.m() * _fraction + p2.m() * fraction : std::numeric_limits<double>::quiet_NaN() );
1409 }
1410 
1411 QgsPointXY QgsGeometryUtils::interpolatePointOnLine( const double x1, const double y1, const double x2, const double y2, const double fraction )
1412 {
1413  const double deltaX = ( x2 - x1 ) * fraction;
1414  const double deltaY = ( y2 - y1 ) * fraction;
1415  return QgsPointXY( x1 + deltaX, y1 + deltaY );
1416 }
1417 
1418 QgsPointXY QgsGeometryUtils::interpolatePointOnLineByValue( const double x1, const double y1, const double v1, const double x2, const double y2, const double v2, const double value )
1419 {
1420  if ( qgsDoubleNear( v1, v2 ) )
1421  return QgsPointXY( x1, y1 );
1422 
1423  const double fraction = ( value - v1 ) / ( v2 - v1 );
1424  return interpolatePointOnLine( x1, y1, x2, y2, fraction );
1425 }
1426 
1427 double QgsGeometryUtils::gradient( const QgsPoint &pt1, const QgsPoint &pt2 )
1428 {
1429  double delta_x = pt2.x() - pt1.x();
1430  double delta_y = pt2.y() - pt1.y();
1431  if ( qgsDoubleNear( delta_x, 0.0 ) )
1432  {
1433  return INFINITY;
1434  }
1435 
1436  return delta_y / delta_x;
1437 }
1438 
1439 void QgsGeometryUtils::coefficients( const QgsPoint &pt1, const QgsPoint &pt2, double &a, double &b, double &c )
1440 {
1441  if ( qgsDoubleNear( pt1.x(), pt2.x() ) )
1442  {
1443  a = 1;
1444  b = 0;
1445  c = -pt1.x();
1446  }
1447  else if ( qgsDoubleNear( pt1.y(), pt2.y() ) )
1448  {
1449  a = 0;
1450  b = 1;
1451  c = -pt1.y();
1452  }
1453  else
1454  {
1455  a = pt1.y() - pt2.y();
1456  b = pt2.x() - pt1.x();
1457  c = pt1.x() * pt2.y() - pt1.y() * pt2.x();
1458  }
1459 
1460 }
1461 
1463 {
1464  QgsLineString line;
1465  QgsPoint p2;
1466 
1467  if ( ( p == s1 ) || ( p == s2 ) )
1468  {
1469  return line;
1470  }
1471 
1472  double a, b, c;
1473  coefficients( s1, s2, a, b, c );
1474 
1475  if ( qgsDoubleNear( a, 0 ) )
1476  {
1477  p2 = QgsPoint( p.x(), s1.y() );
1478  }
1479  else if ( qgsDoubleNear( b, 0 ) )
1480  {
1481  p2 = QgsPoint( s1.x(), p.y() );
1482  }
1483  else
1484  {
1485  double y = ( -c - a * p.x() ) / b;
1486  double m = gradient( s1, s2 );
1487  double d2 = 1 + m * m;
1488  double H = p.y() - y;
1489  double dx = m * H / d2;
1490  double dy = m * dx;
1491  p2 = QgsPoint( p.x() + dx, y + dy );
1492  }
1493 
1494  line.addVertex( p );
1495  line.addVertex( p2 );
1496 
1497  return line;
1498 }
1499 
1500 double QgsGeometryUtils::lineAngle( double x1, double y1, double x2, double y2 )
1501 {
1502  double at = std::atan2( y2 - y1, x2 - x1 );
1503  double a = -at + M_PI_2;
1504  return normalizedAngle( a );
1505 }
1506 
1507 double QgsGeometryUtils::angleBetweenThreePoints( double x1, double y1, double x2, double y2, double x3, double y3 )
1508 {
1509  double angle1 = std::atan2( y1 - y2, x1 - x2 );
1510  double angle2 = std::atan2( y3 - y2, x3 - x2 );
1511  return normalizedAngle( angle1 - angle2 );
1512 }
1513 
1514 double QgsGeometryUtils::linePerpendicularAngle( double x1, double y1, double x2, double y2 )
1515 {
1516  double a = lineAngle( x1, y1, x2, y2 );
1517  a += M_PI_2;
1518  return normalizedAngle( a );
1519 }
1520 
1521 double QgsGeometryUtils::averageAngle( double x1, double y1, double x2, double y2, double x3, double y3 )
1522 {
1523  // calc average angle between the previous and next point
1524  double a1 = lineAngle( x1, y1, x2, y2 );
1525  double a2 = lineAngle( x2, y2, x3, y3 );
1526  return averageAngle( a1, a2 );
1527 }
1528 
1529 double QgsGeometryUtils::averageAngle( double a1, double a2 )
1530 {
1531  a1 = normalizedAngle( a1 );
1532  a2 = normalizedAngle( a2 );
1533  double clockwiseDiff = 0.0;
1534  if ( a2 >= a1 )
1535  {
1536  clockwiseDiff = a2 - a1;
1537  }
1538  else
1539  {
1540  clockwiseDiff = a2 + ( 2 * M_PI - a1 );
1541  }
1542  double counterClockwiseDiff = 2 * M_PI - clockwiseDiff;
1543 
1544  double resultAngle = 0;
1545  if ( clockwiseDiff <= counterClockwiseDiff )
1546  {
1547  resultAngle = a1 + clockwiseDiff / 2.0;
1548  }
1549  else
1550  {
1551  resultAngle = a1 - counterClockwiseDiff / 2.0;
1552  }
1553  return normalizedAngle( resultAngle );
1554 }
1555 
1557  const QgsVector3D &P2, const QgsVector3D &P22 )
1558 {
1559  QgsVector3D u1 = P12 - P1;
1560  QgsVector3D u2 = P22 - P2;
1561  QgsVector3D u3 = QgsVector3D::crossProduct( u1, u2 );
1562  if ( u3.length() == 0 ) return 1;
1563  u3.normalize();
1564  QgsVector3D dir = P1 - P2;
1565  return std::fabs( ( QgsVector3D::dotProduct( dir, u3 ) ) ); // u3 is already normalized
1566 }
1567 
1569  const QgsVector3D &P2, const QgsVector3D &P22,
1570  QgsVector3D &X1, double epsilon )
1571 {
1572  QgsVector3D d = P2 - P1;
1573  QgsVector3D u1 = P12 - P1;
1574  u1.normalize();
1575  QgsVector3D u2 = P22 - P2;
1576  u2.normalize();
1577  QgsVector3D u3 = QgsVector3D::crossProduct( u1, u2 );
1578 
1579  if ( std::fabs( u3.x() ) <= epsilon &&
1580  std::fabs( u3.y() ) <= epsilon &&
1581  std::fabs( u3.z() ) <= epsilon )
1582  {
1583  // The rays are almost parallel.
1584  return false;
1585  }
1586 
1587  // X1 and X2 are the closest points on lines
1588  // we want to find X1 (lies on u1)
1589  // solving the linear equation in r1 and r2: Xi = Pi + ri*ui
1590  // we are only interested in X1 so we only solve for r1.
1591  float a1 = QgsVector3D::dotProduct( u1, u1 ), b1 = QgsVector3D::dotProduct( u1, u2 ), c1 = QgsVector3D::dotProduct( u1, d );
1592  float a2 = QgsVector3D::dotProduct( u1, u2 ), b2 = QgsVector3D::dotProduct( u2, u2 ), c2 = QgsVector3D::dotProduct( u2, d );
1593  if ( !( std::fabs( b1 ) > epsilon ) )
1594  {
1595  // Denominator is close to zero.
1596  return false;
1597  }
1598  if ( !( a2 != -1 && a2 != 1 ) )
1599  {
1600  // Lines are parallel
1601  return false;
1602  }
1603 
1604  double r1 = ( c2 - b2 * c1 / b1 ) / ( a2 - b2 * a1 / b1 );
1605  X1 = P1 + u1 * r1;
1606 
1607  return true;
1608 }
1609 
1611  const QgsVector3D &Lb1, const QgsVector3D &Lb2,
1612  QgsVector3D &intersection )
1613 {
1614 
1615  // if all Vector are on the same plane (have the same Z), use the 2D intersection
1616  // else return a false result
1617  if ( qgsDoubleNear( La1.z(), La2.z() ) && qgsDoubleNear( La1.z(), Lb1.z() ) && qgsDoubleNear( La1.z(), Lb2.z() ) )
1618  {
1619  QgsPoint ptInter;
1620  bool isIntersection;
1621  segmentIntersection( QgsPoint( La1.x(), La1.y() ),
1622  QgsPoint( La2.x(), La2.y() ),
1623  QgsPoint( Lb1.x(), Lb1.y() ),
1624  QgsPoint( Lb2.x(), Lb2.y() ),
1625  ptInter,
1626  isIntersection,
1627  1e-8,
1628  true );
1629  intersection.set( ptInter.x(), ptInter.y(), La1.z() );
1630  return true;
1631  }
1632 
1633  // first check if lines have an exact intersection point
1634  // do it by checking if the shortest distance is exactly 0
1635  float distance = skewLinesDistance( La1, La2, Lb1, Lb2 );
1636  if ( qgsDoubleNear( distance, 0.0 ) )
1637  {
1638  // 3d lines have exact intersection point.
1639  QgsVector3D C = La2;
1640  QgsVector3D D = Lb2;
1641  QgsVector3D e = La1 - La2;
1642  QgsVector3D f = Lb1 - Lb2;
1643  QgsVector3D g = D - C;
1644  if ( qgsDoubleNear( ( QgsVector3D::crossProduct( f, g ) ).length(), 0.0 ) || qgsDoubleNear( ( QgsVector3D::crossProduct( f, e ) ).length(), 0.0 ) )
1645  {
1646  // Lines have no intersection, are they parallel?
1647  return false;
1648  }
1649 
1650  QgsVector3D fgn = QgsVector3D::crossProduct( f, g );
1651  fgn.normalize();
1652 
1653  QgsVector3D fen = QgsVector3D::crossProduct( f, e );
1654  fen.normalize();
1655 
1656  int di = -1;
1657  if ( fgn == fen ) // same direction?
1658  di *= -1;
1659 
1660  intersection = C + e * di * ( QgsVector3D::crossProduct( f, g ).length() / QgsVector3D::crossProduct( f, e ).length() );
1661  return true;
1662  }
1663 
1664  // try to calculate the approximate intersection point
1665  QgsVector3D X1, X2;
1666  bool firstIsDone = skewLinesProjection( La1, La2, Lb1, Lb2, X1 );
1667  bool secondIsDone = skewLinesProjection( Lb1, Lb2, La1, La2, X2 );
1668 
1669  if ( !firstIsDone || !secondIsDone )
1670  {
1671  // Could not obtain projection point.
1672  return false;
1673  }
1674 
1675  intersection = ( X1 + X2 ) / 2.0;
1676  return true;
1677 }
1678 
1679 double QgsGeometryUtils::triangleArea( double aX, double aY, double bX, double bY, double cX, double cY )
1680 {
1681  return 0.5 * std::abs( ( aX - cX ) * ( bY - aY ) - ( aX - bX ) * ( cY - aY ) );
1682 }
1683 
1684 void QgsGeometryUtils::weightedPointInTriangle( const double aX, const double aY, const double bX, const double bY, const double cX, const double cY,
1685  double weightB, double weightC, double &pointX, double &pointY )
1686 {
1687  // if point will be outside of the triangle, invert weights
1688  if ( weightB + weightC > 1 )
1689  {
1690  weightB = 1 - weightB;
1691  weightC = 1 - weightC;
1692  }
1693 
1694  const double rBx = weightB * ( bX - aX );
1695  const double rBy = weightB * ( bY - aY );
1696  const double rCx = weightC * ( cX - aX );
1697  const double rCy = weightC * ( cY - aY );
1698 
1699  pointX = rBx + rCx + aX;
1700  pointY = rBy + rCy + aY;
1701 }
1702 
1704 {
1705  bool rc = false;
1706 
1707  for ( const QgsPoint &pt : points )
1708  {
1709  if ( pt.is3D() )
1710  {
1711  point.convertTo( QgsWkbTypes::addZ( point.wkbType() ) );
1712  point.setZ( pt.z() );
1713  rc = true;
1714  break;
1715  }
1716  }
1717 
1718  return rc;
1719 }
double sqrDist(double x, double y) const SIP_HOLDGIL
Returns the squared distance between this point a specified x, y coordinate.
Definition: qgspointxy.h:175
Maximum distance between an arbitrary point on the original curve and closest point on its approximat...
3 Class for storage of 3D vectors similar to QVector3D, with the difference that it uses double preci...
Definition: qgsvector3d.h:31
static QString pointsToJSON(const QgsPointSequence &points, int precision)
Returns a geoJSON coordinates string.
int precision
static bool circleClockwise(double angle1, double angle2, double angle3) SIP_HOLDGIL
Returns true if the circle defined by three angles is ordered clockwise.
double y
Definition: qgspoint.h:42
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.
virtual bool isEmpty() const
Returns true if the geometry is empty.
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 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 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, cY).
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, angle2 and angle3.
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 QVector< SelfIntersection > selfIntersections(const QgsAbstractGeometry *geom, int part, int ring, double tolerance)
Find self intersections in a polyline.
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, y1) to (x2, y2).
bool addZValue(double zValue=0) override
Adds a z-dimension to the geometry, initialized to a preset value.
Definition: qgspoint.cpp:541
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 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 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 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.
double length() const
Returns the length of the vector.
Definition: qgsvector3d.h:112
double y
Definition: qgspointxy.h:48
A class to represent a 2D point.
Definition: qgspointxy.h:43
bool qgsDoubleNear(double a, double b, double epsilon=4 *std::numeric_limits< double >::epsilon())
Compare two doubles (but allow some difference)
Definition: qgis.h:316
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 Type parseType(const QString &wktStr)
Attempts to extract the WKB type from a WKT string.
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.
Curve polygon geometry type.
double x() const SIP_HOLDGIL
Returns the vector&#39;s x-component.
Definition: qgsvector.h:147
static QgsPoint midpoint(const QgsPoint &pt1, const QgsPoint &pt2) SIP_HOLDGIL
Returns a middle point between points pt1 and pt2.
QgsPoint project(double distance, double azimuth, double inclination=90.0) const SIP_HOLDGIL
Returns a new point which correspond to this point projected by a specified distance with specified a...
Definition: qgspoint.cpp:715
static Type addZ(Type type) SIP_HOLDGIL
Adds the z dimension to a WKB type and returns the new type.
Definition: qgswkbtypes.h:1146
static QgsLineString perpendicularSegment(const QgsPoint &p, const QgsPoint &s1, const QgsPoint &s2) SIP_HOLDGIL
Create a perpendicular line segment from p to segment [s1, s2].
SegmentationToleranceType
Segmentation tolerance as maximum angle or maximum difference between approximation and circle...
static QgsPoint closestVertex(const QgsAbstractGeometry &geom, const QgsPoint &pt, QgsVertexId &id)
Returns the closest vertex to a geometry for a specified point.
virtual bool nextVertex(QgsVertexId &id, QgsPoint &vertex) const =0
Returns next vertex id and coordinates.
double y() const
Returns Y coordinate.
Definition: qgsvector3d.h:51
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
static QStringList wktGetChildBlocks(const QString &wkt, const QString &defaultType=QString())
Parses a WKT string and returns of list of blocks contained in the WKT.
int vertex
Vertex number.
void set(double x, double y, double z)
Sets vector coordinates.
Definition: qgsvector3d.h:56
static double dotProduct(const QgsVector3D &v1, const QgsVector3D &v2)
Returns the dot product of two vectors.
Definition: qgsvector3d.h:98
static QgsVector3D crossProduct(const QgsVector3D &v1, const QgsVector3D &v2)
Returns the cross product of two vectors.
Definition: qgsvector3d.h:104
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.
int part
Part number.
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 Type addM(Type type) SIP_HOLDGIL
Adds the m dimension to a WKB type and returns the new type.
Definition: qgswkbtypes.h:1171
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, y1) to (x2, y2) and (x2, y2) to (x3, y3).
void normalize()
Normalizes the current vector in place.
Definition: qgsvector3d.h:118
bool isEmpty() const override SIP_HOLDGIL
Returns true if the geometry is empty.
Definition: qgspoint.cpp:742
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.
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
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 bool hasM(Type type) SIP_HOLDGIL
Tests whether a WKB type contains m values.
Definition: qgswkbtypes.h:1100
virtual double length() const
Returns the planar, 2-dimensional length of the geometry.
double z() const
Returns Z coordinate.
Definition: qgsvector3d.h:53
static QString pointsToWKT(const QgsPointSequence &points, int precision, bool is3D, bool isMeasure)
Returns a WKT coordinate list.
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.
Type
The WKB type describes the number of dimensions a geometry has.
Definition: qgswkbtypes.h:69
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...
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.
static bool setZValueFromPoints(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 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) ...
virtual double segmentLength(QgsVertexId startVertex) const =0
Returns the length of the segment of the geometry which begins at startVertex.
void addVertex(const QgsPoint &pt)
Adds a new vertex to the end of the line string.
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)...
Utility class for identifying a unique vertex within a geometry.
const double DEFAULT_SEGMENT_EPSILON
Default snapping tolerance for segments.
Definition: qgis.h:757
Geometry collection.
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).
QString qgsDoubleToString(double a, int precision=17)
Returns a string representation of a double.
Definition: qgis.h:276
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...
void setZ(double z) SIP_HOLDGIL
Sets the point&#39;s z-coordinate.
Definition: qgspoint.h:293
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 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.
Abstract base class for curved geometry type.
Definition: qgscurve.h:35
void set(double x, double y) SIP_HOLDGIL
Sets the x and y value of the point.
Definition: qgspointxy.h:124
Abstract base class for all geometries.
static QVector< QgsLineString * > extractLineStrings(const QgsAbstractGeometry *geom)
Returns list of linestrings extracted from the passed geometry.
Point geometry type, with support for z-dimension and m-values.
Definition: qgspoint.h:37
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.
AxisOrder
Axis order for GML generation.
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 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...
double x
Definition: qgspointxy.h:47
A class to represent a vector.
Definition: qgsvector.h:29
static double sqrDistance2D(const QgsPoint &pt1, const QgsPoint &pt2) SIP_HOLDGIL
Returns the squared 2D distance between two points.
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 json pointsToJson(const QgsPointSequence &points, int precision)
Returns coordinates as json object.
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...
QVector< QgsPoint > QgsPointSequence
static bool hasZ(Type type) SIP_HOLDGIL
Tests whether a WKB type contains the z-dimension.
Definition: qgswkbtypes.h:1050
static double normalizedAngle(double angle) SIP_HOLDGIL
Ensures that an angle is in the range 0 <= angle < 2 pi.
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 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...
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:332
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...
bool isMeasure() const SIP_HOLDGIL
Returns true if the geometry contains m values.
double qgsRound(double number, int places)
Returns a double number, rounded (as close as possible) to the specified number of places...
Definition: qgis.h:364
double distance(double x, double y) const SIP_HOLDGIL
Returns the distance between this point and a specified x, y coordinate.
Definition: qgspointxy.h:196
bool is3D() const SIP_HOLDGIL
Returns true if the geometry is 3D and contains a z-value.
static QgsPointSequence pointsFromWKT(const QString &wktCoordinateList, bool is3D, bool isMeasure)
Returns a list of points contained in a WKT string.
VertexType type
Vertex type.
Line string geometry type, with support for z-dimension and m-values.
Definition: qgslinestring.h:43
static QgsPoint pointOnLineWithDistance(const QgsPoint &startPoint, const QgsPoint &directionPoint, double distance) SIP_HOLDGIL
Returns a point a specified distance toward a second point.
virtual int vertexCount(int part=0, int ring=0) const =0
Returns the number of vertices of which this geometry is built.
double y() const SIP_HOLDGIL
Returns the vector&#39;s y-component.
Definition: qgsvector.h:156
double length() const SIP_HOLDGIL
Returns the length of the vector.
Definition: qgsvector.h:128
bool convertTo(QgsWkbTypes::Type type) override
Converts the geometry to a specified type.
Definition: qgspoint.cpp:610
QgsWkbTypes::Type wkbType() const SIP_HOLDGIL
Returns the WKB type of the geometry.
bool addMValue(double mValue=0) override
Adds a measure to the geometry, initialized to a preset value.
Definition: qgspoint.cpp:552
static void pointsToWKB(QgsWkbPtr &wkb, const QgsPointSequence &points, bool is3D, bool isMeasure)
Returns a LinearRing { uint32 numPoints; Point points[numPoints]; }.
static int circleCircleIntersections(QgsPointXY center1, double radius1, QgsPointXY center2, double radius2, QgsPointXY &intersection1, QgsPointXY &intersection2) SIP_HOLDGIL
Calculates the intersections points between the circle with center center1 and radius radius1 and the...
static bool lineIntersection(const QgsPoint &p1, QgsVector v1, const QgsPoint &p2, QgsVector v2, QgsPoint &intersection) SIP_HOLDGIL
Computes the intersection between two lines.
static double gradient(const QgsPoint &pt1, const QgsPoint &pt2) SIP_HOLDGIL
Returns the gradient of a line defined by points pt1 and pt2.
double z
Definition: qgspoint.h:43
static QgsPoint closestPoint(const QgsAbstractGeometry &geometry, const QgsPoint &point)
Returns the nearest point on a segment of a geometry for the specified point.
static double linePerpendicularAngle(double x1, double y1, double x2, double y2) SIP_HOLDGIL
Calculates the perpendicular angle to a line joining two points.
virtual QgsPoint vertexAt(QgsVertexId id) const =0
Returns the point corresponding to a specified vertex id.
static double distanceToVertex(const QgsAbstractGeometry &geom, QgsVertexId id)
Returns the distance along a geometry from its first vertex to the specified vertex.
int ring
Ring number.
double x() const
Returns X coordinate.
Definition: qgsvector3d.h:49
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.
double m
Definition: qgspoint.h:44
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.
double x
Definition: qgspoint.h:41
bool isValid() const SIP_HOLDGIL
Returns true if the vertex id is valid.