SAGA API v9.10
Loading...
Searching...
No Matches
tin_elements.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// tin_elements.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 "tin.h"
54
55
57// //
58// //
59// //
61
62//---------------------------------------------------------
63CSG_TIN_Node::CSG_TIN_Node(CSG_TIN *pOwner, sLong Index)
64 : CSG_Table_Record(pOwner, Index)
65{
66 m_Point.x = m_Point.y = 0.;
67
68 m_Neighbors = NULL;
69 m_nNeighbors = 0;
70
71 m_Triangles = NULL;
72 m_nTriangles = 0;
73}
74
75//---------------------------------------------------------
76CSG_TIN_Node::~CSG_TIN_Node(void)
77{
78 _Del_Relations();
79}
80
81//---------------------------------------------------------
82bool CSG_TIN_Node::_Add_Triangle(CSG_TIN_Triangle *pTriangle)
83{
84 for(int i=0; i<m_nTriangles; i++)
85 {
86 if( m_Triangles[i] == pTriangle )
87 {
88 return( false );
89 }
90 }
91
92 m_Triangles = (CSG_TIN_Triangle **)SG_Realloc(m_Triangles, ((uLong)m_nTriangles + 1) * sizeof(CSG_TIN_Triangle *));
93 m_Triangles[m_nTriangles++] = pTriangle;
94
95// _Add_Neighbor(pTriangle->Get_Point());
96// _Add_Neighbor(pTriangle->Get_Point(1));
97// _Add_Neighbor(pTriangle->Get_Point(2));
98
99 return( true );
100}
101
102//---------------------------------------------------------
103bool CSG_TIN_Node::_Add_Neighbor(CSG_TIN_Node *pNeighbor)
104{
105 if( pNeighbor == this )
106 {
107 return( false );
108 }
109
110 for(int i=0; i<m_nNeighbors; i++)
111 {
112 if( m_Neighbors[i] == pNeighbor )
113 {
114 return( false );
115 }
116 }
117
118 m_Neighbors = (CSG_TIN_Node **)SG_Realloc(m_Neighbors, ((uLong)m_nNeighbors + 1) * sizeof(CSG_TIN_Node *));
119 m_Neighbors[m_nNeighbors++] = pNeighbor;
120
121 return( true );
122}
123
124//---------------------------------------------------------
125bool CSG_TIN_Node::_Del_Relations(void)
126{
127 if( m_nTriangles > 0 )
128 {
129 SG_Free(m_Triangles);
130 m_Triangles = NULL;
131 m_nTriangles = 0;
132 }
133
134 if( m_nNeighbors > 0 )
135 {
136 SG_Free(m_Neighbors);
137 m_Neighbors = NULL;
138 m_nNeighbors = 0;
139 }
140
141 return( true );
142}
143
144//---------------------------------------------------------
145double CSG_TIN_Node::Get_Gradient(int iNeighbor, int iField)
146{
147 double dx, dy, dz;
148 CSG_TIN_Node *pNeighbor;
149
150 if( (pNeighbor = Get_Neighbor(iNeighbor)) != NULL )
151 {
152 dx = Get_X() - pNeighbor->Get_X();
153 dy = Get_Y() - pNeighbor->Get_Y();
154 dz = asDouble(iField) - pNeighbor->asDouble(iField);
155
156 if( (dx = sqrt(dx*dx + dy*dy)) > 0.0 )
157 {
158 return( dz / dx );
159 }
160 }
161
162 return( 0.0 );
163}
164
165//---------------------------------------------------------
166int SG_TIN_Compare_Triangle_Center(const void *pz1, const void *pz2)
167{
168 double z1 = ((TSG_Point_3D *)pz1)->z;
169 double z2 = ((TSG_Point_3D *)pz2)->z;
170
171 return( z1 < z2 ? -1 : z1 > z2 ? 1 : 0 );
172}
173
174//---------------------------------------------------------
176{
177 Points.Clear();
178
179 if( m_nTriangles >= 3 )
180 {
181 int i; CSG_Points_3D p;
182
183 for(i=0; i<m_nTriangles; i++)
184 {
185 TSG_Point c = m_Triangles[i]->Get_CircumCircle_Point();
186
187 p.Add(c.x, c.y, atan2(m_Point.y - c.y, m_Point.x - c.x));
188 }
189
190 qsort(&(p[0]), p.Get_Count(), sizeof(TSG_Point_3D), SG_TIN_Compare_Triangle_Center);
191
192 for(i=0; i<m_nTriangles; i++)
193 {
194 Points.Add(p[i].x, p[i].y);
195 }
196
197 return( true );
198 }
199
200 return( false );
201}
202
203//---------------------------------------------------------
205{
206 CSG_Points Points;
207
208 if( Get_Polygon(Points) )
209 {
210 return( SG_Get_Polygon_Area(Points) );
211 }
212
213 return( 0.0 );
214}
215
216
218// //
219// //
220// //
222
223//---------------------------------------------------------
224CSG_TIN_Edge::CSG_TIN_Edge(CSG_TIN_Node *a, CSG_TIN_Node *b)
225{
226 m_Nodes[0] = a;
227 m_Nodes[1] = b;
228}
229
230//---------------------------------------------------------
231CSG_TIN_Edge::~CSG_TIN_Edge(void)
232{
233}
234
235
237// //
238// //
239// //
241
242//---------------------------------------------------------
243CSG_TIN_Triangle::CSG_TIN_Triangle(CSG_TIN_Node *a, CSG_TIN_Node *b, CSG_TIN_Node *c)
244{
245 m_Nodes[0] = a;
246 m_Nodes[1] = b;
247 m_Nodes[2] = c;
248
249 //-----------------------------------------------------
250 double xMin, yMin, xMax, yMax;
251
252 xMin = xMax = a->Get_X();
253 yMin = yMax = a->Get_Y();
254
255 if( xMin > b->Get_X() )
256 xMin = b->Get_X();
257 else if( xMax < b->Get_X() )
258 xMax = b->Get_X();
259
260 if( yMin > b->Get_Y() )
261 yMin = b->Get_Y();
262 else if( yMax < b->Get_Y() )
263 yMax = b->Get_Y();
264
265 if( xMin > c->Get_X() )
266 xMin = c->Get_X();
267 else if( xMax < c->Get_X() )
268 xMax = c->Get_X();
269
270 if( yMin > c->Get_Y() )
271 yMin = c->Get_Y();
272 else if( yMax < c->Get_Y() )
273 yMax = c->Get_Y();
274
275 m_Extent.Assign(xMin, yMin, xMax, yMax);
276
277 //-----------------------------------------------------
278 m_Area = fabs( a->Get_X() * (b->Get_Y() - c->Get_Y())
279 + b->Get_X() * (c->Get_Y() - a->Get_Y())
280 + c->Get_X() * (a->Get_Y() - b->Get_Y()) ) / 2.0;
281
282 //-----------------------------------------------------
283 TSG_Point Points[3];
284
285 Points[0] = m_Nodes[0]->Get_Point();
286 Points[1] = m_Nodes[1]->Get_Point();
287 Points[2] = m_Nodes[2]->Get_Point();
288
289 SG_Get_Triangle_CircumCircle(Points, m_Center, m_Radius);
290}
291
292//---------------------------------------------------------
293CSG_TIN_Triangle::~CSG_TIN_Triangle(void)
294{}
295
296
297//---------------------------------------------------------
299{
300 return( is_Containing(Point.x, Point.y) );
301}
302
303//---------------------------------------------------------
304#define IS_ONLINE(A, B) (A.y == B.y && ((A.x <= x && x <= B.x) || (B.x <= x && x <= A.x)))
305
306//---------------------------------------------------------
307bool CSG_TIN_Triangle::is_Containing(double x, double y)
308{
309 if( m_Extent.Contains(x, y) )
310 {
311 int nCrossings;
312 TSG_Point A, B, C;
313
314 if( (x == m_Nodes[0]->Get_Point().x && y == m_Nodes[0]->Get_Point().y)
315 || (x == m_Nodes[1]->Get_Point().x && y == m_Nodes[1]->Get_Point().y)
316 || (x == m_Nodes[2]->Get_Point().x && y == m_Nodes[2]->Get_Point().y) )
317 return( true );
318
319 if( y == m_Extent.Get_YMin() || y == m_Extent.Get_YMax() )
320 {
321 if( IS_ONLINE(m_Nodes[0]->Get_Point(), m_Nodes[1]->Get_Point())
322 || IS_ONLINE(m_Nodes[1]->Get_Point(), m_Nodes[2]->Get_Point())
323 || IS_ONLINE(m_Nodes[2]->Get_Point(), m_Nodes[0]->Get_Point()) )
324 return( true );
325 }
326
327 nCrossings = 0;
328
329 if( (y == m_Nodes[0]->Get_Point().y && x > m_Nodes[0]->Get_Point().x)
330 || (y == m_Nodes[1]->Get_Point().y && x > m_Nodes[1]->Get_Point().x)
331 || (y == m_Nodes[2]->Get_Point().y && x > m_Nodes[2]->Get_Point().x) )
332 nCrossings = -1;
333
334 A.x = m_Extent.xMin - 1.0;
335 B.x = x;
336 A.y = B.y = y;
337
338 if( SG_Get_Crossing(C, m_Nodes[0]->Get_Point(), m_Nodes[1]->Get_Point(), A, B) )
339 nCrossings++;
340
341 if( SG_Get_Crossing(C, m_Nodes[1]->Get_Point(), m_Nodes[2]->Get_Point(), A, B) )
342 nCrossings++;
343
344 if( SG_Get_Crossing(C, m_Nodes[2]->Get_Point(), m_Nodes[0]->Get_Point(), A, B) )
345 nCrossings++;
346
347 return( nCrossings == 1 );
348 }
349
350 return( false );
351}
352
353//---------------------------------------------------------
354bool CSG_TIN_Triangle::Get_Value(int zField, const TSG_Point &p, double &z)
355{
356 return( Get_Value(zField, p.x, p.y, z) );
357}
358
359bool CSG_TIN_Triangle::Get_Value(int zField, double x, double y, double &z)
360{
361 CSG_Vector B, Z(3);
362 CSG_Matrix M(3, 3), Mt;
363
364 for(int i=0; i<3; i++)
365 {
366 M[i][0] = 1.0;
367 M[i][1] = m_Nodes[i]->Get_X();
368 M[i][2] = m_Nodes[i]->Get_Y();
369 Z[i] = m_Nodes[i]->asDouble(zField);
370 }
371
372 Mt = M.Get_Transpose();
373 B = (Mt * M).Get_Inverse() * (Mt * Z);
374
375 z = B[0] + B[1] * x + B[2] * y;
376
377 return( true );
378}
379
380//---------------------------------------------------------
381bool CSG_TIN_Triangle::Get_Gradient(int zField, double &Decline, double &Azimuth)
382{
383 int i;
384 double x[3], y[3], z[3], A, B, C;
385
386 for(i=0; i<3; i++)
387 {
388 x[i] = m_Nodes[i]->Get_X();
389 y[i] = m_Nodes[i]->Get_Y();
390 z[i] = m_Nodes[i]->asDouble(zField);
391 }
392
393 A = z[0] * (x[1] - x[2]) + z[1] * (x[2] - x[0]) + z[2] * (x[0] - x[1]);
394 B = y[0] * (z[1] - z[2]) + y[1] * (z[2] - z[0]) + y[2] * (z[0] - z[1]);
395 C = x[0] * (y[1] - y[2]) + x[1] * (y[2] - y[0]) + x[2] * (y[0] - y[1]);
396
397 if( C != 0.0 )
398 {
399 A = - A / C;
400 B = - B / C;
401
402 Decline = atan(sqrt(A*A + B*B));
403
404 if( A != 0.0 )
405 Azimuth = M_PI_180 + atan2(B, A);
406 else
407 Azimuth = B > 0.0 ? M_PI_270 : (B < 0.0 ? M_PI_090 : -1.0);
408
409 return( true );
410 }
411
412 Decline = -1.0;
413 Azimuth = -1.0;
414
415 return( false );
416}
417
418
420// //
421// //
422// //
424
425//---------------------------------------------------------
unsigned long long uLong
Definition api_core.h:159
signed long long sLong
Definition api_core.h:158
SAGA_API_DLL_EXPORT void SG_Free(void *memblock)
SAGA_API_DLL_EXPORT void * SG_Realloc(void *memblock, size_t size)
CSG_Matrix Get_Transpose(void) const
bool Add(double x, double y, double z)
sLong Get_Count(void) const
Definition geo_tools.h:338
bool Clear(void)
Definition geo_tools.h:191
bool Add(double x, double y)
bool Get_Polygon(CSG_Points &Points)
CSG_TIN_Node * Get_Neighbor(int iNeighbor)
Definition tin.h:105
double Get_Y(void)
Definition tin.h:102
double Get_Polygon_Area(void)
double Get_X(void)
Definition tin.h:101
double Get_Gradient(int iNeighbor, int iField)
bool Get_Value(int zField, const TSG_Point &p, double &z)
bool is_Containing(const TSG_Point &p)
bool Get_Gradient(int zField, double &Decline, double &Azimuth)
Definition tin.h:222
double asDouble(int Field) const
#define B
double SG_Get_Polygon_Area(TSG_Point *Points, int nPoints)
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)
#define A
bool SG_Get_Triangle_CircumCircle(TSG_Point Triangle[3], TSG_Point &Point, double &Radius)
#define C
struct SSG_Point TSG_Point
struct SSG_Point_3D TSG_Point_3D
#define M_PI_090
Definition mat_tools.h:100
#define M_PI_180
Definition mat_tools.h:102
#define M_PI_270
Definition mat_tools.h:104
double x
Definition geo_tools.h:129
double y
Definition geo_tools.h:129
#define IS_ONLINE(A, B)
int SG_TIN_Compare_Triangle_Center(const void *pz1, const void *pz2)