Last active
June 2, 2023 11:35
-
-
Save elyase/451cbc00152cb99feac6 to your computer and use it in GitHub Desktop.
2D curve curvature
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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 |
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.
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
doesn't seem to work
[Output:]