[Numpy-discussion] broadcasting with numpy.interp

[greg whittier]

Hi all,

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])

Removing that last loop via some sort of broadcasting seems like it
should be possible, but that doesn't seem to work with interp.  While
interp allows the x-coordinates of interpolation points to be a
multi-dimensional array, it expects the x- and y-coordinates of the
data points to be 1-d arrays.  Any suggestions for speeding this up?

Thanks,
Greg