python - 如何最好地在 python 中创建 Lorenz Attractor 的填充体积以进行体积渲染

Question

我正在尝试创建一个 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)

Answer 1

您的数据在体积中，只是像素化。如果你模糊体积，例如用高斯，你会得到更有用的东西。例如:

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)

为了使音量更平滑，您可能需要减小积分步长。