最常用的三次样条插值函数

'*******************************************************************
' 三次样条插值函数
'*******************************************************************
Sub tspLine(ByVal n As Integer, ByVal ch As Integer, ByVal tx1 As Single, ByVal tx2 As Single, ByVal ty1 As Single, ByVal ty2 As Single)
Dim a(10) As Single, b(10) As Single, c(10) As Single, dx(10) As Single, dy(10) As Single
Dim qx(10) As Single, qy(10) As Single
Dim tt As Single, bx3 As Single, bx4 As Single, by3 As Single, by4 As Single
Dim cx As Single, cy As Single, t(10) As Single, px(10) As Single, py(10) As Single
Dim u(3) As Single, v(3) As Single, I As Integer
num = 0
For I = 1 To n
t(I) = hypot(X(I) - X(I - 1), Y(I) - Y(I - 1))
Next I
Select Case ch
Case 0 '抛物条件
  u(0) = (X(1) - X(0)) / t(1): u(1) = (X(2) - X(1)) / t(2)
  u(2) = (u(1) - u(0)) / (t(2) + t(1))
  tx1 = u(0) - u(2) * t(1)
  u(0) = (Y(1) - Y(0)) / t(1): u(1) = (Y(2) - Y(1)) / t(2)
  u(2) = (u(1) - u(0)) / (t(2) + t(1))
  ty1 = u(0) - u(2) * t(1)
  u(0) = (X(n) - X(n - 1)) / t(n): u(1) = (X(n - 1) - X(n - 2)) / t(n - 1)
  u(2) = (u(0) - u(1)) / (t(n) + t(n - 1))
  tx2 = u(0) + u(2) * t(n)
  u(0) = (Y(n) - Y(n - 1)) / t(n): u(1) = (Y(n - 1) - Y(n - 2)) / t(n - 1)
  u(2) = (u(0) - u(1)) / (t(n) + t(n - 1))
  ty2 = u(0) + u(2) * t(n)
Case 1 '夹持条件
a(0) = 1: c(0) = 0: dx(0) = tx1: dy(0) = ty1
a(n) = 1: b(n) = 0: dx(n) = tx2: dy(n) = ty2
Case 2 '自由条件
a(0) = 2: c(0) = 1
dx(0) = 3 * (X(1) - X(0)) / t(1): dy(0) = 3 * (Y(1) - Y(0)) / t(1)
a(n) = 2: b(n) = 1
dx(n) = 3 * (X(n) - X(n - 1)) / t(n): dy(n) = 3 * (Y(n) - Y(n - 1)) / t(n)
Case 3 '循环条件
a(0) = 2: c(0) = 1
dx(0) = 3 * (X(1) - X(0)) / t(1) - (t(1) * (X(2) - X(1)) / t(2) - X(1) + X(0)) / (t(1) + t(2))
dy(0) = 3 * (Y(1) - Y(0)) / t(1) - (t(1) * (Y(2) - Y(1)) / t(2) - Y(1) + Y(0)) / (t(1) + t(2))
a(n) = 2: b(n) = 1
dx(n) = 3 * (X(n) - X(n - 1)) / t(n)
dx(n) = dx(n) + (X(n) - X(n - 1) - t(n) * (X(n - 1) - X(n - 2)) / t(n - 1)) / (t(n) + t(n - 1))
dy(n) = 3 * (Y(n) - Y(n - 1)) / t(n)
dy(n) = dy(n) + (Y(n) - Y(n - 1) - t(n) * (Y(n - 1) - Y(n - 2)) / t(n - 1)) / (t(n) + t(n - 1))
End Select

'计算方程组系数阵和常数阵
For I = 1 To n - 1
a(I) = 2 * (t(I) + t(I + 1)): b(I) = t(I + 1): c(I) = t(I)
dx(I) = 3 * (t(I) * (X(I + 1) - X(I)) / t(I + 1) + t(I + 1) * (X(I) - X(I - 1)) / t(I))
dy(I) = 3 * (t(I) * (Y(I + 1) - Y(I)) / t(I + 1) + t(I + 1) * (Y(I) - Y(I - 1)) / t(I))
Next I
'采用追赶法解方程组
c(0) = c(0) / a(0)
For I = 1 To n - 1
a(I) = a(I) - b(I) * c(I - 1): c(I) = c(I) / a(I)
Next I
a(n) = a(n) - b(n) * c(I - 1)
qx(0) = dx(0) / a(0): qy(0) = dy(0) / a(0)
For I = 1 To n
qx(I) = (dx(I) - b(I) * qx(I - 1)) / a(I)
qy(I) = (dy(I) - b(I) * qy(I - 1)) / a(I)
Next I
px(n) = qx(n): py(n) = qy(n)
For I = n - 1 To 0 Step -1
px(I) = qx(I) - c(I) * px(I + 1)
py(I) = qy(I) - c(I) * py(I + 1)
Next I
'计算曲线上点的坐标
For I = 0 To n - 1
bx3 = (3 * (X(I + 1) - X(I)) / t(I + 1) - 2 * px(I) - px(I + 1)) / t(I + 1)
bx4 = ((2 * (X(I) - X(I + 1)) / t(I + 1) + px(I) + px(I + 1)) / t(I + 1)) / t(I + 1)
by3 = (3 * (Y(I + 1) - Y(I)) / t(I + 1) - 2 * py(I) - py(I + 1)) / t(I + 1)
by4 = ((2 * (Y(I) - Y(I + 1)) / t(I + 1) + py(I) + py(I + 1)) / t(I + 1)) / t(I + 1)
tt = 0
While (tt <= t(I + 1))
cx = X(I) + (px(I) + (bx3 + bx4 * tt) * tt) * tt
cy = Y(I) + (py(I) + (by3 + by4 * tt) * tt) * tt
u1(num) = cx: v1(num) = cy: num = num + 1: tt = tt + 0.5
Wend
If I = 0 Then num1 = num
If I = n - 1 Then num2 = num
u1(num) = X(I + 1): v1(num) = Y(I + 1): num = num + 1
Next I
End Sub