Splines in VB6

Preview:

Option Explicit
'
Public Type P_Type
    X As Double         ' Coordinate x dei punti.
    Y As Double         ' Coordinate y dei punti.
    'z As Double         ' Coordinate z dei punti.
End Type
Public Sub Bezier_C(Pi() As P_Type, Pc() As P_Type)
'
'   Ritorna, nel vettore Pc(), i valori della curva di Bezier calcolata
'   al valore u (0 <= u <= 1). La curva e' calcolata in modo
'   parametrico con il valore 0 di u corrispondente a Pc(0)
'   ed il valore 1 corrispondente a Pc(NPC_1).
'   Questo algoritmo ricalca la forma classica del polinomio
'   di Bernstein.
'
   Dim I&, K&, NPI_1&, NPC_1&, NF&, u#, BF#
'
    NPI_1 = UBound(Pi)
    NPC_1 = UBound(Pc)
    'NF = Prodotto(NPI_1)
'
    For I = 0 To NPC_1
        u = CDbl(I) / CDbl(NPC_1)
        Pc(I).X = 0#
        Pc(I).Y = 0#
        'Pc(I).z = 0#
        For K = 0 To NPI_1
            'BF = NF * (u ^ K) * ((1 - u) ^ (NPI_1 - K)) _
                / (Prodotto(K) * Prodotto(NPI_1 - K))
            BF = Prodotto(NPI_1, K + 1) * (u ^ K) * ((1# - u) ^ (NPI_1 - K)) _
               / Prodotto(NPI_1 - K)
            Pc(I).X = Pc(I).X + Pi(K).X * BF
            Pc(I).Y = Pc(I).Y + Pi(K).Y * BF
            'Pc(I).z = Pc(I).z + Pi(K).z * BF
        Next K
    Next I
'
'
'
End Sub
Public Sub Bezier_1(Pi() As P_Type, Pc() As P_Type)
'
'   Ritorna, nel vettore Pc(), i valori della curva di Bezier.
'   La curva e' calcolata in modo parametrico (0 <= u < 1)
'   con il valore 0 di u corrispondente a Pc(0);
'   Attenzione: il punto Pc(NPC_1), corrispondente al valore u = 1,
'               non puo' essere calcolato.
'
'   Parametri:
'       Pi(0 to NPI - 1):   Vettore dei punti, dati, da
'                           approssimare.
'       Pc(0 to NPC - 1):   Vettore dei punti, calcolati,
'                           della curva approssimante.
'
   Dim I&, K&, NPI_1&, NPC_1&
   Dim u#, u_1#, ue#, u_1e#, BF#
'
    NPI_1 = UBound(Pi) ' N. di punti da approssimare - 1.
    NPC_1 = UBound(Pc) ' N. di punti sulla curva - 1.
'
    ' La curva inizia sempre da Pi(0) -> u = 0:
    Pc(0).X = Pi(0).X
    Pc(0).Y = Pi(0).Y
    'Pc(0).z = Pi(0).z
'
    For I = 1 To NPC_1 - 1
        u = CDbl(I) / CDbl(NPC_1)
        ue = 1#
        u_1 = 1# - u
        u_1e = u_1 ^ NPI_1
'
        Pc(I).X = 0#
        Pc(I).Y = 0#
        'Pc(I).z = 0#
        For K = 0 To NPI_1
            BF = Prodotto(NPI_1, K + 1) * ue * u_1e / Prodotto(NPI_1 - K)
            Pc(I).X = Pc(I).X + Pi(K).X * BF
            Pc(I).Y = Pc(I).Y + Pi(K).Y * BF
            'Pc(I).z = Pc(I).z + Pi(K).z * BF
'
            ue = ue * u
            u_1e = u_1e / u_1
        Next K
    Next I
'
    ' La curva finisce sempre su Pi(NPI_1) -> u = 1:
    Pc(NPC_1).X = Pi(NPI_1).X
    Pc(NPC_1).Y = Pi(NPI_1).Y
    'Pc(NPC_1).z = Pi(NPI_1).z
'
'
'
End Sub
Public Sub Bezier(Pi() As P_Type, Pc() As P_Type)
'
'   Ritorna, nel vettore Pc(), i valori della curva di Bezier.
'   La curva e' calcolata in modo parametrico (0 <= u < 1)
'   con il valore 0 di u corrispondente a Pc(0);
'
'   Questa versione elimina alcuni problemi di "underflow"
'   e di "overflow" presentati dalla Bezier_1 e dalla Bezier_C.
'
'   Parametri:
'       Pi(0 to NPI - 1):   Vettore dei punti, dati, da
'                           approssimare.
'       Pc(0 to NPC - 1):   Vettore dei punti, calcolati,
'                           della curva approssimante.
'
   Dim I&, K&, NPI_1&, NPC_1&
   Dim u#, u_1#, ue#, BF#
   Static NPI_1_O&, CB_Tav#()
'
    NPI_1 = UBound(Pi) ' N. di punti da approssimare - 1 (deve essere 2 <= NPI_1 <= 1029).
    NPC_1 = UBound(Pc) ' N. di punti sulla curva - 1.
'
    If NPI_1_O <> NPI_1 Then
        ' Prepara la tavola dei coefficienti binomiali:
        ReDim CB_Tav#(0 To NPI_1)
        For K = 0 To NPI_1
            CB_Tav(K) = rncr(NPI_1, K)
        Next K
'
        NPI_1_O = NPI_1
    End If
'
    ' La curva inizia sempre da Pi(0) -> u = 0:
    Pc(0).X = Pi(0).X
    Pc(0).Y = Pi(0).Y
    'Pc(0).z = Pi(0).z
'
    For I = 1 To NPC_1 - 1
        u = CDbl(I) / CDbl(NPC_1)
        ue = 1#
        u_1 = 1# - u
'
        Pc(I).X = 0#
        Pc(I).Y = 0#
        'Pc(I).z = 0#
        For K = 0 To NPI_1
            BF = CB_Tav(K) * ue * u_1 ^ (NPI_1 - K)
'
            Pc(I).X = Pc(I).X + Pi(K).X * BF
            Pc(I).Y = Pc(I).Y + Pi(K).Y * BF
            'Pc(I).z = Pc(I).z + Pi(K).z * BF
'
            ue = ue * u
        Next K
    Next I
'
    ' La curva finisce sempre su Pi(NPI_1) -> u = 1:
    Pc(NPC_1).X = Pi(NPI_1).X
    Pc(NPC_1).Y = Pi(NPI_1).Y
    'Pc(NPC_1).z = Pi(NPI_1).z
'
'
'
End Sub
Public Sub Bezier_P(Pi() As P_Type, Pc() As P_Type)
'
'   Ritorna, nel vettore Pc(), i valori della curva di Bezier calcolata
'   al valore u (0 <= u < 1). La curva e' calcolata in modo
'   parametrico con il valore 0 di u corrispondente a Pc(0);
'   Attenzione: il punto Pc(NPC_1), corrispondente al valore u = 1,
'               non puo' essere calcolato.
'
'   Questo algoritmo (tratto da una pubblicazione di P. Bourke
'   e tradotto dal C) e' particolarmente interessante, in quanto
'   evita l' uso dei fattoriali della forma normale.
'
    Dim K&, I&, KN&, NPI_1&, NPC_1&, NN&, NKN&
    Dim u#, uk#, unk#, Blend#
'
    NPI_1 = UBound(Pi)
    NPC_1 = UBound(Pc)
'
    For I = 0 To NPC_1 - 1
        u = CDbl(I) / CDbl(NPC_1)
        uk = 1#
        unk = (1# - u) ^ NPI_1
'
        Pc(I).X = 0#
        Pc(I).Y = 0#
        'Pc(I).z = 0#
'
        For K = 0 To NPI_1
            NN = NPI_1
            KN = K
            NKN = NPI_1 - K
            Blend = uk * unk
            uk = uk * u
            unk = unk / (1# - u)
            Do While NN >= 1
                Blend = Blend * CDbl(NN)
                NN = NN - 1
                If KN > 1 Then
                    Blend = Blend / CDbl(KN)
                    KN = KN - 1
                End If
                If NKN > 1 Then
                    Blend = Blend / CDbl(NKN)
                    NKN = NKN - 1
                End If
            Loop
'
            Pc(I).X = Pc(I).X + Pi(K).X * Blend
            Pc(I).Y = Pc(I).Y + Pi(K).Y * Blend
            'Pc(I).z = Pc(I).z + Pi(K).z * Blend
        Next K
    Next I
'
    ' La curva finisce sempre su Pi(NPI_1) -> u = 1:
    Pc(NPC_1).X = Pi(NPI_1).X
    Pc(NPC_1).Y = Pi(NPI_1).Y
    'Pc(NPC_1).z = Pi(NPI_1).z
'
'
'
End Sub
Private Function Prodotto(ByVal N2&, Optional ByVal N1& = 2) As Double
'
'   Ritorna il prodotto dei numeri, consecutivi, interi e positivi,
'   da N1 a N2 (0 < N1 <= N2). Se N1 > N2 ritorna 1.
'   Se N1 manca, ritorna il Fattoriale di N2; in questo caso puo'
'   anche essere N2 = 0 perche', per definizione, e' 0! = 1:
'
    Dim Pr#, I&
'
    Pr = 1#
    For I = N1 To N2
        Pr = Pr * CDbl(I)
    Next I
'
    Prodotto = Pr
'
'
'
End Function
Public Sub B_Spline(Pi() As P_Type, ByVal NK&, Pc() As P_Type)
'
'   Ritorna, nel vettore Pc(), i valori della curva B-Spline.
'   La curva e' calcolata in modo parametrico (0 <= u <= 1)
'   con il valore 0 di u corrispondente a Pc(0) ed il valore
'   1 corrispondente a Pc(NPC_1).
'
'   Parametri:
'       Pi(0 to NPI - 1):   Vettore dei punti, dati, da
'                           approssimare.
'       Pc(0 to NPC - 1):   Vettore dei punti, calcolati,
'                           della curva approssimante.
'       NK:                 Numero di nodi della curva
'                           approssimante:
'                           NK = 2    -> segmenti di retta.
'                           NK = 3    -> curve quadratiche.
'                           ..   .       ..................
'                           NK = NPI  -> splines di Bezier.

    Dim NPI_1&, NPC_1&, I&, J&, tmax#, u#, ut#, bn#()
    Const Eps = 0.0000001
'
    NPI_1 = UBound(Pi)  ' N. di punti da approssimare - 1.
    NPC_1 = UBound(Pc)  ' N. di punti sulla curva - 1.
    tmax = NPI_1 - NK + 2
'
    ' La curva inizia sempre da Pi(0) -> u = 0:
    Pc(0).X = Pi(0).X
    Pc(0).Y = Pi(0).Y
'
    For I = 1 To NPC_1 - 1
        u = CDbl(I) / CDbl(NPC_1)
        ut = u * tmax
        If Abs(ut - CDbl(NPI_1 + NK - 2)) <= Eps Then
            Pc(I).X = Pi(NPI_1).X
            Pc(I).Y = Pi(NPI_1).Y
        Else
            Call B_Basis(NPI_1, ut, NK, bn())
            Pc(I).X = 0#
            Pc(I).Y = 0#
            For J = 0 To NPI_1
                Pc(I).X = Pc(I).X + bn(J) * Pi(J).X
                Pc(I).Y = Pc(I).Y + bn(J) * Pi(J).Y
            Next J
        End If
    Next I
'
    ' La curva finisce sempre su Pi(NPI_1) -> u = 1:
    Pc(NPC_1).X = Pi(NPI_1).X
    Pc(NPC_1).Y = Pi(NPI_1).Y
'
'
'
End Sub
Private Sub B_Basis(ByVal NPI_1&, ByVal ut#, ByVal K&, bn#())
'
'   Compute the basis function (also called weight)
'   for the B-Spline approximation curve:
'
    Dim NT&, I&, J&
    Dim b0#, b1#, bl0#, bl1#, bu0#, bu1#
    ReDim bn#(0 To NPI_1 + 1), bn0#(0 To NPI_1 + 1), t#(0 To NPI_1 + K + 1)
'
    NT = NPI_1 + K + 1
    For I = 0 To NT
        If (I < K) Then t(I) = 0#
        If ((I >= K) And (I <= NPI_1)) Then t(I) = CDbl(I - K + 1)
        If (I > NPI_1) Then t(I) = CDbl(NPI_1 - K + 2)
    Next I
    For I = 0 To NPI_1
        bn0(I) = 0#
        If ((ut >= t(I)) And (ut < t(I + 1))) Then bn0(I) = 1#
        If ((t(I) = 0#) And (t(I + 1) = 0#)) Then bn0(I) = 0#
    Next I
'
    For J = 2 To K
        For I = 0 To NPI_1
            bu0 = (ut - t(I)) * bn0(I)
            bl0 = t(I + J - 1) - t(I)
            If (bl0 = 0#) Then
                b0 = 0#
            Else
                b0 = bu0 / bl0
            End If
            bu1 = (t(I + J) - ut) * bn0(I + 1)
            bl1 = t(I + J) - t(I + 1)
            If (bl1 = 0#) Then
                b1 = 0#
            Else
                b1 = bu1 / bl1
            End If
            bn(I) = b0 + b1
        Next I
        For I = 0 To NPI_1
            bn0(I) = bn(I)
        Next I
    Next J
'
'
'
End Sub
Public Sub C_Spline(Pi() As P_Type, Pc() As P_Type)
'
'   Ritorna, nel vettore Pc(), i valori della curva C-Spline.
'   La curva e' calcolata in modo parametrico (0 <= u <= 1)
'   con il valore 0 di u corrispondente a Pc(0) ed il valore
'   1 corrispondente a Pc(NPC_1).
'
'   Parametri:
'       Pi(0 to NPI - 1):   Vettore dei punti, dati, da
'                           interpolare.
'       Pc(0 to NPC - 1):   Vettore dei punti, calcolati,
'                           della curva interpolante.
'
    Dim NPI_1&, NPC_1&, I&, J&
    Dim u#, ui#, uui#
    Dim cof() As P_Type
'
    NPI_1 = UBound(Pi)      ' N. di punti da interpolare - 1.
    NPC_1 = UBound(Pc)      ' N. di punti sulla curva - 1.
'
    Call Find_CCof(Pi(), NPI_1 + 1, cof())
'
    ' La curva inizia sempre da Pi(0) -> u = 0:
    Pc(0).X = Pi(0).X
    Pc(0).Y = Pi(0).Y
'
    For I = 1 To NPC_1 - 1
        u = CDbl(I) / CDbl(NPC_1)
        J = Int(u * CDbl(NPI_1)) + 1
        If (J > (NPI_1)) Then J = NPI_1
'
        ui = CDbl(J - 1) / CDbl(NPI_1)
        uui = u - ui
'
        Pc(I).X = cof(4, J).X * uui ^ 3 + cof(3, J).X * uui ^ 2 + cof(2, J).X * uui + cof(1, J).X
        Pc(I).Y = cof(4, J).Y * uui ^ 3 + cof(3, J).Y * uui ^ 2 + cof(2, J).Y * uui + cof(1, J).Y
    Next I
'
    ' La curva finisce sempre su Pi(NPI_1) -> u = 1:
    Pc(NPC_1).X = Pi(NPI_1).X
    Pc(NPC_1).Y = Pi(NPI_1).Y
'
'
'
End Sub
Private Function rncr(ByVal N&, ByVal K&) As Double
'
'   Calcola i coefficienti binomiali Cn,k come:
'    rncr = N! / (K! * (N - K)!)
'
'   Nota: La funzione ha senso solo per 0 < N, K <= N
'         e 0 <= K.  Nessun errore viene segnalato.
'
    Dim I&, rncr_T#
'
    If ((N < 1) Or (K < 1) Or (N = K)) Then
        rncr = 1#
'
    Else
        rncr_T = 1#
        For I = 1 To N - K
            rncr_T = rncr_T * (1# + CDbl(K) / CDbl(I))
        Next I
'
        rncr = rncr_T
    End If
'
'
'
End Function
Public Sub T_Spline(Pi() As P_Type, ByVal VZ&, Pc() As P_Type)
'
'   Ritorna, nel vettore Pc(), i valori della curva T-Spline.
'   La curva e' calcolata in modo parametrico (0 <= u <= 1)
'   con il valore 0 di u corrispondente a Pc(0) ed il valore
'   1 corrispondente a Pc(NPC_1).
'
'   Parametri:
'       Pi(0 to NPI - 1):   Vettore dei punti, dati, da
'                           interpolare.
'       Pc(0 to NPC - 1):   Vettore dei punti, calcolati,
'                           della curva interpolante.
'       VZ:                 Valore di tensione della curva
'                           (1 <= VZ <= 100): valori grandi
'                           di VZ appiattiscono la curva.
'
    Dim NPI_1&, NPC_1&, I&, J&
    Dim H#, z#, z2i#, szh#, u#, u0#, u1#, du1#, du0#
    Dim s() As P_Type
'
    NPI_1 = UBound(Pi)      ' N. di punti da interpolare - 1.
    NPC_1 = UBound(Pc)      ' N. di punti sulla curva - 1.
    z = CDbl(VZ)
'
    Call Find_TCof(Pi(), NPI_1 + 1, s(), z)
'
    ' La curva inizia sempre da Pi(0) -> u = 0:
    Pc(0).X = Pi(0).X
    Pc(0).Y = Pi(0).Y
'
    H = 1# / CDbl(NPI_1)
    szh = Sinh(z * H)
    z2i = 1# / z / z
    For I = 1 To NPC_1 - 1
        u = CDbl(I) / CDbl(NPC_1)
        J = Int(u * CDbl(NPI_1)) + 1
        If (J > (NPI_1)) Then J = NPI_1
'
        u0 = CDbl(J - 1) / CDbl(NPI_1)
        u1 = CDbl(J) / CDbl(NPI_1)
        du1 = u1 - u
        du0 = u - u0
'
        Pc(I).X = s(J).X * z2i * Sinh(z * du1) / szh + (Pi(J - 1).X - s(J).X * z2i) * du1 / H
        Pc(I).X = Pc(I).X + s(J + 1).X * z2i * Sinh(z * du0) / szh + (Pi(J).X - s(J + 1).X * z2i) * du0 / H
    
        Pc(I).Y = s(J).Y * z2i * Sinh(z * du1) / szh + (Pi(J - 1).Y - s(J).Y * z2i) * du1 / H
        Pc(I).Y = Pc(I).Y + s(J + 1).Y * z2i * Sinh(z * du0) / szh + (Pi(J).Y - s(J + 1).Y * z2i) * du0 / H
    Next I
'
    ' La curva finisce sempre su Pi(NPI_1) -> u = 1:
    Pc(NPC_1).X = Pi(NPI_1).X
    Pc(NPC_1).Y = Pi(NPI_1).Y
'
'
'
End Sub
Private Sub Find_TCof(Pi() As P_Type, ByVal NPI&, s() As P_Type, ByVal z#)
'
'   Find the coefficients of the T-Spline
'   using constant interval:
'
    Dim I&, H#, a0#, b0#, zh#, z2#
'
    ReDim s(1 To NPI) As P_Type, f(1 To NPI) As P_Type
    ReDim a(1 To NPI) As Double, B(1 To NPI) As Double, C(1 To NPI) As Double
'
    H = 1# / CDbl(NPI - 1)
    zh = z * H
    a0 = 1# / H - z / Sinh(zh)
    b0 = z * 2# * Cosh(zh) / Sinh(zh) - 2# / H
    For I = 1 To NPI - 2
        a(I) = a0
        B(I) = b0
        C(I) = a0
    Next I
'
    z2 = z * z / H
    For I = 1 To NPI - 2
        f(I).X = (Pi(I + 1).X - 2# * Pi(I).X + Pi(I - 1).X) * z2
        f(I).Y = (Pi(I + 1).Y - 2# * Pi(I).Y + Pi(I - 1).Y) * z2
    Next I
'
    Call TRIDAG(a(), B(), C(), f(), s(), NPI - 2)
    For I = 1 To NPI - 2
        s(NPI - I).X = s(NPI - 1 - I).X
        s(NPI - I).Y = s(NPI - 1 - I).Y
    Next I
'
    s(1).X = 0#
    s(NPI).X = 0#
    s(1).Y = 0#
    s(NPI).Y = 0#
'
'
'
End Sub
Private Sub Find_CCof(Pi() As P_Type, ByVal NPI&, cof() As P_Type)
'
'   Find the coefficients of the cubic spline
'   using constant interval parameterization:
'
    Dim I&, H#
'
    ReDim s(1 To NPI) As P_Type, f(1 To NPI) As P_Type, cof(1 To 4, 1 To NPI) As P_Type
    ReDim a(1 To NPI) As Double, B(1 To NPI) As Double, C(1 To NPI) As Double
'
    H = 1# / CDbl(NPI - 1)
    For I = 1 To NPI - 2
        a(I) = 1#
        B(I) = 4#
        C(I) = 1#
    Next I
'
    For I = 1 To NPI - 2
        f(I).X = 6# * (Pi(I + 1).X - 2# * Pi(I).X + Pi(I - 1).X) / H / H
        f(I).Y = 6# * (Pi(I + 1).Y - 2# * Pi(I).Y + Pi(I - 1).Y) / H / H
    Next I
'
    Call TRIDAG(a(), B(), C(), f(), s(), NPI - 2)
    For I = 1 To NPI - 2
        s(NPI - I).X = s(NPI - 1 - I).X
        s(NPI - I).Y = s(NPI - 1 - I).Y
    Next I
'
    s(1).X = 0#
    s(NPI).X = 0#
    s(1).Y = 0#
    s(NPI).Y = 0#
    For I = 1 To NPI - 1
        cof(4, I).X = (s(I + 1).X - s(I).X) / 6# / H
        cof(4, I).Y = (s(I + 1).Y - s(I).Y) / 6# / H
        cof(3, I).X = s(I).X / 2#
        cof(3, I).Y = s(I).Y / 2#
        cof(2, I).X = (Pi(I).X - Pi(I - 1).X) / H - (2# * s(I).X + s(I + 1).X) * H / 6#
        cof(2, I).Y = (Pi(I).Y - Pi(I - 1).Y) / H - (2# * s(I).Y + s(I + 1).Y) * H / 6#
        cof(1, I).X = Pi(I - 1).X
        cof(1, I).Y = Pi(I - 1).Y
    Next I
'
'
'
End Sub
Private Sub TRIDAG(a#(), B#(), C#(), f() As P_Type, s() As P_Type, ByVal NPI_2&)
'
'   Solves the tridiagonal linear system of equations:
'
    Dim J&, bet#
    ReDim gam#(1 To NPI_2)
'
    If B(1) = 0 Then Exit Sub
'
    bet = B(1)
    s(1).X = f(1).X / bet
    s(1).Y = f(1).Y / bet
    For J = 2 To NPI_2
        gam(J) = C(J - 1) / bet
        bet = B(J) - a(J) * gam(J)
        If (bet = 0) Then Exit Sub
        s(J).X = (f(J).X - a(J) * s(J - 1).X) / bet
        s(J).Y = (f(J).Y - a(J) * s(J - 1).Y) / bet
    Next J
'
    For J = NPI_2 - 1 To 1 Step -1
        s(J).X = s(J).X - gam(J + 1) * s(J + 1).X
        s(J).Y = s(J).Y - gam(J + 1) * s(J + 1).Y
    Next J
'
'
'
End Sub
Private Function Cosh(ByVal z As Double) As Double
'
'   Ritorna il coseno iperbolico di z#:
'
    Cosh = (Exp(z) + Exp(-z)) / 2#
'
'
'
End Function
Private Function Sinh(ByVal z As Double) As Double
'
'   Ritorna il seno iperbolico di z#:
'
    Sinh = (Exp(z) - Exp(-z)) / 2#
'
'
'
End Function