Coverage for tests/test_detectAndMeasure.py: 99%
1230 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-13 09:35 +0000
1# This file is part of ip_diffim.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
22import numpy as np
23import os
24import unittest
25import uuid
26from unittest import mock
27import requests
29import astropy.units
31import lsst.afw.detection as afwDetection
32import lsst.afw.geom as afwGeom
33import lsst.afw.image as afwImage
34import lsst.afw.table as afwTable
35import lsst.afw.math as afwMath
36import lsst.geom
37from lsst.images import (Box, DifferenceImage, DifferenceImageTemplateInfo, Polygon,
38 get_legacy_difference_image_mask_planes, get_legacy_optional_mask_planes)
39from lsst.images.convolution_kernels import ImageBasisConvolutionKernel
40from lsst.images.fields import ChebyshevField
41import lsst.images.psfs
42import lsst.images.serialization
43from lsst.images.tests import compare_masked_image_to_legacy
44from lsst.ip.diffim import detectAndMeasure, subtractImages
45from lsst.afw.table import IdFactory
46from lsst.afw.cameraGeom.testUtils import CameraWrapper, DetectorWrapper
47import lsst.meas.algorithms as measAlg
48import lsst.pipe.base
49import lsst.pipe.base.testUtils
50from lsst.pipe.base import InvalidQuantumError, UpstreamFailureNoWorkFound, AlgorithmError
51import lsst.utils.tests
52import lsst.meas.base.tests
53import lsst.daf.base as dafBase
54from lsst.afw.coord import Observatory, Weather
55import lsst.geom as geom
56import lsst.pex.config as pexConfig
58from utils import makeTestExposureRecord, makeTestImage, makeTestVisitInfo, checkMask
61class DetectAndMeasureTestBase:
63 def _check_diaSource(self, refSources, diaSource, refIds=None,
64 matchDistance=1., scale=1., usePsfFlux=True,
65 rtol=0.025, atol=None):
66 """Match a diaSource with a source in a reference catalog
67 and compare properties.
69 Parameters
70 ----------
71 refSources : `lsst.afw.table.SourceCatalog`
72 The reference catalog.
73 diaSource : `lsst.afw.table.SourceRecord`
74 The new diaSource to match to the reference catalog.
75 refIds : `list` of `int`, optional
76 Source IDs of previously associated diaSources.
77 matchDistance : `float`, optional
78 Maximum distance allowed between the detected and reference source
79 locations, in pixels.
80 scale : `float`, optional
81 Optional factor to scale the flux by before performing the test.
82 usePsfFlux : `bool`, optional
83 If set, test the PsfInstFlux field, otherwise use ApInstFlux.
84 rtol : `float`, optional
85 Relative tolerance of the flux value test.
86 atol : `float`, optional
87 Absolute tolerance of the flux value test.
88 """
89 if diaSource["sky_source"]:
90 return
91 distance = np.sqrt((diaSource.getX() - refSources.getX())**2
92 + (diaSource.getY() - refSources.getY())**2)
93 self.assertLess(min(distance), matchDistance)
94 src = refSources[np.argmin(distance)]
95 if refIds is not None:
96 # Check that the same source was not previously associated
97 self.assertNotIn(src.getId(), refIds)
98 refIds.append(src.getId())
99 if atol is None:
100 atol = rtol*src.getPsfInstFlux() if usePsfFlux else rtol*src.getApInstFlux()
101 if usePsfFlux:
102 self.assertFloatsAlmostEqual(src.getPsfInstFlux()*scale, diaSource.getPsfInstFlux(),
103 rtol=rtol, atol=atol)
104 else:
105 self.assertFloatsAlmostEqual(src.getApInstFlux()*scale, diaSource.getApInstFlux(),
106 rtol=rtol, atol=atol)
108 def _check_values(self, values, minValue=None, maxValue=None):
109 """Verify that an array has finite values, and optionally that they are
110 within specified minimum and maximum bounds.
112 Parameters
113 ----------
114 values : `numpy.ndarray`
115 Array of values to check.
116 minValue : `float`, optional
117 Minimum allowable value.
118 maxValue : `float`, optional
119 Maximum allowable value.
120 """
121 self.assertTrue(np.all(np.isfinite(values)))
122 if minValue is not None:
123 self.assertTrue(np.all(values >= minValue))
124 if maxValue is not None:
125 self.assertTrue(np.all(values <= maxValue))
127 def _setup_detection(self, doSkySources=True, nSkySources=5,
128 doSubtractBackground=False, run_sattle=False,
129 badSourceFlags=None, **kwargs):
130 """Setup and configure the detection and measurement PipelineTask.
132 Parameters
133 ----------
134 doSkySources : `bool`, optional
135 Generate sky sources.
136 nSkySources : `int`, optional
137 The number of sky sources to add in isolated background regions.
138 **kwargs
139 Any additional config parameters to set.
141 Returns
142 -------
143 `lsst.pipe.base.PipelineTask`
144 The configured Task to use for detection and measurement.
145 """
146 config = self.detectionTask.ConfigClass()
147 config.doSkySources = doSkySources
148 if doSkySources:
149 config.skySources.nSources = nSkySources
150 config.update(**kwargs)
152 config.run_sattle = run_sattle
154 # Make a realistic id generator so that output catalog ids are useful.
155 dataId = lsst.daf.butler.DataCoordinate.standardize(
156 instrument="I",
157 visit=42,
158 detector=12,
159 universe=lsst.daf.butler.DimensionUniverse(),
160 )
161 if badSourceFlags is None:
162 badSourceFlags = ["base_PixelFlags_flag_offimage", ]
163 config.idGenerator.packer.name = "observation"
164 config.idGenerator.packer["observation"].n_observations = 10000
165 config.idGenerator.packer["observation"].n_detectors = 99
166 config.idGenerator.n_releases = 8
167 config.idGenerator.release_id = 2
168 config.doSubtractBackground = doSubtractBackground
169 config.badSourceFlags = badSourceFlags
170 self.idGenerator = config.idGenerator.apply(dataId)
172 return self.detectionTask(config=config)
175class DetectAndMeasureTest(DetectAndMeasureTestBase, lsst.utils.tests.TestCase):
176 detectionTask = detectAndMeasure.DetectAndMeasureTask
178 def test_detection_xy0(self):
179 """Basic functionality test with non-zero x0 and y0.
180 """
181 # Set up the simulated images
182 noiseLevel = 1.
183 staticSeed = 1
184 fluxLevel = 500
185 psfSize = 2.4
186 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel, "x0": 12345, "y0": 67890}
187 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
188 science.getInfo().setVisitInfo(makeVisitInfo())
189 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
190 difference = science.clone()
192 # Configure the detection Task
193 # Set the subtraction residual metric threshold high, since the template
194 # is not actually subtracted from the science image.
195 detectionTask = self._setup_detection(doDeblend=True, badSubtractionRatioThreshold=1.,
196 doSkySources=False)
198 # Run detection and check the results
199 output = detectionTask.run(science, matchedTemplate, difference, sources,
200 idFactory=self.idGenerator.make_table_id_factory())
201 subtractedMeasuredExposure = output.subtractedMeasuredExposure
203 # Catalog ids should be very large from this id generator.
204 self.assertTrue(all(output.diaSources['id'] > 1000000000))
205 self.assertImagesEqual(subtractedMeasuredExposure.image, difference.image)
207 # all of the sources should have been detected
208 self.assertEqual(len(output.diaSources), len(sources))
209 # no sources should be flagged as negative
210 self.assertEqual(len(~output.diaSources["is_negative"]), len(output.diaSources))
211 refIds = []
212 for source in sources:
213 self._check_diaSource(output.diaSources, source, refIds=refIds)
215 def test_raise_bad_psf(self):
216 """Detection should raise if the PSF width is NaN
217 """
218 # Set up the simulated images
219 noiseLevel = 1.
220 staticSeed = 1
221 fluxLevel = 500
222 psfSize = 2.4
223 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel, "x0": 12345, "y0": 67890}
224 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
225 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
226 difference = science.clone()
228 psfShape = difference.getPsf().computeBBox(lsst.geom.Point2D(0, 0)).getDimensions()
229 psfcI = afwImage.ImageD(lsst.geom.Extent2I(psfShape[1], psfShape[0]))
230 psfNew = np.zeros(psfShape)
231 psfNew[0, :] = 1
232 psfcI.array = psfNew
233 psfcK = afwMath.FixedKernel(psfcI)
234 psf = measAlg.KernelPsf(psfcK)
235 difference.setPsf(psf)
237 # Configure the detection Task
238 detectionTask = self._setup_detection()
240 # Run detection and check the results
241 with self.assertRaises(UpstreamFailureNoWorkFound):
242 detectionTask.run(science, matchedTemplate, difference, sources)
244 def test_measurements_finite(self):
245 """Measured fluxes and centroids should always be finite.
246 """
247 columnNames = ["coord_ra", "coord_dec", "ip_diffim_forced_PsfFlux_instFlux",
248 "ip_diffim_forced_template_PsfFlux_instFlux"]
250 # Set up the simulated images
251 noiseLevel = 1.
252 staticSeed = 1
253 transientSeed = 6
254 xSize = 256
255 ySize = 256
256 psfSize = 2.4
257 kwargs = {"psfSize": psfSize, "x0": 0, "y0": 0,
258 "xSize": xSize, "ySize": ySize}
259 science, sources = makeTestImage(seed=staticSeed, noiseLevel=noiseLevel, noiseSeed=6,
260 nSrc=1, **kwargs)
261 science.getInfo().setVisitInfo(makeVisitInfo())
262 matchedTemplate, _ = makeTestImage(seed=staticSeed, noiseLevel=noiseLevel/4, noiseSeed=7,
263 nSrc=1, **kwargs)
264 rng = np.random.RandomState(3)
265 xLoc = np.arange(-5, xSize+5, 10)
266 rng.shuffle(xLoc)
267 yLoc = np.arange(-5, ySize+5, 10)
268 rng.shuffle(yLoc)
269 transients, transientSources = makeTestImage(seed=transientSeed,
270 nSrc=len(xLoc), fluxLevel=1000.,
271 noiseLevel=noiseLevel, noiseSeed=8,
272 xLoc=xLoc, yLoc=yLoc,
273 **kwargs)
274 difference = science.clone()
275 difference.maskedImage -= matchedTemplate.maskedImage
276 difference.maskedImage += transients.maskedImage
278 # Configure the detection Task
279 detectionTask = self._setup_detection(doForcedMeasurement=True)
281 # Run detection and check the results
282 output = detectionTask.run(science, matchedTemplate, difference, sources)
284 for column in columnNames:
285 self._check_values(output.diaSources[column])
286 self._check_values(output.diaSources.getX(), minValue=0, maxValue=xSize)
287 self._check_values(output.diaSources.getY(), minValue=0, maxValue=ySize)
288 self._check_values(output.diaSources.getPsfInstFlux())
290 def test_raise_config_schema_mismatch(self):
291 """Check that sources with specified flags are removed from the catalog.
292 """
293 # Configure the detection Task, and and set a config that is not in the schema
294 with self.assertRaises(InvalidQuantumError):
295 self._setup_detection(badSourceFlags=["Bogus_flag_42"])
297 def test_remove_unphysical(self):
298 """Check that sources with specified flags are removed from the catalog.
299 """
300 # Set up the simulated images
301 noiseLevel = 1.
302 staticSeed = 1
303 xSize = 256
304 ySize = 256
305 psfSize = 2.4
306 kwargs = {"psfSize": psfSize, "xSize": xSize, "ySize": ySize}
307 science, sources = makeTestImage(seed=staticSeed, noiseLevel=noiseLevel, noiseSeed=6,
308 nSrc=1, **kwargs)
309 science.getInfo().setVisitInfo(makeVisitInfo())
310 matchedTemplate, _ = makeTestImage(seed=staticSeed, noiseLevel=noiseLevel/4, noiseSeed=7,
311 nSrc=1, **kwargs)
312 transients, transientSources = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=8, nSrc=1, **kwargs)
313 science.maskedImage += transients.maskedImage
314 difference = science.clone()
315 bbox = difference.getBBox()
316 difference.maskedImage -= matchedTemplate.maskedImage
318 # Configure the detection Task, and remove unphysical sources
319 detectionTask = self._setup_detection(doForcedMeasurement=False, doSkySources=True, nSkySources=20,
320 badSourceFlags=["base_PixelFlags_flag_offimage", ])
322 # Run detection and check the results
323 diaSources = detectionTask.run(science, matchedTemplate, difference, sources).diaSources
324 badDiaSrcDoRemove = ~bbox.contains(diaSources.getX(), diaSources.getY())
325 nBadDoRemove = np.count_nonzero(badDiaSrcDoRemove)
326 # Verify that all sources are physical
327 self.assertEqual(nBadDoRemove, 0)
328 # Set a few centroids outside the image bounding box
329 nSetBad = 5
330 for src in diaSources[0: nSetBad]:
331 src["slot_Centroid_x"] += xSize
332 src["slot_Centroid_y"] += ySize
333 src["base_PixelFlags_flag_offimage"] = True
334 # Verify that these sources are outside the image
335 badDiaSrc = ~bbox.contains(diaSources.getX(), diaSources.getY())
336 nBad = np.count_nonzero(badDiaSrc)
337 self.assertEqual(nBad, nSetBad)
338 diaSourcesNoBad = detectionTask._removeBadSources(diaSources)
339 badDiaSrcNoBad = ~bbox.contains(diaSourcesNoBad.getX(), diaSourcesNoBad.getY())
341 # Verify that no sources outside the image bounding box remain
342 self.assertEqual(np.count_nonzero(badDiaSrcNoBad), 0)
343 self.assertEqual(len(diaSourcesNoBad), len(diaSources) - nSetBad)
345 def test_detect_transients(self):
346 """Run detection on a difference image containing transients.
347 """
348 # Set up the simulated images
349 noiseLevel = 1.
350 staticSeed = 1
351 transientSeed = 6
352 nTransients = 10
353 transientFlux = 1000
354 fluxLevel = 500
355 psfSize = 2.4
356 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel}
357 scienceBase, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
358 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
360 # Configure the detection Task
361 detectionTask = self._setup_detection(doMerge=False, doSkySources=True)
362 kwargs["seed"] = transientSeed
363 kwargs["nSrc"] = nTransients
364 kwargs["fluxLevel"] = transientFlux
366 # Run detection and check the results
367 def _run_and_check_detections(positive=True):
368 transients, transientSources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, **kwargs)
369 science = scienceBase.clone()
370 if positive:
371 science.maskedImage += transients.maskedImage
372 else:
373 science.maskedImage -= transients.maskedImage
374 difference = science.clone()
375 difference.maskedImage -= matchedTemplate.maskedImage
376 # NOTE: NoiseReplacer (run by forcedMeasurement) can modify the
377 # science image if we've e.g. removed parents post-deblending.
378 # Pass a clone of the science image, so that it doesn't disrupt
379 # later tests.
380 output = detectionTask.run(science, matchedTemplate, difference, sources)
381 refIds = []
382 scale = 1. if positive else -1.
383 for diaSource in output.diaSources:
384 self._check_diaSource(transientSources, diaSource, refIds=refIds, scale=scale)
385 _run_and_check_detections(positive=True)
386 _run_and_check_detections(positive=False)
388 def test_detect_transients_with_background(self):
389 """Run detection on a difference image containing transients and a background.
390 """
391 # Set up the simulated images
392 noiseLevel = 1.
393 staticSeed = 1
394 transientSeed = 6
395 nTransients = 10
396 transientFlux = 1000
397 fluxLevel = 500
398 xSize = 512
399 ySize = 512
400 x0 = 123
401 y0 = 456
402 psfSize = 2.4
403 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel,
404 "xSize": xSize, "ySize": ySize, "x0": x0, "y0": y0}
405 params = [2.2, 2.1, 2.0, 1.2, 1.1, 1.0]
407 bbox2D = lsst.geom.Box2D(lsst.geom.Point2D(x0, y0), lsst.geom.Extent2D(xSize, ySize))
408 background_model = afwMath.Chebyshev1Function2D(params, bbox2D)
409 scienceBase, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6,
410 background=background_model, **kwargs)
411 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
413 # Configure the detection Task
414 detectionTask = self._setup_detection(doMerge=False, doSubtractBackground=True)
415 kwargs["seed"] = transientSeed
416 kwargs["nSrc"] = nTransients
417 kwargs["fluxLevel"] = transientFlux
419 # Run detection and check the results
420 def _run_and_check_detections(positive=True):
421 transients, transientSources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, **kwargs)
422 science = scienceBase.clone()
423 if positive:
424 science.maskedImage += transients.maskedImage
425 else:
426 science.maskedImage -= transients.maskedImage
427 difference = science.clone()
428 difference.maskedImage -= matchedTemplate.maskedImage
429 # NOTE: NoiseReplacer (run by forcedMeasurement) can modify the
430 # science image if we've e.g. removed parents post-deblending.
431 # Pass a clone of the science image, so that it doesn't disrupt
432 # later tests.
433 output = detectionTask.run(science, matchedTemplate, difference, sources)
434 refIds = []
435 scale = 1. if positive else -1.
436 for transient in transientSources:
437 self._check_diaSource(output.diaSources, transient, refIds=refIds, scale=scale)
438 _run_and_check_detections(positive=True)
439 _run_and_check_detections(positive=False)
441 def test_mask_cosmic_rays(self):
442 """Run detection on a difference image containing a cosmic ray.
443 """
444 # Set up the simulated images
445 noiseLevel = 1.
446 staticSeed = 1
447 fluxLevel = 500
448 xSize = 400
449 ySize = 400
450 psfSize = 2.4
451 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel,
452 "xSize": xSize, "ySize": ySize}
453 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
454 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
455 crMask = science.mask.getPlaneBitMask("CR")
457 # Configure the detection Task
458 detectionTask = self._setup_detection(doMerge=False, doSkySources=True)
460 # Test that no CRs are detected
461 transients, transientSources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, nSrc=10, **kwargs)
462 science.maskedImage += transients.maskedImage
463 difference = science.clone()
464 difference.maskedImage -= matchedTemplate.maskedImage
465 output = detectionTask.run(science, matchedTemplate, difference, sources)
466 crMaskSet = (output.subtractedMeasuredExposure.mask.array & crMask) > 0
467 self.assertTrue(np.all(crMaskSet == 0))
469 crX0 = round(sources.getX()[0] - science.getX0())
470 crY0 = round(sources.getY()[0] - science.getY0())
471 crX1 = round(sources.getX()[1] - science.getX0())
472 crY1 = round(sources.getY()[1] - science.getY0())
473 # Add CR-like shape and check that CR is detected
474 # Pick two locations on top of sources, since that is what is likely to
475 # be missed in the first stage of CR rejection.
476 crMaskSetInput = np.zeros(difference.image.array.shape, bool)
477 crMaskSetInput[crX0:crX0+1, crY0:crY0+5] = True
478 crMaskSetInput[crX1:crX1+5, crY1:crY1+1] = True
479 difference.image.array[crMaskSetInput] += 1234
480 output = detectionTask.run(science, matchedTemplate, difference)
481 crMaskSet = (output.subtractedMeasuredExposure.mask.array & crMask) > 0
482 self.assertFalse(np.any(crMaskSet[~crMaskSetInput]))
483 self.assertTrue(np.all(crMaskSet[crMaskSetInput]))
485 def test_missing_mask_planes(self):
486 """Check that detection runs with missing mask planes.
487 """
488 # Set up the simulated images
489 noiseLevel = 1.
490 fluxLevel = 500
491 psfSize = 2.4
492 kwargs = {"psfSize": psfSize, "fluxLevel": fluxLevel, "addMaskPlanes": []}
493 # Use different seeds for the science and template so every source is a diaSource
494 science, sources = makeTestImage(seed=5, noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
495 science.getInfo().setVisitInfo(makeVisitInfo())
496 matchedTemplate, _ = makeTestImage(seed=6, noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
498 difference = science.clone()
499 difference.maskedImage -= matchedTemplate.maskedImage
500 detectionTask = self._setup_detection(raiseOnBadSubtractionRatio=False, raiseOnNoDiaSources=False)
502 # Verify that detection runs without errors
503 detectionTask.run(science, matchedTemplate, difference, sources)
505 def test_detect_dipoles(self):
506 """Run detection on a difference image containing dipoles.
507 """
508 # Set up the simulated images
509 noiseLevel = 1.
510 staticSeed = 1
511 fluxLevel = 1000
512 fluxRange = 1.5
513 nSources = 10
514 offset = 1
515 xSize = 300
516 ySize = 300
517 kernelSize = 31
518 psfSize = 2.4
519 # Avoid placing sources near the edge for this test, so that we can
520 # easily check that the correct number of sources are detected.
521 templateBorderSize = kernelSize//2
522 dipoleFlag = "ip_diffim_DipoleFit_classification"
523 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel, "fluxRange": fluxRange,
524 "nSrc": nSources, "templateBorderSize": templateBorderSize, "kernelSize": kernelSize,
525 "xSize": xSize, "ySize": ySize}
526 dipoleFlag = "ip_diffim_DipoleFit_classification"
527 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
528 science.getInfo().setVisitInfo(makeVisitInfo())
529 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
530 difference = science.clone()
531 matchedTemplate.image.array[...] = np.roll(matchedTemplate.image.array[...], offset, axis=0)
532 matchedTemplate.variance.array[...] = np.roll(matchedTemplate.variance.array[...], offset, axis=0)
533 matchedTemplate.mask.array[...] = np.roll(matchedTemplate.mask.array[...], offset, axis=0)
534 difference.maskedImage -= matchedTemplate.maskedImage[science.getBBox()]
536 # Need a higher residual metric threshold, since we are purposely
537 # creating poor subtractions.
538 detectionTask = self._setup_detection(doMerge=True, badSubtractionRatioThreshold=0.3)
539 output = detectionTask.run(science, matchedTemplate, difference, sources)
540 diaSources = output.diaSources[~output.diaSources["sky_source"]].copy(deep=True)
541 self.assertEqual(len(diaSources), len(sources))
542 # no sources should be flagged as negative
543 self.assertEqual(len(~diaSources["is_negative"]), len(diaSources))
544 refIds = []
545 for diaSource in diaSources:
546 if diaSource[dipoleFlag]: 546 ↛ 553line 546 didn't jump to line 553 because the condition on line 546 was always true
547 self._check_diaSource(sources, diaSource, refIds=refIds, scale=0,
548 rtol=0.05, atol=None, usePsfFlux=False)
549 self.assertFloatsAlmostEqual(diaSource["ip_diffim_DipoleFit_orientation"],
550 -np.pi / 2, atol=2.)
551 self.assertFloatsAlmostEqual(diaSource["ip_diffim_DipoleFit_separation"], offset, rtol=0.1)
552 else:
553 raise ValueError("DiaSource with ID %s is not a dipole!", diaSource.getId())
555 def test_sky_sources(self):
556 """Add sky sources and check that they are sufficiently far from other
557 sources and have negligible flux.
558 """
559 # Set up the simulated images
560 noiseLevel = 1.
561 staticSeed = 1
562 transientSeed = 6
563 transientFluxLevel = 1000.
564 transientFluxRange = 1.5
565 fluxLevel = 500
566 psfSize = 2.4
567 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel}
568 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
569 science.getInfo().setVisitInfo(makeVisitInfo())
570 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
571 transients, transientSources = makeTestImage(seed=transientSeed, psfSize=2.4,
572 nSrc=10, fluxLevel=transientFluxLevel,
573 fluxRange=transientFluxRange,
574 noiseLevel=noiseLevel, noiseSeed=8)
575 difference = science.clone()
576 difference.maskedImage -= matchedTemplate.maskedImage
577 difference.maskedImage += transients.maskedImage
578 kernelWidth = np.max(science.psf.getKernel().getDimensions())//2
580 # Configure the detection Task
581 detectionTask = self._setup_detection(doSkySources=True)
583 # Run detection and check the results
584 output = detectionTask.run(science, matchedTemplate, difference, sources,
585 idFactory=self.idGenerator.make_table_id_factory())
586 skySources = output.diaSources[output.diaSources["sky_source"]]
587 self.assertEqual(len(skySources), detectionTask.config.skySources.nSources)
588 for skySource in skySources:
589 # Disable the sky_source flag to enable checks.
590 skySource["sky_source"] = False
591 # The sky sources should not be close to any other source
592 with self.assertRaises(AssertionError):
593 self._check_diaSource(transientSources, skySource, matchDistance=kernelWidth)
594 with self.assertRaises(AssertionError):
595 self._check_diaSource(sources, skySource, matchDistance=kernelWidth)
596 # The sky sources should have low flux levels.
597 self._check_diaSource(transientSources, skySource, matchDistance=1000, scale=0.,
598 atol=np.sqrt(transientFluxRange*transientFluxLevel))
600 # Catalog ids should be very large from this id generator.
601 self.assertTrue(all(output.diaSources['id'] > 1000000000))
603 def test_exclude_mask_detections(self):
604 """Sources with certain bad mask planes set should not be detected.
605 """
606 # Set up the simulated images
607 noiseLevel = 1.
608 staticSeed = 1
609 transientSeed = 6
610 fluxLevel = 500
611 radius = 2
612 psfSize = 2.4
613 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel}
614 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
615 science.getInfo().setVisitInfo(makeVisitInfo())
616 detector = DetectorWrapper(numAmps=1).detector
617 science.setDetector(detector)
618 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
620 # Configure the detection Task
621 badSourceFlags = ["base_PixelFlags_flag_offimage",
622 "base_PixelFlags_flag_edgeCenterAll",
623 "base_PixelFlags_flag_badCenterAll",
624 "base_PixelFlags_flag_saturatedCenterAll", ]
625 detectionTask = self._setup_detection(nSkySources=0, badSourceFlags=badSourceFlags)
626 excludeMaskPlanes = ["EDGE", "BAD", "SAT"]
627 nBad = len(excludeMaskPlanes)
628 self.assertGreater(nBad, 0)
629 kwargs["seed"] = transientSeed
630 kwargs["nSrc"] = nBad
631 kwargs["fluxLevel"] = 1000
633 # Run detection and check the results
634 def _run_and_check_detections(setFlags=True):
635 transients, transientSources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, **kwargs)
636 difference = science.clone()
637 difference.maskedImage -= matchedTemplate.maskedImage
638 difference.maskedImage += transients.maskedImage
639 if setFlags:
640 for src, badMask in zip(transientSources, excludeMaskPlanes):
641 srcX = int(src.getX())
642 srcY = int(src.getY())
643 srcBbox = lsst.geom.Box2I(lsst.geom.Point2I(srcX - radius, srcY - radius),
644 lsst.geom.Extent2I(2*radius + 1, 2*radius + 1))
645 difference[srcBbox].mask.array |= lsst.afw.image.Mask.getPlaneBitMask(badMask)
646 with self.assertRaises(AlgorithmError):
647 output = detectionTask.run(science, matchedTemplate, difference, sources)
648 else:
649 output = detectionTask.run(science, matchedTemplate, difference, sources)
650 refIds = []
651 goodSrcFlags = checkMask(difference.mask, transientSources, excludeMaskPlanes)
652 self.assertEqual(np.sum(~goodSrcFlags), 0)
653 for diaSource, goodSrcFlag in zip(output.diaSources, goodSrcFlags):
654 if ~goodSrcFlag: 654 ↛ 655line 654 didn't jump to line 655 because the condition on line 654 was never true
655 with self.assertRaises(AssertionError):
656 self._check_diaSource(transientSources, diaSource, refIds=refIds)
657 else:
658 self._check_diaSource(transientSources, diaSource, refIds=refIds)
659 _run_and_check_detections(setFlags=False)
660 _run_and_check_detections(setFlags=True)
662 def test_fake_mask_plane_propagation(self):
663 """Test that we have the mask planes related to fakes in diffim images.
664 This is testing method called updateMasks
665 """
666 xSize = 256
667 ySize = 256
668 science, sources = makeTestImage(psfSize=2.4, xSize=xSize, ySize=ySize, doApplyCalibration=True)
669 science_fake_img, science_fake_sources = makeTestImage(
670 psfSize=2.4, xSize=xSize, ySize=ySize, seed=5, nSrc=3, noiseLevel=0.25, fluxRange=1
671 )
672 template, _ = makeTestImage(psfSize=2.4, xSize=xSize, ySize=ySize, doApplyCalibration=True)
673 tmplt_fake_img, tmplt_fake_sources = makeTestImage(
674 psfSize=2.4, xSize=xSize, ySize=ySize, seed=9, nSrc=3, noiseLevel=0.25, fluxRange=1
675 )
676 # created fakes and added them to the images
677 science.image += science_fake_img.image
678 template.image += tmplt_fake_img.image
680 # TODO: DM-40796 update to INJECTED names when source injection gets refactored
681 # adding mask planes to both science and template images
682 science.mask.addMaskPlane("FAKE")
683 science_fake_bitmask = science.mask.getPlaneBitMask("FAKE")
684 template.mask.addMaskPlane("FAKE")
685 template_fake_bitmask = template.mask.getPlaneBitMask("FAKE")
687 # makeTestImage sets the DETECTED plane on the sources; we can use
688 # that to set the FAKE plane on the science and template images.
689 detected = science_fake_img.mask.getPlaneBitMask("DETECTED")
690 fake_pixels = (science_fake_img.mask.array & detected).nonzero()
691 science.mask.array[fake_pixels] |= science_fake_bitmask
692 detected = tmplt_fake_img.mask.getPlaneBitMask("DETECTED")
693 fake_pixels = (tmplt_fake_img.mask.array & detected).nonzero()
694 template.mask.array[fake_pixels] |= science_fake_bitmask
696 science_fake_masked = (science.mask.array & science_fake_bitmask) > 0
697 template_fake_masked = (template.mask.array & template_fake_bitmask) > 0
699 subtractConfig = subtractImages.AlardLuptonSubtractTask.ConfigClass()
700 subtractConfig.sourceSelector.signalToNoise.fluxField = "truth_instFlux"
701 subtractConfig.sourceSelector.signalToNoise.errField = "truth_instFluxErr"
702 subtractTask = subtractImages.AlardLuptonSubtractTask(config=subtractConfig)
703 subtraction = subtractTask.run(template, science, sources)
705 # check subtraction mask plane is set where we set the previous masks
706 diff_mask = subtraction.difference.mask
708 # science mask should be now in INJECTED
709 inj_masked = (diff_mask.array & diff_mask.getPlaneBitMask("INJECTED")) > 0
711 # template mask should be now in INJECTED_TEMPLATE
712 injTmplt_masked = (diff_mask.array & diff_mask.getPlaneBitMask("INJECTED_TEMPLATE")) > 0
714 self.assertFloatsEqual(inj_masked.astype(int), science_fake_masked.astype(int))
715 # The template is convolved, so the INJECTED_TEMPLATE mask plane may
716 # include more pixels than the FAKE mask plane
717 injTmplt_masked &= template_fake_masked
718 self.assertFloatsEqual(injTmplt_masked.astype(int), template_fake_masked.astype(int))
720 # Now check that detection of fakes have the correct flag for injections
721 detectionTask = self._setup_detection()
723 output = detectionTask.run(subtraction.matchedScience,
724 subtraction.matchedTemplate,
725 subtraction.difference,
726 subtraction.kernelSources)
727 diaSources = output.diaSources[~output.diaSources["sky_source"]].copy(deep=True)
729 sci_refIds = []
730 tmpl_refIds = []
731 for diaSrc in diaSources:
732 if diaSrc['base_PsfFlux_instFlux'] > 0:
733 self._check_diaSource(science_fake_sources, diaSrc, scale=1, refIds=sci_refIds)
734 self.assertTrue(diaSrc['base_PixelFlags_flag_injected'])
735 self.assertTrue(diaSrc['base_PixelFlags_flag_injectedCenter'])
736 self.assertFalse(diaSrc['base_PixelFlags_flag_injected_template'])
737 self.assertFalse(diaSrc['base_PixelFlags_flag_injected_templateCenter'])
738 else:
739 self._check_diaSource(tmplt_fake_sources, diaSrc, scale=-1, refIds=tmpl_refIds)
740 self.assertTrue(diaSrc['base_PixelFlags_flag_injected_template'])
741 self.assertTrue(diaSrc['base_PixelFlags_flag_injected_templateCenter'])
742 self.assertFalse(diaSrc['base_PixelFlags_flag_injected'])
743 self.assertFalse(diaSrc['base_PixelFlags_flag_injectedCenter'])
745 def test_mask_streaks(self):
746 """Run detection on a difference image containing a streak.
747 """
748 # Set up the simulated images
749 noiseLevel = 1.
750 staticSeed = 1
751 fluxLevel = 500
752 xSize = 400
753 ySize = 400
754 psfSize = 2.4
755 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel,
756 "xSize": xSize, "ySize": ySize}
757 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
758 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
760 # Configure the detection Task
761 detectionTask = self._setup_detection(doMerge=False, doMaskStreaks=True,
762 raiseOnNoDiaSources=False, raiseOnBadSubtractionRatio=False)
764 # Test that no streaks are detected
765 difference = science.clone()
766 difference.maskedImage -= matchedTemplate.maskedImage
767 output = detectionTask.run(science, matchedTemplate, difference, sources)
768 outMask = output.subtractedMeasuredExposure.mask.array
769 streakMask = output.subtractedMeasuredExposure.mask.getPlaneBitMask("STREAK")
770 streakMaskSet = (outMask & streakMask) > 0
771 self.assertTrue(np.all(streakMaskSet == 0))
773 # Add streak-like shape and check that streak is detected
774 difference.image.array[20:23, 40:200] += 50
775 output = detectionTask.run(science, matchedTemplate, difference, sources)
776 outMask = output.subtractedMeasuredExposure.mask.array
777 streakMask = output.subtractedMeasuredExposure.mask.getPlaneBitMask("STREAK")
778 streakMaskSet = (outMask & streakMask) > 0
779 self.assertTrue(np.all(streakMaskSet[20:23, 40:200]))
781 # Add a more trapezoid shaped streak across an image that is
782 # fainter and check that it is also detected
783 bbox = science.getBBox()
784 difference = science.clone()
785 difference.maskedImage -= matchedTemplate.maskedImage
786 width = 100
787 amplitude = 5
788 x0 = -100 + bbox.getBeginX()
789 y0 = -100 + bbox.getBeginY()
790 x1 = xSize + x0 + 100
791 y1 = ySize/2
792 corner_coords = [lsst.geom.Point2D(x0, y0),
793 lsst.geom.Point2D(x0, y0 + width),
794 lsst.geom.Point2D(x1, y1),
795 lsst.geom.Point2D(x1 + width, y1)]
796 streak_trapezoid = afwGeom.Polygon(corner_coords)
797 streak_image = streak_trapezoid.createImage(bbox)
798 streak_image *= amplitude
799 difference.image.array += streak_image.array
800 output = detectionTask.run(science, matchedTemplate, difference, sources)
801 outMask = output.subtractedMeasuredExposure.mask.array
802 streakMask = output.subtractedMeasuredExposure.mask.getPlaneBitMask("STREAK")
803 streakMaskSet = (outMask & streakMask) > 0
804 # Check that pixel values in streak_image equal the streak amplitude
805 streak_check = np.where(streak_image.array == amplitude)
806 self.assertTrue(np.all(streakMaskSet[streak_check]))
807 # Check that the entire image was not masked STREAK
808 self.assertFalse(np.all(streakMaskSet))
810 def _setup_sattle_tests(self):
811 noiseLevel = 1.
812 staticSeed = 1
813 fluxLevel = 500
814 psfSize = 2.4
815 shared_kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel,
816 "x0": 12345, "y0": 67890}
817 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6,
818 **shared_kwargs)
819 science.getInfo().setVisitInfo(makeVisitInfo(id=2))
820 detector = DetectorWrapper(numAmps=1).detector
821 science.setDetector(detector)
822 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel / 4,
823 noiseSeed=7, **shared_kwargs)
824 difference = science.clone()
826 detectionTask = self._setup_detection(doDeblend=True,
827 badSubtractionRatioThreshold=1.,
828 doSkySources=False, run_sattle=True)
830 return science, matchedTemplate, difference, sources, detectionTask
832 @mock.patch.dict(os.environ, {"SATTLE_URI_BASE": "fake_host:1234"})
833 def test_sattle_not_available(self):
834 science, matchedTemplate, difference, sources, detectionTask = self._setup_sattle_tests()
836 response = MockResponse({"allow_list": []}, 500, "sattle internal error")
837 with mock.patch('requests.put', return_value=response):
838 with self.assertRaises(requests.exceptions.HTTPError):
839 detectionTask.run(science, matchedTemplate, difference, sources,
840 idFactory=IdFactory.makeSimple())
842 @mock.patch.dict(os.environ, {"SATTLE_URI_BASE": "fake_host:1234"})
843 def test_visit_id_not_in_sattle(self):
844 science, matchedTemplate, difference, sources, detectionTask = self._setup_sattle_tests()
846 response = MockResponse({"allow_list": []}, 404, "missing visit cache")
847 # visit id not in sattle raises
848 with self.assertRaises(requests.exceptions.HTTPError):
849 with mock.patch('lsst.ip.diffim.detectAndMeasure.requests.put',
850 return_value=response):
851 with mock.patch('lsst.ip.diffim.utils.populate_sattle_visit_cache'):
852 detectionTask.run(science, matchedTemplate, difference, sources,
853 idFactory=IdFactory.makeSimple())
855 @mock.patch.dict(os.environ, {"SATTLE_URI_BASE": "fake_host:1234"})
856 def test_filter_satellites_some_allowed(self):
857 science, matchedTemplate, difference, sources, detectionTask = self._setup_sattle_tests()
859 allowed_ids = [1, 5]
860 response = MockResponse({"allow_list": allowed_ids}, 200, "some allowed")
861 with mock.patch('requests.put', return_value=response):
862 output = detectionTask.run(science, matchedTemplate, difference, sources,
863 idFactory=IdFactory.makeSimple())
865 self.assertEqual(len(output.diaSources), 2)
867 # Output should be sources 1 and 5 allowed out of 20
868 self.assertEqual(set(output.diaSources['id']), set(allowed_ids))
870 @mock.patch.dict(os.environ, {"SATTLE_URI_BASE": "fake_host:1234"})
871 def test_filter_satellites_all_allowed(self):
872 science, matchedTemplate, difference, sources, detectionTask = self._setup_sattle_tests()
874 allowed_ids = list(range(1, 21))
875 response = MockResponse({"allow_list": allowed_ids}, 200, "all allowed")
876 # Run detection and check the results
877 with mock.patch('requests.put', return_value=response):
878 output = detectionTask.run(science, matchedTemplate, difference, sources,
879 idFactory=IdFactory.makeSimple())
881 # Output should be all sources that went in. 20 go in, 20 should come out
882 self.assertEqual(len(output.diaSources), 20)
884 self.assertEqual(set(output.diaSources['id']), set(allowed_ids))
886 @mock.patch.dict(os.environ, {"SATTLE_URI_BASE": "fake_host:1234"})
887 def test_filter_satellites_none_allowed(self):
888 science, matchedTemplate, difference, sources, detectionTask = self._setup_sattle_tests()
890 response = MockResponse({"allow_list": []}, 200, "none allowed")
891 # Run detection and confirm it raises for no diasources
892 with self.assertRaises(detectAndMeasure.NoDiaSourcesError):
893 with mock.patch('requests.put', return_value=response):
894 detectionTask.run(science, matchedTemplate, difference, sources,
895 idFactory=IdFactory.makeSimple())
897 @mock.patch.dict(os.environ, {"SATTLE_URI_BASE": ""})
898 def test_fail_on_sattle_misconfiguration(self):
899 with self.assertRaises(pexConfig.FieldValidationError):
900 self._setup_detection(run_sattle=True)
902 def test_trailed_glints(self):
903 """Test that the glint_trail column works, and that
904 the trailed_glints output contains the expected information.
905 """
906 noiseLevel = 1.
907 staticSeed = 1
908 diffim, diaSources = makeTestImage(seed=staticSeed, noiseLevel=noiseLevel, noiseSeed=6)
909 self._check_values(diaSources['glint_trail'])
911 # Run detection and return the output Struct so we can check it
912 def _detection_wrapper(diffim, diaSources):
913 detectionTask = self._setup_detection()
914 scienceBase, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6)
915 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7)
916 science = scienceBase.clone()
917 science.maskedImage -= diffim.maskedImage
918 difference = science.clone()
919 difference.maskedImage -= matchedTemplate.maskedImage
920 output = detectionTask.run(science, matchedTemplate, difference, sources)
921 return output
923 output = _detection_wrapper(diffim, diaSources)
924 self.assertTrue('slopes' in output.glintTrailInfo)
925 self.assertTrue('intercepts' in output.glintTrailInfo)
926 self.assertTrue('stderrs' in output.glintTrailInfo)
927 self.assertTrue('lengths' in output.glintTrailInfo)
928 self.assertTrue('angles' in output.glintTrailInfo)
931class DetectAndMeasureScoreTest(DetectAndMeasureTestBase, lsst.utils.tests.TestCase):
932 detectionTask = detectAndMeasure.DetectAndMeasureScoreTask
934 def test_detection_xy0(self):
935 """Basic functionality test with non-zero x0 and y0.
936 """
937 # Set up the simulated images
938 noiseLevel = 1.
939 staticSeed = 1
940 fluxLevel = 500
941 kernelSize = 31
942 psfSize = 2.4
943 # Buffer source positions away from the score image's EDGE strip.
944 # Preconvolution adds ~kernelSize//2 pixels of EDGE on top of the
945 # EDGE already set by detection smoothing on the science image.
946 templateBorderSize = kernelSize//2
947 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel, "x0": 12345, "y0": 67890,
948 "kernelSize": kernelSize, "templateBorderSize": templateBorderSize}
949 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
950 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
951 difference = science.clone()
952 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
953 scienceKernel = science.psf.getKernel()
954 score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
956 # Configure the detection Task
957 # Set the subtraction residual metric threshold high, since the template
958 # is not actually subtracted from the science image.
959 detectionTask = self._setup_detection(doDeblend=True, badSubtractionRatioThreshold=1.,
960 doSkySources=False)
962 # Run detection and check the results
963 output = detectionTask.run(science, matchedTemplate, difference, score, sources,
964 idFactory=self.idGenerator.make_table_id_factory())
966 # Catalog ids should be very large from this id generator.
967 self.assertTrue(all(output.diaSources['id'] > 1000000000))
968 subtractedMeasuredExposure = output.subtractedMeasuredExposure
970 self.assertImagesEqual(subtractedMeasuredExposure.image, difference.image)
972 self.assertEqual(len(output.diaSources), len(sources))
973 # no sources should be flagged as negative
974 self.assertEqual(len(~output.diaSources["is_negative"]), len(output.diaSources))
975 refIds = []
976 for source in sources:
977 self._check_diaSource(output.diaSources, source, refIds=refIds)
979 def test_measurements_finite(self):
980 """Measured fluxes and centroids should always be finite.
981 """
982 columnNames = ["coord_ra", "coord_dec", "ip_diffim_forced_PsfFlux_instFlux",
983 "ip_diffim_forced_template_PsfFlux_instFlux"]
985 # Set up the simulated images
986 noiseLevel = 1.
987 staticSeed = 1
988 transientSeed = 6
989 xSize = 256
990 ySize = 256
991 psfSize = 2.4
992 kwargs = {"psfSize": psfSize, "x0": 0, "y0": 0,
993 "xSize": xSize, "ySize": ySize}
994 science, sources = makeTestImage(seed=staticSeed, noiseLevel=noiseLevel, noiseSeed=6,
995 nSrc=1, **kwargs)
996 matchedTemplate, _ = makeTestImage(seed=staticSeed, noiseLevel=noiseLevel/4, noiseSeed=7,
997 nSrc=1, **kwargs)
998 rng = np.random.RandomState(3)
999 xLoc = np.arange(-5, xSize+5, 10)
1000 rng.shuffle(xLoc)
1001 yLoc = np.arange(-5, ySize+5, 10)
1002 rng.shuffle(yLoc)
1003 transients, transientSources = makeTestImage(seed=transientSeed,
1004 nSrc=len(xLoc), fluxLevel=1000.,
1005 noiseLevel=noiseLevel, noiseSeed=8,
1006 xLoc=xLoc, yLoc=yLoc,
1007 **kwargs)
1008 difference = science.clone()
1009 difference.maskedImage -= matchedTemplate.maskedImage
1010 difference.maskedImage += transients.maskedImage
1011 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
1012 scienceKernel = science.psf.getKernel()
1013 score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
1015 # Configure the detection Task
1016 detectionTask = self._setup_detection(doForcedMeasurement=True)
1018 # Run detection and check the results
1019 output = detectionTask.run(science, matchedTemplate, difference, score, sources)
1021 for column in columnNames:
1022 self._check_values(output.diaSources[column])
1023 self._check_values(output.diaSources.getX(), minValue=0, maxValue=xSize)
1024 self._check_values(output.diaSources.getY(), minValue=0, maxValue=ySize)
1025 self._check_values(output.diaSources.getPsfInstFlux())
1027 def test_detect_transients(self):
1028 """Run detection on a difference image containing transients.
1029 """
1030 # Set up the simulated images
1031 noiseLevel = 1.
1032 staticSeed = 1
1033 transientSeed = 6
1034 nTransients = 10
1035 transientFlux = 1000
1036 fluxLevel = 500
1037 psfSize = 2.4
1038 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel}
1039 scienceBase, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
1040 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
1041 scienceKernel = scienceBase.psf.getKernel()
1042 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
1044 # Configure the detection Task
1045 detectionTask = self._setup_detection(doMerge=False)
1046 kwargs["seed"] = transientSeed
1047 kwargs["nSrc"] = nTransients
1048 kwargs["fluxLevel"] = transientFlux
1050 # Run detection and check the results
1051 def _run_and_check_detections(positive=True):
1052 """Simulate positive or negative transients and run detection.
1054 Parameters
1055 ----------
1056 positive : `bool`, optional
1057 If set, use positive transient sources.
1058 """
1060 transients, transientSources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, **kwargs)
1061 science = scienceBase.clone()
1062 if positive:
1063 science.maskedImage += transients.maskedImage
1064 else:
1065 science.maskedImage -= transients.maskedImage
1066 difference = science.clone()
1067 difference.maskedImage -= matchedTemplate.maskedImage
1068 score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
1069 # NOTE: NoiseReplacer (run by forcedMeasurement) can modify the
1070 # science image if we've e.g. removed parents post-deblending.
1071 # Pass a clone of the science image, so that it doesn't disrupt
1072 # later tests.
1073 output = detectionTask.run(science, matchedTemplate, difference, score, sources)
1074 refIds = []
1075 scale = 1. if positive else -1.
1076 # sources near the edge may have untrustworthy centroids
1077 goodSrcFlags = ~output.diaSources['base_PixelFlags_flag_edge']
1078 for diaSource, goodSrcFlag in zip(output.diaSources, goodSrcFlags):
1079 if goodSrcFlag:
1080 self._check_diaSource(transientSources, diaSource, refIds=refIds, scale=scale)
1081 _run_and_check_detections(positive=True)
1082 _run_and_check_detections(positive=False)
1084 def test_detect_dipoles(self):
1085 """Run detection on a difference image containing dipoles.
1086 """
1087 # Set up the simulated images
1088 noiseLevel = 1.
1089 staticSeed = 1
1090 fluxLevel = 1000
1091 fluxRange = 1.5
1092 nSources = 10
1093 offset = 1
1094 xSize = 300
1095 ySize = 300
1096 kernelSize = 31
1097 psfSize = 2.4
1098 # Avoid placing sources near the edge for this test, so that we can
1099 # easily check that the correct number of sources are detected.
1100 templateBorderSize = kernelSize//2
1101 dipoleFlag = "ip_diffim_DipoleFit_classification"
1102 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel, "fluxRange": fluxRange,
1103 "nSrc": nSources, "templateBorderSize": templateBorderSize, "kernelSize": kernelSize,
1104 "xSize": xSize, "ySize": ySize}
1105 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
1106 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
1107 difference = science.clone()
1108 # Shift the template by a pixel in order to make dipoles in the difference image.
1109 matchedTemplate.image.array[...] = np.roll(matchedTemplate.image.array[...], offset, axis=0)
1110 matchedTemplate.variance.array[...] = np.roll(matchedTemplate.variance.array[...], offset, axis=0)
1111 matchedTemplate.mask.array[...] = np.roll(matchedTemplate.mask.array[...], offset, axis=0)
1112 difference.maskedImage -= matchedTemplate.maskedImage[science.getBBox()]
1113 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
1114 scienceKernel = science.psf.getKernel()
1115 score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
1117 detectionTask = self._setup_detection(badSubtractionRatioThreshold=0.3)
1118 output = detectionTask.run(science, matchedTemplate, difference, score, sources)
1119 diaSources = output.diaSources[~output.diaSources["sky_source"]].copy(deep=True)
1120 self.assertEqual(len(diaSources), len(sources))
1121 refIds = []
1122 for diaSource in diaSources:
1123 if diaSource[dipoleFlag]: 1123 ↛ 1130line 1123 didn't jump to line 1130 because the condition on line 1123 was always true
1124 self._check_diaSource(sources, diaSource, refIds=refIds, scale=0,
1125 rtol=0.05, atol=None, usePsfFlux=False)
1126 self.assertFloatsAlmostEqual(diaSource["ip_diffim_DipoleFit_orientation"],
1127 -np.pi / 2, atol=2.)
1128 self.assertFloatsAlmostEqual(diaSource["ip_diffim_DipoleFit_separation"], offset, rtol=0.1)
1129 else:
1130 raise ValueError("DiaSource with ID %s is not a dipole!", diaSource.getId())
1132 def test_sky_sources(self):
1133 """Add sky sources and check that they are sufficiently far from other
1134 sources and have negligible flux.
1135 """
1136 # Set up the simulated images
1137 noiseLevel = 1.
1138 staticSeed = 1
1139 transientSeed = 6
1140 transientFluxLevel = 1000.
1141 transientFluxRange = 1.5
1142 fluxLevel = 500
1143 psfSize = 2.4
1144 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel}
1145 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
1146 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
1147 transients, transientSources = makeTestImage(seed=transientSeed, psfSize=2.4,
1148 nSrc=10, fluxLevel=transientFluxLevel,
1149 fluxRange=transientFluxRange,
1150 noiseLevel=noiseLevel, noiseSeed=8)
1151 difference = science.clone()
1152 difference.maskedImage -= matchedTemplate.maskedImage
1153 difference.maskedImage += transients.maskedImage
1154 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
1155 scienceKernel = science.psf.getKernel()
1156 kernelWidth = np.max(scienceKernel.getDimensions())//2
1157 score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
1159 # Configure the detection Task
1160 detectionTask = self._setup_detection(doSkySources=True)
1162 # Run detection and check the results
1163 output = detectionTask.run(science, matchedTemplate, difference, score, sources,
1164 idFactory=self.idGenerator.make_table_id_factory())
1165 nSkySourcesGenerated = detectionTask.metadata["n_skySources"]
1166 skySources = output.diaSources[output.diaSources["sky_source"]]
1167 self.assertEqual(len(skySources), nSkySourcesGenerated)
1168 for skySource in skySources:
1169 # Disable the sky_source flag to enable checks.
1170 skySource["sky_source"] = False
1171 # The sky sources should not be close to any other source
1172 with self.assertRaises(AssertionError):
1173 self._check_diaSource(transientSources, skySource, matchDistance=kernelWidth)
1174 with self.assertRaises(AssertionError):
1175 self._check_diaSource(sources, skySource, matchDistance=kernelWidth)
1176 # The sky sources should have low flux levels.
1177 self._check_diaSource(transientSources, skySource, matchDistance=1000, scale=0.,
1178 atol=np.sqrt(transientFluxRange*transientFluxLevel))
1180 # Catalog ids should be very large from this id generator.
1181 self.assertTrue(all(output.diaSources['id'] > 1000000000))
1183 def test_exclude_mask_detections(self):
1184 """Sources with certain bad mask planes set should not be detected.
1185 """
1186 # Set up the simulated images
1187 noiseLevel = 1.
1188 staticSeed = 1
1189 transientSeed = 6
1190 fluxLevel = 500
1191 radius = 2
1192 psfSize = 2.4
1193 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel}
1194 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
1195 science.getInfo().setVisitInfo(makeVisitInfo())
1196 detector = DetectorWrapper(numAmps=1).detector
1197 science.setDetector(detector)
1198 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
1200 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
1201 scienceKernel = science.psf.getKernel()
1202 # Configure the detection Task
1203 badSourceFlags = ["base_PixelFlags_flag_offimage",
1204 "base_PixelFlags_flag_edgeCenterAll",
1205 "base_PixelFlags_flag_badCenterAll",
1206 "base_PixelFlags_flag_saturatedCenterAll", ]
1207 detectionTask = self._setup_detection(nSkySources=0, badSourceFlags=badSourceFlags)
1208 excludeMaskPlanes = ["EDGE", "BAD", "SAT"]
1209 nBad = len(excludeMaskPlanes)
1210 self.assertGreater(nBad, 0)
1211 kwargs["seed"] = transientSeed
1212 kwargs["nSrc"] = nBad
1213 kwargs["fluxLevel"] = 1000
1215 # Run detection and check the results
1216 def _run_and_check_detections(setFlags=True):
1217 transients, transientSources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, **kwargs)
1218 difference = science.clone()
1219 difference.maskedImage -= matchedTemplate.maskedImage
1220 difference.maskedImage += transients.maskedImage
1221 if setFlags:
1222 for src, badMask in zip(transientSources, excludeMaskPlanes):
1223 srcX = int(src.getX())
1224 srcY = int(src.getY())
1225 srcBbox = lsst.geom.Box2I(lsst.geom.Point2I(srcX - radius, srcY - radius),
1226 lsst.geom.Extent2I(2*radius + 1, 2*radius + 1))
1227 difference[srcBbox].mask.array |= lsst.afw.image.Mask.getPlaneBitMask(badMask)
1228 score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
1229 if setFlags:
1230 with self.assertRaises(AlgorithmError):
1231 output = detectionTask.run(science, matchedTemplate, difference, score, sources)
1232 else:
1233 output = detectionTask.run(science, matchedTemplate, difference, score, sources)
1234 refIds = []
1235 goodSrcFlags = checkMask(difference.mask, transientSources, excludeMaskPlanes)
1236 self.assertEqual(np.sum(~goodSrcFlags), 0)
1237 for diaSource, goodSrcFlag in zip(output.diaSources, goodSrcFlags):
1238 if ~goodSrcFlag: 1238 ↛ 1239line 1238 didn't jump to line 1239 because the condition on line 1238 was never true
1239 with self.assertRaises(AssertionError):
1240 self._check_diaSource(transientSources, diaSource, refIds=refIds)
1241 else:
1242 self._check_diaSource(transientSources, diaSource, refIds=refIds)
1243 _run_and_check_detections(setFlags=False)
1244 _run_and_check_detections(setFlags=True)
1246 def test_detect_transients_with_background(self):
1247 """Background measured on the difference image should be subtracted
1248 from the score image so that transients are still detected cleanly.
1249 """
1250 # Set up the simulated images
1251 noiseLevel = 1.
1252 staticSeed = 1
1253 transientSeed = 6
1254 nTransients = 10
1255 transientFlux = 1000
1256 fluxLevel = 500
1257 xSize = 512
1258 ySize = 512
1259 x0 = 123
1260 y0 = 456
1261 kernelSize = 31
1262 psfSize = 2.4
1263 # Avoid placing sources near the edge for this test, so that we can
1264 # easily check that the correct number of sources are detected.
1265 templateBorderSize = kernelSize//2
1266 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel,
1267 "xSize": xSize, "ySize": ySize, "x0": x0, "y0": y0,
1268 "kernelSize": kernelSize, "templateBorderSize": templateBorderSize}
1269 chebyshev_params = [2.2, 2.1, 2.0, 1.2, 1.1, 1.0]
1271 # The background model must cover the grown image bbox, otherwise it
1272 # would be extrapolated outside its declared domain.
1273 bbox2D = lsst.geom.Box2D(
1274 lsst.geom.Point2D(x0 - templateBorderSize, y0 - templateBorderSize),
1275 lsst.geom.Extent2D(xSize + 2*templateBorderSize, ySize + 2*templateBorderSize))
1276 background_model = afwMath.Chebyshev1Function2D(chebyshev_params, bbox2D)
1277 scienceBase, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6,
1278 background=background_model, **kwargs)
1279 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
1280 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
1281 scienceKernel = scienceBase.psf.getKernel()
1283 # Configure the detection Task
1284 detectionTask = self._setup_detection(doMerge=False, doSubtractBackground=True)
1285 kwargs["seed"] = transientSeed
1286 kwargs["nSrc"] = nTransients
1287 kwargs["fluxLevel"] = transientFlux
1289 def _run_and_check_detections(positive=True):
1290 transients, transientSources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, **kwargs)
1291 science = scienceBase.clone()
1292 if positive:
1293 science.maskedImage += transients.maskedImage
1294 else:
1295 science.maskedImage -= transients.maskedImage
1296 difference = science.clone()
1297 difference.maskedImage -= matchedTemplate.maskedImage
1298 score = subtractTask._convolveExposure(difference, scienceKernel,
1299 subtractTask.convolutionControl)
1300 # Record the score image metric before detection so we can verify
1301 # that the background was actually subtracted from it.
1302 originalScoreMedian = np.median(score.image.array[np.isfinite(score.image.array)])
1303 originalDiffimMedian = np.median(difference.image.array[np.isfinite(difference.image.array)])
1304 output = detectionTask.run(science, matchedTemplate, difference, score, sources)
1306 # The background subtraction should leave a non-empty BackgroundList.
1307 self.assertGreater(len(output.differenceBackground), 0)
1308 # The difference image should have been modified in place by the
1309 # background subtraction, so its median should shift noticeably.
1310 subtractedScoreMedian = np.median(score.image.array[np.isfinite(score.image.array)])
1311 self.assertLess(np.abs(subtractedScoreMedian), np.abs(originalScoreMedian))
1312 subtractedDiffimMedian = np.median(difference.image.array[np.isfinite(difference.image.array)])
1313 self.assertLess(np.abs(subtractedDiffimMedian), np.abs(originalDiffimMedian))
1315 refIds = []
1316 scale = 1. if positive else -1.
1317 for transient in transientSources:
1318 self._check_diaSource(output.diaSources, transient, refIds=refIds, scale=scale)
1319 _run_and_check_detections(positive=True)
1320 _run_and_check_detections(positive=False)
1322 def test_background_detection_mask_not_kept(self):
1323 """The detection masks of the difference image should match those of
1324 the score image, without the detections used to fit the background.
1325 """
1326 kwargs = {"seed": 1, "psfSize": 2.4, "fluxLevel": 500, "kernelSize": 31,
1327 "templateBorderSize": 15}
1328 science, sources = makeTestImage(noiseLevel=1., noiseSeed=6, **kwargs)
1329 matchedTemplate, _ = makeTestImage(noiseLevel=0.25, noiseSeed=7, **kwargs)
1330 difference = science.clone()
1331 # Small-scale structure that the first-pass background fit removes
1332 # but the final fit does not, so the two passes detect different
1333 # pixels.
1334 difference.image.array += 0.3*np.sin(np.arange(difference.getWidth())/20.)
1335 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
1336 score = subtractTask._convolveExposure(difference, science.psf.getKernel(),
1337 subtractTask.convolutionControl)
1338 detectionTask = self._setup_detection(doSubtractBackground=True, doSkySources=False,
1339 badSubtractionRatioThreshold=1.)
1341 output = detectionTask.run(science, matchedTemplate, difference, score, sources)
1343 differenceMask = output.subtractedMeasuredExposure.mask
1344 scoreMask = output.scoreMeasuredExposure.mask
1345 for plane in ("DETECTED", "DETECTED_NEGATIVE"):
1346 with self.subTest(plane=plane):
1347 differencePix = differenceMask.array & differenceMask.getPlaneBitMask(plane) > 0
1348 scorePix = scoreMask.array & scoreMask.getPlaneBitMask(plane) > 0
1349 self.assertGreater(scorePix.sum(), 0)
1350 np.testing.assert_array_equal(differencePix, scorePix)
1352 def test_mask_cosmic_rays(self):
1353 """Cosmic rays detected on the difference image should propagate
1354 to the mask of the returned (measured) exposure. This version creates
1355 and uses a Score image for detection.
1356 """
1357 # Set up the simulated images
1358 noiseLevel = 1.
1359 staticSeed = 1
1360 fluxLevel = 500
1361 xSize = 400
1362 ySize = 400
1363 psfSize = 2.4
1364 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel,
1365 "xSize": xSize, "ySize": ySize}
1366 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
1367 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
1368 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
1369 scienceKernel = science.psf.getKernel()
1370 crMask = science.mask.getPlaneBitMask("CR")
1372 # Configure the detection Task
1373 detectionTask = self._setup_detection(doMerge=False, doSkySources=True)
1375 # Test that no CRs are detected when none are present.
1376 transients, _ = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, nSrc=10, **kwargs)
1377 science.maskedImage += transients.maskedImage
1378 difference = science.clone()
1379 difference.maskedImage -= matchedTemplate.maskedImage
1380 score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
1381 output = detectionTask.run(science, matchedTemplate, difference, score, sources)
1382 self.assertFalse(np.any((output.subtractedMeasuredExposure.mask.array & crMask) > 0))
1384 crX0 = round(sources.getX()[0] - science.getX0())
1385 crY0 = round(sources.getY()[0] - science.getY0())
1386 crX1 = round(sources.getX()[1] - science.getX0())
1387 crY1 = round(sources.getY()[1] - science.getY0())
1388 # Inject CR-like shapes into the difference image. CR detection runs
1389 # on the difference, and the mask should propagate to the measured
1390 # exposure returned by the task.
1391 crMaskSetInput = np.zeros(difference.image.array.shape, bool)
1392 crMaskSetInput[crX0:crX0+1, crY0:crY0+5] = True
1393 crMaskSetInput[crX1:crX1+5, crY1:crY1+1] = True
1394 difference.image.array[crMaskSetInput] += 1234
1395 score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
1396 output = detectionTask.run(science, matchedTemplate, difference, score, sources)
1397 crMaskSet = (output.subtractedMeasuredExposure.mask.array & crMask) > 0
1398 self.assertFalse(np.any(crMaskSet[~crMaskSetInput]))
1399 self.assertTrue(np.all(crMaskSet[crMaskSetInput]))
1402class TestNegativePeaks(lsst.utils.tests.TestCase):
1403 """Tests of deblending and merging negative peaks, to test fixes for the
1404 various problems found on DM-48596.
1405 """
1407 def testDeblendNegatives(self):
1408 """Test that negative peaks get deblended and not destroyed: DM-48704.
1409 This is only a test of deblending, not of merging.
1410 """
1411 # Make a science image with one blend of two positive sources, to
1412 # subtract from an empty template, resulting in a negative diffim.
1413 bbox = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(100, 100))
1414 dataset = lsst.meas.base.tests.TestDataset(bbox)
1415 delta = 10
1416 with dataset.addBlend() as family:
1417 family.addChild(instFlux=2E5, centroid=lsst.geom.Point2D(50, 72))
1418 family.addChild(instFlux=2.5E5, centroid=lsst.geom.Point2D(50+delta, 74))
1419 science, catalog = dataset.realize(noise=1.0,
1420 schema=lsst.meas.base.tests.TestDataset.makeMinimalSchema())
1421 dataset = lsst.meas.base.tests.TestDataset(bbox)
1422 template, _ = dataset.realize(noise=1.0,
1423 schema=lsst.meas.base.tests.TestDataset.makeMinimalSchema())
1424 difference = template.clone()
1425 difference.image -= science.image
1427 config = detectAndMeasure.DetectAndMeasureTask.ConfigClass()
1428 config.doDeblend = True
1429 task = detectAndMeasure.DetectAndMeasureTask(config=config)
1430 # prelude steps taken from `detectAndMeasure.run`
1431 task._prepareInputs(difference)
1432 table = lsst.afw.table.SourceTable.make(task.schema)
1433 results = task.detection.run(table=table, exposure=difference, doSmooth=True)
1434 # Just run the deblend step so we can check the footprints independently.
1435 sources, positives, negatives = task._deblend(difference, results.positive, results.negative)
1437 # DM-48704 fixed a problem where the peaks were in the footprints, but
1438 # the spans were empty.
1439 self.assertEqual(len(negatives), 2)
1440 self.assertEqual(len(positives), 0)
1441 self.assertGreater(negatives[0].getFootprint().getSpans().getArea(), 0)
1442 self.assertGreater(negatives[1].getFootprint().getSpans().getArea(), 0)
1443 # Deblended children are HeavyFootprints; we have to make sure the
1444 # pixel values in those are correct (though DetectAndMeasureTask
1445 # doesn't use the fact that they're Heavy).
1446 # The sources are positive in the science image, and negative in the
1447 # diffim, so the minimum value in the deblended negative footprint is
1448 # the maximum value in the science catalog footprint (ignoring the
1449 # noise in the science and template, hence rtol).
1450 # (catalog[0] is the parent; we want the children)
1451 self.assertFloatsAlmostEqual(negatives[0].getFootprint().getImageArray().min(),
1452 -catalog[1].getFootprint().getImageArray().max(), rtol=1e-4)
1453 self.assertFloatsAlmostEqual(negatives[1].getFootprint().getImageArray().min(),
1454 -catalog[2].getFootprint().getImageArray().max(), rtol=1e-4)
1456 self.assertEqual(sources["is_negative"].sum(), 2)
1458 def testMergeFootprints(self):
1459 """Test that merging footprints does not cause negative ones to
1460 disappear (e.g. get merged into non-connected footprints).
1462 As implemented, the diffim will have three positive sources (one a
1463 blended pair), and 7 negative sources (two a blended pair).
1464 """
1465 # Make a template image with multiple blends, designed to trigger the
1466 # negative-footprint-related bug in `FootprintSet.merge`.
1467 bbox = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(400, 200))
1468 dataset = lsst.meas.base.tests.TestDataset(bbox)
1469 # Because these sources are on the template, they will be negative on
1470 # the diffim, unless the pixels are explicitly made negative via
1471 # `template.image.subset` below.
1472 # positive isolated source on diffim
1473 dataset.addSource(instFlux=.5E5, centroid=lsst.geom.Point2D(25, 26))
1474 # negative isolated source on diffim
1475 dataset.addSource(instFlux=.7E5, centroid=lsst.geom.Point2D(75, 24),
1476 shape=lsst.afw.geom.Quadrupole(12, 7, 2))
1477 delta = 10
1478 # negative blended pair on diffim
1479 with dataset.addBlend() as family:
1480 family.addChild(instFlux=1E5, centroid=lsst.geom.Point2D(150, 72))
1481 family.addChild(instFlux=1.5E5, centroid=lsst.geom.Point2D(150+delta, 74))
1482 # positive blended pair on diffim
1483 with dataset.addBlend() as family:
1484 family.addChild(instFlux=2E5, centroid=lsst.geom.Point2D(250, 72))
1485 family.addChild(instFlux=2.5E5, centroid=lsst.geom.Point2D(250+delta, 74))
1486 # negative blended pair on diffim
1487 with dataset.addBlend() as family:
1488 family.addChild(instFlux=3E5, centroid=lsst.geom.Point2D(350, 72))
1489 family.addChild(instFlux=3.5E5, centroid=lsst.geom.Point2D(350+delta, 74))
1490 # negative isolated source on diffim
1491 dataset.addSource(instFlux=4E5, centroid=lsst.geom.Point2D(75, 124))
1492 # negative isolated source on diffim
1493 dataset.addSource(instFlux=5E5, centroid=lsst.geom.Point2D(175, 124))
1494 schema = lsst.meas.base.tests.TestDataset.makeMinimalSchema()
1495 schema.addField("sky_source", "Flag", "Sky source.")
1496 template, catalog = dataset.realize(noise=100.0, schema=schema)
1497 # These will be positive sources on the diffim (as noted above).
1498 template.image.subset(lsst.geom.Box2I(lsst.geom.Point2I(15, 15),
1499 lsst.geom.Point2I(35, 35))).array *= -1
1500 template.image.subset(lsst.geom.Box2I(lsst.geom.Point2I(230, 60),
1501 lsst.geom.Point2I(270, 85))).array *= -1
1502 dataset = lsst.meas.base.tests.TestDataset(bbox)
1503 science, _ = dataset.realize(noise=100.0, schema=schema)
1504 difference = science.clone()
1505 difference.image -= template.image
1507 config = detectAndMeasure.DetectAndMeasureTask.ConfigClass()
1508 config.doDeblend = True
1509 config.raiseOnBadSubtractionRatio = False
1510 config.raiseOnNoDiaSources = False
1511 task = detectAndMeasure.DetectAndMeasureTask(config=config)
1512 result = task.run(science, template, difference, catalog)
1514 # The original bug manifested as the (175,124) source disappaering in
1515 # the merge step, and being included in the peaks of a duplicated
1516 # (75,124) source, even though their footprints are not connected.
1517 mask = np.isclose(result.diaSources["slot_Centroid_x"], 75) & \
1518 np.isclose(result.diaSources["slot_Centroid_y"], 124)
1519 self.assertEqual(mask.sum(), 1)
1520 peaks = result.diaSources[mask][0].getFootprint().peaks
1521 self.assertEqual(len(peaks), 1)
1522 self.assertEqual(peaks[0].getIx(), 75)
1523 self.assertEqual(peaks[0].getIy(), 124)
1525 mask = np.isclose(result.diaSources["slot_Centroid_x"], 175) & \
1526 np.isclose(result.diaSources["slot_Centroid_y"], 124)
1527 self.assertEqual(mask.sum(), 1)
1528 peaks = result.diaSources[mask][0].getFootprint().peaks
1529 self.assertEqual(len(peaks), 1)
1530 self.assertEqual(peaks[0].getIx(), 175)
1531 self.assertEqual(peaks[0].getIy(), 124)
1533 # This checks that all the returned footprints have exactly one peak;
1534 # it's a more general test than the above, but I'm keeping that for
1535 # the more explicit test of the original bugged sources.
1536 peak_x = np.array([peak['f_x'] for src in result.diaSources for peak in src.getFootprint().peaks])
1537 peak_y = np.array([peak['f_y'] for src in result.diaSources for peak in src.getFootprint().peaks])
1538 peaks = np.column_stack((peak_x, peak_y))
1539 unique_peaks, counts = np.unique(peaks, axis=0, return_counts=True)
1540 self.assertEqual(np.sum(counts > 1), 0)
1542 # There are three positive sources that should not have `is_negative`
1543 # set, independent of how deblending/merging of negative sources is
1544 # handled.
1545 self.assertEqual((~result.diaSources["is_negative"]).sum(), 3)
1547 def testReorderPeaksBySignificance(self):
1548 """A merged footprint whose most significant peak is negative should
1549 have that peak first after re-ordering, even though the merge lists
1550 positive peaks before negative ones.
1551 """
1552 config = detectAndMeasure.DetectAndMeasureTask.ConfigClass()
1553 task = detectAndMeasure.DetectAndMeasureTask(config=config)
1555 # Peak schema with a significance field, as SourceDetectionTask
1556 # produces for pixel_stdev thresholds.
1557 mapper = afwTable.SchemaMapper(afwDetection.PeakTable.makeMinimalSchema())
1558 mapper.addMinimalSchema(afwDetection.PeakTable.makeMinimalSchema())
1559 mapper.addOutputField("significance", type=float, doc="")
1560 peakSchema = mapper.getOutputSchema()
1562 diaSources = afwTable.SourceCatalog(task.schema)
1563 source = diaSources.addNew()
1564 spans = afwGeom.SpanSet(geom.Box2I(geom.Point2I(0, 0), geom.Extent2I(10, 10)))
1565 footprint = afwDetection.Footprint(spans, peakSchema)
1566 # Merge order is [positive, negative]: a weak positive peak first,
1567 # then a more significant negative peak.
1568 footprint.addPeak(1, 1, 30.0)
1569 footprint.peaks[-1]["significance"] = 5.0
1570 footprint.addPeak(2, 2, -500.0)
1571 footprint.peaks[-1]["significance"] = 40.0
1572 source.setFootprint(footprint)
1574 task._reorderPeaksBySignificance(diaSources)
1576 peaks = source.getFootprint().peaks
1577 self.assertEqual(peaks[0].getPeakValue(), -500.0)
1578 self.assertEqual(peaks[0]["significance"], 40.0)
1579 self.assertEqual(peaks[1].getPeakValue(), 30.0)
1580 self.assertEqual(peaks[1]["significance"], 5.0)
1583def makeVisitInfo(id=1):
1584 """Return a non-NaN visitInfo."""
1585 return afwImage.VisitInfo(id=id,
1586 exposureTime=10.01,
1587 darkTime=11.02,
1588 date=dafBase.DateTime(65321.1, dafBase.DateTime.MJD, dafBase.DateTime.TAI),
1589 ut1=12345.1,
1590 era=45.1*geom.degrees,
1591 boresightRaDec=geom.SpherePoint(23.1, 73.2, geom.degrees),
1592 boresightAzAlt=geom.SpherePoint(134.5, 33.3, geom.degrees),
1593 boresightAirmass=1.73,
1594 boresightRotAngle=73.2*geom.degrees,
1595 rotType=afwImage.RotType.SKY,
1596 observatory=Observatory(
1597 11.1*geom.degrees, 22.2*geom.degrees, 0.333),
1598 weather=Weather(1.1, 2.2, 34.5),
1599 )
1602class DetectAndMeasureOutputFormatTest(DetectAndMeasureTestBase, lsst.utils.tests.TestCase):
1603 """Tests of the conversion of the output images to `lsst.images` types.
1604 """
1605 detectionTask = detectAndMeasure.DetectAndMeasureTask
1607 @classmethod
1608 def setUpClass(cls):
1609 super().setUpClass()
1610 # `DifferenceImage.from_legacy` needs a detector with per-amplifier
1611 # raw geometry and a FIELD_ANGLE transform; the simpler
1612 # `DetectorWrapper` detector used elsewhere in this file does not
1613 # have either, and fails to convert.
1614 cls.detector = list(CameraWrapper().camera)[0]
1616 def setUp(self):
1617 self.dataId = lsst.daf.butler.DataCoordinate.standardize(
1618 instrument="I",
1619 visit=42,
1620 detector=12,
1621 universe=lsst.daf.butler.DimensionUniverse(),
1622 )
1623 # The exposure record that `runQuantum` rebuilds from the
1624 # observation metadata of the input difference image.
1625 self.exposureRecord = makeTestExposureRecord(self.dataId.universe,
1626 instrument=self.dataId["instrument"],
1627 visit=self.dataId["visit"])
1629 @staticmethod
1630 def _makeTemplateInfo(patch=42):
1631 """Return one record of a coadd that went into a template.
1633 Parameters
1634 ----------
1635 patch : `int`, optional
1636 Identifier of the patch the record describes.
1638 Returns
1639 -------
1640 info : `lsst.images.DifferenceImageTemplateInfo`
1641 The record.
1642 """
1643 return DifferenceImageTemplateInfo(
1644 skymap="skymap", tract=9813, patch=patch, dataset_id=uuid.uuid4(),
1645 dataset_run="a/template/run", bounds=Polygon.from_box(Box.factory[0:50, 0:50]),
1646 psf_shape_xx=4.0, psf_shape_yy=4.0, psf_shape_xy=0.0, psf_shape_flag=False)
1648 @staticmethod
1649 def _makeKernel(scale=1.0):
1650 """Return a PSF matching kernel of the one form that can be attached
1651 to a difference image.
1653 Parameters
1654 ----------
1655 scale : `float`, optional
1656 Coefficient of the first basis kernel, to tell two kernels apart.
1658 Returns
1659 -------
1660 kernel : `lsst.images.convolution_kernels.ImageBasisConvolutionKernel`
1661 The kernel.
1662 """
1663 basis = []
1664 for sigma in (1.0, 2.0):
1665 image = afwImage.ImageD(5, 5)
1666 y, x = np.mgrid[-2:3, -2:3]
1667 image.array[:] = np.exp(-(x**2 + y**2)/(2*sigma**2))
1668 basis.append(afwMath.FixedKernel(image))
1669 bbox = lsst.geom.Box2D(lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Extent2I(100, 100)))
1670 legacy = afwMath.LinearCombinationKernel(basis, afwMath.Chebyshev1Function2D(1, bbox))
1671 legacy.setSpatialParameters([[scale, 0.0, 0.0], [0.5, 0.0, 0.0]])
1672 return ImageBasisConvolutionKernel.from_legacy(legacy)
1674 @staticmethod
1675 def _makePhotometricScaling(exposure):
1676 """Return a constant photometric scaling for an input image.
1678 Parameters
1679 ----------
1680 exposure : `lsst.afw.image.ExposureF`
1681 Image from `_make_images`.
1683 Returns
1684 -------
1685 photometricScaling : `lsst.images.fields.ChebyshevField`
1686 The scaling.
1687 """
1688 return ChebyshevField(Box.from_legacy(exposure.getBBox()), np.array([[0.25]]),
1689 unit=astropy.units.nJy/astropy.units.electron)
1691 def _makeFutureInput(self, exposure, templateInfo=None, kernel=None, photometricScaling=None):
1692 """Convert a legacy image to a `lsst.images.DifferenceImage`.
1694 Parameters
1695 ----------
1696 exposure : `lsst.afw.image.ExposureF`
1697 Image from `_make_images`.
1698 templateInfo : `list` [`lsst.images.DifferenceImageTemplateInfo`],
1699 optional
1700 Record of the coadds that went into the template.
1701 kernel : `lsst.images.convolution_kernels.ImageBasisConvolutionKernel`,
1702 optional
1703 PSF matching kernel.
1704 photometricScaling : `lsst.images.fields.BaseField`, optional
1705 Photometric scaling of the science image.
1707 Returns
1708 -------
1709 image : `lsst.images.DifferenceImage`
1710 The converted image.
1711 """
1712 image = DifferenceImage.from_legacy(
1713 exposure, unit=astropy.units.nJy, exposure_record=self.exposureRecord)
1714 if templateInfo is not None:
1715 image.templates = templateInfo
1716 if kernel is not None:
1717 image.kernel = kernel
1718 if photometricScaling is not None:
1719 image.photometric_scaling = photometricScaling
1720 return image
1722 def _make_images(self, withScore=False):
1723 """Make a science/template/difference image set for detection.
1725 The template is not actually subtracted from the science image, so
1726 that there are plenty of detections; configure the detection task
1727 with ``badSubtractionRatioThreshold=1.`` to match.
1729 Parameters
1730 ----------
1731 withScore : `bool`, optional
1732 Also make and return a score (maximum likelihood) image?
1734 Returns
1735 -------
1736 images : `tuple`
1737 The science exposure, the source catalog, the matched template,
1738 the difference image, and (if ``withScore``) the score image.
1739 """
1740 kernelSize = 31
1741 kwargs = {"seed": 1, "psfSize": 2.4, "fluxLevel": 500, "x0": 12345, "y0": 67890,
1742 "kernelSize": kernelSize, "templateBorderSize": kernelSize//2}
1743 science, sources = makeTestImage(noiseLevel=1., noiseSeed=6, **kwargs)
1744 matchedTemplate, _ = makeTestImage(noiseLevel=0.25, noiseSeed=7, **kwargs)
1745 # Add components to the basic test image
1746 science.setDetector(self.detector)
1747 matchedTemplate.setDetector(self.detector)
1748 for exposure in (science, matchedTemplate):
1749 exposure.info.setVisitInfo(makeTestVisitInfo(self.dataId["visit"]))
1750 difference = science.clone()
1751 if not withScore:
1752 return science, sources, matchedTemplate, difference
1753 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
1754 score = subtractTask._convolveExposure(difference, science.psf.getKernel(),
1755 subtractTask.convolutionControl)
1756 score.setDetector(self.detector)
1757 return science, sources, matchedTemplate, difference, score
1759 def _setup_task(self, **kwargs):
1760 """Configure a detection task to go with `_make_images`.
1762 Parameters
1763 ----------
1764 **kwargs
1765 Any additional config parameters to set.
1767 Returns
1768 -------
1769 task : `lsst.pipe.base.PipelineTask`
1770 The configured task to use for detection and measurement.
1771 """
1772 # The template is not actually subtracted from the science image in
1773 # `_make_images`, so the subtraction residual metric threshold has to
1774 # be set high.
1775 kwargs.setdefault("badSubtractionRatioThreshold", 1.)
1776 return self._setup_detection(doSkySources=False, **kwargs)
1778 def _check_converted(self, image, legacy):
1779 """Check the properties common to every converted output image.
1781 Parameters
1782 ----------
1783 image : `lsst.images.DifferenceImage`
1784 The converted image.
1785 legacy : `lsst.afw.image.ExposureF`
1786 The legacy output exposure it was converted from.
1787 """
1788 self.assertIsInstance(image, DifferenceImage)
1789 self.assertEqual(image.unit, astropy.units.nJy)
1790 self.assertIsNotNone(image.detector)
1791 # `_makeFutureInput` was not given a photometric scaling.
1792 self.assertIsNone(image.photometric_scaling)
1793 compare_masked_image_to_legacy(image, legacy.maskedImage,
1794 plane_map=get_legacy_difference_image_mask_planes())
1795 self.assertTrue(np.any(image.mask.get("DETECTED")))
1797 def test_connections(self):
1798 """Check that the output storage classes follow ``image_type``,
1799 and that the background is a separate output only in legacy mode.
1800 """
1801 for configClass, connectionsClass in (
1802 (detectAndMeasure.DetectAndMeasureConfig, detectAndMeasure.DetectAndMeasureConnections),
1803 (detectAndMeasure.DetectAndMeasureScoreConfig,
1804 detectAndMeasure.DetectAndMeasureScoreConnections),
1805 ):
1806 for imageType, storageClass in (("legacy", "ExposureF"), ("future", "DifferenceImage")):
1807 with self.subTest(configClass=configClass.__name__, imageType=imageType):
1808 config = configClass()
1809 config.image_type = imageType
1810 config.doSubtractBackground = True
1811 with self.assertWarns(FutureWarning):
1812 config.doWriteBackground = True
1813 connections = connectionsClass(config=config)
1814 self.assertEqual(connections.subtractedMeasuredExposure.storageClass, storageClass)
1815 if connectionsClass is detectAndMeasure.DetectAndMeasureScoreConnections:
1816 self.assertEqual(connections.scoreMeasuredExposure.storageClass, storageClass)
1817 if imageType == "legacy":
1818 self.assertIn("differenceBackground", connections.outputs)
1819 self.assertEqual(connections.differenceBackground.storageClass, "Background")
1820 else:
1821 # The background is recorded on the output images.
1822 self.assertNotIn("differenceBackground", connections.outputs)
1824 def test_run_returns_legacy_types(self):
1825 """``run`` always returns legacy types, even in future mode."""
1826 science, sources, matchedTemplate, difference = self._make_images()
1827 task = self._setup_task(image_type="future")
1828 output = task.run(science, matchedTemplate, difference, sources)
1829 self.assertIsInstance(output.subtractedMeasuredExposure, afwImage.ExposureF)
1831 def test_legacy_mode_does_not_convert(self):
1832 """Legacy mode leaves the outputs unconverted."""
1833 science, sources, matchedTemplate, difference = self._make_images()
1834 task = self._setup_task(doSubtractBackground=True, image_type="legacy")
1835 output = task.run(science, matchedTemplate, difference, sources)
1836 self.assertIsInstance(output.subtractedMeasuredExposure, afwImage.ExposureF)
1837 self.assertIsInstance(output.differenceBackground, afwMath.BackgroundList)
1838 lsst.pipe.base.testUtils.assertValidOutput(task, output)
1840 def test_convert_outputs_to_future(self):
1841 """Convert the outputs of ``DetectAndMeasureTask``."""
1842 science, sources, matchedTemplate, difference = self._make_images()
1843 task = self._setup_task(doSubtractBackground=True, image_type="future")
1844 inputImage = self._makeFutureInput(difference)
1845 output = task.run(science, matchedTemplate, inputImage.to_legacy(), sources)
1846 legacy = output.subtractedMeasuredExposure
1847 background = output.differenceBackground
1848 task.convert_outputs_to_future(output, {"subtractedMeasuredExposure": inputImage})
1850 self.assertIs(output.subtractedMeasuredExposure, inputImage)
1851 self._check_converted(inputImage, legacy)
1852 self.assertEqual(list(inputImage.backgrounds), ["subtracted"])
1853 # The separate background dataset is still written in future mode.
1854 self.assertIs(output.differenceBackground, background)
1855 self.assertIsInstance(output.differenceBackground, afwMath.BackgroundList)
1857 def test_convert_outputs_to_future_copied_pixels(self):
1858 """The output pixels are copied into the input image when the legacy
1859 image that was measured does not share them.
1860 """
1861 science, sources, matchedTemplate, difference = self._make_images()
1862 task = self._setup_task(doSubtractBackground=True, image_type="future")
1863 inputImage = self._makeFutureInput(difference)
1864 original = inputImage.image.array.copy()
1865 output = task.run(science, matchedTemplate, inputImage.to_legacy(copy=True), sources)
1866 legacy = output.subtractedMeasuredExposure
1867 # Background subtraction changed the pixels of the copy only.
1868 self.assertFalse(np.array_equal(legacy.image.array, original))
1869 np.testing.assert_array_equal(inputImage.image.array, original)
1871 task.convert_outputs_to_future(output, {"subtractedMeasuredExposure": inputImage})
1873 self._check_converted(inputImage, legacy)
1875 def test_convert_outputs_to_future_source_injection(self):
1876 """Convert a difference image from a pipeline that injects sources.
1878 The injection planes are optional, so the conversion adds them for
1879 an image that has pixels set in them instead of raising.
1880 """
1881 planes = ("INJECTED", "INJECTED_CORE", "INJECTED_TEMPLATE", "INJECTED_CORE_TEMPLATE")
1882 science, sources, matchedTemplate, difference = self._make_images()
1883 task = self._setup_task(image_type="future")
1884 inputImage = self._makeFutureInput(difference)
1885 output = task.run(science, matchedTemplate, inputImage.to_legacy(), sources)
1886 legacy = output.subtractedMeasuredExposure
1887 for n, plane in enumerate(planes):
1888 legacy.mask.addMaskPlane(plane)
1889 legacy.mask.array[0, n] |= legacy.mask.getPlaneBitMask(plane)
1891 task.convert_outputs_to_future(output, {"subtractedMeasuredExposure": inputImage})
1893 image = output.subtractedMeasuredExposure
1894 self.assertIsInstance(image, DifferenceImage)
1895 # Name the optional planes in the map, so that dropping one fails
1896 # the comparison instead of skipping it.
1897 planeMap = get_legacy_difference_image_mask_planes() | get_legacy_optional_mask_planes()
1898 compare_masked_image_to_legacy(image, legacy.maskedImage, plane_map=planeMap)
1900 def test_convert_outputs_to_future_score(self):
1901 """Convert the outputs of ``DetectAndMeasureScoreTask``."""
1902 self.detectionTask = detectAndMeasure.DetectAndMeasureScoreTask
1903 science, sources, matchedTemplate, difference, score = self._make_images(withScore=True)
1904 task = self._setup_task(doSubtractBackground=True, image_type="future")
1905 inputImages = {"subtractedMeasuredExposure": self._makeFutureInput(difference),
1906 "scoreMeasuredExposure": self._makeFutureInput(score)}
1907 output = task.run(science, matchedTemplate,
1908 inputImages["subtractedMeasuredExposure"].to_legacy(),
1909 inputImages["scoreMeasuredExposure"].to_legacy(), sources)
1910 legacyOutputs = {name: getattr(output, name) for name in inputImages}
1911 task.convert_outputs_to_future(output, inputImages)
1913 for name, legacy in legacyOutputs.items():
1914 with self.subTest(name=name):
1915 self.assertIs(getattr(output, name), inputImages[name])
1916 self._check_converted(getattr(output, name), legacy)
1917 # The same background is subtracted from both images.
1918 self.assertIn("subtracted", getattr(output, name).backgrounds)
1920 def test_convert_outputs_to_future_empty_background(self):
1921 """An empty background list cannot be converted to a field, and must
1922 not break the conversion of the images.
1923 """
1924 science, sources, matchedTemplate, difference = self._make_images()
1925 task = self._setup_task(doSubtractBackground=False, image_type="future")
1926 inputImage = self._makeFutureInput(difference)
1927 output = task.run(science, matchedTemplate, inputImage.to_legacy(), sources)
1928 self.assertEqual(len(output.differenceBackground), 0)
1929 legacy = output.subtractedMeasuredExposure
1930 task.convert_outputs_to_future(output, {"subtractedMeasuredExposure": inputImage})
1931 self._check_converted(output.subtractedMeasuredExposure, legacy)
1932 self.assertEqual(len(output.subtractedMeasuredExposure.backgrounds), 0)
1934 def test_stray_struct_field_is_ignored(self):
1935 """``DetectAndMeasureTask`` seeds ``scoreMeasuredExposure`` in its
1936 result struct even though it has no such output connection; check
1937 that such a stray field is ignored.
1938 """
1939 science, sources, matchedTemplate, difference = self._make_images()
1940 task = self._setup_task(image_type="legacy")
1941 output = task.run(science, matchedTemplate, difference, sources)
1942 output.scoreMeasuredExposure = difference
1943 connections = task.config.ConnectionsClass(config=task.config)
1944 self.assertNotIn("scoreMeasuredExposure", connections.outputs)
1945 lsst.pipe.base.testUtils.assertValidOutput(task, output)
1947 def test_run_quantum(self):
1948 """``runQuantum`` converts before putting in future mode, and leaves
1949 the legacy types alone otherwise.
1950 """
1951 templateInfo = [self._makeTemplateInfo()]
1952 kernel = self._makeKernel()
1953 for imageType, expectedType in (("legacy", afwImage.ExposureF),
1954 ("future", DifferenceImage)):
1955 with self.subTest(imageType=imageType):
1956 science, sources, matchedTemplate, difference = self._make_images()
1957 photometricScaling = self._makePhotometricScaling(difference)
1958 if imageType == "future":
1959 difference = self._makeFutureInput(difference, templateInfo, kernel,
1960 photometricScaling)
1961 task = self._setup_task(doSubtractBackground=True, image_type=imageType)
1962 butlerQC = _RecordingQuantumContext(self.dataId)
1963 task.runQuantum(butlerQC,
1964 _FakeRefs(science=science, matchedTemplate=matchedTemplate,
1965 difference=difference, kernelSources=sources),
1966 _FakeRefs())
1967 out = butlerQC.put_values.subtractedMeasuredExposure
1968 self.assertIsInstance(out, expectedType)
1969 if imageType == "legacy":
1970 self.assertIsInstance(butlerQC.put_values.differenceBackground,
1971 afwMath.BackgroundList)
1972 else:
1973 self.assertIsNotNone(out.backgrounds.subtracted)
1974 if imageType == "future":
1975 # The provenance of the image this was measured from
1976 # survives, rather than being worked out again.
1977 self.assertEqual(out.templates, templateInfo)
1978 self.assertEqual(out.kernel, kernel)
1979 self.assertIs(out.photometric_scaling, photometricScaling)
1981 def test_run_quantum_without_provenance(self):
1982 """An input difference image with no kernel or template records
1983 raises an error that names what is missing, and writes nothing.
1984 """
1985 templateInfo = [self._makeTemplateInfo()]
1986 kernel = self._makeKernel()
1987 for missing, templates, inputKernel in (("PSF matching kernel", templateInfo, None),
1988 ("template coadd records", None, kernel)):
1989 with self.subTest(missing=missing):
1990 science, sources, matchedTemplate, difference = self._make_images()
1991 task = self._setup_task(image_type="future")
1992 butlerQC = _RecordingQuantumContext(self.dataId)
1993 refs = _FakeRefs(science=science, matchedTemplate=matchedTemplate,
1994 difference=self._makeFutureInput(difference, templates, inputKernel),
1995 kernelSources=sources)
1997 with self.assertRaisesRegex(lsst.pipe.base.InvalidQuantumError, f"'difference'.*{missing}"):
1998 task.runQuantum(butlerQC, refs, _FakeRefs())
1999 self.assertIsNone(butlerQC.put_values)
2001 def test_run_quantum_score_task_provenance(self):
2002 """Each output takes the kernel of the image it is measured from,
2003 which for the score task are two different kernels.
2004 """
2005 self.detectionTask = detectAndMeasure.DetectAndMeasureScoreTask
2006 science, sources, matchedTemplate, difference, score = self._make_images(withScore=True)
2007 templateInfo = [self._makeTemplateInfo()]
2008 differenceKernel = self._makeKernel(scale=1.0)
2009 scoreKernel = self._makeKernel(scale=3.0)
2010 self.assertNotEqual(differenceKernel, scoreKernel)
2011 photometricScaling = self._makePhotometricScaling(difference)
2012 task = self._setup_task(image_type="future")
2013 butlerQC = _RecordingQuantumContext(self.dataId)
2015 task.runQuantum(butlerQC,
2016 _FakeRefs(science=science, matchedTemplate=matchedTemplate,
2017 difference=self._makeFutureInput(difference, templateInfo,
2018 differenceKernel, photometricScaling),
2019 scoreExposure=self._makeFutureInput(score, templateInfo, scoreKernel,
2020 photometricScaling),
2021 kernelSources=sources),
2022 _FakeRefs())
2024 out = butlerQC.put_values.subtractedMeasuredExposure
2025 scoreOut = butlerQC.put_values.scoreMeasuredExposure
2026 self.assertEqual(out.templates, templateInfo)
2027 self.assertEqual(scoreOut.templates, templateInfo)
2028 self.assertEqual(out.kernel, differenceKernel)
2029 self.assertEqual(scoreOut.kernel, scoreKernel)
2030 self.assertIs(out.photometric_scaling, photometricScaling)
2031 self.assertIs(scoreOut.photometric_scaling, photometricScaling)
2033 def test_connections_future_reads_difference_image(self):
2034 """The inputs whose provenance is copied ask for the type that
2035 carries it, and only in future mode.
2036 """
2037 for imageType, expected in (("legacy", "ExposureF"), ("future", "DifferenceImage")):
2038 with self.subTest(imageType=imageType):
2039 config = detectAndMeasure.DetectAndMeasureScoreTask.ConfigClass()
2040 config.image_type = imageType
2041 connections = config.connections.ConnectionsClass(config=config)
2042 self.assertEqual(connections.difference.storageClass, expected)
2043 self.assertEqual(connections.scoreExposure.storageClass, expected)
2044 # Nothing is copied off these, so the butler converts them.
2045 self.assertEqual(connections.science.storageClass, "ExposureF")
2046 self.assertEqual(connections.matchedTemplate.storageClass, "ExposureF")
2049class _FakeRefs:
2050 """Stand-in for an input or output ``QuantizedConnection``."""
2051 def __init__(self, **kwargs):
2052 self.__dict__.update(kwargs)
2055class _RecordingQuantumContext:
2056 """Minimal `~lsst.pipe.base.QuantumContext` that records what was put.
2058 Parameters
2059 ----------
2060 dataId : `lsst.daf.butler.DataCoordinate`
2061 Data ID to report as the quantum data ID.
2062 """
2063 def __init__(self, dataId):
2064 self.quantum = _FakeRefs(dataId=dataId)
2065 self.put_values = None
2067 def get(self, refs):
2068 return dict(refs.__dict__)
2070 def put(self, values, refs):
2071 self.put_values = values
2074class MockResponse:
2075 """Provide a mock for requests.put calls"""
2076 def __init__(self, json_data, status_code, text):
2077 self.json_data = json_data
2078 self.status_code = status_code
2079 self.text = text
2081 def json(self):
2082 return self.json_data
2084 def raise_for_status(self):
2085 if self.status_code != 200:
2086 raise requests.exceptions.HTTPError
2089def setup_module(module):
2090 lsst.utils.tests.init()
2093class MemoryTestCase(lsst.utils.tests.MemoryTestCase):
2094 pass
2097if __name__ == "__main__": 2097 ↛ 2098line 2097 didn't jump to line 2098 because the condition on line 2097 was never true
2098 lsst.utils.tests.init()
2099 unittest.main()