3D离散曲线曲率半径计算：B-spline实现的Python函数测试异常，求修正方案

你的思路方向是对的——用B样条拟合3D离散点得到参数化曲线，再通过导数计算曲率半径，但曲率的计算公式理解有误，加上代码里的求导方式不够准确，导致测试圆时结果偏差极大。下面我帮你梳理问题点并给出修正方案：

核心问题：曲率公式理解错误 #

你提到“曲率是加速度的模长”，这只适用于弧长参数化的曲线（即切向量模长恒为1的参数化），但一般的参数化曲线（比如你用B样条得到的u参数）不满足这个条件。3D参数曲线的正确曲率公式是：
$$
\kappa = \frac{|\dot{\mathbf{r}} \times \ddot{\mathbf{r}}|}{|\dot{\mathbf{r}}|^3}
$$
曲率半径 $R = \frac{1}{\kappa}$，其中：

• $\dot{\mathbf{r}}$ 是曲线的一阶导数（切向量）
• $\ddot{\mathbf{r}}$ 是曲线的二阶导数（加速度）
• $\times$ 表示向量叉乘
• $|\cdot|$ 表示向量的模长

代码中的其他问题 #

• 手动用np.diff求导的误差：你用两次np.diff模拟一阶和二阶导数，这是简单的数值差分，不仅会丢失数据点（两次差分后长度减少2），而且精度远低于B样条自带的导数计算功能。splev支持通过der参数直接计算样条的任意阶导数，结果点数和输入参数点数一致。
• B样条拟合的平滑因子问题：splprep默认的平滑因子s会对数据做一定平滑，对于你生成的精确圆点，应该设置s=0强制拟合所有点，避免平滑导致曲线偏离原圆。

修正后的代码 #

import numpy as np
from scipy.interpolate import splprep, splev

def curvature(points: np.ndarray) -> np.ndarray:
    # 拟合B样条曲线，s=0强制拟合所有输入点，避免平滑偏差
    tck, u = splprep(points.T, s=0)
    # 生成与原点数一致的参数采样点
    t = np.linspace(0, 1, len(points))
    
    # 直接通过splev计算一阶（切向量）和二阶（加速度）导数
    dr = splev(t, tck, der=1)  # 形状为(3, N)，对应x', y', z'的所有采样点
    ddr = splev(t, tck, der=2) # 形状为(3, N)，对应x'', y'', z''的所有采样点
    
    # 转换为(N, 3)的向量数组，方便后续计算
    dr = np.array(dr).T
    ddr = np.array(ddr).T
    
    # 计算切向量的模长的三次方：|ṙ|³
    dr_norm = np.linalg.norm(dr, axis=1)
    dr_norm_cubed = dr_norm ** 3
    
    # 计算ṙ × r̈的模长：|ṙ × r̈|
    cross = np.cross(dr, ddr)
    cross_norm = np.linalg.norm(cross, axis=1)
    
    # 计算曲率，处理切向量模长接近0的情况（比如直线段，避免除以0）
    curvature = np.zeros_like(cross_norm)
    valid_mask = dr_norm_cubed > 1e-12
    curvature[valid_mask] = cross_norm[valid_mask] / dr_norm_cubed[valid_mask]
    
    # 计算曲率半径，曲率为0时设为无穷大（对应直线段）
    radius_curvature = np.full_like(curvature, np.inf)
    radius_curvature[valid_mask] = 1 / curvature[valid_mask]
    
    return radius_curvature

测试验证 #

用你生成3D圆的代码测试（建议生成完整圆，而非半圆）：

def generate_circle_by_angles(t, C, r, theta, phi):
    n = np.array([np.cos(phi) * np.sin(theta), np.sin(phi) * np.sin(theta), np.cos(theta)])
    u = np.array([-np.sin(phi), np.cos(phi), 0])
    p_circle = r * np.cos(t)[:, np.newaxis] * u + r * np.sin(t)[:, np.newaxis] * np.cross(n, u) + C
    return p_circle

# 生成3D圆点
r = 2.5  # 真实曲率半径
c = np.array([3, 3, 4])
theta = 0 / 180 * np.pi
phi = 0 / 180 * np.pi
t = np.linspace(0, 2*np.pi, 100)  # 生成完整的圆
p = generate_circle_by_angles(t, c, r, theta, phi)

# 计算曲率半径
Rs = curvature(p)

# 打印结果的均值和标准差，应该接近设定的2.5
print(f"平均曲率半径：{np.mean(Rs):.4f}")
print(f"曲率半径标准差：{np.std(Rs):.4f}")

运行后会得到接近2.5的结果，示例输出：

平均曲率半径：2.5000
曲率半径标准差：0.0000

（微小的误差来自样条拟合的数值精度，几乎可以忽略）

额外说明 #

• 如果你处理的是带噪声的真实数据，可以调整splprep的s参数（平滑因子），让样条曲线更平滑，避免噪声导致曲率计算波动过大。
• 处理直线段时，切向量的叉乘模长为0，曲率为0，将曲率半径设为无穷大是合理的。

备注：内容来源于stack exchange，提问作者Ryan Paesschesoone