Coverage for tests / test_detectAndMeasure.py: 8%
829 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-05-06 08:49 +0000
1# This file is part of ip_diffim.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
22import numpy as np
23import os
24import unittest
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_mask_cosmic_rays(self):
418 """Run detection on a difference image containing a cosmic ray.
419 """
420 # Set up the simulated images
421 noiseLevel = 1.
422 staticSeed = 1
423 fluxLevel = 500
424 xSize = 400
425 ySize = 400
426 kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel, "xSize": xSize, "ySize": ySize}
427 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
428 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
429 crMask = science.mask.getPlaneBitMask("CR")
431 # Configure the detection Task
432 detectionTask = self._setup_detection(doMerge=False, doSkySources=True)
434 # Test that no CRs are detected
435 transients, transientSources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, nSrc=10, **kwargs)
436 science.maskedImage += transients.maskedImage
437 difference = science.clone()
438 difference.maskedImage -= matchedTemplate.maskedImage
439 output = detectionTask.run(science, matchedTemplate, difference, sources)
440 crMaskSet = (output.subtractedMeasuredExposure.mask.array & crMask) > 0
441 self.assertTrue(np.all(crMaskSet == 0))
443 crX0 = round(sources.getX()[0] - science.getX0())
444 crY0 = round(sources.getY()[0] - science.getY0())
445 crX1 = round(sources.getX()[1] - science.getX0())
446 crY1 = round(sources.getY()[1] - science.getY0())
447 # Add CR-like shape and check that CR is detected
448 # Pick two locations on top of sources, since that is what is likely to
449 # be missed in the first stage of CR rejection.
450 difference.image.array[crX0:crX0+1, crY0:crY0+5] += 1234 # Arbitrary but large flux for the CRs
451 difference.image.array[crX1:crX1+5, crY1:crY1+1] += 1234
452 output = detectionTask.run(science, matchedTemplate, difference)
453 crMaskSet = (output.subtractedMeasuredExposure.mask.array & crMask) > 0
454 self.assertTrue(np.all(crMaskSet[crX0:crX0+1, crY0:crY0+5]))
455 self.assertTrue(np.all(crMaskSet[crX1:crX1+5, crY1:crY1+1]))
457 def test_missing_mask_planes(self):
458 """Check that detection runs with missing mask planes.
459 """
460 # Set up the simulated images
461 noiseLevel = 1.
462 fluxLevel = 500
463 kwargs = {"psfSize": 2.4, "fluxLevel": fluxLevel, "addMaskPlanes": []}
464 # Use different seeds for the science and template so every source is a diaSource
465 science, sources = makeTestImage(seed=5, noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
466 science.getInfo().setVisitInfo(makeVisitInfo())
467 matchedTemplate, _ = makeTestImage(seed=6, noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
469 difference = science.clone()
470 difference.maskedImage -= matchedTemplate.maskedImage
471 detectionTask = self._setup_detection(raiseOnBadSubtractionRatio=False, raiseOnNoDiaSources=False)
473 # Verify that detection runs without errors
474 detectionTask.run(science, matchedTemplate, difference, sources)
476 def test_detect_dipoles(self):
477 """Run detection on a difference image containing dipoles.
478 """
479 # Set up the simulated images
480 noiseLevel = 1.
481 staticSeed = 1
482 fluxLevel = 1000
483 fluxRange = 1.5
484 nSources = 10
485 offset = 1
486 xSize = 300
487 ySize = 300
488 kernelSize = 32
489 # Avoid placing sources near the edge for this test, so that we can
490 # easily check that the correct number of sources are detected.
491 templateBorderSize = kernelSize//2
492 dipoleFlag = "ip_diffim_DipoleFit_classification"
493 kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel, "fluxRange": fluxRange,
494 "nSrc": nSources, "templateBorderSize": templateBorderSize, "kernelSize": kernelSize,
495 "xSize": xSize, "ySize": ySize}
496 dipoleFlag = "ip_diffim_DipoleFit_classification"
497 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
498 science.getInfo().setVisitInfo(makeVisitInfo())
499 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
500 difference = science.clone()
501 matchedTemplate.image.array[...] = np.roll(matchedTemplate.image.array[...], offset, axis=0)
502 matchedTemplate.variance.array[...] = np.roll(matchedTemplate.variance.array[...], offset, axis=0)
503 matchedTemplate.mask.array[...] = np.roll(matchedTemplate.mask.array[...], offset, axis=0)
504 difference.maskedImage -= matchedTemplate.maskedImage[science.getBBox()]
506 # Need a higher residual metric threshold, since we are purposely
507 # creating poor subtractions.
508 detectionTask = self._setup_detection(doMerge=True, badSubtractionRatioThreshold=0.3)
509 output = detectionTask.run(science, matchedTemplate, difference, sources)
510 diaSources = output.diaSources[~output.diaSources["sky_source"]].copy(deep=True)
511 self.assertEqual(len(diaSources), len(sources))
512 # no sources should be flagged as negative
513 self.assertEqual(len(~diaSources["is_negative"]), len(diaSources))
514 refIds = []
515 for diaSource in diaSources:
516 if diaSource[dipoleFlag]:
517 self._check_diaSource(sources, diaSource, refIds=refIds, scale=0,
518 rtol=0.05, atol=None, usePsfFlux=False)
519 self.assertFloatsAlmostEqual(diaSource["ip_diffim_DipoleFit_orientation"],
520 -np.pi / 2, atol=2.)
521 self.assertFloatsAlmostEqual(diaSource["ip_diffim_DipoleFit_separation"], offset, rtol=0.1)
522 else:
523 raise ValueError("DiaSource with ID %s is not a dipole!", diaSource.getId())
525 def test_sky_sources(self):
526 """Add sky sources and check that they are sufficiently far from other
527 sources and have negligible flux.
528 """
529 # Set up the simulated images
530 noiseLevel = 1.
531 staticSeed = 1
532 transientSeed = 6
533 transientFluxLevel = 1000.
534 transientFluxRange = 1.5
535 fluxLevel = 500
536 kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel}
537 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
538 science.getInfo().setVisitInfo(makeVisitInfo())
539 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
540 transients, transientSources = makeTestImage(seed=transientSeed, psfSize=2.4,
541 nSrc=10, fluxLevel=transientFluxLevel,
542 fluxRange=transientFluxRange,
543 noiseLevel=noiseLevel, noiseSeed=8)
544 difference = science.clone()
545 difference.maskedImage -= matchedTemplate.maskedImage
546 difference.maskedImage += transients.maskedImage
547 kernelWidth = np.max(science.psf.getKernel().getDimensions())//2
549 # Configure the detection Task
550 detectionTask = self._setup_detection(doSkySources=True)
552 # Run detection and check the results
553 output = detectionTask.run(science, matchedTemplate, difference, sources,
554 idFactory=self.idGenerator.make_table_id_factory())
555 skySources = output.diaSources[output.diaSources["sky_source"]]
556 self.assertEqual(len(skySources), detectionTask.config.skySources.nSources)
557 for skySource in skySources:
558 # Disable the sky_source flag to enable checks.
559 skySource["sky_source"] = False
560 # The sky sources should not be close to any other source
561 with self.assertRaises(AssertionError):
562 self._check_diaSource(transientSources, skySource, matchDistance=kernelWidth)
563 with self.assertRaises(AssertionError):
564 self._check_diaSource(sources, skySource, matchDistance=kernelWidth)
565 # The sky sources should have low flux levels.
566 self._check_diaSource(transientSources, skySource, matchDistance=1000, scale=0.,
567 atol=np.sqrt(transientFluxRange*transientFluxLevel))
569 # Catalog ids should be very large from this id generator.
570 self.assertTrue(all(output.diaSources['id'] > 1000000000))
572 def test_exclude_mask_detections(self):
573 """Sources with certain bad mask planes set should not be detected.
574 """
575 # Set up the simulated images
576 noiseLevel = 1.
577 staticSeed = 1
578 transientSeed = 6
579 fluxLevel = 500
580 radius = 2
581 kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel}
582 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
583 science.getInfo().setVisitInfo(makeVisitInfo())
584 detector = DetectorWrapper(numAmps=1).detector
585 science.setDetector(detector)
586 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
588 # Configure the detection Task
589 badSourceFlags = ["base_PixelFlags_flag_offimage",
590 "base_PixelFlags_flag_edgeCenterAll",
591 "base_PixelFlags_flag_badCenterAll",
592 "base_PixelFlags_flag_saturatedCenterAll", ]
593 detectionTask = self._setup_detection(nSkySources=0, badSourceFlags=badSourceFlags)
594 excludeMaskPlanes = ["EDGE", "BAD", "SAT"]
595 nBad = len(excludeMaskPlanes)
596 self.assertGreater(nBad, 0)
597 kwargs["seed"] = transientSeed
598 kwargs["nSrc"] = nBad
599 kwargs["fluxLevel"] = 1000
601 # Run detection and check the results
602 def _detection_wrapper(setFlags=True):
603 transients, transientSources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, **kwargs)
604 difference = science.clone()
605 difference.maskedImage -= matchedTemplate.maskedImage
606 difference.maskedImage += transients.maskedImage
607 if setFlags:
608 for src, badMask in zip(transientSources, excludeMaskPlanes):
609 srcX = int(src.getX())
610 srcY = int(src.getY())
611 srcBbox = lsst.geom.Box2I(lsst.geom.Point2I(srcX - radius, srcY - radius),
612 lsst.geom.Extent2I(2*radius + 1, 2*radius + 1))
613 difference[srcBbox].mask.array |= lsst.afw.image.Mask.getPlaneBitMask(badMask)
614 with self.assertRaises(AlgorithmError):
615 output = detectionTask.run(science, matchedTemplate, difference, sources)
616 else:
617 output = detectionTask.run(science, matchedTemplate, difference, sources)
618 refIds = []
619 goodSrcFlags = checkMask(difference.mask, transientSources, excludeMaskPlanes)
620 self.assertEqual(np.sum(~goodSrcFlags), 0)
621 for diaSource, goodSrcFlag in zip(output.diaSources, goodSrcFlags):
622 if ~goodSrcFlag:
623 with self.assertRaises(AssertionError):
624 self._check_diaSource(transientSources, diaSource, refIds=refIds)
625 else:
626 self._check_diaSource(transientSources, diaSource, refIds=refIds)
627 _detection_wrapper(setFlags=False)
628 _detection_wrapper(setFlags=True)
630 def test_fake_mask_plane_propagation(self):
631 """Test that we have the mask planes related to fakes in diffim images.
632 This is testing method called updateMasks
633 """
634 xSize = 256
635 ySize = 256
636 science, sources = makeTestImage(psfSize=2.4, xSize=xSize, ySize=ySize, doApplyCalibration=True)
637 science_fake_img, science_fake_sources = makeTestImage(
638 psfSize=2.4, xSize=xSize, ySize=ySize, seed=5, nSrc=3, noiseLevel=0.25, fluxRange=1
639 )
640 template, _ = makeTestImage(psfSize=2.4, xSize=xSize, ySize=ySize, doApplyCalibration=True)
641 tmplt_fake_img, tmplt_fake_sources = makeTestImage(
642 psfSize=2.4, xSize=xSize, ySize=ySize, seed=9, nSrc=3, noiseLevel=0.25, fluxRange=1
643 )
644 # created fakes and added them to the images
645 science.image += science_fake_img.image
646 template.image += tmplt_fake_img.image
648 # TODO: DM-40796 update to INJECTED names when source injection gets refactored
649 # adding mask planes to both science and template images
650 science.mask.addMaskPlane("FAKE")
651 science_fake_bitmask = science.mask.getPlaneBitMask("FAKE")
652 template.mask.addMaskPlane("FAKE")
653 template_fake_bitmask = template.mask.getPlaneBitMask("FAKE")
655 # makeTestImage sets the DETECTED plane on the sources; we can use
656 # that to set the FAKE plane on the science and template images.
657 detected = science_fake_img.mask.getPlaneBitMask("DETECTED")
658 fake_pixels = (science_fake_img.mask.array & detected).nonzero()
659 science.mask.array[fake_pixels] |= science_fake_bitmask
660 detected = tmplt_fake_img.mask.getPlaneBitMask("DETECTED")
661 fake_pixels = (tmplt_fake_img.mask.array & detected).nonzero()
662 template.mask.array[fake_pixels] |= science_fake_bitmask
664 science_fake_masked = (science.mask.array & science_fake_bitmask) > 0
665 template_fake_masked = (template.mask.array & template_fake_bitmask) > 0
667 subtractConfig = subtractImages.AlardLuptonSubtractTask.ConfigClass()
668 subtractConfig.sourceSelector.signalToNoise.fluxField = "truth_instFlux"
669 subtractConfig.sourceSelector.signalToNoise.errField = "truth_instFluxErr"
670 subtractTask = subtractImages.AlardLuptonSubtractTask(config=subtractConfig)
671 subtraction = subtractTask.run(template, science, sources)
673 # check subtraction mask plane is set where we set the previous masks
674 diff_mask = subtraction.difference.mask
676 # science mask should be now in INJECTED
677 inj_masked = (diff_mask.array & diff_mask.getPlaneBitMask("INJECTED")) > 0
679 # template mask should be now in INJECTED_TEMPLATE
680 injTmplt_masked = (diff_mask.array & diff_mask.getPlaneBitMask("INJECTED_TEMPLATE")) > 0
682 self.assertFloatsEqual(inj_masked.astype(int), science_fake_masked.astype(int))
683 # The template is convolved, so the INJECTED_TEMPLATE mask plane may
684 # include more pixels than the FAKE mask plane
685 injTmplt_masked &= template_fake_masked
686 self.assertFloatsEqual(injTmplt_masked.astype(int), template_fake_masked.astype(int))
688 # Now check that detection of fakes have the correct flag for injections
689 detectionTask = self._setup_detection()
691 output = detectionTask.run(subtraction.matchedScience,
692 subtraction.matchedTemplate,
693 subtraction.difference,
694 subtraction.kernelSources)
695 diaSources = output.diaSources[~output.diaSources["sky_source"]].copy(deep=True)
697 sci_refIds = []
698 tmpl_refIds = []
699 for diaSrc in diaSources:
700 if diaSrc['base_PsfFlux_instFlux'] > 0:
701 self._check_diaSource(science_fake_sources, diaSrc, scale=1, refIds=sci_refIds)
702 self.assertTrue(diaSrc['base_PixelFlags_flag_injected'])
703 self.assertTrue(diaSrc['base_PixelFlags_flag_injectedCenter'])
704 self.assertFalse(diaSrc['base_PixelFlags_flag_injected_template'])
705 self.assertFalse(diaSrc['base_PixelFlags_flag_injected_templateCenter'])
706 else:
707 self._check_diaSource(tmplt_fake_sources, diaSrc, scale=-1, refIds=tmpl_refIds)
708 self.assertTrue(diaSrc['base_PixelFlags_flag_injected_template'])
709 self.assertTrue(diaSrc['base_PixelFlags_flag_injected_templateCenter'])
710 self.assertFalse(diaSrc['base_PixelFlags_flag_injected'])
711 self.assertFalse(diaSrc['base_PixelFlags_flag_injectedCenter'])
713 def test_mask_streaks(self):
714 """Run detection on a difference image containing a streak.
715 """
716 # Set up the simulated images
717 noiseLevel = 1.
718 staticSeed = 1
719 fluxLevel = 500
720 xSize = 400
721 ySize = 400
722 kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel, "xSize": xSize, "ySize": ySize}
723 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
724 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
726 # Configure the detection Task
727 detectionTask = self._setup_detection(doMerge=False, doMaskStreaks=True,
728 raiseOnNoDiaSources=False, raiseOnBadSubtractionRatio=False)
730 # Test that no streaks are detected
731 difference = science.clone()
732 difference.maskedImage -= matchedTemplate.maskedImage
733 output = detectionTask.run(science, matchedTemplate, difference, sources)
734 outMask = output.subtractedMeasuredExposure.mask.array
735 streakMask = output.subtractedMeasuredExposure.mask.getPlaneBitMask("STREAK")
736 streakMaskSet = (outMask & streakMask) > 0
737 self.assertTrue(np.all(streakMaskSet == 0))
739 # Add streak-like shape and check that streak is detected
740 difference.image.array[20:23, 40:200] += 50
741 output = detectionTask.run(science, matchedTemplate, difference, sources)
742 outMask = output.subtractedMeasuredExposure.mask.array
743 streakMask = output.subtractedMeasuredExposure.mask.getPlaneBitMask("STREAK")
744 streakMaskSet = (outMask & streakMask) > 0
745 self.assertTrue(np.all(streakMaskSet[20:23, 40:200]))
747 # Add a more trapezoid shaped streak across an image that is
748 # fainter and check that it is also detected
749 bbox = science.getBBox()
750 difference = science.clone()
751 difference.maskedImage -= matchedTemplate.maskedImage
752 width = 100
753 amplitude = 5
754 x0 = -100 + bbox.getBeginX()
755 y0 = -100 + bbox.getBeginY()
756 x1 = xSize + x0 + 100
757 y1 = ySize/2
758 corner_coords = [lsst.geom.Point2D(x0, y0),
759 lsst.geom.Point2D(x0, y0 + width),
760 lsst.geom.Point2D(x1, y1),
761 lsst.geom.Point2D(x1 + width, y1)]
762 streak_trapezoid = afwGeom.Polygon(corner_coords)
763 streak_image = streak_trapezoid.createImage(bbox)
764 streak_image *= amplitude
765 difference.image.array += streak_image.array
766 output = detectionTask.run(science, matchedTemplate, difference, sources)
767 outMask = output.subtractedMeasuredExposure.mask.array
768 streakMask = output.subtractedMeasuredExposure.mask.getPlaneBitMask("STREAK")
769 streakMaskSet = (outMask & streakMask) > 0
770 # Check that pixel values in streak_image equal the streak amplitude
771 streak_check = np.where(streak_image.array == amplitude)
772 self.assertTrue(np.all(streakMaskSet[streak_check]))
773 # Check that the entire image was not masked STREAK
774 self.assertFalse(np.all(streakMaskSet))
776 def _setup_sattle_tests(self):
777 noiseLevel = 1.
778 staticSeed = 1
779 fluxLevel = 500
780 shared_kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel,
781 "x0": 12345, "y0": 67890}
782 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6,
783 **shared_kwargs)
784 science.getInfo().setVisitInfo(makeVisitInfo(id=2))
785 detector = DetectorWrapper(numAmps=1).detector
786 science.setDetector(detector)
787 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel / 4,
788 noiseSeed=7, **shared_kwargs)
789 difference = science.clone()
791 detectionTask = self._setup_detection(doDeblend=True,
792 badSubtractionRatioThreshold=1.,
793 doSkySources=False, run_sattle=True)
795 return science, matchedTemplate, difference, sources, detectionTask
797 @mock.patch.dict(os.environ, {"SATTLE_URI_BASE": "fake_host:1234"})
798 def test_sattle_not_available(self):
799 science, matchedTemplate, difference, sources, detectionTask = self._setup_sattle_tests()
801 response = MockResponse({"allow_list": []}, 500, "sattle internal error")
802 with mock.patch('requests.put', return_value=response):
803 with self.assertRaises(requests.exceptions.HTTPError):
804 detectionTask.run(science, matchedTemplate, difference, sources,
805 idFactory=IdFactory.makeSimple())
807 @mock.patch.dict(os.environ, {"SATTLE_URI_BASE": "fake_host:1234"})
808 def test_visit_id_not_in_sattle(self):
809 science, matchedTemplate, difference, sources, detectionTask = self._setup_sattle_tests()
811 response = MockResponse({"allow_list": []}, 404, "missing visit cache")
812 # visit id not in sattle raises
813 with self.assertRaises(requests.exceptions.HTTPError):
814 with mock.patch('lsst.ip.diffim.detectAndMeasure.requests.put',
815 return_value=response):
816 with mock.patch('lsst.ip.diffim.utils.populate_sattle_visit_cache'):
817 detectionTask.run(science, matchedTemplate, difference, sources,
818 idFactory=IdFactory.makeSimple())
820 @mock.patch.dict(os.environ, {"SATTLE_URI_BASE": "fake_host:1234"})
821 def test_filter_satellites_some_allowed(self):
822 science, matchedTemplate, difference, sources, detectionTask = self._setup_sattle_tests()
824 allowed_ids = [1, 5]
825 response = MockResponse({"allow_list": allowed_ids}, 200, "some allowed")
826 with mock.patch('requests.put', return_value=response):
827 output = detectionTask.run(science, matchedTemplate, difference, sources,
828 idFactory=IdFactory.makeSimple())
830 self.assertEqual(len(output.diaSources), 2)
832 # Output should be sources 1 and 5 allowed out of 20
833 self.assertEqual(set(output.diaSources['id']), set(allowed_ids))
835 @mock.patch.dict(os.environ, {"SATTLE_URI_BASE": "fake_host:1234"})
836 def test_filter_satellites_all_allowed(self):
837 science, matchedTemplate, difference, sources, detectionTask = self._setup_sattle_tests()
839 allowed_ids = list(range(1, 21))
840 response = MockResponse({"allow_list": allowed_ids}, 200, "all allowed")
841 # Run detection and check the results
842 with mock.patch('requests.put', return_value=response):
843 output = detectionTask.run(science, matchedTemplate, difference, sources,
844 idFactory=IdFactory.makeSimple())
846 # Output should be all sources that went in. 20 go in, 20 should come out
847 self.assertEqual(len(output.diaSources), 20)
849 self.assertEqual(set(output.diaSources['id']), set(allowed_ids))
851 @mock.patch.dict(os.environ, {"SATTLE_URI_BASE": "fake_host:1234"})
852 def test_filter_satellites_none_allowed(self):
853 science, matchedTemplate, difference, sources, detectionTask = self._setup_sattle_tests()
855 response = MockResponse({"allow_list": []}, 200, "none allowed")
856 # Run detection and confirm it raises for no diasources
857 with self.assertRaises(detectAndMeasure.NoDiaSourcesError):
858 with mock.patch('requests.put', return_value=response):
859 detectionTask.run(science, matchedTemplate, difference, sources,
860 idFactory=IdFactory.makeSimple())
862 @mock.patch.dict(os.environ, {"SATTLE_URI_BASE": ""})
863 def test_fail_on_sattle_misconfiguration(self):
864 with self.assertRaises(pexConfig.FieldValidationError):
865 self._setup_detection(run_sattle=True)
867 def test_trailed_glints(self):
868 """Test that the glint_trail column works, and that
869 the trailed_glints output contains the expected information.
870 """
871 noiseLevel = 1.
872 staticSeed = 1
873 diffim, diaSources = makeTestImage(seed=staticSeed, noiseLevel=noiseLevel, noiseSeed=6)
874 self._check_values(diaSources['glint_trail'])
876 # Run detection and return the output Struct so we can check it
877 def _detection_wrapper(diffim, diaSources):
878 detectionTask = self._setup_detection()
879 scienceBase, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6)
880 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7)
881 science = scienceBase.clone()
882 science.maskedImage -= diffim.maskedImage
883 difference = science.clone()
884 difference.maskedImage -= matchedTemplate.maskedImage
885 output = detectionTask.run(science, matchedTemplate, difference, sources)
886 return output
888 output = _detection_wrapper(diffim, diaSources)
889 self.assertTrue('slopes' in output.glintTrailInfo)
890 self.assertTrue('intercepts' in output.glintTrailInfo)
891 self.assertTrue('stderrs' in output.glintTrailInfo)
892 self.assertTrue('lengths' in output.glintTrailInfo)
893 self.assertTrue('angles' in output.glintTrailInfo)
896class DetectAndMeasureScoreTest(DetectAndMeasureTestBase, lsst.utils.tests.TestCase):
897 detectionTask = detectAndMeasure.DetectAndMeasureScoreTask
899 def test_detection_xy0(self):
900 """Basic functionality test with non-zero x0 and y0.
901 """
902 # Set up the simulated images
903 noiseLevel = 1.
904 staticSeed = 1
905 fluxLevel = 500
906 kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel, "x0": 12345, "y0": 67890}
907 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
908 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
909 difference = science.clone()
910 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
911 scienceKernel = science.psf.getKernel()
912 score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
914 # Configure the detection Task
915 # Set the subtraction residual metric threshold high, since the template
916 # is not actually subtracted from the science image.
917 detectionTask = self._setup_detection(doDeblend=True, badSubtractionRatioThreshold=1.,
918 doSkySources=False)
920 # Run detection and check the results
921 output = detectionTask.run(science, matchedTemplate, difference, score, sources,
922 idFactory=self.idGenerator.make_table_id_factory())
924 # Catalog ids should be very large from this id generator.
925 self.assertTrue(all(output.diaSources['id'] > 1000000000))
926 subtractedMeasuredExposure = output.subtractedMeasuredExposure
928 self.assertImagesEqual(subtractedMeasuredExposure.image, difference.image)
930 self.assertEqual(len(output.diaSources), len(sources))
931 # no sources should be flagged as negative
932 self.assertEqual(len(~output.diaSources["is_negative"]), len(output.diaSources))
933 # TODO DM-41496: restore this block once we handle detections on edge
934 # pixels better; at least one of these sources currently has a bad
935 # centroid because most of the source is rejected as EDGE.
936 # refIds = []
937 # for source in sources:
938 # self._check_diaSource(output.diaSources, source, refIds=refIds)
940 def test_measurements_finite(self):
941 """Measured fluxes and centroids should always be finite.
942 """
943 columnNames = ["coord_ra", "coord_dec", "ip_diffim_forced_PsfFlux_instFlux",
944 "ip_diffim_forced_template_PsfFlux_instFlux"]
946 # Set up the simulated images
947 noiseLevel = 1.
948 staticSeed = 1
949 transientSeed = 6
950 xSize = 256
951 ySize = 256
952 kwargs = {"psfSize": 2.4, "x0": 0, "y0": 0,
953 "xSize": xSize, "ySize": ySize}
954 science, sources = makeTestImage(seed=staticSeed, noiseLevel=noiseLevel, noiseSeed=6,
955 nSrc=1, **kwargs)
956 matchedTemplate, _ = makeTestImage(seed=staticSeed, noiseLevel=noiseLevel/4, noiseSeed=7,
957 nSrc=1, **kwargs)
958 rng = np.random.RandomState(3)
959 xLoc = np.arange(-5, xSize+5, 10)
960 rng.shuffle(xLoc)
961 yLoc = np.arange(-5, ySize+5, 10)
962 rng.shuffle(yLoc)
963 transients, transientSources = makeTestImage(seed=transientSeed,
964 nSrc=len(xLoc), fluxLevel=1000.,
965 noiseLevel=noiseLevel, noiseSeed=8,
966 xLoc=xLoc, yLoc=yLoc,
967 **kwargs)
968 difference = science.clone()
969 difference.maskedImage -= matchedTemplate.maskedImage
970 difference.maskedImage += transients.maskedImage
971 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
972 scienceKernel = science.psf.getKernel()
973 score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
975 # Configure the detection Task
976 detectionTask = self._setup_detection(doForcedMeasurement=True)
978 # Run detection and check the results
979 output = detectionTask.run(science, matchedTemplate, difference, score, sources)
981 for column in columnNames:
982 self._check_values(output.diaSources[column])
983 self._check_values(output.diaSources.getX(), minValue=0, maxValue=xSize)
984 self._check_values(output.diaSources.getY(), minValue=0, maxValue=ySize)
985 self._check_values(output.diaSources.getPsfInstFlux())
987 def test_detect_transients(self):
988 """Run detection on a difference image containing transients.
989 """
990 # Set up the simulated images
991 noiseLevel = 1.
992 staticSeed = 1
993 transientSeed = 6
994 fluxLevel = 500
995 kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel}
996 scienceBase, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
997 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
998 scienceKernel = scienceBase.psf.getKernel()
999 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
1001 # Configure the detection Task
1002 detectionTask = self._setup_detection(doMerge=False)
1003 kwargs["seed"] = transientSeed
1004 kwargs["nSrc"] = 10
1005 kwargs["fluxLevel"] = 1000
1007 # Run detection and check the results
1008 def _detection_wrapper(positive=True):
1009 """Simulate positive or negative transients and run detection.
1011 Parameters
1012 ----------
1013 positive : `bool`, optional
1014 If set, use positive transient sources.
1015 """
1017 transients, transientSources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, **kwargs)
1018 science = scienceBase.clone()
1019 if positive:
1020 science.maskedImage += transients.maskedImage
1021 else:
1022 science.maskedImage -= transients.maskedImage
1023 difference = science.clone()
1024 difference.maskedImage -= matchedTemplate.maskedImage
1025 score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
1026 # NOTE: NoiseReplacer (run by forcedMeasurement) can modify the
1027 # science image if we've e.g. removed parents post-deblending.
1028 # Pass a clone of the science image, so that it doesn't disrupt
1029 # later tests.
1030 output = detectionTask.run(science, matchedTemplate, difference, score, sources)
1031 refIds = []
1032 scale = 1. if positive else -1.
1033 # sources near the edge may have untrustworthy centroids
1034 goodSrcFlags = ~output.diaSources['base_PixelFlags_flag_edge']
1035 for diaSource, goodSrcFlag in zip(output.diaSources, goodSrcFlags):
1036 if goodSrcFlag:
1037 self._check_diaSource(transientSources, diaSource, refIds=refIds, scale=scale)
1038 _detection_wrapper(positive=True)
1039 _detection_wrapper(positive=False)
1041 def test_detect_dipoles(self):
1042 """Run detection on a difference image containing dipoles.
1043 """
1044 # Set up the simulated images
1045 noiseLevel = 1.
1046 staticSeed = 1
1047 fluxLevel = 1000
1048 fluxRange = 1.5
1049 nSources = 10
1050 offset = 1
1051 xSize = 300
1052 ySize = 300
1053 kernelSize = 32
1054 # Avoid placing sources near the edge for this test, so that we can
1055 # easily check that the correct number of sources are detected.
1056 templateBorderSize = kernelSize//2
1057 dipoleFlag = "ip_diffim_DipoleFit_classification"
1058 kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel, "fluxRange": fluxRange,
1059 "nSrc": nSources, "templateBorderSize": templateBorderSize, "kernelSize": kernelSize,
1060 "xSize": xSize, "ySize": ySize}
1061 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
1062 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
1063 difference = science.clone()
1064 # Shift the template by a pixel in order to make dipoles in the difference image.
1065 matchedTemplate.image.array[...] = np.roll(matchedTemplate.image.array[...], offset, axis=0)
1066 matchedTemplate.variance.array[...] = np.roll(matchedTemplate.variance.array[...], offset, axis=0)
1067 matchedTemplate.mask.array[...] = np.roll(matchedTemplate.mask.array[...], offset, axis=0)
1068 difference.maskedImage -= matchedTemplate.maskedImage[science.getBBox()]
1069 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
1070 scienceKernel = science.psf.getKernel()
1071 score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
1073 detectionTask = self._setup_detection(badSubtractionRatioThreshold=0.3)
1074 output = detectionTask.run(science, matchedTemplate, difference, score, sources)
1075 diaSources = output.diaSources[~output.diaSources["sky_source"]].copy(deep=True)
1076 self.assertEqual(len(diaSources), len(sources))
1077 refIds = []
1078 for diaSource in diaSources:
1079 if diaSource[dipoleFlag]:
1080 self._check_diaSource(sources, diaSource, refIds=refIds, scale=0,
1081 rtol=0.05, atol=None, usePsfFlux=False)
1082 self.assertFloatsAlmostEqual(diaSource["ip_diffim_DipoleFit_orientation"],
1083 -np.pi / 2, atol=2.)
1084 self.assertFloatsAlmostEqual(diaSource["ip_diffim_DipoleFit_separation"], offset, rtol=0.1)
1085 else:
1086 raise ValueError("DiaSource with ID %s is not a dipole!", diaSource.getId())
1088 def test_sky_sources(self):
1089 """Add sky sources and check that they are sufficiently far from other
1090 sources and have negligible flux.
1091 """
1092 # Set up the simulated images
1093 noiseLevel = 1.
1094 staticSeed = 1
1095 transientSeed = 6
1096 transientFluxLevel = 1000.
1097 transientFluxRange = 1.5
1098 fluxLevel = 500
1099 kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel}
1100 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
1101 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
1102 transients, transientSources = makeTestImage(seed=transientSeed, psfSize=2.4,
1103 nSrc=10, fluxLevel=transientFluxLevel,
1104 fluxRange=transientFluxRange,
1105 noiseLevel=noiseLevel, noiseSeed=8)
1106 difference = science.clone()
1107 difference.maskedImage -= matchedTemplate.maskedImage
1108 difference.maskedImage += transients.maskedImage
1109 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
1110 scienceKernel = science.psf.getKernel()
1111 kernelWidth = np.max(scienceKernel.getDimensions())//2
1112 score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
1114 # Configure the detection Task
1115 detectionTask = self._setup_detection(doSkySources=True)
1117 # Run detection and check the results
1118 output = detectionTask.run(science, matchedTemplate, difference, score, sources,
1119 idFactory=self.idGenerator.make_table_id_factory())
1120 nSkySourcesGenerated = detectionTask.metadata["n_skySources"]
1121 skySources = output.diaSources[output.diaSources["sky_source"]]
1122 self.assertEqual(len(skySources), nSkySourcesGenerated)
1123 for skySource in skySources:
1124 # Disable the sky_source flag to enable checks.
1125 skySource["sky_source"] = False
1126 # The sky sources should not be close to any other source
1127 with self.assertRaises(AssertionError):
1128 self._check_diaSource(transientSources, skySource, matchDistance=kernelWidth)
1129 with self.assertRaises(AssertionError):
1130 self._check_diaSource(sources, skySource, matchDistance=kernelWidth)
1131 # The sky sources should have low flux levels.
1132 self._check_diaSource(transientSources, skySource, matchDistance=1000, scale=0.,
1133 atol=np.sqrt(transientFluxRange*transientFluxLevel))
1135 # Catalog ids should be very large from this id generator.
1136 self.assertTrue(all(output.diaSources['id'] > 1000000000))
1138 def test_exclude_mask_detections(self):
1139 """Sources with certain bad mask planes set should not be detected.
1140 """
1141 # Set up the simulated images
1142 noiseLevel = 1.
1143 staticSeed = 1
1144 transientSeed = 6
1145 fluxLevel = 500
1146 radius = 2
1147 kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel}
1148 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
1149 science.getInfo().setVisitInfo(makeVisitInfo())
1150 detector = DetectorWrapper(numAmps=1).detector
1151 science.setDetector(detector)
1152 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
1154 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
1155 scienceKernel = science.psf.getKernel()
1156 # Configure the detection Task
1157 badSourceFlags = ["base_PixelFlags_flag_offimage",
1158 "base_PixelFlags_flag_edgeCenterAll",
1159 "base_PixelFlags_flag_badCenterAll",
1160 "base_PixelFlags_flag_saturatedCenterAll", ]
1161 detectionTask = self._setup_detection(nSkySources=0, badSourceFlags=badSourceFlags)
1162 excludeMaskPlanes = ["EDGE", "BAD", "SAT"]
1163 nBad = len(excludeMaskPlanes)
1164 self.assertGreater(nBad, 0)
1165 kwargs["seed"] = transientSeed
1166 kwargs["nSrc"] = nBad
1167 kwargs["fluxLevel"] = 1000
1169 # Run detection and check the results
1170 def _detection_wrapper(setFlags=True):
1171 transients, transientSources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, **kwargs)
1172 difference = science.clone()
1173 difference.maskedImage -= matchedTemplate.maskedImage
1174 difference.maskedImage += transients.maskedImage
1175 if setFlags:
1176 for src, badMask in zip(transientSources, excludeMaskPlanes):
1177 srcX = int(src.getX())
1178 srcY = int(src.getY())
1179 srcBbox = lsst.geom.Box2I(lsst.geom.Point2I(srcX - radius, srcY - radius),
1180 lsst.geom.Extent2I(2*radius + 1, 2*radius + 1))
1181 difference[srcBbox].mask.array |= lsst.afw.image.Mask.getPlaneBitMask(badMask)
1182 score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
1183 if setFlags:
1184 with self.assertRaises(AlgorithmError):
1185 output = detectionTask.run(science, matchedTemplate, difference, score, sources)
1186 else:
1187 output = detectionTask.run(science, matchedTemplate, difference, score, sources)
1188 refIds = []
1189 goodSrcFlags = checkMask(difference.mask, transientSources, excludeMaskPlanes)
1190 self.assertEqual(np.sum(~goodSrcFlags), 0)
1191 for diaSource, goodSrcFlag in zip(output.diaSources, goodSrcFlags):
1192 if ~goodSrcFlag:
1193 with self.assertRaises(AssertionError):
1194 self._check_diaSource(transientSources, diaSource, refIds=refIds)
1195 else:
1196 self._check_diaSource(transientSources, diaSource, refIds=refIds)
1197 _detection_wrapper(setFlags=False)
1198 _detection_wrapper(setFlags=True)
1201class TestNegativePeaks(lsst.utils.tests.TestCase):
1202 """Tests of deblending and merging negative peaks, to test fixes for the
1203 various problems found on DM-48596.
1204 """
1206 def testDeblendNegatives(self):
1207 """Test that negative peaks get deblended and not destroyed: DM-48704.
1208 This is only a test of deblending, not of merging.
1209 """
1210 # Make a science image with one blend of two positive sources, to
1211 # subtract from an empty template, resulting in a negative diffim.
1212 bbox = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(100, 100))
1213 dataset = lsst.meas.base.tests.TestDataset(bbox)
1214 delta = 10
1215 with dataset.addBlend() as family:
1216 family.addChild(instFlux=2E5, centroid=lsst.geom.Point2D(50, 72))
1217 family.addChild(instFlux=2.5E5, centroid=lsst.geom.Point2D(50+delta, 74))
1218 science, catalog = dataset.realize(noise=1.0,
1219 schema=lsst.meas.base.tests.TestDataset.makeMinimalSchema())
1220 dataset = lsst.meas.base.tests.TestDataset(bbox)
1221 template, _ = dataset.realize(noise=1.0,
1222 schema=lsst.meas.base.tests.TestDataset.makeMinimalSchema())
1223 difference = template.clone()
1224 difference.image -= science.image
1226 config = detectAndMeasure.DetectAndMeasureTask.ConfigClass()
1227 config.doDeblend = True
1228 task = detectAndMeasure.DetectAndMeasureTask(config=config)
1229 # prelude steps taken from `detectAndMeasure.run`
1230 task._prepareInputs(difference)
1231 table = lsst.afw.table.SourceTable.make(task.schema)
1232 results = task.detection.run(table=table, exposure=difference, doSmooth=True)
1233 # Just run the deblend step so we can check the footprints independently.
1234 sources, positives, negatives = task._deblend(difference, results.positive, results.negative)
1236 # DM-48704 fixed a problem where the peaks were in the footprints, but
1237 # the spans were empty.
1238 self.assertEqual(len(negatives), 2)
1239 self.assertEqual(len(positives), 0)
1240 self.assertGreater(negatives[0].getFootprint().getSpans().getArea(), 0)
1241 self.assertGreater(negatives[1].getFootprint().getSpans().getArea(), 0)
1242 # Deblended children are HeavyFootprints; we have to make sure the
1243 # pixel values in those are correct (though DetectAndMeasureTask
1244 # doesn't use the fact that they're Heavy).
1245 # The sources are positive in the science image, and negative in the
1246 # diffim, so the minimum value in the deblended negative footprint is
1247 # the maximum value in the science catalog footprint (ignoring the
1248 # noise in the science and template, hence rtol).
1249 # (catalog[0] is the parent; we want the children)
1250 self.assertFloatsAlmostEqual(negatives[0].getFootprint().getImageArray().min(),
1251 -catalog[1].getFootprint().getImageArray().max(), rtol=1e-4)
1252 self.assertFloatsAlmostEqual(negatives[1].getFootprint().getImageArray().min(),
1253 -catalog[2].getFootprint().getImageArray().max(), rtol=1e-4)
1255 self.assertEqual(sources["is_negative"].sum(), 2)
1257 def testMergeFootprints(self):
1258 """Test that merging footprints does not cause negative ones to
1259 disappear (e.g. get merged into non-connected footprints).
1261 As implemented, the diffim will have three positive sources (one a
1262 blended pair), and 7 negative sources (two a blended pair).
1263 """
1264 # Make a template image with multiple blends, designed to trigger the
1265 # negative-footprint-related bug in `FootprintSet.merge`.
1266 bbox = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(400, 200))
1267 dataset = lsst.meas.base.tests.TestDataset(bbox)
1268 # Because these sources are on the template, they will be negative on
1269 # the diffim, unless the pixels are explicitly made negative via
1270 # `template.image.subset` below.
1271 # positive isolated source on diffim
1272 dataset.addSource(instFlux=.5E5, centroid=lsst.geom.Point2D(25, 26))
1273 # negative isolated source on diffim
1274 dataset.addSource(instFlux=.7E5, centroid=lsst.geom.Point2D(75, 24),
1275 shape=lsst.afw.geom.Quadrupole(12, 7, 2))
1276 delta = 10
1277 # negative blended pair on diffim
1278 with dataset.addBlend() as family:
1279 family.addChild(instFlux=1E5, centroid=lsst.geom.Point2D(150, 72))
1280 family.addChild(instFlux=1.5E5, centroid=lsst.geom.Point2D(150+delta, 74))
1281 # positive blended pair on diffim
1282 with dataset.addBlend() as family:
1283 family.addChild(instFlux=2E5, centroid=lsst.geom.Point2D(250, 72))
1284 family.addChild(instFlux=2.5E5, centroid=lsst.geom.Point2D(250+delta, 74))
1285 # negative blended pair on diffim
1286 with dataset.addBlend() as family:
1287 family.addChild(instFlux=3E5, centroid=lsst.geom.Point2D(350, 72))
1288 family.addChild(instFlux=3.5E5, centroid=lsst.geom.Point2D(350+delta, 74))
1289 # negative isolated source on diffim
1290 dataset.addSource(instFlux=4E5, centroid=lsst.geom.Point2D(75, 124))
1291 # negative isolated source on diffim
1292 dataset.addSource(instFlux=5E5, centroid=lsst.geom.Point2D(175, 124))
1293 schema = lsst.meas.base.tests.TestDataset.makeMinimalSchema()
1294 schema.addField("sky_source", "Flag", "Sky source.")
1295 template, catalog = dataset.realize(noise=100.0, schema=schema)
1296 # These will be positive sources on the diffim (as noted above).
1297 template.image.subset(lsst.geom.Box2I(lsst.geom.Point2I(15, 15),
1298 lsst.geom.Point2I(35, 35))).array *= -1
1299 template.image.subset(lsst.geom.Box2I(lsst.geom.Point2I(230, 60),
1300 lsst.geom.Point2I(270, 85))).array *= -1
1301 dataset = lsst.meas.base.tests.TestDataset(bbox)
1302 science, _ = dataset.realize(noise=100.0, schema=schema)
1303 difference = science.clone()
1304 difference.image -= template.image
1306 config = detectAndMeasure.DetectAndMeasureTask.ConfigClass()
1307 config.doDeblend = True
1308 config.raiseOnBadSubtractionRatio = False
1309 config.raiseOnNoDiaSources = False
1310 task = detectAndMeasure.DetectAndMeasureTask(config=config)
1311 result = task.run(science, template, difference, catalog)
1313 # The original bug manifested as the (175,124) source disappaering in
1314 # the merge step, and being included in the peaks of a duplicated
1315 # (75,124) source, even though their footprints are not connected.
1316 mask = np.isclose(result.diaSources["slot_Centroid_x"], 75) & \
1317 np.isclose(result.diaSources["slot_Centroid_y"], 124)
1318 self.assertEqual(mask.sum(), 1)
1319 peaks = result.diaSources[mask][0].getFootprint().peaks
1320 self.assertEqual(len(peaks), 1)
1321 self.assertEqual(peaks[0].getIx(), 75)
1322 self.assertEqual(peaks[0].getIy(), 124)
1324 mask = np.isclose(result.diaSources["slot_Centroid_x"], 175) & \
1325 np.isclose(result.diaSources["slot_Centroid_y"], 124)
1326 self.assertEqual(mask.sum(), 1)
1327 peaks = result.diaSources[mask][0].getFootprint().peaks
1328 self.assertEqual(len(peaks), 1)
1329 self.assertEqual(peaks[0].getIx(), 175)
1330 self.assertEqual(peaks[0].getIy(), 124)
1332 # This checks that all the returned footprints have exactly one peak;
1333 # it's a more general test than the above, but I'm keeping that for
1334 # the more explicit test of the original bugged sources.
1335 peak_x = np.array([peak['f_x'] for src in result.diaSources for peak in src.getFootprint().peaks])
1336 peak_y = np.array([peak['f_y'] for src in result.diaSources for peak in src.getFootprint().peaks])
1337 peaks = np.column_stack((peak_x, peak_y))
1338 unique_peaks, counts = np.unique(peaks, axis=0, return_counts=True)
1339 self.assertEqual(np.sum(counts > 1), 0)
1341 # There are three positive sources that should not have `is_negative`
1342 # set, independent of how deblending/merging of negative sources is
1343 # handled.
1344 self.assertEqual((~result.diaSources["is_negative"]).sum(), 3)
1347def makeVisitInfo(id=1):
1348 """Return a non-NaN visitInfo."""
1349 return afwImage.VisitInfo(id=id,
1350 exposureTime=10.01,
1351 darkTime=11.02,
1352 date=dafBase.DateTime(65321.1, dafBase.DateTime.MJD, dafBase.DateTime.TAI),
1353 ut1=12345.1,
1354 era=45.1*geom.degrees,
1355 boresightRaDec=geom.SpherePoint(23.1, 73.2, geom.degrees),
1356 boresightAzAlt=geom.SpherePoint(134.5, 33.3, geom.degrees),
1357 boresightAirmass=1.73,
1358 boresightRotAngle=73.2*geom.degrees,
1359 rotType=afwImage.RotType.SKY,
1360 observatory=Observatory(
1361 11.1*geom.degrees, 22.2*geom.degrees, 0.333),
1362 weather=Weather(1.1, 2.2, 34.5),
1363 )
1366class MockResponse:
1367 """Provide a mock for requests.put calls"""
1368 def __init__(self, json_data, status_code, text):
1369 self.json_data = json_data
1370 self.status_code = status_code
1371 self.text = text
1373 def json(self):
1374 return self.json_data
1376 def raise_for_status(self):
1377 if self.status_code != 200:
1378 raise requests.exceptions.HTTPError
1381def setup_module(module):
1382 lsst.utils.tests.init()
1385class MemoryTestCase(lsst.utils.tests.MemoryTestCase):
1386 pass
1389if __name__ == "__main__": 1389 ↛ 1390line 1389 didn't jump to line 1390 because the condition on line 1389 was never true
1390 lsst.utils.tests.init()
1391 unittest.main()