diff --git a/.github/workflows/test_installation.yml b/.github/workflows/python-package.yml similarity index 95% rename from .github/workflows/test_installation.yml rename to .github/workflows/python-package.yml index 4c89794..e62236a 100644 --- a/.github/workflows/test_installation.yml +++ b/.github/workflows/python-package.yml @@ -18,7 +18,7 @@ jobs: runs-on: ubuntu-latest strategy: matrix: - python-version: ["3.8", "3.9", "3.10"] + python-version: ["3.9", "3.10", "3.11", "3.12"] steps: - uses: actions/checkout@v2 @@ -64,4 +64,4 @@ jobs: - name: Test command run: | - poetry run breizorro -h + poetry run breizorro --help diff --git a/MANIFEST.in b/MANIFEST.in index a308077..a012c54 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -1,3 +1,4 @@ include LICENSE include README.rst include requirements.txt +include breizorro/breizorro.yaml diff --git a/README.rst b/README.rst index bf96a67..390605a 100644 --- a/README.rst +++ b/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 diff --git a/_config.yml b/_config.yml new file mode 100644 index 0000000..ee9483e --- /dev/null +++ b/_config.yml @@ -0,0 +1,5 @@ +title: Breizorro +description: Image masking and cataloguing suite +show_downloads: true +google_analytics: +theme: jekyll-theme-midnight diff --git a/breizorro.html b/breizorro.html new file mode 100644 index 0000000..e43f14a --- /dev/null +++ b/breizorro.html @@ -0,0 +1,63 @@ + + + + + Mask Editor + + + + + + + +
+ + + + + \ No newline at end of file diff --git a/breizorro/__init__.py b/breizorro/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/breizorro/breizorro.py b/breizorro/breizorro.py index 3c32455..32f0aa4 100644 --- a/breizorro/breizorro.py +++ b/breizorro/breizorro.py @@ -6,10 +6,17 @@ 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 @@ -17,6 +24,9 @@ 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,12 +76,13 @@ 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) - elif ',' not in isl_spec: + elif ':' not in isl_spec: raise ValueError(f"invalid island specification: {isl_spec}") - c = SkyCoord(*isl_spec.split(',', 1)) + c = SkyCoord(*isl_spec.split(':', 1)) x, y = wcs.world_to_pixel(c) x = round(float(x)) y = round(float(y)) @@ -104,95 +95,44 @@ 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'): reg = reg.to_pixel(wcs) mask_image[reg.to_mask().to_image(mask_image.shape) != 0] = 0 - -def main(): +def main(restored_image, mask_image, threshold, boxsize, savenoise, merge, subtract, + number_islands, remove_islands, ignore_missing_islands, extract_islands, + minimum_size, make_binary, invert, dilate, erode, fill_holes, sum_peak, + ncpu, beam_size, gui, outfile, outcatalog, outregion): LOGGER.info("Welcome to breizorro") # Get version from pkg_resources import get_distribution _version = get_distribution('breizorro').version LOGGER.info(f"Version: {_version}") - LOGGER.info("Usage: breizorro --help") - parser = ArgumentParser(description='breizorro [options] --restored-image restored_image') - parser.add_argument('-r', '--restored-image', dest="imagename", metavar="IMAGE", - help="Restored image file from which to build mask") - parser.add_argument('-m', '--mask-image', dest="maskname", metavar="MASK", - help="Input mask file(s). Either --restored-image or --mask-image must be specfied.") - parser.add_argument('-t', '--threshold', dest='threshold', default=6.5, - help='Sigma threshold for masking (default = 6.5)') - parser.add_argument('-b', '--boxsize', dest='boxsize', default=50, - help='Box size over which to compute stats (default = 50)') - parser.add_argument('--savenoise', dest='savenoise', action='store_true', default=False, - help='Enable to export noise image as FITS file (default=do not save noise image)') - - parser.add_argument('--merge', dest='merge', metavar="MASK(s)|REG(s)", nargs='+', - help='Merge in one or more masks or region files') - parser.add_argument('--subtract', dest='subtract', metavar="MASK(s)|REG(s)", nargs='+', - help='Subract one or more masks or region files') - - parser.add_argument('--number-islands', dest='islands', action='store_true', default=False, - help='Number the islands detected (default=do not number islands)') - parser.add_argument('--remove-islands', dest='remove_isl', metavar='N|COORD', type=str, nargs='+', - help='List of islands to remove from input mask. e.g. --remove-islands 1 18 20 20h10m13s,14d15m20s') - parser.add_argument('--ignore-missing-islands', action='store_true', - help='If an island specified by coordinates does not exist, do not throw an error') - parser.add_argument('--extract-islands', dest='extract_isl', metavar='N|COORD', type=str, nargs='+', - help='List of islands to extract from input mask. e.g. --extract-islands 1 18 20 20h10m13s,14d15m20s') - parser.add_argument('--minimum-size', dest='minsize', type=int, - help='Remove islands that have areas fewer than or equal to the specified number of pixels') - parser.add_argument('--make-binary', action="store_true", - 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, - help='Apply N iterations of erosion') - parser.add_argument('--fill-holes', dest='fill_holes', action='store_true', - help='Fill holes (i.e. entirely closed regions) in mask') - - 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('-o', '--outfile', dest='outfile', default='', - help='Suffix for mask image (default based on input name') - - parser.add_argument('--gui', dest='gui', action='store_true', default=False, - help='Open mask in gui.') - args = parser.parse_args() - threshold = float(args.threshold) - boxsize = int(args.boxsize) - dilate = int(args.dilate) - savenoise = args.savenoise - outfile = args.outfile - - if args.imagename and args.maskname: - parser.error("Either --restored-image or --mask-image must be specified, but not both") - elif not args.imagename and not args.maskname: - parser.error("Either --restored-image or --mask-image must be specified") + + if restored_image and mask_image: + LOGGER.error("Either --restored-image or --mask-image must be specified, but not both") + elif not restored_image and not mask_image: + LOGGER.error("Either --restored-image or --mask-image must be specified") # define input file, and get its name and extension - input_file = args.imagename or args.maskname + input_file = restored_image or mask_image name = '.'.join(input_file.split('.')[:-1]) ext = input_file.split('.')[-1] # first, load or generate mask - if args.imagename: - input_image, input_header = get_image(input_file) + if restored_image: + 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) @@ -211,15 +151,15 @@ def main(): mask_header = input_header mask_header['BUNIT'] = 'mask' - out_mask_fits = args.outfile or f"{name}.mask.fits" + out_mask_fits = outfile or f"{name}.mask.fits" - elif args.maskname: - mask_image, mask_header = get_image(args.maskname) + elif mask_image: + mask_image, mask_header = get_image_data(mask_image) LOGGER.info(f"Input mask loaded") - out_mask_fits = args.outfile or f"{name}.out.{ext}" + out_mask_fits = outfile or f"{name}.out.{ext}" else: - parser.error("Either --restored-image or --mask-image must be specified") + LOGGER.error("Either --restored-image or --mask-image must be specified") sys.exit(1) wcs = WCS(mask_header) @@ -231,41 +171,41 @@ 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) except: - LOGGER.error(f"{merge} is neither a FITS file not a regions file") - raise + msg = f"{merge} is neither a FITS file not a regions file" + LOGGER.error(msg) + raise(msg) return fits, regs - - if args.merge: - for merge in args.merge: - fits, regs = load_fits_or_region(merge) + if merge: + for _merge in merge: + fits, regs = load_fits_or_region(_merge) if fits: - LOGGER.info(f"Treating {merge} as a FITS mask") + LOGGER.info(f"Treating {_merge} as a FITS mask") mask_image += fits[0] LOGGER.info("Merged into mask") else: - LOGGER.info(f"Merging in {len(regs)} regions from {merge}") + LOGGER.info(f"Merging in {len(regs)} regions from {_merge}") add_regions(mask_image, regs, wcs) mask_image = mask_image != 0 mask_header['BUNIT'] = 'mask' - if args.subtract: - for subtract in args.subtract: - fits, regs = load_fits_or_region(subtract) + if subtract: + for _subtract in subtract: + fits, regs = load_fits_or_region(_subtract) if fits: - LOGGER.info(f"treating {subtract} as a FITS mask") + LOGGER.info(f"treating {_subtract} as a FITS mask") mask_image[fits[0] != 0] = 0 LOGGER.info("Subtracted from mask") else: - LOGGER.info(f"Subtracting {len(regs)} regions from {subtract}") + LOGGER.info(f"Subtracting {len(regs)} regions from {_subtract}") remove_regions(mask_image, regs, wcs) - if args.islands: + if number_islands: LOGGER.info(f"(Re)numbering islands") mask_image = mask_image != 0 # mask_image = mask_image.byteswap().newbyteorder() @@ -273,58 +213,58 @@ def load_fits_or_region(filename): mask_header['BUNIT'] = 'Source_ID' LOGGER.info(f"Number of islands: {num_features}") - if args.remove_isl: - LOGGER.info(f"Removing islands: {args.remove_isl}") - for isl_spec in args.remove_isl: - isl = resolve_island(isl_spec, mask_image, wcs, ignore_missing=args.ignore_missing_islands) + if remove_islands: + LOGGER.info(f"Removing islands: {remove_islands}") + for isl_spec in remove_islands: + isl = resolve_island(isl_spec, mask_image, wcs, ignore_missing=ignore_missing_islands) if isl != None: mask_image[mask_image == isl] = 0 - if args.extract_isl: - LOGGER.info(f"Extracting islands: {args.extract_isl}") + if extract_islands: + LOGGER.info(f"Extracting islands: {extract_islands}") new_mask_image = np.zeros_like(mask_image) - for isl_spec in args.extract_isl: + for isl_spec in extract_islands: isl = resolve_island(isl_spec, mask_image, wcs) new_mask_image[mask_image == isl] = isl mask_image = new_mask_image - if args.minsize: - LOGGER.info(f"Removing islands that occupy fewer than or equal to {args.minsize} pixels") + if minimum_size: + LOGGER.info(f"Removing islands that occupy fewer than or equal to {minimum_size} pixels") mask_image = mask_image != 0 island_labels, num_features = label(mask_image) island_areas = numpy.array(scipy.ndimage.sum(mask_image,island_labels, numpy.arange(island_labels.max()+1))) - min_mask = island_areas >= args.minsize + min_mask = island_areas >= minimum_size mask_image = min_mask[island_labels.ravel()].reshape(island_labels.shape) - if args.make_binary: + if make_binary: LOGGER.info(f"Converting mask to binary") mask_image = mask_image!=0 mask_header['BUNIT'] = 'mask' - if args.invert: + if invert: LOGGER.info(f"Inverting mask") mask_image = mask_image==0 - if args.dilate: - LOGGER.info(f"Dilating mask using a ball of R={args.dilate}pix") - R = args.dilate + if dilate: + LOGGER.info(f"Dilating mask using a ball of R={dilate}pix") + R = dilate r = np.arange(-R, R+1) struct = np.sqrt(r[:, np.newaxis]**2 + r[np.newaxis,:]**2) <= R mask_image = binary_dilation(input=mask_image, structure=struct) - if args.erode: - LOGGER.info(f"Applying {args.erode} iteration(s) of erosion") - N = args.erode + if erode: + LOGGER.info(f"Applying {erode} iteration(s) of erosion") + N = erode mask_image = binary_erosion(input=mask_image, iterations=N) - if args.fill_holes: + if fill_holes: LOGGER.info(f"Filling closed regions") mask_image = binary_fill_holes(mask_image) - if args.sum_peak: + if sum_peak: # This mainly to produce an image that mask out super extended sources (via sum-to-peak flux ratio) # This is useful to allow source finder to detect mainly point-like sources for cross-matching purposes only. - LOGGER.info(f"Including only flux islands with a sum-peak ratio below: {args.sum_peak}") + LOGGER.info(f"Including only flux islands with a sum-peak ratio below: {sum_peak}") extended_islands = [] mask_image_label, num_features = label(mask_image) island_objects = find_objects(mask_image_label.astype(int)) @@ -332,7 +272,7 @@ def load_fits_or_region(filename): isl_sum = (input_image[island] * mask_image[island]).sum() isl_peak = (input_image[island] * mask_image[island]).max() isl_sum_peak = isl_sum / isl_peak - if isl_sum_peak > float(args.sum_peak): + if isl_sum_peak > float(sum_peak): extended_islands.append(island) new_mask_image = np.zeros_like(mask_image) new_mask_image = new_mask_image == 0 @@ -342,36 +282,97 @@ 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(0) - if args.gui: + if outcatalog or 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 outregion: + regions.Regions(polygon_regions).write(outregion, format='ds9') + LOGGER.info(f"Saving regions in {outregion}") + + if outcatalog and restored_image: + try: + import warnings + # Suppress FittingWarnings from Astropy + # WARNING: The fit may be unsuccessful; check fit_info['message'] for more information. [astropy.modeling.fitting] + # Use context manager for handling warnings + with warnings.catch_warnings(): + warnings.resetwarnings() + warnings.filterwarnings('ignore', category=UserWarning, append=True) + from breizorro.catalog import multiprocess_contours + except ModuleNotFoundError: + msg = "Running breizorro source detector requires optional dependencies, please re-install with: pip install breizorro[all]" + LOGGER.error(msg) + raise(msg) + source_list = [] + image_data, hdu_header = get_image_data(restored_image) + fitsinfo = fitsInfo(restored_image) + mean_beam = None # Use beam info from the image header by default + if beam_size: + mean_beam = beam_size + 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(outcatalog, 'w') + catalog_out = f'# processing fits image: {restored_image} \n' + f.write(catalog_out) + 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 for cataloguing. This might take a while...') + source_list = multiprocess_contours(contours, image_data, fitsinfo, noise, 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: {outcatalog}') + + if 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) + display(input_file, mask_image, outcatalog, source_list) LOGGER.info(f"Enforcing that mask to binary") mask_image = mask_image!=0 diff --git a/breizorro/breizorro.yaml b/breizorro/breizorro.yaml new file mode 100644 index 0000000..ab9a7d4 --- /dev/null +++ b/breizorro/breizorro.yaml @@ -0,0 +1,111 @@ +cabs: + breizorro: + info: "Mask creation and manipulation for radio astronomy images (https://github.com/ratt-ru/breizorro)" + command: breizorro + inputs: + restored-image: + info: "Restored image file from which to build the mask" + dtype: File + metavar: IMAGE + abbreviation: r + mask-image: + info: "Input mask file(s). Either restored-image or mask-image must be specified." + dtype: File + abbreviation: m + metavar: MASK + threshold: + info: "Sigma threshold for masking (default = 6.5)" + dtype: float + default: 6.5 + abbreviation: t + boxsize: + info: "Box size over which to compute stats (default = 50)" + dtype: int + default: 50 + abbreviation: b + savenoise: + info: "Export noise image as FITS file" + dtype: bool + default: false + merge: + info: "Merge one or more masks or region files" + dtype: List[File] + metavar: "MASK(s)|REG(s)" + policies: + repeat: '[]' + subtract: + info: "Subtract one or more masks or region files" + dtype: List[File] + metavar: "MASK(s)|REG(s)" + policies: + repeat: '[]' + number-islands: + info: "Number the islands detected" + dtype: bool + default: false + remove-islands: + info: "Remove islands from input mask (list by number or coordinates) e.g. --remove-islands 1,18,20,20h10m13s:14d15m20s" + dtype: List[Union[int, str]] + metavar: "N|COORD" + policies: + repeat: '[]' + ignore-missing-islands: + info: "Do not throw an error if an island specified by coordinates does not exist" + dtype: bool + default: false + extract-islands: + info: "Extract islands from input mask (list by number or coordinates) e.g. --extract-islands 1,18,20,20h10m13s:14d15m20s" + dtype: List[Union[int, str]] + metavar: "N|COORD" + policies: + repeat: ' ' + minimum-size: + info: "Remove islands with areas fewer than or equal to the specified number of pixels" + dtype: int + make-binary: + info: "Replace all island numbers with 1" + dtype: bool + default: false + invert: + info: "Invert the mask" + dtype: bool + default: false + dilate: + info: "Apply dilation with a radius of R pixels" + dtype: int + metavar: "R" + default: 0 + erode: + info: "Apply N iterations of erosion" + dtype: int + metavar: "N" + default: 0 + fill-holes: + info: "Fill holes (closed regions) in the mask" + dtype: bool + default: false + sum-peak: + info: "Sum-to-peak ratio of flux islands to mask in original image" + dtype: float + ncpu: + info: "Number of processors to use for cataloging" + dtype: int + abbreviation: j + beam-size: + info: "Average beam size in arcsec if missing in the image header" + dtype: float + gui: + info: "Open mask in bokeh html gui" + dtype: bool + default: false + outputs: + outfile: + info: "Suffix for the mask image (default based on input name)" + dtype: File + abbreviation: o + outcatalog: + info: "Generate a catalog based on the region mask" + dtype: File + outregion: + info: "Generate polygon regions from the mask" + dtype: File diff --git a/breizorro/catalog.py b/breizorro/catalog.py new file mode 100644 index 0000000..0b2e784 --- /dev/null +++ b/breizorro/catalog.py @@ -0,0 +1,108 @@ +import numpy as np +from photutils import centroids +from multiprocessing import Process, Queue +from operator import itemgetter, attrgetter +from regions import PolygonSkyRegion, PolygonPixelRegion + +from breizorro.utils import get_source_size, deg2ra, deg2dec +from breizorro.utils import fitsInfo, calculate_beam_area + +def process_contour(contour, image_data, fitsinfo, noise_out): + use_max = 0 + lon = 0 + pix_size = fitsinfo['ddec'] * 3600.0 + bmaj, bmin,_ = np.array(fitsinfo['b_size']) * 3600.0 + mean_beam = 0.5 * (bmaj + bmin) + pix_beam = calculate_beam_area(bmaj/2, bmin/2, pix_size) + wcs = fitsinfo['wcs'] + while len(wcs.array_shape) > 2: + wcs = wcs.dropaxis(len(wcs.array_shape) - 1) + + contour_sky = wcs.pixel_to_world(contour[:, 1], contour[:, 0]) + polygon_region = PolygonSkyRegion(vertices=contour_sky) + pix_region = polygon_region.to_pixel(wcs) + mask = pix_region.to_mask().to_image(image_data.shape[-2:]) + # Calculate the number of pixels in the source region (where mask > 0) + source_area_pix = np.sum(mask > 0) # Count of pixels in the masked source region + # Now calculate the number of beams covering the source + source_beams = source_area_pix / pix_beam # Number of beams covering the source + try: + data = mask * image_data + nndata = data # np.flip(data, axis=0) + #nndata = nndata[~np.isnan(nndata)] + total_flux = np.sum(nndata[nndata != -0.0])/pix_beam + peak_flux = nndata.max()/pix_beam + except: + peak_flux = 0.0 + if total_flux: + total_peak_ratio = np.abs((total_flux - peak_flux) / total_flux) + # Flux density error estimation + ten_pc_error = 0.1 * total_flux # a 10% error term as an additional conservative estimate + beam_error = np.sqrt(source_beams) * noise_out + flux_density_error = np.sqrt(ten_pc_error**2 + beam_error**2) # combined error + # Calculate weighted centroid + _centroids = centroids.centroid_2dg(data) + centroid_x, centroid_y = _centroids + ra, dec = wcs.all_pix2world(centroid_x, centroid_y, 0) + # Ensure RA is positive + if ra < 0: + ra += 360 + source_flux = (round(total_flux, 5), round(flux_density_error, 5)) + source_size = get_source_size(contour, pix_size, mean_beam, image_data, total_peak_ratio, _centroids) + #source_pos = format_source_coordinates(ra, dec) + source = (ra, dec) + source_flux + source_size + catalog_out = ' '.join(str(src_prop) for src_prop in source) + else: + # Dummy source to be eliminated + lon = -np.inf + catalog_out = '' + return (ra, catalog_out, use_max) + + +def multiprocess_contours(contours, image_data, fitsinfo, noise_out, ncpu=None): + + def contour_worker(input, output): + for func, args in iter(input.get, 'STOP'): + result = func(*args) + output.put(result) + + source_list = [] + # Start worker processes + if not ncpu: + try: + import multiprocessing + ncpu = multiprocessing.cpu_count() + except: + pass + TASKS = [] + for i in range(len(contours)): + contour = contours[i] + if len(contour) > 2: + x = [] + y = [] + for j in range(len(contour)): + x.append(contour[j][0]) + y.append(contour[j][1]) + TASKS.append((process_contour,(contour, image_data, fitsinfo, noise_out))) + task_queue = Queue() + done_queue = Queue() + # Submit tasks + for task in TASKS: + task_queue.put(task) + for i in range(ncpu): + Process(target=contour_worker, args=(task_queue, done_queue)).start() + + num_max = 0 + # Get the results from parallel processing + for i in range(len(TASKS)): + catalog_out = done_queue.get(timeout=1800) + if catalog_out[0] > -np.inf: + source_list.append(catalog_out) + num_max += catalog_out[2] + # Tell child processes to stop + for i in range(ncpu): + task_queue.put('STOP') + + ra_sorted_list = sorted(source_list, key = itemgetter(0)) + + return ra_sorted_list diff --git a/breizorro/gui.py b/breizorro/gui.py new file mode 100644 index 0000000..448078a --- /dev/null +++ b/breizorro/gui.py @@ -0,0 +1,192 @@ +from breizorro.utils import fitsInfo, calculate_beam_area, format_source_coordinates + +import numpy as np +from bokeh.models import Label, BoxEditTool, FreehandDrawTool +from bokeh.models import LabelSet, DataTable, TableColumn +from bokeh.plotting import figure, ColumnDataSource, output_file, show +from bokeh.io import curdoc, export_png +from bokeh.models import Toggle, HoverTool +from bokeh.models import LabelSet, CustomJS +from bokeh.layouts import row +curdoc().theme = 'caliber' + + +def display(imagename, mask_image, outcatalog, source_list): + """ + Display the image with overlaid source positions and associated catalog information. + + Parameters + ---------- + imagename : str + The name or path of the FITS file to be displayed. + + outcatalog : bool + If True, overlay the detected sources on the image using the provided source list. + + source_list : list of tuples + A list containing information about detected sources, where each entry is a tuple + with the following format: + (RA, 'RA DEC I I_err Emaj_s Emin_s PA_d', flag). + + - RA : float + Right Ascension in degrees. + - DEC : float + Declination in degrees. + - I : float + Intensity or flux measurement. + - I_err : float + Error in the intensity measurement. + - Emaj_s : float + Major axis error (in arcseconds). + - Emin_s : float + Minor axis error (in arcseconds). + - PA_d : float + Position angle (in degrees). + - flag : int + A flag indicating some property of the source (e.g., detection confidence). + + Returns + ------- + None + Displays/Saves an image with overlaid catalog scatter points and mask image. + """ + # Get fits information + fitsinfo = fitsInfo(imagename) + # Origin coordinates + origin_ra, origin_dec = fitsinfo['centre'] + # Pixel width in degrees + pixel_width = fitsinfo['dra'] + pixel_height = fitsinfo['ddec'] + # Calculate the extent of the image in degrees + # We assume a square image for simplicity + extent_x = fitsinfo['numPix'] * pixel_width # Assume equal pixels in each dimension + extent_y = fitsinfo['numPix'] * pixel_height # Assume equal pixels in each dimension + # Ensure RA is always positive (in degrees, 0 to 360) + origin_ra = origin_ra % 360 + # Specify the coordinates for the image + x_range = (origin_ra - extent_x/2.0, origin_ra + extent_x/2.0) + y_range = (origin_dec - extent_y/2.0, origin_dec + extent_y/2.0) + if outcatalog: + # Extracting data from source_list + x_coords = [float(d[1].split(' ')[0]) for d in source_list] + y_coords = [float(d[1].split(' ')[1]) for d in source_list] + labels = [f"{format_source_coordinates(float(d[1].split(' ')[0]), float(d[1].split(' ')[1]))}" + for d in source_list] + + #Create data + source = ColumnDataSource(data=dict(x=x_coords, y=y_coords, label=labels)) + + # Assuming `source_list` is already populated with your data + # Example source_list: [(ra, 'ra dec i i_err emaj_s emin_s pa_d', flag)] + # Parse the source_list and split each string into its components + data = { + 'name': [f'src{i}' for i in range(len(source_list))], + 'ra_deg': [float(d[1].split(' ')[0]) for d in source_list], + 'dec_deg': [float(d[1].split(' ')[1]) for d in source_list], + 'i': [float(d[1].split(' ')[2]) for d in source_list], + 'error': [float(d[1].split(' ')[3]) for d in source_list], + 'emaj_s': [float(d[1].split(' ')[4]) for d in source_list], + 'emin_s': [float(d[1].split(' ')[5]) for d in source_list], + 'pa_d': [float(d[1].split(' ')[6]) for d in source_list], + } + + # Format RA and DEC to hh:mm:ss and dd:mm:ss + formatted_coords = [format_source_coordinates(ra, dec) for ra, dec in zip(data['ra_deg'], data['dec_deg'])] + formatted_RA = [coord[0] for coord in formatted_coords] + formatted_DEC = [coord[1] for coord in formatted_coords] + + # Add formatted coordinates to the data source + data['RA'] = formatted_RA + data['DEC'] = formatted_DEC + + # Create a ColumnDataSource for the table + table_source = ColumnDataSource(data=data) + from bokeh.models import NumberFormatter + # Define the number formatter for decimal and scientific notation + decimal_formatter = NumberFormatter(format="0.000") + scientific_formatter = NumberFormatter(format="0.000e") + # Define the columns for the DataTable + columns = [ + TableColumn(field="name", title="name"), + TableColumn(field="ra_deg", title="ra (deg)", formatter=decimal_formatter), + TableColumn(field="dec_deg", title="dec (deg)", formatter=decimal_formatter), + TableColumn(field="i", title="i (Jy)", formatter=scientific_formatter), + TableColumn(field="error", title="i_err (Jy)", formatter=scientific_formatter), + TableColumn(field="emaj_s", title="emaj_s (arcsec)"), + TableColumn(field="emin_s", title="emin_s (arcsec)"), + TableColumn(field="pa_d", title="pa_d (deg)"), + ] + + # Create the DataTable widget + data_table = DataTable(source=table_source, columns=columns, width=600, height=600) + + # Create the plot + p = figure(title="Breizorro Source Catalog", + x_axis_label="Right Ascension", + y_axis_label="Declination", + #x_range=(max(x_coords), min(x_coords)), + y_range=(min(y_coords), max(y_coords)), + match_aspect=True, + tooltips=[("x", "$x"), ("y", "$y"), ("value", "@image")]) + # Plot the image + image_renderer = p.image(image=[np.flip(mask_image, axis=1)], + x=x_range[0], + y=y_range[0], + dw=fitsinfo['numPix']*fitsinfo['dra'], + dh=fitsinfo['numPix']*fitsinfo['ddec'], + palette="Greys256", level="image") + scatter_renderer = p.scatter('ra_deg', 'dec_deg', + size=6, color="red", + source=table_source, + legend_label="Detection", + level="overlay") + + # Add labels to the scatter points (optional, can hide later as needed) + labels = LabelSet(x='RA', y='DEC', text='Name', source=table_source, text_font_size="10pt", text_baseline="middle", text_align="center") + + # Add hover tool for scatter points + hover = HoverTool() + hover.tooltips = [("Name", "@name"), ("RA", "@RA"), + ("DEC", "@DEC"), ("Flux", "@i")] + hover.renderers = [scatter_renderer] + p.add_tools(hover) + # Enable legend click to hide/show scatter points + p.legend.click_policy = "hide" + p.legend.label_text_color = "white" # Set the legend text color to white + p.x_range.flipped = True + p.title.align = "center" + p.add_layout(labels) + p.grid.grid_line_width = 0.5 + #from bokeh.models import DatetimeTickFormatter + #p.xaxis.formatter = DatetimeTickFormatter(hours=["%Hh%Mm%Ss"]) + #p.yaxis.formatter = DatetimeTickFormatter(minutes=["%Dd%Mm%Ss"]) + from bokeh.models import CustomJSTickFormatter + # Formatter for RA axis (converting degrees to hh:mm:ss) + p.xaxis.formatter = CustomJSTickFormatter(code=""" + // Convert RA from degrees to hours, minutes, and seconds + let hours = Math.floor(tick / 15); // 15 degrees per hour + let minutes = Math.floor((tick % 15) * 4); + let seconds = Math.floor((((tick % 15) * 4) - minutes) * 60); + return hours + "h " + minutes + "m " + seconds + "s"; + """) + + # Formatter for Dec axis (converting degrees to dd:mm:ss) + p.yaxis.formatter = CustomJSTickFormatter(code=""" + // Convert Dec from degrees to degrees, minutes, and seconds + let degrees = Math.floor(Math.abs(tick)); + let minutes = Math.floor((Math.abs(tick) - degrees) * 60); + let seconds = Math.floor((((Math.abs(tick) - degrees) * 60) - minutes) * 60); + let sign = tick < 0 ? '-' : ''; + return sign + degrees + "° " + minutes + "' " + seconds + '"'; + """) + # JavaScript callback to select the corresponding row in the table when hovering over a point + callback = CustomJS(args=dict(source=table_source), code=""" + var indices = cb_data.index['1d'].indices; // Get the index of the hovered point + source.selected.indices = indices; // Select the corresponding row in the table + """) + # Connect the hover tool to the callback + hover.callback = callback + layout = row(p, data_table, sizing_mode="stretch_both") + output_file("breizorro.html", title="Mask Editor") + #export_png(p, filename="breizorro.png") + show(layout) diff --git a/breizorro/main.py b/breizorro/main.py index 58717c1..8047838 100644 --- a/breizorro/main.py +++ b/breizorro/main.py @@ -1,4 +1,12 @@ +import os +import click +from scabha.schema_utils import clickify_parameters +from omegaconf import OmegaConf from breizorro.breizorro import main -def driver(): - main() +schemas = OmegaConf.load(os.path.join(os.path.dirname(__file__), "breizorro.yaml")) + +@click.command("breizorro") +@clickify_parameters(schemas.cabs.get("breizorro")) +def driver(**kw): + main(**kw) diff --git a/breizorro/utils.py b/breizorro/utils.py new file mode 100644 index 0000000..d23e0a8 --- /dev/null +++ b/breizorro/utils.py @@ -0,0 +1,310 @@ +import sys +import numpy +import shutil +import logging +import argparse +import os.path +import re +import numpy as np + +import astropy.units as u +from astropy.io import fits +from astropy.wcs import WCS +from astropy.coordinates import Angle +from astropy.coordinates import SkyCoord + +from regions import PixCoord +from regions import PolygonSkyRegion, PolygonPixelRegion + +def deg_to_hms(ra_deg): + ra_hours = ra_deg / 15 # 360 degrees = 24 hours + hours = int(ra_hours) + minutes = int((ra_hours - hours) * 60) + seconds = (ra_hours - hours - minutes / 60) * 3600 + return hours, minutes, seconds + + +def deg_to_dms(dec_deg): + degrees = int(dec_deg) + dec_minutes = abs(dec_deg - degrees) * 60 + minutes = int(dec_minutes) + seconds = (dec_minutes - minutes) * 60 + return degrees, minutes, seconds + + +def format_source_coordinates(ra_deg, dec_deg): + coord = SkyCoord(ra=ra_deg*u.deg, dec=dec_deg*u.deg, frame='icrs') + ra_str = coord.ra.to_string(unit=u.hour, sep=':', precision=2) + dec_str = coord.dec.to_string(unit=u.deg, sep=':', precision=2) + return ra_str, dec_str + +def format_source_coordinates0(coord_ra_deg, coord_dec_deg): + h,m,s = deg_to_hms(coord_ra_deg) + if h < 10: + h = '0' + str(h) + else: + h = str(h) + if m < 10: + m = '0' + str(m) + else: + m = str(m) + s = round(s,2) + if s < 10: + s = '0' + str(s) + else: + s = str(s) + if len(s) < 5: + s = s + '0' + h_m_s = h + ':' + m + ':' + s + + d,m,s = deg_to_dms(coord_dec_deg) + if d >= 0 and d < 10: + d = '0' + str(d) + elif d < 0 and abs(d) < 10: + d = '-0' + str(d) + else: + d = str(d) + if m < 10: + m = '0' + str(m) + else: + m = str(m) + s = round(s,2) + if s < 10: + s = '0' + str(s) + else: + s = str(s) + if len(s) < 5: + s = s + '0' + d_m_s = d + ':' + m + ':' + s + + src_pos = (h_m_s, d_m_s) + return src_pos + +def deg2dec(dec_deg, deci=2): + """Converts declination in degrees to dms coordinates + + Parameters + ---------- + dec_deg : float + Declination in degrees + dec: int + Decimal places in float format + + Returns + ------- + dms : str + Declination in degrees:arcmin:arcsec format + """ + DD = int(dec_deg) + dec_deg_abs = np.abs(dec_deg) + DD_abs = np.abs(DD) + MM = int((dec_deg_abs - DD_abs)*60) + SS = round((((dec_deg_abs - DD_abs)*60)-MM), deci) + return "%s:%s:%s"%(DD,MM,SS) + + +def deg2ra(ra_deg, deci=2): + """Converts right ascension in hms coordinates to degrees + + Parameters + ---------- + ra_deg : float + ra in degrees format + + Returns + ------- + HH:MM:SS : str + + """ + if ra_deg < 0: + ra_deg = 360 + ra_deg + HH = int((ra_deg*24)/360.) + MM = int((((ra_deg*24)/360.)-HH)*60) + SS = round(((((((ra_deg*24)/360.)-HH)*60)-MM)*60), deci) + return "%s:%s:%s"%(HH,MM,SS) + + +def calculate_area(bmaj, bmin, pix_size): + """ + Calculate the area of an ellipse represented by its major and minor axes, + given the pixel size. + + Parameters: + bmaj (float): Major axis raduis of the ellipse in arcseconds. + bmin (float): Minor axis radius of the ellipse in arcseconds. + pix_size (float): Pixel size in arcseconds. + + Returns: + float: Calculated area of the ellipse in square pixels. + """ + # Calculate the semi-major and semi-minor axes in pixels + a_pixels = bmaj / pix_size + b_pixels = bmin / pix_size + + # Calculate the area of the ellipse using the formula: π * a * b + area = np.pi * a_pixels * b_pixels + + return area + + +def fitsInfo(fitsname=None): + """Get fits header info. + + Parameters + ---------- + fitsname : fits file + Restored image (cube) + + Returns + ------- + fitsinfo : dict + Dictionary of fits information + e.g. {'wcs': wcs, 'ra': ra, 'dec': dec, + 'dra': dra, 'ddec': ddec, 'raPix': raPix, + 'decPix': decPix, 'b_size': beam_size, + 'numPix': numPix, 'centre': centre, + 'skyArea': skyArea, 'naxis': naxis} + + """ + with fits.open(fitsname) as hdu: + hdr = hdu[0].header + ra = hdr['CRVAL1'] + dra = abs(hdr['CDELT1']) + raPix = hdr['CRPIX1'] + dec = hdr['CRVAL2'] + ddec = abs(hdr['CDELT2']) + decPix = hdr['CRPIX2'] + wcs = WCS(hdr) + numPix = hdr['NAXIS1'] + naxis = hdr['NAXIS'] + try: + beam_size = (hdr['BMAJ'], hdr['BMIN'], hdr['BPA']) + except: + beam_size = None + try: + centre = (hdr['CRVAL1'], hdr['CRVAL2']) + except: + centre = None + try: + freq0=None + for i in range(1, hdr['NAXIS']+1): + if hdr['CTYPE{0:d}'.format(i)].startswith('FREQ'): + freq0 = hdr['CRVAL{0:d}'.format(i)] + except: + freq0=None + + skyArea = (numPix * ddec) ** 2 + fitsinfo = {'wcs': wcs, 'ra': ra, 'dec': dec, 'naxis': naxis, + 'dra': dra, 'ddec': ddec, 'raPix': raPix, + 'decPix': decPix, 'b_size': beam_size, + 'numPix': numPix, 'centre': centre, + 'skyArea': skyArea, 'freq0': freq0} + return fitsinfo + + +def get_image_data(fitsfile): + """ + Reads a FITS file and returns a tuple of the image array and the header. + + Parameters: + fitsfile (str): The path to the FITS file. + + Returns: + tuple: A tuple containing the image array and the header. + """ + with fits.open(fitsfile) as input_hdu: + # Check the dimensionality of the data and extract the appropriate slice + if len(input_hdu[0].data.shape) == 2: + image = np.array(input_hdu[0].data[:, :]) + elif len(input_hdu[0].data.shape) == 3: + image = np.array(input_hdu[0].data[0, :, :]) + else: + image = np.array(input_hdu[0].data[0, 0, :, :]) + + header = input_hdu[0].header + + return image, header + + +def maxDist(contour, pixel_size, x_centroid, y_centroid): + """Calculate maximum extent and position angle of a contour. + + Parameters: + contour : list of [x, y] pairs + List of coordinates defining the contour. + pixel_size : float + Size of a pixel in the image (e.g., arcseconds per pixel). + x_centroid: int + X-coordinate of the centroid. + y_centroid: int + Y-coordinate of the centroid. + + Returns: + e_maj : float + Major axis of the contour in angular units (e.g., arcseconds). + e_min : float + Minor axis of the contour in angular units (e.g., arcseconds). + pos_angle : float + Position angle of the contour (in degrees). + """ + contour_array = np.array(contour) + distances = np.linalg.norm(contour_array - [y_centroid, x_centroid], axis=1) + + # Find the indices of the points with maximum and minimum distances + max_idx = np.argmax(distances) + min_idx = np.argmin(distances) + + # Calculate major and minor axes + e_maj = np.max(distances) * pixel_size + e_min = np.min(distances) * pixel_size + + # Calculate position angle + dx, dy = contour_array[max_idx] - [x_centroid, y_centroid] + pos_angle = np.degrees(np.arctan2(dy, dx)) + + return e_maj, e_min, pos_angle + +def calculate_beam_area(bmaj, bmin, pix_size): + """ + Calculate the area of an ellipse represented by its major and minor axes, + given the pixel size. + + Parameters: + bmaj (float): Major axis of the ellipse in arcseconds. + bmin (float): Minor axis of the ellipse in arcseconds. + pix_size (float): Pixel size in arcseconds. + + Returns: + area (float): Calculated area of the ellipse in square pixels. + """ + # Calculate the semi-major and semi-minor axes in pixels + a_pixels = bmaj / pix_size + b_pixels = bmin / pix_size + + # Calculate the area of the ellipse using the formula: π * a * b + area = np.pi * a_pixels * b_pixels + + return area + + +def get_source_size(contour, pixel_size, mean_beam, image, int_peak_ratio, centroids): + result = maxDist(contour,pixel_size,*centroids) + src_angle = result[:-1] + pos_angle = result[-1] + contour_pixels = PixCoord([c[0] for c in contour], [c[1] for c in contour]) + p = PolygonPixelRegion(vertices=contour_pixels, meta={'label': 'Region'}) + source_beam_ratio = p.area / mean_beam + # first test for point source + point_source = False + if (int_peak_ratio <= 0.2) or (src_angle[0] <= mean_beam): + point_source = True + if source_beam_ratio <= 1.0: + point_source = True + if point_source: + src_size = (0.0, 0.0, 0.0) + else: + emaj = round(src_angle[0],2) + emin = round(src_angle[-1],2) + pa = round(pos_angle,2) + src_size = (emaj, emin, pa) + return src_size diff --git a/index.md b/index.md new file mode 100644 index 0000000..8ba624e --- /dev/null +++ b/index.md @@ -0,0 +1,127 @@ +

+ Breizorro Logo +

+ +Breizorro is a flexible software program made to simplify image analysis tasks, including identifying emission islands and generating and modifying image masks, which are frequently used in radio interferometric imaging. + + +# Parameter definition + +``` +Usage: breizorro [OPTIONS] + +Options: + --restored-image PATH Restored image file from which to build the + mask + --mask-image PATH Input mask file(s). Either restored-image or + mask-image must be specified. + --threshold FLOAT Sigma threshold for masking (default = 6.5) + --boxsize INTEGER Box size over which to compute stats + (default = 50) + --savenoise / --no-savenoise Export noise image as FITS file + --merge TEXT Merge one or more masks or region files + --subtract TEXT Subtract one or more masks or region files + --number-islands / --no-number-islands + Number the islands detected + --remove-islands TEXT Remove islands from input mask (list by + number or coordinates) + --ignore-missing-islands / --no-ignore-missing-islands + Do not throw an error if an island specified + by coordinates does not exist + --extract-islands TEXT Extract islands from input mask (list by + number or coordinates) + --minimum-size INTEGER Remove islands with areas fewer than or + equal to the specified number of pixels + --make-binary / --no-make-binary + Replace all island numbers with 1 + --invert / --no-invert Invert the mask + --dilate INTEGER Apply dilation with a radius of R pixels + --erode INTEGER Apply N iterations of erosion + --fill-holes / --no-fill-holes Fill holes (closed regions) in the mask + --sum-peak FLOAT Sum-to-peak ratio of flux islands to mask in + original image + --ncpu INTEGER Number of processors to use for cataloging + --beam-size FLOAT Average beam size in arcsec if missing in + the image header + --gui / --no-gui Open mask in bokeh html gui + --outfile PATH Suffix for the mask image (default based on + input name) + --outcatalog PATH Generate a catalog based on the region mask + --outregion PATH Generate polygon regions from the mask + --help Show this message and exit. +``` + + +# Image Masking and Transformation + +Breizorro uses the minimal filter to generate binary masks by replacing each pixel's value with the lowest value found in its near vicinity. This is incredibly effective at reducing noise and highlighting or smoothing particular regions of an image, such as the regions surrounding bright or small sources. Users can specify a window (or kernel) of a specific size that moves over the image as well as a sigma threshold for masking. + +``` +breizorro -t 6.5 -b 50 -r circinus-MFS-image.fits +``` + +![mypipelinerun_circinus_p3_3-MFS-image fits-mypipelinerun_circinus_p3_3-MFS-image mask fits-image-2024-09-11-10-35-43](https://github.com/user-attachments/assets/cfd2f918-340a-4148-96a2-c00ca41b33d0) + + +``` +breizorro -r circinus-MFS-image.fits --sum-peak 500 +``` + +![mypipelinerun_circinus_p3_3-MFS-image mask fits-mypipelinerun_circinus_p3_3-MFS-image mask fits-image-2024-09-11-13-23-14](https://github.com/user-attachments/assets/0ff50068-ec8a-42bf-8539-9f68f15a1ea9) + +# Region Generation and Manipulation + +Breizorro makes it easier to create and work with regions using image masks. It includes labelling, eliminating, extracting, and filtering regions (islands) based on user-specified criteria. Users can refine their regions of interest using techniques such as erosion, dilation, hole-filling, binary masking, and inversion. + +``` +breizorro -r circinus-MFS-image.fits --outregion circinus.reg +``` + +![mypipelinerun_circinus_p3_3-MFS-image fits-image-2024-09-11-10-38-15](https://github.com/user-attachments/assets/14f435e1-6234-4515-9597-c3002a644975) + +``` +breizorro -r circinus-MFS-image.fits --merge west.reg --dilate 1 --fill-holes +``` + +![mypipelinerun_circinus_p3_3-MFS-image mask fits-mypipelinerun_circinus_p3_3-MFS-image mask fits-image-2024-09-11-13-59-39](https://github.com/user-attachments/assets/2308c7b7-2ec0-4895-b93b-5d96f3d99337) + + +# Cataloguing and Visualization + +Breizorro enables catalogue generation from extracted regions, saving source properties to an ASCII/text file. +This is particularly useful for analyzing fields dominated by point sources. +By efficiently parameterizing and cataloguing compact sources, Breizorro enables rapid cross-matching. + +``` +breizorro -r deep2-MFS-image.fits --outcatalog deep2.txt +``` + +``` +# processing fits image: deep2-MFS-image.fits +# mean beam size (arcsec): 6.31 +# original image peak flux (Jy/beam): 0.023107271641492844 +# noise out (µJy/beam): 51.62 +# cutt-off flux (mJy/beam): 0.52 +# freq0 (Hz): 1049833007.8125 +# number of sources detected: 100 +# +#format: name ra_d dec_d i i_err emaj_s emin_s pa_d +src0 60.64975237181635 -79.86348735299585 0.00369 0.0001 11.66 0.0 120.96 +src1 60.67140679629887 -80.36126947232378 7e-05 0.0001 0.0 0.0 0.0 +src2 60.877334392876136 -79.76691180042988 0.00113 0.0001 8.25 0.0 104.04 +src3 60.894387589282964 -79.73476060235502 0.00023 0.0001 0.0 0.0 0.0 +src4 60.95196895741357 -79.68021884337942 0.00054 0.0001 0.0 0.0 0.0 +src5 61.00657438220518 -80.17745581694626 0.00018 0.0001 0.0 0.0 0.0 +src6 61.04136845645677 -80.48097816446368 0.00032 0.0001 0.0 0.0 0.0 +src7 61.061502553348895 -79.93907167920088 0.01435 0.0001 21.54 0.0 111.8 +src8 61.18650068905602 -80.40952881341589 0.00027 0.0001 0.0 0.0 0.0 +src9 61.32594143681846 -79.91488530476678 7e-05 0.0001 0.0 0.0 0.0 +src10 61.47462421089555 -79.96912396162651 9e-05 0.0001 0.0 0.0 0.0 +src11 61.56419377531346 -80.18249455745902 0.0022 0.0001 10.77 0.0 158.2 +``` + +# Contributors + +Thank you to the people who have contributed to this project. + +[![Contributors](https://contrib.rocks/image?repo=ratt-ru/breizorro)](https://github.com/ratt-ru/breizorro/graphs/contributors) diff --git a/pyproject.toml b/pyproject.toml index 9ee5e6d..6ada5d9 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [tool.poetry] name = "breizorro" -version = "0.1.2" +version = "0.1.3" description = "Creates a binary mask given a FITS image" authors = ["Ian Heywood & RATT "] readme = "README.rst" @@ -24,22 +24,36 @@ packages = [ breizorro = "breizorro.main:driver" [tool.poetry.dependencies] -python = "^3.8" -astropy = "*" -numpy = "*" +python = ">=3.8, <3.13" +astropy = [ + { version = ">=6.1.0", python = ">=3.12"}, + { version = "<6.1.0", python = ">=3.8, <3.12"} +] +numpy = [ + { version = ">=2.0.0", python = ">=3.12"}, + { version = "<2.0.0", python = ">=3.8, <3.12"} +] +omegaconf = "*" regions = "*" -scipy = "*" +scipy = [ + { version = ">=1.13.0", python = ">=3.12"}, + { version = "<1.13.0", python = ">=3.8, <3.12"} +] +stimela = ">=2.0" # Optional dependencies start here bokeh = { version = "*", optional = true} +photutils = { version = "*", optional = true} pytest = { version = "*", optional=true } pytest-flake8 = { version = "*", optional=true } +scikit-image = { version = "*", optional = true} [build-system] requires = ["setuptools", "poetry-core"] build-backend = "poetry.core.masonry.api" [tool.poetry.extras] +all = ["bokeh", "photutils", "scikit-image"] gui = ["bokeh"] testing = ["pytest", "pytest-flake8"]