4

I'm trying to emulate Excel's

Insert>Scatter>Scatter with smooth lines and markers

command in Matplotlib

The scipy function interpolate creates a similar effect, with some nice examples of how to simply implement this here: How to draw cubic spline in matplotlib

However Excel's spline algorithm is also able to generate a smooth curve through just three points (e.g. x = [0,1,2] y = [4,2,1]); and it isn't possible to do this with cubic splines.

I have seen discussions that suggest that the Excel algorithm uses Catmull-Rom splines; but don't really understand these, or how they could be adapted to Matplotlib: http://answers.microsoft.com/en-us/office/forum/office_2007-excel/how-does-excel-plot-smooth-curves/c751e8ff-9f99-4ac7-a74a-fba41ac80300

Is there a simple way of modifying the above examples to achieve smooth curves through three or more points using the interpolate library?

Many thanks

Community
  • 1
  • 1
Dave
  • 515
  • 1
  • 8
  • 17

2 Answers2

6

By now you may have found the Wikipedia page for the Centripetal Catmull-Rom spline, but in case you haven't, it includes this sample code:

import numpy
import matplotlib.pyplot as plt

def CatmullRomSpline(P0, P1, P2, P3, nPoints=100):
  """
  P0, P1, P2, and P3 should be (x,y) point pairs that define the
  Catmull-Rom spline.
  nPoints is the number of points to include in this curve segment.
  """
  # Convert the points to numpy so that we can do array multiplication
  P0, P1, P2, P3 = map(numpy.array, [P0, P1, P2, P3])

  # Calculate t0 to t4
  alpha = 0.5
  def tj(ti, Pi, Pj):
    xi, yi = Pi
    xj, yj = Pj
    return ( ( (xj-xi)**2 + (yj-yi)**2 )**0.5 )**alpha + ti

  t0 = 0
  t1 = tj(t0, P0, P1)
  t2 = tj(t1, P1, P2)
  t3 = tj(t2, P2, P3)

  # Only calculate points between P1 and P2
  t = numpy.linspace(t1,t2,nPoints)

  # Reshape so that we can multiply by the points P0 to P3
  # and get a point for each value of t.
  t = t.reshape(len(t),1)

  A1 = (t1-t)/(t1-t0)*P0 + (t-t0)/(t1-t0)*P1
  A2 = (t2-t)/(t2-t1)*P1 + (t-t1)/(t2-t1)*P2
  A3 = (t3-t)/(t3-t2)*P2 + (t-t2)/(t3-t2)*P3

  B1 = (t2-t)/(t2-t0)*A1 + (t-t0)/(t2-t0)*A2
  B2 = (t3-t)/(t3-t1)*A2 + (t-t1)/(t3-t1)*A3

  C  = (t2-t)/(t2-t1)*B1 + (t-t1)/(t2-t1)*B2
  return C

def CatmullRomChain(P):
  """
  Calculate Catmull Rom for a chain of points and return the combined curve.
  """
  sz = len(P)

  # The curve C will contain an array of (x,y) points.
  C = []
  for i in range(sz-3):
    c = CatmullRomSpline(P[i], P[i+1], P[i+2], P[i+3])
    C.extend(c)

  return C

which nicely computes the interpolation for n >= 4 points like so:

points = [[0,1.5],[2,2],[3,1],[4,0.5],[5,1],[6,2],[7,3]]
c = CatmullRomChain(points)
px, py = zip(*points)
x, y = zip(*c)

plt.plot(x, y)
plt.plot(px, py, 'or')

resulting in this matplotlib image:

enter image description here

Update:

Alternatively, there is a scipy.interpolate function for BarycentricInterpolator that appears to do what you're looking for. It is rather straightforward to use and works for cases in which you have only 3 data points.

from scipy.interpolate import BarycentricInterpolator

# create some data points
points1 = [[0, 2], [1, 4], [2, -2], [3, 6], [4, 2]]
points2 = [[1, 1], [2, 5], [3, -1]]

# put data into x, y tuples
x1, y1 =zip(*points1)
x2, y2 = zip(*points2)

# create the interpolator
bci1 = BarycentricInterpolator(x1, y1)
bci2 = BarycentricInterpolator(x2, y2)

# define dense x-axis for interpolating over
x1_new = np.linspace(min(x1), max(x1), 1000)
x2_new = np.linspace(min(x2), max(x2), 1000)

# plot it all
plt.plot(x1, y1, 'o')
plt.plot(x2, y2, 'o')
plt.plot(x1_new, bci1(x1_new))
plt.plot(x2_new, bci2(x2_new))
plt.xlim(-1, 5)

enter image description here

Update 2

Another option within scipy is akima interpolation via Akima1DInterpolator. It is as easy to implement as Barycentric, but has the advantage that it avoids large oscillations at the edge of a data set. Here's a few test cases that exhibit all the criteria you've asked for so far.

from scipy.interpolate import Akima1DInterpolator

x1, y1 = np.arange(13), np.random.randint(-10, 10, 13)
x2, y2 = [0,2,3,6,12], [100,50,30,18,14]
x3, y3 = [4, 6, 8], [60, 80, 40]

akima1 = Akima1DInterpolator(x1, y1)
akima2 = Akima1DInterpolator(x2, y2)
akima3 = Akima1DInterpolator(x3, y3)

x1_new = np.linspace(min(x1), max(x1), 1000)
x2_new = np.linspace(min(x2), max(x2), 1000)
x3_new = np.linspace(min(x3), max(x3), 1000)

plt.plot(x1, y1, 'bo')
plt.plot(x2, y2, 'go')
plt.plot(x3, y3, 'ro')
plt.plot(x1_new, akima1(x1_new), 'b', label='random points')
plt.plot(x2_new, akima2(x2_new), 'g', label='exponential')
plt.plot(x3_new, akima3(x3_new), 'r', label='3 points')
plt.xlim(-1, 15)
plt.ylim(-10, 110)
plt.legend(loc='best')

enter image description here

lanery
  • 5,222
  • 3
  • 29
  • 43
  • Thanks for the prompt response and example. This does not answer my original question though: 1) The spline does not go through the first and last points 2) It does not work for n=3 . Also, does this mean that there is no way of emulating Excel's spline function using the interpolate library? I'd be grateful for any further suggestions. – Dave May 14 '16 at 10:06
  • Referring back to the wikipedia page, Catmull-Rom "is a type of interpolating spline **defined by four control points, `p0`, `p1`, `p2`, and `p3`, with the curve drawn only drawn from `p1` to `p2`**," so the spline will never go through the first and last points. My guess is that Excel makes some sort of reasonable extrapolation when it does the spline fit such that the splines typically look like they fit well. – lanery May 14 '16 at 10:28
  • And to answer your second question, I couldn't find anything in `scipy.interpolate` that looks like its using Catmull-Rom. That said, there's another spline method, [`BarycentricInterpolator`](http://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.BarycentricInterpolator.html#scipy.interpolate.BarycentricInterpolator) that looks like it may do what you want. Give me a minute to update my answer. – lanery May 14 '16 at 10:44
  • See also: http://stackoverflow.com/questions/1251438/catmull-rom-splines-in-python – lanery May 14 '16 at 11:05
  • Many thanks - had seen http://stackoverflow.com/questions/1251438 but couldn't work out how to apply it to get the behaviour I wanted. I wasn't aware of BarycentricInterpolator - this is much easier to use and does exactly what I needed. 1+ – Dave May 16 '16 at 05:49
  • Subsiduary related question: Barycentric gives nice result for n=3, but w data points for an exponential curve eg x = [0,2,3,6,12], y = [100,50,30,18,14], Barycentric gives ugly results. In contrast, Catmull-Rom gives a good fit for exponentials, but doesn't make the spline pass through the first and last points. Excel produces nice exponential fit splines from just 3 points, so I guess there must be more to the algorithm than just Catmull-Rom splines... Is there any way to tweak the above examples, or another different approach, which will allow me to emulate Excel's behaviour to do this? – Dave May 16 '16 at 10:29
  • @Dave, 3rd times the charm! While I think the akima interpolator suits your needs, I couldn't find any evidence to support that it is what Excel uses, so I hope that's not a problem. My guess is that Excel is using Catmull-Rom plus some sophisticated extrapolation methods. I've also seen Bezier curves mentioned. Either way, I hope my new update solves your problem. – lanery May 17 '16 at 03:03
  • 1
    BarycentricInterpolator fits a *single (high-degree)* polynomial to each coordinate of the points, using Lagrange interpolation. This is known to be ill-conditioned for many (say more than 8) points. In contrast, other schemes use splines which are *piecewise low-degree* polynomials. I certainly would not recommend BarycentricInterpolator for the general case. – Hugues Mar 27 '21 at 18:35
0

@Lanery: Re: Update 2: The best just got better!

Had to redefine the lists x2,y2,x3,y3 as numpy arrays to get your example to work on my system (Spyder / Python 2.7) :

x2 = np.array([0,2,3,6,12])
y2 = np.array([100,50,30,18,14])
x3 = np.array([4, 6, 8])
y3 = np.array([60, 80, 40])

But now works like a dream! Many thanks for your expertise and clear explanations.

DaveW
  • 185
  • 1
  • 8