3D Visualization¶
Here I discuss my learning about making 3D visualizations of scientific simulation data. There are two goals: interactive visualization to understand the data, and production of beautiful movies and plots for publication/press coverage etc. $\newcommand{\d}{\mathrm{d}}\newcommand{\vect}[1]{\vec{#1}}\newcommand{\abs}[1]{\lvert#1\rvert}\newcommand{\I}{\mathrm{i}}\newcommand{\braket}[1]{\langle#1\rangle}\newcommand{\ket}[1]{\left|#1\right\rangle}\newcommand{\bra}[1]{\left\langle#1\right|}\newcommand{\pdiff}[2]{\frac{\partial #1}{\partial #2}}\newcommand{\diff}[3][]{\frac{\d^{#1}#2}{\d{#3}{}^{#1}}}$
Tools¶
Here is a short list of some relevant tools:
-
matplotlib: Although the 3D visualization aspects are very primative and slow, one can often get a good understanding of a 3D data-set using the 2D primatives. For example, a collection of 2D projections or slices can often give a good idea of how a simulation is evolving. Contour plots or colour can be used to visualize a third dimension. -
MayaVi2: This is a VTK-based visualization toolkit that can be used for interactive 3D visualization. There is a GUI where you can manually edit and interact with the scene, as well as a rudimentary command-line system for automatically generating visualizations. One problem with all VTK visualizations is that they expect the aspect ratio of the various coordinates to be the same (it was designed to visualize objects in real space, like a CAT-scan of a head). I have found it challenging to visualize data with very different scales in a natural way. -
Blender: Supposedly a very good piece of software for constructing scenes i.e. adding lighting sources, textures, etc. I am hoping to be able to insert data into a scene to generate good-looking visualizations. -
YafaRay,Luxrender: These are ray-traceing engines that integrate with Blender for rendering a scene.
Sample Problem: Vortex Ring in a Tube¶
For a sample problem, we shall try visualizing the motion of a vortex ring in a harmonic trap based on the Thomas-Fermi (TF) and Logarithmic approximations of Pitaevskii. We take the GPE to have the form:
$$ \hbar\I\dot{\psi} = \frac{-\hbar^2\nabla^2}{2M}\psi + (V-\mu[\rho])\psi $$where $M$ is the effective particle mass ($M=2m$ where $m$ is ther fermionic mas for bosonic dimers in the UFG, or $M=m$ for real bosons) and $\mu[\rho]$ is the equation of state, which depends non-linearly on the density $\rho$.
The motion of a vortex ring is based on the following relationships for the energy and momentum of a vortex ring with radius $R$ at position $Z$ in a trap with background density $\rho$:
$$ P_R \approx \frac{2\pi\hbar}{M}\int_0^{R} \rho(r, Z) 2\pi r \d{r}\\ E_R \approx \frac{2\pi^2\hbar^2}{M^2}R\rho(R,Z)\ln\frac{R}{h}, $$and $h=k_F^{-1} = \hbar/\sqrt{2M\mu[\rho]}$ is the healing length, which may be deduced by equating the kinetic energy with the potential energy:
$$ \frac{\hbar^2 k_F}{2M} = \mu[\rho]. $$This sets the length-scale for the size of the vortex core. In our phenomenological model here, we shall assume that a vortex core has the following approximate structure (see e.g. Fetter:2009)
$$ \psi = \frac{x + \I y}{\sqrt{2 + k_F^2(x^2+y^2)}}. $$The TF approximation is to approximate the density $\rho$ with the TF density
$$ \rho_{TF}(r, z) = \rho[\mu - V(r, z)] = \rho\left[ \mu\left(1 - \frac{r^2}{R_\perp^2} - \frac{z^2}{R_z^2}\right)\right] $$where we have specialized to a Harmonic trap with TF radii $R_\perp$ and $R_z$ and where $\mu[\rho] = \mathcal{E}'(\rho)$ is the (inverse) equation of state. The Logarithmic approximation assumes that one can treat $L=\ln R/\xi$ as a constant.
The velocity of the ring is then given by:
$$ V = \dot{Z} = \left.\pdiff{E_R}{P_R}\right|_{Z} = \frac{(\partial E_R/\partial R)_Z} {(\partial P_R/\partial R)_Z} = \frac{\hbar L}{2M}\left(\frac{1}{R} + \partial_R\ln\rho\right) $$We now introduce the dimensionless coordinates:
$$ X = \frac{Z}{R_z}, \qquad Y = \frac{R}{R_\perp}, \qquad n = \frac{\rho}{\rho_0}, \qquad \tau = \frac{t}{T_0}, \qquad f = \frac{E}{E_0},\qquad \tilde{\mu} = \frac{\mu}{\mu_0}\\ \rho_0 = \rho[\mu], \qquad T_0 = \frac{2MR_\perp R_z}{\hbar L}, \qquad E_0 = \frac{2\pi^2\hbar^2R_\perp \rho_0 L}{M^2}, \qquad \mu_0 = \mu[\rho_0] = \frac{\hbar^2k_F^2}{2M} $$so that the equations become:
$$ \dot{X} = \frac{1}{Y} + \partial_{Y} \ln n, \qquad f(X,Y) = Y n(X, Y) = Y n(1-X^2-Y^2) = C = \text{const}. $$From these, the following system can be derived:
$$ \dot{X} = \frac{1}{Y} - \frac{2Y^2n'}{C},\qquad \dot{Y} = \frac{2XYn'}{C}. $$The equation of state for the systems of interest are
$$ n(1-X^2-Y^2) = \begin{cases} (1-X^2-Y^2) & \text{GPE}\\ (1-X^2-Y^2)^{3/2} & \text{UFG}. \end{cases} $$Physically, these correspond to equations of state:
$$ \rho(\mu) = \begin{cases} \frac{\mu}{g} & \text{GPE}\\ \frac{k_F^3}{3\pi^2} = \frac{(2m\mu)^{3/2}}{3\pi^2\hbar^3} & \text{UFG} \end{cases} $$%%file vortex_ring.py
import math
import numpy as np
import scipy.integrate
class VortexRing(object):
# Units terms of k_F
kF_0 = 1.0
hbar = 1.0
m = 1.0
micron = 2.8334800990033500879
def __init__(self, lam=6.0, # Aspect ratio
model='UFG', a=1.0, xi=0.37,
x0=0.0, y0=0.2):
self.R_z = 265.214532348242562 * self.micron
self.R_perp = self.R_z/lam
self.L = 5.0 # Logarithmic factor
if model == 'GPE':
self.M = self.m
self.n = lambda tmu:tmu
self.dn = lambda tmu:1.0
self.tmu = lambda n: n
self.a = a
self.g = 4*np.pi*a
self.mu_0 = (self.hbar*self.kF_0)**2 / 2.0/self.m
self.rho_0 = self.mu_0 / self.g
elif model == 'UFG':
self.M = 2.0*self.m
self.n = lambda tmu: tmu**(3./2.)
self.dn = lambda tmu: 3./2.*np.sqrt(tmu)
self.tmu = lambda n: n**(2./3.)
self.rho_0 = self.kF_0**3 / (3.0*np.pi**2)
eF_0 = (self.hbar*self.kF_0)**2/2/self.m
self.xi = xi
self.mu_0 = xi * eF_0
else:
raise NotImplementedError
self.T_0 = 2*self.M * self.R_perp*self.R_z / (self.hbar * self.L)
self.E_0 = 2 * (np.pi*self.hbar/self.M)**2 * self.R_perp * self.rho_0 * self.L
self.xy0 = [x0, y0]
self._init_ode()
def f(self, t, y):
x, y = y
dn = self.dn(1-x**2-y**2)
return [1/y - 2*y**2*dn/self._c,
2*x*y*dn/self._c]
def kF(self, rho):
mu = self.mu_0*self.tmu(rho/self.rho_0)
return np.sqrt(2*self.m*mu)/self.hbar
def get_zr(self, t):
if t == self._ode.t:
y = self._ode.y
else:
y = self._ode.integrate(t)
self._ode.set_initial_value(y, t=t)
return y * [self.R_z, self.R_perp]
def _init_ode(self):
x0, y0 = self.xy0
self._c = y0*self.n(1-y0**2-x0**2)
self._ode = scipy.integrate.ode(self.f)
self._ode.set_initial_value(self.xy0, t=0.0)
def get_trajectory(self, t, N=100):
ts = np.linspace(0, t, N)
zrs = np.array([self.get_zr(_t) for _t in ts])
return ts, zrs
def psi(self, z, r, zr=None):
"""We just use the TF density with a vortex ring imprinted
at position ``."""
degen = self.M/self.m
X = z/self.R_z
R = r/self.R_perp
tmu = (1 - R**2 - X**2)
n = np.where(tmu > 0, self.n(tmu), 0.0)
rho = self.rho_0 * n
psi = np.sqrt(rho)/degen
if zr is not None:
z_0, r_0 = zr
w = z - z_0 + (r - r_0)*1j
kF2 = self.kF(rho)**2 + np.finfo(float).tiny
psi = psi * (w/np.sqrt(2/kF2 + abs(w)**2))
return psi
def grid3(self, Nx):
"""Return a 3D grid."""
# Note, MayaVi needs equally spaced data or else one
# needs to futz with extents which is a pain.
x_max = self.R_z
r_max = self.R_perp
dx = 2*x_max/Nx
Nr = 2*r_max/dx
x, y, z = np.meshgrid(np.linspace(-x_max, x_max, Nx),
np.linspace(-r_max, r_max, Nr),
np.linspace(-r_max, r_max, Nr),
sparse=True, indexing='ij')
return x, y, z
def get_data(self, t, Nx=100):
x, y, z = self.grid3(Nx=Nx)
r = np.sqrt(y**2+z**2)
psi = self.psi(x, r, zr=self.get_zr(t))
ks = [2*np.pi * np.fft.fftfreq(_x.size, d=np.diff(_x.ravel())[0]).reshape(_x.shape)
for _x in [x, y, z]]
vs = [(psi.conj()*np.fft.ifft(_k*np.fft.fft(psi, axis=_n), axis=_n)).imag
for _n, _k in enumerate(ks)]
return x, y, z, psi, vs
def get_psi(self, N, zr=None):
z_max = self.R_z
z = np.linspace(-z_max, z_max, N)[:,None]
r = np.linspace(0, self.R_perp, N//2)[None,:]
psi = self.psi(z=z, r=r, zr=zr)
return z, r, psi
%pylab inline
import vortex_ring;reload(vortex_ring)
from vortex_ring import VortexRing
v = VortexRing(model='UFG', lam=2)
ts, zrs = v.get_trajectory(t=1.0, N=10000)
plot(zrs[:,0], zrs[:,1])
Here is an example of plotting slices and projections of the density:
from mmf.utils import mmf_plot as mmfplt
def plot_2d(x, r, psi):
plt.subplot(221)
mmfplt.imcontourf(x, r, abs(psi)**2, aspect=1)
plt.subplot(221)
mmfplt.imcontourf(x, r, abs(psi)**2, aspect=1)
z, r, psi = v.get_psi(100, zr=[30,40])
z, r, psi0 = v.get_psi(100)
figure(figsize=(20,10))
plot_2d(z, r, psi)
figure(figsize=(20,10))
plot_2d(z, r, abs(psi0)**2*abs(abs(psi0)**2-abs(psi)**2))
Gridded Data Format (GDF)¶
The gridded data format (GDF) is a simple HDF5 format for specifying volumetric data:
#f = h5py.File('RD0005.gdf')
f = h5py.File('ring.hd5')
f["/simulation_parameters"].keys()
f.close()
import h5py
!rm -f ring.hd5
rho = abs(psi.T)**2
with h5py.File('ring.hd5') as f:
f['/gridded_data_format/format_version'] = 1.0
f['/gridded_data_format/data_software'] = "Michael's Blog"
f['/gridded_data_format/data_software_version'] = 0
f['/gridded_data_format/data_comment'] = "Sample data generated for blog"
sp = f.create_group('/simulation_parameters').attrs
sp['refine_by'] = 0
sp['dimensionality'] = len(rho.shape)
sp['current_time'] = 0
sp['domain_dimensions'] = rho.shape
sp['domain_left_edge'] = (r.min(), z.min())
sp['domain_right_edge'] = (r.max(), z.max())
sp['unique_identifier'] = 'ring'
sp['cosmological_simulation'] = False
sp['num_ghost_zones'] = 0
sp['field_ordering'] = rho.flags.c_contiguous
sp['boundary_conditions'] = [1,0,0,0] # periodic
sp['geometry'] = 2 # cylindrical
f.create_group('/field_types/density')
f['/particle_types'] = "Sample data generated for blog"
f['/grid_left_index'] = [(0, 0)]
f['/grid_dimensions'] = [rho.shape]
f['/grid_level'] = [0] # Allows for grid refinement
f['/grid_parent_id'] = [-1] # ????
f['/grid_particle_count'] = [0, 0]
f['/data/grid_%010i/density' % (0,)] = rho
import yt
ds = yt.load('ring.hd5')
#ds = yt.load('RD0005.gdf')
#ds._parse_parameter_file()
ds.print_stats()
ds.parameters
yt.ProjectionPlot(ds, "z", "density")
Blender¶
Here we attempt to output the core structure in a voxel format following the suggestion here.
import struct
from vortex_ring import VortexRing
v = VortexRing(model='UFG', lam=2)
N = 64
x, y, z = v.grid3(N)
r = np.sqrt(y**2 + z**2)
psi0 = v.psi(x, r)
psi = v.psi(x, r, zr=[30,40])
rho0 = 2*abs(psi0)**2
rho = 2*abs(psi)**2
core = abs(rho - rho0)
xcent = 0.5
ycent = 0.5
zcent = 0.5
radius = 0.45
filename = "voxelsphere.raw"
data = ((core/core.max())*255).astype(int)
data = ((rho/rho.max())*255).astype(int)
print data.shape
with open(filename, "wb") as file:
for _d in reversed(data.ravel()):
file.write(struct.pack('B', _d))
MayaVi¶
MayaVi is a nice VTK-based visualization environment (similar to Visit and Paraview). There is a scripting mlab interface to this that allows you to generate plots. Once the plots are generated, you can then interact with them via the MayaVi GUI which will launch when you call mayavi.mlab.show.
One problem with this workflow is that interacting with MayaVi locks up the notebook: you need to close the MayaVi GUI window to continue. It would be nice to launch the MayaVi process in another python process and (maybe in a parallel IPython engine?) but I have not yet figured out how to do this yet.
A second issue is that you must run ipython with a Framework build of python so it can access the screen, otherwise you get the following message:
This program needs access to the screen.
Please run with a Framework build of python, and only when you are
logged in on the main display of your Mac.
My solution is presently to create a shell script called ipythonw that I include in my path (I put it in ~/usr/local/bin/ipythonw). Here is what I have in ipythonw:
#!/data/apps/anaconda/1.3.1/python.app/Contents/MacOS/python
if __name__ == '__main__':
from IPython import start_ipython
start_ipython()
I found out about this framework build of python by looking at the mayavi2 script installed by conda.
Single Image¶
We start by processing a single frame to show some of the things that can be done. Later we shall build an app to visualize all of the data.
%pylab wx
import mayavi
import mayavi.mlab as mlab
import vortex_ring;reload(vortex_ring);from vortex_ring import VortexRing
v = VortexRing(model='UFG', lam=2)
x, y, z, psi, (vx, vy, vz) = v.get_data(t=0.1, Nx=100)
r = np.sqrt(y**2 + z**2)
psi0 = v.psi(x, r)
rho0 = 2*abs(psi0)**2
rho = 2*abs(psi)**2
Here is a simple example of using MayaVi 2 to display the vortex ring. We plot the iso-density contours, and make them transparent so we can see into the object. The scale of the object is set by manipulating the extents.
extent = array([(_x.min(), _x.max()) for _x in x, y, z]).ravel()
mlab.clf()
contours = mlab.contour3d(rho, contours=20,
opacity=0.1,
#extent=extent, # This works here but is a pain elsewhere
)
mlab.flow(vx, vy, vz,
seedtype='plane', seed_resolution=20)
mlab.show(stop=True)
#x.shape
mlab.clf()
mlab.show(stop=True)
In order to make the vortex ring stand out, we consider the difference between the TF density and the vortex solution (this will need some work with a real simulation - perhaps using a high-pass filter?):
from mayavi import mlab
import mayavi.tools
psi = psis[0]
core = abs(psi0)**2 * (abs(abs(psi0)**2 - abs(psi)**2))
core /= core.max()
mlab.clf()
#mlab.contour3d(rho, contours=20, opacity=0.1)
pipeline = mayavi.mlab.pipeline.volume(mayavi.mlab.pipeline.scalar_field(rhos[58]), color=(1.0, 1.0, 1.0))
pipeline = mayavi.mlab.pipeline.volume(mayavi.mlab.pipeline.scalar_field(core), color=(0.0, 0.0, 1.0))
mayavi.mlab.show(stop=True)
Interactive Application¶
This is all better done in GUI where we can control the frame etc. The typical use-case will be to save all of the frames to an HDF5 file, and then load that up. Our application makes use of this, wrapping the data so that the entire data-set need not be loaded into memory at once. First we generate the data.
%%time
import os.path
import numpy as np
import h5py
import vortex_ring;reload(vortex_ring);from vortex_ring import VortexRing
v = VortexRing(model='UFG', lam=2)
Nx = 100
x, y, z, psi0, (vx, vy, vz) = v.get_data(t=0.0, Nx=Nx)
filename = "vortex_ring.hd5"
if os.path.exists(filename):
# One subtlety about hdf5 files is that you cannot
# insert an array if it is already there. Hence we
# ensure we delete the file first.
os.remove(filename)
with h5py.File(filename) as f:
# Create an extensible dataset that can grow.
f.create_dataset('psis', shape=(1,) + psi0.shape,
maxshape=(None,) + psi0.shape,
dtype=psi0.dtype)
f['ts'] = 0
ts = []
for t in np.linspace(0, 1, 100):
ts.append(t)
x, y, z, psi, (vx, vy, vz) = v.get_data(t=t, Nx=Nx)
with h5py.File(filename) as f:
psis = f['psis']
psis.resize(len(ts), axis=0)
psis[len(ts)-1] = psi
del f['ts']
f['ts'] = ts
if 'x' not in f.keys():
f['x'] = x
f['y'] = y
f['z'] = z
import numpy as np
import h5py
import traits.api as api
import traitsui.api as uapi
from traits.api import HasTraits, Int, Str, Array, Range, Instance, on_trait_change
from traitsui.api import View, Item, Group
from mayavi.core.api import PipelineBase
from mayavi.core.ui.api import MayaviScene, SceneEditor, MlabSceneModel
import pyface.timer.api
class MyModel(HasTraits):
rho0 = Array()
filename = Str("vortex_ring.hd5")
frames = Int(0)
frame = Range(0, 100, 0, exclude_high=True, mode='slider')
scene = Instance(MlabSceneModel, ())
rho_plot = Instance(PipelineBase)
core_plot = Instance(PipelineBase)
rho_contour = Range(0.0, 1.0, 0.5, mode='slider')
core_contour = Range(0.0, 1.0, 0.5, mode='slider')
play = api.Bool(False)
delay = api.Range(1, 1000, value=40)
timer = api.Instance(pyface.timer.api.Timer)
@property
def rho(self):
"""Return rho for the current frame"""
with h5py.File(self.filename, 'r') as f:
# Need to convert this to an array to actually read the data
psi = np.asarray(f['psis'][self.frame])
rho = abs(psi)**2
return rho
@property
def ts(self):
"""Return rho for the current frame"""
with h5py.File(self.filename, 'r') as f:
# Need to convert this to an array to actually read the data
ts = np.asarray(f['ts'])
return ts
def __init__(self, **traits):
self.frames = len(self.ts)
self.rho0 = 0*self.rho
for self.frame in xrange(self.frames):
self.rho0 += self.rho
self.rho0 /= self.frames
api.HasTraits.__init__(self, **traits)
self.start_timer()
def start_timer(self):
self.timer = pyface.timer.api.Timer(self.delay, self.onTimer)
def onTimer(self, *v):
r"""Callback responding to each timer tick and updating the frame."""
if self.play:
self.frame = ((self.frame + 1) % self.frames)
# When the scene is activated, or when the parameters are changed, we
# update the plot.
@on_trait_change('frame,rho_contour,core_contour,scene.activated')
def update_plot(self):
rho = self.rho
core = self.rho0 * abs(rho - self.rho0)
core /= core.max()
rho_contour = self.rho_contour * rho.max()
core_contour = self.core_contour * core.max()
if self.rho_plot is None or self.core_plot is None:
#self.rho_plot = self.scene.mlab.pipeline.volume(
# self.scene.mlab.pipeline.scalar_field(rho),
# color=(1.0, 1.0, 1.0))
#self.core_plot = self.scene.mlab.pipeline.volume(
# self.scene.mlab.pipeline.scalar_field(core),
# color=(0.0, 0.0, 1.0))
self.rho_plot = self.scene.mlab.pipeline.iso_surface(
self.scene.mlab.pipeline.scalar_field(rho),
contours=[rho_contour],
opacity=0.1,
color=(1.0, 1.0, 1.0))
self.core_plot = self.scene.mlab.pipeline.iso_surface(
self.scene.mlab.pipeline.scalar_field(core),
contours=[core_contour],
color=(0.0, 0.0, 1.0))
else:
self.rho_plot.mlab_source.set(scalars=rho)
self.rho_plot.contour.contours = [rho_contour]
self.core_plot.mlab_source.set(scalars=core)
self.core_plot.contour.contours = [core_contour]
def _delay_changed(self, name, old, new):
if old != new:
self.timer.Stop()
del self.timer
self.start_timer()
# The layout of the dialog created
view = uapi.View(
uapi.Item('scene', editor=SceneEditor(scene_class=MayaviScene),
height=250, width=300, show_label=False),
uapi.Group(
'_',
'frame',
'rho_contour',
'core_contour',
uapi.HGroup(
uapi.Item(name='play'),
uapi.Item(name='delay',
editor=uapi.DefaultOverride(mode='slider')))),
resizable=True,
)
my_model = MyModel()
my_model.configure_traits()
MayaVi in an Engine¶
If you have started a cluster, you can use the following code to run MayaVi on one of those engines. Here we define a %pxpush magic that allows us to push data to the engine.
from IPython import get_ipython
from IPython.core.magic import Magics, line_magic, magics_class, register_line_magic
@register_line_magic
def pxpush(line):
vars = line.replace(',', ' ').split()
ip = get_ipython()
view = ip.magics_manager.registry['ParallelMagics'].view
view.push({_v: ip.ev(_v) for _v in vars})
Now try to connect to the engine (you must start this separately from the Dashboard) and then set the default target to the zeroth process.
import IPython.parallel
c = IPython.parallel.Client()
%pxconfig -t 0
Compute the data, then push it to the client.
x, y, z = v.grid3(100)
r = np.sqrt(y**2 + z**2)
psi = v.psi(x, r, xy=[30,40])
%pxpush x, y, z, psi
Here we do the plotting. By specifying stop=True, a dialog will pop up allowing you to suspend the interaction. This will keep the plot window open (and in the position you like).
%%px
%pylab wx
import mayavi
import mayavi.mlab
extent = array([(_x.min(), _x.max()) for _x in x, y, z]).ravel()
mayavi.mlab.contour3d(2*abs(psi)**2, contours=20,
opacity=0.1,
extent=extent)
mayavi.mlab.show(stop=True)
psi = v.psi(x, r, xy=[30,100])
%pxpush psi
%%px
mayavi.mlab.show(stop=True)
mayavi.mlab.contour3d(2*abs(psi)**2, contours=20,
opacity=0.1,
extent=extent)
mayavi.mlab.show(stop=True)
POV-Ray¶
POV-Ray is a command-line ray-tracing program. It has a simple markup language for describing and rendering scenes. Here is an example from the tutorial.
from IPython.core.magic import register_cell_magic
from IPython.core.display import Image
@register_cell_magic
def povray(name, cell):
ini_name = None
if name == '':
raise ValueError("Must provide a filename")
elif name.endswith('.ini'):
ini_name = name
name = name[:-4]
if not os.path.exists(ini_name):
raise valueError("Ini file %s must exist" % (ini_name, ))
pov_name = "%s.pov" % (name,)
image_name = "%s.png" % (name,)
if os.path.exists(image_name):
os.remove(image_name)
with open(pov_name, 'w') as _f:
_f.write(cell)
if ini_name:
res = !povray $ini_name
else:
res = !povray +L/data/apps/povray/3.7/POV-Ray3_7_Mac_Unofficial/include/ $pov_name
if os.path.exists(image_name):
return Image(image_name)
else:
print(res.n)
%%povray demo
#include "colors.inc"
background { color Cyan }
camera {
location <0, 2, -3>
look_at <0, 1, 2>
}
sphere {
<0, 1, 2>, 2
texture {
pigment { color Yellow }
}
}
light_source { <2, 4, -3> color White}
As with Blender, we can store the data as a voxel format. The header is 6 bytes with the three-dimensional volume, then the voxel cells follow
import numpy as np
import struct
import vortex_ring;reload(vortex_ring)
from vortex_ring import VortexRing
v = VortexRing(model='UFG', lam=2)
N = 256
x, y, z = v.grid3(N)
r = np.sqrt(y**2 + z**2)
psi0 = v.psi(x, r)
psi = v.psi(x, r, zr=[30,40])
rho0 = 2*abs(psi0)**2
rho = 2*abs(psi)**2
core = abs(rho - rho0)
def write(density, name):
data = ((density/density.max())*255).astype(np.uint8)
filename = "%s.df3" % (name,)
print data.shape
with open(filename, "wb") as file:
file.write(np.array(data.shape).astype(np.uint16).byteswap(True))
file.write(data.T.ravel().astype(np.uint8))
write(rho, 'rho')
write(core, 'core')
%%povray demo
#include "colors.inc"
background { color Black }
camera {
location <0, 0, 3>
look_at <0, 0, 0>
}
#declare boxtexture = texture {
pigment {
rgbf 1
}
}
#declare boxinterior = interior {
media {
intervals 100
samples 1,20
emission <1,0,0>
absorption <0,0,0>
scattering { 1,<0,0,0> }
confidence 0.9999
confidence 0.5
variance 1/1000
variance 1/16
density {
density_file df3 "demo.df3" interpolate 1
}
}
}
box {
<0, 0, 0>, <1, 1, 1>
texture { boxtexture }
interior { boxinterior }
translate <0, 0, 0>
scale <1, 1, 2>
hollow
}
light_source { <2, 4, -3> color White}
Nx, Ny, Nz = data.shape
print Nx, Ny, Nz
%%file demo.ini
Input_File_Name=demo.pov ; input file name(s)
Library_Path=/data/apps/povray/3.7/POV-Ray3_7_Mac_Unofficial/include/
+W400 +H400 ; resolution: W(idth) times H(ight)
Display=on
Quality=9 ; 0 .. 11 (best), default: 9
Bounding=on
;Bounding_Threshold=25
Light_Buffer=on ; (requires bounding on)
Vista_Buffer=on ; (requires bounding on)
Antialias=on
Sampling_Method=1 ; 1 or 2 (both adaptive, 2 recursive)
Antialias_Threshold=0.1 ;
Antialias_Depth=6 ; 1 .. 9 (best), default: 6
Jitter=on ;
Jitter_Amount=1.0 ; 0.0 .. 1.0 (best)
Field_Render=off
Output_to_File=on
Output_File_Type=N ; N for PNG
Output_Alpha=off
Bits_per_Color=16 ; 8 .. 16 for PNG
;Output_File_Name=
%%povray demo.ini
#include "colors.inc"
#declare NX = 256;
#declare NY = 128;
#declare NZ = 128;
#declare DIAG = <NX,NY,NZ>;
global_settings
{ ambient_light <1,1,1>
assumed_gamma 1
}
camera
{ location <0,-7/4*NY,2/3*NZ>
up z
right x // default: 4/3*x
sky <0,0,1>
look_at <0,0,0>
}
light_source
{ <2*NX,-NY,2*NZ>
color rgb <1,1,1>
media_interaction on
media_attenuation on
shadowless
}
#declare DENS = interior
{ media
{ intervals 100 // number of ray-intervals
ratio 0.5
samples 3,3 // maximum,minimum of samples per voxel
method 3 // 1, 2 or 3 (3 requires samples 3,3)
emission 3*<1,1,1>/1000
absorption <0,0,0>/1000
scattering { 1, <0,0,0> }
confidence 0.999 // default: 0.9
variance 1/10000 // default: 1/128
density
{ density_file df3 "rho.df3"
interpolate 1
color_map // colour map with (smooth) linear transition(s)
{ [0.0 rgb <0.0,0.0,0.0>] // 0.0 ~ 'black'
[1.0 rgb <1.0,1.0,1.0>] // 0.2 ~ 'white'
}
}
}
media
{ intervals 100 // number of ray-intervals
ratio 0.5
samples 3,3 // maximum,minimum of samples per voxel
method 3 // 1, 2 or 3 (3 requires samples 3,3)
emission 3*<1,1,1>/10
absorption <0,0,0>/1000
scattering { 1, <0,0,0> }
confidence 0.999 // default: 0.9
variance 1/10000 // default: 1/128
density
{ density_file df3 "core.df3"
interpolate 1
color_map // colour map with (smooth) linear transition(s)
{ [0.0 rgb <0.0,0.0,0.0>] // 0.0 ~ 'black'
[1.0 rgb <0.0,0.0,1.0>] // 0.2 ~ 'blue'
}
}
}
}
box
{ <0,0,0>, <1,1,1>
pigment { rgbt <0,0,0,1> }
hollow
interior { DENS }
scale DIAG
translate -DIAG/2
rotate <0,0,45>
//rotate <0,0,360*clock> // rotation around z-axis
}
