贝塞尔路径加宽

Question

我有一条贝塞尔曲线 B，其中有点 S、C1、C2、E，以及代表宽度的正数 w。有没有一种方法可以快速计算两条贝塞尔曲线B1、B2的控制点，使得B1和B2之间的东西就是B表示的加宽路径？

更正式地说：计算 B1、B2 的良好 Bezier 近似的控制点，其中 B1 = {(x,y) + N(x,y)(w/2) | C 中的 (x,y)}
B2 = {(x,y) - N(x,y)(w/2) | B2 = {(x,y) - N(x,y)(w/2) | C 中的 (x,y)},
其中 N(x,y) 是法线 C 在 (x,y) 处。

我说好的近似值是因为 B1、B2 可能不是多项式曲线（我不确定它们是否是）。

Answer 1

从数学角度来看，贝塞尔曲线的精确平行相当难看（它需要 10 次多项式）。

容易做的是从贝塞尔曲线的多边形近似计算加宽（即从贝塞尔曲线计算线段，然后沿着曲线两侧的法线移动点）。

如果你的厚度与曲率相比不是太大，这会给出很好的结果......“远平行”本身就是一个怪物（并且甚至很难找到什么是开放曲线的平行的定义这会让每个人都高兴）。

一旦两侧有了两条多段线，如果您需要该表示，您可以做的是为这些路径找到最佳近似贝塞尔曲线。我再次认为，对于“正常情况”（即相当细的线），即使两侧各有一条贝塞尔曲线也应该非常准确（误差应该比线的粗细小得多）。

编辑：实际上，即使对于相当正常的情况，使用单个贝塞尔曲线看起来也比我预期的要糟糕得多。我还尝试在每一侧使用两个贝塞尔弧，结果更好，但仍然不完美。当然，误差远小于线条的粗细，因此除非线条非常粗，否则它可能是一个合理的选择。在下图中，它显示了加厚的贝塞尔曲线（每点加厚）、每侧使用单个贝塞尔曲线的近似值以及每侧使用两个贝塞尔曲线的近似值。

编辑2：根据要求，我添加了用于获取图片的代码;它是用 python 编写的，只需要 Qt。这段代码并不适合其他人阅读，因此我使用了一些可能不会在实际生产代码中使用的技巧。该算法也非常低效，但我并不关心速度（这本来是一个一次性程序，看看这个想法是否有效）。

#
# 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