[Yt-svn] yt-commit r1674 - in trunk/yt: extensions/image_panner raven

mturk at wrangler.dreamhost.com mturk at wrangler.dreamhost.com
Sun Mar 28 09:20:53 PDT 2010


Author: mturk
Date: Sun Mar 28 09:20:52 2010
New Revision: 1674
URL: http://yt.enzotools.org/changeset/1674

Log:
Backport from hg:

 * Optimizations to the pixelization step, which should increase performance by about 40% at most for some pixelizations.
 * Adding an image panner, to be better documented, but which provides us with
the ability to control tile display walls, use Chaco for a "google maps" style
interface to plots, and simplify our overall interaction with variable mesh
plots as a whole.

I'm pretty excited about this image panner, because it's a very cool and fun
way to interact with 2D AMR data.  And it's able to operate remotely, so the
client can run on a local machine (where Chaco may be easier to install, for
intance) getting remote images passed back.  It can also handle "windowing", so
that the pixelization jobs can be distributed across machines -- for instance,
for tile displays!



Added:
   trunk/yt/extensions/image_panner/
   trunk/yt/extensions/image_panner/__init__.py
   trunk/yt/extensions/image_panner/pan_and_scan_widget.py
   trunk/yt/extensions/image_panner/vm_panner.py
Modified:
   trunk/yt/raven/PlotTypes.py
   trunk/yt/raven/_MPL.c

Added: trunk/yt/extensions/image_panner/__init__.py
==============================================================================
--- (empty file)
+++ trunk/yt/extensions/image_panner/__init__.py	Sun Mar 28 09:20:52 2010
@@ -0,0 +1,26 @@
+"""
+Author: Matthew Turk <matthewturk at gmail.com>
+Affiliation: KIPAC/SLAC/Stanford
+Homepage: http://yt.enzotools.org/
+License:
+  Copyright (C) 2010 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 vm_panner import VariableMeshPanner, WindowedVariableMeshPanner, \
+                MultipleWindowVariableMeshPanner, ImageSaver, \
+                PanningCeleritasStreamer

Added: trunk/yt/extensions/image_panner/pan_and_scan_widget.py
==============================================================================
--- (empty file)
+++ trunk/yt/extensions/image_panner/pan_and_scan_widget.py	Sun Mar 28 09:20:52 2010
@@ -0,0 +1,200 @@
+"""
+Author: Matthew Turk <matthewturk at gmail.com>
+Affiliation: KIPAC/SLAC/Stanford
+Homepage: http://yt.enzotools.org/
+License:
+  Copyright (C) 2010 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/>.
+"""
+
+# Major library imports
+from vm_panner import VariableMeshPanner
+from numpy import linspace, meshgrid, pi, sin, mgrid, zeros
+
+# Enthought library imports
+from enthought.enable.api import Component, ComponentEditor, Window
+from enthought.traits.api import HasTraits, Instance, Button, Any, Callable, \
+        on_trait_change, Bool
+from enthought.traits.ui.api import Item, Group, View
+
+# Chaco imports
+from enthought.chaco.api import ArrayPlotData, jet, Plot, HPlotContainer, \
+        ColorBar, DataRange1D, DataRange2D, LinearMapper, ImageData, \
+        CMapImagePlot
+from enthought.chaco.tools.api import PanTool, ZoomTool, RangeSelection, \
+        RangeSelectionOverlay, RangeSelection
+from zoom_overlay import ZoomOverlay
+from enthought.chaco.tools.image_inspector_tool import ImageInspectorTool, \
+     ImageInspectorOverlay
+
+if not hasattr(DataRange2D, "_subranges_updated"):
+    print "You'll need to add _subranges updated to enthought/chaco/data_range_2d.py"
+    print 'Add this at the correct indentation level:'
+    print
+    print '    @on_trait_change("_xrange.updated,_yrange.updated")'
+    print '    def _subranges_updated(self):'
+    print '        self.updated = True'
+    print
+    raise RuntimeError
+
+class FunctionImageData(ImageData):
+    # The function to call with the low and high values of the range.
+    # It should return an array of values.
+    func = Callable
+
+    # A reference to a datarange
+    data_range = Instance(DataRange2D)
+
+    def __init__(self, **kw):
+        # Explicitly call the AbstractDataSource constructor because
+        # the ArrayDataSource ctor wants a data array
+        ImageData.__init__(self, **kw)
+        self.recalculate()
+
+    @on_trait_change('data_range.updated')
+    def recalculate(self):
+        if self.func is not None and self.data_range is not None:
+            newarray = self.func(self.data_range.low, self.data_range.high)
+            ImageData.set_data(self, newarray)
+        else:
+            self._data = zeros((512,512),dtype=float)
+
+    def set_data(self, *args, **kw):
+        raise RuntimeError("Cannot set numerical data on a FunctionDataSource")
+
+    def set_mask(self, mask):
+        raise NotImplementedError
+
+    def remove_mask(self):
+        raise NotImplementedError
+
+class ImagePixelizerHelper(object):
+    index = None
+    def __init__(self, panner):
+        self.panner = panner
+
+    def __call__(self, low, high):
+        b = self.panner.set_low_high(low, high)
+        if self.index is not None:
+            num_x_ticks = b.shape[0] + 1
+            num_y_ticks = b.shape[1] + 1
+            xs = mgrid[low[0]:high[0]:num_x_ticks*1j]
+            ys = mgrid[low[1]:high[1]:num_y_ticks*1j]
+            self.index.set_data( xs, ys )
+        return b
+
+class VMImagePlot(HasTraits):
+    plot = Instance(Plot)
+    fid = Instance(FunctionImageData)
+    img_plot = Instance(CMapImagePlot)
+    panner = Instance(VariableMeshPanner)
+    helper = Instance(ImagePixelizerHelper)
+
+    def _plot_default(self):
+        pd = ArrayPlotData()
+        plot = Plot(pd)
+        self.fid._data = self.panner.buffer
+
+        pd.set_data("imagedata", self.fid)
+
+        img_plot = plot.img_plot("imagedata", colormap=jet,
+                                 interpolation='nearest',
+                                 xbounds=(0.0, 1.0),
+                                 ybounds=(0.0, 1.0))[0]
+        self.fid.data_range = plot.range2d
+        self.helper.index = img_plot.index
+        self.img_plot = img_plot
+        return plot
+
+    def _fid_default(self):
+        return FunctionImageData(func = self.helper)
+
+    def _helper_default(self):
+        return ImagePixelizerHelper(self.panner)
+
+class VariableMeshPannerView(HasTraits):
+
+    plot = Instance(Plot)
+    spawn_zoom = Button
+    vm_plot = Instance(VMImagePlot)
+    use_tools = Bool(True)
+    
+    traits_view = View(
+                    Group(
+                        Item('container', editor=ComponentEditor(size=(512,512)), 
+                             show_label=False),
+                        Item('spawn_zoom', show_label=False),
+                        orientation = "vertical"),
+                    width = 800, height=800,
+                    resizable=True, title="Pan and Scan",
+                    )
+
+    def _vm_plot_default(self):
+        return VMImagePlot(panner=self.panner)
+    
+    def __init__(self, **kwargs):
+        super(VariableMeshPannerView, self).__init__(**kwargs)
+        # Create the plot
+
+        plot = self.vm_plot.plot
+        img_plot = self.vm_plot.img_plot
+
+        if self.use_tools:
+            plot.tools.append(PanTool(img_plot))
+            zoom = ZoomTool(component=img_plot, tool_mode="box", always_on=False)
+            plot.overlays.append(zoom)
+            imgtool = ImageInspectorTool(img_plot)
+            img_plot.tools.append(imgtool)
+            overlay = ImageInspectorOverlay(component=img_plot, image_inspector=imgtool,
+                                            bgcolor="white", border_visible=True)
+            img_plot.overlays.append(overlay)
+
+
+        image_value_range = DataRange1D(self.vm_plot.fid)
+        cbar_index_mapper = LinearMapper(range=image_value_range)
+        self.colorbar = ColorBar(index_mapper=cbar_index_mapper,
+                                 plot=img_plot,
+                                 padding_right=40,
+                                 resizable='v',
+                                 width=30)
+
+        self.colorbar.tools.append(
+            PanTool(self.colorbar, constrain_direction="y", constrain=True))
+        zoom_overlay = ZoomTool(self.colorbar, axis="index", tool_mode="range",
+                                always_on=True, drag_button="right")
+        self.colorbar.overlays.append(zoom_overlay)
+
+        # create a range selection for the colorbar
+        range_selection = RangeSelection(component=self.colorbar)
+        self.colorbar.tools.append(range_selection)
+        self.colorbar.overlays.append(
+                RangeSelectionOverlay(component=self.colorbar,
+                                      border_color="white",
+                                      alpha=0.8, fill_color="lightgray"))
+
+        # we also want to the range selection to inform the cmap plot of
+        # the selection, so set that up as well
+        range_selection.listeners.append(img_plot)
+
+        self.container = HPlotContainer(padding=30)
+        self.container.add(self.colorbar)
+        self.container.add(self.vm_plot.plot)
+
+    def _spawn_zoom_fired(self):
+        np = self.panner.source.pf.h.image_panner(
+                self.panner.source, self.panner.size, self.panner.field)
+        new_window = VariableMeshPannerView(panner = np)

Added: trunk/yt/extensions/image_panner/vm_panner.py
==============================================================================
--- (empty file)
+++ trunk/yt/extensions/image_panner/vm_panner.py	Sun Mar 28 09:20:52 2010
@@ -0,0 +1,329 @@
+"""
+Author: Matthew Turk <matthewturk at gmail.com>
+Affiliation: KIPAC/SLAC/Stanford
+Homepage: http://yt.enzotools.org/
+License:
+  Copyright (C) 2010 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/>.
+"""
+
+import numpy as na
+import types, os
+from yt.raven import FixedResolutionBuffer, ObliqueFixedResolutionBuffer
+from yt.lagos import data_object_registry, AMRProjBase, AMRSliceBase, \
+                     x_dict, y_dict
+
+class VariableMeshPanner(object):
+    _buffer = None
+
+    def __init__(self, source, size, field, callback = None,
+                 viewport_callback = None):
+        if not isinstance(source, (AMRProjBase, AMRSliceBase)):
+            raise RuntimeError
+        if callback is None:
+            callback = lambda a: None
+        self.callback = callback
+        if viewport_callback is None:
+            viewport_callback = lambda a, b: None
+        self.viewport_callback = viewport_callback
+        self.size = size
+        self.source = source
+        self.field = field
+        self.xlim, self.ylim = self.bounds
+
+    def _run_callbacks(self):
+        self.callback(self.buffer)
+        self.viewport_callback(self.xlim, self.ylim)
+
+    @property
+    def bounds(self):
+        DLE, DRE = self.pf["DomainLeftEdge"], self.pf["DomainRightEdge"]
+        ax = self.source.axis
+        xax, yax = x_dict[ax], y_dict[ax]
+        xbounds = DLE[xax], DRE[xax]
+        ybounds = DLE[yax], DRE[yax]
+        return (xbounds, ybounds)
+
+    @property
+    def width(self):
+        Wx = self.xlim[1] - self.xlim[0]
+        Wy = self.ylim[1] - self.ylim[0]
+        return (Wx, Wy)
+
+    def zoom(self, factor):
+        Wx, Wy = self.width
+        centerx = self.xlim[0] + Wx*0.5
+        centery = self.ylim[0] + Wy*0.5
+        nWx, nWy = Wx/factor, Wy/factor
+        self.xlim = (centerx - nWx*0.5, centerx + nWx*0.5)
+        self.ylim = (centery - nWy*0.5, centery + nWy*0.5)
+        self._run_callbacks()
+
+    def pan(self, deltas):
+        self.xlim = (self.xlim[0] + deltas[0], self.xlim[1] + deltas[0])
+        self.ylim = (self.ylim[0] + deltas[1], self.ylim[1] + deltas[1])
+        self._run_callbacks()
+
+    def pan_x(self, delta):
+        self.pan( (delta, 0.0) )
+
+    def pan_y(self, delta):
+        self.pan( (0.0, delta) )
+
+    def pan_rel(self, deltas):
+        Wx, Wy = self.width
+        px = deltas[0] * Wx
+        py = deltas[1] * Wy
+        self.pan( (px, py) )
+
+    def pan_rel_x(self, delta):
+        self.pan_rel( (delta, 0.0) )
+
+    def pan_rel_y(self, delta):
+        self.pan_rel( (0.0, delta) )
+
+    @property
+    def buffer(self):
+        my_bounds = self.xlim + self.ylim # Tuples concatenate
+        # If the buffer is None, we just kill it regardless.
+        if self._buffer is None:
+            self._regenerate_buffer()
+        elif self._buffer.bounds != my_bounds:
+            self._regenerate_buffer()
+        return self._buffer[self.field]
+
+    def _regenerate_buffer(self):
+        new_buffer = FixedResolutionBuffer(
+            self.source, self.xlim + self.ylim,
+            self.size)
+        self._buffer = new_buffer
+
+    def set_low_high(self, low, high):
+        self.xlim = (low[0], high[0])
+        self.ylim = (low[1], high[1])
+        return na.log10(self.buffer)
+
+    def set_limits(self, xlim, ylim):
+        self.xlim = xlim
+        self.ylim = ylim
+        self._run_callbacks()
+
+data_object_registry["image_panner"] = VariableMeshPanner
+
+class WindowedVariableMeshPanner(VariableMeshPanner):
+
+    def __init__(self, source, full_size, my_size, start_indices,
+                 field, callback = None):
+        self.my_size = my_size
+        self.start_indices = start_indices
+        VariableMeshPanner.__init__(self, source, full_size, field, callback)
+
+    def _regenerate_buffer(self):
+        dx = (self.xlim[1] - self.xlim[0])/self.size[0]
+        dy = (self.ylim[1] - self.ylim[0])/self.size[1]
+        my_lim = (self.xlim[0] + dx*self.start_indices[0],
+                  self.xlim[0] + dx*(self.start_indices[0] + self.my_size[0]),
+                  self.ylim[0] + dx*self.start_indices[1],
+                  self.ylim[0] + dx*(self.start_indices[1] + self.my_size[1]))
+        new_buffer = FixedResolutionBuffer(self.source, my_lim, self.my_size)
+        self._buffer = new_buffer
+
+data_object_registry["windowed_image_panner"] = WindowedVariableMeshPanner
+
+class MultipleWindowVariableMeshPanner(object):
+    def __init__(self, windows):
+        self.windows = windows
+
+    def zoom(self, factor):
+        for w in self.windows: w.zoom(factor)
+
+    def pan(self, deltas):
+        for w in self.windows: w.pan(factor)
+
+    def pan_x(self, delta):
+        for w in self.windows: w.pan_x(delta)
+
+    def pan_y(self, delta):
+        for w in self.windows: w.pan_y(delta)
+
+    def pan_rel(self, deltas):
+        for w in self.windows: w.pan_rel(deltas)
+
+    def pan_rel_x(self, delta):
+        for w in self.windows: w.pan_rel_x(delta)
+
+    def pan_rel_y(self, delta):
+        for w in self.windows: w.pan_rel_y(delta)
+
+    def set_limits(self, xlim, ylim):
+        for w in self.windows: w.set_limits(xlim, ylim)
+
+class RemoteWindowedVariableMeshController(MultipleWindowVariableMeshPanner):
+    def __init__(self, source, mec = None):
+        if mec is None:
+            from IPython.kernel.client import get_multiengine_client
+            mec = get_multiengine_client()
+        self.mec = mec
+        self.mec.execute("import yt.extensions.image_panner")
+        self._var_name = "_image_panner_%s" % (id(self))
+        self._pf_name = "_pf_%s" % (id(self))
+        self._source_name = "_source_%s" % (id(self))
+        self.source = source
+        self.mec.execute("from yt.mods import *")
+        self.mec.execute("from yt.funcs import iterable")
+        self.mec.push({self._pf_name: self.pf})
+        self.mec.execute("%s.h" % self._pf_name)
+        self.mec.push({self._source_name: self.source})
+        # Now, because the double pickling tosses a PF hash reference inside
+        # the unpickled object, we work around it a little
+        self.mec.execute("while iterable(%s): %s = %s[1]" % (
+            self._source_name, self._source_name, self._source_name))
+        self.windows = []
+
+    def add_window(self, *args, **kwargs):
+        engine_id = len(self.windows)
+        an = "_args_%s" % id(self)
+        kn = "_kwargs_%s" % id(self)
+        if 'callback' not in kwargs:
+            kwargs['callback'] = ImageSaver(engine_id)
+        self.mec.push({an: args, kn: kwargs}, engine_id)
+        exec_string = "%s = %s.h.windowed_image_panner(%s, *%s, **%s)" % (
+            self._var_name, self._pf_name, self._source_name, an, kn)
+        self.mec.execute(exec_string, engine_id)
+        self.windows.append(WindowedVariableMeshPannerProxy(
+            self.mec, engine_id, self._var_name, id(self)))
+
+data_object_registry["remote_image_panner"] = RemoteWindowedVariableMeshController
+
+_wrapped_methods = ["zoom", "pan", "pan_x", "pan_y", "pan_rel",
+                     "pan_rel_x", "pan_rel_y", "set_limits"]
+
+class WindowedVariableMeshPannerProxy(object):
+    class __metaclass__(type):
+        def __new__(cls, name, b, d):
+            # We add on a bunch of proxy functions
+            def return_proxy(fname):
+                def func(self, *args, **kwargs):
+                    vn = "_ret_%s" % self._cid
+                    an = "_args_%s" % self._cid
+                    kn = "_kwargs_%s" % self._cid
+                    self.mec.push({an: args, kn: kwargs}, self.engine_id)
+                    exec_string = "%s = %s.%s(*%s, **%s)" % (
+                        vn, self._var_name, fname, an, kn)
+                    print "Executing %s on %s" % (exec_string, self.engine_id)
+                    self.mec.execute(exec_string, self.engine_id)
+                    return self.mec.pull(vn, self.engine_id)
+                return func
+            new_dict = {}
+            new_dict.update(d)
+            for f in _wrapped_methods:
+                print "Constructing proxy for", f
+                new_dict[f] = return_proxy(f)
+            return type.__new__(cls, name, b, new_dict)
+
+    def __init__(self, mec, engine_id, var_name, cid):
+        # mec here is, optionally, an instance of MultiEngineClient
+        self._var_name = var_name
+        self._cid = cid
+        self.engine_id = engine_id
+        self.mec = mec
+
+    @property
+    def bounds(self):
+        vn = "_ret_%s" % self._cid
+        self.mec.execute("%s = %s.bounds" % (vn, self._var_name),
+                         self.engine_id)
+        return self.mec.pull(vn, self.engine_id)
+
+    @property
+    def width(self):
+        vn = "_ret_%s" % self._cid
+        self.mec.execute("%s = %s.width" % (vn, self._var_name),
+                         self.engine_id)
+        return self.mec.pull(vn, self.engine_id)
+
+    @property
+    def buffer(self):
+        vn = "_ret_%s" % self._cid
+        self.mec.execute("%s = %s.buffer" % (vn, self._var_name),
+                         self.engine_id)
+        return self.mec.pull(vn, self.engine_id)
+
+    def _regenerate_buffer(self):
+        return
+
+    def _run_callback(self):
+        self.mec.execute("%s._regenerate_buffer()" % self._var_name,
+                         self.engine_id)
+        self.mec.execute("%s.callback(%s.buffer)" % (
+            self._var_name, self._var_name), self.engine_id)
+
+class ImageSaver(object):
+    def __init__(self, tile_id):
+        self.tile_id = tile_id
+
+        import matplotlib;matplotlib.use("Agg");import pylab
+        self.pylab = pylab
+        self.pylab.clf()
+        fig = pylab.gcf()
+        fig.subplots_adjust(left=0.0, bottom=0.0,
+                            right=1.0, top=1.0,
+                            wspace=0.0, hspace=0.0)
+        fig.set_dpi(100.0)
+        fig.set_size_inches((5.12, 5.12))
+
+    def __call__(self, val):
+        self.pylab.clf()
+        self.pylab.imshow(na.log10(val), interpolation='nearest')
+        self.pylab.savefig("wimage_%03i.png" % self.tile_id)
+
+class PanningCeleritasStreamer(object):
+    _initialized = False
+    def __init__(self, tile_id, cmap = "algae", port = 9988,
+                 zlim = (0.0, 1.0), take_log = True):
+        self.tile_id = tile_id
+        self._port = port
+        self.cmap = cmap
+        self.zlim = zlim
+        self.take_log = True
+
+    def initialize(self, shape):
+        if isinstance(self.cmap, types.StringTypes):
+            import matplotlib.cm
+            self.cmap = matplotlib.cm.get_cmap(self.cmap)
+
+        import celeritas_streamer
+        self.cs = celeritas_streamer.CeleritasStream()
+        #print "Setting shape: %s and port: %s in %s" % (
+        #    shape, self._port, os.getpid())
+        self.cs.setSize(*shape)
+        self.cs.setLocalPort(self._port)
+        self.cs.initialize()
+        self._initialized = True
+
+    def __call__(self, val):
+        if not self._initialized: self.initialize(val.shape)
+        if self.take_log:
+            vv = na.log10(val)
+        else:
+            vv = val.copy()
+        na.subtract(vv, self.zlim[0], vv)
+        na.divide(vv, (self.zlim[1]-self.zlim[0]), vv)
+        new_buf = self.cmap(vv)[:,:,:3]
+        na.multiply(new_buf, 255.0, new_buf)
+        new_buf = new_buf.astype('uint8')
+        self.cs.readFromRGBMemAndSend(new_buf)

Modified: trunk/yt/raven/PlotTypes.py
==============================================================================
--- trunk/yt/raven/PlotTypes.py	(original)
+++ trunk/yt/raven/PlotTypes.py	Sun Mar 28 09:20:52 2010
@@ -280,6 +280,7 @@
             self.ymin = DLE[yax] - DD[yax]
             self.ymax = DRE[yax] + DD[yax]
             self._period = (DD[xax], DD[yax])
+            self._edges = ( (DLE[xax], DRE[xax]), (DLE[yax], DRE[yax]) )
         else:
             # Not quite sure how to deal with this, particularly
             # in the Orion case.  Cutting planes are tricky.
@@ -296,6 +297,7 @@
 class VMPlot(RavenPlot):
     _antialias = True
     _period = (0.0, 0.0)
+    _edges = True
     def __init__(self, data, field, figure = None, axes = None,
                  use_colorbar = True, size=None, periodic = False):
         fields = ['X', 'Y', field, 'X width', 'Y width']
@@ -360,13 +362,20 @@
         # 'px' == pixel x, or x in the plane of the slice
         # 'x' == actual x
         aa = int(self._antialias)
+        if self._edges is True or \
+           x0 < self._edges[0][0] or x1 > self._edges[0][1] or \
+           y0 < self._edges[1][0] or y1 > self._edges[1][1]:
+            check_period = 1
+        else:
+            check_period = 0
         buff = _MPL.Pixelize(self.data['px'],
                             self.data['py'],
                             self.data['pdx'],
                             self.data['pdy'],
                             self[self.axis_names["Z"]],
                             int(height), int(width),
-                            (x0, x1, y0, y1),aa,self._period).transpose()
+                            (x0, x1, y0, y1),aa,self._period,
+                            check_period).transpose()
         return buff
 
     def _redraw_image(self, *args):

Modified: trunk/yt/raven/_MPL.c
==============================================================================
--- trunk/yt/raven/_MPL.c	(original)
+++ trunk/yt/raven/_MPL.c	Sun Mar 28 09:20:52 2010
@@ -58,17 +58,20 @@
   double x_min, x_max, y_min, y_max;
   double period_x, period_y;
   period_x = period_y = 0;
+  int check_period = 1;
 
-  if (!PyArg_ParseTuple(args, "OOOOOII(dddd)|i(dd)",
+  if (!PyArg_ParseTuple(args, "OOOOOII(dddd)|i(dd)i",
       &xp, &yp, &dxp, &dyp, &dp, &cols, &rows,
       &x_min, &x_max, &y_min, &y_max,
-      &antialias, &period_x, &period_y))
+      &antialias, &period_x, &period_y, &check_period))
       return PyErr_Format(_pixelizeError, "Pixelize: Invalid Parameters.");
 
   double width = x_max - x_min;
   double height = y_max - y_min;
   double px_dx = width / ((double) rows);
   double px_dy = height / ((double) cols);
+  double ipx_dx = 1.0 / px_dx;
+  double ipx_dy = 1.0 / px_dy;
 
   // Check we have something to output to
   if (rows == 0 || cols ==0)
@@ -121,6 +124,9 @@
   double lypx, rypx, lxpx, rxpx, overlap1, overlap2;
   npy_float64 oxsp, oysp, xsp, ysp, dxsp, dysp, dsp;
   int xiter[2], yiter[2];
+  double xiterv[2], yiterv[2];
+
+  
 
   npy_intp dims[] = {rows, cols};
   PyArrayObject *my_array =
@@ -129,6 +135,7 @@
   //npy_float64 *gridded = (npy_float64 *) my_array->data;
 
   xiter[0] = yiter[0] = 0;
+  xiterv[0] = yiterv[0] = 0.0;
 
   for(i=0;i<rows;i++)for(j=0;j<cols;j++)
       *(npy_float64*) PyArray_GETPTR2(my_array, i, j) = 0.0;
@@ -142,33 +149,35 @@
     dysp = *((npy_float64 *)PyArray_GETPTR1(dy, p));
     dsp = *((npy_float64 *)PyArray_GETPTR1(d, p));
     xiter[1] = yiter[1] = 999;
-    if (oxsp < x_min) {xiter[1] = +1;}
-    else if (oxsp > x_max) {xiter[1] = -1;}
-    if (oysp < y_min) {yiter[1] = +1;}
-    else if (oysp > y_max) {yiter[1] = -1;}
+    if(check_period == 1) {
+      if (oxsp < x_min) {xiter[1] = +1; xiterv[1] = period_x;}
+      else if (oxsp > x_max) {xiter[1] = -1; xiterv[1] = -period_x;}
+      if (oysp < y_min) {yiter[1] = +1; yiterv[1] = period_y;}
+      else if (oysp > y_max) {yiter[1] = -1; yiterv[1] = -period_y;}
+    }
     for(xi = 0; xi < 2; xi++) {
       if(xiter[xi] == 999)continue;
-      xsp = oxsp + xiter[xi]*period_x;
+      xsp = oxsp + xiterv[xi];
+      if((xsp+dxsp<x_min) || (xsp-dxsp>x_max)) continue;
       for(yi = 0; yi < 2; yi++) {
         if(yiter[yi] == 999)continue;
-        ysp = oysp + yiter[yi]*period_y;
-        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);
+        ysp = oysp + yiterv[yi];
+        if((ysp+dysp<y_min) || (ysp-dysp>y_max)) continue;
+        lc = max(((xsp-dxsp-x_min)*ipx_dx),0);
+        lr = max(((ysp-dysp-y_min)*ipx_dy),0);
+        rc = min(((xsp+dxsp-x_min)*ipx_dx), rows);
+        rr = min(((ysp+dysp-y_min)*ipx_dy), cols);
         for (i=lr;i<rr;i++) {
           lypx = px_dy * i + y_min;
           rypx = px_dy * (i+1) + y_min;
-          overlap2 = ((min(rypx, ysp+dysp) - max(lypx, (ysp-dysp)))/px_dy);
+          overlap2 = dsp*((min(rypx, ysp+dysp) - max(lypx, (ysp-dysp)))*ipx_dy);
           for (j=lc;j<rc;j++) {
             lxpx = px_dx * j + x_min;
             rxpx = px_dx * (j+1) + x_min;
-            overlap1 = ((min(rxpx, xsp+dxsp) - max(lxpx, (xsp-dxsp)))/px_dx);
+            overlap1 = ((min(rxpx, xsp+dxsp) - max(lxpx, (xsp-dxsp)))*ipx_dx);
             if (overlap1 < 0.0 || overlap2 < 0.0) continue;
             if (antialias == 1)
-              *(npy_float64*) PyArray_GETPTR2(my_array, j, i) += (dsp*overlap1)*overlap2;
+              *(npy_float64*) PyArray_GETPTR2(my_array, j, i) += overlap1*overlap2;
             else *(npy_float64*) PyArray_GETPTR2(my_array, j, i) = dsp;
           }
         }



More information about the yt-svn mailing list