[SciPy-User] [SciPy-user] 3D interpolation

denis denis-bz-gg@t-online...
Mon Dec 28 12:26:31 CST 2009


Vincent, you're right, splines can overshoot, though not by a lot:
0 1 1 0 overshoots by 20 % in the code snippet below.
Burak, does order=1, just averaging 8 neighbors, work ?

There are many kinds of cubic splines, including
-- global cubic, which map_coordinates seems to use
-- local B-spline: does not interpolate, go through the input points,
    but can't overshoot, stays in the convex hull of the input points
-- local Catmull-Rom: interpolates, but can overshoot (over shoots and
leaves :)
    ("Local" means that e.g. in 1d, out[ 3 to 4 ] is calculated from
the 4 nearest input points at 2 3 4 5,
    or equivalently that in[3] affects only out[ 1 to 5 ]: a desirable
property.)

Different corners of Scipy have different spline routines, some with
doc.
To see exactly what they do, you can only plot their impulse response
like so.

cheers
  -- denis

""" plot map_coordinates( 1d spike )
"""

from __future__ import division
import sys
import numpy as np
from scipy.ndimage import map_coordinates
import pylab as pl

N = 10
exec( "\n".join( sys.argv[1:] ))  # N= ...
np.set_printoptions( 2, threshold=100, suppress=True )  # .2f

a = np.r_[ 5*[0], 1., 5*[0], 1,1, 5*[0], 5*[1] ]  # 0 1 1 0  -> 1.2
na = len(a) - 1
x = np.linspace( 0, na, na*N + 1 )
print "a:", a.astype(int)

for order in (3,):
    z = map_coordinates( a, [x], order=order )
    print "z[%d] halfint:" % order, z[N//2::N]
    pl.plot( x, z )

pl.show()


More information about the SciPy-User mailing list