[Numpy-discussion] broadcasting with numpy.interp

[Friedrich Romstedt]

2010/11/16 greg whittier <gregwh at gmail.com>:
> I'd like to be able to speed up the following code.
>
> def replace_dead(cube, dead):
>   # cube.shape == (320, 640, 1200)
>   # dead.shape == (320, 640)
>   # cube[i,j,:] are bad points to be replaced via interpolation if
> dead[i,j] == True
>
>    bands = np.arange(0, cube.shape[0])
>    for line in range(cube.shape[1]):
>        dead_bands = bands[dead[:, line] == True]
>        good_bands = bands[dead[:, line] == False]
>        for sample in range(cube.shape[2]):
>            # interp returns fp[0] for x < xp[0] and fp[-1] for x > xp[-1]
>            cube[dead_bands, line, sample] = \
>                np.interp(dead_bands,
>                          good_bands,
>                          cube[good_bands, line, sample])

I assume you just need *some* interpolation, not that specific one?
In that case, I'd suggest the following:

1)  Use a 2d interpolation, taking into account all nearest neighbours.
2)  For this, use a looped interpolation in this nearest-neighbour sense:
    a)  Generate sums of all unmasked nearest-neighbour values
    b)  Generate counts for the nearest neighbours present
    c)  Replace the bad values by the sums divided by the count.
    d)  Continue at (a) if there are bad values left

Bad values which are neighbouring each other (>= 3) need multiple
passes through the loop.  It should be pretty fast.

If this is what you have in mind, maybe we (or I) can make up some code.

Friedrich