diff --git a/python/lsst/pipe/tasks/calibrateImage.py b/python/lsst/pipe/tasks/calibrateImage.py index 95961b961..2fd99f04c 100644 --- a/python/lsst/pipe/tasks/calibrateImage.py +++ b/python/lsst/pipe/tasks/calibrateImage.py @@ -26,13 +26,19 @@ import requests import os +from lsst.afw.geom import SpanSet import lsst.afw.table as afwTable import lsst.afw.image as afwImage +import lsst.afw.math as afwMath from lsst.ip.diffim.utils import evaluateMaskFraction, populate_sattle_visit_cache import lsst.meas.algorithms import lsst.meas.algorithms.installGaussianPsf import lsst.meas.algorithms.measureApCorr import lsst.meas.algorithms.setPrimaryFlags +from lsst.meas.algorithms.adaptive_thresholds import ( + AdaptiveThresholdDetectionTask, + AdaptiveThresholdBackgroundTask, +) import lsst.meas.base import lsst.meas.astrom import lsst.meas.deblender @@ -283,9 +289,14 @@ class CalibrateImageConfig(pipeBase.PipelineTaskConfig, pipelineConnections=Cali doc="Task to perform intial background subtraction, before first detection pass.", ) psf_detection = pexConfig.ConfigurableField( - target=lsst.meas.algorithms.SourceDetectionTask, + target=AdaptiveThresholdDetectionTask, doc="Task to detect sources for PSF determination." ) + do_adaptive_threshold_detection = pexConfig.Field( + dtype=bool, + default=True, + doc="Implement the adaptive detection thresholding approach?", + ) psf_source_measurement = pexConfig.ConfigurableField( target=lsst.meas.base.SingleFrameMeasurementTask, doc="Task to measure sources to be used for psf estimation." @@ -300,14 +311,18 @@ class CalibrateImageConfig(pipeBase.PipelineTaskConfig, pipelineConnections=Cali "for the bright stars used for psf estimation.", ) - # TODO DM-39203: we can remove aperture correction from this task once we are - # using the shape-based star/galaxy code. + # TODO DM-39203: we can remove aperture correction from this task once we + # are using the shape-based star/galaxy code. measure_aperture_correction = pexConfig.ConfigurableField( target=lsst.meas.algorithms.measureApCorr.MeasureApCorrTask, doc="Task to compute the aperture correction from the bright stars." ) # subtasks used during star measurement + star_background = pexConfig.ConfigurableField( + target=AdaptiveThresholdBackgroundTask, + doc="Task to perform final background subtraction, just before photoCal.", + ) star_detection = pexConfig.ConfigurableField( target=lsst.meas.algorithms.SourceDetectionTask, doc="Task to detect stars to return in the output catalog." @@ -451,25 +466,19 @@ def setDefaults(self): self.install_simple_psf.fwhm = 4 # S/N>=50 sources for PSF determination, but detection to S/N=10. - # The thresholdValue sets the minimum flux in a pixel to be included in the - # footprint, while peaks are only detected when they are above + # The thresholdValue sets the minimum flux in a pixel to be included + # in the footprint, while peaks are only detected when they are above # thresholdValue * includeThresholdMultiplier. The low thresholdValue # ensures that the footprints are large enough for the noise replacer # to mask out faint undetected neighbors that are not to be measured. - self.psf_detection.thresholdValue = 10.0 - self.psf_detection.includeThresholdMultiplier = 5.0 - # TODO investigation: Probably want False here, but that may require - # tweaking the background spatial scale, to make it small enough to - # prevent extra peaks in the wings of bright objects. - self.psf_detection.doTempLocalBackground = False - # NOTE: we do want reEstimateBackground=True in psf_detection, so that - # each measurement step is done with the best background available. + self.psf_detection.baseline.thresholdValue = 10.0 + self.psf_detection.baseline.includeThresholdMultiplier = 5.0 # Minimal measurement plugins for PSF determination. # TODO DM-39203: We can drop GaussianFlux and PsfFlux, if we use # shapeHSM/moments for star/galaxy separation. - # TODO DM-39203: we can remove aperture correction from this task once - # we are using the shape-based star/galaxy code. + # TODO DM-39203: we can remove aperture correction from this task + # once we are using the shape-based star/galaxy code. self.psf_source_measurement.plugins = ["base_PixelFlags", "base_SdssCentroid", "ext_shapeHSM_HsmSourceMoments", @@ -494,6 +503,7 @@ def setDefaults(self): # Detection for good S/N for astrometry/photometry and other # downstream tasks; detection mask to S/N>=5, but S/N>=10 peaks. + self.star_detection.reEstimateBackground = False self.star_detection.thresholdValue = 5.0 self.star_detection.includeThresholdMultiplier = 2.0 self.star_measurement.plugins = ["base_PixelFlags", @@ -512,9 +522,9 @@ def setDefaults(self): self.star_measurement.plugins["base_CircularApertureFlux"].radii = [12.0] self.star_measurement.plugins["base_CompensatedTophatFlux"].apertures = [12] - # We measure and apply the normalization aperture correction with the - # psf_normalized_calibration_flux task, and we only apply the normalization - # aperture correction for the full list of stars. + # We measure and apply the normalization aperture correction with + # the psf_normalized_calibration_flux task, and we only apply the + # normalization aperture correction for the full list of stars. self.star_normalized_calibration_flux.do_measure_ap_corr = False # Select stars with reliable measurements and no bad flags. @@ -536,8 +546,8 @@ def setDefaults(self): # Only reject sky sources; we already selected good stars. self.astrometry.sourceSelector["science"].doFlags = True - self.astrometry.sourceSelector["science"].flags.good = ["calib_psf_candidate"] - self.astrometry.sourceSelector["science"].flags.bad = [] + self.astrometry.sourceSelector["science"].flags.good = [] # ["calib_psf_candidate"] + # self.astrometry.sourceSelector["science"].flags.bad = [] self.astrometry.sourceSelector["science"].doUnresolved = False self.astrometry.sourceSelector["science"].doIsolated = False self.astrometry.sourceSelector["science"].doRequirePrimary = False @@ -606,7 +616,7 @@ def validate(self): "CalibrateImageTask.psf_subtract_background must be configured with " "doApplyFlatBackgroundRatio=True if do_illumination_correction=True." ) - if self.psf_detection.reEstimateBackground: + if getattr(self.psf_detection, "reEstimateBackground", False): if not self.psf_detection.doApplyFlatBackgroundRatio: raise pexConfig.FieldValidationError( CalibrateImageConfig.psf_detection, @@ -614,7 +624,7 @@ def validate(self): "CalibrateImageTask.psf_detection background must be configured with " "doApplyFlatBackgroundRatio=True if do_illumination_correction=True." ) - if self.star_detection.reEstimateBackground: + if getattr(self.star_detection, "reEstimateBackground", False): if not self.star_detection.doApplyFlatBackgroundRatio: raise pexConfig.FieldValidationError( CalibrateImageConfig.star_detection, @@ -627,6 +637,26 @@ def validate(self): if not os.getenv("SATTLE_URI_BASE"): raise pexConfig.FieldValidationError(CalibrateImageConfig.run_sattle, self, "Sattle requested but URI environment variable not set.") + if ( + issubclass(self.psf_detection.target, AdaptiveThresholdDetectionTask) + != self.do_adaptive_threshold_detection + ): + raise pexConfig.FieldValidationError( + CalibrateImageConfig.psf_detection, self, + "AdaptiveThresholdDetectionTask must be used for the psf_detection subtask " + "if and only if do_adaptive_threshold_detection is True." + ) + + if not self.do_adaptive_threshold_detection: + if not self.star_detection.reEstimateBackground: + raise pexConfig.FieldValidationError( + CalibrateImageConfig.psf_detection, + self, + "If not using the adaptive threshold detection implementation " + "(i.e. config.do_adaptive_threshold_detection = False), " + "CalibrateImageTask.star_detection background must be configured " + "with reEstimateBackground = True to maintain original behavior." + ) class CalibrateImageTask(pipeBase.PipelineTask): @@ -673,7 +703,8 @@ def __init__(self, initial_stars_schema=None, **kwargs): # astrometric fitting, and aperture correction calculations. self.psf_fields = ("calib_psf_candidate", "calib_psf_used", "calib_psf_reserved", "calib_astrometry_used", - # TODO DM-39203: these can be removed once apcorr is gone. + # TODO DM-39203: + # these can be removed once apcorr is gone. "apcorr_slot_CalibFlux_used", "apcorr_base_GaussianFlux_used", "apcorr_base_PsfFlux_used",) for field in self.psf_fields: @@ -689,6 +720,7 @@ def __init__(self, initial_stars_schema=None, **kwargs): doc="Maximum value in the star image used to train PSF.") afwTable.CoordKey.addErrorFields(initial_stars_schema) + self.makeSubtask("star_background") self.makeSubtask("star_detection", schema=initial_stars_schema) self.makeSubtask("star_sky_sources", schema=initial_stars_schema) self.makeSubtask("star_deblend", schema=initial_stars_schema) @@ -806,7 +838,8 @@ def run( Parameters ---------- - exposures : `lsst.afw.image.Exposure` or `list` [`lsst.afw.image.Exposure`] + exposures : `lsst.afw.image.Exposure` or \ + `list` [`lsst.afw.image.Exposure`] Post-ISR exposure(s), with an initial WCS, VisitInfo, and Filter. Modified in-place during processing if only one is passed. If two exposures are passed, treat them as snaps and combine @@ -854,10 +887,12 @@ def run( This is `None` if ``config.do_calibrate_pixels`` is `False`. ``astrometry_matches`` Reference catalog stars matches used in the astrometric fit. - (`list` [`lsst.afw.table.ReferenceMatch`] or `lsst.afw.table.BaseCatalog`) + (`list` [`lsst.afw.table.ReferenceMatch`] or + `lsst.afw.table.BaseCatalog`). ``photometry_matches`` Reference catalog stars matches used in the photometric fit. - (`list` [`lsst.afw.table.ReferenceMatch`] or `lsst.afw.table.BaseCatalog`) + (`list` [`lsst.afw.table.ReferenceMatch`] or + `lsst.afw.table.BaseCatalog`). ``mask`` Copy of the mask plane of `exposure`. (`lsst.afw.image.Mask`) @@ -906,11 +941,12 @@ def run( illumination_correction, ) - result.psf_stars_footprints, result.background, _ = self._compute_psf( + psf_detections, result.background, _ = self._compute_psf( result.exposure, id_generator, background_to_photometric_ratio=result.background_to_photometric_ratio, ) + result.psf_stars_footprints = psf_detections.sources have_fit_psf = True # Check if all centroids have been flagged. This should happen @@ -926,32 +962,50 @@ def run( psf_size=psf_shape.getDeterminantRadius(), ) + if result.background is None: + result.background = afwMath.BackgroundList() + self._measure_aperture_correction(result.exposure, result.psf_stars_footprints) result.psf_stars = result.psf_stars_footprints.asAstropy() - # Run astrometry using PSF candidate stars + # Run astrometry using PSF candidate stars. + # Update "the psf_stars" source cooordinates with the current wcs. + afwTable.updateSourceCoords( + result.exposure.wcs, + sourceList=result.psf_stars_footprints, + include_covariance=self.config.do_include_astrometric_errors, + ) astrometry_matches, astrometry_meta = self._fit_astrometry( result.exposure, result.psf_stars_footprints ) - self.metadata["astrometry_matches_count"] = len(astrometry_matches) + num_astrometry_matches = len(astrometry_matches) + self.metadata["astrometry_matches_count"] = num_astrometry_matches if "astrometry_matches" in self.config.optional_outputs: result.astrometry_matches = lsst.meas.astrom.denormalizeMatches(astrometry_matches, astrometry_meta) result.psf_stars = result.psf_stars_footprints.asAstropy() + if self.config.do_adaptive_threshold_detection: + self._remeasure_star_background(result) + # Run the stars_detection subtask for the photometric calibration. result.stars_footprints = self._find_stars( result.exposure, result.background, id_generator, background_to_photometric_ratio=result.background_to_photometric_ratio, + psf_detections=psf_detections, + num_astrometry_matches=num_astrometry_matches, ) - self._match_psf_stars(result.psf_stars_footprints, result.stars_footprints) + psf = result.exposure.getPsf() + psfSigma = psf.computeShape(result.exposure.getBBox().getCenter()).getDeterminantRadius() + self._match_psf_stars(result.psf_stars_footprints, result.stars_footprints, + psfSigma=psfSigma) - # Update the source cooordinates with the current wcs. + # Update the "stars" source cooordinates with the current wcs. afwTable.updateSourceCoords( result.exposure.wcs, sourceList=result.stars_footprints, - include_covariance=self.config.do_include_astrometric_errors + include_covariance=self.config.do_include_astrometric_errors, ) summary_stat_catalog = result.stars_footprints @@ -969,7 +1023,8 @@ def run( photometry_meta, photo_calib = self._fit_photometry(result.exposure, result.stars_footprints) have_fit_photometry = True self.metadata["photometry_matches_count"] = len(photometry_matches) - # fit_photometry returns a new catalog, so we need a new astropy table view. + # fit_photometry returns a new catalog, so we need a new astropy + # table view. result.stars = result.stars_footprints.asAstropy() # summary stats don't make use of the calibrated fluxes, but we # might as well use the best catalog we've got in case that @@ -989,12 +1044,12 @@ def run( result.exposure.setWcs(None) if not have_fit_photometry: result.exposure.setPhotoCalib(None) - # Summary stat calculations can handle missing components gracefully, - # but we want to run them as late as possible (but still before we - # calibrate pixels, if we do that at all). - # So we run them after we succeed or if we get an AlgorithmError. We - # intentionally don't use 'finally' here because we don't want to run - # them if we get some other kind of error. + # Summary stat calculations can handle missing components + # gracefully, but we want to run them as late as possible (but + # still before wecalibrate pixels, if we do that at all). + # So we run them after we succeed or if we get an AlgorithmError. + # We intentionally don't use 'finally' here because we don't + # want to run them if we get some other kind of error. self._summarize(result.exposure, summary_stat_catalog, result.background) raise else: @@ -1032,7 +1087,8 @@ def _apply_illumination_correction(self, exposure, background_flat, illumination background_flat : `lsst.afw.image.Exposure` Flat image that had previously been applied to exposure. illumination_correction : `lsst.afw.image.Exposure` - Illumination correction image to convert to photometric-flattened image. + Illumination correction image to convert to photometric-flattened + image. Returns ------- @@ -1085,8 +1141,10 @@ def _compute_psf(self, exposure, id_generator, background_to_photometric_ratio=N Returns ------- - sources : `lsst.afw.table.SourceCatalog` - Catalog of detected bright sources. + detections : `lsst.pipe.base.Struct` + Struct returned by the ``psf_detection`` subtask, with its + ``sources`` attribute further updated by measurement and PSF + modeling. background : `lsst.afw.math.BackgroundList` Background that was fit to the exposure during detection. cell_set : `lsst.afw.math.SpatialCellSet` @@ -1099,16 +1157,19 @@ def log_psf(msg, addToMetadata=False): Parameters ---------- - msg : `str` - Message to prepend the log info with. - addToMetadata : `bool`, optional - Whether to add the final psf sigma value to the task metadata - (the default is False). + msg : `str` + Message to prepend the log info with. + addToMetadata : `bool`, optional + Whether to add the final psf sigma value to the task + metadata (the default is False). """ position = exposure.psf.getAveragePosition() sigma = exposure.psf.computeShape(position).getDeterminantRadius() dimensions = exposure.psf.computeImage(position).getDimensions() - median_background = np.median(background.getImage().array) + if background is not None: + median_background = np.median(background.getImage().array) + else: + median_background = float("nan") self.log.info("%s sigma=%0.4f, dimensions=%s; median background=%0.2f", msg, sigma, dimensions, median_background) if addToMetadata: @@ -1126,48 +1187,51 @@ def log_psf(msg, addToMetadata=False): self.psf_repair.run(exposure=exposure, keepCRs=True) table = afwTable.SourceTable.make(self.psf_schema, id_generator.make_table_id_factory()) - # Re-estimate the background during this detection step, so that - # measurement uses the most accurate background-subtraction. detections = self.psf_detection.run( table=table, exposure=exposure, background=background, backgroundToPhotometricRatio=background_to_photometric_ratio, ) + self.metadata["initial_psf_positive_footprint_count"] = detections.numPos self.metadata["initial_psf_negative_footprint_count"] = detections.numNeg self.metadata["initial_psf_positive_peak_count"] = detections.numPosPeaks self.metadata["initial_psf_negative_peak_count"] = detections.numNegPeaks self.psf_source_measurement.run(detections.sources, exposure) psf_result = self.psf_measure_psf.run(exposure=exposure, sources=detections.sources) - # Replace the initial PSF with something simpler for the second - # repair/detect/measure/measure_psf step: this can help it converge. - self.install_simple_psf.run(exposure=exposure) - log_psf("Rerunning with simple PSF:") - # TODO investigation: Should we only re-run repair here, to use the - # new PSF? Maybe we *do* need to re-run measurement with PsfFlux, to - # use the fitted PSF? - # TODO investigation: do we need a separate measurement task here - # for the post-psf_measure_psf step, since we only want to do PsfFlux - # and GaussianFlux *after* we have a PSF? Maybe that's not relevant - # once DM-39203 is merged? - self.psf_repair.run(exposure=exposure, keepCRs=True) - # Re-estimate the background during this detection step, so that - # measurement uses the most accurate background-subtraction. - detections = self.psf_detection.run( - table=table, - exposure=exposure, - background=background, - backgroundToPhotometricRatio=background_to_photometric_ratio, - ) + # DO WE REALLY NEED THE 2nd ROUND OF PSF FITTING?? + if not self.config.do_adaptive_threshold_detection: + # Replace the initial PSF with something simpler for the second + # repair/detect/measure/measure_psf step: this can help it + # converge. + self.install_simple_psf.run(exposure=exposure) + + log_psf("Rerunning with simple PSF:") + # TODO investigation: Should we only re-run repair here, to use the + # new PSF? Maybe we *do* need to re-run measurement with PsfFlux, + # to use the fitted PSF? + # TODO investigation: do we need a separate measurement task here + # for the post-psf_measure_psf step, since we only want to do + # PsfFlux and GaussianFlux *after* we have a PSF? Maybe that's not + # relevant once DM-39203 is merged? + self.psf_repair.run(exposure=exposure, keepCRs=True) + # Re-estimate the background during this detection step, so that + # measurement uses the most accurate background-subtraction. + detections = self.psf_detection.run( + table=table, + exposure=exposure, + background=background, + backgroundToPhotometricRatio=background_to_photometric_ratio, + ) + self.psf_source_measurement.run(detections.sources, exposure) + psf_result = self.psf_measure_psf.run(exposure=exposure, sources=detections.sources) + self.metadata["simple_psf_positive_footprint_count"] = detections.numPos self.metadata["simple_psf_negative_footprint_count"] = detections.numNeg self.metadata["simple_psf_positive_peak_count"] = detections.numPosPeaks self.metadata["simple_psf_negative_peak_count"] = detections.numNegPeaks - self.psf_source_measurement.run(detections.sources, exposure) - psf_result = self.psf_measure_psf.run(exposure=exposure, sources=detections.sources) - log_psf("Final PSF:", addToMetadata=True) # Final repair with final PSF, removing cosmic rays this time. @@ -1177,7 +1241,7 @@ def log_psf(msg, addToMetadata=False): # PSF is set on exposure; candidates are returned to use for # calibration flux normalization and aperture corrections. - return detections.sources, background, psf_result.cellSet + return detections, background, psf_result.cellSet def _measure_aperture_correction(self, exposure, bright_sources): """Measure and set the ApCorrMap on the Exposure, using @@ -1193,8 +1257,8 @@ def _measure_aperture_correction(self, exposure, bright_sources): Exposure to set the ApCorrMap on. bright_sources : `lsst.afw.table.SourceCatalog` Catalog of detected bright sources; modified to include columns - necessary for point source determination for the aperture correction - calculation. + necessary for point source determination for the aperture + correction calculation. """ norm_ap_corr_map = self.psf_normalized_calibration_flux.run( exposure=exposure, @@ -1209,7 +1273,8 @@ def _measure_aperture_correction(self, exposure, bright_sources): exposure.info.setApCorrMap(ap_corr_map) - def _find_stars(self, exposure, background, id_generator, background_to_photometric_ratio=None): + def _find_stars(self, exposure, background, id_generator, background_to_photometric_ratio=None, *, + psf_detections, num_astrometry_matches): """Detect stars on an exposure that has a PSF model, and measure their PSF, circular aperture, compensated gaussian fluxes. @@ -1225,6 +1290,10 @@ def _find_stars(self, exposure, background, id_generator, background_to_photomet background_to_photometric_ratio : `lsst.afw.image.Image`, optional Image to convert photometric-flattened image to background-flattened image. + psf_detections : `lsst.pipe.base.Struct` + Struct returned from the ``psf_detection`` subtask. + num_astrometry_matches : `int` + Number of astrometry matches. Returns ------- @@ -1234,14 +1303,68 @@ def _find_stars(self, exposure, background, id_generator, background_to_photomet """ table = afwTable.SourceTable.make(self.initial_stars_schema.schema, id_generator.make_table_id_factory()) - # Re-estimate the background during this detection step, so that - # measurement uses the most accurate background-subtraction. - detections = self.star_detection.run( - table=table, - exposure=exposure, - background=background, - backgroundToPhotometricRatio=background_to_photometric_ratio, - ) + + inAdaptiveDet = True + maxAdaptiveDetIter = 8 + adaptiveDetIter = 0 + threshFactor = 0.2 + if self.config.do_adaptive_threshold_detection: + while inAdaptiveDet and adaptiveDetIter < maxAdaptiveDetIter: + inAdaptiveDet = False + adaptiveDetectionConfig = lsst.meas.algorithms.SourceDetectionConfig() + adaptiveDetectionConfig.reEstimateBackground = False + adaptiveDetectionConfig.includeThresholdMultiplier = 2.0 + psfThreshold = ( + psf_detections.thresholdValue*psf_detections.includeThresholdMultiplier + ) + adaptiveDetectionConfig.thresholdValue = max( + self.config.star_detection.thresholdValue, + threshFactor*psfThreshold/adaptiveDetectionConfig.includeThresholdMultiplier + ) + self.log.info("Using adaptive threshold detection (nIter = %d) with " + "thresholdValue = %.2f and multiplier = %.1f", + adaptiveDetIter, adaptiveDetectionConfig.thresholdValue, + adaptiveDetectionConfig.includeThresholdMultiplier) + adaptiveDetectionTask = lsst.meas.algorithms.SourceDetectionTask( + config=adaptiveDetectionConfig + ) + detections = adaptiveDetectionTask.run( + table=table, + exposure=exposure, + background=background, + backgroundToPhotometricRatio=background_to_photometric_ratio, + ) + nFootprint = len(detections.sources) + nPeak = 0 + nIsolated = 0 + for src in detections.sources: + nPeakSrc = len(src.getFootprint().getPeaks()) + if nPeakSrc == 1: + nIsolated += 1 + nPeak += nPeakSrc + minIsolated = min(400, max(3, 0.005*nPeak, 0.6*num_astrometry_matches)) + if nFootprint > 0: + self.log.info("Adaptive threshold detection _find_stars nIter: %d; " + "nPeak/nFootprint = %.2f (max is 800); nIsolated = %d (min is %.1f).", + adaptiveDetIter, nPeak/nFootprint, nIsolated, minIsolated) + if nPeak/nFootprint > 800 or nIsolated < minIsolated: + threshFactor = max(0.01, 1.5*threshFactor) + inAdaptiveDet = True + self.log.warning("nPeak/nFootprint = %.2f (max is 800); nIsolated = %d " + "(min is %.1f).", nPeak/nFootprint, nIsolated, minIsolated) + else: + self.log.warning("No footprints detected on image and rerunning...") + inAdaptiveDet = True + adaptiveDetIter += 1 + else: + # Re-estimate the background during this detection step, so that + # measurement uses the most accurate background-subtraction. + detections = self.star_detection.run( + table=table, + exposure=exposure, + background=background, + backgroundToPhotometricRatio=background_to_photometric_ratio, + ) sources = detections.sources self.star_sky_sources.run(exposure.mask, id_generator.catalog_id, sources) @@ -1272,7 +1395,8 @@ def _find_stars(self, exposure, background, id_generator, background_to_photomet sel[indices] = True sources = sources[sel] - # TODO investigation: Could this deblender throw away blends of non-PSF sources? + # TODO investigation: Could this deblender throw away blends of + # non-PSF sources? self.star_deblend.run(exposure=exposure, sources=sources) # The deblender may not produce a contiguous catalog; ensure # contiguity for subsequent tasks. @@ -1300,7 +1424,7 @@ def _find_stars(self, exposure, background, id_generator, background_to_photomet else: return result.sourceCat - def _match_psf_stars(self, psf_stars, stars): + def _match_psf_stars(self, psf_stars, stars, psfSigma=None): """Match calibration stars to psf stars, to identify which were psf candidates, and which were used or reserved during psf measurement and the astrometric fit. @@ -1320,9 +1444,10 @@ def _match_psf_stars(self, psf_stars, stars): This code was adapted from CalibrateTask.copyIcSourceFields(). """ control = afwTable.MatchControl() + matchRadius = 3.0 if psfSigma is None else max(3.0, psfSigma) # in pixels # Return all matched objects, to separate blends. control.findOnlyClosest = False - matches = afwTable.matchXy(psf_stars, stars, 3.0, control) + matches = afwTable.matchXy(psf_stars, stars, matchRadius, control) deblend_key = stars.schema["deblend_nChild"].asKey() matches = [m for m in matches if m[1].get(deblend_key) == 0] @@ -1445,7 +1570,8 @@ def _apply_photometry(self, exposure, background, background_to_photometric_rati "Background calibration assumes a constant PhotoCalib; PhotoCalTask should always return that." for bg in background: - # The statsImage is a view, but we can't assign to a function call in python. + # The statsImage is a view, but we can't assign to a function call + # in python. binned_image = bg[0].getStatsImage() binned_image *= photo_calib.getCalibrationMean() @@ -1501,3 +1627,62 @@ def _update_wcs_with_camera_model(self, exposure, cameraModel): detector = cameraModel[exposure.detector.getId()] newWcs = createInitialSkyWcs(exposure.visitInfo, detector) exposure.setWcs(newWcs) + + def _remeasure_star_background(self, result): + """Remeasure the exposure's background with iterative adaptive + threshold detection. + + Parameters + ---------- + result : `lsst.pipe.base.Struct` + The current state of the result Strut from the run method of + CalibrateImageTask. Will be modified in place. + + Returns + ------- + result : `lsst.pipe.base.Struct` + The modified result Struct with the new background subtracted. + """ + result.background = self.star_background.run(result.exposure, background=result.background).background + + # Perform one more round of background subtraction that is just an + # overall pedestal (order = 0). This is intended to account for + # any potential gross oversubtraction imposed by the higher-order + # subtraction. + # Dilate DETECTED mask a bit more if it's below 50% detected. + nPixToDilate = 2 + if detected_fraction < 0.5: + dilatedMask = result.exposure.mask.clone() + for maskName in detected_mask_planes: + # Compute the grown detection mask plane using SpanSet + detectedMaskBit = dilatedMask.getPlaneBitMask(maskName) + detectedMaskSpanSet = SpanSet.fromMask(dilatedMask, detectedMaskBit) + detectedMaskSpanSet = detectedMaskSpanSet.dilated(nPixToDilate) + detectedMaskSpanSet = detectedMaskSpanSet.clippedTo(dilatedMask.getBBox()) + # Clear the detected mask plane + detectedMask = dilatedMask.getMaskPlane(maskName) + dilatedMask.clearMaskPlane(detectedMask) + # Set the mask plane to the dilated one + detectedMaskSpanSet.setMask(dilatedMask, detectedMaskBit) + + detected_fraction_dilated = self._compute_mask_fraction(dilatedMask, + detected_mask_planes, + bad_mask_planes) + result.exposure.mask = dilatedMask + self.log.debug("detected_fraction_orig = %.3f detected_fraction_dilated = %.3f", + detected_fraction_orig, detected_fraction_dilated) + + pedestalBackgroundConfig = lsst.meas.algorithms.SubtractBackgroundConfig() + pedestalBackgroundConfig.statisticsProperty = "MEDIAN" + pedestalBackgroundConfig.approxOrderX = 0 + pedestalBackgroundConfig.binSize = 64 + pedestalBackgroundTask = lsst.meas.algorithms.SubtractBackgroundTask(config=pedestalBackgroundConfig) + pedestalBackgroundList = pedestalBackgroundTask.run( + exposure=result.exposure, background=result.background).background + # Isolate the final pedestal background to log the computed value + pedestalBackground = afwMath.BackgroundList() + pedestalBackground.append(pedestalBackgroundList[1]) + pedestalBackgroundLevel = pedestalBackground.getImage().array[0, 0] + self.log.warning("Subtracted pedestalBackgroundLevel = %.4f", pedestalBackgroundLevel) + + return result diff --git a/python/lsst/pipe/tasks/multiBand.py b/python/lsst/pipe/tasks/multiBand.py index ed3b6ac35..be286d9a6 100644 --- a/python/lsst/pipe/tasks/multiBand.py +++ b/python/lsst/pipe/tasks/multiBand.py @@ -53,8 +53,9 @@ from lsst.meas.extensions.scarlet.io import updateCatalogFootprints from lsst.meas.astrom import DirectMatchTask, denormalizeMatches from lsst.pipe.tasks.propagateSourceFlags import PropagateSourceFlagsTask -import lsst.afw.table as afwTable +import lsst.afw.image as afwImage import lsst.afw.math as afwMath +import lsst.afw.table as afwTable from lsst.daf.base import PropertyList from lsst.skymap import BaseSkyMap @@ -315,6 +316,17 @@ def run(self, exposure, idFactory, expId, patchInfo=None): if hasattr(detections, "background") and detections.background: for bg in detections.background: backgrounds.append(bg) + if len(backgrounds) == 0: + # Persist a constant background with value of 0.0 to get around + # inability to persist empty BackgroundList. + bgLevel = 0.0 + bgStats = afwImage.MaskedImageF(1, 1) + bgStats.set(bgLevel, 0, bgLevel) + bg = afwMath.BackgroundMI(exposure.getBBox(), bgStats) + bgData = (bg, afwMath.Interpolate.LINEAR, afwMath.REDUCE_INTERP_ORDER, + afwMath.ApproximateControl.UNKNOWN, 0, 0, False) + backgrounds.append(bgData) + return Struct(outputSources=sources, outputBackgrounds=backgrounds, outputExposure=exposure) def _cropToExactBinning(self, exposure, patchInfo): diff --git a/python/lsst/pipe/tasks/photoCal.py b/python/lsst/pipe/tasks/photoCal.py index 60366b503..7c7538271 100644 --- a/python/lsst/pipe/tasks/photoCal.py +++ b/python/lsst/pipe/tasks/photoCal.py @@ -404,7 +404,13 @@ def run(self, exposure, sourceCat, expId=0): filterLabel = exposure.getFilter() # Match sources - matchResults = self.match.run(sourceCat, filterLabel.bandLabel) + if exposure.visitInfo is not None: + epoch = exposure.visitInfo.date.toAstropy() + else: + epoch = None + self.log.warning("visitInfo is None for exposure %d. Setting epoch to None.", expId) + + matchResults = self.match.run(sourceCat, filterLabel.bandLabel, epoch=epoch) matches = matchResults.matches reserveResults = self.reserve.run([mm.second for mm in matches], expId=expId) diff --git a/tests/test_calibrateImage.py b/tests/test_calibrateImage.py index e656f7276..2a9f8a65f 100644 --- a/tests/test_calibrateImage.py +++ b/tests/test_calibrateImage.py @@ -39,6 +39,7 @@ import lsst.daf.butler.tests as butlerTests import lsst.geom import lsst.meas.algorithms +from lsst.meas.algorithms.adaptive_thresholds import AdaptiveThresholdDetectionTask from lsst.meas.algorithms import testUtils import lsst.meas.extensions.psfex import lsst.meas.base @@ -115,6 +116,14 @@ def setUp(self): # Test-specific configuration: self.config = CalibrateImageTask.ConfigClass() + # Maintain original, no adaptive threshold detection, configs values. + self.config.do_adaptive_threshold_detection = False + self.config.psf_detection.retarget(lsst.meas.algorithms.SourceDetectionTask) + self.config.psf_detection.thresholdValue = 10 + self.config.psf_detection.includeThresholdMultiplier = 5 + self.config.psf_detection.doTempLocalBackground = False + self.config.psf_detection.reEstimateBackground = True + self.config.star_detection.reEstimateBackground = True # We don't have many sources, so have to fit simpler models. self.config.psf_detection.background.approxOrderX = 1 self.config.star_detection.background.approxOrderX = 1 @@ -123,6 +132,7 @@ def setUp(self): # Use PCA psf fitter, as psfex fails if there are only 4 stars. self.config.psf_measure_psf.psfDeterminer = 'pca' # We don't have many test points, so can't match on complicated shapes. + # self.config.astrometry.sourceSelector["science"].flags.good = ["calib_psf_candidate"] self.config.astrometry.sourceSelector["science"].flags.good = [] self.config.astrometry.matcher.numPointsForShape = 3 self.config.run_sattle = False @@ -154,7 +164,8 @@ def setUp(self): # Something about this test dataset prefers a larger threshold here. self.config.star_selector["science"].unresolved.maximum = 0.2 - def _check_run(self, calibrate, result, expect_calibrated_pixels: bool = True): + def _check_run(self, calibrate, result, expect_calibrated_pixels: bool = True, + expect_n_background: int = 4): """Test the result of CalibrateImage.run(). Parameters @@ -168,7 +179,7 @@ def _check_run(self, calibrate, result, expect_calibrated_pixels: bool = True): """ # Background should have 4 elements: 3 from compute_psf and one from # re-estimation during source detection. - self.assertEqual(len(result.background), 4) + self.assertEqual(len(result.background), expect_n_background) # Both afw and astropy psf_stars catalogs should be populated. self.assertEqual(result.psf_stars["calib_psf_used"].sum(), 3) @@ -217,13 +228,13 @@ def _check_run(self, calibrate, result, expect_calibrated_pixels: bool = True): # Check that the psf_stars cross match worked correctly. matches = esutil.numpy_util.match(result.psf_stars["id"], result.stars["psf_id"]) self.assertFloatsAlmostEqual(result.psf_stars["slot_Centroid_x"][matches[0]], - result.stars["slot_Centroid_x"][matches[1]], atol=2e-5) + result.stars["slot_Centroid_x"][matches[1]], atol=3e-4) if "astrometry_matches" in self.config.optional_outputs: matches = esutil.numpy_util.match(result.astrometry_matches["src_id"], result.photometry_matches["src_psf_id"]) self.assertFloatsAlmostEqual(result.astrometry_matches["src_slot_Centroid_x"][matches[0]], result.photometry_matches["src_slot_Centroid_x"][matches[1]], - atol=2e-5) + atol=3e-4) def test_run(self): """Test that run() returns reasonable values to be butler put. @@ -235,6 +246,29 @@ def test_run(self): self._check_run(calibrate, result) + def test_run_adaptive_threshold_detection(self): + """Test that run() runs with adaptive threshold detection turned on. + """ + # Set the adaptive threshold detection, config values... + self.config.do_adaptive_threshold_detection = True + self.config.psf_detection.retarget(AdaptiveThresholdDetectionTask) + self.config.psf_detection.baseline.thresholdValue = 10 + self.config.psf_detection.baseline.includeThresholdMultiplier = 5 + self.config.psf_detection.minFootprint = 4 + self.config.psf_detection.minIsolated = 4 + self.config.psf_detection.sufficientIsolated = 4 + self.config.star_detection.reEstimateBackground = False + + calibrate = CalibrateImageTask(config=self.config) + calibrate.astrometry.setRefObjLoader(self.ref_loader) + calibrate.photometry.match.setRefObjLoader(self.ref_loader) + with self.assertLogs("lsst.calibrateImage", level="INFO") as cm: + result = calibrate.run(exposures=self.exposure) + subString = "Using adaptive threshold detection " + self.assertTrue(any(subString in s for s in cm.output)) + + self._check_run(calibrate, result, expect_n_background=2) + def test_run_downsample(self): """Test that run() runs with downsample. """ @@ -313,7 +347,8 @@ def test_compute_psf(self): that a PSF is assigned to the expopsure. """ calibrate = CalibrateImageTask(config=self.config) - psf_stars, background, candidates = calibrate._compute_psf(self.exposure, self.id_generator) + psf_detections, background, _ = calibrate._compute_psf(self.exposure, self.id_generator) + psf_stars = psf_detections.sources # Catalog ids should be very large from this id generator. self.assertTrue(all(psf_stars['id'] > 1000000000)) @@ -385,7 +420,8 @@ def test_measure_aperture_correction(self): exposure. """ calibrate = CalibrateImageTask(config=self.config) - psf_stars, background, candidates = calibrate._compute_psf(self.exposure, self.id_generator) + psf_detections, _, _ = calibrate._compute_psf(self.exposure, self.id_generator) + psf_stars = psf_detections.sources # First check that the exposure doesn't have an ApCorrMap. self.assertIsNone(self.exposure.apCorrMap) @@ -400,10 +436,12 @@ def test_find_stars(self): in the image and returns them in the output catalog. """ calibrate = CalibrateImageTask(config=self.config) - psf_stars, background, candidates = calibrate._compute_psf(self.exposure, self.id_generator) + psf_detections, background, _ = calibrate._compute_psf(self.exposure, self.id_generator) + psf_stars = psf_detections.sources calibrate._measure_aperture_correction(self.exposure, psf_stars) - stars = calibrate._find_stars(self.exposure, background, self.id_generator) + stars = calibrate._find_stars(self.exposure, background, self.id_generator, + psf_detections=psf_detections, num_astrometry_matches=0) # Catalog ids should be very large from this id generator. self.assertTrue(all(stars['id'] > 1000000000)) @@ -428,15 +466,13 @@ def test_astrometry(self): """ calibrate = CalibrateImageTask(config=self.config) calibrate.astrometry.setRefObjLoader(self.ref_loader) - psf_stars, background, candidates = calibrate._compute_psf(self.exposure, self.id_generator) + psf_detections, background, _ = calibrate._compute_psf(self.exposure, self.id_generator) + psf_stars = psf_detections.sources calibrate._measure_aperture_correction(self.exposure, psf_stars) - stars = calibrate._find_stars(self.exposure, background, self.id_generator) - - calibrate._fit_astrometry(self.exposure, stars) + calibrate._fit_astrometry(self.exposure, psf_stars) # Check that we got reliable matches with the truth coordinates. - sky = stars["sky_source"] - fitted = SkyCoord(stars[~sky]['coord_ra'], stars[~sky]['coord_dec'], unit="radian") + fitted = SkyCoord(psf_stars['coord_ra'], psf_stars['coord_dec'], unit="radian") truth = SkyCoord(self.truth_cat['coord_ra'], self.truth_cat['coord_dec'], unit="radian") idx, d2d, _ = fitted.match_to_catalog_sky(truth) np.testing.assert_array_less(d2d.to_value(u.milliarcsecond), 35.0) @@ -448,9 +484,11 @@ def test_photometry(self): calibrate = CalibrateImageTask(config=self.config) calibrate.astrometry.setRefObjLoader(self.ref_loader) calibrate.photometry.match.setRefObjLoader(self.ref_loader) - psf_stars, background, candidates = calibrate._compute_psf(self.exposure, self.id_generator) + psf_detections, background, _ = calibrate._compute_psf(self.exposure, self.id_generator) + psf_stars = psf_detections.sources calibrate._measure_aperture_correction(self.exposure, psf_stars) - stars = calibrate._find_stars(self.exposure, background, self.id_generator) + stars = calibrate._find_stars(self.exposure, background, self.id_generator, + psf_detections=psf_detections, num_astrometry_matches=0) calibrate._fit_astrometry(self.exposure, stars) stars, matches, meta, photoCalib = calibrate._fit_photometry(self.exposure, stars) @@ -489,9 +527,11 @@ def test_match_psf_stars(self): and candidates. """ calibrate = CalibrateImageTask(config=self.config) - psf_stars, background, candidates = calibrate._compute_psf(self.exposure, self.id_generator) + psf_detections, background, _ = calibrate._compute_psf(self.exposure, self.id_generator) + psf_stars = psf_detections.sources calibrate._measure_aperture_correction(self.exposure, psf_stars) - stars = calibrate._find_stars(self.exposure, background, self.id_generator) + stars = calibrate._find_stars(self.exposure, background, self.id_generator, + psf_detections=psf_detections, num_astrometry_matches=0) # There should be no psf-related flags set at first. self.assertEqual(stars["calib_psf_candidate"].sum(), 0) @@ -803,7 +843,6 @@ def test_runQuantum_illumination_correction(self): config = CalibrateImageTask.ConfigClass() config.do_illumination_correction = True config.psf_subtract_background.doApplyFlatBackgroundRatio = True - config.psf_detection.doApplyFlatBackgroundRatio = True config.star_detection.doApplyFlatBackgroundRatio = True task = CalibrateImageTask(config=config) lsst.pipe.base.testUtils.assertValidInitOutput(task)