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**

`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**

`scipy.ndimage.gaussian_filter1d(data, window, axis=2)`

?
**#4**

`(x, y, bnd-wind:bnd+wind)`

and after fitting, return H[:,:,bnd]
**#5**

`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**