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