From 882235914c4f98e096d850f18d6557587904dec9 Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Wed, 21 Oct 2026 15:04:58 -0700 Subject: [PATCH 1/4] Correct the template variance for warping kernel correlation Warping with a kernel normalized to unit sum moves noise power off the variance plane's diagonal and into covariance between neighboring pixels, which a variance plane cannot represent, leaving the variance too small for anything measured over more than one pixel. --- python/lsst/ip/diffim/getTemplate.py | 111 ++++++++++++++++++++++++--- 1 file changed, 99 insertions(+), 12 deletions(-) diff --git a/python/lsst/ip/diffim/getTemplate.py b/python/lsst/ip/diffim/getTemplate.py index 0bf168c4..aa9920df 100644 --- a/python/lsst/ip/diffim/getTemplate.py +++ b/python/lsst/ip/diffim/getTemplate.py @@ -20,6 +20,7 @@ # along with this program. If not, see . import collections import dataclasses +import functools import astropy.units from astropy.stats import gaussian_sigma_to_fwhm @@ -65,6 +66,47 @@ """ +@functools.cache +def _warpKernelSumSq(kernelName, nPhase=100): + """Compute the phase-averaged sum of squared values of a warping kernel. + + Warping with a kernel normalized to unit sum multiplies the variance by + this sum, which is less than one. The power the warp took off the + diagonal became covariance between neighboring pixels, which a variance + plane cannot represent, so dividing by it restores the noise power: the + variance of a flux measured over a region and not a single pixel. + + Parameters + ---------- + kernelName : `str` + Name of the warping kernel, as `lsst.afw.math.Warper` takes it. + nPhase : `int`, optional + Number of sub-pixel phases to average over. + + Returns + ------- + sumSq : `float` + Sum of the squared 2-D kernel values, averaged over sub-pixel + phase: 0.798 for ``lanczos3``, 4/9 for ``bilinear``, and 1 for no + interpolation. + + Notes + ----- + The kernels are separable and the phases of the two axes independent, + so this is the square of the 1-D average. + """ + # Disable the cache since we are evaluating fewer phase points than the + # default cache size. + kernel = afwMath.Warper(kernelName, cacheSize=0).getWarpingKernel() + image = afwImage.ImageD(kernel.getDimensions()) + total = 0.0 + for phase in (np.arange(nPhase) + 0.5)/nPhase: + kernel.setKernelParameters((float(phase), 0.0)) + kernel.computeImage(image, True) + total += float((image.array.sum(axis=0)**2).sum()) + return (total/nPhase)**2 + + class GetTemplateConnections( pipeBase.PipelineTaskConnections, dimensions=("instrument", "visit", "detector"), @@ -149,6 +191,29 @@ class GetTemplateConfig( dtype=afwMath.Warper.ConfigClass, doc="warper configuration", ) + doCorrectVarianceWarpCorrelation = pexConfig.Field( + dtype=bool, + default=True, + doc="Correct the template variance for the correlation that warping" + " introduced between neighboring pixels? This covers both warps: the" + " one that built the coadds, named by coaddWarpKernel, and the one" + " this task performs, named by warp.warpingKernelName. It is a single" + " scalar factor, fixed entirely by those two kernels and independent" + " of the pixel scales involved.", + ) + coaddWarpKernel = pexConfig.ChoiceField( + dtype=str, + default="lanczos3", + allowed={ + "bilinear": "bilinear interpolation", + "lanczos3": "Lanczos kernel of order 3", + "lanczos4": "Lanczos kernel of order 4", + "lanczos5": "Lanczos kernel of order 5", + }, + doc="Warping kernel that was used to build the input coadds. This" + " cannot be recovered from the coadds themselves, and is used only if" + " doCorrectVarianceWarpCorrelation is True.", + ) coaddPsf = pexConfig.ConfigField( doc="Configuration for CoaddPsf", dtype=CoaddPsfConfig, @@ -650,14 +715,20 @@ def run(self, *, coaddExposureHandles, bbox, wcs, dataIds, physical_filter, visi def _correctVariance(self, template, plateScaleFactor): """Correct the template variance plane for the effects of warping. - Apply a multiplicative correction factor to the template variance if - the plate scale of the coadd's constituent images is different than - the science image the template is being constructed for. This should - only be necessary if the coadd images were from a different instrument - than the science image. If the plate scale of the science instrument is - smaller than the plate scale of the coadd instrument, then the template - pixels will be correlated and the true variance will be higher than the - image pixel noise level would suggest. + Warping introduces a covariance between neighboring pixels, while the + input variance plane captures only the diagonal term. Compute a + multiplicative correction factor to apply to the final variance so that + the detection significance on the convolved difference image is correct + on average. + + Separately, also apply a scalar multiplicative correction factor to the + template variance if the plate scale of the coadd's constituent images + is different than the science image the template is being constructed + for. This should only be necessary if the coadd images were from a + different instrument than the science image. If the plate scale of the + science instrument is smaller than the plate scale of the coadd + instrument, then the template pixels will be correlated and the true + variance will be higher than the image pixel noise level would suggest. Parameters ---------- @@ -677,14 +748,30 @@ def _correctVariance(self, template, plateScaleFactor): " so the template variance cannot be corrected for the change in pixel area." ) self.metadata["variancePlateScaleFactor"] = plateScaleFactor - if plateScaleFactor == 1.0: + correlationFactor = 1.0 + if self.config.doCorrectVarianceWarpCorrelation: + # Both warps correlate the noise: the one that built the coadds, + # and the one this task just performed. + correlationFactor = 1/(_warpKernelSumSq(self.config.coaddWarpKernel) + * _warpKernelSumSq(self.config.warp.warpingKernelName)) + self.metadata["varianceWarpCorrelationFactor"] = correlationFactor + scale = correlationFactor*plateScaleFactor + if scale == 1.0: return - template.variance.array *= plateScaleFactor - self.metadata["templateVarianceCorrectionFactor"] = plateScaleFactor + template.variance.array *= scale + self.metadata["templateVarianceCorrectionFactor"] = scale self.log.info( - "Corrected the template variance plane for the change in pixel area by a factor of %.4f.", + "Variance correction factor due to plate scale changes: %.4f", plateScaleFactor, ) + self.log.info( + "Variance correction factor due to warping correlations: %.4f", + correlationFactor, + ) + self.log.info( + "Overall template variance correction factor: %.4f", + scale, + ) def _plateScaleFactor(self, coaddInputs, wcs, bbox): """Compute the variance correction for the total change in pixel area From 7e7eacf628fa854125ed99ac0331342f02711c8c Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Wed, 21 Oct 2026 17:33:00 -0700 Subject: [PATCH 2/4] Turn off doScaleVariance of the template by default The correction factors should now result in a better calculation of the variance, and doScaleVariance will actually make the result worse. --- python/lsst/ip/diffim/getTemplate.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/python/lsst/ip/diffim/getTemplate.py b/python/lsst/ip/diffim/getTemplate.py index aa9920df..98219892 100644 --- a/python/lsst/ip/diffim/getTemplate.py +++ b/python/lsst/ip/diffim/getTemplate.py @@ -248,8 +248,10 @@ class GetTemplateConfig( ) doScaleVariance = pexConfig.Field( dtype=bool, - default=True, - doc="Scale variance of the template image?" + default=False, + doc="Scale the variance of the template image to match its per-pixel" + " scatter? Warping correlates the noise, so when combined with" + " doCorrectVarianceWarpCorrelation this overcounts the correction.", ) scaleVariance = pexConfig.ConfigurableField( target=ScaleVarianceTask, From 69fc8b256586ac441da121439c7ff18da25b00eb Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Wed, 21 Oct 2026 17:38:21 -0700 Subject: [PATCH 3/4] Test the template warping correlation correction --- tests/test_getTemplate.py | 259 +++++++++++++++++++++++++++++++------- 1 file changed, 216 insertions(+), 43 deletions(-) diff --git a/tests/test_getTemplate.py b/tests/test_getTemplate.py index 9ae6bee0..fea1af70 100644 --- a/tests/test_getTemplate.py +++ b/tests/test_getTemplate.py @@ -58,6 +58,15 @@ display = lsst.afw.display.Display() display.frame = 1 +# Sum of squared warping kernel values along one axis, with the kernel +# normalized to unit sum at each sub-pixel phase, averaged over phase. +KERNEL_SUM_SQ_1D = { + "bilinear": 2/3, + "lanczos3": 0.8933699084, + "lanczos4": 0.9141249535, + "lanczos5": 0.9342245513, +} + def _showTemplate(box, template): """Show the corners of the template we made in this test.""" @@ -80,6 +89,8 @@ class GetTemplateTaskTestCase(lsst.utils.tests.TestCase): """ def setUp(self): self.scale = 0.2 # arcsec/pixel + # Use 5th order to minimize artifacts on the templates. + self.coaddWarpKernel = "lanczos5" self.skymap = self._makeSkymap() self.patches = collections.defaultdict(list) self.dataIds = collections.defaultdict(list) @@ -151,9 +162,8 @@ def _makePatches(self, tract): # Show the tract center, colored by tract id. display.dot("x", point.x, point.y, ctype=color, size=30) - # Use 5th order to minimize artifacts on the templates. config = lsst.afw.math.Warper.ConfigClass() - config.warpingKernelName = "lanczos5" + config.warpingKernelName = self.coaddWarpKernel warper = lsst.afw.math.Warper.fromConfig(config) for patchId in range(tract.num_patches.x*tract.num_patches.y): patch = tract.getPatchInfo(patchId) @@ -218,6 +228,14 @@ def _checkMetadata(self, template, config, box, wcs, nPsfs): self.assertTrue(template.getInfo().hasCoaddInputs()) self.assertEqual(len(template.getInfo().getCoaddInputs().ccds), nPsfs) + def _makeConfig(self): + """Make a GetTemplateTask config whose coaddWarpKernel is the kernel + the test patches were made with. + """ + config = lsst.ip.diffim.GetTemplateTask.ConfigClass() + config.coaddWarpKernel = self.coaddWarpKernel + return config + def _makeCoaddPsf(self, template, config): """Return the CoaddPsf that ``run`` approximated with a Gaussian, rebuilt from the coadd inputs recorded on the template. @@ -246,15 +264,13 @@ def _checkPixels(self, template, config, box): self.assertImagesAlmostEqual(template.image, expected.image, rtol=.1, atol=4*calibration) # Variance plane ==4 in the original image (realize() takes a noise # sigma). Warping sets the level from the pixel areas and the two - # warping kernels, which `_correctVariance` corrects up to the - # difference between the configured coaddWarpKernel and the lanczos5 - # `_makePatches` really used. A per-pixel ripple that no scalar - # correction can remove remains on top of that, so check the level - # tightly and allow for the ripple around it. + # warping kernels, which `_correctVariance` corrects. A per-pixel + # ripple that no scalar correction can remove remains on top of that, + # so check the level tightly and allow for the ripple around it. variance = template.variance.array[np.isfinite(template.variance.array)] median = np.median(variance) self.assertFloatsAlmostEqual(median, np.median(expected.variance.array), - rtol=0.35, msg="variance level differs") + rtol=0.02, msg="variance level differs") self.assertLess(np.percentile(variance, 99)/median, 1.6, msg="variance ripple too large") self.assertGreater(np.percentile(variance, 1)/median, 0.6, msg="variance ripple too large") # Not checking the mask, as warping changes the sizes of the masks. @@ -265,7 +281,7 @@ def testRunOneTractInput(self): tract input correctly. """ box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180)) - task = lsst.ip.diffim.GetTemplateTask() + task = lsst.ip.diffim.GetTemplateTask(config=self._makeConfig()) # Restrict to tract 0, since the box fits in just that tract. result = task.run(coaddExposureHandles={0: self.patches[0]}, bbox=box, @@ -283,7 +299,7 @@ def testRunOneTractMultipleInputs(self): mostly NaN after warping are merged correctly in the output. """ box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180)) - task = lsst.ip.diffim.GetTemplateTask() + task = lsst.ip.diffim.GetTemplateTask(config=self._makeConfig()) result = task.run(coaddExposureHandles=self.patches, bbox=box, wcs=self.exposure.wcs, @@ -298,7 +314,7 @@ def testRunTwoTracts(self): """Test a bounding box that crosses tract boundaries. """ box = lsst.geom.Box2I(lsst.geom.Point2I(200, 200), lsst.geom.Point2I(600, 600)) - task = lsst.ip.diffim.GetTemplateTask() + task = lsst.ip.diffim.GetTemplateTask(config=self._makeConfig()) result = task.run(coaddExposureHandles=self.patches, bbox=box, wcs=self.exposure.wcs, @@ -316,7 +332,7 @@ def testRunCalibratesToNanojansky(self): box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180)) calibration = self.exposure.photoCalib.getCalibrationMean() self.assertNotEqual(calibration, 1.0) - task = lsst.ip.diffim.GetTemplateTask() + task = lsst.ip.diffim.GetTemplateTask(config=self._makeConfig()) templates = {} for photoCalib in (self.exposure.photoCalib, lsst.afw.image.PhotoCalib(1.0)): self._remakePatches(photoCalib) @@ -348,7 +364,7 @@ def testGaussianPsfGridFallback(self): Gaussian width is averaged over a grid of positions instead. """ box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180)) - task = lsst.ip.diffim.GetTemplateTask() + task = lsst.ip.diffim.GetTemplateTask(config=self._makeConfig()) result = task.run(coaddExposureHandles={0: self.patches[0]}, bbox=box, wcs=self.exposure.wcs, @@ -376,7 +392,7 @@ def testRunNoTemplate(self): """A bounding box that doesn't overlap the patches will raise. """ box = lsst.geom.Box2I(lsst.geom.Point2I(1200, 1200), lsst.geom.Point2I(1600, 1600)) - task = lsst.ip.diffim.GetTemplateTask() + task = lsst.ip.diffim.GetTemplateTask(config=self._makeConfig()) with self.assertRaisesRegex(lsst.pipe.base.NoWorkFound, "No patches found"): task.run(coaddExposureHandles=self.patches, bbox=box, @@ -393,7 +409,7 @@ def testMissingPatches(self): # tract=0, patch=1 is the lower-left corner, as displayed in DS9. self.patches[0].pop(1) box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180)) - task = lsst.ip.diffim.GetTemplateTask() + task = lsst.ip.diffim.GetTemplateTask(config=self._makeConfig()) result = task.run(coaddExposureHandles=self.patches, bbox=box, wcs=self.exposure.wcs, @@ -424,7 +440,7 @@ def testNanInputs(self, box=None, nInput=None): patchCoadd.variance[bbox].array *= np.nan box = lsst.geom.Box2I(lsst.geom.Point2I(200, 200), lsst.geom.Point2I(600, 600)) - task = lsst.ip.diffim.GetTemplateTask() + task = lsst.ip.diffim.GetTemplateTask(config=self._makeConfig()) result = task.run(coaddExposureHandles=self.patches, bbox=box, wcs=self.exposure.wcs, @@ -437,16 +453,22 @@ 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 _runCorrection(self, handles=None, instrument="SciCam"): - """Run the task on tract 0, with a science image from ``instrument`` - and variance scaling turned off. + def _runCorrection(self, doCorrectVarianceWarpCorrelation=False, doScaleVariance=False, + handles=None, instrument="SciCam", coaddWarpKernel=None, + templateWarpKernel="lanczos3"): + """Run the task on tract 0, with a science image from ``instrument``. The plate scale correction is applied only when ``instrument`` is not the instrument of the template coadds, which is ``"TestCam"`` unless - ``handles`` says otherwise. + ``handles`` says otherwise. ``coaddWarpKernel`` defaults to the kernel + the patches were made with, so the reconstructed factors are exact. """ - config = lsst.ip.diffim.GetTemplateTask.ConfigClass() - config.doScaleVariance = False + config = self._makeConfig() + config.doCorrectVarianceWarpCorrelation = doCorrectVarianceWarpCorrelation + config.doScaleVariance = doScaleVariance + if coaddWarpKernel is not None: + config.coaddWarpKernel = coaddWarpKernel + config.warp.warpingKernelName = templateWarpKernel task = lsst.ip.diffim.GetTemplateTask(config=config) result = task.run(coaddExposureHandles={0: handles or self.patches[0]}, bbox=lsst.geom.Box2I(self.varianceBox), @@ -509,6 +531,38 @@ def _patchHandles(self, tract, records, instrument="TestCam"): coadd, storageClass="ExposureF", copy=True, dataId=ref.dataId)) return handles + def testCorrectVarianceWarpCorrelation(self): + """The correlation correction covers both warps: the one that built + the coadds and the one this task performs. + """ + _, off = self._runCorrection(instrument="TestCam") + for coaddKernel, templateKernel in (("lanczos3", "lanczos3"), + ("lanczos5", "lanczos3"), + ("lanczos5", "lanczos5"), + ("bilinear", "lanczos3")): + with self.subTest(coaddWarpKernel=coaddKernel, templateWarpKernel=templateKernel): + task, on = self._runCorrection(doCorrectVarianceWarpCorrelation=True, + instrument="TestCam", + coaddWarpKernel=coaddKernel, + templateWarpKernel=templateKernel) + factor = task.metadata["varianceWarpCorrelationFactor"] + # A kernel spreads power, so this always makes the plane + # larger. + self.assertGreater(factor, 1.0) + # The phase average converges slowly for bilinear, because + # its sum of squares has a kink where the phase wraps. + rtol = 1e-4 if "bilinear" in (coaddKernel, templateKernel) else 1e-5 + self.assertFloatsAlmostEqual( + factor, + 1/(KERNEL_SUM_SQ_1D[coaddKernel]*KERNEL_SUM_SQ_1D[templateKernel])**2, + rtol=rtol) + self.assertEqual(task.metadata["variancePlateScaleFactor"], 1.0) + if templateKernel == "lanczos3": + # A pure scalar, with no per-pixel structure of its own; + # only comparable when `off` used the same template warp. + self.assertFloatsAlmostEqual(on.variance.array, off.variance.array*factor, + rtol=1e-5, ignoreNaNs=True) + def testCorrectVariancePlateScale(self): """The plate scale correction is the total pixel area change from the images the coadds were built from to the science image. @@ -523,6 +577,7 @@ def testCorrectVariancePlateScale(self): # net change in pixel area even when the instruments differ. task, on = self._runCorrection() self.assertFloatsAlmostEqual(task.metadata["variancePlateScaleFactor"], 1.0, rtol=1e-6) + self.assertNotIn("varianceWarpCorrelationFactor", task.metadata) # Coarser original pixels than science pixels, the DECam-template # case: the correction is the ratio of their areas. @@ -536,6 +591,49 @@ def testCorrectVariancePlateScale(self): self.assertFloatsAlmostEqual(on.variance.array, off.variance.array*factor, rtol=1e-5, ignoreNaNs=True) + def testVarianceCorrectionsCompose(self): + """The two corrections are independent and multiply together. + + The patches are warped from the exposure and the template is warped + back to its pixel grid, so the two warps cancel geometrically, leaving + only the two kernels' correlation to undo. The substituted coadd + inputs claim coarser original pixels than the exposure really has, so + with both corrections applied the variance must come back to its + original value times the claimed change in pixel area. + """ + plateScale = 1.315 + handles = self._patchHandles(0, [(self._makeScaledWcs(plateScale), self.exposure.getBBox())]) + # A science image from the coadds' instrument applies no correction + # at all; one from another instrument applies both. + _, off = self._runCorrection(handles=handles, instrument="TestCam") + task, on = self._runCorrection(doCorrectVarianceWarpCorrelation=True, handles=handles) + + # Both together is exactly the product of the two factors. + self.assertFloatsAlmostEqual(task.metadata["variancePlateScaleFactor"], plateScale**2, rtol=1e-6) + scale = (task.metadata["varianceWarpCorrelationFactor"] + * task.metadata["variancePlateScaleFactor"]) + self.assertFloatsAlmostEqual(task.metadata["templateVarianceCorrectionFactor"], + scale, rtol=1e-12) + self.assertFloatsAlmostEqual(on.variance.array, off.variance.array*scale, + rtol=1e-5, ignoreNaNs=True) + + # Only the variance plane is touched. Nothing about the correction + # feeds back into how the patches are combined. + self.assertImagesEqual(on.image, off.image) + self.assertMasksEqual(on.mask, off.mask) + + # realize() takes a noise sigma, so the original variance is 4, not 2, + # before calibrating to nJy. + calibration = self.exposure.photoCalib.getCalibrationMean() + original = calibration**2*np.median(self.exposure.variance.array) + self.assertFloatsAlmostEqual(np.median(on.variance.array), original*plateScale**2, rtol=0.02) + + # Which is a real correction: uncorrected, the plane is too small by + # the two warping kernels' phase-averaged sums of squares. + self.assertFloatsAlmostEqual( + np.median(off.variance.array)/original, + (KERNEL_SUM_SQ_1D[self.coaddWarpKernel]*KERNEL_SUM_SQ_1D["lanczos3"])**2, rtol=0.02) + def testCorrectVariancePlateScaleUsesOneRecord(self): """Only the first usable coadd input is read. @@ -567,7 +665,7 @@ def testCorrectVariancePlateScaleNextTract(self): """If the first coadd kept in a tract cannot give the plate scale, the first coadd kept in the next tract is used. """ - task = lsst.ip.diffim.GetTemplateTask() + task = lsst.ip.diffim.GetTemplateTask(config=self._makeConfig()) handles = {0: self._patchHandles(0, []), 1: self._patchHandles(1, [(self._makeScaledWcs(1.315), self.exposure.getBBox())])} task.run(coaddExposureHandles=handles, @@ -590,7 +688,7 @@ def _runTract0(self, handles): """Build a template from the given tract 0 patch handles. """ box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180)) - task = lsst.ip.diffim.GetTemplateTask() + task = lsst.ip.diffim.GetTemplateTask(config=self._makeConfig()) return task.run(coaddExposureHandles={0: handles}, bbox=box, wcs=self.exposure.wcs, @@ -654,7 +752,7 @@ def _runLogged(self, handles, dataIds=None, box=None, instrument=None): dataIds = {tract: self.dataIds[tract] for tract in handles} if box is None: box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180)) - task = lsst.ip.diffim.GetTemplateTask() + task = lsst.ip.diffim.GetTemplateTask(config=self._makeConfig()) with self.assertLogs(task.log.name, level="INFO") as cm: result = task.run(coaddExposureHandles=handles, bbox=box, @@ -770,7 +868,7 @@ 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 = self._makeConfig() config.doScaleVariance = doScaleVariance task = lsst.ip.diffim.GetTemplateTask(config=config) result = task.run(coaddExposureHandles={0: self._scaleInputVariance(0, varianceScale)}, @@ -782,25 +880,19 @@ def _configureAndRunTask(doScaleVariance, varianceScale=1.): # With scaling disabled the subtask is never constructed, and nothing # is recorded in the metadata. - taskOff, templateOff = _configureAndRunTask(False) + taskOff, templateOff = _configureAndRunTask(doScaleVariance=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) + # The image is correlated after warping, so ``doScaleVariance`` will + # apply a correction. The value here is empirical from a warped image + # with white noise and no sources. + taskOn, templateOn = _configureAndRunTask(doScaleVariance=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.") + self.assertFloatsAlmostEqual(factor, 1.145, atol=0.01, + msg="Measured template variance scaling changed; check the" + " warping kernels used by _makePatches and the task," + " and the order of the variance corrections.") # The only difference from the unscaled template is the constant # factor applied to the variance plane. self.assertFloatsAlmostEqual(templateOn.variance.array, @@ -812,7 +904,7 @@ def _configureAndRunTask(doScaleVariance, varianceScale=1.): # 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) + taskLow, templateLow = _configureAndRunTask(doScaleVariance=True, varianceScale=1/scaleFactor) self.assertFloatsAlmostEqual(taskLow.metadata["scaleTemplateVarianceFactor"], factor*scaleFactor, rtol=1e-5) self.assertImagesAlmostEqual(templateLow.variance, templateOn.variance, rtol=1e-5) @@ -838,7 +930,7 @@ def _runLegacyForFuture(self, raiseOnUndefinedMaskMap=True, handles=None): The bounding box the template was requested on, before the task grew it by the template border. """ - config = lsst.ip.diffim.GetTemplateTask.ConfigClass() + config = self._makeConfig() config.output_image_type = "future" config.raiseOnUndefinedMaskMap = raiseOnUndefinedMaskMap task = lsst.ip.diffim.GetTemplateTask(config=config) @@ -1082,6 +1174,87 @@ def testConvertOutputsToFuturePsfToLegacy(self): self.assertFloatsAlmostEqual(psf.getSigma(), result.template.psf.sigma) +class WarpCorrelationFactorTestCase(lsst.utils.tests.TestCase): + """Test the warping correlation correction on pure noise, with no other + variance correction. + """ + + def setUp(self): + self.size = 1024 + self.bbox = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Extent2I(self.size, self.size)) + + def _makeWcs(self, angle, offset): + """Make a WCS with a fixed plate scale, rotated by ``angle`` degrees + and with its reference pixel moved by ``offset`` pixels. + """ + center = lsst.geom.Point2D(self.size/2 + offset[0], self.size/2 + offset[1]) + cdMatrix = lsst.afw.geom.makeCdMatrix(0.2*lsst.geom.arcseconds, angle*lsst.geom.degrees) + return lsst.afw.geom.makeSkyWcs(center, lsst.geom.SpherePoint(45, 30, lsst.geom.degrees), + cdMatrix) + + def testCorrelationFactorRestoresNoisePower(self): + """Warp pure noise to a rotated, offset frame with the coadd kernel, + then to a third frame with the task's own warper. + + Warping leaves the variance of a flux summed over a region (the noise + power) unchanged, but lowers the per-pixel scatter and, by more, the + variance plane. The correction must bring the variance plane back to + the noise power of the original image. + """ + sigma = 2.0 + rng = np.random.default_rng(42) + original = lsst.afw.image.ExposureF(self.bbox, self._makeWcs(0, (0, 0))) + original.image.array[:] = rng.normal(0, sigma, original.image.array.shape) + original.variance.array[:] = sigma**2 + + config = lsst.ip.diffim.GetTemplateTask.ConfigClass() + config.doScaleVariance = False + config.coaddWarpKernel = "lanczos3" + task = lsst.ip.diffim.GetTemplateTask(config=config) + coadd = lsst.afw.math.Warper(config.coaddWarpKernel).warpExposure( + self._makeWcs(37, (0.31, 0.72)), original, destBBox=self.bbox) + template = task.warper.warpExposure(self._makeWcs(-21, (0.58, 0.23)), coadd, destBBox=self.bbox) + + # The central quarter is covered in every frame. + inner = lsst.geom.Box2I(lsst.geom.Point2I(self.size//4, self.size//4), + lsst.geom.Extent2I(self.size//2, self.size//2)) + originalStd = original.image[inner].array.std() + image = template.image[inner].array.astype(float) + pixelStd = image.std() + # Each warp moves flux between neighboring pixels but doesn't create + # or destroy any, so the noise in a flux summed over a region is the + # same before and after. Summing the image over blocks and measuring + # how much those sums scatter gives that quantity directly, without + # depending on the variance plane. Every placement of the block is + # used, from a summed-area table, to reduce the noise. + blockSize = 8 + table = np.pad(image.cumsum(axis=0).cumsum(axis=1), ((1, 0), (1, 0))) + sums = (table[blockSize:, blockSize:] - table[:-blockSize, blockSize:] + - table[blockSize:, :-blockSize] + table[:-blockSize, :-blockSize]) + blockStd = sums.std()/blockSize + planeStdBefore = np.sqrt(template.variance[inner].array.mean()) + task._correctVariance(template, plateScaleFactor=1.0) + planeStdAfter = np.sqrt(template.variance[inner].array.mean()) + factor = task.metadata["varianceWarpCorrelationFactor"] + + # afw propagates the variance as if each warp's input pixels were + # independent, so the plane drops by exactly the factor's reciprocal. + self.assertFloatsAlmostEqual(planeStdBefore, sigma/np.sqrt(factor), rtol=1e-5) + self.assertFloatsAlmostEqual(planeStdAfter, sigma, rtol=1e-5) + + # The warps conserve the noise power, which the corrected plane + # matches. + self.assertFloatsAlmostEqual(blockStd, originalStd, rtol=0.03) + self.assertFloatsAlmostEqual(blockStd, planeStdAfter, rtol=0.03) + + # The per-pixel scatter drops, though by less than the uncorrected + # plane, because the second warp receives correlated input. The + # corrected plane deliberately exceeds it. + self.assertLess(pixelStd, 0.9*originalStd) + self.assertGreater(pixelStd, planeStdBefore) + self.assertLess(pixelStd, planeStdAfter) + + class GetTemplateConnectionsTestCase(lsst.utils.tests.TestCase): """Test the connections that ``output_image_type`` switches on. From e1ca2e30146dfdf660867c9d2780513a05a23266 Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Wed, 21 Oct 2026 17:38:26 -0700 Subject: [PATCH 4/4] Build the getTemplate test patches once per class Warping all 16 patches took a couple seconds per test, and testRunCalibratesToNanojansky warped them twice more; it now relabels the patches' calibration instead. The getTemplate tests are about 6x faster with these changes. --- tests/test_getTemplate.py | 95 ++++++++++++++++++++++----------------- 1 file changed, 54 insertions(+), 41 deletions(-) diff --git a/tests/test_getTemplate.py b/tests/test_getTemplate.py index fea1af70..b073635d 100644 --- a/tests/test_getTemplate.py +++ b/tests/test_getTemplate.py @@ -87,25 +87,40 @@ class GetTemplateTaskTestCase(lsst.utils.tests.TestCase): and tests GetTemplateTask's output against boxes that overlap various combinations of one or multiple tracts. """ - def setUp(self): - self.scale = 0.2 # arcsec/pixel - # Use 5th order to minimize artifacts on the templates. - self.coaddWarpKernel = "lanczos5" - self.skymap = self._makeSkymap() - self.patches = collections.defaultdict(list) - self.dataIds = collections.defaultdict(list) - self.coaddRefs = {} - self.exposure = self._makeExposure() - self.varianceBox = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180)) + scale = 0.2 # arcsec/pixel + # Use 5th order to minimize artifacts on the templates. + coaddWarpKernel = "lanczos5" + + @classmethod + def setUpClass(cls): + # Warping the patches is the slowest part of every test, so it is done + # once. Tests must not modify the exposure; the patch handles return + # copies, and setUp gives each test its own lists of them. + cls.skymap = cls._makeSkymap() + cls.patches = collections.defaultdict(list) + cls.dataIds = collections.defaultdict(list) + cls.exposure = cls._makeExposure() if debug: - display.image(self.exposure, "base exposure") + display.image(cls.exposure, "base exposure") for tract_id in range(4): - tract = self.skymap.generateTract(tract_id) - self._makePatches(tract) + tract = cls.skymap.generateTract(tract_id) + cls._makePatches(tract) - def _makeSkymap(self): + @classmethod + def tearDownClass(cls): + del cls.skymap, cls.patches, cls.dataIds, cls.exposure + + def setUp(self): + cls = type(self) + self.patches = collections.defaultdict(list, {k: list(v) for k, v in cls.patches.items()}) + self.dataIds = collections.defaultdict(list, {k: list(v) for k, v in cls.dataIds.items()}) + self.coaddRefs = {} + self.varianceBox = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180)) + + @classmethod + def _makeSkymap(cls): """Make a Skymap with 4 tracts with 4 patches each. """ tractScale = 0.02 # degrees @@ -121,14 +136,15 @@ def _makeSkymap(self): # Half the tract center step size, to keep the tract overlap small. config.radiusList = [tractScale/2 for c in coords] config.projection = "TAN" - config.pixelScale = self.scale + config.pixelScale = cls.scale config.tractOverlap = 0.0005 config.tractBuilder = "legacy" config.tractBuilder["legacy"].patchInnerDimensions = (300, 300) config.tractBuilder["legacy"].patchBorder = 10 return lsst.skymap.DiscreteSkyMap(config=config) - def _makeExposure(self): + @classmethod + def _makeExposure(cls): """Create a large image to break up into tracts and patches. The image will have a source every 100 pixels in x and y, and a WCS @@ -141,7 +157,7 @@ def _makeExposure(self): # quadrant of the image, and the other tracts fill the rest; slight # extra rotation as a check on the final warp layout, scaled by 5% # from the patch pixel scale. - cd_matrix = lsst.afw.geom.makeCdMatrix(1.05*self.scale*lsst.geom.arcseconds, 93*lsst.geom.degrees) + cd_matrix = lsst.afw.geom.makeCdMatrix(1.05*cls.scale*lsst.geom.arcseconds, 93*lsst.geom.degrees) wcs = lsst.afw.geom.makeSkyWcs(lsst.geom.Point2D(120, 150), lsst.geom.SpherePoint(0, 0, lsst.geom.radians), cd_matrix) @@ -152,18 +168,19 @@ def _makeExposure(self): exposure.setFilter(lsst.afw.image.FilterLabel("a", "a_test")) return exposure - def _makePatches(self, tract): + @classmethod + def _makePatches(cls, tract): """Populate the patches and dataId dicts, keyed on tract id, with the warps of the main exposure and minimal dataIds, respectively. """ if debug: color = ['red', 'green', 'cyan', 'yellow'][tract.tract_id] - point = self.exposure.wcs.skyToPixel(tract.ctr_coord) + point = cls.exposure.wcs.skyToPixel(tract.ctr_coord) # Show the tract center, colored by tract id. display.dot("x", point.x, point.y, ctype=color, size=30) config = lsst.afw.math.Warper.ConfigClass() - config.warpingKernelName = self.coaddWarpKernel + config.warpingKernelName = cls.coaddWarpKernel warper = lsst.afw.math.Warper.fromConfig(config) for patchId in range(tract.num_patches.x*tract.num_patches.y): patch = tract.getPatchInfo(patchId) @@ -171,20 +188,20 @@ def _makePatches(self, tract): if debug: # Show the patch corners as patch ids, colored by tract id. - points = self.exposure.wcs.skyToPixel(patch.wcs.pixelToSky([lsst.geom.Point2D(x) + points = cls.exposure.wcs.skyToPixel(patch.wcs.pixelToSky([lsst.geom.Point2D(x) for x in box.getCorners()])) for p in points: display.dot(patchId, p.x, p.y, ctype=color) # This is mostly taken from drp_tasks makePsfMatchedWarp, but # ip_diffim cannot depend on drp_tasks. - xyTransform = lsst.afw.geom.makeWcsPairTransform(self.exposure.wcs, patch.wcs) - warpedPsf = lsst.meas.algorithms.WarpedPsf(self.exposure.psf, xyTransform) - warped = warper.warpExposure(patch.wcs, self.exposure, destBBox=box) + xyTransform = lsst.afw.geom.makeWcsPairTransform(cls.exposure.wcs, patch.wcs) + warpedPsf = lsst.meas.algorithms.WarpedPsf(cls.exposure.psf, xyTransform) + warped = warper.warpExposure(patch.wcs, cls.exposure, destBBox=box) warped.setPsf(warpedPsf) warped.getInfo().setCoaddInputs( - self._makeCoaddInputs([(self.exposure.wcs, self.exposure.getBBox())])) + cls._makeCoaddInputs([(cls.exposure.wcs, cls.exposure.getBBox())])) dataRef = pipeBase.InMemoryDatasetHandle( warped, storageClass="ExposureF", @@ -194,13 +211,13 @@ def _makePatches(self, tract): patch=patch, ) ) - self.patches[tract.tract_id].append(dataRef) + cls.patches[tract.tract_id].append(dataRef) dataCoordinate = DataCoordinate.standardize({"tract": tract.tract_id, "patch": patchId, "band": "a", "skymap": "skymap"}, universe=DimensionUniverse()) - self.dataIds[tract.tract_id].append(dataCoordinate) + cls.dataIds[tract.tract_id].append(dataCoordinate) def _checkMetadata(self, template, config, box, wcs, nPsfs): """Check that the various metadata components were set correctly. @@ -332,33 +349,29 @@ def testRunCalibratesToNanojansky(self): box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180)) calibration = self.exposure.photoCalib.getCalibrationMean() self.assertNotEqual(calibration, 1.0) + # The same coadd pixels, labelled as already in nJy. + uncalibrated = [] + for handle in self.patches[0]: + coadd = handle.get() + coadd.setPhotoCalib(lsst.afw.image.PhotoCalib(1.0)) + uncalibrated.append(pipeBase.InMemoryDatasetHandle(coadd, storageClass="ExposureF", copy=True, + dataId=handle.dataId)) task = lsst.ip.diffim.GetTemplateTask(config=self._makeConfig()) templates = {} - for photoCalib in (self.exposure.photoCalib, lsst.afw.image.PhotoCalib(1.0)): - self._remakePatches(photoCalib) - result = task.run(coaddExposureHandles={0: self.patches[0]}, + for key, handles in ((calibration, self.patches[0]), (1.0, uncalibrated)): + result = task.run(coaddExposureHandles={0: handles}, bbox=box, wcs=self.exposure.wcs, dataIds={0: self.dataIds[0]}, physical_filter="a_test") self.assertEqual(result.template.photoCalib, lsst.afw.image.PhotoCalib(1.0)) - templates[photoCalib.getCalibrationMean()] = result.template + templates[key] = result.template self.assertFloatsAlmostEqual(templates[calibration].image.array, calibration*templates[1.0].image.array, rtol=1e-6) self.assertFloatsAlmostEqual(templates[calibration].variance.array, calibration**2*templates[1.0].variance.array, rtol=1e-6) - def _remakePatches(self, photoCalib): - """Rebuild the coadd patches from the main exposure with a different - photometric calibration. - """ - self.exposure.setPhotoCalib(photoCalib) - self.patches.clear() - self.dataIds.clear() - for tract_id in range(4): - self._makePatches(self.skymap.generateTract(tract_id)) - def testGaussianPsfGridFallback(self): """If the CoaddPsf cannot be evaluated at its average position, the Gaussian width is averaged over a grid of positions instead.