SAGA API  v9.5
geo_functions.cpp
Go to the documentation of this file.
1 
3 // //
4 // SAGA //
5 // //
6 // System for Automated Geoscientific Analyses //
7 // //
8 // Application Programming Interface //
9 // //
10 // Library: SAGA_API //
11 // //
12 //-------------------------------------------------------//
13 // //
14 // geo_functions.cpp //
15 // //
16 // Copyright (C) 2005 by Olaf Conrad //
17 // //
18 //-------------------------------------------------------//
19 // //
20 // This file is part of 'SAGA - System for Automated //
21 // Geoscientific Analyses'. //
22 // //
23 // This library is free software; you can redistribute //
24 // it and/or modify it under the terms of the GNU Lesser //
25 // General Public License as published by the Free //
26 // Software Foundation, either version 2.1 of the //
27 // License, or (at your option) any later version. //
28 // //
29 // This library is distributed in the hope that it will //
30 // be useful, but WITHOUT ANY WARRANTY; without even the //
31 // implied warranty of MERCHANTABILITY or FITNESS FOR A //
32 // PARTICULAR PURPOSE. See the GNU Lesser General Public //
33 // License for more details. //
34 // //
35 // You should have received a copy of the GNU Lesser //
36 // General Public License along with this program; if //
37 // not, see <http://www.gnu.org/licenses/>. //
38 // //
39 //-------------------------------------------------------//
40 // //
41 // contact: Olaf Conrad //
42 // Institute of Geography //
43 // University of Goettingen //
44 // Goldschmidtstr. 5 //
45 // 37077 Goettingen //
46 // Germany //
47 // //
48 // e-mail: oconrad@saga-gis.org //
49 // //
51 
52 //---------------------------------------------------------
53 #include "geo_tools.h"
54 #include "mat_tools.h"
55 
56 
58 // //
59 // //
60 // //
62 
63 //---------------------------------------------------------
64 bool SG_Is_Equal(double a, double b, double epsilon)
65 {
66  return( fabs(a - b) <= epsilon );
67 }
68 
69 //---------------------------------------------------------
70 bool SG_Is_Equal(const TSG_Point &A, const TSG_Point &B, double epsilon)
71 {
72  return( SG_Is_Equal(A.x, B.x, epsilon)
73  && SG_Is_Equal(A.y, B.y, epsilon) );
74 }
75 
76 //---------------------------------------------------------
77 bool SG_Is_Between(double x, double a, double b, double epsilon)
78 {
79  return( (a - epsilon <= x && x <= b + epsilon)
80  || (b - epsilon <= x && x <= a + epsilon) );
81 }
82 
83 bool SG_Is_Between(const TSG_Point &Point, const TSG_Point &Corner_A, const TSG_Point &Corner_B, double epsilon)
84 {
85  return( SG_Is_Between(Point.x, Corner_A.x, Corner_B.x, epsilon)
86  && SG_Is_Between(Point.y, Corner_A.y, Corner_B.y, epsilon) );
87 }
88 
89 
91 // //
92 // //
93 // //
95 
96 //---------------------------------------------------------
97 double SG_Get_Length(double dx, double dy)
98 {
99  return( sqrt(dx*dx + dy*dy) );
100 }
101 
102 //---------------------------------------------------------
103 double SG_Get_Distance(double ax, double ay, double bx, double by, bool bPolar)
104 {
105  return( bPolar ? SG_Get_Distance_Polar(ax, ay, bx, by) : SG_Get_Distance(ax, ay, bx, by) );
106 }
107 
108 //---------------------------------------------------------
109 double SG_Get_Distance(const TSG_Point &A, const TSG_Point &B, bool bPolar)
110 {
111  return( bPolar ? SG_Get_Distance_Polar(A, B) : SG_Get_Distance(A, B) );
112 }
113 
114 //---------------------------------------------------------
115 double SG_Get_Distance(double ax, double ay, double bx, double by)
116 {
117  ax -= bx;
118  ay -= by;
119 
120  return( sqrt(ax*ax + ay*ay) );
121 }
122 
123 //---------------------------------------------------------
124 double SG_Get_Distance(const TSG_Point &A, const TSG_Point &B)
125 {
126  double dx, dy;
127 
128  dx = B.x - A.x;
129  dy = B.y - A.y;
130 
131  return( sqrt(dx*dx + dy*dy) );
132 }
133 
134 //---------------------------------------------------------
135 double SG_Get_Distance(double ax, double ay, double az, double bx, double by, double bz)
136 {
137  ax -= bx;
138  ay -= by;
139  az -= bz;
140 
141  return( sqrt(ax*ax + ay*ay + az*az) );
142 }
143 
144 //---------------------------------------------------------
145 double SG_Get_Distance(const TSG_Point_3D &A, const TSG_Point_3D &B)
146 {
147  double dx, dy, dz;
148 
149  dx = B.x - A.x;
150  dy = B.y - A.y;
151  dz = B.z - A.z;
152 
153  return( sqrt(dx*dx + dy*dy + dz*dz) );
154 }
155 
156 //---------------------------------------------------------
157 double SG_Get_Distance_Polar(double aLon, double aLat, double bLon, double bLat, double a, double e, bool bDegree)
158 {
159  if( bDegree )
160  {
161  aLon *= M_DEG_TO_RAD;
162  aLat *= M_DEG_TO_RAD;
163  bLon *= M_DEG_TO_RAD;
164  bLat *= M_DEG_TO_RAD;
165  }
166 
167  if( e <= 0. )
168  {
169  return( a * acos(sin(aLat) * sin(bLat) + cos(aLat) * cos(bLat) * cos(bLon - aLon)) );
170  }
171  else
172  {
173  double F = (aLat + bLat) / 2.;
174  double G = (aLat - bLat) / 2.;
175  double l = (aLon - bLon) / 2.;
176 
177  double sin2_F = SG_Get_Square(sin(F));
178  double cos2_F = SG_Get_Square(cos(F));
179  double sin2_G = SG_Get_Square(sin(G));
180  double cos2_G = SG_Get_Square(cos(G));
181  double sin2_l = SG_Get_Square(sin(l));
182  double cos2_l = SG_Get_Square(cos(l));
183 
184  double S = sin2_G * cos2_l + cos2_F * sin2_l;
185  double C = cos2_G * cos2_l + sin2_F * sin2_l;
186 
187  double w = atan(sqrt(S / C));
188  double D = 2. * w * a;
189 
190  double R = sqrt(S * C) / w;
191  double H1 = (3. * R - 1.) / (2. * C);
192  double H2 = (3. * R + 1.) / (2. * S);
193 
194  double f = 1. / e;
195 
196  double d = D * (1. + f * H1 * sin2_F * cos2_G - f * H2 * cos2_F * sin2_G);
197 
198  return( d );
199  }
200 }
201 
202 //---------------------------------------------------------
203 double SG_Get_Distance_Polar(const TSG_Point &A, const TSG_Point &B, double a, double e, bool bDegree)
204 {
205  return( SG_Get_Distance_Polar(A.x, A.y, B.x, B.y, a, e, bDegree) );
206 }
207 
208 
210 // //
211 // //
212 // //
214 
215 //---------------------------------------------------------
216 double SG_Get_Angle_Of_Direction(double dx, double dy)
217 {
218  if( dx == 0. )
219  {
220  return( dy > 0. ? 0. : dy < 0. ? M_PI_180 : -M_PI_360 );
221  }
222 
223  dx = M_PI_090 - atan2(dy, dx);
224 
225  return( dx < 0. ? M_PI_360 + dx : dx );
226 }
227 
228 //---------------------------------------------------------
229 double SG_Get_Angle_Of_Direction(double ax, double ay, double bx, double by)
230 {
231  return( SG_Get_Angle_Of_Direction(bx - ax, by - ay) );
232 }
233 
234 //---------------------------------------------------------
235 double SG_Get_Angle_Of_Direction(const TSG_Point &A)
236 {
237  return( SG_Get_Angle_Of_Direction(A.x, A.y) );
238 }
239 
240 //---------------------------------------------------------
241 double SG_Get_Angle_Of_Direction(const TSG_Point &A, const TSG_Point &B)
242 {
243  return( SG_Get_Angle_Of_Direction(B.x - A.x, B.y - A.y) );
244 }
245 
246 //---------------------------------------------------------
247 double SG_Get_Angle_Difference(double a, double b)
248 {
249  double d = fmod(b - a, M_PI_360);
250 
251  if( d < 0. ) d += M_PI_360;
252 
253  return( d > M_PI_180 ? d - M_PI_180 : d );
254 }
255 
256 //---------------------------------------------------------
257 bool SG_is_Angle_Between(double Angle, double Angle_Min, double Angle_Max, bool bCheckRange)
258 {
259  if( bCheckRange )
260  {
261  Angle = fmod(Angle , M_PI_360); if( Angle < 0. ) Angle += M_PI_360;
262  Angle_Min = fmod(Angle_Min, M_PI_360); if( Angle_Min < 0. ) Angle_Min += M_PI_360;
263  Angle_Max = fmod(Angle_Max, M_PI_360); if( Angle_Max < 0. ) Angle_Max += M_PI_360;
264  }
265 
266  return( Angle_Min <= Angle_Max
267  ? Angle_Min <= Angle && Angle <= Angle_Max
268  : Angle_Min <= Angle || Angle <= Angle_Max
269  );
270 }
271 
272 
274 // //
275 // //
276 // //
278 
279 //---------------------------------------------------------
280 bool SG_Get_Crossing(TSG_Point &Crossing, const TSG_Point &a1, const TSG_Point &a2, const TSG_Point &b1, const TSG_Point &b2, bool bExactMatch)
281 {
282  //-----------------------------------------------------
283  if( bExactMatch
284  && ( (M_GET_MAX(a1.x, a2.x) < M_GET_MIN(b1.x, b2.x))
285  || (M_GET_MIN(a1.x, a2.x) > M_GET_MAX(b1.x, b2.x))
286  || (M_GET_MAX(a1.y, a2.y) < M_GET_MIN(b1.y, b2.y))
287  || (M_GET_MIN(a1.y, a2.y) > M_GET_MAX(b1.y, b2.y)) ) )
288  {
289  return( false );
290  }
291 
292  //-----------------------------------------------------
293  if( (a1.x == b1.x && a1.y == b1.y) || (a1.x == b2.x && a1.y == b2.y) )
294  {
295  Crossing = a1;
296 
297  return( true );
298  }
299 
300  if( (a2.x == b1.x && a2.y == b1.y) || (a2.x == b2.x && a2.y == b2.y) )
301  {
302  Crossing = a2;
303 
304  return( true );
305  }
306 
307  //-----------------------------------------------------
308  double lambda, div, a_dx, a_dy, b_dx, b_dy;
309 
310  a_dx = a2.x - a1.x;
311  a_dy = a2.y - a1.y;
312 
313  b_dx = b2.x - b1.x;
314  b_dy = b2.y - b1.y;
315 
316  if( (div = a_dx * b_dy - b_dx * a_dy) != 0. )
317  {
318  lambda = ((b1.x - a1.x) * b_dy - b_dx * (b1.y - a1.y)) / div;
319 
320  Crossing.x = a1.x + lambda * a_dx;
321  Crossing.y = a1.y + lambda * a_dy;
322 
323  if( !bExactMatch )
324  {
325  return( true );
326  }
327  else if( 0. <= lambda && lambda <= 1. )
328  {
329  lambda = ((b1.x - a1.x) * a_dy - a_dx * (b1.y - a1.y)) / div;
330 
331  if( 0. <= lambda && lambda <= 1. )
332  {
333  return( true );
334  }
335  }
336  }
337 
338  return( false );
339 }
340 
341 //---------------------------------------------------------
342 bool SG_Get_Crossing_InRegion(TSG_Point &Crossing, const TSG_Point &a, const TSG_Point &b, const TSG_Rect &Region)
343 {
344  TSG_Point ra, rb;
345 
346  //-----------------------------------------------------
347  ra.y = Region.yMin;
348  rb.y = Region.yMax;
349 
350  ra.x = rb.x = Region.xMin;
351 
352  if( SG_Get_Crossing(Crossing, a, b, ra, rb, true) )
353  {
354  return( true );
355  }
356 
357  //-----------------------------------------------------
358  ra.x = rb.x = Region.xMax;
359 
360  if( SG_Get_Crossing(Crossing, a, b, ra, rb, true) )
361  {
362  return( true );
363  }
364 
365  //-----------------------------------------------------
366  ra.x = Region.xMin;
367  ra.y = Region.yMax;
368 
369  if( SG_Get_Crossing(Crossing, a, b, ra, rb, true) )
370  {
371  return( true );
372  }
373 
374  //-----------------------------------------------------
375  ra.y = rb.y = Region.yMin;
376 
377  if( SG_Get_Crossing(Crossing, a, b, ra, rb, true) )
378  {
379  return( true );
380  }
381 
382  //-----------------------------------------------------
383  return( false );
384 }
385 
386 
388 // //
389 // //
390 // //
392 
393 //---------------------------------------------------------
394 bool SG_Is_Point_On_Line(const TSG_Point &Point, const TSG_Point &Line_A, const TSG_Point &Line_B, bool bExactMatch, double Epsilon)
395 {
396  if( SG_Is_Equal(Line_B.x, Line_A.x, Epsilon) ) // vertical line
397  {
398  return( SG_Is_Between(Point.y, Line_A.y, Line_B.y, Epsilon) && (!bExactMatch || SG_Is_Between(Point.x, Line_A.x, Line_B.x, Epsilon)) );
399  }
400 
401  if( bExactMatch && !SG_Is_Between(Point, Line_A, Line_B, Epsilon) )
402  {
403  return( false );
404  }
405 
406  double b = (Line_B.y - Line_A.y) / (Line_B.x - Line_A.x);
407  double a = Line_A.y - b * Line_A.x;
408 
409  return( SG_Is_Equal(Point.y, a + b * Point.x, Epsilon) );
410 }
411 
412 //---------------------------------------------------------
413 double SG_Get_Distance_To_Line(const TSG_Point &Point, const TSG_Point &Line_A, const TSG_Point &Line_B, bool bExactMatch)
414 {
415  TSG_Point Line_Point;
416 
417  return( SG_Get_Nearest_Point_On_Line(Point, Line_A, Line_B, Line_Point, bExactMatch) );
418 }
419 
420 //---------------------------------------------------------
421 double SG_Get_Nearest_Point_On_Line(const TSG_Point &Point, const TSG_Point &Line_A, const TSG_Point &Line_B, TSG_Point &Line_Point, bool bExactMatch)
422 {
423  CSG_Point Point_Ortho(
424  Point.x - (Line_B.y - Line_A.y),
425  Point.y + (Line_B.x - Line_A.x)
426  );
427 
428  if( !SG_Get_Crossing(Line_Point, Line_A, Line_B, Point, Point_Ortho, false) )
429  {
430  return( -1. );
431  }
432 
433  if( !bExactMatch || (bExactMatch
434  && SG_IS_BETWEEN(Line_A.x, Line_Point.x, Line_B.x)
435  && SG_IS_BETWEEN(Line_A.y, Line_Point.y, Line_B.y)) )
436  {
437  return( SG_Get_Distance(Point, Line_Point) );
438  }
439 
440  double Distance_A = SG_Get_Distance(Point, Line_A);
441  double Distance_B = SG_Get_Distance(Point, Line_B);
442 
443  if( Distance_A < Distance_B )
444  {
445  Line_Point = Line_A;
446 
447  return( Distance_A );
448  }
449  else
450  {
451  Line_Point = Line_B;
452 
453  return( Distance_B );
454  }
455 }
456 
457 
459 // //
460 // //
461 // //
463 
464 //---------------------------------------------------------
465 bool SG_Get_Triangle_CircumCircle(TSG_Point Triangle[3], TSG_Point &Point, double &Radius)
466 {
467  #define A Triangle[0]
468  #define B Triangle[1]
469  #define C Triangle[2]
470 
471  //-----------------------------------------------------
472  TSG_Point AB, AC, AB_M, AC_M, AB_N, AC_N;
473 
474  AB.x = B.x - A.x;
475  AB.y = B.y - A.y;
476  AB_M.x = A.x + AB.x / 2.;
477  AB_M.y = A.y + AB.y / 2.;
478  AB_N.x = AB_M.x - AB.y;
479  AB_N.y = AB_M.y + AB.x;
480 
481  AC.x = C.x - A.x;
482  AC.y = C.y - A.y;
483  AC_M.x = A.x + AC.x / 2.;
484  AC_M.y = A.y + AC.y / 2.;
485  AC_N.x = AC_M.x - AC.y;
486  AC_N.y = AC_M.y + AC.x;
487 
488  if( SG_Get_Crossing(Point, AB_M, AB_N, AC_M, AC_N, false) )
489  {
490  AB.x = A.x - Point.x;
491  AB.y = A.y - Point.y;
492 
493  Radius = sqrt(AB.x*AB.x + AB.y*AB.y);
494 
495  return( true );
496  }
497 
498  return( false );
499 
500  //-----------------------------------------------------
501  #undef A
502  #undef B
503  #undef C
504 }
505 
506 
508 // //
509 // //
510 // //
512 
513 //---------------------------------------------------------
514 double SG_Get_Polygon_Area(TSG_Point *Points, int nPoints)
515 {
516  double Area = 0.;
517 
518  if( nPoints >= 3 )
519  {
520  TSG_Point *iP = Points, *jP = Points + nPoints - 1;
521 
522  for(int i=0; i<nPoints; i++, jP=iP++)
523  {
524  Area += (jP->x * iP->y) - (iP->x * jP->y);
525  }
526 
527  Area /= 2.;
528  }
529 
530  return( Area );
531 }
532 
533 //---------------------------------------------------------
534 double SG_Get_Polygon_Area(const CSG_Points &Points)
535 {
536  double Area = 0.;
537 
538  if( Points.Get_Count() >= 3 )
539  {
540  for(sLong i=0, j=Points.Get_Count()-1; i<Points.Get_Count(); j=i++)
541  {
542  Area += (Points[j].x * Points[i].y)
543  - (Points[i].x * Points[j].y);
544  }
545 
546  Area /= 2.;
547  }
548 
549  return( Area );
550 }
551 
552 
554 // //
555 // //
556 // //
558 
559 //---------------------------------------------------------
M_DEG_TO_RAD
#define M_DEG_TO_RAD
Definition: mat_tools.h:111
SG_Get_Triangle_CircumCircle
bool SG_Get_Triangle_CircumCircle(TSG_Point Triangle[3], TSG_Point &Point, double &Radius)
Definition: geo_functions.cpp:465
SG_Get_Crossing_InRegion
bool SG_Get_Crossing_InRegion(TSG_Point &Crossing, const TSG_Point &a, const TSG_Point &b, const TSG_Rect &Region)
Definition: geo_functions.cpp:342
A
#define A
SG_Get_Polygon_Area
double SG_Get_Polygon_Area(TSG_Point *Points, int nPoints)
Definition: geo_functions.cpp:514
SG_Get_Distance_To_Line
double SG_Get_Distance_To_Line(const TSG_Point &Point, const TSG_Point &Line_A, const TSG_Point &Line_B, bool bExactMatch)
Definition: geo_functions.cpp:413
C
#define C
SSG_Point_3D
Definition: geo_tools.h:264
SSG_Rect::xMax
double xMax
Definition: geo_tools.h:465
SG_Get_Angle_Of_Direction
double SG_Get_Angle_Of_Direction(double dx, double dy)
Definition: geo_functions.cpp:216
SSG_Point
Definition: geo_tools.h:128
SG_Get_Nearest_Point_On_Line
double SG_Get_Nearest_Point_On_Line(const TSG_Point &Point, const TSG_Point &Line_A, const TSG_Point &Line_B, TSG_Point &Line_Point, bool bExactMatch)
Definition: geo_functions.cpp:421
SSG_Rect::xMin
double xMin
Definition: geo_tools.h:465
SSG_Rect
Definition: geo_tools.h:464
mat_tools.h
CSG_Point
Definition: geo_tools.h:135
CSG_Points::Get_Count
sLong Get_Count(void) const
Definition: geo_tools.h:200
SG_Get_Distance_Polar
double SG_Get_Distance_Polar(double aLon, double aLat, double bLon, double bLat, double a, double e, bool bDegree)
Definition: geo_functions.cpp:157
SG_Get_Distance
double SG_Get_Distance(double ax, double ay, double bx, double by, bool bPolar)
Definition: geo_functions.cpp:103
sLong
signed long long sLong
Definition: api_core.h:158
SG_Get_Length
double SG_Get_Length(double dx, double dy)
Definition: geo_functions.cpp:97
SG_Get_Square
double SG_Get_Square(double Value)
Definition: mat_tools.cpp:70
M_PI_180
#define M_PI_180
Definition: mat_tools.h:104
SG_is_Angle_Between
bool SG_is_Angle_Between(double Angle, double Angle_Min, double Angle_Max, bool bCheckRange)
Definition: geo_functions.cpp:257
SSG_Rect::yMax
double yMax
Definition: geo_tools.h:465
SG_Is_Equal
bool SG_Is_Equal(double a, double b, double epsilon)
Definition: geo_functions.cpp:64
CSG_Points
Definition: geo_tools.h:184
SG_Get_Angle_Difference
double SG_Get_Angle_Difference(double a, double b)
Definition: geo_functions.cpp:247
M_GET_MIN
#define M_GET_MIN(a, b)
Definition: mat_tools.h:143
SG_Is_Point_On_Line
bool SG_Is_Point_On_Line(const TSG_Point &Point, const TSG_Point &Line_A, const TSG_Point &Line_B, bool bExactMatch, double Epsilon)
Definition: geo_functions.cpp:394
M_GET_MAX
#define M_GET_MAX(a, b)
Definition: mat_tools.h:144
B
#define B
M_PI_090
#define M_PI_090
Definition: mat_tools.h:102
SSG_Point::x
double x
Definition: geo_tools.h:129
SSG_Rect::yMin
double yMin
Definition: geo_tools.h:465
M_PI_360
#define M_PI_360
Definition: mat_tools.h:108
SSG_Point::y
double y
Definition: geo_tools.h:129
SG_IS_BETWEEN
#define SG_IS_BETWEEN(a, x, b)
Definition: geo_tools.h:90
SG_Is_Between
bool SG_Is_Between(double x, double a, double b, double epsilon)
Definition: geo_functions.cpp:77
SG_Get_Crossing
bool SG_Get_Crossing(TSG_Point &Crossing, const TSG_Point &a1, const TSG_Point &a2, const TSG_Point &b1, const TSG_Point &b2, bool bExactMatch)
Definition: geo_functions.cpp:280
geo_tools.h