Loading [MathJax]/extensions/tex2jax.js
LSST Applications 28.0.0,g1653933729+a8ce1bb630,g1a997c3884+a8ce1bb630,g28da252d5a+5bd70b7e6d,g2bbee38e9b+638fca75ac,g2bc492864f+638fca75ac,g3156d2b45e+07302053f8,g347aa1857d+638fca75ac,g35bb328faa+a8ce1bb630,g3a166c0a6a+638fca75ac,g3e281a1b8c+7bbb0b2507,g4005a62e65+17cd334064,g414038480c+5b5cd4fff3,g41af890bb2+4ffae9de63,g4e1a3235cc+0f1912dca3,g6249c6f860+3c3976f90c,g80478fca09+46aba80bd6,g82479be7b0+77990446f6,g858d7b2824+78ba4d1ce1,g89c8672015+f667a5183b,g9125e01d80+a8ce1bb630,ga5288a1d22+2a6264e9ca,gae0086650b+a8ce1bb630,gb58c049af0+d64f4d3760,gc22bb204ba+78ba4d1ce1,gc28159a63d+638fca75ac,gcf0d15dbbd+32ddb6096f,gd6b7c0dfd1+3e339405e9,gda3e153d99+78ba4d1ce1,gda6a2b7d83+32ddb6096f,gdaeeff99f8+1711a396fd,gdd5a9049c5+b18c39e5e3,ge2409df99d+a5e4577cdc,ge33fd446bb+78ba4d1ce1,ge79ae78c31+638fca75ac,gf0baf85859+64e8883e75,gf5289d68f6+e1b046a8d7,gfa443fc69c+91d9ed1ecf,gfda6b12a05+8419469a56
LSST Data Management Base Package
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
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
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
41_dimensions = ("instrument", "visit", "detector")
42_defaultTemplates = {"coaddName": "deep", "fakesType": ""}
43
44
45class SubtractInputConnections(lsst.pipe.base.PipelineTaskConnections,
46 dimensions=_dimensions,
47 defaultTemplates=_defaultTemplates):
48 template = connectionTypes.Input(
49 doc="Input warped template to subtract.",
50 dimensions=("instrument", "visit", "detector"),
51 storageClass="ExposureF",
52 name="{fakesType}{coaddName}Diff_templateExp"
53 )
54 science = connectionTypes.Input(
55 doc="Input science exposure to subtract from.",
56 dimensions=("instrument", "visit", "detector"),
57 storageClass="ExposureF",
58 name="{fakesType}calexp"
59 )
60 sources = connectionTypes.Input(
61 doc="Sources measured on the science exposure; "
62 "used to select sources for making the matching kernel.",
63 dimensions=("instrument", "visit", "detector"),
64 storageClass="SourceCatalog",
65 name="{fakesType}src"
66 )
67 visitSummary = connectionTypes.Input(
68 doc=("Per-visit catalog with final calibration objects. "
69 "These catalogs use the detector id for the catalog id, "
70 "sorted on id for fast lookup."),
71 dimensions=("instrument", "visit"),
72 storageClass="ExposureCatalog",
73 name="finalVisitSummary",
74 )
75
76 def __init__(self, *, config=None):
77 super().__init__(config=config)
78 if not config.doApplyExternalCalibrations:
79 del self.visitSummary
80
81
82class SubtractImageOutputConnections(lsst.pipe.base.PipelineTaskConnections,
83 dimensions=_dimensions,
84 defaultTemplates=_defaultTemplates):
85 difference = connectionTypes.Output(
86 doc="Result of subtracting convolved template from science image.",
87 dimensions=("instrument", "visit", "detector"),
88 storageClass="ExposureF",
89 name="{fakesType}{coaddName}Diff_differenceTempExp",
90 )
91 matchedTemplate = connectionTypes.Output(
92 doc="Warped and PSF-matched template used to create `subtractedExposure`.",
93 dimensions=("instrument", "visit", "detector"),
94 storageClass="ExposureF",
95 name="{fakesType}{coaddName}Diff_matchedExp",
96 )
97 psfMatchingKernel = connectionTypes.Output(
98 doc="Kernel used to PSF match the science and template images.",
99 dimensions=("instrument", "visit", "detector"),
100 storageClass="MatchingKernel",
101 name="{fakesType}{coaddName}Diff_psfMatchKernel",
102 )
103
104
105class SubtractScoreOutputConnections(lsst.pipe.base.PipelineTaskConnections,
106 dimensions=_dimensions,
107 defaultTemplates=_defaultTemplates):
108 scoreExposure = connectionTypes.Output(
109 doc="The maximum likelihood image, used for the detection of diaSources.",
110 dimensions=("instrument", "visit", "detector"),
111 storageClass="ExposureF",
112 name="{fakesType}{coaddName}Diff_scoreExp",
113 )
114 psfMatchingKernel = connectionTypes.Output(
115 doc="Kernel used to PSF match the science and template images.",
116 dimensions=("instrument", "visit", "detector"),
117 storageClass="MatchingKernel",
118 name="{fakesType}{coaddName}Diff_psfScoreMatchKernel",
119 )
120
121
128 target=MakeKernelTask,
129 doc="Task to construct a matching kernel for convolution.",
130 )
131 doDecorrelation = lsst.pex.config.Field(
132 dtype=bool,
133 default=True,
134 doc="Perform diffim decorrelation to undo pixel correlation due to A&L "
135 "kernel convolution? If True, also update the diffim PSF."
136 )
138 target=DecorrelateALKernelTask,
139 doc="Task to decorrelate the image difference.",
140 )
141 requiredTemplateFraction = lsst.pex.config.Field(
142 dtype=float,
143 default=0.1,
144 doc="Raise NoWorkFound and do not attempt image subtraction if template covers less than this "
145 " fraction of pixels. Setting to 0 will always attempt image subtraction."
146 )
147 minTemplateFractionForExpectedSuccess = lsst.pex.config.Field(
148 dtype=float,
149 default=0.2,
150 doc="Raise NoWorkFound if PSF-matching fails and template covers less than this fraction of pixels."
151 " If the fraction of pixels covered by the template is less than this value (and greater than"
152 " requiredTemplateFraction) this task is attempted but failure is anticipated and tolerated."
153 )
154 doScaleVariance = lsst.pex.config.Field(
155 dtype=bool,
156 default=True,
157 doc="Scale variance of the image difference?"
158 )
160 target=ScaleVarianceTask,
161 doc="Subtask to rescale the variance of the template to the statistically expected level."
162 )
163 doSubtractBackground = lsst.pex.config.Field(
164 doc="Subtract the background fit when solving the kernel?",
165 dtype=bool,
166 default=True,
167 )
168 doApplyExternalCalibrations = lsst.pex.config.Field(
169 doc=(
170 "Replace science Exposure's calibration objects with those"
171 " in visitSummary. Ignored if `doApplyFinalizedPsf is True."
172 ),
173 dtype=bool,
174 default=False,
175 )
177 target=ScienceSourceSelectorTask,
178 doc="Task to select sources to be used for PSF matching.",
179 )
180 detectionThreshold = lsst.pex.config.Field(
181 dtype=float,
182 default=10,
183 doc="Minimum signal to noise ratio of detected sources "
184 "to use for calculating the PSF matching kernel."
185 )
186 detectionThresholdMax = lsst.pex.config.Field(
187 dtype=float,
188 default=500,
189 doc="Maximum signal to noise ratio of detected sources "
190 "to use for calculating the PSF matching kernel."
191 )
192 maxKernelSources = lsst.pex.config.Field(
193 dtype=int,
194 default=1000,
195 doc="Maximum number of sources to use for calculating the PSF matching kernel."
196 "Set to -1 to disable."
197 )
198 minKernelSources = lsst.pex.config.Field(
199 dtype=int,
200 default=3,
201 doc="Minimum number of sources needed for calculating the PSF matching kernel."
202 )
204 dtype=str,
205 doc="Flags that, if set, the associated source should not "
206 "be used to determine the PSF matching kernel.",
207 default=("sky_source", "slot_Centroid_flag",
208 "slot_ApFlux_flag", "slot_PsfFlux_flag",
209 "base_PixelFlags_flag_interpolated",
210 "base_PixelFlags_flag_saturated",
211 "base_PixelFlags_flag_bad",
212 ),
213 deprecated="This field is no longer used and will be removed after v28."
214 "Set the equivalent field for the sourceSelector subtask instead.",
215 )
216 excludeMaskPlanes = lsst.pex.config.ListField(
217 dtype=str,
218 default=("NO_DATA", "BAD", "SAT", "EDGE", "FAKE"),
219 doc="Mask planes to exclude when selecting sources for PSF matching."
220 )
222 dtype=str,
223 default=("NO_DATA", "BAD", "SAT", "EDGE"),
224 doc="Mask planes to interpolate over."
225 )
226 preserveTemplateMask = lsst.pex.config.ListField(
227 dtype=str,
228 default=("NO_DATA", "BAD",),
229 doc="Mask planes from the template to propagate to the image difference."
230 )
231 renameTemplateMask = lsst.pex.config.ListField(
232 dtype=str,
233 default=("SAT", "INJECTED", "INJECTED_CORE",),
234 doc="Mask planes from the template to propagate to the image difference"
235 "with '_TEMPLATE' appended to the name."
236 )
237 allowKernelSourceDetection = lsst.pex.config.Field(
238 dtype=bool,
239 default=False,
240 doc="Re-run source detection for kernel candidates if an error is"
241 " encountered while calculating the matching kernel."
242 )
243
244 def setDefaults(self):
245 self.makeKernel.kernel.name = "AL"
246 self.makeKernel.kernel.active.fitForBackground = self.doSubtractBackground
247 self.makeKernel.kernel.active.spatialKernelOrder = 1
248 self.makeKernel.kernel.active.spatialBgOrder = 2
249 self.sourceSelector.doUnresolved = True # apply star-galaxy separation
250 self.sourceSelector.doIsolated = True # apply isolated star selection
251 self.sourceSelector.doRequirePrimary = True # apply primary flag selection
252 self.sourceSelector.doSkySources = False # Do not include sky sources
253 self.sourceSelector.doSignalToNoise = True # apply signal to noise filter
254 self.sourceSelector.signalToNoise.minimum = 10
255 self.sourceSelector.signalToNoise.maximum = 500
256
257
258class AlardLuptonSubtractConfig(AlardLuptonSubtractBaseConfig, lsst.pipe.base.PipelineTaskConfig,
259 pipelineConnections=AlardLuptonSubtractConnections):
261 dtype=str,
262 default="convolveTemplate",
263 allowed={"auto": "Choose which image to convolve at runtime.",
264 "convolveScience": "Only convolve the science image.",
265 "convolveTemplate": "Only convolve the template image."},
266 doc="Choose which image to convolve at runtime, or require that a specific image is convolved."
267 )
268
269
270class AlardLuptonSubtractTask(lsst.pipe.base.PipelineTask):
271 """Compute the image difference of a science and template image using
272 the Alard & Lupton (1998) algorithm.
273 """
274 ConfigClass = AlardLuptonSubtractConfig
275 _DefaultName = "alardLuptonSubtract"
276
277 def __init__(self, **kwargs):
278 super().__init__(**kwargs)
279 self.makeSubtask("decorrelate")
280 self.makeSubtask("makeKernel")
281 self.makeSubtask("sourceSelector")
282 if self.config.doScaleVariance:
283 self.makeSubtask("scaleVariance")
284
286 # Normalization is an extra, unnecessary, calculation and will result
287 # in mis-subtraction of the images if there are calibration errors.
288 self.convolutionControl.setDoNormalize(False)
289 self.convolutionControl.setDoCopyEdge(True)
290
291 def _applyExternalCalibrations(self, exposure, visitSummary):
292 """Replace calibrations (psf, and ApCorrMap) on this exposure with
293 external ones.".
294
295 Parameters
296 ----------
297 exposure : `lsst.afw.image.exposure.Exposure`
298 Input exposure to adjust calibrations.
299 visitSummary : `lsst.afw.table.ExposureCatalog`
300 Exposure catalog with external calibrations to be applied. Catalog
301 uses the detector id for the catalog id, sorted on id for fast
302 lookup.
303
304 Returns
305 -------
306 exposure : `lsst.afw.image.exposure.Exposure`
307 Exposure with adjusted calibrations.
308 """
309 detectorId = exposure.info.getDetector().getId()
310
311 row = visitSummary.find(detectorId)
312 if row is None:
313 self.log.warning("Detector id %s not found in external calibrations catalog; "
314 "Using original calibrations.", detectorId)
315 else:
316 psf = row.getPsf()
317 apCorrMap = row.getApCorrMap()
318 if psf is None:
319 self.log.warning("Detector id %s has None for psf in "
320 "external calibrations catalog; Using original psf and aperture correction.",
321 detectorId)
322 elif apCorrMap is None:
323 self.log.warning("Detector id %s has None for apCorrMap in "
324 "external calibrations catalog; Using original psf and aperture correction.",
325 detectorId)
326 else:
327 exposure.setPsf(psf)
328 exposure.info.setApCorrMap(apCorrMap)
329
330 return exposure
331
332 @timeMethod
333 def run(self, template, science, sources, visitSummary=None):
334 """PSF match, subtract, and decorrelate two images.
335
336 Parameters
337 ----------
338 template : `lsst.afw.image.ExposureF`
339 Template exposure, warped to match the science exposure.
340 science : `lsst.afw.image.ExposureF`
341 Science exposure to subtract from the template.
342 sources : `lsst.afw.table.SourceCatalog`
343 Identified sources on the science exposure. This catalog is used to
344 select sources in order to perform the AL PSF matching on stamp
345 images around them.
346 visitSummary : `lsst.afw.table.ExposureCatalog`, optional
347 Exposure catalog with external calibrations to be applied. Catalog
348 uses the detector id for the catalog id, sorted on id for fast
349 lookup.
350
351 Returns
352 -------
353 results : `lsst.pipe.base.Struct`
354 ``difference`` : `lsst.afw.image.ExposureF`
355 Result of subtracting template and science.
356 ``matchedTemplate`` : `lsst.afw.image.ExposureF`
357 Warped and PSF-matched template exposure.
358 ``backgroundModel`` : `lsst.afw.math.Function2D`
359 Background model that was fit while solving for the
360 PSF-matching kernel
361 ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
362 Kernel used to PSF-match the convolved image.
363
364 Raises
365 ------
366 RuntimeError
367 If an unsupported convolution mode is supplied.
368 RuntimeError
369 If there are too few sources to calculate the PSF matching kernel.
370 lsst.pipe.base.NoWorkFound
371 Raised if fraction of good pixels, defined as not having NO_DATA
372 set, is less then the configured requiredTemplateFraction
373 """
374
375 self._prepareInputs(template, science, visitSummary=visitSummary)
376
377 # In the event that getPsfFwhm fails, evaluate the PSF on a grid.
378 fwhmExposureBuffer = self.config.makeKernel.fwhmExposureBuffer
379 fwhmExposureGrid = self.config.makeKernel.fwhmExposureGrid
380
381 # Calling getPsfFwhm on template.psf fails on some rare occasions when
382 # the template has no input exposures at the average position of the
383 # stars. So we try getPsfFwhm first on template, and if that fails we
384 # evaluate the PSF on a grid specified by fwhmExposure* fields.
385 # To keep consistent definitions for PSF size on the template and
386 # science images, we use the same method for both.
387 # In the try block below, we catch two exceptions:
388 # 1. InvalidParameterError, in case the point where we are evaluating
389 # the PSF lands in a gap in the template.
390 # 2. RangeError, in case the template coverage is so poor that we end
391 # up near a region with no data.
392 try:
393 templatePsfSize = getPsfFwhm(template.psf)
394 sciencePsfSize = getPsfFwhm(science.psf)
396 self.log.info("Unable to evaluate PSF at the average position. "
397 "Evaluting PSF on a grid of points."
398 )
399 templatePsfSize = evaluateMeanPsfFwhm(template,
400 fwhmExposureBuffer=fwhmExposureBuffer,
401 fwhmExposureGrid=fwhmExposureGrid
402 )
403 sciencePsfSize = evaluateMeanPsfFwhm(science,
404 fwhmExposureBuffer=fwhmExposureBuffer,
405 fwhmExposureGrid=fwhmExposureGrid
406 )
407 self.log.info("Science PSF FWHM: %f pixels", sciencePsfSize)
408 self.log.info("Template PSF FWHM: %f pixels", templatePsfSize)
409 self.metadata.add("sciencePsfSize", sciencePsfSize)
410 self.metadata.add("templatePsfSize", templatePsfSize)
411
412 # Calculate estimated image depths, i.e., limiting magnitudes
413 maglim_science = self._calculateMagLim(science, fallbackPsfSize=sciencePsfSize)
414 fluxlim_science = (maglim_science*u.ABmag).to_value(u.nJy)
415 maglim_template = self._calculateMagLim(template, fallbackPsfSize=templatePsfSize)
416 if np.isnan(maglim_template):
417 self.log.info("Cannot evaluate template limiting mag; adopting science limiting mag for diffim")
418 maglim_diffim = maglim_science
419 else:
420 fluxlim_template = (maglim_template*u.ABmag).to_value(u.nJy)
421 maglim_diffim = (np.sqrt(fluxlim_science**2 + fluxlim_template**2)*u.nJy).to(u.ABmag).value
422 self.metadata.add("scienceLimitingMagnitude", maglim_science)
423 self.metadata.add("templateLimitingMagnitude", maglim_template)
424 self.metadata.add("diffimLimitingMagnitude", maglim_diffim)
425
426 if self.config.mode == "auto":
427 convolveTemplate = _shapeTest(template,
428 science,
429 fwhmExposureBuffer=fwhmExposureBuffer,
430 fwhmExposureGrid=fwhmExposureGrid)
431 if convolveTemplate:
432 if sciencePsfSize < templatePsfSize:
433 self.log.info("Average template PSF size is greater, "
434 "but science PSF greater in one dimension: convolving template image.")
435 else:
436 self.log.info("Science PSF size is greater: convolving template image.")
437 else:
438 self.log.info("Template PSF size is greater: convolving science image.")
439 elif self.config.mode == "convolveTemplate":
440 self.log.info("`convolveTemplate` is set: convolving template image.")
441 convolveTemplate = True
442 elif self.config.mode == "convolveScience":
443 self.log.info("`convolveScience` is set: convolving science image.")
444 convolveTemplate = False
445 else:
446 raise RuntimeError("Cannot handle AlardLuptonSubtract mode: %s", self.config.mode)
447
448 try:
449 sourceMask = science.mask.clone()
450 sourceMask.array |= template[science.getBBox()].mask.array
451 selectSources = self._sourceSelector(sources, sourceMask)
452 if convolveTemplate:
453 self.metadata.add("convolvedExposure", "Template")
454 subtractResults = self.runConvolveTemplate(template, science, selectSources)
455 else:
456 self.metadata.add("convolvedExposure", "Science")
457 subtractResults = self.runConvolveScience(template, science, selectSources)
458
459 except (RuntimeError, lsst.pex.exceptions.Exception) as e:
460 self.log.warning("Failed to match template. Checking coverage")
461 # Raise NoWorkFound if template fraction is insufficient
462 checkTemplateIsSufficient(template[science.getBBox()], self.log,
463 self.config.minTemplateFractionForExpectedSuccess,
464 exceptionMessage="Template coverage lower than expected to succeed."
465 f" Failure is tolerable: {e}")
466 # checkTemplateIsSufficient did not raise NoWorkFound, so raise original exception
467 raise e
468
469 return subtractResults
470
471 def runConvolveTemplate(self, template, science, selectSources):
472 """Convolve the template image with a PSF-matching kernel and subtract
473 from the science image.
474
475 Parameters
476 ----------
477 template : `lsst.afw.image.ExposureF`
478 Template exposure, warped to match the science exposure.
479 science : `lsst.afw.image.ExposureF`
480 Science exposure to subtract from the template.
481 selectSources : `lsst.afw.table.SourceCatalog`
482 Identified sources on the science exposure. This catalog is used to
483 select sources in order to perform the AL PSF matching on stamp
484 images around them.
485
486 Returns
487 -------
488 results : `lsst.pipe.base.Struct`
489
490 ``difference`` : `lsst.afw.image.ExposureF`
491 Result of subtracting template and science.
492 ``matchedTemplate`` : `lsst.afw.image.ExposureF`
493 Warped and PSF-matched template exposure.
494 ``backgroundModel`` : `lsst.afw.math.Function2D`
495 Background model that was fit while solving for the PSF-matching kernel
496 ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
497 Kernel used to PSF-match the template to the science image.
498 """
499 try:
500 kernelSources = self.makeKernel.selectKernelSources(template, science,
501 candidateList=selectSources,
502 preconvolved=False)
503 kernelResult = self.makeKernel.run(template, science, kernelSources,
504 preconvolved=False)
505 except Exception as e:
506 if self.config.allowKernelSourceDetection:
507 self.log.warning("Error encountered trying to construct the matching kernel"
508 f" Running source detection and retrying. {e}")
509 kernelSources = self.makeKernel.selectKernelSources(template, science,
510 candidateList=None,
511 preconvolved=False)
512 kernelResult = self.makeKernel.run(template, science, kernelSources,
513 preconvolved=False)
514 else:
515 raise e
516
517 matchedTemplate = self._convolveExposure(template, kernelResult.psfMatchingKernel,
519 bbox=science.getBBox(),
520 psf=science.psf,
521 photoCalib=science.photoCalib)
522
523 difference = _subtractImages(science, matchedTemplate,
524 backgroundModel=(kernelResult.backgroundModel
525 if self.config.doSubtractBackground else None))
526 correctedExposure = self.finalize(template, science, difference,
527 kernelResult.psfMatchingKernel,
528 templateMatched=True)
529
530 return lsst.pipe.base.Struct(difference=correctedExposure,
531 matchedTemplate=matchedTemplate,
532 matchedScience=science,
533 backgroundModel=kernelResult.backgroundModel,
534 psfMatchingKernel=kernelResult.psfMatchingKernel)
535
536 def runConvolveScience(self, template, science, selectSources):
537 """Convolve the science image with a PSF-matching kernel and subtract
538 the template image.
539
540 Parameters
541 ----------
542 template : `lsst.afw.image.ExposureF`
543 Template exposure, warped to match the science exposure.
544 science : `lsst.afw.image.ExposureF`
545 Science exposure to subtract from the template.
546 selectSources : `lsst.afw.table.SourceCatalog`
547 Identified sources on the science exposure. This catalog is used to
548 select sources in order to perform the AL PSF matching on stamp
549 images around them.
550
551 Returns
552 -------
553 results : `lsst.pipe.base.Struct`
554
555 ``difference`` : `lsst.afw.image.ExposureF`
556 Result of subtracting template and science.
557 ``matchedTemplate`` : `lsst.afw.image.ExposureF`
558 Warped template exposure. Note that in this case, the template
559 is not PSF-matched to the science image.
560 ``backgroundModel`` : `lsst.afw.math.Function2D`
561 Background model that was fit while solving for the PSF-matching kernel
562 ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
563 Kernel used to PSF-match the science image to the template.
564 """
565 bbox = science.getBBox()
566 kernelSources = self.makeKernel.selectKernelSources(science, template,
567 candidateList=selectSources,
568 preconvolved=False)
569 kernelResult = self.makeKernel.run(science, template, kernelSources,
570 preconvolved=False)
571 modelParams = kernelResult.backgroundModel.getParameters()
572 # We must invert the background model if the matching kernel is solved for the science image.
573 kernelResult.backgroundModel.setParameters([-p for p in modelParams])
574
575 kernelImage = lsst.afw.image.ImageD(kernelResult.psfMatchingKernel.getDimensions())
576 norm = kernelResult.psfMatchingKernel.computeImage(kernelImage, doNormalize=False)
577
578 matchedScience = self._convolveExposure(science, kernelResult.psfMatchingKernel,
580 psf=template.psf)
581
582 # Place back on native photometric scale
583 matchedScience.maskedImage /= norm
584 matchedTemplate = template.clone()[bbox]
585 matchedTemplate.maskedImage /= norm
586 matchedTemplate.setPhotoCalib(science.photoCalib)
587
588 difference = _subtractImages(matchedScience, matchedTemplate,
589 backgroundModel=(kernelResult.backgroundModel
590 if self.config.doSubtractBackground else None))
591
592 correctedExposure = self.finalize(template, science, difference,
593 kernelResult.psfMatchingKernel,
594 templateMatched=False)
595
596 return lsst.pipe.base.Struct(difference=correctedExposure,
597 matchedTemplate=matchedTemplate,
598 matchedScience=matchedScience,
599 backgroundModel=kernelResult.backgroundModel,
600 psfMatchingKernel=kernelResult.psfMatchingKernel,)
601
602 def finalize(self, template, science, difference, kernel,
603 templateMatched=True,
604 preConvMode=False,
605 preConvKernel=None,
606 spatiallyVarying=False):
607 """Decorrelate the difference image to undo the noise correlations
608 caused by convolution.
609
610 Parameters
611 ----------
612 template : `lsst.afw.image.ExposureF`
613 Template exposure, warped to match the science exposure.
614 science : `lsst.afw.image.ExposureF`
615 Science exposure to subtract from the template.
616 difference : `lsst.afw.image.ExposureF`
617 Result of subtracting template and science.
618 kernel : `lsst.afw.math.Kernel`
619 An (optionally spatially-varying) PSF matching kernel
620 templateMatched : `bool`, optional
621 Was the template PSF-matched to the science image?
622 preConvMode : `bool`, optional
623 Was the science image preconvolved with its own PSF
624 before PSF matching the template?
625 preConvKernel : `lsst.afw.detection.Psf`, optional
626 If not `None`, then the science image was pre-convolved with
627 (the reflection of) this kernel. Must be normalized to sum to 1.
628 spatiallyVarying : `bool`, optional
629 Compute the decorrelation kernel spatially varying across the image?
630
631 Returns
632 -------
633 correctedExposure : `lsst.afw.image.ExposureF`
634 The decorrelated image difference.
635 """
636 if self.config.doDecorrelation:
637 self.log.info("Decorrelating image difference.")
638 # We have cleared the template mask plane, so copy the mask plane of
639 # the image difference so that we can calculate correct statistics
640 # during decorrelation
641 correctedExposure = self.decorrelate.run(science, template[science.getBBox()], difference, kernel,
642 templateMatched=templateMatched,
643 preConvMode=preConvMode,
644 preConvKernel=preConvKernel,
645 spatiallyVarying=spatiallyVarying).correctedExposure
646 else:
647 self.log.info("NOT decorrelating image difference.")
648 correctedExposure = difference
649 return correctedExposure
650
651 def _calculateMagLim(self, exposure, nsigma=5.0, fallbackPsfSize=None):
652 """Calculate an exposure's limiting magnitude.
653
654 This method uses the photometric zeropoint together with the
655 PSF size from the average position of the exposure.
656
657 Parameters
658 ----------
659 exposure : `lsst.afw.image.Exposure`
660 The target exposure to calculate the limiting magnitude for.
661 nsigma : `float`, optional
662 The detection threshold in sigma.
663 fallbackPsfSize : `float`, optional
664 PSF FWHM to use in the event the exposure PSF cannot be retrieved.
665
666 Returns
667 -------
668 maglim : `astropy.units.Quantity`
669 The limiting magnitude of the exposure, or np.nan.
670 """
671 try:
672 psf = exposure.getPsf()
673 psf_shape = psf.computeShape(psf.getAveragePosition())
674 except (lsst.pex.exceptions.InvalidParameterError, afwDetection.InvalidPsfError):
675 if fallbackPsfSize is not None:
676 self.log.info("Unable to evaluate PSF, using fallback FWHM %f", fallbackPsfSize)
677 psf_area = np.pi*(fallbackPsfSize/2)**2
678 zeropoint = exposure.getPhotoCalib().instFluxToMagnitude(1)
679 maglim = zeropoint - 2.5*np.log10(nsigma*np.sqrt(psf_area))
680 else:
681 self.log.info("Unable to evaluate PSF, setting maglim to nan")
682 maglim = np.nan
683 else:
684 # Get a more accurate area than `psf_shape.getArea()` via moments
685 psf_area = np.pi*np.sqrt(psf_shape.getIxx()*psf_shape.getIyy())
686 zeropoint = exposure.getPhotoCalib().instFluxToMagnitude(1)
687 maglim = zeropoint - 2.5*np.log10(nsigma*np.sqrt(psf_area))
688 finally:
689 return maglim
690
691 @staticmethod
692 def _validateExposures(template, science):
693 """Check that the WCS of the two Exposures match, the template bbox
694 contains the science bbox, and that the bands match.
695
696 Parameters
697 ----------
698 template : `lsst.afw.image.ExposureF`
699 Template exposure, warped to match the science exposure.
700 science : `lsst.afw.image.ExposureF`
701 Science exposure to subtract from the template.
702
703 Raises
704 ------
705 AssertionError
706 Raised if the WCS of the template is not equal to the science WCS,
707 if the science image is not fully contained in the template
708 bounding box, or if the bands do not match.
709 """
710 assert template.wcs == science.wcs, \
711 "Template and science exposure WCS are not identical."
712 templateBBox = template.getBBox()
713 scienceBBox = science.getBBox()
714 assert science.filter.bandLabel == template.filter.bandLabel, \
715 "Science and template exposures have different bands: %s, %s" % \
716 (science.filter, template.filter)
717
718 assert templateBBox.contains(scienceBBox), \
719 "Template bbox does not contain all of the science image."
720
721 def _convolveExposure(self, exposure, kernel, convolutionControl,
722 bbox=None,
723 psf=None,
724 photoCalib=None,
725 interpolateBadMaskPlanes=False,
726 ):
727 """Convolve an exposure with the given kernel.
728
729 Parameters
730 ----------
731 exposure : `lsst.afw.Exposure`
732 exposure to convolve.
733 kernel : `lsst.afw.math.LinearCombinationKernel`
734 PSF matching kernel computed in the ``makeKernel`` subtask.
735 convolutionControl : `lsst.afw.math.ConvolutionControl`
736 Configuration for convolve algorithm.
737 bbox : `lsst.geom.Box2I`, optional
738 Bounding box to trim the convolved exposure to.
739 psf : `lsst.afw.detection.Psf`, optional
740 Point spread function (PSF) to set for the convolved exposure.
741 photoCalib : `lsst.afw.image.PhotoCalib`, optional
742 Photometric calibration of the convolved exposure.
743
744 Returns
745 -------
746 convolvedExp : `lsst.afw.Exposure`
747 The convolved image.
748 """
749 convolvedExposure = exposure.clone()
750 if psf is not None:
751 convolvedExposure.setPsf(psf)
752 if photoCalib is not None:
753 convolvedExposure.setPhotoCalib(photoCalib)
754 if interpolateBadMaskPlanes and self.config.badMaskPlanes is not None:
755 nInterp = _interpolateImage(convolvedExposure.maskedImage,
756 self.config.badMaskPlanes)
757 self.metadata.add("nInterpolated", nInterp)
758 convolvedImage = lsst.afw.image.MaskedImageF(convolvedExposure.getBBox())
759 lsst.afw.math.convolve(convolvedImage, convolvedExposure.maskedImage, kernel, convolutionControl)
760 convolvedExposure.setMaskedImage(convolvedImage)
761 if bbox is None:
762 return convolvedExposure
763 else:
764 return convolvedExposure[bbox]
765
766 def _sourceSelector(self, sources, mask):
767 """Select sources from a catalog that meet the selection criteria.
768
769 Parameters
770 ----------
771 sources : `lsst.afw.table.SourceCatalog`
772 Input source catalog to select sources from.
773 mask : `lsst.afw.image.Mask`
774 The image mask plane to use to reject sources
775 based on their location on the ccd.
776
777 Returns
778 -------
779 selectSources : `lsst.afw.table.SourceCatalog`
780 The input source catalog, with flagged and low signal-to-noise
781 sources removed.
782
783 Raises
784 ------
785 RuntimeError
786 If there are too few sources to compute the PSF matching kernel
787 remaining after source selection.
788 """
789
790 selected = self.sourceSelector.selectSources(sources).selected
791 nInitialSelected = np.count_nonzero(selected)
792 selected *= self._checkMask(mask, sources, self.config.excludeMaskPlanes)
793 nSelected = np.count_nonzero(selected)
794 self.log.info("Rejecting %i candidate sources: an excluded template mask plane is set.",
795 nInitialSelected - nSelected)
796 selectSources = sources[selected].copy(deep=True)
797 # Trim selectSources if they exceed ``maxKernelSources``.
798 # Keep the highest signal-to-noise sources of those selected.
799 if (len(selectSources) > self.config.maxKernelSources) & (self.config.maxKernelSources > 0):
800 signalToNoise = selectSources.getPsfInstFlux()/selectSources.getPsfInstFluxErr()
801 indices = np.argsort(signalToNoise)
802 indices = indices[-self.config.maxKernelSources:]
803 selected = np.zeros(len(selectSources), dtype=bool)
804 selected[indices] = True
805 selectSources = selectSources[selected].copy(deep=True)
806
807 self.log.info("%i/%i=%.1f%% of sources selected for PSF matching from the input catalog",
808 len(selectSources), len(sources), 100*len(selectSources)/len(sources))
809 if len(selectSources) < self.config.minKernelSources:
810 self.log.error("Too few sources to calculate the PSF matching kernel: "
811 "%i selected but %i needed for the calculation.",
812 len(selectSources), self.config.minKernelSources)
813 raise RuntimeError("Cannot compute PSF matching kernel: too few sources selected.")
814 self.metadata.add("nPsfSources", len(selectSources))
815
816 return selectSources
817
818 @staticmethod
819 def _checkMask(mask, sources, excludeMaskPlanes, checkAdjacent=True):
820 """Exclude sources that are located on or adjacent to masked pixels.
821
822 Parameters
823 ----------
824 mask : `lsst.afw.image.Mask`
825 The image mask plane to use to reject sources
826 based on the location of their centroid on the ccd.
827 sources : `lsst.afw.table.SourceCatalog`
828 The source catalog to evaluate.
829 excludeMaskPlanes : `list` of `str`
830 List of the names of the mask planes to exclude.
831
832 Returns
833 -------
834 flags : `numpy.ndarray` of `bool`
835 Array indicating whether each source in the catalog should be
836 kept (True) or rejected (False) based on the value of the
837 mask plane at its location.
838 """
839 setExcludeMaskPlanes = [
840 maskPlane for maskPlane in excludeMaskPlanes if maskPlane in mask.getMaskPlaneDict()
841 ]
842
843 excludePixelMask = mask.getPlaneBitMask(setExcludeMaskPlanes)
844
845 xv = (np.rint(sources.getX() - mask.getX0())).astype(int)
846 yv = (np.rint(sources.getY() - mask.getY0())).astype(int)
847
848 flags = np.ones(len(sources), dtype=bool)
849 if checkAdjacent:
850 pixRange = (0, -1, 1)
851 else:
852 pixRange = (0,)
853 for j in pixRange:
854 for i in pixRange:
855 mv = mask.array[yv + j, xv + i]
856 flags *= np.bitwise_and(mv, excludePixelMask) == 0
857 return flags
858
859 def _prepareInputs(self, template, science, visitSummary=None):
860 """Perform preparatory calculations common to all Alard&Lupton Tasks.
861
862 Parameters
863 ----------
864 template : `lsst.afw.image.ExposureF`
865 Template exposure, warped to match the science exposure. The
866 variance plane of the template image is modified in place.
867 science : `lsst.afw.image.ExposureF`
868 Science exposure to subtract from the template. The variance plane
869 of the science image is modified in place.
870 visitSummary : `lsst.afw.table.ExposureCatalog`, optional
871 Exposure catalog with external calibrations to be applied. Catalog
872 uses the detector id for the catalog id, sorted on id for fast
873 lookup.
874 """
875 self._validateExposures(template, science)
876 if visitSummary is not None:
877 self._applyExternalCalibrations(science, visitSummary=visitSummary)
878 templateCoverageFraction = checkTemplateIsSufficient(
879 template[science.getBBox()], self.log,
880 requiredTemplateFraction=self.config.requiredTemplateFraction,
881 exceptionMessage="Not attempting subtraction. To force subtraction,"
882 " set config requiredTemplateFraction=0"
883 )
884 self.metadata.add("templateCoveragePercent", 100*templateCoverageFraction)
885
886 if self.config.doScaleVariance:
887 # Scale the variance of the template and science images before
888 # convolution, subtraction, or decorrelation so that they have the
889 # correct ratio.
890 templateVarFactor = self.scaleVariance.run(template.maskedImage)
891 sciVarFactor = self.scaleVariance.run(science.maskedImage)
892 self.log.info("Template variance scaling factor: %.2f", templateVarFactor)
893 self.metadata.add("scaleTemplateVarianceFactor", templateVarFactor)
894 self.log.info("Science variance scaling factor: %.2f", sciVarFactor)
895 self.metadata.add("scaleScienceVarianceFactor", sciVarFactor)
896
897 # Erase existing detection mask planes.
898 # We don't want the detection mask from the science image
899 self.updateMasks(template, science)
900
901 def updateMasks(self, template, science):
902 """Update the science and template mask planes before differencing.
903
904 Parameters
905 ----------
906 template : `lsst.afw.image.Exposure`
907 Template exposure, warped to match the science exposure.
908 The template mask planes will be erased, except for a few specified
909 in the task config.
910 science : `lsst.afw.image.Exposure`
911 Science exposure to subtract from the template.
912 The DETECTED and DETECTED_NEGATIVE mask planes of the science image
913 will be erased.
914 """
915 self._clearMask(science.mask, clearMaskPlanes=["DETECTED", "DETECTED_NEGATIVE"])
916
917 # We will clear ALL template mask planes, except for those specified
918 # via the `preserveTemplateMask` config. Mask planes specified via
919 # the `renameTemplateMask` config will be copied to new planes with
920 # "_TEMPLATE" appended to their names, and the original mask plane will
921 # be cleared.
922 clearMaskPlanes = [mp for mp in template.mask.getMaskPlaneDict().keys()
923 if mp not in self.config.preserveTemplateMask]
924 renameMaskPlanes = [mp for mp in self.config.renameTemplateMask
925 if mp in template.mask.getMaskPlaneDict().keys()]
926
927 # propagate the mask plane related to Fake source injection
928 # NOTE: the fake source injection sets FAKE plane, but it should be INJECTED
929 # NOTE: This can be removed in DM-40796
930 if "FAKE" in science.mask.getMaskPlaneDict().keys():
931 self.log.info("Adding injected mask plane to science image")
932 self._renameMaskPlanes(science.mask, "FAKE", "INJECTED")
933 if "FAKE" in template.mask.getMaskPlaneDict().keys():
934 self.log.info("Adding injected mask plane to template image")
935 self._renameMaskPlanes(template.mask, "FAKE", "INJECTED_TEMPLATE")
936 if "INJECTED" in renameMaskPlanes:
937 renameMaskPlanes.remove("INJECTED")
938 if "INJECTED_TEMPLATE" in clearMaskPlanes:
939 clearMaskPlanes.remove("INJECTED_TEMPLATE")
940
941 for maskPlane in renameMaskPlanes:
942 self._renameMaskPlanes(template.mask, maskPlane, maskPlane + "_TEMPLATE")
943 self._clearMask(template.mask, clearMaskPlanes=clearMaskPlanes)
944
945 @staticmethod
946 def _renameMaskPlanes(mask, maskPlane, newMaskPlane):
947 """Rename a mask plane by adding the new name and copying the data.
948
949 Parameters
950 ----------
951 mask : `lsst.afw.image.Mask`
952 The mask image to update in place.
953 maskPlane : `str`
954 The name of the existing mask plane to copy.
955 newMaskPlane : `str`
956 The new name of the mask plane that will be added.
957 If the mask plane already exists, it will be updated in place.
958 """
959 mask.addMaskPlane(newMaskPlane)
960 originBitMask = mask.getPlaneBitMask(maskPlane)
961 destinationBitMask = mask.getPlaneBitMask(newMaskPlane)
962 mask.array |= ((mask.array & originBitMask) > 0)*destinationBitMask
963
964 def _clearMask(self, mask, clearMaskPlanes=None):
965 """Clear the mask plane of an exposure.
966
967 Parameters
968 ----------
969 mask : `lsst.afw.image.Mask`
970 The mask plane to erase, which will be modified in place.
971 clearMaskPlanes : `list` of `str`, optional
972 Erase the specified mask planes.
973 If not supplied, the entire mask will be erased.
974 """
975 if clearMaskPlanes is None:
976 clearMaskPlanes = list(mask.getMaskPlaneDict().keys())
977
978 bitMaskToClear = mask.getPlaneBitMask(clearMaskPlanes)
979 mask &= ~bitMaskToClear
980
981
983 SubtractScoreOutputConnections):
984 pass
985
986
988 pipelineConnections=AlardLuptonPreconvolveSubtractConnections):
989 pass
990
991
993 """Subtract a template from a science image, convolving the science image
994 before computing the kernel, and also convolving the template before
995 subtraction.
996 """
997 ConfigClass = AlardLuptonPreconvolveSubtractConfig
998 _DefaultName = "alardLuptonPreconvolveSubtract"
999
1000 def run(self, template, science, sources, visitSummary=None):
1001 """Preconvolve the science image with its own PSF,
1002 convolve the template image with a PSF-matching kernel and subtract
1003 from the preconvolved science image.
1004
1005 Parameters
1006 ----------
1007 template : `lsst.afw.image.ExposureF`
1008 The template image, which has previously been warped to the science
1009 image. The template bbox will be padded by a few pixels compared to
1010 the science bbox.
1011 science : `lsst.afw.image.ExposureF`
1012 The science exposure.
1013 sources : `lsst.afw.table.SourceCatalog`
1014 Identified sources on the science exposure. This catalog is used to
1015 select sources in order to perform the AL PSF matching on stamp
1016 images around them.
1017 visitSummary : `lsst.afw.table.ExposureCatalog`, optional
1018 Exposure catalog with complete external calibrations. Catalog uses
1019 the detector id for the catalog id, sorted on id for fast lookup.
1020
1021 Returns
1022 -------
1023 results : `lsst.pipe.base.Struct`
1024 ``scoreExposure`` : `lsst.afw.image.ExposureF`
1025 Result of subtracting the convolved template and science
1026 images. Attached PSF is that of the original science image.
1027 ``matchedTemplate`` : `lsst.afw.image.ExposureF`
1028 Warped and PSF-matched template exposure. Attached PSF is that
1029 of the original science image.
1030 ``matchedScience`` : `lsst.afw.image.ExposureF`
1031 The science exposure after convolving with its own PSF.
1032 Attached PSF is that of the original science image.
1033 ``backgroundModel`` : `lsst.afw.math.Function2D`
1034 Background model that was fit while solving for the
1035 PSF-matching kernel
1036 ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
1037 Final kernel used to PSF-match the template to the science
1038 image.
1039 """
1040 self._prepareInputs(template, science, visitSummary=visitSummary)
1041
1042 # TODO: DM-37212 we need to mirror the kernel in order to get correct cross correlation
1043 scienceKernel = science.psf.getKernel()
1044 matchedScience = self._convolveExposure(science, scienceKernel, self.convolutionControlconvolutionControl,
1045 interpolateBadMaskPlanes=True)
1046 self.metadata.add("convolvedExposure", "Preconvolution")
1047 try:
1048 selectSources = self._sourceSelector(sources, matchedScience.mask)
1049 subtractResults = self.runPreconvolve(template, science, matchedScience,
1050 selectSources, scienceKernel)
1051
1052 except (RuntimeError, lsst.pex.exceptions.Exception) as e:
1053 self.loglog.warning("Failed to match template. Checking coverage")
1054 # Raise NoWorkFound if template fraction is insufficient
1055 checkTemplateIsSufficient(template[science.getBBox()], self.loglog,
1056 self.config.minTemplateFractionForExpectedSuccess,
1057 exceptionMessage="Template coverage lower than expected to succeed."
1058 f" Failure is tolerable: {e}")
1059 # checkTemplateIsSufficient did not raise NoWorkFound, so raise original exception
1060 raise e
1061
1062 return subtractResults
1063
1064 def runPreconvolve(self, template, science, matchedScience, selectSources, preConvKernel):
1065 """Convolve the science image with its own PSF, then convolve the
1066 template with a matching kernel and subtract to form the Score
1067 exposure.
1068
1069 Parameters
1070 ----------
1071 template : `lsst.afw.image.ExposureF`
1072 Template exposure, warped to match the science exposure.
1073 science : `lsst.afw.image.ExposureF`
1074 Science exposure to subtract from the template.
1075 matchedScience : `lsst.afw.image.ExposureF`
1076 The science exposure, convolved with the reflection of its own PSF.
1077 selectSources : `lsst.afw.table.SourceCatalog`
1078 Identified sources on the science exposure. This catalog is used to
1079 select sources in order to perform the AL PSF matching on stamp
1080 images around them.
1081 preConvKernel : `lsst.afw.math.Kernel`
1082 The reflection of the kernel that was used to preconvolve the
1083 `science` exposure. Must be normalized to sum to 1.
1084
1085 Returns
1086 -------
1087 results : `lsst.pipe.base.Struct`
1088
1089 ``scoreExposure`` : `lsst.afw.image.ExposureF`
1090 Result of subtracting the convolved template and science
1091 images. Attached PSF is that of the original science image.
1092 ``matchedTemplate`` : `lsst.afw.image.ExposureF`
1093 Warped and PSF-matched template exposure. Attached PSF is that
1094 of the original science image.
1095 ``matchedScience`` : `lsst.afw.image.ExposureF`
1096 The science exposure after convolving with its own PSF.
1097 Attached PSF is that of the original science image.
1098 ``backgroundModel`` : `lsst.afw.math.Function2D`
1099 Background model that was fit while solving for the
1100 PSF-matching kernel
1101 ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
1102 Final kernel used to PSF-match the template to the science
1103 image.
1104 """
1105 bbox = science.getBBox()
1106 innerBBox = preConvKernel.shrinkBBox(bbox)
1107
1108 kernelSources = self.makeKernel.selectKernelSources(template[innerBBox], matchedScience[innerBBox],
1109 candidateList=selectSources,
1110 preconvolved=True)
1111 kernelResult = self.makeKernel.run(template[innerBBox], matchedScience[innerBBox], kernelSources,
1112 preconvolved=True)
1113
1114 matchedTemplate = self._convolveExposure(template, kernelResult.psfMatchingKernel,
1116 bbox=bbox,
1117 psf=science.psf,
1118 interpolateBadMaskPlanes=True,
1119 photoCalib=science.photoCalib)
1120 score = _subtractImages(matchedScience, matchedTemplate,
1121 backgroundModel=(kernelResult.backgroundModel
1122 if self.config.doSubtractBackground else None))
1123 correctedScore = self.finalize(template[bbox], science, score,
1124 kernelResult.psfMatchingKernel,
1125 templateMatched=True, preConvMode=True,
1126 preConvKernel=preConvKernel)
1127
1128 return lsst.pipe.base.Struct(scoreExposure=correctedScore,
1129 matchedTemplate=matchedTemplate,
1130 matchedScience=matchedScience,
1131 backgroundModel=kernelResult.backgroundModel,
1132 psfMatchingKernel=kernelResult.psfMatchingKernel)
1133
1134
1135def checkTemplateIsSufficient(templateExposure, logger, requiredTemplateFraction=0.,
1136 exceptionMessage=""):
1137 """Raise NoWorkFound if template coverage < requiredTemplateFraction
1138
1139 Parameters
1140 ----------
1141 templateExposure : `lsst.afw.image.ExposureF`
1142 The template exposure to check
1143 logger : `logging.Logger`
1144 Logger for printing output.
1145 requiredTemplateFraction : `float`, optional
1146 Fraction of pixels of the science image required to have coverage
1147 in the template.
1148 exceptionMessage : `str`, optional
1149 Message to include in the exception raised if the template coverage
1150 is insufficient.
1151
1152 Returns
1153 -------
1154 templateCoverageFraction: `float`
1155 Fraction of pixels in the template with data.
1156
1157 Raises
1158 ------
1159 lsst.pipe.base.NoWorkFound
1160 Raised if fraction of good pixels, defined as not having NO_DATA
1161 set, is less than the requiredTemplateFraction
1162 """
1163 # Count the number of pixels with the NO_DATA mask bit set
1164 # counting NaN pixels is insufficient because pixels without data are often intepolated over)
1165 pixNoData = np.count_nonzero(templateExposure.mask.array
1166 & templateExposure.mask.getPlaneBitMask('NO_DATA'))
1167 pixGood = templateExposure.getBBox().getArea() - pixNoData
1168 templateCoverageFraction = pixGood/templateExposure.getBBox().getArea()
1169 logger.info("template has %d good pixels (%.1f%%)", pixGood, 100*templateCoverageFraction)
1170
1171 if templateCoverageFraction < requiredTemplateFraction:
1172 message = ("Insufficient Template Coverage. (%.1f%% < %.1f%%)" % (
1173 100*templateCoverageFraction,
1174 100*requiredTemplateFraction))
1175 raise lsst.pipe.base.NoWorkFound(message + " " + exceptionMessage)
1176 return templateCoverageFraction
1177
1178
1179def _subtractImages(science, template, backgroundModel=None):
1180 """Subtract template from science, propagating relevant metadata.
1181
1182 Parameters
1183 ----------
1184 science : `lsst.afw.Exposure`
1185 The input science image.
1186 template : `lsst.afw.Exposure`
1187 The template to subtract from the science image.
1188 backgroundModel : `lsst.afw.MaskedImage`, optional
1189 Differential background model
1190
1191 Returns
1192 -------
1193 difference : `lsst.afw.Exposure`
1194 The subtracted image.
1195 """
1196 difference = science.clone()
1197 if backgroundModel is not None:
1198 difference.maskedImage -= backgroundModel
1199 difference.maskedImage -= template.maskedImage
1200 return difference
1201
1202
1203def _shapeTest(exp1, exp2, fwhmExposureBuffer, fwhmExposureGrid):
1204 """Determine that the PSF of ``exp1`` is not wider than that of ``exp2``.
1205
1206 Parameters
1207 ----------
1208 exp1 : `~lsst.afw.image.Exposure`
1209 Exposure with the reference point spread function (PSF) to evaluate.
1210 exp2 : `~lsst.afw.image.Exposure`
1211 Exposure with a candidate point spread function (PSF) to evaluate.
1212 fwhmExposureBuffer : `float`
1213 Fractional buffer margin to be left out of all sides of the image
1214 during the construction of the grid to compute mean PSF FWHM in an
1215 exposure, if the PSF is not available at its average position.
1216 fwhmExposureGrid : `int`
1217 Grid size to compute the mean FWHM in an exposure, if the PSF is not
1218 available at its average position.
1219 Returns
1220 -------
1221 result : `bool`
1222 True if ``exp1`` has a PSF that is not wider than that of ``exp2`` in
1223 either dimension.
1224 """
1225 try:
1226 shape1 = getPsfFwhm(exp1.psf, average=False)
1227 shape2 = getPsfFwhm(exp2.psf, average=False)
1229 shape1 = evaluateMeanPsfFwhm(exp1,
1230 fwhmExposureBuffer=fwhmExposureBuffer,
1231 fwhmExposureGrid=fwhmExposureGrid
1232 )
1233 shape2 = evaluateMeanPsfFwhm(exp2,
1234 fwhmExposureBuffer=fwhmExposureBuffer,
1235 fwhmExposureGrid=fwhmExposureGrid
1236 )
1237 return shape1 <= shape2
1238
1239 # Results from getPsfFwhm is a tuple of two values, one for each dimension.
1240 xTest = shape1[0] <= shape2[0]
1241 yTest = shape1[1] <= shape2[1]
1242 return xTest | yTest
1243
1244
1245def _interpolateImage(maskedImage, badMaskPlanes, fallbackValue=None):
1246 """Replace masked image pixels with interpolated values.
1247
1248 Parameters
1249 ----------
1250 maskedImage : `lsst.afw.image.MaskedImage`
1251 Image on which to perform interpolation.
1252 badMaskPlanes : `list` of `str`
1253 List of mask planes to interpolate over.
1254 fallbackValue : `float`, optional
1255 Value to set when interpolation fails.
1256
1257 Returns
1258 -------
1259 result: `float`
1260 The number of masked pixels that were replaced.
1261 """
1262 imgBadMaskPlanes = [
1263 maskPlane for maskPlane in badMaskPlanes if maskPlane in maskedImage.mask.getMaskPlaneDict()
1264 ]
1265
1266 image = maskedImage.image.array
1267 badPixels = (maskedImage.mask.array & maskedImage.mask.getPlaneBitMask(imgBadMaskPlanes)) > 0
1268 image[badPixels] = np.nan
1269 if fallbackValue is None:
1270 fallbackValue = np.nanmedian(image)
1271 # For this initial implementation, skip the interpolation and just fill with
1272 # the median value.
1273 image[badPixels] = fallbackValue
1274 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)
run(self, template, science, sources, visitSummary=None)
_calculateMagLim(self, exposure, nsigma=5.0, fallbackPsfSize=None)
runConvolveTemplate(self, template, science, selectSources)
_checkMask(mask, sources, excludeMaskPlanes, checkAdjacent=True)
_convolveExposure(self, exposure, kernel, convolutionControl, bbox=None, psf=None, photoCalib=None, interpolateBadMaskPlanes=False)
runConvolveScience(self, template, science, selectSources)
finalize(self, template, science, difference, kernel, templateMatched=True, preConvMode=False, preConvKernel=None, spatiallyVarying=False)
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)
_interpolateImage(maskedImage, badMaskPlanes, fallbackValue=None)
checkTemplateIsSufficient(templateExposure, logger, requiredTemplateFraction=0., exceptionMessage="")
_shapeTest(exp1, exp2, fwhmExposureBuffer, fwhmExposureGrid)