B_Spline_Approximation

B样条曲线的拟合主要是一个LSQ（least squares) 拟合问题，主要思想也是最小二乘法的思想，这与B-Spline曲线插值不同，拟合的曲线是尽量接近数据点，而不是完全通过。主要的方法可以参考cs3621

这里我定义了一个BS_curve类，类中的方法包括数据的参数化(parameterization)，节点(knots)的生成，计算系数 N i , p N_{i,p} Ni,p，De_Boor算法以及最小二乘拟合(approximation),完整代码如下：

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import mathclass BS_curve(object):def __init__(self,n,p,cp=None,knots=None):self.n = n # n+1 control points >>> p0,p1,,,pnself.p = pif cp:self.cp = cpself.u = knotsself.m = knots.shape[0]-1 # m+1 knots >>> u0,u1,,,nmelse:self.cp = Noneself.u = Noneself.m = Noneself.paras = Nonedef check(self):if self.m == self.n + self.p + 1:return 1else:return 0def coeffs(self,uq):# n+1 control points >>> p0,p1,,,pn# m+1 knots >>> u0,u1,,,nm# algorithm is from https://pages.mtu.edu/~shene/COURSES/cs3621/NOTES/spline/B-spline/bspline-curve-coef.html#N[] holds all intermediate and the final results# in fact N is longer than control points,this is just to hold the intermediate value# at last, we juest extract a part of N,that is N[0:n+1]N = np.zeros(self.m+1,dtype=np.float64) # rule out special cases Important Properties of clamped B-spline curveif uq == self.u[0]:N[0] = 1.0return N[0:self.n+1]elif uq == self.u[self.m]:N[self.n] = 1.0return N[0:self.n+1]# now u is between u0 and um# first find k uq in span [uk,uk+1)check = uq - self.uind = check >=0k = np.max(np.nonzero(ind))# sk >>> multiplicity of u[k]sk = np.sum(self.u==self.u[k])N[k] = 1.0 # degree 0# degree d goes from 1 to pfor d in range(1,self.p+1):r_max = self.m - d - 1 # the maximum subscript value of N in degree d,the minimum is 0if k-d >=0:if self.u[k+1]-self.u[k-d+1]:N[k-d] = (self.u[k+1]-uq)/(self.u[k+1]-self.u[k-d+1])*N[k-d+1] #right (south-west corner) term onlyelse:N[k-d] = (self.u[k+1]-uq)/1*N[k-d+1] #right (south-west corner) term onlyfor i in range(k-d+1,(k-1)+1):if i>=0 and i<=r_max:Denominator1 = self.u[i+d]-self.u[i]Denominator2 = self.u[i+d+1]-self.u[i+1]# 0/0=0if Denominator1 == 0:Denominator1 = 1if Denominator2 == 0:Denominator2 = 1N[i] = (uq-self.u[i])/(Denominator1)*N[i]+(self.u[i+d+1]-uq)/(Denominator2)*N[i+1]if k <= r_max:if self.u[k+d]-self.u[k]:N[k] = (uq-self.u[k])/(self.u[k+d]-self.u[k])*N[k]else:N[k] = (uq-self.u[k])/1*N[k]return N[0:self.n+1]def De_Boor(self,uq):# Input: a value u# Output: the point on the curve, C(u)# first find k uq in span [uk,uk+1)check = uq - self.uind = check >=0k = np.max(np.nonzero(ind))# inserting uq h timesif uq in self.u:# sk >>> multiplicity of u[k]sk = np.sum(self.u==self.u[k])h = self.p - skelse:sk = 0h = self.p# rule out special casesif h == -1:if k == self.p:return np.array(self.cp[0])elif k == self.m:return np.array(self.cp[-1])# initial values of P(affected control points) >>> Pk-s,0 Pk-s-1,0 ... Pk-p+1,0P = self.cp[k-self.p:k-sk+1]P = P.copy()dis = k-self.p # the index distance between storage loaction and varibale i# 1-hfor r in range(1,h+1):# k-p >> k-sktemp = [] # uesd for Storing variables of the current stagefor i in range(k-self.p+r,k-sk+1):a_ir = (uq-self.u[i])/(self.u[i+self.p-r+1]-self.u[i])temp.append((1-a_ir)*P[i-dis-1]+a_ir*P[i-dis])P[k-self.p+r-dis:k-sk+1-dis] = np.array(temp)# the last value is what we wantreturn P[-1]def bs(self,us):y = []for x in us:y.append(self.De_Boor(x))y = np.array(y)return ydef estimate_parameters(self,data_points,method="centripetal"):pts = data_points.copy()N = pts.shape[0]w = pts.shape[1]Li = []for i in range(1,N):Li.append(np.sum([pts[i,j]**2 for j in range(w)])**0.5)L = np.sum(Li)t= [0]for i in range(len(Li)):Lki = 0for j in range(i+1):Lki += Li[j]t.append(Lki/L)t = np.array(t)self.paras = tind = t>1.0t[ind] = 1.0return tdef get_knots(self,method="average"):knots = np.zeros(self.p+1).tolist()paras_temp = self.paras.copy()# m = n+p+1self.m = self.n + self.p + 1# we only need m+1 knots# so we just select m+1-(p+1)-(p+1)+(p-1)+1+1  paras to averagenum = self.m - self.p  # select n+1 parasind = np.linspace(0,paras_temp.shape[0]-1,num)ind = ind.astype(int)paras_knots = paras_temp[ind]for j in range(1,self.n-self.p+1):k_temp = 0# the maximun of variable i is n-1for i in range(j,j+self.p-1+1):k_temp += paras_knots[i]k_temp /= self.pknots.append(k_temp)add = np.ones(self.p+1).tolist()knots = knots + addknots = np.array(knots)self.u = knotsself.m = knots.shape[0]-1return knotsdef set_paras(self,parameters):self.paras = parametersdef set_knots(self,knots):self.u = knotsdef approximation(self,pts):## Obtain a set of parameters t0, ..., tn#pts_paras = self.estimate_parameters(pts)## knot vector U;#knots = self.get_knots()num = pts.shape[0]-1 # (num+1) is the number of data pointsP = np.zeros((self.n+1,pts.shape[1]),dtype=np.float64) # n+1 control pointsP[0] = pts[0]P[-1] = pts[-1]# compute NN = []for uq in self.paras:N_temp = self.coeffs(uq)N.append(N_temp)N = np.array(N)Q = [0] # hold the locationfor k in range(1,num-1+1):Q_temp = pts[k] - N[k,0]*pts[0] - N[k,self.n]*pts[-1]Q.append(Q_temp)b = [0]for i in range(1,self.n-1+1):b_temp = 0for k in range(1,num-1+1):b_temp += N[k,i]*Q[k]b.append(b_temp)b = b[1::]b = np.array(b)N = N[:,1:(self.n-1)+1]A = np.dot(N.T,N)cpm = np.linalg.solve(A,b)P[1:self.n] = cpmself.cp = Preturn Pif __name__ =="__main__":bs = BS_curve(8,3)xx = np.linspace(0,4*np.pi,101)yy = np.sin(xx)+0.6*np.random.random(101)fig = plt.figure(figsize=(10,5))ax = fig.add_subplot(111)ax.scatter(xx,yy)data = np.array([xx,yy]).Tparas = bs.estimate_parameters(data)knots = bs.get_knots()if bs.check():cp = bs.approximation(data)uq = np.linspace(0,1,101)y = bs.bs(uq)ax.plot(y[:,0],y[:,1],'-r')ax.plot(cp[:,0],cp[:,1],'-b*')plt.show()

上面代码的结果如下图所示：

B样条曲线拟合(B_Spline_Approximation)相关推荐

ITK:将样条曲线拟合到点集
ITK:将样条曲线拟合到点集内容提要 C++实现代码内容提要将样条曲线拟合到点集. C++实现代码 #include "itkBSplineScatteredDataPointSetT ...
Matlab光滑曲线多项式拟合与样条曲线拟合的两个案例
%多项式曲线拟合 figure(1) matrix2=[]; %新建空矩阵 h1=polyfit(matrix1(:,1),matrix1(:,2),3); %计算多项式拟合系数,3-拟合次数 mat ...
matlab 拟合光滑曲线图,Matlab光滑曲线多项式拟合与样条曲线拟合的两个案例
%多项式曲线拟合 figure(1) matrix2=[]; %新建空矩阵 h1=polyfit(matrix1(:,1),matrix1(:,2),3); %计算多项式拟合系数,3-拟合次数 mat ...
轨迹绕圈算法_基于三次B样条曲线拟合的智能车轨迹跟踪算法
收稿日期:2017-10-30; 修回日期:2017-12-10; 录用日期:2017-12-19. 基金项目: 国家自然科学基金资助项目( 91420202,61372088) . 作者简介: 张永 ...
三次B样条曲线拟合算法
1 三次B样条曲线方程 B样条曲线分为近似拟合和插值拟合,所谓近似拟合就是不过特征点,而插值拟合就是通过特征点,但是插值拟合需要经过反算得到控制点再拟合出过特征点的B样条曲线方程.这里会一次介绍两种拟 ...
Matlab中对离散数据点进行B样条曲线拟合
1.拟合出的曲线通过离散的路径点 x= [0;0.0128205128205128;0.0256410256410256;0.0384615384615385;0.0512820512820513;0 ...
PCL B样条曲线拟合（2d/3d）
文章目录一.简介 1.1定义 1.2最小二乘拟合二.实现代码(PCL) 三.实现效果参考资料一.简介 1.1定义一条B样条曲线可以被定义成 n + 1 n+1 n+
开源项目推荐：Bezier曲线、B-Spline和NURBS的区别与《THE NURBS BOOK 2nd》简介，曲线拟合可视化工具
一.基本概念 B-Spline:B样条曲线 NURBS(Non Uniform Rational B-Spline):非均匀有理B样条曲线 B样条曲线有三种类型: 当起始点和终止点的重复度为最高次数加 ...
三次Beizer曲线拟合算法
1 三次Beizer曲线方程介绍 Beizer曲线的一些特性这里不再赘述,大家可以去网上查看一些资料,很详细.最近用到轮廓拟合,所以用三次Beizer曲线效果还可以,有插值和近似拟合(插值就是曲线过点 ...

B样条曲线拟合(B_Spline_Approximation)

B_Spline_Approximation

B样条曲线拟合(B_Spline_Approximation)相关推荐

最新文章

热门文章