gts/iso.c

   1 /* GTS - Library for the manipulation of triangulated surfaces
   2  * Copyright (C) 1999 Stéphane Popinet
   3  *
   4  * This library is free software; you can redistribute it and/or
   5  * modify it under the terms of the GNU Library General Public
   6  * License as published by the Free Software Foundation; either
   7  * version 2 of the License, or (at your option) any later version.
   8  *
   9  * This library is distributed in the hope that it will be useful,
  10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12  * Library General Public License for more details.
  13  *
  14  * You should have received a copy of the GNU Library General Public
  15  * License along with this library; if not, write to the
  16  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  17  * Boston, MA 02111-1307, USA.
  18  */
  19
  20 #include "gts.h"
  21
  22 typedef enum { LEFT = 0, RIGHT = 1 } Orientation;
  23
  24 typedef struct {
  25   GtsVertex * v;
  26   Orientation orientation;
  27 } OrientedVertex;
  28
  29 struct _GtsIsoSlice {
  30   OrientedVertex *** vertices;
  31   guint nx, ny;
  32 };
  33
  34 /* coordinates of the edges of the cube (see doc/isocube.fig) */
  35 static guint c[12][4] = {
  36   {0, 0, 0, 0}, {0, 0, 0, 1}, {0, 0, 1, 1}, {0, 0, 1, 0},
  37   {1, 0, 0, 0}, {1, 0, 0, 1}, {1, 1, 0, 1}, {1, 1, 0, 0},
  38   {2, 0, 0, 0}, {2, 1, 0, 0}, {2, 1, 1, 0}, {2, 0, 1, 0}};
  39
  40 /* first index is the edge number, second index is the edge orientation
  41    (RIGHT or LEFT), third index are the edges which this edge may connect to
  42    in order */
  43 static guint edge[12][2][3] = {
  44   {{9, 1, 8}, {4, 3, 7}},   /* 0 */
  45   {{6, 2, 5}, {8, 0, 9}},   /* 1 */
  46   {{10, 3, 11}, {5, 1, 6}}, /* 2 */
  47   {{7, 0, 4}, {11, 2, 10}}, /* 3 */
  48   {{3, 7, 0}, {8, 5, 11}},  /* 4 */
  49   {{11, 4, 8}, {1, 6, 2}},  /* 5 */
  50   {{2, 5, 1}, {9, 7, 10}},  /* 6 */
  51   {{10, 6, 9}, {0, 4, 3}},  /* 7 */
  52   {{5, 11, 4}, {0, 9, 1}},  /* 8 */
  53   {{1, 8, 0}, {7, 10, 6}},  /* 9 */
  54   {{6, 9, 7}, {3, 11, 2}},  /* 10 */
  55   {{2, 10, 3}, {4, 8, 5}}   /* 11 */
  56 };
  57
  58 static void ** malloc2D (guint nx, guint ny, gulong size)
  59 {
  60   void ** m = g_malloc (nx*sizeof (void *));
  61   guint i;
  62
  63   for (i = 0; i < nx; i++)
  64     m[i] = g_malloc0 (ny*size);
  65
  66   return m;
  67 }
  68
  69 static void free2D (void ** m, guint nx)
  70 {
  71   guint i;
  72
  73   g_return_if_fail (m != NULL);
  74
  75   for (i = 0; i < nx; i++)
  76     g_free (m[i]);
  77   g_free (m);
  78 }
  79
  80 /**
  81  * gts_grid_plane_new:
  82  * @nx:
  83  * @ny:
  84  *
  85  * Returns:
  86  */
  87 GtsGridPlane * gts_grid_plane_new (guint nx, guint ny)
  88 {
  89   GtsGridPlane * g = g_malloc (sizeof (GtsGridPlane));
  90
  91   g->p = (GtsPoint **) malloc2D (nx, ny, sizeof (GtsPoint));
  92   g->nx = nx;
  93   g->ny = ny;
  94
  95   return g;
  96 }
  97
  98 /**
  99  * gts_grid_plane_destroy:
 100  * @g:
 101  *
 102  */
 103 void gts_grid_plane_destroy (GtsGridPlane * g)
 104 {
 105   g_return_if_fail (g != NULL);
 106
 107   free2D ((void **) g->p, g->nx);
 108   g_free (g);
 109 }
 110
 111 /**
 112  * gts_iso_slice_new:
 113  * @nx: number of vertices in the x direction.
 114  * @ny: number of vertices in the y direction.
 115  *
 116  * Returns: a new #GtsIsoSlice.
 117  */
 118 GtsIsoSlice * gts_iso_slice_new (guint nx, guint ny)
 119 {
 120   GtsIsoSlice * slice;
 121
 122   g_return_val_if_fail (nx > 1, NULL);
 123   g_return_val_if_fail (ny > 1, NULL);
 124
 125   slice = g_malloc (sizeof (GtsIsoSlice));
 126
 127   slice->vertices = g_malloc (3*sizeof (OrientedVertex **));
 128   slice->vertices[0] =
 129     (OrientedVertex **) malloc2D (nx, ny, sizeof (OrientedVertex));
 130   slice->vertices[1] =
 131     (OrientedVertex **) malloc2D (nx - 1, ny, sizeof (OrientedVertex));
 132   slice->vertices[2] =
 133     (OrientedVertex **) malloc2D (nx, ny - 1, sizeof (OrientedVertex));
 134   slice->nx = nx;
 135   slice->ny = ny;
 136
 137   return slice;
 138 }
 139
 140 /**
 141  * gts_iso_slice_fill:
 142  * @slice: a #GtsIsoSlice.
 143  * @plane1: a #GtsGridPlane.
 144  * @plane2: another #GtsGridPlane.
 145  * @f1: values of the function corresponding to @plane1.
 146  * @f2: values of the function corresponding to @plane2.
 147  * @iso: isosurface value.
 148  * @klass: a #GtsVertexClass or one of its descendant to be used for the
 149  * new vertices.
 150  *
 151  * Fill @slice with the coordinates of the vertices defined by
 152  * f1 (x,y,z) = @iso and f2 (x, y, z) = @iso.
 153  */
 154 void gts_iso_slice_fill (GtsIsoSlice * slice,
 155                          GtsGridPlane * plane1,
 156                          GtsGridPlane * plane2,
 157                          gdouble ** f1,
 158                          gdouble ** f2,
 159                          gdouble iso,
 160                          GtsVertexClass * klass)
 161 {
 162   OrientedVertex *** vertices;
 163   GtsPoint ** p1, ** p2 = NULL;
 164   guint i, j, nx, ny;
 165
 166   g_return_if_fail (slice != NULL);
 167   g_return_if_fail (plane1 != NULL);
 168   g_return_if_fail (f1 != NULL);
 169   g_return_if_fail (f2 == NULL || plane2 != NULL);
 170
 171   p1 = plane1->p;
 172   if (plane2)
 173     p2 = plane2->p;
 174   vertices = slice->vertices;
 175   nx = slice->nx;
 176   ny = slice->ny;
 177
 178   if (f2)
 179     for (i = 0; i < nx; i++)
 180       for (j = 0; j < ny; j++) {
 181         gdouble v1 = f1[i][j] - iso;
 182         gdouble v2 = f2[i][j] - iso;
 183         if ((v1 >= 0. && v2 < 0.) || (v1 < 0. && v2 >= 0.)) {
 184           gdouble c2 = v1/(v1 - v2), c1 = 1. - c2;
 185           vertices[0][i][j].v =
 186             gts_vertex_new (klass,
 187                             c1*p1[i][j].x + c2*p2[i][j].x,
 188                             c1*p1[i][j].y + c2*p2[i][j].y,
 189                             c1*p1[i][j].z + c2*p2[i][j].z);
 190           vertices[0][i][j].orientation = v2 >= 0. ? RIGHT : LEFT;
 191         }
 192         else
 193           vertices[0][i][j].v = NULL;
 194       }
 195   for (i = 0; i < nx - 1; i++)
 196     for (j = 0; j < ny; j++) {
 197       gdouble v1 = f1[i][j] - iso;
 198       gdouble v2 = f1[i+1][j] - iso;
 199       if ((v1 >= 0. && v2 < 0.) || (v1 < 0. && v2 >= 0.)) {
 200         gdouble c2 = v1/(v1 - v2), c1 = 1. - c2;
 201         vertices[1][i][j].v =
 202           gts_vertex_new (klass,
 203                           c1*p1[i][j].x + c2*p1[i+1][j].x,
 204                           c1*p1[i][j].y + c2*p1[i+1][j].y,
 205                           c1*p1[i][j].z + c2*p1[i+1][j].z);
 206         vertices[1][i][j].orientation = v2 >= 0. ? RIGHT : LEFT;
 207       }
 208       else
 209         vertices[1][i][j].v = NULL;
 210     }
 211   for (i = 0; i < nx; i++)
 212     for (j = 0; j < ny - 1; j++) {
 213       gdouble v1 = f1[i][j] - iso;
 214       gdouble v2 = f1[i][j+1] - iso;
 215       if ((v1 >= 0. && v2 < 0.) || (v1 < 0. && v2 >= 0.)) {
 216         gdouble c2 = v1/(v1 - v2), c1 = 1. - c2;
 217         vertices[2][i][j].v =
 218           gts_vertex_new (klass,
 219                           c1*p1[i][j].x + c2*p1[i][j+1].x,
 220                           c1*p1[i][j].y + c2*p1[i][j+1].y,
 221                           c1*p1[i][j].z + c2*p1[i][j+1].z);
 222         vertices[2][i][j].orientation = v2 >= 0. ? RIGHT : LEFT;
 223       }
 224       else
 225         vertices[2][i][j].v = NULL;
 226     }
 227 }
 228
 229 /**
 230  * gts_iso_slice_fill_cartesian:
 231  * @slice: a #GtsIsoSlice.
 232  * @g: a #GtsCartesianGrid.
 233  * @f1: values of the function for plane z = @g.z.
 234  * @f2: values of the function for plane z = @g.z + @g.dz.
 235  * @iso: isosurface value.
 236  * @klass: a #GtsVertexClass.
 237  *
 238  * Fill @slice with the coordinates of the vertices defined by
 239  * f1 (x,y,z) = @iso and f2 (x, y, z) = @iso.
 240  */
 241 void gts_iso_slice_fill_cartesian (GtsIsoSlice * slice,
 242                                    GtsCartesianGrid g,
 243                                    gdouble ** f1,
 244                                    gdouble ** f2,
 245                                    gdouble iso,
 246                                    GtsVertexClass * klass)
 247 {
 248   OrientedVertex *** vertices;
 249   guint i, j;
 250   gdouble x, y;
 251
 252   g_return_if_fail (slice != NULL);
 253   g_return_if_fail (f1 != NULL);
 254
 255   vertices = slice->vertices;
 256
 257   if (f2)
 258     for (i = 0, x = g.x; i < g.nx; i++, x += g.dx)
 259       for (j = 0, y = g.y; j < g.ny; j++, y += g.dy) {
 260         gdouble v1 = f1[i][j] - iso;
 261         gdouble v2 = f2[i][j] - iso;
 262         if ((v1 >= 0. && v2 < 0.) || (v1 < 0. && v2 >= 0.)) {
 263           vertices[0][i][j].v =
 264             gts_vertex_new (klass,
 265                             x, y, g.z + g.dz*v1/(v1 - v2));
 266           vertices[0][i][j].orientation = v2 >= 0. ? RIGHT : LEFT;
 267         }
 268         else
 269           vertices[0][i][j].v = NULL;
 270       }
 271   for (i = 0, x = g.x; i < g.nx - 1; i++, x += g.dx)
 272     for (j = 0, y = g.y; j < g.ny; j++, y += g.dy) {
 273       gdouble v1 = f1[i][j] - iso;
 274       gdouble v2 = f1[i+1][j] - iso;
 275       if ((v1 >= 0. && v2 < 0.) || (v1 < 0. && v2 >= 0.)) {
 276         vertices[1][i][j].v =
 277           gts_vertex_new (klass, x + g.dx*v1/(v1 - v2), y, g.z);
 278         vertices[1][i][j].orientation = v2 >= 0. ? RIGHT : LEFT;
 279       }
 280       else
 281         vertices[1][i][j].v = NULL;
 282     }
 283   for (i = 0, x = g.x; i < g.nx; i++, x += g.dx)
 284     for (j = 0, y = g.y; j < g.ny - 1; j++, y += g.dy) {
 285       gdouble v1 = f1[i][j] - iso;
 286       gdouble v2 = f1[i][j+1] - iso;
 287       if ((v1 >= 0. && v2 < 0.) || (v1 < 0. && v2 >= 0.)) {
 288         vertices[2][i][j].v =
 289           gts_vertex_new (klass, x, y + g.dy*v1/(v1 - v2), g.z);
 290         vertices[2][i][j].orientation = v2 >= 0. ? RIGHT : LEFT;
 291       }
 292       else
 293         vertices[2][i][j].v = NULL;
 294     }
 295 }
 296
 297 /**
 298  * gts_iso_slice_destroy:
 299  * @slice: a #GtsIsoSlice.
 300  *
 301  * Free all memory allocated for @slice.
 302  */
 303 void gts_iso_slice_destroy (GtsIsoSlice * slice)
 304 {
 305   g_return_if_fail (slice != NULL);
 306
 307   free2D ((void **) slice->vertices[0], slice->nx);
 308   free2D ((void **) slice->vertices[1], slice->nx - 1);
 309   free2D ((void **) slice->vertices[2], slice->nx);
 310   g_free (slice->vertices);
 311   g_free (slice);
 312 }
 313
 314 /**
 315  * gts_isosurface_slice:
 316  * @slice1: a #GtsIsoSlice.
 317  * @slice2: another #GtsIsoSlice.
 318  * @surface: a #GtsSurface.
 319  *
 320  * Given two successive slices @slice1 and @slice2 link their vertices with
 321  * segments and triangles which are added to @surface.
 322  */
 323 void gts_isosurface_slice (GtsIsoSlice * slice1,
 324                            GtsIsoSlice * slice2,
 325                            GtsSurface * surface)
 326 {
 327   guint j, k, l, nx, ny;
 328   OrientedVertex *** vertices[2];
 329   GtsVertex * va[12];
 330
 331   g_return_if_fail (slice1 != NULL);
 332   g_return_if_fail (slice2 != NULL);
 333   g_return_if_fail (surface != NULL);
 334   g_return_if_fail (slice1->nx == slice2->nx && slice1->ny == slice2->ny);
 335
 336   vertices[0] = slice1->vertices;
 337   vertices[1] = slice2->vertices;
 338   nx = slice1->nx;
 339   ny = slice1->ny;
 340
 341   /* link vertices with segments and triangles */
 342   for (j = 0; j < nx - 1; j++)
 343     for (k = 0; k < ny - 1; k++) {
 344       gboolean cube_is_cut = FALSE;
 345       for (l = 0; l < 12; l++) {
 346         guint nv = 0, e = l;
 347         OrientedVertex ov =
 348           vertices[c[e][1]][c[e][0]][j + c[e][2]][k + c[e][3]];
 349         while (ov.v && !GTS_OBJECT (ov.v)->reserved) {
 350           guint m = 0, * ne = edge[e][ov.orientation];
 351           va[nv++] = ov.v;
 352           GTS_OBJECT (ov.v)->reserved = surface;
 353           ov.v = NULL;
 354           while (m < 3 && !ov.v) {
 355             e = ne[m++];
 356             ov = vertices[c[e][1]][c[e][0]][j + c[e][2]][k + c[e][3]];
 357           }
 358         }
 359         /* create edges and faces */
 360         if (nv > 2) {
 361           GtsEdge * e1, * e2, * e3;
 362           guint m;
 363           if (!(e1 = GTS_EDGE (gts_vertices_are_connected (va[0], va[1]))))
 364             e1 = gts_edge_new (surface->edge_class, va[0], va[1]);
 365           for (m = 1; m < nv - 1; m++) {
 366             if (!(e2 = GTS_EDGE (gts_vertices_are_connected (va[m], va[m+1]))))
 367               e2 = gts_edge_new (surface->edge_class, va[m], va[m+1]);
 368             if (!(e3 = GTS_EDGE (gts_vertices_are_connected (va[m+1], va[0]))))
 369               e3 = gts_edge_new (surface->edge_class, va[m+1], va[0]);
 370             gts_surface_add_face (surface,
 371                                   gts_face_new (surface->face_class,
 372                                                 e1, e2, e3));
 373             e1 = e3;
 374           }
 375         }
 376         if (nv > 0)
 377           cube_is_cut = TRUE;
 378       }
 379       if (cube_is_cut)
 380         for (l = 0; l < 12; l++) {
 381           GtsVertex * v =
 382             vertices[c[l][1]][c[l][0]][j + c[l][2]][k + c[l][3]].v;
 383           if (v)
 384             GTS_OBJECT (v)->reserved = NULL;
 385         }
 386     }
 387 }
 388
 389 #define SWAP(s1, s2, tmp) (tmp = s1, s1 = s2, s2 = tmp)
 390
 391 /**
 392  * gts_isosurface_cartesian:
 393  * @surface: a #GtsSurface.
 394  * @g: a #GtsCartesianGrid.
 395  * @f: a #GtsIsoCartesianFunc.
 396  * @data: user data to be passed to @f.
 397  * @iso: isosurface value.
 398  *
 399  * Adds to @surface new faces defining the isosurface f(x,y,z) = @iso. By
 400  * convention, the normals to the surface are pointing toward the positive
 401  * values of f(x,y,z) - @iso.
 402  *
 403  * The user function @f is called successively for each value of the z
 404  * coordinate defined by @g. It must fill the corresponding (x,y) plane with
 405  * the values of the function for which the isosurface is to be computed.
 406  */
 407 void gts_isosurface_cartesian (GtsSurface * surface,
 408                                GtsCartesianGrid g,
 409                                GtsIsoCartesianFunc f,
 410                                gpointer data,
 411                                gdouble iso)
 412 {
 413   void * tmp;
 414   gdouble ** f1, ** f2;
 415   GtsIsoSlice * slice1, * slice2;
 416   guint i;
 417
 418   g_return_if_fail (surface != NULL);
 419   g_return_if_fail (f != NULL);
 420   g_return_if_fail (g.nx > 1);
 421   g_return_if_fail (g.ny > 1);
 422   g_return_if_fail (g.nz > 1);
 423
 424   slice1 = gts_iso_slice_new (g.nx, g.ny);
 425   slice2 = gts_iso_slice_new (g.nx, g.ny);
 426   f1 = (gdouble **) malloc2D (g.nx, g.ny, sizeof (gdouble));
 427   f2 = (gdouble **) malloc2D (g.nx, g.ny, sizeof (gdouble));
 428
 429   (*f) (f1, g, 0, data);
 430   g.z += g.dz;
 431   (*f) (f2, g, 1, data);
 432   g.z -= g.dz;
 433   gts_iso_slice_fill_cartesian (slice1, g, f1, f2, iso,
 434                                 surface->vertex_class);
 435   g.z += g.dz;
 436   for (i = 2; i < g.nz; i++) {
 437     g.z += g.dz;
 438     (*f) (f1, g, i, data);
 439     SWAP (f1, f2, tmp);
 440     g.z -= g.dz;
 441     gts_iso_slice_fill_cartesian (slice2, g, f1, f2, iso,
 442                                   surface->vertex_class);
 443     g.z += g.dz;
 444     gts_isosurface_slice (slice1, slice2, surface);
 445     SWAP (slice1, slice2, tmp);
 446   }
 447   gts_iso_slice_fill_cartesian (slice2, g, f2, NULL, iso,
 448                                 surface->vertex_class);
 449   gts_isosurface_slice (slice1, slice2, surface);
 450
 451   gts_iso_slice_destroy (slice1);
 452   gts_iso_slice_destroy (slice2);
 453   free2D ((void **) f1, g.nx);
 454   free2D ((void **) f2, g.nx);
 455 }