diff --git a/python/lsst/ip/diffim/getTemplate.py b/python/lsst/ip/diffim/getTemplate.py index 02ba35d5..c8f3059c 100644 --- a/python/lsst/ip/diffim/getTemplate.py +++ b/python/lsst/ip/diffim/getTemplate.py @@ -34,7 +34,7 @@ from lsst.skymap import BaseSkyMap from lsst.ip.diffim.dcrModel import DcrModel -from lsst.meas.algorithms import CoaddPsf, CoaddPsfConfig, SubtractBackgroundTask +from lsst.meas.algorithms import CoaddPsf, CoaddPsfConfig, SubtractBackgroundTask, ScaleVarianceTask from lsst.utils.timer import timeMethod __all__ = [ @@ -120,6 +120,15 @@ class GetTemplateConfig( doc="Minimum fraction of unmasked pixels needed to set the" " HIGH_VARIANCE mask plane.", ) + doScaleVariance = pexConfig.Field( + dtype=bool, + default=True, + doc="Scale variance of the template image?" + ) + scaleVariance = pexConfig.ConfigurableField( + target=ScaleVarianceTask, + doc="Subtask to rescale the variance of the template to the statistically expected level." + ) def setDefaults(self): # Use a smaller cache: per SeparableKernel.computeCache, this should @@ -152,6 +161,8 @@ class GetTemplateTask(pipeBase.PipelineTask): def __init__(self, *args, **kwargs): super().__init__(*args, **kwargs) + if self.config.doScaleVariance: + self.makeSubtask("scaleVariance") self.warper = afwMath.Warper.fromConfig(self.config.warp) self.schema = afwTable.ExposureTable.makeMinimalSchema() self.schema.addField( @@ -381,6 +392,14 @@ def run(self, *, coaddExposureHandles, bbox, wcs, dataIds, physical_filter, visi if count == 0: raise pipeBase.NoWorkFound("No valid pixels in warped template.") + if self.config.doScaleVariance: + # Scale the variance of the template image before subtraction, if + # needed. Note that the science variance is scaled + # independently in ``AlardLuptonSubtractTask``. + varianceFactor = self.scaleVariance.run(template.maskedImage) + self.log.info("Template variance scaling factor: %.2f", varianceFactor) + self.metadata["scaleTemplateVarianceFactor"] = varianceFactor + # Make a single catalog containing all the inputs that were accepted. catalog = afwTable.ExposureCatalog(self.schema) catalog.reserve(sum([len(c) for c in catalogs])) @@ -394,6 +413,7 @@ def run(self, *, coaddExposureHandles, bbox, wcs, dataIds, physical_filter, visi template.setFilter(afwImage.FilterLabel(band, physical_filter)) template.setPhotoCalib(photoCalib) template.setPsf(self._makePsf(template, catalog, wcs)) + # Record the input coadd patches as the template's coadd inputs. coaddInputs = afwImage.CoaddInputs(afwTable.ExposureTable.makeMinimalSchema(), self.schema) coaddInputs.ccds.extend(catalog, deep=True) diff --git a/python/lsst/ip/diffim/subtractImages.py b/python/lsst/ip/diffim/subtractImages.py index 10455c1c..09b48301 100644 --- a/python/lsst/ip/diffim/subtractImages.py +++ b/python/lsst/ip/diffim/subtractImages.py @@ -207,7 +207,8 @@ class AlardLuptonSubtractBaseConfig(lsst.pex.config.Config): doScaleVariance = lsst.pex.config.Field( dtype=bool, default=True, - doc="Scale variance of the image difference?" + doc="Scale variance of the science image? Note that the template variance is NOT scaled" + " here. The template variance may be scaled independently in ``GetTemplateTask``." ) scaleVariance = lsst.pex.config.ConfigurableField( target=ScaleVarianceTask, @@ -1044,13 +1045,11 @@ def _prepareInputs(self, template, science, visitSummary=None): self.metadata["templateCoveragePercent"] = 100*templateCoverageFraction if self.config.doScaleVariance: - # Scale the variance of the template and science images before - # convolution, subtraction, or decorrelation so that they have the - # correct ratio. - templateVarFactor = self.scaleVariance.run(template.maskedImage) + # Scale the variance of the science image before + # convolution, subtraction, or decorrelation so that it has the + # correct ratio. Note that the template variance is scaled + # independently in ``GetTemplateTask``. sciVarFactor = self.scaleVariance.run(science.maskedImage) - self.log.info("Template variance scaling factor: %.2f", templateVarFactor) - self.metadata["scaleTemplateVarianceFactor"] = templateVarFactor self.log.info("Science variance scaling factor: %.2f", sciVarFactor) self.metadata["scaleScienceVarianceFactor"] = sciVarFactor diff --git a/tests/test_getTemplate.py b/tests/test_getTemplate.py index bad896f5..efcb2df4 100644 --- a/tests/test_getTemplate.py +++ b/tests/test_getTemplate.py @@ -337,6 +337,94 @@ def testNanInputs(self, box=None, nInput=None): # in the template are closer to the original anymore. self.assertTrue(np.isfinite(result.template.image.array).all()) + def _scaleInputVariance(self, tract, factor): + """Return fresh handles for one tract's patches, with their variance + planes multiplied by ``factor``. + + Parameters + ---------- + tract : `int` + Id of the tract whose patches should be copied. + factor : `float` + Factor to multiply the input variance planes by. + + Returns + ------- + handles : `list` [`lsst.pipe.base.InMemoryDatasetHandle`] + Handles to the modified patches. + """ + handles = [] + for handle in self.patches[tract]: + # ``copy=True`` on the original handles means this is a copy, so + # the patches shared with the other tests are left untouched. + patch = handle.get() + patch.variance.array *= factor + handles.append(pipeBase.InMemoryDatasetHandle(patch, + storageClass="ExposureF", + copy=True, + dataId=handle.dataId)) + return handles + + def testScaleVariance(self): + """Test that the template variance plane is rescaled to match the + empirical pixel noise, and that the factor used is recorded in the + task metadata. + """ + scaleFactor = 1.345 + box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180)) + + def _configureAndRunTask(doScaleVariance, varianceScale=1.): + """Build a template from tract 0, optionally rescaling the input + variance planes by ``varianceScale`` first. + """ + config = lsst.ip.diffim.GetTemplateTask.ConfigClass() + config.doScaleVariance = doScaleVariance + task = lsst.ip.diffim.GetTemplateTask(config=config) + # Task modifies the input bbox, so pass a copy. + result = task.run(coaddExposureHandles={0: self._scaleInputVariance(0, varianceScale)}, + bbox=lsst.geom.Box2I(box), + wcs=self.exposure.wcs, + dataIds={0: self.dataIds[0]}, + physical_filter="a_test") + return task, result.template + + # With scaling disabled the subtask is never constructed, and nothing + # is recorded in the metadata. + taskOff, templateOff = _configureAndRunTask(False) + self.assertFalse(hasattr(taskOff, "scaleVariance")) + self.assertNotIn("scaleTemplateVarianceFactor", taskOff.metadata) + + # Both warps -- lanczos5 in ``_makePatches`` and lanczos3 in the + # task -- correlate the noise. The variance plane tracks only the + # per-pixel diagonal, which the second warp leaves too low, so + # ``scaleVariance`` measures a factor well above 1 even though the + # input variance planes are correct. + # + taskOn, templateOn = _configureAndRunTask(True) + factor = taskOn.metadata["scaleTemplateVarianceFactor"] + # TODO DM-55879: this value is pinned on purpose. The lanczos warping + # kernels introduce small correlations that artificially suppress the + # image pixel stddev and inflate the variance scaling factor. This + # should be changed to 1.0 after DM-55879 is merged. + self.assertFloatsAlmostEqual(factor, 1.1465, atol=0.01, + msg="Measured template variance scaling changed; see the" + " comment above if the correlation correction landed.") + # The only difference from the unscaled template is the constant + # factor applied to the variance plane. + self.assertFloatsAlmostEqual(templateOn.variance.array, + templateOff.variance.array*factor, rtol=1e-5) + # Tolerance here is float32 round-off: repeated runs of the task are + # not bitwise identical. + self.assertImagesAlmostEqual(templateOn.image, templateOff.image, rtol=1e-5, atol=1e-5) + + # If the input variance planes under-estimate the noise by a known + # factor, the measured factor grows by that amount and the same + # output variance plane is recovered. + taskLow, templateLow = _configureAndRunTask(True, varianceScale=1/scaleFactor) + self.assertFloatsAlmostEqual(taskLow.metadata["scaleTemplateVarianceFactor"], + factor*scaleFactor, rtol=1e-5) + self.assertImagesAlmostEqual(templateLow.variance, templateOn.variance, rtol=1e-5) + def setup_module(module): lsst.utils.tests.init() diff --git a/tests/test_subtractTask.py b/tests/test_subtractTask.py index fa30f736..3847aa86 100644 --- a/tests/test_subtractTask.py +++ b/tests/test_subtractTask.py @@ -349,7 +349,6 @@ def _run_and_check_images(statsCtrl, statsCtrlDetect, scienceNoiseLevel, templat templateBorderSize=20, doApplyCalibration=True) task = self._setup_subtraction(mode="convolveScience") output = task.run(template, science, sources) - self.assertFloatsAlmostEqual(task.metadata["scaleTemplateVarianceFactor"], 1., atol=.05) self.assertFloatsAlmostEqual(task.metadata["scaleScienceVarianceFactor"], 1., atol=.05) # Mean of difference image should be close to zero. nGoodPix = np.sum(np.isfinite(output.difference.image.array)) @@ -385,7 +384,6 @@ def _run_and_check_images(statsCtrl, statsCtrlDetect, scienceNoiseLevel, templat templateBorderSize=20, doApplyCalibration=True) task = self._setup_subtraction() output = task.run(template, science, sources) - self.assertFloatsAlmostEqual(task.metadata["scaleTemplateVarianceFactor"], 1., atol=.05) self.assertFloatsAlmostEqual(task.metadata["scaleScienceVarianceFactor"], 1., atol=.05) # There should be no NaNs in the image if we convolve the template with a buffer self.assertTrue(np.all(np.isfinite(output.difference.image.array))) @@ -632,8 +630,6 @@ def _run_and_check_images(science, template, sources, statsCtrl, ) output = task.run(template.clone(), science.clone(), sources) if doScaleVariance: - self.assertFloatsAlmostEqual(task.metadata["scaleTemplateVarianceFactor"], - scaleFactor, atol=0.05) self.assertFloatsAlmostEqual(task.metadata["scaleScienceVarianceFactor"], scaleFactor, atol=0.05) @@ -646,7 +642,8 @@ def _run_and_check_images(science, template, sources, statsCtrl, statsCtrl) if doScaleVariance: - templateNoise *= scaleFactor + # Only the science variance is scaled here. The template + # variance is scaled independently in ``GetTemplateTask``. scienceNoise *= scaleFactor varMean = computeRobustStatistics(output.difference.variance, output.difference.mask, statsCtrl) self.assertFloatsAlmostEqual(varMean, scienceNoise + templateNoise, rtol=0.1) @@ -666,8 +663,7 @@ def _run_and_check_images(science, template, sources, statsCtrl, doDecorrelation=False, doScaleVariance=False) # Verify that the variance plane of the difference image is correct - # when the template variance plane is incorrect - template.variance.array /= scaleFactor + # when the input science variance plane is incorrect science.variance.array /= scaleFactor _run_and_check_images(science, template, sources, statsCtrl, doDecorrelation=True, doScaleVariance=True, scaleFactor=scaleFactor) @@ -700,8 +696,6 @@ def _run_and_check_images(science, template, sources, statsCtrl, ) output = task.run(template.clone(), science.clone(), sources) if doScaleVariance: - self.assertFloatsAlmostEqual(task.metadata["scaleTemplateVarianceFactor"], - scaleFactor, atol=0.05) self.assertFloatsAlmostEqual(task.metadata["scaleScienceVarianceFactor"], scaleFactor, atol=0.05) @@ -714,7 +708,8 @@ def _run_and_check_images(science, template, sources, statsCtrl, statsCtrl) if doScaleVariance: - templateNoise *= scaleFactor + # Only the science variance is scaled here. The template + # variance is scaled independently in ``GetTemplateTask``. scienceNoise *= scaleFactor varMean = computeRobustStatistics(output.difference.variance, output.difference.mask, statsCtrl) @@ -735,9 +730,8 @@ def _run_and_check_images(science, template, sources, statsCtrl, doDecorrelation=False, doScaleVariance=False) # Verify that the variance plane of the difference image is correct - # when the template and science variance planes are incorrect + # when the input science variance plane is incorrect science.variance.array /= scaleFactor - template.variance.array /= scaleFactor _run_and_check_images(science, template, sources, statsCtrl, doDecorrelation=True, doScaleVariance=True, scaleFactor=scaleFactor) _run_and_check_images(science, template, sources, statsCtrl, @@ -1259,8 +1253,6 @@ def _run_and_check_images(science, template, sources, statsCtrl, ) output = task.run(template.clone(), science.clone(), sources) if doScaleVariance: - self.assertFloatsAlmostEqual(task.metadata["scaleTemplateVarianceFactor"], - scaleFactor, atol=0.05) self.assertFloatsAlmostEqual(task.metadata["scaleScienceVarianceFactor"], scaleFactor, atol=0.05) @@ -1278,7 +1270,8 @@ def _run_and_check_images(science, template, sources, statsCtrl, output.matchedTemplate.mask, statsCtrl) if doScaleVariance: - templateNoise *= scaleFactor + # Only the science variance is scaled here. The template + # variance is scaled independently in ``GetTemplateTask``. scienceNoise *= scaleFactor varMean = computeRobustStatistics(output.scoreExposure.variance, output.scoreExposure.mask, @@ -1302,8 +1295,7 @@ def _run_and_check_images(science, template, sources, statsCtrl, doDecorrelation=False, doScaleVariance=False) # Verify that the variance plane of the Score image is correct - # when the template variance plane is incorrect - template.variance.array /= scaleFactor + # when the input science variance plane is incorrect science.variance.array /= scaleFactor _run_and_check_images(science, template, sources, statsCtrl, doDecorrelation=True, doScaleVariance=True, scaleFactor=scaleFactor)