准均匀B样条拟合3d点

参考：B样条曲线讲解视频_哔哩哔哩_bilibili

一些概念

• 待拟合的点称之为样本点，也是控制点，控制点组成的多边形称之为控制多边形。
• 基函数不为0的区间称之为支撑区间。

b样条相比于贝塞尔曲线优点

• 可以局部拟合，不会牵一发而动全身。这是由其基函数的支撑区间是局部的决定的，贝塞尔基函数的支撑区间是全局的。
• 支撑区间和b样条阶次有关，阶次越大，支撑区间越大。平滑效果越好，但更容易牵一发而动全身。
• 由于支持局部拟合，因此可以对控制多边形的逼近程度更好。

b样条拟合性质

• 如果有n个待拟合的点，想用k次B样拟合，则节点knots的数量为：n+k+1
• 准均匀knots的分布是

  k+1个0

  n-1-k个 (0，1）均匀分布

  k+1个1

    k <= n-1：可以从节点表中间可以看出，k>n-1后，中间部分大小为负数。
• 基函数的生成与控制点的数量有关系，与控制点具体的内容（取值大小）无关。
• 最终b样条的拟合点 是控制点（样本点）的加权，权值来自与基函数。
  • 以3阶样条为例子，t从0到1，t = t0处，与几个基函数会对应的几个交点解释：
  • 
    • 是对应样本点的权值Fi,k(t)，对样本点加权求和就是b样条拟合点。
      • 对应理论公式
      • 
    • 几个交点值的求和=1，也就是权值和为1
    • 代码里，Hnm中，每一列就是一个基函数。Hnm（m是控制点数量） * points ，就是拟合的点结果。注意，b样条基函数对控制点加权时，是对各个维度独立加权，因此把3d点的拟合代码修改成2d的，不需要修改基函数Hnm，只需要修改points即可。
  • 一个测试，6个蓝点是控制点，红线是3次b样条上的100个拟合点。
  • 

import numpy as np
import matplotlib.pyplot as plt

#  基于节点表knots, 第i个k阶基函数，在t处的取值。
def basefunction(t, i, k, knots):
    # i = base_funtion_idx
    # k是权重基函数的次数。
    # t 从 0-1 取值，是拟合曲线上的采样间隔。
    # knots：生成控制点权重基函数的节点系数。
    if k==0:
        return 1.0 if t>=knots[i] and t<knots[i+1] else 0
    
    if (t<knots[i] or t>knots[i+k+1]):
        return 0
    else:
        denominator1 = knots[i+k] - knots[i]
        coeff1 = 0 if denominator1 == 0 else (t-knots[i]) / denominator1
        denominator2 = knots[i+k+1] - knots[i+1]
        coeff2 = 0 if denominator2 == 0 else (knots[i+k+1] - t) / denominator2
        return coeff1*basefunction(t, i, k-1, knots) +\
              coeff2 * basefunction(t, i+1, k-1, knots)
    
# 每一列为一个基函数，共有control_point_number个。
def GetHnm(k, control_point_number, fitpoints_number):
    # n = fitpoints_number
    # m = 控制点 = control_point_number
    knots = [0.0] * k + list(np.linspace(0, 1, (control_point_number - k + 1))) + [1.0] *k
   
    basefunction_idx = [i for i in range(0, len(knots) - k -1)]
    ts = np.linspace(0, 1, fitpoints_number)
    Hnm = []
    
    for i in basefunction_idx:
        basefunction_i_value = []
        for t in ts:
            value = basefunction(t, i, k, knots)
            basefunction_i_value.append(value)
        # 最后一个base function的最后一位设置为1
        if  i == basefunction_idx[-1]:
            basefunction_i_value[-1] = 1
        Hnm.append(basefunction_i_value)
       
    return np.transpose(np.array(Hnm))


def drawbasefunction(k, control_point_number, fitpoints_number):
    for i in range(0, k):
        Hnm = GetHnm(i, control_point_number, fitpoints_number)
        ts = np.linspace(0, 1, fitpoints_number)
        for col in range(0, Hnm.shape[1]):
            ys = Hnm[:, col]
            plt.subplot(k, 1, i+1)
            plt.plot(ts, ys)
    plt.show()


def drawfit(origin_points, fit_points):
    ax = plt.axes(projection='3d')
    ax.scatter(origin_points[:, 0], origin_points[:, 1], origin_points[:, 2], label = "origin points")
    ax.plot(fit_points[:, 0], fit_points[:, 1], fit_points[:, 2], label = "fit points", color = "red")
    plt.legend()
    plt.show()


points = np.array([[1, 1.1, 1.2], [2, 2.1, 1.9], [1.9, 2.9, 2.9],  [3.9, 4, 4.2], [5.9, 6.3,6], [3.9, 6.6, 6.8]])
k = 3
fitpoints_number = 100

drawbasefunction(k, len(points), 100)

Hnm = GetHnm(k, len(points), fitpoints_number)
fit_points = np.matmul(Hnm, points)

drawfit(points, fit_points)