贝塞尔曲线及插值全解

1、贝塞尔曲线及插值这里主要讲一下如何在excel及vb中实现贝塞尔曲线插值，程序来源于互联网（程序作者：海底眼（Mr.DragonPan在excel中用宏实现），本文作为少量修改，方便在vb中调用，经运行证明是没错的，下面程序可作成一个模块放到vb或vba中调用：Excel的平滑线散点图，可以根据两组分别代表X-Y坐标的散点数值产生曲线图但是，却没有提供这个曲线图的公式，所以无法查找曲线上的点坐标后来我在以下这个网页找到了详细的说明和示例程序 HYPERLINK /Smooth_curve_bezier_example_file.zip /Smooth_curve_bezierexample_f

2、ile.zip根据其中采用的算法，进一步增添根据X坐标求Y坐标，或根据Y坐标求X坐标，更切合实际需求这个自定义函数按照Excel的曲线算法（三次贝塞尔分段插值），计算平滑曲线上任意一点的点坐标Excel的平滑曲线的大致算法是:给出了两组X-Y数值以后，每一对X-Y坐标称为节点，然后在每两个节点之间画出三次贝塞尔曲线（下面简称曲线）贝塞尔曲线的算法网上有很多资源，这里不介绍了，只作简单说明每条曲线都由四个节点开始，计算出四个贝塞尔控制点，然后根据控制点画出唯一一条曲线假设曲线的源数据是节点1,节点2,节点3,节点4（Dotl,Dot2,Dot3,Dot4）那么贝塞尔控制点的计算如下Dot2是第一

3、个控制点,也是曲点的起点，Dot3是第四个控制点也是曲线的终点k、八、第二个控制点的位置是:过第一个控制点（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,

4、终点）与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|=1and=0and=1ConstError10=Error:known_valueisnotonthecurve(definedbygivenknown_xandknown_y)ConstNoRoot=NoRootConstMaxErr=0.00000001ConstM

5、axLoop=1000DimSizeX,SizeYAsLongDimDot1AsVector塞尔控制点的四个节点DimDot2AsVectorDimDot3AsVectorDimDot4AsVectorDimBezierPt1AsVector制点DimBezierPt2AsVectorDimBezierPt3AsVectorDimBezierPt4AsVectorDimOffsetTo2AsVector终点的距离关系DimOffsetTo3AsVector输入区域的大小输入区域里面，用作计算贝生成贝塞尔曲线的四个贝塞尔控第二,三个贝塞尔控制点跟起点，DimValueTypeAsVariant输

6、入待查数值的类型、”X代表输入的是X坐标，求对应的Y坐标DimInterpol_hereAsBoolean当前分段曲线是否包含待查数值Dimkey_value,a,b,c,dAsDouble贝塞尔曲线插值多项式的系数Dimt1,t2,t3AsVariant贝塞尔曲线插值多项式的根Dima3,a2,a1,a0AsDoubleDimsize%PublicSubbefit(ByRefknown_x()AsDouble,ByRefknown_y()AsDouble,sizeAsInteger,known_valueAsDouble,result()AsVariant,OptionalStartKnot

7、AsLong=1,Optionalknown_value_typeAsVariant=”x”)子过程方便VB中调用主程序开始，至少要输入五个参数，第一个是X坐标系列，然后是Y坐标系列,第三个是坐标点数，第四个是待查数值，第五个是返回值第六个参数是从哪一段曲线开始查找，如果曲线可以返回多个值，那么分别指定起始节点就可以找出全部合要求的点第七个参数是待查数值的类型，”x代表输入x坐标求对应y坐标，”y”则相反,t”是直接输入贝塞尔插值多项式的参数DimjAsLongDimx1Value，y1Value，x2Value，y2Value，x3Value，y3ValueAsVariantDimError

8、MsgAsVariant待查数值的类型转化为小写，ValueType=LCase(known_value_type)并赋值到全局变量ValueTypekey_value=known_value待查数值赋值到全局变量key_valueErrorMsg=ErrorCheck(known_x,known_y,StartKnot)检查输入错误IfErrorMsgNoErrorThen有错误就返回错误信息，退出程序result=Array(ErrorMsg,ErrorMsg,ErrorMsg,ErrorMsg,ErrorMsg,ErrorMsg)ExitSubEndIfSizeX=UBound(know

9、n_x)Forj=StartKnotTosizeSizeX-1从指定的节点开始，没有指定节点就从1开始CallFindFourDots(known_x,known_y,j)找出输入X-Y点坐标里面，应该用于计算的四个结点CallFindFourBezierPoints(Dot1,Dot2,Dot3,Dot4)根据四个结点计算四个贝塞尔控制点CallFindABCD根据待查数值的类型，和贝塞尔控制点,计算贝塞尔插值多项式的系数CallFind_t检查贝塞尔曲线是否包含待查数值IfInterpol_here=TrueThenExitForNextjIfInterpol_here=TrueThen计

10、算点坐标，并返回以下是由四个贝塞尔控制点决定的,贝塞尔曲线的参数方程xlValue=(1-tl)八3*BezierPtl.x+3*tl*(1-tl)八2*BezierPt2.x+3*tl八2*(1-tl)*BezierPt3.x+tl八3*BezierPt4.xylValue=(l-tl)八3*BezierPtl.y+3*tl*(l-tl)八2*BezierPt2.y+3*tl八2*(l-tl)*BezierPt3.y+tl八3*BezierPt4.yx2Value=(l-t2)八3*BezierPtl.x+3*t2*(l-t2)八2*BezierPt2.x+3*t2八2*(l-t2)*Bez

11、ierPt3.x+t2八3*BezierPt4.xy2Value=(l-t2)八3*BezierPtl.y+3*t2*(l-t2)八2*BezierPt2.y+3*t2八2*(l-t2)*BezierPt3.y+t2八3*BezierPt4.yx3Value=(l-t3)八3*BezierPtl.x+3*t3*(l-t3)八2*BezierPt2.x+3*t3八2*(l-t3)*BezierPt3.x+t3八3*BezierPt4.xy3Value=(l-t3)八3*BezierPtl.y+3*t3*(l-t3)八2*BezierPt2.y+3*t3八2*(l-t3)*BezierPt3.y+

12、t3八3*BezierPt4.yresult=Array(x1Value,y1Value,x2Value,y2Value,x3Value,y3Value)Elseresult=Array(Error10,Error10,Error10,Error10,Error10,Error10)EndIfEndSub/*FunctionErrorCheck(ByRefknown_x()AsDouble,ByRefknown_y()AsDouble,StartKnot)AsVariantErrorCheck=NoErrorSizeX=UBound(known_x)known_x.CountSizeY=UBo

13、und(known_y)known_y.CountIfSizeXSizeYThen假如输入的X坐标数目不等于Y坐标数目ErrorCheck=Error1ExitFunctionEndIfIfSizeX3Then输入的X-Y坐标对少于三个ErrorCheck=Error2ExitFunctionEndIfIf(StartKnot=SizeX)Then指定的起始节点超出范围ErrorCheck=Error3ExitFunctionEndIfIf(ValueTypexAndValueTypeyAndValueTypet)Then输入的待查数值类型不是x,y,tErrorCheck=Error4Exi

14、tFunctionEndIfIf(ValueType=tAndkey_value1)Or(ValueType=tAndkey_value0)Thent类型的范围是0-1ErrorCheck=Error5ExitFunctionEndIfEndFunction/*SubFindFourDots(ByRefknown_x()AsDouble,ByRefknown_y()AsDouble,j)根据X-Y数值，及起始节点，找出用于计算的四个结点坐标Ifj=1Then第一个结点Dot2=Dot1Dot1.x=known_x(1)Dot1.y=known_y(1)ElseDot1.x=known_x(j-

15、1)Dot1.y=known_y(j-1)EndIfDot2.x=known_x(j)Dot2.y=known_y(j)Dot3.x=known_x(j+1)Dot3.y=known_y(j+1)Ifj=SizeX-1Then最后一个结点Dot4=Dot3Dot4.x=Dot3.xDot4.y=Dot3.yElseDot4.x=known_x(j+2)Dot4.y=known_y(j+2)EndIfEndSub/*SubFindFourBezierPoints(Dot1AsVector,Dot2AsVector,Dot3AsVector,Dot4AsVector)Dimd12,d23,d34,

16、d13,d14,d24AsDoubled12=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/6d23/2)And(d24/6d23/2)ThenIf(Dot1.xDot2

17、.xOrDot1.yDot2.y)ThenOffsetTo2=AmultiF(OffsetTo2,1/6)If(Dot1.x=Dot2.xAndDot1.y=Dot2.y)ThenOffsetTo2=AmultiF(OffsetTo2,1/3)If(Dot3.xDot4.xOrDot3.yDot4.y)ThenOffsetTo3=AmultiF(OffsetTo3,1/6)If(Dot3.x=Dot4.xAndDot3.y=Dot4.y)ThenOffsetTo3=AmultiF(OffsetTo3,1/3)ElseIf(d13/6=d23/2)And(d24/6=d23/2)ThenOffs

18、etTo2=AmultiF(OffsetTo2,d23/12)OffsetTo3=AmultiF(OffsetTo3,d23/12)ElseIf(d13/6=d23/2)ThenOffsetTo2=AmultiF(OffsetTo2,d23/2/d13)OffsetTo3=AmultiF(OffsetTo3,d23/2/d13)ElseIf(d24/6=d23/2)ThenOffsetTo2=AmultiF(OffsetTo2,d23/2/d24)OffsetTo3=AmultiF(OffsetTo3,d23/2/d24)EndIfBezierPt2=AaddB(BezierPt1,Offse

19、tTo2)向量加法BezierPt3=AaddB(BezierPt4,OffsetTo3)EndSub/*FunctionDistAtoB(dotaAsVector,dotbAsVector)AsDoubleDistAtoB=(dota.xdotb.x)八2+(dota.ydotb.y)八2)八0.5EndFunctionFunctionAaddB(dotaAsVector,dotbAsVector)AsVectorAaddB.x=dota.x+dotb.xAaddB.y=dota.y+dotb.yEndFunctionFunctionAsubB(dotaAsVector,dotbAsVect

20、or)AsVectorAsubB.x=dota.xdotb.xAsubB.y=dota.ydotb.yEndFunctionFunctionAmultiF(dotaAsVector,MultiFactorAsDouble)AsVectorAmultiF.x=dota.x*MultiFactorAmultiF.y=dota.y*MultiFactorEndFunction/*SubFindABCD()IfValueType=xThen参数类型是x,需要解参数方程f(t)=x，这里设定参数方程的系数a=-BezierPt1.x+3*BezierPt2.x-3*BezierPt3.x+BezierP

21、t4.xb=3*BezierPt1.x-6*BezierPt2.x+3*BezierPt3.xc=-3*BezierPt1.x+3*BezierPt2.xd=BezierPt1.x-key_valueEndIfIfValueType=yThen参数类型是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_v

22、alueEndIfEndSub/*SubFind_t()计算当f(t)=待查数值时,t应该是什么数值DimtArrAsVariantInterpol_here=TrueIfValueType=tThen待查数值类型为t,那么无需计算t1=key_valuet2=key_valuet3=key_valueExitSubEndIf否则，解三次贝塞尔参数tArr=Solve_Order3_Equation(a,b,c,d)方程f(t)=待查数值t1=tArr(1)解得方程的三个根t2=tArr(2)t3=tArr(3)If(t11Ort11Ort21Ort30)Thent3=NoRootEndIfI

23、f(IsNumeric(t1)=FalseAndIsNumeric(t2)=FalseAndIsNumeric(t3)=False)ThenInterpol_here=FalseEndIf三个根都不合要求，代表曲线上没有包含待查数值的点If(t1=NoRootAndt2NoRoot)Then至少有一个根，则用它代替NoRoot的结果,方便Excel画图t1=t2EndIfIf(t1=NoRootAndt3NoRoot)Thent1=t3EndIfIf(t2=NoRoot)Thent2=t1If(t3=NoRoot)Thent3=t1EndSub/*牛顿法解三次方程，先求解方程的导函数，得到方程

24、的拐点（导数等于0的点）然后分三段用迭代法分别求三个根PublicFunctionSolve_Order3_Equation(p3,p2,p1,P0,OptionalStartingAsDouble=-10000000000#,OptionalEndingAsDouble=10000000000#)AsVariantDimTwo_X,TurningPoint,x1,x2,x3AsVariantDimxAsDoublea3=p3a2=p2a1=p1a0=P0 x1=NoRootx2=NoRootx3=NoRootx1=Newton_Solve(Starting)Ifa3=0Then如果三次方程没

25、有三次项Two_X=Solve_Order2_Equation(a2,a1,a0)解释法直接求二次方程的解x1=Two_X(1)x2=Two_X(2)ElseTurningPoint=Solve_Order2_Equation(3*a3,2*a2,1*a1)求解f(t)=0If(TurningPoint(1)=NoRootAndTurningPoint(2)=NoRoot)Then分段求根x=0 x1=Newton_Solve(x)ElseIf(TurningPoint(1)NoRootAndTurningPoint(2)=NoRoot)ThenIff_x(Starting)*f_x(Turn

26、ingPoint(1)0Thenx=(Starting+TurningPoint(1)/2x1=Newton_Solve(x)EndIfIff_x(TurningPoint(2)*f_x(Ending)0Thenx=(TurningPoint(2)+Ending)/2x3=Newton_Solve(x)EndIfElseIf(TurningPoint(1)NoRootAndTurningPoint(2)NoRoot)ThenIff_x(Starting)*f_x(TurningPoint(1)0Thenx=(Starting+TurningPoint(1)/2x1=Newton_Solve(x)EndIfIff_x(TurningPoint(1)*f_x(TurningPoint(2)0Thenx=(TurningPoint(1)+TurningPoint(2)/2x2=Newton_Solve(x)EndIfIff_x(TurningPoint(2)*f_x(Ending)0Thenx=(TurningPoint(2)+Ending)/2x3=New