ParametricSpline_example2

import numpy as np
import scipy.interpolate as si
import matplotlib.pyplot as plt
from string import split, join

def readXY(filepath="./XYdata.txt"):

 data=open(filepath, "r")
 X=[]
 Y=[]
 for line in data:
  splitline=split(line)
  try:
   X.append(float(splitline[0]))
   Y.append(float(splitline[1]))
  except:
   pass
 data.close()
 return X, Y

def getControlPoints(knots, k):

 n = len(knots)-1-k
 cx = np.zeros(n)
 for i in range(n):
  tsum = 0
  for j in range(1,k+1):
   tsum+=knots[i+j]
  cx[i] = float(tsum)/k
 return cx

def plotSplineFunction(title = "Spline Function", offset_x = 0.5, offset_y = 0.05):

 plt.plot(coeffs_x, coeffs_y, 'go-')
 plt.plot(x, y, 'ro')
 plt.plot(xP, yP, lw=2)
 plt.plot(knots, knotsy, 'gs') # ,markersize=10
 plt.axis([xmin-offset_x, xmax+offset_x, ymin-offset_y, ymax+offset_y])
 plt.grid(True)
 plt.title(title)
 plt.show()

def findLSQSplineRep():
 global knots, knots_full, coeffs_y, coeffs_x, yP, knotsy

 lsqspline = si.LSQUnivariateSpline( x, y, knots, k=k, w = w)
 knots = lsqspline.get_knots()
 coeffs_y = lsqspline.get_coeffs()
 yP = lsqspline(xP)
 knotsy = lsqspline(knots)
 knots_full = np.concatenate(([knots[0]]*k, knots, [knots[-1]]*k ))
 coeffs_x = getControlPoints(knots_full, k)

def findSmoothedSplineRep(s=0.005):
 global knots, knots_full, coeffs_y, coeffs_x, yP, knotsy

 spline = si.UnivariateSpline( x, y, k=k , s=s, w=w)
 knots = spline.get_knots()
 coeffs_y = spline.get_coeffs()
 yP = spline(xP)
 knotsy = spline(knots)
 knots_full = np.concatenate(([knots[0]]*k, knots, [knots[-1]]*k ))
 coeffs_x = getControlPoints(knots_full, k)


##################################### LOAD DATA POINTS #######################################

x, y = readXY("./data_charts/cy_f004.txt")

# print "x:", x
# print "y:", y

### set variables
nsample = 100
nknot = 3
k=3
num_points=len(x)
xmax=max(x)
xmin=min(x)
ymax=max(y)
ymin=min(y)
xP = np.linspace( x[0], x[-1], nsample )
# xP = np.linspace( x[0]-2, x[-1]+2, nsample )

plt.plot(x, y, 'ro')
plt.axis([xmin-0.5, xmax+0.5, ymin-0.05, ymax+0.05])
plt.grid(True)
plt.title("Data points")
plt.show()

### define weight vector to push further approximations stick to end points
wend = 3
w = [wend] + [1]*(num_points-2) + [wend]


################################ FIND LSQ SPLINE REPRESENTATION  ##############################

### try with uniform knot vector

knot_offset=(xmax-xmin)/(nknot+1)
knots = np.linspace( knot_offset, xmax-knot_offset, nknot )

findLSQSplineRep()
plotSplineFunction(title = "LSQ spline function with uniform knot vector", offset_x = 1.)

### manually select best fitting non-uniform knot vector

knots = [1.2, 1.85, 2.1, 2.4, 2.6] # !!
knots = [1.2, 1.85, 2.3, 2.6, 3.2] # !
knots = [1.2, 1.85, 2.0, 2.7, 3.2] # !!!
knots = [1.2, 1.85, 2.0, 2.8, 3.5] # !!!!

findLSQSplineRep()
plotSplineFunction(title = "LSQ spline function with manual selected non-uniform knot vector", offset_x = 1.)

############################        FIND SMOOTHED SPLINE REPRESENTATION   ###########################

findSmoothedSplineRep(s=0.1)
plotSplineFunction(title = "Smoothed spline function s=0.1")
# print "s=0.1, knots:", knots
# findSmoothedSplineRep(s=0)
# plotSplineFunction(title = "Polynomial spline interpolation")
# print "s=0, knots:", knots
findSmoothedSplineRep()
plotSplineFunction(title = "Smoothed spline function s=0.005")
# print "s=0.05, knots:", knots

############################    MANUALLY SELECT BEST_FITTING CONTROL POINTS   #############################

### try to comment next two rows to look at the original smoothed spline curve after parametrization
coeffs_x = [ x[0], 1.1,   1.7,  2.16,  3.1,   3.0,   x[-1]]
coeffs_y = [ y[0], 0.11,  0.17,  0.32,  0.27,  0.19,  y[-1]]

############################        PARAMETRIC SPLINE REPRESENTATION       #############################

### norm knot vector to number of points

num_knots = len(knots_full)
ka = (knots_full[-1] - knots_full[0])/(num_points)
knotsp = np.zeros(num_knots)
for i in range(num_knots):
 knotsp[i] = num_points-((knots_full[-1] - knots_full[i]))/ka
# print "knotsp:", knotsp

### find parametric spline representation based on control points coordinates (coeffs_x, coeffs_y)

tckX = knotsp , coeffs_x, k
tckY = knotsp , coeffs_y, k

splineX = si.UnivariateSpline._from_tck(tckX)
splineY = si.UnivariateSpline._from_tck(tckY)

coeffs_p = getControlPoints(knotsp, k)
tP = np.linspace(0,  num_points, 100)
# tP = np.linspace(-20,  num_points+20, 100)
xP = splineX(tP)
yP = splineY(tP)

################################## PLOT PARAMETRIC SPLINE CURVE #####################################

offset_x = (xmax-xmin)*0.05
offset_y = (ymax-ymin)*0.1
ct = 0.5
tadd = 0.5

# offset_x = (xmax-xmin)*0.2
# offset_y = (ymax-ymin)*0.2
# ct = 0.7
# tadd = 1.0

knotpoints_y = [ymin-offset_y*ct]*len(knotsp)
knotpoints_x = [xmin-offset_x*ct]*len(knotsp)

fig = plt.figure()

ax = fig.add_subplot(224)
ax.grid(True)
plt.plot(coeffs_x, coeffs_p, '-og')
plt.plot(xP, tP, 'r', lw=2)
plt.plot(knotpoints_x, knotsp, '>', ms=6, color='black')
plt.ylim([-tadd, num_points+tadd])
plt.xlim([xmin - offset_x, xmax + offset_x])
plt.ylabel('t', rotation=0,  labelpad=20, fontweight='bold', fontsize=14)
plt.xlabel('x', labelpad=5, fontweight='bold', fontsize=14)
ax.invert_yaxis()
plt.title('Spline function x(t)')

ax = fig.add_subplot(221)
ax.grid(True)
plt.plot(coeffs_p, coeffs_y, '-og')
plt.plot(tP, yP, 'r', lw=2)
plt.plot(knotsp,knotpoints_y, '^', ms=6, color='black')
plt.xlim([-tadd, num_points+tadd])
plt.ylim([ymin - offset_y, ymax + offset_y])
plt.ylabel('y', labelpad=20, rotation=0, fontweight='bold', fontsize=14)
plt.xlabel('t', labelpad=10, fontweight='bold', fontsize=14)
ax.invert_xaxis()
plt.title('Spline function y(t)')

ax = fig.add_subplot(222)
ax.grid(True)
plt.plot(x, y, 'ro')
plt.plot(coeffs_x, coeffs_y, '-og')
plt.plot(xP, yP, 'b', lw=2.5)
plt.xlim([xmin - offset_x, xmax + offset_x])
plt.ylim([ymin - offset_y, ymax + offset_y])
plt.title('Spline curve f(x(t), y(t))')

plt.show()