使用 numpy/scipy 从 3D 数组计算等值面

Question

我有一个 3D numpy 数组，其中包含给定函数的值。我想计算一个 2D 等值面，或一组表示此函数的某些值的等值面。

在这种特殊情况下，可以独立处理 3D 数组的每个 1D 列 (column = myarray[i, j, :])。所以我想知道的是函数等于某个值的最后一个索引位置（二维数组），比如myvalue。

一些（慢）代码来举例说明：

# myarray = 3D ndarray
import numpy as np
from scipy import interpolate

result = np.zeros(nx, ny)
z_values = np.arange(nz)

for i in range(nx):
    for j in range(ny):
        f = interpolate.interp1d(my_array[i, j], z_values)
        result[i, j] = f(myvalue)

我知道这可以通过np.ndenumerate 和其他技巧来加快速度，但想知道是否已经有一种更简单的方法来制作这种等值面。我在ndimage 或其他库中找不到任何东西。我知道 mayavi2 和 vtk 有很多工具可以处理等值面，但我的目标不是可视化——我想对这些等值面值进行计算，而不是显示它们。另外，vtk 的很多等值面方法似乎都涉及到多边形等，而我需要的只是每个等值面值的二维位置数组。

Answer 1

仅使用numpy，您可以使用argsort、sort、take 和适当的数组操作获得一个很好的解决方案。下面的函数使用加权平均来计算等值面：

def calc_iso_surface(my_array, my_value, zs, interp_order=6, power_parameter=0.5):
    if interp_order < 1: interp_order = 1
    from numpy import argsort, take, clip, zeros
    dist = (my_array - my_value)**2
    arg = argsort(dist,axis=2)
    dist.sort(axis=2)
    w_total = 0.
    z = zeros(my_array.shape[:2], dtype=float)
    for i in xrange(int(interp_order)):
        zi = take(zs, arg[:,:,i])
        valuei = dist[:,:,i]
        wi = 1/valuei
        clip(wi, 0, 1.e6, out=wi) # avoiding overflows
        w_total += wi**power_parameter
        z += zi*wi**power_parameter
    z /= w_total
    return z

此解决方案不处理有多个z 对应于my_value 的情况。以下代码中给出了构建等值面的应用示例：

from numpy import meshgrid, sin, cos, pi, linspace
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
dx = 100; dy =  50; dz = 25
nx = 200; ny = 100; nz = 100
xs = linspace(0,dx,nx)
ys = linspace(0,dy,ny)
zs = linspace(0,dz,nz)
X,Y,Z = meshgrid( xs, ys, zs, dtype=float)
my_array = sin(0.3*pi+0.4*pi*X/dx)*sin(0.3*pi+0.4*pi*Y/dy)*(Z/dz)

fig = plt.figure()
ax = fig.gca(projection='3d')

z = calc_iso_surface( my_array, my_value=0.1, zs=zs, interp_order=6 )
ax.plot_surface(X[:,:,0], Y[:,:,0], z, cstride=4, rstride=4, color='g')

z = calc_iso_surface( my_array, my_value=0.2, zs=zs, interp_order=6 )
ax.plot_surface(X[:,:,0], Y[:,:,0], z, cstride=4, rstride=4, color='y')

z = calc_iso_surface( my_array, my_value=0.3, zs=zs, interp_order=6 )
ax.plot_surface(X[:,:,0], Y[:,:,0], z, cstride=4, rstride=4, color='b')

plt.ion()
plt.show()

您还可以使用不同的插值函数。请参阅下面的一个示例，该示例取两个最接近的zs 的平均值：

def calc_iso_surface_2(my_array, my_value, zs):
    '''Takes the average of the two closest zs
    '''
    from numpy import argsort, take
    dist = (my_array - my_value)**2
    arg = argsort(dist,axis=2)
    z0 = take(zs, arg[:,:,0])
    z1 = take(zs, arg[:,:,1])
    z = (z0+z1)/2
    return z