[Yt-svn] yt-commit r424 - trunk/yt/raven
mturk at wrangler.dreamhost.com
mturk at wrangler.dreamhost.com
Thu May 1 14:34:21 PDT 2008
Author: mturk
Date: Thu May 1 14:34:21 2008
New Revision: 424
URL: http://yt.spacepope.org/changeset/424
Log:
Moved the cutting plane pixelization into C.
Modified:
trunk/yt/raven/PlotTypes.py
trunk/yt/raven/_MPL.c
Modified: trunk/yt/raven/PlotTypes.py
==============================================================================
--- trunk/yt/raven/PlotTypes.py (original)
+++ trunk/yt/raven/PlotTypes.py Thu May 1 14:34:21 2008
@@ -401,23 +401,14 @@
px_min, px_max = self.xlim
py_min, py_max = self.ylim
l, b, width, height = self._axes.bbox.get_bounds()
- pxs, pys, pzs = self.data['px'], self.data['py'], self.data['pz']
- xs, ys, zs = self.data['x'], self.data['y'], self.data['z']
- dxs, dys, dzs = self.data['pdx'], self.data['pdy'], self.data['pdz']
- field = self.axis_names['Z']
- ds = self.data[field]
- indices = na.argsort(dxs)[::-1]
- nx = indices.size
- inv_mat = self.data._inv_mat
- center = na.array(self.data.center)
- buff = na.zeros((width,height), dtype='float64')
- count = na.zeros((width,height), dtype='float64')
- weave.inline(_pixelize_cp,
- ['pxs','pys','pzs','xs','ys','zs','dxs','dys','dzs',
- 'buff','ds','nx','inv_mat','width','height','count',
- 'px_min','px_max','py_min','py_max', 'center', 'indices'],
- compiler='gcc', type_converters=converters.blitz,
- auto_downcast = 0, verbose=2)
+ indices = na.argsort(self.data['dx'])[::-1]
+ buff = _MPL.CPixelize( self.data['x'], self.data['y'], self.data['z'],
+ self.data['px'], self.data['py'],
+ self.data['pdx'], self.data['pdy'], self.data['pdz'],
+ self.data.center, self.data._inv_mat, indices,
+ self.data[self.axis_names['Z']],
+ int(width), int(height),
+ (px_min, px_max, py_min, py_max))
return buff
def _refresh_display_width(self, width=None):
@@ -799,51 +790,3 @@
plot._axes.set_ylim(yy0,yy1)
plot._axes.hold(False)
return runCallback
-
-_pixelize_cp = r"""
-
-long double md, cxpx, cypx;
-long double cx, cy, cz;
-long double lrx, lry, lrz;
-long double rrx, rry, rrz;
-int lc, lr, rc, rr, p;
-
-long double px_dx, px_dy, px_dz, overlap1, overlap2, overlap3;
-px_dx = (px_max-px_min)/height;
-px_dy = (py_max-py_min)/width;
-px_dz = sqrt(0.5 * (px_dy*px_dy + px_dx*px_dx));
-
-#define min(X,Y) ((X) < (Y) ? (X) : (Y))
-#define max(X,Y) ((X) > (Y) ? (X) : (Y))
-using namespace std;
-
-//for(int i=0; i<width; i++) for(int j=0; j<height; j++) count(i,j)=buff(i,j)=0.0;
-
-for(int pp=0; pp<nx; pp++)
-{
- p = indices(pp);
- // Any point we want to plot is at most this far from the center
- md = 2.0*sqrtl(dxs(p)*dxs(p) + dys(p)*dys(p) + dzs(p)*dzs(p));
- if(((pxs(p)+md<px_min) ||
- (pxs(p)-md>px_max)) ||
- ((pys(p)+md<py_min) ||
- (pys(p)-md>py_max))) continue;
- lc = max(floorl((pxs(p)-md-px_min)/px_dx),0);
- lr = max(floorl((pys(p)-md-py_min)/px_dy),0);
- rc = min(ceill((pxs(p)+md-px_min)/px_dx),height);
- rr = min(ceill((pys(p)+md-py_min)/px_dy),width);
- for (int i=lr;i<rr;i++) {
- cypx = px_dy * (i+0.5) + py_min;
- for (int j=lc;j<rc;j++) {
- cxpx = px_dx * (j+0.5) + px_min;
- cx = inv_mat(0,0)*cxpx + inv_mat(0,1)*cypx + center(0);
- cy = inv_mat(1,0)*cxpx + inv_mat(1,1)*cypx + center(1);
- cz = inv_mat(2,0)*cxpx + inv_mat(2,1)*cypx + center(2);
- if( ((xs(p)-cx)>1.01*dxs(p)) || ((xs(p)-cx)<(-1.01*dxs(p)))
- || ((ys(p)-cy)>1.01*dys(p)) || ((ys(p)-cy)<(-1.01*dys(p)))
- || ((zs(p)-cz)>1.01*dzs(p)) || ((zs(p)-cz)<(-1.01*dzs(p))) ) continue;
- buff(i,j) = ds(p);
- }
- }
-}
-"""
Modified: trunk/yt/raven/_MPL.c
==============================================================================
--- trunk/yt/raven/_MPL.c (original)
+++ trunk/yt/raven/_MPL.c Thu May 1 14:34:21 2008
@@ -181,9 +181,217 @@
}
+static PyObject* Py_CPixelize(PyObject *obj, PyObject *args) {
+
+ PyObject *xp, *yp, *zp, *pxp, *pyp,
+ *dxp, *dyp, *dzp, *dp,
+ *centerp, *inv_matp, *indicesp;
+ unsigned int rows, cols;
+ double px_min, px_max, py_min, py_max;
+
+ if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOII(dddd)",
+ &xp, &yp, &zp, &pxp, &pyp, &dxp, &dyp, &dzp, ¢erp, &inv_matp,
+ &indicesp, &dp, &cols, &rows, &px_min, &px_max, &py_min, &py_max))
+ return PyErr_Format(_pixelizeError, "CPixelize: Invalid Parameters.");
+
+ double width = px_max - px_min;
+ double height = py_max - py_min;
+ long double px_dx = width / ((double) rows);
+ long double px_dy = height / ((double) cols);
+
+ // Check we have something to output to
+ if (rows == 0 || cols ==0)
+ PyErr_Format( _pixelizeError, "Cannot scale to zero size.");
+
+ // Get numeric arrays
+ PyArrayObject *x = (PyArrayObject *) PyArray_FromAny(xp,
+ PyArray_DescrFromType(NPY_FLOAT64), 1, 1, NPY_C_CONTIGUOUS, NULL);
+ if (x == NULL) {
+ PyErr_Format( _pixelizeError, "x is of incorrect type (wanted 1D float)");
+ goto _fail;
+ }
+
+ PyArrayObject *y = (PyArrayObject *) PyArray_FromAny(yp,
+ PyArray_DescrFromType(NPY_FLOAT64), 1, 1, NPY_C_CONTIGUOUS, NULL);
+ if ((y == NULL) || (PyArray_SIZE(y) != PyArray_SIZE(x))) {
+ PyErr_Format( _pixelizeError, "y is of incorrect type (wanted 1D float)");
+ goto _fail;
+ }
+
+ PyArrayObject *z = (PyArrayObject *) PyArray_FromAny(zp,
+ PyArray_DescrFromType(NPY_FLOAT64), 1, 1, NPY_C_CONTIGUOUS, NULL);
+ if ((z == NULL) || (PyArray_SIZE(y) != PyArray_SIZE(x))) {
+ PyErr_Format( _pixelizeError, "z is of incorrect type (wanted 1D float)");
+ goto _fail;
+ }
+
+ PyArrayObject *px = (PyArrayObject *) PyArray_FromAny(pxp,
+ PyArray_DescrFromType(NPY_FLOAT64), 1, 1, NPY_C_CONTIGUOUS, NULL);
+ if ((px == NULL) || (PyArray_SIZE(y) != PyArray_SIZE(x))) {
+ PyErr_Format( _pixelizeError, "px is of incorrect type (wanted 1D float)");
+ goto _fail;
+ }
+
+ PyArrayObject *py = (PyArrayObject *) PyArray_FromAny(pyp,
+ PyArray_DescrFromType(NPY_FLOAT64), 1, 1, NPY_C_CONTIGUOUS, NULL);
+ if ((py == NULL) || (PyArray_SIZE(y) != PyArray_SIZE(x))) {
+ PyErr_Format( _pixelizeError, "py is of incorrect type (wanted 1D float)");
+ goto _fail;
+ }
+
+ PyArrayObject *d = (PyArrayObject *) PyArray_FromAny(dp,
+ PyArray_DescrFromType(NPY_FLOAT64), 1, 1, NPY_C_CONTIGUOUS, NULL);
+ if ((d == NULL) || (PyArray_SIZE(d) != PyArray_SIZE(x))) {
+ PyErr_Format( _pixelizeError, "data is of incorrect type (wanted 1D float)");
+ goto _fail;
+ }
+
+ PyArrayObject *dx = (PyArrayObject *) PyArray_FromAny(dxp,
+ PyArray_DescrFromType(NPY_FLOAT64), 1, 1, NPY_C_CONTIGUOUS, NULL);
+ if ((dx == NULL) || (PyArray_SIZE(dx) != PyArray_SIZE(x))) {
+ PyErr_Format( _pixelizeError, "dx is of incorrect type (wanted 1D float)");
+ goto _fail;
+ }
+ PyArrayObject *dy = (PyArrayObject *) PyArray_FromAny(dyp,
+ PyArray_DescrFromType(NPY_FLOAT64), 1, 1, NPY_C_CONTIGUOUS, NULL);
+ if ((dy == NULL) || (PyArray_SIZE(dy) != PyArray_SIZE(x))) {
+ PyErr_Format( _pixelizeError, "dy is of incorrect type (wanted 1D float)");
+ goto _fail;
+ }
+ PyArrayObject *dz = (PyArrayObject *) PyArray_FromAny(dzp,
+ PyArray_DescrFromType(NPY_FLOAT64), 1, 1, NPY_C_CONTIGUOUS, NULL);
+ if ((dz == NULL) || (PyArray_SIZE(dz) != PyArray_SIZE(x))) {
+ PyErr_Format( _pixelizeError, "dz is of incorrect type (wanted 1D float)");
+ goto _fail;
+ }
+ PyArrayObject *center = (PyArrayObject *) PyArray_FromAny(centerp,
+ PyArray_DescrFromType(NPY_FLOAT64), 1, 1, NPY_C_CONTIGUOUS, NULL);
+ if ((dz == NULL) || (PyArray_SIZE(center) != 3)) {
+ PyErr_Format( _pixelizeError, "Center must have three points");
+ goto _fail;
+ }
+ PyArrayObject *inv_mat = (PyArrayObject *) PyArray_FromAny(inv_matp,
+ PyArray_DescrFromType(NPY_FLOAT64), 2, 2, NPY_C_CONTIGUOUS, NULL);
+ if ((inv_mat == NULL) || (PyArray_SIZE(inv_mat) != 9)) {
+ PyErr_Format( _pixelizeError, "inv_mat must be three by three");
+ goto _fail;
+ }
+ PyArrayObject *indices = (PyArrayObject *) PyArray_FromAny(indicesp,
+ PyArray_DescrFromType(NPY_INT64), 1, 1, NPY_C_CONTIGUOUS, NULL);
+ if ((indices == NULL) || (PyArray_SIZE(indices) != PyArray_SIZE(dx))) {
+ PyErr_Format( _pixelizeError, "indices must be same length as dx");
+ goto _fail;
+ }
+
+ // Check dimensions match
+ int nx = x->dimensions[0];
+
+ npy_float64 *gridded;
+ gridded = malloc(sizeof(npy_float64)*rows*cols);
+ if (gridded == NULL) {
+ PyErr_Format(_pixelizeError, "Could not allocate memory for image");
+ goto _fail;
+ }
+
+ // Calculate the pointer arrays to map input x to output x
+ int i, j, p;
+ int lc, lr, rc, rr;
+ long double md, cxpx, cypx;
+ long double cx, cy, cz;
+
+ npy_float64 *xs = (npy_float64 *) PyArray_GETPTR1(x, 0);
+ npy_float64 *ys = (npy_float64 *) PyArray_GETPTR1(y, 0);
+ npy_float64 *zs = (npy_float64 *) PyArray_GETPTR1(z, 0);
+ npy_float64 *pxs = (npy_float64 *) PyArray_GETPTR1(px, 0);
+ npy_float64 *pys = (npy_float64 *) PyArray_GETPTR1(py, 0);
+ npy_float64 *dxs = (npy_float64 *) PyArray_GETPTR1(dx, 0);
+ npy_float64 *dys = (npy_float64 *) PyArray_GETPTR1(dy, 0);
+ npy_float64 *dzs = (npy_float64 *) PyArray_GETPTR1(dz, 0);
+ npy_float64 *ds = (npy_float64 *) PyArray_GETPTR1(d, 0); // We check this above
+ npy_float64 *centers = (npy_float64 *) PyArray_GETPTR1(center,0);
+ npy_int64 *indicess = (npy_int64 *) PyArray_GETPTR1(indices,0);
+
+ npy_intp dims[] = {rows, cols};
+ PyArrayObject *my_array =
+ (PyArrayObject *) PyArray_SimpleNewFromDescr(2, dims,
+ PyArray_DescrFromType(NPY_FLOAT64));
+ gridded = (npy_float64 *) my_array->data;
+
+ npy_float64 inv_mats[3][3];
+ for(i=0;i<3;i++)for(j=0;j<3;j++)
+ inv_mats[i][j]=*(npy_float64*)PyArray_GETPTR2(inv_mat,i,j);
+
+ int pp;
+ for(p=0;p<cols*rows;p++)gridded[p]=0.0;
+ for(pp=0; pp<nx; pp++)
+ {
+ p = indicess[pp];
+ // Any point we want to plot is at most this far from the center
+ md = 2.0*sqrtl(dxs[p]*dxs[p] + dys[p]*dys[p] + dzs[p]*dzs[p]);
+ if(((pxs[p]+md<px_min) ||
+ (pxs[p]-md>px_max)) ||
+ ((pys[p]+md<py_min) ||
+ (pys[p]-md>py_max))) continue;
+ lc = max(floorl((pxs[p]-md-px_min)/px_dx),0);
+ lr = max(floorl((pys[p]-md-py_min)/px_dy),0);
+ rc = min(ceill((pxs[p]+md-px_min)/px_dx),rows);
+ rr = min(ceill((pys[p]+md-py_min)/px_dy),cols);
+ for (i=lr;i<rr;i++) {
+ cypx = px_dy * (i+0.5) + py_min;
+ for (j=lc;j<rc;j++) {
+ cxpx = px_dx * (j+0.5) + px_min;
+ cx = inv_mats[0][0]*cxpx + inv_mats[0][1]*cypx + centers[0];
+ cy = inv_mats[1][0]*cxpx + inv_mats[1][1]*cypx + centers[1];
+ cz = inv_mats[2][0]*cxpx + inv_mats[2][1]*cypx + centers[2];
+ if( (fabs(xs[p]-cx)>dxs[p]) ||
+ (fabs(ys[p]-cy)>dys[p]) ||
+ (fabs(zs[p]-cz)>dzs[p])) continue;
+ gridded[j*cols+i] += ds[p];
+ }
+ }
+ }
+
+ // Attatch output buffer to output buffer
+
+ Py_DECREF(x);
+ Py_DECREF(y);
+ Py_DECREF(z);
+ Py_DECREF(px);
+ Py_DECREF(py);
+ Py_DECREF(d);
+ Py_DECREF(dx);
+ Py_DECREF(dy);
+ Py_DECREF(dz);
+ Py_DECREF(center);
+ Py_DECREF(indices);
+ Py_DECREF(inv_mat);
+
+ PyObject *return_value = Py_BuildValue("N", my_array);
+
+ return return_value;
+
+ _fail:
+
+ Py_XDECREF(x);
+ Py_XDECREF(y);
+ Py_XDECREF(z);
+ Py_XDECREF(px);
+ Py_XDECREF(py);
+ Py_XDECREF(d);
+ Py_XDECREF(dx);
+ Py_XDECREF(dy);
+ Py_XDECREF(dz);
+ Py_XDECREF(center);
+ Py_XDECREF(indices);
+ Py_XDECREF(inv_mat);
+
+ return NULL;
+
+}
+
static PyMethodDef __MPLMethods[] = {
{"Pixelize", Py_Pixelize, METH_VARARGS, _pixelizeDocstring},
- //{"CPixelize", Py_CPixelize, METH_VARARGS, NULL},
+ {"CPixelize", Py_CPixelize, METH_VARARGS, NULL},
{NULL, NULL} /* Sentinel */
};
More information about the yt-svn
mailing list