diff --git a/python/lsst/ip/diffim/detectAndMeasure.py b/python/lsst/ip/diffim/detectAndMeasure.py index b5d8cbf8d..9fe66e0d9 100644 --- a/python/lsst/ip/diffim/detectAndMeasure.py +++ b/python/lsst/ip/diffim/detectAndMeasure.py @@ -27,7 +27,7 @@ import lsst.afw.table as afwTable import lsst.daf.base as dafBase import lsst.geom -from lsst.ip.diffim.utils import evaluateMaskFraction +from lsst.ip.diffim.utils import evaluateMaskFraction, computeDifferenceImageMetrics from lsst.meas.algorithms import SkyObjectsTask, SourceDetectionTask, SetPrimaryFlagsTask, MaskStreaksTask from lsst.meas.base import ForcedMeasurementTask, ApplyApCorrTask, DetectorVisitIdGeneratorConfig import lsst.meas.deblender @@ -68,6 +68,12 @@ class DetectAndMeasureConnections(pipeBase.PipelineTaskConnections, storageClass="ExposureF", name="{fakesType}{coaddName}Diff_differenceTempExp", ) + kernelSources = pipeBase.connectionTypes.Input( + doc="Final selection of sources used for psf matching.", + dimensions=("instrument", "visit", "detector"), + storageClass="SourceCatalog", + name="{fakesType}{coaddName}Diff_psfMatchSources" + ) outputSchema = pipeBase.connectionTypes.InitOutput( doc="Schema (as an example catalog) for output DIASource catalog.", storageClass="SourceCatalog", @@ -390,7 +396,7 @@ def runQuantum(self, butlerQC: pipeBase.QuantumContext, butlerQC.put(outputs, outputRefs) @timeMethod - def run(self, science, matchedTemplate, difference, + def run(self, science, matchedTemplate, difference, kernelSources, idFactory=None): """Detect and measure sources on a difference image. @@ -411,6 +417,8 @@ def run(self, science, matchedTemplate, difference, difference image. difference : `lsst.afw.image.ExposureF` Result of subtracting template from the science image. + kernelSources : `lsst.afw.table.SourceCatalog` + Final selection of sources that was used for psf matching. idFactory : `lsst.afw.table.IdFactory`, optional Generator object used to assign ids to detected sources in the difference image. Ids from this generator are not set until after @@ -473,7 +481,7 @@ def run(self, science, matchedTemplate, difference, results.negative) result = self.processResults(science, matchedTemplate, difference, - sources, idFactory, + sources, idFactory, kernelSources, positives=positives, negatives=negatives) result.differenceBackground = background @@ -483,7 +491,7 @@ def run(self, science, matchedTemplate, difference, negatives = afwTable.SourceCatalog(self.schema) results.negative.makeSources(negatives) result = self.processResults(science, matchedTemplate, difference, - results.sources, idFactory, + results.sources, idFactory, kernelSources, positives=positives, negatives=negatives) result.differenceBackground = background @@ -514,7 +522,7 @@ def _prepareInputs(self, difference): mask.addMaskPlane(mp) mask &= ~mask.getPlaneBitMask(self.config.clearMaskPlanes) - def processResults(self, science, matchedTemplate, difference, sources, idFactory, + def processResults(self, science, matchedTemplate, difference, sources, idFactory, kernelSources, positives=None, negatives=None,): """Measure and process the results of source detection. @@ -532,6 +540,8 @@ def processResults(self, science, matchedTemplate, difference, sources, idFactor idFactory : `lsst.afw.table.IdFactory` Generator object used to assign ids to detected sources in the difference image. + kernelSources : `lsst.afw.table.SourceCatalog` + Final selection of sources that was used for psf matching. positives : `lsst.afw.table.SourceCatalog`, optional Positive polarity footprints. negatives : `lsst.afw.table.SourceCatalog`, optional @@ -597,7 +607,7 @@ def processResults(self, science, matchedTemplate, difference, sources, idFactor if self.config.doForcedMeasurement: self.measureForcedSources(diaSources, science, difference.getWcs()) - self.calculateMetrics(science, difference, diaSources) + self.calculateMetrics(science, difference, diaSources, kernelSources) measurementResults = pipeBase.Struct( subtractedMeasuredExposure=difference, @@ -787,7 +797,7 @@ def measureForcedSources(self, diaSources, science, wcs): for diaSource, forcedSource in zip(diaSources, forcedSources): diaSource.assign(forcedSource, mapper) - def calculateMetrics(self, science, difference, diaSources): + def calculateMetrics(self, science, difference, diaSources, kernelSources): """Add difference image QA metrics to the Task metadata. This may be used to produce corresponding metrics (see @@ -795,12 +805,14 @@ def calculateMetrics(self, science, difference, diaSources): Parameters ---------- - science : `lsst.afw.image.ExposureF` + science : `lsst.afw.image.ExposureF` Science exposure that was subtracted. difference : `lsst.afw.image.Exposure` The target difference image to calculate metrics for. diaSources : `lsst.afw.table.SourceCatalog` The catalog of detected sources. + kernelSources : `lsst.afw.table.SourceCatalog` + Final selection of sources that was used for psf matching. """ mask = difference.mask badPix = (mask.array & mask.getPlaneBitMask(self.config.detection.excludeMaskPlanes)) > 0 @@ -823,75 +835,20 @@ def calculateMetrics(self, science, difference, diaSources): self.metadata["%s_mask_fraction"%maskPlane.lower()] = -1 self.log.info("Unable to calculate metrics for mask plane %s: not in image"%maskPlane) - residualInfo = self._computeDifferenceResidual(science, difference, diaSources) - self.metadata["residualFootprintRatioMean"] = residualInfo.ratio - self.metadata["residualFootprintRatioStdev"] = residualInfo.ratio_std + if self.config.doSkySources: + skySources = diaSources[diaSources["sky_source"]] + else: + skySources = None + metrics = computeDifferenceImageMetrics(science, difference, kernelSources, sky_sources=skySources) + + self.metadata["residualFootprintRatioMean"] = metrics.differenceFootprintRatioMean + self.metadata["residualFootprintRatioStdev"] = metrics.differenceFootprintRatioStdev + self.metadata["differenceFootprintSkyRatioMean"] = metrics.differenceFootprintSkyRatioMean + self.metadata["differenceFootprintSkyRatioStdev"] = metrics.differenceFootprintSkyRatioStdev self.log.info("Mean, stdev of ratio of difference to science " "pixels in star footprints: %5.4f, %5.4f", - residualInfo.ratio, residualInfo.ratio_std) - - def _computeDifferenceResidual(self, science, difference, diaSources): - """Compute quality metrics (saved to the task metadata) on the - difference image, at the locations of detected stars on the science - image (determined by the DETECTED mask plane). This restricts the metric - to locations that should be well-subtracted. - - Parameters - ---------- - science : `lsst.afw.image.ExposureF` - Science exposure that was subtracted. - difference : `lsst.afw.image.ExposureF` - Result of subtracting template and science. - diaSources : `lsst.afw.table.SourceCatalog` - The catalog of detected sources. - - Returns - ------- - results: `lsst.pipe.base.Struct` - ``ratio`` : `np.ndarray` - Mean of the ratio of the absolute value of the difference image - (with the mean absolute value of the sky regions on the - difference image removed) to the science image, computed in the - footprints of stars detected on the science image. - ``ratio_std`` : `np.ndarray` - Standard Deviation across footprints of the above ratio. - """ - good = science.mask.array & science.mask.getPlaneBitMask('DETECTED') > 0 - badPix = (difference.mask.array - & difference.mask.getPlaneBitMask(self.config.detection.excludeMaskPlanes)) > 0 - good &= ~badPix - if not self.config.doSkySources: - # Place sky sources to sample the background - skySourceFootprints = self.skySources.run(mask=science.mask, seed=difference.info.id, - catalog=afwTable.SourceCatalog(self.schema)) - else: - skySourceFootprints = [sky.getFootprint() for sky in diaSources[diaSources["sky_source"]]] - n = len(skySourceFootprints) - sky_science = np.zeros(n) - sky_difference = np.zeros(n) - for i, footprint in enumerate(skySourceFootprints): - heavy = lsst.afw.detection.makeHeavyFootprint(footprint, science.maskedImage) - heavy_diff = lsst.afw.detection.makeHeavyFootprint(footprint, difference.maskedImage) - sky_science[i] = abs(heavy.getImageArray()).mean() - sky_difference[i] = abs(heavy_diff.getImageArray()).mean() - - # Calculate the difference between the absolute value of the science - # image within detected footprints and the sky sources. - # Then do the same for the difference image - if np.sum(good) > 1: - science_mean = np.abs(np.mean(abs(science.image.array[good])) - sky_science.mean()) - science_std = np.std(abs(science.image.array[good])) - difference_mean = np.abs(np.mean(abs(difference.image.array[good])) - sky_difference.mean()) - difference_std = np.std(abs(difference.image.array[good])) - # Calculate stddev of the science and difference image values and - # propagate to the stddev of the quotient. - # This avoids potential divide by zero errors. - ratio = difference_mean/science_mean - ratio_std = ratio*np.sqrt((difference_std/difference_mean)**2 + (science_std/science_mean)**2) - else: - ratio = np.nan - ratio_std = np.nan - return pipeBase.Struct(ratio=ratio, ratio_std=ratio_std) + self.metadata["residualFootprintRatioMean"], + self.metadata["residualFootprintRatioStdev"]) def _runStreakMasking(self, difference): """Do streak masking and optionally save the resulting streak @@ -995,7 +952,7 @@ class DetectAndMeasureScoreTask(DetectAndMeasureTask): _DefaultName = "detectAndMeasureScore" @timeMethod - def run(self, science, matchedTemplate, difference, scoreExposure, + def run(self, science, matchedTemplate, difference, scoreExposure, kernelSources, idFactory=None): """Detect and measure sources on a score image. @@ -1010,6 +967,8 @@ def run(self, science, matchedTemplate, difference, scoreExposure, Result of subtracting template from the science image. scoreExposure : `lsst.afw.image.ExposureF` Score or maximum likelihood difference image + kernelSources : `lsst.afw.table.SourceCatalog` + Final selection of sources that was used for psf matching. idFactory : `lsst.afw.table.IdFactory`, optional Generator object used to assign ids to detected sources in the difference image. Ids from this generator are not set until after @@ -1047,7 +1006,7 @@ def run(self, science, matchedTemplate, difference, scoreExposure, results.negative) return self.processResults(science, matchedTemplate, difference, - sources, idFactory, + sources, idFactory, kernelSources, positives=positives, negatives=negatives) @@ -1057,6 +1016,6 @@ def run(self, science, matchedTemplate, difference, scoreExposure, negatives = afwTable.SourceCatalog(self.schema) results.negative.makeSources(negatives) return self.processResults(science, matchedTemplate, difference, - results.sources, idFactory, + results.sources, idFactory, kernelSources, positives=positives, negatives=negatives) diff --git a/python/lsst/ip/diffim/subtractImages.py b/python/lsst/ip/diffim/subtractImages.py index f9d916cd7..2a42395df 100644 --- a/python/lsst/ip/diffim/subtractImages.py +++ b/python/lsst/ip/diffim/subtractImages.py @@ -26,7 +26,7 @@ import lsst.afw.image import lsst.afw.math import lsst.geom -from lsst.ip.diffim.utils import evaluateMeanPsfFwhm, getPsfFwhm +from lsst.ip.diffim.utils import evaluateMeanPsfFwhm, getPsfFwhm, computeDifferenceImageMetrics from lsst.meas.algorithms import ScaleVarianceTask, ScienceSourceSelectorTask import lsst.pex.config import lsst.pipe.base @@ -433,88 +433,19 @@ def run(self, template, science, sources, visitSummary=None): # checkTemplateIsSufficient did not raise NoWorkFound, so raise original exception raise e - self.computeImageMetrics(science, subtractResults.difference, sources) + metrics = computeDifferenceImageMetrics(science, subtractResults.difference, sources) - return subtractResults - - def computeImageMetrics(self, science, difference, stars): - r"""Compute quality metrics (saved to the task metadata) on the - difference image, at the locations of detected stars on the science - image. This restricts the metric to locations that should be - well-subtracted. - - Parameters - ---------- - science : `lsst.afw.image.ExposureF` - Science exposure that was subtracted. - difference : `lsst.afw.image.ExposureF` - Result of subtracting template and science. - stars : `lsst.afw.table.SourceCatalog` - Good calibration sources detected on science image; these - footprints are what the metrics are computed on. - - Notes - ----- - The task metadata does not include docstrings, so descriptions of the - computed metrics are given here: - - differenceFootprintRatioMean - Mean of the ratio of the absolute value of the difference image - (with the mean absolute value of the sky regions on the difference - image removed) to the science image, computed in the footprints - of stars detected on the science image (the sums below are of the - pixels in each star or sky footprint): - :math:`\mathrm{mean}_{footprints}((\sum |difference| - - \mathrm{mean}(\sum |difference_{sky}|)) / \sum science)` - differenceFootprintRatioStdev - Standard Deviation across footprints of the above ratio. - differenceFootprintSkyRatioMean - Mean of the ratio of the absolute value of sky source regions on - the difference image to the science image (the sum below is of the - pixels in each sky source footprint): - :math:`\mathrm{mean}_{footprints}(\sum |difference_{sky}| / \sum science_{sky})` - differenceFootprintSkyRatioStdev - Standard Deivation across footprints of the above sky ratio. - """ - def footprint_mean(sources, sky=0): - """Compute ratio of the absolute value of the diffim to the science - image, within each source footprint, subtracting the sky from the - diffim values if provided. - """ - n = len(sources) - science_footprints = np.zeros(n) - difference_footprints = np.zeros(n) - ratio = np.zeros(n) - for i, record in enumerate(sources): - footprint = record.getFootprint() - heavy = lsst.afw.detection.makeHeavyFootprint(footprint, science.maskedImage) - heavy_diff = lsst.afw.detection.makeHeavyFootprint(footprint, difference.maskedImage) - science_footprints[i] = abs(heavy.getImageArray()).mean() - difference_footprints[i] = abs(heavy_diff.getImageArray()).mean() - ratio[i] = abs((difference_footprints[i] - sky)/science_footprints[i]) - return science_footprints, difference_footprints, ratio - - sky = stars["sky_source"] - if np.count_nonzero(sky) > 0: - sky_science, sky_difference, sky_ratio = footprint_mean(stars[sky]) - sky_mean = sky_ratio.mean() - sky_std = sky_ratio.std() - sky_difference = sky_difference.mean() - else: - sky_mean = np.nan - sky_std = np.nan - sky_difference = 0 - science_footprints, difference_footprints, ratio = footprint_mean(stars[~sky], sky_difference) - - self.metadata["differenceFootprintRatioMean"] = ratio.mean() - self.metadata["differenceFootprintRatioStdev"] = ratio.std() - self.metadata["differenceFootprintSkyRatioMean"] = sky_mean - self.metadata["differenceFootprintSkyRatioStdev"] = sky_std + self.metadata["differenceFootprintRatioMean"] = metrics.differenceFootprintRatioMean + self.metadata["differenceFootprintRatioStdev"] = metrics.differenceFootprintRatioStdev + self.metadata["differenceFootprintSkyRatioMean"] = metrics.differenceFootprintSkyRatioMean + self.metadata["differenceFootprintSkyRatioStdev"] = metrics.differenceFootprintSkyRatioStdev self.log.info("Mean, stdev of ratio of difference to science " "pixels in star footprints: %5.4f, %5.4f", self.metadata["differenceFootprintRatioMean"], self.metadata["differenceFootprintRatioStdev"]) + return subtractResults + def runConvolveTemplate(self, template, science, selectSources): """Convolve the template image with a PSF-matching kernel and subtract from the science image. diff --git a/python/lsst/ip/diffim/utils.py b/python/lsst/ip/diffim/utils.py index e52abd9eb..76f251640 100644 --- a/python/lsst/ip/diffim/utils.py +++ b/python/lsst/ip/diffim/utils.py @@ -23,14 +23,17 @@ __all__ = ["evaluateMeanPsfFwhm", "getPsfFwhm", "getKernelCenterDisplacement", + "computeDifferenceImageMetrics", ] import itertools import numpy as np import lsst.geom as geom +import lsst.afw.detection as afwDetection import lsst.afw.image as afwImage import lsst.afw.math as afwMath from lsst.pex.exceptions import InvalidParameterError, RangeError +import lsst.pipe.base from lsst.utils.logging import getLogger @@ -331,3 +334,83 @@ def evaluateMaskFraction(mask, maskPlane): """ nMaskSet = np.count_nonzero((mask.array & mask.getPlaneBitMask(maskPlane))) return nMaskSet/mask.array.size + + +def computeDifferenceImageMetrics(science, difference, stars, sky_sources=None): + r"""Compute quality metrics (saved to the task metadata) on the + difference image, at the locations of detected stars on the science + image. This restricts the metric to locations that should be + well-subtracted. + + Parameters + ---------- + science : `lsst.afw.image.ExposureF` + Science exposure that was subtracted. + difference : `lsst.afw.image.ExposureF` + Result of subtracting template and science. + stars : `lsst.afw.table.SourceCatalog` + Good calibration sources detected on science image; these + footprints are what the metrics are computed on. + + Returns + ------- + metrics : `lsst.pipe.base.Struct` + + ``differenceFootprintRatioMean`` : `float` + Mean of the ratio of the absolute value of the difference image + (with the mean absolute value of the sky regions on the difference + image removed) to the science image, computed in the footprints + of stars detected on the science image (the sums below are of the + pixels in each star or sky footprint): + :math:`\mathrm{mean}_{footprints}((\sum |difference| - + \mathrm{mean}(\sum |difference_{sky}|)) / \sum science)` + ``differenceFootprintRatioStdev`` : `float` + Standard Deviation across footprints of the above ratio. + ``differenceFootprintSkyRatioMean`` : `float` + Mean of the ratio of the absolute value of sky source regions on + the difference image to the science image (the sum below is of the + pixels in each sky source footprint): + :math:`\mathrm{mean}_{footprints}(\sum |difference_{sky}| / \sum science_{sky})` + ``differenceFootprintSkyRatioStdev`` : `float` + Standard Deivation across footprints of the above sky ratio. + """ + def footprint_mean(sources, sky=0): + """Compute ratio of the absolute value of the diffim to the science + image, within each source footprint, subtracting the sky from the + diffim values if provided. + """ + n = len(sources) + science_footprints = np.zeros(n) + difference_footprints = np.zeros(n) + ratio = np.zeros(n) + for i, record in enumerate(sources): + footprint = record.getFootprint() + heavy = afwDetection.makeHeavyFootprint(footprint, science.maskedImage) + heavy_diff = afwDetection.makeHeavyFootprint(footprint, difference.maskedImage) + science_footprints[i] = abs(heavy.getImageArray()).mean() + difference_footprints[i] = abs(heavy_diff.getImageArray()).mean() + ratio[i] = abs((difference_footprints[i] - sky)/science_footprints[i]) + return science_footprints, difference_footprints, ratio + + if "sky_source" in stars.schema: + sky = stars["sky_source"] + selectStars = stars[~sky] + if sky_sources is None: + sky_sources = stars[sky] + else: + selectStars = stars + if sky_sources is not None: + sky_science, sky_difference, sky_ratio = footprint_mean(sky_sources) + sky_mean = sky_ratio.mean() + sky_std = sky_ratio.std() + sky_difference = sky_difference.mean() + else: + sky_mean = np.nan + sky_std = np.nan + sky_difference = 0 + science_footprints, difference_footprints, ratio = footprint_mean(selectStars, sky_difference) + return lsst.pipe.base.Struct(differenceFootprintRatioMean=ratio.mean(), + differenceFootprintRatioStdev=ratio.std(), + differenceFootprintSkyRatioMean=sky_mean, + differenceFootprintSkyRatioStdev=sky_std, + ) diff --git a/tests/test_detectAndMeasure.py b/tests/test_detectAndMeasure.py index 449ef9825..e25d8e0ff 100644 --- a/tests/test_detectAndMeasure.py +++ b/tests/test_detectAndMeasure.py @@ -63,6 +63,8 @@ def _check_diaSource(self, refSources, diaSource, refIds=None, atol : `float`, optional Absolute tolerance of the flux value test. """ + if diaSource["sky_source"]: + return distance = np.sqrt((diaSource.getX() - refSources.getX())**2 + (diaSource.getY() - refSources.getY())**2) self.assertLess(min(distance), matchDistance) @@ -159,7 +161,7 @@ def test_detection_xy0(self): detectionTask = self._setup_detection(doDeblend=True) # Run detection and check the results - output = detectionTask.run(science, matchedTemplate, difference, + output = detectionTask.run(science, matchedTemplate, difference, sources, idFactory=self.idGenerator.make_table_id_factory()) subtractedMeasuredExposure = output.subtractedMeasuredExposure @@ -201,7 +203,7 @@ def test_raise_bad_psf(self): # Run detection and check the results with self.assertRaises(UpstreamFailureNoWorkFound): - detectionTask.run(science, matchedTemplate, difference) + detectionTask.run(science, matchedTemplate, difference, sources) def test_measurements_finite(self): """Measured fluxes and centroids should always be finite. @@ -238,7 +240,7 @@ def test_measurements_finite(self): detectionTask = self._setup_detection(doForcedMeasurement=True) # Run detection and check the results - output = detectionTask.run(science, matchedTemplate, difference) + output = detectionTask.run(science, matchedTemplate, difference, sources) for column in columnNames: self._check_values(output.diaSources[column]) @@ -275,7 +277,7 @@ def test_remove_unphysical(self): badSourceFlags=["base_PixelFlags_flag_offimage", ]) # Run detection and check the results - diaSources = detectionTask.run(science, matchedTemplate, difference).diaSources + diaSources = detectionTask.run(science, matchedTemplate, difference, sources).diaSources badDiaSrcDoRemove = ~bbox.contains(diaSources.getX(), diaSources.getY()) nBadDoRemove = np.count_nonzero(badDiaSrcDoRemove) # Verify that all sources are physical @@ -310,7 +312,7 @@ def test_detect_transients(self): matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs) # Configure the detection Task - detectionTask = self._setup_detection(doMerge=False) + detectionTask = self._setup_detection(doMerge=False, doSkySources=True) kwargs["seed"] = transientSeed kwargs["nSrc"] = 10 kwargs["fluxLevel"] = 1000 @@ -328,7 +330,7 @@ def _detection_wrapper(positive=True): # 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) + output = detectionTask.run(science.clone(), matchedTemplate, difference, sources) refIds = [] scale = 1. if positive else -1. for diaSource in output.diaSources: @@ -377,7 +379,7 @@ def _detection_wrapper(positive=True): # 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) + output = detectionTask.run(science.clone(), matchedTemplate, difference, sources) refIds = [] scale = 1. if positive else -1. for transient in transientSources: @@ -401,7 +403,7 @@ def test_missing_mask_planes(self): detectionTask = self._setup_detection() # Verify that detection runs without errors - detectionTask.run(science, matchedTemplate, difference) + detectionTask.run(science, matchedTemplate, difference, sources) def test_detect_dipoles(self): """Run detection on a difference image containing dipoles. @@ -433,7 +435,7 @@ def test_detect_dipoles(self): difference.maskedImage -= matchedTemplate.maskedImage[science.getBBox()] detectionTask = self._setup_detection(doMerge=True) - output = detectionTask.run(science, matchedTemplate, difference) + output = detectionTask.run(science, matchedTemplate, difference, sources) self.assertEqual(len(output.diaSources), len(sources)) # no sources should be flagged as negative self.assertEqual(len(~output.diaSources["is_negative"]), len(output.diaSources)) @@ -475,11 +477,13 @@ def test_sky_sources(self): detectionTask = self._setup_detection(doSkySources=True) # Run detection and check the results - output = detectionTask.run(science, matchedTemplate, difference, + output = detectionTask.run(science, matchedTemplate, difference, sources, idFactory=self.idGenerator.make_table_id_factory()) skySources = output.diaSources[output.diaSources["sky_source"]] self.assertEqual(len(skySources), detectionTask.config.skySources.nSources) for skySource in skySources: + # Disable the sky_source flag to enable checks. + skySource["sky_source"] = False # The sky sources should not be close to any other source with self.assertRaises(AssertionError): self._check_diaSource(transientSources, skySource, matchDistance=kernelWidth) @@ -527,7 +531,7 @@ 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) + output = detectionTask.run(science, matchedTemplate, difference, sources) refIds = [] goodSrcFlags = checkMask(difference.mask, transientSources, excludeMaskPlanes) if setFlags: @@ -607,7 +611,8 @@ def test_fake_mask_plane_propagation(self): output = detectionTask.run(subtraction.matchedScience, subtraction.matchedTemplate, - subtraction.difference) + subtraction.difference, + subtraction.kernelSources) sci_refIds = [] tmpl_refIds = [] @@ -644,7 +649,7 @@ def test_mask_streaks(self): # Test that no streaks are detected difference = science.clone() difference.maskedImage -= matchedTemplate.maskedImage - output = detectionTask.run(science, matchedTemplate, difference) + output = detectionTask.run(science, matchedTemplate, difference, sources) outMask = output.subtractedMeasuredExposure.mask.array streakMask = output.subtractedMeasuredExposure.mask.getPlaneBitMask("STREAK") streakMaskSet = (outMask & streakMask) > 0 @@ -652,7 +657,7 @@ def test_mask_streaks(self): # Add streak-like shape and check that streak is detected difference.image.array[20:23, 40:200] += 50 - output = detectionTask.run(science, matchedTemplate, difference) + output = detectionTask.run(science, matchedTemplate, difference, sources) outMask = output.subtractedMeasuredExposure.mask.array streakMask = output.subtractedMeasuredExposure.mask.getPlaneBitMask("STREAK") streakMaskSet = (outMask & streakMask) > 0 @@ -677,7 +682,7 @@ def test_mask_streaks(self): streak_image = streak_trapezoid.createImage(bbox) streak_image *= amplitude difference.image.array += streak_image.array - output = detectionTask.run(science, matchedTemplate, difference) + output = detectionTask.run(science, matchedTemplate, difference, sources) outMask = output.subtractedMeasuredExposure.mask.array streakMask = output.subtractedMeasuredExposure.mask.getPlaneBitMask("STREAK") streakMaskSet = (outMask & streakMask) > 0 @@ -710,7 +715,7 @@ def test_detection_xy0(self): detectionTask = self._setup_detection(doDeblend=True) # Run detection and check the results - output = detectionTask.run(science, matchedTemplate, difference, score, + output = detectionTask.run(science, matchedTemplate, difference, score, sources, idFactory=self.idGenerator.make_table_id_factory()) # Catalog ids should be very large from this id generator. @@ -769,7 +774,7 @@ def test_measurements_finite(self): detectionTask = self._setup_detection(doForcedMeasurement=True) # Run detection and check the results - output = detectionTask.run(science, matchedTemplate, difference, score) + output = detectionTask.run(science, matchedTemplate, difference, score, sources) for column in columnNames: self._check_values(output.diaSources[column]) @@ -819,7 +824,7 @@ def _detection_wrapper(positive=True): # 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) + output = detectionTask.run(science.clone(), matchedTemplate, difference, score, sources) refIds = [] scale = 1. if positive else -1. # sources near the edge may have untrustworthy centroids @@ -863,7 +868,7 @@ def test_detect_dipoles(self): score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl) detectionTask = self._setup_detection() - output = detectionTask.run(science, matchedTemplate, difference, score) + output = detectionTask.run(science, matchedTemplate, difference, score, sources) self.assertEqual(len(output.diaSources), len(sources)) refIds = [] for diaSource in output.diaSources: @@ -906,12 +911,14 @@ def test_sky_sources(self): detectionTask = self._setup_detection(doSkySources=True) # Run detection and check the results - output = detectionTask.run(science, matchedTemplate, difference, score, + output = detectionTask.run(science, matchedTemplate, difference, score, sources, idFactory=self.idGenerator.make_table_id_factory()) nSkySourcesGenerated = detectionTask.metadata["n_skySources"] skySources = output.diaSources[output.diaSources["sky_source"]] self.assertEqual(len(skySources), nSkySourcesGenerated) for skySource in skySources: + # Disable the sky_source flag to enable checks. + skySource["sky_source"] = False # The sky sources should not be close to any other source with self.assertRaises(AssertionError): self._check_diaSource(transientSources, skySource, matchDistance=kernelWidth) @@ -962,7 +969,7 @@ 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) + output = detectionTask.run(science, matchedTemplate, difference, score, sources) refIds = [] goodSrcFlags = checkMask(difference.mask, transientSources, excludeMaskPlanes) if setFlags: @@ -1072,23 +1079,23 @@ def testMergeFootprints(self): dataset.addSource(instFlux=4E5, centroid=lsst.geom.Point2D(75, 124)) # negative isolated source on diffim dataset.addSource(instFlux=5E5, centroid=lsst.geom.Point2D(175, 124)) - template, catalog = dataset.realize(noise=100.0, - schema=lsst.meas.base.tests.TestDataset.makeMinimalSchema()) + schema = lsst.meas.base.tests.TestDataset.makeMinimalSchema() + schema.addField("sky_source", "Flag", "Sky source.") + template, catalog = dataset.realize(noise=100.0, schema=schema) # These will be positive sources on the diffim (as noted above). template.image.subset(lsst.geom.Box2I(lsst.geom.Point2I(15, 15), lsst.geom.Point2I(35, 35))).array *= -1 template.image.subset(lsst.geom.Box2I(lsst.geom.Point2I(230, 60), lsst.geom.Point2I(270, 85))).array *= -1 dataset = lsst.meas.base.tests.TestDataset(bbox) - science, _ = dataset.realize(noise=100.0, - schema=lsst.meas.base.tests.TestDataset.makeMinimalSchema()) + science, _ = dataset.realize(noise=100.0, schema=schema) difference = science.clone() difference.image -= template.image config = detectAndMeasure.DetectAndMeasureTask.ConfigClass() config.doDeblend = True task = detectAndMeasure.DetectAndMeasureTask(config=config) - result = task.run(science, template, difference) + result = task.run(science, template, difference, catalog) # The original bug manifested as the (175,124) source disappaering in # the merge step, and being included in the peaks of a duplicated