Skip to content

Commit

Permalink
Deblend difference image sources
Browse files Browse the repository at this point in the history 
We were detecting footprints of blended sources, but our measurements
were of the blend, not the individual sources.

* Slightly increase test threshold; we're now comparing with the truth
catalog, so might expect a slightly worse match.
* Don't re-use the science exposure when running detectAndMeasure:
noiseReplacer doesn't handle our "no parents" deblended catalog well
when restoring the image.
* Reject sources near the edge: the existing test was too fragile to
small changes in source centroids.
  • Loading branch information
@parejkoj
parejkoj committed Mar 7, 2024
1 parent a9087ad commit 21b8f0c
Show file tree
Hide file tree
Showing 2 changed files with 104 additions and 17 deletions.
98 changes: 90 additions & 8 deletions python/lsst/ip/diffim/detectAndMeasure.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,9 +21,10 @@

import numpy as np

import lsst.afw.detection as afwDetection
import lsst.afw.table as afwTable
import lsst.daf.base as dafBase
from lsst.meas.algorithms import SkyObjectsTask, SourceDetectionTask
from lsst.meas.algorithms import SkyObjectsTask, SourceDetectionTask, SetPrimaryFlagsTask
from lsst.meas.base import ForcedMeasurementTask, ApplyApCorrTask, DetectorVisitIdGeneratorConfig
import lsst.meas.extensions.trailedSources # noqa: F401
import lsst.meas.extensions.shapeHSM
Expand Down Expand Up @@ -103,6 +104,10 @@ class DetectAndMeasureConfig(pipeBase.PipelineTaskConfig,
target=SourceDetectionTask,
doc="Final source detection for diaSource measurement",
)
deblend = pexConfig.ConfigurableField(
target=lsst.meas.deblender.SourceDeblendTask,
doc="Task to split blended sources into their components."
)
measurement = pexConfig.ConfigurableField(
target=DipoleFitTask,
doc="Task to measure sources on the difference image.",
Expand Down Expand Up @@ -139,6 +144,10 @@ class DetectAndMeasureConfig(pipeBase.PipelineTaskConfig,
target=SkyObjectsTask,
doc="Generate sky sources",
)
setPrimaryFlags = pexConfig.ConfigurableField(
target=SetPrimaryFlagsTask,
doc="Task to add isPrimary and deblending-related flags to the catalog."
)
badSourceFlags = lsst.pex.config.ListField(
dtype=str,
doc="Sources with any of these flags set are removed before writing the output catalog.",
Expand Down Expand Up @@ -198,6 +207,8 @@ def __init__(self, **kwargs):

self.algMetadata = dafBase.PropertyList()
self.makeSubtask("detection", schema=self.schema)
self.makeSubtask("deblend", schema=self.schema)
self.makeSubtask("setPrimaryFlags", schema=self.schema, isSingleFrame=True)
self.makeSubtask("measurement", schema=self.schema,
algMetadata=self.algMetadata)
if self.config.doApCorr:
Expand Down Expand Up @@ -283,20 +294,29 @@ def run(self, science, matchedTemplate, difference,
``diaSources`` : `lsst.afw.table.SourceCatalog`
The catalog of detected sources.
"""
if idFactory is None:
idFactory = lsst.meas.base.IdGenerator().make_table_id_factory()

# Ensure that we start with an empty detection mask.
mask = difference.mask
mask &= ~(mask.getPlaneBitMask("DETECTED") | mask.getPlaneBitMask("DETECTED_NEGATIVE"))

table = afwTable.SourceTable.make(self.schema, idFactory)
table = afwTable.SourceTable.make(self.schema)
table.setMetadata(self.algMetadata)
results = self.detection.run(
table=table,
exposure=difference,
doSmooth=True,
)

return self.processResults(science, matchedTemplate, difference, results.sources, table,
positiveFootprints=results.positive, negativeFootprints=results.negative)
sources, positives, negatives = self._deblend(difference,
results.positive,
results.negative,
idFactory)

return self.processResults(science, matchedTemplate, difference, sources, table,
positiveFootprints=positives,
negativeFootprints=negatives)

def processResults(self, science, matchedTemplate, difference, sources, table,
positiveFootprints=None, negativeFootprints=None,):
Expand Down Expand Up @@ -358,6 +378,64 @@ def processResults(self, science, matchedTemplate, difference, sources, table,

return measurementResults

def _deblend(self, difference, positiveFootprints, negativeFootprints, idFactory):
"""Deblend the positive and negative footprints and return a catalog
containing just the children, and the deblended footprints.
Parameters
----------
difference : `lsst.afw.image.Exposure`
Result of subtracting template from the science image.
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 to assign ids to detected sources.
Returns
-------
sources : `lsst.afw.table.SourceCatalog`
Positive and negative deblended children.
positives, negatives : `lsst.afw.detection.FootprintSet`
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

table = afwTable.SourceTable.make(self.schema, idFactory)
table.setMetadata(self.algMetadata)

def deblend(footprints):
"""Deblend a positive or negative footprint set,
and return the deblended children.
"""
sources = afwTable.SourceCatalog(table)
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)

# TODO: id generation here might be a problem; we are reusing the same
# id generator, which might just keep adding to the next id, instead
# of starting from scratch for the final catalog.
# TODO solution?: just use idFactory post-deblending, not in initial detection?
# TODO: also do that in CalibrateImage?
sources = afwTable.SourceCatalog(table)
sources.reserve(len(positives) + len(negatives))
sources.extend(positives, deep=True)
sources.extend(negatives, deep=True)
return sources, makeFootprints(positives), makeFootprints(negatives)

def _removeBadSources(self, diaSources):
"""Remove bad diaSources from the catalog.
Expand Down Expand Up @@ -447,8 +525,7 @@ def measureForcedSources(self, diaSources, science, wcs):
"""
# Run forced psf photometry on the PVI at the diaSource locations.
# Copy the measured flux and error into the diaSource.
forcedSources = self.forcedMeasurement.generateMeasCat(
science, diaSources, wcs)
forcedSources = self.forcedMeasurement.generateMeasCat(science, diaSources, wcs)
self.forcedMeasurement.run(forcedSources, science, diaSources, wcs)
mapper = afwTable.SchemaMapper(forcedSources.schema, diaSources.schema)
mapper.addMapping(forcedSources.schema.find("base_PsfFlux_instFlux")[0],
Expand Down Expand Up @@ -559,5 +636,10 @@ 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())

return self.processResults(science, matchedTemplate, difference, results.sources, table,
positiveFootprints=results.positive, negativeFootprints=results.negative)
sources, positives, negatives = self._deblend(difference,
results.positive,
results.negative,
idFactory)

return self.processResults(science, matchedTemplate, difference, sources, table,
positiveFootprints=positives, negativeFootprints=negatives)
23 changes: 14 additions & 9 deletions tests/test_detectAndMeasure.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,7 @@ class DetectAndMeasureTestBase:

def _check_diaSource(self, refSources, diaSource, refIds=None,
matchDistance=1., scale=1., usePsfFlux=True,
rtol=0.02, atol=None):
rtol=0.021, atol=None):
"""Match a diaSource with a source in a reference catalog
and compare properties.
Expand Down Expand Up @@ -274,7 +274,11 @@ def _detection_wrapper(positive=True):
difference.maskedImage += transients.maskedImage
else:
difference.maskedImage -= transients.maskedImage
output = detectionTask.run(science, matchedTemplate, difference)
# 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)
refIds = []
scale = 1. if positive else -1.
for diaSource in output.diaSources:
Expand Down Expand Up @@ -622,16 +626,17 @@ def _detection_wrapper(positive=True):
else:
difference.maskedImage -= transients.maskedImage
score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
output = detectionTask.run(science, matchedTemplate, difference, score)
# 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)
refIds = []
scale = 1. if positive else -1.
goodSrcFlags = subtractTask._checkMask(score.mask, transientSources,
subtractTask.config.badMaskPlanes)
# sources near the edge may have untrustworthy centroids
goodSrcFlags = ~output.diaSources['base_PixelFlags_flag_edge']
for diaSource, goodSrcFlag in zip(output.diaSources, goodSrcFlags):
if ~goodSrcFlag:
with self.assertRaises(AssertionError):
self._check_diaSource(transientSources, diaSource, refIds=refIds, scale=scale)
else:
if goodSrcFlag:
self._check_diaSource(transientSources, diaSource, refIds=refIds, scale=scale)
_detection_wrapper(positive=True)
_detection_wrapper(positive=False)
Expand Down

0 comments on commit 21b8f0c

Please sign in to comment.