From 1fd215e0739fb81f6366bb2ca2bf26937e472194 Mon Sep 17 00:00:00 2001 From: Christoph Deil Date: Fri, 23 May 2014 11:37:35 +0200 Subject: [PATCH] Work on BoundingBox code --- gammapy/image/measure.py | 131 +++++++++++++++++++++++++++++++++------ gammapy/image/utils.py | 46 +++++++------- 2 files changed, 136 insertions(+), 41 deletions(-) diff --git a/gammapy/image/measure.py b/gammapy/image/measure.py index 6bee1d38adc..2d9117d28d7 100644 --- a/gammapy/image/measure.py +++ b/gammapy/image/measure.py @@ -19,47 +19,142 @@ class BoundingBox(object): """Rectangular bounding box. + TODO: Link to bounding box docs. + Parameters ---------- - TODO + x_start, y_start : float + Start pixel coordinates (0-based indexing, inclusive) + x_stop, y_stop : float + Stop pixel coordinates (0-based indexing, exclusive) """ - def __init__(self, xmin, xmax, ymin, ymax): - self.xmin = xmin - self.xmax = xmax - self.ymin = ymin - self.ymax = ymax + def __init__(self, x_start, y_start, x_stop, y_stop): + self.x_start = x_start + self.y_start = y_start + self.x_stop = x_stop + self.y_stop = y_stop @staticmethod - def for_circle(x, y, radius, nx, ny): + def for_circle(x, y, radius): """Compute bounding box for a circle. Parameters ---------- + x, y : float + Circle center position + radius : float + Circle radius (pix) + + Returns + ------- + bounding_box : BoundingBox + Bounding box + """ + x, y, radius = int(x), int(y), int(radius) + x_start = x - radius + y_start = y - radius + x_stop = x + radius + y_stop = y + radius + return BoundingBox(x_start, y_start, x_stop, y_stop) + + + def __str__(self): + ss = 'x = ({x_start}, {x_stop}), '.format(**self) + ss += 'y = ({y_start}, {y_stop}) '.format(**self) + return ss + + def intersection(self, other): + """Compute intersection of two bounding boxes. + + Parameters + ---------- + other : `BoundingBox` + Other bounding box. + + Returns + ------- + intersection : `BoundingBox` + New intersection bounding box + + Examples + -------- + >>> b1 = BoundingBox(1, 5, 10, 11) + >>> b2 = BoundingBox(2, 4, 13, 9) + >>> b1.intersection(b2) TODO + """ + x_start = max(self.x_start, other.x_start) + x_stop = min(self.x_stop, other.x_stop) + y_start = max(self.y_start, other.y_start) + y_stop = min(self.x_stop, other.x_stop) + return BoundingBox(x_start, y_start, x_stop, y_stop) + + + def overlaps(self, other): + """Does this bounding box overlap the other bounding box? + + Parameters + ---------- + other : `BoundingBox` + Other bounding box Returns ------- + overlaps : bool + True of there is overlap + + Examples + -------- TODO """ - x, y, radius = int(x), int(y), int(radius) - xmin = max(x - radius, 0) - xmax = min(x + radius, nx) - ymin = max(y - radius, 0) - ymax = min(y + radius, ny) - bbox = BoundingBox(xmin, xmax, ymin, ymax) + return self.intersection(other).is_empty + - return bbox + def __eq__(self, other): + return self.slice == other.slice + + @property + def ftcopy_string(self): + """Bounding box in ftcopy string format. - def to_ftcopy_string(self): - """Create bounding box string in ftcopy format. + Examples + -------- + >>> from gammapy.image import BoundingBox + >>> bounding_box = BoundingBox(7, 9, 43, 44) + >>> bounding_box.ftcopy_string + TODO + """ + return '[{x_start}:{x_stop},{y_start}:{y_stop}]'.format(**self) + + + @property + def slice(self): + """Bounding box in slice format. Examples -------- + >>> from gammapy.image import BoundingBox + >>> bounding_box = BoundingBox(7, 9, 43, 44) + >>> bounding_box.slice TODO - """ - return '[{xmin}:{xmax},{ymin}:{ymax}]'.format(**self) + """ + return (self.y_slice, self.x_slice) + + @property + def x_slice(self): + return slice(self.x_start, self.x_stop) + + @property + def y_slice(self): + return slice(self.y_start, self.y_stop) + + @property + def is_empty(self): + x_is_empty = (self.x_start >= self.x_stop) + y_is_empty = (self.y_start >= self.y_stop) + return x_is_empty or y_is_empty def bbox(mask, margin, binsz): diff --git a/gammapy/image/utils.py b/gammapy/image/utils.py index 2f00de896a3..6f71e769dd4 100644 --- a/gammapy/image/utils.py +++ b/gammapy/image/utils.py @@ -14,7 +14,7 @@ 'binary_opening_circle', 'binary_ring', 'calc_footprint', 'contains', 'coordinates', 'cube_to_image', 'cube_to_spec', - 'cut_out', + 'crop_image', 'disk_correlate', 'exclusion_distance', 'image_groupby', 'images_to_cube', 'make_empty_image', 'make_header', @@ -697,37 +697,32 @@ def make_empty_image(nxpix=100, nypix=100, binsz=0.1, xref=0, yref=0, fill=0, data = fill * np.ones(shape, dtype=dtype) return fits.ImageHDU(data, header) -def cut_out(image, center, fov=[2, 2]): - """Cut out a part of an image. + +def crop_image(image, bounding_box): + """Crop an image (cut out a rectangular part). Parameters ---------- image : `astropy.io.fits.ImageHDU` - Input image - center : [glon, glat] - Cutout center position - fov : [glon_fov, glat_fov] - Cutout field of view dimensions + Image + bounding_box : `~gammapy.image.BoundingBox` + Bounding box Returns ------- - TODO + new_image : `astropy.io.fits.ImageHDU` + Cropped image + + See Also + -------- + paste_cutout_into_image """ - raise NotImplementedError - - # Unpack center and fov - glon, glat = center - glon_fov, glat_fov = fov + data = image.data[bounding_box.slice] + header = image.header.copy() - # Calculate image limits - glon_lim = [glon - glon_fov, glon + glon_fov] - glat_lim = [glat - glat_fov, glat + glat_fov] + # TODO: fix header keywords and test against ftcopy - # Cut out the requested part - xlim = image.proj.topixel((glon_lim, [0, 0]))[0] - xlim = [xlim[1], xlim[0]] # longitude axis runs backwards - ylim = image.proj.topixel(([0, 0], glat_lim))[1] - image.set_limits(pxlim=xlim, pylim=ylim) + return fits.ImageHDU(data=data, header=header) def cube_to_image(cube, slicepos=None): @@ -825,13 +820,16 @@ def paste_cutout_into_image(total, cutout, method='sum'): ------- total : `astropy.io.fits.ImageHDU` A reference to the total input HDU that was modified in-place. + + See Also + -------- + crop_image """ # find offset lon, lat = WCS(cutout.header).wcs_pix2world(0, 0, 0) x, y = WCS(total.header).wcs_world2pix(lon, lat, 0) x, y = int(np.round(x)), int(np.round(y)) dy, dx = cutout.shape - # import IPython; IPython.embed(); 1/0 if method == 'sum': total.data[y : y + dy, x : x + dx] += cutout.data @@ -840,6 +838,8 @@ def paste_cutout_into_image(total, cutout, method='sum'): else: raise ValueError('Invalid method: {0}'.format(method)) + return total + def block_reduce_hdu(input_hdu, block_size, func, cval=0): """Provides block reduce functionality for image HDUs.