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