diff options
Diffstat (limited to '')
| -rw-r--r-- | rtgrpblib_draw_cubic_bezier.c | 205 |
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 |