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)

标签： pythonnumpyscipyvolume-renderingstrange-attractor

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)

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