• B样条曲线拟合(B_Spline_Approximation)
千次阅读 热门讨论
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} ，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)

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.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样条曲线拟合。实现简单高效，项目中实际使用的代码。
• b样条曲线拟合，用着很好用
• 代码是基于C++的三次B样条曲线拟合代码，包含插值拟合，近似拟合就不放代码了，较简单，我的博客中有相关论文链接。http://blog.csdn.net/liumangmao1314/article/details/54588155
• 本代码为三次B样条曲线拟合10个点的程序，并利用OpenGL对最后拟合的曲线进行绘制，注意：需要安装EIGEN矩阵库和OpenGL才能运行。
• 解压密码为：hur.cn 主要采用C++编程实现，过控制点的三次B样条曲线拟合，可以用于各种高级的曲线拟合方面。
• 在工程上往往需要进行三维曲线拟合，该Matlab代码可以对三组离散数据进行三维B样条曲线拟合
• 以渐开线为例，实现b样条曲线插值拟合，可以任意改变曲线类型
• 三次B样条曲线拟合C++ B样条曲线的方程：P=∑i=0nPiFi,k(t)\sum_{i=0}^nP_iF_{i,k}(t)∑i=0n​Pi​Fi,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)
其中 F i , k ( t ) F_{i,k}(t) 为基函数，三次B样条的基函数为：
F 0 , 3 ( t ) = 1 6 ( 1 − t ) 3 F_{0,3}(t)=\displaystyle{\frac{1}{6}{(1-t)}^3}

F 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)}

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)}

F 3 , 3 ( t ) = 1 6 t 3 F_{3,3}(t)=\displaystyle{\frac{1}{6}t^3}

所以，三次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 = w 0 + w 1 t + w 2 t 2 + w 3 t 3 P=w_0+w_1t+w_2t^2+w_3t^3 (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)}

w 1 = − 1 2 ( P 0 − P 2 ) w_1=\displaystyle{-\frac{1}{2}(P_0-P_2)}

w 2 = 1 2 ( P 0 − 2 P 1 + P 2 ) w_2=\displaystyle{\frac{1}{2}(P_0-2P_1+P_2)}

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)}

因此，每四个离散点就可拟合一段曲线，比如 P 0 , P 1 , P 2 , P 3 P_0,P_1,P_2,P_3 可以拟合一段光滑曲线， P 1 , P 2 , P 3 , P 4 P_1,P_2,P_3,P_4 可以拟合下一段，相邻两段曲线是平滑过渡的，以此类推N个点可以拟合出N-3段平滑相接的曲线。
以上是非闭合曲线的拟合，闭合曲线只需离散点集首尾相连，也就是说，还需用 P N − 1 , P 0 , P 1 , P 2 P_{N-1},P_0,P_1,P_2 拟合一段曲线。
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;
}


非闭合拟合效果：

闭合曲线拟合效果：

展开全文
• 我刚查了下，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样条(NURBS)曲线的算法，实现了满足一定精度要求的数据点云拟合以及控制点的计算。该方法首先通过点云外形特征提取主特征点，把主特征点...
• 基于B样条曲线拟合出现的问题和困难,提出了一种新的B样条曲线拟合方法.该方法成功地避免了数据点参数化的问题,并使得逼近曲线具有较好的形状和接近弧长参数化的节点向量.
• %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样条曲线拟合问题中，将节点作为自由变量可大幅提高拟合精度，但这就使曲线拟合问题转化为求解困难的连续多峰值、多变量非线性优化问题，当待拟合的曲线是不连续、有尖点情况，就更为困难。针对这一问题，基于混沌...
• ## 三次B样条曲线拟合算法

万次阅读 多人点赞 2017-01-17 22:10:28
三次B样条曲线方程B样条曲线分为近似拟合和插值拟合，所谓近似拟合就是不过特征点，而插值拟合就是通过特征点，但是插值拟合需要经过反算得到控制点再拟合出过特征点的B样条曲线方程。这里会一次介绍两种拟合算法。...
• 基于三次B样条曲线拟合的智能车轨迹跟踪算法
• 绘制B样条曲线，可以修改参数，给出控制点，拟合b样条
• b样条曲线（B-spline curve）是指在数学的子学科数值分析里的一种特殊的表示形式。它是B-样条基曲线的线性组合。由Isaac Jacob Schoenberg创造。 B-样条是贝兹曲线(又称贝塞尔曲线)的一种一般化，可以进一步推广为非...
• 由于物体的结构约束导致通过一个横截面获得的点云通常不会是规则的圆形，那么求这个维度就不能直接用圆形周长公式，因此这个时候可以利用pcl库中的B样条曲线进行曲线拟合。 具体曲线拟合的实现在pcl利用的是ON_...
• 例程： 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样条曲线。 ...
• % ref: 闭合 B 样条曲线控制点的快速求解算法及应用 % http://www.doc88.com/p-5714423317458.html % https://blog.csdn.net/liumangmao1314/article/details/54588155 ========================...
• %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
• 基于粒子群算法的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:43
B样条曲线是在Bezier 曲线基础上发展起来的一类曲线，它克服了Bezier 曲线整体控制性所带来的不便，最常用的是二次和三次B样条曲线。 2.二次B样条 2.1 参数方程 已知三个平面离散点P0、P1、P2，由这三点可以定义...
• 图像二值化之后，可以用cv2.CHAIN_APPROX_NONE或者cv2....再对上面保存的csv文件读取，进行B样条曲线拟合。可以修改拟合点的数量。 import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d
• 为解决数控系统进行连续微线段加工时加减速频繁、运行速度缓慢、加工路径不连续等问题,提出了最小二乘3次B样条曲线逼近拟合算法,采用该算法实现对连续微线段的逼近.文中通过分析连续微线段加工路径的几何特性,提出了...
• 基于三次B样条曲线拟合的列车定位方法研究

...