Test and fix rtgrpblib_draw_quadratic_bezier for perfect lines

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

1 files changed, 194 insertions, 12 deletions

diff --git a/rtgrpblib_draw_quadratic_bezier.c b/rtgrpblib_draw_quadratic_bezier.c
index 8aae11e..e2202c6 100644
--- a/rtgrpblib_draw_quadratic_bezier.c
+++ b/rtgrpblib_draw_quadratic_bezier.c
@@ -90,23 +90,22 @@ rtgrpblib_draw_quadratic_bezier(RASTER *restrict raster, double x1, double y1,
 	size_t nts = 0;
 	size_t i;
 
-#ifdef TODO /* untested */
 	/* Can we downgrade the curve to a linear Bézier curve?
 	 *
 	 *  (y2 - y1)/(x2 - x1) = (y3 - y2)/(x3 - x2) = (y3 - y1)/(x3 - x1)
 	 *  (y2 - y1)/(x2 - x1) = (y3 - y2)/(x3 - x2)
 	 *  (y2 - y1)(x3 - x2) = (y3 - y2)(x2 - x1)
 	 */
-	if ((y3 - y1) * (x3 - x2) == (y3 - y2) * (x2 - x1)) {
-		if (x1 <= x2 && x2 <= x3 && y1 <= y2 && y2 <= y3)
-			rtgrpblib_draw_linear_bezier(raster, x1, y1, x3, y3);
-		else if (x1 <= x3 && x3 <= x2 && y1 <= y3 && y3 <= y2)
-			rtgrpblib_draw_linear_bezier(raster, x1, y1, x2, y2);
-		else
-			rtgrpblib_draw_linear_bezier(raster, x2, y2, x3, y3);
+	if ((y2 - y1) * (x3 - x2) == (y3 - y2) * (x2 - x1)) {
+		/* Drawing (x1, y1)--(x2, y2) plus (x2, y2)--(x3, y3)
+		 * is equivalent to drawing (x1, y1)--(x3, y3) because
+		 * the result is signed and if (x2, y2) is outside the
+		 * (x1, y1)--(x3, y3) line (x2, y2)--(x3, y3) erases
+		 * the part that (x1, y1)--(x2, y2) extends to
+		 * (x1, y1)--(x3, y3) */
+		rtgrpblib_draw_linear_bezier(raster, x1, y1, x3, y3);
 		return;
 	}
-#endif
 
 	/* Beginning and end of curve */
 	ts[nts++] = 0.0;
@@ -126,12 +125,12 @@ rtgrpblib_draw_quadratic_bezier(RASTER *restrict raster, double x1, double y1,
 
 		/* Remove any segments above or below the raster */
 		y = evaluate_quadratic_bezier(t, y1, y2, y3);
-		if (y < 0 || y >= (double)raster->height)
+		if (y < 0 || y > (double)raster->height)
 			continue;
 
 		/* If the segment is inside the raster, draw it, */
 		x = evaluate_quadratic_bezier(t, x1, x2, x3);
-		if (0 <= x && x < (double)raster->width) {
+		if (0 <= x && x <= (double)raster->width) {
 			draw_bounded_quadratic_bezier(raster, x1, y1, x2, y2, x3, y3, t1, t2);
 			continue;
 		}
@@ -157,10 +156,193 @@ rtgrpblib_draw_quadratic_bezier(RASTER *restrict raster, double x1, double y1,
 #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_line(double x1, double y1, double x2, double y2, double x3, double y3, int exact)
+{
+	size_t i;
+
+	alarm(5);
+	rtgrpblib_draw_quadratic_bezier(raster, x1, y1, x2, y2, x3, y3);
+
+	alarm(5);
+	rtgrpblib_draw_linear_bezier(refraster, x1, y1, x3, y3);
+
+	alarm(0);
+
+	if (exact) {
+		ASSERT(!memcmp(raster->cells, refraster->cells, 100 * sizeof(*raster->cells)));
+	} else {
+		for (i = 0; i < 100; i++) {
+			if (!eq(raster->cells[i].cell_coverage, refraster->cells[i].cell_coverage))
+				fprintf(stderr, "c%02zu: %la ref: %la; diff: %lg=%la (tol: %la)\n", i, raster->cells[i].cell_coverage, refraster->cells[i].cell_coverage, refraster->cells[i].cell_coverage - raster->cells[i].cell_coverage, refraster->cells[i].cell_coverage - raster->cells[i].cell_coverage, TOLERANCE);
+			if (!eq(raster->cells[i].opposite_coverage, refraster->cells[i].opposite_coverage))
+				fprintf(stderr, "o%02zu: %la ref: %la; diff: %lg=%la (tol: %la)\n", i, raster->cells[i].opposite_coverage, refraster->cells[i].opposite_coverage, raster->cells[i].opposite_coverage - refraster->cells[i].opposite_coverage, raster->cells[i].opposite_coverage - refraster->cells[i].opposite_coverage, TOLERANCE);
+		}
+		for (i = 0; i < 100; i++) {
+			if (!eq(raster->cells[i].cell_coverage, refraster->cells[i].cell_coverage) ||
+			    !eq(raster->cells[i].opposite_coverage, refraster->cells[i].opposite_coverage)) {
+				print_raster(refraster);
+				print_raster(raster);
+			}
+			ASSERT(eq(raster->cells[i].cell_coverage, refraster->cells[i].cell_coverage));
+			ASSERT(eq(raster->cells[i].opposite_coverage, refraster->cells[i].opposite_coverage));
+		}
+	}
+	rtgrpblib_reset_raster(raster, 10, 10);
+	rtgrpblib_reset_raster(refraster, 10, 10);
+}
+
+
+static void
+check_simple_lines(void)
+{
+	check_line(1, 1, 2, 2, 3, 3, 1);
+	check_line(1, 1, 3, 3, 2, 2, 1);
+	check_line(2, 2, 1, 1, 3, 3, 1);
+	check_line(2, 2, 3, 3, 1, 1, 1);
+	check_line(3, 3, 1, 1, 2, 2, 1);
+	check_line(3, 3, 2, 2, 1, 1, 1);
+	check_line(1, 1, 3, 3, 3, 3, 1);
+	check_line(1, 1, 1, 1, 3, 3, 1);
+
+	check_line(1, 4, 2, 4, 3, 4, 1);
+	check_line(1, 4, 3, 4, 2, 4, 1);
+	check_line(2, 4, 1, 4, 3, 4, 1);
+	check_line(2, 4, 3, 4, 1, 4, 1);
+	check_line(3, 4, 1, 4, 2, 4, 1);
+	check_line(3, 4, 2, 4, 1, 4, 1);
+	check_line(1, 4, 3, 4, 3, 4, 1);
+	check_line(1, 4, 1, 4, 3, 4, 1);
+
+	check_line(4, 1, 4, 2, 4, 3, 1);
+	check_line(4, 1, 4, 3, 4, 2, 1);
+	check_line(4, 2, 4, 1, 4, 3, 1);
+	check_line(4, 2, 4, 3, 4, 1, 1);
+	check_line(4, 3, 4, 1, 4, 2, 1);
+	check_line(4, 3, 4, 2, 4, 1, 1);
+	check_line(4, 1, 4, 3, 4, 3, 1);
+	check_line(4, 1, 4, 1, 4, 3, 1);
+}
+
+
+static void
+check_random_bounded_lines(int exact)
+{
+	double x1, x2, x3;
+	double y1, y2, y3;
+	double mid;
+	size_t i;
+
+	for (i = 0; i < 100000UL; i++) {
+		x1 = floor(frand(0, 11));
+		y1 = floor(frand(0, 11));
+		x3 = floor(frand(0, 11));
+		y3 = floor(frand(0, 11));
+		x1 -= (double)(x1 > 10);
+		x2 -= (double)(x2 > 10);
+		y1 -= (double)(y1 > 10);
+		y2 -= (double)(y2 > 10);
+		mid = frand(0, 1);
+
+		if (x3 - x1) {
+			x2 = floor(mid * (x3 - x1) + x1);
+			y2 = y1 + (x2 - x1) * (y3 - y1) / (x3 - x1);
+		} else if (y3 - y1) {
+			x2 = x1;
+			y2 = floor(mid * (y3 - y1) + y1);
+		} else {
+			x2 = x1;
+			y2 = y1;
+		}
+
+		if (!exact)
+			check_line(x1, y1, x2, y2, x3, y3, x2 == floor(x2) && y2 == floor(y2));
+		else if (x2 == floor(x2) && y2 == floor(y2))
+			check_line(x1, y1, x2, y2, x3, y3, 1);
+	}
+}
+
+
+static void
+check_random_unbounded_lines(int exact)
+{
+	double x1, x2, x3;
+	double y1, y2, y3;
+	double mid;
+	size_t i;
+
+	for (i = 0; i < 100000UL; i++) {
+		x1 = floor(frand(-5, 15));
+		y1 = floor(frand(-5, 15));
+		x3 = floor(frand(-5, 15));
+		y3 = floor(frand(-5, 15));
+		mid = frand(-2, 2);
+
+		if (x3 - x1) {
+			x2 = floor(mid * (x3 - x1) + x1);
+			y2 = y1 + (x2 - x1) * (y3 - y1) / (x3 - x1);
+		} else if (y3 - y1) {
+			x2 = x1;
+			y2 = floor(mid * (y3 - y1) + y1);
+		} else {
+			x2 = x1;
+			y2 = y1;
+		}
+
+		if (!exact)
+			check_line(x1, y1, x2, y2, x3, y3, x2 == floor(x2) && y2 == floor(y2));
+		else if (x2 == floor(x2) && y2 == floor(y2))
+			check_line(x1, y1, x2, y2, x3, y3, 1);
+	}
+}
+
+
 int
 main(void)
 {
-	return 0; /* TODO add test */
+	raster = rtgrpblib_create_raster(10, 10);
+	ASSERT(raster);
+	refraster = rtgrpblib_create_raster(10, 10);
+	ASSERT(refraster);
+	srand((unsigned)(uintptr_t)raster);
+
+	check_simple_lines();
+	check_random_bounded_lines(1);
+	check_random_unbounded_lines(1);
+	//TODO check_random_bounded_lines(0);
+	//TODO check_random_unbounded_lines(0);
+	raster->draftness = -1;
+	refraster->draftness = -1;
+	check_simple_lines();
+	check_random_bounded_lines(1);
+	check_random_unbounded_lines(1);
+#ifdef TODO_NOT_IMPLEMENTED
+	check_random_bounded_lines(0);
+	check_random_unbounded_lines(0);
+#endif
+	raster->draftness = 0.5;
+	refraster->draftness = 0.5;
+
+	/* TODO add tests */
+
+	free(raster);
+	free(refraster);
+	return 0;
 }