Bresenham's Algorithm Anti-aliased quadratic Bézier curve with line thickness

Question

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

Answer 1

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