/*
 * area.c - Area fill
 *
 * Written 2012 by Werner Almesberger
 * Copyright 2012 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.
 */

/*
 * We use the following requirement to simplify toolpath generation: the
 * outlines must be designed such that the tool can pass along all the
 * outlines without cutting into anything it's not supposed to.
 */

#include <stdio.h>
#include <stddef.h>
#include <math.h>
#include <assert.h>

#include "util.h"
#include "path.h"
#include "area.h"


#define	EPSILON	0.0001


static int bbox(const struct path *path,
    double *xa, double *ya, double *xb, double *yb)
{
	const struct point *p = path->first;

	if (!p)
		return 0;
	*xa = *xb = p->x;
	*ya = *yb = p->y;
	while (p) {
		if (p->x < *xa)
			*xa = p->x;
		if (p->x > *xb)
			*xb = p->x;
		if (p->y < *ya)
			*ya = p->y;
		if (p->y > *yb)
			*yb = p->y;
		p = p->next;
	}
	return 1;
}


/*
 * @@@ this is a bit too simple. E.g., it would report A as being inside B
 * in this case:
 *
 *          +---+
 *   +---+  |   |
 *   | A |  |   |
 *   +---+  |   |
 *          | B |
 * +--------+   |
 * |            |
 * +------------+
 */

static int is_inside(const struct path *a, const struct path *b)
{
    	double xa, ya, xb, yb;
	const struct point *p;

	if (!bbox(b, &xa, &ya, &xb, &yb))
		return 0;
	for (p = a->first; p; p = p->next)
		if (p->x < xa || p->x > xb ||
		    p->y < ya || p->y > yb)
			return 0;
	return 1;
}


/*
 * Solve
 *
 * ax+na*bx = cx+nb*dx
 * ay+na*by = cy+nb*dy
 *
 * which is
 *
 * na*bx + nb*-dx = cx-ax
 * na*by + nb*-dy = cy-ay
 *
 * which we the solve with Cramer's rule:
 * http://en.wikipedia.org/wiki/Cramer's_rule
 */

static int intersect(double ax, double ay, double bx, double by,
    double cx, double cy, double dx, double dy, double *na, double *nb)
{
	double det;

	det = dx*by-bx*dy;
	if (fabs(det) < EPSILON)
		return 0;

	*na = (dx*(cy-ay)-dy*(cx-ax))/det;
	*nb = (bx*(cy-ay)-by*(cx-ax))/det;
	return 1;
}


/*
 * Solve
 *
 * (ax+n*bx-cx)^2+(ay+n*by-cy)^2 = r^2 for n
 *
 * http://en.wikipedia.org/wiki/Quadratic_equation
 */

static int touch(double ax, double ay, double bx, double by,
    double cx, double cy, double r, double *n)
{
	double dx = cx-ax;
	double dy = cy-ay;
	double a = bx*bx+by*by;
	double b = -2*bx*dx-2*by*dy;
	double c = dx*dx+dy*dy-r*r;
	double d, n0, n1;

	d = b*b-4*a*c;
	if (d < 0)
		return 0;
	d = sqrt(d);
	n0 = (-b-d)/2/a;
	n1 = (-b+d)/2/a;

	if (n0 > 0) {
		*n = n0;
		return 1;
	}
	if (n1 > 0) {
		*n = n1;
		return 1;
	}
	return 0;
}


static int hit_segment(double fx, double fy, double tx, double ty,
    const struct point *a, const struct point *b, double r, int enter,
    double *n)
{
	double dx, dy, d;
	double px, py;
	double na, nb;

printf("      seg (%g,%g)+(%g,%g) -> (%g,%g)-(%g,%g)\n",
    fx, fy, tx, ty, a->x, a->y, b->x, b->y);
	tx -= fx;
	ty -= fy;

	dx = b->x-a->x;
	dy = b->y-a->y;

	if (enter && dy < 0)
		return 0;
	if (!enter && dy > 0)
		return 0;

	d = hypot(dx, dy);

	px = a->x-dy/d*r;
	py = a->y+dx/d*r;

	if (!intersect(fx, fy, tx, ty, px, py, dx, dy, &na, &nb))
		return 0;
printf("\tna %g (%g) nb %g (%g)\n", na, fx+tx*na, nb, fx+tx*nb);
	if (nb <= 0) {
		if (!touch(fx, fy, tx, ty,  a->x, a->y, r, &na))
			return 0;
	}
	if (nb >= 1) {
		if (!touch(fx, fy, tx, ty,  b->x, b->y, r, &na))
			return 0;
	}
	if (na <= 0 || na >= 1)
		return 0;
	*n = na;
	return 1;
}


static int hit_path(double fx, double fy, double tx, double ty,
    const struct path *path, int inside, int enter, double r, double *x)
{
	const struct point *p;
	int left;
	double nx, tmp;
	int found = 0;

	left = path_tool_is_left(path);
	if (inside)
		left = !left;
	for (p = path->first; p != path->last; p = p->next) {
		if (hit_segment(fx, fy, tx, ty,
		    left ? p : p->next, left ? p->next : p, r, enter, &tmp)) {
			if (!found || nx > tmp)
				nx = tmp;
			found = 1;
		}
	}
	if (found)
		*x = fx+nx*(tx-fx);
	return found;
}


static const struct path **subordinates(const struct path *paths,
    const struct path *path)
{
	const struct path **sub, **w, **a, **b;;
	const struct path *p;
	int n = 0;

	for (p = paths; p; p = p->next)
		n++;
	sub = alloc_size(sizeof(struct path *)*n);
	w = sub;
	for (p = paths; p; p = p->next)
		if (p != path && is_inside(p, path) && !is_inside(path, p))
			*w++ = p;
	*w = NULL;
	for (a = sub; a != w; a++)
		for (b = sub; b != w; b++)
			if (a != b && is_inside(*a, *b)) {
				*a = *w--;
				*w = NULL;
				a--;
				break;
			}
	return sub;
}


static void do_line(const struct path *path, const struct path **sub,
    double xa, double xb, double y, double r_tool, double overlap,
    struct path **res)
{
	const struct path *last = path;
	const struct path **s;
	struct path *new;
	double x, next;

printf("  y=%g\n", y);
	if (!hit_path(xa-3*r_tool, y, xb, y, last, 1, 0, r_tool, &x))
		return;
	while (1) {
printf("    x=%g\n", x);
		next = xb;
		last = NULL;
		if (hit_path(x, y, xb, y, path, 1, 1, r_tool, &next))
			last = path;
		for (s = sub; *s; s++)
			if (hit_path(x, y, next, y, *s, 0, 1, r_tool, &next))
				last = *s;
		if (next-x > 2*r_tool-2*overlap) {
			new = path_new(r_tool, "");
			path_add(new, x, y, path->first->z);
			path_add(new, next, y, path->first->z);
			new->next = *res;
			*res = new;
		}
		if (!last)
			return;
		if (!hit_path(next+EPSILON, y, xb, y, last, last == path, 0,
		    r_tool, &x))
			return;
	}
}


static void add_outline(const struct path *path, int inside, struct path **res)
{
	struct path *new;
	int left;

	left = path_tool_is_left(path);
	new = path_offset(path, inside ? !left : left, 0);
	new->next = *res;
	*res = new;
}


static void fill_path(const struct path *paths, const struct path *path,
    double r_tool, double overlap, struct path **res)
{
	const struct path **sub, **s;
	const struct path **sub2, **s2;
    	double xa, ya, xb, yb;
	int n, i;

	if (!bbox(path, &xa, &ya, &xb, &yb))
		return;
	sub = subordinates(paths, path);
	xa += r_tool;
	ya += 3*r_tool-overlap;
	xb -= r_tool;
	yb -= 3*r_tool-overlap;
	n = ceil((yb-ya)/(2*r_tool-overlap));
printf("x[%g:%g] y[%g:%g] n=%d\n", xa, xb, ya, yb, n);
	for (i = 0; i <= n; i++)
		do_line(path, sub, xa, xb, ya+(yb-ya)*((double) i/n),
		    r_tool, overlap, res);
	for (s = sub; *s; s++) {
		sub2 = subordinates(paths, *s);
		for (s2 = sub2; *s2; s2++)
			fill_path(paths, *s2, r_tool, overlap, res);
		free(sub2);
		add_outline(*s, 0, res);
	}
	free(sub);
	add_outline(path, 1, res);
}


struct path *area(const struct path *path, double r_tool, double overlap)
{
	struct path *res = NULL;

	if (!path)
		return NULL;
	fill_path(path, path_find_leftmost(path), r_tool, overlap, &res);
	return res;
}