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..

diff --git a/src/grits.h b/src/grits.h
index 94cc1ab0e0b27ffdfde48bab531d965d3587eae3..688504210c9aea0a6119502575f6a799bf4c12f3 100644 (file)
--- a/src/grits.h
+++ b/src/grits.h
@@ -34,6 +34,7 @@
 #include <objects/grits-tile.h>
 #include <objects/grits-marker.h>
 #include <objects/grits-callback.h>
+#include <objects/grits-volume.h>
 
 /* Plugins */
 #include <grits-plugin.h>

diff --git a/src/objects/Makefile.am b/src/objects/Makefile.am
index 19a3b843f3da3bcf2be6b109c91c3d89be4e4466..bb1a2a4b017f98be94eb753e8e751aaadb5d26c4 100644 (file)
--- 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 (file)
index 0000000..be08279
--- /dev/null
+++ b/src/objects/grits-volume.c
@@ -0,0 +1,175 @@
+/*
+ * Copyright (C) 2009-2010 Andy Spencer <andy753421@gmail.com>
+ *
+ * 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 <http://www.gnu.org/licenses/>.
+ */
+
+/**
+ * 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 <config.h>
+#include <math.h>
+#include <glib.h>
+#include <GL/gl.h>
+#include <gdk/gdkgl.h>
+#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 (file)
index 0000000..6eb4e83
--- /dev/null
+++ b/src/objects/grits-volume.h
@@ -0,0 +1,76 @@
+/*
+ * Copyright (C) 2009-2010 Andy Spencer <andy753421@gmail.com>
+ *
+ * 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 <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef __GRITS_VOLUME_H__
+#define __GRITS_VOLUME_H__
+
+#include <glib.h>
+#include <glib-object.h>
+#include <cairo.h>
+#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 (file)
index 0000000..30a4e3b
--- /dev/null
+++ b/src/objects/marching.c
@@ -0,0 +1,556 @@
+/**
+ * Copyright (C) 2002 Jamie Zawinski <jwz@jwz.org>
+ * Copyright (C) 2009-2010 Andy Spencer <andy753421@gmail.com>
+ *
+ * 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 <pbourke@swin.edu.au>
+ * http://astronomy.swin.edu.au/~pbourke/modelling/polygonise/
+ */
+
+#include <math.h>
+#include <glib.h>
+#include <string.h>
+
+#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 <pbourke@swin.edu.au>
+ * Modified by Andy Spencer <andy753421@gmail.com>
+ */
+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 (file)
index 0000000..98ee729
--- /dev/null
+++ b/src/objects/marching.h
@@ -0,0 +1,82 @@
+/*
+ * Copyright (C) 2002 Jamie Zawinski <jwz@jwz.org>
+ * Copyright (C) 2009-2010 Andy Spencer <andy753421@gmail.com>
+ *
+ * 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 <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef __VOLUME_H__
+#define __VOLUME_H__
+
+#include <glib.h>
+
+/* 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