贝塞尔曲线及插值.doc

贝塞尔曲线及插值 2010-07-01 21:41 贝塞尔曲线介绍可参考中文维基百科，图文并茂，这里就不啰嗦了 /zh-cn/%E8%B2%9D%E8%8C%B2%E6%9B%B2%E7%B7%9A 这里主要讲一下如何在excel及vb中实现贝塞尔曲线插值，程序来源于互联网（程序作者: 海底眼(Mr. Dragon Pan在excel中用宏实现)，本文作为少量修改，方便在vb中调用，经运行证明是没错的，下面程序可作成一个模块放到vb或vba中调用： - Excel的平滑线散点图，可以根据两组分别代表X-Y坐标的散点数值产生曲线图 但是，却没有提供这个曲线图的公式，所以无法查找曲线上的点坐标 后来我在以下这个网页找到了详细的说明和示例程序 . /Smooth_curve_bezier_example_file.zip . 根据其中采用的算法，进一步增添根据X坐标求Y坐标，或根据Y坐标求X坐标，更切合实际需求 这个自定义函数按照Excel的曲线算法(三次贝塞尔分段插值),计算平滑曲线上任意一点的点坐标 Excel的平滑曲线的大致算法是: 给出了两组X-Y数值以后，每一对X-Y坐标称为节点，然后在每两个节点之间画出三次贝塞尔曲线(下面简称曲线) 贝塞尔曲线的算法网上有很多资源，这里不介绍了，只作简单说明 每条曲线都由四个节点开始，计算出四个贝塞尔控制点，然后根据控制点画出唯一一条曲线 假设曲线的源数据是节点1,节点2,节点3,节点4(Dot1,Dot2,Dot3,Dot4) 那么贝塞尔控制点的计算如下 Dot2是第一个控制点,也是曲点的起点，Dot3是第四个控制点也是曲线的终点 第二个控制点的位置是: 过第一个控制点(Dot2,起点)，与Dot1, Dot3的连线平行，且与Dot2距离为 1/6 * 线段Dot1_Dot3的长度 假如是图形的第一段曲线，取节点1,1,2,3进行计算,即 Dot2 = Dot1 且第二个控制点与第一控制点距离取 1/3 * |Dot1_Dot3|,而不是1/6 * |Dot1_Dot3| 假如 1/2 * |Dot2_Dot3| 1/6 * |Dot1_Dot3| 那么第二个控制点与第一控制点距离取 1/2 * |Dot2_Dot3|,而不是1/6 * |Dot1_Dot3| 第三个控制点的位置是: 过第四个控制点(Dot3,终点)，与Dot2, Dot4的连线平行，且与Dot3距离为 1/6 * |Dot2_Dot4| 假如是图形的最后一段曲线，取节点Last-2,Last-1,Last,Last进行计算,即 Dot4 = Dot3 且第三个控制点与第四控制点距离取 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_x and known_y)Const NoRoot = No RootConst MaxErr = 0.00000001Const MaxLoop = 1000Dim SizeX, SizeY As Long 输入区域的大小Dim Dot1 As Vector 输入区域里面，用作计算贝塞尔控制点的四个节点Dim Dot2 As VectorDim Dot3 As VectorDim Dot4 As VectorDim BezierPt1 As Vector 生成贝塞尔曲线的四个贝塞尔控制点Dim BezierPt2 As VectorDim BezierPt3 As VectorDim BezierPt4 As VectorDim OffsetTo2 As Vector 第二,三个贝塞尔控制点跟起点，终点的距离关系Dim OffsetTo3 As VectorDim ValueType As Variant 输入待查数值的类型,x代表输入的是X坐标，求对应的Y坐标Dim Interpol_here As Boolean 当前分段曲线是否包含待查数值Dim key_value, a, b, c, d As Double 贝塞尔曲线插值多项式的系数Dim t1, t2, t3 As Variant 贝塞尔曲线插值多项式的根Dim a3, a2, a1, a0 As DoubleDim size%Public 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中调用-主程序开始，至少要输入五个参数，第一个是X坐标系列，然后是Y坐标系列，第三个是坐标点数，第四个是待查数值，第五个是返回值第六个参数是从哪一段曲线开始查找，如果曲线可以返回多个值，那么分别指定起始节点就可以找出全部合要求的点第七个参数是待查数值的类型,x代表输入x坐标求对应y坐标，y则相反,t是直接输入贝塞尔插值多项式的参数-Dim j As LongDim x1Value, y1Value, x2Value, y2Value, x3Value, y3Value As VariantDim ErrorMsg As VariantValueType = LCase(known_value_type) 待查数值的类型转化为小写，并赋值到全局变量ValueTypekey_value = known_value 待查数值赋值到全局变量key_valueErrorMsg = ErrorCheck(known_x, known_y, StartKnot) 检查输入错误If ErrorMsg NoError Then 有错误就返回错误信息，退出程序 result = Array(ErrorMsg, ErrorMsg, ErrorMsg, ErrorMsg, ErrorMsg, ErrorMsg) Exit SubEnd IfSizeX = UBound(known_x)For j = StartKnot To size SizeX - 1 从指定的节点开始，没有指定节点就从1开始 Call FindFourDots(known_x, known_y, j) 找出输入X-Y点坐标里面，应该用于计算的四个结点 Call FindFourBezierPoints(Dot1, Dot2, Dot3, Dot4) 根据四个结点计算四个贝塞尔控制点 Call FindABCD 根据待查数值的类型，和贝塞尔控制点,计算贝塞尔插值多项式的系数 Call Find_t 检查贝塞尔曲线是否包含待查数值 If Interpol_here = True Then Exit ForNext jIf Interpol_here = True Then 计算点坐标，并返回 以下是由四个贝塞尔控制点决定的,贝塞尔曲线的参数方程 x1Value = (1 - t1) 3 * BezierPt1.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 + 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 * 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 SubFunction 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坐标数目不等于Y坐标数目ErrorCheck = Error1Exit FunctionEnd IfIf SizeX 3 Then 输入的X-Y坐标对少于三个ErrorCheck = Error2Exit FunctionEnd IfIf (StartKnot = SizeX) Then 指定的起始节点超出范围ErrorCheck = Error3Exit FunctionEnd IfIf (ValueType x And ValueType y And ValueType t) Then 输入的待查数值类型不是x, y, tErrorCheck = Error4Exit FunctionEnd IfIf (ValueType = t And key_value 1) Or (ValueType = t And key_value 0) Then t 类型的范围是0-1ErrorCheck = Error5Exit FunctionEnd IfEnd FunctionSub FindFourDots(ByRef known_x() As Double, ByRef known_y() As Double, j) 根据X-Y数值，及起始节点，找出用于计算的四个结点坐标 If j = 1 Then 第一个结点 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 = known_x(j) Dot2.y = known_y(j) Dot3.x = known_x(j + 1) Dot3.y = known_y(j + 1) If j = SizeX - 1 Then 最后一个结点 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 SubSub FindFourBezierPoints(Dot1 As Vector, Dot2 As Vector, Dot3 As Vector, Dot4 As Vector)Dim d12, d23, d34, d13, d14, d24 As Doubled12 = DistAtoB(Dot1, Dot2) 计算平面坐标系上的两点距离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(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) 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) OffsetTo3 = 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 SubFunction DistAtoB(dota As Vector, dotb 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 = dota.y - dotb.yEnd FunctionFunction AmultiF(dota As Vector, MultiFactor As Double) As VectorAmultiF.x = dota.x * MultiFactorAmultiF.y = dota.y * MultiFactorEnd FunctionSub FindABCD()If ValueType = x Then 参数类型是x, 需要解参数方程 f(t) = x，这里设定参数方程的系数a = -BezierPt1.x + 3 * BezierPt2.x - 3 * BezierPt3.x + BezierPt4.xb = 3 * BezierPt1.x - 6 * BezierPt2.x + 3 * BezierPt3.xc = -3 * BezierPt1.x + 3 * BezierPt2.xd = BezierPt1.x - key_valueEnd IfIf ValueType = y Then 参数类型是x, 需要解参数方程 f(t) = y，这里设定参数方程的系数a = -BezierPt1.y + 3 * BezierPt2.y - 3 * BezierPt3.y + BezierPt4.yb = 3 * BezierPt1.y - 6 * BezierPt2.y + 3 * BezierPt3.yc = -3 * BezierPt1.y + 3 * BezierPt2.yd = BezierPt1.y - key_valueEnd IfEnd SubSub Find_t() 计算当 f(t) = 待查数值时, t应该是什么数值Dim tArr As VariantInterpol_here = TrueIf ValueType = t Then 待查数值类型为t,那么无需计算 t1 = key_value t2 = key_value t3 = key_value Exit SubEnd IftArr = Solve_Order3_Equation(a, b, c, d) 否则，解三次贝塞尔参数方程 f(t) = 待查数值t1 = tArr(1) 解得方程的三个根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 三个根都不合要求，代表曲线上没有包含待查数值的点If (t1 = NoRoot And t2 NoRoot) Then 至少有一个根，则用它代替NoRoot的结果,方便Excel画图 t1 = t2End IfIf (t1 = NoRoot And t3 NoRoot) Then t1 = t3End IfIf (t2 = NoRoot) Then t2 = t1If (t3 = NoRoot) Then t3 = t1End Sub. 牛顿法解三次方程，先求解方程的导函数，得到方程的拐点(导数等于0的点) 然后分三段用迭代法分别求三个根.Public Function Solve_Order3_Equation(p3, p2, p1, P0, Optional Starting As Double = -10000000000#, Optional Ending As Double = 10000000000#) 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 如果三次方程没有三次项 Two_X = Solve_Order2_Equation(a2, a1, a0) 解释法直接求二次方程的解 x1 = Two_X(1) x2 = 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 = (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 En