Performance improvements

glue/core/roi.py

@@ -3,11 +3,12 @@
 import copy
 
 import numpy as np
+
 from matplotlib.patches import Ellipse, Polygon, Rectangle, Path as MplPath, PathPatch
 from matplotlib.transforms import IdentityTransform, blended_transform_factory
 
 from glue.core.exceptions import UndefinedROI
-from glue.utils import points_inside_poly, iterate_chunks
+from glue.utils import points_inside_poly, polygon_mask, iterate_chunks, broadcast_to
 
 
 np.seterr(all='ignore')
@@ -61,6 +62,18 @@ class Roi(object):  # pragma: no cover
     method to test whether a collection of 2D points lies inside the region.
     """
 
+    def mask(self, shape):
+        """
+        Return a mask which is the projection of the ROI onto pixel space.
+        The mask is True if the center of the pixel is inside the ROI.
+
+        Parameters
+        ----------
+        shape : tuple
+            The shape of the array to project the ROI onto
+        """
+        raise NotImplementedError()
+
     def contains(self, x, y):
         """Return true/false for each x/y pair.
 
@@ -135,8 +148,11 @@ def __init__(self, x=None, y=None):
         self.x = x
         self.y = y
 
+    def mask(self, shape):
+        return broadcast_to(False, shape)
+
     def contains(self, x, y):
-        return False
+        return broadcast_to(False, np.shape(x))
 
     def move_to(self, x, y):
         self.x = x
@@ -208,6 +224,15 @@ def width(self):
     def height(self):
         return self.ymax - self.ymin
 
+    def mask(self, shape):
+        mask = np.zeros(shape, dtype=bool)
+        imin = int(np.ceil(self.xmin))
+        imax = int(np.floor(self.xmax))
+        jmin = int(np.ceil(self.ymin))
+        jmax = int(np.floor(self.ymax))
+        mask[jmin:jmax + 1, imin:imax + 1] = True
+        return mask
+
     def contains(self, x, y):
         """
         Test whether a set of (x,y) points falls within
@@ -225,8 +250,8 @@ def contains(self, x, y):
         if not self.defined():
             raise UndefinedROI
 
-        return (x > self.xmin) & (x < self.xmax) & \
-               (y > self.ymin) & (y < self.ymax)
+        return (x >= self.xmin) & (x <= self.xmax) & \
+               (y >= self.ymin) & (y <= self.ymax)
 
     def update_limits(self, xmin, ymin, xmax, ymax):
         """
@@ -317,12 +342,25 @@ def move_to(self, center):
         self.min += delta
         self.max += delta
 
+    def mask(self, shape):
+        mask = np.zeros(shape, dtype=bool)
+        if self.ori == 'x':
+            imin = int(np.ceil(self.min))
+            imax = int(np.floor(self.max))
+            mask[:, imin:imax + 1] = True
+        else:
+            jmin = int(np.ceil(self.min))
+            jmax = int(np.floor(self.max))
+            print(jmin, jmax)
+            mask[jmin:jmax + 1, :] = True
+        return mask
+
     def contains(self, x, y):
         if not self.defined():
             raise UndefinedROI()
 
         coord = x if self.ori == 'x' else y
-        return (coord > self.min) & (coord < self.max)
+        return (coord >= self.min) & (coord <= self.max)
 
     def transformed(self, xfunc=None, yfunc=None):
         if self.ori == 'x':
@@ -381,6 +419,16 @@ def __init__(self, xc=None, yc=None, radius=None):
         self.yc = yc
         self.radius = radius
 
+    def mask(self, shape):
+        mask = np.zeros(shape, dtype=bool)
+        imin = int(np.ceil(self.xc - self.radius))
+        imax = int(np.floor(self.xc + self.radius))
+        jmin = int(np.ceil(self.yc - self.radius))
+        jmax = int(np.floor(self.yc + self.radius))
+        x, y = np.meshgrid(np.arange(imin, imax + 1), np.arange(jmin, jmax + 1))
+        mask[jmin:jmax + 1, imin:imax + 1] = self.contains(x, y)
+        return mask
+
     def contains(self, x, y):
         """
         Test whether a set of (x,y) points falls within
@@ -557,6 +605,9 @@ def __str__(self):
         result += ')'
         return result
 
+    def mask(self, shape):
+        return polygon_mask(shape, self.vx, self.vy)
+
     def contains(self, x, y):
         """
         Test whether a set of (x,y) points falls within

glue/core/subset.py

@@ -570,20 +570,21 @@ def to_mask(self, data, view=None):
             # Note that we can only do this if the view (if present) preserved
             # the dimensionality, which is why we checked that x.ndim == data.ndim
 
-            subset = []
+            final_shape = []
             for i in range(data.ndim):
                 if i == self.xatt.axis or i == self.yatt.axis:
-                    subset.append(slice(None))
+                    final_shape.append(data.shape[i])
                 else:
-                    subset.append(slice(0, 1))
-
-            x_slice = x[subset]
-            y_slice = y[subset]
+                    final_shape.append(1)
 
             if self.roi.defined():
-                result = self.roi.contains(x_slice, y_slice)
+                slice_shape = data.shape[self.yatt.axis], data.shape[self.xatt.axis]
+                result = self.roi.mask(slice_shape)
+                if self.xatt.axis < self.yatt.axis:
+                    result = result.T
+                result = result.reshape(final_shape)
             else:
-                result = np.zeros(x_slice.shape, dtype=bool)
+                result = False
 
             result = broadcast_to(result, x.shape)
 

glue/core/tests/test_roi.py

@@ -6,7 +6,7 @@
 import numpy as np
 from mock import MagicMock
 from matplotlib.figure import Figure
-from numpy.testing import assert_almost_equal
+from numpy.testing import assert_almost_equal, assert_equal
 
 from .. import roi as r
 from ..component import CategoricalComponent
@@ -151,7 +151,7 @@ def test_contains(self):
         x = np.array([0, 1, 2, 3])
         y = np.array([-np.inf, 100, 200, 0])
         np.testing.assert_array_equal(roi.contains(x, y),
-                                      [False, False, True, False])
+                                      [False, True, True, True])
 
     def test_contains_undefined(self):
         roi = XRangeROI()
@@ -192,7 +192,7 @@ def test_contains(self):
         y = np.array([0, 1, 2, 3])
         x = np.array([-np.inf, 100, 200, 0])
         np.testing.assert_array_equal(roi.contains(x, y),
-                                      [False, False, True, False])
+                                      [False, True, True, True])
 
     def test_contains_undefined(self):
         roi = YRangeROI()
@@ -998,3 +998,28 @@ def test_data_to_pixel(self):
 
 
 del TestMpl  # prevents unittest discovery from running abstract base class
+
+
+ROIS = [PointROI(4.1, 4.3),
+        RectangularROI(12.1, 14.3, 5.5, 10),
+        XRangeROI(10.2, 30),
+        YRangeROI(10, 22.2),
+        CircularROI(5.3, 8.3, 3.8),
+        PolygonalROI([4.4, 8.9, 8.8, 4.3], [3.4, 3.3, 9.9, 9.7])]
+
+
+@pytest.mark.parametrize('roi', ROIS)
+def test_mask(roi):
+
+    shape = (100, 100)
+    y, x = np.indices((100, 100))
+
+    mask1 = roi.contains(x, y)
+    mask2 = roi.mask(shape)
+
+    from astropy.io import fits
+
+    fits.writeto('mask1_{0}.fits'.format(roi.__class__.__name__), mask1.astype(int), overwrite=True)
+    fits.writeto('mask2_{0}.fits'.format(roi.__class__.__name__), mask2.astype(int), overwrite=True)
+
+    assert_equal(mask1, mask2)

glue/utils/array.py

@@ -449,10 +449,10 @@ def compute_statistic(statistic, data, mask=None, axis=None, finite=True,
 
     if (finite or positive or mask is not None) and data.dtype.kind != 'M':
 
-        keep = np.ones(data.shape, dtype=bool)
-
         if finite:
-            keep &= np.isfinite(data)
+            keep = np.isfinite(data)
+        else:
+            keep = np.ones(data.shape, dtype=bool)
 
         if positive:
             keep &= data > 0

glue/utils/geometry.py

@@ -4,7 +4,33 @@
 
 from glue.utils import unbroadcast, broadcast_to
 
-__all__ = ['points_inside_poly', 'polygon_line_intersections', 'floodfill']
+__all__ = ['polygon_mask', 'points_inside_poly', 'polygon_line_intersections', 'floodfill']
+
+
+def polygon_mask(shape, vx, vy):
+    """
+    Given an array shape and polygon vertices, return a boolean mask indicating
+    which values are inside the polygon.
+    """
+
+    vx = np.asanyarray(vx)
+    vy = np.asanyarray(vy)
+
+    xmin = int(np.floor(np.min(vx)))
+    xmax = int(np.ceil(np.max(vx)))
+    ymin = int(np.floor(np.min(vy)))
+    ymax = int(np.ceil(np.max(vy)))
+
+    ix = np.arange(xmin, xmax + 1, dtype=int)
+    iy = np.arange(ymin, ymax + 1, dtype=int)
+
+    x, y = np.meshgrid(ix, iy)
+
+    inside = np.zeros(shape, bool)
+
+    inside[ymin:ymax + 1, xmin:xmax + 1] = points_inside_poly(x, y, vx, vy)
+
+    return inside
 
 
 def points_inside_poly(x, y, vx, vy):