First commit

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

1 files changed, 205 insertions, 0 deletions

diff --git a/rtgrpblib_draw_cubic_bezier.c b/rtgrpblib_draw_cubic_bezier.c
new file mode 100644
index 0000000..733e750
--- /dev/null
+++ b/rtgrpblib_draw_cubic_bezier.c
@@ -0,0 +1,205 @@
+/* See LICENSE file for copyright and license details. */
+#include "common.h"
+#ifndef TEST
+
+
+static size_t
+solve_cubic_bezier(double *restrict ts, double a, double b, double c, double d, double z)
+{
+	/*
+	 * B(t) = t(q - p) + p
+	 *   where
+	 *     p = t(v - u) + u
+	 *     q = t(w - v) + v
+	 *       where
+	 *          u = t(b - a) + a
+	 *          v = t(c - b) + b
+	 *          w = t(d - c) + c
+	 * 
+	 * B(t) = t³(d - 3c + 3b - a) + t²(3c - 6b + 3a) + t(3b - 3a) + a
+	 * 
+	 * z = t³(d - 3c + 3b - a) + t²(3c - 6b + 3a) + t(3b - 3a) + a
+	 * 
+	 * t³(d - 3c + 3b - a) + t²(3c - 6b + 3a) + t(3b - 3a) + a - z = 0
+	 *
+	 * A := d - 3c + 3b - a
+	 * B := 3c - 6b + 3a
+	 * C := 3b - 3a
+	 * D := a - z
+	 * 
+	 * t³A + t²B + tC + D = 0
+	 */
+	double c3 = 3.0 * c;
+	double b3 = 3.0 * b;
+	double a3 = 3.0 * a;
+	double A = (d - c3) + (b3 - a);
+	double B = (c3 - b3) - (b3 - a3);
+	double C = b3 - a3;
+	double D = a - z;
+
+	if (A)
+		return solve_cubic(ts, A, B, C, D);
+	else if (B)
+		return solve_quadratic(ts, B, C, D);
+	else if (C)
+		return solve_linear(ts, C, D);
+	else
+		return 0;
+}
+
+
+static double
+evaluate_cubic_bezier(double t, double a, double b, double c, double d)
+{
+	double u = fma(t, b - a, a);
+	double v = fma(t, c - b, b);
+	double w = fma(t, d - c, c);
+	double p = fma(t, v - u, u);
+	double q = fma(t, w - v, v);
+	return fma(t, q - p, p);
+}
+
+
+static double
+evaluate_cubic_bezier_derivative(double t, double a, double b, double c, double d)
+{
+	double u = fma(t, b - a, a);
+	double v = fma(t, c - b, b);
+	double w = fma(t, d - c, c);
+	return fma(t, (w - v) + (u - v), v - u);
+}
+
+
+static void
+draw_bounded_cubic_bezier(RASTER *restrict raster, double x1, double y1, double x2, double y2,
+                          double x3, double y3, double x4, double y4, double t1, double t2)
+{
+	double t = t1, x0, y0, dx, dy, x, y;
+	double w = (double)raster->width - 0.0001;
+	double h = (double)raster->height - 0.0001;
+
+	x0 = evaluate_cubic_bezier(t, x1, x2, x3, x4);
+	y0 = evaluate_cubic_bezier(t, y1, y2, y3, y4);
+	x0 = fmin(fmax(0, x0), w);
+	y0 = fmin(fmax(0, y0), h);
+
+	for (; t < t2; x0 = x, y0 = y) {
+		dx = evaluate_cubic_bezier_derivative(t, x1, x2, x3, x4);
+		dy = evaluate_cubic_bezier_derivative(t, y1, y2, y3, y4);
+
+		t += raster->draftness / hypot(dx, dy);
+		t = fmin(t, t2);
+		x = evaluate_cubic_bezier(t, x1, x2, x3, x4);
+		y = evaluate_cubic_bezier(t, y1, y2, y3, y4);
+		x = fmin(fmax(0, x), w);
+		y = fmin(fmax(0, y), h);
+
+		draw_bounded_line(raster, x0, y0, x, y);
+	}
+}
+
+
+void
+rtgrpblib_draw_cubic_bezier(RASTER *restrict raster, double x1, double y1, double x2, double y2,
+                            double x3, double y3, double x4, double y4)
+{
+	double x, y, yt1, yt2, dy;
+	double ts[2+4*3], t, t1, t2;
+	size_t nts = 0;
+	size_t i;
+
+#ifdef TODO /* untested */
+	/* Can we downgrade the curve to a quadratic bézier curve?
+	 * 
+	 * B(t) = (1-t)²Q₁ + 2(1-t)tQ₂ + t²Q₃
+	 * B(t) = (1-t)B(t) + tB(t) =
+	 *      = (1-t)((1-t)²Q₁ + 2(1-t)tQ₂ + t²Q₃) + t((1-t)²Q₁ + 2(1-t)tQ₂ + t²Q₃)
+	 *      = (1-t)³Q₁ + 2(1-t)²tQ₂ + (1-t)t²Q₃ + (1-t)²tQ₁ + 2(1-t)t²Q₂ + t³Q₃
+	 *      = (1-t)³Q₁ + 2(1-t)²tQ₂ + (1-t)²tQ₁ + 2(1-t)t²Q₂ + (1-t)t²Q₃ + t³Q₃
+	 *      = (1-t)³Q₁ + (1-t)²t(Q₁ + 2Q₂) + (1-t)t²(2Q₂ + Q₃) + t³Q₃
+	 *      = (1-t)³Q₁ + 3(1-t)²t(Q₁ + 2Q₂)/3 + 3(1-t)t²(2Q₂ + Q₃)/3 + t³Q₃
+	 * 
+	 * B(t) = (1-t)³C₁ + 3(1-t)²tC₂ + 3(1-t)t²C₃ + t³C₄
+	 * 
+	 * C₁ = Q₁, C₄ = Q₃
+	 * 
+	 * B(t) = (1-t)³Q₁ + 3(1-t)²t(Q₁ + 2Q₂)/3 + 3(1-t)t²(2Q₂ + Q₃)/3 + t³Q₃
+	 * B(t) = (1-t)³C₁ + 3(1-t)²tC₂ + 3(1-t)t²C₃ + t³C₄
+	 * 
+	 * (1-t)³Q₁ + 3(1-t)²t(Q₁ + 2Q₂)/3 + 3(1-t)t²(2Q₂ + Q₃)/3 + t³Q₃ = (1-t)³C₁ + 3(1-t)²tC₂ + 3(1-t)t²C₃ + t³C₄w
+	 * 3(1-t)²t(Q₁ + 2Q₂)/3 + 3(1-t)t²(2Q₂ + Q₃)/3 = 3(1-t)²tC₂ + 3(1-t)t²C₃
+	 * (1-t)²t(Q₁ + 2Q₂)/3 + (1-t)t²(2Q₂ + Q₃)/3 = (1-t)²tC₂ + (1-t)t²C₃
+	 * (1-t)t(Q₁ + 2Q₂)/3 + t²(2Q₂ + Q₃)/3 = (1-t)tC₂ + t²C₃
+	 * (1-t)(Q₁ + 2Q₂)/3 + t(2Q₂ + Q₃)/3 = (1-t)C₂ + tC₃
+	 * s(Q₁ + 2Q₂)/3 + t(2Q₂ + Q₃)/3 = sC₂ + tC₃
+	 * (Q₁ + 2Q₂)/3 = C₂, (2Q₂ + Q₃)/3 = C₃
+	 * Q₁ + 2Q₂ = 3C₂, 2Q₂ + Q₃ = 3C₃
+	 * C₁ + 2Q₂ = 3C₂, 2Q₂ + C₄ = 3C₃
+	 * 2Q₂ = 3C₂ - C₁, 2Q₂ = 3C₃ - C₄
+	 * 2Q₂ = 3C₂ - C₁ = 3C₃ - C₄
+	 */
+	if (x2 + x2 + x2 - x1 == x3 + x3 + x3 - x4 && y2 + y2 + y2 - y1 == y3 + y3 + y3 - y4) {
+		rtgrpblib_draw_quadratic_bezier(raster, x1, y1, x2 + (x2 - x1) / 2.0, y2 + (y2 - y1) / 2.0, x4, y4);
+		return;
+	}
+#endif
+
+	/* Beginning and end of curve */
+	ts[nts++] = 0.0;
+	ts[nts++] = 1.0;
+	/* Get points where the end intersects one of the edges of the raster */
+	nts += solve_cubic_bezier(&ts[nts], y1, y2, y3, y4, 0);
+	nts += solve_cubic_bezier(&ts[nts], y1, y2, y3, y4, (double)raster->height);
+	nts += solve_cubic_bezier(&ts[nts], x1, x2, x3, x4, 0);
+	nts += solve_cubic_bezier(&ts[nts], x1, x2, x3, x4, (double)raster->width);
+	/* Sort all points of interest */
+	qsort(ts, nts, sizeof(*ts), doublepcmp);
+
+	for (i = 0; i < nts - 1; i++) {
+		t1 = ts[i];
+		t2 = ts[i + 1];
+		if (t1 == t2)
+			continue;
+		t = 0.5 * (t1 + t2);
+
+		/* Remove any segments above or below the raster */
+		y = evaluate_cubic_bezier(t, y1, y2, y3, y4);
+		if (y < 0 || y >= (double)raster->height)
+			continue;
+
+		/* If the segment is inside the raster, draw it, */
+		x = evaluate_cubic_bezier(t, x1, x2, x3, x4);
+		if (0 <= x && x < (double)raster->width) {
+			draw_bounded_cubic_bezier(raster, x1, y1, x2, y2, x3, y3, x4, y4, t1, t2);
+			continue;
+		}
+
+		/* otherwise draw it flattened to the edge; however,
+		 * we will reduce it to a from the start to the end
+		 * and draw its project on the edge, which works
+		 * because if the curve bends it will cancel it self
+		 * out except within this project; this project will
+		 * also preserv the direction of part of the curve
+		 * that is not cancelled out */
+		yt1 = evaluate_cubic_bezier(t1, y1, y2, y3, y4);
+		yt2 = evaluate_cubic_bezier(t2, y1, y2, y3, y4);
+		dy = yt2 - yt1;
+		if (x < 0)
+			draw_vertical_line(raster, 0, yt1, yt2, SIGNUM(dy));
+		else
+			draw_vertical_line_opposite_only(raster, raster->width - 1, yt1, yt2, SIGNUM(dy));
+	}
+}
+
+
+#else
+
+
+int
+main(void)
+{
+	return 0; /* TODO add test */
+}
+
+
+#endif