#
# pychip.py
# Michalski
# 20090818
# Piecewise =
cubic Hermite =
interpolation (monotonic...) in Python
#
# References:
#
# =
Wikipedia: Monotone cubic =
interpolation
# =
Cubic Hermite spline
#
# A cubic =
Hermte spline is a =
third degree spline with each polynomial of the spline
# in Hermite form. The Hermite form consists of two control points and two =
control
# tangents =
for each polynomial. Each interpolation is performed on one =
sub-interval
# at a time =
(piece-wise). A monotone cubic interpolation is a variant of =
cubic
# =
interpolation that preserves monotonicity of the data to be interpolated =
(in other
# words, it =
controls overshoot). Monotonicity is preserved by linear =
interpolation
# but not by =
cubic interpolation.
#
# Use:
# There are =
two separate calls, the first call, pchip_init(), computes the slopes =
that
# the interpolator needs. If =
there are a large number of points to compute,
# it is more efficient to compute the =
slopes once, rather than for every point
# =
pchip_init() along with X, Y, and a vector of desired "xnew"s and computes a =
vector
# of "ynew"s. If only a =
handful of points is needed, pchip() is a third function
# which combines a call to pchip_init() =
followed by pchip_eval().
#
import pylab as =
P
def pchip(x, y, xnew):
# Compute the slopes used by the piecewise =
cubic Hermite interpolator
# Use these slopes (along =
with the Hermite basis =
function) to interpolate
ynew =3D pchip_eval(x, =
y, xnew)
return =
ynew
def x_is_okay(x,xvec):
# Make sure "x" and "xvec" satisfy the conditions =
for
# running the pchip interpolator
n =3D len(x)
=
# Make sure "x" is in sorted order (brute force, but =
works...)
xx =3D x.copy()
total_matches =3D (xx =
=3D=3D x).sum()
if total_matches !=3D n:
=
print "*" * 50
=
print "x_is_okay()"
=
print "x values weren't in sorted order --- =
aborting"
return False
=
# Make sure 'x' doesn't have any repeated values
1:] - =
x[:-1]
if (delta =3D=3D =
0.0).any():
print "*" * 50
print =
"x_is_okay()"
print "x values =
weren't monotonic--- aborting"
return False
# Check for in-range xvec values (beyond upper =
edge)
check =3D xvec > x[-1]
if check.any():
print "*" * 50
=
print "x_is_okay()"
=
print "Certain 'xvec' values are beyond the upper end of =
'x'"
print "x_max =3D =
", x[-1]
print "out-of-range =
xvec's =3D ", xvec[indices]
print =
"out-of-range xvec indices =3D ", =
indices
return False
=
# Second - check for in-range xvec values (beyond lower edge)
0]
if check.any():
print "*" * 50
=
print "x_is_okay()"
=
print "Certain 'xvec' values are beyond the lower end of =
'x'"
print "x_min =3D =
", x[0]
print "out-of-range =
xvec's =3D ", xvec[indices]
print =
"out-of-range xvec indices =3D ", =
indices
return False
=
return True
def pchip_eval(x, y, m, xvec):
# Evaluate the piecewise =
cubic Hermite interpolant with monoticity =
preserved
#
# x =3D array containing the =
x-data
# y =3D array containing the =
y-data
# m =3D slopes at each (x,y) =
point [computed to preserve monotonicity]
# =
xnew =3D new "x" value =
where the interpolation is desired
#
# x must be sorted low to =
high... (no repeats)
# y can have repeated =
values
#
# This works with either a scalar or =
vector of "xvec"
mm =3D =
len(xvec)
=
############################
# Make sure there aren't =
problems with the input data
=
############################
if not =
x_is_okay(x, xvec):
print =
"pchip_eval2() - ill formed 'x' vector!!!!!!!!!!!!!"
# Cause a =
hard crash...
=
STOP_pchip_eval2
# =
Find the indices "k" such that x[k] < xvec < x[k+1]
=
# Create "copies" of "x" as rows in a mxn 2-dimensional =
vector
xx =3D =
P.resize(x,(mm,n)).transpose()
xxx =3D xx > =
xvec
=
# Compute column by column differences
1,:] - =
xxx[1:,:]
# Collapse over =
rows...
k =3D z.argmax(axis=3D0)
=
# Create the Hermite coefficients
h =3D x[k+1] - x[k]
t =3D (xvec - =
x[k]) / h[k]
=
# Hermite basis =
functions
h00 =3D (2 * t**3) - (3 * t**2) + 1
3 - (2 * =
t**2) + t
2* =
t**3) + (3 * t**2)
h11 =3D t**3 - t**2
=
# Compute the interpolated value of "y"
1] + h11*h*m[k+1]
return =
ynew
def pchip_init(x,y):
=
# Evaluate the piecewise cubic Hermite interpolant with monoticity preserved
#
# x =3D array =
containing the x-data
# y =3D array containing the =
y-data
#
# x must be sorted low to =
high... (no repeats)
# y can have repeated =
values
#
# x input conditioning is =
assumed but not checked
=
# Compute the slopes of the secant lines between successive points
1:] - =
y[:-1]) / (x[1:] - x[:-1])
=
# Initialize the tangents at every points as the average of the =
secants
'd')
=
# At the endpoints - use one-sided differences
0] =3D =
delta[0]
1] =3D =
delta[-1]
# In the middle - use the =
average of the secants
m[1:-1] =3D =
(delta[:-1] + delta[1:]) / 2.0
=
# Special case: intervals where y[k] =3D=3D y[k+1]
# Setting these slopes to =
zero guarantees the spline connecting
# these points =
will be flat which preserves monotonicity
indices_to_fix =
=3D P.compress((delta =3D=3D 0.0), =
range(n))
# print "zero slope indices to fix =3D =
", indices_to_fix
for ii in indices_to_fix:
0.0
m[ii+1] =3D 0.0
1]/delta
beta =3D m[1:]/delta
dist =3D =
alpha**2 + beta**2
tau =3D 3.0 / P.sqrt(dist)
# To prevent overshoot or =
undershoot, restrict =
the position vector
# (alpha, beta) to a circle of radius =
3. If (alpha**2 + beta**2)>9,
# then set m[k] =3D =
tau[k]alpha[k]delta[k] and m[k+1] =3D tau[k]beta[b]delta[k]
=
# where tau =3D 3/sqrt(alpha**2 + beta**2).
# Find the indices that need =
adjustment
over =3D (dist > 9.0)
indices_to_fix =3D =
P.compress(over, range(n))
# print "overshoot indices to fix... =3D =
", indices_to_fix
for ii in indices_to_fix:
m[ii+1] =3D tau[ii] * beta[ii] * =
delta[ii]
return m
def CubicHermiteSpline(x, y, =
x_new):
=
# Piecewise Cubic Hermite Interpolation using Catmull-Rom
=
# method for computing the slopes.
#
# Note - this only works =
if delta-x is uniform?
=
# Find the two points which "bracket" "x_new"
False
for ii in range(len(x)-1):
if (x[ii] <=3D x_new) and (x[ii+1] > =
x_new):
found_it =
=3D True
break
if not found_it:
print
print "requested =
x=3D<%f> outside X range[%f,%f]" % =
(x_new, x[0], x[-1])
# Starting and ending data =
points
x0 =3D x[ii]
1]
y1 =3D y[ii+1]
=
# Starting and ending tangents (using Catmull-Rom spline method)
# Handle special cases =
(hit one of the endpoints...)
if ii =3D=3D 0:
=
# Hit lower endpoint
m0 =3D =
(y[1] - y[0])
m1 =3D (y[2] - y[0]) / 2.0
elif ii =3D=3D =
(len(x) - 2):
# Hit upper =
endpoints
1] - y[ii-1]) / =
2.0
=
m1 =3D (y[ii+1] - y[ii])
else:
=
# Inside the field...
m0 =3D =
(y[ii+1] - y[ii-1])/ 2.0
m1 =3D (y[ii+2] - y[ii]) / 2.0
=
# Normalize to x_new to [0,1] interval
t =3D (x_new - x0) / =
h
# Compute the four Hermite basis =
functions
h00 =3D ( 2.0 * t**3) - =
(3.0 * t**2) + 1.0
1.0 * =
t**3) - (2.0 * t**2) + =
t
2.0 * =
t**3) + (3.0 * t**2)
h11 =3D ( 1.0 * t**3) - =
(1.0 * t**2)
1
y_new =3D (h00 * y0) + =
(h10 * h * m0) + (h01 * y1) + (h11 * h * m1)
return y_new
def main():
=
############################################################
<=
div style=3D"margin-top: 0px; margin-right: 0px; margin-bottom: 0px; =
margin-left: 0px; font: normal normal normal 11px/normal Monaco; color: =
rgb(192, 35, 30); "> =
# Sine wave test