On 9 January 2013 16:02, <raph...@mameghani.de> wrote:
> Hi,
>
> I want to interpolate (with quadratic splines) a stack of 2D-arrays/matrices
> y1, y2, y3, ... in a third dimension (which I call x) e.g. for crossfading
> images. I already have a working code which unfortunately still contains two
> explicit loops over the rows and colums of the matrices. Inside these loops I
> simply use 'interp1d' from scipy suitable for 1D-interpolations. Is anybody
> here aware of a better, more efficient solution of my problem? Maybe
> somewhere out there a compiled routine for my problem already exists in a
> python library... :-)
It's possible. I wouldn't be surprised if there wasn't any existing
code ready for you to use.
>
> My code:
>
> -----============================================-----
> from scipy.interpolate import interp1d
> from numpy import array, empty_like, dstack
>
> x = [0.0, 0.25, 0.5, 0.75, 1.0]
>
> y1 = array([[1, 10, 100, 1000], [1, 10, 100, 1000]], float)
> y2 = array([[2, 20, 200, 2000], [2, 20, 200, 2000]], float)
> y3 = array([[3, 30, 300, 3000], [4, 40, 400, 4000]], float)
> y4 = array([[4, 40, 400, 4000], [8, 80, 800, 8000]], float)
> y5 = array([[5, 50, 500, 5000], [16, 160, 1600, 16000]], float)
>
> y = dstack((y1, y2, y3, y4, y5))
>
> y_interpol = empty_like(y[:, :, 0])
> i_range, j_range = y.shape[:2]
>
> for i in xrange(i_range):
> for j in xrange(j_range):
> # interpolated value for x = 0.2
> y_interpol[i,j] = interp1d(x, y[i, j,:], kind='quadratic')(0.2)
>
> print y_interpol
> -----============================================-----
Since numpy arrays make it so easy to form linear combinations of
arrays without loops I would probably eliminate the loops and just
form the appropriate combinations of the image arrays. For example, to
use linear interpolation you could do:
def interp_frames_linear(times, frames, t):
'''times is a vector of floats
frames is a 3D array whose nth page is the image for time t[n]
t is the time to interpolate for
'''
# Find the two frames to interpolate between
# Probably a better way of doing this
for n in range(len(t)-1):
if times[n] <= t < times[n+1]:
break
else:
raise OutOfBoundsError
# Interpolate between the two images
alpha = (t - times[n]) / (times[n+1] - times[n])
return (1 - alpha) * frames[:, :, n] + alpha * frames[:, :, n+1]
I'm not really sure how quadratic interpolation is supposed to work
(I've only ever used linear and cubic) but you should be able to do
the same sort of thing.
Oscar
--
http://mail.python.org/mailman/listinfo/python-list
