基于骨架提取与傅里叶描述符的曲线相似度计算

一、问题背景

在计算机视觉领域，曲线相似度计算是形状匹配、手写识别、运动轨迹分析等任务的核心技术。本文介绍通过骨架提取结合傅里叶描述符的方法，实现两条曲线的形状相似度量化。

二、方法流程

1. 骨架提取

def skeletonize_image(binary_image):
    skeleton = skeletonize(binary_image // 255)
    return (skeleton * 255).astype(np.uint8)

• 使用skimage的骨架化算法

• 将二值图像转换为单像素宽度的骨架

• 优势：消除线宽差异对比较的影响

2. 轮廓提取与处理

contours, _ = cv2.findContours(...)
filtered = sorted(contours, key=len, reverse=True)

• 使用OpenCV提取骨架线段

• 过滤短线段并按长度降序排列

• 取最长两条曲线作为比较对象（示例）

3. 曲线对齐（PCA）

centroid = np.mean(curve, axis=0)
rotation_matrix = eigvecs[...]
aligned_curve = curve_centered @ rotation_matrix

• 通过主成分分析对齐曲线方向

• 消除空间位置和旋转带来的差异

4. 曲线重采样

interp_x = interp1d(...)
interp_y = interp1d(...)
resampled_curve = np.vstack(...)

• 统一曲线采样点数量（默认100点）

• 保证不同长度曲线的可比性

5. 傅里叶描述符

complex_curve = curve[:,0] + 1j*curve[:,1]
dft = np.fft.fft(complex_curve)
magnitude = np.abs(dft[:10])

• 将曲线转换为复数序列

• 取前10个傅里叶系数的模值

• 特征优势：平移/旋转不变性

6.相似度计算

fourier_similarity = np.linalg.norm(des1 - des2)

• 计算傅里叶描述符的欧氏距离
• 值越小表示相似度越高

三、完整示例

"""
曲线相似度计算核心算法
环境要求：Python 3.8+，numpy, opencv, scikit-image, scipy
"""
import numpy as np
import cv2
from scipy.interpolate import interp1d
from skimage.morphology import skeletonize

def skeletonize_image(binary_img):
    """二值图像骨架提取"""
    skeleton = skeletonize(binary_img // 255)
    return (skeleton * 255).astype(np.uint8)

def process_contours(skeleton_img, min_length=20):
    """轮廓提取与处理"""
    contours, _ = cv2.findContours(skeleton_img, 
                                 cv2.RETR_EXTERNAL,
                                 cv2.CHAIN_APPROX_NONE)
    contours = [c[:,0,:] for c in contours if len(c) >= min_length]
    return sorted(contours, key=len, reverse=True)

def pca_alignment(curve):
    """基于主成分分析的曲线对齐"""
    centroid = np.mean(curve, axis=0)
    cov = np.cov((curve - centroid).T)
    _, eigvecs = np.linalg.eig(cov)
    return (curve - centroid) @ eigvecs

def uniform_resample(curve, n_points=100):
    """均匀重采样"""
    dist = np.cumsum(np.sqrt(np.sum(np.diff(curve, axis=0)**2, axis=1)))
    dist = np.insert(dist, 0, 0)
    return np.column_stack([
        interp1d(dist, curve[:,0])(np.linspace(0, dist[-1], n_points)),
        interp1d(dist, curve[:,1])(np.linspace(0, dist[-1], n_points))
    ])

def fourier_descriptor(curve, n_coeff=10):
    """傅里叶描述符计算"""
    complex_curve = curve[:,0] + 1j*curve[:,1]
    return np.abs(np.fft.fft(complex_curve)[:n_coeff])

def curve_similarity(img_path, n_points=100, n_coeff=10):
    """完整计算流程"""
    # 图像预处理
    img = cv2.imread(img_path, cv2.IMREAD_GRAYSCALE)
    _, binary = cv2.threshold(img, 127, 255, cv2.THRESH_BINARY)
    
    # 骨架提取
    skeleton = skeletonize_image(binary)
    
    # 轮廓处理
    contours = process_contours(skeleton)
    if len(contours) < 2: 
        raise ValueError("需至少检测到两条曲线")
    
    # 处理前两条曲线
    curve1 = uniform_resample(pca_alignment(contours[0]), n_points)
    curve2 = uniform_resample(pca_alignment(contours[1]), n_points)
    
    # 计算相似度
    fd1 = fourier_descriptor(curve1, n_coeff)
    fd2 = fourier_descriptor(curve2, n_coeff)
    return np.linalg.norm(fd1 - fd2)

if __name__ == "__main__":
    # 示例用法
    similarity = curve_similarity("test.png")
    print(f"曲线相似度指标: {similarity:.2f}")