How best to create a filled volume of Lorenz Attractor in python for volume rendering

I am trying to create a 3D array to then perform volume rendering on (in other software or volume rendering packages) of strange attractor like Lorenz Attractor. It is easy enough to plot the attractor from data points and provide a value to assign a color and visualize in matplotlib for example.

However I would like a filled volume array. I have tried interpolation methods like griddata but it doesn't give the desired result. What I am envisioning is something like:

Which is from the wikipedia page.

Here is what I have tried but if you open the result in a simple viewer it doesn't look great. I am thinking instead maybe doing a interpolation only between the points that make up the x,y,z array... I am a little lost after playing with this for several hours. What I think I need is to take the points and do some sort of interpolation or filling into an array, here I am calling interp_im. This can then be viewed in volume rendering. Any help is greatly appreciated on this!

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)

Your data is in the volume, it is just pixelated. If you blur the volume, for example with a gaussian, you get something much more usable. For example:

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)

I would then use something like napari to view the data in 3D:

import napari

with napari.gui_qt():
    napari.view_image(v)

To make the volume smoother you may want to reduce your integration step size.