>>>>> "Massimo" == DiPierro, Massimo <[EMAIL PROTECTED]> writes:
Massimo> www.fermiqcd.net <http://www.fermiqcd.net/> and I am
Massimo> starting a new project for visualization in lattice QCD.
I've attached some code which you may find useful as a starter. Note
that it will NOT be good for you right away because it is quite old,
written for Mayavi1, plain VTK and Numeric.
But still, some of what's there may provide a first cut. There are
classes in there to load Topological Densities in the format produced
by the version 5 MILC code in pure gauge configurations (from the
binary lattice files). My use for it was to slice the field
configurations (I also did some fermion density plots using the same
code) in any direction and generate HTML tables of isosurfaces for
easy inspection of a large number of lattices. I found it easier to
look at the data this way by dumping 'contact sheets' of the field
config than to manually load hundreds of configs.
If this of any use to you, do with it as you wish; I don't use this
code at all anymore.
Cheers,
f
"""Miscellaneous utilities for VTK programming."""
__author__ = 'Fernando Perez'
__email__ = '[EMAIL PROTECTED]'
#--------------------------------------------------------------------------
# Needed modules
import types,copy
from mayavi import ivtk
from vtkpython import *
#--------------------------------------------------------------------------
# Functions
def __make_globals():
global opacity_fns,color_fns
class table:pass # dummy class
# Opacity functions
op = table()
op.flat = [(0,1),(1,1)] # minimal valid list without symmetrization
# various v-shaped functions which hide mid-range values
op.simple_v = {'points':[(0,1),(0.5,0)],'symmetric':1}
op.v1 = {'points':[(0,1),(.1,.9),(.3,.2),(.5,0)],'symmetric':1}
op.v2 = {'points':[(0,1),(.05,.8),(.15,.5),(.3,0.05),(.4,0.01),(.5,0)],
'symmetric':1}
op.v3 = {'points':[(0,1),(.05,.8),(.15,.5),(.3,0.01),(.4,0.0),(.5,0)],
'symmetric':1}
op.v4 = {'points':[(0.0, 0.94948854441504127),
(0.035714285714285712, 0.80424645547447116),
(0.071428571428571425, 0.65109382296701002),
(0.10714285714285714, 0.49696568131514102),
(0.14285714285714285, 0.35009836810263129),
(0.17857142857142855, 0.2193496591388204),
(0.21428571428571427, 0.11333665732781846),
(0.25, 0.039504140754667641),
(0.2857142857142857, 0.0032619764428993394),
(0.3214285714285714, 0.0),
(0.5, 0.0)],
'symmetric':1}
op.v5 = {'points':[(0.0, 0.99762935250668083),
(0.035714285714285712, 0.9184236784755615),
(0.071428571428571425, 0.82533229455592594),
(0.10714285714285714, 0.71980451982071758),
(0.14285714285714285, 0.60472644671805775),
(0.17857142857142855, 0.48445139475713539),
(0.21428571428571427, 0.36463421932744544),
(0.25, 0.25185472655077878),
(0.2857142857142857, 0.15305251748275239),
(0.3214285714285714, 0.074834654304759413),
(0.3571428571428571, 0.02274961903772077),
(0.39285714285714285, 0.00063759243707800506),
(0.42857142857142855, 0.0),
(0.5, 0.0)],
'symmetric':1}
op.v_gauss = {'points':[(0.0, 0.99326205300091452),
(0.035714285714285712, 0.98658288958329921),
(0.071428571428571425, 0.97461176465766319),
(0.10714285714285714, 0.95434922494481678),
(0.14285714285714285, 0.92199796853146232),
(0.17857142857142855, 0.87334982624593316),
(0.21428571428571427, 0.80458955009027378),
(0.25, 0.71349520313980985),
(0.2857142857142857, 0.60082978415651644),
(0.3214285714285714, 0.47152252861281918),
(0.3571428571428571, 0.33512968029560419),
(0.39285714285714285, 0.2051420318745103),
(0.42857142857142855, 0.097007305927970022),
(0.46428571428571425, 0.025187568150640605),
(0.5, 0.0)],
'symmetric':1}
[ 0.99326205, 0.98173503, 0.95568107, 0.90374171, 0.81286037, 0.67433423,
0.49270894, 0.29267172, 0.11719747, 0.01375494, 0.01375494,
0.11719747, 0.29267172, 0.49270894, 0.67433423, 0.81286037,
0.90374171, 0.95568107, 0.98173503, 0.99326205]
op.sq = {'points':[(0,1),(.25,1),(.26,0),(.5,0)],'symmetric':1}
# Color functions
cf = table()
cf.red_turq = {'points':[(0,1,0,0),(.5,1,1,0),(1,0,1,1)],'mode':'RGB'}
cf.red_blue = {'points':[(0,0,1,1),(1,0.667,1,1)],'mode':'HSV'}
cf.blue_red = {'points':[(0,0.667,1,1),(1,0,1,1)],'mode':'HSV'}
# set globals
opacity_fns = op
color_fns = cf
__make_globals()
def symmetrize(lst):
"""Symmetrize a list of (x,y) pairs about the endpoint.
Returns a list with 2n-1 elements (n=len(lst)) which is now
symmetric about its midpoint.
The original list is unchanged.
Note. The algorithm is simple and inefficient. This is meant to be
used for simple cases such as transfer functions in volume rendering,
so zero effort went into optimizing it. Don't use this to manipulate
huge lists."""
mirror = copy.deepcopy(lst)
mirror.reverse()
sym = copy.deepcopy(lst)
x_end = 2*lst[-1][0]
midpt = len(sym)
sym.extend(mirror[1:])
for n in range(midpt,len(sym)):
sym[n] = (x_end - sym[n][0],sym[n][1])
return sym
def build_PiecewiseFunction(pwf,nmin,nmax,pts):
"""Build a vtkPiecewiseFunction given a list of points.
build_PiecewiseFunction(pwf,nmin,nmax,pts):
-pwf: vtkPiecewiseFunction object, modified in-place.
-nmin,nmax: range of input values for the function (0,255 typically)
-pts: a list or tuple of (x,y) pairs mapped to the [0,1] interval. The
first pair must be (0,y0) and the last must be (1,y_fin) to ensure the
function is defined over the complete (nmin,nmax) range unless the data is
to be symmetrized.
In order to symmetrize a list, pts must instead be a dict with the form:
pts = {'points':[(0,y0),...,(0.5,y_mid)],'symmetric':1}
-symmetric: if 1, the last pair must be (0.5,y_fin) and the list willl
be symmetrized about its midpoint."""
if type(pts) == types.DictType:
if pts['symmetric']:
pts = symmetrize(pts['points'])
else:
raise ValueError,'dicts can only be used to symmetrize lists'
if pts[0][0] != 0 or pts[-1][0] != 1:
raise ValueError,'first,last pair must be (0,y0) and (1,y_fin)'
for pt in pts:
pwf.AddPoint(nmin+int(pt[0]*nmax),pt[1])
def build_ColorTransferFunction(ctf,nmin,nmax,pts):
"""Build a vtkColorTransferFunction given a list of points.
Similar to build_PiecewiseFunction(), except it doesn't have the
symmetrization options.
pts is a dict of the form:
{'points':[(0,v1_1,v1_2,v1_3),...,(1,vN_1,vN_2,vN_3)],
'mode':'RGB'}
Use mode 'RGB' (default) or 'HSV' to specify how the above vi_j are to be
interpreted."""
if pts['points'][0][0] != 0 or pts['points'][-1][0] != 1:
raise ValueError,'first,last pair MUST be (0,y0) and (1,y_fin)'
if pts['mode'] == 'RGB':
add_segment = vtkColorTransferFunction.AddRGBSegment
elif pts['mode'] == 'HSV':
add_segment = vtkColorTransferFunction.AddHSVSegment
else:
raise ValueError,"mode can only be 'RGB' or 'HSV'"
# We must do it with segments to be able to build tables which go
# 'backwards' in colorspace (doesn't work with AddPoint)
pt = pts['points']
for n in range(len(pt)-1):
ptn = pt[n]
add_segment(ctf,nmin+int(ptn[0]*nmax),ptn[1],ptn[2],ptn[3],
nmin+int(pt[n+1][0]*nmax),*pt[n+1][1:4])
def volume_render(vtk_file,
mapper='texture',
dsample = 1,
opacity_function = opacity_fns.v3,
color_function = color_fns.blue_red,
background_color = (0.93,0.93,0.84),
outline_color = (0.1,0.1,0.1) ):
"""Open a VTK volume renderer for the given file.
Details for some parameters (the others are self-explanatory):
- mapper: 'ray' or 'texture' (default). Ray tracing is MUCH slower.
- dsample: distance between samples for ray tracing
interpolation. Extremely slow if < 1, but gives very nice images.
"""
# Build the vtk pipeline
# Data reader
reader = vtkStructuredPointsReader()
reader.SetFileName(vtk_file)
# Texture mappers are fairly fast
textureMapper = vtkVolumeTextureMapper2D()
textureMapper.SetInput(reader.GetOutput())
# We can also use a much slower ray tracing mapper
rayMapper = vtkVolumeRayCastMapper()
rayMapper.SetInput(reader.GetOutput())
compositeFunction = vtkVolumeRayCastCompositeFunction()
rayMapper.SetVolumeRayCastFunction(compositeFunction)
rayMapper.SetSampleDistance(dsample) # VERY SLOW if dsample < 1
# Transfer functions for opacity and color
nmin = 0; nmax=255 # volumetric data is in bytes
opacity = vtkPiecewiseFunction()
build_PiecewiseFunction(opacity,nmin,nmax,opacity_function)
color = vtkColorTransferFunction()
build_ColorTransferFunction(color,nmin,nmax,color_function)
# Fix in the volume properties (opacity+color)
volumeProperty = vtkVolumeProperty()
volumeProperty.SetColor(color)
volumeProperty.SetScalarOpacity(opacity)
volumeProperty.ShadeOn()
volumeProperty.SetInterpolationTypeToLinear()
# Create the volume actor
dataVolume = vtkVolume()
dataVolume.SetProperty(volumeProperty)
# Choose the mapper to use. Default to texture (faster)
if mapper=='ray':
dataVolume.SetMapper(rayMapper)
else:
dataVolume.SetMapper(textureMapper)
# Put an outline around it
outline = vtkOutlineFilter ()
outline.SetInput (reader.GetOutput ())
outline_map = vtkPolyDataMapper ()
outline_map.SetInput(outline.GetOutput())
outline_act = vtkActor ()
outline_act.SetMapper(outline_map)
outline_act.GetProperty().SetColor(outline_color)
# Let's use Prabhu's great ivtk to see the thing
iv = ivtk.viewer()
camera = iv.renwin.GetActiveCamera()
camera.SetPosition (-25,-50,30)
camera.SetViewUp (0,.3,.95)
iv.renwin.set_background(background_color)
iv.AddActors([dataVolume,outline_act])
#******************** End of file <vtkutils.py> ****************************
"""Classes for visualization of 4d scalar data using MayaVi.
"""
__author__ = "Fernando Perez <[EMAIL PROTECTED]>"
#*****************************************************************************
# Needed modules
import sys,struct,types,time,shutil,os
try:
import HTMLgen as html
except ImportError:
print "Warning: no HTML generation available"
import mayavi,gracePlot,pyvtk
from Numeric import *
from vtkpython import *
from Scientific.Statistics.Histogram import Histogram
# my own generic stuff
from IPython.genutils import *
from IPython.numutils import *
import vtkutils
# stuff from modules
ivtk = mayavi.ivtk
#*****************************************************************************
class Data4d:
"""Class for visualizing 4-dimensional scalar data.
It is meant to be subclassed by classes which implement at least the
constructor with appropriate routines for reading the data files.
"""
CameraPosition = (-25,-50,30)
CameraViewUp = (0,.3,.95)
def __init__(self):
raise NotImplementedError
def cut_slice(self,dir,N):
"""Slice a data set along direction dir at N.
Note that vtk datasets use the same (odd) convention as topological
density files: data is indexed 'in reverse':
data(x,y,z,t) = data_arr[t,z,y,x].
For this reason the slicing below looks 'backward'."""
if dir == 'x':
self.slice = self.data[:,:,:,N]
elif dir == 'y':
self.slice = self.data[:,:,N,:]
elif dir == 'z':
self.slice = self.data[:,N,:,:]
elif dir == 't':
self.slice = self.data[N,:,:,:]
else:
raise ValueError,'direction must be one of x,y,z,t.'
self.slice_dir = dir
self.slice_val = N
self.axes_labels = qw('x y z t')
self.axes_labels.remove(dir)
# compute the histogram for the 3d slice
self.histo3d = Histogram(ravel(self.slice),100)
def write_text_file(self,fname):
"""Dump the 4d data to a simple text file."""
outfile = file(fname,'w')
nx,ny,nz,nt = self.nx,self.ny,self.nz,self.nt
print >> outfile, 'nx',nx
print >> outfile, 'ny',ny
print >> outfile, 'nz',nz
print >> outfile, 'nt',nt
for it in range (nt):
for iz in range (nz):
for iy in range (ny):
for ix in range (nx):
print >> outfile, ix,iy,iz,it,self.data[it,iz,iy,ix]
def reset(self):
"""Reset the histogram clipping"""
self.cut_slice(self.slice_dir,self.slice_val)
def histo3_plot(self,plotter=None):
"""Print min/max data and plot a histogram of the current 3-slice."""
histo = self.histo3d
if self.clipped:
type = 'Clipped'
else:
type = 'Unclipped'
if plotter == None:
plotter = gracePlot.gracePlot()
plotter.histoPlot(histo)
plotter.title('%s=%s %s 3-d slice' %
(self.slice_dir,self.slice_val,type) )
print 'Min:',histo.min,'\nMax:',histo.max
def histo4_make(self,force=0):
"""Make 4-d histogram."""
if force:
try:
del self.histo4
except AttributeError: pass
try:
self.histo4
except AttributeError:
nbins = 250
print 'Building histograms of full dataset, this may take a moment...',
sys.stdout.flush()
self.histo4 = Histogram(ravel(self.data),nbins)
print
def histo4_plot(self,plotter=None):
"""Plot the 4-d histograms.
If they haven't been made, it builds them automatically."""
self.histo4_make()
plotter = gracePlot.gracePlot()
plotter.title('4-d data set')
plotter.histoPlot(self.histo4)
def clip_slice(self,min_c=None,max_c=None):
self.min_c = min_c or self.histo3d.min
self.max_c = max_c or self.histo3d.max
self.slice = clip(self.slice,self.min_c,self.max_c)
self.histo3d = Histogram(ravel(self.slice),100)
self.clipped = 1
def clip_all(self,min_c=None,max_c=None,fraction=None):
"""Clip all data"""
self.histo4_make()
if fraction is None:
min_c = min_c or self.histo4.min
max_c = max_c or self.histo4.max
self.data = clip(self.data,min_c,max_c)
else:
# this should find the cuts in the histogram where the data is a
# certain fraction of the total and then do the actual clipping
raise NotImplementedError
hist_bins = self.histo4.array[:,0]
hist_counts = self.histo4.array[:,1]
count_max = max(self.histo4.array[:,1])
count_min = fraction * count_max
max_c = self.histo4.max
min_c = None
self.histo4_make(force=1)
def make_vtk(self):
"""Build a vtk data obj from the clipped 3d slice."""
dims = list(self.slice.shape)
# We need to give the dims as (x,y,z), not (z,y,x):
dims.reverse()
grid = pyvtk.StructuredPoints(dims)
# RectilinearGrid is another option for this [use (x,y,z) order]:
#grid = pyvtk.RectilinearGrid(range(dims[2]),range(dims[1]),range(dims[0]))
if self.slice.iscontiguous():
dat = self.slice.flat
else:
dat = ravel(self.slice)
point_data = pyvtk.PointData(pyvtk.Scalars(dat,
'Data Density: '+self.slice_dir+
'='+`self.slice_val`,'default'))
header = 'Data density slice %s=%s' % (self.slice_dir,self.slice_val)
return pyvtk.VtkData(grid,header,point_data)
def make_vtk_vol(self):
"""Build a vtk data obj from the clipped 3d slice."""
# Define function which can scale density to range nmin, nmax
densmax=self.histo3d.min
densmin=self.histo3d.max
nmin=0
nmax=255
dims = list(self.slice.shape)
# We need to give the dims as (x,y,z), not (z,y,x):
dims.reverse()
grid = pyvtk.StructuredPoints(dims)
# RectilinearGrid is another option for this [use (x,y,z) order]:
#grid = pyvtk.RectilinearGrid(range(dims[2]),range(dims[1]),range(dims[0]))
if self.slice.iscontiguous():
dat = self.slice.flat
else:
dat = ravel(self.slice)
# Data must be stored as a byte array for VTK's volume rendering to
# work. Is the map(int()) really needed?
dat = array(map(int,(dat-densmin)/(densmax-densmin)*(nmax-nmin)+nmin),
'b')
point_data = pyvtk.PointData(pyvtk.Scalars(dat,
'Data Density: '+self.slice_dir+
'='+`self.slice_val`,'default'))
header = 'Data density slice %s=%s' % (self.slice_dir,self.slice_val)
return pyvtk.VtkData(grid,header,point_data)
def frac2val(self,frac):
"""Return the actual value of a fraction (0-1) of the data's range.
"""
return self.histo3d.min + frac*(self.histo3d.max - self.histo3d.min)
def load_isosurface(self,isoval,mode='relative'):
"""Load an isosurface at val (as a fraction of data range)"""
isosurface = self.mv.load_module('IsoSurface',0)
if mode=='relative':
isoval = self.frac2val(isoval) # turn the fraction into a data value
else:
if isoval < self.histo3d.min:
isoval = self.frac2val(0.01)
elif isoval > self.histo3d.max:
isoval = self.frac2val(0.99)
#print 'ISO val=%s, mode=%s' % (isoval,mode) # dbg
isosurface.slider_pos = isoval
isosurface.normals_on.set(1)
isosurface.do_normals()
isosurface.cont_fil.SetValue (0, isoval)
# call this explicitly (it was already called at construction time) to
# update the gui with our changes:
isosurface._gui_init()
def getLutIndex(self,l, val):
c = l.GetColor(val)
for i in range(l.GetNumberOfColors()):
if l.GetTableValue(i)[:-1] == c:
return i
def square_window(self,lut,min_val,max_val,invert=0):
"""In-place modify a lut with a square window.
It sets alpha outside of min/max values (given as data values) to
0 and 1 inside. If invert=1, the alpha table is the opposite.
"""
imin = self.getLutIndex(lut,min_val)
imax = self.getLutIndex(lut,max_val)
# Set up table of transparencies
lut_nvalues = lut.GetNumberOfTableValues()
alpha = invert * ones(lut_nvalues,'f')
alpha[imin:imax] = 1-invert
# Change the lut according to the alpha table
for i in range(0,lut_nvalues):
val = lut.GetTableValue (i)
lut.SetTableValue(i, val[0], val[1], val[2], alpha[i])
def load_ContourGridPlane(self,axis=0,pos=0,contours=0):
"""Load a ContourGridPlane module at given axis/position.
- axis: use the convention: x,y,z == 0,1,2
- pos can be given as -1 for the end of an axis.
- contours: Total # of contours per plane. If contours is negative,
then colormaps are shown instead of individual contour lines. This is
much faster but gives colormaps which 'bleed' out of the isosurfaces,
since they are 1x1 squares and not smoothed to follow the isosurface
contours.
"""
cplane = self.mv.load_module('ContourGridPlane',0)
cplane.axis_var.set(axis)
if pos == -1:
pos = cplane.dim[axis]
cplane.slider_pos = pos
if contours > 0:
cplane.contour_on.set (1)
cplane.n_cnt.set (contours)
cplane.do_contour()
cplane.config_extents() # needed so the plane actually shows
def box_countour_planes(self,contours,invert=0,mode='relative'):
"""Put contour planes around the 3d box with a clipped colormap.
They are loaded in a separate module manager.
Parameters:
- contours: Total # of contours per plane. If contours is negative,
then colormaps are shown instead of individual contour lines. This is
much faster but gives colormaps which 'bleed' out of the isosurfaces,
since they are 1x1 squares and not smoothed to follow the isosurface
contours.
- invert: whether to invert the lookup table.
"""
self.add_cutoff_lut_mod_manager(self.iso_val0,self.iso_val1,
invert,mode)
if contours>0:
print 'Building contour box. This will take a while...'
for dir in range(3):
self.load_ContourGridPlane(dir,0,contours)
self.load_ContourGridPlane(dir,-1,contours)
def add_cutoff_lut_mod_manager(self,min_val,max_val,invert=0,
mode='relative'):
"""Add a module manager with a modified lut for transparency.
"""
#print 'cutoff lut. min=%s, max=%s, invert=%s' % \
# (min_val,max_val,invert) # dbg
# new module manager.
self.mv.load_mm('new')
dvm = self.mv.get_current_dvm()
mm = dvm.get_current_module_mgr()
lh = mm.scalar_lut_h # get the scalar_lut_handler (LutHandler)
lut = lh.lut # get the vtkLookupTable
if mode=='relative':
min_val = self.frac2val(min_val)
max_val = self.frac2val(max_val)
#print 'min %s max %s ' % (min_val, max_val) # dbg
self.square_window(lut,min_val,max_val,invert)
def view(self,iso_val0=0.25,iso_val1=None,check_histo=1,contours=-1,
invert=1,iso_mode='relative',legend=1):
"""Visualize in MayaVi the currently active slice.
Parameters:
- iso_val0: lower value of the isosurface cut plane.
- iso_val1: higher value of the isosurface cut plane.
- check_histo: show histogram to check for data clipping or not.
- contours: how many contours to show (0 disables).
- invert: whether the contour planes should use an inverted LUT.
- iso_mode: 'relative'/'absolute'. Whether isosurfaces are specified
as [0,1] relative values or as absolute data values. """
# relative_3 and relative are the same for this function (which
# doesn't know about the 4d data)
if iso_mode.startswith('relative'):
iso_mode = 'relative'
self.iso_val0 = iso_val0
if iso_val1 is None:
if iso_mode == 'relative':
iso_val1 = 1 - iso_val0
else:
raise ValueError,'in absolute mode you must give max value'
self.iso_val1 = iso_val1
if check_histo:
plotter = gracePlot.gracePlot()
while check_histo:
self.histo3_plot(plotter)
do_clip = raw_input('\nClip values [N]? ')
if do_clip.lower() != 'y':
break
low = float(raw_input('Low ['+`self.histo3d.min`+'] = ') or
self.histo3d.min)
high = float(raw_input('High ['+`self.histo3d.max`+'] = ') or
self.histo3d.max)
self.clip_slice(low,high)
fname3d = '%s.%s%s.vtk' % (self.fname,self.slice_dir,self.slice_val)
self.fname3d = fname3d
self.vtkdata = self.make_vtk()
self.vtkdata.tofile(fname3d)
# The config option turns on/off showing a GUI control
self.mv = mayavi.mayavi()
self.mv.open_vtk(fname3d,config=0)
# Let's not clutter up with all these temp files
os.unlink(fname3d)
# Now load modules and configure some by hand:
# Mayavi bug: load the isosurfaces *first*, otherwise window repainting
# breaks
self.load_isosurface(iso_val0,iso_mode)
self.load_isosurface(iso_val1,iso_mode)
self.mv.load_module('Outline',0)
axes = self.mv.load_module('Axes',0)
for axis,label,i in zip(qw('X Y Z'),self.axes_labels,range(3)):
apply (eval ("axes.axes.Set%sLabel"%(axis)), (label,))
axes.axis_txt_var[i].set(label)
loc = self.mv.load_module('Locator',0)
loc.axes.SetOrigin([0,0,0])
loc.coord_var.set([0,0,0])
# Rotate view to a nice position and angle
view = Data4d.CameraPosition + Data4d.CameraViewUp # join tuples
self.mv.renwin.update_view(*view)
# add a legend
if legend:
dvm = self.mv.get_current_dvm()
mm = dvm.get_current_module_mgr ()
slh = mm.get_scalar_lut_handler ()
# turn legend on
slh.legend_on.set(1) # set to 1/0 for on/off
slh.legend_on_off () # this _must_ be called to apply change
# Contour planes. Expensive if contours > 0, cheap if <0.
if contours:
self.box_countour_planes(contours,invert,iso_mode)
def save_png(self,dir='',fname=None,size=None):
"""Save snapshot as a png file. Give name without extension.
If no filename is given, the name of the 3-d slice is used.
If size is given, convert is called with 'convert -geometry size
fname' (so it does an in-place rewrite of the file to the given size
spec. """
if fname is None:
fname = self.fname3d.replace('vtk','png')
elif fname.endswith('.png'):
pass
else:
fname = fname +'.png'
fname = os.path.join(dir,fname)
#print '***png file',fname # dbg
self.mv.renwin.save_png(fname)
if size is not None:
xsys('convert -geometry %s %s %s' % (size,fname,fname))
def save_eps(self,dir='',fname=None,zip=0):
"""Save snapshot as a PostScript file. Give name without extension.
By default the file is gzipped, use zip=0 to prevent this."""
if fname is None:
fname = self.fname3d.replace('vtk','eps')
elif fname.endswith('.eps'):
pass
else:
fname = fname +'.eps'
fname = os.path.join(dir,fname)
print '***ps file',fname # dbg
self.mv.renwin.save_ps(fname)
# fix this file and zip it
tmp_file = os.path.join(dir,'__tmp.eps')
xsys('eps2eps %s %s' % (fname,tmp_file),verbose=1)
os.rename(tmp_file,fname)
if zip:
xsys('gzip -f %s' % fname)
def volume(self,mapper='texture',dsample=1,opacity=None):
"""Show volumetric density.
Call with mapper='ray' for using ray tracer instead of texture mapper.
dsample: distance between samples for ray tracer. If <1 it gets slow.
opacity: list of (x,y) pairs for opacity function (see vtkutils.py).
If unspecified, it defaults to the instance's opacity table
(self.opacity). Absent that, to a builtin default.
"""
# Defaults
if opacity is not None:
opacity_function = opacity
elif hasattr(self,'opacity'):
opacity_function = self.opacity
else:
opacity_function = vtkutils.opacity_fns.v3
fname3d_vol = '%s.%s%s.vol.vtk' % (self.fname,self.slice_dir,self.slice_val)
self.vtkdata_vol = self.make_vtk_vol()
self.vtkdata_vol.tofile(fname3d_vol)
vtkutils.volume_render(fname3d_vol,mapper=mapper,
opacity_function=opacity_function,
color_function=vtkutils.color_fns.blue_red,
background_color=(1,1,1))
def slice_table(self,dir='x',step=1,fname=None,
iso_val0=0.25,iso_val1=None,iso_mode='relative',
contours=-1,invert=1,
out_dir='',image_size=310,make_eps=0):
"""Make an HTML page with a table of 3-d slices.
Use iso_mode = 'relative_3' for each 3d slice to have its own
isosurfaces computed with values relative to the its histogram only
(and not to the full 4d histogram).
"""
if fname is None:
fname = os.path.basename(self.fname)
html_dir = os.path.join(out_dir,fname+'.slices')
try:
os.mkdir(html_dir)
except:
pass
# for the full table, reset the iso_values so that they are relative
# to the data range ends, but measured on the full 4d dataset. Once
# the values are determined, then the mode should become absolute
if iso_mode=='relative':
self.histo4_make()
h4 = self.histo4
if iso_val1 is None:
iso_val1 = 1 - iso_val0
iso_val0 = h4.min + iso_val0*(h4.max-h4.min)
iso_val1 = h4.min + iso_val1*(h4.max-h4.min)
iso_mode = 'absolute'
# Build the slice table, which can fail for various reasons (bad or
# zero data, typically).
images = []
dim = dict(zip(qw('t z y x'),self.data.shape))
start = time.time()
for slice in range(0,dim[dir],step):
try:
self.cut_slice(dir,slice)
except:
warn('Slice %s could not be generated. \n'
'Removing HTML dir %s and exiting slice_table().'
% (slice, html_dir))
shutil.rmtree(html_dir)
return
else:
self.view(iso_val0=iso_val0,iso_val1=iso_val1,iso_mode=iso_mode,
check_histo=0,contours=contours,invert=invert)
pngfile = os.path.basename(self.fname3d.replace('.vtk','.png'))
if make_eps:
epsfile = pngfile.replace('.png','.eps')
self.save_eps(html_dir,epsfile)
self.save_png(html_dir,pngfile,image_size)
images.append(['%s=%s'%(dir,slice),pngfile])
self.mv.quit() # don't litter with open mayavi instances
rtime = time.time() - start
print 'Rendering time: %4.2f seconds for %s slices.' % \
(rtime ,len(images) )
# If we didn't fail in the above, proceed with the HTML generation
doc = html.SimpleDocument(title=fname)
doc.append(html.Heading(1,fname,align='center'))
# now make a table with the images
tab = html.Table('Slices along %s direction' % dir,
column1_align='center',cell_align='center')
tab.body = []
col = 1; l1 = []; l2 = []
for img in images:
#print 'col',col,'image',img # dbg
l1.append(html.Image(img[1]))
l2.append(img[0])
if col%3:
col += 1
else:
#print 'finished row' # dbg
tab.body.extend([l1,l2])
col = 1; l1 = []; l2 = []
if col>1:
tab.body.extend([l1,l2]) # don't forget last (incomplete) row
doc.append(tab)
doc.write(os.path.join(html_dir,'index.html'))
#---------------------------------------------------------------------------
class Eigenmode(Data4d):
def __init__(self,fname,mode_num,slice_dir='x',slice_val=0):
"""Store and visualize fermionic eigenmode information.
The main data attributes are:
data_p1: <psi_bar*(gamma_5+1)*psi>
data_p2: <psi_bar*(gamma_5-1)*psi>
data_chi: (p1+p2)/2 = <psi_bar*gamma_5*psi>
data_psi: (p1-p2)/2 = <psi_bar*psi>
The default data attribute is initialized to data_chi.
"""
self.read(fname,mode_num)
self.data = self.data_chi # by default plot the chiral density
self.fname = fname
self.mode_num = mode_num
self.cut_slice(slice_dir,slice_val) # creates an actual 3d slice
self.clipped = 0
self.opacity = [(0,1),(0.3,0.8),(1,0)] # for volume rendering
def view(self,type=None,slice_dir=None,slice_val=None,
iso_val0=0.25,iso_val1=None,check_histo=1,
contours=-1,invert=0,iso_mode='relative_3',legend=1):
"""View either chirality(+/-), chiral density or fermion density.
Specify type as 'p1','p2','chi' or 'psi'.
If no type is specified, the currently active slice is viewed."""
new_slice = 0
if slice_dir is None:
slice_dir = self.slice_dir
new_slice = 1
if slice_val is None:
slice_val = self.slice_val
new_slice = 1
if type is not None:
data = 'data_' + type
self.data = getattr(self,data)
new_slice = 1
if new_slice:
self.cut_slice(slice_dir,slice_val)
Data4d.view(self,iso_val0,iso_val1,check_histo,contours,invert,
iso_mode,legend)
def read(self,fname,mode_num):
"""Read in a text density file.
Assumes nx=ny=nz (but not necessarily =nt)."""
for line in file(fname):
if line.startswith('nx'):
nx = ny = nz = int(line.strip().split()[1])
if line.startswith('nt'):
nt = int(line.strip().split()[1])
break
print 'Reading data with nx=%s nt=%s' % (nx,nt)
# a -> g5+1
cmd = "awk '/^%sPG5P/ { print $6}' %s" % (mode_num,fname)
tfile = os.popen(cmd)
data_p1_flat = array(map(float,tfile))
data_p1_flat.shape = (nt,nz,ny,nx)
# b -> g5-1
cmd = "awk '/^%sPG5P/ { print $7}' %s" % (mode_num,fname)
tfile = os.popen(cmd)
data_p2_flat = array(map(float,tfile))
data_p2_flat.shape = (nt,nz,ny,nx)
with(self,nx=nx,ny=ny,nz=nz,nt=nt,
data_p1=data_p1_flat,
data_p2=data_p2_flat,
data_chi = 0.5*(data_p1_flat+data_p2_flat),
data_psi = 0.5*(data_p1_flat-data_p2_flat),)
data_size = nx*ny*nz*nt
print 'Read',data_size,'data points. Shape:',(nx,ny,nz,nt)
def histo4_all_make(self):
"""Make histograms of the full 4d data.
Makes histo4_<type>, where <type> is in [p1,p2,chi,psi].
You can then view them with histo4_plot."""
try:
self.histo4_all_made
except AttributeError:
pass
else:
return
nbins = 150
print 'Building histograms of full dataset, this may take a moment...',
sys.stdout.flush()
try:
self.histo4_p1 = Histogram(ravel(self.data_p1),nbins)
except:
warn("Error making p1 histogram")
try:
self.histo4_p2 = Histogram(ravel(self.data_p2),nbins)
except:
warn("Error making p2 histogram")
try:
self.histo4_chi = Histogram(ravel(self.data_chi),nbins)
except:
warn("Error making chi histogram")
try:
self.histo4_psi = Histogram(ravel(self.data_psi),nbins)
except:
warn("Error making psi histogram")
self.histo4_all_made = 1
print 'Done.'
def histo4_all_plot(self,plotter=None):
"""Plot the 4-d histograms.
If they haven't been made, it builds them automatically."""
self.histo4_all_make()
plotter = gracePlot.gracePlot()
plotter.multi(2,2,0.1,0.2,0.4) # 2x2 plot table
plotter.title('p1: <psi_b(g5+1)psi>')
if hasattr(self,'histo4_p1'):
plotter.histoPlot(self.histo4_p1)
plotter.focus(0,1)
plotter.title('p2: <psi_b(g5-1)psi>')
if hasattr(self,'histo4_p2'):
plotter.histoPlot(self.histo4_p2)
plotter.focus(1,0)
plotter.title('chi: <psi_b*g5*psi>')
if hasattr(self,'histo4_chi'):
plotter.histoPlot(self.histo4_chi)
plotter.focus(1,1)
plotter.title('psi: <psi_b*g5*psi>')
if hasattr(self,'histo4_psi'):
plotter.histoPlot(self.histo4_psi)
def slice_table(self,type='chi',dir='x',step=1,
iso_val0=0.25,iso_val1=None,iso_mode='relative_3',
contours=-1,invert=0,
out_dir='',image_size=350,make_eps=0):
"""Make an HTML page with a table of 3-d slices.
"""
data = 'data_' + type
self.data = getattr(self,data)
fname = os.path.basename(self.fname) + \
'_mode_%s_%s' % (self.mode_num,type)
Data4d.slice_table(self,dir,step,fname,iso_val0,iso_val1,iso_mode,
contours,invert,
out_dir,image_size,make_eps)
#---------------------------------------------------------------------------
class TopologicalDensity(Data4d):
def __init__(self,fname,slice_dir='x',slice_val=0):
self.read(fname)
self.fname = fname
self.cut_slice(slice_dir,slice_val)
self.clipped = 0
def read(self,fname):
"""Read in a topological density file.
Does NOT do byte-reversal.
"""
# MILC code:
TOPO_VERSION = 66051
# of bytes to skip in file *from the beginning* before data starts
offset = 108
topo = open(fname)
topo_version,nx,ny,nz,nt = struct.unpack('i'*5,topo.read(5*4))
if topo_version != TOPO_VERSION:
warn('Wrong topological density version number (or '
'byte-reversed file)!\n'
'Topo_version = %s, should be %s' %
(topo_version,TOPO_VERSION) , level=4)
topo.seek(offset)
data_size = nx*ny*nz*nt*4 # these files use C float, not double
topo_flat = fromstring(topo.read(data_size),Float32)
topo_flat.shape = (nt,nz,ny,nx)
with(self,nx=nx,ny=ny,nz=nz,nt=nt,data=topo_flat)
print 'Read',data_size/4,'data points. Shape:',(nx,ny,nz,nt)
-------------------------------------------------------------------------
Take Surveys. Earn Cash. Influence the Future of IT
Join SourceForge.net's Techsay panel and you'll get the chance to share your
opinions on IT & business topics through brief surveys - and earn cash
http://www.techsay.com/default.php?page=join.php&p=sourceforge&CID=DEVDEV
_______________________________________________
MayaVi-users mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/mayavi-users
