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()

AboAmmar · Answer 1 · 2023-02-11T12:11:42.937

You should first identify the corners; calculate the angle between each segment of the curve, and if the angle is larger than some threshold, you can treat the point as a corner. Then, you can split the curve at the corners, compute the distances between the points, and interpolate each segment individually.

This is how I do it.

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)

# Angle threshold
threshold = np.pi / 6

# Identify the corners
corners = [0, len(x)-1]
for i in range(len(x) - 2):
    dx1, dy1 = x[i+1]-x[i], y[i+1]-y[i]
    dx2, dy2 = x[i+2]-x[i+1], y[i+2]-y[i+1]
    angle = np.arctan2(dy2, dx2) - np.arctan2(dy1, dx1)
    if abs(angle) > threshold:
        corners.append(i+1)
corners.sort()

# Split the curve into segments
segments = [(start, i) for start, i in zip([0] + corners, corners)]

# Interpolate each segment
xinter_list, yinter_list = [], []
for start, end in segments:
    xseg, yseg = x[start:end+1], y[start:end+1]
    
    # Compute the distances
    dx, dy = xseg[+1:]-xseg[:-1],  yseg[+1:]-yseg[:-1]
    ds = np.array((0, *np.sqrt(dx*dx+dy*dy)))

    # Compute the total distance from the 1st point
    s = np.cumsum(ds)
    
    # Interpolate  
    xinter_list.append(np.interp(np.linspace(0,s[-1], 30), s, xseg))
    yinter_list.append(np.interp(np.linspace(0,s[-1], 30), s, yseg))

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

JohanC · Accepted Answer · 2023-02-11T01:57:52.757

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()

thank you for the message. It would work for this four point square example, however I need to be able to distinguish sharp points by the adjacent line segments' angles. Some points in my primary shapes will have no angle. — Dr. Pontchartrain, Feb 11 '23 at 00:16
My apologies. I was using a square to keep it simple. I have updated the example with a polygon. — Dr. Pontchartrain, Feb 11 '23 at 00:29

gboffi · Answer 3 · 2023-02-11T14:25:05.077

Let's say you have a list of "corners", e.g.

corners = 0, 12,43,48

(note that corners comprises the indices of the first and the last points, i.e. the length of x in this example should be 49).

Next you compute sand subsequently a list of corner positions

c_pos = [s[c] for c in corners]

and you can eventually compute the list of points were you want to evaluate the interpolants

ds = (s[-1]-s[0])/npoints # e.g., npoints = 200
arrays = [np.linspace(c_pos0, c_pos1, (c_pos1-c_pos0/ds), endpoint=c_pos1==s[-1])
              for c_pos0, c_pos1 in zip(c_pos[0:], c_pos[1:])
s_new = np.concatenate(arrays)

xinterp = np.interp(s_new, s, x)
...

A few words on

np.linspace(c_pos0, c_pos1, (c_pos1-c_pos0/ds), endpoint=c_pos1==s[-1])

that is meaning: the arrays in the list arrays are "open interval", [c0, c0+Δ,..., c1), except the last one, because we want to draw a closed curve.

This is the easy part of the answer to your question, the hard one is to find the corners...

If someone thinks they have found a robust method to determine the corners, please feel free to incorporate this method in their answers. — gboffi, Feb 11 '23 at 00:27
I've finally got a function that is finding the corners & their indices in a not so elegant way. How may I incorporate your code into it? https://trinket.io/python3/23a153c64e — Dr. Pontchartrain, Feb 11 '23 at 01:17

score 0 · Answer 4 · answered Feb 10 '23 at 23:44

0

It wont be possible to make sharp corners if xinter and yinter steps are not a divider of your side lenght. In this case you would need a stepsize of 25 for example since your side lenght is 100.

Try this:

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

answered Feb 10 '23 at 23:44

tetris programming

809
1
2
13

Thank you for the comment. It's good to know that this is how to get corners with the existing code. I don't mind putting points in addition to the `xinter` & `yinter` if I can identify the corners somehow. – Dr. Pontchartrain Feb 10 '23 at 23:48

Maintaining Sharp Corners in a Numpy Interpolation

4 Answers4

Finding sharp angles

Answer to original question

Linked