Skip to content

Instantly share code, notes, and snippets.

@elyase
Last active June 2, 2023 11:35
Show Gist options
  • Star 10 You must be signed in to star a gist
  • Fork 1 You must be signed in to fork a gist
  • Save elyase/451cbc00152cb99feac6 to your computer and use it in GitHub Desktop.
Save elyase/451cbc00152cb99feac6 to your computer and use it in GitHub Desktop.
2D curve curvature
from scipy.interpolate import UnivariateSpline
import numpy as np
def curvature_splines(x, y=None, error=0.1):
"""Calculate the signed curvature of a 2D curve at each point
using interpolating splines.
Parameters
----------
x,y: numpy.array(dtype=float) shape (n_points, )
or
y=None and
x is a numpy.array(dtype=complex) shape (n_points, )
In the second case the curve is represented as a np.array
of complex numbers.
error : float
The admisible error when interpolating the splines
Returns
-------
curvature: numpy.array shape (n_points, )
Note: This is 2-3x slower (1.8 ms for 2000 points) than `curvature_gradient`
but more accurate, especially at the borders.
"""
# handle list of complex case
if y is None:
x, y = x.real, x.imag
t = np.arange(x.shape[0])
std = error * np.ones_like(x)
fx = UnivariateSpline(t, x, k=4, w=1 / np.sqrt(std))
fy = UnivariateSpline(t, y, k=4, w=1 / np.sqrt(std))
xˈ = fx.derivative(1)(t)
xˈˈ = fx.derivative(2)(t)
yˈ = fy.derivative(1)(t)
yˈˈ = fy.derivative(2)(t)
curvature = (xˈ* yˈˈ - yˈ* xˈˈ) / np.power(xˈ** 2 + yˈ** 2, 3 / 2)
return curvature
@lzl200102109
Copy link

Line 39 to 43 is wrong. Please get rid of the prime symbols and make the power 1.5, instead of (3 / 2). Because it will give an integer 1 instead.

@lorenadocarmoinpe
Copy link

Please, is there somebody can explain this code? I am new in Python and I have to calculate the curvature and did not get the idea correctly.

Copy link

ghost commented Oct 10, 2019

doesn't seem to work

x = np.cos(np.arange(20) / 20 * 2 * np.pi)
y = np.sin(np.arange(20) / 20 * 2 * np.pi)
curvature_splines(x, y)

[Output:]

array([0.65568335, 1.18395371, 1.44758424, 1.21313059, 0.93134681,
       0.79791282, 0.79862571, 0.90228408, 1.0675021 , 1.20413488,
       1.20413488, 1.0675021 , 0.90228408, 0.79862571, 0.79791282,
       0.93134681, 1.21313059, 1.44758424, 1.18395371, 0.65568335])

@salem-harvard
Copy link

doesn't seem to work

x = np.cos(np.arange(20) / 20 * 2 * np.pi)
y = np.sin(np.arange(20) / 20 * 2 * np.pi)
curvature_splines(x, y)

[Output:]

array([0.65568335, 1.18395371, 1.44758424, 1.21313059, 0.93134681,
       0.79791282, 0.79862571, 0.90228408, 1.0675021 , 1.20413488,
       1.20413488, 1.0675021 , 0.90228408, 0.79862571, 0.79791282,
       0.93134681, 1.21313059, 1.44758424, 1.18395371, 0.65568335])

It works if you reduce the error to something very small. Like error=0.000000001.

@tjusxh
Copy link

tjusxh commented Nov 4, 2020

Please, what's the formula? I want to get the curvature in the three-dimension. How to do it. Thank you very much.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment