c++calgorithmbezierbresenham
Bresenham's Algorithm Anti-aliased quadratic Bézier curve with line thickness
In the document found at http://members.chello.at/easyfilter/bresenham.html it walks through the process of creating curves based arround the Bresenham's Algorithm, and it ends with an algorithm to create anti-aliased thick lines, but how can this be aplied to the Quadratic Bezier Curves (the Anti Aliased version found on the same site)?
I tried change the plotQuadBezierSegAA function's last line to use the algorithm of the thick line, but obviously it did not worked as it is computing the other pixels one by one. (I changed from plotLineAA(x0,y0, x2,y2); to plotLineWidth(x0, y0, x2, y2, wd);)
I also tried to draw more curves slightly shifted until it had the wanted thickness, but it creates problems with the anti aliasing colors. (A for loop that shifted by the x and y step (xx and yy variables) and recursively called the plotQuadBezierSegWidth).
None of this tries actualy worked, so could please someone help me acomplish the thickness in this curves. (The algorithm so far is the one from plotQuadBezierSegAA found on that site).
Code for the shifting:
void plotQuadBezierSegAA(int x0, int y0, int x1, int y1, int x2, int y2, int wd)
{
int sx = x2-x1, sy = y2-y1;
long xx = x0-x1, yy = y0-y1, xy; /* relative values for checks */
double dx, dy, err, ed, cur = xx*sy-yy*sx; /* curvature */
assert(xx*sx >= 0 && yy*sy >= 0); /* sign of gradient must not change */
if (sx*(long)sx+sy*(long)sy > xx*xx+yy*yy) { /* begin with longer part */
x2 = x0; x0 = sx+x1; y2 = y0; y0 = sy+y1; cur = -cur; /* swap P0 P2 */
}
if (cur != 0)
{ /* no straight line */
xx += sx; xx *= sx = x0 < x2 ? 1 : -1; /* x step direction */
yy += sy; yy *= sy = y0 < y2 ? 1 : -1; /* y step direction */
// SHIFTING HERE
plotQuadBezierSegAA(x0 + xx, y0 + yy, x1 + xx, y1 + yy, x2 + xx, y2 + yy, wd - 1);
xy = 2*xx*yy; xx *= xx; yy *= yy; /* differences 2nd degree */
if (cur*sx*sy < 0) { /* negated curvature? */
xx = -xx; yy = -yy; xy = -xy; cur = -cur;
}
dx = 4.0*sy*(x1-x0)*cur+xx-xy; /* differences 1st degree */
dy = 4.0*sx*(y0-y1)*cur+yy-xy;
xx += xx; yy += yy; err = dx+dy+xy; /* error 1st step */
do {
cur = fmin(dx+xy,-xy-dy);
ed = fmax(dx+xy,-xy-dy); /* approximate error distance */
ed = 255/(ed+2*ed*cur*cur/(4.*ed*ed+cur*cur));
setPixelAA(x0,y0, ed*fabs(err-dx-dy-xy)); /* plot curve */
if (x0 == x2 && y0 == y2) return;/* last pixel -> curve finished */
x1 = x0; cur = dx-err; y1 = 2*err+dy < 0;
if (2*err+dx > 0) { /* x step */
if (err-dy < ed) setPixelAA(x0,y0+sy, ed*fabs(err-dy));
x0 += sx; dx -= xy; err += dy += yy;
}
if (y1) { /* y step */
if (cur < ed) setPixelAA(x1+sx,y0, ed*fabs(cur));
y0 += sy; dy -= xy; err += dx += xx;
}
} while (dy < dx); /* gradient negates -> close curves */
}
plotLineAA(x0,y0, x2,y2); /* plot remaining needle to end */
}
Solution
In the scratchpad (js file here: http://members.chello.at/easyfilter/bresenham.js lines 909 to 1049) it uses the plotQuadRationalBezierWidth with the w set as 1 to compute a normal Bezier curve with a specified amount of width for the thickness of the line (the wd parameter).
The ported c code from the file:
void plotQuadRationalBezierWidthSeg(
int x0, int y0, int x1, int y1, int x2, int y2, int w, float th) { /* plot a limited rational Bezier segment
of thickness th, squared weight */
int sx = x2 - x1, sy = y2 - y1; /* relative values for checks */
int dx = x0 - x2, dy = y0 - y2, xx = x0 - x1, yy = y0 - y1;
double xy = xx * sy + yy * sx, cur = xx * sy - yy * sx, err, e2,
ed; /* curvature */
bool fullBreak = false;
if (cur != 0.0 && w > 0.0) { /* no straight line */
if (sx * sx + sy * sy >
xx * xx + yy * yy) { /* begin with longer part */
x2 = x0;
x0 -= dx;
y2 = y0;
y0 -= dy;
cur = -cur; /* swap P0 P2 */
}
xx = 2.0 * (4.0 * w * sx * xx + dx * dx); /* differences 2nd degree */
yy = 2.0 * (4.0 * w * sy * yy + dy * dy);
sx = x0 < x2 ? 1 : -1; /* x step direction */
sy = y0 < y2 ? 1 : -1; /* y step direction */
xy = -2.0 * sx * sy * (2.0 * w * xy + dx * dy);
if (cur * sx * sy < 0) { /* negated curvature? */
xx = -xx;
yy = -yy;
cur = -cur;
xy = -xy;
}
dx = 4.0 * w * (x1 - x0) * sy * cur +
xx / 2.0; /* differences 1st degree */
dy = 4.0 * w * (y0 - y1) * sx * cur + yy / 2.0;
if (w < 0.5 &&
(dx + xx <= 0 || dy + yy >= 0)) { /* flat ellipse, algo fails */
cur = (w + 1.0) / 2.0;
w = std::sqrtf(w);
xy = 1.0 / (w + 1.0);
sx = std::floor((x0 + 2.0 * w * x1 + x2) * xy / 2.0 +
0.5); /* subdivide curve */
sy = std::floor((y0 + 2.0 * w * y1 + y2) * xy / 2.0 +
0.5); /* plot separately */
dx = std::floor((w * x1 + x0) * xy + 0.5);
dy = std::floor((y1 * w + y0) * xy + 0.5);
plotQuadRationalBezierWidthSeg(x0, y0, dx, dy, sx, sy, cur, th);
dx = std::floor((w * x1 + x2) * xy + 0.5);
dy = std::floor((y1 * w + y2) * xy + 0.5);
return plotQuadRationalBezierWidthSeg(sx, sy, dx, dy, x2, y2, cur,
th);
}
for (err = 0;
dy + 2 * yy < 0 && dx + 2 * xx > 0;) /* loop of steep/flat curve */
if (dx + dy + xy < 0) { /* steep curve */
do {
ed = -dy -
2 * dy * dx * dx /
(4. * dy * dy + dx * dx); /* approximate sqrt */
w = (th - 1) * ed; /* scale line width */
x1 = std::floor((err - ed - w / 2) / dy); /* start offset */
e2 = err - x1 * dy - w / 2; /* error value at offset */
x1 = x0 - x1 * sx; /* start point */
setPixel(x1, y0, 255 * e2 / ed); /* aliasing pre-pixel */
for (e2 = -w - dy - e2; e2 - dy < ed; e2 -= dy)
setPixel(x1 += sx, y0); /* pixel on thick line */
setPixel(x1 + sx, y0,
255 * e2 / ed); /* aliasing post-pixel */
if (y0 == y2) return; /* last pixel -> curve finished */
y0 += sy;
dy += xy;
err += dx;
dx += xx; /* y step */
if (2 * err + dy > 0) { /* e_x+e_xy > 0 */
x0 += sx;
dx += xy;
err += dy;
dy += yy; /* x step */
}
if (x0 != x2 && (dx + 2 * xx <= 0 || dy + 2 * yy >= 0))
if (std::abs(y2 - y0) > std::abs(x2 - x0)) {
fullBreak = true;
break;
} else
break; /* other curve near */
} while (dx + dy + xy < 0); /* gradient still steep? */
if (fullBreak) break;
/* change from steep to flat curve */
for (cur = err - dy - w / 2, y1 = y0; cur < ed;
y1 += sy, cur += dx) {
for (e2 = cur, x1 = x0; e2 - dy < ed; e2 -= dy)
setPixel(x1 -= sx, y1); /* pixel on thick line */
setPixel(x1 - sx, y1,
255 * e2 / ed); /* aliasing post-pixel */
}
} else { /* flat curve */
do {
ed = dx +
2 * dx * dy * dy /
(4. * dx * dx + dy * dy); /* approximate sqrt */
w = (th - 1) * ed; /* scale line width */
y1 = std::floor((err + ed + w / 2) / dx); /* start offset */
e2 = y1 * dx - w / 2 - err; /* error value at offset */
y1 = y0 - y1 * sy; /* start point */
setPixel(x0, y1, 255 * e2 / ed); /* aliasing pre-pixel */
for (e2 = dx - e2 - w; e2 + dx < ed; e2 += dx)
setPixel(x0, y1 += sy); /* pixel on thick line */
setPixel(x0, y1 + sy,
255 * e2 / ed); /* aliasing post-pixel */
if (x0 == x2) return; /* last pixel -> curve finished */
x0 += sx;
dx += xy;
err += dy;
dy += yy; /* x step */
if (2 * err + dx < 0) { /* e_y+e_xy < 0 */
y0 += sy;
dy += xy;
err += dx;
dx += xx; /* y step */
}
if (y0 != y2 && (dx + 2 * xx <= 0 || dy + 2 * yy >= 0))
if (std::abs(y2 - y0) <= std::abs(x2 - x0)) {
fullBreak = true;
break;
} else
break; /* other curve near */
} while (dx + dy + xy >= 0); /* gradient still flat? */
if (fullBreak) break;
/* change from flat to steep curve */
for (cur = -err + dx - w / 2, x1 = x0; cur < ed;
x1 += sx, cur -= dy) {
for (e2 = cur, y1 = y0; e2 + dx < ed; e2 += dx)
setPixel(x1, y1 -= sy); /* pixel on thick line */
setPixel(x1, y1 - sy,
255 * e2 / ed); /* aliasing post-pixel */
}
}
}
setLine(x0, y0, x2, y2, th); /* confusing error values */
}
void plotQuadRationalBezierWidth(
int x0, int y0, int x1, int y1, int x2, int y2, int w,
float th) { /* plot any anti-aliased quadratic rational Bezier curve */
int x = x0 - 2 * x1 + x2, y = y0 - 2 * y1 + y2;
double xx = x0 - x1, yy = y0 - y1, ww, t, q;
if (xx * (x2 - x1) > 0) { /* horizontal cut at P4? */
if (yy * (y2 - y1) > 0) /* vertical cut at P6 too? */
if (std::abs(xx * y) > std::abs(yy * x)) { /* which first? */
x0 = x2;
x2 = xx + x1;
y0 = y2;
y2 = yy + y1; /* swap points */
} /* now horizontal cut at P4 comes first */
if (x0 == x2 || w == 1.0)
t = (x0 - x1) / x;
else { /* non-rational or rational case */
q = std::sqrt(4.0 * w * w * (x0 - x1) * (x2 - x1) +
(x2 - x0) * (x2 - x0));
if (x1 < x0) q = -q;
t = (2.0 * w * (x0 - x1) - x0 + x2 + q) /
(2.0 * (1.0 - w) * (x2 - x0)); /* t at P4 */
}
q = 1.0 / (2.0 * t * (1.0 - t) * (w - 1.0) + 1.0); /* sub-divide at t */
xx = (t * t * (x0 - 2.0 * w * x1 + x2) + 2.0 * t * (w * x1 - x0) + x0) *
q; /* = P4 */
yy = (t * t * (y0 - 2.0 * w * y1 + y2) + 2.0 * t * (w * y1 - y0) + y0) *
q;
ww = t * (w - 1.0) + 1.0;
ww *= ww * q; /* squared weight P3 */
w = ((1.0 - t) * (w - 1.0) + 1.0) * std::sqrt(q); /* weight P8 */
x = std::floor(xx + 0.5);
y = std::floor(yy + 0.5); /* P4 */
yy = (xx - x0) * (y1 - y0) / (x1 - x0) + y0; /* intersect P3 | P0 P1 */
plotQuadRationalBezierWidthSeg(x0, y0, x, std::floor(yy + 0.5), x, y,
ww, th);
yy = (xx - x2) * (y1 - y2) / (x1 - x2) + y2; /* intersect P4 | P1 P2 */
y1 = std::floor(yy + 0.5);
x0 = x1 = x;
y0 = y; /* P0 = P4, P1 = P8 */
}
if ((y0 - y1) * (y2 - y1) > 0) { /* vertical cut at P6? */
if (y0 == y2 || w == 1.0)
t = (y0 - y1) / (y0 - 2.0 * y1 + y2);
else { /* non-rational or rational case */
q = std::sqrt(4.0 * w * w * (y0 - y1) * (y2 - y1) +
(y2 - y0) * (y2 - y0));
if (y1 < y0) q = -q;
t = (2.0 * w * (y0 - y1) - y0 + y2 + q) /
(2.0 * (1.0 - w) * (y2 - y0)); /* t at P6 */
}
q = 1.0 / (2.0 * t * (1.0 - t) * (w - 1.0) + 1.0); /* sub-divide at t */
xx = (t * t * (x0 - 2.0 * w * x1 + x2) + 2.0 * t * (w * x1 - x0) + x0) *
q; /* = P6 */
yy = (t * t * (y0 - 2.0 * w * y1 + y2) + 2.0 * t * (w * y1 - y0) + y0) *
q;
ww = t * (w - 1.0) + 1.0;
ww *= ww * q; /* squared weight P5 */
w = ((1.0 - t) * (w - 1.0) + 1.0) * std::sqrt(q); /* weight P7 */
x = std::floor(xx + 0.5);
y = std::floor(yy + 0.5); /* P6 */
xx = (x1 - x0) * (yy - y0) / (y1 - y0) + x0; /* intersect P6 | P0 P1 */
plotQuadRationalBezierWidthSeg(x0, y0, std::floor(xx + 0.5), y, x, y,
ww, th);
xx = (x1 - x2) * (yy - y2) / (y1 - y2) + x2; /* intersect P7 | P1 P2 */
x1 = std::floor(xx + 0.5);
x0 = x;
y0 = y1 = y; /* P0 = P6, P1 = P7 */
}
plotQuadRationalBezierWidthSeg(x0, y0, x1, y1, x2, y2, w * w, th);
}
void plotQuadBezierWidth(int x0, int y0, int x1, int y1, int x2, int y2, int wd) {
plotQuadRationalBezierWidth(x0, y0, x1, y1, x2, y2, 1, wd);
}
Example (Thickness 5px):
plotQuadBezierWidth(10, 10, 10, 300, 100, 500, 5);
Example (Thickness 1px):
plotQuadBezierWidth(10, 10, 10, 300, 100, 500, 1);
- Is this declaration UB?
- return type defaults to 'int' [-Wimplicit-int]
- What is the scope of `fesetround()`?
- Which specific optimization flag causes libm functions to be treated as pure?
- How to render text in SDL2?
- Why would you use 'extern "C++"'?
- Strange Behavior Compiler Ignoring NULL Check Unless I Print Something in the if Statement
- Fast inverse square root using fixed point instead of floating point
- What is the const qualifier attached to in C: the memory area or the pointer?
- GCC options for strictest C code?
- How to do an explicit fall-through in C
- How do compilers treat CONST qualifier when the pointer points to a memory location obtained with malloc()?
- C: cmocka headers - how to unittest?
- Why in C when I print a double with a one decimal it round it to the next number
- Android C to Java SWIG unable to compile: incompatible types: byte cannot be converted to SWIGTYPE_p_uint8_t
- GNU Make in Ubuntu giving fatal error: rpc/types.h: No such file or directory
- How can I exclude non-numeric keys? CS50 Caesar Pset2
- How change every struct in an array of pointers?
- Optimized 2x2 matrix multiplication: Slow assembly versus fast SIMD
- Simple frame by frame video decoder library
- GCC no longer implements <varargs.h>
- Contents of IO buffer unknown == unsafe?
- Avoiding strcpy overflow destination warning
- Sort program not working, not sure why
- Fast & accurate atan/arctan approximation algorithm
- What's the difference between strtok_r and strtok_s in C?
- How memory address for pointer to arrays is same as an element in 2D array?
- Which is the best way to suppress "unused variable" warning
- How to use ellipsis in c's case statement?
- Fast ceiling of an integer division in C / C++