form1.frm

三次样条插值函数

VERSION 5.00
Begin VB.Form Form1 
   Caption         =   "Form1"
   ClientHeight    =   8865
   ClientLeft      =   60
   ClientTop       =   450
   ClientWidth     =   11535
   LinkTopic       =   "Form1"
   ScaleHeight     =   8865
   ScaleWidth      =   11535
   StartUpPosition =   3  '窗口缺省
   Begin VB.PictureBox Picture1 
      Height          =   5415
      Left            =   1680
      ScaleHeight     =   5355
      ScaleWidth      =   7155
      TabIndex        =   0
      Top             =   1320
      Width           =   7215
   End
   Begin VB.Label Label1 
      Caption         =   "Label1"
      Height          =   495
      Left            =   2400
      TabIndex        =   1
      Top             =   480
      Width           =   1935
   End
   Begin VB.Line Line1 
      X1              =   1080
      X2              =   1080
      Y1              =   2400
      Y2              =   4320
   End
End
Attribute VB_Name = "Form1"
Attribute VB_GlobalNameSpace = False
Attribute VB_Creatable = False
Attribute VB_PredeclaredId = True
Attribute VB_Exposed = False
Option Explicit

Dim Xi() As Single, Yi() As Single, u1(40000) As Single, v1(40000) As Single
Dim a(100) As Single, b(100) As Single, c(100) As Single, dx(100) As Single, dy(100) As Single
Dim qx(100) As Single, qy(100) As Single
Dim Num As Long, Nn As Integer
Dim ii As Long
Const PP = 100
Function hypot(ByVal X As Single, ByVal Y As Single)
    hypot = Sqr(X * X + Y * Y)
End Function

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)
   '三次样条方程系数计算
   '参数意义:
   'n: 给定点的个数,ch:边界条件类型,(tx1,ty1)  (tx2,ty2) :还不知道什么意思,是边界条件吧,呵呵
   
   On Error Resume Next
    Dim tt As Single, bx3 As Single, bx4 As Single, by3 As Single, by4 As Single
    Dim cx As Single, cy As Single, t(100) As Single, px(100) As Single, py(100) As Single
    Dim u(3) As Single, v(3) As Single, i As Integer
    
    Num = 0
    For i = 1 To n
        t(i) = hypot(Xi(i) - Xi(i - 1), Yi(i) - Yi(i - 1))
    Next i
    Select Case ch
    Case 0 '抛物条件
        u(0) = (Xi(1) - Xi(0)) / t(1): u(1) = (Xi(2) - Xi(1)) / t(2)
        u(2) = (u(1) - u(0)) / (t(2) + t(1))
        tx1 = u(0) - u(2) * t(1)
        u(0) = (Yi(1) - Yi(0)) / t(1): u(1) = (Yi(2) - Yi(1)) / t(2)
        u(2) = (u(1) - u(0)) / (t(2) + t(1))
        ty1 = u(0) - u(2) * t(1)
        u(0) = (Xi(n) - Xi(n - 1)) / t(n): u(1) = (Xi(n - 1) - Xi(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) = (Yi(n) - Yi(n - 1)) / t(n): u(1) = (Yi(n - 1) - Yi(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 * (Xi(1) - Xi(0)) / t(1): dy(0) = 3 * (Yi(1) - Yi(0)) / t(1)
        a(n) = 2: b(n) = 1
        dx(n) = 3 * (Xi(n) - Xi(n - 1)) / t(n): dy(n) = 3 * (Yi(n) - Yi(n - 1)) / t(n)
    Case 3 '循环条件  周期样条
        a(0) = 2: c(0) = 1
        dx(0) = 3 * (Xi(1) - Xi(0)) / t(1) - (t(1) * (Xi(2) - Xi(1)) / t(2) - Xi(1) + Xi(0)) / (t(1) + t(2))
        dy(0) = 3 * (Yi(1) - Yi(0)) / t(1) - (t(1) * (Yi(2) - Yi(1)) / t(2) - Yi(1) + Yi(0)) / (t(1) + t(2))
        a(n) = 2: b(n) = 1
        dx(n) = 3 * (Xi(n) - Xi(n - 1)) / t(n)
        dx(n) = dx(n) + (Xi(n) - Xi(n - 1) - t(n) * (Xi(n - 1) - Xi(n - 2)) / t(n - 1)) / (t(n) + t(n - 1))
        dy(n) = 3 * (Yi(n) - Yi(n - 1)) / t(n)
        dy(n) = dy(n) + (Yi(n) - Yi(n - 1) - t(n) * (Yi(n - 1) - Yi(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) * (Xi(i + 1) - Xi(i)) / t(i + 1) + t(i + 1) * (Xi(i) - Xi(i - 1)) / t(i))
        dy(i) = 3 * (t(i) * (Yi(i + 1) - Yi(i)) / t(i + 1) + t(i + 1) * (Yi(i) - Yi(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 * (Xi(i + 1) - Xi(i)) / t(i + 1) - 2 * px(i) - px(i + 1)) / t(i + 1)
        bx4 = ((2 * (Xi(i) - Xi(i + 1)) / t(i + 1) + px(i) + px(i + 1)) / t(i + 1)) / t(i + 1)
        by3 = (3 * (Yi(i + 1) - Yi(i)) / t(i + 1) - 2 * py(i) - py(i + 1)) / t(i + 1)
        by4 = ((2 * (Yi(i) - Yi(i + 1)) / t(i + 1) + py(i) + py(i + 1)) / t(i + 1)) / t(i + 1)
        tt = 0
        While (tt <= t(i + 1))
         cx = Xi(i) + (px(i) + (bx3 + bx4 * tt) * tt) * tt
         cy = Yi(i) + (py(i) + (by3 + by4 * tt) * tt) * tt
         u1(Num) = cx: v1(Num) = cy: Num = Num + 1: tt = tt + 0.5
        Wend
        u1(Num) = Xi(i + 1): v1(Num) = Yi(i + 1): Num = Num + 1
    Next i
End Sub

Private Sub 画样条()
On Error Resume Next

    Dim i As Long
    
    Picture1.Cls

    Picture1.ForeColor = RGB(200, 200, 200)
    Picture1.Line (Picture1.ScaleLeft, 0)-(Picture1.ScaleLeft + Picture1.ScaleWidth, 0)
    Picture1.Line (0, Picture1.ScaleTop)-(0, Picture1.ScaleTop + Picture1.ScaleHeight)
    Picture1.ForeColor = RGB(255, 0, 0)

    Picture1.DrawWidth = 4
    For i = 0 To Nn - 1
        Picture1.PSet (Xi(i), Yi(i)), vbBlue    '画点
    Next i
    tspLine Nn - 1, 2, 0, 0, 0, 0 '求点坐标,
    Picture1.DrawWidth = 1
    Picture1.PSet (u1(0), v1(0))
    For i = 1 To Num - 1
        Picture1.Line -(u1(i), v1(i))   '画线
    Next i
End Sub
Private Sub Form_Load()
    Picture1.Scale (-500, 500)-(500, -500)
    Picture1.ForeColor = RGB(200, 200, 200)
    Picture1.Line (Picture1.ScaleLeft, 0)-(Picture1.ScaleLeft + Picture1.ScaleWidth, 0)
    Picture1.Line (0, Picture1.ScaleTop)-(0, Picture1.ScaleTop + Picture1.ScaleHeight)
    
End Sub

Private Sub Picture1_MouseDown(Button As Integer, Shift As Integer, X As Single, Y As Single)
    If Button = 1 Then
'        ReDim Preserve Xi(Nn)
'        ReDim Preserve Yi(Nn)
        If Nn = 0 Then
            Picture1.Cls
            Picture1.ForeColor = RGB(200, 200, 200)
            Picture1.Line (Picture1.ScaleLeft, 0)-(Picture1.ScaleLeft + Picture1.ScaleWidth, 0)
            Picture1.Line (0, Picture1.ScaleTop)-(0, Picture1.ScaleTop + Picture1.ScaleHeight)
        End If
        Xi(Nn) = X
        Yi(Nn) = Y
        Picture1.DrawWidth = 4
        Picture1.PSet (Xi(Nn), Yi(Nn)), vbBlue    '画点
        Picture1.DrawWidth = 1
        Nn = Nn + 1
        ReDim Preserve Xi(Nn)
        ReDim Preserve Yi(Nn)
    End If
    If Button = 2 Then   '右键结束画点
        Nn = 0
        Erase Xi(), Yi()   '释放数组
        Line1.Visible = False
    End If
    
End Sub

Private Sub Picture1_MouseMove(Button As Integer, Shift As Integer, X As Single, Y As Single)
    If Nn = 0 Then
        ReDim Xi(Nn)
        ReDim Yi(Nn)
        
      Else
        Line1.X1 = Xi(Nn - 1)
        Line1.Y1 = Yi(Nn - 1)
        Line1.X2 = X
        Line1.Y2 = Y
        Line1.Visible = True
    End If
'    Xi(Nn) = X
'    Yi(Nn) = Y
'    If Nn > 1 Then 画样条
    Label1.Caption = "X=" & X & " Y=" & Y
End Sub

Private Sub Picture1_MouseUp(Button As Integer, Shift As Integer, X As Single, Y As Single)
    If Nn > 1 Then 画样条
End Sub