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