CODEWITHSUNDEEP
pythonscipyinterpolation
megaprimatus
megaprimatus

Reputation: 83

scipy.interpolate smooth closed curve

I have this code, which gives me this closed line thru sets of points. When I change the degree k=3 Im expecting to see the same line but more like ellipse shaped and not these bumpy curves like shown below in a second image. Is there any other settings needs to be put into "splprep". Thanks

import numpy as np
from scipy.interpolate import interp1d, splprep, splev
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(), u.max(), 1000)
x_new, y_new = splev(u_new, tck, der=0)

plt.plot(pts[:,0], pts[:,1], 'ro')
plt.plot(x_new, y_new, 'b--')
plt.show()

enter image description here

enter image description here

Upvotes: 0

Views: 1079

Answers (1)

Michael Kopp
Michael Kopp

Reputation: 1711

What you could do, is, do do an interpolation between two consecutive points. This makes sure, that the final curve goes through all points. You can then additionally constraint the derivative of this interpolation. This would make sure, that two adjacent segments are "smoothly" connected to each other. Taking your example and using CubicHermiteSpline could for example result in this smooth curve (red) in comparison to the piecewise linear (blue) interpolation:

enter image description here

There is some degree of freedom, in how to select the constraint for the derivative. In the example I've decided that for point i, the derivative should be similar to the slope of the vector connecting points i-1 and i+1. You could chose other conditions here, too.

The code to represent the figure is

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)
# 0.45 is chosen to get roughly the part not covered by the piecewise spline demo
u_new = np.linspace(u.min() + 0.45, u.max(), 1000)
x_new, y_new = splev(u_new, tck, der=0)


fix, axs = plt.subplots()

axs.plot(x_new, y_new, "b-")


def derivative_via_neighbors(index) -> float:
    delta = pts[index + 1] - pts[index - 1]
    return delta[1] / delta[0]


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(pts[:, 0], pts[:, 1], "ro")
# axs.set_aspect("equal", "box")
plt.show()

Note, that you need to take care to handle the points at y=0 (where the derivative might become infinity). Those cases are not handled by the code above (that's why `range(5) is chosen conveniently for me).

Upvotes: 2

Related Questions