1 /* GTS - Library for the manipulation of triangulated surfaces
2 * Copyright (C) 1999 St�phane Popinet
3 *
4 * This library is free software; you can redistribute it and/or
5 * modify it under the terms of the GNU Library General Public
6 * License as published by the Free Software Foundation; either
7 * version 2 of the License, or (at your option) any later version.
8 *
9 * This library is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 * Library General Public License for more details.
13 *
14 * You should have received a copy of the GNU Library General Public
15 * License along with this library; if not, write to the
16 * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
17 * Boston, MA 02111-1307, USA.
18 */
19
20 #include <stdlib.h>
21 #include "config.h"
22 #ifdef HAVE_GETOPT_H
23 # include <getopt.h>
24 #endif /* HAVE_GETOPT_H */
25 #ifdef HAVE_UNISTD_H
26 # include <unistd.h>
27 #endif /* HAVE_UNISTD_H */
28 #include "gts.h"
29
30 typedef struct _CartesianGrid CartesianGrid;
31
32 struct _CartesianGrid {
33 GtsVertex *** vertices;
34 guint nx, ny;
35 gdouble xmin, xmax, ymin, ymax;
36 };
37
malloc2D(guint nx,guint ny,gulong size)38 static void ** malloc2D (guint nx, guint ny, gulong size)
39 {
40 void ** m = g_malloc (nx*sizeof (void *));
41 guint i;
42
43 for (i = 0; i < nx; i++)
44 m[i] = g_malloc0 (ny*size);
45
46 return m;
47 }
48
free2D(void ** m,guint nx)49 static void free2D (void ** m, guint nx)
50 {
51 guint i;
52
53 g_return_if_fail (m != NULL);
54
55 for (i = 0; i < nx; i++)
56 g_free (m[i]);
57 g_free (m);
58 }
59
cartesian_grid_new(guint nx,guint ny)60 static CartesianGrid * cartesian_grid_new (guint nx, guint ny)
61 {
62 CartesianGrid * grid;
63
64 grid = g_malloc (sizeof (CartesianGrid));
65 grid->vertices = (GtsVertex ***) malloc2D (nx, ny, sizeof (GtsVertex *));
66 grid->nx = nx;
67 grid->ny = ny;
68 grid->xmin = G_MAXDOUBLE;
69 grid->xmax = - G_MAXDOUBLE;
70 grid->ymin = G_MAXDOUBLE;
71 grid->ymax = - G_MAXDOUBLE;
72
73 return grid;
74 }
75
cartesian_grid_destroy(CartesianGrid * g,gboolean destroy_vertices)76 static void cartesian_grid_destroy (CartesianGrid * g,
77 gboolean destroy_vertices)
78 {
79 g_return_if_fail (g != NULL);
80
81 if (destroy_vertices) {
82 guint i, j;
83
84 gts_allow_floating_vertices = TRUE;
85 for (i = 0; i < g->nx; i++)
86 for (j = 0; j < g->ny; j++)
87 if (g->vertices[i][j])
88 gts_object_destroy (GTS_OBJECT (g->vertices[i][j]));
89 gts_allow_floating_vertices = FALSE;
90 }
91
92 free2D ((void **) g->vertices, g->nx);
93 g_free (g);
94 }
95
cartesian_grid_read(GtsVertexClass * klass,guint * line)96 static CartesianGrid * cartesian_grid_read (GtsVertexClass * klass,
97 guint * line)
98 {
99 CartesianGrid * grid;
100 guint nx, ny, line_number = 1, i, j;
101
102 g_return_val_if_fail (klass != NULL, NULL);
103
104 if (scanf ("%u %u", &nx, &ny) != 2) {
105 if (line)
106 *line = line_number;
107 return NULL;
108 }
109 line_number++;
110
111 grid = cartesian_grid_new (nx, ny);
112 for (i = 0; i < nx; i++)
113 for (j = 0; j < ny; j++) {
114 gdouble x, y, z;
115 if (scanf ("%lf %lf %lf", &x, &y, &z) != 3) {
116 cartesian_grid_destroy (grid, TRUE);
117 if (line)
118 *line = line_number;
119 return NULL;
120 }
121 if (z != 0.)
122 grid->vertices[i][j] = gts_vertex_new (klass, x, y, z);
123 if (x > grid->xmax) grid->xmax = x;
124 if (x < grid->xmin) grid->xmin = x;
125 if (y > grid->ymax) grid->ymax = y;
126 if (y < grid->ymin) grid->ymin = y;
127 line_number++;
128 }
129
130 return grid;
131 }
132
new_edge(GtsVertex * v1,GtsVertex * v2)133 static GtsEdge * new_edge (GtsVertex * v1, GtsVertex * v2)
134 {
135 GtsSegment * s = gts_vertices_are_connected (v1, v2);
136 return s == NULL ?
137 gts_edge_new (GTS_EDGE_CLASS (gts_constraint_class ()), v1, v2) :
138 GTS_EDGE (s);
139 }
140
cartesian_grid_triangulate(CartesianGrid * g,GtsSurface * s)141 static void cartesian_grid_triangulate (CartesianGrid * g,
142 GtsSurface * s)
143 {
144 gint i, j;
145 GtsVertex *** v;
146
147 g_return_if_fail (g != NULL);
148 g_return_if_fail (s != NULL);
149
150 v = g->vertices;
151 for (i = 0; i < g->nx - 1; i++)
152 for (j = 0; j < g->ny - 1; j++)
153 if (v[i][j]) {
154 if (v[i][j+1]) {
155 if (v[i+1][j+1]) {
156 GtsEdge * e1 = new_edge (v[i][j+1], v[i][j]);
157 GtsEdge * e2 =
158 gts_edge_new (GTS_EDGE_CLASS (gts_constraint_class ()),
159 v[i][j], v[i+1][j+1]);
160 GtsEdge * e3 = new_edge (v[i+1][j+1], v[i][j+1]);
161 gts_surface_add_face (s, gts_face_new (s->face_class, e1, e2, e3));
162 if (v[i+1][j]) {
163 e1 = new_edge (v[i+1][j], v[i+1][j+1]);
164 e3 = new_edge (v[i][j], v[i+1][j]);
165 gts_surface_add_face (s,
166 gts_face_new (s->face_class, e1, e2, e3));
167 }
168 }
169 else if (v[i+1][j]) {
170 GtsEdge * e1 = new_edge (v[i][j+1], v[i][j]);
171 GtsEdge * e2 = new_edge (v[i][j], v[i+1][j]);
172 GtsEdge * e3 =
173 gts_edge_new (GTS_EDGE_CLASS (gts_constraint_class ()),
174 v[i+1][j], v[i][j+1]);
175 gts_surface_add_face (s, gts_face_new (s->face_class, e1, e2, e3));
176 }
177 }
178 else if (v[i+1][j] && v[i+1][j+1]) {
179 GtsEdge * e1 = new_edge (v[i][j], v[i+1][j]);
180 GtsEdge * e2 = new_edge (v[i+1][j], v[i+1][j+1]);
181 GtsEdge * e3 =
182 gts_edge_new (GTS_EDGE_CLASS (gts_constraint_class ()),
183 v[i+1][j+1], v[i][j]);
184 gts_surface_add_face (s, gts_face_new (s->face_class, e1, e2, e3));
185 }
186 }
187 else if (v[i][j+1] && v[i+1][j+1] && v[i+1][j]) {
188 GtsEdge * e1 = new_edge (v[i+1][j], v[i+1][j+1]);
189 GtsEdge * e2 = new_edge (v[i+1][j+1], v[i][j+1]);
190 GtsEdge * e3 =
191 gts_edge_new (GTS_EDGE_CLASS (gts_constraint_class ()),
192 v[i][j+1], v[i+1][j]);
193 gts_surface_add_face (s, gts_face_new (s->face_class, e1, e2, e3));
194 }
195 }
196
triangle_is_hole(GtsTriangle * t)197 static gboolean triangle_is_hole (GtsTriangle * t)
198 {
199 GtsEdge * e1, * e2, * e3;
200 GtsVertex * v1, * v2, * v3;
201
202 gts_triangle_vertices_edges (t, NULL, &v1, &v2, &v3, &e1, &e2, &e3);
203
204 if ((GTS_IS_CONSTRAINT (e1) && GTS_SEGMENT (e1)->v2 != v1) ||
205 (GTS_IS_CONSTRAINT (e2) && GTS_SEGMENT (e2)->v2 != v2) ||
206 (GTS_IS_CONSTRAINT (e3) && GTS_SEGMENT (e3)->v2 != v3))
207 return TRUE;
208 return FALSE;
209 }
210
mark_as_hole(GtsFace * f,GtsSurface * s)211 static void mark_as_hole (GtsFace * f, GtsSurface * s)
212 {
213 GtsEdge * e1, * e2, * e3;
214
215 if (GTS_OBJECT (f)->reserved == f)
216 return;
217
218 GTS_OBJECT (f)->reserved = f;
219 e1 = GTS_TRIANGLE (f)->e1;
220 e2 = GTS_TRIANGLE (f)->e2;
221 e3 = GTS_TRIANGLE (f)->e3;
222
223 if (!GTS_IS_CONSTRAINT (e1)) {
224 GSList * i = e1->triangles;
225 while (i) {
226 GtsFace * f1 = i->data;
227 if (f1 != f && GTS_IS_FACE (f1) && gts_face_has_parent_surface (f1, s))
228 mark_as_hole (f1, s);
229 i = i->next;
230 }
231 }
232 if (!GTS_IS_CONSTRAINT (e2)) {
233 GSList * i = e2->triangles;
234 while (i) {
235 GtsFace * f1 = i->data;
236 if (f1 != f && GTS_IS_FACE (f1) && gts_face_has_parent_surface (f1, s))
237 mark_as_hole (f1, s);
238 i = i->next;
239 }
240 }
241 if (!GTS_IS_CONSTRAINT (e3)) {
242 GSList * i = e3->triangles;
243 while (i) {
244 GtsFace * f1 = i->data;
245 if (f1 != f && GTS_IS_FACE (f1) && gts_face_has_parent_surface (f1, s))
246 mark_as_hole (f1, s);
247 i = i->next;
248 }
249 }
250 }
251
edge_mark_as_hole(GtsEdge * e,GtsSurface * s)252 static void edge_mark_as_hole (GtsEdge * e, GtsSurface * s)
253 {
254 GSList * i = e->triangles;
255
256 while (i) {
257 GtsFace * f = i->data;
258 if (GTS_IS_FACE (f) &&
259 gts_face_has_parent_surface (f, s) &&
260 triangle_is_hole (GTS_TRIANGLE (f)))
261 mark_as_hole (f, s);
262 i = i->next;
263 }
264 }
265
face_is_marked(GtsObject * o)266 static gboolean face_is_marked (GtsObject * o)
267 {
268 if (o->reserved == o || gts_triangle_is_duplicate (GTS_TRIANGLE (o)))
269 return TRUE;
270 return FALSE;
271 }
272
build_constraint_list(GtsEdge * e,gpointer * data)273 static void build_constraint_list (GtsEdge * e, gpointer * data)
274 {
275 GtsSurface * s = data[0];
276 GSList ** constraints = data[1];
277
278 if (gts_edge_is_boundary (e, s))
279 *constraints = g_slist_prepend (*constraints, e);
280 }
281
cartesian_grid_triangulate_holes(CartesianGrid * grid,GtsSurface * s)282 static GtsSurface * cartesian_grid_triangulate_holes (CartesianGrid * grid,
283 GtsSurface * s)
284 {
285 GtsVertex * v1, * v2, * v3, * v4;
286 GtsEdge * e1, * e2, * e3, * e4, * e5;
287 gdouble w, h;
288 GtsSurface * box;
289 GSList * constraints = NULL, * vertices = NULL, * i;
290 gpointer data[2];
291
292 g_return_val_if_fail (grid != NULL, NULL);
293 g_return_val_if_fail (s != NULL, NULL);
294
295 /* build enclosing box */
296 w = grid->xmax - grid->xmin;
297 h = grid->ymax - grid->ymin;
298 v1 = gts_vertex_new (s->vertex_class, grid->xmin - w, grid->ymin - h, 0.);
299 v2 = gts_vertex_new (s->vertex_class, grid->xmax + w, grid->ymin - h, 0.);
300 v3 = gts_vertex_new (s->vertex_class, grid->xmax + w, grid->ymax + h, 0.);
301 v4 = gts_vertex_new (s->vertex_class, grid->xmin - w, grid->ymax + h, 0.);
302
303 e1 = gts_edge_new (s->edge_class, v1, v2);
304 e2 = gts_edge_new (s->edge_class, v2, v3);
305 e3 = gts_edge_new (s->edge_class, v3, v4);
306 e4 = gts_edge_new (s->edge_class, v4, v1);
307 e5 = gts_edge_new (s->edge_class, v1, v3);
308
309 box = gts_surface_new (GTS_SURFACE_CLASS (GTS_OBJECT (s)->klass),
310 s->face_class,
311 s->edge_class,
312 s->vertex_class);
313 gts_surface_add_face (box, gts_face_new (s->face_class, e1, e2, e5));
314 gts_surface_add_face (box, gts_face_new (s->face_class, e3, e4, e5));
315
316 /* build vertex and constraint list from the boundaries of the input
317 surface s */
318 data[0] = s;
319 data[1] = &constraints;
320 gts_surface_foreach_edge (s, (GtsFunc) build_constraint_list, data);
321 vertices = gts_vertices_from_segments (constraints);
322
323 /* triangulate holes */
324 i = vertices;
325 while (i) {
326 g_assert (!gts_delaunay_add_vertex (box, i->data, NULL));
327 i = i->next;
328 }
329 g_slist_free (vertices);
330
331 i = constraints;
332 while (i) {
333 g_assert (!gts_delaunay_add_constraint (box, i->data));
334 i = i->next;
335 }
336
337 /* destroy corners of the enclosing box */
338 gts_allow_floating_vertices = TRUE;
339 gts_object_destroy (GTS_OBJECT (v1));
340 gts_object_destroy (GTS_OBJECT (v2));
341 gts_object_destroy (GTS_OBJECT (v3));
342 gts_object_destroy (GTS_OBJECT (v4));
343 gts_allow_floating_vertices = FALSE;
344
345 /* remove parts of the mesh which are not holes */
346 i = constraints;
347 while (i) {
348 edge_mark_as_hole (i->data, box);
349 i = i->next;
350 }
351 g_slist_free (constraints);
352
353 /* remove marked and duplicate faces */
354 gts_surface_foreach_face_remove (box, (GtsFunc) face_is_marked, NULL);
355
356 /* box now contains only the triangulated holes. */
357 return box;
358 }
359
main(int argc,char * argv[])360 int main (int argc, char * argv[])
361 {
362 gboolean verbose = FALSE;
363 gboolean triangulate_holes = TRUE;
364 gboolean write_holes = FALSE;
365 int c = 0;
366 CartesianGrid * grid;
367 guint line = 0;
368 GtsSurface * s;
369 GTimer * timer;
370
371 /* parse options using getopt */
372 while (c != EOF) {
373 #ifdef HAVE_GETOPT_LONG
374 static struct option long_options[] = {
375 {"holes", no_argument, NULL, 'H'},
376 {"keep", no_argument, NULL, 'k'},
377 {"help", no_argument, NULL, 'h'},
378 {"verbose", no_argument, NULL, 'v'}
379 };
380 int option_index = 0;
381 switch ((c = getopt_long (argc, argv, "hvHk",
382 long_options, &option_index))) {
383 #else /* not HAVE_GETOPT_LONG */
384 switch ((c = getopt (argc, argv, "hvHk"))) {
385 #endif /* not HAVE_GETOPT_LONG */
386 case 'H': /* holes */
387 write_holes = TRUE;
388 break;
389 case 'k': /* keep */
390 triangulate_holes = FALSE;
391 break;
392 case 'v': /* verbose */
393 verbose = TRUE;
394 break;
395 case 'h': /* help */
396 fprintf (stderr,
397 "Usage: cartesian [OPTION] < FILE\n"
398 "Triangulates vertices of a regular cartesian mesh\n"
399 "(possibly containing holes)\n"
400 "\n"
401 " -k --keep keep holes\n"
402 " -H --holes write holes only\n"
403 " -v --verbose print statistics about the surface\n"
404 " -h --help display this help and exit\n"
405 "\n"
406 "Reports bugs to %s\n",
407 GTS_MAINTAINER);
408 return 0; /* success */
409 break;
410 case '?': /* wrong options */
411 fprintf (stderr, "Try `cartesian --help' for more information.\n");
412 return 1; /* failure */
413 }
414 }
415
416 grid = cartesian_grid_read (gts_vertex_class (), &line);
417 if (grid == NULL) {
418 fprintf (stderr,
419 "cartesian: file on standard input is not a valid cartesian grid\n"
420 "error at line %u\n", line);
421 return 1;
422 }
423
424 timer = g_timer_new ();
425 g_timer_start (timer);
426 s = gts_surface_new (gts_surface_class (),
427 gts_face_class (),
428 gts_edge_class (),
429 gts_vertex_class ());
430 cartesian_grid_triangulate (grid, s);
431 if (triangulate_holes) {
432 GtsSurface * holes = cartesian_grid_triangulate_holes (grid, s);
433 if (write_holes) {
434 gts_object_destroy (GTS_OBJECT (s));
435 s = holes;
436 }
437 else
438 gts_surface_merge (s, holes);
439 }
440 g_timer_stop (timer);
441
442 if (verbose) {
443 gts_surface_print_stats (s, stderr);
444 fprintf (stderr, "# Triangulation time: %g s speed: %.0f vertex/s\n",
445 g_timer_elapsed (timer, NULL),
446 gts_surface_vertex_number (s)/g_timer_elapsed (timer, NULL));
447 }
448
449 gts_surface_write (s, stdout);
450
451 return 0; /* success */
452 }
453