From 6d2d2bbd0a6c2cada801abdf40037d48eb9ed53e Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Fri, 1 Mar 2024 09:34:03 -0800 Subject: [PATCH 1/4] Move mask updating code to function --- python/lsst/ip/diffim/subtractImages.py | 125 ++++++++++++++++-------- 1 file changed, 82 insertions(+), 43 deletions(-) diff --git a/python/lsst/ip/diffim/subtractImages.py b/python/lsst/ip/diffim/subtractImages.py index b4d35983..61d65155 100644 --- a/python/lsst/ip/diffim/subtractImages.py +++ b/python/lsst/ip/diffim/subtractImages.py @@ -233,9 +233,15 @@ class AlardLuptonSubtractBaseConfig(lsst.pex.config.Config): ) preserveTemplateMask = lsst.pex.config.ListField( dtype=str, - default=("NO_DATA", "BAD", "SAT", "FAKE", "INJECTED", "INJECTED_CORE"), + default=("NO_DATA", "BAD",), doc="Mask planes from the template to propagate to the image difference." ) + renameTemplateMask = lsst.pex.config.ListField( + dtype=str, + default=("SAT", "INJECTED", "INJECTED_CORE",), + doc="Mask planes from the template to propagate to the image difference" + "with '_TEMPLATE' appended to the name." + ) allowKernelSourceDetection = lsst.pex.config.Field( dtype=bool, default=False, @@ -624,11 +630,6 @@ def finalize(self, template, science, difference, kernel, correctedExposure : `lsst.afw.image.ExposureF` The decorrelated image difference. """ - # Erase existing detection mask planes. - # We don't want the detection mask from the science image - - self.updateMasks(template, science, difference) - if self.config.doDecorrelation: self.log.info("Decorrelating image difference.") # We have cleared the template mask plane, so copy the mask plane of @@ -644,35 +645,6 @@ def finalize(self, template, science, difference, kernel, correctedExposure = difference return correctedExposure - def updateMasks(self, template, science, difference): - """Update the mask planes on images for finalizing.""" - - bbox = science.getBBox() - mask = difference.mask - mask &= ~(mask.getPlaneBitMask("DETECTED") | mask.getPlaneBitMask("DETECTED_NEGATIVE")) - - self.log.info("Adding injected mask planes") - mask.addMaskPlane("INJECTED") - mask.addMaskPlane("INJECTED_TEMPLATE") - - if "FAKE" in science.mask.getMaskPlaneDict().keys(): - # propagate the mask plane related to Fake source injection - # NOTE: the fake source injection sets FAKE plane, but it should be INJECTED - # NOTE: This can be removed in DM-40796 - diffInjectedBitMask = mask.getPlaneBitMask("INJECTED") - diffInjTmpltBitMask = mask.getPlaneBitMask("INJECTED_TEMPLATE") - - scienceFakeBitMask = science.mask.getPlaneBitMask('FAKE') - tmpltFakeBitMask = template[bbox].mask.getPlaneBitMask('FAKE') - - injScienceMaskArray = ((science.mask.array & scienceFakeBitMask) > 0) * diffInjectedBitMask - injTemplateMaskArray = ((template[bbox].mask.array & tmpltFakeBitMask) > 0) * diffInjTmpltBitMask - - mask.array |= injScienceMaskArray - mask.array |= injTemplateMaskArray - - template[bbox].mask.array[...] = difference.mask.array[...] - @staticmethod def _validateExposures(template, science): """Check that the WCS of the two Exposures match, and the template bbox @@ -868,20 +840,87 @@ def _prepareInputs(self, template, science, visitSummary=None): self.metadata.add("scaleTemplateVarianceFactor", templateVarFactor) self.log.info("Science variance scaling factor: %.2f", sciVarFactor) self.metadata.add("scaleScienceVarianceFactor", sciVarFactor) - self._clearMask(template) - def _clearMask(self, template): - """Clear the mask plane of the template. + # Erase existing detection mask planes. + # We don't want the detection mask from the science image + self.updateMasks(template, science) + + def updateMasks(self, template, science): + """Update the science and template mask planes before differencing. Parameters ---------- - template : `lsst.afw.image.ExposureF` + template : `lsst.afw.image.Exposure` Template exposure, warped to match the science exposure. - The mask plane will be modified in place. + The template mask planes will be erased, except for a few specified + in the task config. + science : `lsst.afw.image.Exposure` + Science exposure to subtract from the template. + The DETECTED and DETECTED_NEGATIVE mask planes of the science image + will be erased. + """ + self._clearMask(science.mask, clearMaskPlanes=["DETECTED", "DETECTED_NEGATIVE"]) + + # We will clear ALL template mask planes, except for those specified + # via the `preserveTemplateMask` config. Mask planes specified via + # the `renameTemplateMask` config will be copied to new planes with + # "_TEMPLATE" appended to their names, and the original mask plane will + # be cleared. + clearMaskPlanes = [mp for mp in template.mask.getMaskPlaneDict().keys() + if mp not in self.config.preserveTemplateMask] + renameMaskPlanes = [mp for mp in self.config.renameTemplateMask + if mp in template.mask.getMaskPlaneDict().keys()] + + # propagate the mask plane related to Fake source injection + # NOTE: the fake source injection sets FAKE plane, but it should be INJECTED + # NOTE: This can be removed in DM-40796 + if "FAKE" in science.mask.getMaskPlaneDict().keys(): + self.log.info("Adding injected mask plane to science image") + self._renameMaskPlanes(science.mask, "FAKE", "INJECTED") + if "FAKE" in template.mask.getMaskPlaneDict().keys(): + self.log.info("Adding injected mask plane to template image") + self._renameMaskPlanes(template.mask, "FAKE", "INJECTED_TEMPLATE") + if "INJECTED" in renameMaskPlanes: + renameMaskPlanes.remove("INJECTED") + if "INJECTED_TEMPLATE" in clearMaskPlanes: + clearMaskPlanes.remove("INJECTED_TEMPLATE") + + for maskPlane in renameMaskPlanes: + self._renameMaskPlanes(template.mask, maskPlane, maskPlane + "_TEMPLATE") + self._clearMask(template.mask, clearMaskPlanes=clearMaskPlanes) + + @staticmethod + def _renameMaskPlanes(mask, maskPlane, newMaskPlane): + """Rename a mask plane by adding the new name and copying the data. + + Parameters + ---------- + mask : `lsst.afw.image.Mask` + The mask image to update in place. + maskPlane : `str` + The name of the existing mask plane to copy. + newMaskPlane : `str` + The new name of the mask plane that will be added. + If the mask plane already exists, it will be updated in place. + """ + mask.addMaskPlane(newMaskPlane) + originBitMask = mask.getPlaneBitMask(maskPlane) + destinationBitMask = mask.getPlaneBitMask(newMaskPlane) + mask.array |= ((mask.array & originBitMask) > 0)*destinationBitMask + + def _clearMask(self, mask, clearMaskPlanes=None): + """Clear the mask plane of the template. + + Parameters + ---------- + mask : `lsst.afw.image.Mask` + The mask plane to erase, which will be modified in place. + clearMaskPlanes : `list` of `str`, optional + Erase the specified mask planes. + If not supplied, the entire mask will be erased. """ - mask = template.mask - clearMaskPlanes = [maskplane for maskplane in mask.getMaskPlaneDict().keys() - if maskplane not in self.config.preserveTemplateMask] + if clearMaskPlanes is None: + clearMaskPlanes = list(mask.getMaskPlaneDict().keys()) bitMaskToClear = mask.getPlaneBitMask(clearMaskPlanes) mask &= ~bitMaskToClear From 86e5c15bd6190977f961383a6bc725c5235604ed Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Fri, 1 Mar 2024 17:42:52 -0800 Subject: [PATCH 2/4] Fix unit test --- tests/test_detectAndMeasure.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/tests/test_detectAndMeasure.py b/tests/test_detectAndMeasure.py index c6961c15..31e705f6 100644 --- a/tests/test_detectAndMeasure.py +++ b/tests/test_detectAndMeasure.py @@ -481,6 +481,9 @@ def test_fake_mask_plane_propagation(self): injTmplt_masked = (diff_mask.array & diff_mask.getPlaneBitMask("INJECTED_TEMPLATE")) > 0 self.assertFloatsEqual(inj_masked.astype(int), science_fake_masked.astype(int)) + # The template is convolved, so the INJECTED_TEMPLATE mask plane may + # include more pixels than the FAKE mask plane + injTmplt_masked &= template_fake_masked self.assertFloatsEqual(injTmplt_masked.astype(int), template_fake_masked.astype(int)) # Now check that detection of fakes have the correct flag for injections From 4821d48b3d9348a26daea980a2a341ad9540176d Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Fri, 1 Mar 2024 23:15:55 -0800 Subject: [PATCH 3/4] Support missing STREAK and INJECTED mask planes --- python/lsst/ip/diffim/detectAndMeasure.py | 3 +++ python/lsst/ip/diffim/utils.py | 16 +++++++++----- tests/test_detectAndMeasure.py | 18 ++++++++++++++++ tests/test_subtractTask.py | 26 +++++++++++++++++++++++ 4 files changed, 58 insertions(+), 5 deletions(-) diff --git a/python/lsst/ip/diffim/detectAndMeasure.py b/python/lsst/ip/diffim/detectAndMeasure.py index 73d2a7c5..6f4ab9a9 100644 --- a/python/lsst/ip/diffim/detectAndMeasure.py +++ b/python/lsst/ip/diffim/detectAndMeasure.py @@ -420,6 +420,9 @@ def measureDiaSources(self, diaSources, science, difference, matchedTemplate): Warped and PSF-matched template that was used produce the difference image. """ + # Ensure that the required mask planes are present + for mp in self.config.measurement.plugins["base_PixelFlags"].masksFpAnywhere: + difference.mask.addMaskPlane(mp) # Note that this may not be correct if we convolved the science image. # In the future we may wish to persist the matchedScience image. self.measurement.run(diaSources, difference, science, matchedTemplate) diff --git a/python/lsst/ip/diffim/utils.py b/python/lsst/ip/diffim/utils.py index 8c204319..2d0df3d8 100644 --- a/python/lsst/ip/diffim/utils.py +++ b/python/lsst/ip/diffim/utils.py @@ -1018,7 +1018,7 @@ def computeAveragePsf(exposure: afwImage.Exposure, return psfImage -def detectTestSources(exposure): +def detectTestSources(exposure, addMaskPlanes=None): """Minimal source detection wrapper suitable for unit tests. Parameters @@ -1032,6 +1032,10 @@ def detectTestSources(exposure): selectSources : Source catalog containing candidates """ + if addMaskPlanes is None: + # add empty streak mask plane in lieu of maskStreaksTask + # And add empty INJECTED and INJECTED_TEMPLATE mask planes + addMaskPlanes = ["STREAK", "INJECTED", "INJECTED_TEMPLATE"] schema = afwTable.SourceTable.makeMinimalSchema() selectDetection = measAlg.SourceDetectionTask(schema=schema) @@ -1044,9 +1048,8 @@ def detectTestSources(exposure): sigma=None, # The appropriate sigma is calculated from the PSF doSmooth=True ) - exposure.mask.addMaskPlane("STREAK") # add empty streak mask plane in lieu of maskStreaksTask - exposure.mask.addMaskPlane("INJECTED") # add empty injected mask plane - exposure.mask.addMaskPlane("INJECTED_TEMPLATE") # add empty injected template mask plane + for mp in addMaskPlanes: + exposure.mask.addMaskPlane(mp) selectSources = detRet.sources selectMeasurement.run(measCat=selectSources, exposure=exposure) @@ -1078,6 +1081,7 @@ def makeTestImage(seed=5, nSrc=20, psfSize=2., noiseLevel=5., yLoc=None, flux=None, clearEdgeMask=False, + addMaskPlanes=None, ): """Make a reproduceable PSF-convolved exposure for testing. @@ -1119,6 +1123,8 @@ def makeTestImage(seed=5, nSrc=20, psfSize=2., noiseLevel=5., If specified, must have length equal to ``nSrc`` clearEdgeMask : `bool`, optional Clear the "EDGE" mask plane after source detection. + addMaskPlanes : `list` of `str`, optional + Mask plane names to add to the image. Returns ------- @@ -1172,7 +1178,7 @@ def makeTestImage(seed=5, nSrc=20, psfSize=2., noiseLevel=5., modelExposure.image.array += noise # Run source detection to set up the mask plane - sourceCat = detectTestSources(modelExposure) + sourceCat = detectTestSources(modelExposure, addMaskPlanes=addMaskPlanes) if clearEdgeMask: modelExposure.mask &= ~modelExposure.mask.getPlaneBitMask("EDGE") modelExposure.setPhotoCalib(afwImage.PhotoCalib(calibration, 0., bbox)) diff --git a/tests/test_detectAndMeasure.py b/tests/test_detectAndMeasure.py index 31e705f6..2fbd8c63 100644 --- a/tests/test_detectAndMeasure.py +++ b/tests/test_detectAndMeasure.py @@ -282,6 +282,24 @@ def _detection_wrapper(positive=True): _detection_wrapper(positive=True) _detection_wrapper(positive=False) + def test_missing_mask_planes(self): + """Check that detection runs with missing mask planes. + """ + # Set up the simulated images + noiseLevel = 1. + fluxLevel = 500 + kwargs = {"psfSize": 2.4, "fluxLevel": fluxLevel, "addMaskPlanes": []} + # Use different seeds for the science and template so every source is a diaSource + science, sources = makeTestImage(seed=5, noiseLevel=noiseLevel, noiseSeed=6, **kwargs) + matchedTemplate, _ = makeTestImage(seed=6, noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs) + + difference = science.clone() + difference.maskedImage -= matchedTemplate.maskedImage + detectionTask = self._setup_detection() + + # Verify that detection runs without errors + detectionTask.run(science, matchedTemplate, difference) + def test_detect_dipoles(self): """Run detection on a difference image containing dipoles. """ diff --git a/tests/test_subtractTask.py b/tests/test_subtractTask.py index a8deab16..4266938d 100644 --- a/tests/test_subtractTask.py +++ b/tests/test_subtractTask.py @@ -173,6 +173,32 @@ def test_equal_images(self): makeStats(), statistic=afwMath.STDEV) self.assertFloatsAlmostEqual(differenceStd, np.sqrt(2)*noiseLevel, rtol=0.1) + def test_equal_images_missing_mask_planes(self): + """Test that running with enough sources produces reasonable output, + with the same size psf in the template and science and with missing + mask planes. + """ + noiseLevel = 1. + science, sources = makeTestImage(psfSize=2.4, noiseLevel=noiseLevel, noiseSeed=6, addMaskPlanes=[]) + template, _ = makeTestImage(psfSize=2.4, noiseLevel=noiseLevel, noiseSeed=7, + templateBorderSize=20, doApplyCalibration=True, addMaskPlanes=[]) + config = subtractImages.AlardLuptonSubtractTask.ConfigClass() + config.doSubtractBackground = False + task = subtractImages.AlardLuptonSubtractTask(config=config) + output = task.run(template, science, sources) + # There shoud be no NaN values in the difference image + self.assertTrue(np.all(np.isfinite(output.difference.image.array))) + # Mean of difference image should be close to zero. + meanError = noiseLevel/np.sqrt(output.difference.image.array.size) + # Make sure to include pixels with the DETECTED mask bit set. + statsCtrl = makeStats(badMaskPlanes=("EDGE", "BAD", "NO_DATA", "DETECTED", "DETECTED_NEGATIVE")) + differenceMean = computeRobustStatistics(output.difference.image, output.difference.mask, statsCtrl) + self.assertFloatsAlmostEqual(differenceMean, 0, atol=5*meanError) + # stddev of difference image should be close to expected value. + differenceStd = computeRobustStatistics(output.difference.image, output.difference.mask, + makeStats(), statistic=afwMath.STDEV) + self.assertFloatsAlmostEqual(differenceStd, np.sqrt(2)*noiseLevel, rtol=0.1) + def test_psf_size(self): """Test that the image subtract task runs without failing, if fwhmExposureBuffer and fwhmExposureGrid parameters are set. From bb107ec5e25aa52475370e7392fc42f5b692ce26 Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Mon, 11 Mar 2024 10:35:12 -0700 Subject: [PATCH 4/4] Clean up docstrings --- python/lsst/ip/diffim/detectAndMeasure.py | 16 ++++++++-------- python/lsst/ip/diffim/utils.py | 4 +++- 2 files changed, 11 insertions(+), 9 deletions(-) diff --git a/python/lsst/ip/diffim/detectAndMeasure.py b/python/lsst/ip/diffim/detectAndMeasure.py index 6f4ab9a9..2ba558e1 100644 --- a/python/lsst/ip/diffim/detectAndMeasure.py +++ b/python/lsst/ip/diffim/detectAndMeasure.py @@ -304,14 +304,6 @@ def processResults(self, science, matchedTemplate, difference, sources, table, Parameters ---------- - sources : `lsst.afw.table.SourceCatalog` - Detected sources on the difference exposure. - positiveFootprints : `lsst.afw.detection.FootprintSet`, optional - Positive polarity footprints. - negativeFootprints : `lsst.afw.detection.FootprintSet`, optional - Negative polarity footprints. - table : `lsst.afw.table.SourceTable` - Table object that will be used to create the SourceCatalog. science : `lsst.afw.image.ExposureF` Science exposure that the template was subtracted from. matchedTemplate : `lsst.afw.image.ExposureF` @@ -319,6 +311,14 @@ def processResults(self, science, matchedTemplate, difference, sources, table, difference image. difference : `lsst.afw.image.ExposureF` Result of subtracting template from the science image. + sources : `lsst.afw.table.SourceCatalog` + Detected sources on the difference exposure. + table : `lsst.afw.table.SourceTable` + Table object that will be used to create the SourceCatalog. + positiveFootprints : `lsst.afw.detection.FootprintSet`, optional + Positive polarity footprints. + negativeFootprints : `lsst.afw.detection.FootprintSet`, optional + Negative polarity footprints. Returns ------- diff --git a/python/lsst/ip/diffim/utils.py b/python/lsst/ip/diffim/utils.py index 2d0df3d8..27ddc6af 100644 --- a/python/lsst/ip/diffim/utils.py +++ b/python/lsst/ip/diffim/utils.py @@ -1026,10 +1026,12 @@ def detectTestSources(exposure, addMaskPlanes=None): exposure : `lsst.afw.image.Exposure` Exposure on which to run detection/measurement The exposure is modified in place to set the 'DETECTED' mask plane. + addMaskPlanes : `list` of `str`, optional + Additional mask planes to add to the maskedImage of the exposure. Returns ------- - selectSources : + selectSources Source catalog containing candidates """ if addMaskPlanes is None: