diff --git a/python/lsst/ip/diffim/detectAndMeasure.py b/python/lsst/ip/diffim/detectAndMeasure.py index 9fe66e0d9..bcba5a2c7 100644 --- a/python/lsst/ip/diffim/detectAndMeasure.py +++ b/python/lsst/ip/diffim/detectAndMeasure.py @@ -45,6 +45,38 @@ "DetectAndMeasureScoreConfig", "DetectAndMeasureScoreTask"] +class BadSubtractionError(pipeBase.AlgorithmError): + """Raised when the residuals in footprints of stars used to compute the + psf-matching kernel exceeds the configured maximum. + """ + def __init__(self, *, ratio, threshold): + msg = ("The ratio of residual power in source footprints on the" + " difference image to the power in the footprints on the" + f" science image was {ratio}, which exceeds the maximum" + f" threshold of {threshold}") + super().__init__(msg) + self.ratio = ratio + self.threshold = threshold + + @property + def metadata(self): + return {"ratio": self.ratio, + "threshold": self.threshold + } + + +class NoDiaSourcesError(pipeBase.AlgorithmError): + """Raised when there are no diaSources detected on an image difference. + """ + def __init__(self): + msg = ("No diaSources detected!") + super().__init__(msg) + + @property + def metadata(self): + return {} + + class DetectAndMeasureConnections(pipeBase.PipelineTaskConnections, dimensions=("instrument", "visit", "detector"), defaultTemplates={"coaddName": "deep", @@ -245,6 +277,32 @@ class DetectAndMeasureConfig(pipeBase.PipelineTaskConfig, doc="Mask planes to clear before running detection.", default=("DETECTED", "DETECTED_NEGATIVE", "NOT_DEBLENDED", "STREAK"), ) + raiseOnBadSubtractionRatio = pexConfig.Field( + dtype=bool, + default=True, + doc="Raise an error if the ratio of power in detected footprints" + " on the difference image to the power in footprints on the science" + " image exceeds ``badSubtractionRatioThreshold``", + ) + badSubtractionRatioThreshold = pexConfig.Field( + dtype=float, + default=0.2, + doc="Maximum ratio of power in footprints on the difference image to" + " the same footprints on the science image." + "Only used if ``raiseOnBadSubtractionRatio`` is set", + ) + badSubtractionVariationThreshold = pexConfig.Field( + dtype=float, + default=0.4, + doc="Maximum standard deviation of the ratio of power in footprints on" + " the difference image to the same footprints on the science image." + "Only used if ``raiseOnBadSubtractionRatio`` is set", + ) + raiseOnNoDiaSources = pexConfig.Field( + dtype=bool, + default=True, + doc="Raise an algorithm error if no diaSources are detected.", + ) idGenerator = DetectorVisitIdGeneratorConfig.make_field() def setDefaults(self): @@ -392,12 +450,32 @@ def runQuantum(self, butlerQC: pipeBase.QuantumContext, inputs = butlerQC.get(inputRefs) idGenerator = self.config.idGenerator.apply(butlerQC.quantum.dataId) idFactory = idGenerator.make_table_id_factory() - outputs = self.run(**inputs, idFactory=idFactory) - butlerQC.put(outputs, outputRefs) + # Specify the fields that `annotate` needs below, to ensure they + # exist, even as None. + measurementResults = pipeBase.Struct( + subtractedMeasuredExposure=None, + diaSources=None, + maskedStreaks=None, + differenceBackground=None, + ) + try: + self.run(**inputs, idFactory=idFactory, measurementResults=measurementResults) + except pipeBase.AlgorithmError as e: + error = pipeBase.AnnotatedPartialOutputsError.annotate( + e, + self, + measurementResults.subtractedMeasuredExposure, + measurementResults.diaSources, + measurementResults.maskedStreaks, + log=self.log + ) + butlerQC.put(measurementResults, outputRefs) + raise error from e + butlerQC.put(measurementResults, outputRefs) @timeMethod def run(self, science, matchedTemplate, difference, kernelSources, - idFactory=None): + idFactory=None, measurementResults=None): """Detect and measure sources on a difference image. The difference image will be convolved with a gaussian approximation of @@ -423,6 +501,10 @@ def run(self, science, matchedTemplate, difference, kernelSources, Generator object used to assign ids to detected sources in the difference image. Ids from this generator are not set until after deblending and merging positive/negative peaks. + measurementResults : `lsst.pipe.base.Struct`, optional + Result struct that is modified to allow saving of partial outputs + for some failure conditions. If the task completes successfully, + this is also returned. Returns ------- @@ -435,6 +517,8 @@ def run(self, science, matchedTemplate, difference, kernelSources, ``differenceBackground`` : `lsst.afw.math.BackgroundList` Background that was subtracted from the difference image. """ + if measurementResults is None: + measurementResults = pipeBase.Struct() if idFactory is None: idFactory = lsst.meas.base.IdGenerator().make_table_id_factory() @@ -474,28 +558,26 @@ def run(self, science, matchedTemplate, difference, kernelSources, doSmooth=True, background=background ) + measurementResults.differenceBackground = background if self.config.doDeblend: sources, positives, negatives = self._deblend(difference, results.positive, results.negative) - result = self.processResults(science, matchedTemplate, difference, - sources, idFactory, kernelSources, - positives=positives, - negatives=negatives) - result.differenceBackground = background else: positives = afwTable.SourceCatalog(self.schema) results.positive.makeSources(positives) negatives = afwTable.SourceCatalog(self.schema) results.negative.makeSources(negatives) - result = self.processResults(science, matchedTemplate, difference, - results.sources, idFactory, kernelSources, - positives=positives, - negatives=negatives) - result.differenceBackground = background - return result + sources = results.sources + + self.processResults(science, matchedTemplate, difference, + sources, idFactory, kernelSources, + positives=positives, + negatives=negatives, + measurementResults=measurementResults) + return measurementResults def _prepareInputs(self, difference): """Ensure that we start with an empty detection and deblended mask. @@ -523,7 +605,7 @@ def _prepareInputs(self, difference): mask &= ~mask.getPlaneBitMask(self.config.clearMaskPlanes) def processResults(self, science, matchedTemplate, difference, sources, idFactory, kernelSources, - positives=None, negatives=None,): + positives=None, negatives=None, measurementResults=None): """Measure and process the results of source detection. Parameters @@ -546,6 +628,10 @@ def processResults(self, science, matchedTemplate, difference, sources, idFactor Positive polarity footprints. negatives : `lsst.afw.table.SourceCatalog`, optional Negative polarity footprints. + measurementResults : `lsst.pipe.base.Struct`, optional + Result struct that is modified to allow saving of partial outputs + for some failure conditions. If the task completes successfully, + this is also returned. Returns ------- @@ -556,6 +642,8 @@ def processResults(self, science, matchedTemplate, difference, sources, idFactor ``diaSources`` : `lsst.afw.table.SourceCatalog` The catalog of detected sources. """ + if measurementResults is None: + measurementResults = pipeBase.Struct() self.metadata["nUnmergedDiaSources"] = len(sources) if self.config.doMerge: # preserve peak schema, if there are any footprints @@ -607,18 +695,22 @@ def processResults(self, science, matchedTemplate, difference, sources, idFactor if self.config.doForcedMeasurement: self.measureForcedSources(diaSources, science, difference.getWcs()) - self.calculateMetrics(science, difference, diaSources, kernelSources) - - measurementResults = pipeBase.Struct( - subtractedMeasuredExposure=difference, - ) - - if len(diaSources) > 0: - measurementResults.diaSources = diaSources + measurementResults.subtractedMeasuredExposure = difference if self.config.doMaskStreaks and self.config.writeStreakInfo: measurementResults.mergeItems(streakInfo, 'maskedStreaks') + self.calculateMetrics(science, difference, diaSources, kernelSources) + + if np.count_nonzero(~diaSources["sky_source"]) > 0: + measurementResults.diaSources = diaSources + elif self.config.raiseOnNoDiaSources: + raise NoDiaSourcesError() + elif len(diaSources) > 0: + # This option allows returning sky sources, + # even if there are no diaSources + measurementResults.diaSources = diaSources + return measurementResults def _deblend(self, difference, positiveFootprints, negativeFootprints): @@ -849,6 +941,13 @@ def calculateMetrics(self, science, difference, diaSources, kernelSources): "pixels in star footprints: %5.4f, %5.4f", self.metadata["residualFootprintRatioMean"], self.metadata["residualFootprintRatioStdev"]) + if self.config.raiseOnBadSubtractionRatio: + if metrics.differenceFootprintRatioMean > self.config.badSubtractionRatioThreshold: + raise BadSubtractionError(ratio=metrics.differenceFootprintRatioMean, + threshold=self.config.badSubtractionRatioThreshold) + if metrics.differenceFootprintRatioStdev > self.config.badSubtractionVariationThreshold: + raise BadSubtractionError(ratio=metrics.differenceFootprintRatioStdev, + threshold=self.config.badSubtractionVariationThreshold) def _runStreakMasking(self, difference): """Do streak masking and optionally save the resulting streak diff --git a/python/lsst/ip/diffim/utils.py b/python/lsst/ip/diffim/utils.py index 76f251640..9648c8f53 100644 --- a/python/lsst/ip/diffim/utils.py +++ b/python/lsst/ip/diffim/utils.py @@ -399,7 +399,8 @@ def footprint_mean(sources, sky=0): sky_sources = stars[sky] else: selectStars = stars - if sky_sources is not None: + # Note that the len() below is only evaluated if sky_sources is not None + if sky_sources is not None and len(sky_sources) > 0: sky_science, sky_difference, sky_ratio = footprint_mean(sky_sources) sky_mean = sky_ratio.mean() sky_std = sky_ratio.std() diff --git a/tests/test_detectAndMeasure.py b/tests/test_detectAndMeasure.py index e25d8e0ff..43af77ff1 100644 --- a/tests/test_detectAndMeasure.py +++ b/tests/test_detectAndMeasure.py @@ -28,7 +28,7 @@ import lsst.geom from lsst.ip.diffim import detectAndMeasure, subtractImages import lsst.meas.algorithms as measAlg -from lsst.pipe.base import InvalidQuantumError, UpstreamFailureNoWorkFound +from lsst.pipe.base import InvalidQuantumError, UpstreamFailureNoWorkFound, AlgorithmError import lsst.utils.tests import lsst.meas.base.tests @@ -101,7 +101,7 @@ def _check_values(self, values, minValue=None, maxValue=None): if maxValue is not None: self.assertTrue(np.all(values <= maxValue)) - def _setup_detection(self, doSkySources=False, nSkySources=5, doSubtractBackground=False, **kwargs): + def _setup_detection(self, doSkySources=True, nSkySources=5, doSubtractBackground=False, **kwargs): """Setup and configure the detection and measurement PipelineTask. Parameters @@ -158,7 +158,10 @@ def test_detection_xy0(self): difference = science.clone() # Configure the detection Task - detectionTask = self._setup_detection(doDeblend=True) + # Set the subtraction residual metric threshold high, since the template + # is not actually subtracted from the science image. + detectionTask = self._setup_detection(doDeblend=True, badSubtractionRatioThreshold=1., + doSkySources=False) # Run detection and check the results output = detectionTask.run(science, matchedTemplate, difference, sources, @@ -268,6 +271,8 @@ def test_remove_unphysical(self): nSrc=1, **kwargs) matchedTemplate, _ = makeTestImage(seed=staticSeed, noiseLevel=noiseLevel/4, noiseSeed=7, nSrc=1, **kwargs) + transients, transientSources = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=8, nSrc=1, **kwargs) + science.maskedImage += transients.maskedImage difference = science.clone() bbox = difference.getBBox() difference.maskedImage -= matchedTemplate.maskedImage @@ -308,7 +313,7 @@ def test_detect_transients(self): transientSeed = 6 fluxLevel = 500 kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel} - science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs) + scienceBase, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs) matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs) # Configure the detection Task @@ -320,17 +325,18 @@ def test_detect_transients(self): # Run detection and check the results def _detection_wrapper(positive=True): transients, transientSources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, **kwargs) - difference = science.clone() - difference.maskedImage -= matchedTemplate.maskedImage + science = scienceBase.clone() if positive: - difference.maskedImage += transients.maskedImage + science.maskedImage += transients.maskedImage else: - difference.maskedImage -= transients.maskedImage + science.maskedImage -= transients.maskedImage + difference = science.clone() + difference.maskedImage -= matchedTemplate.maskedImage # NOTE: NoiseReplacer (run by forcedMeasurement) can modify the # science image if we've e.g. removed parents post-deblending. # Pass a clone of the science image, so that it doesn't disrupt # later tests. - output = detectionTask.run(science.clone(), matchedTemplate, difference, sources) + output = detectionTask.run(science, matchedTemplate, difference, sources) refIds = [] scale = 1. if positive else -1. for diaSource in output.diaSources: @@ -356,8 +362,8 @@ def test_detect_transients_with_background(self): bbox2D = lsst.geom.Box2D(lsst.geom.Point2D(x0, y0), lsst.geom.Extent2D(xSize, ySize)) background_model = afwMath.Chebyshev1Function2D(params, bbox2D) - science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, - background=background_model, **kwargs) + scienceBase, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, + background=background_model, **kwargs) matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs) # Configure the detection Task @@ -369,17 +375,18 @@ def test_detect_transients_with_background(self): # Run detection and check the results def _detection_wrapper(positive=True): transients, transientSources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, **kwargs) - difference = science.clone() - difference.maskedImage -= matchedTemplate.maskedImage + science = scienceBase.clone() if positive: - difference.maskedImage += transients.maskedImage + science.maskedImage += transients.maskedImage else: - difference.maskedImage -= transients.maskedImage + science.maskedImage -= transients.maskedImage + difference = science.clone() + difference.maskedImage -= matchedTemplate.maskedImage # NOTE: NoiseReplacer (run by forcedMeasurement) can modify the # science image if we've e.g. removed parents post-deblending. # Pass a clone of the science image, so that it doesn't disrupt # later tests. - output = detectionTask.run(science.clone(), matchedTemplate, difference, sources) + output = detectionTask.run(science, matchedTemplate, difference, sources) refIds = [] scale = 1. if positive else -1. for transient in transientSources: @@ -400,7 +407,7 @@ def test_missing_mask_planes(self): difference = science.clone() difference.maskedImage -= matchedTemplate.maskedImage - detectionTask = self._setup_detection() + detectionTask = self._setup_detection(raiseOnBadSubtractionRatio=False, raiseOnNoDiaSources=False) # Verify that detection runs without errors detectionTask.run(science, matchedTemplate, difference, sources) @@ -434,13 +441,16 @@ def test_detect_dipoles(self): matchedTemplate.mask.array[...] = np.roll(matchedTemplate.mask.array[...], offset, axis=0) difference.maskedImage -= matchedTemplate.maskedImage[science.getBBox()] - detectionTask = self._setup_detection(doMerge=True) + # Need a higher residual metric threshold, since we are purposely + # creating poor subtractions. + detectionTask = self._setup_detection(doMerge=True, badSubtractionRatioThreshold=0.3) output = detectionTask.run(science, matchedTemplate, difference, sources) - self.assertEqual(len(output.diaSources), len(sources)) + diaSources = output.diaSources[~output.diaSources["sky_source"]].copy(deep=True) + self.assertEqual(len(diaSources), len(sources)) # no sources should be flagged as negative - self.assertEqual(len(~output.diaSources["is_negative"]), len(output.diaSources)) + self.assertEqual(len(~diaSources["is_negative"]), len(diaSources)) refIds = [] - for diaSource in output.diaSources: + for diaSource in diaSources: if diaSource[dipoleFlag]: self._check_diaSource(sources, diaSource, refIds=refIds, scale=0, rtol=0.05, atol=None, usePsfFlux=False) @@ -531,13 +541,14 @@ def _detection_wrapper(setFlags=True): srcBbox = lsst.geom.Box2I(lsst.geom.Point2I(srcX - radius, srcY - radius), lsst.geom.Extent2I(2*radius + 1, 2*radius + 1)) difference[srcBbox].mask.array |= lsst.afw.image.Mask.getPlaneBitMask(badMask) - output = detectionTask.run(science, matchedTemplate, difference, sources) - refIds = [] - goodSrcFlags = checkMask(difference.mask, transientSources, excludeMaskPlanes) if setFlags: - self.assertEqual(np.sum(~goodSrcFlags), nBad) - self.assertFalse(hasattr(output, "diaSources")) + with self.assertRaises(AlgorithmError): + output = detectionTask.run(science, matchedTemplate, difference, sources) + return else: + output = detectionTask.run(science, matchedTemplate, difference, sources) + refIds = [] + goodSrcFlags = checkMask(difference.mask, transientSources, excludeMaskPlanes) self.assertEqual(np.sum(~goodSrcFlags), 0) for diaSource, goodSrcFlag in zip(output.diaSources, goodSrcFlags): if ~goodSrcFlag: @@ -613,10 +624,11 @@ def test_fake_mask_plane_propagation(self): subtraction.matchedTemplate, subtraction.difference, subtraction.kernelSources) + diaSources = output.diaSources[~output.diaSources["sky_source"]].copy(deep=True) sci_refIds = [] tmpl_refIds = [] - for diaSrc in output.diaSources: + for diaSrc in diaSources: if diaSrc['base_PsfFlux_instFlux'] > 0: self._check_diaSource(science_fake_sources, diaSrc, scale=1, refIds=sci_refIds) self.assertTrue(diaSrc['base_PixelFlags_flag_injected']) @@ -644,7 +656,8 @@ def test_mask_streaks(self): matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs) # Configure the detection Task - detectionTask = self._setup_detection(doMerge=False, doMaskStreaks=True) + detectionTask = self._setup_detection(doMerge=False, doMaskStreaks=True, + raiseOnNoDiaSources=False, raiseOnBadSubtractionRatio=False) # Test that no streaks are detected difference = science.clone() @@ -712,7 +725,10 @@ def test_detection_xy0(self): score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl) # Configure the detection Task - detectionTask = self._setup_detection(doDeblend=True) + # Set the subtraction residual metric threshold high, since the template + # is not actually subtracted from the science image. + detectionTask = self._setup_detection(doDeblend=True, badSubtractionRatioThreshold=1., + doSkySources=False) # Run detection and check the results output = detectionTask.run(science, matchedTemplate, difference, score, sources, @@ -791,9 +807,9 @@ def test_detect_transients(self): transientSeed = 6 fluxLevel = 500 kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel} - science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs) + scienceBase, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs) matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs) - scienceKernel = science.psf.getKernel() + scienceKernel = scienceBase.psf.getKernel() subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask() # Configure the detection Task @@ -813,18 +829,19 @@ def _detection_wrapper(positive=True): """ transients, transientSources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, **kwargs) - difference = science.clone() - difference.maskedImage -= matchedTemplate.maskedImage + science = scienceBase.clone() if positive: - difference.maskedImage += transients.maskedImage + science.maskedImage += transients.maskedImage else: - difference.maskedImage -= transients.maskedImage + science.maskedImage -= transients.maskedImage + difference = science.clone() + difference.maskedImage -= matchedTemplate.maskedImage score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl) # NOTE: NoiseReplacer (run by forcedMeasurement) can modify the # science image if we've e.g. removed parents post-deblending. # Pass a clone of the science image, so that it doesn't disrupt # later tests. - output = detectionTask.run(science.clone(), matchedTemplate, difference, score, sources) + output = detectionTask.run(science, matchedTemplate, difference, score, sources) refIds = [] scale = 1. if positive else -1. # sources near the edge may have untrustworthy centroids @@ -867,11 +884,12 @@ def test_detect_dipoles(self): scienceKernel = science.psf.getKernel() score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl) - detectionTask = self._setup_detection() + detectionTask = self._setup_detection(badSubtractionRatioThreshold=0.3) output = detectionTask.run(science, matchedTemplate, difference, score, sources) - self.assertEqual(len(output.diaSources), len(sources)) + diaSources = output.diaSources[~output.diaSources["sky_source"]].copy(deep=True) + self.assertEqual(len(diaSources), len(sources)) refIds = [] - for diaSource in output.diaSources: + for diaSource in diaSources: if diaSource[dipoleFlag]: self._check_diaSource(sources, diaSource, refIds=refIds, scale=0, rtol=0.05, atol=None, usePsfFlux=False) @@ -969,13 +987,14 @@ def _detection_wrapper(setFlags=True): lsst.geom.Extent2I(2*radius + 1, 2*radius + 1)) difference[srcBbox].mask.array |= lsst.afw.image.Mask.getPlaneBitMask(badMask) score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl) - output = detectionTask.run(science, matchedTemplate, difference, score, sources) - refIds = [] - goodSrcFlags = checkMask(difference.mask, transientSources, excludeMaskPlanes) if setFlags: - self.assertEqual(np.sum(~goodSrcFlags), nBad) - self.assertFalse(hasattr(output, "diaSources")) + with self.assertRaises(AlgorithmError): + output = detectionTask.run(science, matchedTemplate, difference, score, sources) + return else: + output = detectionTask.run(science, matchedTemplate, difference, score, sources) + refIds = [] + goodSrcFlags = checkMask(difference.mask, transientSources, excludeMaskPlanes) self.assertEqual(np.sum(~goodSrcFlags), 0) for diaSource, goodSrcFlag in zip(output.diaSources, goodSrcFlags): if ~goodSrcFlag: @@ -1094,6 +1113,8 @@ def testMergeFootprints(self): config = detectAndMeasure.DetectAndMeasureTask.ConfigClass() config.doDeblend = True + config.raiseOnBadSubtractionRatio = False + config.raiseOnNoDiaSources = False task = detectAndMeasure.DetectAndMeasureTask(config=config) result = task.run(science, template, difference, catalog)