fix(display): rewrite line routine to use implicit function

The old approach was using distance based error, the implicit function
approach yields slightly different lines which are more agreeable with
the WIP triangle branch.

Author: ulahello

display/rendering/canvas.c

@@ -127,88 +127,54 @@ DEFN_RENDER(circle)
 
 DEFN_RENDER(line)
 {
-	const struct color color = line->c;
-
-	// I think this is essentially Bresenham's line algorithm, as in, that
-	// is where the ideas came from, but with the caveat that I didn't want
-	// a separate vertical and horizontal drawing function.
-	//
-	// ## Terms
-	//
-	// Since this is octant-agnostic, we generalize the notion of a step
-	// between pixels.
+	// This approach is guided by https://zingl.github.io/Bresenham.pdf
+	// (A Rasterizing Algorithm for Drawing Curves, by Alois Zingl).
+	// The math is based on the implicit function for a line,
 	//
-	// Considering the vector (dx, dy), one component will have greater or
-	// equal magnitude than the other, so at every iteration, the pixel
-	// position will be incremented along that component. I refer to this
-	// increment vector as the step.
+	// > f(x, y) = (y - y0)(x1 - x0) - (x - x0)(y1 - y0)
 	//
-	// The pixel position may also be incremented along the component with
-	// lesser magnitude, but this happens a fraction of the time (a fraction
-	// in the range [0, 1]). I refer to this as the lesser step.
+	// which for all points on the line equals zero. Since the slope of a
+	// line remains constant, we can test the "error" (evaluating the
+	// function) at the next possible points and use a cheap comparison to
+	// decide which to move to with the least error.
 	//
-	// ## Process for deriving
+	// For lines, the errors for directly adjacent pixels can be expressed
+	// in terms of the error for the diagonal pixel and the total distances
+	// `dy` and `dx`. This means we only need one error variable.
 	//
-	// 1. Determine whether to step in the lesser direction.
-	//    - Given the two possible future points (either cur + step or cur +
-	//      step + lstep), and the point on the line which lies between
-	//      them, which future point is closer to the line? This calculation
-	//      is made with real numbers, approximated with floats.
-	// 2. Scale the two distances such that the comparison is still
-	//    meaningful, while eliminating floating point division. Now, we can
-	//    use integers!
-	// 3. Avoid recomputing the distances from scratch each iteration.
-	//    - Determine the initial value, and the difference between the
-	//      (i+1)th and ith values for each branch.
-
-	const int32_t dx = (int32_t)line->x1 - (int32_t)line->x0;
-	const int32_t dy = (int32_t)line->y1 - (int32_t)line->y0;
-
-	const uint16_t steps = (uint16_t)max(abs(dx), abs(dy));
-
-	if (steps == 0) {
-		draw_point(c, line->x0, line->y0, color);
-		return;
-	}
-
-	// step_x and step_y are both in the range [-1, 1], and
-	// at least one is -1 or 1
-	const int32_t step_x = dx / steps;
-	const int32_t step_y = dy / steps;
+	// The next optimization is to track the initial value and difference
+	// per iteration, rather than recomputing the error each time.
 
-	const int32_t lstep_x = step_x == 0 ? (dx != 0 ? dx / abs(dx) : 0) : 0;
-	const int32_t lstep_y = step_y == 0 ? (dy != 0 ? dy / abs(dy) : 0) : 0;
+	const int32_t dx = abs(line->x1 - line->x0);
+	const int32_t dy = abs(line->y1 - line->y0);
+	const int32_t sx = line->x0 < line->x1 ? 1 : -1;
+	const int32_t sy = line->y0 < line->y1 ? 1 : -1;
 
-	uint16_t x = line->x0;
-	uint16_t y = line->y0;
+	int32_t x = line->x0;
+	int32_t y = line->y0;
 
-	int32_t d0_x = steps * step_x;
-	int32_t d0_y = steps * step_y;
-	int32_t d1_x = -steps * (step_x + lstep_x);
-	int32_t d1_y = -steps * (step_y + lstep_y);
-
-	for (int32_t i = 0; i <= steps; i++) {
-		draw_point(c, x, y, color);
-
-		const int32_t px = abs(d0_x) - abs(d1_x);
-		const int32_t py = abs(d0_y) - abs(d1_y);
-
-		x += step_x;
-		d0_x += steps * step_x - dx;
-		d1_x -= steps * step_x - dx;
-		if (px >= 0) {
-			x += lstep_x;
-			d0_x += steps * lstep_x;
-			d1_x -= steps * lstep_x;
+	// We track the error for the next diagonal pixel (x + sx, y + sy).
+	//
+	// This relies on the assumption that the slope is positive, but we took
+	// the absolute value for `dy` and `dx`, and this lie is self contained
+	// in the mathy state and accounted for by `sx` and `sy`.
+	int32_t e = dx - dy;
+
+	while (true) {
+		draw_point(c, x, y, line->c);
+
+		// Check if we reached the other end.
+		if (x == line->x1 && y == line->y1)
+			break;
+
+		const int32_t test = 2 * e;
+		if (test > -dy) {
+			x += sx;
+			e -= dy;
 		}
-
-		y += step_y;
-		d0_y += steps * step_y - dy;
-		d1_y -= steps * step_y - dy;
-		if (py >= 0) {
-			y += lstep_y;
-			d0_y += steps * lstep_y;
-			d1_y -= steps * lstep_y;
+		if (test < dx) {
+			y += sy;
+			e += dx;
 		}
 	}
 }

display/testing.c

@@ -33,6 +33,7 @@ static void run_these_tests(
 static void test_rects(struct canvas *c);
 static void test_circles(struct canvas *c);
 static void test_lines_burst(struct canvas *c);
+static void test_lines_array(struct canvas *c);
 static void test_copy_rect(struct canvas *c);
 static void test_bezier2(struct canvas *c);
 
@@ -60,6 +61,13 @@ void run_tests(const char *dump_dir)
 		    .width = 32,
 		    .height = 32,
 		},
+		{
+		    .fill_color = BG,
+		    .draw_fn = test_lines_array,
+		    .output_path = "lines-array.data",
+		    .width = 128,
+		    .height = 128,
+		},
 		{
 		    .fill_color = BG,
 		    .draw_fn = test_copy_rect,
@@ -229,6 +237,45 @@ static void test_lines_burst(struct canvas *c)
 	}
 }
 
+static void test_lines_array(struct canvas *c)
+{
+	const int32_t cols = 7;
+	const int32_t rows = 13;
+	const int32_t gap = 1;
+
+	const int32_t col_len = c->width / cols;
+	const int32_t row_len = c->height / rows;
+
+	for (int32_t col = 0; col < cols; col++) {
+		for (int32_t row = 0; row < rows; row++) {
+			int32_t col_start = col * c->width / cols;
+			int32_t row_start = row * c->height / rows;
+			int32_t cx = col_start + col_len / 2;
+			int32_t cy = row_start + row_len / 2;
+
+			float p = (float)(cols * row + col) /
+			    ((float)cols * (float)rows);
+
+			rendering_draw_line(c,
+			    &(struct line) {
+				.x0 = cx +
+				    roundf(cosf(p * TAU) *
+					((float)col_len / 2 - gap)),
+				.y0 = cy +
+				    roundf(sinf(p * TAU) *
+					((float)row_len / 2 - gap)),
+				.x1 = cx +
+				    roundf(cosf(p * TAU + PI) *
+					((float)col_len / 2 - gap)),
+				.y1 = cy +
+				    roundf(sinf(p * TAU + PI) *
+					((float)row_len / 2 - gap)),
+				.c = FG,
+			    });
+		}
+	}
+}
+
 static void test_copy_rect(struct canvas *c)
 {
 	// Fill the screen with something.
@@ -308,8 +355,8 @@ static void test_bezier2(struct canvas *c)
 
 	for (int32_t col = 0; col < cols; col++) {
 		for (int32_t row = 0; row < rows; row++) {
-			int32_t col_start = col * col_len;
-			int32_t row_start = row * row_len;
+			int32_t col_start = col * c->width / cols;
+			int32_t row_start = row * c->height / rows;
 			int32_t dx = col * col_len / cols;
 			int32_t dy = row * row_len / rows;