python求表面积代码_用于计算python中的体积或表面积的良好算法

我正在尝试计算3D numpy数组的体积(或表面积).在许多情况下,体素是各向异性的,并且我在每个方向上具有像素到厘米的转换因子.

有没有人知道找到工具包或包来做上述的好地方？

现在,我有一些内部代码,但我希望在准确性方面升级到更具工业实力的东西.

编辑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 %