QGIS API Documentation  3.0.2-Girona (307d082)
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 *
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 
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( "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( "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( "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( "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( "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 
201 double MathUtils::distPointFromLine( QgsPoint *thepoint, QgsPoint *p1, QgsPoint *p2 )
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( "warning, null pointer" );
217  return 0;
218  }
219 }
220 
221 int MathUtils::faculty( int n )
222 {
223  return std::tgamma( n + 1 );
224 }
225 
226 bool MathUtils::inCircle( QgsPoint *testp, QgsPoint *p1, QgsPoint *p2, QgsPoint *p3 )
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( "warning, null pointer" );
262  return false;
263  }
264 }
265 
266 bool MathUtils::inDiametral( QgsPoint *p1, QgsPoint *p2, QgsPoint *point )
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( "warning, null pointer" );
307  return 0;
308  }
309 }
310 
311 bool MathUtils::lineIntersection( QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, QgsPoint *p4 )
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( "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( "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 
419 double MathUtils::triArea( QgsPoint *pa, QgsPoint *pb, QgsPoint *pc )
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( "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( "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( "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( "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( "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( "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( "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 
632 void MathUtils::normalFromPoints( QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, Vector3D *vec )
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( "warning: vectors are parallel" );
668  return 0;
669  }
670  }
671 
672 
673  else//null pointer
674  {
675  QgsDebugMsg( "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( "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( "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( "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( "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( "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( "warning, null pointer" );
798  return 0;
799  }
800 }
Vector3D::setX
void setX(double x)
Sets the x-component of the vector.
Definition: Vector3D.h:103
QgsPoint::y
double y
Definition: qgspoint.h:42
MathUtils::normalRight
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
QgsPoint::setZ
void setZ(double z)
Sets the point&#39;s z-coordinate.
Definition: qgspoint.h:216
MathUtils::cFDerBernsteinPoly
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
QgsDebugMsg
#define QgsDebugMsg(str)
Definition: qgslogger.h:38
MathUtils::lineIntersection
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
MathUtils::triArea
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
MathUtils::BarycentricToXY
bool ANALYSIS_EXPORT BarycentricToXY(double u, double v, double w, QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, QgsPoint *result)
Definition: MathUtils.cpp:50
MathUtils::calcBarycentricCoordinates
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
MathUtils::inDiametral
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
MathUtils::angle
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
MathUtils::pointInsideTriangle
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
Vector3D::getZ
double getZ() const
Returns the z-component of the vector.
Definition: Vector3D.h:98
MathUtils::normalMinDistance
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
Vector3D
Class Vector3D represents a 3D-Vector, capable to store x-,y- and z-coordinates in double values...
Definition: Vector3D.h:33
Vector3D::setZ
void setZ(double z)
Sets the z-component of the vector.
Definition: Vector3D.h:113
MathUtils::normalLeft
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
MathUtils::circumcenter
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
MathUtils::distPointFromLine
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
QgsPoint
Point geometry type, with support for z-dimension and m-values.
Definition: qgspoint.h:37
MathUtils::planeTest
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
QgsPoint::setX
void setX(double x)
Sets the point&#39;s x-coordinate.
Definition: qgspoint.h:192
QgsPoint::setY
void setY(double y)
Sets the point&#39;s y-coordinate.
Definition: qgspoint.h:203
Vector3D::getY
double getY() const
Returns the y-component of the vector.
Definition: Vector3D.h:93
MathUtils::inCircle
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
Vector3D::setY
void setY(double y)
Sets the y-component of the vector.
Definition: Vector3D.h:108
MathUtils::crossVec
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
MathUtils::calcBernsteinPoly
double ANALYSIS_EXPORT calcBernsteinPoly(int n, int i, double t)
Calculates the value of a Bernstein polynomial.
Definition: MathUtils.cpp:101
QgsPoint::z
double z
Definition: qgspoint.h:43
Vector3D::getX
double getX() const
Returns the x-component of the vector.
Definition: Vector3D.h:88
MathUtils::normalFromPoints
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
MathUtils::lower
int ANALYSIS_EXPORT lower(int n, int i)
Lower function.
Definition: MathUtils.cpp:407
MathUtils::cFDerCubicHermitePoly
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
MathUtils::leftOf
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
MathUtils::calcCubicHermitePoly
double ANALYSIS_EXPORT calcCubicHermitePoly(int n, int i, double t)
Calculates the value of a cubic Hermite polynomial.
Definition: MathUtils.cpp:433
MathUtils::derVec
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
MathUtils::faculty
int ANALYSIS_EXPORT faculty(int n)
Faculty function.
Definition: MathUtils.cpp:221
QgsPoint::x
double x
Definition: qgspoint.h:41
Generated on Sat Apr 21 2018 11:45:00 for QGIS API Documentation by   doxygen 1.8.13