计算 3d 绘图的体积

Question

数据来自测量。 plotted data的图片

我尝试使用 trapz 两次，但得到错误代码：“ValueError：操作数无法与形状一起广播 (1,255) (256,531)” x 有 256 个点，y 有 532 个点，Z 也是一个 256 x 532 长度的二维数组。代码如下：

import numpy as np

 img=np.loadtxt('focus_x.txt')
 m=0
 m=np.max(img)
 Z=img/m 

 X=np.loadtxt("pixelx.txt",float)
 Y=np.loadtxt("pixely.txt",float)

 [X, Y] = np.meshgrid(X, Y)
 volume=np.trapz(X,np.trapz(Y,Z))

Answer 1

The docs 声明 trapz 应该这样使用

intermediate = np.trapz(Z, x)
result = np.trapz(intermediate, y)

trapz 正在减少其操作数的维数（默认情况下在最后一个轴上），使用可选一维横坐标数组来确定积分的子区间；它没有使用网格进行操作。

一个完整的例子。

首先，我们使用sympy 计算一个简单双线性函数在矩形域(0, 5) × (0, 7) 上的积分

In [1]: import sympy as sp, numpy as np
In [2]: x, y = sp.symbols('x y')
In [3]: f = 1 + 2*x + y + x*y
In [4]: f.integrate((x, 0, 5)).integrate((y, 0, 7))
Out[4]: 2555/4

现在我们计算积分的梯形逼近（碰巧，逼近对于双线性函数是精确的）——我们需要坐标数组

In [5]: x, y = np.linspace(0, 5, 11), np.linspace(0, 7, 22)

（注意采样在两个方向上是不同的，并且与trapz使用的默认值不同）——我们需要一个网格来计算被积函数，我们需要计算被积函数

In [6]: X, Y = np.meshgrid(x, y)
In [7]: z = 1 + 2*X + Y + X*Y

最终我们计算积分

In [8]: 4*np.trapz(np.trapz(z, x), y)
Out[8]: 2555.0