-
2020-07-06 23:44:24
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 math class BS_curve(object): def __init__(self,n,p,cp=None,knots=None): self.n = n # n+1 control points >>> p0,p1,,,pn self.p = p if cp: self.cp = cp self.u = knots self.m = knots.shape[0]-1 # m+1 knots >>> u0,u1,,,nm else: self.cp = None self.u = None self.m = None self.paras = None def check(self): if self.m == self.n + self.p + 1: return 1 else: return 0 def 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 curve if uq == self.u[0]: N[0] = 1.0 return N[0:self.n+1] elif uq == self.u[self.m]: N[self.n] = 1.0 return N[0:self.n+1] # now u is between u0 and um # first find k uq in span [uk,uk+1) check = uq - self.u ind = check >=0 k = 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 p for 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 0 if 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 only else: N[k-d] = (self.u[k+1]-uq)/1*N[k-d+1] #right (south-west corner) term only for 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=0 if Denominator1 == 0: Denominator1 = 1 if Denominator2 == 0: Denominator2 = 1 N[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.u ind = check >=0 k = np.max(np.nonzero(ind)) # inserting uq h times if uq in self.u: # sk >>> multiplicity of u[k] sk = np.sum(self.u==self.u[k]) h = self.p - sk else: sk = 0 h = self.p # rule out special cases if 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,0 P = 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-h for r in range(1,h+1): # k-p >> k-sk temp = [] # uesd for Storing variables of the current stage for 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 want return P[-1] def bs(self,us): y = [] for x in us: y.append(self.De_Boor(x)) y = np.array(y) return y def 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 = 0 for j in range(i+1): Lki += Li[j] t.append(Lki/L) t = np.array(t) self.paras = t ind = t>1.0 t[ind] = 1.0 return t def get_knots(self,method="average"): knots = np.zeros(self.p+1).tolist() paras_temp = self.paras.copy() # m = n+p+1 self.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 average num = self.m - self.p # select n+1 paras ind = 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-1 for i in range(j,j+self.p-1+1): k_temp += paras_knots[i] k_temp /= self.p knots.append(k_temp) add = np.ones(self.p+1).tolist() knots = knots + add knots = np.array(knots) self.u = knots self.m = knots.shape[0]-1 return knots def set_paras(self,parameters): self.paras = parameters def set_knots(self,knots): self.u = knots def 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 points P = np.zeros((self.n+1,pts.shape[1]),dtype=np.float64) # n+1 control points P[0] = pts[0] P[-1] = pts[-1] # compute N N = [] for uq in self.paras: N_temp = self.coeffs(uq) N.append(N_temp) N = np.array(N) Q = [0] # hold the location for 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 = 0 for 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] = cpm self.cp = P return P if __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]).T paras = 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样条曲线拟合
2019-01-15 12:07:57B样条曲线拟合。实现简单高效,项目中实际使用的代码。 -
b样条曲线拟合
2019-03-31 15:53:54b样条曲线拟合,用着很好用 -
基于C++的三次B样条曲线拟合代码
2017-06-08 22:48:24代码是基于C++的三次B样条曲线拟合代码,包含插值拟合,近似拟合就不放代码了,较简单,我的博客中有相关论文链接。http://blog.csdn.net/liumangmao1314/article/details/54588155 -
三次B样条曲线拟合点集程序(C++)
2018-01-04 17:54:08本代码为三次B样条曲线拟合10个点的程序,并利用OpenGL对最后拟合的曲线进行绘制,注意:需要安装EIGEN矩阵库和OpenGL才能运行。 -
C++ 过控制点的三次B样条曲线拟合
2017-12-20 10:41:13解压密码为:hur.cn 主要采用C++编程实现,过控制点的三次B样条曲线拟合,可以用于各种高级的曲线拟合方面。 -
三维B样条曲线拟合Matlab程序
2015-04-23 19:52:57在工程上往往需要进行三维曲线拟合,该Matlab代码可以对三组离散数据进行三维B样条曲线拟合 -
jiankaixian--B-yangtiao_B样条_渐开线_B样条曲线拟合_
2021-10-01 17:07:28以渐开线为例,实现b样条曲线插值拟合,可以任意改变曲线类型 -
三次B样条曲线拟合C++
2020-08-21 11:20:17三次B样条曲线拟合C++ B样条曲线的方程:P=∑i=0nPiFi,k(t)\sum_{i=0}^nP_iF_{i,k}(t)∑i=0nPiFi,k(t) 其中Fi,k(t)F_{i,k}(t)Fi,k(t)为基函数,三次B样条的基函数为: F0,3(t)=16(1−t)3F_{0,3}(t)=\...B样条曲线的方程:P= ∑ i = 0 n P i F i , k ( t ) \sum_{i=0}^nP_iF_{i,k}(t) ∑i=0nPiFi,k(t)
其中 F i , k ( t ) F_{i,k}(t) Fi,k(t)为基函数,三次B样条的基函数为:
F 0 , 3 ( t ) = 1 6 ( 1 − t ) 3 F_{0,3}(t)=\displaystyle{\frac{1}{6}{(1-t)}^3} F0,3(t)=61(1−t)3F 1 , 3 ( t ) = 1 6 ( 3 t 3 − 6 t 2 + 4 ) F_{1,3}(t)=\displaystyle{\frac{1}{6}(3t^3-6t^2+4)} F1,3(t)=61(3t3−6t2+4)
F 2 , 3 ( t ) = 1 6 ( − 3 t 3 + 3 t 2 + 3 t + 1 ) F_{2,3}(t)=\displaystyle{\frac{1}{6}(-3t^3+3t^2+3t+1)} F2,3(t)=61(−3t3+3t2+3t+1)
F 3 , 3 ( t ) = 1 6 t 3 F_{3,3}(t)=\displaystyle{\frac{1}{6}t^3} F3,3(t)=61t3
所以,三次B样条的方程式为:
P = P 0 F 0 , 3 ( t ) + P 1 F 1 , 3 ( t ) + P 2 F 2 , 3 ( t ) + P 3 F 3 , 3 ( t ) P=P_0F_{0,3}(t)+P_1F_{1,3}(t)+P_2F_{2,3}(t)+P_3F_{3,3}(t) P=P0F0,3(t)+P1F1,3(t)+P2F2,3(t)+P3F3,3(t)
把基函数代入可以简化为:
P = w 0 + w 1 t + w 2 t 2 + w 3 t 3 P=w_0+w_1t+w_2t^2+w_3t^3 P=w0+w1t+w2t2+w3t3 (0≤t<≤1)
其中, w 0 = 1 6 ( P 0 + 4 P 1 + P 2 ) w_0=\displaystyle{\frac{1}{6}(P_0+4P_1+P_2)} w0=61(P0+4P1+P2)w 1 = − 1 2 ( P 0 − P 2 ) w_1=\displaystyle{-\frac{1}{2}(P_0-P_2)} w1=−21(P0−P2)
w 2 = 1 2 ( P 0 − 2 P 1 + P 2 ) w_2=\displaystyle{\frac{1}{2}(P_0-2P_1+P_2)} w2=21(P0−2P1+P2)
w 3 = − 1 6 ( P 0 − 3 P 1 + 3 P 2 − P 3 ) w_3=\displaystyle{-\frac{1}{6}(P_0-3P_1+3P_2-P_3)} w3=−61(P0−3P1+3P2−P3)
因此,每四个离散点就可拟合一段曲线,比如 P 0 , P 1 , P 2 , P 3 P_0,P_1,P_2,P_3 P0,P1,P2,P3可以拟合一段光滑曲线, P 1 , P 2 , P 3 , P 4 P_1,P_2,P_3,P_4 P1,P2,P3,P4可以拟合下一段,相邻两段曲线是平滑过渡的,以此类推N个点可以拟合出N-3段平滑相接的曲线。
以上是非闭合曲线的拟合,闭合曲线只需离散点集首尾相连,也就是说,还需用 P N − 1 , P 0 , P 1 , P 2 P_{N-1},P_0,P_1,P_2 PN−1,P0,P1,P2拟合一段曲线。
c++代码如下:/*B样条曲线拟合 @return 返回拟合得到的曲线 @discretePoints 输入的离散点,至少4个点 @closed 是否拟合闭合曲线,true表示闭合,false不闭合 @stride 拟合精度 */ vector<Point2f> BSplineFit(vector<Point2f> discretePoints, bool closed, double stride = 0.01) { vector<Point2f> fittingPoints; for (int i = 0; i < (closed ? discretePoints.size() : discretePoints.size() - 1); i++) { Point2f xy[4]; xy[0] = (discretePoints[i] + 4 * discretePoints[(i + 1) % discretePoints.size()] + discretePoints[(i + 2) % discretePoints.size()]) / 6; xy[1] = -(discretePoints[i] - discretePoints[(i + 2) % discretePoints.size()]) / 2; xy[2] = (discretePoints[i] - 2 * discretePoints[(i + 1) % discretePoints.size()] + discretePoints[(i + 2) % discretePoints.size()]) / 2; xy[3] = -(discretePoints[i] - 3 * discretePoints[(i + 1) % discretePoints.size()] + 3 * discretePoints[(i + 2) % discretePoints.size()] - discretePoints[(i + 3) % discretePoints.size()]) / 6; for (double t = 0; t <= 1; t += stride) { Point2f totalPoints = Point2f(0, 0); for (int j = 0; j < 4; j++) { totalPoints += xy[j] * pow(t, j); } fittingPoints.push_back(totalPoints); } } return fittingPoints; }
非闭合拟合效果:
闭合曲线拟合效果:
-
请matlab高手过来看看 怎么用b样条曲线拟合离散点
2021-03-03 17:08:50我刚查了下,b样条曲线拟合就是拟合成光滑曲线。这里可以尝试Matlab的polyfit命令,我尝试了好几个,发现在5阶的时候已经非常接近了,当然如果你需要更高精度,可以继续提高阶次。代码:x=[1:20];y=[42 45 47 49 52 ...我刚查了下,b样条曲线拟合就是拟合成光滑曲线。这里可以尝试Matlab的polyfit命令,我尝试了好几个,发现在5阶的时候已经非常接近了,当然如果你需要更高精度,可以继续提高阶次。
代码:
x=[1:20];
y=[42 45 47 49 52 59 66 74 85 98 111 125 136 147 157 162 164 167 168 168];
plot(x,y,'r')
hold on
p=polyfit(x,y,5)
z=p(1)*x.^5+p(2)*x.^4+p(3)*x.^3+p(4)*x.^2+p(5)*x+p(6);
plot(x,z,'b')
legend('红色原来数据曲线','蓝色直接模拟曲线')
输出结果:
p =
0.0006 -0.0315 0.5628 -3.4653 10.5082 34.1178
所以拟合结果是:
这是个人愚见,希望对你有帮助,有疑问请追问,若满意还望采纳,祝生活愉快!
-
基于最少控制点的非均匀有理B样条曲线拟合 (2008年)
2021-05-17 11:38:44针对叶片型线的优化设计,提出采用自适应方法提取合适的节点来插值非均匀有理B样条(NURBS)曲线的算法,实现了满足一定精度要求的数据点云拟合以及控制点的计算。该方法首先通过点云外形特征提取主特征点,把主特征点... -
B样条曲线拟合相关资料
2014-09-13 09:36:04基于B样条曲线拟合出现的问题和困难,提出了一种新的B样条曲线拟合方法.该方法成功地避免了数据点参数化的问题,并使得逼近曲线具有较好的形状和接近弧长参数化的节点向量. -
求解分段方程的Matlab三次B样条曲线拟合算法,matlab
2020-12-24 11:35:08%hg为x和y拟合的系数,16行, %第1,2行分别为第一段x,y的系数,3,4为第二段,类推 plot(x,y) hold on plot([xx(i+1),xx(i+2)],[yy(i+1),yy(i+2)]) hold on end 结果 方程参数为hg 第一行所有为第一个方程横坐标x...参考
代码如下:
clc
clear
xx=[6.852,5.934,5.317,4.617,3.924,3.232,2.525,1.882,0.999];
yy=[1.399,1.399,1.226,0.859,0.212,0.339,-0.657,-0.892,-0.892];
xx=[xx(1)-xx(2)+xx(1),xx,xx(end)-xx(end-1)+xx(end)];
yy=[yy(1)-yy(2)+yy(1),yy,yy(end)-yy(end-1)+yy(end)];
hg=[];
for i=1:8
t=(0:0.001:1);
x0=xx(i);x1=xx(i+1);x2=xx(i+2);x3=xx(i+3);
y0=yy(i);y1=yy(i+1);y2=yy(i+2);y3=yy(i+3);
a0=(x0+4*x1+x2)/6;a1=-(x0-x2)/2;a2=(x0-2*x1+x2)/2;a3=-(x0-3*x1+3*x2-x3)/6;
b0=(y0+4*y1+y2)/6;b1=-(y0-y2)/2;b2=(y0-2*y1+y2)/2;b3=-(y0-3*y1+3*y2-y3)/6;
x=a0+a1*t+a2*t.^2+a3*t.^3;
y=b0+b1*t+b2*t.^2+b3*t.^3;
hg=[hg;a0,a1,a2,a3;b0,b1,b2,b3];%hg为x和y拟合的系数,16行,
%第1,2行分别为第一段x,y的系数,3,4为第二段,类推
plot(x,y)
hold on
plot([xx(i+1),xx(i+2)],[yy(i+1),yy(i+2)])
hold on
end
结果
方程参数为hg
第一行所有为第一个方程横坐标x的常数项、一次项、二次项、三次项。
第二行所有为第一个方程纵坐标y的常数项、一次项、二次项、三次项。
3、4为第二个方程。
直到最后
-
论文研究-混沌蚂蚁群优化求解自由节点B样条曲线拟合.pdf
2019-09-10 14:23:29B样条曲线拟合问题中,将节点作为自由变量可大幅提高拟合精度,但这就使曲线拟合问题转化为求解困难的连续多峰值、多变量非线性优化问题,当待拟合的曲线是不连续、有尖点情况,就更为困难。针对这一问题,基于混沌... -
三次B样条曲线拟合算法
2017-01-17 22:10:28三次B样条曲线方程B样条曲线分为近似拟合和插值拟合,所谓近似拟合就是不过特征点,而插值拟合就是通过特征点,但是插值拟合需要经过反算得到控制点再拟合出过特征点的B样条曲线方程。这里会一次介绍两种拟合算法。... -
基于三次B样条曲线拟合的智能车轨迹跟踪算法
2021-12-02 22:20:21基于三次B样条曲线拟合的智能车轨迹跟踪算法 -
B样条曲线_B样条_B样条曲线拟合_B样条曲线_点拟合
2016-12-06 21:30:15绘制B样条曲线,可以修改参数,给出控制点,拟合b样条 -
B样条曲线拟合数据
2022-05-11 17:59:40b样条曲线(B-spline curve)是指在数学的子学科数值分析里的一种特殊的表示形式。它是B-样条基曲线的线性组合。由Isaac Jacob Schoenberg创造。 B-样条是贝兹曲线(又称贝塞尔曲线)的一种一般化,可以进一步推广为非... -
pcl实现三次B样条曲线拟合不规则圆并计算所得拟合曲线长度
2022-04-05 17:20:22由于物体的结构约束导致通过一个横截面获得的点云通常不会是规则的圆形,那么求这个维度就不能直接用圆形周长公式,因此这个时候可以利用pcl库中的B样条曲线进行曲线拟合。 具体曲线拟合的实现在pcl利用的是ON_... -
Matlab中对离散数据点进行B样条曲线拟合
2020-09-09 17:44:50例程: x = [1:20]; y = [42 45 47 49 52 59 66 74 85 98 ...其中:(x,y)为离散的数据点,spapi(3,x,y)表示用3阶B样条曲线对离散的数据点进行拟合,并且要经过给定的离散数据点,fnplt(sp)代表画出该B样条曲线。 ... -
三维闭合B样条曲线拟合算法Matlab代码
2019-05-23 22:26:13% ref: 闭合 B 样条曲线控制点的快速求解算法及应用 % http://www.doc88.com/p-5714423317458.html % https://blog.csdn.net/liumangmao1314/article/details/54588155 ========================... -
matlab三次B样条曲线拟合算法求分段方程
2020-11-29 23:08:16%hg为x和y拟合的系数,16行, %第1,2行分别为第一段x,y的系数,3,4为第二段,类推 plot(x,y) hold on plot([xx(i+1),xx(i+2)],[yy(i+1),yy(i+2)]) hold on end 结果 方程参数为hg 第一行所有为第一个方程横坐标x... -
三维航迹的B样条曲线拟合算法.pdf
2021-10-11 04:32:22三维航迹的B样条曲线拟合算法.pdf -
基于粒子群算法的B样条曲线拟合.pdf
2021-09-29 10:08:34基于粒子群算法的B样条曲线拟合.pdf -
PCL学习:平面点云B样条曲线拟合
2019-07-10 17:58:23//初始化曲线拟合对象 fit.assemble (curve_params); //装配曲线参数 fit.solve (); //拟合 VisualizeCurve (fit.m_nurbs, 1.0, 0.0, 0.0, false); //可视化拟合曲线 } // visualize viewer.setSize ... -
B样条曲线拟合原理
2017-01-13 18:51:43B样条曲线是在Bezier 曲线基础上发展起来的一类曲线,它克服了Bezier 曲线整体控制性所带来的不便,最常用的是二次和三次B样条曲线。 2.二次B样条 2.1 参数方程 已知三个平面离散点P0、P1、P2,由这三点可以定义... -
得到图像轮廓点的像素坐标,用B样条曲线拟合,并求曲率
2021-02-05 11:00:04图像二值化之后,可以用cv2.CHAIN_APPROX_NONE或者cv2....再对上面保存的csv文件读取,进行B样条曲线拟合。可以修改拟合点的数量。 import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d -
数控加工中连续微线段轨迹的B样条曲线拟合 (2012年)
2021-06-15 12:23:21为解决数控系统进行连续微线段加工时加减速频繁、运行速度缓慢、加工路径不连续等问题,提出了最小二乘3次B样条曲线逼近拟合算法,采用该算法实现对连续微线段的逼近.文中通过分析连续微线段加工路径的几何特性,提出了... -
基于三次B样条曲线拟合的列车定位方法研究
2021-04-05 06:40:37基于三次B样条曲线拟合的列车定位方法研究