/* 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);
/* TODO deal with unchanged t and oversized steps (e.g. locally zero derivative) */
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
