Coverage for tests/test_detectAndMeasure.py: 99%
968 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-07-26 09:34 +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.detection as afwDetection
29import lsst.afw.geom as afwGeom
30import lsst.afw.image as afwImage
31import lsst.afw.table as afwTable
32import lsst.afw.math as afwMath
33import lsst.geom
34from lsst.ip.diffim import detectAndMeasure, subtractImages
35from lsst.afw.table import IdFactory
36from lsst.afw.cameraGeom.testUtils import DetectorWrapper
37import lsst.meas.algorithms as measAlg
38from lsst.pipe.base import InvalidQuantumError, UpstreamFailureNoWorkFound, AlgorithmError
39import lsst.utils.tests
40import lsst.meas.base.tests
41import lsst.daf.base as dafBase
42from lsst.afw.coord import Observatory, Weather
43import lsst.geom as geom
44import lsst.pex.config as pexConfig
46from utils import makeTestImage, checkMask
49class DetectAndMeasureTestBase:
51 def _check_diaSource(self, refSources, diaSource, refIds=None,
52 matchDistance=1., scale=1., usePsfFlux=True,
53 rtol=0.025, atol=None):
54 """Match a diaSource with a source in a reference catalog
55 and compare properties.
57 Parameters
58 ----------
59 refSources : `lsst.afw.table.SourceCatalog`
60 The reference catalog.
61 diaSource : `lsst.afw.table.SourceRecord`
62 The new diaSource to match to the reference catalog.
63 refIds : `list` of `int`, optional
64 Source IDs of previously associated diaSources.
65 matchDistance : `float`, optional
66 Maximum distance allowed between the detected and reference source
67 locations, in pixels.
68 scale : `float`, optional
69 Optional factor to scale the flux by before performing the test.
70 usePsfFlux : `bool`, optional
71 If set, test the PsfInstFlux field, otherwise use ApInstFlux.
72 rtol : `float`, optional
73 Relative tolerance of the flux value test.
74 atol : `float`, optional
75 Absolute tolerance of the flux value test.
76 """
77 if diaSource["sky_source"]:
78 return
79 distance = np.sqrt((diaSource.getX() - refSources.getX())**2
80 + (diaSource.getY() - refSources.getY())**2)
81 self.assertLess(min(distance), matchDistance)
82 src = refSources[np.argmin(distance)]
83 if refIds is not None:
84 # Check that the same source was not previously associated
85 self.assertNotIn(src.getId(), refIds)
86 refIds.append(src.getId())
87 if atol is None:
88 atol = rtol*src.getPsfInstFlux() if usePsfFlux else rtol*src.getApInstFlux()
89 if usePsfFlux:
90 self.assertFloatsAlmostEqual(src.getPsfInstFlux()*scale, diaSource.getPsfInstFlux(),
91 rtol=rtol, atol=atol)
92 else:
93 self.assertFloatsAlmostEqual(src.getApInstFlux()*scale, diaSource.getApInstFlux(),
94 rtol=rtol, atol=atol)
96 def _check_values(self, values, minValue=None, maxValue=None):
97 """Verify that an array has finite values, and optionally that they are
98 within specified minimum and maximum bounds.
100 Parameters
101 ----------
102 values : `numpy.ndarray`
103 Array of values to check.
104 minValue : `float`, optional
105 Minimum allowable value.
106 maxValue : `float`, optional
107 Maximum allowable value.
108 """
109 self.assertTrue(np.all(np.isfinite(values)))
110 if minValue is not None:
111 self.assertTrue(np.all(values >= minValue))
112 if maxValue is not None:
113 self.assertTrue(np.all(values <= maxValue))
115 def _setup_detection(self, doSkySources=True, nSkySources=5,
116 doSubtractBackground=False, run_sattle=False,
117 badSourceFlags=None, **kwargs):
118 """Setup and configure the detection and measurement PipelineTask.
120 Parameters
121 ----------
122 doSkySources : `bool`, optional
123 Generate sky sources.
124 nSkySources : `int`, optional
125 The number of sky sources to add in isolated background regions.
126 **kwargs
127 Any additional config parameters to set.
129 Returns
130 -------
131 `lsst.pipe.base.PipelineTask`
132 The configured Task to use for detection and measurement.
133 """
134 config = self.detectionTask.ConfigClass()
135 config.doSkySources = doSkySources
136 if doSkySources:
137 config.skySources.nSources = nSkySources
138 config.update(**kwargs)
140 config.run_sattle = run_sattle
142 # Make a realistic id generator so that output catalog ids are useful.
143 dataId = lsst.daf.butler.DataCoordinate.standardize(
144 instrument="I",
145 visit=42,
146 detector=12,
147 universe=lsst.daf.butler.DimensionUniverse(),
148 )
149 if badSourceFlags is None:
150 badSourceFlags = ["base_PixelFlags_flag_offimage", ]
151 config.idGenerator.packer.name = "observation"
152 config.idGenerator.packer["observation"].n_observations = 10000
153 config.idGenerator.packer["observation"].n_detectors = 99
154 config.idGenerator.n_releases = 8
155 config.idGenerator.release_id = 2
156 config.doSubtractBackground = doSubtractBackground
157 config.badSourceFlags = badSourceFlags
158 self.idGenerator = config.idGenerator.apply(dataId)
160 return self.detectionTask(config=config)
163class DetectAndMeasureTest(DetectAndMeasureTestBase, lsst.utils.tests.TestCase):
164 detectionTask = detectAndMeasure.DetectAndMeasureTask
166 def test_detection_xy0(self):
167 """Basic functionality test with non-zero x0 and y0.
168 """
169 # Set up the simulated images
170 noiseLevel = 1.
171 staticSeed = 1
172 fluxLevel = 500
173 psfSize = 2.4
174 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel, "x0": 12345, "y0": 67890}
175 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
176 science.getInfo().setVisitInfo(makeVisitInfo())
177 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
178 difference = science.clone()
180 # Configure the detection Task
181 # Set the subtraction residual metric threshold high, since the template
182 # is not actually subtracted from the science image.
183 detectionTask = self._setup_detection(doDeblend=True, badSubtractionRatioThreshold=1.,
184 doSkySources=False)
186 # Run detection and check the results
187 output = detectionTask.run(science, matchedTemplate, difference, sources,
188 idFactory=self.idGenerator.make_table_id_factory())
189 subtractedMeasuredExposure = output.subtractedMeasuredExposure
191 # Catalog ids should be very large from this id generator.
192 self.assertTrue(all(output.diaSources['id'] > 1000000000))
193 self.assertImagesEqual(subtractedMeasuredExposure.image, difference.image)
195 # all of the sources should have been detected
196 self.assertEqual(len(output.diaSources), len(sources))
197 # no sources should be flagged as negative
198 self.assertEqual(len(~output.diaSources["is_negative"]), len(output.diaSources))
199 refIds = []
200 for source in sources:
201 self._check_diaSource(output.diaSources, source, refIds=refIds)
203 def test_raise_bad_psf(self):
204 """Detection should raise if the PSF width is NaN
205 """
206 # Set up the simulated images
207 noiseLevel = 1.
208 staticSeed = 1
209 fluxLevel = 500
210 psfSize = 2.4
211 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel, "x0": 12345, "y0": 67890}
212 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
213 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
214 difference = science.clone()
216 psfShape = difference.getPsf().computeBBox(lsst.geom.Point2D(0, 0)).getDimensions()
217 psfcI = afwImage.ImageD(lsst.geom.Extent2I(psfShape[1], psfShape[0]))
218 psfNew = np.zeros(psfShape)
219 psfNew[0, :] = 1
220 psfcI.array = psfNew
221 psfcK = afwMath.FixedKernel(psfcI)
222 psf = measAlg.KernelPsf(psfcK)
223 difference.setPsf(psf)
225 # Configure the detection Task
226 detectionTask = self._setup_detection()
228 # Run detection and check the results
229 with self.assertRaises(UpstreamFailureNoWorkFound):
230 detectionTask.run(science, matchedTemplate, difference, sources)
232 def test_measurements_finite(self):
233 """Measured fluxes and centroids should always be finite.
234 """
235 columnNames = ["coord_ra", "coord_dec", "ip_diffim_forced_PsfFlux_instFlux",
236 "ip_diffim_forced_template_PsfFlux_instFlux"]
238 # Set up the simulated images
239 noiseLevel = 1.
240 staticSeed = 1
241 transientSeed = 6
242 xSize = 256
243 ySize = 256
244 psfSize = 2.4
245 kwargs = {"psfSize": psfSize, "x0": 0, "y0": 0,
246 "xSize": xSize, "ySize": ySize}
247 science, sources = makeTestImage(seed=staticSeed, noiseLevel=noiseLevel, noiseSeed=6,
248 nSrc=1, **kwargs)
249 science.getInfo().setVisitInfo(makeVisitInfo())
250 matchedTemplate, _ = makeTestImage(seed=staticSeed, noiseLevel=noiseLevel/4, noiseSeed=7,
251 nSrc=1, **kwargs)
252 rng = np.random.RandomState(3)
253 xLoc = np.arange(-5, xSize+5, 10)
254 rng.shuffle(xLoc)
255 yLoc = np.arange(-5, ySize+5, 10)
256 rng.shuffle(yLoc)
257 transients, transientSources = makeTestImage(seed=transientSeed,
258 nSrc=len(xLoc), fluxLevel=1000.,
259 noiseLevel=noiseLevel, noiseSeed=8,
260 xLoc=xLoc, yLoc=yLoc,
261 **kwargs)
262 difference = science.clone()
263 difference.maskedImage -= matchedTemplate.maskedImage
264 difference.maskedImage += transients.maskedImage
266 # Configure the detection Task
267 detectionTask = self._setup_detection(doForcedMeasurement=True)
269 # Run detection and check the results
270 output = detectionTask.run(science, matchedTemplate, difference, sources)
272 for column in columnNames:
273 self._check_values(output.diaSources[column])
274 self._check_values(output.diaSources.getX(), minValue=0, maxValue=xSize)
275 self._check_values(output.diaSources.getY(), minValue=0, maxValue=ySize)
276 self._check_values(output.diaSources.getPsfInstFlux())
278 def test_raise_config_schema_mismatch(self):
279 """Check that sources with specified flags are removed from the catalog.
280 """
281 # Configure the detection Task, and and set a config that is not in the schema
282 with self.assertRaises(InvalidQuantumError):
283 self._setup_detection(badSourceFlags=["Bogus_flag_42"])
285 def test_remove_unphysical(self):
286 """Check that sources with specified flags are removed from the catalog.
287 """
288 # Set up the simulated images
289 noiseLevel = 1.
290 staticSeed = 1
291 xSize = 256
292 ySize = 256
293 psfSize = 2.4
294 kwargs = {"psfSize": psfSize, "xSize": xSize, "ySize": ySize}
295 science, sources = makeTestImage(seed=staticSeed, noiseLevel=noiseLevel, noiseSeed=6,
296 nSrc=1, **kwargs)
297 science.getInfo().setVisitInfo(makeVisitInfo())
298 matchedTemplate, _ = makeTestImage(seed=staticSeed, noiseLevel=noiseLevel/4, noiseSeed=7,
299 nSrc=1, **kwargs)
300 transients, transientSources = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=8, nSrc=1, **kwargs)
301 science.maskedImage += transients.maskedImage
302 difference = science.clone()
303 bbox = difference.getBBox()
304 difference.maskedImage -= matchedTemplate.maskedImage
306 # Configure the detection Task, and remove unphysical sources
307 detectionTask = self._setup_detection(doForcedMeasurement=False, doSkySources=True, nSkySources=20,
308 badSourceFlags=["base_PixelFlags_flag_offimage", ])
310 # Run detection and check the results
311 diaSources = detectionTask.run(science, matchedTemplate, difference, sources).diaSources
312 badDiaSrcDoRemove = ~bbox.contains(diaSources.getX(), diaSources.getY())
313 nBadDoRemove = np.count_nonzero(badDiaSrcDoRemove)
314 # Verify that all sources are physical
315 self.assertEqual(nBadDoRemove, 0)
316 # Set a few centroids outside the image bounding box
317 nSetBad = 5
318 for src in diaSources[0: nSetBad]:
319 src["slot_Centroid_x"] += xSize
320 src["slot_Centroid_y"] += ySize
321 src["base_PixelFlags_flag_offimage"] = True
322 # Verify that these sources are outside the image
323 badDiaSrc = ~bbox.contains(diaSources.getX(), diaSources.getY())
324 nBad = np.count_nonzero(badDiaSrc)
325 self.assertEqual(nBad, nSetBad)
326 diaSourcesNoBad = detectionTask._removeBadSources(diaSources)
327 badDiaSrcNoBad = ~bbox.contains(diaSourcesNoBad.getX(), diaSourcesNoBad.getY())
329 # Verify that no sources outside the image bounding box remain
330 self.assertEqual(np.count_nonzero(badDiaSrcNoBad), 0)
331 self.assertEqual(len(diaSourcesNoBad), len(diaSources) - nSetBad)
333 def test_detect_transients(self):
334 """Run detection on a difference image containing transients.
335 """
336 # Set up the simulated images
337 noiseLevel = 1.
338 staticSeed = 1
339 transientSeed = 6
340 nTransients = 10
341 transientFlux = 1000
342 fluxLevel = 500
343 psfSize = 2.4
344 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel}
345 scienceBase, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
346 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
348 # Configure the detection Task
349 detectionTask = self._setup_detection(doMerge=False, doSkySources=True)
350 kwargs["seed"] = transientSeed
351 kwargs["nSrc"] = nTransients
352 kwargs["fluxLevel"] = transientFlux
354 # Run detection and check the results
355 def _run_and_check_detections(positive=True):
356 transients, transientSources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, **kwargs)
357 science = scienceBase.clone()
358 if positive:
359 science.maskedImage += transients.maskedImage
360 else:
361 science.maskedImage -= transients.maskedImage
362 difference = science.clone()
363 difference.maskedImage -= matchedTemplate.maskedImage
364 # NOTE: NoiseReplacer (run by forcedMeasurement) can modify the
365 # science image if we've e.g. removed parents post-deblending.
366 # Pass a clone of the science image, so that it doesn't disrupt
367 # later tests.
368 output = detectionTask.run(science, matchedTemplate, difference, sources)
369 refIds = []
370 scale = 1. if positive else -1.
371 for diaSource in output.diaSources:
372 self._check_diaSource(transientSources, diaSource, refIds=refIds, scale=scale)
373 _run_and_check_detections(positive=True)
374 _run_and_check_detections(positive=False)
376 def test_detect_transients_with_background(self):
377 """Run detection on a difference image containing transients and a background.
378 """
379 # Set up the simulated images
380 noiseLevel = 1.
381 staticSeed = 1
382 transientSeed = 6
383 nTransients = 10
384 transientFlux = 1000
385 fluxLevel = 500
386 xSize = 512
387 ySize = 512
388 x0 = 123
389 y0 = 456
390 psfSize = 2.4
391 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel,
392 "xSize": xSize, "ySize": ySize, "x0": x0, "y0": y0}
393 params = [2.2, 2.1, 2.0, 1.2, 1.1, 1.0]
395 bbox2D = lsst.geom.Box2D(lsst.geom.Point2D(x0, y0), lsst.geom.Extent2D(xSize, ySize))
396 background_model = afwMath.Chebyshev1Function2D(params, bbox2D)
397 scienceBase, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6,
398 background=background_model, **kwargs)
399 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
401 # Configure the detection Task
402 detectionTask = self._setup_detection(doMerge=False, doSubtractBackground=True)
403 kwargs["seed"] = transientSeed
404 kwargs["nSrc"] = nTransients
405 kwargs["fluxLevel"] = transientFlux
407 # Run detection and check the results
408 def _run_and_check_detections(positive=True):
409 transients, transientSources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, **kwargs)
410 science = scienceBase.clone()
411 if positive:
412 science.maskedImage += transients.maskedImage
413 else:
414 science.maskedImage -= transients.maskedImage
415 difference = science.clone()
416 difference.maskedImage -= matchedTemplate.maskedImage
417 # NOTE: NoiseReplacer (run by forcedMeasurement) can modify the
418 # science image if we've e.g. removed parents post-deblending.
419 # Pass a clone of the science image, so that it doesn't disrupt
420 # later tests.
421 output = detectionTask.run(science, matchedTemplate, difference, sources)
422 refIds = []
423 scale = 1. if positive else -1.
424 for transient in transientSources:
425 self._check_diaSource(output.diaSources, transient, refIds=refIds, scale=scale)
426 _run_and_check_detections(positive=True)
427 _run_and_check_detections(positive=False)
429 def test_mask_cosmic_rays(self):
430 """Run detection on a difference image containing a cosmic ray.
431 """
432 # Set up the simulated images
433 noiseLevel = 1.
434 staticSeed = 1
435 fluxLevel = 500
436 xSize = 400
437 ySize = 400
438 psfSize = 2.4
439 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel,
440 "xSize": xSize, "ySize": ySize}
441 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
442 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
443 crMask = science.mask.getPlaneBitMask("CR")
445 # Configure the detection Task
446 detectionTask = self._setup_detection(doMerge=False, doSkySources=True)
448 # Test that no CRs are detected
449 transients, transientSources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, nSrc=10, **kwargs)
450 science.maskedImage += transients.maskedImage
451 difference = science.clone()
452 difference.maskedImage -= matchedTemplate.maskedImage
453 output = detectionTask.run(science, matchedTemplate, difference, sources)
454 crMaskSet = (output.subtractedMeasuredExposure.mask.array & crMask) > 0
455 self.assertTrue(np.all(crMaskSet == 0))
457 crX0 = round(sources.getX()[0] - science.getX0())
458 crY0 = round(sources.getY()[0] - science.getY0())
459 crX1 = round(sources.getX()[1] - science.getX0())
460 crY1 = round(sources.getY()[1] - science.getY0())
461 # Add CR-like shape and check that CR is detected
462 # Pick two locations on top of sources, since that is what is likely to
463 # be missed in the first stage of CR rejection.
464 crMaskSetInput = np.zeros(difference.image.array.shape, bool)
465 crMaskSetInput[crX0:crX0+1, crY0:crY0+5] = True
466 crMaskSetInput[crX1:crX1+5, crY1:crY1+1] = True
467 difference.image.array[crMaskSetInput] += 1234
468 output = detectionTask.run(science, matchedTemplate, difference)
469 crMaskSet = (output.subtractedMeasuredExposure.mask.array & crMask) > 0
470 self.assertFalse(np.any(crMaskSet[~crMaskSetInput]))
471 self.assertTrue(np.all(crMaskSet[crMaskSetInput]))
473 def test_missing_mask_planes(self):
474 """Check that detection runs with missing mask planes.
475 """
476 # Set up the simulated images
477 noiseLevel = 1.
478 fluxLevel = 500
479 psfSize = 2.4
480 kwargs = {"psfSize": psfSize, "fluxLevel": fluxLevel, "addMaskPlanes": []}
481 # Use different seeds for the science and template so every source is a diaSource
482 science, sources = makeTestImage(seed=5, noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
483 science.getInfo().setVisitInfo(makeVisitInfo())
484 matchedTemplate, _ = makeTestImage(seed=6, noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
486 difference = science.clone()
487 difference.maskedImage -= matchedTemplate.maskedImage
488 detectionTask = self._setup_detection(raiseOnBadSubtractionRatio=False, raiseOnNoDiaSources=False)
490 # Verify that detection runs without errors
491 detectionTask.run(science, matchedTemplate, difference, sources)
493 def test_detect_dipoles(self):
494 """Run detection on a difference image containing dipoles.
495 """
496 # Set up the simulated images
497 noiseLevel = 1.
498 staticSeed = 1
499 fluxLevel = 1000
500 fluxRange = 1.5
501 nSources = 10
502 offset = 1
503 xSize = 300
504 ySize = 300
505 kernelSize = 31
506 psfSize = 2.4
507 # Avoid placing sources near the edge for this test, so that we can
508 # easily check that the correct number of sources are detected.
509 templateBorderSize = kernelSize//2
510 dipoleFlag = "ip_diffim_DipoleFit_classification"
511 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel, "fluxRange": fluxRange,
512 "nSrc": nSources, "templateBorderSize": templateBorderSize, "kernelSize": kernelSize,
513 "xSize": xSize, "ySize": ySize}
514 dipoleFlag = "ip_diffim_DipoleFit_classification"
515 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
516 science.getInfo().setVisitInfo(makeVisitInfo())
517 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
518 difference = science.clone()
519 matchedTemplate.image.array[...] = np.roll(matchedTemplate.image.array[...], offset, axis=0)
520 matchedTemplate.variance.array[...] = np.roll(matchedTemplate.variance.array[...], offset, axis=0)
521 matchedTemplate.mask.array[...] = np.roll(matchedTemplate.mask.array[...], offset, axis=0)
522 difference.maskedImage -= matchedTemplate.maskedImage[science.getBBox()]
524 # Need a higher residual metric threshold, since we are purposely
525 # creating poor subtractions.
526 detectionTask = self._setup_detection(doMerge=True, badSubtractionRatioThreshold=0.3)
527 output = detectionTask.run(science, matchedTemplate, difference, sources)
528 diaSources = output.diaSources[~output.diaSources["sky_source"]].copy(deep=True)
529 self.assertEqual(len(diaSources), len(sources))
530 # no sources should be flagged as negative
531 self.assertEqual(len(~diaSources["is_negative"]), len(diaSources))
532 refIds = []
533 for diaSource in diaSources:
534 if diaSource[dipoleFlag]: 534 ↛ 541line 534 didn't jump to line 541 because the condition on line 534 was always true
535 self._check_diaSource(sources, diaSource, refIds=refIds, scale=0,
536 rtol=0.05, atol=None, usePsfFlux=False)
537 self.assertFloatsAlmostEqual(diaSource["ip_diffim_DipoleFit_orientation"],
538 -np.pi / 2, atol=2.)
539 self.assertFloatsAlmostEqual(diaSource["ip_diffim_DipoleFit_separation"], offset, rtol=0.1)
540 else:
541 raise ValueError("DiaSource with ID %s is not a dipole!", diaSource.getId())
543 def test_sky_sources(self):
544 """Add sky sources and check that they are sufficiently far from other
545 sources and have negligible flux.
546 """
547 # Set up the simulated images
548 noiseLevel = 1.
549 staticSeed = 1
550 transientSeed = 6
551 transientFluxLevel = 1000.
552 transientFluxRange = 1.5
553 fluxLevel = 500
554 psfSize = 2.4
555 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel}
556 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
557 science.getInfo().setVisitInfo(makeVisitInfo())
558 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
559 transients, transientSources = makeTestImage(seed=transientSeed, psfSize=2.4,
560 nSrc=10, fluxLevel=transientFluxLevel,
561 fluxRange=transientFluxRange,
562 noiseLevel=noiseLevel, noiseSeed=8)
563 difference = science.clone()
564 difference.maskedImage -= matchedTemplate.maskedImage
565 difference.maskedImage += transients.maskedImage
566 kernelWidth = np.max(science.psf.getKernel().getDimensions())//2
568 # Configure the detection Task
569 detectionTask = self._setup_detection(doSkySources=True)
571 # Run detection and check the results
572 output = detectionTask.run(science, matchedTemplate, difference, sources,
573 idFactory=self.idGenerator.make_table_id_factory())
574 skySources = output.diaSources[output.diaSources["sky_source"]]
575 self.assertEqual(len(skySources), detectionTask.config.skySources.nSources)
576 for skySource in skySources:
577 # Disable the sky_source flag to enable checks.
578 skySource["sky_source"] = False
579 # The sky sources should not be close to any other source
580 with self.assertRaises(AssertionError):
581 self._check_diaSource(transientSources, skySource, matchDistance=kernelWidth)
582 with self.assertRaises(AssertionError):
583 self._check_diaSource(sources, skySource, matchDistance=kernelWidth)
584 # The sky sources should have low flux levels.
585 self._check_diaSource(transientSources, skySource, matchDistance=1000, scale=0.,
586 atol=np.sqrt(transientFluxRange*transientFluxLevel))
588 # Catalog ids should be very large from this id generator.
589 self.assertTrue(all(output.diaSources['id'] > 1000000000))
591 def test_exclude_mask_detections(self):
592 """Sources with certain bad mask planes set should not be detected.
593 """
594 # Set up the simulated images
595 noiseLevel = 1.
596 staticSeed = 1
597 transientSeed = 6
598 fluxLevel = 500
599 radius = 2
600 psfSize = 2.4
601 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel}
602 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
603 science.getInfo().setVisitInfo(makeVisitInfo())
604 detector = DetectorWrapper(numAmps=1).detector
605 science.setDetector(detector)
606 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
608 # Configure the detection Task
609 badSourceFlags = ["base_PixelFlags_flag_offimage",
610 "base_PixelFlags_flag_edgeCenterAll",
611 "base_PixelFlags_flag_badCenterAll",
612 "base_PixelFlags_flag_saturatedCenterAll", ]
613 detectionTask = self._setup_detection(nSkySources=0, badSourceFlags=badSourceFlags)
614 excludeMaskPlanes = ["EDGE", "BAD", "SAT"]
615 nBad = len(excludeMaskPlanes)
616 self.assertGreater(nBad, 0)
617 kwargs["seed"] = transientSeed
618 kwargs["nSrc"] = nBad
619 kwargs["fluxLevel"] = 1000
621 # Run detection and check the results
622 def _run_and_check_detections(setFlags=True):
623 transients, transientSources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, **kwargs)
624 difference = science.clone()
625 difference.maskedImage -= matchedTemplate.maskedImage
626 difference.maskedImage += transients.maskedImage
627 if setFlags:
628 for src, badMask in zip(transientSources, excludeMaskPlanes):
629 srcX = int(src.getX())
630 srcY = int(src.getY())
631 srcBbox = lsst.geom.Box2I(lsst.geom.Point2I(srcX - radius, srcY - radius),
632 lsst.geom.Extent2I(2*radius + 1, 2*radius + 1))
633 difference[srcBbox].mask.array |= lsst.afw.image.Mask.getPlaneBitMask(badMask)
634 with self.assertRaises(AlgorithmError):
635 output = detectionTask.run(science, matchedTemplate, difference, sources)
636 else:
637 output = detectionTask.run(science, matchedTemplate, difference, sources)
638 refIds = []
639 goodSrcFlags = checkMask(difference.mask, transientSources, excludeMaskPlanes)
640 self.assertEqual(np.sum(~goodSrcFlags), 0)
641 for diaSource, goodSrcFlag in zip(output.diaSources, goodSrcFlags):
642 if ~goodSrcFlag: 642 ↛ 643line 642 didn't jump to line 643 because the condition on line 642 was never true
643 with self.assertRaises(AssertionError):
644 self._check_diaSource(transientSources, diaSource, refIds=refIds)
645 else:
646 self._check_diaSource(transientSources, diaSource, refIds=refIds)
647 _run_and_check_detections(setFlags=False)
648 _run_and_check_detections(setFlags=True)
650 def test_fake_mask_plane_propagation(self):
651 """Test that we have the mask planes related to fakes in diffim images.
652 This is testing method called updateMasks
653 """
654 xSize = 256
655 ySize = 256
656 science, sources = makeTestImage(psfSize=2.4, xSize=xSize, ySize=ySize, doApplyCalibration=True)
657 science_fake_img, science_fake_sources = makeTestImage(
658 psfSize=2.4, xSize=xSize, ySize=ySize, seed=5, nSrc=3, noiseLevel=0.25, fluxRange=1
659 )
660 template, _ = makeTestImage(psfSize=2.4, xSize=xSize, ySize=ySize, doApplyCalibration=True)
661 tmplt_fake_img, tmplt_fake_sources = makeTestImage(
662 psfSize=2.4, xSize=xSize, ySize=ySize, seed=9, nSrc=3, noiseLevel=0.25, fluxRange=1
663 )
664 # created fakes and added them to the images
665 science.image += science_fake_img.image
666 template.image += tmplt_fake_img.image
668 # TODO: DM-40796 update to INJECTED names when source injection gets refactored
669 # adding mask planes to both science and template images
670 science.mask.addMaskPlane("FAKE")
671 science_fake_bitmask = science.mask.getPlaneBitMask("FAKE")
672 template.mask.addMaskPlane("FAKE")
673 template_fake_bitmask = template.mask.getPlaneBitMask("FAKE")
675 # makeTestImage sets the DETECTED plane on the sources; we can use
676 # that to set the FAKE plane on the science and template images.
677 detected = science_fake_img.mask.getPlaneBitMask("DETECTED")
678 fake_pixels = (science_fake_img.mask.array & detected).nonzero()
679 science.mask.array[fake_pixels] |= science_fake_bitmask
680 detected = tmplt_fake_img.mask.getPlaneBitMask("DETECTED")
681 fake_pixels = (tmplt_fake_img.mask.array & detected).nonzero()
682 template.mask.array[fake_pixels] |= science_fake_bitmask
684 science_fake_masked = (science.mask.array & science_fake_bitmask) > 0
685 template_fake_masked = (template.mask.array & template_fake_bitmask) > 0
687 subtractConfig = subtractImages.AlardLuptonSubtractTask.ConfigClass()
688 subtractConfig.sourceSelector.signalToNoise.fluxField = "truth_instFlux"
689 subtractConfig.sourceSelector.signalToNoise.errField = "truth_instFluxErr"
690 subtractTask = subtractImages.AlardLuptonSubtractTask(config=subtractConfig)
691 subtraction = subtractTask.run(template, science, sources)
693 # check subtraction mask plane is set where we set the previous masks
694 diff_mask = subtraction.difference.mask
696 # science mask should be now in INJECTED
697 inj_masked = (diff_mask.array & diff_mask.getPlaneBitMask("INJECTED")) > 0
699 # template mask should be now in INJECTED_TEMPLATE
700 injTmplt_masked = (diff_mask.array & diff_mask.getPlaneBitMask("INJECTED_TEMPLATE")) > 0
702 self.assertFloatsEqual(inj_masked.astype(int), science_fake_masked.astype(int))
703 # The template is convolved, so the INJECTED_TEMPLATE mask plane may
704 # include more pixels than the FAKE mask plane
705 injTmplt_masked &= template_fake_masked
706 self.assertFloatsEqual(injTmplt_masked.astype(int), template_fake_masked.astype(int))
708 # Now check that detection of fakes have the correct flag for injections
709 detectionTask = self._setup_detection()
711 output = detectionTask.run(subtraction.matchedScience,
712 subtraction.matchedTemplate,
713 subtraction.difference,
714 subtraction.kernelSources)
715 diaSources = output.diaSources[~output.diaSources["sky_source"]].copy(deep=True)
717 sci_refIds = []
718 tmpl_refIds = []
719 for diaSrc in diaSources:
720 if diaSrc['base_PsfFlux_instFlux'] > 0:
721 self._check_diaSource(science_fake_sources, diaSrc, scale=1, refIds=sci_refIds)
722 self.assertTrue(diaSrc['base_PixelFlags_flag_injected'])
723 self.assertTrue(diaSrc['base_PixelFlags_flag_injectedCenter'])
724 self.assertFalse(diaSrc['base_PixelFlags_flag_injected_template'])
725 self.assertFalse(diaSrc['base_PixelFlags_flag_injected_templateCenter'])
726 else:
727 self._check_diaSource(tmplt_fake_sources, diaSrc, scale=-1, refIds=tmpl_refIds)
728 self.assertTrue(diaSrc['base_PixelFlags_flag_injected_template'])
729 self.assertTrue(diaSrc['base_PixelFlags_flag_injected_templateCenter'])
730 self.assertFalse(diaSrc['base_PixelFlags_flag_injected'])
731 self.assertFalse(diaSrc['base_PixelFlags_flag_injectedCenter'])
733 def test_mask_streaks(self):
734 """Run detection on a difference image containing a streak.
735 """
736 # Set up the simulated images
737 noiseLevel = 1.
738 staticSeed = 1
739 fluxLevel = 500
740 xSize = 400
741 ySize = 400
742 psfSize = 2.4
743 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel,
744 "xSize": xSize, "ySize": ySize}
745 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
746 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
748 # Configure the detection Task
749 detectionTask = self._setup_detection(doMerge=False, doMaskStreaks=True,
750 raiseOnNoDiaSources=False, raiseOnBadSubtractionRatio=False)
752 # Test that no streaks are detected
753 difference = science.clone()
754 difference.maskedImage -= matchedTemplate.maskedImage
755 output = detectionTask.run(science, matchedTemplate, difference, sources)
756 outMask = output.subtractedMeasuredExposure.mask.array
757 streakMask = output.subtractedMeasuredExposure.mask.getPlaneBitMask("STREAK")
758 streakMaskSet = (outMask & streakMask) > 0
759 self.assertTrue(np.all(streakMaskSet == 0))
761 # Add streak-like shape and check that streak is detected
762 difference.image.array[20:23, 40:200] += 50
763 output = detectionTask.run(science, matchedTemplate, difference, sources)
764 outMask = output.subtractedMeasuredExposure.mask.array
765 streakMask = output.subtractedMeasuredExposure.mask.getPlaneBitMask("STREAK")
766 streakMaskSet = (outMask & streakMask) > 0
767 self.assertTrue(np.all(streakMaskSet[20:23, 40:200]))
769 # Add a more trapezoid shaped streak across an image that is
770 # fainter and check that it is also detected
771 bbox = science.getBBox()
772 difference = science.clone()
773 difference.maskedImage -= matchedTemplate.maskedImage
774 width = 100
775 amplitude = 5
776 x0 = -100 + bbox.getBeginX()
777 y0 = -100 + bbox.getBeginY()
778 x1 = xSize + x0 + 100
779 y1 = ySize/2
780 corner_coords = [lsst.geom.Point2D(x0, y0),
781 lsst.geom.Point2D(x0, y0 + width),
782 lsst.geom.Point2D(x1, y1),
783 lsst.geom.Point2D(x1 + width, y1)]
784 streak_trapezoid = afwGeom.Polygon(corner_coords)
785 streak_image = streak_trapezoid.createImage(bbox)
786 streak_image *= amplitude
787 difference.image.array += streak_image.array
788 output = detectionTask.run(science, matchedTemplate, difference, sources)
789 outMask = output.subtractedMeasuredExposure.mask.array
790 streakMask = output.subtractedMeasuredExposure.mask.getPlaneBitMask("STREAK")
791 streakMaskSet = (outMask & streakMask) > 0
792 # Check that pixel values in streak_image equal the streak amplitude
793 streak_check = np.where(streak_image.array == amplitude)
794 self.assertTrue(np.all(streakMaskSet[streak_check]))
795 # Check that the entire image was not masked STREAK
796 self.assertFalse(np.all(streakMaskSet))
798 def _setup_sattle_tests(self):
799 noiseLevel = 1.
800 staticSeed = 1
801 fluxLevel = 500
802 psfSize = 2.4
803 shared_kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel,
804 "x0": 12345, "y0": 67890}
805 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6,
806 **shared_kwargs)
807 science.getInfo().setVisitInfo(makeVisitInfo(id=2))
808 detector = DetectorWrapper(numAmps=1).detector
809 science.setDetector(detector)
810 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel / 4,
811 noiseSeed=7, **shared_kwargs)
812 difference = science.clone()
814 detectionTask = self._setup_detection(doDeblend=True,
815 badSubtractionRatioThreshold=1.,
816 doSkySources=False, run_sattle=True)
818 return science, matchedTemplate, difference, sources, detectionTask
820 @mock.patch.dict(os.environ, {"SATTLE_URI_BASE": "fake_host:1234"})
821 def test_sattle_not_available(self):
822 science, matchedTemplate, difference, sources, detectionTask = self._setup_sattle_tests()
824 response = MockResponse({"allow_list": []}, 500, "sattle internal error")
825 with mock.patch('requests.put', return_value=response):
826 with self.assertRaises(requests.exceptions.HTTPError):
827 detectionTask.run(science, matchedTemplate, difference, sources,
828 idFactory=IdFactory.makeSimple())
830 @mock.patch.dict(os.environ, {"SATTLE_URI_BASE": "fake_host:1234"})
831 def test_visit_id_not_in_sattle(self):
832 science, matchedTemplate, difference, sources, detectionTask = self._setup_sattle_tests()
834 response = MockResponse({"allow_list": []}, 404, "missing visit cache")
835 # visit id not in sattle raises
836 with self.assertRaises(requests.exceptions.HTTPError):
837 with mock.patch('lsst.ip.diffim.detectAndMeasure.requests.put',
838 return_value=response):
839 with mock.patch('lsst.ip.diffim.utils.populate_sattle_visit_cache'):
840 detectionTask.run(science, matchedTemplate, difference, sources,
841 idFactory=IdFactory.makeSimple())
843 @mock.patch.dict(os.environ, {"SATTLE_URI_BASE": "fake_host:1234"})
844 def test_filter_satellites_some_allowed(self):
845 science, matchedTemplate, difference, sources, detectionTask = self._setup_sattle_tests()
847 allowed_ids = [1, 5]
848 response = MockResponse({"allow_list": allowed_ids}, 200, "some allowed")
849 with mock.patch('requests.put', return_value=response):
850 output = detectionTask.run(science, matchedTemplate, difference, sources,
851 idFactory=IdFactory.makeSimple())
853 self.assertEqual(len(output.diaSources), 2)
855 # Output should be sources 1 and 5 allowed out of 20
856 self.assertEqual(set(output.diaSources['id']), set(allowed_ids))
858 @mock.patch.dict(os.environ, {"SATTLE_URI_BASE": "fake_host:1234"})
859 def test_filter_satellites_all_allowed(self):
860 science, matchedTemplate, difference, sources, detectionTask = self._setup_sattle_tests()
862 allowed_ids = list(range(1, 21))
863 response = MockResponse({"allow_list": allowed_ids}, 200, "all allowed")
864 # Run detection and check the results
865 with mock.patch('requests.put', return_value=response):
866 output = detectionTask.run(science, matchedTemplate, difference, sources,
867 idFactory=IdFactory.makeSimple())
869 # Output should be all sources that went in. 20 go in, 20 should come out
870 self.assertEqual(len(output.diaSources), 20)
872 self.assertEqual(set(output.diaSources['id']), set(allowed_ids))
874 @mock.patch.dict(os.environ, {"SATTLE_URI_BASE": "fake_host:1234"})
875 def test_filter_satellites_none_allowed(self):
876 science, matchedTemplate, difference, sources, detectionTask = self._setup_sattle_tests()
878 response = MockResponse({"allow_list": []}, 200, "none allowed")
879 # Run detection and confirm it raises for no diasources
880 with self.assertRaises(detectAndMeasure.NoDiaSourcesError):
881 with mock.patch('requests.put', return_value=response):
882 detectionTask.run(science, matchedTemplate, difference, sources,
883 idFactory=IdFactory.makeSimple())
885 @mock.patch.dict(os.environ, {"SATTLE_URI_BASE": ""})
886 def test_fail_on_sattle_misconfiguration(self):
887 with self.assertRaises(pexConfig.FieldValidationError):
888 self._setup_detection(run_sattle=True)
890 def test_trailed_glints(self):
891 """Test that the glint_trail column works, and that
892 the trailed_glints output contains the expected information.
893 """
894 noiseLevel = 1.
895 staticSeed = 1
896 diffim, diaSources = makeTestImage(seed=staticSeed, noiseLevel=noiseLevel, noiseSeed=6)
897 self._check_values(diaSources['glint_trail'])
899 # Run detection and return the output Struct so we can check it
900 def _detection_wrapper(diffim, diaSources):
901 detectionTask = self._setup_detection()
902 scienceBase, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6)
903 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7)
904 science = scienceBase.clone()
905 science.maskedImage -= diffim.maskedImage
906 difference = science.clone()
907 difference.maskedImage -= matchedTemplate.maskedImage
908 output = detectionTask.run(science, matchedTemplate, difference, sources)
909 return output
911 output = _detection_wrapper(diffim, diaSources)
912 self.assertTrue('slopes' in output.glintTrailInfo)
913 self.assertTrue('intercepts' in output.glintTrailInfo)
914 self.assertTrue('stderrs' in output.glintTrailInfo)
915 self.assertTrue('lengths' in output.glintTrailInfo)
916 self.assertTrue('angles' in output.glintTrailInfo)
919class DetectAndMeasureScoreTest(DetectAndMeasureTestBase, lsst.utils.tests.TestCase):
920 detectionTask = detectAndMeasure.DetectAndMeasureScoreTask
922 def test_detection_xy0(self):
923 """Basic functionality test with non-zero x0 and y0.
924 """
925 # Set up the simulated images
926 noiseLevel = 1.
927 staticSeed = 1
928 fluxLevel = 500
929 kernelSize = 31
930 psfSize = 2.4
931 # Buffer source positions away from the score image's EDGE strip.
932 # Preconvolution adds ~kernelSize//2 pixels of EDGE on top of the
933 # EDGE already set by detection smoothing on the science image.
934 templateBorderSize = kernelSize//2
935 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel, "x0": 12345, "y0": 67890,
936 "kernelSize": kernelSize, "templateBorderSize": templateBorderSize}
937 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
938 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
939 difference = science.clone()
940 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
941 scienceKernel = science.psf.getKernel()
942 score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
944 # Configure the detection Task
945 # Set the subtraction residual metric threshold high, since the template
946 # is not actually subtracted from the science image.
947 detectionTask = self._setup_detection(doDeblend=True, badSubtractionRatioThreshold=1.,
948 doSkySources=False)
950 # Run detection and check the results
951 output = detectionTask.run(science, matchedTemplate, difference, score, sources,
952 idFactory=self.idGenerator.make_table_id_factory())
954 # Catalog ids should be very large from this id generator.
955 self.assertTrue(all(output.diaSources['id'] > 1000000000))
956 subtractedMeasuredExposure = output.subtractedMeasuredExposure
958 self.assertImagesEqual(subtractedMeasuredExposure.image, difference.image)
960 self.assertEqual(len(output.diaSources), len(sources))
961 # no sources should be flagged as negative
962 self.assertEqual(len(~output.diaSources["is_negative"]), len(output.diaSources))
963 refIds = []
964 for source in sources:
965 self._check_diaSource(output.diaSources, source, refIds=refIds)
967 def test_measurements_finite(self):
968 """Measured fluxes and centroids should always be finite.
969 """
970 columnNames = ["coord_ra", "coord_dec", "ip_diffim_forced_PsfFlux_instFlux",
971 "ip_diffim_forced_template_PsfFlux_instFlux"]
973 # Set up the simulated images
974 noiseLevel = 1.
975 staticSeed = 1
976 transientSeed = 6
977 xSize = 256
978 ySize = 256
979 psfSize = 2.4
980 kwargs = {"psfSize": psfSize, "x0": 0, "y0": 0,
981 "xSize": xSize, "ySize": ySize}
982 science, sources = makeTestImage(seed=staticSeed, noiseLevel=noiseLevel, noiseSeed=6,
983 nSrc=1, **kwargs)
984 matchedTemplate, _ = makeTestImage(seed=staticSeed, noiseLevel=noiseLevel/4, noiseSeed=7,
985 nSrc=1, **kwargs)
986 rng = np.random.RandomState(3)
987 xLoc = np.arange(-5, xSize+5, 10)
988 rng.shuffle(xLoc)
989 yLoc = np.arange(-5, ySize+5, 10)
990 rng.shuffle(yLoc)
991 transients, transientSources = makeTestImage(seed=transientSeed,
992 nSrc=len(xLoc), fluxLevel=1000.,
993 noiseLevel=noiseLevel, noiseSeed=8,
994 xLoc=xLoc, yLoc=yLoc,
995 **kwargs)
996 difference = science.clone()
997 difference.maskedImage -= matchedTemplate.maskedImage
998 difference.maskedImage += transients.maskedImage
999 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
1000 scienceKernel = science.psf.getKernel()
1001 score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
1003 # Configure the detection Task
1004 detectionTask = self._setup_detection(doForcedMeasurement=True)
1006 # Run detection and check the results
1007 output = detectionTask.run(science, matchedTemplate, difference, score, sources)
1009 for column in columnNames:
1010 self._check_values(output.diaSources[column])
1011 self._check_values(output.diaSources.getX(), minValue=0, maxValue=xSize)
1012 self._check_values(output.diaSources.getY(), minValue=0, maxValue=ySize)
1013 self._check_values(output.diaSources.getPsfInstFlux())
1015 def test_detect_transients(self):
1016 """Run detection on a difference image containing transients.
1017 """
1018 # Set up the simulated images
1019 noiseLevel = 1.
1020 staticSeed = 1
1021 transientSeed = 6
1022 nTransients = 10
1023 transientFlux = 1000
1024 fluxLevel = 500
1025 psfSize = 2.4
1026 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel}
1027 scienceBase, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
1028 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
1029 scienceKernel = scienceBase.psf.getKernel()
1030 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
1032 # Configure the detection Task
1033 detectionTask = self._setup_detection(doMerge=False)
1034 kwargs["seed"] = transientSeed
1035 kwargs["nSrc"] = nTransients
1036 kwargs["fluxLevel"] = transientFlux
1038 # Run detection and check the results
1039 def _run_and_check_detections(positive=True):
1040 """Simulate positive or negative transients and run detection.
1042 Parameters
1043 ----------
1044 positive : `bool`, optional
1045 If set, use positive transient sources.
1046 """
1048 transients, transientSources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, **kwargs)
1049 science = scienceBase.clone()
1050 if positive:
1051 science.maskedImage += transients.maskedImage
1052 else:
1053 science.maskedImage -= transients.maskedImage
1054 difference = science.clone()
1055 difference.maskedImage -= matchedTemplate.maskedImage
1056 score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
1057 # NOTE: NoiseReplacer (run by forcedMeasurement) can modify the
1058 # science image if we've e.g. removed parents post-deblending.
1059 # Pass a clone of the science image, so that it doesn't disrupt
1060 # later tests.
1061 output = detectionTask.run(science, matchedTemplate, difference, score, sources)
1062 refIds = []
1063 scale = 1. if positive else -1.
1064 # sources near the edge may have untrustworthy centroids
1065 goodSrcFlags = ~output.diaSources['base_PixelFlags_flag_edge']
1066 for diaSource, goodSrcFlag in zip(output.diaSources, goodSrcFlags):
1067 if goodSrcFlag:
1068 self._check_diaSource(transientSources, diaSource, refIds=refIds, scale=scale)
1069 _run_and_check_detections(positive=True)
1070 _run_and_check_detections(positive=False)
1072 def test_detect_dipoles(self):
1073 """Run detection on a difference image containing dipoles.
1074 """
1075 # Set up the simulated images
1076 noiseLevel = 1.
1077 staticSeed = 1
1078 fluxLevel = 1000
1079 fluxRange = 1.5
1080 nSources = 10
1081 offset = 1
1082 xSize = 300
1083 ySize = 300
1084 kernelSize = 31
1085 psfSize = 2.4
1086 # Avoid placing sources near the edge for this test, so that we can
1087 # easily check that the correct number of sources are detected.
1088 templateBorderSize = kernelSize//2
1089 dipoleFlag = "ip_diffim_DipoleFit_classification"
1090 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel, "fluxRange": fluxRange,
1091 "nSrc": nSources, "templateBorderSize": templateBorderSize, "kernelSize": kernelSize,
1092 "xSize": xSize, "ySize": ySize}
1093 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
1094 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
1095 difference = science.clone()
1096 # Shift the template by a pixel in order to make dipoles in the difference image.
1097 matchedTemplate.image.array[...] = np.roll(matchedTemplate.image.array[...], offset, axis=0)
1098 matchedTemplate.variance.array[...] = np.roll(matchedTemplate.variance.array[...], offset, axis=0)
1099 matchedTemplate.mask.array[...] = np.roll(matchedTemplate.mask.array[...], offset, axis=0)
1100 difference.maskedImage -= matchedTemplate.maskedImage[science.getBBox()]
1101 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
1102 scienceKernel = science.psf.getKernel()
1103 score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
1105 detectionTask = self._setup_detection(badSubtractionRatioThreshold=0.3)
1106 output = detectionTask.run(science, matchedTemplate, difference, score, sources)
1107 diaSources = output.diaSources[~output.diaSources["sky_source"]].copy(deep=True)
1108 self.assertEqual(len(diaSources), len(sources))
1109 refIds = []
1110 for diaSource in diaSources:
1111 if diaSource[dipoleFlag]: 1111 ↛ 1118line 1111 didn't jump to line 1118 because the condition on line 1111 was always true
1112 self._check_diaSource(sources, diaSource, refIds=refIds, scale=0,
1113 rtol=0.05, atol=None, usePsfFlux=False)
1114 self.assertFloatsAlmostEqual(diaSource["ip_diffim_DipoleFit_orientation"],
1115 -np.pi / 2, atol=2.)
1116 self.assertFloatsAlmostEqual(diaSource["ip_diffim_DipoleFit_separation"], offset, rtol=0.1)
1117 else:
1118 raise ValueError("DiaSource with ID %s is not a dipole!", diaSource.getId())
1120 def test_sky_sources(self):
1121 """Add sky sources and check that they are sufficiently far from other
1122 sources and have negligible flux.
1123 """
1124 # Set up the simulated images
1125 noiseLevel = 1.
1126 staticSeed = 1
1127 transientSeed = 6
1128 transientFluxLevel = 1000.
1129 transientFluxRange = 1.5
1130 fluxLevel = 500
1131 psfSize = 2.4
1132 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel}
1133 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
1134 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
1135 transients, transientSources = makeTestImage(seed=transientSeed, psfSize=2.4,
1136 nSrc=10, fluxLevel=transientFluxLevel,
1137 fluxRange=transientFluxRange,
1138 noiseLevel=noiseLevel, noiseSeed=8)
1139 difference = science.clone()
1140 difference.maskedImage -= matchedTemplate.maskedImage
1141 difference.maskedImage += transients.maskedImage
1142 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
1143 scienceKernel = science.psf.getKernel()
1144 kernelWidth = np.max(scienceKernel.getDimensions())//2
1145 score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
1147 # Configure the detection Task
1148 detectionTask = self._setup_detection(doSkySources=True)
1150 # Run detection and check the results
1151 output = detectionTask.run(science, matchedTemplate, difference, score, sources,
1152 idFactory=self.idGenerator.make_table_id_factory())
1153 nSkySourcesGenerated = detectionTask.metadata["n_skySources"]
1154 skySources = output.diaSources[output.diaSources["sky_source"]]
1155 self.assertEqual(len(skySources), nSkySourcesGenerated)
1156 for skySource in skySources:
1157 # Disable the sky_source flag to enable checks.
1158 skySource["sky_source"] = False
1159 # The sky sources should not be close to any other source
1160 with self.assertRaises(AssertionError):
1161 self._check_diaSource(transientSources, skySource, matchDistance=kernelWidth)
1162 with self.assertRaises(AssertionError):
1163 self._check_diaSource(sources, skySource, matchDistance=kernelWidth)
1164 # The sky sources should have low flux levels.
1165 self._check_diaSource(transientSources, skySource, matchDistance=1000, scale=0.,
1166 atol=np.sqrt(transientFluxRange*transientFluxLevel))
1168 # Catalog ids should be very large from this id generator.
1169 self.assertTrue(all(output.diaSources['id'] > 1000000000))
1171 def test_exclude_mask_detections(self):
1172 """Sources with certain bad mask planes set should not be detected.
1173 """
1174 # Set up the simulated images
1175 noiseLevel = 1.
1176 staticSeed = 1
1177 transientSeed = 6
1178 fluxLevel = 500
1179 radius = 2
1180 psfSize = 2.4
1181 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel}
1182 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
1183 science.getInfo().setVisitInfo(makeVisitInfo())
1184 detector = DetectorWrapper(numAmps=1).detector
1185 science.setDetector(detector)
1186 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
1188 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
1189 scienceKernel = science.psf.getKernel()
1190 # Configure the detection Task
1191 badSourceFlags = ["base_PixelFlags_flag_offimage",
1192 "base_PixelFlags_flag_edgeCenterAll",
1193 "base_PixelFlags_flag_badCenterAll",
1194 "base_PixelFlags_flag_saturatedCenterAll", ]
1195 detectionTask = self._setup_detection(nSkySources=0, badSourceFlags=badSourceFlags)
1196 excludeMaskPlanes = ["EDGE", "BAD", "SAT"]
1197 nBad = len(excludeMaskPlanes)
1198 self.assertGreater(nBad, 0)
1199 kwargs["seed"] = transientSeed
1200 kwargs["nSrc"] = nBad
1201 kwargs["fluxLevel"] = 1000
1203 # Run detection and check the results
1204 def _run_and_check_detections(setFlags=True):
1205 transients, transientSources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, **kwargs)
1206 difference = science.clone()
1207 difference.maskedImage -= matchedTemplate.maskedImage
1208 difference.maskedImage += transients.maskedImage
1209 if setFlags:
1210 for src, badMask in zip(transientSources, excludeMaskPlanes):
1211 srcX = int(src.getX())
1212 srcY = int(src.getY())
1213 srcBbox = lsst.geom.Box2I(lsst.geom.Point2I(srcX - radius, srcY - radius),
1214 lsst.geom.Extent2I(2*radius + 1, 2*radius + 1))
1215 difference[srcBbox].mask.array |= lsst.afw.image.Mask.getPlaneBitMask(badMask)
1216 score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
1217 if setFlags:
1218 with self.assertRaises(AlgorithmError):
1219 output = detectionTask.run(science, matchedTemplate, difference, score, sources)
1220 else:
1221 output = detectionTask.run(science, matchedTemplate, difference, score, sources)
1222 refIds = []
1223 goodSrcFlags = checkMask(difference.mask, transientSources, excludeMaskPlanes)
1224 self.assertEqual(np.sum(~goodSrcFlags), 0)
1225 for diaSource, goodSrcFlag in zip(output.diaSources, goodSrcFlags):
1226 if ~goodSrcFlag: 1226 ↛ 1227line 1226 didn't jump to line 1227 because the condition on line 1226 was never true
1227 with self.assertRaises(AssertionError):
1228 self._check_diaSource(transientSources, diaSource, refIds=refIds)
1229 else:
1230 self._check_diaSource(transientSources, diaSource, refIds=refIds)
1231 _run_and_check_detections(setFlags=False)
1232 _run_and_check_detections(setFlags=True)
1234 def test_detect_transients_with_background(self):
1235 """Background measured on the difference image should be subtracted
1236 from the score image so that transients are still detected cleanly.
1237 """
1238 # Set up the simulated images
1239 noiseLevel = 1.
1240 staticSeed = 1
1241 transientSeed = 6
1242 nTransients = 10
1243 transientFlux = 1000
1244 fluxLevel = 500
1245 xSize = 512
1246 ySize = 512
1247 x0 = 123
1248 y0 = 456
1249 kernelSize = 31
1250 psfSize = 2.4
1251 # Avoid placing sources near the edge for this test, so that we can
1252 # easily check that the correct number of sources are detected.
1253 templateBorderSize = kernelSize//2
1254 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel,
1255 "xSize": xSize, "ySize": ySize, "x0": x0, "y0": y0,
1256 "kernelSize": kernelSize, "templateBorderSize": templateBorderSize}
1257 chebyshev_params = [2.2, 2.1, 2.0, 1.2, 1.1, 1.0]
1259 # The background model must cover the grown image bbox, otherwise it
1260 # would be extrapolated outside its declared domain.
1261 bbox2D = lsst.geom.Box2D(
1262 lsst.geom.Point2D(x0 - templateBorderSize, y0 - templateBorderSize),
1263 lsst.geom.Extent2D(xSize + 2*templateBorderSize, ySize + 2*templateBorderSize))
1264 background_model = afwMath.Chebyshev1Function2D(chebyshev_params, bbox2D)
1265 scienceBase, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6,
1266 background=background_model, **kwargs)
1267 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
1268 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
1269 scienceKernel = scienceBase.psf.getKernel()
1271 # Configure the detection Task
1272 detectionTask = self._setup_detection(doMerge=False, doSubtractBackground=True)
1273 kwargs["seed"] = transientSeed
1274 kwargs["nSrc"] = nTransients
1275 kwargs["fluxLevel"] = transientFlux
1277 def _run_and_check_detections(positive=True):
1278 transients, transientSources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, **kwargs)
1279 science = scienceBase.clone()
1280 if positive:
1281 science.maskedImage += transients.maskedImage
1282 else:
1283 science.maskedImage -= transients.maskedImage
1284 difference = science.clone()
1285 difference.maskedImage -= matchedTemplate.maskedImage
1286 score = subtractTask._convolveExposure(difference, scienceKernel,
1287 subtractTask.convolutionControl)
1288 # Record the score image metric before detection so we can verify
1289 # that the background was actually subtracted from it.
1290 originalScoreMedian = np.median(score.image.array[np.isfinite(score.image.array)])
1291 originalDiffimMedian = np.median(difference.image.array[np.isfinite(difference.image.array)])
1292 output = detectionTask.run(science, matchedTemplate, difference, score, sources)
1294 # The background subtraction should leave a non-empty BackgroundList.
1295 self.assertGreater(len(output.differenceBackground), 0)
1296 # The difference image should have been modified in place by the
1297 # background subtraction, so its median should shift noticeably.
1298 subtractedScoreMedian = np.median(score.image.array[np.isfinite(score.image.array)])
1299 self.assertLess(np.abs(subtractedScoreMedian), np.abs(originalScoreMedian))
1300 subtractedDiffimMedian = np.median(difference.image.array[np.isfinite(difference.image.array)])
1301 self.assertLess(np.abs(subtractedDiffimMedian), np.abs(originalDiffimMedian))
1303 refIds = []
1304 scale = 1. if positive else -1.
1305 for transient in transientSources:
1306 self._check_diaSource(output.diaSources, transient, refIds=refIds, scale=scale)
1307 _run_and_check_detections(positive=True)
1308 _run_and_check_detections(positive=False)
1310 def test_mask_cosmic_rays(self):
1311 """Cosmic rays detected on the difference image should propagate
1312 to the mask of the returned (measured) exposure. This version creates
1313 and uses a Score image for detection.
1314 """
1315 # Set up the simulated images
1316 noiseLevel = 1.
1317 staticSeed = 1
1318 fluxLevel = 500
1319 xSize = 400
1320 ySize = 400
1321 psfSize = 2.4
1322 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel,
1323 "xSize": xSize, "ySize": ySize}
1324 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
1325 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
1326 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
1327 scienceKernel = science.psf.getKernel()
1328 crMask = science.mask.getPlaneBitMask("CR")
1330 # Configure the detection Task
1331 detectionTask = self._setup_detection(doMerge=False, doSkySources=True)
1333 # Test that no CRs are detected when none are present.
1334 transients, _ = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, nSrc=10, **kwargs)
1335 science.maskedImage += transients.maskedImage
1336 difference = science.clone()
1337 difference.maskedImage -= matchedTemplate.maskedImage
1338 score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
1339 output = detectionTask.run(science, matchedTemplate, difference, score, sources)
1340 self.assertFalse(np.any((output.subtractedMeasuredExposure.mask.array & crMask) > 0))
1342 crX0 = round(sources.getX()[0] - science.getX0())
1343 crY0 = round(sources.getY()[0] - science.getY0())
1344 crX1 = round(sources.getX()[1] - science.getX0())
1345 crY1 = round(sources.getY()[1] - science.getY0())
1346 # Inject CR-like shapes into the difference image. CR detection runs
1347 # on the difference, and the mask should propagate to the measured
1348 # exposure returned by the task.
1349 crMaskSetInput = np.zeros(difference.image.array.shape, bool)
1350 crMaskSetInput[crX0:crX0+1, crY0:crY0+5] = True
1351 crMaskSetInput[crX1:crX1+5, crY1:crY1+1] = True
1352 difference.image.array[crMaskSetInput] += 1234
1353 score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
1354 output = detectionTask.run(science, matchedTemplate, difference, score, sources)
1355 crMaskSet = (output.subtractedMeasuredExposure.mask.array & crMask) > 0
1356 self.assertFalse(np.any(crMaskSet[~crMaskSetInput]))
1357 self.assertTrue(np.all(crMaskSet[crMaskSetInput]))
1360class TestNegativePeaks(lsst.utils.tests.TestCase):
1361 """Tests of deblending and merging negative peaks, to test fixes for the
1362 various problems found on DM-48596.
1363 """
1365 def testDeblendNegatives(self):
1366 """Test that negative peaks get deblended and not destroyed: DM-48704.
1367 This is only a test of deblending, not of merging.
1368 """
1369 # Make a science image with one blend of two positive sources, to
1370 # subtract from an empty template, resulting in a negative diffim.
1371 bbox = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(100, 100))
1372 dataset = lsst.meas.base.tests.TestDataset(bbox)
1373 delta = 10
1374 with dataset.addBlend() as family:
1375 family.addChild(instFlux=2E5, centroid=lsst.geom.Point2D(50, 72))
1376 family.addChild(instFlux=2.5E5, centroid=lsst.geom.Point2D(50+delta, 74))
1377 science, catalog = dataset.realize(noise=1.0,
1378 schema=lsst.meas.base.tests.TestDataset.makeMinimalSchema())
1379 dataset = lsst.meas.base.tests.TestDataset(bbox)
1380 template, _ = dataset.realize(noise=1.0,
1381 schema=lsst.meas.base.tests.TestDataset.makeMinimalSchema())
1382 difference = template.clone()
1383 difference.image -= science.image
1385 config = detectAndMeasure.DetectAndMeasureTask.ConfigClass()
1386 config.doDeblend = True
1387 task = detectAndMeasure.DetectAndMeasureTask(config=config)
1388 # prelude steps taken from `detectAndMeasure.run`
1389 task._prepareInputs(difference)
1390 table = lsst.afw.table.SourceTable.make(task.schema)
1391 results = task.detection.run(table=table, exposure=difference, doSmooth=True)
1392 # Just run the deblend step so we can check the footprints independently.
1393 sources, positives, negatives = task._deblend(difference, results.positive, results.negative)
1395 # DM-48704 fixed a problem where the peaks were in the footprints, but
1396 # the spans were empty.
1397 self.assertEqual(len(negatives), 2)
1398 self.assertEqual(len(positives), 0)
1399 self.assertGreater(negatives[0].getFootprint().getSpans().getArea(), 0)
1400 self.assertGreater(negatives[1].getFootprint().getSpans().getArea(), 0)
1401 # Deblended children are HeavyFootprints; we have to make sure the
1402 # pixel values in those are correct (though DetectAndMeasureTask
1403 # doesn't use the fact that they're Heavy).
1404 # The sources are positive in the science image, and negative in the
1405 # diffim, so the minimum value in the deblended negative footprint is
1406 # the maximum value in the science catalog footprint (ignoring the
1407 # noise in the science and template, hence rtol).
1408 # (catalog[0] is the parent; we want the children)
1409 self.assertFloatsAlmostEqual(negatives[0].getFootprint().getImageArray().min(),
1410 -catalog[1].getFootprint().getImageArray().max(), rtol=1e-4)
1411 self.assertFloatsAlmostEqual(negatives[1].getFootprint().getImageArray().min(),
1412 -catalog[2].getFootprint().getImageArray().max(), rtol=1e-4)
1414 self.assertEqual(sources["is_negative"].sum(), 2)
1416 def testMergeFootprints(self):
1417 """Test that merging footprints does not cause negative ones to
1418 disappear (e.g. get merged into non-connected footprints).
1420 As implemented, the diffim will have three positive sources (one a
1421 blended pair), and 7 negative sources (two a blended pair).
1422 """
1423 # Make a template image with multiple blends, designed to trigger the
1424 # negative-footprint-related bug in `FootprintSet.merge`.
1425 bbox = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(400, 200))
1426 dataset = lsst.meas.base.tests.TestDataset(bbox)
1427 # Because these sources are on the template, they will be negative on
1428 # the diffim, unless the pixels are explicitly made negative via
1429 # `template.image.subset` below.
1430 # positive isolated source on diffim
1431 dataset.addSource(instFlux=.5E5, centroid=lsst.geom.Point2D(25, 26))
1432 # negative isolated source on diffim
1433 dataset.addSource(instFlux=.7E5, centroid=lsst.geom.Point2D(75, 24),
1434 shape=lsst.afw.geom.Quadrupole(12, 7, 2))
1435 delta = 10
1436 # negative blended pair on diffim
1437 with dataset.addBlend() as family:
1438 family.addChild(instFlux=1E5, centroid=lsst.geom.Point2D(150, 72))
1439 family.addChild(instFlux=1.5E5, centroid=lsst.geom.Point2D(150+delta, 74))
1440 # positive blended pair on diffim
1441 with dataset.addBlend() as family:
1442 family.addChild(instFlux=2E5, centroid=lsst.geom.Point2D(250, 72))
1443 family.addChild(instFlux=2.5E5, centroid=lsst.geom.Point2D(250+delta, 74))
1444 # negative blended pair on diffim
1445 with dataset.addBlend() as family:
1446 family.addChild(instFlux=3E5, centroid=lsst.geom.Point2D(350, 72))
1447 family.addChild(instFlux=3.5E5, centroid=lsst.geom.Point2D(350+delta, 74))
1448 # negative isolated source on diffim
1449 dataset.addSource(instFlux=4E5, centroid=lsst.geom.Point2D(75, 124))
1450 # negative isolated source on diffim
1451 dataset.addSource(instFlux=5E5, centroid=lsst.geom.Point2D(175, 124))
1452 schema = lsst.meas.base.tests.TestDataset.makeMinimalSchema()
1453 schema.addField("sky_source", "Flag", "Sky source.")
1454 template, catalog = dataset.realize(noise=100.0, schema=schema)
1455 # These will be positive sources on the diffim (as noted above).
1456 template.image.subset(lsst.geom.Box2I(lsst.geom.Point2I(15, 15),
1457 lsst.geom.Point2I(35, 35))).array *= -1
1458 template.image.subset(lsst.geom.Box2I(lsst.geom.Point2I(230, 60),
1459 lsst.geom.Point2I(270, 85))).array *= -1
1460 dataset = lsst.meas.base.tests.TestDataset(bbox)
1461 science, _ = dataset.realize(noise=100.0, schema=schema)
1462 difference = science.clone()
1463 difference.image -= template.image
1465 config = detectAndMeasure.DetectAndMeasureTask.ConfigClass()
1466 config.doDeblend = True
1467 config.raiseOnBadSubtractionRatio = False
1468 config.raiseOnNoDiaSources = False
1469 task = detectAndMeasure.DetectAndMeasureTask(config=config)
1470 result = task.run(science, template, difference, catalog)
1472 # The original bug manifested as the (175,124) source disappaering in
1473 # the merge step, and being included in the peaks of a duplicated
1474 # (75,124) source, even though their footprints are not connected.
1475 mask = np.isclose(result.diaSources["slot_Centroid_x"], 75) & \
1476 np.isclose(result.diaSources["slot_Centroid_y"], 124)
1477 self.assertEqual(mask.sum(), 1)
1478 peaks = result.diaSources[mask][0].getFootprint().peaks
1479 self.assertEqual(len(peaks), 1)
1480 self.assertEqual(peaks[0].getIx(), 75)
1481 self.assertEqual(peaks[0].getIy(), 124)
1483 mask = np.isclose(result.diaSources["slot_Centroid_x"], 175) & \
1484 np.isclose(result.diaSources["slot_Centroid_y"], 124)
1485 self.assertEqual(mask.sum(), 1)
1486 peaks = result.diaSources[mask][0].getFootprint().peaks
1487 self.assertEqual(len(peaks), 1)
1488 self.assertEqual(peaks[0].getIx(), 175)
1489 self.assertEqual(peaks[0].getIy(), 124)
1491 # This checks that all the returned footprints have exactly one peak;
1492 # it's a more general test than the above, but I'm keeping that for
1493 # the more explicit test of the original bugged sources.
1494 peak_x = np.array([peak['f_x'] for src in result.diaSources for peak in src.getFootprint().peaks])
1495 peak_y = np.array([peak['f_y'] for src in result.diaSources for peak in src.getFootprint().peaks])
1496 peaks = np.column_stack((peak_x, peak_y))
1497 unique_peaks, counts = np.unique(peaks, axis=0, return_counts=True)
1498 self.assertEqual(np.sum(counts > 1), 0)
1500 # There are three positive sources that should not have `is_negative`
1501 # set, independent of how deblending/merging of negative sources is
1502 # handled.
1503 self.assertEqual((~result.diaSources["is_negative"]).sum(), 3)
1505 def testReorderPeaksBySignificance(self):
1506 """A merged footprint whose most significant peak is negative should
1507 have that peak first after re-ordering, even though the merge lists
1508 positive peaks before negative ones.
1509 """
1510 config = detectAndMeasure.DetectAndMeasureTask.ConfigClass()
1511 task = detectAndMeasure.DetectAndMeasureTask(config=config)
1513 # Peak schema with a significance field, as SourceDetectionTask
1514 # produces for pixel_stdev thresholds.
1515 mapper = afwTable.SchemaMapper(afwDetection.PeakTable.makeMinimalSchema())
1516 mapper.addMinimalSchema(afwDetection.PeakTable.makeMinimalSchema())
1517 mapper.addOutputField("significance", type=float, doc="")
1518 peakSchema = mapper.getOutputSchema()
1520 diaSources = afwTable.SourceCatalog(task.schema)
1521 source = diaSources.addNew()
1522 spans = afwGeom.SpanSet(geom.Box2I(geom.Point2I(0, 0), geom.Extent2I(10, 10)))
1523 footprint = afwDetection.Footprint(spans, peakSchema)
1524 # Merge order is [positive, negative]: a weak positive peak first,
1525 # then a more significant negative peak.
1526 footprint.addPeak(1, 1, 30.0)
1527 footprint.peaks[-1]["significance"] = 5.0
1528 footprint.addPeak(2, 2, -500.0)
1529 footprint.peaks[-1]["significance"] = 40.0
1530 source.setFootprint(footprint)
1532 task._reorderPeaksBySignificance(diaSources)
1534 peaks = source.getFootprint().peaks
1535 self.assertEqual(peaks[0].getPeakValue(), -500.0)
1536 self.assertEqual(peaks[0]["significance"], 40.0)
1537 self.assertEqual(peaks[1].getPeakValue(), 30.0)
1538 self.assertEqual(peaks[1]["significance"], 5.0)
1541def makeVisitInfo(id=1):
1542 """Return a non-NaN visitInfo."""
1543 return afwImage.VisitInfo(id=id,
1544 exposureTime=10.01,
1545 darkTime=11.02,
1546 date=dafBase.DateTime(65321.1, dafBase.DateTime.MJD, dafBase.DateTime.TAI),
1547 ut1=12345.1,
1548 era=45.1*geom.degrees,
1549 boresightRaDec=geom.SpherePoint(23.1, 73.2, geom.degrees),
1550 boresightAzAlt=geom.SpherePoint(134.5, 33.3, geom.degrees),
1551 boresightAirmass=1.73,
1552 boresightRotAngle=73.2*geom.degrees,
1553 rotType=afwImage.RotType.SKY,
1554 observatory=Observatory(
1555 11.1*geom.degrees, 22.2*geom.degrees, 0.333),
1556 weather=Weather(1.1, 2.2, 34.5),
1557 )
1560class MockResponse:
1561 """Provide a mock for requests.put calls"""
1562 def __init__(self, json_data, status_code, text):
1563 self.json_data = json_data
1564 self.status_code = status_code
1565 self.text = text
1567 def json(self):
1568 return self.json_data
1570 def raise_for_status(self):
1571 if self.status_code != 200:
1572 raise requests.exceptions.HTTPError
1575def setup_module(module):
1576 lsst.utils.tests.init()
1579class MemoryTestCase(lsst.utils.tests.MemoryTestCase):
1580 pass
1583if __name__ == "__main__": 1583 ↛ 1584line 1583 didn't jump to line 1584 because the condition on line 1583 was never true
1584 lsst.utils.tests.init()
1585 unittest.main()