Skip to content
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

Merged
merged 11 commits into from
Nov 30, 2016

Conversation

astrofrog
Copy link
Member

@astrofrog astrofrog commented Nov 30, 2016

This implements to_mask for Polygon in line with other shapes.

It does not use code from ginga for example - I did look at that code, but the issue is we need to be able to efficiently compute the sub-pixel overlap without making large arrays of pixels (the ginga code operates with vectorized operations on Numpy arrays).

However, the code I'm adding here for the point/polygon stuff is pretty lightweight. Example usage:

import numpy as np
import matplotlib.pyplot as plt
from regions import PolygonPixelRegion, PixCoord

vx = np.array([-3, 5, 2, 0.5, -1]).astype(float)
vy = np.array([-2., -3., 4., 2., 3.]).astype(float)

vertices = PixCoord(vx, vy)

region = PolygonPixelRegion(vertices)

mask = region.to_mask(mode='subpixels', subpixels=10)

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.imshow(mask, cmap=plt.cm.viridis, origin='lower',
          interpolation='nearest',
          extent=mask.bbox.extent)
ax.add_patch(region.as_patch(edgecolor='white', facecolor='none'))
fig.savefig('polygon.png')

polygon

This needs tests of course, but waiting first to see if people are open to this.

@cdeil cdeil added this to the 0.2 milestone Nov 30, 2016
@cdeil cdeil self-assigned this Nov 30, 2016
@cdeil
Copy link
Member

cdeil commented Nov 30, 2016

@astrofrog - Awesome!

I have to run now, but I'll try to review tonight.

What's the origin of the code in regions/_geometry/pnpoly.pyx? Which project? Do you have a URL where you got it?

# The code in this file was adapted from code written by Greg von Winckel and
# which was released under the following license:
#
# ----------------------------------------------------------------------------
#
# The MIT License (MIT)
#
# Copyright (c) 2014 Greg von Winckel

I think we should also acknowledge them in the docs or changelog.

The ideal situation would be if they agree to contribute the code under our BSD-3 license via a comment here or by email, just so that we don't complicate the license situation of Astropy.
Alternatively, having some MIT licensed code in Astropy is fine of course as well.

@astrofrog
Copy link
Member Author

astrofrog commented Nov 30, 2016

The ideal situation would be if they agree to contribute the code under our BSD-3 license via a comment here or by email, just so that we don't complicate the license situation of Astropy.
Alternatively, having some MIT licensed code in Astropy is fine of course as well.

Just to be clear, when we take MIT code and we put it in Astropy, we are re-licensing it as BSD but just acknowledging the original license. So the code is not currently MIT, we are releasing it in a BSD package, but we're just including the original license for reference.

@astrofrog
Copy link
Member Author

I added the URLs where the original code was adapted from

@astrofrog
Copy link
Member Author

Another note: I think there are some inefficiencies in the Cython code at the moment that could be sped up, but I'm not too concerned about performance right now. I think maybe once we have a lot of the functionality in place we could do a performance sprint where we try and speed everything up where needed?

@cdeil
Copy link
Member

cdeil commented Nov 30, 2016

Just to be clear, when we take MIT code and we put it in Astropy, we are re-licensing it as BSD but just acknowledging the original license. So the code is not currently MIT, we are releasing it in a BSD package, but we're just including the original license for reference.

Thanks. The origin / license is now clearer in regions/_geometry/pnpoly.pyx.

Still, it doesn't say what you write in that comment. If we're re-licensing as BSD, then I'd suggest the first line in the file to be our usual

# Licensed under a 3-clause BSD style license - see LICENSE.rst

and then after that to explain the origin and keep the MIT license.

In the end it's just ~ 15 lines of code. We could just do a BSD licensed clean-room implementation before putting it in Astropy core, i.e. one person writes down or links to a description of the algorithm, and a second writes the code from that only. I've always wanted to do that. But not now, I'm OK to merge this and will shut about license now and have a look at the code. :-)

@@ -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),
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There's this error:
https://travis-ci.org/astropy/regions/jobs/180103544#L819

E       AttributeError: 'int' object has no attribute 'astype'

Besides fixing the error, may I suggest you change the code to do the type conversions on separate lines before the points_in_polygon call and put them in temp variables?
The efficiency is the same, and the error message for invalid data like here will improve.
At the moment from the traceback it's not clear which variable / line is the problematic one.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should be fixed

@cdeil
Copy link
Member

cdeil commented Nov 30, 2016

In your first comment on this PR, you have a plot with a polygon and a mask. Very nice!
Please add it to the docs, either the polygon class docstring, or the high-level getting started docs section.

@astrofrog
Copy link
Member Author

@cdeil - the 'getting started' documentation is growing quite long - what do you think about splitting it up into multiple pages?

@@ -0,0 +1,5 @@
cimport numpy as np
Copy link
Member

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.)

Copy link
Member Author

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' ;)

Copy link
Member Author

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)

Copy link
Member

@cdeil cdeil left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've left two more small inline comments.

I didn't review the algorithm / math.

@astrofrog - Should we merge this when travis-ci passes, or do you want more review here?
Maybe we leave it open for a day or two, and if people want to contribute tests or docs, they can make pull requests against this branch?
Or we merge in master now, to increase the chance this gets some usage before the regions 0.2 release?

# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

cimport cython
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this line can be removed!?

cimport cython

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done


for i in range(n):
j = (i+n-1) % n
if(((vy[i] > y) != (vy[j] > y)) and \
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This if statement spanning over 4 lines is very hard to read.
Can you compute the parts on separate lines and store them in temp variables, and then have a one-line conditional?

Copy link
Member Author

@astrofrog astrofrog Nov 30, 2016

Choose a reason for hiding this comment

The 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.

Copy link
Member

Choose a reason for hiding this comment

The 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.

@astrofrog
Copy link
Member Author

@cdeil - I'll add an arraydiff test now so there is at least some form of testing. But we can merge before the docs. I think for 0.2 it'd be nice to split up the docs into more pages - at the moment 'getting started' is a bit monolithic - would you have time to work on that?

@cdeil
Copy link
Member

cdeil commented Nov 30, 2016

the 'getting started' documentation is growing quite long - what do you think about splitting it up into multiple pages?

I was thinking about this today, too.
Yes, splitting it onto multiple pages is fine with me.

The other idea I had was that we create an IPython notebook with some nice demos to get people's attention. We could make it a binder, so that they can try it out online.
The only problem with that is that we grow the repo size, probably by a few MB with images, and I don't think this is something we want to end up in the Astropy core repo at the end.

@astrofrog
Copy link
Member Author

astrofrog commented Nov 30, 2016

The only problem with that is that we grow the repo size, probably by a few MB with images, and I don't think this is something we want to end up in the Astropy core repo at the end.

What about putting that in the astropy-tutorials repo? It's already on binder that way.

@cdeil
Copy link
Member

cdeil commented Nov 30, 2016

would you have time to work on that?

Very little, but yes, I'll try to put some more time into docs before making the release.

@cdeil
Copy link
Member

cdeil commented Nov 30, 2016

What about putting that in the astropy-tutorials repo?

Yes, probably that would be the best solution.
Made an issue as a reminder and for further discussion on the notebook: #88

@cdeil
Copy link
Member

cdeil commented Nov 30, 2016

I've split out the discussion / work on high-level docs re-structuring into a separate issue:
#89

@astrofrog
Copy link
Member Author

@cdeil - I implemented your comments. If the tests pass, I think we can merge this and add documentation after.

@cdeil
Copy link
Member

cdeil commented Nov 30, 2016

For the polygon code, only the polygonal_overlap_single_subpixel docstring hints at the method that's implemented.

@astrofrog - If you know which algorithm is implemented in the polygon functions, could you mention it by name and link e.g. to wikipedia?

I mean you are linking to https://github.com/gregvw/pnpoly/blob/master/pnpoly.pyx and https://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html, but from looking at the code and those pages for 1 minute it's not clear to me which algorithm we have.

Is it one of the algorithms mentioned here?
https://en.wikipedia.org/wiki/Point_in_polygon

@cdeil
Copy link
Member

cdeil commented Nov 30, 2016

+1 to merge and add more tests / docs in the coming days.

@astrofrog
Copy link
Member Author

@cdeil - the algorithm is the even-odd rule

@astrofrog astrofrog merged commit f456f49 into astropy:master Nov 30, 2016
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

2 participants