Account for finite pixel size in PointPixelRegion.contains

[adonath]

Currently the PointPixelRegion.contains method always returns False. However this is not correct for the cases of sampled image data, because the size of the pixel is finite, namely 1 x 1 pix. This mean the .contains() method should return true if a position is contained in the finite support of the pixel.

[keflavich]

I'm not sure I agree - isn't PointPixelRegion definitionally a point, therefore it can't contain anything? The Pixel part just refers to the coordinate system.

I can certainly imagine situations where it would be useful to have a pixel-sized region, though. I think it would be best to have that be a RectanglePixelRegion with dimensions and orientation equal to a pixel size.

[adonath]

Thanks @keflavich! I guess the main question is what should happen in the following cases:

from regions import PointPixelRegion, PointSkyRegion
from astropy.coordinates import SkyCoord
from astropy.wcs import WCS
from astropy.io import fits
from astropy.utils.data import get_pkg_data_filename


center = PixCoord(0, 0)
point = PointPixelRegion(center=center)
print(point.contains(center))


fn = get_pkg_data_filename("galactic_center/gc_msx_e.fits", package='astropy.wcs.tests')
header = fits.getheader(fn)
wcs = WCS(header)

center = SkyCoord("0d", "0d", frame="galactic")
point = PointSkyRegion(center=center)
print(point.contains(center, wcs=wcs))

My expectation was, that both cases would return True, so the exact center should be contained in the region. I expected the behavior of a delta function, where, when evaluated at the exact position, the value is non-zero. When a wcs is given I think it can only be within the numerical precision and pixel size of the WCS.

This also allows to actually implement the PointPixelRegion.to_mask() method, which is currently undefined.

[larrybradley]

I agree with @keflavich here that a Point shouldn't contain anything. In your use case, if you want to check the center, you could use point.center == mycenter.

[adonath]

Thanks for the feedback @keflavich and @larrybradley, however I don't agree with you both here. Let me try to illustrate with a few use-cases, why I think the behavior should be the way I requested in my initial message. I think the main argument is consistent behavior with other regions and behavior in context of integration:

Creating Masks

from matplotlib import pyplot as plt
from regions import CirclePixelRegion, PointPixelRegion, PixCoord

circle = CirclePixelRegion(PixCoord(7.3, 6.8), radius=5)
mask = circle.to_mask(mode="exact").to_image((15, 15))
plt.imshow(mask);

Gives:

Now for the PointPixelRegion

import numpy as np
from matplotlib import pyplot as plt
from regions import PointPixelRegion, PixCoord

point = PointPixelRegion(PixCoord(7.3, 6.8))

# This currently raises NotImplementedError
mask = circle.to_mask().to_image((15, 15))

# However it is actually well defined in pixel space
y, x = np.indices((15, 15))
dx = np.abs(x - point.center.x)
dx = np.where(dx < 1, 1 - dx, 0)

dy = np.abs(y - point.center.y)
dy = np.where(dy < 1, 1 - dy, 0)

mask = dx * dy

plt.imshow(mask)

Should give:

The computation of the mask is numerically very well defined. There is no reason to not implement it.

Checking Overlap of Region with Mask Image

from matplotlib import pyplot as plt
from regions import CirclePixelRegion, PointPixelRegion, PixCoord

circle_mask = CirclePixelRegion(PixCoord(51, 51), radius=20)
mask = circle_mask.to_mask().to_image((101, 101)).astype(bool)

y, x = np.indices((101, 101))
coords_mask = PixCoord(x[mask], y[mask]) 

circle = CirclePixelRegion(PixCoord(37, 38), radius=10)
print(circle.contains(coords_mask).any())

point = PointPixelRegion(PixCoord(37, 38))
print(point.contains(coords_mask).any())

ax = plt.subplot()
ax.imshow(mask)
circle.plot(ax=ax, color="r")
point.plot(ax=ax, color="b")

Gives:

In this case the statement for the red circle clearly gives the right result, however the point region will always return False, even if it overlaps with the region defined by the mask image. This does not make sense to me.

Integration

Often regions are used to define integration boundaries. This is related to the mask behavior I mentioned above. See e.g.:

from matplotlib import pyplot as plt
from regions import CirclePixelRegion, PointPixelRegion, PixCoord

data = Gaussian2DKernel(x_stddev=5).array

circle = CirclePixelRegion(PixCoord(20, 20), radius=5)
weights = circle.to_mask(mode='exact').to_image(data.shape)
print((data * weights).sum())

point = PointPixelRegion(PixCoord(20, 20))

# Fails, but should give the value at the center of the Gauss,
# thats how a delta function behave when integrating
weights = point.to_mask(mode='exact').to_image(data.shape)
print((data * weights).sum())

ax = plt.subplot()
ax.imshow(data)
circle.plot(ax=ax, color="r")
point.plot(ax=ax, color="b")

And the illustration:

These are all use-cases where the behavior for the PointPixelRegion is very well defined, but current behavior in regions is either incorrect or not implemented / defined.

I'm not sure I agree - isn't PointPixelRegion definitionally a point, therefore it can't contain anything? The Pixel part just refers to the coordinate system.

I think this is only half-true. To me the pixel coordinate system implies, that the data space is discretized / represented as a pixelized image. This is also the behavior currently implemented for other regions, especially with the .to_mask() method or the definition of discrete bounding boxes. I can see your point for the .contains() method, but the behavior of e.g. point_region.to_mask(), where one clearly "enters" discretized image space, is uncontroversial to me.

[keflavich]

I agree that you've shown cases where a region that contains a single pixel is useful. The question before was mostly one of terminology; in a mathematically purist view, a point can't contain anything. But, I think you're making a good case that a PointRegion is more useful if it can be defined to contain something, or at least, if it can be defined to be analogous to the Dirac delta (if used in our discrete space).

Is there any possible downside to treating PointPixelRegion.contains as a check of whether the test pixel is equal to PointPixelRegion.center? I don't see one immediately but I want to at least voice the thought.

I also see the use of defining a to_mask method that returns the fractional overlap between each pixel and the pixel region. However, we should make quite clear what the different modes do - exact should be as in your example, where partial overlap is included. Should center be "true iff the pixel center is the region center" or should it be "true if the region center is within the bounds of the pixel"?

I'd be in favor of a PR that implements this functionality. Maybe we'd call this just a PixelRegion instead of a PointPixelRegion, though?

I suggest that this is a critical case to document very explicitly. For most of the regions in this package, a small difference in interpretation of what an edge is will not result in major differences in the outcome of a project. For these pixel regions, when we're potentially talking about a boolean difference between 'included' and 'not included', the consequences could be more pronounced.

[adonath]

I guess the general issue is more conceptional:

• For which region methods / properties is the pixel space used "just as coordinate system"?
• For which region methods / properties are regions actually "digitized" or "discretized"?

Indeed .contains() just uses the pixel coordinate system and does not discretize for extended regions. In this case I would agree that one should maybe just check for the equivalency of the position and the point region center and not handle it as a 1x1 pixel region. However for all methods where the region is actually discretized e.g. .to_mask() it should be handled as a 1x1 pixel region. So my proposal would be:

• Adapt PointPixelRegion.contains() to check for the equivalency of the position to the region center
• Actually implement PointPixelRegion.to_mask() and PointPixelRegion.bounding_box with the behavior I proposed above of handling it as a 1x1 pixel region

One more thing I found: I think it's non ideal to return a Line2D patch in PixelSkyRegion instead of a normal matplotlib Patch. Because this will always require an extra check when looping over or plotting multiple regions. Because the Line2D plotting arguments are different than the ones for Patch. E.g. facecolor and edgecolor will work for Patch, but the corresponding argument for Line2D are markeredgecolor and markerfacecolor. So one typically ends up with code like:

for region_pix in regions_pix:
    
    if isinstance(region, PointSkyRegion):
        artist = region_pix.as_artist(**kwargs_line2d)
    else:
        artist = region_pix.as_artist(**kwargs)

    ax.add_artist(artist)

Can this behavior be unified as well? Or is it intended?

[keflavich]

I agree, I like your proposed points.

Regarding the inconsistency of the Line2D object: I agree, it would be better to have these be consistent. Is there a matplotlib object type that fits this? Maybe PathPatch? (could you re-raise this in a different Issue, though?)