aboutsummaryrefslogtreecommitdiffstats
path: root/src/contour.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/contour.c')
-rw-r--r--src/contour.c383
1 files changed, 383 insertions, 0 deletions
diff --git a/src/contour.c b/src/contour.c
new file mode 100644
index 0000000..8022adb
--- /dev/null
+++ b/src/contour.c
@@ -0,0 +1,383 @@
+#include "project.h"
+
+#ifndef PI
+#define PI 3.141592654
+#endif
+
+
+typedef struct {
+ int x, y;
+} Point;
+
+typedef struct contour_struct {
+ float *data;
+ float **dcol;
+ unsigned char *flags;
+ unsigned char **fcol;
+ int x, y;
+ int x1, y1;
+ Jwgline current_line;
+ int current_line_size;
+ int lp;
+ Handle h;
+} *Contour;
+
+
+#define INRANGE(l1,l2,c) ( ((c)<(l2)) && ((c)>=(l1)) )
+#define RANGEX(c,x) ( ((x)>0) ? ( ((x)<((c)->x1)) ? (x):((c)->x1)) : 1)
+#define RANGEY(c,y) (((y)>0) ? ( ((y)<((c)->y1)) ? (y):((c)->y1)) : 1)
+
+#define FLAGS(c,p) ((c->fcol[p.y])[p.x])
+#define DATA(c,p) ((c->dcol[p.y])[p.x])
+
+#define DEF_LINE_SIZE 8192
+
+
+
+
+
+void
+cont_reset_point(Contour c)
+{
+ c->current_line.npts = 0;
+}
+
+void
+cont_add_point(Contour c, Jwgpos pos)
+{
+ Jwgpos *ptr;
+ int n;
+
+ if (c->current_line.npts == c->current_line_size) {
+ /* uh-oh we're about to over flow */
+
+ ptr = c->current_line.data;
+
+
+ n = c->current_line_size;
+ printf("JWG: CONTOUR: Increasing maximum line size from %d to %d\n",
+ n, n * 2);
+
+ c->current_line_size = n * 2;
+
+ c->current_line.data = malloc(sizeof(Jwgpos) * n * 2);
+ bcopy(ptr, c->current_line.data, sizeof(Jwgpos) * n);
+
+ free(ptr);
+
+ }
+ jwg_xform(c->h, &pos);
+
+ c->current_line.data[c->current_line.npts] = pos;
+ c->current_line.npts++;
+
+}
+
+unsigned char
+find_direction(Point p1, Point p2)
+{
+ p2.x -= p1.x;
+ p2.y -= p1.y;
+
+ switch (p2.y) {
+ case -1:
+ switch (p2.x) {
+ case -1:
+ return (1);
+ case 0:
+ return (2);
+ case 1:
+ return (4);
+ }
+ case 0:
+ switch (p2.x) {
+ case -1:
+ return (128);
+ case 1:
+ return (8);
+ }
+ case 1:
+ switch (p2.x) {
+ case -1:
+ return (64);
+ case 0:
+ return (32);
+ case 1:
+ return (16);
+ }
+ }
+
+
+ fprintf(stderr, "Contouring cockup type 3 ... (bye bye)\n");
+ abort();
+
+}
+
+Point
+add_direction(Point p, int d)
+{
+ Point dirs[8] = {{-1, -1}, {0, -1}, {1, -1}, {1, 0}, {1, 1}, {0, 1}, {-1, 1}, {-1, 0}};
+
+ p.x = dirs[d].x;
+ p.y = dirs[d].y;
+
+ return (p);
+}
+
+
+int
+points_marked(Contour c, Point p1, Point p2)
+{
+ return (c->fcol[p1.y])[p1.x] & find_direction(p1, p2);
+}
+
+void
+mark_points(Contour c, Point p1, Point p2)
+{
+ (c->fcol[p1.y])[p1.x] |= find_direction(p1, p2);
+ (c->fcol[p2.y])[p2.x] |= find_direction(p2, p1);
+}
+
+
+
+trace(Contour c, Point p1, Point p2, Point pl, float level)
+{
+ Point pn;
+ float v1, v2, vn;
+ Jwgpos pos;
+ float r;
+
+printf(".");
+fflush(stdout);
+
+ if (!c->lp)
+ {
+ pos.x = 1;
+ pos.y = 1;
+
+ cont_add_point(c, pos);
+ }
+
+ v1 = DATA(c, p1);
+ v2 = DATA(c, p2);
+
+ while (1) {
+
+ /*
+ * printf("(%d,%d,%f) (%d,%d,%f) (%d,%d,%f)\n", p1.x, p1.y, DATA(c, p1),
+ * p2.x, p2.y, DATA(c, p2), pl.x, pl.y, DATA(c, pl));
+ */
+
+ r = (v2 - v1);
+ if (r == 0.0) {
+ r = 0.5;
+ } else {
+ r = (level - v1) / r;
+ }
+
+ pos.x = r * (float) RANGEX(c, p2.x);
+ pos.y = r * (float) RANGEY(c, p2.y);
+
+ r = 1.0 - r;
+
+ pos.x += r * (float) RANGEX(c, p1.x);
+ pos.y += r * (float) RANGEY(c, p1.y);
+
+ cont_add_point(c, pos);
+
+
+
+
+
+
+
+ if (DATA(c, p1) > level) {
+ printf("Contouring cockup type 1 -- bye bye\n");
+ abort();
+ }
+ if (DATA(c, p2) <= level) {
+ printf("Contouring cockup type 2 -- bye bye\n");
+ abort();
+ }
+ if (points_marked(c, p1, p2)) {
+ if (c->lp) {
+ jwg_write_line(c->h, &c->current_line);
+ c->current_line.npts = 0;
+ } else {
+ pos.x = 1;
+ pos.y = 1;
+
+ cont_add_point(c, pos);
+ }
+
+
+ return;
+ }
+ mark_points(c, p1, p2);
+
+ pn.x = p1.x + p2.x - pl.x;
+ pn.y = p1.y + p2.y - pl.y;
+
+ vn = DATA(c, pn);
+
+ if (vn <= level) {
+ pl = p1;
+ p1 = pn;
+ v1 = vn;
+ } else {
+ pl = p2;
+ p2 = pn;
+ v2 = vn;
+ }
+
+
+
+
+ }
+
+}
+
+
+
+
+
+possible(Contour c, int x1, int y1, int x2, int y2, int xl, int yl, float level)
+{
+ Point p1, p2, pl;
+
+ p1.x = x1;
+ p1.y = y1;
+ p2.x = x2;
+ p2.y = y2;
+ pl.x = xl;
+ pl.y = yl;
+
+ if (points_marked(c, p1, p2))
+ return;
+
+ trace(c, p1, p2, pl, level);
+
+
+}
+
+
+contour(Contour c, float level)
+{
+ int x, y;
+ float **c0, **c1;
+ float *r00, *r10, *r11, *r01;
+
+ c1 = c0 = c->dcol;
+ c1--;
+
+ for (y = 0; y < c->y; ++y) {
+ r00 = r01 = *(c0);
+ r01--;
+ if (y) {
+ r10 = r11 = *(c1);
+ r11--;
+ }
+ for (x = 0; x < c->x; ++x) {
+
+ if (x) {
+ if (INRANGE(*r00, *r01, level))
+ possible(c, x, y, x - 1, y, x, y - 1, level);
+ if (INRANGE(*r01, *r00, level))
+ possible(c, x - 1, y, x, y, x - 1, y + 1, level);
+ }
+ if (y) {
+ if (INRANGE(*r00, *r10, level))
+ possible(c, x, y - 1, x, y, x - 1, y, level);
+ if (INRANGE(*r10, *r00, level))
+ possible(c, x, y, x, y - 1, x + 1, y - 1, level);
+ }
+ r00++;
+ r10++;
+ r11++;
+ r01++;
+ }
+ c0++;
+ c1++;
+ }
+
+}
+
+
+
+void
+do_contour(Handle h, float level, int side, int lp)
+{
+ Contour c = (Contour) malloc(sizeof(struct contour_struct));
+ int i;
+ float hug = 1.0e10 * (float) side;
+
+
+ c->data = h->data->data;
+ c->x1 = c->x = h->data->w;
+ c->y1 = c->y = h->data->h;
+
+
+ c->dcol = (float **) malloc(sizeof(float *) * c->y);
+
+
+
+ c->dcol[0] = c->data;
+ for (i = 1; i < c->y; ++i)
+ c->dcol[i] = c->dcol[i - 1] + c->x;
+
+ c->flags = (unsigned char *) malloc(c->x * c->y);
+ bzero(c->flags, c->x * c->y);
+
+ c->fcol = (unsigned char **) malloc(c->y * sizeof(unsigned char *));
+
+ c->fcol[0] = c->flags;
+ for (i = 1; i < c->y; ++i)
+ c->fcol[i] = c->fcol[i - 1] + c->x;
+
+ c->x1--;
+ c->y1--;
+
+ for (i = 0; i < c->y; ++i) {
+ (c->dcol[i])[0] = hug;
+ (c->dcol[i])[c->x1] = hug;
+ }
+ for (i = 0; i < c->x; ++i) {
+ (c->dcol[0])[i] = hug;
+ (c->dcol[c->y1])[i] = hug;
+ }
+
+ c->x1--;
+ c->y1--;
+
+ /*
+ * { Point p;
+ *
+ * for (p.y=0;p.y<c->y;++p.y) { for (p.x=0;p.x<c->x;++p.x) {
+ *
+ * printf("%.1f ",DATA(c,p));
+ *
+ * } printf("\n");
+ *
+ * }
+ * } */
+
+ c->current_line.npts = 0;
+ c->current_line_size = DEF_LINE_SIZE;
+ c->current_line.data = (Jwgpos *) malloc(sizeof(Jwgpos) * DEF_LINE_SIZE);
+
+ c->h = h;
+ c->lp = lp;
+
+ contour(c, level);
+ printf("\n");
+
+
+ if (!c->lp)
+ jwg_write_polygon(c->h, &c->current_line);
+
+ free(c->flags);
+ free(c->fcol);
+ free(c->dcol);
+ free(c->current_line.data);
+
+}
generated by cgit v1.2.3 (git 2.25.1) at 2024-05-19 02:34:32 +0000