1
2 ///////////////////////////////////////////////////////////
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 // //
50 ///////////////////////////////////////////////////////////
51
52 //---------------------------------------------------------
53 #include "geo_tools.h"
54 #include "mat_tools.h"
55
56
57 ///////////////////////////////////////////////////////////
58 // //
59 // //
60 // //
61 ///////////////////////////////////////////////////////////
62
63 //---------------------------------------------------------
SG_Is_Equal(double a,double b,double epsilon)64 bool SG_Is_Equal(double a, double b, double epsilon)
65 {
66 return( fabs(a - b) <= epsilon );
67 }
68
69 //---------------------------------------------------------
SG_Is_Equal(const TSG_Point & A,const TSG_Point & B,double epsilon)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 //---------------------------------------------------------
SG_Is_Between(double x,double a,double b,double epsilon)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
SG_Is_Between(const TSG_Point & Point,const TSG_Point & Corner_A,const TSG_Point & Corner_B,double epsilon)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
90 ///////////////////////////////////////////////////////////
91 // //
92 // //
93 // //
94 ///////////////////////////////////////////////////////////
95
96 //---------------------------------------------------------
SG_Get_Length(double dx,double dy)97 double SG_Get_Length(double dx, double dy)
98 {
99 return( sqrt(dx*dx + dy*dy) );
100 }
101
102 //---------------------------------------------------------
SG_Get_Distance(double ax,double ay,double bx,double by,bool bPolar)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 //---------------------------------------------------------
SG_Get_Distance(const TSG_Point & A,const TSG_Point & B,bool bPolar)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 //---------------------------------------------------------
SG_Get_Distance(double ax,double ay,double bx,double by)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 //---------------------------------------------------------
SG_Get_Distance(const TSG_Point & A,const TSG_Point & B)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 //---------------------------------------------------------
SG_Get_Distance(double ax,double ay,double az,double bx,double by,double bz)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 //---------------------------------------------------------
SG_Get_Distance(const TSG_Point_Z & A,const TSG_Point_Z & B)145 double SG_Get_Distance(const TSG_Point_Z &A, const TSG_Point_Z &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 //---------------------------------------------------------
SG_Get_Distance_Polar(double aLon,double aLat,double bLon,double bLat,double a,double e,bool bDegree)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 //---------------------------------------------------------
SG_Get_Distance_Polar(const TSG_Point & A,const TSG_Point & B,double a,double e,bool bDegree)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
209 ///////////////////////////////////////////////////////////
210 // //
211 // //
212 // //
213 ///////////////////////////////////////////////////////////
214
215 //---------------------------------------------------------
SG_Get_Angle_Of_Direction(double dx,double dy)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 //---------------------------------------------------------
SG_Get_Angle_Of_Direction(double ax,double ay,double bx,double by)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 //---------------------------------------------------------
SG_Get_Angle_Of_Direction(const TSG_Point & A)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 //---------------------------------------------------------
SG_Get_Angle_Of_Direction(const TSG_Point & A,const TSG_Point & B)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 //---------------------------------------------------------
SG_Get_Angle_Difference(double a,double b)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 //---------------------------------------------------------
SG_is_Angle_Between(double Angle,double Angle_Min,double Angle_Max,bool bCheckRange)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
273 ///////////////////////////////////////////////////////////
274 // //
275 // //
276 // //
277 ///////////////////////////////////////////////////////////
278
279 //---------------------------------------------------------
SG_Get_Crossing(TSG_Point & Crossing,const TSG_Point & a1,const TSG_Point & a2,const TSG_Point & b1,const TSG_Point & b2,bool bExactMatch)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 //---------------------------------------------------------
SG_Get_Crossing_InRegion(TSG_Point & Crossing,const TSG_Point & a,const TSG_Point & b,const TSG_Rect & Region)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
387 ///////////////////////////////////////////////////////////
388 // //
389 // //
390 // //
391 ///////////////////////////////////////////////////////////
392
393 //---------------------------------------------------------
SG_Is_Point_On_Line(const TSG_Point & Point,const TSG_Point & Ln_A,const TSG_Point & Ln_B,bool bExactMatch,double Epsilon)394 bool SG_Is_Point_On_Line(const TSG_Point &Point, const TSG_Point &Ln_A, const TSG_Point &Ln_B, bool bExactMatch, double Epsilon)
395 {
396 if( SG_Is_Equal(Ln_B.x, Ln_A.x, Epsilon) ) // vertical line
397 {
398 return( SG_Is_Between(Point.y, Ln_A.y, Ln_B.y, Epsilon) && (!bExactMatch || SG_Is_Between(Point.x, Ln_A.x, Ln_B.x, Epsilon)) );
399 }
400
401 if( bExactMatch && !SG_Is_Between(Point, Ln_A, Ln_B, Epsilon) )
402 {
403 return( false );
404 }
405
406 double b = (Ln_B.y - Ln_A.y) / (Ln_B.x - Ln_A.x);
407 double a = Ln_A.y - b * Ln_A.x;
408
409 return( SG_Is_Equal(Point.y, a + b * Point.x, Epsilon) );
410 }
411
412 //---------------------------------------------------------
SG_Get_Nearest_Point_On_Line(const TSG_Point & Point,const TSG_Point & Ln_A,const TSG_Point & Ln_B,TSG_Point & Ln_Point,bool bExactMatch)413 double SG_Get_Nearest_Point_On_Line(const TSG_Point &Point, const TSG_Point &Ln_A, const TSG_Point &Ln_B, TSG_Point &Ln_Point, bool bExactMatch)
414 {
415 double dx, dy, Distance, d;
416 TSG_Point Point_B;
417
418 Point_B.x = Point.x - (Ln_B.y - Ln_A.y);
419 Point_B.y = Point.y + (Ln_B.x - Ln_A.x);
420
421 if( SG_Get_Crossing(Ln_Point, Ln_A, Ln_B, Point, Point_B, false) )
422 {
423 if( !bExactMatch || (bExactMatch && SG_IS_BETWEEN(Ln_A.x, Ln_Point.x, Ln_B.x) && SG_IS_BETWEEN(Ln_A.y, Ln_Point.y, Ln_B.y)) )
424 {
425 dx = Point.x - Ln_Point.x;
426 dy = Point.y - Ln_Point.y;
427 Distance = sqrt(dx*dx + dy*dy);
428 }
429 else
430 {
431 dx = Point.x - Ln_A.x;
432 dy = Point.y - Ln_A.y;
433 d = sqrt(dx*dx + dy*dy);
434
435 dx = Point.x - Ln_B.x;
436 dy = Point.y - Ln_B.y;
437 Distance = sqrt(dx*dx + dy*dy);
438
439 if( d < Distance )
440 {
441 Distance = d;
442 Ln_Point = Ln_A;
443 }
444 else
445 {
446 Ln_Point = Ln_B;
447 }
448 }
449
450 return( Distance );
451 }
452
453 return( -1. );
454 }
455
456
457 ///////////////////////////////////////////////////////////
458 // //
459 // //
460 // //
461 ///////////////////////////////////////////////////////////
462
463 //---------------------------------------------------------
SG_Get_Triangle_CircumCircle(TSG_Point Triangle[3],TSG_Point & Point,double & Radius)464 bool SG_Get_Triangle_CircumCircle(TSG_Point Triangle[3], TSG_Point &Point, double &Radius)
465 {
466 #define A Triangle[0]
467 #define B Triangle[1]
468 #define C Triangle[2]
469
470 //-----------------------------------------------------
471 TSG_Point AB, AC, AB_M, AC_M, AB_N, AC_N;
472
473 AB.x = B.x - A.x;
474 AB.y = B.y - A.y;
475 AB_M.x = A.x + AB.x / 2.;
476 AB_M.y = A.y + AB.y / 2.;
477 AB_N.x = AB_M.x - AB.y;
478 AB_N.y = AB_M.y + AB.x;
479
480 AC.x = C.x - A.x;
481 AC.y = C.y - A.y;
482 AC_M.x = A.x + AC.x / 2.;
483 AC_M.y = A.y + AC.y / 2.;
484 AC_N.x = AC_M.x - AC.y;
485 AC_N.y = AC_M.y + AC.x;
486
487 if( SG_Get_Crossing(Point, AB_M, AB_N, AC_M, AC_N, false) )
488 {
489 AB.x = A.x - Point.x;
490 AB.y = A.y - Point.y;
491
492 Radius = sqrt(AB.x*AB.x + AB.y*AB.y);
493
494 return( true );
495 }
496
497 return( false );
498
499 //-----------------------------------------------------
500 #undef A
501 #undef B
502 #undef C
503 }
504
505
506 ///////////////////////////////////////////////////////////
507 // //
508 // //
509 // //
510 ///////////////////////////////////////////////////////////
511
512 //---------------------------------------------------------
SG_Get_Polygon_Area(TSG_Point * Points,int nPoints)513 double SG_Get_Polygon_Area(TSG_Point *Points, int nPoints)
514 {
515 double Area = 0.;
516
517 if( nPoints >= 3 )
518 {
519 int i;
520 TSG_Point *jP, *iP;
521
522 for(i=0, iP=Points, jP=Points+nPoints-1; 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 //---------------------------------------------------------
SG_Get_Polygon_Area(const CSG_Points & Points)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(int i=0, j=Points.Get_Count()-1; i<Points.Get_Count(); j=i++)
541 {
542 Area += (Points.Get_X(j) * Points.Get_Y(i))
543 - (Points.Get_X(i) * Points.Get_Y(j));
544 }
545
546 Area /= 2.;
547 }
548
549 return( Area );
550 }
551
552
553 ///////////////////////////////////////////////////////////
554 // //
555 // //
556 // //
557 ///////////////////////////////////////////////////////////
558
559 //---------------------------------------------------------
560