bezier path widening

Question

I have a bezier curve B with points S, C1, C2, E, and a positive number w representing width. Is there a way of quickly computing the control points of two bezier curves B1, B2 such that the stuff between B1 and B2 is the widened path represented by B?

More formally: compute the control points of good Bezier approximations to B1, B2, where B1 = {(x,y) + N(x,y)(w/2) | (x,y) in C}
B2 = {(x,y) - N(x,y)(w/2) | (x,y) in C},
where N(x,y) is the normal of C at (x,y).

I say good approximations because B1, B2 might not be polynomial curves (I'm not sure if they are).

Answer 1

The exact parallel of a bezier curve is quite ugly from a mathematical point of view (it requires 10th-degree polynomials).

What is easy to do is compute a widening from a polygonal approximation of the bezier (that is you compute line segments from the bezier and then move the points along the normals on the two sides of the curve).

This gives good results if your thickness isn't too big compared to the curvature... a "far parallel" instead is a monster on its own (and it's not even easy to find a definition of what is a parallel of an open curve that would make everyone happy).

Once you have two polylines for the two sides what you can do is finding a best approximating bezier for those paths if you need that representation. Once again I think that for "normal cases" (that is reasonably thin lines) even just a single bezier arc for each of the two sides should be quite accurate (the error should be much smaller than the thickness of the line).

EDIT: Indeed using a single bezier arc looks much worse than I would have expected even for reasonably normal cases. I tried also using two bezier arcs for each side and the result are better but still not perfect. The error is of course much smaller than the thickness of the line so unless lines are very thick it could be a reasonable option. In the following picture it's shown a thickened bezier (with per-point thickening), an approximation using a single bezier arc for each side and an approximation using two bezier arcs for each side.

EDIT 2: As requested I add the code I used to get the pictures; it's in python and requires only Qt. This code wasn't meant to be read by others so I used some tricks that probably I wouldn't use in real production code. The algorithm is also very inefficient but I didn't care about speed (this was meant to be a one-shot program to see if the idea works).

#
# This code has been written during an ego-pumping session on
# www.stackoverflow.com, while trying to reply to an interesting
# question. Do whatever you want with it but don't blame me if
# doesn't do what *you* think it should do or even if doesn't do
# what *I* say it should do.
#
# Comments of course are welcome...
#
# Andrea "6502" Griffini
#
# Requirements: Qt and PyQt
#
import sys
from PyQt4.Qt import *

QW = QWidget

bezlevels = 5

def avg(a, b):
    """Average of two (x, y) points"""
    xa, ya = a
    xb, yb = b
    return ((xa + xb)*0.5, (ya + yb)*0.5)

def bez3split(p0, p1, p2,p3):
    """
    Given the control points of a bezier cubic arc computes the
    control points of first and second half
    """
    p01 = avg(p0, p1)
    p12 = avg(p1, p2)
    p23 = avg(p2, p3)
    p012 = avg(p01, p12)
    p123 = avg(p12, p23)
    p0123 = avg(p012, p123)
    return [(p0, p01, p012, p0123),
            (p0123, p123, p23, p3)]

def bez3(p0, p1, p2, p3, levels=bezlevels):
    """
    Builds a bezier cubic arc approximation using a fixed
    number of half subdivisions.
    """
    if levels <= 0:
        return [p0, p3]
    else:
        (a0, a1, a2, a3), (b0, b1, b2, b3) = bez3split(p0, p1, p2, p3)
        return (bez3(a0, a1, a2, a3, levels-1) +
                bez3(b0, b1, b2, b3, levels-1)[1:])

def thickPath(pts, d):
    """
    Given a polyline and a distance computes an approximation
    of the two one-sided offset curves and returns it as two
    polylines with the same number of vertices as input.

    NOTE: Quick and dirty approach, just uses a "normal" for every
          vertex computed as the perpendicular to the segment joining
          the previous and next vertex.
          No checks for self-intersections (those happens when the
          distance is too big for the local curvature), and no check
          for degenerate input (e.g. multiple points).
    """
    l1 = []
    l2 = []
    for i in xrange(len(pts)):
        i0 = max(0, i - 1)             # previous index
        i1 = min(len(pts) - 1, i + 1)  # next index
        x, y = pts[i]
        x0, y0 = pts[i0]
        x1, y1 = pts[i1]
        dx = x1 - x0
        dy = y1 - y0
        L = (dx**2 + dy**2) ** 0.5
        nx = - d*dy / L
        ny = d*dx / L
        l1.append((x - nx, y - ny))
        l2.append((x + nx, y + ny))
    return l1, l2

def dist2(x0, y0, x1, y1):
    "Squared distance between two points"
    return (x1 - x0)**2 + (y1 - y0)**2

def dist(x0, y0, x1, y1):
    "Distance between two points"
    return ((x1 - x0)**2 + (y1 - y0)**2) ** 0.5

def ibez(pts, levels=bezlevels):
    """
    Inverse-bezier computation.
    Given a list of points computes the control points of a
    cubic bezier arc that approximates them.
    """
    #
    # NOTE:
    #
    # This is a very specific routine that only works
    # if the input has been obtained from the computation
    # of a bezier arc with "levels" levels of subdivisions
    # because computes the distance as the maximum of the
    # distances of *corresponding points*.
    # Note that for "big" changes in the input from the
    # original bezier I dont't think is even true that the
    # best parameters for a curve-curve match would also
    # minimize the maximum distance between corresponding
    # points. For a more general input a more general
    # path-path error estimation is needed.
    #
    # The minimizing algorithm is a step descent on the two
    # middle control points starting with a step of about
    # 1/10 of the lenght of the input to about 1/1000.
    # It's slow and ugly but required no dependencies and
    # is just a bunch of lines of code, so I used that.
    #
    # Note that there is a closed form solution for finding
    # the best bezier approximation given starting and
    # ending points and a list of intermediate parameter
    # values and points, and this formula also could be
    # used to implement a much faster and accurate
    # inverse-bezier in the general case.
    # If you care about the problem of inverse-bezier then
    # I'm pretty sure there are way smarter methods around.
    #
    # The minimization used here is very specific, slow
    # and not so accurate. It's not production-quality code.
    # You have been warned.
    #

    # Start with a straight line bezier arc (surely not
    # the best choice but this is just a toy).
    x0, y0 = pts[0]
    x3, y3 = pts[-1]
    x1, y1 = (x0*3 + x3) / 4.0, (y0*3 + y3) / 4.0
    x2, y2 = (x0 + x3*3) / 4.0, (y0 + y3*3) / 4.0
    L = sum(dist(*(pts[i] + pts[i-1])) for i in xrange(len(pts) - 1))
    step = L / 10
    limit = step / 100

    # Function to minimize = max((a[i] - b[i])**2)
    def err(x0, y0, x1, y1, x2, y2, x3, y3):
        return max(dist2(*(x+p)) for x, p in zip(pts, bez3((x0, y0), (x1, y1),
                                                           (x2, y2), (x3, y3),
                                                           levels)))
    while step > limit:
        best = None
        for dx1 in (-step, 0,  step):
            for dy1 in (-step, 0, step):
                for dx2 in (-step, 0, step):
                    for dy2 in (-step, 0, step):
                        e = err(x0, y0,
                                x1+dx1, y1+dy1,
                                x2+dx2, y2+dy2,
                                x3, y3)
                        if best is None or e < best[0] * 0.9999:
                            best = e, dx1, dy1, dx2, dy2
        e, dx1, dy1, dx2, dy2 = best
        if (dx1, dy1, dx2, dy2) == (0, 0, 0, 0):
            # We got to a minimum for this step => refine
            step *= 0.5
        else:
            # We're still moving
            x1 += dx1
            y1 += dy1
            x2 += dx2
            y2 += dy2

    return [(x0, y0), (x1, y1), (x2, y2), (x3, y3)]

def poly(pts):
    "Converts a list of (x, y) points to a QPolygonF)"
    return QPolygonF(map(lambda p: QPointF(*p), pts))

class Viewer(QW):
    def __init__(self, parent):
        QW.__init__(self, parent)
        self.pts = [(100, 100), (200, 100), (200, 200), (100, 200)]
        self.tracking = None    # Mouse dragging callback
        self.ibez = 0           # Thickening algorithm selector

    def sizeHint(self):
        return QSize(900, 700)

    def wheelEvent(self, e):
        # Moving the wheel changes between
        # - original polygonal thickening
        # - single-arc thickening
        # - double-arc thickening
        self.ibez = (self.ibez + 1) % 3
        self.update()

    def paintEvent(self, e):
        dc = QPainter(self)
        dc.setRenderHints(QPainter.Antialiasing)

        # First build the curve and the polygonal thickening
        pts = bez3(*self.pts)
        l1, l2 = thickPath(pts, 15)

        # Apply inverse bezier computation if requested
        if self.ibez == 1:
            # Single arc
            l1 = bez3(*ibez(l1))
            l2 = bez3(*ibez(l2))
        elif self.ibez == 2:
            # Double arc
            l1 = (bez3(*ibez(l1[:len(l1)/2+1], bezlevels-1)) +
                  bez3(*ibez(l1[len(l1)/2:], bezlevels-1))[1:])
            l2 = (bez3(*ibez(l2[:len(l2)/2+1], bezlevels-1)) +
                  bez3(*ibez(l2[len(l2)/2:], bezlevels-1))[1:])

        # Draw results
        dc.setBrush(QBrush(QColor(0, 255, 0)))
        dc.drawPolygon(poly(l1 + l2[::-1]))
        dc.drawPolyline(poly(pts))
        dc.drawPolyline(poly(self.pts))

        # Draw control points
        dc.setBrush(QBrush(QColor(255, 0, 0)))
        dc.setPen(QPen(Qt.NoPen))
        for x, y in self.pts:
            dc.drawEllipse(QRectF(x-3, y-3, 6, 6))

        # Display the algorithm that has been used
        dc.setPen(QPen(QColor(0, 0, 0)))
        dc.drawText(20, 20,
                    ["Polygonal", "Single-arc", "Double-arc"][self.ibez])

    def mousePressEvent(self, e):
        # Find closest control point
        i = min(range(len(self.pts)),
                key=lambda i: (e.x() - self.pts[i][0])**2 +
                              (e.y() - self.pts[i][1])**2)

        # Setup a callback for mouse dragging
        self.tracking = lambda p: self.pts.__setitem__(i, p)

    def mouseMoveEvent(self, e):
        if self.tracking:
            self.tracking((e.x(), e.y()))
            self.update()

    def mouseReleaseEvent(self, e):
        self.tracking = None

# Qt boilerplate
class MyDialog(QDialog):
    def __init__(self, parent):
        QDialog.__init__(self, parent)
        self.ws = Viewer(self)
        L = QVBoxLayout(self)
        L.addWidget(self.ws)
        self.setModal(True)
        self.show()

app = QApplication([])
aa = MyDialog(None)
aa.exec_()
aa = None