diff --git a/CHANGES.rst b/CHANGES.rst index d7beee957..d92674d01 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -81,6 +81,11 @@ tweakreg - Fix a bug due to which source catalog may contain sources outside of the bounding box. [#947] +source_detection +---------------- + +- Support for PSF fitting (optional) for accurate centroids. [#841] + 0.12.0 (2023-08-18) =================== diff --git a/pyproject.toml b/pyproject.toml index 5bce18789..7168ce445 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -21,7 +21,7 @@ dependencies = [ 'numpy >=1.22', 'photutils @ git+https://github.com/astropy/photutils.git', 'pyparsing >=2.4.7', - 'requests >=2.22', + 'requests >=2.26', 'rad >= 0.17.1', # 'rad @ git+https://github.com/spacetelescope/rad.git@main', 'roman_datamodels >= 0.17.1', @@ -34,7 +34,7 @@ dependencies = [ 'spherical-geometry >= 1.2.22', 'stsci.imagestats >= 1.6.3', 'drizzle >= 1.13.7', - 'webbpsf == 1.1.1', + 'webbpsf == 1.2.1', ] dynamic = ['version'] diff --git a/romancal/lib/psf.py b/romancal/lib/psf.py index bd2ca74b3..99e4733a7 100644 --- a/romancal/lib/psf.py +++ b/romancal/lib/psf.py @@ -31,26 +31,26 @@ # Phase C central wavelengths [micron], released by Goddard (Jan 2023): # https://roman.ipac.caltech.edu/sims/Param_db.html#wfi_filters filter_central_wavelengths = { - "WFI_Filter_F062_Center": 0.620, - "WFI_Filter_F087_Center": 0.869, - "WFI_Filter_F106_Center": 1.060, - "WFI_Filter_F129_Center": 1.293, - "WFI_Filter_F146_Center": 1.464, - "WFI_Filter_F158_Center": 1.577, - "WFI_Filter_F184_Center": 1.842, - "WFI_Filter_F213_Center": 2.125, + "F062": 0.620, + "F087": 0.869, + "F106": 1.060, + "F129": 1.293, + "F146": 1.464, + "F158": 1.577, + "F184": 1.842, + "F213": 2.125, } default_finder = DAOStarFinder( # these defaults extracted from the # romancal SourceDetectionStep - fwhm=2.0, - threshold=2.0, + fwhm=1.0, + threshold=0.0, sharplo=0.0, sharphi=1.0, roundlo=-1.0, roundhi=1.0, - peakmax=1000.0, + peakmax=None, ) @@ -58,9 +58,9 @@ def create_gridded_psf_model( path_prefix, filt, detector, - oversample=12, - fov_pixels=12, - sqrt_n_psfs=4, + oversample=11, + fov_pixels=9, + sqrt_n_psfs=2, overwrite=False, buffer_pixels=100, instrument_options=None, @@ -82,7 +82,8 @@ def create_gridded_psf_model( Computed gridded PSF model for this SCA. Examples include: `"SCA01"` or `"SCA18"`. oversample : int, optional - Oversample factor, default is 12. See WebbPSF docs for details [1]_. + Oversample factor, default is 11. See WebbPSF docs for details [1]_. + Choosing an odd number makes the pixel convolution more accurate. fov_pixels : int, optional Field of view width [pixels]. Default is 12. See WebbPSF docs for details [1]_. @@ -150,17 +151,12 @@ def create_gridded_psf_model( if instrument_options is not None: wfi.options.update(instrument_options) - central_wavelength_meters = ( - filter_central_wavelengths[f"WFI_Filter_{filt}_Center"] * 1e-6 * u.m - ) - # Initialize the PSF library inst = gridded_library.CreatePSFLibrary( instrument=wfi, filter_name=filt, detectors=detector.upper(), num_psfs=n_psfs, - monochromatic=central_wavelength_meters, oversample=oversample, fov_pixels=fov_pixels, add_distortion=False, @@ -279,7 +275,7 @@ def fit_psf_to_image_model( """ if grouper is None: # minimum separation before sources are fit simultaneously: - grouper = SourceGrouper(min_separation=20) # [pix] + grouper = SourceGrouper(min_separation=5) # [pix] if fitter is None: fitter = LevMarLSQFitter(calc_uncertainties=True) diff --git a/romancal/lib/tests/test_psf.py b/romancal/lib/tests/test_psf.py index 9fb9108eb..4cc2cb0bf 100644 --- a/romancal/lib/tests/test_psf.py +++ b/romancal/lib/tests/test_psf.py @@ -4,177 +4,175 @@ import os import tempfile +from copy import deepcopy import numpy as np import pytest from astropy import units as u -from astropy.nddata import overlap_slices from photutils.psf import PSFPhotometry from roman_datamodels import maker_utils as testutil from roman_datamodels.datamodels import ImageModel from romancal.lib.psf import create_gridded_psf_model, fit_psf_to_image_model -n_sources = 10 -image_model_shape = (100, 100) +n_trials = 15 +image_model_shape = (50, 50) rng = np.random.default_rng(0) -@pytest.fixture -def setup_inputs(): - def _setup( - nrows=image_model_shape[0], ncols=image_model_shape[1], noise=1.0, seed=None - ): - """ - Return ImageModel of level 2 image. - """ - shape = (nrows, ncols) - wfi_image = testutil.mk_level2_image(shape=shape) +def setup_inputs( + shape=image_model_shape, + noise=1.0, + seed=None, +): + """ + Return ImageModel of level 2 image. + """ + wfi_image = testutil.mk_level2_image(shape=shape) + wfi_image.data = u.Quantity( + np.ones(shape, dtype=np.float32), u.electron / u.s, dtype=np.float32 + ) + wfi_image.meta.filename = "filename" + wfi_image.meta.instrument["optical_element"] = "F087" + + # add noise to data + if noise is not None: + setup_rng = np.random.default_rng(seed or 19) wfi_image.data = u.Quantity( - np.ones(shape, dtype=np.float32), u.electron / u.s, dtype=np.float32 + setup_rng.normal(scale=noise, size=shape), + u.electron / u.s, + dtype=np.float32, ) - wfi_image.meta.filename = "filename" - - # add noise to data - if noise is not None: - setup_rng = np.random.default_rng(seed or 19) - wfi_image.data = u.Quantity( - setup_rng.normal(scale=noise, size=shape), - u.electron / u.s, - dtype=np.float32, - ) - wfi_image.err = noise * np.ones(shape, dtype=np.float32) * u.electron / u.s + wfi_image.err = noise * np.ones(shape, dtype=np.float32) * u.electron / u.s - # add dq array - wfi_image.dq = np.zeros(shape, dtype=np.uint32) + # add dq array + wfi_image.dq = np.zeros(shape, dtype=np.uint32) - # construct ImageModel - mod = ImageModel(wfi_image) + # construct ImageModel + mod = ImageModel(wfi_image) - return mod + filt = mod.meta.instrument["optical_element"] + detector = mod.meta["instrument"]["detector"].replace("WFI", "SCA") - return _setup - - -def add_synthetic_sources( - image_model, - psf_model, - true_x, - true_y, - true_amp, - oversample, - xname="x_0", - yname="y_0", -): - fit_models = [] - - # ensure truths are arrays: - true_x, true_y, true_amp = ( - np.atleast_1d(truth) for truth in [true_x, true_y, true_amp] - ) - - for x, y, amp in zip(true_x, true_y, true_amp): - psf = psf_model.copy() - psf.parameters = [amp, x, y] - fit_models.append(psf) - - synth_image = image_model.data - synth_err = image_model.err - psf_shape = np.array(psf_model.data.shape[1:]) // oversample - - for fit_model in fit_models: - x0 = getattr(fit_model, xname).value - y0 = getattr(fit_model, yname).value - slc_lg, _ = overlap_slices(synth_image.shape, psf_shape, (y0, x0), mode="trim") - yy, xx = np.mgrid[slc_lg] - model_data = fit_model(xx, yy) * image_model.data.unit - model_err = np.sqrt(model_data.value) * model_data.unit - synth_image[slc_lg] += ( - rng.normal( - model_data.to_value(image_model.data.unit), - model_err.to_value(image_model.data.unit), - size=model_data.shape, - ) - * image_model.data.unit - ) - synth_err[slc_lg] = np.sqrt(synth_err[slc_lg] ** 2 + model_err**2) - - -@pytest.mark.webbpsf -@pytest.mark.parametrize( - "dx, dy, true_amp", - zip( - rng.uniform(-1, 1, n_sources), - rng.uniform(-1, 1, n_sources), - np.geomspace(10, 10_000, n_sources), - ), -) -def test_psf_fit(setup_inputs, dx, dy, true_amp, seed=42): # input parameters for PSF model: - filt = "F087" - detector = "SCA01" - oversample = 12 - fov_pixels = 15 + webbpsf_config = dict( + filt=filt, + detector=detector, + oversample=12, + fov_pixels=15, + ) dir_path = tempfile.gettempdir() - filename_prefix = f"psf_model_{filt}" + filename_prefix = f"psf_model_{webbpsf_config['filt']}" file_path = os.path.join(dir_path, filename_prefix) # compute gridded PSF model: psf_model, centroids = create_gridded_psf_model( file_path, - filt, - detector, - oversample=oversample, - fov_pixels=fov_pixels, - overwrite=False, + webbpsf_config["filt"], + webbpsf_config["detector"], + oversample=webbpsf_config["oversample"], + fov_pixels=webbpsf_config["fov_pixels"], + overwrite=True, logging_level="ERROR", ) - # generate an ImageModel - image_model = setup_inputs(seed=seed) - init_data_stddev = np.std(image_model.data.value) + return mod, webbpsf_config, psf_model + - # add synthetic sources to the ImageModel: - true_x = image_model_shape[0] / 2 + dx - true_y = image_model_shape[1] / 2 + dy - add_synthetic_sources( - image_model, psf_model, true_x, true_y, true_amp, oversample=oversample +def add_sources(image_model, psf_model, x_true, y_true, amp_true, background=10): + shape = image_model.data.shape + + x_true, y_true, amp_true = ( + np.atleast_1d(val) for val in [x_true, y_true, amp_true] ) - if fov_pixels % 2 == 0: - fit_shape = (fov_pixels + 1, fov_pixels + 1) - else: - fit_shape = (fov_pixels, fov_pixels) - - # fit the PSF to the ImageModel: - results_table, photometry = fit_psf_to_image_model( - image_model=image_model, - photometry_cls=PSFPhotometry, - psf_model=psf_model, - x_init=true_x, - y_init=true_y, - fit_shape=fit_shape, + fit_models = [] + for amp, x, y in zip(amp_true, x_true, y_true): + model = deepcopy(psf_model) + model.parameters = [amp, x, y] + fit_models.append(model) + + photometry = PSFPhotometry(psf_model, (15, 15)) + photometry.fit_results = dict(local_bkg=np.ones(len(x_true)) * 0) + photometry._fit_models = fit_models + + image = photometry.make_model_image(shape, (19, 19)) + + image += rng.normal(background, 1, size=shape) + image_model.data = image * np.ones_like(image_model.data) + image_model.err = background * np.ones_like(image_model.err) + + +class TestPSFFitting: + def setup_method(self): + self.image_model, self.webbpsf_config, self.psf_model = setup_inputs( + shape=image_model_shape, + ) + + @pytest.mark.webbpsf + @pytest.mark.parametrize( + "dx, dy, true_amp", + zip( + rng.uniform(-1, 1, n_trials), + rng.uniform(-1, 1, n_trials), + np.geomspace(1_000, 100_000, n_trials), + ), ) + def test_psf_fit(self, dx, dy, true_amp): + # generate an ImageModel + image_model = deepcopy(self.image_model) + init_data_stddev = np.std(image_model.data.value) + + # add synthetic sources to the ImageModel: + true_x = image_model_shape[0] / 2 + dx + true_y = image_model_shape[1] / 2 + dy + add_sources(image_model, self.psf_model, true_x, true_y, true_amp) + + if self.webbpsf_config["fov_pixels"] % 2 == 0: + fit_shape = ( + self.webbpsf_config["fov_pixels"] + 1, + self.webbpsf_config["fov_pixels"] + 1, + ) + else: + fit_shape = ( + self.webbpsf_config["fov_pixels"], + self.webbpsf_config["fov_pixels"], + ) - # difference between input and output, normalized by the - # uncertainty. Has units of sigma: - delta_x = np.abs(true_x - results_table["x_fit"]) / results_table["x_err"] - delta_y = np.abs(true_y - results_table["y_fit"]) / results_table["y_err"] + # fit the PSF to the ImageModel: + results_table, photometry = fit_psf_to_image_model( + image_model=image_model, + photometry_cls=PSFPhotometry, + psf_model=self.psf_model, + x_init=true_x, + y_init=true_y, + fit_shape=fit_shape, + ) - sigma_threshold = 3.5 - assert np.all(delta_x < sigma_threshold) - assert np.all(delta_y < sigma_threshold) + # difference between input and output, normalized by the + # uncertainty. Has units of sigma: + delta_x = np.abs(true_x - results_table["x_fit"]) / results_table["x_err"] + delta_y = np.abs(true_y - results_table["y_fit"]) / results_table["y_err"] - # now check that the uncertainties aren't way too large, which could cause - # the above test to pass even when the fits are bad. Use overly-simple approximation - # that astrometric uncertainty be proportional to the PSF's FWHM / SNR: - approx_snr = true_amp / init_data_stddev - approx_fwhm = 1 - approx_centroid_err = approx_fwhm / approx_snr + sigma_threshold = 3.5 + assert np.all(delta_x < sigma_threshold) + assert np.all(delta_y < sigma_threshold) - # centroid err heuristic above is an underestimate, so we scale it up: - scale_factor_approx = 100 + # now check that the uncertainties aren't way too large, which could cause + # the above test to pass even when the fits are bad. Use overly-simple + # approximation that astrometric uncertainty be proportional to the + # PSF's FWHM / SNR: + approx_snr = true_amp / init_data_stddev + approx_fwhm = 1 + approx_centroid_err = approx_fwhm / approx_snr - assert np.all(results_table["x_err"] < scale_factor_approx * approx_centroid_err) - assert np.all(results_table["y_err"] < scale_factor_approx * approx_centroid_err) + # centroid err heuristic above is an underestimate, so we scale it up: + scale_factor_approx = 100 + + assert np.all( + results_table["x_err"] < scale_factor_approx * approx_centroid_err + ) + assert np.all( + results_table["y_err"] < scale_factor_approx * approx_centroid_err + ) diff --git a/romancal/source_detection/source_detection_step.py b/romancal/source_detection/source_detection_step.py index 324b0d822..4d1f40b35 100644 --- a/romancal/source_detection/source_detection_step.py +++ b/romancal/source_detection/source_detection_step.py @@ -5,6 +5,7 @@ import logging +import astropy.units as u import numpy as np from asdf import AsdfFile from astropy.stats import SigmaClip @@ -19,7 +20,7 @@ from roman_datamodels import datamodels as rdm from roman_datamodels import maker_utils -from romancal.lib import dqflags +from romancal.lib import dqflags, psf from romancal.stpipe import RomanStep log = logging.getLogger(__name__) @@ -35,7 +36,7 @@ class SourceDetectionStep(RomanStep): """ spec = """ - kernel_fwhm = float(default=2.) # DAOStarFinder:Size of Gaussian kernel, + kernel_fwhm = float(default=None) # DAOStarFinder:Size of Gaussian kernel, # in pixels. sharplo = float(default=0.) # DAOStarFinder: Lower bound for sharpness. # Typical values of sharpness range from 0 (flat) to 1 (delta function). @@ -49,13 +50,13 @@ class SourceDetectionStep(RomanStep): scalar_threshold = float(default=None) # Detection threshold, to # be used for entire image. Assumed to be in same units as data, and is # an absolute threshold over background. - calc_threshold = boolean(default=True) # Calculate a single absoulte + calc_threshold = boolean(default=True) # Calculate a single absolute # detection threshold from image based on background. snr_threshold = float(default=3.0) # if calc_threshold_img, # the SNR for the threshold image. bkg_estimator = string(default='median') # if calc_threshold_img, # choice of mean, median, or mode. - bkg_boxsize = integer(default=3) # if calc_threshold_img, + bkg_boxsize = integer(default=9) # if calc_threshold_img, # size of box in pixels for 2D background. bkg_sigma = float(default=2.0) # if calc_threshold_img, # n sigma for sigma clipping bkgrnd. @@ -65,6 +66,7 @@ class SourceDetectionStep(RomanStep): # Will overwrite an existing catalog of the same name. output_cat_filetype = option('asdf', 'ecsv', default='asdf') # Used if #save_catalogs=True - file type of output catalog. + fit_psf = boolean(default=False) # fit source PSFs for accurate astrometry """ def process(self, input): @@ -88,6 +90,8 @@ def process(self, input): (dqflags.pixel["DO_NOT_USE"]) & input_model.dq ).astype(bool) + filt = input_model.meta.instrument["optical_element"] + # if a pre-determined threshold value for detection for the whole # image is provided, use this if self.scalar_threshold is not None: @@ -98,14 +102,28 @@ def process(self, input): # image by calculating a 2D background image, using this to create # a 2d threshold image, and using the median of the 2d threshold # image as the scalar detection threshold for the whole image - elif self.calc_threshold is not None: + elif self.calc_threshold: log.info("Determining detection threshold from image.") bkg = self._calc_2D_background() threshold_img = bkg.background + self.snr_threshold * bkg.background_rms threshold = np.median(threshold_img) log.info(f"Calculated a detection threshold of {threshold} from image.") - - log.info("Detecting sources with DAOFind, using entire image array.") + else: + threshold = np.median(self.data) + + if self.kernel_fwhm is None: + # estimate the FWHM of the PSF using rough estimate + # from the diffraction limit: + central_wavelength = psf.filter_central_wavelengths[filt] * u.um + diffraction_limit = float(1.22 * central_wavelength / (2.4 * u.m)) + assume_pixel_scale = 0.11 * u.arcsec / u.pix + expect_psf_pix = np.sqrt( + ((diffraction_limit * u.rad).to(u.arcsec) / assume_pixel_scale) ** 2 + + 1 * u.pix**2 + ).to_value(u.pix) + self.kernel_fwhm = expect_psf_pix + + log.info("Detecting sources with DAOStarFinder, using entire image array.") daofind = DAOStarFinder( fwhm=self.kernel_fwhm, threshold=threshold, @@ -122,8 +140,10 @@ def process(self, input): sources = daofind(self.data, mask=self.coverage_mask) elif self.calc_threshold is not None: - # subtrack background from data if calculating abs. threshold + # subtract background from data if calculating abs. threshold sources = daofind(self.data - bkg.background, mask=self.coverage_mask) + else: + sources = daofind(self.data, mask=self.coverage_mask) # reduce table to minimal number of columns, just source ID, # positions, and fluxes @@ -140,17 +160,54 @@ def process(self, input): dtype=(int, np.float64, np.float64, np.float64), ) - # attach source catalog to output model as array in meta.source_detecion - # the table will be stored as a 1D array with the four columns - # concatenated, in order, with units attached - catalog_as_array = np.array( - [ - catalog["id"].value, - catalog["xcentroid"].value, - catalog["ycentroid"].value, - catalog["flux"].value, - ] - ) + if self.fit_psf and len(sources): + # refine astrometry by fitting PSF models to the detected sources + log.info("Constructing a gridded PSF model.") + detector = input_model.meta.instrument["detector"].replace("WFI", "SCA") + gridded_psf_model, _ = psf.create_gridded_psf_model( + path_prefix="tmp", + filt=filt, + detector=detector, + overwrite=True, + logging_level="ERROR", + ) + + log.info( + "Fitting a PSF model to sources for improved astrometric precision." + ) + psf_photometry_table, photometry = psf.fit_psf_to_image_model( + image_model=input_model, + psf_model=gridded_psf_model, + x_init=catalog["xcentroid"], + y_init=catalog["ycentroid"], + exclude_out_of_bounds=True, + ) + + # attach source catalog to output model as array in + # meta.source_detection the table will be stored as a + # 1D array with the four columns concatenated, in order, + # with units attached + catalog_as_array = np.array( + [ + psf_photometry_table["id"].value, + psf_photometry_table["x_fit"].value, + psf_photometry_table["y_fit"].value, + psf_photometry_table["flux_fit"].value, + ] + ) + else: + # attach source catalog to output model as array in + # meta.source_detection the table will be stored as a + # 1D array with the four columns concatenated, in order, + # with units attached + catalog_as_array = np.array( + [ + catalog["id"].value, + catalog["xcentroid"].value, + catalog["ycentroid"].value, + catalog["flux"].value, + ] + ) # create meta.source detection section in file # if save_catalogs is True, this will be updated with the @@ -185,6 +242,11 @@ def process(self, input): # and save_catalogs is false, since it is not in the schema input_model.meta.source_detection["tweakreg_catalog"] = catalog_as_array + if self.fit_psf: + input_model.meta.source_detection[ + "psf_catalog" + ] = psf_photometry_table + input_model.meta.cal_step["source_detection"] = "COMPLETE" # just pass input model to next step - catalog is stored in meta diff --git a/romancal/source_detection/tests/test_source_detection_step.py b/romancal/source_detection/tests/test_source_detection_step.py index 8ae2317dc..d03280c90 100644 --- a/romancal/source_detection/tests/test_source_detection_step.py +++ b/romancal/source_detection/tests/test_source_detection_step.py @@ -2,220 +2,62 @@ Unit tests for the Roman source detection step code """ -import os +from copy import deepcopy import numpy as np import pytest from astropy import units as u -from astropy.convolution import Gaussian2DKernel -from roman_datamodels import maker_utils as testutil -from roman_datamodels.datamodels import ImageModel +from romancal.lib.tests.test_psf import add_sources, setup_inputs from romancal.source_detection import SourceDetectionStep +n_sources = 10 +image_model_shape = (100, 100) -@pytest.fixture -def setup_inputs(): - def _setup(nrows=100, ncols=100, noise=1.0, seed=None): - """Return ImageModel of lvl 2 image""" - shape = (100, 100) # size of test image - wfi_image = testutil.mk_level2_image(shape=shape) - wfi_image.data = u.Quantity( - np.ones(shape, dtype=np.float32), u.electron / u.s, dtype=np.float32 +class TestSourceDetection: + def setup_method(self): + self.image_model, self.webbpsf_config, self.psf_model = setup_inputs( + shape=image_model_shape, ) - wfi_image.meta.filename = "filename" - # add noise to data - if noise is not None: - rng = np.random.default_rng(seed or 19) - wfi_image.data += u.Quantity( - noise * rng.uniform(size=shape), u.electron / u.s, dtype=np.float32 - ) + @pytest.mark.webbpsf + def test_dao_vs_psf(self, seed=0): + rng = np.random.default_rng(seed) + image_model = deepcopy(self.image_model) - # add dq array + n_models = 10 + amp_true = rng.normal(1e3, 100, n_models) + along_diag = np.arange(100, 950, 90) / 1000 * image_model_shape[0] + x_true = along_diag + rng.normal(scale=0.5, size=n_models) + y_true = along_diag + rng.normal(scale=0.5, size=n_models) - wfi_image.dq = np.zeros(shape, dtype=np.uint32) - # construct ImageModel - mod = ImageModel(wfi_image) + add_sources(image_model, self.psf_model, x_true, y_true, amp_true) - return mod + source_detect = SourceDetectionStep() + source_detect.scalar_threshold = 100 + source_detect.peakmax = None + dao_result = source_detect.process(image_model) + idx, x_dao, y_dao, amp_dao = dao_result.meta.source_detection.tweakreg_catalog - return _setup + # check that all injected targets are found by DAO: + assert len(x_dao) == len(x_true) + source_detect.fit_psf = True + psf_result = source_detect.process(image_model) + psf_catalog = psf_result.meta.source_detection.psf_catalog -def add_random_gauss( - arr, x_positions, y_positions, min_amp=200, max_amp=500, seed=None -): - """Add random 2D Gaussians to `arr` at specified positions, - with random amplitudes from `min_amp` to `max_amp`. Assumes - units of e-/s.""" + extract_columns = ["x_fit", "x_err", "y_fit", "y_err", "flux_fit"] + x_psf, x_err, y_psf, y_err, amp_psf = psf_catalog[extract_columns].itercols() - rng = np.random.default_rng(seed or 29) + # check that typical PSF centroids are more accurate than DAO centroids: + assert np.median(np.abs(x_dao - x_true)) > np.median(np.abs(x_psf - x_true)) - for i, x in enumerate(x_positions): - y = y_positions[i] - gauss = Gaussian2DKernel(2, x_size=21, y_size=21).array - amp = rng.integers(200, 700) - arr[y - 10 : y + 11, x - 10 : x + 11] += ( - u.Quantity(gauss, u.electron / u.s, dtype=np.float32) * amp - ) - - -def test_source_detection_defaults(setup_inputs): - """Test SourceDetectionStep with its default parameters. The detection - threshold will be chosen based on the image's background level.""" - - model = setup_inputs(seed=0) - - # add in 12 sources, roughly evenly distributed - # sort by true_x so they can be matched up to output - - true_x = np.array([20, 25, 30, 35, 40, 45, 50, 55, 60, 70, 80, 82]) - true_y = np.array([26, 80, 44, 19, 66, 39, 29, 72, 54, 29, 80, 62]) - - # at each position, place a 2d gaussian - # randomly vary flux from 100 to 500 for each source - add_random_gauss(model.data, true_x, true_y, seed=0) - - # call SourceDetectionStep with default parameters - sd = SourceDetectionStep() - res = sd.process(model) - - # unpack output catalog array - _, xcentroid, ycentroid, flux = res.meta.source_detection.tweakreg_catalog - - # sort based on x coordinate, like input - ycentroid = [x for y, x in sorted(zip(xcentroid, ycentroid))] - flux = [x for y, x in sorted(zip(xcentroid, flux))] - xcentroid = sorted(xcentroid) - - # check that the number of input and output sources are the same - assert len(flux) == len(true_x) - - # check that their locations agree - # atol=0.2 seems to be the lowest safe value for this right now. - - assert np.allclose(np.abs(xcentroid - true_x), 0.0, atol=0.25) - assert np.allclose(np.abs(ycentroid - true_y), 0.0, atol=0.25) - - -def test_source_detection_scalar_threshold(setup_inputs): - """Test SourceDetectionStep using the option to choose a detection - threshold for entire image.""" - - model = setup_inputs(seed=0) - - # add in 12 sources, roughly evenly distributed - # sort by true_x so they can be matched up to output - - true_x = np.array([20, 25, 30, 35, 40, 45, 50, 55, 60, 70, 80, 82]) - true_y = np.array([26, 80, 44, 19, 66, 39, 29, 72, 54, 29, 80, 62]) - - # at each position, place a 2d gaussian - # randomly vary flux from 100 to 500 for each source - add_random_gauss(model.data, true_x, true_y, seed=0) - - # call SourceDetectionStep with default parameters - sd = SourceDetectionStep() - sd.scalar_threshold = 2.0 - res = sd.process(model) - - # unpack output catalog array - _, xcentroid, ycentroid, flux = res.meta.source_detection.tweakreg_catalog - - # sort based on x coordinate, like input - ycentroid = [x for y, x in sorted(zip(xcentroid, ycentroid))] - flux = [x for y, x in sorted(zip(xcentroid, flux))] - xcentroid = sorted(xcentroid) - - # check that the number of input and output sources are the same - assert len(flux) == len(true_x) - - # check that their locations agree - # atol=0.2 seems to be the lowest safe value for this right now. - - assert np.allclose(np.abs(xcentroid - true_x), 0.0, atol=0.25) - assert np.allclose(np.abs(ycentroid - true_y), 0.0, atol=0.25) - - -@pytest.mark.skipif( - os.environ.get("CI") == "true", - reason="Roman CRDS servers are not currently available outside the internal " - "network", -) -def test_outputs(tmp_path, setup_inputs): - """Make sure `save_catalogs` and `output_cat_filetype` work correctly.""" - - try: - # The contextlib.chdir() context manager was added in Python 3.11 - # so once we drop support for Python 3.10 we can remove this try/except - from contextlib import chdir - - context = chdir(tmp_path) - except ImportError: - # This reproduces enough of the behavior of contextlib.chdir() - # for Python < 3.11 (this part can be removed at a later date) - from contextlib import contextmanager - from pathlib import Path - - @contextmanager - def ctx(): - cwd = Path.cwd() - os.chdir(tmp_path) - try: - yield - finally: - os.chdir(cwd) - - context = ctx() - - with context: - model = setup_inputs(seed=0) - - # add a single source to image so a non-empty catalog is produced - add_random_gauss(model.data, [50], [50], seed=0) - - # run step and direct it to save catalog. default format should be asdf - sd = SourceDetectionStep() - sd.save_catalogs = True - sd.process(model) - # make sure file exists - expected_output_path = os.path.join(tmp_path, "filename_tweakreg_catalog.asdf") - assert os.path.isfile(expected_output_path) - - # run again, specifying that catalog should be saved in ecsv format - sd = SourceDetectionStep() - sd.save_catalogs = True - sd.output_cat_filetype = "ecsv" - sd.process(model) - expected_output_path = os.path.join(tmp_path, "filename_tweakreg_catalog.ecsv") - assert os.path.isfile(expected_output_path) - - -def test_limiting_catalog_size(setup_inputs): - """Test to make sure setting `max_sources` limits the size of the - output catalog to contain only the N brightest sources""" - - model = setup_inputs(seed=0) - - amps = [200, 300, 400] # flux - pos = [20, 50, 80] # 3 sources in a line - for i in range(3): - xy = pos[i] - gauss = Gaussian2DKernel(2, x_size=20, y_size=20).array - model.data[xy - 10 : xy + 10, xy - 10 : xy + 10] += ( - u.Quantity(gauss, u.electron / u.s, dtype=np.float32) * amps[i] - ) - - sd = SourceDetectionStep() - sd.max_sources = 2 - res = sd.process(model) - - _, xcentroid, ycentroid, flux = res.meta.source_detection.tweakreg_catalog - - # make sure only 2 of the three sources are returned in output catalog - assert len(flux) == 2 + # check that the typical/worst PSF centroid is still within some tolerance: + pixel_scale = 0.11 * u.arcsec / u.pix + centroid_residuals = np.abs(x_psf - x_true) * u.pix * pixel_scale + assert np.max(centroid_residuals) < 11 * u.mas + assert np.median(centroid_residuals) < 3 * u.mas - # and make sure one of them is not the first, dimmest source - assert np.all(xcentroid > 22) # dimmest position is really 20, give a + # check that typical residuals are consistent with their errors: + assert np.median(np.abs(x_psf - x_true) / x_err) < 2