First commit

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

20 files changed, 1893 insertions, 0 deletions

diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..d8a83f0
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,15 @@
+*\#*
+*~
+*.o
+*.a
+*.lo
+*.su
+*.so
+*.so.*
+*.dll
+*.dylib
+*.gch
+*.gcov
+*.gcno
+*.gcda
+*.test
diff --git a/LICENSE b/LICENSE
new file mode 100644
index 0000000..0be2ccf
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,15 @@
+ISC License
+
+© 2023 Mattias Andrée <maandree@kth.se>
+
+Permission to use, copy, modify, and/or distribute this software for any
+purpose with or without fee is hereby granted, provided that the above
+copyright notice and this permission notice appear in all copies.
+
+THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
+OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
diff --git a/Makefile b/Makefile
new file mode 100644
index 0000000..1db2138
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,88 @@
+.POSIX:
+
+CONFIGFILE = config.mk
+include $(CONFIGFILE)
+
+OS = linux
+# Linux:   linux
+# Mac OS:  macos
+# Windows: windows
+include mk/$(OS).mk
+
+
+LIB_MAJOR = 1
+LIB_MINOR = 0
+LIB_VERSION = $(LIB_MAJOR).$(LIB_MINOR)
+LIB_NAME = rifunktionsteckensnittsglyfrasteriseringsprogrambiblioteket
+
+
+OBJ =\
+	equations.o\
+	lines.o\
+	sorting.o\
+	rtgrpblib_create_raster.o\
+	rtgrpblib_draw_linear_bezier.o\
+	rtgrpblib_draw_quadratic_bezier.o\
+	rtgrpblib_draw_cubic_bezier.o\
+	rtgrpblib_fill_shapes.o\
+	rtgrpblib_reset_raster.o\
+	rtgrpblib_set_draftness.o
+
+HDR =\
+	lib$(LIB_NAME).h\
+	common.h
+
+LOBJ = $(OBJ:.o=.lo)
+TESTS = $(OBJ:.o=.test)
+
+
+all: lib$(LIB_NAME).a lib$(LIB_NAME).$(LIBEXT) $(TESTS)
+$(OBJ): $(HDR)
+$(LOBJ): $(HDR)
+$(TESTS): $(HDR) lib$(LIB_NAME).a
+
+.c.o:
+	$(CC) -c -o $@ $< $(CFLAGS) $(CPPFLAGS)
+
+.c.lo:
+	$(CC) -fPIC -c -o $@ $< $(CFLAGS) $(CPPFLAGS)
+
+.c.test:
+	$(CC) -o $@ $< lib$(LIB_NAME).a $(CFLAGS) $(CPPFLAGS) -DTEST $(LDFLAGS)
+
+lib$(LIB_NAME).a: $(OBJ)
+	@rm -f -- $@
+	$(AR) rc $@ $(OBJ)
+	$(AR) ts $@ > /dev/null
+
+lib$(LIB_NAME).$(LIBEXT): $(LOBJ)
+	$(CC) $(LIBFLAGS) -o $@ $(LOBJ) $(LDFLAGS)
+
+check: $(TESTS)
+	@for t in $(TESTS); do printf './%s\n' $$t; ./$$t || exit 1; done
+
+install: lib$(LIB_NAME).a lib$(LIB_NAME).$(LIBEXT)
+	mkdir -p -- "$(DESTDIR)$(PREFIX)/lib"
+	mkdir -p -- "$(DESTDIR)$(PREFIX)/include"
+	cp -- lib$(LIB_NAME).a "$(DESTDIR)$(PREFIX)/lib/"
+	cp -- lib$(LIB_NAME).$(LIBEXT) "$(DESTDIR)$(PREFIX)/lib/lib$(LIB_NAME).$(LIBMINOREXT)"
+	$(FIX_INSTALL_NAME) "$(DESTDIR)$(PREFIX)/lib/lib$(LIB_NAME).$(LIBMINOREXT)"
+	ln -sf -- lib$(LIB_NAME).$(LIBMINOREXT) "$(DESTDIR)$(PREFIX)/lib/lib$(LIB_NAME).$(LIBMAJOREXT)"
+	ln -sf -- lib$(LIB_NAME).$(LIBMAJOREXT) "$(DESTDIR)$(PREFIX)/lib/lib$(LIB_NAME).$(LIBEXT)"
+	cp -- lib$(LIB_NAME).h "$(DESTDIR)$(PREFIX)/include/"
+
+uninstall:
+	-rm -f -- "$(DESTDIR)$(PREFIX)/lib/lib$(LIB_NAME).a"
+	-rm -f -- "$(DESTDIR)$(PREFIX)/lib/lib$(LIB_NAME).$(LIBMAJOREXT)"
+	-rm -f -- "$(DESTDIR)$(PREFIX)/lib/lib$(LIB_NAME).$(LIBMINOREXT)"
+	-rm -f -- "$(DESTDIR)$(PREFIX)/lib/lib$(LIB_NAME).$(LIBEXT)"
+	-rm -f -- "$(DESTDIR)$(PREFIX)/include/lib$(LIB_NAME).h"
+
+clean:
+	-rm -f -- *.o *.a *.lo *.su *.so *.so.* *.dll *.dylib *.test
+	-rm -f -- *.gch *.gcov *.gcno *.gcda *.$(LIBEXT)
+
+.SUFFIXES:
+.SUFFIXES: .lo .o .c .test
+
+.PHONY: all check install uninstall clean
diff --git a/common.h b/common.h
new file mode 100644
index 0000000..ae9a1d7
--- /dev/null
+++ b/common.h
@@ -0,0 +1,147 @@
+/* See LICENSE file for copyright and license details. */
+#include "librifunktionsteckensnittsglyfrasteriseringsprogrambiblioteket.h"
+#include <errno.h>
+#include <math.h>
+#include <stddef.h>
+#include <stdio.h>
+#include <stdint.h>
+#include <stdlib.h>
+#include <string.h>
+
+
+#if defined(__GNUC__) && !defined(__clang__)
+# pragma GCC diagnostic ignored "-Wunsuffixed-float-constants"
+# pragma GCC diagnostic ignored "-Wempty-body"
+# pragma GCC diagnostic ignored "-Wfloat-equal"
+# define PURE_FUNCTION __attribute__((__pure__))
+# define CONST_FUNCTION __attribute__((__const__))
+#elif defined(__clang__)
+# pragma clang diagnostic ignored "-Wcomma"
+# pragma clang diagnostic ignored "-Wfloat-equal"
+# pragma clang diagnostic ignored "-Wfloat-conversion"
+# define PURE_FUNCTION
+# define CONST_FUNCTION
+#else
+# define PURE_FUNCTION
+# define CONST_FUNCTION
+#endif
+
+
+#define SIGNUM(X) (((X) > 0) - ((X) < 0))
+
+#define draw_vertical_line_opposite_only(...)
+
+
+typedef struct rtgrpblib_cell {
+	/**
+	 * The sum all contribution lines make to
+	 * the cell's ink level
+	 * 
+	 * For each line, this is calculated by
+	 * the signed area of the right triangle,
+	 * whose catheti are parallel to the edges
+	 * and whose hypotenuse is the line's segment
+	 * in the cell, plus the area of the rectangle
+	 * whose left edge is the opposite cathetus
+	 * of the right triangle, and whose right edge
+	 * intersects with the right cell edge
+	 */
+	double cell_coverage; /* ink here */
+
+	/**
+	 * The sum all contribution lines make to
+	 * the shadow on the cell's right edge:
+	 * the ink level on the next cell
+	 * 
+	 * For each line, this is calculated by the
+	 * signed size of the opposite (right-side)
+	 * cathetus of the right triangle, whose
+	 * catheti are parallel to the edges and
+	 * whose hypotenuse is the line's segment
+	 * in the cell
+	 */
+	double opposite_coverage; /* ink there */
+} CELL;
+
+typedef struct rtgrpblib_raster {
+	size_t height;
+	size_t width;
+	double draftness; /* TODO support 0 */
+	CELL cells[];
+} RASTER;
+
+
+
+#define TOLERANCE 0.000001
+
+CONST_FUNCTION
+static inline int
+iszeroish(double x)
+{
+	return x <= TOLERANCE && x >= -TOLERANCE;
+}
+
+
+
+/* equations.c */
+#define solve_cubic rtgrpblib_solve_cubic__
+#define solve_quadratic rtgrpblib_solve_quadratic__
+#define solve_linear rtgrpblib_solve_linear__
+size_t solve_cubic(double *ts, double a, double b, double c, double d);
+size_t solve_quadratic(double *ts, double a, double b, double c);
+size_t solve_linear(double *ts, double a, double b);
+
+
+/* lines.c */
+#ifndef draw_vertical_line_opposite_only
+# define draw_vertical_line_opposite_only rtgrpblib_draw_vertical_line_opposite_only__
+void draw_vertical_line_opposite_only(RASTER *restrict raster, size_t x, double y1, double y2, int ydir);
+#else
+# define ROWWISE_RESET_INKLEVEL
+#endif
+#define draw_vertical_line rtgrpblib_draw_vertical_line__
+#define draw_diagonal_line rtgrpblib_draw_diagonal_line__
+#define draw_bounded_line rtgrpblib_draw_bounded_line__
+void draw_vertical_line(RASTER *restrict raster, double x1, double y1, double y2, int ydir);
+void draw_diagonal_line(RASTER *restrict raster, double x1, double y1, double x2, double y2,
+                        double dx, double dy, int xdir, int ydir);
+void draw_bounded_line(RASTER *restrict raster, double x1, double y1, double x2, double y2);
+
+
+/* sorting.c */
+#define doublepcmp rtgrpblib_doublepcmp__
+PURE_FUNCTION int doublepcmp(const void *avp, const void *bvp);
+
+
+
+#ifdef TEST
+
+#if defined(__GNUC__) && !defined(__clang__)
+# pragma GCC diagnostic ignored "-Wunsuffixed-float-constants"
+# pragma GCC diagnostic ignored "-Wfloat-equal"
+# pragma GCC diagnostic ignored "-Wsuggest-attribute=const"
+#elif defined(__clang__)
+# pragma clang diagnostic ignored "-Wfloat-equal"
+#endif
+
+static inline int
+eq(double a, double b)
+{
+	return iszeroish(a - b);
+}
+
+static inline int
+tolerant_eq(double a, double b)
+{
+	return iszeroish((a - b) / 100.0);
+}
+
+# define ASSERT(ASSERTION)\
+	do {\
+		if (!(ASSERTION)) {\
+			fprintf(stderr, "Failed assertion at line %i: %s\n", __LINE__, #ASSERTION);\
+			exit(1);\
+		}\
+	} while (0)
+
+#endif
diff --git a/config.mk b/config.mk
new file mode 100644
index 0000000..eb89423
--- /dev/null
+++ b/config.mk
@@ -0,0 +1,10 @@
+PREFIX    = /usr
+MANPREFIX = $(PREFIX)/share/man
+
+CC = cc
+
+CPPFLAGS = -D_DEFAULT_SOURCE -D_BSD_SOURCE -D_XOPEN_SOURCE=700 -D_GNU_SOURCE
+CFLAGS   = -std=c99 -Wall
+LDFLAGS  = -lm
+
+# Do NOT use -ffast-math
diff --git a/demo.c b/demo.c
new file mode 100644
index 0000000..76473e7
--- /dev/null
+++ b/demo.c
@@ -0,0 +1,91 @@
+/* See LICENSE file for copyright and license details. */
+#include "common.h"
+
+
+int
+main(void)
+{
+#define P(I) (points[I].x - 0), (points[I].y - 0)
+	/* TODO test drawing outside of raster */
+
+	RTGRPBLIB_RASTER *raster;
+	double *image;
+	size_t x, y, i, j, n;
+	size_t width = 40;
+	size_t height = 30;
+	struct {
+		double x, y;
+	} points[] = {
+#if 0
+		{20, 5},
+		{30, 15},
+		{20, 25},
+		{10, 15}
+#elif 0
+		{20, 5},
+		{35, 15},
+		{20, 25},
+		{5, 15}
+#else
+		{5.25, 5.25},
+		{35.75, 5.25},
+		{35.75, 25.75},
+		{5.25, 25.75}
+#endif
+	};
+
+	raster = rtgrpblib_create_raster(width, height);
+	image  = calloc(width * height, sizeof(*image));
+
+#if 0
+	n = sizeof(points) / sizeof(*points);
+	for (i = 0; i < n; i++) {
+		j = (i + 1) % n;
+		rtgrpblib_draw_linear_bezier(raster, P(i), P(j));
+	}
+#elif 0
+	rtgrpblib_draw_quadratic_bezier(raster, P(2), P(3), P(0));
+	rtgrpblib_draw_linear_bezier(raster, P(0), P(2));
+#else
+	rtgrpblib_draw_cubic_bezier(raster, P(0), P(1), P(2), P(3));
+	rtgrpblib_draw_linear_bezier(raster, P(3), P(0));
+#endif
+
+	rtgrpblib_fill_shapes(image, width, raster);
+
+#if 1
+	printf("\nOutline (area)\n");
+	for (y = 0; y < height; y++) {
+		for (x = 0; x < width; x++)
+			printf(raster->cells[y * width + x].cell_coverage ? "%+.1lf " : " 0   ",
+			       raster->cells[y * width + x].cell_coverage);
+		printf("\n");
+	}
+	printf("\nOutline (shadow)\n");
+	for (y = 0; y < height; y++) {
+		for (x = 0; x < width; x++)
+			printf(raster->cells[y * width + x].opposite_coverage ? "%+.1lf " : " 0   ",
+			       raster->cells[y * width + x].opposite_coverage);
+		printf("\n");
+	}
+	printf("\nFill\n");
+	for (y = 0; y < height; y++) {
+		for (x = 0; x < width; x++) {
+			unsigned val = (unsigned)(image[y * width + x] * 255 + 0.5);
+			printf("%02X ", val - (val == 128));
+		}
+		printf("\n");
+	}
+#else
+	printf("P2\n# %s\n%zu %zu\n255\n", "Note: no gamma correction applied", width, height);
+	for (y = 0; y < height; y++) {
+		for (x = 0; x < width; x++)
+			printf("%i ", (int)(image->cells[y * width + x] * 255 + 0.5));
+		printf("\n");
+	}
+#endif
+
+	free(raster);
+	free(image);
+	return 0;
+}
diff --git a/equations.c b/equations.c
new file mode 100644
index 0000000..571fc26
--- /dev/null
+++ b/equations.c
@@ -0,0 +1,609 @@
+/* See LICENSE file for copyright and license details. */
+#include "common.h"
+#ifndef TEST
+
+
+CONST_FUNCTION
+static double
+refine_cubic_solution(double t, double a, double b, double c, double d)
+{
+	size_t i;
+	double x;
+	for (i = 0; x = t*t*t*a + t*t*b + t*c + d, !iszeroish(x) && i < 10; i++)
+		t -= x / (t*t*(a*3) + t*(b*2) + c);
+	return t;
+}
+
+
+static size_t
+solve_cubic_with_1_known_root(double *ts, double a, double b, double t1)
+{
+	/*
+	 * (t³ + t²a + tb + c) / (t - t₁) = t² + tp + q
+	 * t³ + t²a + tb + c = (t² + tp + q)(t - t₁)
+	 * t³ + t²a + tb + c = t²(t - t₁) + tp(t - t₁) + q(t - t₁)
+	 * t³ + t²a + tb + c = t³ - t²t₁ + t²p - tpt₁ + tq - qt₁
+	 * t³ + t²a + tb + c = t³ + t²(p - t₁) + t(q - pt₁) - qt₁
+	 *
+	 * ⎧a = p - t₁
+	 * ⎨b = q - pt₁
+	 * ⎩c = -qt₁
+	 *
+	 * ⎧a + t₁ = p
+	 * ⎨b + pt₁ = q
+	 * ⎩c = -qt₁
+	 *
+	 * ⎰a + t₁ = p
+	 * ⎱b + pt₁ = q
+	 *
+	 * ⎰p = a + t₁
+	 * ⎱q = b + pt₁
+	 */
+	double p = a + t1;
+	double q = b + p * t1;
+
+	return solve_quadratic(ts, 1, p, q);
+}
+
+
+static size_t
+solve_cubic_with_2_known_roots(double *ts, double a, double b, double c, double t1, double t2)
+{
+	/*
+	 * (t - t₁)(t - t₂)(t - t₃)
+	 * (t - t₁)t(t - t₃) - (t - t₁)t₂(t - t₃)
+	 * (t - t₁)(t² - t₃t) - (t - t₁)(t₂t - t₂t₃)
+	 * (t - t₁)t² - (t - t₁)t₃t - (t - t₁)t₂t + (t - t₁)t₂t₃
+	 * t³ - t₁t² - t₃t² + t₁t₃t - t₂t² + t₁t₂t + t₂t₃t - t₁t₂t₃
+	 * t³ - t²(t₁ + t₂ + t₃) + t(t₁t₂ + t₁t₃ + t₂t₃) - t₁t₂t₃
+	 * 
+	 * ⎧a = -t₁ - t₂ - t₃
+	 * ⎨b = t₁t₂ + t₁t₃ + t₂t₃
+	 * ⎩c = -t₁t₂t₃
+	 * 
+	 * ⎧t₃ = -t₁ - t₂ - a
+	 * ⎨t₃ = (b - t₁t₂)/(t₁ + t₂)
+	 * ⎩t₃ = -c/t₁t₂
+	 */
+	double t12sum  = t1 + t2;
+	double t12prod = t1 * t2;
+	double t3a = -a - t12sum;
+	double t3b = iszeroish(t12sum) ? t3a : (b - t12prod) / t12sum;
+	double t3c = iszeroish(t12prod) ? t3a : -c / t12prod;
+	if (iszeroish(t3a - t3b) && iszeroish(t3a - t3c)) {
+		ts[0] = (t3a + t3b + t3c) / 3.0;
+		return 0.0 <= ts[0] && ts[0] <= 1.0;
+	}
+	t1 = 0.5 * (t1 + t2);
+	t1 = refine_cubic_solution(t1, 1.0, a, b, c);
+	return solve_cubic_with_1_known_root(ts, a, b, t1);
+}
+
+
+size_t
+solve_cubic(double *ts, double a, double b, double c, double d)
+{
+	double t1, t2, t3;
+	unsigned int solutions = 0;
+	size_t n = 0;
+
+	double u = 1 / (a * -3);
+
+	double D0 = b*b - 3*(a*c);
+	double D1 = b*b*b - (4.5*b)*(a*c) + (13.5*a)*(a*d);
+
+	double kk = D1*D1 - D0*D0*D0;
+
+	double r, i;
+
+	double Kr;
+	if (kk >= 0) {
+		Kr = D1 < 0 ? cbrt(D1 + sqrt(kk)) : cbrt(D1 - sqrt(kk));
+	} else {
+		double ki = sqrt(-kk);
+		double h = pow(D1*D1 + ki*ki, 1.0 / 6.0);
+		double v = (0.1 / 0.3) * atan2(ki, D1);
+
+		Kr = h * cos(v);
+
+		if (1 /* || fabs(cos(v)) > 1.0 - TOLERANCE */) {
+			double Ki = h * sin(v);
+			double t = 1.0 / (Kr * Kr + Ki * Ki);
+			double m = t * D0 + 1;
+			if (iszeroish(Ki - D0 * Ki * t)) {
+				ts[n] = u * (b + Kr + D0 * Kr * t);
+				ts[n] = t1 = refine_cubic_solution(ts[n], a, b, c, d);
+				solutions |= 1;
+				n += 0.0 <= ts[n] && ts[n] <= 1.0;
+			}
+			if (iszeroish(Ki - sqrt(3.0) * Kr)) {
+				ts[n] = u * (b - (0.5 * Kr + sqrt(0.75) * Ki) * m);
+				ts[n] = t2 = refine_cubic_solution(ts[n], a, b, c, d);
+				solutions |= 2;
+				n += 0.0 <= ts[n] && ts[n] <= 1.0;
+			}
+			if (iszeroish(Ki + sqrt(3.0) * Kr)) {
+				ts[n] = u * (b - (0.5 * Kr - sqrt(0.75) * Ki) * m);
+				ts[n] = t3 = refine_cubic_solution(ts[n], a, b, c, d);
+				solutions |= 4;
+				n += 0.0 <= ts[n] && ts[n] <= 1.0;
+			}
+
+			/* Sometimes the maths doesn't work well enough */
+			if (solutions != 7) {
+#if defined(__GNUC__) && !defined(__clang__)
+# pragma GCC diagnostic push
+# pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
+#elif defined(__clang__)
+# pragma clang diagnostic push
+# pragma clang diagnostic ignored "-Wconditional-uninitialized"
+#endif
+				switch (solutions) {
+				case 1:
+					n += solve_cubic_with_1_known_root(&ts[n], b / a, c / a, t1);
+					break;
+				case 2:
+					n += solve_cubic_with_1_known_root(&ts[n], b / a, c / a, t2);
+					break;
+				case 4:
+					n += solve_cubic_with_1_known_root(&ts[n], b / a, c / a, t3);
+					break;
+				case 7 - 1:
+					n += solve_cubic_with_2_known_roots(&ts[n], b / a, c / a, d / a, t2, t3);
+					break;
+				case 7 - 2:
+					n += solve_cubic_with_2_known_roots(&ts[n], b / a, c / a, d / a, t1, t3);
+					break;
+				case 7 - 4:
+					n += solve_cubic_with_2_known_roots(&ts[n], b / a, c / a, d / a, t1, t2);
+					break;
+				}
+#if defined(__GNUC__) && !defined(__clang__)
+# pragma GCC diagnostic pop
+#elif defined(__clang__)
+# pragma clang diagnostic pop
+#endif
+			}
+			return n;
+		}
+	}
+
+	if (!iszeroish(Kr)) {
+		r = Kr + D0 / Kr;
+		i = u * sqrt(0.75) * (Kr - D0 / Kr);
+
+	real_case:
+		ts[n] = u * (b + r);
+		ts[n] = t1 = refine_cubic_solution(ts[n], a, b, c, d);
+		n += 0.0 <= ts[n] && ts[n] <= 1.0;
+		if (iszeroish(i)) {
+			ts[n] = u * (b - 0.5 * r);
+			ts[n] = refine_cubic_solution(ts[n], a, b, c, d);
+			n += 0.0 <= ts[n] && ts[n] <= 1.0;
+			return n;
+		} else {
+			n += solve_cubic_with_1_known_root(&ts[n], b / a, c / a, t1);
+			return n;
+		}
+	} else {
+		double D0divKr = D0 / Kr;
+		r = Kr + D0divKr;
+		i = u * sqrt(0.75) * (Kr - D0divKr);
+		if (isfinite(D0divKr)) /* NOT -ffast-math compatible */
+			goto real_case;
+
+		ts[n] = u * b;
+		ts[n] = t1 = refine_cubic_solution(ts[n], a, b, c, d);
+		n += 0.0 <= ts[n] && ts[n] <= 1.0;
+		n += solve_cubic_with_1_known_root(&ts[n], b / a, c / a, t1);
+		return n;
+	}
+}
+
+
+size_t
+solve_quadratic(double *ts, double a, double b, double c)
+{
+	size_t n = 0;
+
+	double negahalf_p = -0.5 * (b / a);
+	double q = c / a;
+
+	double D = negahalf_p * negahalf_p - q;
+
+	if (D < 0) {
+		if (D < -TOLERANCE)
+			return 0;
+
+		ts[n] = negahalf_p;
+		n += (0.0 <= ts[n] && ts[n] <= 1.0);
+
+		return n;
+	} else {
+		ts[n] = negahalf_p + sqrt(D);
+		n += (0.0 <= ts[n] && ts[n] <= 1.0);
+
+		ts[n] = negahalf_p - sqrt(D);
+		n += (0.0 <= ts[n] && ts[n] <= 1.0);
+
+		return n;
+	}
+}
+
+
+size_t
+solve_linear(double *ts, double a, double b)
+{
+	ts[0] = -b / a;
+	return (0.0 <= ts[0] && ts[0] <= 1.0);
+}
+
+
+#else
+
+
+#define Re(X) X, 0
+#define Im(X) 0, X
+#define Cx(R, I) R, I
+
+static void
+check_cubic(double r1, double i1, double r2, double i2, double r3, double i3, double k, int strict_real)
+{
+	size_t nexpected = 0, ngot, i, j;
+	double expected[3], got[3];
+	double error, maxerror = 0;
+	double ra, rb, rc;
+	double ia, ib, ic;
+	int found[3], correct;
+
+	if (iszeroish(i1))
+		expected[nexpected++] = r1;
+	if (iszeroish(i2))
+		expected[nexpected++] = r2;
+	if (iszeroish(i3))
+		expected[nexpected++] = r3;
+
+	/*
+	 * k(t₁ - t)(t₂ - t)(t₃ - t) = 0
+	 * k(t₁ - t)t₂(t₃ - t) - k(t₁ - t)t(t₃ - t) = 0
+	 * kt₂(t₁ - t)(t₃ - t) - kt(t₁ - t)(t₃ - t) = 0
+	 * kt₂(t₁(t₃ - t) - t(t₃ - t)) - kt(t₁(t₃ - t) - t(t₃ - t)) = 0
+	 * kt₂(t₁t₃ - t₁t - tt₃ + t²) - kt(t₁t₃ - t₁t - tt₃ + t²) = 0
+	 * kt₂(t₁t₃ - t₁t - tt₃ + t²) + kt(-t₁t₃ + t₁t + tt₃ - t²) = 0
+	 * kt₂t₁t₃ - kt₂t₁t - kt₂tt₃ + kt₂t² - ktt₁t₃ + ktt₁t + kt²t₃ - kt³ = 0
+	 * kt₂t₁t₃ - kt₂t₁t¹ - kt₂t¹t₃ + kt₂t² - kt¹t₁t₃ + kt²t₁ + kt²t₃ - kt³ = 0
+	 * -kt³ + kt₂t² + kt²t₁ + kt²t₃ - kt₂t₁t¹ - kt₂t¹t₃ - kt¹t₁t₃ + kt₂t₁t₃ = 0
+	 * kt³ - kt₂t² - kt²t₁ - kt²t₃ + kt₂t₁t¹ + kt₂t¹t₃ + kt¹t₁t₃ - kt₂t₁t₃ = 0
+	 * t³k - t²k(t₁ + t₂ + t₃) + t¹k(t₂t₁ + t₂t₃ + t₁t₃) - kt₂t₁t₃ = 0
+	 * ⎧a := t₁ + t₂ + t₃      ⎫
+	 * ⎨b := t₁t₂ + t₁t₃ + t₂t₃⎬
+	 * ⎩c := t₁t₂t₃            ⎭
+	 * t³k - t²ka + t¹kb - kc = 0
+	 */
+
+	ra = r1 + r2 + r3;
+	ia = i1 + i2 + i3;
+	rb = r1*r2 - i1*i2 + r1*r3 - i1*i3 + r2*r3 - i2*i3;
+	ib = r1*i2 + i1*r2 + r1*i3 + i1*r3 + r2*i3 + i2*r3;
+	rc = r1*r2*r3 - r1*i2*i3 - i1*r2*i3 - i1*i2*r3;
+	ic = r1*r2*i3 + r1*i2*r3 + i1*r2*r3 - i1*i2*i3;
+
+	if (strict_real) {
+		ASSERT(iszeroish(ia));
+		ASSERT(iszeroish(ib));
+		ASSERT(iszeroish(ic));
+	} else {
+		if (!iszeroish(ia) || !iszeroish(ib) || !iszeroish(ic))
+			return;
+	}
+
+	ngot = solve_cubic(got, k, -k*ra, k*rb, -k*rc);
+
+	for (i = 0; i < sizeof(found) / sizeof(*found); i++)
+		found[i] = (i >= nexpected || expected[i] < 0.0 + TOLERANCE || expected[i] > 1.0 - TOLERANCE);
+	if (!nexpected) {
+		ASSERT(!ngot);
+		return;
+	}
+	for (i = 0; i < ngot; i++) {
+		correct = (got[i] < 0.0 + TOLERANCE || got[i] > 1.0 - TOLERANCE);
+		error = fabs(got[i] - expected[0]);
+		for (j = 0; j < nexpected; j++) {
+			error = fmin(error, fabs(got[i] - expected[j]));
+			if (tolerant_eq(got[i], expected[j]))
+				found[j] = correct = 1;
+		}
+		if (!correct)
+			maxerror = fmax(maxerror, error);
+	}
+	if (!found[0] || !found[1] || !found[2] || !tolerant_eq(maxerror, 0)) {
+		fprintf(stderr, "\n%lg, %lg, %lg, %lg, %lg, %lg, %lg\n", r1, i1, r2, i2, r3, i3, k);
+		fprintf(stderr, "\t(%la, %la, %la, %la, %la, %la, %la)\n", r1, i1, r2, i2, r3, i3, k);
+		fprintf(stderr, "\tmaxerror: %lf (%la)\n", maxerror, maxerror);
+		fprintf(stderr, "\tfound: %i, %i, %i\n", found[0], found[1], found[2]);
+		fprintf(stderr, "\tgot: %zu\n", ngot);
+		for (i = 0; i < ngot; i++)
+			fprintf(stderr, "\t\t%lf (%la)\n", got[i], got[i]);
+		fprintf(stderr, "\texpected: %zu\n", nexpected);
+		for (i = 0; i < nexpected; i++)
+			fprintf(stderr, "\t\t%lf (%la)\n", expected[i], expected[i]);
+	}
+	ASSERT(found[0]);
+	ASSERT(found[1]);
+	ASSERT(found[2]);
+	ASSERT(tolerant_eq(maxerror, 0));
+}
+
+
+static void
+check_quadratic(double r1, double i1, double r2, double i2, double k, int strict_real)
+{
+	size_t nexpected = 0, ngot, i, j;
+	double expected[2], got[2];
+	double error, maxerror = 0;
+	double ra, rb;
+	double ia, ib;
+	int found[2], correct;
+
+	if (iszeroish(i1))
+		expected[nexpected++] = r1;
+	if (iszeroish(i2))
+		expected[nexpected++] = r2;
+
+	/*
+	 * k(t₁ - t)(t₂ - t) = 0
+	 * k(t₁(t₂ - t) - t(t₂ - t)) = 0
+	 * k(t₁t₂ - t¹t₁ - t¹t₂ + t²) = 0
+	 * kt₁t₂ - t¹kt₁ - t¹kt₂ + kt² = 0
+	 * t²k - t¹kt₁ - t¹kt₂ + kt₁t₂ = 0
+	 * t²k - t¹k(t₁ + kt₂) + kt₁t₂ = 0
+	 * ⎰a := t₁ + t₂⎱
+	 * ⎱b := t₁t₂   ⎰
+	 * t²k - t¹ka + kb = 0
+	 */
+
+	ra = r1 + r2;
+	ia = i1 + i2;
+	rb = r1*r2 - i1*i2;
+	ib = r1*i2 + i1*r2;
+
+	if (strict_real) {
+		ASSERT(iszeroish(ia));
+		ASSERT(iszeroish(ib));
+	} else {
+		if (!iszeroish(ia) || !iszeroish(ib))
+			return;
+	}
+
+	ngot = solve_quadratic(got, k, -k*ra, k*rb);
+
+	for (i = 0; i < sizeof(found) / sizeof(*found); i++)
+		found[i] = (i >= nexpected || expected[i] < 0.0 + TOLERANCE || expected[i] > 1.0 - TOLERANCE);
+	if (!nexpected) {
+		ASSERT(!ngot);
+		return;
+	}
+	for (i = 0; i < ngot; i++) {
+		correct = (got[i] < 0.0 + TOLERANCE || got[i] > 1.0 - TOLERANCE);
+		error = fabs(got[i] - expected[0]);
+		for (j = 0; j < nexpected; j++) {
+			error = fmin(error, fabs(got[i] - expected[j]));
+			if (tolerant_eq(got[i], expected[j]))
+				found[j] = correct = 1;
+		}
+		if (!correct)
+			maxerror = fmax(maxerror, error);
+	}
+	if (!found[0] || !found[1] || !tolerant_eq(maxerror, 0)) {
+		fprintf(stderr, "\n%lg, %lg, %lg, %lg, %lg\n", r1, i1, r2, i2, k);
+		fprintf(stderr, "\t(%la, %la, %la, %la, %la)\n", r1, i1, r2, i2, k);
+		fprintf(stderr, "\tmaxerror: %lf (%la)\n", maxerror, maxerror);
+		fprintf(stderr, "\tfound: %i, %i\n", found[0], found[1]);
+		fprintf(stderr, "\tgot: %zu\n", ngot);
+		for (i = 0; i < ngot; i++)
+			fprintf(stderr, "\t\t%lf (%la)\n", got[i], got[i]);
+		fprintf(stderr, "\texpected: %zu\n", nexpected);
+		for (i = 0; i < nexpected; i++)
+			fprintf(stderr, "\t\t%lf (%la)\n", expected[i], expected[i]);
+	}
+	ASSERT(found[0]);
+	ASSERT(found[1]);
+	ASSERT(tolerant_eq(maxerror, 0));
+}
+
+
+static void
+check_linear(double r1, double i1, double k, int strict_real)
+{
+	size_t nexpected = 0, ngot, i;
+	double expected[1] = {999999}, got[1];
+	double ra;
+	double ia;
+
+	if (iszeroish(i1) && 0.0 <= r1 && r1 <= 1.0)
+		expected[nexpected++] = r1;
+
+	/*
+	 * k(t - t₁) = 0
+	 * kt - kt₁ = 0
+	 */
+
+	ra = r1;
+	ia = i1;
+
+	if (strict_real) {
+		ASSERT(iszeroish(ia));
+	} else {
+		if (!iszeroish(ia))
+			return;
+	}
+
+	ngot = solve_linear(got, k, -k*ra);
+
+	ASSERT(ngot == nexpected);
+	for (i = 0; i < ngot; i++)
+		ASSERT(eq(got[i], expected[i]));
+}
+
+
+static void
+check_cubic_randomly(void)
+{
+	double r1 = rand() / (double)RAND_MAX;
+	double i1 = rand() / (double)RAND_MAX;
+	double r2 = rand() / (double)RAND_MAX;
+	double i2 = rand() / (double)RAND_MAX;
+	double r3 = rand() / (double)RAND_MAX;
+	double i3 = rand() / (double)RAND_MAX;
+	double k = rand() / (double)RAND_MAX;
+	r1 = r1 * 3 - 1;
+	i1 = i1 * 3 - 1;
+	r2 = r2 * 3 - 1;
+	i2 = i2 * 3 - 1;
+	r3 = r3 * 3 - 1;
+	i3 = i3 * 3 - 1;
+	k = k * 3 - 1;
+	k += iszeroish(k);
+
+	check_cubic(Cx(r1, i1), Cx(r2, i2), Cx(r3, i3), k, 0);
+	check_cubic(Im(i1), Cx(r2, i2), Cx(r3, i3), k, 0);
+	check_cubic(Re(r1), Cx(r2, i2), Cx(r3, i3), k, 0);
+	check_cubic(Cx(r1, i1), Im(i2), Cx(r3, i3), k, 0);
+	check_cubic(Im(i1), Im(i2), Cx(r3, i3), k, 0);
+	check_cubic(Re(r1), Im(i2), Cx(r3, i3), k, 0);
+	check_cubic(Cx(r1, i1), Re(r2), Cx(r3, i3), k, 0);
+	check_cubic(Im(i1), Re(r2), Cx(r3, i3), k, 0);
+	check_cubic(Re(r1), Re(r2), Cx(r3, i3), k, 0);
+	check_cubic(Cx(r1, i1), Cx(r1, i1), Cx(r3, i3), k, 0);
+	check_cubic(Im(i1), Im(i1), Cx(r3, i3), k, 0);
+	check_cubic(Re(r1), Re(r1), Cx(r3, i3), k, 0);
+	check_cubic(Cx(r1, i1), Cx(r3, i3), Cx(r3, i3), k, 0);
+	check_cubic(Im(i1), Cx(r3, i3), Cx(r3, i3), k, 0);
+	check_cubic(Re(r1), Cx(r3, i3), Cx(r3, i3), k, 0);
+
+	check_cubic(Cx(r1, i1), Cx(r2, i2), Im(i3), k, 0);
+	check_cubic(Im(i1), Cx(r2, i2), Im(i3), k, 0);
+	check_cubic(Re(r1), Cx(r2, i2), Im(i3), k, 0);
+	check_cubic(Cx(r1, i1), Im(i2), Im(i3), k, 0);
+	check_cubic(Im(i1), Im(i2), Im(i3), k, 0);
+	check_cubic(Re(r1), Im(i2), Im(i3), k, 0);
+	check_cubic(Cx(r1, i1), Re(r2), Im(i3), k, 0);
+	check_cubic(Im(i1), Re(r2), Im(i3), k, 0);
+	check_cubic(Re(r1), Re(r2), Im(i3), k, 0);
+	check_cubic(Cx(r1, i1), Cx(r1, i1), Im(i3), k, 0);
+	check_cubic(Im(i1), Im(i1), Im(i3), k, 0);
+	check_cubic(Re(r1), Re(r1), Im(i3), k, 0);
+	check_cubic(Cx(r1, i1), Im(i3), Im(i3), k, 0);
+	check_cubic(Im(i1), Im(i3), Im(i3), k, 0);
+	check_cubic(Re(r1), Im(i3), Im(i3), k, 0);
+
+	check_cubic(Cx(r1, i1), Cx(r2, i2), Re(r3), k, 0);
+	check_cubic(Im(i1), Cx(r2, i2), Re(r3), k, 0);
+	check_cubic(Re(r1), Cx(r2, i2), Re(r3), k, 0);
+	check_cubic(Cx(r1, i1), Im(i2), Re(r3), k, 0);
+	check_cubic(Im(i1), Im(i2), Re(r3), k, 0);
+	check_cubic(Re(r1), Im(i2), Re(r3), k, 0);
+	check_cubic(Cx(r1, i1), Re(r2), Re(r3), k, 0);
+	check_cubic(Im(i1), Re(r2), Re(r3), k, 0);
+	check_cubic(Re(r1), Re(r2), Re(r3), k, 1);
+	check_cubic(Cx(r1, i1), Cx(r1, i1), Re(r3), k, 0);
+	check_cubic(Im(i1), Im(i1), Re(r3), k, 0);
+	check_cubic(Re(r1), Re(r1), Re(r3), k, 1);
+	check_cubic(Cx(r1, i1), Re(r3), Re(r3), k, 0);
+	check_cubic(Im(i1), Re(r3), Re(r3), k, 0);
+	check_cubic(Re(r1), Re(r3), Re(r3), k, 1);
+
+	check_cubic(Cx(r1, i1), Cx(r1, i1), Cx(r1, i1), k, 0);
+	check_cubic(Im(i1), Im(i1), Im(i1), k, 0);
+	check_cubic(Re(r1), Re(r1), Re(r1), k, 1);
+}
+
+
+static void
+check_quadratic_randomly(void)
+{
+	double r1 = rand() / (double)RAND_MAX;
+	double i1 = rand() / (double)RAND_MAX;
+	double r2 = rand() / (double)RAND_MAX;
+	double i2 = rand() / (double)RAND_MAX;
+	double k = rand() / (double)RAND_MAX;
+	r1 = r1 * 3 - 1;
+	i1 = i1 * 3 - 1;
+	r2 = r2 * 3 - 1;
+	i2 = i2 * 3 - 1;
+	k = k * 3 - 1;
+	k += iszeroish(k);
+
+	check_quadratic(Cx(r1, i1), Cx(r2, i2), k, 0);
+	check_quadratic(Im(i1), Cx(r2, i2), k, 0);
+	check_quadratic(Re(r1), Cx(r2, i2), k, 0);
+	check_quadratic(Cx(r1, i1), Im(i2), k, 0);
+	check_quadratic(Im(i1), Im(i2), k, 0);
+	check_quadratic(Re(r1), Im(i2), k, 0);
+	check_quadratic(Cx(r1, i1), Re(r2), k, 0);
+	check_quadratic(Im(i1), Re(r2), k, 0);
+	check_quadratic(Re(r1), Re(r2), k, 1);
+	check_quadratic(Cx(r1, i1), Cx(r1, i1), k, 0);
+	check_quadratic(Im(i1), Im(i1), k, 0);
+	check_quadratic(Re(r1), Re(r1), k, 1);
+}
+
+
+static void
+check_linear_randomly(void)
+{
+	double r1 = rand() / (double)RAND_MAX;
+	double i1 = rand() / (double)RAND_MAX;
+	double k = rand() / (double)RAND_MAX;
+	r1 = r1 * 3 - 1;
+	i1 = i1 * 3 - 1;
+	k = k * 3 - 1;
+	k += iszeroish(k);
+
+	check_linear(Cx(r1, i1), k, 0);
+	check_linear(Im(i1), k, 0);
+	check_linear(Re(r1), k, 0);
+}
+
+
+int
+main(void)
+{
+	size_t i;
+
+	void *r = malloc(1);
+	srand((unsigned)(uintptr_t)r);
+	free(r);
+
+	check_linear(Re(0), 1.0, 1);
+	check_linear(Re(0), 1.5, 1);
+	check_linear(Re(0), 0.5, 1);
+	check_linear(Re(0), -0.5, 1);
+	check_linear(Re(0), -1.0, 1);
+	check_linear(Re(0), -1.5, 1);
+	for (i = 0; i < 100000UL; i++)
+		check_linear_randomly();
+
+	check_quadratic(Re(0), Re(0), 1.0, 1);
+	check_quadratic(Re(0), Re(0), 1.5, 1);
+	check_quadratic(Re(0), Re(0), 0.5, 1);
+	check_quadratic(Re(0), Re(0), -0.5, 1);
+	check_quadratic(Re(0), Re(0), -1.0, 1);
+	check_quadratic(Re(0), Re(0), -1.5, 1);
+	for (i = 0; i < 10000000UL; i++)
+		check_quadratic_randomly();
+
+	check_cubic(Re(0), Re(0), Re(0), 1.0, 1);
+	check_cubic(Re(0), Re(0), Re(0), 1.5, 1);
+	check_cubic(Re(0), Re(0), Re(0), 0.5, 1);
+	check_cubic(Re(0), Re(0), Re(0), -0.5, 1);
+	check_cubic(Re(0), Re(0), Re(0), -1.0, 1);
+	check_cubic(Re(0), Re(0), Re(0), -1.5, 1);
+	for (i = 0; i < 1000000UL; i++)
+		check_cubic_randomly();
+
+	return 0;
+}
+
+
+#endif
diff --git a/librifunktionsteckensnittsglyfrasteriseringsprogrambiblioteket.h b/librifunktionsteckensnittsglyfrasteriseringsprogrambiblioteket.h
new file mode 100644
index 0000000..0fceece
--- /dev/null
+++ b/librifunktionsteckensnittsglyfrasteriseringsprogrambiblioteket.h
@@ -0,0 +1,33 @@
+/* See LICENSE file for copyright and license details. */
+#ifndef LIBRIFUNKTIONSTECKENSNITTSGLYFRASTERISERINGSPROGRAMBIBLIOTEKET_H
+#define LIBRIFUNKTIONSTECKENSNITTSGLYFRASTERISERINGSPROGRAMBIBLIOTEKET_H
+
+#include <stddef.h>
+
+
+typedef struct rtgrpblib_raster RTGRPBLIB_RASTER;
+
+
+RTGRPBLIB_RASTER *rtgrpblib_create_raster(size_t width, size_t height);
+int rtgrpblib_reset_raster(RTGRPBLIB_RASTER *raster, size_t width, size_t height);
+
+void rtgrpblib_set_draftness(RTGRPBLIB_RASTER *raster, double draftness);
+
+void rtgrpblib_fill_shapes(double *restrict image, size_t rowsize, const RTGRPBLIB_RASTER *raster);
+
+void rtgrpblib_draw_linear_bezier(RTGRPBLIB_RASTER *restrict raster,
+                                  double x1, double y1,
+                                  double x2, double y2);
+
+void rtgrpblib_draw_quadratic_bezier(RTGRPBLIB_RASTER *restrict raster,
+                                     double x1, double y1,
+                                     double x2, double y2,
+                                     double x3, double y3);
+
+void rtgrpblib_draw_cubic_bezier(RTGRPBLIB_RASTER *restrict raster,
+                                 double x1, double y1,
+                                 double x2, double y2,
+                                 double x3, double y3,
+                                 double x4, double y4);
+
+#endif
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
diff --git a/mk/linux.mk b/mk/linux.mk
new file mode 100644
index 0000000..ad58f69
--- /dev/null
+++ b/mk/linux.mk
@@ -0,0 +1,6 @@
+LIBEXT      = so
+LIBFLAGS    = -shared -Wl,-soname,lib$(LIB_NAME).$(LIBEXT).$(LIB_MAJOR)
+LIBMAJOREXT = $(LIBEXT).$(LIB_MAJOR)
+LIBMINOREXT = $(LIBEXT).$(LIB_VERSION)
+
+FIX_INSTALL_NAME = :
diff --git a/mk/macos.mk b/mk/macos.mk
new file mode 100644
index 0000000..265b653
--- /dev/null
+++ b/mk/macos.mk
@@ -0,0 +1,6 @@
+LIBEXT      = dylib
+LIBFLAGS    = -dynamiclib -Wl,-compatibility_version,$(LIB_MAJOR) -Wl,-current_version,$(LIB_VERSION)
+LIBMAJOREXT = $(LIB_MAJOR).$(LIBEXT)
+LIBMINOREXT = $(LIB_VERSION).$(LIBEXT)
+
+FIX_INSTALL_NAME = install_name_tool -id "$(PREFIX)/lib/librifunktionsteckensnittsglyfrasteriseringsprogrambiblioteket.$(LIBMAJOREXT)"
diff --git a/mk/windows.mk b/mk/windows.mk
new file mode 100644
index 0000000..ed5ec8d
--- /dev/null
+++ b/mk/windows.mk
@@ -0,0 +1,6 @@
+LIBEXT      = dll
+LIBFLAGS    = -shared
+LIBMAJOREXT = $(LIB_MAJOR).$(LIBEXT)
+LIBMINOREXT = $(LIB_VERSION).$(LIBEXT)
+
+FIX_INSTALL_NAME = :
diff --git a/rtgrpblib_create_raster.c b/rtgrpblib_create_raster.c
new file mode 100644
index 0000000..89ebad6
--- /dev/null
+++ b/rtgrpblib_create_raster.c
@@ -0,0 +1,40 @@
+/* See LICENSE file for copyright and license details. */
+#include "common.h"
+#ifndef TEST
+
+
+RASTER *
+rtgrpblib_create_raster(size_t width, size_t height)
+{
+	RASTER *raster;
+
+	if (!width || !height) {
+		errno = EINVAL;
+		return NULL;
+	}
+
+	if (width > (SIZE_MAX - offsetof(RASTER, cells)) / sizeof(*raster->cells) / height)
+		goto enomem;
+	raster = calloc(1, offsetof(RASTER, cells) + height * width * sizeof(*raster->cells));
+	if (!raster) {
+	enomem:
+		errno = ENOMEM;
+		return NULL;
+	}
+
+	raster->width = width;
+	raster->height = height;
+	raster->draftness = 0.5;
+	return raster;
+}
+
+
+#else
+
+int
+main(void)
+{
+	return 0; /* TODO add test */
+}
+
+#endif
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
diff --git a/rtgrpblib_draw_linear_bezier.c b/rtgrpblib_draw_linear_bezier.c
new file mode 100644
index 0000000..5c83591
--- /dev/null
+++ b/rtgrpblib_draw_linear_bezier.c
@@ -0,0 +1,116 @@
+/* See LICENSE file for copyright and license details. */
+#include "common.h"
+#ifndef TEST
+
+
+void
+rtgrpblib_draw_linear_bezier(RASTER *restrict raster, double x1, double y1, double x2, double y2)
+{
+	double w, h, dx, dy, x, y;
+	int xdir, ydir;
+
+	dx = x2 - x1;
+	dy = y2 - y1;
+	xdir = SIGNUM(dx);
+	ydir = SIGNUM(dy);
+
+	if (!ydir) {
+		/* For horizontal lines, it is enough that we have the
+		 * corners mapped onto the raster, which we must have
+		 * since a glyph cannot just contain a line, but outlines:
+		 * in `fill`, the line will not make a contribution because
+		 * `.opposite_coverage` for each cell will be zero, and
+		 * consequentially, `.cell_coverage` will also be zero */
+		return;
+	}
+
+	/* We cut of everything above and below the raster,
+	 * as these do not contribute to the result */
+	h = (double)raster->height;
+	if (y1 < 0 && y2 < 0)
+		return;
+	if (y1 >= h && y2 >= h)
+		return;
+	if (y1 < 0 && y2 >= 0) {
+		x1 = x1 + (0 - y1) * dx / dy;
+		y1 = 0;
+	} else if (y2 < 0 && y1 >= 0) {
+		x2 = x1 + (0 - y1) * dx / dy;
+		y2 = 0;
+	}
+	if (y1 < h && y2 >= h) {
+		x2 = x1 + (h - y1) * dx / dy;
+		y2 = h;
+	} else if (y2 < h && y1 >= h) {
+		x1 = x1 + (h - y1) * dx / dy;
+		y1 = h;
+	}
+
+	/* Dealing with the left and the right section
+	 * is not as simple as `fill` assumes that the
+	 * sum of each line's `.opposite_coverage` is 0.
+	 * Therefore must must identify the parts of the
+	 * line that are left of the raster, right of
+	 * the raster, and inside the raster. The parts
+	 * that are outside the raster are move to the
+	 * edge, or rather their show on the edges are
+	 * used. */
+	w = (double)raster->width;
+	if (x1 <= 0 && x2 <= 0) {
+		draw_vertical_line(raster, 0, y1, y2, ydir);
+		return;
+	}
+	if (x1 >= w && x2 >= w) {
+		draw_vertical_line_opposite_only(raster, raster->width - 1, y1, y2, ydir);
+		return;
+	}
+	if (x1 < 0 && x2 >= 0) {
+		y = y1 + (0 - x1) * dy / dx;
+		x = 0;
+		draw_vertical_line(raster, 0, y1, y, ydir);
+		x1 = x;
+		y1 = y;
+	} else if (x2 < 0 && x1 >= 0) {
+		y = y1 + (0 - x1) * dy / dx;
+		x = 0;
+		draw_vertical_line(raster, 0, y, y2, ydir);
+		x2 = x;
+		y2 = y;
+	}
+	if (x1 < w && x2 >= w) {
+		y = y1 + (w - x1) * dy / dx;
+		x = w;
+		draw_vertical_line_opposite_only(raster, raster->width - 1, y, y2, ydir);
+		x2 = x;
+		y2 = y;
+	} else if (x2 < w && x1 >= w) {
+		y = y1 + (w - x1) * dy / dx;
+		x = w;
+		draw_vertical_line_opposite_only(raster, raster->width - 1, y1, y, ydir);
+		x1 = x;
+		y1 = y;
+	}
+
+	/* Now we can finally draw the part of the
+	 * line that is inside the raster */
+	if (!xdir) {
+		/* Optimisation for vertical lines. It also serves
+		 * to illustrate how the algorithm is designed. */
+		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
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
diff --git a/rtgrpblib_fill_shapes.c b/rtgrpblib_fill_shapes.c
new file mode 100644
index 0000000..285a1fd
--- /dev/null
+++ b/rtgrpblib_fill_shapes.c
@@ -0,0 +1,38 @@
+/* See LICENSE file for copyright and license details. */
+#include "common.h"
+#ifndef TEST
+
+
+void
+rtgrpblib_fill_shapes(double *restrict image, size_t rowsize, const RASTER *raster)
+{
+	const CELL *cells = raster->cells;
+	double inklevel;
+	size_t x, y;
+
+#ifndef ROWWISE_RESET_INKLEVEL
+	inklevel = 0;
+#endif
+	for (y = 0; y < raster->height; y++) {
+#ifdef ROWWISE_RESET_INKLEVEL
+		inklevel = 0;
+#endif
+		for (x = 0; x < raster->width; x++) {
+			image[x] = fmin(fabs(inklevel + cells[x].cell_coverage), 1.0);
+			inklevel += cells[x].opposite_coverage;
+		}
+		image = &image[rowsize];
+		cells = &cells[raster->width];
+	}
+}
+
+
+#else
+
+int
+main(void)
+{
+	return 0; /* TODO add test */
+}
+
+#endif
diff --git a/rtgrpblib_reset_raster.c b/rtgrpblib_reset_raster.c
new file mode 100644
index 0000000..2cfef42
--- /dev/null
+++ b/rtgrpblib_reset_raster.c
@@ -0,0 +1,29 @@
+/* See LICENSE file for copyright and license details. */
+#include "common.h"
+#ifndef TEST
+
+
+int
+rtgrpblib_reset_raster(RASTER *raster, size_t width, size_t height)
+{
+	if (!width || !height || width > raster->width * raster->height / height) {
+		errno = EINVAL;
+		return -1;
+	}
+
+	raster->width = width;
+	raster->height = height;
+	memset(raster->cells, 0, width * height * sizeof(*raster->cells));
+	return 0;
+}
+
+
+#else
+
+int
+main(void)
+{
+	return 0; /* TODO add test */
+}
+
+#endif
diff --git a/rtgrpblib_set_draftness.c b/rtgrpblib_set_draftness.c
new file mode 100644
index 0000000..dbdeed0
--- /dev/null
+++ b/rtgrpblib_set_draftness.c
@@ -0,0 +1,21 @@
+/* See LICENSE file for copyright and license details. */
+#include "common.h"
+#ifndef TEST
+
+
+void
+rtgrpblib_set_draftness(RTGRPBLIB_RASTER *raster, double draftness)
+{
+	raster->draftness = draftness;
+}
+
+
+#else
+
+int
+main(void)
+{
+	return 0; /* TODO add test */
+}
+
+#endif
diff --git a/sorting.c b/sorting.c
new file mode 100644
index 0000000..e4bfdaf
--- /dev/null
+++ b/sorting.c
@@ -0,0 +1,25 @@
+/* See LICENSE file for copyright and license details. */
+#include "common.h"
+#ifndef TEST
+
+
+int
+doublepcmp(const void *avp, const void *bvp)
+{
+	const double *ap = avp;
+	const double *bp = bvp;
+	return (*ap > *bp) - (*ap < *bp);
+}
+
+
+#else
+
+
+int
+main(void)
+{
+	return 0; /* TODO add test */
+}
+
+
+#endif