I'm trying to re-implement this paper, which requires me to fit b-spline to a set of 2D points. Later I need to sample the curvature of these points at same arc length intervals. Is there any way to do this with scipy?
I already managed to fit a b-spline to the data. The red crosses indicate the original data and the curvature is indicated by the color of the curve. The blue dots should be at equal arc length intervals which is not working at the moment.
I know that one problem is that I compute the arc length only as an approximation, which is really bad in my case, since I only have a couple of points. One idea I had was to first fit the discrete points with a B-spline and then extract the arc length at many more points. I could then create a new B-spline representation with this info. Does this work? It seems unnecessarily complicated.
import numpy as np
import scipy as sp
import scipy.interpolate
import matplotlib.pyplot as plt
f, (ax1, ax2) = plt.subplots(1, 2, figsize=(20, 10))
# x,y coordinates of contour points, not monotonically increasing
x = np.array([87., 98., 100., 95., 100., 108., 110., 118., 120., 117., 105., 100., 92., 90.])
y = np.array([42., 35., 32., 25., 18., 20., 27., 27., 35., 46., 45., 48., 55., 51. ])
# build a spline representation of the contour
# Summed euclidian distance is good approximation https://stackoverflow.com/questions/36795447/length-measurement-along-a-b-spline
squared_diff_x, squared_diff_y = np.diff([x, y]) ** 2
arclen = np.sqrt(squared_diff_x + squared_diff_y)
arclen_cum = np.concatenate(([0], arclen.cumsum()))
spline, u = sp.interpolate.splprep([x, y], u=arclen_cum, s=0)
eval_points = np.linspace(0, arclen_cum[-1], 500)
interp = sp.interpolate.splev(eval_points, spline)
derv2 = sp.interpolate.splev(eval_points, spline, der=2)
# https://en.wikipedia.org/wiki/Curvature
curv = np.sqrt(derv2[0]**2 + derv2[1]**2)
ax1.scatter(interp[0], interp[1], c=curv, cmap='magma')
# Check how sampling at equal arc length intervals looks like
a, b = sp.interpolate.splev(np.linspace(0, arclen_cum[-1], 50), spline)
ax1.scatter(a, b, s=100)
ax1.scatter(x, y, s=200, color='r', marker='X')
ax2.plot(curv)
plt.title('Curvature of function')
plt.xlabel('x')
plt.ylabel('Curvature')