From 2099758314bf2d5e03d59173af0aebb52159e446 Mon Sep 17 00:00:00 2001 From: Lauren MacArthur Date: Fri, 28 Feb 2025 11:40:29 -0800 Subject: [PATCH 01/12] Add epoch for photoCal refCat loading This is to allow for proper motion corrections. --- python/lsst/pipe/tasks/photoCal.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) 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) From d44331f7aa08a8cd031251e83575c749acefdfc2 Mon Sep 17 00:00:00 2001 From: Lauren MacArthur Date: Sun, 24 Aug 2025 17:40:45 -0700 Subject: [PATCH 02/12] Enable adaptive thresholding detection strategy This implementation is to allow the single frame processing calibration to handle a wide range of "scenes", e.g. crowded and/or high nebulosity ields in addition to the more "typical" regions, by implementing an adaptive and iterative thresholding based on the current detection masks and tuning the thresholds to adapt to the current scene. In general, higher detection thresholds are requird for non-standard fields to avoid large contiguous detection footprints which leave too few good isolated sources for calibration. To achieve the previous behaviour, the adaptive scheme can be turned off by setting the do_adaptive_threshold_detection config to False. --- python/lsst/pipe/tasks/calibrateImage.py | 489 ++++++++++++++++++++--- 1 file changed, 435 insertions(+), 54 deletions(-) diff --git a/python/lsst/pipe/tasks/calibrateImage.py b/python/lsst/pipe/tasks/calibrateImage.py index 95961b961..e9d24f5db 100644 --- a/python/lsst/pipe/tasks/calibrateImage.py +++ b/python/lsst/pipe/tasks/calibrateImage.py @@ -26,8 +26,10 @@ 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 @@ -286,6 +288,15 @@ class CalibrateImageConfig(pipeBase.PipelineTaskConfig, pipelineConnections=Cali target=lsst.meas.algorithms.SourceDetectionTask, 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_adaptive_threshold_detection = pexConfig.ConfigurableField( + target=lsst.meas.algorithms.AdaptiveThresholdDetectionTask, + doc="Task to adaptively detect sources for PSF determination.", + ) psf_source_measurement = pexConfig.ConfigurableField( target=lsst.meas.base.SingleFrameMeasurementTask, doc="Task to measure sources to be used for psf estimation." @@ -308,6 +319,10 @@ class CalibrateImageConfig(pipeBase.PipelineTaskConfig, pipelineConnections=Cali ) # subtasks used during star measurement + star_background = pexConfig.ConfigurableField( + target=lsst.meas.algorithms.SubtractBackgroundTask, + 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." @@ -462,8 +477,7 @@ def setDefaults(self): # 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.reEstimateBackground = False # Minimal measurement plugins for PSF determination. # TODO DM-39203: We can drop GaussianFlux and PsfFlux, if we use @@ -494,6 +508,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", @@ -536,8 +551,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 @@ -628,6 +643,26 @@ def validate(self): raise pexConfig.FieldValidationError(CalibrateImageConfig.run_sattle, self, "Sattle requested but URI environment variable not set.") + if not self.do_adaptive_threshold_detection: + if not self.psf_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.psf_detection " + "background must be configured with " + "reEstimateBackground = True to maintain original behavior." + ) + 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): """Compute the PSF, aperture corrections, astrometric and photometric @@ -659,6 +694,7 @@ def __init__(self, initial_stars_schema=None, **kwargs): afwTable.CoordKey.addErrorFields(self.psf_schema) self.makeSubtask("psf_detection", schema=self.psf_schema) self.makeSubtask("psf_source_measurement", schema=self.psf_schema) + self.makeSubtask("psf_adaptive_threshold_detection") self.makeSubtask("psf_measure_psf", schema=self.psf_schema) self.makeSubtask("psf_normalized_calibration_flux", schema=self.psf_schema) @@ -689,6 +725,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) @@ -906,7 +943,7 @@ def run( illumination_correction, ) - result.psf_stars_footprints, result.background, _ = self._compute_psf( + result.psf_stars_footprints, result.background, _, adaptiveDetResStruct = self._compute_psf( result.exposure, id_generator, background_to_photometric_ratio=result.background_to_photometric_ratio, @@ -926,32 +963,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, + adaptiveDetResStruct=adaptiveDetResStruct, + 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 @@ -1108,7 +1163,10 @@ def log_psf(msg, addToMetadata=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 +1184,65 @@ 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, - ) + if not self.config.do_adaptive_threshold_detection: + adaptiveDetResStruct = None + # 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, + ) + else: + initialThreshold = self.config.psf_detection.thresholdValue + initialThresholdMultiplier = self.config.psf_detection.includeThresholdMultiplier + adaptiveDetResStruct = self.psf_adaptive_threshold_detection.run( + table, exposure, + initialThreshold=initialThreshold, + initialThresholdMultiplier=initialThresholdMultiplier, + doReEstimageBackgroud=False, + ) + detections = adaptiveDetResStruct.detections + 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 +1252,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.sources, background, psf_result.cellSet, adaptiveDetResStruct def _measure_aperture_correction(self, exposure, bright_sources): """Measure and set the ApCorrMap on the Exposure, using @@ -1209,7 +1284,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, + adaptiveDetResStruct=None, num_astrometry_matches=None): """Detect stars on an exposure that has a PSF model, and measure their PSF, circular aperture, compensated gaussian fluxes. @@ -1234,14 +1310,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 adaptiveDetResStruct is not None: + while inAdaptiveDet and adaptiveDetIter < maxAdaptiveDetIter: + inAdaptiveDet = False + adaptiveDetectionConfig = lsst.meas.algorithms.SourceDetectionConfig() + adaptiveDetectionConfig.reEstimateBackground = False + adaptiveDetectionConfig.includeThresholdMultiplier = 2.0 + psfThreshold = ( + adaptiveDetResStruct.thresholdValue*adaptiveDetResStruct.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) @@ -1300,7 +1430,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 +1450,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] @@ -1501,3 +1632,253 @@ def _update_wcs_with_camera_model(self, exposure, cameraModel): detector = cameraModel[exposure.detector.getId()] newWcs = createInitialSkyWcs(exposure.visitInfo, detector) exposure.setWcs(newWcs) + + def _compute_mask_fraction(self, mask, mask_planes, bad_mask_planes): + """Evaluate the fraction of masked pixels in a (set of) mask plane(s). + + Parameters + ---------- + mask : `lsst.afw.image.Mask` + The mask on which to evaluate the fraction. + mask_planes : `list`, `str` + The mask planes for which to evaluate the fraction. + bad_mask_planes : `list`, `str` + The mask planes to exclude from the computation. + + Returns + ------- + detected_fraction : `float` + The calculated fraction of masked pixels + """ + bad_pixel_mask = afwImage.Mask.getPlaneBitMask(bad_mask_planes) + n_good_pix = np.sum(mask.array & bad_pixel_mask == 0) + if n_good_pix == 0: + detected_fraction = float("nan") + return detected_fraction + detected_pixel_mask = afwImage.Mask.getPlaneBitMask(mask_planes) + n_detected_pix = np.sum((mask.array & detected_pixel_mask != 0) + & (mask.array & bad_pixel_mask == 0)) + detected_fraction = n_detected_pix/n_good_pix + return detected_fraction + + def _compute_per_amp_fraction(self, exposure, detected_fraction, mask_planes, bad_mask_planes): + """Evaluate the maximum per-amplifier fraction of masked pixels. + + Parameters + ---------- + exposure : `lsst.afw.image.ExposureF` + The exposure on which to compute the per-amp masked fraction. + detected_fraction : `float` + The current detected_fraction of the ``mask_planes`` for the + full detector. + mask_planes : `list`, `str` + The mask planes for which to evaluate the fraction. + bad_mask_planes : `list`, `str` + The mask planes to exclude from the computation. + + Returns + ------- + n_above_max_per_amp : `int` + The number of amplifiers with masked fractions above a maximum + value (set by the current full-detector ``detected_fraction``). + highest_detected_fraction_per_amp : `float` + The highest value of the per-amplifier fraction of masked pixels. + no_zero_det_amps : `bool` + A boolean representing whether any of the amplifiers has zero + masked pixels. + """ + highest_detected_fraction_per_amp = -9.99 + n_above_max_per_amp = 0 + n_no_zero_det_amps = 0 + no_zero_det_amps = True + amps = exposure.detector.getAmplifiers() + if amps is not None: + for ia, amp in enumerate(amps): + amp_bbox = amp.getBBox() + exp_bbox = exposure.getBBox() + if not exp_bbox.contains(amp_bbox): + self.log.info("Bounding box of amplifier (%s) does not fit in exposure's " + "bounding box (%s). Skipping...", amp_bbox, exp_bbox) + continue + sub_image = exposure.subset(amp.getBBox()) + detected_fraction_amp = self._compute_mask_fraction(sub_image.mask, + mask_planes, + bad_mask_planes) + self.log.debug("Current detected fraction for amplifier %s = %.3f", + amp.getName(), detected_fraction_amp) + if detected_fraction_amp < 0.002: + n_no_zero_det_amps += 1 + if n_no_zero_det_amps > 2: + no_zero_det_amps = False + break + highest_detected_fraction_per_amp = max(detected_fraction_amp, + highest_detected_fraction_per_amp) + if highest_detected_fraction_per_amp > min(0.998, max(0.8, 3.0*detected_fraction)): + n_above_max_per_amp += 1 + if n_above_max_per_amp > 2: + break + else: + self.log.info("No amplifier object for detector %d, so skipping per-amp " + "detection fraction checks.", exposure.detector.getId()) + return n_above_max_per_amp, highest_detected_fraction_per_amp, no_zero_det_amps + + 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. + """ + # Restore the previously measured backgroud and remeasure it + # using an adaptive threshold detection iteration to ensure a + # "Goldilocks Zone" for the fraction of detected pixels. + backgroundOrig = result.background.clone() + median_background = np.median(backgroundOrig.getImage().array) + self.log.warning("Original median_background = %.2f", median_background) + result.exposure.image.array += result.background.getImage().array + result.background = afwMath.BackgroundList() + + origMask = result.exposure.mask.clone() + bad_mask_planes = ["BAD", "EDGE", "NO_DATA"] + detected_mask_planes = ["DETECTED", "DETECTED_NEGATIVE"] + nPixToDilate = 10 + detected_fraction_orig = self._compute_mask_fraction(result.exposure.mask, + detected_mask_planes, + bad_mask_planes) + minDetFracForFinalBg = 0.02 + maxDetFracForFinalBg = 0.93 + doDilatedOrigDetectionMask = True + if doDilatedOrigDetectionMask: + # Dilate the current detected mask planes and don't clear + # them in the detection step. + inDilating = True + while inDilating: + dilatedMask = origMask.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) + if detected_fraction_dilated < maxDetFracForFinalBg or nPixToDilate == 1: + inDilating = False + else: + nPixToDilate -= 1 + result.exposure.mask = dilatedMask + self.log.warning("detected_fraction_orig = %.3f detected_fraction_dilated = %.3f", + detected_fraction_orig, detected_fraction_dilated) + n_above_max_per_amp = -99 + highest_detected_fraction_per_amp = float("nan") + doCheckPerAmpDetFraction = True + if doCheckPerAmpDetFraction: # detected_fraction < maxDetFracForFinalBg: + n_above_max_per_amp, highest_detected_fraction_per_amp, no_zero_det_amps = \ + self._compute_per_amp_fraction(result.exposure, detected_fraction_dilated, + detected_mask_planes, bad_mask_planes) + self.log.warning("Dilated mask: n_above_max_per_amp = %d, " + "highest_detected_fraction_per_amp = %.3f", + n_above_max_per_amp, highest_detected_fraction_per_amp) + + detected_fraction = 1.0 + inBackgroundDet = True + detected_fraction = 1.0 if inBackgroundDet else detected_fraction_dilated + maxIter = 40 + nIter = 0 + nFootprintTemp = 1e12 + starBackgroundDetectionConfig = lsst.meas.algorithms.SourceDetectionConfig() + starBackgroundDetectionConfig.doTempLocalBackground = False + starBackgroundDetectionConfig.nSigmaToGrow = 70.0 + starBackgroundDetectionConfig.reEstimateBackground = False + starBackgroundDetectionConfig.includeThresholdMultiplier = 1.0 + starBackgroundDetectionConfig.thresholdValue = max(2.0, 0.2*median_background) + starBackgroundDetectionConfig.thresholdType = "pixel_stdev" # "stdev" + + n_above_max_per_amp = -99 + highest_detected_fraction_per_amp = float("nan") + doCheckPerAmpDetFraction = True + + while inBackgroundDet: + currentThresh = starBackgroundDetectionConfig.thresholdValue + if detected_fraction > maxDetFracForFinalBg: + starBackgroundDetectionConfig.thresholdValue = 1.07*currentThresh + if nFootprintTemp < 3 and detected_fraction > 0.9*maxDetFracForFinalBg: + starBackgroundDetectionConfig.thresholdValue = 1.2*currentThresh + if n_above_max_per_amp > 1: + starBackgroundDetectionConfig.thresholdValue = 1.1*currentThresh + if detected_fraction < minDetFracForFinalBg: + starBackgroundDetectionConfig.thresholdValue = 0.8*currentThresh + starBackgroundDetectionTask = lsst.meas.algorithms.SourceDetectionTask( + config=starBackgroundDetectionConfig) + table = afwTable.SourceTable.make(self.initial_stars_schema.schema) + tempDetections = starBackgroundDetectionTask.run( + table=table, exposure=result.exposure, clearMask=True) + result.exposure.mask |= dilatedMask + nFootprintTemp = len(tempDetections.sources) + detected_fraction = self._compute_mask_fraction(result.exposure.mask, detected_mask_planes, + bad_mask_planes) + self.log.info("nIter = %d, thresh = %.2f: Fraction of pixels marked as DETECTED or " + "DETECTED_NEGATIVE in star_background_detection = %.3f " + "(max is %.3f; min is %.3f)", + nIter, starBackgroundDetectionConfig.thresholdValue, + detected_fraction, maxDetFracForFinalBg, minDetFracForFinalBg) + + n_amp = len(result.exposure.detector.getAmplifiers()) + if doCheckPerAmpDetFraction: # detected_fraction < maxDetFracForFinalBg: + n_above_max_per_amp, highest_detected_fraction_per_amp, no_zero_det_amps = \ + self._compute_per_amp_fraction(result.exposure, detected_fraction, + detected_mask_planes, bad_mask_planes) + + if not no_zero_det_amps: + starBackgroundDetectionConfig.thresholdValue = 0.95*currentThresh + nIter += 1 + if nIter > maxIter: + inBackgroundDet = False + + if (detected_fraction < maxDetFracForFinalBg and detected_fraction > minDetFracForFinalBg + and n_above_max_per_amp < int(0.75*n_amp) + and no_zero_det_amps): + if (n_above_max_per_amp < max(1, int(0.15*n_amp)) + or detected_fraction < 0.85*maxDetFracForFinalBg): + inBackgroundDet = False + else: + self.log.warning("Making small tweak....") + starBackgroundDetectionConfig.thresholdValue = 1.05*currentThresh + self.log.warning("n_above_max_per_amp = %d (abs max is %d)", n_above_max_per_amp, int(0.75*n_amp)) + + self.log.info("Fraction of pixels marked as DETECTED or DETECTED_NEGATIVE is now %.5f " + "(highest per amp section = %.5f)", + detected_fraction, highest_detected_fraction_per_amp) + + if detected_fraction > maxDetFracForFinalBg: + result.exposure.mask = dilatedMask + self.log.warning("Final fraction of pixels marked as DETECTED or DETECTED_NEGATIVE " + "was too large in star_background_detection = %.3f (max = %.3f). " + "Reverting to dilated mask from PSF detection...", + detected_fraction, maxDetFracForFinalBg) + star_background = self.star_background.run(exposure=result.exposure).background + result.background.append(star_background[0]) + + # Clear detected mask plane before final round of detection + mask = result.exposure.mask + for mp in detected_mask_planes: + if mp not in mask.getMaskPlaneDict(): + mask.addMaskPlane(mp) + mask &= ~mask.getPlaneBitMask(detected_mask_planes) + + return result From 50a61a751c24fa0aab2aa0574b6e26a62cd1537e Mon Sep 17 00:00:00 2001 From: Lauren MacArthur Date: Mon, 25 Aug 2025 16:39:57 -0700 Subject: [PATCH 03/12] Update test_calibrateImage for adaptive detection --- tests/test_calibrateImage.py | 48 ++++++++++++++++++++++++++++-------- 1 file changed, 38 insertions(+), 10 deletions(-) diff --git a/tests/test_calibrateImage.py b/tests/test_calibrateImage.py index e656f7276..8fa33841e 100644 --- a/tests/test_calibrateImage.py +++ b/tests/test_calibrateImage.py @@ -123,9 +123,14 @@ 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 + # Maintain original, no adaptive threshold detection, configs values. + self.config.do_adaptive_threshold_detection = False + self.config.psf_detection.reEstimateBackground = True + self.config.star_detection.reEstimateBackground = True # ApFlux has more noise than PsfFlux (the latter unrealistically small # in this test data), so we need to do magnitude rejection at higher # sigma, otherwise we can lose otherwise good sources. @@ -154,7 +159,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 +174,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 +223,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 +241,28 @@ def test_run(self): self._check_run(calibrate, result) + def test_run_adaptive_threshold_deteection(self): + """Test that run() runs with adaptive threshold detection turned on. + """ + config = copy.copy(self.config) + # Set the adaptive threshold detection, config values... + config.do_adaptive_threshold_detection = True + config.psf_adaptive_threshold_detection.minFootprint = 4 + config.psf_adaptive_threshold_detection.minIsolated = 4 + config.psf_adaptive_threshold_detection.sufficientIsolated = 4 + config.psf_detection.reEstimateBackground = False + config.star_detection.reEstimateBackground = False + + calibrate = CalibrateImageTask(config=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=1) + def test_run_downsample(self): """Test that run() runs with downsample. """ @@ -313,7 +341,7 @@ 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_stars, background, candidates, _ = calibrate._compute_psf(self.exposure, self.id_generator) # Catalog ids should be very large from this id generator. self.assertTrue(all(psf_stars['id'] > 1000000000)) @@ -385,7 +413,7 @@ 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_stars, background, candidates, _ = calibrate._compute_psf(self.exposure, self.id_generator) # First check that the exposure doesn't have an ApCorrMap. self.assertIsNone(self.exposure.apCorrMap) @@ -400,7 +428,7 @@ 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_stars, background, candidates, _ = calibrate._compute_psf(self.exposure, self.id_generator) calibrate._measure_aperture_correction(self.exposure, psf_stars) stars = calibrate._find_stars(self.exposure, background, self.id_generator) @@ -428,7 +456,7 @@ 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_stars, background, candidates, _ = calibrate._compute_psf(self.exposure, self.id_generator) calibrate._measure_aperture_correction(self.exposure, psf_stars) stars = calibrate._find_stars(self.exposure, background, self.id_generator) @@ -448,7 +476,7 @@ 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_stars, background, candidates, _ = calibrate._compute_psf(self.exposure, self.id_generator) calibrate._measure_aperture_correction(self.exposure, psf_stars) stars = calibrate._find_stars(self.exposure, background, self.id_generator) calibrate._fit_astrometry(self.exposure, stars) @@ -489,7 +517,7 @@ 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_stars, background, candidates, _ = calibrate._compute_psf(self.exposure, self.id_generator) calibrate._measure_aperture_correction(self.exposure, psf_stars) stars = calibrate._find_stars(self.exposure, background, self.id_generator) From c8f5abaae4e83fc03d28750f12953d0d7ba4ba1c Mon Sep 17 00:00:00 2001 From: Lauren MacArthur Date: Wed, 10 Sep 2025 14:43:48 -0700 Subject: [PATCH 04/12] Flake8 fixes. --- python/lsst/pipe/tasks/calibrateImage.py | 62 +++++++++++++----------- 1 file changed, 35 insertions(+), 27 deletions(-) diff --git a/python/lsst/pipe/tasks/calibrateImage.py b/python/lsst/pipe/tasks/calibrateImage.py index e9d24f5db..a3db15b93 100644 --- a/python/lsst/pipe/tasks/calibrateImage.py +++ b/python/lsst/pipe/tasks/calibrateImage.py @@ -311,8 +311,8 @@ 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." @@ -466,8 +466,8 @@ 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. @@ -482,8 +482,8 @@ def setDefaults(self): # 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", @@ -527,9 +527,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. @@ -709,7 +709,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: @@ -843,7 +844,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 @@ -891,10 +893,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`) @@ -1024,7 +1028,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 @@ -1044,12 +1049,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: @@ -1087,7 +1092,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 ------- @@ -1157,8 +1163,8 @@ def log_psf(msg, addToMetadata=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). + 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() @@ -1268,8 +1274,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, @@ -1402,7 +1408,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. @@ -1576,7 +1583,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() From fb78891c9fe6413f6ebb4cd1590e19bb609cf427 Mon Sep 17 00:00:00 2001 From: Lauren MacArthur Date: Thu, 23 Oct 2025 16:46:18 -0700 Subject: [PATCH 05/12] Slight hack for unpersistable empty BackgroundList --- python/lsst/pipe/tasks/multiBand.py | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) 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): From 2aa0488baf5a00b8613d5723157754e393d8d935 Mon Sep 17 00:00:00 2001 From: Lauren MacArthur Date: Mon, 3 Nov 2025 14:11:18 -0800 Subject: [PATCH 06/12] Add final pedetstal backgrond in calibrateImage --- python/lsst/pipe/tasks/calibrateImage.py | 40 ++++++++++++++++++++++++ tests/test_calibrateImage.py | 2 +- 2 files changed, 41 insertions(+), 1 deletion(-) diff --git a/python/lsst/pipe/tasks/calibrateImage.py b/python/lsst/pipe/tasks/calibrateImage.py index a3db15b93..ed8ff7589 100644 --- a/python/lsst/pipe/tasks/calibrateImage.py +++ b/python/lsst/pipe/tasks/calibrateImage.py @@ -1882,6 +1882,46 @@ def _remeasure_star_background(self, result): star_background = self.star_background.run(exposure=result.exposure).background result.background.append(star_background[0]) + # 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) + # Clear detected mask plane before final round of detection mask = result.exposure.mask for mp in detected_mask_planes: diff --git a/tests/test_calibrateImage.py b/tests/test_calibrateImage.py index 8fa33841e..e8a92ebea 100644 --- a/tests/test_calibrateImage.py +++ b/tests/test_calibrateImage.py @@ -261,7 +261,7 @@ def test_run_adaptive_threshold_deteection(self): subString = "Using adaptive threshold detection " self.assertTrue(any(subString in s for s in cm.output)) - self._check_run(calibrate, result, expect_n_background=1) + self._check_run(calibrate, result, expect_n_background=2) def test_run_downsample(self): """Test that run() runs with downsample. From 2a7188f21f592fb0cf27fec76d389e3c334b0410 Mon Sep 17 00:00:00 2001 From: Jim Bosch Date: Fri, 31 Oct 2025 10:49:29 -0400 Subject: [PATCH 07/12] Move new adaptive threshold detection task to its own module. --- python/lsst/pipe/tasks/calibrateImage.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/python/lsst/pipe/tasks/calibrateImage.py b/python/lsst/pipe/tasks/calibrateImage.py index ed8ff7589..80bc40424 100644 --- a/python/lsst/pipe/tasks/calibrateImage.py +++ b/python/lsst/pipe/tasks/calibrateImage.py @@ -35,6 +35,7 @@ import lsst.meas.algorithms.installGaussianPsf import lsst.meas.algorithms.measureApCorr import lsst.meas.algorithms.setPrimaryFlags +from lsst.meas.algorithms.adaptive_thresholds import AdaptiveThresholdDetectionTask import lsst.meas.base import lsst.meas.astrom import lsst.meas.deblender @@ -294,7 +295,7 @@ class CalibrateImageConfig(pipeBase.PipelineTaskConfig, pipelineConnections=Cali doc="Implement the adaptive detection thresholding approach?", ) psf_adaptive_threshold_detection = pexConfig.ConfigurableField( - target=lsst.meas.algorithms.AdaptiveThresholdDetectionTask, + target=AdaptiveThresholdDetectionTask, doc="Task to adaptively detect sources for PSF determination.", ) psf_source_measurement = pexConfig.ConfigurableField( From 7d7a8fed36f13719a890be1ccc5f114ea867e574 Mon Sep 17 00:00:00 2001 From: Jim Bosch Date: Fri, 31 Oct 2025 11:08:43 -0400 Subject: [PATCH 08/12] Drop support for bg reestimation in adaptive det task. --- python/lsst/pipe/tasks/calibrateImage.py | 1 - 1 file changed, 1 deletion(-) diff --git a/python/lsst/pipe/tasks/calibrateImage.py b/python/lsst/pipe/tasks/calibrateImage.py index 80bc40424..4956b13f4 100644 --- a/python/lsst/pipe/tasks/calibrateImage.py +++ b/python/lsst/pipe/tasks/calibrateImage.py @@ -1208,7 +1208,6 @@ def log_psf(msg, addToMetadata=False): table, exposure, initialThreshold=initialThreshold, initialThresholdMultiplier=initialThresholdMultiplier, - doReEstimageBackgroud=False, ) detections = adaptiveDetResStruct.detections From 1146cdd9144f9587ff1748bde5eeb5539dfd337e Mon Sep 17 00:00:00 2001 From: Jim Bosch Date: Sat, 1 Nov 2025 08:45:02 -0400 Subject: [PATCH 09/12] Fix doc indent. --- python/lsst/pipe/tasks/calibrateImage.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/python/lsst/pipe/tasks/calibrateImage.py b/python/lsst/pipe/tasks/calibrateImage.py index 4956b13f4..8cf1f4ef4 100644 --- a/python/lsst/pipe/tasks/calibrateImage.py +++ b/python/lsst/pipe/tasks/calibrateImage.py @@ -1161,11 +1161,11 @@ 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() From 07d5d00f7188e957c544fd58627bcbf511fef9fb Mon Sep 17 00:00:00 2001 From: Jim Bosch Date: Sun, 2 Nov 2025 08:37:07 -0500 Subject: [PATCH 10/12] Make CalibrateImageTask astrometry test match production. The unit test was still fitting astrometry to the second-stage star catalog, even though the task itself had switched to using the PSF star catalog. --- tests/test_calibrateImage.py | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/tests/test_calibrateImage.py b/tests/test_calibrateImage.py index e8a92ebea..188c4d304 100644 --- a/tests/test_calibrateImage.py +++ b/tests/test_calibrateImage.py @@ -458,13 +458,11 @@ def test_astrometry(self): calibrate.astrometry.setRefObjLoader(self.ref_loader) psf_stars, background, candidates, _ = calibrate._compute_psf(self.exposure, self.id_generator) 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") + sky = psf_stars["sky_source"] + fitted = SkyCoord(psf_stars[~sky]['coord_ra'], psf_stars[~sky]['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) From bc7a85a5651f573eb578cfe3046b77cfc62dbbd1 Mon Sep 17 00:00:00 2001 From: Jim Bosch Date: Sun, 2 Nov 2025 09:56:31 -0500 Subject: [PATCH 11/12] Rework adaptive detection into a replacement for regular detection. --- python/lsst/pipe/tasks/calibrateImage.py | 88 ++++++++++-------------- tests/test_calibrateImage.py | 63 ++++++++++------- 2 files changed, 74 insertions(+), 77 deletions(-) diff --git a/python/lsst/pipe/tasks/calibrateImage.py b/python/lsst/pipe/tasks/calibrateImage.py index 8cf1f4ef4..67a92fe6b 100644 --- a/python/lsst/pipe/tasks/calibrateImage.py +++ b/python/lsst/pipe/tasks/calibrateImage.py @@ -286,7 +286,7 @@ 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( @@ -294,10 +294,6 @@ class CalibrateImageConfig(pipeBase.PipelineTaskConfig, pipelineConnections=Cali default=True, doc="Implement the adaptive detection thresholding approach?", ) - psf_adaptive_threshold_detection = pexConfig.ConfigurableField( - target=AdaptiveThresholdDetectionTask, - doc="Task to adaptively detect sources for PSF determination.", - ) psf_source_measurement = pexConfig.ConfigurableField( target=lsst.meas.base.SingleFrameMeasurementTask, doc="Task to measure sources to be used for psf estimation." @@ -472,13 +468,8 @@ def setDefaults(self): # 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 - self.psf_detection.reEstimateBackground = False + 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 @@ -622,7 +613,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, @@ -630,7 +621,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, @@ -643,17 +634,17 @@ 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.psf_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.psf_detection " - "background must be configured with " - "reEstimateBackground = True to maintain original behavior." - ) if not self.star_detection.reEstimateBackground: raise pexConfig.FieldValidationError( CalibrateImageConfig.psf_detection, @@ -695,7 +686,6 @@ def __init__(self, initial_stars_schema=None, **kwargs): afwTable.CoordKey.addErrorFields(self.psf_schema) self.makeSubtask("psf_detection", schema=self.psf_schema) self.makeSubtask("psf_source_measurement", schema=self.psf_schema) - self.makeSubtask("psf_adaptive_threshold_detection") self.makeSubtask("psf_measure_psf", schema=self.psf_schema) self.makeSubtask("psf_normalized_calibration_flux", schema=self.psf_schema) @@ -948,11 +938,12 @@ def run( illumination_correction, ) - result.psf_stars_footprints, result.background, _, adaptiveDetResStruct = 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 @@ -999,7 +990,7 @@ def run( result.background, id_generator, background_to_photometric_ratio=result.background_to_photometric_ratio, - adaptiveDetResStruct=adaptiveDetResStruct, + psf_detections=psf_detections, num_astrometry_matches=num_astrometry_matches, ) psf = result.exposure.getPsf() @@ -1147,8 +1138,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` @@ -1191,25 +1184,12 @@ 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()) - if not self.config.do_adaptive_threshold_detection: - adaptiveDetResStruct = None - # 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, - ) - else: - initialThreshold = self.config.psf_detection.thresholdValue - initialThresholdMultiplier = self.config.psf_detection.includeThresholdMultiplier - adaptiveDetResStruct = self.psf_adaptive_threshold_detection.run( - table, exposure, - initialThreshold=initialThreshold, - initialThresholdMultiplier=initialThresholdMultiplier, - ) - detections = adaptiveDetResStruct.detections + 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 @@ -1258,7 +1238,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, adaptiveDetResStruct + return detections, background, psf_result.cellSet def _measure_aperture_correction(self, exposure, bright_sources): """Measure and set the ApCorrMap on the Exposure, using @@ -1290,8 +1270,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, - adaptiveDetResStruct=None, num_astrometry_matches=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. @@ -1307,6 +1287,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 ------- @@ -1321,14 +1305,14 @@ def _find_stars(self, exposure, background, id_generator, background_to_photomet maxAdaptiveDetIter = 8 adaptiveDetIter = 0 threshFactor = 0.2 - if adaptiveDetResStruct is not None: + 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 = ( - adaptiveDetResStruct.thresholdValue*adaptiveDetResStruct.includeThresholdMultiplier + psf_detections.thresholdValue*psf_detections.includeThresholdMultiplier ) adaptiveDetectionConfig.thresholdValue = max( self.config.star_detection.thresholdValue, diff --git a/tests/test_calibrateImage.py b/tests/test_calibrateImage.py index 188c4d304..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 @@ -127,10 +136,6 @@ def setUp(self): self.config.astrometry.sourceSelector["science"].flags.good = [] self.config.astrometry.matcher.numPointsForShape = 3 self.config.run_sattle = False - # Maintain original, no adaptive threshold detection, configs values. - self.config.do_adaptive_threshold_detection = False - self.config.psf_detection.reEstimateBackground = True - self.config.star_detection.reEstimateBackground = True # ApFlux has more noise than PsfFlux (the latter unrealistically small # in this test data), so we need to do magnitude rejection at higher # sigma, otherwise we can lose otherwise good sources. @@ -241,19 +246,20 @@ def test_run(self): self._check_run(calibrate, result) - def test_run_adaptive_threshold_deteection(self): + def test_run_adaptive_threshold_detection(self): """Test that run() runs with adaptive threshold detection turned on. """ - config = copy.copy(self.config) # Set the adaptive threshold detection, config values... - config.do_adaptive_threshold_detection = True - config.psf_adaptive_threshold_detection.minFootprint = 4 - config.psf_adaptive_threshold_detection.minIsolated = 4 - config.psf_adaptive_threshold_detection.sufficientIsolated = 4 - config.psf_detection.reEstimateBackground = False - config.star_detection.reEstimateBackground = False + 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=config) + 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: @@ -341,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)) @@ -413,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) @@ -428,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)) @@ -456,13 +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) calibrate._fit_astrometry(self.exposure, psf_stars) # Check that we got reliable matches with the truth coordinates. - sky = psf_stars["sky_source"] - fitted = SkyCoord(psf_stars[~sky]['coord_ra'], psf_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) @@ -474,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) @@ -515,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) @@ -829,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) From ff6d1b154d6528c27a36bba598c53bf62e94ac45 Mon Sep 17 00:00:00 2001 From: Jim Bosch Date: Mon, 3 Nov 2025 09:50:05 -0500 Subject: [PATCH 12/12] Move adaptive background step to a new subtask. --- python/lsst/pipe/tasks/calibrateImage.py | 240 +---------------------- 1 file changed, 6 insertions(+), 234 deletions(-) diff --git a/python/lsst/pipe/tasks/calibrateImage.py b/python/lsst/pipe/tasks/calibrateImage.py index 67a92fe6b..2fd99f04c 100644 --- a/python/lsst/pipe/tasks/calibrateImage.py +++ b/python/lsst/pipe/tasks/calibrateImage.py @@ -35,7 +35,10 @@ import lsst.meas.algorithms.installGaussianPsf import lsst.meas.algorithms.measureApCorr import lsst.meas.algorithms.setPrimaryFlags -from lsst.meas.algorithms.adaptive_thresholds import AdaptiveThresholdDetectionTask +from lsst.meas.algorithms.adaptive_thresholds import ( + AdaptiveThresholdDetectionTask, + AdaptiveThresholdBackgroundTask, +) import lsst.meas.base import lsst.meas.astrom import lsst.meas.deblender @@ -317,7 +320,7 @@ class CalibrateImageConfig(pipeBase.PipelineTaskConfig, pipelineConnections=Cali # subtasks used during star measurement star_background = pexConfig.ConfigurableField( - target=lsst.meas.algorithms.SubtractBackgroundTask, + target=AdaptiveThresholdBackgroundTask, doc="Task to perform final background subtraction, just before photoCal.", ) star_detection = pexConfig.ConfigurableField( @@ -1625,95 +1628,6 @@ def _update_wcs_with_camera_model(self, exposure, cameraModel): newWcs = createInitialSkyWcs(exposure.visitInfo, detector) exposure.setWcs(newWcs) - def _compute_mask_fraction(self, mask, mask_planes, bad_mask_planes): - """Evaluate the fraction of masked pixels in a (set of) mask plane(s). - - Parameters - ---------- - mask : `lsst.afw.image.Mask` - The mask on which to evaluate the fraction. - mask_planes : `list`, `str` - The mask planes for which to evaluate the fraction. - bad_mask_planes : `list`, `str` - The mask planes to exclude from the computation. - - Returns - ------- - detected_fraction : `float` - The calculated fraction of masked pixels - """ - bad_pixel_mask = afwImage.Mask.getPlaneBitMask(bad_mask_planes) - n_good_pix = np.sum(mask.array & bad_pixel_mask == 0) - if n_good_pix == 0: - detected_fraction = float("nan") - return detected_fraction - detected_pixel_mask = afwImage.Mask.getPlaneBitMask(mask_planes) - n_detected_pix = np.sum((mask.array & detected_pixel_mask != 0) - & (mask.array & bad_pixel_mask == 0)) - detected_fraction = n_detected_pix/n_good_pix - return detected_fraction - - def _compute_per_amp_fraction(self, exposure, detected_fraction, mask_planes, bad_mask_planes): - """Evaluate the maximum per-amplifier fraction of masked pixels. - - Parameters - ---------- - exposure : `lsst.afw.image.ExposureF` - The exposure on which to compute the per-amp masked fraction. - detected_fraction : `float` - The current detected_fraction of the ``mask_planes`` for the - full detector. - mask_planes : `list`, `str` - The mask planes for which to evaluate the fraction. - bad_mask_planes : `list`, `str` - The mask planes to exclude from the computation. - - Returns - ------- - n_above_max_per_amp : `int` - The number of amplifiers with masked fractions above a maximum - value (set by the current full-detector ``detected_fraction``). - highest_detected_fraction_per_amp : `float` - The highest value of the per-amplifier fraction of masked pixels. - no_zero_det_amps : `bool` - A boolean representing whether any of the amplifiers has zero - masked pixels. - """ - highest_detected_fraction_per_amp = -9.99 - n_above_max_per_amp = 0 - n_no_zero_det_amps = 0 - no_zero_det_amps = True - amps = exposure.detector.getAmplifiers() - if amps is not None: - for ia, amp in enumerate(amps): - amp_bbox = amp.getBBox() - exp_bbox = exposure.getBBox() - if not exp_bbox.contains(amp_bbox): - self.log.info("Bounding box of amplifier (%s) does not fit in exposure's " - "bounding box (%s). Skipping...", amp_bbox, exp_bbox) - continue - sub_image = exposure.subset(amp.getBBox()) - detected_fraction_amp = self._compute_mask_fraction(sub_image.mask, - mask_planes, - bad_mask_planes) - self.log.debug("Current detected fraction for amplifier %s = %.3f", - amp.getName(), detected_fraction_amp) - if detected_fraction_amp < 0.002: - n_no_zero_det_amps += 1 - if n_no_zero_det_amps > 2: - no_zero_det_amps = False - break - highest_detected_fraction_per_amp = max(detected_fraction_amp, - highest_detected_fraction_per_amp) - if highest_detected_fraction_per_amp > min(0.998, max(0.8, 3.0*detected_fraction)): - n_above_max_per_amp += 1 - if n_above_max_per_amp > 2: - break - else: - self.log.info("No amplifier object for detector %d, so skipping per-amp " - "detection fraction checks.", exposure.detector.getId()) - return n_above_max_per_amp, highest_detected_fraction_per_amp, no_zero_det_amps - def _remeasure_star_background(self, result): """Remeasure the exposure's background with iterative adaptive threshold detection. @@ -1729,142 +1643,7 @@ def _remeasure_star_background(self, result): result : `lsst.pipe.base.Struct` The modified result Struct with the new background subtracted. """ - # Restore the previously measured backgroud and remeasure it - # using an adaptive threshold detection iteration to ensure a - # "Goldilocks Zone" for the fraction of detected pixels. - backgroundOrig = result.background.clone() - median_background = np.median(backgroundOrig.getImage().array) - self.log.warning("Original median_background = %.2f", median_background) - result.exposure.image.array += result.background.getImage().array - result.background = afwMath.BackgroundList() - - origMask = result.exposure.mask.clone() - bad_mask_planes = ["BAD", "EDGE", "NO_DATA"] - detected_mask_planes = ["DETECTED", "DETECTED_NEGATIVE"] - nPixToDilate = 10 - detected_fraction_orig = self._compute_mask_fraction(result.exposure.mask, - detected_mask_planes, - bad_mask_planes) - minDetFracForFinalBg = 0.02 - maxDetFracForFinalBg = 0.93 - doDilatedOrigDetectionMask = True - if doDilatedOrigDetectionMask: - # Dilate the current detected mask planes and don't clear - # them in the detection step. - inDilating = True - while inDilating: - dilatedMask = origMask.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) - if detected_fraction_dilated < maxDetFracForFinalBg or nPixToDilate == 1: - inDilating = False - else: - nPixToDilate -= 1 - result.exposure.mask = dilatedMask - self.log.warning("detected_fraction_orig = %.3f detected_fraction_dilated = %.3f", - detected_fraction_orig, detected_fraction_dilated) - n_above_max_per_amp = -99 - highest_detected_fraction_per_amp = float("nan") - doCheckPerAmpDetFraction = True - if doCheckPerAmpDetFraction: # detected_fraction < maxDetFracForFinalBg: - n_above_max_per_amp, highest_detected_fraction_per_amp, no_zero_det_amps = \ - self._compute_per_amp_fraction(result.exposure, detected_fraction_dilated, - detected_mask_planes, bad_mask_planes) - self.log.warning("Dilated mask: n_above_max_per_amp = %d, " - "highest_detected_fraction_per_amp = %.3f", - n_above_max_per_amp, highest_detected_fraction_per_amp) - - detected_fraction = 1.0 - inBackgroundDet = True - detected_fraction = 1.0 if inBackgroundDet else detected_fraction_dilated - maxIter = 40 - nIter = 0 - nFootprintTemp = 1e12 - starBackgroundDetectionConfig = lsst.meas.algorithms.SourceDetectionConfig() - starBackgroundDetectionConfig.doTempLocalBackground = False - starBackgroundDetectionConfig.nSigmaToGrow = 70.0 - starBackgroundDetectionConfig.reEstimateBackground = False - starBackgroundDetectionConfig.includeThresholdMultiplier = 1.0 - starBackgroundDetectionConfig.thresholdValue = max(2.0, 0.2*median_background) - starBackgroundDetectionConfig.thresholdType = "pixel_stdev" # "stdev" - - n_above_max_per_amp = -99 - highest_detected_fraction_per_amp = float("nan") - doCheckPerAmpDetFraction = True - - while inBackgroundDet: - currentThresh = starBackgroundDetectionConfig.thresholdValue - if detected_fraction > maxDetFracForFinalBg: - starBackgroundDetectionConfig.thresholdValue = 1.07*currentThresh - if nFootprintTemp < 3 and detected_fraction > 0.9*maxDetFracForFinalBg: - starBackgroundDetectionConfig.thresholdValue = 1.2*currentThresh - if n_above_max_per_amp > 1: - starBackgroundDetectionConfig.thresholdValue = 1.1*currentThresh - if detected_fraction < minDetFracForFinalBg: - starBackgroundDetectionConfig.thresholdValue = 0.8*currentThresh - starBackgroundDetectionTask = lsst.meas.algorithms.SourceDetectionTask( - config=starBackgroundDetectionConfig) - table = afwTable.SourceTable.make(self.initial_stars_schema.schema) - tempDetections = starBackgroundDetectionTask.run( - table=table, exposure=result.exposure, clearMask=True) - result.exposure.mask |= dilatedMask - nFootprintTemp = len(tempDetections.sources) - detected_fraction = self._compute_mask_fraction(result.exposure.mask, detected_mask_planes, - bad_mask_planes) - self.log.info("nIter = %d, thresh = %.2f: Fraction of pixels marked as DETECTED or " - "DETECTED_NEGATIVE in star_background_detection = %.3f " - "(max is %.3f; min is %.3f)", - nIter, starBackgroundDetectionConfig.thresholdValue, - detected_fraction, maxDetFracForFinalBg, minDetFracForFinalBg) - - n_amp = len(result.exposure.detector.getAmplifiers()) - if doCheckPerAmpDetFraction: # detected_fraction < maxDetFracForFinalBg: - n_above_max_per_amp, highest_detected_fraction_per_amp, no_zero_det_amps = \ - self._compute_per_amp_fraction(result.exposure, detected_fraction, - detected_mask_planes, bad_mask_planes) - - if not no_zero_det_amps: - starBackgroundDetectionConfig.thresholdValue = 0.95*currentThresh - nIter += 1 - if nIter > maxIter: - inBackgroundDet = False - - if (detected_fraction < maxDetFracForFinalBg and detected_fraction > minDetFracForFinalBg - and n_above_max_per_amp < int(0.75*n_amp) - and no_zero_det_amps): - if (n_above_max_per_amp < max(1, int(0.15*n_amp)) - or detected_fraction < 0.85*maxDetFracForFinalBg): - inBackgroundDet = False - else: - self.log.warning("Making small tweak....") - starBackgroundDetectionConfig.thresholdValue = 1.05*currentThresh - self.log.warning("n_above_max_per_amp = %d (abs max is %d)", n_above_max_per_amp, int(0.75*n_amp)) - - self.log.info("Fraction of pixels marked as DETECTED or DETECTED_NEGATIVE is now %.5f " - "(highest per amp section = %.5f)", - detected_fraction, highest_detected_fraction_per_amp) - - if detected_fraction > maxDetFracForFinalBg: - result.exposure.mask = dilatedMask - self.log.warning("Final fraction of pixels marked as DETECTED or DETECTED_NEGATIVE " - "was too large in star_background_detection = %.3f (max = %.3f). " - "Reverting to dilated mask from PSF detection...", - detected_fraction, maxDetFracForFinalBg) - star_background = self.star_background.run(exposure=result.exposure).background - result.background.append(star_background[0]) + 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 @@ -1906,11 +1685,4 @@ def _remeasure_star_background(self, result): pedestalBackgroundLevel = pedestalBackground.getImage().array[0, 0] self.log.warning("Subtracted pedestalBackgroundLevel = %.4f", pedestalBackgroundLevel) - # Clear detected mask plane before final round of detection - mask = result.exposure.mask - for mp in detected_mask_planes: - if mp not in mask.getMaskPlaneDict(): - mask.addMaskPlane(mp) - mask &= ~mask.getPlaneBitMask(detected_mask_planes) - return result