【Python】 NURBS/BSpline（B样条） 与geomdl 实现效果一致，纯干货（效果图+Python源代码）（动态规划）

上行为geomdl库实现效果，下行为本文的方法实现效果

左列为open模式，中间为clame模式，右侧为close模式

说明：当W不指定的时候NURBS等价于BSpline

from matplotlib import pyplot as plt 
import numpy as np
from geomdl import NURBS 
from geomdl.visualization import VisMPL

class NURBSCurve:
    def __init__(self, T, P, W=None):
        self.reset(T, P, W)

    def reset(self, T, P, W=None): 
        self.T = T
        self.P = P
        if W is None:
            self.W = np.ones_like(P)
        else:
            self.W = W
        self.m = len(T)
        self.n = len(P)
        self.k = self.m - self.n - 1 
        self.dpN = np.zeros(self.m) # 动态规划求N的
        self.P2 = None

    def caculate_dpN(self, t):
        for k in range(self.k+1): # 动态规划
            for i in range(self.n): 
                if k == 0:
                    self.dpN[i] = 1 * (self.T[i] <= t < self.T[i+1])
                else:
                    w1 = w2 = 0
                    if self.T[i+k] != self.T[i]:
                        w1 = (t-self.T[i])/(self.T[i+k]-self.T[i])
                    if self.T[i+k+1] != self.T[i+1]:
                        w2 = (self.T[i+k+1] - t)/(self.T[i+k+1]-self.T[i+1])  
                    self.dpN[i] = w1*self.dpN[i] + w2*self.dpN[i+1]
        
    def __call__(self, t:int):
        c1 = 1e-30 # 分子 float 范围 1.7e38
        c2 = 1e-30 # 分母 float 范围 1.7e38
        self.caculate_dpN(t)
        for i in range(self.n):  # 0 - n
            c1 += self.P[i] * self.dpN[i] * self.W[i]
            c2 += self.dpN[i] * self.W[i]
        return c1/c2
    
    def resample(self, sample:int): 
        if self.P2 != None:
            if self.P2.shape[0] == sample:
                return self.P2 
        first = self.T[self.k]
        gama = 1e-9
        last = self.T[-self.k-1] * (1-gama) + gama * self.T[-self.k-2] # 
        ts = np.linspace(first, last, sample)
        self.P2 = np.array([self(t) for t in ts])
        return self.P2

    def plot(self, sample):
        P2 = self.resample(sample)
        plt.figure(figsize=[6, 8], dpi=96)
        plt.plot(self.P.T[0], self.P.T[1], 'b-.')
        plt.plot(P2.T[0], P2.T[1], color = 'red', marker='+', label='curve')
        plt.scatter(self.P.T[0],self.P.T[1], color='b', label='control points') 
        plt.gca().set_aspect('equal')
        plt.legend()
          
   
        plt.xlabel('x')
        plt.ylabel('y')
        plt.tight_layout() 
 
     

def main():
    P = np.array([[0, 1,],
                [1, 0],
                [0, 2],
                [3, 0],
                [3, 1],
                [2, 4],
                ],) 
    W = [1] * (len(P) - 1) + [1,] 
    p = 3  # 2阶
    n = P.shape[0] - 1  # 从 0 - n, P的长度是n+1
    m = n + p + 1

 

    # clame
    T1 = [0]*p + np.linspace(0, 1, m + 1 - 2 * p).tolist() +  [1] * p  # 从0-m, 个数是 m+1个 
    curve = NURBS.Curve() 
    curve.degree = p
    curve.ctrlpts = P 
    curve.weights = W
    curve.knotvector = T1
    curve.delta = 0.01
    curve.vis = VisMPL.VisCurve2D() 
    curve.render(filename="nurbs-clame.png", plot=False)   
    plt.close()
    c = NURBSCurve(T1, P, curve.weights)
    c.plot(int(1/curve.delta))
    plt.savefig("nurbs-clame(our).png")
    plt.close()


 
    # open
    T2 = np.linspace(0, 1, m + 1).tolist() # 从0-m, 个数是 m+1个 
    curve = NURBS.Curve() 
    curve.degree = p
    curve.ctrlpts = P
    curve.weights = W
    curve.knotvector = T2
    curve.delta = 0.01
    curve.vis = VisMPL.VisCurve2D() 
    curve.render(filename="nurbs-open.png", plot=False)   
    plt.figure(figsize=[10, 8], dpi=96) 
    plt.close()
    c = NURBSCurve(T2, P, curve.weights)
    c.plot(int(1/curve.delta))
    plt.savefig("nurbs-open(our).png")
    plt.close()


 
    # close
    P = np.concatenate((P, P[:p+1]), axis=0)
    W = np.concatenate((W, W[:p+1]), axis=0)
    n = P.shape[0] - 1  # 从 0 - n, P的长度是n+1 
    m = n + p + 1
    T2 = np.linspace(0, 1, m + 1).tolist() # 从0-m, 个数是 m+1个 
    # T2[3] = T2[2]
    curve = NURBS.Curve()
    curve.degree = p
    curve.ctrlpts = P
    curve.weights = W
    curve.knotvector = T2
    curve.delta = 0.01
    curve.vis = VisMPL.VisCurve2D()  
    curve.render(filename="nurbs-close.png", plot=False)  
    plt.close()
    c = NURBSCurve(T2, P, curve.weights)
    c.plot(int(1/curve.delta))
    plt.savefig("nurbs-close(our).png")
    plt.close()


if __name__ == '__main__': 
    main()