/* 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;
start -= (size_t)(start == raster->height);
end -= (size_t)(end == raster->height);
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
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_against_reference(double x1, double y1, double x2, double y2)
{
size_t i;
alarm(5);
draw_linear_bezier_reference(refraster, x1, y1, x2, y2);
alarm(0);
for (i = 0; i < 100; i++) {
ASSERT(eq(raster->cells[i].cell_coverage, refraster->cells[i].cell_coverage));
ASSERT(eq(raster->cells[i].opposite_coverage, refraster->cells[i].opposite_coverage));
}
}
static void
check_vertical_line(double x, double y1, double y2)
{
alarm(5);
draw_vertical_line(raster, x, y1, y2, SIGNUM(y2 - y1));
alarm(0);
check_against_reference(x, y1, x, y2);
}
static void
check_diagonal_line(double x1, double y1, double x2, double y2)
{
double dx = x2 - x1;
double dy = y2 - y1;
alarm(5);
draw_diagonal_line(raster, x1, y1, x2, y2, dx, dy, SIGNUM(dx), SIGNUM(dy));
alarm(0);
check_against_reference(x1, y1, x2, y2);
}
static void
check_bounded_line(double x1, double y1, double x2, double y2)
{
alarm(5);
draw_bounded_line(raster, x1, y1, x2, y2);
alarm(0);
check_against_reference(x1, y1, x2, y2);
}
#ifndef ROWWISE_RESET_INKLEVEL
static void
check_vertical_line_opposite_only(double y1, double y2)
{
alarm(5);
draw_vertical_line_opposite_only(raster, 9, y1, y2, SIGNUM(y2 - y1));
alarm(0);
check_against_reference(20, y1, 20, y2);
}
#endif
static void
reset(void)
{
rtgrpblib_reset_raster(raster, 10, 10);
rtgrpblib_reset_raster(refraster, 10, 10);
}
int
main(void)
{
double x1, x2;
double y1, y2;
size_t i;
raster = rtgrpblib_create_raster(10, 10);
ASSERT(raster);
refraster = rtgrpblib_create_raster(10, 10);
ASSERT(refraster);
srand((unsigned)(uintptr_t)raster);
for (i = 0; i < 10000UL; i++) {
x1 = frand(0, 10);
y1 = frand(0, 10);
x2 = frand(0, 10);
y2 = frand(0, 10);
check_vertical_line(x1, y1, y2);
reset();
if (x1 != x2) {
check_diagonal_line(x1, y1, x2, y2);
reset();
}
check_bounded_line(x1, y1, x1, y2);
reset();
check_bounded_line(x1, y1, x2, y2);
reset();
#ifndef ROWWISE_RESET_INKLEVEL
check_vertical_line_opposite_only(y1, y2);
reset();
#endif
}
free(raster);
free(refraster);
return 0;
}
#endif
