Evaluate Spline (y = f(x)) in Excel

Question

I have a NURBs spline I manually fit over some data:

Control Points (x, y):
(0, 5)
(0, 3)
(1, 2)
(10, 2)
(100, 5)

Knot Vector [0, 0, 0, 1, 2, 2, 2]

In Excel, I need a y = f(x) function that would solve like so:

f(0) = 5
f(0.4) = 3.148
f(1) = 2.679
f(10) = 2.155

My goal is to deploy a standalone Excel file to users. I know how to get these values with a library like Rhino3dm. I'm hoping to find something native to Excel, like a VBA, Python, or Excel Labs library that can do this.

I've found libraries that can evaluate a curve along the domain, but not at x. And I've found libraries that can evaluate x, but only on curves created by interpolation.

Answer 1

The way to get y = f(x) on a spline is with a binary search over x = f(t) (De Boor's algorithm). I wrote it in VBA for use in Excel and am sharing it here in case anyone else needs it.

Here's a quick overview of the code:

• CalcKnots() - calculates a knot vector specific to DeBoor's (terminated with degree + 1 instead of degree multiples). It's used to simplify things during recursion, causing the CV index to align with the knot index. I also scale the knot vector to increase finite precision.
• findKnotSpan() - figures out the proper knot index, Tk, for DeBoor's
• deBoor() - recursively evaluates a spline at position t and returns the physical location of x (or y).
• binSearchDeBoor() - binary search over deBoor(t) to find a t within tolerance of your target x
• EvalSpline() - This is what you use in Excel to get y = f(x). It calls binSearchDeBoor(x) to find t, then calls deBoor(t) to solve for y

And the VBA:

Private Function CalcKnots(Degree As Integer, cvX() As Double, cvY() As Double) As Double()
    
    Dim cvLen As Integer, spans As Integer, knotLen As Integer, i As Integer
    Dim xyDistance As Double, knotDelta As Double, acc As Double
    
    cvLen = UBound(cvX) + 1
    spans = cvLen - Degree
    
    ' DeBoors requires degree + 1 knot multiples at the ends
    knotLen = Degree + cvLen + 1
    Dim knots() As Double
    ReDim knots(knotLen - 1)

    ' scale the knot vector with the spline length to minimize issues with precision
    ' since it's a rough scaling, use the x distance + y distance, no need to add a sqrt()
    xyDistance = Abs((cvX(cvLen - 1) - cvX(0))) + Abs((cvY(cvLen - 1) - cvY(0)))
    knotDelta = xyDistance / (spans)
    acc = knotDelta
    
    For i = 0 To Degree
        knots(i) = 0
    Next i
    
    For i = Degree + 1 To Degree + spans - 1
        knots(i) = acc
        acc = acc + knotDelta
    Next i
    
    For i = Degree + spans To knotLen - 1
        knots(i) = acc
    Next i
    
    CalcKnots = knots
    
End Function

Private Function findKnotSpan(t As Double, knots() As Double) As Integer
    ' find knot span index Tk, where Tk < t < Tk+1
    Dim k As Integer
    For k = 0 To UBound(knots)
        If t < knots(k) Then
            findKnotSpan = k - 1
            Exit Function
        End If
    Next k
End Function

Public Function deBoor(idx As Integer, Degree As Integer, Tk As Integer, t As Double, knots() As Double, cvX() As Double) As Double
    ' idx = recursion index, initially the degree
    ' Tk = knot span index, where Tk < t < Tk+1
    ' t = domain parameter
    ' knots = array of knot values
    ' cvX = array of 1d control points
    
    Dim alpha As Double, left As Double, right As Double
    If idx = 0 Then
        deBoor = cvX(Tk)
    Else
        alpha = (t - knots(Tk)) / (knots(Tk + Degree + 1 - idx) - knots(Tk))
        left = deBoor(idx - 1, Degree, Tk - 1, t, knots, cvX) * (1 - alpha)
        right = deBoor(idx - 1, Degree, Tk, t, knots, cvX) * alpha
        deBoor = left + right
    End If
End Function

Private Function binSearchDeBoor(Degree As Integer, TargetValue As Double, Tolerance As Double, knots() As Double, cvX() As Double) As Double
    Dim knotsLen As Integer, Tk As Integer
    knotsLen = UBound(knots) + 1
    
    Dim delta As Double, dMin As Double, dMax As Double, t As Double, tMin As Double, tMax As Double
    dMin = 0 - Tolerance
    dMax = 0 + Tolerance
    tMin = 0
    tMax = knots(knotsLen - 1)
    t = tMax / 2
    
    Do
        Tk = findKnotSpan(t, knots)
        delta = deBoor(Degree, Degree, Tk, t, knots, cvX) - TargetValue
        Select Case delta
            Case Is < dMin
                tMin = t
                t = tMin + ((tMax - tMin) / 2)
            Case Is > dMax
                tMax = t
                t = tMin + ((tMax - tMin) / 2)
            Case Else
                binSearchDeBoor = t
                Exit Function
        End Select
    Loop
End Function

Public Function EvalSpline(Degree as Integer, TargetValue as Double, SearchRange As Range, ResultRange As Range) As Double
    If SearchRange.Count <> ResultRange.Count Then
        Err.Raise vbObjectError + 1000, "EvalSpline", "Both ranges must be the same length"
    End If
    Dim rangeLen As Integer
    rangeLen = SearchRange.Count
    Dim cvS() As Double
    Dim cvR() As Double
    ReDim cvS(rangeLen - 1)
    ReDim cvR(rangeLen - 1)
    For i = 0 To rangeLen - 1
        cvS(i) = SearchRange(i + 1)
        cvR(i) = ResultRange(i + 1)
    Next i
    
    Dim knots() As Double
    knots = CalcKnots(Degree, cvS, cvR)
    
    Dim t As Double, Tk As Integer
    t = binSearchDeBoor(Degree, TargetValue, Tolerance, knots, cvS)
    Tk = findKnotSpan(t, knots)
    EvalSpline = deBoor(Degree, Degree, Tk, t, knots, cvR)
    
End Function