Ошибка вычисления кривизны с numpy

Question

Я пытаюсь вычислить кривизну 2D-кривой в каждой точке, используя 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 не равномерно распределены.

EDIT: Я использую 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 этого не делает.