Coverage for python/lsst/ip/diffim/subtractImages.py: 87%
498 statements
« prev ^ index » next coverage.py v7.16.0, created at 2026-09-17 09:18 +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/>.
22from astropy import units as u
23from astropy.stats import gaussian_fwhm_to_sigma
24import numpy as np
26import lsst.afw.detection as afwDetection
27import lsst.afw.image
28import lsst.afw.math
29import lsst.geom
30from lsst.ip.diffim.utils import (evaluateMeanPsfFwhm, getPsfFwhm,
31 computeDifferenceImageMetrics,
32 checkMask, setSourceFootprints)
33from lsst.meas.algorithms import ScaleVarianceTask, ScienceSourceSelectorTask
34import lsst.pex.config
35import lsst.pipe.base
36import lsst.pex.exceptions
37from lsst.pipe.base import connectionTypes
38from . import MakeKernelTask, DecorrelateALKernelTask
39from lsst.utils.timer import timeMethod
41__all__ = ["AlardLuptonSubtractConfig", "AlardLuptonSubtractTask",
42 "AlardLuptonPreconvolveSubtractConfig", "AlardLuptonPreconvolveSubtractTask",
43 "SimplifiedSubtractConfig", "SimplifiedSubtractTask",
44 "InsufficientKernelSourcesError"]
46_dimensions = ("instrument", "visit", "detector")
47_defaultTemplates = {"coaddName": "deep", "fakesType": ""}
50class InsufficientKernelSourcesError(lsst.pipe.base.AlgorithmError):
51 """Raised when there are too few sources to calculate the PSF matching
52 kernel.
53 """
54 def __init__(self, *, nSources, nRequired):
55 msg = (f"Only {nSources} sources were selected for PSF matching,"
56 f" but {nRequired} are required.")
57 super().__init__(msg)
58 self.nSources = nSources
59 self.nRequired = nRequired
61 @property
62 def metadata(self):
63 return {"nSources": self.nSources,
64 "nRequired": self.nRequired
65 }
68class SubtractInputConnections(lsst.pipe.base.PipelineTaskConnections,
69 dimensions=_dimensions,
70 defaultTemplates=_defaultTemplates):
71 template = connectionTypes.Input(
72 doc="Input warped template to subtract.",
73 dimensions=("instrument", "visit", "detector"),
74 storageClass="ExposureF",
75 name="{fakesType}{coaddName}Diff_templateExp"
76 )
77 science = connectionTypes.Input(
78 doc="Input science exposure to subtract from.",
79 dimensions=("instrument", "visit", "detector"),
80 storageClass="ExposureF",
81 name="{fakesType}calexp"
82 )
83 sources = connectionTypes.Input(
84 doc="Sources measured on the science exposure; "
85 "used to select sources for making the matching kernel.",
86 dimensions=("instrument", "visit", "detector"),
87 storageClass="SourceCatalog",
88 name="{fakesType}src"
89 )
90 visitSummary = connectionTypes.Input(
91 doc=("Per-visit catalog with final calibration objects. "
92 "These catalogs use the detector id for the catalog id, "
93 "sorted on id for fast lookup."),
94 dimensions=("instrument", "visit"),
95 storageClass="ExposureCatalog",
96 name="finalVisitSummary",
97 )
99 def __init__(self, *, config=None):
100 super().__init__(config=config)
101 if not config.doApplyExternalCalibrations:
102 del self.visitSummary
105class SubtractImageOutputConnections(lsst.pipe.base.PipelineTaskConnections,
106 dimensions=_dimensions,
107 defaultTemplates=_defaultTemplates):
108 difference = connectionTypes.Output(
109 doc="Result of subtracting convolved template from science image.",
110 dimensions=("instrument", "visit", "detector"),
111 storageClass="ExposureF",
112 name="{fakesType}{coaddName}Diff_differenceTempExp",
113 )
114 matchedTemplate = connectionTypes.Output(
115 doc="Warped and PSF-matched template used to create `subtractedExposure`.",
116 dimensions=("instrument", "visit", "detector"),
117 storageClass="ExposureF",
118 name="{fakesType}{coaddName}Diff_matchedExp",
119 )
120 psfMatchingKernel = connectionTypes.Output(
121 doc="Kernel used to PSF match the science and template images.",
122 dimensions=("instrument", "visit", "detector"),
123 storageClass="MatchingKernel",
124 name="{fakesType}{coaddName}Diff_psfMatchKernel",
125 )
126 kernelSources = connectionTypes.Output(
127 doc="Final selection of sources used for psf matching.",
128 dimensions=("instrument", "visit", "detector"),
129 storageClass="SourceCatalog",
130 name="{fakesType}{coaddName}Diff_psfMatchSources"
131 )
134class SubtractScoreOutputConnections(lsst.pipe.base.PipelineTaskConnections,
135 dimensions=_dimensions,
136 defaultTemplates=_defaultTemplates):
137 scoreExposure = connectionTypes.Output(
138 doc="The maximum likelihood image, used for the detection of diaSources.",
139 dimensions=("instrument", "visit", "detector"),
140 storageClass="ExposureF",
141 name="{fakesType}{coaddName}Diff_scoreTempExp",
142 )
143 psfMatchingKernel = connectionTypes.Output(
144 doc="Kernel used to PSF match the science and template images.",
145 dimensions=("instrument", "visit", "detector"),
146 storageClass="MatchingKernel",
147 name="{fakesType}{coaddName}Diff_psfScoreMatchKernel",
148 )
149 kernelSources = connectionTypes.Output(
150 doc="Final selection of sources used for psf matching.",
151 dimensions=("instrument", "visit", "detector"),
152 storageClass="SourceCatalog",
153 name="{fakesType}{coaddName}Diff_psfScoreMatchSources"
154 )
157class AlardLuptonSubtractConnections(SubtractInputConnections, SubtractImageOutputConnections):
158 pass
161class SimplifiedSubtractConnections(SubtractInputConnections, SubtractImageOutputConnections):
162 inputPsfMatchingKernel = connectionTypes.Input(
163 doc="Kernel used to PSF match the science and template images.",
164 dimensions=("instrument", "visit", "detector"),
165 storageClass="MatchingKernel",
166 name="{fakesType}{coaddName}Diff_psfMatchKernel",
167 )
169 def __init__(self, *, config=None):
170 super().__init__(config=config)
171 del self.sources
172 if config.useExistingKernel:
173 del self.psfMatchingKernel
174 del self.kernelSources
175 else:
176 del self.inputPsfMatchingKernel
179class AlardLuptonSubtractBaseConfig(lsst.pex.config.Config):
180 makeKernel = lsst.pex.config.ConfigurableField(
181 target=MakeKernelTask,
182 doc="Task to construct a matching kernel for convolution.",
183 )
184 doDecorrelation = lsst.pex.config.Field(
185 dtype=bool,
186 default=True,
187 doc="Perform diffim decorrelation to undo pixel correlation due to A&L "
188 "kernel convolution? If True, also update the diffim PSF."
189 )
190 decorrelate = lsst.pex.config.ConfigurableField(
191 target=DecorrelateALKernelTask,
192 doc="Task to decorrelate the image difference.",
193 )
194 requiredTemplateFraction = lsst.pex.config.Field(
195 dtype=float,
196 default=0.1,
197 doc="Raise NoWorkFound and do not attempt image subtraction if template covers less than this "
198 " fraction of pixels. Setting to 0 will always attempt image subtraction."
199 )
200 minTemplateFractionForExpectedSuccess = lsst.pex.config.Field(
201 dtype=float,
202 default=0.2,
203 doc="Raise NoWorkFound if PSF-matching fails and template covers less than this fraction of pixels."
204 " If the fraction of pixels covered by the template is less than this value (and greater than"
205 " requiredTemplateFraction) this task is attempted but failure is anticipated and tolerated."
206 )
207 doScaleVariance = lsst.pex.config.Field(
208 dtype=bool,
209 default=True,
210 doc="Scale variance of the science image? Note that the template variance is NOT scaled"
211 " here. The template variance may be scaled independently in ``GetTemplateTask``."
212 )
213 scaleVariance = lsst.pex.config.ConfigurableField(
214 target=ScaleVarianceTask,
215 doc="Subtask to rescale the variance of the template to the statistically expected level."
216 )
217 doSubtractBackground = lsst.pex.config.Field(
218 doc="Subtract the background fit when solving the kernel? "
219 "It is generally better to instead subtract the background in detectAndMeasure.",
220 dtype=bool,
221 default=False,
222 )
223 doApplyExternalCalibrations = lsst.pex.config.Field(
224 doc=(
225 "Replace science Exposure's calibration objects with those"
226 " in visitSummary. Ignored if `doApplyFinalizedPsf is True."
227 ),
228 dtype=bool,
229 default=False,
230 )
231 sourceSelector = lsst.pex.config.ConfigurableField(
232 target=ScienceSourceSelectorTask,
233 doc="Task to select sources to be used for PSF matching.",
234 )
235 fallbackSourceSelector = lsst.pex.config.ConfigurableField(
236 target=ScienceSourceSelectorTask,
237 doc="Task to select sources to be used for PSF matching."
238 "Used only if the kernel calculation fails and"
239 "`allowKernelSourceDetection` is set. The fallback source detection"
240 " will not include all of the same plugins as the original source "
241 " detection, so not all of the same flags can be used.",
242 )
243 detectionThreshold = lsst.pex.config.Field(
244 dtype=float,
245 default=10,
246 doc="Minimum signal to noise ratio of detected sources "
247 "to use for calculating the PSF matching kernel.",
248 deprecated="No longer used. Will be removed after v30"
249 )
250 detectionThresholdMax = lsst.pex.config.Field(
251 dtype=float,
252 default=500,
253 doc="Maximum signal to noise ratio of detected sources "
254 "to use for calculating the PSF matching kernel.",
255 deprecated="No longer used. Will be removed after v30"
256 )
257 restrictKernelEdgeSources = lsst.pex.config.Field(
258 dtype=bool,
259 default=True,
260 doc="Exclude sources close to the edge from the kernel calculation?"
261 )
262 maxKernelSources = lsst.pex.config.Field(
263 dtype=int,
264 default=1000,
265 doc="Maximum number of sources to use for calculating the PSF matching kernel."
266 "Set to -1 to disable."
267 )
268 minKernelSources = lsst.pex.config.Field(
269 dtype=int,
270 default=3,
271 doc="Minimum number of sources needed for calculating the PSF matching kernel."
272 )
273 excludeMaskPlanes = lsst.pex.config.ListField(
274 dtype=str,
275 default=("NO_DATA", "BAD", "SAT", "EDGE", "FAKE", "HIGH_VARIANCE"),
276 doc="Template mask planes to exclude when selecting sources for PSF matching.",
277 )
278 badMaskPlanes = lsst.pex.config.ListField(
279 dtype=str,
280 default=("NO_DATA", "BAD", "SAT", "EDGE"),
281 doc="Mask planes to interpolate over."
282 )
283 preserveTemplateMask = lsst.pex.config.ListField(
284 dtype=str,
285 default=("NO_DATA", "BAD", "HIGH_VARIANCE"),
286 doc="Mask planes from the template to propagate to the image difference."
287 )
288 renameTemplateMask = lsst.pex.config.ListField(
289 dtype=str,
290 default=("SAT", "INJECTED", "INJECTED_CORE",),
291 doc="Mask planes from the template to propagate to the image difference"
292 "with '_TEMPLATE' appended to the name."
293 )
294 preserveMaskPlanes = lsst.pex.config.ListField(
295 dtype=str,
296 default=("INJECTED", "INJECTED_CORE", "INJECTED_TEMPLATE", "INJECTED_CORE_TEMPLATE"),
297 doc="Mask planes to preserve without dilation when convolving the image.",
298 )
299 allowKernelSourceDetection = lsst.pex.config.Field(
300 dtype=bool,
301 default=False,
302 doc="Re-run source detection for kernel candidates if an error is"
303 " encountered while calculating the matching kernel."
304 )
306 def setDefaults(self):
307 self.makeKernel.kernel.name = "AL"
308 # Always include background fitting in the kernel fit,
309 # even if it is not subtracted
310 self.makeKernel.kernel.active.fitForBackground = True
311 self.makeKernel.kernel.active.spatialKernelOrder = 1
312 self.makeKernel.kernel.active.spatialBgOrder = 2
313 # Shared source selector settings
314 doSkySources = False # Do not include sky sources
315 doSignalToNoise = True # apply signal to noise filter
316 doUnresolved = True # apply star-galaxy separation
317 signalToNoiseMinimum = 10
318 signalToNoiseMaximum = 500
319 self.sourceSelector.doIsolated = True # apply isolated star selection
320 self.sourceSelector.doRequirePrimary = True # apply primary flag selection
321 self.sourceSelector.doUnresolved = doUnresolved
322 self.sourceSelector.doSkySources = doSkySources
323 self.sourceSelector.doSignalToNoise = doSignalToNoise
324 self.sourceSelector.signalToNoise.minimum = signalToNoiseMinimum
325 self.sourceSelector.signalToNoise.maximum = signalToNoiseMaximum
326 # The following two configs should not be necessary to be turned on for
327 # PSF-matching, and the fallback kernel source selection will fail if
328 # they are set since it does not run deblending.
329 self.fallbackSourceSelector.doIsolated = False # Do not apply isolated star selection
330 self.fallbackSourceSelector.doRequirePrimary = False # Do not apply primary flag selection
331 self.fallbackSourceSelector.doUnresolved = doUnresolved
332 self.fallbackSourceSelector.doSkySources = doSkySources
333 self.fallbackSourceSelector.doSignalToNoise = doSignalToNoise
334 self.fallbackSourceSelector.signalToNoise.minimum = signalToNoiseMinimum
335 self.fallbackSourceSelector.signalToNoise.maximum = signalToNoiseMaximum
338class AlardLuptonSubtractConfig(AlardLuptonSubtractBaseConfig, lsst.pipe.base.PipelineTaskConfig,
339 pipelineConnections=AlardLuptonSubtractConnections):
340 mode = lsst.pex.config.ChoiceField(
341 dtype=str,
342 default="convolveTemplate",
343 allowed={"auto": "Choose which image to convolve at runtime.",
344 "convolveScience": "Only convolve the science image.",
345 "convolveTemplate": "Only convolve the template image."},
346 doc="Choose which image to convolve at runtime, or require that a specific image is convolved."
347 )
350class AlardLuptonSubtractTask(lsst.pipe.base.PipelineTask):
351 """Compute the image difference of a science and template image using
352 the Alard & Lupton (1998) algorithm.
353 """
354 ConfigClass = AlardLuptonSubtractConfig
355 _DefaultName = "alardLuptonSubtract"
356 usePreconvolution = False
357 """Whether this task preconvolves the science image with its own PSF
358 before kernel-matching. Subclasses that preconvolve override this to
359 `True`."""
361 def __init__(self, **kwargs):
362 super().__init__(**kwargs)
363 self.makeSubtask("decorrelate")
364 self.makeSubtask("makeKernel")
365 self.makeSubtask("sourceSelector")
366 self.makeSubtask("fallbackSourceSelector")
367 if self.config.doScaleVariance:
368 self.makeSubtask("scaleVariance")
370 self.convolutionControl = lsst.afw.math.ConvolutionControl()
371 # Normalization is an extra, unnecessary, calculation and will result
372 # in mis-subtraction of the images if there are calibration errors.
373 self.convolutionControl.setDoNormalize(False)
374 self.convolutionControl.setDoCopyEdge(True)
376 def _applyExternalCalibrations(self, exposure, visitSummary):
377 """Replace calibrations (psf, and ApCorrMap) on this exposure with
378 external ones.".
380 Parameters
381 ----------
382 exposure : `lsst.afw.image.exposure.Exposure`
383 Input exposure to adjust calibrations.
384 visitSummary : `lsst.afw.table.ExposureCatalog`
385 Exposure catalog with external calibrations to be applied. Catalog
386 uses the detector id for the catalog id, sorted on id for fast
387 lookup.
389 Returns
390 -------
391 exposure : `lsst.afw.image.exposure.Exposure`
392 Exposure with adjusted calibrations.
393 """
394 detectorId = exposure.info.getDetector().getId()
396 row = visitSummary.find(detectorId)
397 if row is None:
398 self.log.warning("Detector id %s not found in external calibrations catalog; "
399 "Using original calibrations.", detectorId)
400 else:
401 psf = row.getPsf()
402 apCorrMap = row.getApCorrMap()
403 if psf is None:
404 self.log.warning("Detector id %s has None for psf in "
405 "external calibrations catalog; Using original psf and aperture correction.",
406 detectorId)
407 elif apCorrMap is None:
408 self.log.warning("Detector id %s has None for apCorrMap in "
409 "external calibrations catalog; Using original psf and aperture correction.",
410 detectorId)
411 else:
412 exposure.setPsf(psf)
413 exposure.info.setApCorrMap(apCorrMap)
415 return exposure
417 def runQuantum(self, butlerQC, inputRefs, outputRefs):
418 inputs = butlerQC.get(inputRefs)
420 try:
421 results = self.run(**inputs)
422 except lsst.pipe.base.AlgorithmError as e:
423 error = lsst.pipe.base.AnnotatedPartialOutputsError.annotate(e, self, log=self.log)
424 # No partial outputs for butler to put
425 raise error from e
427 butlerQC.put(results, outputRefs)
429 @timeMethod
430 def run(self, template, science, sources, visitSummary=None):
431 """PSF match, subtract, and decorrelate two images.
433 Parameters
434 ----------
435 template : `lsst.afw.image.ExposureF`
436 Template exposure, warped to match the science exposure.
437 science : `lsst.afw.image.ExposureF`
438 Science exposure to subtract from the template.
439 sources : `lsst.afw.table.SourceCatalog`
440 Identified sources on the science exposure. This catalog is used to
441 select sources in order to perform the AL PSF matching on stamp
442 images around them.
443 visitSummary : `lsst.afw.table.ExposureCatalog`, optional
444 Exposure catalog with external calibrations to be applied. Catalog
445 uses the detector id for the catalog id, sorted on id for fast
446 lookup.
448 Returns
449 -------
450 results : `lsst.pipe.base.Struct`
451 ``difference`` : `lsst.afw.image.ExposureF`
452 Result of subtracting template and science.
453 ``matchedTemplate`` : `lsst.afw.image.ExposureF`
454 Warped and PSF-matched template exposure.
455 ``backgroundModel`` : `lsst.afw.math.Function2D`
456 Background model that was fit while solving for the
457 PSF-matching kernel
458 ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
459 Kernel used to PSF-match the convolved image.
460 ``kernelSources` : `lsst.afw.table.SourceCatalog`
461 Sources from the input catalog that were used to construct the
462 PSF-matching kernel.
463 """
464 self._prepareInputs(template, science, visitSummary=visitSummary)
466 convolveTemplate = self.chooseConvolutionMethod(template, science)
467 self.matchedPsfSize = self.sciencePsfSize if convolveTemplate else self.templatePsfSize
469 kernelResult = self.runMakeKernel(template, science, sources=sources,
470 convolveTemplate=convolveTemplate,
471 runSourceDetection=False)
473 if self.config.doSubtractBackground:
474 backgroundModel = kernelResult.backgroundModel
475 else:
476 backgroundModel = None
477 if convolveTemplate:
478 subtractResults = self.runConvolveTemplate(template, science, kernelResult.psfMatchingKernel,
479 backgroundModel=backgroundModel)
480 else:
481 subtractResults = self.runConvolveScience(template, science, kernelResult.psfMatchingKernel,
482 backgroundModel=backgroundModel)
483 subtractResults.kernelSources = kernelResult.kernelSources
485 metrics = computeDifferenceImageMetrics(science, subtractResults.difference, sources)
487 self.metadata["differenceFootprintRatioMean"] = metrics.differenceFootprintRatioMean
488 self.metadata["differenceFootprintRatioStdev"] = metrics.differenceFootprintRatioStdev
489 self.metadata["differenceFootprintSkyRatioMean"] = metrics.differenceFootprintSkyRatioMean
490 self.metadata["differenceFootprintSkyRatioStdev"] = metrics.differenceFootprintSkyRatioStdev
491 self.log.info("Mean, stdev of ratio of difference to science "
492 "pixels in star footprints: %5.4f, %5.4f",
493 self.metadata["differenceFootprintRatioMean"],
494 self.metadata["differenceFootprintRatioStdev"])
496 return subtractResults
498 def chooseConvolutionMethod(self, template, science):
499 """Determine whether the template should be convolved with the PSF
500 matching kernel.
502 Parameters
503 ----------
504 template : `lsst.afw.image.ExposureF`
505 Template exposure, warped to match the science exposure.
506 science : `lsst.afw.image.ExposureF`
507 Science exposure to subtract from the template.
509 Returns
510 -------
511 convolveTemplate : `bool`
512 Convolve the template to match the two images?
514 Raises
515 ------
516 RuntimeError
517 If an unsupported convolution mode is supplied.
518 """
519 if self.usePreconvolution: 519 ↛ 520line 519 didn't jump to line 520 because the condition on line 519 was never true
520 raise RuntimeError("Choosing a convolution method is incompatible with preconvolution!")
521 if self.config.mode == "auto":
522 convolveTemplate = _shapeTest(template,
523 science,
524 fwhmExposureBuffer=self.config.makeKernel.fwhmExposureBuffer,
525 fwhmExposureGrid=self.config.makeKernel.fwhmExposureGrid)
526 if convolveTemplate:
527 if self.sciencePsfSize < self.templatePsfSize: 527 ↛ 528line 527 didn't jump to line 528 because the condition on line 527 was never true
528 self.log.info("Average template PSF size is greater, "
529 "but science PSF greater in one dimension: convolving template image.")
530 else:
531 self.log.info("Science PSF size is greater: convolving template image.")
532 else:
533 self.log.info("Template PSF size is greater: convolving science image.")
534 elif self.config.mode == "convolveTemplate":
535 self.log.info("`convolveTemplate` is set: convolving template image.")
536 convolveTemplate = True
537 elif self.config.mode == "convolveScience": 537 ↛ 541line 537 didn't jump to line 541 because the condition on line 537 was always true
538 self.log.info("`convolveScience` is set: convolving science image.")
539 convolveTemplate = False
540 else:
541 raise RuntimeError(f"Cannot handle AlardLuptonSubtract mode: {self.config.mode}")
542 return convolveTemplate
544 def runMakeKernel(self, template, science, sources=None, convolveTemplate=True, runSourceDetection=False):
545 """Construct the PSF-matching kernel. Not used for preconvolution.
547 Parameters
548 ----------
549 template : `lsst.afw.image.ExposureF`
550 Template exposure, warped to match the science exposure.
551 science : `lsst.afw.image.ExposureF`
552 Science exposure to subtract from the template.
553 sources : `lsst.afw.table.SourceCatalog`
554 Identified sources on the science exposure. This catalog is used to
555 select sources in order to perform the AL PSF matching on stamp
556 images around them.
557 Not used if ``runSourceDetection`` is set.
558 convolveTemplate : `bool`, optional
559 Construct the matching kernel to convolve the template?
560 runSourceDetection : `bool`, optional
561 Run a minimal version of source detection to determine kernel
562 candidates? If False, a source list to select kernel candidates
563 from must be supplied.
565 Returns
566 -------
567 results : `lsst.pipe.base.Struct`
568 ``backgroundModel`` : `lsst.afw.math.Function2D`
569 Background model that was fit while solving for the
570 PSF-matching kernel
571 ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
572 Kernel used to PSF-match the convolved image.
573 ``kernelSources` : `lsst.afw.table.SourceCatalog`
574 Sources from the input catalog that were used to construct the
575 PSF-matching kernel.
576 """
577 if self.usePreconvolution: 577 ↛ 578line 577 didn't jump to line 578 because the condition on line 577 was never true
578 raise RuntimeError("Incorrect matching kernel calculation configured. "
579 "`runMakeKernel` can't be called if `usePreconvolution` is set.")
580 if convolveTemplate:
581 reference = template
582 target = science
583 referenceFwhmPix = self.templatePsfSize
584 targetFwhmPix = self.sciencePsfSize
585 else:
586 reference = science
587 target = template
588 referenceFwhmPix = self.sciencePsfSize
589 targetFwhmPix = self.templatePsfSize
590 try:
591 if runSourceDetection:
592 kernelSources = self.runKernelSourceDetection(template, science)
593 else:
594 kernelSources = self._sourceSelector(template, science, sources)
595 kernelResult = self.makeKernel.run(reference, target, kernelSources,
596 preconvolved=False,
597 templateFwhmPix=referenceFwhmPix,
598 scienceFwhmPix=targetFwhmPix)
599 except (RuntimeError, lsst.pex.exceptions.Exception) as e:
600 self.log.warning("Failed to match template. Checking coverage")
601 # Raise NoWorkFound if template fraction is insufficient
602 checkTemplateIsSufficient(template[science.getBBox()], science, self.log,
603 self.config.minTemplateFractionForExpectedSuccess,
604 exceptionMessage="Template coverage lower than expected to succeed."
605 f" Failure is tolerable: {e}")
606 # checkTemplateIsSufficient did not raise NoWorkFound, so raise original exception
607 raise e
609 return lsst.pipe.base.Struct(backgroundModel=kernelResult.backgroundModel,
610 psfMatchingKernel=kernelResult.psfMatchingKernel,
611 kernelSources=kernelSources)
613 def runKernelSourceDetection(self, template, science):
614 """Run detection on the science image and use the template mask plane
615 to reject candidate sources.
617 Parameters
618 ----------
619 template : `lsst.afw.image.ExposureF`
620 Template exposure, warped to match the science exposure.
621 science : `lsst.afw.image.ExposureF`
622 Science exposure to subtract from the template.
624 Returns
625 -------
626 kernelSources : `lsst.afw.table.SourceCatalog`
627 Sources from the input catalog to use to construct the
628 PSF-matching kernel.
629 """
630 kernelSize = self.makeKernel.makeKernelBasisList(
631 self.templatePsfSize, self.matchedPsfSize)[0].getWidth()
632 sources = self.makeKernel.makeCandidateList(template, science, kernelSize,
633 candidateList=None,
634 sigma=gaussian_fwhm_to_sigma*self.sciencePsfSize)
635 return self._sourceSelector(template, science, sources, fallback=True)
637 def runConvolveTemplate(self, template, science, psfMatchingKernel, backgroundModel=None):
638 """Convolve the template image with a PSF-matching kernel and subtract
639 from the science image.
641 Parameters
642 ----------
643 template : `lsst.afw.image.ExposureF`
644 Template exposure, warped to match the science exposure.
645 science : `lsst.afw.image.ExposureF`
646 Science exposure to subtract from the template.
647 psfMatchingKernel : `lsst.afw.math.Kernel`
648 Kernel to be used to PSF-match the science image to the template.
649 backgroundModel : `lsst.afw.math.Function2D`, optional
650 Background model that was fit while solving for the PSF-matching
651 kernel.
653 Returns
654 -------
655 results : `lsst.pipe.base.Struct`
657 ``difference`` : `lsst.afw.image.ExposureF`
658 Result of subtracting template and science.
659 ``matchedTemplate`` : `lsst.afw.image.ExposureF`
660 Warped and PSF-matched template exposure.
661 ``backgroundModel`` : `lsst.afw.math.Function2D`
662 Background model that was fit while solving for the PSF-matching kernel
663 ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
664 Kernel used to PSF-match the template to the science image.
665 """
666 self.metadata["convolvedExposure"] = "Template"
668 matchedTemplate = self._convolveExposure(template, psfMatchingKernel,
669 self.convolutionControl,
670 bbox=science.getBBox(),
671 psf=science.psf,
672 photoCalib=science.photoCalib)
674 difference = _subtractImages(science, matchedTemplate, backgroundModel=backgroundModel)
675 correctedExposure = self.finalize(template, science, difference,
676 psfMatchingKernel,
677 templateMatched=True)
679 return lsst.pipe.base.Struct(difference=correctedExposure,
680 matchedTemplate=matchedTemplate,
681 matchedScience=science,
682 backgroundModel=backgroundModel,
683 psfMatchingKernel=psfMatchingKernel)
685 def runConvolveScience(self, template, science, psfMatchingKernel, backgroundModel=None):
686 """Convolve the science image with a PSF-matching kernel and subtract
687 the template image.
689 Parameters
690 ----------
691 template : `lsst.afw.image.ExposureF`
692 Template exposure, warped to match the science exposure.
693 science : `lsst.afw.image.ExposureF`
694 Science exposure to subtract from the template.
695 psfMatchingKernel : `lsst.afw.math.Kernel`
696 Kernel to be used to PSF-match the science image to the template.
697 backgroundModel : `lsst.afw.math.Function2D`, optional
698 Background model that was fit while solving for the PSF-matching
699 kernel.
701 Returns
702 -------
703 results : `lsst.pipe.base.Struct`
705 ``difference`` : `lsst.afw.image.ExposureF`
706 Result of subtracting template and science.
707 ``matchedTemplate`` : `lsst.afw.image.ExposureF`
708 Warped template exposure. Note that in this case, the template
709 is not PSF-matched to the science image.
710 ``backgroundModel`` : `lsst.afw.math.Function2D`
711 Background model that was fit while solving for the PSF-matching kernel
712 ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
713 Kernel used to PSF-match the science image to the template.
714 """
715 self.metadata["convolvedExposure"] = "Science"
716 bbox = science.getBBox()
718 kernelImage = lsst.afw.image.ImageD(psfMatchingKernel.getDimensions())
719 xcen, ycen = bbox.getCenter()
720 norm = psfMatchingKernel.computeImage(kernelImage, doNormalize=False, x=xcen, y=ycen)
722 matchedScience = self._convolveExposure(science, psfMatchingKernel,
723 self.convolutionControl,
724 psf=template.psf)
726 # Place back on native photometric scale
727 matchedScience.maskedImage /= norm
728 matchedTemplate = template.clone()[bbox]
729 matchedTemplate.setPhotoCalib(science.photoCalib)
731 if backgroundModel is not None:
732 # We must invert the background model if the matching kernel is solved for the science image.
733 invertedBackground = backgroundModel.clone()
734 invertedBackground.setParameters([-p for p in backgroundModel.getParameters()])
735 backgroundModel = invertedBackground
737 difference = _subtractImages(matchedScience, matchedTemplate, backgroundModel=backgroundModel)
739 correctedExposure = self.finalize(template, science, difference,
740 psfMatchingKernel,
741 templateMatched=False)
743 return lsst.pipe.base.Struct(difference=correctedExposure,
744 matchedTemplate=matchedTemplate,
745 matchedScience=matchedScience,
746 backgroundModel=backgroundModel,
747 psfMatchingKernel=psfMatchingKernel)
749 def finalize(self, template, science, difference, kernel,
750 templateMatched=True,
751 preConvMode=False,
752 preConvKernel=None,
753 spatiallyVarying=False):
754 """Decorrelate the difference image to undo the noise correlations
755 caused by convolution.
757 Parameters
758 ----------
759 template : `lsst.afw.image.ExposureF`
760 Template exposure, warped to match the science exposure.
761 science : `lsst.afw.image.ExposureF`
762 Science exposure to subtract from the template.
763 difference : `lsst.afw.image.ExposureF`
764 Result of subtracting template and science.
765 kernel : `lsst.afw.math.Kernel`
766 An (optionally spatially-varying) PSF matching kernel
767 templateMatched : `bool`, optional
768 Was the template PSF-matched to the science image?
769 preConvMode : `bool`, optional
770 Was the science image preconvolved with its own PSF
771 before PSF matching the template?
772 preConvKernel : `lsst.afw.detection.Psf`, optional
773 If not `None`, then the science image was pre-convolved with
774 (the reflection of) this kernel. Must be normalized to sum to 1.
775 spatiallyVarying : `bool`, optional
776 Compute the decorrelation kernel spatially varying across the image?
778 Returns
779 -------
780 correctedExposure : `lsst.afw.image.ExposureF`
781 The decorrelated image difference.
782 """
783 if self.config.doDecorrelation:
784 self.log.info("Decorrelating image difference.")
785 # We have cleared the template mask plane, so copy the mask plane of
786 # the image difference so that we can calculate correct statistics
787 # during decorrelation
788 correctedExposure = self.decorrelate.run(science, template[science.getBBox()], difference, kernel,
789 templateMatched=templateMatched,
790 preConvMode=preConvMode,
791 preConvKernel=preConvKernel,
792 spatiallyVarying=spatiallyVarying).correctedExposure
793 else:
794 self.log.info("NOT decorrelating image difference.")
795 correctedExposure = difference
796 return correctedExposure
798 def _calculateMagLim(self, exposure, nsigma=5.0, fallbackPsfSize=None):
799 """Calculate an exposure's limiting magnitude.
801 This method uses the photometric zeropoint together with the
802 PSF size from the average position of the exposure.
804 Parameters
805 ----------
806 exposure : `lsst.afw.image.Exposure`
807 The target exposure to calculate the limiting magnitude for.
808 nsigma : `float`, optional
809 The detection threshold in sigma.
810 fallbackPsfSize : `float`, optional
811 PSF FWHM to use in the event the exposure PSF cannot be retrieved.
813 Returns
814 -------
815 maglim : `astropy.units.Quantity`
816 The limiting magnitude of the exposure, or np.nan.
817 """
818 if exposure.photoCalib is None: 818 ↛ 819line 818 didn't jump to line 819 because the condition on line 818 was never true
819 return np.nan
820 try:
821 psf = exposure.getPsf()
822 psf_shape = psf.computeShape(psf.getAveragePosition())
823 except (lsst.pex.exceptions.InvalidParameterError,
824 afwDetection.InvalidPsfError,
825 lsst.pex.exceptions.RangeError):
826 if fallbackPsfSize is None:
827 self.log.info("Unable to evaluate PSF, setting maglim to nan")
828 return np.nan
829 self.log.info("Unable to evaluate PSF, using fallback FWHM %f", fallbackPsfSize)
830 psf_area = np.pi*(fallbackPsfSize/2)**2
831 else:
832 # Get a more accurate area than `psf_shape.getArea()` via moments
833 psf_area = np.pi*np.sqrt(psf_shape.getIxx()*psf_shape.getIyy())
835 zeropoint = exposure.photoCalib.instFluxToMagnitude(1)
836 return zeropoint - 2.5*np.log10(nsigma*np.sqrt(psf_area))
838 @staticmethod
839 def _validateExposures(template, science):
840 """Check that the WCS of the two Exposures match, the template bbox
841 contains the science bbox, and that the bands match.
843 Parameters
844 ----------
845 template : `lsst.afw.image.ExposureF`
846 Template exposure, warped to match the science exposure.
847 science : `lsst.afw.image.ExposureF`
848 Science exposure to subtract from the template.
850 Raises
851 ------
852 AssertionError
853 Raised if the WCS of the template is not equal to the science WCS,
854 if the science image is not fully contained in the template
855 bounding box, or if the bands do not match.
856 """
857 assert template.wcs == science.wcs, \
858 "Template and science exposure WCS are not identical."
859 templateBBox = template.getBBox()
860 scienceBBox = science.getBBox()
861 assert science.filter.bandLabel == template.filter.bandLabel, \
862 "Science and template exposures have different bands: %s, %s" % \
863 (science.filter, template.filter)
865 assert templateBBox.contains(scienceBBox), \
866 "Template bbox does not contain all of the science image."
868 def _convolveExposure(self, exposure, kernel, convolutionControl,
869 bbox=None,
870 psf=None,
871 photoCalib=None,
872 interpolateBadMaskPlanes=False,
873 ):
874 """Convolve an exposure with the given kernel.
876 Parameters
877 ----------
878 exposure : `lsst.afw.Exposure`
879 exposure to convolve.
880 kernel : `lsst.afw.math.LinearCombinationKernel`
881 PSF matching kernel computed in the ``makeKernel`` subtask.
882 convolutionControl : `lsst.afw.math.ConvolutionControl`
883 Configuration for convolve algorithm.
884 bbox : `lsst.geom.Box2I`, optional
885 Bounding box to trim the convolved exposure to.
886 psf : `lsst.afw.detection.Psf`, optional
887 Point spread function (PSF) to set for the convolved exposure.
888 photoCalib : `lsst.afw.image.PhotoCalib`, optional
889 Photometric calibration of the convolved exposure.
890 interpolateBadMaskPlanes : `bool`, optional
891 If set, interpolate over mask planes specified in
892 ``config.badMaskPlanes`` before convolving the image.
894 Returns
895 -------
896 convolvedExp : `lsst.afw.Exposure`
897 The convolved image.
898 """
899 convolvedExposure = exposure.clone()
900 if psf is not None:
901 convolvedExposure.setPsf(psf)
902 if photoCalib is not None:
903 convolvedExposure.setPhotoCalib(photoCalib)
904 if interpolateBadMaskPlanes and self.config.badMaskPlanes is not None:
905 nInterp = _interpolateImage(convolvedExposure.maskedImage,
906 self.config.badMaskPlanes)
907 self.metadata["nInterpolated"] = nInterp
909 # Snapshot the footprints of mask planes that must not be dilated by
910 # the convolution.
911 preservePlanes = [mp for mp in self.config.preserveMaskPlanes
912 if mp in convolvedExposure.mask.getMaskPlaneDict()]
913 maskResetDict = {
914 mp: (convolvedExposure.mask.array
915 & convolvedExposure.mask.getPlaneBitMask(mp)) > 0
916 for mp in preservePlanes
917 }
919 convolvedImage = lsst.afw.image.MaskedImageF(convolvedExposure.getBBox())
920 lsst.afw.math.convolve(convolvedImage, convolvedExposure.maskedImage, kernel, convolutionControl)
921 convolvedExposure.setMaskedImage(convolvedImage)
923 # Undo the convolution's dilation of the preserved planes: clear the
924 # dilated bits, then restore each mask plane for the pixels that were
925 # previously set.
926 self._clearMask(convolvedExposure.mask, clearMaskPlanes=preservePlanes)
927 for maskPlane, maskSetPixels in maskResetDict.items():
928 bit = convolvedExposure.mask.getPlaneBitMask(maskPlane)
929 convolvedExposure.mask.array[maskSetPixels] |= bit
931 if bbox is None:
932 return convolvedExposure
933 else:
934 return convolvedExposure[bbox]
936 def _sourceSelector(self, template, science, sources, fallback=False):
937 """Select sources from a catalog that meet the selection criteria.
938 The selection criteria include any configured parameters of the
939 `sourceSelector` subtask, as well as checking the science and template
940 mask planes.
942 Parameters
943 ----------
944 template : `lsst.afw.image.ExposureF`
945 Template exposure, warped to match the science exposure.
946 science : `lsst.afw.image.ExposureF`
947 Science exposure to subtract from the template.
948 sources : `lsst.afw.table.SourceCatalog`
949 Input source catalog to select sources from.
950 fallback : `bool`, optional
951 Switch indicating the source selector is being called after
952 running the fallback source detection subtask, which does not run a
953 full set of measurement plugins and can't use the same settings for
954 the source selector.
956 Returns
957 -------
958 kernelSources : `lsst.afw.table.SourceCatalog`
959 The input source catalog, with flagged and low signal-to-noise
960 sources removed and footprints added.
962 Raises
963 ------
964 InsufficientKernelSourcesError
965 An AlgorithmError that is raised if there are not enough PSF
966 candidates to construct the PSF matching kernel.
967 """
968 if fallback:
969 selected = self.fallbackSourceSelector.selectSources(sources).selected
970 else:
971 selected = self.sourceSelector.selectSources(sources).selected
972 # It is OK to use just self.matchedPsfSize for the science PSF here,
973 # since we are just using it to calculate the size of the matching
974 # kernel.
975 kSize = self.makeKernel.makeKernelBasisList(self.templatePsfSize, self.matchedPsfSize)[0].getWidth()
976 selectSources = sources[selected].copy(deep=True)
977 # Set the footprints, to be used in `makeKernel` and `checkMask`.
978 kernelSources = setSourceFootprints(selectSources, kernelSize=kSize)
979 bbox = science.getBBox()
980 if self.usePreconvolution:
981 # Exclude a wider buffer around the edge of the image to
982 # account for an extra convolution.
983 bbox.grow(-kSize)
984 if self.config.restrictKernelEdgeSources:
985 bbox.grow(-kSize)
986 # Remove sources that land on masked pixels
987 scienceSelected = checkMask(science.mask[bbox], kernelSources, self.config.excludeMaskPlanes)
988 templateSelected = checkMask(template.mask[bbox], kernelSources, self.config.excludeMaskPlanes)
989 maskSelected = scienceSelected & templateSelected
990 kernelSources = kernelSources[maskSelected].copy(deep=True)
991 # Trim kernelSources if they exceed ``maxKernelSources``.
992 # Keep the highest signal-to-noise sources of those selected.
993 if (len(kernelSources) > self.config.maxKernelSources) & (self.config.maxKernelSources > 0):
994 signalToNoise = kernelSources.getPsfInstFlux()/kernelSources.getPsfInstFluxErr()
995 indices = np.argsort(signalToNoise)
996 indices = indices[-self.config.maxKernelSources:]
997 selected = np.zeros(len(kernelSources), dtype=bool)
998 selected[indices] = True
999 kernelSources = kernelSources[selected].copy(deep=True)
1001 self.log.info("%i/%i=%.1f%% of sources selected for PSF matching from the input catalog",
1002 len(kernelSources), len(sources), 100*len(kernelSources)/len(sources))
1003 if len(kernelSources) < self.config.minKernelSources:
1004 self.log.error("Too few sources to calculate the PSF matching kernel: "
1005 "%i selected but %i needed for the calculation.",
1006 len(kernelSources), self.config.minKernelSources)
1007 if self.config.allowKernelSourceDetection and not fallback: 1007 ↛ 1011line 1007 didn't jump to line 1011 because the condition on line 1007 was never true
1008 # The fallback source detection pipeline calls this method, so
1009 # allowing source detection in that case would create an endless
1010 # loop
1011 kernelSources = self.runKernelSourceDetection(template, science)
1012 else:
1013 raise InsufficientKernelSourcesError(nSources=len(kernelSources),
1014 nRequired=self.config.minKernelSources)
1016 self.metadata["nPsfSources"] = len(kernelSources)
1018 return kernelSources
1020 def _prepareInputs(self, template, science, visitSummary=None):
1021 """Perform preparatory calculations common to all Alard&Lupton Tasks.
1023 Parameters
1024 ----------
1025 template : `lsst.afw.image.ExposureF`
1026 Template exposure, warped to match the science exposure. The
1027 variance plane of the template image is modified in place.
1028 science : `lsst.afw.image.ExposureF`
1029 Science exposure to subtract from the template. The variance plane
1030 of the science image is modified in place.
1031 visitSummary : `lsst.afw.table.ExposureCatalog`, optional
1032 Exposure catalog with external calibrations to be applied. Catalog
1033 uses the detector id for the catalog id, sorted on id for fast
1034 lookup.
1035 """
1036 self._validateExposures(template, science)
1037 if visitSummary is not None: 1037 ↛ 1038line 1037 didn't jump to line 1038 because the condition on line 1037 was never true
1038 self._applyExternalCalibrations(science, visitSummary=visitSummary)
1039 templateCoverageFraction = checkTemplateIsSufficient(
1040 template[science.getBBox()], science, self.log,
1041 requiredTemplateFraction=self.config.requiredTemplateFraction,
1042 exceptionMessage="Not attempting subtraction. To force subtraction,"
1043 " set config requiredTemplateFraction=0"
1044 )
1045 self.metadata["templateCoveragePercent"] = 100*templateCoverageFraction
1047 if self.config.doScaleVariance:
1048 # Scale the variance of the science image before
1049 # convolution, subtraction, or decorrelation so that it has the
1050 # correct ratio. Note that the template variance is scaled
1051 # independently in ``GetTemplateTask``.
1052 sciVarFactor = self.scaleVariance.run(science.maskedImage)
1053 self.log.info("Science variance scaling factor: %.2f", sciVarFactor)
1054 self.metadata["scaleScienceVarianceFactor"] = sciVarFactor
1056 # Erase existing detection mask planes.
1057 # We don't want the detection mask from the science image
1058 self.updateMasks(template, science)
1060 # Calling getPsfFwhm on template.psf fails on some rare occasions when
1061 # the template has no input exposures at the average position of the
1062 # stars. So we try getPsfFwhm first on template, and if that fails we
1063 # evaluate the PSF on a grid specified by fwhmExposure* fields.
1064 # To keep consistent definitions for PSF size on the template and
1065 # science images, we use the same method for both.
1066 # In the try block below, we catch two exceptions:
1067 # 1. InvalidParameterError, in case the point where we are evaluating
1068 # the PSF lands in a gap in the template.
1069 # 2. RangeError, in case the template coverage is so poor that we end
1070 # up near a region with no data.
1071 try:
1072 self.templatePsfSize = getPsfFwhm(template.psf)
1073 self.sciencePsfSize = getPsfFwhm(science.psf)
1074 except lsst.pex.exceptions.Exception:
1075 # Catch a broad range of exceptions, since some are C++ only
1076 # Catching:
1077 # - lsst::geom::SingularTransformException
1078 # - lsst.pex.exceptions.InvalidParameterError
1079 # - lsst.pex.exceptions.RangeError
1080 self.log.info("Unable to evaluate PSF at the average position. "
1081 "Evaluting PSF on a grid of points."
1082 )
1083 self.templatePsfSize = evaluateMeanPsfFwhm(
1084 template,
1085 fwhmExposureBuffer=self.config.makeKernel.fwhmExposureBuffer,
1086 fwhmExposureGrid=self.config.makeKernel.fwhmExposureGrid
1087 )
1088 self.sciencePsfSize = evaluateMeanPsfFwhm(
1089 science,
1090 fwhmExposureBuffer=self.config.makeKernel.fwhmExposureBuffer,
1091 fwhmExposureGrid=self.config.makeKernel.fwhmExposureGrid
1092 )
1093 self.log.info("Science PSF FWHM: %f pixels", self.sciencePsfSize)
1094 self.log.info("Template PSF FWHM: %f pixels", self.templatePsfSize)
1095 self.metadata["sciencePsfSize"] = self.sciencePsfSize
1096 self.metadata["templatePsfSize"] = self.templatePsfSize
1098 # Calculate estimated image depths, i.e., limiting magnitudes
1099 maglim_science = self._calculateMagLim(science, fallbackPsfSize=self.sciencePsfSize)
1100 if np.isnan(maglim_science): 1100 ↛ 1101line 1100 didn't jump to line 1101 because the condition on line 1100 was never true
1101 self.log.warning("Limiting magnitude of the science image is NaN!")
1102 fluxlim_science = (maglim_science*u.ABmag).to_value(u.nJy)
1103 maglim_template = self._calculateMagLim(template, fallbackPsfSize=self.templatePsfSize)
1104 if np.isnan(maglim_template): 1104 ↛ 1105line 1104 didn't jump to line 1105 because the condition on line 1104 was never true
1105 self.log.info("Cannot evaluate template limiting mag; adopting science limiting mag for diffim")
1106 maglim_diffim = maglim_science
1107 else:
1108 fluxlim_template = (maglim_template*u.ABmag).to_value(u.nJy)
1109 maglim_diffim = (np.sqrt(fluxlim_science**2 + fluxlim_template**2)*u.nJy).to(u.ABmag).value
1110 self.metadata["scienceLimitingMagnitude"] = maglim_science
1111 self.metadata["templateLimitingMagnitude"] = maglim_template
1112 self.metadata["diffimLimitingMagnitude"] = maglim_diffim
1114 def updateMasks(self, template, science):
1115 """Update the science and template mask planes before differencing.
1117 Parameters
1118 ----------
1119 template : `lsst.afw.image.Exposure`
1120 Template exposure, warped to match the science exposure.
1121 The template mask planes will be erased, except for a few specified
1122 in the task config.
1123 science : `lsst.afw.image.Exposure`
1124 Science exposure to subtract from the template.
1125 The DETECTED and DETECTED_NEGATIVE mask planes of the science image
1126 will be erased.
1127 """
1128 self._clearMask(science.mask, clearMaskPlanes=["DETECTED", "DETECTED_NEGATIVE"])
1130 # We will clear ALL template mask planes, except for those specified
1131 # via the `preserveTemplateMask` config. Mask planes specified via
1132 # the `renameTemplateMask` config will be copied to new planes with
1133 # "_TEMPLATE" appended to their names, and the original mask plane will
1134 # be cleared.
1135 clearMaskPlanes = [mp for mp in template.mask.getMaskPlaneDict().keys()
1136 if mp not in self.config.preserveTemplateMask]
1137 renameMaskPlanes = [mp for mp in self.config.renameTemplateMask
1138 if mp in template.mask.getMaskPlaneDict().keys()]
1140 # propagate the mask plane related to Fake source injection
1141 # NOTE: the fake source injection sets FAKE plane, but it should be INJECTED
1142 # NOTE: This can be removed in DM-40796
1143 if "FAKE" in science.mask.getMaskPlaneDict().keys():
1144 self.log.info("Adding injected mask plane to science image")
1145 self._renameMaskPlanes(science.mask, "FAKE", "INJECTED")
1146 if "FAKE" in template.mask.getMaskPlaneDict().keys():
1147 self.log.info("Adding injected mask plane to template image")
1148 self._renameMaskPlanes(template.mask, "FAKE", "INJECTED_TEMPLATE")
1149 if "INJECTED" in renameMaskPlanes: 1149 ↛ 1151line 1149 didn't jump to line 1151 because the condition on line 1149 was always true
1150 renameMaskPlanes.remove("INJECTED")
1151 if "INJECTED_TEMPLATE" in clearMaskPlanes: 1151 ↛ 1154line 1151 didn't jump to line 1154 because the condition on line 1151 was always true
1152 clearMaskPlanes.remove("INJECTED_TEMPLATE")
1154 for maskPlane in renameMaskPlanes:
1155 self._renameMaskPlanes(template.mask, maskPlane, maskPlane + "_TEMPLATE")
1156 self._clearMask(template.mask, clearMaskPlanes=clearMaskPlanes)
1158 @staticmethod
1159 def _renameMaskPlanes(mask, maskPlane, newMaskPlane):
1160 """Rename a mask plane by adding the new name and copying the data.
1162 Parameters
1163 ----------
1164 mask : `lsst.afw.image.Mask`
1165 The mask image to update in place.
1166 maskPlane : `str`
1167 The name of the existing mask plane to copy.
1168 newMaskPlane : `str`
1169 The new name of the mask plane that will be added.
1170 If the mask plane already exists, it will be updated in place.
1171 """
1172 mask.addMaskPlane(newMaskPlane)
1173 originBitMask = mask.getPlaneBitMask(maskPlane)
1174 destinationBitMask = mask.getPlaneBitMask(newMaskPlane)
1175 mask.array |= ((mask.array & originBitMask) > 0)*destinationBitMask
1177 def _clearMask(self, mask, clearMaskPlanes=None):
1178 """Clear the mask plane of an exposure.
1180 Parameters
1181 ----------
1182 mask : `lsst.afw.image.Mask`
1183 The mask plane to erase, which will be modified in place.
1184 clearMaskPlanes : `list` of `str`, optional
1185 Erase the specified mask planes.
1186 If not supplied, the entire mask will be erased.
1187 """
1188 if clearMaskPlanes is None: 1188 ↛ 1189line 1188 didn't jump to line 1189 because the condition on line 1188 was never true
1189 clearMaskPlanes = list(mask.getMaskPlaneDict().keys())
1191 bitMaskToClear = mask.getPlaneBitMask(clearMaskPlanes)
1192 mask &= ~bitMaskToClear
1195class AlardLuptonPreconvolveSubtractConnections(SubtractInputConnections,
1196 SubtractScoreOutputConnections):
1197 pass
1200class AlardLuptonPreconvolveSubtractConfig(AlardLuptonSubtractBaseConfig, lsst.pipe.base.PipelineTaskConfig,
1201 pipelineConnections=AlardLuptonPreconvolveSubtractConnections):
1202 pass
1205class AlardLuptonPreconvolveSubtractTask(AlardLuptonSubtractTask):
1206 """Subtract a template from a science image, convolving the science image
1207 before computing the kernel, and also convolving the template before
1208 subtraction.
1209 """
1210 ConfigClass = AlardLuptonPreconvolveSubtractConfig
1211 _DefaultName = "alardLuptonPreconvolveSubtract"
1212 usePreconvolution = True
1214 def run(self, template, science, sources, visitSummary=None):
1215 """Preconvolve the science image with its own PSF,
1216 convolve the template image with a PSF-matching kernel and subtract
1217 from the preconvolved science image.
1219 Parameters
1220 ----------
1221 template : `lsst.afw.image.ExposureF`
1222 The template image, which has previously been warped to the science
1223 image. The template bbox will be padded by a few pixels compared to
1224 the science bbox.
1225 science : `lsst.afw.image.ExposureF`
1226 The science exposure.
1227 sources : `lsst.afw.table.SourceCatalog`
1228 Identified sources on the science exposure. This catalog is used to
1229 select sources in order to perform the AL PSF matching on stamp
1230 images around them.
1231 visitSummary : `lsst.afw.table.ExposureCatalog`, optional
1232 Exposure catalog with complete external calibrations. Catalog uses
1233 the detector id for the catalog id, sorted on id for fast lookup.
1235 Returns
1236 -------
1237 results : `lsst.pipe.base.Struct`
1238 ``scoreExposure`` : `lsst.afw.image.ExposureF`
1239 Result of subtracting the convolved template and science
1240 images. Attached PSF is that of the original science image.
1241 ``matchedTemplate`` : `lsst.afw.image.ExposureF`
1242 Warped and PSF-matched template exposure. Attached PSF is that
1243 of the original science image.
1244 ``matchedScience`` : `lsst.afw.image.ExposureF`
1245 The science exposure after convolving with its own PSF.
1246 Attached PSF is that of the original science image.
1247 ``backgroundModel`` : `lsst.afw.math.Function2D`
1248 Background model that was fit while solving for the
1249 PSF-matching kernel
1250 ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
1251 Final kernel used to PSF-match the template to the science
1252 image.
1253 """
1254 self._prepareInputs(template, science, visitSummary=visitSummary)
1256 convolutionKernel = self._makePreconvolutionKernel(science.psf)
1257 matchedScience = self._convolveExposure(science, convolutionKernel, self.convolutionControl,
1258 interpolateBadMaskPlanes=True)
1259 self.metadata["convolvedExposure"] = "Preconvolution"
1261 self.matchedPsfSize = self.sciencePsfSize*np.sqrt(2)
1262 self.log.info("Preconvolved science PSF FWHM: %f pixels", self.matchedPsfSize)
1263 self.metadata["preconvolvedSciencePsfSize"] = self.matchedPsfSize
1264 try:
1265 kernelSources = self._sourceSelector(template, matchedScience, sources)
1266 subtractResults = self.runPreconvolve(template, science, matchedScience,
1267 kernelSources, convolutionKernel)
1269 except (RuntimeError, lsst.pex.exceptions.Exception) as e:
1270 self.log.warning("Failed to match template. Checking coverage")
1271 # Raise NoWorkFound if template fraction is insufficient
1272 checkTemplateIsSufficient(template[science.getBBox()], science, self.log,
1273 self.config.minTemplateFractionForExpectedSuccess,
1274 exceptionMessage="Template coverage lower than expected to succeed."
1275 f" Failure is tolerable: {e}")
1276 # checkTemplateIsSufficient did not raise NoWorkFound, so raise original exception
1277 raise e
1279 return subtractResults
1281 @staticmethod
1282 def _flagScoreEdge(mask, innerBBox):
1283 """Set the EDGE mask bit on pixels outside a known-valid region.
1285 Parameters
1286 ----------
1287 mask : `~lsst.afw.image.Mask`
1288 Exposure mask that will be modified in place. Must have
1289 an ``EDGE`` mask plane.
1290 innerBBox : `~lsst.geom.Box2I`
1291 The valid inner region. Pixels
1292 outside this bbox will have their ``EDGE`` bit set.
1293 """
1294 bbox = mask.getBBox()
1295 edgeBit = mask.getPlaneBitMask("EDGE")
1296 dx0 = innerBBox.getMinX() - bbox.getMinX()
1297 dx1 = bbox.getMaxX() - innerBBox.getMaxX()
1298 dy0 = innerBBox.getMinY() - bbox.getMinY()
1299 dy1 = bbox.getMaxY() - innerBBox.getMaxY()
1300 if dy0 > 0: 1300 ↛ 1302line 1300 didn't jump to line 1302 because the condition on line 1300 was always true
1301 mask.array[:dy0, :] |= edgeBit
1302 if dy1 > 0: 1302 ↛ 1304line 1302 didn't jump to line 1304 because the condition on line 1302 was always true
1303 mask.array[-dy1:, :] |= edgeBit
1304 if dx0 > 0: 1304 ↛ 1306line 1304 didn't jump to line 1306 because the condition on line 1304 was always true
1305 mask.array[:, :dx0] |= edgeBit
1306 if dx1 > 0: 1306 ↛ exitline 1306 didn't return from function '_flagScoreEdge' because the condition on line 1306 was always true
1307 mask.array[:, -dx1:] |= edgeBit
1309 @staticmethod
1310 def _makePreconvolutionKernel(psf):
1311 """Build a normalized, reflected matched-filter kernel from a PSF.
1313 Convolving an image with this kernel is equivalent to correlating
1314 the image with the PSF, so peaks in the output align with the PSF's
1315 centroid — even for asymmetric PSFs. The kernel is evaluated at the
1316 PSF's average position and returned as a constant
1317 `~lsst.afw.math.Kernel`.
1319 Parameters
1320 ----------
1321 psf : `~lsst.afw.detection.Psf`
1322 The PSF to derive the preconvolution kernel from.
1324 Returns
1325 -------
1326 kernel : `~lsst.afw.math.Kernel`
1327 The PSF reflected about both axes, normalized to sum to one.
1329 Raises
1330 ------
1331 ValueError
1332 Raised if the PSF kernel has an even size along either axis.
1333 It's not possible to center an even-sized kernel.
1334 """
1335 avgPos = psf.getAveragePosition()
1336 localKernel = psf.getLocalKernel(avgPos)
1337 dims = localKernel.getDimensions()
1338 if dims.x % 2 == 0 or dims.y % 2 == 0: 1338 ↛ 1339line 1338 didn't jump to line 1339 because the condition on line 1338 was never true
1339 raise ValueError(
1340 f"Preconvolution requires an odd-sized PSF kernel, got {dims.x}x{dims.y}. "
1341 )
1342 kimg = lsst.afw.image.ImageD(dims)
1343 localKernel.computeImage(kimg, doNormalize=True) # normalize to unit sum
1344 # Reflect about the kernel center. PSF kernels are odd-sized,
1345 # so ``[::-1, ::-1]`` places the peak at the same pixel.
1346 kimg.array[...] = kimg.array[::-1, ::-1]
1347 return lsst.afw.math.FixedKernel(kimg)
1349 def runPreconvolve(self, template, science, matchedScience, kernelSources, preConvKernel):
1350 """Convolve the science image with its own PSF, then convolve the
1351 template with a matching kernel and subtract to form the Score
1352 exposure.
1354 Parameters
1355 ----------
1356 template : `lsst.afw.image.ExposureF`
1357 Template exposure, warped to match the science exposure.
1358 science : `lsst.afw.image.ExposureF`
1359 Science exposure to subtract from the template.
1360 matchedScience : `lsst.afw.image.ExposureF`
1361 The science exposure, convolved with the reflection of its own PSF.
1362 kernelSources : `lsst.afw.table.SourceCatalog`
1363 Identified sources on the science exposure. This catalog is used to
1364 select sources in order to perform the AL PSF matching on stamp
1365 images around them.
1366 preConvKernel : `lsst.afw.math.Kernel`
1367 The kernel that was used to preconvolve the ``science``
1368 exposure. Must be normalized to sum to 1.
1370 Returns
1371 -------
1372 results : `lsst.pipe.base.Struct`
1374 ``scoreExposure`` : `lsst.afw.image.ExposureF`
1375 Result of subtracting the convolved template and science
1376 images. Attached PSF is that of the original science image.
1377 ``matchedTemplate`` : `lsst.afw.image.ExposureF`
1378 Warped and PSF-matched template exposure. Attached PSF is that
1379 of the original science image.
1380 ``matchedScience`` : `lsst.afw.image.ExposureF`
1381 The science exposure after convolving with its own PSF.
1382 Attached PSF is that of the original science image.
1383 ``backgroundModel`` : `lsst.afw.math.Function2D`
1384 Background model that was fit while solving for the
1385 PSF-matching kernel
1386 ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
1387 Final kernel used to PSF-match the template to the science
1388 image.
1389 """
1390 bbox = science.getBBox()
1391 innerBBox = preConvKernel.shrinkBBox(bbox)
1393 kernelResult = self.makeKernel.run(template[innerBBox], matchedScience[innerBBox], kernelSources,
1394 preconvolved=True,
1395 templateFwhmPix=self.templatePsfSize,
1396 scienceFwhmPix=self.matchedPsfSize)
1398 matchedTemplate = self._convolveExposure(template, kernelResult.psfMatchingKernel,
1399 self.convolutionControl,
1400 bbox=bbox,
1401 psf=science.psf,
1402 interpolateBadMaskPlanes=True,
1403 photoCalib=science.photoCalib)
1404 score = _subtractImages(matchedScience, matchedTemplate,
1405 backgroundModel=(kernelResult.backgroundModel
1406 if self.config.doSubtractBackground else None))
1407 correctedScore = self.finalize(template[bbox], science, score,
1408 kernelResult.psfMatchingKernel,
1409 templateMatched=True, preConvMode=True,
1410 preConvKernel=preConvKernel)
1412 # Flag the outer ``preConvKernel/2``-wide border as EDGE.
1413 self._flagScoreEdge(correctedScore.mask, innerBBox)
1415 return lsst.pipe.base.Struct(scoreExposure=correctedScore,
1416 matchedTemplate=matchedTemplate,
1417 matchedScience=matchedScience,
1418 backgroundModel=kernelResult.backgroundModel,
1419 psfMatchingKernel=kernelResult.psfMatchingKernel,
1420 kernelSources=kernelSources)
1423def checkTemplateIsSufficient(templateExposure, scienceExposure, logger, requiredTemplateFraction=0.,
1424 exceptionMessage=""):
1425 """Raise NoWorkFound if template coverage < requiredTemplateFraction
1427 Parameters
1428 ----------
1429 templateExposure : `lsst.afw.image.ExposureF`
1430 The template exposure to check
1431 logger : `logging.Logger`
1432 Logger for printing output.
1433 requiredTemplateFraction : `float`, optional
1434 Fraction of pixels of the science image required to have coverage
1435 in the template.
1436 exceptionMessage : `str`, optional
1437 Message to include in the exception raised if the template coverage
1438 is insufficient.
1440 Returns
1441 -------
1442 templateCoverageFraction: `float`
1443 Fraction of pixels in the template with data.
1445 Raises
1446 ------
1447 lsst.pipe.base.NoWorkFound
1448 Raised if fraction of good pixels, defined as not having NO_DATA
1449 set, is less than the requiredTemplateFraction
1450 """
1451 # Count the number of pixels with the NO_DATA mask bit set
1452 # counting NaN pixels is insufficient because pixels without data are often intepolated over)
1453 noTemplate = templateExposure.mask.array & templateExposure.mask.getPlaneBitMask('NO_DATA')
1454 # Also need to account for missing data in the science image,
1455 # because template coverage there doesn't help
1456 noScience = scienceExposure.mask.array & scienceExposure.mask.getPlaneBitMask('NO_DATA')
1457 pixNoData = np.count_nonzero(noTemplate | noScience)
1458 pixGood = templateExposure.getBBox().getArea() - pixNoData
1459 templateCoverageFraction = pixGood/templateExposure.getBBox().getArea()
1460 logger.info("template has %d good pixels (%.1f%%)", pixGood, 100*templateCoverageFraction)
1462 if templateCoverageFraction < requiredTemplateFraction:
1463 message = ("Insufficient Template Coverage. (%.1f%% < %.1f%%)" % (
1464 100*templateCoverageFraction,
1465 100*requiredTemplateFraction))
1466 raise lsst.pipe.base.NoWorkFound(message + " " + exceptionMessage)
1467 return templateCoverageFraction
1470def _subtractImages(science, template, backgroundModel=None):
1471 """Subtract template from science, propagating relevant metadata.
1473 Parameters
1474 ----------
1475 science : `lsst.afw.Exposure`
1476 The input science image.
1477 template : `lsst.afw.Exposure`
1478 The template to subtract from the science image.
1479 backgroundModel : `lsst.afw.MaskedImage`, optional
1480 Differential background model
1482 Returns
1483 -------
1484 difference : `lsst.afw.Exposure`
1485 The subtracted image.
1486 """
1487 difference = science.clone()
1488 if backgroundModel is not None:
1489 difference.maskedImage -= backgroundModel
1490 difference.maskedImage -= template.maskedImage
1491 return difference
1494def _shapeTest(exp1, exp2, fwhmExposureBuffer, fwhmExposureGrid):
1495 """Determine that the PSF of ``exp1`` is not wider than that of ``exp2``.
1497 Parameters
1498 ----------
1499 exp1 : `~lsst.afw.image.Exposure`
1500 Exposure with the reference point spread function (PSF) to evaluate.
1501 exp2 : `~lsst.afw.image.Exposure`
1502 Exposure with a candidate point spread function (PSF) to evaluate.
1503 fwhmExposureBuffer : `float`
1504 Fractional buffer margin to be left out of all sides of the image
1505 during the construction of the grid to compute mean PSF FWHM in an
1506 exposure, if the PSF is not available at its average position.
1507 fwhmExposureGrid : `int`
1508 Grid size to compute the mean FWHM in an exposure, if the PSF is not
1509 available at its average position.
1510 Returns
1511 -------
1512 result : `bool`
1513 True if ``exp1`` has a PSF that is not wider than that of ``exp2`` in
1514 either dimension.
1515 """
1516 try:
1517 shape1 = getPsfFwhm(exp1.psf, average=False)
1518 shape2 = getPsfFwhm(exp2.psf, average=False)
1519 except (lsst.pex.exceptions.InvalidParameterError, lsst.pex.exceptions.RangeError):
1520 shape1 = evaluateMeanPsfFwhm(exp1,
1521 fwhmExposureBuffer=fwhmExposureBuffer,
1522 fwhmExposureGrid=fwhmExposureGrid
1523 )
1524 shape2 = evaluateMeanPsfFwhm(exp2,
1525 fwhmExposureBuffer=fwhmExposureBuffer,
1526 fwhmExposureGrid=fwhmExposureGrid
1527 )
1528 return shape1 <= shape2
1530 # Results from getPsfFwhm is a tuple of two values, one for each dimension.
1531 xTest = shape1[0] <= shape2[0]
1532 yTest = shape1[1] <= shape2[1]
1533 return xTest | yTest
1536class SimplifiedSubtractConfig(AlardLuptonSubtractBaseConfig, lsst.pipe.base.PipelineTaskConfig,
1537 pipelineConnections=SimplifiedSubtractConnections):
1538 mode = lsst.pex.config.ChoiceField(
1539 dtype=str,
1540 default="convolveTemplate",
1541 allowed={"auto": "Choose which image to convolve at runtime.",
1542 "convolveScience": "Only convolve the science image.",
1543 "convolveTemplate": "Only convolve the template image."},
1544 doc="Choose which image to convolve at runtime, or require that a specific image is convolved."
1545 )
1546 useExistingKernel = lsst.pex.config.Field(
1547 dtype=bool,
1548 default=True,
1549 doc="Use a pre-existing PSF matching kernel?"
1550 "If False, source detection and measurement will be run."
1551 )
1554class SimplifiedSubtractTask(AlardLuptonSubtractTask):
1555 """Compute the image difference of a science and template image using
1556 the Alard & Lupton (1998) algorithm.
1557 """
1558 ConfigClass = SimplifiedSubtractConfig
1559 _DefaultName = "simplifiedSubtract"
1561 @timeMethod
1562 def run(self, template, science, visitSummary=None, inputPsfMatchingKernel=None):
1563 """PSF match, subtract, and decorrelate two images.
1565 Parameters
1566 ----------
1567 template : `lsst.afw.image.ExposureF`
1568 Template exposure, warped to match the science exposure.
1569 science : `lsst.afw.image.ExposureF`
1570 Science exposure to subtract from the template.
1571 visitSummary : `lsst.afw.table.ExposureCatalog`, optional
1572 Exposure catalog with external calibrations to be applied. Catalog
1573 uses the detector id for the catalog id, sorted on id for fast
1574 lookup.
1575 inputPsfMatchingKernel : `lsst.afw.math.Kernel`, optional
1576 Pre-existing PSF matching kernel to use for convolution.
1577 Required, and only used, if ``config.useExistingKernel`` is set.
1579 Returns
1580 -------
1581 results : `lsst.pipe.base.Struct`
1582 ``difference`` : `lsst.afw.image.ExposureF`
1583 Result of subtracting template and science.
1584 ``matchedTemplate`` : `lsst.afw.image.ExposureF`
1585 Warped and PSF-matched template exposure.
1586 ``backgroundModel`` : `lsst.afw.math.Function2D`
1587 Background model that was fit while solving for the
1588 PSF-matching kernel
1589 ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
1590 Kernel used to PSF-match the convolved image.
1591 ``kernelSources` : `lsst.afw.table.SourceCatalog`
1592 Sources detected on the science image that were used to
1593 construct the PSF-matching kernel.
1595 Raises
1596 ------
1597 lsst.pipe.base.NoWorkFound
1598 Raised if fraction of good pixels, defined as not having NO_DATA
1599 set, is less then the configured requiredTemplateFraction
1600 """
1601 self._prepareInputs(template, science, visitSummary=visitSummary)
1603 convolveTemplate = self.chooseConvolutionMethod(template, science)
1604 self.matchedPsfSize = self.sciencePsfSize if convolveTemplate else self.templatePsfSize
1606 if self.config.useExistingKernel:
1607 psfMatchingKernel = inputPsfMatchingKernel
1608 backgroundModel = None
1609 kernelSources = None
1610 else:
1611 kernelResult = self.runMakeKernel(template, science, convolveTemplate=convolveTemplate,
1612 runSourceDetection=True)
1613 psfMatchingKernel = kernelResult.psfMatchingKernel
1614 kernelSources = kernelResult.kernelSources
1615 if self.config.doSubtractBackground: 1615 ↛ 1616line 1615 didn't jump to line 1616 because the condition on line 1615 was never true
1616 backgroundModel = kernelResult.backgroundModel
1617 else:
1618 backgroundModel = None
1619 if convolveTemplate: 1619 ↛ 1623line 1619 didn't jump to line 1623 because the condition on line 1619 was always true
1620 subtractResults = self.runConvolveTemplate(template, science, psfMatchingKernel,
1621 backgroundModel=backgroundModel)
1622 else:
1623 subtractResults = self.runConvolveScience(template, science, psfMatchingKernel,
1624 backgroundModel=backgroundModel)
1625 if kernelSources is not None:
1626 subtractResults.kernelSources = kernelSources
1627 return subtractResults
1630def _interpolateImage(maskedImage, badMaskPlanes, fallbackValue=None):
1631 """Replace masked image pixels with interpolated values.
1633 Parameters
1634 ----------
1635 maskedImage : `lsst.afw.image.MaskedImage`
1636 Image on which to perform interpolation.
1637 badMaskPlanes : `list` of `str`
1638 List of mask planes to interpolate over.
1639 fallbackValue : `float`, optional
1640 Value to set when interpolation fails.
1642 Returns
1643 -------
1644 result: `float`
1645 The number of masked pixels that were replaced.
1646 """
1647 imgBadMaskPlanes = [
1648 maskPlane for maskPlane in badMaskPlanes if maskPlane in maskedImage.mask.getMaskPlaneDict()
1649 ]
1651 image = maskedImage.image.array
1652 badPixels = (maskedImage.mask.array & maskedImage.mask.getPlaneBitMask(imgBadMaskPlanes)) > 0
1653 image[badPixels] = np.nan
1654 if fallbackValue is None: 1654 ↛ 1658line 1654 didn't jump to line 1658 because the condition on line 1654 was always true
1655 fallbackValue = np.nanmedian(image)
1656 # For this initial implementation, skip the interpolation and just fill with
1657 # the median value.
1658 image[badPixels] = fallbackValue
1659 return np.sum(badPixels)