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

/*
 * The implementation in this fill is not supposed to be fast,
 * it is supposed to be absolutely correct. The implementation
 * shall be as clear and obvious as possible. It is important
 * that everything written to cells is as precise as possible
 * and that the function behavious correctly when outside the
 * raster. This implementation shall be usable to generate an
 * as perfect as possible rasterisation and shall be usable to
 * test the correctness and accuracy of the fast implementation.
 */


void
draw_linear_bezier_reference(RASTER *restrict raster, double x1, double y1, double x2, double y2)
{
	double dx = x2 - x1;
	double dy = y2 - y1;
	int xdir = SIGNUM(dx);
	int ydir = SIGNUM(dy);
	double kx = dx / dy;
	double ky = dy / dx;
	double prevX, prevY, x, y, cx, cy, h2x, h2y;
	double xfloor, yfloor;
	ssize_t cellX, cellY;
	size_t cell;

	prevX = x1;
	prevY = y1;

	for (; prevX != x2 || prevY != y2; prevX = x, prevY = y) {
		/* Find next vertical and next hozitonal grid line */
		x = xdir < 0 ? floor(prevX) : xdir > 0 ? ceil(prevX) : prevX;
		y = ydir < 0 ? floor(prevY) : ydir > 0 ? ceil(prevY) : prevY;
		x = iszeroish(x - prevX) ? x + xdir : x;
		y = iszeroish(y - prevY) ? y + ydir : y;

		/* Limit to the extend of the drawn line segment */
		if ((xdir < 0 && x < x2) || (xdir > 0 && x > x2))
			x = x2;
		if ((ydir < 0 && y < y2) || (ydir > 0 && y > y2))
			y = y2;

		/* Find next cell edge */
		if (xdir && ydir) {
			cy = y1 + (x - x1) * ky;
			cx = x1 + (y - y1) * kx;
			h2x = (cy - prevY) * (cy - prevY) + (x - prevX) * (x - prevX);
			h2y = (cx - prevX) * (cx - prevX) + (y - prevY) * (y - prevY);
			if (h2x < h2y)
				y = cy;
			else
				x = cx;

			/* Make sure loop will exit when we reach the end */
			if ((xdir < 0 && x <= x2) || (xdir > 0 && x >= x2) ||
			    (ydir < 0 && y <= y2) || (ydir > 0 && y >= y2)) {
				x = x2;
				y = y2;
			}
		}

		/* Select cell to draw in */
		xfloor = xdir >= 0 ? floor(prevX) : floor(x);
		yfloor = ydir >= 0 ? floor(prevY) : floor(y);
		cellX = (ssize_t)xfloor;
		cellY = (ssize_t)yfloor;

		/* Do not draw if above or below the raster */
		if (cellY < 0 || (size_t)cellY >= raster->height)
			continue;

		if (cellX < 0) {
			/* Draw on first column in raster of outside of raster on the left side,
			 * with full horizontal coverage, `.opposite_coverage` one the actual
			 * cell (outside of the raster) will contribute to the first column's
			 * `.cell_coverage` */
			cell = (size_t)cellY * raster->width;
			raster->cells[cell].opposite_coverage += y - prevY;
			raster->cells[cell].cell_coverage += y - prevY;

		} else if ((size_t)cellX >= raster->width) {
			/* If outside of the raster, on the right side, add to the last
			 * columns contributions to the next column, so that the ink level
			 * balances out to zero */
			cell = (size_t)cellY * raster->width + (raster->width - 1);
			raster->cells[cell].opposite_coverage += y - prevY;

		} else {
			/* If inside the raster, just draw the line */
			cell = (size_t)cellY * raster->width + (size_t)cellX;
			raster->cells[cell].opposite_coverage += y - prevY;
			raster->cells[cell].cell_coverage += (y - prevY) * fabs(xfloor + 1.0 - (x + prevX) / 2.0);
		}
	}
}


#else


#define CHECK_ZEROED()\
	do {\
		size_t i__;\
		for (i__ = 0; i__ < 100; i__++) {\
			ASSERT(!raster->cells[i__].cell_coverage);\
			ASSERT(!raster->cells[i__].opposite_coverage);\
		}\
	} while (0)

#define CHECK_CELL(X, Y, AREA, SHADOW)\
	do {\
		size_t i__ = (size_t)(Y) * raster->width + (size_t)(X);\
		ASSERT(raster->cells[i__].cell_coverage == (AREA));\
		ASSERT(raster->cells[i__].opposite_coverage == (SHADOW));\
		raster->cells[i__].cell_coverage = 0;\
		raster->cells[i__].opposite_coverage = 0;\
	} while (0)

#define SUB_CELL(X, Y, AREA, SHADOW)\
	do {\
		size_t i__ = (size_t)(Y) * raster->width + (size_t)(X);\
		raster->cells[i__].cell_coverage -= (AREA);\
		raster->cells[i__].opposite_coverage -= (SHADOW);\
	} while (0)


static RASTER *raster;


static void
draw(double x1, double y1, double x2, double y2)
{
	alarm(5);
	draw_linear_bezier_reference(raster, x1, y1, x2, y2);
	alarm(0);
}


static void
check_horizontal(void)
{
	ssize_t y, x1, x2;
	double fy, fx1, fx2;

	for (y = -2; y < 12; y++) {
		for (x1 = -2; x1 < 12; x1++) {
			for (x2 = -2; x2 < 12; x2++) {
				fy = (double)y;
				fx1 = (double)x1;
				fx2 = (double)x2;

				draw(fx1, fy, fx2, fy);
				CHECK_ZEROED();

				draw(fx1 + 0.25, fy, fx2 + 0.25, fy);
				CHECK_ZEROED();

				draw(fx1, fy + 0.25, fx2, fy + 0.25);
				CHECK_ZEROED();

				draw(fx1 + 0.25, fy + 0.25, fx2 + 0.25, fy + 0.25);
				CHECK_ZEROED();
			}
		}
	}
}


static void
check_vertical(void)
{
	ssize_t x, y1, y2, i, ydir;
	double yhit;
	double fx, fy1, fy2;

	for (x = -2; x < 12; x++) {
		for (y1 = -2; y1 < 12; y1++) {
			for (y2 = -2; y2 < 12; y2++) {
				ydir = SIGNUM(y2 - y1);
				fx = (double)x;
				fy1 = (double)y1;
				fy2 = (double)y2;

				draw(fx, fy1, fx, fy2);
				for (i = 0; i < 10; i++) {
					if (i < y1 && i < y2)
						continue;
					if (i >= y1 && i >= y2)
						continue;
					CHECK_CELL(x < 0 ? 0 : x < 10 ? x : 9,
					           i,
						   x < 0 ? (double)ydir : x < 10 ? (double)ydir : 0,
					           (double)ydir);
				}
				CHECK_ZEROED();

				draw(fx + 0.25, fy1, fx + 0.25, fy2);
				for (i = 0; i < 10; i++) {
					if (i < y1 && i < y2)
						continue;
					if (i >= y1 && i >= y2)
						continue;
					CHECK_CELL(x < 0 ? 0 : x < 10 ? x : 9,
					           i,
						   x < 0 ? (double)ydir : x < 10 ? (double)ydir * 0.75 : 0,
					           (double)ydir);
				}
				CHECK_ZEROED();

				draw(fx, fy1 + 0.25, fx, fy2 + 0.25);
				for (i = 0; i < 10; i++) {
					if (i < y1 + (y1 > y2) && i < y2 + (y2 > y1))
						continue;
					if (i >= y1 + (y1 > y2) && i >= y2 + (y2 > y1))
						continue;
					yhit = i == (y1 < y2 ? y1 : y2) ? 0.75 : i == (y2 > y1 ? y2 : y1) ? 0.25 : 1.0;
					CHECK_CELL(x < 0 ? 0 : x < 10 ? x : 9,
					           i,
						   (x < 0 ? (double)ydir : x < 10 ? (double)ydir : 0) * yhit,
					           (double)ydir * yhit);
				}
				CHECK_ZEROED();

				draw(fx + 0.25, fy1 + 0.25, fx + 0.25, fy2 + 0.25);
				for (i = 0; i < 10; i++) {
					if (i < y1 + (y1 > y2) && i < y2 + (y2 > y1))
						continue;
					if (i >= y1 + (y1 > y2) && i >= y2 + (y2 > y1))
						continue;
					yhit = i == (y1 < y2 ? y1 : y2) ? 0.75 : i == (y2 > y1 ? y2 : y1) ? 0.25 : 1.0;
					CHECK_CELL(x < 0 ? 0 : x < 10 ? x : 9,
					           i,
						   (x < 0 ? (double)ydir : x < 10 ? (double)ydir * 0.75 : 0) * yhit,
					           (double)ydir * yhit);
				}
				CHECK_ZEROED();
			}
		}
	}
}


static void
check_11_diagonal(ssize_t xshift, ssize_t yshift)
{
	ssize_t i1, i2, i, imin, imax, x, y;
	double fi1, fi2, idir;
	double fxshift = (double)xshift;
	double fyshift = (double)yshift;

	for (i1 = -2; i1 < 12; i1++) {
		for (i2 = -2; i2 < 12; i2++) {
			imin = i1 < i2 ? i1 : i2;
			imax = i1 > i2 ? i1 : i2;
			fi1 = (double)i1;
			fi2 = (double)i2;
			idir = SIGNUM(fi2 - fi1);

			draw(fi1 + fxshift, fi1 + fyshift, fi2 + fxshift, fi2 + fyshift);
			for (i = imin; i < imax; i++) {
				x = i + xshift;
				y = i + yshift;
				if (y < 0 || y >= 10)
					continue;
				if (x >= 10)
					CHECK_CELL(9, y, 0, idir);
				else if (x >= 0)
					CHECK_CELL(x, y, 0.5 * idir, idir);
				else
					CHECK_CELL(0, y, idir, idir);
			}
			CHECK_ZEROED();

			draw(fi1 + fxshift, fi1 + fyshift + 0.25, fi2 + fxshift, fi2 + fyshift + 0.25);
			for (i = imin; i < imax; i++) {
				x = i + xshift;
				y = i + yshift;
				if (0 <= y && y < 10) {
					if (x >= 10)
						SUB_CELL(9, y, 0, idir * 0.75);
					else if (x >= 0)
						SUB_CELL(x, y, idir * (0.75 - 0.75 * 0.75 / 2), idir * 0.75);
					else
						SUB_CELL(0, y, idir * 0.75, idir * 0.75);
				}
				x = i + xshift;
				y = i + yshift + 1;
				if (0 <= y && y < 10) {
					if (x >= 10)
						SUB_CELL(9, y, 0, idir * 0.25);
					else if (x >= 0)
						SUB_CELL(x, y, idir * 0.25 * 0.25 / 2, idir * 0.25);
					else
						SUB_CELL(0, y, idir * 0.25, idir * 0.25);
				}
			}
			CHECK_ZEROED();

			draw(fi1 + fxshift + 0.25, fi1 + fyshift, fi2 + fxshift + 0.25, fi2 + fyshift);
			for (i = imin; i < imax; i++) {
				x = i + xshift;
				y = i + yshift;
				if (0 <= y && y < 10) {
					if (x >= 10)
						SUB_CELL(9, y, 0, idir * 0.75);
					else if (x >= 0)
						SUB_CELL(x, y, idir * 0.75 * 0.75 / 2, idir * 0.75);
					else
						SUB_CELL(0, y, idir * 0.75, idir * 0.75);
				}
				x = i + xshift + 1;
				y = i + yshift;
				if (0 <= y && y < 10) {
					if (x >= 10)
						SUB_CELL(9, y, 0, idir * 0.25);
					else if (x >= 0)
						SUB_CELL(x, y, idir * (0.25 - 0.25 * 0.25 / 2), idir * 0.25);
					else
						SUB_CELL(0, y, idir * 0.25, idir * 0.25);
				}
			}
			CHECK_ZEROED();

			draw(fi1 + fxshift + 0.25, fi1 + fyshift + 0.25, fi2 + fxshift + 0.25, fi2 + fyshift + 0.25);
			for (i = imin; i < imax; i++) {
				x = i + xshift;
				y = i + yshift;
				if (0 <= y && y < 10) {
					if (x >= 10)
						SUB_CELL(9, y, 0, idir * 0.75);
					else if (x >= 0)
						SUB_CELL(x, y, idir * 0.75 * 0.75 / 2, idir * 0.75);
					else
						SUB_CELL(0, y, idir * 0.75, idir * 0.75);
				}
				x = i + xshift + 1;
				y = i + yshift + 1;
				if (0 <= y && y < 10) {
					if (x >= 10)
						SUB_CELL(9, y, 0, idir * 0.25);
					else if (x >= 0)
						SUB_CELL(x, y, idir * (0.25 - 0.25 * 0.25 / 2), idir * 0.25);
					else
						SUB_CELL(0, y, idir * 0.25, idir * 0.25);
				}
			}
			CHECK_ZEROED();
		}
	}
}


static void
check_11_antidiagonal(ssize_t xshift, ssize_t yshift)
{
	ssize_t i1, i2, i, imin, imax, x, y;
	double fi1, fi2, ydir;
	double fxshift = (double)xshift;
	double fyshift = (double)yshift;

	for (i1 = -2; i1 < 12; i1++) {
		for (i2 = -2; i2 < 12; i2++) {
			imin = i1 < i2 ? i1 : i2;
			imax = i1 > i2 ? i1 : i2;
			fi1 = (double)i1;
			fi2 = (double)i2;
			ydir = -SIGNUM(fi2 - fi1);

			draw(fi1 + fxshift, fi2 + fyshift, fi2 + fxshift, fi1 + fyshift);
			for (i = imin; i < imax; i++) {
				x = i + xshift;
				y = imax - i + imin + yshift - 1;
				if (y < 0 || y >= 10)
					continue;
				if (x >= 10)
					CHECK_CELL(9, y, 0, ydir);
				else if (x >= 0)
					CHECK_CELL(x, y, 0.5 * ydir, ydir);
				else
					CHECK_CELL(0, y, ydir, ydir);
			}
			CHECK_ZEROED();

			draw(fi1 + fxshift, fi2 + fyshift + 0.25, fi2 + fxshift, fi1 + fyshift + 0.25);
			for (i = imin; i < imax; i++) {
				x = i + xshift;
				y = imax - i + imin + yshift - 1;
				if (0 <= y && y < 10) {
					if (x >= 10)
						SUB_CELL(9, y, 0, ydir * 0.75);
					else if (x >= 0)
						SUB_CELL(x, y, ydir * 0.75 * 0.75 / 2, ydir * 0.75);
					else
						SUB_CELL(0, y, ydir * 0.75, ydir * 0.75);
				}
				x = i + xshift;
				y = imax - i + imin + yshift - 1 + 1;
				if (0 <= y && y < 10) {
					if (x >= 10)
						SUB_CELL(9, y, 0, ydir * 0.25);
					else if (x >= 0)
						SUB_CELL(x, y, ydir * (0.25 - 0.25 * 0.25 / 2), ydir * 0.25);
					else
						SUB_CELL(0, y, ydir * 0.25, ydir * 0.25);
				}
			}
			CHECK_ZEROED();

			draw(fi1 + fxshift + 0.25, fi2 + fyshift, fi2 + fxshift + 0.25, fi1 + fyshift);
			for (i = imin; i < imax; i++) {
				x = i + xshift;
				y = imax - i + imin + yshift - 1;
				if (0 <= y && y < 10) {
					if (x >= 10)
						SUB_CELL(9, y, 0, ydir * 0.75);
					else if (x >= 0)
						SUB_CELL(x, y, ydir * 0.75 * 0.75 / 2, ydir * 0.75);
					else
						SUB_CELL(0, y, ydir * 0.75, ydir * 0.75);
				}
				x = i + xshift + 1;
				y = imax - i + imin + yshift - 1;
				if (0 <= y && y < 10) {
					if (x >= 10)
						SUB_CELL(9, y, 0, ydir * 0.25);
					else if (x >= 0)
						SUB_CELL(x, y, ydir * (0.25 - 0.25 * 0.25 / 2), ydir * 0.25);
					else
						SUB_CELL(0, y, ydir * 0.25, ydir * 0.25);
				}
			}
			CHECK_ZEROED();

			draw(fi1 + fxshift + 0.25, fi2 + fyshift + 0.25, fi2 + fxshift + 0.25, fi1 + fyshift + 0.25);
			for (i = imin; i < imax; i++) {
				x = i + xshift;
				y = imax - i + imin + yshift - 1;
				if (0 <= y && y < 10) {
					if (x >= 10)
						SUB_CELL(9, y, 0, ydir * 0.5);
					else if (x >= 0)
						SUB_CELL(x, y, ydir / 8.0, ydir * 0.5);
					else
						SUB_CELL(0, y, ydir * 0.5, ydir * 0.5);
				}
				x = i + xshift;
				y = imax - i + imin + yshift - 1 + 1;
				if (0 <= y && y < 10) {
					if (x >= 10)
						SUB_CELL(9, y, 0, ydir * 0.25);
					else if (x >= 0)
						SUB_CELL(x, y, ydir * (2.5 / 16.0), ydir * 0.25);
					else
						SUB_CELL(0, y, ydir * 0.25, ydir * 0.25);
				}
				x = i + xshift + 1;
				y = imax - i + imin + yshift - 1;
				if (0 <= y && y < 10) {
					if (x >= 10)
						SUB_CELL(9, y, 0, ydir * 0.25);
					else if (x >= 0)
						SUB_CELL(x, y, ydir * (3.5 / 16.0), ydir * 0.25);
					else
						SUB_CELL(0, y, ydir * 0.25, ydir * 0.25);
				}
			}
			CHECK_ZEROED();
		}
	}
}


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
try_random_line(double xmin, double xmax, double ymin, double ymax)
{
	double x1 = frand(xmin, xmax);
	double y1 = frand(ymin, ymax);
	double x2 = frand(xmin, xmax);
	double y2 = frand(ymin, ymax);
	draw(x1, y1, x2, y2);
}


static void
check_random_line(double rand_xmin, double rand_xmax, double rand_ymin, double rand_ymax)
{
	double x1 = frand(rand_xmin, rand_xmax);
	double y1 = frand(rand_ymin, rand_ymax);
	double x2 = frand(rand_xmin, rand_xmax);
	double y2 = frand(rand_ymin, rand_ymax);
	double dx = x2 - x1;
	double dy = y2 - y1;
	double ymin = fmin(y1, y2);
	double ymax = fmax(y1, y2);
	double xmin = fmin(x1, x2);
	double xmax = fmax(x1, x2);
	double ydir = SIGNUM(y2 - y1);
	double inklevel, area, expected_area, inkdir, dist;
	size_t y, x, jumps;
	int hit;

	draw(x1, y1, x2, y2);

	area = 0;
	for (y = 0; y < 10; y++) {
		inklevel = 0;
		jumps = 0;
		inkdir = 0;
		for (x = 0; x < 10; x++) {
			area += raster->cells[y * 10 + x].cell_coverage + inklevel;
			inklevel += raster->cells[y * 10 + x].opposite_coverage;

			if (inkdir == 0) {
				inkdir = raster->cells[y * 10 + x].opposite_coverage;
			} else if (inkdir < 0 && raster->cells[y * 10 + x].opposite_coverage > 0) {
				ASSERT(!jumps);
				jumps += 1;
				inkdir = raster->cells[y * 10 + x].opposite_coverage;
			} else if (inkdir > 0 && raster->cells[y * 10 + x].opposite_coverage < 0) {
				ASSERT(!jumps);
				jumps += 1;
				inkdir = raster->cells[y * 10 + x].opposite_coverage;
			}

			if (!dy)
				continue;

			hit  = (double)x >= floor(xmin) && (double)x <= ceil(xmax);
			hit |= !x && xmin < 0;
			hit &= (double)y >= floor(ymin) && (double)y <= ceil(ymax);
			ASSERT(hit || !raster->cells[y * 10 + x].cell_coverage);

			hit  = (double)x >= floor(xmin) && (double)x <= ceil(xmax);
			hit &= (double)y >= floor(ymin) && (double)y <= ceil(ymax);
			if (!hit || !dx || !x)
				continue;

			dist = ((double)x+0.5)*dy/dx - ((double)y+0.5) + (y1 - x1*dy/dx);
			dist /= hypot(dy/dx, 1);
			dist = fabs(dist);
			if (dist > 0.708)
				ASSERT(!raster->cells[y * 10 + x].cell_coverage);
		}

		if ((double)y < floor(ymin) || (double)y >= ceil(ymax))
			ASSERT(inklevel == 0.0);
		else if ((double)y == floor(ymin))
			ASSERT(eq(inklevel, ydir * (fmin(ymax, (double)y + 1.0) - ymin)));
		else if ((double)y == floor(ymax))
			ASSERT(eq(inklevel, ydir * fmod(ymax, 1.0)));
		else
			ASSERT(eq(inklevel, ydir));
	}

	if (ymax <= 0 || ymin >= 10 || !dy || fmin(x1, x2) >= 10) {
		expected_area = 0;
	} else if (fmax(x1, x2) < 0) {
		expected_area = ydir * 10 * (fmin(10, ymax) - fmax(ymin, 0));
	} else if (!dx) {
		expected_area = ydir * (10 - fmin(fmax(0, x1), 10)) * (fmin(10, ymax) - fmax(ymin, 0));
	} else {
#define y(x) (y1 + ((x) - x1) * k)
#define x(y) (x1 + ((y) - y1) / k)
		double k = dy / dx;
		double a = xmin;
		double b = xmax;
		double mid;

		a = y(a) > 10 ? x(10) : y(a) < 0 ? x(0) : a;
		b = y(b) > 10 ? x(10) : y(b) < 0 ? x(0) : b;
		b = fmin(b, 10);
		a = fmin(a, 10);
		mid = (a + b) / 2;

		if (mid < 0) {
			expected_area = (10 - 0) * fabs(y(b) - y(a));
			if (b > 0)
				expected_area -= fabs(y(b) - y(0)) * (b - 0) / 2;
		} else {
			expected_area = (10 - mid) * fabs(y(b) - y(a));
			if (a < 0)
				expected_area -= fabs(y(0) - y(a)) * (0 - a) / 2;
		}
		expected_area *= ydir;
#undef y
#undef x
	}
	ASSERT(eq(area, expected_area));

	rtgrpblib_reset_raster(raster, 10, 10);
}


int
main(void)
{
	size_t i;

	raster = rtgrpblib_create_raster(10, 10);
	ASSERT(raster);
	srand((unsigned)(uintptr_t)raster);

	check_horizontal();
	check_vertical();

	check_11_diagonal(0, 0);
	check_11_diagonal(0, 2);
	check_11_diagonal(0, -2);
	check_11_diagonal(2, 0);
	check_11_diagonal(2, 2);
	check_11_diagonal(2, -2);
	check_11_diagonal(-2, 0);
	check_11_diagonal(-2, 2);
	check_11_diagonal(-2, -2);

	check_11_antidiagonal(0, 0);
	check_11_antidiagonal(0, 2);
	check_11_antidiagonal(0, -2);
	check_11_antidiagonal(2, 0);
	check_11_antidiagonal(2, 2);
	check_11_antidiagonal(2, -2);
	check_11_antidiagonal(-2, 0);
	check_11_antidiagonal(-2, 2);
	check_11_antidiagonal(-2, -2);

	for (i = 0; i < 100000UL; i++) {
		try_random_line(0., 1., 0., 1.);
		try_random_line(-5., 15., -5., 15.);
	}
	rtgrpblib_reset_raster(raster, 10, 10);

	for (i = 0; i < 10000UL; i++) {
		check_random_line(0., 1., 0., 1.);
		check_random_line(-5., 15., -5., 15.);
	}

	free(raster);
	return 0;
}


#endif