/* See LICENSE file for copyright and license details. */
#include "common.h"
#ifndef TEST
/*
* The implementation in this fill is not supposed to be fast,
* it is supposed to be absolutely correct. The implementation
* shall be as clear and obvious as possible. It is important
* that everything written to cells is as precise as possible
* and that the function behavious correctly when outside the
* raster. This implementation shall be usable to generate an
* as perfect as possible rasterisation and shall be usable to
* test the correctness and accuracy of the fast implementation.
*/
void
draw_linear_bezier_reference(RASTER *restrict raster, double x1, double y1, double x2, double y2)
{
double dx = x2 - x1;
double dy = y2 - y1;
int xdir = SIGNUM(dx);
int ydir = SIGNUM(dy);
double kx = dx / dy;
double ky = dy / dx;
double prevX, prevY, x, y, cx, cy, h2x, h2y;
double xfloor, yfloor;
ssize_t cellX, cellY;
size_t cell;
prevX = x1;
prevY = y1;
for (; prevX != x2 || prevY != y2; prevX = x, prevY = y) {
/* Find next vertical and next hozitonal grid line */
x = xdir < 0 ? floor(prevX) : xdir > 0 ? ceil(prevX) : prevX;
y = ydir < 0 ? floor(prevY) : ydir > 0 ? ceil(prevY) : prevY;
x = x == prevX ? x + xdir : x;
y = y == prevY ? y + ydir : y;
/* Limit to the extend of the drawn line segment */
if ((xdir < 0 && x < x2) || (xdir > 0 && x > x2))
x = x2;
if ((ydir < 0 && y < y2) || (ydir > 0 && y > y2))
y = y2;
/* Find next cell edge */
if (xdir && ydir) {
cy = y1 + (x - x1) * ky;
cx = x1 + (y - y1) * kx;
h2x = (cy - prevY) * (cy - prevY) + (x - prevX) * (x - prevX);
h2y = (cx - prevX) * (cx - prevX) + (y - prevY) * (y - prevY);
if (h2x < h2y)
y = cy;
else
x = cx;
/* Make sure loop will exit when we reach the end */
if ((xdir < 0 && x <= x2) || (xdir > 0 && x >= x2) ||
(ydir < 0 && y <= y2) || (ydir > 0 && y >= y2)) {
x = x2;
y = y2;
}
}
/* Select cell to draw in */
xfloor = xdir >= 0 ? floor(prevX) : floor(x);
yfloor = ydir >= 0 ? floor(prevY) : floor(y);
cellX = (ssize_t)xfloor;
cellY = (ssize_t)yfloor;
/* Do not draw if above or below the raster */
if (cellY < 0 || (size_t)cellY >= raster->height)
continue;
if (cellX < 0) {
/* Draw on first column in raster of outside of raster on the left side,
* with full horizontal coverage, `.opposite_coverage` one the actual
* cell (outside of the raster) will contribute to the first column's
* `.cell_coverage` */
cell = (size_t)cellY * raster->width;
raster->cells[cell].opposite_coverage += y - prevY;
raster->cells[cell].cell_coverage += y - prevY;
} else if ((size_t)cellX >= raster->width) {
/* If outside of the raster, on the right side, add to the last
* columns contributions to the next column, so that the ink level
* balances out to zero */
cell = (size_t)cellY * raster->width + (raster->width - 1);
raster->cells[cell].opposite_coverage += y - prevY;
} else {
/* If inside the raster, just draw the line */
cell = (size_t)cellY * raster->width + (size_t)cellX;
raster->cells[cell].opposite_coverage += y - prevY;
raster->cells[cell].cell_coverage += (y - prevY) * fabs(xfloor + 1.0 - (x + prevX) / 2.0);
}
}
}
#else
#define CHECK_ZEROED()\
do {\
size_t i__;\
for (i__ = 0; i__ < 100; i__++) {\
ASSERT(!raster->cells[i__].cell_coverage);\
ASSERT(!raster->cells[i__].opposite_coverage);\
}\
} while (0)
#define CHECK_CELL(X, Y, AREA, SHADOW)\
do {\
size_t i__ = (size_t)(Y) * raster->width + (size_t)(X);\
ASSERT(raster->cells[i__].cell_coverage == (AREA));\
ASSERT(raster->cells[i__].opposite_coverage == (SHADOW));\
raster->cells[i__].cell_coverage = 0;\
raster->cells[i__].opposite_coverage = 0;\
} while (0)
#define SUB_CELL(X, Y, AREA, SHADOW)\
do {\
size_t i__ = (size_t)(Y) * raster->width + (size_t)(X);\
raster->cells[i__].cell_coverage -= (AREA);\
raster->cells[i__].opposite_coverage -= (SHADOW);\
} while (0)
static RASTER *raster;
static void
draw(double x1, double y1, double x2, double y2)
{
alarm(5);
draw_linear_bezier_reference(raster, x1, y1, x2, y2);
alarm(0);
}
static void
check_horizontal(void)
{
ssize_t y, x1, x2;
double fy, fx1, fx2;
for (y = -2; y < 12; y++) {
for (x1 = -2; x1 < 12; x1++) {
for (x2 = -2; x2 < 12; x2++) {
fy = (double)y;
fx1 = (double)x1;
fx2 = (double)x2;
draw(fx1, fy, fx2, fy);
CHECK_ZEROED();
draw(fx1 + 0.25, fy, fx2 + 0.25, fy);
CHECK_ZEROED();
draw(fx1, fy + 0.25, fx2, fy + 0.25);
CHECK_ZEROED();
draw(fx1 + 0.25, fy + 0.25, fx2 + 0.25, fy + 0.25);
CHECK_ZEROED();
}
}
}
}
static void
check_vertical(void)
{
ssize_t x, y1, y2, i, ydir;
double yhit;
double fx, fy1, fy2;
for (x = -2; x < 12; x++) {
for (y1 = -2; y1 < 12; y1++) {
for (y2 = -2; y2 < 12; y2++) {
ydir = SIGNUM(y2 - y1);
fx = (double)x;
fy1 = (double)y1;
fy2 = (double)y2;
draw(fx, fy1, fx, fy2);
for (i = 0; i < 10; i++) {
if (i < y1 && i < y2)
continue;
if (i >= y1 && i >= y2)
continue;
CHECK_CELL(x < 0 ? 0 : x < 10 ? x : 9,
i,
x < 0 ? (double)ydir : x < 10 ? (double)ydir : 0,
(double)ydir);
}
CHECK_ZEROED();
draw(fx + 0.25, fy1, fx + 0.25, fy2);
for (i = 0; i < 10; i++) {
if (i < y1 && i < y2)
continue;
if (i >= y1 && i >= y2)
continue;
CHECK_CELL(x < 0 ? 0 : x < 10 ? x : 9,
i,
x < 0 ? (double)ydir : x < 10 ? (double)ydir * 0.75 : 0,
(double)ydir);
}
CHECK_ZEROED();
draw(fx, fy1 + 0.25, fx, fy2 + 0.25);
for (i = 0; i < 10; i++) {
if (i < y1 + (y1 > y2) && i < y2 + (y2 > y1))
continue;
if (i >= y1 + (y1 > y2) && i >= y2 + (y2 > y1))
continue;
yhit = i == (y1 < y2 ? y1 : y2) ? 0.75 : i == (y2 > y1 ? y2 : y1) ? 0.25 : 1.0;
CHECK_CELL(x < 0 ? 0 : x < 10 ? x : 9,
i,
(x < 0 ? (double)ydir : x < 10 ? (double)ydir : 0) * yhit,
(double)ydir * yhit);
}
CHECK_ZEROED();
draw(fx + 0.25, fy1 + 0.25, fx + 0.25, fy2 + 0.25);
for (i = 0; i < 10; i++) {
if (i < y1 + (y1 > y2) && i < y2 + (y2 > y1))
continue;
if (i >= y1 + (y1 > y2) && i >= y2 + (y2 > y1))
continue;
yhit = i == (y1 < y2 ? y1 : y2) ? 0.75 : i == (y2 > y1 ? y2 : y1) ? 0.25 : 1.0;
CHECK_CELL(x < 0 ? 0 : x < 10 ? x : 9,
i,
(x < 0 ? (double)ydir : x < 10 ? (double)ydir * 0.75 : 0) * yhit,
(double)ydir * yhit);
}
CHECK_ZEROED();
}
}
}
}
static void
check_11_diagonal(ssize_t xshift, ssize_t yshift)
{
ssize_t i1, i2, i, imin, imax, x, y;
double fi1, fi2, idir;
double fxshift = (double)xshift;
double fyshift = (double)yshift;
for (i1 = -2; i1 < 12; i1++) {
for (i2 = -2; i2 < 12; i2++) {
imin = i1 < i2 ? i1 : i2;
imax = i1 > i2 ? i1 : i2;
fi1 = (double)i1;
fi2 = (double)i2;
idir = SIGNUM(fi2 - fi1);
draw(fi1 + fxshift, fi1 + fyshift, fi2 + fxshift, fi2 + fyshift);
for (i = imin; i < imax; i++) {
x = i + xshift;
y = i + yshift;
if (y < 0 || y >= 10)
continue;
if (x >= 10)
CHECK_CELL(9, y, 0, idir);
else if (x >= 0)
CHECK_CELL(x, y, 0.5 * idir, idir);
else
CHECK_CELL(0, y, idir, idir);
}
CHECK_ZEROED();
draw(fi1 + fxshift, fi1 + fyshift + 0.25, fi2 + fxshift, fi2 + fyshift + 0.25);
for (i = imin; i < imax; i++) {
x = i + xshift;
y = i + yshift;
if (0 <= y && y < 10) {
if (x >= 10)
SUB_CELL(9, y, 0, idir * 0.75);
else if (x >= 0)
SUB_CELL(x, y, idir * (0.75 - 0.75 * 0.75 / 2), idir * 0.75);
else
SUB_CELL(0, y, idir * 0.75, idir * 0.75);
}
x = i + xshift;
y = i + yshift + 1;
if (0 <= y && y < 10) {
if (x >= 10)
SUB_CELL(9, y, 0, idir * 0.25);
else if (x >= 0)
SUB_CELL(x, y, idir * 0.25 * 0.25 / 2, idir * 0.25);
else
SUB_CELL(0, y, idir * 0.25, idir * 0.25);
}
}
CHECK_ZEROED();
draw(fi1 + fxshift + 0.25, fi1 + fyshift, fi2 + fxshift + 0.25, fi2 + fyshift);
for (i = imin; i < imax; i++) {
x = i + xshift;
y = i + yshift;
if (0 <= y && y < 10) {
if (x >= 10)
SUB_CELL(9, y, 0, idir * 0.75);
else if (x >= 0)
SUB_CELL(x, y, idir * 0.75 * 0.75 / 2, idir * 0.75);
else
SUB_CELL(0, y, idir * 0.75, idir * 0.75);
}
x = i + xshift + 1;
y = i + yshift;
if (0 <= y && y < 10) {
if (x >= 10)
SUB_CELL(9, y, 0, idir * 0.25);
else if (x >= 0)
SUB_CELL(x, y, idir * (0.25 - 0.25 * 0.25 / 2), idir * 0.25);
else
SUB_CELL(0, y, idir * 0.25, idir * 0.25);
}
}
CHECK_ZEROED();
draw(fi1 + fxshift + 0.25, fi1 + fyshift + 0.25, fi2 + fxshift + 0.25, fi2 + fyshift + 0.25);
for (i = imin; i < imax; i++) {
x = i + xshift;
y = i + yshift;
if (0 <= y && y < 10) {
if (x >= 10)
SUB_CELL(9, y, 0, idir * 0.75);
else if (x >= 0)
SUB_CELL(x, y, idir * 0.75 * 0.75 / 2, idir * 0.75);
else
SUB_CELL(0, y, idir * 0.75, idir * 0.75);
}
x = i + xshift + 1;
y = i + yshift + 1;
if (0 <= y && y < 10) {
if (x >= 10)
SUB_CELL(9, y, 0, idir * 0.25);
else if (x >= 0)
SUB_CELL(x, y, idir * (0.25 - 0.25 * 0.25 / 2), idir * 0.25);
else
SUB_CELL(0, y, idir * 0.25, idir * 0.25);
}
}
CHECK_ZEROED();
}
}
}
static void
check_11_antidiagonal(ssize_t xshift, ssize_t yshift)
{
ssize_t i1, i2, i, imin, imax, x, y;
double fi1, fi2, ydir;
double fxshift = (double)xshift;
double fyshift = (double)yshift;
for (i1 = -2; i1 < 12; i1++) {
for (i2 = -2; i2 < 12; i2++) {
imin = i1 < i2 ? i1 : i2;
imax = i1 > i2 ? i1 : i2;
fi1 = (double)i1;
fi2 = (double)i2;
ydir = -SIGNUM(fi2 - fi1);
draw(fi1 + fxshift, fi2 + fyshift, fi2 + fxshift, fi1 + fyshift);
for (i = imin; i < imax; i++) {
x = i + xshift;
y = imax - i + imin + yshift - 1;
if (y < 0 || y >= 10)
continue;
if (x >= 10)
CHECK_CELL(9, y, 0, ydir);
else if (x >= 0)
CHECK_CELL(x, y, 0.5 * ydir, ydir);
else
CHECK_CELL(0, y, ydir, ydir);
}
CHECK_ZEROED();
draw(fi1 + fxshift, fi2 + fyshift + 0.25, fi2 + fxshift, fi1 + fyshift + 0.25);
for (i = imin; i < imax; i++) {
x = i + xshift;
y = imax - i + imin + yshift - 1;
if (0 <= y && y < 10) {
if (x >= 10)
SUB_CELL(9, y, 0, ydir * 0.75);
else if (x >= 0)
SUB_CELL(x, y, ydir * 0.75 * 0.75 / 2, ydir * 0.75);
else
SUB_CELL(0, y, ydir * 0.75, ydir * 0.75);
}
x = i + xshift;
y = imax - i + imin + yshift - 1 + 1;
if (0 <= y && y < 10) {
if (x >= 10)
SUB_CELL(9, y, 0, ydir * 0.25);
else if (x >= 0)
SUB_CELL(x, y, ydir * (0.25 - 0.25 * 0.25 / 2), ydir * 0.25);
else
SUB_CELL(0, y, ydir * 0.25, ydir * 0.25);
}
}
CHECK_ZEROED();
draw(fi1 + fxshift + 0.25, fi2 + fyshift, fi2 + fxshift + 0.25, fi1 + fyshift);
for (i = imin; i < imax; i++) {
x = i + xshift;
y = imax - i + imin + yshift - 1;
if (0 <= y && y < 10) {
if (x >= 10)
SUB_CELL(9, y, 0, ydir * 0.75);
else if (x >= 0)
SUB_CELL(x, y, ydir * 0.75 * 0.75 / 2, ydir * 0.75);
else
SUB_CELL(0, y, ydir * 0.75, ydir * 0.75);
}
x = i + xshift + 1;
y = imax - i + imin + yshift - 1;
if (0 <= y && y < 10) {
if (x >= 10)
SUB_CELL(9, y, 0, ydir * 0.25);
else if (x >= 0)
SUB_CELL(x, y, ydir * (0.25 - 0.25 * 0.25 / 2), ydir * 0.25);
else
SUB_CELL(0, y, ydir * 0.25, ydir * 0.25);
}
}
CHECK_ZEROED();
draw(fi1 + fxshift + 0.25, fi2 + fyshift + 0.25, fi2 + fxshift + 0.25, fi1 + fyshift + 0.25);
for (i = imin; i < imax; i++) {
x = i + xshift;
y = imax - i + imin + yshift - 1;
if (0 <= y && y < 10) {
if (x >= 10)
SUB_CELL(9, y, 0, ydir * 0.5);
else if (x >= 0)
SUB_CELL(x, y, ydir / 8.0, ydir * 0.5);
else
SUB_CELL(0, y, ydir * 0.5, ydir * 0.5);
}
x = i + xshift;
y = imax - i + imin + yshift - 1 + 1;
if (0 <= y && y < 10) {
if (x >= 10)
SUB_CELL(9, y, 0, ydir * 0.25);
else if (x >= 0)
SUB_CELL(x, y, ydir * (2.5 / 16.0), ydir * 0.25);
else
SUB_CELL(0, y, ydir * 0.25, ydir * 0.25);
}
x = i + xshift + 1;
y = imax - i + imin + yshift - 1;
if (0 <= y && y < 10) {
if (x >= 10)
SUB_CELL(9, y, 0, ydir * 0.25);
else if (x >= 0)
SUB_CELL(x, y, ydir * (3.5 / 16.0), ydir * 0.25);
else
SUB_CELL(0, y, ydir * 0.25, ydir * 0.25);
}
}
CHECK_ZEROED();
}
}
}
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
try_random_line(double xmin, double xmax, double ymin, double ymax)
{
double x1 = frand(xmin, xmax);
double y1 = frand(ymin, ymax);
double x2 = frand(xmin, xmax);
double y2 = frand(ymin, ymax);
draw(x1, y1, x2, y2);
}
static void
check_random_line(double rand_xmin, double rand_xmax, double rand_ymin, double rand_ymax)
{
double x1 = frand(rand_xmin, rand_xmax);
double y1 = frand(rand_ymin, rand_ymax);
double x2 = frand(rand_xmin, rand_xmax);
double y2 = frand(rand_ymin, rand_ymax);
double dx = x2 - x1;
double dy = y2 - y1;
double ymin = fmin(y1, y2);
double ymax = fmax(y1, y2);
double xmin = fmin(x1, x2);
double xmax = fmax(x1, x2);
double ydir = SIGNUM(y2 - y1);
double inklevel, area, expected_area, inkdir, dist;
size_t y, x, jumps;
int hit;
draw(x1, y1, x2, y2);
area = 0;
for (y = 0; y < 10; y++) {
inklevel = 0;
jumps = 0;
inkdir = 0;
for (x = 0; x < 10; x++) {
area += raster->cells[y * 10 + x].cell_coverage + inklevel;
inklevel += raster->cells[y * 10 + x].opposite_coverage;
if (inkdir == 0) {
inkdir = raster->cells[y * 10 + x].opposite_coverage;
} else if (inkdir < 0 && raster->cells[y * 10 + x].opposite_coverage > 0) {
ASSERT(!jumps);
jumps += 1;
inkdir = raster->cells[y * 10 + x].opposite_coverage;
} else if (inkdir > 0 && raster->cells[y * 10 + x].opposite_coverage < 0) {
ASSERT(!jumps);
jumps += 1;
inkdir = raster->cells[y * 10 + x].opposite_coverage;
}
if (!dy)
continue;
hit = x >= floor(xmin) && x <= ceil(xmax);
hit |= !x && xmin < 0;
hit &= y >= floor(ymin) && y <= ceil(ymax);
ASSERT(hit || !raster->cells[y * 10 + x].cell_coverage);
hit = x >= floor(xmin) && x <= ceil(xmax);
hit &= y >= floor(ymin) && y <= ceil(ymax);
if (!hit || !dx || !x)
continue;
dist = ((double)x+0.5)*dy/dx - ((double)y+0.5) + (y1 - x1*dy/dx);
dist /= hypot(dy/dx, 1);
dist = fabs(dist);
if (dist > 0.708)
ASSERT(!raster->cells[y * 10 + x].cell_coverage);
}
if ((double)y < floor(ymin) || (double)y >= ceil(ymax))
ASSERT(inklevel == 0.0);
else if ((double)y == floor(ymin))
ASSERT(eq(inklevel, ydir * (fmin(ymax, (double)y + 1.0) - ymin)));
else if ((double)y == floor(ymax))
ASSERT(eq(inklevel, ydir * fmod(ymax, 1.0)));
else
ASSERT(eq(inklevel, ydir));
}
if (ymax <= 0 || ymin >= 10 || !dy || fmin(x1, x2) >= 10) {
expected_area = 0;
} else if (fmax(x1, x2) < 0) {
expected_area = ydir * 10 * (fmin(10, ymax) - fmax(ymin, 0));
} else if (!dx) {
expected_area = ydir * (10 - fmin(fmax(0, x1), 10)) * (fmin(10, ymax) - fmax(ymin, 0));
} else {
#define y(x) (y1 + ((x) - x1) * k)
#define x(y) (x1 + ((y) - y1) / k)
double k = dy / dx;
double a = xmin;
double b = xmax;
double mid;
a = y(a) > 10 ? x(10) : y(a) < 0 ? x(0) : a;
b = y(b) > 10 ? x(10) : y(b) < 0 ? x(0) : b;
b = fmin(b, 10);
a = fmin(a, 10);
mid = (a + b) / 2;
if (mid < 0) {
expected_area = (10 - 0) * fabs(y(b) - y(a));
if (b > 0)
expected_area -= fabs(y(b) - y(0)) * (b - 0) / 2;
} else {
expected_area = (10 - mid) * fabs(y(b) - y(a));
if (a < 0)
expected_area -= fabs(y(0) - y(a)) * (0 - a) / 2;
}
expected_area *= ydir;
#undef y
#undef x
}
ASSERT(eq(area, expected_area));
rtgrpblib_reset_raster(raster, 10, 10);
}
int
main(void)
{
size_t i;
raster = rtgrpblib_create_raster(10, 10);
ASSERT(raster);
srand((unsigned)(uintptr_t)raster);
check_horizontal();
check_vertical();
check_11_diagonal(0, 0);
check_11_diagonal(0, 2);
check_11_diagonal(0, -2);
check_11_diagonal(2, 0);
check_11_diagonal(2, 2);
check_11_diagonal(2, -2);
check_11_diagonal(-2, 0);
check_11_diagonal(-2, 2);
check_11_diagonal(-2, -2);
check_11_antidiagonal(0, 0);
check_11_antidiagonal(0, 2);
check_11_antidiagonal(0, -2);
check_11_antidiagonal(2, 0);
check_11_antidiagonal(2, 2);
check_11_antidiagonal(2, -2);
check_11_antidiagonal(-2, 0);
check_11_antidiagonal(-2, 2);
check_11_antidiagonal(-2, -2);
for (i = 0; i < 100000UL; i++) {
try_random_line(0., 1., 0., 1.);
try_random_line(-5., 15., -5., 15.);
}
rtgrpblib_reset_raster(raster, 10, 10);
for (i = 0; i < 10000UL; i++) {
check_random_line(0., 1., 0., 1.);
check_random_line(-5., 15., -5., 15.);
}
free(raster);
return 0;
}
#endif
