Re: [Numpy-discussion] broadcasting with numpy.interp

On Wed, Nov 24, 2010 at 3:16 PM, Friedrich Romstedt < friedrichromstedt@gmail.com> wrote:

2010/11/16 greg whittier :

...

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

Thanks so much for the response! Sorry I didn't respond earlier. I put it aside until I found time to try and understand part 2 of your response and forgot about it. I'm not really looking for 2d interpolation at the moment, but I can see needing it in the future. Right now, I just want to interpolate along one of the three axes. I think what you're suggesting might work for 1d or 2d depending on how you find the nearest neighbors. What routine would you use? Also, when you say "unmasked" do you mean literally using masked arrays? Thanks, Greg