贝塞尔曲线及插值

贝塞尔曲线及插值这里主要讲一下如何在excel及vb中实现贝塞尔曲线插值，程序来源于互联网（程序作者:海底眼(Mr.DragonPan在excel中用宏实现)，本文作为少量修改，方便在vb中调用，经运行证明是没错的，下面程序可作成一个模块放到vb或vba中调用：'Excel的平滑线散点图，可以根据两组分别代表X-Y坐标的散点数值产生曲线图'但是，却没有提供这个曲线图的公式，所以无法查找曲线上的点坐标'后来我在以下这个网页找到了详细的说明和示例程序'..............................................................................'根据其中采用的算法，进一步增添根据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/6*|Dot2_Dot4|'那么第二个控制点与第一控制点距离取1/2*|Dot2_Dot4|,而不是1/6*|Dot2_Dot4|'...............................................................................................'这个自定义函数的计算流程是'Step1:检查输入的X-Y数值是否有错误，如(输入不够三个点，X-Y的数量不一致,起始搜索节点超过范围等等)'Step2:从参数指定的节点开始，计算出四个贝塞尔控制点，得到贝塞尔插值多项式方程，'然后代入已知的待求数值，看它能不能满足f(t)=0有解(即曲线包含待查数值)'Step3:如果f(t)=0有解，根据解出来的t值计算X-Y坐标，退出程序,否则继续检查下一段曲线'Step4:如果所有分段曲线都不包含待查数值，退出程序'...............................................................................................OptionExplicitOptionBase1'所有数组的第一个元素编号为1(默认为0)TypeVector'自定义数据结构(用二维向量代表坐标系里面的点坐标)xAsDoubleyAsDoubleEndTypeConstNoError=Noerror'错误提示信息ConstError1=Error:Thesizeofknown_xmustequaltosizeofknown_yConstError2=Error:Thesizeofknown_xmustequaltoorgreaterthan3ConstError3=Error:StartKnotmustbe=1and=count(known_x)-1ConstError4=Error:known_value_typemustbex,y,ortConstError5=Error:Whenknown_value_typeist,known_valuemust=0and=1ConstError10=Error:known_valueisnotonthecurve(definedbygivenknown_xandknown_y)ConstNoRoot=NoRootConstMaxErr=0.00000001ConstMaxLoop=1000DimSizeX,SizeYAsLong'输入区域的大小DimDot1AsVector'输入区域里面，用作计算贝塞尔控制点的四个节点DimDot2AsVectorDimDot3AsVectorDimDot4AsVectorDimBezierPt1AsVector'生成贝塞尔曲线的四个贝塞尔控制点DimBezierPt2AsVectorDimBezierPt3AsVectorDimBezierPt4AsVectorDimOffsetTo2AsVector'第二,三个贝塞尔控制点跟起点，终点的距离关系DimOffsetTo3AsVectorDimValueTypeAsVariant'输入待查数值的类型,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,OptionalStartKnotAsLong=1,Optionalknown_value_typeAsVariant=x)''--------------------------------------子过程方便VB中调用-----------------------------------------------------------'主程序开始，至少要输入五个参数，第一个是X坐标系列，然后是Y坐标系列，第三个是坐标点数，第四个是待查数值，第五个是返回值'第六个参数是从哪一段曲线开始查找，如果曲线可以返回多个值，那么分别指定起始节点就可以找出全部合要求的点'第七个参数是待查数值的类型,x代表输入x坐标求对应y坐标，y则相反,t是直接输入贝塞尔插值多项式的参数'-------------------------------------------------------------------------------------------------DimjAsLongDimx1Value,y1Value,x2Value,y2Value,x3Value,y3ValueAsVariantDimErrorMsgAsVariantValueType=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)ExitSubEndIf'SizeX=UBound(known_x)Forj=StartKnotTosize'SizeX-1'从指定的节点开始，没有指定节点就从1开始CallFindFourDots(known_x,known_y,j)'找出输入X-Y点坐标里面，应该用于计算的四个结点CallFindFourBezierPoints(Dot1,Dot2,Dot3,Dot4)'根据四个结点计算四个贝塞尔控制点CallFindABCD'根据待查数值的类型，和贝塞尔控制点,计算贝塞尔插值多项式的系数CallFind_t'检查贝塞尔曲线是否包含待查数值IfInterpol_here=TrueThenExitForNextjIfInterpol_here=TrueThen'计算点坐标，并返回'以下是由四个贝塞尔控制点决定的,贝塞尔曲线的参数方程x1Value=(1-t1)^3*BezierPt1.x+3*t1*(1-t1)^2*BezierPt2.x+3*t1^2*(1-t1)*BezierPt3.x+t1^3*BezierPt4.xy1Value=(1-t1)^3*BezierPt1.y+3*t1*(1-t1)^2*BezierPt2.y+3*t1^2*(1-t1)*BezierPt3.y+t1^3*BezierPt4.yx2Value=(1-t2)^3*BezierPt1.x+3*t2*(1-t2)^2*BezierPt2.x+3*t2^2*(1-t2)*BezierPt3.x+t2^3*BezierPt4.xy2Value=(1-t2)^3*BezierPt1.y+3*t2*(1-t2)^2*BezierPt2.y+3*t2^2*(1-t2)*BezierPt3.y+t2^3*BezierPt4.yx3Value=(1-t3)^3*BezierPt1.x+3*t3*(1-t3)^2*BezierPt2.x+3*t3^2*(1-t3)*BezierPt3.x+t3^3*BezierPt4.xy3Value=(1-t3)^3*BezierPt1.y+3*t3*(1-t3)^2*BezierPt2.y+3*t3^2*(1-t3)*BezierPt3.y+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=UBound(known_y)'known_y.CountIfSizeXSizeY