python - 在 Python 中绘制 3d 数组的最有效方法是什么？

Question

在 Python 中绘制 3d 数组的最有效方法是什么？

例如:

volume = np.random.rand(512, 512, 512)

其中数组项表示每个像素的灰度颜色。

---

以下代码运行速度太慢:

import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import matplotlib.pyplot as plt

fig = plt.figure()
ax = fig.gca(projection='3d')
volume = np.random.rand(20, 20, 20)
for x in range(len(volume[:, 0, 0])):
    for y in range(len(volume[0, :, 0])):
        for z in range(len(volume[0, 0, :])):
            ax.scatter(x, y, z, c = tuple([volume[x, y, z], volume[x, y, z], volume[x, y, z], 1]))
plt.show()

Answer 1

为了获得更好的性能，请尽可能避免多次调用 ax.scatter。相反，将所有 x、y、z 坐标和颜色打包到 1D 数组(或列表)，然后调用 ax.scatter 一次:

ax.scatter(x, y, z, c=volume.ravel())

---

问题(在 CPU 时间和内存方面)随 size**3 增长，其中 size 是立方体的边长。

此外，ax.scatter 将尝试渲染所有 size**3 点而不考虑事实上，大多数这些点都被外面的那些点遮住了外壳。

这将有助于减少 volume 中的点数——也许是以某种方式对其进行总结或重采样/插值——在渲染之前。

我们还可以将所需的 CPU 和内存从 O(size**3) 减少到 O(size**2)通过仅绘制外壳:

import functools
import itertools as IT
import numpy as np
import scipy.ndimage as ndimage
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def cartesian_product_broadcasted(*arrays):
    """
    http://stackoverflow.com/a/11146645/190597 (senderle)
    """
    broadcastable = np.ix_(*arrays)
    broadcasted = np.broadcast_arrays(*broadcastable)
    dtype = np.result_type(*arrays)
    rows, cols = functools.reduce(np.multiply, broadcasted[0].shape), len(broadcasted)
    out = np.empty(rows * cols, dtype=dtype)
    start, end = 0, rows
    for a in broadcasted:
        out[start:end] = a.reshape(-1)
        start, end = end, end + rows
    return out.reshape(cols, rows).T

# @profile  # used with `python -m memory_profiler script.py` to measure memory usage
def main():
    fig = plt.figure()
    ax = fig.add_subplot(1, 1, 1, projection='3d')

    size = 512
    volume = np.random.rand(size, size, size)
    x, y, z = cartesian_product_broadcasted(*[np.arange(size, dtype='int16')]*3).T
    mask = ((x == 0) | (x == size-1) 
            | (y == 0) | (y == size-1) 
            | (z == 0) | (z == size-1))
    x = x[mask]
    y = y[mask]
    z = z[mask]
    volume = volume.ravel()[mask]

    ax.scatter(x, y, z, c=volume, cmap=plt.get_cmap('Greys'))
    plt.show()

if __name__ == '__main__':
    main()

但请注意，即使只绘制外壳，也可以使用size=512 我们仍然需要大约 1.3 GiB 的内存。还要注意，即使你有足够的总内存，但由于缺少 RAM，程序使用交换空间，程序的整体速度也会急剧减速。如果您发现自己处于这种情况，那么唯一的解决方案是找到一种更聪明的方法来渲染可接受的图像使用更少的点数，或者购买更多的 RAM。