Created
January 17, 2022 20:32
-
-
Save kopp/54052bfa9258f3c81dc9d0f7b5d0853a to your computer and use it in GitHub Desktop.
Interpolate between two points with splines
This file contains 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
import numpy as np | |
from scipy.interpolate import interp1d, splprep, splev, CubicHermiteSpline | |
import matplotlib.pyplot as plt | |
pts = np.array( | |
[ | |
[-846724, 0], | |
[-423362, 10029], | |
[0, 13942], | |
[289000, 14733], | |
[547558, 13942], | |
[730000, 11948], | |
[846746, 9015], | |
[846746, 0], | |
[423373, -8311], | |
[0, -12759], | |
[-486000, -14733], | |
[-846724, -12759], | |
[-846724, 0], | |
] | |
) | |
tck, u = splprep(pts.T, u=None, s=0, k=1, per=1) | |
u_new = np.linspace(u.min() + 0.45, u.max(), 1000) | |
x_new, y_new = splev(u_new, tck, der=0) | |
center = np.mean(pts, axis=0) | |
def derivative_at_point(p) -> float: | |
radial_vector = p - center | |
return -(radial_vector[0] / radial_vector[1]) / np.linalg.norm(radial_vector) | |
def derivative_via_neighbors(index) -> float: | |
delta = pts[index + 1] - pts[index - 1] | |
return delta[1] / delta[0] | |
fix, axs = plt.subplots() | |
axs.plot(x_new, y_new, "b-") | |
for i in range(6): | |
spline = CubicHermiteSpline( | |
pts[i : i + 2, 0], | |
pts[i : i + 2, 1], | |
[derivative_at_point(pts[i]), derivative_at_point(pts[i + 1])], | |
) | |
spline_x = np.linspace(*spline.x) | |
spline_y = spline(spline_x) | |
# axs.plot(spline_x, spline_y, "g-") | |
for i in range(5): | |
spline = CubicHermiteSpline( | |
pts[i : i + 2, 0], | |
pts[i : i + 2, 1], | |
[derivative_via_neighbors(i), derivative_via_neighbors(i + 1)], | |
) | |
spline_x = np.linspace(*spline.x) | |
spline_y = spline(spline_x) | |
axs.plot(spline_x, spline_y, "r--") | |
# axs.plot(center[0], center[1], "r+") | |
axs.plot(pts[:, 0], pts[:, 1], "ro") | |
# axs.set_aspect("equal", "box") | |
plt.show() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment