用numpy计算曲率时出错 [英] Error calculating curvature with numpy

问题描述

我正在尝试使用此处的公式.我遇到的问题是，尽管我获得了应有的恒定值，但该值不是正确的值.这是我的代码:

I am trying to calculate the curvature of a 2D curve at each point using the formula here. The problem I am having is that while I am getting a constant value as it should be, this value is not the correct one. Here is my code:

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点不是均匀分布的.

I was expecting to get something close to a 5. Can somebody help me spot the error or suggest a more robust way to do this? In practice my x,y points are not evenly spaced.

我使用gaussian_filter1d而不是np.gradient作为衍生工具，因为在此处是一种更可靠的方法，尤其是对于二阶导数.

I am using gaussian_filter1d instead of np.gradient for the derivative because it was shown here that this is a more robust method, especially for the second derivative.

推荐答案

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

The formula for curvature depends on the first and second derivatives of x and y.

您的代码假定gaussian_filter1d与x的一阶导数相同.不是.

Your code is assuming that the gaussian_filter1d is the same as the first derivative of x. It is not.

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

编辑如果您想通过gaussian_filter1d，应该没问题，但是二阶导数的计算没有达到您的期望.这是一些工作代码，在这里我完成了两个一阶导数以获得x2和y2:

edit If you want to go through gaussian_filter1d, you should be alright but the calculation of the second derivative is not doing what you expect. Here is some working code where I've done 2 first derivatives to get x2 and 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的过滤器却无法工作.

After a lot of careful checking I saw that y2 wasn't looking much like -y and similarly for x2. The change I made from your code is that now y2 and x2 come from y1 and x1 with gaussian_filter1d having order=1. I don't know enough about what the filter is doing to be able to say why two passes through the filter with order=1 seems to work but a single pass with order=2 doesn't.

这篇关于用numpy计算曲率时出错的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！