首页
学习
活动
专区
圈层
工具
发布
首页
学习
活动
专区
圈层
工具
MCP广场
社区首页 >问答首页 >执行numpy数组逐点内插的更快方法?

执行numpy数组逐点内插的更快方法?
EN

Stack Overflow用户
提问于 2011-04-08 05:23:22
回答 1查看 336关注 0票数 0

我有一个3D datacube,它有两个空间维度,第三个维度是2D图像每个点的多波段光谱。

代码语言:javascript
运行
复制
H[x, y, bands]

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

我还想在光谱上“平滑”图像,以对抗光谱中的弱光噪声。也就是说,对于频带bnd,我选择了一个大小为wind的窗口,并将n次多项式拟合到该窗口中的频谱。使用polyfit和polyval,我可以找到频带bnd在该点的拟合光谱值。

现在,如果我想要来自拟合值的bnd的整个图像,那么我必须在图像的每个(i,j)上执行这个窗口拟合。我还想要bnd的二阶导数图像,即拟合光谱在每个点的二阶导数的值。

在这些点上运行,我可以对每个x*y光谱进行多项式拟合。虽然这是有效的,但这是一个逐点操作。有没有什么方法可以更快地做到这一点?

EN

回答 1

Stack Overflow用户

发布于 2011-04-08 05:33:52

如果对固定dx集的点(x+dxi,yi)进行最小二乘多项式拟合,然后计算x处的结果多项式,则结果是yi的(固定)线性组合。多项式的导数也是如此。所以你只需要切片的线性组合。查找"Savitzky-Golay filters“。

编辑以添加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,这意味着(取偏导数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

或者,简化一下,

  • 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在矩阵的主对角线上。

因此,您可以对每个空间位置执行矩阵乘法;但是,由于每个位置的矩阵都是相同的,因此只需调用一次tensordot即可完成此操作。因此,如果S包含我刚才描述的矩阵(呃,等等,不,我刚刚描述的矩阵的转置),并且A是您的3维数据立方体,那么您的光谱平滑数据立方体就是numpy.tensordot(A,S)

这将是重复我的警告的好地方:我没有检查上面几段中的任何细节,这只是为了给出它是如何工作的,以及为什么你可以在一个单一的线性代数操作中完成所有的事情。

票数 1
EN
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接:

https://stackoverflow.com/questions/5587823

复制
相关文章

相似问题

领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档