Coverage for tests/test_detectAndMeasure.py: 99%
1232 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-31 02:11 -0700
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)
929 self.assertTrue('center_xs' in output.glintTrailInfo)
930 self.assertTrue('center_ys' in output.glintTrailInfo)
933class DetectAndMeasureScoreTest(DetectAndMeasureTestBase, lsst.utils.tests.TestCase):
934 detectionTask = detectAndMeasure.DetectAndMeasureScoreTask
936 def test_detection_xy0(self):
937 """Basic functionality test with non-zero x0 and y0.
938 """
939 # Set up the simulated images
940 noiseLevel = 1.
941 staticSeed = 1
942 fluxLevel = 500
943 kernelSize = 31
944 psfSize = 2.4
945 # Buffer source positions away from the score image's EDGE strip.
946 # Preconvolution adds ~kernelSize//2 pixels of EDGE on top of the
947 # EDGE already set by detection smoothing on the science image.
948 templateBorderSize = kernelSize//2
949 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel, "x0": 12345, "y0": 67890,
950 "kernelSize": kernelSize, "templateBorderSize": templateBorderSize}
951 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
952 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
953 difference = science.clone()
954 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
955 scienceKernel = science.psf.getKernel()
956 score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
958 # Configure the detection Task
959 # Set the subtraction residual metric threshold high, since the template
960 # is not actually subtracted from the science image.
961 detectionTask = self._setup_detection(doDeblend=True, badSubtractionRatioThreshold=1.,
962 doSkySources=False)
964 # Run detection and check the results
965 output = detectionTask.run(science, matchedTemplate, difference, score, sources,
966 idFactory=self.idGenerator.make_table_id_factory())
968 # Catalog ids should be very large from this id generator.
969 self.assertTrue(all(output.diaSources['id'] > 1000000000))
970 subtractedMeasuredExposure = output.subtractedMeasuredExposure
972 self.assertImagesEqual(subtractedMeasuredExposure.image, difference.image)
974 self.assertEqual(len(output.diaSources), len(sources))
975 # no sources should be flagged as negative
976 self.assertEqual(len(~output.diaSources["is_negative"]), len(output.diaSources))
977 refIds = []
978 for source in sources:
979 self._check_diaSource(output.diaSources, source, refIds=refIds)
981 def test_measurements_finite(self):
982 """Measured fluxes and centroids should always be finite.
983 """
984 columnNames = ["coord_ra", "coord_dec", "ip_diffim_forced_PsfFlux_instFlux",
985 "ip_diffim_forced_template_PsfFlux_instFlux"]
987 # Set up the simulated images
988 noiseLevel = 1.
989 staticSeed = 1
990 transientSeed = 6
991 xSize = 256
992 ySize = 256
993 psfSize = 2.4
994 kwargs = {"psfSize": psfSize, "x0": 0, "y0": 0,
995 "xSize": xSize, "ySize": ySize}
996 science, sources = makeTestImage(seed=staticSeed, noiseLevel=noiseLevel, noiseSeed=6,
997 nSrc=1, **kwargs)
998 matchedTemplate, _ = makeTestImage(seed=staticSeed, noiseLevel=noiseLevel/4, noiseSeed=7,
999 nSrc=1, **kwargs)
1000 rng = np.random.RandomState(3)
1001 xLoc = np.arange(-5, xSize+5, 10)
1002 rng.shuffle(xLoc)
1003 yLoc = np.arange(-5, ySize+5, 10)
1004 rng.shuffle(yLoc)
1005 transients, transientSources = makeTestImage(seed=transientSeed,
1006 nSrc=len(xLoc), fluxLevel=1000.,
1007 noiseLevel=noiseLevel, noiseSeed=8,
1008 xLoc=xLoc, yLoc=yLoc,
1009 **kwargs)
1010 difference = science.clone()
1011 difference.maskedImage -= matchedTemplate.maskedImage
1012 difference.maskedImage += transients.maskedImage
1013 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
1014 scienceKernel = science.psf.getKernel()
1015 score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
1017 # Configure the detection Task
1018 detectionTask = self._setup_detection(doForcedMeasurement=True)
1020 # Run detection and check the results
1021 output = detectionTask.run(science, matchedTemplate, difference, score, sources)
1023 for column in columnNames:
1024 self._check_values(output.diaSources[column])
1025 self._check_values(output.diaSources.getX(), minValue=0, maxValue=xSize)
1026 self._check_values(output.diaSources.getY(), minValue=0, maxValue=ySize)
1027 self._check_values(output.diaSources.getPsfInstFlux())
1029 def test_detect_transients(self):
1030 """Run detection on a difference image containing transients.
1031 """
1032 # Set up the simulated images
1033 noiseLevel = 1.
1034 staticSeed = 1
1035 transientSeed = 6
1036 nTransients = 10
1037 transientFlux = 1000
1038 fluxLevel = 500
1039 psfSize = 2.4
1040 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel}
1041 scienceBase, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
1042 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
1043 scienceKernel = scienceBase.psf.getKernel()
1044 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
1046 # Configure the detection Task
1047 detectionTask = self._setup_detection(doMerge=False)
1048 kwargs["seed"] = transientSeed
1049 kwargs["nSrc"] = nTransients
1050 kwargs["fluxLevel"] = transientFlux
1052 # Run detection and check the results
1053 def _run_and_check_detections(positive=True):
1054 """Simulate positive or negative transients and run detection.
1056 Parameters
1057 ----------
1058 positive : `bool`, optional
1059 If set, use positive transient sources.
1060 """
1062 transients, transientSources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, **kwargs)
1063 science = scienceBase.clone()
1064 if positive:
1065 science.maskedImage += transients.maskedImage
1066 else:
1067 science.maskedImage -= transients.maskedImage
1068 difference = science.clone()
1069 difference.maskedImage -= matchedTemplate.maskedImage
1070 score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
1071 # NOTE: NoiseReplacer (run by forcedMeasurement) can modify the
1072 # science image if we've e.g. removed parents post-deblending.
1073 # Pass a clone of the science image, so that it doesn't disrupt
1074 # later tests.
1075 output = detectionTask.run(science, matchedTemplate, difference, score, sources)
1076 refIds = []
1077 scale = 1. if positive else -1.
1078 # sources near the edge may have untrustworthy centroids
1079 goodSrcFlags = ~output.diaSources['base_PixelFlags_flag_edge']
1080 for diaSource, goodSrcFlag in zip(output.diaSources, goodSrcFlags):
1081 if goodSrcFlag:
1082 self._check_diaSource(transientSources, diaSource, refIds=refIds, scale=scale)
1083 _run_and_check_detections(positive=True)
1084 _run_and_check_detections(positive=False)
1086 def test_detect_dipoles(self):
1087 """Run detection on a difference image containing dipoles.
1088 """
1089 # Set up the simulated images
1090 noiseLevel = 1.
1091 staticSeed = 1
1092 fluxLevel = 1000
1093 fluxRange = 1.5
1094 nSources = 10
1095 offset = 1
1096 xSize = 300
1097 ySize = 300
1098 kernelSize = 31
1099 psfSize = 2.4
1100 # Avoid placing sources near the edge for this test, so that we can
1101 # easily check that the correct number of sources are detected.
1102 templateBorderSize = kernelSize//2
1103 dipoleFlag = "ip_diffim_DipoleFit_classification"
1104 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel, "fluxRange": fluxRange,
1105 "nSrc": nSources, "templateBorderSize": templateBorderSize, "kernelSize": kernelSize,
1106 "xSize": xSize, "ySize": ySize}
1107 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
1108 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
1109 difference = science.clone()
1110 # Shift the template by a pixel in order to make dipoles in the difference image.
1111 matchedTemplate.image.array[...] = np.roll(matchedTemplate.image.array[...], offset, axis=0)
1112 matchedTemplate.variance.array[...] = np.roll(matchedTemplate.variance.array[...], offset, axis=0)
1113 matchedTemplate.mask.array[...] = np.roll(matchedTemplate.mask.array[...], offset, axis=0)
1114 difference.maskedImage -= matchedTemplate.maskedImage[science.getBBox()]
1115 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
1116 scienceKernel = science.psf.getKernel()
1117 score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
1119 detectionTask = self._setup_detection(badSubtractionRatioThreshold=0.3)
1120 output = detectionTask.run(science, matchedTemplate, difference, score, sources)
1121 diaSources = output.diaSources[~output.diaSources["sky_source"]].copy(deep=True)
1122 self.assertEqual(len(diaSources), len(sources))
1123 refIds = []
1124 for diaSource in diaSources:
1125 if diaSource[dipoleFlag]: 1125 ↛ 1132line 1125 didn't jump to line 1132 because the condition on line 1125 was always true
1126 self._check_diaSource(sources, diaSource, refIds=refIds, scale=0,
1127 rtol=0.05, atol=None, usePsfFlux=False)
1128 self.assertFloatsAlmostEqual(diaSource["ip_diffim_DipoleFit_orientation"],
1129 -np.pi / 2, atol=2.)
1130 self.assertFloatsAlmostEqual(diaSource["ip_diffim_DipoleFit_separation"], offset, rtol=0.1)
1131 else:
1132 raise ValueError("DiaSource with ID %s is not a dipole!", diaSource.getId())
1134 def test_sky_sources(self):
1135 """Add sky sources and check that they are sufficiently far from other
1136 sources and have negligible flux.
1137 """
1138 # Set up the simulated images
1139 noiseLevel = 1.
1140 staticSeed = 1
1141 transientSeed = 6
1142 transientFluxLevel = 1000.
1143 transientFluxRange = 1.5
1144 fluxLevel = 500
1145 psfSize = 2.4
1146 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel}
1147 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
1148 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
1149 transients, transientSources = makeTestImage(seed=transientSeed, psfSize=2.4,
1150 nSrc=10, fluxLevel=transientFluxLevel,
1151 fluxRange=transientFluxRange,
1152 noiseLevel=noiseLevel, noiseSeed=8)
1153 difference = science.clone()
1154 difference.maskedImage -= matchedTemplate.maskedImage
1155 difference.maskedImage += transients.maskedImage
1156 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
1157 scienceKernel = science.psf.getKernel()
1158 kernelWidth = np.max(scienceKernel.getDimensions())//2
1159 score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
1161 # Configure the detection Task
1162 detectionTask = self._setup_detection(doSkySources=True)
1164 # Run detection and check the results
1165 output = detectionTask.run(science, matchedTemplate, difference, score, sources,
1166 idFactory=self.idGenerator.make_table_id_factory())
1167 nSkySourcesGenerated = detectionTask.metadata["n_skySources"]
1168 skySources = output.diaSources[output.diaSources["sky_source"]]
1169 self.assertEqual(len(skySources), nSkySourcesGenerated)
1170 for skySource in skySources:
1171 # Disable the sky_source flag to enable checks.
1172 skySource["sky_source"] = False
1173 # The sky sources should not be close to any other source
1174 with self.assertRaises(AssertionError):
1175 self._check_diaSource(transientSources, skySource, matchDistance=kernelWidth)
1176 with self.assertRaises(AssertionError):
1177 self._check_diaSource(sources, skySource, matchDistance=kernelWidth)
1178 # The sky sources should have low flux levels.
1179 self._check_diaSource(transientSources, skySource, matchDistance=1000, scale=0.,
1180 atol=np.sqrt(transientFluxRange*transientFluxLevel))
1182 # Catalog ids should be very large from this id generator.
1183 self.assertTrue(all(output.diaSources['id'] > 1000000000))
1185 def test_exclude_mask_detections(self):
1186 """Sources with certain bad mask planes set should not be detected.
1187 """
1188 # Set up the simulated images
1189 noiseLevel = 1.
1190 staticSeed = 1
1191 transientSeed = 6
1192 fluxLevel = 500
1193 radius = 2
1194 psfSize = 2.4
1195 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel}
1196 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
1197 science.getInfo().setVisitInfo(makeVisitInfo())
1198 detector = DetectorWrapper(numAmps=1).detector
1199 science.setDetector(detector)
1200 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
1202 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
1203 scienceKernel = science.psf.getKernel()
1204 # Configure the detection Task
1205 badSourceFlags = ["base_PixelFlags_flag_offimage",
1206 "base_PixelFlags_flag_edgeCenterAll",
1207 "base_PixelFlags_flag_badCenterAll",
1208 "base_PixelFlags_flag_saturatedCenterAll", ]
1209 detectionTask = self._setup_detection(nSkySources=0, badSourceFlags=badSourceFlags)
1210 excludeMaskPlanes = ["EDGE", "BAD", "SAT"]
1211 nBad = len(excludeMaskPlanes)
1212 self.assertGreater(nBad, 0)
1213 kwargs["seed"] = transientSeed
1214 kwargs["nSrc"] = nBad
1215 kwargs["fluxLevel"] = 1000
1217 # Run detection and check the results
1218 def _run_and_check_detections(setFlags=True):
1219 transients, transientSources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, **kwargs)
1220 difference = science.clone()
1221 difference.maskedImage -= matchedTemplate.maskedImage
1222 difference.maskedImage += transients.maskedImage
1223 if setFlags:
1224 for src, badMask in zip(transientSources, excludeMaskPlanes):
1225 srcX = int(src.getX())
1226 srcY = int(src.getY())
1227 srcBbox = lsst.geom.Box2I(lsst.geom.Point2I(srcX - radius, srcY - radius),
1228 lsst.geom.Extent2I(2*radius + 1, 2*radius + 1))
1229 difference[srcBbox].mask.array |= lsst.afw.image.Mask.getPlaneBitMask(badMask)
1230 score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
1231 if setFlags:
1232 with self.assertRaises(AlgorithmError):
1233 output = detectionTask.run(science, matchedTemplate, difference, score, sources)
1234 else:
1235 output = detectionTask.run(science, matchedTemplate, difference, score, sources)
1236 refIds = []
1237 goodSrcFlags = checkMask(difference.mask, transientSources, excludeMaskPlanes)
1238 self.assertEqual(np.sum(~goodSrcFlags), 0)
1239 for diaSource, goodSrcFlag in zip(output.diaSources, goodSrcFlags):
1240 if ~goodSrcFlag: 1240 ↛ 1241line 1240 didn't jump to line 1241 because the condition on line 1240 was never true
1241 with self.assertRaises(AssertionError):
1242 self._check_diaSource(transientSources, diaSource, refIds=refIds)
1243 else:
1244 self._check_diaSource(transientSources, diaSource, refIds=refIds)
1245 _run_and_check_detections(setFlags=False)
1246 _run_and_check_detections(setFlags=True)
1248 def test_detect_transients_with_background(self):
1249 """Background measured on the difference image should be subtracted
1250 from the score image so that transients are still detected cleanly.
1251 """
1252 # Set up the simulated images
1253 noiseLevel = 1.
1254 staticSeed = 1
1255 transientSeed = 6
1256 nTransients = 10
1257 transientFlux = 1000
1258 fluxLevel = 500
1259 xSize = 512
1260 ySize = 512
1261 x0 = 123
1262 y0 = 456
1263 kernelSize = 31
1264 psfSize = 2.4
1265 # Avoid placing sources near the edge for this test, so that we can
1266 # easily check that the correct number of sources are detected.
1267 templateBorderSize = kernelSize//2
1268 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel,
1269 "xSize": xSize, "ySize": ySize, "x0": x0, "y0": y0,
1270 "kernelSize": kernelSize, "templateBorderSize": templateBorderSize}
1271 chebyshev_params = [2.2, 2.1, 2.0, 1.2, 1.1, 1.0]
1273 # The background model must cover the grown image bbox, otherwise it
1274 # would be extrapolated outside its declared domain.
1275 bbox2D = lsst.geom.Box2D(
1276 lsst.geom.Point2D(x0 - templateBorderSize, y0 - templateBorderSize),
1277 lsst.geom.Extent2D(xSize + 2*templateBorderSize, ySize + 2*templateBorderSize))
1278 background_model = afwMath.Chebyshev1Function2D(chebyshev_params, bbox2D)
1279 scienceBase, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6,
1280 background=background_model, **kwargs)
1281 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
1282 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
1283 scienceKernel = scienceBase.psf.getKernel()
1285 # Configure the detection Task
1286 detectionTask = self._setup_detection(doMerge=False, doSubtractBackground=True)
1287 kwargs["seed"] = transientSeed
1288 kwargs["nSrc"] = nTransients
1289 kwargs["fluxLevel"] = transientFlux
1291 def _run_and_check_detections(positive=True):
1292 transients, transientSources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, **kwargs)
1293 science = scienceBase.clone()
1294 if positive:
1295 science.maskedImage += transients.maskedImage
1296 else:
1297 science.maskedImage -= transients.maskedImage
1298 difference = science.clone()
1299 difference.maskedImage -= matchedTemplate.maskedImage
1300 score = subtractTask._convolveExposure(difference, scienceKernel,
1301 subtractTask.convolutionControl)
1302 # Record the score image metric before detection so we can verify
1303 # that the background was actually subtracted from it.
1304 originalScoreMedian = np.median(score.image.array[np.isfinite(score.image.array)])
1305 originalDiffimMedian = np.median(difference.image.array[np.isfinite(difference.image.array)])
1306 output = detectionTask.run(science, matchedTemplate, difference, score, sources)
1308 # The background subtraction should leave a non-empty BackgroundList.
1309 self.assertGreater(len(output.differenceBackground), 0)
1310 # The difference image should have been modified in place by the
1311 # background subtraction, so its median should shift noticeably.
1312 subtractedScoreMedian = np.median(score.image.array[np.isfinite(score.image.array)])
1313 self.assertLess(np.abs(subtractedScoreMedian), np.abs(originalScoreMedian))
1314 subtractedDiffimMedian = np.median(difference.image.array[np.isfinite(difference.image.array)])
1315 self.assertLess(np.abs(subtractedDiffimMedian), np.abs(originalDiffimMedian))
1317 refIds = []
1318 scale = 1. if positive else -1.
1319 for transient in transientSources:
1320 self._check_diaSource(output.diaSources, transient, refIds=refIds, scale=scale)
1321 _run_and_check_detections(positive=True)
1322 _run_and_check_detections(positive=False)
1324 def test_background_detection_mask_not_kept(self):
1325 """The detection masks of the difference image should match those of
1326 the score image, without the detections used to fit the background.
1327 """
1328 kwargs = {"seed": 1, "psfSize": 2.4, "fluxLevel": 500, "kernelSize": 31,
1329 "templateBorderSize": 15}
1330 science, sources = makeTestImage(noiseLevel=1., noiseSeed=6, **kwargs)
1331 matchedTemplate, _ = makeTestImage(noiseLevel=0.25, noiseSeed=7, **kwargs)
1332 difference = science.clone()
1333 # Small-scale structure that the first-pass background fit removes
1334 # but the final fit does not, so the two passes detect different
1335 # pixels.
1336 difference.image.array += 0.3*np.sin(np.arange(difference.getWidth())/20.)
1337 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
1338 score = subtractTask._convolveExposure(difference, science.psf.getKernel(),
1339 subtractTask.convolutionControl)
1340 detectionTask = self._setup_detection(doSubtractBackground=True, doSkySources=False,
1341 badSubtractionRatioThreshold=1.)
1343 output = detectionTask.run(science, matchedTemplate, difference, score, sources)
1345 differenceMask = output.subtractedMeasuredExposure.mask
1346 scoreMask = output.scoreMeasuredExposure.mask
1347 for plane in ("DETECTED", "DETECTED_NEGATIVE"):
1348 with self.subTest(plane=plane):
1349 differencePix = differenceMask.array & differenceMask.getPlaneBitMask(plane) > 0
1350 scorePix = scoreMask.array & scoreMask.getPlaneBitMask(plane) > 0
1351 self.assertGreater(scorePix.sum(), 0)
1352 np.testing.assert_array_equal(differencePix, scorePix)
1354 def test_mask_cosmic_rays(self):
1355 """Cosmic rays detected on the difference image should propagate
1356 to the mask of the returned (measured) exposure. This version creates
1357 and uses a Score image for detection.
1358 """
1359 # Set up the simulated images
1360 noiseLevel = 1.
1361 staticSeed = 1
1362 fluxLevel = 500
1363 xSize = 400
1364 ySize = 400
1365 psfSize = 2.4
1366 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel,
1367 "xSize": xSize, "ySize": ySize}
1368 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
1369 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
1370 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
1371 scienceKernel = science.psf.getKernel()
1372 crMask = science.mask.getPlaneBitMask("CR")
1374 # Configure the detection Task
1375 detectionTask = self._setup_detection(doMerge=False, doSkySources=True)
1377 # Test that no CRs are detected when none are present.
1378 transients, _ = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, nSrc=10, **kwargs)
1379 science.maskedImage += transients.maskedImage
1380 difference = science.clone()
1381 difference.maskedImage -= matchedTemplate.maskedImage
1382 score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
1383 output = detectionTask.run(science, matchedTemplate, difference, score, sources)
1384 self.assertFalse(np.any((output.subtractedMeasuredExposure.mask.array & crMask) > 0))
1386 crX0 = round(sources.getX()[0] - science.getX0())
1387 crY0 = round(sources.getY()[0] - science.getY0())
1388 crX1 = round(sources.getX()[1] - science.getX0())
1389 crY1 = round(sources.getY()[1] - science.getY0())
1390 # Inject CR-like shapes into the difference image. CR detection runs
1391 # on the difference, and the mask should propagate to the measured
1392 # exposure returned by the task.
1393 crMaskSetInput = np.zeros(difference.image.array.shape, bool)
1394 crMaskSetInput[crX0:crX0+1, crY0:crY0+5] = True
1395 crMaskSetInput[crX1:crX1+5, crY1:crY1+1] = True
1396 difference.image.array[crMaskSetInput] += 1234
1397 score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
1398 output = detectionTask.run(science, matchedTemplate, difference, score, sources)
1399 crMaskSet = (output.subtractedMeasuredExposure.mask.array & crMask) > 0
1400 self.assertFalse(np.any(crMaskSet[~crMaskSetInput]))
1401 self.assertTrue(np.all(crMaskSet[crMaskSetInput]))
1404class TestNegativePeaks(lsst.utils.tests.TestCase):
1405 """Tests of deblending and merging negative peaks, to test fixes for the
1406 various problems found on DM-48596.
1407 """
1409 def testDeblendNegatives(self):
1410 """Test that negative peaks get deblended and not destroyed: DM-48704.
1411 This is only a test of deblending, not of merging.
1412 """
1413 # Make a science image with one blend of two positive sources, to
1414 # subtract from an empty template, resulting in a negative diffim.
1415 bbox = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(100, 100))
1416 dataset = lsst.meas.base.tests.TestDataset(bbox)
1417 delta = 10
1418 with dataset.addBlend() as family:
1419 family.addChild(instFlux=2E5, centroid=lsst.geom.Point2D(50, 72))
1420 family.addChild(instFlux=2.5E5, centroid=lsst.geom.Point2D(50+delta, 74))
1421 science, catalog = dataset.realize(noise=1.0,
1422 schema=lsst.meas.base.tests.TestDataset.makeMinimalSchema())
1423 dataset = lsst.meas.base.tests.TestDataset(bbox)
1424 template, _ = dataset.realize(noise=1.0,
1425 schema=lsst.meas.base.tests.TestDataset.makeMinimalSchema())
1426 difference = template.clone()
1427 difference.image -= science.image
1429 config = detectAndMeasure.DetectAndMeasureTask.ConfigClass()
1430 config.doDeblend = True
1431 task = detectAndMeasure.DetectAndMeasureTask(config=config)
1432 # prelude steps taken from `detectAndMeasure.run`
1433 task._prepareInputs(difference)
1434 table = lsst.afw.table.SourceTable.make(task.schema)
1435 results = task.detection.run(table=table, exposure=difference, doSmooth=True)
1436 # Just run the deblend step so we can check the footprints independently.
1437 sources, positives, negatives = task._deblend(difference, results.positive, results.negative)
1439 # DM-48704 fixed a problem where the peaks were in the footprints, but
1440 # the spans were empty.
1441 self.assertEqual(len(negatives), 2)
1442 self.assertEqual(len(positives), 0)
1443 self.assertGreater(negatives[0].getFootprint().getSpans().getArea(), 0)
1444 self.assertGreater(negatives[1].getFootprint().getSpans().getArea(), 0)
1445 # Deblended children are HeavyFootprints; we have to make sure the
1446 # pixel values in those are correct (though DetectAndMeasureTask
1447 # doesn't use the fact that they're Heavy).
1448 # The sources are positive in the science image, and negative in the
1449 # diffim, so the minimum value in the deblended negative footprint is
1450 # the maximum value in the science catalog footprint (ignoring the
1451 # noise in the science and template, hence rtol).
1452 # (catalog[0] is the parent; we want the children)
1453 self.assertFloatsAlmostEqual(negatives[0].getFootprint().getImageArray().min(),
1454 -catalog[1].getFootprint().getImageArray().max(), rtol=1e-4)
1455 self.assertFloatsAlmostEqual(negatives[1].getFootprint().getImageArray().min(),
1456 -catalog[2].getFootprint().getImageArray().max(), rtol=1e-4)
1458 self.assertEqual(sources["is_negative"].sum(), 2)
1460 def testMergeFootprints(self):
1461 """Test that merging footprints does not cause negative ones to
1462 disappear (e.g. get merged into non-connected footprints).
1464 As implemented, the diffim will have three positive sources (one a
1465 blended pair), and 7 negative sources (two a blended pair).
1466 """
1467 # Make a template image with multiple blends, designed to trigger the
1468 # negative-footprint-related bug in `FootprintSet.merge`.
1469 bbox = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(400, 200))
1470 dataset = lsst.meas.base.tests.TestDataset(bbox)
1471 # Because these sources are on the template, they will be negative on
1472 # the diffim, unless the pixels are explicitly made negative via
1473 # `template.image.subset` below.
1474 # positive isolated source on diffim
1475 dataset.addSource(instFlux=.5E5, centroid=lsst.geom.Point2D(25, 26))
1476 # negative isolated source on diffim
1477 dataset.addSource(instFlux=.7E5, centroid=lsst.geom.Point2D(75, 24),
1478 shape=lsst.afw.geom.Quadrupole(12, 7, 2))
1479 delta = 10
1480 # negative blended pair on diffim
1481 with dataset.addBlend() as family:
1482 family.addChild(instFlux=1E5, centroid=lsst.geom.Point2D(150, 72))
1483 family.addChild(instFlux=1.5E5, centroid=lsst.geom.Point2D(150+delta, 74))
1484 # positive blended pair on diffim
1485 with dataset.addBlend() as family:
1486 family.addChild(instFlux=2E5, centroid=lsst.geom.Point2D(250, 72))
1487 family.addChild(instFlux=2.5E5, centroid=lsst.geom.Point2D(250+delta, 74))
1488 # negative blended pair on diffim
1489 with dataset.addBlend() as family:
1490 family.addChild(instFlux=3E5, centroid=lsst.geom.Point2D(350, 72))
1491 family.addChild(instFlux=3.5E5, centroid=lsst.geom.Point2D(350+delta, 74))
1492 # negative isolated source on diffim
1493 dataset.addSource(instFlux=4E5, centroid=lsst.geom.Point2D(75, 124))
1494 # negative isolated source on diffim
1495 dataset.addSource(instFlux=5E5, centroid=lsst.geom.Point2D(175, 124))
1496 schema = lsst.meas.base.tests.TestDataset.makeMinimalSchema()
1497 schema.addField("sky_source", "Flag", "Sky source.")
1498 template, catalog = dataset.realize(noise=100.0, schema=schema)
1499 # These will be positive sources on the diffim (as noted above).
1500 template.image.subset(lsst.geom.Box2I(lsst.geom.Point2I(15, 15),
1501 lsst.geom.Point2I(35, 35))).array *= -1
1502 template.image.subset(lsst.geom.Box2I(lsst.geom.Point2I(230, 60),
1503 lsst.geom.Point2I(270, 85))).array *= -1
1504 dataset = lsst.meas.base.tests.TestDataset(bbox)
1505 science, _ = dataset.realize(noise=100.0, schema=schema)
1506 difference = science.clone()
1507 difference.image -= template.image
1509 config = detectAndMeasure.DetectAndMeasureTask.ConfigClass()
1510 config.doDeblend = True
1511 config.raiseOnBadSubtractionRatio = False
1512 config.raiseOnNoDiaSources = False
1513 task = detectAndMeasure.DetectAndMeasureTask(config=config)
1514 result = task.run(science, template, difference, catalog)
1516 # The original bug manifested as the (175,124) source disappaering in
1517 # the merge step, and being included in the peaks of a duplicated
1518 # (75,124) source, even though their footprints are not connected.
1519 mask = np.isclose(result.diaSources["slot_Centroid_x"], 75) & \
1520 np.isclose(result.diaSources["slot_Centroid_y"], 124)
1521 self.assertEqual(mask.sum(), 1)
1522 peaks = result.diaSources[mask][0].getFootprint().peaks
1523 self.assertEqual(len(peaks), 1)
1524 self.assertEqual(peaks[0].getIx(), 75)
1525 self.assertEqual(peaks[0].getIy(), 124)
1527 mask = np.isclose(result.diaSources["slot_Centroid_x"], 175) & \
1528 np.isclose(result.diaSources["slot_Centroid_y"], 124)
1529 self.assertEqual(mask.sum(), 1)
1530 peaks = result.diaSources[mask][0].getFootprint().peaks
1531 self.assertEqual(len(peaks), 1)
1532 self.assertEqual(peaks[0].getIx(), 175)
1533 self.assertEqual(peaks[0].getIy(), 124)
1535 # This checks that all the returned footprints have exactly one peak;
1536 # it's a more general test than the above, but I'm keeping that for
1537 # the more explicit test of the original bugged sources.
1538 peak_x = np.array([peak['f_x'] for src in result.diaSources for peak in src.getFootprint().peaks])
1539 peak_y = np.array([peak['f_y'] for src in result.diaSources for peak in src.getFootprint().peaks])
1540 peaks = np.column_stack((peak_x, peak_y))
1541 unique_peaks, counts = np.unique(peaks, axis=0, return_counts=True)
1542 self.assertEqual(np.sum(counts > 1), 0)
1544 # There are three positive sources that should not have `is_negative`
1545 # set, independent of how deblending/merging of negative sources is
1546 # handled.
1547 self.assertEqual((~result.diaSources["is_negative"]).sum(), 3)
1549 def testReorderPeaksBySignificance(self):
1550 """A merged footprint whose most significant peak is negative should
1551 have that peak first after re-ordering, even though the merge lists
1552 positive peaks before negative ones.
1553 """
1554 config = detectAndMeasure.DetectAndMeasureTask.ConfigClass()
1555 task = detectAndMeasure.DetectAndMeasureTask(config=config)
1557 # Peak schema with a significance field, as SourceDetectionTask
1558 # produces for pixel_stdev thresholds.
1559 mapper = afwTable.SchemaMapper(afwDetection.PeakTable.makeMinimalSchema())
1560 mapper.addMinimalSchema(afwDetection.PeakTable.makeMinimalSchema())
1561 mapper.addOutputField("significance", type=float, doc="")
1562 peakSchema = mapper.getOutputSchema()
1564 diaSources = afwTable.SourceCatalog(task.schema)
1565 source = diaSources.addNew()
1566 spans = afwGeom.SpanSet(geom.Box2I(geom.Point2I(0, 0), geom.Extent2I(10, 10)))
1567 footprint = afwDetection.Footprint(spans, peakSchema)
1568 # Merge order is [positive, negative]: a weak positive peak first,
1569 # then a more significant negative peak.
1570 footprint.addPeak(1, 1, 30.0)
1571 footprint.peaks[-1]["significance"] = 5.0
1572 footprint.addPeak(2, 2, -500.0)
1573 footprint.peaks[-1]["significance"] = 40.0
1574 source.setFootprint(footprint)
1576 task._reorderPeaksBySignificance(diaSources)
1578 peaks = source.getFootprint().peaks
1579 self.assertEqual(peaks[0].getPeakValue(), -500.0)
1580 self.assertEqual(peaks[0]["significance"], 40.0)
1581 self.assertEqual(peaks[1].getPeakValue(), 30.0)
1582 self.assertEqual(peaks[1]["significance"], 5.0)
1585def makeVisitInfo(id=1):
1586 """Return a non-NaN visitInfo."""
1587 return afwImage.VisitInfo(id=id,
1588 exposureTime=10.01,
1589 darkTime=11.02,
1590 date=dafBase.DateTime(65321.1, dafBase.DateTime.MJD, dafBase.DateTime.TAI),
1591 ut1=12345.1,
1592 era=45.1*geom.degrees,
1593 boresightRaDec=geom.SpherePoint(23.1, 73.2, geom.degrees),
1594 boresightAzAlt=geom.SpherePoint(134.5, 33.3, geom.degrees),
1595 boresightAirmass=1.73,
1596 boresightRotAngle=73.2*geom.degrees,
1597 rotType=afwImage.RotType.SKY,
1598 observatory=Observatory(
1599 11.1*geom.degrees, 22.2*geom.degrees, 0.333),
1600 weather=Weather(1.1, 2.2, 34.5),
1601 )
1604class DetectAndMeasureOutputFormatTest(DetectAndMeasureTestBase, lsst.utils.tests.TestCase):
1605 """Tests of the conversion of the output images to `lsst.images` types.
1606 """
1607 detectionTask = detectAndMeasure.DetectAndMeasureTask
1609 @classmethod
1610 def setUpClass(cls):
1611 super().setUpClass()
1612 # `DifferenceImage.from_legacy` needs a detector with per-amplifier
1613 # raw geometry and a FIELD_ANGLE transform; the simpler
1614 # `DetectorWrapper` detector used elsewhere in this file does not
1615 # have either, and fails to convert.
1616 cls.detector = list(CameraWrapper().camera)[0]
1618 def setUp(self):
1619 self.dataId = lsst.daf.butler.DataCoordinate.standardize(
1620 instrument="I",
1621 visit=42,
1622 detector=12,
1623 universe=lsst.daf.butler.DimensionUniverse(),
1624 )
1625 # The exposure record that `runQuantum` rebuilds from the
1626 # observation metadata of the input difference image.
1627 self.exposureRecord = makeTestExposureRecord(self.dataId.universe,
1628 instrument=self.dataId["instrument"],
1629 visit=self.dataId["visit"])
1631 @staticmethod
1632 def _makeTemplateInfo(patch=42):
1633 """Return one record of a coadd that went into a template.
1635 Parameters
1636 ----------
1637 patch : `int`, optional
1638 Identifier of the patch the record describes.
1640 Returns
1641 -------
1642 info : `lsst.images.DifferenceImageTemplateInfo`
1643 The record.
1644 """
1645 return DifferenceImageTemplateInfo(
1646 skymap="skymap", tract=9813, patch=patch, dataset_id=uuid.uuid4(),
1647 dataset_run="a/template/run", bounds=Polygon.from_box(Box.factory[0:50, 0:50]),
1648 psf_shape_xx=4.0, psf_shape_yy=4.0, psf_shape_xy=0.0, psf_shape_flag=False)
1650 @staticmethod
1651 def _makeKernel(scale=1.0):
1652 """Return a PSF matching kernel of the one form that can be attached
1653 to a difference image.
1655 Parameters
1656 ----------
1657 scale : `float`, optional
1658 Coefficient of the first basis kernel, to tell two kernels apart.
1660 Returns
1661 -------
1662 kernel : `lsst.images.convolution_kernels.ImageBasisConvolutionKernel`
1663 The kernel.
1664 """
1665 basis = []
1666 for sigma in (1.0, 2.0):
1667 image = afwImage.ImageD(5, 5)
1668 y, x = np.mgrid[-2:3, -2:3]
1669 image.array[:] = np.exp(-(x**2 + y**2)/(2*sigma**2))
1670 basis.append(afwMath.FixedKernel(image))
1671 bbox = lsst.geom.Box2D(lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Extent2I(100, 100)))
1672 legacy = afwMath.LinearCombinationKernel(basis, afwMath.Chebyshev1Function2D(1, bbox))
1673 legacy.setSpatialParameters([[scale, 0.0, 0.0], [0.5, 0.0, 0.0]])
1674 return ImageBasisConvolutionKernel.from_legacy(legacy)
1676 @staticmethod
1677 def _makePhotometricScaling(exposure):
1678 """Return a constant photometric scaling for an input image.
1680 Parameters
1681 ----------
1682 exposure : `lsst.afw.image.ExposureF`
1683 Image from `_make_images`.
1685 Returns
1686 -------
1687 photometricScaling : `lsst.images.fields.ChebyshevField`
1688 The scaling.
1689 """
1690 return ChebyshevField(Box.from_legacy(exposure.getBBox()), np.array([[0.25]]),
1691 unit=astropy.units.nJy/astropy.units.electron)
1693 def _makeFutureInput(self, exposure, templateInfo=None, kernel=None, photometricScaling=None):
1694 """Convert a legacy image to a `lsst.images.DifferenceImage`.
1696 Parameters
1697 ----------
1698 exposure : `lsst.afw.image.ExposureF`
1699 Image from `_make_images`.
1700 templateInfo : `list` [`lsst.images.DifferenceImageTemplateInfo`],
1701 optional
1702 Record of the coadds that went into the template.
1703 kernel : `lsst.images.convolution_kernels.ImageBasisConvolutionKernel`,
1704 optional
1705 PSF matching kernel.
1706 photometricScaling : `lsst.images.fields.BaseField`, optional
1707 Photometric scaling of the science image.
1709 Returns
1710 -------
1711 image : `lsst.images.DifferenceImage`
1712 The converted image.
1713 """
1714 image = DifferenceImage.from_legacy(
1715 exposure, unit=astropy.units.nJy, exposure_record=self.exposureRecord)
1716 if templateInfo is not None:
1717 image.templates = templateInfo
1718 if kernel is not None:
1719 image.kernel = kernel
1720 if photometricScaling is not None:
1721 image.photometric_scaling = photometricScaling
1722 return image
1724 def _make_images(self, withScore=False):
1725 """Make a science/template/difference image set for detection.
1727 The template is not actually subtracted from the science image, so
1728 that there are plenty of detections; configure the detection task
1729 with ``badSubtractionRatioThreshold=1.`` to match.
1731 Parameters
1732 ----------
1733 withScore : `bool`, optional
1734 Also make and return a score (maximum likelihood) image?
1736 Returns
1737 -------
1738 images : `tuple`
1739 The science exposure, the source catalog, the matched template,
1740 the difference image, and (if ``withScore``) the score image.
1741 """
1742 kernelSize = 31
1743 kwargs = {"seed": 1, "psfSize": 2.4, "fluxLevel": 500, "x0": 12345, "y0": 67890,
1744 "kernelSize": kernelSize, "templateBorderSize": kernelSize//2}
1745 science, sources = makeTestImage(noiseLevel=1., noiseSeed=6, **kwargs)
1746 matchedTemplate, _ = makeTestImage(noiseLevel=0.25, noiseSeed=7, **kwargs)
1747 # Add components to the basic test image
1748 science.setDetector(self.detector)
1749 matchedTemplate.setDetector(self.detector)
1750 for exposure in (science, matchedTemplate):
1751 exposure.info.setVisitInfo(makeTestVisitInfo(self.dataId["visit"]))
1752 difference = science.clone()
1753 if not withScore:
1754 return science, sources, matchedTemplate, difference
1755 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
1756 score = subtractTask._convolveExposure(difference, science.psf.getKernel(),
1757 subtractTask.convolutionControl)
1758 score.setDetector(self.detector)
1759 return science, sources, matchedTemplate, difference, score
1761 def _setup_task(self, **kwargs):
1762 """Configure a detection task to go with `_make_images`.
1764 Parameters
1765 ----------
1766 **kwargs
1767 Any additional config parameters to set.
1769 Returns
1770 -------
1771 task : `lsst.pipe.base.PipelineTask`
1772 The configured task to use for detection and measurement.
1773 """
1774 # The template is not actually subtracted from the science image in
1775 # `_make_images`, so the subtraction residual metric threshold has to
1776 # be set high.
1777 kwargs.setdefault("badSubtractionRatioThreshold", 1.)
1778 return self._setup_detection(doSkySources=False, **kwargs)
1780 def _check_converted(self, image, legacy):
1781 """Check the properties common to every converted output image.
1783 Parameters
1784 ----------
1785 image : `lsst.images.DifferenceImage`
1786 The converted image.
1787 legacy : `lsst.afw.image.ExposureF`
1788 The legacy output exposure it was converted from.
1789 """
1790 self.assertIsInstance(image, DifferenceImage)
1791 self.assertEqual(image.unit, astropy.units.nJy)
1792 self.assertIsNotNone(image.detector)
1793 # `_makeFutureInput` was not given a photometric scaling.
1794 self.assertIsNone(image.photometric_scaling)
1795 compare_masked_image_to_legacy(image, legacy.maskedImage,
1796 plane_map=get_legacy_difference_image_mask_planes())
1797 self.assertTrue(np.any(image.mask.get("DETECTED")))
1799 def test_connections(self):
1800 """Check that the output storage classes follow ``image_type``,
1801 and that the background is a separate output only in legacy mode.
1802 """
1803 for configClass, connectionsClass in (
1804 (detectAndMeasure.DetectAndMeasureConfig, detectAndMeasure.DetectAndMeasureConnections),
1805 (detectAndMeasure.DetectAndMeasureScoreConfig,
1806 detectAndMeasure.DetectAndMeasureScoreConnections),
1807 ):
1808 for imageType, storageClass in (("legacy", "ExposureF"), ("future", "DifferenceImage")):
1809 with self.subTest(configClass=configClass.__name__, imageType=imageType):
1810 config = configClass()
1811 config.image_type = imageType
1812 config.doSubtractBackground = True
1813 with self.assertWarns(FutureWarning):
1814 config.doWriteBackground = True
1815 connections = connectionsClass(config=config)
1816 self.assertEqual(connections.subtractedMeasuredExposure.storageClass, storageClass)
1817 if connectionsClass is detectAndMeasure.DetectAndMeasureScoreConnections:
1818 self.assertEqual(connections.scoreMeasuredExposure.storageClass, storageClass)
1819 if imageType == "legacy":
1820 self.assertIn("differenceBackground", connections.outputs)
1821 self.assertEqual(connections.differenceBackground.storageClass, "Background")
1822 else:
1823 # The background is recorded on the output images.
1824 self.assertNotIn("differenceBackground", connections.outputs)
1826 def test_run_returns_legacy_types(self):
1827 """``run`` always returns legacy types, even in future mode."""
1828 science, sources, matchedTemplate, difference = self._make_images()
1829 task = self._setup_task(image_type="future")
1830 output = task.run(science, matchedTemplate, difference, sources)
1831 self.assertIsInstance(output.subtractedMeasuredExposure, afwImage.ExposureF)
1833 def test_legacy_mode_does_not_convert(self):
1834 """Legacy mode leaves the outputs unconverted."""
1835 science, sources, matchedTemplate, difference = self._make_images()
1836 task = self._setup_task(doSubtractBackground=True, image_type="legacy")
1837 output = task.run(science, matchedTemplate, difference, sources)
1838 self.assertIsInstance(output.subtractedMeasuredExposure, afwImage.ExposureF)
1839 self.assertIsInstance(output.differenceBackground, afwMath.BackgroundList)
1840 lsst.pipe.base.testUtils.assertValidOutput(task, output)
1842 def test_convert_outputs_to_future(self):
1843 """Convert the outputs of ``DetectAndMeasureTask``."""
1844 science, sources, matchedTemplate, difference = self._make_images()
1845 task = self._setup_task(doSubtractBackground=True, image_type="future")
1846 inputImage = self._makeFutureInput(difference)
1847 output = task.run(science, matchedTemplate, inputImage.to_legacy(), sources)
1848 legacy = output.subtractedMeasuredExposure
1849 background = output.differenceBackground
1850 task.convert_outputs_to_future(output, {"subtractedMeasuredExposure": inputImage})
1852 self.assertIs(output.subtractedMeasuredExposure, inputImage)
1853 self._check_converted(inputImage, legacy)
1854 self.assertEqual(list(inputImage.backgrounds), ["subtracted"])
1855 # The separate background dataset is still written in future mode.
1856 self.assertIs(output.differenceBackground, background)
1857 self.assertIsInstance(output.differenceBackground, afwMath.BackgroundList)
1859 def test_convert_outputs_to_future_copied_pixels(self):
1860 """The output pixels are copied into the input image when the legacy
1861 image that was measured does not share them.
1862 """
1863 science, sources, matchedTemplate, difference = self._make_images()
1864 task = self._setup_task(doSubtractBackground=True, image_type="future")
1865 inputImage = self._makeFutureInput(difference)
1866 original = inputImage.image.array.copy()
1867 output = task.run(science, matchedTemplate, inputImage.to_legacy(copy=True), sources)
1868 legacy = output.subtractedMeasuredExposure
1869 # Background subtraction changed the pixels of the copy only.
1870 self.assertFalse(np.array_equal(legacy.image.array, original))
1871 np.testing.assert_array_equal(inputImage.image.array, original)
1873 task.convert_outputs_to_future(output, {"subtractedMeasuredExposure": inputImage})
1875 self._check_converted(inputImage, legacy)
1877 def test_convert_outputs_to_future_source_injection(self):
1878 """Convert a difference image from a pipeline that injects sources.
1880 The injection planes are optional, so the conversion adds them for
1881 an image that has pixels set in them instead of raising.
1882 """
1883 planes = ("INJECTED", "INJECTED_CORE", "INJECTED_TEMPLATE", "INJECTED_CORE_TEMPLATE")
1884 science, sources, matchedTemplate, difference = self._make_images()
1885 task = self._setup_task(image_type="future")
1886 inputImage = self._makeFutureInput(difference)
1887 output = task.run(science, matchedTemplate, inputImage.to_legacy(), sources)
1888 legacy = output.subtractedMeasuredExposure
1889 for n, plane in enumerate(planes):
1890 legacy.mask.addMaskPlane(plane)
1891 legacy.mask.array[0, n] |= legacy.mask.getPlaneBitMask(plane)
1893 task.convert_outputs_to_future(output, {"subtractedMeasuredExposure": inputImage})
1895 image = output.subtractedMeasuredExposure
1896 self.assertIsInstance(image, DifferenceImage)
1897 # Name the optional planes in the map, so that dropping one fails
1898 # the comparison instead of skipping it.
1899 planeMap = get_legacy_difference_image_mask_planes() | get_legacy_optional_mask_planes()
1900 compare_masked_image_to_legacy(image, legacy.maskedImage, plane_map=planeMap)
1902 def test_convert_outputs_to_future_score(self):
1903 """Convert the outputs of ``DetectAndMeasureScoreTask``."""
1904 self.detectionTask = detectAndMeasure.DetectAndMeasureScoreTask
1905 science, sources, matchedTemplate, difference, score = self._make_images(withScore=True)
1906 task = self._setup_task(doSubtractBackground=True, image_type="future")
1907 inputImages = {"subtractedMeasuredExposure": self._makeFutureInput(difference),
1908 "scoreMeasuredExposure": self._makeFutureInput(score)}
1909 output = task.run(science, matchedTemplate,
1910 inputImages["subtractedMeasuredExposure"].to_legacy(),
1911 inputImages["scoreMeasuredExposure"].to_legacy(), sources)
1912 legacyOutputs = {name: getattr(output, name) for name in inputImages}
1913 task.convert_outputs_to_future(output, inputImages)
1915 for name, legacy in legacyOutputs.items():
1916 with self.subTest(name=name):
1917 self.assertIs(getattr(output, name), inputImages[name])
1918 self._check_converted(getattr(output, name), legacy)
1919 # The same background is subtracted from both images.
1920 self.assertIn("subtracted", getattr(output, name).backgrounds)
1922 def test_convert_outputs_to_future_empty_background(self):
1923 """An empty background list cannot be converted to a field, and must
1924 not break the conversion of the images.
1925 """
1926 science, sources, matchedTemplate, difference = self._make_images()
1927 task = self._setup_task(doSubtractBackground=False, image_type="future")
1928 inputImage = self._makeFutureInput(difference)
1929 output = task.run(science, matchedTemplate, inputImage.to_legacy(), sources)
1930 self.assertEqual(len(output.differenceBackground), 0)
1931 legacy = output.subtractedMeasuredExposure
1932 task.convert_outputs_to_future(output, {"subtractedMeasuredExposure": inputImage})
1933 self._check_converted(output.subtractedMeasuredExposure, legacy)
1934 self.assertEqual(len(output.subtractedMeasuredExposure.backgrounds), 0)
1936 def test_stray_struct_field_is_ignored(self):
1937 """``DetectAndMeasureTask`` seeds ``scoreMeasuredExposure`` in its
1938 result struct even though it has no such output connection; check
1939 that such a stray field is ignored.
1940 """
1941 science, sources, matchedTemplate, difference = self._make_images()
1942 task = self._setup_task(image_type="legacy")
1943 output = task.run(science, matchedTemplate, difference, sources)
1944 output.scoreMeasuredExposure = difference
1945 connections = task.config.ConnectionsClass(config=task.config)
1946 self.assertNotIn("scoreMeasuredExposure", connections.outputs)
1947 lsst.pipe.base.testUtils.assertValidOutput(task, output)
1949 def test_run_quantum(self):
1950 """``runQuantum`` converts before putting in future mode, and leaves
1951 the legacy types alone otherwise.
1952 """
1953 templateInfo = [self._makeTemplateInfo()]
1954 kernel = self._makeKernel()
1955 for imageType, expectedType in (("legacy", afwImage.ExposureF),
1956 ("future", DifferenceImage)):
1957 with self.subTest(imageType=imageType):
1958 science, sources, matchedTemplate, difference = self._make_images()
1959 photometricScaling = self._makePhotometricScaling(difference)
1960 if imageType == "future":
1961 difference = self._makeFutureInput(difference, templateInfo, kernel,
1962 photometricScaling)
1963 task = self._setup_task(doSubtractBackground=True, image_type=imageType)
1964 butlerQC = _RecordingQuantumContext(self.dataId)
1965 task.runQuantum(butlerQC,
1966 _FakeRefs(science=science, matchedTemplate=matchedTemplate,
1967 difference=difference, kernelSources=sources),
1968 _FakeRefs())
1969 out = butlerQC.put_values.subtractedMeasuredExposure
1970 self.assertIsInstance(out, expectedType)
1971 if imageType == "legacy":
1972 self.assertIsInstance(butlerQC.put_values.differenceBackground,
1973 afwMath.BackgroundList)
1974 else:
1975 self.assertIsNotNone(out.backgrounds.subtracted)
1976 if imageType == "future":
1977 # The provenance of the image this was measured from
1978 # survives, rather than being worked out again.
1979 self.assertEqual(out.templates, templateInfo)
1980 self.assertEqual(out.kernel, kernel)
1981 self.assertIs(out.photometric_scaling, photometricScaling)
1983 def test_run_quantum_without_provenance(self):
1984 """An input difference image with no kernel or template records
1985 raises an error that names what is missing, and writes nothing.
1986 """
1987 templateInfo = [self._makeTemplateInfo()]
1988 kernel = self._makeKernel()
1989 for missing, templates, inputKernel in (("PSF matching kernel", templateInfo, None),
1990 ("template coadd records", None, kernel)):
1991 with self.subTest(missing=missing):
1992 science, sources, matchedTemplate, difference = self._make_images()
1993 task = self._setup_task(image_type="future")
1994 butlerQC = _RecordingQuantumContext(self.dataId)
1995 refs = _FakeRefs(science=science, matchedTemplate=matchedTemplate,
1996 difference=self._makeFutureInput(difference, templates, inputKernel),
1997 kernelSources=sources)
1999 with self.assertRaisesRegex(lsst.pipe.base.InvalidQuantumError, f"'difference'.*{missing}"):
2000 task.runQuantum(butlerQC, refs, _FakeRefs())
2001 self.assertIsNone(butlerQC.put_values)
2003 def test_run_quantum_score_task_provenance(self):
2004 """Each output takes the kernel of the image it is measured from,
2005 which for the score task are two different kernels.
2006 """
2007 self.detectionTask = detectAndMeasure.DetectAndMeasureScoreTask
2008 science, sources, matchedTemplate, difference, score = self._make_images(withScore=True)
2009 templateInfo = [self._makeTemplateInfo()]
2010 differenceKernel = self._makeKernel(scale=1.0)
2011 scoreKernel = self._makeKernel(scale=3.0)
2012 self.assertNotEqual(differenceKernel, scoreKernel)
2013 photometricScaling = self._makePhotometricScaling(difference)
2014 task = self._setup_task(image_type="future")
2015 butlerQC = _RecordingQuantumContext(self.dataId)
2017 task.runQuantum(butlerQC,
2018 _FakeRefs(science=science, matchedTemplate=matchedTemplate,
2019 difference=self._makeFutureInput(difference, templateInfo,
2020 differenceKernel, photometricScaling),
2021 scoreExposure=self._makeFutureInput(score, templateInfo, scoreKernel,
2022 photometricScaling),
2023 kernelSources=sources),
2024 _FakeRefs())
2026 out = butlerQC.put_values.subtractedMeasuredExposure
2027 scoreOut = butlerQC.put_values.scoreMeasuredExposure
2028 self.assertEqual(out.templates, templateInfo)
2029 self.assertEqual(scoreOut.templates, templateInfo)
2030 self.assertEqual(out.kernel, differenceKernel)
2031 self.assertEqual(scoreOut.kernel, scoreKernel)
2032 self.assertIs(out.photometric_scaling, photometricScaling)
2033 self.assertIs(scoreOut.photometric_scaling, photometricScaling)
2035 def test_connections_future_reads_difference_image(self):
2036 """The inputs whose provenance is copied ask for the type that
2037 carries it, and only in future mode.
2038 """
2039 for imageType, expected in (("legacy", "ExposureF"), ("future", "DifferenceImage")):
2040 with self.subTest(imageType=imageType):
2041 config = detectAndMeasure.DetectAndMeasureScoreTask.ConfigClass()
2042 config.image_type = imageType
2043 connections = config.connections.ConnectionsClass(config=config)
2044 self.assertEqual(connections.difference.storageClass, expected)
2045 self.assertEqual(connections.scoreExposure.storageClass, expected)
2046 # Nothing is copied off these, so the butler converts them.
2047 self.assertEqual(connections.science.storageClass, "ExposureF")
2048 self.assertEqual(connections.matchedTemplate.storageClass, "ExposureF")
2051class _FakeRefs:
2052 """Stand-in for an input or output ``QuantizedConnection``."""
2053 def __init__(self, **kwargs):
2054 self.__dict__.update(kwargs)
2057class _RecordingQuantumContext:
2058 """Minimal `~lsst.pipe.base.QuantumContext` that records what was put.
2060 Parameters
2061 ----------
2062 dataId : `lsst.daf.butler.DataCoordinate`
2063 Data ID to report as the quantum data ID.
2064 """
2065 def __init__(self, dataId):
2066 self.quantum = _FakeRefs(dataId=dataId)
2067 self.put_values = None
2069 def get(self, refs):
2070 return dict(refs.__dict__)
2072 def put(self, values, refs):
2073 self.put_values = values
2076class MockResponse:
2077 """Provide a mock for requests.put calls"""
2078 def __init__(self, json_data, status_code, text):
2079 self.json_data = json_data
2080 self.status_code = status_code
2081 self.text = text
2083 def json(self):
2084 return self.json_data
2086 def raise_for_status(self):
2087 if self.status_code != 200:
2088 raise requests.exceptions.HTTPError
2091def setup_module(module):
2092 lsst.utils.tests.init()
2095class MemoryTestCase(lsst.utils.tests.MemoryTestCase):
2096 pass
2099if __name__ == "__main__": 2099 ↛ 2100line 2099 didn't jump to line 2100 because the condition on line 2099 was never true
2100 lsst.utils.tests.init()
2101 unittest.main()