/* 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;
/* 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 ((y2 - y1) * (x3 - x2) == (y3 - y2) * (x2 - x1)) {
/* Drawing (x1, y1)--(x2, y2) plus (x2, y2)--(x3, y3)
* is equivalent to drawing (x1, y1)--(x3, y3) because
* the result is signed and if (x2, y2) is outside the
* (x1, y1)--(x3, y3) line (x2, y2)--(x3, y3) erases
* the part that (x1, y1)--(x2, y2) extends to
* (x1, y1)--(x3, y3) */
rtgrpblib_draw_linear_bezier(raster, x1, y1, x3, y3);
return;
}
/* 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
static RASTER *raster;
static RASTER *refraster;
static double
frand(double a, double b)
{
int i = rand();
double f = (double)i / (double)RAND_MAX;
double max = fmax(a, b);
double min = fmin(a, b);
return fma(f, max - min, min);
}
static void
check_line(double x1, double y1, double x2, double y2, double x3, double y3, int exact)
{
size_t i;
alarm(5);
rtgrpblib_draw_quadratic_bezier(raster, x1, y1, x2, y2, x3, y3);
alarm(5);
rtgrpblib_draw_linear_bezier(refraster, x1, y1, x3, y3);
alarm(0);
if (exact) {
ASSERT(!memcmp(raster->cells, refraster->cells, 100 * sizeof(*raster->cells)));
} else {
for (i = 0; i < 100; i++) {
if (!eq(raster->cells[i].cell_coverage, refraster->cells[i].cell_coverage))
fprintf(stderr, "c%02zu: %la ref: %la; diff: %lg=%la (tol: %la)\n", i, raster->cells[i].cell_coverage, refraster->cells[i].cell_coverage, refraster->cells[i].cell_coverage - raster->cells[i].cell_coverage, refraster->cells[i].cell_coverage - raster->cells[i].cell_coverage, TOLERANCE);
if (!eq(raster->cells[i].opposite_coverage, refraster->cells[i].opposite_coverage))
fprintf(stderr, "o%02zu: %la ref: %la; diff: %lg=%la (tol: %la)\n", i, raster->cells[i].opposite_coverage, refraster->cells[i].opposite_coverage, raster->cells[i].opposite_coverage - refraster->cells[i].opposite_coverage, raster->cells[i].opposite_coverage - refraster->cells[i].opposite_coverage, TOLERANCE);
}
for (i = 0; i < 100; i++) {
if (!eq(raster->cells[i].cell_coverage, refraster->cells[i].cell_coverage) ||
!eq(raster->cells[i].opposite_coverage, refraster->cells[i].opposite_coverage)) {
print_raster(refraster);
print_raster(raster);
}
ASSERT(eq(raster->cells[i].cell_coverage, refraster->cells[i].cell_coverage));
ASSERT(eq(raster->cells[i].opposite_coverage, refraster->cells[i].opposite_coverage));
}
}
rtgrpblib_reset_raster(raster, 10, 10);
rtgrpblib_reset_raster(refraster, 10, 10);
}
static void
check_simple_lines(void)
{
check_line(1, 1, 2, 2, 3, 3, 1);
check_line(1, 1, 3, 3, 2, 2, 1);
check_line(2, 2, 1, 1, 3, 3, 1);
check_line(2, 2, 3, 3, 1, 1, 1);
check_line(3, 3, 1, 1, 2, 2, 1);
check_line(3, 3, 2, 2, 1, 1, 1);
check_line(1, 1, 3, 3, 3, 3, 1);
check_line(1, 1, 1, 1, 3, 3, 1);
check_line(1, 4, 2, 4, 3, 4, 1);
check_line(1, 4, 3, 4, 2, 4, 1);
check_line(2, 4, 1, 4, 3, 4, 1);
check_line(2, 4, 3, 4, 1, 4, 1);
check_line(3, 4, 1, 4, 2, 4, 1);
check_line(3, 4, 2, 4, 1, 4, 1);
check_line(1, 4, 3, 4, 3, 4, 1);
check_line(1, 4, 1, 4, 3, 4, 1);
check_line(4, 1, 4, 2, 4, 3, 1);
check_line(4, 1, 4, 3, 4, 2, 1);
check_line(4, 2, 4, 1, 4, 3, 1);
check_line(4, 2, 4, 3, 4, 1, 1);
check_line(4, 3, 4, 1, 4, 2, 1);
check_line(4, 3, 4, 2, 4, 1, 1);
check_line(4, 1, 4, 3, 4, 3, 1);
check_line(4, 1, 4, 1, 4, 3, 1);
}
static void
check_random_bounded_lines(int exact)
{
double x1, x2, x3;
double y1, y2, y3;
double mid;
size_t i;
for (i = 0; i < 100000UL; i++) {
x1 = floor(frand(0, 11));
y1 = floor(frand(0, 11));
x3 = floor(frand(0, 11));
y3 = floor(frand(0, 11));
x1 -= (double)(x1 > 10);
x2 -= (double)(x2 > 10);
y1 -= (double)(y1 > 10);
y2 -= (double)(y2 > 10);
mid = frand(0, 1);
if (x3 - x1) {
x2 = floor(mid * (x3 - x1) + x1);
y2 = y1 + (x2 - x1) * (y3 - y1) / (x3 - x1);
} else if (y3 - y1) {
x2 = x1;
y2 = floor(mid * (y3 - y1) + y1);
} else {
x2 = x1;
y2 = y1;
}
if (!exact)
check_line(x1, y1, x2, y2, x3, y3, x2 == floor(x2) && y2 == floor(y2));
else if (x2 == floor(x2) && y2 == floor(y2))
check_line(x1, y1, x2, y2, x3, y3, 1);
}
}
static void
check_random_unbounded_lines(int exact)
{
double x1, x2, x3;
double y1, y2, y3;
double mid;
size_t i;
for (i = 0; i < 100000UL; i++) {
x1 = floor(frand(-5, 15));
y1 = floor(frand(-5, 15));
x3 = floor(frand(-5, 15));
y3 = floor(frand(-5, 15));
mid = frand(-2, 2);
if (x3 - x1) {
x2 = floor(mid * (x3 - x1) + x1);
y2 = y1 + (x2 - x1) * (y3 - y1) / (x3 - x1);
} else if (y3 - y1) {
x2 = x1;
y2 = floor(mid * (y3 - y1) + y1);
} else {
x2 = x1;
y2 = y1;
}
if (!exact)
check_line(x1, y1, x2, y2, x3, y3, x2 == floor(x2) && y2 == floor(y2));
else if (x2 == floor(x2) && y2 == floor(y2))
check_line(x1, y1, x2, y2, x3, y3, 1);
}
}
int
main(void)
{
raster = rtgrpblib_create_raster(10, 10);
ASSERT(raster);
refraster = rtgrpblib_create_raster(10, 10);
ASSERT(refraster);
srand((unsigned)(uintptr_t)raster);
check_simple_lines();
check_random_bounded_lines(1);
check_random_unbounded_lines(1);
//TODO check_random_bounded_lines(0);
//TODO check_random_unbounded_lines(0);
raster->draftness = -1;
refraster->draftness = -1;
check_simple_lines();
check_random_bounded_lines(1);
check_random_unbounded_lines(1);
#ifdef TODO_NOT_IMPLEMENTED
check_random_bounded_lines(0);
check_random_unbounded_lines(0);
#endif
raster->draftness = 0.5;
refraster->draftness = 0.5;
/* TODO add tests */
free(raster);
free(refraster);
return 0;
}
#endif
