Bspline样条曲线和曲面

附件1 B样条曲线代码

import numpy as np
import matplotlib.pyplot as plt


def bspline_basis(i, k, t, knot_vector):
    """计算 B 样条基函数 N_{i,k}(t)"""
    if k == 0:
        return 1.0 if knot_vector[i] <= t <= knot_vector[i + 1] else 0.0
    else:
        denom1 = knot_vector[i + k] - knot_vector[i]
        denom2 = knot_vector[i + k + 1] - knot_vector[i + 1]

        term1 = 0.0
        term2 = 0.0

        if denom1 != 0:
            term1 = (t - knot_vector[i]) / denom1 * bspline_basis(i, k - 1, t, knot_vector)

        if denom2 != 0:
            term2 = (knot_vector[i + k + 1] - t) / denom2 * bspline_basis(i + 1, k - 1, t, knot_vector)

        return term1 + term2


def generate_open_bspline(control_points, degree, num_points=100):
    """生成开放 B 样条曲线"""
    n = len(control_points) - 1  # 控制点数量 - 1
    k = degree  # B 样条的阶数
    m = n + k + 2  # 节点向量的长度

    # 生成开放节点向量
    knot_vector = np.zeros(m)
    knot_vector[k:m - k] = np.linspace(0, 1, m - 2 * k)
    knot_vector[m - k:] = 1
    print(knot_vector)

    # 生成曲线
    curve = np.zeros((num_points, 2))
    t_values = np.linspace(0, 1, num_points)

    for i, t in enumerate(t_values):
        for j in range(n + 1):
            b = bspline_basis(j, k, t, knot_vector)
            curve[i] += b * control_points[j]

    return curve, knot_vector


def main():
    # 示例控制点
    control_points = np.array([
        [0.0, 0.0], [1.0, 2.0], [2.0, 2.0], [3.0, 0.0], [4.0, 3.0], [5.0, 1.0], [8.0, 0.0]
    ])

    degree = 3  # B 样条的阶数

    # 生成开放的 B 样条曲线
    curve, knot_vector = generate_open_bspline(control_points, degree)

    # 绘制 B 样条曲线及其控制点
    plt.figure(figsize=(8, 6))
    plt.plot(curve[:, 0], curve[:, 1], label='B-Spline Curve', color='blue')
    plt.plot(control_points[:, 0], control_points[:, 1], 'ro--', label='Control Points')

    plt.title('Manually Generated Open B-Spline Curve')
    plt.legend()
    plt.grid(True)
    plt.show()


if __name__ == "__main__":
    main()

 运行结果：

附件2 B样条曲面代码

import numpy as np
import matplotlib.pyplot as plt
import utils as utl
from scipy.interpolate import BSpline


# 样条基函数的递归定义
def bspline_basis(i, k, t, knot_vector):
    """计算 B 样条基函数 N_{i,k}(t)"""
    if k == 0:
        return 1.0 if knot_vector[i] <= t <= knot_vector[i + 1] else 0.0
    else:
        denom1 = knot_vector[i + k] - knot_vector[i]
        denom2 = knot_vector[i + k + 1] - knot_vector[i + 1]

        term1 = 0.0
        term2 = 0.0

        if denom1 != 0:
            term1 = (t - knot_vector[i]) / denom1 * bspline_basis(i, k - 1, t, knot_vector)

        if denom2 != 0:
            term2 = (knot_vector[i + k + 1] - t) / denom2 * bspline_basis(i + 1, k - 1, t, knot_vector)

        return term1 + term2


# B样条曲面函数
def bspline_surface(u_values, v_values, control_points, p, q, knot_u, knot_v):
    yn = control_points.shape[0]
    xm = control_points.shape[1]

    # 计算曲面上的每个点
    num_points_u = len(u_values)
    num_points_v = len(v_values)
    surface_points = np.zeros((num_points_u, num_points_v, 3))
    for i, u in enumerate(u_values):
        for j, v in enumerate(v_values):
            sum_pt = np.zeros(3)
            for ui in range(yn):
                for vi in range(xm):
                    N_ui = bspline_basis(ui, p, u, knot_u)
                    M_vi = bspline_basis(vi, q, v, knot_v)
                    sum_pt += N_ui * M_vi * control_points[ui, vi]
                    # print(f'i={i},j={j},N_ui={N_ui},M_vi={M_vi},sum_pt={sum_pt}')
            surface_points[i, j] = sum_pt
            
    return surface_points


def calc_knots(n, k):
    nk = n - 1  # 控制点数量 - 1
    m = nk + k + 2
    knot = np.zeros(m)  # 节点向量的长度
    knot[k:m - k] = np.linspace(0, 1, m - 2 * k)
    knot[m - k:] = 1

    # knot = np.linspace(0, m-1, m)
    return knot


def data_tst(x_max=8, y_max=5, xn=8, yn=8):
    xp = np.linspace(0, x_max, xn)
    yp = np.linspace(0, y_max, yn)
    xq, yq = np.meshgrid(xp, yp)
    zq = np.sin(xq) * np.cos(yq)  # 示例函数值
    pts = np.zeros((yn, xn, 3))
    pts[:, :, 0] = xq
    pts[:, :, 1] = yq
    pts[:, :, 2] = zq
    return pts


if __name__=='__main__':
    # 定义控制点 (3D坐标)
    control_pts = data_tst()
    x_max = np.max(control_pts[:, :, 0])
    y_max = np.max(control_pts[:, :, 1])

    # 定义样条的阶数
    p = 3  # U方向的B样条阶数
    q = 3  # V方向的B样条阶数

    # 定义U和V方向的节点向量， n 是U方向的控制点数减1，U方向节点向量的长度为 n + p + 2
    # m 是V方向的控制点数减1，v方向节点向量的长度为 m + q + 2， p、q 是B样条的阶数
    # 节点向量的首尾重复次数应当等于B样条的阶数加1
    knot_y = calc_knots(control_pts.shape[0], p)
    knot_x = calc_knots(control_pts.shape[1], q)
    knot_x = knot_x * x_max
    knot_y = knot_y * y_max
    print(f'knot_x= {knot_x}')
    print(f'knot_x= {knot_y}')

    # 定义曲面上的点数
    num_points_x = 100
    num_points_y = 100
    # 生成参数 u 和 v 的网格

    y_values = np.linspace(0, y_max, num_points_y)
    x_values = np.linspace(0, x_max, num_points_x)

    surf_pts = bspline_surface(y_values, x_values, control_pts, p, q, knot_y, knot_x)

    # 如果需要，可以将生成的点可视化（例如使用matplotlib）
    utl.plt_show(control_pts[:, :, 2], 1, 'control_pts')
    utl.plt_show(surf_pts[:, :, 2], 1, 'surf_pts')

    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    ax.scatter(control_pts[:, :, 0], control_pts[:, :, 1], control_pts[:, :, 2], color='r', s=5, label='Control Points')
    ax.plot_surface(surf_pts[:, :, 0], surf_pts[:, :, 1], surf_pts[:, :, 2], cmap='viridis', alpha=1.0)

    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_zlabel('Z')
    plt.show()

运行结果：