Support computing area for compound regions

[adonath]

The following example currently raises a NotImplementedError:

from regions import CirclePixelRegion, PixCoord

circle_1 = CirclePixelRegion(
    center=PixCoord(0, 0),
    radius=10
)

circle_2 = CirclePixelRegion(
    center=PixCoord(10, 0),
    radius=10
)

compound = circle_1 & circle_2
print(compound.area)

It would be great to have this supported!

[larrybradley]

This would be a very nice feature! However, it would be very challenging to implement because the area attribute represents the exact analytical area of the shape. To calculate this generally for compound shape combinations is not trivial. Even for your very simple example above with two circles, the overlapping area is (https://mathworld.wolfram.com/Circle-CircleIntersection.html): 👀

If you are looking for something based on the region pixel masks (as an approximation or as the real "pixel" area), you can use the RegionMask object. It requires using the mode='center' method to give "full" (not fractional) pixel weights:

from regions import CirclePixelRegion, PixCoord

circle1 = CirclePixelRegion(PixCoord(15, 20), radius=10)
circle2 = CirclePixelRegion(PixCoord(25, 20), radius=10)
compound = circle1 & circle2
mask = compound.to_mask(mode='center')
np.count_nonzero(mask.data)
# 117

The exact overlap area for this "simple" case is:

def circle_overlap_area(radius1, radius2, distance):
    f1 = radius1**2 * np.arccos((distance**2 + radius1**2 - radius2**2) / (2 * distance * radius1))
    f2 = radius2**2 * np.arccos((distance**2 + radius2**2 - radius1**2) / (2 * distance * radius2))
    a = radius1 + radius2 - distance
    b = distance + radius1 - radius2
    c = distance - radius1 + radius2
    d = distance + radius1 + radius2
    f3 = 0.5 * np.sqrt(a * b * c * d)
    return f1 + f2 - f3

circle_overlap_area(10, 10, 10)
# 122.83696986087567

maskimg = mask.to_image((40, 40))
plt.imshow(maskimg)
circle1.plot(color='blue', lw=3)
circle2.plot(color='green', lw=3)

[adonath]

Thanks @larrybradley for this nice illustration! I agree that even the case of two circle is far from simple. Especially because compound regions are not only intersections, so the different cases of operators need to be handled as well.

However I think even supporting a few standard cases will still be useful. I guess the complexity arises, once regions are overlapping. So even raising an error only in this case could be fine. BTW: I think a method like PixelRegion.overlaps(other_region=) could be useful as well (I might open a separate issue for this).

The numerical approach seems like a good solution as well. With a little workaround the exact case can be supported as well (not sure how this integrates in regions though....):

import numpy as np
from regions import CirclePixelRegion, PixCoord

circle1 = CirclePixelRegion(PixCoord(15, 20), radius=10)
circle2 = CirclePixelRegion(PixCoord(25, 20), radius=10)
mask_1 = circle1.to_mask(mode='exact')
mask_2 = circle2.to_mask(mode='exact')

shape = (50, 50)

m_1 = mask_1.to_image(shape=shape)
m_2 = mask_2.to_image(shape=shape)
m_all = m_1 * m_2

print(m_all.sum())
#122.80909547049568

[larrybradley]

The "exact" mode is a lot better here, but still an approximation because of the combination of the fractional-pixel masks (also see #73).