貝塞爾曲線及插值

1、貝塞爾曲線及插值 這里主要講一下如何在excel及vb中實(shí)現(xiàn)貝塞爾曲線插值，程序來源于互聯(lián)網(wǎng)（程序作者: 海底眼(Mr. Dragon Pan在excel中用宏實(shí)現(xiàn))，本文作為少量修改，方便在vb中調(diào)用，經(jīng)運(yùn)行證明是沒錯(cuò)的，下面程序可作成一個(gè)模塊放到vb或vba中調(diào)用： Excel的平滑線散點(diǎn)圖，可以根據(jù)兩組分別代表X-Y坐標(biāo)的散點(diǎn)數(shù)值產(chǎn)生曲線圖 但是，卻沒有提供這個(gè)曲線圖的公式，所以無法查找曲線上的點(diǎn)坐標(biāo) 后來我在以下這個(gè)網(wǎng)頁找到了詳細(xì)的說明和示例程序 . /Smooth_curve_bezier_example_file.zip . 根據(jù)其中采用的

2、算法，進(jìn)一步增添根據(jù)X坐標(biāo)求Y坐標(biāo)，或根據(jù)Y坐標(biāo)求X坐標(biāo)，更切合實(shí)際需求 這個(gè)自定義函數(shù)按照Excel的曲線算法(三次貝塞爾分段插值),計(jì)算平滑曲線上任意一點(diǎn)的點(diǎn)坐標(biāo) Excel的平滑曲線的大致算法是: 給出了兩組X-Y數(shù)值以后，每一對X-Y坐標(biāo)稱為節(jié)點(diǎn)，然后在每兩個(gè)節(jié)點(diǎn)之間畫出三次貝塞爾曲線(下面簡稱曲線) 貝塞爾曲線的算法網(wǎng)上有很多資源，這里不介紹了，只作簡單說明 每條曲線都由四個(gè)節(jié)點(diǎn)開始，計(jì)算出四個(gè)貝塞爾控制點(diǎn)，然后根據(jù)控制點(diǎn)畫出唯一一條曲線 假設(shè)曲線的源數(shù)據(jù)是節(jié)點(diǎn)1,節(jié)點(diǎn)2,節(jié)點(diǎn)3,節(jié)點(diǎn)4(Dot1,Dot2,Dot3,Dot4) 那么貝塞爾控制點(diǎn)的計(jì)算如下 Dot2是第一個(gè)控制點(diǎn),也

3、是曲點(diǎn)的起點(diǎn)，Dot3是第四個(gè)控制點(diǎn)也是曲線的終點(diǎn) 第二個(gè)控制點(diǎn)的位置是: 過第一個(gè)控制點(diǎn)(Dot2,起點(diǎn))，與Dot1, Dot3的連線平行，且與Dot2距離為 1/6 * 線段Dot1_Dot3的長度 假如是圖形的第一段曲線，取節(jié)點(diǎn)1,1,2,3進(jìn)行計(jì)算,即 Dot2 = Dot1 且第二個(gè)控制點(diǎn)與第一控制點(diǎn)距離取 1/3 * |Dot1_Dot3|,而不是1/6 * |Dot1_Dot3| 假如 1/2 * |Dot2_Dot3| 1/6 * |Dot1_Dot3| 那么第二個(gè)控制點(diǎn)與第一控制點(diǎn)距離取 1/2 * |Dot2_Dot3|,而不是1/6 * |Dot1_Dot3| 第三個(gè)控

4、制點(diǎn)的位置是: 過第四個(gè)控制點(diǎn)(Dot3,終點(diǎn))，與Dot2, Dot4的連線平行，且與Dot3距離為 1/6 * |Dot2_Dot4| 假如是圖形的最后一段曲線，取節(jié)點(diǎn)Last-2,Last-1,Last,Last進(jìn)行計(jì)算,即 Dot4 = Dot3 且第三個(gè)控制點(diǎn)與第四控制點(diǎn)距離取 1/3 * |Dot2_Dot4|,而不是1/6 * |Dot2_Dot4| 假如 1/2 * |Dot2_Dot3| =1 and =0 and =1Const Error10 = Error: known_value is not on the curve (defined by given known_

5、x and known_y)Const NoRoot = No RootConst MaxErr = 0.Const MaxLoop = 1000Dim SizeX, SizeY As Long 輸入?yún)^(qū)域的大小Dim Dot1 As Vector 輸入?yún)^(qū)域里面，用作計(jì)算貝塞爾控制點(diǎn)的四個(gè)節(jié)點(diǎn)Dim Dot2 As VectorDim Dot3 As VectorDim Dot4 As VectorDim BezierPt1 As Vector 生成貝塞爾曲線的四個(gè)貝塞爾控制點(diǎn)Dim BezierPt2 As VectorDim BezierPt3 As VectorDim BezierPt4

6、As VectorDim OffsetTo2 As Vector 第二,三個(gè)貝塞爾控制點(diǎn)跟起點(diǎn)，終點(diǎn)的距離關(guān)系Dim OffsetTo3 As VectorDim ValueType As Variant 輸入待查數(shù)值的類型,x代表輸入的是X坐標(biāo)，求對應(yīng)的Y坐標(biāo)Dim Interpol_here As Boolean 當(dāng)前分段曲線是否包含待查數(shù)值Dim key_value, a, b, c, d As Double 貝塞爾曲線插值多項(xiàng)式的系數(shù)Dim t1, t2, t3 As Variant 貝塞爾曲線插值多項(xiàng)式的根Dim a3, a2, a1, a0 As DoubleDim size%Pu

7、blic Sub befit(ByRef known_x() As Double, ByRef known_y() As Double, size As Integer, known_value As Double, result() As Variant, Optional StartKnot As Long = 1, Optional known_value_type As Variant = x)-子過程方便VB中調(diào)用-主程序開始，至少要輸入五個(gè)參數(shù)，第一個(gè)是X坐標(biāo)系列，然后是Y坐標(biāo)系列，第三個(gè)是坐標(biāo)點(diǎn)數(shù)，第四個(gè)是待查數(shù)值，第五個(gè)是返回值第六個(gè)參數(shù)是從哪一段曲線開始查找，如果曲線可以返回

8、多個(gè)值，那么分別指定起始節(jié)點(diǎn)就可以找出全部合要求的點(diǎn)第七個(gè)參數(shù)是待查數(shù)值的類型,x代表輸入x坐標(biāo)求對應(yīng)y坐標(biāo)，y則相反,t是直接輸入貝塞爾插值多項(xiàng)式的參數(shù)-Dim j As LongDim x1Value, y1Value, x2Value, y2Value, x3Value, y3Value As VariantDim ErrorMsg As VariantValueType = LCase(known_value_type) 待查數(shù)值的類型轉(zhuǎn)化為小寫，并賦值到全局變量ValueTypekey_value = known_value 待查數(shù)值賦值到全局變量key_valueErrorMsg

9、= ErrorCheck(known_x, known_y, StartKnot) 檢查輸入錯(cuò)誤If ErrorMsg NoError Then 有錯(cuò)誤就返回錯(cuò)誤信息，退出程序 result = Array(ErrorMsg, ErrorMsg, ErrorMsg, ErrorMsg, ErrorMsg, ErrorMsg) Exit SubEnd IfSizeX = UBound(known_x)For j = StartKnot To size SizeX - 1 從指定的節(jié)點(diǎn)開始，沒有指定節(jié)點(diǎn)就從1開始 Call FindFourDots(known_x, known_y, j) 找出輸

10、入X-Y點(diǎn)坐標(biāo)里面，應(yīng)該用于計(jì)算的四個(gè)結(jié)點(diǎn) Call FindFourBezierPoints(Dot1, Dot2, Dot3, Dot4) 根據(jù)四個(gè)結(jié)點(diǎn)計(jì)算四個(gè)貝塞爾控制點(diǎn) Call FindABCD 根據(jù)待查數(shù)值的類型，和貝塞爾控制點(diǎn),計(jì)算貝塞爾插值多項(xiàng)式的系數(shù) Call Find_t 檢查貝塞爾曲線是否包含待查數(shù)值 If Interpol_here = True Then Exit ForNext jIf Interpol_here = True Then 計(jì)算點(diǎn)坐標(biāo)，并返回 以下是由四個(gè)貝塞爾控制點(diǎn)決定的,貝塞爾曲線的參數(shù)方程 x1Value = (1 - t1) 3 * Bezie

11、rPt1.x + 3 * t1 * (1 - t1) 2 * BezierPt2.x + 3 * t1 2 * (1 - t1) * BezierPt3.x + t1 3 * BezierPt4.x y1Value = (1 - t1) 3 * BezierPt1.y + 3 * t1 * (1 - t1) 2 * BezierPt2.y + 3 * t1 2 * (1 - t1) * BezierPt3.y + t1 3 * BezierPt4.y x2Value = (1 - t2) 3 * BezierPt1.x + 3 * t2 * (1 - t2) 2 * BezierPt2.x +

12、 3 * t2 2 * (1 - t2) * BezierPt3.x + t2 3 * BezierPt4.x y2Value = (1 - t2) 3 * BezierPt1.y + 3 * t2 * (1 - t2) 2 * BezierPt2.y + 3 * t2 2 * (1 - t2) * BezierPt3.y + t2 3 * BezierPt4.y x3Value = (1 - t3) 3 * BezierPt1.x + 3 * t3 * (1 - t3) 2 * BezierPt2.x + 3 * t3 2 * (1 - t3) * BezierPt3.x + t3 3 *

13、BezierPt4.x y3Value = (1 - t3) 3 * BezierPt1.y + 3 * t3 * (1 - t3) 2 * BezierPt2.y + 3 * t3 2 * (1 - t3) * BezierPt3.y + t3 3 * BezierPt4.y result = Array(x1Value, y1Value, x2Value, y2Value, x3Value, y3Value)Else result = Array(Error10, Error10, Error10, Error10, Error10, Error10)End IfEnd Sub/*Func

14、tion ErrorCheck(ByRef known_x() As Double, ByRef known_y() As Double, StartKnot) As VariantErrorCheck = NoErrorSizeX = UBound(known_x) known_x.CountSizeY = UBound(known_y) known_y.CountIf SizeX SizeY Then 假如輸入的X坐標(biāo)數(shù)目不等于Y坐標(biāo)數(shù)目ErrorCheck = Error1Exit FunctionEnd IfIf SizeX 3 Then 輸入的X-Y坐標(biāo)對少于三個(gè)ErrorCheck

15、 = Error2Exit FunctionEnd IfIf (StartKnot = SizeX) Then 指定的起始節(jié)點(diǎn)超出范圍ErrorCheck = Error3Exit FunctionEnd IfIf (ValueType x And ValueType y And ValueType t) Then 輸入的待查數(shù)值類型不是x, y, tErrorCheck = Error4Exit FunctionEnd IfIf (ValueType = t And key_value 1) Or (ValueType = t And key_value 0) Then t 類型的范圍是0-

16、1ErrorCheck = Error5Exit FunctionEnd IfEnd Function/*Sub FindFourDots(ByRef known_x() As Double, ByRef known_y() As Double, j) 根據(jù)X-Y數(shù)值，及起始節(jié)點(diǎn)，找出用于計(jì)算的四個(gè)結(jié)點(diǎn)坐標(biāo) If j = 1 Then 第一個(gè)結(jié)點(diǎn) Dot2 = Dot1 Dot1.x = known_x(1) Dot1.y = known_y(1) Else Dot1.x = known_x(j - 1) Dot1.y = known_y(j - 1) End If Dot2.x = know

17、n_x(j) Dot2.y = known_y(j) Dot3.x = known_x(j + 1) Dot3.y = known_y(j + 1) If j = SizeX - 1 Then 最后一個(gè)結(jié)點(diǎn) Dot4 = Dot3 Dot4.x = Dot3.x Dot4.y = Dot3.y Else Dot4.x = known_x(j + 2) Dot4.y = known_y(j + 2) End IfEnd Sub/*Sub FindFourBezierPoints(Dot1 As Vector, Dot2 As Vector, Dot3 As Vector, Dot4 As Vec

18、tor)Dim d12, d23, d34, d13, d14, d24 As Doubled12 = DistAtoB(Dot1, Dot2) 計(jì)算平面坐標(biāo)系上的兩點(diǎn)距離d23 = DistAtoB(Dot2, Dot3)d34 = DistAtoB(Dot3, Dot4)d13 = DistAtoB(Dot1, Dot3)d14 = DistAtoB(Dot1, Dot4)d24 = DistAtoB(Dot2, Dot4)BezierPt1 = Dot2BezierPt4 = Dot3OffsetTo2 = AsubB(Dot3, Dot1) 向量減法OffsetTo3 = AsubB(

19、Dot2, Dot4)If (d13 / 6 d23 / 2) And (d24 / 6 d23 / 2) Then If (Dot1.x Dot2.x Or Dot1.y Dot2.y) Then OffsetTo2 = AmultiF(OffsetTo2, 1 / 6) If (Dot1.x = Dot2.x And Dot1.y = Dot2.y) Then OffsetTo2 = AmultiF(OffsetTo2, 1 / 3) If (Dot3.x Dot4.x Or Dot3.y Dot4.y) Then OffsetTo3 = AmultiF(OffsetTo3, 1 / 6)

20、 If (Dot3.x = Dot4.x And Dot3.y = Dot4.y) Then OffsetTo3 = AmultiF(OffsetTo3, 1 / 3)ElseIf (d13 / 6 = d23 / 2) And (d24 / 6 = d23 / 2) Then OffsetTo2 = AmultiF(OffsetTo2, d23 / 12) OffsetTo3 = AmultiF(OffsetTo3, d23 / 12)ElseIf (d13 / 6 = d23 / 2) Then OffsetTo2 = AmultiF(OffsetTo2, d23 / 2 / d13) O

21、ffsetTo3 = AmultiF(OffsetTo3, d23 / 2 / d13)ElseIf (d24 / 6 = d23 / 2) Then OffsetTo2 = AmultiF(OffsetTo2, d23 / 2 / d24) OffsetTo3 = AmultiF(OffsetTo3, d23 / 2 / d24)End IfBezierPt2 = AaddB(BezierPt1, OffsetTo2) 向量加法BezierPt3 = AaddB(BezierPt4, OffsetTo3)End Sub/*Function DistAtoB(dota As Vector, d

22、otb As Vector) As DoubleDistAtoB = (dota.x - dotb.x) 2 + (dota.y - dotb.y) 2) 0.5End FunctionFunction AaddB(dota As Vector, dotb As Vector) As VectorAaddB.x = dota.x + dotb.xAaddB.y = dota.y + dotb.yEnd FunctionFunction AsubB(dota As Vector, dotb As Vector) As VectorAsubB.x = dota.x - dotb.xAsubB.y

23、= dota.y - dotb.yEnd FunctionFunction AmultiF(dota As Vector, MultiFactor As Double) As VectorAmultiF.x = dota.x * MultiFactorAmultiF.y = dota.y * MultiFactorEnd Function/*Sub FindABCD()If ValueType = x Then 參數(shù)類型是x, 需要解參數(shù)方程 f(t) = x，這里設(shè)定參數(shù)方程的系數(shù)a = -BezierPt1.x + 3 * BezierPt2.x - 3 * BezierPt3.x + B

24、ezierPt4.xb = 3 * BezierPt1.x - 6 * BezierPt2.x + 3 * BezierPt3.xc = -3 * BezierPt1.x + 3 * BezierPt2.xd = BezierPt1.x - key_valueEnd IfIf ValueType = y Then 參數(shù)類型是x, 需要解參數(shù)方程 f(t) = y，這里設(shè)定參數(shù)方程的系數(shù)a = -BezierPt1.y + 3 * BezierPt2.y - 3 * BezierPt3.y + BezierPt4.yb = 3 * BezierPt1.y - 6 * BezierPt2.y +

25、3 * BezierPt3.yc = -3 * BezierPt1.y + 3 * BezierPt2.yd = BezierPt1.y - key_valueEnd IfEnd Sub/*Sub Find_t() 計(jì)算當(dāng) f(t) = 待查數(shù)值時(shí), t應(yīng)該是什么數(shù)值Dim tArr As VariantInterpol_here = TrueIf ValueType = t Then 待查數(shù)值類型為t,那么無需計(jì)算 t1 = key_value t2 = key_value t3 = key_value Exit SubEnd IftArr = Solve_Order3_Equation(a

26、, b, c, d) 否則，解三次貝塞爾參數(shù)方程 f(t) = 待查數(shù)值t1 = tArr(1) 解得方程的三個(gè)根t2 = tArr(2)t3 = tArr(3)If (t1 1 Or t1 1 Or t2 1 Or t3 0) Then t3 = NoRootEnd IfIf (IsNumeric(t1) = False And IsNumeric(t2) = False And IsNumeric(t3) = False) Then Interpol_here = FalseEnd If 三個(gè)根都不合要求，代表曲線上沒有包含待查數(shù)值的點(diǎn)If (t1 = NoRoot And t2 NoRo

27、ot) Then 至少有一個(gè)根，則用它代替NoRoot的結(jié)果,方便Excel畫圖 t1 = t2End IfIf (t1 = NoRoot And t3 NoRoot) Then t1 = t3End IfIf (t2 = NoRoot) Then t2 = t1If (t3 = NoRoot) Then t3 = t1End Sub/*. 牛頓法解三次方程，先求解方程的導(dǎo)函數(shù)，得到方程的拐點(diǎn)(導(dǎo)數(shù)等于0的點(diǎn)) 然后分三段用迭代法分別求三個(gè)根.Public Function Solve_Order3_Equation(p3, p2, p1, P0, Optional Starting As D

28、ouble = -#, Optional Ending As Double = #) As VariantDim Two_X, TurningPoint, x1, x2, x3 As VariantDim x As Doublea3 = p3a2 = p2a1 = p1a0 = P0x1 = NoRootx2 = NoRootx3 = NoRootx1 = Newton_Solve(Starting)If a3 = 0 Then 如果三次方程沒有三次項(xiàng) Two_X = Solve_Order2_Equation(a2, a1, a0) 解釋法直接求二次方程的解 x1 = Two_X(1) x2

29、 = Two_X(2)Else TurningPoint = Solve_Order2_Equation(3 * a3, 2 * a2, 1 * a1) 求解 f(t) = 0 If (TurningPoint(1) = NoRoot And TurningPoint(2) = NoRoot) Then 分段求根 x = 0 x1 = Newton_Solve(x) ElseIf (TurningPoint(1) NoRoot And TurningPoint(2) = NoRoot) Then If f_x(Starting) * f_x(TurningPoint(1) 0 Then x =

30、 (Starting + TurningPoint(1) / 2 x1 = Newton_Solve(x) End If If f_x(TurningPoint(2) * f_x(Ending) 0 Then x = (TurningPoint(2) + Ending) / 2 x3 = Newton_Solve(x) End If ElseIf (TurningPoint(1) NoRoot And TurningPoint(2) NoRoot) Then If f_x(Starting) * f_x(TurningPoint(1) 0 Then x = (Starting + TurningPoint(1) / 2 x1 = Newton_Solve(x) End If If f_x(TurningPoint(1) * f_x(TurningPoint(2) 0 Then x = (TurningPoint(1) + TurningPoint(2) / 2 x2 = Newton_Solve(x) End If If f_x(TurningPoint(2) * f_x(Ending) 0 Then x = (TurningPoint(2) + Ending) / 2 x3 = Newton_Solve(x) End If End IfEnd IfSolve