diff --git a/flint/configuration.py b/flint/configuration.py index 69a45fb8..a30c9537 100644 --- a/flint/configuration.py +++ b/flint/configuration.py @@ -36,7 +36,6 @@ def get_selfcal_options_from_yaml(input_yaml: Optional[Path] = None) -> Dict: 2: {"solint": "30s", "calmode": "p", "uvrange": ">235m", "nspw": 4}, 3: {"solint": "60s", "calmode": "ap", "uvrange": ">235m", "nspw": 4}, 4: {"solint": "30s", "calmode": "ap", "uvrange": ">235m", "nspw": 4}, - 5: {"solint": "30s", "calmode": "ap", "uvrange": ">235m", "nspw": 1}, } diff --git a/flint/masking.py b/flint/masking.py index 132aabfc..caea8cf0 100644 --- a/flint/masking.py +++ b/flint/masking.py @@ -2,9 +2,11 @@ thought being towards FITS images. """ +from __future__ import annotations + from argparse import ArgumentParser from pathlib import Path -from typing import Optional, Tuple +from typing import NamedTuple, Optional import numpy as np from astropy.io import fits @@ -17,13 +19,46 @@ binary_erosion as scipy_binary_erosion, # Rename to distinguish from skimage ) from scipy.ndimage import label -from skimage.filters import butterworth -from skimage.morphology import binary_erosion from flint.logging import logger from flint.naming import FITSMaskNames, create_fits_mask_names +class MaskingOptions(NamedTuple): + """Contains options for the creation of clean masks from some subject + image. Clipping levels specified are in units of RMS (or sigma). They + are NOT in absolute units. + """ + + 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. """ + flood_fill_positive_seed_clip: float = 4.5 + """The clipping level to seed islands that will be grown to lower SNR""" + flood_fill_positive_flood_clip: float = 1.5 + """Clipping level used to grow seeded islands down to""" + suppress_artefacts: bool = True + """Whether to attempt artefacts based on the presence of sigificant negatives""" + suppress_artefacts_negative_seed_clip: Optional[float] = 5 + """The significance level of a negative island for the sidelobe suppresion to be activated. This should be a positive number (the signal map is internally inverted)""" + suppress_artefacts_guard_negative_dilation: float = 40 + """The minimum positive signifance pixels should have to be guarded when attempting to suppress artefacts around bright sources""" + grow_low_snr_islands: bool = True + """Whether to attempt to grow a mask to capture islands of low SNR (e.g. diffuse emission)""" + grow_low_snr_clip: float = 1.75 + """The minimum signifance 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""" + + def with_options(self, **kwargs) -> MaskingOptions: + """Return a new instance of the MaskingOptions""" + _dict = self._asdict() + _dict.update(**kwargs) + + return MaskingOptions(**_dict) + + def extract_beam_mask_from_mosaic( fits_beam_image_path: Path, fits_mosaic_mask_names: FITSMaskNames ) -> FITSMaskNames: @@ -165,7 +200,8 @@ def reverse_negative_flood_fill( signal: Optional[np.ndarray] = None, positive_seed_clip: float = 4, positive_flood_clip: float = 2, - negative_seed_clip: Optional[float] = 5, + suppress_artefacts: bool = False, + negative_seed_clip: float = 5, guard_negative_dilation: float = 50, grow_low_snr: Optional[float] = 2, grow_low_island_size: int = 512, @@ -205,6 +241,7 @@ def reverse_negative_flood_fill( signal(Optional[np.ndarray], optional): A signal map. Defaults to None. positive_seed_clip (float, optional): Initial clip of the mask before islands are grown. Defaults to 4. positive_flood_clip (float, optional): Pixels above `positive_seed_clip` are dilated to this threshold. Defaults to 2. + suppress_artefacts (boo, optional): Attempt to suppress regions around presumed artefacts. Defaults to False. negative_seed_clip (Optional[float], optional): Initial clip of negative pixels. This operation is on the inverted signal mask (so this value should be a positive number). If None this second operation is not performed. Defaults to 5. guard_negative_dilation (float, optional): Positive pixels from the computed signal mask will be above this threshold to be protect from the negative island mask dilation. Defaults to 50. grow_low__snr (Optional[float], optional): Attempt to grow islands of contigous pixels above thius low SNR ration. If None this is not performed. Defaults to 2. @@ -239,14 +276,16 @@ def reverse_negative_flood_fill( structure=np.ones((3, 3)), ) + # TODO: This function should be divided up + negative_dilated_mask = None + # Now do the same but on negative islands. The assumption here is that: # - no genuine source of negative sky emission # - negative islands are around bright sources with deconvolution/calibration errors # - if there are brightish negative islands there is also positive brightish arteefact islands nearby # For this reason the guard mask should be sufficently high to protect the # main source but nuke the fask positive islands - negative_dilated_mask = None - if negative_seed_clip: + if suppress_artefacts: negative_mask = negative_signal > negative_seed_clip negative_dilated_mask = scipy_binary_dilation( input=negative_mask, @@ -270,110 +309,12 @@ def reverse_negative_flood_fill( return positive_dilated_mask.astype(np.int32) -def create_snr_mask_wbutter_from_fits( - fits_image_path: Path, - fits_rms_path: Path, - fits_bkg_path: Path, - create_signal_fits: bool = False, - min_snr: float = 5, - connectivity_shape: Tuple[int, int] = (4, 4), - 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: - > signal = (image - background) / rms - - Before deriving a signal map the image is first smoothed using a butterworth filter, and a - crude rescaling factor is applied based on the ratio of the maximum pixel values before and - after the smoothing is applied. - - This is done in a staged manner to minimise the number of (potentially large) images - held in memory. - - 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. - - Once the signal map as been computed, all pixels below ``min_snr`` are flagged. The resulting - islands then have a binary erosion applied to contract the resultingn islands. - - 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`` - create_signal_fits (bool, optional): Create an output signal map. Defaults to False. - min_snr (float, optional): Minimum signal-to-noise ratio for the masking to include a pixel. Defaults to 3.5. - connectivity_shape (Tuple[int, int], optional): The connectivity matrix used in the scikit-image binary erosion applied to the mask. Defaults to (4, 4). - overwrite (bool): Passed to `fits.writeto`, and will overwrite files should they exist. Defaults to True. - - Returns: - FITSMaskNames: Container describing the signal and mask FITS image paths. If ``create_signal_path`` is None, then the ``signal_fits`` attribute will be None. - """ - logger.info(f"Creating a mask image with SNR>{min_snr:.2f}") - mask_names = create_fits_mask_names( - fits_image=fits_image_path, include_signal_path=create_signal_fits - ) - - with fits.open(fits_image_path) as fits_image: - fits_header = fits_image[0].header - - image_max = np.nanmax(fits_image[0].data) - image_butter = butterworth( - np.nan_to_num(np.squeeze(fits_image[0].data)), 0.045, high_pass=False - ) - - butter_max = np.nanmax(image_butter) - scale_ratio = image_max / butter_max - logger.info(f"Scaling smoothed image by {scale_ratio:.4f}") - image_butter *= image_max / butter_max - - with fits.open(fits_bkg_path) as fits_bkg: - logger.info("Subtracting background") - signal_data = image_butter - np.squeeze(fits_bkg[0].data) - - with fits.open(fits_rms_path) as fits_rms: - logger.info("Dividing by RMS") - signal_data /= np.squeeze(fits_rms[0].data) - - 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, - ) - - # Following the help in wsclean: - # WSClean accepts masks in CASA format and in fits file format. A mask is a - # normal, single polarization image file, where all zero values are interpreted - # as being not masked, and all non-zero values are interpreted as masked. In the - # case of a fits file, the file may either contain a single frequency or it may - # contain a cube of images. - logger.info(f"Clipping using a {min_snr=}") - mask_data = (signal_data > min_snr).astype(np.int32) - - logger.info(f"Applying binary erosion with {connectivity_shape=}") - mask_data = binary_erosion(mask_data, np.ones(connectivity_shape)) - - logger.info(f"Writing {mask_names.mask_fits}") - fits.writeto( - filename=mask_names.mask_fits, - data=mask_data.astype(np.int32), - header=fits_header, - overwrite=overwrite, - ) - - return mask_names - - def create_snr_mask_from_fits( fits_image_path: Path, fits_rms_path: Path, fits_bkg_path: Path, + masking_options: MaskingOptions, create_signal_fits: bool = False, - min_snr: float = 3.5, - attempt_reverse_nergative_flood_fill: bool = True, 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. @@ -393,15 +334,13 @@ def create_snr_mask_from_fits( 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. create_signal_fits (bool, optional): Create an output signal map. Defaults to False. - min_snr (float, optional): Minimum signal-to-noise ratio for the masking to include a pixel. Defaults to 3.5. - attempt_negative_flood_fill (bool): Attempt to filter out negative sidelobes from the bask. See `reverse_negative_flood_fill`. Defaults to True. overwrite (bool): Passed to `fits.writeto`, and will overwrite files should they exist. Defaults to True. Returns: FITSMaskNames: Container describing the signal and mask FITS image paths. If ``create_signal_path`` is None, then the ``signal_fits`` attribute will be None. """ - logger.info(f"Creating a mask image with SNR>{min_snr:.2f}") mask_names = create_fits_mask_names( fits_image=fits_image_path, include_signal_path=create_signal_fits ) @@ -431,20 +370,20 @@ def create_snr_mask_from_fits( # as being not masked, and all non-zero values are interpreted as masked. In the # case of a fits file, the file may either contain a single frequency or it may # contain a cube of images. - if attempt_reverse_nergative_flood_fill: + if masking_options.flood_fill: mask_data = reverse_negative_flood_fill( signal=np.squeeze(signal_data), - positive_seed_clip=4, - positive_flood_clip=1.5, - negative_seed_clip=4, - guard_negative_dilation=30, - grow_low_snr=1.75, - grow_low_island_size=768, + positive_seed_clip=masking_options.flood_fill_positive_seed_clip, + positive_flood_clip=masking_options.flood_fill_positive_flood_clip, + negative_seed_clip=masking_options.suppress_artefacts_negative_seed_clip, + guard_negative_dilation=masking_options.suppress_artefacts_guard_negative_dilation, + grow_low_snr=masking_options.grow_low_snr_clip, + grow_low_island_size=masking_options.grow_low_snr_island_size, ) mask_data = mask_data.reshape(signal_data.shape) else: - logger.info(f"Clipping using a {min_snr=}") - mask_data = (signal_data > min_snr).astype(np.int32) + logger.info(f"Clipping using a {masking_options.base_snr_clip=}") + mask_data = (signal_data > masking_options.base_snr_clip).astype(np.int32) logger.info(f"Writing {mask_names.mask_fits}") fits.writeto( @@ -529,23 +468,14 @@ def cli(): args = parser.parse_args() if args.mode == "snrmask": - if args.user_butterworth: - create_snr_mask_wbutter_from_fits( - fits_image_path=args.image, - fits_rms_path=args.rms, - fits_bkg_path=args.bkg, - create_signal_fits=args.save_signal, - min_snr=args.min_snr, - connectivity_shape=tuple(args.connectivity_shape), - ) - else: - create_snr_mask_from_fits( - fits_image_path=args.image, - fits_rms_path=args.rms, - fits_bkg_path=args.bkg, - create_signal_fits=args.save_signal, - min_snr=args.min_snr, - ) + masking_options = MaskingOptions(base_snr_clip=args.min_snr) + create_snr_mask_from_fits( + fits_image_path=args.image, + fits_rms_path=args.rms, + fits_bkg_path=args.bkg, + create_signal_fits=args.save_signal, + masking_options=masking_options, + ) elif args.mode == "extractmask": extract_beam_mask_from_mosaic( diff --git a/flint/options.py b/flint/options.py index 39db2cf6..b6c0d6a8 100644 --- a/flint/options.py +++ b/flint/options.py @@ -90,11 +90,9 @@ class FieldOptions(NamedTuple): """Primary beam attentuation cutoff to use during linmos""" use_preflagger: bool = True """Whether to apply (or search for solutions with) bandpass solutions that have gone through the preflagging operations""" - use_smoothed: bool = True + use_smoothed: bool = False """Whether to apply (or search for solutions with) a bandpass smoothing operation applied""" use_beam_masks: bool = True """Construct beam masks from MFS images to use for the next round of imaging. """ use_beam_masks_from: int = 2 """If `use_beam_masks` is True, start using them from this round of self-calibration""" - use_beam_mask_wbutterworth: bool = False - """If `use_beam_masks` is True, this will specify whether a Butterworth filter is used to smooth the image before the S/N clip is applied""" diff --git a/flint/prefect/common/imaging.py b/flint/prefect/common/imaging.py index a5be6b88..9a4eedca 100644 --- a/flint/prefect/common/imaging.py +++ b/flint/prefect/common/imaging.py @@ -23,8 +23,8 @@ from flint.imager.wsclean import ImageSet, WSCleanCommand, wsclean_imager from flint.logging import logger from flint.masking import ( + MaskingOptions, create_snr_mask_from_fits, - create_snr_mask_wbutter_from_fits, extract_beam_mask_from_mosaic, ) from flint.ms import MS, preprocess_askap_ms, rename_column_in_ms, split_by_field @@ -462,7 +462,7 @@ def task_create_image_mask_model( image: Union[LinmosCommand, ImageSet, WSCleanCommand], image_products: AegeanOutputs, min_snr: Optional[float] = 3.5, - with_butterworth: bool = False, + 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 to an appropriate imager. This is derived using a simple signal to noise cut. @@ -471,7 +471,8 @@ def task_create_image_mask_model( 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. - with_butterworth (bool, optional): whether to taper the input image with a Butterworth filter before masking. + update_masking_options (Optional[Dict[str,Any]], optional): Updated options supplied to the default MaskingOptions. Defaults to None. + Raises: ValueError: Raised when ``image_products`` are not known @@ -502,23 +503,17 @@ def task_create_image_mask_model( logger.info(f"Using {source_rms=}") logger.info(f"Using {source_bkg=}") - if with_butterworth: - mask_names = create_snr_mask_wbutter_from_fits( - fits_image_path=source_image, - fits_bkg_path=source_bkg, - fits_rms_path=source_rms, - create_signal_fits=True, - min_snr=min_snr, - connectivity_shape=(3, 3), - ) - else: - mask_names = create_snr_mask_from_fits( - fits_image_path=source_image, - fits_bkg_path=source_bkg, - fits_rms_path=source_rms, - create_signal_fits=True, - min_snr=min_snr, - ) + masking_options = MaskingOptions() + if update_masking_options: + masking_options = masking_options.with_options(**update_masking_options) + + mask_names = create_snr_mask_from_fits( + fits_image_path=source_image, + fits_bkg_path=source_bkg, + fits_rms_path=source_rms, + masking_options=masking_options, + create_signal_fits=True, + ) logger.info(f"Created {mask_names.mask_fits}") diff --git a/flint/prefect/flows/continuum_pipeline.py b/flint/prefect/flows/continuum_pipeline.py index d9505deb..a7d38107 100644 --- a/flint/prefect/flows/continuum_pipeline.py +++ b/flint/prefect/flows/continuum_pipeline.py @@ -263,7 +263,6 @@ def process_science_fields( image=wsclean_cmds, image_products=beam_aegean_outputs, min_snr=3.5, - with_butterworth=field_options.use_beam_mask_wbutterworth, ) wsclean_options["auto_mask"] = 1.25 wsclean_options["auto_threshold"] = 1.0 @@ -482,12 +481,6 @@ def get_parser() -> ArgumentParser: default=False, help="Whether to use (or search for solutions with) the preflagger operations applied to the bandpass gain solutions", ) - parser.add_argument( - "--use-smoothed", - action="store_true", - default=False, - help="Whether to use (or search for solutions with) the smoothing operations applied to the bandpass gain solutions", - ) parser.add_argument( "--use-beam-masks", default=False, @@ -500,12 +493,6 @@ def get_parser() -> ArgumentParser: type=int, help="If --use-beam-masks is provided, this option specifies from which round of self-calibration the masking operation will be used onwards from. ", ) - parser.add_argument( - "--use-beam-masks-wbutterworth", - default=False, - action="store_true", - help="If --use-beam-masks is provided, this will specify whether a Butterworth filter is first used to smooth an image before applying the S/N cut", - ) return parser @@ -537,10 +524,8 @@ def cli() -> None: beam_cutoff=args.beam_cutoff, pb_cutoff=args.pb_cutoff, use_preflagger=args.use_preflagger, - use_smoothed=args.use_smoothed, use_beam_masks=args.use_beam_masks, use_beam_masks_from=args.use_beam_masks_from, - use_beam_mask_wbutterworth=args.use_beam_masks_wbutterworth, ) setup_run_process_science_field( diff --git a/tests/test_masking.py b/tests/test_masking.py index aa3186dc..0582ee7d 100644 --- a/tests/test_masking.py +++ b/tests/test_masking.py @@ -4,7 +4,7 @@ import pytest from astropy.io import fits -from flint.masking import create_snr_mask_from_fits +from flint.masking import MaskingOptions, create_snr_mask_from_fits from flint.naming import FITSMaskNames SHAPE = (100, 100) @@ -28,11 +28,24 @@ def fits_dir(tmpdir): return fits_dir +def test_make_masking_options(): + """Just a dump test to make sure the options structure is ok""" + + masking_options = MaskingOptions() + + assert masking_options.base_snr_clip != -1 + + masking_options = masking_options.with_options(base_snr_clip=-1) + assert masking_options.base_snr_clip == -1 + + def test_fits_masking(fits_dir): + masking_options = MaskingOptions(flood_fill=False) names = create_snr_mask_from_fits( fits_image_path=fits_dir / "image.fits", fits_rms_path=fits_dir / "rms.fits", fits_bkg_path=fits_dir / "bkg.fits", + masking_options=masking_options, ) assert isinstance(names, FITSMaskNames) @@ -45,10 +58,12 @@ def test_fits_masking(fits_dir): def test_fits_masking_with_signal(fits_dir): + masking_options = MaskingOptions(flood_fill=False) names = create_snr_mask_from_fits( fits_image_path=fits_dir / "image.fits", fits_rms_path=fits_dir / "rms.fits", fits_bkg_path=fits_dir / "bkg.fits", + masking_options=masking_options, create_signal_fits=True, ) diff --git a/tests/test_options.py b/tests/test_options.py index 03452dd3..b04009d0 100644 --- a/tests/test_options.py +++ b/tests/test_options.py @@ -2,6 +2,7 @@ somewhat tracked, especially when using an argparse object to create it """ + from pathlib import Path from flint.options import FieldOptions @@ -46,12 +47,10 @@ def test_create_field_options(): beam_cutoff=args.beam_cutoff, pb_cutoff=args.pb_cutoff, use_preflagger=args.use_preflagger, - use_smoothed=args.use_smoothed, ) assert isinstance(field_options, FieldOptions) assert field_options.use_preflagger is False - assert field_options.use_smoothed is False assert field_options.zip_ms is True assert field_options.linmos_residuals is True assert field_options.rounds == 2 @@ -75,7 +74,6 @@ def test_create_field_options2(): --aegean-container '/scratch3/gal16b/containers/aegean.sif' --reference-catalogue-directory '/scratch3/gal16b/reference_catalogues/' --use-preflagger - --use-smoothed """.split() ) @@ -96,12 +94,10 @@ def test_create_field_options2(): beam_cutoff=args.beam_cutoff, pb_cutoff=args.pb_cutoff, use_preflagger=args.use_preflagger, - use_smoothed=args.use_smoothed, ) assert isinstance(field_options, FieldOptions) assert field_options.use_preflagger is True - assert field_options.use_smoothed is True assert field_options.zip_ms is False assert field_options.linmos_residuals is False assert field_options.rounds == 2