在python中计算体积或表面积的好算法
我想要计算3D numpy数组的体积(或表面积)。在很多情况下,体素是各向异性的,我在每个方向都有像素到厘米的转换因子。在python中计算体积或表面积的好算法
有没有人知道一个好的地方找到一个工具包或包做上述?? ??
现在,我有一些内部代码,但我期待升级到更精确的工业实力。
编辑1:这是一些(差)样品data。这比典型的球体小得多。我可以在生成它时添加更好的数据!它在(自身)肿瘤脑肿瘤。
一个选项是使用VTK
。 (我将在这里使用tvtk
python绑定...)
至少在某些情况下,获取等值面内的区域会更精确一些。
另外,就表面积而言,tvtk.MassProperties
也计算表面积。它是mass.surface_area
(下面的代码中有mass
对象)。
import numpy as np
from tvtk.api import tvtk
def main():
# Generate some data with anisotropic cells...
# x,y,and z will range from -2 to 2, but with a
# different (20, 15, and 5 for x, y, and z) number of steps
x,y,z = np.mgrid[-2:2:20j, -2:2:15j, -2:2:5j]
r = np.sqrt(x**2 + y**2 + z**2)
dx, dy, dz = [np.diff(it, axis=a)[0,0,0] for it, a in zip((x,y,z),(0,1,2))]
# Your actual data is a binary (logical) array
max_radius = 1.5
data = (r <= max_radius).astype(np.int8)
ideal_volume = 4.0/3 * max_radius**3 * np.pi
coarse_volume = data.sum() * dx * dy * dz
est_volume = vtk_volume(data, (dx, dy, dz), (x.min(), y.min(), z.min()))
coarse_error = 100 * (coarse_volume - ideal_volume)/ideal_volume
vtk_error = 100 * (est_volume - ideal_volume)/ideal_volume
print 'Ideal volume', ideal_volume
print 'Coarse approximation', coarse_volume, 'Error', coarse_error, '%'
print 'VTK approximation', est_volume, 'Error', vtk_error, '%'
def vtk_volume(data, spacing=(1,1,1), origin=(0,0,0)):
data[data == 0] = -1
grid = tvtk.ImageData(spacing=spacing, origin=origin)
grid.point_data.scalars = data.T.ravel() # It wants fortran order???
grid.point_data.scalars.name = 'scalars'
grid.dimensions = data.shape
iso = tvtk.ImageMarchingCubes(input=grid)
mass = tvtk.MassProperties(input=iso.output)
return mass.volume
main()
这产生了:
Ideal volume 14.1371669412
Coarse approximation 14.7969924812 Error 4.66731094565 %
VTK approximation 14.1954890878 Error 0.412544796894 %
这看起来很棒!当我的老板回来时我会试一试(我甚至没有管理员权限,所以我不能安装vtk :()。 – tylerthemiler 2012-01-19 00:08:41
体素将是相当简单,规则的多面体,不是吗?计算每一个的数量并对它们进行求和。
主要问题是z方向太粗糙。我们采用一种贝壳方法来平均切片,但这不够好。 – tylerthemiler 2012-01-18 00:54:52
你需要定义你的意思有点更清晰。如果你的数据实际上是一个3D数组,那么整个网格所占据的体积是'nx * dx * ny * dy * nz * dz',但我敢肯定你并不是这个意思......你想要体积内的等值面? – 2012-01-18 00:56:49
我认为你是对的。它是一个X×Y×Z的二进制数组,我想要包含在其二进制部分周边的所有内容。它通常(但不总是)球形。 – tylerthemiler 2012-01-18 01:01:25
这听起来很有趣,你可以发布一些示例数据的链接吗?用'pickle.dump'保存numpy数组 – wim 2012-01-18 01:14:56