/* See LICENSE file for copyright and license details. */
#include "common.h"
#ifndef TEST


static size_t
solve_quadratic_bezier(double *restrict ts, double a, double b, double c, double z)
{
	/*
	 * B(t) = t(q - p) + p
	 *   where
	 *     p = t(b - a) + a
	 *     q = t(c - b) + b
	 * 
	 * B(t) = t²(c - 2b + a) + t(2b) + a
	 *
	 * z = t²(c - 2b + a) + t(2b) + a
	 *
	 * t²(c - 2b + a) + t(2b) + a - z = 0
	 * 
	 * A := c - 2b + a
	 * B := 2b
	 * C := a - z
	 * 
	 * t²A + tB + C = 0
	 */

	double A = (c - a) + (b + b);
	double B = b + b;
	double C = a - z;

	if (A)
		return solve_quadratic(ts, A, B, C);
	else if (B)
		return solve_linear(ts, B, C);
	else
		return 0;
}


static double
evaluate_quadratic_bezier(double t, double a, double b, double c)
{
	double p = fma(t, b - a, a);
	double q = fma(t, c - b, b);
	return fma(t, q - p, p);
}


static void
draw_bounded_quadratic_bezier(RASTER *restrict raster, double x1, double y1, double x2, double y2,
                              double x3, double y3, double t1, double t2)
{
	double t = t1, x0, y0, dx, dy, x, y;
	double w = (double)raster->width - 0.0001;
	double h = (double)raster->height - 0.0001;

	double kx = (x3 - x2) + (x1 - x2), mx = x2 - x1;
	double ky = (y3 - y2) + (y1 - y2), my = y2 - y1;

	x0 = evaluate_quadratic_bezier(t, x1, x2, x3);
	y0 = evaluate_quadratic_bezier(t, y1, y2, y3);
	x0 = fmin(fmax(0, x0), w);
	y0 = fmin(fmax(0, y0), h);

	for (; t < t2; x0 = x, y0 = y) {
		dx = fma(t, kx, mx);
		dy = fma(t, ky, my);

		t += raster->draftness / hypot(dx, dy);
		/* TODO deal with unchanged t and oversized steps (e.g. locally zero derivative) */
		t = fmin(t, t2);
		x = evaluate_quadratic_bezier(t, x1, x2, x3);
		y = evaluate_quadratic_bezier(t, y1, y2, y3);

		x = fmin(fmax(0, x), w);
		y = fmin(fmax(0, y), h);

		draw_bounded_line(raster, x0, y0, x, y);
	}
}


void
rtgrpblib_draw_quadratic_bezier(RASTER *restrict raster, double x1, double y1,
                                double x2, double y2, double x3, double y3)
{
	double x, y, yt1, yt2, dy;
	double ts[2+2*4], t, t1, t2;
	size_t nts = 0;
	size_t i;

	/* Can we downgrade the curve to a linear Bézier curve?
	 *
	 *  (y2 - y1)/(x2 - x1) = (y3 - y2)/(x3 - x2) = (y3 - y1)/(x3 - x1)
	 *  (y2 - y1)/(x2 - x1) = (y3 - y2)/(x3 - x2)
	 *  (y2 - y1)(x3 - x2) = (y3 - y2)(x2 - x1)
	 */
	if ((y2 - y1) * (x3 - x2) == (y3 - y2) * (x2 - x1)) {
		/* Drawing (x1, y1)--(x2, y2) plus (x2, y2)--(x3, y3)
		 * is equivalent to drawing (x1, y1)--(x3, y3) because
		 * the result is signed and if (x2, y2) is outside the
		 * (x1, y1)--(x3, y3) line (x2, y2)--(x3, y3) erases
		 * the part that (x1, y1)--(x2, y2) extends to
		 * (x1, y1)--(x3, y3) */
		rtgrpblib_draw_linear_bezier(raster, x1, y1, x3, y3);
		return;
	}

	/* Beginning and end of curve */
	ts[nts++] = 0.0;
	ts[nts++] = 1.0;
	/* Get points where the end intersects one of the edges of the raster */
	nts += solve_quadratic_bezier(&ts[nts], y1, y2, y3, (double)raster->height);
	nts += solve_quadratic_bezier(&ts[nts], x1, x2, x3, (double)raster->width);
	/* Sort all points of interest */
	qsort(ts, nts, sizeof(*ts), doublepcmp);

	for (i = 0; i < nts - 1; i++) {
		t1 = ts[i];
		t2 = ts[i + 1];
		if (t1 == t2)
			continue;
		t = 0.5 * (t1 + t2);

		/* Remove any segments above or below the raster */
		y = evaluate_quadratic_bezier(t, y1, y2, y3);
		if (y < 0 || y > (double)raster->height)
			continue;

		/* If the segment is inside the raster, draw it, */
		x = evaluate_quadratic_bezier(t, x1, x2, x3);
		if (0 <= x && x <= (double)raster->width) {
			draw_bounded_quadratic_bezier(raster, x1, y1, x2, y2, x3, y3, t1, t2);
			continue;
		}

		/* otherwise draw it flattened to the edge; however,
		 * we will reduce it to a from the start to the end
		 * and draw its project on the edge, which works
		 * because if the curve bends it will cancel it self
		 * out except within this project; this project will
		 * also preserv the direction of part of the curve
		 * that is not cancelled out */
		yt1 = evaluate_quadratic_bezier(t1, y1, y2, y3);
		yt2 = evaluate_quadratic_bezier(t2, y1, y2, y3);
		dy = yt2 - yt1;
		if (x < 0)
			draw_vertical_line(raster, 0, yt1, yt2, SIGNUM(dy));
		else
			draw_vertical_line_opposite_only(raster, raster->width - 1, yt1, yt2, SIGNUM(dy));
	}
}


#else


static RASTER *raster;
static RASTER *refraster;


static double
frand(double a, double b)
{
	int i = rand();
	double f = (double)i / (double)RAND_MAX;
	double max = fmax(a, b);
	double min = fmin(a, b);
	return fma(f, max - min, min);
}


static void
check_line(double x1, double y1, double x2, double y2, double x3, double y3, int exact)
{
	size_t i;

	alarm(5);
	rtgrpblib_draw_quadratic_bezier(raster, x1, y1, x2, y2, x3, y3);

	alarm(5);
	rtgrpblib_draw_linear_bezier(refraster, x1, y1, x3, y3);

	alarm(0);

	if (exact) {
		ASSERT(!memcmp(raster->cells, refraster->cells, 100 * sizeof(*raster->cells)));
	} else {
		for (i = 0; i < 100; i++) {
			if (!eq(raster->cells[i].cell_coverage, refraster->cells[i].cell_coverage))
				fprintf(stderr, "c%02zu: %la ref: %la; diff: %lg=%la (tol: %la)\n", i, raster->cells[i].cell_coverage, refraster->cells[i].cell_coverage, refraster->cells[i].cell_coverage - raster->cells[i].cell_coverage, refraster->cells[i].cell_coverage - raster->cells[i].cell_coverage, TOLERANCE);
			if (!eq(raster->cells[i].opposite_coverage, refraster->cells[i].opposite_coverage))
				fprintf(stderr, "o%02zu: %la ref: %la; diff: %lg=%la (tol: %la)\n", i, raster->cells[i].opposite_coverage, refraster->cells[i].opposite_coverage, raster->cells[i].opposite_coverage - refraster->cells[i].opposite_coverage, raster->cells[i].opposite_coverage - refraster->cells[i].opposite_coverage, TOLERANCE);
		}
		for (i = 0; i < 100; i++) {
			if (!eq(raster->cells[i].cell_coverage, refraster->cells[i].cell_coverage) ||
			    !eq(raster->cells[i].opposite_coverage, refraster->cells[i].opposite_coverage)) {
				print_raster(refraster);
				print_raster(raster);
			}
			ASSERT(eq(raster->cells[i].cell_coverage, refraster->cells[i].cell_coverage));
			ASSERT(eq(raster->cells[i].opposite_coverage, refraster->cells[i].opposite_coverage));
		}
	}
	rtgrpblib_reset_raster(raster, 10, 10);
	rtgrpblib_reset_raster(refraster, 10, 10);
}


static void
check_simple_lines(void)
{
	check_line(1, 1, 2, 2, 3, 3, 1);
	check_line(1, 1, 3, 3, 2, 2, 1);
	check_line(2, 2, 1, 1, 3, 3, 1);
	check_line(2, 2, 3, 3, 1, 1, 1);
	check_line(3, 3, 1, 1, 2, 2, 1);
	check_line(3, 3, 2, 2, 1, 1, 1);
	check_line(1, 1, 3, 3, 3, 3, 1);
	check_line(1, 1, 1, 1, 3, 3, 1);

	check_line(1, 4, 2, 4, 3, 4, 1);
	check_line(1, 4, 3, 4, 2, 4, 1);
	check_line(2, 4, 1, 4, 3, 4, 1);
	check_line(2, 4, 3, 4, 1, 4, 1);
	check_line(3, 4, 1, 4, 2, 4, 1);
	check_line(3, 4, 2, 4, 1, 4, 1);
	check_line(1, 4, 3, 4, 3, 4, 1);
	check_line(1, 4, 1, 4, 3, 4, 1);

	check_line(4, 1, 4, 2, 4, 3, 1);
	check_line(4, 1, 4, 3, 4, 2, 1);
	check_line(4, 2, 4, 1, 4, 3, 1);
	check_line(4, 2, 4, 3, 4, 1, 1);
	check_line(4, 3, 4, 1, 4, 2, 1);
	check_line(4, 3, 4, 2, 4, 1, 1);
	check_line(4, 1, 4, 3, 4, 3, 1);
	check_line(4, 1, 4, 1, 4, 3, 1);
}


static void
check_random_bounded_lines(int exact)
{
	double x1, x2, x3;
	double y1, y2, y3;
	double mid;
	size_t i;

	for (i = 0; i < 100000UL; i++) {
		x1 = floor(frand(0, 11));
		y1 = floor(frand(0, 11));
		x3 = floor(frand(0, 11));
		y3 = floor(frand(0, 11));
		x1 -= (double)(x1 > 10);
		x2 -= (double)(x2 > 10);
		y1 -= (double)(y1 > 10);
		y2 -= (double)(y2 > 10);
		mid = frand(0, 1);

		if (x3 - x1) {
			x2 = floor(mid * (x3 - x1) + x1);
			y2 = y1 + (x2 - x1) * (y3 - y1) / (x3 - x1);
		} else if (y3 - y1) {
			x2 = x1;
			y2 = floor(mid * (y3 - y1) + y1);
		} else {
			x2 = x1;
			y2 = y1;
		}

		if (!exact)
			check_line(x1, y1, x2, y2, x3, y3, x2 == floor(x2) && y2 == floor(y2));
		else if (x2 == floor(x2) && y2 == floor(y2))
			check_line(x1, y1, x2, y2, x3, y3, 1);
	}
}


static void
check_random_unbounded_lines(int exact)
{
	double x1, x2, x3;
	double y1, y2, y3;
	double mid;
	size_t i;

	for (i = 0; i < 100000UL; i++) {
		x1 = floor(frand(-5, 15));
		y1 = floor(frand(-5, 15));
		x3 = floor(frand(-5, 15));
		y3 = floor(frand(-5, 15));
		mid = frand(-2, 2);

		if (x3 - x1) {
			x2 = floor(mid * (x3 - x1) + x1);
			y2 = y1 + (x2 - x1) * (y3 - y1) / (x3 - x1);
		} else if (y3 - y1) {
			x2 = x1;
			y2 = floor(mid * (y3 - y1) + y1);
		} else {
			x2 = x1;
			y2 = y1;
		}

		if (!exact)
			check_line(x1, y1, x2, y2, x3, y3, x2 == floor(x2) && y2 == floor(y2));
		else if (x2 == floor(x2) && y2 == floor(y2))
			check_line(x1, y1, x2, y2, x3, y3, 1);
	}
}


int
main(void)
{
	raster = rtgrpblib_create_raster(10, 10);
	ASSERT(raster);
	refraster = rtgrpblib_create_raster(10, 10);
	ASSERT(refraster);
	srand((unsigned)(uintptr_t)raster);

	check_simple_lines();
	check_random_bounded_lines(1);
	check_random_unbounded_lines(1);
	//TODO check_random_bounded_lines(0);
	//TODO check_random_unbounded_lines(0);
	raster->draftness = -1;
	refraster->draftness = -1;
	check_simple_lines();
	check_random_bounded_lines(1);
	check_random_unbounded_lines(1);
#ifdef TODO_NOT_IMPLEMENTED
	check_random_bounded_lines(0);
	check_random_unbounded_lines(0);
#endif
	raster->draftness = 0.5;
	refraster->draftness = 0.5;

	/* TODO add tests */

	free(raster);
	free(refraster);
	return 0;
}


#endif