Skip to content

Commit

Permalink
Move exclusion mask to maker / run method
Browse files Browse the repository at this point in the history
  • Loading branch information
maxnoe committed Jan 27, 2022
1 parent 4b57311 commit e4536a9
Show file tree
Hide file tree
Showing 3 changed files with 51 additions and 25 deletions.
46 changes: 30 additions & 16 deletions gammapy/makers/background/reflected.py
Expand Up @@ -22,7 +22,7 @@ class RegionsFinder(metaclass=ABCMeta):
'''Baseclass for regions finders'''

@abstractmethod
def run(self, region, center):
def run(self, region, center, exclusion_mask=None):
"""Find regions to calculate background counts.
Parameters
Expand All @@ -31,6 +31,8 @@ def run(self, region, center):
Region to rotate
center : `~astropy.coordinates.SkyCoord`
Rotation point
exclusion_mask : `~gammapy.maps.WcsNDMap`, optional
Exclusion mask
Returns
-------
Expand Down Expand Up @@ -65,8 +67,6 @@ class ReflectedRegionsFinder(RegionsFinder):
Minimal distance from input region
max_region_number : int, optional
Maximum number of regions to use
exclusion_mask : `~gammapy.maps.WcsNDMap`, optional
Exclusion mask
binsz : `~astropy.coordinates.Angle`
Bin size of the reference map used for region finding.
Expand Down Expand Up @@ -94,7 +94,6 @@ def __init__(
min_distance="0 rad",
min_distance_input="0.1 rad",
max_region_number=10000,
exclusion_mask=None,
binsz="0.01 deg",
):

Expand All @@ -106,10 +105,6 @@ def __init__(
self.min_distance = Angle(min_distance)
self.min_distance_input = Angle(min_distance_input)

if exclusion_mask and not exclusion_mask.is_mask:
raise ValueError("Exclusion mask must contain boolean values")

self.exclusion_mask = exclusion_mask
self.max_region_number = max_region_number
self.binsz = Angle(binsz)

Expand Down Expand Up @@ -172,18 +167,22 @@ def _region_angular_size(region, reference_geom, center_pix):

return angular_size

def _exclusion_mask_ref(self, reference_geom):
@staticmethod
def _exclusion_mask_ref(reference_geom, exclusion_mask):
"""Exclusion mask reprojected"""
if self.exclusion_mask:
mask = self.exclusion_mask.interp_to_geom(reference_geom, fill_value=True)
if exclusion_mask:
mask = exclusion_mask.interp_to_geom(reference_geom, fill_value=True)
else:
mask = WcsNDMap.from_geom(geom=reference_geom, data=True)
return mask

def _get_excluded_pixels(self, reference_geom):
@staticmethod
def _get_excluded_pixels(reference_geom, exclusion_mask):
"""Excluded pix coords"""
# find excluded PixCoords
exclusion_mask = self._exclusion_mask_ref(reference_geom)
exclusion_mask = ReflectedRegionsFinder._exclusion_mask_ref(
reference_geom, exclusion_mask,
)
pix_y, pix_x = np.nonzero(~exclusion_mask.data)
return PixCoord(pix_x, pix_y)

Expand All @@ -198,7 +197,7 @@ def _get_angle_range(self, region, reference_geom, center_pix):
angle_max = FULL_CIRCLE - angle_min - self.min_distance_input
return angle_min, angle_max

def run(self, region, center):
def run(self, region, center, exclusion_mask=None):
"""Find reflected regions.
Parameters
Expand All @@ -224,7 +223,7 @@ def run(self, region, center):
)

region_pix = self._get_region_pixels(region, reference_geom)
excluded_pixels = self._get_excluded_pixels(reference_geom)
excluded_pixels = self._get_excluded_pixels(reference_geom, exclusion_mask)

angle = angle_min + self.min_distance_input
while angle < angle_max:
Expand All @@ -247,21 +246,34 @@ def run(self, region, center):
class ReflectedRegionsBackgroundMaker(Maker):
"""Reflected regions background maker.
Parameters
Attributes
----------
regions_finder: RegionsFinder
if not given, a `ReflectedRegionsFinder` will be created and
any of the ``**kwargs`` will be forwarded to the `ReflectedRegionsFinder`.
exclusion_mask : `~gammapy.maps.WcsNDMap`, optional
Exclusion mask
"""

tag = "ReflectedRegionsBackgroundMaker"

def __init__(
self,
region_finder=None,
exclusion_mask=None,
**kwargs,
):

if exclusion_mask and not exclusion_mask.is_mask:
raise ValueError("Exclusion mask must contain boolean values")

self.exclusion_mask = exclusion_mask

if region_finder is None:
self.region_finder = ReflectedRegionsFinder(**kwargs)
else:
if len(kwargs) != 0:
raise ValueError('No kwargs can be given if providing a region_finder')
self.region_finder = region_finder

def make_counts_off(self, dataset, observation):
Expand Down Expand Up @@ -296,12 +308,14 @@ def make_counts_off(self, dataset, observation):
rad_max=observation.rad_max,
events=observation.events,
region_finder=self.region_finder,
exclusion_mask=self.exclusion_mask,
)

else:
regions, wcs = self.region_finder.run(
center=observation.pointing_radec,
region=dataset.counts.geom.region,
exclusion_mask=self.exclusion_mask,
)

energy_axis = dataset.counts.geom.axes["energy"]
Expand Down
28 changes: 19 additions & 9 deletions gammapy/makers/background/tests/test_reflected.py
Expand Up @@ -54,8 +54,11 @@ def observations():

@pytest.fixture()
def reflected_bkg_maker(exclusion_mask):
finder = ReflectedRegionsFinder(exclusion_mask=exclusion_mask)
return ReflectedRegionsBackgroundMaker(region_finder=finder)
finder = ReflectedRegionsFinder()
return ReflectedRegionsBackgroundMaker(
region_finder=finder,
exclusion_mask=exclusion_mask,
)


region_finder_param = [
Expand All @@ -74,27 +77,33 @@ def test_find_reflected_regions(
):
pointing = pointing_pos
finder = ReflectedRegionsFinder(
exclusion_mask=exclusion_mask,
min_distance_input="0 deg",
)
regions, _ = finder.run(center=pointing, region=on_region)
regions, _ = finder.run(
center=pointing, region=on_region,
exclusion_mask=exclusion_mask,
)
assert len(regions) == nreg1
assert_quantity_allclose(regions[3].center.icrs.ra, reg3_ra, rtol=1e-2)

# Test without exclusion
finder.exclusion_mask = None
regions, _ = finder.run(center=pointing, region=on_region)
assert len(regions) == nreg2

# Test with too small exclusion
small_mask = exclusion_mask.cutout(pointing, Angle("0.1 deg"))
finder.exclusion_mask = small_mask
regions, _ = finder.run(center=pointing, region=on_region)
regions, _ = finder.run(
center=pointing, region=on_region,
exclusion_mask=small_mask,
)
assert len(regions) == nreg3

# Test with maximum number of regions
finder.max_region_number = 5
regions, _ = finder.run(center=pointing, region=on_region)
regions, _ = finder.run(
center=pointing, region=on_region,
exclusion_mask=small_mask,
)
assert len(regions) == 5

# Test with an other type of region
Expand All @@ -109,6 +118,7 @@ def test_find_reflected_regions(
regions, _ = finder.run(
region=on_ellipse_annulus,
center=pointing,
exclusion_mask=small_mask,
)
assert len(regions) == 5

Expand Down Expand Up @@ -139,12 +149,12 @@ def test_non_circular_regions(region, nreg):
def test_bad_on_region(exclusion_mask, on_region):
pointing = SkyCoord(83.63, 22.01, unit="deg", frame="icrs")
finder = ReflectedRegionsFinder(
exclusion_mask=exclusion_mask,
min_distance_input="0 deg",
)
regions, _ = finder.run(
center=pointing,
region=on_region,
exclusion_mask=exclusion_mask,
)
assert len(regions) == 0

Expand Down
2 changes: 2 additions & 0 deletions gammapy/makers/utils.py
Expand Up @@ -495,6 +495,7 @@ def make_counts_off_rad_max(
rad_max,
events,
region_finder,
exclusion_mask=None,
):
"""Extract the OFF counts and the ON / OFF acceptance considering for the
sizes of the ON and OFF regions the values in the `RAD_MAX_2D` table.
Expand Down Expand Up @@ -532,6 +533,7 @@ def make_counts_off_rad_max(
regions, wcs = region_finder.run(
center=events.pointing_radec,
region=on_region,
exclusion_mask=exclusion_mask,
)

if len(regions) > 0:
Expand Down

0 comments on commit e4536a9

Please sign in to comment.