First commit

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

1 files changed, 166 insertions, 0 deletions

diff --git a/rtgrpblib_draw_quadratic_bezier.c b/rtgrpblib_draw_quadratic_bezier.c
new file mode 100644
index 0000000..dd4a484
--- /dev/null
+++ b/rtgrpblib_draw_quadratic_bezier.c
@@ -0,0 +1,166 @@
+/* See LICENSE file for copyright and license details. */
+#include "common.h"
+#ifndef TEST
+
+
+static size_t
+solve_quadratic_bezier(double *restrict ts, double a, double b, double c, double z)
+{
+	/*
+	 * B(t) = t(q - p) + p
+	 *   where
+	 *     p = t(b - a) + a
+	 *     q = t(c - b) + b
+	 * 
+	 * B(t) = t²(c - 2b + a) + t(2b) + a
+	 *
+	 * z = t²(c - 2b + a) + t(2b) + a
+	 *
+	 * t²(c - 2b + a) + t(2b) + a - z = 0
+	 * 
+	 * A := c - 2b + a
+	 * B := 2b
+	 * C := a - z
+	 * 
+	 * t²A + tB + C = 0
+	 */
+
+	double A = (c - a) + (b + b);
+	double B = b + b;
+	double C = a - z;
+
+	if (A)
+		return solve_quadratic(ts, A, B, C);
+	else if (B)
+		return solve_linear(ts, B, C);
+	else
+		return 0;
+}
+
+
+static double
+evaluate_quadratic_bezier(double t, double a, double b, double c)
+{
+	double p = fma(t, b - a, a);
+	double q = fma(t, c - b, b);
+	return fma(t, q - p, p);
+}
+
+
+static void
+draw_bounded_quadratic_bezier(RASTER *restrict raster, double x1, double y1, double x2, double y2,
+                              double x3, double y3, 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;
+
+	double kx = (x3 - x2) + (x1 - x2), mx = x2 - x1;
+	double ky = (y3 - y2) + (y1 - y2), my = y2 - y1;
+
+	x0 = evaluate_quadratic_bezier(t, x1, x2, x3);
+	y0 = evaluate_quadratic_bezier(t, y1, y2, y3);
+	x0 = fmin(fmax(0, x0), w);
+	y0 = fmin(fmax(0, y0), h);
+
+	for (; t < t2; x0 = x, y0 = y) {
+		dx = fma(t, kx, mx);
+		dy = fma(t, ky, my);
+
+		t += raster->draftness / hypot(dx, dy);
+		t = fmin(t, t2);
+		x = evaluate_quadratic_bezier(t, x1, x2, x3);
+		y = evaluate_quadratic_bezier(t, y1, y2, y3);
+
+		x = fmin(fmax(0, x), w);
+		y = fmin(fmax(0, y), h);
+
+		draw_bounded_line(raster, x0, y0, x, y);
+	}
+}
+
+
+void
+rtgrpblib_draw_quadratic_bezier(RASTER *restrict raster, double x1, double y1,
+                                double x2, double y2, double x3, double y3)
+{
+	double x, y, yt1, yt2, dy;
+	double ts[2+2*4], t, t1, t2;
+	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);
+		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_quadratic_bezier(&ts[nts], y1, y2, y3, (double)raster->height);
+	nts += solve_quadratic_bezier(&ts[nts], x1, x2, x3, (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_quadratic_bezier(t, y1, y2, y3);
+		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) {
+			draw_bounded_quadratic_bezier(raster, x1, y1, x2, y2, x3, y3, 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_quadratic_bezier(t1, y1, y2, y3);
+		yt2 = evaluate_quadratic_bezier(t2, y1, y2, y3);
+		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