1 /**********************************************************
2 * Version $Id$
3 *********************************************************/
4
5 ///////////////////////////////////////////////////////////
6 // //
7 // SAGA //
8 // //
9 // System for Automated Geoscientific Analyses //
10 // //
11 // Application Programming Interface //
12 // //
13 // Library: SAGA_API //
14 // //
15 //-------------------------------------------------------//
16 // //
17 // tin_elements.cpp //
18 // //
19 // Copyright (C) 2005 by Olaf Conrad //
20 // //
21 //-------------------------------------------------------//
22 // //
23 // This file is part of 'SAGA - System for Automated //
24 // Geoscientific Analyses'. //
25 // //
26 // This library is free software; you can redistribute //
27 // it and/or modify it under the terms of the GNU Lesser //
28 // General Public License as published by the Free //
29 // Software Foundation, either version 2.1 of the //
30 // License, or (at your option) any later version. //
31 // //
32 // This library is distributed in the hope that it will //
33 // be useful, but WITHOUT ANY WARRANTY; without even the //
34 // implied warranty of MERCHANTABILITY or FITNESS FOR A //
35 // PARTICULAR PURPOSE. See the GNU Lesser General Public //
36 // License for more details. //
37 // //
38 // You should have received a copy of the GNU Lesser //
39 // General Public License along with this program; if //
40 // not, see <http://www.gnu.org/licenses/>. //
41 // //
42 //-------------------------------------------------------//
43 // //
44 // contact: Olaf Conrad //
45 // Institute of Geography //
46 // University of Goettingen //
47 // Goldschmidtstr. 5 //
48 // 37077 Goettingen //
49 // Germany //
50 // //
51 // e-mail: oconrad@saga-gis.org //
52 // //
53 ///////////////////////////////////////////////////////////
54
55 //---------------------------------------------------------
56
57
58 ///////////////////////////////////////////////////////////
59 // //
60 // //
61 // //
62 ///////////////////////////////////////////////////////////
63
64 //---------------------------------------------------------
65 #include "tin.h"
66
67
68 ///////////////////////////////////////////////////////////
69 // //
70 // //
71 // //
72 ///////////////////////////////////////////////////////////
73
74 //---------------------------------------------------------
CSG_TIN_Node(CSG_TIN * pOwner,int Index)75 CSG_TIN_Node::CSG_TIN_Node(CSG_TIN *pOwner, int Index)
76 : CSG_Table_Record(pOwner, Index)
77 {
78 m_Point.x = m_Point.y = 0.0;
79
80 m_Neighbors = NULL;
81 m_nNeighbors = 0;
82
83 m_Triangles = NULL;
84 m_nTriangles = 0;
85 }
86
87 //---------------------------------------------------------
~CSG_TIN_Node(void)88 CSG_TIN_Node::~CSG_TIN_Node(void)
89 {
90 _Del_Relations();
91 }
92
93 //---------------------------------------------------------
_Add_Triangle(CSG_TIN_Triangle * pTriangle)94 bool CSG_TIN_Node::_Add_Triangle(CSG_TIN_Triangle *pTriangle)
95 {
96 for(int i=0; i<m_nTriangles; i++)
97 {
98 if( m_Triangles[i] == pTriangle )
99 {
100 return( false );
101 }
102 }
103
104 m_Triangles = (CSG_TIN_Triangle **)SG_Realloc(m_Triangles, (m_nTriangles + 1) * sizeof(CSG_TIN_Triangle *));
105 m_Triangles[m_nTriangles++] = pTriangle;
106
107 // _Add_Neighbor(pTriangle->Get_Point(0));
108 // _Add_Neighbor(pTriangle->Get_Point(1));
109 // _Add_Neighbor(pTriangle->Get_Point(2));
110
111 return( true );
112 }
113
114 //---------------------------------------------------------
_Add_Neighbor(CSG_TIN_Node * pNeighbor)115 bool CSG_TIN_Node::_Add_Neighbor(CSG_TIN_Node *pNeighbor)
116 {
117 if( pNeighbor == this )
118 {
119 return( false );
120 }
121
122 for(int i=0; i<m_nNeighbors; i++)
123 {
124 if( m_Neighbors[i] == pNeighbor )
125 {
126 return( false );
127 }
128 }
129
130 m_Neighbors = (CSG_TIN_Node **)SG_Realloc(m_Neighbors, (m_nNeighbors + 1) * sizeof(CSG_TIN_Node *));
131 m_Neighbors[m_nNeighbors++] = pNeighbor;
132
133 return( true );
134 }
135
136 //---------------------------------------------------------
_Del_Relations(void)137 bool CSG_TIN_Node::_Del_Relations(void)
138 {
139 if( m_nTriangles > 0 )
140 {
141 SG_Free(m_Triangles);
142 m_Triangles = NULL;
143 m_nTriangles = 0;
144 }
145
146 if( m_nNeighbors > 0 )
147 {
148 SG_Free(m_Neighbors);
149 m_Neighbors = NULL;
150 m_nNeighbors = 0;
151 }
152
153 return( true );
154 }
155
156 //---------------------------------------------------------
Get_Gradient(int iNeighbor,int iField)157 double CSG_TIN_Node::Get_Gradient(int iNeighbor, int iField)
158 {
159 double dx, dy, dz;
160 CSG_TIN_Node *pNeighbor;
161
162 if( (pNeighbor = Get_Neighbor(iNeighbor)) != NULL )
163 {
164 dx = Get_X() - pNeighbor->Get_X();
165 dy = Get_Y() - pNeighbor->Get_Y();
166 dz = asDouble(iField) - pNeighbor->asDouble(iField);
167
168 if( (dx = sqrt(dx*dx + dy*dy)) > 0.0 )
169 {
170 return( dz / dx );
171 }
172 }
173
174 return( 0.0 );
175 }
176
177 //---------------------------------------------------------
SG_TIN_Compare_Triangle_Center(const void * pz1,const void * pz2)178 int SG_TIN_Compare_Triangle_Center(const void *pz1, const void *pz2)
179 {
180 double z1 = ((TSG_Point_Z *)pz1)->z;
181 double z2 = ((TSG_Point_Z *)pz2)->z;
182
183 return( z1 < z2 ? -1 : z1 > z2 ? 1 : 0 );
184 }
185
186 //---------------------------------------------------------
Get_Polygon(CSG_Points & Points)187 bool CSG_TIN_Node::Get_Polygon(CSG_Points &Points)
188 {
189 Points.Clear();
190
191 if( m_nTriangles >= 3 )
192 {
193 int i; CSG_Points_Z p;
194
195 for(i=0; i<m_nTriangles; i++)
196 {
197 TSG_Point c = m_Triangles[i]->Get_CircumCircle_Point();
198
199 p.Add(c.x, c.y, atan2(m_Point.y - c.y, m_Point.x - c.x));
200 }
201
202 qsort(&(p[0]), p.Get_Count(), sizeof(TSG_Point_Z), SG_TIN_Compare_Triangle_Center);
203
204 for(i=0; i<m_nTriangles; i++)
205 {
206 Points.Add(p[i].x, p[i].y);
207 }
208
209 return( true );
210 }
211
212 return( false );
213 }
214
215 //---------------------------------------------------------
Get_Polygon_Area(void)216 double CSG_TIN_Node::Get_Polygon_Area(void)
217 {
218 CSG_Points Points;
219
220 if( Get_Polygon(Points) )
221 {
222 return( SG_Get_Polygon_Area(Points) );
223 }
224
225 return( 0.0 );
226 }
227
228
229 ///////////////////////////////////////////////////////////
230 // //
231 // //
232 // //
233 ///////////////////////////////////////////////////////////
234
235 //---------------------------------------------------------
CSG_TIN_Edge(CSG_TIN_Node * a,CSG_TIN_Node * b)236 CSG_TIN_Edge::CSG_TIN_Edge(CSG_TIN_Node *a, CSG_TIN_Node *b)
237 {
238 m_Nodes[0] = a;
239 m_Nodes[1] = b;
240 }
241
242 //---------------------------------------------------------
~CSG_TIN_Edge(void)243 CSG_TIN_Edge::~CSG_TIN_Edge(void)
244 {
245 }
246
247
248 ///////////////////////////////////////////////////////////
249 // //
250 // //
251 // //
252 ///////////////////////////////////////////////////////////
253
254 //---------------------------------------------------------
CSG_TIN_Triangle(CSG_TIN_Node * a,CSG_TIN_Node * b,CSG_TIN_Node * c)255 CSG_TIN_Triangle::CSG_TIN_Triangle(CSG_TIN_Node *a, CSG_TIN_Node *b, CSG_TIN_Node *c)
256 {
257 m_Nodes[0] = a;
258 m_Nodes[1] = b;
259 m_Nodes[2] = c;
260
261 //-----------------------------------------------------
262 double xMin, yMin, xMax, yMax;
263
264 xMin = xMax = a->Get_X();
265 yMin = yMax = a->Get_Y();
266
267 if( xMin > b->Get_X() )
268 xMin = b->Get_X();
269 else if( xMax < b->Get_X() )
270 xMax = b->Get_X();
271
272 if( yMin > b->Get_Y() )
273 yMin = b->Get_Y();
274 else if( yMax < b->Get_Y() )
275 yMax = b->Get_Y();
276
277 if( xMin > c->Get_X() )
278 xMin = c->Get_X();
279 else if( xMax < c->Get_X() )
280 xMax = c->Get_X();
281
282 if( yMin > c->Get_Y() )
283 yMin = c->Get_Y();
284 else if( yMax < c->Get_Y() )
285 yMax = c->Get_Y();
286
287 m_Extent.Assign(xMin, yMin, xMax, yMax);
288
289 //-----------------------------------------------------
290 m_Area = fabs( a->Get_X() * (b->Get_Y() - c->Get_Y())
291 + b->Get_X() * (c->Get_Y() - a->Get_Y())
292 + c->Get_X() * (a->Get_Y() - b->Get_Y()) ) / 2.0;
293
294 //-----------------------------------------------------
295 TSG_Point Points[3];
296
297 Points[0] = m_Nodes[0]->Get_Point();
298 Points[1] = m_Nodes[1]->Get_Point();
299 Points[2] = m_Nodes[2]->Get_Point();
300
301 SG_Get_Triangle_CircumCircle(Points, m_Center, m_Radius);
302 }
303
304 //---------------------------------------------------------
~CSG_TIN_Triangle(void)305 CSG_TIN_Triangle::~CSG_TIN_Triangle(void)
306 {}
307
308
309 //---------------------------------------------------------
is_Containing(const TSG_Point & Point)310 bool CSG_TIN_Triangle::is_Containing(const TSG_Point &Point)
311 {
312 return( is_Containing(Point.x, Point.y) );
313 }
314
315 //---------------------------------------------------------
316 #define IS_ONLINE(A, B) (A.y == B.y && ((A.x <= x && x <= B.x) || (B.x <= x && x <= A.x)))
317
318 //---------------------------------------------------------
is_Containing(double x,double y)319 bool CSG_TIN_Triangle::is_Containing(double x, double y)
320 {
321 if( m_Extent.Contains(x, y) )
322 {
323 int nCrossings;
324 TSG_Point A, B, C;
325
326 if( (x == m_Nodes[0]->Get_Point().x && y == m_Nodes[0]->Get_Point().y)
327 || (x == m_Nodes[1]->Get_Point().x && y == m_Nodes[1]->Get_Point().y)
328 || (x == m_Nodes[2]->Get_Point().x && y == m_Nodes[2]->Get_Point().y) )
329 return( true );
330
331 if( y == m_Extent.Get_YMin() || y == m_Extent.Get_YMax() )
332 {
333 if( IS_ONLINE(m_Nodes[0]->Get_Point(), m_Nodes[1]->Get_Point())
334 || IS_ONLINE(m_Nodes[1]->Get_Point(), m_Nodes[2]->Get_Point())
335 || IS_ONLINE(m_Nodes[2]->Get_Point(), m_Nodes[0]->Get_Point()) )
336 return( true );
337 }
338
339 nCrossings = 0;
340
341 if( (y == m_Nodes[0]->Get_Point().y && x > m_Nodes[0]->Get_Point().x)
342 || (y == m_Nodes[1]->Get_Point().y && x > m_Nodes[1]->Get_Point().x)
343 || (y == m_Nodes[2]->Get_Point().y && x > m_Nodes[2]->Get_Point().x) )
344 nCrossings = -1;
345
346 A.x = m_Extent.m_rect.xMin - 1.0;
347 B.x = x;
348 A.y = B.y = y;
349
350 if( SG_Get_Crossing(C, m_Nodes[0]->Get_Point(), m_Nodes[1]->Get_Point(), A, B) )
351 nCrossings++;
352
353 if( SG_Get_Crossing(C, m_Nodes[1]->Get_Point(), m_Nodes[2]->Get_Point(), A, B) )
354 nCrossings++;
355
356 if( SG_Get_Crossing(C, m_Nodes[2]->Get_Point(), m_Nodes[0]->Get_Point(), A, B) )
357 nCrossings++;
358
359 return( nCrossings == 1 );
360 }
361
362 return( false );
363 }
364
365 //---------------------------------------------------------
Get_Value(int zField,const TSG_Point & p,double & z)366 bool CSG_TIN_Triangle::Get_Value(int zField, const TSG_Point &p, double &z)
367 {
368 return( Get_Value(zField, p.x, p.y, z) );
369 }
370
Get_Value(int zField,double x,double y,double & z)371 bool CSG_TIN_Triangle::Get_Value(int zField, double x, double y, double &z)
372 {
373 CSG_Vector B, Z(3);
374 CSG_Matrix M(3, 3), Mt;
375
376 for(int i=0; i<3; i++)
377 {
378 M[i][0] = 1.0;
379 M[i][1] = m_Nodes[i]->Get_X();
380 M[i][2] = m_Nodes[i]->Get_Y();
381 Z[i] = m_Nodes[i]->asDouble(zField);
382 }
383
384 Mt = M.Get_Transpose();
385 B = (Mt * M).Get_Inverse() * (Mt * Z);
386
387 z = B[0] + B[1] * x + B[2] * y;
388
389 return( true );
390 }
391
392 //---------------------------------------------------------
Get_Gradient(int zField,double & Decline,double & Azimuth)393 bool CSG_TIN_Triangle::Get_Gradient(int zField, double &Decline, double &Azimuth)
394 {
395 int i;
396 double x[3], y[3], z[3], A, B, C;
397
398 for(i=0; i<3; i++)
399 {
400 x[i] = m_Nodes[i]->Get_X();
401 y[i] = m_Nodes[i]->Get_Y();
402 z[i] = m_Nodes[i]->asDouble(zField);
403 }
404
405 A = z[0] * (x[1] - x[2]) + z[1] * (x[2] - x[0]) + z[2] * (x[0] - x[1]);
406 B = y[0] * (z[1] - z[2]) + y[1] * (z[2] - z[0]) + y[2] * (z[0] - z[1]);
407 C = x[0] * (y[1] - y[2]) + x[1] * (y[2] - y[0]) + x[2] * (y[0] - y[1]);
408
409 if( C != 0.0 )
410 {
411 A = - A / C;
412 B = - B / C;
413
414 Decline = atan(sqrt(A*A + B*B));
415
416 if( A != 0.0 )
417 Azimuth = M_PI_180 + atan2(B, A);
418 else
419 Azimuth = B > 0.0 ? M_PI_270 : (B < 0.0 ? M_PI_090 : -1.0);
420
421 return( true );
422 }
423
424 Decline = -1.0;
425 Azimuth = -1.0;
426
427 return( false );
428 }
429
430
431 ///////////////////////////////////////////////////////////
432 // //
433 // //
434 // //
435 ///////////////////////////////////////////////////////////
436
437 //---------------------------------------------------------
438