From b162a2714584370fcf843d43cc0d6ac33d8c04cf Mon Sep 17 00:00:00 2001 From: tgalvin Date: Mon, 7 Oct 2024 17:24:48 +0800 Subject: [PATCH 01/11] Beginning of a toy --- flint/masking.py | 41 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 41 insertions(+) diff --git a/flint/masking.py b/flint/masking.py index 7008f6e2..bea25ba4 100644 --- a/flint/masking.py +++ b/flint/masking.py @@ -8,6 +8,7 @@ from pathlib import Path from typing import Collection, NamedTuple, Optional, Union +import astropy.units as u import numpy as np from astropy.io import fits from astropy.wcs import WCS @@ -19,6 +20,7 @@ binary_erosion as scipy_binary_erosion, # Rename to distinguish from skimage ) from scipy.ndimage import label, minimum_filter +from radio_beam import Beam from flint.logging import logger from flint.naming import FITSMaskNames, create_fits_mask_names @@ -107,6 +109,45 @@ def consider_beam_mask_round( ) # type: ignore +def make_beam_msk_kernel( + fits_header: fits.Header, kernel_size=100, minimum_response: float = 0.6 +) -> np.ndarray: + """Make a mask using the shape of a beam in a FITS Header object. The + beam properties in the ehader are used to generate the two-dimensional + Gaussian main lobe, from which a cut is made based on the minimum + power. + + Args: + fits_header (fits.Header): The FITS header to create beam from + kernel_size (int, optional): Size of the output kernel in pixels. Will be a square. Defaults to 100. + minimum_response (float, optional): Minimum response of the beam shape for the mask to be constructed from. Defaults to 0.6. + + Raises: + KeyError: Raised if CDELT1 and CDELT2 missing + + Returns: + np.ndarray: Boolean mask of the kernel shape + """ + + POSITION_KEYS = ("CDELT1", "CDELT2") + if all([key in fits_header for key in POSITION_KEYS]): + raise KeyError(f"{POSITION_KEYS=} all need to be present") + + beam = Beam.from_fits_header(fits_header) + assert isinstance(beam, Beam) + + cdelt1, cdelt2 = np.abs(fits_header["CDELT1"]), np.abs(fits_header["CDELT2"]) # type: ignore + assert np.isclose( + cdelt1, cdelt2 + ), f"Pixel scales {cdelt1=} {cdelt2=}, but must be equal" + + k = beam.as_kernel( + pixscale=cdelt1 * u.Unit("deg"), x_size=kernel_size, y_size=kernel_size + ) + + return k.array > (np.max(k.array) * minimum_response) + + def extract_beam_mask_from_mosaic( fits_beam_image_path: Path, fits_mosaic_mask_names: FITSMaskNames ) -> FITSMaskNames: From 3b1386054e7df287fb50f3bd9076f639ab3911e8 Mon Sep 17 00:00:00 2001 From: tgalvin Date: Mon, 7 Oct 2024 17:26:05 +0800 Subject: [PATCH 02/11] branch and typo fix --- flint/masking.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/flint/masking.py b/flint/masking.py index bea25ba4..49fc3f4a 100644 --- a/flint/masking.py +++ b/flint/masking.py @@ -109,7 +109,7 @@ def consider_beam_mask_round( ) # type: ignore -def make_beam_msk_kernel( +def make_beam_mask_kernel( fits_header: fits.Header, kernel_size=100, minimum_response: float = 0.6 ) -> np.ndarray: """Make a mask using the shape of a beam in a FITS Header object. The From 3ffabb50b92677959c75e6bb06ce941c304468a5 Mon Sep 17 00:00:00 2001 From: tgalvin Date: Mon, 7 Oct 2024 22:07:06 +0800 Subject: [PATCH 03/11] added initial toy beam shape erode mask / tests --- flint/masking.py | 46 +++++++++++++++++++++++++-- tests/test_masking.py | 73 +++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 117 insertions(+), 2 deletions(-) diff --git a/flint/masking.py b/flint/masking.py index 49fc3f4a..3ee7c991 100644 --- a/flint/masking.py +++ b/flint/masking.py @@ -65,6 +65,10 @@ class MaskingOptions(NamedTuple): """Size of the boxcar filter""" minimum_boxcar_increase_factor: float = 1.2 """The factor used to construct minimum positive signal threshold for an island """ + beam_shape_erode: bool = False + """Erode the mask using the shape of the restoring beam""" + beam_shape_erode_minimum_response: float = 0.6 + """The minimum response of the beam that is used to form t he erode structure shape""" def with_options(self, **kwargs) -> MaskingOptions: """Return a new instance of the MaskingOptions""" @@ -109,7 +113,7 @@ def consider_beam_mask_round( ) # type: ignore -def make_beam_mask_kernel( +def create_beam_mask_kernel( fits_header: fits.Header, kernel_size=100, minimum_response: float = 0.6 ) -> np.ndarray: """Make a mask using the shape of a beam in a FITS Header object. The @@ -128,9 +132,10 @@ def make_beam_mask_kernel( Returns: np.ndarray: Boolean mask of the kernel shape """ + assert minimum_response > 0, f"{minimum_response=}, should be positive" POSITION_KEYS = ("CDELT1", "CDELT2") - if all([key in fits_header for key in POSITION_KEYS]): + if not all([key in fits_header for key in POSITION_KEYS]): raise KeyError(f"{POSITION_KEYS=} all need to be present") beam = Beam.from_fits_header(fits_header) @@ -148,6 +153,43 @@ def make_beam_mask_kernel( return k.array > (np.max(k.array) * minimum_response) +def beam_shape_erode( + mask: np.ndarray, fits_header: fits.Header, minimum_response: float = 0.6 +) -> np.ndarray: + """Construct a kernel representing the shape of the restoring beam at + a particular level, and use it as the basis of a binary erosion of the + input mask. + + The ``fits_header`` is used to construct the beam shape that matches the + same pixel size + + Args: + mask (np.ndarray): The current mask that will be eroded based on the beam shape + fits_header (fits.Header): The fits header of the mask used to generate the beam kernel shape + minimum_response (float, optional): The minimum response of the main restoring beam to craft the shape from. Defaults to 0.6. + + Returns: + np.ndarray: The eroded beam shape + """ + + if not all([key in fits_header for key in ["BMAJ", "BMIN", "BPA"]]): + logger.warning( + "Beam parameters missing. Not performing the beam shape erosion. " + ) + return mask + + logger.info("Eroding the mask using the beam shape") + beam_mask_kernel = create_beam_mask_kernel( + fits_header=fits_header, minimum_response=minimum_response + ) + + erode_mask = scipy_binary_erosion( + input=mask, iterations=1, structure=beam_mask_kernel + ) + + return erode_mask + + def extract_beam_mask_from_mosaic( fits_beam_image_path: Path, fits_mosaic_mask_names: FITSMaskNames ) -> FITSMaskNames: diff --git a/tests/test_masking.py b/tests/test_masking.py index d5991037..4beb54e1 100644 --- a/tests/test_masking.py +++ b/tests/test_masking.py @@ -7,7 +7,9 @@ from flint.masking import ( MaskingOptions, _verify_set_positive_seed_clip, + beam_shape_erode, consider_beam_mask_round, + create_beam_mask_kernel, create_snr_mask_from_fits, minimum_boxcar_artefact_mask, ) @@ -16,6 +18,77 @@ SHAPE = (100, 100) +def test_create_beam_mask_kernel(): + """See whether the beam kernel creation mask can return a correct mask""" + fits_header = fits.Header( + dict( + CDELT1=-0.000694444444444444, + CDELT2=0.000694444444444444, + BMAJ=0.00340540107886635, + BMIN=0.00283268735470751, + BPA=74.6618858613889, + ) + ) + mask_1 = create_beam_mask_kernel(fits_header=fits_header) + assert mask_1.shape == (100, 100) + assert np.sum(mask_1) == 12 + + mask_2 = create_beam_mask_kernel(fits_header=fits_header, minimum_response=0.1) + assert mask_2.shape == (100, 100) + assert np.sum(mask_2) == 52 + + with pytest.raises(AssertionError): + fits_header = fits.Header( + dict( + CDELT1=-0.000694444444444444, + CDELT2=0.2000694444444444444, + BMAJ=0.00340540107886635, + BMIN=0.00283268735470751, + BPA=74.6618858613889, + ) + ) + create_beam_mask_kernel(fits_header=fits_header) + + with pytest.raises(KeyError): + fits_header = fits.Header( + dict( + CDELT1=-0.000694444444444444, + BMAJ=0.00340540107886635, + BMIN=0.00283268735470751, + BPA=74.6618858613889, + ) + ) + create_beam_mask_kernel(fits_header=fits_header) + + +def test_beam_shape_erode(): + """Ensure that the beam shape erosion approach works. This should drop out pixels + should the beam shape structure connectivity be statisifed""" + fits_header = fits.Header( + dict( + CDELT1=-0.000694444444444444, + CDELT2=0.000694444444444444, + BMAJ=0.00340540107886635, + BMIN=0.00283268735470751, + BPA=74.6618858613889, + ) + ) + + mask = np.zeros((500, 500)).astype(bool) + + mask[300, 300] = True + assert np.sum(mask) == 1 + + new_mask = beam_shape_erode(mask=mask, fits_header=fits_header) + assert np.sum(new_mask) == 0 + + mask = np.zeros((200, 200)).astype(bool) + mask[100:130, 100:130] = True + assert np.sum(mask) == 900 + erode_mask = beam_shape_erode(mask=mask, fits_header=fits_header) + assert np.sum(erode_mask) == 729 + + def test_consider_beam_masking_round(): """Test to ensure the beam mask consideration log is correct""" lower = ("all", "ALL", "aLl") From 81ce006f8f2b3cb68f944b11b94321af7d2b6556 Mon Sep 17 00:00:00 2001 From: tgalvin Date: Mon, 7 Oct 2024 22:09:44 +0800 Subject: [PATCH 04/11] another test --- tests/test_masking.py | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/tests/test_masking.py b/tests/test_masking.py index 4beb54e1..9cd1e53f 100644 --- a/tests/test_masking.py +++ b/tests/test_masking.py @@ -89,6 +89,26 @@ def test_beam_shape_erode(): assert np.sum(erode_mask) == 729 +def test_beam_shape_erode_nobeam(): + """Ensure that the beam shape erosion approach works. This should drop out pixels + should the beam shape structure connectivity be statisifed. This should simply return + the input array since there is no beam information""" + fits_header = fits.Header( + dict( + CDELT1=-0.000694444444444444, + CDELT2=0.000694444444444444, + ) + ) + + mask = np.zeros((500, 500)).astype(bool) + + mask[300, 300] = True + assert np.sum(mask) == 1 + new_mask = beam_shape_erode(mask=mask, fits_header=fits_header) + assert new_mask is mask + assert np.sum(new_mask) == 1 + + def test_consider_beam_masking_round(): """Test to ensure the beam mask consideration log is correct""" lower = ("all", "ALL", "aLl") From d295ff0e1b944ad429160aa4c4478c4f98345506 Mon Sep 17 00:00:00 2001 From: tgalvin Date: Tue, 8 Oct 2024 16:36:16 +0800 Subject: [PATCH 05/11] cleaning up CLI / added some tests --- flint/masking.py | 99 ++++++++++++++++++++++++++++++++++--------- tests/test_masking.py | 30 +++++++++++-- 2 files changed, 106 insertions(+), 23 deletions(-) diff --git a/flint/masking.py b/flint/masking.py index 3ee7c991..84d116d1 100644 --- a/flint/masking.py +++ b/flint/masking.py @@ -4,7 +4,7 @@ from __future__ import annotations -from argparse import ArgumentParser +from argparse import ArgumentParser, Namespace from pathlib import Path from typing import Collection, NamedTuple, Optional, Union @@ -26,6 +26,8 @@ from flint.naming import FITSMaskNames, create_fits_mask_names from flint.utils import get_pixels_per_beam +# TODO: Need to remove a fair amount of old approaches, and deprecate some of the toy functions + class MaskingOptions(NamedTuple): """Contains options for the creation of clean masks from some subject @@ -36,16 +38,16 @@ class MaskingOptions(NamedTuple): base_snr_clip: float = 4 """A base clipping level to be used should other options not be activated""" flood_fill: bool = True - """Whether to attempt to flood fill when constructing a mask. This should be `True` for `grow_low_snr_islands` and `suppress_artefacts to have an effect. """ + """Whether to attempt to flood fill when constructing a mask. This should be `True` for ``grow_low_snr_islands`` and ``suppress_artefacts`` to have an effect. """ flood_fill_positive_seed_clip: float = 4.5 - """The clipping level to seed islands that will be grown to lower SNR""" + """The clipping level to seed islands that will be grown to lower signal metric""" flood_fill_positive_flood_clip: float = 1.5 """Clipping level used to grow seeded islands down to""" flood_fill_use_mbc: bool = False """If True, the clipping levels are used as the `increase_factor` when using a minimum absolute clip""" flood_fill_use_mbc_box_size: int = 75 """The size of the mbc box size should mbc be used""" - suppress_artefacts: bool = True + suppress_artefacts: bool = False """Whether to attempt artefacts based on the presence of significant negatives""" suppress_artefacts_negative_seed_clip: float = 5 """The significance level of a negative island for the sidelobe suppression to be activated. This should be a positive number (the signal map is internally inverted)""" @@ -59,7 +61,7 @@ class MaskingOptions(NamedTuple): """The minimum significance levels of pixels to be to seed low SNR islands for consideration""" grow_low_snr_island_size: int = 768 """The number of pixels an island has to be for it to be accepted""" - minimum_boxcar: bool = True + minimum_boxcar: bool = False """Use the boxcar minimum threshold to compare to remove artefacts""" minimum_boxcar_size: int = 100 """Size of the boxcar filter""" @@ -346,6 +348,7 @@ def minimum_boxcar_artefact_mask( island region is larger to the maximum signal then that island is omitted. + This is planned to be deprecated. Args: signal (np.ndarray): The input signl to use @@ -356,6 +359,10 @@ def minimum_boxcar_artefact_mask( Returns: np.ndarray: _description_ """ + import warnings + + warnings.warn("minimum_boxcar_artefacts to be removed. ", DeprecationWarning) + return island_mask logger.info( f"Running boxcar minimum island clip with {boxcar_size=} {increase_factor=}" @@ -407,7 +414,7 @@ def suppress_artefact_mask( pixels_per_beam: Optional[float] = None, large_island_threshold: float = 1.0, ) -> np.ndarray: - """Attempt to grow mask that sepresses artefacts around bright sources. Small islands + """Attempt to grow mask that suppresses artefacts around bright sources. Small islands of negative emission seed pixels, which then grow out. Bright positive pixels are not allowed to change (which presumably are the source of negative artetfacts). @@ -419,6 +426,8 @@ def suppress_artefact_mask( For this reason the guard mask should be sufficiently high to protect the main source but nuke the fask positive islands + Consider using ``minimum_absolute_clip``, which has a simpler set of options and has + largely the same intended result. Args: signal (np.ndarray): The signal mask, @@ -662,14 +671,14 @@ def create_snr_mask_from_fits( ) with fits.open(fits_image_path) as fits_image: - fits_header = fits_image[0].header + fits_header = fits_image[0].header # type: ignore with fits.open(fits_bkg_path) as fits_bkg: logger.info("Subtracting background") - signal_data = fits_image[0].data - fits_bkg[0].data + signal_data = fits_image[0].data - fits_bkg[0].data # type: ignore with fits.open(fits_rms_path) as fits_rms: logger.info("Dividing by RMS") - signal_data /= fits_rms[0].data + signal_data /= fits_rms[0].data # type: ignore if create_signal_fits: logger.info(f"Writing {mask_names.signal_fits}") @@ -701,6 +710,13 @@ def create_snr_mask_from_fits( logger.info(f"Clipping using a {masking_options.base_snr_clip=}") mask_data = (signal_data > masking_options.base_snr_clip).astype(np.int32) + if masking_options.beam_shape_erode: + mask_data = beam_shape_erode( + mask=mask_data, + fits_header=fits_header, + minimum_response=masking_options.beam_shape_erode_minimum_response, + ) + logger.info(f"Writing {mask_names.mask_fits}") fits.writeto( filename=mask_names.mask_fits, @@ -722,7 +738,7 @@ def get_parser() -> ArgumentParser: ) fits_parser = subparser.add_parser( - "snrmask", + "mask", help="Create a mask for an image, using its RMS and BKG images (e.g. outputs from BANE). Output FITS image will default to the image with a mask suffix.", ) fits_parser.add_argument("image", type=Path, help="Path to the input image. ") @@ -740,22 +756,52 @@ def get_parser() -> ArgumentParser: help="Save the signal map. Defaults to the same as image with a signal suffix. ", ) fits_parser.add_argument( - "--min-snr", + "--base-snr-clip", type=float, default=4, - help="The minimum SNR required for a pixel to be marked as valid. ", + help="A base clipping level to be used should other options not be activated", + ) + fits_parser.add_argument( + "--flood-fill", + action="store_true", + default=False, + help="Whether to attempt to flood fill when constructing a mask. This should be `True` for `grow_low_snr_islands` and `suppress_artefacts to have an effect. ", + ) + fits_parser.add_argument( + "--flood-fill-positive-seed-clip", + type=float, + default=4.5, + help="The clipping level to seed islands that will be grown to lower signal metric", + ) + fits_parser.add_argument( + "--flood-fill-positive-flood-clip", + type=float, + default=1.5, + help="Clipping level used to grow seeded islands down to", ) fits_parser.add_argument( - "--use-butterworth", + "--flood-fill-use-mbc", action="store_true", - help="Apply a butterworth filter to smooth the total intensity image before computing the signal map. ", + default=False, + help="If True, the clipping levels are used as the `increase_factor` when using a minimum absolute clip. ", ) fits_parser.add_argument( - "--connectivity-shape", - default=(4, 4), - nargs=2, + "--flood-fill-use-mbc-box-size", type=int, - help="The connectivity matrix to use when applying a binary erosion. Only used when using the butterworth smoothing filter. ", + default=75, + help="The size of the mbc box size should mbc be used", + ) + fits_parser.add_argument( + "--beam-shape-erode", + action="store_true", + default=False, + help="Erode the mask using the shape of the restoring beam", + ) + fits_parser.add_argument( + "--beam-shape-erode-minimum-response", + type=float, + default=0.6, + help="The minimum response of the beam that is used to form t he erode structure shape", ) extract_parser = subparser.add_parser( @@ -778,13 +824,28 @@ def get_parser() -> ArgumentParser: return parser +def _args_to_mask_options(args: Namespace) -> MaskingOptions: + """Convert the args namespace to a MaskingOptions""" + masking_options = MaskingOptions( + base_snr_clip=args.base_snr_clip, + flood_fill=args.flood_fill, + flood_fill_positive_seed_clip=args.flood_fill_positive_seed_clip, + flood_fill_positive_flood_clip=args.flood_fill_positive_flood_clip, + flood_fill_use_mbc=args.flood_fill_use_mbc, + flood_fill_use_mbc_box_size=args.flood_fill_use_mbc_box_size, + beam_shape_erode=args.beam_shape_erode, + beam_shape_erode_minimum_response=args.beam_shape_erode_minimum_response, + ) + return masking_options + + def cli(): parser = get_parser() args = parser.parse_args() if args.mode == "snrmask": - masking_options = MaskingOptions(base_snr_clip=args.min_snr) + masking_options = _args_to_mask_options(args=args) create_snr_mask_from_fits( fits_image_path=args.image, fits_rms_path=args.rms, diff --git a/tests/test_masking.py b/tests/test_masking.py index 9cd1e53f..968ccdec 100644 --- a/tests/test_masking.py +++ b/tests/test_masking.py @@ -6,11 +6,13 @@ from flint.masking import ( MaskingOptions, + _args_to_mask_options, _verify_set_positive_seed_clip, beam_shape_erode, consider_beam_mask_round, create_beam_mask_kernel, create_snr_mask_from_fits, + get_parser, minimum_boxcar_artefact_mask, ) from flint.naming import FITSMaskNames @@ -18,6 +20,21 @@ SHAPE = (100, 100) +def test_arg_parser_cli_and_masking_options(): + """See if the CLI parser is constructed with correct set of + options which are properly converted to a MaskingOptions object""" + parser = get_parser() + args = parser.parse_args( + args="mask img rms bkg --flood-fill --flood-fill-positive-seed-clip 10 --flood-fill-positive-flood-clip 1. --flood-fill-use-mbc --flood-fill-use-mbc-box-size 100".split() + ) + masking_options = _args_to_mask_options(args=args) + assert isinstance(masking_options, MaskingOptions) + assert masking_options.flood_fill + assert masking_options.flood_fill_use_mbc + assert masking_options.flood_fill_positive_seed_clip == 10.0 + assert masking_options.flood_fill_positive_flood_clip == 1.0 + + def test_create_beam_mask_kernel(): """See whether the beam kernel creation mask can return a correct mask""" fits_header = fits.Header( @@ -148,13 +165,16 @@ def test_minimum_boxcar_artefact(): signal=img, island_mask=img_mask, boxcar_size=10 ) assert np.all(img_mask == out_mask) - assert img_mask is not out_mask + assert out_mask is img_mask # minimum boxcar artefact mask to be deprecated + # assert img_mask is not out_mask img[41:45, 30:40] = -20 out_mask = minimum_boxcar_artefact_mask( signal=img, island_mask=img_mask, boxcar_size=10 ) - assert not np.all(img_mask == out_mask) + assert out_mask is img_mask # minimum boxcar artefact mask to be deprecated + + # assert not np.all(img_mask == out_mask) def test_minimum_boxcar_artefact_blanked(): @@ -171,8 +191,10 @@ def test_minimum_boxcar_artefact_blanked(): out_mask = minimum_boxcar_artefact_mask( signal=img, island_mask=img_mask, boxcar_size=10, increase_factor=1000 ) - assert out_mask is not img_mask - assert np.all(out_mask[30:40, 30:40] == False) # noqa + assert out_mask is img_mask # minimum boxcar artefact mask will be deprecated + + # assert out_mask is not img_mask + # assert np.all(out_mask[30:40, 30:40] == False) # noqa def test_minimum_boxcar_large_bright_island(): From fc7ead23659188f4535202eb882020c900b1ad88 Mon Sep 17 00:00:00 2001 From: tgalvin Date: Tue, 8 Oct 2024 16:39:02 +0800 Subject: [PATCH 06/11] added notes to change log --- CHANGELOG.md | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 16f07f6b..30792f40 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,11 @@ - cubes are supported when computing the yandasoft linmos weights and trimming - `--coadd-cubes` option added to co-add cubes on the final imaging round together to form a single field spectral cube +- Cleaning up the `flint_masking` CLI: + - added more options + - removed references to butterworth filter + - marked `minimum_boxcar)artefact_mask` as deprecated and to be removed +- Added initial `beam_shape_erode` to masking operations # 0.2.6 From 0b04ec87fae3a5439e3b432033c9dc64a50f4965 Mon Sep 17 00:00:00 2001 From: tgalvin Date: Tue, 8 Oct 2024 21:23:32 +0800 Subject: [PATCH 07/11] added a couple of todos, rearranging --- flint/masking.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/flint/masking.py b/flint/masking.py index 84d116d1..a323569a 100644 --- a/flint/masking.py +++ b/flint/masking.py @@ -670,11 +670,14 @@ def create_snr_mask_from_fits( fits_image=fits_image_path, include_signal_path=create_signal_fits ) + # TODOL Make the bkg and rms images optional. Don't need to load if mbc is usede with fits.open(fits_image_path) as fits_image: fits_header = fits_image[0].header # type: ignore - with fits.open(fits_bkg_path) as fits_bkg: - logger.info("Subtracting background") - signal_data = fits_image[0].data - fits_bkg[0].data # type: ignore + signal_data = fits_image[0].data # type: ignore + + with fits.open(fits_bkg_path) as fits_bkg: + logger.info("Subtracting background") + signal_data -= fits_bkg[0].data # type: ignore with fits.open(fits_rms_path) as fits_rms: logger.info("Dividing by RMS") @@ -698,7 +701,8 @@ def create_snr_mask_from_fits( # case of a fits file, the file may either contain a single frequency or it may # contain a cube of images. if masking_options.flood_fill: - # TODO: This function should really just accept a MaskingOptions directly + # TODO: The image and signal masks both don't need to be inputs. Image is only used + # if mbc = True mask_data = reverse_negative_flood_fill( signal=np.squeeze(signal_data), image=np.squeeze(fits.getdata(fits_image_path)), From b5a49e8fc86d5d7900668a525ecd2badafc6298d Mon Sep 17 00:00:00 2001 From: tgalvin Date: Wed, 9 Oct 2024 15:54:07 +0800 Subject: [PATCH 08/11] initial tweaks to better clean API --- flint/masking.py | 129 ++++++++++++++++++++------------ flint/prefect/common/imaging.py | 2 - tests/test_masking.py | 46 +++++++++++- 3 files changed, 125 insertions(+), 52 deletions(-) diff --git a/flint/masking.py b/flint/masking.py index a323569a..35a06e25 100644 --- a/flint/masking.py +++ b/flint/masking.py @@ -517,11 +517,8 @@ def _verify_set_positive_seed_clip( def reverse_negative_flood_fill( + base_image: np.ndarray, masking_options: MaskingOptions, - image: Optional[np.ndarray] = None, - rms: Optional[np.ndarray] = None, - background: Optional[np.ndarray] = None, - signal: Optional[np.ndarray] = None, pixels_per_beam: Optional[float] = None, ) -> np.ndarray: """Attempt to: @@ -549,15 +546,12 @@ def reverse_negative_flood_fill( * if there are bright negative artefacts there are likely bright positive artefacts nearby * such negative pixels are ~10% level artefacts from a presumed bright sources - Optionally, the `growlow_snr_mask` may also be considered via the `grow_low_snr` and `grow_low_island_size` + Optionally, the `grow_low_snr_mask` may also be considered via the `grow_low_snr` and `grow_low_island_size` parameters. Args: + base_image (np.ndarray): The base image or signal map that is used throughout the fill procedure. masking_options (MaskingOptions): Options to carry out masking. - image (Optional[np.ndarray], optional): The total intensity pixels to have the mask for. Defaults to None. - rms (Optional[np.ndarray], optional): The noise across the image. Defaults to None. - background (Optional[np.ndarray], optional): The background acros the image. If None, zeros are assumed. Defaults to None. - signal(Optional[np.ndarray], optional): A signal map. Defaults to None. pixels_per_beam (Optional[float], optional): The number of pixels that cover a beam. Defaults to None. Returns: @@ -567,18 +561,14 @@ def reverse_negative_flood_fill( logger.info("Will be reversing flood filling") logger.info(f"{masking_options=} ") - signal = _get_signal_image( - image=image, rms=rms, background=background, signal=signal - ) - - if masking_options.flood_fill_use_mbc and image is not None: + if masking_options.flood_fill_use_mbc: positive_mask = minimum_absolute_clip( - image=image, + image=base_image, increase_factor=masking_options.flood_fill_positive_seed_clip, box_size=masking_options.flood_fill_use_mbc_box_size, ) flood_floor_mask = minimum_absolute_clip( - image=image, + image=base_image, increase_factor=masking_options.flood_fill_positive_flood_clip, box_size=masking_options.flood_fill_use_mbc_box_size, ) @@ -586,13 +576,13 @@ def reverse_negative_flood_fill( # Sanity check the upper clip level, you rotten seadog positive_seed_clip = _verify_set_positive_seed_clip( positive_seed_clip=masking_options.flood_fill_positive_seed_clip, - signal=signal, + signal=base_image, ) # Here we create the mask image that will start the binary dilation # process, and we will ensure only pixels above the `positive_flood_clip` # are allowed to be dilated. In other words we are growing the mask - positive_mask = signal >= positive_seed_clip - flood_floor_mask = signal > masking_options.flood_fill_positive_flood_clip + positive_mask = base_image >= positive_seed_clip + flood_floor_mask = base_image > masking_options.flood_fill_positive_flood_clip positive_dilated_mask = scipy_binary_dilation( input=positive_mask, @@ -603,7 +593,7 @@ def reverse_negative_flood_fill( if masking_options.minimum_boxcar: positive_dilated_mask = minimum_boxcar_artefact_mask( - signal=signal, + signal=base_image, island_mask=positive_dilated_mask, boxcar_size=masking_options.minimum_boxcar_size, increase_factor=masking_options.minimum_boxcar_increase_factor, @@ -612,7 +602,7 @@ def reverse_negative_flood_fill( negative_dilated_mask = None if masking_options.suppress_artefacts: negative_dilated_mask = suppress_artefact_mask( - signal=signal, + signal=base_image, negative_seed_clip=masking_options.suppress_artefacts_negative_seed_clip, guard_negative_dilation=masking_options.suppress_artefacts_guard_negative_dilation, pixels_per_beam=pixels_per_beam, @@ -624,7 +614,7 @@ def reverse_negative_flood_fill( if masking_options.grow_low_snr_island: low_snr_mask = grow_low_snr_mask( - signal=signal, + signal=base_image, grow_low_snr=masking_options.grow_low_snr_island_clip, grow_low_island_size=masking_options.grow_low_snr_island_size, region_mask=negative_dilated_mask, @@ -634,17 +624,57 @@ def reverse_negative_flood_fill( return positive_dilated_mask.astype(np.int32) +def _create_signal_from_rmsbkg( + image: Union[Path, np.ndarray], + rms: Union[Path, np.ndarray], + bkg: Union[Path, np.ndarray], +) -> np.ndarray: + logger.info("Creating signal image") + + if isinstance(image, Path): + with fits.open(image) as in_fits: + logger.info(f"Loading {image}") + image = in_fits[0].data # type: ignore + + assert isinstance( + image, np.ndarray + ), f"Expected the image to be a numpy array by now, instead have {type(image)}" + + if isinstance(bkg, Path): + with fits.open(bkg) as in_fits: + logger.info(f"Loading {bkg=}") + bkg = in_fits[0].data # type: ignore + + logger.info("Subtracting background") + image -= bkg + + if isinstance(rms, Path): + with fits.open(rms) as in_fits: + logger.info(f"Loading {rms=}") + rms = in_fits[0].data # type: ignore + + logger.info("Dividing by rms") + image /= rms + + return np.array(image) + + +def _need_to_make_signal(masking_options: MaskingOptions) -> bool: + """Isolated functions to consider whether a signal image is needed""" + return not masking_options.flood_fill_use_mbc + + def create_snr_mask_from_fits( fits_image_path: Path, - fits_rms_path: Path, - fits_bkg_path: Path, masking_options: MaskingOptions, + fits_rms_path: Optional[Path], + fits_bkg_path: Optional[Path], create_signal_fits: bool = False, overwrite: bool = True, ) -> FITSMaskNames: """Create a mask for an input FITS image based on a signal to noise given a corresponding pair of RMS and background FITS images. - Internally the signal image is computed as something akin to: + Internally should a signal image be needed it is computed as something akin to: > signal = (image - background) / rms This is done in a staged manner to minimise the number of (potentially large) images @@ -653,13 +683,15 @@ def create_snr_mask_from_fits( Each of the input images needs to share the same shape. This means that compression features offered by some tooling (e.g. BANE --compress) can not be used. + Depending on the `MaksingOptions` used the signal image may not be needed. + Once the signal map as been computed, all pixels below ``min_snr`` are flagged. Args: fits_image_path (Path): Path to the FITS file containing an image - fits_rms_path (Path): Path to the FITS file with an RMS image corresponding to ``fits_image_path`` - fits_bkg_path (Path): Path to the FITS file with an baclground image corresponding to ``fits_image_path`` masking_options (MaskingOptions): Configurables on the masking operation procedure. + fits_rms_path (Optional[Path], optional): Path to the FITS file with an RMS image corresponding to ``fits_image_path``. Defaults to None. + fits_bkg_path (Optional[Path], optional): Path to the FITS file with an background image corresponding to ``fits_image_path``. Defaults to None. create_signal_fits (bool, optional): Create an output signal map. Defaults to False. overwrite (bool): Passed to `fits.writeto`, and will overwrite files should they exist. Defaults to True. @@ -675,23 +707,23 @@ def create_snr_mask_from_fits( fits_header = fits_image[0].header # type: ignore signal_data = fits_image[0].data # type: ignore - with fits.open(fits_bkg_path) as fits_bkg: - logger.info("Subtracting background") - signal_data -= fits_bkg[0].data # type: ignore - - with fits.open(fits_rms_path) as fits_rms: - logger.info("Dividing by RMS") - signal_data /= fits_rms[0].data # type: ignore - - if create_signal_fits: - logger.info(f"Writing {mask_names.signal_fits}") - fits.writeto( - filename=mask_names.signal_fits, - data=signal_data, - header=fits_header, - overwrite=overwrite, + if _need_to_make_signal(masking_options=masking_options): + assert isinstance(fits_rms_path, Path) and isinstance( + fits_bkg_path, Path + ), "Expected paths for input RMS and bkg FITS files" + signal_data = _create_signal_from_rmsbkg( + image=signal_data, rms=fits_rms_path, bkg=fits_bkg_path ) + if create_signal_fits: + logger.info(f"Writing {mask_names.signal_fits}") + fits.writeto( + filename=mask_names.signal_fits, + data=signal_data, + header=fits_header, + overwrite=overwrite, + ) + pixels_per_beam = get_pixels_per_beam(fits_path=fits_image_path) # Following the help in wsclean: @@ -704,8 +736,7 @@ def create_snr_mask_from_fits( # TODO: The image and signal masks both don't need to be inputs. Image is only used # if mbc = True mask_data = reverse_negative_flood_fill( - signal=np.squeeze(signal_data), - image=np.squeeze(fits.getdata(fits_image_path)), + base_image=np.squeeze(signal_data), masking_options=masking_options, pixels_per_beam=pixels_per_beam, ) @@ -743,14 +774,14 @@ def get_parser() -> ArgumentParser: fits_parser = subparser.add_parser( "mask", - help="Create a mask for an image, using its RMS and BKG images (e.g. outputs from BANE). Output FITS image will default to the image with a mask suffix.", + help="Create a mask for an image, potentially using its RMS and BKG images (e.g. outputs from BANE). Output FITS image will default to the image with a mask suffix.", ) fits_parser.add_argument("image", type=Path, help="Path to the input image. ") fits_parser.add_argument( - "rms", type=Path, help="Path to the RMS of the input image. " + "--rms-fits", type=Path, help="Path to the RMS of the input image. " ) fits_parser.add_argument( - "bkg", type=Path, help="Path to the BKG of the input image. " + "--bkg-fits", type=Path, help="Path to the BKG of the input image. " ) fits_parser.add_argument( @@ -852,8 +883,8 @@ def cli(): masking_options = _args_to_mask_options(args=args) create_snr_mask_from_fits( fits_image_path=args.image, - fits_rms_path=args.rms, - fits_bkg_path=args.bkg, + fits_rms_path=args.rms_fits, + fits_bkg_path=args.bkg_fits, create_signal_fits=args.save_signal, masking_options=masking_options, ) diff --git a/flint/prefect/common/imaging.py b/flint/prefect/common/imaging.py index 60359dd8..5f265b17 100644 --- a/flint/prefect/common/imaging.py +++ b/flint/prefect/common/imaging.py @@ -816,7 +816,6 @@ def _create_convolve_linmos_cubes( def task_create_image_mask_model( image: Union[LinmosCommand, ImageSet, WSCleanCommand], image_products: AegeanOutputs, - min_snr: Optional[float] = 3.5, update_masking_options: Optional[Dict[str, Any]] = None, ) -> FITSMaskNames: """Create a mask from a linmos image, with the intention of providing it as a clean mask @@ -825,7 +824,6 @@ def task_create_image_mask_model( Args: linmos_parset (LinmosCommand): Linmos command and associated meta-data image_products (AegeanOutputs): Images of the RMS and BKG - min_snr (float, optional): The minimum S/N a pixel should be for it to be included in the clean mask. update_masking_options (Optional[Dict[str,Any]], optional): Updated options supplied to the default MaskingOptions. Defaults to None. diff --git a/tests/test_masking.py b/tests/test_masking.py index 968ccdec..6e46ed96 100644 --- a/tests/test_masking.py +++ b/tests/test_masking.py @@ -6,7 +6,9 @@ from flint.masking import ( MaskingOptions, + _create_signal_from_rmsbkg, _args_to_mask_options, + _need_to_make_signal, _verify_set_positive_seed_clip, beam_shape_erode, consider_beam_mask_round, @@ -20,12 +22,54 @@ SHAPE = (100, 100) +def test_create_signal_from_rmsbkg(): + """Make sure the operations around the creation of a signal image from + rms / bkg images follows and handles paths vs numpy arrays""" + shape = (100, 100) + image = np.zeros(shape) + 10 + bkg = np.ones(shape) + rms = np.zeros(shape) + 3 + + signal = _create_signal_from_rmsbkg(image=image, rms=rms, bkg=bkg) + assert np.all(signal == 3.0) + + +def test_create_signal_from_rmsbkg_with_fits(tmpdir): + """Make sure the operations around the creation of a signal image from + rms / bkg images follows and handles paths vs numpy arrays. Unlike the + above this uses fits for some inputs""" + out_path = Path(tmpdir) / "fits" + out_path.mkdir(parents=True, exist_ok=True) + + shape = (100, 100) + image = np.zeros(shape) + 10 + image_fits = Path(out_path) / "image.fits" + fits.writeto(image_fits, data=image) + + bkg = np.ones(shape) + rms = np.zeros(shape) + 3 + rms_fits = Path(out_path) / "rms.fits" + fits.writeto(rms_fits, data=rms) + + signal = _create_signal_from_rmsbkg(image=image_fits, rms=rms_fits, bkg=bkg) + assert np.all(signal == 3.0) + + +def test_need_to_make_signal(): + """Ensure the conditions around creating a signal image make sense""" + masking_options = MaskingOptions() + assert _need_to_make_signal(masking_options=masking_options) + + masking_options = MaskingOptions(flood_fill_use_mbc=True) + assert not _need_to_make_signal(masking_options=masking_options) + + def test_arg_parser_cli_and_masking_options(): """See if the CLI parser is constructed with correct set of options which are properly converted to a MaskingOptions object""" parser = get_parser() args = parser.parse_args( - args="mask img rms bkg --flood-fill --flood-fill-positive-seed-clip 10 --flood-fill-positive-flood-clip 1. --flood-fill-use-mbc --flood-fill-use-mbc-box-size 100".split() + args="mask img --rms-fits rms --bkg-fits bkg --flood-fill --flood-fill-positive-seed-clip 10 --flood-fill-positive-flood-clip 1. --flood-fill-use-mbc --flood-fill-use-mbc-box-size 100".split() ) masking_options = _args_to_mask_options(args=args) assert isinstance(masking_options, MaskingOptions) From 8a8118d04be2618eb217c56843a92552bb88738b Mon Sep 17 00:00:00 2001 From: tgalvin Date: Wed, 9 Oct 2024 16:34:07 +0800 Subject: [PATCH 09/11] matching dimensionality in beam shape erode --- flint/masking.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/flint/masking.py b/flint/masking.py index 35a06e25..c59ce2e7 100644 --- a/flint/masking.py +++ b/flint/masking.py @@ -185,11 +185,16 @@ def beam_shape_erode( fits_header=fits_header, minimum_response=minimum_response ) + # This handles any unsqueezed dimensions + beam_mask_kernel = beam_mask_kernel.reshape( + mask.shape[:-2] + beam_mask_kernel.shape + ) + erode_mask = scipy_binary_erosion( input=mask, iterations=1, structure=beam_mask_kernel ) - return erode_mask + return erode_mask.astype(mask.dtype) def extract_beam_mask_from_mosaic( @@ -879,7 +884,7 @@ def cli(): args = parser.parse_args() - if args.mode == "snrmask": + if args.mode == "mask": masking_options = _args_to_mask_options(args=args) create_snr_mask_from_fits( fits_image_path=args.image, From 469e7b6a0605abcea98e8d993fb07a4f3b00e331 Mon Sep 17 00:00:00 2001 From: tgalvin Date: Wed, 9 Oct 2024 16:56:10 +0800 Subject: [PATCH 10/11] added doc --- flint/masking.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/flint/masking.py b/flint/masking.py index c59ce2e7..ffc726a7 100644 --- a/flint/masking.py +++ b/flint/masking.py @@ -134,7 +134,9 @@ def create_beam_mask_kernel( Returns: np.ndarray: Boolean mask of the kernel shape """ - assert minimum_response > 0, f"{minimum_response=}, should be positive" + assert ( + 0.0 < minimum_response < 1.0 + ), f"{minimum_response=}, should be between 0 to 1 (exclusive)" POSITION_KEYS = ("CDELT1", "CDELT2") if not all([key in fits_header for key in POSITION_KEYS]): @@ -180,7 +182,7 @@ def beam_shape_erode( ) return mask - logger.info("Eroding the mask using the beam shape") + logger.info(f"Eroding the mask using the beam shape with {minimum_response=}") beam_mask_kernel = create_beam_mask_kernel( fits_header=fits_header, minimum_response=minimum_response ) @@ -841,7 +843,7 @@ def get_parser() -> ArgumentParser: "--beam-shape-erode-minimum-response", type=float, default=0.6, - help="The minimum response of the beam that is used to form t he erode structure shape", + help="The minimum response of the beam that is used to form t he erode structure shape. Smaller numbers correspond to a larger shape which means islands are more aggressively removed", ) extract_parser = subparser.add_parser( From a5b8161990b9bcf6eeb39613ea8486728769acf4 Mon Sep 17 00:00:00 2001 From: tgalvin Date: Wed, 9 Oct 2024 22:34:32 +0800 Subject: [PATCH 11/11] removed old min_snr argument --- flint/prefect/flows/continuum_pipeline.py | 1 - 1 file changed, 1 deletion(-) diff --git a/flint/prefect/flows/continuum_pipeline.py b/flint/prefect/flows/continuum_pipeline.py index 2b31a4fe..b196d6e2 100644 --- a/flint/prefect/flows/continuum_pipeline.py +++ b/flint/prefect/flows/continuum_pipeline.py @@ -400,7 +400,6 @@ def process_science_fields( fits_beam_masks = task_create_image_mask_model.map( image=wsclean_cmds, image_products=beam_aegean_outputs, - min_snr=3.5, update_masking_options=unmapped(masking_options), )