/* 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.
 */


static void
draw_cell(RASTER *restrict raster, double x1, double y1, double x2, double y2, double dx, double dy,
          double cxmin, double cxmax, double rxmin, double rxmax, double rymin, double rymax, ssize_t x, size_t cellY)
{
	double cy1, cy2, ydir = (double)SIGNUM(dy);
	double rightedge, opposite, midadjacent;
	size_t cell;

	rightedge = cxmax;
	cxmin = fmax(cxmin, rxmin);
	cxmax = fmin(cxmax, rxmax);
	cy1 = dx ? y1 + (cxmin - x1) * dy / dx : y1;
	cy2 = dx ? y2 + (cxmax - x2) * dy / dx : y2;
	cy1 = fmin(fmax(cy1, rymin), rymax);
	cy2 = fmin(fmax(cy2, rymin), rymax);

	opposite = ydir * fabs(cy2 - cy1);
	midadjacent = (cxmin + cxmax) / 2.0;

	if (x < 0) {
		raster->cells[cellY].opposite_coverage += opposite;
		raster->cells[cellY].cell_coverage += opposite;
	} else if ((size_t)x >= raster->width) {
		cell = cellY + (raster->width - 1);
		raster->cells[cell].opposite_coverage += opposite;
	} else {
		cell = cellY + (size_t)x;
		raster->cells[cell].opposite_coverage += opposite;
		raster->cells[cell].cell_coverage += opposite * fabs(rightedge - midadjacent);
	}
}

void
draw_linear_bezier_reference(RASTER *restrict raster, double x1, double y1, double x2, double y2)
{
	ssize_t iwidth = (ssize_t)raster->width;
	double width = (double)raster->width;
	double height = (double)raster->height;
	double dx = x2 - x1;
	double dy = y2 - y1;
	double fymin = fmin(fmax(fmin(y1, y2), 0.0), height);
	double fymax = fmin(fmax(fmax(y1, y2), 0.0), height);
	double floor_fymin = floor(fymin);
	double ceil_fymax = ceil(fymax);
	size_t ymin = (size_t)floor_fymin;
	size_t ymax = (size_t)ceil_fymax;
	size_t y, cellY;
	double rowx1, rowx2;
	double rxmin, rxmax, floor_rxmin;
	double rymax, rymin, ceil_rxmax;
	ssize_t x, xmin, xmax, xmin_cropped, xmax_cropped;
	double cxmin, cxmax;

	ymax += (size_t)(ymax == ymin);
	for (y = ymin; y < ymax; y++) {
		rymin = (double)y;
		rymax = (double)y + 1.0;
		rymin = fmax(rymin, fymin);
		rymax = fmin(rymax, fymax);
		rowx1 = dy ? x1 + (rymin - y1) * dx / dy : x1;
		rowx2 = dy ? x2 + (rymax - y2) * dx / dy : x2;
		rxmin = fmin(rowx1, rowx2);
		rxmax = fmax(rowx1, rowx2);
		cellY = y * raster->width;

		floor_rxmin = floor(rxmin);
		ceil_rxmax = ceil(rxmax);
		xmin = (ssize_t)floor_rxmin;
		xmax = (ssize_t)ceil_rxmax;
		xmax += (ssize_t)(xmax == xmin);
		xmin_cropped = xmin;
		xmax_cropped = xmax;
		if (xmin < 0.0) {
			cxmin = (double)xmin;
			cxmax = xmax < 0 ? (double)(xmax + 1) : 0.0;
			draw_cell(raster, x1, y1, x2, y2, dx, dy, cxmin, cxmax, rxmin, rxmax, rymin, rymax, xmin, cellY);
			xmin_cropped = 0;
		}
		if (xmax > iwidth) {
			cxmin = fmax(floor_rxmin, width);
			cxmax = (double)xmax;
			draw_cell(raster, x1, y1, x2, y2, dx, dy, cxmin, cxmax, rxmin, rxmax, rymin, rymax, (ssize_t)cxmin, cellY);
			xmax_cropped = iwidth;
		}
		for (x = xmin_cropped; x < xmax_cropped; x++) {
			cxmin = (double)x;
			cxmax = (double)x + 1.0;
			draw_cell(raster, x1, y1, x2, y2, dx, dy, cxmin, cxmax, rxmin, rxmax, rymin, rymax, x, cellY);
		}
	}
}


#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