Path: csiph.com!v102.xanadu-bbs.net!xanadu-bbs.net!goblin2!goblin.stu.neva.ru!newsfeed.xs4all.nl!newsfeed3.news.xs4all.nl!xs4all!post.news.xs4all.nl!not-for-mail
MIME-Version: 1.0
In-Reply-To: <5b5b7e09-7482-4609-b774-d143de366615@googlegroups.com>
References: <5b5b7e09-7482-4609-b774-d143de366615@googlegroups.com>
Date: Wed, 9 Jan 2013 22:59:56 +0000
Subject: Re: Interpolating/crossfading a stack of matrices
From: Oscar Benjamin <oscar.j.benjamin@gmail.com>
To: raphael@mameghani.de
Content-Type: text/plain; charset=ISO-8859-1
Content-Transfer-Encoding: quoted-printable
Cc: Python List <python-list@python.org>
Precedence: list
Newsgroups: comp.lang.python
Message-ID: <mailman.343.1357772410.2939.python-list@python.org>
Lines: 74
NNTP-Posting-Host: 2001:888:2000:d::a6
Xref: csiph.com comp.lang.python:36524

On 9 January 2013 16:02,  <raphael@mameghani.de> wrote:
> Hi,
>
> I want to interpolate (with quadratic splines) a stack of 2D-arrays/matri=
ces y1, y2, y3, ... in a third dimension (which I call x) e.g. for crossfad=
ing images. I already have a working code which unfortunately still contain=
s two explicit loops over the rows and colums of the matrices. Inside these=
 loops I simply use 'interp1d' from scipy suitable for 1D-interpolations. I=
s anybody here aware of a better, more efficient solution of my problem? Ma=
ybe 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:
>
> -----=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D-----
> from scipy.interpolate import interp1d
> from numpy import array, empty_like, dstack
>
> x =3D [0.0, 0.25, 0.5, 0.75, 1.0]
>
> y1 =3D array([[1, 10, 100, 1000], [1,  10,   100,  1000]], float)
> y2 =3D array([[2, 20, 200, 2000], [2,  20,   200,  2000]], float)
> y3 =3D array([[3, 30, 300, 3000], [4,  40,   400,  4000]], float)
> y4 =3D array([[4, 40, 400, 4000], [8,  80,   800,  8000]], float)
> y5 =3D array([[5, 50, 500, 5000], [16, 160, 1600, 16000]], float)
>
> y =3D dstack((y1, y2, y3, y4, y5))
>
> y_interpol =3D empty_like(y[:, :, 0])
> i_range, j_range =3D y.shape[:2]
>
> for i in xrange(i_range):
>     for j in xrange(j_range):
>         # interpolated value for x =3D 0.2
>         y_interpol[i,j] =3D interp1d(x, y[i, j,:], kind=3D'quadratic')(0.=
2)
>
> print y_interpol
> -----=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D-----

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] <=3D t < times[n+1]:
            break
    else:
        raise OutOfBoundsError

    # Interpolate between the two images
    alpha =3D (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