SAGA API v9.10
Loading...
Searching...
No Matches
tin_triangulation.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_triangulation.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//
54// The Delaunay Triangulation algorithm used here is based
55// on Paul Bourke's C source codes, which can be found at:
56//
57// http://astronomy.swin.edu.au/~pbourke/
58//
59//---------------------------------------------------------
60
61
62//---------------------------------------------------------
63#include "tin.h"
64
65
67// //
68// //
69// //
71
72//---------------------------------------------------------
73int SG_TIN_Compare(const void *pp1, const void *pp2)
74{
75 CSG_TIN_Node *p1 = *((CSG_TIN_Node **)pp1);
76 CSG_TIN_Node *p2 = *((CSG_TIN_Node **)pp2);
77
78 if( p1->Get_X() < p2->Get_X() ) { return( -1 ); }
79 if( p1->Get_X() > p2->Get_X() ) { return( 1 ); }
80 if( p1->Get_Y() < p2->Get_Y() ) { return( -1 ); }
81 if( p1->Get_Y() > p2->Get_Y() ) { return( 1 ); }
82
83 return( 0 );
84}
85
86
88// //
89// //
90// //
92
93//---------------------------------------------------------
94bool CSG_TIN::_Triangulate(void)
95{
96 _Destroy_Edges(); _Destroy_Triangles();
97
98 //-----------------------------------------------------
99 CSG_TIN_Node **Nodes = (CSG_TIN_Node **)SG_Malloc((Get_Node_Count() + 3l) * sizeof(CSG_TIN_Node *));
100
101 for(sLong i=0; i<Get_Node_Count(); i++)
102 {
103 Nodes[i] = Get_Node(i); Nodes[i]->_Del_Relations();
104 }
105
106 //-----------------------------------------------------
107 qsort(Nodes, Get_Node_Count(), sizeof(CSG_TIN_Node *), SG_TIN_Compare);
108
109 for(sLong i=0, j=0, n=Get_Node_Count(); j<n; i++) // remove duplicates
110 {
111 Nodes[i] = Nodes[j++];
112
113 while( j < n
114 && Nodes[i]->Get_X() == Nodes[j]->Get_X()
115 && Nodes[i]->Get_Y() == Nodes[j]->Get_Y() )
116 {
117 Del_Node(Nodes[j++]->Get_Index(), false);
118 }
119 }
120
121 //-----------------------------------------------------
122 for(sLong i=Get_Node_Count(); i<Get_Node_Count()+3l; i++)
123 {
124 Nodes[i] = new CSG_TIN_Node(this, 0);
125 }
126
127 //-----------------------------------------------------
128 TTIN_Triangle *Triangles = (TTIN_Triangle *)SG_Malloc(Get_Node_Count() * 3l * sizeof(TTIN_Triangle));
129
130 int nTriangles; bool bResult = _Triangulate(Nodes, (int)Get_Node_Count(), Triangles, nTriangles);
131
132 if( bResult )
133 {
134 for(int i=0; i<nTriangles && SG_UI_Process_Set_Progress(i, nTriangles); i++)
135 {
136 _Add_Triangle(Nodes[Triangles[i].p1], Nodes[Triangles[i].p2], Nodes[Triangles[i].p3]);
137 }
138 }
139
140 SG_Free(Triangles);
141
142 //-----------------------------------------------------
143 for(sLong i=Get_Node_Count(); i<Get_Node_Count()+3l; i++)
144 {
145 delete(Nodes[i]);
146 }
147
148 SG_Free(Nodes);
149
150 SG_UI_Process_Set_Ready();
151
152 return( bResult );
153}
154
155
157// //
159
160//---------------------------------------------------------
161bool CSG_TIN::_Triangulate(CSG_TIN_Node **Points, int nPoints, TTIN_Triangle *Triangles, int &nTriangles)
162{
163 if( nPoints >= 3 ) // Update extent...
164 {
165 m_Extent.Assign(Points[0]->Get_Point(), Points[0]->Get_Point());
166
167 for(int i=1; i<nPoints; i++)
168 {
169 m_Extent.Union(Points[i]->Get_Point());
170 }
171 }
172 else
173 {
174 return( false );
175 }
176
177 int status = 0, emax = 200, trimax = 4 * nPoints;
178
179 //-----------------------------------------------------
180 int *complete = (int *)SG_Malloc(trimax * sizeof(int)); // Allocate memory for the completeness list, flag for each triangle
181
182 if( complete == NULL )
183 {
184 return( false );
185 }
186
187 //-----------------------------------------------------
188 TTIN_Edge *edges = (TTIN_Edge *)SG_Malloc(emax * sizeof(TTIN_Edge)); // Allocate memory for the edge list
189
190 if( edges == NULL )
191 {
192 SG_Free(complete); return( false );
193 }
194
195 //-----------------------------------------------------
196 // Find the maximum and minimum vertex bounds.
197 // This is to allow calculation of the bounding triangle
198 // Set up the supertriangle
199 // This is a triangle which encompasses all the sample points.
200 // The supertriangle coordinates are added to the end of the
201 // vertex list. The supertriangle is the first triangle in
202 // the triangle list.
203
204 double dmax = m_Extent.Get_XRange() > m_Extent.Get_YRange() ? m_Extent.Get_XRange() : m_Extent.Get_YRange();
205
206 Points[nPoints + 0]->m_Point.x = m_Extent.Get_XCenter() - 20 * dmax;
207 Points[nPoints + 1]->m_Point.x = m_Extent.Get_XCenter();
208 Points[nPoints + 2]->m_Point.x = m_Extent.Get_XCenter() + 20 * dmax;
209
210 Points[nPoints + 0]->m_Point.y = m_Extent.Get_YCenter() - dmax;
211 Points[nPoints + 1]->m_Point.y = m_Extent.Get_YCenter() + 20 * dmax;
212 Points[nPoints + 2]->m_Point.y = m_Extent.Get_YCenter() - dmax;
213
214 Triangles[0].p1 = nPoints;
215 Triangles[0].p2 = nPoints + 1;
216 Triangles[0].p3 = nPoints + 2;
217
218 complete [0] = false;
219
220 nTriangles = 1;
221
222 //-----------------------------------------------------
223 // Include each point one at a time into the existing mesh
224 for(int i=0; i<nPoints && SG_UI_Process_Set_Progress(i, nPoints); i++)
225 {
226 double xp = Points[i]->Get_X();
227 double yp = Points[i]->Get_Y();
228
229 int nedge = 0;
230
231 //-------------------------------------------------
232 // Set up the edge buffer.
233 // If the point (xp,yp) lies inside the circumcircle then the
234 // three edges of that triangle are added to the edge buffer
235 // and that triangle is removed.
236 for(int j=0; j<nTriangles; j++)
237 {
238 if( complete[j] )
239 {
240 continue;
241 }
242
243 double x1 = Points[Triangles[j].p1]->Get_X();
244 double y1 = Points[Triangles[j].p1]->Get_Y();
245 double x2 = Points[Triangles[j].p2]->Get_X();
246 double y2 = Points[Triangles[j].p2]->Get_Y();
247 double x3 = Points[Triangles[j].p3]->Get_X();
248 double y3 = Points[Triangles[j].p3]->Get_Y();
249
250 double xc, yc, r; int inside = _CircumCircle(xp, yp, x1, y1, x2, y2, x3, y3, &xc, &yc, &r);
251
252 if( xc + r < xp )
253 {
254 complete[j] = true;
255 }
256
257 if( inside )
258 {
259 if( nedge + 3 >= emax ) // Check that we haven't exceeded the edge list size
260 {
261 emax += 100; TTIN_Edge *_edges = (TTIN_Edge *)SG_Realloc(edges, emax * sizeof(TTIN_Edge));
262
263 if( _edges == NULL )
264 {
265 SG_Free(edges); SG_Free(complete); return( false );
266 }
267
268 edges = _edges;
269 }
270
271 edges[nedge + 0].p1 = Triangles[j].p1;
272 edges[nedge + 0].p2 = Triangles[j].p2;
273 edges[nedge + 1].p1 = Triangles[j].p2;
274 edges[nedge + 1].p2 = Triangles[j].p3;
275 edges[nedge + 2].p1 = Triangles[j].p3;
276 edges[nedge + 2].p2 = Triangles[j].p1;
277
278 nedge += 3;
279
280 Triangles[j] = Triangles[nTriangles - 1];
281 complete [j] = complete [nTriangles - 1];
282
283 nTriangles--; j--;
284 }
285 }
286
287 //-------------------------------------------------
288 // Tag multiple edges
289 // Note: if all triangles are specified anticlockwise then all
290 // interior edges are opposite pointing in direction.
291 for(int j=0; j<nedge-1; j++)
292 {
293 for(int k=j+1; k<nedge; k++)
294 {
295 if( (edges[j].p1 == edges[k].p2) && (edges[j].p2 == edges[k].p1) )
296 {
297 edges[j].p1 = -1;
298 edges[j].p2 = -1;
299 edges[k].p1 = -1;
300 edges[k].p2 = -1;
301 }
302
303 // Shouldn't need the following, see note above
304 if( (edges[j].p1 == edges[k].p1) && (edges[j].p2 == edges[k].p2) )
305 {
306 edges[j].p1 = -1;
307 edges[j].p2 = -1;
308 edges[k].p1 = -1;
309 edges[k].p2 = -1;
310 }
311 }
312 }
313
314 //-------------------------------------------------
315 // Form new triangles for the current point
316 // Skipping over any tagged edges.
317 // All edges are arranged in clockwise order.
318 for(int j=0; j<nedge; j++)
319 {
320 if( edges[j].p1 < 0 || edges[j].p2 < 0 )
321 {
322 continue;
323 }
324
325 if( nTriangles >= trimax )
326 {
327 SG_Free(edges); SG_Free(complete); return( false );
328 }
329
330 Triangles[nTriangles].p1 = edges[j].p1;
331 Triangles[nTriangles].p2 = edges[j].p2;
332 Triangles[nTriangles].p3 = i;
333 complete [nTriangles] = false;
334 nTriangles++;
335 }
336 }
337
338 //-----------------------------------------------------
339 // Remove triangles with supertriangle vertices
340 // These are triangles which have a vertex number greater than nPoints
341 for(int i=0; i<nTriangles; i++)
342 {
343 if( Triangles[i].p1 >= nPoints
344 || Triangles[i].p2 >= nPoints
345 || Triangles[i].p3 >= nPoints )
346 {
347 Triangles[i] = Triangles[nTriangles - 1];
348 nTriangles--;
349 i--;
350 }
351 }
352
353 //-----------------------------------------------------
354 return( true );
355}
356
357
359// //
361
362//---------------------------------------------------------
363// Return true if a point (xp,yp) is inside the circumcircle made up
364// of the points (x1,y1), (x2,y2), (x3,y3)
365// The circumcircle centre is returned in (xc,yc) and the radius r
366// NOTE: A point on the edge is inside the circumcircle
367//
368
369//---------------------------------------------------------
370#define IS_IDENTICAL(a, b) (a == b) // #define IS_IDENTICAL(a, b) (fabs(a - b) < 0.0001)
371
372//---------------------------------------------------------
373int CSG_TIN::_CircumCircle(double xp, double yp, double x1, double y1, double x2, double y2, double x3, double y3, double *xc, double *yc, double *r)
374{
375 double m1, m2, mx1, mx2, my1, my2,
376 dx, dy, rsqr, drsqr;
377
378 // Check for coincident points
379 if( IS_IDENTICAL(y1, y2) && IS_IDENTICAL(y2, y3) )
380 {
381 return( false );
382 }
383
384 //-----------------------------------------------------
385 if( IS_IDENTICAL(y2, y1) )
386 {
387 m2 = -(x3 - x2) / (y3 - y2);
388 mx2 = (x2 + x3) / 2.0;
389 my2 = (y2 + y3) / 2.0;
390 *xc = (x2 + x1) / 2.0;
391
392 *yc = m2 * (*xc - mx2) + my2;
393 }
394 else if( IS_IDENTICAL(y3, y2) )
395 {
396 m1 = -(x2 - x1) / (y2 - y1);
397 mx1 = (x1 + x2) / 2.0;
398 my1 = (y1 + y2) / 2.0;
399 *xc = (x3 + x2) / 2.0;
400
401 *yc = m1 * (*xc - mx1) + my1;
402 }
403 else
404 {
405 m1 = -(x2 - x1) / (y2 - y1);
406 m2 = -(x3 - x2) / (y3 - y2);
407 mx1 = (x1 + x2) / 2.0;
408 mx2 = (x2 + x3) / 2.0;
409 my1 = (y1 + y2) / 2.0;
410 my2 = (y2 + y3) / 2.0;
411
412 *xc = (m1 * mx1 - m2 * mx2 + my2 - my1) / (m1 - m2);
413 *yc = m1 * (*xc - mx1) + my1;
414 }
415
416 //-----------------------------------------------------
417 dx = x2 - *xc;
418 dy = y2 - *yc;
419 rsqr = dx*dx + dy*dy;
420 *r = sqrt(rsqr);
421
422 dx = xp - *xc;
423 dy = yp - *yc;
424 drsqr = dx*dx + dy*dy;
425
426 return( drsqr <= rsqr ? 1 : 0 );
427}
428
429
431// //
432// //
433// //
435
436//---------------------------------------------------------
SG_UI_Process_Set_Ready
bool SG_UI_Process_Set_Ready(void)
Definition api_callback.cpp:305
SG_UI_Process_Set_Progress
bool SG_UI_Process_Set_Progress(int Position, int Range)
Definition api_callback.cpp:255
SG_Malloc
SAGA_API_DLL_EXPORT void * SG_Malloc(size_t size)
Definition api_memory.cpp:65
sLong
signed long long sLong
Definition api_core.h:158
SG_Free
SAGA_API_DLL_EXPORT void SG_Free(void *memblock)
Definition api_memory.cpp:83
SG_Realloc
SAGA_API_DLL_EXPORT void * SG_Realloc(void *memblock, size_t size)
Definition api_memory.cpp:77
CSG_TIN_Node
Definition tin.h:93
CSG_TIN_Node::Get_Y
double Get_Y(void)
Definition tin.h:102
CSG_TIN_Node::Get_X
double Get_X(void)
Definition tin.h:101
CSG_TIN::_Add_Triangle
CSG_TIN_Triangle * _Add_Triangle(CSG_TIN_Node *a, CSG_TIN_Node *b, CSG_TIN_Node *c)
Definition tin.cpp:481
CSG_TIN::m_Extent
CSG_Rect m_Extent
Definition tin.h:295
CSG_TIN::Get_Node_Count
sLong Get_Node_Count(void) const
Definition tin.h:262
CSG_TIN::Del_Node
bool Del_Node(sLong Index, bool bUpdateNow)
Definition tin.cpp:433
CSG_TIN::Get_Node
CSG_TIN_Node * Get_Node(sLong Index) const
Definition tin.h:263
CSG_TIN::_Destroy_Edges
bool _Destroy_Edges(void)
Definition tin.cpp:258
CSG_TIN::_Destroy_Triangles
bool _Destroy_Triangles(void)
Definition tin.cpp:274
CSG_TIN::_CircumCircle
int _CircumCircle(double xp, double yp, double x1, double y1, double x2, double y2, double x3, double y3, double *xc, double *yc, double *r)
Definition tin_triangulation.cpp:373
CSG_TIN::_Triangulate
bool _Triangulate(void)
Definition tin_triangulation.cpp:94
CSG_Table::Get_Index
sLong Get_Index(sLong Index) const
Definition table.h:405
CSG_TIN::TTIN_Edge
Definition tin.h:277
CSG_TIN::TTIN_Edge::p1
int p1
Definition tin.h:278
CSG_TIN::TTIN_Edge::p2
int p2
Definition tin.h:278
CSG_TIN::TTIN_Triangle
Definition tin.h:283
CSG_TIN::TTIN_Triangle::p2
int p2
Definition tin.h:284
CSG_TIN::TTIN_Triangle::p3
int p3
Definition tin.h:284
CSG_TIN::TTIN_Triangle::p1
int p1
Definition tin.h:284
SSG_Point::x
double x
Definition geo_tools.h:129
SSG_Point::y
double y
Definition geo_tools.h:129
SG_TIN_Compare
int SG_TIN_Compare(const void *pp1, const void *pp2)
Definition tin_triangulation.cpp:73
IS_IDENTICAL
#define IS_IDENTICAL(a, b)
Definition tin_triangulation.cpp:370