summaryrefslogtreecommitdiff
path: root/meme/slice.c
diff options
context:
space:
mode:
authorWerner Almesberger <werner@almesberger.net>2015-03-14 11:10:08 (GMT)
committerWerner Almesberger <werner@almesberger.net>2015-03-14 11:10:08 (GMT)
commit5063fcc259c3f5c485ade53ea2a1b1882bec2d30 (patch)
tree579fa96ebac862b54651c8f0aad219c43ae0a4c4 /meme/slice.c
parent7d70682e30d87204e1e441c7f6c7ecbd68fe1a07 (diff)
downloadmisc-5063fcc259c3f5c485ade53ea2a1b1882bec2d30.zip
misc-5063fcc259c3f5c485ade53ea2a1b1882bec2d30.tar.gz
misc-5063fcc259c3f5c485ade53ea2a1b1882bec2d30.tar.bz2
meme/: new option -s slice.gp for intersect the mesh with a plane (WIP)
Diffstat (limited to 'meme/slice.c')
-rw-r--r--meme/slice.c194
1 files changed, 194 insertions, 0 deletions
diff --git a/meme/slice.c b/meme/slice.c
new file mode 100644
index 0000000..8bfbd1e
--- /dev/null
+++ b/meme/slice.c
@@ -0,0 +1,194 @@
+/*
+ * slice.c - Intersect the mesh with a plane
+ *
+ * Written 2015 by Werner Almesberger
+ * Copyright 2015 Werner Almesberger
+ *
+ * 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 2 of the License, or
+ * (at your option) any later version.
+ */
+
+
+#include <stdbool.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+
+#include "util.h"
+#include "mesh.h"
+#include "slice.h"
+
+
+#define EPSILON 1e-6
+
+
+struct point {
+ double x, y, z;
+ struct point *next;
+};
+
+static struct point *points = NULL;
+
+
+/* ----- Slicing plane ----------------------------------------------------- */
+
+
+static double z0, cx, cy; /* z = cx * x + cy * y + z0 */
+
+
+
+static void slice_plane(void)
+{
+ /* @@@ */
+ cx = cy = 0;
+ z0 = points->z;
+}
+
+
+/* ----- Slice the mesh ---------------------------------------------------- */
+
+
+static bool above(const struct vertex *v)
+{
+ return v->z > z0 + cx * v->x + cy * v->y;
+}
+
+
+static bool solve(const struct vertex *va, const struct vertex *vb, double *res)
+{
+ double dx, dy, dz;
+ double m, a;
+
+ dx = vb->x - va->x;
+ dy = vb->y - va->y;
+ dz = vb->z - va->z;
+
+ /*
+ * We solve this equation for t:
+ *
+ * x(t) = t * dx + va->x;
+ * y(t) = t * dy + va->y;
+ * z(t) = t * dz + va->z;
+ * cx * x(t) + cy * y(t) + z0 = z(t)
+ *
+ * Calculate m and a such that
+ * m * t + a = 0
+ */
+
+ m = cx * dx + cx * dy - dz;
+ a = cx * va->x + cy * va->y - va->z + z0;
+
+ if (fabs(m) < EPSILON)
+ return 0;
+
+ *res = -a / m;
+ return 1;
+}
+
+static void calc(void (*fn)(double x, double y, double z, bool last),
+ const struct vertex *va, const struct vertex *vb, double t,
+ bool needle, bool last)
+{
+ double dx, dy, dz;
+
+ dx = vb->x - va->x;
+ dy = vb->y - va->y;
+ dz = vb->z - va->z;
+
+ /* @@@ needle correction */
+
+ fn((va->x + t * dx) / 1000, (va->y + t * dy) / 1000,
+ (va->z + t * dz) / 1000, last);
+}
+
+
+static void segment(void (*fn)(double x, double y, double z, bool last),
+ const struct vertex *va, const struct vertex *vb, const struct vertex *vc,
+ bool needle)
+{
+ double tb, tc;
+
+ if (solve(va, vb, &tb) && solve(va, vc, &tc)) {
+ calc(fn, va, vb, tb, needle, 0);
+ calc(fn, va, vc, tc, needle, 1);
+ }
+}
+
+
+void slice_run(void (*fn)(double x, double y, double z, bool last), bool needle)
+{
+ const struct facet *f;
+ bool a, b, c;
+
+ slice_plane();
+ for (f = facets; f; f = f->next) {
+ a = above(f->v[0]);
+ b = above(f->v[1]);
+ c = above(f->v[2]);
+ if (a && b && c)
+ continue;
+ if (!a && !b && !c)
+ continue;
+ if ((b && c) || (!b && !c))
+ segment(fn, f->v[0], f->v[1], f->v[2], needle);
+ if ((a && c) || (!a && !c))
+ segment(fn, f->v[1], f->v[0], f->v[2], needle);
+ if ((a && b) || (!a && !b))
+ segment(fn, f->v[2], f->v[0], f->v[1], needle);
+ }
+}
+
+
+/* ----- Points that define the slicing plane ------------------------------ */
+
+
+static void slice_load_file(FILE *file)
+{
+ struct point *p;
+ int lineno = 0;
+ char buf[1024];
+ int n;
+ double zoff;
+
+ p = alloc_type(struct point);
+ while (fgets(buf, sizeof(buf), file)) {
+ lineno++;
+ if (*buf == '#')
+ continue;
+ n = sscanf(buf, "%lf %lf %lf %lf\n",
+ &p->x, &p->y, &p->z, &zoff);
+ switch (n) {
+ case -1:
+ continue;
+ case 4:
+ p->z += zoff;
+ /* fall through */
+ case 3:
+ break;
+ default:
+ fprintf(stderr, "invalid data at line %d\n", lineno);
+ exit(1);
+ }
+
+ p->next = points;
+ points = p;
+ p = alloc_type(struct point);
+ }
+ free(p);
+}
+
+
+void slice_load(const char *name)
+{
+ FILE *file;
+
+ file = fopen(name, "r");
+ if (!file) {
+ perror(name);
+ exit(1);
+ }
+ slice_load_file(file);
+ fclose(file);
+}