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);
