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