[SciPy-User] Efficient 3d smoothing/interpolation

[Jose Gomez-Dans]

Hi,

A two part question , I guess! I have a set of images (2d datasets), and I
want to interpolate from one timestep to the next for each pixel when the
data is invalid. Since I can get invalid data on the edges (missing first
few time steps), I want to be able to set these values to the closest (in
time) valid values for those pixels.

Let's say I have a 3D dataset: a 2d (2000x2000) dataset with 180 time
frames. Let's call it A
>>> print A.shape
(180, 2000, 2000)

Some of the values in the dataset will be missing (flagged by a known
value). I would like to interpolate over them (and maybe do some smoothing,
but that's easy with splines) in time. Additionally, since the ends of my
time series for a given "pixel" may be missing, I want to set these
endpoints to the closest valid value in the time series (this is required to
be able to interpolate over these points). My first attempt uses loops, but
is obviously quite slow, and is quite probably the ugliest piece of code
I've seen in a long while! So please bear with me for a minute ;)

t_axis = numpy.arange(A.shape[0])
for ix in xrange(2000):
    for iy in xrange(2000):
        # get the time series for the current pixel
        t_series = A[:, ix, iy ]
        # Locate valid samples
        valid = t_series <> INVALID_VALUE
        # if the edges are invalid, fill in with the closest valid value.
        if not valid [0]:
            t_series[0] = t_series[valid][0]
            valid[0] = True
        if not valid [-1]:
            t_series[-1] = t_series[valid][-1]
            valid[-1] = True
        # do 1D interpolation on t_series, returning the interpolated values
where valid isn't True. Use eg scipy.interpolate.Univariate.Spline
        sp = scipy.interpolate.UnivariateSpline ( t_axis[valid],
t_series[valid] )
        it_series = numpy.where ( valid, t_series, sp ( t_axis ) )
        # Pack the data back into a 3D array

I don't know whether standard scipy (or ndimage) has stuff to solve this
problem in a reasonable way. Any ideas on how to improve on this much
welcomed!

Thanks!
Jose
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20100414/b67cf964/attachment.html>