python - 在 python 中计算体积或表面积的好算法

Question

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

有谁知道哪里可以找到执行上述操作的工具包或包？

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

Edit1:这是一些(较差的)样本 data .这比典型的球体小得多。当我可以生成它时，我会添加更好的数据!它位于(自身)tumorBrain.tumor 中。

Answer 1

一种选择是使用 VTK。 (我将在这里使用 tvtk python 绑定(bind)...)

至少在某些情况下，获取等值面内的面积会更准确一些。

此外，就表面积而言，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 %