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