Python 中计算体积或表面积的良好算法

Question

我正在尝试计算三维numpy数组的体积（或表面积）。在很多情况下，体素是各向异性的，并且我有每个方向上的像素到厘米转换因子。

请问是否有人知道一个好的工具包或软件来完成上述任务？

目前我有一些内部代码，但我正在寻找更准确、更强大的工具。

编辑1：这里有一些（糟糕的）data 示例数据。这比典型的球小得多。当我能够生成更好的数据时，我将添加它！它在(self.)tumorBrain.tumor 中。

- tylerthemiler

---

你需要更清楚地定义你的意思。如果你的数据实际上是一个三维数组，那么整个网格占据的体积为 nx * dx * ny * dy * nz * dz，但我很确定你不是这个意思... 你是想要等值面内部的体积吗？ - Joe Kington

我认为你是正确的。这是一个 X x Y x Z 二进制数组，我想要计算出它二进制部分周围的所有内容的体积。它通常（但不总是）呈球形状。 - tylerthemiler

这听起来很有趣，你能发一些示例数据的链接吗？只需使用pickle.dump保存NumPy数组即可。 - wim

我为你添加了一些示例数据。它们不是很好，如果你愿意的话，我稍后会添加更多！ - tylerthemiler

@tylerthemiller - 只是提醒一下：如果你要pickle numpy数组，请确保指定除默认ascii协议以外的其他协议。否则，你会得到一个巨大的文件（就像你的例子一样）。例如：pickle.dump(data, f, pickle.HIGHEST_PROTOCOL) 或者，你可以直接使用numpy.save。此外，我们无法加载你的数据，因为你pickled某种dataStructures对象，而没有你的代码，我们无法unpickle它。尝试只pickle numpy数组和一个(dx, dy, dz)元组。 - Joe Kington

由于您正在处理医学图像，因此您应该查看invesalius - 一款用Python编写的免费3D医学成像软件（由巴西政府赞助）- 它可能具有您需要的功能，或者允许您进行编码。http://www.softwarepublico.gov.br/ver-comunidade?community_id=626732 - jsbueno

Answer 1

有一个选择是使用VTK。(这里我将使用tvtk的Python绑定...)

至少在某些情况下，获取等值面内部的区域会更加准确。

另外，就表面积而言，tvtk.MassProperties也可以计算表面积。代码中的mass.surface_area即为其表面积。

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 %

- Joe Kington

---

这看起来很棒！等我的老板回来我会试一下（我甚至没有管理员权限，所以无法安装vtk :( ）。 - tylerthemiler

如果你遇到了 traits.trait_errors.TraitError: The 'input' trait of an ImageMarchingCubes instance is 'read only'. 错误，请查看我下面的答案更新（针对 tvtk 更改的代码更新）。 - SkyLeach

Answer 2

如果您尝试使用上面Joe的答案，您会得到：

traits.trait_errors.TraitError: ImageMarchingCubes实例的“input”特性为“只读”。

这里是所需更改以及显示如何修复它的差异。

import numpy as np
from tvtk.api import tvtk
from tvtk.common import configure_input


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()
    configure_input(iso, grid)  # <== will work
    # iso = tvtk.ImageMarchingCubes(input=grid)
    mass = tvtk.MassProperties()
    configure_input(mass, iso)
    # mass = tvtk.MassProperties(input=iso.output)
    return mass.volume


if __name__ == '__main__':
    main()

2a3,4
> from tvtk.common import configure_input
> 
6c8
<     # x,y,and z will range from -2 to 2, but with a 
---
>     # x,y,and z will range from -2 to 2, but with a
8c10
<     x,y,z = np.mgrid[-2:2:20j, -2:2:15j, -2:2:5j]
---
>     x, y, z = np.mgrid[-2:2:20j, -2:2:15j, -2:2:5j]
11c13,14
<     dx, dy, dz = [np.diff(it, axis=a)[0,0,0] for it, a in zip((x,y,z),(0,1,2))]
---
>     dx, dy, dz = [np.diff(it, axis=a)[0, 0, 0] for it, a in zip(
>         (x, y, z), (0, 1, 2))]
24,26c27,30
<     print 'Ideal volume', ideal_volume
<     print 'Coarse approximation', coarse_volume, 'Error', coarse_error, '%'
<     print 'VTK approximation', est_volume, 'Error', vtk_error, '%'
---
>     print('Ideal volume', ideal_volume)
>     print('Coarse approximation', coarse_volume, 'Error', coarse_error, '%')
>     print('VTK approximation', est_volume, 'Error', vtk_error, '%')
> 
28c32
< def vtk_volume(data, spacing=(1,1,1), origin=(0,0,0)):
---
> def vtk_volume(data, spacing=(1, 1, 1), origin=(0, 0, 0)):
31c35
<     grid.point_data.scalars = data.T.ravel() # It wants fortran order???
---
>     grid.point_data.scalars = data.T.ravel()  # It wants fortran order???
35,36c39,44
<     iso = tvtk.ImageMarchingCubes(input=grid)
<     mass = tvtk.MassProperties(input=iso.output)
---
>     iso = tvtk.ImageMarchingCubes()
>     configure_input(iso, grid)  # <== will work
>     # iso = tvtk.ImageMarchingCubes(input=grid)
>     mass = tvtk.MassProperties()
>     configure_input(mass, iso)
>     # mass = tvtk.MassProperties(input=iso.output)
39c47,49
< main()
---
> 
> if __name__ == '__main__':
>     main()

• 运行2to3以便在python 3中运行
• 根据autopep8的要求修复代码以符合PEP8标准（包括语法、长度和间距更改）
• 由于TVTK中的这些更改（Github代码更改），需要导入configure_imput
• 修改ImageMarchingCubes构造函数中的input= kwarg
• 修改MassProperties构造函数中的input= kwargs
• 将对main()的调用封装在直接调用中（以防止被导入时执行）#最佳实践

- SkyLeach

Answer 3

这些体素会是相当简单的，规则的多面体，对吧？计算每个体素的体积并将它们相加。

- Charlie Martin

---

主要问题在于z方向太粗糙了。我们采用一种外壳方法来对切片之间进行平均，但这还不够好。 - tylerthemiler