Error calculating curvature with numpy

Question

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,
         ...

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.

EDIT: 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.

Answer 1

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

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

Look up np.gradient(x,dalpha) where dalpha is the step size.

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)

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.