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