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


static size_t
solve_cubic_bezier(double *restrict ts, double a, double b, double c, double d, double z)
{
	/*
	 * B(t) = t(q - p) + p
	 *   where
	 *     p = t(v - u) + u
	 *     q = t(w - v) + v
	 *       where
	 *          u = t(b - a) + a
	 *          v = t(c - b) + b
	 *          w = t(d - c) + c
	 * 
	 * B(t) = t³(d - 3c + 3b - a) + t²(3c - 6b + 3a) + t(3b - 3a) + a
	 * 
	 * z = t³(d - 3c + 3b - a) + t²(3c - 6b + 3a) + t(3b - 3a) + a
	 * 
	 * t³(d - 3c + 3b - a) + t²(3c - 6b + 3a) + t(3b - 3a) + a - z = 0
	 *
	 * A := d - 3c + 3b - a
	 * B := 3c - 6b + 3a
	 * C := 3b - 3a
	 * D := a - z
	 * 
	 * t³A + t²B + tC + D = 0
	 */
	double c3 = 3.0 * c;
	double b3 = 3.0 * b;
	double a3 = 3.0 * a;
	double A = (d - c3) + (b3 - a);
	double B = (c3 - b3) - (b3 - a3);
	double C = b3 - a3;
	double D = a - z;

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


static double
evaluate_cubic_bezier(double t, double a, double b, double c, double d)
{
	double u = fma(t, b - a, a);
	double v = fma(t, c - b, b);
	double w = fma(t, d - c, c);
	double p = fma(t, v - u, u);
	double q = fma(t, w - v, v);
	return fma(t, q - p, p);
}


static double
evaluate_cubic_bezier_derivative(double t, double a, double b, double c, double d)
{
	double u = fma(t, b - a, a);
	double v = fma(t, c - b, b);
	double w = fma(t, d - c, c);
	return fma(t, (w - v) + (u - v), v - u);
}


static void
draw_bounded_cubic_bezier(RASTER *restrict raster, double x1, double y1, double x2, double y2,
                          double x3, double y3, double x4, double y4, 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;

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

	for (; t < t2; x0 = x, y0 = y) {
		dx = evaluate_cubic_bezier_derivative(t, x1, x2, x3, x4);
		dy = evaluate_cubic_bezier_derivative(t, y1, y2, y3, y4);

		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_cubic_bezier(t, x1, x2, x3, x4);
		y = evaluate_cubic_bezier(t, y1, y2, y3, y4);
		x = fmin(fmax(0, x), w);
		y = fmin(fmax(0, y), h);

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


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

#ifdef TODO /* untested */
	/* Can we downgrade the curve to a quadratic Bézier curve?
	 * 
	 * B(t) = (1-t)²Q₁ + 2(1-t)tQ₂ + t²Q₃
	 * B(t) = (1-t)B(t) + tB(t) =
	 *      = (1-t)((1-t)²Q₁ + 2(1-t)tQ₂ + t²Q₃) + t((1-t)²Q₁ + 2(1-t)tQ₂ + t²Q₃)
	 *      = (1-t)³Q₁ + 2(1-t)²tQ₂ + (1-t)t²Q₃ + (1-t)²tQ₁ + 2(1-t)t²Q₂ + t³Q₃
	 *      = (1-t)³Q₁ + 2(1-t)²tQ₂ + (1-t)²tQ₁ + 2(1-t)t²Q₂ + (1-t)t²Q₃ + t³Q₃
	 *      = (1-t)³Q₁ + (1-t)²t(Q₁ + 2Q₂) + (1-t)t²(2Q₂ + Q₃) + t³Q₃
	 *      = (1-t)³Q₁ + 3(1-t)²t(Q₁ + 2Q₂)/3 + 3(1-t)t²(2Q₂ + Q₃)/3 + t³Q₃
	 * 
	 * B(t) = (1-t)³C₁ + 3(1-t)²tC₂ + 3(1-t)t²C₃ + t³C₄
	 * 
	 * C₁ = Q₁, C₄ = Q₃
	 * 
	 * B(t) = (1-t)³Q₁ + 3(1-t)²t(Q₁ + 2Q₂)/3 + 3(1-t)t²(2Q₂ + Q₃)/3 + t³Q₃
	 * B(t) = (1-t)³C₁ + 3(1-t)²tC₂ + 3(1-t)t²C₃ + t³C₄
	 * 
	 * (1-t)³Q₁ + 3(1-t)²t(Q₁ + 2Q₂)/3 + 3(1-t)t²(2Q₂ + Q₃)/3 + t³Q₃ = (1-t)³C₁ + 3(1-t)²tC₂ + 3(1-t)t²C₃ + t³C₄w
	 * 3(1-t)²t(Q₁ + 2Q₂)/3 + 3(1-t)t²(2Q₂ + Q₃)/3 = 3(1-t)²tC₂ + 3(1-t)t²C₃
	 * (1-t)²t(Q₁ + 2Q₂)/3 + (1-t)t²(2Q₂ + Q₃)/3 = (1-t)²tC₂ + (1-t)t²C₃
	 * (1-t)t(Q₁ + 2Q₂)/3 + t²(2Q₂ + Q₃)/3 = (1-t)tC₂ + t²C₃
	 * (1-t)(Q₁ + 2Q₂)/3 + t(2Q₂ + Q₃)/3 = (1-t)C₂ + tC₃
	 * s(Q₁ + 2Q₂)/3 + t(2Q₂ + Q₃)/3 = sC₂ + tC₃
	 * (Q₁ + 2Q₂)/3 = C₂, (2Q₂ + Q₃)/3 = C₃
	 * Q₁ + 2Q₂ = 3C₂, 2Q₂ + Q₃ = 3C₃
	 * C₁ + 2Q₂ = 3C₂, 2Q₂ + C₄ = 3C₃
	 * 2Q₂ = 3C₂ - C₁, 2Q₂ = 3C₃ - C₄
	 * 2Q₂ = 3C₂ - C₁ = 3C₃ - C₄
	 */
	if (x2 + x2 + x2 - x1 == x3 + x3 + x3 - x4 && y2 + y2 + y2 - y1 == y3 + y3 + y3 - y4) {
		rtgrpblib_draw_quadratic_bezier(raster, x1, y1, x2 + (x2 - x1) / 2.0, y2 + (y2 - y1) / 2.0, x4, y4);
		return;
	}
#endif

	/* 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_cubic_bezier(&ts[nts], y1, y2, y3, y4, 0);
	nts += solve_cubic_bezier(&ts[nts], y1, y2, y3, y4, (double)raster->height);
	nts += solve_cubic_bezier(&ts[nts], x1, x2, x3, x4, 0);
	nts += solve_cubic_bezier(&ts[nts], x1, x2, x3, x4, (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_cubic_bezier(t, y1, y2, y3, y4);
		if (y < 0 || y > (double)raster->height)
			continue;

		/* If the segment is inside the raster, draw it, */
		x = evaluate_cubic_bezier(t, x1, x2, x3, x4);
		if (0 <= x && x <= (double)raster->width) {
			draw_bounded_cubic_bezier(raster, x1, y1, x2, y2, x3, y3, x4, y4, 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_cubic_bezier(t1, y1, y2, y3, y4);
		yt2 = evaluate_cubic_bezier(t2, y1, y2, y3, y4);
		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


int
main(void)
{
	return 0; /* TODO add test */
}


#endif