Maintaining Sharp Corners in a Numpy Interpolation

Question

I am interpolating a shape with numpy's linspace and interp using gboffi's magnificent code from this post (included, below).

This works well, however, the corners sometimes get missed and the resulting softened shape is undesired.

I'd like to maintain the sharp corners of my shapes with an angle threshold parameter. Is there any way to keep the corners of the shape's interpolation if the next line segment is of a sharp enough angle? Thank you!

from matplotlib import pyplot as plt
import numpy as np

x = np.array([815.9, 693.2, 570.4, 462.4, 354.4, 469.9, 585.4, 700.6, 815.9])
y = np.array([529.9, 637.9, 746, 623.2, 500.5, 326.9, 153.3, 341.6, 529.9])

fig, ax = plt.subplots(1)
ax.set_aspect('equal')
ax.scatter(x, y, s=40, zorder=3, alpha=0.3)

# compute the distances, ds, between points
dx, dy = x[+1:]-x[:-1],  y[+1:]-y[:-1]
ds = np.array((0, *np.sqrt(dx*dx+dy*dy)))

# compute the total distance from the 1st point, measured on the curve
s = np.cumsum(ds)

# interpolate
xinter = np.interp(np.linspace(0,s[-1], 30), s, x)
yinter = np.interp(np.linspace(0,s[-1], 30), s, y)

# plot the interpolated points
ax.plot(xinter, yinter)
ax.scatter(xinter, yinter, s=5, zorder=4)
plt.show()

Answer 1

Finding sharp angles

To find all point where the angle is larger than some cut-off, you could calculate the arc cosines of the dot products of the normalized difference vectors.

# cosines of the angles are the normalized dotproduct of the difference vectors, ignore first and last point
cosines = np.array( [1, *((dx[:-1]*dx[1:] +  dy[:-1]*dy[1:]) / ds[1:-1] /ds[2:]),1])

An alternative formula is the atan2 of the cross product and the dot product, which avoids divisions.

from matplotlib import pyplot as plt
import numpy as np

x = np.array([815.9, 693.2, 570.4, 462.4, 354.4, 469.9, 585.4, 700.6, 815.9])
y = np.array([529.9, 637.9, 746, 623.2, 500.5, 326.9, 153.3, 341.6, 529.9])

fig, ax = plt.subplots(1)
ax.set_aspect('equal')
ax.scatter(x, y, s=40, zorder=3, alpha=0.3)

# compute the distances, ds, between points
dx, dy = x[+1:] - x[:-1], y[+1:] - y[:-1]
ds = np.append(0, np.sqrt(dx * dx + dy * dy))

# compute the total distance from the 1st point, measured on the curve
s = np.cumsum(ds)

# angle: atan2 of cross product and dot product
cut_off_angle = 10  # angles larger than this will be considered sharp
angles_rad = np.arctan2(dx[:-1] * dy[1:] - dy[:-1] * dx[1:],
                        dx[:-1] * dx[1:] + dy[:-1] * dy[1:])
# convert to degrees, and pad with zeros for first and last point
angles = np.pad(np.degrees(angles_rad), 1)

# interpolates
s_long = np.sort(np.append(np.linspace(0, s[-1], 30),
                           s[np.abs(angles) > cut_off_angle]))
xinter = np.interp(s_long, s, x)
yinter = np.interp(s_long, s, y)

# plot the interpolated points
ax.plot(xinter, yinter)
ax.scatter(xinter, yinter, s=5, zorder=4)
plt.show()

Answer to original question

A simple extension is to insert the points of s into the array used for the intermediate points (np.linspace(0,s[-1], 30)).

from matplotlib import pyplot as plt
import numpy as np

x = np.array([100, 200, 200, 100, 100])
y = np.array([100, 100, 200, 200, 100])

fig, ax = plt.subplots(1)
ax.set_aspect('equal')
ax.scatter(x, y, s=40, zorder=3, alpha=0.3)

# compute the distances, ds, between points
dx, dy = x[+1:] - x[:-1], y[+1:] - y[:-1]
ds = np.array((0, *np.sqrt(dx * dx + dy * dy)))

# compute the total distance from the 1st point, measured on the curve
s = np.cumsum(ds)

# interpolates
s_long = np.sort(np.append(np.linspace(0, s[-1], 30)[1:-1], s))
xinter = np.interp(s_long, s, x)
yinter = np.interp(s_long, s, y)

# plot the interpolated points
ax.plot(xinter, yinter)
ax.scatter(xinter, yinter, s=5, zorder=4)
plt.show()