From b0f80483e3dc40f18c261edb9b38135ae29607ee Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Thu, 14 Mar 2024 16:07:35 +0200 Subject: [PATCH 01/48] Add cataloging routiine --- breizorro/breizorro.py | 455 ++++++++++++++++++++++++++++++++++++++++- 1 file changed, 447 insertions(+), 8 deletions(-) diff --git a/breizorro/breizorro.py b/breizorro/breizorro.py index 3c32455..5148aa9 100644 --- a/breizorro/breizorro.py +++ b/breizorro/breizorro.py @@ -6,10 +6,16 @@ import os.path import re import numpy as np + 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 +23,12 @@ import scipy.special import scipy.ndimage +from skimage.draw import polygon as skimage_polygon +from skimage.measure import find_contours + +from multiprocessing import Process, Queue + +from operator import itemgetter, attrgetter def create_logger(): """Create a console logger""" @@ -86,6 +98,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 +117,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'): @@ -117,6 +132,362 @@ def remove_regions(mask_image, regs, wcs): mask_image[reg.to_mask().to_image(mask_image.shape) != 0] = 0 +def contour_worker(input, output): + for func, args in iter(input.get, 'STOP'): + result = func(*args) + output.put(result) + + +def calculate_weighted_centroid(x, y, flux_values): + # Calculate the total flux within the region + total_flux = np.sum(flux_values) + # Initialize variables for weighted sums + weighted_sum_x = 0 + weighted_sum_y = 0 + # Loop through all pixels within the region + for xi, yi, flux in zip(x, y, flux_values): + # Add the weighted contribution of each pixel to the centroid + weighted_sum_x += xi * flux + weighted_sum_y += yi * flux + # Calculate the centroid coordinates + centroid_x = weighted_sum_x / total_flux + centroid_y = weighted_sum_y / total_flux + return centroid_x, centroid_y + + +def process_contour(x, y, image_data, fitsinfo, flux_cutoff, noise_out, beam): + rr, cc = skimage_polygon(x,y) + data_result = image_data[rr, cc] + pixel_size = fitsinfo['ddec'] * 3600.0 + if fitsinfo['b_size']: + bmaj,bmin,_ = np.array(fitsinfo['b_size']) * 3600.0 + mean_beam = 0.5 * (bmaj + bmin) + else: + mean_beam = beam + pixels_beam = calculate_area(bmaj, bmin, pixel_size) + total_flux = data_result.sum() / pixels_beam + contained_points = len(rr) # needed for estimating flux density error + use_max = 0 + lon = 0 + wcs = fitsinfo['wcs'] + while len(wcs.array_shape) > 2: + wcs = wcs.dropaxis(len(wcs.array_shape) - 1) + try: + peak_flux = data_result.max() + except: + LOGGER.warn('Failure to get maximum within contour') + LOGGER.info('Probably a contour at edge of image - skipping') + peak_flux = 0.0 + if peak_flux >= flux_cutoff: + total_peak_ratio = np.abs((total_flux - peak_flux) / total_flux) + beam_error = contained_points/pixels_beam * noise_out + ten_pc_error = 0.1 * total_flux + flux_density_error = np.sqrt(ten_pc_error * ten_pc_error + beam_error * beam_error) + contour = [] + for i in range(len(x)): + contour.append([x[i],y[i]]) + centroid = calculate_weighted_centroid(x,y, data_result) + pix_centroid = PixCoord(centroid[0], centroid[1]) + contour_pixels = PixCoord(np.array(x), np.array(y)) + p = PolygonPixelRegion(vertices=contour_pixels, meta={'label': 'Region'}) + if p.contains(pix_centroid)[0]: + ra, dec = wcs.all_pix2world(pix_centroid.x, pix_centroid.y,0) + else: + use_max = 1 + LOGGER.warn('centroid lies outside polygon - using peak position') + location = np.unravel_index(np.argmax(data_result, axis=None), data_result.shape) + x_pos = rr[location] + y_pos = cc[location] + data_max = image_data[x_pos,y_pos] + pos_pixels = PixCoord(x_pos, y_pos) + ra, dec = wcs.all_pix2world(pos_pixels.x, pos_pixels.y, 0) + + source_flux = (round(total_flux * 1000, 3), round(flux_density_error * 1000, 4)) + source_size = get_source_size(contour, pixel_size, mean_beam, image_data, total_peak_ratio) + source_pos = format_source_coordinates(ra, dec) + source = source_pos + (ra, dec) + (total_flux, flux_density_error) + 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 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(coord_x_deg, coord_y_deg): + h,m,s = deg_to_hms(coord_x_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' + + d,m,s = deg_to_dms(coord_y_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' + + h_m_s = h + ':' + m + ':' + s + d_m_s = d + ':' + m + ':' + s + src_pos = (h_m_s, d_m_s) + return src_pos + +def maxDist(contour, pixel_size): + """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). + + Returns: + ang_size : float + Maximum extent of the contour in angular units (e.g., arcseconds). + pos_angle : float + Position angle of the contour (in degrees). + """ + src_size = 0 + pos_angle = None + + # Convert the contour to a numpy array for easier calculations + contour_array = np.array(contour) + + # Calculate pairwise distances between all points in the contour + for i in range(len(contour_array)): + for j in range(i+1, len(contour_array)): + # Calculate Euclidean distance between points i and j + distance = np.linalg.norm(contour_array[i] - contour_array[j]) * pixel_size + + # Calculate positional angle between points i and j + dx, dy = contour_array[j] - contour_array[i] + angle = np.degrees(np.arctan2(dy, dx)) + + # Update max_distance, max_points, and pos_angle if the calculated distance is greater + if distance > src_size: + src_size = distance + pos_angle = angle + + return src_size, pos_angle + + +def get_source_size(contour, pixel_size, mean_beam, image, int_peak_ratio): + result = maxDist(contour,pixel_size) + src_angle = result[0] + 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 <= mean_beam): + point_source = True + if source_beam_ratio <= 1.0: + point_source = True + if point_source: + src_size = (0.0, 0.0) + print(f"Point source because {int_peak_ratio} <= 0.2 and {src_angle} <= {mean_beam}") + else: + ang = round(src_angle,2) + pa = round(pos_angle,2) + src_size = (ang, pa) + return src_size + + +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} + + """ + hdu = fits.open(fitsname) + 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 calculate_area(b_major, b_minor, pixel_size): + """ + Calculate the area of an ellipse represented by its major and minor axes, + given the pixel size. + + Parameters: + b_major (float): Major axis of the ellipse in arcseconds. + b_minor (float): Minor axis of the ellipse in arcseconds. + pixel_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 = b_major / pixel_size + b_pixels = b_minor / pixel_size + + # Calculate the area of the ellipse using the formula: π * a * b + area = np.pi * a_pixels * b_pixels + + return area + + +def generate_source_list(filename, outfile, mask_image, threshold_value, noise, num_processors, mean_beam): + image_data, hdu_header = get_image(filename) + fitsinfo = fitsInfo(filename) + f = open(outfile, 'w') + catalog_out = f'# processing fits image: {filename} \n' + f.write(catalog_out) + incoming_dimensions = fitsinfo['naxis'] + pixel_size = fitsinfo['ddec'] * 3600.0 + 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 5') + 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_value + catalog_out = f'# limiting_flux (mJy/beam): {round(limiting_flux*1000,2)} \n' + f.write(catalog_out) + + contours = find_contours(mask_image, 0.5) + # Start worker processes + if not num_processors: + try: + import multiprocessing + num_processors = multiprocessing.cpu_count() + LOGGER.info(f'Setting number of processors to {num_processors}') + 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,(x,y, image_data, fitsinfo, limiting_flux, noise, mean_beam))) + task_queue = Queue() + done_queue = Queue() + # Submit tasks + LOGGER.info('Submitting distributed tasks. This might take a while...') + for task in TASKS: + task_queue.put(task) + for i in range(num_processors): + Process(target=contour_worker, args=(task_queue, done_queue)).start() + + source_list = [] + num_max = 0 + # Get the results from parallel processing + for i in range(len(TASKS)): + catalog_out = done_queue.get(timeout=300) + 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(num_processors): + task_queue.put('STOP') + + ra_sorted_list = sorted(source_list, key = itemgetter(0)) + LOGGER.info(f'Number of sources detected {len(ra_sorted_list)}') + catalog_out = f'# number of sources detected {len(ra_sorted_list)} \n' + f.write(catalog_out) + catalog_out = f'# number of souces using peak for source position: {num_max} \n' + f.write(catalog_out) + catalog_out = '#\n# source ra_hms dec_dms ra(deg) dec(deg) flux(mJy) error ang_size_(arcsec) pos_angle_(deg)\n' + f.write(catalog_out) + for i in range(len(ra_sorted_list)): + output = str(i) + ', ' + ra_sorted_list[i][1] + '\n' + f.write(output) + f.close() + return source_list + def main(): LOGGER.info("Welcome to breizorro") # Get version @@ -167,9 +538,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, + 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.') @@ -343,23 +723,81 @@ 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.outcatalog: + num_proc = None # Assign based on the number of processors on a system + if args.ncpu: + num_proc = int(args.ncpu) + mean_beam = None # Use beam info from the image header by default + if args.beam: + mean_beam = float(args.beam) + noise = np.median(noise_image) + source_list = generate_source_list(args.imagename, args.outcatalog, mask_image, + threshold, noise, num_proc, mean_beam) + LOGGER.info(f'Source catalog saved: {args.outcatalog}') + if args.gui: try: - from bokeh.models import BoxEditTool, ColumnDataSource, FreehandDrawTool + from bokeh.models import Label, BoxEditTool, ColumnDataSource, FreehandDrawTool from bokeh.plotting import figure, output_file, show + from bokeh.io import curdoc, export_png from bokeh.io import curdoc curdoc().theme = 'caliber' 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 + fitsinfo = fitsInfo(args.imagename or args.maskname) + + # Origin coordinates + origin_ra, origin_dec = fitsinfo['centre'] + + # Pixel width in degrees + pixel_width = fitsinfo['ddec'] + + # Calculate the extent of the image in degrees + # We assume a square image for simplicity + extent = fitsinfo['numPix'] * pixel_width # Assume equal pixels in each dimension + + # Specify the coordinates for the image + xy_origin = (origin_ra - extent/2.0, origin_dec - extent/2.0) + + # Define the plot p = figure(tooltips=[("x", "$x"), ("y", "$y"), ("value", "@image")]) - p.x_range.range_padding = p.y_range.range_padding = 0 + p.x_range.range_padding = 0 + p.x_range.flipped = True 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") + # Plot the image using degree coordinates + p.image(image=[np.flip(mask_image, axis=1)], x=xy_origin[0], y=xy_origin[1], dw=extent, dh=extent, palette="Greys256", level="image") + + + # Extracting data from source_list + x_coords = [float(d[1].split(', ')[2]) for d in source_list] + y_coords = [float(d[1].split(', ')[3]) for d in source_list] + labels = [f"({d[1].split(', ')[0]}, {d[1].split(', ')[1]})" for d in source_list] + # Plotting points + p.circle(x_coords, y_coords, size=1, color="red", alpha=0.5) + #for i, (x, y, label) in enumerate(zip(x_coords, y_coords, labels)): + # p.add_layout(Label(x=x, y=y, text=label, text_baseline="middle", text_align="center", text_font_size="10pt")) p.grid.grid_line_width = 0.5 src1 = ColumnDataSource({'x': [], 'y': [], 'width': [], 'height': [], 'alpha': []}) src2 = ColumnDataSource({'xs': [], 'ys': [], 'alpha': []}) @@ -371,11 +809,12 @@ def load_fits_or_region(filename): p.add_tools(draw_tool2) p.toolbar.active_drag = draw_tool1 output_file("breizorro.html", title="Mask Editor") + export_png(p, filename="breizorro.png") show(p) - LOGGER.info(f"Enforcing that mask to binary") - mask_image = mask_image!=0 - mask_header['BUNIT'] = 'mask' + 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) From b3998b8026ed6a220b72960d82bb2a1d0cb0db1f Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Sun, 17 Mar 2024 21:18:11 +0200 Subject: [PATCH 02/48] updates --- breizorro/breizorro.py | 39 +++++++++++++++++++++++---------------- 1 file changed, 23 insertions(+), 16 deletions(-) diff --git a/breizorro/breizorro.py b/breizorro/breizorro.py index 5148aa9..33395b6 100644 --- a/breizorro/breizorro.py +++ b/breizorro/breizorro.py @@ -166,6 +166,7 @@ def process_contour(x, y, image_data, fitsinfo, flux_cutoff, noise_out, beam): mean_beam = beam pixels_beam = calculate_area(bmaj, bmin, pixel_size) total_flux = data_result.sum() / pixels_beam + print(total_flux) contained_points = len(rr) # needed for estimating flux density error use_max = 0 lon = 0 @@ -186,21 +187,26 @@ def process_contour(x, y, image_data, fitsinfo, flux_cutoff, noise_out, beam): contour = [] for i in range(len(x)): contour.append([x[i],y[i]]) - centroid = calculate_weighted_centroid(x,y, data_result) - pix_centroid = PixCoord(centroid[0], centroid[1]) - contour_pixels = PixCoord(np.array(x), np.array(y)) - p = PolygonPixelRegion(vertices=contour_pixels, meta={'label': 'Region'}) - if p.contains(pix_centroid)[0]: - ra, dec = wcs.all_pix2world(pix_centroid.x, pix_centroid.y,0) - else: + #centroid = calculate_weighted_centroid(x,y, data_result) + #pix_centroid = PixCoord(centroid[0], centroid[1]) + #contour_pixels = PixCoord(np.array(x), np.array(y)) + #p = PolygonPixelRegion(vertices=contour_pixels, meta={'label': 'Region'}) + #ra, dec = wcs.all_pix2world(pix_centroid.x, pix_centroid.y,0) + #if p.contains(pix_centroid)[0]: + # ra, dec = wcs.all_pix2world(pix_centroid.x, pix_centroid.y,0) + #else: + if True: use_max = 1 LOGGER.warn('centroid lies outside polygon - using peak position') location = np.unravel_index(np.argmax(data_result, axis=None), data_result.shape) x_pos = rr[location] y_pos = cc[location] + data_max = image_data[x_pos,y_pos] / pixels_beam data_max = image_data[x_pos,y_pos] pos_pixels = PixCoord(x_pos, y_pos) ra, dec = wcs.all_pix2world(pos_pixels.x, pos_pixels.y, 0) + print(data_max) + print(f'{ra},{dec}') source_flux = (round(total_flux * 1000, 3), round(flux_density_error * 1000, 4)) source_size = get_source_size(contour, pixel_size, mean_beam, image_data, total_peak_ratio) @@ -244,6 +250,7 @@ def format_source_coordinates(coord_x_deg, coord_y_deg): s = str(s) if len(s) < 5: s = s + '0' + h_m_s = h + ':' + m + ':' + s d,m,s = deg_to_dms(coord_y_deg) if d >= 0 and d < 10: @@ -263,9 +270,8 @@ def format_source_coordinates(coord_x_deg, coord_y_deg): s = str(s) if len(s) < 5: s = s + '0' - - h_m_s = h + ':' + m + ':' + s d_m_s = d + ':' + m + ':' + s + src_pos = (h_m_s, d_m_s) return src_pos @@ -426,7 +432,7 @@ def generate_source_list(filename, outfile, mask_image, threshold_value, noise, raise('No beam information found. Specify mean beam in arcsec: --beam-size 5') 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' + 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) @@ -777,17 +783,18 @@ def load_fits_or_region(filename): extent = fitsinfo['numPix'] * pixel_width # Assume equal pixels in each dimension # Specify the coordinates for the image - xy_origin = (origin_ra - extent/2.0, origin_dec - extent/2.0) + x_range = (origin_ra - extent/2.0, origin_ra + extent/2.0) + y_range = (origin_dec - extent/2.0, origin_dec + extent/2.0) # Define the plot - p = figure(tooltips=[("x", "$x"), ("y", "$y"), ("value", "@image")]) + p = figure(tooltips=[("x", "$x"), ("y", "$y"), ("value", "@image")], y_range=y_range) p.x_range.range_padding = 0 p.x_range.flipped = True p.title.text = out_mask_fits # must give a vector of image data for image parameter # Plot the image using degree coordinates - p.image(image=[np.flip(mask_image, axis=1)], x=xy_origin[0], y=xy_origin[1], dw=extent, dh=extent, palette="Greys256", level="image") + p.image(image=[np.flip(mask_image, axis=1)], x=x_range[0], y=y_range[0], dw=extent, dh=extent, palette="Greys256", level="image") # Extracting data from source_list @@ -795,9 +802,9 @@ def load_fits_or_region(filename): y_coords = [float(d[1].split(', ')[3]) for d in source_list] labels = [f"({d[1].split(', ')[0]}, {d[1].split(', ')[1]})" for d in source_list] # Plotting points - p.circle(x_coords, y_coords, size=1, color="red", alpha=0.5) - #for i, (x, y, label) in enumerate(zip(x_coords, y_coords, labels)): - # p.add_layout(Label(x=x, y=y, text=label, text_baseline="middle", text_align="center", text_font_size="10pt")) + p.circle(x_coords, y_coords, size=3, color="red", alpha=0.5) + for i, (x, y, label) in enumerate(zip(x_coords, y_coords, labels)): + p.add_layout(Label(x=x, y=y, text=label, text_baseline="middle", text_align="center", text_font_size="10pt")) p.grid.grid_line_width = 0.5 src1 = ColumnDataSource({'x': [], 'y': [], 'width': [], 'height': [], 'alpha': []}) src2 = ColumnDataSource({'xs': [], 'ys': [], 'alpha': []}) From f904e458fd5554040fe36341ac3f66b07fa24e92 Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Wed, 20 Mar 2024 11:53:48 +0200 Subject: [PATCH 03/48] Restructuring and clean-ups --- breizorro/__init__.py | 0 breizorro/breizorro.py | 348 ++++++++--------------------------------- breizorro/utils.py | 250 +++++++++++++++++++++++++++++ 3 files changed, 319 insertions(+), 279 deletions(-) create mode 100644 breizorro/__init__.py create mode 100644 breizorro/utils.py 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 33395b6..c7670ca 100644 --- a/breizorro/breizorro.py +++ b/breizorro/breizorro.py @@ -30,6 +30,10 @@ from operator import itemgetter, attrgetter +from breizorro.utils import get_source_size, format_source_coordinates +from breizorro.utils import fitsInfo + + def create_logger(): """Create a console logger""" log = logging.getLogger(__name__) @@ -132,30 +136,8 @@ def remove_regions(mask_image, regs, wcs): mask_image[reg.to_mask().to_image(mask_image.shape) != 0] = 0 -def contour_worker(input, output): - for func, args in iter(input.get, 'STOP'): - result = func(*args) - output.put(result) - - -def calculate_weighted_centroid(x, y, flux_values): - # Calculate the total flux within the region - total_flux = np.sum(flux_values) - # Initialize variables for weighted sums - weighted_sum_x = 0 - weighted_sum_y = 0 - # Loop through all pixels within the region - for xi, yi, flux in zip(x, y, flux_values): - # Add the weighted contribution of each pixel to the centroid - weighted_sum_x += xi * flux - weighted_sum_y += yi * flux - # Calculate the centroid coordinates - centroid_x = weighted_sum_x / total_flux - centroid_y = weighted_sum_y / total_flux - return centroid_x, centroid_y - +def process_contour(contour, image_data, fitsinfo, flux_cutoff, noise_out, beam): -def process_contour(x, y, image_data, fitsinfo, flux_cutoff, noise_out, beam): rr, cc = skimage_polygon(x,y) data_result = image_data[rr, cc] pixel_size = fitsinfo['ddec'] * 3600.0 @@ -219,233 +201,20 @@ def process_contour(x, y, image_data, fitsinfo, flux_cutoff, noise_out, beam): catalog_out = '' return (ra, catalog_out, use_max) -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(coord_x_deg, coord_y_deg): - h,m,s = deg_to_hms(coord_x_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_y_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 maxDist(contour, pixel_size): - """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). - - Returns: - ang_size : float - Maximum extent of the contour in angular units (e.g., arcseconds). - pos_angle : float - Position angle of the contour (in degrees). - """ - src_size = 0 - pos_angle = None - - # Convert the contour to a numpy array for easier calculations - contour_array = np.array(contour) - - # Calculate pairwise distances between all points in the contour - for i in range(len(contour_array)): - for j in range(i+1, len(contour_array)): - # Calculate Euclidean distance between points i and j - distance = np.linalg.norm(contour_array[i] - contour_array[j]) * pixel_size - - # Calculate positional angle between points i and j - dx, dy = contour_array[j] - contour_array[i] - angle = np.degrees(np.arctan2(dy, dx)) - - # Update max_distance, max_points, and pos_angle if the calculated distance is greater - if distance > src_size: - src_size = distance - pos_angle = angle - - return src_size, pos_angle - - -def get_source_size(contour, pixel_size, mean_beam, image, int_peak_ratio): - result = maxDist(contour,pixel_size) - src_angle = result[0] - 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 <= mean_beam): - point_source = True - if source_beam_ratio <= 1.0: - point_source = True - if point_source: - src_size = (0.0, 0.0) - print(f"Point source because {int_peak_ratio} <= 0.2 and {src_angle} <= {mean_beam}") - else: - ang = round(src_angle,2) - pa = round(pos_angle,2) - src_size = (ang, pa) - return src_size - - -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} - - """ - hdu = fits.open(fitsname) - 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 calculate_area(b_major, b_minor, pixel_size): - """ - Calculate the area of an ellipse represented by its major and minor axes, - given the pixel size. - - Parameters: - b_major (float): Major axis of the ellipse in arcseconds. - b_minor (float): Minor axis of the ellipse in arcseconds. - pixel_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 = b_major / pixel_size - b_pixels = b_minor / pixel_size - - # Calculate the area of the ellipse using the formula: π * a * b - area = np.pi * a_pixels * b_pixels - - return area +def multiprocess_contour(contours, ncpu=None): + def contour_worker(input, output): + for func, args in iter(input.get, 'STOP'): + result = func(*args) + output.put(result) -def generate_source_list(filename, outfile, mask_image, threshold_value, noise, num_processors, mean_beam): - image_data, hdu_header = get_image(filename) - fitsinfo = fitsInfo(filename) - f = open(outfile, 'w') - catalog_out = f'# processing fits image: {filename} \n' - f.write(catalog_out) - incoming_dimensions = fitsinfo['naxis'] - pixel_size = fitsinfo['ddec'] * 3600.0 - 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 5') - 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_value - catalog_out = f'# limiting_flux (mJy/beam): {round(limiting_flux*1000,2)} \n' - f.write(catalog_out) - - contours = find_contours(mask_image, 0.5) + source_list = [] # Start worker processes - if not num_processors: + if not ncpu: try: import multiprocessing - num_processors = multiprocessing.cpu_count() + ncpu = multiprocessing.cpu_count() LOGGER.info(f'Setting number of processors to {num_processors}') except: pass @@ -458,17 +227,16 @@ def generate_source_list(filename, outfile, mask_image, threshold_value, noise, for j in range(len(contour)): x.append(contour[j][0]) y.append(contour[j][1]) - TASKS.append((process_contour,(x,y, image_data, fitsinfo, limiting_flux, noise, mean_beam))) + TASKS.append((process_contour,(contour, image_data, fitsinfo, limiting_flux, noise, mean_beam))) task_queue = Queue() done_queue = Queue() # Submit tasks LOGGER.info('Submitting distributed tasks. This might take a while...') for task in TASKS: task_queue.put(task) - for i in range(num_processors): + for i in range(ncpu): Process(target=contour_worker, args=(task_queue, done_queue)).start() - source_list = [] num_max = 0 # Get the results from parallel processing for i in range(len(TASKS)): @@ -481,18 +249,9 @@ def generate_source_list(filename, outfile, mask_image, threshold_value, noise, task_queue.put('STOP') ra_sorted_list = sorted(source_list, key = itemgetter(0)) - LOGGER.info(f'Number of sources detected {len(ra_sorted_list)}') - catalog_out = f'# number of sources detected {len(ra_sorted_list)} \n' - f.write(catalog_out) - catalog_out = f'# number of souces using peak for source position: {num_max} \n' - f.write(catalog_out) - catalog_out = '#\n# source ra_hms dec_dms ra(deg) dec(deg) flux(mJy) error ang_size_(arcsec) pos_angle_(deg)\n' - f.write(catalog_out) - for i in range(len(ra_sorted_list)): - output = str(i) + ', ' + ra_sorted_list[i][1] + '\n' - f.write(output) - f.close() - return source_list + + return ra_sorted_list + def main(): LOGGER.info("Welcome to breizorro") @@ -729,7 +488,7 @@ 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: + if args.outcatalog or args.outregion: contours = find_contours(mask_image, 0.5) polygon_regions = [] for contour in contours: @@ -741,21 +500,52 @@ def load_fits_or_region(filename): 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.outcatalog: - num_proc = None # Assign based on the number of processors on a system - if args.ncpu: - num_proc = int(args.ncpu) + 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: + source_list = [] + image_data, hdu_header = get_image(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) noise = np.median(noise_image) - source_list = generate_source_list(args.imagename, args.outcatalog, mask_image, - threshold, noise, num_proc, mean_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') + + 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'# limiting_flux (mJy/beam): {round(limiting_flux*1000,2)} \n' + f.write(catalog_out) + source_list = multiprocess_contours(contours, args.ncpu) + catalog_out = f'# number of sources detected {len(source_list)} \n' + f.write(catalog_out) + catalog_out = f'# number of souces using peak for source position: {num_max} \n' + f.write(catalog_out) + catalog_out = '#\n# source ra_hms dec_dms ra(deg) dec(deg) flux(mJy) error ang_size_(arcsec) pos_angle_(deg)\n' + f.write(catalog_out) + for i in range(len(source_list)): + output = str(i) + ', ' + source_list[i][1] + '\n' + f.write(output) + f.close() LOGGER.info(f'Source catalog saved: {args.outcatalog}') if args.gui: @@ -796,15 +586,15 @@ def load_fits_or_region(filename): # Plot the image using degree coordinates p.image(image=[np.flip(mask_image, axis=1)], x=x_range[0], y=y_range[0], dw=extent, dh=extent, palette="Greys256", level="image") - - # Extracting data from source_list - x_coords = [float(d[1].split(', ')[2]) for d in source_list] - y_coords = [float(d[1].split(', ')[3]) for d in source_list] - labels = [f"({d[1].split(', ')[0]}, {d[1].split(', ')[1]})" for d in source_list] - # Plotting points - p.circle(x_coords, y_coords, size=3, color="red", alpha=0.5) - for i, (x, y, label) in enumerate(zip(x_coords, y_coords, labels)): - p.add_layout(Label(x=x, y=y, text=label, text_baseline="middle", text_align="center", text_font_size="10pt")) + if args.outcatalog: + # Extracting data from source_list + x_coords = [float(d[1].split(', ')[2]) for d in source_list] + y_coords = [float(d[1].split(', ')[3]) for d in source_list] + labels = [f"({d[1].split(', ')[0]}, {d[1].split(', ')[1]})" for d in source_list] + # Plotting points + p.circle(x_coords, y_coords, size=3, color="red", alpha=0.5) + for i, (x, y, label) in enumerate(zip(x_coords, y_coords, labels)): + p.add_layout(Label(x=x, y=y, text=label, text_baseline="middle", text_align="center", text_font_size="10pt")) p.grid.grid_line_width = 0.5 src1 = ColumnDataSource({'x': [], 'y': [], 'width': [], 'height': [], 'alpha': []}) src2 = ColumnDataSource({'xs': [], 'ys': [], 'alpha': []}) diff --git a/breizorro/utils.py b/breizorro/utils.py new file mode 100644 index 0000000..f7bab51 --- /dev/null +++ b/breizorro/utils.py @@ -0,0 +1,250 @@ +import sys +import numpy +import shutil +import logging +import argparse +import os.path +import re +import numpy as np + +from astropy.io import fits +from astropy.wcs import WCS +from astropy.coordinates import Angle +from astropy.coordinates import SkyCoord + + +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(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 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 of the ellipse in arcseconds. + bmin (float): Minor axis 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 calculate_weighted_centroid(x, y, flux_values): + # Calculate the total flux within the region + total_flux = np.sum(flux_values) + # Initialize variables for weighted sums + weighted_sum_x = 0 + weighted_sum_y = 0 + # Loop through all pixels within the region + for xi, yi, flux in zip(x, y, flux_values): + # Add the weighted contribution of each pixel to the centroid + weighted_sum_x += xi * flux + weighted_sum_y += yi * flux + # Calculate the centroid coordinates + centroid_x = weighted_sum_x / total_flux + centroid_y = weighted_sum_y / total_flux + return centroid_x, centroid_y + + +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} + + """ + hdu = fits.open(fitsname) + 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 maxDist(contour, pixel_size): + """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). + + Returns: + ang_size : float + Maximum extent of the contour in angular units (e.g., arcseconds). + pos_angle : float + Position angle of the contour (in degrees). + """ + src_size = 0 + pos_angle = None + + # Convert the contour to a numpy array for easier calculations + contour_array = np.array(contour) + + # Calculate pairwise distances between all points in the contour + for i in range(len(contour_array)): + for j in range(i+1, len(contour_array)): + # Calculate Euclidean distance between points i and j + distance = np.linalg.norm(contour_array[i] - contour_array[j]) * pixel_size + + # Calculate positional angle between points i and j + dx, dy = contour_array[j] - contour_array[i] + angle = np.degrees(np.arctan2(dy, dx)) + + # Update max_distance, max_points, and pos_angle if the calculated distance is greater + if distance > src_size: + src_size = distance + pos_angle = angle + + return src_size, 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 = b_major / pixel_size + b_pixels = b_minor / pixel_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): + result = maxDist(contour,pixel_size) + src_angle = result[0] + 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 <= mean_beam): + point_source = True + if source_beam_ratio <= 1.0: + point_source = True + if point_source: + src_size = (0.0, 0.0) + print(f"Point source because {int_peak_ratio} <= 0.2 and {src_angle} <= {mean_beam}") + else: + ang = round(src_angle,2) + pa = round(pos_angle,2) + src_size = (ang, pa) + return src_size From 4b488189161a2f86cdd2e871226a1f3ae415495a Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Fri, 22 Mar 2024 10:06:29 +0200 Subject: [PATCH 04/48] Add utils and init file --- breizorro/utils.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/breizorro/utils.py b/breizorro/utils.py index f7bab51..0a5e07e 100644 --- a/breizorro/utils.py +++ b/breizorro/utils.py @@ -12,6 +12,8 @@ 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 @@ -90,7 +92,7 @@ def calculate_area(bmaj, bmin, 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 + area = a_pixels * b_pixels return area From 3ad1d583c584a0a7a4425233a1023e367e0a47a1 Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Fri, 22 Mar 2024 10:07:00 +0200 Subject: [PATCH 05/48] source detection updates --- breizorro/breizorro.py | 82 +++++++++++++++++++++--------------------- 1 file changed, 41 insertions(+), 41 deletions(-) diff --git a/breizorro/breizorro.py b/breizorro/breizorro.py index c7670ca..04da87a 100644 --- a/breizorro/breizorro.py +++ b/breizorro/breizorro.py @@ -31,7 +31,7 @@ from operator import itemgetter, attrgetter from breizorro.utils import get_source_size, format_source_coordinates -from breizorro.utils import fitsInfo +from breizorro.utils import fitsInfo, calculate_area def create_logger(): @@ -136,39 +136,40 @@ def remove_regions(mask_image, regs, wcs): mask_image[reg.to_mask().to_image(mask_image.shape) != 0] = 0 -def process_contour(contour, image_data, fitsinfo, flux_cutoff, noise_out, beam): - - rr, cc = skimage_polygon(x,y) - data_result = image_data[rr, cc] - pixel_size = fitsinfo['ddec'] * 3600.0 - if fitsinfo['b_size']: - bmaj,bmin,_ = np.array(fitsinfo['b_size']) * 3600.0 - mean_beam = 0.5 * (bmaj + bmin) - else: - mean_beam = beam - pixels_beam = calculate_area(bmaj, bmin, pixel_size) - total_flux = data_result.sum() / pixels_beam - print(total_flux) - contained_points = len(rr) # needed for estimating flux density error +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_area(bmaj, bmin, 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:]) try: - peak_flux = data_result.max() + data = mask * image_data + nndata = np.flip(data, axis=0) + nndata = nndata[~np.isnan(nndata)] + total_flux = np.sum(nndata[nndata != -0.0])/pix_beam + peak_flux = data.max()/pix_beam except: LOGGER.warn('Failure to get maximum within contour') LOGGER.info('Probably a contour at edge of image - skipping') peak_flux = 0.0 - if peak_flux >= flux_cutoff: + if peak_flux: total_peak_ratio = np.abs((total_flux - peak_flux) / total_flux) - beam_error = contained_points/pixels_beam * noise_out - ten_pc_error = 0.1 * total_flux - flux_density_error = np.sqrt(ten_pc_error * ten_pc_error + beam_error * beam_error) - contour = [] - for i in range(len(x)): - contour.append([x[i],y[i]]) + flux_density_error = 0.001 + #beam_error = contained_points/pixels_beam * noise_out + #ten_pc_error = 0.1 * total_flux + #flux_density_error = np.sqrt(ten_pc_error * ten_pc_error + beam_error * beam_error) + #contour = [] + #for i in range(len(x)): + # contour.append([x[i],y[i]]) #centroid = calculate_weighted_centroid(x,y, data_result) #pix_centroid = PixCoord(centroid[0], centroid[1]) #contour_pixels = PixCoord(np.array(x), np.array(y)) @@ -180,18 +181,17 @@ def process_contour(contour, image_data, fitsinfo, flux_cutoff, noise_out, beam) if True: use_max = 1 LOGGER.warn('centroid lies outside polygon - using peak position') - location = np.unravel_index(np.argmax(data_result, axis=None), data_result.shape) - x_pos = rr[location] - y_pos = cc[location] - data_max = image_data[x_pos,y_pos] / pixels_beam - data_max = image_data[x_pos,y_pos] - pos_pixels = PixCoord(x_pos, y_pos) + location = np.unravel_index(np.argmax(data, axis=None), data.shape) + #x_pos = rr[location] + #y_pos = cc[location] + #data_max = image_data[x_pos,y_pos] / pixels_beam + #data_max = image_data[x_pos,y_pos] + pos_pixels = PixCoord(*location) ra, dec = wcs.all_pix2world(pos_pixels.x, pos_pixels.y, 0) - print(data_max) - print(f'{ra},{dec}') + ra, dec = float(ra), float(dec) source_flux = (round(total_flux * 1000, 3), round(flux_density_error * 1000, 4)) - source_size = get_source_size(contour, pixel_size, mean_beam, image_data, total_peak_ratio) + source_size = get_source_size(contour, pix_size, mean_beam, image_data, total_peak_ratio) source_pos = format_source_coordinates(ra, dec) source = source_pos + (ra, dec) + (total_flux, flux_density_error) + source_size catalog_out = ', '.join(str(src_prop) for src_prop in source) @@ -202,7 +202,7 @@ def process_contour(contour, image_data, fitsinfo, flux_cutoff, noise_out, beam) return (ra, catalog_out, use_max) -def multiprocess_contour(contours, ncpu=None): +def multiprocess_contours(contours, image_data, fitsinfo, mean_beam, ncpu=None): def contour_worker(input, output): for func, args in iter(input.get, 'STOP'): @@ -215,7 +215,7 @@ def contour_worker(input, output): try: import multiprocessing ncpu = multiprocessing.cpu_count() - LOGGER.info(f'Setting number of processors to {num_processors}') + LOGGER.info(f'Setting number of processors to {ncpu}') except: pass TASKS = [] @@ -227,7 +227,7 @@ def contour_worker(input, output): for j in range(len(contour)): x.append(contour[j][0]) y.append(contour[j][1]) - TASKS.append((process_contour,(contour, image_data, fitsinfo, limiting_flux, noise, mean_beam))) + TASKS.append((process_contour,(contour, image_data, fitsinfo, mean_beam))) task_queue = Queue() done_queue = Queue() # Submit tasks @@ -245,7 +245,7 @@ def contour_worker(input, output): source_list.append(catalog_out) num_max += catalog_out[2] # Tell child processes to stop - for i in range(num_processors): + for i in range(ncpu): task_queue.put('STOP') ra_sorted_list = sorted(source_list, key = itemgetter(0)) @@ -535,16 +535,16 @@ def load_fits_or_region(filename): limiting_flux = noise * threshold catalog_out = f'# limiting_flux (mJy/beam): {round(limiting_flux*1000,2)} \n' f.write(catalog_out) - source_list = multiprocess_contours(contours, args.ncpu) + source_list = multiprocess_contours(contours, image_data, fitsinfo, mean_beam, args.ncpu) catalog_out = f'# number of sources detected {len(source_list)} \n' f.write(catalog_out) - catalog_out = f'# number of souces using peak for source position: {num_max} \n' - f.write(catalog_out) - catalog_out = '#\n# source ra_hms dec_dms ra(deg) dec(deg) flux(mJy) error ang_size_(arcsec) pos_angle_(deg)\n' + #catalog_out = f'# number of souces using peak for source position: {num_max} \n' + #f.write(catalog_out) + catalog_out = '#\n# source ra_hms dec_dms ra(deg) dec(deg) flux(mJy) error ang_size(arcsec) pos_angle(deg)\n' f.write(catalog_out) for i in range(len(source_list)): output = str(i) + ', ' + source_list[i][1] + '\n' - f.write(output) + f.write(output) f.close() LOGGER.info(f'Source catalog saved: {args.outcatalog}') From 9b2bfb77b2aaacbbc3830e90af41e04341d60e8f Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Fri, 22 Mar 2024 10:07:28 +0200 Subject: [PATCH 06/48] Bump version --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 9ee5e6d..69530b1 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" From cf1e8e8fdb89a488537139ca955f62e5e291cdf7 Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Sun, 28 Apr 2024 21:58:03 +0200 Subject: [PATCH 07/48] positionns --- breizorro/breizorro.py | 33 +++++++++++++++---------- breizorro/utils.py | 55 +++++++++++++++++++++++++++++++++++++----- 2 files changed, 69 insertions(+), 19 deletions(-) diff --git a/breizorro/breizorro.py b/breizorro/breizorro.py index 04da87a..80105d7 100644 --- a/breizorro/breizorro.py +++ b/breizorro/breizorro.py @@ -30,8 +30,8 @@ from operator import itemgetter, attrgetter -from breizorro.utils import get_source_size, format_source_coordinates -from breizorro.utils import fitsInfo, calculate_area +from breizorro.utils import get_source_size, format_source_coordinates, deg2ra, deg2dec +from breizorro.utils import fitsInfo, calculate_beam_area def create_logger(): @@ -142,7 +142,7 @@ def process_contour(contour, image_data, fitsinfo, noise_out): pix_size = fitsinfo['ddec'] * 3600.0 bmaj, bmin,_ = np.array(fitsinfo['b_size']) * 3600.0 mean_beam = 0.5 * (bmaj + bmin) - pix_beam = calculate_area(bmaj, bmin, pix_size) + 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) @@ -153,15 +153,15 @@ def process_contour(contour, image_data, fitsinfo, noise_out): mask = pix_region.to_mask().to_image(image_data.shape[-2:]) try: data = mask * image_data - nndata = np.flip(data, axis=0) - nndata = nndata[~np.isnan(nndata)] + nndata = data # np.flip(data, axis=0) + #nndata = nndata[~np.isnan(nndata)] total_flux = np.sum(nndata[nndata != -0.0])/pix_beam - peak_flux = data.max()/pix_beam + peak_flux = nndata.max()/pix_beam except: LOGGER.warn('Failure to get maximum within contour') LOGGER.info('Probably a contour at edge of image - skipping') peak_flux = 0.0 - if peak_flux: + if total_flux: total_peak_ratio = np.abs((total_flux - peak_flux) / total_flux) flux_density_error = 0.001 #beam_error = contained_points/pixels_beam * noise_out @@ -181,7 +181,9 @@ def process_contour(contour, image_data, fitsinfo, noise_out): if True: use_max = 1 LOGGER.warn('centroid lies outside polygon - using peak position') - location = np.unravel_index(np.argmax(data, axis=None), data.shape) + location = list(np.unravel_index(np.argmax(nndata, axis=None), data.shape)) + location.reverse() + print(location) #x_pos = rr[location] #y_pos = cc[location] #data_max = image_data[x_pos,y_pos] / pixels_beam @@ -189,11 +191,14 @@ def process_contour(contour, image_data, fitsinfo, noise_out): pos_pixels = PixCoord(*location) ra, dec = wcs.all_pix2world(pos_pixels.x, pos_pixels.y, 0) ra, dec = float(ra), float(dec) + #import IPython; IPython.embed() source_flux = (round(total_flux * 1000, 3), round(flux_density_error * 1000, 4)) source_size = get_source_size(contour, pix_size, mean_beam, image_data, total_peak_ratio) source_pos = format_source_coordinates(ra, dec) - source = source_pos + (ra, dec) + (total_flux, flux_density_error) + source_size + print(source_pos) + print(total_flux) + source = source_pos + (ra, dec) + source_flux + source_size catalog_out = ', '.join(str(src_prop) for src_prop in source) else: # Dummy source to be eliminated @@ -219,6 +224,7 @@ def contour_worker(input, output): except: pass TASKS = [] + #import IPython; IPython.embed() for i in range(len(contours)): contour = contours[i] if len(contour) > 2: @@ -533,10 +539,10 @@ def load_fits_or_region(filename): catalog_out = f'# noise out (µJy/beam): {round(noise*1000000,2)} \n' f.write(catalog_out) limiting_flux = noise * threshold - catalog_out = f'# limiting_flux (mJy/beam): {round(limiting_flux*1000,2)} \n' + catalog_out = f'# cutt-off flux (mJy/beam): {round(limiting_flux*1000,2)} \n' f.write(catalog_out) source_list = multiprocess_contours(contours, image_data, fitsinfo, mean_beam, args.ncpu) - catalog_out = f'# number of sources detected {len(source_list)} \n' + catalog_out = f'# number of sources detected: {len(source_list)} \n' f.write(catalog_out) #catalog_out = f'# number of souces using peak for source position: {num_max} \n' #f.write(catalog_out) @@ -578,7 +584,8 @@ def load_fits_or_region(filename): # Define the plot p = figure(tooltips=[("x", "$x"), ("y", "$y"), ("value", "@image")], y_range=y_range) - p.x_range.range_padding = 0 + #p.x_range.range_padding = 0 + #p.y_range.range_padding = 0 p.x_range.flipped = True p.title.text = out_mask_fits @@ -606,7 +613,7 @@ def load_fits_or_region(filename): p.add_tools(draw_tool2) p.toolbar.active_drag = draw_tool1 output_file("breizorro.html", title="Mask Editor") - export_png(p, filename="breizorro.png") + export_png(p, filename="breizorro1.png") show(p) LOGGER.info(f"Enforcing that mask to binary") diff --git a/breizorro/utils.py b/breizorro/utils.py index 0a5e07e..a8c68c7 100644 --- a/breizorro/utils.py +++ b/breizorro/utils.py @@ -73,6 +73,49 @@ def format_source_coordinates(coord_ra_deg, coord_dec_deg): 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): """ @@ -80,8 +123,8 @@ def calculate_area(bmaj, bmin, pix_size): given the pixel size. Parameters: - bmaj (float): Major axis of the ellipse in arcseconds. - bmin (float): Minor axis of the ellipse in arcseconds. + 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: @@ -92,7 +135,7 @@ def calculate_area(bmaj, bmin, pix_size): b_pixels = bmin / pix_size # Calculate the area of the ellipse using the formula: π * a * b - area = a_pixels * b_pixels + area = np.pi * a_pixels * b_pixels return area @@ -220,11 +263,11 @@ def calculate_beam_area(bmaj, bmin, pix_size): area (float): Calculated area of the ellipse in square pixels. """ # Calculate the semi-major and semi-minor axes in pixels - a_pixels = b_major / pixel_size - b_pixels = b_minor / pixel_size + 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 + area = np.pi * a_pixels * b_pixels return area From d574bf3c5b28af056c97ab21e41ace1b0e17f21f Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Thu, 12 Sep 2024 14:13:52 +0200 Subject: [PATCH 08/48] Refactor routines --- breizorro/breizorro.py | 263 +++++++---------------------------------- breizorro/catalog.py | 104 ++++++++++++++++ breizorro/gui.py | 169 ++++++++++++++++++++++++++ breizorro/utils.py | 90 ++++++++------ 4 files changed, 369 insertions(+), 257 deletions(-) create mode 100644 breizorro/catalog.py create mode 100644 breizorro/gui.py diff --git a/breizorro/breizorro.py b/breizorro/breizorro.py index 80105d7..d8a5c8f 100644 --- a/breizorro/breizorro.py +++ b/breizorro/breizorro.py @@ -7,6 +7,7 @@ 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 @@ -26,12 +27,8 @@ from skimage.draw import polygon as skimage_polygon from skimage.measure import find_contours -from multiprocessing import Process, Queue - -from operator import itemgetter, attrgetter - from breizorro.utils import get_source_size, format_source_coordinates, deg2ra, deg2dec -from breizorro.utils import fitsInfo, calculate_beam_area +from breizorro.utils import get_image_data, fitsInfo, calculate_beam_area def create_logger(): @@ -47,40 +44,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): @@ -136,129 +113,6 @@ def remove_regions(mask_image, regs, wcs): mask_image[reg.to_mask().to_image(mask_image.shape) != 0] = 0 -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:]) - 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: - LOGGER.warn('Failure to get maximum within contour') - LOGGER.info('Probably a contour at edge of image - skipping') - peak_flux = 0.0 - if total_flux: - total_peak_ratio = np.abs((total_flux - peak_flux) / total_flux) - flux_density_error = 0.001 - #beam_error = contained_points/pixels_beam * noise_out - #ten_pc_error = 0.1 * total_flux - #flux_density_error = np.sqrt(ten_pc_error * ten_pc_error + beam_error * beam_error) - #contour = [] - #for i in range(len(x)): - # contour.append([x[i],y[i]]) - #centroid = calculate_weighted_centroid(x,y, data_result) - #pix_centroid = PixCoord(centroid[0], centroid[1]) - #contour_pixels = PixCoord(np.array(x), np.array(y)) - #p = PolygonPixelRegion(vertices=contour_pixels, meta={'label': 'Region'}) - #ra, dec = wcs.all_pix2world(pix_centroid.x, pix_centroid.y,0) - #if p.contains(pix_centroid)[0]: - # ra, dec = wcs.all_pix2world(pix_centroid.x, pix_centroid.y,0) - #else: - if True: - use_max = 1 - LOGGER.warn('centroid lies outside polygon - using peak position') - location = list(np.unravel_index(np.argmax(nndata, axis=None), data.shape)) - location.reverse() - print(location) - #x_pos = rr[location] - #y_pos = cc[location] - #data_max = image_data[x_pos,y_pos] / pixels_beam - #data_max = image_data[x_pos,y_pos] - pos_pixels = PixCoord(*location) - ra, dec = wcs.all_pix2world(pos_pixels.x, pos_pixels.y, 0) - ra, dec = float(ra), float(dec) - #import IPython; IPython.embed() - - source_flux = (round(total_flux * 1000, 3), round(flux_density_error * 1000, 4)) - source_size = get_source_size(contour, pix_size, mean_beam, image_data, total_peak_ratio) - source_pos = format_source_coordinates(ra, dec) - print(source_pos) - print(total_flux) - source = source_pos + (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, mean_beam, 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() - LOGGER.info(f'Setting number of processors to {ncpu}') - except: - pass - TASKS = [] - #import IPython; IPython.embed() - 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, mean_beam))) - task_queue = Queue() - done_queue = Queue() - # Submit tasks - LOGGER.info('Submitting distributed tasks. This might take a while...') - 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=300) - 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 - - def main(): LOGGER.info("Welcome to breizorro") # Get version @@ -309,7 +163,7 @@ 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, + 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.') @@ -343,7 +197,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) @@ -365,7 +219,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}" @@ -382,7 +236,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) @@ -493,6 +347,10 @@ 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: contours = find_contours(mask_image, 0.5) @@ -512,13 +370,22 @@ def load_fits_or_region(filename): 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(args.imagename) + 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) - noise = np.median(noise_image) if mean_beam: LOGGER.info(f'Using user provided size: {mean_beam}') elif fitsinfo['b_size']: @@ -527,6 +394,7 @@ def load_fits_or_region(filename): 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) @@ -541,80 +409,29 @@ def load_fits_or_region(filename): 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, mean_beam, 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 = f'# number of souces using peak for source position: {num_max} \n' - #f.write(catalog_out) - catalog_out = '#\n# source ra_hms dec_dms ra(deg) dec(deg) flux(mJy) error ang_size(arcsec) pos_angle(deg)\n' + 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 = str(i) + ', ' + source_list[i][1] + '\n' + 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 Label, BoxEditTool, ColumnDataSource, FreehandDrawTool - from bokeh.plotting import figure, output_file, show - from bokeh.io import curdoc, export_png - 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 ...") - fitsinfo = fitsInfo(args.imagename or args.maskname) - - # Origin coordinates - origin_ra, origin_dec = fitsinfo['centre'] - - # Pixel width in degrees - pixel_width = fitsinfo['ddec'] - - # Calculate the extent of the image in degrees - # We assume a square image for simplicity - extent = fitsinfo['numPix'] * pixel_width # Assume equal pixels in each dimension - - # Specify the coordinates for the image - x_range = (origin_ra - extent/2.0, origin_ra + extent/2.0) - y_range = (origin_dec - extent/2.0, origin_dec + extent/2.0) - - # Define the plot - p = figure(tooltips=[("x", "$x"), ("y", "$y"), ("value", "@image")], y_range=y_range) - #p.x_range.range_padding = 0 - #p.y_range.range_padding = 0 - p.x_range.flipped = True - p.title.text = out_mask_fits - - # must give a vector of image data for image parameter - # Plot the image using degree coordinates - p.image(image=[np.flip(mask_image, axis=1)], x=x_range[0], y=y_range[0], dw=extent, dh=extent, palette="Greys256", level="image") - - if args.outcatalog: - # Extracting data from source_list - x_coords = [float(d[1].split(', ')[2]) for d in source_list] - y_coords = [float(d[1].split(', ')[3]) for d in source_list] - labels = [f"({d[1].split(', ')[0]}, {d[1].split(', ')[1]})" for d in source_list] - # Plotting points - p.circle(x_coords, y_coords, size=3, color="red", alpha=0.5) - for i, (x, y, label) in enumerate(zip(x_coords, y_coords, labels)): - p.add_layout(Label(x=x, y=y, text=label, text_baseline="middle", text_align="center", text_font_size="10pt")) - 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") - export_png(p, filename="breizorro1.png") - show(p) + 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 diff --git a/breizorro/catalog.py b/breizorro/catalog.py new file mode 100644 index 0000000..a100e80 --- /dev/null +++ b/breizorro/catalog.py @@ -0,0 +1,104 @@ +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:]) + #contained_points = len(polygon_region.vertices.y) #len(rr) # needed for estimating flux density error + 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 = 0.0001 + #beam_error = contained_points/pixels_beam * noise_out + #ten_pc_error = 0.1 * total_flux + #flux_density_error = np.sqrt(ten_pc_error * ten_pc_error + beam_error * beam_error) + # Calculate weighted centroid + centroid_x, centroid_y = centroids.centroid_2dg(data) + 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) + #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, mean_beam, 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, mean_beam))) + 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..d7c4a36 --- /dev/null +++ b/breizorro/gui.py @@ -0,0 +1,169 @@ +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['ddec'] + # Calculate the extent of the image in degrees + # We assume a square image for simplicity + extent = fitsinfo['numPix'] * pixel_width # 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/2.0, origin_ra + extent/2.0) + y_range = (origin_dec - extent/2.0, origin_dec + extent/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="RA (deg)", + y_axis_label="DEC (deg)", + y_range=(min(y_coords)-0.1, max(y_coords)+0.1), # Add margin + 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.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"]) + # 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) + output_file("breizorro.html", title="Mask Editor") + #export_png(p, filename="breizorro.png") + show(layout) diff --git a/breizorro/utils.py b/breizorro/utils.py index a8c68c7..221584c 100644 --- a/breizorro/utils.py +++ b/breizorro/utils.py @@ -7,6 +7,7 @@ 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 @@ -31,7 +32,13 @@ def deg_to_dms(dec_deg): return degrees, minutes, seconds -def format_source_coordinates(coord_ra_deg, coord_dec_deg): +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) @@ -139,22 +146,6 @@ def calculate_area(bmaj, bmin, pix_size): return area -def calculate_weighted_centroid(x, y, flux_values): - # Calculate the total flux within the region - total_flux = np.sum(flux_values) - # Initialize variables for weighted sums - weighted_sum_x = 0 - weighted_sum_y = 0 - # Loop through all pixels within the region - for xi, yi, flux in zip(x, y, flux_values): - # Add the weighted contribution of each pixel to the centroid - weighted_sum_x += xi * flux - weighted_sum_y += yi * flux - # Calculate the centroid coordinates - centroid_x = weighted_sum_x / total_flux - centroid_y = weighted_sum_y / total_flux - return centroid_x, centroid_y - def fitsInfo(fitsname=None): """Get fits header info. @@ -175,8 +166,8 @@ def fitsInfo(fitsname=None): 'skyArea': skyArea, 'naxis': naxis} """ - hdu = fits.open(fitsname) - hdr = hdu[0].header + with fits.open(fitsname) as hdu: + hdr = hdu[0].header ra = hdr['CRVAL1'] dra = abs(hdr['CDELT1']) raPix = hdr['CRPIX1'] @@ -210,6 +201,30 @@ def fitsInfo(fitsname=None): '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): """Calculate maximum extent and position angle of a contour. @@ -220,18 +235,21 @@ def maxDist(contour, pixel_size): Size of a pixel in the image (e.g., arcseconds per pixel). Returns: - ang_size : float - Maximum extent of the contour in angular units (e.g., arcseconds). + 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). """ - src_size = 0 + e_maj = 0 + e_min = np.inf # Start with a very large value for the minor axis pos_angle = None # Convert the contour to a numpy array for easier calculations contour_array = np.array(contour) - # Calculate pairwise distances between all points in the contour + # Step 1: Loop through all pairs of points to find major and minor axes for i in range(len(contour_array)): for j in range(i+1, len(contour_array)): # Calculate Euclidean distance between points i and j @@ -241,12 +259,16 @@ def maxDist(contour, pixel_size): dx, dy = contour_array[j] - contour_array[i] angle = np.degrees(np.arctan2(dy, dx)) - # Update max_distance, max_points, and pos_angle if the calculated distance is greater - if distance > src_size: - src_size = distance + # Update e_maj and pos_angle if this distance is the largest + if distance > e_maj: + e_maj = distance pos_angle = angle - return src_size, pos_angle + # Update e_min if this distance is the smallest + if distance < e_min: + e_min = distance + + return e_maj, e_min, pos_angle def calculate_beam_area(bmaj, bmin, pix_size): @@ -274,22 +296,22 @@ def calculate_beam_area(bmaj, bmin, pix_size): def get_source_size(contour, pixel_size, mean_beam, image, int_peak_ratio): result = maxDist(contour,pixel_size) - src_angle = result[0] - pos_angle = result[1] + 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 <= mean_beam): + 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) - print(f"Point source because {int_peak_ratio} <= 0.2 and {src_angle} <= {mean_beam}") + src_size = (0.0, 0.0, 0.0) else: - ang = round(src_angle,2) + emaj = round(src_angle[0],2) + emin = round(src_angle[-1],2) pa = round(pos_angle,2) - src_size = (ang, pa) + src_size = (emaj, emin, pa) return src_size From ece367a7e8618ec4f8cfc96e6132f143b18ecea2 Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Sat, 19 Oct 2024 11:06:51 +0200 Subject: [PATCH 09/48] Rename test_installation.yml to python-package.yml --- .github/workflows/{test_installation.yml => python-package.yml} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename .github/workflows/{test_installation.yml => python-package.yml} (100%) diff --git a/.github/workflows/test_installation.yml b/.github/workflows/python-package.yml similarity index 100% rename from .github/workflows/test_installation.yml rename to .github/workflows/python-package.yml From e8e6075ff28edcf6d3745d8fe155d80b429ceb52 Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Sun, 20 Oct 2024 12:03:04 +0200 Subject: [PATCH 10/48] Create _config.yml --- _config.yml | 5 +++++ 1 file changed, 5 insertions(+) create mode 100644 _config.yml diff --git a/_config.yml b/_config.yml new file mode 100644 index 0000000..7a7df3b --- /dev/null +++ b/_config.yml @@ -0,0 +1,5 @@ +title: Midnight theme +description: Midnight is a theme for GitHub Pages. +show_downloads: true +google_analytics: +theme: jekyll-theme-midnight From d6a9beba2fa9af71ed44bb16991d47d3c6d07d3a Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Sun, 20 Oct 2024 12:12:29 +0200 Subject: [PATCH 11/48] Update _config.yml --- _config.yml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/_config.yml b/_config.yml index 7a7df3b..17baf7c 100644 --- a/_config.yml +++ b/_config.yml @@ -1,5 +1,5 @@ -title: Midnight theme -description: Midnight is a theme for GitHub Pages. +title: Breizorro +description: Breizorro is an image masking and cataloguing show_downloads: true -google_analytics: +google_analytics: theme: jekyll-theme-midnight From e5798d22d41b097f7bf05daeb9f3fcfd568609b8 Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Sun, 20 Oct 2024 12:27:48 +0200 Subject: [PATCH 12/48] Create index.md --- index.md | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) create mode 100644 index.md diff --git a/index.md b/index.md new file mode 100644 index 0000000..0d8b3ad --- /dev/null +++ b/index.md @@ -0,0 +1,25 @@ + +--- +layout: default +--- + +![Octocat](https://github.githubassets.com/images/icons/emoji/octocat.png) + +Breizorro is a flexible software program made to simplify image analysis tasks, including identifying emission islands and generating and modifying picture masks, which are frequently used in radio interferometric imaging. + +# 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 a picture, 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. + +# Region Generation and Manipulation + +Breizorro makes it easier to create and work with regions using image masks. Labeling, eliminating, extracting, and filtering regions (islands) based on user-specified criteria are all included in this. Users can employ techniques including erosion, dilation, hole-filling, binary masking, and inversion to refine their regions of interest. + +# Cataloguing and Visualization + +Breizorro enables catalog generation from extracted regions, saving source properties to an ASCII/text +file. + +``` +breizorro -r circinus-MFS-image.fits --save-catalog circinus.txt +``` From 1d182ae849bdd0b8caf10c4d7514bdc9c3776a57 Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Sun, 20 Oct 2024 12:42:26 +0200 Subject: [PATCH 13/48] Update index.md --- index.md | 48 +++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 45 insertions(+), 3 deletions(-) diff --git a/index.md b/index.md index 0d8b3ad..071e418 100644 --- a/index.md +++ b/index.md @@ -1,7 +1,4 @@ ---- -layout: default ---- ![Octocat](https://github.githubassets.com/images/icons/emoji/octocat.png) @@ -11,15 +8,60 @@ Breizorro is a flexible software program made to simplify image analysis tasks, 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 a picture, 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. Labeling, eliminating, extracting, and filtering regions (islands) based on user-specified criteria are all included in this. Users can employ techniques including erosion, dilation, hole-filling, binary masking, and inversion to refine their regions of interest. +```breizorro -r circinus-MFS-image.fits --save-regions 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 catalog generation from extracted regions, saving source properties to an ASCII/text file. +```breizorro -r circinus-MFS-image.fits --save-catalog circinus.txt``` + +``` +# processing fits image: circinus-MFS-image.fits +# mean beam size (arcsec): 37.97 +# original image peak flux (Jy/beam): 0.768315315246582 +# noise out (μJy/beam): 737.33 +# cutt-off flux (mJy/beam): 4.79 +# freq0 (Hz): 1419944335.9375 +# number of sources detected: 35 +# +#format: name ra_d dec_d i i_err emaj_s emin_s pa_d +src0 -147.88904620760084 -65.52043870236054 0.00468 0.0001 50.04 0.0 177.71 +src1 -147.84162114219356 -65.44040921224196 0.00047 0.0001 0.0 0.0 0.0 +src2 -147.7980562104931 -65.39402001087845 0.00089 0.0001 0.0 0.0 0.0 +src3 -147.74069959466678 -65.08592423258469 0.0082 0.0001 51.61 0.0 144.46 +src4 -147.6583486798732 -65.28753961408073 0.09922 0.0001 86.58 0.0 -173.37 +src5 -147.59829620974554 -65.28104296056243 0.00629 0.0001 47.07 0.0 167.74 +src6 -147.51788583319714 -65.72316001301358 0.06639 0.0001 71.47 0.0 162.07 +src7 -147.49179712006742 -64.88584736668952 0.00092 0.0001 0.0 0.0 0.0 +src8 -147.42970139783463 -65.65050232604096 0.00069 0.0001 0.0 0.0 0.0 +src9 -147.4293961973031 -65.44158144271519 0.00134 0.0001 0.0 0.0 0.0 +src10 -147.28370646739054 -65.42936506202037 9e-05 0.0001 0.0 0.0 0.0 +``` + ``` breizorro -r circinus-MFS-image.fits --save-catalog circinus.txt ``` + +![Screenshot 2024-09-11 100116](https://github.com/user-attachments/assets/79d3ece6-5d96-48a1-a06e-fbae2987d333) From dd51336320bd019eef896955488574b7593be12a Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Sun, 20 Oct 2024 12:48:28 +0200 Subject: [PATCH 14/48] Update README.rst --- README.rst | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/README.rst b/README.rst index bf96a67..839468d 100644 --- a/README.rst +++ b/README.rst @@ -88,6 +88,12 @@ To show help message and exit 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 +122,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 From e46f763ae870e97502fb7309257945ba4d1b094a Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Sun, 20 Oct 2024 12:49:44 +0200 Subject: [PATCH 15/48] Update index.md --- index.md | 22 +++++++++++++++------- 1 file changed, 15 insertions(+), 7 deletions(-) diff --git a/index.md b/index.md index 071e418..2651dac 100644 --- a/index.md +++ b/index.md @@ -1,5 +1,3 @@ - - ![Octocat](https://github.githubassets.com/images/icons/emoji/octocat.png) Breizorro is a flexible software program made to simplify image analysis tasks, including identifying emission islands and generating and modifying picture masks, which are frequently used in radio interferometric imaging. @@ -8,12 +6,16 @@ Breizorro is a flexible software program made to simplify image analysis tasks, 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 a picture, 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``` +``` +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``` +``` +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) @@ -21,11 +23,15 @@ Breizorro uses the minimal filter to generate binary masks by replacing each pix Breizorro makes it easier to create and work with regions using image masks. Labeling, eliminating, extracting, and filtering regions (islands) based on user-specified criteria are all included in this. Users can employ techniques including erosion, dilation, hole-filling, binary masking, and inversion to refine their regions of interest. -```breizorro -r circinus-MFS-image.fits --save-regions circinus.reg``` +``` +breizorro -r circinus-MFS-image.fits --save-regions 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``` +``` +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) @@ -35,7 +41,9 @@ Breizorro makes it easier to create and work with regions using image masks. Lab Breizorro enables catalog generation from extracted regions, saving source properties to an ASCII/text file. -```breizorro -r circinus-MFS-image.fits --save-catalog circinus.txt``` +``` +breizorro -r circinus-MFS-image.fits --save-catalog circinus.txt +``` ``` # processing fits image: circinus-MFS-image.fits From 53e7408614b264c27c03ef07315dad3ed7b953db Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Sun, 20 Oct 2024 13:26:02 +0200 Subject: [PATCH 16/48] Update index.md --- index.md | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/index.md b/index.md index 2651dac..5516f7c 100644 --- a/index.md +++ b/index.md @@ -1,4 +1,4 @@ -![Octocat](https://github.githubassets.com/images/icons/emoji/octocat.png) +![breiz](https://github.com/user-attachments/assets/9f2376eb-37a1-4f73-8675-3c92b5da1e99) Breizorro is a flexible software program made to simplify image analysis tasks, including identifying emission islands and generating and modifying picture masks, which are frequently used in radio interferometric imaging. @@ -21,7 +21,8 @@ breizorro -r circinus-MFS-image.fits --sum-peak 500 # Region Generation and Manipulation -Breizorro makes it easier to create and work with regions using image masks. Labeling, eliminating, extracting, and filtering regions (islands) based on user-specified criteria are all included in this. Users can employ techniques including erosion, dilation, hole-filling, binary masking, and inversion to refine their regions of interest. +Breizorro makes it easier to create and work with regions using image masks. Labeling, eliminating, extracting, and filtering regions (islands) based on user-specified criteria are all included in this. Us![breiz](https://github.com/user-attachments/assets/7b1f1a68-d0df-45fe-88eb-61800af3bf19) +ers can employ techniques including erosion, dilation, hole-filling, binary masking, and inversion to refine their regions of interest. ``` breizorro -r circinus-MFS-image.fits --save-regions circinus.reg From 896ce23c8515c1068e99d9023f9f2667a31d9f36 Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Sun, 20 Oct 2024 13:26:51 +0200 Subject: [PATCH 17/48] Update index.md --- index.md | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/index.md b/index.md index 5516f7c..db91f47 100644 --- a/index.md +++ b/index.md @@ -21,8 +21,7 @@ breizorro -r circinus-MFS-image.fits --sum-peak 500 # Region Generation and Manipulation -Breizorro makes it easier to create and work with regions using image masks. Labeling, eliminating, extracting, and filtering regions (islands) based on user-specified criteria are all included in this. Us![breiz](https://github.com/user-attachments/assets/7b1f1a68-d0df-45fe-88eb-61800af3bf19) -ers can employ techniques including erosion, dilation, hole-filling, binary masking, and inversion to refine their regions of interest. +Breizorro makes it easier to create and work with regions using image masks. It includes labeling, 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 --save-regions circinus.reg From 015dd709b39ed67577dae16538dd75552524fb45 Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Sun, 20 Oct 2024 13:33:51 +0200 Subject: [PATCH 18/48] Update index.md --- index.md | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/index.md b/index.md index db91f47..bcd32f1 100644 --- a/index.md +++ b/index.md @@ -1,4 +1,6 @@ -![breiz](https://github.com/user-attachments/assets/9f2376eb-37a1-4f73-8675-3c92b5da1e99) +

+ Breizorro Logo +

Breizorro is a flexible software program made to simplify image analysis tasks, including identifying emission islands and generating and modifying picture masks, which are frequently used in radio interferometric imaging. From 1b6a99d3460db099d3b23573b26836c943d0fcd9 Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Sun, 20 Oct 2024 14:10:40 +0200 Subject: [PATCH 19/48] Update README.rst --- README.rst | 58 ------------------------------------------------------ 1 file changed, 58 deletions(-) diff --git a/README.rst b/README.rst index 839468d..390605a 100644 --- a/README.rst +++ b/README.rst @@ -30,64 +30,6 @@ 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 ============= From 64cbe402806923a7e5c340c896255bf938a0dbf1 Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Sun, 20 Oct 2024 14:13:43 +0200 Subject: [PATCH 20/48] Update index.md --- index.md | 51 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 51 insertions(+) diff --git a/index.md b/index.md index bcd32f1..5ebb1fb 100644 --- a/index.md +++ b/index.md @@ -4,6 +4,57 @@ Breizorro is a flexible software program made to simplify image analysis tasks, including identifying emission islands and generating and modifying picture masks, which are frequently used in radio interferometric imaging. + +# Parameter definition + +``` + 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. +``` + + # 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 a picture, 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. From f0daab2f048823b8d64472875a1ce83942b147ce Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Sun, 20 Oct 2024 14:23:56 +0200 Subject: [PATCH 21/48] Update index.md --- index.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/index.md b/index.md index 5ebb1fb..7cf5b5d 100644 --- a/index.md +++ b/index.md @@ -2,7 +2,7 @@ Breizorro Logo

-Breizorro is a flexible software program made to simplify image analysis tasks, including identifying emission islands and generating and modifying picture masks, which are frequently used in radio interferometric imaging. +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 From 975ff99cf2308b4415bc12134d517c56190948e2 Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Sun, 20 Oct 2024 15:36:36 +0200 Subject: [PATCH 22/48] Update _config.yml --- _config.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/_config.yml b/_config.yml index 17baf7c..ee9483e 100644 --- a/_config.yml +++ b/_config.yml @@ -1,5 +1,5 @@ title: Breizorro -description: Breizorro is an image masking and cataloguing +description: Image masking and cataloguing suite show_downloads: true google_analytics: theme: jekyll-theme-midnight From 36a4e78d9b925b5e43877f675d90f70cd4a5001c Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Sun, 20 Oct 2024 20:33:04 +0200 Subject: [PATCH 23/48] Update index.md --- index.md | 88 ++++++++++++++++++++++++++++---------------------------- 1 file changed, 44 insertions(+), 44 deletions(-) diff --git a/index.md b/index.md index 7cf5b5d..9cdd6ce 100644 --- a/index.md +++ b/index.md @@ -8,50 +8,50 @@ Breizorro is a flexible software program made to simplify image analysis tasks, # Parameter definition ``` - 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. +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. ``` From b36c3b37c7c3cd41026aaa1295172f58744bef53 Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Sun, 20 Oct 2024 20:56:42 +0200 Subject: [PATCH 24/48] Update index.md --- index.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/index.md b/index.md index 9cdd6ce..a21ff29 100644 --- a/index.md +++ b/index.md @@ -57,7 +57,7 @@ optional arguments: # 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 a picture, 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 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 From 660f5b15bdade3a4b996eab35184c8f839037c37 Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Sun, 20 Oct 2024 21:06:08 +0200 Subject: [PATCH 25/48] Update index.md --- index.md | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/index.md b/index.md index a21ff29..c158eb9 100644 --- a/index.md +++ b/index.md @@ -126,3 +126,10 @@ breizorro -r circinus-MFS-image.fits --save-catalog circinus.txt ``` ![Screenshot 2024-09-11 100116](https://github.com/user-attachments/assets/79d3ece6-5d96-48a1-a06e-fbae2987d333) + + +# 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) From 83acfa1e5f22535d0f5a69d08d57adb2aafc6c14 Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Tue, 22 Oct 2024 11:50:23 +0200 Subject: [PATCH 26/48] legend color --- breizorro/gui.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/breizorro/gui.py b/breizorro/gui.py index d7c4a36..49d3edd 100644 --- a/breizorro/gui.py +++ b/breizorro/gui.py @@ -122,7 +122,7 @@ def display(imagename, mask_image, outcatalog, source_list): p = figure(title="Breizorro Source Catalog", x_axis_label="RA (deg)", y_axis_label="DEC (deg)", - y_range=(min(y_coords)-0.1, max(y_coords)+0.1), # Add margin + y_range=(min(y_coords), max(y_coords)), match_aspect=True, tooltips=[("x", "$x"), ("y", "$y"), ("value", "@image")]) # Plot the image @@ -149,6 +149,7 @@ def display(imagename, mask_image, outcatalog, source_list): 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) From f1aa77e14de005003a0c30d6e488c1edfa9fac1d Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Tue, 22 Oct 2024 13:05:13 +0200 Subject: [PATCH 27/48] Add files via upload --- breizorro.html | 63 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 63 insertions(+) create mode 100644 breizorro.html 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 From bf9d06d604372761478587725911e1f9e72bc1a1 Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Tue, 22 Oct 2024 13:07:49 +0200 Subject: [PATCH 28/48] Update index.md --- index.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/index.md b/index.md index c158eb9..bb44aa2 100644 --- a/index.md +++ b/index.md @@ -122,11 +122,12 @@ src10 -147.28370646739054 -65.42936506202037 9e-05 0.0001 0.0 0.0 0.0 ``` ``` -breizorro -r circinus-MFS-image.fits --save-catalog circinus.txt +breizorro -r circinus-MFS-image.fits --save-catalog circinus.txt --gui ``` ![Screenshot 2024-09-11 100116](https://github.com/user-attachments/assets/79d3ece6-5d96-48a1-a06e-fbae2987d333) +For a live demo of the viewer, visit: [breiz-viewer](breizorro.html) # Contributors From 337f6923fdc152b6c88c9e4c843a5e5145bd3cec Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Sun, 27 Oct 2024 17:44:45 +0200 Subject: [PATCH 29/48] plot extent and sizing mode --- breizorro/gui.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/breizorro/gui.py b/breizorro/gui.py index 49d3edd..8fe2ca5 100644 --- a/breizorro/gui.py +++ b/breizorro/gui.py @@ -55,15 +55,17 @@ def display(imagename, mask_image, outcatalog, source_list): # Origin coordinates origin_ra, origin_dec = fitsinfo['centre'] # Pixel width in degrees - pixel_width = fitsinfo['ddec'] + pixel_width = fitsinfo['dra'] + pixel_height = fitsinfo['ddec'] # Calculate the extent of the image in degrees # We assume a square image for simplicity - extent = fitsinfo['numPix'] * pixel_width # Assume equal pixels in each dimension + 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/2.0, origin_ra + extent/2.0) - y_range = (origin_dec - extent/2.0, origin_dec + extent/2.0) + 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] @@ -164,7 +166,7 @@ def display(imagename, mask_image, outcatalog, source_list): """) # Connect the hover tool to the callback hover.callback = callback - layout = row(p, data_table) + layout = row(p, data_table, sizing_mode="stretch_both") output_file("breizorro.html", title="Mask Editor") #export_png(p, filename="breizorro.png") show(layout) From 9f6b8cad27570e381515fe5fd30b0311a1d0c34b Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Sun, 27 Oct 2024 18:49:50 +0200 Subject: [PATCH 30/48] Add photutils as deps for cataloging --- pyproject.toml | 3 +++ 1 file changed, 3 insertions(+) diff --git a/pyproject.toml b/pyproject.toml index 69530b1..ab32ee9 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -32,6 +32,7 @@ scipy = "*" # Optional dependencies start here bokeh = { version = "*", optional = true} +photutils = { version = "*", optional = true} pytest = { version = "*", optional=true } pytest-flake8 = { version = "*", optional=true } @@ -40,7 +41,9 @@ requires = ["setuptools", "poetry-core"] build-backend = "poetry.core.masonry.api" [tool.poetry.extras] +all = ["bokeh", "photutils"] gui = ["bokeh"] +catalog = ["photutils"] testing = ["pytest", "pytest-flake8"] [tool.poetry.group.tests] From c80f3bb007913bc5d5e9099b79973aa030f509eb Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Wed, 30 Oct 2024 12:51:40 +0200 Subject: [PATCH 31/48] format axis --- breizorro/gui.py | 24 ++++++++++++++++++++++-- 1 file changed, 22 insertions(+), 2 deletions(-) diff --git a/breizorro/gui.py b/breizorro/gui.py index 8fe2ca5..448078a 100644 --- a/breizorro/gui.py +++ b/breizorro/gui.py @@ -122,8 +122,9 @@ def display(imagename, mask_image, outcatalog, source_list): # Create the plot p = figure(title="Breizorro Source Catalog", - x_axis_label="RA (deg)", - y_axis_label="DEC (deg)", + 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")]) @@ -159,6 +160,25 @@ def display(imagename, mask_image, outcatalog, source_list): #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 From 3b3eca615d208ae1f2a6ff6f6204aad173a4c99f Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Thu, 7 Nov 2024 19:54:39 +0200 Subject: [PATCH 32/48] fit source shape --- breizorro/breizorro.py | 2 +- breizorro/catalog.py | 22 ++++++++++++--------- breizorro/utils.py | 43 ++++++++++++++++++------------------------ 3 files changed, 32 insertions(+), 35 deletions(-) diff --git a/breizorro/breizorro.py b/breizorro/breizorro.py index d8a5c8f..05097eb 100644 --- a/breizorro/breizorro.py +++ b/breizorro/breizorro.py @@ -410,7 +410,7 @@ def load_fits_or_region(filename): 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, mean_beam, args.ncpu) + 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' diff --git a/breizorro/catalog.py b/breizorro/catalog.py index a100e80..0b2e784 100644 --- a/breizorro/catalog.py +++ b/breizorro/catalog.py @@ -22,7 +22,10 @@ def process_contour(contour, image_data, fitsinfo, noise_out): polygon_region = PolygonSkyRegion(vertices=contour_sky) pix_region = polygon_region.to_pixel(wcs) mask = pix_region.to_mask().to_image(image_data.shape[-2:]) - #contained_points = len(polygon_region.vertices.y) #len(rr) # needed for estimating flux density error + # 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) @@ -33,18 +36,19 @@ def process_contour(contour, image_data, fitsinfo, noise_out): peak_flux = 0.0 if total_flux: total_peak_ratio = np.abs((total_flux - peak_flux) / total_flux) - flux_density_error = 0.0001 - #beam_error = contained_points/pixels_beam * noise_out - #ten_pc_error = 0.1 * total_flux - #flux_density_error = np.sqrt(ten_pc_error * ten_pc_error + beam_error * beam_error) + # 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 - centroid_x, centroid_y = centroids.centroid_2dg(data) + _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) + 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) @@ -55,7 +59,7 @@ def process_contour(contour, image_data, fitsinfo, noise_out): return (ra, catalog_out, use_max) -def multiprocess_contours(contours, image_data, fitsinfo, mean_beam, ncpu=None): +def multiprocess_contours(contours, image_data, fitsinfo, noise_out, ncpu=None): def contour_worker(input, output): for func, args in iter(input.get, 'STOP'): @@ -79,7 +83,7 @@ def contour_worker(input, output): for j in range(len(contour)): x.append(contour[j][0]) y.append(contour[j][1]) - TASKS.append((process_contour,(contour, image_data, fitsinfo, mean_beam))) + TASKS.append((process_contour,(contour, image_data, fitsinfo, noise_out))) task_queue = Queue() done_queue = Queue() # Submit tasks diff --git a/breizorro/utils.py b/breizorro/utils.py index 221584c..d23e0a8 100644 --- a/breizorro/utils.py +++ b/breizorro/utils.py @@ -225,7 +225,8 @@ def get_image_data(fitsfile): return image, header -def maxDist(contour, pixel_size): + +def maxDist(contour, pixel_size, x_centroid, y_centroid): """Calculate maximum extent and position angle of a contour. Parameters: @@ -233,6 +234,10 @@ def maxDist(contour, pixel_size): 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 @@ -242,35 +247,23 @@ def maxDist(contour, pixel_size): pos_angle : float Position angle of the contour (in degrees). """ - e_maj = 0 - e_min = np.inf # Start with a very large value for the minor axis - pos_angle = None - - # Convert the contour to a numpy array for easier calculations contour_array = np.array(contour) + distances = np.linalg.norm(contour_array - [y_centroid, x_centroid], axis=1) - # Step 1: Loop through all pairs of points to find major and minor axes - for i in range(len(contour_array)): - for j in range(i+1, len(contour_array)): - # Calculate Euclidean distance between points i and j - distance = np.linalg.norm(contour_array[i] - contour_array[j]) * pixel_size + # Find the indices of the points with maximum and minimum distances + max_idx = np.argmax(distances) + min_idx = np.argmin(distances) - # Calculate positional angle between points i and j - dx, dy = contour_array[j] - contour_array[i] - angle = np.degrees(np.arctan2(dy, dx)) + # Calculate major and minor axes + e_maj = np.max(distances) * pixel_size + e_min = np.min(distances) * pixel_size - # Update e_maj and pos_angle if this distance is the largest - if distance > e_maj: - e_maj = distance - pos_angle = angle - - # Update e_min if this distance is the smallest - if distance < e_min: - e_min = distance + # 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, @@ -294,8 +287,8 @@ def calculate_beam_area(bmaj, bmin, pix_size): return area -def get_source_size(contour, pixel_size, mean_beam, image, int_peak_ratio): - result = maxDist(contour,pixel_size) +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]) From 1fc13d9494330547ced8d1098e94f51d8546c972 Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Fri, 8 Nov 2024 17:22:42 +0200 Subject: [PATCH 33/48] add scikit-image to find contours --- pyproject.toml | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index ab32ee9..58e029c 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -41,9 +41,8 @@ requires = ["setuptools", "poetry-core"] build-backend = "poetry.core.masonry.api" [tool.poetry.extras] -all = ["bokeh", "photutils"] +all = ["bokeh", "photutils", "scikit-image"] gui = ["bokeh"] -catalog = ["photutils"] testing = ["pytest", "pytest-flake8"] [tool.poetry.group.tests] From aa79eb31fce7cee03748f245b73e402a85c332ef Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Fri, 8 Nov 2024 17:22:53 +0200 Subject: [PATCH 34/48] add scikit-image to find contours --- pyproject.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/pyproject.toml b/pyproject.toml index 58e029c..c96da35 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -35,6 +35,7 @@ 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"] From ead9bbb6eeeefb2d1f51145c2e3ea3de56a55372 Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Fri, 8 Nov 2024 18:23:49 +0200 Subject: [PATCH 35/48] Update index.md --- index.md | 54 ++++++++++++++++++++++++------------------------------ 1 file changed, 24 insertions(+), 30 deletions(-) diff --git a/index.md b/index.md index bb44aa2..df122df 100644 --- a/index.md +++ b/index.md @@ -74,7 +74,7 @@ breizorro -r circinus-MFS-image.fits --sum-peak 500 # Region Generation and Manipulation -Breizorro makes it easier to create and work with regions using image masks. It includes labeling, 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 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 --save-regions circinus.reg @@ -91,44 +91,38 @@ breizorro -r circinus-MFS-image.fits --merge west.reg --dilate 1 --fill-holes # Cataloguing and Visualization -Breizorro enables catalog generation from extracted regions, saving source properties to an ASCII/text -file. +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 circinus-MFS-image.fits --save-catalog circinus.txt +breizorro -r deep2-MFS-image.fits --save-catalog deep2.txt ``` ``` -# processing fits image: circinus-MFS-image.fits -# mean beam size (arcsec): 37.97 -# original image peak flux (Jy/beam): 0.768315315246582 -# noise out (μJy/beam): 737.33 -# cutt-off flux (mJy/beam): 4.79 -# freq0 (Hz): 1419944335.9375 -# number of sources detected: 35 +# 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 -147.88904620760084 -65.52043870236054 0.00468 0.0001 50.04 0.0 177.71 -src1 -147.84162114219356 -65.44040921224196 0.00047 0.0001 0.0 0.0 0.0 -src2 -147.7980562104931 -65.39402001087845 0.00089 0.0001 0.0 0.0 0.0 -src3 -147.74069959466678 -65.08592423258469 0.0082 0.0001 51.61 0.0 144.46 -src4 -147.6583486798732 -65.28753961408073 0.09922 0.0001 86.58 0.0 -173.37 -src5 -147.59829620974554 -65.28104296056243 0.00629 0.0001 47.07 0.0 167.74 -src6 -147.51788583319714 -65.72316001301358 0.06639 0.0001 71.47 0.0 162.07 -src7 -147.49179712006742 -64.88584736668952 0.00092 0.0001 0.0 0.0 0.0 -src8 -147.42970139783463 -65.65050232604096 0.00069 0.0001 0.0 0.0 0.0 -src9 -147.4293961973031 -65.44158144271519 0.00134 0.0001 0.0 0.0 0.0 -src10 -147.28370646739054 -65.42936506202037 9e-05 0.0001 0.0 0.0 0.0 +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 ``` -``` -breizorro -r circinus-MFS-image.fits --save-catalog circinus.txt --gui -``` - -![Screenshot 2024-09-11 100116](https://github.com/user-attachments/assets/79d3ece6-5d96-48a1-a06e-fbae2987d333) - -For a live demo of the viewer, visit: [breiz-viewer](breizorro.html) - # Contributors Thank you to the people who have contributed to this project. From e35df02fbf10dd54175ca34f316c685b89710b41 Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Fri, 8 Nov 2024 18:57:44 +0200 Subject: [PATCH 36/48] Optional deps --- breizorro/breizorro.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/breizorro/breizorro.py b/breizorro/breizorro.py index 05097eb..01a4cb5 100644 --- a/breizorro/breizorro.py +++ b/breizorro/breizorro.py @@ -24,9 +24,6 @@ import scipy.special import scipy.ndimage -from skimage.draw import polygon as skimage_polygon -from skimage.measure import find_contours - from breizorro.utils import get_source_size, format_source_coordinates, deg2ra, deg2dec from breizorro.utils import get_image_data, fitsInfo, calculate_beam_area @@ -353,6 +350,11 @@ def load_fits_or_region(filename): 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: From db804fd1d1be9f90af53267296985da0fd25019b Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Fri, 13 Dec 2024 21:25:26 +0200 Subject: [PATCH 37/48] Update deps --- .github/workflows/python-package.yml | 2 +- pyproject.toml | 4 +++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml index 4c89794..bdb4545 100644 --- a/.github/workflows/python-package.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 diff --git a/pyproject.toml b/pyproject.toml index c96da35..ac1105b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -24,11 +24,13 @@ packages = [ breizorro = "breizorro.main:driver" [tool.poetry.dependencies] -python = "^3.8" +python = ">=3.8, <3.13" astropy = "*" numpy = "*" +omegaconf = "*" regions = "*" scipy = "*" +stimela = ">=2.0" # Optional dependencies start here bokeh = { version = "*", optional = true} From 56ffc9687dafbeb3f30c1479eb6b1998f7a39689 Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Fri, 13 Dec 2024 21:27:00 +0200 Subject: [PATCH 38/48] use clickify_parameters --- breizorro/breizorro.py | 232 ++++++++++++++------------------ breizorro/config/breizorro.yaml | 120 +++++++++++++++++ breizorro/main.py | 12 +- index.md | 89 ++++++------ 4 files changed, 271 insertions(+), 182 deletions(-) create mode 100644 breizorro/config/breizorro.yaml diff --git a/breizorro/breizorro.py b/breizorro/breizorro.py index 01a4cb5..600c499 100644 --- a/breizorro/breizorro.py +++ b/breizorro/breizorro.py @@ -109,8 +109,10 @@ def remove_regions(mask_image, regs, wcs): 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 @@ -118,82 +120,44 @@ def main(): 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('-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.') - 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") + + # Example usage of the parameters + LOGGER.info(f"Restored Image: {restored_image}") + LOGGER.info(f"Mask Image: {mask_image}") + LOGGER.info(f"Threshold: {threshold}") + LOGGER.info(f"Box Size: {boxsize}") + LOGGER.info(f"Save Noise: {savenoise}") + LOGGER.info(f"Merge: {merge}") + LOGGER.info(f"Subtract: {subtract}") + LOGGER.info(f"Number Islands: {number_islands}") + LOGGER.info(f"Remove Islands: {remove_islands}") + LOGGER.info(f"Ignore Missing Islands: {ignore_missing_islands}") + LOGGER.info(f"Extract Islands: {extract_islands}") + LOGGER.info(f"Minimum Size: {minimum_size}") + LOGGER.info(f"Make Binary: {make_binary}") + LOGGER.info(f"Invert: {invert}") + LOGGER.info(f"Dilate: {dilate}") + LOGGER.info(f"Erode: {erode}") + LOGGER.info(f"Fill Holes: {fill_holes}") + LOGGER.info(f"Sum Peak: {sum_peak}") + LOGGER.info(f"NCPU: {ncpu}") + LOGGER.info(f"Beam Size: {beam_size}") + LOGGER.info(f"Outfile: {outfile}") + LOGGER.info(f"Save Catalog: {outcatalog}") + LOGGER.info(f"Save Regions: {outregion}") + + 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: + if restored_image: input_image, input_header = get_image_data(input_file) LOGGER.info(f"Generating mask using threshold {threshold}") @@ -213,15 +177,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_data(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) @@ -242,32 +206,31 @@ def load_fits_or_region(filename): raise 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() @@ -275,58 +238,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)) @@ -334,7 +297,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 @@ -349,7 +312,7 @@ def load_fits_or_region(filename): LOGGER.info("Done") sys.exit(1) - if args.outcatalog or args.outregion: + if outcatalog or outregion: try: from skimage.measure import find_contours except ImportError: @@ -367,27 +330,30 @@ def load_fits_or_region(filename): # 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 outregion: + regions.Regions(polygon_regions).write(outregion, format='ds9') + LOGGER.info(f"Saving regions in {outregion}") - if args.outcatalog and args.imagename: + 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] - warnings.resetwarnings() - warnings.filterwarnings('ignore', category=UserWarning, append=True) + # 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: - LOGGER.error("Running breizorro source detector requires optional dependencies, please re-install with: pip install breizorro[catalog]") - raise('Missing cataloguing dependencies') + 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(args.imagename) - fitsinfo = fitsInfo(args.imagename) + 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 args.beam: - mean_beam = float(args.beam) + if beam_size: + mean_beam = beam_size if mean_beam: LOGGER.info(f'Using user provided size: {mean_beam}') elif fitsinfo['b_size']: @@ -397,11 +363,9 @@ def load_fits_or_region(filename): 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 = open(outcatalog, 'w') + catalog_out = f'# processing fits image: {restored_image} \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' @@ -411,8 +375,8 @@ def load_fits_or_region(filename): 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) + 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' @@ -423,9 +387,9 @@ def load_fits_or_region(filename): output = 'src' + str(i) + ' ' + source_list[i][1] + '\n' f.write(output) f.close() - LOGGER.info(f'Source catalog saved: {args.outcatalog}') + LOGGER.info(f'Source catalog saved: {outcatalog}') - if args.gui: + if gui: try: from breizorro.gui import display except ModuleNotFoundError: @@ -433,7 +397,7 @@ def load_fits_or_region(filename): raise('Missing GUI dependencies') LOGGER.info("Loading Gui ...") - display(args.imagename or args.maskname, mask_image, args.outcatalog, source_list) + display(input_file, mask_image, outcatalog, source_list) LOGGER.info(f"Enforcing that mask to binary") mask_image = mask_image!=0 @@ -441,4 +405,4 @@ def load_fits_or_region(filename): shutil.copyfile(input_file, out_mask_fits) # to provide a template flush_fits(mask_image, out_mask_fits, mask_header) - LOGGER.info("Done") + LOGGER.info("Done") \ No newline at end of file diff --git a/breizorro/config/breizorro.yaml b/breizorro/config/breizorro.yaml new file mode 100644 index 0000000..310f89f --- /dev/null +++ b/breizorro/config/breizorro.yaml @@ -0,0 +1,120 @@ +opts: + log: + dir: logs/log-{config.run.datetime} + name: log-{info.taskname} + nest: 10 + symlink: log + +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 + required: false + mask_image: + info: "Input mask file(s). Either restored-image or mask-image must be specified." + dtype: File + required: false + threshold: + info: "Sigma threshold for masking (default = 6.5)" + dtype: float + required: false + default: 6.5 + boxsize: + info: "Box size over which to compute stats (default = 50)" + dtype: int + required: false + default: 50 + savenoise: + info: "Export noise image as FITS file" + dtype: bool + required: false + default: false + merge: + info: "Merge one or more masks or region files" + dtype: list[File] + required: false + subtract: + info: "Subtract one or more masks or region files" + dtype: list[File] + required: false + number_islands: + info: "Number the islands detected" + dtype: bool + required: false + default: false + remove_islands: + info: "Remove islands from input mask (list by number or coordinates)" + dtype: list[str] + required: false + ignore_missing_islands: + info: "Do not throw an error if an island specified by coordinates does not exist" + dtype: bool + required: false + default: false + extract_islands: + info: "Extract islands from input mask (list by number or coordinates)" + dtype: list[str] + required: false + minimum_size: + info: "Remove islands with areas fewer than or equal to the specified number of pixels" + dtype: int + required: false + make_binary: + info: "Replace all island numbers with 1" + dtype: bool + required: false + default: false + invert: + info: "Invert the mask" + dtype: bool + required: false + default: false + dilate: + info: "Apply dilation with a radius of R pixels" + dtype: int + required: false + default: 0 + erode: + info: "Apply N iterations of erosion" + dtype: int + required: false + default: 0 + fill_holes: + info: "Fill holes (closed regions) in the mask" + dtype: bool + required: false + default: false + sum_peak: + info: "Sum-to-peak ratio of flux islands to mask in original image" + dtype: float + required: false + ncpu: + info: "Number of processors to use for cataloging" + dtype: int + required: false + beam_size: + info: "Average beam size in arcsec if missing in the image header" + dtype: float + required: false + gui: + info: "Open mask in bokeh html gui" + dtype: bool + required: false + default: false + outputs: + outfile: + info: "Suffix for the mask image (default based on input name)" + dtype: File + required: false + outcatalog: + info: "Generate a catalog based on the region mask" + dtype: File + required: false + outregion: + info: "Generate polygon regions from the mask" + dtype: File + required: false diff --git a/breizorro/main.py b/breizorro/main.py index 58717c1..9bead95 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__), "config/breizorro.yaml")) + +@click.command("breizorro") +@clickify_parameters(schemas.cabs.get("breizorro")) +def driver(**kw): + main(**kw) diff --git a/index.md b/index.md index df122df..8ba624e 100644 --- a/index.md +++ b/index.md @@ -8,50 +8,47 @@ Breizorro is a flexible software program made to simplify image analysis tasks, # Parameter definition ``` -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. +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. ``` @@ -77,7 +74,7 @@ breizorro -r circinus-MFS-image.fits --sum-peak 500 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 --save-regions circinus.reg +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) @@ -96,7 +93,7 @@ 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 --save-catalog deep2.txt +breizorro -r deep2-MFS-image.fits --outcatalog deep2.txt ``` ``` From 19dbbec7440a3e566b479a4f033c3ea6ed43db67 Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Fri, 13 Dec 2024 21:37:28 +0200 Subject: [PATCH 39/48] Pin numpy to < 2 for older Python --- pyproject.toml | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index ac1105b..ce737a8 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -26,7 +26,10 @@ breizorro = "breizorro.main:driver" [tool.poetry.dependencies] python = ">=3.8, <3.13" astropy = "*" -numpy = "*" +numpy = [ + { version = ">=2.0.0", python = ">=3.12"}, + { version = "<2.0.0", python = ">=3.8, <3.12"} +] omegaconf = "*" regions = "*" scipy = "*" From 438e8bdff95fe2468d4130e699e6fcfdd0f6892f Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Fri, 13 Dec 2024 21:41:40 +0200 Subject: [PATCH 40/48] test help command --- .github/workflows/python-package.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml index bdb4545..e62236a 100644 --- a/.github/workflows/python-package.yml +++ b/.github/workflows/python-package.yml @@ -64,4 +64,4 @@ jobs: - name: Test command run: | - poetry run breizorro -h + poetry run breizorro --help From 7104796d1bb914a6b1ab769ad0fc15d776c59f63 Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Fri, 13 Dec 2024 21:49:03 +0200 Subject: [PATCH 41/48] update deps --- pyproject.toml | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index ce737a8..6ada5d9 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -25,14 +25,20 @@ breizorro = "breizorro.main:driver" [tool.poetry.dependencies] python = ">=3.8, <3.13" -astropy = "*" +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 From 53fef7e9a9d26a4023c20817921d1d055b77e038 Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Fri, 13 Dec 2024 21:56:32 +0200 Subject: [PATCH 42/48] exit 0 --- breizorro/breizorro.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/breizorro/breizorro.py b/breizorro/breizorro.py index 600c499..a96255b 100644 --- a/breizorro/breizorro.py +++ b/breizorro/breizorro.py @@ -310,7 +310,7 @@ def load_fits_or_region(filename): 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) + sys.exit(0) if outcatalog or outregion: try: From 285d453ce52374b7db9929abce2511f38f976f68 Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Sat, 14 Dec 2024 15:43:22 +0200 Subject: [PATCH 43/48] Upate config structure --- MANIFEST.in | 1 + breizorro/{config => }/breizorro.yaml | 53 +++++++++------------------ 2 files changed, 19 insertions(+), 35 deletions(-) rename breizorro/{config => }/breizorro.yaml (77%) 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/breizorro/config/breizorro.yaml b/breizorro/breizorro.yaml similarity index 77% rename from breizorro/config/breizorro.yaml rename to breizorro/breizorro.yaml index 310f89f..65f30de 100644 --- a/breizorro/config/breizorro.yaml +++ b/breizorro/breizorro.yaml @@ -1,10 +1,3 @@ -opts: - log: - dir: logs/log-{config.run.datetime} - name: log-{info.taskname} - nest: 10 - symlink: log - cabs: breizorro: info: "Mask creation and manipulation for radio astronomy images (https://github.com/ratt-ru/breizorro)" @@ -13,108 +6,98 @@ cabs: restored_image: info: "Restored image file from which to build the mask" dtype: File - required: false + abbreviation: r mask_image: info: "Input mask file(s). Either restored-image or mask-image must be specified." dtype: File - required: false + abbreviation: m threshold: info: "Sigma threshold for masking (default = 6.5)" dtype: float - required: false default: 6.5 + abbreviation: t boxsize: info: "Box size over which to compute stats (default = 50)" dtype: int - required: false default: 50 + abbreviation: b savenoise: info: "Export noise image as FITS file" dtype: bool - required: false default: false merge: info: "Merge one or more masks or region files" - dtype: list[File] - required: false + dtype: List[File] + policies: + repeat: ' ' subtract: info: "Subtract one or more masks or region files" - dtype: list[File] - required: false + dtype: List[File] + policies: + repeat: list number_islands: info: "Number the islands detected" dtype: bool - required: false default: false remove_islands: info: "Remove islands from input mask (list by number or coordinates)" - dtype: list[str] - required: false + dtype: List[str] + policies: + repeat: list ignore_missing_islands: info: "Do not throw an error if an island specified by coordinates does not exist" dtype: bool - required: false default: false extract_islands: info: "Extract islands from input mask (list by number or coordinates)" - dtype: list[str] - required: false + dtype: List[str] + policies: + repeat: list minimum_size: info: "Remove islands with areas fewer than or equal to the specified number of pixels" dtype: int - required: false make_binary: info: "Replace all island numbers with 1" dtype: bool - required: false default: false invert: info: "Invert the mask" dtype: bool - required: false default: false dilate: info: "Apply dilation with a radius of R pixels" dtype: int - required: false default: 0 erode: info: "Apply N iterations of erosion" dtype: int - required: false default: 0 fill_holes: info: "Fill holes (closed regions) in the mask" dtype: bool - required: false default: false sum_peak: info: "Sum-to-peak ratio of flux islands to mask in original image" dtype: float - required: false ncpu: info: "Number of processors to use for cataloging" dtype: int - required: false + abbreviation: j beam_size: info: "Average beam size in arcsec if missing in the image header" dtype: float - required: false gui: info: "Open mask in bokeh html gui" dtype: bool - required: false default: false outputs: outfile: info: "Suffix for the mask image (default based on input name)" dtype: File - required: false + abbreviation: o outcatalog: info: "Generate a catalog based on the region mask" dtype: File - required: false outregion: info: "Generate polygon regions from the mask" dtype: File - required: false From d32db93948f30055ed2907f3cf9f3780e3127205 Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Sat, 14 Dec 2024 15:44:25 +0200 Subject: [PATCH 44/48] Clean up --- breizorro/breizorro.py | 33 ++++----------------------------- breizorro/main.py | 2 +- 2 files changed, 5 insertions(+), 30 deletions(-) diff --git a/breizorro/breizorro.py b/breizorro/breizorro.py index a96255b..3d64d22 100644 --- a/breizorro/breizorro.py +++ b/breizorro/breizorro.py @@ -118,34 +118,8 @@ def main(restored_image, mask_image, threshold, boxsize, savenoise, merge, subtr from pkg_resources import get_distribution _version = get_distribution('breizorro').version LOGGER.info(f"Version: {_version}") - LOGGER.info("Usage: breizorro --help") - # Example usage of the parameters - LOGGER.info(f"Restored Image: {restored_image}") - LOGGER.info(f"Mask Image: {mask_image}") - LOGGER.info(f"Threshold: {threshold}") - LOGGER.info(f"Box Size: {boxsize}") - LOGGER.info(f"Save Noise: {savenoise}") - LOGGER.info(f"Merge: {merge}") - LOGGER.info(f"Subtract: {subtract}") - LOGGER.info(f"Number Islands: {number_islands}") - LOGGER.info(f"Remove Islands: {remove_islands}") - LOGGER.info(f"Ignore Missing Islands: {ignore_missing_islands}") - LOGGER.info(f"Extract Islands: {extract_islands}") - LOGGER.info(f"Minimum Size: {minimum_size}") - LOGGER.info(f"Make Binary: {make_binary}") - LOGGER.info(f"Invert: {invert}") - LOGGER.info(f"Dilate: {dilate}") - LOGGER.info(f"Erode: {erode}") - LOGGER.info(f"Fill Holes: {fill_holes}") - LOGGER.info(f"Sum Peak: {sum_peak}") - LOGGER.info(f"NCPU: {ncpu}") - LOGGER.info(f"Beam Size: {beam_size}") - LOGGER.info(f"Outfile: {outfile}") - LOGGER.info(f"Save Catalog: {outcatalog}") - LOGGER.info(f"Save Regions: {outregion}") - 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: @@ -202,8 +176,9 @@ def load_fits_or_region(filename): 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 merge: @@ -405,4 +380,4 @@ def load_fits_or_region(filename): shutil.copyfile(input_file, out_mask_fits) # to provide a template flush_fits(mask_image, out_mask_fits, mask_header) - LOGGER.info("Done") \ No newline at end of file + LOGGER.info("Done") diff --git a/breizorro/main.py b/breizorro/main.py index 9bead95..8047838 100644 --- a/breizorro/main.py +++ b/breizorro/main.py @@ -4,7 +4,7 @@ from omegaconf import OmegaConf from breizorro.breizorro import main -schemas = OmegaConf.load(os.path.join(os.path.dirname(__file__), "config/breizorro.yaml")) +schemas = OmegaConf.load(os.path.join(os.path.dirname(__file__), "breizorro.yaml")) @click.command("breizorro") @clickify_parameters(schemas.cabs.get("breizorro")) From 3d2685285a8d7da610f75b627c34679bee76f5d0 Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Sat, 14 Dec 2024 23:43:02 +0200 Subject: [PATCH 45/48] add metavar --- breizorro/breizorro.yaml | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/breizorro/breizorro.yaml b/breizorro/breizorro.yaml index 65f30de..bc0514d 100644 --- a/breizorro/breizorro.yaml +++ b/breizorro/breizorro.yaml @@ -6,11 +6,13 @@ cabs: 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 @@ -28,16 +30,19 @@ cabs: 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: list number_islands: info: "Number the islands detected" dtype: bool + metavar: "N|COORD" default: false remove_islands: info: "Remove islands from input mask (list by number or coordinates)" @@ -51,6 +56,7 @@ cabs: extract_islands: info: "Extract islands from input mask (list by number or coordinates)" dtype: List[str] + metavar: "N|COORD" policies: repeat: list minimum_size: @@ -67,10 +73,12 @@ cabs: 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" From 866e08d812ccb5cdb3ecc37c0098fa0d6f07c265 Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Thu, 19 Dec 2024 00:26:40 +0200 Subject: [PATCH 46/48] Update repeat policies --- breizorro/breizorro.py | 6 +++--- breizorro/breizorro.yaml | 8 ++++---- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/breizorro/breizorro.py b/breizorro/breizorro.py index 3d64d22..48cc97b 100644 --- a/breizorro/breizorro.py +++ b/breizorro/breizorro.py @@ -374,9 +374,9 @@ def load_fits_or_region(filename): LOGGER.info("Loading Gui ...") display(input_file, mask_image, outcatalog, source_list) - LOGGER.info(f"Enforcing that mask to binary") - mask_image = mask_image!=0 - mask_header['BUNIT'] = 'mask' + 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) diff --git a/breizorro/breizorro.yaml b/breizorro/breizorro.yaml index bc0514d..26ed1ad 100644 --- a/breizorro/breizorro.yaml +++ b/breizorro/breizorro.yaml @@ -32,13 +32,13 @@ cabs: dtype: List[File] metavar: "MASK(s)|REG(s)" policies: - repeat: ' ' + repeat: '[]' subtract: info: "Subtract one or more masks or region files" dtype: List[File] metavar: "MASK(s)|REG(s)" policies: - repeat: list + repeat: '[]' number_islands: info: "Number the islands detected" dtype: bool @@ -48,7 +48,7 @@ cabs: info: "Remove islands from input mask (list by number or coordinates)" dtype: List[str] policies: - repeat: list + repeat: ' ' ignore_missing_islands: info: "Do not throw an error if an island specified by coordinates does not exist" dtype: bool @@ -58,7 +58,7 @@ cabs: dtype: List[str] metavar: "N|COORD" policies: - repeat: list + repeat: ' ' minimum_size: info: "Remove islands with areas fewer than or equal to the specified number of pixels" dtype: int From 11d9211dc832ff44c62b33b677b0f37674fed273 Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Fri, 20 Dec 2024 00:58:49 +0200 Subject: [PATCH 47/48] Update dtype and add examples --- breizorro/breizorro.py | 4 ++-- breizorro/breizorro.yaml | 12 ++++++------ 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/breizorro/breizorro.py b/breizorro/breizorro.py index 48cc97b..32f0aa4 100644 --- a/breizorro/breizorro.py +++ b/breizorro/breizorro.py @@ -80,9 +80,9 @@ def make_noise_map(restored_image, boxsize): 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)) diff --git a/breizorro/breizorro.yaml b/breizorro/breizorro.yaml index 26ed1ad..34b55a0 100644 --- a/breizorro/breizorro.yaml +++ b/breizorro/breizorro.yaml @@ -42,20 +42,20 @@ cabs: number_islands: info: "Number the islands detected" dtype: bool - metavar: "N|COORD" default: false remove_islands: - info: "Remove islands from input mask (list by number or coordinates)" - dtype: List[str] + 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: ' ' + 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)" - dtype: List[str] + 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: ' ' From 4d0d6092cde6a8ff78c16e6e71d861f4c1648f40 Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Mon, 13 Jan 2025 16:27:46 +0200 Subject: [PATCH 48/48] Use rather than --- breizorro/breizorro.yaml | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/breizorro/breizorro.yaml b/breizorro/breizorro.yaml index 34b55a0..ab9a7d4 100644 --- a/breizorro/breizorro.yaml +++ b/breizorro/breizorro.yaml @@ -3,12 +3,12 @@ cabs: info: "Mask creation and manipulation for radio astronomy images (https://github.com/ratt-ru/breizorro)" command: breizorro inputs: - restored_image: + restored-image: info: "Restored image file from which to build the mask" dtype: File metavar: IMAGE abbreviation: r - mask_image: + mask-image: info: "Input mask file(s). Either restored-image or mask-image must be specified." dtype: File abbreviation: m @@ -39,30 +39,30 @@ cabs: metavar: "MASK(s)|REG(s)" policies: repeat: '[]' - number_islands: + number-islands: info: "Number the islands detected" dtype: bool default: false - remove_islands: + 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: + ignore-missing-islands: info: "Do not throw an error if an island specified by coordinates does not exist" dtype: bool default: false - extract_islands: + 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: + minimum-size: info: "Remove islands with areas fewer than or equal to the specified number of pixels" dtype: int - make_binary: + make-binary: info: "Replace all island numbers with 1" dtype: bool default: false @@ -80,18 +80,18 @@ cabs: dtype: int metavar: "N" default: 0 - fill_holes: + fill-holes: info: "Fill holes (closed regions) in the mask" dtype: bool default: false - sum_peak: + 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: + beam-size: info: "Average beam size in arcsec if missing in the image header" dtype: float gui: