diff --git a/python/lsst/ip/diffim/detectAndMeasure.py b/python/lsst/ip/diffim/detectAndMeasure.py index 61113a9d7..535e03836 100644 --- a/python/lsst/ip/diffim/detectAndMeasure.py +++ b/python/lsst/ip/diffim/detectAndMeasure.py @@ -21,7 +21,6 @@ import numpy as np -import lsst.afw.detection as afwDetection import lsst.afw.table as afwTable import lsst.daf.base as dafBase import lsst.geom @@ -77,6 +76,12 @@ class DetectAndMeasureConnections(pipeBase.PipelineTaskConnections, storageClass="SourceCatalog", name="{fakesType}{coaddName}Diff_diaSrc", ) + removedSources = pipeBase.connectionTypes.Output( + doc="Sources and blends removed from the diaSources catalog.", + dimensions=("instrument", "visit", "detector"), + storageClass="SourceCatalog", + name="{fakesType}{coaddName}Diff_diaSrc_removed", + ) subtractedMeasuredExposure = pipeBase.connectionTypes.Output( doc="Difference image with detection mask plane filled in.", dimensions=("instrument", "visit", "detector"), @@ -104,7 +109,8 @@ class DetectAndMeasureConfig(pipeBase.PipelineTaskConfig, dtype=bool, default=True, doc="Merge positive and negative diaSources with grow radius " - "set by growFootprint" + "set by growFootprint", + deprecated="This field is no longer used and will be removed after v28." ) doForcedMeasurement = pexConfig.Field( dtype=bool, @@ -143,7 +149,8 @@ class DetectAndMeasureConfig(pipeBase.PipelineTaskConfig, growFootprint = pexConfig.Field( dtype=int, default=2, - doc="Grow positive and negative footprints by this many pixels before merging" + doc="Grow positive and negative footprints by this many pixels before merging", + deprecated="This field is no longer used and will be removed after v28." ) diaSourceMatchRadius = pexConfig.Field( dtype=float, @@ -350,16 +357,11 @@ def run(self, science, matchedTemplate, difference, doSmooth=True, ) - sources, positives, negatives = self._deblend(difference, - results.positive, - results.negative) + sources = self._buildCatalogAndDeblend(difference, results.positive, results.negative, idFactory) - return self.processResults(science, matchedTemplate, difference, sources, idFactory, - positiveFootprints=positives, - negativeFootprints=negatives) + return self._measureSources(science, matchedTemplate, difference, sources) - def processResults(self, science, matchedTemplate, difference, sources, idFactory, - positiveFootprints=None, negativeFootprints=None,): + def _measureSources(self, science, matchedTemplate, difference, initialDiaSources): """Measure and process the results of source detection. Parameters @@ -371,15 +373,8 @@ def processResults(self, science, matchedTemplate, difference, sources, idFactor difference image. difference : `lsst.afw.image.ExposureF` Result of subtracting template from the science image. - sources : `lsst.afw.table.SourceCatalog` + initialDiaSources : `lsst.afw.table.SourceCatalog` Detected sources on the difference exposure. - idFactory : `lsst.afw.table.IdFactory` - Generator object used to assign ids to detected sources in the - difference image. - positiveFootprints : `lsst.afw.detection.FootprintSet`, optional - Positive polarity footprints. - negativeFootprints : `lsst.afw.detection.FootprintSet`, optional - Negative polarity footprints. Returns ------- @@ -390,38 +385,14 @@ def processResults(self, science, matchedTemplate, difference, sources, idFactor ``diaSources`` : `lsst.afw.table.SourceCatalog` The catalog of detected sources. """ - self.metadata.add("nUnmergedDiaSources", len(sources)) - if self.config.doMerge: - fpSet = positiveFootprints - fpSet.merge(negativeFootprints, self.config.growFootprint, - self.config.growFootprint, False) - initialDiaSources = afwTable.SourceCatalog(self.schema) - fpSet.makeSources(initialDiaSources) - self.log.info("Merging detections into %d sources", len(initialDiaSources)) - else: - initialDiaSources = sources - - # Assign source ids at the end: deblend/merge mean that we don't keep - # track of parents and children, we only care about the final ids. - for source in initialDiaSources: - source.setId(idFactory()) - # Ensure sources added after this get correct ids. - initialDiaSources.getTable().setIdFactory(idFactory) - initialDiaSources.setMetadata(self.algMetadata) - self.metadata.add("nMergedDiaSources", len(initialDiaSources)) + initialDiaSources.setMetadata(self.algMetadata) if self.config.doMaskStreaks: streakInfo = self._runStreakMasking(difference.maskedImage) - if self.config.doSkySources: - self.addSkySources(initialDiaSources, difference.mask, difference.info.id) - - if not initialDiaSources.isContiguous(): - initialDiaSources = initialDiaSources.copy(deep=True) - self.measureDiaSources(initialDiaSources, science, difference, matchedTemplate) - diaSources = self._removeBadSources(initialDiaSources) + diaSources, removedSources = self._removeBadSources(initialDiaSources) if self.config.doForcedMeasurement: self.measureForcedSources(diaSources, science, difference.getWcs()) @@ -431,15 +402,21 @@ def processResults(self, science, matchedTemplate, difference, sources, idFactor measurementResults = pipeBase.Struct( subtractedMeasuredExposure=difference, diaSources=diaSources, + removedSources=removedSources, ) if self.config.doMaskStreaks and self.config.writeStreakInfo: measurementResults.mergeItems(streakInfo, 'maskedStreaks') return measurementResults - def _deblend(self, difference, positiveFootprints, negativeFootprints): - """Deblend the positive and negative footprints and return a catalog - containing just the children, and the deblended footprints. + def _buildCatalogAndDeblend(self, difference, positiveFootprints, negativeFootprints, idFactory): + """Create a source catalog and deblend when possible. + + This method creates a source catalog from the positive and negative + footprints, and then deblends the footprints that have only positive + peaks. This is because `lsst.meas.deblender.SourceDeblendTask` is not + designed to handle footprints that have negative flux, resulting in + garbage output. Parameters ---------- @@ -448,6 +425,9 @@ def _deblend(self, difference, positiveFootprints, negativeFootprints): positiveFootprints, negativeFootprints : `lsst.afw.detection.FootprintSet` Positive and negative polarity footprints measured on ``difference`` to be deblended separately. + idFactory : `lsst.afw.table.IdFactory` + Generator object used to assign ids to detected sources in the + difference image. Returns ------- @@ -457,33 +437,67 @@ def _deblend(self, difference, positiveFootprints, negativeFootprints): Deblended positive and negative polarity footprints measured on ``difference``. """ - def makeFootprints(sources): - footprints = afwDetection.FootprintSet(difference.getBBox()) - footprints.setFootprints([src.getFootprint() for src in sources]) - return footprints - - def deblend(footprints): - """Deblend a positive or negative footprint set, - and return the deblended children. - """ - sources = afwTable.SourceCatalog(self.schema) - footprints.makeSources(sources) - self.deblend.run(exposure=difference, sources=sources) - self.setPrimaryFlags.run(sources) - children = sources["detect_isDeblendedSource"] == 1 - sources = sources[children].copy(deep=True) - # Clear parents, so that measurement plugins behave correctly. - sources['parent'] = 0 - return sources.copy(deep=True) - - positives = deblend(positiveFootprints) - negatives = deblend(negativeFootprints) - - sources = afwTable.SourceCatalog(self.schema) - sources.reserve(len(positives) + len(negatives)) - sources.extend(positives, deep=True) - sources.extend(negatives, deep=True) - return sources, makeFootprints(positives), makeFootprints(negatives) + self.metadata.add( + "nUnmergedDiaSources", + len(positiveFootprints.getFootprints()) + len(negativeFootprints.getFootprints()) + ) + # Merge the positive and negative footprints. + # The original detection FootprintSets already grew each detection, + # so there is no need to grow them again. + merged_footprints = positiveFootprints + merged_footprints.merge(negativeFootprints, 0, 0, False) + + # Create a source catalog from the footprints. + table = afwTable.SourceTable.make(self.schema, idFactory) + sources = afwTable.SourceCatalog(table) + merged_footprints.makeSources(sources) + + # Sky sources must be added before deblending, otherwise the + # sources with parent == 0 will be out of order and + # downstream measurement tasks cannot run. + if self.config.doSkySources: + self.addSkySources(sources, difference.mask, difference.info.id) + + # Find the footprints with only positive peaks and no negative peaks. + footprints = [src.getFootprint() for src in sources] + nPeaks = np.array([len(fp.peaks) for fp in footprints]) + blend_ids = [] + skipped_ids = [] + for src in sources[nPeaks > 1]: + peaks = src.getFootprint().peaks + peak_flux = np.array([peak.getPeakValue() for peak in peaks]) + if np.all(peak_flux > 0) or np.all(peak_flux < 0): + blend_ids.append(src.getId()) + else: + skipped_ids.append(src.getId()) + + # Deblend the footprints that can be deblended with SourceDeblendTask. + temp_cat = afwTable.SourceCatalog(sources.table) + for sid in blend_ids: + temp_cat.append(sources.find(sid)) + first_child_index = len(temp_cat) + self.deblend.run(exposure=difference, sources=temp_cat) + + # Update the deblended parents and add their children + # to the sources catalog. + sources.extend(temp_cat[first_child_index:], deep=False) + + # Set deblending flags. + # Since SourceDeblendTask already has a method for setting + # flags for skipped parents, we just call that method here + # instead of setting the flags manually. + for sid in skipped_ids: + src = sources.find(sid) + self.deblend.skipParent(src, difference.mask) + + # Set detection and primary flags + self.setPrimaryFlags.run(sources) + + table = afwTable.SourceTable.make(self.schema, idFactory) + _sources = afwTable.SourceCatalog(table) + _sources.extend(sources, deep=True) + + return _sources def _removeBadSources(self, diaSources): """Remove unphysical diaSources from the catalog. @@ -499,17 +513,33 @@ def _removeBadSources(self, diaSources): The updated catalog of detected sources, with any source that has a flag in ``config.badSourceFlags`` set removed. """ - selector = np.ones(len(diaSources), dtype=bool) + bad_sources = np.zeros(len(diaSources), dtype=bool) for flag in self.config.badSourceFlags: flags = diaSources[flag] nBad = np.count_nonzero(flags) if nBad > 0: self.log.debug("Found %d unphysical sources with flag %s.", nBad, flag) - selector &= ~flags - nBadTotal = np.count_nonzero(~selector) + bad_sources |= flags + # Remove parents of bad children, and children of bad parents. + # This is due to the measurement plugins, where it is assumed that + # parent blends have all of their children in the catalog, and vice versa. + parents = ( + (diaSources["parent"] == 0) + & (diaSources["deblend_nChild"] > 1) + & (~diaSources["deblend_skipped"]) + ) + children = (diaSources["parent"] != 0) + parents_to_remove = np.unique(diaSources["parent"][children & bad_sources]) + bad_parents = parents & bad_sources + bad_parents |= np.in1d(diaSources["id"], parents_to_remove) + children_to_remove = np.in1d(diaSources["parent"], diaSources["id"][bad_parents]) + bad_children = (bad_sources & children) | children_to_remove + bad_sources |= bad_parents | bad_children + + nBadTotal = np.count_nonzero(bad_sources) self.metadata.add("nRemovedBadFlaggedSources", nBadTotal) - self.log.info("Removed %d unphysical sources.", nBadTotal) - return diaSources[selector].copy(deep=True) + self.log.info("Removed %d unphysical sources and blends.", nBadTotal) + return diaSources[~bad_sources].copy(deep=True), diaSources[bad_sources] def addSkySources(self, diaSources, mask, seed, subtask=None): @@ -731,9 +761,6 @@ def run(self, science, matchedTemplate, difference, scoreExposure, # Copy the detection mask from the Score image to the difference image difference.mask.assign(scoreExposure.mask, scoreExposure.getBBox()) - sources, positives, negatives = self._deblend(difference, - results.positive, - results.negative) + sources = self._buildCatalogAndDeblend(difference, results.positive, results.negative, idFactory) - return self.processResults(science, matchedTemplate, difference, sources, idFactory, - positiveFootprints=positives, negativeFootprints=negatives) + return self._measureSources(science, matchedTemplate, difference, sources) diff --git a/tests/test_detectAndMeasure.py b/tests/test_detectAndMeasure.py index 48e46f7e9..0932641b8 100644 --- a/tests/test_detectAndMeasure.py +++ b/tests/test_detectAndMeasure.py @@ -247,8 +247,9 @@ def test_remove_unphysical(self): badDiaSrc = ~bbox.contains(diaSources.getX(), diaSources.getY()) nBad = np.count_nonzero(badDiaSrc) self.assertEqual(nBad, nSetBad) - diaSourcesNoBad = detectionTask._removeBadSources(diaSources) + diaSourcesNoBad, removedSources = detectionTask._removeBadSources(diaSources) badDiaSrcNoBad = ~bbox.contains(diaSourcesNoBad.getX(), diaSourcesNoBad.getY()) + self.assertEqual(len(removedSources), nSetBad) # Verify that no sources outside the image bounding box remain self.assertEqual(np.count_nonzero(badDiaSrcNoBad), 0) @@ -288,7 +289,9 @@ def _detection_wrapper(positive=True): output = detectionTask.run(science.clone(), matchedTemplate, difference) refIds = [] scale = 1. if positive else -1. - for diaSource in output.diaSources: + # Deblending provides multiple copies of the main parent source, + # so we only check the parents. Deblending is checked in another test. + for diaSource in output.diaSources[output.diaSources["parent"] == 0]: self._check_diaSource(transientSources, diaSource, refIds=refIds, scale=scale) _detection_wrapper(positive=True) _detection_wrapper(positive=False) @@ -517,7 +520,7 @@ def test_fake_mask_plane_propagation(self): sci_refIds = [] tmpl_refIds = [] - for diaSrc in output.diaSources: + for diaSrc in output.diaSources[output.diaSources["parent"] == 0]: 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']) @@ -687,7 +690,10 @@ def _detection_wrapper(positive=True): scale = 1. if positive else -1. # sources near the edge may have untrustworthy centroids goodSrcFlags = ~output.diaSources['base_PixelFlags_flag_edge'] - for diaSource, goodSrcFlag in zip(output.diaSources, goodSrcFlags): + # Deblending provides multiple copies of the main parent source, + # so we only check the parents. Deblending is checked in another test. + for diaSource, goodSrcFlag in zip(output.diaSources[output.diaSources["parent"] == 0], + goodSrcFlags): if goodSrcFlag: self._check_diaSource(transientSources, diaSource, refIds=refIds, scale=scale) _detection_wrapper(positive=True)