# Cubic spline interpolation results

I have a set of data points on which i am trying to do cubic spline interpolation. Below is the snapshot of the curve with the input data points marked in green color. And the red color marked point is the test point.

The steps i am following to interpolate is

• Form piecewise spline equations between points . cubic equation : ax^3 + bx^2 + cx + d = P(x)
• Form first derivative equations at interior points
• Form Second derivative equations at interior points
• Set Second derivative equations at exterior points to zero

Get the equations created above in matrix form and solve for the coefficients(a,b,c,d) of each spline.

But for any given test point (X) other than input points.... the spline equation does not yield proper (Y) .

As a test case in the above curve.....if i use X = 152.73 (red colour marked point on the curve) in the first spline equation ... then,

Y = (966.3375) * (152.73^3) + (-443577.902) * (152.73^2) + (67872238.9) * (152.73) + (-3461748870) = 613.07 which is wrong . Y should be returned as 699.92

can anyone explain what is going wrong in this? I tried Quadratic spline interpolation as well, it also returns same kind of results. I hope i am following the correct steps .If not, Please direct me in doing so. I am very new to this stuff. Please help.

Sorry for my bad english and math terminologies

• You give the form of the "cubic" as ax^3 + bx^2 + c^x + d = P when it should be ax^3 + bx^2 + cx + d = P. Was that just a typo here or did you also make that mistake in your calculations. – user247327 May 17 '16 at 15:41
• @user247327 Typo here – 911Rapid May 17 '16 at 15:42
• The lay of the data dictates that you should let $x=ay^3+by^2+cy+d$ or else use a different family of splines like Bezier. – user5713492 May 17 '16 at 17:49
• @user5713492 I dont understand how it matters. Changing X to Y and vice-versa. But anways, i will try that. – 911Rapid May 17 '16 at 18:15
• @user5713492 Btw, I am not totally sure....But aren't Beizer curves are for approximation as opposed to interpolation ?? And , i think they can't be used for very precise calculations and stuff in like engineering drawings. Correct me, if i am wrong – 911Rapid May 17 '16 at 18:18

The problem was that the data looked kind of like this:

For this data, when an attempt was made to fit $y=a_ix^3+b_ix^2+c_ix+d_i$ between the spline points and it didn't work well because the data can't be expressed as a function of $x$; instead there can be multiple values of $y$ for each value of $x$. There are spline algorithms that can handle this, but if it is desired to use an equation as noted above, in this case we can instead fit $x=a_iy^3+b_iy^2+c_iy+d_i$ in each interval. For these data we anticipate good results just by looking at the data in that the graph seems to depict $x$ as a smooth function of $y$.

If you are set on keeping the data as-is and not flip x/y as others have pointed out, I might suggest using Catmull-Rom, maybe with a chordial (alpha=1) parameterization if you like keeping the "humps" large as in a cubic spline (though I prefer centripedal, though uniform is simplest to work with).

Here's an example from python i made using those points: (http://imgur.com/CwMMou4) (not enough rep to post images directly :/)

Here is another example with some more extreme points to show how the 3 parameterizations behave differently: View post on imgur.com

Notice 3 lines for different alpha values: Green line = uniform (a=0) Blue = Centripedal (a=0.5, no cusps or loops) Red = Chordial (a=1)

Notice though that I had to arbitrarily create endpoints (I added custom endpoints that maintained slope between first and last 2 points), otherwise catmull rom won't be able to interpolate first and last two points, without knowing which tangent you want to use at the start/end.

My Python Code if you'd like:

# coding: utf-8

# In[1]:

import numpy
import matplotlib.pyplot as plt
get_ipython().magic(u'pylab inline')

# In[2]:

def CatmullRomSpline(P0, P1, P2, P3, a, 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 = a
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,alpha):
"""
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],alpha)
C.extend(c)

return C

# In[15]:

# Define a set of points for curve to go through
#Points = numpy.random.rand(12,2)
Points=array([array([153.01,722.67]),array([152.73,699.92]),array([152.91,683.04]),array([154.6,643.45]),
array([158.07,603.97])])
x1=Points[0][0]
x2=Points[1][0]
y1=Points[0][1]
y2=Points[1][1]
x3=Points[-2][0]
x4=Points[-1][0]
y3=Points[-2][1]
y4=Points[-1][1]
dom=max(Points[:,0])-min(Points[:,0])
rng=max(Points[:,1])-min(Points[:,1])
prex=x1+sign(x1-x2)*dom*0.01
prey=(y1-y2)/(x1-x2)*(prex-x1)+y1
endx=x4+sign(x4-x3)*dom*0.01
endy=(y4-y3)/(x4-x3)*(endx-x4)+y4
print len(Points)
Points=list(Points)
Points.insert(0,array([prex,prey]))
Points.append(array([endx,endy]))
print len(Points)

# In[16]:

#Define alpha
a=0.

# Calculate the Catmull-Rom splines through the points
c = CatmullRomChain(Points,a)

# Convert the Catmull-Rom curve points into x and y arrays and plot
x,y = zip(*c)
plt.plot(x,y,c='green',zorder=10)

# Plot the control points
px, py = zip(*Points)
plt.plot(px,py,'or')

a=0.5
c = CatmullRomChain(Points,a)
x,y = zip(*c)
plt.plot(x,y,c='blue')

a=1.
c = CatmullRomChain(Points,a)
x,y = zip(*c)
plt.plot(x,y,c='red')

plt.grid(b=True)
plt.show()

# In[14]:

Points