Make regions

README.rst

@@ -41,7 +41,7 @@ To show help message and exit
                       [--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]
+                      [-o OUTFILE] [--save-regions OUTREGION] [--gui]
 
      breizorro [options] --restored-image restored_image
 
@@ -86,7 +86,9 @@ To show help message and exit
                                 ratio above 100
           -o OUTFILE, --outfile OUTFILE
                                 Suffix for mask image (default based on input name
-          --gui                 Open mask in gui.
+          --save-regions OUTREGION
+                                Generate polygon regions from the mask
+          --gui                 Open mask in gui
 
 =======
 License

breizorro/breizorro.py

@@ -3,20 +3,22 @@
 import shutil
 import logging
 import argparse
+import regions
 import os.path
 import re
 import numpy as np
 from astropy.io import fits
 from astropy.coordinates import SkyCoord
 from astropy.wcs import WCS
-import regions
 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 regions import PolygonSkyRegion
+from skimage.measure import find_contours
 
 def create_logger():
     """Create a console logger"""
@@ -156,7 +158,7 @@ def main():
                          help='Replace all island numbers with 1')
     parser.add_argument('--invert', action="store_true",
                          help='Invert the mask')
-    
+
     parser.add_argument('--dilate', dest='dilate', metavar="R", type=int, default=0,
                         help='Apply dilation with a radius of R pixels')
     parser.add_argument('--erode', dest='erode', metavar="N", type=int, default=0,
@@ -170,6 +172,8 @@ def main():
 
     parser.add_argument('-o', '--outfile', dest='outfile', default='',
                         help='Suffix for mask image (default based on input name')
+    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.')
@@ -343,6 +347,23 @@ def load_fits_or_region(filename):
         mask_image = input_image * new_mask_image
         LOGGER.info(f"Number of extended islands found: {len(extended_islands)}")
 
+    if args.outregion:
+        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)
+        regions.Regions(polygon_regions).write(args.outregion, format='ds9')
+        LOGGER.info(f"Number of regions found: {len(polygon_regions)}")
+        LOGGER.info(f"Saving regions in {args.outregion}")
+
+
     if args.gui:
         try:
             from bokeh.models import BoxEditTool, ColumnDataSource, FreehandDrawTool

pyproject.toml

@@ -29,6 +29,7 @@ astropy = "*"
 numpy = "*"
 regions = "*"
 scipy = "*"
+scikit-image = "*"
 
 # Optional dependencies start here
 bokeh = { version = "*", optional = true}
@@ -54,4 +55,4 @@ optional = true
 [tool.poetry.group.docs.dependencies]
 Sphinx = "^5.3.0"
 sphinx-copybutton = "^0.5.0"
-furo = "^2022.9.15"
+furo = "^2022.9.15"