1 /* GTS-Library conform marching tetrahedra algorithm
2 * Copyright (C) 2002 Gert Wollny
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.
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.
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.
25 # define M_SQRT2 1.41421356237309504880
26 #endif /* NATIVE_WIN32 */
39 /* this helper is a lookup table for vertices */
42 GtsVertex
** vtop
, ** vmid
, **vbot
;
46 GHashTable
* vbot
, * vtop
;
50 static helper_t
* init_helper (int nx
, int ny
)
53 helper_t
*retval
= g_malloc0 (sizeof (helper_t
));
57 retval
->vtop
= g_malloc0 (sizeof (GtsVertex
*)*nxy
);
58 retval
->vmid
= g_malloc0 (sizeof (GtsVertex
*)*nxy
);
59 retval
->vbot
= g_malloc0 (sizeof (GtsVertex
*)*nxy
);
63 static helper_bcl
* init_helper_bcl (void)
65 helper_bcl
*retval
= g_malloc0 (sizeof (helper_bcl
));
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
);
72 static void free_helper (helper_t
* h
)
80 static void free_helper_bcl (helper_bcl
* h
)
82 g_hash_table_destroy (h
->vtop
);
83 g_hash_table_destroy (h
->vbot
);
87 /* move the vertices in the bottom slice to the top, and clear the
88 other slices in the lookup tables */
89 static void helper_advance (helper_t
* h
)
91 GtsVertex
** help
= h
->vbot
;
95 memset (h
->vmid
, 0, 4*sizeof(GtsVertex
*) * h
->nx
* h
->ny
);
96 memset (h
->vbot
, 0, 4*sizeof(GtsVertex
*) * h
->nx
* h
->ny
);
99 static void helper_advance_bcl (helper_bcl
* h
)
101 GHashTable
* help
= g_hash_table_new (g_str_hash
, g_str_equal
);
103 g_hash_table_destroy (h
->vbot
);
108 /* find the zero-crossing of line through v1 and v2 and return the
109 corresponding GtsVertex */
110 static GtsVertex
* get_vertex (gint mz
,
111 const tetra_vertex_t
* v1
,
112 const tetra_vertex_t
* v2
,
114 GtsCartesianGrid
* g
,
115 GtsVertexClass
* klass
)
118 gint x
, y
, index
, idx2
, z
;
119 gdouble dx
, dy
, dz
, d
;
121 g_assert (v1
->d
- v2
->d
!= 0.);
124 d
= v1
->d
/(v1
->d
- v2
->d
);
128 if (v1
->x
!= v2
->x
) {
133 if (v1
->y
!= v2
->y
) {
138 if (v1
->z
!= v2
->z
) {
143 if (v1
->x
> v2
->x
) { x
= v2
->x
; dx
= 1.0 - dx
; }
146 if (v1
->y
> v2
->y
) { y
= v2
->y
; dy
= 1.0 - dy
;}
149 if (v1
->z
> v2
->z
) { z
= v2
->z
; dz
= 1.0 - dz
;}
151 idx2
= 4 * ( x
+ y
* help
->nx
) + index
;
154 vertex
= (mz
== z
) ? &help
->vtop
[idx2
] : &help
->vbot
[idx2
];
156 vertex
= &help
->vmid
[idx2
];
158 if (mz
!= z
&& dz
!= 0.0) {
159 fprintf(stderr
, "%f \n", dz
);
162 /* if vertex is not yet created, do it now */
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
);
172 static GtsVertex
* get_vertex_bcl (gint mz
,
173 const tetra_vertex_t
* v1
,
174 const tetra_vertex_t
* v2
,
176 GtsCartesianGrid
* g
,
177 GtsVertexClass
* klass
)
181 gchar
* s1
, * s2
, * hash
;
182 gdouble x1
, x2
, y1
, y2
, z1
, z2
, d
;
184 g_assert (v1
->d
- v2
->d
!= 0.);
186 /* first find correct hash table */
187 if ((v1
->z
> mz
) && (v2
->z
> mz
))
192 d
= v1
->d
/ (v1
->d
- v2
->d
);
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
);
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
);
203 /* return existing vertex or make a new one */
204 v
= g_hash_table_lookup (table
, hash
);
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
;
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
);
218 g_hash_table_insert (table
, g_strdup(hash
), v
);
227 /* create an edge connecting the zero crossings of lines through a
228 pair of vertices, or return an existing one */
229 static GtsEdge
* get_edge (GtsVertex
* v1
, GtsVertex
* v2
,
230 GtsEdgeClass
* klass
)
238 s
= gts_vertices_are_connected (v1
, v2
);
243 edge
= gts_edge_new (klass
, v1
, v2
);
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
)
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
);
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
);
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
);
274 t
= gts_face_new (surface
->face_class
, e1
, e2
, e3
);
275 gts_surface_add_face (surface
, t
);
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
)
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
);
294 if (v1
== v2
|| v2
== v3
|| v1
== v3
)
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
);
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
);
307 t
= gts_face_new (surface
->face_class
, e1
, e2
, e3
);
308 gts_surface_add_face (surface
, t
);
311 /* create a new slice of site nx \times ny */
312 static slice_t
* new_slice (gint nx
, gint ny
)
315 slice_t
* retval
= g_malloc (sizeof (slice_t
));
317 retval
->data
= g_malloc (nx
*sizeof(gdouble
*));
320 for (x
= 0; x
< nx
; x
++)
321 retval
->data
[x
] = g_malloc (ny
*sizeof (gdouble
));
325 /* initialize a slice with inival */
326 static void slice_init (slice_t
* slice
, gdouble inival
)
332 for (x
= 0; x
< slice
->nx
; x
++)
333 for (y
= 0; y
< slice
->ny
; y
++)
334 slice
->data
[x
][y
] = inival
;
337 /* free the memory of a slice */
338 static void free_slice (slice_t
* slice
)
342 g_return_if_fail (slice
!= NULL
);
344 for (x
= 0; x
< slice
->nx
; x
++)
345 g_free (slice
->data
[x
]);
346 g_free (slice
->data
);
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
,
356 gint z
, GtsCartesianGrid
* g
)
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;
368 case 0: return; /* all inside or outside */
370 case 14:rev
= !parity
;
371 case 1:add_face (surface
, a
, b
, a
, d
, a
, c
, rev
, help
, z
, g
);
373 case 13:rev
= ! parity
;
374 case 2:add_face (surface
, a
, b
, b
, c
, b
, d
, rev
, help
, z
, g
);
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
);
380 case 11:rev
= !parity
;
381 case 4:add_face (surface
, a
, c
, c
, d
, b
, c
, rev
, help
, z
, g
);
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
);
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
);
391 case 7:rev
= !parity
;
392 case 8:add_face (surface
, a
, d
, b
, d
, c
, d
, rev
, help
, z
, 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
,
403 gint z
, GtsCartesianGrid
* g
)
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;
415 case 0: return; /* all inside or outside */
418 case 1:add_face_bcl (surface
, a
, b
, a
, d
, a
, c
, rev
, help
, z
, g
);
421 case 2:add_face_bcl (surface
, a
, b
, b
, c
, b
, d
, rev
, help
, z
, g
);
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
);
428 case 4:add_face_bcl (surface
, a
, c
, c
, d
, b
, c
, rev
, help
, z
, g
);
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
);
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
);
439 case 8:add_face_bcl (surface
, a
, d
, b
, d
, c
, d
, rev
, help
, z
, g
);
444 static void iso_slice_evaluate (slice_t
* s1
, slice_t
* s2
,
446 gint z
, GtsSurface
* surface
, helper_t
* help
)
449 tetra_vertex_t v0
, v1
, v2
, v3
, v4
, v5
, v6
, v7
;
450 gdouble
** s1p
= s1
->data
;
451 gdouble
** s2p
= s2
->data
;
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);
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];
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
);
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
);
482 static void iso_slice_evaluate_bcl (slice_t
* s1
, slice_t
* s2
, slice_t
* s3
,
484 gint z
, GtsSurface
* surface
,
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
;
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;
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;
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;
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;
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
];
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
);
542 /* copy src into dest by stripping off the iso value and leave out
543 the boundary (which should be G_MINDOUBLE) */
544 static void copy_to_bounded (slice_t
* dest
, slice_t
* src
,
545 gdouble iso
, gdouble fill
)
549 gdouble
* dest_ptr
= dest
->data
[0];
551 g_assert(dest
->ny
== src
->ny
+ 2);
552 g_assert(dest
->nx
== src
->nx
+ 2);
554 for (y
= 0; y
< dest
->ny
; ++y
, ++dest_ptr
)
557 for (x
= 1; x
< src
->nx
- 1; ++x
) {
558 dest_ptr
= dest
->data
[x
];
559 src_ptr
= src
->data
[x
-1];
561 for (y
= 0; y
< src
->ny
; ++y
, ++dest_ptr
, ++src_ptr
)
562 *dest_ptr
= *src_ptr
- iso
;
566 dest_ptr
= dest
->data
[y
];
568 for (y
= 0; y
< dest
->ny
; ++y
, ++dest_ptr
)
572 static void iso_sub (slice_t
* s
, gdouble iso
)
576 for (x
= 0; x
< s
->nx
; ++x
) {
577 gdouble
*ptr
= s
->data
[x
];
579 for (y
= 0; y
< s
->ny
; ++y
, ++ptr
)
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.
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
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
603 void gts_isosurface_tetra_bounded (GtsSurface
* surface
,
605 GtsIsoCartesianFunc f
,
609 slice_t
*slice1
, *slice2
, *transfer_slice
;
610 GtsCartesianGrid g_intern
= g
;
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);
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);
624 /* initialize the first slice as OUTSIDE */
625 slice_init (slice1
, -1.0);
627 /* create a slice of the original image size */
628 transfer_slice
= new_slice (g
.nx
, g
.ny
);
630 /* adapt the parameters to our enlarged image */
634 g_intern
.nx
= g
.nx
+ 2;
635 g_intern
.ny
= g
.ny
+ 2;
638 /* create the helper for vertex-lookup */
639 helper
= init_helper (g_intern
.nx
, g_intern
.ny
);
641 /* go slicewise through the data */
647 f (transfer_slice
->data
, g
, z
, data
);
650 /* copy slice in enlarged image and mark the border as OUTSIDE */
651 copy_to_bounded (slice2
, transfer_slice
, iso
, -1.);
654 iso_slice_evaluate (slice1
, slice2
, g_intern
, z
, surface
, helper
);
656 /* switch the input slices */
657 hs
= slice1
; slice1
= slice2
; slice2
= hs
;
659 /* switch the vertex lookup tables */
660 helper_advance(helper
);
664 /* initialize the last slice as OUTSIDE */
665 slice_init (slice2
, - 1.0);
667 /* close the object */
668 iso_slice_evaluate(slice1
, slice2
, g_intern
, z
, surface
, helper
);
670 free_helper (helper
);
673 free_slice (transfer_slice
);
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.
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.
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
693 void gts_isosurface_tetra (GtsSurface
* surface
,
695 GtsIsoCartesianFunc f
,
699 slice_t
*slice1
, *slice2
;
702 GtsCartesianGrid g_internal
;
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);
710 memcpy (&g_internal
, &g
, sizeof (GtsCartesianGrid
));
712 /* create the helper slices */
713 slice1
= new_slice (g
.nx
, g
.ny
);
714 slice2
= new_slice (g
.nx
, g
.ny
);
716 /* create the helper for vertex-lookup */
717 helper
= init_helper (g
.nx
, g
.ny
);
720 f (slice1
->data
, g
, z
, data
);
721 iso_sub (slice1
, iso
);
726 /* go slicewise through the data */
731 f (slice2
->data
, g
, z
, data
);
732 iso_sub (slice2
, iso
);
737 iso_slice_evaluate (slice1
, slice2
, g_internal
, z
-1, surface
, helper
);
739 /* switch the input slices */
740 hs
= slice1
; slice1
= slice2
; slice2
= hs
;
742 /* switch the vertex lookup tables */
743 helper_advance (helper
);
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.
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.
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
770 * This version produces the dual "body-centered" faces relative to
771 * the faces produced by gts_isosurface_tetra().
773 void gts_isosurface_tetra_bcl (GtsSurface
* surface
,
775 GtsIsoCartesianFunc f
,
779 slice_t
*slice1
, *slice2
, *slice3
;
782 GtsCartesianGrid g_internal
;
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);
790 memcpy (&g_internal
, &g
, sizeof (GtsCartesianGrid
));
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
);
797 /* create the helper for vertex-lookup */
798 helper
= init_helper_bcl ();
801 f (slice1
->data
, g
, z
, data
);
802 iso_sub (slice1
, iso
);
807 f (slice2
->data
, g
, z
, data
);
808 iso_sub (slice1
, iso
);
813 /* go slicewise through the data */
818 f (slice3
->data
, g
, z
, data
);
819 iso_sub (slice3
, iso
);
824 iso_slice_evaluate_bcl (slice1
, slice2
, slice3
, g_internal
, z
-2,
827 /* switch the input slices */
828 hs
= slice1
; slice1
= slice2
; slice2
= slice3
; slice3
= hs
;
830 /* switch the vertex lookup tables */
831 helper_advance_bcl (helper
);
836 free_helper_bcl(helper
);