First commit

Signed-off-by: Mattias Andrée <maandree@kth.se>

1 files changed, 227 insertions, 0 deletions

diff --git a/lines.c b/lines.c
new file mode 100644
index 0000000..dfe2448
--- /dev/null
+++ b/lines.c
@@ -0,0 +1,227 @@
+/* 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;
+
+	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
+
+
+int
+main(void)
+{
+	return 0; /* TODO add test */
+}
+
+
+#endif