path: root/lines.c

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


#ifndef ROWWISE_RESET_INKLEVEL
void
draw_vertical_line_opposite_only(RASTER *restrict raster, size_t x, double y1, double y2, int ydir)
{
	size_t y, start, end;
	CELL *cells;

	cells = &raster->cells[x];

	start = (size_t)y1;
	end = (size_t)y2;

	y = start;
	if (ydir > 0) {
		cells[end * raster->width].opposite_coverage -= 1.0 - (y2 - (double)end);
		cells[start * raster->width].opposite_coverage -= y1 - (double)start;
		cells = &cells[start * raster->width];
		do {
			cells->opposite_coverage += 1.0;
			cells += raster->width;
		} while (y++ != end);
	} else {
		cells[end * raster->width].opposite_coverage += y2 - (double)end;
		cells[start * raster->width].opposite_coverage += 1.0 - (y1 - (double)start);
		cells = &cells[start * raster->width];
		do {
			cells->opposite_coverage -= 1.0;
			cells -= raster->width;
		} while (y-- != end);
	}
}
#endif


void
draw_vertical_line(RASTER *restrict raster, double x1, double y1, double y2, int ydir)
{
	size_t x, y, start, end;
	double cell_coverage;
	CELL *cells;

	x = (size_t)x1;
	cell_coverage = 1.0 - (x1 - (double)x);
	cells = &raster->cells[x];

	start = (size_t)y1;
	end = (size_t)y2;
	start -= (size_t)(start == raster->height);
	end -= (size_t)(end == raster->height);

	y = start;
	if (ydir > 0) {
		cells[end * raster->width].opposite_coverage -= 1.0 - (y2 - (double)end);
		cells[end * raster->width].cell_coverage -= (1.0 - (y2 - (double)end)) * cell_coverage;
		cells[start * raster->width].opposite_coverage -= y1 - (double)start;
		cells[start * raster->width].cell_coverage -= (y1 - (double)start) * cell_coverage;
		cells = &cells[start * raster->width];
		do {
			cells->opposite_coverage += 1.0;
			cells->cell_coverage += cell_coverage;
			cells += raster->width;
		} while (y++ != end);
	} else {
		cells[end * raster->width].opposite_coverage += y2 - (double)end;
		cells[end * raster->width].cell_coverage += (y2 - (double)end) * cell_coverage;
		cells[start * raster->width].opposite_coverage += 1.0 - (y1 - (double)start);
		cells[start * raster->width].cell_coverage += (1.0 - (y1 - (double)start)) * cell_coverage;
		cells = &cells[start * raster->width];
		do {
			cells->opposite_coverage -= 1.0;
			cells->cell_coverage -= cell_coverage;
			cells -= raster->width;
		} while (y-- != end);
	}
}


void
draw_diagonal_line(RASTER *restrict raster, double x1, double y1, double x2, double y2,
                   double dx, double dy, int xdir, int ydir)
{
	/* This code is taken, with rewrites and added documentation,
	 * from Thomas Oltmann's ISC licensed libschrift. */

	ssize_t rowsize;
	CELL *cell;
	size_t startCellX, startCellY, step, nsteps = 0;
	double ydiff, halfdx, one_minus_x1_plus_column;
	double nextProgress, prevProgress;
	double nextXProgress, nextYProgress;
	double xDistanceToProgress, yDistanceToProgress;
	int alongX;

	/* See just beneath */
	xDistanceToProgress = fabs(1.0 / dx);
	yDistanceToProgress = fabs(1.0 / dy);

	/* Get the column (x1, y1) is in. This is nothing
	 * more complicated then the floor of `x1`, however
	 * if `x1 > x2` and `x1` is on a cell edge, we say
	 * that it is in the left cell rather than the
	 * right cell; hence `ceil(x1) - 1`.
	 * Then we figure out how many steps will take
	 * horizontally when drawing the line. If (x1, y1)
	 * and (x2, y2) are in the same column, this is 0.
	 * Hence this the difference between the floor of
	 * min(x1, x2) and ceiling of max(x1, x2), less 1;
	 * less 1 because if (x1, y1) and (x2, y2) are in
	 * the same column the difference is 1: it is how
	 * many cells wide are covered, which is the 1
	 * greater than the number of edges between the
	 * cells.
	 * Third, we calculate (`nextXProgress`) where the
	 * line next (from (x1, y1)) intersections with a
	 * vertical line in the cell grid, and convert
	 * it the horizontal progression towards (x2, y2)
	 * by dividing by the unsigned projection of the
	 * line onto the horizontal raster axis.
	 */
	if (xdir > 0) {
		double x1floor = floor(x1), x2ceil = ceil(x2);
		startCellX = (size_t)x1floor;
		nsteps += (size_t)x2ceil - (size_t)x1floor - 1;
		nextXProgress = (x1floor + 1.0 - x1) * xDistanceToProgress;
	} else {
		double x2floor = floor(x2), x1ceil = ceil(x1);
		startCellX = (size_t)x1ceil - 1;
		nsteps += (size_t)x1ceil - (size_t)x2floor - 1;
		nextXProgress = (x1 + 1.0 - x1ceil) * xDistanceToProgress;
	}

	/* Same as above, but on the vertical axis */
	if (ydir > 0) {
		double y1floor = floor(y1), y2ceil = ceil(y2);
		startCellY = (size_t)y1floor;
		nsteps += (size_t)y2ceil - (size_t)y1floor - 1;
		nextYProgress = (y1floor + 1.0 - y1) * yDistanceToProgress;
	} else {
		double y2floor = floor(y2), y1ceil = ceil(y1);
		startCellY = (size_t)y1ceil - 1;
		nsteps += (size_t)y1ceil - (size_t)y2floor - 1;
		nextYProgress = (y1 + 1.0 - y1ceil) * yDistanceToProgress;
	}

	/* These values are in [0, 1], where 0 corresponds to (x1, y1)
	 * and 1 corresponds to (x2, y2). `prevProgress` is the amount of
	 * work done at the beginning of for-loop body, and `nextProgress`
	 * is the amount of work done at the end of the for-loop body. */
	prevProgress = 0.0;
	nextProgress = fmin(nextXProgress, nextYProgress);

	/* These are just tricks to reduce the number of operations made */
	halfdx = 0.5 * dx;
	one_minus_x1_plus_column = 1.0 - x1 + (double)startCellX;
	cell = &raster->cells[startCellY * raster->width + startCellX];
	rowsize = ydir * (ssize_t)raster->width;

	/* Visit all cells the line intersect (and one extra for every time
	 * the line intersects a cell corner); `nsteps` is the number of
	 * horizontal steps plus the number of vertical steps (Manhattan walk). */
	for (step = 0; step < nsteps; step++) {
		/* `nextProgress` and `prevProgress` are [0, 1] values
		 * measuring the progression from (x1, y1) to (x2, y2),
		 * hence `nextProgress - prevProgress` the the progress
		 * this interation makes, and it multiplied by `dy` is
		 * how much we travel vertically. This is how much
		 * additional ink shall be applied to all cells on this
		 * cells right (in the same row). */
		cell->opposite_coverage += ydiff = (nextProgress - prevProgress) * dy;
		/* This is the mount of ink we need in in this cell
		 * (minus all the sum of `.opposite_coverage` of
		 * all cells on the left in the same row). This is
		 * the area of a right trapezoid made by cells right
		 * edge, the line `y = prevProgress * dy`, the line
		 * `y = nextProgress * dy` and the draw line. */
		cell->cell_coverage += ydiff * (one_minus_x1_plus_column - halfdx * (prevProgress + nextProgress));

		/* Step either vertically or horizontally, whichever is shortest */
		prevProgress = nextProgress;
		alongX = nextXProgress < nextYProgress;
		nextXProgress += alongX ? xDistanceToProgress : 0;
		nextYProgress += alongX ? 0 : yDistanceToProgress;
		cell += alongX ? xdir : 0;
		cell += alongX ? 0 : rowsize;
		one_minus_x1_plus_column += alongX ? xdir : 0;
		nextProgress = fmin(nextXProgress, nextYProgress);
	}

	/* In case we get beyond (x2, y2), that is, if `prevProgress > 1`. */
	cell->opposite_coverage += ydiff = (1.0 - prevProgress) * dy;
	cell->cell_coverage += ydiff * (one_minus_x1_plus_column - halfdx * (prevProgress + 1.0));
}


void
draw_bounded_line(RASTER *restrict raster, double x1, double y1, double x2, double y2)
{
	double dx, dy = y2 - y1;
	int xdir, ydir = SIGNUM(dy);

	if (!ydir)
		return;

	dx = x2 - x1;
	xdir = SIGNUM(dx);

	if (!xdir)
		draw_vertical_line(raster, x1, y1, y2, ydir);
	else
		draw_diagonal_line(raster, x1, y1, x2, y2, dx, dy, xdir, ydir);
}


#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_against_reference(double x1, double y1, double x2, double y2)
{
	size_t i;
	alarm(5);
	draw_linear_bezier_reference(refraster, x1, y1, x2, y2);
	alarm(0);
	for (i = 0; i < 100; i++) {
		ASSERT(intolerant_eq(raster->cells[i].cell_coverage, refraster->cells[i].cell_coverage));
		ASSERT(intolerant_eq(raster->cells[i].opposite_coverage, refraster->cells[i].opposite_coverage));
	}
}


static void
check_vertical_line(double x, double y1, double y2)
{
	alarm(5);
	draw_vertical_line(raster, x, y1, y2, SIGNUM(y2 - y1));
	alarm(0);
	check_against_reference(x, y1, x, y2);
}


static void
check_diagonal_line(double x1, double y1, double x2, double y2)
{
	double dx = x2 - x1;
	double dy = y2 - y1;
	alarm(5);
	draw_diagonal_line(raster, x1, y1, x2, y2, dx, dy, SIGNUM(dx), SIGNUM(dy));
	alarm(0);
	check_against_reference(x1, y1, x2, y2);
}


static void
check_bounded_line(double x1, double y1, double x2, double y2)
{
	alarm(5);
	draw_bounded_line(raster, x1, y1, x2, y2);
	alarm(0);
	check_against_reference(x1, y1, x2, y2);
}


#ifndef ROWWISE_RESET_INKLEVEL
static void
check_vertical_line_opposite_only(double y1, double y2)
{
	alarm(5);
	draw_vertical_line_opposite_only(raster, 9, y1, y2, SIGNUM(y2 - y1));
	alarm(0);
	check_against_reference(20, y1, 20, y2);
}
#endif


static void
reset(void)
{
	rtgrpblib_reset_raster(raster, 10, 10);
	rtgrpblib_reset_raster(refraster, 10, 10);
}


int
main(void)
{
	double x1, x2;
	double y1, y2;
	size_t i;

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

	for (i = 0; i < 10000UL; i++) {
		x1 = frand(0, 10);
		y1 = frand(0, 10);
		x2 = frand(0, 10);
		y2 = frand(0, 10);

		check_vertical_line(x1, y1, y2);
		reset();

		if (x1 != x2) {
			check_diagonal_line(x1, y1, x2, y2);
			reset();
		}

		check_bounded_line(x1, y1, x1, y2);
		reset();

		check_bounded_line(x1, y1, x2, y2);
		reset();

#ifndef ROWWISE_RESET_INKLEVEL
		check_vertical_line_opposite_only(y1, y2);
		reset();
#endif
	}

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


#endif