[Numpy-discussion] broadcasting with numpy.interp
Wed Dec 1 13:16:13 CST 2010
On Wed, Nov 24, 2010 at 3:16 PM, Friedrich Romstedt <
> 2010/11/16 greg whittier <email@example.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)
> > for line in range(cube.shape):
> > dead_bands = bands[dead[:, line] == True]
> > good_bands = bands[dead[:, line] == False]
> > for sample in range(cube.shape):
> > # interp returns fp for x < xp 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.
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?
-------------- next part --------------
An HTML attachment was scrubbed...
More information about the NumPy-Discussion