SAGA API v9.10
Loading...
Searching...
No Matches
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//---------------------------------------------------------
64bool SG_Is_Equal(double a, double b, double epsilon)
65{
66 return( fabs(a - b) <= epsilon );
67}
68
69//---------------------------------------------------------
70bool 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//---------------------------------------------------------
77bool 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
83bool 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//---------------------------------------------------------
97double SG_Get_Length(double dx, double dy)
98{
99 return( sqrt(dx*dx + dy*dy) );
100}
101
102//---------------------------------------------------------
103double 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//---------------------------------------------------------
109double 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//---------------------------------------------------------
115double 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//---------------------------------------------------------
124double 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//---------------------------------------------------------
135double 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//---------------------------------------------------------
145double 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//---------------------------------------------------------
157double 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//---------------------------------------------------------
203double 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//---------------------------------------------------------
216double 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//---------------------------------------------------------
229double 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//---------------------------------------------------------
235double SG_Get_Angle_Of_Direction(const TSG_Point &A)
236{
237 return( SG_Get_Angle_Of_Direction(A.x, A.y) );
238}
239
240//---------------------------------------------------------
241double 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//---------------------------------------------------------
247double 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//---------------------------------------------------------
257bool 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//---------------------------------------------------------
280bool 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//---------------------------------------------------------
342bool 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//---------------------------------------------------------
394bool 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//---------------------------------------------------------
413double 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//---------------------------------------------------------
421double 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//---------------------------------------------------------
465bool 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//---------------------------------------------------------
514double 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//---------------------------------------------------------
534double 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//---------------------------------------------------------
sLong
signed long long sLong
Definition api_core.h:158
CSG_Point
Definition geo_tools.h:135
CSG_Points
Definition geo_tools.h:185
CSG_Points::Get_Count
sLong Get_Count(void) const
Definition geo_tools.h:201
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
B
#define B
SG_Get_Polygon_Area
double SG_Get_Polygon_Area(TSG_Point *Points, int nPoints)
Definition geo_functions.cpp:514
SG_Get_Angle_Difference
double SG_Get_Angle_Difference(double a, double b)
Definition geo_functions.cpp:247
SG_Get_Distance
double SG_Get_Distance(double ax, double ay, double bx, double by, bool bPolar)
Definition geo_functions.cpp:103
SG_Get_Length
double SG_Get_Length(double dx, double dy)
Definition geo_functions.cpp:97
SG_is_Angle_Between
bool SG_is_Angle_Between(double Angle, double Angle_Min, double Angle_Max, bool bCheckRange)
Definition geo_functions.cpp:257
SG_Get_Angle_Of_Direction
double SG_Get_Angle_Of_Direction(double dx, double dy)
Definition geo_functions.cpp:216
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
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
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
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
A
#define A
SG_Is_Between
bool SG_Is_Between(double x, double a, double b, double epsilon)
Definition geo_functions.cpp:77
SG_Get_Triangle_CircumCircle
bool SG_Get_Triangle_CircumCircle(TSG_Point Triangle[3], TSG_Point &Point, double &Radius)
Definition geo_functions.cpp:465
C
#define C
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_Is_Equal
bool SG_Is_Equal(double a, double b, double epsilon)
Definition geo_functions.cpp:64
geo_tools.h
TSG_Point
struct SSG_Point TSG_Point
TSG_Point_3D
struct SSG_Point_3D TSG_Point_3D
TSG_Rect
struct SSG_Rect TSG_Rect
SG_IS_BETWEEN
#define SG_IS_BETWEEN(a, x, b)
Definition geo_tools.h:90
SG_Get_Square
double SG_Get_Square(double Value)
Definition mat_tools.cpp:70
mat_tools.h
M_PI_360
#define M_PI_360
Definition mat_tools.h:106
M_PI_090
#define M_PI_090
Definition mat_tools.h:100
M_GET_MAX
#define M_GET_MAX(a, b)
Definition mat_tools.h:142
M_PI_180
#define M_PI_180
Definition mat_tools.h:102
M_GET_MIN
#define M_GET_MIN(a, b)
Definition mat_tools.h:141
M_DEG_TO_RAD
#define M_DEG_TO_RAD
Definition mat_tools.h:109
SSG_Point::x
double x
Definition geo_tools.h:129
SSG_Point::y
double y
Definition geo_tools.h:129
SSG_Rect::xMin
double xMin
Definition geo_tools.h:468
SSG_Rect::xMax
double xMax
Definition geo_tools.h:468
SSG_Rect::yMin
double yMin
Definition geo_tools.h:468
SSG_Rect::yMax
double yMax
Definition geo_tools.h:468