1 三次B样条曲线方程
B样条曲线分为近似拟合和插值拟合，所谓近似拟合就是不过特征点，而插值拟合就是通过特征点，但是插值拟合需要经过反算得到控制点再拟合出过特征点的B样条曲线方程。这里会一次介绍两种拟合算法。首先介绍B样条的曲线方程。
B样条曲线的总方程为：P(t)=∑i=0nPiFi,k(t)P(t)=\sum_{i=0}^{n} P_{i}F_{i,k}(t)P(t)=∑i=0n​Pi​Fi,k​(t)    (1)

其中PiP_iPi​是控制曲线的特征点，Fi,k(u)F_{i,k}(u)Fi,k​(u)则是K阶B样条基函数。
1.1 三次B样条曲线方程中基函数为：
Fi,k(t)=1k!∑m=0k−i(−1)m(mk+1)(t+k−m−j)kF_{i,k}(t)=\frac{1}{k!}\sum_{m=0}^{k-i}(-1)^{m}\binom{m}{k+1}(t+k-m-j)^kFi,k​(t)=k!1​∑m=0k−i​(−1)m(k+1m​)(t+k−m−j)k  (2)

其中(mk+1)\binom{m}{k+1}(k+1m​)表示阶乘，化成看的明白的式子就是：

（


将图片上的基函数代入到方程（1）中，就是：
P(t)=P0∗F0,3(t)+P1∗F1,3(t)+P2∗F2,3(t)+P3∗F3,3(t)P(t)= 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(t)=P0​∗F0,3​(t)+P1​∗F1,3​(t)+P2​∗F2,3​(t)+P3​∗F3,3​(t)  (3)

方程（3）就是三次B样条曲线方程。

2019-04-18 更新

有小伙伴提到上式(2)的j是什么意思，其实j就是控制点的索引值。这里我把书上的公式摘抄下来，可能看起来更为清晰。

参考书籍：《计算机图形学 第3版》何援军 第13章
三次B样条曲线计算
F0,3(t)=13!∑j=03(−1)jC4j)(t+3−0−j)3F_{0,3}(t)=\frac{1}{3!}\sum_{j=0}^{3}(-1)^{j}C^{j}_{4})(t+3-0-j)^3F0,3​(t)=3!1​∑j=03​(−1)jC4j​)(t+3−0−j)3
=16[(−1)0C40(t+3)3+(−1)1C41(t+2)3+(−1)2C42(t+1)3+(−1)3C43t3]=\frac{1}{6}[(-1)^{0}C^{0}_{4}(t+3)^{3}+(-1)^{1}C^{1}_{4}(t+2)^{3}+(-1)^{2} C^{2}_{4}(t+1)^{3}+(-1)^{3}C^{3}_{4}t^{3}]=61​[(−1)0C40​(t+3)3+(−1)1C41​(t+2)3+(−1)2C42​(t+1)3+(−1)3C43​t3]
=16(−t3+3t2−3t+1)=16(1−t)3=\frac{1}{6}(-t^{3}+3t^{2}-3t+1)=\frac{1}{6}(1-t)^{3}=61​(−t3+3t2−3t+1)=61​(1−t)3
F1,3(t)=13!∑j=02(−1)jC4j)(t+3−1−j)3F_{1,3}(t)=\frac{1}{3!}\sum_{j=0}^{2}(-1)^{j}C^{j}_{4})(t+3-1-j)^3F1,3​(t)=3!1​∑j=02​(−1)jC4j​)(t+3−1−j)3
=16[(−1)0C40(t+2)3+(−1)1C41(t+1)3+(−1)2C42t3]=\frac{1}{6}[(-1)^{0}C^{0}_{4}(t+2)^{3}+(-1)^{1}C^{1}_{4}(t+1)^{3}+(-1)^{2} C^{2}_{4}t^{3}]=61​[(−1)0C40​(t+2)3+(−1)1C41​(t+1)3+(−1)2C42​t3]
=16(3t3−6t2+4)=\frac{1}{6}(3t^{3}-6t^{2}+4)=61​(3t3−6t2+4)
F2,3(t)=13!∑j=01(−1)jC4j)(t+3−2−j)3F_{2,3}(t)=\frac{1}{3!}\sum_{j=0}^{1}(-1)^{j}C^{j}_{4})(t+3-2-j)^3F2,3​(t)=3!1​∑j=01​(−1)jC4j​)(t+3−2−j)3
=16[(−1)0C40(t+1)3+(−1)1C41t3=\frac{1}{6}[(-1)^{0}C^{0}_{4}(t+1)^{3}+(-1)^{1}C^{1}_{4}t^{3}=61​[(−1)0C40​(t+1)3+(−1)1C41​t3
=16(−3t3+3t2+3t+1)=\frac{1}{6}(-3t^{3}+3t^{2}+3t+1)=61​(−3t3+3t2+3t+1)
F3,3(t)=13!∑j=00(−1)jC4j)(t+3−3−j)3F_{3,3}(t)=\frac{1}{3!}\sum_{j=0}^{0}(-1)^{j}C^{j}_{4})(t+3-3-j)^3F3,3​(t)=3!1​∑j=00​(−1)jC4j​)(t+3−3−j)3
=16[(−1)0C40t3=\frac{1}{6}[(-1)^{0}C^{0}_{4}t^{3}=61​[(−1)0C40​t3
=16t3=\frac{1}{6}t^{3}=61​t3
2 三次B样条曲线近似拟合
近似拟合很简单。不需要求控制点，求得Fi,k(t)F_{i,k}(t)Fi,k​(t)，由上述方程（3），代入P0,P1,P2,P3P_0,P_1,P_2,P_3P0​,P1​,P2​,P3​就可以得到由这四个点近似拟合的一段三次B样条曲线，起始点在P0P_0P0​，终点在P1P_1P1​，对于闭合轮廓，最后一段可以取前两点做辅助，拟合实验结果我最后一块给出。这种近似拟合曲线光滑，但是最大不足就是不过特征点，也就是不过PiP_iPi​，需要过点需要反求控制点再拟合。
3 三次B样条插值拟合
插值拟合较为复杂。其实也不算是很复杂，找资料过程和理解过程是一个复杂的过程。不过有了前面大神做工作，我们只是借用别人的成果写代码就好了。我给大家看一篇论文，大家可以百度或者去知网搜索，闭合 B 样条曲线控制点的快速求解算法及应用。文章讲解了反求控制点的具体步骤，写的非常详细，基本上贴近代码的那种。大家可以根据这篇论文反求控制点，拟合出来的三次B样条曲线是经过PiP_iPi​的。代码就不放了，很多，可以根据我给的那篇论文直接编写相应代码，有问题可以私信我，知无不言。

##4 拟合结果
    原轮廓
   近似拟合轮廓。可以看到没过黑色特征点，只是近似拟合
  插值拟。可以看到曲线经过黑色特征点，不过有一些不足之处。
##5 总结

三次B样条曲线拟合轮廓效果还是可以，较之Beizer（可以参考我博客三次Beizer曲线拟合算法），B样条将一些细节描述的很好，很多细节之处都贴近原轮廓，但是有一些不足之处，可以看到对直线拟合效果不是很好。两篇博客都是关于闭合轮廓的拟合，对于非闭合或者只是一段曲线拟合，还有一种曲线是很好的，《数值分析》提到过，叫三次样条插值拟合，拟合效果很好，我做过拟合一元三次方程曲线，拟合效果跟原曲线非常贴近，不过过程中需要用到追赶法，而追赶法需要满足一个条件，对于闭合曲线三次样条插值是不满足这个条件的，所以我没去深研究，大家可以去试一试。谢谢大家！
-------------------------------------------------2018-10-30--------------------------------------------

真的很感谢大家的支持，这一年都比较忙，找实习and找工作，而且这个东西是研一上学期搞的，现在看都毫无印象，对自己说一句：牛逼。我把代码放在网盘了（能运行，但现在没效果，大家可以检查下，之前是OK的），本来想上传到CSDN，好像要C币，而且要审核，太墨迹了。

链接: https://pan.baidu.com/s/1mSQMmvL71gwEAqgiT6O9Gg 提取码: xv5f

B样条曲线的方程：P=$\sum_{i=0}^nP_iF_{i,k}(t)$

其中$F_{i,k}(t)$为基函数，三次B样条的基函数为：

$F_{0,3}(t)=\displaystyle{\frac{1}{6}{(1-t)}^3}$
$F_{1,3}(t)=\displaystyle{\frac{1}{6}(3t^3-6t^2+4)}$
$F_{2,3}(t)=\displaystyle{\frac{1}{6}(-3t^3+3t^2+3t+1)}$
$F_{3,3}(t)=\displaystyle{\frac{1}{6}t^3}$
所以，三次B样条的方程式为:

$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_1t+w_2t^2+w_3t^3$         (0≤t<≤1)

其中，$w_0=\displaystyle{\frac{1}{6}(P_0+4P_1+P_2)}$
$w_1=\displaystyle{-\frac{1}{2}(P_0-P_2)}$
$w_2=\displaystyle{\frac{1}{2}(P_0-2P_1+P_2)}$
$w_3=\displaystyle{-\frac{1}{6}(P_0-3P_1+3P_2-P_3)}$
因此，每四个离散点就可拟合一段曲线，比如$P_0,P_1,P_2,P_3$可以拟合一段光滑曲线，$P_1,P_2,P_3,P_4$可以拟合下一段，相邻两段曲线是平滑过渡的，以此类推N个点可以拟合出N-3段平滑相接的曲线。

以上是非闭合曲线的拟合，闭合曲线只需离散点集首尾相连，也就是说，还需用$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_Spline_Approximation
B样条曲线的拟合主要是一个LSQ（least squares) 拟合问题，主要思想也是最小二乘法的思想，这与B-Spline曲线插值不同，拟合的曲线是尽量接近数据点，而不是完全通过。主要的方法可以参考cs3621
这里我定义了一个BS_curve类，类中的方法包括数据的参数化(parameterization)，节点(knots)的生成，计算系数$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)

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

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