Hey there! 作为图像处理新手，你遇到的这个单像素骨架曲线曲率计算问题其实挺典型的——离散阶梯状的曲线直接算梯度确实会导致曲率失真，我给你整理几个实用的平滑方案，搭配标准的曲率计算方法，都是业界常用的思路，上手不难～

一、先理清核心痛点 #

单像素骨架是离散的阶梯状曲线，每个点的邻域只有8个方向可选，导致斜率只能是0、1或无穷大，直接用梯度算曲率自然会得到一堆无意义的数值。解决思路是先把离散的坐标点序列转换成连续平滑的曲线，再计算曲率，而不是直接对图像做卷积（因为图像卷积会模糊骨架的单像素结构）。

二、常用的曲线平滑方法 #

这些方法都是针对骨架的坐标点序列（不是原始二值图像）操作，第一步需要先提取并排序骨架点：

1. 高斯平滑（最常用的基础方法） #

高斯平滑是用高斯核对坐标序列做一维滤波，能有效消除阶梯噪声，同时保留整体弯曲趋势，适合大多数场景。

• 实现步骤：
  1. 提取骨架的(x,y)坐标：用cv2.findNonZero或skimage.measure.find_contours提取所有骨架点，注意要把点转换成有序的序列（不然平滑会混乱）。
  2. 对x、y坐标分别做高斯滤波：用scipy.ndimage.gaussian_filter1d调整平滑程度。
• 代码示例：

import numpy as np
from scipy.ndimage import gaussian_filter1d
import cv2

# 假设skeleton是你的二值骨架数组（0=背景，255=骨架）
skeleton = ... 

# 提取骨架坐标并排序（关键步骤：确保点沿曲线顺序排列）
points = cv2.findNonZero(skeleton).squeeze()  # 形状(N,2)，N为点的数量

# 简单的排序函数（如果提取的点已经有序可跳过）
def sort_skeleton_points(points):
    sorted_points = [points[0]]
    remaining = points[1:].tolist()
    while remaining:
        last_point = sorted_points[-1]
        # 找到剩余点中距离最近的点
        distances = np.linalg.norm(np.array(remaining) - last_point, axis=1)
        closest_idx = np.argmin(distances)
        sorted_points.append(remaining.pop(closest_idx))
    return np.array(sorted_points)

sorted_points = sort_skeleton_points(points)

# 高斯平滑：sigma越大越平滑，建议在1-3之间调整
sigma = 2
x_smoothed = gaussian_filter1d(sorted_points[:, 0], sigma=sigma)
y_smoothed = gaussian_filter1d(sorted_points[:, 1], sigma=sigma)

2. Savitzky-Golay平滑（滑动窗口多项式拟合） #

本质是用滑动窗口内的点拟合低次多项式，既能平滑噪声，又能更好保留局部弯曲特征，适合弯曲细节较多的曲线。

• 代码示例：

from scipy.signal import savgol_filter

# window_length必须是奇数，polyorder是多项式次数（建议2或3）
window_length = 7
polyorder = 2
x_smoothed = savgol_filter(sorted_points[:, 0], window_length, polyorder)
y_smoothed = savgol_filter(sorted_points[:, 1], window_length, polyorder)

3. 三次样条/贝塞尔曲线拟合（连续曲线建模） #

用参数化的连续曲线拟合离散点，得到平滑的曲线表达式，适合整体形状规整的曲线，计算曲率更精准。

• 代码示例：

from scipy.interpolate import splprep, splev

# 三次样条插值：s是平滑因子，s越大越平滑
tck, u_param = splprep([sorted_points[:,0], sorted_points[:,1]], s=10, k=3)
# 生成更多平滑点（比如1000个，方便计算曲率）
u_new = np.linspace(u_param.min(), u_param.max(), 1000)
x_smoothed, y_smoothed = splev(u_new, tck)

三、平滑后计算曲率的标准公式 #

得到平滑的参数化曲线(x(t), y(t))后，用经典的参数化曲线曲率公式计算：
$$
k = \frac{|x'(t)y''(t) - x''(t)y'(t)|}{(x'(t)^2 + y'(t)²⁾{3/2}}
$$
用numpy的梯度函数就能轻松实现：

# 计算一阶导数
dx_dt = np.gradient(x_smoothed)
dy_dt = np.gradient(y_smoothed)
# 计算二阶导数
d2x_dt2 = np.gradient(dx_dt)
d2y_dt2 = np.gradient(dy_dt)
# 计算曲率
curvature = np.abs(dx_dt * d2y_dt2 - d2x_dt2 * dy_dt) / (dx_dt**2 + dy_dt**2)**1.5

高曲率的区域就是你要找的弯曲/弯折处啦～

四、新手注意事项 #

• 排序骨架点是关键：如果点是乱序的，平滑会把不相邻的点混在一起，结果完全错误。
• 调整平滑参数：根据你的曲线噪声程度和细节需求调整sigma、window_length等参数，多试几次找到最优值。
• 可选的预处理：如果骨架有小毛刺，可以先用cv2.morphologyEx做一次小的开运算去除噪声，再提取骨架。

内容的提问来源于stack exchange，提问作者Abhra Basak