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 <math.h>
21 #include <stdlib.h>
22 #include <locale.h>
23 #include <string.h>
24 #include <pgm.h>
25 #include "config.h"
26 #ifdef HAVE_GETOPT_H
27 # include <getopt.h>
28 #endif /* HAVE_GETOPT_H */
29 #ifdef HAVE_UNISTD_H
30 # include <unistd.h>
31 #endif /* HAVE_UNISTD_H */
32 #include "gts.h"
33
34 typedef enum { NUMBER, COST } StopOptions;
35
stop_number(gdouble cost,guint number,guint * max)36 static gboolean stop_number (gdouble cost, guint number, guint * max)
37 {
38 if (number >= *max)
39 return TRUE;
40 return FALSE;
41 }
42
stop_number_log_cost(gdouble cost,guint number,guint * max)43 static gboolean stop_number_log_cost (gdouble cost, guint number, guint * max)
44 {
45 if (number >= *max)
46 return TRUE;
47 fprintf (stderr, "%d %g\n", number, cost);
48 return FALSE;
49 }
50
stop_number_verbose(gdouble cost,guint number,guint * max)51 static gboolean stop_number_verbose (gdouble cost, guint number, guint * max)
52 {
53 static guint nmax = 0, nold = 0;
54 static GTimer * timer = NULL, * total_timer = NULL;
55
56 if (timer == NULL) {
57 nmax = nold = number;
58 timer = g_timer_new ();
59 total_timer = g_timer_new ();
60 g_timer_start (total_timer);
61 }
62
63 if (number != nold && number % 121 == 0 &&
64 number > nmax && nmax <= *max) {
65 gdouble total_elapsed = g_timer_elapsed (total_timer, NULL);
66 gdouble remaining;
67 gdouble hours, mins, secs;
68 gdouble hours1, mins1, secs1;
69
70 g_timer_stop (timer);
71
72 hours = floor (total_elapsed/3600.);
73 mins = floor ((total_elapsed - 3600.*hours)/60.);
74 secs = floor (total_elapsed - 3600.*hours - 60.*mins);
75
76 remaining = total_elapsed*((*max - nmax)/(gdouble) (number - nmax) - 1.);
77 hours1 = floor (remaining/3600.);
78 mins1 = floor ((remaining - 3600.*hours1)/60.);
79 secs1 = floor (remaining - 3600.*hours1 - 60.*mins1);
80
81 fprintf (stderr,
82 "\rVertices: %8u %3.0f%% %6.0f vertices/s "
83 "Elapsed: %02.0f:%02.0f:%02.0f "
84 "Remaining: %02.0f:%02.0f:%02.0f ",
85 number,
86 100.*(number - nmax)/(*max - nmax),
87 (number - nold)/g_timer_elapsed (timer, NULL),
88 hours, mins, secs,
89 hours1, mins1, secs1);
90 fflush (stderr);
91
92 nold = number;
93 g_timer_start (timer);
94 }
95 if (number >= *max) {
96 g_timer_destroy (timer);
97 g_timer_destroy (total_timer);
98 return TRUE;
99 }
100 return FALSE;
101 }
102
stop_cost(gdouble cost,guint number,gdouble * min)103 static gboolean stop_cost (gdouble cost, guint number, gdouble * min)
104 {
105 if (cost < *min)
106 return TRUE;
107 return FALSE;
108 }
109
stop_cost_verbose(gdouble cost,guint number,gdouble * min)110 static gboolean stop_cost_verbose (gdouble cost, guint number, gdouble * min)
111 {
112 if (number % 121 == 0) {
113 fprintf (stderr, "\rVertices: %10u Cost: %10.2f ", number, cost);
114 fflush (stderr);
115 }
116 if (cost < *min)
117 return TRUE;
118 return FALSE;
119 }
120
121 /* ListFace: Header */
122
123 typedef struct _ListFace ListFace;
124
125 struct _ListFace {
126 /*< private >*/
127 GtsListFace parent;
128
129 /*< public >*/
130 GtsEHeap * heap;
131 GtsEHeapPair * pair;
132 GtsVertex * best;
133 GtsVector p;
134 };
135
136 typedef struct _ListFaceClass ListFaceClass;
137
138 struct _ListFaceClass {
139 /*< private >*/
140 GtsFaceClass parent_class;
141
142 /*< public >*/
143 void (* cost_init) (ListFace *);
144 };
145
146 #define LIST_FACE(obj) GTS_OBJECT_CAST (obj,\
147 ListFace,\
148 list_face_class ())
149 #define LIST_FACE_CLASS(klass) GTS_OBJECT_CLASS_CAST (klass,\
150 ListFaceClass,\
151 list_face_class())
152 #define IS_LIST_FACE(obj) (gts_object_is_from_class (obj,\
153 list_face_class ()))
154
155 static ListFaceClass * list_face_class (void);
156
157 /* ListFace: Object */
158
list_face_destroy(GtsObject * object)159 static void list_face_destroy (GtsObject * object)
160 {
161 ListFace * f = LIST_FACE (object);
162
163 if (f->heap) {
164 gts_eheap_remove (f->heap, f->pair);
165 f->heap = NULL;
166 }
167
168 (* GTS_OBJECT_CLASS (list_face_class ())->parent_class->destroy)
169 (object);
170 }
171
triangle_plane(GtsTriangle * f,GtsVector p)172 static void triangle_plane (GtsTriangle * f, GtsVector p)
173 {
174 GtsPoint * v1, * v2, * v3;
175 gdouble x1, x2, y1, y2, det;
176
177 v1 = GTS_POINT (GTS_SEGMENT (f->e1)->v1);
178 v2 = GTS_POINT (GTS_SEGMENT (f->e1)->v2);
179 v3 = GTS_POINT (gts_triangle_vertex (f));
180
181 x1 = v2->x - v1->x;
182 y1 = v2->y - v1->y;
183 x2 = v3->x - v1->x;
184 y2 = v3->y - v1->y;
185 det = x1*y2 - x2*y1;
186 g_assert (det != 0.);
187
188 p[0] = (y2*(v2->z - v1->z) - y1*(v3->z - v1->z))/det;
189 p[1] = (-x2*(v2->z - v1->z) + x1*(v3->z - v1->z))/det;
190 p[2] = ((- v1->x*y2 + v1->y*x2)*(v2->z - v1->z) +
191 (- v1->y*x1 + v1->x*y1)*(v3->z - v1->z))/det + v1->z;
192 }
193
list_face_cost_init(ListFace * f)194 static void list_face_cost_init (ListFace * f)
195 {
196 triangle_plane (GTS_TRIANGLE (f), f->p);
197 }
198
list_face_class_init(ListFaceClass * klass)199 static void list_face_class_init (ListFaceClass * klass)
200 {
201 GTS_OBJECT_CLASS (klass)->destroy = list_face_destroy;
202 klass->cost_init = list_face_cost_init;
203 }
204
list_face_class(void)205 static ListFaceClass * list_face_class (void)
206 {
207 static ListFaceClass * klass = NULL;
208
209 if (klass == NULL) {
210 GtsObjectClassInfo list_face_info = {
211 "ListFace",
212 sizeof (ListFace),
213 sizeof (ListFaceClass),
214 (GtsObjectClassInitFunc) list_face_class_init,
215 (GtsObjectInitFunc) NULL,
216 (GtsArgSetFunc) NULL,
217 (GtsArgGetFunc) NULL
218 };
219 klass = gts_object_class_new (GTS_OBJECT_CLASS (gts_list_face_class ()),
220 &list_face_info);
221 }
222
223 return klass;
224 }
225
226 typedef gdouble (*CostFunc) (ListFace *, GtsPoint *, gpointer);
227
list_face_update(ListFace * f,gpointer * data)228 static void list_face_update (ListFace * f, gpointer * data)
229 {
230 if (IS_LIST_FACE (f)) {
231 GSList * i = GTS_LIST_FACE (f)->points;
232 GtsEHeap * heap = data[0];
233 CostFunc cost = data[1];
234 gpointer cost_data = data[2];
235
236 if (i) {
237 gdouble maxerr = 0.;
238 ListFaceClass * klass = LIST_FACE_CLASS (GTS_OBJECT (f)->klass);
239
240 (* klass->cost_init) (f);
241 f->best = NULL;
242 while (i) {
243 GtsPoint * v = i->data;
244 gdouble e = (* cost) (f, v, cost_data);
245
246 if (e > maxerr) {
247 maxerr = e;
248 f->best = GTS_VERTEX (v);
249 }
250 i = i->next;
251 }
252 if (f->heap) {
253 gts_eheap_remove (f->heap, f->pair);
254 f->heap = NULL;
255 }
256 if (maxerr > 0.) {
257 f->pair = gts_eheap_insert_with_key (heap, f, - maxerr);
258 f->heap = heap;
259 }
260 }
261 }
262 }
263
list_face_clear_heap(ListFace * f)264 static void list_face_clear_heap (ListFace * f)
265 {
266 f->heap = NULL;
267 }
268
surface_hf_refine(GtsSurface * s,CostFunc cost_func,gpointer cost_data,GtsStopFunc stop_func,gpointer stop_data)269 static void surface_hf_refine (GtsSurface * s,
270 CostFunc cost_func,
271 gpointer cost_data,
272 GtsStopFunc stop_func,
273 gpointer stop_data)
274 {
275 GtsEHeap * heap;
276 gdouble top_cost;
277 guint nv = 4;
278 GtsListFace * f;
279 gpointer data[3];
280
281 g_return_if_fail (s != NULL);
282 g_return_if_fail (cost_func != NULL);
283 g_return_if_fail (stop_func != NULL);
284
285 data[0] = heap = gts_eheap_new (NULL, NULL);
286 data[1] = cost_func;
287 data[2] = cost_data;
288 gts_surface_foreach_face (s, (GtsFunc) list_face_update, data);
289 while ((f = gts_eheap_remove_top (heap, &top_cost)) &&
290 !(*stop_func) (- top_cost, nv, stop_data)) {
291 GtsVertex * v = LIST_FACE (f)->best;
292 GSList * t, * i;
293
294 LIST_FACE (f)->heap = NULL;
295 gts_delaunay_add_vertex_to_face (s, v, GTS_FACE (f));
296 i = t = gts_vertex_triangles (v, NULL);
297 while (i) {
298 list_face_update (i->data, data);
299 i = i->next;
300 }
301 g_slist_free (t);
302 nv++;
303 }
304
305 if (f)
306 LIST_FACE (f)->heap = NULL;
307
308 gts_eheap_foreach (heap, (GFunc) list_face_clear_heap, NULL);
309 gts_eheap_destroy (heap);
310 }
311
prepend(GtsListFace * f,gray ** g,guint i,guint j)312 static void prepend (GtsListFace * f,
313 gray ** g, guint i, guint j)
314 {
315 f->points = g_slist_prepend (f->points,
316 gts_vertex_new (gts_vertex_class (),
317 i, j, g[j][i]));
318 }
319
destroy_unused(GtsListFace * f)320 static void destroy_unused (GtsListFace * f)
321 {
322 g_slist_foreach (f->points, (GFunc) gts_object_destroy, NULL);
323 g_slist_free (f->points);
324 f->points = NULL;
325 }
326
happrox(gray ** g,gint width,gint height,CostFunc cost_func,gpointer cost_data,GtsStopFunc stop_func,gpointer stop_data)327 static GtsSurface * happrox (gray ** g, gint width, gint height,
328 CostFunc cost_func,
329 gpointer cost_data,
330 GtsStopFunc stop_func,
331 gpointer stop_data)
332 {
333 GtsSurface * s = gts_surface_new (gts_surface_class (),
334 GTS_FACE_CLASS (list_face_class ()),
335 gts_edge_class (),
336 gts_vertex_class ());
337 GtsVertex * v1 = gts_vertex_new (s->vertex_class, 0., 0., g[0][0]);
338 GtsVertex * v2 = gts_vertex_new (s->vertex_class,
339 0., height - 1, g[height - 1][0]);
340 GtsVertex * v3 = gts_vertex_new (s->vertex_class,
341 width - 1, 0., g[0][width - 1]);
342 GtsVertex * v4 = gts_vertex_new (s->vertex_class,
343 width - 1, height - 1,
344 g[height - 1][width - 1]);
345 guint i, j;
346 GSList * corners = NULL;
347 GtsTriangle * t;
348 GtsVertex * w1, * w2, * w3;
349 GtsListFace * f;
350
351 /* creates enclosing triangle */
352 corners = g_slist_prepend (corners, v1);
353 corners = g_slist_prepend (corners, v2);
354 corners = g_slist_prepend (corners, v3);
355 corners = g_slist_prepend (corners, v4);
356 t = gts_triangle_enclosing (gts_triangle_class (), corners, 100.);
357 g_slist_free (corners);
358 gts_triangle_vertices (t, &w1, &w2, &w3);
359
360 f = GTS_LIST_FACE (gts_face_new (s->face_class, t->e1, t->e2, t->e3));
361 gts_surface_add_face (s, GTS_FACE (f));
362
363 /* add PGM vertices (corners excepted) to point list of f */
364 for (i = 1; i < width - 1; i++) {
365 for (j = 1; j < height - 1; j++)
366 prepend (f, g, i, j);
367 prepend (f, g, i, 0);
368 prepend (f, g, i, height - 1);
369 }
370 for (j = 1; j < height - 1; j++) {
371 prepend (f, g, 0, j);
372 prepend (f, g, width - 1, j);
373 }
374 pgm_freearray (g, height);
375
376 /* add four corners to initial triangulation */
377 g_assert (gts_delaunay_add_vertex_to_face (s, v1, GTS_FACE (f)) == NULL);
378 f = GTS_LIST_FACE (gts_point_locate (GTS_POINT (v2), s, NULL));
379 g_assert (gts_delaunay_add_vertex_to_face (s, v2, GTS_FACE (f)) == NULL);
380 f = GTS_LIST_FACE (gts_point_locate (GTS_POINT (v3), s, NULL));
381 g_assert (gts_delaunay_add_vertex_to_face (s, v3, GTS_FACE (f)) == NULL);
382 f = GTS_LIST_FACE (gts_point_locate (GTS_POINT (v4), s, NULL));
383 g_assert (gts_delaunay_add_vertex_to_face (s, v4, GTS_FACE (f)) == NULL);
384
385 /* refine surface */
386 surface_hf_refine (s, cost_func, cost_data, stop_func, stop_data);
387
388 /* destroy unused vertices */
389 gts_surface_foreach_face (s, (GtsFunc) destroy_unused, NULL);
390
391 /* destroy enclosing triangle */
392 gts_allow_floating_vertices = TRUE;
393 gts_object_destroy (GTS_OBJECT (w1));
394 gts_object_destroy (GTS_OBJECT (w2));
395 gts_object_destroy (GTS_OBJECT (w3));
396 gts_allow_floating_vertices = FALSE;
397
398 return s;
399 }
400
happrox_list(GSList * points,gboolean keep_enclosing,gboolean closed,CostFunc cost_func,gpointer cost_data,GtsStopFunc stop_func,gpointer stop_data)401 static GtsSurface * happrox_list (GSList * points,
402 gboolean keep_enclosing,
403 gboolean closed,
404 CostFunc cost_func,
405 gpointer cost_data,
406 GtsStopFunc stop_func,
407 gpointer stop_data)
408 {
409 GtsSurface * s = gts_surface_new (gts_surface_class (),
410 GTS_FACE_CLASS (list_face_class ()),
411 gts_edge_class (),
412 gts_vertex_class ());
413 GtsTriangle * t;
414 GtsVertex * w1, * w2, * w3;
415 GtsListFace * f;
416
417 /* creates enclosing triangle */
418 t = gts_triangle_enclosing (gts_triangle_class (), points, 10.);
419 gts_triangle_vertices (t, &w1, &w2, &w3);
420 GTS_POINT (w1)->z = GTS_POINT (w2)->z = GTS_POINT (w3)->z =
421 keep_enclosing ? -10. : -1e30;
422
423 f = GTS_LIST_FACE (gts_face_new (s->face_class, t->e1, t->e2, t->e3));
424 gts_surface_add_face (s, GTS_FACE (f));
425 f->points = points;
426
427 /* refine surface */
428 surface_hf_refine (s, cost_func, cost_data, stop_func, stop_data);
429
430 /* destroy unused vertices */
431 gts_surface_foreach_face (s, (GtsFunc) destroy_unused, NULL);
432
433 /* destroy enclosing triangle */
434 if (!keep_enclosing) {
435 gts_allow_floating_vertices = TRUE;
436 gts_object_destroy (GTS_OBJECT (w1));
437 gts_object_destroy (GTS_OBJECT (w2));
438 gts_object_destroy (GTS_OBJECT (w3));
439 gts_allow_floating_vertices = FALSE;
440 }
441 else if (closed) {
442 GSList * l = gts_surface_boundary (s);
443 GtsFace * f;
444
445 g_assert (g_slist_length (l) == 3);
446 f = gts_face_new (s->face_class, l->data, l->next->data, l->next->next->data);
447 gts_surface_add_face (s, f);
448 if (!gts_face_is_compatible (f, s))
449 gts_triangle_revert (GTS_TRIANGLE (f));
450 g_slist_free (l);
451 gts_object_destroy (GTS_OBJECT (t));
452 }
453 else
454 gts_object_destroy (GTS_OBJECT (t));
455
456 return s;
457 }
458
height_cost(ListFace * f,GtsPoint * p)459 static gdouble height_cost (ListFace * f, GtsPoint * p)
460 {
461 return fabs (p->x*f->p[0] + p->y*f->p[1] + f->p[2] - p->z);
462 }
463
relative_height_cost(ListFace * f,GtsPoint * p,gdouble * z)464 static gdouble relative_height_cost (ListFace * f, GtsPoint * p, gdouble * z)
465 {
466 gdouble dz = p->x*f->p[0] + p->y*f->p[1] + f->p[2] - p->z;
467
468 return fabs (p->z) < *z ? fabs (dz/(*z)) : fabs (dz/p->z);
469 }
470
main(int argc,char * argv[])471 int main (int argc, char * argv[])
472 {
473 int c = 0;
474 gboolean verbose = FALSE;
475 GtsSurface * s;
476 CostFunc cost_func = (CostFunc) height_cost;
477 gpointer cost_data = NULL;
478 GtsStopFunc stop_func = NULL;
479 gpointer stop_data = NULL;
480 StopOptions stop = NUMBER;
481 guint number = 0;
482 gdouble cmin = 0.0;
483 gboolean logcost = FALSE;
484 gboolean flat = FALSE;
485 gdouble relative = 0.;
486 gboolean keep_enclosing = FALSE;
487 gboolean closed = FALSE;
488
489 if (!setlocale (LC_ALL, "POSIX"))
490 g_warning ("cannot set locale to POSIX");
491
492 while (c != EOF) {
493 #ifdef HAVE_GETOPT_LONG
494 static struct option long_options[] = {
495 {"closed", no_argument, NULL, 'C'},
496 {"keep", no_argument, NULL, 'k'},
497 {"relative", required_argument, NULL, 'r'},
498 {"flat", no_argument, NULL, 'f'},
499 {"log", no_argument, NULL, 'l'},
500 {"number", required_argument, NULL, 'n'},
501 {"cost", required_argument, NULL, 'c'},
502 {"help", no_argument, NULL, 'h'},
503 {"verbose", no_argument, NULL, 'v'}
504 };
505 int option_index = 0;
506 switch ((c = getopt_long (argc, argv, "hvn:c:lfr:kC",
507 long_options, &option_index))) {
508 #else /* not HAVE_GETOPT_LONG */
509 switch ((c = getopt (argc, argv, "hvn:c:lfr:kC"))) {
510 #endif /* not HAVE_GETOPT_LONG */
511 case 'C': /* closed */
512 closed = TRUE;
513 keep_enclosing = TRUE;
514 break;
515 case 'k': /* keep */
516 keep_enclosing = TRUE;
517 break;
518 case 'r': /* relative */
519 cost_func = (CostFunc) relative_height_cost;
520 relative = atof (optarg);
521 cost_data = &relative;
522 if (relative <= 0.) {
523 fprintf (stderr, "happrox: argument for -r option must be strictly positive\n");
524 return 1;
525 }
526 break;
527 case 'f': /* flat file */
528 flat = TRUE;
529 break;
530 case 'l': /* log cost */
531 logcost = TRUE;
532 break;
533 case 'n': /* stop by number */
534 stop = NUMBER;
535 number = atoi (optarg);
536 break;
537 case 'c': /* stop by cost */
538 stop = COST;
539 cmin = atof (optarg);
540 break;
541 case 'h': /* help */
542 fprintf (stderr,
543 "Usage: happrox [OPTION]... < [input.pgm|input] > output.gts\n"
544 "Returns a simplified triangulation of a set of points using\n"
545 "algorithm III of Garland and Heckbert (1995).\n"
546 "\n"
547 " -n N, --number=N stop the refinement process if the number of\n"
548 " vertices is larger than N\n"
549 " -c C, --cost=C stop the refinement process if the cost of insertion\n"
550 " of a vertex is smaller than C\n"
551 " -f --flat input is a flat file with three x,y,z columns\n"
552 " (default is PGM file)\n"
553 " -r Z --relative=Z use relative height cost for all heights larger than Z\n"
554 " -k --keep keep enclosing triangle\n"
555 " -C --closed close the surface\n"
556 " -l --log log evolution of cost\n"
557 " -v, --verbose display surface statistics\n"
558 " -h, --help display this help and exit\n"
559 "\n"
560 "Report bugs to %s\n",
561 GTS_MAINTAINER);
562 return 0;
563 break;
564 case 'v':
565 verbose = TRUE;
566 break;
567 case '?': /* wrong options */
568 fprintf (stderr, "Try `happrox --help' for more information.\n");
569 return 1;
570 }
571 }
572
573 switch (stop) {
574 case NUMBER:
575 if (verbose)
576 stop_func = (GtsStopFunc) stop_number_verbose;
577 else
578 stop_func = (GtsStopFunc) stop_number;
579 if (logcost)
580 stop_func = (GtsStopFunc) stop_number_log_cost;
581 stop_data = &number;
582 break;
583 case COST:
584 if (verbose)
585 stop_func = (GtsStopFunc) stop_cost_verbose;
586 else
587 stop_func = (GtsStopFunc) stop_cost;
588 stop_data = &cmin;
589 break;
590 default:
591 g_assert_not_reached ();
592 }
593
594 if (flat) {
595 GSList * points = NULL;
596 guint n = 0;
597 gdouble x, y, z;
598
599 while (scanf ("%lf %lf %lf", &x, &y, &z) == 3) {
600 points = g_slist_prepend (points, gts_vertex_new (gts_vertex_class (), x, y, z));
601 n++;
602 }
603 if (verbose)
604 fprintf (stderr, "happrox: %d vertices\n", n);
605
606 s = happrox_list (points, keep_enclosing, closed, cost_func, cost_data, stop_func, stop_data);
607 }
608 else {
609 gint width, height;
610 gray ** g, maxgray;
611
612 g = pgm_readpgm (stdin, &width, &height, &maxgray);
613 if (verbose)
614 fprintf (stderr, "happrox: width: %d height: %d maxgray: %d\n",
615 width, height, maxgray);
616
617 s = happrox (g, width, height, cost_func, cost_data, stop_func, stop_data);
618 }
619
620 if (verbose) {
621 fputc ('\n', stderr);
622 gts_surface_print_stats (s, stderr);
623 }
624 gts_surface_write (s, stdout);
625
626 return 0;
627 }
628