Coverage for tests / test_detectAndMeasure.py: 9%
801 statements
« prev ^ index » next coverage.py v7.13.4, created at 2026-03-05 23:15 -0800
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
25from unittest import mock
26import requests
28import lsst.afw.geom as afwGeom
29import lsst.afw.image as afwImage
30import lsst.afw.math as afwMath
31import lsst.geom
32from lsst.ip.diffim import detectAndMeasure, subtractImages
33from lsst.afw.table import IdFactory
34from lsst.afw.cameraGeom.testUtils import DetectorWrapper
35import lsst.meas.algorithms as measAlg
36from lsst.pipe.base import InvalidQuantumError, UpstreamFailureNoWorkFound, AlgorithmError
37import lsst.utils.tests
38import lsst.meas.base.tests
39import lsst.daf.base as dafBase
40from lsst.afw.coord import Observatory, Weather
41import lsst.geom as geom
42import lsst.pex.config as pexConfig
44from utils import makeTestImage, checkMask
47class DetectAndMeasureTestBase:
49 def _check_diaSource(self, refSources, diaSource, refIds=None,
50 matchDistance=1., scale=1., usePsfFlux=True,
51 rtol=0.025, atol=None):
52 """Match a diaSource with a source in a reference catalog
53 and compare properties.
55 Parameters
56 ----------
57 refSources : `lsst.afw.table.SourceCatalog`
58 The reference catalog.
59 diaSource : `lsst.afw.table.SourceRecord`
60 The new diaSource to match to the reference catalog.
61 refIds : `list` of `int`, optional
62 Source IDs of previously associated diaSources.
63 matchDistance : `float`, optional
64 Maximum distance allowed between the detected and reference source
65 locations, in pixels.
66 scale : `float`, optional
67 Optional factor to scale the flux by before performing the test.
68 usePsfFlux : `bool`, optional
69 If set, test the PsfInstFlux field, otherwise use ApInstFlux.
70 rtol : `float`, optional
71 Relative tolerance of the flux value test.
72 atol : `float`, optional
73 Absolute tolerance of the flux value test.
74 """
75 if diaSource["sky_source"]:
76 return
77 distance = np.sqrt((diaSource.getX() - refSources.getX())**2
78 + (diaSource.getY() - refSources.getY())**2)
79 self.assertLess(min(distance), matchDistance)
80 src = refSources[np.argmin(distance)]
81 if refIds is not None:
82 # Check that the same source was not previously associated
83 self.assertNotIn(src.getId(), refIds)
84 refIds.append(src.getId())
85 if atol is None:
86 atol = rtol*src.getPsfInstFlux() if usePsfFlux else rtol*src.getApInstFlux()
87 if usePsfFlux:
88 self.assertFloatsAlmostEqual(src.getPsfInstFlux()*scale, diaSource.getPsfInstFlux(),
89 rtol=rtol, atol=atol)
90 else:
91 self.assertFloatsAlmostEqual(src.getApInstFlux()*scale, diaSource.getApInstFlux(),
92 rtol=rtol, atol=atol)
94 def _check_values(self, values, minValue=None, maxValue=None):
95 """Verify that an array has finite values, and optionally that they are
96 within specified minimum and maximum bounds.
98 Parameters
99 ----------
100 values : `numpy.ndarray`
101 Array of values to check.
102 minValue : `float`, optional
103 Minimum allowable value.
104 maxValue : `float`, optional
105 Maximum allowable value.
106 """
107 self.assertTrue(np.all(np.isfinite(values)))
108 if minValue is not None:
109 self.assertTrue(np.all(values >= minValue))
110 if maxValue is not None:
111 self.assertTrue(np.all(values <= maxValue))
113 def _setup_detection(self, doSkySources=True, nSkySources=5,
114 doSubtractBackground=False, run_sattle=False,
115 badSourceFlags=None, **kwargs):
116 """Setup and configure the detection and measurement PipelineTask.
118 Parameters
119 ----------
120 doSkySources : `bool`, optional
121 Generate sky sources.
122 nSkySources : `int`, optional
123 The number of sky sources to add in isolated background regions.
124 **kwargs
125 Any additional config parameters to set.
127 Returns
128 -------
129 `lsst.pipe.base.PipelineTask`
130 The configured Task to use for detection and measurement.
131 """
132 config = self.detectionTask.ConfigClass()
133 config.doSkySources = doSkySources
134 if doSkySources:
135 config.skySources.nSources = nSkySources
136 config.update(**kwargs)
138 config.run_sattle = run_sattle
140 # Make a realistic id generator so that output catalog ids are useful.
141 dataId = lsst.daf.butler.DataCoordinate.standardize(
142 instrument="I",
143 visit=42,
144 detector=12,
145 universe=lsst.daf.butler.DimensionUniverse(),
146 )
147 if badSourceFlags is None:
148 badSourceFlags = ["base_PixelFlags_flag_offimage", ]
149 config.idGenerator.packer.name = "observation"
150 config.idGenerator.packer["observation"].n_observations = 10000
151 config.idGenerator.packer["observation"].n_detectors = 99
152 config.idGenerator.n_releases = 8
153 config.idGenerator.release_id = 2
154 config.doSubtractBackground = doSubtractBackground
155 config.badSourceFlags = badSourceFlags
156 self.idGenerator = config.idGenerator.apply(dataId)
158 return self.detectionTask(config=config)
161class DetectAndMeasureTest(DetectAndMeasureTestBase, lsst.utils.tests.TestCase):
162 detectionTask = detectAndMeasure.DetectAndMeasureTask
164 def test_detection_xy0(self):
165 """Basic functionality test with non-zero x0 and y0.
166 """
167 # Set up the simulated images
168 noiseLevel = 1.
169 staticSeed = 1
170 fluxLevel = 500
171 kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel, "x0": 12345, "y0": 67890}
172 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
173 science.getInfo().setVisitInfo(makeVisitInfo())
174 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
175 difference = science.clone()
177 # Configure the detection Task
178 # Set the subtraction residual metric threshold high, since the template
179 # is not actually subtracted from the science image.
180 detectionTask = self._setup_detection(doDeblend=True, badSubtractionRatioThreshold=1.,
181 doSkySources=False)
183 # Run detection and check the results
184 output = detectionTask.run(science, matchedTemplate, difference, sources,
185 idFactory=self.idGenerator.make_table_id_factory())
186 subtractedMeasuredExposure = output.subtractedMeasuredExposure
188 # Catalog ids should be very large from this id generator.
189 self.assertTrue(all(output.diaSources['id'] > 1000000000))
190 self.assertImagesEqual(subtractedMeasuredExposure.image, difference.image)
192 # all of the sources should have been detected
193 self.assertEqual(len(output.diaSources), len(sources))
194 # no sources should be flagged as negative
195 self.assertEqual(len(~output.diaSources["is_negative"]), len(output.diaSources))
196 refIds = []
197 for source in sources:
198 self._check_diaSource(output.diaSources, source, refIds=refIds)
200 def test_raise_bad_psf(self):
201 """Detection should raise if the PSF width is NaN
202 """
203 # Set up the simulated images
204 noiseLevel = 1.
205 staticSeed = 1
206 fluxLevel = 500
207 kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel, "x0": 12345, "y0": 67890}
208 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
209 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
210 difference = science.clone()
212 psfShape = difference.getPsf().computeBBox(lsst.geom.Point2D(0, 0)).getDimensions()
213 psfcI = afwImage.ImageD(lsst.geom.Extent2I(psfShape[1], psfShape[0]))
214 psfNew = np.zeros(psfShape)
215 psfNew[0, :] = 1
216 psfcI.array = psfNew
217 psfcK = afwMath.FixedKernel(psfcI)
218 psf = measAlg.KernelPsf(psfcK)
219 difference.setPsf(psf)
221 # Configure the detection Task
222 detectionTask = self._setup_detection()
224 # Run detection and check the results
225 with self.assertRaises(UpstreamFailureNoWorkFound):
226 detectionTask.run(science, matchedTemplate, difference, sources)
228 def test_measurements_finite(self):
229 """Measured fluxes and centroids should always be finite.
230 """
231 columnNames = ["coord_ra", "coord_dec", "ip_diffim_forced_PsfFlux_instFlux",
232 "ip_diffim_forced_template_PsfFlux_instFlux"]
234 # Set up the simulated images
235 noiseLevel = 1.
236 staticSeed = 1
237 transientSeed = 6
238 xSize = 256
239 ySize = 256
240 kwargs = {"psfSize": 2.4, "x0": 0, "y0": 0,
241 "xSize": xSize, "ySize": ySize}
242 science, sources = makeTestImage(seed=staticSeed, noiseLevel=noiseLevel, noiseSeed=6,
243 nSrc=1, **kwargs)
244 science.getInfo().setVisitInfo(makeVisitInfo())
245 matchedTemplate, _ = makeTestImage(seed=staticSeed, noiseLevel=noiseLevel/4, noiseSeed=7,
246 nSrc=1, **kwargs)
247 rng = np.random.RandomState(3)
248 xLoc = np.arange(-5, xSize+5, 10)
249 rng.shuffle(xLoc)
250 yLoc = np.arange(-5, ySize+5, 10)
251 rng.shuffle(yLoc)
252 transients, transientSources = makeTestImage(seed=transientSeed,
253 nSrc=len(xLoc), fluxLevel=1000.,
254 noiseLevel=noiseLevel, noiseSeed=8,
255 xLoc=xLoc, yLoc=yLoc,
256 **kwargs)
257 difference = science.clone()
258 difference.maskedImage -= matchedTemplate.maskedImage
259 difference.maskedImage += transients.maskedImage
261 # Configure the detection Task
262 detectionTask = self._setup_detection(doForcedMeasurement=True)
264 # Run detection and check the results
265 output = detectionTask.run(science, matchedTemplate, difference, sources)
267 for column in columnNames:
268 self._check_values(output.diaSources[column])
269 self._check_values(output.diaSources.getX(), minValue=0, maxValue=xSize)
270 self._check_values(output.diaSources.getY(), minValue=0, maxValue=ySize)
271 self._check_values(output.diaSources.getPsfInstFlux())
273 def test_raise_config_schema_mismatch(self):
274 """Check that sources with specified flags are removed from the catalog.
275 """
276 # Configure the detection Task, and and set a config that is not in the schema
277 with self.assertRaises(InvalidQuantumError):
278 self._setup_detection(badSourceFlags=["Bogus_flag_42"])
280 def test_remove_unphysical(self):
281 """Check that sources with specified flags are removed from the catalog.
282 """
283 # Set up the simulated images
284 noiseLevel = 1.
285 staticSeed = 1
286 xSize = 256
287 ySize = 256
288 kwargs = {"psfSize": 2.4, "xSize": xSize, "ySize": ySize}
289 science, sources = makeTestImage(seed=staticSeed, noiseLevel=noiseLevel, noiseSeed=6,
290 nSrc=1, **kwargs)
291 science.getInfo().setVisitInfo(makeVisitInfo())
292 matchedTemplate, _ = makeTestImage(seed=staticSeed, noiseLevel=noiseLevel/4, noiseSeed=7,
293 nSrc=1, **kwargs)
294 transients, transientSources = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=8, nSrc=1, **kwargs)
295 science.maskedImage += transients.maskedImage
296 difference = science.clone()
297 bbox = difference.getBBox()
298 difference.maskedImage -= matchedTemplate.maskedImage
300 # Configure the detection Task, and remove unphysical sources
301 detectionTask = self._setup_detection(doForcedMeasurement=False, doSkySources=True, nSkySources=20,
302 badSourceFlags=["base_PixelFlags_flag_offimage", ])
304 # Run detection and check the results
305 diaSources = detectionTask.run(science, matchedTemplate, difference, sources).diaSources
306 badDiaSrcDoRemove = ~bbox.contains(diaSources.getX(), diaSources.getY())
307 nBadDoRemove = np.count_nonzero(badDiaSrcDoRemove)
308 # Verify that all sources are physical
309 self.assertEqual(nBadDoRemove, 0)
310 # Set a few centroids outside the image bounding box
311 nSetBad = 5
312 for src in diaSources[0: nSetBad]:
313 src["slot_Centroid_x"] += xSize
314 src["slot_Centroid_y"] += ySize
315 src["base_PixelFlags_flag_offimage"] = True
316 # Verify that these sources are outside the image
317 badDiaSrc = ~bbox.contains(diaSources.getX(), diaSources.getY())
318 nBad = np.count_nonzero(badDiaSrc)
319 self.assertEqual(nBad, nSetBad)
320 diaSourcesNoBad = detectionTask._removeBadSources(diaSources)
321 badDiaSrcNoBad = ~bbox.contains(diaSourcesNoBad.getX(), diaSourcesNoBad.getY())
323 # Verify that no sources outside the image bounding box remain
324 self.assertEqual(np.count_nonzero(badDiaSrcNoBad), 0)
325 self.assertEqual(len(diaSourcesNoBad), len(diaSources) - nSetBad)
327 def test_detect_transients(self):
328 """Run detection on a difference image containing transients.
329 """
330 # Set up the simulated images
331 noiseLevel = 1.
332 staticSeed = 1
333 transientSeed = 6
334 fluxLevel = 500
335 kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel}
336 scienceBase, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
337 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
339 # Configure the detection Task
340 detectionTask = self._setup_detection(doMerge=False, doSkySources=True)
341 kwargs["seed"] = transientSeed
342 kwargs["nSrc"] = 10
343 kwargs["fluxLevel"] = 1000
345 # Run detection and check the results
346 def _detection_wrapper(positive=True):
347 transients, transientSources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, **kwargs)
348 science = scienceBase.clone()
349 if positive:
350 science.maskedImage += transients.maskedImage
351 else:
352 science.maskedImage -= transients.maskedImage
353 difference = science.clone()
354 difference.maskedImage -= matchedTemplate.maskedImage
355 # NOTE: NoiseReplacer (run by forcedMeasurement) can modify the
356 # science image if we've e.g. removed parents post-deblending.
357 # Pass a clone of the science image, so that it doesn't disrupt
358 # later tests.
359 output = detectionTask.run(science, matchedTemplate, difference, sources)
360 refIds = []
361 scale = 1. if positive else -1.
362 for diaSource in output.diaSources:
363 self._check_diaSource(transientSources, diaSource, refIds=refIds, scale=scale)
364 _detection_wrapper(positive=True)
365 _detection_wrapper(positive=False)
367 def test_detect_transients_with_background(self):
368 """Run detection on a difference image containing transients and a background.
369 """
370 # Set up the simulated images
371 noiseLevel = 1.
372 staticSeed = 1
373 transientSeed = 6
374 fluxLevel = 500
375 xSize = 512
376 ySize = 512
377 x0 = 123
378 y0 = 456
379 kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel,
380 "xSize": xSize, "ySize": ySize, "x0": x0, "y0": y0}
381 params = [2.2, 2.1, 2.0, 1.2, 1.1, 1.0]
383 bbox2D = lsst.geom.Box2D(lsst.geom.Point2D(x0, y0), lsst.geom.Extent2D(xSize, ySize))
384 background_model = afwMath.Chebyshev1Function2D(params, bbox2D)
385 scienceBase, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6,
386 background=background_model, **kwargs)
387 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
389 # Configure the detection Task
390 detectionTask = self._setup_detection(doMerge=False, doSubtractBackground=True)
391 kwargs["seed"] = transientSeed
392 kwargs["nSrc"] = 10
393 kwargs["fluxLevel"] = 1000
395 # Run detection and check the results
396 def _detection_wrapper(positive=True):
397 transients, transientSources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, **kwargs)
398 science = scienceBase.clone()
399 if positive:
400 science.maskedImage += transients.maskedImage
401 else:
402 science.maskedImage -= transients.maskedImage
403 difference = science.clone()
404 difference.maskedImage -= matchedTemplate.maskedImage
405 # NOTE: NoiseReplacer (run by forcedMeasurement) can modify the
406 # science image if we've e.g. removed parents post-deblending.
407 # Pass a clone of the science image, so that it doesn't disrupt
408 # later tests.
409 output = detectionTask.run(science, matchedTemplate, difference, sources)
410 refIds = []
411 scale = 1. if positive else -1.
412 for transient in transientSources:
413 self._check_diaSource(output.diaSources, transient, refIds=refIds, scale=scale)
414 _detection_wrapper(positive=True)
415 _detection_wrapper(positive=False)
417 def test_missing_mask_planes(self):
418 """Check that detection runs with missing mask planes.
419 """
420 # Set up the simulated images
421 noiseLevel = 1.
422 fluxLevel = 500
423 kwargs = {"psfSize": 2.4, "fluxLevel": fluxLevel, "addMaskPlanes": []}
424 # Use different seeds for the science and template so every source is a diaSource
425 science, sources = makeTestImage(seed=5, noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
426 science.getInfo().setVisitInfo(makeVisitInfo())
427 matchedTemplate, _ = makeTestImage(seed=6, noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
429 difference = science.clone()
430 difference.maskedImage -= matchedTemplate.maskedImage
431 detectionTask = self._setup_detection(raiseOnBadSubtractionRatio=False, raiseOnNoDiaSources=False)
433 # Verify that detection runs without errors
434 detectionTask.run(science, matchedTemplate, difference, sources)
436 def test_detect_dipoles(self):
437 """Run detection on a difference image containing dipoles.
438 """
439 # Set up the simulated images
440 noiseLevel = 1.
441 staticSeed = 1
442 fluxLevel = 1000
443 fluxRange = 1.5
444 nSources = 10
445 offset = 1
446 xSize = 300
447 ySize = 300
448 kernelSize = 32
449 # Avoid placing sources near the edge for this test, so that we can
450 # easily check that the correct number of sources are detected.
451 templateBorderSize = kernelSize//2
452 dipoleFlag = "ip_diffim_DipoleFit_classification"
453 kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel, "fluxRange": fluxRange,
454 "nSrc": nSources, "templateBorderSize": templateBorderSize, "kernelSize": kernelSize,
455 "xSize": xSize, "ySize": ySize}
456 dipoleFlag = "ip_diffim_DipoleFit_classification"
457 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
458 science.getInfo().setVisitInfo(makeVisitInfo())
459 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
460 difference = science.clone()
461 matchedTemplate.image.array[...] = np.roll(matchedTemplate.image.array[...], offset, axis=0)
462 matchedTemplate.variance.array[...] = np.roll(matchedTemplate.variance.array[...], offset, axis=0)
463 matchedTemplate.mask.array[...] = np.roll(matchedTemplate.mask.array[...], offset, axis=0)
464 difference.maskedImage -= matchedTemplate.maskedImage[science.getBBox()]
466 # Need a higher residual metric threshold, since we are purposely
467 # creating poor subtractions.
468 detectionTask = self._setup_detection(doMerge=True, badSubtractionRatioThreshold=0.3)
469 output = detectionTask.run(science, matchedTemplate, difference, sources)
470 diaSources = output.diaSources[~output.diaSources["sky_source"]].copy(deep=True)
471 self.assertEqual(len(diaSources), len(sources))
472 # no sources should be flagged as negative
473 self.assertEqual(len(~diaSources["is_negative"]), len(diaSources))
474 refIds = []
475 for diaSource in diaSources:
476 if diaSource[dipoleFlag]:
477 self._check_diaSource(sources, diaSource, refIds=refIds, scale=0,
478 rtol=0.05, atol=None, usePsfFlux=False)
479 self.assertFloatsAlmostEqual(diaSource["ip_diffim_DipoleFit_orientation"],
480 -np.pi / 2, atol=2.)
481 self.assertFloatsAlmostEqual(diaSource["ip_diffim_DipoleFit_separation"], offset, rtol=0.1)
482 else:
483 raise ValueError("DiaSource with ID %s is not a dipole!", diaSource.getId())
485 def test_sky_sources(self):
486 """Add sky sources and check that they are sufficiently far from other
487 sources and have negligible flux.
488 """
489 # Set up the simulated images
490 noiseLevel = 1.
491 staticSeed = 1
492 transientSeed = 6
493 transientFluxLevel = 1000.
494 transientFluxRange = 1.5
495 fluxLevel = 500
496 kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel}
497 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
498 science.getInfo().setVisitInfo(makeVisitInfo())
499 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
500 transients, transientSources = makeTestImage(seed=transientSeed, psfSize=2.4,
501 nSrc=10, fluxLevel=transientFluxLevel,
502 fluxRange=transientFluxRange,
503 noiseLevel=noiseLevel, noiseSeed=8)
504 difference = science.clone()
505 difference.maskedImage -= matchedTemplate.maskedImage
506 difference.maskedImage += transients.maskedImage
507 kernelWidth = np.max(science.psf.getKernel().getDimensions())//2
509 # Configure the detection Task
510 detectionTask = self._setup_detection(doSkySources=True)
512 # Run detection and check the results
513 output = detectionTask.run(science, matchedTemplate, difference, sources,
514 idFactory=self.idGenerator.make_table_id_factory())
515 skySources = output.diaSources[output.diaSources["sky_source"]]
516 self.assertEqual(len(skySources), detectionTask.config.skySources.nSources)
517 for skySource in skySources:
518 # Disable the sky_source flag to enable checks.
519 skySource["sky_source"] = False
520 # The sky sources should not be close to any other source
521 with self.assertRaises(AssertionError):
522 self._check_diaSource(transientSources, skySource, matchDistance=kernelWidth)
523 with self.assertRaises(AssertionError):
524 self._check_diaSource(sources, skySource, matchDistance=kernelWidth)
525 # The sky sources should have low flux levels.
526 self._check_diaSource(transientSources, skySource, matchDistance=1000, scale=0.,
527 atol=np.sqrt(transientFluxRange*transientFluxLevel))
529 # Catalog ids should be very large from this id generator.
530 self.assertTrue(all(output.diaSources['id'] > 1000000000))
532 def test_exclude_mask_detections(self):
533 """Sources with certain bad mask planes set should not be detected.
534 """
535 # Set up the simulated images
536 noiseLevel = 1.
537 staticSeed = 1
538 transientSeed = 6
539 fluxLevel = 500
540 radius = 2
541 kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel}
542 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
543 science.getInfo().setVisitInfo(makeVisitInfo())
544 detector = DetectorWrapper(numAmps=1).detector
545 science.setDetector(detector)
546 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
548 # Configure the detection Task
549 badSourceFlags = ["base_PixelFlags_flag_offimage",
550 "base_PixelFlags_flag_edgeCenterAll",
551 "base_PixelFlags_flag_badCenterAll",
552 "base_PixelFlags_flag_saturatedCenterAll", ]
553 detectionTask = self._setup_detection(nSkySources=0, badSourceFlags=badSourceFlags)
554 excludeMaskPlanes = ["EDGE", "BAD", "SAT"]
555 nBad = len(excludeMaskPlanes)
556 self.assertGreater(nBad, 0)
557 kwargs["seed"] = transientSeed
558 kwargs["nSrc"] = nBad
559 kwargs["fluxLevel"] = 1000
561 # Run detection and check the results
562 def _detection_wrapper(setFlags=True):
563 transients, transientSources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, **kwargs)
564 difference = science.clone()
565 difference.maskedImage -= matchedTemplate.maskedImage
566 difference.maskedImage += transients.maskedImage
567 if setFlags:
568 for src, badMask in zip(transientSources, excludeMaskPlanes):
569 srcX = int(src.getX())
570 srcY = int(src.getY())
571 srcBbox = lsst.geom.Box2I(lsst.geom.Point2I(srcX - radius, srcY - radius),
572 lsst.geom.Extent2I(2*radius + 1, 2*radius + 1))
573 difference[srcBbox].mask.array |= lsst.afw.image.Mask.getPlaneBitMask(badMask)
574 with self.assertRaises(AlgorithmError):
575 output = detectionTask.run(science, matchedTemplate, difference, sources)
576 else:
577 output = detectionTask.run(science, matchedTemplate, difference, sources)
578 refIds = []
579 goodSrcFlags = checkMask(difference.mask, transientSources, excludeMaskPlanes)
580 self.assertEqual(np.sum(~goodSrcFlags), 0)
581 for diaSource, goodSrcFlag in zip(output.diaSources, goodSrcFlags):
582 if ~goodSrcFlag:
583 with self.assertRaises(AssertionError):
584 self._check_diaSource(transientSources, diaSource, refIds=refIds)
585 else:
586 self._check_diaSource(transientSources, diaSource, refIds=refIds)
587 _detection_wrapper(setFlags=False)
588 _detection_wrapper(setFlags=True)
590 def test_fake_mask_plane_propagation(self):
591 """Test that we have the mask planes related to fakes in diffim images.
592 This is testing method called updateMasks
593 """
594 xSize = 256
595 ySize = 256
596 science, sources = makeTestImage(psfSize=2.4, xSize=xSize, ySize=ySize, doApplyCalibration=True)
597 science_fake_img, science_fake_sources = makeTestImage(
598 psfSize=2.4, xSize=xSize, ySize=ySize, seed=5, nSrc=3, noiseLevel=0.25, fluxRange=1
599 )
600 template, _ = makeTestImage(psfSize=2.4, xSize=xSize, ySize=ySize, doApplyCalibration=True)
601 tmplt_fake_img, tmplt_fake_sources = makeTestImage(
602 psfSize=2.4, xSize=xSize, ySize=ySize, seed=9, nSrc=3, noiseLevel=0.25, fluxRange=1
603 )
604 # created fakes and added them to the images
605 science.image += science_fake_img.image
606 template.image += tmplt_fake_img.image
608 # TODO: DM-40796 update to INJECTED names when source injection gets refactored
609 # adding mask planes to both science and template images
610 science.mask.addMaskPlane("FAKE")
611 science_fake_bitmask = science.mask.getPlaneBitMask("FAKE")
612 template.mask.addMaskPlane("FAKE")
613 template_fake_bitmask = template.mask.getPlaneBitMask("FAKE")
615 # makeTestImage sets the DETECTED plane on the sources; we can use
616 # that to set the FAKE plane on the science and template images.
617 detected = science_fake_img.mask.getPlaneBitMask("DETECTED")
618 fake_pixels = (science_fake_img.mask.array & detected).nonzero()
619 science.mask.array[fake_pixels] |= science_fake_bitmask
620 detected = tmplt_fake_img.mask.getPlaneBitMask("DETECTED")
621 fake_pixels = (tmplt_fake_img.mask.array & detected).nonzero()
622 template.mask.array[fake_pixels] |= science_fake_bitmask
624 science_fake_masked = (science.mask.array & science_fake_bitmask) > 0
625 template_fake_masked = (template.mask.array & template_fake_bitmask) > 0
627 subtractConfig = subtractImages.AlardLuptonSubtractTask.ConfigClass()
628 subtractConfig.sourceSelector.signalToNoise.fluxField = "truth_instFlux"
629 subtractConfig.sourceSelector.signalToNoise.errField = "truth_instFluxErr"
630 subtractTask = subtractImages.AlardLuptonSubtractTask(config=subtractConfig)
631 subtraction = subtractTask.run(template, science, sources)
633 # check subtraction mask plane is set where we set the previous masks
634 diff_mask = subtraction.difference.mask
636 # science mask should be now in INJECTED
637 inj_masked = (diff_mask.array & diff_mask.getPlaneBitMask("INJECTED")) > 0
639 # template mask should be now in INJECTED_TEMPLATE
640 injTmplt_masked = (diff_mask.array & diff_mask.getPlaneBitMask("INJECTED_TEMPLATE")) > 0
642 self.assertFloatsEqual(inj_masked.astype(int), science_fake_masked.astype(int))
643 # The template is convolved, so the INJECTED_TEMPLATE mask plane may
644 # include more pixels than the FAKE mask plane
645 injTmplt_masked &= template_fake_masked
646 self.assertFloatsEqual(injTmplt_masked.astype(int), template_fake_masked.astype(int))
648 # Now check that detection of fakes have the correct flag for injections
649 detectionTask = self._setup_detection()
651 output = detectionTask.run(subtraction.matchedScience,
652 subtraction.matchedTemplate,
653 subtraction.difference,
654 subtraction.kernelSources)
655 diaSources = output.diaSources[~output.diaSources["sky_source"]].copy(deep=True)
657 sci_refIds = []
658 tmpl_refIds = []
659 for diaSrc in diaSources:
660 if diaSrc['base_PsfFlux_instFlux'] > 0:
661 self._check_diaSource(science_fake_sources, diaSrc, scale=1, refIds=sci_refIds)
662 self.assertTrue(diaSrc['base_PixelFlags_flag_injected'])
663 self.assertTrue(diaSrc['base_PixelFlags_flag_injectedCenter'])
664 self.assertFalse(diaSrc['base_PixelFlags_flag_injected_template'])
665 self.assertFalse(diaSrc['base_PixelFlags_flag_injected_templateCenter'])
666 else:
667 self._check_diaSource(tmplt_fake_sources, diaSrc, scale=-1, refIds=tmpl_refIds)
668 self.assertTrue(diaSrc['base_PixelFlags_flag_injected_template'])
669 self.assertTrue(diaSrc['base_PixelFlags_flag_injected_templateCenter'])
670 self.assertFalse(diaSrc['base_PixelFlags_flag_injected'])
671 self.assertFalse(diaSrc['base_PixelFlags_flag_injectedCenter'])
673 def test_mask_streaks(self):
674 """Run detection on a difference image containing a streak.
675 """
676 # Set up the simulated images
677 noiseLevel = 1.
678 staticSeed = 1
679 fluxLevel = 500
680 xSize = 400
681 ySize = 400
682 kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel, "xSize": xSize, "ySize": ySize}
683 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
684 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
686 # Configure the detection Task
687 detectionTask = self._setup_detection(doMerge=False, doMaskStreaks=True,
688 raiseOnNoDiaSources=False, raiseOnBadSubtractionRatio=False)
690 # Test that no streaks are detected
691 difference = science.clone()
692 difference.maskedImage -= matchedTemplate.maskedImage
693 output = detectionTask.run(science, matchedTemplate, difference, sources)
694 outMask = output.subtractedMeasuredExposure.mask.array
695 streakMask = output.subtractedMeasuredExposure.mask.getPlaneBitMask("STREAK")
696 streakMaskSet = (outMask & streakMask) > 0
697 self.assertTrue(np.all(streakMaskSet == 0))
699 # Add streak-like shape and check that streak is detected
700 difference.image.array[20:23, 40:200] += 50
701 output = detectionTask.run(science, matchedTemplate, difference, sources)
702 outMask = output.subtractedMeasuredExposure.mask.array
703 streakMask = output.subtractedMeasuredExposure.mask.getPlaneBitMask("STREAK")
704 streakMaskSet = (outMask & streakMask) > 0
705 self.assertTrue(np.all(streakMaskSet[20:23, 40:200]))
707 # Add a more trapezoid shaped streak across an image that is
708 # fainter and check that it is also detected
709 bbox = science.getBBox()
710 difference = science.clone()
711 difference.maskedImage -= matchedTemplate.maskedImage
712 width = 100
713 amplitude = 5
714 x0 = -100 + bbox.getBeginX()
715 y0 = -100 + bbox.getBeginY()
716 x1 = xSize + x0 + 100
717 y1 = ySize/2
718 corner_coords = [lsst.geom.Point2D(x0, y0),
719 lsst.geom.Point2D(x0, y0 + width),
720 lsst.geom.Point2D(x1, y1),
721 lsst.geom.Point2D(x1 + width, y1)]
722 streak_trapezoid = afwGeom.Polygon(corner_coords)
723 streak_image = streak_trapezoid.createImage(bbox)
724 streak_image *= amplitude
725 difference.image.array += streak_image.array
726 output = detectionTask.run(science, matchedTemplate, difference, sources)
727 outMask = output.subtractedMeasuredExposure.mask.array
728 streakMask = output.subtractedMeasuredExposure.mask.getPlaneBitMask("STREAK")
729 streakMaskSet = (outMask & streakMask) > 0
730 # Check that pixel values in streak_image equal the streak amplitude
731 streak_check = np.where(streak_image.array == amplitude)
732 self.assertTrue(np.all(streakMaskSet[streak_check]))
733 # Check that the entire image was not masked STREAK
734 self.assertFalse(np.all(streakMaskSet))
736 def _setup_sattle_tests(self):
737 noiseLevel = 1.
738 staticSeed = 1
739 fluxLevel = 500
740 shared_kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel,
741 "x0": 12345, "y0": 67890}
742 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6,
743 **shared_kwargs)
744 science.getInfo().setVisitInfo(makeVisitInfo(id=2))
745 detector = DetectorWrapper(numAmps=1).detector
746 science.setDetector(detector)
747 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel / 4,
748 noiseSeed=7, **shared_kwargs)
749 difference = science.clone()
751 detectionTask = self._setup_detection(doDeblend=True,
752 badSubtractionRatioThreshold=1.,
753 doSkySources=False, run_sattle=True)
755 return science, matchedTemplate, difference, sources, detectionTask
757 @mock.patch.dict(os.environ, {"SATTLE_URI_BASE": "fake_host:1234"})
758 def test_sattle_not_available(self):
759 science, matchedTemplate, difference, sources, detectionTask = self._setup_sattle_tests()
761 response = MockResponse({"allow_list": []}, 500, "sattle internal error")
762 with mock.patch('requests.put', return_value=response):
763 with self.assertRaises(requests.exceptions.HTTPError):
764 detectionTask.run(science, matchedTemplate, difference, sources,
765 idFactory=IdFactory.makeSimple())
767 @mock.patch.dict(os.environ, {"SATTLE_URI_BASE": "fake_host:1234"})
768 def test_visit_id_not_in_sattle(self):
769 science, matchedTemplate, difference, sources, detectionTask = self._setup_sattle_tests()
771 response = MockResponse({"allow_list": []}, 404, "missing visit cache")
772 # visit id not in sattle raises
773 with self.assertRaises(requests.exceptions.HTTPError):
774 with mock.patch('lsst.ip.diffim.detectAndMeasure.requests.put',
775 return_value=response):
776 with mock.patch('lsst.ip.diffim.utils.populate_sattle_visit_cache'):
777 detectionTask.run(science, matchedTemplate, difference, sources,
778 idFactory=IdFactory.makeSimple())
780 @mock.patch.dict(os.environ, {"SATTLE_URI_BASE": "fake_host:1234"})
781 def test_filter_satellites_some_allowed(self):
782 science, matchedTemplate, difference, sources, detectionTask = self._setup_sattle_tests()
784 allowed_ids = [1, 5]
785 response = MockResponse({"allow_list": allowed_ids}, 200, "some allowed")
786 with mock.patch('requests.put', return_value=response):
787 output = detectionTask.run(science, matchedTemplate, difference, sources,
788 idFactory=IdFactory.makeSimple())
790 self.assertEqual(len(output.diaSources), 2)
792 # Output should be sources 1 and 5 allowed out of 20
793 self.assertEqual(set(output.diaSources['id']), set(allowed_ids))
795 @mock.patch.dict(os.environ, {"SATTLE_URI_BASE": "fake_host:1234"})
796 def test_filter_satellites_all_allowed(self):
797 science, matchedTemplate, difference, sources, detectionTask = self._setup_sattle_tests()
799 allowed_ids = list(range(1, 21))
800 response = MockResponse({"allow_list": allowed_ids}, 200, "all allowed")
801 # Run detection and check the results
802 with mock.patch('requests.put', return_value=response):
803 output = detectionTask.run(science, matchedTemplate, difference, sources,
804 idFactory=IdFactory.makeSimple())
806 # Output should be all sources that went in. 20 go in, 20 should come out
807 self.assertEqual(len(output.diaSources), 20)
809 self.assertEqual(set(output.diaSources['id']), set(allowed_ids))
811 @mock.patch.dict(os.environ, {"SATTLE_URI_BASE": "fake_host:1234"})
812 def test_filter_satellites_none_allowed(self):
813 science, matchedTemplate, difference, sources, detectionTask = self._setup_sattle_tests()
815 response = MockResponse({"allow_list": []}, 200, "none allowed")
816 # Run detection and confirm it raises for no diasources
817 with self.assertRaises(detectAndMeasure.NoDiaSourcesError):
818 with mock.patch('requests.put', return_value=response):
819 detectionTask.run(science, matchedTemplate, difference, sources,
820 idFactory=IdFactory.makeSimple())
822 @mock.patch.dict(os.environ, {"SATTLE_URI_BASE": ""})
823 def test_fail_on_sattle_misconfiguration(self):
824 with self.assertRaises(pexConfig.FieldValidationError):
825 self._setup_detection(run_sattle=True)
827 def test_trailed_glints(self):
828 """Test that the glint_trail column works, and that
829 the trailed_glints output contains the expected information.
830 """
831 noiseLevel = 1.
832 staticSeed = 1
833 diffim, diaSources = makeTestImage(seed=staticSeed, noiseLevel=noiseLevel, noiseSeed=6)
834 self._check_values(diaSources['glint_trail'])
836 # Run detection and return the output Struct so we can check it
837 def _detection_wrapper(diffim, diaSources):
838 detectionTask = self._setup_detection()
839 scienceBase, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6)
840 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7)
841 science = scienceBase.clone()
842 science.maskedImage -= diffim.maskedImage
843 difference = science.clone()
844 difference.maskedImage -= matchedTemplate.maskedImage
845 output = detectionTask.run(science, matchedTemplate, difference, sources)
846 return output
848 output = _detection_wrapper(diffim, diaSources)
849 self.assertTrue('slopes' in output.glintTrailInfo)
850 self.assertTrue('intercepts' in output.glintTrailInfo)
851 self.assertTrue('stderrs' in output.glintTrailInfo)
852 self.assertTrue('lengths' in output.glintTrailInfo)
853 self.assertTrue('angles' in output.glintTrailInfo)
856class DetectAndMeasureScoreTest(DetectAndMeasureTestBase, lsst.utils.tests.TestCase):
857 detectionTask = detectAndMeasure.DetectAndMeasureScoreTask
859 def test_detection_xy0(self):
860 """Basic functionality test with non-zero x0 and y0.
861 """
862 # Set up the simulated images
863 noiseLevel = 1.
864 staticSeed = 1
865 fluxLevel = 500
866 kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel, "x0": 12345, "y0": 67890}
867 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
868 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
869 difference = science.clone()
870 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
871 scienceKernel = science.psf.getKernel()
872 score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
874 # Configure the detection Task
875 # Set the subtraction residual metric threshold high, since the template
876 # is not actually subtracted from the science image.
877 detectionTask = self._setup_detection(doDeblend=True, badSubtractionRatioThreshold=1.,
878 doSkySources=False)
880 # Run detection and check the results
881 output = detectionTask.run(science, matchedTemplate, difference, score, sources,
882 idFactory=self.idGenerator.make_table_id_factory())
884 # Catalog ids should be very large from this id generator.
885 self.assertTrue(all(output.diaSources['id'] > 1000000000))
886 subtractedMeasuredExposure = output.subtractedMeasuredExposure
888 self.assertImagesEqual(subtractedMeasuredExposure.image, difference.image)
890 self.assertEqual(len(output.diaSources), len(sources))
891 # no sources should be flagged as negative
892 self.assertEqual(len(~output.diaSources["is_negative"]), len(output.diaSources))
893 # TODO DM-41496: restore this block once we handle detections on edge
894 # pixels better; at least one of these sources currently has a bad
895 # centroid because most of the source is rejected as EDGE.
896 # refIds = []
897 # for source in sources:
898 # self._check_diaSource(output.diaSources, source, refIds=refIds)
900 def test_measurements_finite(self):
901 """Measured fluxes and centroids should always be finite.
902 """
903 columnNames = ["coord_ra", "coord_dec", "ip_diffim_forced_PsfFlux_instFlux",
904 "ip_diffim_forced_template_PsfFlux_instFlux"]
906 # Set up the simulated images
907 noiseLevel = 1.
908 staticSeed = 1
909 transientSeed = 6
910 xSize = 256
911 ySize = 256
912 kwargs = {"psfSize": 2.4, "x0": 0, "y0": 0,
913 "xSize": xSize, "ySize": ySize}
914 science, sources = makeTestImage(seed=staticSeed, noiseLevel=noiseLevel, noiseSeed=6,
915 nSrc=1, **kwargs)
916 matchedTemplate, _ = makeTestImage(seed=staticSeed, noiseLevel=noiseLevel/4, noiseSeed=7,
917 nSrc=1, **kwargs)
918 rng = np.random.RandomState(3)
919 xLoc = np.arange(-5, xSize+5, 10)
920 rng.shuffle(xLoc)
921 yLoc = np.arange(-5, ySize+5, 10)
922 rng.shuffle(yLoc)
923 transients, transientSources = makeTestImage(seed=transientSeed,
924 nSrc=len(xLoc), fluxLevel=1000.,
925 noiseLevel=noiseLevel, noiseSeed=8,
926 xLoc=xLoc, yLoc=yLoc,
927 **kwargs)
928 difference = science.clone()
929 difference.maskedImage -= matchedTemplate.maskedImage
930 difference.maskedImage += transients.maskedImage
931 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
932 scienceKernel = science.psf.getKernel()
933 score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
935 # Configure the detection Task
936 detectionTask = self._setup_detection(doForcedMeasurement=True)
938 # Run detection and check the results
939 output = detectionTask.run(science, matchedTemplate, difference, score, sources)
941 for column in columnNames:
942 self._check_values(output.diaSources[column])
943 self._check_values(output.diaSources.getX(), minValue=0, maxValue=xSize)
944 self._check_values(output.diaSources.getY(), minValue=0, maxValue=ySize)
945 self._check_values(output.diaSources.getPsfInstFlux())
947 def test_detect_transients(self):
948 """Run detection on a difference image containing transients.
949 """
950 # Set up the simulated images
951 noiseLevel = 1.
952 staticSeed = 1
953 transientSeed = 6
954 fluxLevel = 500
955 kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel}
956 scienceBase, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
957 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
958 scienceKernel = scienceBase.psf.getKernel()
959 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
961 # Configure the detection Task
962 detectionTask = self._setup_detection(doMerge=False)
963 kwargs["seed"] = transientSeed
964 kwargs["nSrc"] = 10
965 kwargs["fluxLevel"] = 1000
967 # Run detection and check the results
968 def _detection_wrapper(positive=True):
969 """Simulate positive or negative transients and run detection.
971 Parameters
972 ----------
973 positive : `bool`, optional
974 If set, use positive transient sources.
975 """
977 transients, transientSources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, **kwargs)
978 science = scienceBase.clone()
979 if positive:
980 science.maskedImage += transients.maskedImage
981 else:
982 science.maskedImage -= transients.maskedImage
983 difference = science.clone()
984 difference.maskedImage -= matchedTemplate.maskedImage
985 score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
986 # NOTE: NoiseReplacer (run by forcedMeasurement) can modify the
987 # science image if we've e.g. removed parents post-deblending.
988 # Pass a clone of the science image, so that it doesn't disrupt
989 # later tests.
990 output = detectionTask.run(science, matchedTemplate, difference, score, sources)
991 refIds = []
992 scale = 1. if positive else -1.
993 # sources near the edge may have untrustworthy centroids
994 goodSrcFlags = ~output.diaSources['base_PixelFlags_flag_edge']
995 for diaSource, goodSrcFlag in zip(output.diaSources, goodSrcFlags):
996 if goodSrcFlag:
997 self._check_diaSource(transientSources, diaSource, refIds=refIds, scale=scale)
998 _detection_wrapper(positive=True)
999 _detection_wrapper(positive=False)
1001 def test_detect_dipoles(self):
1002 """Run detection on a difference image containing dipoles.
1003 """
1004 # Set up the simulated images
1005 noiseLevel = 1.
1006 staticSeed = 1
1007 fluxLevel = 1000
1008 fluxRange = 1.5
1009 nSources = 10
1010 offset = 1
1011 xSize = 300
1012 ySize = 300
1013 kernelSize = 32
1014 # Avoid placing sources near the edge for this test, so that we can
1015 # easily check that the correct number of sources are detected.
1016 templateBorderSize = kernelSize//2
1017 dipoleFlag = "ip_diffim_DipoleFit_classification"
1018 kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel, "fluxRange": fluxRange,
1019 "nSrc": nSources, "templateBorderSize": templateBorderSize, "kernelSize": kernelSize,
1020 "xSize": xSize, "ySize": ySize}
1021 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
1022 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
1023 difference = science.clone()
1024 # Shift the template by a pixel in order to make dipoles in the difference image.
1025 matchedTemplate.image.array[...] = np.roll(matchedTemplate.image.array[...], offset, axis=0)
1026 matchedTemplate.variance.array[...] = np.roll(matchedTemplate.variance.array[...], offset, axis=0)
1027 matchedTemplate.mask.array[...] = np.roll(matchedTemplate.mask.array[...], offset, axis=0)
1028 difference.maskedImage -= matchedTemplate.maskedImage[science.getBBox()]
1029 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
1030 scienceKernel = science.psf.getKernel()
1031 score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
1033 detectionTask = self._setup_detection(badSubtractionRatioThreshold=0.3)
1034 output = detectionTask.run(science, matchedTemplate, difference, score, sources)
1035 diaSources = output.diaSources[~output.diaSources["sky_source"]].copy(deep=True)
1036 self.assertEqual(len(diaSources), len(sources))
1037 refIds = []
1038 for diaSource in diaSources:
1039 if diaSource[dipoleFlag]:
1040 self._check_diaSource(sources, diaSource, refIds=refIds, scale=0,
1041 rtol=0.05, atol=None, usePsfFlux=False)
1042 self.assertFloatsAlmostEqual(diaSource["ip_diffim_DipoleFit_orientation"],
1043 -np.pi / 2, atol=2.)
1044 self.assertFloatsAlmostEqual(diaSource["ip_diffim_DipoleFit_separation"], offset, rtol=0.1)
1045 else:
1046 raise ValueError("DiaSource with ID %s is not a dipole!", diaSource.getId())
1048 def test_sky_sources(self):
1049 """Add sky sources and check that they are sufficiently far from other
1050 sources and have negligible flux.
1051 """
1052 # Set up the simulated images
1053 noiseLevel = 1.
1054 staticSeed = 1
1055 transientSeed = 6
1056 transientFluxLevel = 1000.
1057 transientFluxRange = 1.5
1058 fluxLevel = 500
1059 kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel}
1060 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
1061 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
1062 transients, transientSources = makeTestImage(seed=transientSeed, psfSize=2.4,
1063 nSrc=10, fluxLevel=transientFluxLevel,
1064 fluxRange=transientFluxRange,
1065 noiseLevel=noiseLevel, noiseSeed=8)
1066 difference = science.clone()
1067 difference.maskedImage -= matchedTemplate.maskedImage
1068 difference.maskedImage += transients.maskedImage
1069 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
1070 scienceKernel = science.psf.getKernel()
1071 kernelWidth = np.max(scienceKernel.getDimensions())//2
1072 score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
1074 # Configure the detection Task
1075 detectionTask = self._setup_detection(doSkySources=True)
1077 # Run detection and check the results
1078 output = detectionTask.run(science, matchedTemplate, difference, score, sources,
1079 idFactory=self.idGenerator.make_table_id_factory())
1080 nSkySourcesGenerated = detectionTask.metadata["n_skySources"]
1081 skySources = output.diaSources[output.diaSources["sky_source"]]
1082 self.assertEqual(len(skySources), nSkySourcesGenerated)
1083 for skySource in skySources:
1084 # Disable the sky_source flag to enable checks.
1085 skySource["sky_source"] = False
1086 # The sky sources should not be close to any other source
1087 with self.assertRaises(AssertionError):
1088 self._check_diaSource(transientSources, skySource, matchDistance=kernelWidth)
1089 with self.assertRaises(AssertionError):
1090 self._check_diaSource(sources, skySource, matchDistance=kernelWidth)
1091 # The sky sources should have low flux levels.
1092 self._check_diaSource(transientSources, skySource, matchDistance=1000, scale=0.,
1093 atol=np.sqrt(transientFluxRange*transientFluxLevel))
1095 # Catalog ids should be very large from this id generator.
1096 self.assertTrue(all(output.diaSources['id'] > 1000000000))
1098 def test_exclude_mask_detections(self):
1099 """Sources with certain bad mask planes set should not be detected.
1100 """
1101 # Set up the simulated images
1102 noiseLevel = 1.
1103 staticSeed = 1
1104 transientSeed = 6
1105 fluxLevel = 500
1106 radius = 2
1107 kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel}
1108 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
1109 science.getInfo().setVisitInfo(makeVisitInfo())
1110 detector = DetectorWrapper(numAmps=1).detector
1111 science.setDetector(detector)
1112 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
1114 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
1115 scienceKernel = science.psf.getKernel()
1116 # Configure the detection Task
1117 badSourceFlags = ["base_PixelFlags_flag_offimage",
1118 "base_PixelFlags_flag_edgeCenterAll",
1119 "base_PixelFlags_flag_badCenterAll",
1120 "base_PixelFlags_flag_saturatedCenterAll", ]
1121 detectionTask = self._setup_detection(nSkySources=0, badSourceFlags=badSourceFlags)
1122 excludeMaskPlanes = ["EDGE", "BAD", "SAT"]
1123 nBad = len(excludeMaskPlanes)
1124 self.assertGreater(nBad, 0)
1125 kwargs["seed"] = transientSeed
1126 kwargs["nSrc"] = nBad
1127 kwargs["fluxLevel"] = 1000
1129 # Run detection and check the results
1130 def _detection_wrapper(setFlags=True):
1131 transients, transientSources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, **kwargs)
1132 difference = science.clone()
1133 difference.maskedImage -= matchedTemplate.maskedImage
1134 difference.maskedImage += transients.maskedImage
1135 if setFlags:
1136 for src, badMask in zip(transientSources, excludeMaskPlanes):
1137 srcX = int(src.getX())
1138 srcY = int(src.getY())
1139 srcBbox = lsst.geom.Box2I(lsst.geom.Point2I(srcX - radius, srcY - radius),
1140 lsst.geom.Extent2I(2*radius + 1, 2*radius + 1))
1141 difference[srcBbox].mask.array |= lsst.afw.image.Mask.getPlaneBitMask(badMask)
1142 score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
1143 if setFlags:
1144 with self.assertRaises(AlgorithmError):
1145 output = detectionTask.run(science, matchedTemplate, difference, score, sources)
1146 else:
1147 output = detectionTask.run(science, matchedTemplate, difference, score, sources)
1148 refIds = []
1149 goodSrcFlags = checkMask(difference.mask, transientSources, excludeMaskPlanes)
1150 self.assertEqual(np.sum(~goodSrcFlags), 0)
1151 for diaSource, goodSrcFlag in zip(output.diaSources, goodSrcFlags):
1152 if ~goodSrcFlag:
1153 with self.assertRaises(AssertionError):
1154 self._check_diaSource(transientSources, diaSource, refIds=refIds)
1155 else:
1156 self._check_diaSource(transientSources, diaSource, refIds=refIds)
1157 _detection_wrapper(setFlags=False)
1158 _detection_wrapper(setFlags=True)
1161class TestNegativePeaks(lsst.utils.tests.TestCase):
1162 """Tests of deblending and merging negative peaks, to test fixes for the
1163 various problems found on DM-48596.
1164 """
1166 def testDeblendNegatives(self):
1167 """Test that negative peaks get deblended and not destroyed: DM-48704.
1168 This is only a test of deblending, not of merging.
1169 """
1170 # Make a science image with one blend of two positive sources, to
1171 # subtract from an empty template, resulting in a negative diffim.
1172 bbox = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(100, 100))
1173 dataset = lsst.meas.base.tests.TestDataset(bbox)
1174 delta = 10
1175 with dataset.addBlend() as family:
1176 family.addChild(instFlux=2E5, centroid=lsst.geom.Point2D(50, 72))
1177 family.addChild(instFlux=2.5E5, centroid=lsst.geom.Point2D(50+delta, 74))
1178 science, catalog = dataset.realize(noise=1.0,
1179 schema=lsst.meas.base.tests.TestDataset.makeMinimalSchema())
1180 dataset = lsst.meas.base.tests.TestDataset(bbox)
1181 template, _ = dataset.realize(noise=1.0,
1182 schema=lsst.meas.base.tests.TestDataset.makeMinimalSchema())
1183 difference = template.clone()
1184 difference.image -= science.image
1186 config = detectAndMeasure.DetectAndMeasureTask.ConfigClass()
1187 config.doDeblend = True
1188 task = detectAndMeasure.DetectAndMeasureTask(config=config)
1189 # prelude steps taken from `detectAndMeasure.run`
1190 task._prepareInputs(difference)
1191 table = lsst.afw.table.SourceTable.make(task.schema)
1192 results = task.detection.run(table=table, exposure=difference, doSmooth=True)
1193 # Just run the deblend step so we can check the footprints independently.
1194 sources, positives, negatives = task._deblend(difference, results.positive, results.negative)
1196 # DM-48704 fixed a problem where the peaks were in the footprints, but
1197 # the spans were empty.
1198 self.assertEqual(len(negatives), 2)
1199 self.assertEqual(len(positives), 0)
1200 self.assertGreater(negatives[0].getFootprint().getSpans().getArea(), 0)
1201 self.assertGreater(negatives[1].getFootprint().getSpans().getArea(), 0)
1202 # Deblended children are HeavyFootprints; we have to make sure the
1203 # pixel values in those are correct (though DetectAndMeasureTask
1204 # doesn't use the fact that they're Heavy).
1205 # The sources are positive in the science image, and negative in the
1206 # diffim, so the minimum value in the deblended negative footprint is
1207 # the maximum value in the science catalog footprint (ignoring the
1208 # noise in the science and template, hence rtol).
1209 # (catalog[0] is the parent; we want the children)
1210 self.assertFloatsAlmostEqual(negatives[0].getFootprint().getImageArray().min(),
1211 -catalog[1].getFootprint().getImageArray().max(), rtol=1e-4)
1212 self.assertFloatsAlmostEqual(negatives[1].getFootprint().getImageArray().min(),
1213 -catalog[2].getFootprint().getImageArray().max(), rtol=1e-4)
1215 self.assertEqual(sources["is_negative"].sum(), 2)
1217 def testMergeFootprints(self):
1218 """Test that merging footprints does not cause negative ones to
1219 disappear (e.g. get merged into non-connected footprints).
1221 As implemented, the diffim will have three positive sources (one a
1222 blended pair), and 7 negative sources (two a blended pair).
1223 """
1224 # Make a template image with multiple blends, designed to trigger the
1225 # negative-footprint-related bug in `FootprintSet.merge`.
1226 bbox = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(400, 200))
1227 dataset = lsst.meas.base.tests.TestDataset(bbox)
1228 # Because these sources are on the template, they will be negative on
1229 # the diffim, unless the pixels are explicitly made negative via
1230 # `template.image.subset` below.
1231 # positive isolated source on diffim
1232 dataset.addSource(instFlux=.5E5, centroid=lsst.geom.Point2D(25, 26))
1233 # negative isolated source on diffim
1234 dataset.addSource(instFlux=.7E5, centroid=lsst.geom.Point2D(75, 24),
1235 shape=lsst.afw.geom.Quadrupole(12, 7, 2))
1236 delta = 10
1237 # negative blended pair on diffim
1238 with dataset.addBlend() as family:
1239 family.addChild(instFlux=1E5, centroid=lsst.geom.Point2D(150, 72))
1240 family.addChild(instFlux=1.5E5, centroid=lsst.geom.Point2D(150+delta, 74))
1241 # positive blended pair on diffim
1242 with dataset.addBlend() as family:
1243 family.addChild(instFlux=2E5, centroid=lsst.geom.Point2D(250, 72))
1244 family.addChild(instFlux=2.5E5, centroid=lsst.geom.Point2D(250+delta, 74))
1245 # negative blended pair on diffim
1246 with dataset.addBlend() as family:
1247 family.addChild(instFlux=3E5, centroid=lsst.geom.Point2D(350, 72))
1248 family.addChild(instFlux=3.5E5, centroid=lsst.geom.Point2D(350+delta, 74))
1249 # negative isolated source on diffim
1250 dataset.addSource(instFlux=4E5, centroid=lsst.geom.Point2D(75, 124))
1251 # negative isolated source on diffim
1252 dataset.addSource(instFlux=5E5, centroid=lsst.geom.Point2D(175, 124))
1253 schema = lsst.meas.base.tests.TestDataset.makeMinimalSchema()
1254 schema.addField("sky_source", "Flag", "Sky source.")
1255 template, catalog = dataset.realize(noise=100.0, schema=schema)
1256 # These will be positive sources on the diffim (as noted above).
1257 template.image.subset(lsst.geom.Box2I(lsst.geom.Point2I(15, 15),
1258 lsst.geom.Point2I(35, 35))).array *= -1
1259 template.image.subset(lsst.geom.Box2I(lsst.geom.Point2I(230, 60),
1260 lsst.geom.Point2I(270, 85))).array *= -1
1261 dataset = lsst.meas.base.tests.TestDataset(bbox)
1262 science, _ = dataset.realize(noise=100.0, schema=schema)
1263 difference = science.clone()
1264 difference.image -= template.image
1266 config = detectAndMeasure.DetectAndMeasureTask.ConfigClass()
1267 config.doDeblend = True
1268 config.raiseOnBadSubtractionRatio = False
1269 config.raiseOnNoDiaSources = False
1270 task = detectAndMeasure.DetectAndMeasureTask(config=config)
1271 result = task.run(science, template, difference, catalog)
1273 # The original bug manifested as the (175,124) source disappaering in
1274 # the merge step, and being included in the peaks of a duplicated
1275 # (75,124) source, even though their footprints are not connected.
1276 mask = np.isclose(result.diaSources["slot_Centroid_x"], 75) & \
1277 np.isclose(result.diaSources["slot_Centroid_y"], 124)
1278 self.assertEqual(mask.sum(), 1)
1279 peaks = result.diaSources[mask][0].getFootprint().peaks
1280 self.assertEqual(len(peaks), 1)
1281 self.assertEqual(peaks[0].getIx(), 75)
1282 self.assertEqual(peaks[0].getIy(), 124)
1284 mask = np.isclose(result.diaSources["slot_Centroid_x"], 175) & \
1285 np.isclose(result.diaSources["slot_Centroid_y"], 124)
1286 self.assertEqual(mask.sum(), 1)
1287 peaks = result.diaSources[mask][0].getFootprint().peaks
1288 self.assertEqual(len(peaks), 1)
1289 self.assertEqual(peaks[0].getIx(), 175)
1290 self.assertEqual(peaks[0].getIy(), 124)
1292 # This checks that all the returned footprints have exactly one peak;
1293 # it's a more general test than the above, but I'm keeping that for
1294 # the more explicit test of the original bugged sources.
1295 peak_x = np.array([peak['f_x'] for src in result.diaSources for peak in src.getFootprint().peaks])
1296 peak_y = np.array([peak['f_y'] for src in result.diaSources for peak in src.getFootprint().peaks])
1297 peaks = np.column_stack((peak_x, peak_y))
1298 unique_peaks, counts = np.unique(peaks, axis=0, return_counts=True)
1299 self.assertEqual(np.sum(counts > 1), 0)
1301 # There are three positive sources that should not have `is_negative`
1302 # set, independent of how deblending/merging of negative sources is
1303 # handled.
1304 self.assertEqual((~result.diaSources["is_negative"]).sum(), 3)
1307def makeVisitInfo(id=1):
1308 """Return a non-NaN visitInfo."""
1309 return afwImage.VisitInfo(id=id,
1310 exposureTime=10.01,
1311 darkTime=11.02,
1312 date=dafBase.DateTime(65321.1, dafBase.DateTime.MJD, dafBase.DateTime.TAI),
1313 ut1=12345.1,
1314 era=45.1*geom.degrees,
1315 boresightRaDec=geom.SpherePoint(23.1, 73.2, geom.degrees),
1316 boresightAzAlt=geom.SpherePoint(134.5, 33.3, geom.degrees),
1317 boresightAirmass=1.73,
1318 boresightRotAngle=73.2*geom.degrees,
1319 rotType=afwImage.RotType.SKY,
1320 observatory=Observatory(
1321 11.1*geom.degrees, 22.2*geom.degrees, 0.333),
1322 weather=Weather(1.1, 2.2, 34.5),
1323 )
1326class MockResponse:
1327 """Provide a mock for requests.put calls"""
1328 def __init__(self, json_data, status_code, text):
1329 self.json_data = json_data
1330 self.status_code = status_code
1331 self.text = text
1333 def json(self):
1334 return self.json_data
1336 def raise_for_status(self):
1337 if self.status_code != 200:
1338 raise requests.exceptions.HTTPError
1341def setup_module(module):
1342 lsst.utils.tests.init()
1345class MemoryTestCase(lsst.utils.tests.MemoryTestCase):
1346 pass
1349if __name__ == "__main__": 1349 ↛ 1350line 1349 didn't jump to line 1350 because the condition on line 1349 was never true
1350 lsst.utils.tests.init()
1351 unittest.main()