/* 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);
		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;

#ifdef TODO /* untested */
	/* 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 ((y3 - y1) * (x3 - x2) == (y3 - y2) * (x2 - x1)) {
		if (x1 <= x2 && x2 <= x3 && y1 <= y2 && y2 <= y3)
			rtgrpblib_draw_linear_bezier(raster, x1, y1, x3, y3);
		else if (x1 <= x3 && x3 <= x2 && y1 <= y3 && y3 <= y2)
			rtgrpblib_draw_linear_bezier(raster, x1, y1, x2, y2);
		else
			rtgrpblib_draw_linear_bezier(raster, x2, y2, x3, y3);
		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_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


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


#endif