Pull regions code from photutils, Ginga or LSST?

[cdeil]

At ADASS 2016 @ejeschke and I had a chat about regions and I wanted to leave some notes here.

There's a pretty complete set of regions implemented in Ginga (list). They have a contains_arr method with containment implemented using Numpy array expressions. See polygon as example.

Ginga only has pixel regions, no sky regions, and the regions code there contains GUI stuff. @ejeschke might be interested in re-using astropy.regions from Ginga, if it's done in an extensible way so that the Ginga GUI code can be added e.g. by holding astropy.region data members or by subclassing. Also @ejeschke isn't a fan of the double class hierarchy we have in regions and photutils where the data members are complex and non-uniform objects like SkyCoord or Angle. I think in Ginga he has something more uniform where every region can be expressed as a list of points. Maybe it's similar to how pyregion.Shape represents regions? @ejeschke, could you elaborate on this point and how the double class hierarchy could be avoided and then we'll try to do a pro-con of the options?

There's also photutils.aperture where so far the plan was to pull the code for ellipse, rectangle from into this regions package and then in the far future photutils would switch to astropy.regions. photutils is Cython-backed code and has sub-pixel stuff, which Ginga doesn't have. Do we want sub-pixel contains fractions in astropy.regions? It's extra work to develop / test for all regions.

@timj has written up (in Section 4.2) some info about he LSST regions code and thoughts on how it could be added or interfaced to astropy.regions here: https://www.lsst.org/scientists/publications/investigating-interoperability-lsst-data-management-software-stack-astropy

Personally I'm inclined to start by pulling the pixel polygon code from Ginga, because that's what we need for data analysis in Gammapy. But if we find "the right way" forward, I'll spend dev time to help implement whatever solution we decide to go with.

@astropy/regions-developers - Thoughts on what we should do? Anyone have time to work on this in the coming months?

[astrofrog]

Do we want sub-pixel contains fractions in astropy.regions? It's extra work to develop / test for all regions.

Yes, I think we do (as I had originally hinted at in the original API: https://github.com/astropy/regions/blob/master/regions/core/core.py#L149)

However, I don't think all regions need to implement all modes for the package to be useful (it's ok to return NotImplementedError). I think using the center of the pixel should be implemented for all regions, and then the exact mode could be added when available. The 'all'/'any' modes can be easily implemented by using 'exact' and checking if the answer is =1 (for all) or >0 (for any).

Note that Matplotlib also has a fast point-in-polygon test written in C that we should investigate and compare to the ginga version.

[larrybradley]

I also think we want to keep the subpixel mode, at least for photometry.

@astrofrog The line you refer to doesn't actually include a 'subpixels' mode (we should add it). Also, there would need to be an additional subpixels keyword to specify the number of subpixels for this mode (like what is currently in photutils).

[larrybradley]

I submitted a PR for adding a subpixels mode and keyword here: #70

[ejeschke]

I think in Ginga he has something more uniform where every region can be expressed as a list of points. Maybe it's similar to how pyregion.Shape represents regions? @ejeschke, could you elaborate on this point and how the double class hierarchy could be avoided and then we'll try to do a pro-con of the options?

First to the bit about being able to access every shape via an attribute that reveals an array of points/coordinates. Such a feature is extremely useful in code where you want to manipulate a shape "from the outside" without really understanding the details of the shape. An example is where the user is editing a point graphically, by moving a single point within a an object that is a point, line, or polygon. If you have a setter/getter (or properties) defined then the code for changing the point simply needs to set the index and value of the point as it is moved and the shape is updated--there is no need for special code to understand that in this kind of object I need to access/set x1, y1 in this object and x, y in this object, etc. It's just point[i] (or coord[i]). The same thing can be exploited for generic code for rotating an object, getting the bounding box, etc. As I understand the region code so far, this should be fairly easy to support simply with a property defined.

As for the double class hierarchy, I'm not really opposed so much as wondering whether the alternatives have been explored thoroughly. From the users API perspective, isn't it reasonable to say I want to make a Polygon and then indicate with a parameter (or the type of coordinates passed) that it is a pixel based polygon or a sky polygon? Currently, the hierarchy looks like Region => PixelRegion => PolygonPixelRegion, etc. and a separate path of Region => SkyRegion => PolygonSkyRegion. Why not Region => PathType => Polygon with a mixin class of SkySurface or PixelSurface? (paths and other objects share many features of polygons, thus the PathType intermediate class)

@ejeschke might be interested in re-using astropy.regions from Ginga, if it's done in an extensible way so that the Ginga GUI code can be added e.g. by holding astropy.region data members or by subclassing.

This will be true for anyone wanting to use astropy regions with e.g. matplotlib or their own GUI. What I want to avoid with ginga is having to maintain a parallel mapping of points/coords so that two sets of points do not need to be updated when manipulating a shape, either in the GUI or programatically for developers. To illustrate this point, think only about how you would maintain a matplotlib polygon and an astropy regions polygon and keep them in sync.

[ejeschke]

Oh, BTW, @cdeil, is the PixelCoordinate (sorry if I am not remembering the name correctly) that you described to me planned to be merged into astropy regions? It would be nice if Cartesian or pixel coordinates can be used with these shapes using a somewhat common way of manipulating them compared to other coordinates.

[TallJimbo]

LSST has some high-quality C++ code for difficult-to-implement spherical geometry operations ( https://github.com/lsst/sphgeom) that can be called from Python, but I don't think we have any particularly unique insight to offer on the APIs for generic spherical and pixel geometry objects.

I do think we have some API design insight (and wrapped C++ code) for regions that represent discrete pixel regions (i.e. regions defined to contain only complete pixels) using run-length-encoding and integer bounding boxes. It doesn't sound like that's on the table here yet, though.

As a side note, you might want to consider "Cartesian" and "Spherical" instead of "Pixel" and "Sky" in terms of naming, unless there are features of these unique to sky coordinate systems and pixel coordinate systems.

[pllim]

@ejeschke started something over on Ginga side, see ejeschke/ginga#787