-
-
Notifications
You must be signed in to change notification settings - Fork 54
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Added initial implementation of Polygon to_mask and as_patch #87
Changes from 4 commits
91c56b7
12df57b
1063e9f
8b232cf
bf3f349
d90426d
9f542c1
b085a2a
c6cef1e
4f7a97f
e02710d
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,5 @@ | ||
cimport numpy as np | ||
|
||
ctypedef np.float64_t DTYPE_t | ||
|
||
cdef int point_in_polygon(double x, double y, np.ndarray[DTYPE_t, ndim=1] vx, np.ndarray[DTYPE_t, ndim=1] vy) |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,112 @@ | ||
# The code in this file was adapted from code written by Greg von Winckel: | ||
# | ||
# https://github.com/gregvw/pnpoly | ||
# | ||
# and which was released under the following license: | ||
# | ||
# ---------------------------------------------------------------------------- | ||
# | ||
# The MIT License (MIT) | ||
# | ||
# Copyright (c) 2014 Greg von Winckel | ||
# | ||
# Permission is hereby granted, free of charge, to any person obtaining a copy of | ||
# this software and associated documentation files (the "Software"), to deal in | ||
# the Software without restriction, including without limitation the rights to | ||
# use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of | ||
# the Software, and to permit persons to whom the Software is furnished to do so, | ||
# subject to the following conditions: | ||
# | ||
# The above copyright notice and this permission notice shall be included in all | ||
# copies or substantial portions of the Software. | ||
# | ||
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR | ||
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS | ||
# FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR | ||
# COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER | ||
# IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN | ||
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. | ||
# | ||
# ---------------------------------------------------------------------------- | ||
# | ||
# This code was itself adapted from code written by W. Randolph Franklin: | ||
# | ||
# http://www.ecse.rpi.edu/~wrf/Research/Short_Notes/pnpoly.html | ||
# | ||
# and released under the following license: | ||
# | ||
# Copyright (c) 1970-2003, Wm. Randolph Franklin | ||
# | ||
# Permission is hereby granted, free of charge, to any person obtaining a copy | ||
# of this software and associated documentation files (the "Software"), to deal | ||
# in the Software without restriction, including without limitation the rights | ||
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell | ||
# copies of the Software, and to permit persons to whom the Software is | ||
# furnished to do so, subject to the following conditions: | ||
# | ||
# Redistributions of source code must retain the above copyright notice, this | ||
# list of conditions and the following disclaimers. | ||
# | ||
# Redistributions in binary form must reproduce the above copyright notice in | ||
# the documentation and/or other materials provided with the distribution. | ||
# | ||
# The name of W. Randolph Franklin may not be used to endorse or promote | ||
# products derived from this Software without specific prior written permission. | ||
# | ||
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR | ||
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, | ||
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE | ||
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER | ||
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, | ||
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE | ||
# SOFTWARE. | ||
|
||
cimport cython | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I think this line can be removed!?
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Done |
||
import numpy as np | ||
cimport numpy as np | ||
|
||
DTYPE_BOOL = np.bool | ||
ctypedef np.uint8_t DTYPE_BOOL_t | ||
|
||
|
||
def points_in_polygon(np.ndarray[DTYPE_t, ndim=1] x, | ||
np.ndarray[DTYPE_t, ndim=1] y, | ||
np.ndarray[DTYPE_t, ndim=1] vx, | ||
np.ndarray[DTYPE_t, ndim=1] vy): | ||
|
||
cdef int i, n | ||
cdef np.ndarray[np.uint8_t, ndim=1, cast=True] result | ||
|
||
n = x.shape[0] | ||
|
||
result = np.zeros(n, DTYPE_BOOL) | ||
|
||
for i in range(n): | ||
result[i] = point_in_polygon(x[i], y[i], vx, vy) | ||
|
||
return result | ||
|
||
|
||
cdef int point_in_polygon(double x, double y, | ||
np.ndarray[DTYPE_t, ndim=1] vx, | ||
np.ndarray[DTYPE_t, ndim=1] vy): | ||
""" | ||
Determine whether a test point (x, y) is within a polygon defined by a set | ||
of vertices (vx, vy) | ||
""" | ||
|
||
cdef int i, j, k, m, n | ||
cdef int result | ||
|
||
n = vx.shape[0] | ||
|
||
result = 0 | ||
|
||
for i in range(n): | ||
j = (i+n-1) % n | ||
if(((vy[i] > y) != (vy[j] > y)) and \ | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. This if statement spanning over 4 lines is very hard to read. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. The 'and' short-circuits the evaluation of the second part of the if statement, so I don't want to pre-compute the second part - but I'll clean it up. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. PyCharm suggests to remove the backslashes, and I agree, I think generally just parens to group multi-line expressions are the preferred style over backslashes. Not sure if that's PEP8 or not. |
||
(x < (vx[j]-vx[i]) * (y-vy[i])/ \ | ||
(vy[j]-vy[i]) + vx[i]) ): | ||
result += 1 | ||
|
||
return result % 2 |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,117 @@ | ||
# Licensed under a 3-clause BSD style license - see LICENSE.rst | ||
|
||
from __future__ import (absolute_import, division, print_function, | ||
unicode_literals) | ||
|
||
import numpy as np | ||
cimport numpy as np | ||
|
||
|
||
__all__ = ['polygonal_overlap_grid'] | ||
|
||
|
||
DTYPE = np.float64 | ||
ctypedef np.float64_t DTYPE_t | ||
|
||
from .pnpoly cimport point_in_polygon | ||
|
||
__all__ = ['polygonal_overlap_grid'] | ||
|
||
|
||
def polygonal_overlap_grid(double xmin, double xmax, double ymin, double ymax, | ||
int nx, int ny, | ||
np.ndarray[DTYPE_t, ndim=1] vx, | ||
np.ndarray[DTYPE_t, ndim=1] vy, | ||
int use_exact, | ||
int subpixels): | ||
""" | ||
polygonal_overlap_grid(xmin, xmax, ymin, ymax, nx, ny, r, | ||
use_exact, subpixels) | ||
|
||
Area of overlap between a polygon and a pixel grid. | ||
|
||
Parameters | ||
---------- | ||
xmin, xmax, ymin, ymax : float | ||
Extent of the grid in the x and y direction. | ||
nx, ny : int | ||
Grid dimensions. | ||
vx, vy : `numpy.ndarray` | ||
The vertices of the polygon | ||
use_exact : 0 or 1 | ||
If ``1`` calculates exact overlap, if ``0`` uses ``subpixel`` number | ||
of subpixels to calculate the overlap. | ||
subpixels : int | ||
Each pixel resampled by this factor in each dimension, thus each | ||
pixel is divided into ``subpixels ** 2`` subpixels. | ||
|
||
Returns | ||
------- | ||
frac : `~numpy.ndarray` (float) | ||
2-d array of shape (ny, nx) giving the fraction of the overlap. | ||
""" | ||
|
||
cdef unsigned int i, j | ||
cdef double x, y, dx, dy, d, pixel_radius | ||
cdef double bxmin, bxmax, bymin, bymax | ||
cdef double pxmin, pxcen, pxmax, pymin, pycen, pymax | ||
|
||
# Define output array | ||
cdef np.ndarray[DTYPE_t, ndim=2] frac = np.zeros([ny, nx], dtype=DTYPE) | ||
|
||
if use_exact == 1: | ||
raise NotImplementedError("Exact mode has not been implemented for " | ||
"polygonal apertures") | ||
|
||
# Find the width of each element in x and y | ||
dx = (xmax - xmin) / nx | ||
dy = (ymax - ymin) / ny | ||
|
||
# Define bounding box | ||
bxmin = vx.min() | ||
bxmax = vx.max() | ||
bymin = vy.min() | ||
bymax = vy.max() | ||
|
||
for i in range(nx): | ||
pxmin = xmin + i * dx # lower end of pixel | ||
pxmax = pxmin + dx # upper end of pixel | ||
if pxmax > bxmin and pxmin < bxmax: | ||
for j in range(ny): | ||
pymin = ymin + j * dy | ||
pymax = pymin + dy | ||
if pymax > bymin and pymin < bymax: | ||
frac[j, i] = polygonal_overlap_single_subpixel(pxmin, pymin, | ||
pxmax, pymax, | ||
vx, vy, subpixels) | ||
|
||
return frac | ||
|
||
|
||
cdef double polygonal_overlap_single_subpixel(double x0, double y0, | ||
double x1, double y1, | ||
np.ndarray[DTYPE_t, ndim=1] vx, | ||
np.ndarray[DTYPE_t, ndim=1] vy, | ||
int subpixels): | ||
""" | ||
Return the fraction of overlap between a polygon and a single pixel | ||
with given extent, using a sub-pixel sampling method. | ||
""" | ||
|
||
cdef unsigned int i, j | ||
cdef double x, y, dx, dy | ||
cdef double frac = 0. # Accumulator. | ||
|
||
dx = (x1 - x0) / subpixels | ||
dy = (y1 - y0) / subpixels | ||
|
||
x = x0 - 0.5 * dx | ||
for i in range(subpixels): | ||
x += dx | ||
y = y0 - 0.5 * dy | ||
for j in range(subpixels): | ||
y += dy | ||
if point_in_polygon(x, y, vx, vy) == 1: | ||
frac += 1. | ||
|
||
return frac / (subpixels * subpixels) |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -1,6 +1,11 @@ | ||
# Licensed under a 3-clause BSD style license - see LICENSE.rst | ||
from __future__ import absolute_import, division, print_function, unicode_literals | ||
from ..core import PixelRegion, SkyRegion | ||
|
||
import numpy as np | ||
|
||
from ..core import PixelRegion, SkyRegion, Mask, BoundingBox | ||
from .._geometry import polygonal_overlap_grid | ||
from .._geometry.pnpoly import points_in_polygon | ||
|
||
__all__ = ['PolygonPixelRegion', 'PolygonSkyRegion'] | ||
|
||
|
@@ -35,8 +40,10 @@ def area(self): | |
raise NotImplementedError | ||
|
||
def contains(self, pixcoord): | ||
# TODO: needs to be implemented | ||
raise NotImplementedError | ||
return points_in_polygon(pixcoord.x.astype(float), | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. There's this error:
Besides fixing the error, may I suggest you change the code to do the type conversions on separate lines before the There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Should be fixed |
||
pixcoord.y.astype(float), | ||
self.vertices.x.astype(float), | ||
self.vertices.y.astype(float)).astype(bool) | ||
|
||
def to_shapely(self): | ||
# TODO: needs to be implemented | ||
|
@@ -48,16 +55,50 @@ def to_sky(self, wcs, mode='local', tolerance=None): | |
|
||
@property | ||
def bounding_box(self, mode='center'): | ||
# TODO: needs to be implemented | ||
raise NotImplementedError | ||
xmin = self.vertices.x.min() | ||
xmax = self.vertices.x.max() | ||
ymin = self.vertices.y.min() | ||
ymax = self.vertices.y.max() | ||
return BoundingBox._from_float(xmin, xmax, ymin, ymax) | ||
|
||
def to_mask(self, mode='center'): | ||
# TODO: needs to be implemented | ||
raise NotImplementedError | ||
def to_mask(self, mode='center', subpixels=5): | ||
|
||
self._validate_mode(mode, subpixels) | ||
|
||
if mode == 'center': | ||
mode = 'subpixels' | ||
subpixels = 1 | ||
|
||
# Find bounding box and mask size | ||
bbox = self.bounding_box | ||
ny, nx = bbox.shape | ||
|
||
# Find position of pixel edges and recenter so that circle is at origin | ||
xmin = float(bbox.ixmin) - 0.5 | ||
xmax = float(bbox.ixmax) - 0.5 | ||
ymin = float(bbox.iymin) - 0.5 | ||
ymax = float(bbox.iymax) - 0.5 | ||
|
||
if mode == 'subpixels': | ||
use_exact = 0 | ||
else: | ||
use_exact = 1 | ||
|
||
fraction = polygonal_overlap_grid( | ||
xmin, xmax, ymin, ymax, nx, ny, | ||
self.vertices.x.astype(float), | ||
self.vertices.y.astype(float), | ||
use_exact, subpixels) | ||
|
||
return Mask(fraction, bbox=bbox) | ||
|
||
def as_patch(self, **kwargs): | ||
# TODO: needs to be implemented | ||
raise NotImplementedError | ||
""" | ||
Matplotlib patch object for this region (`matplotlib.patches.Polygon`). | ||
""" | ||
from matplotlib.patches import Polygon | ||
xy = np.vstack([self.vertices.x, self.vertices.y]).transpose() | ||
return Polygon(xy=xy, **kwargs) | ||
|
||
|
||
class PolygonSkyRegion(SkyRegion): | ||
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Add standard license line at the top, even if it's a very small file.
Can you add a docstring with a sentence saying why we need a
pxd
file at all?(my Cython is always rusty because I use it so rarely. It's the same for most Python programmers.)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
My knowledge of Cython is 'because Cython complains if it's not there' ;)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Done (it's because it's needed for cimport to work)