#ifndef   ISOCURVES_H
#define   ISOCURVES_H
#include <stdlib.h>	/* For size_t */

typedef struct point	 point_t;
struct point {
	double		 x;
	double		 y;
#ifdef USE_3D
	double		 z;
#endif
};

typedef struct tri	 tri_t;
struct tri {
	size_t		 p[3];
};

typedef struct path	 path_t;
struct path {
	struct path	*next;
	size_t		 points;
	struct point	 point[];
};

path_t *isocurves(const size_t	points,
                  const point_t	point[],
                  const double	value[],
                  const size_t	tris,
                  const tri_t	tri[],
                  const double	level);

static inline point_t between(const point_t p0, const point_t p1, const double t)
{
	if (t <= 0.0)
		return p0;

	if (t >= 1.0)
		return p1;

#ifdef USE_3D
	return (struct point){ .x = p0.x + t * (p1.x - p0.x),
	                       .y = p0.y + t * (p1.y - p0.y),
	                       .z = p0.z + t * (p1.z - p0.z) };
#else
	return (struct point){ .x = p0.x + t * (p1.x - p0.x),
	                       .y = p0.y + t * (p1.y - p0.y) };
#endif
}

#endif /* ISOCURVES_H */

#include <stdlib.h>
#include <errno.h>
#include "isocurves.h"

struct edge {
	point_t		 c;		/* Intersected point */
	size_t		 p[2];		/* Point indices identifying this edge */
	size_t		 e;		/* Next edge along the contour curve */
	signed char	 s;		/* Status: 3 = start, 2 = tracing, 1 = unused, 0 = used. */
};

path_t *isocurves(const size_t	 points,
                  const point_t	 point[],
                  const double	 value[],
                  const size_t	 tris,
                  const tri_t	 tri[],
                  const double	 level)
{
	size_t		 edges_max = 0;
	size_t		 edges = 0;
	struct edge	*edge = NULL;

	path_t		*root = NULL, *path;

	size_t		 t, e, i, n;


	/* We need at least one triangle and three points. */
	if (points < (size_t)3 || tris < (size_t)1) {
		errno = EINVAL;
		return NULL;
	}

	/* Check that all triangle points are valid. */
	t = tris; while (t--)
		if (tri[t].p[0] >= points ||
		    tri[t].p[1] >= points ||
		    tri[t].p[2] >= points ||
		    tri[t].p[0] == tri[t].p[1] ||
		    tri[t].p[0] == tri[t].p[2] ||
		    tri[t].p[1] == tri[t].p[2]) {
			errno = EDOM;
			return NULL;
		}

	/* Find triangles bisected by the isocurve. */
	t = tris; while (t--) {
		const size_t	p[3] = { tri[t].p[0], tri[t].p[1], tri[t].p[2] };

		const int	type = ( (value[p[0]] >= level) ? 1 : 0 )
		         	     | ( (value[p[1]] >= level) ? 2 : 0 )
		         	     | ( (value[p[2]] >= level) ? 4 : 0 );

		size_t		curr_p[2], curr_e;
		size_t		next_p[2], next_e;

		/* Nothing in this triangle? */
		if (type < 1 || type > 6)
			continue;

		/* We know there is a contour line bisecting this triangle,
		 * since the three points straddle the contour level.
		 * Orient the line segment so that the upper half is on the right.
		*/

		/* Determine which edge starts the contour line segment
		 * passing through this triangle. For efficiency,
		 * pick the points so that curr_p[0] < curr_p[1].
		*/
		switch (type) {
		case 1: case 5:
			if (p[0] < p[1]) {
				curr_p[0] = p[0];
				curr_p[1] = p[1];
			} else {
				curr_p[0] = p[1];
				curr_p[1] = p[0];
			}
			break;
		case 2: case 3:
			if (p[1] < p[2]) {
				curr_p[0] = p[1];
				curr_p[1] = p[2];
			} else {
				curr_p[0] = p[2];
				curr_p[1] = p[1];
			}
			break;
		case 4: case 6:
			if (p[0] < p[2]) {
				curr_p[0] = p[0];
				curr_p[1] = p[2];
			} else {
				curr_p[0] = p[2];
				curr_p[1] = p[0];
			}
			break;
		}

		/* Determine which ends starts the contour line segment
		 * passing through this triangle. For efficiency,
		 * pick the points so that next_p[0] < next_p[1].
		*/
		switch (type) {
		case 1: case 3:
			if (p[0] < p[2]) {
				next_p[0] = p[0];
				next_p[1] = p[2];
			} else {
				next_p[0] = p[2];
				next_p[1] = p[0];
			}
			break;
		case 2: case 6:
			if (p[0] < p[1]) {
				next_p[0] = p[0];
				next_p[1] = p[1];
			} else {
				next_p[0] = p[1];
				next_p[1] = p[0];
			}
			break;
		case 4: case 5:
			if (p[1] < p[2]) {
				next_p[0] = p[1];
				next_p[1] = p[2];
			} else {
				next_p[0] = p[2];
				next_p[1] = p[1];
			}
			break;
		}

		/* Increase edge buffer allocation, if necessary.
		 * we may add up to two edges at once.
		*/
		if (edges + (size_t)2 >= edges_max) {
			const size_t	 max = edges + (size_t)4096;
			struct edge	*temp;

			temp = realloc(edge, max * sizeof(struct edge));
			if (!temp) {
				if (edge)
					free(edge);
				errno = ENOMEM;
				return NULL;
			}

			edge = temp;
			edges_max = max;
		}

		/* Find the index of the edge where the bisecting line segment ends.
		 * If we create a new edge, calculate the point where the contour
		 * curve crosses the edge, too. (Mark the new edge active, "s=1", for later.)
		*/
		next_e = 0;
		while (next_e < edges && (edge[next_e].p[0] != next_p[0] || edge[next_e].p[1] != next_p[1]))
			next_e++;
		if (next_e >= edges) {
			/* Create new edge. */
			next_e = edges++;
			edge[next_e].p[0] = next_p[0];
			edge[next_e].p[1] = next_p[1];
			edge[next_e].e = next_e;
			edge[next_e].c = between(point[next_p[0]], point[next_p[1]], (level - value[next_p[0]]) / (value[next_p[1]] - value[next_p[0]]));
			edge[next_e].s = 1;
		}

		/* Find the index of the edge where the bisecting line segment starts.
		 * If we create a new edge, calculate the point where the contour
		 * curve crosses the edge, too. (Mark the new edge active, "s=1", for later.)
		*/
		curr_e = 0;
		while (curr_e < edges && (edge[curr_e].p[0] != curr_p[0] || edge[curr_e].p[1] != curr_p[1]))
			curr_e++;
		if (curr_e < edges) {
			/* Found; link. (This should always work.) */
			if (edge[curr_e].e == curr_e)
				edge[curr_e].e = next_e;
		} else {
			/* Create new edge. */
			curr_e = edges++;
			edge[curr_e].p[0] = curr_p[0];
			edge[curr_e].p[1] = curr_p[1];
			edge[curr_e].e = next_e;
			edge[curr_e].c = between(point[curr_p[0]], point[curr_p[1]], (level - value[curr_p[0]]) / (value[curr_p[1]] - value[curr_p[0]]));
			edge[curr_e].s = 1;
		}
	}

	/* Scan thorough all active edges. .s has values:
	 * 0=used, 1=unused/active, 2=tracing, 3=starting point.
	*/
	e = edges; while (e--) if (edge[e].s == 1) {

		/* Mark starting point. */
		edge[e].s = 3;

		/* Trace the curve; we need to know the length. */
		n = 2; /* Extra points to allocate */
		i = edge[e].e;
		while (edge[i].s == 1) {
			n++;
			edge[i].s = 2;
			i = edge[i].e;
		}

		/* Ignore points. */
		if (n < (size_t)3) {
			edge[e].s = 0;
			continue;
		}

		/* Allocate a new sub-path. */
		path = malloc(sizeof(struct path) + n * sizeof(struct point));
		if (!path) {
			/* Oops. Out of memory. */
			if (edge)
				free(edge);
			while (root) {
				path = root;
				root = root->next;
				free(path);
			}
			errno = ENOMEM;
			return NULL;
		}

		/* Prepend to the singly-linked list of sub-paths. */
		path->next = root;
		root = path;

		/* Retrace the curve, but this time copy and mark it too. */
		n = 0;
		i = e;
		edge[i].s = 2;
		while (edge[i].s == 2) {
			path->point[n] = edge[i].c;
			n++;
			edge[i].s = 0;
			i = edge[i].e;
		}

		/* If it is a closed curve, duplicate the starting point. */
		if (i == e) {
			path->point[n] = edge[e].c;
			n++;
		}
		path->points = n;
	}

	/* The edges are no longer needed. */
	if (edge)
		free(edge);

	/* Done. */
	errno = 0;
	return root;
}

#include <stdio.h>
#include <string.h>
#include <errno.h>
#include "isocurves.h"


int main(int argc, char *argv[])
{
	size_t	 points = 0;
	point_t	*point = NULL;
	double	*value = NULL;

	size_t	 max_tris = 0;
	size_t	 tris = 0;
	tri_t	*tri = NULL;

	path_t	*path, *curr;

	size_t	 offset, stride, i, n;

	double	 min, val;

	int	 arg, o;
	long	 nx, ny, x, y;
	char	 dummy;

	if (argc < 4) {
		fprintf(stderr, "\n"
		                "Usage: %s xpoints ypoints level:color ...\n"
		                "\n"
		                "This program will read xpoints*ypoints scalar values from\n"
		                "standard input, and output an SVG image, with specified\n"
		                "levels drawn as contour plots, and filled with the desired color.\n"
		                "Example level and color definitions:\n"
		                "\t3.141:#ffffff\tWhite contour at 3.141\n"
		                "\t-2.542#ff0000\tRed contour at -2.542\n"
		                "\t0.0012#00ff00\tGreen contour at 0.0012\n"
		                "\n"
		                "The resulting contours contain a lot of redundant information.\n"
		                "This program does NOT try to optimize the contour paths.\n"
		                "\n"
		              , argv[0]);
		return 1;
	}

	/* Parse grid size. */
	if (sscanf(argv[1], "%ld %c", &nx, &dummy) != 1 || nx < 2L) {
		fprintf(stderr, "%s: Invalid number of points along the x axis.\n", argv[1]);
		return 1;
	}
	if (sscanf(argv[2], "%ld %c", &ny, &dummy) != 1 || ny < 2L) {
		fprintf(stderr, "%s: Invalid number of points along the y axis.\n", argv[2]);
		return 1;
	}

	/* Reserve one grid point as a border. */
	stride = (size_t)(nx + 2L);
	offset = stride + (size_t)1;

	/* Allocate memory for the points, data, and triangles. */
	points = stride * (size_t)(ny + 2L);
	max_tris = (size_t)2 * points;
	point = malloc(points * sizeof(point_t));
	value = malloc(points * sizeof(double));
	tri = malloc(max_tris * sizeof(tri_t));
	if (!point || !value || !tri) {
		fprintf(stderr, "Cannot allocate enough memory.\n");
		return 1;
	}

	/* Initialize minimum to some ridiculous value. */
	min = 1.0e307;

	/* Read samples. */
	fprintf(stderr, "Reading samples for a %ldx%ld grid: ", nx, ny);
	fflush(stderr);
	for (y = 0L; y < ny; y++)
		for (x = 0L; x < nx; x++)
			if (fscanf(stdin, "%lf", &val) > 0) {
				if (min > val) min = val;
				value[offset + x + stride * y] = val;
				point[offset + x + stride * y].x = (double)x;
				point[offset + x + stride * y].y = (double)y;
			} else {
				fprintf(stderr, "Error reading sample %ldx%ld.\n", x + 1L, y + 1L);
				return 1;
			}

	/* Draw border with the minimum value, so we get full cycles. */
	for (y = 0L; y < ny; y++) {
		value[offset - 1 + stride * y] = min;
		point[offset - 1 + stride * y].x = -1.0;
		point[offset - 1 + stride * y].y = (double)y;

		value[offset + nx + stride * y] = min;
		point[offset + nx + stride * y].x = (double)nx;
		point[offset + nx + stride * y].y = (double)y;
	}
	for (x = -1L; x <= nx; x++) {
		value[offset + x - stride] = min;
		point[offset + x - stride].x = (double)x;
		point[offset + x - stride].y = -1.0;

		value[offset + x + stride * ny] = min;
		point[offset + x + stride * ny].x = (double)x;
		point[offset + x + stride * ny].y = (double)ny;
	}

	fprintf(stderr, "Done.\nGenerating triangular mesh: ");
	fflush(stderr);

	for (y = -1L; y < ny; y++)
		for (x = -1L; x < nx; x++) {
			const size_t	ul = offset + x + stride * y;
			const size_t	ur = ul + 1;
			const size_t	dl = ul + stride;
			const size_t	dr = dl + 1;
			if ((x ^ y) & 1) {
				tri[tris].p[0] = ul;
				tri[tris].p[1] = ur;
				tri[tris].p[2] = dl;
				tris++;
				tri[tris].p[0] = ur;
				tri[tris].p[1] = dr;
				tri[tris].p[2] = dl;
				tris++;				
			} else {
				tri[tris].p[0] = ul;
				tri[tris].p[1] = dr;
				tri[tris].p[2] = dl;
				tris++;
				tri[tris].p[0] = ul;
				tri[tris].p[1] = ur;
				tri[tris].p[2] = dr;
				tris++;
			}
		}

	fprintf(stderr, "%ld triangles.\n", (long)tris);
	fflush(stderr);

	printf("<?xml version=\"1.0\" standalone=\"no\" ?>\n"
	       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\" \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n"
	       "<svg width=\"20cm\" height=\"%.1fcm\" viewBox=\"0 0 %ld %ld\" xmlns=\"http://www.w3.org/2000/svg\" version=\"1.1\">\n",
	       20.0*ny/nx, nx, ny);

	for (arg = 3; arg < argc; arg++) {
		if (sscanf(argv[arg], "%lf%n", &val, &o) < 1) {
			fprintf(stderr, "%s: Cannot parse sample level. Skipping.\n", argv[arg]);
			fflush(stderr);
			continue;
		}
		while ((argv[arg][o] > 0 && argv[arg][o] < 33) ||
		       (argv[arg][o] == ':' || argv[arg][o] == '+' || argv[arg][o] == '/'))
			o++;
		if (!argv[arg][o]) {
			fprintf(stderr, "%s: No fill color specified. Skipping.\n", argv[arg]);
			fflush(stderr);
			continue;
		}

		fprintf(stderr, "Generating contour %.3g (%s): ", val, (const char *)&argv[arg][o]);
		fflush(stderr);

		path = isocurves(points, point, value, tris, tri, val);
		if (!path) {
			fprintf(stderr, "%s.\n", strerror(errno));
			fflush(stderr);
			continue;
		}

		printf("<path d=\"");

		while (path) {
			curr = path;
			path = path->next;

			n = curr->points - (size_t)1;
			printf("M %.5g,%.5g L", curr->point[0].x, curr->point[0].y);
			for (i = 1; i < n ; i++)
				printf(" %.5g,%.5g", curr->point[i].x, curr->point[i].y);

			if (curr->point[0].x == curr->point[n].x && curr->point[0].y == curr->point[n].y)
				printf(" Z");
			else
				printf(" %.5g,%.5g", curr->point[n].x, curr->point[n].y);

			if (path)
				putchar('\n');

			free(curr);
		}
		printf("\"       fill-rule=\"evenodd\" stroke=\"none\" fill=\"%s\" title=\"%.5g\" />\n", (const char *)&argv[arg][o], val);
		fflush(stdout);

		fprintf(stderr, "Done.\n");
		fflush(stderr);
	}

	printf("</svg>\n");
	return 0;
}

gcc -Wall -pedantic -std=c99 -O3 -o test main.c isocurves.c

N=50

awk -v n=$N 'BEGIN {
	u=8.0/n
	for (y=0;y<n;y++)
		for (x=0;x<n;x++)
			printf("%.5g\n", cos((x*u-4.0)**2)*cos((y*u-4.0)**2))
}' | ./test $N $N '-0.5#900' '0.0#f00' '0.5#ff0' '0.75#fff' > test.svg