SAGA API  v9.5
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 //---------------------------------------------------------
73 int 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 //---------------------------------------------------------
94 bool 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 //---------------------------------------------------------
161 bool 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 //---------------------------------------------------------
373 int 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 //---------------------------------------------------------
CSG_Rect::Assign
CSG_Rect & Assign(double xMin, double yMin, double xMax, double yMax)
Definition: geo_classes.cpp:727
CSG_TIN::TTIN_Edge::p1
int p1
Definition: tin.h:278
CSG_Table::Get_Index
sLong Get_Index(sLong Index) const
Definition: table.h:397
CSG_TIN::TTIN_Triangle
Definition: tin.h:283
CSG_TIN_Node::Get_X
double Get_X(void)
Definition: tin.h:101
CSG_TIN::Del_Node
bool Del_Node(sLong Index, bool bUpdateNow)
Definition: tin.cpp:428
CSG_TIN::TTIN_Triangle::p2
int p2
Definition: tin.h:284
CSG_TIN::Get_Node
CSG_TIN_Node * Get_Node(sLong Index) const
Definition: tin.h:263
SG_Malloc
SAGA_API_DLL_EXPORT void * SG_Malloc(size_t size)
Definition: api_memory.cpp:65
CSG_TIN::Get_Node_Count
sLong Get_Node_Count(void) const
Definition: tin.h:262
SG_Free
SAGA_API_DLL_EXPORT void SG_Free(void *memblock)
Definition: api_memory.cpp:83
tin.h
CSG_TIN::_Destroy_Triangles
bool _Destroy_Triangles(void)
Definition: tin.cpp:269
CSG_TIN::TTIN_Edge
Definition: tin.h:277
CSG_Rect::Get_XRange
double Get_XRange(void) const
Definition: geo_tools.h:507
CSG_TIN::TTIN_Triangle::p1
int p1
Definition: tin.h:284
CSG_TIN::TTIN_Triangle::p3
int p3
Definition: tin.h:284
sLong
signed long long sLong
Definition: api_core.h:158
CSG_TIN_Node::Get_Y
double Get_Y(void)
Definition: tin.h:102
SG_TIN_Compare
int SG_TIN_Compare(const void *pp1, const void *pp2)
Definition: tin_triangulation.cpp:73
CSG_Rect::Union
CSG_Rect & Union(double x, double y)
Definition: geo_classes.cpp:824
IS_IDENTICAL
#define IS_IDENTICAL(a, b)
Definition: tin_triangulation.cpp:370
SG_UI_Process_Set_Progress
bool SG_UI_Process_Set_Progress(int Position, int Range)
Definition: api_callback.cpp:256
SSG_Point::x
double x
Definition: geo_tools.h:129
CSG_TIN_Node
Definition: tin.h:93
CSG_Rect::Get_XCenter
double Get_XCenter(void) const
Definition: geo_tools.h:517
CSG_TIN::m_Extent
CSG_Rect m_Extent
Definition: tin.h:295
CSG_TIN::_Add_Triangle
CSG_TIN_Triangle * _Add_Triangle(CSG_TIN_Node *a, CSG_TIN_Node *b, CSG_TIN_Node *c)
Definition: tin.cpp:476
SSG_Point::y
double y
Definition: geo_tools.h:129
CSG_Rect::Get_YCenter
double Get_YCenter(void) const
Definition: geo_tools.h:518
SG_UI_Process_Set_Ready
bool SG_UI_Process_Set_Ready(void)
Definition: api_callback.cpp:306
SG_Realloc
SAGA_API_DLL_EXPORT void * SG_Realloc(void *memblock, size_t size)
Definition: api_memory.cpp:77
CSG_TIN::_Triangulate
bool _Triangulate(void)
Definition: tin_triangulation.cpp:94
CSG_Rect::Get_YRange
double Get_YRange(void) const
Definition: geo_tools.h:508
CSG_TIN::_Destroy_Edges
bool _Destroy_Edges(void)
Definition: tin.cpp:253
CSG_TIN::TTIN_Edge::p2
int p2
Definition: tin.h:278
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