Coverage for python/lsst/ip/diffim/subtractImages.py: 90%
566 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-14 19:10 +0000
1# This file is part of ip_diffim.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
22import dataclasses
24from astropy import units as u
25from astropy.stats import gaussian_fwhm_to_sigma
26import numpy as np
28import lsst.afw.detection as afwDetection
29import lsst.afw.image
30import lsst.afw.math
31import lsst.geom
32from lsst.images import DifferenceImage
33from lsst.images.convolution_kernels import ImageBasisConvolutionKernel
34from lsst.ip.diffim.utils import (evaluateMeanPsfFwhm, getPsfFwhm,
35 computeDifferenceImageMetrics, computePSFNoiseEquivalentArea,
36 checkMask, record_from_obs_info, setSourceFootprints,
37 get_difference_image_provenance)
38from lsst.meas.algorithms import ScaleVarianceTask, ScienceSourceSelectorTask
39import lsst.pex.config
40import lsst.pipe.base
41import lsst.pex.exceptions
42from lsst.pipe.base import connectionTypes
43from . import MakeKernelTask, DecorrelateALKernelTask
44from lsst.utils.timer import timeMethod
46__all__ = ["AlardLuptonSubtractConfig", "AlardLuptonSubtractTask",
47 "AlardLuptonPreconvolveSubtractConfig", "AlardLuptonPreconvolveSubtractTask",
48 "SimplifiedSubtractConfig", "SimplifiedSubtractTask",
49 "InsufficientKernelSourcesError"]
51_dimensions = ("instrument", "visit", "detector")
52_defaultTemplates = {"coaddName": "deep", "fakesType": ""}
55class InsufficientKernelSourcesError(lsst.pipe.base.AlgorithmError):
56 """Raised when there are too few sources to calculate the PSF matching
57 kernel.
58 """
59 def __init__(self, *, nSources, nRequired):
60 msg = (f"Only {nSources} sources were selected for PSF matching,"
61 f" but {nRequired} are required.")
62 super().__init__(msg)
63 self.nSources = nSources
64 self.nRequired = nRequired
66 @property
67 def metadata(self):
68 return {"nSources": self.nSources,
69 "nRequired": self.nRequired
70 }
73class SubtractInputConnections(lsst.pipe.base.PipelineTaskConnections,
74 dimensions=_dimensions,
75 defaultTemplates=_defaultTemplates):
76 template = connectionTypes.Input(
77 doc="Input warped template to subtract."
78 " With image_type='future' this is read as a lsst.images.DifferenceImage, so that the"
79 " record of the coadds that went into it can be passed on to the outputs.",
80 dimensions=("instrument", "visit", "detector"),
81 storageClass="ExposureF",
82 name="{fakesType}{coaddName}Diff_templateExp"
83 )
84 science = connectionTypes.Input(
85 doc="Input science exposure to subtract from."
86 " With image_type='future' this is read as a lsst.images.VisitImage, so that its"
87 " photometric scaling can be passed on to the outputs.",
88 dimensions=("instrument", "visit", "detector"),
89 storageClass="ExposureF",
90 name="{fakesType}calexp"
91 )
92 sources = connectionTypes.Input(
93 doc="Sources measured on the science exposure; "
94 "used to select sources for making the matching kernel.",
95 dimensions=("instrument", "visit", "detector"),
96 storageClass="SourceCatalog",
97 name="{fakesType}src"
98 )
99 visitSummary = connectionTypes.Input(
100 doc=("Per-visit catalog with final calibration objects. "
101 "These catalogs use the detector id for the catalog id, "
102 "sorted on id for fast lookup."),
103 dimensions=("instrument", "visit"),
104 storageClass="ExposureCatalog",
105 name="finalVisitSummary",
106 )
108 def __init__(self, *, config=None):
109 super().__init__(config=config)
110 if not config.doApplyExternalCalibrations: 110 ↛ 112line 110 didn't jump to line 112 because the condition on line 110 was always true
111 del self.visitSummary
112 if config.image_type == "future":
113 self.template = dataclasses.replace(self.template, storageClass="DifferenceImage")
114 self.science = dataclasses.replace(self.science, storageClass="VisitImage")
117class SubtractImageOutputConnections(lsst.pipe.base.PipelineTaskConnections,
118 dimensions=_dimensions,
119 defaultTemplates=_defaultTemplates):
120 difference = connectionTypes.Output(
121 doc="Result of subtracting convolved template from science image.",
122 dimensions=("instrument", "visit", "detector"),
123 storageClass="ExposureF",
124 name="{fakesType}{coaddName}Diff_differenceTempExp",
125 )
126 matchedTemplate = connectionTypes.Output(
127 doc="Warped and PSF-matched template used to create `subtractedExposure`.",
128 dimensions=("instrument", "visit", "detector"),
129 storageClass="ExposureF",
130 name="{fakesType}{coaddName}Diff_matchedExp",
131 )
132 psfMatchingKernel = connectionTypes.Output(
133 doc="Kernel used to PSF match the science and template images.",
134 dimensions=("instrument", "visit", "detector"),
135 storageClass="MatchingKernel",
136 name="{fakesType}{coaddName}Diff_psfMatchKernel",
137 )
138 kernelSources = connectionTypes.Output(
139 doc="Final selection of sources used for psf matching.",
140 dimensions=("instrument", "visit", "detector"),
141 storageClass="SourceCatalog",
142 name="{fakesType}{coaddName}Diff_psfMatchSources"
143 )
145 def __init__(self, *, config=None):
146 super().__init__(config=config)
147 if config.image_type == "future":
148 self.difference = dataclasses.replace(self.difference, storageClass="DifferenceImage")
149 self.matchedTemplate = dataclasses.replace(self.matchedTemplate,
150 storageClass="DifferenceImage")
151 self.psfMatchingKernel = dataclasses.replace(self.psfMatchingKernel,
152 storageClass="ConvolutionKernel")
155class SubtractScoreOutputConnections(lsst.pipe.base.PipelineTaskConnections,
156 dimensions=_dimensions,
157 defaultTemplates=_defaultTemplates):
158 scoreExposure = connectionTypes.Output(
159 doc="The maximum likelihood image, used for the detection of diaSources.",
160 dimensions=("instrument", "visit", "detector"),
161 storageClass="ExposureF",
162 name="{fakesType}{coaddName}Diff_scoreTempExp",
163 )
164 psfMatchingKernel = connectionTypes.Output(
165 doc="Kernel used to PSF match the science and template images.",
166 dimensions=("instrument", "visit", "detector"),
167 storageClass="MatchingKernel",
168 name="{fakesType}{coaddName}Diff_psfScoreMatchKernel",
169 )
170 kernelSources = connectionTypes.Output(
171 doc="Final selection of sources used for psf matching.",
172 dimensions=("instrument", "visit", "detector"),
173 storageClass="SourceCatalog",
174 name="{fakesType}{coaddName}Diff_psfScoreMatchSources"
175 )
177 def __init__(self, *, config=None):
178 super().__init__(config=config)
179 if config.image_type == "future":
180 self.scoreExposure = dataclasses.replace(self.scoreExposure, storageClass="DifferenceImage")
181 self.psfMatchingKernel = dataclasses.replace(self.psfMatchingKernel,
182 storageClass="ConvolutionKernel")
185class AlardLuptonSubtractConnections(SubtractInputConnections, SubtractImageOutputConnections):
186 pass
189class SimplifiedSubtractConnections(SubtractInputConnections, SubtractImageOutputConnections):
190 inputPsfMatchingKernel = connectionTypes.Input(
191 doc="Kernel used to PSF match the science and template images.",
192 dimensions=("instrument", "visit", "detector"),
193 storageClass="MatchingKernel",
194 name="{fakesType}{coaddName}Diff_psfMatchKernel",
195 )
197 def __init__(self, *, config=None):
198 super().__init__(config=config)
199 del self.sources
200 if config.useExistingKernel: 200 ↛ 207line 200 didn't jump to line 207 because the condition on line 200 was always true
201 del self.psfMatchingKernel
202 del self.kernelSources
203 if config.image_type == "future":
204 self.inputPsfMatchingKernel = dataclasses.replace(self.inputPsfMatchingKernel,
205 storageClass="ConvolutionKernel")
206 else:
207 del self.inputPsfMatchingKernel
210class AlardLuptonSubtractBaseConfig(lsst.pex.config.Config):
211 makeKernel = lsst.pex.config.ConfigurableField(
212 target=MakeKernelTask,
213 doc="Task to construct a matching kernel for convolution.",
214 )
215 doDecorrelation = lsst.pex.config.Field(
216 dtype=bool,
217 default=True,
218 doc="Perform diffim decorrelation to undo pixel correlation due to A&L "
219 "kernel convolution? If True, also update the diffim PSF."
220 )
221 decorrelate = lsst.pex.config.ConfigurableField(
222 target=DecorrelateALKernelTask,
223 doc="Task to decorrelate the image difference.",
224 )
225 requiredTemplateFraction = lsst.pex.config.Field(
226 dtype=float,
227 default=0.1,
228 doc="Raise NoWorkFound and do not attempt image subtraction if template covers less than this "
229 " fraction of pixels. Setting to 0 will always attempt image subtraction."
230 )
231 minTemplateFractionForExpectedSuccess = lsst.pex.config.Field(
232 dtype=float,
233 default=0.2,
234 doc="Raise NoWorkFound if PSF-matching fails and template covers less than this fraction of pixels."
235 " If the fraction of pixels covered by the template is less than this value (and greater than"
236 " requiredTemplateFraction) this task is attempted but failure is anticipated and tolerated."
237 )
238 doScaleVariance = lsst.pex.config.Field(
239 dtype=bool,
240 default=True,
241 doc="Scale variance of the science image? Note that the template variance is NOT scaled"
242 " here. The template variance may be scaled independently in ``GetTemplateTask``."
243 )
244 scaleVariance = lsst.pex.config.ConfigurableField(
245 target=ScaleVarianceTask,
246 doc="Subtask to rescale the variance of the template to the statistically expected level."
247 )
248 doSubtractBackground = lsst.pex.config.Field(
249 doc="Subtract the background fit when solving the kernel? "
250 "It is generally better to instead subtract the background in detectAndMeasure.",
251 dtype=bool,
252 default=False,
253 deprecated="Subtract the background in detectAndMeasure instead. Will be removed after v31.",
254 )
255 doApplyExternalCalibrations = lsst.pex.config.Field(
256 doc=(
257 "Replace science Exposure's calibration objects with those"
258 " in visitSummary. Ignored if `doApplyFinalizedPsf is True."
259 ),
260 dtype=bool,
261 default=False,
262 )
263 sourceSelector = lsst.pex.config.ConfigurableField(
264 target=ScienceSourceSelectorTask,
265 doc="Task to select sources to be used for PSF matching.",
266 )
267 fallbackSourceSelector = lsst.pex.config.ConfigurableField(
268 target=ScienceSourceSelectorTask,
269 doc="Task to select sources to be used for PSF matching."
270 "Used only if the kernel calculation fails and"
271 "`allowKernelSourceDetection` is set. The fallback source detection"
272 " will not include all of the same plugins as the original source "
273 " detection, so not all of the same flags can be used.",
274 )
275 detectionThreshold = lsst.pex.config.Field(
276 dtype=float,
277 default=10,
278 doc="Minimum signal to noise ratio of detected sources "
279 "to use for calculating the PSF matching kernel.",
280 deprecated="No longer used. Will be removed after v30"
281 )
282 detectionThresholdMax = lsst.pex.config.Field(
283 dtype=float,
284 default=500,
285 doc="Maximum signal to noise ratio of detected sources "
286 "to use for calculating the PSF matching kernel.",
287 deprecated="No longer used. Will be removed after v30"
288 )
289 restrictKernelEdgeSources = lsst.pex.config.Field(
290 dtype=bool,
291 default=True,
292 doc="Exclude sources close to the edge from the kernel calculation?"
293 )
294 maxKernelSources = lsst.pex.config.Field(
295 dtype=int,
296 default=1000,
297 doc="Maximum number of sources to use for calculating the PSF matching kernel."
298 "Set to -1 to disable."
299 )
300 minKernelSources = lsst.pex.config.Field(
301 dtype=int,
302 default=3,
303 doc="Minimum number of sources needed for calculating the PSF matching kernel."
304 )
305 excludeMaskPlanes = lsst.pex.config.ListField(
306 dtype=str,
307 default=("NO_DATA", "BAD", "SAT", "EDGE", "FAKE", "HIGH_VARIANCE"),
308 doc="Template mask planes to exclude when selecting sources for PSF matching.",
309 )
310 badMaskPlanes = lsst.pex.config.ListField(
311 dtype=str,
312 default=("NO_DATA", "BAD", "SAT", "EDGE"),
313 doc="Mask planes to interpolate over."
314 )
315 preserveTemplateMask = lsst.pex.config.ListField(
316 dtype=str,
317 default=("NO_DATA", "BAD", "HIGH_VARIANCE"),
318 doc="Mask planes from the template to propagate to the image difference."
319 )
320 renameTemplateMask = lsst.pex.config.ListField(
321 dtype=str,
322 default=("SAT", "INJECTED", "INJECTED_CORE",),
323 doc="Mask planes from the template to propagate to the image difference"
324 "with '_TEMPLATE' appended to the name."
325 )
326 preserveMaskPlanes = lsst.pex.config.ListField(
327 dtype=str,
328 default=("INJECTED", "INJECTED_CORE", "INJECTED_TEMPLATE", "INJECTED_CORE_TEMPLATE"),
329 doc="Mask planes to preserve without dilation when convolving the image.",
330 )
331 allowKernelSourceDetection = lsst.pex.config.Field(
332 dtype=bool,
333 default=False,
334 doc="Re-run source detection for kernel candidates if an error is"
335 " encountered while calculating the matching kernel."
336 )
337 image_type = lsst.pex.config.ChoiceField[str](
338 doc="Image type of the input template, the image outputs of this"
339 " task (difference and matchedTemplate, or scoreExposure), and the PSF"
340 " matching kernel.",
341 allowed={
342 "legacy": "Read and write lsst.afw.image.ExposureF and"
343 " lsst.afw.math.Kernel.",
344 "future": "Read and write lsst.images.DifferenceImage and"
345 " lsst.images.convolution_kernels.ConvolutionKernel.",
346 },
347 optional=False,
348 default="legacy",
349 )
351 def setDefaults(self):
352 self.makeKernel.kernel.name = "AL"
353 # Always include background fitting in the kernel fit,
354 # even if it is not subtracted
355 self.makeKernel.kernel.active.fitForBackground = True
356 self.makeKernel.kernel.active.spatialKernelOrder = 1
357 self.makeKernel.kernel.active.spatialBgOrder = 2
358 # Shared source selector settings
359 doSkySources = False # Do not include sky sources
360 doSignalToNoise = True # apply signal to noise filter
361 doUnresolved = True # apply star-galaxy separation
362 signalToNoiseMinimum = 10
363 signalToNoiseMaximum = 500
364 self.sourceSelector.doIsolated = True # apply isolated star selection
365 self.sourceSelector.doRequirePrimary = True # apply primary flag selection
366 self.sourceSelector.doUnresolved = doUnresolved
367 self.sourceSelector.doSkySources = doSkySources
368 self.sourceSelector.doSignalToNoise = doSignalToNoise
369 self.sourceSelector.signalToNoise.minimum = signalToNoiseMinimum
370 self.sourceSelector.signalToNoise.maximum = signalToNoiseMaximum
371 # The following two configs should not be necessary to be turned on for
372 # PSF-matching, and the fallback kernel source selection will fail if
373 # they are set since it does not run deblending.
374 self.fallbackSourceSelector.doIsolated = False # Do not apply isolated star selection
375 self.fallbackSourceSelector.doRequirePrimary = False # Do not apply primary flag selection
376 self.fallbackSourceSelector.doUnresolved = doUnresolved
377 self.fallbackSourceSelector.doSkySources = doSkySources
378 self.fallbackSourceSelector.doSignalToNoise = doSignalToNoise
379 self.fallbackSourceSelector.signalToNoise.minimum = signalToNoiseMinimum
380 self.fallbackSourceSelector.signalToNoise.maximum = signalToNoiseMaximum
382 def validate(self):
383 super().validate()
384 if self.doSubtractBackground and self.image_type == "future": 384 ↛ 385line 384 didn't jump to line 385 because the condition on line 384 was never true
385 raise lsst.pex.config.FieldValidationError(
386 AlardLuptonSubtractBaseConfig.doSubtractBackground, self,
387 "The background subtracted when solving the kernel is not recorded on"
388 " image_type='future' outputs; subtract the background in detectAndMeasure instead."
389 )
392class AlardLuptonSubtractConfig(AlardLuptonSubtractBaseConfig, lsst.pipe.base.PipelineTaskConfig,
393 pipelineConnections=AlardLuptonSubtractConnections):
394 mode = lsst.pex.config.ChoiceField(
395 dtype=str,
396 default="convolveTemplate",
397 allowed={"auto": "Choose which image to convolve at runtime.",
398 "convolveScience": "Only convolve the science image.",
399 "convolveTemplate": "Only convolve the template image."},
400 doc="Choose which image to convolve at runtime, or require that a specific image is convolved."
401 )
404class AlardLuptonSubtractTask(lsst.pipe.base.PipelineTask):
405 """Compute the image difference of a science and template image using
406 the Alard & Lupton (1998) algorithm.
407 """
408 ConfigClass = AlardLuptonSubtractConfig
409 _DefaultName = "alardLuptonSubtract"
410 usePreconvolution = False
411 """Whether this task preconvolves the science image with its own PSF
412 before kernel-matching. Subclasses that preconvolve override this to
413 `True`."""
415 futureImageOutputs = ("difference", "matchedTemplate")
416 """Names of the images in the results struct that are converted in future
417 mode (`tuple` [`str`]). Subclasses with different image outputs override
418 this."""
420 def __init__(self, **kwargs):
421 super().__init__(**kwargs)
422 self.makeSubtask("decorrelate")
423 self.makeSubtask("makeKernel")
424 self.makeSubtask("sourceSelector")
425 self.makeSubtask("fallbackSourceSelector")
426 if self.config.doScaleVariance:
427 self.makeSubtask("scaleVariance")
429 self.convolutionControl = lsst.afw.math.ConvolutionControl()
430 # Normalization is an extra, unnecessary, calculation and will result
431 # in mis-subtraction of the images if there are calibration errors.
432 self.convolutionControl.setDoNormalize(False)
433 self.convolutionControl.setDoCopyEdge(True)
435 def _applyExternalCalibrations(self, exposure, visitSummary):
436 """Replace calibrations (psf, and ApCorrMap) on this exposure with
437 external ones.".
439 Parameters
440 ----------
441 exposure : `lsst.afw.image.exposure.Exposure`
442 Input exposure to adjust calibrations.
443 visitSummary : `lsst.afw.table.ExposureCatalog`
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 exposure : `lsst.afw.image.exposure.Exposure`
451 Exposure with adjusted calibrations.
452 """
453 detectorId = exposure.info.getDetector().getId()
455 row = visitSummary.find(detectorId)
456 if row is None:
457 self.log.warning("Detector id %s not found in external calibrations catalog; "
458 "Using original calibrations.", detectorId)
459 else:
460 psf = row.getPsf()
461 apCorrMap = row.getApCorrMap()
462 if psf is None:
463 self.log.warning("Detector id %s has None for psf in "
464 "external calibrations catalog; Using original psf and aperture correction.",
465 detectorId)
466 elif apCorrMap is None:
467 self.log.warning("Detector id %s has None for apCorrMap in "
468 "external calibrations catalog; Using original psf and aperture correction.",
469 detectorId)
470 else:
471 exposure.setPsf(psf)
472 exposure.info.setApCorrMap(apCorrMap)
474 return exposure
476 def runQuantum(self, butlerQC, inputRefs, outputRefs):
477 inputs = butlerQC.get(inputRefs)
478 templateInfo = None
479 exposureRecord = None
480 photometricScaling = None
481 if self.config.image_type == "future":
482 template = inputs["template"]
483 dataId = butlerQC.quantum.dataId
484 # Read the `TemplateInfo` and the observation metadata from the
485 # template while it is still in the `DifferenceImage` format.
486 templateInfo = get_difference_image_provenance(template, "templates", "template", dataId)
487 exposureRecord = record_from_obs_info(template.obs_info, dataId["instrument"],
488 dataId["visit"], dataId.universe)
489 inputs["template"] = template.to_legacy()
490 photometricScaling = inputs["science"].photometric_scaling
491 inputs["science"] = inputs["science"].to_legacy()
492 if "inputPsfMatchingKernel" in inputs:
493 inputs["inputPsfMatchingKernel"] = inputs["inputPsfMatchingKernel"].to_legacy()
495 try:
496 results = self.run(**inputs)
497 except lsst.pipe.base.AlgorithmError as e:
498 error = lsst.pipe.base.AnnotatedPartialOutputsError.annotate(e, self, log=self.log)
499 # No partial outputs for butler to put, so nothing to convert.
500 raise error from e
502 if self.config.image_type == "future":
503 self.convert_outputs_to_future(results, exposureRecord, templateInfo=templateInfo,
504 photometricScaling=photometricScaling)
506 butlerQC.put(results, outputRefs)
508 def convert_outputs_to_future(self, results, exposureRecord, templateInfo=None,
509 photometricScaling=None):
510 """Convert the image outputs to `lsst.images` types.
512 Each image named in `futureImageOutputs` that is present on
513 ``results`` is replaced by an `lsst.images.DifferenceImage`, with the
514 PSF matching kernel attached to it. ``psfMatchingKernel`` is replaced
515 by the same kernel as an
516 `lsst.images.convolution_kernels.ImageBasisConvolutionKernel`.
518 Parameters
519 ----------
520 results : `lsst.pipe.base.Struct`
521 Output struct to read and modify in place.
522 exposureRecord : `lsst.daf.butler.DimensionRecord`
523 The ``exposure`` record of the observation, which supplies the
524 observation metadata recorded on each converted output.
525 templateInfo : `list` [`lsst.images.DifferenceImageTemplateInfo`], \
526 optional
527 Record of the coadds that went into the template, taken from the
528 template this task subtracted.
529 photometricScaling : `lsst.images.fields.BaseField`, optional
530 Photometric scaling of the science image.
532 Raises
533 ------
534 TypeError
535 Raised if the PSF matching kernel cannot be stored on an
536 `lsst.images.DifferenceImage`.
537 RuntimeError
538 Raised if an output has pixels set in a mask plane that
539 `lsst.images.get_legacy_difference_image_mask_planes` does not
540 map.
542 Notes
543 -----
544 This uses the default legacy mask plane map of
545 `lsst.images.get_legacy_difference_image_mask_planes` and adds any
546 planes from `~lsst.images.get_legacy_optional_mask_planes` that are set
547 for any pixels in the image.
548 """
549 # The science image and every output built from it share the science
550 # pixel grid, so its detector applies to all of them.
551 science = results.matchedScience
552 try:
553 kernel = ImageBasisConvolutionKernel.from_legacy(results.psfMatchingKernel)
554 except (TypeError, ValueError, AttributeError) as e:
555 raise TypeError("The supplied PSF matching kernel cannot be attached to an"
556 f" lsst.images.DifferenceImage: {e}") from e
557 for name in self.futureImageOutputs:
558 exposure = getattr(results, name, None)
559 if exposure is None:
560 continue
561 if exposure.getDetector() is None:
562 # Fall back on the science detector for images without a
563 # detector of their own (e.g. matchedTemplate)
564 exposure.setDetector(science.getDetector())
565 image = DifferenceImage.from_legacy(
566 exposure,
567 unit=u.nJy,
568 exposure_record=exposureRecord,
569 )
570 if kernel is not None: 570 ↛ 572line 570 didn't jump to line 572 because the condition on line 570 was always true
571 image.kernel = kernel
572 if templateInfo:
573 image.templates = templateInfo
574 if photometricScaling is not None:
575 image.photometric_scaling = photometricScaling
576 setattr(results, name, image)
577 results.psfMatchingKernel = kernel
579 @timeMethod
580 def run(self, template, science, sources, visitSummary=None):
581 """PSF match, subtract, and decorrelate two images.
583 Parameters
584 ----------
585 template : `lsst.afw.image.ExposureF`
586 Template exposure, warped to match the science exposure.
587 science : `lsst.afw.image.ExposureF`
588 Science exposure to subtract from the template.
589 sources : `lsst.afw.table.SourceCatalog`
590 Identified sources on the science exposure. This catalog is used to
591 select sources in order to perform the AL PSF matching on stamp
592 images around them.
593 visitSummary : `lsst.afw.table.ExposureCatalog`, optional
594 Exposure catalog with external calibrations to be applied. Catalog
595 uses the detector id for the catalog id, sorted on id for fast
596 lookup.
598 Returns
599 -------
600 results : `lsst.pipe.base.Struct`
601 ``difference`` : `lsst.afw.image.ExposureF`
602 Result of subtracting template and science.
603 ``matchedTemplate`` : `lsst.afw.image.ExposureF`
604 Warped and PSF-matched template exposure.
605 ``backgroundModel`` : `lsst.afw.math.Function2D`
606 Background model that was fit while solving for the
607 PSF-matching kernel
608 ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
609 Kernel used to PSF-match the convolved image.
610 ``kernelSources` : `lsst.afw.table.SourceCatalog`
611 Sources from the input catalog that were used to construct the
612 PSF-matching kernel.
613 """
614 self._prepareInputs(template, science, visitSummary=visitSummary)
616 convolveTemplate = self.chooseConvolutionMethod(template, science)
617 self.matchedPsfSize = self.sciencePsfSize if convolveTemplate else self.templatePsfSize
619 kernelResult = self.runMakeKernel(template, science, sources=sources,
620 convolveTemplate=convolveTemplate,
621 runSourceDetection=False)
623 if self.config.doSubtractBackground:
624 backgroundModel = kernelResult.backgroundModel
625 else:
626 backgroundModel = None
627 if convolveTemplate:
628 subtractResults = self.runConvolveTemplate(template, science, kernelResult.psfMatchingKernel,
629 backgroundModel=backgroundModel)
630 else:
631 subtractResults = self.runConvolveScience(template, science, kernelResult.psfMatchingKernel,
632 backgroundModel=backgroundModel)
633 subtractResults.kernelSources = kernelResult.kernelSources
635 metrics = computeDifferenceImageMetrics(science, subtractResults.difference, sources)
637 self.metadata["differenceFootprintRatioMean"] = metrics.differenceFootprintRatioMean
638 self.metadata["differenceFootprintRatioStdev"] = metrics.differenceFootprintRatioStdev
639 self.metadata["differenceFootprintSkyRatioMean"] = metrics.differenceFootprintSkyRatioMean
640 self.metadata["differenceFootprintSkyRatioStdev"] = metrics.differenceFootprintSkyRatioStdev
641 self.log.info("Mean, stdev of ratio of difference to science "
642 "pixels in star footprints: %5.4f, %5.4f",
643 self.metadata["differenceFootprintRatioMean"],
644 self.metadata["differenceFootprintRatioStdev"])
646 return subtractResults
648 def chooseConvolutionMethod(self, template, science):
649 """Determine whether the template should be convolved with the PSF
650 matching kernel.
652 Parameters
653 ----------
654 template : `lsst.afw.image.ExposureF`
655 Template exposure, warped to match the science exposure.
656 science : `lsst.afw.image.ExposureF`
657 Science exposure to subtract from the template.
659 Returns
660 -------
661 convolveTemplate : `bool`
662 Convolve the template to match the two images?
664 Raises
665 ------
666 RuntimeError
667 If an unsupported convolution mode is supplied.
668 """
669 if self.usePreconvolution: 669 ↛ 670line 669 didn't jump to line 670 because the condition on line 669 was never true
670 raise RuntimeError("Choosing a convolution method is incompatible with preconvolution!")
671 if self.config.mode == "auto":
672 convolveTemplate = _shapeTest(template,
673 science,
674 fwhmExposureBuffer=self.config.makeKernel.fwhmExposureBuffer,
675 fwhmExposureGrid=self.config.makeKernel.fwhmExposureGrid)
676 if convolveTemplate:
677 if self.sciencePsfSize < self.templatePsfSize: 677 ↛ 678line 677 didn't jump to line 678 because the condition on line 677 was never true
678 self.log.info("Average template PSF size is greater, "
679 "but science PSF greater in one dimension: convolving template image.")
680 else:
681 self.log.info("Science PSF size is greater: convolving template image.")
682 else:
683 self.log.info("Template PSF size is greater: convolving science image.")
684 elif self.config.mode == "convolveTemplate":
685 self.log.info("`convolveTemplate` is set: convolving template image.")
686 convolveTemplate = True
687 elif self.config.mode == "convolveScience": 687 ↛ 691line 687 didn't jump to line 691 because the condition on line 687 was always true
688 self.log.info("`convolveScience` is set: convolving science image.")
689 convolveTemplate = False
690 else:
691 raise RuntimeError(f"Cannot handle AlardLuptonSubtract mode: {self.config.mode}")
692 return convolveTemplate
694 def runMakeKernel(self, template, science, sources=None, convolveTemplate=True, runSourceDetection=False):
695 """Construct the PSF-matching kernel. Not used for preconvolution.
697 Parameters
698 ----------
699 template : `lsst.afw.image.ExposureF`
700 Template exposure, warped to match the science exposure.
701 science : `lsst.afw.image.ExposureF`
702 Science exposure to subtract from the template.
703 sources : `lsst.afw.table.SourceCatalog`
704 Identified sources on the science exposure. This catalog is used to
705 select sources in order to perform the AL PSF matching on stamp
706 images around them.
707 Not used if ``runSourceDetection`` is set.
708 convolveTemplate : `bool`, optional
709 Construct the matching kernel to convolve the template?
710 runSourceDetection : `bool`, optional
711 Run a minimal version of source detection to determine kernel
712 candidates? If False, a source list to select kernel candidates
713 from must be supplied.
715 Returns
716 -------
717 results : `lsst.pipe.base.Struct`
718 ``backgroundModel`` : `lsst.afw.math.Function2D`
719 Background model that was fit while solving for the
720 PSF-matching kernel
721 ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
722 Kernel used to PSF-match the convolved image.
723 ``kernelSources` : `lsst.afw.table.SourceCatalog`
724 Sources from the input catalog that were used to construct the
725 PSF-matching kernel.
726 """
727 if self.usePreconvolution: 727 ↛ 728line 727 didn't jump to line 728 because the condition on line 727 was never true
728 raise RuntimeError("Incorrect matching kernel calculation configured. "
729 "`runMakeKernel` can't be called if `usePreconvolution` is set.")
730 if convolveTemplate:
731 reference = template
732 target = science
733 referenceFwhmPix = self.templatePsfSize
734 targetFwhmPix = self.sciencePsfSize
735 else:
736 reference = science
737 target = template
738 referenceFwhmPix = self.sciencePsfSize
739 targetFwhmPix = self.templatePsfSize
740 try:
741 if runSourceDetection:
742 kernelSources = self.runKernelSourceDetection(template, science)
743 else:
744 kernelSources = self._sourceSelector(template, science, sources)
745 kernelResult = self.makeKernel.run(reference, target, kernelSources,
746 preconvolved=False,
747 templateFwhmPix=referenceFwhmPix,
748 scienceFwhmPix=targetFwhmPix)
749 except (RuntimeError, lsst.pex.exceptions.Exception) as e:
750 self.log.warning("Failed to match template. Checking coverage")
751 # Raise NoWorkFound if template fraction is insufficient
752 checkTemplateIsSufficient(template[science.getBBox()], science, self.log,
753 self.config.minTemplateFractionForExpectedSuccess,
754 exceptionMessage="Template coverage lower than expected to succeed."
755 f" Failure is tolerable: {e}")
756 # checkTemplateIsSufficient did not raise NoWorkFound, so raise original exception
757 raise e
759 return lsst.pipe.base.Struct(backgroundModel=kernelResult.backgroundModel,
760 psfMatchingKernel=kernelResult.psfMatchingKernel,
761 kernelSources=kernelSources)
763 def runKernelSourceDetection(self, template, science):
764 """Run detection on the science image and use the template mask plane
765 to reject candidate sources.
767 Parameters
768 ----------
769 template : `lsst.afw.image.ExposureF`
770 Template exposure, warped to match the science exposure.
771 science : `lsst.afw.image.ExposureF`
772 Science exposure to subtract from the template.
774 Returns
775 -------
776 kernelSources : `lsst.afw.table.SourceCatalog`
777 Sources from the input catalog to use to construct the
778 PSF-matching kernel.
779 """
780 kernelSize = self.makeKernel.makeKernelBasisList(
781 self.templatePsfSize, self.matchedPsfSize)[0].getWidth()
782 sources = self.makeKernel.makeCandidateList(template, science, kernelSize,
783 candidateList=None,
784 sigma=gaussian_fwhm_to_sigma*self.sciencePsfSize)
785 return self._sourceSelector(template, science, sources, fallback=True)
787 def runConvolveTemplate(self, template, science, psfMatchingKernel, backgroundModel=None):
788 """Convolve the template image with a PSF-matching kernel and subtract
789 from the science image.
791 Parameters
792 ----------
793 template : `lsst.afw.image.ExposureF`
794 Template exposure, warped to match the science exposure.
795 science : `lsst.afw.image.ExposureF`
796 Science exposure to subtract from the template.
797 psfMatchingKernel : `lsst.afw.math.Kernel`
798 Kernel to be used to PSF-match the science image to the template.
799 backgroundModel : `lsst.afw.math.Function2D`, optional
800 Background model that was fit while solving for the PSF-matching
801 kernel.
803 Returns
804 -------
805 results : `lsst.pipe.base.Struct`
807 ``difference`` : `lsst.afw.image.ExposureF`
808 Result of subtracting template and science.
809 ``matchedTemplate`` : `lsst.afw.image.ExposureF`
810 Warped and PSF-matched template exposure.
811 ``backgroundModel`` : `lsst.afw.math.Function2D`
812 Background model that was fit while solving for the PSF-matching kernel
813 ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
814 Kernel used to PSF-match the template to the science image.
815 """
816 self.metadata["convolvedExposure"] = "Template"
818 matchedTemplate = self._convolveExposure(template, psfMatchingKernel,
819 self.convolutionControl,
820 bbox=science.getBBox(),
821 psf=science.psf,
822 photoCalib=science.photoCalib)
824 difference = _subtractImages(science, matchedTemplate, backgroundModel=backgroundModel)
825 correctedExposure = self.finalize(template, science, difference,
826 psfMatchingKernel,
827 templateMatched=True)
829 return lsst.pipe.base.Struct(difference=correctedExposure,
830 matchedTemplate=matchedTemplate,
831 matchedScience=science,
832 backgroundModel=backgroundModel,
833 psfMatchingKernel=psfMatchingKernel)
835 def runConvolveScience(self, template, science, psfMatchingKernel, backgroundModel=None):
836 """Convolve the science image with a PSF-matching kernel and subtract
837 the template image.
839 Parameters
840 ----------
841 template : `lsst.afw.image.ExposureF`
842 Template exposure, warped to match the science exposure.
843 science : `lsst.afw.image.ExposureF`
844 Science exposure to subtract from the template.
845 psfMatchingKernel : `lsst.afw.math.Kernel`
846 Kernel to be used to PSF-match the science image to the template.
847 backgroundModel : `lsst.afw.math.Function2D`, optional
848 Background model that was fit while solving for the PSF-matching
849 kernel.
851 Returns
852 -------
853 results : `lsst.pipe.base.Struct`
855 ``difference`` : `lsst.afw.image.ExposureF`
856 Result of subtracting template and science.
857 ``matchedTemplate`` : `lsst.afw.image.ExposureF`
858 Warped template exposure. Note that in this case, the template
859 is not PSF-matched to the science image.
860 ``backgroundModel`` : `lsst.afw.math.Function2D`
861 Background model that was fit while solving for the PSF-matching kernel
862 ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
863 Kernel used to PSF-match the science image to the template.
864 """
865 self.metadata["convolvedExposure"] = "Science"
866 bbox = science.getBBox()
868 kernelImage = lsst.afw.image.ImageD(psfMatchingKernel.getDimensions())
869 xcen, ycen = bbox.getCenter()
870 norm = psfMatchingKernel.computeImage(kernelImage, doNormalize=False, x=xcen, y=ycen)
872 matchedScience = self._convolveExposure(science, psfMatchingKernel,
873 self.convolutionControl,
874 psf=template.psf)
876 # Place back on native photometric scale
877 matchedScience.maskedImage /= norm
878 matchedTemplate = template.clone()[bbox]
879 matchedTemplate.setPhotoCalib(science.photoCalib)
881 if backgroundModel is not None:
882 # We must invert the background model if the matching kernel is solved for the science image.
883 invertedBackground = backgroundModel.clone()
884 invertedBackground.setParameters([-p for p in backgroundModel.getParameters()])
885 backgroundModel = invertedBackground
887 difference = _subtractImages(matchedScience, matchedTemplate, backgroundModel=backgroundModel)
889 correctedExposure = self.finalize(template, science, difference,
890 psfMatchingKernel,
891 templateMatched=False)
893 return lsst.pipe.base.Struct(difference=correctedExposure,
894 matchedTemplate=matchedTemplate,
895 matchedScience=matchedScience,
896 backgroundModel=backgroundModel,
897 psfMatchingKernel=psfMatchingKernel)
899 def finalize(self, template, science, difference, kernel,
900 templateMatched=True,
901 preConvMode=False,
902 preConvKernel=None,
903 spatiallyVarying=False):
904 """Decorrelate the difference image to undo the noise correlations
905 caused by convolution.
907 Parameters
908 ----------
909 template : `lsst.afw.image.ExposureF`
910 Template exposure, warped to match the science exposure.
911 science : `lsst.afw.image.ExposureF`
912 Science exposure to subtract from the template.
913 difference : `lsst.afw.image.ExposureF`
914 Result of subtracting template and science.
915 kernel : `lsst.afw.math.Kernel`
916 An (optionally spatially-varying) PSF matching kernel
917 templateMatched : `bool`, optional
918 Was the template PSF-matched to the science image?
919 preConvMode : `bool`, optional
920 Was the science image preconvolved with its own PSF
921 before PSF matching the template?
922 preConvKernel : `lsst.afw.detection.Psf`, optional
923 If not `None`, then the science image was pre-convolved with
924 (the reflection of) this kernel. Must be normalized to sum to 1.
925 spatiallyVarying : `bool`, optional
926 Compute the decorrelation kernel spatially varying across the image?
928 Returns
929 -------
930 correctedExposure : `lsst.afw.image.ExposureF`
931 The decorrelated image difference.
932 """
933 if self.config.doDecorrelation:
934 self.log.info("Decorrelating image difference.")
935 # We have cleared the template mask plane, so copy the mask plane of
936 # the image difference so that we can calculate correct statistics
937 # during decorrelation
938 correctedExposure = self.decorrelate.run(science, template[science.getBBox()], difference, kernel,
939 templateMatched=templateMatched,
940 preConvMode=preConvMode,
941 preConvKernel=preConvKernel,
942 spatiallyVarying=spatiallyVarying).correctedExposure
943 else:
944 self.log.info("NOT decorrelating image difference.")
945 correctedExposure = difference
946 return correctedExposure
948 def _calculateMagLim(self, exposure, nsigma=5.0, fallbackPsfSize=None):
949 """Calculate the point-source limiting magnitude of an exposure.
951 The limit is the flux of a point source detected at ``nsigma`` in
952 the median pixel noise, integrated over the effective area of the
953 PSF, so it does not depend on the units of the pixels.
955 Parameters
956 ----------
957 exposure : `lsst.afw.image.Exposure`
958 The exposure to calculate the limiting magnitude for.
959 nsigma : `float`, optional
960 The detection threshold in sigma.
961 fallbackPsfSize : `float`, optional
962 PSF FWHM in pixels to use if the exposure PSF cannot be
963 evaluated.
965 Returns
966 -------
967 maglim : `float`
968 The limiting AB magnitude, or `numpy.nan` if the exposure has no
969 photometric calibration, no usable pixels, or no PSF.
970 """
971 if exposure.photoCalib is None: 971 ↛ 972line 971 didn't jump to line 972 because the condition on line 971 was never true
972 return np.nan
973 try:
974 psfArea = computePSFNoiseEquivalentArea(exposure.getPsf())
975 except (lsst.pex.exceptions.InvalidParameterError,
976 afwDetection.InvalidPsfError,
977 lsst.pex.exceptions.RangeError):
978 if fallbackPsfSize is None:
979 self.log.info("Unable to evaluate PSF, setting maglim to nan")
980 return np.nan
981 self.log.info("Unable to evaluate PSF, using fallback FWHM %f", fallbackPsfSize)
982 # Effective area of a Gaussian PSF with this FWHM.
983 psfArea = 4*np.pi*(fallbackPsfSize/np.sqrt(8*np.log(2)))**2
984 mask = exposure.mask
985 good = (mask.array & mask.getPlaneBitMask(self.config.badMaskPlanes)) == 0
986 variance = exposure.variance.array[good]
987 variance = variance[np.isfinite(variance) & (variance > 0)]
988 if variance.size == 0:
989 return np.nan
990 fluxLim = nsigma*np.sqrt(np.median(variance)*psfArea)
991 return exposure.photoCalib.instFluxToMagnitude(fluxLim)
993 @staticmethod
994 def _validateExposures(template, science):
995 """Check that the WCS of the two Exposures match, the template bbox
996 contains the science bbox, and that the bands match.
998 Parameters
999 ----------
1000 template : `lsst.afw.image.ExposureF`
1001 Template exposure, warped to match the science exposure.
1002 science : `lsst.afw.image.ExposureF`
1003 Science exposure to subtract from the template.
1005 Raises
1006 ------
1007 AssertionError
1008 Raised if the WCS of the template is not equal to the science WCS,
1009 if the science image is not fully contained in the template
1010 bounding box, or if the bands do not match.
1011 """
1012 assert template.wcs == science.wcs, \
1013 "Template and science exposure WCS are not identical."
1014 templateBBox = template.getBBox()
1015 scienceBBox = science.getBBox()
1016 assert science.filter.bandLabel == template.filter.bandLabel, \
1017 "Science and template exposures have different bands: %s, %s" % \
1018 (science.filter, template.filter)
1020 assert templateBBox.contains(scienceBBox), \
1021 "Template bbox does not contain all of the science image."
1023 def _convolveExposure(self, exposure, kernel, convolutionControl,
1024 bbox=None,
1025 psf=None,
1026 photoCalib=None,
1027 interpolateBadMaskPlanes=False,
1028 ):
1029 """Convolve an exposure with the given kernel.
1031 Parameters
1032 ----------
1033 exposure : `lsst.afw.Exposure`
1034 exposure to convolve.
1035 kernel : `lsst.afw.math.LinearCombinationKernel`
1036 PSF matching kernel computed in the ``makeKernel`` subtask.
1037 convolutionControl : `lsst.afw.math.ConvolutionControl`
1038 Configuration for convolve algorithm.
1039 bbox : `lsst.geom.Box2I`, optional
1040 Bounding box to trim the convolved exposure to.
1041 psf : `lsst.afw.detection.Psf`, optional
1042 Point spread function (PSF) to set for the convolved exposure.
1043 photoCalib : `lsst.afw.image.PhotoCalib`, optional
1044 Photometric calibration of the convolved exposure.
1045 interpolateBadMaskPlanes : `bool`, optional
1046 If set, interpolate over mask planes specified in
1047 ``config.badMaskPlanes`` before convolving the image.
1049 Returns
1050 -------
1051 convolvedExp : `lsst.afw.Exposure`
1052 The convolved image.
1053 """
1054 convolvedExposure = exposure.clone()
1055 if psf is not None:
1056 convolvedExposure.setPsf(psf)
1057 if photoCalib is not None:
1058 convolvedExposure.setPhotoCalib(photoCalib)
1059 if interpolateBadMaskPlanes and self.config.badMaskPlanes is not None:
1060 nInterp = _interpolateImage(convolvedExposure.maskedImage,
1061 self.config.badMaskPlanes)
1062 self.metadata["nInterpolated"] = nInterp
1064 # Snapshot the footprints of mask planes that must not be dilated by
1065 # the convolution.
1066 preservePlanes = [mp for mp in self.config.preserveMaskPlanes
1067 if mp in convolvedExposure.mask.getMaskPlaneDict()]
1068 maskResetDict = {
1069 mp: (convolvedExposure.mask.array
1070 & convolvedExposure.mask.getPlaneBitMask(mp)) > 0
1071 for mp in preservePlanes
1072 }
1074 convolvedImage = lsst.afw.image.MaskedImageF(convolvedExposure.getBBox())
1075 lsst.afw.math.convolve(convolvedImage, convolvedExposure.maskedImage, kernel, convolutionControl)
1076 convolvedExposure.setMaskedImage(convolvedImage)
1078 # Undo the convolution's dilation of the preserved planes: clear the
1079 # dilated bits, then restore each mask plane for the pixels that were
1080 # previously set.
1081 self._clearMask(convolvedExposure.mask, clearMaskPlanes=preservePlanes)
1082 for maskPlane, maskSetPixels in maskResetDict.items():
1083 bit = convolvedExposure.mask.getPlaneBitMask(maskPlane)
1084 convolvedExposure.mask.array[maskSetPixels] |= bit
1086 if bbox is None:
1087 return convolvedExposure
1088 else:
1089 return convolvedExposure[bbox]
1091 def _sourceSelector(self, template, science, sources, fallback=False):
1092 """Select sources from a catalog that meet the selection criteria.
1093 The selection criteria include any configured parameters of the
1094 `sourceSelector` subtask, as well as checking the science and template
1095 mask planes.
1097 Parameters
1098 ----------
1099 template : `lsst.afw.image.ExposureF`
1100 Template exposure, warped to match the science exposure.
1101 science : `lsst.afw.image.ExposureF`
1102 Science exposure to subtract from the template.
1103 sources : `lsst.afw.table.SourceCatalog`
1104 Input source catalog to select sources from.
1105 fallback : `bool`, optional
1106 Switch indicating the source selector is being called after
1107 running the fallback source detection subtask, which does not run a
1108 full set of measurement plugins and can't use the same settings for
1109 the source selector.
1111 Returns
1112 -------
1113 kernelSources : `lsst.afw.table.SourceCatalog`
1114 The input source catalog, with flagged and low signal-to-noise
1115 sources removed and footprints added.
1117 Raises
1118 ------
1119 InsufficientKernelSourcesError
1120 An AlgorithmError that is raised if there are not enough PSF
1121 candidates to construct the PSF matching kernel.
1122 """
1123 if fallback:
1124 selected = self.fallbackSourceSelector.selectSources(sources).selected
1125 else:
1126 selected = self.sourceSelector.selectSources(sources).selected
1127 # It is OK to use just self.matchedPsfSize for the science PSF here,
1128 # since we are just using it to calculate the size of the matching
1129 # kernel.
1130 kSize = self.makeKernel.makeKernelBasisList(self.templatePsfSize, self.matchedPsfSize)[0].getWidth()
1131 selectSources = sources[selected].copy(deep=True)
1132 # Set the footprints, to be used in `makeKernel` and `checkMask`.
1133 kernelSources = setSourceFootprints(selectSources, kernelSize=kSize)
1134 bbox = science.getBBox()
1135 if self.usePreconvolution:
1136 # Exclude a wider buffer around the edge of the image to
1137 # account for an extra convolution.
1138 bbox.grow(-kSize)
1139 if self.config.restrictKernelEdgeSources:
1140 bbox.grow(-kSize)
1141 # Remove sources that land on masked pixels
1142 scienceSelected = checkMask(science.mask[bbox], kernelSources, self.config.excludeMaskPlanes)
1143 templateSelected = checkMask(template.mask[bbox], kernelSources, self.config.excludeMaskPlanes)
1144 maskSelected = scienceSelected & templateSelected
1145 kernelSources = kernelSources[maskSelected].copy(deep=True)
1146 # Trim kernelSources if they exceed ``maxKernelSources``.
1147 # Keep the highest signal-to-noise sources of those selected.
1148 if (len(kernelSources) > self.config.maxKernelSources) & (self.config.maxKernelSources > 0):
1149 signalToNoise = kernelSources.getPsfInstFlux()/kernelSources.getPsfInstFluxErr()
1150 indices = np.argsort(signalToNoise)
1151 indices = indices[-self.config.maxKernelSources:]
1152 selected = np.zeros(len(kernelSources), dtype=bool)
1153 selected[indices] = True
1154 kernelSources = kernelSources[selected].copy(deep=True)
1156 self.log.info("%i/%i=%.1f%% of sources selected for PSF matching from the input catalog",
1157 len(kernelSources), len(sources), 100*len(kernelSources)/len(sources))
1158 if len(kernelSources) < self.config.minKernelSources:
1159 self.log.error("Too few sources to calculate the PSF matching kernel: "
1160 "%i selected but %i needed for the calculation.",
1161 len(kernelSources), self.config.minKernelSources)
1162 if self.config.allowKernelSourceDetection and not fallback: 1162 ↛ 1166line 1162 didn't jump to line 1166 because the condition on line 1162 was never true
1163 # The fallback source detection pipeline calls this method, so
1164 # allowing source detection in that case would create an endless
1165 # loop
1166 kernelSources = self.runKernelSourceDetection(template, science)
1167 else:
1168 raise InsufficientKernelSourcesError(nSources=len(kernelSources),
1169 nRequired=self.config.minKernelSources)
1171 self.metadata["nPsfSources"] = len(kernelSources)
1173 return kernelSources
1175 def _prepareInputs(self, template, science, visitSummary=None):
1176 """Perform preparatory calculations common to all Alard&Lupton Tasks.
1178 Parameters
1179 ----------
1180 template : `lsst.afw.image.ExposureF`
1181 Template exposure, warped to match the science exposure. The
1182 variance plane of the template image is modified in place.
1183 science : `lsst.afw.image.ExposureF`
1184 Science exposure to subtract from the template. The variance plane
1185 of the science image is modified in place.
1186 visitSummary : `lsst.afw.table.ExposureCatalog`, optional
1187 Exposure catalog with external calibrations to be applied. Catalog
1188 uses the detector id for the catalog id, sorted on id for fast
1189 lookup.
1190 """
1191 self._validateExposures(template, science)
1192 if visitSummary is not None: 1192 ↛ 1193line 1192 didn't jump to line 1193 because the condition on line 1192 was never true
1193 self._applyExternalCalibrations(science, visitSummary=visitSummary)
1194 templateCoverageFraction = checkTemplateIsSufficient(
1195 template[science.getBBox()], science, self.log,
1196 requiredTemplateFraction=self.config.requiredTemplateFraction,
1197 exceptionMessage="Not attempting subtraction. To force subtraction,"
1198 " set config requiredTemplateFraction=0"
1199 )
1200 self.metadata["templateCoveragePercent"] = 100*templateCoverageFraction
1202 if self.config.doScaleVariance:
1203 # Scale the variance of the science image before
1204 # convolution, subtraction, or decorrelation so that it has the
1205 # correct ratio. Note that the template variance is scaled
1206 # independently in ``GetTemplateTask``.
1207 sciVarFactor = self.scaleVariance.run(science.maskedImage)
1208 self.log.info("Science variance scaling factor: %.2f", sciVarFactor)
1209 self.metadata["scaleScienceVarianceFactor"] = sciVarFactor
1211 # Erase existing detection mask planes.
1212 # We don't want the detection mask from the science image
1213 self.updateMasks(template, science)
1215 # Calling getPsfFwhm on template.psf fails on some rare occasions when
1216 # the template has no input exposures at the average position of the
1217 # stars. So we try getPsfFwhm first on template, and if that fails we
1218 # evaluate the PSF on a grid specified by fwhmExposure* fields.
1219 # To keep consistent definitions for PSF size on the template and
1220 # science images, we use the same method for both.
1221 # In the try block below, we catch two exceptions:
1222 # 1. InvalidParameterError, in case the point where we are evaluating
1223 # the PSF lands in a gap in the template.
1224 # 2. RangeError, in case the template coverage is so poor that we end
1225 # up near a region with no data.
1226 try:
1227 self.templatePsfSize = getPsfFwhm(template.psf)
1228 self.sciencePsfSize = getPsfFwhm(science.psf)
1229 except lsst.pex.exceptions.Exception:
1230 # Catch a broad range of exceptions, since some are C++ only
1231 # Catching:
1232 # - lsst::geom::SingularTransformException
1233 # - lsst.pex.exceptions.InvalidParameterError
1234 # - lsst.pex.exceptions.RangeError
1235 self.log.info("Unable to evaluate PSF at the average position. "
1236 "Evaluting PSF on a grid of points."
1237 )
1238 self.templatePsfSize = evaluateMeanPsfFwhm(
1239 template,
1240 fwhmExposureBuffer=self.config.makeKernel.fwhmExposureBuffer,
1241 fwhmExposureGrid=self.config.makeKernel.fwhmExposureGrid
1242 )
1243 self.sciencePsfSize = evaluateMeanPsfFwhm(
1244 science,
1245 fwhmExposureBuffer=self.config.makeKernel.fwhmExposureBuffer,
1246 fwhmExposureGrid=self.config.makeKernel.fwhmExposureGrid
1247 )
1248 self.log.info("Science PSF FWHM: %f pixels", self.sciencePsfSize)
1249 self.log.info("Template PSF FWHM: %f pixels", self.templatePsfSize)
1250 self.metadata["sciencePsfSize"] = self.sciencePsfSize
1251 self.metadata["templatePsfSize"] = self.templatePsfSize
1253 # Calculate estimated image depths, i.e., limiting magnitudes
1254 maglim_science = self._calculateMagLim(science, fallbackPsfSize=self.sciencePsfSize)
1255 if np.isnan(maglim_science): 1255 ↛ 1256line 1255 didn't jump to line 1256 because the condition on line 1255 was never true
1256 self.log.warning("Limiting magnitude of the science image is NaN!")
1257 fluxlim_science = (maglim_science*u.ABmag).to_value(u.nJy)
1258 maglim_template = self._calculateMagLim(template, fallbackPsfSize=self.templatePsfSize)
1259 if np.isnan(maglim_template): 1259 ↛ 1260line 1259 didn't jump to line 1260 because the condition on line 1259 was never true
1260 self.log.info("Cannot evaluate template limiting mag; adopting science limiting mag for diffim")
1261 maglim_diffim = maglim_science
1262 else:
1263 fluxlim_template = (maglim_template*u.ABmag).to_value(u.nJy)
1264 maglim_diffim = (np.sqrt(fluxlim_science**2 + fluxlim_template**2)*u.nJy).to(u.ABmag).value
1265 self.metadata["scienceLimitingMagnitude"] = maglim_science
1266 self.metadata["templateLimitingMagnitude"] = maglim_template
1267 self.metadata["diffimLimitingMagnitude"] = maglim_diffim
1269 def updateMasks(self, template, science):
1270 """Update the science and template mask planes before differencing.
1272 Parameters
1273 ----------
1274 template : `lsst.afw.image.Exposure`
1275 Template exposure, warped to match the science exposure.
1276 The template mask planes will be erased, except for a few specified
1277 in the task config.
1278 science : `lsst.afw.image.Exposure`
1279 Science exposure to subtract from the template.
1280 The DETECTED and DETECTED_NEGATIVE mask planes of the science image
1281 will be erased.
1282 """
1283 self._clearMask(science.mask, clearMaskPlanes=["DETECTED", "DETECTED_NEGATIVE"])
1285 # We will clear ALL template mask planes, except for those specified
1286 # via the `preserveTemplateMask` config. Mask planes specified via
1287 # the `renameTemplateMask` config will be copied to new planes with
1288 # "_TEMPLATE" appended to their names, and the original mask plane will
1289 # be cleared.
1290 renameMaskPlanes = [mp for mp in self.config.renameTemplateMask
1291 if mp in template.mask.getMaskPlaneDict().keys()]
1292 # afw's mask planes are global, so if the "*_TEMPLATE" mask planes were
1293 # defined anywhere previously they would now also be cleared unless we
1294 # exclude them.
1295 renamedMaskPlanes = {maskPlane + "_TEMPLATE" for maskPlane in renameMaskPlanes}
1296 clearMaskPlanes = [mp for mp in template.mask.getMaskPlaneDict().keys()
1297 if mp not in self.config.preserveTemplateMask
1298 and mp not in renamedMaskPlanes]
1300 # propagate the mask plane related to Fake source injection
1301 # NOTE: the fake source injection sets FAKE plane, but it should be INJECTED
1302 # NOTE: This can be removed in DM-40796
1303 if "FAKE" in science.mask.getMaskPlaneDict().keys():
1304 self.log.info("Adding injected mask plane to science image")
1305 self._renameMaskPlanes(science.mask, "FAKE", "INJECTED")
1306 if "FAKE" in template.mask.getMaskPlaneDict().keys():
1307 self.log.info("Adding injected mask plane to template image")
1308 self._renameMaskPlanes(template.mask, "FAKE", "INJECTED_TEMPLATE")
1309 if "INJECTED" in renameMaskPlanes: 1309 ↛ 1311line 1309 didn't jump to line 1311 because the condition on line 1309 was always true
1310 renameMaskPlanes.remove("INJECTED")
1311 if "INJECTED_TEMPLATE" in clearMaskPlanes: 1311 ↛ 1312line 1311 didn't jump to line 1312 because the condition on line 1311 was never true
1312 clearMaskPlanes.remove("INJECTED_TEMPLATE")
1314 for maskPlane in renameMaskPlanes:
1315 self._renameMaskPlanes(template.mask, maskPlane, maskPlane + "_TEMPLATE")
1316 self._clearMask(template.mask, clearMaskPlanes=clearMaskPlanes)
1318 @staticmethod
1319 def _renameMaskPlanes(mask, maskPlane, newMaskPlane):
1320 """Rename a mask plane by adding the new name and copying the data.
1322 Parameters
1323 ----------
1324 mask : `lsst.afw.image.Mask`
1325 The mask image to update in place.
1326 maskPlane : `str`
1327 The name of the existing mask plane to copy.
1328 newMaskPlane : `str`
1329 The new name of the mask plane that will be added.
1330 If the mask plane already exists, it will be updated in place.
1331 """
1332 mask.addMaskPlane(newMaskPlane)
1333 originBitMask = mask.getPlaneBitMask(maskPlane)
1334 destinationBitMask = mask.getPlaneBitMask(newMaskPlane)
1335 mask.array |= ((mask.array & originBitMask) > 0)*destinationBitMask
1337 def _clearMask(self, mask, clearMaskPlanes=None):
1338 """Clear the mask plane of an exposure.
1340 Parameters
1341 ----------
1342 mask : `lsst.afw.image.Mask`
1343 The mask plane to erase, which will be modified in place.
1344 clearMaskPlanes : `list` of `str`, optional
1345 Erase the specified mask planes.
1346 If not supplied, the entire mask will be erased.
1347 """
1348 if clearMaskPlanes is None: 1348 ↛ 1349line 1348 didn't jump to line 1349 because the condition on line 1348 was never true
1349 clearMaskPlanes = list(mask.getMaskPlaneDict().keys())
1351 bitMaskToClear = mask.getPlaneBitMask(clearMaskPlanes)
1352 mask &= ~bitMaskToClear
1355class AlardLuptonPreconvolveSubtractConnections(SubtractInputConnections,
1356 SubtractScoreOutputConnections):
1357 pass
1360class AlardLuptonPreconvolveSubtractConfig(AlardLuptonSubtractBaseConfig, lsst.pipe.base.PipelineTaskConfig,
1361 pipelineConnections=AlardLuptonPreconvolveSubtractConnections):
1362 pass
1365class AlardLuptonPreconvolveSubtractTask(AlardLuptonSubtractTask):
1366 """Subtract a template from a science image, convolving the science image
1367 before computing the kernel, and also convolving the template before
1368 subtraction.
1369 """
1370 ConfigClass = AlardLuptonPreconvolveSubtractConfig
1371 _DefaultName = "alardLuptonPreconvolveSubtract"
1372 usePreconvolution = True
1373 # This task just writes the score image instead of the difference image and
1374 # matchedTemplate
1375 futureImageOutputs = ("scoreExposure",)
1377 def run(self, template, science, sources, visitSummary=None):
1378 """Preconvolve the science image with its own PSF,
1379 convolve the template image with a PSF-matching kernel and subtract
1380 from the preconvolved science image.
1382 Parameters
1383 ----------
1384 template : `lsst.afw.image.ExposureF`
1385 The template image, which has previously been warped to the science
1386 image. The template bbox will be padded by a few pixels compared to
1387 the science bbox.
1388 science : `lsst.afw.image.ExposureF`
1389 The science exposure.
1390 sources : `lsst.afw.table.SourceCatalog`
1391 Identified sources on the science exposure. This catalog is used to
1392 select sources in order to perform the AL PSF matching on stamp
1393 images around them.
1394 visitSummary : `lsst.afw.table.ExposureCatalog`, optional
1395 Exposure catalog with complete external calibrations. Catalog uses
1396 the detector id for the catalog id, sorted on id for fast lookup.
1398 Returns
1399 -------
1400 results : `lsst.pipe.base.Struct`
1401 ``scoreExposure`` : `lsst.afw.image.ExposureF`
1402 Result of subtracting the convolved template and science
1403 images. Attached PSF is that of the original science image.
1404 ``matchedTemplate`` : `lsst.afw.image.ExposureF`
1405 Warped and PSF-matched template exposure. Attached PSF is that
1406 of the original science image.
1407 ``matchedScience`` : `lsst.afw.image.ExposureF`
1408 The science exposure after convolving with its own PSF.
1409 Attached PSF is that of the original science image.
1410 ``backgroundModel`` : `lsst.afw.math.Function2D`
1411 Background model that was fit while solving for the
1412 PSF-matching kernel
1413 ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
1414 Final kernel used to PSF-match the template to the science
1415 image.
1416 """
1417 self._prepareInputs(template, science, visitSummary=visitSummary)
1419 convolutionKernel = self._makePreconvolutionKernel(science.psf)
1420 matchedScience = self._convolveExposure(science, convolutionKernel, self.convolutionControl,
1421 interpolateBadMaskPlanes=True)
1422 self.metadata["convolvedExposure"] = "Preconvolution"
1424 self.matchedPsfSize = self.sciencePsfSize*np.sqrt(2)
1425 self.log.info("Preconvolved science PSF FWHM: %f pixels", self.matchedPsfSize)
1426 self.metadata["preconvolvedSciencePsfSize"] = self.matchedPsfSize
1427 try:
1428 kernelSources = self._sourceSelector(template, matchedScience, sources)
1429 subtractResults = self.runPreconvolve(template, science, matchedScience,
1430 kernelSources, convolutionKernel)
1432 except (RuntimeError, lsst.pex.exceptions.Exception) as e:
1433 self.log.warning("Failed to match template. Checking coverage")
1434 # Raise NoWorkFound if template fraction is insufficient
1435 checkTemplateIsSufficient(template[science.getBBox()], science, self.log,
1436 self.config.minTemplateFractionForExpectedSuccess,
1437 exceptionMessage="Template coverage lower than expected to succeed."
1438 f" Failure is tolerable: {e}")
1439 # checkTemplateIsSufficient did not raise NoWorkFound, so raise original exception
1440 raise e
1442 return subtractResults
1444 @staticmethod
1445 def _flagScoreEdge(mask, innerBBox):
1446 """Set the EDGE mask bit on pixels outside a known-valid region.
1448 Parameters
1449 ----------
1450 mask : `~lsst.afw.image.Mask`
1451 Exposure mask that will be modified in place. Must have
1452 an ``EDGE`` mask plane.
1453 innerBBox : `~lsst.geom.Box2I`
1454 The valid inner region. Pixels
1455 outside this bbox will have their ``EDGE`` bit set.
1456 """
1457 bbox = mask.getBBox()
1458 edgeBit = mask.getPlaneBitMask("EDGE")
1459 dx0 = innerBBox.getMinX() - bbox.getMinX()
1460 dx1 = bbox.getMaxX() - innerBBox.getMaxX()
1461 dy0 = innerBBox.getMinY() - bbox.getMinY()
1462 dy1 = bbox.getMaxY() - innerBBox.getMaxY()
1463 if dy0 > 0: 1463 ↛ 1465line 1463 didn't jump to line 1465 because the condition on line 1463 was always true
1464 mask.array[:dy0, :] |= edgeBit
1465 if dy1 > 0: 1465 ↛ 1467line 1465 didn't jump to line 1467 because the condition on line 1465 was always true
1466 mask.array[-dy1:, :] |= edgeBit
1467 if dx0 > 0: 1467 ↛ 1469line 1467 didn't jump to line 1469 because the condition on line 1467 was always true
1468 mask.array[:, :dx0] |= edgeBit
1469 if dx1 > 0: 1469 ↛ exitline 1469 didn't return from function '_flagScoreEdge' because the condition on line 1469 was always true
1470 mask.array[:, -dx1:] |= edgeBit
1472 @staticmethod
1473 def _makePreconvolutionKernel(psf):
1474 """Build a normalized, reflected matched-filter kernel from a PSF.
1476 Convolving an image with this kernel is equivalent to correlating
1477 the image with the PSF, so peaks in the output align with the PSF's
1478 centroid — even for asymmetric PSFs. The kernel is evaluated at the
1479 PSF's average position and returned as a constant
1480 `~lsst.afw.math.Kernel`.
1482 Parameters
1483 ----------
1484 psf : `~lsst.afw.detection.Psf`
1485 The PSF to derive the preconvolution kernel from.
1487 Returns
1488 -------
1489 kernel : `~lsst.afw.math.Kernel`
1490 The PSF reflected about both axes, normalized to sum to one.
1492 Raises
1493 ------
1494 ValueError
1495 Raised if the PSF kernel has an even size along either axis.
1496 It's not possible to center an even-sized kernel.
1497 """
1498 avgPos = psf.getAveragePosition()
1499 localKernel = psf.getLocalKernel(avgPos)
1500 dims = localKernel.getDimensions()
1501 if dims.x % 2 == 0 or dims.y % 2 == 0: 1501 ↛ 1502line 1501 didn't jump to line 1502 because the condition on line 1501 was never true
1502 raise ValueError(
1503 f"Preconvolution requires an odd-sized PSF kernel, got {dims.x}x{dims.y}. "
1504 )
1505 kimg = lsst.afw.image.ImageD(dims)
1506 localKernel.computeImage(kimg, doNormalize=True) # normalize to unit sum
1507 # Reflect about the kernel center. PSF kernels are odd-sized,
1508 # so ``[::-1, ::-1]`` places the peak at the same pixel.
1509 kimg.array[...] = kimg.array[::-1, ::-1]
1510 return lsst.afw.math.FixedKernel(kimg)
1512 def runPreconvolve(self, template, science, matchedScience, kernelSources, preConvKernel):
1513 """Convolve the science image with its own PSF, then convolve the
1514 template with a matching kernel and subtract to form the Score
1515 exposure.
1517 Parameters
1518 ----------
1519 template : `lsst.afw.image.ExposureF`
1520 Template exposure, warped to match the science exposure.
1521 science : `lsst.afw.image.ExposureF`
1522 Science exposure to subtract from the template.
1523 matchedScience : `lsst.afw.image.ExposureF`
1524 The science exposure, convolved with the reflection of its own PSF.
1525 kernelSources : `lsst.afw.table.SourceCatalog`
1526 Identified sources on the science exposure. This catalog is used to
1527 select sources in order to perform the AL PSF matching on stamp
1528 images around them.
1529 preConvKernel : `lsst.afw.math.Kernel`
1530 The kernel that was used to preconvolve the ``science``
1531 exposure. Must be normalized to sum to 1.
1533 Returns
1534 -------
1535 results : `lsst.pipe.base.Struct`
1537 ``scoreExposure`` : `lsst.afw.image.ExposureF`
1538 Result of subtracting the convolved template and science
1539 images. Attached PSF is that of the original science image.
1540 ``matchedTemplate`` : `lsst.afw.image.ExposureF`
1541 Warped and PSF-matched template exposure. Attached PSF is that
1542 of the original science image.
1543 ``matchedScience`` : `lsst.afw.image.ExposureF`
1544 The science exposure after convolving with its own PSF.
1545 Attached PSF is that of the original science image.
1546 ``backgroundModel`` : `lsst.afw.math.Function2D`
1547 Background model that was fit while solving for the
1548 PSF-matching kernel
1549 ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
1550 Final kernel used to PSF-match the template to the science
1551 image.
1552 """
1553 bbox = science.getBBox()
1554 innerBBox = preConvKernel.shrinkBBox(bbox)
1556 kernelResult = self.makeKernel.run(template[innerBBox], matchedScience[innerBBox], kernelSources,
1557 preconvolved=True,
1558 templateFwhmPix=self.templatePsfSize,
1559 scienceFwhmPix=self.matchedPsfSize)
1561 matchedTemplate = self._convolveExposure(template, kernelResult.psfMatchingKernel,
1562 self.convolutionControl,
1563 bbox=bbox,
1564 psf=science.psf,
1565 interpolateBadMaskPlanes=True,
1566 photoCalib=science.photoCalib)
1567 score = _subtractImages(matchedScience, matchedTemplate,
1568 backgroundModel=(kernelResult.backgroundModel
1569 if self.config.doSubtractBackground else None))
1570 correctedScore = self.finalize(template[bbox], science, score,
1571 kernelResult.psfMatchingKernel,
1572 templateMatched=True, preConvMode=True,
1573 preConvKernel=preConvKernel)
1575 # Flag the outer ``preConvKernel/2``-wide border as EDGE.
1576 self._flagScoreEdge(correctedScore.mask, innerBBox)
1578 return lsst.pipe.base.Struct(scoreExposure=correctedScore,
1579 matchedTemplate=matchedTemplate,
1580 matchedScience=matchedScience,
1581 backgroundModel=kernelResult.backgroundModel,
1582 psfMatchingKernel=kernelResult.psfMatchingKernel,
1583 kernelSources=kernelSources)
1586def checkTemplateIsSufficient(templateExposure, scienceExposure, logger, requiredTemplateFraction=0.,
1587 exceptionMessage=""):
1588 """Raise NoWorkFound if template coverage < requiredTemplateFraction
1590 Parameters
1591 ----------
1592 templateExposure : `lsst.afw.image.ExposureF`
1593 The template exposure to check
1594 logger : `logging.Logger`
1595 Logger for printing output.
1596 requiredTemplateFraction : `float`, optional
1597 Fraction of pixels of the science image required to have coverage
1598 in the template.
1599 exceptionMessage : `str`, optional
1600 Message to include in the exception raised if the template coverage
1601 is insufficient.
1603 Returns
1604 -------
1605 templateCoverageFraction: `float`
1606 Fraction of pixels in the template with data.
1608 Raises
1609 ------
1610 lsst.pipe.base.NoWorkFound
1611 Raised if fraction of good pixels, defined as not having NO_DATA
1612 set, is less than the requiredTemplateFraction
1613 """
1614 # Count the number of pixels with the NO_DATA mask bit set
1615 # counting NaN pixels is insufficient because pixels without data are often intepolated over)
1616 noTemplate = templateExposure.mask.array & templateExposure.mask.getPlaneBitMask('NO_DATA')
1617 # Also need to account for missing data in the science image,
1618 # because template coverage there doesn't help
1619 noScience = scienceExposure.mask.array & scienceExposure.mask.getPlaneBitMask('NO_DATA')
1620 pixNoData = np.count_nonzero(noTemplate | noScience)
1621 pixGood = templateExposure.getBBox().getArea() - pixNoData
1622 templateCoverageFraction = pixGood/templateExposure.getBBox().getArea()
1623 logger.info("template has %d good pixels (%.1f%%)", pixGood, 100*templateCoverageFraction)
1625 if templateCoverageFraction < requiredTemplateFraction:
1626 message = ("Insufficient Template Coverage. (%.1f%% < %.1f%%)" % (
1627 100*templateCoverageFraction,
1628 100*requiredTemplateFraction))
1629 raise lsst.pipe.base.NoWorkFound(message + " " + exceptionMessage)
1630 return templateCoverageFraction
1633def _subtractImages(science, template, backgroundModel=None):
1634 """Subtract template from science, propagating relevant metadata.
1636 Parameters
1637 ----------
1638 science : `lsst.afw.Exposure`
1639 The input science image.
1640 template : `lsst.afw.Exposure`
1641 The template to subtract from the science image.
1642 backgroundModel : `lsst.afw.MaskedImage`, optional
1643 Differential background model
1645 Returns
1646 -------
1647 difference : `lsst.afw.Exposure`
1648 The subtracted image.
1649 """
1650 difference = science.clone()
1651 if backgroundModel is not None:
1652 difference.maskedImage -= backgroundModel
1653 difference.maskedImage -= template.maskedImage
1654 return difference
1657def _shapeTest(exp1, exp2, fwhmExposureBuffer, fwhmExposureGrid):
1658 """Determine that the PSF of ``exp1`` is not wider than that of ``exp2``.
1660 Parameters
1661 ----------
1662 exp1 : `~lsst.afw.image.Exposure`
1663 Exposure with the reference point spread function (PSF) to evaluate.
1664 exp2 : `~lsst.afw.image.Exposure`
1665 Exposure with a candidate point spread function (PSF) to evaluate.
1666 fwhmExposureBuffer : `float`
1667 Fractional buffer margin to be left out of all sides of the image
1668 during the construction of the grid to compute mean PSF FWHM in an
1669 exposure, if the PSF is not available at its average position.
1670 fwhmExposureGrid : `int`
1671 Grid size to compute the mean FWHM in an exposure, if the PSF is not
1672 available at its average position.
1673 Returns
1674 -------
1675 result : `bool`
1676 True if ``exp1`` has a PSF that is not wider than that of ``exp2`` in
1677 either dimension.
1678 """
1679 try:
1680 shape1 = getPsfFwhm(exp1.psf, average=False)
1681 shape2 = getPsfFwhm(exp2.psf, average=False)
1682 except (lsst.pex.exceptions.InvalidParameterError, lsst.pex.exceptions.RangeError):
1683 shape1 = evaluateMeanPsfFwhm(exp1,
1684 fwhmExposureBuffer=fwhmExposureBuffer,
1685 fwhmExposureGrid=fwhmExposureGrid
1686 )
1687 shape2 = evaluateMeanPsfFwhm(exp2,
1688 fwhmExposureBuffer=fwhmExposureBuffer,
1689 fwhmExposureGrid=fwhmExposureGrid
1690 )
1691 return shape1 <= shape2
1693 # Results from getPsfFwhm is a tuple of two values, one for each dimension.
1694 xTest = shape1[0] <= shape2[0]
1695 yTest = shape1[1] <= shape2[1]
1696 return xTest | yTest
1699class SimplifiedSubtractConfig(AlardLuptonSubtractBaseConfig, lsst.pipe.base.PipelineTaskConfig,
1700 pipelineConnections=SimplifiedSubtractConnections):
1701 mode = lsst.pex.config.ChoiceField(
1702 dtype=str,
1703 default="convolveTemplate",
1704 allowed={"auto": "Choose which image to convolve at runtime.",
1705 "convolveScience": "Only convolve the science image.",
1706 "convolveTemplate": "Only convolve the template image."},
1707 doc="Choose which image to convolve at runtime, or require that a specific image is convolved."
1708 )
1709 useExistingKernel = lsst.pex.config.Field(
1710 dtype=bool,
1711 default=True,
1712 doc="Use a pre-existing PSF matching kernel?"
1713 "If False, source detection and measurement will be run."
1714 )
1717class SimplifiedSubtractTask(AlardLuptonSubtractTask):
1718 """Compute the image difference of a science and template image using
1719 the Alard & Lupton (1998) algorithm.
1720 """
1721 ConfigClass = SimplifiedSubtractConfig
1722 _DefaultName = "simplifiedSubtract"
1724 @timeMethod
1725 def run(self, template, science, visitSummary=None, inputPsfMatchingKernel=None):
1726 """PSF match, subtract, and decorrelate two images.
1728 Parameters
1729 ----------
1730 template : `lsst.afw.image.ExposureF`
1731 Template exposure, warped to match the science exposure.
1732 science : `lsst.afw.image.ExposureF`
1733 Science exposure to subtract from the template.
1734 visitSummary : `lsst.afw.table.ExposureCatalog`, optional
1735 Exposure catalog with external calibrations to be applied. Catalog
1736 uses the detector id for the catalog id, sorted on id for fast
1737 lookup.
1738 inputPsfMatchingKernel : `lsst.afw.math.Kernel`, optional
1739 Pre-existing PSF matching kernel to use for convolution.
1740 Required, and only used, if ``config.useExistingKernel`` is set.
1742 Returns
1743 -------
1744 results : `lsst.pipe.base.Struct`
1745 ``difference`` : `lsst.afw.image.ExposureF`
1746 Result of subtracting template and science.
1747 ``matchedTemplate`` : `lsst.afw.image.ExposureF`
1748 Warped and PSF-matched template exposure.
1749 ``backgroundModel`` : `lsst.afw.math.Function2D`
1750 Background model that was fit while solving for the
1751 PSF-matching kernel
1752 ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
1753 Kernel used to PSF-match the convolved image.
1754 ``kernelSources` : `lsst.afw.table.SourceCatalog`
1755 Sources detected on the science image that were used to
1756 construct the PSF-matching kernel.
1758 Raises
1759 ------
1760 lsst.pipe.base.NoWorkFound
1761 Raised if fraction of good pixels, defined as not having NO_DATA
1762 set, is less then the configured requiredTemplateFraction
1763 """
1764 self._prepareInputs(template, science, visitSummary=visitSummary)
1766 convolveTemplate = self.chooseConvolutionMethod(template, science)
1767 self.matchedPsfSize = self.sciencePsfSize if convolveTemplate else self.templatePsfSize
1769 if self.config.useExistingKernel:
1770 psfMatchingKernel = inputPsfMatchingKernel
1771 backgroundModel = None
1772 kernelSources = None
1773 else:
1774 kernelResult = self.runMakeKernel(template, science, convolveTemplate=convolveTemplate,
1775 runSourceDetection=True)
1776 psfMatchingKernel = kernelResult.psfMatchingKernel
1777 kernelSources = kernelResult.kernelSources
1778 if self.config.doSubtractBackground: 1778 ↛ 1779line 1778 didn't jump to line 1779 because the condition on line 1778 was never true
1779 backgroundModel = kernelResult.backgroundModel
1780 else:
1781 backgroundModel = None
1782 if convolveTemplate: 1782 ↛ 1786line 1782 didn't jump to line 1786 because the condition on line 1782 was always true
1783 subtractResults = self.runConvolveTemplate(template, science, psfMatchingKernel,
1784 backgroundModel=backgroundModel)
1785 else:
1786 subtractResults = self.runConvolveScience(template, science, psfMatchingKernel,
1787 backgroundModel=backgroundModel)
1788 if kernelSources is not None:
1789 subtractResults.kernelSources = kernelSources
1790 return subtractResults
1793def _interpolateImage(maskedImage, badMaskPlanes, fallbackValue=None):
1794 """Replace masked image pixels with interpolated values.
1796 Parameters
1797 ----------
1798 maskedImage : `lsst.afw.image.MaskedImage`
1799 Image on which to perform interpolation.
1800 badMaskPlanes : `list` of `str`
1801 List of mask planes to interpolate over.
1802 fallbackValue : `float`, optional
1803 Value to set when interpolation fails.
1805 Returns
1806 -------
1807 result: `float`
1808 The number of masked pixels that were replaced.
1809 """
1810 imgBadMaskPlanes = [
1811 maskPlane for maskPlane in badMaskPlanes if maskPlane in maskedImage.mask.getMaskPlaneDict()
1812 ]
1814 image = maskedImage.image.array
1815 badPixels = (maskedImage.mask.array & maskedImage.mask.getPlaneBitMask(imgBadMaskPlanes)) > 0
1816 image[badPixels] = np.nan
1817 if fallbackValue is None: 1817 ↛ 1821line 1817 didn't jump to line 1821 because the condition on line 1817 was always true
1818 fallbackValue = np.nanmedian(image)
1819 # For this initial implementation, skip the interpolation and just fill with
1820 # the median value.
1821 image[badPixels] = fallbackValue
1822 return np.sum(badPixels)