Coverage for tests/test_subtractTask.py: 99%
804 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-10 09:33 +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.detection as afwDetection
27import lsst.afw.image as afwImage
28import lsst.afw.math as afwMath
29import lsst.afw.table as afwTable
30import lsst.geom
31import lsst.meas.algorithms as measAlg
32from lsst.ip.diffim import subtractImages, InsufficientKernelSourcesError
33from lsst.pex.config import FieldValidationError
34from lsst.pipe.base import NoWorkFound
35import lsst.utils.tests
36import numpy as np
37from lsst.ip.diffim.utils import (computeRobustStatistics, computePSFNoiseEquivalentArea,
38 evaluateMeanPsfFwhm, getPsfFwhm)
39from lsst.pex.exceptions import InvalidParameterError
41from utils import makeStats, makeTestImage, CustomCoaddPsf
44class AlardLuptonSubtractTestBase:
45 goodPsfSize = 2.0
46 midPsfSize = 2.4
47 badPsfSize = 2.8
49 def _setup_subtraction(self, fluxField="truth_instFlux", errField="truth_instFluxErr", **kwargs):
50 """Setup and configure the image subtraction PipelineTask.
52 Parameters
53 ----------
54 fluxField : `str`, optional
55 Name of the flux field in the source catalog.
56 errField : `str`, optional
57 Name of the flux error field in the source catalog.
58 **kwargs
59 Any additional config parameters to set.
61 Returns
62 -------
63 `lsst.pipe.base.PipelineTask`
64 The configured Task to use for detection and measurement.
65 """
66 config = self.subtractTask.ConfigClass()
67 config.doSubtractBackground = False
68 config.restrictKernelEdgeSources = False
69 config.sourceSelector.signalToNoise.fluxField = fluxField
70 config.sourceSelector.signalToNoise.errField = errField
71 config.sourceSelector.doUnresolved = True
72 config.sourceSelector.doIsolated = True
73 config.sourceSelector.doRequirePrimary = True
74 config.sourceSelector.doFlags = True
75 config.sourceSelector.doSignalToNoise = True
76 config.sourceSelector.flags.bad = ["base_PsfFlux_flag", ]
77 config.update(**kwargs)
79 return self.subtractTask(config=config)
82class AlardLuptonSubtractTest(AlardLuptonSubtractTestBase, lsst.utils.tests.TestCase):
83 subtractTask = subtractImages.AlardLuptonSubtractTask
85 def test_allowed_config_modes(self):
86 """Verify the allowable modes for convolution.
87 """
88 config = subtractImages.AlardLuptonSubtractTask.ConfigClass()
89 config.mode = 'auto'
90 config.mode = 'convolveScience'
91 config.mode = 'convolveTemplate'
93 with self.assertRaises(FieldValidationError):
94 config.mode = 'aotu'
96 def test_mismatched_template(self):
97 """Test that an error is raised if the template
98 does not fully contain the science image.
99 """
100 xSize = 200
101 ySize = 200
102 science, sources = makeTestImage(psfSize=self.midPsfSize, xSize=xSize + 20, ySize=ySize + 20)
103 template, _ = makeTestImage(psfSize=self.midPsfSize, xSize=xSize, ySize=ySize,
104 doApplyCalibration=True)
105 task = self._setup_subtraction()
106 with self.assertRaises(AssertionError):
107 task.run(template, science, sources)
109 def test_mismatched_filter(self):
110 """Test that an error is raised if the science and template have
111 different bands.
112 """
113 xSize = 200
114 ySize = 200
115 science, sources = makeTestImage(psfSize=self.midPsfSize, xSize=xSize + 20, ySize=ySize + 20,
116 band="g", physicalFilter="g noCamera")
117 template, _ = makeTestImage(psfSize=self.midPsfSize, xSize=xSize, ySize=ySize,
118 doApplyCalibration=True, band="not-g", physicalFilter="not-g noCamera")
119 task = self._setup_subtraction()
120 with self.assertRaises(AssertionError):
121 task.run(template, science, sources)
123 def test_incomplete_template_coverage(self):
124 noiseLevel = 1.
125 border = 20
126 xSize = 400
127 ySize = 400
128 nSources = 80
129 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6,
130 nSrc=nSources, xSize=xSize, ySize=ySize)
131 template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
132 nSrc=nSources, templateBorderSize=border, doApplyCalibration=True,
133 xSize=xSize, ySize=ySize)
135 science_height = science.getBBox().getDimensions().getY()
137 def _run_and_check_coverage(template_coverage,
138 requiredTemplateFraction=0.1,
139 minTemplateFractionForExpectedSuccess=0.2):
140 template_cut = template.clone()
141 template_height = int(science_height*template_coverage + border)
142 template_cut.image.array[:, template_height:] = 0
143 template_cut.mask.array[:, template_height:] = template_cut.mask.getPlaneBitMask('NO_DATA')
144 task = self._setup_subtraction(
145 requiredTemplateFraction=requiredTemplateFraction,
146 minTemplateFractionForExpectedSuccess=minTemplateFractionForExpectedSuccess
147 )
148 if template_coverage < requiredTemplateFraction:
149 doRaise = True
150 elif template_coverage < minTemplateFractionForExpectedSuccess:
151 doRaise = True
152 else:
153 doRaise = False
154 if doRaise:
155 with self.assertRaises(NoWorkFound):
156 task.run(template_cut, science.clone(), sources.copy(deep=True))
157 else:
158 task.run(template_cut, science.clone(), sources.copy(deep=True))
159 _run_and_check_coverage(template_coverage=0.09)
160 _run_and_check_coverage(template_coverage=0.15)
161 _run_and_check_coverage(template_coverage=0.7)
163 def test_clear_template_mask(self):
164 noiseLevel = 1.
165 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6)
166 template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
167 templateBorderSize=20, doApplyCalibration=True)
168 diffimEmptyMaskPlanes = ["DETECTED", "DETECTED_NEGATIVE"]
169 task = self._setup_subtraction(mode="convolveTemplate")
170 # Ensure that each each mask plane is set for some pixels
171 mask = template.mask
172 x0 = 50
173 x1 = 75
174 y0 = 150
175 y1 = 175
176 scienceMaskCheck = {}
177 for maskPlane in mask.getMaskPlaneDict().keys():
178 scienceMaskCheck[maskPlane] = np.sum(science.mask.array & mask.getPlaneBitMask(maskPlane) > 0)
179 mask.array[x0: x1, y0: y1] |= mask.getPlaneBitMask(maskPlane)
180 self.assertTrue(np.sum(mask.array & mask.getPlaneBitMask(maskPlane) > 0))
182 output = task.run(template, science, sources)
183 # Verify that the template mask has been modified in place
184 for maskPlane in mask.getMaskPlaneDict().keys():
185 if maskPlane in diffimEmptyMaskPlanes:
186 self.assertTrue(np.sum(mask.array & mask.getPlaneBitMask(maskPlane) == 0))
187 elif maskPlane in task.config.preserveTemplateMask:
188 self.assertTrue(np.sum(mask.array & mask.getPlaneBitMask(maskPlane) > 0))
189 else:
190 self.assertTrue(np.sum(mask.array & mask.getPlaneBitMask(maskPlane) == 0))
191 # Mask planes set in the science image should also be set in the difference
192 # Except the "DETECTED" planes should have been cleared
193 diffimMask = output.difference.mask
194 for maskPlane, scienceSum in scienceMaskCheck.items():
195 diffimSum = np.sum(diffimMask.array & mask.getPlaneBitMask(maskPlane) > 0)
196 if maskPlane in diffimEmptyMaskPlanes:
197 self.assertEqual(diffimSum, 0)
198 else:
199 self.assertTrue(diffimSum >= scienceSum)
201 def test_equal_images(self):
202 """Test that running with enough sources produces reasonable output,
203 with the same size psf in the template and science.
204 """
205 noiseLevel = 1.
206 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6)
207 template, _ = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
208 templateBorderSize=20, doApplyCalibration=True)
209 task = self._setup_subtraction()
210 output = task.run(template, science, sources)
211 # There shoud be no NaN values in the difference image
212 self.assertTrue(np.all(np.isfinite(output.difference.image.array)))
213 # Mean of difference image should be close to zero.
214 meanError = noiseLevel/np.sqrt(output.difference.image.array.size)
215 # Make sure to include pixels with the DETECTED mask bit set.
216 statsCtrl = makeStats(badMaskPlanes=("EDGE", "BAD", "NO_DATA", "DETECTED", "DETECTED_NEGATIVE"))
217 differenceMean = computeRobustStatistics(output.difference.image, output.difference.mask, statsCtrl)
218 self.assertFloatsAlmostEqual(differenceMean, 0, atol=5*meanError)
219 # stddev of difference image should be close to expected value.
220 differenceStd = computeRobustStatistics(output.difference.image, output.difference.mask,
221 makeStats(), statistic=afwMath.STDEV)
222 self.assertFloatsAlmostEqual(differenceStd, np.sqrt(2)*noiseLevel, rtol=0.1)
224 def test_equal_images_missing_mask_planes(self):
225 """Test that running with enough sources produces reasonable output,
226 with the same size psf in the template and science and with missing
227 mask planes.
228 """
229 noiseLevel = 1.
230 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6,
231 addMaskPlanes=[])
232 template, _ = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
233 templateBorderSize=20, doApplyCalibration=True, addMaskPlanes=[])
234 task = self._setup_subtraction()
235 output = task.run(template, science, sources)
236 # There shoud be no NaN values in the difference image
237 self.assertTrue(np.all(np.isfinite(output.difference.image.array)))
238 # Mean of difference image should be close to zero.
239 meanError = noiseLevel/np.sqrt(output.difference.image.array.size)
240 # Make sure to include pixels with the DETECTED mask bit set.
241 statsCtrl = makeStats(badMaskPlanes=("EDGE", "BAD", "NO_DATA", "DETECTED", "DETECTED_NEGATIVE"))
242 differenceMean = computeRobustStatistics(output.difference.image, output.difference.mask, statsCtrl)
243 self.assertFloatsAlmostEqual(differenceMean, 0, atol=5*meanError)
244 # stddev of difference image should be close to expected value.
245 differenceStd = computeRobustStatistics(output.difference.image, output.difference.mask,
246 makeStats(), statistic=afwMath.STDEV)
247 self.assertFloatsAlmostEqual(differenceStd, np.sqrt(2)*noiseLevel, rtol=0.1)
249 def test_psf_size(self):
250 """Test that the image subtract task runs without failing, if
251 fwhmExposureBuffer and fwhmExposureGrid parameters are set.
252 """
253 noiseLevel = 1.
254 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6)
255 template, _ = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
256 templateBorderSize=20, doApplyCalibration=True)
258 schema = afwTable.ExposureTable.makeMinimalSchema()
259 weightKey = schema.addField("weight", type="D", doc="Coadd weight")
260 exposureCatalog = afwTable.ExposureCatalog(schema)
261 kernel = measAlg.DoubleGaussianPsf(7, 7, 2.0).getKernel()
262 psf = measAlg.KernelPsf(kernel, template.getBBox().getCenter())
264 record = exposureCatalog.addNew()
265 record.setPsf(psf)
266 record.setWcs(template.wcs)
267 record.setD(weightKey, 1.0)
268 record.setBBox(template.getBBox())
270 customPsf = CustomCoaddPsf(exposureCatalog, template.wcs)
271 template.setPsf(customPsf)
273 # Test that we get an exception if we simply get the FWHM at center.
274 with self.assertRaises(InvalidParameterError):
275 getPsfFwhm(template.psf, True)
277 with self.assertRaises(InvalidParameterError):
278 getPsfFwhm(template.psf, False)
280 # Test that evaluateMeanPsfFwhm runs successfully on the template.
281 evaluateMeanPsfFwhm(template, fwhmExposureBuffer=0.05, fwhmExposureGrid=10)
283 # Since the PSF is spatially invariant, the FWHM should be the same at
284 # all points in the science image.
285 fwhm1 = getPsfFwhm(science.psf, False)
286 fwhm2 = evaluateMeanPsfFwhm(science, fwhmExposureBuffer=0.05, fwhmExposureGrid=10)
287 self.assertAlmostEqual(fwhm1[0], fwhm2, places=13)
288 self.assertAlmostEqual(fwhm1[1], fwhm2, places=13)
290 self.assertAlmostEqual(evaluateMeanPsfFwhm(science, fwhmExposureBuffer=0.05,
291 fwhmExposureGrid=10),
292 getPsfFwhm(science.psf, True), places=7
293 )
295 # Test that the image subtraction task runs successfully.
296 task = self._setup_subtraction()
298 # Test that the task runs if we take the mean FWHM on a grid.
299 with self.assertLogs(level="INFO") as cm:
300 task.run(template, science, sources)
302 # Check that evaluateMeanPsfFwhm was called.
303 # This tests that getPsfFwhm failed raising InvalidParameterError,
304 # that is caught and handled appropriately.
305 logMessage = ("INFO:lsst.alardLuptonSubtract:Unable to evaluate PSF at the average position. "
306 "Evaluting PSF on a grid of points."
307 )
308 self.assertIn(logMessage, cm.output)
310 def test_auto_convolveTemplate(self):
311 """Test that auto mode gives the same result as convolveTemplate when
312 the template psf is the smaller.
313 """
314 noiseLevel = 1.
315 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6)
316 template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
317 templateBorderSize=20, doApplyCalibration=True)
318 task = self._setup_subtraction(mode="convolveTemplate")
319 output = task.run(template.clone(), science.clone(), sources)
321 task = self._setup_subtraction(mode="auto")
322 outputAuto = task.run(template, science, sources)
323 self.assertMaskedImagesEqual(output.difference.maskedImage, outputAuto.difference.maskedImage)
325 def test_auto_convolveScience(self):
326 """Test that auto mode gives the same result as convolveScience when
327 the science psf is the smaller.
328 """
329 noiseLevel = 1.
330 science, sources = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=6)
331 template, _ = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
332 templateBorderSize=20, doApplyCalibration=True)
333 task = self._setup_subtraction(mode="convolveScience")
334 output = task.run(template.clone(), science.clone(), sources)
336 task = self._setup_subtraction(mode="auto")
337 outputAuto = task.run(template, science, sources)
338 self.assertMaskedImagesEqual(output.difference.maskedImage, outputAuto.difference.maskedImage)
340 def test_science_better(self):
341 """Test that running with enough sources produces reasonable output,
342 with the science psf being smaller than the template.
343 """
344 statsCtrl = makeStats()
345 statsCtrlDetect = makeStats(badMaskPlanes=("EDGE", "BAD", "NO_DATA"))
347 def _run_and_check_images(statsCtrl, statsCtrlDetect, scienceNoiseLevel, templateNoiseLevel):
348 science, sources = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=scienceNoiseLevel,
349 noiseSeed=6)
350 template, _ = makeTestImage(psfSize=self.midPsfSize, noiseLevel=templateNoiseLevel, noiseSeed=7,
351 templateBorderSize=20, doApplyCalibration=True)
352 task = self._setup_subtraction(mode="convolveScience")
353 output = task.run(template, science, sources)
354 self.assertFloatsAlmostEqual(task.metadata["scaleScienceVarianceFactor"], 1., atol=.05)
355 # Mean of difference image should be close to zero.
356 nGoodPix = np.sum(np.isfinite(output.difference.image.array))
357 meanError = (scienceNoiseLevel + templateNoiseLevel)/np.sqrt(nGoodPix)
358 diffimMean = computeRobustStatistics(output.difference.image, output.difference.mask,
359 statsCtrlDetect)
361 self.assertFloatsAlmostEqual(diffimMean, 0, atol=5*meanError)
362 # stddev of difference image should be close to expected value.
363 noiseLevel = np.sqrt(scienceNoiseLevel**2 + templateNoiseLevel**2)
364 varianceMean = computeRobustStatistics(output.difference.variance, output.difference.mask,
365 statsCtrl)
366 diffimStd = computeRobustStatistics(output.difference.image, output.difference.mask,
367 statsCtrl, statistic=afwMath.STDEV)
368 self.assertFloatsAlmostEqual(varianceMean, noiseLevel**2, rtol=0.1)
369 self.assertFloatsAlmostEqual(diffimStd, noiseLevel, rtol=0.1)
371 _run_and_check_images(statsCtrl, statsCtrlDetect, scienceNoiseLevel=1., templateNoiseLevel=1.)
372 _run_and_check_images(statsCtrl, statsCtrlDetect, scienceNoiseLevel=1., templateNoiseLevel=.1)
373 _run_and_check_images(statsCtrl, statsCtrlDetect, scienceNoiseLevel=.1, templateNoiseLevel=.1)
375 def test_template_better(self):
376 """Test that running with enough sources produces reasonable output,
377 with the template psf being smaller than the science.
378 """
379 statsCtrl = makeStats()
380 statsCtrlDetect = makeStats(badMaskPlanes=("EDGE", "BAD", "NO_DATA"))
382 def _run_and_check_images(statsCtrl, statsCtrlDetect, scienceNoiseLevel, templateNoiseLevel):
383 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=scienceNoiseLevel,
384 noiseSeed=6)
385 template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=templateNoiseLevel, noiseSeed=7,
386 templateBorderSize=20, doApplyCalibration=True)
387 task = self._setup_subtraction()
388 output = task.run(template, science, sources)
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["scaleScienceVarianceFactor"],
636 scaleFactor, atol=0.05)
638 scienceNoise = computeRobustStatistics(science.variance, science.mask, statsCtrl)
639 if doDecorrelation:
640 templateNoise = computeRobustStatistics(template.variance, template.mask, statsCtrl)
641 else:
642 templateNoise = computeRobustStatistics(output.matchedTemplate.variance,
643 output.matchedTemplate.mask,
644 statsCtrl)
646 if doScaleVariance:
647 # Only the science variance is scaled here. The template
648 # variance is scaled independently in ``GetTemplateTask``.
649 scienceNoise *= scaleFactor
650 varMean = computeRobustStatistics(output.difference.variance, output.difference.mask, statsCtrl)
651 self.assertFloatsAlmostEqual(varMean, scienceNoise + templateNoise, rtol=0.1)
653 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=scienceNoiseLevel, noiseSeed=6)
654 template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=templateNoiseLevel, noiseSeed=7,
655 templateBorderSize=20, doApplyCalibration=True)
656 # Verify that the variance plane of the difference image is correct
657 # when the template and science variance planes are correct
658 _run_and_check_images(science, template, sources, statsCtrl,
659 doDecorrelation=True, doScaleVariance=True)
660 _run_and_check_images(science, template, sources, statsCtrl,
661 doDecorrelation=True, doScaleVariance=False)
662 _run_and_check_images(science, template, sources, statsCtrl,
663 doDecorrelation=False, doScaleVariance=True)
664 _run_and_check_images(science, template, sources, statsCtrl,
665 doDecorrelation=False, doScaleVariance=False)
667 # Verify that the variance plane of the difference image is correct
668 # when the input science variance plane is incorrect
669 science.variance.array /= scaleFactor
670 _run_and_check_images(science, template, sources, statsCtrl,
671 doDecorrelation=True, doScaleVariance=True, scaleFactor=scaleFactor)
672 _run_and_check_images(science, template, sources, statsCtrl,
673 doDecorrelation=True, doScaleVariance=False, scaleFactor=scaleFactor)
674 _run_and_check_images(science, template, sources, statsCtrl,
675 doDecorrelation=False, doScaleVariance=True, scaleFactor=scaleFactor)
676 _run_and_check_images(science, template, sources, statsCtrl,
677 doDecorrelation=False, doScaleVariance=False, scaleFactor=scaleFactor)
679 def test_scale_variance_convolve_science(self):
680 """Check variance scaling of the image difference.
681 """
682 scienceNoiseLevel = 4.
683 templateNoiseLevel = 2.
684 scaleFactor = 1.345
685 # Make sure to include pixels with the DETECTED mask bit set.
686 statsCtrl = makeStats(badMaskPlanes=("EDGE", "BAD", "NO_DATA"))
688 def _run_and_check_images(science, template, sources, statsCtrl,
689 doDecorrelation, doScaleVariance, scaleFactor=1.):
690 """Check that the variance plane matches the expected value for
691 different configurations of ``doDecorrelation`` and ``doScaleVariance``.
692 """
694 task = self._setup_subtraction(mode="convolveScience",
695 doDecorrelation=doDecorrelation,
696 doScaleVariance=doScaleVariance,
697 restrictKernelEdgeSources=False,
698 )
699 output = task.run(template.clone(), science.clone(), sources)
700 if doScaleVariance:
701 self.assertFloatsAlmostEqual(task.metadata["scaleScienceVarianceFactor"],
702 scaleFactor, atol=0.05)
704 templateNoise = computeRobustStatistics(template.variance, template.mask, statsCtrl)
705 if doDecorrelation:
706 scienceNoise = computeRobustStatistics(science.variance, science.mask, statsCtrl)
707 else:
708 scienceNoise = computeRobustStatistics(output.matchedScience.variance,
709 output.matchedScience.mask,
710 statsCtrl)
712 if doScaleVariance:
713 # Only the science variance is scaled here. The template
714 # variance is scaled independently in ``GetTemplateTask``.
715 scienceNoise *= scaleFactor
717 varMean = computeRobustStatistics(output.difference.variance, output.difference.mask, statsCtrl)
718 self.assertFloatsAlmostEqual(varMean, scienceNoise + templateNoise, rtol=0.1)
720 science, sources = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=scienceNoiseLevel, noiseSeed=6)
721 template, _ = makeTestImage(psfSize=self.midPsfSize, noiseLevel=templateNoiseLevel, noiseSeed=7,
722 templateBorderSize=20, doApplyCalibration=True)
723 # Verify that the variance plane of the difference image is correct
724 # when the template and science variance planes are correct
725 _run_and_check_images(science, template, sources, statsCtrl,
726 doDecorrelation=True, doScaleVariance=True)
727 _run_and_check_images(science, template, sources, statsCtrl,
728 doDecorrelation=True, doScaleVariance=False)
729 _run_and_check_images(science, template, sources, statsCtrl,
730 doDecorrelation=False, doScaleVariance=True)
731 _run_and_check_images(science, template, sources, statsCtrl,
732 doDecorrelation=False, doScaleVariance=False)
734 # Verify that the variance plane of the difference image is correct
735 # when the input science variance plane is incorrect
736 science.variance.array /= scaleFactor
737 _run_and_check_images(science, template, sources, statsCtrl,
738 doDecorrelation=True, doScaleVariance=True, scaleFactor=scaleFactor)
739 _run_and_check_images(science, template, sources, statsCtrl,
740 doDecorrelation=True, doScaleVariance=False, scaleFactor=scaleFactor)
741 _run_and_check_images(science, template, sources, statsCtrl,
742 doDecorrelation=False, doScaleVariance=True, scaleFactor=scaleFactor)
743 _run_and_check_images(science, template, sources, statsCtrl,
744 doDecorrelation=False, doScaleVariance=False, scaleFactor=scaleFactor)
746 def test_exposure_properties_convolve_template(self):
747 """Check that all necessary exposure metadata is included
748 when the template is convolved.
749 """
750 noiseLevel = 1.
751 seed = 37
752 rng = np.random.RandomState(seed)
753 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6)
754 psf = science.psf
755 psfAvgPos = psf.getAveragePosition()
756 psfSize = getPsfFwhm(science.psf)
757 psfImg = psf.computeKernelImage(psfAvgPos)
758 template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
759 templateBorderSize=20, doApplyCalibration=True)
761 # Generate a random aperture correction map
762 apCorrMap = lsst.afw.image.ApCorrMap()
763 for name in ("a", "b", "c"):
764 apCorrMap.set(name, lsst.afw.math.ChebyshevBoundedField(science.getBBox(), rng.randn(3, 3)))
765 science.info.setApCorrMap(apCorrMap)
767 def _run_and_check_images(doDecorrelation):
768 """Check that the metadata is correct with or without decorrelation.
769 """
770 task = self._setup_subtraction(mode="convolveTemplate",
771 doDecorrelation=doDecorrelation,
772 )
773 output = task.run(template.clone(), science.clone(), sources)
774 psfOut = output.difference.psf
775 psfAvgPos = psfOut.getAveragePosition()
776 if doDecorrelation:
777 # Decorrelation requires recalculating the PSF,
778 # so it will not be the same as the input
779 psfOutSize = getPsfFwhm(science.psf)
780 self.assertFloatsAlmostEqual(psfSize, psfOutSize)
781 else:
782 psfOutImg = psfOut.computeKernelImage(psfAvgPos)
783 self.assertImagesAlmostEqual(psfImg, psfOutImg)
785 # check PSF, WCS, bbox, filterLabel, photoCalib, aperture correction
786 self._compare_apCorrMaps(apCorrMap, output.difference.info.getApCorrMap())
787 self.assertWcsAlmostEqualOverBBox(science.wcs, output.difference.wcs, science.getBBox())
788 self.assertEqual(science.filter, output.difference.filter)
789 self.assertEqual(science.photoCalib, output.difference.photoCalib)
790 _run_and_check_images(doDecorrelation=True)
791 _run_and_check_images(doDecorrelation=False)
793 def test_exposure_properties_convolve_science(self):
794 """Check that all necessary exposure metadata is included
795 when the science image is convolved.
796 """
797 noiseLevel = 1.
798 seed = 37
799 rng = np.random.RandomState(seed)
800 science, sources = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=6)
801 template, _ = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
802 templateBorderSize=20, doApplyCalibration=True)
803 psf = template.psf
804 psfAvgPos = psf.getAveragePosition()
805 psfSize = getPsfFwhm(template.psf)
806 psfImg = psf.computeKernelImage(psfAvgPos)
808 # Generate a random aperture correction map
809 apCorrMap = lsst.afw.image.ApCorrMap()
810 for name in ("a", "b", "c"):
811 apCorrMap.set(name, lsst.afw.math.ChebyshevBoundedField(science.getBBox(), rng.randn(3, 3)))
812 science.info.setApCorrMap(apCorrMap)
814 def _run_and_check_images(doDecorrelation):
815 """Check that the metadata is correct with or without decorrelation.
816 """
817 task = self._setup_subtraction(mode="convolveScience",
818 doDecorrelation=doDecorrelation,
819 )
820 output = task.run(template.clone(), science.clone(), sources)
821 if doDecorrelation:
822 # Decorrelation requires recalculating the PSF,
823 # so it will not be the same as the input
824 psfOutSize = getPsfFwhm(template.psf)
825 self.assertFloatsAlmostEqual(psfSize, psfOutSize)
826 else:
827 psfOut = output.difference.psf
828 psfAvgPos = psfOut.getAveragePosition()
829 psfOutImg = psfOut.computeKernelImage(psfAvgPos)
830 self.assertImagesAlmostEqual(psfImg, psfOutImg)
832 # check PSF, WCS, bbox, filterLabel, photoCalib, aperture correction
833 self._compare_apCorrMaps(apCorrMap, output.difference.info.getApCorrMap())
834 self.assertWcsAlmostEqualOverBBox(science.wcs, output.difference.wcs, science.getBBox())
835 self.assertEqual(science.filter, output.difference.filter)
836 self.assertEqual(science.photoCalib, output.difference.photoCalib)
838 _run_and_check_images(doDecorrelation=True)
839 _run_and_check_images(doDecorrelation=False)
841 def _compare_apCorrMaps(self, a, b):
842 """Compare two ApCorrMaps for equality, without assuming that their BoundedFields have the
843 same addresses (i.e. so we can compare after serialization).
845 This function is taken from ``ApCorrMapTestCase`` in afw/tests/.
847 Parameters
848 ----------
849 a, b : `lsst.afw.image.ApCorrMap`
850 The two aperture correction maps to compare.
851 """
852 self.assertEqual(len(a), len(b))
853 for name, value in list(a.items()):
854 value2 = b.get(name)
855 self.assertIsNotNone(value2)
856 self.assertEqual(value.getBBox(), value2.getBBox())
857 self.assertFloatsAlmostEqual(
858 value.getCoefficients(), value2.getCoefficients(), rtol=0.0)
860 def test_fake_mask_plane_propagation(self):
861 """Test that we have the mask planes related to fakes in diffim images.
862 This is testing method called updateMasks
863 """
864 xSize = 200
865 ySize = 200
866 science, sources = makeTestImage(psfSize=self.midPsfSize, xSize=xSize, ySize=ySize)
867 science_fake_img, science_fake_sources = makeTestImage(
868 psfSize=self.midPsfSize, xSize=xSize, ySize=ySize, seed=7, nSrc=2, noiseLevel=0.25, fluxRange=1
869 )
870 template, _ = makeTestImage(psfSize=self.midPsfSize, xSize=xSize, ySize=ySize,
871 doApplyCalibration=True)
872 tmplt_fake_img, tmplt_fake_sources = makeTestImage(
873 psfSize=self.midPsfSize, xSize=xSize, ySize=ySize, seed=9, nSrc=2, noiseLevel=0.25, fluxRange=1
874 )
875 # created fakes and added them to the images
876 science.image.array += science_fake_img.image.array
877 template.image.array += tmplt_fake_img.image.array
879 # TODO: DM-40796 update to INJECTED names when source injection gets refactored
880 # adding mask planes to both science and template images
881 science_mask_planes = science.mask.addMaskPlane("FAKE")
882 template_mask_planes = template.mask.addMaskPlane("FAKE")
884 for a_science_source in science_fake_sources:
885 # 3 x 3 masking of the source locations is fine
886 bbox = lsst.geom.Box2I(
887 lsst.geom.Point2I(a_science_source.getX(), a_science_source.getY()), lsst.geom.Extent2I(3, 3)
888 )
889 science[bbox].mask.array |= science_mask_planes
891 for a_template_source in tmplt_fake_sources:
892 # 3 x 3 masking of the source locations is fine
893 bbox = lsst.geom.Box2I(
894 lsst.geom.Point2I(a_template_source.getX(), a_template_source.getY()),
895 lsst.geom.Extent2I(3, 3)
896 )
897 template[bbox].mask.array |= template_mask_planes
899 science_fake_masked = (science.mask.array & science.mask.getPlaneBitMask("FAKE")) > 0
900 template_fake_masked = (template.mask.array & template.mask.getPlaneBitMask("FAKE")) > 0
902 task = self._setup_subtraction()
903 subtraction = task.run(template, science, sources)
905 # check subtraction mask plane is set where we set the previous masks
906 diff_mask = subtraction.difference.mask
908 # science mask should be now in INJECTED
909 inj_masked = (diff_mask.array & diff_mask.getPlaneBitMask("INJECTED")) > 0
911 # template mask should be now in INJECTED_TEMPLATE
912 injTmplt_masked = (diff_mask.array & diff_mask.getPlaneBitMask("INJECTED_TEMPLATE")) > 0
914 self.assertEqual(np.sum(inj_masked.astype(int)-science_fake_masked.astype(int)), 0)
915 self.assertEqual(np.sum(injTmplt_masked.astype(int)-template_fake_masked.astype(int)), 0)
917 def test_calculateMagLim(self):
918 """Test that the limiting magnitude matches the analytic value for a
919 Gaussian PSF and flat noise, and does not depend on the pixel units
920 or on masked pixels.
921 """
922 sigma = 1.5
923 psfSigma = 2.0
924 calibration = 3.0
925 bbox = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Extent2I(200, 200))
926 exposure = afwImage.ExposureF(bbox)
927 exposure.variance.array[:, :] = sigma**2
928 exposure.setPsf(afwDetection.GaussianPsf(41, 41, psfSigma))
929 exposure.setPhotoCalib(afwImage.PhotoCalib(calibration))
930 task = self._setup_subtraction()
932 maglim = task._calculateMagLim(exposure)
933 fluxLim = 5*sigma*np.sqrt(4*np.pi*psfSigma**2)
934 self.assertFloatsAlmostEqual(maglim, afwImage.PhotoCalib(calibration).instFluxToMagnitude(fluxLim),
935 atol=1e-3)
937 # Pixels in other units, with a matching calibration.
938 scaled = exposure.clone()
939 scaled.maskedImage *= 10.0
940 scaled.setPhotoCalib(afwImage.PhotoCalib(calibration/10.0))
941 self.assertFloatsAlmostEqual(task._calculateMagLim(scaled), maglim, atol=1e-6)
943 # Masked pixels do not count toward the noise.
944 masked = exposure.clone()
945 masked.variance.array[:, :100] = 1e6
946 masked.mask.array[:, :100] |= masked.mask.getPlaneBitMask("NO_DATA")
947 self.assertFloatsAlmostEqual(task._calculateMagLim(masked), maglim, atol=1e-6)
949 masked.mask.array[:, :] |= masked.mask.getPlaneBitMask("NO_DATA")
950 self.assertTrue(np.isnan(task._calculateMagLim(masked)))
952 def test_metadata_metrics(self):
953 """Verify fields are added to metadata when subtraction is run, and
954 that the difference image limiting magnitude is calculated correctly,
955 both with a "good" and "bad" seeing template.
956 """
957 science, sources = makeTestImage(psfSize=self.badPsfSize, noiseLevel=1)
958 template_good, _ = makeTestImage(psfSize=self.midPsfSize, doApplyCalibration=True, noiseLevel=0.25,
959 templateBorderSize=20)
960 template_bad, _ = makeTestImage(psfSize=9.5, doApplyCalibration=True, noiseLevel=0.25,
961 templateBorderSize=20)
963 # Add a few sky objects; sky footprints are needed for some metrics.
964 config = measAlg.SkyObjectsTask.ConfigClass()
965 config.nSources = 3
966 skyTask = measAlg.SkyObjectsTask(config=config, name="skySources")
967 skyTask.skySourceKey = sources.schema["sky_source"].asKey()
968 skyTask.run(science.mask, 10, catalog=sources)
969 sources = sources.copy(deep=True)
970 # Add centroids, since these sources were added post-measurement.
971 for record in sources[sources["sky_source"]]:
972 record["truth_x"] = record.getFootprint().getPeaks()[0].getFx()
973 record["truth_y"] = record.getFootprint().getPeaks()[0].getFy()
975 # The metadata fields are attached to the subtractTask, so we do
976 # need to run that; run it for both "good" and "bad" seeing templates
978 subtractTask_good = self._setup_subtraction()
979 _ = subtractTask_good.run(template_good.clone(), science.clone(), sources)
980 subtractTask_bad = self._setup_subtraction()
981 _ = subtractTask_bad.run(template_bad.clone(), science.clone(), sources)
983 # Test that the diffim limiting magnitude combines the science and
984 # template limits in flux. The limits depend on the variance and mask
985 # planes, which `run` modifies, so take them from the metadata.
986 def combineMagLims(metadata):
987 fluxlim_science = (metadata['scienceLimitingMagnitude']*u.ABmag).to_value(u.nJy)
988 fluxlim_template = (metadata['templateLimitingMagnitude']*u.ABmag).to_value(u.nJy)
989 return (np.hypot(fluxlim_science, fluxlim_template)*u.nJy).to(u.ABmag).value
991 for task in (subtractTask_good, subtractTask_bad):
992 self.assertFloatsAlmostEqual(task.metadata['diffimLimitingMagnitude'],
993 combineMagLims(task.metadata), atol=1e-6)
994 # The template with the wider PSF is shallower.
995 self.assertLess(subtractTask_bad.metadata['templateLimitingMagnitude'],
996 subtractTask_good.metadata['templateLimitingMagnitude'])
998 # Create a template with a PSF that is not defined at the image center.
999 # First, make an exposure catalog so we can force the template to have
1000 # a bad (off-image) PSF. It must have a record with a weight field
1001 # and a BBox in order to let us set the PSF manually.
1002 template_offimage, _ = makeTestImage()
1003 schema = afwTable.ExposureTable.makeMinimalSchema()
1004 weightKey = schema.addField("weight", type="D", doc="Coadd weight")
1005 exposureCatalog = afwTable.ExposureCatalog(schema)
1006 record = exposureCatalog.addNew()
1007 record.setD(weightKey, 1.0)
1008 record.setBBox(template_offimage.getBBox())
1009 kernel = measAlg.DoubleGaussianPsf(7, 7, 2.0).getKernel()
1010 psf = measAlg.KernelPsf(kernel, template_offimage.getBBox().getCenter())
1011 record.setPsf(psf)
1012 record.setWcs(template_offimage.wcs)
1013 custom_offimage_psf = CustomCoaddPsf(exposureCatalog, template_offimage.wcs)
1014 template_offimage.setPsf(custom_offimage_psf)
1016 # Test that template PSF size edge cases are handled correctly.
1017 subtractTask_offimage = self._setup_subtraction()
1018 _ = subtractTask_offimage.run(template_offimage.clone(), science.clone(), sources)
1019 # Test that providing no fallbackPsfSize results in a nan template
1020 # limiting magnitude.
1021 maglim_template_offimage = subtractTask_offimage._calculateMagLim(template_offimage)
1022 self.assertTrue(np.isnan(maglim_template_offimage))
1023 # Test that the task falls back to the provided fallbackPsfSize, so
1024 # the template and diffim limiting magnitudes are finite.
1025 self.assertTrue(np.isfinite(subtractTask_offimage.metadata['templateLimitingMagnitude']))
1026 self.assertFloatsAlmostEqual(subtractTask_offimage.metadata['diffimLimitingMagnitude'],
1027 combineMagLims(subtractTask_offimage.metadata), atol=1e-6)
1029 # Test that several other expected metadata metrics exist
1030 self.assertIn('scienceLimitingMagnitude', subtractTask_good.metadata)
1031 self.assertIn('templateLimitingMagnitude', subtractTask_good.metadata)
1033 # The mean ratio metric should be much worse on the "bad" subtraction.
1034 self.assertLess(subtractTask_good.metadata['differenceFootprintRatioMean'], 0.02)
1035 self.assertGreater(subtractTask_bad.metadata['differenceFootprintRatioMean'], 0.12)
1038class AlardLuptonPreconvolveSubtractTest(AlardLuptonSubtractTestBase, lsst.utils.tests.TestCase):
1039 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask
1041 def test_mismatched_template(self):
1042 """Test that an error is raised if the template
1043 does not fully contain the science image.
1044 """
1045 xSize = 200
1046 ySize = 200
1047 science, sources = makeTestImage(psfSize=self.midPsfSize, xSize=xSize + 20, ySize=ySize + 20)
1048 template, _ = makeTestImage(psfSize=self.midPsfSize, xSize=xSize, ySize=ySize,
1049 doApplyCalibration=True)
1050 task = self._setup_subtraction()
1051 with self.assertRaises(AssertionError):
1052 task.run(template, science, sources)
1054 def test_equal_images(self):
1055 """Test that running with enough sources produces reasonable output,
1056 with the same size psf in the template and science.
1057 """
1058 noiseLevel = 1.
1059 xSize = 400
1060 ySize = 400
1061 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6,
1062 xSize=xSize, ySize=ySize)
1063 template, _ = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
1064 templateBorderSize=20, doApplyCalibration=True,
1065 xSize=xSize, ySize=ySize)
1066 task = self._setup_subtraction()
1067 output = task.run(template, science, sources)
1068 # There shoud be no NaN values in the Score image
1069 self.assertTrue(np.all(np.isfinite(output.scoreExposure.image.array)))
1070 # Mean of Score image should be close to zero.
1071 meanError = noiseLevel/np.sqrt(output.scoreExposure.image.array.size)
1072 # Make sure to include pixels with the DETECTED mask bit set.
1073 statsCtrl = makeStats(badMaskPlanes=("EDGE", "BAD", "NO_DATA"))
1074 scoreMean = computeRobustStatistics(output.scoreExposure.image,
1075 output.scoreExposure.mask,
1076 statsCtrl)
1077 self.assertFloatsAlmostEqual(scoreMean, 0, atol=5*meanError)
1078 nea = computePSFNoiseEquivalentArea(science.psf)
1079 # stddev of Score image should be close to expected value.
1080 scoreStd = computeRobustStatistics(output.scoreExposure.image, output.scoreExposure.mask,
1081 statsCtrl=statsCtrl, statistic=afwMath.STDEV)
1082 self.assertFloatsAlmostEqual(scoreStd, np.sqrt(2)*noiseLevel/np.sqrt(nea), rtol=0.1)
1084 def test_incomplete_template_coverage(self):
1085 noiseLevel = 1.
1086 border = 20
1087 xSize = 400
1088 ySize = 400
1089 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6,
1090 xSize=xSize, ySize=ySize)
1091 template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
1092 templateBorderSize=border, doApplyCalibration=True,
1093 xSize=xSize, ySize=ySize)
1095 science_height = science.getBBox().getDimensions().getY()
1097 def _run_and_check_coverage(template_coverage,
1098 requiredTemplateFraction=0.1,
1099 minTemplateFractionForExpectedSuccess=0.2):
1100 template_cut = template.clone()
1101 template_height = int(science_height*template_coverage + border)
1102 template_cut.image.array[:, template_height:] = 0
1103 template_cut.mask.array[:, template_height:] = template_cut.mask.getPlaneBitMask('NO_DATA')
1104 task = self._setup_subtraction(
1105 requiredTemplateFraction=requiredTemplateFraction,
1106 minTemplateFractionForExpectedSuccess=minTemplateFractionForExpectedSuccess
1107 )
1108 if template_coverage < task.config.requiredTemplateFraction:
1109 doRaise = True
1110 elif template_coverage < task.config.minTemplateFractionForExpectedSuccess:
1111 doRaise = True
1112 else:
1113 doRaise = False
1114 if doRaise:
1115 with self.assertRaises(NoWorkFound):
1116 task.run(template_cut, science.clone(), sources.copy(deep=True))
1117 else:
1118 task.run(template_cut, science.clone(), sources.copy(deep=True))
1119 _run_and_check_coverage(template_coverage=0.09)
1120 _run_and_check_coverage(template_coverage=0.15)
1121 _run_and_check_coverage(template_coverage=.7)
1123 def test_clear_template_mask(self):
1124 noiseLevel = 1.
1125 xSize = 400
1126 ySize = 400
1127 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6,
1128 xSize=xSize, ySize=ySize)
1129 template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
1130 templateBorderSize=20, doApplyCalibration=True,
1131 xSize=xSize, ySize=ySize)
1132 diffimEmptyMaskPlanes = ["DETECTED", "DETECTED_NEGATIVE"]
1133 task = self._setup_subtraction()
1134 # Ensure that each each mask plane is set for some pixels
1135 mask = template.mask
1136 x0 = 50
1137 x1 = 75
1138 y0 = 150
1139 y1 = 175
1140 scienceMaskCheck = {}
1141 for maskPlane in mask.getMaskPlaneDict().keys():
1142 scienceMaskCheck[maskPlane] = np.sum(science.mask.array & mask.getPlaneBitMask(maskPlane) > 0)
1143 mask.array[x0: x1, y0: y1] |= mask.getPlaneBitMask(maskPlane)
1144 self.assertTrue(np.sum(mask.array & mask.getPlaneBitMask(maskPlane) > 0))
1146 output = task.run(template, science, sources)
1147 # Verify that the template mask has been modified in place
1148 for maskPlane in mask.getMaskPlaneDict().keys():
1149 if maskPlane in diffimEmptyMaskPlanes:
1150 self.assertTrue(np.sum(mask.array & mask.getPlaneBitMask(maskPlane) == 0))
1151 elif maskPlane in task.config.preserveTemplateMask:
1152 self.assertTrue(np.sum(mask.array & mask.getPlaneBitMask(maskPlane) > 0))
1153 else:
1154 self.assertTrue(np.sum(mask.array & mask.getPlaneBitMask(maskPlane) == 0))
1155 # Mask planes set in the science image should also be set in the difference
1156 # Except the "DETECTED" planes should have been cleared
1157 diffimMask = output.scoreExposure.mask
1158 for maskPlane, scienceSum in scienceMaskCheck.items():
1159 diffimSum = np.sum(diffimMask.array & mask.getPlaneBitMask(maskPlane) > 0)
1160 if maskPlane in diffimEmptyMaskPlanes:
1161 self.assertEqual(diffimSum, 0)
1162 else:
1163 self.assertTrue(diffimSum >= scienceSum)
1165 def test_agnostic_template_psf(self):
1166 """Test that the Score image is the same whether the template PSF is
1167 larger or smaller than the science image PSF.
1168 """
1169 noiseLevel = .3
1170 xSize = 400
1171 ySize = 400
1172 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel,
1173 noiseSeed=6, templateBorderSize=0,
1174 xSize=xSize, ySize=ySize)
1175 template1, _ = makeTestImage(psfSize=self.badPsfSize, noiseLevel=noiseLevel,
1176 noiseSeed=7, doApplyCalibration=True,
1177 xSize=xSize, ySize=ySize)
1178 template2, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel,
1179 noiseSeed=8, doApplyCalibration=True,
1180 xSize=xSize, ySize=ySize)
1181 task = self._setup_subtraction()
1183 science_better = task.run(template1, science.clone(), sources)
1184 template_better = task.run(template2, science, sources)
1185 bbox = science_better.scoreExposure.getBBox().clippedTo(template_better.scoreExposure.getBBox())
1187 delta = template_better.scoreExposure[bbox].clone()
1188 delta.image -= science_better.scoreExposure[bbox].image
1189 delta.variance -= science_better.scoreExposure[bbox].variance
1190 delta.mask.array &= science_better.scoreExposure[bbox].mask.array
1192 statsCtrl = makeStats(badMaskPlanes=("EDGE", "BAD", "NO_DATA"))
1193 # Mean of delta should be very close to zero.
1194 nGoodPix = np.sum(np.isfinite(delta.image.array))
1195 meanError = 2*noiseLevel/np.sqrt(nGoodPix)
1196 deltaMean = computeRobustStatistics(delta.image, delta.mask, statsCtrl)
1197 deltaStd = computeRobustStatistics(delta.image, delta.mask, statsCtrl,
1198 statistic=afwMath.STDEV)
1199 self.assertFloatsAlmostEqual(deltaMean, 0, atol=5*meanError)
1200 nea = computePSFNoiseEquivalentArea(science.psf)
1201 # stddev of Score image should be close to expected value
1202 self.assertFloatsAlmostEqual(deltaStd, np.sqrt(2)*noiseLevel/np.sqrt(nea), rtol=.1)
1204 def test_few_sources(self):
1205 """Test with only 1 source, to check that we get a useful error.
1206 """
1207 xSize = 256
1208 ySize = 256
1209 science, sources = makeTestImage(psfSize=self.midPsfSize, nSrc=10, xSize=xSize, ySize=ySize)
1210 template, _ = makeTestImage(psfSize=self.goodPsfSize, nSrc=10, xSize=xSize, ySize=ySize,
1211 doApplyCalibration=True)
1212 task = self._setup_subtraction()
1213 sources = sources[0:1]
1214 with self.assertRaises(InsufficientKernelSourcesError):
1215 task.run(template, science, sources)
1217 def test_background_subtraction(self):
1218 """Check that we can recover the background,
1219 and that it is subtracted correctly in the Score image.
1220 """
1221 noiseLevel = 1.
1222 xSize = 512
1223 ySize = 512
1224 x0 = 123
1225 y0 = 456
1226 template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
1227 templateBorderSize=20,
1228 xSize=xSize, ySize=ySize, x0=x0, y0=y0,
1229 doApplyCalibration=True)
1230 params = [2.2, 2.1, 2.0, 1.2, 1.1, 1.0]
1232 bbox2D = lsst.geom.Box2D(lsst.geom.Point2D(x0, y0), lsst.geom.Extent2D(xSize, ySize))
1233 background_model = afwMath.Chebyshev1Function2D(params, bbox2D)
1234 science, sources = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=6,
1235 background=background_model,
1236 xSize=xSize, ySize=ySize, x0=x0, y0=y0)
1237 # Don't use ``self._setup_subtraction()`` here.
1238 # Modifying the config of a subtask is messy.
1239 config = subtractImages.AlardLuptonPreconvolveSubtractTask.ConfigClass()
1241 config.sourceSelector.signalToNoise.fluxField = "truth_instFlux"
1242 config.sourceSelector.signalToNoise.errField = "truth_instFluxErr"
1243 config.doSubtractBackground = True
1245 config.makeKernel.kernel.name = "AL"
1246 config.makeKernel.kernel.active.fitForBackground = True
1247 config.makeKernel.kernel.active.spatialKernelOrder = 1
1248 config.makeKernel.kernel.active.spatialBgOrder = 2
1249 statsCtrl = makeStats(badMaskPlanes=("EDGE", "BAD", "NO_DATA"))
1251 task = subtractImages.AlardLuptonPreconvolveSubtractTask(config=config)
1252 output = task.run(template.clone(), science.clone(), sources)
1254 # We should be fitting the same number of parameters as were in the input model
1255 self.assertEqual(output.backgroundModel.getNParameters(), background_model.getNParameters())
1257 # The parameters of the background fit should be close to the input model
1258 self.assertFloatsAlmostEqual(np.array(output.backgroundModel.getParameters()),
1259 np.array(params), rtol=0.2)
1261 # stddev of Score image should be close to expected value.
1262 # This will fail if we have mis-subtracted the background.
1263 stdVal = computeRobustStatistics(output.scoreExposure.image, output.scoreExposure.mask,
1264 statsCtrl, statistic=afwMath.STDEV)
1265 # get the img psf Noise Equivalent Area value
1266 nea = computePSFNoiseEquivalentArea(science.psf)
1267 self.assertFloatsAlmostEqual(stdVal, np.sqrt(2)*noiseLevel/np.sqrt(nea), rtol=0.12)
1269 def test_scale_variance(self):
1270 """Check variance scaling of the Score image.
1271 """
1272 scienceNoiseLevel = 4.
1273 templateNoiseLevel = 2.
1274 scaleFactor = 1.345
1275 xSize = 400
1276 ySize = 400
1277 # Make sure to include pixels with the DETECTED mask bit set.
1278 statsCtrl = makeStats(badMaskPlanes=("EDGE", "BAD", "NO_DATA"))
1280 def _run_and_check_images(science, template, sources, statsCtrl,
1281 doDecorrelation, doScaleVariance, scaleFactor=1.):
1282 """Check that the variance plane matches the expected value for
1283 different configurations of ``doDecorrelation`` and ``doScaleVariance``.
1284 """
1286 task = self._setup_subtraction(doDecorrelation=doDecorrelation,
1287 doScaleVariance=doScaleVariance,
1288 )
1289 output = task.run(template.clone(), science.clone(), sources)
1290 if doScaleVariance:
1291 self.assertFloatsAlmostEqual(task.metadata["scaleScienceVarianceFactor"],
1292 scaleFactor, atol=0.05)
1294 scienceNoise = computeRobustStatistics(science.variance, science.mask, statsCtrl)
1295 # get the img psf Noise Equivalent Area value
1296 nea = computePSFNoiseEquivalentArea(science.psf)
1297 scienceNoise /= nea
1298 if doDecorrelation:
1299 templateNoise = computeRobustStatistics(template.variance, template.mask, statsCtrl)
1300 templateNoise /= nea
1301 else:
1302 # Don't divide by NEA in this case, since the template is convolved
1303 # and in the same units as the Score exposure.
1304 templateNoise = computeRobustStatistics(output.matchedTemplate.variance,
1305 output.matchedTemplate.mask,
1306 statsCtrl)
1307 if doScaleVariance:
1308 # Only the science variance is scaled here. The template
1309 # variance is scaled independently in ``GetTemplateTask``.
1310 scienceNoise *= scaleFactor
1311 varMean = computeRobustStatistics(output.scoreExposure.variance,
1312 output.scoreExposure.mask,
1313 statsCtrl)
1314 self.assertFloatsAlmostEqual(varMean, scienceNoise + templateNoise, rtol=0.1)
1316 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=scienceNoiseLevel, noiseSeed=6,
1317 xSize=xSize, ySize=ySize)
1318 template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=templateNoiseLevel, noiseSeed=7,
1319 templateBorderSize=20, doApplyCalibration=True,
1320 xSize=xSize, ySize=ySize)
1321 # Verify that the variance plane of the Score image is correct
1322 # when the template and science variance planes are correct
1323 _run_and_check_images(science, template, sources, statsCtrl,
1324 doDecorrelation=True, doScaleVariance=True)
1325 _run_and_check_images(science, template, sources, statsCtrl,
1326 doDecorrelation=True, doScaleVariance=False)
1327 _run_and_check_images(science, template, sources, statsCtrl,
1328 doDecorrelation=False, doScaleVariance=True)
1329 _run_and_check_images(science, template, sources, statsCtrl,
1330 doDecorrelation=False, doScaleVariance=False)
1332 # Verify that the variance plane of the Score image is correct
1333 # when the input science variance plane is incorrect
1334 science.variance.array /= scaleFactor
1335 _run_and_check_images(science, template, sources, statsCtrl,
1336 doDecorrelation=True, doScaleVariance=True, scaleFactor=scaleFactor)
1337 _run_and_check_images(science, template, sources, statsCtrl,
1338 doDecorrelation=True, doScaleVariance=False, scaleFactor=scaleFactor)
1339 _run_and_check_images(science, template, sources, statsCtrl,
1340 doDecorrelation=False, doScaleVariance=True, scaleFactor=scaleFactor)
1341 _run_and_check_images(science, template, sources, statsCtrl,
1342 doDecorrelation=False, doScaleVariance=False, scaleFactor=scaleFactor)
1344 def test_exposure_properties(self):
1345 """Check that all necessary exposure metadata is included
1346 with the Score image.
1347 """
1348 noiseLevel = 1.
1349 xSize = 400
1350 ySize = 400
1351 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6,
1352 xSize=xSize, ySize=ySize)
1353 psf = science.psf
1354 psfAvgPos = psf.getAveragePosition()
1355 psfSize = getPsfFwhm(science.psf)
1356 psfImg = psf.computeKernelImage(psfAvgPos)
1357 template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
1358 templateBorderSize=20, doApplyCalibration=True,
1359 xSize=xSize, ySize=ySize)
1361 def _run_and_check_images(doDecorrelation):
1362 """Check that the metadata is correct with or without decorrelation.
1363 """
1364 task = self._setup_subtraction(doDecorrelation=doDecorrelation)
1365 output = task.run(template.clone(), science.clone(), sources)
1366 psfOut = output.scoreExposure.psf
1367 psfAvgPos = psfOut.getAveragePosition()
1368 if doDecorrelation:
1369 # Decorrelation requires recalculating the PSF,
1370 # so it will not be the same as the input
1371 psfOutSize = getPsfFwhm(science.psf)
1372 self.assertFloatsAlmostEqual(psfSize, psfOutSize)
1373 else:
1374 psfOutImg = psfOut.computeKernelImage(psfAvgPos)
1375 self.assertImagesAlmostEqual(psfImg, psfOutImg)
1377 # check PSF, WCS, bbox, filterLabel, photoCalib
1378 self.assertWcsAlmostEqualOverBBox(science.wcs, output.scoreExposure.wcs, science.getBBox())
1379 self.assertEqual(science.filter, output.scoreExposure.filter)
1380 self.assertEqual(science.photoCalib, output.scoreExposure.photoCalib)
1381 _run_and_check_images(doDecorrelation=True)
1382 _run_and_check_images(doDecorrelation=False)
1385class SimplifiedSubtractTest(AlardLuptonSubtractTestBase, lsst.utils.tests.TestCase):
1386 subtractTask = subtractImages.SimplifiedSubtractTask
1388 def test_runSimplifiedTaskWithExistingKernel(self):
1389 """Test that the simplified task produces the same output as
1390 `AlardLuptonSubtractTask` if it uses the AL kernel.
1391 """
1392 noiseLevel = 1.
1393 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6)
1394 template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
1395 templateBorderSize=20, doApplyCalibration=True)
1396 alTask = AlardLuptonSubtractTest._setup_subtraction(AlardLuptonSubtractTest())
1397 task = self._setup_subtraction(useExistingKernel=True)
1399 alResults = alTask.run(template.clone(), science.clone(), sources)
1400 results = task.run(template.clone(), science.clone(),
1401 inputPsfMatchingKernel=alResults.psfMatchingKernel)
1403 self.assertMaskedImagesEqual(alResults.difference, results.difference)
1405 def test_runSimplifiedTaskWithSourceDetection(self):
1406 """Test that the simplified task with source detection produces
1407 reasonable output.
1408 """
1409 noiseLevel = 1.
1410 science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6)
1411 template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
1412 templateBorderSize=20, doApplyCalibration=True)
1413 task = self._setup_subtraction(useExistingKernel=False,
1414 fluxField="base_PsfFlux_instFlux",
1415 errField="base_PsfFlux_instFluxErr",
1416 )
1418 output = task.run(template, science)
1420 # There shoud be no NaN values in the difference image
1421 self.assertTrue(np.all(np.isfinite(output.difference.image.array)))
1422 # Mean of difference image should be close to zero.
1423 meanError = noiseLevel/np.sqrt(output.difference.image.array.size)
1424 # Make sure to include pixels with the DETECTED mask bit set.
1425 statsCtrl = makeStats(badMaskPlanes=("EDGE", "BAD", "NO_DATA", "DETECTED", "DETECTED_NEGATIVE"))
1426 differenceMean = computeRobustStatistics(output.difference.image, output.difference.mask, statsCtrl)
1427 self.assertFloatsAlmostEqual(differenceMean, 0, atol=5*meanError)
1428 # stddev of difference image should be close to expected value.
1429 differenceStd = computeRobustStatistics(output.difference.image, output.difference.mask,
1430 makeStats(), statistic=afwMath.STDEV)
1431 self.assertFloatsAlmostEqual(differenceStd, np.sqrt(2)*noiseLevel, rtol=0.1)
1434def setup_module(module):
1435 lsst.utils.tests.init()
1438class MemoryTestCase(lsst.utils.tests.MemoryTestCase):
1439 pass
1442if __name__ == "__main__": 1442 ↛ 1443line 1442 didn't jump to line 1443 because the condition on line 1442 was never true
1443 lsst.utils.tests.init()
1444 unittest.main()