QGIS API Documentation 3.37.0-Master (614ebf5f6b2)
Loading...
Searching...
No Matches
qgsgeometryutils_base.cpp
Go to the documentation of this file.
1/***************************************************************************
2 qgsgeometryutils_base.cpp
3 -------------------------------------------------------------------
4Date : 14 september 2023
5Copyright : (C) 2023 by Loïc Bartoletti
6email : loic dot bartoletti at oslandia 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
17#include "qgsvector3d.h"
18#include "qgsvector.h"
19
20double QgsGeometryUtilsBase::sqrDistToLine( double ptX, double ptY, double x1, double y1, double x2, double y2, double &minDistX, double &minDistY, double epsilon )
21{
22 minDistX = x1;
23 minDistY = y1;
24
25 double dx = x2 - x1;
26 double dy = y2 - y1;
27
28 if ( !qgsDoubleNear( dx, 0.0 ) || !qgsDoubleNear( dy, 0.0 ) )
29 {
30 const double t = ( ( ptX - x1 ) * dx + ( ptY - y1 ) * dy ) / ( dx * dx + dy * dy );
31 if ( t > 1 )
32 {
33 minDistX = x2;
34 minDistY = y2;
35 }
36 else if ( t > 0 )
37 {
38 minDistX += dx * t;
39 minDistY += dy * t;
40 }
41 }
42
43 dx = ptX - minDistX;
44 dy = ptY - minDistY;
45
46 const double dist = dx * dx + dy * dy;
47
48 //prevent rounding errors if the point is directly on the segment
49 if ( qgsDoubleNear( dist, 0.0, epsilon ) )
50 {
51 minDistX = ptX;
52 minDistY = ptY;
53 return 0.0;
54 }
55
56 return dist;
57}
58
59int QgsGeometryUtilsBase::leftOfLine( const double x, const double y, const double x1, const double y1, const double x2, const double y2 )
60{
61 const double f1 = x - x1;
62 const double f2 = y2 - y1;
63 const double f3 = y - y1;
64 const double f4 = x2 - x1;
65 const double test = ( f1 * f2 - f3 * f4 );
66 // return -1, 0, or 1
67 return qgsDoubleNear( test, 0.0 ) ? 0 : ( test < 0 ? -1 : 1 );
68}
69
70void QgsGeometryUtilsBase::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 )
71{
72 const double dx = x2 - x1;
73 const double dy = y2 - y1;
74 const double length = std::sqrt( dx * dx + dy * dy );
75
76 if ( qgsDoubleNear( length, 0.0 ) )
77 {
78 x = x1;
79 y = y1;
80 if ( z && z1 )
81 *z = *z1;
82 if ( m && m1 )
83 *m = *m1;
84 }
85 else
86 {
87 const double scaleFactor = distance / length;
88 x = x1 + dx * scaleFactor;
89 y = y1 + dy * scaleFactor;
90 if ( z && z1 && z2 )
91 *z = *z1 + ( *z2 - *z1 ) * scaleFactor;
92 if ( m && m1 && m2 )
93 *m = *m1 + ( *m2 - *m1 ) * scaleFactor;
94 }
95}
96
97void QgsGeometryUtilsBase::perpendicularOffsetPointAlongSegment( double x1, double y1, double x2, double y2, double proportion, double offset, double *x, double *y )
98{
99 // calculate point along segment
100 const double mX = x1 + ( x2 - x1 ) * proportion;
101 const double mY = y1 + ( y2 - y1 ) * proportion;
102 const double pX = x1 - x2;
103 const double pY = y1 - y2;
104 double normalX = -pY;
105 double normalY = pX; //#spellok
106 const double normalLength = sqrt( ( normalX * normalX ) + ( normalY * normalY ) ); //#spellok
107 normalX /= normalLength;
108 normalY /= normalLength; //#spellok
109
110 *x = mX + offset * normalX;
111 *y = mY + offset * normalY; //#spellok
112}
113
114double QgsGeometryUtilsBase::ccwAngle( double dy, double dx )
115{
116 const double angle = std::atan2( dy, dx ) * 180 / M_PI;
117 if ( angle < 0 )
118 {
119 return 360 + angle;
120 }
121 else if ( angle > 360 )
122 {
123 return 360 - angle;
124 }
125 return angle;
126}
127
128bool QgsGeometryUtilsBase::circleClockwise( double angle1, double angle2, double angle3 )
129{
130 if ( angle3 >= angle1 )
131 {
132 return !( angle2 > angle1 && angle2 < angle3 );
133 }
134 else
135 {
136 return !( angle2 > angle1 || angle2 < angle3 );
137 }
138}
139
140bool QgsGeometryUtilsBase::circleAngleBetween( double angle, double angle1, double angle2, bool clockwise )
141{
142 if ( clockwise )
143 {
144 if ( angle2 < angle1 )
145 {
146 return ( angle <= angle1 && angle >= angle2 );
147 }
148 else
149 {
150 return ( angle <= angle1 || angle >= angle2 );
151 }
152 }
153 else
154 {
155 if ( angle2 > angle1 )
156 {
157 return ( angle >= angle1 && angle <= angle2 );
158 }
159 else
160 {
161 return ( angle >= angle1 || angle <= angle2 );
162 }
163 }
164}
165
166bool QgsGeometryUtilsBase::angleOnCircle( double angle, double angle1, double angle2, double angle3 )
167{
168 const bool clockwise = circleClockwise( angle1, angle2, angle3 );
169 return circleAngleBetween( angle, angle1, angle3, clockwise );
170}
171
172void QgsGeometryUtilsBase::circleCenterRadius( double x1, double y1, double x2, double y2, double x3, double y3, double &radius, double &centerX, double &centerY )
173{
174 double dx21, dy21, dx31, dy31, h21, h31, d;
175
176 //closed circle
177 if ( qgsDoubleNear( x1, x3 ) && qgsDoubleNear( y1, y3 ) )
178 {
179 centerX = ( x1 + x2 ) / 2.0;
180 centerY = ( y1 + y2 ) / 2.0;
181 radius = std::sqrt( std::pow( centerX - x1, 2.0 ) + std::pow( centerY - y1, 2.0 ) );
182 return;
183 }
184
185 // Using Cartesian circumcenter eguations from page https://en.wikipedia.org/wiki/Circumscribed_circle
186 dx21 = x2 - x1;
187 dy21 = y2 - y1;
188 dx31 = x3 - x1;
189 dy31 = y3 - y1;
190
191 h21 = std::pow( dx21, 2.0 ) + std::pow( dy21, 2.0 );
192 h31 = std::pow( dx31, 2.0 ) + std::pow( dy31, 2.0 );
193
194 // 2*Cross product, d<0 means clockwise and d>0 counterclockwise sweeping angle
195 d = 2 * ( dx21 * dy31 - dx31 * dy21 );
196
197 // Check colinearity, Cross product = 0
198 if ( qgsDoubleNear( std::fabs( d ), 0.0, 0.00000000001 ) )
199 {
200 radius = -1.0;
201 return;
202 }
203
204 // Calculate centroid coordinates and radius
205 centerX = x1 + ( h21 * dy31 - h31 * dy21 ) / d;
206 centerY = y1 - ( h21 * dx31 - h31 * dx21 ) / d;
207 radius = std::sqrt( std::pow( centerX - x1, 2.0 ) + std::pow( centerY - y1, 2.0 ) );
208}
209
210double QgsGeometryUtilsBase::circleLength( double x1, double y1, double x2, double y2, double x3, double y3 )
211{
212 double centerX{0.0};
213 double centerY{0.0};
214 double radius{0.0};
215 circleCenterRadius( x1, y1, x2, y2, x3, y3, radius, centerX, centerY );
216 double length = M_PI / 180.0 * radius * sweepAngle( centerX, centerY, x1, y1, x2, y2, x3, y3 );
217 if ( length < 0 )
218 {
219 length = -length;
220 }
221 return length;
222}
223
224double QgsGeometryUtilsBase::sweepAngle( double centerX, double centerY, double x1, double y1, double x2, double y2, double x3, double y3 )
225{
226 const double p1Angle = QgsGeometryUtilsBase::ccwAngle( y1 - centerY, x1 - centerX );
227 const double p2Angle = QgsGeometryUtilsBase::ccwAngle( y2 - centerY, x2 - centerX );
228 const double p3Angle = QgsGeometryUtilsBase::ccwAngle( y3 - centerY, x3 - centerX );
229
230 if ( p3Angle >= p1Angle )
231 {
232 if ( p2Angle > p1Angle && p2Angle < p3Angle )
233 {
234 return ( p3Angle - p1Angle );
235 }
236 else
237 {
238 return ( - ( p1Angle + ( 360 - p3Angle ) ) );
239 }
240 }
241 else
242 {
243 if ( p2Angle < p1Angle && p2Angle > p3Angle )
244 {
245 return ( -( p1Angle - p3Angle ) );
246 }
247 else
248 {
249 return ( p3Angle + ( 360 - p1Angle ) );
250 }
251 }
252}
253
254double QgsGeometryUtilsBase::interpolateArcValue( double angle, double a1, double a2, double a3, double zm1, double zm2, double zm3 )
255{
256 /* Counter-clockwise sweep */
257 if ( a1 < a2 )
258 {
259 if ( angle <= a2 )
260 return zm1 + ( zm2 - zm1 ) * ( angle - a1 ) / ( a2 - a1 );
261 else
262 return zm2 + ( zm3 - zm2 ) * ( angle - a2 ) / ( a3 - a2 );
263 }
264 /* Clockwise sweep */
265 else
266 {
267 if ( angle >= a2 )
268 return zm1 + ( zm2 - zm1 ) * ( a1 - angle ) / ( a1 - a2 );
269 else
270 return zm2 + ( zm3 - zm2 ) * ( a2 - angle ) / ( a2 - a3 );
271 }
272}
274{
275 double clippedAngle = angle;
276 if ( clippedAngle >= M_PI * 2 || clippedAngle <= -2 * M_PI )
277 {
278 clippedAngle = std::fmod( clippedAngle, 2 * M_PI );
279 }
280 if ( clippedAngle < 0.0 )
281 {
282 clippedAngle += 2 * M_PI;
283 }
284 return clippedAngle;
285}
286
287
288int QgsGeometryUtilsBase::closestSideOfRectangle( double right, double bottom, double left, double top, double x, double y )
289{
290 // point outside rectangle
291 if ( x <= left && y <= bottom )
292 {
293 const double dx = left - x;
294 const double dy = bottom - y;
295 if ( qgsDoubleNear( dx, dy ) )
296 return 6;
297 else if ( dx < dy )
298 return 5;
299 else
300 return 7;
301 }
302 else if ( x >= right && y >= top )
303 {
304 const double dx = x - right;
305 const double dy = y - top;
306 if ( qgsDoubleNear( dx, dy ) )
307 return 2;
308 else if ( dx < dy )
309 return 1;
310 else
311 return 3;
312 }
313 else if ( x >= right && y <= bottom )
314 {
315 const double dx = x - right;
316 const double dy = bottom - y;
317 if ( qgsDoubleNear( dx, dy ) )
318 return 4;
319 else if ( dx < dy )
320 return 5;
321 else
322 return 3;
323 }
324 else if ( x <= left && y >= top )
325 {
326 const double dx = left - x;
327 const double dy = y - top;
328 if ( qgsDoubleNear( dx, dy ) )
329 return 8;
330 else if ( dx < dy )
331 return 1;
332 else
333 return 7;
334 }
335 else if ( x <= left )
336 return 7;
337 else if ( x >= right )
338 return 3;
339 else if ( y <= bottom )
340 return 5;
341 else if ( y >= top )
342 return 1;
343
344 // point is inside rectangle
345 const double smallestX = std::min( right - x, x - left );
346 const double smallestY = std::min( top - y, y - bottom );
347 if ( smallestX < smallestY )
348 {
349 // closer to left/right side
350 if ( right - x < x - left )
351 return 3; // closest to right side
352 else
353 return 7;
354 }
355 else
356 {
357 // closer to top/bottom side
358 if ( top - y < y - bottom )
359 return 1; // closest to top side
360 else
361 return 5;
362 }
363}
364
365void QgsGeometryUtilsBase::perpendicularCenterSegment( double pointx, double pointy, double segmentPoint1x, double segmentPoint1y, double segmentPoint2x, double segmentPoint2y, double &perpendicularSegmentPoint1x, double &perpendicularSegmentPoint1y, double &perpendicularSegmentPoint2x, double &perpendicularSegmentPoint2y, double desiredSegmentLength )
366{
367 QgsVector segmentVector = QgsVector( segmentPoint2x - segmentPoint1x, segmentPoint2y - segmentPoint1y );
368 QgsVector perpendicularVector = segmentVector.perpVector();
369 if ( desiredSegmentLength != 0 )
370 {
371 perpendicularVector = perpendicularVector.normalized() * ( desiredSegmentLength ) / 2;
372 }
373
374 perpendicularSegmentPoint1x = pointx - perpendicularVector.x();
375 perpendicularSegmentPoint1y = pointy - perpendicularVector.y();
376 perpendicularSegmentPoint2x = pointx + perpendicularVector.x();
377 perpendicularSegmentPoint2y = pointy + perpendicularVector.y();
378}
379
380double QgsGeometryUtilsBase::lineAngle( double x1, double y1, double x2, double y2 )
381{
382 const double at = std::atan2( y2 - y1, x2 - x1 );
383 const double a = -at + M_PI_2;
384 return normalizedAngle( a );
385}
386
387double QgsGeometryUtilsBase::angleBetweenThreePoints( double x1, double y1, double x2, double y2, double x3, double y3 )
388{
389 const double angle1 = std::atan2( y1 - y2, x1 - x2 );
390 const double angle2 = std::atan2( y3 - y2, x3 - x2 );
391 return normalizedAngle( angle1 - angle2 );
392}
393
394double QgsGeometryUtilsBase::linePerpendicularAngle( double x1, double y1, double x2, double y2 )
395{
396 double a = lineAngle( x1, y1, x2, y2 );
397 a += M_PI_2;
398 return normalizedAngle( a );
399}
400
401double QgsGeometryUtilsBase::averageAngle( double x1, double y1, double x2, double y2, double x3, double y3 )
402{
403 // calc average angle between the previous and next point
404 const double a1 = lineAngle( x1, y1, x2, y2 );
405 const double a2 = lineAngle( x2, y2, x3, y3 );
406 return averageAngle( a1, a2 );
407}
408
409double QgsGeometryUtilsBase::averageAngle( double a1, double a2 )
410{
411 a1 = normalizedAngle( a1 );
412 a2 = normalizedAngle( a2 );
413 double clockwiseDiff = 0.0;
414 if ( a2 >= a1 )
415 {
416 clockwiseDiff = a2 - a1;
417 }
418 else
419 {
420 clockwiseDiff = a2 + ( 2 * M_PI - a1 );
421 }
422 const double counterClockwiseDiff = 2 * M_PI - clockwiseDiff;
423
424 double resultAngle = 0;
425 if ( clockwiseDiff <= counterClockwiseDiff )
426 {
427 resultAngle = a1 + clockwiseDiff / 2.0;
428 }
429 else
430 {
431 resultAngle = a1 - counterClockwiseDiff / 2.0;
432 }
433 return normalizedAngle( resultAngle );
434}
435
437 const QgsVector3D &P2, const QgsVector3D &P22 )
438{
439 const QgsVector3D u1 = P12 - P1;
440 const QgsVector3D u2 = P22 - P2;
442 if ( u3.length() == 0 ) return 1;
443 u3.normalize();
444 const QgsVector3D dir = P1 - P2;
445 return std::fabs( ( QgsVector3D::dotProduct( dir, u3 ) ) ); // u3 is already normalized
446}
447
449 const QgsVector3D &P2, const QgsVector3D &P22,
450 QgsVector3D &X1, double epsilon )
451{
452 const QgsVector3D d = P2 - P1;
453 QgsVector3D u1 = P12 - P1;
454 u1.normalize();
455 QgsVector3D u2 = P22 - P2;
456 u2.normalize();
457 const QgsVector3D u3 = QgsVector3D::crossProduct( u1, u2 );
458
459 if ( std::fabs( u3.x() ) <= epsilon &&
460 std::fabs( u3.y() ) <= epsilon &&
461 std::fabs( u3.z() ) <= epsilon )
462 {
463 // The rays are almost parallel.
464 return false;
465 }
466
467 // X1 and X2 are the closest points on lines
468 // we want to find X1 (lies on u1)
469 // solving the linear equation in r1 and r2: Xi = Pi + ri*ui
470 // we are only interested in X1 so we only solve for r1.
471 const double a1 = QgsVector3D::dotProduct( u1, u1 ), b1 = QgsVector3D::dotProduct( u1, u2 ), c1 = QgsVector3D::dotProduct( u1, d );
472 const double a2 = QgsVector3D::dotProduct( u1, u2 ), b2 = QgsVector3D::dotProduct( u2, u2 ), c2 = QgsVector3D::dotProduct( u2, d );
473 if ( !( std::fabs( b1 ) > epsilon ) )
474 {
475 // Denominator is close to zero.
476 return false;
477 }
478 if ( !( a2 != -1 && a2 != 1 ) )
479 {
480 // Lines are parallel
481 return false;
482 }
483
484 const double r1 = ( c2 - b2 * c1 / b1 ) / ( a2 - b2 * a1 / b1 );
485 X1 = P1 + u1 * r1;
486
487 return true;
488}
489
490bool QgsGeometryUtilsBase::lineIntersection( double p1x, double p1y, QgsVector v1, double p2x, double p2y, QgsVector v2, double &intersectionX, double &intersectionY )
491{
492 const double d = v1.y() * v2.x() - v1.x() * v2.y();
493
494 if ( qgsDoubleNear( d, 0 ) )
495 return false;
496
497 const double dx = p2x - p1x;
498 const double dy = p2y - p1y;
499 const double k = ( dy * v2.x() - dx * v2.y() ) / d;
500
501 intersectionX = p1x + v1.x() * k;
502 intersectionY = p1y + v1.y() * k;
503
504 return true;
505}
506
507static bool equals( double p1x, double p1y, double p2x, double p2y, double epsilon = 1e-8 )
508{
509 return qgsDoubleNear( p1x, p2x, epsilon ) && qgsDoubleNear( p1y, p2y, epsilon );
510}
511
512bool QgsGeometryUtilsBase::segmentIntersection( double p1x, double p1y, double p2x, double p2y, double q1x, double q1y, double q2x, double q2y, double &intersectionPointX, double &intersectionPointY, bool &isIntersection, double tolerance, bool acceptImproperIntersection )
513{
514 isIntersection = false;
515 intersectionPointX = intersectionPointY = std::numeric_limits<double>::quiet_NaN();
516
517 QgsVector v( p2x - p1x, p2y - p1y );
518 QgsVector w( q2x - q1x, q2y - q1y );
519 const double vl = v.length();
520 const double wl = w.length();
521
522 if ( qgsDoubleNear( vl, 0.0, tolerance ) || qgsDoubleNear( wl, 0.0, tolerance ) )
523 {
524 return false;
525 }
526 v = v / vl;
527 w = w / wl;
528
529 if ( !lineIntersection( p1x, p1y, v, q1x, q1y, w, intersectionPointX, intersectionPointY ) )
530 {
531 return false;
532 }
533
534 isIntersection = true;
535 if ( acceptImproperIntersection )
536 {
537 if ( ( equals( p1x, p1y, q1x, q1y ) ) || ( equals( p1x, p1y, q2x, q2y ) ) )
538 {
539 intersectionPointX = p1x;
540 intersectionPointY = p1y;
541 return true;
542 }
543 else if ( ( equals( p1x, p1y, q2x, q2y ) ) || ( equals( p2x, p2y, q2x, q2y ) ) )
544 {
545 intersectionPointX = p2x;
546 intersectionPointY = p2y;
547 return true;
548 }
549
550 double x, y;
551 if (
552 // intersectionPoint = p1
553 qgsDoubleNear( sqrDistToLine( p1x, p1y, q1x, q1y, q2x, q2y, x, y, tolerance ), 0.0, tolerance ) ||
554 // intersectionPoint = p2
555 qgsDoubleNear( sqrDistToLine( p2x, p2y, q1x, q1y, q2x, q2y, x, y, tolerance ), 0.0, tolerance ) ||
556 // intersectionPoint = q1
557 qgsDoubleNear( sqrDistToLine( q1x, q1y, p1x, p1y, p2x, p2y, x, y, tolerance ), 0.0, tolerance ) ||
558 // intersectionPoint = q2
559 qgsDoubleNear( sqrDistToLine( q2x, q2y, p1x, p1y, p2x, p2y, x, y, tolerance ), 0.0, tolerance )
560 )
561 {
562 return true;
563 }
564 }
565
566 const double lambdav = QgsVector( intersectionPointX - p1x, intersectionPointY - p1y ) * v;
567 if ( lambdav < 0. + tolerance || lambdav > vl - tolerance )
568 return false;
569
570 const double lambdaw = QgsVector( intersectionPointX - q1x, intersectionPointY - q1y ) * w;
571 return !( lambdaw < 0. + tolerance || lambdaw >= wl - tolerance );
572
573}
574
575
577 const QgsVector3D &Lb1, const QgsVector3D &Lb2,
578 QgsVector3D &intersection )
579{
580
581 // if all Vector are on the same plane (have the same Z), use the 2D intersection
582 // else return a false result
583 if ( qgsDoubleNear( La1.z(), La2.z() ) && qgsDoubleNear( La1.z(), Lb1.z() ) && qgsDoubleNear( La1.z(), Lb2.z() ) )
584 {
585 double ptInterX = 0.0, ptInterY = 0.0;
586 bool isIntersection = false;
587 segmentIntersection( La1.x(), La1.y(),
588 La2.x(), La2.y(),
589 Lb1.x(), Lb1.y(),
590 Lb2.x(), Lb2.y(),
591 ptInterX, ptInterY,
592 isIntersection,
593 1e-8,
594 true );
595 intersection.set( ptInterX, ptInterY, La1.z() );
596 return true;
597 }
598
599 // first check if lines have an exact intersection point
600 // do it by checking if the shortest distance is exactly 0
601 const double distance = skewLinesDistance( La1, La2, Lb1, Lb2 );
602 if ( qgsDoubleNear( distance, 0.0 ) )
603 {
604 // 3d lines have exact intersection point.
605 const QgsVector3D C = La2;
606 const QgsVector3D D = Lb2;
607 const QgsVector3D e = La1 - La2;
608 const QgsVector3D f = Lb1 - Lb2;
609 const QgsVector3D g = D - C;
610 if ( qgsDoubleNear( ( QgsVector3D::crossProduct( f, g ) ).length(), 0.0 ) || qgsDoubleNear( ( QgsVector3D::crossProduct( f, e ) ).length(), 0.0 ) )
611 {
612 // Lines have no intersection, are they parallel?
613 return false;
614 }
615
617 fgn.normalize();
618
620 fen.normalize();
621
622 int di = -1;
623 if ( fgn == fen ) // same direction?
624 di *= -1;
625
626 intersection = C + e * di * ( QgsVector3D::crossProduct( f, g ).length() / QgsVector3D::crossProduct( f, e ).length() );
627 return true;
628 }
629
630 // try to calculate the approximate intersection point
631 QgsVector3D X1, X2;
632 const bool firstIsDone = skewLinesProjection( La1, La2, Lb1, Lb2, X1 );
633 const bool secondIsDone = skewLinesProjection( Lb1, Lb2, La1, La2, X2 );
634
635 if ( !firstIsDone || !secondIsDone )
636 {
637 // Could not obtain projection point.
638 return false;
639 }
640
641 intersection = ( X1 + X2 ) / 2.0;
642 return true;
643}
644
645double QgsGeometryUtilsBase::triangleArea( double aX, double aY, double bX, double bY, double cX, double cY )
646{
647 return 0.5 * std::abs( ( aX - cX ) * ( bY - aY ) - ( aX - bX ) * ( cY - aY ) );
648}
649
650double QgsGeometryUtilsBase::pointFractionAlongLine( double x1, double y1, double x2, double y2, double px, double py )
651{
652 const double dxp = px - x1;
653 const double dyp = py - y1;
654
655 const double dxl = x2 - x1;
656 const double dyl = y2 - y1;
657
658 return std::sqrt( ( dxp * dxp ) + ( dyp * dyp ) ) / std::sqrt( ( dxl * dxl ) + ( dyl * dyl ) );
659}
660
661void QgsGeometryUtilsBase::weightedPointInTriangle( const double aX, const double aY, const double bX, const double bY, const double cX, const double cY,
662 double weightB, double weightC, double &pointX, double &pointY )
663{
664 // if point will be outside of the triangle, invert weights
665 if ( weightB + weightC > 1 )
666 {
667 weightB = 1 - weightB;
668 weightC = 1 - weightC;
669 }
670
671 const double rBx = weightB * ( bX - aX );
672 const double rBy = weightB * ( bY - aY );
673 const double rCx = weightC * ( cX - aX );
674 const double rCy = weightC * ( cY - aY );
675
676 pointX = rBx + rCx + aX;
677 pointY = rBy + rCy + aY;
678}
679
680bool QgsGeometryUtilsBase::pointsAreCollinear( double x1, double y1, double x2, double y2, double x3, double y3, double epsilon )
681{
682 return qgsDoubleNear( x1 * ( y2 - y3 ) + x2 * ( y3 - y1 ) + x3 * ( y1 - y2 ), 0, epsilon );
683};
684
685double QgsGeometryUtilsBase::azimuth( double x1, double y1, double x2, double y2 )
686{
687 const double dx = x2 - x1;
688 const double dy = y2 - y1;
689 return ( std::atan2( dx, dy ) * 180.0 / M_PI );
690}
691
692bool QgsGeometryUtilsBase::angleBisector( double aX, double aY, double bX, double bY, double cX, double cY, double dX, double dY,
693 double &pointX, double &pointY, double &angle )
694{
695 angle = ( azimuth( aX, aY, bX, bY ) + azimuth( cX, cY, dX, dY ) ) / 2.0;
696
697 bool intersection = false;
698 QgsGeometryUtilsBase::segmentIntersection( aX, aY, bX, bY, cX, cY, dX, dY, pointX, pointY, intersection );
699
700 return intersection;
701}
702
703void QgsGeometryUtilsBase::project( double aX, double aY, double aZ, double distance, double azimuth, double inclination, double &resultX, double &resultY, double &resultZ )
704{
705 const double radsXy = azimuth * M_PI / 180.0;
706 double dx = 0.0, dy = 0.0, dz = 0.0;
707
708 inclination = std::fmod( inclination, 360.0 );
709
710 if ( std::isnan( aZ ) && qgsDoubleNear( inclination, 90.0 ) )
711 {
712 dx = distance * std::sin( radsXy );
713 dy = distance * std::cos( radsXy );
714 }
715 else
716 {
717 const double radsZ = inclination * M_PI / 180.0;
718 dx = distance * std::sin( radsZ ) * std::sin( radsXy );
719 dy = distance * std::sin( radsZ ) * std::cos( radsXy );
720 dz = distance * std::cos( radsZ );
721 }
722
723 resultX = aX + dx;
724 resultY = aY + dy;
725 resultZ = aZ + dz;
726}
727
728bool QgsGeometryUtilsBase::bisector( double aX, double aY, double bX, double bY, double cX, double cY,
729 double &pointX, double &pointY )
730{
731 const double angle = ( azimuth( aX, aY, bX, bY ) + azimuth( aX, aY, cX, cY ) ) / 2.0;
732
733 bool intersection = false;
734 double dX = 0.0, dY = 0.0, dZ = 0.0;
735 project( aX, aY, std::numeric_limits<double>::quiet_NaN(), 1.0, angle, 90.0, dX, dY, dZ );
736 segmentIntersection( bX, bY, cX, cY, aX, aY, dX, dY, pointX, pointY, intersection );
737
738 return intersection;
739}
static double circleLength(double x1, double y1, double x2, double y2, double x3, double y3)
Length of a circular string segment defined by pt1, pt2, pt3.
static void pointOnLineWithDistance(double x1, double y1, double x2, double y2, double distance, double &x, double &y, double *z1=nullptr, double *z2=nullptr, double *z=nullptr, double *m1=nullptr, double *m2=nullptr, double *m=nullptr)
Calculates the point a specified distance from (x1, y1) toward a second point (x2,...
static double angleBetweenThreePoints(double x1, double y1, double x2, double y2, double x3, double y3)
Calculates the angle between the lines AB and BC, where AB and BC described by points a,...
static double ccwAngle(double dy, double dx)
Returns the counter clockwise angle between a line with components dx, dy and the line with dx > 0 an...
static bool circleClockwise(double angle1, double angle2, double angle3)
Returns true if the circle defined by three angles is ordered clockwise.
static double sweepAngle(double centerX, double centerY, double x1, double y1, double x2, double y2, double x3, double y3)
Calculates angle of a circular string part defined by pt1, pt2, pt3.
static double linePerpendicularAngle(double x1, double y1, double x2, double y2)
Calculates the perpendicular angle to a line joining two points.
static bool bisector(double aX, double aY, double bX, double bY, double cX, double cY, double &pointX, double &pointY)
Returns the point (pointX, pointY) forming the bisector from point (aX, aY) to the segment (bX,...
static double lineAngle(double x1, double y1, double x2, double y2)
Calculates the direction of line joining two points in radians, clockwise from the north direction.
static double averageAngle(double x1, double y1, double x2, double y2, double x3, double y3)
Calculates the average angle (in radians) between the two linear segments from (x1,...
static bool skewLinesProjection(const QgsVector3D &P1, const QgsVector3D &P12, const QgsVector3D &P2, const QgsVector3D &P22, QgsVector3D &X1, double epsilon=0.0001)
A method to project one skew line onto another.
static int closestSideOfRectangle(double right, double bottom, double left, double top, double x, double y)
Returns a number representing the closest side of a rectangle defined by /a right,...
static void project(double aX, double aY, double aZ, double distance, double azimuth, double inclination, double &resultX, double &resultY, double &resultZ)
Returns coordinates of a point which corresponds to this point projected by a specified distance with...
static double interpolateArcValue(double angle, double a1, double a2, double a3, double zm1, double zm2, double zm3)
Interpolate a value at given angle on circular arc given values (zm1, zm2, zm3) at three different an...
static void perpendicularOffsetPointAlongSegment(double x1, double y1, double x2, double y2, double proportion, double offset, double *x, double *y)
Calculates a point a certain proportion of the way along the segment from (x1, y1) to (x2,...
static double skewLinesDistance(const QgsVector3D &P1, const QgsVector3D &P12, const QgsVector3D &P2, const QgsVector3D &P22)
An algorithm to calculate the shortest distance between two skew lines.
static double pointFractionAlongLine(double x1, double y1, double x2, double y2, double px, double py)
Given the line (x1, y1) to (x2, y2) and a point (px, py) returns the fraction of the line length at w...
static double normalizedAngle(double angle)
Ensures that an angle is in the range 0 <= angle < 2 pi.
static void weightedPointInTriangle(double aX, double aY, double bX, double bY, double cX, double cY, double weightB, double weightC, double &pointX, double &pointY)
Returns a weighted point inside the triangle denoted by the points (aX, aY), (bX, bY) and (cX,...
static void perpendicularCenterSegment(double centerPointX, double centerPointY, double segmentPoint1x, double segmentPoint1y, double segmentPoint2x, double segmentPoint2y, double &perpendicularSegmentPoint1x, double &perpendicularSegmentPoint1y, double &perpendicularSegmentPoint2x, double &perpendicularSegmentPoint2y, double segmentLength=0)
Create a perpendicular line segment to a given segment [segmentPoint1,segmentPoint2] with its center ...
static bool circleAngleBetween(double angle, double angle1, double angle2, bool clockwise)
Returns true if, in a circle, angle is between angle1 and angle2.
static double triangleArea(double aX, double aY, double bX, double bY, double cX, double cY)
Returns the area of the triangle denoted by the points (aX, aY), (bX, bY) and (cX,...
static bool angleBisector(double aX, double aY, double bX, double bY, double cX, double cY, double dX, double dY, double &pointX, double &pointY, double &angle)
Returns the point (pointX, pointY) forming the bisector from segment (aX aY) (bX bY) and segment (bX,...
static double sqrDistToLine(double ptX, double ptY, double x1, double y1, double x2, double y2, double &minDistX, double &minDistY, double epsilon)
Returns the squared distance between a point and a line.
static bool angleOnCircle(double angle, double angle1, double angle2, double angle3)
Returns true if an angle is between angle1 and angle3 on a circle described by angle1,...
static bool segmentIntersection(double p1x, double p1y, double p2x, double p2y, double q1x, double q1y, double q2x, double q2y, double &intersectionPointX, double &intersectionPointY, bool &isIntersection, double tolerance=1e-8, bool acceptImproperIntersection=false)
Compute the intersection between two segments.
static int leftOfLine(const double x, const double y, const double x1, const double y1, const double x2, const double y2)
Returns a value < 0 if the point (x, y) is left of the line from (x1, y1) -> (x2, y2).
static bool linesIntersection3D(const QgsVector3D &La1, const QgsVector3D &La2, const QgsVector3D &Lb1, const QgsVector3D &Lb2, QgsVector3D &intersection)
An algorithm to calculate an (approximate) intersection of two lines in 3D.
static bool pointsAreCollinear(double x1, double y1, double x2, double y2, double x3, double y3, double epsilon)
Given the points (x1, y1), (x2, y2) and (x3, y3) returns true if these points can be considered colli...
static void circleCenterRadius(double x1, double y1, double x2, double y2, double x3, double y3, double &radius, double &centerX, double &centerY)
Returns radius and center of the circle through (x1 y1), (x2 y2), (x3 y3)
static bool lineIntersection(double p1x, double p1y, QgsVector v1, double p2x, double p2y, QgsVector v2, double &intersectionX, double &intersectionY)
Computes the intersection between two lines.
static double azimuth(double x1, double y1, double x2, double y2)
Calculates Cartesian azimuth between points (x1, y1) and (x2, y2) (clockwise in degree,...
Class for storage of 3D vectors similar to QVector3D, with the difference that it uses double precisi...
Definition qgsvector3d.h:31
double y() const
Returns Y coordinate.
Definition qgsvector3d.h:50
double z() const
Returns Z coordinate.
Definition qgsvector3d.h:52
static double dotProduct(const QgsVector3D &v1, const QgsVector3D &v2)
Returns the dot product of two vectors.
double x() const
Returns X coordinate.
Definition qgsvector3d.h:48
void normalize()
Normalizes the current vector in place.
static QgsVector3D crossProduct(const QgsVector3D &v1, const QgsVector3D &v2)
Returns the cross product of two vectors.
void set(double x, double y, double z)
Sets vector coordinates.
Definition qgsvector3d.h:73
double length() const
Returns the length of the vector.
A class to represent a vector.
Definition qgsvector.h:30
double y() const
Returns the vector's y-component.
Definition qgsvector.h:152
QgsVector normalized() const
Returns the vector's normalized (or "unit") vector (ie same angle but length of 1....
Definition qgsvector.cpp:28
QgsVector perpVector() const
Returns the perpendicular vector to this vector (rotated 90 degrees counter-clockwise)
Definition qgsvector.h:160
double x() const
Returns the vector's x-component.
Definition qgsvector.h:143
double length() const
Returns the length of the vector.
Definition qgsvector.h:124
bool qgsDoubleNear(double a, double b, double epsilon=4 *std::numeric_limits< double >::epsilon())
Compare two doubles (but allow some difference)
Definition qgis.h:5445