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