[Yt-svn] yt-commit r538 - in trunk/yt: lagos raven
mturk at wrangler.dreamhost.com
mturk at wrangler.dreamhost.com
Thu Jun 5 00:10:53 PDT 2008
Author: mturk
Date: Thu Jun 5 00:10:52 2008
New Revision: 538
URL: http://yt.spacepope.org/changeset/538
Log:
Added support for 3D profiles, and displaying of these profiles through S2PLOT
in volume rendering format.
Axes are labeled, but it doesn't yet have a colorbar.
See: http://yt.enzotools.org/attachment/wiki/Screenshots/3DBinnedProfile.png
for an example.
The API is still in flux, and I'm still trying to figure out the best way to
present this functionality to the user. However, it's currently already
useful...
This closes #115.
Modified:
trunk/yt/lagos/PointCombine.c
trunk/yt/lagos/Profiles.py
trunk/yt/raven/Plot3DInterface.py
Modified: trunk/yt/lagos/PointCombine.c
==============================================================================
--- trunk/yt/lagos/PointCombine.c (original)
+++ trunk/yt/lagos/PointCombine.c Thu Jun 5 00:10:52 2008
@@ -383,6 +383,137 @@
}
+static PyObject *_profile3DError;
+
+static PyObject *Py_Bin3DProfile(PyObject *obj, PyObject *args)
+{
+ int i, j, k;
+ PyObject *obins_x, *obins_y, *obins_z, *owsource, *obsource, *owresult,
+ *obresult, *oused;
+ PyArrayObject *bins_x, *bins_y, *bins_z, *wsource, *bsource, *wresult,
+ *bresult, *used;
+ bins_x = bins_y = bins_z = wsource = bsource = wresult = bresult = used = NULL;
+
+ if (!PyArg_ParseTuple(args, "OOOOOOOO",
+ &obins_x, &obins_y, &obins_z, &owsource, &obsource,
+ &owresult, &obresult, &oused))
+ return PyErr_Format(_profile3DError,
+ "Bin3DProfile: Invalid parameters.");
+ i = 0;
+
+ bins_x = (PyArrayObject *) PyArray_FromAny(obins_x,
+ PyArray_DescrFromType(NPY_INT64), 1, 1,
+ NPY_IN_ARRAY, NULL);
+ if(bins_x==NULL) {
+ PyErr_Format(_profile3DError,
+ "Bin3DProfile: One dimension required for bins_x.");
+ goto _fail;
+ }
+
+ bins_y = (PyArrayObject *) PyArray_FromAny(obins_y,
+ PyArray_DescrFromType(NPY_INT64), 1, 1,
+ NPY_IN_ARRAY, NULL);
+ if((bins_y==NULL) || (PyArray_SIZE(bins_x) != PyArray_SIZE(bins_y))) {
+ PyErr_Format(_profile3DError,
+ "Bin3DProfile: One dimension required for bins_y, same size as bins_x.");
+ goto _fail;
+ }
+
+ bins_z = (PyArrayObject *) PyArray_FromAny(obins_z,
+ PyArray_DescrFromType(NPY_INT64), 1, 1,
+ NPY_IN_ARRAY, NULL);
+ if((bins_z==NULL) || (PyArray_SIZE(bins_x) != PyArray_SIZE(bins_z))) {
+ PyErr_Format(_profile3DError,
+ "Bin3DProfile: One dimension required for bins_z, same size as bins_x.");
+ goto _fail;
+ }
+
+ wsource = (PyArrayObject *) PyArray_FromAny(owsource,
+ PyArray_DescrFromType(NPY_FLOAT64), 1, 1,
+ NPY_IN_ARRAY, NULL);
+ if((wsource==NULL) || (PyArray_SIZE(bins_x) != PyArray_SIZE(wsource))) {
+ PyErr_Format(_profile3DError,
+ "Bin3DProfile: One dimension required for wsource, same size as bins_x.");
+ goto _fail;
+ }
+
+ bsource = (PyArrayObject *) PyArray_FromAny(obsource,
+ PyArray_DescrFromType(NPY_FLOAT64), 1, 1,
+ NPY_IN_ARRAY, NULL);
+ if((bsource==NULL) || (PyArray_SIZE(bins_x) != PyArray_SIZE(bsource))) {
+ PyErr_Format(_profile3DError,
+ "Bin3DProfile: One dimension required for bsource, same size as bins_x.");
+ goto _fail;
+ }
+
+ wresult = (PyArrayObject *) PyArray_FromAny(owresult,
+ PyArray_DescrFromType(NPY_FLOAT64), 3,3,
+ NPY_INOUT_ARRAY | NPY_UPDATEIFCOPY, NULL);
+ if(wresult==NULL){
+ PyErr_Format(_profile3DError,
+ "Bin3DProfile: Two dimensions required for wresult.");
+ goto _fail;
+ }
+
+ bresult = (PyArrayObject *) PyArray_FromAny(obresult,
+ PyArray_DescrFromType(NPY_FLOAT64), 3,3,
+ NPY_INOUT_ARRAY | NPY_UPDATEIFCOPY, NULL);
+ if((bresult==NULL) ||(PyArray_SIZE(wresult) != PyArray_SIZE(bresult))
+ || (PyArray_DIM(bresult,0) != PyArray_DIM(wresult,0))){
+ PyErr_Format(_profile3DError,
+ "Bin3DProfile: Two dimensions required for bresult, same shape as wresult.");
+ goto _fail;
+ }
+
+ used = (PyArrayObject *) PyArray_FromAny(oused,
+ PyArray_DescrFromType(NPY_FLOAT64), 3,3,
+ NPY_INOUT_ARRAY | NPY_UPDATEIFCOPY, NULL);
+ if((used==NULL) ||(PyArray_SIZE(used) != PyArray_SIZE(wresult))
+ || (PyArray_DIM(used,0) != PyArray_DIM(wresult,0))){
+ PyErr_Format(_profile3DError,
+ "Bin3DProfile: Two dimensions required for used, same shape as wresult.");
+ goto _fail;
+ }
+
+ npy_float64 wval, bval;
+ int n;
+
+ for(n=0; n<bins_x->dimensions[0]; n++) {
+ i = *(npy_int64*)PyArray_GETPTR1(bins_x, n);
+ j = *(npy_int64*)PyArray_GETPTR1(bins_y, n);
+ k = *(npy_int64*)PyArray_GETPTR1(bins_z, n);
+ bval = *(npy_float64*)PyArray_GETPTR1(bsource, n);
+ wval = *(npy_float64*)PyArray_GETPTR1(wsource, n);
+ *(npy_float64*)PyArray_GETPTR3(wresult, i, j, k) += wval;
+ *(npy_float64*)PyArray_GETPTR3(bresult, i, j, k) += wval*bval;
+ *(npy_float64*)PyArray_GETPTR3(used, i, j, k) = 1.0;
+ }
+
+ Py_DECREF(bins_x);
+ Py_DECREF(bins_y);
+ Py_DECREF(bins_z);
+ Py_DECREF(wsource);
+ Py_DECREF(bsource);
+ Py_DECREF(wresult);
+ Py_DECREF(bresult);
+ Py_DECREF(used);
+
+ PyObject *onum_found = PyInt_FromLong((long)1);
+ return onum_found;
+
+ _fail:
+ Py_XDECREF(bins_x);
+ Py_XDECREF(bins_y);
+ Py_XDECREF(bins_z);
+ Py_XDECREF(wsource);
+ Py_XDECREF(bsource);
+ Py_XDECREF(wresult);
+ Py_XDECREF(bresult);
+ Py_XDECREF(used);
+ return NULL;
+
+}
+
static PyObject *_dataCubeError;
static PyObject *DataCubeGeneric(PyObject *obj, PyObject *args,
@@ -952,6 +1083,7 @@
{"DataCubeRefine", Py_DataCubeRefine, METH_VARARGS},
{"DataCubeReplace", Py_DataCubeReplace, METH_VARARGS},
{"Bin2DProfile", Py_Bin2DProfile, METH_VARARGS},
+ {"Bin3DProfile", Py_Bin3DProfile, METH_VARARGS},
{"FindContours", Py_FindContours, METH_VARARGS},
{"FindBindingEnergy", Py_FindBindingEnergy, METH_VARARGS},
{NULL, NULL} /* Sentinel */
@@ -975,6 +1107,8 @@
PyDict_SetItemString(d, "error", _dataCubeError);
_profile2DError = PyErr_NewException("PointCombine.Profile2DError", NULL, NULL);
PyDict_SetItemString(d, "error", _profile2DError);
+ _profile3DError = PyErr_NewException("PointCombine.Profile3DError", NULL, NULL);
+ PyDict_SetItemString(d, "error", _profile3DError);
_findContoursError = PyErr_NewException("PointCombine.FindContoursError", NULL, NULL);
PyDict_SetItemString(d, "error", _findContoursError);
import_array();
Modified: trunk/yt/lagos/Profiles.py
==============================================================================
--- trunk/yt/lagos/Profiles.py (original)
+++ trunk/yt/lagos/Profiles.py Thu Jun 5 00:10:52 2008
@@ -104,7 +104,6 @@
def __setitem__(self, key, value):
self._data[key] = value
-
# @todo: Fix accumulation with overriding
class BinnedProfile1D(BinnedProfile):
def __init__(self, data_source, n_bins, bin_field,
@@ -186,26 +185,20 @@
x_n_bins, x_bin_field, x_lower_bound, x_upper_bound, x_log,
y_n_bins, y_bin_field, y_lower_bound, y_upper_bound, y_log,
lazy_reader=False):
-
- self.total_cells = 0
BinnedProfile.__init__(self, data_source, lazy_reader)
self.x_bin_field = x_bin_field
self.y_bin_field = y_bin_field
self._x_log = x_log
self._y_log = y_log
- if x_log:
- self[x_bin_field] = na.logspace(na.log10(x_lower_bound*0.99),
- na.log10(x_upper_bound*1.01),
- x_n_bins)
- else:
- self[x_bin_field] = na.linspace(
+ if x_log: self[x_bin_field] = na.logspace(na.log10(x_lower_bound*0.99),
+ na.log10(x_upper_bound*1.01),
+ x_n_bins)
+ else: self[x_bin_field] = na.linspace(
x_lower_bound*0.99, x_upper_bound*1.01, x_n_bins)
- if y_log:
- self[y_bin_field] = na.logspace(na.log10(y_lower_bound*0.99),
- na.log10(y_upper_bound*1.01),
- y_n_bins)
- else:
- self[y_bin_field] = na.linspace(
+ if y_log: self[y_bin_field] = na.logspace(na.log10(y_lower_bound*0.99),
+ na.log10(y_upper_bound*1.01),
+ y_n_bins)
+ else: self[y_bin_field] = na.linspace(
y_lower_bound*0.99, y_upper_bound*1.01, y_n_bins)
if na.any(na.isnan(self[x_bin_field])) \
or na.any(na.isnan(self[y_bin_field])):
@@ -240,7 +233,6 @@
bin_indices_y = args[2].ravel().astype('int64')
source_data = source_data[mi]
weight_data = weight_data[mi]
- self.total_cells += bin_indices_x.size
nx = bin_indices_x.size
#mylog.debug("Binning %s / %s times", source_data.size, nx)
PointCombine.Bin2DProfile(bin_indices_x, bin_indices_y, weight_data, source_data,
@@ -293,3 +285,115 @@
field_data[:,line].tofile(fid, sep="\t", format=format)
fid.write("\n")
fid.close()
+
+class BinnedProfile3D(BinnedProfile):
+ def __init__(self, data_source,
+ x_n_bins, x_bin_field, x_lower_bound, x_upper_bound, x_log,
+ y_n_bins, y_bin_field, y_lower_bound, y_upper_bound, y_log,
+ z_n_bins, z_bin_field, z_lower_bound, z_upper_bound, z_log,
+ lazy_reader=False):
+ BinnedProfile.__init__(self, data_source, lazy_reader)
+ self.x_bin_field = x_bin_field
+ self.y_bin_field = y_bin_field
+ self.z_bin_field = z_bin_field
+ self._x_log = x_log
+ self._y_log = y_log
+ self._z_log = z_log
+ if x_log: self[x_bin_field] = na.logspace(na.log10(x_lower_bound*0.99),
+ na.log10(x_upper_bound*1.01),
+ x_n_bins)
+ else: self[x_bin_field] = na.linspace(
+ x_lower_bound*0.99, x_upper_bound*1.01, x_n_bins)
+ if y_log: self[y_bin_field] = na.logspace(na.log10(y_lower_bound*0.99),
+ na.log10(y_upper_bound*1.01),
+ y_n_bins)
+ else: self[y_bin_field] = na.linspace(
+ y_lower_bound*0.99, y_upper_bound*1.01, y_n_bins)
+ if z_log: self[z_bin_field] = na.logspace(na.log10(z_lower_bound*0.99),
+ na.log10(z_upper_bound*1.01),
+ z_n_bins)
+ else: self[z_bin_field] = na.linspace(
+ z_lower_bound*0.99, z_upper_bound*1.01, z_n_bins)
+ if na.any(na.isnan(self[x_bin_field])) \
+ or na.any(na.isnan(self[y_bin_field])) \
+ or na.any(na.isnan(self[z_bin_field])):
+ mylog.error("Your min/max values for x, y or z have given me a nan.")
+ mylog.error("Usually this means you are asking for log, with a zero bound.")
+ raise ValueError
+ if not lazy_reader:
+ self._args = self._get_bins(data_source)
+
+ def _get_empty_field(self):
+ return na.zeros((self[self.x_bin_field].size,
+ self[self.y_bin_field].size,
+ self[self.z_bin_field].size), dtype='float64')
+
+ @preserve_source_parameters
+ def _bin_field(self, source, field, weight, accumulation,
+ args, check_cut=False):
+ if check_cut:
+ pointI = self._data_source._get_point_indices(source)
+ source_data = source[field][pointI].ravel().astype('float64')
+ weight_data = na.ones(source_data.shape).astype('float64')
+ if weight: weight_data = source[weight][pointI].ravel().astype('float64')
+ else:
+ source_data = source[field].ravel().astype('float64')
+ weight_data = na.ones(source_data.shape).astype('float64')
+ if weight: weight_data = source[weight].ravel().astype('float64')
+ self.total_stuff = source_data.sum()
+ binned_field = self._get_empty_field()
+ weight_field = self._get_empty_field()
+ used_field = self._get_empty_field()
+ mi = args[0]
+ bin_indices_x = args[1].ravel().astype('int64')
+ bin_indices_y = args[2].ravel().astype('int64')
+ bin_indices_z = args[3].ravel().astype('int64')
+ source_data = source_data[mi]
+ weight_data = weight_data[mi]
+ PointCombine.Bin3DProfile(
+ bin_indices_x, bin_indices_y, bin_indices_z,
+ weight_data, source_data,
+ weight_field, binned_field, used_field)
+ if accumulation: # Fix for laziness
+ if not iterable(accumulation):
+ raise SyntaxError("Accumulation needs to have length 2")
+ if accumulation[0]:
+ binned_field = na.add.accumulate(binned_field, axis=0)
+ if accumulation[1]:
+ binned_field = na.add.accumulate(binned_field, axis=1)
+ if accumulation[2]:
+ binned_field = na.add.accumulate(binned_field, axis=2)
+ return binned_field, weight_field, used_field.astype('bool')
+
+ @preserve_source_parameters
+ def _get_bins(self, source, check_cut=False):
+ if check_cut:
+ cm = self._data_source._get_point_indices(source)
+ source_data_x = source[self.x_bin_field][cm]
+ source_data_y = source[self.y_bin_field][cm]
+ source_data_z = source[self.z_bin_field][cm]
+ else:
+ source_data_x = source[self.x_bin_field]
+ source_data_y = source[self.y_bin_field]
+ source_data_z = source[self.z_bin_field]
+ if source_data_x.size == 0:
+ return
+ mi = na.where( (source_data_x > self[self.x_bin_field].min())
+ & (source_data_x < self[self.x_bin_field].max())
+ & (source_data_y > self[self.y_bin_field].min())
+ & (source_data_y < self[self.y_bin_field].max())
+ & (source_data_z > self[self.z_bin_field].min())
+ & (source_data_z < self[self.z_bin_field].max()))
+ sd_x = source_data_x[mi]
+ sd_y = source_data_y[mi]
+ sd_z = source_data_z[mi]
+ if sd_x.size == 0 or sd_y.size == 0 or sd_z.size == 0:
+ return
+ bin_indices_x = na.digitize(sd_x, self[self.x_bin_field])
+ bin_indices_y = na.digitize(sd_y, self[self.y_bin_field])
+ bin_indices_z = na.digitize(sd_z, self[self.z_bin_field])
+ # Now we set up our inverse bin indices
+ return (mi, bin_indices_x, bin_indices_y, bin_indices_z)
+
+ def write_out(self, filename, format="%0.16e"):
+ pass # Will eventually dump HDF5
Modified: trunk/yt/raven/Plot3DInterface.py
==============================================================================
--- trunk/yt/raven/Plot3DInterface.py (original)
+++ trunk/yt/raven/Plot3DInterface.py Thu Jun 5 00:10:52 2008
@@ -38,38 +38,36 @@
func(obj, *args, **kwargs)
return check_started
-class VolumeRenderingCube(object):
- def __init__(self, pf, center=None, width=1, unit='1',
- field='Density', dims=128, take_log=True,
+class VolumeRendering(object):
+ def __init__(self, data, take_log=True,
window_opts="/S2MONO", cmap="rainbow",
- amin=0.01, amax=0.1):
- self.pf = pf
- self.width = width/pf[unit]
- if center is None: center = pf.h.find_max("Density")[1]
- self.center = center
- self.field = field
- self.dims = dims
- dx = self.width / dims
- self.max_level = na.unique(pf.h.gridDxs[pf.h.gridDxs>=dx]).argmax()+1
- self.__setup_data(take_log)
+ amin=0.0, amax=0.1, bounds = None):
+ if bounds is None:
+ self.x0,self.x1, self.y0,self.y1, self.z0,self.z1 = \
+ 0.0,1.0, 0.0,1.0, 0.0,1.0
+ else:
+ self.x0,self.x1, self.y0,self.y1, self.z0,self.z1 = \
+ bounds
+ self.__setup_data(data, take_log)
self.vrid = None
self.isoids = []
self.window_opts = window_opts
self.cmap = cmap
self._amin, self._amax = amin, amax
- self.tr = na.array([ 0.0, (1.0-0.0)/(self.dims-1.0), 0, 0,
- 0.0, 0, (1.0-0.0)/(self.dims-1.0), 0,
- 0.0, 0, 0, (1.0-0.0)/(self.dims-1.0) ])
+ nx, ny, nz = self.dims
+ self.tr = na.array([
+ self.x0, (self.x1-self.x0)/(nx-1.0), 0 , 0,
+ self.y0, 0 , (self.y1-self.y0)/(ny-1.0), 0,
+ self.z0, 0 , 0 , (self.z1-self.z0)/(nz-1.0)
+ ])
self.started = False
- def __setup_data(self, take_log):
- self.data_grid = self.pf.h.covering_grid(self.max_level,
- self.center - self.width/2.0,
- self.center + self.width/2.0,
- [self.dims]*3, fields=[self.field])
- if take_log: self.data=na.log10(self.data_grid[self.field])
- else: self.data=self.data_grid[self.field]
+ def __setup_data(self, data, take_log):
+ self.dims = na.array(data.shape)
+ if take_log: self.data=na.log10(data).copy()
+ else: self.data=data.copy()
self._dmin = self.data.min()
+ if take_log: self._dmin=na.log10(data[data>0]).min()
self._dmax = self.data.max()
@must_have_s2plot
@@ -95,27 +93,36 @@
def __add_isosurface(self, val, alpha, cm):
scaled_val = ((val-self._dmin)/(self._dmax-self._dmin))
r,g,b,a = cm(scaled_val)
- nx = self.dims-1
+ nx,ny,nz = self.dims
print "Adding",val,scaled_val,alpha,r,g,b
- return s2plot.ns2cis(self.data, self.dims, self.dims, self.dims,
- 0, nx, 0, nx, 0, nx, self.tr, val,
+ return s2plot.ns2cis(self.data, nx, ny, nz,
+ 0, nx-1, 0, ny-1, 0, ny-1,
+ self.tr, val,
1, 't', alpha, r,g,b)
def run(self, pre_call=None):
self.__setup_s2plot()
self.__setup_volrendering()
self.__register_callbacks()
+ self._setup_labels()
if pre_call is not None: pre_call(self)
- self.__start_rendering()
+ s2plot.s2disp(-1, 1)
def restart(self):
- self.__start_rendering()
+ s2plot.s2disp(-1, 0)
+
+ def _setup_labels(self):
+ pass
def __setup_s2plot(self):
- dx = 1.0/(self.dims-1.0)
+ dx,dy,dz = 1.0/(self.dims-1.0)
s2plot.s2opendo(self.window_opts)
- s2plot.s2swin(-dx,1.0+dx,-dx,1.0+dx,-dx,1.0+dx) # We mandate cube
- s2plot.s2box("BCDE",0,0,"BCDE",0,0,"BCDE",0,0)
+ s2plot.s2swin(-dx/1.0,1.0+dx/1.0,
+ -dy/1.0,1.0+dy/1.0,
+ -dz/1.0,1.0+dz/1.0)
+ #opts = "BCDETMNOPQ"
+ opts = "BCDE"
+ s2plot.s2box(opts,0,0,opts,0,0,opts,0,0)
s2plot.s2scir(1000,2000) # Set colour range
s2plot.s2icm(self.cmap,1000,2000) # Install colour map
amb = {'r':0.8, 'g':0.8, 'b':0.8} # ambient light
@@ -124,9 +131,9 @@
self.started = True
def __setup_volrendering(self):
- nd = self.dims
- self.vrid = s2plot.ns2cvr(self.data, nd, nd, nd,
- 0, nd-1, 0, nd-1, 0, nd-1,
+ ndx,ndy,ndz = self.dims
+ self.vrid = s2plot.ns2cvr(self.data, ndx, ndy, ndz,
+ 0, ndx-1, 0, ndy-1, 0, ndz-1,
self.tr, 's',
self._dmin, self._dmax, self._amin, self._amax)
@@ -134,8 +141,50 @@
# More should go here for functional changes to the object
s2plot.cs2scb(self.__my_callback) # Install a dynamic callback
- def __start_rendering(self):
- s2plot.s2disp(-1, 1)
-
def __my_callback(self, t, kc):
s2plot.ds2dvr(self.vrid, 0)
+
+class VolumeRenderingDataCube(VolumeRendering):
+ def __init__(self, pf, center=None, width=1, unit='1',
+ field='Density', dims=128, **kwargs):
+ self.pf = pf
+ self.width = width/pf[unit]
+ if center is None: center = pf.h.find_max("Density")[1]
+ self.center = center
+ self.field = field
+ self.dims = dims
+ dx = self.width / dims
+ self.max_level = na.unique(pf.h.gridDxs[pf.h.gridDxs>=dx]).argmax()+1
+ self.data_grid = self.__get_data()
+ VolumeRendering.__init__(self, self.data_grid[field], **kwargs)
+
+ def __get_data(self):
+ data_grid = self.pf.h.covering_grid(self.max_level,
+ self.center - self.width/2.0,
+ self.center + self.width/2.0,
+ [self.dims]*3, fields=[self.field])
+ return data_grid
+
+class VolumeRendering3DProfile(VolumeRendering):
+ def __init__(self, profile, field, **kwargs):
+ self.profile = profile
+ self.field = field
+ if 'bounds' not in kwargs:
+ kwargs['bounds'] = self._setup_bounds()
+ VolumeRendering.__init__(self, self.profile[field], **kwargs)
+
+ def _setup_bounds(self):
+ nolog = lambda a: a
+ self.bf = [[nolog, self.profile.x_bin_field], \
+ [nolog, self.profile.y_bin_field], \
+ [nolog, self.profile.z_bin_field]]
+ if self.profile._x_log: self.bf[0][0] = na.log10
+ if self.profile._y_log: self.bf[1][0] = na.log10
+ if self.profile._z_log: self.bf[2][0] = na.log10
+ xmi, ymi, zmi = [f(self.profile[fn]).min() for f,fn in self.bf]
+ xma, yma, zma = [f(self.profile[fn]).max() for f,fn in self.bf]
+ return (xmi,xma,ymi,yma,zmi,zma)
+
+ def _setup_labels(self):
+ s2plot.s2env(self.x0,self.x1, self.y0,self.y1, self.z0,self.z1, 0, 1)
+ s2plot.s2lab(self.bf[0][1],self.bf[1][1],self.bf[2][1],self.field)
More information about the yt-svn
mailing list