Search code examples
pythonrecursionmathbezierbspline

Python ZeroDivisionError in Open Uniform B-Spline Curve

So I was working with implementing B-spline curves in python (and I am aware there are existing libraries, but I wanted to do it myself) and it worked well for non-open uniform B-spline curves as shown here:

enter image description here

The graph on the right showing the Basis Function (Cox de Boor's recursion formula) calculations correspond really well with what is shown in this video with the only difference being the number of control points. (https://www.youtube.com/watch?v=qhQrRCJ-mVg&t=2136):

enter image description here

As soon as I try to make this an open uniform B-Spline curve. I expect a graph like this (https://www.youtube.com/watch?v=qhQrRCJ-mVg&t=2501):

enter image description here

But instead, I get a ZeroDivisionError:

a = ((t - knotVector[i]) / (knotVector[i + j] - knotVector[i]) * CoxDeBoorRecursion(i, j - 1, t, knotVector))
ZeroDivisionError: float division by zero

where the knot vector is defined as [0, 0, 0, 1, 2, 3, 3, 3] like shown in the video at https://www.youtube.com/watch?v=qhQrRCJ-mVg&t=2460.

Here is my code:

def UniformBSpline(t, controlPoints, open=False):
    sumX, sumY = 0, 0
    degree = 2
    knotVector = [x for x in range(len(controlPoints) + 2 + degree)]

    if open:
        knotVector = knotVector[0:len(controlPoints)]
        for _ in range(degree):
            knotVector.insert(0, knotVector[0])
            knotVector.insert(-1, knotVector[-1])

    for i in range(len(controlPoints)):
        sumX += CoxDeBoorRecursion(i, degree, t, knotVector) * controlPoints[i].x
        sumY += CoxDeBoorRecursion(i, degree, t, knotVector) * controlPoints[i].y

    return sumX, sumY

def CoxDeBoorRecursion(i, j, t, knotVector):
    if j == 0:
        return 1 if knotVector[i] <= t < knotVector[i+1] else 0

    a = ((t - knotVector[i]) / (knotVector[i + j] - knotVector[i]) * CoxDeBoorRecursion(i, j - 1, t, knotVector))
    b = (((knotVector[i + j + 1] - t)/(knotVector[i + j + 1] - knotVector[i + 1])) * CoxDeBoorRecursion(i + 1, j - 1, t, knotVector))
    return a + b

This seems to be an issue in the function I am using which is defined as:

enter image description here

How can I fix this? and is this perhaps a better question for the Math StackExchange?

Solution

  • There are two errors in your code:

    The knot vector for uniform closed B-Splines is supposed to count up to the number of segments, not the number of control points. The number of segments is the number of control points minus the degree. I.e., if you have three points for a degree-2 spline (i.e., a simple Bezier curve), you get a single segment. For each additional point, you get one more segment. Hence, the knot vector has to be:

    if open:
    knotVector = [x for x in range(len(controlPoints) - degree + 1)]    
    for _ in range(degree):
        knotVector.insert(0, knotVector[0])
        knotVector.insert(-1, knotVector[-1])

    For four control points, this will give you

    [0, 0, 0, 1, 2, 2, 2]

    You have a similar problem for the closed case. That's why your spline has that corner at (0, 0). Haven't checked exactly, though.

    The second problem: The weights can have a divisor of zero if knots repeat. If this is the case, just ignore that part of the weight:

    d1 = (knotVector[i + j] - knotVector[i])
a = ((t - knotVector[i]) / d1 * CoxDeBoorRecursion(i, j - 1, t, knotVector)) if d1 > 0 else 0

d2 = (knotVector[i + j + 1] - knotVector[i + 1])
b = (((knotVector[i + j + 1] - t)/d2) * CoxDeBoorRecursion(i + 1, j - 1, t, knotVector)) if d2 > 0 else 0

    You can also extract that part into a separate function like the Wikipedia article does (function ω).