Coverage for tests / test_detectAndMeasure.py: 8%
944 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-05-07 08:35 +0000
1# This file is part of ip_diffim.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
22import numpy as np
23import os
24import unittest
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 psfSize = 2.4
172 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel, "x0": 12345, "y0": 67890}
173 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
174 science.getInfo().setVisitInfo(makeVisitInfo())
175 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
176 difference = science.clone()
178 # Configure the detection Task
179 # Set the subtraction residual metric threshold high, since the template
180 # is not actually subtracted from the science image.
181 detectionTask = self._setup_detection(doDeblend=True, badSubtractionRatioThreshold=1.,
182 doSkySources=False)
184 # Run detection and check the results
185 output = detectionTask.run(science, matchedTemplate, difference, sources,
186 idFactory=self.idGenerator.make_table_id_factory())
187 subtractedMeasuredExposure = output.subtractedMeasuredExposure
189 # Catalog ids should be very large from this id generator.
190 self.assertTrue(all(output.diaSources['id'] > 1000000000))
191 self.assertImagesEqual(subtractedMeasuredExposure.image, difference.image)
193 # all of the sources should have been detected
194 self.assertEqual(len(output.diaSources), len(sources))
195 # no sources should be flagged as negative
196 self.assertEqual(len(~output.diaSources["is_negative"]), len(output.diaSources))
197 refIds = []
198 for source in sources:
199 self._check_diaSource(output.diaSources, source, refIds=refIds)
201 def test_raise_bad_psf(self):
202 """Detection should raise if the PSF width is NaN
203 """
204 # Set up the simulated images
205 noiseLevel = 1.
206 staticSeed = 1
207 fluxLevel = 500
208 psfSize = 2.4
209 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel, "x0": 12345, "y0": 67890}
210 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
211 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
212 difference = science.clone()
214 psfShape = difference.getPsf().computeBBox(lsst.geom.Point2D(0, 0)).getDimensions()
215 psfcI = afwImage.ImageD(lsst.geom.Extent2I(psfShape[1], psfShape[0]))
216 psfNew = np.zeros(psfShape)
217 psfNew[0, :] = 1
218 psfcI.array = psfNew
219 psfcK = afwMath.FixedKernel(psfcI)
220 psf = measAlg.KernelPsf(psfcK)
221 difference.setPsf(psf)
223 # Configure the detection Task
224 detectionTask = self._setup_detection()
226 # Run detection and check the results
227 with self.assertRaises(UpstreamFailureNoWorkFound):
228 detectionTask.run(science, matchedTemplate, difference, sources)
230 def test_measurements_finite(self):
231 """Measured fluxes and centroids should always be finite.
232 """
233 columnNames = ["coord_ra", "coord_dec", "ip_diffim_forced_PsfFlux_instFlux",
234 "ip_diffim_forced_template_PsfFlux_instFlux"]
236 # Set up the simulated images
237 noiseLevel = 1.
238 staticSeed = 1
239 transientSeed = 6
240 xSize = 256
241 ySize = 256
242 psfSize = 2.4
243 kwargs = {"psfSize": psfSize, "x0": 0, "y0": 0,
244 "xSize": xSize, "ySize": ySize}
245 science, sources = makeTestImage(seed=staticSeed, noiseLevel=noiseLevel, noiseSeed=6,
246 nSrc=1, **kwargs)
247 science.getInfo().setVisitInfo(makeVisitInfo())
248 matchedTemplate, _ = makeTestImage(seed=staticSeed, noiseLevel=noiseLevel/4, noiseSeed=7,
249 nSrc=1, **kwargs)
250 rng = np.random.RandomState(3)
251 xLoc = np.arange(-5, xSize+5, 10)
252 rng.shuffle(xLoc)
253 yLoc = np.arange(-5, ySize+5, 10)
254 rng.shuffle(yLoc)
255 transients, transientSources = makeTestImage(seed=transientSeed,
256 nSrc=len(xLoc), fluxLevel=1000.,
257 noiseLevel=noiseLevel, noiseSeed=8,
258 xLoc=xLoc, yLoc=yLoc,
259 **kwargs)
260 difference = science.clone()
261 difference.maskedImage -= matchedTemplate.maskedImage
262 difference.maskedImage += transients.maskedImage
264 # Configure the detection Task
265 detectionTask = self._setup_detection(doForcedMeasurement=True)
267 # Run detection and check the results
268 output = detectionTask.run(science, matchedTemplate, difference, sources)
270 for column in columnNames:
271 self._check_values(output.diaSources[column])
272 self._check_values(output.diaSources.getX(), minValue=0, maxValue=xSize)
273 self._check_values(output.diaSources.getY(), minValue=0, maxValue=ySize)
274 self._check_values(output.diaSources.getPsfInstFlux())
276 def test_raise_config_schema_mismatch(self):
277 """Check that sources with specified flags are removed from the catalog.
278 """
279 # Configure the detection Task, and and set a config that is not in the schema
280 with self.assertRaises(InvalidQuantumError):
281 self._setup_detection(badSourceFlags=["Bogus_flag_42"])
283 def test_remove_unphysical(self):
284 """Check that sources with specified flags are removed from the catalog.
285 """
286 # Set up the simulated images
287 noiseLevel = 1.
288 staticSeed = 1
289 xSize = 256
290 ySize = 256
291 psfSize = 2.4
292 kwargs = {"psfSize": psfSize, "xSize": xSize, "ySize": ySize}
293 science, sources = makeTestImage(seed=staticSeed, noiseLevel=noiseLevel, noiseSeed=6,
294 nSrc=1, **kwargs)
295 science.getInfo().setVisitInfo(makeVisitInfo())
296 matchedTemplate, _ = makeTestImage(seed=staticSeed, noiseLevel=noiseLevel/4, noiseSeed=7,
297 nSrc=1, **kwargs)
298 transients, transientSources = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=8, nSrc=1, **kwargs)
299 science.maskedImage += transients.maskedImage
300 difference = science.clone()
301 bbox = difference.getBBox()
302 difference.maskedImage -= matchedTemplate.maskedImage
304 # Configure the detection Task, and remove unphysical sources
305 detectionTask = self._setup_detection(doForcedMeasurement=False, doSkySources=True, nSkySources=20,
306 badSourceFlags=["base_PixelFlags_flag_offimage", ])
308 # Run detection and check the results
309 diaSources = detectionTask.run(science, matchedTemplate, difference, sources).diaSources
310 badDiaSrcDoRemove = ~bbox.contains(diaSources.getX(), diaSources.getY())
311 nBadDoRemove = np.count_nonzero(badDiaSrcDoRemove)
312 # Verify that all sources are physical
313 self.assertEqual(nBadDoRemove, 0)
314 # Set a few centroids outside the image bounding box
315 nSetBad = 5
316 for src in diaSources[0: nSetBad]:
317 src["slot_Centroid_x"] += xSize
318 src["slot_Centroid_y"] += ySize
319 src["base_PixelFlags_flag_offimage"] = True
320 # Verify that these sources are outside the image
321 badDiaSrc = ~bbox.contains(diaSources.getX(), diaSources.getY())
322 nBad = np.count_nonzero(badDiaSrc)
323 self.assertEqual(nBad, nSetBad)
324 diaSourcesNoBad = detectionTask._removeBadSources(diaSources)
325 badDiaSrcNoBad = ~bbox.contains(diaSourcesNoBad.getX(), diaSourcesNoBad.getY())
327 # Verify that no sources outside the image bounding box remain
328 self.assertEqual(np.count_nonzero(badDiaSrcNoBad), 0)
329 self.assertEqual(len(diaSourcesNoBad), len(diaSources) - nSetBad)
331 def test_detect_transients(self):
332 """Run detection on a difference image containing transients.
333 """
334 # Set up the simulated images
335 noiseLevel = 1.
336 staticSeed = 1
337 transientSeed = 6
338 nTransients = 10
339 transientFlux = 1000
340 fluxLevel = 500
341 psfSize = 2.4
342 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel}
343 scienceBase, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
344 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
346 # Configure the detection Task
347 detectionTask = self._setup_detection(doMerge=False, doSkySources=True)
348 kwargs["seed"] = transientSeed
349 kwargs["nSrc"] = nTransients
350 kwargs["fluxLevel"] = transientFlux
352 # Run detection and check the results
353 def _run_and_check_detections(positive=True):
354 transients, transientSources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, **kwargs)
355 science = scienceBase.clone()
356 if positive:
357 science.maskedImage += transients.maskedImage
358 else:
359 science.maskedImage -= transients.maskedImage
360 difference = science.clone()
361 difference.maskedImage -= matchedTemplate.maskedImage
362 # NOTE: NoiseReplacer (run by forcedMeasurement) can modify the
363 # science image if we've e.g. removed parents post-deblending.
364 # Pass a clone of the science image, so that it doesn't disrupt
365 # later tests.
366 output = detectionTask.run(science, matchedTemplate, difference, sources)
367 refIds = []
368 scale = 1. if positive else -1.
369 for diaSource in output.diaSources:
370 self._check_diaSource(transientSources, diaSource, refIds=refIds, scale=scale)
371 _run_and_check_detections(positive=True)
372 _run_and_check_detections(positive=False)
374 def test_detect_transients_with_background(self):
375 """Run detection on a difference image containing transients and a background.
376 """
377 # Set up the simulated images
378 noiseLevel = 1.
379 staticSeed = 1
380 transientSeed = 6
381 nTransients = 10
382 transientFlux = 1000
383 fluxLevel = 500
384 xSize = 512
385 ySize = 512
386 x0 = 123
387 y0 = 456
388 psfSize = 2.4
389 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel,
390 "xSize": xSize, "ySize": ySize, "x0": x0, "y0": y0}
391 params = [2.2, 2.1, 2.0, 1.2, 1.1, 1.0]
393 bbox2D = lsst.geom.Box2D(lsst.geom.Point2D(x0, y0), lsst.geom.Extent2D(xSize, ySize))
394 background_model = afwMath.Chebyshev1Function2D(params, bbox2D)
395 scienceBase, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6,
396 background=background_model, **kwargs)
397 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
399 # Configure the detection Task
400 detectionTask = self._setup_detection(doMerge=False, doSubtractBackground=True)
401 kwargs["seed"] = transientSeed
402 kwargs["nSrc"] = nTransients
403 kwargs["fluxLevel"] = transientFlux
405 # Run detection and check the results
406 def _run_and_check_detections(positive=True):
407 transients, transientSources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, **kwargs)
408 science = scienceBase.clone()
409 if positive:
410 science.maskedImage += transients.maskedImage
411 else:
412 science.maskedImage -= transients.maskedImage
413 difference = science.clone()
414 difference.maskedImage -= matchedTemplate.maskedImage
415 # NOTE: NoiseReplacer (run by forcedMeasurement) can modify the
416 # science image if we've e.g. removed parents post-deblending.
417 # Pass a clone of the science image, so that it doesn't disrupt
418 # later tests.
419 output = detectionTask.run(science, matchedTemplate, difference, sources)
420 refIds = []
421 scale = 1. if positive else -1.
422 for transient in transientSources:
423 self._check_diaSource(output.diaSources, transient, refIds=refIds, scale=scale)
424 _run_and_check_detections(positive=True)
425 _run_and_check_detections(positive=False)
427 def test_mask_cosmic_rays(self):
428 """Run detection on a difference image containing a cosmic ray.
429 """
430 # Set up the simulated images
431 noiseLevel = 1.
432 staticSeed = 1
433 fluxLevel = 500
434 xSize = 400
435 ySize = 400
436 psfSize = 2.4
437 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel,
438 "xSize": xSize, "ySize": ySize}
439 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
440 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
441 crMask = science.mask.getPlaneBitMask("CR")
443 # Configure the detection Task
444 detectionTask = self._setup_detection(doMerge=False, doSkySources=True)
446 # Test that no CRs are detected
447 transients, transientSources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, nSrc=10, **kwargs)
448 science.maskedImage += transients.maskedImage
449 difference = science.clone()
450 difference.maskedImage -= matchedTemplate.maskedImage
451 output = detectionTask.run(science, matchedTemplate, difference, sources)
452 crMaskSet = (output.subtractedMeasuredExposure.mask.array & crMask) > 0
453 self.assertTrue(np.all(crMaskSet == 0))
455 crX0 = round(sources.getX()[0] - science.getX0())
456 crY0 = round(sources.getY()[0] - science.getY0())
457 crX1 = round(sources.getX()[1] - science.getX0())
458 crY1 = round(sources.getY()[1] - science.getY0())
459 # Add CR-like shape and check that CR is detected
460 # Pick two locations on top of sources, since that is what is likely to
461 # be missed in the first stage of CR rejection.
462 crMaskSetInput = np.zeros(difference.image.array.shape, bool)
463 crMaskSetInput[crX0:crX0+1, crY0:crY0+5] = True
464 crMaskSetInput[crX1:crX1+5, crY1:crY1+1] = True
465 difference.image.array[crMaskSetInput] += 1234
466 output = detectionTask.run(science, matchedTemplate, difference)
467 crMaskSet = (output.subtractedMeasuredExposure.mask.array & crMask) > 0
468 self.assertFalse(np.any(crMaskSet[~crMaskSetInput]))
469 self.assertTrue(np.all(crMaskSet[crMaskSetInput]))
471 def test_missing_mask_planes(self):
472 """Check that detection runs with missing mask planes.
473 """
474 # Set up the simulated images
475 noiseLevel = 1.
476 fluxLevel = 500
477 psfSize = 2.4
478 kwargs = {"psfSize": psfSize, "fluxLevel": fluxLevel, "addMaskPlanes": []}
479 # Use different seeds for the science and template so every source is a diaSource
480 science, sources = makeTestImage(seed=5, noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
481 science.getInfo().setVisitInfo(makeVisitInfo())
482 matchedTemplate, _ = makeTestImage(seed=6, noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
484 difference = science.clone()
485 difference.maskedImage -= matchedTemplate.maskedImage
486 detectionTask = self._setup_detection(raiseOnBadSubtractionRatio=False, raiseOnNoDiaSources=False)
488 # Verify that detection runs without errors
489 detectionTask.run(science, matchedTemplate, difference, sources)
491 def test_detect_dipoles(self):
492 """Run detection on a difference image containing dipoles.
493 """
494 # Set up the simulated images
495 noiseLevel = 1.
496 staticSeed = 1
497 fluxLevel = 1000
498 fluxRange = 1.5
499 nSources = 10
500 offset = 1
501 xSize = 300
502 ySize = 300
503 kernelSize = 31
504 psfSize = 2.4
505 # Avoid placing sources near the edge for this test, so that we can
506 # easily check that the correct number of sources are detected.
507 templateBorderSize = kernelSize//2
508 dipoleFlag = "ip_diffim_DipoleFit_classification"
509 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel, "fluxRange": fluxRange,
510 "nSrc": nSources, "templateBorderSize": templateBorderSize, "kernelSize": kernelSize,
511 "xSize": xSize, "ySize": ySize}
512 dipoleFlag = "ip_diffim_DipoleFit_classification"
513 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
514 science.getInfo().setVisitInfo(makeVisitInfo())
515 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
516 difference = science.clone()
517 matchedTemplate.image.array[...] = np.roll(matchedTemplate.image.array[...], offset, axis=0)
518 matchedTemplate.variance.array[...] = np.roll(matchedTemplate.variance.array[...], offset, axis=0)
519 matchedTemplate.mask.array[...] = np.roll(matchedTemplate.mask.array[...], offset, axis=0)
520 difference.maskedImage -= matchedTemplate.maskedImage[science.getBBox()]
522 # Need a higher residual metric threshold, since we are purposely
523 # creating poor subtractions.
524 detectionTask = self._setup_detection(doMerge=True, badSubtractionRatioThreshold=0.3)
525 output = detectionTask.run(science, matchedTemplate, difference, sources)
526 diaSources = output.diaSources[~output.diaSources["sky_source"]].copy(deep=True)
527 self.assertEqual(len(diaSources), len(sources))
528 # no sources should be flagged as negative
529 self.assertEqual(len(~diaSources["is_negative"]), len(diaSources))
530 refIds = []
531 for diaSource in diaSources:
532 if diaSource[dipoleFlag]:
533 self._check_diaSource(sources, diaSource, refIds=refIds, scale=0,
534 rtol=0.05, atol=None, usePsfFlux=False)
535 self.assertFloatsAlmostEqual(diaSource["ip_diffim_DipoleFit_orientation"],
536 -np.pi / 2, atol=2.)
537 self.assertFloatsAlmostEqual(diaSource["ip_diffim_DipoleFit_separation"], offset, rtol=0.1)
538 else:
539 raise ValueError("DiaSource with ID %s is not a dipole!", diaSource.getId())
541 def test_sky_sources(self):
542 """Add sky sources and check that they are sufficiently far from other
543 sources and have negligible flux.
544 """
545 # Set up the simulated images
546 noiseLevel = 1.
547 staticSeed = 1
548 transientSeed = 6
549 transientFluxLevel = 1000.
550 transientFluxRange = 1.5
551 fluxLevel = 500
552 psfSize = 2.4
553 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel}
554 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
555 science.getInfo().setVisitInfo(makeVisitInfo())
556 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
557 transients, transientSources = makeTestImage(seed=transientSeed, psfSize=2.4,
558 nSrc=10, fluxLevel=transientFluxLevel,
559 fluxRange=transientFluxRange,
560 noiseLevel=noiseLevel, noiseSeed=8)
561 difference = science.clone()
562 difference.maskedImage -= matchedTemplate.maskedImage
563 difference.maskedImage += transients.maskedImage
564 kernelWidth = np.max(science.psf.getKernel().getDimensions())//2
566 # Configure the detection Task
567 detectionTask = self._setup_detection(doSkySources=True)
569 # Run detection and check the results
570 output = detectionTask.run(science, matchedTemplate, difference, sources,
571 idFactory=self.idGenerator.make_table_id_factory())
572 skySources = output.diaSources[output.diaSources["sky_source"]]
573 self.assertEqual(len(skySources), detectionTask.config.skySources.nSources)
574 for skySource in skySources:
575 # Disable the sky_source flag to enable checks.
576 skySource["sky_source"] = False
577 # The sky sources should not be close to any other source
578 with self.assertRaises(AssertionError):
579 self._check_diaSource(transientSources, skySource, matchDistance=kernelWidth)
580 with self.assertRaises(AssertionError):
581 self._check_diaSource(sources, skySource, matchDistance=kernelWidth)
582 # The sky sources should have low flux levels.
583 self._check_diaSource(transientSources, skySource, matchDistance=1000, scale=0.,
584 atol=np.sqrt(transientFluxRange*transientFluxLevel))
586 # Catalog ids should be very large from this id generator.
587 self.assertTrue(all(output.diaSources['id'] > 1000000000))
589 def test_exclude_mask_detections(self):
590 """Sources with certain bad mask planes set should not be detected.
591 """
592 # Set up the simulated images
593 noiseLevel = 1.
594 staticSeed = 1
595 transientSeed = 6
596 fluxLevel = 500
597 radius = 2
598 psfSize = 2.4
599 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel}
600 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
601 science.getInfo().setVisitInfo(makeVisitInfo())
602 detector = DetectorWrapper(numAmps=1).detector
603 science.setDetector(detector)
604 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
606 # Configure the detection Task
607 badSourceFlags = ["base_PixelFlags_flag_offimage",
608 "base_PixelFlags_flag_edgeCenterAll",
609 "base_PixelFlags_flag_badCenterAll",
610 "base_PixelFlags_flag_saturatedCenterAll", ]
611 detectionTask = self._setup_detection(nSkySources=0, badSourceFlags=badSourceFlags)
612 excludeMaskPlanes = ["EDGE", "BAD", "SAT"]
613 nBad = len(excludeMaskPlanes)
614 self.assertGreater(nBad, 0)
615 kwargs["seed"] = transientSeed
616 kwargs["nSrc"] = nBad
617 kwargs["fluxLevel"] = 1000
619 # Run detection and check the results
620 def _run_and_check_detections(setFlags=True):
621 transients, transientSources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, **kwargs)
622 difference = science.clone()
623 difference.maskedImage -= matchedTemplate.maskedImage
624 difference.maskedImage += transients.maskedImage
625 if setFlags:
626 for src, badMask in zip(transientSources, excludeMaskPlanes):
627 srcX = int(src.getX())
628 srcY = int(src.getY())
629 srcBbox = lsst.geom.Box2I(lsst.geom.Point2I(srcX - radius, srcY - radius),
630 lsst.geom.Extent2I(2*radius + 1, 2*radius + 1))
631 difference[srcBbox].mask.array |= lsst.afw.image.Mask.getPlaneBitMask(badMask)
632 with self.assertRaises(AlgorithmError):
633 output = detectionTask.run(science, matchedTemplate, difference, sources)
634 else:
635 output = detectionTask.run(science, matchedTemplate, difference, sources)
636 refIds = []
637 goodSrcFlags = checkMask(difference.mask, transientSources, excludeMaskPlanes)
638 self.assertEqual(np.sum(~goodSrcFlags), 0)
639 for diaSource, goodSrcFlag in zip(output.diaSources, goodSrcFlags):
640 if ~goodSrcFlag:
641 with self.assertRaises(AssertionError):
642 self._check_diaSource(transientSources, diaSource, refIds=refIds)
643 else:
644 self._check_diaSource(transientSources, diaSource, refIds=refIds)
645 _run_and_check_detections(setFlags=False)
646 _run_and_check_detections(setFlags=True)
648 def test_fake_mask_plane_propagation(self):
649 """Test that we have the mask planes related to fakes in diffim images.
650 This is testing method called updateMasks
651 """
652 xSize = 256
653 ySize = 256
654 science, sources = makeTestImage(psfSize=2.4, xSize=xSize, ySize=ySize, doApplyCalibration=True)
655 science_fake_img, science_fake_sources = makeTestImage(
656 psfSize=2.4, xSize=xSize, ySize=ySize, seed=5, nSrc=3, noiseLevel=0.25, fluxRange=1
657 )
658 template, _ = makeTestImage(psfSize=2.4, xSize=xSize, ySize=ySize, doApplyCalibration=True)
659 tmplt_fake_img, tmplt_fake_sources = makeTestImage(
660 psfSize=2.4, xSize=xSize, ySize=ySize, seed=9, nSrc=3, noiseLevel=0.25, fluxRange=1
661 )
662 # created fakes and added them to the images
663 science.image += science_fake_img.image
664 template.image += tmplt_fake_img.image
666 # TODO: DM-40796 update to INJECTED names when source injection gets refactored
667 # adding mask planes to both science and template images
668 science.mask.addMaskPlane("FAKE")
669 science_fake_bitmask = science.mask.getPlaneBitMask("FAKE")
670 template.mask.addMaskPlane("FAKE")
671 template_fake_bitmask = template.mask.getPlaneBitMask("FAKE")
673 # makeTestImage sets the DETECTED plane on the sources; we can use
674 # that to set the FAKE plane on the science and template images.
675 detected = science_fake_img.mask.getPlaneBitMask("DETECTED")
676 fake_pixels = (science_fake_img.mask.array & detected).nonzero()
677 science.mask.array[fake_pixels] |= science_fake_bitmask
678 detected = tmplt_fake_img.mask.getPlaneBitMask("DETECTED")
679 fake_pixels = (tmplt_fake_img.mask.array & detected).nonzero()
680 template.mask.array[fake_pixels] |= science_fake_bitmask
682 science_fake_masked = (science.mask.array & science_fake_bitmask) > 0
683 template_fake_masked = (template.mask.array & template_fake_bitmask) > 0
685 subtractConfig = subtractImages.AlardLuptonSubtractTask.ConfigClass()
686 subtractConfig.sourceSelector.signalToNoise.fluxField = "truth_instFlux"
687 subtractConfig.sourceSelector.signalToNoise.errField = "truth_instFluxErr"
688 subtractTask = subtractImages.AlardLuptonSubtractTask(config=subtractConfig)
689 subtraction = subtractTask.run(template, science, sources)
691 # check subtraction mask plane is set where we set the previous masks
692 diff_mask = subtraction.difference.mask
694 # science mask should be now in INJECTED
695 inj_masked = (diff_mask.array & diff_mask.getPlaneBitMask("INJECTED")) > 0
697 # template mask should be now in INJECTED_TEMPLATE
698 injTmplt_masked = (diff_mask.array & diff_mask.getPlaneBitMask("INJECTED_TEMPLATE")) > 0
700 self.assertFloatsEqual(inj_masked.astype(int), science_fake_masked.astype(int))
701 # The template is convolved, so the INJECTED_TEMPLATE mask plane may
702 # include more pixels than the FAKE mask plane
703 injTmplt_masked &= template_fake_masked
704 self.assertFloatsEqual(injTmplt_masked.astype(int), template_fake_masked.astype(int))
706 # Now check that detection of fakes have the correct flag for injections
707 detectionTask = self._setup_detection()
709 output = detectionTask.run(subtraction.matchedScience,
710 subtraction.matchedTemplate,
711 subtraction.difference,
712 subtraction.kernelSources)
713 diaSources = output.diaSources[~output.diaSources["sky_source"]].copy(deep=True)
715 sci_refIds = []
716 tmpl_refIds = []
717 for diaSrc in diaSources:
718 if diaSrc['base_PsfFlux_instFlux'] > 0:
719 self._check_diaSource(science_fake_sources, diaSrc, scale=1, refIds=sci_refIds)
720 self.assertTrue(diaSrc['base_PixelFlags_flag_injected'])
721 self.assertTrue(diaSrc['base_PixelFlags_flag_injectedCenter'])
722 self.assertFalse(diaSrc['base_PixelFlags_flag_injected_template'])
723 self.assertFalse(diaSrc['base_PixelFlags_flag_injected_templateCenter'])
724 else:
725 self._check_diaSource(tmplt_fake_sources, diaSrc, scale=-1, refIds=tmpl_refIds)
726 self.assertTrue(diaSrc['base_PixelFlags_flag_injected_template'])
727 self.assertTrue(diaSrc['base_PixelFlags_flag_injected_templateCenter'])
728 self.assertFalse(diaSrc['base_PixelFlags_flag_injected'])
729 self.assertFalse(diaSrc['base_PixelFlags_flag_injectedCenter'])
731 def test_mask_streaks(self):
732 """Run detection on a difference image containing a streak.
733 """
734 # Set up the simulated images
735 noiseLevel = 1.
736 staticSeed = 1
737 fluxLevel = 500
738 xSize = 400
739 ySize = 400
740 psfSize = 2.4
741 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel,
742 "xSize": xSize, "ySize": ySize}
743 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
744 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
746 # Configure the detection Task
747 detectionTask = self._setup_detection(doMerge=False, doMaskStreaks=True,
748 raiseOnNoDiaSources=False, raiseOnBadSubtractionRatio=False)
750 # Test that no streaks are detected
751 difference = science.clone()
752 difference.maskedImage -= matchedTemplate.maskedImage
753 output = detectionTask.run(science, matchedTemplate, difference, sources)
754 outMask = output.subtractedMeasuredExposure.mask.array
755 streakMask = output.subtractedMeasuredExposure.mask.getPlaneBitMask("STREAK")
756 streakMaskSet = (outMask & streakMask) > 0
757 self.assertTrue(np.all(streakMaskSet == 0))
759 # Add streak-like shape and check that streak is detected
760 difference.image.array[20:23, 40:200] += 50
761 output = detectionTask.run(science, matchedTemplate, difference, sources)
762 outMask = output.subtractedMeasuredExposure.mask.array
763 streakMask = output.subtractedMeasuredExposure.mask.getPlaneBitMask("STREAK")
764 streakMaskSet = (outMask & streakMask) > 0
765 self.assertTrue(np.all(streakMaskSet[20:23, 40:200]))
767 # Add a more trapezoid shaped streak across an image that is
768 # fainter and check that it is also detected
769 bbox = science.getBBox()
770 difference = science.clone()
771 difference.maskedImage -= matchedTemplate.maskedImage
772 width = 100
773 amplitude = 5
774 x0 = -100 + bbox.getBeginX()
775 y0 = -100 + bbox.getBeginY()
776 x1 = xSize + x0 + 100
777 y1 = ySize/2
778 corner_coords = [lsst.geom.Point2D(x0, y0),
779 lsst.geom.Point2D(x0, y0 + width),
780 lsst.geom.Point2D(x1, y1),
781 lsst.geom.Point2D(x1 + width, y1)]
782 streak_trapezoid = afwGeom.Polygon(corner_coords)
783 streak_image = streak_trapezoid.createImage(bbox)
784 streak_image *= amplitude
785 difference.image.array += streak_image.array
786 output = detectionTask.run(science, matchedTemplate, difference, sources)
787 outMask = output.subtractedMeasuredExposure.mask.array
788 streakMask = output.subtractedMeasuredExposure.mask.getPlaneBitMask("STREAK")
789 streakMaskSet = (outMask & streakMask) > 0
790 # Check that pixel values in streak_image equal the streak amplitude
791 streak_check = np.where(streak_image.array == amplitude)
792 self.assertTrue(np.all(streakMaskSet[streak_check]))
793 # Check that the entire image was not masked STREAK
794 self.assertFalse(np.all(streakMaskSet))
796 def _setup_sattle_tests(self):
797 noiseLevel = 1.
798 staticSeed = 1
799 fluxLevel = 500
800 psfSize = 2.4
801 shared_kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel,
802 "x0": 12345, "y0": 67890}
803 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6,
804 **shared_kwargs)
805 science.getInfo().setVisitInfo(makeVisitInfo(id=2))
806 detector = DetectorWrapper(numAmps=1).detector
807 science.setDetector(detector)
808 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel / 4,
809 noiseSeed=7, **shared_kwargs)
810 difference = science.clone()
812 detectionTask = self._setup_detection(doDeblend=True,
813 badSubtractionRatioThreshold=1.,
814 doSkySources=False, run_sattle=True)
816 return science, matchedTemplate, difference, sources, detectionTask
818 @mock.patch.dict(os.environ, {"SATTLE_URI_BASE": "fake_host:1234"})
819 def test_sattle_not_available(self):
820 science, matchedTemplate, difference, sources, detectionTask = self._setup_sattle_tests()
822 response = MockResponse({"allow_list": []}, 500, "sattle internal error")
823 with mock.patch('requests.put', return_value=response):
824 with self.assertRaises(requests.exceptions.HTTPError):
825 detectionTask.run(science, matchedTemplate, difference, sources,
826 idFactory=IdFactory.makeSimple())
828 @mock.patch.dict(os.environ, {"SATTLE_URI_BASE": "fake_host:1234"})
829 def test_visit_id_not_in_sattle(self):
830 science, matchedTemplate, difference, sources, detectionTask = self._setup_sattle_tests()
832 response = MockResponse({"allow_list": []}, 404, "missing visit cache")
833 # visit id not in sattle raises
834 with self.assertRaises(requests.exceptions.HTTPError):
835 with mock.patch('lsst.ip.diffim.detectAndMeasure.requests.put',
836 return_value=response):
837 with mock.patch('lsst.ip.diffim.utils.populate_sattle_visit_cache'):
838 detectionTask.run(science, matchedTemplate, difference, sources,
839 idFactory=IdFactory.makeSimple())
841 @mock.patch.dict(os.environ, {"SATTLE_URI_BASE": "fake_host:1234"})
842 def test_filter_satellites_some_allowed(self):
843 science, matchedTemplate, difference, sources, detectionTask = self._setup_sattle_tests()
845 allowed_ids = [1, 5]
846 response = MockResponse({"allow_list": allowed_ids}, 200, "some allowed")
847 with mock.patch('requests.put', return_value=response):
848 output = detectionTask.run(science, matchedTemplate, difference, sources,
849 idFactory=IdFactory.makeSimple())
851 self.assertEqual(len(output.diaSources), 2)
853 # Output should be sources 1 and 5 allowed out of 20
854 self.assertEqual(set(output.diaSources['id']), set(allowed_ids))
856 @mock.patch.dict(os.environ, {"SATTLE_URI_BASE": "fake_host:1234"})
857 def test_filter_satellites_all_allowed(self):
858 science, matchedTemplate, difference, sources, detectionTask = self._setup_sattle_tests()
860 allowed_ids = list(range(1, 21))
861 response = MockResponse({"allow_list": allowed_ids}, 200, "all allowed")
862 # Run detection and check the results
863 with mock.patch('requests.put', return_value=response):
864 output = detectionTask.run(science, matchedTemplate, difference, sources,
865 idFactory=IdFactory.makeSimple())
867 # Output should be all sources that went in. 20 go in, 20 should come out
868 self.assertEqual(len(output.diaSources), 20)
870 self.assertEqual(set(output.diaSources['id']), set(allowed_ids))
872 @mock.patch.dict(os.environ, {"SATTLE_URI_BASE": "fake_host:1234"})
873 def test_filter_satellites_none_allowed(self):
874 science, matchedTemplate, difference, sources, detectionTask = self._setup_sattle_tests()
876 response = MockResponse({"allow_list": []}, 200, "none allowed")
877 # Run detection and confirm it raises for no diasources
878 with self.assertRaises(detectAndMeasure.NoDiaSourcesError):
879 with mock.patch('requests.put', return_value=response):
880 detectionTask.run(science, matchedTemplate, difference, sources,
881 idFactory=IdFactory.makeSimple())
883 @mock.patch.dict(os.environ, {"SATTLE_URI_BASE": ""})
884 def test_fail_on_sattle_misconfiguration(self):
885 with self.assertRaises(pexConfig.FieldValidationError):
886 self._setup_detection(run_sattle=True)
888 def test_trailed_glints(self):
889 """Test that the glint_trail column works, and that
890 the trailed_glints output contains the expected information.
891 """
892 noiseLevel = 1.
893 staticSeed = 1
894 diffim, diaSources = makeTestImage(seed=staticSeed, noiseLevel=noiseLevel, noiseSeed=6)
895 self._check_values(diaSources['glint_trail'])
897 # Run detection and return the output Struct so we can check it
898 def _detection_wrapper(diffim, diaSources):
899 detectionTask = self._setup_detection()
900 scienceBase, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6)
901 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7)
902 science = scienceBase.clone()
903 science.maskedImage -= diffim.maskedImage
904 difference = science.clone()
905 difference.maskedImage -= matchedTemplate.maskedImage
906 output = detectionTask.run(science, matchedTemplate, difference, sources)
907 return output
909 output = _detection_wrapper(diffim, diaSources)
910 self.assertTrue('slopes' in output.glintTrailInfo)
911 self.assertTrue('intercepts' in output.glintTrailInfo)
912 self.assertTrue('stderrs' in output.glintTrailInfo)
913 self.assertTrue('lengths' in output.glintTrailInfo)
914 self.assertTrue('angles' in output.glintTrailInfo)
917class DetectAndMeasureScoreTest(DetectAndMeasureTestBase, lsst.utils.tests.TestCase):
918 detectionTask = detectAndMeasure.DetectAndMeasureScoreTask
920 def test_detection_xy0(self):
921 """Basic functionality test with non-zero x0 and y0.
922 """
923 # Set up the simulated images
924 noiseLevel = 1.
925 staticSeed = 1
926 fluxLevel = 500
927 kernelSize = 31
928 psfSize = 2.4
929 # Buffer source positions away from the score image's EDGE strip.
930 # Preconvolution adds ~kernelSize//2 pixels of EDGE on top of the
931 # EDGE already set by detection smoothing on the science image.
932 templateBorderSize = kernelSize//2
933 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel, "x0": 12345, "y0": 67890,
934 "kernelSize": kernelSize, "templateBorderSize": templateBorderSize}
935 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
936 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
937 difference = science.clone()
938 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
939 scienceKernel = science.psf.getKernel()
940 score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
942 # Configure the detection Task
943 # Set the subtraction residual metric threshold high, since the template
944 # is not actually subtracted from the science image.
945 detectionTask = self._setup_detection(doDeblend=True, badSubtractionRatioThreshold=1.,
946 doSkySources=False)
948 # Run detection and check the results
949 output = detectionTask.run(science, matchedTemplate, difference, score, sources,
950 idFactory=self.idGenerator.make_table_id_factory())
952 # Catalog ids should be very large from this id generator.
953 self.assertTrue(all(output.diaSources['id'] > 1000000000))
954 subtractedMeasuredExposure = output.subtractedMeasuredExposure
956 self.assertImagesEqual(subtractedMeasuredExposure.image, difference.image)
958 self.assertEqual(len(output.diaSources), len(sources))
959 # no sources should be flagged as negative
960 self.assertEqual(len(~output.diaSources["is_negative"]), len(output.diaSources))
961 refIds = []
962 for source in sources:
963 self._check_diaSource(output.diaSources, source, refIds=refIds)
965 def test_measurements_finite(self):
966 """Measured fluxes and centroids should always be finite.
967 """
968 columnNames = ["coord_ra", "coord_dec", "ip_diffim_forced_PsfFlux_instFlux",
969 "ip_diffim_forced_template_PsfFlux_instFlux"]
971 # Set up the simulated images
972 noiseLevel = 1.
973 staticSeed = 1
974 transientSeed = 6
975 xSize = 256
976 ySize = 256
977 psfSize = 2.4
978 kwargs = {"psfSize": psfSize, "x0": 0, "y0": 0,
979 "xSize": xSize, "ySize": ySize}
980 science, sources = makeTestImage(seed=staticSeed, noiseLevel=noiseLevel, noiseSeed=6,
981 nSrc=1, **kwargs)
982 matchedTemplate, _ = makeTestImage(seed=staticSeed, noiseLevel=noiseLevel/4, noiseSeed=7,
983 nSrc=1, **kwargs)
984 rng = np.random.RandomState(3)
985 xLoc = np.arange(-5, xSize+5, 10)
986 rng.shuffle(xLoc)
987 yLoc = np.arange(-5, ySize+5, 10)
988 rng.shuffle(yLoc)
989 transients, transientSources = makeTestImage(seed=transientSeed,
990 nSrc=len(xLoc), fluxLevel=1000.,
991 noiseLevel=noiseLevel, noiseSeed=8,
992 xLoc=xLoc, yLoc=yLoc,
993 **kwargs)
994 difference = science.clone()
995 difference.maskedImage -= matchedTemplate.maskedImage
996 difference.maskedImage += transients.maskedImage
997 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
998 scienceKernel = science.psf.getKernel()
999 score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
1001 # Configure the detection Task
1002 detectionTask = self._setup_detection(doForcedMeasurement=True)
1004 # Run detection and check the results
1005 output = detectionTask.run(science, matchedTemplate, difference, score, sources)
1007 for column in columnNames:
1008 self._check_values(output.diaSources[column])
1009 self._check_values(output.diaSources.getX(), minValue=0, maxValue=xSize)
1010 self._check_values(output.diaSources.getY(), minValue=0, maxValue=ySize)
1011 self._check_values(output.diaSources.getPsfInstFlux())
1013 def test_detect_transients(self):
1014 """Run detection on a difference image containing transients.
1015 """
1016 # Set up the simulated images
1017 noiseLevel = 1.
1018 staticSeed = 1
1019 transientSeed = 6
1020 nTransients = 10
1021 transientFlux = 1000
1022 fluxLevel = 500
1023 psfSize = 2.4
1024 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel}
1025 scienceBase, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
1026 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
1027 scienceKernel = scienceBase.psf.getKernel()
1028 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
1030 # Configure the detection Task
1031 detectionTask = self._setup_detection(doMerge=False)
1032 kwargs["seed"] = transientSeed
1033 kwargs["nSrc"] = nTransients
1034 kwargs["fluxLevel"] = transientFlux
1036 # Run detection and check the results
1037 def _run_and_check_detections(positive=True):
1038 """Simulate positive or negative transients and run detection.
1040 Parameters
1041 ----------
1042 positive : `bool`, optional
1043 If set, use positive transient sources.
1044 """
1046 transients, transientSources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, **kwargs)
1047 science = scienceBase.clone()
1048 if positive:
1049 science.maskedImage += transients.maskedImage
1050 else:
1051 science.maskedImage -= transients.maskedImage
1052 difference = science.clone()
1053 difference.maskedImage -= matchedTemplate.maskedImage
1054 score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
1055 # NOTE: NoiseReplacer (run by forcedMeasurement) can modify the
1056 # science image if we've e.g. removed parents post-deblending.
1057 # Pass a clone of the science image, so that it doesn't disrupt
1058 # later tests.
1059 output = detectionTask.run(science, matchedTemplate, difference, score, sources)
1060 refIds = []
1061 scale = 1. if positive else -1.
1062 # sources near the edge may have untrustworthy centroids
1063 goodSrcFlags = ~output.diaSources['base_PixelFlags_flag_edge']
1064 for diaSource, goodSrcFlag in zip(output.diaSources, goodSrcFlags):
1065 if goodSrcFlag:
1066 self._check_diaSource(transientSources, diaSource, refIds=refIds, scale=scale)
1067 _run_and_check_detections(positive=True)
1068 _run_and_check_detections(positive=False)
1070 def test_detect_dipoles(self):
1071 """Run detection on a difference image containing dipoles.
1072 """
1073 # Set up the simulated images
1074 noiseLevel = 1.
1075 staticSeed = 1
1076 fluxLevel = 1000
1077 fluxRange = 1.5
1078 nSources = 10
1079 offset = 1
1080 xSize = 300
1081 ySize = 300
1082 kernelSize = 31
1083 psfSize = 2.4
1084 # Avoid placing sources near the edge for this test, so that we can
1085 # easily check that the correct number of sources are detected.
1086 templateBorderSize = kernelSize//2
1087 dipoleFlag = "ip_diffim_DipoleFit_classification"
1088 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel, "fluxRange": fluxRange,
1089 "nSrc": nSources, "templateBorderSize": templateBorderSize, "kernelSize": kernelSize,
1090 "xSize": xSize, "ySize": ySize}
1091 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
1092 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
1093 difference = science.clone()
1094 # Shift the template by a pixel in order to make dipoles in the difference image.
1095 matchedTemplate.image.array[...] = np.roll(matchedTemplate.image.array[...], offset, axis=0)
1096 matchedTemplate.variance.array[...] = np.roll(matchedTemplate.variance.array[...], offset, axis=0)
1097 matchedTemplate.mask.array[...] = np.roll(matchedTemplate.mask.array[...], offset, axis=0)
1098 difference.maskedImage -= matchedTemplate.maskedImage[science.getBBox()]
1099 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
1100 scienceKernel = science.psf.getKernel()
1101 score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
1103 detectionTask = self._setup_detection(badSubtractionRatioThreshold=0.3)
1104 output = detectionTask.run(science, matchedTemplate, difference, score, sources)
1105 diaSources = output.diaSources[~output.diaSources["sky_source"]].copy(deep=True)
1106 self.assertEqual(len(diaSources), len(sources))
1107 refIds = []
1108 for diaSource in diaSources:
1109 if diaSource[dipoleFlag]:
1110 self._check_diaSource(sources, diaSource, refIds=refIds, scale=0,
1111 rtol=0.05, atol=None, usePsfFlux=False)
1112 self.assertFloatsAlmostEqual(diaSource["ip_diffim_DipoleFit_orientation"],
1113 -np.pi / 2, atol=2.)
1114 self.assertFloatsAlmostEqual(diaSource["ip_diffim_DipoleFit_separation"], offset, rtol=0.1)
1115 else:
1116 raise ValueError("DiaSource with ID %s is not a dipole!", diaSource.getId())
1118 def test_sky_sources(self):
1119 """Add sky sources and check that they are sufficiently far from other
1120 sources and have negligible flux.
1121 """
1122 # Set up the simulated images
1123 noiseLevel = 1.
1124 staticSeed = 1
1125 transientSeed = 6
1126 transientFluxLevel = 1000.
1127 transientFluxRange = 1.5
1128 fluxLevel = 500
1129 psfSize = 2.4
1130 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel}
1131 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
1132 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
1133 transients, transientSources = makeTestImage(seed=transientSeed, psfSize=2.4,
1134 nSrc=10, fluxLevel=transientFluxLevel,
1135 fluxRange=transientFluxRange,
1136 noiseLevel=noiseLevel, noiseSeed=8)
1137 difference = science.clone()
1138 difference.maskedImage -= matchedTemplate.maskedImage
1139 difference.maskedImage += transients.maskedImage
1140 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
1141 scienceKernel = science.psf.getKernel()
1142 kernelWidth = np.max(scienceKernel.getDimensions())//2
1143 score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
1145 # Configure the detection Task
1146 detectionTask = self._setup_detection(doSkySources=True)
1148 # Run detection and check the results
1149 output = detectionTask.run(science, matchedTemplate, difference, score, sources,
1150 idFactory=self.idGenerator.make_table_id_factory())
1151 nSkySourcesGenerated = detectionTask.metadata["n_skySources"]
1152 skySources = output.diaSources[output.diaSources["sky_source"]]
1153 self.assertEqual(len(skySources), nSkySourcesGenerated)
1154 for skySource in skySources:
1155 # Disable the sky_source flag to enable checks.
1156 skySource["sky_source"] = False
1157 # The sky sources should not be close to any other source
1158 with self.assertRaises(AssertionError):
1159 self._check_diaSource(transientSources, skySource, matchDistance=kernelWidth)
1160 with self.assertRaises(AssertionError):
1161 self._check_diaSource(sources, skySource, matchDistance=kernelWidth)
1162 # The sky sources should have low flux levels.
1163 self._check_diaSource(transientSources, skySource, matchDistance=1000, scale=0.,
1164 atol=np.sqrt(transientFluxRange*transientFluxLevel))
1166 # Catalog ids should be very large from this id generator.
1167 self.assertTrue(all(output.diaSources['id'] > 1000000000))
1169 def test_exclude_mask_detections(self):
1170 """Sources with certain bad mask planes set should not be detected.
1171 """
1172 # Set up the simulated images
1173 noiseLevel = 1.
1174 staticSeed = 1
1175 transientSeed = 6
1176 fluxLevel = 500
1177 radius = 2
1178 psfSize = 2.4
1179 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel}
1180 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
1181 science.getInfo().setVisitInfo(makeVisitInfo())
1182 detector = DetectorWrapper(numAmps=1).detector
1183 science.setDetector(detector)
1184 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
1186 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
1187 scienceKernel = science.psf.getKernel()
1188 # Configure the detection Task
1189 badSourceFlags = ["base_PixelFlags_flag_offimage",
1190 "base_PixelFlags_flag_edgeCenterAll",
1191 "base_PixelFlags_flag_badCenterAll",
1192 "base_PixelFlags_flag_saturatedCenterAll", ]
1193 detectionTask = self._setup_detection(nSkySources=0, badSourceFlags=badSourceFlags)
1194 excludeMaskPlanes = ["EDGE", "BAD", "SAT"]
1195 nBad = len(excludeMaskPlanes)
1196 self.assertGreater(nBad, 0)
1197 kwargs["seed"] = transientSeed
1198 kwargs["nSrc"] = nBad
1199 kwargs["fluxLevel"] = 1000
1201 # Run detection and check the results
1202 def _run_and_check_detections(setFlags=True):
1203 transients, transientSources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, **kwargs)
1204 difference = science.clone()
1205 difference.maskedImage -= matchedTemplate.maskedImage
1206 difference.maskedImage += transients.maskedImage
1207 if setFlags:
1208 for src, badMask in zip(transientSources, excludeMaskPlanes):
1209 srcX = int(src.getX())
1210 srcY = int(src.getY())
1211 srcBbox = lsst.geom.Box2I(lsst.geom.Point2I(srcX - radius, srcY - radius),
1212 lsst.geom.Extent2I(2*radius + 1, 2*radius + 1))
1213 difference[srcBbox].mask.array |= lsst.afw.image.Mask.getPlaneBitMask(badMask)
1214 score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
1215 if setFlags:
1216 with self.assertRaises(AlgorithmError):
1217 output = detectionTask.run(science, matchedTemplate, difference, score, sources)
1218 else:
1219 output = detectionTask.run(science, matchedTemplate, difference, score, sources)
1220 refIds = []
1221 goodSrcFlags = checkMask(difference.mask, transientSources, excludeMaskPlanes)
1222 self.assertEqual(np.sum(~goodSrcFlags), 0)
1223 for diaSource, goodSrcFlag in zip(output.diaSources, goodSrcFlags):
1224 if ~goodSrcFlag:
1225 with self.assertRaises(AssertionError):
1226 self._check_diaSource(transientSources, diaSource, refIds=refIds)
1227 else:
1228 self._check_diaSource(transientSources, diaSource, refIds=refIds)
1229 _run_and_check_detections(setFlags=False)
1230 _run_and_check_detections(setFlags=True)
1232 def test_detect_transients_with_background(self):
1233 """Background measured on the difference image should be subtracted
1234 from the score image so that transients are still detected cleanly.
1235 """
1236 # Set up the simulated images
1237 noiseLevel = 1.
1238 staticSeed = 1
1239 transientSeed = 6
1240 nTransients = 10
1241 transientFlux = 1000
1242 fluxLevel = 500
1243 xSize = 512
1244 ySize = 512
1245 x0 = 123
1246 y0 = 456
1247 kernelSize = 31
1248 psfSize = 2.4
1249 # Avoid placing sources near the edge for this test, so that we can
1250 # easily check that the correct number of sources are detected.
1251 templateBorderSize = kernelSize//2
1252 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel,
1253 "xSize": xSize, "ySize": ySize, "x0": x0, "y0": y0,
1254 "kernelSize": kernelSize, "templateBorderSize": templateBorderSize}
1255 chebyshev_params = [2.2, 2.1, 2.0, 1.2, 1.1, 1.0]
1257 # The background model must cover the grown image bbox, otherwise it
1258 # would be extrapolated outside its declared domain.
1259 bbox2D = lsst.geom.Box2D(
1260 lsst.geom.Point2D(x0 - templateBorderSize, y0 - templateBorderSize),
1261 lsst.geom.Extent2D(xSize + 2*templateBorderSize, ySize + 2*templateBorderSize))
1262 background_model = afwMath.Chebyshev1Function2D(chebyshev_params, bbox2D)
1263 scienceBase, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6,
1264 background=background_model, **kwargs)
1265 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
1266 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
1267 scienceKernel = scienceBase.psf.getKernel()
1269 # Configure the detection Task
1270 detectionTask = self._setup_detection(doMerge=False, doSubtractBackground=True)
1271 kwargs["seed"] = transientSeed
1272 kwargs["nSrc"] = nTransients
1273 kwargs["fluxLevel"] = transientFlux
1275 def _run_and_check_detections(positive=True):
1276 transients, transientSources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, **kwargs)
1277 science = scienceBase.clone()
1278 if positive:
1279 science.maskedImage += transients.maskedImage
1280 else:
1281 science.maskedImage -= transients.maskedImage
1282 difference = science.clone()
1283 difference.maskedImage -= matchedTemplate.maskedImage
1284 score = subtractTask._convolveExposure(difference, scienceKernel,
1285 subtractTask.convolutionControl)
1286 # Record the score image metric before detection so we can verify
1287 # that the background was actually subtracted from it.
1288 originalScoreMedian = np.median(score.image.array[np.isfinite(score.image.array)])
1289 originalDiffimMedian = np.median(difference.image.array[np.isfinite(difference.image.array)])
1290 output = detectionTask.run(science, matchedTemplate, difference, score, sources)
1292 # The background subtraction should leave a non-empty BackgroundList.
1293 self.assertGreater(len(output.differenceBackground), 0)
1294 # The difference image should have been modified in place by the
1295 # background subtraction, so its median should shift noticeably.
1296 subtractedScoreMedian = np.median(score.image.array[np.isfinite(score.image.array)])
1297 self.assertLess(np.abs(subtractedScoreMedian), np.abs(originalScoreMedian))
1298 subtractedDiffimMedian = np.median(difference.image.array[np.isfinite(difference.image.array)])
1299 self.assertLess(np.abs(subtractedDiffimMedian), np.abs(originalDiffimMedian))
1301 refIds = []
1302 scale = 1. if positive else -1.
1303 for transient in transientSources:
1304 self._check_diaSource(output.diaSources, transient, refIds=refIds, scale=scale)
1305 _run_and_check_detections(positive=True)
1306 _run_and_check_detections(positive=False)
1308 def test_mask_cosmic_rays(self):
1309 """Cosmic rays detected on the difference image should propagate
1310 to the mask of the returned (measured) exposure. This version creates
1311 and uses a Score image for detection.
1312 """
1313 # Set up the simulated images
1314 noiseLevel = 1.
1315 staticSeed = 1
1316 fluxLevel = 500
1317 xSize = 400
1318 ySize = 400
1319 psfSize = 2.4
1320 kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel,
1321 "xSize": xSize, "ySize": ySize}
1322 science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
1323 matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
1324 subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
1325 scienceKernel = science.psf.getKernel()
1326 crMask = science.mask.getPlaneBitMask("CR")
1328 # Configure the detection Task
1329 detectionTask = self._setup_detection(doMerge=False, doSkySources=True)
1331 # Test that no CRs are detected when none are present.
1332 transients, _ = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, nSrc=10, **kwargs)
1333 science.maskedImage += transients.maskedImage
1334 difference = science.clone()
1335 difference.maskedImage -= matchedTemplate.maskedImage
1336 score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
1337 output = detectionTask.run(science, matchedTemplate, difference, score, sources)
1338 self.assertFalse(np.any((output.subtractedMeasuredExposure.mask.array & crMask) > 0))
1340 crX0 = round(sources.getX()[0] - science.getX0())
1341 crY0 = round(sources.getY()[0] - science.getY0())
1342 crX1 = round(sources.getX()[1] - science.getX0())
1343 crY1 = round(sources.getY()[1] - science.getY0())
1344 # Inject CR-like shapes into the difference image. CR detection runs
1345 # on the difference, and the mask should propagate to the measured
1346 # exposure returned by the task.
1347 crMaskSetInput = np.zeros(difference.image.array.shape, bool)
1348 crMaskSetInput[crX0:crX0+1, crY0:crY0+5] = True
1349 crMaskSetInput[crX1:crX1+5, crY1:crY1+1] = True
1350 difference.image.array[crMaskSetInput] += 1234
1351 score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
1352 output = detectionTask.run(science, matchedTemplate, difference, score, sources)
1353 crMaskSet = (output.subtractedMeasuredExposure.mask.array & crMask) > 0
1354 self.assertFalse(np.any(crMaskSet[~crMaskSetInput]))
1355 self.assertTrue(np.all(crMaskSet[crMaskSetInput]))
1358class TestNegativePeaks(lsst.utils.tests.TestCase):
1359 """Tests of deblending and merging negative peaks, to test fixes for the
1360 various problems found on DM-48596.
1361 """
1363 def testDeblendNegatives(self):
1364 """Test that negative peaks get deblended and not destroyed: DM-48704.
1365 This is only a test of deblending, not of merging.
1366 """
1367 # Make a science image with one blend of two positive sources, to
1368 # subtract from an empty template, resulting in a negative diffim.
1369 bbox = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(100, 100))
1370 dataset = lsst.meas.base.tests.TestDataset(bbox)
1371 delta = 10
1372 with dataset.addBlend() as family:
1373 family.addChild(instFlux=2E5, centroid=lsst.geom.Point2D(50, 72))
1374 family.addChild(instFlux=2.5E5, centroid=lsst.geom.Point2D(50+delta, 74))
1375 science, catalog = dataset.realize(noise=1.0,
1376 schema=lsst.meas.base.tests.TestDataset.makeMinimalSchema())
1377 dataset = lsst.meas.base.tests.TestDataset(bbox)
1378 template, _ = dataset.realize(noise=1.0,
1379 schema=lsst.meas.base.tests.TestDataset.makeMinimalSchema())
1380 difference = template.clone()
1381 difference.image -= science.image
1383 config = detectAndMeasure.DetectAndMeasureTask.ConfigClass()
1384 config.doDeblend = True
1385 task = detectAndMeasure.DetectAndMeasureTask(config=config)
1386 # prelude steps taken from `detectAndMeasure.run`
1387 task._prepareInputs(difference)
1388 table = lsst.afw.table.SourceTable.make(task.schema)
1389 results = task.detection.run(table=table, exposure=difference, doSmooth=True)
1390 # Just run the deblend step so we can check the footprints independently.
1391 sources, positives, negatives = task._deblend(difference, results.positive, results.negative)
1393 # DM-48704 fixed a problem where the peaks were in the footprints, but
1394 # the spans were empty.
1395 self.assertEqual(len(negatives), 2)
1396 self.assertEqual(len(positives), 0)
1397 self.assertGreater(negatives[0].getFootprint().getSpans().getArea(), 0)
1398 self.assertGreater(negatives[1].getFootprint().getSpans().getArea(), 0)
1399 # Deblended children are HeavyFootprints; we have to make sure the
1400 # pixel values in those are correct (though DetectAndMeasureTask
1401 # doesn't use the fact that they're Heavy).
1402 # The sources are positive in the science image, and negative in the
1403 # diffim, so the minimum value in the deblended negative footprint is
1404 # the maximum value in the science catalog footprint (ignoring the
1405 # noise in the science and template, hence rtol).
1406 # (catalog[0] is the parent; we want the children)
1407 self.assertFloatsAlmostEqual(negatives[0].getFootprint().getImageArray().min(),
1408 -catalog[1].getFootprint().getImageArray().max(), rtol=1e-4)
1409 self.assertFloatsAlmostEqual(negatives[1].getFootprint().getImageArray().min(),
1410 -catalog[2].getFootprint().getImageArray().max(), rtol=1e-4)
1412 self.assertEqual(sources["is_negative"].sum(), 2)
1414 def testMergeFootprints(self):
1415 """Test that merging footprints does not cause negative ones to
1416 disappear (e.g. get merged into non-connected footprints).
1418 As implemented, the diffim will have three positive sources (one a
1419 blended pair), and 7 negative sources (two a blended pair).
1420 """
1421 # Make a template image with multiple blends, designed to trigger the
1422 # negative-footprint-related bug in `FootprintSet.merge`.
1423 bbox = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(400, 200))
1424 dataset = lsst.meas.base.tests.TestDataset(bbox)
1425 # Because these sources are on the template, they will be negative on
1426 # the diffim, unless the pixels are explicitly made negative via
1427 # `template.image.subset` below.
1428 # positive isolated source on diffim
1429 dataset.addSource(instFlux=.5E5, centroid=lsst.geom.Point2D(25, 26))
1430 # negative isolated source on diffim
1431 dataset.addSource(instFlux=.7E5, centroid=lsst.geom.Point2D(75, 24),
1432 shape=lsst.afw.geom.Quadrupole(12, 7, 2))
1433 delta = 10
1434 # negative blended pair on diffim
1435 with dataset.addBlend() as family:
1436 family.addChild(instFlux=1E5, centroid=lsst.geom.Point2D(150, 72))
1437 family.addChild(instFlux=1.5E5, centroid=lsst.geom.Point2D(150+delta, 74))
1438 # positive blended pair on diffim
1439 with dataset.addBlend() as family:
1440 family.addChild(instFlux=2E5, centroid=lsst.geom.Point2D(250, 72))
1441 family.addChild(instFlux=2.5E5, centroid=lsst.geom.Point2D(250+delta, 74))
1442 # negative blended pair on diffim
1443 with dataset.addBlend() as family:
1444 family.addChild(instFlux=3E5, centroid=lsst.geom.Point2D(350, 72))
1445 family.addChild(instFlux=3.5E5, centroid=lsst.geom.Point2D(350+delta, 74))
1446 # negative isolated source on diffim
1447 dataset.addSource(instFlux=4E5, centroid=lsst.geom.Point2D(75, 124))
1448 # negative isolated source on diffim
1449 dataset.addSource(instFlux=5E5, centroid=lsst.geom.Point2D(175, 124))
1450 schema = lsst.meas.base.tests.TestDataset.makeMinimalSchema()
1451 schema.addField("sky_source", "Flag", "Sky source.")
1452 template, catalog = dataset.realize(noise=100.0, schema=schema)
1453 # These will be positive sources on the diffim (as noted above).
1454 template.image.subset(lsst.geom.Box2I(lsst.geom.Point2I(15, 15),
1455 lsst.geom.Point2I(35, 35))).array *= -1
1456 template.image.subset(lsst.geom.Box2I(lsst.geom.Point2I(230, 60),
1457 lsst.geom.Point2I(270, 85))).array *= -1
1458 dataset = lsst.meas.base.tests.TestDataset(bbox)
1459 science, _ = dataset.realize(noise=100.0, schema=schema)
1460 difference = science.clone()
1461 difference.image -= template.image
1463 config = detectAndMeasure.DetectAndMeasureTask.ConfigClass()
1464 config.doDeblend = True
1465 config.raiseOnBadSubtractionRatio = False
1466 config.raiseOnNoDiaSources = False
1467 task = detectAndMeasure.DetectAndMeasureTask(config=config)
1468 result = task.run(science, template, difference, catalog)
1470 # The original bug manifested as the (175,124) source disappaering in
1471 # the merge step, and being included in the peaks of a duplicated
1472 # (75,124) source, even though their footprints are not connected.
1473 mask = np.isclose(result.diaSources["slot_Centroid_x"], 75) & \
1474 np.isclose(result.diaSources["slot_Centroid_y"], 124)
1475 self.assertEqual(mask.sum(), 1)
1476 peaks = result.diaSources[mask][0].getFootprint().peaks
1477 self.assertEqual(len(peaks), 1)
1478 self.assertEqual(peaks[0].getIx(), 75)
1479 self.assertEqual(peaks[0].getIy(), 124)
1481 mask = np.isclose(result.diaSources["slot_Centroid_x"], 175) & \
1482 np.isclose(result.diaSources["slot_Centroid_y"], 124)
1483 self.assertEqual(mask.sum(), 1)
1484 peaks = result.diaSources[mask][0].getFootprint().peaks
1485 self.assertEqual(len(peaks), 1)
1486 self.assertEqual(peaks[0].getIx(), 175)
1487 self.assertEqual(peaks[0].getIy(), 124)
1489 # This checks that all the returned footprints have exactly one peak;
1490 # it's a more general test than the above, but I'm keeping that for
1491 # the more explicit test of the original bugged sources.
1492 peak_x = np.array([peak['f_x'] for src in result.diaSources for peak in src.getFootprint().peaks])
1493 peak_y = np.array([peak['f_y'] for src in result.diaSources for peak in src.getFootprint().peaks])
1494 peaks = np.column_stack((peak_x, peak_y))
1495 unique_peaks, counts = np.unique(peaks, axis=0, return_counts=True)
1496 self.assertEqual(np.sum(counts > 1), 0)
1498 # There are three positive sources that should not have `is_negative`
1499 # set, independent of how deblending/merging of negative sources is
1500 # handled.
1501 self.assertEqual((~result.diaSources["is_negative"]).sum(), 3)
1504def makeVisitInfo(id=1):
1505 """Return a non-NaN visitInfo."""
1506 return afwImage.VisitInfo(id=id,
1507 exposureTime=10.01,
1508 darkTime=11.02,
1509 date=dafBase.DateTime(65321.1, dafBase.DateTime.MJD, dafBase.DateTime.TAI),
1510 ut1=12345.1,
1511 era=45.1*geom.degrees,
1512 boresightRaDec=geom.SpherePoint(23.1, 73.2, geom.degrees),
1513 boresightAzAlt=geom.SpherePoint(134.5, 33.3, geom.degrees),
1514 boresightAirmass=1.73,
1515 boresightRotAngle=73.2*geom.degrees,
1516 rotType=afwImage.RotType.SKY,
1517 observatory=Observatory(
1518 11.1*geom.degrees, 22.2*geom.degrees, 0.333),
1519 weather=Weather(1.1, 2.2, 34.5),
1520 )
1523class MockResponse:
1524 """Provide a mock for requests.put calls"""
1525 def __init__(self, json_data, status_code, text):
1526 self.json_data = json_data
1527 self.status_code = status_code
1528 self.text = text
1530 def json(self):
1531 return self.json_data
1533 def raise_for_status(self):
1534 if self.status_code != 200:
1535 raise requests.exceptions.HTTPError
1538def setup_module(module):
1539 lsst.utils.tests.init()
1542class MemoryTestCase(lsst.utils.tests.MemoryTestCase):
1543 pass
1546if __name__ == "__main__": 1546 ↛ 1547line 1546 didn't jump to line 1547 because the condition on line 1546 was never true
1547 lsst.utils.tests.init()
1548 unittest.main()