Nurbs样条线算法推导及python实现

Nurbs样条线算法推导及python实现注：在现在的建模软件中Nurbs曲线还引入一个跨度的概念，在这里并不体现，这里是有几个控制点就是控制点数-1阶，最后好像到36阶曲线就发散了耶三个控制点：a=(ax,ay)b=(bx,by)c=(cx,cy)0≤t≤1e为ab线上一点则：ae⃑⃑⃑⃑=t∗ab⃑⃑⃑⃑e=(ex,ey)=(t*(bx-ax)+ax,t*(by-ay)+ay)同理，同一t值，f为bc线上一点为：f=(fx,fy)=(t*(cx-bx)+bx,t*(cy-by)+by)g点为动线ef上的动点为：g=(t*(fx-ex)+ex,t*(fy-ey)+ey)=(gx,gy)gx=t^2*(cx-2*bx+ax)+2*t*(bx-ax)+axgy=t^2*(cy-2*by+ay)+2*t*(by-ay)+ay动点g的运动方程入上图所示0:(x0,y0);1:(x1,y1);2:(x3,y3);x=t2∙(x2−2∙x1+x0)+t∙2∙(x1−x0)+x0;y=t2∙(y2−2∙y1+y0)+t∙2∙(y1−y0)+y0;abc四个控制点0≤t≤1a=(ax,ay)(x0,y0)b=(bx,by)(x1,y1)c=(cx,cy)(x2,y2)d=(dx,dy)(x3,y3)e=(ex,ey)=(t*(bx-ax)+ax,t*(by-ay)+ay)t∙(x1−x0)+x0,…f=(fx,fy)=(t*(cx-bx)+bx,t*(cy-by)+by)t∙(x2−x1)+x1,…g=(gx,gy)=(t*(dx-cx)+cx,t*(dy-cy)+cy)t∙(x3−x2)+x2,…h=(hx,hy)=(t*(fx-ex)+ex,t*(fy-ey)+ey)t2∙(x2−2∙x1+x0)+2∙t∙(x1−x0)+x0,…I=(ix,iy)=(t*(gx-fx)+fx,t*(gy-fy)+fy)t2∙(x3−2∙x2+x1)+2∙t∙(x2−x1)+x1,…j=(jx,jy)=(t*(ix-hx)+hx,t*(iy-hy)+hy)abcdt3∙(x3−3∙x2+3∙x1−x0)+3∙t2∙(x2−2∙x1+x0)+3∙t∙(x1−x0)+x0,…t3∙(x3−3∙(x2−x1)−x0)+3∙t2∙(x2−2∙x1+x0)+3∙t∙(x1−x0)+x0,…N个控制点：defcv(N,t):n=[Point(t*(N[i+1].x-N[i].x)+N[i].x,t*(N[i+1].y-N[i].y)+N[i].yforiinrange(len(N)-1))]iflen(n)==1:returnn[0]else:returncv(n,t)0≤t≤1N=[p0,p1,p2…pn]printcv(N,t)python代码实现：（按分段数计算样条线上的点）defcv(N,seg):step=1/segsP=N[0]eP=N[-1]p=[sP]oldN=Nforiinrange(step,1,step):t=i*stepcount=len(N)whilecount1:n=[Point(t*(N[i+1].x-N[i].x)+N[i].x,t*(N[i+1].y-N[i].y)+N[i].y)foriinrange(count-1)]N=ncount-=1N=oldNp.append(n[0])p.append(eP)returnp因为y和x的运算过程相同的故下面只求x用多项式的方法优化运算：各项常数确定：Cx,Cy(x0,y0)(x1,y1)(x2,y2)(x3,y3)(x4,y4)⇊t∙(x1−x0)+x0,…t∙(x2−x1)+x1,…t∙(x3−x2)+x2,…t∙(x4−x3)+x3,…⇊t2∙((x2−x1)−(x1−x0))+2∙t∙(x1−x0)+x0t2∙((x3−x2)−(x2−x1))+2∙t∙(x2−x1)+x1t2∙((x4−x3)−(x3−x2))+2∙t∙(x3−x2)+x2⇊t3∙(((x3−x2)−(x2−x1))−((x2−x1)−(x1−x0)))+3∙t2∙((x2−x1)−(x1−x0))+3∙t∙(x1−x0)+x0t3∙(((x4−x3)−(x3−x2))−((x3−x2)−(x2−x1)))+3∙t2∙((x3−x2)−(x2−x1))+3∙t∙(x2−x1)+x1⇊t4∙((((x4−x3)−(x3−x2))−((x3−x2)−(x2−x1)))−(((x3−x2)−(x2−x1))−((x2−x1)−(x1−x0))))+4∙t3∙(((x3−x2)−(x2−x1))−((x2−x1)−(x1−x0)))+6∙t2∙((x2−x1)−(x1−x0))+4∙t∙(x1−x0)+x0;01234C0=x04∙C1=4∙(x1−C0)6∙C2=6∙((x2−x1)−C1)4∙C3=4∙(((x3−x2)−(x2−x1))−C2)C4=((((x4−x3)−(x3−x2))−((x3−x2)−(x2−x1)))−C3)⇊D1=x1D2=(x2−D1)D3=((x3−x2)−D2)D4=(((x4−x3)−(x3−x2))−D3)⇊E2=x2E3=(x3−E2)E4=((x4−x3)−E3)⇊F3=x3F4=(x4−F3)⇊G4=x4⇊F4=(G4−F3)⟺G4=F3+F4E4=(F4−E3)⟺F4=E3+E4E3=(F3−E2)⟺F3=E2+E3D4=(E4−D3)⟺E4=D3+D4D3=(E3−D2)⟺E3=D2+D3D2=(E2−D1)⟺E2=D1+D2C4=(D4−C3)⟺D4=C3+C4C3=(D3−C2)⟺D3=C2+C3C2=(D2−C1)⟺D2=C1+C2C1=(D1−C0)⟺D1=C0+C1G4F3,F4E2,E3,E4D1,D2,D3,D4C0,C1,C2,C3,C4⇕x4x3,F4x2,E3,E4x1,D2,D3,D4x0,C1,C2,C3,C4各项系数增量：scale1,11,2,11,3,3,11,4,6,4,1各项系数增量scale:scalen,i=scalen−1,i−1+scalen−1,i；(scalen,0=1,scalen,n=1)scalei=scalen,idefgetScales(n):'''各项常系数的缩放值1,1n=11,2,1n=21,3,3,1n=31,4,6,4,1n=4......返回第n行数据n为控制点数目'''ifn2:returnNonen=n-2s1=[1,1]scales=[1,1]foriinrange(n):scales=[s1[j]+s1[j+1]forjinrange(len(s1)-1)]scales.insert(0,1)scales.append(1)s1=scalesreturnscales各项常系数CCn,i=Cn−1,i−1−Cn,i−1；(Cn,0=x0或y0或z0)Ci=scalei∙Cn,idefgetCs(scales,points):'''各项常系数x_4n=0x_3,F_4n=11234x_2,E_3,E_4n=2x_1,D_2,D_3,D_4n=3x_0,C_1,C_2,C_3,C_4n=5......points为控制点（三角最左边对应各控制点同一轴向坐标，顺序自下向上对应控制点输入顺序）返回各轴轴常系数-第n行数据（杨辉三角底行）'''points=list(points)points.reverse()axises=zip(*points)Cs=[]foraxisinaxises:C1=[axis[0]]foriinrange(1,len(axis)):C=[]C.append(axis[i])forjinrange(len(C1)):C.append(C1[j]-C[j])C1=CC=tuple(map(lambdac,s:c*s,C,scales))Cs.append(C)returnCs各分段各项ti值：defgetTs(n,seg):'''n控制点数目，seg分段数按分段数获取曲线上所有分段点对应0=t=1范围内t**i的所有数据'''step=1/segTs=[]foriinrange(seg+1):t1=[]s=i*stepforjinrange(n):t1.append(s**j)Ts.append(t1)Ts.append([1]*n)returnTs样条线多项式：X=∑ti∙Cin−1i=0；(n为控制点数目)defcv(Ts,Cs):'''返回样条线上所有分段点坐标'''n=len(Cs[0])nAxis=len(Cs)ps=[]foriinrange(len(Ts)):sumR=[0foriinrange(nAxis)]forjinrange(n):forkinrange(nAxis):sumR[k]+=Ts[i][j]*Cs[k][j]ps.append(sumR)returnps样条线长度：未完待续…判断一点是否在样条线内：未完待续..Python代码：#nurbs#-*-coding:utf-8-*-#fromcopyimportdeepcopydefgetScales(self,n):'''各项常系数的缩放值1,1n=11,2,1n=21,3,3,1n=31,4,6,4,1n=4......返回第n行数据n为控制点数目'''ifn2:returnNonen=n-2s1=[1,1]scales=[1,1]foriinrange(n):scales=[s1[j]+s1[j+1]forjinrange(len(s1)-1)]scales.insert(0,1)scales.append(1)s1=scalesreturnscalesdefgetCs(self,scales,points):'''各项常系数x_4n=0x_3,F_4n=1x_2,E_3,E_4n=2x_1,D_2,D_3,D_4n=3x_0,C_1,C_2,C_3,C_4n=5......points为控制点（三角最左边对应各控制点同一轴向坐标，顺序自下向上对应控制点输入顺序）返回各轴轴常系数-第n行数据（杨辉三角底行）'''ifscalesisNone:returnNone#points=list(points)#points.reverse()#axises=zip(*points)axises=zip(*reversed(points))Cs=[]foraxisinaxises:C1=[axis[0]]foriinrange(1,len(axis)):C=[]C.append(axis[i])forjinrange(len(C1)):C.append(C1[j]-C[j])C1=C#print(t)C=tuple(map(lambdac,s:c*s,C,scales))Cs.append(C)returnCsdefgetTs(self,n,seg):'''n控制点数目，seg分段数按分段数获取曲线上所有分段点对应0=t=1范围内t**i的所有数据'''ifn1:returnNonestep=1/segTs=[]foriinrange(seg+1):s=i*step#t1=[]#forjinrange(n):#t1.append(s**j)t1=[s**jforjinrange(n)]Ts.append(t1)Ts.append([1]*n)returnTsdefgetCvps(self,Ts,Cs):'''返回样条线上所有分段点坐标'''ifCsisNone:returnNonen=len(Cs[0])nAxis=len(Cs)ps=[]foriinrange(len(Ts)-1):sumR=[0foriinrange(nAxis)]forjinrange(n):forkinrange(nAxis):sumR[k]+=Ts[i][j]*Cs[k][j]ps.append