From 830c81681d646b04a5d8d8379f94cb757556f818 Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Wed, 13 Jun 2018 09:39:53 +0100 Subject: [PATCH 1/2] Performance improvement in compute_statistic --- glue/utils/array.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/glue/utils/array.py b/glue/utils/array.py index 5492ba9f6..d21e43f4d 100644 --- a/glue/utils/array.py +++ b/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 From 7d06f1e6f7856da4b79d803a8060e60126b4a59d Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Wed, 13 Jun 2018 09:41:01 +0100 Subject: [PATCH 2/2] Improve performance of RoiSubsetState for whole datasets when selecting along pixel coordinates, by adding a way for ROIs to return masks rather than treating each pixel as an individual point --- glue/core/roi.py | 61 ++++++++++++++++++++++++++++++++++--- glue/core/subset.py | 17 ++++++----- glue/core/tests/test_roi.py | 31 +++++++++++++++++-- glue/utils/geometry.py | 28 ++++++++++++++++- 4 files changed, 120 insertions(+), 17 deletions(-) diff --git a/glue/core/roi.py b/glue/core/roi.py index 4cdeffede..cb76f5ce7 100644 --- a/glue/core/roi.py +++ b/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 diff --git a/glue/core/subset.py b/glue/core/subset.py index 893a4cbbe..b4fc744f5 100644 --- a/glue/core/subset.py +++ b/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) diff --git a/glue/core/tests/test_roi.py b/glue/core/tests/test_roi.py index d631388ee..16991fc95 100644 --- a/glue/core/tests/test_roi.py +++ b/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) diff --git a/glue/utils/geometry.py b/glue/utils/geometry.py index 9b11b7c07..c2671a07c 100644 --- a/glue/utils/geometry.py +++ b/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):