Write regions and catalogs

README.rst

@@ -30,63 +30,11 @@ To show help message and exit
 
    $ breizorro --help
 
-     breizorro.breizorro - 2022-08-24 11:07:39,311 INFO - Welcome to breizorro
-     breizorro.breizorro - 2022-08-24 11:07:39,375 INFO - Version: 0.1.1
-     breizorro.breizorro - 2022-08-24 11:07:39,375 INFO - Usage: breizorro --help
-     usage: breizorro [-h] [-r IMAGE] [-m MASK] [-t THRESHOLD] [-b BOXSIZE]
-                      [--savenoise] [--merge MASKs|REGs) [MASK(s|REGs) ...]]
-                      [--subtract MASK(s|REGs) [MASK(s|REGs ...]]
-                      [--number-islands] [--remove-islands N|COORD [N|COORD ...]]
-                      [--ignore-missing-islands]
-                      [--extract-islands N|COORD [N|COORD ...]]
-                      [--minimum-size MINSIZE] [--make-binary] [--invert]
-                      [--dilate R] [--erode N] [--fill-holes] [--sum-peak SUM_PEAK]
-                      [-o OUTFILE] [--gui]
-
-     breizorro [options] --restored-image restored_image
-
-     optional arguments:
-          -h, --help            show this help message and exit
-          -r IMAGE, --restored-image IMAGE
-                                Restored image file from which to build mask
-          -m MASK, --mask-image MASK
-                                Input mask file(s). Either --restored-image or --mask-
-                                image must be specfied.
-          -t THRESHOLD, --threshold THRESHOLD
-                                Sigma threshold for masking (default = 6.5)
-          -b BOXSIZE, --boxsize BOXSIZE
-                                Box size over which to compute stats (default = 50)
-          --savenoise           Enable to export noise image as FITS file (default=do
-                                not save noise image)
-          --merge MASK(s)|REG(s) [MASK(s)|REG(s) ...]
-                                Merge in one or more masks or region files
-          --subtract MASK(s)|REG(s) [MASK(s)|REG(s) ...]
-                                Subract one or more masks or region files
-          --number-islands      Number the islands detected (default=do not number
-                                islands)
-          --remove-islands N|COORD [N|COORD ...]
-                                List of islands to remove from input mask. e.g.
-                                --remove-islands 1 18 20 20h10m13s,14d15m20s
-          --ignore-missing-islands
-                                If an island specified by coordinates does not exist,
-                                do not throw an error
-          --extract-islands N|COORD [N|COORD ...]
-                                List of islands to extract from input mask. e.g.
-                                --extract-islands 1 18 20 20h10m13s,14d15m20s
-          --minimum-size MINSIZE
-                                Remove islands that have areas fewer than or equal to
-                                the specified number of pixels
-          --make-binary         Replace all island numbers with 1
-          --invert              Invert the mask
-          --dilate R            Apply dilation with a radius of R pixels
-          --erode N             Apply N iterations of erosion
-          --fill-holes          Fill holes (i.e. entirely closed regions) in mask
-          --sum-peak SUM_PEAK   Sum to peak ratio of flux islands to mask in original
-                                image.e.g. --sum-peak 100 will mask everything with a
-                                ratio above 100
-          -o OUTFILE, --outfile OUTFILE
-                                Suffix for mask image (default based on input name
-          --gui                 Open mask in gui.
+=============
+Documentation
+=============
+
+Documentation is available on the GitHub page_.
 
 =======
 License
@@ -116,3 +64,4 @@ standards pep8_.
 .. _source: https://github.com/ratt-ru/breizorro
 .. _license: https://github.com/ratt-ru/breizorro/blob/main/LICENSE
 .. _pep8: https://www.python.org/dev/peps/pep-0008
+.. _page: https://ratt-ru.github.io/breizorro

_config.yml

@@ -0,0 +1,5 @@
+title: Breizorro
+description: Image masking and cataloguing suite
+show_downloads: true
+google_analytics: 
+theme: jekyll-theme-midnight

breizorro/breizorro.py

@@ -6,17 +6,27 @@
 import os.path
 import re
 import numpy as np
+
+import astropy.units as u
 from astropy.io import fits
-from astropy.coordinates import SkyCoord
 from astropy.wcs import WCS
+from astropy.coordinates import Angle
+from astropy.coordinates import SkyCoord
+
 import regions
+from regions import PixCoord
+from regions import PolygonSkyRegion, PolygonPixelRegion
+
 from argparse import ArgumentParser
 
 from scipy.ndimage.morphology import binary_dilation, binary_erosion, binary_fill_holes
 from scipy.ndimage.measurements import label, find_objects
 import scipy.special
 import scipy.ndimage
 
+from breizorro.utils import get_source_size, format_source_coordinates, deg2ra, deg2dec
+from breizorro.utils import get_image_data, fitsInfo, calculate_beam_area
+
 
 def create_logger():
     """Create a console logger"""
@@ -31,40 +41,20 @@ def create_logger():
 
 LOGGER = create_logger()
 
-def get_image(fitsfile):
-    """
-    Reads FITS file, returns tuple of image_array, header
-    """
-    LOGGER.info(f"Reading {fitsfile} data")
-    input_hdu = fits.open(fitsfile)[0]
-    if len(input_hdu.data.shape) == 2:
-        image = numpy.array(input_hdu.data[:, :])
-    elif len(input_hdu.data.shape) == 3:
-        image = numpy.array(input_hdu.data[0, :, :])
-    else:
-        image = numpy.array(input_hdu.data[0, 0, :, :])
-    return image, input_hdu.header
-
-
-def get_image_header(fitsfile):
-    LOGGER.info(f"Reading {fitsfile} header")
-    input_hdu = fits.open(fitsfile)[0]
-    return input_hdu.header
-
 
 def flush_fits(newimage, fitsfile, header=None):
     LOGGER.info(f"Writing {fitsfile}")
-    f = fits.open(fitsfile, mode='update')
-    input_hdu = f[0]
-    if len(input_hdu.data.shape) == 2:
-        input_hdu.data[:, :] = newimage
-    elif len(input_hdu.data.shape) == 3:
-        input_hdu.data[0, :, :] = newimage
-    else:
-        input_hdu.data[0, 0, :, :] = newimage
-    if header:
-        input_hdu.header = header
-    f.flush()
+    with fits.open(fitsfile, mode='update') as f:
+        input_hdu = f[0]
+        if len(input_hdu.data.shape) == 2:
+            input_hdu.data[:, :] = newimage
+        elif len(input_hdu.data.shape) == 3:
+            input_hdu.data[0, :, :] = newimage
+        else:
+            input_hdu.data[0, 0, :, :] = newimage
+        if header:
+            input_hdu.header = header
+        f.flush()
 
 
 def make_noise_map(restored_image, boxsize):
@@ -86,6 +76,7 @@ def make_noise_map(restored_image, boxsize):
     LOGGER.info(f"Median noise value is {median_noise}")
     return noise
 
+
 def resolve_island(isl_spec, mask_image, wcs, ignore_missing=False):
     if re.match("^\d+$", isl_spec):
         return int(isl_spec)
@@ -104,12 +95,14 @@ def resolve_island(isl_spec, mask_image, wcs, ignore_missing=False):
             raise ValueError(f"coordinates {c} do not select a valid island")
     return value
 
+
 def add_regions(mask_image, regs, wcs):
     for reg in regs:
         if hasattr(reg, 'to_pixel'):
             reg = reg.to_pixel(wcs)
         mask_image += reg.to_mask().to_image(mask_image.shape)
 
+
 def remove_regions(mask_image, regs, wcs):
     for reg in regs:
         if hasattr(reg, 'to_pixel'):
@@ -167,9 +160,18 @@ def main():
     parser.add_argument('--sum-peak', dest='sum_peak', default=None,
                         help='Sum to peak ratio of flux islands to mask in original image.'
                              'e.g. --sum-peak 100 will mask everything with a ratio above 100')
+    parser.add_argument('-ncpu', '--ncpu', dest='ncpu', default=None, type=int,
+                        help='Number of processors to use for cataloging.')
+    parser.add_argument('-beam', '--beam-size', dest='beam', default=None,
+                        help='Average beam size in arcesc incase beam info is missing in image header.')
 
     parser.add_argument('-o', '--outfile', dest='outfile', default='',
                         help='Suffix for mask image (default based on input name')
+    parser.add_argument('--save-catalog', dest='outcatalog', default='',
+                        help='Generate catalog based on region mask')
+    parser.add_argument('--save-regions', dest='outregion', default='',
+                         help='Generate polygon regions from the mask')
+
 
     parser.add_argument('--gui', dest='gui', action='store_true', default=False,
                          help='Open mask in gui.')
@@ -192,7 +194,7 @@ def main():
 
     # first, load or generate mask
     if args.imagename:
-        input_image, input_header = get_image(input_file)
+        input_image, input_header = get_image_data(input_file)
         LOGGER.info(f"Generating mask using threshold {threshold}")
 
         noise_image = make_noise_map(input_image, boxsize)
@@ -214,7 +216,7 @@ def main():
         out_mask_fits = args.outfile or f"{name}.mask.fits"
 
     elif args.maskname:
-        mask_image, mask_header = get_image(args.maskname)
+        mask_image, mask_header = get_image_data(args.maskname)
         LOGGER.info(f"Input mask loaded")
 
         out_mask_fits = args.outfile or f"{name}.out.{ext}"
@@ -231,7 +233,7 @@ def load_fits_or_region(filename):
         fits = regs = None
         # read as FITS or region
         try:
-            fits = get_image(filename)
+            fits = get_image_data(filename)
         except OSError:
             try:
                 regs = regions.Regions.read(filename)
@@ -342,40 +344,100 @@ def load_fits_or_region(filename):
         mask_header['BUNIT'] = 'Jy/beam'
         mask_image = input_image * new_mask_image
         LOGGER.info(f"Number of extended islands found: {len(extended_islands)}")
+        shutil.copyfile(input_file, out_mask_fits)  # to provide a template
+        flush_fits(mask_image, out_mask_fits, mask_header)
+        LOGGER.info("Done")
+        sys.exit(1)
+
+    if args.outcatalog or args.outregion:
+        try:
+            from skimage.measure import find_contours
+        except ImportError:
+            LOGGER.info('pip install breizorro[all] to use cataloguing feature.')
+            exit(1)
+        contours = find_contours(mask_image, 0.5)
+        polygon_regions = []
+        for contour in contours:
+            # Convert the contour points to pixel coordinates
+            contour_pixels = contour
+            # Convert the pixel coordinates to Sky coordinates
+            contour_sky = wcs.pixel_to_world(contour_pixels[:, 1], contour_pixels[:, 0])
+            # Create a Polygon region from the Sky coordinates
+            polygon_region = PolygonSkyRegion(vertices=contour_sky, meta={'label': 'Region'})
+            # Add the polygon region to the list
+            polygon_regions.append(polygon_region)
+        LOGGER.info(f"Number of regions found: {len(polygon_regions)}")
+        if args.outregion:
+            regions.Regions(polygon_regions).write(args.outregion, format='ds9')
+            LOGGER.info(f"Saving regions in {args.outregion}")
+
+    if args.outcatalog and args.imagename:
+        try:
+            import warnings
+            # Suppress FittingWarnings from Astropy
+            # WARNING: The fit may be unsuccessful; check fit_info['message'] for more information. [astropy.modeling.fitting]
+            warnings.resetwarnings()
+            warnings.filterwarnings('ignore', category=UserWarning, append=True)
+            from breizorro.catalog import multiprocess_contours
+        except ModuleNotFoundError:
+            LOGGER.error("Running breizorro source detector requires optional dependencies, please re-install with: pip install breizorro[catalog]")
+            raise('Missing cataloguing dependencies')
+        source_list = []
+        image_data, hdu_header = get_image_data(args.imagename)
+        fitsinfo = fitsInfo(args.imagename)
+        mean_beam = None # Use beam info from the image header by default
+        if args.beam:
+            mean_beam = float(args.beam)
+        if mean_beam:
+            LOGGER.info(f'Using user provided size: {mean_beam}')
+        elif fitsinfo['b_size']:
+            bmaj,bmin,_ = np.array(fitsinfo['b_size']) * 3600.0
+            mean_beam = 0.5 * (bmaj + bmin)
+        else:
+            raise('No beam information found. Specify mean beam in arcsec: --beam-size 6.5')
+
+        noise = np.median(noise_image)
+        f = open(args.outcatalog, 'w')
+        catalog_out = f'# processing fits image: {args.imagename}  \n'
+        f.write(catalog_out)
+        image_dimensions = fitsinfo['naxis']
+        pixel_size = fitsinfo['ddec'] * 3600.0
+        catalog_out = f'# mean beam size (arcsec): {round(mean_beam,2)} \n' 
+        f.write(catalog_out)
+        catalog_out = f'# original image peak flux (Jy/beam): {image_data.max()} \n'
+        f.write(catalog_out)
+        catalog_out = f'# noise out (µJy/beam): {round(noise*1000000,2)} \n'
+        f.write(catalog_out)
+        limiting_flux = noise * threshold
+        catalog_out = f'# cutt-off flux  (mJy/beam): {round(limiting_flux*1000,2)} \n'
+        f.write(catalog_out)
+        LOGGER.info('Submitting distributed tasks. This might take a while...') 
+        source_list = multiprocess_contours(contours, image_data, fitsinfo, noise, args.ncpu)
+        catalog_out = f"# freq0 (Hz): {fitsinfo['freq0']} \n"
+        f.write(catalog_out)
+        catalog_out = f'# number of sources detected: {len(source_list)} \n'
+        f.write(catalog_out)
+        catalog_out = '#\n#format: name ra_d dec_d i i_err emaj_s emin_s pa_d\n'
+        f.write(catalog_out)
+        for i in range(len(source_list)):
+            output = 'src' + str(i) + ' ' + source_list[i][1] + '\n'
+            f.write(output)
+        f.close()
+        LOGGER.info(f'Source catalog saved: {args.outcatalog}')
 
     if args.gui:
         try:
-            from bokeh.models import BoxEditTool, ColumnDataSource, FreehandDrawTool
-            from bokeh.plotting import figure, output_file, show
-            from bokeh.io import curdoc
-            curdoc().theme = 'caliber'
+            from breizorro.gui import display
         except ModuleNotFoundError:
             LOGGER.error("Running breizorro gui requires optional dependencies, please re-install with: pip install breizorro[gui]")
+            raise('Missing GUI dependencies')
 
         LOGGER.info("Loading Gui ...")
-        d = mask_image
-        p = figure(tooltips=[("x", "$x"), ("y", "$y"), ("value", "@image")])
-        p.x_range.range_padding = p.y_range.range_padding = 0
-        p.title.text = out_mask_fits
-
-        # must give a vector of image data for image parameter
-        p.image(image=[d], x=0, y=0, dw=10, dh=10, palette="Greys256", level="image")
-        p.grid.grid_line_width = 0.5
-        src1 = ColumnDataSource({'x': [], 'y': [], 'width': [], 'height': [], 'alpha': []})
-        src2 = ColumnDataSource({'xs': [], 'ys': [], 'alpha': []})
-        renderer1 = p.rect('x', 'y', 'width', 'height', source=src1, alpha='alpha')
-        renderer2 = p.multi_line('xs', 'ys', source=src2, alpha='alpha')
-        draw_tool1 = BoxEditTool(renderers=[renderer1], empty_value=1)
-        draw_tool2 = FreehandDrawTool(renderers=[renderer2])
-        p.add_tools(draw_tool1)
-        p.add_tools(draw_tool2)
-        p.toolbar.active_drag = draw_tool1
-        output_file("breizorro.html", title="Mask Editor")
-        show(p)
-
-        LOGGER.info(f"Enforcing that mask to binary")
-        mask_image = mask_image!=0
-        mask_header['BUNIT'] = 'mask'
+        display(args.imagename or args.maskname, mask_image, args.outcatalog, source_list)
+
+    LOGGER.info(f"Enforcing that mask to binary")
+    mask_image = mask_image!=0
+    mask_header['BUNIT'] = 'mask'
 
     shutil.copyfile(input_file, out_mask_fits)  # to provide a template
     flush_fits(mask_image, out_mask_fits, mask_header)