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