python - 使用 numpy 计算曲率时出错

Question

我正在尝试使用 formula here 计算二维曲线在每个点的曲率.我遇到的问题是，虽然我得到了一个应该的常量值，但这个值不正确。这是我的代码:

from scipy.ndimage import gaussian_filter1d
import numpy as np

def curvature(x, y):
    #first and second derivative
    x1 = gaussian_filter1d(x, sigma=1, order=1, mode='wrap')
    x2 = gaussian_filter1d(x, sigma=1, order=2, mode='wrap')
    y1 = gaussian_filter1d(y, sigma=1, order=1, mode='wrap')
    y2 = gaussian_filter1d(y, sigma=1, order=2, mode='wrap')
    return np.abs(x1*y2 - y1*x2) / np.power(x1**2 + y1**2, 3/2)

# make circle data
alpha = np.linspace(-np.pi/2,np.pi/2, 1000)
R = 5
x = R*np.cos(alpha)
y = R*np.sin(alpha)

>>> 1 / curvature(x, y)
array([  9.60e+02,   5.65e+01,   4.56e-01,   1.41e-02,   6.04e-01,
         6.04e-01,   6.04e-01,   6.04e-01,   6.04e-01,   6.04e-01,
         6.04e-01,   6.04e-01,   6.04e-01,   6.04e-01,   6.04e-01,
         ...

我原本希望得到接近 5 的分数。有人可以帮我找出错误或建议更可靠的方法吗？实际上，我的 x,y 点不是均匀分布的。

编辑:我使用 gaussian_filter1d 而不是 np.gradient 作为导数，因为它显示为 here这是一种更稳健的方法，尤其是对于二阶导数。

Answer 1

曲率公式取决于 x 和 y 的一阶和二阶导数。

您的代码假设 gaussian_filter1d 与 x 的一阶导数相同。它不是。

查找 np.gradient(x,dalpha)，其中 dalpha 是步长。

编辑 如果你想通过 gaussian_filter1d，你应该没问题，但二阶导数的计算没有达到你的预期。这是一些工作代码，我已经完成了 2 个一阶导数以获得 x2 和 y2:

import numpy as np
def curvature(x, y):
    #first and second derivative
    dalpha = np.pi/1000
    x1 = gaussian_filter1d(x, sigma=1, order=1, mode='wrap')
    x2 = gaussian_filter1d(x1, sigma=1, order=1, mode='wrap')
    y1 = gaussian_filter1d(y, sigma=1, order=1, mode='wrap')
    y2 = gaussian_filter1d(y1, sigma=1, order=1, mode='wrap')
    return np.abs(x1*y2 - y1*x2) / np.power(x1**2 + y1**2, 3./2)

# make circle data
alpha = np.linspace(-np.pi/2,np.pi/2, 1000)
R = 5
x = R*np.cos(alpha)
y = R*np.sin(alpha)

print 1/curvature(x,y)

经过大量仔细检查后，我发现 y2 看起来不太像 -y，x2 也很相似。我对您的代码所做的更改是，现在 y2 和 x2 来自 y1 和 x1 以及 gaussian_filter1d 具有 order=1。我不太了解过滤器在做什么，无法说明为什么使用 order=1 两次通过过滤器似乎有效，但一次通过 order=2 没有。