Skip to content

Commit

Permalink
Add score image detection and measurement tests
Browse files Browse the repository at this point in the history
  • Loading branch information
@isullivan
isullivan committed Jan 10, 2023
1 parent 4606bd0 commit 60ce4e1
Showing 1 changed file with 256 additions and 1 deletion.
257 changes: 256 additions & 1 deletion tests/test_detectAndMeasure.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@
import unittest

import lsst.geom
from lsst.ip.diffim import detectAndMeasure
from lsst.ip.diffim import detectAndMeasure, subtractImages
from lsst.ip.diffim.utils import makeTestImage
import lsst.utils.tests

Expand Down Expand Up @@ -265,6 +265,261 @@ def _check_diaSource(self, refSources, diaSource, refIds=None,
rtol=rtol, atol=atol)


class DetectAndMeasureScoreTest(lsst.utils.tests.TestCase):

def test_detection_runs(self):
"""Basic smoke test.
"""
noiseLevel = 1.
staticSeed = 1
fluxLevel = 500
kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel, "x0": 0, "y0": 0}
science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
difference = science.clone()
subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
scienceKernel = science.psf.getKernel()
score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
config = detectAndMeasure.DetectAndMeasureScoreTask.ConfigClass()
config.doApCorr = False
task = detectAndMeasure.DetectAndMeasureScoreTask(config=config)
output = task.run(science, matchedTemplate, difference, score)
subtractedMeasuredExposure = output.subtractedMeasuredExposure
self.assertImagesEqual(subtractedMeasuredExposure.image, difference.image)

def test_detection_xy0(self):
"""Basic smoke test with non-zero x0 and y0.
"""
noiseLevel = 1.
staticSeed = 1
fluxLevel = 500
kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel, "x0": 12345, "y0": 67890}
science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
difference = science.clone()
subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
scienceKernel = science.psf.getKernel()
score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
config = detectAndMeasure.DetectAndMeasureScoreTask.ConfigClass()
config.doApCorr = False
config.doMerge = False
config.doSkySources = False
config.doForcedMeasurement = False
task = detectAndMeasure.DetectAndMeasureScoreTask(config=config)
output = task.run(science, matchedTemplate, difference, score)
subtractedMeasuredExposure = output.subtractedMeasuredExposure

self.assertImagesEqual(subtractedMeasuredExposure.image, difference.image)

def test_measurements_finite(self):
"""Measured fluxes and centroids should always be finite.
"""
columnNames = ["coord_ra", "coord_dec", "ip_diffim_forced_PsfFlux_instFlux"]
noiseLevel = 1.
staticSeed = 1
transientSeed = 6
xSize = 256
ySize = 256
kwargs = {"psfSize": 2.4, "x0": 0, "y0": 0,
"xSize": xSize, "ySize": ySize}
science, sources = makeTestImage(seed=staticSeed, noiseLevel=noiseLevel, noiseSeed=6,
nSrc=1, **kwargs)
matchedTemplate, _ = makeTestImage(seed=staticSeed, noiseLevel=noiseLevel/4, noiseSeed=7,
nSrc=1, **kwargs)
rng = np.random.RandomState(3)
xLoc = np.arange(-5, xSize+5, 10)
rng.shuffle(xLoc)
yLoc = np.arange(-5, ySize+5, 10)
rng.shuffle(yLoc)
transients, transientSources = makeTestImage(seed=transientSeed,
nSrc=len(xLoc), fluxLevel=1000.,
noiseLevel=noiseLevel, noiseSeed=8,
xLoc=xLoc, yLoc=yLoc,
**kwargs)
difference = science.clone()
difference.maskedImage -= matchedTemplate.maskedImage
difference.maskedImage += transients.maskedImage
subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
scienceKernel = science.psf.getKernel()
score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
config = detectAndMeasure.DetectAndMeasureScoreTask.ConfigClass()
config.doApCorr = False
config.doMerge = False
config.doSkySources = False
config.doForcedMeasurement = True
task = detectAndMeasure.DetectAndMeasureScoreTask(config=config)
output = task.run(science, matchedTemplate, difference, score)

def _check_values(values, minValue=None, maxValue=None):
self.assertTrue(np.all(np.isfinite(values)))
if minValue is not None:
self.assertTrue(np.all(values > minValue))
if maxValue is not None:
self.assertTrue(np.all(values < maxValue))
for column in columnNames:
_check_values(output.diaSources[column])
_check_values(output.diaSources.getX(), minValue=0, maxValue=xSize)
_check_values(output.diaSources.getY(), minValue=0, maxValue=ySize)
_check_values(output.diaSources.getPsfInstFlux())

def test_detect_transients(self):
"""Run detection on a difference image containing transients.
"""
noiseLevel = 1.
staticSeed = 1
transientSeed = 6
fluxLevel = 500
kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel}
science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
config = detectAndMeasure.DetectAndMeasureScoreTask.ConfigClass()
config.doApCorr = False
config.doMerge = False
config.doSkySources = False
config.doForcedMeasurement = False
detectionTask = detectAndMeasure.DetectAndMeasureScoreTask(config=config)
subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
scienceKernel = science.psf.getKernel()

def _detection_wrapper(polarity=1):
transients, transientSources = makeTestImage(seed=transientSeed, psfSize=2.4,
nSrc=10, fluxLevel=1000.,
noiseLevel=noiseLevel, noiseSeed=8)
difference = science.clone()
difference.maskedImage -= matchedTemplate.maskedImage
if polarity > 0:
difference.maskedImage += transients.maskedImage
else:
difference.maskedImage -= transients.maskedImage
score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
output = detectionTask.run(science, matchedTemplate, difference, score)
refIds = []
for diaSource in output.diaSources:
self._check_diaSource(transientSources, diaSource, refIds=refIds, scale=polarity)
_detection_wrapper(polarity=1)
_detection_wrapper(polarity=-1)

def test_detect_dipoles(self):
"""Run detection on a difference image containing dipoles.
"""
noiseLevel = 1.
staticSeed = 1
fluxLevel = 1000
fluxRange = 1.5
nSources = 10
offset = 1
dipoleFlag = "ip_diffim_DipoleFit_flag_classification"
kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel, "fluxRange": fluxRange,
"nSrc": nSources}
science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
difference = science.clone()
matchedTemplate.image.array[...] = np.roll(matchedTemplate.image.array[...], offset, axis=0)
matchedTemplate.variance.array[...] = np.roll(matchedTemplate.variance.array[...], offset, axis=0)
matchedTemplate.mask.array[...] = np.roll(matchedTemplate.mask.array[...], offset, axis=0)
difference.maskedImage -= matchedTemplate.maskedImage[science.getBBox()]
subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
scienceKernel = science.psf.getKernel()
score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
config = detectAndMeasure.DetectAndMeasureScoreTask.ConfigClass()
config.doApCorr = False
config.doMerge = False
config.doSkySources = False
config.doForcedMeasurement = False
detectionTask = detectAndMeasure.DetectAndMeasureScoreTask(config=config)
output = detectionTask.run(science, matchedTemplate, difference, score)
self.assertIn(dipoleFlag, output.diaSources.schema.getNames())
nSourcesDet = len(sources)
self.assertEqual(len(output.diaSources), 2*nSourcesDet)
refIds = []
# The diaSource check should fail if we don't merge positive and negative footprints
for diaSource in output.diaSources:
with self.assertRaises(AssertionError):
self._check_diaSource(sources, diaSource, refIds=refIds, scale=0,
atol=np.sqrt(fluxRange*fluxLevel))

config.doMerge = True
detectionTask2 = detectAndMeasure.DetectAndMeasureScoreTask(config=config)
output2 = detectionTask2.run(science, matchedTemplate, difference, score)
self.assertEqual(len(output2.diaSources), nSourcesDet)
refIds = []
for diaSource in output2.diaSources:
if diaSource[dipoleFlag]:
self._check_diaSource(sources, diaSource, refIds=refIds, scale=0,
rtol=0.05, atol=None, usePsfFlux=False)
self.assertFloatsAlmostEqual(diaSource["ip_diffim_DipoleFit_orientation"], -90., atol=2.)
self.assertFloatsAlmostEqual(diaSource["ip_diffim_DipoleFit_separation"], offset, rtol=0.1)
else:
raise ValueError("DiaSource with ID %s is not a dipole!", diaSource.getId())

def test_sky_sources(self):
"""Add sky sources and check that they are sufficiently far from other
sources and have negligible flux.
"""
noiseLevel = 1.
staticSeed = 1
transientSeed = 6
transientFluxLevel = 1000.
transientFluxRange = 1.5
fluxLevel = 500
kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel}
science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
config = detectAndMeasure.DetectAndMeasureScoreTask.ConfigClass()
config.doApCorr = False
config.doMerge = False
config.doSkySources = True
config.doForcedMeasurement = False
config.skySources.nSources = 5
kernelWidth = np.max(science.psf.getKernel().getDimensions())//2
detectionTask = detectAndMeasure.DetectAndMeasureScoreTask(config=config)
transients, transientSources = makeTestImage(seed=transientSeed, psfSize=2.4,
nSrc=10, fluxLevel=transientFluxLevel,
fluxRange=transientFluxRange,
noiseLevel=noiseLevel, noiseSeed=8)
difference = science.clone()
difference.maskedImage -= matchedTemplate.maskedImage
difference.maskedImage += transients.maskedImage
subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
scienceKernel = science.psf.getKernel()
score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
output = detectionTask.run(science, matchedTemplate, difference, score)
skySources = output.diaSources[output.diaSources["sky_source"]]
self.assertEqual(len(skySources), config.skySources.nSources)
for skySource in skySources:
# The sky sources should not be close to any other source
with self.assertRaises(AssertionError):
self._check_diaSource(transientSources, skySource, matchDistance=kernelWidth)
with self.assertRaises(AssertionError):
self._check_diaSource(sources, skySource, matchDistance=kernelWidth)
# The sky sources should have low flux levels.
self._check_diaSource(transientSources, skySource, matchDistance=1000, scale=0.,
atol=np.sqrt(transientFluxRange*transientFluxLevel))

def _check_diaSource(self, refSources, diaSource, refIds=None,
matchDistance=1., scale=1., usePsfFlux=True,
rtol=0.02, atol=None):
"""Match a diaSource with a source in a reference catalog
and compare properties.
"""
distance = np.sqrt((diaSource.getX() - refSources.getX())**2
+ (diaSource.getY() - refSources.getY())**2)
self.assertLess(min(distance), matchDistance)
src = refSources[np.argmin(distance)]
if refIds is not None:
# Check that the same source was not previously associated
self.assertNotIn(src.getId(), refIds)
refIds.append(src.getId())
if atol is None:
atol = rtol*src.getPsfInstFlux() if usePsfFlux else rtol*src.getApInstFlux()
if usePsfFlux:
self.assertFloatsAlmostEqual(src.getPsfInstFlux()*scale, diaSource.getPsfInstFlux(),
rtol=rtol, atol=atol)
else:
self.assertFloatsAlmostEqual(src.getApInstFlux()*scale, diaSource.getApInstFlux(),
rtol=rtol, atol=atol)


def setup_module(module):
lsst.utils.tests.init()

Expand Down

0 comments on commit 60ce4e1

Please sign in to comment.