QGIS API Documentation 3.39.0-Master (d0dedde5474)
Searching...
No Matches
MathUtils.cpp
Go to the documentation of this file.
1/***************************************************************************
2 MathUtils.cpp
3 -------------
4 copyright : (C) 2004 by Marco Hugentobler
6 ***************************************************************************/
7
8/***************************************************************************
9 * *
10 * This program is free software; you can redistribute it and/or modify *
12 * the Free Software Foundation; either version 2 of the License, or *
13 * (at your option) any later version. *
14 * *
15 ***************************************************************************/
16
17#include "MathUtils.h"
18#include "qgslogger.h"
19#include "qgspoint.h"
20#include "Vector3D.h"
22
23bool MathUtils::calcBarycentricCoordinates( double x, double y, QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, QgsPoint *result )
24{
25 if ( p1 && p2 && p3 && result )
26 {
27 QgsPoint p( x, y, 0 );
28 const double area = triArea( p1, p2, p3 );
29 if ( area == 0 )//p1, p2, p3 are in a line
30 {
31 return false;
32 }
33 const double area1 = triArea( &p, p2, p3 );
34 const double area2 = triArea( p1, &p, p3 );
35 const double area3 = triArea( p1, p2, &p );
36 const double u = area1 / area;
37 const double v = area2 / area;
38 const double w = area3 / area;
39 result->setX( u );
40 result->setY( v );
41 result->setZ( w );
42 return true;
43 }
44 else//null pointer
45 {
46 QgsDebugError( QStringLiteral( "warning, null pointer" ) );
47 return false;
48 }
49}
50
51bool MathUtils::BarycentricToXY( double u, double v, double w, QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, QgsPoint *result )//this is wrong at the moment. Furthermore, the case, where the denominators are 0 have to be treated (two ways of calculating px and py)
52{
53 Q_UNUSED( w )
54
55 double px, py;
56 if ( p1 && p2 && p3 && result )
57 {
58 const double area = triArea( p1, p2, p3 );
59
60 if ( area == 0 )
61 {
62 QgsDebugError( QStringLiteral( "warning, p1, p2 and p3 are in a line" ) );
63 return false;
64 }
65
66 const double denominator = ( ( p2->y() - p3->y() ) * ( p1->x() - p3->x() ) - ( p3->y() - p1->y() ) * ( p3->x() - p2->x() ) );
67 if ( denominator != 0 )//drop out py in the two equations
68 {
69 px = ( 2 * u * area * ( p1->x() - p3->x() ) - 2 * v * area * ( p3->x() - p2->x() ) - p2->x() * p3->y() * ( p1->x() - p3->x() ) + p3->x() * p1->y() * ( p3->x() - p2->x() ) + p3->x() * p2->y() * ( p1->x() - p3->x() ) - p1->x() * p3->y() * ( p3->x() - p2->x() ) ) / denominator;
70 if ( ( p3->x() - p2->x() ) != 0 )
71 {
72 py = ( 2 * u * area - px * ( p2->y() - p3->y() ) - p2->x() * p3->y() + p3->x() * p2->y() ) / ( p3->x() - p2->x() );
73 }
74 else
75 {
76 py = ( 2 * v * area - px * ( p3->y() - p1->y() ) - p3->x() * p1->y() + p1->x() * p3->y() ) / ( p1->x() - p3->x() );
77 }
78 }
79 else//drop out px in the two equations(maybe this possibility occurs only, if p1, p2 and p3 are coplanar
80 {
81 py = ( 2 * u * area * ( p3->y() - p1->y() ) - 2 * v * area * ( p2->y() - p3->y() ) - p2->x() * p3->y() * ( p3->y() - p1->y() ) + p3->x() * p1->y() * ( p2->y() - p3->y() ) + p3->x() * p2->y() * ( p3->y() - p1->y() ) - p1->x() * p3->y() * ( p2->y() - p3->y() ) ) / ( ( p3->x() - p2->x() ) * ( p3->y() - p1->y() ) - ( p1->x() - p3->x() ) * ( p2->y() - p3->y() ) );
82 if ( ( p2->y() - p3->y() ) != 0 )
83 {
84 px = ( 2 * u * area - py * ( p3->x() - p2->x() ) - p2->x() * p3->y() + p3->x() * p2->y() ) / ( p2->y() - p3->y() );
85 }
86 else
87 {
88 px = ( 2 * v * area - py * ( p1->x() - p3->x() ) - p3->x() * p1->y() + p1->x() * p3->y() ) / ( p3->y() - p1->y() );
89 }
90 }
91 result->setX( px );
92 result->setY( py );
93 return true;
94 }
95 else//null pointer
96 {
97 QgsDebugError( QStringLiteral( "warning, null pointer" ) );
98 return false;
99 }
100}
101
102double MathUtils::calcBernsteinPoly( int n, int i, double t )
103{
104 if ( i < 0 )
105 {
106 return 0;
107 }
108
109 return lower( n, i ) * std::pow( t, i ) * std::pow( ( 1 - t ), ( n - i ) );
110}
111
112double MathUtils::cFDerBernsteinPoly( int n, int i, double t )
113{
114 return n * ( calcBernsteinPoly( n - 1, i - 1, t ) - calcBernsteinPoly( n - 1, i, t ) );
115}
116
117bool MathUtils::circumcenter( QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, QgsPoint *result )//version using the property that the distances from p1, p2, p3 to the circumcenter have to be equal. Possibly there is a bug
118{
119 if ( p1 && p2 && p3 && result )
120 {
121 const double distp1p2 = p1->distance( *p2 );
122 const double distp2p3 = p2->distance( *p3 );
123 if ( distp1p2 > distp2p3 )
124 {
125 //swap p1 and p3 to avoid round-off errors
126 QgsPoint *temp = p1;
127 p1 = p3;
128 p3 = temp;
129 }
130 const double denominator = -p3->x() * p2->y() + p3->x() * p1->y() + p1->x() * p2->y() + p2->x() * p3->y() - p2->x() * p1->y() - p1->x() * p3->y();
131 //if one of the denominator is zero we will have problems
132 if ( denominator == 0 )
133 {
134 QgsDebugError( QStringLiteral( "error: the three points are in a line" ) );
135 return false;
136 }
137 else
138 {
139 result->setX( 0.5 * ( p1->x()*p1->x()*p2->y() - p1->x()*p1->x()*p3->y() - p3->x()*p3->x()*p2->y() - p1->y()*p2->x()*p2->x() - p1->y()*p1->y()*p3->y() - p3->y()*p3->y()*p2->y() + p1->y()*p1->y()*p2->y() + p3->y()*p2->x()*p2->x() - p1->y()*p2->y()*p2->y() + p1->y()*p3->y()*p3->y() + p1->y()*p3->x()*p3->x() + p3->y()*p2->y()*p2->y() ) / denominator );
140 result->setY( -0.5 * ( p3->x()*p2->x()*p2->x() + p2->x()*p1->y()*p1->y() + p3->x()*p2->y()*p2->y() - p3->x()*p1->x()*p1->x() + p1->x()*p3->y()*p3->y() - p3->x()*p1->y()*p1->y() - p1->x()*p2->x()*p2->x() - p2->x()*p3->y()*p3->y() - p1->x()*p2->y()*p2->y() - p2->x()*p3->x()*p3->x() + p1->x()*p3->x()*p3->x() + p2->x()*p1->x()*p1->x() ) / denominator );
141
142 return true;
143 }
144 }
145 else
146 {
147 QgsDebugError( QStringLiteral( "warning, null pointer" ) );
148 return false;
149 }
150}
151
153{
154 if ( thepoint && p1 && p2 )
155 {
156 Vector3D normal( 0, 0, 0 );
157 Vector3D line( p2->x() - p1->x(), p2->y() - p1->y(), 0 );
158 MathUtils::normalLeft( &line, &normal, 1 );
159 const double a = normal.getX();
160 const double b = normal.getY();
161 const double c = -( normal.getX() * p2->x() + normal.getY() * p2->y() );
162 const double distance = std::fabs( ( a * thepoint->x() + b * thepoint->y() + c ) / ( std::sqrt( a * a + b * b ) ) );
163 return distance;
164 }
165 else
166 {
167 QgsDebugError( QStringLiteral( "warning, null pointer" ) );
168 return 0;
169 }
170}
171
173{
174 return std::tgamma( n + 1 );
175}
176
178{
179 const double tolerance = 0.0001;//if the amount of aValue is below this, testp is approximately on the circle and we have to define another criterion to tell, if it is inside or outside.
180
181 if ( testp && p1 && p2 && p3 )
182 {
183 double ax = p1->x();
184 double ay = p1->y();
185 double bx = p2->x();
186 double by = p2->y();
187 double cx = p3->x();
188 double cy = p3->y();
189 double px = testp->x();
190 double py = testp->y();
191
192 const double xmin = std::min( std::min( ax, px ), std::min( bx, cx ) );
193 const double ymin = std::min( std::min( ay, py ), std::min( by, cy ) );
194 ax -= xmin;
195 bx -= xmin;
196 cx -= xmin;
197 px -= xmin;
198 ay -= ymin;
199 by -= ymin;
200 cy -= ymin;
201 py -= ymin;
202
203 double aValue = ( ax * ax + ay * ay ) * triArea( p2, p3, testp );
204 aValue = aValue - ( ( bx * bx + by * by ) * triArea( p1, p3, testp ) );
205 aValue = aValue + ( ( cx * cx + cy * cy ) * triArea( p1, p2, testp ) );
206 aValue = aValue - ( ( px * px + py * py ) * triArea( p1, p2, p3 ) );
207
208 return aValue > tolerance;
209 }
210 else
211 {
212 QgsDebugError( QStringLiteral( "warning, null pointer" ) );
213 return false;
214 }
215}
216
218{
219 return angle( p1, point, p2, point ) > 90;
220}
221
222double MathUtils::leftOf( const QgsPoint &thepoint, const QgsPoint *p1, const QgsPoint *p2 )
223{
224 if ( p1 && p2 )
225 {
226 //just for debugging
227 const double f1 = thepoint.x() - p1->x();
228 const double f2 = p2->y() - p1->y();
229 const double f3 = thepoint.y() - p1->y();
230 const double f4 = p2->x() - p1->x();
231 return f1 * f2 - f3 * f4;
232 //return thepoint->getX()-p1->getX())*(p2->getY()-p1->getY())-(thepoint->getY()-p1->getY())*(p2->getX()-p1->getX();//calculating the vectorproduct
233 }
234 else
235 {
236 QgsDebugError( QStringLiteral( "warning, null pointer" ) );
237 return 0;
238 }
239}
240
242{
243 if ( p1 && p2 && p3 && p4 )
244 {
245 double t1, t2;
246 const Vector3D p1p2( p2->x() - p1->x(), p2->y() - p1->y(), 0 );
247 const Vector3D p3p4( p4->x() - p3->x(), p4->y() - p3->y(), 0 );
248 if ( ( p3p4.getX()*p1p2.getY() - p3p4.getY()*p1p2.getX() ) != 0 && p1p2.getX() != 0 ) //avoid division through zero
249 {
250 t2 = ( p1->x() * p1p2.getY() - p1->y() * p1p2.getX() + p3->y() * p1p2.getX() - p3->x() * p1p2.getY() ) / ( p3p4.getX() * p1p2.getY() - p3p4.getY() * p1p2.getX() );
251 t1 = ( p3->x() - p1->x() + t2 * p3p4.getX() ) / p1p2.getX();
252 }
253 else if ( ( p1p2.getX()*p3p4.getY() - p1p2.getY()*p3p4.getX() ) != 0 && p3p4.getX() != 0 )
254 {
255 t1 = ( p3->x() * p3p4.getY() - p3->y() * p3p4.getX() - p1->x() * p3p4.getY() + p1->y() * p3p4.getX() ) / ( p1p2.getX() * p3p4.getY() - p1p2.getY() * p3p4.getX() );
256 t2 = ( p1->x() + t1 * p1p2.getX() - p3->x() ) / p3p4.getX();
257 }
258 else//the lines are parallel
259 {
260 return false;
261 }
262
263 if ( t1 > 0 && t1 < 1 && t2 > 0 && t2 < 1 )
264 {
265 if ( ( *p1 ) == ( *p3 ) || ( *p1 ) == ( *p4 ) || ( *p2 ) == ( *p3 ) || ( *p2 ) == ( *p4 ) ) //the lines touch each other, so they do not cross
266 {
267 return false;
268 }
269 return true;
270 }
271 else
272 {
273 return false;
274 }
275 }
276
277 else
278 {
279 QgsDebugError( QStringLiteral( "warning, null pointer" ) );
280 return false;
281 }
282}
283
284bool MathUtils::lineIntersection( QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, QgsPoint *p4, QgsPoint *intersection_point )
285{
286 if ( p1 && p2 && p3 && p4 )
287 {
288 double t1, t2;
289 const Vector3D p1p2( p2->x() - p1->x(), p2->y() - p1->y(), 0 );
290 const Vector3D p3p4( p4->x() - p3->x(), p4->y() - p3->y(), 0 );
291 if ( ( p3p4.getX()*p1p2.getY() - p3p4.getY()*p1p2.getX() ) != 0 && p1p2.getX() != 0 ) //avoid division through zero
292 {
293 t2 = ( p1->x() * p1p2.getY() - p1->y() * p1p2.getX() + p3->y() * p1p2.getX() - p3->x() * p1p2.getY() ) / ( p3p4.getX() * p1p2.getY() - p3p4.getY() * p1p2.getX() );
294 t1 = ( p3->x() - p1->x() + t2 * p3p4.getX() ) / p1p2.getX();
295 }
296 else if ( ( p1p2.getX()*p3p4.getY() - p1p2.getY()*p3p4.getX() ) != 0 && p3p4.getX() != 0 )
297 {
298 t1 = ( p3->x() * p3p4.getY() - p3->y() * p3p4.getX() - p1->x() * p3p4.getY() + p1->y() * p3p4.getX() ) / ( p1p2.getX() * p3p4.getY() - p1p2.getY() * p3p4.getX() );
299 t2 = ( p1->x() + t1 * p1p2.getX() - p3->x() ) / p3p4.getX();
300 }
301 else//the lines are parallel
302 {
303 intersection_point->setX( 0 );
304 intersection_point->setY( 0 );
305 intersection_point->setZ( 0 );
306 return false;
307 }
308
309 if ( t1 > 0 && t1 < 1 && t2 > 0 && t2 < 1 )
310 {
311 if ( ( *p1 ) == ( *p3 ) || ( *p1 ) == ( *p4 ) || ( *p2 ) == ( *p3 ) || ( *p2 ) == ( *p4 ) ) //the lines touch each other, so they do not cross
312 {
313 intersection_point->setX( 0 );
314 intersection_point->setY( 0 );
315 intersection_point->setZ( 0 );
316 return false;
317 }
318 //calculate the intersection point
319 intersection_point->setX( p1->x() * ( 1 - t1 ) + p2->x()*t1 );
320 intersection_point->setY( p1->y() * ( 1 - t1 ) + p2->y()*t1 );
321 intersection_point->setZ( 0 );
322 return true;
323 }
324 else
325 {
326 return false;
327 }
328 }
329
330 else
331 {
332 QgsDebugError( QStringLiteral( "warning, null pointer" ) );
333 return false;
334 }
335}
336
337int MathUtils::lower( int n, int i )
338{
339 if ( i >= 0 && i <= n )
340 {
341 return faculty( n ) / ( faculty( i ) * faculty( n - i ) );
342 }
343 else
344 {
345 return 0;
346 }
347}
348
350{
351 if ( pa && pb && pc )
352 {
353 const double deter = ( pa->x() * pb->y() + pb->x() * pc->y() + pc->x() * pa->y() - pa->x() * pc->y() - pb->x() * pa->y() - pc->x() * pb->y() );
354 return 0.5 * deter;
355 }
356 else//null pointer
357 {
358 QgsDebugError( QStringLiteral( "warning, null pointer" ) );
359 return 0;
360 }
361}
362
363double MathUtils::calcCubicHermitePoly( int n, int i, double t )
364{
365 if ( n != 3 || i > n )
366 {
367 QgsDebugError( QStringLiteral( "error, can't calculate hermite polynom" ) );
368 }
369
370 if ( n == 3 && i == 0 )
371 {
372 return ( calcBernsteinPoly( 3, 0, t ) + calcBernsteinPoly( 3, 1, t ) );
373 }
374
375 if ( n == 3 && i == 1 )
376 {
377 return ( calcBernsteinPoly( 3, 1, t ) * 0.33333333 );
378 }
379
380 if ( n == 3 && i == 2 )
381 {
382 return ( calcBernsteinPoly( 3, 2, t ) * -0.33333333 );
383 }
384
385 if ( n == 3 && i == 3 )
386 {
387 return ( calcBernsteinPoly( 3, 2, t ) + calcBernsteinPoly( 3, 3, t ) );
388 }
389 else//something went wrong
390 {
391 QgsDebugError( QStringLiteral( "unexpected error" ) );
392 return 0;
393 }
394}
395
396double MathUtils::cFDerCubicHermitePoly( int n, int i, double t )
397{
398 if ( n != 3 || i > n )
399 {
400 QgsDebugError( QStringLiteral( "error, can't calculate hermite polynom" ) );
401 }
402
403 if ( n == 3 && i == 0 )
404 {
405 //cout << "cFDerBernsteinPoly(3,0,t): " << cFDerBernsteinPoly(3,0,t) << endl << flush;
406 //cout << "cFDerBernsteinPoly(3,1,t): " << cFDerBernsteinPoly(3,1,t) << endl << flush;
407 return ( cFDerBernsteinPoly( 3, 0, t ) + cFDerBernsteinPoly( 3, 1, t ) );
408 }
409
410 if ( n == 3 && i == 1 )
411 {
412 //cout << "cFDerBernsteinPoly(3,1,t)/3: " << cFDerBernsteinPoly(3,1,t)/3 << endl << flush;
413 return ( cFDerBernsteinPoly( 3, 1, t ) / 3 );
414 }
415
416 if ( n == 3 && i == 2 )
417 {
418 //cout << "cFDerBernsteinPoly(3,2,t)/3: " << cFDerBernsteinPoly(3,2,t)/3 << endl << flush;
419 return ( -cFDerBernsteinPoly( 3, 2, t ) / 3 );
420 }
421
422 if ( n == 3 && i == 3 )
423 {
424 //cout << "cFDerBernsteinPoly(3,2,t): " << cFDerBernsteinPoly(3,2,t) << endl << flush;
425 //cout << "cFDerBernsteinPoly(3,3,t): " << cFDerBernsteinPoly(3,3,t) << endl << flush;
426 return ( cFDerBernsteinPoly( 3, 2, t ) + cFDerBernsteinPoly( 3, 3, t ) );
427 }
428 else
429 {
430 QgsDebugError( QStringLiteral( "unexpected error" ) );
431 return 0;
432 }
433}
434
435bool MathUtils::derVec( const Vector3D *v1, const Vector3D *v2, Vector3D *result, double x, double y )
436{
437 if ( v1 && v2 && result )//no null pointers
438 {
439 const double u = ( x * v2->getY() - y * v2->getX() ) / ( v1->getX() * v2->getY() - v1->getY() * v2->getX() );
440 const double v = ( x * v1->getY() - y * v1->getX() ) / ( v2->getX() * v1->getY() - v2->getY() * v1->getX() );
441 result->setX( x );
442 result->setY( y );
443 result->setZ( u * v1->getZ() + v * v2->getZ() );
444 return true;
445 }
446 return false;
447}
448
449
450bool MathUtils::normalLeft( Vector3D *v1, Vector3D *result, double length )
451{
452 if ( v1 && result )//we don't like null pointers
453 {
454 if ( v1->getY() == 0 )//this would cause a division with zero
455 {
456 result->setX( 0 );
457
458 if ( v1->getX() < 0 )
459 {
460 result->setY( -length );
461 }
462
463 else
464 {
465 result->setY( length );
466 }
467 return true;
468 }
469
470 //coefficients for the quadratic equation
471 const double a = 1 + ( v1->getX() * v1->getX() ) / ( v1->getY() * v1->getY() );
472 const double b = 0;
473 const double c = -( length * length );
474 const double d = b * b - 4 * a * c;
475
476 if ( d < 0 )//no solution in R
477 {
478 QgsDebugError( QStringLiteral( "Determinant Error" ) );
479 return false;
480 }
481
482 result->setX( ( -b + std::sqrt( d ) ) / ( 2 * a ) ); //take one of the two solutions of the quadratic equation
483 result->setY( ( -result->getX()*v1->getX() ) / v1->getY() );
484
485 const QgsPoint point1( 0, 0, 0 );
486 const QgsPoint point2( v1->getX(), v1->getY(), 0 );
487 const QgsPoint point3( result->getX(), result->getY(), 0 );
488
489 if ( !( leftOf( point1, &point2, &point3 ) < 0 ) )//if we took the solution on the right side, change the sign of the components of the result
490 {
491 result->setX( -result->getX() );
492 result->setY( -result->getY() );
493 }
494
495 return true;
496 }
497
498 else
499 {
500 return false;
501 }
502}
503
504
505
506bool MathUtils::normalRight( Vector3D *v1, Vector3D *result, double length )
507{
508 if ( v1 && result )//we don't like null pointers
509 {
510
511 if ( v1->getY() == 0 )//this would cause a division with zero
512 {
513 result->setX( 0 );
514
515 if ( v1->getX() < 0 )
516 {
517 result->setY( length );
518 }
519
520 else
521 {
522 result->setY( -length );
523 }
524 return true;
525 }
526
527 //coefficients for the quadratic equation
528 const double a = 1 + ( v1->getX() * v1->getX() ) / ( v1->getY() * v1->getY() );
529 const double b = 0;
530 const double c = -( length * length );
531 const double d = b * b - 4 * a * c;
532
533 if ( d < 0 )//no solution in R
534 {
535 QgsDebugError( QStringLiteral( "Determinant Error" ) );
536 return false;
537 }
538
539 result->setX( ( -b + std::sqrt( d ) ) / ( 2 * a ) ); //take one of the two solutions of the quadratic equation
540 result->setY( ( -result->getX()*v1->getX() ) / v1->getY() );
541
542 const QgsPoint point1( 0, 0, 0 );
543 const QgsPoint point2( v1->getX(), v1->getY(), 0 );
544 const QgsPoint point3( result->getX(), result->getY(), 0 );
545
546 if ( ( leftOf( point1, &point2, &point3 ) < 0 ) ) //if we took the solution on the right side, change the sign of the components of the result
547 {
548 result->setX( -result->getX() );
549 result->setY( -result->getY() );
550 }
551
552 return true;
553 }
554
555 else
556 {
557 return false;
558 }
559}
560
561
563{
564 if ( p1 && p2 && p3 && vec )//no null pointers
565 {
566 const double ax = p2->x() - p1->x();
567 const double ay = p2->y() - p1->y();
568 const double az = p2->z() - p1->z();
569 const double bx = p3->x() - p1->x();
570 const double by = p3->y() - p1->y();
571 const double bz = p3->z() - p1->z();
572
573 vec->setX( ay * bz - az * by );
574 vec->setY( az * bx - ax * bz );
575 vec->setZ( ax * by - ay * bx );
576 }
577
578}
579
580double MathUtils::crossVec( QgsPoint *first, Vector3D *vec1, QgsPoint *second, Vector3D *vec2 )
581{
582 if ( first && vec1 && second && vec2 )
583 {
584 if ( ( vec2->getX()*vec1->getY() - vec2->getY()*vec1->getX() ) != 0 )
585 {
586 /*cout << "first: " << first->getX() << "//" << first->getY() << "//" << first->getZ() << endl << flush;
587 cout << "vec1: " << vec1->getX() << "//" << vec1->getY() << "//" << vec1->getZ() << endl << flush;
588 cout << "second: " << second->getX() << "//" << second->getY() << "//" << second->getZ() << endl << flush;
589 cout << "vec2: " << vec2->getX() << "//" << vec2->getY() << "//" << vec2->getZ() << endl << flush;
590 cout << "t2: " << ((first->getX()*vec1->getY()-first->getY()*vec1->getX()-second->getX()*vec1->getY()+second->getY()*vec1->getX())/(vec2->getX()*vec1->getY()-vec2->getY()*vec1->getX())) << endl << flush;*/
591
592 return ( ( first->x() * vec1->getY() - first->y() * vec1->getX() - second->x() * vec1->getY() + second->y() * vec1->getX() ) / ( vec2->getX() * vec1->getY() - vec2->getY() * vec1->getX() ) );
593
594 }
595 else//if a division by zero would occur
596 {
597 QgsDebugError( QStringLiteral( "warning: vectors are parallel" ) );
598 return 0;
599 }
600 }
601
602
603 else//null pointer
604 {
605 QgsDebugError( QStringLiteral( "warning, null pointer" ) );
606 return 0;
607 }
608}
609
610bool MathUtils::pointInsideTriangle( double x, double y, QgsPoint *p1, QgsPoint *p2, QgsPoint *p3 )
611{
612 const QgsPoint thepoint( x, y, 0 );
613 if ( MathUtils::leftOf( thepoint, p1, p2 ) > 0 )
614 {
615 return false;
616 }
617 if ( MathUtils::leftOf( thepoint, p2, p3 ) > 0 )
618 {
619 return false;
620 }
621 if ( MathUtils::leftOf( thepoint, p3, p1 ) > 0 )
622 {
623 return false;
624 }
625 return true;
626}
627
628bool MathUtils::normalMinDistance( Vector3D *tangent, Vector3D *target, Vector3D *result )
629{
630 if ( tangent && target && result )
631 {
632 const double xt = tangent->getX();
633 const double yt = tangent->getY();
634 const double zt = tangent->getZ();
635
636 const double xw = target->getX();
637 const double yw = target->getY();
638 const double zw = target->getZ();
639
640 double xg1, yg1, zg1;//the coordinates of the first result
641 double xg2, yg2, zg2;//the coordinates of the second result
642
643 //calculate xg
644 const double xgalpha1 = 1 / ( 2 * xt * xt * yw * yw * zt * zt - 2 * zt * zt * zt * xt * zw * xw + yt * yt * yt * yt * zw * zw + yt * yt * zw * zw * zt * zt + xt * xt * yt * yt * xw * xw + xt * xt * yw * yw * yt * yt - 2 * xt * xt * xt * zt * zw * xw + yt * yt * yt * yt * xw * xw + yt * yt * yw * yw * zt * zt + 2 * xt * xt * yt * yt * zw * zw - 2 * yt * yt * yt * yw * zt * zw + zt * zt * xt * xt * zw * zw + zt * zt * zt * zt * xw * xw + xt * xt * zt * zt * xw * xw + 2 * zt * zt * xw * xw * yt * yt - 2 * xt * xt * yw * zt * yt * zw - 2 * xt * yt * yt * yt * xw * yw - 2 * xt * xt * xt * yw * yt * xw - 2 * xt * zt * zt * xw * yt * yw - 2 * xt * zt * xw * yt * yt * zw - 2 * yw * zt * zt * zt * yt * zw + xt * xt * xt * xt * yw * yw + yw * yw * zt * zt * zt * zt + xt * xt * xt * xt * zw * zw );
645 if ( xgalpha1 < 0 )
646 {
647 QgsDebugError( QStringLiteral( "warning, only complex solution of xg" ) );
648 return false;
649 }
650 xg1 = std::sqrt( xgalpha1 ) * ( -yt * yw * xt + yt * yt * xw + xw * zt * zt - zt * xt * zw );
651 xg2 = -sqrt( xgalpha1 ) * ( -yt * yw * xt + yt * yt * xw + xw * zt * zt - zt * xt * zw );
652
653 //calculate yg
654 const double ygalpha1 = 1 / ( 2 * xt * xt * yw * yw * zt * zt - 2 * zt * zt * zt * xt * zw * xw + yt * yt * yt * yt * zw * zw + yt * yt * zw * zw * zt * zt + xt * xt * yt * yt * xw * xw + xt * xt * yw * yw * yt * yt - 2 * xt * xt * xt * zt * zw * xw + yt * yt * yt * yt * xw * xw + yt * yt * yw * yw * zt * zt + 2 * xt * xt * yt * yt * zw * zw - 2 * yt * yt * yt * yw * zt * zw + zt * zt * xt * xt * zw * zw + zt * zt * zt * zt * xw * xw + xt * xt * zt * zt * xw * xw + 2 * zt * zt * xw * xw * yt * yt - 2 * xt * xt * yw * zt * yt * zw - 2 * xt * yt * yt * yt * xw * yw - 2 * xt * xt * xt * yw * yt * xw - 2 * xt * zt * zt * xw * yt * yw - 2 * xt * zt * xw * yt * yt * zw - 2 * yw * zt * zt * zt * yt * zw + xt * xt * xt * xt * yw * yw + yw * yw * zt * zt * zt * zt + xt * xt * xt * xt * zw * zw );
655 if ( ygalpha1 < 0 )
656 {
657 QgsDebugError( QStringLiteral( "warning, only complex solution of yg" ) );
658 return false;
659 }
660 yg1 = -sqrt( ygalpha1 ) * ( -yw * xt * xt - zt * zt * yw + zt * yt * zw + yt * xw * xt );
661 yg2 = std::sqrt( ygalpha1 ) * ( -yw * xt * xt - zt * zt * yw + zt * yt * zw + yt * xw * xt );
662
663 //calculate zg
664 const double zgalpha1 = 1 / ( 2 * xt * xt * yw * yw * zt * zt - 2 * zt * zt * zt * xt * zw * xw + yt * yt * yt * yt * zw * zw + yt * yt * zw * zw * zt * zt + xt * xt * yt * yt * xw * xw + xt * xt * yw * yw * yt * yt - 2 * xt * xt * xt * zt * zw * xw + yt * yt * yt * yt * xw * xw + yt * yt * yw * yw * zt * zt + 2 * xt * xt * yt * yt * zw * zw - 2 * yt * yt * yt * yw * zt * zw + zt * zt * xt * xt * zw * zw + zt * zt * zt * zt * xw * xw + xt * xt * zt * zt * xw * xw + 2 * zt * zt * xw * xw * yt * yt - 2 * xt * xt * yw * zt * yt * zw - 2 * xt * yt * yt * yt * xw * yw - 2 * xt * xt * xt * yw * yt * xw - 2 * xt * zt * zt * xw * yt * yw - 2 * xt * zt * xw * yt * yt * zw - 2 * yw * zt * zt * zt * yt * zw + xt * xt * xt * xt * yw * yw + yw * yw * zt * zt * zt * zt + xt * xt * xt * xt * zw * zw );
665 if ( zgalpha1 < 0 )
666 {
667 QgsDebugError( QStringLiteral( "warning, only complex solution of zg" ) );
668 return false;
669 }
670 zg1 = -sqrt( zgalpha1 ) * ( yt * yw * zt - yt * yt * zw + xw * zt * xt - xt * xt * zw );
671 zg2 = std::sqrt( zgalpha1 ) * ( yt * yw * zt - yt * yt * zw + xw * zt * xt - xt * xt * zw );
672
673 const double distance1 = QgsGeometryUtilsBase::distance2D( xw, yw, xg1, yg1 );
674 const double distance2 = QgsGeometryUtilsBase::distance2D( xw, yw, xg2, yg2 );
675
676 if ( distance1 <= distance2 )//find out, which solution is the maximum and which the minimum
677 {
678 result->setX( xg1 );
679 result->setY( yg1 );
680 result->setZ( zg1 );
681 }
682 else
683 {
684 result->setX( xg2 );
685 result->setY( yg2 );
686 result->setZ( zg2 );
687 }
688 return true;
689 }
690
691 else
692 {
693 QgsDebugError( QStringLiteral( "warning, null pointer" ) );
694 return false;
695 }
696}
697
698
699double MathUtils::planeTest( QgsPoint *test, QgsPoint *pt1, QgsPoint *pt2, QgsPoint *pt3 )
700{
701 if ( test && pt1 && pt2 && pt3 )
702 {
703 const double a = ( pt1->z() * ( pt2->y() - pt3->y() ) + pt2->z() * ( pt3->y() - pt1->y() ) + pt3->z() * ( pt1->y() - pt2->y() ) ) / ( ( pt1->x() - pt2->x() ) * ( pt2->y() - pt3->y() ) - ( pt2->x() - pt3->x() ) * ( pt1->y() - pt2->y() ) );
704 const double b = ( pt1->z() * ( pt2->x() - pt3->x() ) + pt2->z() * ( pt3->x() - pt1->x() ) + pt3->z() * ( pt1->x() - pt2->x() ) ) / ( ( pt1->y() - pt2->y() ) * ( pt2->x() - pt3->x() ) - ( pt2->y() - pt3->y() ) * ( pt1->x() - pt2->x() ) );
705 const double c = pt1->z() - a * pt1->x() - b * pt1->y();
706 const double zpredicted = test->x() * a + test->y() * b + c;
707 return ( test->z() - zpredicted );
708 }
709 else
710 {
711 QgsDebugError( QStringLiteral( "warning, null pointer" ) );
712 return 0;
713 }
714}
715
717{
718 if ( p1 && p2 && p3 && p4 )
719 {
720 const Vector3D v1( p2->x() - p1->x(), p2->y() - p1->y(), 0 );
721 const Vector3D v2( p4->x() - p3->x(), p4->y() - p3->y(), 0 );
722 const double value = acos( ( v1.getX() * v2.getX() + v1.getY() * v2.getY() ) / ( v1.getLength() * v2.getLength() ) ) * 180 / M_PI;
723 return value;
724 }
725 else
726 {
727 QgsDebugError( QStringLiteral( "warning, null pointer" ) );
728 return 0;
729 }
730}
static double distance2D(double x1, double y1, double x2, double y2)
Returns the 2D distance between (x1, y1) and (x2, y2).
Point geometry type, with support for z-dimension and m-values.
Definition qgspoint.h:49
void setY(double y)
Sets the point's y-coordinate.
Definition qgspoint.h:343
void setX(double x)
Sets the point's x-coordinate.
Definition qgspoint.h:332
double z
Definition qgspoint.h:54
double x
Definition qgspoint.h:52
double distance(double x, double y) const
Returns the Cartesian 2D distance between this point and a specified x, y coordinate.
Definition qgspoint.h:393
void setZ(double z)
Sets the point's z-coordinate.
Definition qgspoint.h:356
double y
Definition qgspoint.h:53
Class Vector3D represents a 3D-Vector, capable to store x-,y- and z-coordinates in double values.
Definition Vector3D.h:36
void setX(double x)
Sets the x-component of the vector.
Definition Vector3D.h:106
double getY() const
Returns the y-component of the vector.
Definition Vector3D.h:96
double getX() const
Returns the x-component of the vector.
Definition Vector3D.h:91
void setY(double y)
Sets the y-component of the vector.
Definition Vector3D.h:111
double getZ() const
Returns the z-component of the vector.
Definition Vector3D.h:101
double getLength() const
Returns the length of the vector.
Definition Vector3D.cpp:19
void setZ(double z)
Sets the z-component of the vector.
Definition Vector3D.h:116
bool ANALYSIS_EXPORT normalLeft(Vector3D *v1, Vector3D *result, double length)
Assigns the vector 'result', which is normal to the vector 'v1', on the left side of v1 and has lengt...
bool ANALYSIS_EXPORT inDiametral(QgsPoint *p1, QgsPoint *p2, QgsPoint *point)
Tests, whether 'point' is inside the diametral circle through 'p1' and 'p2'.
bool ANALYSIS_EXPORT derVec(const Vector3D *v1, const Vector3D *v2, Vector3D *result, double x, double y)
Calculates the z-component of a vector with coordinates 'x' and 'y'which is in the same tangent plane...
double ANALYSIS_EXPORT calcBernsteinPoly(int n, int i, double t)
Calculates the value of a Bernstein polynomial.
int ANALYSIS_EXPORT faculty(int n)
Faculty function.
bool ANALYSIS_EXPORT BarycentricToXY(double u, double v, double w, QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, QgsPoint *result)
Definition MathUtils.cpp:51
double ANALYSIS_EXPORT angle(QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, QgsPoint *p4)
Calculates the angle between two segments (in 2 dimension, z-values are ignored)
double ANALYSIS_EXPORT triArea(QgsPoint *pa, QgsPoint *pb, QgsPoint *pc)
Returns the area of a triangle. If the points are ordered counterclockwise, the value will be positiv...
bool ANALYSIS_EXPORT pointInsideTriangle(double x, double y, QgsPoint *p1, QgsPoint *p2, QgsPoint *p3)
Returns true, if the point with coordinates x and y is inside (or at the edge) of the triangle p1,...
bool ANALYSIS_EXPORT normalRight(Vector3D *v1, Vector3D *result, double length)
Assigns the vector 'result', which is normal to the vector 'v1', on the right side of v1 and has leng...
double ANALYSIS_EXPORT crossVec(QgsPoint *first, Vector3D *vec1, QgsPoint *second, Vector3D *vec2)
Calculates the intersection of the two vectors vec1 and vec2, which start at first(vec1) and second(v...
bool ANALYSIS_EXPORT lineIntersection(QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, QgsPoint *p4)
Returns true, if line1 (p1 to p2) and line2 (p3 to p4) intersect. If the lines have an endpoint in co...
int ANALYSIS_EXPORT lower(int n, int i)
Lower function.
double ANALYSIS_EXPORT planeTest(QgsPoint *test, QgsPoint *pt1, QgsPoint *pt2, QgsPoint *pt3)
Tests, if 'test' is in the same plane as 'p1', 'p2' and 'p3' and returns the z-difference from the pl...
bool ANALYSIS_EXPORT circumcenter(QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, QgsPoint *result)
Calculates the center of the circle passing through p1, p2 and p3. Returns true in case of success an...
double ANALYSIS_EXPORT cFDerBernsteinPoly(int n, int i, double t)
Calculates the first derivative of a Bernstein polynomial with respect to the parameter t.
double ANALYSIS_EXPORT cFDerCubicHermitePoly(int n, int i, double t)
Calculates the first derivative of a cubic Hermite polynomial with respect to the parameter t.
double ANALYSIS_EXPORT leftOf(const QgsPoint &thepoint, const QgsPoint *p1, const QgsPoint *p2)
Returns whether 'thepoint' is left or right of the line from 'p1' to 'p2'. Negative values mean left ...
double ANALYSIS_EXPORT distPointFromLine(QgsPoint *thepoint, QgsPoint *p1, QgsPoint *p2)
Calculates the (2 dimensional) distance from 'thepoint' to the line defined by p1 and p2.
bool ANALYSIS_EXPORT calcBarycentricCoordinates(double x, double y, QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, QgsPoint *result)
Calculates the barycentric coordinates of a point (x,y) with respect to p1, p2, p3 and stores the thr...
Definition MathUtils.cpp:23
double ANALYSIS_EXPORT calcCubicHermitePoly(int n, int i, double t)
Calculates the value of a cubic Hermite polynomial.
bool ANALYSIS_EXPORT inCircle(QgsPoint *testp, QgsPoint *p1, QgsPoint *p2, QgsPoint *p3)
Tests, whether 'testp' is inside the circle through 'p1', 'p2' and 'p3'.
void ANALYSIS_EXPORT normalFromPoints(QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, Vector3D *vec)
Calculates the normal vector of the plane through the points p1, p2 and p3 and assigns the result to ...
bool ANALYSIS_EXPORT normalMinDistance(Vector3D *tangent, Vector3D *target, Vector3D *result)
Calculates a Vector orthogonal to 'tangent' with length 1 and closest possible to result....
As part of the API refactoring and improvements which landed in the Processing API was substantially reworked from the x version This was done in order to allow much of the underlying Processing framework to be ported into c
#define QgsDebugError(str)
Definition qgslogger.h:38