From 605bfb16a23716d80c42c8ab8d9abb5bed6ae45c Mon Sep 17 00:00:00 2001 From: wirawan Date: Fri, 22 Jan 2010 18:48:19 +0000 Subject: [PATCH] * Simplified spline interpolation tools for 2-D curves (single curve or piecewise curve). --- math/spline_2d.py | 91 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 91 insertions(+) create mode 100644 math/spline_2d.py diff --git a/math/spline_2d.py b/math/spline_2d.py new file mode 100644 index 0000000..431ed66 --- /dev/null +++ b/math/spline_2d.py @@ -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)