Bezier Pixel Curve

Question

The quadratic / cubic Bezier curve code that I find through google basically works by dividing the line into a series of points and connecting them with straight lines. Rasterization occurs in a linear algorithm, and not in one of them. Algorithms such as Bresenham work pixel by pixel to rasterize the line and can be optimized (see Po-Han Lin solution ).

What is a quadratic bezier curve algorithm that works on pixel pixel algorithms, and not by constructing a series of points?

Answer 1

A variant of the Bresenem algorithm works with quadratic functions, such as circles, ellipses, and parabolas, so it should work with quadratic Bezier curves.

I tried to implement the implementation, but then found one on the Internet: http://members.chello.at/~easyfilter/bresenham.html .

If you need more detailed or additional examples, the page mentioned above has a link to a 100-page PDF document on the method: http://members.chello.at/~easyfilter/Bresenham.pdf .

Here is the code of the Alois Zinl site for constructing any quadratic Bezier curve. The first procedure divides the curve into changes in horizontal and vertical gradients:

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 */ } 

The second procedure actually displays a segment of the Bezier curve (one without gradient changes):

 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 */ } 

Code for smoothing is also available on the site.

Relevant Zingl site features for cubic Bezier curves

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

and

 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 */ } 

As Mike Pomax Kamermans noted, the solution for cubic Bezier curves on the site is not complete; in particular, there are problems with smoothing Bezier cubic curves, and the discussion of rational Bezier cubic curves is incomplete.

Answer 2

You can use the De Castellau algorithm to divide the curve into enough parts, each of which is a pixel.

This is the equation for finding the point [x, y] on a quadratic curve in the interval T:

 // 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 starting control point T==0.00. // At the curve 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; 

For practical use of this equation, you can enter about 1000 T values ​​from 0.00 to 1.00. This results in a set of 1000 points guaranteed along a quadratic curve.

Computing 1000 points along the curve is probably over-sampling (some calculated points will have the same pixel coordinate), so you will want to deduplicate 1000 points until the set displays unique pixel coordinates along the curve.

There is a similar equation for Kubica Bezier curves.

Here is an example of code that displays a quadratic curve as a set of calculated pixels:

 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> 

Answer 3

The fact is that "line segments" when they are created small enough are equivalent to pixels. Bezier curves are not linearly intersecting curves, so we cannot easily "go to the next pixel" in one step, as we can for lines or circular arcs.

You could, of course, take the tangent at any time for t that you already have, and then guess which next value of t' will lie further in the pixel. However, it usually happens that you guess and think it is wrong, because the curve does not behave linearly, then you check how your assumption is “turned off”, correct your guess, and then check again. Repeat until you get to the next pixel: this is much slower than just smoothing the curve to a large number of line segments, which is a quick operation.

If you select the number of segments so that they correspond to the length of the curve, then, given the displayed display, no one can say that you have smoothed the curve.

There are methods for reparametrizing Bezier curves, but they are expensive, and different canonical curves require different reparametrization, so it really is not faster. Most preferable for discrete displays is building a LUT (lookup table) for your curve with a length that works for the size of the curve on the display, and then using that LUT as the base data for drawing, detecting intersections, etc. etc.

Answer 4

First of all, I would like to say that the fastest and most reliable way to visualize Bezier curves is to approximate them with a polyline through an adaptive unit, and then display the polyline. @MarkE's approach of drawing many points selected on a curve is pretty quick, but it can skip pixels. Here I describe another approach that is closest to line spacing (although it is slow and hard to implement reliably).

Usually I will consider the curve parameter as time. Here is the pseudo code:

• Place the cursor on the first control point, find the surrounding pixel.
• For each side of the pixel (four in total), check when the Bezier curves cross their line, allowing quadratic equations.
• Among all calculated lateral intersection points, choose the one that will be executed strictly in the future, but as soon as possible.
• Move to an adjacent pixel depending on which side was better.
• Set the current time to the time of best lateral intersection.
• Repeat step 2.

This algorithm works until the time parameter exceeds one. Also note that he has serious problems with curves that touch the side of the pixel. I believe this is solvable with a special check.

Here is the main code:

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

Full code is available here . It uses loadbmp.h, here .

Answer 5

There is a Bresenham Bezier curve algorithm that works for some Bezier curves. In particular, those curves where the control point is inside the bounding box of the start and end points. That is, we can guarantee the absence of gradient shifts, that if we move in the + x or + y direction, we will ALWAYS move in this direction.

http://members.chello.at/~easyfilter/bresenham.html

However, its licensing code is completely ambiguous.

But his method is to calculate the curvature and perform integer math to find the point at which the curvature will have a difference in pixels, then subtract the bit of curvature again and continue. It is this Bresenham method that is applied to the curve algorithm.

Again, I cannot confirm the licensing of the source code for the following code.

 void plotQuadBezierSeg(int x0, int y0, int x1, int y1, int x2, int y2) { int dx_p1_p2 = x2 - x1; int dy_p1_p2 = y2 - y1; long dx_p0_p1 = x0 - x1; long dy_p0_p1 = y0 - y1; long xy; /* relative values for checks */ double dx, dy, err; double curvature = dx_p0_p1 * dy_p1_p2 - dy_p0_p1 * dx_p1_p2; /* curvature */ assert (dx_p0_p1 * dx_p1_p2 <= 0 && dy_p0_p1 * dy_p1_p2 <= 0); /* sign of gradient must not change */ if (dx_p1_p2 * (long) dx_p1_p2 + dy_p1_p2 * (long) dy_p1_p2 > dx_p0_p1 * dx_p0_p1 + dy_p0_p1 * dy_p0_p1) { /* begin with longer part */ x2 = x0; x0 = dx_p1_p2 + x1; y2 = y0; y0 = dy_p1_p2 + y1; curvature = -curvature; /* swap P0 P2 */ } if (curvature != 0) { /* no straight line */ dx_p0_p1 += dx_p1_p2; dx_p0_p1 *= dx_p1_p2 = x0 < x2 ? 1 : -1; /* x step direction */ dy_p0_p1 += dy_p1_p2; dy_p0_p1 *= dy_p1_p2 = y0 < y2 ? 1 : -1; /* y step direction */ xy = 2 * dx_p0_p1 * dy_p0_p1; dx_p0_p1 *= dx_p0_p1; dy_p0_p1 *= dy_p0_p1; /* differences 2nd degree */ if (curvature * dx_p1_p2 * dy_p1_p2 < 0) { /* negated curvature? */ dx_p0_p1 = -dx_p0_p1; dy_p0_p1 = -dy_p0_p1; xy = -xy; curvature = -curvature; } dx = 4.0 * dy_p1_p2 * curvature * (x1 - x0) + dx_p0_p1 - xy; /* differences 1st degree */ dy = 4.0 * dx_p1_p2 * curvature * (y0 - y1) + dy_p0_p1 - xy; dx_p0_p1 += dx_p0_p1; dy_p0_p1 += dy_p0_p1; 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) ? 1 : 0; /* save value for test of y step */ if (2 * err > dy) { x0 += dx_p1_p2; dx -= xy; err += dy += dy_p0_p1; } /* x step */ if (y1 == 1) { y0 += dy_p1_p2; dy -= xy; err += dx += dx_p0_p1; } /* y step */ } while (dy < dx); /* gradient negates -> algorithm fails */ } drawLine(x0, y0, x2, y2); /* plot remaining part to end */ } 

drawLine there refers to the Bresenham line drawing algorithm and is applied there only to the end.

It also has a fantastic demo canvas in javascript:

http://members.chello.at/~easyfilter/canvas.html