Skip to content

Commit

Permalink
Merge pull request #235 from lsst/tickets/DM-36717
Browse files Browse the repository at this point in the history 
DM-36717: Diffim bug fixes
  • Loading branch information
@isullivan
isullivan committed Nov 4, 2022
2 parents 8a851bd + 82e78da commit fd94035
Show file tree
Hide file tree
Showing 3 changed files with 125 additions and 29 deletions.
118 changes: 100 additions & 18 deletions python/lsst/ip/diffim/subtractImages.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -145,6 +145,19 @@ class AlardLuptonSubtractConfig(lsst.pipe.base.PipelineTaskConfig,
dtype=bool,
default=False,
)
detectionThreshold = lsst.pex.config.Field(
dtype=float,
default=10,
doc="Minimum signal to noise ration of detected sources "
"to use for calculating the PSF matching kernel."
)
badSourceFlags = lsst.pex.config.ListField(
dtype=str,
doc="Flags that, if set, the associated source should not "
"be used to determine the PSF matching kernel.",
default=("sky_source", "slot_Centroid_flag",
"slot_ApFlux_flag", "slot_PsfFlux_flag", ),
)

forceCompatibility = lsst.pex.config.Field(
dtype=bool,
Expand Down Expand Up @@ -263,6 +276,8 @@ def run(self, template, science, sources, finalizedPsfApCorrCatalog=None):
------
RuntimeError
If an unsupported convolution mode is supplied.
RuntimeError
If there are too few sources to calculate the PSF matching kernel.
lsst.pipe.base.NoWorkFound
Raised if fraction of good pixels, defined as not having NO_DATA
set, is less then the configured requiredTemplateFraction
Expand All @@ -280,15 +295,18 @@ def run(self, template, science, sources, finalizedPsfApCorrCatalog=None):
sources = None
sciencePsfSize = getPsfFwhm(science.psf)
templatePsfSize = getPsfFwhm(template.psf)
self.log.info("Science PSF size: %f", sciencePsfSize)
self.log.info("Template PSF size: %f", templatePsfSize)
self.log.info("Science PSF FWHM: %f pixels", sciencePsfSize)
self.log.info("Template PSF FWHM: %f pixels", templatePsfSize)
if self.config.mode == "auto":
if sciencePsfSize < templatePsfSize:
self.log.info("Template PSF size is greater: convolving science image.")
convolveTemplate = False
convolveTemplate = _shapeTest(template.psf, science.psf)
if convolveTemplate:
if sciencePsfSize < templatePsfSize:
self.log.info("Average template PSF size is greater, "
"but science PSF greater in one dimension: convolving template image.")
else:
self.log.info("Science PSF size is greater: convolving template image.")
else:
self.log.info("Science PSF size is greater: convolving template image.")
convolveTemplate = True
self.log.info("Template PSF size is greater: convolving science image.")
elif self.config.mode == "convolveTemplate":
self.log.info("`convolveTemplate` is set: convolving template image.")
convolveTemplate = True
Expand All @@ -297,6 +315,10 @@ def run(self, template, science, sources, finalizedPsfApCorrCatalog=None):
convolveTemplate = False
else:
raise RuntimeError("Cannot handle AlardLuptonSubtract mode: %s", self.config.mode)
# put the template on the same photometric scale as the science image
photoCalib = template.getPhotoCalib()
self.log.info("Applying photometric calibration to template: %f", photoCalib.getCalibrationMean())
template.maskedImage = photoCalib.calibrateImage(template.maskedImage)

if self.config.doScaleVariance and not self.config.forceCompatibility:
# Scale the variance of the template and science images before
Expand All @@ -309,13 +331,17 @@ def run(self, template, science, sources, finalizedPsfApCorrCatalog=None):
self.log.info("Science variance scaling factor: %.2f", sciVarFactor)
self.metadata.add("scaleScienceVarianceFactor", sciVarFactor)

kernelSources = self.makeKernel.selectKernelSources(template, science,
candidateList=sources,
preconvolved=False)
selectSources = self._sourceSelector(sources)
self.log.info("%i sources used out of %i from the input catalog", len(selectSources), len(sources))
if len(selectSources) < self.config.makeKernel.nStarPerCell:
self.log.warning("Too few sources to calculate the PSF matching kernel: "
"%i selected but %i needed for the calculation.",
len(selectSources), self.config.makeKernel.nStarPerCell)
raise RuntimeError("Cannot compute PSF matching kernel: too few sources selected.")
if convolveTemplate:
subtractResults = self.runConvolveTemplate(template, science, kernelSources)
subtractResults = self.runConvolveTemplate(template, science, selectSources)
else:
subtractResults = self.runConvolveScience(template, science, kernelSources)
subtractResults = self.runConvolveScience(template, science, selectSources)

if self.config.doScaleVariance and self.config.forceCompatibility:
# The old behavior scaled the variance of the final image difference.
Expand All @@ -325,7 +351,7 @@ def run(self, template, science, sources, finalizedPsfApCorrCatalog=None):

return subtractResults

def runConvolveTemplate(self, template, science, sources):
def runConvolveTemplate(self, template, science, selectSources):
"""Convolve the template image with a PSF-matching kernel and subtract
from the science image.
Expand All @@ -335,7 +361,7 @@ def runConvolveTemplate(self, template, science, sources):
Template exposure, warped to match the science exposure.
science : `lsst.afw.image.ExposureF`
Science exposure to subtract from the template.
sources : `lsst.afw.table.SourceCatalog`
selectSources : `lsst.afw.table.SourceCatalog`
Identified sources on the science exposure. This catalog is used to
select sources in order to perform the AL PSF matching on stamp
images around them.
Expand All @@ -357,7 +383,12 @@ def runConvolveTemplate(self, template, science, sources):
# Compatibility option to maintain old behavior
# This should be removed in the future!
template = template[science.getBBox()]
kernelResult = self.makeKernel.run(template, science, sources, preconvolved=False)

kernelSources = self.makeKernel.selectKernelSources(template, science,
candidateList=selectSources,
preconvolved=False)
kernelResult = self.makeKernel.run(template, science, kernelSources,
preconvolved=False)

matchedTemplate = self._convolveExposure(template, kernelResult.psfMatchingKernel,
self.convolutionControl,
Expand All @@ -376,7 +407,7 @@ def runConvolveTemplate(self, template, science, sources):
backgroundModel=kernelResult.backgroundModel,
psfMatchingKernel=kernelResult.psfMatchingKernel)

def runConvolveScience(self, template, science, sources):
def runConvolveScience(self, template, science, selectSources):
"""Convolve the science image with a PSF-matching kernel and subtract the template image.
Parameters
Expand All @@ -385,7 +416,7 @@ def runConvolveScience(self, template, science, sources):
Template exposure, warped to match the science exposure.
science : `lsst.afw.image.ExposureF`
Science exposure to subtract from the template.
sources : `lsst.afw.table.SourceCatalog`
selectSources : `lsst.afw.table.SourceCatalog`
Identified sources on the science exposure. This catalog is used to
select sources in order to perform the AL PSF matching on stamp
images around them.
Expand All @@ -408,7 +439,11 @@ def runConvolveScience(self, template, science, sources):
# Compatibility option to maintain old behavior
# This should be removed in the future!
template = template[science.getBBox()]
kernelResult = self.makeKernel.run(science, template, sources, preconvolved=False)
kernelSources = self.makeKernel.selectKernelSources(science, template,
candidateList=selectSources,
preconvolved=False)
kernelResult = self.makeKernel.run(science, template, kernelSources,
preconvolved=False)
modelParams = kernelResult.backgroundModel.getParameters()
# We must invert the background model if the matching kernel is solved for the science image.
kernelResult.backgroundModel.setParameters([-p for p in modelParams])
Expand Down Expand Up @@ -557,6 +592,31 @@ def _convolveExposure(exposure, kernel, convolutionControl,
else:
return convolvedExposure[bbox]

def _sourceSelector(self, sources):
"""Select sources from a catalog that meet the selection criteria.
Parameters
----------
sources : `lsst.afw.table.SourceCatalog`
Input source catalog to select sources from.
Returns
-------
`lsst.afw.table.SourceCatalog`
The source catalog filtered to include only the selected sources.
"""
flags = [True, ]*len(sources)
for flag in self.config.badSourceFlags:
try:
flags *= ~sources[flag]
except Exception as e:
self.log.warning("Could not apply source flag: %s", e)
sToNFlag = (sources.getPsfInstFlux()/sources.getPsfInstFluxErr()) > self.config.detectionThreshold
flags *= sToNFlag
selectSources = sources[flags]

return selectSources.copy(deep=True)


def checkTemplateIsSufficient(templateExposure, logger, requiredTemplateFraction=0.):
"""Raise NoWorkFound if template coverage < requiredTemplateFraction
Expand Down Expand Up @@ -615,3 +675,25 @@ def _subtractImages(science, template, backgroundModel=None):
difference.maskedImage -= backgroundModel
difference.maskedImage -= template.maskedImage
return difference


def _shapeTest(psf1, psf2):
"""Determine whether psf1 is narrower in either dimension than psf2.
Parameters
----------
psf1 : `lsst.afw.detection.Psf`
Reference point spread function (PSF) to evaluate.
psf2 : `lsst.afw.detection.Psf`
Candidate point spread function (PSF) to evaluate.
Returns
-------
`bool`
Returns True if psf1 is narrower than psf2 in either dimension.
"""
shape1 = getPsfFwhm(psf1, average=False)
shape2 = getPsfFwhm(psf2, average=False)
xTest = shape1[0] < shape2[0]
yTest = shape1[1] < shape2[1]
return xTest | yTest
27 changes: 20 additions & 7 deletions python/lsst/ip/diffim/utils.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1080,20 +1080,33 @@ def detectDipoleSources(self, doMerge=True, diffim=None, detectSigma=5.5, grow=3
return detectTask, schema


def getPsfFwhm(psf):
"""Calculate the FWHM in pixels of a supplied PSF.
def getPsfFwhm(psf, average=True):
"""Directly calculate the horizontal and vertical widths
of a PSF at half its maximum value.
Parameters
----------
psf : `lsst.afw.detection.Psf`
Point spread function (PSF) to evaluate.
average : `bool`, optional
Set to return the average width.
Returns
-------
psfSize : `float`
psfSize : `float`, or `tuple` of `float`
The FWHM of the PSF computed at its average position.
Returns the widths along the Y and X axes,
or the average of the two if `average` is set.
"""
sigma2fwhm = 2.*np.sqrt(2.*np.log(2.))
psfAvgPos = psf.getAveragePosition()
psfSize = psf.computeShape(psfAvgPos).getDeterminantRadius()*sigma2fwhm
return psfSize
pos = psf.getAveragePosition()
image = psf.computeKernelImage(pos).array
peak = psf.computePeak(pos)
peakLocs = np.unravel_index(np.argmax(image), image.shape)

def sliceWidth(image, threshold, peaks, axis):
vec = image.take(peaks[1 - axis], axis=axis)
low = np.interp(threshold, vec[:peaks[axis] + 1], np.arange(peaks[axis] + 1))
high = np.interp(threshold, vec[:peaks[axis] - 1:-1], np.arange(len(vec) - 1, peaks[axis] - 1, -1))
return high - low
width = (sliceWidth(image, peak/2., peakLocs, axis=0), sliceWidth(image, peak/2., peakLocs, axis=1))
return np.mean(width) if average else width
9 changes: 5 additions & 4 deletions tests/test_subtractTask.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -365,12 +365,13 @@ def test_few_sources(self):
"""
xSize = 256
ySize = 256
science, sources = makeTestImage(psfSize=2.4, nSrc=1, xSize=xSize, ySize=ySize)
template, _ = makeTestImage(psfSize=2.0, nSrc=1, xSize=xSize, ySize=ySize, doApplyCalibration=True)
science, sources = makeTestImage(psfSize=2.4, nSrc=10, xSize=xSize, ySize=ySize)
template, _ = makeTestImage(psfSize=2.0, nSrc=10, xSize=xSize, ySize=ySize, doApplyCalibration=True)
config = subtractImages.AlardLuptonSubtractTask.ConfigClass()
task = subtractImages.AlardLuptonSubtractTask(config=config)
with self.assertRaisesRegex(lsst.pex.exceptions.Exception,
'Unable to determine kernel sum; 0 candidates'):
sources = sources[0:1]
with self.assertRaisesRegex(RuntimeError,
"Cannot compute PSF matching kernel: too few sources selected."):
task.run(template, science, sources)

def test_order_equal_images(self):
Expand Down

0 comments on commit fd94035

Please sign in to comment.