iso.c

Go to the documentation of this file.

/* GTS - Library for the manipulation of triangulated surfaces
 * Copyright (C) 1999 Stéphane Popinet
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Library General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Library General Public License for more details.
 *
 * You should have received a copy of the GNU Library General Public
 * License along with this library; if not, write to the
 * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
 * Boston, MA 02111-1307, USA.
 */

#include "gts.h"

typedef enum { LEFT = 0, RIGHT = 1 } Orientation;

typedef struct {
  GtsVertex * v;
  Orientation orientation;
} OrientedVertex;

struct _GtsIsoSlice {
  OrientedVertex *** vertices;
  guint nx, ny;
};

/* coordinates of the edges of the cube (see doc/isocube.fig) */
static guint c[12][4] = {
  {0, 0, 0, 0}, {0, 0, 0, 1}, {0, 0, 1, 1}, {0, 0, 1, 0},
  {1, 0, 0, 0}, {1, 0, 0, 1}, {1, 1, 0, 1}, {1, 1, 0, 0},
  {2, 0, 0, 0}, {2, 1, 0, 0}, {2, 1, 1, 0}, {2, 0, 1, 0}};

/* first index is the edge number, second index is the edge orientation 
   (RIGHT or LEFT), third index are the edges which this edge may connect to
   in order */
static guint edge[12][2][3] = {
  {{9, 1, 8}, {4, 3, 7}},   /* 0 */
  {{6, 2, 5}, {8, 0, 9}},   /* 1 */
  {{10, 3, 11}, {5, 1, 6}}, /* 2 */
  {{7, 0, 4}, {11, 2, 10}}, /* 3 */
  {{3, 7, 0}, {8, 5, 11}},  /* 4 */
  {{11, 4, 8}, {1, 6, 2}},  /* 5 */
  {{2, 5, 1}, {9, 7, 10}},  /* 6 */
  {{10, 6, 9}, {0, 4, 3}},  /* 7 */
  {{5, 11, 4}, {0, 9, 1}},  /* 8 */
  {{1, 8, 0}, {7, 10, 6}},  /* 9 */
  {{6, 9, 7}, {3, 11, 2}},  /* 10 */
  {{2, 10, 3}, {4, 8, 5}}   /* 11 */
};

static void ** malloc2D (guint nx, guint ny, gulong size)
{
  void ** m = g_malloc (nx*sizeof (void *));
  guint i;

  for (i = 0; i < nx; i++)
    m[i] = g_malloc0 (ny*size);

  return m;
}

static void free2D (void ** m, guint nx)
{
  guint i;

  g_return_if_fail (m != NULL);

  for (i = 0; i < nx; i++)
    g_free (m[i]);
  g_free (m);
}

GtsGridPlane * gts_grid_plane_new (guint nx, guint ny)
{
  GtsGridPlane * g = g_malloc (sizeof (GtsGridPlane));

  g->p = (GtsPoint **) malloc2D (nx, ny, sizeof (GtsPoint));
  g->nx = nx;
  g->ny = ny;
  
  return g;
}

void gts_grid_plane_destroy (GtsGridPlane * g)
{
  g_return_if_fail (g != NULL);

  free2D ((void **) g->p, g->nx);
  g_free (g);
}

GtsIsoSlice * gts_iso_slice_new (guint nx, guint ny)
{
  GtsIsoSlice * slice;

  g_return_val_if_fail (nx > 1, NULL);
  g_return_val_if_fail (ny > 1, NULL);

  slice = g_malloc (sizeof (GtsIsoSlice));

  slice->vertices = g_malloc (3*sizeof (OrientedVertex **));
  slice->vertices[0] = 
    (OrientedVertex **) malloc2D (nx, ny, sizeof (OrientedVertex));
  slice->vertices[1] = 
    (OrientedVertex **) malloc2D (nx - 1, ny, sizeof (OrientedVertex));
  slice->vertices[2] = 
    (OrientedVertex **) malloc2D (nx, ny - 1, sizeof (OrientedVertex));
  slice->nx = nx;
  slice->ny = ny;

  return slice;
}

void gts_iso_slice_fill (GtsIsoSlice * slice,
                         GtsGridPlane * plane1,
                         GtsGridPlane * plane2,
                         gdouble ** f1,
                         gdouble ** f2,
                         gdouble iso,
                         GtsVertexClass * klass)
{
  OrientedVertex *** vertices;
  GtsPoint ** p1, ** p2 = NULL;
  guint i, j, nx, ny;

  g_return_if_fail (slice != NULL);
  g_return_if_fail (plane1 != NULL);
  g_return_if_fail (f1 != NULL);
  g_return_if_fail (f2 == NULL || plane2 != NULL);

  p1 = plane1->p;
  if (plane2) 
    p2 = plane2->p;
  vertices = slice->vertices;
  nx = slice->nx;
  ny = slice->ny;

  if (f2)
    for (i = 0; i < nx; i++)
      for (j = 0; j < ny; j++) {
        gdouble v1 = f1[i][j] - iso;
        gdouble v2 = f2[i][j] - iso;
        if ((v1 >= 0. && v2 < 0.) || (v1 < 0. && v2 >= 0.)) {
          gdouble c2 = v1/(v1 - v2), c1 = 1. - c2;
          vertices[0][i][j].v = 
            gts_vertex_new (klass,
                            c1*p1[i][j].x + c2*p2[i][j].x,
                            c1*p1[i][j].y + c2*p2[i][j].y,
                            c1*p1[i][j].z + c2*p2[i][j].z);
          vertices[0][i][j].orientation = v2 >= 0. ? RIGHT : LEFT;
        }
        else
          vertices[0][i][j].v = NULL;
      }
  for (i = 0; i < nx - 1; i++)
    for (j = 0; j < ny; j++) {
      gdouble v1 = f1[i][j] - iso;
      gdouble v2 = f1[i+1][j] - iso;
      if ((v1 >= 0. && v2 < 0.) || (v1 < 0. && v2 >= 0.)) {
        gdouble c2 = v1/(v1 - v2), c1 = 1. - c2;
        vertices[1][i][j].v = 
          gts_vertex_new (klass,
                          c1*p1[i][j].x + c2*p1[i+1][j].x,
                          c1*p1[i][j].y + c2*p1[i+1][j].y,
                          c1*p1[i][j].z + c2*p1[i+1][j].z);
        vertices[1][i][j].orientation = v2 >= 0. ? RIGHT : LEFT;
      }
      else
        vertices[1][i][j].v = NULL;
    }
  for (i = 0; i < nx; i++)
    for (j = 0; j < ny - 1; j++) {
      gdouble v1 = f1[i][j] - iso;
      gdouble v2 = f1[i][j+1] - iso;
      if ((v1 >= 0. && v2 < 0.) || (v1 < 0. && v2 >= 0.)) {
        gdouble c2 = v1/(v1 - v2), c1 = 1. - c2;
        vertices[2][i][j].v = 
          gts_vertex_new (klass,
                          c1*p1[i][j].x + c2*p1[i][j+1].x,
                          c1*p1[i][j].y + c2*p1[i][j+1].y,
                          c1*p1[i][j].z + c2*p1[i][j+1].z);
        vertices[2][i][j].orientation = v2 >= 0. ? RIGHT : LEFT;
      }
      else
        vertices[2][i][j].v = NULL;
    }
}
 
void gts_iso_slice_fill_cartesian (GtsIsoSlice * slice,
                                   GtsCartesianGrid g,
                                   gdouble ** f1,
                                   gdouble ** f2,
                                   gdouble iso,
                                   GtsVertexClass * klass)
{
  OrientedVertex *** vertices;
  guint i, j;
  gdouble x, y;

  g_return_if_fail (slice != NULL);
  g_return_if_fail (f1 != NULL);

  vertices = slice->vertices;

  if (f2)
    for (i = 0, x = g.x; i < g.nx; i++, x += g.dx)
      for (j = 0, y = g.y; j < g.ny; j++, y += g.dy) {
        gdouble v1 = f1[i][j] - iso;
        gdouble v2 = f2[i][j] - iso;
        if ((v1 >= 0. && v2 < 0.) || (v1 < 0. && v2 >= 0.)) {
          vertices[0][i][j].v = 
            gts_vertex_new (klass,
                            x, y, g.z + g.dz*v1/(v1 - v2));
          vertices[0][i][j].orientation = v2 >= 0. ? RIGHT : LEFT;
        }
        else
          vertices[0][i][j].v = NULL;
      }
  for (i = 0, x = g.x; i < g.nx - 1; i++, x += g.dx)
    for (j = 0, y = g.y; j < g.ny; j++, y += g.dy) {
      gdouble v1 = f1[i][j] - iso;
      gdouble v2 = f1[i+1][j] - iso;
      if ((v1 >= 0. && v2 < 0.) || (v1 < 0. && v2 >= 0.)) {
        vertices[1][i][j].v = 
          gts_vertex_new (klass, x + g.dx*v1/(v1 - v2), y, g.z);
        vertices[1][i][j].orientation = v2 >= 0. ? RIGHT : LEFT;
      }
      else
        vertices[1][i][j].v = NULL;
    }
  for (i = 0, x = g.x; i < g.nx; i++, x += g.dx)
    for (j = 0, y = g.y; j < g.ny - 1; j++, y += g.dy) {
      gdouble v1 = f1[i][j] - iso;
      gdouble v2 = f1[i][j+1] - iso;
      if ((v1 >= 0. && v2 < 0.) || (v1 < 0. && v2 >= 0.)) {
        vertices[2][i][j].v = 
          gts_vertex_new (klass, x, y + g.dy*v1/(v1 - v2), g.z);
        vertices[2][i][j].orientation = v2 >= 0. ? RIGHT : LEFT;
      }
      else
        vertices[2][i][j].v = NULL;
    }
}

void gts_iso_slice_destroy (GtsIsoSlice * slice)
{
  g_return_if_fail (slice != NULL);

  free2D ((void **) slice->vertices[0], slice->nx);
  free2D ((void **) slice->vertices[1], slice->nx - 1);
  free2D ((void **) slice->vertices[2], slice->nx);  
  g_free (slice->vertices);
  g_free (slice);
}

void gts_isosurface_slice (GtsIsoSlice * slice1,
                           GtsIsoSlice * slice2,
                           GtsSurface * surface)
{
  guint j, k, l, nx, ny;
  OrientedVertex *** vertices[2];
  GtsVertex * va[12];

  g_return_if_fail (slice1 != NULL);
  g_return_if_fail (slice2 != NULL);
  g_return_if_fail (surface != NULL);
  g_return_if_fail (slice1->nx == slice2->nx && slice1->ny == slice2->ny);

  vertices[0] = slice1->vertices;
  vertices[1] = slice2->vertices;
  nx = slice1->nx;
  ny = slice1->ny;

  /* link vertices with segments and triangles */
  for (j = 0; j < nx - 1; j++)
    for (k = 0; k < ny - 1; k++) {
      gboolean cube_is_cut = FALSE;
      for (l = 0; l < 12; l++) {
        guint nv = 0, e = l;
        OrientedVertex ov = 
          vertices[c[e][1]][c[e][0]][j + c[e][2]][k + c[e][3]];
        while (ov.v && !GTS_OBJECT (ov.v)->reserved) {
          guint m = 0, * ne = edge[e][ov.orientation];
          va[nv++] = ov.v;
          GTS_OBJECT (ov.v)->reserved = surface;
          ov.v = NULL;
          while (m < 3 && !ov.v) {
            e = ne[m++];
            ov = vertices[c[e][1]][c[e][0]][j + c[e][2]][k + c[e][3]];
          }
        }
        /* create edges and faces */
        if (nv > 2) {
          GtsEdge * e1, * e2, * e3;
          guint m;
          if (!(e1 = GTS_EDGE (gts_vertices_are_connected (va[0], va[1]))))
            e1 = gts_edge_new (surface->edge_class, va[0], va[1]);
          for (m = 1; m < nv - 1; m++) {
            if (!(e2 = GTS_EDGE (gts_vertices_are_connected (va[m], va[m+1]))))
              e2 = gts_edge_new (surface->edge_class, va[m], va[m+1]);
            if (!(e3 = GTS_EDGE (gts_vertices_are_connected (va[m+1], va[0]))))
              e3 = gts_edge_new (surface->edge_class, va[m+1], va[0]);
            gts_surface_add_face (surface, 
                                  gts_face_new (surface->face_class,
                                                e1, e2, e3));
            e1 = e3;
          }
        }
        if (nv > 0)
          cube_is_cut = TRUE;
      }
      if (cube_is_cut)
        for (l = 0; l < 12; l++) {
          GtsVertex * v = 
            vertices[c[l][1]][c[l][0]][j + c[l][2]][k + c[l][3]].v;
          if (v)
            GTS_OBJECT (v)->reserved = NULL;
        }
    }
}

#define SWAP(s1, s2, tmp) (tmp = s1, s1 = s2, s2 = tmp)

void gts_isosurface_cartesian (GtsSurface * surface,
                               GtsCartesianGrid g,
                               GtsIsoCartesianFunc f,
                               gpointer data,
                               gdouble iso)
{
  void * tmp;
  gdouble ** f1, ** f2;
  GtsIsoSlice * slice1, * slice2;
  guint i;

  g_return_if_fail (surface != NULL);
  g_return_if_fail (f != NULL);
  g_return_if_fail (g.nx > 1);
  g_return_if_fail (g.ny > 1);
  g_return_if_fail (g.nz > 1);

  slice1 = gts_iso_slice_new (g.nx, g.ny);
  slice2 = gts_iso_slice_new (g.nx, g.ny);
  f1 = (gdouble **) malloc2D (g.nx, g.ny, sizeof (gdouble));
  f2 = (gdouble **) malloc2D (g.nx, g.ny, sizeof (gdouble));

  (*f) (f1, g, 0, data);
  g.z += g.dz;
  (*f) (f2, g, 1, data);
  g.z -= g.dz;
  gts_iso_slice_fill_cartesian (slice1, g, f1, f2, iso, 
                                surface->vertex_class);
  g.z += g.dz;
  for (i = 2; i < g.nz; i++) {
    g.z += g.dz;
    (*f) (f1, g, i, data);
    SWAP (f1, f2, tmp);
    g.z -= g.dz;
    gts_iso_slice_fill_cartesian (slice2, g, f1, f2, iso, 
                                  surface->vertex_class);
    g.z += g.dz;
    gts_isosurface_slice (slice1, slice2, surface);
    SWAP (slice1, slice2, tmp);
  }
  gts_iso_slice_fill_cartesian (slice2, g, f2, NULL, iso,
                                surface->vertex_class);
  gts_isosurface_slice (slice1, slice2, surface);

  gts_iso_slice_destroy (slice1);
  gts_iso_slice_destroy (slice2);
  free2D ((void **) f1, g.nx);
  free2D ((void **) f2, g.nx);
}