pcb 4.1.1
An interactive printed circuit board layout editor.

isotetra.c

Go to the documentation of this file.
00001 /* GTS-Library conform marching tetrahedra algorithm 
00002  * Copyright (C) 2002 Gert Wollny
00003  *
00004  * This library is free software; you can redistribute it and/or
00005  * modify it under the terms of the GNU Library General Public
00006  * License as published by the Free Software Foundation; either
00007  * version 2 of the License, or (at your option) any later version.
00008  *
00009  * This library is distributed in the hope that it will be useful,
00010  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00011  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00012  * Library General Public License for more details.
00013  *
00014  * You should have received a copy of the GNU Library General Public
00015  * License along with this library; if not, write to the
00016  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
00017  * Boston, MA 02111-1307, USA.
00018  */
00019 
00020 #include <math.h>
00021 #include <string.h>
00022 #include <gts.h>
00023 #ifdef NATIVE_WIN32
00024 # include <memory.h>
00025 # define M_SQRT2                1.41421356237309504880
00026 #endif /* NATIVE_WIN32 */
00027 
00028 typedef struct {
00029   gint nx, ny; 
00030   gdouble ** data; 
00031 } slice_t;
00032 
00033 typedef struct {
00034   gint x, y, z;
00035   gboolean mid;
00036   gdouble d; 
00037 } tetra_vertex_t; 
00038 
00039 /* this helper is a lookup table for vertices */
00040 typedef struct {
00041   gint nx, ny; 
00042   GtsVertex ** vtop, ** vmid, **vbot;
00043 } helper_t ;
00044 
00045 typedef struct {
00046   GHashTable * vbot, * vtop;
00047 } helper_bcl ;
00048 
00049 
00050 static helper_t * init_helper (int nx, int ny) 
00051 {
00052   gint nxy = 4*nx*ny; 
00053   helper_t *retval = g_malloc0 (sizeof (helper_t));
00054 
00055   retval->nx = nx; 
00056   retval->ny = ny; 
00057   retval->vtop = g_malloc0 (sizeof (GtsVertex *)*nxy);
00058   retval->vmid = g_malloc0 (sizeof (GtsVertex *)*nxy);
00059   retval->vbot = g_malloc0 (sizeof (GtsVertex *)*nxy);
00060   return retval;
00061 }
00062 
00063 static helper_bcl * init_helper_bcl (void)
00064 {
00065   helper_bcl *retval = g_malloc0 (sizeof (helper_bcl));
00066 
00067   retval->vtop = g_hash_table_new (g_str_hash, g_str_equal);
00068   retval->vbot = g_hash_table_new (g_str_hash, g_str_equal);
00069   return retval;
00070 }
00071 
00072 static void free_helper (helper_t * h) 
00073 {
00074   g_free (h->vtop);
00075   g_free (h->vmid);
00076   g_free (h->vbot);
00077   g_free (h);
00078 }
00079 
00080 static void free_helper_bcl (helper_bcl * h) 
00081 {
00082   g_hash_table_destroy (h->vtop);
00083   g_hash_table_destroy (h->vbot);
00084   g_free (h);
00085 }
00086 
00087 /* move the vertices in the bottom slice to the top, and clear the
00088    other slices in the lookup tables */
00089 static void helper_advance (helper_t * h) 
00090 {
00091   GtsVertex ** help = h->vbot;
00092   h->vbot = h->vtop; 
00093   h->vtop = help;
00094   
00095   memset (h->vmid, 0, 4*sizeof(GtsVertex *) * h->nx * h->ny);
00096   memset (h->vbot, 0, 4*sizeof(GtsVertex *) * h->nx * h->ny);
00097 }
00098 
00099 static void helper_advance_bcl (helper_bcl * h) 
00100 {
00101   GHashTable * help = g_hash_table_new (g_str_hash, g_str_equal);
00102 
00103   g_hash_table_destroy (h->vbot);
00104   h->vbot = h->vtop;
00105   h->vtop = help;
00106 }
00107 
00108 /* find the zero-crossing of line through v1 and v2 and return the
00109    corresponding GtsVertex */
00110 static GtsVertex * get_vertex (gint mz, 
00111                                const tetra_vertex_t * v1, 
00112                                const tetra_vertex_t * v2, 
00113                                helper_t * help, 
00114                                GtsCartesianGrid * g,
00115                                GtsVertexClass * klass)
00116 {
00117   GtsVertex ** vertex; 
00118   gint x, y, index, idx2, z; 
00119   gdouble dx, dy, dz, d; 
00120 
00121   g_assert (v1->d - v2->d != 0.);
00122   
00123   dx = dy = dz = 0.0;
00124   d = v1->d/(v1->d - v2->d);
00125 
00126   index = 0;
00127   
00128   if (v1->x != v2->x) {
00129     index |= 1;
00130     dx = d;
00131   }
00132   
00133   if (v1->y != v2->y) {
00134     index |= 2;
00135     dy = d;
00136   }
00137   
00138   if (v1->z != v2->z) {
00139     dz = d;
00140   }
00141 
00142   x = v1->x;
00143   if (v1->x > v2->x) {  x = v2->x; dx = 1.0 - dx; }
00144   
00145   y = v1->y;
00146   if (v1->y > v2->y) {  y = v2->y; dy = 1.0 - dy;}
00147   
00148   z = v1->z;
00149   if (v1->z > v2->z) {  z = v2->z; dz = 1.0 - dz;}
00150   
00151   idx2 = 4 * ( x + y * help->nx ) + index;
00152   
00153   if (v1->z == v2->z)
00154     vertex = (mz == z) ? &help->vtop[idx2] : &help->vbot[idx2];
00155   else
00156     vertex = &help->vmid[idx2];
00157   
00158   if (mz != z && dz != 0.0) {
00159     fprintf(stderr, "%f \n", dz);
00160   }
00161   
00162   /* if vertex is not yet created, do it now */
00163   if (!*vertex)
00164     *vertex = gts_vertex_new (klass,
00165                               g->dx * ( x + dx) + g->x, 
00166                               g->dy * ( y + dy) + g->y, 
00167                               g->dz * ( z + dz) + g->z);
00168   
00169   return *vertex;
00170 }
00171 
00172 static GtsVertex * get_vertex_bcl (gint mz, 
00173                                    const tetra_vertex_t * v1, 
00174                                    const tetra_vertex_t * v2, 
00175                                    helper_bcl * help, 
00176                                    GtsCartesianGrid * g,
00177                                    GtsVertexClass * klass)
00178 {
00179   GtsVertex * v;
00180   GHashTable * table;
00181   gchar * s1, * s2, * hash;
00182   gdouble x1, x2, y1, y2, z1, z2, d;
00183 
00184   g_assert (v1->d - v2->d != 0.);
00185 
00186   /* first find correct hash table */  
00187   if ((v1->z > mz) && (v2->z > mz))
00188     table = help->vtop;
00189   else
00190     table = help->vbot;
00191 
00192   d = v1->d / (v1->d - v2->d);
00193 
00194   /* sort vertices */
00195   s1 = g_strdup_printf ("%d %d %d %d", v1->x, v1->y, v1->z, v1->mid);
00196   s2 = g_strdup_printf ("%d %d %d %d", v2->x, v2->y, v2->z, v2->mid);
00197 
00198   hash = (d == 0.0) ? g_strdup (s1) :
00199     (d == 1.0) ? g_strdup (s2) :
00200     (strcmp (s1, s2) < 0) ? g_strjoin (" ", s1, s2, NULL) :
00201     g_strjoin (" ", s2, s1, NULL);
00202 
00203   /* return existing vertex or make a new one */
00204   v = g_hash_table_lookup (table, hash);
00205   if (!v){
00206 
00207     x1 = g->dx * (v1->x + (v1->mid / 2.0)) + g->x;
00208     x2 = g->dx * (v2->x + (v2->mid / 2.0)) + g->x;
00209     y1 = g->dy * (v1->y + (v1->mid / 2.0)) + g->y;
00210     y2 = g->dy * (v2->y + (v2->mid / 2.0)) + g->y;
00211     z1 = g->dz * (v1->z + (v1->mid / 2.0)) + g->z;
00212     z2 = g->dz * (v2->z + (v2->mid / 2.0)) + g->z;
00213 
00214     v = gts_vertex_new (klass, x1 * (1.0 - d) + d * x2,
00215                         y1 * (1.0 - d) + d * y2,
00216                         z1 * (1.0 - d) + d * z2);
00217 
00218     g_hash_table_insert (table, g_strdup(hash), v);
00219   }
00220   g_free (s1);
00221   g_free (s2);
00222   g_free (hash);
00223 
00224   return v;
00225 }
00226 
00227 /* create an edge connecting the zero crossings of lines through a
00228    pair of vertices, or return an existing one */
00229 static GtsEdge * get_edge (GtsVertex * v1, GtsVertex * v2,
00230                            GtsEdgeClass * klass)
00231 {
00232   GtsSegment *s;
00233   GtsEdge *edge; 
00234   
00235   g_assert (v1);
00236   g_assert (v2);
00237   
00238   s = gts_vertices_are_connected (v1, v2);
00239   
00240   if (GTS_IS_EDGE (s))
00241     edge = GTS_EDGE(s);
00242   else
00243     edge = gts_edge_new (klass, v1, v2);
00244   return edge; 
00245 }
00246 
00247 static void add_face (GtsSurface * surface, 
00248                       const tetra_vertex_t * a1, const tetra_vertex_t * a2, 
00249                       const tetra_vertex_t * b1, const tetra_vertex_t * b2, 
00250                       const tetra_vertex_t * c1, const tetra_vertex_t * c2, 
00251                       gint rev, helper_t * help, 
00252                       gint z, GtsCartesianGrid * g)
00253 {
00254   GtsFace * t; 
00255   GtsEdge * e1, * e2, * e3;     
00256   GtsVertex * v1 = get_vertex (z, a1, a2, help, g, surface->vertex_class);
00257   GtsVertex * v2 = get_vertex (z, b1, b2, help, g, surface->vertex_class);
00258   GtsVertex * v3 = get_vertex (z, c1, c2, help, g, surface->vertex_class);
00259 
00260   g_assert (v1 != v2);
00261   g_assert (v2 != v3);
00262   g_assert (v1 != v3);
00263 
00264   if (!rev) {
00265     e1 = get_edge (v1, v2, surface->edge_class);
00266     e2 = get_edge (v2, v3, surface->edge_class);
00267     e3 = get_edge (v1, v3, surface->edge_class);
00268   } else {
00269     e1 = get_edge (v1, v3, surface->edge_class);
00270     e2 = get_edge (v2, v3, surface->edge_class);
00271     e3 = get_edge (v1, v2, surface->edge_class);        
00272   }
00273   
00274   t = gts_face_new (surface->face_class, e1, e2, e3);   
00275   gts_surface_add_face (surface, t);
00276 }
00277 
00278 static void add_face_bcl (GtsSurface * surface, 
00279                           const tetra_vertex_t * a1, 
00280                           const tetra_vertex_t * a2, 
00281                           const tetra_vertex_t * b1, 
00282                           const tetra_vertex_t * b2, 
00283                           const tetra_vertex_t * c1, 
00284                           const tetra_vertex_t * c2, 
00285                           gint rev, helper_bcl * help, 
00286                           gint z, GtsCartesianGrid * g)
00287 {
00288   GtsFace * t; 
00289   GtsEdge * e1, * e2, * e3;     
00290   GtsVertex * v1 = get_vertex_bcl (z, a1, a2, help, g, surface->vertex_class);
00291   GtsVertex * v2 = get_vertex_bcl (z, b1, b2, help, g, surface->vertex_class);
00292   GtsVertex * v3 = get_vertex_bcl (z, c1, c2, help, g, surface->vertex_class);
00293 
00294   if (v1 == v2 || v2 == v3 || v1 == v3)
00295     return;
00296 
00297   if (!rev) {
00298     e1 = get_edge (v1, v2, surface->edge_class);
00299     e2 = get_edge (v2, v3, surface->edge_class);
00300     e3 = get_edge (v1, v3, surface->edge_class);
00301   } else {
00302     e1 = get_edge (v1, v3, surface->edge_class);
00303     e2 = get_edge (v2, v3, surface->edge_class);
00304     e3 = get_edge (v1, v2, surface->edge_class);        
00305   }
00306 
00307   t = gts_face_new (surface->face_class, e1, e2, e3);   
00308   gts_surface_add_face (surface, t);
00309 }
00310 
00311 /* create a new slice of site nx \times ny */
00312 static slice_t * new_slice (gint nx, gint ny) 
00313 {
00314   gint x; 
00315   slice_t * retval = g_malloc (sizeof (slice_t));
00316 
00317   retval->data = g_malloc (nx*sizeof(gdouble *));
00318   retval->nx = nx;
00319   retval->ny = ny;  
00320   for (x = 0; x < nx; x++) 
00321     retval->data[x] = g_malloc (ny*sizeof (gdouble));
00322   return retval; 
00323 }
00324 
00325 /* initialize a slice with inival */
00326 static void slice_init (slice_t * slice, gdouble inival)
00327 {
00328   gint x, y; 
00329   
00330   g_assert (slice);
00331         
00332   for (x = 0; x < slice->nx; x++) 
00333     for (y = 0; y < slice->ny; y++)
00334       slice->data[x][y] = inival; 
00335 }
00336 
00337 /* free the memory of a slice */
00338 static void free_slice (slice_t * slice) 
00339 {
00340   gint x; 
00341         
00342   g_return_if_fail (slice != NULL);
00343         
00344   for (x = 0; x < slice->nx; x++) 
00345     g_free (slice->data[x]);  
00346   g_free (slice->data);
00347   g_free (slice);
00348 }
00349 
00350 static void analyze_tetrahedra (const tetra_vertex_t * a, 
00351                                 const tetra_vertex_t * b, 
00352                                 const tetra_vertex_t * c, 
00353                                 const tetra_vertex_t * d, 
00354                                 gint parity, GtsSurface * surface, 
00355                                 helper_t * help, 
00356                                 gint z, GtsCartesianGrid * g)
00357 {
00358   gint rev = parity; 
00359   gint code = 0; 
00360         
00361   if (a->d >= 0.) code |= 1; 
00362   if (b->d >= 0.) code |= 2; 
00363   if (c->d >= 0.) code |= 4; 
00364   if (d->d >= 0.) code |= 8;
00365                 
00366   switch (code) {
00367   case 15:
00368   case 0: return; /* all inside or outside */           
00369   
00370   case 14:rev = !parity;
00371   case  1:add_face (surface, a, b, a, d, a, c, rev, help, z, g);
00372           break;
00373   case 13:rev = ! parity;  
00374   case  2:add_face (surface, a, b, b, c, b, d, rev, help, z, g);
00375           break;
00376   case 12:rev = !parity;          
00377   case  3:add_face (surface, a, d, a, c, b, c, rev, help, z, g);
00378           add_face (surface, a, d, b, c, b, d, rev, help, z, g);
00379           break;
00380   case 11:rev = !parity;          
00381   case  4:add_face (surface, a, c, c, d, b, c, rev, help, z, g);
00382           break;
00383   case 10:rev = !parity;          
00384   case 5: add_face (surface, a, b, a, d, c, d, rev, help, z, g);
00385           add_face (surface, a, b, c, d, b, c, rev, help, z, g);
00386           break;        
00387   case  9:rev = !parity; 
00388   case  6:add_face (surface, a, b, a, c, c, d, rev, help, z, g);
00389           add_face (surface, a, b, c, d, b, d, rev, help, z, g);
00390           break;
00391   case  7:rev = !parity;
00392   case  8:add_face (surface, a, d, b, d, c, d, rev, help, z, g);
00393     break; 
00394   }
00395 }
00396 
00397 static void analyze_tetrahedra_bcl (const tetra_vertex_t * a, 
00398                                     const tetra_vertex_t * b, 
00399                                     const tetra_vertex_t * c, 
00400                                     const tetra_vertex_t * d, 
00401                                     GtsSurface * surface, 
00402                                     helper_bcl * help, 
00403                                     gint z, GtsCartesianGrid * g)
00404 {
00405   gint rev = 0;
00406   gint code = 0; 
00407         
00408   if (a->d >= 0.) code |= 1; 
00409   if (b->d >= 0.) code |= 2; 
00410   if (c->d >= 0.) code |= 4; 
00411   if (d->d >= 0.) code |= 8;
00412 
00413   switch (code) {
00414   case 15:
00415   case 0: return; /* all inside or outside */           
00416 
00417   case 14:rev = !rev;
00418   case  1:add_face_bcl (surface, a, b, a, d, a, c, rev, help, z, g);
00419           break;
00420   case 13:rev = !rev;  
00421   case  2:add_face_bcl (surface, a, b, b, c, b, d, rev, help, z, g);
00422           break;
00423   case 12:rev = !rev;     
00424   case  3:add_face_bcl (surface, a, d, a, c, b, c, rev, help, z, g);
00425           add_face_bcl (surface, a, d, b, c, b, d, rev, help, z, g);
00426           break;
00427   case 11:rev = !rev;     
00428   case  4:add_face_bcl (surface, a, c, c, d, b, c, rev, help, z, g);
00429           break;
00430   case 10:rev = !rev;     
00431   case 5: add_face_bcl (surface, a, b, a, d, c, d, rev, help, z, g);
00432           add_face_bcl (surface, a, b, c, d, b, c, rev, help, z, g);
00433           break;        
00434   case  9:rev = !rev; 
00435   case  6:add_face_bcl (surface, a, b, a, c, c, d, rev, help, z, g);
00436           add_face_bcl (surface, a, b, c, d, b, d, rev, help, z, g);
00437           break;
00438   case  7:rev = !rev;
00439   case  8:add_face_bcl (surface, a, d, b, d, c, d, rev, help, z, g);
00440     break;
00441   }
00442 }
00443 
00444 static void  iso_slice_evaluate (slice_t * s1, slice_t * s2, 
00445                                  GtsCartesianGrid g, 
00446                                  gint z, GtsSurface * surface, helper_t * help)
00447 {
00448   gint x,y; 
00449   tetra_vertex_t v0, v1, v2, v3, v4, v5, v6, v7; 
00450   gdouble ** s1p = s1->data; 
00451   gdouble ** s2p = s2->data; 
00452         
00453   for (y = 0; y < g.ny-1; y++)
00454     for (x = 0; x < g.nx-1; x++) {
00455       gint parity = (((x ^ y) ^ z) & 1);
00456       
00457       v0.x = x  ; v0.y = y  ; v0.z = z  ; v0.mid = FALSE; v0.d = s1p[x  ][y  ];
00458       v1.x = x  ; v1.y = y+1; v1.z = z  ; v1.mid = FALSE; v1.d = s1p[x  ][y+1];
00459       v2.x = x+1; v2.y = y  ; v2.z = z  ; v2.mid = FALSE; v2.d = s1p[x+1][y  ];
00460       v3.x = x+1; v3.y = y+1; v3.z = z  ; v3.mid = FALSE; v3.d = s1p[x+1][y+1];
00461       v4.x = x  ; v4.y = y  ; v4.z = z+1; v4.mid = FALSE; v4.d = s2p[x  ][y  ];
00462       v5.x = x  ; v5.y = y+1; v5.z = z+1; v5.mid = FALSE; v5.d = s2p[x  ][y+1];
00463       v6.x = x+1; v6.y = y  ; v6.z = z+1; v6.mid = FALSE; v6.d = s2p[x+1][y  ];
00464       v7.x = x+1; v7.y = y+1; v7.z = z+1; v7.mid = FALSE; v7.d = s2p[x+1][y+1];
00465       
00466       if (parity == 0) {
00467         analyze_tetrahedra (&v0, &v1, &v2, &v4, parity, surface, help, z, &g);
00468         analyze_tetrahedra (&v7, &v1, &v4, &v2, parity, surface, help, z, &g);
00469         analyze_tetrahedra (&v1, &v7, &v3, &v2, parity, surface, help, z, &g);
00470         analyze_tetrahedra (&v1, &v7, &v4, &v5, parity, surface, help, z, &g);
00471         analyze_tetrahedra (&v2, &v6, &v4, &v7, parity, surface, help, z, &g);
00472       }else{
00473         analyze_tetrahedra (&v4, &v5, &v6, &v0, parity, surface, help, z, &g);
00474         analyze_tetrahedra (&v3, &v5, &v0, &v6, parity, surface, help, z, &g);
00475         analyze_tetrahedra (&v5, &v3, &v7, &v6, parity, surface, help, z, &g);
00476         analyze_tetrahedra (&v5, &v3, &v0, &v1, parity, surface, help, z, &g);
00477         analyze_tetrahedra (&v6, &v2, &v0, &v3, parity, surface, help, z, &g);
00478       }
00479     }
00480 }
00481 
00482 static void  iso_slice_evaluate_bcl (slice_t * s1, slice_t * s2, slice_t * s3,
00483                                      GtsCartesianGrid g, 
00484                                      gint z, GtsSurface * surface, 
00485                                      helper_bcl * help)
00486 {
00487   gint x,y; 
00488   tetra_vertex_t v0, v1, v2, v3, v4, v5, v6, v7, v8, v9, w0; 
00489   gdouble ** s1p = s1->data;
00490   gdouble ** s2p = s2->data;
00491   gdouble ** s3p = s3->data;
00492         
00493   for (y = 0; y < g.ny-2; y++)
00494     for (x = 0; x < g.nx-2; x++) {
00495       v0.x = x  ; v0.y = y  ; v0.z = z  ; v0.mid = TRUE;
00496       v0.d = (s1p[x  ][y  ] + s2p[x  ][y  ] +
00497               s1p[x  ][y+1] + s2p[x  ][y+1] +
00498               s1p[x+1][y  ] + s2p[x+1][y  ] +
00499               s1p[x+1][y+1] + s2p[x+1][y+1])/8.0;
00500 
00501       v1.x = x+1; v1.y = y  ; v1.z = z  ; v1.mid = TRUE;
00502       v1.d = (s1p[x+1][y  ] + s2p[x+1][y  ] +
00503               s1p[x+1][y+1] + s2p[x+1][y+1] +
00504               s1p[x+2][y  ] + s2p[x+2][y  ] +
00505               s1p[x+2][y+1] + s2p[x+2][y+1])/8.0;
00506 
00507       v2.x = x  ; v2.y = y+1; v2.z = z  ; v2.mid = TRUE;
00508       v2.d = (s1p[x  ][y+1] + s2p[x  ][y+1] +
00509               s1p[x  ][y+2] + s2p[x  ][y+2] +
00510               s1p[x+1][y+1] + s2p[x+1][y+1] +
00511               s1p[x+1][y+2] + s2p[x+1][y+2])/8.0;
00512 
00513       v3.x = x  ; v3.y = y  ; v3.z = z+1; v3.mid = TRUE;
00514       v3.d = (s2p[x  ][y  ] + s3p[x  ][y  ] +
00515               s2p[x  ][y+1] + s3p[x  ][y+1] +
00516               s2p[x+1][y  ] + s3p[x+1][y  ] +
00517               s2p[x+1][y+1] + s3p[x+1][y+1])/8.0;
00518 
00519       v4.x = x+1; v4.y = y  ; v4.z = z  ; v4.mid = FALSE; v4.d = s1p[x+1][y  ];
00520       v5.x = x  ; v5.y = y+1; v5.z = z  ; v5.mid = FALSE; v5.d = s1p[x  ][y+1];
00521       v6.x = x+1; v6.y = y+1; v6.z = z  ; v6.mid = FALSE; v6.d = s1p[x+1][y+1];
00522       v7.x = x+1; v7.y = y  ; v7.z = z+1; v7.mid = FALSE; v7.d = s2p[x+1][y  ];
00523       v8.x = x  ; v8.y = y+1; v8.z = z+1; v8.mid = FALSE; v8.d = s2p[x  ][y+1];
00524       v9.x = x+1; v9.y = y+1; v9.z = z+1; v9.mid = FALSE; v9.d = s2p[x+1][y+1];
00525       w0.x = x  ; w0.y = y  ; w0.z = z+1; w0.mid = FALSE; w0.d = s2p[x  ][y  ];
00526 
00527       analyze_tetrahedra_bcl (&v0, &v9, &v6, &v1, surface, help, z, &g);
00528       analyze_tetrahedra_bcl (&v0, &v6, &v4, &v1, surface, help, z, &g);
00529       analyze_tetrahedra_bcl (&v0, &v4, &v7, &v1, surface, help, z, &g);
00530       analyze_tetrahedra_bcl (&v0, &v7, &v9, &v1, surface, help, z, &g);
00531       analyze_tetrahedra_bcl (&v0, &v5, &v6, &v2, surface, help, z, &g);
00532       analyze_tetrahedra_bcl (&v0, &v6, &v9, &v2, surface, help, z, &g);
00533       analyze_tetrahedra_bcl (&v0, &v9, &v8, &v2, surface, help, z, &g);
00534       analyze_tetrahedra_bcl (&v0, &v8, &v5, &v2, surface, help, z, &g);
00535       analyze_tetrahedra_bcl (&v0, &v8, &v9, &v3, surface, help, z, &g);
00536       analyze_tetrahedra_bcl (&v0, &v9, &v7, &v3, surface, help, z, &g);
00537       analyze_tetrahedra_bcl (&v0, &v7, &w0, &v3, surface, help, z, &g);
00538       analyze_tetrahedra_bcl (&v0, &w0, &v8, &v3, surface, help, z, &g);
00539     }
00540 }
00541 
00542 /*  copy src into dest by stripping off the iso value and leave out
00543     the boundary (which should be G_MINDOUBLE) */
00544 static void copy_to_bounded (slice_t * dest, slice_t * src, 
00545                              gdouble iso, gdouble fill)
00546 {
00547   gint x,y; 
00548   gdouble * src_ptr;
00549   gdouble * dest_ptr = dest->data[0];
00550   
00551   g_assert(dest->ny == src->ny + 2);
00552   g_assert(dest->nx == src->nx + 2);
00553         
00554   for (y = 0; y < dest->ny; ++y, ++dest_ptr)
00555     *dest_ptr = fill; 
00556         
00557   for (x = 1; x < src->nx - 1; ++x) {
00558     dest_ptr = dest->data[x];
00559     src_ptr = src->data[x-1];
00560     *dest_ptr++ = fill;
00561     for (y = 0; y < src->ny; ++y, ++dest_ptr, ++src_ptr)
00562       *dest_ptr  = *src_ptr - iso;
00563     *dest_ptr++ = fill; 
00564   }
00565   
00566   dest_ptr = dest->data[y];
00567   
00568   for (y = 0; y < dest->ny; ++y, ++dest_ptr)
00569     *dest_ptr = fill; 
00570 }
00571 
00572 static void iso_sub (slice_t * s, gdouble iso)
00573 {
00574   gint x,y; 
00575 
00576   for (x = 0; x < s->nx; ++x) {
00577     gdouble *ptr = s->data[x];
00578 
00579     for (y = 0; y < s->ny; ++y, ++ptr)
00580       *ptr -= iso; 
00581   }
00582 }
00583 
00584 
00603 void gts_isosurface_tetra_bounded (GtsSurface * surface,
00604                                    GtsCartesianGrid g,
00605                                    GtsIsoCartesianFunc f,
00606                                    gpointer data,
00607                                    gdouble iso)
00608 {
00609   slice_t *slice1, *slice2, *transfer_slice; 
00610   GtsCartesianGrid g_intern = g; 
00611   helper_t *helper;
00612   gint z; 
00613         
00614   g_return_if_fail (surface != NULL);
00615   g_return_if_fail (f != NULL);
00616   g_return_if_fail (g.nx > 1);
00617   g_return_if_fail (g.ny > 1);
00618   g_return_if_fail (g.nz > 1);
00619 
00620   /* create the helper slices */
00621   slice1 = new_slice (g.nx + 2, g.ny + 2);
00622   slice2 = new_slice (g.nx + 2, g.ny + 2);
00623         
00624   /*  initialize the first slice as OUTSIDE */
00625   slice_init (slice1, -1.0);
00626         
00627   /* create a slice of the original image size */
00628   transfer_slice = new_slice (g.nx, g.ny);
00629         
00630   /* adapt the parameters to our enlarged image */
00631   g_intern.x -= g.dx;
00632   g_intern.y -= g.dy; 
00633   g_intern.z -= g.dz;
00634   g_intern.nx = g.nx + 2;
00635   g_intern.ny = g.ny + 2;       
00636   g_intern.nz = g.nz;
00637         
00638   /* create the helper for vertex-lookup */
00639   helper = init_helper (g_intern.nx, g_intern.ny);
00640         
00641   /* go slicewise through the data */
00642   z = 0; 
00643   while (z < g.nz) {
00644     slice_t * hs; 
00645     
00646     /* request slice */
00647     f (transfer_slice->data, g, z, data);
00648     g.z += g.dz; 
00649     
00650     /* copy slice in enlarged image and mark the border as OUTSIDE */
00651     copy_to_bounded (slice2, transfer_slice, iso, -1.);
00652     
00653     /* triangulate */
00654     iso_slice_evaluate (slice1, slice2, g_intern, z, surface, helper);
00655     
00656     /* switch the input slices */
00657     hs = slice1; slice1 = slice2; slice2 = hs; 
00658     
00659     /* switch the vertex lookup tables */
00660     helper_advance(helper);
00661     ++z; 
00662   }
00663   
00664   /* initialize the last slice as OUTSIDE */
00665   slice_init (slice2, - 1.0);
00666                 
00667   /* close the object */
00668   iso_slice_evaluate(slice1, slice2, g_intern, z, surface, helper);
00669   
00670   free_helper (helper);
00671   free_slice (slice1);
00672   free_slice (slice2);
00673   free_slice (transfer_slice);  
00674 }
00675 
00693 void gts_isosurface_tetra (GtsSurface * surface,
00694                            GtsCartesianGrid g,
00695                            GtsIsoCartesianFunc f,
00696                            gpointer data,
00697                            gdouble iso)
00698 {
00699   slice_t *slice1, *slice2; 
00700   helper_t *helper;
00701   gint z; 
00702   GtsCartesianGrid g_internal;
00703   
00704   g_return_if_fail (surface != NULL);
00705   g_return_if_fail (f != NULL);
00706   g_return_if_fail (g.nx > 1);
00707   g_return_if_fail (g.ny > 1);
00708   g_return_if_fail (g.nz > 1);
00709 
00710   memcpy (&g_internal, &g, sizeof (GtsCartesianGrid));
00711         
00712   /* create the helper slices */
00713   slice1 = new_slice (g.nx, g.ny);
00714   slice2 = new_slice (g.nx, g.ny);
00715   
00716   /* create the helper for vertex-lookup */
00717   helper = init_helper (g.nx, g.ny);
00718         
00719   z = 0;
00720   f (slice1->data, g, z, data);
00721   iso_sub (slice1, iso); 
00722   
00723   z++; 
00724   g.z += g.dz;
00725   
00726   /* go slicewise through the data */
00727   while (z < g.nz) {
00728     slice_t * hs; 
00729     
00730     /* request slice */
00731     f (slice2->data, g, z, data);
00732     iso_sub (slice2, iso);
00733      
00734     g.z += g.dz;
00735     
00736     /* triangulate */
00737     iso_slice_evaluate (slice1, slice2, g_internal, z-1, surface, helper);
00738     
00739     /* switch the input slices */
00740     hs = slice1; slice1 = slice2; slice2 = hs; 
00741     
00742     /* switch the vertex lookup tables */
00743     helper_advance (helper);
00744     
00745     ++z; 
00746   }
00747 
00748   free_helper(helper);
00749   free_slice(slice1);
00750   free_slice(slice2);   
00751 }
00752 
00773 void gts_isosurface_tetra_bcl (GtsSurface * surface,
00774                                GtsCartesianGrid g,
00775                                GtsIsoCartesianFunc f,
00776                                gpointer data,
00777                                gdouble iso)
00778 {
00779   slice_t *slice1, *slice2, *slice3;
00780   helper_bcl *helper;
00781   gint z; 
00782   GtsCartesianGrid g_internal;
00783   
00784   g_return_if_fail (surface != NULL);
00785   g_return_if_fail (f != NULL);
00786   g_return_if_fail (g.nx > 1);
00787   g_return_if_fail (g.ny > 1);
00788   g_return_if_fail (g.nz > 1);
00789 
00790   memcpy (&g_internal, &g, sizeof (GtsCartesianGrid));
00791         
00792   /* create the helper slices */
00793   slice1 = new_slice (g.nx, g.ny);
00794   slice2 = new_slice (g.nx, g.ny);
00795   slice3 = new_slice (g.nx, g.ny);
00796   
00797   /* create the helper for vertex-lookup */
00798   helper = init_helper_bcl ();
00799         
00800   z = 0;
00801   f (slice1->data, g, z, data);
00802   iso_sub (slice1, iso); 
00803 
00804   z++; 
00805   g.z += g.dz;
00806 
00807   f (slice2->data, g, z, data);
00808   iso_sub (slice1, iso); 
00809   
00810   z++; 
00811   g.z += g.dz;
00812   
00813   /* go slicewise through the data */
00814   while (z < g.nz) {
00815     slice_t * hs; 
00816     
00817     /* request slice */
00818     f (slice3->data, g, z, data);
00819     iso_sub (slice3, iso);
00820      
00821     g.z += g.dz;
00822     
00823     /* triangulate */
00824     iso_slice_evaluate_bcl (slice1, slice2, slice3, g_internal, z-2, 
00825                             surface, helper);
00826     
00827     /* switch the input slices */
00828     hs = slice1; slice1 = slice2; slice2 = slice3; slice3 = hs;
00829     
00830     /* switch the vertex lookup tables */
00831     helper_advance_bcl (helper);
00832     
00833     ++z; 
00834   }
00835 
00836   free_helper_bcl(helper);
00837   free_slice(slice1);
00838   free_slice(slice2);   
00839   free_slice(slice3);   
00840 }