parent
5af79a8000
commit
605bfb16a2
1 changed files with 91 additions and 0 deletions
@ -0,0 +1,91 @@ |
||||
# $Id: spline_2d.py,v 1.1 2010-01-22 18:48:19 wirawan Exp $ |
||||
# |
||||
# wpylib.math.spline_2d module |
||||
# Created: 20091204 |
||||
# Wirawan Purwanto |
||||
# |
||||
|
||||
import numpy |
||||
import scipy.interpolate |
||||
|
||||
class spline_2d: |
||||
"""Simple interpolation for two-dimensional curve. |
||||
Created to handle the quirks of current SciPy interpolation |
||||
routines.""" |
||||
# Important notes on spline CAVEATS: |
||||
# - for some reason we HAVE to make a copy of the 'x' array |
||||
# (to make it contiguous, probably?) |
||||
# - also, the x values better be sorted in ascending order, or else |
||||
# the routine would return nonsense (i.e. NaN's). |
||||
def __init__(self, x, y): |
||||
self.init(x,y) |
||||
|
||||
def init(self, x, y): |
||||
# First, the x must be sorted, so we make a private copy of |
||||
# the data: |
||||
self.data = numpy.array(zip(x, y), dtype=[('x', float), ('y', float)]) |
||||
# Quirk 1: The x axis data must be sorted ascending |
||||
self.data.sort(order=['x']) |
||||
self.x = self.data['x'] |
||||
self.y = self.data['y'] |
||||
# Quirk 2: the x data for spline function must be contiguous |
||||
# See below in init_spline_params() |
||||
#self.x_copy = self.x.copy() |
||||
|
||||
try: |
||||
del self.spline_params |
||||
except: |
||||
pass |
||||
|
||||
def init_spline_params(self): |
||||
"""Initialize spline params with default params. |
||||
You can call something to initialize the spline params before |
||||
calling the first spline function if you want different, non-default |
||||
parameters.""" |
||||
self.spline_params \ |
||||
= scipy.interpolate.splmake(self.x.copy(), self.y) |
||||
|
||||
def spline(self, xnew): |
||||
try: |
||||
params = self.spline_params |
||||
except: |
||||
self.init_spline_params() |
||||
return scipy.interpolate.spleval(self.spline_params, xnew) |
||||
|
||||
|
||||
|
||||
class spline_2d_piecewise: |
||||
"""Simple spline_2d interpolator with piecewise datasets. |
||||
Interpolation is possible only in the ranges defined by the piecewise |
||||
datasets. |
||||
No checking is done whether the pieces are overlapping, discontinuous, etc. |
||||
The first piece found enclosing the coordinate will be taken for |
||||
interpolation.""" |
||||
def __init__(self, *datasets): |
||||
self.init(*datasets) |
||||
|
||||
def init(self, *datasets): |
||||
#if len(dsets) % 2: |
||||
# raise ValueError, "The input datasets must be given in x, y pairs |
||||
self.pieces = [] |
||||
for dset in datasets: |
||||
x = dset[0] |
||||
y = dset[1] |
||||
xmin = numpy.min(x) |
||||
xmax = numpy.max(x) |
||||
piece = spline_2d(x, y) |
||||
piece.xmin = xmin |
||||
piece.xmax = xmax |
||||
self.pieces.append(piece) |
||||
|
||||
def in_range(self, piece, x): |
||||
return piece.xmin <= x and x <= piece.xmax |
||||
|
||||
def get_piece(self, x): |
||||
for p in self.pieces: |
||||
if self.in_range(p, x): |
||||
return p |
||||
raise ValueError, "Out-of-range x value = %g" % x |
||||
|
||||
def spline(self, x): |
||||
return self.get_piece(x).spline(x) |
Loading…
Reference in new issue