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