Faster way to perform point-wise interplation of numpy array?

Faster way to perform point-wise interplation of numpy array?,numpy,performance,curve-fitting,smoothing,Numpy,Performance,Curve Fitting,Smoothing,I have a 3D datacube, with two spatial dimensions and the third being a multi-band spectrum at each point of the 2D image. H[x, y, bands] Given a wavelength (or band number), I would like to extract the 2D image corresponding to that wavelength. This would be simply an array slice like H[:,:,bnd]. Similarly, given a spatial location (i,j) the spectrum at that location is H[i,j]. I would also like to 'smooth' the image spectrally, to counter low-light noise in the spectra. That is for band bn

I have a 3D datacube, with two spatial dimensions and the third being a multi-band spectrum at each point of the 2D image.

H[x, y, bands]

Given a wavelength (or band number), I would like to extract the 2D image corresponding to that wavelength. This would be simply an array slice like

H[:,:,bnd]
. Similarly, given a spatial location (i,j) the spectrum at that location is
H[i,j]
.

I would also like to 'smooth' the image spectrally, to counter low-light noise in the spectra. That is for band

bnd
, I choose a window of size
wind
and fit a n-degree polynomial to the spectrum in that window. With polyfit and polyval I can find the fitted spectral value at that point for band
bnd
.

Now, if I want the whole image of

bnd
from the fitted value, then I have to perform this windowed-fitting at each
(i,j)
of the image. I also want the 2nd-derivative image of
bnd
, that is, the value of the 2nd-derivative of the fitted spectrum at each point.

Running over the points, I could polyfit-polyval-polyder each of the

x*y
spectra. While this works, this is a point-wise operation. Is there some pytho-numponic way to do this faster?


#1

If you do least-squares polynomial fitting to points (x+dx[i],y[i]) for a fixed set of dx and then evaluate the resulting polynomial at x, the result is a (fixed) linear combination of the y[i]. The same is true for the derivatives of the polynomial. So you just need a linear combination of the slices. Look up "Savitzky-Golay filters".

EDITED to add a brief example of how S-G filters work. I haven't checked any of the details and you should therefore not rely on it to be correct.

So, suppose you take a filter of width 5 and degree 2. That is, for each band (ignoring, for the moment, ones at the start and end) we'll take that one and the two on either side, fit a quadratic curve, and look at its value in the middle.

So, if f(x) ~= ax^2+bx+c and f(-2),f(-1),f(0),f(1),f(2) = p,q,r,s,t then we want 4a-2b+c ~= p, a-b+c ~= q, etc. Least-squares fitting means minimizing (4a-2b+c-p)^2 + (a-b+c-q)^2 + (c-r)^2 + (a+b+c-s)^2 + (4a+2b+c-t)^2, which means (taking partial derivatives w.r.t. a,b,c):

  • 4(4a-2b+c-p)+(a-b+c-q)+(a+b+c-s)+4(4a+2b+c-t)=0
  • -2(4a-2b+c-p)-(a-b+c-q)+(a+b+c-s)+2(4a+2b+c-t)=0
  • (4a-2b+c-p)+(a-b+c-q)+(c-r)+(a+b+c-s)+(4a+2b+c-t)=0

or, simplifying,

  • 22a+10c = 4p+q+s+4t
  • 10b = -2p-q+s+2t
  • 10a+5c = p+q+r+s+t

so a,b,c = p-q/2-r-s/2+t, (2(t-p)+(s-q))/10, (p+q+r+s+t)/5-(2p-q-2r-s+2t).

And of course c is the value of the fitted polynomial at 0, and therefore is the smoothed value we want. So for each spatial position, we have a vector of input spectral data, from which we compute the smoothed spectral data by multiplying by a matrix whose rows (apart from the first and last couple) look like [0 ... 0 -9/5 4/5 11/5 4/5 -9/5 0 ... 0], with the central 11/5 on the main diagonal of the matrix.

So you could do a matrix multiplication for each spatial position; but since it's the same matrix everywhere you can do it with a single call to tensordot. So if S contains the matrix I just described (er, wait, no, the transpose of the matrix I just described) and A is your 3-dimensional data cube, your spectrally-smoothed data cube would be numpy.tensordot(A,S).

This would be a good point at which to repeat my warning: I haven't checked any of the details in the few paragraphs above, which are just meant to give an indication of how it all works and why you can do the whole thing in a single linear-algebra operation.


#2

On playing the RTFM game, I found that the np.polyfit function can take an array of y-columns to fit for x[i], y[i]. This speeds up my operation by one dimension because I can process the cube line-by-line. This semi-answers my own question -- but is there a better way?

#3

It sounds like you're just smoothing over one dimension... Why not just call scipy.ndimage.gaussian_filter1d(data, window, axis=2)?

#4

Thanks for that info. The SG filter is indeed a great tool to get rid of data noise. However, this does not really answer my question here because this would still be a point-wise calculation of each spectrum. The savitzky_golay filter takes a 1D array as input. I want to take a cube of shape (x, y, bnd-wind:bnd+wind) and after fitting, return H[:,:,bnd]

#5

I wasn't suggesting that you call anyone's savitzky_golay function. I was suggesting that you compute the relevant coefficients, and then do the whole calculation plane by plane. You can probably do it with a single call to numpy.tensordot, in fact.

#6

I'm afraid I do not understand the suggestion. Could you please elaborate, perhaps with an example? I'm using no matrix operations as of now, just array slicing.