From aae93be6e35c5c760d87f08e7993194adf1f6e44 Mon Sep 17 00:00:00 2001 From: Andy Spencer Date: Mon, 24 Jan 2011 03:30:34 +0000 Subject: [PATCH] Add GritsVolume for 3D volume rendering This is based on the Marching Cubes algorithm. Be warned, the API for this is not quite entirely stable.. --- src/grits.h | 1 + src/objects/Makefile.am | 8 +- src/objects/grits-volume.c | 175 ++++++++++++ src/objects/grits-volume.h | 76 +++++ src/objects/marching.c | 556 +++++++++++++++++++++++++++++++++++++ src/objects/marching.h | 82 ++++++ 6 files changed, 896 insertions(+), 2 deletions(-) create mode 100644 src/objects/grits-volume.c create mode 100644 src/objects/grits-volume.h create mode 100644 src/objects/marching.c create mode 100644 src/objects/marching.h diff --git a/src/grits.h b/src/grits.h index 94cc1ab..6885042 100644 --- a/src/grits.h +++ b/src/grits.h @@ -34,6 +34,7 @@ #include #include #include +#include /* Plugins */ #include diff --git a/src/objects/Makefile.am b/src/objects/Makefile.am index 19a3b84..bb1a2a4 100644 --- a/src/objects/Makefile.am +++ b/src/objects/Makefile.am @@ -9,14 +9,18 @@ grits_objects_include_HEADERS = \ grits-object.h \ grits-marker.h \ grits-callback.h \ - grits-tile.h + grits-tile.h \ + grits-volume.h \ + marching.h noinst_LTLIBRARIES = libgrits-objects.la libgrits_objects_la_SOURCES = \ grits-object.c grits-object.h \ grits-marker.c grits-marker.h \ grits-callback.c grits-callback.h \ - grits-tile.c grits-tile.h + grits-tile.c grits-tile.h \ + grits-volume.c grits-volume.h \ + marching.c marching.h libgrits_objects_la_LDFLAGS = -static MAINTAINERCLEANFILES = Makefile.in diff --git a/src/objects/grits-volume.c b/src/objects/grits-volume.c new file mode 100644 index 0000000..be08279 --- /dev/null +++ b/src/objects/grits-volume.c @@ -0,0 +1,175 @@ +/* + * Copyright (C) 2009-2010 Andy Spencer + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program 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 General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +/** + * SECTION:grits-volume + * @short_description: 3-D gridded vovlume + * + * Each #GritsVolume consistes of a 3-dimentional grid of data points each + * consisting of a value. + * + * Currently iso-surfaces are extracted and displayed when rendering the volume + * data. + */ + +#include +#include +#include +#include +#include +#include "grits-volume.h" + +/* Drawing */ +static void draw_points(VolGrid *grid, gdouble level) +{ + glPointSize(200./grid->xs); + glBegin(GL_POINTS); + for (int x = 0; x < grid->xs; x++) + for (int y = 0; y < grid->ys; y++) + for (int z = 0; z < grid->zs; z++) { + VolPoint *pt = vol_grid_get(grid, x, y, z); + if (pt->value < level) + continue; + + g_debug("(%d,%d,%d) value=%f", x, y, z, pt->value); + + glColor4f(1.0, 1.0, 1.0, pt->value/100); + glVertex3dv((double*)&pt->c); + } + glEnd(); +} + +static void draw_iso(GList *tris) +{ + g_debug("GritsVolume: draw_iso"); + glDisable(GL_CULL_FACE); + glBegin(GL_TRIANGLES); + for (GList *cur = tris; cur; cur = cur->next) { + VolTriangle *tri = cur->data; + VolCoord c[3] = {tri->v[0]->c, tri->v[1]->c, tri->v[2]->c }; + VolCoord n[3] = {tri->v[0]->norm, tri->v[1]->norm, tri->v[2]->norm}; + + /* Normalize normal vector */ + for (int i = 0; i < 3; i++) { + double total = sqrt( + n[i].x * n[i].x + + n[i].y * n[i].y + + n[i].z * n[i].z); + if (total == 0) + continue; // wtf + n[i].x = n[i].x / total; + n[i].y = n[i].y / total; + n[i].z = n[i].z / total; + } + + //for (int i = 0; i < 3; i++) + // g_debug("norm=%f,%f,%f", + // n[i].x, n[i].y, n[i].z); + //glNormal3dv((double*)&tri->norm); + glNormal3dv((double*)&n[0]); glVertex3dv((double*)&c[0]); + glNormal3dv((double*)&n[1]); glVertex3dv((double*)&c[1]); + glNormal3dv((double*)&n[2]); glVertex3dv((double*)&c[2]); + } + glEnd(); +} + +static void draw(GritsObject *_volume, GritsOpenGL *opengl) +{ + g_debug("GritsVolume: draw"); + GritsVolume *volume = GRITS_VOLUME(_volume); + + gfloat amb[4] = {}; + gfloat dif[4] = {}; + switch (volume->disp) { + case GRITS_VOLUME_SURFACE: + glDisable(GL_COLOR_MATERIAL); + amb[0] = (double)volume->color[0] / 0xff; + amb[1] = (double)volume->color[1] / 0xff; + amb[2] = (double)volume->color[2] / 0xff; + amb[3] = 1; + dif[0] = (double)volume->color[0] / 0xff; + dif[1] = (double)volume->color[1] / 0xff; + dif[2] = (double)volume->color[2] / 0xff; + dif[3] = 1; + glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT, amb); + glMaterialfv(GL_FRONT_AND_BACK, GL_DIFFUSE, dif); + draw_iso(volume->tris); + break; + case GRITS_VOLUME_POINTS: + draw_points(volume->grid, volume->level); + break; + } +} + +/* Idle update */ +gboolean update_iso(gpointer _volume) +{ + GritsVolume *volume = _volume; + if (volume->tris) { + g_list_foreach(volume->tris, (GFunc)vol_triangle_free, NULL); + g_list_free(volume->tris); + } + volume->tris = marching_cubes(volume->grid, volume->level); + volume->update_id = 0; + grits_object_queue_draw(GRITS_OBJECT(volume)); + return FALSE; +} + +/*********** + * Methods * + ***********/ +void grits_volume_set_level(GritsVolume *volume, gdouble level) +{ + volume->level = level; + if (!volume->update_id) + volume->update_id = g_idle_add(update_iso, volume); +} + +GritsVolume *grits_volume_new(VolGrid *grid) +{ + g_debug("GritsVolume: new - %p[%d][%d][%d]", + grid, grid->xs, grid->ys, grid->zs); + GritsVolume *volume = g_object_new(GRITS_TYPE_VOLUME, NULL); + volume->grid = grid; + return volume; +} + +/* GObject code */ +G_DEFINE_TYPE(GritsVolume, grits_volume, GRITS_TYPE_OBJECT); +static void grits_volume_init(GritsVolume *volume) +{ +} + +static void grits_volume_finalize(GObject *_volume) +{ + GritsVolume *volume = GRITS_VOLUME(_volume); + volume->color[0] = 1; + volume->color[1] = 1; + volume->color[2] = 1; + volume->color[3] = 1; + //g_debug("GritsVolume: finalize - %s", volume->label); +} + +static void grits_volume_class_init(GritsVolumeClass *klass) +{ + g_debug("GritsVolume: class_init"); + GObjectClass *gobject_class = G_OBJECT_CLASS(klass); + gobject_class->finalize = grits_volume_finalize; + + GritsObjectClass *object_class = GRITS_OBJECT_CLASS(klass); + object_class->draw = draw; +} diff --git a/src/objects/grits-volume.h b/src/objects/grits-volume.h new file mode 100644 index 0000000..6eb4e83 --- /dev/null +++ b/src/objects/grits-volume.h @@ -0,0 +1,76 @@ +/* + * Copyright (C) 2009-2010 Andy Spencer + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program 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 General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#ifndef __GRITS_VOLUME_H__ +#define __GRITS_VOLUME_H__ + +#include +#include +#include +#include "grits-object.h" +#include "marching.h" + +/*************** + * GritsVolume * + ***************/ +#define GRITS_TYPE_VOLUME (grits_volume_get_type()) +#define GRITS_VOLUME(obj) (G_TYPE_CHECK_INSTANCE_CAST((obj), GRITS_TYPE_VOLUME, GritsVolume)) +#define GRITS_IS_VOLUME(obj) (G_TYPE_CHECK_INSTANCE_TYPE((obj), GRITS_TYPE_VOLUME)) +#define GRITS_VOLUME_CLASS(klass) (G_TYPE_CHECK_CLASS_CAST ((klass), GRITS_TYPE_VOLUME, GritsVolumeClass)) +#define GRITS_IS_VOLUME_CLASS(klass) (G_TYPE_CHECK_CLASS_TYPE ((klass), GRITS_TYPE_VOLUME)) +#define GRITS_VOLUME_GET_CLASS(obj) (G_TYPE_INSTANCE_GET_CLASS ((obj), GRITS_TYPE_VOLUME, GritsVolumeClass)) + +typedef enum { + GRITS_VOLUME_CARTESIAN, + GRITS_VOLUME_SPHERICAL, +} GritsVolumeProj; + +typedef enum { + GRITS_VOLUME_SURFACE, + GRITS_VOLUME_POINTS, +} GritsVolumeDisp; + +typedef struct _GritsVolume GritsVolume; +typedef struct _GritsVolumeClass GritsVolumeClass; + +struct _GritsVolume { + GritsObject parent_instance; + + /* Instance members */ + GritsVolumeProj proj; // projection + GritsVolumeDisp disp; // display mode + + /* Internal */ + VolGrid *grid; + GList *tris; + gdouble level; + guint8 color[4]; + gint update_id; +}; + +struct _GritsVolumeClass { + GritsObjectClass parent_class; +}; + +GType grits_volume_get_type(void); + +/* Methods */ +void grits_volume_set_level(GritsVolume *volume, gdouble level); + +GritsVolume *grits_volume_new(VolGrid *grid); + +#endif diff --git a/src/objects/marching.c b/src/objects/marching.c new file mode 100644 index 0000000..30a4e3b --- /dev/null +++ b/src/objects/marching.c @@ -0,0 +1,556 @@ +/** + * Copyright (C) 2002 Jamie Zawinski + * Copyright (C) 2009-2010 Andy Spencer + * + * Utility functions to create "marching cubes" meshes from 3d fields. + * + * Permission to use, copy, modify, distribute, and sell this software and its + * documentation for any purpose is hereby granted without fee, provided that + * the above copyright notice appear in all copies and that both that + * copyright notice and this permission notice appear in supporting + * documentation. No representations are made about the suitability of this + * software for any purpose. It is provided "as is" without express or + * implied warranty. + * + * Marching cubes implementation by Paul Bourke + * http://astronomy.swin.edu.au/~pbourke/modelling/polygonise/ + */ + +#include +#include +#include + +#include "grits-util.h" +#include "marching.h" + +/** + * Indexing convention: + * + * Vertices: Edges: + * + * 4 ______________ 5 ______________ + * /| /| /| 4 /| + * / | 6 / | 7 / |8 5 / | + * 7 /_____________/ | /______________/ | 9 + * | | | | | | 6 | | + * | 0 |_________|___| 1 | |_________|10_| + * | / | / 11 | 3/ 0 | / + * | / | / | / | / 1 + * 3 |/____________|/ 2 |/____________|/ + * 2 + */ + +static const int edge_table[256] = { + 0x0 , 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c, + 0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00, + 0x190, 0x99 , 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c, + 0x99c, 0x895, 0xb9f, 0xa96, 0xd9a, 0xc93, 0xf99, 0xe90, + 0x230, 0x339, 0x33 , 0x13a, 0x636, 0x73f, 0x435, 0x53c, + 0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30, + 0x3a0, 0x2a9, 0x1a3, 0xaa , 0x7a6, 0x6af, 0x5a5, 0x4ac, + 0xbac, 0xaa5, 0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0, + 0x460, 0x569, 0x663, 0x76a, 0x66 , 0x16f, 0x265, 0x36c, + 0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963, 0xa69, 0xb60, + 0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0xff , 0x3f5, 0x2fc, + 0xdfc, 0xcf5, 0xfff, 0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0, + 0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x55 , 0x15c, + 0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950, + 0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6, 0x2cf, 0x1c5, 0xcc , + 0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0, + 0x8c0, 0x9c9, 0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc, + 0xcc , 0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9, 0x7c0, + 0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c, + 0x15c, 0x55 , 0x35f, 0x256, 0x55a, 0x453, 0x759, 0x650, + 0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc, + 0x2fc, 0x3f5, 0xff , 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0, + 0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f, 0xd65, 0xc6c, + 0x36c, 0x265, 0x16f, 0x66 , 0x76a, 0x663, 0x569, 0x460, + 0xca0, 0xda9, 0xea3, 0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac, + 0x4ac, 0x5a5, 0x6af, 0x7a6, 0xaa , 0x1a3, 0x2a9, 0x3a0, + 0xd30, 0xc39, 0xf33, 0xe3a, 0x936, 0x83f, 0xb35, 0xa3c, + 0x53c, 0x435, 0x73f, 0x636, 0x13a, 0x33 , 0x339, 0x230, + 0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c, + 0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x99 , 0x190, + 0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c, + 0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x0 +}; + +static const int tri_table[256][16] = { + {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 0, 8, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 9, 2, 10, 0, 2, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 2, 8, 3, 2, 10, 8, 10, 9, 8, -1, -1, -1, -1, -1, -1, -1}, + { 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 0, 11, 2, 8, 11, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 1, 11, 2, 1, 9, 11, 9, 8, 11, -1, -1, -1, -1, -1, -1, -1}, + { 3, 10, 1, 11, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 0, 10, 1, 0, 8, 10, 8, 11, 10, -1, -1, -1, -1, -1, -1, -1}, + { 3, 9, 0, 3, 11, 9, 11, 10, 9, -1, -1, -1, -1, -1, -1, -1}, + { 9, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 4, 3, 0, 7, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 0, 1, 9, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 4, 1, 9, 4, 7, 1, 7, 3, 1, -1, -1, -1, -1, -1, -1, -1}, + { 1, 2, 10, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 3, 4, 7, 3, 0, 4, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1}, + { 9, 2, 10, 9, 0, 2, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1}, + { 2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1, -1, -1, -1}, + { 8, 4, 7, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {11, 4, 7, 11, 2, 4, 2, 0, 4, -1, -1, -1, -1, -1, -1, -1}, + { 9, 0, 1, 8, 4, 7, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1}, + { 4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, -1, -1, -1, -1}, + { 3, 10, 1, 3, 11, 10, 7, 8, 4, -1, -1, -1, -1, -1, -1, -1}, + { 1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, -1, -1, -1, -1}, + { 4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, -1, -1, -1, -1}, + { 4, 7, 11, 4, 11, 9, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1}, + { 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 9, 5, 4, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 0, 5, 4, 1, 5, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 8, 5, 4, 8, 3, 5, 3, 1, 5, -1, -1, -1, -1, -1, -1, -1}, + { 1, 2, 10, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 3, 0, 8, 1, 2, 10, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1}, + { 5, 2, 10, 5, 4, 2, 4, 0, 2, -1, -1, -1, -1, -1, -1, -1}, + { 2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, -1, -1, -1, -1}, + { 9, 5, 4, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 0, 11, 2, 0, 8, 11, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1}, + { 0, 5, 4, 0, 1, 5, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1}, + { 2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, -1, -1, -1, -1}, + {10, 3, 11, 10, 1, 3, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1}, + { 4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, -1, -1, -1, -1}, + { 5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, -1, -1, -1, -1}, + { 5, 4, 8, 5, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1}, + { 9, 7, 8, 5, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 9, 3, 0, 9, 5, 3, 5, 7, 3, -1, -1, -1, -1, -1, -1, -1}, + { 0, 7, 8, 0, 1, 7, 1, 5, 7, -1, -1, -1, -1, -1, -1, -1}, + { 1, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 9, 7, 8, 9, 5, 7, 10, 1, 2, -1, -1, -1, -1, -1, -1, -1}, + {10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, -1, -1, -1, -1}, + { 8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, -1, -1, -1, -1}, + { 2, 10, 5, 2, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1}, + { 7, 9, 5, 7, 8, 9, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1}, + { 9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, -1, -1, -1, -1}, + { 2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, -1, -1, -1, -1}, + {11, 2, 1, 11, 1, 7, 7, 1, 5, -1, -1, -1, -1, -1, -1, -1}, + { 9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, -1, -1, -1, -1}, + { 5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0, -1}, + {11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0, -1}, + {11, 10, 5, 7, 11, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 0, 8, 3, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 9, 0, 1, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 1, 8, 3, 1, 9, 8, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1}, + { 1, 6, 5, 2, 6, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 1, 6, 5, 1, 2, 6, 3, 0, 8, -1, -1, -1, -1, -1, -1, -1}, + { 9, 6, 5, 9, 0, 6, 0, 2, 6, -1, -1, -1, -1, -1, -1, -1}, + { 5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, -1, -1, -1, -1}, + { 2, 3, 11, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {11, 0, 8, 11, 2, 0, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1}, + { 0, 1, 9, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1}, + { 5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, -1, -1, -1, -1}, + { 6, 3, 11, 6, 5, 3, 5, 1, 3, -1, -1, -1, -1, -1, -1, -1}, + { 0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, -1, -1, -1, -1}, + { 3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, -1, -1, -1, -1}, + { 6, 5, 9, 6, 9, 11, 11, 9, 8, -1, -1, -1, -1, -1, -1, -1}, + { 5, 10, 6, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 4, 3, 0, 4, 7, 3, 6, 5, 10, -1, -1, -1, -1, -1, -1, -1}, + { 1, 9, 0, 5, 10, 6, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1}, + {10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, -1, -1, -1, -1}, + { 6, 1, 2, 6, 5, 1, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1}, + { 1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, -1, -1, -1, -1}, + { 8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, -1, -1, -1, -1}, + { 7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9, -1}, + { 3, 11, 2, 7, 8, 4, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1}, + { 5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, -1, -1, -1, -1}, + { 0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1}, + { 9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6, -1}, + { 8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, -1, -1, -1, -1}, + { 5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11, -1}, + { 0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7, -1}, + { 6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9, -1, -1, -1, -1}, + {10, 4, 9, 6, 4, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 4, 10, 6, 4, 9, 10, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1}, + {10, 0, 1, 10, 6, 0, 6, 4, 0, -1, -1, -1, -1, -1, -1, -1}, + { 8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, -1, -1, -1, -1}, + { 1, 4, 9, 1, 2, 4, 2, 6, 4, -1, -1, -1, -1, -1, -1, -1}, + { 3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, -1, -1, -1, -1}, + { 0, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 8, 3, 2, 8, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1}, + {10, 4, 9, 10, 6, 4, 11, 2, 3, -1, -1, -1, -1, -1, -1, -1}, + { 0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, -1, -1, -1, -1}, + { 3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, -1, -1, -1, -1}, + { 6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1, -1}, + { 9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, -1, -1, -1, -1}, + { 8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1, -1}, + { 3, 11, 6, 3, 6, 0, 0, 6, 4, -1, -1, -1, -1, -1, -1, -1}, + { 6, 4, 8, 11, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 7, 10, 6, 7, 8, 10, 8, 9, 10, -1, -1, -1, -1, -1, -1, -1}, + { 0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, -1, -1, -1, -1}, + {10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, -1, -1, -1, -1}, + {10, 6, 7, 10, 7, 1, 1, 7, 3, -1, -1, -1, -1, -1, -1, -1}, + { 1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, -1, -1, -1, -1}, + { 2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9, -1}, + { 7, 8, 0, 7, 0, 6, 6, 0, 2, -1, -1, -1, -1, -1, -1, -1}, + { 7, 3, 2, 6, 7, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, -1, -1, -1, -1}, + { 2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7, -1}, + { 1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11, -1}, + {11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1, -1, -1, -1, -1}, + { 8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6, -1}, + { 0, 9, 1, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, -1, -1, -1, -1}, + { 7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 3, 0, 8, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 0, 1, 9, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 8, 1, 9, 8, 3, 1, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1}, + {10, 1, 2, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 1, 2, 10, 3, 0, 8, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1}, + { 2, 9, 0, 2, 10, 9, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1}, + { 6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8, -1, -1, -1, -1}, + { 7, 2, 3, 6, 2, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 7, 0, 8, 7, 6, 0, 6, 2, 0, -1, -1, -1, -1, -1, -1, -1}, + { 2, 7, 6, 2, 3, 7, 0, 1, 9, -1, -1, -1, -1, -1, -1, -1}, + { 1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, -1, -1, -1, -1}, + {10, 7, 6, 10, 1, 7, 1, 3, 7, -1, -1, -1, -1, -1, -1, -1}, + {10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, -1, -1, -1, -1}, + { 0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, -1, -1, -1, -1}, + { 7, 6, 10, 7, 10, 8, 8, 10, 9, -1, -1, -1, -1, -1, -1, -1}, + { 6, 8, 4, 11, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 3, 6, 11, 3, 0, 6, 0, 4, 6, -1, -1, -1, -1, -1, -1, -1}, + { 8, 6, 11, 8, 4, 6, 9, 0, 1, -1, -1, -1, -1, -1, -1, -1}, + { 9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, -1, -1, -1, -1}, + { 6, 8, 4, 6, 11, 8, 2, 10, 1, -1, -1, -1, -1, -1, -1, -1}, + { 1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6, -1, -1, -1, -1}, + { 4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, -1, -1, -1, -1}, + {10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3, -1}, + { 8, 2, 3, 8, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1}, + { 0, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8, -1, -1, -1, -1}, + { 1, 9, 4, 1, 4, 2, 2, 4, 6, -1, -1, -1, -1, -1, -1, -1}, + { 8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, -1, -1, -1, -1}, + {10, 1, 0, 10, 0, 6, 6, 0, 4, -1, -1, -1, -1, -1, -1, -1}, + { 4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3, -1}, + {10, 9, 4, 6, 10, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 4, 9, 5, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 0, 8, 3, 4, 9, 5, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1}, + { 5, 0, 1, 5, 4, 0, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1}, + {11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5, -1, -1, -1, -1}, + { 9, 5, 4, 10, 1, 2, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1}, + { 6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5, -1, -1, -1, -1}, + { 7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2, -1, -1, -1, -1}, + { 3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6, -1}, + { 7, 2, 3, 7, 6, 2, 5, 4, 9, -1, -1, -1, -1, -1, -1, -1}, + { 9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7, -1, -1, -1, -1}, + { 3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0, -1, -1, -1, -1}, + { 6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8, -1}, + { 9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7, -1, -1, -1, -1}, + { 1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4, -1}, + { 4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10, -1}, + { 7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10, -1, -1, -1, -1}, + { 6, 9, 5, 6, 11, 9, 11, 8, 9, -1, -1, -1, -1, -1, -1, -1}, + { 3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, -1, -1, -1, -1}, + { 0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, -1, -1, -1, -1}, + { 6, 11, 3, 6, 3, 5, 5, 3, 1, -1, -1, -1, -1, -1, -1, -1}, + { 1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, -1, -1, -1, -1}, + { 0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10, -1}, + {11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5, -1}, + { 6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3, -1, -1, -1, -1}, + { 5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, -1, -1, -1, -1}, + { 9, 5, 6, 9, 6, 0, 0, 6, 2, -1, -1, -1, -1, -1, -1, -1}, + { 1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8, -1}, + { 1, 5, 6, 2, 1, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6, -1}, + {10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0, -1, -1, -1, -1}, + { 0, 3, 8, 5, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {10, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {11, 5, 10, 7, 5, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {11, 5, 10, 11, 7, 5, 8, 3, 0, -1, -1, -1, -1, -1, -1, -1}, + { 5, 11, 7, 5, 10, 11, 1, 9, 0, -1, -1, -1, -1, -1, -1, -1}, + {10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1, -1, -1, -1, -1}, + {11, 1, 2, 11, 7, 1, 7, 5, 1, -1, -1, -1, -1, -1, -1, -1}, + { 0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11, -1, -1, -1, -1}, + { 9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, -1, -1, -1, -1}, + { 7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2, -1}, + { 2, 5, 10, 2, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1}, + { 8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, -1, -1, -1, -1}, + { 9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2, -1, -1, -1, -1}, + { 9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2, -1}, + { 1, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 0, 8, 7, 0, 7, 1, 1, 7, 5, -1, -1, -1, -1, -1, -1, -1}, + { 9, 0, 3, 9, 3, 5, 5, 3, 7, -1, -1, -1, -1, -1, -1, -1}, + { 9, 8, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 5, 8, 4, 5, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1}, + { 5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, -1, -1, -1, -1}, + { 0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, -1, -1, -1, -1}, + {10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4, -1}, + { 2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, -1, -1, -1, -1}, + { 0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11, -1}, + { 0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5, -1}, + { 9, 4, 5, 2, 11, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, -1, -1, -1, -1}, + { 5, 10, 2, 5, 2, 4, 4, 2, 0, -1, -1, -1, -1, -1, -1, -1}, + { 3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9, -1}, + { 5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2, -1, -1, -1, -1}, + { 8, 4, 5, 8, 5, 3, 3, 5, 1, -1, -1, -1, -1, -1, -1, -1}, + { 0, 4, 5, 1, 0, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5, -1, -1, -1, -1}, + { 9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 4, 11, 7, 4, 9, 11, 9, 10, 11, -1, -1, -1, -1, -1, -1, -1}, + { 0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11, -1, -1, -1, -1}, + { 1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, -1, -1, -1, -1}, + { 3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4, -1}, + { 4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2, -1, -1, -1, -1}, + { 9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3, -1}, + {11, 7, 4, 11, 4, 2, 2, 4, 0, -1, -1, -1, -1, -1, -1, -1}, + {11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4, -1, -1, -1, -1}, + { 2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, -1, -1, -1, -1}, + { 9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7, -1}, + { 3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10, -1}, + { 1, 10, 2, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 4, 9, 1, 4, 1, 7, 7, 1, 3, -1, -1, -1, -1, -1, -1, -1}, + { 4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1, -1, -1, -1, -1}, + { 4, 0, 3, 7, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 9, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 3, 0, 9, 3, 9, 11, 11, 9, 10, -1, -1, -1, -1, -1, -1, -1}, + { 0, 1, 10, 0, 10, 8, 8, 10, 11, -1, -1, -1, -1, -1, -1, -1}, + { 3, 1, 10, 11, 3, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 1, 2, 11, 1, 11, 9, 9, 11, 8, -1, -1, -1, -1, -1, -1, -1}, + { 3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9, -1, -1, -1, -1}, + { 0, 2, 11, 8, 0, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 3, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 2, 3, 8, 2, 8, 10, 10, 8, 9, -1, -1, -1, -1, -1, -1, -1}, + { 9, 10, 2, 0, 9, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8, -1, -1, -1, -1}, + { 1, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 1, 3, 8, 9, 1, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 0, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + { 0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1} +}; + +static void do_normals(VolTriangle *tri) +{ + /* Triangle normal */ + crossd3((gdouble*)tri->v[0], + (gdouble*)tri->v[1], + (gdouble*)tri->v[2], + (gdouble*)&tri->norm); + + double mag = sqrt(tri->norm.x*tri->norm.x + + tri->norm.y*tri->norm.y + + tri->norm.z*tri->norm.z); + if (mag == 0) + return; // wtf + + /* Vertex normal */ + for (int vi = 0; vi < 3; vi++) { + double *tri_norm = (double*)&tri->norm; + double *vert_norm = (double*)&tri->v[vi]->norm; + for (int i = 0; i < 3; i++) + vert_norm[i] += tri_norm[i]/mag/mag; + } +} + +/* Linearly interpolate the position where an isosurface cuts + * an edge between two vertices, each with their own scalar value + */ +gdouble dist = 0; +static VolVertex *interp_vertex( + VolPoint *p1, VolVertex *v, VolPoint *p2, double level) +{ + dist = 0; + + /* caching */ + if (v) return v; + + /* Invalid values */ + if (isnan(p1->value)) level = p2->value; + if (isnan(p2->value)) level = p1->value; + + /* Debug */ + dist = distd((gdouble*)p1, (gdouble*)p2); + + double mu = (level - p1->value) / (p2->value - p1->value); + v = g_new0(VolVertex, 1); + //v->c.x = (p2->c.x + p1->c.x)/2; + //v->c.y = (p2->c.y + p1->c.y)/2; + //v->c.z = (p2->c.z + p1->c.z)/2; + v->c.x = p1->c.x + mu * (p2->c.x - p1->c.x); + v->c.y = p1->c.y + mu * (p2->c.y - p1->c.y); + v->c.z = p1->c.z + mu * (p2->c.z - p1->c.z); + return v; +} + + +/** + * By Paul Bourke + * Modified by Andy Spencer + */ +GList *march_one_cube(VolCell *cell, double level) +{ + VolPoint **c = cell->corner; + VolVertex **e = cell->edge; + + /* Determine the index into the edge table which + * tells us which vertices are inside of the surface */ + int ci = 0; + for (int i = 0; i < 8; i++) + if (c[i]->value < level || isnan(c[i]->value)) + ci |= 1 << i; + + /* Cube is entirely in/out of the surface */ + if (edge_table[ci] == 0) + return 0; + + /* Find the vertices where the surface intersects the cube */ + if (edge_table[ci] & 1<<0 ) e[0 ] = interp_vertex(c[0], e[0 ], c[1], level); + if (edge_table[ci] & 1<<1 ) e[1 ] = interp_vertex(c[1], e[1 ], c[2], level); + if (edge_table[ci] & 1<<2 ) e[2 ] = interp_vertex(c[2], e[2 ], c[3], level); + if (edge_table[ci] & 1<<3 ) e[3 ] = interp_vertex(c[3], e[3 ], c[0], level); + if (edge_table[ci] & 1<<4 ) e[4 ] = interp_vertex(c[4], e[4 ], c[5], level); + if (edge_table[ci] & 1<<5 ) e[5 ] = interp_vertex(c[5], e[5 ], c[6], level); + if (edge_table[ci] & 1<<6 ) e[6 ] = interp_vertex(c[6], e[6 ], c[7], level); + if (edge_table[ci] & 1<<7 ) e[7 ] = interp_vertex(c[7], e[7 ], c[4], level); + if (edge_table[ci] & 1<<8 ) e[8 ] = interp_vertex(c[0], e[8 ], c[4], level); + if (edge_table[ci] & 1<<9 ) e[9 ] = interp_vertex(c[1], e[9 ], c[5], level); + if (edge_table[ci] & 1<<10) e[10] = interp_vertex(c[2], e[10], c[6], level); + if (edge_table[ci] & 1<<11) e[11] = interp_vertex(c[3], e[11], c[7], level); + + /* Create the triangle */ + GList *tris = NULL; + for (int i = 0; tri_table[ci][i] != -1; i += 3) { + VolTriangle *tri = g_new0(VolTriangle, 1); + /* Use z,y,z */ + tri->v[2] = e[tri_table[ci][i+0]]; + tri->v[1] = e[tri_table[ci][i+1]]; + tri->v[0] = e[tri_table[ci][i+2]]; + tri->v[2]->tris++; + tri->v[1]->tris++; + tri->v[0]->tris++; + do_normals(tri); + tris = g_list_prepend(tris, tri); + } + return tris; +} + +/* Find triangles along an isosurface in a grid of points */ +//GList *marching_cubes_simple(VolGrid *grid, double level) +//{ +// GList *tris = NULL; +// for (int x = 0; x+1 < grid->xs; x++) +// for (int y = 0; y+1 < grid->ys; y++) +// for (int z = 0; z+1 < grid->zs; z++) { +// VolCell cell = {}; +// cell.corner[0] = &grid->points[x ][y ][z ]; +// cell.corner[1] = &grid->points[x ][y+1][z ]; +// cell.corner[2] = &grid->points[x+1][y+1][z ]; +// cell.corner[3] = &grid->points[x+1][y ][z ]; +// cell.corner[4] = &grid->points[x ][y ][z+1]; +// cell.corner[5] = &grid->points[x ][y+1][z+1]; +// cell.corner[6] = &grid->points[x+1][y+1][z+1]; +// cell.corner[7] = &grid->points[x+1][y ][z+1]; +// GList *these = march_one_cube(&cell, level); +// tris = g_list_concat(these, tris); +// } +// return tris; +//} + +/** + * Indexing convention: + * + * Vertices: Edges: + * + * 4 ______________ 5 ______________ + * /| /| /| 4 /| | 2 + * / | 6 / | 7 / |8 5 / | | + * 7 /_____________/ | /______________/ | 9 | 1 + * | | | | | | 6 | | X---------- + * | 0 |_________|___| 1 | |_________|10_| / + * | / | / 11 | 3/ 0 | / / + * | / | / | / | / 1 / 0 + * 3 |/____________|/ 2 |/____________|/ + * 2 + */ + +GList *marching_cubes(VolGrid *grid, double level) +{ + int xs = grid->xs; + int ys = grid->ys; + int zs = grid->zs; + + GList *tris = NULL; + VolVertex **verts = g_new0(VolVertex*, xs*ys*2*3); +#define EDGE(x, y, z, num) (num + (x)*3 + (y)*xs*3 + ((z)%2)*ys*xs*3) + for (int z = 0; z+1 < zs; z++) { + memset(&verts[EDGE(0,0,z+1,0)], 0, xs*ys*3*sizeof(VolVertex*)); + for (int y = 0; y+1 < ys; y++) + for (int x = 0; x+1 < xs; x++) { + VolCell cell = {}; + cell.corner[0] = vol_grid_get(grid, x , y , z ); + cell.corner[1] = vol_grid_get(grid, x , y+1, z ); + cell.corner[2] = vol_grid_get(grid, x+1, y+1, z ); + cell.corner[3] = vol_grid_get(grid, x+1, y , z ); + cell.corner[4] = vol_grid_get(grid, x , y , z+1); + cell.corner[5] = vol_grid_get(grid, x , y+1, z+1); + cell.corner[6] = vol_grid_get(grid, x+1, y+1, z+1); + cell.corner[7] = vol_grid_get(grid, x+1, y , z+1); + cell.edge[0 ] = verts[EDGE(x ,y ,z ,1)]; + cell.edge[1 ] = verts[EDGE(x ,y+1,z ,0)]; + cell.edge[2 ] = verts[EDGE(x+1,y ,z ,1)]; + cell.edge[3 ] = verts[EDGE(x ,y ,z ,0)]; + cell.edge[4 ] = verts[EDGE(x ,y ,z+1,1)]; + cell.edge[7 ] = verts[EDGE(x ,y ,z+1,0)]; + cell.edge[8 ] = verts[EDGE(x ,y ,z ,2)]; + cell.edge[9 ] = verts[EDGE(x ,y+1,z ,2)]; + cell.edge[10] = verts[EDGE(x+1,y+1,z ,2)]; + cell.edge[11] = verts[EDGE(x+1,y ,z ,2)]; + GList *these = march_one_cube(&cell, level); + if (these) { + verts[EDGE(x ,y+1,z ,0)] = cell.edge[1 ]; + verts[EDGE(x+1,y ,z ,1)] = cell.edge[2 ]; + verts[EDGE(x ,y ,z+1,1)] = cell.edge[4 ]; + verts[EDGE(x ,y+1,z+1,0)] = cell.edge[5 ]; + verts[EDGE(x+1,y ,z+1,1)] = cell.edge[6 ]; + verts[EDGE(x ,y ,z+1,0)] = cell.edge[7 ]; + verts[EDGE(x ,y+1,z ,2)] = cell.edge[9 ]; + verts[EDGE(x+1,y+1,z ,2)] = cell.edge[10]; + verts[EDGE(x+1,y ,z ,2)] = cell.edge[11]; + tris = g_list_concat(these, tris); + } + } } +#undef EDGE + g_free(verts); + g_debug("Marching: generated %d triangles", g_list_length(tris)); + return tris; +} + +void vol_triangle_free(VolTriangle *tri) +{ + for (int i = 0; i < 3; i++) + if (--tri->v[i]->tris == 0) + g_free(tri->v[i]); + g_free(tri); +} + +VolGrid *vol_grid_new(int xs, int ys, int zs) +{ + g_debug("VolGrid: new - %d,%d,%d", xs, ys, zs); + VolGrid *grid = g_new0(VolGrid, 1); + grid->data = g_new0(VolPoint, xs*ys*zs); + grid->xs = xs; + grid->ys = ys; + grid->zs = zs; + return grid; +} + +void vol_grid_free(VolGrid *grid) +{ + g_free(grid->data); + g_free(grid); +} diff --git a/src/objects/marching.h b/src/objects/marching.h new file mode 100644 index 0000000..98ee729 --- /dev/null +++ b/src/objects/marching.h @@ -0,0 +1,82 @@ +/* + * Copyright (C) 2002 Jamie Zawinski + * Copyright (C) 2009-2010 Andy Spencer + * + * Marching cubes implementation based on code from from the xscreensaver + * package. + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program 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 General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#ifndef __VOLUME_H__ +#define __VOLUME_H__ + +#include + +/* VolCoord */ +typedef struct { + double x, y, z; +} VolCoord; + +/* VolPoint */ +typedef struct { + VolCoord c; /* Coordinate of the point */ + double value; /* Function value at point */ +} VolPoint; + +/* VolVertex */ +typedef struct { + VolCoord c; /* Coordinate of the vertex */ + gint tris; /* Count of associated triangles */ + VolCoord norm; /* Vertex normal */ +} VolVertex; + +/* VolTriangle */ +typedef struct { + VolVertex *v[3]; /* Vertices */ + VolCoord norm; /* Surface normal */ +} VolTriangle; + +/* VolCell */ +typedef struct { + VolPoint *corner[8]; /* Corners */ + VolVertex *edge[12]; /* Points along the edges (for caching) */ +} VolCell; + +/* VolGrid */ +typedef struct { + int xs, ys, zs; /* Dimensions of points */ + gpointer data; /* 3-D grid of points */ +} VolGrid; + +/* Find triangles on an isosurface that pass through the cell */ +GList *march_one_cube(VolCell *cell, double level); + +/* Find triangles along an isosurface in a grid of points */ +GList *marching_cubes(VolGrid *grid, double level); + +/* Free/unref data for triangle */ +void vol_triangle_free(VolTriangle *tri); + +/* Grid functions */ +static inline VolPoint *vol_grid_get(VolGrid *grid, int x, int y, int z) +{ + VolPoint (*points)[grid->ys][grid->zs] = grid->data; + return &points[x][y][z]; +} +VolGrid *vol_grid_new(int xs, int ys, int zs); +void vol_grid_free(VolGrid *grid); + + +#endif -- 2.43.2