Work on BoundingBox code

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):

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.