在 Excel 中计算样条曲线 (y = f(x))

Question

我有一个 NURB 样条线，我手动拟合一些数据：

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

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

在 Excel 中，我需要一个

y = f(x)

函数来解决如下问题：

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

我的目标是向用户部署独立的 Excel 文件。 我知道如何使用像 Rhino3dm 这样的库来获取这些值。 我希望找到 Excel 原生的东西，比如可以做到这一点的 VBA、Python 或 Excel Labs 库。

我发现一些库可以 评估沿域的曲线，但不能评估 x 处的曲线。 我发现了可以评估 x 的库，但只能在 通过插值创建的曲线上。

Answer 1

在样条线上获取

y = f(x)

 的方法是对 

x = f(t)

 进行二分搜索（

De Boor 算法）。 我用 VBA 编写了它，以便在 Excel 中使用，并在这里分享它，以防其他人需要它。

以下是代码的快速概述：

CalcKnots() - 计算特定于 DeBoor 的结向量（以
• degree + 1

   而不是 

  degree

   倍数终止）。  它用于简化递归过程中的事情，使 CV 索引与结索引对齐。  我还缩放结向量以提高有限精度。

findKnotSpan() - 计算出 DeBoor 的正确结索引
• Tk
deBoor() - 递归地计算位置
• t

   处的样条线并返回 

  x

  （或 

  y

  ）的物理位置。

binSearchDeBoor() - 对“deBoor(t)”进行二分搜索，以在目标
• t 的

  tolerance

  内找到

  x

EvalSpline() - 这是您在 Excel 中用来获取
• y = f(x)

   的方法。  它调用 

  binSearchDeBoor(x)

   来查找 

  t

  ，然后调用 

  deBoor(t)

   来求解 

  y

还有 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