[Yt-svn] yt-commit r744 - in trunk/yt: lagos raven
mturk at wrangler.dreamhost.com
mturk at wrangler.dreamhost.com
Thu Aug 28 11:05:48 PDT 2008
Author: mturk
Date: Thu Aug 28 11:05:46 2008
New Revision: 744
URL: http://yt.spacepope.org/changeset/744
Log:
Some minor changes:
* fixed some old method names, and re-set up projections to autoproject,
but there is probably still a bug there
* changed _MPL.Pixelize to avoid re-copying the arrays, and instead just do
PyArray_GETPTR1 calls
* Fixed calls to min() and max() to work with data that potentially has NaNs
in it.
Bigger changes:
* Added a spectral frequency integrator, which accepts Cloudy data outputs and
makes frequency bins -- **not-well-tested** but seems to work and give good
results compared to older tools.
* FixedResolutionBuffers, which basically respect the appropriate protocol for
data access, but do not derived from EnzoData. Provides back an image from a
2D data source. (Not Cutting Plane, though.)
* AnnuliProfiler, for dealing with FRB's and looking at total output in a
region, as well as a profile of the output in a region.
Added:
trunk/yt/lagos/SpectralIntegrator.py
trunk/yt/raven/FixedResolution.py
Modified:
trunk/yt/lagos/BaseDataTypes.py
trunk/yt/raven/PlotCollection.py
trunk/yt/raven/PlotTypes.py
trunk/yt/raven/_MPL.c
Modified: trunk/yt/lagos/BaseDataTypes.py
==============================================================================
--- trunk/yt/lagos/BaseDataTypes.py (original)
+++ trunk/yt/lagos/BaseDataTypes.py Thu Aug 28 11:05:46 2008
@@ -516,7 +516,7 @@
self.coord += dx * val
else:
raise ValueError(val)
- self.refreshData()
+ self._refresh_data()
def _generate_coords(self):
points = []
@@ -925,7 +925,8 @@
#@time_execution
def get_data(self, fields = None):
- if fields is None: fields = ensure_list(self.fields)
+ if fields is None: fields = self.fields
+ fields = ensure_list(fields)
coord_data = []
field_data = []
dxs = []
Added: trunk/yt/lagos/SpectralIntegrator.py
==============================================================================
--- (empty file)
+++ trunk/yt/lagos/SpectralIntegrator.py Thu Aug 28 11:05:46 2008
@@ -0,0 +1,88 @@
+"""
+Integrator classes to deal with interpolation and integration of input spectral
+bins. Currently only supports Cloudy-style data.
+
+Author: Matthew Turk <matthewturk at gmail.com>
+Affiliation: KIPAC/SLAC/Stanford
+Homepage: http://yt.enzotools.org/
+License:
+ Copyright (C) 2007-2008 Matthew Turk. All Rights Reserved.
+
+ This file is part of yt.
+
+ yt is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+"""
+
+from yt.lagos import *
+
+class SpectralFrequencyIntegrator(object):
+ def __init__(self, table, field_names,
+ bounds, ev_bounds):
+ """
+ From a table, interpolate over field_names to get resultant luminosity.
+ Table must be of the style such that it is ordered by
+ [field_names[0], field_names[1], ev]
+ """
+ self.table = table
+ self.field_names = field_names
+
+ self.bounds = bounds
+ self.ev_bounds = ev_bounds
+ self.ev_vals = na.logspace(ev_bounds[0], ev_bounds[1], table.shape[-1])
+
+ def _get_interpolator(self, ev_min, ev_max):
+ """
+ Integrates from ev_min to ev_max and returns an interpolator.
+ """
+ e_is, e_ie = na.digitize([ev_min, ev_max], self.ev_vals)
+ bin_table = na.trapz(self.table[...,e_is-1:e_ie],
+ 2.41799e17*
+ (self.ev_vals[e_is:e_ie+1]-self.ev_vals[e_is-1:e_is]),
+ axis=-1)
+ bin_table = na.log10(bin_table.clip(1e-80,bin_table.max()))
+ return BilinearFieldInterpolator(
+ bin_table, self.bounds, self.field_names[:],
+ truncate=True)
+
+ def add_frequency_bin_field(self, ev_min, ev_max):
+ """
+ Add a new field to the global fieldInfo dictionary, which is an
+ integrated bin from *ev_min* to *ev_max*.
+
+ Returns the name of the new field.
+ """
+ interp = self._get_interpolator(ev_min, ev_max)
+ name = "XRay_%s_%s" % (ev_min, ev_max)
+ def frequency_bin_field(field, data):
+ dd = {'NumberDensity' : na.log10(data["NumberDensity"]),
+ 'Temperature' : na.log10(data["Temperature"])}
+ return 10**interp(dd)*(data.dx**2.0)
+ def _conv(data):
+ return (data.convert("cm")**3.0)
+ add_field(name, convert_function=_conv, function=frequency_bin_field,
+ projection_conversion="1")
+ return name
+
+def create_table_from_textfiles(pattern, rho_spec, e_spec, T_spec):
+ rho_n_bins, rho_min, rho_max = rho_spec
+ e_n_bins, e_min, e_max = e_spec
+ T_n_bins, T_min, T_max = T_spec
+ # The second one is the fast-varying one
+ rho_is, e_is = na.mgrid[0:rho_n_bins,0:e_n_bins]
+ table = na.zeros((rho_n_bins, T_n_bins, e_n_bins), dtype='float64')
+ mylog.info("Parsing Cloudy files")
+ for i,ri,ei in zip(range(rho_n_bins*e_n_bins), rho_is.ravel(), e_is.ravel()):
+ table[ri,:,ei] = [float(l.split()[-1]) for l in open(pattern%(i+1)) if l[0] != "#"]
+ return table
+
Added: trunk/yt/raven/FixedResolution.py
==============================================================================
--- (empty file)
+++ trunk/yt/raven/FixedResolution.py Thu Aug 28 11:05:46 2008
@@ -0,0 +1,96 @@
+"""
+Fixed resolution buffer support, along with a primitive image analysis tool.
+
+Author: Matthew Turk <matthewturk at gmail.com>
+Affiliation: KIPAC/SLAC/Stanford
+Homepage: http://yt.enzotools.org/
+License:
+ Copyright (C) 2008 Matthew Turk. All Rights Reserved.
+
+ This file is part of yt.
+
+ yt is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+"""
+
+from yt.raven import *
+
+import _MPL
+
+class FixedResolutionBuffer(object):
+ def __init__(self, data_source, bounds, buff_size, antialias = True):
+ """
+ Accepts a 2D data object, such as a Projection or Slice, and implements
+ a protocol for generating a pixelized, fixed-resolution buffer.
+ *bounds* is (px_min,px_max,py_min,py_max), *buff_size* is
+ (width, height), and antialias is a boolean referring to whether or not
+ the buffer should have pixel boundary antialiasing.
+ """
+ self.data_source = data_source
+ self.bounds = bounds
+ self.buff_size = buff_size
+ self.antialias = antialias
+ self.data = {}
+
+ def __getitem__(self, item):
+ if item in self.data: return self.data[item]
+ buff = _MPL.Pixelize(self.data_source['px'],
+ self.data_source['py'],
+ self.data_source['pdx'],
+ self.data_source['pdy'],
+ self.data_source[item],
+ self.buff_size[0], self.buff_size[1],
+ self.bounds, int(self.antialias))
+ self[item] = buff
+ return buff
+
+ def __setitem__(self, item, val):
+ self.data[item] = val
+
+class AnnuliProfiler(object):
+ def __init__(self, fixed_buffer, center, num_bins, min_radius, max_radius):
+ """
+ This is a very simple class, principally used to sum up total values
+ inside annuli in a fixed resolution buffer. It accepts *fixed_buffer*,
+ which should be a FixedResolutionBuffer, *center*, which is in pixel
+ coordinates. *num_bins*, *min_radius* and *max_radius* all refer to
+ the binning properties for the annuli. Note that these are all in
+ pixel values.
+ """
+ self.fixed_buffer = fixed_buffer
+ self.center = center
+ self.num_bins = num_bins
+ self.min_radius = min_radius
+ self.max_radius = max_radius
+ self.bins = na.linspace(min_radius, max_radius, num_bins)
+ self.data = {}
+ self.radii = self._setup_bins()
+
+ def _setup_bins(self):
+ new_x = na.arange(self.fixed_buffer.buff_size[0]) - self.center[0]
+ new_y = na.arange(self.fixed_buffer.buff_size[1]) - self.center[1]
+ radii = na.sqrt((new_x**2.0)[None,:] +
+ (new_y**2.0)[:,None])
+ self.bin_indices = na.digitize(radii.ravel(), self.bins)
+
+ def __getitem__(self, item):
+ if item in self.data: return self.data[item]
+ binned = na.zeros(self.num_bins, dtype='float64')
+ for i in range(self.num_bins):
+ binned[i] = na.sum(self.fixed_buffer[item].ravel()[self.bin_indices==i])
+ self[item] = binned
+ return binned
+
+ def __setitem__(self, item, val):
+ self.data[item] = val
+
Modified: trunk/yt/raven/PlotCollection.py
==============================================================================
--- trunk/yt/raven/PlotCollection.py (original)
+++ trunk/yt/raven/PlotCollection.py Thu Aug 28 11:05:46 2008
@@ -39,7 +39,7 @@
self._run_id = deliverator_id
self.pf = pf
if center == None:
- v,self.c = pf.h.findMax("Density") # @todo: ensure no caching
+ v,self.c = pf.h.find_max("Density") # @todo: ensure no caching
else:
self.c = center
if deliverator_id > 0:
Modified: trunk/yt/raven/PlotTypes.py
==============================================================================
--- trunk/yt/raven/PlotTypes.py (original)
+++ trunk/yt/raven/PlotTypes.py Thu Aug 28 11:05:46 2008
@@ -167,6 +167,8 @@
"""
Set the z boundaries of this plot.
"""
+ # This next call fixes some things, but is slower...
+ #self._redraw_image()
self.norm.autoscale(na.array([zmin,zmax]))
self.image.changed()
if self.colorbar is not None:
@@ -295,25 +297,21 @@
self._axes.clear() # To help out the colorbar
buff = self._get_buff()
mylog.debug("Received buffer of min %s and max %s (%s %s)",
- buff.min(), buff.max(),
+ na.nanmin(buff), na.nanmax(buff),
self[self.axis_names["Z"]].min(),
self[self.axis_names["Z"]].max())
if self.log_field:
bI = na.where(buff > 0)
- newmin = buff[bI].min()
- newmax = buff[bI].max()
+ newmin = na.nanmin(buff[bI])
+ newmax = na.nanmax(buff[bI])
else:
- newmin = buff.min()
- newmax = buff.max()
+ newmin = na.nanmin(buff)
+ newmax = na.nanmax(buff)
if self.do_autoscale:
self.norm.autoscale(na.array((newmin,newmax)))
self.image = \
self._axes.imshow(buff, interpolation='nearest', norm = self.norm,
aspect=1.0, picker=True, origin='lower')
- if self.cmap is not None:
- self.image.set_cmap(self.cmap)
- if self.colorbar is not None:
- self.colorbar.set_cmap(self.cmap)
self._reset_image_parameters()
self._run_callbacks()
@@ -327,6 +325,7 @@
if self.colorbar != None:
self.image.set_norm(self.norm)
self.colorbar.set_norm(self.norm)
+ if self.cmap: self.colorbar.set_cmap(self.cmap)
if self.do_autoscale: _notify(self.image, self.colorbar)
self.autoset_label()
@@ -370,8 +369,8 @@
self._redraw_image()
def autoscale(self):
- zmin = self._axes.images[-1]._A.min()
- zmax = self._axes.images[-1]._A.max()
+ zmin = na.nanmin(self._axes.images[-1]._A)
+ zmax = na.nanmax(self._axes.images[-1]._A)
self.set_zlim(zmin, zmax)
def switch_y(self, *args, **kwargs):
Modified: trunk/yt/raven/_MPL.c
==============================================================================
--- trunk/yt/raven/_MPL.c (original)
+++ trunk/yt/raven/_MPL.c Thu Aug 28 11:05:46 2008
@@ -72,34 +72,34 @@
// Get numeric arrays
PyArrayObject *x = (PyArrayObject *) PyArray_FromAny(xp,
- PyArray_DescrFromType(NPY_FLOAT64), 1, 1, NPY_C_CONTIGUOUS, NULL);
+ PyArray_DescrFromType(NPY_FLOAT64), 1, 1, 0, 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);
+ PyArray_DescrFromType(NPY_FLOAT64), 1, 1, 0, NULL);
if ((y == NULL) || (PyArray_SIZE(y) != PyArray_SIZE(x))) {
PyErr_Format( _pixelizeError, "y 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);
+ PyArray_DescrFromType(NPY_FLOAT64), 1, 1, 0, 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);
+ PyArray_DescrFromType(NPY_FLOAT64), 1, 1, 0, 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);
+ PyArray_DescrFromType(NPY_FLOAT64), 1, 1, 0, NULL);
if ((dy == NULL) || (PyArray_SIZE(dy) != PyArray_SIZE(x))) {
PyErr_Format( _pixelizeError, "dy is of incorrect type (wanted 1D float)");
goto _fail;
@@ -115,40 +115,43 @@
int i, j, p;
double lc, lr, rc, rr;
double lypx, rypx, lxpx, rxpx, overlap1, overlap2;
- npy_float64 *xs = (npy_float64 *) PyArray_GETPTR1(x, 0);
- npy_float64 *ys = (npy_float64 *) PyArray_GETPTR1(y, 0);
- npy_float64 *dxs = (npy_float64 *) PyArray_GETPTR1(dx, 0);
- npy_float64 *dys = (npy_float64 *) PyArray_GETPTR1(dy, 0);
- npy_float64 *ds = (npy_float64 *) PyArray_GETPTR1(d, 0); // We check this above
+ npy_float64 xsp, ysp, dxsp, dysp, dsp;
npy_intp dims[] = {rows, cols};
PyArrayObject *my_array =
(PyArrayObject *) PyArray_SimpleNewFromDescr(2, dims,
PyArray_DescrFromType(NPY_FLOAT64));
- npy_float64 *gridded = (npy_float64 *) my_array->data;
+ //npy_float64 *gridded = (npy_float64 *) my_array->data;
- for(p=0;p<cols*rows;p++)gridded[p]=0.0;
+ for(i=0;i<rows;i++)for(j=0;j<cols;j++)
+ *(npy_float64*) PyArray_GETPTR2(my_array, j, i) = 0.0;
for(p=0;p<nx;p++)
{
- if(((xs[p]+dxs[p]<x_min) ||
- (xs[p]-dxs[p]>x_max)) ||
- ((ys[p]+dys[p]<y_min) ||
- (ys[p]-dys[p]>y_max))) continue;
- lc = max(((xs[p]-dxs[p]-x_min)/px_dx),0);
- lr = max(((ys[p]-dys[p]-y_min)/px_dy),0);
- rc = min(((xs[p]+dxs[p]-x_min)/px_dx), rows);
- rr = min(((ys[p]+dys[p]-y_min)/px_dy), cols);
+ xsp = *((npy_float64 *)PyArray_GETPTR1(x, p));
+ ysp = *((npy_float64 *)PyArray_GETPTR1(y, p));
+ dxsp = *((npy_float64 *)PyArray_GETPTR1(dx, p));
+ dysp = *((npy_float64 *)PyArray_GETPTR1(dy, p));
+ dsp = *((npy_float64 *)PyArray_GETPTR1(d, p));
+ if(((xsp+dxsp<x_min) ||
+ (xsp-dxsp>x_max)) ||
+ ((ysp+dysp<y_min) ||
+ (ysp-dysp>y_max))) continue;
+ lc = max(((xsp-dxsp-x_min)/px_dx),0);
+ lr = max(((ysp-dysp-y_min)/px_dy),0);
+ rc = min(((xsp+dxsp-x_min)/px_dx), rows);
+ rr = min(((ysp+dysp-y_min)/px_dy), cols);
for (i=lr;i<rr;i++) {
lypx = px_dy * i + y_min;
rypx = px_dy * (i+1) + y_min;
- overlap2 = ((min(rypx, ys[p]+dys[p]) - max(lypx, (ys[p]-dys[p])))/px_dy);
+ overlap2 = ((min(rypx, ysp+dysp) - max(lypx, (ysp-dysp)))/px_dy);
for (j=lc;j<rc;j++) {
lxpx = px_dx * j + x_min;
rxpx = px_dx * (j+1) + x_min;
- overlap1 = ((min(rxpx, xs[p]+dxs[p]) - max(lxpx, (xs[p]-dxs[p])))/px_dx);
+ overlap1 = ((min(rxpx, xsp+dxsp) - max(lxpx, (xsp-dxsp)))/px_dx);
if (overlap1 < 0.0 || overlap2 < 0.0) continue;
- if (antialias == 1) gridded[j*cols+i] += (ds[p]*overlap1)*overlap2;
- else gridded[j*cols+i] = ds[p];
+ if (antialias == 1)
+ *(npy_float64*) PyArray_GETPTR2(my_array, j, i) += (dsp*overlap1)*overlap2;
+ else *(npy_float64*) PyArray_GETPTR2(my_array, j, i) = dsp;
}
}
}
@@ -200,61 +203,61 @@
// Get numeric arrays
PyArrayObject *x = (PyArrayObject *) PyArray_FromAny(xp,
- PyArray_DescrFromType(NPY_FLOAT64), 1, 1, NPY_C_CONTIGUOUS, NULL);
+ PyArray_DescrFromType(NPY_FLOAT64), 1, 1, 0, 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);
+ PyArray_DescrFromType(NPY_FLOAT64), 1, 1, 0, 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);
+ PyArray_DescrFromType(NPY_FLOAT64), 1, 1, 0, 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);
+ PyArray_DescrFromType(NPY_FLOAT64), 1, 1, 0, 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);
+ PyArray_DescrFromType(NPY_FLOAT64), 1, 1, 0, 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);
+ PyArray_DescrFromType(NPY_FLOAT64), 1, 1, 0, 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);
+ PyArray_DescrFromType(NPY_FLOAT64), 1, 1, 0, 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);
+ PyArray_DescrFromType(NPY_FLOAT64), 1, 1, 0, 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);
+ PyArray_DescrFromType(NPY_FLOAT64), 1, 1, 0, NULL);
if ((dz == NULL) || (PyArray_SIZE(dz) != PyArray_SIZE(x))) {
PyErr_Format( _pixelizeError, "dz is of incorrect type (wanted 1D float)");
goto _fail;
@@ -266,13 +269,13 @@
goto _fail;
}
PyArrayObject *inv_mat = (PyArrayObject *) PyArray_FromAny(inv_matp,
- PyArray_DescrFromType(NPY_FLOAT64), 2, 2, NPY_C_CONTIGUOUS, NULL);
+ PyArray_DescrFromType(NPY_FLOAT64), 2, 2, 0, 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);
+ PyArray_DescrFromType(NPY_INT64), 1, 1, 0, NULL);
if ((indices == NULL) || (PyArray_SIZE(indices) != PyArray_SIZE(dx))) {
PyErr_Format( _pixelizeError, "indices must be same length as dx");
goto _fail;
@@ -287,17 +290,8 @@
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 xsp, ysp, zsp, pxsp, pysp, dxsp, dysp, dzsp, dsp;
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 =
@@ -315,17 +309,26 @@
for(p=0;p<cols*rows;p++)gridded[p]=mask[p]=0.0;
for(pp=0; pp<nx; pp++)
{
- p = indicess[pp];
+ p = *((npy_int64 *) PyArray_GETPTR1(indices, pp));
+ npy_float64 xsp = *((npy_float64 *) PyArray_GETPTR1(x, p));
+ npy_float64 ysp = *((npy_float64 *) PyArray_GETPTR1(y, p));
+ npy_float64 zsp = *((npy_float64 *) PyArray_GETPTR1(z, p));
+ npy_float64 pxsp = *((npy_float64 *) PyArray_GETPTR1(px, p));
+ npy_float64 pysp = *((npy_float64 *) PyArray_GETPTR1(py, p));
+ npy_float64 dxsp = *((npy_float64 *) PyArray_GETPTR1(dx, p));
+ npy_float64 dysp = *((npy_float64 *) PyArray_GETPTR1(dy, p));
+ npy_float64 dzsp = *((npy_float64 *) PyArray_GETPTR1(dz, p));
+ npy_float64 dsp = *((npy_float64 *) PyArray_GETPTR1(d, p)); // We check this above
// 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);
+ md = 2.0*sqrtl(dxsp*dxsp + dysp*dysp + dzsp*dzsp);
+ if(((pxsp+md<px_min) ||
+ (pxsp-md>px_max)) ||
+ ((pysp+md<py_min) ||
+ (pysp-md>py_max))) continue;
+ lc = max(floorl((pxsp-md-px_min)/px_dx),0);
+ lr = max(floorl((pysp-md-py_min)/px_dy),0);
+ rc = min(ceill((pxsp+md-px_min)/px_dx),rows);
+ rr = min(ceill((pysp+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++) {
@@ -333,11 +336,11 @@
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)*0.95>dxs[p]) ||
- (fabs(ys[p]-cy)*0.95>dys[p]) ||
- (fabs(zs[p]-cz)*0.95>dzs[p])) continue;
+ if( (fabs(xsp-cx)*0.95>dxsp) ||
+ (fabs(ysp-cy)*0.95>dysp) ||
+ (fabs(zsp-cz)*0.95>dzsp)) continue;
mask[j*cols+i] += 1;
- gridded[j*cols+i] += ds[p];
+ gridded[j*cols+i] += dsp;
}
}
}
More information about the yt-svn
mailing list