gts/isotetra.c

   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
  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
  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
  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
  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 */
  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
  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 */
 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
 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 */
 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
 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
 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 */
 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 */
 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 */
 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
 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
 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
 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
 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) */
 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
 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  */
 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  */
 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  */
 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 }