How do I smooth out the edges of a closed line similar to d3's curveCardinal method implementation?

Question

I have a few data points that I am connecting using a closed line plot, and I want the line to have smooth edges similar to how the curveCardinal methods in d3 do it. Link Here

Here's a minimal example of what I want to do:

import numpy as np
from matplotlib import pyplot as plt

x = np.array([0.5, 0.13, 0.4, 0.5, 0.6, 0.7, 0.5])
y = np.array([1.0, 0.7, 0.5, 0.2, 0.4, 0.6, 1.0])

fig, ax = plt.subplots()

ax.plot(x, y)
ax.scatter(x, y)

Now, I'd like to smooth out/interpolate the line similar to d3's curveCardinal methods. Here are a few things that I've tried.

from scipy import interpolate

tck, u = interpolate.splprep([x, y], s=0, per=True)  
xi, yi = interpolate.splev(np.linspace(0, 1, 100), tck)

fig, ax = plt.subplots(1, 1)

ax.plot(xi, yi, '-b')
ax.plot(x, y, 'k')
ax.scatter(x[:2], y[:2], s=200)
ax.scatter(x, y)

The result of the above code is not bad, but I was hoping that the curve would stay closer to the line when the data points are far apart (I increased the size of two such data points above to highlight this). Essentially, have the curve stay close to the line.

Using interp1d (has the same problem as the code above):

from scipy.interpolate import interp1d

x = [0.5, 0.13, 0.4, 0.5, 0.6, 0.7, 0.5]
y = [1.0, 0.7, 0.5, 0.2, 0.4, 0.6, 1.0]

orig_len = len(x)
x = x[-3:-1] + x + x[1:3]
y = y[-3:-1] + y + y[1:3]

t = np.arange(len(x))
ti = np.linspace(2, orig_len + 1, 10 * orig_len)

kind='cubic'
xi = interp1d(t, x, kind=kind)(ti)
yi = interp1d(t, y, kind=kind)(ti)

fig, ax = plt.subplots()
ax.plot(xi, yi, 'g')
ax.plot(x, y, 'k')
ax.scatter(x, y)

I also looked at the Chaikins Corner Cutting algorithm, but I don't like the result.

def chaikins_corner_cutting(coords, refinements=5):
    coords = np.array(coords)

    for _ in range(refinements):
        L = coords.repeat(2, axis=0)
        R = np.empty_like(L)
        R[0] = L[0]
        R[2::2] = L[1:-1:2]
        R[1:-1:2] = L[2::2]
        R[-1] = L[-1]
        coords = L * 0.75 + R * 0.25

    return coords


fig, ax = plt.subplots()

ax.plot(x, y, 'k', linewidth=1)
ax.plot(chaikins_corner_cutting(x, 4), chaikins_corner_cutting(y, 4))

I also, superficially, looked at Bezier curves, matplotlibs PathPatch, and Fancy box implementations, but I couldn't get any satisfactory results.

Suggestions are greatly appreciated.

Answer 1

So, here's how I ended up doing it. I decided to introduce new points between every two existing data points. The following image shows how I am adding these new points. Red are data that I have. Using a convex hull I calculate the geometric center of the data points and draw lines to it from each point (shown with blue lines). Divide these lines twice in half and connect the resulting points (green line). The center of the green line is the new point added.

Here are the functions that accomplish this:

def midpoint(p1, p2, sf=1):
    """Calculate the midpoint, with an optional 
    scaling-factor (sf)"""
    xm = ((p1[0]+p2[0])/2) * sf
    ym = ((p1[1]+p2[1])/2) * sf
    return (xm, ym)

def star_curv(old_x, old_y):
    """ Interpolates every point by a star-shaped curve. It does so by adding
    "fake" data points in-between every two data points, and pushes these "fake"
    points towards the center of the graph (roughly 1/4 of the way).
    """

    try:
        points = np.array([old_x, old_y]).reshape(7, 2)
        hull = ConvexHull(points)
        x_mid = np.mean(hull.points[hull.vertices,0])
        y_mid = np.mean(hull.points[hull.vertices,1])
    except:
        x_mid = 0.5
        y_mid = 0.5

    c=1
    x, y = [], []
    for i, j in zip(old_x, old_y):
        x.append(i)
        y.append(j)
        try:
            xm_i, ym_i = midpoint((i, j),
                                  midpoint((i, j), (x_mid, y_mid)))

            xm_j, ym_j = midpoint((old_x[c], old_y[c]),
                                  midpoint((old_x[c], old_y[c]), (x_mid, y_mid)))

            xm, ym = midpoint((xm_i, ym_i), (xm_j, ym_j))
            x.append(xm)
            y.append(ym)
            c += 1
        except IndexError:
            break


    orig_len = len(x)
    x = x[-3:-1] + x + x[1:3]
    y = y[-3:-1] + y + y[1:3]


    t = np.arange(len(x))
    ti = np.linspace(2, orig_len + 1, 10 * orig_len)

    kind='quadratic'
    xi = interp1d(t, x, kind=kind)(ti)
    yi = interp1d(t, y, kind=kind)(ti)

    return xi, yi

Here's how it looks:

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d
from scipy.spatial import ConvexHull

x = [0.5, 0.13, 0.4, 0.5, 0.6, 0.7, 0.5]
y = [1.0, 0.7, 0.5, 0.2, 0.4, 0.6, 1.0]

xi, yi = star_curv(x, y)

fig, ax = plt.subplots()
ax.plot(xi, yi, 'g')

ax.plot(x, y, 'k', alpha=0.5)
ax.scatter(x, y, color='r')

The result is especially noticeable when the data points are more symmetric, for example the following x, y values give the results in the image below:

x = [0.5, 0.32, 0.34, 0.5, 0.66, 0.65, 0.5]
y = [0.71, 0.6, 0.41, 0.3, 0.41, 0.59, 0.71]

Comparison between the interpolation presented here, with the default interp1d interpolation.

Answer 2

I would create another array with the vertices extended in/out or up/down by about 5%. So if a point is lower than the average of the neighbouring points, make it a bit lower still.

Then do a linear interpolation between the new points, say 10 points/edge. Finally do a spline between the second last point per edge and the actual vertex. If you use Bezier curves, you can make the spline come in at the same angle on each side.

It's a bit of work, but of course you can use this anywhere.