MathUtils.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  MathUtils.cpp
3  -------------
4  copyright : (C) 2004 by Marco Hugentobler
5  email : [email protected]
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"
21
22 bool 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  double area = triArea( p1, p2, p3 );
28  if ( area == 0 )//p1, p2, p3 are in a line
29  {
30  return false;
31  }
32  double area1 = triArea( &p, p2, p3 );
33  double area2 = triArea( p1, &p, p3 );
34  double area3 = triArea( p1, p2, &p );
35  double u = area1 / area;
36  double v = area2 / area;
37  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
50 bool 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  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  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//dorp 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
101 double 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
111 double 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
116 bool 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  double distp1p2 = std::sqrt( ( p1->x() - p2->x() ) * ( p1->x() - p2->x() ) + ( p1->y() - p2->y() ) * ( p1->y() - p2->y() ) );
121  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  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
164 bool 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  double a = normal.getX();
209  double b = normal.getY();
210  double c = -( normal.getX() * p2->x() + normal.getY() * p2->y() );
211  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
221 int MathUtils::faculty( int n )
222 {
223  return std::tgamma( n + 1 );
224 }
225
227 {
228  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  double xmin = std::min( std::min( ax, px ), std::min( bx, cx ) );
242  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
272 bool 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
292 double MathUtils::leftOf( const QgsPoint &thepoint, const QgsPoint *p1, const QgsPoint *p2 )
293 {
294  if ( p1 && p2 )
295  {
296  //just for debugging
297  double f1 = thepoint.x() - p1->x();
298  double f2 = p2->y() - p1->y();
299  double f3 = thepoint.y() - p1->y();
300  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  Vector3D p1p2( p2->x() - p1->x(), p2->y() - p1->y(), 0 );
317  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
354 bool 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  Vector3D p1p2( p2->x() - p1->x(), p2->y() - p1->y(), 0 );
360  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
407 int 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  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
433 double 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
466 double 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
505 bool 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  double u = ( x * v2->getY() - y * v2->getX() ) / ( v1->getX() * v2->getY() - v1->getY() * v2->getX() );
510  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
520 bool 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  double a = 1 + ( v1->getX() * v1->getX() ) / ( v1->getY() * v1->getY() );
542  double b = 0;
543  double c = -( length * length );
544  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  QgsPoint point1( 0, 0, 0 );
556  QgsPoint point2( v1->getX(), v1->getY(), 0 );
557  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
576 bool 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  double a = 1 + ( v1->getX() * v1->getX() ) / ( v1->getY() * v1->getY() );
599  double b = 0;
600  double c = -( length * length );
601  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  QgsPoint point1( 0, 0, 0 );
613  QgsPoint point2( v1->getX(), v1->getY(), 0 );
614  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  double ax = p2->x() - p1->x();
637  double ay = p2->y() - p1->y();
638  double az = p2->z() - p1->z();
639  double bx = p3->x() - p1->x();
640  double by = p3->y() - p1->y();
641  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
650 double 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
680 bool MathUtils::pointInsideTriangle( double x, double y, QgsPoint *p1, QgsPoint *p2, QgsPoint *p3 )
681 {
682  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
698 bool MathUtils::normalMinDistance( Vector3D *tangent, Vector3D *target, Vector3D *result )
699 {
700  if ( tangent && target && result )
701  {
702  double xt = tangent->getX();
703  double yt = tangent->getY();
704  double zt = tangent->getZ();
705
706  double xw = target->getX();
707  double yw = target->getY();
708  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  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  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  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  double distance1 = std::sqrt( ( xw - xg1 ) * ( xw - xg1 ) + ( yw - yg1 ) * ( yw - yg1 ) + ( zw - zg1 ) * ( zw - zg1 ) );
744  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
769 double MathUtils::planeTest( QgsPoint *test, QgsPoint *pt1, QgsPoint *pt2, QgsPoint *pt3 )
770 {
771  if ( test && pt1 && pt2 && pt3 )
772  {
773  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  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  double c = pt1->z() - a * pt1->x() - b * pt1->y();
776  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
786 double MathUtils::angle( QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, QgsPoint *p4 )
787 {
788  if ( p1 && p2 && p3 && p4 )
789  {
790  Vector3D v1( p2->x() - p1->x(), p2->y() - p1->y(), 0 );
791  Vector3D v2( p4->x() - p3->x(), p4->y() - p3->y(), 0 );
792  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 }
void setX(double x)
Sets the x-component of the vector.
Definition: Vector3D.h:103
double y
Definition: qgspoint.h:42
bool ANALYSIS_EXPORT normalRight(Vector3D *v1, Vector3D *result, double length)
Assigns the vector &#39;result&#39;, which is normal to the vector &#39;v1&#39;, on the right side of v1 and has leng...
Definition: MathUtils.cpp:576
void setZ(double z)
Sets the point&#39;s z-coordinate.
Definition: qgspoint.h:237
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
#define QgsDebugMsg(str)
Definition: qgslogger.h:38
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
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 BarycentricToXY(double u, double v, double w, QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, QgsPoint *result)
Definition: MathUtils.cpp:50
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
bool ANALYSIS_EXPORT inDiametral(QgsPoint *p1, QgsPoint *p2, QgsPoint *point)
Tests, whether &#39;point&#39; is inside the diametral circle through &#39;p1&#39; and &#39;p2&#39;.
Definition: MathUtils.cpp:266
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
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
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
bool ANALYSIS_EXPORT normalMinDistance(Vector3D *tangent, Vector3D *target, Vector3D *result)
Calculates a Vector orthogonal to &#39;tangent&#39; with length 1 and closest possible to result...
Definition: MathUtils.cpp:698
Class Vector3D represents a 3D-Vector, capable to store x-,y- and z-coordinates in double values...
Definition: Vector3D.h:33
void setZ(double z)
Sets the z-component of the vector.
Definition: Vector3D.h:113
bool ANALYSIS_EXPORT normalLeft(Vector3D *v1, Vector3D *result, double length)
Assigns the vector &#39;result&#39;, which is normal to the vector &#39;v1&#39;, on the left side of v1 and has lengt...
Definition: MathUtils.cpp:520
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 distPointFromLine(QgsPoint *thepoint, QgsPoint *p1, QgsPoint *p2)
Calculates the (2 dimensional) distance from &#39;thepoint&#39; to the line defined by p1 and p2...
Definition: MathUtils.cpp:201
Point geometry type, with support for z-dimension and m-values.
Definition: qgspoint.h:37
double ANALYSIS_EXPORT planeTest(QgsPoint *test, QgsPoint *pt1, QgsPoint *pt2, QgsPoint *pt3)
Tests, if &#39;test&#39; is in the same plane as &#39;p1&#39;, &#39;p2&#39; and &#39;p3&#39; and returns the z-difference from the pl...
Definition: MathUtils.cpp:769
double getY() const
Returns the y-component of the vector.
Definition: Vector3D.h:93
double getX() const
Returns the x-component of the vector.
Definition: Vector3D.h:88
void setX(double x)
Sets the point&#39;s x-coordinate.
Definition: qgspoint.h:213
void setY(double y)
Sets the point&#39;s y-coordinate.
Definition: qgspoint.h:224
bool ANALYSIS_EXPORT inCircle(QgsPoint *testp, QgsPoint *p1, QgsPoint *p2, QgsPoint *p3)
Tests, whether &#39;testp&#39; is inside the circle through &#39;p1&#39;, &#39;p2&#39; and &#39;p3&#39;.
Definition: MathUtils.cpp:226
void setY(double y)
Sets the y-component of the vector.
Definition: Vector3D.h:108
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
double ANALYSIS_EXPORT calcBernsteinPoly(int n, int i, double t)
Calculates the value of a Bernstein polynomial.
Definition: MathUtils.cpp:101
double z
Definition: qgspoint.h:43
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
double getZ() const
Returns the z-component of the vector.
Definition: Vector3D.h:98
int ANALYSIS_EXPORT lower(int n, int i)
Lower function.
Definition: MathUtils.cpp:407
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 &#39;thepoint&#39; is left or right of the line from &#39;p1&#39; to &#39;p2&#39;. Negativ values mean left a...
Definition: MathUtils.cpp:292
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 derVec(const Vector3D *v1, const Vector3D *v2, Vector3D *result, double x, double y)
Calculates the z-component of a vector with coordinates &#39;x&#39; and &#39;y&#39;which is in the same tangent plane...
Definition: MathUtils.cpp:505
int ANALYSIS_EXPORT faculty(int n)
Faculty function.
Definition: MathUtils.cpp:221
double x
Definition: qgspoint.h:41