1 /* GTS-Library conform marching tetrahedra algorithm
2 * Copyright (C) 2002 Gert Wollny
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 <string.h>
22 #include <gts.h>
23 #ifdef NATIVE_WIN32
24 # include <memory.h>
25 # define M_SQRT2 1.41421356237309504880
26 #endif /* NATIVE_WIN32 */
27
28 typedef struct {
29 gint nx, ny;
30 gdouble ** data;
31 } slice_t;
32
33 typedef struct {
34 gint x, y, z;
35 gboolean mid;
36 gdouble d;
37 } tetra_vertex_t;
38
39 /* this helper is a lookup table for vertices */
40 typedef struct {
41 gint nx, ny;
42 GtsVertex ** vtop, ** vmid, **vbot;
43 } helper_t ;
44
45 typedef struct {
46 GHashTable * vbot, * vtop;
47 } helper_bcl ;
48
49
init_helper(int nx,int ny)50 static helper_t * init_helper (int nx, int ny)
51 {
52 gint nxy = 4*nx*ny;
53 helper_t *retval = g_malloc0 (sizeof (helper_t));
54
55 retval->nx = nx;
56 retval->ny = ny;
57 retval->vtop = g_malloc0 (sizeof (GtsVertex *)*nxy);
58 retval->vmid = g_malloc0 (sizeof (GtsVertex *)*nxy);
59 retval->vbot = g_malloc0 (sizeof (GtsVertex *)*nxy);
60 return retval;
61 }
62
init_helper_bcl(void)63 static helper_bcl * init_helper_bcl (void)
64 {
65 helper_bcl *retval = g_malloc0 (sizeof (helper_bcl));
66
67 retval->vtop = g_hash_table_new (g_str_hash, g_str_equal);
68 retval->vbot = g_hash_table_new (g_str_hash, g_str_equal);
69 return retval;
70 }
71
free_helper(helper_t * h)72 static void free_helper (helper_t * h)
73 {
74 g_free (h->vtop);
75 g_free (h->vmid);
76 g_free (h->vbot);
77 g_free (h);
78 }
79
free_helper_bcl(helper_bcl * h)80 static void free_helper_bcl (helper_bcl * h)
81 {
82 g_hash_table_destroy (h->vtop);
83 g_hash_table_destroy (h->vbot);
84 g_free (h);
85 }
86
87 /* move the vertices in the bottom slice to the top, and clear the
88 other slices in the lookup tables */
helper_advance(helper_t * h)89 static void helper_advance (helper_t * h)
90 {
91 GtsVertex ** help = h->vbot;
92 h->vbot = h->vtop;
93 h->vtop = help;
94
95 memset (h->vmid, 0, 4*sizeof(GtsVertex *) * h->nx * h->ny);
96 memset (h->vbot, 0, 4*sizeof(GtsVertex *) * h->nx * h->ny);
97 }
98
helper_advance_bcl(helper_bcl * h)99 static void helper_advance_bcl (helper_bcl * h)
100 {
101 GHashTable * help = g_hash_table_new (g_str_hash, g_str_equal);
102
103 g_hash_table_destroy (h->vbot);
104 h->vbot = h->vtop;
105 h->vtop = help;
106 }
107
108 /* find the zero-crossing of line through v1 and v2 and return the
109 corresponding GtsVertex */
get_vertex(gint mz,const tetra_vertex_t * v1,const tetra_vertex_t * v2,helper_t * help,GtsCartesianGrid * g,GtsVertexClass * klass)110 static GtsVertex * get_vertex (gint mz,
111 const tetra_vertex_t * v1,
112 const tetra_vertex_t * v2,
113 helper_t * help,
114 GtsCartesianGrid * g,
115 GtsVertexClass * klass)
116 {
117 GtsVertex ** vertex;
118 gint x, y, index, idx2, z;
119 gdouble dx, dy, dz, d;
120
121 g_assert (v1->d - v2->d != 0.);
122
123 dx = dy = dz = 0.0;
124 d = v1->d/(v1->d - v2->d);
125
126 index = 0;
127
128 if (v1->x != v2->x) {
129 index |= 1;
130 dx = d;
131 }
132
133 if (v1->y != v2->y) {
134 index |= 2;
135 dy = d;
136 }
137
138 if (v1->z != v2->z) {
139 dz = d;
140 }
141
142 x = v1->x;
143 if (v1->x > v2->x) { x = v2->x; dx = 1.0 - dx; }
144
145 y = v1->y;
146 if (v1->y > v2->y) { y = v2->y; dy = 1.0 - dy;}
147
148 z = v1->z;
149 if (v1->z > v2->z) { z = v2->z; dz = 1.0 - dz;}
150
151 idx2 = 4 * ( x + y * help->nx ) + index;
152
153 if (v1->z == v2->z)
154 vertex = (mz == z) ? &help->vtop[idx2] : &help->vbot[idx2];
155 else
156 vertex = &help->vmid[idx2];
157
158 if (mz != z && dz != 0.0) {
159 fprintf(stderr, "%f \n", dz);
160 }
161
162 /* if vertex is not yet created, do it now */
163 if (!*vertex)
164 *vertex = gts_vertex_new (klass,
165 g->dx * ( x + dx) + g->x,
166 g->dy * ( y + dy) + g->y,
167 g->dz * ( z + dz) + g->z);
168
169 return *vertex;
170 }
171
get_vertex_bcl(gint mz,const tetra_vertex_t * v1,const tetra_vertex_t * v2,helper_bcl * help,GtsCartesianGrid * g,GtsVertexClass * klass)172 static GtsVertex * get_vertex_bcl (gint mz,
173 const tetra_vertex_t * v1,
174 const tetra_vertex_t * v2,
175 helper_bcl * help,
176 GtsCartesianGrid * g,
177 GtsVertexClass * klass)
178 {
179 GtsVertex * v;
180 GHashTable * table;
181 gchar * s1, * s2, * hash;
182 gdouble x1, x2, y1, y2, z1, z2, d;
183
184 g_assert (v1->d - v2->d != 0.);
185
186 /* first find correct hash table */
187 if ((v1->z > mz) && (v2->z > mz))
188 table = help->vtop;
189 else
190 table = help->vbot;
191
192 d = v1->d / (v1->d - v2->d);
193
194 /* sort vertices */
195 s1 = g_strdup_printf ("%d %d %d %d", v1->x, v1->y, v1->z, v1->mid);
196 s2 = g_strdup_printf ("%d %d %d %d", v2->x, v2->y, v2->z, v2->mid);
197
198 hash = (d == 0.0) ? g_strdup (s1) :
199 (d == 1.0) ? g_strdup (s2) :
200 (strcmp (s1, s2) < 0) ? g_strjoin (" ", s1, s2, NULL) :
201 g_strjoin (" ", s2, s1, NULL);
202
203 /* return existing vertex or make a new one */
204 v = g_hash_table_lookup (table, hash);
205 if (!v){
206
207 x1 = g->dx * (v1->x + (v1->mid / 2.0)) + g->x;
208 x2 = g->dx * (v2->x + (v2->mid / 2.0)) + g->x;
209 y1 = g->dy * (v1->y + (v1->mid / 2.0)) + g->y;
210 y2 = g->dy * (v2->y + (v2->mid / 2.0)) + g->y;
211 z1 = g->dz * (v1->z + (v1->mid / 2.0)) + g->z;
212 z2 = g->dz * (v2->z + (v2->mid / 2.0)) + g->z;
213
214 v = gts_vertex_new (klass, x1 * (1.0 - d) + d * x2,
215 y1 * (1.0 - d) + d * y2,
216 z1 * (1.0 - d) + d * z2);
217
218 g_hash_table_insert (table, g_strdup(hash), v);
219 }
220 g_free (s1);
221 g_free (s2);
222 g_free (hash);
223
224 return v;
225 }
226
227 /* create an edge connecting the zero crossings of lines through a
228 pair of vertices, or return an existing one */
get_edge(GtsVertex * v1,GtsVertex * v2,GtsEdgeClass * klass)229 static GtsEdge * get_edge (GtsVertex * v1, GtsVertex * v2,
230 GtsEdgeClass * klass)
231 {
232 GtsSegment *s;
233 GtsEdge *edge;
234
235 g_assert (v1);
236 g_assert (v2);
237
238 s = gts_vertices_are_connected (v1, v2);
239
240 if (GTS_IS_EDGE (s))
241 edge = GTS_EDGE(s);
242 else
243 edge = gts_edge_new (klass, v1, v2);
244 return edge;
245 }
246
add_face(GtsSurface * surface,const tetra_vertex_t * a1,const tetra_vertex_t * a2,const tetra_vertex_t * b1,const tetra_vertex_t * b2,const tetra_vertex_t * c1,const tetra_vertex_t * c2,gint rev,helper_t * help,gint z,GtsCartesianGrid * g)247 static void add_face (GtsSurface * surface,
248 const tetra_vertex_t * a1, const tetra_vertex_t * a2,
249 const tetra_vertex_t * b1, const tetra_vertex_t * b2,
250 const tetra_vertex_t * c1, const tetra_vertex_t * c2,
251 gint rev, helper_t * help,
252 gint z, GtsCartesianGrid * g)
253 {
254 GtsFace * t;
255 GtsEdge * e1, * e2, * e3;
256 GtsVertex * v1 = get_vertex (z, a1, a2, help, g, surface->vertex_class);
257 GtsVertex * v2 = get_vertex (z, b1, b2, help, g, surface->vertex_class);
258 GtsVertex * v3 = get_vertex (z, c1, c2, help, g, surface->vertex_class);
259
260 g_assert (v1 != v2);
261 g_assert (v2 != v3);
262 g_assert (v1 != v3);
263
264 if (!rev) {
265 e1 = get_edge (v1, v2, surface->edge_class);
266 e2 = get_edge (v2, v3, surface->edge_class);
267 e3 = get_edge (v1, v3, surface->edge_class);
268 } else {
269 e1 = get_edge (v1, v3, surface->edge_class);
270 e2 = get_edge (v2, v3, surface->edge_class);
271 e3 = get_edge (v1, v2, surface->edge_class);
272 }
273
274 t = gts_face_new (surface->face_class, e1, e2, e3);
275 gts_surface_add_face (surface, t);
276 }
277
add_face_bcl(GtsSurface * surface,const tetra_vertex_t * a1,const tetra_vertex_t * a2,const tetra_vertex_t * b1,const tetra_vertex_t * b2,const tetra_vertex_t * c1,const tetra_vertex_t * c2,gint rev,helper_bcl * help,gint z,GtsCartesianGrid * g)278 static void add_face_bcl (GtsSurface * surface,
279 const tetra_vertex_t * a1,
280 const tetra_vertex_t * a2,
281 const tetra_vertex_t * b1,
282 const tetra_vertex_t * b2,
283 const tetra_vertex_t * c1,
284 const tetra_vertex_t * c2,
285 gint rev, helper_bcl * help,
286 gint z, GtsCartesianGrid * g)
287 {
288 GtsFace * t;
289 GtsEdge * e1, * e2, * e3;
290 GtsVertex * v1 = get_vertex_bcl (z, a1, a2, help, g, surface->vertex_class);
291 GtsVertex * v2 = get_vertex_bcl (z, b1, b2, help, g, surface->vertex_class);
292 GtsVertex * v3 = get_vertex_bcl (z, c1, c2, help, g, surface->vertex_class);
293
294 if (v1 == v2 || v2 == v3 || v1 == v3)
295 return;
296
297 if (!rev) {
298 e1 = get_edge (v1, v2, surface->edge_class);
299 e2 = get_edge (v2, v3, surface->edge_class);
300 e3 = get_edge (v1, v3, surface->edge_class);
301 } else {
302 e1 = get_edge (v1, v3, surface->edge_class);
303 e2 = get_edge (v2, v3, surface->edge_class);
304 e3 = get_edge (v1, v2, surface->edge_class);
305 }
306
307 t = gts_face_new (surface->face_class, e1, e2, e3);
308 gts_surface_add_face (surface, t);
309 }
310
311 /* create a new slice of site nx \times ny */
new_slice(gint nx,gint ny)312 static slice_t * new_slice (gint nx, gint ny)
313 {
314 gint x;
315 slice_t * retval = g_malloc (sizeof (slice_t));
316
317 retval->data = g_malloc (nx*sizeof(gdouble *));
318 retval->nx = nx;
319 retval->ny = ny;
320 for (x = 0; x < nx; x++)
321 retval->data[x] = g_malloc (ny*sizeof (gdouble));
322 return retval;
323 }
324
325 /* initialize a slice with inival */
slice_init(slice_t * slice,gdouble inival)326 static void slice_init (slice_t * slice, gdouble inival)
327 {
328 gint x, y;
329
330 g_assert (slice);
331
332 for (x = 0; x < slice->nx; x++)
333 for (y = 0; y < slice->ny; y++)
334 slice->data[x][y] = inival;
335 }
336
337 /* free the memory of a slice */
free_slice(slice_t * slice)338 static void free_slice (slice_t * slice)
339 {
340 gint x;
341
342 g_return_if_fail (slice != NULL);
343
344 for (x = 0; x < slice->nx; x++)
345 g_free (slice->data[x]);
346 g_free (slice->data);
347 g_free (slice);
348 }
349
analyze_tetrahedra(const tetra_vertex_t * a,const tetra_vertex_t * b,const tetra_vertex_t * c,const tetra_vertex_t * d,gint parity,GtsSurface * surface,helper_t * help,gint z,GtsCartesianGrid * g)350 static void analyze_tetrahedra (const tetra_vertex_t * a,
351 const tetra_vertex_t * b,
352 const tetra_vertex_t * c,
353 const tetra_vertex_t * d,
354 gint parity, GtsSurface * surface,
355 helper_t * help,
356 gint z, GtsCartesianGrid * g)
357 {
358 gint rev = parity;
359 gint code = 0;
360
361 if (a->d >= 0.) code |= 1;
362 if (b->d >= 0.) code |= 2;
363 if (c->d >= 0.) code |= 4;
364 if (d->d >= 0.) code |= 8;
365
366 switch (code) {
367 case 15:
368 case 0: return; /* all inside or outside */
369
370 case 14:rev = !parity;
371 case 1:add_face (surface, a, b, a, d, a, c, rev, help, z, g);
372 break;
373 case 13:rev = ! parity;
374 case 2:add_face (surface, a, b, b, c, b, d, rev, help, z, g);
375 break;
376 case 12:rev = !parity;
377 case 3:add_face (surface, a, d, a, c, b, c, rev, help, z, g);
378 add_face (surface, a, d, b, c, b, d, rev, help, z, g);
379 break;
380 case 11:rev = !parity;
381 case 4:add_face (surface, a, c, c, d, b, c, rev, help, z, g);
382 break;
383 case 10:rev = !parity;
384 case 5: add_face (surface, a, b, a, d, c, d, rev, help, z, g);
385 add_face (surface, a, b, c, d, b, c, rev, help, z, g);
386 break;
387 case 9:rev = !parity;
388 case 6:add_face (surface, a, b, a, c, c, d, rev, help, z, g);
389 add_face (surface, a, b, c, d, b, d, rev, help, z, g);
390 break;
391 case 7:rev = !parity;
392 case 8:add_face (surface, a, d, b, d, c, d, rev, help, z, g);
393 break;
394 }
395 }
396
analyze_tetrahedra_bcl(const tetra_vertex_t * a,const tetra_vertex_t * b,const tetra_vertex_t * c,const tetra_vertex_t * d,GtsSurface * surface,helper_bcl * help,gint z,GtsCartesianGrid * g)397 static void analyze_tetrahedra_bcl (const tetra_vertex_t * a,
398 const tetra_vertex_t * b,
399 const tetra_vertex_t * c,
400 const tetra_vertex_t * d,
401 GtsSurface * surface,
402 helper_bcl * help,
403 gint z, GtsCartesianGrid * g)
404 {
405 gint rev = 0;
406 gint code = 0;
407
408 if (a->d >= 0.) code |= 1;
409 if (b->d >= 0.) code |= 2;
410 if (c->d >= 0.) code |= 4;
411 if (d->d >= 0.) code |= 8;
412
413 switch (code) {
414 case 15:
415 case 0: return; /* all inside or outside */
416
417 case 14:rev = !rev;
418 case 1:add_face_bcl (surface, a, b, a, d, a, c, rev, help, z, g);
419 break;
420 case 13:rev = !rev;
421 case 2:add_face_bcl (surface, a, b, b, c, b, d, rev, help, z, g);
422 break;
423 case 12:rev = !rev;
424 case 3:add_face_bcl (surface, a, d, a, c, b, c, rev, help, z, g);
425 add_face_bcl (surface, a, d, b, c, b, d, rev, help, z, g);
426 break;
427 case 11:rev = !rev;
428 case 4:add_face_bcl (surface, a, c, c, d, b, c, rev, help, z, g);
429 break;
430 case 10:rev = !rev;
431 case 5: add_face_bcl (surface, a, b, a, d, c, d, rev, help, z, g);
432 add_face_bcl (surface, a, b, c, d, b, c, rev, help, z, g);
433 break;
434 case 9:rev = !rev;
435 case 6:add_face_bcl (surface, a, b, a, c, c, d, rev, help, z, g);
436 add_face_bcl (surface, a, b, c, d, b, d, rev, help, z, g);
437 break;
438 case 7:rev = !rev;
439 case 8:add_face_bcl (surface, a, d, b, d, c, d, rev, help, z, g);
440 break;
441 }
442 }
443
iso_slice_evaluate(slice_t * s1,slice_t * s2,GtsCartesianGrid g,gint z,GtsSurface * surface,helper_t * help)444 static void iso_slice_evaluate (slice_t * s1, slice_t * s2,
445 GtsCartesianGrid g,
446 gint z, GtsSurface * surface, helper_t * help)
447 {
448 gint x,y;
449 tetra_vertex_t v0, v1, v2, v3, v4, v5, v6, v7;
450 gdouble ** s1p = s1->data;
451 gdouble ** s2p = s2->data;
452
453 for (y = 0; y < g.ny-1; y++)
454 for (x = 0; x < g.nx-1; x++) {
455 gint parity = (((x ^ y) ^ z) & 1);
456
457 v0.x = x ; v0.y = y ; v0.z = z ; v0.mid = FALSE; v0.d = s1p[x ][y ];
458 v1.x = x ; v1.y = y+1; v1.z = z ; v1.mid = FALSE; v1.d = s1p[x ][y+1];
459 v2.x = x+1; v2.y = y ; v2.z = z ; v2.mid = FALSE; v2.d = s1p[x+1][y ];
460 v3.x = x+1; v3.y = y+1; v3.z = z ; v3.mid = FALSE; v3.d = s1p[x+1][y+1];
461 v4.x = x ; v4.y = y ; v4.z = z+1; v4.mid = FALSE; v4.d = s2p[x ][y ];
462 v5.x = x ; v5.y = y+1; v5.z = z+1; v5.mid = FALSE; v5.d = s2p[x ][y+1];
463 v6.x = x+1; v6.y = y ; v6.z = z+1; v6.mid = FALSE; v6.d = s2p[x+1][y ];
464 v7.x = x+1; v7.y = y+1; v7.z = z+1; v7.mid = FALSE; v7.d = s2p[x+1][y+1];
465
466 if (parity == 0) {
467 analyze_tetrahedra (&v0, &v1, &v2, &v4, parity, surface, help, z, &g);
468 analyze_tetrahedra (&v7, &v1, &v4, &v2, parity, surface, help, z, &g);
469 analyze_tetrahedra (&v1, &v7, &v3, &v2, parity, surface, help, z, &g);
470 analyze_tetrahedra (&v1, &v7, &v4, &v5, parity, surface, help, z, &g);
471 analyze_tetrahedra (&v2, &v6, &v4, &v7, parity, surface, help, z, &g);
472 }else{
473 analyze_tetrahedra (&v4, &v5, &v6, &v0, parity, surface, help, z, &g);
474 analyze_tetrahedra (&v3, &v5, &v0, &v6, parity, surface, help, z, &g);
475 analyze_tetrahedra (&v5, &v3, &v7, &v6, parity, surface, help, z, &g);
476 analyze_tetrahedra (&v5, &v3, &v0, &v1, parity, surface, help, z, &g);
477 analyze_tetrahedra (&v6, &v2, &v0, &v3, parity, surface, help, z, &g);
478 }
479 }
480 }
481
iso_slice_evaluate_bcl(slice_t * s1,slice_t * s2,slice_t * s3,GtsCartesianGrid g,gint z,GtsSurface * surface,helper_bcl * help)482 static void iso_slice_evaluate_bcl (slice_t * s1, slice_t * s2, slice_t * s3,
483 GtsCartesianGrid g,
484 gint z, GtsSurface * surface,
485 helper_bcl * help)
486 {
487 gint x,y;
488 tetra_vertex_t v0, v1, v2, v3, v4, v5, v6, v7, v8, v9, w0;
489 gdouble ** s1p = s1->data;
490 gdouble ** s2p = s2->data;
491 gdouble ** s3p = s3->data;
492
493 for (y = 0; y < g.ny-2; y++)
494 for (x = 0; x < g.nx-2; x++) {
495 v0.x = x ; v0.y = y ; v0.z = z ; v0.mid = TRUE;
496 v0.d = (s1p[x ][y ] + s2p[x ][y ] +
497 s1p[x ][y+1] + s2p[x ][y+1] +
498 s1p[x+1][y ] + s2p[x+1][y ] +
499 s1p[x+1][y+1] + s2p[x+1][y+1])/8.0;
500
501 v1.x = x+1; v1.y = y ; v1.z = z ; v1.mid = TRUE;
502 v1.d = (s1p[x+1][y ] + s2p[x+1][y ] +
503 s1p[x+1][y+1] + s2p[x+1][y+1] +
504 s1p[x+2][y ] + s2p[x+2][y ] +
505 s1p[x+2][y+1] + s2p[x+2][y+1])/8.0;
506
507 v2.x = x ; v2.y = y+1; v2.z = z ; v2.mid = TRUE;
508 v2.d = (s1p[x ][y+1] + s2p[x ][y+1] +
509 s1p[x ][y+2] + s2p[x ][y+2] +
510 s1p[x+1][y+1] + s2p[x+1][y+1] +
511 s1p[x+1][y+2] + s2p[x+1][y+2])/8.0;
512
513 v3.x = x ; v3.y = y ; v3.z = z+1; v3.mid = TRUE;
514 v3.d = (s2p[x ][y ] + s3p[x ][y ] +
515 s2p[x ][y+1] + s3p[x ][y+1] +
516 s2p[x+1][y ] + s3p[x+1][y ] +
517 s2p[x+1][y+1] + s3p[x+1][y+1])/8.0;
518
519 v4.x = x+1; v4.y = y ; v4.z = z ; v4.mid = FALSE; v4.d = s1p[x+1][y ];
520 v5.x = x ; v5.y = y+1; v5.z = z ; v5.mid = FALSE; v5.d = s1p[x ][y+1];
521 v6.x = x+1; v6.y = y+1; v6.z = z ; v6.mid = FALSE; v6.d = s1p[x+1][y+1];
522 v7.x = x+1; v7.y = y ; v7.z = z+1; v7.mid = FALSE; v7.d = s2p[x+1][y ];
523 v8.x = x ; v8.y = y+1; v8.z = z+1; v8.mid = FALSE; v8.d = s2p[x ][y+1];
524 v9.x = x+1; v9.y = y+1; v9.z = z+1; v9.mid = FALSE; v9.d = s2p[x+1][y+1];
525 w0.x = x ; w0.y = y ; w0.z = z+1; w0.mid = FALSE; w0.d = s2p[x ][y ];
526
527 analyze_tetrahedra_bcl (&v0, &v9, &v6, &v1, surface, help, z, &g);
528 analyze_tetrahedra_bcl (&v0, &v6, &v4, &v1, surface, help, z, &g);
529 analyze_tetrahedra_bcl (&v0, &v4, &v7, &v1, surface, help, z, &g);
530 analyze_tetrahedra_bcl (&v0, &v7, &v9, &v1, surface, help, z, &g);
531 analyze_tetrahedra_bcl (&v0, &v5, &v6, &v2, surface, help, z, &g);
532 analyze_tetrahedra_bcl (&v0, &v6, &v9, &v2, surface, help, z, &g);
533 analyze_tetrahedra_bcl (&v0, &v9, &v8, &v2, surface, help, z, &g);
534 analyze_tetrahedra_bcl (&v0, &v8, &v5, &v2, surface, help, z, &g);
535 analyze_tetrahedra_bcl (&v0, &v8, &v9, &v3, surface, help, z, &g);
536 analyze_tetrahedra_bcl (&v0, &v9, &v7, &v3, surface, help, z, &g);
537 analyze_tetrahedra_bcl (&v0, &v7, &w0, &v3, surface, help, z, &g);
538 analyze_tetrahedra_bcl (&v0, &w0, &v8, &v3, surface, help, z, &g);
539 }
540 }
541
542 /* copy src into dest by stripping off the iso value and leave out
543 the boundary (which should be G_MINDOUBLE) */
copy_to_bounded(slice_t * dest,slice_t * src,gdouble iso,gdouble fill)544 static void copy_to_bounded (slice_t * dest, slice_t * src,
545 gdouble iso, gdouble fill)
546 {
547 gint x,y;
548 gdouble * src_ptr;
549 gdouble * dest_ptr = dest->data[0];
550
551 g_assert(dest->ny == src->ny + 2);
552 g_assert(dest->nx == src->nx + 2);
553
554 for (y = 0; y < dest->ny; ++y, ++dest_ptr)
555 *dest_ptr = fill;
556
557 for (x = 1; x < src->nx - 1; ++x) {
558 dest_ptr = dest->data[x];
559 src_ptr = src->data[x-1];
560 *dest_ptr++ = fill;
561 for (y = 0; y < src->ny; ++y, ++dest_ptr, ++src_ptr)
562 *dest_ptr = *src_ptr - iso;
563 *dest_ptr++ = fill;
564 }
565
566 dest_ptr = dest->data[y];
567
568 for (y = 0; y < dest->ny; ++y, ++dest_ptr)
569 *dest_ptr = fill;
570 }
571
iso_sub(slice_t * s,gdouble iso)572 static void iso_sub (slice_t * s, gdouble iso)
573 {
574 gint x,y;
575
576 for (x = 0; x < s->nx; ++x) {
577 gdouble *ptr = s->data[x];
578
579 for (y = 0; y < s->ny; ++y, ++ptr)
580 *ptr -= iso;
581 }
582 }
583
584
585 /**
586 * gts_isosurface_tetra_bounded:
587 * @surface: a #GtsSurface.
588 * @g: a #GtsCartesianGrid.
589 * @f: a #GtsIsoCartesianFunc.
590 * @data: user data to be passed to @f.
591 * @iso: isosurface value.
592 *
593 * Adds to @surface new faces defining the isosurface f(x,y,z) =
594 * @iso. By convention, the normals to the surface are pointing toward
595 * the positive values of f(x,y,z) - @iso. To ensure a closed object,
596 * a boundary of G_MINDOUBLE is added around the domain
597 *
598 * The user function @f is called successively for each value of the z
599 * coordinate defined by @g. It must fill the corresponding (x,y)
600 * plane with the values of the function for which the isosurface is
601 * to be computed.
602 */
gts_isosurface_tetra_bounded(GtsSurface * surface,GtsCartesianGrid g,GtsIsoCartesianFunc f,gpointer data,gdouble iso)603 void gts_isosurface_tetra_bounded (GtsSurface * surface,
604 GtsCartesianGrid g,
605 GtsIsoCartesianFunc f,
606 gpointer data,
607 gdouble iso)
608 {
609 slice_t *slice1, *slice2, *transfer_slice;
610 GtsCartesianGrid g_intern = g;
611 helper_t *helper;
612 gint z;
613
614 g_return_if_fail (surface != NULL);
615 g_return_if_fail (f != NULL);
616 g_return_if_fail (g.nx > 1);
617 g_return_if_fail (g.ny > 1);
618 g_return_if_fail (g.nz > 1);
619
620 /* create the helper slices */
621 slice1 = new_slice (g.nx + 2, g.ny + 2);
622 slice2 = new_slice (g.nx + 2, g.ny + 2);
623
624 /* initialize the first slice as OUTSIDE */
625 slice_init (slice1, -1.0);
626
627 /* create a slice of the original image size */
628 transfer_slice = new_slice (g.nx, g.ny);
629
630 /* adapt the parameters to our enlarged image */
631 g_intern.x -= g.dx;
632 g_intern.y -= g.dy;
633 g_intern.z -= g.dz;
634 g_intern.nx = g.nx + 2;
635 g_intern.ny = g.ny + 2;
636 g_intern.nz = g.nz;
637
638 /* create the helper for vertex-lookup */
639 helper = init_helper (g_intern.nx, g_intern.ny);
640
641 /* go slicewise through the data */
642 z = 0;
643 while (z < g.nz) {
644 slice_t * hs;
645
646 /* request slice */
647 f (transfer_slice->data, g, z, data);
648 g.z += g.dz;
649
650 /* copy slice in enlarged image and mark the border as OUTSIDE */
651 copy_to_bounded (slice2, transfer_slice, iso, -1.);
652
653 /* triangulate */
654 iso_slice_evaluate (slice1, slice2, g_intern, z, surface, helper);
655
656 /* switch the input slices */
657 hs = slice1; slice1 = slice2; slice2 = hs;
658
659 /* switch the vertex lookup tables */
660 helper_advance(helper);
661 ++z;
662 }
663
664 /* initialize the last slice as OUTSIDE */
665 slice_init (slice2, - 1.0);
666
667 /* close the object */
668 iso_slice_evaluate(slice1, slice2, g_intern, z, surface, helper);
669
670 free_helper (helper);
671 free_slice (slice1);
672 free_slice (slice2);
673 free_slice (transfer_slice);
674 }
675
676 /**
677 * gts_isosurface_tetra:
678 * @surface: a #GtsSurface.
679 * @g: a #GtsCartesianGrid.
680 * @f: a #GtsIsoCartesianFunc.
681 * @data: user data to be passed to @f.
682 * @iso: isosurface value.
683 *
684 * Adds to @surface new faces defining the isosurface f(x,y,z) =
685 * @iso. By convention, the normals to the surface are pointing toward
686 * the positive values of f(x,y,z) - @iso.
687 *
688 * The user function @f is called successively for each value of the z
689 * coordinate defined by @g. It must fill the corresponding (x,y)
690 * plane with the values of the function for which the isosurface is
691 * to be computed.
692 */
gts_isosurface_tetra(GtsSurface * surface,GtsCartesianGrid g,GtsIsoCartesianFunc f,gpointer data,gdouble iso)693 void gts_isosurface_tetra (GtsSurface * surface,
694 GtsCartesianGrid g,
695 GtsIsoCartesianFunc f,
696 gpointer data,
697 gdouble iso)
698 {
699 slice_t *slice1, *slice2;
700 helper_t *helper;
701 gint z;
702 GtsCartesianGrid g_internal;
703
704 g_return_if_fail (surface != NULL);
705 g_return_if_fail (f != NULL);
706 g_return_if_fail (g.nx > 1);
707 g_return_if_fail (g.ny > 1);
708 g_return_if_fail (g.nz > 1);
709
710 memcpy (&g_internal, &g, sizeof (GtsCartesianGrid));
711
712 /* create the helper slices */
713 slice1 = new_slice (g.nx, g.ny);
714 slice2 = new_slice (g.nx, g.ny);
715
716 /* create the helper for vertex-lookup */
717 helper = init_helper (g.nx, g.ny);
718
719 z = 0;
720 f (slice1->data, g, z, data);
721 iso_sub (slice1, iso);
722
723 z++;
724 g.z += g.dz;
725
726 /* go slicewise through the data */
727 while (z < g.nz) {
728 slice_t * hs;
729
730 /* request slice */
731 f (slice2->data, g, z, data);
732 iso_sub (slice2, iso);
733
734 g.z += g.dz;
735
736 /* triangulate */
737 iso_slice_evaluate (slice1, slice2, g_internal, z-1, surface, helper);
738
739 /* switch the input slices */
740 hs = slice1; slice1 = slice2; slice2 = hs;
741
742 /* switch the vertex lookup tables */
743 helper_advance (helper);
744
745 ++z;
746 }
747
748 free_helper(helper);
749 free_slice(slice1);
750 free_slice(slice2);
751 }
752
753 /**
754 * gts_isosurface_tetra_bcl:
755 * @surface: a #GtsSurface.
756 * @g: a #GtsCartesianGrid.
757 * @f: a #GtsIsoCartesianFunc.
758 * @data: user data to be passed to @f.
759 * @iso: isosurface value.
760 *
761 * Adds to @surface new faces defining the isosurface f(x,y,z) =
762 * @iso. By convention, the normals to the surface are pointing toward
763 * the positive values of f(x,y,z) - @iso.
764 *
765 * The user function @f is called successively for each value of the z
766 * coordinate defined by @g. It must fill the corresponding (x,y)
767 * plane with the values of the function for which the isosurface is
768 * to be computed.
769 *
770 * This version produces the dual "body-centered" faces relative to
771 * the faces produced by gts_isosurface_tetra().
772 */
gts_isosurface_tetra_bcl(GtsSurface * surface,GtsCartesianGrid g,GtsIsoCartesianFunc f,gpointer data,gdouble iso)773 void gts_isosurface_tetra_bcl (GtsSurface * surface,
774 GtsCartesianGrid g,
775 GtsIsoCartesianFunc f,
776 gpointer data,
777 gdouble iso)
778 {
779 slice_t *slice1, *slice2, *slice3;
780 helper_bcl *helper;
781 gint z;
782 GtsCartesianGrid g_internal;
783
784 g_return_if_fail (surface != NULL);
785 g_return_if_fail (f != NULL);
786 g_return_if_fail (g.nx > 1);
787 g_return_if_fail (g.ny > 1);
788 g_return_if_fail (g.nz > 1);
789
790 memcpy (&g_internal, &g, sizeof (GtsCartesianGrid));
791
792 /* create the helper slices */
793 slice1 = new_slice (g.nx, g.ny);
794 slice2 = new_slice (g.nx, g.ny);
795 slice3 = new_slice (g.nx, g.ny);
796
797 /* create the helper for vertex-lookup */
798 helper = init_helper_bcl ();
799
800 z = 0;
801 f (slice1->data, g, z, data);
802 iso_sub (slice1, iso);
803
804 z++;
805 g.z += g.dz;
806
807 f (slice2->data, g, z, data);
808 iso_sub (slice1, iso);
809
810 z++;
811 g.z += g.dz;
812
813 /* go slicewise through the data */
814 while (z < g.nz) {
815 slice_t * hs;
816
817 /* request slice */
818 f (slice3->data, g, z, data);
819 iso_sub (slice3, iso);
820
821 g.z += g.dz;
822
823 /* triangulate */
824 iso_slice_evaluate_bcl (slice1, slice2, slice3, g_internal, z-2,
825 surface, helper);
826
827 /* switch the input slices */
828 hs = slice1; slice1 = slice2; slice2 = slice3; slice3 = hs;
829
830 /* switch the vertex lookup tables */
831 helper_advance_bcl (helper);
832
833 ++z;
834 }
835
836 free_helper_bcl(helper);
837 free_slice(slice1);
838 free_slice(slice2);
839 free_slice(slice3);
840 }
841