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

Nurbs样条线算法推导及python实现注：在现在的建模软件中Nurbs曲线还引入一个跨度的概念，在这里并不体现，这里是有几个控制点就是控制点数-1阶，最后好像到36阶曲线就发散了耶三个控制点：cbaa=(ax, ay)b=(bx, by)c=(cx, cy)0t1e为ab线上一点则：ae=t*abe=(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= t2*(cx-2*bx+ax)+2*t*(bx-ax)+axgy= t2*(cy-2*by+ay)+2*t*(by-ay)+ay动点g的运动方程入上图所示0:x0, y0;1:x1, y1;2:x3, y3;x=t2x2-2x1+x0+t2x1-x0+x0;y=t2y2-2y1+y0+t2y1-y0+y0;四个控制点cdba0t1a=(ax, ay) x0, y0b=(bx, by) x1, y1c=(cx, cy) x2, y2d=(dx, dy) x3, y3e=(ex, ey)=(t*(bx-ax)+ax, t*(by-ay)+ay) tx1-x0+x0,f = (fx , fy)=(t*(cx-bx)+bx, t*(cy-by)+by)tx2-x1+x1,g= (gx, gy)=(t*(dx-cx)+cx, t*(dy-cy)+cy)tx3-x2+x2,h=(hx, hy)=(t*(fx-ex)+ex, t*(fy-ey)+ey)t2x2-2x1+x0+2tx1-x0+x0,I = (ix , iy )=(t*(gx-fx)+fx, t*(gy-fy)+fy)t2x3-2x2+x1+2tx2-x1+x1,j = (jx, jy )=(t*(ix-hx)+hx, t*(iy-hy)+hy)t3x3-3x2+3x1-x0+3t2x2-2x1+x0+3tx1-x0+x0,t3x3-3x2-x1-x0+3t2x2-2x1+x0+3tx1-x0+x0,N个控制点：def cv(N, t):n= Point(t*( Ni+1.x- Ni.x)+Ni.x, t*( Ni+1.y- Ni.y)+Ni.yfor i in range(len(N)-1) )if len(n)=1: return n0else: return cv(n, t)0t1N=p0, p1, p2pnprint cv(N, t)python代码实现：（按分段数计算样条线上的点）def cv(N, seg):step=1/segsP=N0eP=N-1p=sPoldN=Nfor i in range(step, 1, step): t=i*stepcount= len(N)while count1:n= Point(t*( Ni+1.x- Ni.x)+Ni.x , t*( Ni+1.y- Ni.y)+Ni.y) for i in range(count-1) N=ncount-=1 N=oldN p.append(n0)p.append(eP)return p用多项式的方法优化运算：各项常数确定：Cx, Cy因为y和x的运算过程相同的故下面只求xx0, y0x1, y1x2, y2x3, y3x4, y4tx1-x0+x0,tx2-x1+x1,tx3-x2+x2,tx4-x3+x3,t2x2-x1-x1-x0+2tx1-x0+x0t2x3-x2-x2-x1+2tx2-x1+x1t2x4-x3-x3-x2+2tx3-x2+x2t3x3-x2-x2-x1- x2-x1-x1-x0+3t2x2-x1-x1-x0+3tx1-x0+x0t3x4-x3-x3-x2- x3-x2-x2-x1+3t2x3-x2-x2-x1+3tx2-x1+x1t4 x4-x3-x3-x2- x3-x2-x2-x1 - x3-x2-x2-x1- x2-x1-x1-x0+4t3x3-x2-x2-x1- x2-x1-x1-x0+6t2x2-x1-x1-x0+4tx1-x0+x0;C0=x04C1=4x1-C06C2=6x2-x1-C14C3=4x3-x2-x2-x1-C2C4= x4-x3-x3-x2- x3-x2-x2-x1-C3D1=x1D2=x2-D1D3=x3-x2-D2D4= x4-x3-x3-x2- D3E2=x2E3=x3-E2E4=x4-x3-E3F3=x3F4=x4-F3G4=x4F4=G4-F3G4=F3+F4E4=F4-E3F4=E3+E4E3=F3-E2F3=E2+E3D4=E4-D3E4=D3+D4D3=E3-D2E3=D2+D3D2=E2-D1E2=D1+D2C4= D4-C3D4=C3+C4C3= D3-C2D3=C2+C3C2= D2-C1D2=C1+C2C1= D1-C0D1=C0+C1G4F3,F4E2,E3, E4D1, D2, D3, D4C0, C1, C2, C3, C401234x4x3,F4x2,E3, E4x1, D2, D3, D4x0, C1, C2, C3, C41234各项系数增量：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=1scalei=scalen, idef getScales(n): 各项常系数的缩放值 1,1 n=1 1,2,1 n=2 1,3,3,1 n=3 1,4,6,4,1 n=4 . . 返回第n行数据 n为控制点数目 if n2: return None n=n-2 s1=1,1 scales=1,1 for i in range(n): scales=s1j+s1j+1 for j in range(len(s1)-1) scales.insert(0,1) scales.append(1) s1=scales return scales各项常系数C Cn, i=Cn-1, i-1-Cn, i-1；Cn, 0=x0或y0或z0C i=scaleiCn, idef getCs(scales,points): 各项常系数 x_4 n=0 x_3,F_4 n=1 x_2,E_3,E_4 n=2 x_1,D_2,D_3,D_4 n=3 x_0,C_1,C_2,C_3,C_4 n=5 . . points为控制点 （三角最左边对应各控制点同一轴向坐标， 顺序自下向上对应控制点输入顺序） 返回各轴轴常系数-第n行数据（杨辉三角底行） points=list(points) points.reverse() axises=zip(*points) Cs= for axis in axises: C1=axis0 for i in range(1,len(axis): C= C.append(axisi) for j in range(len(C1): C.append(C1j-Cj) C1=C C=tuple(map(lambda c,s:c*s, C, scales) Cs.append(C) return Cs各分段各项ti值：def getTs(n, seg): n控制点数目，seg分段数 按分段数获取曲线上所有分段点对应 0=t=1范围内t*i的所有数据 step=1/seg Ts= for i in range(seg+1): t1= s=i*step for j in range(n): t1.append(s*j) Ts.append(t1) Ts.append(1*n) return Ts样条线多项式：X=i=0n-1tiCi；(n为控制点数目)def cv(Ts, Cs): 返回样条线上所有分段点坐标 n=len(Cs0) nAxis=len(Cs) ps= for i in range(len(Ts): sumR=0 for i in range(nAxis) for j in range(n): for k in range(nAxis): sumRk+=Tsij*Cskj ps.append(sumR) return ps样条线长度：未完待续判断一点是否在样条线内：未完待续.Python代码：#nurbs# -*- coding: utf-8 -*-#from copy import deepcopydef getScales(self, n): 各项常系数的缩放值 1,1 n=1 1,2,1 n=2 1,3,3,1 n=3 1,4,6,4,1 n=4 . . 返回第n行数据 n为控制点数目 if n2: return None n=n-2 s1=1,1 scales=1,1 for i in range(n): scales=s1j+s1j+1 for j in range(len(s1)-1) scales.insert(0,1) scales.append(1) s1=scales return scalesdef getCs(self, scales, points): 各项常系数 x_4 n=0 x_3,F_4 n=1 x_2,E_3,E_4 n=2 x_1,D_2,D_3,D_4 n=3 x_0,C_1,C_2,C_3,C_4 n=5 . . points为控制点 （三角最左边对应各控制点同一轴向坐标， 顺序自下向上对应控制点输入顺序） 返回各轴轴常系数-第n行数据（杨辉三角底行） if scales is None: return None #points=list(points)#points.reverse()#axises=zip(*points) axises=zip(*reversed(points) Cs= for axis in axises: C1=axis0 for i in range(1,len(axis): C= C.append(axisi) for j in range(len(C1): C.append(C1j-Cj) C1=C #print (t) C=tuple(map(lambda c,s:c*s, C, scales) Cs.append(C) return Csdef getTs(self, n, seg): n控制点数目，seg分段数 按分段数获取曲线上所有分段点对应 0=t=1范围内t*i的所有数据 if n1: return None step=1/seg Ts= for i in range(seg+1): s=i*step #t1=#for j in range(n):#t1.append(s*j) t1=s*j for j in range(n) Ts.append(t1) Ts.append(1*n) return Ts def getCvps(self, Ts, Cs): 返回样条线上所有分段点坐标 if Cs is None: return None n=len(Cs0) nAxis=len(Cs) ps= for i in range(len(Ts)-1): sumR=0 for i in range(nAxis) for j in range(n): for k in range(nAxis): sumRk+=Tsij*Cskj ps.append(sumR) return ps class Nurbs(object): def _init_(self, points=(), segs=100): self._seg=segs self._points=points self._n=len(self._points) self._scales=getScales(self, self._n) self._Cs=getCs(self, self._scales, self._points) self._Ts=getTs(self, self._n, self._seg) self._cvps=getCvps(self, self._Ts, self.Cs) #get or set ctrPoints property def ctrPoints(self): return deepcopy(self._points) ctrPoints.setter def ctrPoints(self, points): if points!=self._points: self._points=deepcopy(points) n=len(self._points) if n!=self._n: self._n=n self._scales=getScales(self, self._n) self._Ts=getTs(self, self._n, self._seg) self._Cs=getCs(self, self._scales, self._points) self._cvps=getCvps(self, self._Ts, self.Cs) #get or set ctrPoints property def seg(self): return self._seg seg.setter def seg(self, segs): if segs!=self._seg: self._seg=segs self._Ts=getTs(self, self._n, self._seg) self._cvps=getCvps(self, self._Ts, self.Cs) property def scales(self): return deepcopy(self._scales) property def Cs(self): return deepcopy(self._Cs) property def Ts(self): return deepcopy(self._Ts) property def cvps(self): return deepcopy(self._cvps) #get ctrPoints amount def _len_(self): return self._n #tkinter交互程序#-*- coding:utf-8 -*-#import tkinter, time, nurbsclass Cv(tkinter.Canvas): def getObjs(self, tag, p, r): objs=self.find_overlapping(p0-r, p1-r, p0+r, p1+r) return tuple(filter(lambda obj:tag in self.gettags(obj), objs) def createPoint(self, pos, radius, color=None, lineColor=None): return self.create_oval(pos0-radius, pos1-radius, pos0+radius, pos1+radius, fill=color, outline=lineColor, tags=ctrp) def modifyPoint(self, id, pos=None, color=None, lineColor=None): if pos: x0, y0, x1, y1=self.coords(id) wm, hm=(x1-x0)/2, (y1-y0)/2 self.coords(id, pos0-wm,pos1-hm, pos0+wm, pos1+hm) if color: self.itemconfig(id, fill=color) if lineColor: self.itemconfig(id, outline=lineColor) def unpackCoords(self, l): return j for i in l for j in i def sortCtrPs(self, dic, l): return dici for i in l def _init_(self, root, size=(640, 480): super()._init_(root, bg=white, width=size0, height=size1) self.lineDict=#id:nb self.curLineId=None#cvId self.dotLineId=None#ctrpCvId self.curCtrPsDict=#ctrpId:pos self.curCtrPsInd=#ctrpId. self.curCtrPId=None#ctrpId self.need2update=False self.editState=new self.running=1 root.bind(, self.eQuit) self.bind(, self.eLeftClick) self.bind(, self.eLeftDrag) self.bind(, self.eLeftRelease) self.bind(, self.eMidClick) self.bind(, self.eRightClick) self.colorNormal=darkorange self.colorPick=green self.colorEnd=blue self.colorDashLine=orange self.colorCurve=darkblue self.stateText=self.create_text(100, 100, text=self.editState) self.pack(side=left) self.segScale=tkinter.Scale(root, command=self.setSeg, width=10, length=480, sliderlength=15, resolution=1)#滑块 self.segScale.pack(side=right) def setSeg(self, seg):#设置分段数 if self.curLineId: self.lineDictself.curLineId.seg=max(1,int(seg) self.need2update=True #print (event) def eQuit(self, event): self.running=0 def eLeftClick(self, event): pos=event.x, event.y if self.editState=new: getCtrps=self.getObjs(ctrp, pos, 3) if len(getCtrps): self.curCtrPId=getCtrps0 self.modifyPoint(self.curCtrPId, lineColor=self.colorPick) else: self.curCtrPId=self.createPoint(pos, 3, lineColor=self.colorPick) self.curCtrPsDictself.curCtrPId=pos self.curCtrPsInd.append(self.curCtrPId) if self.curLineId is None and len(self.curCtrPsInd)1: ctrps=self.sortCtrPs(self.curCtrPsDict, self.curCtrPsInd) curve=nurbs.Nurbs(ctrps) self.curLineId=self.create_line(curve.cvps, tags=curve, fill=self.colorCurve) self.lineDictself.curLineId=curve self.dotLineId=self.create_line(ctrps, dash=1, fill=self.colorDashLine) elif self.editState=select: getCurve=self.getObjs(curve, pos, 3) if len(getCurve): self.curLineId=getCurve0 #self.curCtrPsDict.clear() #self.curCtrPsInd.clear() ctrps=self.lineDictself.curLineId.ctrPoints for pos in ctrps: point=self.createPoint(pos, 3, lineColor=self.colorNormal) self.curCtrPsDictpoint=pos self.curCtrPsInd.append(point) self.modifyPoint(self.curCtrPsInd-1, lineColor=self.colorEnd) self.dotLineId=self.create_line(ctrps, dash=1, fill=self.colorDashLine) self.need2update=True self.segScale.set(self.lineDictself.curLineId.seg)#获取分段数 self.editState=new self.itemconfig(self.stateText, text=self.editState) def eLeftDrag(self, event): pos=event.x, event.y if not self.curCtrPId is None: self.curCtrPsDictself.curCtrPId=pos self.modifyPoint(self.curCtrPId, pos) self.need2update=True def eLeftRelease(self, event): pos=event.x, event.y if not self.curCtrPId is None: print (self.curCtrPId) self.curCtrPsDictself.curCtrPId=pos self.modifyPoint(self.curCtrPId, pos) if self.curCtrPId=self.curCtrPsInd-1: if len(self.curCtrPsInd)1: self.modifyPoint(self.curCtrPsInd-2, lineColor=self.colorNormal) self.modifyPoint(self.curCtrPId, lineColor=self.colorEnd) else: self.modifyPoint(self.curCtrPId, lineColor=self.colorNormal) self.curCtrPId=None self.need2update=True def eMidClick(self, event): pos=event.x, event.y if self.curLineId: self.delete(*self.curCtrPsInd) self.delete(self.dotLineId) self.curLineI