使用numpy/scipy从3D阵列计算等值面 [英] Using numpy/scipy to calculate iso-surface from 3D array
问题描述
我有一个3D numpy数组,其中包含给定函数的值.我想计算一个二维等值面,或一组代表该函数某些值的等值面.
I have a 3D numpy array that contains the values of a given function. I want to calculate a 2D iso-surface, or a set of iso-surfaces that represent certain values of this function.
在这种特定情况下,可以独立处理3D数组的每个1D列(column = myarray[i, j, :]
).因此,我想知道的是函数等于某个值(例如myvalue
)的最后一个索引位置(二维数组).
In this particular case, each 1D column (column = myarray[i, j, :]
) of the 3D array can be treated independently. So what I would like to know are the last index positions (2D array) where the function is equal to a certain value, say myvalue
.
一些(慢速)代码示例:
Some (slow) code to exemplify:
# 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的许多等值面方法似乎都涉及多边形等,而我需要的只是每个等值面值的二维数组.
I know this can be sped up a bit with np.ndenumerate
and other tricks, but was wondering if there is already a simpler way of doing this kind of iso-surface. I couldn't find anything in ndimage
or other libraries. I know that mayavi2 and vtk have a lot of tools to deal with iso-surfaces, but my aim here is not visualisation -- I want to perform calculations on those iso-surface values, not display them. Plus, a lot of the iso-surface methods of vtk seem to involve polygons and the like, and what I need is just a 2D array of positions for each iso-surface value.
推荐答案
仅使用numpy
,您可以使用argsort
,sort
,take
和适当的数组操作来获得良好的解决方案.下面的函数使用加权平均值来计算等值面:
Using only numpy
you can get a good solution using argsort
, sort
, take
and the proper array manipulation. The function below uses a weighted average to compute the iso-surface:
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
此解决方案无法处理与my_value
对应的多个z
的情况.以下代码给出了构建下面的等值面的应用示例:
This solution does not handle situations where there is more than one z
corresponding to my_value
. An application example to build the iso-surfaces below is given in the following code:
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
的平均值:
You can also play with different interpolation functions. See below one example that takes the average of the two closest 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
这篇关于使用numpy/scipy从3D阵列计算等值面的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!