更快的方法来执行 numpy 数组的逐点插值？

Question

我有一个 3D 数据立方体，具有两个空间维度，第三个是 2D 图像每个点的多波段光谱。

H[x, y, bands]

给定一个波长（或波段编号），我想提取与该波长对应的二维图像。这只是一个像H[:,:,bnd] 这样的数组切片。同样，给定一个空间位置 (i,j)，该位置的频谱为H[i,j]。

我还想对图像进行光谱“平滑”，以消除光谱中的低光噪声。对于波段bnd，我选择了一个大小为wind 的窗口，并将n 次多项式拟合到该窗口中的频谱。使用 polyfit 和 polyval，我可以找到频带 bnd 的拟合光谱值。

现在，如果我想从拟合值中获得bnd 的整个图像，那么我必须在图像的每个(i,j) 处执行此窗口拟合。我还要bnd的二阶导数图像，也就是每个点的拟合光谱二阶导数的值。

遍历这些点，我可以对每个x*y 光谱进行 polyfit-polyval-polyder。虽然这有效，但这是一个逐点操作。有没有一些 pytho-numponic 方法可以更快地做到这一点？

Answer 1

如果您对一组固定 dx 的点 (x+dx[i],y[i]) 进行最小二乘多项式拟合，然后在 x 处计算结果多项式，则结果是（固定）线性组合y[i]。多项式的导数也是如此。所以你只需要切片的线性组合。查找“Savitzky-Golay 过滤器”。

已编辑添加 S-G 过滤器如何工作的简短示例。我没有检查任何细节，因此你不应该依赖它是正确的。

因此，假设您采用宽度为 5 和度数为 2 的过滤器。也就是说，对于每个波段（暂时忽略开始和结束的波段），我们将在两侧分别取一个和两个，拟合一条二次曲线，看中间的值。

所以，如果 f(x) ~= ax^2+bx+c 并且 f(-2),f(-1),f(0),f(1),f(2) = p,q ,r,s,t 然后我们想要 4a-2b+c ~= p, a-b+c ~= q 等。最小二乘拟合意味着最小化 (4a-2b+c-p)^2 + (a-b+ c-q)^2 + (c-r)^2 + (a+b+c-s)^2 + (4a+2b+c-t)^2，意思是（对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

或者，简化，

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

所以 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)。

当然 c 是拟合多项式在 0 处的值，因此是我们想要的平滑值。所以对于每个空间位置，我们有一个输入光谱数据的向量，我们通过乘以一个矩阵来计算平滑的光谱数据，该矩阵的行（除了第一对和最后一对）看起来像 [0 ... 0 -9/ 5 4/5 11/5 4/5 -9/5 0 ... 0]，中心 11/5 在矩阵的主对角线上。

所以你可以对每个空间位置进行矩阵乘法；但由于它在任何地方都是 same 矩阵，因此只需调用tensordot 即可。所以如果S 包含我刚刚描述的矩阵（呃，等等，不，我刚刚描述的矩阵的转置）并且A 是你的3维数据立方体，你的光谱平滑数据立方体将是numpy.tensordot(A,S)。

这将是重复我的警告的好点：我没有检查上面几段中的任何细节，这只是为了说明它是如何工作的以及为什么你可以这样做整个事情都在一个线性代数运算中。