精华内容
下载资源
问答
  • 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} 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()
    

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

    更多相关内容
  • B样条曲线拟合

    2019-01-15 12:07:57
    B样条曲线拟合。实现简单高效,项目中实际使用的代码。
  • b样条曲线拟合

    2019-03-31 15:53:54
    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++

    千次阅读 2020-08-21 11:20:17
    三次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) 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(1t)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)} F1,3(t)=61(3t36t2+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(P0P2)

    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(P02P1+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(P03P1+3P2P3)

    因此,每四个离散点就可拟合一段曲线,比如 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 PN1,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;
    }
    

    非闭合拟合效果:

    在这里插入图片描述
    闭合曲线拟合效果:
    在这里插入图片描述
    在这里插入图片描述

    展开全文
  • 我刚查了下,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

    5c63fef9ef01e4cbaf5860e8b52e029a.png

    所以拟合结果是:

    4a0a5bf7509123ed7a453a4bba9d716d.png

    这是个人愚见,希望对你有帮助,有疑问请追问,若满意还望采纳,祝生活愉快!

    5557fc66944a353f92513d4afeb39ece.png

    b8bc992455007ec622ee9d37f0fe3a85.png

    展开全文
  • 针对叶片型线的优化设计,提出采用自适应方法提取合适的节点来插值非均匀有理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样条曲线拟合数据

    2022-05-11 17:59:40
    b样条曲线(B-spline curve)是指在数学的子学科数值分析里的一种特殊的表示形式。它是B-样条基曲线的线性组合。由Isaac Jacob Schoenberg创造。 B-样条是贝兹曲线(又称贝塞尔曲线)的一种一般化,可以进一步推广为非...
  • 由于物体的结构约束导致通过一个横截面获得的点云通常不会是规则的圆形,那么求这个维度就不能直接用圆形周长公式,因此这个时候可以利用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样条曲线。 ...
  • % 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样条曲线拟合的列车定位方法研究

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 5,255
精华内容 2,102
关键字:

b样条曲线拟合