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 <math.h>
22 #include "config.h"
23 #ifdef HAVE_GETOPT_H
24 # include <getopt.h>
25 #endif /* HAVE_GETOPT_H */
26 #ifdef HAVE_UNISTD_H
27 # include <unistd.h>
28 #endif /* HAVE_UNISTD_H */
29 #include "gts.h"
30
31 #define HEAP_INSERT_EDGE(h, e) (GTS_OBJECT (e)->reserved = gts_eheap_insert (h, e))
32 #define HEAP_REMOVE_EDGE(h, e) (gts_eheap_remove (h, GTS_OBJECT (e)->reserved),\
33 GTS_OBJECT (e)->reserved = NULL)
34
triangles_angle(GtsPoint * p1,GtsPoint * p2,GtsPoint * p3,GtsPoint * p4)35 static gdouble triangles_angle (GtsPoint * p1, GtsPoint * p2,
36 GtsPoint * p3, GtsPoint * p4)
37 {
38 gdouble x1, y1, z1, x2, y2, z2;
39 gdouble nx1, ny1, nz1, nx2, ny2, nz2;
40 gdouble pvx, pvy, pvz;
41 gdouble theta;
42
43 x1 = p2->x - p1->x;
44 y1 = p2->y - p1->y;
45 z1 = p2->z - p1->z;
46
47 x2 = p3->x - p1->x;
48 y2 = p3->y - p1->y;
49 z2 = p3->z - p1->z;
50
51 nx1 = y1*z2 - z1*y2;
52 ny1 = z1*x2 - x1*z2;
53 nz1 = x1*y2 - y1*x2;
54
55 x1 = p1->x - p2->x;
56 y1 = p1->y - p2->y;
57 z1 = p1->z - p2->z;
58
59 x2 = p4->x - p2->x;
60 y2 = p4->y - p2->y;
61 z2 = p4->z - p2->z;
62
63 nx2 = y1*z2 - z1*y2;
64 ny2 = z1*x2 - x1*z2;
65 nz2 = x1*y2 - y1*x2;
66
67 pvx = ny1*nz2 - nz1*ny2;
68 pvy = nz1*nx2 - nx1*nz2;
69 pvz = nx1*ny2 - ny1*nx2;
70
71 theta = atan2 (sqrt (pvx*pvx + pvy*pvy + pvz*pvz),
72 nx1*nx2 + ny1*ny2 + nz1*nz2) - M_PI;
73 return theta < - M_PI ? theta + 2.*M_PI : theta;
74
75 }
76
edge_swap_cost(GtsEdge * e)77 static gdouble edge_swap_cost (GtsEdge * e)
78 {
79 GSList * i;
80 GtsTriangle * t1 = NULL, * t2 = NULL;
81 GtsVertex * v1, * v2, * v3, * v4;
82 GtsEdge * e1, * e2, * e3, * e4;
83 gdouble ab, aa;
84
85 i = e->triangles;
86 while (i) {
87 if (GTS_IS_FACE (i->data)) {
88 if (!t1) t1 = i->data;
89 else if (!t2) t2 = i->data;
90 else return G_MAXDOUBLE;
91 }
92 i = i->next;
93 }
94 if (!t1 || !t2)
95 return G_MAXDOUBLE;
96
97 gts_triangle_vertices_edges (t1, e, &v1, &v2, &v3, &e, &e3, &e4);
98 gts_triangle_vertices_edges (t2, e, &v2, &v1, &v4, &e, &e1, &e2);
99
100 ab = triangles_angle (GTS_POINT (v1), GTS_POINT (v2),
101 GTS_POINT (v3), GTS_POINT (v4));
102 aa = triangles_angle (GTS_POINT (v3), GTS_POINT (v4),
103 GTS_POINT (v2), GTS_POINT (v1));
104 return fabs (ab) - fabs (aa);
105 }
106
edge_swap(GtsEdge * e,GtsSurface * s,GtsEHeap * heap)107 static void edge_swap (GtsEdge * e, GtsSurface * s, GtsEHeap * heap)
108 {
109 GSList * i;
110 GtsTriangle * t1 = NULL, * t2 = NULL;
111 GtsVertex * v1, * v2, * v3, * v4;
112 GtsEdge * e1, * e2, * e3, * e4;
113
114 i = e->triangles;
115 while (i) {
116 if (GTS_IS_FACE (i->data)) {
117 if (!t1) t1 = i->data;
118 else if (!t2) t2 = i->data;
119 else g_assert_not_reached ();
120 }
121 i = i->next;
122 }
123 g_assert (t1 && t2);
124
125 gts_triangle_vertices_edges (t1, e, &v1, &v2, &v3, &e, &e3, &e4);
126 gts_triangle_vertices_edges (t2, e, &v2, &v1, &v4, &e, &e1, &e2);
127
128 gts_object_destroy (GTS_OBJECT (e));
129 e = gts_edge_new (s->edge_class, v3, v4);
130 gts_surface_add_face (s, gts_face_new (s->face_class, e, e4, e1));
131 gts_surface_add_face (s, gts_face_new (s->face_class, e, e2, e3));
132
133 HEAP_INSERT_EDGE (heap, e);
134 HEAP_REMOVE_EDGE (heap, e1);
135 HEAP_INSERT_EDGE (heap, e1);
136 HEAP_REMOVE_EDGE (heap, e2);
137 HEAP_INSERT_EDGE (heap, e2);
138 HEAP_REMOVE_EDGE (heap, e3);
139 HEAP_INSERT_EDGE (heap, e3);
140 HEAP_REMOVE_EDGE (heap, e4);
141 HEAP_INSERT_EDGE (heap, e4);
142 }
143
create_heap_optimize(GtsEdge * e,GtsEHeap * heap)144 static void create_heap_optimize (GtsEdge * e, GtsEHeap * heap)
145 {
146 HEAP_INSERT_EDGE (heap, e);
147 }
148
surface_optimize(GtsSurface * surface,gdouble max_cost)149 static void surface_optimize (GtsSurface * surface,
150 gdouble max_cost)
151 {
152 GtsEHeap * heap;
153 GtsEdge * e;
154 gdouble top_cost;
155
156 heap = gts_eheap_new ((GtsKeyFunc) edge_swap_cost, NULL);
157 gts_eheap_freeze (heap);
158 gts_surface_foreach_edge (surface, (GtsFunc) create_heap_optimize, heap);
159 gts_eheap_thaw (heap);
160
161 gts_allow_floating_edges = TRUE;
162 while ((e = gts_eheap_remove_top (heap, &top_cost)) &&
163 top_cost < max_cost)
164 edge_swap (e, surface, heap);
165 gts_allow_floating_edges = FALSE;
166
167 if (e) GTS_OBJECT (e)->reserved = NULL;
168 gts_eheap_foreach (heap, (GFunc) gts_object_reset_reserved, NULL);
169
170 gts_eheap_destroy (heap);
171 }
172
angle_stats(GtsEdge * e,GtsRange * angle)173 static void angle_stats (GtsEdge * e, GtsRange * angle)
174 {
175 GSList * i;
176 GtsTriangle * t1 = NULL, * t2 = NULL;
177
178 i = e->triangles;
179 while (i) {
180 if (GTS_IS_FACE (i->data)) {
181 if (!t1) t1 = i->data;
182 else if (!t2) t2 = i->data;
183 else return;
184 }
185 i = i->next;
186 }
187 if (!t1 || !t2)
188 return;
189
190 gts_range_add_value (angle, fabs (gts_triangles_angle (t1, t2)));
191 }
192
main(int argc,char * argv[])193 int main (int argc, char * argv[])
194 {
195 GtsSurface * s;
196 gboolean verbose = FALSE;
197 gdouble threshold;
198 int c = 0;
199 GtsFile * fp;
200 GtsRange angle;
201
202 /* parse options using getopt */
203 while (c != EOF) {
204 #ifdef HAVE_GETOPT_LONG
205 static struct option long_options[] = {
206 {"help", no_argument, NULL, 'h'},
207 {"verbose", no_argument, NULL, 'v'}
208 };
209 int option_index = 0;
210 switch ((c = getopt_long (argc, argv, "hv",
211 long_options, &option_index))) {
212 #else /* not HAVE_GETOPT_LONG */
213 switch ((c = getopt (argc, argv, "hv"))) {
214 #endif /* not HAVE_GETOPT_LONG */
215 case 'v': /* verbose */
216 verbose = TRUE;
217 break;
218 case 'h': /* help */
219 fprintf (stderr,
220 "Usage: optimize [OPTION] THRESHOLD < FILE\n"
221 "\n"
222 " -v --verbose print statistics about the surface\n"
223 " -h --help display this help and exit\n"
224 "\n"
225 "Report bugs to %s\n",
226 GTS_MAINTAINER);
227 return 0; /* success */
228 break;
229 case '?': /* wrong options */
230 fprintf (stderr, "Try `optimize --help' for more information.\n");
231 return 1; /* failure */
232 }
233 }
234
235 if (optind >= argc) { /* missing threshold */
236 fprintf (stderr,
237 "optimize: missing THRESHOLD\n"
238 "Try `optimize --help' for more information.\n");
239 return 1; /* failure */
240 }
241 threshold = atof (argv[optind]);
242
243 if (threshold < 0.0) { /* threshold must be positive */
244 fprintf (stderr,
245 "optimize: THRESHOLD must be >= 0.0\n"
246 "Try `optimize --help' for more information.\n");
247 return 1; /* failure */
248 }
249
250 /* read surface in */
251 s = gts_surface_new (gts_surface_class (),
252 gts_face_class (),
253 gts_edge_class (),
254 gts_vertex_class ());
255 fp = gts_file_new (stdin);
256 if (gts_surface_read (s, fp)) {
257 fputs ("optimize: file on standard input is not a valid GTS file\n",
258 stderr);
259 fprintf (stderr, "stdin:%d:%d: %s\n", fp->line, fp->pos, fp->error);
260 return 1; /* failure */
261 }
262
263 /* if verbose on print stats */
264 if (verbose) {
265 gts_surface_print_stats (s, stderr);
266 gts_range_init (&angle);
267 gts_surface_foreach_edge (s, (GtsFunc) angle_stats, &angle);
268 gts_range_update (&angle);
269 fputs ("# angle : ", stderr);
270 gts_range_print (&angle, stderr);
271 fputc ('\n', stderr);
272 }
273
274 surface_optimize (s, -threshold);
275
276 /* if verbose on print stats */
277 if (verbose) {
278 gts_surface_print_stats (s, stderr);
279 gts_range_init (&angle);
280 gts_surface_foreach_edge (s, (GtsFunc) angle_stats, &angle);
281 gts_range_update (&angle);
282 fputs ("# angle : ", stderr);
283 gts_range_print (&angle, stderr);
284 fputc ('\n', stderr);
285 }
286
287 /* write surface */
288 gts_surface_write (s, stdout);
289
290 return 0; /* success */
291 }
292