Merge pull request #87 from astrofrog/polygon

Added initial implementation of Polygon to_mask and as_patch

regions/_geometry/__init__.py

@@ -6,7 +6,8 @@
 from .circular_overlap import *
 from .elliptical_overlap import *
 from .rectangular_overlap import *
+from .polygonal_overlap import *
 
 
 __all__ = ['circular_overlap_grid', 'elliptical_overlap_grid',
-           'rectangular_overlap_grid']
+           'rectangular_overlap_grid', 'polygonal_overlap_grid']

regions/_geometry/core.pxd

@@ -1,3 +1,7 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+
+# This file is needed in order to be able to cimport functions into other Cython files
+
 cdef double distance(double x1, double y1, double x2, double y2)
 cdef double area_arc(double x1, double y1, double x2, double y2, double R)
 cdef double area_triangle(double x1, double y1, double x2, double y2, double x3, double y3)

regions/_geometry/pnpoly.pxd

@@ -0,0 +1,9 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+
+# This file is needed in order to be able to cimport functions into other Cython files
+
+cimport numpy as np
+
+ctypedef np.float64_t DTYPE_t
+
+cdef int point_in_polygon(double x, double y, np.ndarray[DTYPE_t, ndim=1] vx, np.ndarray[DTYPE_t, ndim=1] vy)

regions/_geometry/pnpoly.pyx

@@ -0,0 +1,116 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+#
+# The code in this file was adapted from code written by Greg von Winckel:
+#
+# https://github.com/gregvw/pnpoly
+#
+# and which was released under the following license:
+#
+# ----------------------------------------------------------------------------
+#
+# The MIT License (MIT)
+#
+# Copyright (c) 2014 Greg von Winckel
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy of
+# this software and associated documentation files (the "Software"), to deal in
+# the Software without restriction, including without limitation the rights to
+# use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of
+# the Software, and to permit persons to whom the Software is furnished to do so,
+# subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in all
+# copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
+# FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
+# COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
+# IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
+# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+#
+# ----------------------------------------------------------------------------
+#
+# This code was itself adapted from code written by W. Randolph Franklin:
+#
+# http://www.ecse.rpi.edu/~wrf/Research/Short_Notes/pnpoly.html
+#
+# and released under the following license:
+#
+# Copyright (c) 1970-2003, Wm. Randolph Franklin
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# Redistributions of source code must retain the above copyright notice, this
+# list of conditions and the following disclaimers.
+#
+# Redistributions in binary form must reproduce the above copyright notice in
+# the documentation and/or other materials provided with the distribution.
+#
+# The name of W. Randolph Franklin may not be used to endorse or promote
+# products derived from this Software without specific prior written permission.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+# SOFTWARE.
+
+import numpy as np
+cimport numpy as np
+
+DTYPE_BOOL = np.bool
+ctypedef np.uint8_t DTYPE_BOOL_t
+
+
+def points_in_polygon(np.ndarray[DTYPE_t, ndim=1] x,
+                      np.ndarray[DTYPE_t, ndim=1] y,
+                      np.ndarray[DTYPE_t, ndim=1] vx,
+                      np.ndarray[DTYPE_t, ndim=1] vy):
+
+    cdef int i, n
+    cdef np.ndarray[np.uint8_t, ndim=1, cast=True] result
+
+    n = x.shape[0]
+
+    result = np.zeros(n, DTYPE_BOOL)
+
+    for i in range(n):
+        result[i] = point_in_polygon(x[i], y[i], vx, vy)
+
+    return result
+
+
+cdef int point_in_polygon(double x, double y,
+                          np.ndarray[DTYPE_t, ndim=1] vx,
+                          np.ndarray[DTYPE_t, ndim=1] vy):
+    """
+    Determine whether a test point (x, y) is within a polygon defined by a set
+    of vertices (vx, vy).
+
+    This uses the even-odd rule, as described here:
+
+        https://en.wikipedia.org/wiki/Even–odd_rule
+    """
+
+    cdef int i, j, k, m, n
+    cdef int result
+
+    n = vx.shape[0]
+
+    result = 0
+
+    for i in range(n):
+        j = (i + n - 1) % n
+        if(((vy[i] > y) != (vy[j] > y)) and
+            (x < (vx[j] - vx[i]) * (y - vy[i]) / (vy[j] - vy[i]) + vx[i])):
+            result += 1
+
+    return result % 2

regions/_geometry/polygonal_overlap.pyx

@@ -0,0 +1,117 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+
+from __future__ import (absolute_import, division, print_function,
+                        unicode_literals)
+
+import numpy as np
+cimport numpy as np
+
+
+__all__ = ['polygonal_overlap_grid']
+
+
+DTYPE = np.float64
+ctypedef np.float64_t DTYPE_t
+
+from .pnpoly cimport point_in_polygon
+
+__all__ = ['polygonal_overlap_grid']
+
+
+def polygonal_overlap_grid(double xmin, double xmax, double ymin, double ymax,
+                         int nx, int ny,
+                         np.ndarray[DTYPE_t, ndim=1] vx,
+                         np.ndarray[DTYPE_t, ndim=1] vy,
+                         int use_exact,
+                         int subpixels):
+    """
+    polygonal_overlap_grid(xmin, xmax, ymin, ymax, nx, ny, r,
+                         use_exact, subpixels)
+
+    Area of overlap between a polygon and a pixel grid.
+
+    Parameters
+    ----------
+    xmin, xmax, ymin, ymax : float
+        Extent of the grid in the x and y direction.
+    nx, ny : int
+        Grid dimensions.
+    vx, vy : `numpy.ndarray`
+        The vertices of the polygon
+    use_exact : 0 or 1
+        If ``1`` calculates exact overlap, if ``0`` uses ``subpixel`` number
+        of subpixels to calculate the overlap.
+    subpixels : int
+        Each pixel resampled by this factor in each dimension, thus each
+        pixel is divided into ``subpixels ** 2`` subpixels.
+
+    Returns
+    -------
+    frac : `~numpy.ndarray` (float)
+        2-d array of shape (ny, nx) giving the fraction of the overlap.
+    """
+
+    cdef unsigned int i, j
+    cdef double x, y, dx, dy, d, pixel_radius
+    cdef double bxmin, bxmax, bymin, bymax
+    cdef double pxmin, pxcen, pxmax, pymin, pycen, pymax
+
+    # Define output array
+    cdef np.ndarray[DTYPE_t, ndim=2] frac = np.zeros([ny, nx], dtype=DTYPE)
+
+    if use_exact == 1:
+        raise NotImplementedError("Exact mode has not been implemented for "
+                                  "polygonal apertures")
+
+    # Find the width of each element in x and y
+    dx = (xmax - xmin) / nx
+    dy = (ymax - ymin) / ny
+
+    # Define bounding box
+    bxmin = vx.min()
+    bxmax = vx.max()
+    bymin = vy.min()
+    bymax = vy.max()
+
+    for i in range(nx):
+        pxmin = xmin + i * dx  # lower end of pixel
+        pxmax = pxmin + dx  # upper end of pixel
+        if pxmax > bxmin and pxmin < bxmax:
+            for j in range(ny):
+                pymin = ymin + j * dy
+                pymax = pymin + dy
+                if pymax > bymin and pymin < bymax:
+                    frac[j, i] = polygonal_overlap_single_subpixel(pxmin, pymin,
+                                                                 pxmax, pymax,
+                                                                 vx, vy, subpixels)
+
+    return frac
+
+
+cdef double polygonal_overlap_single_subpixel(double x0, double y0,
+                                            double x1, double y1,
+                                            np.ndarray[DTYPE_t, ndim=1] vx,
+                                            np.ndarray[DTYPE_t, ndim=1] vy,
+                                            int subpixels):
+    """
+    Return the fraction of overlap between a polygon and a single pixel
+    with given extent, using a sub-pixel sampling method.
+    """
+
+    cdef unsigned int i, j
+    cdef double x, y, dx, dy
+    cdef double frac = 0.  # Accumulator.
+
+    dx = (x1 - x0) / subpixels
+    dy = (y1 - y0) / subpixels
+
+    x = x0 - 0.5 * dx
+    for i in range(subpixels):
+        x += dx
+        y = y0 - 0.5 * dy
+        for j in range(subpixels):
+            y += dy
+            if point_in_polygon(x, y, vx, vy) == 1:
+                frac += 1.
+
+    return frac / (subpixels * subpixels)

regions/shapes/polygon.py

@@ -1,6 +1,12 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
+
 from __future__ import absolute_import, division, print_function, unicode_literals
-from ..core import PixelRegion, SkyRegion
+
+import numpy as np
+
+from ..core import PixelRegion, SkyRegion, Mask, BoundingBox
+from .._geometry import polygonal_overlap_grid
+from .._geometry.pnpoly import points_in_polygon
 
 __all__ = ['PolygonPixelRegion', 'PolygonSkyRegion']
 
@@ -35,8 +41,15 @@ def area(self):
         raise NotImplementedError
 
     def contains(self, pixcoord):
-        # TODO: needs to be implemented
-        raise NotImplementedError
+        if pixcoord.isscalar:
+            x = np.array([pixcoord.x], dtype=float)
+            y = np.array([pixcoord.y], dtype=float)
+        else:
+            x = pixcoord.x.astype(float)
+            y = pixcoord.y.astype(float)
+        return points_in_polygon(x, y,
+                                 self.vertices.x.astype(float),
+                                 self.vertices.y.astype(float)).astype(bool)
 
     def to_shapely(self):
         # TODO: needs to be implemented
@@ -48,16 +61,50 @@ def to_sky(self, wcs, mode='local', tolerance=None):
 
     @property
     def bounding_box(self, mode='center'):
-        # TODO: needs to be implemented
-        raise NotImplementedError
+        xmin = self.vertices.x.min()
+        xmax = self.vertices.x.max()
+        ymin = self.vertices.y.min()
+        ymax = self.vertices.y.max()
+        return BoundingBox._from_float(xmin, xmax, ymin, ymax)
 
-    def to_mask(self, mode='center'):
-        # TODO: needs to be implemented
-        raise NotImplementedError
+    def to_mask(self, mode='center', subpixels=5):
+
+        self._validate_mode(mode, subpixels)
+
+        if mode == 'center':
+            mode = 'subpixels'
+            subpixels = 1
+
+        # Find bounding box and mask size
+        bbox = self.bounding_box
+        ny, nx = bbox.shape
+
+        # Find position of pixel edges and recenter so that circle is at origin
+        xmin = float(bbox.ixmin) - 0.5
+        xmax = float(bbox.ixmax) - 0.5
+        ymin = float(bbox.iymin) - 0.5
+        ymax = float(bbox.iymax) - 0.5
+
+        if mode == 'subpixels':
+            use_exact = 0
+        else:
+            use_exact = 1
+
+        fraction = polygonal_overlap_grid(
+            xmin, xmax, ymin, ymax, nx, ny,
+            self.vertices.x.astype(float),
+            self.vertices.y.astype(float),
+            use_exact, subpixels)
+
+        return Mask(fraction, bbox=bbox)
 
     def as_patch(self, **kwargs):
-        # TODO: needs to be implemented
-        raise NotImplementedError
+        """
+        Matplotlib patch object for this region (`matplotlib.patches.Polygon`).
+        """
+        from matplotlib.patches import Polygon
+        xy = np.vstack([self.vertices.x, self.vertices.y]).transpose()
+        return Polygon(xy=xy, **kwargs)
 
 
 class PolygonSkyRegion(SkyRegion):

regions/shapes/tests/reference/test_to_mask_poly-mode_center.txt

@@ -0,0 +1,7 @@
+1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
+1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00
+0.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00
+0.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00
+0.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00
+0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
+0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00

regions/shapes/tests/reference/test_to_mask_poly-mode_subpixels-subpixels_1.txt

@@ -0,0 +1,7 @@
+1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
+1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00
+0.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00
+0.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00
+0.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00
+0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
+0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00

regions/shapes/tests/reference/test_to_mask_poly-mode_subpixels-subpixels_18.txt

@@ -0,0 +1,7 @@
+4.19753086e-01 3.98148148e-01 2.31481481e-01 6.17283951e-02 0.00000000e+00 0.00000000e+00 0.00000000e+00
+6.04938272e-01 1.00000000e+00 1.00000000e+00 1.00000000e+00 8.88888889e-01 6.94444444e-01 1.85185185e-02
+4.07407407e-01 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 5.83333333e-01 0.00000000e+00
+1.94444444e-01 1.00000000e+00 1.00000000e+00 1.00000000e+00 9.72222222e-01 1.11111111e-01 0.00000000e+00
+2.16049383e-02 9.72222222e-01 1.00000000e+00 1.00000000e+00 6.01851852e-01 0.00000000e+00 0.00000000e+00
+0.00000000e+00 3.91975309e-01 8.88888889e-01 9.81481481e-01 1.29629630e-01 0.00000000e+00 0.00000000e+00
+0.00000000e+00 0.00000000e+00 2.77777778e-02 2.56172840e-01 0.00000000e+00 0.00000000e+00 0.00000000e+00

regions/shapes/tests/reference/test_to_mask_poly-mode_subpixels-subpixels_5.txt

@@ -0,0 +1,7 @@
+4.40000000e-01 4.00000000e-01 2.00000000e-01 4.00000000e-02 0.00000000e+00 0.00000000e+00 0.00000000e+00
+6.00000000e-01 1.00000000e+00 1.00000000e+00 1.00000000e+00 8.80000000e-01 6.80000000e-01 0.00000000e+00
+4.00000000e-01 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 5.60000000e-01 0.00000000e+00
+2.00000000e-01 1.00000000e+00 1.00000000e+00 1.00000000e+00 9.60000000e-01 8.00000000e-02 0.00000000e+00
+0.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 6.00000000e-01 0.00000000e+00 0.00000000e+00
+0.00000000e+00 4.00000000e-01 9.20000000e-01 1.00000000e+00 1.60000000e-01 0.00000000e+00 0.00000000e+00
+0.00000000e+00 0.00000000e+00 4.00000000e-02 2.80000000e-01 0.00000000e+00 0.00000000e+00 0.00000000e+00