add quadtree and basic world model

11 files changed, 666 insertions, 0 deletions

diff --git a/experimental/tectonics/Makefile b/experimental/tectonics/Makefile
new file mode 100644
index 0000000..dd04325
--- /dev/null
+++ b/experimental/tectonics/Makefile
@@ -0,0 +1,12 @@
+CC=gcc
+CFLAGS=-lm `pkg-config --cflags --libs cairo` -ggdb
+DEPS=tectonics.h
+
+%.o: %.c $(DEPS)
+	$(CC) -c -o $@ $< $(CFLAGS)
+
+all: world.o util.o tectonics.o quadtree.o
+		$(CC) -o tectonics world.o util.o tectonics.o quadtree.o $(CFLAGS)
+
+test: all quadtree.test.o tests.o
+		$(CC) -o tests world.o util.o quadtree.o quadtree.test.o tests.o $(CFLAGS)
diff --git a/experimental/tectonics/minunit.h b/experimental/tectonics/minunit.h
new file mode 100644
index 0000000..d62d5e5
--- /dev/null
+++ b/experimental/tectonics/minunit.h
@@ -0,0 +1,44 @@
+#pragma once
+
+#include <stdio.h>
+#include <string.h>
+
+#define UNIT_TEST
+
+#define STR_IMPL(x) #x
+#define STR(x) STR_IMPL(x)
+
+#define MU_INDENT "   "
+
+/* minunit testing macros from /www.jera.com/techinfo/jtns/jtn002.html */
+#define mu_assert_msg(test, message) do { if (!(test)) return message "\n" MU_INDENT " [" __FILE__ ":" STR(__LINE__) "]"; } while (0)
+#define mu_assert(test) mu_assert_msg((test), "assert failed: " #test)
+#define mu_assert_equal(a, b) mu_assert_msg(a == b, "'" #a "' is not equal to '" #b "'")
+#define mu_assert_unequal(a, b) mu_assert_msg(a != b, "'" #a "' is equal to '" #b "'")
+#define mu_assert_streq(a, b) mu_assert_msg(strcmp(a, b) == 0, "'" #a "' is not equal to '" #b "'")
+
+#define mu_run_test(name, test) do {					\
+      const char *message = test();					\
+      tests_run++;							\
+      if (message) {							\
+	 printf(MU_INDENT "test '%s' failed: %s\n", name, message);	\
+	 tests_failed++;						\
+      }									\
+   } while (0)
+#define mu_run_suite(suite) do {					\
+      int run_old = tests_run;						\
+      int failed_old = tests_failed;					\
+      printf("suite: " #suite "\n");					\
+      suite();								\
+      printf(MU_INDENT "ran %d tests, %d failed\n\n",			\
+	     tests_run - run_old,					\
+	     tests_failed - failed_old);				\
+   } while(0)
+#define mu_tests_finished() do {					\
+      printf("ran %d tests, %d failed\n", tests_run, tests_failed);	\
+      return tests_failed;						\
+   } while(0)
+
+#define mu_test const char *
+
+extern int tests_run, tests_failed;
diff --git a/experimental/tectonics/output.png b/experimental/tectonics/output.png
new file mode 100644
index 0000000..041699a
--- /dev/null
+++ b/experimental/tectonics/output.png
diff --git a/experimental/tectonics/quadtree.c b/experimental/tectonics/quadtree.c
new file mode 100644
index 0000000..dead9f2
--- /dev/null
+++ b/experimental/tectonics/quadtree.c
@@ -0,0 +1,135 @@
+#include <math.h>
+#include "tectonics.h"
+
+bool contains_point(struct quad_region_t region, struct point_t pt)
+{
+   if (-region.half_dim < (region.center.x - pt.x) <= region.half_dim)
+      return false;
+
+   if (-region.half_dim < (region.center.y - pt.y) <= region.half_dim)
+      return false;
+
+   return true;
+}
+
+
+struct quadtree_node_t new_node(struct point_t center, double half_dim)
+{
+   struct quadtree_node_t node;
+
+   node.region = (struct quad_region_t) { center, half_dim };
+   node.id = -1;
+   node.nw = NULL;
+   node.ne = NULL;
+   node.sw = NULL;
+   node.se = NULL;
+
+   return node;
+}
+
+
+void subdivide(struct quadtree_node_t *node, struct point_t *points)
+{
+   node->nw = malloc(sizeof(struct quadtree_node_t));
+   node->ne = malloc(sizeof(struct quadtree_node_t));
+   node->sw = malloc(sizeof(struct quadtree_node_t));
+   node->se = malloc(sizeof(struct quadtree_node_t));
+
+   struct point_t center = node->region.center;
+   double quarter_dim = node->region.half_dim / 2;
+
+   struct point_t c_nw =
+      { center.x - quarter_dim, center.y - quarter_dim };
+   *(node->nw) = new_node(c_nw, quarter_dim);
+
+   struct point_t c_ne =
+      { center.x + quarter_dim, center.y - quarter_dim };
+   *(node->ne) = new_node(c_ne, quarter_dim);
+
+   struct point_t c_sw =
+      { center.x - quarter_dim, center.y + quarter_dim };
+   *(node->sw) = new_node(c_sw, quarter_dim);
+
+   struct point_t c_se =
+      { center.x + quarter_dim, center.y + quarter_dim };
+   *(node->se) = new_node(c_se, quarter_dim);
+
+   if (node->id != -1) {
+   /* move point to a child */
+      if (contains_point(node->nw->region, points[node->id]))
+	 node->nw->id = node->id;
+      else if (contains_point(node->ne->region, points[node->id]))
+	 node->ne->id = node->id;
+      else if (contains_point(node->sw->region, points[node->id]))
+	 node->sw->id = node->id;
+      else
+	 node->se->id = node->id;
+      node->id = -1;
+   }
+}
+
+bool insert(struct quadtree_node_t *node,
+	    struct point_t *points, int id)
+{
+   if (!contains_point(node->region, points[id]))
+      return false;
+
+   if (node->id == -1 && node->nw == NULL) {
+      node->id = id;
+      return true;
+   }
+
+   if (node->nw == NULL)
+      subdivide(node, points);
+
+   if (insert(node->nw, points, id)) return true;
+   if (insert(node->ne, points, id)) return true;
+   if (insert(node->sw, points, id)) return true;
+   if (insert(node->se, points, id)) return true;
+
+   return false;
+}
+
+
+int get_closest(struct quadtree_node_t *node,
+		struct point_t pt)
+{
+   if (!contains_point(node->region, pt))
+      return -1;
+
+   if (node->id != -1)
+      return node->id;
+
+   int closest;
+
+   if ((closest = get_closest(node->nw, pt)) != -1) return closest;
+   if ((closest = get_closest(node->ne, pt)) != -1) return closest;
+   if ((closest = get_closest(node->sw, pt)) != -1) return closest;
+   if ((closest = get_closest(node->se, pt)) != -1) return closest;
+
+   return -1;
+}
+
+
+void draw_quadtree(cairo_t *cr, struct quadtree_node_t *node)
+{
+   int width, height;
+   get_cairo_size(cr, &width, &height);
+
+   cairo_set_source_rgba(cr, 0, 0, 1, 1);
+   struct point_t center = node->region.center;
+   double half_dim = node->region.half_dim;
+
+   double x = (center.x - half_dim) * width;
+   double y = (center.y - half_dim) * height;
+   
+   cairo_rectangle(cr, x, y, 2*half_dim * width, 2*half_dim * height);
+   cairo_stroke(cr);
+
+   if (node->nw != NULL) {
+      draw_quadtree(cr, node->nw);
+      draw_quadtree(cr, node->ne);
+      draw_quadtree(cr, node->sw);
+      draw_quadtree(cr, node->se);
+   }
+}
diff --git a/experimental/tectonics/quadtree.test.c b/experimental/tectonics/quadtree.test.c
new file mode 100644
index 0000000..273087b
--- /dev/null
+++ b/experimental/tectonics/quadtree.test.c
@@ -0,0 +1,225 @@
+#include "minunit.h"
+#include "tests.h"
+#include "tectonics.h"
+
+mu_test test_contains_point();
+mu_test test_subdivide();
+mu_test test_insert();
+mu_test test_insert_two();
+mu_test test_get_closest();
+
+void quadtree_tests()
+{
+   mu_run_test("check if contains_point() works",
+	       test_contains_point);
+   mu_run_test("subdivide quadtree node",
+	       test_subdivide);
+   mu_run_test("insert second point",
+	       test_insert);
+   mu_run_test("insert two more points",
+	       test_insert_two);
+   mu_run_test("look up nearest point",
+	       test_get_closest);
+}
+
+
+/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
+
+mu_test test_contains_point()
+{
+   struct quad_region_t region;
+   region.center = (struct point_t){ 0.5, 0.5 };
+   region.half_dim = 0.5;
+
+   mu_assert(contains_point(region, (struct point_t){ 0.5, 0.5 }));
+   mu_assert(contains_point(region, (struct point_t){ 0.25, 0.5 }));
+   mu_assert(contains_point(region, (struct point_t){ 0.5, 0.25 }));
+   mu_assert(contains_point(region, (struct point_t){ 0, 0 }));
+
+
+   mu_assert(!contains_point(region, (struct point_t){ 1.25, 0.5 }));
+   mu_assert(!contains_point(region, (struct point_t){ 0.25, 1.5 }));
+   mu_assert(!contains_point(region, (struct point_t){ 1, 1 }));
+
+   return 0;
+}
+
+
+mu_test test_subdivide()
+{
+   struct quadtree_node_t node;
+   struct point_t pts[1] = { { 0.75, 0.75 } };
+
+   node.region = (struct quad_region_t){ { 0.5, 0.5 }, 0.5 };
+   node.id = 0;
+   node.nw = NULL;
+   node.ne = NULL;
+   node.sw = NULL;
+   node.se = NULL;
+
+   subdivide(&node, pts);
+
+   mu_assert_equal(node.id, -1);
+
+   mu_assert_equal(node.nw->region.center.x, 0.25);
+   mu_assert_equal(node.nw->region.center.y, 0.25);
+   mu_assert_equal(node.nw->region.half_dim, 0.25);
+   mu_assert_equal(node.nw->nw, NULL);
+   mu_assert_equal(node.nw->id, -1);
+
+   mu_assert_equal(node.ne->region.center.x, 0.75);
+   mu_assert_equal(node.ne->region.center.y, 0.25);
+   mu_assert_equal(node.ne->region.half_dim, 0.25);
+   mu_assert_equal(node.ne->nw, NULL);
+   mu_assert_equal(node.ne->id, -1);
+
+   mu_assert_equal(node.sw->region.center.x, 0.25);
+   mu_assert_equal(node.sw->region.center.y, 0.75);
+   mu_assert_equal(node.sw->region.half_dim, 0.25);
+   mu_assert_equal(node.sw->nw, NULL);
+   mu_assert_equal(node.sw->id, -1);
+
+   mu_assert_equal(node.se->region.center.x, 0.75);
+   mu_assert_equal(node.se->region.center.y, 0.75);
+   mu_assert_equal(node.se->region.half_dim, 0.25);
+   mu_assert_equal(node.se->nw, NULL);
+   mu_assert_equal(node.se->id, 0);
+
+   return 0;
+}
+
+
+mu_test test_insert()
+{
+   struct quadtree_node_t node;
+   struct point_t pts[2] = { {0.75, 0.75}, {0.25, 0.25} };
+
+   node.region = (struct quad_region_t){ { 0.5, 0.5 }, 0.5 };
+   node.id = 0;
+   node.nw = NULL;
+   node.ne = NULL;
+   node.sw = NULL;
+   node.se = NULL;
+
+   insert(&node, pts, 1);
+
+   mu_assert_equal(node.id, -1);
+
+   mu_assert_equal(node.nw->region.center.x, 0.25);
+   mu_assert_equal(node.nw->region.center.y, 0.25);
+   mu_assert_equal(node.nw->region.half_dim, 0.25);
+   mu_assert_equal(node.nw->nw, NULL);
+   mu_assert_equal(node.nw->id, 1);
+
+   mu_assert_equal(node.ne->region.center.x, 0.75);
+   mu_assert_equal(node.ne->region.center.y, 0.25);
+   mu_assert_equal(node.ne->region.half_dim, 0.25);
+   mu_assert_equal(node.ne->nw, NULL);
+   mu_assert_equal(node.ne->id, -1);
+
+   mu_assert_equal(node.sw->region.center.x, 0.25);
+   mu_assert_equal(node.sw->region.center.y, 0.75);
+   mu_assert_equal(node.sw->region.half_dim, 0.25);
+   mu_assert_equal(node.sw->nw, NULL);
+   mu_assert_equal(node.sw->id, -1);
+
+   mu_assert_equal(node.se->region.center.x, 0.75);
+   mu_assert_equal(node.se->region.center.y, 0.75);
+   mu_assert_equal(node.se->region.half_dim, 0.25);
+   mu_assert_equal(node.se->nw, NULL);
+   mu_assert_equal(node.se->id, 0);
+
+   return 0;
+}
+
+
+mu_test test_insert_two()
+{
+   struct quadtree_node_t node;
+   struct point_t pts[3] =
+      { {0.75, 0.75}, {0.3, 0.3}, {0.2, 0.2} };
+
+   node.region = (struct quad_region_t){ { 0.5, 0.5 }, 0.5 };
+   node.id = 0;
+   node.nw = NULL;
+   node.ne = NULL;
+   node.sw = NULL;
+   node.se = NULL;
+
+   insert(&node, pts, 1);
+   insert(&node, pts, 2);
+
+   mu_assert_equal(node.id, -1);
+
+   mu_assert_equal(node.nw->region.center.x, 0.25);
+   mu_assert_equal(node.nw->region.center.y, 0.25);
+   mu_assert_equal(node.nw->region.half_dim, 0.25);
+   mu_assert_equal(node.nw->id, -1);
+
+   /*printf("%d, %d, %d, %d\n",
+	  node.nw->nw->id,
+	  node.nw->ne->id,
+	  node.nw->sw->id,
+	  node.nw->se->id);*/
+
+   mu_assert_equal(node.nw->nw->id,  2);
+   mu_assert_equal(node.nw->nw->nw,  NULL);
+   mu_assert_equal(node.nw->ne->id, -1);
+   mu_assert_equal(node.nw->ne->nw,  NULL);
+   mu_assert_equal(node.nw->sw->id, -1);
+   mu_assert_equal(node.nw->sw->nw,  NULL);
+   mu_assert_equal(node.nw->se->id,  1);
+   mu_assert_equal(node.nw->se->nw,  NULL);
+
+   mu_assert_equal(node.ne->region.center.x, 0.75);
+   mu_assert_equal(node.ne->region.center.y, 0.25);
+   mu_assert_equal(node.ne->region.half_dim, 0.25);
+   mu_assert_equal(node.ne->nw, NULL);
+   mu_assert_equal(node.ne->id, -1);
+
+   mu_assert_equal(node.sw->region.center.x, 0.25);
+   mu_assert_equal(node.sw->region.center.y, 0.75);
+   mu_assert_equal(node.sw->region.half_dim, 0.25);
+   mu_assert_equal(node.sw->nw, NULL);
+   mu_assert_equal(node.sw->id, -1);
+
+   mu_assert_equal(node.se->region.center.x, 0.75);
+   mu_assert_equal(node.se->region.center.y, 0.75);
+   mu_assert_equal(node.se->region.half_dim, 0.25);
+   mu_assert_equal(node.se->nw, NULL);
+   mu_assert_equal(node.se->id, 0);
+
+   return 0;
+}
+
+
+mu_test test_get_closest()
+{
+   struct quadtree_node_t node;
+   struct point_t pts[3] =
+      { {0.75, 0.75}, {0.3, 0.3}, {0.2, 0.2} };
+
+   node.region = (struct quad_region_t){ { 0.5, 0.5 }, 0.5 };
+   node.id = 0;
+   node.nw = NULL;
+   node.ne = NULL;
+   node.sw = NULL;
+   node.se = NULL;
+
+   insert(&node, pts, 1);
+   insert(&node, pts, 2);
+
+   int closest_75 = get_closest(&node, (struct point_t){0.75, 0.75});
+   mu_assert_equal(closest_75, 0);
+
+   int closest_33 = get_closest(&node, (struct point_t){0.3, 0.3});
+   mu_assert_equal(closest_33, 1);
+
+   int closest_22 = get_closest(&node, (struct point_t){0.2, 0.2});
+   mu_assert_equal(closest_22, 2);
+
+   int closest_nil = get_closest(&node, (struct point_t){100, 100});
+   mu_assert_equal(closest_nil, -1);
+
+   return 0;
+}
diff --git a/experimental/tectonics/tectonics.c b/experimental/tectonics/tectonics.c
new file mode 100644
index 0000000..f294944
--- /dev/null
+++ b/experimental/tectonics/tectonics.c
@@ -0,0 +1,24 @@
+#include <cairo.h>
+
+#include "tectonics.h"
+
+#define X_RES 512
+#define Y_RES X_RES
+
+int main()
+{
+   cairo_surface_t *surface = cairo_image_surface_create (CAIRO_FORMAT_ARGB32, X_RES, Y_RES);
+   cairo_t *cr = cairo_create(surface);
+
+   struct world_t world;
+   create_world(&world, 1000);
+   if (world.points == NULL)
+      return 1;
+   render_world(cr, &world);
+   free_world(&world);
+   
+   cairo_destroy(cr);
+   cairo_surface_write_to_png(surface, "output.png");
+   cairo_surface_destroy(surface);
+   return 0;
+}
diff --git a/experimental/tectonics/tectonics.h b/experimental/tectonics/tectonics.h
new file mode 100644
index 0000000..5f39ff7
--- /dev/null
+++ b/experimental/tectonics/tectonics.h
@@ -0,0 +1,76 @@
+#ifndef TECTONICS_H
+#define TECTONICS_H
+
+#include <stdlib.h>
+#include <stdbool.h>
+#include <cairo.h>
+
+struct point_t {
+   double x, y;
+   void *data;
+};
+
+
+struct quad_region_t {
+   struct point_t center;
+   double half_dim;
+};
+
+
+struct quadtree_node_t {
+   struct quad_region_t region;
+
+   int id;
+   
+   /* children */
+   struct quadtree_node_t *nw, *ne, *sw, *se;
+};
+
+
+struct world_t {
+   struct point_t *points;
+   struct quadtree_node_t tree;
+   int n_points;
+};
+
+/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ *
+ * quadtree
+ *
+ * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ */
+
+bool contains_point(struct quad_region_t region, struct point_t pt);
+struct quadtree_node_t new_node(struct point_t center, double half_dim);
+void subdivide(struct quadtree_node_t *node, struct point_t *points);
+bool insert(struct quadtree_node_t *node,
+	    struct point_t *points, int id);
+int get_closest(struct quadtree_node_t *root,
+		struct point_t pt);
+void draw_quadtree(cairo_t *cr, struct quadtree_node_t *root);
+
+/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ *
+ * util
+ *
+ * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ */
+
+void get_cairo_size(cairo_t *cr, int *width, int *height);
+double rand01();
+double distance(double x0, double y0, double x1, double y1);
+
+
+/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ *
+ * worlds
+ *
+ * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ */
+
+void create_world(struct world_t *world, int n_points);
+void free_world(struct world_t *world);
+void render_world(cairo_t *cr, struct world_t *world);
+void run_world_step(struct world_t *world);
+
+#endif
diff --git a/experimental/tectonics/tests.c b/experimental/tectonics/tests.c
new file mode 100644
index 0000000..5e00525
--- /dev/null
+++ b/experimental/tectonics/tests.c
@@ -0,0 +1,11 @@
+#include "tests.h"
+#include "minunit.h"
+
+int tests_run = 0;
+int tests_failed = 0;
+
+int main() {
+   RUN_TESTS();
+   mu_tests_finished();
+   return tests_failed;
+}
diff --git a/experimental/tectonics/tests.h b/experimental/tectonics/tests.h
new file mode 100644
index 0000000..c49c46f
--- /dev/null
+++ b/experimental/tectonics/tests.h
@@ -0,0 +1,9 @@
+#ifndef TESTS_H
+#define TESTS_H
+
+void quadtree_tests();
+
+#define RUN_TESTS()				\
+   mu_run_suite(quadtree_tests);
+
+#endif
diff --git a/experimental/tectonics/util.c b/experimental/tectonics/util.c
new file mode 100644
index 0000000..a0a32a3
--- /dev/null
+++ b/experimental/tectonics/util.c
@@ -0,0 +1,25 @@
+#include <math.h>
+
+#include "tectonics.h"
+
+void get_cairo_size(cairo_t *cr, int *width, int *height)
+{
+   cairo_surface_t *surface = cairo_get_target(cr);
+   *width = cairo_image_surface_get_width(surface);
+   *height = cairo_image_surface_get_height(surface);
+}
+
+
+double rand01()
+{
+   return ((double) rand())/RAND_MAX;
+}
+
+
+double distance(double x0, double y0, double x1, double y1)
+{
+   double dx = x1-x0;
+   double dy = y1-y0;
+
+   return sqrt(dx*dx + dy*dy);
+}
diff --git a/experimental/tectonics/world.c b/experimental/tectonics/world.c
new file mode 100644
index 0000000..cfc5d14
--- /dev/null
+++ b/experimental/tectonics/world.c
@@ -0,0 +1,105 @@
+#include <stdio.h>
+#include <math.h>
+
+#include "tectonics.h"
+
+struct centroid_t {
+   double x, y;
+   int updates;
+};
+
+static int closest_point(struct world_t *world, double x, double y)
+{
+   double dist = 1e9;
+   int closest = -1;
+   for (int i=0; i<world->n_points; i++) {
+      struct point_t *pt = world->points + i;
+      double d = distance(x, y, pt->x, pt->y);
+      if (d <= dist) {
+	 dist = d;
+	 closest = i;
+      }
+   }
+   return closest;
+}
+
+static void relax_points(struct world_t *world, int iterations)
+{
+   struct centroid_t *centroid =
+      malloc(sizeof(struct centroid_t) * world->n_points);
+   
+   for (int i=0; i<world->n_points; i++) {
+      struct point_t *pt = world->points + i;
+      struct centroid_t *c = centroid + i;
+      c->x = pt->x;
+      c->y = pt->y;
+      c->updates = 1;
+   }
+
+   for (int i=0; i<iterations; i++) {
+      printf("%02f%%\n", 100*((double)i+1)/iterations);
+      double x = rand01();
+      double y = rand01();
+      int closest = closest_point(world, x, y);
+
+      int t = centroid[closest].updates;
+      centroid[closest].x = (t*centroid[closest].x + x)/(t+1);
+      centroid[closest].y = (t*centroid[closest].y + y)/(t+1);
+      centroid[closest].updates += 1;
+   }
+
+   for (int i=0; i<world->n_points; i++) {
+      struct point_t *pt = world->points + i;
+      struct centroid_t *c = centroid + i;
+      pt->x = c->x;
+      pt->y = c->y;
+   }
+}
+
+void create_world(struct world_t *world, int n_points)
+{
+   world->n_points = n_points;
+   world->points = malloc(sizeof(struct point_t) * n_points);
+   if (world->points == NULL) {
+      fprintf(stderr,
+	      "ERROR: failed to allocate %d points\n",
+	      n_points);
+      return;
+   }
+
+   struct point_t center = { 0.5, 0.5 };
+   world->tree = new_node(center, 0.5);
+
+   for (int i=0; i<n_points; i++) {
+      struct point_t *pt = world->points + i;
+      pt->x = rand01();
+      pt->y = rand01();
+
+      insert(&(world->tree), world->points, i);
+   }
+}
+
+
+void free_world(struct world_t *world)
+{
+   free(world->points);
+}
+
+
+void render_world(cairo_t *cr, struct world_t *world)
+{
+   int width, height;
+   get_cairo_size(cr, &width, &height);
+   double r = ((double)width) / (100*sqrt(world->n_points));
+   cairo_set_source_rgba(cr, 1, 0, 0, 1);
+
+   for (int i=0; i<world->n_points; i++) {
+      struct point_t *pt = world->points + i;
+      double xc = pt->x * width;
+      double yc = pt->y * height;
+      cairo_arc(cr, xc, yc, r, 0, 2*M_PI);
+      cairo_stroke(cr);
+   }
+
+   draw_quadtree(cr, &(world->tree));
+}