From 67526480c575684f71ea456d6cbf28ce31acb593 Mon Sep 17 00:00:00 2001
From: Derek Homeier <dhomeie@gwdg.de>
Date: Sun, 26 Sep 2021 02:36:59 +0200
Subject: [PATCH 01/12] ENH: Add angle of rotation to RectangularROI and
 EllipticalROI

---
 glue/core/roi.py | 136 ++++++++++++++++++++++++++++++++++++++++-------
 1 file changed, 118 insertions(+), 18 deletions(-)

diff --git a/glue/core/roi.py b/glue/core/roi.py
index 3a3656836..573e7a25b 100644
--- a/glue/core/roi.py
+++ b/glue/core/roi.py
@@ -55,6 +55,12 @@ def pixel_to_axes(axes, x, y):
     return axes.transAxes.inverted().transform(xy)
 
 
+def rotation(alpha):
+    """Return rotation matrix for angle alpha around origin.
+    """
+    return np.array([[np.cos(alpha), np.sin(alpha)], [-np.sin(alpha), np.cos(alpha)]])
+
+
 class Roi(object):  # pragma: no cover
 
     """
@@ -163,21 +169,37 @@ class RectangularROI(Roi):
 
     """
     A 2D rectangular region of interest.
+
+    Parameters
+    ----------
+    xmin, xmax :  float, optional
+        x coordinates of left and right border
+    ymin, ymax :  float, optional
+        y coordinates of lower and upper border
+    theta : float, optional
+        Angle of clockwise rotation around center
     """
 
-    def __init__(self, xmin=None, xmax=None, ymin=None, ymax=None):
+    def __init__(self, xmin=None, xmax=None, ymin=None, ymax=None, theta=None):
         super(RectangularROI, self).__init__()
         self.xmin = xmin
         self.xmax = xmax
         self.ymin = ymin
         self.ymax = ymax
+        self.theta = theta
+        if self.theta is None or self.theta == 0.0:
+            self.rotation = np.identity(2)
+        else:
+            self.rotation = rotation(self.theta)
+        self.invrot = self.rotation * [[1, -1], [-1, 1]]
 
     def __str__(self):
-        if self.defined():
-            return "x=[%0.3f, %0.3f], y=[%0.3f, %0.3f]" % (self.xmin,
-                                                           self.xmax,
-                                                           self.ymin,
-                                                           self.ymax)
+        if self.defined() and self.theta is None:
+            return f"x=[{self.xmin:.3f}, {self.xmax:.3f}], y=[{self.ymin:.3f}, {self.ymax:.3f}]"
+        elif self.defined():
+            return (f"center=({self.center()[0]:.3f}, {self.center()[1]:.3f}), "
+                    f"size=({self.width():.3f} x {self.height():.3f}), "
+                    f"theta={self.theta:.3f} radian")
         else:
             return "Undefined Rectangular ROI"
 
@@ -229,8 +251,26 @@ def contains(self, x, y):
         if not self.defined():
             raise UndefinedROI
 
-        return (x > self.xmin) & (x < self.xmax) & \
-               (y > self.ymin) & (y < self.ymax)
+        if self.theta is not None and self.theta != 0.0:
+            if not isinstance(x, np.ndarray):
+                x = np.asarray(x)
+            if not isinstance(y, np.ndarray):
+                y = np.asarray(y)
+
+            inside = np.zeros_like(x, dtype=bool)
+            bounds = self.to_polygon()
+            xc, yc = self.center()
+            keep = ((x >= bounds[0].min()) & (x <= bounds[0].max()) &
+                    (y >= bounds[1].min()) & (y <= bounds[1].max()))
+            x = x[keep] - xc
+            y = y[keep] - yc
+            shape = (2,) + x.shape
+            x, y = (self.invrot @ [x.flatten(), y.flatten()]).reshape(shape)
+            inside[keep] = (abs(x) <= self.width() / 2) & (abs(y) <= self.height() / 2)
+            return inside
+
+        else:
+            return (x > self.xmin) & (x < self.xmax) & (y > self.ymin) & (y < self.ymax)
 
     def update_limits(self, xmin, ymin, xmax, ymax):
         """
@@ -255,8 +295,13 @@ def defined(self):
 
     def to_polygon(self):
         if self.defined():
-            return (np.array([self.xmin, self.xmax, self.xmax, self.xmin, self.xmin]),
-                    np.array([self.ymin, self.ymin, self.ymax, self.ymax, self.ymin]))
+            if self.theta is None or self.theta == 0.0:
+                return (np.array([self.xmin, self.xmax, self.xmax, self.xmin, self.xmin]),
+                        np.array([self.ymin, self.ymin, self.ymax, self.ymax, self.ymin]))
+            else:
+                corners = (np.array([-1, 1, 1, -1, -1]) * self.width() / 2,
+                           np.array([-1, -1, 1, 1, -1]) * self.height() / 2)
+            return tuple((self.rotation @ corners) + np.array(self.center()).reshape((2, 1)))
         else:
             return [], []
 
@@ -471,25 +516,50 @@ def __setgluestate__(cls, rec, context):
 
 class EllipticalROI(Roi):
     """
-    A 2D elliptical region of interest.
+    A 2D elliptical region of interest at (xc, yc) with semimajor/minor axes radius_[xy].
+
+    Parameters
+    ----------
+    xc :  float, optional
+        x coordinate of center
+    yc :  float, optional
+        y coordinate of center
+    radius_x :  float, optional
+        Semiaxis along x axis
+    radius_y :  float, optional
+        Semiaxis along y axis
+    theta : float, optional
+        Angle of clockwise rotation around (xc, yc)
     """
 
     def __init__(self, xc=None, yc=None, radius_x=None, radius_y=None, theta=None):
         super(EllipticalROI, self).__init__()
-        if theta is not None and theta != 0:
-            raise NotImplementedError("Rotated ellipses are not yet supported")
         self.xc = xc
         self.yc = yc
         self.radius_x = radius_x
         self.radius_y = radius_y
+        self.theta = theta
+        if self.theta is None or self.theta == 0.0:
+            self.rotation = np.identity(2)
+        else:
+            self.rotation = rotation(self.theta)
+        self.invrot = self.rotation * [[1, -1], [-1, 1]]
+
+    def __str__(self):
+        if self.defined():
+            return (f"center=({self.xc:.3f}, {self.yc:.3f}), "
+                    f"semiaxes=({self.radius_x:.3f} x {self.radius_y:.3f}), "
+                    f"theta={self.theta:.3f} radian")
+        else:
+            return "Undefined Elliptical ROI"
 
     def contains(self, x, y):
         """
         Test whether a set of (x,y) points falls within
         the region of interest
 
-        :param x: A list of x points
-        :param y: A list of y points
+        :param x: A list of x coordinates
+        :param y: A list of y coordinates
 
         *Returns*
 
@@ -504,6 +574,19 @@ def contains(self, x, y):
             x = np.asarray(x)
         if not isinstance(y, np.ndarray):
             y = np.asarray(y)
+
+        if self.theta is not None:
+            inside = np.zeros_like(x, dtype=bool)
+            bounds = self.bounds()
+            keep = ((x >= bounds[0][0]) & (x <= bounds[0][1]) &
+                    (y >= bounds[1][0]) & (y <= bounds[1][1]))
+            x = x[keep] - self.xc
+            y = y[keep] - self.yc
+            shape = (2,) + x.shape
+            x, y = (self.invrot @ [x.flatten(), y.flatten()]).reshape(shape)
+            inside[keep] = ((x ** 2 / self.radius_x ** 2 + y ** 2 / self.radius_y ** 2) < 1.)
+            return inside
+
         return (((x - self.xc) ** 2 / self.radius_x ** 2 +
                  (y - self.yc) ** 2 / self.radius_y ** 2) < 1.)
 
@@ -523,14 +606,27 @@ def defined(self):
                 self.radius_x is not None and
                 self.radius_y is not None)
 
+    def get_center(self):
+        return self.xc, self.yc
+
     def to_polygon(self):
         """ Returns x, y, where each is a list of points """
         if not self.defined():
             return [], []
         theta = np.linspace(0, 2 * np.pi, num=20)
-        x = self.xc + self.radius_x * np.cos(theta)
-        y = self.yc + self.radius_y * np.sin(theta)
-        return x, y
+        x = self.radius_x * np.cos(theta)
+        y = self.radius_y * np.sin(theta)
+        x, y = self.rotation @ (x, y)
+        return x + self.xc, y + self.yc
+
+    def bounds(self):
+        """ Returns (conservatively estimated) boundary values in x and y """
+        if self.theta is None or self.theta == 0.0:
+            return [[self.xc - self.radius_x, self.xc + self.radius_x],
+                    [self.yc - self.radius_y, self.yc + self.radius_y]]
+        radius = max(self.radius_x, self.radius_y)
+        return [[self.xc - radius, self.xc + radius],
+                [self.yc - radius, self.yc + radius]]
 
     def transformed(self, xfunc=None, yfunc=None):
         return PolygonalROI(*self.to_polygon()).transformed(xfunc=xfunc, yfunc=yfunc)
@@ -674,6 +770,10 @@ def contains(self, x, y):
         result = points_inside_poly(x, y, self.vx, self.vy)
         return result
 
+    # There are several possible definitions of the centre; this is easiest to calculate.
+    def center(self):
+        return np.mean(self.vx), np.mean(self.vy)
+
     def move_to(self, xdelta, ydelta):
         self.vx = list(map(lambda x: x + xdelta, self.vx))
         self.vy = list(map(lambda y: y + ydelta, self.vy))

From 28db2fa603aa20fafc967387af62fdc4f8609aef Mon Sep 17 00:00:00 2001
From: Derek Homeier <dhomeie@gwdg.de>
Date: Sun, 26 Sep 2021 02:37:12 +0200
Subject: [PATCH 02/12] TST: first tests for RectangularROI, EllipticalROI with
 theta

---
 glue/core/tests/test_roi.py    | 13 +++++++++++--
 glue/core/tests/test_subset.py | 16 +++++++++++++++-
 2 files changed, 26 insertions(+), 3 deletions(-)

diff --git a/glue/core/tests/test_roi.py b/glue/core/tests/test_roi.py
index e43b81257..b44dff9c7 100644
--- a/glue/core/tests/test_roi.py
+++ b/glue/core/tests/test_roi.py
@@ -325,15 +325,19 @@ class TestEllipse(object):
     def setup_method(self, method):
         self.roi_empty = EllipticalROI()
         self.roi = EllipticalROI(1, 2, 3, 4)
+        self.roi_rotated = EllipticalROI(1, 2, 3, 0.4, theta=np.pi / 6)
 
     def test_undefined_on_creation(self):
         assert not self.roi_empty.defined()
         assert self.roi.defined()
+        assert self.roi_rotated.defined()
 
     def test_contains_on_undefined_contains_raises(self):
         with pytest.raises(UndefinedROI):
             self.roi_empty.contains(1, 1)
-        self.roi.contains(1, 1)
+        assert self.roi.contains(1, 1)
+        assert self.roi_rotated.contains(0, 2.5)
+        assert not self.roi_rotated.contains(0, 2)
 
     def test_set_center(self):
         assert self.roi.contains(0, 0)
@@ -366,6 +370,11 @@ def test_poly(self):
         assert poly.contains(0, 0)
         assert not poly.contains(10, 0)
 
+        x, y = self.roi_rotated.to_polygon()
+        poly = PolygonalROI(vx=x, vy=y)
+        assert poly.contains(0, 2.5)
+        assert not poly.contains(0, 2)
+
     def test_poly_undefined(self):
         x, y = self.roi_empty.to_polygon()
         assert x == []
@@ -931,7 +940,7 @@ def test_circular_roi_representation():
 
     # Case 2: linear-linear axes with non-square aspect ratio
 
-    ax.set_xlim(-40, 100)
+    ax.set_xlim(10, 100)
     ax.set_ylim(-2, 5)
     ax.set_xscale('linear')
     ax.set_yscale('linear')
diff --git a/glue/core/tests/test_subset.py b/glue/core/tests/test_subset.py
index e953efc0b..9effaed60 100644
--- a/glue/core/tests/test_subset.py
+++ b/glue/core/tests/test_subset.py
@@ -851,7 +851,8 @@ def test_roi_subset_state(self):
         roi = RectangularROI(xmin=0, xmax=3, ymin=1.1, ymax=1.4)
 
         subset = self.data.new_subset()
-        subset.subset_state = RoiSubsetState(xatt=self.data.id['a'], yatt=self.data.id['c'], roi=roi)
+        subset.subset_state = RoiSubsetState(xatt=self.data.id['a'], yatt=self.data.id['c'],
+                                             roi=roi)
         assert_equal(self.data.subsets[0].to_mask(), [0, 1, 0, 0])
 
         data_clone = clone(self.data)
@@ -871,6 +872,19 @@ def test_roi_subset_state_with_pretransform(self):
 
         assert_equal(data_clone.subsets[0].to_mask(), [0, 1, 0, 0])
 
+    def test_roi_subset_state_rotated(self):
+
+        roi = RectangularROI(xmin=0, xmax=3, ymin=1.1, ymax=1.4, theta=np.pi / 3)
+
+        subset = self.data.new_subset()
+        subset.subset_state = RoiSubsetState(xatt=self.data.id['a'], yatt=self.data.id['c'],
+                                             roi=roi)
+        assert_equal(self.data.subsets[0].to_mask(), [0, 0, 0, 1])
+
+        data_clone = clone(self.data)
+
+        assert_equal(data_clone.subsets[0].to_mask(), [0, 1, 0, 0])
+
 
 @requires_scipy
 def test_floodfill_subset_state():

From c48423f8ed2b2e05c2e01808415674632cae21b4 Mon Sep 17 00:00:00 2001
From: Derek Homeier <dhomeie@gwdg.de>
Date: Tue, 28 Sep 2021 15:22:42 +0200
Subject: [PATCH 03/12] Speed up 90-deg rot; add `rotate_[to|by]` methods;
 theta in gluestate

---
 glue/core/roi.py | 156 +++++++++++++++++++++++++++++++++--------------
 1 file changed, 110 insertions(+), 46 deletions(-)

diff --git a/glue/core/roi.py b/glue/core/roi.py
index 573e7a25b..a00586c2c 100644
--- a/glue/core/roi.py
+++ b/glue/core/roi.py
@@ -117,6 +117,14 @@ def move_to(self, x, y):
         """Translate the ROI to a center of (x, y)"""
         raise NotImplementedError()
 
+    def rotate_to(self, theta):
+        """Set rotation angle of ROI around center to theta"""
+        raise NotImplementedError()
+
+    def rotate_by(self, dtheta):
+        """Rotate the ROI around center by angle dtheta"""
+        raise NotImplementedError()
+
     def defined(self):
         """ Returns whether or not the subset is properly defined """
         raise NotImplementedError()
@@ -186,15 +194,15 @@ def __init__(self, xmin=None, xmax=None, ymin=None, ymax=None, theta=None):
         self.xmax = xmax
         self.ymin = ymin
         self.ymax = ymax
-        self.theta = theta
-        if self.theta is None or self.theta == 0.0:
+        self.theta = 0 if theta is None else theta
+        if np.isclose(self.theta % np.pi, 0.0, atol=1e-9):
             self.rotation = np.identity(2)
         else:
             self.rotation = rotation(self.theta)
         self.invrot = self.rotation * [[1, -1], [-1, 1]]
 
     def __str__(self):
-        if self.defined() and self.theta is None:
+        if self.defined() and self.theta == 0:
             return f"x=[{self.xmin:.3f}, {self.xmax:.3f}], y=[{self.ymin:.3f}, {self.ymax:.3f}]"
         elif self.defined():
             return (f"center=({self.center()[0]:.3f}, {self.center()[1]:.3f}), "
@@ -215,6 +223,17 @@ def move_to(self, x, y):
         self.ymin += dy
         self.ymax += dy
 
+    def rotate_to(self, theta):
+        self.theta = 0 if theta is None else theta
+        if np.isclose(self.theta % np.pi, 0.0, atol=1e-9):
+            self.rotation = np.identity(2)
+        else:
+            self.rotation = rotation(self.theta)
+        self.invrot = self.rotation * [[1, -1], [-1, 1]]
+
+    def rotate_by(self, dtheta):
+        self.rotate_to(self.theta + dtheta)
+
     def transpose(self, copy=True):
         if copy:
             new = self.copy()
@@ -239,8 +258,8 @@ def contains(self, x, y):
         Test whether a set of (x,y) points falls within
         the region of interest
 
-        :param x: A scalar or numpy array of x points
-        :param y: A scalar or numpy array of y points
+        :param x: A scalar or iterable of x coordinates
+        :param y: A scalar or iterable of y coordinates
 
         *Returns*
 
@@ -251,26 +270,30 @@ def contains(self, x, y):
         if not self.defined():
             raise UndefinedROI
 
-        if self.theta is not None and self.theta != 0.0:
-            if not isinstance(x, np.ndarray):
-                x = np.asarray(x)
-            if not isinstance(y, np.ndarray):
-                y = np.asarray(y)
-
-            inside = np.zeros_like(x, dtype=bool)
-            bounds = self.to_polygon()
-            xc, yc = self.center()
-            keep = ((x >= bounds[0].min()) & (x <= bounds[0].max()) &
-                    (y >= bounds[1].min()) & (y <= bounds[1].max()))
-            x = x[keep] - xc
-            y = y[keep] - yc
-            shape = (2,) + x.shape
-            x, y = (self.invrot @ [x.flatten(), y.flatten()]).reshape(shape)
-            inside[keep] = (abs(x) <= self.width() / 2) & (abs(y) <= self.height() / 2)
-            return inside
+        if not isinstance(x, np.ndarray):
+            x = np.asarray(x)
+        if not isinstance(y, np.ndarray):
+            y = np.asarray(y)
 
-        else:
+        if np.isclose(self.theta % np.pi, 0.0, atol=1e-9):
             return (x > self.xmin) & (x < self.xmax) & (y > self.ymin) & (y < self.ymax)
+        elif np.isclose(self.theta % (np.pi / 2), 0.0, atol=1e-9):
+            xc, yc = self.center()
+            xext = self.height() / 2
+            yext = self.width() / 2
+            return (x > xc - xext) & (x < xc + xext) & (y > yc - yext) & (y < yc + yext)
+
+        inside = np.zeros_like(x, dtype=bool)
+        bounds = self.to_polygon()
+        xc, yc = self.center()
+        keep = ((x >= bounds[0].min()) & (x <= bounds[0].max()) &
+                (y >= bounds[1].min()) & (y <= bounds[1].max()))
+        x = x[keep] - xc
+        y = y[keep] - yc
+        shape = (2,) + x.shape
+        x, y = (self.invrot @ [x.flatten(), y.flatten()]).reshape(shape)
+        inside[keep] = (abs(x) <= self.width() / 2) & (abs(y) <= self.height() / 2)
+        return inside
 
     def update_limits(self, xmin, ymin, xmax, ymax):
         """
@@ -295,7 +318,7 @@ def defined(self):
 
     def to_polygon(self):
         if self.defined():
-            if self.theta is None or self.theta == 0.0:
+            if np.isclose(self.theta % np.pi, 0.0, atol=1e-9):
                 return (np.array([self.xmin, self.xmax, self.xmax, self.xmin, self.xmin]),
                         np.array([self.ymin, self.ymin, self.ymax, self.ymax, self.ymin]))
             else:
@@ -316,12 +339,14 @@ def __gluestate__(self, context):
         return dict(xmin=context.do(self.xmin),
                     xmax=context.do(self.xmax),
                     ymin=context.do(self.ymin),
-                    ymax=context.do(self.ymax))
+                    ymax=context.do(self.ymax),
+                    theta=context.do(self.theta))
 
     @classmethod
     def __setgluestate__(cls, rec, context):
         return cls(xmin=context.object(rec['xmin']), xmax=context.object(rec['xmax']),
-                   ymin=context.object(rec['ymin']), ymax=context.object(rec['ymax']))
+                   ymin=context.object(rec['ymin']), ymax=context.object(rec['ymax']),
+                   theta=context.object(rec['theta']))
 
 
 class RangeROI(Roi):
@@ -438,8 +463,8 @@ def contains(self, x, y):
         Test whether a set of (x,y) points falls within
         the region of interest
 
-        :param x: A list of x points
-        :param y: A list of y points
+        :param x: A scalar or iterable of x coordinates
+        :param y: A scalar or iterable of y coordinates
 
         *Returns*
 
@@ -538,8 +563,8 @@ def __init__(self, xc=None, yc=None, radius_x=None, radius_y=None, theta=None):
         self.yc = yc
         self.radius_x = radius_x
         self.radius_y = radius_y
-        self.theta = theta
-        if self.theta is None or self.theta == 0.0:
+        self.theta = 0 if theta is None else theta
+        if np.isclose(self.theta % np.pi, 0.0, atol=1e-9):
             self.rotation = np.identity(2)
         else:
             self.rotation = rotation(self.theta)
@@ -558,8 +583,8 @@ def contains(self, x, y):
         Test whether a set of (x,y) points falls within
         the region of interest
 
-        :param x: A list of x coordinates
-        :param y: A list of y coordinates
+        :param x: A scalar or iterable of x coordinates
+        :param y: A scalar or iterable of y coordinates
 
         *Returns*
 
@@ -575,7 +600,16 @@ def contains(self, x, y):
         if not isinstance(y, np.ndarray):
             y = np.asarray(y)
 
-        if self.theta is not None:
+        if np.isclose(self.theta % np.pi, 0.0, atol=1e-9):
+            return (((x - self.xc) ** 2 / self.radius_x ** 2 +
+                     (y - self.yc) ** 2 / self.radius_y ** 2) < 1.)
+        elif np.isclose(self.theta % (np.pi / 2), 0.0, atol=1e-9):
+            return (((x - self.xc) ** 2 / self.radius_y ** 2 +
+                     (y - self.yc) ** 2 / self.radius_x ** 2) < 1.)
+        else:
+            # Pre-select points inside the bounding rectangle. In principle this could be
+            # used to speed up the non-rotated cases above as well, but will only pay off
+            # the overhead for datasets much larger than the region (e.g. < ~1 % inside).
             inside = np.zeros_like(x, dtype=bool)
             bounds = self.bounds()
             keep = ((x >= bounds[0][0]) & (x <= bounds[0][1]) &
@@ -587,9 +621,6 @@ def contains(self, x, y):
             inside[keep] = ((x ** 2 / self.radius_x ** 2 + y ** 2 / self.radius_y ** 2) < 1.)
             return inside
 
-        return (((x - self.xc) ** 2 / self.radius_x ** 2 +
-                 (y - self.yc) ** 2 / self.radius_y ** 2) < 1.)
-
     def reset(self):
         """
         Reset the rectangular region.
@@ -621,12 +652,16 @@ def to_polygon(self):
 
     def bounds(self):
         """ Returns (conservatively estimated) boundary values in x and y """
-        if self.theta is None or self.theta == 0.0:
+        if self.theta is None or np.isclose(self.theta % (np.pi), 0.0, atol=1e-9):
             return [[self.xc - self.radius_x, self.xc + self.radius_x],
                     [self.yc - self.radius_y, self.yc + self.radius_y]]
-        radius = max(self.radius_x, self.radius_y)
-        return [[self.xc - radius, self.xc + radius],
-                [self.yc - radius, self.yc + radius]]
+        elif np.isclose(self.theta % (np.pi / 2), 0.0, atol=1e-9):
+            return [[self.xc - self.radius_y, self.xc + self.radius_y],
+                    [self.yc - self.radius_x, self.yc + self.radius_x]]
+        else:
+            radius = max(self.radius_x, self.radius_y)
+            return [[self.xc - radius, self.xc + radius],
+                    [self.yc - radius, self.yc + radius]]
 
     def transformed(self, xfunc=None, yfunc=None):
         return PolygonalROI(*self.to_polygon()).transformed(xfunc=xfunc, yfunc=yfunc)
@@ -635,16 +670,28 @@ def move_to(self, xdelta, ydelta):
         self.xc += xdelta
         self.yc += ydelta
 
+    def rotate_to(self, theta):
+        self.theta = 0 if theta is None else theta
+        if np.isclose(self.theta % np.pi, 0.0, atol=1e-9):
+            self.rotation = np.identity(2)
+        else:
+            self.rotation = rotation(self.theta)
+        self.invrot = self.rotation * [[1, -1], [-1, 1]]
+
+    def rotate_by(self, dtheta):
+        self.rotate_to(self.theta + dtheta)
+
     def __gluestate__(self, context):
         return dict(xc=context.do(self.xc),
                     yc=context.do(self.yc),
                     radius_x=context.do(self.radius_x),
-                    radius_y=context.do(self.radius_y))
+                    radius_y=context.do(self.radius_y),
+                    theta=context.do(self.theta))
 
     @classmethod
     def __setgluestate__(cls, rec, context):
         return cls(xc=rec['xc'], yc=rec['yc'],
-                   radius_x=rec['radius_x'], radius_y=rec['radius_y'])
+                   radius_x=rec['radius_x'], radius_y=rec['radius_y'], theta=rec['theta'])
 
 
 class VertexROIBase(Roi):
@@ -663,6 +710,7 @@ def __init__(self, vx=None, vy=None):
             self.vx = []
         if self.vy is None:
             self.vy = []
+        self.theta = 0
 
     def transformed(self, xfunc=None, yfunc=None):
         vx = self.vx if xfunc is None else xfunc(np.asarray(self.vx))
@@ -729,8 +777,7 @@ def __gluestate__(self, context):
 
     @classmethod
     def __setgluestate__(cls, rec, context):
-        return cls(vx=context.object(rec['vx']),
-                   vy=context.object(rec['vy']))
+        return cls(vx=context.object(rec['vx']), vy=context.object(rec['vy']))
 
 
 class PolygonalROI(VertexROIBase):
@@ -770,14 +817,31 @@ def contains(self, x, y):
         result = points_inside_poly(x, y, self.vx, self.vy)
         return result
 
-    # There are several possible definitions of the centre; this is easiest to calculate.
+    # There are several possible definitions of the centre; this is easiest to calculate
+    # (and invariant under rotation?) - do not include starting vertex twice!
     def center(self):
-        return np.mean(self.vx), np.mean(self.vy)
+        if not self.defined():
+            raise UndefinedROI
+        if self.vx[-1] == self.vx[0] and self.vy[:-1] == self.vy[0]:
+            return np.mean(self.vx[:-1]), np.mean(self.vy[:-1])
+        else:
+            return np.mean(self.vx), np.mean(self.vy)
 
     def move_to(self, xdelta, ydelta):
         self.vx = list(map(lambda x: x + xdelta, self.vx))
         self.vy = list(map(lambda y: y + ydelta, self.vy))
 
+    def rotate_to(self, theta):
+        theta = 0 if theta is None else theta
+        dtheta = theta - self.theta
+        if self.defined() and not np.isclose(dtheta % np.pi, 0.0, atol=1e-9):
+            dx, dy = np.array([self.vx, self.vy]) - np.array(self.center()).reshape(2, 1)
+            self.vx, self.vy = rotation(dtheta) @ (dx, dy) + np.array(self.center()).reshape(2, 1)
+        self.theta = theta
+
+    def rotate_by(self, dtheta):
+        self.rotate_to(self.theta + dtheta)
+
 
 class Projected3dROI(Roi):
     """"A region of interest defined in screen coordinates.

From 39a6c52c09a3e45a11fd17060765fe898398a147 Mon Sep 17 00:00:00 2001
From: Derek Homeier <dhomeie@gwdg.de>
Date: Tue, 28 Sep 2021 15:35:00 +0200
Subject: [PATCH 04/12] TST: add tests for rotated Rect, Elliptical, Poly ROIs

---
 glue/core/tests/test_roi.py    | 44 ++++++++++++++++++++++++++++++++++
 glue/core/tests/test_subset.py | 38 ++++++++++++++++++++++++++++-
 2 files changed, 81 insertions(+), 1 deletion(-)

diff --git a/glue/core/tests/test_roi.py b/glue/core/tests/test_roi.py
index b44dff9c7..b5b8b0e99 100644
--- a/glue/core/tests/test_roi.py
+++ b/glue/core/tests/test_roi.py
@@ -78,6 +78,17 @@ def test_reset(self):
         with pytest.raises(UndefinedROI):
             self.roi.contains(5, 5)
 
+    def test_set_rotation(self):
+        self.roi.update_limits(0, 0, 10, 4)
+        self.roi.rotate_by(-np.pi / 6)
+        assert self.roi.contains(8, 6)
+        self.roi.rotate_by(-np.pi / 3)
+        assert not self.roi.contains(9, 6)
+        assert self.roi.contains(6, 6.9)
+        self.roi.rotate_to(np.pi / 3)
+        assert not self.roi.contains(5, 6)
+        assert self.roi.contains(6, -3)
+
     def test_empty_to_polygon(self):
         x, y = self.roi.to_polygon()
         assert x == []
@@ -88,6 +99,12 @@ def test_to_polygon(self):
         x, y = self.roi.to_polygon()
         poly = PolygonalROI(vx=x, vy=y)
         assert poly.contains(5, 5)
+        self.roi.update_limits(0, 0, 10, 4)
+        self.roi.rotate_to(np.pi / 3)
+        x, y = self.roi.to_polygon()
+        poly = PolygonalROI(vx=x, vy=y)
+        assert not poly.contains(8, 3)
+        assert poly.contains(4, 6)
 
     def test_ndarray(self):
         self.roi.update_limits(0, 0, 10, 10)
@@ -357,6 +374,15 @@ def test_set_radius(self):
         assert self.roi.contains(0, 0)
         assert self.roi.contains(12, 12)
 
+    def test_set_rotation(self):
+        self.roi_rotated.rotate_to(0.55)
+        assert self.roi_rotated.contains(-1.5, 3.6)
+        self.roi_rotated.rotate_by(np.pi / 3)
+        assert self.roi_rotated.contains(1.2, 4.5)
+        self.roi_rotated.rotate_to(np.pi / 3)
+        assert not self.roi_rotated.contains(1.2, 4.5)
+        assert self.roi_rotated.contains(0.5, 3.5)
+
     def test_contains_many(self):
         x = [0, 0, 0, 0, 0]
         y = [0, 0, 0, 0, 0]
@@ -407,6 +433,12 @@ def test_serialization(self):
         assert_almost_equal(new_roi.radius_x, 3)
         assert_almost_equal(new_roi.radius_y, 4)
 
+    def test_serialize_rotated(self):
+        new_roi = roundtrip_roi(self.roi_rotated)
+        assert new_roi.radius_y == 0.4
+        assert new_roi.theta == np.pi / 6
+        assert new_roi.contains(-1.5, 3.5)
+
 
 class TestPolygon(object):
     def setup_method(self, method):
@@ -503,11 +535,23 @@ def test_str(self):
         """ __str__ returns a string """
         assert type(str(self.roi)) == str
 
+    def test_rotate(self):
+        self.define_as_square()
+        self.roi.rotate_to(np.pi / 4)
+        assert self.roi.contains([.5, .5, 1.2], [1.2, -0.2, .5]).all()
+        assert not self.roi.contains([1.5, 1.5], [0, 0]).any()
+
     def test_serialization(self):
         self.define_as_square()
         new_roi = roundtrip_roi(self.roi)
         assert_almost_equal(new_roi.vx, np.array([0, 0, 1, 1]))
         assert_almost_equal(new_roi.vy, np.array([0, 1, 1, 0]))
+        self.roi.rotate_to(np.pi / 4)
+        new_roi = roundtrip_roi(self.roi)
+        assert new_roi.theta == 0
+        sqh = np.sqrt(0.5)
+        assert_almost_equal(new_roi.vx, np.array([-sqh, 0, sqh, 0]) + 0.5)
+        assert_almost_equal(new_roi.vy, np.array([0, sqh, 0, -sqh]) + 0.5)
 
 
 class TestProjected3dROI(object):
diff --git a/glue/core/tests/test_subset.py b/glue/core/tests/test_subset.py
index 9effaed60..dd9718af1 100644
--- a/glue/core/tests/test_subset.py
+++ b/glue/core/tests/test_subset.py
@@ -872,6 +872,14 @@ def test_roi_subset_state_with_pretransform(self):
 
         assert_equal(data_clone.subsets[0].to_mask(), [0, 1, 0, 0])
 
+        roi.move_to(3, 0)
+        roi.rotate_to(1.5)
+
+        subset = data_clone.new_subset()
+        subset.subset_state = RoiSubsetState(xatt=data_clone.id['a'], yatt=data_clone.id['c'],
+                                             roi=roi, pretransform=example_transform)
+        assert_equal(data_clone.subsets[1].to_mask(), [0, 1, 0, 1])
+
     def test_roi_subset_state_rotated(self):
 
         roi = RectangularROI(xmin=0, xmax=3, ymin=1.1, ymax=1.4, theta=np.pi / 3)
@@ -883,7 +891,15 @@ def test_roi_subset_state_rotated(self):
 
         data_clone = clone(self.data)
 
-        assert_equal(data_clone.subsets[0].to_mask(), [0, 1, 0, 0])
+        assert_equal(data_clone.subsets[0].to_mask(), [0, 0, 0, 1])
+
+        roi = RectangularROI(xmin=0, xmax=3, ymin=1.1, ymax=1.4)
+        roi.rotate_by(0.9)
+
+        subset = data_clone.new_subset()
+        subset.subset_state = RoiSubsetState(xatt=data_clone.id['a'], yatt=data_clone.id['c'],
+                                             roi=roi)
+        assert_equal(data_clone.subsets[1].to_mask(), [0, 0, 0, 1])
 
 
 @requires_scipy
@@ -1059,3 +1075,23 @@ def test_roi_reduction():
     assert_equal(out[:, :, 0, 1], expected_slice)
     assert_equal(out[:, :, 1, 0], expected_slice)
     assert_equal(out[:, :, 1, 1], expected_slice)
+
+    roi.rotate_by(np.pi / 2)
+    state = RoiSubsetState(xatt=data4d.pixel_component_ids[0], yatt=data4d.pixel_component_ids[1], roi=roi)
+    out = state.to_mask(data4d)
+    expected_slice = np.array([[0, 0, 0, 0], [0, 0, 1, 1], [0, 0, 1, 1]])
+    assert_equal(out[:, :, 0, 0], expected_slice)
+    assert_equal(out[:, :, 0, 1], expected_slice)
+    assert_equal(out[:, :, 1, 0], expected_slice)
+    assert_equal(out[:, :, 1, 1], expected_slice)
+
+    roi = RectangularROI(0.5, 2.5, 1.9, 3.1)
+    roi.rotate_to(np.pi / 4)
+    state = RoiSubsetState(xatt=data4d.pixel_component_ids[0], yatt=data4d.pixel_component_ids[1], roi=roi)
+    out = state.to_mask(data4d)
+    expected_slice = np.array([[0, 0, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]])
+    assert_equal(out[:, :, 0, 0], expected_slice)
+    assert_equal(out[:, :, 0, 1], expected_slice)
+    assert_equal(out[:, :, 1, 0], expected_slice)
+    assert_equal(out[:, :, 1, 1], expected_slice)
+

From 7849cec7c86b28cbe247e1d898c3dd87ee448fd7 Mon Sep 17 00:00:00 2001
From: Derek Homeier <dhomeie@gwdg.de>
Date: Wed, 29 Sep 2021 13:14:45 +0200
Subject: [PATCH 05/12] Change sense of rotation to anticlockwise

---
 glue/core/roi.py               | 14 +++++++++-----
 glue/core/tests/test_roi.py    | 30 +++++++++++++++---------------
 glue/core/tests/test_subset.py |  8 ++++----
 3 files changed, 28 insertions(+), 24 deletions(-)

diff --git a/glue/core/roi.py b/glue/core/roi.py
index a00586c2c..16adb6dbb 100644
--- a/glue/core/roi.py
+++ b/glue/core/roi.py
@@ -56,9 +56,9 @@ def pixel_to_axes(axes, x, y):
 
 
 def rotation(alpha):
-    """Return rotation matrix for angle alpha around origin.
+    """Return rotation matrix for angle alpha (increasing anticlockwise) around origin.
     """
-    return np.array([[np.cos(alpha), np.sin(alpha)], [-np.sin(alpha), np.cos(alpha)]])
+    return np.array([[np.cos(alpha), -np.sin(alpha)], [np.sin(alpha), np.cos(alpha)]])
 
 
 class Roi(object):  # pragma: no cover
@@ -185,7 +185,7 @@ class RectangularROI(Roi):
     ymin, ymax :  float, optional
         y coordinates of lower and upper border
     theta : float, optional
-        Angle of clockwise rotation around center
+        Angle of anticlockwise rotation around center
     """
 
     def __init__(self, xmin=None, xmax=None, ymin=None, ymax=None, theta=None):
@@ -224,6 +224,7 @@ def move_to(self, x, y):
         self.ymax += dy
 
     def rotate_to(self, theta):
+        """ Rotate anticlockwise around center to position angle theta (radian) """
         self.theta = 0 if theta is None else theta
         if np.isclose(self.theta % np.pi, 0.0, atol=1e-9):
             self.rotation = np.identity(2)
@@ -317,6 +318,7 @@ def defined(self):
         return self.xmin is not None
 
     def to_polygon(self):
+        """ Returns vertices x, y, where each is an array of coordinates """
         if self.defined():
             if np.isclose(self.theta % np.pi, 0.0, atol=1e-9):
                 return (np.array([self.xmin, self.xmax, self.xmax, self.xmin, self.xmin]),
@@ -554,7 +556,7 @@ class EllipticalROI(Roi):
     radius_y :  float, optional
         Semiaxis along y axis
     theta : float, optional
-        Angle of clockwise rotation around (xc, yc)
+        Angle of anticlockwise rotation around (xc, yc)
     """
 
     def __init__(self, xc=None, yc=None, radius_x=None, radius_y=None, theta=None):
@@ -641,7 +643,7 @@ def get_center(self):
         return self.xc, self.yc
 
     def to_polygon(self):
-        """ Returns x, y, where each is a list of points """
+        """ Returns vertices x, y, where each is an array of coordinates """
         if not self.defined():
             return [], []
         theta = np.linspace(0, 2 * np.pi, num=20)
@@ -671,6 +673,7 @@ def move_to(self, xdelta, ydelta):
         self.yc += ydelta
 
     def rotate_to(self, theta):
+        """ Rotate anticlockwise around center to position angle theta (radian) """
         self.theta = 0 if theta is None else theta
         if np.isclose(self.theta % np.pi, 0.0, atol=1e-9):
             self.rotation = np.identity(2)
@@ -832,6 +835,7 @@ def move_to(self, xdelta, ydelta):
         self.vy = list(map(lambda y: y + ydelta, self.vy))
 
     def rotate_to(self, theta):
+        """ Rotate anticlockwise around center to position angle theta (radian) """
         theta = 0 if theta is None else theta
         dtheta = theta - self.theta
         if self.defined() and not np.isclose(dtheta % np.pi, 0.0, atol=1e-9):
diff --git a/glue/core/tests/test_roi.py b/glue/core/tests/test_roi.py
index b5b8b0e99..05a605691 100644
--- a/glue/core/tests/test_roi.py
+++ b/glue/core/tests/test_roi.py
@@ -80,12 +80,12 @@ def test_reset(self):
 
     def test_set_rotation(self):
         self.roi.update_limits(0, 0, 10, 4)
-        self.roi.rotate_by(-np.pi / 6)
+        self.roi.rotate_by(np.pi / 6)
         assert self.roi.contains(8, 6)
-        self.roi.rotate_by(-np.pi / 3)
+        self.roi.rotate_by(np.pi / 3)
         assert not self.roi.contains(9, 6)
         assert self.roi.contains(6, 6.9)
-        self.roi.rotate_to(np.pi / 3)
+        self.roi.rotate_to(-np.pi / 3)
         assert not self.roi.contains(5, 6)
         assert self.roi.contains(6, -3)
 
@@ -104,7 +104,7 @@ def test_to_polygon(self):
         x, y = self.roi.to_polygon()
         poly = PolygonalROI(vx=x, vy=y)
         assert not poly.contains(8, 3)
-        assert poly.contains(4, 6)
+        assert poly.contains(6, 6)
 
     def test_ndarray(self):
         self.roi.update_limits(0, 0, 10, 10)
@@ -353,8 +353,8 @@ def test_contains_on_undefined_contains_raises(self):
         with pytest.raises(UndefinedROI):
             self.roi_empty.contains(1, 1)
         assert self.roi.contains(1, 1)
-        assert self.roi_rotated.contains(0, 2.5)
-        assert not self.roi_rotated.contains(0, 2)
+        assert self.roi_rotated.contains(2, 2.5)
+        assert not self.roi_rotated.contains(2, 2)
 
     def test_set_center(self):
         assert self.roi.contains(0, 0)
@@ -376,12 +376,12 @@ def test_set_radius(self):
 
     def test_set_rotation(self):
         self.roi_rotated.rotate_to(0.55)
-        assert self.roi_rotated.contains(-1.5, 3.6)
+        assert self.roi_rotated.contains(-1.5, 0.4)
         self.roi_rotated.rotate_by(np.pi / 3)
-        assert self.roi_rotated.contains(1.2, 4.5)
+        assert self.roi_rotated.contains(0.8, 4.5)
         self.roi_rotated.rotate_to(np.pi / 3)
-        assert not self.roi_rotated.contains(1.2, 4.5)
-        assert self.roi_rotated.contains(0.5, 3.5)
+        assert not self.roi_rotated.contains(0.8, 4.5)
+        assert self.roi_rotated.contains(0.5, 0.5)
 
     def test_contains_many(self):
         x = [0, 0, 0, 0, 0]
@@ -398,8 +398,8 @@ def test_poly(self):
 
         x, y = self.roi_rotated.to_polygon()
         poly = PolygonalROI(vx=x, vy=y)
-        assert poly.contains(0, 2.5)
-        assert not poly.contains(0, 2)
+        assert poly.contains(2, 2.5)
+        assert not poly.contains(2, 2)
 
     def test_poly_undefined(self):
         x, y = self.roi_empty.to_polygon()
@@ -437,7 +437,7 @@ def test_serialize_rotated(self):
         new_roi = roundtrip_roi(self.roi_rotated)
         assert new_roi.radius_y == 0.4
         assert new_roi.theta == np.pi / 6
-        assert new_roi.contains(-1.5, 3.5)
+        assert new_roi.contains(-1.5, 0.5)
 
 
 class TestPolygon(object):
@@ -550,8 +550,8 @@ def test_serialization(self):
         new_roi = roundtrip_roi(self.roi)
         assert new_roi.theta == 0
         sqh = np.sqrt(0.5)
-        assert_almost_equal(new_roi.vx, np.array([-sqh, 0, sqh, 0]) + 0.5)
-        assert_almost_equal(new_roi.vy, np.array([0, sqh, 0, -sqh]) + 0.5)
+        assert_almost_equal(new_roi.vx, np.array([0, -sqh, 0, sqh]) + 0.5)
+        assert_almost_equal(new_roi.vy, np.array([-sqh, 0, sqh, 0]) + 0.5)
 
 
 class TestProjected3dROI(object):
diff --git a/glue/core/tests/test_subset.py b/glue/core/tests/test_subset.py
index dd9718af1..7548f18c3 100644
--- a/glue/core/tests/test_subset.py
+++ b/glue/core/tests/test_subset.py
@@ -873,7 +873,7 @@ def test_roi_subset_state_with_pretransform(self):
         assert_equal(data_clone.subsets[0].to_mask(), [0, 1, 0, 0])
 
         roi.move_to(3, 0)
-        roi.rotate_to(1.5)
+        roi.rotate_to(-1.5)
 
         subset = data_clone.new_subset()
         subset.subset_state = RoiSubsetState(xatt=data_clone.id['a'], yatt=data_clone.id['c'],
@@ -882,7 +882,7 @@ def test_roi_subset_state_with_pretransform(self):
 
     def test_roi_subset_state_rotated(self):
 
-        roi = RectangularROI(xmin=0, xmax=3, ymin=1.1, ymax=1.4, theta=np.pi / 3)
+        roi = RectangularROI(xmin=0, xmax=3, ymin=1.1, ymax=1.4, theta=-np.pi / 3)
 
         subset = self.data.new_subset()
         subset.subset_state = RoiSubsetState(xatt=self.data.id['a'], yatt=self.data.id['c'],
@@ -894,7 +894,7 @@ def test_roi_subset_state_rotated(self):
         assert_equal(data_clone.subsets[0].to_mask(), [0, 0, 0, 1])
 
         roi = RectangularROI(xmin=0, xmax=3, ymin=1.1, ymax=1.4)
-        roi.rotate_by(0.9)
+        roi.rotate_by(-0.9)
 
         subset = data_clone.new_subset()
         subset.subset_state = RoiSubsetState(xatt=data_clone.id['a'], yatt=data_clone.id['c'],
@@ -1089,7 +1089,7 @@ def test_roi_reduction():
     roi.rotate_to(np.pi / 4)
     state = RoiSubsetState(xatt=data4d.pixel_component_ids[0], yatt=data4d.pixel_component_ids[1], roi=roi)
     out = state.to_mask(data4d)
-    expected_slice = np.array([[0, 0, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]])
+    expected_slice = np.array([[0, 0, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]])
     assert_equal(out[:, :, 0, 0], expected_slice)
     assert_equal(out[:, :, 0, 1], expected_slice)
     assert_equal(out[:, :, 1, 0], expected_slice)

From 87ec7606c9922bad386f4a1472307ade44bbba42 Mon Sep 17 00:00:00 2001
From: Derek Homeier <dhomeie@gwdg.de>
Date: Wed, 29 Sep 2021 22:14:38 +0200
Subject: [PATCH 06/12] Add `Projected3dROI.roi_2d` rotation; move `rotate_by`
 to base `Roi` class

---
 glue/core/roi.py               | 28 +++++++++++-----------------
 glue/core/tests/test_roi.py    | 13 +++++++++++++
 glue/core/tests/test_subset.py |  1 -
 3 files changed, 24 insertions(+), 18 deletions(-)

diff --git a/glue/core/roi.py b/glue/core/roi.py
index 16adb6dbb..bff0be959 100644
--- a/glue/core/roi.py
+++ b/glue/core/roi.py
@@ -117,14 +117,6 @@ def move_to(self, x, y):
         """Translate the ROI to a center of (x, y)"""
         raise NotImplementedError()
 
-    def rotate_to(self, theta):
-        """Set rotation angle of ROI around center to theta"""
-        raise NotImplementedError()
-
-    def rotate_by(self, dtheta):
-        """Rotate the ROI around center by angle dtheta"""
-        raise NotImplementedError()
-
     def defined(self):
         """ Returns whether or not the subset is properly defined """
         raise NotImplementedError()
@@ -134,6 +126,14 @@ def to_polygon(self):
         as a polygon."""
         raise NotImplementedError()
 
+    def rotate_to(self, theta):
+        """Set rotation angle of ROI around center to theta"""
+        raise NotImplementedError()
+
+    def rotate_by(self, dtheta):
+        """Rotate the ROI around center by angle dtheta"""
+        self.rotate_to(getattr(self, 'theta', 0.0) + dtheta)
+
     def copy(self):
         """
         Return a clone of the ROI
@@ -232,9 +232,6 @@ def rotate_to(self, theta):
             self.rotation = rotation(self.theta)
         self.invrot = self.rotation * [[1, -1], [-1, 1]]
 
-    def rotate_by(self, dtheta):
-        self.rotate_to(self.theta + dtheta)
-
     def transpose(self, copy=True):
         if copy:
             new = self.copy()
@@ -681,9 +678,6 @@ def rotate_to(self, theta):
             self.rotation = rotation(self.theta)
         self.invrot = self.rotation * [[1, -1], [-1, 1]]
 
-    def rotate_by(self, dtheta):
-        self.rotate_to(self.theta + dtheta)
-
     def __gluestate__(self, context):
         return dict(xc=context.do(self.xc),
                     yc=context.do(self.yc),
@@ -843,9 +837,6 @@ def rotate_to(self, theta):
             self.vx, self.vy = rotation(dtheta) @ (dx, dy) + np.array(self.center()).reshape(2, 1)
         self.theta = theta
 
-    def rotate_by(self, dtheta):
-        self.rotate_to(self.theta + dtheta)
-
 
 class Projected3dROI(Roi):
     """"A region of interest defined in screen coordinates.
@@ -925,6 +916,9 @@ def to_polygon(self):
     def transformed(self, xfunc=None, yfunc=None):
         return self.roi_2d.transformed(xfunc, yfunc)
 
+    def rotate_to(self, theta):
+        return self.roi_2d.rotate_to(theta)
+
 
 class Path(VertexROIBase):
 
diff --git a/glue/core/tests/test_roi.py b/glue/core/tests/test_roi.py
index 05a605691..288795fef 100644
--- a/glue/core/tests/test_roi.py
+++ b/glue/core/tests/test_roi.py
@@ -585,6 +585,19 @@ def test_forward(self):
         assert roi.to_polygon() == roi_2d.to_polygon()
         assert roi.defined() == roi_2d.defined()
 
+    def test_rotate2d(self):
+        """ Test rotation of the 2d ROI """
+        roi_2d = PolygonalROI(vx=[1.5, 3.5, 3.5, 1.5], vy=[4.5, 4.5, 5.5, 5.5])
+        roi = Projected3dROI(roi_2d=roi_2d, projection_matrix=self.xyzw2yxzw)
+        assert roi.contains3d(self.x, self.y, self.z).tolist() == [True, False, False]
+        assert roi.contains3d(self.x_nd, self.y_nd, self.z_nd).tolist() == [[True, False], [False, True], [False, False]]
+
+        roi.rotate_by(np.pi / 2)
+        assert_almost_equal(roi.roi_2d.theta, np.pi / 2)
+        assert_almost_equal(roi_2d.theta, np.pi / 2)
+        assert roi.contains3d(self.x, self.y, self.z).tolist() == [True, True, False]
+        assert roi.contains3d(self.x_nd, self.y_nd, self.z_nd).tolist() == [[True, False], [True, True], [False, True]]
+
 
 class TestCategorical(object):
 
diff --git a/glue/core/tests/test_subset.py b/glue/core/tests/test_subset.py
index 7548f18c3..8f6f394e9 100644
--- a/glue/core/tests/test_subset.py
+++ b/glue/core/tests/test_subset.py
@@ -1094,4 +1094,3 @@ def test_roi_reduction():
     assert_equal(out[:, :, 0, 1], expected_slice)
     assert_equal(out[:, :, 1, 0], expected_slice)
     assert_equal(out[:, :, 1, 1], expected_slice)
-

From 7cce26e6b9e15f9d9f954cf5af91200cbeb95ae9 Mon Sep 17 00:00:00 2001
From: Derek Homeier <dhomeie@gwdg.de>
Date: Thu, 30 Sep 2021 23:45:52 +0200
Subject: [PATCH 07/12] Refactoring `rotation[_matrix_2d]`; default `theta` in
 gluestate; docstring updates

---
 glue/core/roi.py       | 130 +++++++++++++++++++----------------------
 glue/utils/geometry.py |  18 +++++-
 2 files changed, 76 insertions(+), 72 deletions(-)

diff --git a/glue/core/roi.py b/glue/core/roi.py
index bff0be959..3c05aa9b1 100644
--- a/glue/core/roi.py
+++ b/glue/core/roi.py
@@ -6,7 +6,7 @@
 
 from glue.core.component import CategoricalComponent
 from glue.core.exceptions import UndefinedROI
-from glue.utils import points_inside_poly, iterate_chunks
+from glue.utils import points_inside_poly, iterate_chunks, rotation_matrix_2d
 
 
 np.seterr(all='ignore')
@@ -55,12 +55,6 @@ def pixel_to_axes(axes, x, y):
     return axes.transAxes.inverted().transform(xy)
 
 
-def rotation(alpha):
-    """Return rotation matrix for angle alpha (increasing anticlockwise) around origin.
-    """
-    return np.array([[np.cos(alpha), -np.sin(alpha)], [np.sin(alpha), np.cos(alpha)]])
-
-
 class Roi(object):  # pragma: no cover
 
     """
@@ -127,11 +121,11 @@ def to_polygon(self):
         raise NotImplementedError()
 
     def rotate_to(self, theta):
-        """Set rotation angle of ROI around center to theta"""
+        """Set rotation angle of ROI around center to theta (radian)"""
         raise NotImplementedError()
 
     def rotate_by(self, dtheta):
-        """Rotate the ROI around center by angle dtheta"""
+        """Rotate the ROI around center by angle dtheta (radian)"""
         self.rotate_to(getattr(self, 'theta', 0.0) + dtheta)
 
     def copy(self):
@@ -180,12 +174,12 @@ class RectangularROI(Roi):
 
     Parameters
     ----------
-    xmin, xmax :  float, optional
-        x coordinates of left and right border
-    ymin, ymax :  float, optional
-        y coordinates of lower and upper border
+    xmin, xmax : float, optional
+        x coordinates of left and right border.
+    ymin, ymax : float, optional
+        y coordinates of lower and upper border.
     theta : float, optional
-        Angle of anticlockwise rotation around center
+        Angle of anticlockwise rotation around center in radian.
     """
 
     def __init__(self, xmin=None, xmax=None, ymin=None, ymax=None, theta=None):
@@ -195,11 +189,6 @@ def __init__(self, xmin=None, xmax=None, ymin=None, ymax=None, theta=None):
         self.ymin = ymin
         self.ymax = ymax
         self.theta = 0 if theta is None else theta
-        if np.isclose(self.theta % np.pi, 0.0, atol=1e-9):
-            self.rotation = np.identity(2)
-        else:
-            self.rotation = rotation(self.theta)
-        self.invrot = self.rotation * [[1, -1], [-1, 1]]
 
     def __str__(self):
         if self.defined() and self.theta == 0:
@@ -226,11 +215,6 @@ def move_to(self, x, y):
     def rotate_to(self, theta):
         """ Rotate anticlockwise around center to position angle theta (radian) """
         self.theta = 0 if theta is None else theta
-        if np.isclose(self.theta % np.pi, 0.0, atol=1e-9):
-            self.rotation = np.identity(2)
-        else:
-            self.rotation = rotation(self.theta)
-        self.invrot = self.rotation * [[1, -1], [-1, 1]]
 
     def transpose(self, copy=True):
         if copy:
@@ -253,16 +237,20 @@ def height(self):
 
     def contains(self, x, y):
         """
-        Test whether a set of (x,y) points falls within
-        the region of interest
+        Test whether a set of (x,y) points falls within the region of interest
 
-        :param x: A scalar or iterable of x coordinates
-        :param y: A scalar or iterable of y coordinates
+        Parameters
+        ----------
+        x : float or array-like
+            x coordinate(s) of point(s).
+        y : float or array-like
+            y coordinate(s) of point(s).
 
-        *Returns*
+        Returns
+        -------
 
-            A list of True/False values, for whether each (x,y)
-            point falls within the ROI
+        inside : bool or `numpy.ndarray`
+            An iterable of True/False values, for whether each (x,y) point falls within the ROI
         """
 
         if not self.defined():
@@ -289,7 +277,7 @@ def contains(self, x, y):
         x = x[keep] - xc
         y = y[keep] - yc
         shape = (2,) + x.shape
-        x, y = (self.invrot @ [x.flatten(), y.flatten()]).reshape(shape)
+        x, y = (rotation_matrix_2d(-self.theta) @ [x.flatten(), y.flatten()]).reshape(shape)
         inside[keep] = (abs(x) <= self.width() / 2) & (abs(y) <= self.height() / 2)
         return inside
 
@@ -323,7 +311,8 @@ def to_polygon(self):
             else:
                 corners = (np.array([-1, 1, 1, -1, -1]) * self.width() / 2,
                            np.array([-1, -1, 1, 1, -1]) * self.height() / 2)
-            return tuple((self.rotation @ corners) + np.array(self.center()).reshape((2, 1)))
+            return tuple((rotation_matrix_2d(self.theta) @ corners) +
+                         np.array(self.center()).reshape((2, 1)))
         else:
             return [], []
 
@@ -345,7 +334,7 @@ def __gluestate__(self, context):
     def __setgluestate__(cls, rec, context):
         return cls(xmin=context.object(rec['xmin']), xmax=context.object(rec['xmax']),
                    ymin=context.object(rec['ymin']), ymax=context.object(rec['ymax']),
-                   theta=context.object(rec['theta']))
+                   theta=context.object(rec.get('theta', 0)))
 
 
 class RangeROI(Roi):
@@ -544,16 +533,16 @@ class EllipticalROI(Roi):
 
     Parameters
     ----------
-    xc :  float, optional
-        x coordinate of center
-    yc :  float, optional
-        y coordinate of center
-    radius_x :  float, optional
-        Semiaxis along x axis
-    radius_y :  float, optional
-        Semiaxis along y axis
+    xc : float, optional
+        x coordinate of center.
+    yc : float, optional
+        y coordinate of center.
+    radius_x : float, optional
+        Semiaxis along x axis.
+    radius_y : float, optional
+        Semiaxis along y axis.
     theta : float, optional
-        Angle of anticlockwise rotation around (xc, yc)
+        Angle of anticlockwise rotation around (xc, yc) in radian.
     """
 
     def __init__(self, xc=None, yc=None, radius_x=None, radius_y=None, theta=None):
@@ -563,11 +552,6 @@ def __init__(self, xc=None, yc=None, radius_x=None, radius_y=None, theta=None):
         self.radius_x = radius_x
         self.radius_y = radius_y
         self.theta = 0 if theta is None else theta
-        if np.isclose(self.theta % np.pi, 0.0, atol=1e-9):
-            self.rotation = np.identity(2)
-        else:
-            self.rotation = rotation(self.theta)
-        self.invrot = self.rotation * [[1, -1], [-1, 1]]
 
     def __str__(self):
         if self.defined():
@@ -582,14 +566,18 @@ def contains(self, x, y):
         Test whether a set of (x,y) points falls within
         the region of interest
 
-        :param x: A scalar or iterable of x coordinates
-        :param y: A scalar or iterable of y coordinates
-
-        *Returns*
+        Parameters
+        ----------
+        x : float or array-like
+            x coordinate(s) of point(s).
+        y : float or array-like
+            y coordinate(s) of point(s).
 
-           A list of True/False values, for whether each (x,y)
-           point falls within the ROI
+        Returns
+        -------
 
+        inside : bool or `numpy.ndarray`
+            An iterable of True/False values, for whether each (x,y) point falls within the ROI
         """
         if not self.defined():
             raise UndefinedROI
@@ -616,7 +604,7 @@ def contains(self, x, y):
             x = x[keep] - self.xc
             y = y[keep] - self.yc
             shape = (2,) + x.shape
-            x, y = (self.invrot @ [x.flatten(), y.flatten()]).reshape(shape)
+            x, y = (rotation_matrix_2d(-self.theta) @ [x.flatten(), y.flatten()]).reshape(shape)
             inside[keep] = ((x ** 2 / self.radius_x ** 2 + y ** 2 / self.radius_y ** 2) < 1.)
             return inside
 
@@ -646,7 +634,7 @@ def to_polygon(self):
         theta = np.linspace(0, 2 * np.pi, num=20)
         x = self.radius_x * np.cos(theta)
         y = self.radius_y * np.sin(theta)
-        x, y = self.rotation @ (x, y)
+        x, y = rotation_matrix_2d(self.theta) @ (x, y)
         return x + self.xc, y + self.yc
 
     def bounds(self):
@@ -672,11 +660,6 @@ def move_to(self, xdelta, ydelta):
     def rotate_to(self, theta):
         """ Rotate anticlockwise around center to position angle theta (radian) """
         self.theta = 0 if theta is None else theta
-        if np.isclose(self.theta % np.pi, 0.0, atol=1e-9):
-            self.rotation = np.identity(2)
-        else:
-            self.rotation = rotation(self.theta)
-        self.invrot = self.rotation * [[1, -1], [-1, 1]]
 
     def __gluestate__(self, context):
         return dict(xc=context.do(self.xc),
@@ -688,7 +671,7 @@ def __gluestate__(self, context):
     @classmethod
     def __setgluestate__(cls, rec, context):
         return cls(xc=rec['xc'], yc=rec['yc'],
-                   radius_x=rec['radius_x'], radius_y=rec['radius_y'], theta=rec['theta'])
+                   radius_x=rec['radius_x'], radius_y=rec['radius_y'], theta=rec.get('theta', 0))
 
 
 class VertexROIBase(Roi):
@@ -792,18 +775,22 @@ def __str__(self):
 
     def contains(self, x, y):
         """
-        Test whether a set of (x,y) points falls within
-        the region of interest
-
-        :param x: A list of x points
-        :param y: A list of y points
+        Test whether a set of (x,y) points falls within the region of interest
 
-        *Returns*
+        Parameters
+        ----------
+        x : float or array-like
+            x coordinate(s) of point(s).
+        y : float or array-like
+            y coordinate(s) of point(s).
 
-           A list of True/False values, for whether each (x,y)
-           point falls within the ROI
+        Returns
+        -------
 
+        inside : bool or `numpy.ndarray`
+            An iterable of True/False values, for whether each (x,y) point falls within the ROI
         """
+
         if not self.defined():
             raise UndefinedROI
         if not isinstance(x, np.ndarray):
@@ -834,7 +821,8 @@ def rotate_to(self, theta):
         dtheta = theta - self.theta
         if self.defined() and not np.isclose(dtheta % np.pi, 0.0, atol=1e-9):
             dx, dy = np.array([self.vx, self.vy]) - np.array(self.center()).reshape(2, 1)
-            self.vx, self.vy = rotation(dtheta) @ (dx, dy) + np.array(self.center()).reshape(2, 1)
+            self.vx, self.vy = (rotation_matrix_2d(dtheta) @ (dx, dy) +
+                                np.array(self.center()).reshape(2, 1))
         self.theta = theta
 
 
diff --git a/glue/utils/geometry.py b/glue/utils/geometry.py
index 07c9de801..26330f8cd 100644
--- a/glue/utils/geometry.py
+++ b/glue/utils/geometry.py
@@ -2,7 +2,23 @@
 
 from glue.utils import unbroadcast, broadcast_to
 
-__all__ = ['points_inside_poly', 'polygon_line_intersections', 'floodfill']
+__all__ = ['points_inside_poly', 'polygon_line_intersections', 'floodfill', 'rotation_matrix_2d']
+
+
+def rotation_matrix_2d(alpha):
+    """
+    Return rotation matrix for angle alpha around origin.
+
+    Parameters
+    ----------
+    alpha : float
+        Rotation angle in radian, increasing for anticlockwise rotation.
+    """
+    if np.asarray(alpha).ndim > 0:
+        # In principle this works on an array as well; would have to return matrix.T then
+        raise ValueError("Only scalar input for angle accepted")
+
+    return np.array([[np.cos(alpha), -np.sin(alpha)], [np.sin(alpha), np.cos(alpha)]])
 
 
 def points_inside_poly(x, y, vx, vy):

From 25168027d2dba55513d777cba183651f699b87e2 Mon Sep 17 00:00:00 2001
From: Derek Homeier <dhomeie@gwdg.de>
Date: Thu, 30 Sep 2021 23:49:43 +0200
Subject: [PATCH 08/12] Docs reorganisation and numpydoc updates; changelog
 entry

---
 CHANGES.md       |   4 +
 glue/core/roi.py | 346 +++++++++++++++++++++++------------------------
 2 files changed, 177 insertions(+), 173 deletions(-)

diff --git a/CHANGES.md b/CHANGES.md
index 78058d3ed..e21905e6f 100644
--- a/CHANGES.md
+++ b/CHANGES.md
@@ -15,6 +15,10 @@ v1.4.0 (2022-05-31)
   is defined, resulting in a broadcast error in arrays large enough to need
   chunking. [#2302]
 
+* Added rotation angle ``theta`` as a property to regions of interest,
+  to be set on instantiation or modified using the ``rotate_to`` and
+  ``rotate_by`` methods. [#2235]
+
 v1.3.0 (2022-04-22)
 -------------------
 
diff --git a/glue/core/roi.py b/glue/core/roi.py
index 3c05aa9b1..08a0b0888 100644
--- a/glue/core/roi.py
+++ b/glue/core/roi.py
@@ -23,9 +23,7 @@
 
 
 def aspect_ratio(axes):
-    """ Returns the pixel height / width of a box that spans 1
-    data unit in x and y
-    """
+    """Returns the pixel height / width of a box that spans 1 data unit in `x` and `y`"""
     width = axes.get_position().width * axes.figure.get_figwidth()
     height = axes.get_position().height * axes.figure.get_figheight()
     xmin, xmax = axes.get_xlim()
@@ -60,41 +58,54 @@ class Roi(object):  # pragma: no cover
     """
     A geometrical 2D region of interest.
 
-    Glue uses Roi's to represent user-drawn regions on plots. There
+    Glue uses ROIs to represent user-drawn regions on plots. There
     are many specific sub-classes of Roi, but they all have a ``contains``
     method to test whether a collection of 2D points lies inside the region.
     """
 
     def contains(self, x, y):
-        """Return true/false for each x/y pair.
+        """
+        Test which of a set of (`x`, `y`) points fall within the region of interest
 
-        :param x: Array of X locations
-        :param y: Array of Y locations
+        Parameters
+        ----------
+        x : float or array-like
+            `x` coordinate(s) of point(s).
+        y : float or array-like
+            `y` coordinate(s) of point(s).
 
-        :returns: A Boolean array, where each element is True
-                  if the corresponding (x,y) tuple is inside the Roi.
+        Returns
+        -------
 
-        :raises: UndefinedROI exception if not defined
+        inside : bool or `numpy.ndarray`
+            An boolean iterable, where each element is `True` if the corresponding
+            (`x`, `y`) tuple is inside the Roi.
+
+        Raises
+        ------
+        UndefinedROI
+            if not defined
         """
         raise NotImplementedError()
 
     def contains3d(self, x, y, z):
-        """Return true/false for each x/y/z pair.
+        """
+        Test which of a set of projected (`x`, `y`, `z`) points fall within the region of interest
 
         Parameters
         ----------
         x : :class:`numpy.ndarray`
-            Array of x locations
+            Array of `x` locations
         y : :class:`numpy.ndarray`
-            Array of y locations
+            Array of `y` locations
         z : :class:`numpy.ndarray`
-            Array of z locations
+            Array of `z` locations
 
         Returns
         -------
         :class:`numpy.ndarray`
             A boolean array, where each element is `True` if the corresponding
-            (x,y,z) tuple is inside the Roi.
+            (`x`, `y`, `z`) tuple is projected inside the linked 2-d Roi.
 
         Raises
         ------
@@ -104,40 +115,38 @@ def contains3d(self, x, y, z):
         raise NotImplementedError()
 
     def center(self):
-        """Return the (x,y) coordinates of the ROI center"""
+        """Return the (`x`, `y`) coordinates of the ROI center"""
         raise NotImplementedError()
 
     def move_to(self, x, y):
-        """Translate the ROI to a center of (x, y)"""
+        """Translate the ROI to a center of (`x`, `y`)"""
         raise NotImplementedError()
 
     def defined(self):
-        """ Returns whether or not the subset is properly defined """
+        """Returns `True` if the ROI is defined"""
         raise NotImplementedError()
 
     def to_polygon(self):
-        """ Returns a tuple of x and y points, approximating the ROI
-        as a polygon."""
+        """
+        Returns vertices `vx`, `vy` of a polygon approximating the Roi,
+        where each is an array of vertex coordinates in `x` and `y`.
+        """
         raise NotImplementedError()
 
     def rotate_to(self, theta):
-        """Set rotation angle of ROI around center to theta (radian)"""
+        """Rotate anticlockwise around center to position angle theta (radian)"""
         raise NotImplementedError()
 
     def rotate_by(self, dtheta):
-        """Rotate the ROI around center by angle dtheta (radian)"""
+        """Rotate the Roi around center by angle dtheta (radian)"""
         self.rotate_to(getattr(self, 'theta', 0.0) + dtheta)
 
     def copy(self):
-        """
-        Return a clone of the ROI
-        """
+        """Return a clone of the Roi"""
         return copy.copy(self)
 
     def transformed(self, xfunc=None, yfunc=None):
-        """
-        A transformed version of the ROI
-        """
+        """A transformed version of the Roi"""
         raise NotImplementedError()
 
 
@@ -175,11 +184,17 @@ class RectangularROI(Roi):
     Parameters
     ----------
     xmin, xmax : float, optional
-        x coordinates of left and right border.
+        `x` coordinates of left and right edge.
     ymin, ymax : float, optional
-        y coordinates of lower and upper border.
+        `y` coordinates of lower and upper edge.
     theta : float, optional
         Angle of anticlockwise rotation around center in radian.
+
+    Notes
+    -----
+        The input and ``update_limits()`` parameters specify the `x` and `y` edge
+        positions *before* any rotation is applied; the size always remains
+        :math:`width` = `xmax` - `xmin`; `height` = `ymax` - `ymin`.
     """
 
     def __init__(self, xmin=None, xmax=None, ymin=None, ymax=None, theta=None):
@@ -213,7 +228,6 @@ def move_to(self, x, y):
         self.ymax += dy
 
     def rotate_to(self, theta):
-        """ Rotate anticlockwise around center to position angle theta (radian) """
         self.theta = 0 if theta is None else theta
 
     def transpose(self, copy=True):
@@ -236,23 +250,6 @@ def height(self):
         return self.ymax - self.ymin
 
     def contains(self, x, y):
-        """
-        Test whether a set of (x,y) points falls within the region of interest
-
-        Parameters
-        ----------
-        x : float or array-like
-            x coordinate(s) of point(s).
-        y : float or array-like
-            y coordinate(s) of point(s).
-
-        Returns
-        -------
-
-        inside : bool or `numpy.ndarray`
-            An iterable of True/False values, for whether each (x,y) point falls within the ROI
-        """
-
         if not self.defined():
             raise UndefinedROI
 
@@ -282,18 +279,14 @@ def contains(self, x, y):
         return inside
 
     def update_limits(self, xmin, ymin, xmax, ymax):
-        """
-        Update the limits of the rectangle
-        """
+        """Update the limits (edge positions) of the rectangle before rotation"""
         self.xmin = min(xmin, xmax)
         self.xmax = max(xmin, xmax)
         self.ymin = min(ymin, ymax)
         self.ymax = max(ymin, ymax)
 
     def reset(self):
-        """
-        Reset the rectangular region.
-        """
+        """Reset the rectangular region"""
         self.xmin = None
         self.xmax = None
         self.ymin = None
@@ -303,7 +296,10 @@ def defined(self):
         return self.xmin is not None
 
     def to_polygon(self):
-        """ Returns vertices x, y, where each is an array of coordinates """
+        """
+        Returns vertices `vx`, `vy` of the rectangular region represented as a polygon,
+        where each is an array of vertex coordinates in `x` and `y`.
+        """
         if self.defined():
             if np.isclose(self.theta % np.pi, 0.0, atol=1e-9):
                 return (np.array([self.xmin, self.xmax, self.xmax, self.xmin, self.xmin]),
@@ -338,9 +334,19 @@ def __setgluestate__(cls, rec, context):
 
 
 class RangeROI(Roi):
+    """
+    A region of interest representing all points within a range in either `x` or `y`.
 
+    Parameters
+    ----------
+    orientation : str
+        One of 'x' or 'y', setting the axis on which to apply the range.
+    min : float, optional
+        Start value of the range.
+    max : float, optional
+        End value of the range.
+    """
     def __init__(self, orientation, min=None, max=None):
-        """:param orientation: 'x' or 'y'. Sets which axis to range"""
         super(RangeROI, self).__init__()
 
         self.min = min
@@ -438,6 +444,15 @@ class CircularROI(Roi):
 
     """
     A 2D circular region of interest.
+
+    Parameters
+    ----------
+    xc : float, optional
+        `x` coordinate of center.
+    yc : float, optional
+        `y` coordinate of center.
+    radius : float, optional
+        Radius of the circle.
     """
 
     def __init__(self, xc=None, yc=None, radius=None):
@@ -447,19 +462,6 @@ def __init__(self, xc=None, yc=None, radius=None):
         self.radius = radius
 
     def contains(self, x, y):
-        """
-        Test whether a set of (x,y) points falls within
-        the region of interest
-
-        :param x: A scalar or iterable of x coordinates
-        :param y: A scalar or iterable of y coordinates
-
-        *Returns*
-
-           A list of True/False values, for whether each (x,y)
-           point falls within the ROI
-
-        """
         if not self.defined():
             raise UndefinedROI
 
@@ -470,16 +472,12 @@ def contains(self, x, y):
         return (x - self.xc) ** 2 + (y - self.yc) ** 2 < self.radius ** 2
 
     def set_center(self, x, y):
-        """
-        Set the center of the circular region
-        """
+        """Set the center of the circular region"""
         self.xc = x
         self.yc = y
 
     def set_radius(self, radius):
-        """
-        Set the radius of the circular region
-        """
+        """Set the radius of the circular region"""
         self.radius = radius
 
     def get_center(self):
@@ -489,20 +487,16 @@ def get_radius(self):
         return self.radius
 
     def reset(self):
-        """
-        Reset the rectangular region.
-        """
+        """Reset the circular region"""
         self.xc = None
         self.yc = None
         self.radius = 0.
 
     def defined(self):
-        """ Returns True if the ROI is defined """
         return self.xc is not None and \
             self.yc is not None and self.radius is not None
 
     def to_polygon(self):
-        """ Returns x, y, where each is a list of points """
         if not self.defined():
             return [], []
         theta = np.linspace(0, 2 * np.pi, num=20)
@@ -529,20 +523,25 @@ def __setgluestate__(cls, rec, context):
 
 class EllipticalROI(Roi):
     """
-    A 2D elliptical region of interest at (xc, yc) with semimajor/minor axes radius_[xy].
+    A 2D elliptical region of interest with semimajor/minor axes `radius_[xy]`.
 
     Parameters
     ----------
     xc : float, optional
-        x coordinate of center.
+        `x` coordinate of center.
     yc : float, optional
-        y coordinate of center.
+        `y` coordinate of center.
     radius_x : float, optional
-        Semiaxis along x axis.
+        Semiaxis along `x` axis.
     radius_y : float, optional
-        Semiaxis along y axis.
+        Semiaxis along `y` axis.
     theta : float, optional
-        Angle of anticlockwise rotation around (xc, yc) in radian.
+        Angle of anticlockwise rotation around (`xc`, `yc`) in radian.
+
+    Notes
+    -----
+        The `radius_x`, `radius_y` properties refer to the semiaxes along the `x` and `y`
+        axes *before* any rotation is applied.
     """
 
     def __init__(self, xc=None, yc=None, radius_x=None, radius_y=None, theta=None):
@@ -562,23 +561,6 @@ def __str__(self):
             return "Undefined Elliptical ROI"
 
     def contains(self, x, y):
-        """
-        Test whether a set of (x,y) points falls within
-        the region of interest
-
-        Parameters
-        ----------
-        x : float or array-like
-            x coordinate(s) of point(s).
-        y : float or array-like
-            y coordinate(s) of point(s).
-
-        Returns
-        -------
-
-        inside : bool or `numpy.ndarray`
-            An iterable of True/False values, for whether each (x,y) point falls within the ROI
-        """
         if not self.defined():
             raise UndefinedROI
 
@@ -609,16 +591,13 @@ def contains(self, x, y):
             return inside
 
     def reset(self):
-        """
-        Reset the rectangular region.
-        """
+        """Reset the rectangular region"""
         self.xc = None
         self.yc = None
         self.radius_x = 0.
         self.radius_y = 0.
 
     def defined(self):
-        """ Returns True if the ROI is defined """
         return (self.xc is not None and
                 self.yc is not None and
                 self.radius_x is not None and
@@ -628,7 +607,6 @@ def get_center(self):
         return self.xc, self.yc
 
     def to_polygon(self):
-        """ Returns vertices x, y, where each is an array of coordinates """
         if not self.defined():
             return [], []
         theta = np.linspace(0, 2 * np.pi, num=20)
@@ -638,7 +616,7 @@ def to_polygon(self):
         return x + self.xc, y + self.yc
 
     def bounds(self):
-        """ Returns (conservatively estimated) boundary values in x and y """
+        """Returns (conservatively estimated) boundary values in `x` and `y`"""
         if self.theta is None or np.isclose(self.theta % (np.pi), 0.0, atol=1e-9):
             return [[self.xc - self.radius_x, self.xc + self.radius_x],
                     [self.yc - self.radius_y, self.yc + self.radius_y]]
@@ -658,7 +636,6 @@ def move_to(self, xdelta, ydelta):
         self.yc += ydelta
 
     def rotate_to(self, theta):
-        """ Rotate anticlockwise around center to position angle theta (radian) """
         self.theta = 0 if theta is None else theta
 
     def __gluestate__(self, context):
@@ -676,13 +653,23 @@ def __setgluestate__(cls, rec, context):
 
 class VertexROIBase(Roi):
 
+    """
+    Class representing a set of vertices e.g. to define a 2D polygon.
+
+    Parameters
+    ----------
+    vx : float or array-like, optional
+        Initial `x` vertices.
+    vy : float or array-like, optional
+        Initial `y` vertices.
+
+    Notes
+    -----
+    This class only comprises the collection of vertices, but does not
+    provide a ``contains`` method.
+    """
+
     def __init__(self, vx=None, vy=None):
-        """
-        :param vx: initial x vertices
-        :type vx: list
-        :param vy: initial y vertices
-        :type vy: list
-        """
         super(VertexROIBase, self).__init__()
         self.vx = vx
         self.vy = vy
@@ -699,18 +686,20 @@ def transformed(self, xfunc=None, yfunc=None):
 
     def add_point(self, x, y):
         """
-        Add another vertex to the ROI
+        Add another vertex to the ROI.
 
-        :param x: The x coordinate
-        :param y: The y coordinate
+        Parameters
+        ----------
+        x : float
+            The `x` coordinate of the point to add.
+        y : float
+            The `y` coordinate of the point to add.
         """
         self.vx.append(x)
         self.vy.append(y)
 
     def reset(self):
-        """
-        Reset the vertex list.
-        """
+        """Reset the vertex lists"""
         self.vx = []
         self.vy = []
 
@@ -720,15 +709,18 @@ def replace_last_point(self, x, y):
             self.vy[-1] = y
 
     def remove_point(self, x, y, thresh=None):
-        """Remove the vertex closest to a reference (xy) point
-
-        :param x: The x coordinate of the reference point
-        :param y: The y coordinate of the reference point
-
-        :param thresh: An optional threshold. If present, the vertex
-                closest to (x,y) will only be removed if the distance
-                is less than thresh
+        """
+        Remove the vertex closest to a reference (`x`, `y`) point.
 
+        Parameters
+        ----------
+        x : float
+            The `x` coordinate of the reference point.
+        y : float
+            The `y` coordinate of the reference point.
+        thresh : float, optional
+            Threshold. If set, the vertex closest to (`x`, `y`) will only be removed
+            if the distance is less than `thresh`.
         """
         if len(self.vx) == 0:
             return
@@ -763,7 +755,14 @@ def __setgluestate__(cls, rec, context):
 class PolygonalROI(VertexROIBase):
 
     """
-    A class to define 2D polygonal regions-of-interest
+    A class to define 2D polygonal regions of interest.
+
+    Parameters
+    ----------
+    vx : float or array-like, optional
+        Initial `x` vertices.
+    vy : float or array-like, optional
+        Initial `y` vertices.
     """
 
     def __str__(self):
@@ -774,23 +773,6 @@ def __str__(self):
         return result
 
     def contains(self, x, y):
-        """
-        Test whether a set of (x,y) points falls within the region of interest
-
-        Parameters
-        ----------
-        x : float or array-like
-            x coordinate(s) of point(s).
-        y : float or array-like
-            y coordinate(s) of point(s).
-
-        Returns
-        -------
-
-        inside : bool or `numpy.ndarray`
-            An iterable of True/False values, for whether each (x,y) point falls within the ROI
-        """
-
         if not self.defined():
             raise UndefinedROI
         if not isinstance(x, np.ndarray):
@@ -816,7 +798,6 @@ def move_to(self, xdelta, ydelta):
         self.vy = list(map(lambda y: y + ydelta, self.vy))
 
     def rotate_to(self, theta):
-        """ Rotate anticlockwise around center to position angle theta (radian) """
         theta = 0 if theta is None else theta
         dtheta = theta - self.theta
         if self.defined() and not np.isclose(dtheta % np.pi, 0.0, atol=1e-9):
@@ -827,12 +808,20 @@ def rotate_to(self, theta):
 
 
 class Projected3dROI(Roi):
-    """"A region of interest defined in screen coordinates.
+    """
+    A region of interest defined in screen coordinates.
 
     The screen coordinates are defined by the projection matrix.
-    The projection matrix converts homogeneous coordinates (x, y, z, w), where
-    w is implicitly 1, to homogeneous screen coordinates (usually the product
-    of the world and projection matrix).
+    The projection matrix converts homogeneous coordinates (`x`, `y`, `z`, `w`),
+    where `w` is implicitly 1, to homogeneous screen coordinates (usually the
+    tensor product of the world coordinate vectors and the projection matrix).
+
+    Parameters
+    ----------
+    2d_roi : `~glue.core.roi.Roi`, optional
+        If specified, this ROI will be used in screen coordinate space.
+    projection_matrix : `~numpy.ndarray`, optional
+        Projection matrix defining the mapping from world onto screen coordinates.
     """
 
     def __init__(self, roi_2d=None, projection_matrix=None):
@@ -841,10 +830,6 @@ def __init__(self, roi_2d=None, projection_matrix=None):
         self.projection_matrix = np.asarray(projection_matrix)
 
     def contains3d(self, x, y, z):
-        """
-        Test whether the projected coordinates are contained in the 2d ROI.
-        """
-
         if not self.defined():
             raise UndefinedROI
 
@@ -1694,8 +1679,14 @@ def _categorical_helper(self, indata):
         """
         A helper function to do the rigamaroll of getting categorical data.
 
-        :param indata: Any type of input data
-        :return: The best guess at the categorical data associated with indata
+        Parameters
+        ----------
+        indata : object
+            Any type of input data
+
+        Returns
+        -------
+            The best guess at the categorical data associated with indata.
         """
 
         try:
@@ -1708,16 +1699,18 @@ def _categorical_helper(self, indata):
 
     def contains(self, x, y):
         """
-        Test whether a set categorical elements fall within
-        the region of interest
+        Test whether a set categorical elements fall within the region of interest.
 
-        :param x: Any array-like object of categories
-                 (includes CategoricalComponenets)
-        :param y: Unused but required for compatibility
-
-        *Returns*
+        Parameters
+        ----------
+        x : array-like
+            An array-like object of categories (includes `CategoricalComponents`).
+        y : object or None
+            Unused but required for compatibility
 
-           A list of True/False values, for whether each x value falls
+        Returns
+        -------
+           A list of True/False values, for whether each `x` value falls
            within the ROI
 
         """
@@ -1733,7 +1726,6 @@ def update_categories(self, categories):
         self.categories = np.unique(self._categorical_helper(categories))
 
     def defined(self):
-        """ Returns True if the ROI is defined """
         return self.categories is not None
 
     def reset(self):
@@ -1742,12 +1734,20 @@ def reset(self):
     @staticmethod
     def from_range(categories, lo, hi):
         """
-        Utility function to help construct the Roi from a range.
+        Utility function to help construct the ROI from a range.
+
+        Parameters
+        ----------
+        categories : object
+            Anything understood by ``._categorical_helper`` ... array, list or component.
+        lo : int or float
+            Lower bound of the range (rounded up to next integer)
+        hi : int or float
+            Upper bound of the range (rounded up to next integer)
 
-        :param cat_comp: Anything understood by ._categorical_helper ... array, list or component
-        :param lo: lower bound of the range
-        :param hi: upper bound of the range
-        :return: CategoricalROI object
+        Returns
+        -------
+            `CategoricalROI` object
         """
 
         # Convert lo and hi to integers. Note that if lo or hi are negative,

From 22f230215bbfde43a0b2defd5f4233e6e7b0aa78 Mon Sep 17 00:00:00 2001
From: Derek Homeier <dhomeie@gwdg.de>
Date: Fri, 1 Oct 2021 14:37:09 +0200
Subject: [PATCH 09/12] DOC: further docstring refinements [ci skip]

---
 glue/core/roi.py | 24 +++++++++++++-----------
 1 file changed, 13 insertions(+), 11 deletions(-)

diff --git a/glue/core/roi.py b/glue/core/roi.py
index 08a0b0888..be00487a0 100644
--- a/glue/core/roi.py
+++ b/glue/core/roi.py
@@ -61,11 +61,13 @@ class Roi(object):  # pragma: no cover
     Glue uses ROIs to represent user-drawn regions on plots. There
     are many specific sub-classes of Roi, but they all have a ``contains``
     method to test whether a collection of 2D points lies inside the region.
+    ROI bounds are generally designed as being exclusive (that is, points
+    situated exactly on a border are considered to lie outside the region).
     """
 
     def contains(self, x, y):
         """
-        Test which of a set of (`x`, `y`) points fall within the region of interest
+        Test which of a set of (`x`, `y`) points fall within the region of interest.
 
         Parameters
         ----------
@@ -76,41 +78,41 @@ def contains(self, x, y):
 
         Returns
         -------
-
-        inside : bool or `numpy.ndarray`
+        inside : bool or `~numpy.ndarray`
             An boolean iterable, where each element is `True` if the corresponding
             (`x`, `y`) tuple is inside the Roi.
 
         Raises
         ------
         UndefinedROI
-            if not defined
+            If not defined.
         """
         raise NotImplementedError()
 
     def contains3d(self, x, y, z):
         """
-        Test which of a set of projected (`x`, `y`, `z`) points fall within the region of interest
+        Test which of a set of projected (`x`, `y`, `z`) points fall within the
+        linked 2D region of interest.
 
         Parameters
         ----------
-        x : :class:`numpy.ndarray`
+        x : :class:`~numpy.ndarray`
             Array of `x` locations
-        y : :class:`numpy.ndarray`
+        y : :class:`~numpy.ndarray`
             Array of `y` locations
-        z : :class:`numpy.ndarray`
+        z : :class:`~numpy.ndarray`
             Array of `z` locations
 
         Returns
         -------
-        :class:`numpy.ndarray`
+        :class:`~numpy.ndarray`
             A boolean array, where each element is `True` if the corresponding
-            (`x`, `y`, `z`) tuple is projected inside the linked 2-d Roi.
+            (`x`, `y`, `z`) tuple is projected inside the associated 2D Roi.
 
         Raises
         ------
         UndefinedROI
-            if not defined
+            If not defined.
         """
         raise NotImplementedError()
 

From b6d9d76fbdfdd555984561ae7be44737768512f0 Mon Sep 17 00:00:00 2001
From: Derek Homeier <dhomeie@gwdg.de>
Date: Tue, 1 Feb 2022 00:00:32 +0100
Subject: [PATCH 10/12] Introduce polygon centroid as default rotation point

---
 glue/core/roi.py            | 88 +++++++++++++++++++++++++++++++------
 glue/core/tests/test_roi.py | 35 ++++++++++++++-
 setup.cfg                   |  4 +-
 3 files changed, 110 insertions(+), 17 deletions(-)

diff --git a/glue/core/roi.py b/glue/core/roi.py
index be00487a0..7a01a32f1 100644
--- a/glue/core/roi.py
+++ b/glue/core/roi.py
@@ -139,9 +139,9 @@ def rotate_to(self, theta):
         """Rotate anticlockwise around center to position angle theta (radian)"""
         raise NotImplementedError()
 
-    def rotate_by(self, dtheta):
+    def rotate_by(self, dtheta, **kwargs):
         """Rotate the Roi around center by angle dtheta (radian)"""
-        self.rotate_to(getattr(self, 'theta', 0.0) + dtheta)
+        self.rotate_to(getattr(self, 'theta', 0.0) + dtheta, **kwargs)
 
     def copy(self):
         """Return a clone of the Roi"""
@@ -673,12 +673,8 @@ class VertexROIBase(Roi):
 
     def __init__(self, vx=None, vy=None):
         super(VertexROIBase, self).__init__()
-        self.vx = vx
-        self.vy = vy
-        if self.vx is None:
-            self.vx = []
-        if self.vy is None:
-            self.vy = []
+        self.vx = [] if vx is None else list(vx)
+        self.vy = [] if vy is None else list(vy)
         self.theta = 0
 
     def transformed(self, xfunc=None, yfunc=None):
@@ -785,27 +781,91 @@ def contains(self, x, y):
         result = points_inside_poly(x, y, self.vx, self.vy)
         return result
 
-    # There are several possible definitions of the centre; this is easiest to calculate
-    # (and invariant under rotation?) - do not include starting vertex twice!
-    def center(self):
+    # There are several possible definitions of the centre; `mean()` is
+    # easiest to calculate, but not robust against adding vertices.
+    def mean(self):
+        """Return arithmetic mean (of vertex positions) of polygon."""
+
         if not self.defined():
             raise UndefinedROI
+        # Do not include starting vertex twice!
         if self.vx[-1] == self.vx[0] and self.vy[:-1] == self.vy[0]:
             return np.mean(self.vx[:-1]), np.mean(self.vy[:-1])
         else:
             return np.mean(self.vx), np.mean(self.vy)
 
+    def area(self, signed=False):
+        """
+        Return area of polygon using the shoelace formula.
+
+        Parameters
+        ----------
+        signed : bool, optional
+            If `True`, return signed area from the cross product calculation,
+            indicating whether vertices are ordered clockwise (negative) or
+            counter clockwise (positive).
+        """
+
+        # Use offsets to improve numerical precision
+        x0, y0 = self.mean()
+        x_ = self.vx - x0
+        y_ = self.vy - y0
+
+        # Shoelace formula; in case where the start vertex is not already duplicated
+        # at the end, final term added manually to avoid an array copy.
+        area_main = np.dot(x_[:-1], y_[1:]) - np.dot(y_[:-1], x_[1:])
+        if not (self.vx[-1] == self.vx[0] and self.vy[:-1] == self.vy[0]):
+            area_main += x_[-1] * y_[0] - y_[-1] * x_[0]
+        if signed:
+            return 0.5 * area_main
+        else:
+            return 0.5 * np.abs(area_main)
+
+    def centroid(self):
+        """Return centroid (centre of mass) of polygon."""
+
+        # See http://paulbourke.net/geometry/polygonmesh/
+        #     https://www.ma.ic.ac.uk/~rn/centroid.pdf
+
+        # Use vertex position offsets from mean to improve numerical precision;
+        # for a triangle the mean already identifies the centroid.
+
+        if len(self.vx) == 3:
+            return self.mean()
+        else:
+            x0, y0 = self.mean()
+
+        if self.vx[-1] == self.vx[0] and self.vy[:-1] == self.vy[0]:
+            x_ = self.vx[:-1] - x0
+            y_ = self.vy[:-1] - y0
+        else:
+            x_ = self.vx - x0
+            y_ = self.vy - y0
+        indices = np.arange(len(x_)) - 1
+
+        xs = x_[indices] + x_
+        ys = y_[indices] + y_
+        dxy = x_[indices] * y_ - y_[indices] * x_
+        scl = 1. / (6 * self.area(signed=True))
+
+        return np.dot(xs, dxy) * scl + x0, np.dot(ys, dxy) * scl + y0
+
     def move_to(self, xdelta, ydelta):
         self.vx = list(map(lambda x: x + xdelta, self.vx))
         self.vy = list(map(lambda y: y + ydelta, self.vy))
 
-    def rotate_to(self, theta):
+    def rotate_to(self, theta, center=None):
+        """
+        Rotate polygon by angle `theta` [radian] around `center` (default centroid).
+        """
+
         theta = 0 if theta is None else theta
+        center = self.centroid() if center is None else center
         dtheta = theta - self.theta
         if self.defined() and not np.isclose(dtheta % np.pi, 0.0, atol=1e-9):
-            dx, dy = np.array([self.vx, self.vy]) - np.array(self.center()).reshape(2, 1)
+            dx, dy = np.array([self.vx, self.vy]) - np.array(center).reshape(2, 1)
             self.vx, self.vy = (rotation_matrix_2d(dtheta) @ (dx, dy) +
-                                np.array(self.center()).reshape(2, 1))
+                                np.array(center).reshape(2, 1)).tolist()
         self.theta = theta
 
 
diff --git a/glue/core/tests/test_roi.py b/glue/core/tests/test_roi.py
index 288795fef..789f15e0a 100644
--- a/glue/core/tests/test_roi.py
+++ b/glue/core/tests/test_roi.py
@@ -536,10 +536,43 @@ def test_str(self):
         assert type(str(self.roi)) == str
 
     def test_rotate(self):
+        """ Test 45 deg rotation of square ROI """
         self.define_as_square()
-        self.roi.rotate_to(np.pi / 4)
+        self.roi.rotate_to(np.pi/4)
         assert self.roi.contains([.5, .5, 1.2], [1.2, -0.2, .5]).all()
         assert not self.roi.contains([1.5, 1.5], [0, 0]).any()
+        assert_almost_equal(self.roi.centroid(), (0.5, 0.5), decimal=12)
+
+    def test_rotate_triangle(self):
+        """ Test incremental rotations of triangular (half-square) ROI """
+        self.define_as_square()
+        assert_almost_equal(self.roi.area(), 1, decimal=12)
+        assert_almost_equal(self.roi.centroid(), (0.5, 0.5), decimal=12)
+        self.roi.remove_point(1, 1)
+        assert_almost_equal(self.roi.area(), 0.5, decimal=12)
+        assert_almost_equal(self.roi.centroid(), (1/3, 1/3), decimal=12)
+        self.roi.rotate_to(np.pi/3)
+        assert_almost_equal(self.roi.centroid(), (1/3, 1/3), decimal=12)
+        self.roi.rotate_by(np.pi/6)
+        assert_almost_equal(self.roi.area(), 0.5, decimal=12)
+        assert_almost_equal(self.roi.centroid(), (1/3, 1/3), decimal=12)
+        assert_almost_equal(self.roi.vx, (2/3, -1/3, 2/3), decimal=12)
+        assert_almost_equal(self.roi.vy, (0, 0, 1), decimal=12)
+
+    def test_append_mock_points(self):
+        """
+        Test that adding points on the side of square ROI conserves area and centroid.
+        """
+        self.define_as_square()
+        assert_almost_equal(self.roi.area(), 1, decimal=12)
+        assert_almost_equal(self.roi.area(signed=True), -1, decimal=12)
+        assert_almost_equal(self.roi.centroid(), (0.5, 0.5), decimal=12)
+        assert_almost_equal(self.roi.mean(), (0.5, 0.5), decimal=12)
+        self.roi.add_point(0.9, 0)
+        self.roi.add_point(0.7, 0)
+        assert_almost_equal(self.roi.area(), 1, decimal=12)
+        assert_almost_equal(self.roi.centroid(), (0.5, 0.5), decimal=12)
+        assert_almost_equal(self.roi.mean(), (0.6, 1/3), decimal=12)
 
     def test_serialization(self):
         self.define_as_square()
diff --git a/setup.cfg b/setup.cfg
index 190001cf5..b2cdda28e 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -106,11 +106,11 @@ glue.viewers.profile.qt.tests = data/*.glu
 glue.viewers.scatter.qt.tests = data/*.glu
 
 [flake8]
-ignore = E501,E731,F841,E127,E741,E402,W504,W605
+ignore = E226,E501,E731,F841,E127,E741,E402,W504,W605
 
 [tool:pytest]
 addopts=-p no:logging
-flake8-ignore = E501,E731,F841,E127,E741,E402,W504,W605
+flake8-ignore = E226,E501,E731,F841,E127,E741,E402,W504,W605
 filterwarnings =
     ignore::PendingDeprecationWarning:xlrd
     ignore:Session._key_changed is deprecated

From bb7ca9dff36c82dcf522c8e2ae8a8be3b984a745 Mon Sep 17 00:00:00 2001
From: Derek Homeier <dhomeie@gwdg.de>
Date: Thu, 19 May 2022 16:30:53 +0200
Subject: [PATCH 11/12] Fix `dtype` error on calling `..._inside_poly` with
 lists

---
 glue/core/roi.py       |  2 +-
 glue/utils/geometry.py | 15 +++++++++++++++
 2 files changed, 16 insertions(+), 1 deletion(-)

diff --git a/glue/core/roi.py b/glue/core/roi.py
index 7a01a32f1..27ece95a7 100644
--- a/glue/core/roi.py
+++ b/glue/core/roi.py
@@ -778,7 +778,7 @@ def contains(self, x, y):
         if not isinstance(y, np.ndarray):
             y = np.asarray(y)
 
-        result = points_inside_poly(x, y, self.vx, self.vy)
+        result = points_inside_poly(x, y, np.asarray(self.vx), np.asarray(self.vy))
         return result
 
     # There are several possible definitions of the centre; `mean()` is
diff --git a/glue/utils/geometry.py b/glue/utils/geometry.py
index 26330f8cd..f53c1ee79 100644
--- a/glue/utils/geometry.py
+++ b/glue/utils/geometry.py
@@ -22,6 +22,21 @@ def rotation_matrix_2d(alpha):
 
 
 def points_inside_poly(x, y, vx, vy):
+    """
+    Test if coordinates ``x``, ``y`` fall inside polygon of vertices ``vx``, ``vy``.
+
+    Parameters
+    ----------
+    x, y : `~numpy.ndarray`
+        Coordinates of the points to test
+    vx, vy : `~numpy.ndarray`
+        The vertices of the polygon
+
+    Returns
+    -------
+    contains : `~numpy.ndarray` of bool
+        Array indicating whether each coordinate pair is inside the polygon.
+    """
 
     if x.dtype.kind == 'M' and vx.dtype.kind == 'M':
         vx = vx.astype(x.dtype).astype(float)

From 1298902d8d96cc9edd6915402354be8c43f593dd Mon Sep 17 00:00:00 2001
From: Derek Homeier <dhomeie@gwdg.de>
Date: Thu, 19 May 2022 23:41:57 +0200
Subject: [PATCH 12/12] Adjust default centre for linear PolygonalROI; docs
 updates

---
 glue/core/roi.py            | 40 ++++++++++++++++++++++++++++++++-----
 glue/core/tests/test_roi.py | 14 +++++++++++++
 2 files changed, 49 insertions(+), 5 deletions(-)

diff --git a/glue/core/roi.py b/glue/core/roi.py
index 27ece95a7..f5f0114bc 100644
--- a/glue/core/roi.py
+++ b/glue/core/roi.py
@@ -136,11 +136,25 @@ def to_polygon(self):
         raise NotImplementedError()
 
     def rotate_to(self, theta):
-        """Rotate anticlockwise around center to position angle theta (radian)"""
+        """
+        Rotate anticlockwise around center to position angle theta.
+
+        Parameters
+        ----------
+        theta : float
+            Angle of anticlockwise rotation around center in radian.
+        """
         raise NotImplementedError()
 
     def rotate_by(self, dtheta, **kwargs):
-        """Rotate the Roi around center by angle dtheta (radian)"""
+        """
+        Rotate the Roi around center by angle dtheta.
+
+        Parameters
+        ----------
+        dtheta : float
+            Change in anticlockwise rotation angle around center in radian.
+        """
         self.rotate_to(getattr(self, 'theta', 0.0) + dtheta, **kwargs)
 
     def copy(self):
@@ -697,9 +711,10 @@ def add_point(self, x, y):
         self.vy.append(y)
 
     def reset(self):
-        """Reset the vertex lists"""
+        """Reset the vertex lists and position angle"""
         self.vx = []
         self.vy = []
+        self.theta = 0
 
     def replace_last_point(self, x, y):
         if len(self.vx) > 0:
@@ -856,12 +871,27 @@ def move_to(self, xdelta, ydelta):
 
     def rotate_to(self, theta, center=None):
         """
-        Rotate polygon by angle `theta` [radian] around `center` (default centroid).
+        Rotate polygon to position angle `theta` around `center`.
+
+        Parameters
+        ----------
+        theta : float
+            Angle of anticlockwise rotation around center in radian.
+        center : pair of float, optional
+            Coordinates of center of rotation. Defaults to
+            :meth:`~glue.core.roi.PolygonalROI.centroid`, for linear
+            "polygons" to :meth:`~glue.core.roi.PolygonalROI.mean`.
         """
 
         theta = 0 if theta is None else theta
-        center = self.centroid() if center is None else center
+        # For linear (1D) "polygons" centroid is not defined.
+        if center is None:
+            if self.area() == 0:
+                center = self.mean()
+            else:
+                center = self.centroid()
         dtheta = theta - self.theta
+
         if self.defined() and not np.isclose(dtheta % np.pi, 0.0, atol=1e-9):
             dx, dy = np.array([self.vx, self.vy]) - np.array(center).reshape(2, 1)
             self.vx, self.vy = (rotation_matrix_2d(dtheta) @ (dx, dy) +
diff --git a/glue/core/tests/test_roi.py b/glue/core/tests/test_roi.py
index 789f15e0a..e5aeac2f7 100644
--- a/glue/core/tests/test_roi.py
+++ b/glue/core/tests/test_roi.py
@@ -559,6 +559,20 @@ def test_rotate_triangle(self):
         assert_almost_equal(self.roi.vx, (2/3, -1/3, 2/3), decimal=12)
         assert_almost_equal(self.roi.vy, (0, 0, 1), decimal=12)
 
+    def test_rotate_polyline(self):
+        """ Test rotation of degenerate (linear) ROI around mean """
+        self.roi.reset()
+        self.roi.add_point(-2, 0)
+        self.roi.add_point(4, 0)
+        assert_almost_equal(self.roi.mean(), (1.0, 0.0), decimal=12)
+        self.roi.add_point(-0.5, 0)
+        self.roi.add_point(-1.5, 0)
+        assert_almost_equal(self.roi.mean(), (0.0, 0.0), decimal=12)
+        assert all(np.isnan(self.roi.centroid()))
+        self.roi.rotate_to(np.pi/2)
+        assert_almost_equal(self.roi.vx, (0, 0, 0, 0), decimal=12)
+        assert_almost_equal(self.roi.vy, (-2, 4, -0.5, -1.5), decimal=12)
+
     def test_append_mock_points(self):
         """
         Test that adding points on the side of square ROI conserves area and centroid.