逐像素贝塞尔曲线
Pixel by pixel Bézier Curve
我通过 google 找到的 quadratic/cubic 贝塞尔曲线代码主要通过将线细分为一系列点并用直线连接它们来工作。光栅化发生在直线算法中，而不是在贝塞尔算法中。像 Bresenham 的算法一样逐个像素地栅格化一条线，并且可以优化（参见 Po-Han Lin's solution）。
什么是二次贝塞尔曲线算法，它像直线算法一样逐像素工作，而不是通过绘制一系列点？
您可以使用 De Casteljau's algorithm 将曲线细分为足够多的片段，每个片段都是一个像素。
这是用于在间隔 T 处的二次曲线上找到 [x,y] 点的方程式：
// Given 3 control points defining the Quadratic curve
// and given T which is an interval between 0.00 and 1.00 along the curve.
// Note:
// At the curve's starting control point T==0.00.
// At the curve's ending control point T==1.00.
var x = Math.pow(1-T,2)*startPt.x + 2 * (1-T) * T * controlPt.x + Math.pow(T,2) * endPt.x;
var y = Math.pow(1-T,2)*startPt.y + 2 * (1-T) * T * controlPt.y + Math.pow(T,2) * endPt.y;
要实际使用此方程式，您可以输入大约 1000 个介于 0.00 和 1.00 之间的 T 值。这导致一组 1000 个点保证沿着二次曲线。
沿曲线计算 1000 个点可能是过采样（一些计算的点将处于相同的像素坐标），因此您需要对 1000 个点进行去重，直到该集表示沿曲线的唯一像素坐标。
三次贝塞尔曲线也有类似的方程。
下面是将二次曲线绘制为一组计算像素的示例代码：
var canvas=document.getElementById("canvas");
var ctx=canvas.getContext("2d");
var points=[];
var lastX=-100;
var lastY=-100;
var startPt={x:50,y:200};
var controlPt={x:150,y:25};
var endPt={x:250,y:100};
for(var t=0;t<1000;t++){
var xyAtT=getQuadraticBezierXYatT(startPt,controlPt,endPt,t/1000);
var x=parseInt(xyAtT.x);
var y=parseInt(xyAtT.y);
if(!(x==lastX && y==lastY)){
points.push(xyAtT);
lastX=x;
lastY=y;
}
}
$('#curve').text('Quadratic Curve made up of '+points.length+' individual points');
ctx.fillStyle='red';
for(var i=0;i<points.length;i++){
var x=points[i].x;
var y=points[i].y;
ctx.fillRect(x,y,1,1);
}
function getQuadraticBezierXYatT(startPt,controlPt,endPt,T) {
var x = Math.pow(1-T,2) * startPt.x + 2 * (1-T) * T * controlPt.x + Math.pow(T,2) * endPt.x;
var y = Math.pow(1-T,2) * startPt.y + 2 * (1-T) * T * controlPt.y + Math.pow(T,2) * endPt.y;
return( {x:x,y:y} );
}
body{ background-color: ivory; }
#canvas{border:1px solid red; margin:0 auto; }
<script src="https://ajax.googleapis.com/ajax/libs/jquery/1.9.1/jquery.min.js"></script>
<h4 id='curve'>Q</h4>
<canvas id="canvas" width=350 height=300></canvas>
Bresenham 算法的一个变体适用于二次函数，如圆、椭圆和抛物线，因此它也适用于二次贝塞尔曲线。
我本来打算尝试一个实现，但后来我在网上找到了一个：http://members.chello.at/~easyfilter/bresenham.html。
如果您需要更多详细信息或其他示例，上面提到的页面有一个 link 到 100 页的 PDF，详细说明了该方法：http://members.chello.at/~easyfilter/Bresenham.pdf.
这是 Alois Zingl 网站上用于绘制任何二次贝塞尔曲线的代码。第一个例程在水平和垂直梯度变化处细分曲线：
void plotQuadBezier(int x0, int y0, int x1, int y1, int x2, int y2)
{ /* plot any quadratic Bezier curve */
int x = x0-x1, y = y0-y1;
double t = x0-2*x1+x2, r;
if ((long)x*(x2-x1) > 0) { /* horizontal cut at P4? */
if ((long)y*(y2-y1) > 0) /* vertical cut at P6 too? */
if (fabs((y0-2*y1+y2)/t*x) > abs(y)) { /* which first? */
x0 = x2; x2 = x+x1; y0 = y2; y2 = y+y1; /* swap points */
} /* now horizontal cut at P4 comes first */
t = (x0-x1)/t;
r = (1-t)*((1-t)*y0+2.0*t*y1)+t*t*y2; /* By(t=P4) */
t = (x0*x2-x1*x1)*t/(x0-x1); /* gradient dP4/dx=0 */
x = floor(t+0.5); y = floor(r+0.5);
r = (y1-y0)*(t-x0)/(x1-x0)+y0; /* intersect P3 | P0 P1 */
plotQuadBezierSeg(x0,y0, x,floor(r+0.5), x,y);
r = (y1-y2)*(t-x2)/(x1-x2)+y2; /* intersect P4 | P1 P2 */
x0 = x1 = x; y0 = y; y1 = floor(r+0.5); /* P0 = P4, P1 = P8 */
}
if ((long)(y0-y1)*(y2-y1) > 0) { /* vertical cut at P6? */
t = y0-2*y1+y2; t = (y0-y1)/t;
r = (1-t)*((1-t)*x0+2.0*t*x1)+t*t*x2; /* Bx(t=P6) */
t = (y0*y2-y1*y1)*t/(y0-y1); /* gradient dP6/dy=0 */
x = floor(r+0.5); y = floor(t+0.5);
r = (x1-x0)*(t-y0)/(y1-y0)+x0; /* intersect P6 | P0 P1 */
plotQuadBezierSeg(x0,y0, floor(r+0.5),y, x,y);
r = (x1-x2)*(t-y2)/(y1-y2)+x2; /* intersect P7 | P1 P2 */
x0 = x; x1 = floor(r+0.5); y0 = y1 = y; /* P0 = P6, P1 = P7 */
}
plotQuadBezierSeg(x0,y0, x1,y1, x2,y2); /* remaining part */
}
第二个例程实际上绘制了一段贝塞尔曲线（没有梯度变化）：
void plotQuadBezierSeg(int x0, int y0, int x1, int y1, int x2, int y2)
{ /* plot a limited quadratic Bezier segment */
int sx = x2-x1, sy = y2-y1;
long xx = x0-x1, yy = y0-y1, xy; /* relative values for checks */
double dx, dy, err, 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 */
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*cur*(x1-x0)+xx-xy; /* differences 1st degree */
dy = 4.0*sx*cur*(y0-y1)+yy-xy;
xx += xx; yy += yy; err = dx+dy+xy; /* error 1st step */
do {
setPixel(x0,y0); /* plot curve */
if (x0 == x2 && y0 == y2) return; /* last pixel -> curve finished */
y1 = 2*err < dx; /* save value for test of y step */
if (2*err > dy) { x0 += sx; dx -= xy; err += dy += yy; } /* x step */
if ( y1 ) { y0 += sy; dy -= xy; err += dx += xx; } /* y step */
} while (dy < 0 && dx > 0); /* gradient negates -> algorithm fails */
}
plotLine(x0,y0, x2,y2); /* plot remaining part to end */
}
站点上还提供了抗锯齿代码。
来自 Zingl 网站的三次贝塞尔曲线的对应函数是
void plotCubicBezier(int x0, int y0, int x1, int y1,
int x2, int y2, int x3, int y3)
{ /* plot any cubic Bezier curve */
int n = 0, i = 0;
long xc = x0+x1-x2-x3, xa = xc-4*(x1-x2);
long xb = x0-x1-x2+x3, xd = xb+4*(x1+x2);
long yc = y0+y1-y2-y3, ya = yc-4*(y1-y2);
long yb = y0-y1-y2+y3, yd = yb+4*(y1+y2);
float fx0 = x0, fx1, fx2, fx3, fy0 = y0, fy1, fy2, fy3;
double t1 = xb*xb-xa*xc, t2, t[5];
/* sub-divide curve at gradient sign changes */
if (xa == 0) { /* horizontal */
if (abs(xc) < 2*abs(xb)) t[n++] = xc/(2.0*xb); /* one change */
} else if (t1 > 0.0) { /* two changes */
t2 = sqrt(t1);
t1 = (xb-t2)/xa; if (fabs(t1) < 1.0) t[n++] = t1;
t1 = (xb+t2)/xa; if (fabs(t1) < 1.0) t[n++] = t1;
}
t1 = yb*yb-ya*yc;
if (ya == 0) { /* vertical */
if (abs(yc) < 2*abs(yb)) t[n++] = yc/(2.0*yb); /* one change */
} else if (t1 > 0.0) { /* two changes */
t2 = sqrt(t1);
t1 = (yb-t2)/ya; if (fabs(t1) < 1.0) t[n++] = t1;
t1 = (yb+t2)/ya; if (fabs(t1) < 1.0) t[n++] = t1;
}
for (i = 1; i < n; i++) /* bubble sort of 4 points */
if ((t1 = t[i-1]) > t[i]) { t[i-1] = t[i]; t[i] = t1; i = 0; }
t1 = -1.0; t[n] = 1.0; /* begin / end point */
for (i = 0; i <= n; i++) { /* plot each segment separately */
t2 = t[i]; /* sub-divide at t[i-1], t[i] */
fx1 = (t1*(t1*xb-2*xc)-t2*(t1*(t1*xa-2*xb)+xc)+xd)/8-fx0;
fy1 = (t1*(t1*yb-2*yc)-t2*(t1*(t1*ya-2*yb)+yc)+yd)/8-fy0;
fx2 = (t2*(t2*xb-2*xc)-t1*(t2*(t2*xa-2*xb)+xc)+xd)/8-fx0;
fy2 = (t2*(t2*yb-2*yc)-t1*(t2*(t2*ya-2*yb)+yc)+yd)/8-fy0;
fx0 -= fx3 = (t2*(t2*(3*xb-t2*xa)-3*xc)+xd)/8;
fy0 -= fy3 = (t2*(t2*(3*yb-t2*ya)-3*yc)+yd)/8;
x3 = floor(fx3+0.5); y3 = floor(fy3+0.5); /* scale bounds to int */
if (fx0 != 0.0) { fx1 *= fx0 = (x0-x3)/fx0; fx2 *= fx0; }
if (fy0 != 0.0) { fy1 *= fy0 = (y0-y3)/fy0; fy2 *= fy0; }
if (x0 != x3 || y0 != y3) /* segment t1 - t2 */
plotCubicBezierSeg(x0,y0, x0+fx1,y0+fy1, x0+fx2,y0+fy2, x3,y3);
x0 = x3; y0 = y3; fx0 = fx3; fy0 = fy3; t1 = t2;
}
}
和
void plotCubicBezierSeg(int x0, int y0, float x1, float y1,
float x2, float y2, int x3, int y3)
{ /* plot limited cubic Bezier segment */
int f, fx, fy, leg = 1;
int sx = x0 < x3 ? 1 : -1, sy = y0 < y3 ? 1 : -1; /* step direction */
float xc = -fabs(x0+x1-x2-x3), xa = xc-4*sx*(x1-x2), xb = sx*(x0-x1-x2+x3);
float yc = -fabs(y0+y1-y2-y3), ya = yc-4*sy*(y1-y2), yb = sy*(y0-y1-y2+y3);
double ab, ac, bc, cb, xx, xy, yy, dx, dy, ex, *pxy, EP = 0.01;
/* check for curve restrains */
/* slope P0-P1 == P2-P3 and (P0-P3 == P1-P2 or no slope change) */
assert((x1-x0)*(x2-x3) < EP && ((x3-x0)*(x1-x2) < EP || xb*xb < xa*xc+EP));
assert((y1-y0)*(y2-y3) < EP && ((y3-y0)*(y1-y2) < EP || yb*yb < ya*yc+EP));
if (xa == 0 && ya == 0) { /* quadratic Bezier */
sx = floor((3*x1-x0+1)/2); sy = floor((3*y1-y0+1)/2); /* new midpoint */
return plotQuadBezierSeg(x0,y0, sx,sy, x3,y3);
}
x1 = (x1-x0)*(x1-x0)+(y1-y0)*(y1-y0)+1; /* line lengths */
x2 = (x2-x3)*(x2-x3)+(y2-y3)*(y2-y3)+1;
do { /* loop over both ends */
ab = xa*yb-xb*ya; ac = xa*yc-xc*ya; bc = xb*yc-xc*yb;
ex = ab*(ab+ac-3*bc)+ac*ac; /* P0 part of self-intersection loop? */
f = ex > 0 ? 1 : sqrt(1+1024/x1); /* calculate resolution */
ab *= f; ac *= f; bc *= f; ex *= f*f; /* increase resolution */
xy = 9*(ab+ac+bc)/8; cb = 8*(xa-ya);/* init differences of 1st degree */
dx = 27*(8*ab*(yb*yb-ya*yc)+ex*(ya+2*yb+yc))/64-ya*ya*(xy-ya);
dy = 27*(8*ab*(xb*xb-xa*xc)-ex*(xa+2*xb+xc))/64-xa*xa*(xy+xa);
/* init differences of 2nd degree */
xx = 3*(3*ab*(3*yb*yb-ya*ya-2*ya*yc)-ya*(3*ac*(ya+yb)+ya*cb))/4;
yy = 3*(3*ab*(3*xb*xb-xa*xa-2*xa*xc)-xa*(3*ac*(xa+xb)+xa*cb))/4;
xy = xa*ya*(6*ab+6*ac-3*bc+cb); ac = ya*ya; cb = xa*xa;
xy = 3*(xy+9*f*(cb*yb*yc-xb*xc*ac)-18*xb*yb*ab)/8;
if (ex < 0) { /* negate values if inside self-intersection loop */
dx = -dx; dy = -dy; xx = -xx; yy = -yy; xy = -xy; ac = -ac; cb = -cb;
} /* init differences of 3rd degree */
ab = 6*ya*ac; ac = -6*xa*ac; bc = 6*ya*cb; cb = -6*xa*cb;
dx += xy; ex = dx+dy; dy += xy; /* error of 1st step */
for (pxy = &xy, fx = fy = f; x0 != x3 && y0 != y3; ) {
setPixel(x0,y0); /* plot curve */
do { /* move sub-steps of one pixel */
if (dx > *pxy || dy < *pxy) goto exit; /* confusing values */
y1 = 2*ex-dy; /* save value for test of y step */
if (2*ex >= dx) { /* x sub-step */
fx--; ex += dx += xx; dy += xy += ac; yy += bc; xx += ab;
}
if (y1 <= 0) { /* y sub-step */
fy--; ex += dy += yy; dx += xy += bc; xx += ac; yy += cb;
}
} while (fx > 0 && fy > 0); /* pixel complete? */
if (2*fx <= f) { x0 += sx; fx += f; } /* x step */
if (2*fy <= f) { y0 += sy; fy += f; } /* y step */
if (pxy == &xy && dx < 0 && dy > 0) pxy = &EP;/* pixel ahead valid */
}
exit: xx = x0; x0 = x3; x3 = xx; sx = -sx; xb = -xb; /* swap legs */
yy = y0; y0 = y3; y3 = yy; sy = -sy; yb = -yb; x1 = x2;
} while (leg--); /* try other end */
plotLine(x0,y0, x3,y3); /* remaining part in case of cusp or crunode */
}
正如 Mike 'Pomax' Kamermans 所指出的，网站上的三次贝塞尔曲线的解决方案并不完整；特别是，抗锯齿三次贝塞尔曲线存在问题，有理三次贝塞尔曲线的讨论不完整。
首先我想说，渲染贝塞尔曲线最快最可靠的方法是通过自适应细分用折线逼近它们，然后渲染折线。 @markE 在曲线上绘制许多采样点的方法相当快，但它可以跳过像素。在这里我描述了另一种方法，它最接近线栅格化（尽管它很慢并且难以稳健地实现）。
我通常将曲线参数视为时间。这是伪代码：
- 将光标放在第一个控制点，找到周围的像素。
- 对于像素的每一侧（总共四个），通过求解二次方程检查贝塞尔曲线何时与其线相交。
- 在所有计算出的侧交叉口时间中，选择一个严格在未来发生的交叉路口时间，但要尽可能早。
- 根据最佳边移动到相邻像素。
- 将当前时间设置为该最佳侧路口的时间。
- 从第 2 步开始重复。
此算法一直有效，直到时间参数超过 1。另请注意，曲线恰好接触像素的一侧存在严重问题。我想它可以通过特殊检查解决。
这里是主要代码：
double WhenEquals(double p0, double p1, double p2, double val, double minp) {
//p0 * (1-t)^2 + p1 * 2t(1 - t) + p2 * t^2 = val
double qa = p0 + p2 - 2 * p1;
double qb = p1 - p0;
double qc = p0 - val;
assert(fabs(qa) > EPS); //singular case must be handled separately
double qd = qb * qb - qa * qc;
if (qd < -EPS)
return INF;
qd = sqrt(max(qd, 0.0));
double t1 = (-qb - qd) / qa;
double t2 = (-qb + qd) / qa;
if (t2 < t1) swap(t1, t2);
if (t1 > minp + EPS)
return t1;
else if (t2 > minp + EPS)
return t2;
return INF;
}
void DrawCurve(const Bezier &curve) {
int cell[2];
for (int c = 0; c < 2; c++)
cell[c] = int(floor(curve.pts[0].a[c]));
DrawPixel(cell[0], cell[1]);
double param = 0.0;
while (1) {
int bc = -1, bs = -1;
double bestTime = 1.0;
for (int c = 0; c < 2; c++)
for (int s = 0; s < 2; s++) {
double crit = WhenEquals(
curve.pts[0].a[c],
curve.pts[1].a[c],
curve.pts[2].a[c],
cell[c] + s, param
);
if (crit < bestTime) {
bestTime = crit;
bc = c, bs = s;
}
}
if (bc < 0)
break;
param = bestTime;
cell[bc] += (2*bs - 1);
DrawPixel(cell[0], cell[1]);
}
}
这里要注意的是 "line segments"，当创建得足够小时，相当于像素。贝塞尔曲线不是可线性遍历的曲线，因此我们不能像直线或圆弧那样轻松地一步 "skip ahead to the next pixel"。
当然，您可以对已有的 t 取任意点的切线，然后猜测下一个值 t' 会比哪个像素更远。但是，通常发生的情况是您猜测并猜错了，因为曲线的行为不是线性的，然后您检查 "off" 您的猜测如何，更正您的猜测，然后再次检查。重复直到收敛到下一个像素：这比将曲线展平为大量线段要慢得多，后者是一种快速操作。
如果您选择适合曲线长度的段数，考虑到它的渲染显示，没有人会告诉您曲线变平了。
有很多方法可以重新参数化贝塞尔曲线，但它们很昂贵，而且不同的典型曲线需要不同的重新参数化，所以这也不是真的更快。对于离散显示器最有用的往往是为您的曲线构建一个 LUT（查找 table），其长度适用于显示器上曲线的大小，然后使用该 LUT 作为您的基础绘图、交点检测等数据等等
我通过 google 找到的 quadratic/cubic 贝塞尔曲线代码主要通过将线细分为一系列点并用直线连接它们来工作。光栅化发生在直线算法中，而不是在贝塞尔算法中。像 Bresenham 的算法一样逐个像素地栅格化一条线，并且可以优化（参见 Po-Han Lin's solution）。
什么是二次贝塞尔曲线算法，它像直线算法一样逐像素工作，而不是通过绘制一系列点？
您可以使用 De Casteljau's algorithm 将曲线细分为足够多的片段，每个片段都是一个像素。
这是用于在间隔 T 处的二次曲线上找到 [x,y] 点的方程式：
// Given 3 control points defining the Quadratic curve
// and given T which is an interval between 0.00 and 1.00 along the curve.
// Note:
// At the curve's starting control point T==0.00.
// At the curve's ending control point T==1.00.
var x = Math.pow(1-T,2)*startPt.x + 2 * (1-T) * T * controlPt.x + Math.pow(T,2) * endPt.x;
var y = Math.pow(1-T,2)*startPt.y + 2 * (1-T) * T * controlPt.y + Math.pow(T,2) * endPt.y;
要实际使用此方程式，您可以输入大约 1000 个介于 0.00 和 1.00 之间的 T 值。这导致一组 1000 个点保证沿着二次曲线。
沿曲线计算 1000 个点可能是过采样（一些计算的点将处于相同的像素坐标），因此您需要对 1000 个点进行去重，直到该集表示沿曲线的唯一像素坐标。
三次贝塞尔曲线也有类似的方程。
下面是将二次曲线绘制为一组计算像素的示例代码：
var canvas=document.getElementById("canvas");
var ctx=canvas.getContext("2d");
var points=[];
var lastX=-100;
var lastY=-100;
var startPt={x:50,y:200};
var controlPt={x:150,y:25};
var endPt={x:250,y:100};
for(var t=0;t<1000;t++){
var xyAtT=getQuadraticBezierXYatT(startPt,controlPt,endPt,t/1000);
var x=parseInt(xyAtT.x);
var y=parseInt(xyAtT.y);
if(!(x==lastX && y==lastY)){
points.push(xyAtT);
lastX=x;
lastY=y;
}
}
$('#curve').text('Quadratic Curve made up of '+points.length+' individual points');
ctx.fillStyle='red';
for(var i=0;i<points.length;i++){
var x=points[i].x;
var y=points[i].y;
ctx.fillRect(x,y,1,1);
}
function getQuadraticBezierXYatT(startPt,controlPt,endPt,T) {
var x = Math.pow(1-T,2) * startPt.x + 2 * (1-T) * T * controlPt.x + Math.pow(T,2) * endPt.x;
var y = Math.pow(1-T,2) * startPt.y + 2 * (1-T) * T * controlPt.y + Math.pow(T,2) * endPt.y;
return( {x:x,y:y} );
}
body{ background-color: ivory; }
#canvas{border:1px solid red; margin:0 auto; }
<script src="https://ajax.googleapis.com/ajax/libs/jquery/1.9.1/jquery.min.js"></script>
<h4 id='curve'>Q</h4>
<canvas id="canvas" width=350 height=300></canvas>
Bresenham 算法的一个变体适用于二次函数，如圆、椭圆和抛物线，因此它也适用于二次贝塞尔曲线。
我本来打算尝试一个实现，但后来我在网上找到了一个：http://members.chello.at/~easyfilter/bresenham.html。
如果您需要更多详细信息或其他示例，上面提到的页面有一个 link 到 100 页的 PDF，详细说明了该方法：http://members.chello.at/~easyfilter/Bresenham.pdf.
这是 Alois Zingl 网站上用于绘制任何二次贝塞尔曲线的代码。第一个例程在水平和垂直梯度变化处细分曲线：
void plotQuadBezier(int x0, int y0, int x1, int y1, int x2, int y2)
{ /* plot any quadratic Bezier curve */
int x = x0-x1, y = y0-y1;
double t = x0-2*x1+x2, r;
if ((long)x*(x2-x1) > 0) { /* horizontal cut at P4? */
if ((long)y*(y2-y1) > 0) /* vertical cut at P6 too? */
if (fabs((y0-2*y1+y2)/t*x) > abs(y)) { /* which first? */
x0 = x2; x2 = x+x1; y0 = y2; y2 = y+y1; /* swap points */
} /* now horizontal cut at P4 comes first */
t = (x0-x1)/t;
r = (1-t)*((1-t)*y0+2.0*t*y1)+t*t*y2; /* By(t=P4) */
t = (x0*x2-x1*x1)*t/(x0-x1); /* gradient dP4/dx=0 */
x = floor(t+0.5); y = floor(r+0.5);
r = (y1-y0)*(t-x0)/(x1-x0)+y0; /* intersect P3 | P0 P1 */
plotQuadBezierSeg(x0,y0, x,floor(r+0.5), x,y);
r = (y1-y2)*(t-x2)/(x1-x2)+y2; /* intersect P4 | P1 P2 */
x0 = x1 = x; y0 = y; y1 = floor(r+0.5); /* P0 = P4, P1 = P8 */
}
if ((long)(y0-y1)*(y2-y1) > 0) { /* vertical cut at P6? */
t = y0-2*y1+y2; t = (y0-y1)/t;
r = (1-t)*((1-t)*x0+2.0*t*x1)+t*t*x2; /* Bx(t=P6) */
t = (y0*y2-y1*y1)*t/(y0-y1); /* gradient dP6/dy=0 */
x = floor(r+0.5); y = floor(t+0.5);
r = (x1-x0)*(t-y0)/(y1-y0)+x0; /* intersect P6 | P0 P1 */
plotQuadBezierSeg(x0,y0, floor(r+0.5),y, x,y);
r = (x1-x2)*(t-y2)/(y1-y2)+x2; /* intersect P7 | P1 P2 */
x0 = x; x1 = floor(r+0.5); y0 = y1 = y; /* P0 = P6, P1 = P7 */
}
plotQuadBezierSeg(x0,y0, x1,y1, x2,y2); /* remaining part */
}
第二个例程实际上绘制了一段贝塞尔曲线（没有梯度变化）：
void plotQuadBezierSeg(int x0, int y0, int x1, int y1, int x2, int y2)
{ /* plot a limited quadratic Bezier segment */
int sx = x2-x1, sy = y2-y1;
long xx = x0-x1, yy = y0-y1, xy; /* relative values for checks */
double dx, dy, err, 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 */
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*cur*(x1-x0)+xx-xy; /* differences 1st degree */
dy = 4.0*sx*cur*(y0-y1)+yy-xy;
xx += xx; yy += yy; err = dx+dy+xy; /* error 1st step */
do {
setPixel(x0,y0); /* plot curve */
if (x0 == x2 && y0 == y2) return; /* last pixel -> curve finished */
y1 = 2*err < dx; /* save value for test of y step */
if (2*err > dy) { x0 += sx; dx -= xy; err += dy += yy; } /* x step */
if ( y1 ) { y0 += sy; dy -= xy; err += dx += xx; } /* y step */
} while (dy < 0 && dx > 0); /* gradient negates -> algorithm fails */
}
plotLine(x0,y0, x2,y2); /* plot remaining part to end */
}
站点上还提供了抗锯齿代码。
来自 Zingl 网站的三次贝塞尔曲线的对应函数是
void plotCubicBezier(int x0, int y0, int x1, int y1,
int x2, int y2, int x3, int y3)
{ /* plot any cubic Bezier curve */
int n = 0, i = 0;
long xc = x0+x1-x2-x3, xa = xc-4*(x1-x2);
long xb = x0-x1-x2+x3, xd = xb+4*(x1+x2);
long yc = y0+y1-y2-y3, ya = yc-4*(y1-y2);
long yb = y0-y1-y2+y3, yd = yb+4*(y1+y2);
float fx0 = x0, fx1, fx2, fx3, fy0 = y0, fy1, fy2, fy3;
double t1 = xb*xb-xa*xc, t2, t[5];
/* sub-divide curve at gradient sign changes */
if (xa == 0) { /* horizontal */
if (abs(xc) < 2*abs(xb)) t[n++] = xc/(2.0*xb); /* one change */
} else if (t1 > 0.0) { /* two changes */
t2 = sqrt(t1);
t1 = (xb-t2)/xa; if (fabs(t1) < 1.0) t[n++] = t1;
t1 = (xb+t2)/xa; if (fabs(t1) < 1.0) t[n++] = t1;
}
t1 = yb*yb-ya*yc;
if (ya == 0) { /* vertical */
if (abs(yc) < 2*abs(yb)) t[n++] = yc/(2.0*yb); /* one change */
} else if (t1 > 0.0) { /* two changes */
t2 = sqrt(t1);
t1 = (yb-t2)/ya; if (fabs(t1) < 1.0) t[n++] = t1;
t1 = (yb+t2)/ya; if (fabs(t1) < 1.0) t[n++] = t1;
}
for (i = 1; i < n; i++) /* bubble sort of 4 points */
if ((t1 = t[i-1]) > t[i]) { t[i-1] = t[i]; t[i] = t1; i = 0; }
t1 = -1.0; t[n] = 1.0; /* begin / end point */
for (i = 0; i <= n; i++) { /* plot each segment separately */
t2 = t[i]; /* sub-divide at t[i-1], t[i] */
fx1 = (t1*(t1*xb-2*xc)-t2*(t1*(t1*xa-2*xb)+xc)+xd)/8-fx0;
fy1 = (t1*(t1*yb-2*yc)-t2*(t1*(t1*ya-2*yb)+yc)+yd)/8-fy0;
fx2 = (t2*(t2*xb-2*xc)-t1*(t2*(t2*xa-2*xb)+xc)+xd)/8-fx0;
fy2 = (t2*(t2*yb-2*yc)-t1*(t2*(t2*ya-2*yb)+yc)+yd)/8-fy0;
fx0 -= fx3 = (t2*(t2*(3*xb-t2*xa)-3*xc)+xd)/8;
fy0 -= fy3 = (t2*(t2*(3*yb-t2*ya)-3*yc)+yd)/8;
x3 = floor(fx3+0.5); y3 = floor(fy3+0.5); /* scale bounds to int */
if (fx0 != 0.0) { fx1 *= fx0 = (x0-x3)/fx0; fx2 *= fx0; }
if (fy0 != 0.0) { fy1 *= fy0 = (y0-y3)/fy0; fy2 *= fy0; }
if (x0 != x3 || y0 != y3) /* segment t1 - t2 */
plotCubicBezierSeg(x0,y0, x0+fx1,y0+fy1, x0+fx2,y0+fy2, x3,y3);
x0 = x3; y0 = y3; fx0 = fx3; fy0 = fy3; t1 = t2;
}
}
和
void plotCubicBezierSeg(int x0, int y0, float x1, float y1,
float x2, float y2, int x3, int y3)
{ /* plot limited cubic Bezier segment */
int f, fx, fy, leg = 1;
int sx = x0 < x3 ? 1 : -1, sy = y0 < y3 ? 1 : -1; /* step direction */
float xc = -fabs(x0+x1-x2-x3), xa = xc-4*sx*(x1-x2), xb = sx*(x0-x1-x2+x3);
float yc = -fabs(y0+y1-y2-y3), ya = yc-4*sy*(y1-y2), yb = sy*(y0-y1-y2+y3);
double ab, ac, bc, cb, xx, xy, yy, dx, dy, ex, *pxy, EP = 0.01;
/* check for curve restrains */
/* slope P0-P1 == P2-P3 and (P0-P3 == P1-P2 or no slope change) */
assert((x1-x0)*(x2-x3) < EP && ((x3-x0)*(x1-x2) < EP || xb*xb < xa*xc+EP));
assert((y1-y0)*(y2-y3) < EP && ((y3-y0)*(y1-y2) < EP || yb*yb < ya*yc+EP));
if (xa == 0 && ya == 0) { /* quadratic Bezier */
sx = floor((3*x1-x0+1)/2); sy = floor((3*y1-y0+1)/2); /* new midpoint */
return plotQuadBezierSeg(x0,y0, sx,sy, x3,y3);
}
x1 = (x1-x0)*(x1-x0)+(y1-y0)*(y1-y0)+1; /* line lengths */
x2 = (x2-x3)*(x2-x3)+(y2-y3)*(y2-y3)+1;
do { /* loop over both ends */
ab = xa*yb-xb*ya; ac = xa*yc-xc*ya; bc = xb*yc-xc*yb;
ex = ab*(ab+ac-3*bc)+ac*ac; /* P0 part of self-intersection loop? */
f = ex > 0 ? 1 : sqrt(1+1024/x1); /* calculate resolution */
ab *= f; ac *= f; bc *= f; ex *= f*f; /* increase resolution */
xy = 9*(ab+ac+bc)/8; cb = 8*(xa-ya);/* init differences of 1st degree */
dx = 27*(8*ab*(yb*yb-ya*yc)+ex*(ya+2*yb+yc))/64-ya*ya*(xy-ya);
dy = 27*(8*ab*(xb*xb-xa*xc)-ex*(xa+2*xb+xc))/64-xa*xa*(xy+xa);
/* init differences of 2nd degree */
xx = 3*(3*ab*(3*yb*yb-ya*ya-2*ya*yc)-ya*(3*ac*(ya+yb)+ya*cb))/4;
yy = 3*(3*ab*(3*xb*xb-xa*xa-2*xa*xc)-xa*(3*ac*(xa+xb)+xa*cb))/4;
xy = xa*ya*(6*ab+6*ac-3*bc+cb); ac = ya*ya; cb = xa*xa;
xy = 3*(xy+9*f*(cb*yb*yc-xb*xc*ac)-18*xb*yb*ab)/8;
if (ex < 0) { /* negate values if inside self-intersection loop */
dx = -dx; dy = -dy; xx = -xx; yy = -yy; xy = -xy; ac = -ac; cb = -cb;
} /* init differences of 3rd degree */
ab = 6*ya*ac; ac = -6*xa*ac; bc = 6*ya*cb; cb = -6*xa*cb;
dx += xy; ex = dx+dy; dy += xy; /* error of 1st step */
for (pxy = &xy, fx = fy = f; x0 != x3 && y0 != y3; ) {
setPixel(x0,y0); /* plot curve */
do { /* move sub-steps of one pixel */
if (dx > *pxy || dy < *pxy) goto exit; /* confusing values */
y1 = 2*ex-dy; /* save value for test of y step */
if (2*ex >= dx) { /* x sub-step */
fx--; ex += dx += xx; dy += xy += ac; yy += bc; xx += ab;
}
if (y1 <= 0) { /* y sub-step */
fy--; ex += dy += yy; dx += xy += bc; xx += ac; yy += cb;
}
} while (fx > 0 && fy > 0); /* pixel complete? */
if (2*fx <= f) { x0 += sx; fx += f; } /* x step */
if (2*fy <= f) { y0 += sy; fy += f; } /* y step */
if (pxy == &xy && dx < 0 && dy > 0) pxy = &EP;/* pixel ahead valid */
}
exit: xx = x0; x0 = x3; x3 = xx; sx = -sx; xb = -xb; /* swap legs */
yy = y0; y0 = y3; y3 = yy; sy = -sy; yb = -yb; x1 = x2;
} while (leg--); /* try other end */
plotLine(x0,y0, x3,y3); /* remaining part in case of cusp or crunode */
}
正如 Mike 'Pomax' Kamermans 所指出的，网站上的三次贝塞尔曲线的解决方案并不完整；特别是，抗锯齿三次贝塞尔曲线存在问题，有理三次贝塞尔曲线的讨论不完整。
首先我想说，渲染贝塞尔曲线最快最可靠的方法是通过自适应细分用折线逼近它们，然后渲染折线。 @markE 在曲线上绘制许多采样点的方法相当快，但它可以跳过像素。在这里我描述了另一种方法，它最接近线栅格化（尽管它很慢并且难以稳健地实现）。
我通常将曲线参数视为时间。这是伪代码：
- 将光标放在第一个控制点，找到周围的像素。
- 对于像素的每一侧（总共四个），通过求解二次方程检查贝塞尔曲线何时与其线相交。
- 在所有计算出的侧交叉口时间中，选择一个严格在未来发生的交叉路口时间，但要尽可能早。
- 根据最佳边移动到相邻像素。
- 将当前时间设置为该最佳侧路口的时间。
- 从第 2 步开始重复。
此算法一直有效，直到时间参数超过 1。另请注意，曲线恰好接触像素的一侧存在严重问题。我想它可以通过特殊检查解决。
这里是主要代码：
double WhenEquals(double p0, double p1, double p2, double val, double minp) {
//p0 * (1-t)^2 + p1 * 2t(1 - t) + p2 * t^2 = val
double qa = p0 + p2 - 2 * p1;
double qb = p1 - p0;
double qc = p0 - val;
assert(fabs(qa) > EPS); //singular case must be handled separately
double qd = qb * qb - qa * qc;
if (qd < -EPS)
return INF;
qd = sqrt(max(qd, 0.0));
double t1 = (-qb - qd) / qa;
double t2 = (-qb + qd) / qa;
if (t2 < t1) swap(t1, t2);
if (t1 > minp + EPS)
return t1;
else if (t2 > minp + EPS)
return t2;
return INF;
}
void DrawCurve(const Bezier &curve) {
int cell[2];
for (int c = 0; c < 2; c++)
cell[c] = int(floor(curve.pts[0].a[c]));
DrawPixel(cell[0], cell[1]);
double param = 0.0;
while (1) {
int bc = -1, bs = -1;
double bestTime = 1.0;
for (int c = 0; c < 2; c++)
for (int s = 0; s < 2; s++) {
double crit = WhenEquals(
curve.pts[0].a[c],
curve.pts[1].a[c],
curve.pts[2].a[c],
cell[c] + s, param
);
if (crit < bestTime) {
bestTime = crit;
bc = c, bs = s;
}
}
if (bc < 0)
break;
param = bestTime;
cell[bc] += (2*bs - 1);
DrawPixel(cell[0], cell[1]);
}
}
这里要注意的是 "line segments"，当创建得足够小时，相当于像素。贝塞尔曲线不是可线性遍历的曲线，因此我们不能像直线或圆弧那样轻松地一步 "skip ahead to the next pixel"。
当然，您可以对已有的 t 取任意点的切线，然后猜测下一个值 t' 会比哪个像素更远。但是，通常发生的情况是您猜测并猜错了，因为曲线的行为不是线性的，然后您检查 "off" 您的猜测如何，更正您的猜测，然后再次检查。重复直到收敛到下一个像素：这比将曲线展平为大量线段要慢得多，后者是一种快速操作。
如果您选择适合曲线长度的段数，考虑到它的渲染显示，没有人会告诉您曲线变平了。
有很多方法可以重新参数化贝塞尔曲线，但它们很昂贵，而且不同的典型曲线需要不同的重新参数化，所以这也不是真的更快。对于离散显示器最有用的往往是为您的曲线构建一个 LUT（查找 table），其长度适用于显示器上曲线的大小，然后使用该 LUT 作为您的基础绘图、交点检测等数据等等