metaballs effect

diff --git a/src/dynarr.c b/src/dynarr.c
new file mode 100644 (file)
index 0000000..f3e9c47
--- /dev/null
+++ b/src/dynarr.c
@@ -0,0 +1,128 @@
+/* dynarr - dynamic resizable C array data structure
+ * author: John Tsiombikas <nuclear@member.fsf.org>
+ * license: public domain
+ */
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "dynarr.h"
+
+/* The array descriptor keeps auxilliary information needed to manipulate
+ * the dynamic array. It's allocated adjacent to the array buffer.
+ */
+struct arrdesc {
+       int nelem, szelem;
+       int max_elem;
+       int bufsz;      /* not including the descriptor */
+};
+
+#define DESC(x)                ((struct arrdesc*)((char*)(x) - sizeof(struct arrdesc)))
+
+void *dynarr_alloc(int elem, int szelem)
+{
+       struct arrdesc *desc;
+
+       if(!(desc = malloc(elem * szelem + sizeof *desc))) {
+               return 0;
+       }
+       desc->nelem = desc->max_elem = elem;
+       desc->szelem = szelem;
+       desc->bufsz = elem * szelem;
+       return (char*)desc + sizeof *desc;
+}
+
+void dynarr_free(void *da)
+{
+       if(da) {
+               free(DESC(da));
+       }
+}
+
+void *dynarr_resize(void *da, int elem)
+{
+       int newsz;
+       void *tmp;
+       struct arrdesc *desc;
+
+       if(!da) return 0;
+       desc = DESC(da);
+
+       newsz = desc->szelem * elem;
+
+       if(!(tmp = realloc(desc, newsz + sizeof *desc))) {
+               return 0;
+       }
+       desc = tmp;
+
+       desc->nelem = desc->max_elem = elem;
+       desc->bufsz = newsz;
+       return (char*)desc + sizeof *desc;
+}
+
+int dynarr_empty(void *da)
+{
+       return DESC(da)->nelem ? 0 : 1;
+}
+
+int dynarr_size(void *da)
+{
+       return DESC(da)->nelem;
+}
+
+
+/* stack semantics */
+void *dynarr_push(void *da, void *item)
+{
+       struct arrdesc *desc;
+       int nelem;
+
+       desc = DESC(da);
+       nelem = desc->nelem;
+
+       if(nelem >= desc->max_elem) {
+               /* need to resize */
+               struct arrdesc *tmp;
+               int newsz = desc->max_elem ? desc->max_elem * 2 : 1;
+
+               if(!(tmp = dynarr_resize(da, newsz))) {
+                       fprintf(stderr, "failed to resize\n");
+                       return da;
+               }
+               da = tmp;
+               desc = DESC(da);
+               desc->nelem = nelem;
+       }
+
+       if(item) {
+               memcpy((char*)da + desc->nelem++ * desc->szelem, item, desc->szelem);
+       }
+       return da;
+}
+
+void *dynarr_pop(void *da)
+{
+       struct arrdesc *desc;
+       int nelem;
+
+       desc = DESC(da);
+       nelem = desc->nelem;
+
+       if(!nelem) return da;
+
+       if(nelem <= desc->max_elem / 3) {
+               /* reclaim space */
+               struct arrdesc *tmp;
+               int newsz = desc->max_elem / 2;
+
+               if(!(tmp = dynarr_resize(da, newsz))) {
+                       fprintf(stderr, "failed to resize\n");
+                       return da;
+               }
+               da = tmp;
+               desc = DESC(da);
+               desc->nelem = nelem;
+       }
+       desc->nelem--;
+
+       return da;
+}

diff --git a/src/dynarr.h b/src/dynarr.h
new file mode 100644 (file)
index 0000000..e528305
--- /dev/null
+++ b/src/dynarr.h
@@ -0,0 +1,38 @@
+/* dynarr - dynamic resizable C array data structure
+ * author: John Tsiombikas <nuclear@member.fsf.org>
+ * license: public domain
+ */
+#ifndef DYNARR_H_
+#define DYNARR_H_
+
+void *dynarr_alloc(int elem, int szelem);
+void dynarr_free(void *da);
+void *dynarr_resize(void *da, int elem);
+
+int dynarr_empty(void *da);
+int dynarr_size(void *da);
+
+/* stack semantics */
+void *dynarr_push(void *da, void *item);
+void *dynarr_pop(void *da);
+
+
+/* usage example:
+ * -------------
+ * int *arr = dynarr_alloc(0, sizeof *arr);
+ *
+ * int x = 10;
+ * arr = dynarr_push(arr, &x);
+ * x = 5;
+ * arr = dynarr_push(arr, &x);
+ * x = 42;
+ * arr = dynarr_push(arr, &x);
+ *
+ * for(i=0; i<dynarr_size(arr); i++) {
+ *     printf("%d\n", arr[i]);
+ *  }
+ *  dynarr_free(arr);
+ */
+
+
+#endif /* DYNARR_H_ */

diff --git a/src/mcubes.h b/src/mcubes.h
new file mode 100644 (file)
index 0000000..2a9034d
--- /dev/null
+++ b/src/mcubes.h
@@ -0,0 +1,294 @@
+static int mc_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 int mc_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}
+};

diff --git a/src/metaball.c b/src/metaball.c
new file mode 100644 (file)
index 0000000..a6dec99
--- /dev/null
+++ b/src/metaball.c
@@ -0,0 +1,296 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <assert.h>
+#include "screen.h"
+#include "demo.h"
+#include "3dgfx.h"
+#include "gfxutil.h"
+#include "util.h"
+#include "metasurf.h"
+#include "dynarr.h"
+
+struct mesh {
+       int prim;
+       struct g3d_vertex *varr;
+       int16_t *iarr;
+       int vcount, icount;
+};
+
+struct metaball {
+       float energy;
+       float pos[3];
+};
+
+static int init(void);
+static void destroy(void);
+static void start(long trans_time);
+static void draw(void);
+static void draw_mesh(struct mesh *mesh);
+static void zsort(struct mesh *m);
+
+static void calc_voxel_field(void);
+static float eval(struct metasurface *ms, float x, float y, float z);
+static void emit_vertex(struct metasurface *ms, float x, float y, float z);
+
+static struct screen scr = {
+       "metaballs",
+       init,
+       destroy,
+       start, 0,
+       draw
+};
+
+static float cam_theta, cam_phi = 25;
+static float cam_dist = 3;
+static struct mesh mmesh;
+
+static struct metasurface *msurf;
+
+#define VOL_SZ 32
+#define VOL_SCALE      10.0f
+#define VOL_HALF_SCALE (VOL_SCALE * 0.5f)
+static float *volume;
+#define VOXEL(x, y, z) (volume[(z) * VOL_SZ * VOL_SZ + (y) * VOL_SZ + (x)])
+
+#define NUM_MBALLS     3
+static struct metaball mball[NUM_MBALLS];
+
+
+struct screen *metaballs_screen(void)
+{
+       return &scr;
+}
+
+static int init(void)
+{
+       if(!(volume = malloc(VOL_SZ * VOL_SZ * VOL_SZ * sizeof *volume))) {
+               fprintf(stderr, "failed to allocate %dx%dx%d voxel field\n", VOL_SZ, VOL_SZ, VOL_SZ);
+               return -1;
+       }
+
+       mball[0].energy = 1.2;
+       mball[1].energy = 0.8;
+       mball[2].energy = 1.0;
+
+       if(!(msurf = msurf_create())) {
+               fprintf(stderr, "failed to initialize metasurf\n");
+               return -1;
+       }
+       msurf_set_resolution(msurf, VOL_SZ, VOL_SZ, VOL_SZ);
+       msurf_set_bounds(msurf, 0, 0, 0, VOL_SCALE, VOL_SCALE, VOL_SCALE);
+       msurf_eval_func(msurf, eval);
+       msurf_set_threshold(msurf, 0.5);
+       msurf_set_inside(msurf, MSURF_GREATER);
+       msurf_vertex_func(msurf, emit_vertex);
+
+       mmesh.prim = G3D_TRIANGLES;
+       mmesh.varr = 0;
+       mmesh.iarr = 0;
+       mmesh.vcount = mmesh.icount = 0;
+
+       return 0;
+}
+
+static void destroy(void)
+{
+       dynarr_free(mmesh.varr);
+}
+
+static void start(long trans_time)
+{
+       g3d_matrix_mode(G3D_PROJECTION);
+       g3d_load_identity();
+       g3d_perspective(50.0, 1.3333333, 0.5, 100.0);
+
+       g3d_enable(G3D_CULL_FACE);
+       g3d_enable(G3D_LIGHTING);
+       g3d_enable(G3D_LIGHT0);
+
+       g3d_polygon_mode(G3D_WIRE);
+}
+
+static void update(void)
+{
+       static int prev_mx, prev_my;
+       static unsigned int prev_bmask;
+
+       if(mouse_bmask) {
+               if((mouse_bmask ^ prev_bmask) == 0) {
+                       int dx = mouse_x - prev_mx;
+                       int dy = mouse_y - prev_my;
+
+                       if(dx || dy) {
+                               if(mouse_bmask & 1) {
+                                       cam_theta += dx * 1.0;
+                                       cam_phi += dy * 1.0;
+
+                                       if(cam_phi < -90) cam_phi = -90;
+                                       if(cam_phi > 90) cam_phi = 90;
+                               }
+                               if(mouse_bmask & 4) {
+                                       cam_dist += dy * 0.5;
+
+                                       if(cam_dist < 0) cam_dist = 0;
+                               }
+                       }
+               }
+       }
+       prev_mx = mouse_x;
+       prev_my = mouse_y;
+       prev_bmask = mouse_bmask;
+
+       {
+               int i, j;
+               float tsec = time_msec / 1000.0f;
+               static float phase[] = {0.0, M_PI / 3.0, M_PI * 0.8};
+               static float speed[] = {0.8, 1.4, 1.0};
+               static float scale[][3] = {{1, 2, 0.8}, {0.5, 1.6, 0.6}, {1.5, 0.7, 0.5}};
+               static float offset[][3] = {{0, 0, 0}, {0.25, 0, 0}, {-0.2, 0.15, 0.2}};
+
+               for(i=0; i<NUM_MBALLS; i++) {
+                       float t = (tsec + phase[i]) * speed[i];
+
+                       for(j=0; j<3; j++) {
+                               float x = sin(t + j * M_PI / 2.0);
+                               mball[i].pos[j] = offset[i][j] + x * scale[i][j];
+                       }
+               }
+       }
+
+       calc_voxel_field();
+
+       dynarr_free(mmesh.varr);
+       mmesh.varr = dynarr_alloc(0, sizeof *mmesh.varr);
+       msurf_polygonize(msurf);
+}
+
+static void draw(void)
+{
+       update();
+
+       memset(fb_pixels, 0, fb_width * fb_height * 2);
+
+       g3d_matrix_mode(G3D_MODELVIEW);
+       g3d_load_identity();
+       g3d_translate(0, 0, -cam_dist);
+       g3d_rotate(cam_phi, 1, 0, 0);
+       g3d_rotate(cam_theta, 0, 1, 0);
+
+       /*g3d_light_pos(0, -10, 10, 20);*/
+
+       /*zsort(&metamesh);*/
+
+       /*g3d_mtl_diffuse(0.6, 0.6, 0.6);*/
+
+       draw_mesh(&mmesh);
+
+       swap_buffers(fb_pixels);
+}
+
+static void draw_mesh(struct mesh *mesh)
+{
+       if(mesh->iarr) {
+               g3d_draw_indexed(mesh->prim, mesh->varr, mesh->vcount, mesh->iarr, mesh->icount);
+       } else {
+               g3d_draw(mesh->prim, mesh->varr, mesh->vcount);
+       }
+}
+
+static struct {
+       struct g3d_vertex *varr;
+       const float *xform;
+} zsort_cls;
+
+static int zsort_cmp(const void *aptr, const void *bptr)
+{
+       const int16_t *a = (const int16_t*)aptr;
+       const int16_t *b = (const int16_t*)bptr;
+
+       const float *m = zsort_cls.xform;
+
+       const struct g3d_vertex *va = zsort_cls.varr + a[0];
+       const struct g3d_vertex *vb = zsort_cls.varr + b[0];
+
+       float za = m[2] * va->x + m[6] * va->y + m[10] * va->z + m[14];
+       float zb = m[2] * vb->x + m[6] * vb->y + m[10] * vb->z + m[14];
+
+       va = zsort_cls.varr + a[2];
+       vb = zsort_cls.varr + b[2];
+
+       za += m[2] * va->x + m[6] * va->y + m[10] * va->z + m[14];
+       zb += m[2] * vb->x + m[6] * vb->y + m[10] * vb->z + m[14];
+
+       return za - zb;
+}
+
+static void zsort(struct mesh *m)
+{
+       int nfaces = m->icount / m->prim;
+
+       zsort_cls.varr = m->varr;
+       zsort_cls.xform = g3d_get_matrix(G3D_MODELVIEW, 0);
+
+       qsort(m->iarr, nfaces, m->prim * sizeof *m->iarr, zsort_cmp);
+}
+
+static void calc_voxel_field(void)
+{
+       int i, j, k, b;
+       float *voxptr = volume;
+
+       for(i=0; i<VOL_SZ; i++) {
+               float z = VOL_SCALE * ((float)i / (float)VOL_SZ - 0.5);
+
+               for(j=0; j<VOL_SZ; j++) {
+                       float y = VOL_SCALE * ((float)j / (float)VOL_SZ - 0.5);
+
+                       for(k=0; k<VOL_SZ; k++) {
+                               float x = VOL_SCALE * ((float)k / (float)VOL_SZ - 0.5);
+
+                               float val = 0.0f;
+                               for(b=0; b<NUM_MBALLS; b++) {
+                                       float dx = mball[b].pos[0] - x;
+                                       float dy = mball[b].pos[1] - y;
+                                       float dz = mball[b].pos[2] - z;
+
+                                       float lensq = dx * dx + dy * dy + dz * dz;
+
+                                       val += lensq == 0.0f ? 1024.0f : mball[b].energy / lensq;
+                               }
+
+                               *voxptr++ = val;
+                       }
+               }
+       }
+}
+
+static float eval(struct metasurface *ms, float x, float y, float z)
+{
+       int xidx = cround64(x);
+       int yidx = cround64(y);
+       int zidx = cround64(z);
+
+       assert(xidx >= 0 && xidx < VOL_SZ);
+       assert(yidx >= 0 && yidx < VOL_SZ);
+       assert(zidx >= 0 && zidx < VOL_SZ);
+
+       return VOXEL(xidx, yidx, zidx);
+}
+
+static void emit_vertex(struct metasurface *ms, float x, float y, float z)
+{
+       struct g3d_vertex v;
+
+       v.x = x - VOL_HALF_SCALE;
+       v.y = y - VOL_HALF_SCALE;
+       v.z = z - VOL_HALF_SCALE;
+       v.r = cround64(255.0 * x / VOL_SCALE);
+       v.g = cround64(255.0 * y / VOL_SCALE);
+       v.b = cround64(255.0 * z / VOL_SCALE);
+
+       mmesh.varr = dynarr_push(mmesh.varr, &v);
+       assert(mmesh.varr);
+       ++mmesh.vcount;
+}

diff --git a/src/metasurf.c b/src/metasurf.c
new file mode 100644 (file)
index 0000000..85aae1c
--- /dev/null
+++ b/src/metasurf.c
@@ -0,0 +1,451 @@
+/*
+metasurf - a library for implicit surface polygonization
+Copyright (C) 2011-2015  John Tsiombikas <nuclear@member.fsf.org>
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU Lesser 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 Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+#include <stdio.h>
+#include <stdlib.h>
+#include "metasurf.h"
+#include "mcubes.h"
+
+#undef USE_MTETRA
+#define USE_MCUBES
+
+#if (defined(USE_MTETRA) && defined(USE_MCUBES)) || (!defined(USE_MTETRA) && !defined(USE_MCUBES))
+#error "pick either USE_MTETRA or USE_MCUBES, not both..."
+#endif
+
+typedef float vec3[3];
+
+struct metasurface {
+       vec3 min, max;
+       int res[3];
+       float thres;
+
+       msurf_eval_func_t eval;
+       msurf_vertex_func_t vertex;
+       msurf_normal_func_t normal;
+       void *udata;
+
+       float dx, dy, dz;
+       int flip;
+
+       vec3 vbuf[3];
+       int nverts;
+};
+
+static int msurf_init(struct metasurface *ms);
+static void process_cell(struct metasurface *ms, vec3 pos, vec3 sz);
+#ifdef USE_MTETRA
+static void process_tetra(struct metasurface *ms, int *idx, vec3 *pos, float *val);
+#endif
+#ifdef USE_MCUBES
+static void process_cube(struct metasurface *ms, vec3 *pos, float *val);
+#endif
+
+
+struct metasurface *msurf_create(void)
+{
+       struct metasurface *ms;
+
+       if(!(ms = malloc(sizeof *ms))) {
+               return 0;
+       }
+       if(msurf_init(ms) == -1) {
+               free(ms);
+       }
+       return ms;
+}
+
+void msurf_free(struct metasurface *ms)
+{
+       free(ms);
+}
+
+static int msurf_init(struct metasurface *ms)
+{
+       ms->thres = 0.0;
+       ms->eval = 0;
+       ms->vertex = 0;
+       ms->normal = 0;
+       ms->udata = 0;
+       ms->min[0] = ms->min[1] = ms->min[2] = -1.0;
+       ms->max[0] = ms->max[1] = ms->max[2] = 1.0;
+       ms->res[0] = ms->res[1] = ms->res[2] = 40;
+       ms->nverts = 0;
+
+       ms->dx = ms->dy = ms->dz = 0.001;
+       ms->flip = 0;
+
+       return 0;
+}
+
+void msurf_set_user_data(struct metasurface *ms, void *udata)
+{
+       ms->udata = udata;
+}
+
+void *msurf_get_user_data(struct metasurface *ms)
+{
+       return ms->udata;
+}
+
+void msurf_set_inside(struct metasurface *ms, int inside)
+{
+       switch(inside) {
+       case MSURF_GREATER:
+               ms->flip = 0;
+               break;
+
+       case MSURF_LESS:
+               ms->flip = 1;
+               break;
+
+       default:
+               fprintf(stderr, "msurf_inside expects MSURF_GREATER or MSURF_LESS\n");
+       }
+}
+
+int msurf_get_inside(struct metasurface *ms)
+{
+       return ms->flip ? MSURF_LESS : MSURF_GREATER;
+}
+
+void msurf_eval_func(struct metasurface *ms, msurf_eval_func_t func)
+{
+       ms->eval = func;
+}
+
+void msurf_vertex_func(struct metasurface *ms, msurf_vertex_func_t func)
+{
+       ms->vertex = func;
+}
+
+void msurf_normal_func(struct metasurface *ms, msurf_normal_func_t func)
+{
+       ms->normal = func;
+}
+
+void msurf_set_bounds(struct metasurface *ms, float xmin, float ymin, float zmin, float xmax, float ymax, float zmax)
+{
+       ms->min[0] = xmin;
+       ms->min[1] = ymin;
+       ms->min[2] = zmin;
+       ms->max[0] = xmax;
+       ms->max[1] = ymax;
+       ms->max[2] = zmax;
+}
+
+void msurf_get_bounds(struct metasurface *ms, float *xmin, float *ymin, float *zmin, float *xmax, float *ymax, float *zmax)
+{
+       *xmin = ms->min[0];
+       *ymin = ms->min[1];
+       *zmin = ms->min[2];
+       *xmax = ms->max[0];
+       *ymax = ms->max[1];
+       *zmax = ms->max[2];
+}
+
+void msurf_set_resolution(struct metasurface *ms, int xres, int yres, int zres)
+{
+       ms->res[0] = xres;
+       ms->res[1] = yres;
+       ms->res[2] = zres;
+}
+
+void msurf_get_resolution(struct metasurface *ms, int *xres, int *yres, int *zres)
+{
+       *xres = ms->res[0];
+       *yres = ms->res[1];
+       *zres = ms->res[2];
+}
+
+void msurf_set_threshold(struct metasurface *ms, float thres)
+{
+       ms->thres = thres;
+}
+
+float msurf_get_threshold(struct metasurface *ms)
+{
+       return ms->thres;
+}
+
+
+int msurf_polygonize(struct metasurface *ms)
+{
+       int i, j, k;
+       vec3 pos, delta;
+
+       if(!ms->eval || !ms->vertex) {
+               fprintf(stderr, "you need to set eval and vertex callbacks before calling msurf_polygonize\n");
+               return -1;
+       }
+
+       for(i=0; i<3; i++) {
+               delta[i] = (ms->max[i] - ms->min[i]) / (float)ms->res[i];
+       }
+
+       pos[0] = ms->min[0];
+       for(i=0; i<ms->res[0] - 1; i++) {
+               pos[1] = ms->min[1];
+               for(j=0; j<ms->res[1] - 1; j++) {
+
+                       pos[2] = ms->min[2];
+                       for(k=0; k<ms->res[2] - 1; k++) {
+
+                               process_cell(ms, pos, delta);
+
+                               pos[2] += delta[2];
+                       }
+                       pos[1] += delta[1];
+               }
+               pos[0] += delta[0];
+       }
+       return 0;
+}
+
+
+static void process_cell(struct metasurface *ms, vec3 pos, vec3 sz)
+{
+       int i;
+       vec3 p[8];
+       float val[8];
+
+#ifdef USE_MTETRA
+       static int tetra[][4] = {
+               {0, 2, 3, 7},
+               {0, 2, 6, 7},
+               {0, 4, 6, 7},
+               {0, 6, 1, 2},
+               {0, 6, 1, 4},
+               {5, 6, 1, 4}
+       };
+#endif
+
+       static const float offs[][3] = {
+               {0.0f, 0.0f, 0.0f},
+               {1.0f, 0.0f, 0.0f},
+               {1.0f, 1.0f, 0.0f},
+               {0.0f, 1.0f, 0.0f},
+               {0.0f, 0.0f, 1.0f},
+               {1.0f, 0.0f, 1.0f},
+               {1.0f, 1.0f, 1.0f},
+               {0.0f, 1.0f, 1.0f}
+       };
+
+       for(i=0; i<8; i++) {
+               p[i][0] = pos[0] + sz[0] * offs[i][2];
+               p[i][1] = pos[1] + sz[1] * offs[i][1];
+               p[i][2] = pos[2] + sz[2] * offs[i][0];
+
+               val[i] = ms->eval(ms, p[i][0], p[i][1], p[i][2]);
+       }
+
+#ifdef USE_MTETRA
+       for(i=0; i<6; i++) {
+               process_tetra(ms, tetra[i], p, val);
+       }
+#endif
+#ifdef USE_MCUBES
+       process_cube(ms, p, val);
+#endif
+}
+
+
+/* ---- marching cubes implementation ---- */
+#ifdef USE_MCUBES
+
+static unsigned int mc_bitcode(float *val, float thres);
+
+static void process_cube(struct metasurface *ms, vec3 *pos, float *val)
+{
+       static const int pidx[12][2] = {
+               {0, 1}, {1, 2}, {2, 3}, {3, 0}, {4, 5}, {5, 6},
+               {6, 7}, {7, 4}, {0, 4}, {1, 5}, {2, 6}, {3, 7}
+       };
+       int i, j;
+       vec3 vert[12];
+       unsigned int code = mc_bitcode(val, ms->thres);
+
+       if(ms->flip) {
+               code = ~code & 0xff;
+       }
+
+       if(mc_edge_table[code] == 0) {
+               return;
+       }
+
+       for(i=0; i<12; i++) {
+               if(mc_edge_table[code] & (1 << i)) {
+                       int p0 = pidx[i][0];
+                       int p1 = pidx[i][1];
+
+                       float t = (ms->thres - val[p0]) / (val[p1] - val[p0]);
+                       vert[i][0] = pos[p0][0] + (pos[p1][0] - pos[p0][0]) * t;
+                       vert[i][1] = pos[p0][1] + (pos[p1][1] - pos[p0][1]) * t;
+                       vert[i][2] = pos[p0][2] + (pos[p1][2] - pos[p0][2]) * t;
+               }
+       }
+
+       for(i=0; mc_tri_table[code][i] != -1; i+=3) {
+               for(j=0; j<3; j++) {
+                       float *v = vert[mc_tri_table[code][i + j]];
+
+                       if(ms->normal) {
+                               float dfdx, dfdy, dfdz;
+                               dfdx = ms->eval(ms, v[0] - ms->dx, v[1], v[2]) - ms->eval(ms, v[0] + ms->dx, v[1], v[2]);
+                               dfdy = ms->eval(ms, v[0], v[1] - ms->dy, v[2]) - ms->eval(ms, v[0], v[1] + ms->dy, v[2]);
+                               dfdz = ms->eval(ms, v[0], v[1], v[2] - ms->dz) - ms->eval(ms, v[0], v[1], v[2] + ms->dz);
+
+                               if(ms->flip) {
+                                       dfdx = -dfdx;
+                                       dfdy = -dfdy;
+                                       dfdz = -dfdz;
+                               }
+                               ms->normal(ms, dfdx, dfdy, dfdz);
+                       }
+
+                       /* TODO multithreadied polygon emmit */
+                       ms->vertex(ms, v[0], v[1], v[2]);
+               }
+       }
+}
+
+static unsigned int mc_bitcode(float *val, float thres)
+{
+       unsigned int i, res = 0;
+
+       for(i=0; i<8; i++) {
+               if(val[i] > thres) {
+                       res |= 1 << i;
+               }
+       }
+       return res;
+}
+#endif /* USE_MCUBES */
+
+
+/* ---- marching tetrahedra implementation (incomplete) ---- */
+#ifdef USE_MTETRA
+
+static unsigned int mt_bitcode(float v0, float v1, float v2, float v3, float thres);
+static void emmit(struct metasurface *ms, float v0, float v1, vec3 p0, vec3 p1, int rev)
+
+
+#define REVBIT(x)      ((x) & 8)
+#define INV(x)         (~(x) & 0xf)
+#define EDGE(a, b)     emmit(ms, val[idx[a]], val[idx[b]], pos[idx[a]], pos[idx[b]], REVBIT(code))
+static void process_tetra(struct metasurface *ms, int *idx, vec3 *pos, float *val)
+{
+       unsigned int code = mt_bitcode(val[idx[0]], val[idx[1]], val[idx[2]], val[idx[3]], ms->thres);
+
+       switch(code) {
+       case 1:
+       case INV(1):
+               EDGE(0, 1);
+               EDGE(0, 2);
+               EDGE(0, 3);
+               break;
+
+       case 2:
+       case INV(2):
+               EDGE(1, 0);
+               EDGE(1, 3);
+               EDGE(1, 2);
+               break;
+
+       case 3:
+       case INV(3):
+               EDGE(0, 3);
+               EDGE(0, 2);
+               EDGE(1, 3);
+
+               EDGE(1, 3);
+               EDGE(1, 2);
+               EDGE(0, 2);
+               break;
+
+       case 4:
+       case INV(4):
+               EDGE(2, 0);
+               EDGE(2, 1);
+               EDGE(2, 3);
+               break;
+
+       case 5:
+       case INV(5):
+               EDGE(0, 1);
+               EDGE(2, 3);
+               EDGE(0, 3);
+
+               EDGE(0, 1);
+               EDGE(1, 2);
+               EDGE(2, 3);
+               break;
+
+       case 6:
+       case INV(6):
+               EDGE(0, 1);
+               EDGE(1, 3);
+               EDGE(2, 3);
+
+               EDGE(0, 1);
+               EDGE(0, 2);
+               EDGE(2, 3);
+               break;
+
+       case 7:
+       case INV(7):
+               EDGE(3, 0);
+               EDGE(3, 2);
+               EDGE(3, 1);
+               break;
+
+       default:
+               break;  /* cases 0 and 15 */
+       }
+}
+
+#define BIT(i) ((v##i > thres) ? (1 << i) : 0)
+static unsigned int mt_bitcode(float v0, float v1, float v2, float v3, float thres)
+{
+       return BIT(0) | BIT(1) | BIT(2) | BIT(3);
+}
+
+static void emmit(struct metasurface *ms, float v0, float v1, vec3 p0, vec3 p1, int rev)
+{
+       int i;
+       float t = (ms->thres - v0) / (v1 - v0);
+
+       vec3 p;
+       for(i=0; i<3; i++) {
+               p[i] = p0[i] + (p1[i] - p0[i]) * t;
+       }
+       ms->vertex(ms, p[0], p[1], p[2]);
+
+       /*for(i=0; i<3; i++) {
+               ms->vbuf[ms->nverts][i] = p0[i] + (p1[i] - p0[i]) * t;
+       }
+
+       if(++ms->nverts >= 3) {
+               ms->nverts = 0;
+
+               for(i=0; i<3; i++) {
+                       int idx = rev ? (2 - i) : i;
+                       ms->vertex(ms, ms->vbuf[idx][0], ms->vbuf[idx][1], ms->vbuf[idx][2]);
+               }
+       }*/
+}
+#endif /* USE_MTETRA */

diff --git a/src/metasurf.h b/src/metasurf.h
new file mode 100644 (file)
index 0000000..7cdd074
--- /dev/null
+++ b/src/metasurf.h
@@ -0,0 +1,79 @@
+/*
+metasurf - a library for implicit surface polygonization
+Copyright (C) 2011-2015  John Tsiombikas <nuclear@member.fsf.org>
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU Lesser 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 Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef METASURF_H_
+#define METASURF_H_
+
+#define MSURF_GREATER  1
+#define MSURF_LESS             0
+
+struct metasurface;
+
+typedef float (*msurf_eval_func_t)(struct metasurface *ms, float, float, float);
+typedef void (*msurf_vertex_func_t)(struct metasurface *ms, float, float, float);
+typedef void (*msurf_normal_func_t)(struct metasurface *ms, float, float, float);
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+struct metasurface *msurf_create(void);
+void msurf_free(struct metasurface *ms);
+
+void msurf_set_user_data(struct metasurface *ms, void *udata);
+void *msurf_get_user_data(struct metasurface *ms);
+
+/* which is inside above or below the threshold */
+void msurf_set_inside(struct metasurface *ms, int inside);
+int msurf_get_inside(struct metasurface *ms);
+
+/* set a scalar field evaluator function */
+void msurf_eval_func(struct metasurface *ms, msurf_eval_func_t func);
+
+/* set a generated vertex callback function */
+void msurf_vertex_func(struct metasurface *ms, msurf_vertex_func_t func);
+
+/* set a generated surface normal callback function (unused yet) */
+void msurf_normal_func(struct metasurface *ms, msurf_normal_func_t func);
+
+/* set the bounding box (default: -1, -1, -1, 1, 1, 1)
+ * keep this as tight as possible to avoid wasting grid resolution
+ */
+void msurf_set_bounds(struct metasurface *ms, float xmin, float ymin, float zmin, float xmax, float ymax, float zmax);
+void msurf_get_bounds(struct metasurface *ms, float *xmin, float *ymin, float *zmin, float *xmax, float *ymax, float *zmax);
+
+/* resolution of the 3D evaluation grid, the bigger, the better, the slower
+ * (default: 40, 40, 40)
+ */
+void msurf_set_resolution(struct metasurface *ms, int xres, int yres, int zres);
+void msurf_get_resolution(struct metasurface *ms, int *xres, int *yres, int *zres);
+
+/* isosurface threshold value (default: 0) */
+void msurf_set_threshold(struct metasurface *ms, float thres);
+float msurf_get_threshold(struct metasurface *ms);
+
+
+/* finally call this to perform the polygonization */
+int msurf_polygonize(struct metasurface *ms);
+
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif /* METASURF_H_ */

diff --git a/src/screen.c b/src/screen.c
index 9aae073..2b551f4 100644 (file)
--- a/src/screen.c
+++ b/src/screen.c
@@ -12,6 +12,7 @@ struct screen *polytest_screen(void);
 struct screen *plasma_screen(void);
 struct screen *bump_screen(void);
 struct screen *thunder_screen(void);
+struct screen *metaballs_screen(void);
 
 #define NUM_SCR 32
 static struct screen *scr[NUM_SCR];
@@ -45,6 +46,9 @@ int scr_init(void)
        if (!(scr[idx++] = thunder_screen())) {
                return -1;
        }
+       if(!(scr[idx++] = metaballs_screen())) {
+               return -1;
+       }
        num_screens = idx;
 
        assert(num_screens <= NUM_SCR);