如何最好地在 python 中创建 Lorenz Attractor 的填充体积以进行体积渲染
How best to create a filled volume of Lorenz Attractor in python for volume rendering
我正在尝试创建一个 3D 阵列,然后对像 Lorenz Attractor 这样的奇怪吸引子(在其他软件或体积渲染包中)执行体积渲染。从数据点绘制吸引子并提供一个值来分配颜色并在 matplotlib 中可视化非常容易。
但是我想要一个填充的体积阵列。我尝试过像 griddata 这样的插值方法,但它没有给出预期的结果。我设想的是这样的:
来自维基百科页面。
这是我尝试过的方法,但如果您在简单的查看器中打开结果,它看起来不太好。我在想,也许只在构成 x、y、z 数组的点之间进行插值……玩了几个小时后我有点迷路了。我认为我需要的是取点并进行某种插值或填充到数组中,这里我调用 interp_im。然后可以在体积渲染中查看。非常感谢任何帮助!
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy.interpolate import griddata
from scipy.interpolate import LinearNDInterpolator
from skimage.external import tifffile
rho = 28.0
sigma = 10.0
beta = 8.0 / 3.0
def f(state, t):
x, y, z = state # Unpack the state vector
return sigma * (y - x), x * (rho - z) - y, x * y - beta * z # Derivatives
state0 = [1.0, 1.0, 1.0]
t = np.arange(0.0, 40.0, 0.01) #t = np.arange(0.0, 40.0, 0.01)
states = odeint(f, state0, t)
# shift x,y,z positions to int for regular image volume
x = states[:, 0]
y = states[:, 1]
z = states[:, 2]
x_min = x.min()
y_min = y.min()
z_min = z.min()
states_int = states + [abs(x_min),abs(y_min),abs(z_min)] + 1
states_int = states_int * 10
states_int = states_int.astype(int)
# values will be in order of tracing for color
values = []
for i,j in enumerate(states_int):
values.append(i*10)
values = np.asarray(values)
fig = plt.figure()
ax = fig.gca(projection='3d')
sc = ax.scatter(states_int[:, 0], states_int[:, 1], states_int[:, 2],c=values)
plt.colorbar(sc)
plt.draw()
plt.show()
#print(x.shape, y.shape, z.shape, values.shape)
#Interpolate for volume rendering
x_ = np.linspace(0,999,500)
y_ = np.linspace(0,999,500)
z_ = np.linspace(0,999,500)
xx,yy,zz = np.meshgrid(x_,y_,z_, sparse = True)
#
# X = states_int.tolist()
#
interp_im = griddata(states_int, values, (xx,yy,zz), method='linear')
interp_im = interp_im.astype(np.uint16)
np.save('interp_im.npy', interp_im)
tifffile.imsave('LorenzAttractor.tif', interp_im)
您的数据在卷中,只是像素化了。如果你模糊体积,例如用高斯,你会得到更有用的东西。例如:
from scipy import ndimage
vol = np.zeros((512, 512, 512), dtype=states_int.dtype)
# add data to vol
vol[tuple(np.split(states_int, vol.ndim, axis=1))] = values[:, np.newaxis]
# apply gaussian filter, sigma=5 in this case
vol = ndimage.gaussian_filter(vol, 5)
然后我会使用类似 napari 的东西来查看 3D 数据:
import napari
with napari.gui_qt():
napari.view_image(v)
为了使音量更平滑,您可能需要减小积分步长。
我正在尝试创建一个 3D 阵列,然后对像 Lorenz Attractor 这样的奇怪吸引子(在其他软件或体积渲染包中)执行体积渲染。从数据点绘制吸引子并提供一个值来分配颜色并在 matplotlib 中可视化非常容易。
但是我想要一个填充的体积阵列。我尝试过像 griddata 这样的插值方法,但它没有给出预期的结果。我设想的是这样的:
来自维基百科页面。
这是我尝试过的方法,但如果您在简单的查看器中打开结果,它看起来不太好。我在想,也许只在构成 x、y、z 数组的点之间进行插值……玩了几个小时后我有点迷路了。我认为我需要的是取点并进行某种插值或填充到数组中,这里我调用 interp_im。然后可以在体积渲染中查看。非常感谢任何帮助!
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy.interpolate import griddata
from scipy.interpolate import LinearNDInterpolator
from skimage.external import tifffile
rho = 28.0
sigma = 10.0
beta = 8.0 / 3.0
def f(state, t):
x, y, z = state # Unpack the state vector
return sigma * (y - x), x * (rho - z) - y, x * y - beta * z # Derivatives
state0 = [1.0, 1.0, 1.0]
t = np.arange(0.0, 40.0, 0.01) #t = np.arange(0.0, 40.0, 0.01)
states = odeint(f, state0, t)
# shift x,y,z positions to int for regular image volume
x = states[:, 0]
y = states[:, 1]
z = states[:, 2]
x_min = x.min()
y_min = y.min()
z_min = z.min()
states_int = states + [abs(x_min),abs(y_min),abs(z_min)] + 1
states_int = states_int * 10
states_int = states_int.astype(int)
# values will be in order of tracing for color
values = []
for i,j in enumerate(states_int):
values.append(i*10)
values = np.asarray(values)
fig = plt.figure()
ax = fig.gca(projection='3d')
sc = ax.scatter(states_int[:, 0], states_int[:, 1], states_int[:, 2],c=values)
plt.colorbar(sc)
plt.draw()
plt.show()
#print(x.shape, y.shape, z.shape, values.shape)
#Interpolate for volume rendering
x_ = np.linspace(0,999,500)
y_ = np.linspace(0,999,500)
z_ = np.linspace(0,999,500)
xx,yy,zz = np.meshgrid(x_,y_,z_, sparse = True)
#
# X = states_int.tolist()
#
interp_im = griddata(states_int, values, (xx,yy,zz), method='linear')
interp_im = interp_im.astype(np.uint16)
np.save('interp_im.npy', interp_im)
tifffile.imsave('LorenzAttractor.tif', interp_im)
您的数据在卷中,只是像素化了。如果你模糊体积,例如用高斯,你会得到更有用的东西。例如:
from scipy import ndimage
vol = np.zeros((512, 512, 512), dtype=states_int.dtype)
# add data to vol
vol[tuple(np.split(states_int, vol.ndim, axis=1))] = values[:, np.newaxis]
# apply gaussian filter, sigma=5 in this case
vol = ndimage.gaussian_filter(vol, 5)
然后我会使用类似 napari 的东西来查看 3D 数据:
import napari
with napari.gui_qt():
napari.view_image(v)
为了使音量更平滑,您可能需要减小积分步长。