diff --git a/.rtd-environment.yml b/.rtd-environment.yml
index 729e4528..2b3070e8 100644
--- a/.rtd-environment.yml
+++ b/.rtd-environment.yml
@@ -5,3 +5,4 @@ dependencies:
 - matplotlib
 - astropy
 - wcsaxes
+- cython
diff --git a/.travis.yml b/.travis.yml
index d82151e3..b78fce42 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -25,7 +25,7 @@ env:
         - ASTROPY_VERSION=stable
         - SETUP_CMD='test'
         - MAIN_CMD='python setup.py'
-        - PIP_DEPENDENCIES=''
+        - PIP_DEPENDENCIES='https://github.com/astrofrog/pytest-arraydiff/archive/master.zip'
         - CONDA_DEPENDENCIES='Cython shapely wcsaxes'
         - CONDA_CHANNELS='astropy'
         - SETUP_XVFB=True
@@ -53,7 +53,7 @@ matrix:
         # Do a test without the optional dependencies
         - python: 3.5
           env: SETUP_CMD='test'
-               CONDA_DEPENDENCIES=''
+               CONDA_DEPENDENCIES='Cython'
 
 
         # Try Astropy development version
@@ -68,10 +68,10 @@ matrix:
         # CONDA_CHANNELS environmental variable).
         - python: 3.3
           env: SETUP_CMD='egg_info'
-               CONDA_DEPENDENCIES=''
+               CONDA_DEPENDENCIES='Cython'
         - python: 3.3
           env: SETUP_CMD='test' NUMPY_VERSION=1.9
-               CONDA_DEPENDENCIES=''
+               CONDA_DEPENDENCIES='Cython'
 
         # Try older numpy versions
         - python: 2.7
diff --git a/appveyor.yml b/appveyor.yml
index 8de3c156..256b601c 100644
--- a/appveyor.yml
+++ b/appveyor.yml
@@ -15,6 +15,7 @@ environment:
       # Cython code, you can set CONDA_DEPENDENCIES=''
       CONDA_CHANNELS: "defaults conda-forge"
       CONDA_DEPENDENCIES: "Cython shapely"
+      PIP_DEPENDENCIES: "https://github.com/astrofrog/pytest-arraydiff/archive/master.zip"
 
   matrix:
 
diff --git a/docs/conf.py b/docs/conf.py
index bd9d64e7..9773af3a 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -163,3 +163,15 @@
 
     edit_on_github_source_root = ""
     edit_on_github_doc_root = "docs"
+
+plot_rcparams = {'savefig.bbox': 'tight',
+                 'savefig.facecolor':'none',
+                 'axes.formatter.useoffset': False,
+                 'xtick.labelsize': 9,
+                 'ytick.labelsize': 9,
+                 'axes.labelsize': 11,
+                 'figure.titlesize': 11,
+                 'figure.subplot.wspace': 0.23,
+                 'figure.subplot.hspace': 0.23}
+
+plot_apply_rcparams = True
diff --git a/docs/getting_started.rst b/docs/getting_started.rst
index 453706d7..5caaa6a0 100644
--- a/docs/getting_started.rst
+++ b/docs/getting_started.rst
@@ -402,17 +402,217 @@ or `regions.PixelRegion.contains` methods:
 Masks
 =====
 
-For aperture photometry, a common operation is to compute, for a given image and region,
-a boolean mask or array of pixel indices defining which pixels (in the whole image or a
-minimal rectangular bounding box) are inside and outside the region.
+For aperture photometry, a common operation is to compute, for a given image and
+region, a mask or array of pixel indices defining which pixels (in the whole
+image or a minimal rectangular bounding box) are inside and outside the region.
+
+All :class:`~regions.PixelRegion` objects have a
+:meth:`~regions.PixelRegion.to_mask` method that returns a
+:class:`~regions.Mask` object that contains information about whether
+pixels are inside the region, and can be used to mask data arrays:
+
+    >>> from regions.core import PixCoord
+    >>> from regions.shapes.circle import CirclePixelRegion
+    >>> center = PixCoord(4., 5.)
+    >>> reg = CirclePixelRegion(center, 2.3411)
+    >>> mask = reg.to_mask()
+    >>> mask
+    <regions.mask.Mask object at 0x10900b5c0>
+    >>> mask.data
+    array([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
+           [ 0.,  0.,  1.,  1.,  1.,  0.,  0.],
+           [ 0.,  1.,  1.,  1.,  1.,  1.,  0.],
+           [ 0.,  1.,  1.,  1.,  1.,  1.,  0.],
+           [ 0.,  1.,  1.,  1.,  1.,  1.,  0.],
+           [ 0.,  0.,  1.,  1.,  1.,  0.,  0.],
+           [ 0.,  0.,  0.,  0.,  0.,  0.,  0.]])
+
+The mask data contains floating point that are between 0 (no overlap) and 1
+(overlap). By default, this is determined by looking only at the central position
+in each pixel, and::
+
+    >>> reg.to_mask()
+
+is equivalent to::
+
+    >>> reg.to_mask(mode='center')
+
+but other modes are available:
+
+* ``mode='exact'``: the overlap is determined using the exact geometrical
+  overlap between pixels and the region. This is slower than using the central
+  position, but allows partial overlap to be treated correctly.
+
+* ``mode='subpixels'``: the overlap is determined by sub-sampling the pixel
+  using a grid of sub-pixels. The number of sub-pixels to use in this mode
+  should be given using the ``subpixels`` argument.
+
+Here are what the different modes look like:
+
+.. plot::
+   :include-source:
+
+    import matplotlib.pyplot as plt
+    from regions.core import PixCoord
+    from regions.shapes.circle import CirclePixelRegion
+
+    center = PixCoord(6.6, 7.2)
+    reg = CirclePixelRegion(center, 5.2)
+
+    plt.figure(figsize=(6, 6))
+
+    mask1 = reg.to_mask(mode='center')
+    plt.subplot(2, 2, 1)
+    plt.title("mode='center'", size=9)
+    plt.imshow(mask1.data, cmap=plt.cm.viridis,
+               interpolation='nearest', origin='lower')
+
+    mask2 = reg.to_mask(mode='exact')
+    plt.subplot(2, 2, 2)
+    plt.title("mode='exact'", size=9)
+    plt.imshow(mask2.data, cmap=plt.cm.viridis,
+               interpolation='nearest', origin='lower')
+
+    mask3 = reg.to_mask(mode='subpixels', subpixels=3)
+    plt.subplot(2, 2, 3)
+    plt.title("mode='subpixels', subpixels=3", size=9)
+    plt.imshow(mask3.data, cmap=plt.cm.viridis,
+               interpolation='nearest', origin='lower')
+
+    mask4 = reg.to_mask(mode='subpixels', subpixels=20)
+    plt.subplot(2, 2, 4)
+    plt.title("mode='subpixels', subpixels=20", size=9)
+    plt.imshow(mask4.data, cmap=plt.cm.viridis,
+               interpolation='nearest', origin='lower')
+
+As we've seen above, the :class:`~regions.Mask` objects have a ``data``
+attribute that contains a Numpy array with the mask values. However, if you
+have for example a circular region with a radius of 3 pixels at a pixel position
+of (1000, 1000), it would be inefficient to store a mask that has a size larger
+than this, so instead we store the mask using the minimal array that contains
+the mask, and the :class:`~regions.Mask` objects also include a ``bbox``
+attribute that is a :class:`~regions.BoundingBox` object used to indicate where
+the mask should be applied in an image.
+
+:class:`~regions.Mask` objects also have a number of methods to make it
+easy to use the masks with data. The :meth:`~regions.Mask.to_image` method
+can be used to obtain an image of the mask in a 2D array of the given shape.
+This places the mask in the correct place in the image and deals properly with
+boundary effects:
+
+.. plot::
+   :include-source:
+
+    import matplotlib.pyplot as plt
+    from regions.core import PixCoord
+    from regions.shapes.circle import CirclePixelRegion
+
+    center = PixCoord(6.6, 7.2)
+    reg = CirclePixelRegion(center, 5.2)
+
+    mask = reg.to_mask(mode='exact')
+    plt.figure(figsize=(4, 4))
+    plt.imshow(mask.to_image((10, 10)), cmap=plt.cm.viridis,
+               interpolation='nearest', origin='lower')
+
+The :meth:`~regions.Mask.cutout` method can be used to create a cutout from
+the input data over the mask bounding box, and the
+:meth:`~regions.Mask.multiply` method can be used to multiply the aperture
+mask with the input data to create a mask-weighted data cutout. All of these
+methods properly handle the cases of partial or no overlap of the aperture mask
+with the data.
 
-To a certain degree, such spatial filtering can be done using the methods described in the previous :ref:`gs-contain`
-section. Apart from that, no high-level functionality for spatial filtering, bounding boxes or aperture photometry
-is available yet.
+These masks can be used as the building blocks for photometry, which we
+demonstrate with a simple example. We start off by getting an example image:
+
+.. plot::
+   :context: reset
+   :include-source:
+   :align: center
+   :nofigs:
+
+    >>> from astropy.io import fits
+    >>> from astropy.utils.data import get_pkg_data_filename
+    >>> filename = get_pkg_data_filename('photometry/M6707HH.fits')
+    >>> hdu = fits.open(filename)[0]
+
+We then define the aperture:
+
+.. plot::
+   :context:
+   :include-source:
+   :align: center
+   :nofigs:
 
-For now, please use `photutils`_ or `pyregion`_.
+    >>> from regions.core import PixCoord
+    >>> from regions.shapes.circle import CirclePixelRegion
+    >>> center = PixCoord(158.5, 1053.5)
+    >>> aperture = CirclePixelRegion(center, 4.)
+
+We convert the aperture to a mask and extract a cutout from the data, as well
+as a cutout with the data multipled by the mask:
+
+.. plot::
+   :context:
+   :include-source:
+   :align: center
+   :nofigs:
+
+    >>> mask = aperture.to_mask(mode='exact')
+    >>> data = mask.cutout(hdu.data)
+    >>> weighted_data = mask.multiply(hdu.data)
+
+We can take a look at the results to make sure the source overlaps with the
+aperture:
+
+.. plot::
+   :context:
+   :include-source:
+   :align: center
+
+    >>> import matplotlib.pyplot as plt
+    >>> plt.subplot(1,3,1)
+    >>> plt.title("Mask", size=9)
+    >>> plt.imshow(mask.data, cmap=plt.cm.viridis,
+    ...            interpolation='nearest', origin='lower',
+    ...            extent=mask.bbox.extent)
+    >>> plt.subplot(1,3,2)
+    >>> plt.title("Data cutout", size=9)
+    >>> plt.imshow(data, cmap=plt.cm.viridis,
+    ...            interpolation='nearest', origin='lower',
+    ...            extent=mask.bbox.extent)
+    >>> plt.subplot(1,3,3)
+    >>> plt.title("Data cutout multiplied by mask", size=9)
+    >>> plt.imshow(weighted_data, cmap=plt.cm.viridis,
+    ...            interpolation='nearest', origin='lower',
+    ...            extent=mask.bbox.extent)
+
+We can also use the ``Mask.bbox`` attribute to look at the extent
+of the mask in the image:
+
+.. plot::
+   :context:
+   :include-source:
+   :align: center
+
+    >>> ax = plt.subplot(1, 1, 1)
+    >>> ax.imshow(hdu.data, cmap=plt.cm.viridis,
+    ...            interpolation='nearest', origin='lower')
+    >>> ax.add_patch(mask.bbox.as_patch(facecolor='none', edgecolor='white'))
+    >>> ax.add_patch(aperture.as_patch(facecolor='none', edgecolor='orange'))
+    >>> ax.set_xlim(120, 180)
+    >>> ax.set_ylim(1020, 1080)
+
+Finally, we can use the mask and data values to compute weighted statistics:
+
+.. plot::
+   :context:
+   :include-source:
+   :align: center
+   :nofigs:
 
-(We plan to merge that functionality from ``photutils`` or ``pyregion`` into this ``regions`` package, or re-implement it.)
+    >>> np.average(data, weights=mask)
+    9364.0126748880211
 
 .. _gs-compound:
 
diff --git a/regions/_geometry/README.md b/regions/_geometry/README.md
new file mode 100644
index 00000000..e1f43834
--- /dev/null
+++ b/regions/_geometry/README.md
@@ -0,0 +1 @@
+# This sub-package contains fast geometrical routines that are then used by different shapes
diff --git a/regions/_geometry/__init__.py b/regions/_geometry/__init__.py
new file mode 100644
index 00000000..d43c2999
--- /dev/null
+++ b/regions/_geometry/__init__.py
@@ -0,0 +1,12 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+"""
+Geometry subpackage for low-level geometry functions.
+"""
+
+from .circular_overlap import *
+from .elliptical_overlap import *
+from .rectangular_overlap import *
+
+
+__all__ = ['circular_overlap_grid', 'elliptical_overlap_grid',
+           'rectangular_overlap_grid']
diff --git a/regions/_geometry/circular_overlap.pyx b/regions/_geometry/circular_overlap.pyx
new file mode 100644
index 00000000..c65e88d4
--- /dev/null
+++ b/regions/_geometry/circular_overlap.pyx
@@ -0,0 +1,231 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+"""
+The functions defined here allow one to determine the exact area of
+overlap of a rectangle and a circle (written by Thomas Robitaille).
+"""
+
+from __future__ import (absolute_import, division, print_function,
+                        unicode_literals)
+
+import numpy as np
+cimport numpy as np
+
+
+__all__ = ['circular_overlap_grid']
+
+
+cdef extern from "math.h":
+
+    double asin(double x)
+    double sin(double x)
+    double sqrt(double x)
+
+
+DTYPE = np.float64
+ctypedef np.float64_t DTYPE_t
+
+# NOTE: Here we need to make sure we use cimport to import the C functions from
+# core (since these were defined with cdef). This also requires the core.pxd
+# file to exist with the function signatures.
+from .core cimport area_arc, area_triangle
+
+
+def circular_overlap_grid(double xmin, double xmax, double ymin, double ymax,
+                          int nx, int ny, double r, int use_exact,
+                          int subpixels):
+    """
+    circular_overlap_grid(xmin, xmax, ymin, ymax, nx, ny, r,
+                             use_exact, subpixels)
+
+    Area of overlap between a circle and a pixel grid. The circle is centered
+    on the origin.
+
+    Parameters
+    ----------
+    xmin, xmax, ymin, ymax : float
+        Extent of the grid in the x and y direction.
+    nx, ny : int
+        Grid dimensions.
+    r : float
+        The radius of the circle.
+    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)
+
+    # Find the width of each element in x and y
+    dx = (xmax - xmin) / nx
+    dy = (ymax - ymin) / ny
+
+    # Find the radius of a single pixel
+    pixel_radius = 0.5 * sqrt(dx * dx + dy * dy)
+
+    # Define bounding box
+    bxmin = -r - 0.5 * dx
+    bxmax = +r + 0.5 * dx
+    bymin = -r - 0.5 * dy
+    bymax = +r + 0.5 * dy
+
+    for i in range(nx):
+        pxmin = xmin + i * dx  # lower end of pixel
+        pxcen = pxmin + dx * 0.5
+        pxmax = pxmin + dx  # upper end of pixel
+        if pxmax > bxmin and pxmin < bxmax:
+            for j in range(ny):
+                pymin = ymin + j * dy
+                pycen = pymin + dy * 0.5
+                pymax = pymin + dy
+                if pymax > bymin and pymin < bymax:
+
+                    # Distance from circle center to pixel center.
+                    d = sqrt(pxcen * pxcen + pycen * pycen)
+
+                    # If pixel center is "well within" circle, count full
+                    # pixel.
+                    if d < r - pixel_radius:
+                        frac[j, i] = 1.
+
+                    # If pixel center is "close" to circle border, find
+                    # overlap.
+                    elif d < r + pixel_radius:
+
+                        # Either do exact calculation or use subpixel
+                        # sampling:
+                        if use_exact:
+                            frac[j, i] = circular_overlap_single_exact(
+                                pxmin, pymin, pxmax, pymax, r) / (dx * dy)
+                        else:
+                            frac[j, i] = circular_overlap_single_subpixel(
+                                pxmin, pymin, pxmax, pymax, r, subpixels)
+
+                    # Otherwise, it is fully outside circle.
+                    # No action needed.
+
+    return frac
+
+
+# NOTE: The following two functions use cdef because they are not
+# intended to be called from the Python code. Using def makes them
+# callable from outside, but also slower. In any case, these aren't useful
+# to call from outside because they only operate on a single pixel.
+
+
+cdef double circular_overlap_single_subpixel(double x0, double y0,
+                                             double x1, double y1,
+                                             double r, int subpixels):
+    """Return the fraction of overlap between a circle and a single pixel
+    with given extent, using a sub-pixel sampling method."""
+
+    cdef unsigned int i, j
+    cdef double x, y, dx, dy, r_squared
+    cdef double frac = 0.  # Accumulator.
+
+    dx = (x1 - x0) / subpixels
+    dy = (y1 - y0) / subpixels
+    r_squared = r ** 2
+
+    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 x * x + y * y < r_squared:
+                frac += 1.
+
+    return frac / (subpixels * subpixels)
+
+
+cdef double circular_overlap_single_exact(double xmin, double ymin,
+                                          double xmax, double ymax,
+                                          double r):
+    """
+    Area of overlap of a rectangle and a circle
+    """
+    if 0. <= xmin:
+        if 0. <= ymin:
+            return circular_overlap_core(xmin, ymin, xmax, ymax, r)
+        elif 0. >= ymax:
+            return circular_overlap_core(-ymax, xmin, -ymin, xmax, r)
+        else:
+            return circular_overlap_single_exact(xmin, ymin, xmax, 0., r) \
+                + circular_overlap_single_exact(xmin, 0., xmax, ymax, r)
+    elif 0. >= xmax:
+        if 0. <= ymin:
+            return circular_overlap_core(-xmax, ymin, -xmin, ymax, r)
+        elif 0. >= ymax:
+            return circular_overlap_core(-xmax, -ymax, -xmin, -ymin, r)
+        else:
+            return circular_overlap_single_exact(xmin, ymin, xmax, 0., r) \
+                + circular_overlap_single_exact(xmin, 0., xmax, ymax, r)
+    else:
+        if 0. <= ymin:
+            return circular_overlap_single_exact(xmin, ymin, 0., ymax, r) \
+                + circular_overlap_single_exact(0., ymin, xmax, ymax, r)
+        if 0. >= ymax:
+            return circular_overlap_single_exact(xmin, ymin, 0., ymax, r) \
+                + circular_overlap_single_exact(0., ymin, xmax, ymax, r)
+        else:
+            return circular_overlap_single_exact(xmin, ymin, 0., 0., r) \
+                + circular_overlap_single_exact(0., ymin, xmax, 0., r) \
+                + circular_overlap_single_exact(xmin, 0., 0., ymax, r) \
+                + circular_overlap_single_exact(0., 0., xmax, ymax, r)
+
+
+def circular_overlap_core(double xmin, double ymin, double xmax, double ymax,
+                          double r):
+    """Assumes that the center of the circle is <= xmin,
+    ymin (can always modify input to conform to this).
+    """
+
+    cdef double area, d1, d2, x1, x2, y1, y2
+
+    if xmin * xmin + ymin * ymin > r * r:
+        area = 0.
+    elif xmax * xmax + ymax * ymax < r * r:
+        area = (xmax - xmin) * (ymax - ymin)
+    else:
+        area = 0.
+        d1 = sqrt(xmax * xmax + ymin * ymin)
+        d2 = sqrt(xmin * xmin + ymax * ymax)
+        if d1 < r and d2 < r:
+            x1, y1 = sqrt(r * r - ymax * ymax), ymax
+            x2, y2 = xmax, sqrt(r * r - xmax * xmax)
+            area = ((xmax - xmin) * (ymax - ymin) -
+                    area_triangle(x1, y1, x2, y2, xmax, ymax) +
+                    area_arc(x1, y1, x2, y2, r))
+        elif d1 < r:
+            x1, y1 = xmin, sqrt(r * r - xmin * xmin)
+            x2, y2 = xmax, sqrt(r * r - xmax * xmax)
+            area = (area_arc(x1, y1, x2, y2, r) +
+                    area_triangle(x1, y1, x1, ymin, xmax, ymin) +
+                    area_triangle(x1, y1, x2, ymin, x2, y2))
+        elif d2 < r:
+            x1, y1 = sqrt(r * r - ymin * ymin), ymin
+            x2, y2 = sqrt(r * r - ymax * ymax), ymax
+            area = (area_arc(x1, y1, x2, y2, r) +
+                    area_triangle(x1, y1, xmin, y1, xmin, ymax) +
+                    area_triangle(x1, y1, xmin, y2, x2, y2))
+        else:
+            x1, y1 = sqrt(r * r - ymin * ymin), ymin
+            x2, y2 = xmin, sqrt(r * r - xmin * xmin)
+            area = (area_arc(x1, y1, x2, y2, r) +
+                    area_triangle(x1, y1, x2, y2, xmin, ymin))
+
+    return area
diff --git a/regions/_geometry/core.pxd b/regions/_geometry/core.pxd
new file mode 100644
index 00000000..d894d08a
--- /dev/null
+++ b/regions/_geometry/core.pxd
@@ -0,0 +1,6 @@
+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)
+cdef double area_arc_unit(double x1, double y1, double x2, double y2)
+cdef int in_triangle(double x, double y, double x1, double y1, double x2, double y2, double x3, double y3)
+cdef double overlap_area_triangle_unit_circle(double x1, double y1, double x2, double y2, double x3, double y3)
diff --git a/regions/_geometry/core.pyx b/regions/_geometry/core.pyx
new file mode 100644
index 00000000..2082eff9
--- /dev/null
+++ b/regions/_geometry/core.pyx
@@ -0,0 +1,377 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+"""The functions here are the core geometry functions."""
+
+from __future__ import (absolute_import, division, print_function,
+                        unicode_literals)
+
+import numpy as np
+cimport numpy as np
+
+
+__all__ = ['elliptical_overlap_grid']
+
+
+cdef extern from "math.h":
+
+    double asin(double x)
+    double sin(double x)
+    double cos(double x)
+    double sqrt(double x)
+    double fabs(double x)
+
+from cpython cimport bool
+
+DTYPE = np.float64
+ctypedef np.float64_t DTYPE_t
+
+cimport cython
+
+
+ctypedef struct point:
+    double x
+    double y
+
+
+ctypedef struct intersections:
+    point p1
+    point p2
+
+
+# NOTE: The following two functions use cdef because they are not intended to be
+# called from the Python code. Using def makes them callable from outside, but
+# also slower. Some functions currently return multiple values, and for those we
+# still use 'def' for now.
+
+
+cdef double distance(double x1, double y1, double x2, double y2):
+    """
+    Distance between two points in two dimensions.
+
+    Parameters
+    ----------
+    x1, y1 : float
+        The coordinates of the first point
+    x2, y2 : float
+        The coordinates of the second point
+
+    Returns
+    -------
+    d : float
+        The Euclidean distance between the two points
+    """
+
+    return sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2)
+
+
+cdef double area_arc(double x1, double y1, double x2, double y2, double r):
+    """
+    Area of a circle arc with radius r between points (x1, y1) and (x2, y2).
+
+    References
+    ----------
+    http://mathworld.wolfram.com/CircularSegment.html
+    """
+
+    cdef double a, theta
+    a = distance(x1, y1, x2, y2)
+    theta = 2. * asin(0.5 * a / r)
+    return 0.5 * r * r * (theta - sin(theta))
+
+
+cdef double area_triangle(double x1, double y1, double x2, double y2, double x3,
+                  double y3):
+    """
+    Area of a triangle defined by three vertices.
+    """
+    return 0.5 * abs(x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2))
+
+
+cdef double area_arc_unit(double x1, double y1, double x2, double y2):
+    """
+    Area of a circle arc with radius R between points (x1, y1) and (x2, y2)
+
+    References
+    ----------
+    http://mathworld.wolfram.com/CircularSegment.html
+    """
+    cdef double a, theta
+    a = distance(x1, y1, x2, y2)
+    theta = 2. * asin(0.5 * a)
+    return 0.5 * (theta - sin(theta))
+
+
+cdef int in_triangle(double x, double y, double x1, double y1, double x2, double y2, double x3, double y3):
+    """
+    Check if a point (x,y) is inside a triangle
+    """
+    cdef int c = 0
+
+    c += ((y1 > y) != (y2 > y) and x < (x2 - x1) * (y - y1) / (y2 - y1) + x1)
+    c += ((y2 > y) != (y3 > y) and x < (x3 - x2) * (y - y2) / (y3 - y2) + x2)
+    c += ((y3 > y) != (y1 > y) and x < (x1 - x3) * (y - y3) / (y1 - y3) + x3)
+
+    return c % 2 == 1
+
+
+cdef intersections circle_line(double x1, double y1, double x2, double y2):
+    """Intersection of a line defined by two points with a unit circle"""
+
+    cdef double a, b, delta, dx, dy
+    cdef double tolerance = 1.e-10
+    cdef intersections inter
+
+    dx = x2 - x1
+    dy = y2 - y1
+
+    if fabs(dx) < tolerance and fabs(dy) < tolerance:
+        inter.p1.x = 2.
+        inter.p1.y = 2.
+        inter.p2.x = 2.
+        inter.p2.y = 2.
+
+    elif fabs(dx) > fabs(dy):
+
+        # Find the slope and intercept of the line
+        a = dy / dx
+        b = y1 - a * x1
+
+        # Find the determinant of the quadratic equation
+        delta = 1. + a * a - b * b
+        if delta > 0.:  # solutions exist
+
+            delta = sqrt(delta)
+
+            inter.p1.x = (- a * b - delta) / (1. + a * a)
+            inter.p1.y = a * inter.p1.x + b
+            inter.p2.x = (- a * b + delta) / (1. + a * a)
+            inter.p2.y = a * inter.p2.x + b
+
+        else:  # no solution, return values > 1
+            inter.p1.x = 2.
+            inter.p1.y = 2.
+            inter.p2.x = 2.
+            inter.p2.y = 2.
+
+    else:
+
+        # Find the slope and intercept of the line
+        a = dx / dy
+        b = x1 - a * y1
+
+        # Find the determinant of the quadratic equation
+        delta = 1. + a * a - b * b
+
+        if delta > 0.:  # solutions exist
+
+            delta = sqrt(delta)
+
+            inter.p1.y = (- a * b - delta) / (1. + a * a)
+            inter.p1.x = a * inter.p1.y + b
+            inter.p2.y = (- a * b + delta) / (1. + a * a)
+            inter.p2.x = a * inter.p2.y + b
+
+        else:  # no solution, return values > 1
+            inter.p1.x = 2.
+            inter.p1.y = 2.
+            inter.p2.x = 2.
+            inter.p2.y = 2.
+
+    return inter
+
+
+cdef point circle_segment_single2(double x1, double y1, double x2, double y2):
+    """
+    The intersection of a line with the unit circle. The intersection the
+    closest to (x2, y2) is chosen.
+    """
+
+    cdef double dx1, dy1, dx2, dy2
+    cdef intersections inter
+    cdef point pt1, pt2, pt
+
+    inter = circle_line(x1, y1, x2, y2)
+
+    pt1 = inter.p1
+    pt2 = inter.p2
+
+    # Can be optimized, but just checking for correctness right now
+    dx1 = fabs(pt1.x - x2)
+    dy1 = fabs(pt1.y - y2)
+    dx2 = fabs(pt2.x - x2)
+    dy2 = fabs(pt2.y - y2)
+
+    if dx1 > dy1:  # compare based on x-axis
+        if dx1 > dx2:
+            pt = pt2
+        else:
+            pt = pt1
+    else:
+        if dy1 > dy2:
+            pt = pt2
+        else:
+            pt = pt1
+
+    return pt
+
+
+cdef intersections circle_segment(double x1, double y1, double x2, double y2):
+    """
+    Intersection(s) of a segment with the unit circle. Discard any
+    solution not on the segment.
+    """
+
+    cdef intersections inter, inter_new
+    cdef point pt1, pt2
+
+    inter = circle_line(x1, y1, x2, y2)
+
+    pt1 = inter.p1
+    pt2 = inter.p2
+
+    if (pt1.x > x1 and pt1.x > x2) or (pt1.x < x1 and pt1.x < x2) or (pt1.y > y1 and pt1.y > y2) or (pt1.y < y1 and pt1.y < y2):
+        pt1.x, pt1.y = 2., 2.
+    if (pt2.x > x1 and pt2.x > x2) or (pt2.x < x1 and pt2.x < x2) or (pt2.y > y1 and pt2.y > y2) or (pt2.y < y1 and pt2.y < y2):
+        pt2.x, pt2.y = 2., 2.
+
+    if pt1.x > 1. and pt2.x < 2.:
+        inter_new.p1 = pt1
+        inter_new.p2 = pt2
+    else:
+        inter_new.p1 = pt2
+        inter_new.p2 = pt1
+
+    return inter_new
+
+
+cdef double overlap_area_triangle_unit_circle(double x1, double y1, double x2, double y2, double x3, double y3):
+    """
+    Given a triangle defined by three points (x1, y1), (x2, y2), and
+    (x3, y3), find the area of overlap with the unit circle.
+    """
+
+    cdef double d1, d2, d3
+    cdef bool in1, in2, in3
+    cdef bool on1, on2, on3
+    cdef double area
+    cdef double PI = np.pi
+    cdef intersections inter
+    cdef point pt1, pt2, pt3, pt4, pt5, pt6, pt_tmp
+
+    # Find distance of all vertices to circle center
+    d1 = x1 * x1 + y1 * y1
+    d2 = x2 * x2 + y2 * y2
+    d3 = x3 * x3 + y3 * y3
+
+    # Order vertices by distance from origin
+    if d1 < d2:
+        if d2 < d3:
+            pass
+        elif d1 < d3:
+            x2, y2, d2, x3, y3, d3 = x3, y3, d3, x2, y2, d2
+        else:
+            x1, y1, d1, x2, y2, d2, x3, y3, d3 = x3, y3, d3, x1, y1, d1, x2, y2, d2
+
+    else:
+        if d1 < d3:
+            x1, y1, d1, x2, y2, d2 = x2, y2, d2, x1, y1, d1
+        elif d2 < d3:
+            x1, y1, d1, x2, y2, d2, x3, y3, d3 = x2, y2, d2, x3, y3, d3, x1, y1, d1
+        else:
+            x1, y1, d1, x2, y2, d2, x3, y3, d3 = x3, y3, d3, x2, y2, d2, x1, y1, d1
+
+    if d1 > d2 or d2 > d3 or d1 > d3:
+        raise Exception("ERROR: vertices did not sort correctly")
+
+    # Determine number of vertices inside circle
+    in1 = d1 < 1
+    in2 = d2 < 1
+    in3 = d3 < 1
+
+    # Determine which vertices are on the circle
+    on1 = fabs(d1 - 1) < 1.e-10
+    on2 = fabs(d2 - 1) < 1.e-10
+    on3 = fabs(d3 - 1) < 1.e-10
+
+    if on3 or in3:  # triangle is completely in circle
+
+        area = area_triangle(x1, y1, x2, y2, x3, y3)
+
+    elif in2 or on2:
+        # If vertex 1 or 2 are on the edge of the circle, then we use the dot
+        # product to vertex 3 to determine whether an intersection takes place.
+        intersect13 = not on1 or x1 * (x3 - x1) + y1 * (y3 - y1) < 0.
+        intersect23 = not on2 or x2 * (x3 - x2) + y2 * (y3 - y2) < 0.
+        if intersect13 and intersect23 and not on2:
+            pt1 = circle_segment_single2(x1, y1, x3, y3)
+            pt2 = circle_segment_single2(x2, y2, x3, y3)
+            area = area_triangle(x1, y1, x2, y2, pt1.x, pt1.y) \
+                 + area_triangle(x2, y2, pt1.x, pt1.y, pt2.x, pt2.y) \
+                 + area_arc_unit(pt1.x, pt1.y, pt2.x, pt2.y)
+        elif intersect13:
+            pt1 = circle_segment_single2(x1, y1, x3, y3)
+            area = area_triangle(x1, y1, x2, y2, pt1.x, pt1.y) \
+                 + area_arc_unit(x2, y2, pt1.x, pt1.y)
+        elif intersect23:
+            pt2 = circle_segment_single2(x2, y2, x3, y3)
+            area = area_triangle(x1, y1, x2, y2, pt2.x, pt2.y) \
+                 + area_arc_unit(x1, y1, pt2.x, pt2.y)
+        else:
+            area = area_arc_unit(x1, y1, x2, y2)
+
+    elif on1:
+        # The triangle is outside the circle
+        area = 0.0
+    elif in1:
+        # Check for intersections of far side with circle
+        inter = circle_segment(x2, y2, x3, y3)
+        pt1 = inter.p1
+        pt2 = inter.p2
+        pt3 = circle_segment_single2(x1, y1, x2, y2)
+        pt4 = circle_segment_single2(x1, y1, x3, y3)
+        if pt1.x > 1.:  # indicates no intersection
+            # Code taken from `sep.h`.
+            # TODO: use `sep` and get rid of this Cython code.
+            if (((0.-pt3.y) * (pt4.x-pt3.x) > (pt4.y-pt3.y) * (0.-pt3.x)) !=
+                ((y1-pt3.y) * (pt4.x-pt3.x) > (pt4.y-pt3.y) * (x1-pt3.x))):
+                area = area_triangle(x1, y1, pt3.x, pt3.y, pt4.x, pt4.y) \
+                     + (PI - area_arc_unit(pt3.x, pt3.y, pt4.x, pt4.y))
+            else:
+                area = area_triangle(x1, y1, pt3.x, pt3.y, pt4.x, pt4.y) \
+                     + area_arc_unit(pt3.x, pt3.y, pt4.x, pt4.y)
+        else:
+            if (pt2.x - x2)**2 + (pt2.y - y2)**2 < (pt1.x - x2)**2 + (pt1.y - y2)**2:
+                pt1, pt2 = pt2, pt1
+            area = area_triangle(x1, y1, pt3.x, pt3.y, pt1.x, pt1.y) \
+                 + area_triangle(x1, y1, pt1.x, pt1.y, pt2.x, pt2.y) \
+                 + area_triangle(x1, y1, pt2.x, pt2.y, pt4.x, pt4.y) \
+                 + area_arc_unit(pt1.x, pt1.y, pt3.x, pt3.y) \
+                 + area_arc_unit(pt2.x, pt2.y, pt4.x, pt4.y)
+    else:
+        inter = circle_segment(x1, y1, x2, y2)
+        pt1 = inter.p1
+        pt2 = inter.p2
+        inter = circle_segment(x2, y2, x3, y3)
+        pt3 = inter.p1
+        pt4 = inter.p2
+        inter = circle_segment(x3, y3, x1, y1)
+        pt5 = inter.p1
+        pt6 = inter.p2
+        if pt1.x <= 1.:
+            xp, yp = 0.5 * (pt1.x + pt2.x), 0.5 * (pt1.y + pt2.y)
+            area = overlap_area_triangle_unit_circle(x1, y1, x3, y3, xp, yp) \
+                 + overlap_area_triangle_unit_circle(x2, y2, x3, y3, xp, yp)
+        elif pt3.x <= 1.:
+            xp, yp = 0.5 * (pt3.x + pt4.x), 0.5 * (pt3.y + pt4.y)
+            area = overlap_area_triangle_unit_circle(x3, y3, x1, y1, xp, yp) \
+                 + overlap_area_triangle_unit_circle(x2, y2, x1, y1, xp, yp)
+        elif pt5.x <= 1.:
+            xp, yp = 0.5 * (pt5.x + pt6.x), 0.5 * (pt5.y + pt6.y)
+            area = overlap_area_triangle_unit_circle(x1, y1, x2, y2, xp, yp) \
+                 + overlap_area_triangle_unit_circle(x3, y3, x2, y2, xp, yp)
+        else:  # no intersections
+            if in_triangle(0., 0., x1, y1, x2, y2, x3, y3):
+                return PI
+            else:
+                return 0.
+
+    return area
diff --git a/regions/_geometry/elliptical_overlap.pyx b/regions/_geometry/elliptical_overlap.pyx
new file mode 100644
index 00000000..4d091982
--- /dev/null
+++ b/regions/_geometry/elliptical_overlap.pyx
@@ -0,0 +1,199 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+"""
+The functions defined here allow one to determine the exact area of
+overlap of an ellipse and a triangle (written by Thomas Robitaille).
+The approach is to divide the rectangle into two triangles, and
+reproject these so that the ellipse is a unit circle, then compute the
+intersection of a triangle with a unit circle.
+"""
+
+from __future__ import (absolute_import, division, print_function,
+                        unicode_literals)
+
+import numpy as np
+cimport numpy as np
+
+
+__all__ = ['elliptical_overlap_grid']
+
+
+cdef extern from "math.h":
+
+    double asin(double x)
+    double sin(double x)
+    double cos(double x)
+    double sqrt(double x)
+
+from cpython cimport bool
+
+DTYPE = np.float64
+ctypedef np.float64_t DTYPE_t
+
+cimport cython
+
+# NOTE: Here we need to make sure we use cimport to import the C functions from
+# core (since these were defined with cdef). This also requires the core.pxd
+# file to exist with the function signatures.
+from .core cimport distance, area_triangle, overlap_area_triangle_unit_circle
+
+
+def elliptical_overlap_grid(double xmin, double xmax, double ymin, double ymax,
+                            int nx, int ny, double rx, double ry, double theta,
+                            int use_exact, int subpixels):
+    """
+    elliptical_overlap_grid(xmin, xmax, ymin, ymax, nx, ny, rx, ry,
+                             use_exact, subpixels)
+
+    Area of overlap between an ellipse and a pixel grid. The ellipse is
+    centered on the origin.
+
+    Parameters
+    ----------
+    xmin, xmax, ymin, ymax : float
+        Extent of the grid in the x and y direction.
+    nx, ny : int
+        Grid dimensions.
+    rx : float
+        The semimajor axis of the ellipse.
+    ry : float
+        The semiminor axis of the ellipse.
+    theta : float
+        The position angle of the semimajor axis in radians (counterclockwise).
+    use_exact : 0 or 1
+        If set to 1, calculates the exact overlap, while if set to 0, uses a
+        subpixel sampling method with ``subpixel`` subpixels in each direction.
+    subpixels : int
+        If ``use_exact`` is 0, each pixel is resampled by this factor in each
+        dimension. Thus, each pixel is divided into ``subpixels ** 2``
+        subpixels.
+
+    Returns
+    -------
+    frac : `~numpy.ndarray`
+        2-d array giving the fraction of the overlap.
+    """
+
+    cdef unsigned int i, j
+    cdef double x, y, dx, dy
+    cdef double bxmin, bxmax, bymin, bymax
+    cdef double pxmin, pxmax, pymin, pymax
+    cdef double norm
+
+    # Define output array
+    cdef np.ndarray[DTYPE_t, ndim=2] frac = np.zeros([ny, nx], dtype=DTYPE)
+
+    # Find the width of each element in x and y
+    dx = (xmax - xmin) / nx
+    dy = (ymax - ymin) / ny
+
+    norm = 1. / (dx * dy)
+
+    # For now we use a bounding circle and then use that to find a bounding box
+    # but of course this is inefficient and could be done better.
+
+    # Find bounding circle radius
+    r = max(rx, ry)
+
+    # Define bounding box
+    bxmin = -r - 0.5 * dx
+    bxmax = +r + 0.5 * dx
+    bymin = -r - 0.5 * dy
+    bymax = +r + 0.5 * dy
+
+    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:
+                    if use_exact:
+                        frac[j, i] = elliptical_overlap_single_exact(
+                            pxmin, pymin, pxmax, pymax, rx, ry, theta) * norm
+                    else:
+                        frac[j, i] = elliptical_overlap_single_subpixel(
+                            pxmin, pymin, pxmax, pymax, rx, ry, theta,
+                            subpixels)
+    return frac
+
+
+# NOTE: The following two functions use cdef because they are not
+# intended to be called from the Python code. Using def makes them
+# callable from outside, but also slower. In any case, these aren't useful
+# to call from outside because they only operate on a single pixel.
+
+
+cdef double elliptical_overlap_single_subpixel(double x0, double y0,
+                                               double x1, double y1,
+                                               double rx, double ry,
+                                               double theta, int subpixels):
+    """
+    Return the fraction of overlap between a ellipse and a single pixel with
+    given extent, using a sub-pixel sampling method.
+    """
+
+    cdef unsigned int i, j
+    cdef double x, y
+    cdef double frac = 0.  # Accumulator.
+    cdef double inv_rx_sq, inv_ry_sq
+    cdef double cos_theta = cos(theta)
+    cdef double sin_theta = sin(theta)
+    cdef double dx, dy
+    cdef double x_tr, y_tr
+
+    dx = (x1 - x0) / subpixels
+    dy = (y1 - y0) / subpixels
+
+    inv_rx_sq = 1. / (rx * rx)
+    inv_ry_sq = 1. / (ry * ry)
+
+    x = x0 - 0.5 * dx
+    for i in range(subpixels):
+        x += dx
+        y = y0 - 0.5 * dy
+        for j in range(subpixels):
+            y += dy
+
+            # Transform into frame of rotated ellipse
+            x_tr = y * sin_theta + x * cos_theta
+            y_tr = y * cos_theta - x * sin_theta
+
+            if x_tr * x_tr * inv_rx_sq + y_tr * y_tr * inv_ry_sq < 1.:
+                frac += 1.
+
+    return frac / (subpixels * subpixels)
+
+
+cdef double elliptical_overlap_single_exact(double xmin, double ymin,
+                                            double xmax, double ymax,
+                                            double rx, double ry,
+                                            double theta):
+    """
+    Given a rectangle defined by (xmin, ymin, xmax, ymax) and an ellipse
+    with major and minor axes rx and ry respectively, position angle theta,
+    and centered at the origin, find the area of overlap.
+    """
+
+    cdef double cos_m_theta = cos(-theta)
+    cdef double sin_m_theta = sin(-theta)
+    cdef double scale
+
+    # Find scale by which the areas will be shrunk
+    scale = rx * ry
+
+    # Reproject rectangle to frame of reference in which ellipse is a
+    # unit circle
+    x1, y1 = ((xmin * cos_m_theta - ymin * sin_m_theta) / rx,
+              (xmin * sin_m_theta + ymin * cos_m_theta) / ry)
+    x2, y2 = ((xmax * cos_m_theta - ymin * sin_m_theta) / rx,
+              (xmax * sin_m_theta + ymin * cos_m_theta) / ry)
+    x3, y3 = ((xmax * cos_m_theta - ymax * sin_m_theta) / rx,
+              (xmax * sin_m_theta + ymax * cos_m_theta) / ry)
+    x4, y4 = ((xmin * cos_m_theta - ymax * sin_m_theta) / rx,
+              (xmin * sin_m_theta + ymax * cos_m_theta) / ry)
+
+    # Divide resulting quadrilateral into two triangles and find
+    # intersection with unit circle
+    return (overlap_area_triangle_unit_circle(x1, y1, x2, y2, x3, y3) +
+            overlap_area_triangle_unit_circle(x1, y1, x4, y4, x3, y3)) * scale
diff --git a/regions/_geometry/rectangular_overlap.pyx b/regions/_geometry/rectangular_overlap.pyx
new file mode 100644
index 00000000..43c493cf
--- /dev/null
+++ b/regions/_geometry/rectangular_overlap.pyx
@@ -0,0 +1,132 @@
+# 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__ = ['rectangular_overlap_grid']
+
+
+cdef extern from "math.h":
+
+    double asin(double x)
+    double sin(double x)
+    double cos(double x)
+    double sqrt(double x)
+    double fabs(double x)
+
+from cpython cimport bool
+
+DTYPE = np.float64
+ctypedef np.float64_t DTYPE_t
+
+cimport cython
+
+
+def rectangular_overlap_grid(double xmin, double xmax, double ymin,
+                             double ymax, int nx, int ny, double width,
+                             double height, double theta, int use_exact,
+                             int subpixels):
+    """
+    rectangular_overlap_grid(xmin, xmax, ymin, ymax, nx, ny, width, height,
+                             use_exact, subpixels)
+
+    Area of overlap between a rectangle and a pixel grid. The rectangle is
+    centered on the origin.
+
+    Parameters
+    ----------
+    xmin, xmax, ymin, ymax : float
+        Extent of the grid in the x and y direction.
+    nx, ny : int
+        Grid dimensions.
+    width : float
+        The width of the rectangle
+    height : float
+        The height of the rectangle
+    theta : float
+        The position angle of the rectangle in radians (counterclockwise).
+    use_exact : 0 or 1
+        If set to 1, calculates the exact overlap, while if set to 0, uses a
+        subpixel sampling method with ``subpixel`` subpixels in each direction.
+    subpixels : int
+        If ``use_exact`` is 0, each pixel is resampled by this factor in each
+        dimension. Thus, each pixel is divided into ``subpixels ** 2``
+        subpixels.
+
+    Returns
+    -------
+    frac : `~numpy.ndarray`
+        2-d array giving the fraction of the overlap.
+    """
+
+    cdef unsigned int i, j
+    cdef double x, y, dx, dy
+    cdef double pxmin, pxmax, pymin, 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 "
+                                  "rectangular apertures")
+
+    # Find the width of each element in x and y
+    dx = (xmax - xmin) / nx
+    dy = (ymax - ymin) / ny
+
+    # TODO: can implement a bounding box here for efficiency (as for the
+    # circular and elliptical aperture photometry)
+
+    for i in range(nx):
+        pxmin = xmin + i * dx  # lower end of pixel
+        pxmax = pxmin + dx  # upper end of pixel
+        for j in range(ny):
+            pymin = ymin + j * dy
+            pymax = pymin + dy
+            frac[j, i] = rectangular_overlap_single_subpixel(
+                pxmin, pymin, pxmax, pymax, width, height, theta,
+                subpixels)
+
+    return frac
+
+
+cdef double rectangular_overlap_single_subpixel(double x0, double y0,
+                                                double x1, double y1,
+                                                double width, double height,
+                                                double theta, int subpixels):
+    """
+    Return the fraction of overlap between a rectangle and a single pixel with
+    given extent, using a sub-pixel sampling method.
+    """
+
+    cdef unsigned int i, j
+    cdef double x, y
+    cdef double frac = 0.  # Accumulator.
+    cdef double cos_theta = cos(theta)
+    cdef double sin_theta = sin(theta)
+    cdef double half_width, half_height
+
+    half_width = width / 2.
+    half_height = height / 2.
+
+    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
+
+            # Transform into frame of rotated rectangle
+            x_tr = y * sin_theta + x * cos_theta
+            y_tr = y * cos_theta - x * sin_theta
+
+            if fabs(x_tr) < half_width and fabs(y_tr) < half_height:
+                frac += 1.
+
+    return frac / (subpixels * subpixels)
diff --git a/regions/_geometry/tests/__init__.py b/regions/_geometry/tests/__init__.py
new file mode 100644
index 00000000..356c0a05
--- /dev/null
+++ b/regions/_geometry/tests/__init__.py
@@ -0,0 +1,4 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+"""
+This package contains affiliated package tests.
+"""
diff --git a/regions/_geometry/tests/test_circular_overlap_grid.py b/regions/_geometry/tests/test_circular_overlap_grid.py
new file mode 100644
index 00000000..76bc4240
--- /dev/null
+++ b/regions/_geometry/tests/test_circular_overlap_grid.py
@@ -0,0 +1,30 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+from __future__ import (absolute_import, division, print_function,
+                        unicode_literals)
+import itertools
+
+from numpy.testing import assert_allclose
+from astropy.tests.helper import pytest
+
+from .. import circular_overlap_grid
+
+
+grid_sizes = [50, 500, 1000]
+circ_sizes = [0.2, 0.4, 0.8]
+use_exact = [0, 1]
+subsamples = [1, 5, 10]
+arg_list = ['grid_size', 'circ_size', 'use_exact', 'subsample']
+
+
+@pytest.mark.parametrize(('grid_size', 'circ_size', 'use_exact', 'subsample'),
+                         list(itertools.product(grid_sizes, circ_sizes,
+                                                use_exact, subsamples)))
+def test_circular_overlap_grid(grid_size, circ_size, use_exact, subsample):
+    """
+    Test normalization of the overlap grid to make sure that a fully
+    enclosed pixel has a value of 1.0.
+    """
+
+    g = circular_overlap_grid(-1.0, 1.0, -1.0, 1.0, grid_size, grid_size,
+                              circ_size, use_exact, subsample)
+    assert_allclose(g.max(), 1.0)
diff --git a/regions/_geometry/tests/test_elliptical_overlap_grid.py b/regions/_geometry/tests/test_elliptical_overlap_grid.py
new file mode 100644
index 00000000..e514d6d4
--- /dev/null
+++ b/regions/_geometry/tests/test_elliptical_overlap_grid.py
@@ -0,0 +1,37 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+from __future__ import (absolute_import, division, print_function,
+                        unicode_literals)
+import itertools
+
+from numpy.testing import assert_allclose
+from astropy.tests.helper import pytest
+
+from .. import elliptical_overlap_grid
+
+
+grid_sizes = [50, 500, 1000]
+maj_sizes = [0.2, 0.4, 0.8]
+min_sizes = [0.2, 0.4, 0.8]
+angles = [0.0, 0.5, 1.0]
+use_exact = [0, 1]
+subsamples = [1, 5, 10]
+arg_list = ['grid_size', 'maj_size', 'min_size', 'angle', 'use_exact',
+            'subsample']
+
+
+@pytest.mark.parametrize(('grid_size', 'maj_size', 'min_size', 'angle',
+                          'use_exact', 'subsample'),
+                         list(itertools.product(grid_sizes, maj_sizes,
+                                                min_sizes, angles, use_exact,
+                                                subsamples)))
+def test_elliptical_overlap_grid(grid_size, maj_size, min_size, angle,
+                                 use_exact, subsample):
+    """
+    Test normalization of the overlap grid to make sure that a fully
+    enclosed pixel has a value of 1.0.
+    """
+
+    g = elliptical_overlap_grid(-1.0, 1.0, -1.0, 1.0, grid_size, grid_size,
+                                maj_size, min_size, angle, use_exact,
+                                subsample)
+    assert_allclose(g.max(), 1.0)
diff --git a/regions/_geometry/tests/test_rectangular_overlap_grid.py b/regions/_geometry/tests/test_rectangular_overlap_grid.py
new file mode 100644
index 00000000..a3dca94e
--- /dev/null
+++ b/regions/_geometry/tests/test_rectangular_overlap_grid.py
@@ -0,0 +1,30 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+from __future__ import (absolute_import, division, print_function,
+                        unicode_literals)
+import itertools
+
+from numpy.testing import assert_allclose
+from astropy.tests.helper import pytest
+
+from .. import rectangular_overlap_grid
+
+
+grid_sizes = [50, 500, 1000]
+rect_sizes = [0.2, 0.4, 0.8]
+angles = [0.0, 0.5, 1.0]
+subsamples = [1, 5, 10]
+arg_list = ['grid_size', 'rect_size', 'angle', 'subsample']
+
+
+@pytest.mark.parametrize(('grid_size', 'rect_size', 'angle', 'subsample'),
+                         list(itertools.product(grid_sizes, rect_sizes,
+                                                angles, subsamples)))
+def test_rectangular_overlap_grid(grid_size, rect_size, angle, subsample):
+    """
+    Test normalization of the overlap grid to make sure that a fully
+    enclosed pixel has a value of 1.0.
+    """
+
+    g = rectangular_overlap_grid(-1.0, 1.0, -1.0, 1.0, grid_size, grid_size,
+                                 rect_size, rect_size, angle, 0, subsample)
+    assert_allclose(g.max(), 1.0)
diff --git a/regions/core/__init__.py b/regions/core/__init__.py
index 619d08a4..892a5c67 100644
--- a/regions/core/__init__.py
+++ b/regions/core/__init__.py
@@ -4,3 +4,5 @@
 from .core import *
 from .pixcoord import *
 from .compound import *
+from .mask import *
+from .bounding_box import *
diff --git a/regions/core/bounding_box.py b/regions/core/bounding_box.py
new file mode 100644
index 00000000..61cae1a8
--- /dev/null
+++ b/regions/core/bounding_box.py
@@ -0,0 +1,50 @@
+class BoundingBox(object):
+    """
+    A simple class used to represent a bounding box in integer pixel
+    coordinates.
+
+    Parameters
+    ----------
+    ixmin, ixmax, iymin, iymax : int
+        The bounding indices. Note that the upper values (iymax and ixmax) are
+        exlusive as for normal slices in Python.
+    """
+
+    def __init__(self, ixmin, ixmax, iymin, iymax):
+        self.ixmin = ixmin
+        self.ixmax = ixmax
+        self.iymin = iymin
+        self.iymax = iymax
+
+    @property
+    def shape(self):
+        """
+        The shape of the bounding box
+        """
+        return self.iymax - self.iymin, self.ixmax - self.ixmin
+
+    @property
+    def slices(self):
+        """
+        The bounding box as a pair of `slice` objects.
+        """
+        return (slice(self.iymin, self.iymax),
+                slice(self.ixmin, self.ixmax))
+
+    @property
+    def extent(self):
+        """
+        The 'extent' of the mask, i.e the bounding box from the bottom left
+        corner of the lower left pixel to the upper right corner of the upper
+        right pixel. This can be used for example when plotting using Matplotlib.
+        """
+        return self.ixmin - 0.5, self.ixmax - 0.5, self.iymin - 0.5, self.iymax - 0.5
+
+    def as_patch(self, **kwargs):
+        """
+        Return a Matplotlib patch that represents the bounding box.
+        """
+        from matplotlib.patches import Rectangle
+        return Rectangle((self.ixmin - 0.5, self.iymin - 0.5),
+                         self.ixmax - self.ixmin,
+                         self.iymax - self.iymin, **kwargs)
diff --git a/regions/core/core.py b/regions/core/core.py
index 2a24f54a..859aaed0 100644
--- a/regions/core/core.py
+++ b/regions/core/core.py
@@ -1,17 +1,22 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
-from __future__ import absolute_import, division, print_function, unicode_literals
+from __future__ import (absolute_import, division, print_function,
+                        unicode_literals)
 import abc
 import operator
+
 from astropy.extern import six
 
+
 __all__ = ['Region', 'PixelRegion', 'SkyRegion']
 
+
 """
 Here we define global variables for the default `origin` and `mode` used
 for WCS transformations throughout the `regions` package.
 
 Their purpose is to simplify achieving uniformity across the codebase.
-They are mainly used as default arguments for methods that do WCS transformations.
+They are mainly used as default arguments for methods that do WCS
+transformations.
 
 They are private (with an underscore), not part of the public API,
 users should not touch them.
@@ -19,6 +24,8 @@
 _DEFAULT_WCS_ORIGIN = 0
 _DEFAULT_WCS_MODE = 'all'
 
+VALID_MASK_MODES = set(['center', 'exact', 'subpixels'])
+
 
 @six.add_metaclass(abc.ABCMeta)
 class Region(object):
@@ -145,6 +152,14 @@ def to_sky(self, wcs, mode='local', tolerance=None):
             The tolerance for the ``'full'`` mode described above.
         """
 
+    @abc.abstractproperty
+    def bounding_box(self):
+        """
+        The minimal bounding box (in integer pixel coordinates) that contains
+        the region.
+        """
+        pass
+
     @abc.abstractmethod
     def to_mask(self, mode='center', subpixels=5):
         """
@@ -152,14 +167,10 @@ def to_mask(self, mode='center', subpixels=5):
 
         Parameters
         ----------
-        mode : { 'center' | 'any' | 'all' | 'exact' | 'subpixels'}
+        mode : { 'center' | 'exact' | 'subpixels'}
             The following modes are available:
                 * ``'center'``: returns 1 for pixels where the center is in
                   the region, and 0 otherwise.
-                * ``'any'``: returns 1 for pixels where any of the pixel is
-                  in the region, and 0 otherwise.
-                * ``'all'``: returns 1 for pixels that are completely inside
-                  the region, 0 otherwise.
                 * ``'exact'``: returns a value between 0 and 1 giving the
                   fractional level of overlap of the pixel with the region.
                 * ``'subpixels'``: A pixel is divided into subpixels and
@@ -176,12 +187,19 @@ def to_mask(self, mode='center', subpixels=5):
 
         Returns
         -------
-        mask : `~numpy.ndarray`
-            A mask indicating whether each pixel is contained in the region.
-        slice_x, slice_y : `slice`
-            Slices for x and y which can be used on an array to extract the
-            same region as the mask.
-        """
+        mask : `~regions.Mask`
+            A region mask object.
+        """
+
+    @staticmethod
+    def _validate_mode(mode, subpixels):
+        if mode not in VALID_MASK_MODES:
+            raise ValueError("Invalid mask mode: {0} (should be one "
+                             "of {1})".format(mode, '/'.join(VALID_MASK_MODES)))
+        if mode == 'subpixels':
+            if not isinstance(subpixels, int) or subpixels <= 0:
+                raise ValueError("Invalid subpixels value: {0} (should be"
+                                 " a strictly positive integer)".format(subpixels))
 
     @abc.abstractmethod
     def to_shapely(self):
diff --git a/regions/core/mask.py b/regions/core/mask.py
new file mode 100644
index 00000000..b328a268
--- /dev/null
+++ b/regions/core/mask.py
@@ -0,0 +1,198 @@
+import numpy as np
+from astropy import units as u
+
+__all__ = ['Mask']
+
+
+class Mask(object):
+    """
+    Class for a region mask.
+
+    Parameters
+    ----------
+    mask : array_like
+        A 2D array of a region mask representing the fractional overlap
+        of the region on the pixel grid.  This should be the full-sized
+        (i.e. not truncated) array that is the direct output of one of
+        the low-level "geometry" functions.
+
+    bbox : `regions.BoundingBox`
+        The bounding box for the region.
+
+    Examples
+    --------
+
+    Usage examples are provided in the :ref:`gs-masks` section of the docs.
+    """
+
+    def __init__(self, mask, bbox):
+        if mask.shape != bbox.shape:
+            raise ValueError("shape of mask and bounding box should match")
+        self.data = np.asanyarray(mask)
+        self.bbox = bbox
+
+    @property
+    def shape(self):
+        """
+        The shape of the mask array.
+        """
+        return self.data.shape
+
+    def __array__(self):
+        """
+        Array representation of the mask array (e.g., for matplotlib).
+        """
+
+        return self.data
+
+    def _overlap_slices(self, shape):
+        """
+        Calculate the slices for the overlapping part of the bounding box
+        and an array of the given shape.
+
+        Parameters
+        ----------
+        shape : tuple of int
+            The ``(ny, nx)`` shape of array where the slices are to be
+            applied.
+
+        Returns
+        -------
+        slices_large : tuple of slices
+            A tuple of slice objects for each axis of the large array,
+            such that ``large_array[slices_large]`` extracts the region
+            of the large array that overlaps with the small array.
+
+        slices_small : slice
+            A tuple of slice objects for each axis of the small array,
+            such that ``small_array[slices_small]`` extracts the region
+            of the small array that is inside the large array.
+        """
+
+        if len(shape) != 2:
+            raise ValueError('input shape must have 2 elements.')
+
+        xmin = self.bbox.ixmin
+        xmax = self.bbox.ixmax
+        ymin = self.bbox.iymin
+        ymax = self.bbox.iymax
+
+        if (xmin >= shape[1] or ymin >= shape[0] or xmax <= 0 or ymax <= 0):
+            # no overlap of the region with the data
+            return None, None
+
+        slices_large = (slice(max(ymin, 0), min(ymax, shape[0])),
+                        slice(max(xmin, 0), min(xmax, shape[1])))
+
+        slices_small = (slice(max(-ymin, 0),
+                              min(ymax - ymin, shape[0] - ymin)),
+                        slice(max(-xmin, 0),
+                              min(xmax - xmin, shape[1] - xmin)))
+
+        return slices_large, slices_small
+
+    def to_image(self, shape):
+        """
+        Return an image of the mask in a 2D array of the given shape,
+        taking any edge effects into account.
+
+        Parameters
+        ----------
+        shape : tuple of int
+            The ``(ny, nx)`` shape of the output array.
+
+        Returns
+        -------
+        result : `~numpy.ndarray`
+            A 2D array of the mask.
+        """
+
+        if len(shape) != 2:
+            raise ValueError('input shape must have 2 elements.')
+
+        mask = np.zeros(shape)
+
+        try:
+            mask[self.bbox.slices] = self.data
+        except ValueError:    # partial or no overlap
+            slices_large, slices_small = self._overlap_slices(shape)
+
+            if slices_small is None:
+                return None    # no overlap
+
+            mask = np.zeros(shape)
+            mask[slices_large] = self.data[slices_small]
+
+        return mask
+
+    def cutout(self, data, fill_value=0.):
+        """
+        Create a cutout from the input data over the mask bounding box,
+        taking any edge effects into account.
+
+        Parameters
+        ----------
+        data : array_like or `~astropy.units.Quantity`
+            A 2D array on which to apply the region mask.
+
+        fill_value : float, optional
+            The value is used to fill pixels where the region mask
+            does not overlap with the input ``data``.  The default is 0.
+
+        Returns
+        -------
+        result : `~numpy.ndarray`
+            A 2D array cut out from the input ``data`` representing the
+            same cutout region as the region mask.  If there is a
+            partial overlap of the region mask with the input data,
+            pixels outside of the data will be assigned to
+            ``fill_value``.  `None` is returned if there is no overlap
+            of the region with the input ``data``.
+        """
+
+        data = np.asanyarray(data)
+        cutout = data[self.bbox.slices]
+
+        if cutout.shape != self.shape:
+            slices_large, slices_small = self._overlap_slices(data.shape)
+
+            if slices_small is None:
+                return None    # no overlap
+
+            cutout = np.zeros(self.shape, dtype=data.dtype)
+            cutout[:] = fill_value
+            cutout[slices_small] = data[slices_large]
+
+            if isinstance(data, u.Quantity):
+                cutout = u.Quantity(cutout, unit=data.unit)
+
+        return cutout
+
+    def multiply(self, data, fill_value=0.):
+        """
+        Apply the region mask to the input data, taking any edge effects
+        into account.
+
+        The result is a mask-weighted cutout from the data.
+
+        Parameters
+        ----------
+        data : array_like or `~astropy.units.Quantity`
+            A 2D array on which to apply the region mask.
+
+        fill_value : float, optional
+            The value is used to fill pixels where the region mask does
+            not overlap with the input ``data``.  The default is 0.
+
+        Returns
+        -------
+        result : `~numpy.ndarray`
+            A 2D mask-weighted cutout from the input ``data``.  If there
+            is a partial overlap of the region mask with the input data,
+            pixels outside of the data will be assigned to
+            ``fill_value`` before being multipled with the mask.  `None`
+            is returned if there is no overlap of the region with the
+            input ``data``.
+        """
+
+        return self.cutout(data, fill_value=fill_value) * self.data
diff --git a/regions/shapes/circle.py b/regions/shapes/circle.py
index eae2f8cd..8b7aff6d 100644
--- a/regions/shapes/circle.py
+++ b/regions/shapes/circle.py
@@ -1,10 +1,16 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 from __future__ import absolute_import, division, print_function, unicode_literals
+
 import math
+
+import numpy as np
 from astropy.coordinates import Angle
 from astropy.wcs.utils import pixel_to_skycoord
-from ..core import PixelRegion, SkyRegion
+
+from ..core import PixelRegion, SkyRegion, Mask, BoundingBox
+from ..utils import skycoord_to_pixel_scale_angle
 from ..utils.wcs_helpers import skycoord_to_pixel_scale_angle
+from .._geometry import circular_overlap_grid
 
 __all__ = ['CirclePixelRegion', 'CircleSkyRegion']
 
@@ -61,9 +67,55 @@ def to_sky(self, wcs, mode='local', tolerance=None):
         radius = Angle(self.radius / scale, 'deg')
         return CircleSkyRegion(center, radius)
 
-    def to_mask(self, mode='center'):
-        # TODO: needs to be implemented
-        raise NotImplementedError
+    @property
+    def bounding_box(self):
+
+        # Find exact bounds
+        xmin = self.center.x - self.radius
+        xmax = self.center.x + self.radius
+        ymin = self.center.y - self.radius
+        ymax = self.center.y + self.radius
+
+        # Find range of pixels. We use round here because if the region extends
+        # to e.g. -0.4, it's enough to set the bounding box lower value to 0
+        # because the 0-th pixel goes from -0.5 to 0.5. At the upper end we add
+        # 1 because the upper limits need to be exlcusive.
+        ixmin = round(xmin)
+        ixmax = round(xmax) + 1
+        iymin = round(ymin)
+        iymax = round(ymax) + 1
+
+        return BoundingBox(ixmin, ixmax, iymin, iymax)
+
+    def to_mask(self, mode='center', subpixels=1):
+
+        # NOTE: assumes this class represents a single circle
+
+        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 - self.center.x
+        xmax = float(bbox.ixmax) - 0.5 - self.center.x
+        ymin = float(bbox.iymin) - 0.5 - self.center.y
+        ymax = float(bbox.iymax) - 0.5 - self.center.y
+
+        if mode == 'subpixels':
+            use_exact = 0
+        else:
+            use_exact = 1
+
+        fraction = circular_overlap_grid(xmin, xmax, ymin, ymax, nx, ny,
+                                         self.radius, use_exact, subpixels)
+
+        return Mask(fraction, bbox=bbox)
 
     def as_patch(self, **kwargs):
         import matplotlib.patches as mpatches
diff --git a/regions/shapes/ellipse.py b/regions/shapes/ellipse.py
index 3f80f07e..8e4d6fa8 100644
--- a/regions/shapes/ellipse.py
+++ b/regions/shapes/ellipse.py
@@ -2,7 +2,8 @@
 from __future__ import absolute_import, division, print_function, unicode_literals
 import math
 from astropy import units as u
-from ..core import PixelRegion, SkyRegion
+from ..core import PixelRegion, SkyRegion, Mask, BoundingBox
+from .._geometry import elliptical_overlap_grid
 
 __all__ = ['EllipsePixelRegion', 'EllipseSkyRegion']
 
@@ -62,9 +63,58 @@ def to_sky(self, wcs, mode='local', tolerance=None):
         # TODO: needs to be implemented
         raise NotImplementedError
 
-    def to_mask(self, mode='center'):
-        # TODO: needs to be implemented
-        raise NotImplementedError
+    @property
+    def bounding_box(self):
+
+        # Find exact bounds
+        # FIXME: this is not the minimal bounding box, and can be optimized
+        xmin = self.center.x - max(self.major, self.minor)
+        xmax = self.center.x + max(self.major, self.minor)
+        ymin = self.center.y - max(self.major, self.minor)
+        ymax = self.center.y + max(self.major, self.minor)
+
+        # Find range of pixels. We use round here because if the region extends
+        # to e.g. -0.4, it's enough to set the bounding box lower value to 0
+        # because the 0-th pixel goes from -0.5 to 0.5. At the upper end we add
+        # 1 because the upper limits need to be exlcusive.
+        ixmin = round(xmin)
+        ixmax = round(xmax) + 1
+        iymin = round(ymin)
+        iymax = round(ymax) + 1
+
+        return BoundingBox(ixmin, ixmax, iymin, iymax)
+
+    def to_mask(self, mode='center', subpixels=5):
+
+        # NOTE: assumes this class represents a single circle
+
+        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 ellipse is at origin
+        xmin = float(bbox.ixmin) - 0.5 - self.center.x
+        xmax = float(bbox.ixmax) - 0.5 - self.center.x
+        ymin = float(bbox.iymin) - 0.5 - self.center.y
+        ymax = float(bbox.iymax) - 0.5 - self.center.y
+
+        if mode == 'subpixels':
+            use_exact = 0
+        else:
+            use_exact = 1
+
+        fraction = elliptical_overlap_grid(xmin, xmax, ymin, ymax, nx, ny,
+                                           self.major, self.minor,
+                                           self.angle.to(u.deg).value,
+                                           use_exact, subpixels)
+
+        return Mask(fraction, bbox=bbox)
 
     def as_patch(self, **kwargs):
         # TODO: needs to be implemented
diff --git a/regions/shapes/point.py b/regions/shapes/point.py
index 67f99eab..b5fc01f6 100644
--- a/regions/shapes/point.py
+++ b/regions/shapes/point.py
@@ -43,6 +43,11 @@ def to_sky(self, wcs, mode='local', tolerance=None):
         # TODO: needs to be implemented
         raise NotImplementedError
 
+    @property
+    def bounding_box(self, mode='center'):
+        # TODO: needs to be implemented
+        raise NotImplementedError
+
     def to_mask(self, mode='center'):
         # TODO: needs to be implemented
         raise NotImplementedError
diff --git a/regions/shapes/polygon.py b/regions/shapes/polygon.py
index 62d6e6ee..7d840206 100644
--- a/regions/shapes/polygon.py
+++ b/regions/shapes/polygon.py
@@ -46,6 +46,11 @@ def to_sky(self, wcs, mode='local', tolerance=None):
         # TODO: needs to be implemented
         raise NotImplementedError
 
+    @property
+    def bounding_box(self, mode='center'):
+        # TODO: needs to be implemented
+        raise NotImplementedError
+
     def to_mask(self, mode='center'):
         # TODO: needs to be implemented
         raise NotImplementedError
diff --git a/regions/shapes/rectangle.py b/regions/shapes/rectangle.py
index 45a694bf..7b853701 100644
--- a/regions/shapes/rectangle.py
+++ b/regions/shapes/rectangle.py
@@ -1,7 +1,11 @@
 # 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
 from astropy import units as u
-from ..core import PixelRegion, SkyRegion
+
+from ..core import PixelRegion, SkyRegion, Mask, BoundingBox
+from .._geometry import rectangular_overlap_grid
 
 __all__ = ['RectanglePixelRegion', 'RectangleSkyRegion']
 
@@ -59,9 +63,59 @@ def to_sky(self, wcs, mode='local', tolerance=None):
         # TODO: needs to be implemented
         raise NotImplementedError
 
-    def to_mask(self, mode='center'):
-        # TODO: needs to be implemented
-        raise NotImplementedError
+    @property
+    def bounding_box(self):
+
+        # Find exact bounds
+        # FIXME: this is not the minimal bounding box, and can be optimized
+        radius = np.hypot(self.width / 2, self.height / 2)
+        xmin = self.center.x - radius
+        xmax = self.center.x + radius
+        ymin = self.center.y - radius
+        ymax = self.center.y + radius
+
+        # Find range of pixels. We use round here because if the region extends
+        # to e.g. -0.4, it's enough to set the bounding box lower value to 0
+        # because the 0-th pixel goes from -0.5 to 0.5. At the upper end we add
+        # 1 because the upper limits need to be exlcusive.
+        ixmin = round(xmin)
+        ixmax = round(xmax) + 1
+        iymin = round(ymin)
+        iymax = round(ymax) + 1
+
+        return BoundingBox(ixmin, ixmax, iymin, iymax)
+
+    def to_mask(self, mode='center', subpixels=5):
+
+        # NOTE: assumes this class represents a single circle
+
+        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 - self.center.x
+        xmax = float(bbox.ixmax) - 0.5 - self.center.x
+        ymin = float(bbox.iymin) - 0.5 - self.center.y
+        ymax = float(bbox.iymax) - 0.5 - self.center.y
+
+        if mode == 'subpixels':
+            use_exact = 0
+        else:
+            use_exact = 1
+
+        fraction = rectangular_overlap_grid(xmin, xmax, ymin, ymax, nx, ny,
+                                            self.width, self.height,
+                                            self.angle.to(u.deg).value,
+                                            use_exact, subpixels)
+
+        return Mask(fraction, bbox=bbox)
 
     def as_patch(self, **kwargs):
         # TODO: needs to be implemented
diff --git a/regions/shapes/setup_package.py b/regions/shapes/setup_package.py
index 48cd980f..0e05da1f 100644
--- a/regions/shapes/setup_package.py
+++ b/regions/shapes/setup_package.py
@@ -1,3 +1,6 @@
+import os
+
 def get_package_data():
-    shape_test = ['data/*.fits']
-    return {'regions.shapes.tests': shape_test}
+    data = [os.path.join('data', '*.fits'),
+            os.path.join('baseline', '*.fits')]
+    return {'regions.shapes.tests': data}
diff --git a/regions/shapes/tests/reference/test_to_mask_circ-mode_center.txt b/regions/shapes/tests/reference/test_to_mask_circ-mode_center.txt
new file mode 100644
index 00000000..5b64916b
--- /dev/null
+++ b/regions/shapes/tests/reference/test_to_mask_circ-mode_center.txt
@@ -0,0 +1,7 @@
+0.00000000e+00 0.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 1.00000000e+00 0.00000000e+00
+1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00
+1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00
+1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.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 0.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00
diff --git a/regions/shapes/tests/reference/test_to_mask_circ-mode_exact.txt b/regions/shapes/tests/reference/test_to_mask_circ-mode_exact.txt
new file mode 100644
index 00000000..e768cb15
--- /dev/null
+++ b/regions/shapes/tests/reference/test_to_mask_circ-mode_exact.txt
@@ -0,0 +1,7 @@
+0.00000000e+00 1.12637824e-01 5.47795293e-01 6.97159714e-01 5.36338125e-01 9.98826171e-02 0.00000000e+00
+1.41612171e-01 9.21332849e-01 1.00000000e+00 1.00000000e+00 1.00000000e+00 9.06385082e-01 1.20533937e-01
+6.46660408e-01 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 6.10634408e-01
+8.43986340e-01 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 8.07960340e-01
+7.30305012e-01 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 6.94279012e-01
+2.73584439e-01 9.91864462e-01 1.00000000e+00 1.00000000e+00 1.00000000e+00 9.86382743e-01 2.43040158e-01
+0.00000000e+00 3.04862210e-01 8.10551293e-01 9.59915714e-01 7.99094125e-01 2.82640956e-01 0.00000000e+00
diff --git a/regions/shapes/tests/reference/test_to_mask_circ-mode_subpixels-subpixels_1.txt b/regions/shapes/tests/reference/test_to_mask_circ-mode_subpixels-subpixels_1.txt
new file mode 100644
index 00000000..5b64916b
--- /dev/null
+++ b/regions/shapes/tests/reference/test_to_mask_circ-mode_subpixels-subpixels_1.txt
@@ -0,0 +1,7 @@
+0.00000000e+00 0.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 1.00000000e+00 0.00000000e+00
+1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00
+1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00
+1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.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 0.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00
diff --git a/regions/shapes/tests/reference/test_to_mask_circ-mode_subpixels-subpixels_18.txt b/regions/shapes/tests/reference/test_to_mask_circ-mode_subpixels-subpixels_18.txt
new file mode 100644
index 00000000..c723a59a
--- /dev/null
+++ b/regions/shapes/tests/reference/test_to_mask_circ-mode_subpixels-subpixels_18.txt
@@ -0,0 +1,7 @@
+0.00000000e+00 1.11111111e-01 5.52469136e-01 7.00617284e-01 5.33950617e-01 9.87654321e-02 0.00000000e+00
+1.41975309e-01 9.16666667e-01 1.00000000e+00 1.00000000e+00 1.00000000e+00 9.10493827e-01 1.20370370e-01
+6.48148148e-01 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 6.04938272e-01
+8.30246914e-01 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 8.17901235e-01
+7.31481481e-01 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 6.94444444e-01
+2.71604938e-01 9.90740741e-01 1.00000000e+00 1.00000000e+00 1.00000000e+00 9.87654321e-01 2.40740741e-01
+0.00000000e+00 3.02469136e-01 8.11728395e-01 9.47530864e-01 7.99382716e-01 2.83950617e-01 0.00000000e+00
diff --git a/regions/shapes/tests/reference/test_to_mask_circ-mode_subpixels-subpixels_5.txt b/regions/shapes/tests/reference/test_to_mask_circ-mode_subpixels-subpixels_5.txt
new file mode 100644
index 00000000..d8485035
--- /dev/null
+++ b/regions/shapes/tests/reference/test_to_mask_circ-mode_subpixels-subpixels_5.txt
@@ -0,0 +1,7 @@
+0.00000000e+00 1.20000000e-01 5.20000000e-01 7.20000000e-01 5.20000000e-01 8.00000000e-02 0.00000000e+00
+1.20000000e-01 9.20000000e-01 1.00000000e+00 1.00000000e+00 1.00000000e+00 8.80000000e-01 1.20000000e-01
+6.40000000e-01 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 6.00000000e-01
+8.00000000e-01 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 8.00000000e-01
+7.20000000e-01 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 7.20000000e-01
+2.80000000e-01 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 2.40000000e-01
+0.00000000e+00 3.20000000e-01 8.00000000e-01 1.00000000e+00 8.00000000e-01 2.80000000e-01 0.00000000e+00
diff --git a/regions/shapes/tests/reference/test_to_mask_elli-mode_center.txt b/regions/shapes/tests/reference/test_to_mask_elli-mode_center.txt
new file mode 100644
index 00000000..70d4ef2b
--- /dev/null
+++ b/regions/shapes/tests/reference/test_to_mask_elli-mode_center.txt
@@ -0,0 +1,7 @@
+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 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 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 0.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
diff --git a/regions/shapes/tests/reference/test_to_mask_elli-mode_exact.txt b/regions/shapes/tests/reference/test_to_mask_elli-mode_exact.txt
new file mode 100644
index 00000000..dd44e405
--- /dev/null
+++ b/regions/shapes/tests/reference/test_to_mask_elli-mode_exact.txt
@@ -0,0 +1,7 @@
+0.00000000e+00 0.00000000e+00 0.00000000e+00 2.01450892e-01 3.89520765e-01 9.56735439e-02 0.00000000e+00
+0.00000000e+00 3.33728284e-03 5.68850724e-01 9.98666422e-01 1.00000000e+00 8.45828552e-01 4.48384112e-03
+0.00000000e+00 4.19946985e-01 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 8.22061396e-02
+6.16413600e-03 8.97539339e-01 1.00000000e+00 1.00000000e+00 1.00000000e+00 9.37056728e-01 1.17070928e-02
+1.11713638e-01 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 5.39575174e-01 0.00000000e+00
+4.03192933e-02 9.57604924e-01 1.00000000e+00 1.00000000e+00 7.70659026e-01 3.79279177e-02 0.00000000e+00
+0.00000000e+00 2.85508104e-01 6.53832220e-01 4.47790886e-01 2.92273877e-02 0.00000000e+00 0.00000000e+00
diff --git a/regions/shapes/tests/reference/test_to_mask_elli-mode_subpixels-subpixels_1.txt b/regions/shapes/tests/reference/test_to_mask_elli-mode_subpixels-subpixels_1.txt
new file mode 100644
index 00000000..70d4ef2b
--- /dev/null
+++ b/regions/shapes/tests/reference/test_to_mask_elli-mode_subpixels-subpixels_1.txt
@@ -0,0 +1,7 @@
+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 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 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 0.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
diff --git a/regions/shapes/tests/reference/test_to_mask_elli-mode_subpixels-subpixels_18.txt b/regions/shapes/tests/reference/test_to_mask_elli-mode_subpixels-subpixels_18.txt
new file mode 100644
index 00000000..46541ba5
--- /dev/null
+++ b/regions/shapes/tests/reference/test_to_mask_elli-mode_subpixels-subpixels_18.txt
@@ -0,0 +1,7 @@
+0.00000000e+00 0.00000000e+00 0.00000000e+00 2.03703704e-01 3.79629630e-01 9.87654321e-02 0.00000000e+00
+0.00000000e+00 3.08641975e-03 5.58641975e-01 1.00000000e+00 1.00000000e+00 8.51851852e-01 3.08641975e-03
+0.00000000e+00 4.19753086e-01 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 8.64197531e-02
+6.17283951e-03 8.98148148e-01 1.00000000e+00 1.00000000e+00 1.00000000e+00 9.38271605e-01 1.23456790e-02
+1.01851852e-01 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 5.43209877e-01 0.00000000e+00
+4.32098765e-02 9.56790123e-01 1.00000000e+00 1.00000000e+00 7.68518519e-01 3.70370370e-02 0.00000000e+00
+0.00000000e+00 2.83950617e-01 6.54320988e-01 4.50617284e-01 2.77777778e-02 0.00000000e+00 0.00000000e+00
diff --git a/regions/shapes/tests/reference/test_to_mask_elli-mode_subpixels-subpixels_5.txt b/regions/shapes/tests/reference/test_to_mask_elli-mode_subpixels-subpixels_5.txt
new file mode 100644
index 00000000..34f4f22d
--- /dev/null
+++ b/regions/shapes/tests/reference/test_to_mask_elli-mode_subpixels-subpixels_5.txt
@@ -0,0 +1,7 @@
+0.00000000e+00 0.00000000e+00 0.00000000e+00 2.00000000e-01 4.00000000e-01 8.00000000e-02 0.00000000e+00
+0.00000000e+00 0.00000000e+00 5.60000000e-01 1.00000000e+00 1.00000000e+00 8.40000000e-01 0.00000000e+00
+0.00000000e+00 4.40000000e-01 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00
+0.00000000e+00 9.20000000e-01 1.00000000e+00 1.00000000e+00 1.00000000e+00 9.20000000e-01 0.00000000e+00
+1.60000000e-01 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 5.60000000e-01 0.00000000e+00
+4.00000000e-02 9.60000000e-01 1.00000000e+00 1.00000000e+00 7.60000000e-01 4.00000000e-02 0.00000000e+00
+0.00000000e+00 3.20000000e-01 6.00000000e-01 4.40000000e-01 4.00000000e-02 0.00000000e+00 0.00000000e+00
diff --git a/regions/shapes/tests/reference/test_to_mask_rect-mode_center.txt b/regions/shapes/tests/reference/test_to_mask_rect-mode_center.txt
new file mode 100644
index 00000000..383eb69f
--- /dev/null
+++ b/regions/shapes/tests/reference/test_to_mask_rect-mode_center.txt
@@ -0,0 +1,8 @@
+0.00000000e+00 0.00000000e+00 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 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 0.00000000e+00
+1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.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 0.00000000e+00 1.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 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
diff --git a/regions/shapes/tests/reference/test_to_mask_rect-mode_subpixels-subpixels_1.txt b/regions/shapes/tests/reference/test_to_mask_rect-mode_subpixels-subpixels_1.txt
new file mode 100644
index 00000000..383eb69f
--- /dev/null
+++ b/regions/shapes/tests/reference/test_to_mask_rect-mode_subpixels-subpixels_1.txt
@@ -0,0 +1,8 @@
+0.00000000e+00 0.00000000e+00 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 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 0.00000000e+00
+1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.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 0.00000000e+00 1.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 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
diff --git a/regions/shapes/tests/reference/test_to_mask_rect-mode_subpixels-subpixels_18.txt b/regions/shapes/tests/reference/test_to_mask_rect-mode_subpixels-subpixels_18.txt
new file mode 100644
index 00000000..67e8bff7
--- /dev/null
+++ b/regions/shapes/tests/reference/test_to_mask_rect-mode_subpixels-subpixels_18.txt
@@ -0,0 +1,8 @@
+0.00000000e+00 0.00000000e+00 3.54938272e-01 6.48148148e-02 0.00000000e+00 0.00000000e+00 0.00000000e+00
+0.00000000e+00 2.46913580e-01 9.93827160e-01 8.98148148e-01 2.96296296e-01 0.00000000e+00 0.00000000e+00
+4.32098765e-02 8.64197531e-01 1.00000000e+00 1.00000000e+00 1.00000000e+00 6.32716049e-01 7.40740741e-02
+5.55555556e-01 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 6.29629630e-01
+2.62345679e-01 8.76543210e-01 1.00000000e+00 1.00000000e+00 1.00000000e+00 9.35185185e-01 1.01851852e-01
+0.00000000e+00 4.93827160e-02 5.86419753e-01 9.93827160e-01 1.00000000e+00 3.70370370e-01 0.00000000e+00
+0.00000000e+00 0.00000000e+00 0.00000000e+00 2.46913580e-01 5.98765432e-01 3.08641975e-03 0.00000000e+00
+0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
diff --git a/regions/shapes/tests/reference/test_to_mask_rect-mode_subpixels-subpixels_5.txt b/regions/shapes/tests/reference/test_to_mask_rect-mode_subpixels-subpixels_5.txt
new file mode 100644
index 00000000..532007b2
--- /dev/null
+++ b/regions/shapes/tests/reference/test_to_mask_rect-mode_subpixels-subpixels_5.txt
@@ -0,0 +1,8 @@
+0.00000000e+00 0.00000000e+00 3.60000000e-01 4.00000000e-02 0.00000000e+00 0.00000000e+00 0.00000000e+00
+0.00000000e+00 2.80000000e-01 1.00000000e+00 8.80000000e-01 2.80000000e-01 0.00000000e+00 0.00000000e+00
+4.00000000e-02 8.80000000e-01 1.00000000e+00 1.00000000e+00 1.00000000e+00 6.80000000e-01 8.00000000e-02
+5.20000000e-01 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00 6.00000000e-01
+2.80000000e-01 8.80000000e-01 1.00000000e+00 1.00000000e+00 1.00000000e+00 9.60000000e-01 1.20000000e-01
+0.00000000e+00 4.00000000e-02 6.00000000e-01 1.00000000e+00 1.00000000e+00 4.00000000e-01 0.00000000e+00
+0.00000000e+00 0.00000000e+00 0.00000000e+00 2.80000000e-01 6.00000000e-01 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
diff --git a/regions/shapes/tests/setup_package.py b/regions/shapes/tests/setup_package.py
new file mode 100644
index 00000000..d48b17f6
--- /dev/null
+++ b/regions/shapes/tests/setup_package.py
@@ -0,0 +1,2 @@
+def get_package_data():
+    return {'regions.shapes.tests': ['reference/*.txt', 'data/*.fits']}
diff --git a/regions/shapes/tests/test_api.py b/regions/shapes/tests/test_api.py
index 28245e27..ad2a4a1b 100644
--- a/regions/shapes/tests/test_api.py
+++ b/regions/shapes/tests/test_api.py
@@ -10,6 +10,8 @@
 from astropy import units as u
 from astropy.coordinates import ICRS
 from astropy.wcs import WCS
+
+from ...core.mask import Mask
 from ...core.core import Region, SkyRegion, PixelRegion
 from ...core.pixcoord import PixCoord
 from ..circle import CirclePixelRegion, CircleSkyRegion
@@ -39,7 +41,7 @@
 ]
 
 SKYPIX_MODES = ['local', 'affine', 'full']
-MASK_MODES = ['center', 'any', 'all', 'exact']
+MASK_MODES = ['center', 'exact', 'subpixels']
 COMMON_WCS = WCS(naxis=2)
 COMMON_WCS.wcs.ctype = 'RA---TAN', 'DEC--TAN'
 
@@ -92,8 +94,7 @@ def test_pix_to_shapely(region):
 def test_pix_to_mask(region, mode):
     try:
         mask = region.to_mask(mode=mode)
-        assert isinstance(mask, np.ndarray)
-        assert mask.ndim == 2
+        assert isinstance(mask, Mask)
     except NotImplementedError:
         pytest.xfail()
 
diff --git a/regions/shapes/tests/test_circle.py b/regions/shapes/tests/test_circle.py
index 1d24916d..a4fb4c02 100644
--- a/regions/shapes/tests/test_circle.py
+++ b/regions/shapes/tests/test_circle.py
@@ -1,6 +1,9 @@
 # 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
 from numpy.testing import assert_allclose
+
 from astropy import units as u
 from astropy.coordinates import SkyCoord
 from astropy.tests.helper import pytest, assert_quantity_allclose
@@ -9,6 +12,7 @@
 from astropy.wcs import WCS
 from ...core import PixCoord
 from ..circle import CirclePixelRegion, CircleSkyRegion
+from .utils import ASTROPY_LT_13
 
 try:
     import matplotlib
@@ -36,7 +40,10 @@ def test_circle_sky():
     center = SkyCoord(3 * u.deg, 4 * u.deg)
     reg = CircleSkyRegion(center, 2 * u.arcsec)
 
-    assert str(reg) == 'CircleSkyRegion\ncenter: <SkyCoord (ICRS): (ra, dec) in deg\n    (3.0, 4.0)>\nradius: 2.0 arcsec'
+    if ASTROPY_LT_13:
+        assert str(reg) == 'CircleSkyRegion\ncenter: <SkyCoord (ICRS): (ra, dec) in deg\n    (3.0, 4.0)>\nradius: 2.0 arcsec'
+    else:
+        assert str(reg) == 'CircleSkyRegion\ncenter: <SkyCoord (ICRS): (ra, dec) in deg\n    ( 3.,  4.)>\nradius: 2.0 arcsec'
 
 
 def test_transformation():
@@ -82,3 +89,11 @@ def test_plot():
     assert p.get_facecolor() == (1, 0, 0, 0.6)
 
     c.plot(ax, alpha=0.6)
+
+
+
+def test_to_mask():
+    center = PixCoord(3, 4)
+    reg = CirclePixelRegion(center, 2)
+    mask = reg.to_mask(mode='exact')
+    assert_allclose(np.sum(mask.data), 4 * np.pi)
diff --git a/regions/shapes/tests/test_ellipse.py b/regions/shapes/tests/test_ellipse.py
index 9b4fcaaf..4c6a7ee8 100644
--- a/regions/shapes/tests/test_ellipse.py
+++ b/regions/shapes/tests/test_ellipse.py
@@ -4,6 +4,7 @@
 from astropy.coordinates import SkyCoord
 from ...core import PixCoord
 from ..ellipse import EllipsePixelRegion, EllipseSkyRegion
+from .utils import ASTROPY_LT_13
 
 
 def test_ellipse_pixel():
@@ -17,6 +18,11 @@ def test_ellipse_sky():
     center = SkyCoord(3, 4, unit='deg')
     reg = EllipseSkyRegion(center, 3 * u.deg, 4 * u.deg, 5 * u.deg)
 
-    expected = ('EllipseSkyRegion\ncenter: <SkyCoord (ICRS): (ra, dec) in deg\n'
-                '    (3.0, 4.0)>\nminor: 3.0 deg\nmajor: 4.0 deg\nangle: 5.0 deg')
+    if ASTROPY_LT_13:
+        expected = ('EllipseSkyRegion\ncenter: <SkyCoord (ICRS): (ra, dec) in deg\n'
+                    '    (3.0, 4.0)>\nminor: 3.0 deg\nmajor: 4.0 deg\nangle: 5.0 deg')
+    else:
+        expected = ('EllipseSkyRegion\ncenter: <SkyCoord (ICRS): (ra, dec) in deg\n'
+                    '    ( 3.,  4.)>\nminor: 3.0 deg\nmajor: 4.0 deg\nangle: 5.0 deg')
+
     assert str(reg) == expected
diff --git a/regions/shapes/tests/test_masks.py b/regions/shapes/tests/test_masks.py
new file mode 100644
index 00000000..37456e4c
--- /dev/null
+++ b/regions/shapes/tests/test_masks.py
@@ -0,0 +1,43 @@
+# This file sets up detailed tests for computing masks with reference images
+
+import itertools
+import pytest
+
+from regions.core import PixCoord
+from regions.shapes.circle import CirclePixelRegion
+from regions.shapes.ellipse import EllipsePixelRegion
+from regions.shapes.rectangle import RectanglePixelRegion
+from astropy.io import fits
+from astropy import units as u
+
+
+REGIONS = [CirclePixelRegion(PixCoord(3.981987, 4.131378), 3.3411),
+           EllipsePixelRegion(PixCoord(3.981987, 4.131378), 3.3411, 2.2233, angle=32 * u.deg),
+           RectanglePixelRegion(PixCoord(3.981987, 4.131378), 4.3411, 5.2233, angle=32 * u.deg)]
+
+MODES = [{'mode': 'center'},
+         {'mode': 'exact'},
+         {'mode': 'subpixels', 'subpixels': 1},
+         {'mode': 'subpixels', 'subpixels': 5},
+         {'mode': 'subpixels', 'subpixels': 18}]
+
+
+def label(value):
+    if isinstance(value, CirclePixelRegion):
+        return 'circ'
+    elif isinstance(value, EllipsePixelRegion):
+        return 'elli'
+    elif isinstance(value, RectanglePixelRegion):
+        return 'rect'
+    else:
+        return '-'.join('{0}_{1}'.format(key, value) for key, value in sorted(value.items()))
+
+
+@pytest.mark.array_compare(fmt='text', write_kwargs={'fmt': '%12.8e'})
+@pytest.mark.parametrize(('region', 'mode'), itertools.product(REGIONS, MODES), ids=label)
+def test_to_mask(region, mode):
+    try:
+        mask = region.to_mask(**mode)
+    except NotImplementedError:
+        pytest.xfail()
+    return mask.data.astype(float)
diff --git a/regions/shapes/tests/test_point.py b/regions/shapes/tests/test_point.py
index 53db943d..6c570c75 100644
--- a/regions/shapes/tests/test_point.py
+++ b/regions/shapes/tests/test_point.py
@@ -3,6 +3,7 @@
 from astropy.coordinates import SkyCoord
 from ...core import PixCoord
 from ..point import PointPixelRegion, PointSkyRegion
+from .utils import ASTROPY_LT_13
 
 
 def test_ellipse_pixel():
@@ -16,4 +17,7 @@ def test_ellipse_sky():
     center = SkyCoord(3, 4, unit='deg')
     reg = PointSkyRegion(center)
 
-    assert str(reg) == 'PointSkyRegion\ncenter: <SkyCoord (ICRS): (ra, dec) in deg\n    (3.0, 4.0)>'
+    if ASTROPY_LT_13:
+        assert str(reg) == 'PointSkyRegion\ncenter: <SkyCoord (ICRS): (ra, dec) in deg\n    (3.0, 4.0)>'
+    else:
+        assert str(reg) == 'PointSkyRegion\ncenter: <SkyCoord (ICRS): (ra, dec) in deg\n    ( 3.,  4.)>'
diff --git a/regions/shapes/tests/test_polygon.py b/regions/shapes/tests/test_polygon.py
index 9520c34e..7f97b2e5 100644
--- a/regions/shapes/tests/test_polygon.py
+++ b/regions/shapes/tests/test_polygon.py
@@ -6,6 +6,7 @@
 from astropy.coordinates import SkyCoord
 from ...core import PixCoord
 from ..polygon import PolygonPixelRegion, PolygonSkyRegion
+from .utils import ASTROPY_LT_13
 
 
 @pytest.fixture
@@ -26,8 +27,12 @@ def test_polygon_pixel(pix_poly):
 
 
 def test_polygon_sky(sky_poly):
-    expected = ('PolygonSkyRegion\nvertices: <SkyCoord (ICRS): (ra, dec) in deg\n'
-                '    [(3.0, 3.0), (4.0, 4.0), (3.0, 4.0)]>')
+    if ASTROPY_LT_13:
+        expected = ('PolygonSkyRegion\nvertices: <SkyCoord (ICRS): (ra, dec) in deg\n'
+                    '    [(3.0, 3.0), (4.0, 4.0), (3.0, 4.0)]>')
+    else:
+        expected = ('PolygonSkyRegion\nvertices: <SkyCoord (ICRS): (ra, dec) in deg\n'
+                    '    [( 3.,  3.), ( 4.,  4.), ( 3.,  4.)]>')
     assert str(sky_poly) == expected
 
 
diff --git a/regions/shapes/tests/test_rectangle.py b/regions/shapes/tests/test_rectangle.py
index f23fb8cb..d334fc6c 100644
--- a/regions/shapes/tests/test_rectangle.py
+++ b/regions/shapes/tests/test_rectangle.py
@@ -4,6 +4,7 @@
 from astropy.coordinates import SkyCoord
 from ...core import PixCoord
 from ..rectangle import RectanglePixelRegion, RectangleSkyRegion
+from .utils import ASTROPY_LT_13
 
 
 def test_ellipse_pixel():
@@ -16,7 +17,10 @@ def test_ellipse_pixel():
 def test_ellipse_sky():
     center = SkyCoord(3, 4, unit='deg')
     reg = RectangleSkyRegion(center, 3 * u.deg, 4 * u.deg, 5 * u.deg)
-
-    expected = ('RectangleSkyRegion\ncenter: <SkyCoord (ICRS): (ra, dec) in deg\n'
-                '    (3.0, 4.0)>\nheight: 3.0 deg\nwidth: 4.0 deg\nangle: 5.0 deg')
+    if ASTROPY_LT_13:
+        expected = ('RectangleSkyRegion\ncenter: <SkyCoord (ICRS): (ra, dec) in deg\n'
+                    '    (3.0, 4.0)>\nheight: 3.0 deg\nwidth: 4.0 deg\nangle: 5.0 deg')
+    else:
+        expected = ('RectangleSkyRegion\ncenter: <SkyCoord (ICRS): (ra, dec) in deg\n'
+                    '    ( 3.,  4.)>\nheight: 3.0 deg\nwidth: 4.0 deg\nangle: 5.0 deg')
     assert str(reg) == expected
diff --git a/regions/shapes/tests/utils.py b/regions/shapes/tests/utils.py
new file mode 100644
index 00000000..40f66172
--- /dev/null
+++ b/regions/shapes/tests/utils.py
@@ -0,0 +1,5 @@
+from distutils.version import LooseVersion
+
+import astropy
+
+ASTROPY_LT_13 = LooseVersion(astropy.__version__) < LooseVersion('1.3')
diff --git a/regions/utils/__init__.py b/regions/utils/__init__.py
index 9c5e6f84..d516d594 100644
--- a/regions/utils/__init__.py
+++ b/regions/utils/__init__.py
@@ -1,3 +1,5 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 """Region utilities.
 """
+
+from .wcs_helpers import *
diff --git a/setup.cfg b/setup.cfg
index 8ed1dfe1..9426778c 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -10,6 +10,7 @@ show-response = 1
 [pytest]
 minversion = 2.2
 norecursedirs = build docs/_build
+addopts = --arraydiff
 # TODO: re-activate doctests.
 # Docs need to be marked up with skip directives to pass.
 # doctest_plus = enabled
@@ -41,4 +42,4 @@ github_project = astropy/regions
 # E901 - SyntaxError or IndentationError
 # E902 - IOError
 select = E101,W191,W291,W292,W293,W391,E111,E112,E113,E901,E902
-exclude = sphinx
\ No newline at end of file
+exclude = sphinx