LSST Applications g0000d66e7c+ce78115f25,g0485b4d2cb+c8d56b10d4,g0fba68d861+fcbc158cd0,g1ec0fe41b4+3e153770da,g1fd858c14a+57ee4e1624,g2440f9efcc+8c5ae1fdc5,g35bb328faa+8c5ae1fdc5,g4d2262a081+1e04cc5a47,g53246c7159+8c5ae1fdc5,g55585698de+7a33f081c8,g56a49b3a55+b9d5cac73f,g60b5630c4e+7a33f081c8,g67b6fd64d1+035c836e50,g78460c75b0+7e33a9eb6d,g786e29fd12+668abc6043,g7ac00fbb6c+b938379438,g8352419a5c+8c5ae1fdc5,g8852436030+5ba78a36c9,g89139ef638+035c836e50,g94187f82dc+7a33f081c8,g989de1cb63+035c836e50,g9d31334357+7a33f081c8,g9f33ca652e+e34120223a,ga815be3f0b+911242149a,gabe3b4be73+8856018cbb,gabf8522325+21619da9f3,gb1101e3267+0b44b44611,gb89ab40317+035c836e50,gc91f06edcd+e59fb3c9bc,gcf25f946ba+5ba78a36c9,gd6cbbdb0b4+958adf5c1f,gde0f65d7ad+6c98dcc924,ge278dab8ac+83c63f4893,ge410e46f29+035c836e50,gf35d7ec915+97dd712d81,gf5e32f922b+8c5ae1fdc5,gf67bdafdda+035c836e50,gf6800124b1+1714c04baa,w.2025.19
LSST Data Management Base Package
Loading...
Searching...
No Matches
subtractImages.py
Go to the documentation of this file.
1# This file is part of ip_diffim.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
21
22from astropy import units as u
23import numpy as np
24
25import lsst.afw.detection as afwDetection
26import lsst.afw.image
27import lsst.afw.math
28import lsst.geom
29from lsst.ip.diffim.utils import evaluateMeanPsfFwhm, getPsfFwhm
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=False,
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 # Always include background fitting in the kernel fit,
246 # even if it is not subtracted
247 self.makeKernel.kernel.active.fitForBackground = True
248 self.makeKernel.kernel.active.spatialKernelOrder = 1
249 self.makeKernel.kernel.active.spatialBgOrder = 2
250 self.sourceSelector.doUnresolved = True # apply star-galaxy separation
251 self.sourceSelector.doIsolated = True # apply isolated star selection
252 self.sourceSelector.doRequirePrimary = True # apply primary flag selection
253 self.sourceSelector.doSkySources = False # Do not include sky sources
254 self.sourceSelector.doSignalToNoise = True # apply signal to noise filter
255 self.sourceSelector.signalToNoise.minimum = 10
256 self.sourceSelector.signalToNoise.maximum = 500
257
258
259class AlardLuptonSubtractConfig(AlardLuptonSubtractBaseConfig, lsst.pipe.base.PipelineTaskConfig,
260 pipelineConnections=AlardLuptonSubtractConnections):
262 dtype=str,
263 default="convolveTemplate",
264 allowed={"auto": "Choose which image to convolve at runtime.",
265 "convolveScience": "Only convolve the science image.",
266 "convolveTemplate": "Only convolve the template image."},
267 doc="Choose which image to convolve at runtime, or require that a specific image is convolved."
268 )
269
270
271class AlardLuptonSubtractTask(lsst.pipe.base.PipelineTask):
272 """Compute the image difference of a science and template image using
273 the Alard & Lupton (1998) algorithm.
274 """
275 ConfigClass = AlardLuptonSubtractConfig
276 _DefaultName = "alardLuptonSubtract"
277
278 def __init__(self, **kwargs):
279 super().__init__(**kwargs)
280 self.makeSubtask("decorrelate")
281 self.makeSubtask("makeKernel")
282 self.makeSubtask("sourceSelector")
283 if self.config.doScaleVariance:
284 self.makeSubtask("scaleVariance")
285
287 # Normalization is an extra, unnecessary, calculation and will result
288 # in mis-subtraction of the images if there are calibration errors.
289 self.convolutionControl.setDoNormalize(False)
290 self.convolutionControl.setDoCopyEdge(True)
291
292 def _applyExternalCalibrations(self, exposure, visitSummary):
293 """Replace calibrations (psf, and ApCorrMap) on this exposure with
294 external ones.".
295
296 Parameters
297 ----------
298 exposure : `lsst.afw.image.exposure.Exposure`
299 Input exposure to adjust calibrations.
300 visitSummary : `lsst.afw.table.ExposureCatalog`
301 Exposure catalog with external calibrations to be applied. Catalog
302 uses the detector id for the catalog id, sorted on id for fast
303 lookup.
304
305 Returns
306 -------
307 exposure : `lsst.afw.image.exposure.Exposure`
308 Exposure with adjusted calibrations.
309 """
310 detectorId = exposure.info.getDetector().getId()
311
312 row = visitSummary.find(detectorId)
313 if row is None:
314 self.log.warning("Detector id %s not found in external calibrations catalog; "
315 "Using original calibrations.", detectorId)
316 else:
317 psf = row.getPsf()
318 apCorrMap = row.getApCorrMap()
319 if psf is None:
320 self.log.warning("Detector id %s has None for psf in "
321 "external calibrations catalog; Using original psf and aperture correction.",
322 detectorId)
323 elif apCorrMap is None:
324 self.log.warning("Detector id %s has None for apCorrMap in "
325 "external calibrations catalog; Using original psf and aperture correction.",
326 detectorId)
327 else:
328 exposure.setPsf(psf)
329 exposure.info.setApCorrMap(apCorrMap)
330
331 return exposure
332
333 @timeMethod
334 def run(self, template, science, sources, visitSummary=None):
335 """PSF match, subtract, and decorrelate two images.
336
337 Parameters
338 ----------
339 template : `lsst.afw.image.ExposureF`
340 Template exposure, warped to match the science exposure.
341 science : `lsst.afw.image.ExposureF`
342 Science exposure to subtract from the template.
343 sources : `lsst.afw.table.SourceCatalog`
344 Identified sources on the science exposure. This catalog is used to
345 select sources in order to perform the AL PSF matching on stamp
346 images around them.
347 visitSummary : `lsst.afw.table.ExposureCatalog`, optional
348 Exposure catalog with external calibrations to be applied. Catalog
349 uses the detector id for the catalog id, sorted on id for fast
350 lookup.
351
352 Returns
353 -------
354 results : `lsst.pipe.base.Struct`
355 ``difference`` : `lsst.afw.image.ExposureF`
356 Result of subtracting template and science.
357 ``matchedTemplate`` : `lsst.afw.image.ExposureF`
358 Warped and PSF-matched template exposure.
359 ``backgroundModel`` : `lsst.afw.math.Function2D`
360 Background model that was fit while solving for the
361 PSF-matching kernel
362 ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
363 Kernel used to PSF-match the convolved image.
364
365 Raises
366 ------
367 RuntimeError
368 If an unsupported convolution mode is supplied.
369 RuntimeError
370 If there are too few sources to calculate the PSF matching kernel.
371 lsst.pipe.base.NoWorkFound
372 Raised if fraction of good pixels, defined as not having NO_DATA
373 set, is less then the configured requiredTemplateFraction
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 self.templatePsfSize = getPsfFwhm(template.psf)
394 self.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 self.templatePsfSize = evaluateMeanPsfFwhm(template,
400 fwhmExposureBuffer=fwhmExposureBuffer,
401 fwhmExposureGrid=fwhmExposureGrid
402 )
403 self.sciencePsfSize = evaluateMeanPsfFwhm(science,
404 fwhmExposureBuffer=fwhmExposureBuffer,
405 fwhmExposureGrid=fwhmExposureGrid
406 )
407 self.log.info("Science PSF FWHM: %f pixels", self.sciencePsfSize)
408 self.log.info("Template PSF FWHM: %f pixels", self.templatePsfSize)
409 self.metadata["sciencePsfSize"] = self.sciencePsfSize
410 self.metadata["templatePsfSize"] = self.templatePsfSize
411
412 # Calculate estimated image depths, i.e., limiting magnitudes
413 maglim_science = self._calculateMagLim(science, fallbackPsfSize=self.sciencePsfSize)
414 fluxlim_science = (maglim_science*u.ABmag).to_value(u.nJy)
415 maglim_template = self._calculateMagLim(template, fallbackPsfSize=self.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["scienceLimitingMagnitude"] = maglim_science
423 self.metadata["templateLimitingMagnitude"] = maglim_template
424 self.metadata["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 self.sciencePsfSize < self.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["convolvedExposure"] = "Template"
454 subtractResults = self.runConvolveTemplate(template, science, selectSources)
455 else:
456 self.metadata["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 self.computeImageMetrics(science, subtractResults.difference, sources)
470
471 return subtractResults
472
473 def computeImageMetrics(self, science, difference, stars):
474 r"""Compute quality metrics (saved to the task metadata) on the
475 difference image, at the locations of detected stars on the science
476 image. This restricts the metric to locations that should be
477 well-subtracted.
478
479 Parameters
480 ----------
481 science : `lsst.afw.image.ExposureF`
482 Science exposure that was subtracted.
483 difference : `lsst.afw.image.ExposureF`
484 Result of subtracting template and science.
485 stars : `lsst.afw.table.SourceCatalog`
486 Good calibration sources detected on science image; these
487 footprints are what the metrics are computed on.
488
489 Notes
490 -----
491 The task metadata does not include docstrings, so descriptions of the
492 computed metrics are given here:
493
494 differenceFootprintRatioMean
495 Mean of the ratio of the absolute value of the difference image
496 (with the mean absolute value of the sky regions on the difference
497 image removed) to the science image, computed in the footprints
498 of stars detected on the science image (the sums below are of the
499 pixels in each star or sky footprint):
500 :math:`\mathrm{mean}_{footprints}((\sum |difference| -
501 \mathrm{mean}(\sum |difference_{sky}|)) / \sum science)`
502 differenceFootprintRatioStdev
503 Standard Deviation across footprints of the above ratio.
504 differenceFootprintSkyRatioMean
505 Mean of the ratio of the absolute value of sky source regions on
506 the difference image to the science image (the sum below is of the
507 pixels in each sky source footprint):
508 :math:`\mathrm{mean}_{footprints}(\sum |difference_{sky}| / \sum science_{sky})`
509 differenceFootprintSkyRatioStdev
510 Standard Deivation across footprints of the above sky ratio.
511 """
512 def footprint_mean(sources, sky=0):
513 """Compute ratio of the absolute value of the diffim to the science
514 image, within each source footprint, subtracting the sky from the
515 diffim values if provided.
516 """
517 n = len(sources)
518 science_footprints = np.zeros(n)
519 difference_footprints = np.zeros(n)
520 ratio = np.zeros(n)
521 for i, record in enumerate(sources):
522 footprint = record.getFootprint()
523 heavy = lsst.afw.detection.makeHeavyFootprint(footprint, science.maskedImage)
524 heavy_diff = lsst.afw.detection.makeHeavyFootprint(footprint, difference.maskedImage)
525 science_footprints[i] = abs(heavy.getImageArray()).mean()
526 difference_footprints[i] = abs(heavy_diff.getImageArray()).mean()
527 ratio[i] = abs((difference_footprints[i] - sky)/science_footprints[i])
528 return science_footprints, difference_footprints, ratio
529
530 sky = stars["sky_source"]
531 sky_science, sky_difference, sky_ratio = footprint_mean(stars[sky])
532 science_footprints, difference_footprints, ratio = footprint_mean(stars[~sky], sky_difference.mean())
533
534 self.metadata["differenceFootprintRatioMean"] = ratio.mean()
535 self.metadata["differenceFootprintRatioStdev"] = ratio.std()
536 self.metadata["differenceFootprintSkyRatioMean"] = sky_ratio.mean()
537 self.metadata["differenceFootprintSkyRatioStdev"] = sky_ratio.std()
538 self.log.info("Mean, stdev of ratio of difference to science "
539 "pixels in star footprints: %5.4f, %5.4f",
540 self.metadata["differenceFootprintRatioMean"],
541 self.metadata["differenceFootprintRatioStdev"])
542
543 def runConvolveTemplate(self, template, science, selectSources):
544 """Convolve the template image with a PSF-matching kernel and subtract
545 from the science image.
546
547 Parameters
548 ----------
549 template : `lsst.afw.image.ExposureF`
550 Template exposure, warped to match the science exposure.
551 science : `lsst.afw.image.ExposureF`
552 Science exposure to subtract from the template.
553 selectSources : `lsst.afw.table.SourceCatalog`
554 Identified sources on the science exposure. This catalog is used to
555 select sources in order to perform the AL PSF matching on stamp
556 images around them.
557
558 Returns
559 -------
560 results : `lsst.pipe.base.Struct`
561
562 ``difference`` : `lsst.afw.image.ExposureF`
563 Result of subtracting template and science.
564 ``matchedTemplate`` : `lsst.afw.image.ExposureF`
565 Warped and PSF-matched template exposure.
566 ``backgroundModel`` : `lsst.afw.math.Function2D`
567 Background model that was fit while solving for the PSF-matching kernel
568 ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
569 Kernel used to PSF-match the template to the science image.
570 """
571 try:
572 kernelSources = self.makeKernel.selectKernelSources(template, science,
573 candidateList=selectSources,
574 preconvolved=False)
575 kernelResult = self.makeKernel.run(template, science, kernelSources,
576 preconvolved=False)
577 except Exception as e:
578 if self.config.allowKernelSourceDetection:
579 self.log.warning("Error encountered trying to construct the matching kernel"
580 f" Running source detection and retrying. {e}")
581 kernelSize = self.makeKernel.makeKernelBasisList(
582 self.templatePsfSize, self.sciencePsfSize)[0].getWidth()
583 sigmaToFwhm = 2*np.log(2*np.sqrt(2))
584 candidateList = self.makeKernel.makeCandidateList(template, science, kernelSize,
585 candidateList=None,
586 sigma=self.sciencePsfSize/sigmaToFwhm)
587 kernelSources = self.makeKernel.selectKernelSources(template, science,
588 candidateList=candidateList,
589 preconvolved=False)
590 kernelResult = self.makeKernel.run(template, science, kernelSources,
591 preconvolved=False)
592 else:
593 raise e
594
595 matchedTemplate = self._convolveExposure(template, kernelResult.psfMatchingKernel,
597 bbox=science.getBBox(),
598 psf=science.psf,
599 photoCalib=science.photoCalib)
600
601 difference = _subtractImages(science, matchedTemplate,
602 backgroundModel=(kernelResult.backgroundModel
603 if self.config.doSubtractBackground else None))
604 correctedExposure = self.finalize(template, science, difference,
605 kernelResult.psfMatchingKernel,
606 templateMatched=True)
607
608 return lsst.pipe.base.Struct(difference=correctedExposure,
609 matchedTemplate=matchedTemplate,
610 matchedScience=science,
611 backgroundModel=kernelResult.backgroundModel,
612 psfMatchingKernel=kernelResult.psfMatchingKernel,
613 kernelSources=kernelSources)
614
615 def runConvolveScience(self, template, science, selectSources):
616 """Convolve the science image with a PSF-matching kernel and subtract
617 the template image.
618
619 Parameters
620 ----------
621 template : `lsst.afw.image.ExposureF`
622 Template exposure, warped to match the science exposure.
623 science : `lsst.afw.image.ExposureF`
624 Science exposure to subtract from the template.
625 selectSources : `lsst.afw.table.SourceCatalog`
626 Identified sources on the science exposure. This catalog is used to
627 select sources in order to perform the AL PSF matching on stamp
628 images around them.
629
630 Returns
631 -------
632 results : `lsst.pipe.base.Struct`
633
634 ``difference`` : `lsst.afw.image.ExposureF`
635 Result of subtracting template and science.
636 ``matchedTemplate`` : `lsst.afw.image.ExposureF`
637 Warped template exposure. Note that in this case, the template
638 is not PSF-matched to the science image.
639 ``backgroundModel`` : `lsst.afw.math.Function2D`
640 Background model that was fit while solving for the PSF-matching kernel
641 ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
642 Kernel used to PSF-match the science image to the template.
643 """
644 bbox = science.getBBox()
645 kernelSources = self.makeKernel.selectKernelSources(science, template,
646 candidateList=selectSources,
647 preconvolved=False)
648 kernelResult = self.makeKernel.run(science, template, kernelSources,
649 preconvolved=False)
650 modelParams = kernelResult.backgroundModel.getParameters()
651 # We must invert the background model if the matching kernel is solved for the science image.
652 kernelResult.backgroundModel.setParameters([-p for p in modelParams])
653
654 kernelImage = lsst.afw.image.ImageD(kernelResult.psfMatchingKernel.getDimensions())
655 norm = kernelResult.psfMatchingKernel.computeImage(kernelImage, doNormalize=False)
656
657 matchedScience = self._convolveExposure(science, kernelResult.psfMatchingKernel,
659 psf=template.psf)
660
661 # Place back on native photometric scale
662 matchedScience.maskedImage /= norm
663 matchedTemplate = template.clone()[bbox]
664 matchedTemplate.maskedImage /= norm
665 matchedTemplate.setPhotoCalib(science.photoCalib)
666
667 difference = _subtractImages(matchedScience, matchedTemplate,
668 backgroundModel=(kernelResult.backgroundModel
669 if self.config.doSubtractBackground else None))
670
671 correctedExposure = self.finalize(template, science, difference,
672 kernelResult.psfMatchingKernel,
673 templateMatched=False)
674
675 return lsst.pipe.base.Struct(difference=correctedExposure,
676 matchedTemplate=matchedTemplate,
677 matchedScience=matchedScience,
678 backgroundModel=kernelResult.backgroundModel,
679 psfMatchingKernel=kernelResult.psfMatchingKernel,
680 kernelSources=kernelSources)
681
682 def finalize(self, template, science, difference, kernel,
683 templateMatched=True,
684 preConvMode=False,
685 preConvKernel=None,
686 spatiallyVarying=False):
687 """Decorrelate the difference image to undo the noise correlations
688 caused by convolution.
689
690 Parameters
691 ----------
692 template : `lsst.afw.image.ExposureF`
693 Template exposure, warped to match the science exposure.
694 science : `lsst.afw.image.ExposureF`
695 Science exposure to subtract from the template.
696 difference : `lsst.afw.image.ExposureF`
697 Result of subtracting template and science.
698 kernel : `lsst.afw.math.Kernel`
699 An (optionally spatially-varying) PSF matching kernel
700 templateMatched : `bool`, optional
701 Was the template PSF-matched to the science image?
702 preConvMode : `bool`, optional
703 Was the science image preconvolved with its own PSF
704 before PSF matching the template?
705 preConvKernel : `lsst.afw.detection.Psf`, optional
706 If not `None`, then the science image was pre-convolved with
707 (the reflection of) this kernel. Must be normalized to sum to 1.
708 spatiallyVarying : `bool`, optional
709 Compute the decorrelation kernel spatially varying across the image?
710
711 Returns
712 -------
713 correctedExposure : `lsst.afw.image.ExposureF`
714 The decorrelated image difference.
715 """
716 if self.config.doDecorrelation:
717 self.log.info("Decorrelating image difference.")
718 # We have cleared the template mask plane, so copy the mask plane of
719 # the image difference so that we can calculate correct statistics
720 # during decorrelation
721 correctedExposure = self.decorrelate.run(science, template[science.getBBox()], difference, kernel,
722 templateMatched=templateMatched,
723 preConvMode=preConvMode,
724 preConvKernel=preConvKernel,
725 spatiallyVarying=spatiallyVarying).correctedExposure
726 else:
727 self.log.info("NOT decorrelating image difference.")
728 correctedExposure = difference
729 return correctedExposure
730
731 def _calculateMagLim(self, exposure, nsigma=5.0, fallbackPsfSize=None):
732 """Calculate an exposure's limiting magnitude.
733
734 This method uses the photometric zeropoint together with the
735 PSF size from the average position of the exposure.
736
737 Parameters
738 ----------
739 exposure : `lsst.afw.image.Exposure`
740 The target exposure to calculate the limiting magnitude for.
741 nsigma : `float`, optional
742 The detection threshold in sigma.
743 fallbackPsfSize : `float`, optional
744 PSF FWHM to use in the event the exposure PSF cannot be retrieved.
745
746 Returns
747 -------
748 maglim : `astropy.units.Quantity`
749 The limiting magnitude of the exposure, or np.nan.
750 """
751 try:
752 psf = exposure.getPsf()
753 psf_shape = psf.computeShape(psf.getAveragePosition())
754 except (lsst.pex.exceptions.InvalidParameterError, afwDetection.InvalidPsfError):
755 if fallbackPsfSize is not None:
756 self.log.info("Unable to evaluate PSF, using fallback FWHM %f", fallbackPsfSize)
757 psf_area = np.pi*(fallbackPsfSize/2)**2
758 zeropoint = exposure.getPhotoCalib().instFluxToMagnitude(1)
759 maglim = zeropoint - 2.5*np.log10(nsigma*np.sqrt(psf_area))
760 else:
761 self.log.info("Unable to evaluate PSF, setting maglim to nan")
762 maglim = np.nan
763 else:
764 # Get a more accurate area than `psf_shape.getArea()` via moments
765 psf_area = np.pi*np.sqrt(psf_shape.getIxx()*psf_shape.getIyy())
766 zeropoint = exposure.getPhotoCalib().instFluxToMagnitude(1)
767 maglim = zeropoint - 2.5*np.log10(nsigma*np.sqrt(psf_area))
768 finally:
769 return maglim
770
771 @staticmethod
772 def _validateExposures(template, science):
773 """Check that the WCS of the two Exposures match, the template bbox
774 contains the science bbox, and that the bands match.
775
776 Parameters
777 ----------
778 template : `lsst.afw.image.ExposureF`
779 Template exposure, warped to match the science exposure.
780 science : `lsst.afw.image.ExposureF`
781 Science exposure to subtract from the template.
782
783 Raises
784 ------
785 AssertionError
786 Raised if the WCS of the template is not equal to the science WCS,
787 if the science image is not fully contained in the template
788 bounding box, or if the bands do not match.
789 """
790 assert template.wcs == science.wcs, \
791 "Template and science exposure WCS are not identical."
792 templateBBox = template.getBBox()
793 scienceBBox = science.getBBox()
794 assert science.filter.bandLabel == template.filter.bandLabel, \
795 "Science and template exposures have different bands: %s, %s" % \
796 (science.filter, template.filter)
797
798 assert templateBBox.contains(scienceBBox), \
799 "Template bbox does not contain all of the science image."
800
801 def _convolveExposure(self, exposure, kernel, convolutionControl,
802 bbox=None,
803 psf=None,
804 photoCalib=None,
805 interpolateBadMaskPlanes=False,
806 ):
807 """Convolve an exposure with the given kernel.
808
809 Parameters
810 ----------
811 exposure : `lsst.afw.Exposure`
812 exposure to convolve.
813 kernel : `lsst.afw.math.LinearCombinationKernel`
814 PSF matching kernel computed in the ``makeKernel`` subtask.
815 convolutionControl : `lsst.afw.math.ConvolutionControl`
816 Configuration for convolve algorithm.
817 bbox : `lsst.geom.Box2I`, optional
818 Bounding box to trim the convolved exposure to.
819 psf : `lsst.afw.detection.Psf`, optional
820 Point spread function (PSF) to set for the convolved exposure.
821 photoCalib : `lsst.afw.image.PhotoCalib`, optional
822 Photometric calibration of the convolved exposure.
823
824 Returns
825 -------
826 convolvedExp : `lsst.afw.Exposure`
827 The convolved image.
828 """
829 convolvedExposure = exposure.clone()
830 if psf is not None:
831 convolvedExposure.setPsf(psf)
832 if photoCalib is not None:
833 convolvedExposure.setPhotoCalib(photoCalib)
834 if interpolateBadMaskPlanes and self.config.badMaskPlanes is not None:
835 nInterp = _interpolateImage(convolvedExposure.maskedImage,
836 self.config.badMaskPlanes)
837 self.metadata["nInterpolated"] = nInterp
838 convolvedImage = lsst.afw.image.MaskedImageF(convolvedExposure.getBBox())
839 lsst.afw.math.convolve(convolvedImage, convolvedExposure.maskedImage, kernel, convolutionControl)
840 convolvedExposure.setMaskedImage(convolvedImage)
841 if bbox is None:
842 return convolvedExposure
843 else:
844 return convolvedExposure[bbox]
845
846 def _sourceSelector(self, sources, mask):
847 """Select sources from a catalog that meet the selection criteria.
848
849 Parameters
850 ----------
851 sources : `lsst.afw.table.SourceCatalog`
852 Input source catalog to select sources from.
853 mask : `lsst.afw.image.Mask`
854 The image mask plane to use to reject sources
855 based on their location on the ccd.
856
857 Returns
858 -------
859 selectSources : `lsst.afw.table.SourceCatalog`
860 The input source catalog, with flagged and low signal-to-noise
861 sources removed.
862
863 Raises
864 ------
865 RuntimeError
866 If there are too few sources to compute the PSF matching kernel
867 remaining after source selection.
868 """
869
870 selected = self.sourceSelector.selectSources(sources).selected
871 nInitialSelected = np.count_nonzero(selected)
872 selected *= self._checkMask(mask, sources, self.config.excludeMaskPlanes)
873 nSelected = np.count_nonzero(selected)
874 self.log.info("Rejecting %i candidate sources: an excluded template mask plane is set.",
875 nInitialSelected - nSelected)
876 selectSources = sources[selected].copy(deep=True)
877 # Trim selectSources if they exceed ``maxKernelSources``.
878 # Keep the highest signal-to-noise sources of those selected.
879 if (len(selectSources) > self.config.maxKernelSources) & (self.config.maxKernelSources > 0):
880 signalToNoise = selectSources.getPsfInstFlux()/selectSources.getPsfInstFluxErr()
881 indices = np.argsort(signalToNoise)
882 indices = indices[-self.config.maxKernelSources:]
883 selected = np.zeros(len(selectSources), dtype=bool)
884 selected[indices] = True
885 selectSources = selectSources[selected].copy(deep=True)
886
887 self.log.info("%i/%i=%.1f%% of sources selected for PSF matching from the input catalog",
888 len(selectSources), len(sources), 100*len(selectSources)/len(sources))
889 if len(selectSources) < self.config.minKernelSources:
890 self.log.error("Too few sources to calculate the PSF matching kernel: "
891 "%i selected but %i needed for the calculation.",
892 len(selectSources), self.config.minKernelSources)
893 if not self.config.allowKernelSourceDetection:
894 raise RuntimeError("Cannot compute PSF matching kernel: too few sources selected.")
895 self.metadata["nPsfSources"] = len(selectSources)
896
897 return selectSources
898
899 @staticmethod
900 def _checkMask(mask, sources, excludeMaskPlanes, checkAdjacent=True):
901 """Exclude sources that are located on or adjacent to masked pixels.
902
903 Parameters
904 ----------
905 mask : `lsst.afw.image.Mask`
906 The image mask plane to use to reject sources
907 based on the location of their centroid on the ccd.
908 sources : `lsst.afw.table.SourceCatalog`
909 The source catalog to evaluate.
910 excludeMaskPlanes : `list` of `str`
911 List of the names of the mask planes to exclude.
912
913 Returns
914 -------
915 flags : `numpy.ndarray` of `bool`
916 Array indicating whether each source in the catalog should be
917 kept (True) or rejected (False) based on the value of the
918 mask plane at its location.
919 """
920 setExcludeMaskPlanes = [
921 maskPlane for maskPlane in excludeMaskPlanes if maskPlane in mask.getMaskPlaneDict()
922 ]
923
924 excludePixelMask = mask.getPlaneBitMask(setExcludeMaskPlanes)
925
926 xv = (np.rint(sources.getX() - mask.getX0())).astype(int)
927 yv = (np.rint(sources.getY() - mask.getY0())).astype(int)
928
929 flags = np.ones(len(sources), dtype=bool)
930 if checkAdjacent:
931 pixRange = (0, -1, 1)
932 else:
933 pixRange = (0,)
934 for j in pixRange:
935 for i in pixRange:
936 mv = mask.array[yv + j, xv + i]
937 flags *= np.bitwise_and(mv, excludePixelMask) == 0
938 return flags
939
940 def _prepareInputs(self, template, science, visitSummary=None):
941 """Perform preparatory calculations common to all Alard&Lupton Tasks.
942
943 Parameters
944 ----------
945 template : `lsst.afw.image.ExposureF`
946 Template exposure, warped to match the science exposure. The
947 variance plane of the template image is modified in place.
948 science : `lsst.afw.image.ExposureF`
949 Science exposure to subtract from the template. The variance plane
950 of the science image is modified in place.
951 visitSummary : `lsst.afw.table.ExposureCatalog`, optional
952 Exposure catalog with external calibrations to be applied. Catalog
953 uses the detector id for the catalog id, sorted on id for fast
954 lookup.
955 """
956 self._validateExposures(template, science)
957 if visitSummary is not None:
958 self._applyExternalCalibrations(science, visitSummary=visitSummary)
959 templateCoverageFraction = checkTemplateIsSufficient(
960 template[science.getBBox()], self.log,
961 requiredTemplateFraction=self.config.requiredTemplateFraction,
962 exceptionMessage="Not attempting subtraction. To force subtraction,"
963 " set config requiredTemplateFraction=0"
964 )
965 self.metadata["templateCoveragePercent"] = 100*templateCoverageFraction
966
967 if self.config.doScaleVariance:
968 # Scale the variance of the template and science images before
969 # convolution, subtraction, or decorrelation so that they have the
970 # correct ratio.
971 templateVarFactor = self.scaleVariance.run(template.maskedImage)
972 sciVarFactor = self.scaleVariance.run(science.maskedImage)
973 self.log.info("Template variance scaling factor: %.2f", templateVarFactor)
974 self.metadata["scaleTemplateVarianceFactor"] = templateVarFactor
975 self.log.info("Science variance scaling factor: %.2f", sciVarFactor)
976 self.metadata["scaleScienceVarianceFactor"] = sciVarFactor
977
978 # Erase existing detection mask planes.
979 # We don't want the detection mask from the science image
980 self.updateMasks(template, science)
981
982 def updateMasks(self, template, science):
983 """Update the science and template mask planes before differencing.
984
985 Parameters
986 ----------
987 template : `lsst.afw.image.Exposure`
988 Template exposure, warped to match the science exposure.
989 The template mask planes will be erased, except for a few specified
990 in the task config.
991 science : `lsst.afw.image.Exposure`
992 Science exposure to subtract from the template.
993 The DETECTED and DETECTED_NEGATIVE mask planes of the science image
994 will be erased.
995 """
996 self._clearMask(science.mask, clearMaskPlanes=["DETECTED", "DETECTED_NEGATIVE"])
997
998 # We will clear ALL template mask planes, except for those specified
999 # via the `preserveTemplateMask` config. Mask planes specified via
1000 # the `renameTemplateMask` config will be copied to new planes with
1001 # "_TEMPLATE" appended to their names, and the original mask plane will
1002 # be cleared.
1003 clearMaskPlanes = [mp for mp in template.mask.getMaskPlaneDict().keys()
1004 if mp not in self.config.preserveTemplateMask]
1005 renameMaskPlanes = [mp for mp in self.config.renameTemplateMask
1006 if mp in template.mask.getMaskPlaneDict().keys()]
1007
1008 # propagate the mask plane related to Fake source injection
1009 # NOTE: the fake source injection sets FAKE plane, but it should be INJECTED
1010 # NOTE: This can be removed in DM-40796
1011 if "FAKE" in science.mask.getMaskPlaneDict().keys():
1012 self.log.info("Adding injected mask plane to science image")
1013 self._renameMaskPlanes(science.mask, "FAKE", "INJECTED")
1014 if "FAKE" in template.mask.getMaskPlaneDict().keys():
1015 self.log.info("Adding injected mask plane to template image")
1016 self._renameMaskPlanes(template.mask, "FAKE", "INJECTED_TEMPLATE")
1017 if "INJECTED" in renameMaskPlanes:
1018 renameMaskPlanes.remove("INJECTED")
1019 if "INJECTED_TEMPLATE" in clearMaskPlanes:
1020 clearMaskPlanes.remove("INJECTED_TEMPLATE")
1021
1022 for maskPlane in renameMaskPlanes:
1023 self._renameMaskPlanes(template.mask, maskPlane, maskPlane + "_TEMPLATE")
1024 self._clearMask(template.mask, clearMaskPlanes=clearMaskPlanes)
1025
1026 @staticmethod
1027 def _renameMaskPlanes(mask, maskPlane, newMaskPlane):
1028 """Rename a mask plane by adding the new name and copying the data.
1029
1030 Parameters
1031 ----------
1032 mask : `lsst.afw.image.Mask`
1033 The mask image to update in place.
1034 maskPlane : `str`
1035 The name of the existing mask plane to copy.
1036 newMaskPlane : `str`
1037 The new name of the mask plane that will be added.
1038 If the mask plane already exists, it will be updated in place.
1039 """
1040 mask.addMaskPlane(newMaskPlane)
1041 originBitMask = mask.getPlaneBitMask(maskPlane)
1042 destinationBitMask = mask.getPlaneBitMask(newMaskPlane)
1043 mask.array |= ((mask.array & originBitMask) > 0)*destinationBitMask
1044
1045 def _clearMask(self, mask, clearMaskPlanes=None):
1046 """Clear the mask plane of an exposure.
1047
1048 Parameters
1049 ----------
1050 mask : `lsst.afw.image.Mask`
1051 The mask plane to erase, which will be modified in place.
1052 clearMaskPlanes : `list` of `str`, optional
1053 Erase the specified mask planes.
1054 If not supplied, the entire mask will be erased.
1055 """
1056 if clearMaskPlanes is None:
1057 clearMaskPlanes = list(mask.getMaskPlaneDict().keys())
1058
1059 bitMaskToClear = mask.getPlaneBitMask(clearMaskPlanes)
1060 mask &= ~bitMaskToClear
1061
1062
1064 SubtractScoreOutputConnections):
1065 pass
1066
1067
1069 pipelineConnections=AlardLuptonPreconvolveSubtractConnections):
1070 pass
1071
1072
1074 """Subtract a template from a science image, convolving the science image
1075 before computing the kernel, and also convolving the template before
1076 subtraction.
1077 """
1078 ConfigClass = AlardLuptonPreconvolveSubtractConfig
1079 _DefaultName = "alardLuptonPreconvolveSubtract"
1080
1081 def run(self, template, science, sources, visitSummary=None):
1082 """Preconvolve the science image with its own PSF,
1083 convolve the template image with a PSF-matching kernel and subtract
1084 from the preconvolved science image.
1085
1086 Parameters
1087 ----------
1088 template : `lsst.afw.image.ExposureF`
1089 The template image, which has previously been warped to the science
1090 image. The template bbox will be padded by a few pixels compared to
1091 the science bbox.
1092 science : `lsst.afw.image.ExposureF`
1093 The science exposure.
1094 sources : `lsst.afw.table.SourceCatalog`
1095 Identified sources on the science exposure. This catalog is used to
1096 select sources in order to perform the AL PSF matching on stamp
1097 images around them.
1098 visitSummary : `lsst.afw.table.ExposureCatalog`, optional
1099 Exposure catalog with complete external calibrations. Catalog uses
1100 the detector id for the catalog id, sorted on id for fast lookup.
1101
1102 Returns
1103 -------
1104 results : `lsst.pipe.base.Struct`
1105 ``scoreExposure`` : `lsst.afw.image.ExposureF`
1106 Result of subtracting the convolved template and science
1107 images. Attached PSF is that of the original science image.
1108 ``matchedTemplate`` : `lsst.afw.image.ExposureF`
1109 Warped and PSF-matched template exposure. Attached PSF is that
1110 of the original science image.
1111 ``matchedScience`` : `lsst.afw.image.ExposureF`
1112 The science exposure after convolving with its own PSF.
1113 Attached PSF is that of the original science image.
1114 ``backgroundModel`` : `lsst.afw.math.Function2D`
1115 Background model that was fit while solving for the
1116 PSF-matching kernel
1117 ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
1118 Final kernel used to PSF-match the template to the science
1119 image.
1120 """
1121 self._prepareInputs(template, science, visitSummary=visitSummary)
1122
1123 # TODO: DM-37212 we need to mirror the kernel in order to get correct cross correlation
1124 scienceKernel = science.psf.getKernel()
1125 matchedScience = self._convolveExposure(science, scienceKernel, self.convolutionControl,
1126 interpolateBadMaskPlanes=True)
1127 self.metadata["convolvedExposure"] = "Preconvolution"
1128 try:
1129 selectSources = self._sourceSelector(sources, matchedScience.mask)
1130 subtractResults = self.runPreconvolve(template, science, matchedScience,
1131 selectSources, scienceKernel)
1132
1133 except (RuntimeError, lsst.pex.exceptions.Exception) as e:
1134 self.log.warning("Failed to match template. Checking coverage")
1135 # Raise NoWorkFound if template fraction is insufficient
1136 checkTemplateIsSufficient(template[science.getBBox()], self.log,
1137 self.config.minTemplateFractionForExpectedSuccess,
1138 exceptionMessage="Template coverage lower than expected to succeed."
1139 f" Failure is tolerable: {e}")
1140 # checkTemplateIsSufficient did not raise NoWorkFound, so raise original exception
1141 raise e
1142
1143 return subtractResults
1144
1145 def runPreconvolve(self, template, science, matchedScience, selectSources, preConvKernel):
1146 """Convolve the science image with its own PSF, then convolve the
1147 template with a matching kernel and subtract to form the Score
1148 exposure.
1149
1150 Parameters
1151 ----------
1152 template : `lsst.afw.image.ExposureF`
1153 Template exposure, warped to match the science exposure.
1154 science : `lsst.afw.image.ExposureF`
1155 Science exposure to subtract from the template.
1156 matchedScience : `lsst.afw.image.ExposureF`
1157 The science exposure, convolved with the reflection of its own PSF.
1158 selectSources : `lsst.afw.table.SourceCatalog`
1159 Identified sources on the science exposure. This catalog is used to
1160 select sources in order to perform the AL PSF matching on stamp
1161 images around them.
1162 preConvKernel : `lsst.afw.math.Kernel`
1163 The reflection of the kernel that was used to preconvolve the
1164 `science` exposure. Must be normalized to sum to 1.
1165
1166 Returns
1167 -------
1168 results : `lsst.pipe.base.Struct`
1169
1170 ``scoreExposure`` : `lsst.afw.image.ExposureF`
1171 Result of subtracting the convolved template and science
1172 images. Attached PSF is that of the original science image.
1173 ``matchedTemplate`` : `lsst.afw.image.ExposureF`
1174 Warped and PSF-matched template exposure. Attached PSF is that
1175 of the original science image.
1176 ``matchedScience`` : `lsst.afw.image.ExposureF`
1177 The science exposure after convolving with its own PSF.
1178 Attached PSF is that of the original science image.
1179 ``backgroundModel`` : `lsst.afw.math.Function2D`
1180 Background model that was fit while solving for the
1181 PSF-matching kernel
1182 ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
1183 Final kernel used to PSF-match the template to the science
1184 image.
1185 """
1186 bbox = science.getBBox()
1187 innerBBox = preConvKernel.shrinkBBox(bbox)
1188
1189 kernelSources = self.makeKernel.selectKernelSources(template[innerBBox], matchedScience[innerBBox],
1190 candidateList=selectSources,
1191 preconvolved=True)
1192 kernelResult = self.makeKernel.run(template[innerBBox], matchedScience[innerBBox], kernelSources,
1193 preconvolved=True)
1194
1195 matchedTemplate = self._convolveExposure(template, kernelResult.psfMatchingKernel,
1196 self.convolutionControl,
1197 bbox=bbox,
1198 psf=science.psf,
1199 interpolateBadMaskPlanes=True,
1200 photoCalib=science.photoCalib)
1201 score = _subtractImages(matchedScience, matchedTemplate,
1202 backgroundModel=(kernelResult.backgroundModel
1203 if self.config.doSubtractBackground else None))
1204 correctedScore = self.finalize(template[bbox], science, score,
1205 kernelResult.psfMatchingKernel,
1206 templateMatched=True, preConvMode=True,
1207 preConvKernel=preConvKernel)
1208
1209 return lsst.pipe.base.Struct(scoreExposure=correctedScore,
1210 matchedTemplate=matchedTemplate,
1211 matchedScience=matchedScience,
1212 backgroundModel=kernelResult.backgroundModel,
1213 psfMatchingKernel=kernelResult.psfMatchingKernel,
1214 kernelSources=kernelSources)
1215
1216
1217def checkTemplateIsSufficient(templateExposure, logger, requiredTemplateFraction=0.,
1218 exceptionMessage=""):
1219 """Raise NoWorkFound if template coverage < requiredTemplateFraction
1220
1221 Parameters
1222 ----------
1223 templateExposure : `lsst.afw.image.ExposureF`
1224 The template exposure to check
1225 logger : `logging.Logger`
1226 Logger for printing output.
1227 requiredTemplateFraction : `float`, optional
1228 Fraction of pixels of the science image required to have coverage
1229 in the template.
1230 exceptionMessage : `str`, optional
1231 Message to include in the exception raised if the template coverage
1232 is insufficient.
1233
1234 Returns
1235 -------
1236 templateCoverageFraction: `float`
1237 Fraction of pixels in the template with data.
1238
1239 Raises
1240 ------
1241 lsst.pipe.base.NoWorkFound
1242 Raised if fraction of good pixels, defined as not having NO_DATA
1243 set, is less than the requiredTemplateFraction
1244 """
1245 # Count the number of pixels with the NO_DATA mask bit set
1246 # counting NaN pixels is insufficient because pixels without data are often intepolated over)
1247 pixNoData = np.count_nonzero(templateExposure.mask.array
1248 & templateExposure.mask.getPlaneBitMask('NO_DATA'))
1249 pixGood = templateExposure.getBBox().getArea() - pixNoData
1250 templateCoverageFraction = pixGood/templateExposure.getBBox().getArea()
1251 logger.info("template has %d good pixels (%.1f%%)", pixGood, 100*templateCoverageFraction)
1252
1253 if templateCoverageFraction < requiredTemplateFraction:
1254 message = ("Insufficient Template Coverage. (%.1f%% < %.1f%%)" % (
1255 100*templateCoverageFraction,
1256 100*requiredTemplateFraction))
1257 raise lsst.pipe.base.NoWorkFound(message + " " + exceptionMessage)
1258 return templateCoverageFraction
1259
1260
1261def _subtractImages(science, template, backgroundModel=None):
1262 """Subtract template from science, propagating relevant metadata.
1263
1264 Parameters
1265 ----------
1266 science : `lsst.afw.Exposure`
1267 The input science image.
1268 template : `lsst.afw.Exposure`
1269 The template to subtract from the science image.
1270 backgroundModel : `lsst.afw.MaskedImage`, optional
1271 Differential background model
1272
1273 Returns
1274 -------
1275 difference : `lsst.afw.Exposure`
1276 The subtracted image.
1277 """
1278 difference = science.clone()
1279 if backgroundModel is not None:
1280 difference.maskedImage -= backgroundModel
1281 difference.maskedImage -= template.maskedImage
1282 return difference
1283
1284
1285def _shapeTest(exp1, exp2, fwhmExposureBuffer, fwhmExposureGrid):
1286 """Determine that the PSF of ``exp1`` is not wider than that of ``exp2``.
1287
1288 Parameters
1289 ----------
1290 exp1 : `~lsst.afw.image.Exposure`
1291 Exposure with the reference point spread function (PSF) to evaluate.
1292 exp2 : `~lsst.afw.image.Exposure`
1293 Exposure with a candidate point spread function (PSF) to evaluate.
1294 fwhmExposureBuffer : `float`
1295 Fractional buffer margin to be left out of all sides of the image
1296 during the construction of the grid to compute mean PSF FWHM in an
1297 exposure, if the PSF is not available at its average position.
1298 fwhmExposureGrid : `int`
1299 Grid size to compute the mean FWHM in an exposure, if the PSF is not
1300 available at its average position.
1301 Returns
1302 -------
1303 result : `bool`
1304 True if ``exp1`` has a PSF that is not wider than that of ``exp2`` in
1305 either dimension.
1306 """
1307 try:
1308 shape1 = getPsfFwhm(exp1.psf, average=False)
1309 shape2 = getPsfFwhm(exp2.psf, average=False)
1311 shape1 = evaluateMeanPsfFwhm(exp1,
1312 fwhmExposureBuffer=fwhmExposureBuffer,
1313 fwhmExposureGrid=fwhmExposureGrid
1314 )
1315 shape2 = evaluateMeanPsfFwhm(exp2,
1316 fwhmExposureBuffer=fwhmExposureBuffer,
1317 fwhmExposureGrid=fwhmExposureGrid
1318 )
1319 return shape1 <= shape2
1320
1321 # Results from getPsfFwhm is a tuple of two values, one for each dimension.
1322 xTest = shape1[0] <= shape2[0]
1323 yTest = shape1[1] <= shape2[1]
1324 return xTest | yTest
1325
1326
1327def _interpolateImage(maskedImage, badMaskPlanes, fallbackValue=None):
1328 """Replace masked image pixels with interpolated values.
1329
1330 Parameters
1331 ----------
1332 maskedImage : `lsst.afw.image.MaskedImage`
1333 Image on which to perform interpolation.
1334 badMaskPlanes : `list` of `str`
1335 List of mask planes to interpolate over.
1336 fallbackValue : `float`, optional
1337 Value to set when interpolation fails.
1338
1339 Returns
1340 -------
1341 result: `float`
1342 The number of masked pixels that were replaced.
1343 """
1344 imgBadMaskPlanes = [
1345 maskPlane for maskPlane in badMaskPlanes if maskPlane in maskedImage.mask.getMaskPlaneDict()
1346 ]
1347
1348 image = maskedImage.image.array
1349 badPixels = (maskedImage.mask.array & maskedImage.mask.getPlaneBitMask(imgBadMaskPlanes)) > 0
1350 image[badPixels] = np.nan
1351 if fallbackValue is None:
1352 fallbackValue = np.nanmedian(image)
1353 # For this initial implementation, skip the interpolation and just fill with
1354 # the median value.
1355 image[badPixels] = fallbackValue
1356 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)