Coverage for tests / test_subtractTask.py: 7%
795 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-05-07 08:35 +0000
1# This file is part of ip_diffim.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
22import unittest
24from astropy import units as u
26import lsst.afw.math as afwMath
27import lsst.afw.table as afwTable
28import lsst.geom
29import lsst.meas.algorithms as measAlg
30from lsst.ip.diffim import subtractImages, InsufficientKernelSourcesError
31from lsst.pex.config import FieldValidationError
32from lsst.pipe.base import NoWorkFound
33import lsst.utils.tests
34import numpy as np
35from lsst.ip.diffim.utils import (computeRobustStatistics, computePSFNoiseEquivalentArea,
36 evaluateMeanPsfFwhm, getPsfFwhm)
37from lsst.pex.exceptions import InvalidParameterError
39from utils import makeStats, makeTestImage, CustomCoaddPsf
42class AlardLuptonSubtractTestBase:
43 goodPsfSize = 2.0
44 midPsfSize = 2.4
45 badPsfSize = 2.8
47 def _setup_subtraction(self, fluxField="truth_instFlux", errField="truth_instFluxErr", **kwargs):
48 """Setup and configure the image subtraction PipelineTask.
50 Parameters
51 ----------
52 fluxField : `str`, optional
53 Name of the flux field in the source catalog.
54 errField : `str`, optional
55 Name of the flux error field in the source catalog.
56 **kwargs
57 Any additional config parameters to set.
59 Returns
60 -------
61 `lsst.pipe.base.PipelineTask`
62 The configured Task to use for detection and measurement.
63 """
64 config = self.subtractTask.ConfigClass()
65 config.doSubtractBackground = False
66 config.restrictKernelEdgeSources = False
67 config.sourceSelector.signalToNoise.fluxField = fluxField
68 config.sourceSelector.signalToNoise.errField = errField
69 config.sourceSelector.doUnresolved = True
70 config.sourceSelector.doIsolated = True
71 config.sourceSelector.doRequirePrimary = True
72 config.sourceSelector.doFlags = True
73 config.sourceSelector.doSignalToNoise = True
74 config.sourceSelector.flags.bad = ["base_PsfFlux_flag", ]
75 config.update(**kwargs)
77 return self.subtractTask(config=config)
80class AlardLuptonSubtractTest(AlardLuptonSubtractTestBase, lsst.utils.tests.TestCase):
81 subtractTask = subtractImages.AlardLuptonSubtractTask
83 def test_allowed_config_modes(self):
84 """Verify the allowable modes for convolution.
85 """
86 config = subtractImages.AlardLuptonSubtractTask.ConfigClass()
87 config.mode = 'auto'
88 config.mode = 'convolveScience'
89 config.mode = 'convolveTemplate'
91 with self.assertRaises(FieldValidationError):
92 config.mode = 'aotu'
94 def test_mismatched_template(self):
95 """Test that an error is raised if the template
96 does not fully contain the science image.
97 """
98 xSize = 200
99 ySize = 200
100 science, sources = makeTestImage(psfSize=self.midPsfSize, xSize=xSize + 20, ySize=ySize + 20)
101 template, _ = makeTestImage(psfSize=self.midPsfSize, xSize=xSize, ySize=ySize,
102 doApplyCalibration=True)
103 task = self._setup_subtraction()
104 with self.assertRaises(AssertionError):
105 task.run(template, science, sources)
107 def test_mismatched_filter(self):
108 """Test that an error is raised if the science and template have
109 different bands.
110 """
111 xSize = 200
112 ySize = 200
113 science, sources = makeTestImage(psfSize=self.midPsfSize, xSize=xSize + 20, ySize=ySize + 20,
114 band="g", physicalFilter="g noCamera")
115 template, _ = makeTestImage(psfSize=self.midPsfSize, xSize=xSize, ySize=ySize,
116 doApplyCalibration=True, band="not-g", physicalFilter="not-g noCamera")
117 task = self._setup_subtraction()
118 with self.assertRaises(AssertionError):
119 task.run(template, science, sources)
121 def test_incomplete_template_coverage(self):
122 noiseLevel = 1.
123 border = 20
124 xSize = 400
125 ySize = 400
126 nSources = 80
127 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6,
128 nSrc=nSources, xSize=xSize, ySize=ySize)
129 template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
130 nSrc=nSources, templateBorderSize=border, doApplyCalibration=True,
131 xSize=xSize, ySize=ySize)
133 science_height = science.getBBox().getDimensions().getY()
135 def _run_and_check_coverage(template_coverage,
136 requiredTemplateFraction=0.1,
137 minTemplateFractionForExpectedSuccess=0.2):
138 template_cut = template.clone()
139 template_height = int(science_height*template_coverage + border)
140 template_cut.image.array[:, template_height:] = 0
141 template_cut.mask.array[:, template_height:] = template_cut.mask.getPlaneBitMask('NO_DATA')
142 task = self._setup_subtraction(
143 requiredTemplateFraction=requiredTemplateFraction,
144 minTemplateFractionForExpectedSuccess=minTemplateFractionForExpectedSuccess
145 )
146 if template_coverage < requiredTemplateFraction:
147 doRaise = True
148 elif template_coverage < minTemplateFractionForExpectedSuccess:
149 doRaise = True
150 else:
151 doRaise = False
152 if doRaise:
153 with self.assertRaises(NoWorkFound):
154 task.run(template_cut, science.clone(), sources.copy(deep=True))
155 else:
156 task.run(template_cut, science.clone(), sources.copy(deep=True))
157 _run_and_check_coverage(template_coverage=0.09)
158 _run_and_check_coverage(template_coverage=0.15)
159 _run_and_check_coverage(template_coverage=0.7)
161 def test_clear_template_mask(self):
162 noiseLevel = 1.
163 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6)
164 template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
165 templateBorderSize=20, doApplyCalibration=True)
166 diffimEmptyMaskPlanes = ["DETECTED", "DETECTED_NEGATIVE"]
167 task = self._setup_subtraction(mode="convolveTemplate")
168 # Ensure that each each mask plane is set for some pixels
169 mask = template.mask
170 x0 = 50
171 x1 = 75
172 y0 = 150
173 y1 = 175
174 scienceMaskCheck = {}
175 for maskPlane in mask.getMaskPlaneDict().keys():
176 scienceMaskCheck[maskPlane] = np.sum(science.mask.array & mask.getPlaneBitMask(maskPlane) > 0)
177 mask.array[x0: x1, y0: y1] |= mask.getPlaneBitMask(maskPlane)
178 self.assertTrue(np.sum(mask.array & mask.getPlaneBitMask(maskPlane) > 0))
180 output = task.run(template, science, sources)
181 # Verify that the template mask has been modified in place
182 for maskPlane in mask.getMaskPlaneDict().keys():
183 if maskPlane in diffimEmptyMaskPlanes:
184 self.assertTrue(np.sum(mask.array & mask.getPlaneBitMask(maskPlane) == 0))
185 elif maskPlane in task.config.preserveTemplateMask:
186 self.assertTrue(np.sum(mask.array & mask.getPlaneBitMask(maskPlane) > 0))
187 else:
188 self.assertTrue(np.sum(mask.array & mask.getPlaneBitMask(maskPlane) == 0))
189 # Mask planes set in the science image should also be set in the difference
190 # Except the "DETECTED" planes should have been cleared
191 diffimMask = output.difference.mask
192 for maskPlane, scienceSum in scienceMaskCheck.items():
193 diffimSum = np.sum(diffimMask.array & mask.getPlaneBitMask(maskPlane) > 0)
194 if maskPlane in diffimEmptyMaskPlanes:
195 self.assertEqual(diffimSum, 0)
196 else:
197 self.assertTrue(diffimSum >= scienceSum)
199 def test_equal_images(self):
200 """Test that running with enough sources produces reasonable output,
201 with the same size psf in the template and science.
202 """
203 noiseLevel = 1.
204 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6)
205 template, _ = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
206 templateBorderSize=20, doApplyCalibration=True)
207 task = self._setup_subtraction()
208 output = task.run(template, science, sources)
209 # There shoud be no NaN values in the difference image
210 self.assertTrue(np.all(np.isfinite(output.difference.image.array)))
211 # Mean of difference image should be close to zero.
212 meanError = noiseLevel/np.sqrt(output.difference.image.array.size)
213 # Make sure to include pixels with the DETECTED mask bit set.
214 statsCtrl = makeStats(badMaskPlanes=("EDGE", "BAD", "NO_DATA", "DETECTED", "DETECTED_NEGATIVE"))
215 differenceMean = computeRobustStatistics(output.difference.image, output.difference.mask, statsCtrl)
216 self.assertFloatsAlmostEqual(differenceMean, 0, atol=5*meanError)
217 # stddev of difference image should be close to expected value.
218 differenceStd = computeRobustStatistics(output.difference.image, output.difference.mask,
219 makeStats(), statistic=afwMath.STDEV)
220 self.assertFloatsAlmostEqual(differenceStd, np.sqrt(2)*noiseLevel, rtol=0.1)
222 def test_equal_images_missing_mask_planes(self):
223 """Test that running with enough sources produces reasonable output,
224 with the same size psf in the template and science and with missing
225 mask planes.
226 """
227 noiseLevel = 1.
228 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6,
229 addMaskPlanes=[])
230 template, _ = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
231 templateBorderSize=20, doApplyCalibration=True, addMaskPlanes=[])
232 task = self._setup_subtraction()
233 output = task.run(template, science, sources)
234 # There shoud be no NaN values in the difference image
235 self.assertTrue(np.all(np.isfinite(output.difference.image.array)))
236 # Mean of difference image should be close to zero.
237 meanError = noiseLevel/np.sqrt(output.difference.image.array.size)
238 # Make sure to include pixels with the DETECTED mask bit set.
239 statsCtrl = makeStats(badMaskPlanes=("EDGE", "BAD", "NO_DATA", "DETECTED", "DETECTED_NEGATIVE"))
240 differenceMean = computeRobustStatistics(output.difference.image, output.difference.mask, statsCtrl)
241 self.assertFloatsAlmostEqual(differenceMean, 0, atol=5*meanError)
242 # stddev of difference image should be close to expected value.
243 differenceStd = computeRobustStatistics(output.difference.image, output.difference.mask,
244 makeStats(), statistic=afwMath.STDEV)
245 self.assertFloatsAlmostEqual(differenceStd, np.sqrt(2)*noiseLevel, rtol=0.1)
247 def test_psf_size(self):
248 """Test that the image subtract task runs without failing, if
249 fwhmExposureBuffer and fwhmExposureGrid parameters are set.
250 """
251 noiseLevel = 1.
252 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6)
253 template, _ = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
254 templateBorderSize=20, doApplyCalibration=True)
256 schema = afwTable.ExposureTable.makeMinimalSchema()
257 weightKey = schema.addField("weight", type="D", doc="Coadd weight")
258 exposureCatalog = afwTable.ExposureCatalog(schema)
259 kernel = measAlg.DoubleGaussianPsf(7, 7, 2.0).getKernel()
260 psf = measAlg.KernelPsf(kernel, template.getBBox().getCenter())
262 record = exposureCatalog.addNew()
263 record.setPsf(psf)
264 record.setWcs(template.wcs)
265 record.setD(weightKey, 1.0)
266 record.setBBox(template.getBBox())
268 customPsf = CustomCoaddPsf(exposureCatalog, template.wcs)
269 template.setPsf(customPsf)
271 # Test that we get an exception if we simply get the FWHM at center.
272 with self.assertRaises(InvalidParameterError):
273 getPsfFwhm(template.psf, True)
275 with self.assertRaises(InvalidParameterError):
276 getPsfFwhm(template.psf, False)
278 # Test that evaluateMeanPsfFwhm runs successfully on the template.
279 evaluateMeanPsfFwhm(template, fwhmExposureBuffer=0.05, fwhmExposureGrid=10)
281 # Since the PSF is spatially invariant, the FWHM should be the same at
282 # all points in the science image.
283 fwhm1 = getPsfFwhm(science.psf, False)
284 fwhm2 = evaluateMeanPsfFwhm(science, fwhmExposureBuffer=0.05, fwhmExposureGrid=10)
285 self.assertAlmostEqual(fwhm1[0], fwhm2, places=13)
286 self.assertAlmostEqual(fwhm1[1], fwhm2, places=13)
288 self.assertAlmostEqual(evaluateMeanPsfFwhm(science, fwhmExposureBuffer=0.05,
289 fwhmExposureGrid=10),
290 getPsfFwhm(science.psf, True), places=7
291 )
293 # Test that the image subtraction task runs successfully.
294 task = self._setup_subtraction()
296 # Test that the task runs if we take the mean FWHM on a grid.
297 with self.assertLogs(level="INFO") as cm:
298 task.run(template, science, sources)
300 # Check that evaluateMeanPsfFwhm was called.
301 # This tests that getPsfFwhm failed raising InvalidParameterError,
302 # that is caught and handled appropriately.
303 logMessage = ("INFO:lsst.alardLuptonSubtract:Unable to evaluate PSF at the average position. "
304 "Evaluting PSF on a grid of points."
305 )
306 self.assertIn(logMessage, cm.output)
308 def test_auto_convolveTemplate(self):
309 """Test that auto mode gives the same result as convolveTemplate when
310 the template psf is the smaller.
311 """
312 noiseLevel = 1.
313 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6)
314 template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
315 templateBorderSize=20, doApplyCalibration=True)
316 task = self._setup_subtraction(mode="convolveTemplate")
317 output = task.run(template.clone(), science.clone(), sources)
319 task = self._setup_subtraction(mode="auto")
320 outputAuto = task.run(template, science, sources)
321 self.assertMaskedImagesEqual(output.difference.maskedImage, outputAuto.difference.maskedImage)
323 def test_auto_convolveScience(self):
324 """Test that auto mode gives the same result as convolveScience when
325 the science psf is the smaller.
326 """
327 noiseLevel = 1.
328 science, sources = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=6)
329 template, _ = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
330 templateBorderSize=20, doApplyCalibration=True)
331 task = self._setup_subtraction(mode="convolveScience")
332 output = task.run(template.clone(), science.clone(), sources)
334 task = self._setup_subtraction(mode="auto")
335 outputAuto = task.run(template, science, sources)
336 self.assertMaskedImagesEqual(output.difference.maskedImage, outputAuto.difference.maskedImage)
338 def test_science_better(self):
339 """Test that running with enough sources produces reasonable output,
340 with the science psf being smaller than the template.
341 """
342 statsCtrl = makeStats()
343 statsCtrlDetect = makeStats(badMaskPlanes=("EDGE", "BAD", "NO_DATA"))
345 def _run_and_check_images(statsCtrl, statsCtrlDetect, scienceNoiseLevel, templateNoiseLevel):
346 science, sources = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=scienceNoiseLevel,
347 noiseSeed=6)
348 template, _ = makeTestImage(psfSize=self.midPsfSize, noiseLevel=templateNoiseLevel, noiseSeed=7,
349 templateBorderSize=20, doApplyCalibration=True)
350 task = self._setup_subtraction(mode="convolveScience")
351 output = task.run(template, science, sources)
352 self.assertFloatsAlmostEqual(task.metadata["scaleTemplateVarianceFactor"], 1., atol=.05)
353 self.assertFloatsAlmostEqual(task.metadata["scaleScienceVarianceFactor"], 1., atol=.05)
354 # Mean of difference image should be close to zero.
355 nGoodPix = np.sum(np.isfinite(output.difference.image.array))
356 meanError = (scienceNoiseLevel + templateNoiseLevel)/np.sqrt(nGoodPix)
357 diffimMean = computeRobustStatistics(output.difference.image, output.difference.mask,
358 statsCtrlDetect)
360 self.assertFloatsAlmostEqual(diffimMean, 0, atol=5*meanError)
361 # stddev of difference image should be close to expected value.
362 noiseLevel = np.sqrt(scienceNoiseLevel**2 + templateNoiseLevel**2)
363 varianceMean = computeRobustStatistics(output.difference.variance, output.difference.mask,
364 statsCtrl)
365 diffimStd = computeRobustStatistics(output.difference.image, output.difference.mask,
366 statsCtrl, statistic=afwMath.STDEV)
367 self.assertFloatsAlmostEqual(varianceMean, noiseLevel**2, rtol=0.1)
368 self.assertFloatsAlmostEqual(diffimStd, noiseLevel, rtol=0.1)
370 _run_and_check_images(statsCtrl, statsCtrlDetect, scienceNoiseLevel=1., templateNoiseLevel=1.)
371 _run_and_check_images(statsCtrl, statsCtrlDetect, scienceNoiseLevel=1., templateNoiseLevel=.1)
372 _run_and_check_images(statsCtrl, statsCtrlDetect, scienceNoiseLevel=.1, templateNoiseLevel=.1)
374 def test_template_better(self):
375 """Test that running with enough sources produces reasonable output,
376 with the template psf being smaller than the science.
377 """
378 statsCtrl = makeStats()
379 statsCtrlDetect = makeStats(badMaskPlanes=("EDGE", "BAD", "NO_DATA"))
381 def _run_and_check_images(statsCtrl, statsCtrlDetect, scienceNoiseLevel, templateNoiseLevel):
382 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=scienceNoiseLevel,
383 noiseSeed=6)
384 template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=templateNoiseLevel, noiseSeed=7,
385 templateBorderSize=20, doApplyCalibration=True)
386 task = self._setup_subtraction()
387 output = task.run(template, science, sources)
388 self.assertFloatsAlmostEqual(task.metadata["scaleTemplateVarianceFactor"], 1., atol=.05)
389 self.assertFloatsAlmostEqual(task.metadata["scaleScienceVarianceFactor"], 1., atol=.05)
390 # There should be no NaNs in the image if we convolve the template with a buffer
391 self.assertTrue(np.all(np.isfinite(output.difference.image.array)))
392 # Mean of difference image should be close to zero.
393 meanError = (scienceNoiseLevel + templateNoiseLevel)/np.sqrt(output.difference.image.array.size)
395 diffimMean = computeRobustStatistics(output.difference.image, output.difference.mask,
396 statsCtrlDetect)
397 self.assertFloatsAlmostEqual(diffimMean, 0, atol=5*meanError)
398 # stddev of difference image should be close to expected value.
399 noiseLevel = np.sqrt(scienceNoiseLevel**2 + templateNoiseLevel**2)
400 varianceMean = computeRobustStatistics(output.difference.variance, output.difference.mask,
401 statsCtrl)
402 diffimStd = computeRobustStatistics(output.difference.image, output.difference.mask,
403 statsCtrl, statistic=afwMath.STDEV)
404 self.assertFloatsAlmostEqual(varianceMean, noiseLevel**2, rtol=0.1)
405 self.assertFloatsAlmostEqual(diffimStd, noiseLevel, rtol=0.1)
407 _run_and_check_images(statsCtrl, statsCtrlDetect, scienceNoiseLevel=1., templateNoiseLevel=1.)
408 _run_and_check_images(statsCtrl, statsCtrlDetect, scienceNoiseLevel=1., templateNoiseLevel=.1)
409 _run_and_check_images(statsCtrl, statsCtrlDetect, scienceNoiseLevel=.1, templateNoiseLevel=.1)
411 def test_symmetry(self):
412 """Test that convolving the science and convolving the template are
413 symmetric: if the psfs are switched between them, the difference image
414 should be nearly the same.
415 """
416 noiseLevel = 1.
417 # Don't include a border for the template, in order to make the results
418 # comparable when we swap which image is treated as the "science" image.
419 science, sources = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel,
420 noiseSeed=6, templateBorderSize=0, doApplyCalibration=True)
421 template, _ = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel,
422 noiseSeed=7, templateBorderSize=0, doApplyCalibration=True)
423 task = self._setup_subtraction(mode='auto')
425 # The science image will be modified in place, so use a copy for the second run.
426 science_better = task.run(template.clone(), science.clone(), sources)
427 template_better = task.run(science, template, sources)
429 delta = template_better.difference.clone()
430 delta.image -= science_better.difference.image
431 delta.variance -= science_better.difference.variance
432 delta.mask.array -= science_better.difference.mask.array
434 statsCtrl = makeStats()
435 # Mean of delta should be very close to zero.
436 nGoodPix = np.sum(np.isfinite(delta.image.array))
437 meanError = 2*noiseLevel/np.sqrt(nGoodPix)
438 deltaMean = computeRobustStatistics(delta.image, delta.mask, statsCtrl)
439 deltaStd = computeRobustStatistics(delta.image, delta.mask, statsCtrl, statistic=afwMath.STDEV)
440 self.assertFloatsAlmostEqual(deltaMean, 0, atol=5*meanError)
441 # stddev of difference image should be close to expected value
442 self.assertFloatsAlmostEqual(deltaStd, 2*np.sqrt(2)*noiseLevel, rtol=.1)
444 def test_few_sources(self):
445 """Test with only 1 source, to check that we get a useful error.
446 """
447 xSize = 256
448 ySize = 256
449 science, sources = makeTestImage(psfSize=self.midPsfSize, nSrc=10, xSize=xSize, ySize=ySize)
450 template, _ = makeTestImage(psfSize=self.goodPsfSize, nSrc=10, xSize=xSize, ySize=ySize,
451 doApplyCalibration=True)
452 task = self._setup_subtraction()
453 sources = sources[0:1]
454 with self.assertRaises(InsufficientKernelSourcesError):
455 task.run(template, science, sources)
457 def test_kernel_source_selector(self):
458 """Check that kernel source selection behaves as expected.
459 """
460 xSize = 256
461 ySize = 256
462 nSourcesSimulated = 20
463 sciencePsfSize = self.midPsfSize
464 templatePsfSize = self.goodPsfSize
465 science, sources = makeTestImage(psfSize=sciencePsfSize, nSrc=nSourcesSimulated,
466 xSize=xSize, ySize=ySize)
467 template, _ = makeTestImage(psfSize=templatePsfSize, nSrc=nSourcesSimulated,
468 xSize=xSize, ySize=ySize, doApplyCalibration=True)
470 def _run_and_check_sources(sourcesIn, maxKernelSources=1000, minKernelSources=3):
471 sources = sourcesIn.copy(deep=True)
473 task = self._setup_subtraction(maxKernelSources=maxKernelSources,
474 minKernelSources=minKernelSources,
475 )
476 task.templatePsfSize = templatePsfSize
477 task.sciencePsfSize = sciencePsfSize
478 task.matchedPsfSize = sciencePsfSize
479 # Verify that source flags are not set in the input catalog
480 # Note that this will use the last flag in the list for the rest of
481 # the test.
482 for badSourceFlag in task.sourceSelector.config.flags.bad:
483 self.assertEqual(np.sum(sources[badSourceFlag]), 0)
484 nSources = len(sources)
485 # Flag a third of the sources
486 sources[0:: 3][badSourceFlag] = True
487 rejectRadius = 2*task.config.makeKernel.kernel.active.kernelSize
488 bbox = science.getBBox()
489 bbox.grow(-rejectRadius)
490 edgeSources = ~bbox.contains(sources.getX(), sources.getY())
491 sources[edgeSources][badSourceFlag] = True
492 nBadSources = np.sum(sources[badSourceFlag])
493 if maxKernelSources > 0:
494 nGoodSources = np.minimum(nSources - nBadSources, maxKernelSources)
495 else:
496 nGoodSources = nSources - nBadSources
498 signalToNoise = sources.getPsfInstFlux()/sources.getPsfInstFluxErr()
499 signalToNoise = signalToNoise[~sources[badSourceFlag]]
500 signalToNoise.sort()
501 selectSources = task._sourceSelector(template, science, sources)
502 self.assertEqual(nGoodSources, len(selectSources))
503 signalToNoiseOut = selectSources.getPsfInstFlux()/selectSources.getPsfInstFluxErr()
504 signalToNoiseOut.sort()
505 self.assertFloatsAlmostEqual(signalToNoise[-nGoodSources:], signalToNoiseOut)
507 _run_and_check_sources(sources)
508 _run_and_check_sources(sources, maxKernelSources=len(sources)//3)
509 _run_and_check_sources(sources, maxKernelSources=-1)
510 with self.assertRaises(RuntimeError):
511 _run_and_check_sources(sources, minKernelSources=1000)
513 def test_order_equal_images(self):
514 """Verify that the result is the same regardless of convolution mode
515 if the images are equivalent.
516 """
517 noiseLevel = .1
518 seed1 = 6
519 seed2 = 7
520 for psfSize in [self.midPsfSize, self.goodPsfSize, self.badPsfSize]:
521 science1, sources1 = makeTestImage(psfSize=psfSize, noiseLevel=noiseLevel, noiseSeed=seed1,
522 clearEdgeMask=True)
523 template1, _ = makeTestImage(psfSize=psfSize, noiseLevel=noiseLevel, noiseSeed=seed2,
524 templateBorderSize=0, doApplyCalibration=True,
525 clearEdgeMask=True)
526 task1 = self._setup_subtraction(mode="convolveTemplate")
527 results_convolveTemplate = task1.run(template1, science1, sources1)
529 science2, sources2 = makeTestImage(psfSize=psfSize, noiseLevel=noiseLevel, noiseSeed=seed1,
530 clearEdgeMask=True)
531 template2, _ = makeTestImage(psfSize=psfSize, noiseLevel=noiseLevel, noiseSeed=seed2,
532 templateBorderSize=0, doApplyCalibration=True,
533 clearEdgeMask=True)
534 task2 = self._setup_subtraction(mode="convolveScience")
535 results_convolveScience = task2.run(template2, science2, sources2)
536 bbox = results_convolveTemplate.difference.getBBox().clippedTo(
537 results_convolveScience.difference.getBBox())
538 diff1 = science1.maskedImage.clone()[bbox]
539 diff1 -= template1.maskedImage[bbox]
540 diff2 = science2.maskedImage.clone()[bbox]
541 diff2 -= template2.maskedImage[bbox]
542 self.assertFloatsAlmostEqual(results_convolveTemplate.difference[bbox].image.array,
543 diff1.image.array,
544 atol=noiseLevel*5.)
545 self.assertFloatsAlmostEqual(results_convolveScience.difference[bbox].image.array,
546 diff2.image.array,
547 atol=noiseLevel*5.)
548 diffErr = noiseLevel*2
549 self.assertMaskedImagesAlmostEqual(results_convolveTemplate.difference[bbox].maskedImage,
550 results_convolveScience.difference[bbox].maskedImage,
551 atol=diffErr*5.)
553 def test_background_subtraction(self):
554 """Check that we can recover the background,
555 and that it is subtracted correctly in the difference image.
557 NOTE: Background subtraction is now turned off by default in
558 subtractImages. It is now run in detectAndMeasure instead, but since the
559 code to run background subtraction is not being removed this test should
560 stay to make sure it continues functioning as intended.
561 """
562 noiseLevel = 1.
563 xSize = 512
564 ySize = 512
565 x0 = 123
566 y0 = 456
567 template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
568 templateBorderSize=20,
569 xSize=xSize, ySize=ySize, x0=x0, y0=y0,
570 doApplyCalibration=True)
571 params = [2.2, 2.1, 2.0, 1.2, 1.1, 1.0]
573 bbox2D = lsst.geom.Box2D(lsst.geom.Point2D(x0, y0), lsst.geom.Extent2D(xSize, ySize))
574 background_model = afwMath.Chebyshev1Function2D(params, bbox2D)
575 science, sources = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=6,
576 background=background_model,
577 xSize=xSize, ySize=ySize, x0=x0, y0=y0)
578 # Don't use ``self._setup_subtraction()`` here.
579 # Modifying the config of a subtask is messy.
580 config = subtractImages.AlardLuptonSubtractTask.ConfigClass()
582 config.sourceSelector.signalToNoise.fluxField = "truth_instFlux"
583 config.sourceSelector.signalToNoise.errField = "truth_instFluxErr"
584 config.doSubtractBackground = True
586 config.makeKernel.kernel.name = "AL"
587 config.makeKernel.kernel.active.fitForBackground = True
588 config.makeKernel.kernel.active.spatialKernelOrder = 1
589 config.makeKernel.kernel.active.spatialBgOrder = 2
590 statsCtrl = makeStats()
592 def _run_and_check_images(config, statsCtrl, mode):
593 """Check that the fit background matches the input model.
594 """
595 config.mode = mode
596 task = subtractImages.AlardLuptonSubtractTask(config=config)
597 output = task.run(template.clone(), science.clone(), sources)
599 # We should be fitting the same number of parameters as were in the input model
600 self.assertEqual(output.backgroundModel.getNParameters(), background_model.getNParameters())
602 # The parameters of the background fit should be close to the input model
603 self.assertFloatsAlmostEqual(np.array(output.backgroundModel.getParameters()),
604 np.array(params), rtol=0.3)
606 # stddev of difference image should be close to expected value.
607 # This will fail if we have mis-subtracted the background.
608 stdVal = computeRobustStatistics(output.difference.image, output.difference.mask,
609 statsCtrl, statistic=afwMath.STDEV)
610 self.assertFloatsAlmostEqual(stdVal, np.sqrt(2)*noiseLevel, rtol=0.1)
612 _run_and_check_images(config, statsCtrl, "convolveTemplate")
613 _run_and_check_images(config, statsCtrl, "convolveScience")
615 def test_scale_variance_convolve_template(self):
616 """Check variance scaling of the image difference.
617 """
618 scienceNoiseLevel = 4.
619 templateNoiseLevel = 2.
620 scaleFactor = 1.345
621 # Make sure to include pixels with the DETECTED mask bit set.
622 statsCtrl = makeStats(badMaskPlanes=("EDGE", "BAD", "NO_DATA"))
624 def _run_and_check_images(science, template, sources, statsCtrl,
625 doDecorrelation, doScaleVariance, scaleFactor=1.):
626 """Check that the variance plane matches the expected value for
627 different configurations of ``doDecorrelation`` and ``doScaleVariance``.
628 """
630 task = self._setup_subtraction(doDecorrelation=doDecorrelation,
631 doScaleVariance=doScaleVariance,
632 )
633 output = task.run(template.clone(), science.clone(), sources)
634 if doScaleVariance:
635 self.assertFloatsAlmostEqual(task.metadata["scaleTemplateVarianceFactor"],
636 scaleFactor, atol=0.05)
637 self.assertFloatsAlmostEqual(task.metadata["scaleScienceVarianceFactor"],
638 scaleFactor, atol=0.05)
640 scienceNoise = computeRobustStatistics(science.variance, science.mask, statsCtrl)
641 if doDecorrelation:
642 templateNoise = computeRobustStatistics(template.variance, template.mask, statsCtrl)
643 else:
644 templateNoise = computeRobustStatistics(output.matchedTemplate.variance,
645 output.matchedTemplate.mask,
646 statsCtrl)
648 if doScaleVariance:
649 templateNoise *= scaleFactor
650 scienceNoise *= scaleFactor
651 varMean = computeRobustStatistics(output.difference.variance, output.difference.mask, statsCtrl)
652 self.assertFloatsAlmostEqual(varMean, scienceNoise + templateNoise, rtol=0.1)
654 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=scienceNoiseLevel, noiseSeed=6)
655 template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=templateNoiseLevel, noiseSeed=7,
656 templateBorderSize=20, doApplyCalibration=True)
657 # Verify that the variance plane of the difference image is correct
658 # when the template and science variance planes are correct
659 _run_and_check_images(science, template, sources, statsCtrl,
660 doDecorrelation=True, doScaleVariance=True)
661 _run_and_check_images(science, template, sources, statsCtrl,
662 doDecorrelation=True, doScaleVariance=False)
663 _run_and_check_images(science, template, sources, statsCtrl,
664 doDecorrelation=False, doScaleVariance=True)
665 _run_and_check_images(science, template, sources, statsCtrl,
666 doDecorrelation=False, doScaleVariance=False)
668 # Verify that the variance plane of the difference image is correct
669 # when the template variance plane is incorrect
670 template.variance.array /= scaleFactor
671 science.variance.array /= scaleFactor
672 _run_and_check_images(science, template, sources, statsCtrl,
673 doDecorrelation=True, doScaleVariance=True, scaleFactor=scaleFactor)
674 _run_and_check_images(science, template, sources, statsCtrl,
675 doDecorrelation=True, doScaleVariance=False, scaleFactor=scaleFactor)
676 _run_and_check_images(science, template, sources, statsCtrl,
677 doDecorrelation=False, doScaleVariance=True, scaleFactor=scaleFactor)
678 _run_and_check_images(science, template, sources, statsCtrl,
679 doDecorrelation=False, doScaleVariance=False, scaleFactor=scaleFactor)
681 def test_scale_variance_convolve_science(self):
682 """Check variance scaling of the image difference.
683 """
684 scienceNoiseLevel = 4.
685 templateNoiseLevel = 2.
686 scaleFactor = 1.345
687 # Make sure to include pixels with the DETECTED mask bit set.
688 statsCtrl = makeStats(badMaskPlanes=("EDGE", "BAD", "NO_DATA"))
690 def _run_and_check_images(science, template, sources, statsCtrl,
691 doDecorrelation, doScaleVariance, scaleFactor=1.):
692 """Check that the variance plane matches the expected value for
693 different configurations of ``doDecorrelation`` and ``doScaleVariance``.
694 """
696 task = self._setup_subtraction(mode="convolveScience",
697 doDecorrelation=doDecorrelation,
698 doScaleVariance=doScaleVariance,
699 restrictKernelEdgeSources=False,
700 )
701 output = task.run(template.clone(), science.clone(), sources)
702 if doScaleVariance:
703 self.assertFloatsAlmostEqual(task.metadata["scaleTemplateVarianceFactor"],
704 scaleFactor, atol=0.05)
705 self.assertFloatsAlmostEqual(task.metadata["scaleScienceVarianceFactor"],
706 scaleFactor, atol=0.05)
708 templateNoise = computeRobustStatistics(template.variance, template.mask, statsCtrl)
709 if doDecorrelation:
710 scienceNoise = computeRobustStatistics(science.variance, science.mask, statsCtrl)
711 else:
712 scienceNoise = computeRobustStatistics(output.matchedScience.variance,
713 output.matchedScience.mask,
714 statsCtrl)
716 if doScaleVariance:
717 templateNoise *= scaleFactor
718 scienceNoise *= scaleFactor
720 varMean = computeRobustStatistics(output.difference.variance, output.difference.mask, statsCtrl)
721 self.assertFloatsAlmostEqual(varMean, scienceNoise + templateNoise, rtol=0.1)
723 science, sources = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=scienceNoiseLevel, noiseSeed=6)
724 template, _ = makeTestImage(psfSize=self.midPsfSize, noiseLevel=templateNoiseLevel, noiseSeed=7,
725 templateBorderSize=20, doApplyCalibration=True)
726 # Verify that the variance plane of the difference image is correct
727 # when the template and science variance planes are correct
728 _run_and_check_images(science, template, sources, statsCtrl,
729 doDecorrelation=True, doScaleVariance=True)
730 _run_and_check_images(science, template, sources, statsCtrl,
731 doDecorrelation=True, doScaleVariance=False)
732 _run_and_check_images(science, template, sources, statsCtrl,
733 doDecorrelation=False, doScaleVariance=True)
734 _run_and_check_images(science, template, sources, statsCtrl,
735 doDecorrelation=False, doScaleVariance=False)
737 # Verify that the variance plane of the difference image is correct
738 # when the template and science variance planes are incorrect
739 science.variance.array /= scaleFactor
740 template.variance.array /= scaleFactor
741 _run_and_check_images(science, template, sources, statsCtrl,
742 doDecorrelation=True, doScaleVariance=True, scaleFactor=scaleFactor)
743 _run_and_check_images(science, template, sources, statsCtrl,
744 doDecorrelation=True, doScaleVariance=False, scaleFactor=scaleFactor)
745 _run_and_check_images(science, template, sources, statsCtrl,
746 doDecorrelation=False, doScaleVariance=True, scaleFactor=scaleFactor)
747 _run_and_check_images(science, template, sources, statsCtrl,
748 doDecorrelation=False, doScaleVariance=False, scaleFactor=scaleFactor)
750 def test_exposure_properties_convolve_template(self):
751 """Check that all necessary exposure metadata is included
752 when the template is convolved.
753 """
754 noiseLevel = 1.
755 seed = 37
756 rng = np.random.RandomState(seed)
757 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6)
758 psf = science.psf
759 psfAvgPos = psf.getAveragePosition()
760 psfSize = getPsfFwhm(science.psf)
761 psfImg = psf.computeKernelImage(psfAvgPos)
762 template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
763 templateBorderSize=20, doApplyCalibration=True)
765 # Generate a random aperture correction map
766 apCorrMap = lsst.afw.image.ApCorrMap()
767 for name in ("a", "b", "c"):
768 apCorrMap.set(name, lsst.afw.math.ChebyshevBoundedField(science.getBBox(), rng.randn(3, 3)))
769 science.info.setApCorrMap(apCorrMap)
771 def _run_and_check_images(doDecorrelation):
772 """Check that the metadata is correct with or without decorrelation.
773 """
774 task = self._setup_subtraction(mode="convolveTemplate",
775 doDecorrelation=doDecorrelation,
776 )
777 output = task.run(template.clone(), science.clone(), sources)
778 psfOut = output.difference.psf
779 psfAvgPos = psfOut.getAveragePosition()
780 if doDecorrelation:
781 # Decorrelation requires recalculating the PSF,
782 # so it will not be the same as the input
783 psfOutSize = getPsfFwhm(science.psf)
784 self.assertFloatsAlmostEqual(psfSize, psfOutSize)
785 else:
786 psfOutImg = psfOut.computeKernelImage(psfAvgPos)
787 self.assertImagesAlmostEqual(psfImg, psfOutImg)
789 # check PSF, WCS, bbox, filterLabel, photoCalib, aperture correction
790 self._compare_apCorrMaps(apCorrMap, output.difference.info.getApCorrMap())
791 self.assertWcsAlmostEqualOverBBox(science.wcs, output.difference.wcs, science.getBBox())
792 self.assertEqual(science.filter, output.difference.filter)
793 self.assertEqual(science.photoCalib, output.difference.photoCalib)
794 _run_and_check_images(doDecorrelation=True)
795 _run_and_check_images(doDecorrelation=False)
797 def test_exposure_properties_convolve_science(self):
798 """Check that all necessary exposure metadata is included
799 when the science image is convolved.
800 """
801 noiseLevel = 1.
802 seed = 37
803 rng = np.random.RandomState(seed)
804 science, sources = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=6)
805 template, _ = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
806 templateBorderSize=20, doApplyCalibration=True)
807 psf = template.psf
808 psfAvgPos = psf.getAveragePosition()
809 psfSize = getPsfFwhm(template.psf)
810 psfImg = psf.computeKernelImage(psfAvgPos)
812 # Generate a random aperture correction map
813 apCorrMap = lsst.afw.image.ApCorrMap()
814 for name in ("a", "b", "c"):
815 apCorrMap.set(name, lsst.afw.math.ChebyshevBoundedField(science.getBBox(), rng.randn(3, 3)))
816 science.info.setApCorrMap(apCorrMap)
818 def _run_and_check_images(doDecorrelation):
819 """Check that the metadata is correct with or without decorrelation.
820 """
821 task = self._setup_subtraction(mode="convolveScience",
822 doDecorrelation=doDecorrelation,
823 )
824 output = task.run(template.clone(), science.clone(), sources)
825 if doDecorrelation:
826 # Decorrelation requires recalculating the PSF,
827 # so it will not be the same as the input
828 psfOutSize = getPsfFwhm(template.psf)
829 self.assertFloatsAlmostEqual(psfSize, psfOutSize)
830 else:
831 psfOut = output.difference.psf
832 psfAvgPos = psfOut.getAveragePosition()
833 psfOutImg = psfOut.computeKernelImage(psfAvgPos)
834 self.assertImagesAlmostEqual(psfImg, psfOutImg)
836 # check PSF, WCS, bbox, filterLabel, photoCalib, aperture correction
837 self._compare_apCorrMaps(apCorrMap, output.difference.info.getApCorrMap())
838 self.assertWcsAlmostEqualOverBBox(science.wcs, output.difference.wcs, science.getBBox())
839 self.assertEqual(science.filter, output.difference.filter)
840 self.assertEqual(science.photoCalib, output.difference.photoCalib)
842 _run_and_check_images(doDecorrelation=True)
843 _run_and_check_images(doDecorrelation=False)
845 def _compare_apCorrMaps(self, a, b):
846 """Compare two ApCorrMaps for equality, without assuming that their BoundedFields have the
847 same addresses (i.e. so we can compare after serialization).
849 This function is taken from ``ApCorrMapTestCase`` in afw/tests/.
851 Parameters
852 ----------
853 a, b : `lsst.afw.image.ApCorrMap`
854 The two aperture correction maps to compare.
855 """
856 self.assertEqual(len(a), len(b))
857 for name, value in list(a.items()):
858 value2 = b.get(name)
859 self.assertIsNotNone(value2)
860 self.assertEqual(value.getBBox(), value2.getBBox())
861 self.assertFloatsAlmostEqual(
862 value.getCoefficients(), value2.getCoefficients(), rtol=0.0)
864 def test_fake_mask_plane_propagation(self):
865 """Test that we have the mask planes related to fakes in diffim images.
866 This is testing method called updateMasks
867 """
868 xSize = 200
869 ySize = 200
870 science, sources = makeTestImage(psfSize=self.midPsfSize, xSize=xSize, ySize=ySize)
871 science_fake_img, science_fake_sources = makeTestImage(
872 psfSize=self.midPsfSize, xSize=xSize, ySize=ySize, seed=7, nSrc=2, noiseLevel=0.25, fluxRange=1
873 )
874 template, _ = makeTestImage(psfSize=self.midPsfSize, xSize=xSize, ySize=ySize,
875 doApplyCalibration=True)
876 tmplt_fake_img, tmplt_fake_sources = makeTestImage(
877 psfSize=self.midPsfSize, xSize=xSize, ySize=ySize, seed=9, nSrc=2, noiseLevel=0.25, fluxRange=1
878 )
879 # created fakes and added them to the images
880 science.image.array += science_fake_img.image.array
881 template.image.array += tmplt_fake_img.image.array
883 # TODO: DM-40796 update to INJECTED names when source injection gets refactored
884 # adding mask planes to both science and template images
885 science_mask_planes = science.mask.addMaskPlane("FAKE")
886 template_mask_planes = template.mask.addMaskPlane("FAKE")
888 for a_science_source in science_fake_sources:
889 # 3 x 3 masking of the source locations is fine
890 bbox = lsst.geom.Box2I(
891 lsst.geom.Point2I(a_science_source.getX(), a_science_source.getY()), lsst.geom.Extent2I(3, 3)
892 )
893 science[bbox].mask.array |= science_mask_planes
895 for a_template_source in tmplt_fake_sources:
896 # 3 x 3 masking of the source locations is fine
897 bbox = lsst.geom.Box2I(
898 lsst.geom.Point2I(a_template_source.getX(), a_template_source.getY()),
899 lsst.geom.Extent2I(3, 3)
900 )
901 template[bbox].mask.array |= template_mask_planes
903 science_fake_masked = (science.mask.array & science.mask.getPlaneBitMask("FAKE")) > 0
904 template_fake_masked = (template.mask.array & template.mask.getPlaneBitMask("FAKE")) > 0
906 task = self._setup_subtraction()
907 subtraction = task.run(template, science, sources)
909 # check subtraction mask plane is set where we set the previous masks
910 diff_mask = subtraction.difference.mask
912 # science mask should be now in INJECTED
913 inj_masked = (diff_mask.array & diff_mask.getPlaneBitMask("INJECTED")) > 0
915 # template mask should be now in INJECTED_TEMPLATE
916 injTmplt_masked = (diff_mask.array & diff_mask.getPlaneBitMask("INJECTED_TEMPLATE")) > 0
918 self.assertEqual(np.sum(inj_masked.astype(int)-science_fake_masked.astype(int)), 0)
919 self.assertEqual(np.sum(injTmplt_masked.astype(int)-template_fake_masked.astype(int)), 0)
921 def test_metadata_metrics(self):
922 """Verify fields are added to metadata when subtraction is run, and
923 that the difference image limiting magnitude is calculated correctly,
924 both with a "good" and "bad" seeing template.
925 """
926 science, sources = makeTestImage(psfSize=self.badPsfSize, noiseLevel=1)
927 template_good, _ = makeTestImage(psfSize=self.midPsfSize, doApplyCalibration=True, noiseLevel=0.25,
928 templateBorderSize=20)
929 template_bad, _ = makeTestImage(psfSize=9.5, doApplyCalibration=True, noiseLevel=0.25,
930 templateBorderSize=20)
932 # Add a few sky objects; sky footprints are needed for some metrics.
933 config = measAlg.SkyObjectsTask.ConfigClass()
934 config.nSources = 3
935 skyTask = measAlg.SkyObjectsTask(config=config, name="skySources")
936 skyTask.skySourceKey = sources.schema["sky_source"].asKey()
937 skyTask.run(science.mask, 10, catalog=sources)
938 sources = sources.copy(deep=True)
939 # Add centroids, since these sources were added post-measurement.
940 for record in sources[sources["sky_source"]]:
941 record["truth_x"] = record.getFootprint().getPeaks()[0].getFx()
942 record["truth_y"] = record.getFootprint().getPeaks()[0].getFy()
944 # The metadata fields are attached to the subtractTask, so we do
945 # need to run that; run it for both "good" and "bad" seeing templates
947 subtractTask_good = self._setup_subtraction()
948 _ = subtractTask_good.run(template_good.clone(), science.clone(), sources)
949 subtractTask_bad = self._setup_subtraction()
950 _ = subtractTask_bad.run(template_bad.clone(), science.clone(), sources)
952 # Test that the diffim limiting magnitudes are computed correctly
953 maglim_science = subtractTask_good._calculateMagLim(science)
954 fluxlim_science = (maglim_science*u.ABmag).to_value(u.nJy)
955 maglim_template_good = subtractTask_good._calculateMagLim(template_good)
956 fluxlim_template_good = (maglim_template_good*u.ABmag).to_value(u.nJy)
957 maglim_template_bad = subtractTask_bad._calculateMagLim(template_bad)
958 fluxlim_template_bad = (maglim_template_bad*u.ABmag).to_value(u.nJy)
960 maglim_good = (np.sqrt(fluxlim_science**2 + fluxlim_template_good**2)*u.nJy).to(u.ABmag).value
961 maglim_bad = (np.sqrt(fluxlim_science**2 + fluxlim_template_bad**2)*u.nJy).to(u.ABmag).value
963 self.assertFloatsAlmostEqual(subtractTask_good.metadata['diffimLimitingMagnitude'],
964 maglim_good, atol=1e-6)
965 self.assertFloatsAlmostEqual(subtractTask_bad.metadata['diffimLimitingMagnitude'],
966 maglim_bad, atol=1e-6)
968 # Create a template with a PSF that is not defined at the image center.
969 # First, make an exposure catalog so we can force the template to have
970 # a bad (off-image) PSF. It must have a record with a weight field
971 # and a BBox in order to let us set the PSF manually.
972 template_offimage, _ = makeTestImage()
973 schema = afwTable.ExposureTable.makeMinimalSchema()
974 weightKey = schema.addField("weight", type="D", doc="Coadd weight")
975 exposureCatalog = afwTable.ExposureCatalog(schema)
976 record = exposureCatalog.addNew()
977 record.setD(weightKey, 1.0)
978 record.setBBox(template_offimage.getBBox())
979 kernel = measAlg.DoubleGaussianPsf(7, 7, 2.0).getKernel()
980 psf = measAlg.KernelPsf(kernel, template_offimage.getBBox().getCenter())
981 record.setPsf(psf)
982 record.setWcs(template_offimage.wcs)
983 custom_offimage_psf = CustomCoaddPsf(exposureCatalog, template_offimage.wcs)
984 template_offimage.setPsf(custom_offimage_psf)
986 # Test that template PSF size edge cases are handled correctly.
987 subtractTask_offimage = self._setup_subtraction()
988 _ = subtractTask_offimage.run(template_offimage.clone(), science.clone(), sources)
989 # Test that providing no fallbackPsfSize results in a nan template
990 # limiting magnitude.
991 maglim_template_offimage = subtractTask_offimage._calculateMagLim(template_offimage)
992 self.assertTrue(np.isnan(maglim_template_offimage))
993 # Test that given the provided fallbackPsfSize, the diffim limiting
994 # magnitude is calculated correctly.
995 maglim_template_offimage = 28.182284789714952
996 fluxlim_template_offimage = (maglim_template_offimage*u.ABmag).to_value(u.nJy)
997 maglim_offimage = (np.sqrt(fluxlim_science**2 + fluxlim_template_offimage**2)*u.nJy).to(u.ABmag).value
998 self.assertEqual(subtractTask_offimage.metadata['diffimLimitingMagnitude'], maglim_offimage)
1000 # Test that several other expected metadata metrics exist
1001 self.assertIn('scienceLimitingMagnitude', subtractTask_good.metadata)
1002 self.assertIn('templateLimitingMagnitude', subtractTask_good.metadata)
1004 # The mean ratio metric should be much worse on the "bad" subtraction.
1005 self.assertLess(subtractTask_good.metadata['differenceFootprintRatioMean'], 0.02)
1006 self.assertGreater(subtractTask_bad.metadata['differenceFootprintRatioMean'], 0.12)
1009class AlardLuptonPreconvolveSubtractTest(AlardLuptonSubtractTestBase, lsst.utils.tests.TestCase):
1010 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask
1012 def test_mismatched_template(self):
1013 """Test that an error is raised if the template
1014 does not fully contain the science image.
1015 """
1016 xSize = 200
1017 ySize = 200
1018 science, sources = makeTestImage(psfSize=self.midPsfSize, xSize=xSize + 20, ySize=ySize + 20)
1019 template, _ = makeTestImage(psfSize=self.midPsfSize, xSize=xSize, ySize=ySize,
1020 doApplyCalibration=True)
1021 task = self._setup_subtraction()
1022 with self.assertRaises(AssertionError):
1023 task.run(template, science, sources)
1025 def test_equal_images(self):
1026 """Test that running with enough sources produces reasonable output,
1027 with the same size psf in the template and science.
1028 """
1029 noiseLevel = 1.
1030 xSize = 400
1031 ySize = 400
1032 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6,
1033 xSize=xSize, ySize=ySize)
1034 template, _ = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
1035 templateBorderSize=20, doApplyCalibration=True,
1036 xSize=xSize, ySize=ySize)
1037 task = self._setup_subtraction()
1038 output = task.run(template, science, sources)
1039 # There shoud be no NaN values in the Score image
1040 self.assertTrue(np.all(np.isfinite(output.scoreExposure.image.array)))
1041 # Mean of Score image should be close to zero.
1042 meanError = noiseLevel/np.sqrt(output.scoreExposure.image.array.size)
1043 # Make sure to include pixels with the DETECTED mask bit set.
1044 statsCtrl = makeStats(badMaskPlanes=("EDGE", "BAD", "NO_DATA"))
1045 scoreMean = computeRobustStatistics(output.scoreExposure.image,
1046 output.scoreExposure.mask,
1047 statsCtrl)
1048 self.assertFloatsAlmostEqual(scoreMean, 0, atol=5*meanError)
1049 nea = computePSFNoiseEquivalentArea(science.psf)
1050 # stddev of Score image should be close to expected value.
1051 scoreStd = computeRobustStatistics(output.scoreExposure.image, output.scoreExposure.mask,
1052 statsCtrl=statsCtrl, statistic=afwMath.STDEV)
1053 self.assertFloatsAlmostEqual(scoreStd, np.sqrt(2)*noiseLevel/np.sqrt(nea), rtol=0.1)
1055 def test_incomplete_template_coverage(self):
1056 noiseLevel = 1.
1057 border = 20
1058 xSize = 400
1059 ySize = 400
1060 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6,
1061 xSize=xSize, ySize=ySize)
1062 template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
1063 templateBorderSize=border, doApplyCalibration=True,
1064 xSize=xSize, ySize=ySize)
1066 science_height = science.getBBox().getDimensions().getY()
1068 def _run_and_check_coverage(template_coverage,
1069 requiredTemplateFraction=0.1,
1070 minTemplateFractionForExpectedSuccess=0.2):
1071 template_cut = template.clone()
1072 template_height = int(science_height*template_coverage + border)
1073 template_cut.image.array[:, template_height:] = 0
1074 template_cut.mask.array[:, template_height:] = template_cut.mask.getPlaneBitMask('NO_DATA')
1075 task = self._setup_subtraction(
1076 requiredTemplateFraction=requiredTemplateFraction,
1077 minTemplateFractionForExpectedSuccess=minTemplateFractionForExpectedSuccess
1078 )
1079 if template_coverage < task.config.requiredTemplateFraction:
1080 doRaise = True
1081 elif template_coverage < task.config.minTemplateFractionForExpectedSuccess:
1082 doRaise = True
1083 else:
1084 doRaise = False
1085 if doRaise:
1086 with self.assertRaises(NoWorkFound):
1087 task.run(template_cut, science.clone(), sources.copy(deep=True))
1088 else:
1089 task.run(template_cut, science.clone(), sources.copy(deep=True))
1090 _run_and_check_coverage(template_coverage=0.09)
1091 _run_and_check_coverage(template_coverage=0.15)
1092 _run_and_check_coverage(template_coverage=.7)
1094 def test_clear_template_mask(self):
1095 noiseLevel = 1.
1096 xSize = 400
1097 ySize = 400
1098 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6,
1099 xSize=xSize, ySize=ySize)
1100 template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
1101 templateBorderSize=20, doApplyCalibration=True,
1102 xSize=xSize, ySize=ySize)
1103 diffimEmptyMaskPlanes = ["DETECTED", "DETECTED_NEGATIVE"]
1104 task = self._setup_subtraction()
1105 # Ensure that each each mask plane is set for some pixels
1106 mask = template.mask
1107 x0 = 50
1108 x1 = 75
1109 y0 = 150
1110 y1 = 175
1111 scienceMaskCheck = {}
1112 for maskPlane in mask.getMaskPlaneDict().keys():
1113 scienceMaskCheck[maskPlane] = np.sum(science.mask.array & mask.getPlaneBitMask(maskPlane) > 0)
1114 mask.array[x0: x1, y0: y1] |= mask.getPlaneBitMask(maskPlane)
1115 self.assertTrue(np.sum(mask.array & mask.getPlaneBitMask(maskPlane) > 0))
1117 output = task.run(template, science, sources)
1118 # Verify that the template mask has been modified in place
1119 for maskPlane in mask.getMaskPlaneDict().keys():
1120 if maskPlane in diffimEmptyMaskPlanes:
1121 self.assertTrue(np.sum(mask.array & mask.getPlaneBitMask(maskPlane) == 0))
1122 elif maskPlane in task.config.preserveTemplateMask:
1123 self.assertTrue(np.sum(mask.array & mask.getPlaneBitMask(maskPlane) > 0))
1124 else:
1125 self.assertTrue(np.sum(mask.array & mask.getPlaneBitMask(maskPlane) == 0))
1126 # Mask planes set in the science image should also be set in the difference
1127 # Except the "DETECTED" planes should have been cleared
1128 diffimMask = output.scoreExposure.mask
1129 for maskPlane, scienceSum in scienceMaskCheck.items():
1130 diffimSum = np.sum(diffimMask.array & mask.getPlaneBitMask(maskPlane) > 0)
1131 if maskPlane in diffimEmptyMaskPlanes:
1132 self.assertEqual(diffimSum, 0)
1133 else:
1134 self.assertTrue(diffimSum >= scienceSum)
1136 def test_agnostic_template_psf(self):
1137 """Test that the Score image is the same whether the template PSF is
1138 larger or smaller than the science image PSF.
1139 """
1140 noiseLevel = .3
1141 xSize = 400
1142 ySize = 400
1143 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel,
1144 noiseSeed=6, templateBorderSize=0,
1145 xSize=xSize, ySize=ySize)
1146 template1, _ = makeTestImage(psfSize=self.badPsfSize, noiseLevel=noiseLevel,
1147 noiseSeed=7, doApplyCalibration=True,
1148 xSize=xSize, ySize=ySize)
1149 template2, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel,
1150 noiseSeed=8, doApplyCalibration=True,
1151 xSize=xSize, ySize=ySize)
1152 task = self._setup_subtraction()
1154 science_better = task.run(template1, science.clone(), sources)
1155 template_better = task.run(template2, science, sources)
1156 bbox = science_better.scoreExposure.getBBox().clippedTo(template_better.scoreExposure.getBBox())
1158 delta = template_better.scoreExposure[bbox].clone()
1159 delta.image -= science_better.scoreExposure[bbox].image
1160 delta.variance -= science_better.scoreExposure[bbox].variance
1161 delta.mask.array &= science_better.scoreExposure[bbox].mask.array
1163 statsCtrl = makeStats(badMaskPlanes=("EDGE", "BAD", "NO_DATA"))
1164 # Mean of delta should be very close to zero.
1165 nGoodPix = np.sum(np.isfinite(delta.image.array))
1166 meanError = 2*noiseLevel/np.sqrt(nGoodPix)
1167 deltaMean = computeRobustStatistics(delta.image, delta.mask, statsCtrl)
1168 deltaStd = computeRobustStatistics(delta.image, delta.mask, statsCtrl,
1169 statistic=afwMath.STDEV)
1170 self.assertFloatsAlmostEqual(deltaMean, 0, atol=5*meanError)
1171 nea = computePSFNoiseEquivalentArea(science.psf)
1172 # stddev of Score image should be close to expected value
1173 self.assertFloatsAlmostEqual(deltaStd, np.sqrt(2)*noiseLevel/np.sqrt(nea), rtol=.1)
1175 def test_few_sources(self):
1176 """Test with only 1 source, to check that we get a useful error.
1177 """
1178 xSize = 256
1179 ySize = 256
1180 science, sources = makeTestImage(psfSize=self.midPsfSize, nSrc=10, xSize=xSize, ySize=ySize)
1181 template, _ = makeTestImage(psfSize=self.goodPsfSize, nSrc=10, xSize=xSize, ySize=ySize,
1182 doApplyCalibration=True)
1183 task = self._setup_subtraction()
1184 sources = sources[0:1]
1185 with self.assertRaises(InsufficientKernelSourcesError):
1186 task.run(template, science, sources)
1188 def test_background_subtraction(self):
1189 """Check that we can recover the background,
1190 and that it is subtracted correctly in the Score image.
1191 """
1192 noiseLevel = 1.
1193 xSize = 512
1194 ySize = 512
1195 x0 = 123
1196 y0 = 456
1197 template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
1198 templateBorderSize=20,
1199 xSize=xSize, ySize=ySize, x0=x0, y0=y0,
1200 doApplyCalibration=True)
1201 params = [2.2, 2.1, 2.0, 1.2, 1.1, 1.0]
1203 bbox2D = lsst.geom.Box2D(lsst.geom.Point2D(x0, y0), lsst.geom.Extent2D(xSize, ySize))
1204 background_model = afwMath.Chebyshev1Function2D(params, bbox2D)
1205 science, sources = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=6,
1206 background=background_model,
1207 xSize=xSize, ySize=ySize, x0=x0, y0=y0)
1208 # Don't use ``self._setup_subtraction()`` here.
1209 # Modifying the config of a subtask is messy.
1210 config = subtractImages.AlardLuptonPreconvolveSubtractTask.ConfigClass()
1212 config.sourceSelector.signalToNoise.fluxField = "truth_instFlux"
1213 config.sourceSelector.signalToNoise.errField = "truth_instFluxErr"
1214 config.doSubtractBackground = True
1216 config.makeKernel.kernel.name = "AL"
1217 config.makeKernel.kernel.active.fitForBackground = True
1218 config.makeKernel.kernel.active.spatialKernelOrder = 1
1219 config.makeKernel.kernel.active.spatialBgOrder = 2
1220 statsCtrl = makeStats(badMaskPlanes=("EDGE", "BAD", "NO_DATA"))
1222 task = subtractImages.AlardLuptonPreconvolveSubtractTask(config=config)
1223 output = task.run(template.clone(), science.clone(), sources)
1225 # We should be fitting the same number of parameters as were in the input model
1226 self.assertEqual(output.backgroundModel.getNParameters(), background_model.getNParameters())
1228 # The parameters of the background fit should be close to the input model
1229 self.assertFloatsAlmostEqual(np.array(output.backgroundModel.getParameters()),
1230 np.array(params), rtol=0.2)
1232 # stddev of Score image should be close to expected value.
1233 # This will fail if we have mis-subtracted the background.
1234 stdVal = computeRobustStatistics(output.scoreExposure.image, output.scoreExposure.mask,
1235 statsCtrl, statistic=afwMath.STDEV)
1236 # get the img psf Noise Equivalent Area value
1237 nea = computePSFNoiseEquivalentArea(science.psf)
1238 self.assertFloatsAlmostEqual(stdVal, np.sqrt(2)*noiseLevel/np.sqrt(nea), rtol=0.12)
1240 def test_scale_variance(self):
1241 """Check variance scaling of the Score image.
1242 """
1243 scienceNoiseLevel = 4.
1244 templateNoiseLevel = 2.
1245 scaleFactor = 1.345
1246 xSize = 400
1247 ySize = 400
1248 # Make sure to include pixels with the DETECTED mask bit set.
1249 statsCtrl = makeStats(badMaskPlanes=("EDGE", "BAD", "NO_DATA"))
1251 def _run_and_check_images(science, template, sources, statsCtrl,
1252 doDecorrelation, doScaleVariance, scaleFactor=1.):
1253 """Check that the variance plane matches the expected value for
1254 different configurations of ``doDecorrelation`` and ``doScaleVariance``.
1255 """
1257 task = self._setup_subtraction(doDecorrelation=doDecorrelation,
1258 doScaleVariance=doScaleVariance,
1259 )
1260 output = task.run(template.clone(), science.clone(), sources)
1261 if doScaleVariance:
1262 self.assertFloatsAlmostEqual(task.metadata["scaleTemplateVarianceFactor"],
1263 scaleFactor, atol=0.05)
1264 self.assertFloatsAlmostEqual(task.metadata["scaleScienceVarianceFactor"],
1265 scaleFactor, atol=0.05)
1267 scienceNoise = computeRobustStatistics(science.variance, science.mask, statsCtrl)
1268 # get the img psf Noise Equivalent Area value
1269 nea = computePSFNoiseEquivalentArea(science.psf)
1270 scienceNoise /= nea
1271 if doDecorrelation:
1272 templateNoise = computeRobustStatistics(template.variance, template.mask, statsCtrl)
1273 templateNoise /= nea
1274 else:
1275 # Don't divide by NEA in this case, since the template is convolved
1276 # and in the same units as the Score exposure.
1277 templateNoise = computeRobustStatistics(output.matchedTemplate.variance,
1278 output.matchedTemplate.mask,
1279 statsCtrl)
1280 if doScaleVariance:
1281 templateNoise *= scaleFactor
1282 scienceNoise *= scaleFactor
1283 varMean = computeRobustStatistics(output.scoreExposure.variance,
1284 output.scoreExposure.mask,
1285 statsCtrl)
1286 self.assertFloatsAlmostEqual(varMean, scienceNoise + templateNoise, rtol=0.1)
1288 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=scienceNoiseLevel, noiseSeed=6,
1289 xSize=xSize, ySize=ySize)
1290 template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=templateNoiseLevel, noiseSeed=7,
1291 templateBorderSize=20, doApplyCalibration=True,
1292 xSize=xSize, ySize=ySize)
1293 # Verify that the variance plane of the Score image is correct
1294 # when the template and science variance planes are correct
1295 _run_and_check_images(science, template, sources, statsCtrl,
1296 doDecorrelation=True, doScaleVariance=True)
1297 _run_and_check_images(science, template, sources, statsCtrl,
1298 doDecorrelation=True, doScaleVariance=False)
1299 _run_and_check_images(science, template, sources, statsCtrl,
1300 doDecorrelation=False, doScaleVariance=True)
1301 _run_and_check_images(science, template, sources, statsCtrl,
1302 doDecorrelation=False, doScaleVariance=False)
1304 # Verify that the variance plane of the Score image is correct
1305 # when the template variance plane is incorrect
1306 template.variance.array /= scaleFactor
1307 science.variance.array /= scaleFactor
1308 _run_and_check_images(science, template, sources, statsCtrl,
1309 doDecorrelation=True, doScaleVariance=True, scaleFactor=scaleFactor)
1310 _run_and_check_images(science, template, sources, statsCtrl,
1311 doDecorrelation=True, doScaleVariance=False, scaleFactor=scaleFactor)
1312 _run_and_check_images(science, template, sources, statsCtrl,
1313 doDecorrelation=False, doScaleVariance=True, scaleFactor=scaleFactor)
1314 _run_and_check_images(science, template, sources, statsCtrl,
1315 doDecorrelation=False, doScaleVariance=False, scaleFactor=scaleFactor)
1317 def test_exposure_properties(self):
1318 """Check that all necessary exposure metadata is included
1319 with the Score image.
1320 """
1321 noiseLevel = 1.
1322 xSize = 400
1323 ySize = 400
1324 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6,
1325 xSize=xSize, ySize=ySize)
1326 psf = science.psf
1327 psfAvgPos = psf.getAveragePosition()
1328 psfSize = getPsfFwhm(science.psf)
1329 psfImg = psf.computeKernelImage(psfAvgPos)
1330 template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
1331 templateBorderSize=20, doApplyCalibration=True,
1332 xSize=xSize, ySize=ySize)
1334 def _run_and_check_images(doDecorrelation):
1335 """Check that the metadata is correct with or without decorrelation.
1336 """
1337 task = self._setup_subtraction(doDecorrelation=doDecorrelation)
1338 output = task.run(template.clone(), science.clone(), sources)
1339 psfOut = output.scoreExposure.psf
1340 psfAvgPos = psfOut.getAveragePosition()
1341 if doDecorrelation:
1342 # Decorrelation requires recalculating the PSF,
1343 # so it will not be the same as the input
1344 psfOutSize = getPsfFwhm(science.psf)
1345 self.assertFloatsAlmostEqual(psfSize, psfOutSize)
1346 else:
1347 psfOutImg = psfOut.computeKernelImage(psfAvgPos)
1348 self.assertImagesAlmostEqual(psfImg, psfOutImg)
1350 # check PSF, WCS, bbox, filterLabel, photoCalib
1351 self.assertWcsAlmostEqualOverBBox(science.wcs, output.scoreExposure.wcs, science.getBBox())
1352 self.assertEqual(science.filter, output.scoreExposure.filter)
1353 self.assertEqual(science.photoCalib, output.scoreExposure.photoCalib)
1354 _run_and_check_images(doDecorrelation=True)
1355 _run_and_check_images(doDecorrelation=False)
1358class SimplifiedSubtractTest(AlardLuptonSubtractTestBase, lsst.utils.tests.TestCase):
1359 subtractTask = subtractImages.SimplifiedSubtractTask
1361 def test_runSimplifiedTaskWithExistingKernel(self):
1362 """Test that the simplified task produces the same output as
1363 `AlardLuptonSubtractTask` if it uses the AL kernel.
1364 """
1365 noiseLevel = 1.
1366 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6)
1367 template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
1368 templateBorderSize=20, doApplyCalibration=True)
1369 alTask = AlardLuptonSubtractTest._setup_subtraction(AlardLuptonSubtractTest())
1370 task = self._setup_subtraction(useExistingKernel=True)
1372 alResults = alTask.run(template.clone(), science.clone(), sources)
1373 results = task.run(template.clone(), science.clone(),
1374 inputPsfMatchingKernel=alResults.psfMatchingKernel)
1376 self.assertMaskedImagesEqual(alResults.difference, results.difference)
1378 def test_runSimplifiedTaskWithSourceDetection(self):
1379 """Test that the simplified task with source detection produces
1380 reasonable output.
1381 """
1382 noiseLevel = 1.
1383 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6)
1384 template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
1385 templateBorderSize=20, doApplyCalibration=True)
1386 task = self._setup_subtraction(useExistingKernel=False,
1387 fluxField="base_PsfFlux_instFlux",
1388 errField="base_PsfFlux_instFluxErr",
1389 )
1391 output = task.run(template, science)
1393 # There shoud be no NaN values in the difference image
1394 self.assertTrue(np.all(np.isfinite(output.difference.image.array)))
1395 # Mean of difference image should be close to zero.
1396 meanError = noiseLevel/np.sqrt(output.difference.image.array.size)
1397 # Make sure to include pixels with the DETECTED mask bit set.
1398 statsCtrl = makeStats(badMaskPlanes=("EDGE", "BAD", "NO_DATA", "DETECTED", "DETECTED_NEGATIVE"))
1399 differenceMean = computeRobustStatistics(output.difference.image, output.difference.mask, statsCtrl)
1400 self.assertFloatsAlmostEqual(differenceMean, 0, atol=5*meanError)
1401 # stddev of difference image should be close to expected value.
1402 differenceStd = computeRobustStatistics(output.difference.image, output.difference.mask,
1403 makeStats(), statistic=afwMath.STDEV)
1404 self.assertFloatsAlmostEqual(differenceStd, np.sqrt(2)*noiseLevel, rtol=0.1)
1407def setup_module(module):
1408 lsst.utils.tests.init()
1411class MemoryTestCase(lsst.utils.tests.MemoryTestCase):
1412 pass
1415if __name__ == "__main__": 1415 ↛ 1416line 1415 didn't jump to line 1416 because the condition on line 1415 was never true
1416 lsst.utils.tests.init()
1417 unittest.main()