24from astropy
import units
as u
31from lsst.utils.introspection
import find_outside_stacklevel
32from lsst.ip.diffim.utils
import evaluateMeanPsfFwhm, getPsfFwhm
38from .
import MakeKernelTask, DecorrelateALKernelTask
39from lsst.utils.timer
import timeMethod
41__all__ = [
"AlardLuptonSubtractConfig",
"AlardLuptonSubtractTask",
42 "AlardLuptonPreconvolveSubtractConfig",
"AlardLuptonPreconvolveSubtractTask"]
44_dimensions = (
"instrument",
"visit",
"detector")
45_defaultTemplates = {
"coaddName":
"deep",
"fakesType":
""}
49 dimensions=_dimensions,
50 defaultTemplates=_defaultTemplates):
51 template = connectionTypes.Input(
52 doc=
"Input warped template to subtract.",
53 dimensions=(
"instrument",
"visit",
"detector"),
54 storageClass=
"ExposureF",
55 name=
"{fakesType}{coaddName}Diff_templateExp"
57 science = connectionTypes.Input(
58 doc=
"Input science exposure to subtract from.",
59 dimensions=(
"instrument",
"visit",
"detector"),
60 storageClass=
"ExposureF",
61 name=
"{fakesType}calexp"
63 sources = connectionTypes.Input(
64 doc=
"Sources measured on the science exposure; "
65 "used to select sources for making the matching kernel.",
66 dimensions=(
"instrument",
"visit",
"detector"),
67 storageClass=
"SourceCatalog",
70 visitSummary = connectionTypes.Input(
71 doc=(
"Per-visit catalog with final calibration objects. "
72 "These catalogs use the detector id for the catalog id, "
73 "sorted on id for fast lookup."),
74 dimensions=(
"instrument",
"visit"),
75 storageClass=
"ExposureCatalog",
76 name=
"finalVisitSummary",
81 if not config.doApplyExternalCalibrations:
86 dimensions=_dimensions,
87 defaultTemplates=_defaultTemplates):
88 difference = connectionTypes.Output(
89 doc=
"Result of subtracting convolved template from science image.",
90 dimensions=(
"instrument",
"visit",
"detector"),
91 storageClass=
"ExposureF",
92 name=
"{fakesType}{coaddName}Diff_differenceTempExp",
94 matchedTemplate = connectionTypes.Output(
95 doc=
"Warped and PSF-matched template used to create `subtractedExposure`.",
96 dimensions=(
"instrument",
"visit",
"detector"),
97 storageClass=
"ExposureF",
98 name=
"{fakesType}{coaddName}Diff_matchedExp",
100 psfMatchingKernel = connectionTypes.Output(
101 doc=
"Kernel used to PSF match the science and template images.",
102 dimensions=(
"instrument",
"visit",
"detector"),
103 storageClass=
"MatchingKernel",
104 name=
"{fakesType}{coaddName}Diff_psfMatchKernel",
109 dimensions=_dimensions,
110 defaultTemplates=_defaultTemplates):
111 scoreExposure = connectionTypes.Output(
112 doc=
"The maximum likelihood image, used for the detection of diaSources.",
113 dimensions=(
"instrument",
"visit",
"detector"),
114 storageClass=
"ExposureF",
115 name=
"{fakesType}{coaddName}Diff_scoreExp",
117 psfMatchingKernel = connectionTypes.Output(
118 doc=
"Kernel used to PSF match the science and template images.",
119 dimensions=(
"instrument",
"visit",
"detector"),
120 storageClass=
"MatchingKernel",
121 name=
"{fakesType}{coaddName}Diff_psfScoreMatchKernel",
131 target=MakeKernelTask,
132 doc=
"Task to construct a matching kernel for convolution.",
137 doc=
"Perform diffim decorrelation to undo pixel correlation due to A&L "
138 "kernel convolution? If True, also update the diffim PSF."
141 target=DecorrelateALKernelTask,
142 doc=
"Task to decorrelate the image difference.",
147 doc=
"Raise NoWorkFound and do not attempt image subtraction if template covers less than this "
148 " fraction of pixels. Setting to 0 will always attempt image subtraction."
153 doc=
"Raise NoWorkFound if PSF-matching fails and template covers less than this fraction of pixels."
154 " If the fraction of pixels covered by the template is less than this value (and greater than"
155 " requiredTemplateFraction) this task is attempted but failure is anticipated and tolerated."
160 doc=
"Scale variance of the image difference?"
163 target=ScaleVarianceTask,
164 doc=
"Subtask to rescale the variance of the template to the statistically expected level."
167 doc=
"Subtract the background fit when solving the kernel?",
173 "Replace science Exposure's calibration objects with those"
174 " in visitSummary. Ignored if `doApplyFinalizedPsf is True."
182 doc=
"Minimum signal to noise ratio of detected sources "
183 "to use for calculating the PSF matching kernel."
188 doc=
"Maximum signal to noise ratio of detected sources "
189 "to use for calculating the PSF matching kernel."
194 doc=
"Maximum number of sources to use for calculating the PSF matching kernel."
195 "Set to -1 to disable."
200 doc=
"Minimum number of sources needed for calculating the PSF matching kernel."
204 doc=
"Flags that, if set, the associated source should not "
205 "be used to determine the PSF matching kernel.",
206 default=(
"sky_source",
"slot_Centroid_flag",
207 "slot_ApFlux_flag",
"slot_PsfFlux_flag",
208 "base_PixelFlags_flag_interpolated",
209 "base_PixelFlags_flag_saturated",
210 "base_PixelFlags_flag_bad",
215 default=(
"NO_DATA",
"BAD",
"SAT",
"EDGE",
"FAKE"),
216 doc=
"Mask planes to exclude when selecting sources for PSF matching."
220 default=(
"NO_DATA",
"BAD",
"SAT",
"EDGE"),
221 doc=
"Mask planes to interpolate over."
225 default=(
"NO_DATA",
"BAD",),
226 doc=
"Mask planes from the template to propagate to the image difference."
230 default=(
"SAT",
"INJECTED",
"INJECTED_CORE",),
231 doc=
"Mask planes from the template to propagate to the image difference"
232 "with '_TEMPLATE' appended to the name."
237 doc=
"Re-run source detection for kernel candidates if an error is"
238 " encountered while calculating the matching kernel."
244 self.
makeKernel.kernel.active.spatialKernelOrder = 1
245 self.
makeKernel.kernel.active.spatialBgOrder = 2
249 pipelineConnections=AlardLuptonSubtractConnections):
252 default=
"convolveTemplate",
253 allowed={
"auto":
"Choose which image to convolve at runtime.",
254 "convolveScience":
"Only convolve the science image.",
255 "convolveTemplate":
"Only convolve the template image."},
256 doc=
"Choose which image to convolve at runtime, or require that a specific image is convolved."
261 """Compute the image difference of a science and template image using
262 the Alard & Lupton (1998) algorithm.
264 ConfigClass = AlardLuptonSubtractConfig
265 _DefaultName =
"alardLuptonSubtract"
269 self.makeSubtask(
"decorrelate")
270 self.makeSubtask(
"makeKernel")
271 if self.config.doScaleVariance:
272 self.makeSubtask(
"scaleVariance")
281 """Replace calibrations (psf, and ApCorrMap) on this exposure with
286 exposure : `lsst.afw.image.exposure.Exposure`
287 Input exposure to adjust calibrations.
288 visitSummary : `lsst.afw.table.ExposureCatalog`
289 Exposure catalog with external calibrations to be applied. Catalog
290 uses the detector id for the catalog id, sorted on id for fast
295 exposure : `lsst.afw.image.exposure.Exposure`
296 Exposure with adjusted calibrations.
298 detectorId = exposure.info.getDetector().getId()
300 row = visitSummary.find(detectorId)
302 self.
log.warning(
"Detector id %s not found in external calibrations catalog; "
303 "Using original calibrations.", detectorId)
306 apCorrMap = row.getApCorrMap()
308 self.
log.warning(
"Detector id %s has None for psf in "
309 "external calibrations catalog; Using original psf and aperture correction.",
311 elif apCorrMap
is None:
312 self.
log.warning(
"Detector id %s has None for apCorrMap in "
313 "external calibrations catalog; Using original psf and aperture correction.",
317 exposure.info.setApCorrMap(apCorrMap)
322 def run(self, template, science, sources, finalizedPsfApCorrCatalog=None,
324 """PSF match, subtract, and decorrelate two images.
328 template : `lsst.afw.image.ExposureF`
329 Template exposure, warped to match the science exposure.
330 science : `lsst.afw.image.ExposureF`
331 Science exposure to subtract from the template.
332 sources : `lsst.afw.table.SourceCatalog`
333 Identified sources on the science exposure. This catalog is used to
334 select sources in order to perform the AL PSF matching on stamp
336 finalizedPsfApCorrCatalog : `lsst.afw.table.ExposureCatalog`, optional
337 Exposure catalog with finalized psf models and aperture correction
338 maps to be applied. Catalog uses the detector id for the catalog
339 id, sorted on id for fast lookup. Deprecated in favor of
340 ``visitSummary``, and will be removed after v26.
341 visitSummary : `lsst.afw.table.ExposureCatalog`, optional
342 Exposure catalog with external calibrations to be applied. Catalog
343 uses the detector id for the catalog id, sorted on id for fast
344 lookup. Ignored (for temporary backwards compatibility) if
345 ``finalizedPsfApCorrCatalog`` is provided.
349 results : `lsst.pipe.base.Struct`
350 ``difference`` : `lsst.afw.image.ExposureF`
351 Result of subtracting template and science.
352 ``matchedTemplate`` : `lsst.afw.image.ExposureF`
353 Warped and PSF-matched template exposure.
354 ``backgroundModel`` : `lsst.afw.math.Function2D`
355 Background model that was fit while solving for the
357 ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
358 Kernel used to PSF-match the convolved image.
363 If an unsupported convolution mode is supplied.
365 If there are too few sources to calculate the PSF matching kernel.
366 lsst.pipe.base.NoWorkFound
367 Raised if fraction of good pixels, defined as not having NO_DATA
368 set, is less then the configured requiredTemplateFraction
371 if finalizedPsfApCorrCatalog
is not None:
373 "The finalizedPsfApCorrCatalog argument is deprecated in favor of the visitSummary "
374 "argument, and will be removed after v26.",
376 stacklevel=find_outside_stacklevel(
"lsst.ip.diffim"),
378 visitSummary = finalizedPsfApCorrCatalog
383 fwhmExposureBuffer = self.config.makeKernel.fwhmExposureBuffer
384 fwhmExposureGrid = self.config.makeKernel.fwhmExposureGrid
398 templatePsfSize = getPsfFwhm(template.psf)
399 sciencePsfSize = getPsfFwhm(science.psf)
401 self.
log.info(
"Unable to evaluate PSF at the average position. "
402 "Evaluting PSF on a grid of points."
404 templatePsfSize = evaluateMeanPsfFwhm(template,
405 fwhmExposureBuffer=fwhmExposureBuffer,
406 fwhmExposureGrid=fwhmExposureGrid
408 sciencePsfSize = evaluateMeanPsfFwhm(science,
409 fwhmExposureBuffer=fwhmExposureBuffer,
410 fwhmExposureGrid=fwhmExposureGrid
412 self.
log.info(
"Science PSF FWHM: %f pixels", sciencePsfSize)
413 self.
log.info(
"Template PSF FWHM: %f pixels", templatePsfSize)
414 self.metadata.add(
"sciencePsfSize", sciencePsfSize)
415 self.metadata.add(
"templatePsfSize", templatePsfSize)
418 maglim_science = self.
_calculateMagLim(science, fallbackPsfSize=sciencePsfSize)
419 fluxlim_science = (maglim_science*u.ABmag).to_value(u.nJy)
420 maglim_template = self.
_calculateMagLim(template, fallbackPsfSize=templatePsfSize)
421 if np.isnan(maglim_template):
422 self.
log.info(
"Cannot evaluate template limiting mag; adopting science limiting mag for diffim")
423 maglim_diffim = maglim_science
425 fluxlim_template = (maglim_template*u.ABmag).to_value(u.nJy)
426 maglim_diffim = (np.sqrt(fluxlim_science**2 + fluxlim_template**2)*u.nJy).to(u.ABmag).value
427 self.metadata.add(
"scienceLimitingMagnitude", maglim_science)
428 self.metadata.add(
"templateLimitingMagnitude", maglim_template)
429 self.metadata.add(
"diffimLimitingMagnitude", maglim_diffim)
431 if self.config.mode ==
"auto":
434 fwhmExposureBuffer=fwhmExposureBuffer,
435 fwhmExposureGrid=fwhmExposureGrid)
437 if sciencePsfSize < templatePsfSize:
438 self.
log.info(
"Average template PSF size is greater, "
439 "but science PSF greater in one dimension: convolving template image.")
441 self.
log.info(
"Science PSF size is greater: convolving template image.")
443 self.
log.info(
"Template PSF size is greater: convolving science image.")
444 elif self.config.mode ==
"convolveTemplate":
445 self.
log.info(
"`convolveTemplate` is set: convolving template image.")
446 convolveTemplate =
True
447 elif self.config.mode ==
"convolveScience":
448 self.
log.info(
"`convolveScience` is set: convolving science image.")
449 convolveTemplate =
False
451 raise RuntimeError(
"Cannot handle AlardLuptonSubtract mode: %s", self.config.mode)
454 sourceMask = science.mask.clone()
455 sourceMask.array |= template[science.getBBox()].mask.array
458 self.metadata.add(
"convolvedExposure",
"Template")
461 self.metadata.add(
"convolvedExposure",
"Science")
465 self.
log.warning(
"Failed to match template. Checking coverage")
468 self.config.minTemplateFractionForExpectedSuccess,
469 exceptionMessage=
"Template coverage lower than expected to succeed."
470 f
" Failure is tolerable: {e}")
474 return subtractResults
477 """Convolve the template image with a PSF-matching kernel and subtract
478 from the science image.
482 template : `lsst.afw.image.ExposureF`
483 Template exposure, warped to match the science exposure.
484 science : `lsst.afw.image.ExposureF`
485 Science exposure to subtract from the template.
486 selectSources : `lsst.afw.table.SourceCatalog`
487 Identified sources on the science exposure. This catalog is used to
488 select sources in order to perform the AL PSF matching on stamp
493 results : `lsst.pipe.base.Struct`
495 ``difference`` : `lsst.afw.image.ExposureF`
496 Result of subtracting template and science.
497 ``matchedTemplate`` : `lsst.afw.image.ExposureF`
498 Warped and PSF-matched template exposure.
499 ``backgroundModel`` : `lsst.afw.math.Function2D`
500 Background model that was fit while solving for the PSF-matching kernel
501 ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
502 Kernel used to PSF-match the template to the science image.
505 kernelSources = self.makeKernel.selectKernelSources(template, science,
506 candidateList=selectSources,
508 kernelResult = self.makeKernel.run(template, science, kernelSources,
510 except Exception
as e:
511 if self.config.allowKernelSourceDetection:
512 self.
log.warning(
"Error encountered trying to construct the matching kernel"
513 f
" Running source detection and retrying. {e}")
514 kernelSources = self.makeKernel.selectKernelSources(template, science,
517 kernelResult = self.makeKernel.run(template, science, kernelSources,
522 matchedTemplate = self.
_convolveExposure(template, kernelResult.psfMatchingKernel,
524 bbox=science.getBBox(),
526 photoCalib=science.photoCalib)
529 backgroundModel=(kernelResult.backgroundModel
530 if self.config.doSubtractBackground
else None))
531 correctedExposure = self.
finalize(template, science, difference,
532 kernelResult.psfMatchingKernel,
533 templateMatched=
True)
535 return lsst.pipe.base.Struct(difference=correctedExposure,
536 matchedTemplate=matchedTemplate,
537 matchedScience=science,
538 backgroundModel=kernelResult.backgroundModel,
539 psfMatchingKernel=kernelResult.psfMatchingKernel)
542 """Convolve the science image with a PSF-matching kernel and subtract
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
558 results : `lsst.pipe.base.Struct`
560 ``difference`` : `lsst.afw.image.ExposureF`
561 Result of subtracting template and science.
562 ``matchedTemplate`` : `lsst.afw.image.ExposureF`
563 Warped template exposure. Note that in this case, the template
564 is not PSF-matched to the science image.
565 ``backgroundModel`` : `lsst.afw.math.Function2D`
566 Background model that was fit while solving for the PSF-matching kernel
567 ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
568 Kernel used to PSF-match the science image to the template.
570 bbox = science.getBBox()
571 kernelSources = self.makeKernel.selectKernelSources(science, template,
572 candidateList=selectSources,
574 kernelResult = self.makeKernel.run(science, template, kernelSources,
576 modelParams = kernelResult.backgroundModel.getParameters()
578 kernelResult.backgroundModel.setParameters([-p
for p
in modelParams])
580 kernelImage = lsst.afw.image.ImageD(kernelResult.psfMatchingKernel.getDimensions())
581 norm = kernelResult.psfMatchingKernel.computeImage(kernelImage, doNormalize=
False)
588 matchedScience.maskedImage /= norm
589 matchedTemplate = template.clone()[bbox]
590 matchedTemplate.maskedImage /= norm
591 matchedTemplate.setPhotoCalib(science.photoCalib)
594 backgroundModel=(kernelResult.backgroundModel
595 if self.config.doSubtractBackground
else None))
597 correctedExposure = self.
finalize(template, science, difference,
598 kernelResult.psfMatchingKernel,
599 templateMatched=
False)
601 return lsst.pipe.base.Struct(difference=correctedExposure,
602 matchedTemplate=matchedTemplate,
603 matchedScience=matchedScience,
604 backgroundModel=kernelResult.backgroundModel,
605 psfMatchingKernel=kernelResult.psfMatchingKernel,)
607 def finalize(self, template, science, difference, kernel,
608 templateMatched=True,
611 spatiallyVarying=False):
612 """Decorrelate the difference image to undo the noise correlations
613 caused by convolution.
617 template : `lsst.afw.image.ExposureF`
618 Template exposure, warped to match the science exposure.
619 science : `lsst.afw.image.ExposureF`
620 Science exposure to subtract from the template.
621 difference : `lsst.afw.image.ExposureF`
622 Result of subtracting template and science.
623 kernel : `lsst.afw.math.Kernel`
624 An (optionally spatially-varying) PSF matching kernel
625 templateMatched : `bool`, optional
626 Was the template PSF-matched to the science image?
627 preConvMode : `bool`, optional
628 Was the science image preconvolved with its own PSF
629 before PSF matching the template?
630 preConvKernel : `lsst.afw.detection.Psf`, optional
631 If not `None`, then the science image was pre-convolved with
632 (the reflection of) this kernel. Must be normalized to sum to 1.
633 spatiallyVarying : `bool`, optional
634 Compute the decorrelation kernel spatially varying across the image?
638 correctedExposure : `lsst.afw.image.ExposureF`
639 The decorrelated image difference.
641 if self.config.doDecorrelation:
642 self.
log.info(
"Decorrelating image difference.")
646 correctedExposure = self.decorrelate.run(science, template[science.getBBox()], difference, kernel,
647 templateMatched=templateMatched,
648 preConvMode=preConvMode,
649 preConvKernel=preConvKernel,
650 spatiallyVarying=spatiallyVarying).correctedExposure
652 self.
log.info(
"NOT decorrelating image difference.")
653 correctedExposure = difference
654 return correctedExposure
657 """Calculate an exposure's limiting magnitude.
659 This method uses the photometric zeropoint together with the
660 PSF size from the average position of the exposure.
664 exposure : `lsst.afw.image.Exposure`
665 The target exposure to calculate the limiting magnitude for.
666 nsigma : `float`, optional
667 The detection threshold in sigma.
668 fallbackPsfSize : `float`, optional
669 PSF FWHM to use in the event the exposure PSF cannot be retrieved.
673 maglim : `astropy.units.Quantity`
674 The limiting magnitude of the exposure, or np.nan.
677 psf = exposure.getPsf()
678 psf_shape = psf.computeShape(psf.getAveragePosition())
680 if fallbackPsfSize
is not None:
681 self.
log.info(
"Unable to evaluate PSF, using fallback FWHM %f", fallbackPsfSize)
682 psf_area = np.pi*(fallbackPsfSize/2)**2
683 zeropoint = exposure.getPhotoCalib().instFluxToMagnitude(1)
684 maglim = zeropoint - 2.5*np.log10(nsigma*np.sqrt(psf_area))
686 self.
log.info(
"Unable to evaluate PSF, setting maglim to nan")
690 psf_area = np.pi*np.sqrt(psf_shape.getIxx()*psf_shape.getIyy())
691 zeropoint = exposure.getPhotoCalib().instFluxToMagnitude(1)
692 maglim = zeropoint - 2.5*np.log10(nsigma*np.sqrt(psf_area))
698 """Check that the WCS of the two Exposures match, and the template bbox
699 contains the science bbox.
703 template : `lsst.afw.image.ExposureF`
704 Template exposure, warped to match the science exposure.
705 science : `lsst.afw.image.ExposureF`
706 Science exposure to subtract from the template.
711 Raised if the WCS of the template is not equal to the science WCS,
712 or if the science image is not fully contained in the template
715 assert template.wcs == science.wcs, \
716 "Template and science exposure WCS are not identical."
717 templateBBox = template.getBBox()
718 scienceBBox = science.getBBox()
720 assert templateBBox.contains(scienceBBox), \
721 "Template bbox does not contain all of the science image."
727 interpolateBadMaskPlanes=False,
729 """Convolve an exposure with the given kernel.
733 exposure : `lsst.afw.Exposure`
734 exposure to convolve.
735 kernel : `lsst.afw.math.LinearCombinationKernel`
736 PSF matching kernel computed in the ``makeKernel`` subtask.
737 convolutionControl : `lsst.afw.math.ConvolutionControl`
738 Configuration for convolve algorithm.
739 bbox : `lsst.geom.Box2I`, optional
740 Bounding box to trim the convolved exposure to.
741 psf : `lsst.afw.detection.Psf`, optional
742 Point spread function (PSF) to set for the convolved exposure.
743 photoCalib : `lsst.afw.image.PhotoCalib`, optional
744 Photometric calibration of the convolved exposure.
748 convolvedExp : `lsst.afw.Exposure`
751 convolvedExposure = exposure.clone()
753 convolvedExposure.setPsf(psf)
754 if photoCalib
is not None:
755 convolvedExposure.setPhotoCalib(photoCalib)
756 if interpolateBadMaskPlanes
and self.config.badMaskPlanes
is not None:
758 self.config.badMaskPlanes)
759 self.metadata.add(
"nInterpolated", nInterp)
760 convolvedImage = lsst.afw.image.MaskedImageF(convolvedExposure.getBBox())
762 convolvedExposure.setMaskedImage(convolvedImage)
764 return convolvedExposure
766 return convolvedExposure[bbox]
769 """Select sources from a catalog that meet the selection criteria.
773 sources : `lsst.afw.table.SourceCatalog`
774 Input source catalog to select sources from.
775 mask : `lsst.afw.image.Mask`
776 The image mask plane to use to reject sources
777 based on their location on the ccd.
781 selectSources : `lsst.afw.table.SourceCatalog`
782 The input source catalog, with flagged and low signal-to-noise
788 If there are too few sources to compute the PSF matching kernel
789 remaining after source selection.
791 flags = np.ones(len(sources), dtype=bool)
792 for flag
in self.config.badSourceFlags:
794 flags *= ~sources[flag]
795 except Exception
as e:
796 self.
log.warning(
"Could not apply source flag: %s", e)
797 signalToNoise = sources.getPsfInstFlux()/sources.getPsfInstFluxErr()
798 sToNFlag = signalToNoise > self.config.detectionThreshold
800 sToNFlagMax = signalToNoise < self.config.detectionThresholdMax
802 flags *= self.
_checkMask(mask, sources, self.config.excludeMaskPlanes)
803 selectSources = sources[flags].copy(deep=
True)
804 if (len(selectSources) > self.config.maxKernelSources) & (self.config.maxKernelSources > 0):
805 signalToNoise = selectSources.getPsfInstFlux()/selectSources.getPsfInstFluxErr()
806 indices = np.argsort(signalToNoise)
807 indices = indices[-self.config.maxKernelSources:]
808 flags = np.zeros(len(selectSources), dtype=bool)
809 flags[indices] =
True
810 selectSources = selectSources[flags].copy(deep=
True)
812 self.
log.info(
"%i/%i=%.1f%% of sources selected for PSF matching from the input catalog",
813 len(selectSources), len(sources), 100*len(selectSources)/len(sources))
814 if len(selectSources) < self.config.minKernelSources:
815 self.
log.error(
"Too few sources to calculate the PSF matching kernel: "
816 "%i selected but %i needed for the calculation.",
817 len(selectSources), self.config.minKernelSources)
818 raise RuntimeError(
"Cannot compute PSF matching kernel: too few sources selected.")
819 self.metadata.add(
"nPsfSources", len(selectSources))
825 """Exclude sources that are located on masked pixels.
829 mask : `lsst.afw.image.Mask`
830 The image mask plane to use to reject sources
831 based on the location of their centroid on the ccd.
832 sources : `lsst.afw.table.SourceCatalog`
833 The source catalog to evaluate.
834 excludeMaskPlanes : `list` of `str`
835 List of the names of the mask planes to exclude.
839 flags : `numpy.ndarray` of `bool`
840 Array indicating whether each source in the catalog should be
841 kept (True) or rejected (False) based on the value of the
842 mask plane at its location.
844 setExcludeMaskPlanes = [
845 maskPlane
for maskPlane
in excludeMaskPlanes
if maskPlane
in mask.getMaskPlaneDict()
848 excludePixelMask = mask.getPlaneBitMask(setExcludeMaskPlanes)
850 xv = np.rint(sources.getX() - mask.getX0())
851 yv = np.rint(sources.getY() - mask.getY0())
853 mv = mask.array[yv.astype(int), xv.astype(int)]
854 flags = np.bitwise_and(mv, excludePixelMask) == 0
858 """Perform preparatory calculations common to all Alard&Lupton Tasks.
862 template : `lsst.afw.image.ExposureF`
863 Template exposure, warped to match the science exposure. The
864 variance plane of the template image is modified in place.
865 science : `lsst.afw.image.ExposureF`
866 Science exposure to subtract from the template. The variance plane
867 of the science image is modified in place.
868 visitSummary : `lsst.afw.table.ExposureCatalog`, optional
869 Exposure catalog with external calibrations to be applied. Catalog
870 uses the detector id for the catalog id, sorted on id for fast
874 if visitSummary
is not None:
877 template[science.getBBox()], self.
log,
878 requiredTemplateFraction=self.config.requiredTemplateFraction,
879 exceptionMessage=
"Not attempting subtraction. To force subtraction,"
880 " set config requiredTemplateFraction=0"
882 self.metadata.add(
"templateCoveragePercent", 100*templateCoverageFraction)
884 if self.config.doScaleVariance:
888 templateVarFactor = self.scaleVariance.run(template.maskedImage)
889 sciVarFactor = self.scaleVariance.run(science.maskedImage)
890 self.
log.info(
"Template variance scaling factor: %.2f", templateVarFactor)
891 self.metadata.add(
"scaleTemplateVarianceFactor", templateVarFactor)
892 self.
log.info(
"Science variance scaling factor: %.2f", sciVarFactor)
893 self.metadata.add(
"scaleScienceVarianceFactor", sciVarFactor)
900 """Update the science and template mask planes before differencing.
904 template : `lsst.afw.image.Exposure`
905 Template exposure, warped to match the science exposure.
906 The template mask planes will be erased, except for a few specified
908 science : `lsst.afw.image.Exposure`
909 Science exposure to subtract from the template.
910 The DETECTED and DETECTED_NEGATIVE mask planes of the science image
913 self.
_clearMask(science.mask, clearMaskPlanes=[
"DETECTED",
"DETECTED_NEGATIVE"])
920 clearMaskPlanes = [mp
for mp
in template.mask.getMaskPlaneDict().keys()
921 if mp
not in self.config.preserveTemplateMask]
922 renameMaskPlanes = [mp
for mp
in self.config.renameTemplateMask
923 if mp
in template.mask.getMaskPlaneDict().keys()]
928 if "FAKE" in science.mask.getMaskPlaneDict().keys():
929 self.
log.info(
"Adding injected mask plane to science image")
931 if "FAKE" in template.mask.getMaskPlaneDict().keys():
932 self.
log.info(
"Adding injected mask plane to template image")
934 if "INJECTED" in renameMaskPlanes:
935 renameMaskPlanes.remove(
"INJECTED")
936 if "INJECTED_TEMPLATE" in clearMaskPlanes:
937 clearMaskPlanes.remove(
"INJECTED_TEMPLATE")
939 for maskPlane
in renameMaskPlanes:
941 self.
_clearMask(template.mask, clearMaskPlanes=clearMaskPlanes)
945 """Rename a mask plane by adding the new name and copying the data.
949 mask : `lsst.afw.image.Mask`
950 The mask image to update in place.
952 The name of the existing mask plane to copy.
954 The new name of the mask plane that will be added.
955 If the mask plane already exists, it will be updated in place.
957 mask.addMaskPlane(newMaskPlane)
958 originBitMask = mask.getPlaneBitMask(maskPlane)
959 destinationBitMask = mask.getPlaneBitMask(newMaskPlane)
960 mask.array |= ((mask.array & originBitMask) > 0)*destinationBitMask
963 """Clear the mask plane of an exposure.
967 mask : `lsst.afw.image.Mask`
968 The mask plane to erase, which will be modified in place.
969 clearMaskPlanes : `list` of `str`, optional
970 Erase the specified mask planes.
971 If not supplied, the entire mask will be erased.
973 if clearMaskPlanes
is None:
974 clearMaskPlanes = list(mask.getMaskPlaneDict().keys())
976 bitMaskToClear = mask.getPlaneBitMask(clearMaskPlanes)
977 mask &= ~bitMaskToClear
981 SubtractScoreOutputConnections):
986 pipelineConnections=AlardLuptonPreconvolveSubtractConnections):
991 """Subtract a template from a science image, convolving the science image
992 before computing the kernel, and also convolving the template before
995 ConfigClass = AlardLuptonPreconvolveSubtractConfig
996 _DefaultName =
"alardLuptonPreconvolveSubtract"
998 def run(self, template, science, sources, visitSummary=None):
999 """Preconvolve the science image with its own PSF,
1000 convolve the template image with a PSF-matching kernel and subtract
1001 from the preconvolved science image.
1005 template : `lsst.afw.image.ExposureF`
1006 The template image, which has previously been warped to the science
1007 image. The template bbox will be padded by a few pixels compared to
1009 science : `lsst.afw.image.ExposureF`
1010 The science exposure.
1011 sources : `lsst.afw.table.SourceCatalog`
1012 Identified sources on the science exposure. This catalog is used to
1013 select sources in order to perform the AL PSF matching on stamp
1015 visitSummary : `lsst.afw.table.ExposureCatalog`, optional
1016 Exposure catalog with complete external calibrations. Catalog uses
1017 the detector id for the catalog id, sorted on id for fast lookup.
1021 results : `lsst.pipe.base.Struct`
1022 ``scoreExposure`` : `lsst.afw.image.ExposureF`
1023 Result of subtracting the convolved template and science
1024 images. Attached PSF is that of the original science image.
1025 ``matchedTemplate`` : `lsst.afw.image.ExposureF`
1026 Warped and PSF-matched template exposure. Attached PSF is that
1027 of the original science image.
1028 ``matchedScience`` : `lsst.afw.image.ExposureF`
1029 The science exposure after convolving with its own PSF.
1030 Attached PSF is that of the original science image.
1031 ``backgroundModel`` : `lsst.afw.math.Function2D`
1032 Background model that was fit while solving for the
1034 ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
1035 Final kernel used to PSF-match the template to the science
1038 self.
_prepareInputs(template, science, visitSummary=visitSummary)
1041 scienceKernel = science.psf.getKernel()
1043 interpolateBadMaskPlanes=
True)
1044 self.metadata.add(
"convolvedExposure",
"Preconvolution")
1047 subtractResults = self.
runPreconvolve(template, science, matchedScience,
1048 selectSources, scienceKernel)
1051 self.
loglog.warning(
"Failed to match template. Checking coverage")
1054 self.config.minTemplateFractionForExpectedSuccess,
1055 exceptionMessage=
"Template coverage lower than expected to succeed."
1056 f
" Failure is tolerable: {e}")
1060 return subtractResults
1062 def runPreconvolve(self, template, science, matchedScience, selectSources, preConvKernel):
1063 """Convolve the science image with its own PSF, then convolve the
1064 template with a matching kernel and subtract to form the Score
1069 template : `lsst.afw.image.ExposureF`
1070 Template exposure, warped to match the science exposure.
1071 science : `lsst.afw.image.ExposureF`
1072 Science exposure to subtract from the template.
1073 matchedScience : `lsst.afw.image.ExposureF`
1074 The science exposure, convolved with the reflection of its own PSF.
1075 selectSources : `lsst.afw.table.SourceCatalog`
1076 Identified sources on the science exposure. This catalog is used to
1077 select sources in order to perform the AL PSF matching on stamp
1079 preConvKernel : `lsst.afw.math.Kernel`
1080 The reflection of the kernel that was used to preconvolve the
1081 `science` exposure. Must be normalized to sum to 1.
1085 results : `lsst.pipe.base.Struct`
1087 ``scoreExposure`` : `lsst.afw.image.ExposureF`
1088 Result of subtracting the convolved template and science
1089 images. Attached PSF is that of the original science image.
1090 ``matchedTemplate`` : `lsst.afw.image.ExposureF`
1091 Warped and PSF-matched template exposure. Attached PSF is that
1092 of the original science image.
1093 ``matchedScience`` : `lsst.afw.image.ExposureF`
1094 The science exposure after convolving with its own PSF.
1095 Attached PSF is that of the original science image.
1096 ``backgroundModel`` : `lsst.afw.math.Function2D`
1097 Background model that was fit while solving for the
1099 ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
1100 Final kernel used to PSF-match the template to the science
1103 bbox = science.getBBox()
1104 innerBBox = preConvKernel.shrinkBBox(bbox)
1106 kernelSources = self.makeKernel.selectKernelSources(template[innerBBox], matchedScience[innerBBox],
1107 candidateList=selectSources,
1109 kernelResult = self.makeKernel.run(template[innerBBox], matchedScience[innerBBox], kernelSources,
1112 matchedTemplate = self.
_convolveExposure(template, kernelResult.psfMatchingKernel,
1116 interpolateBadMaskPlanes=
True,
1117 photoCalib=science.photoCalib)
1119 backgroundModel=(kernelResult.backgroundModel
1120 if self.config.doSubtractBackground
else None))
1121 correctedScore = self.
finalize(template[bbox], science, score,
1122 kernelResult.psfMatchingKernel,
1123 templateMatched=
True, preConvMode=
True,
1124 preConvKernel=preConvKernel)
1126 return lsst.pipe.base.Struct(scoreExposure=correctedScore,
1127 matchedTemplate=matchedTemplate,
1128 matchedScience=matchedScience,
1129 backgroundModel=kernelResult.backgroundModel,
1130 psfMatchingKernel=kernelResult.psfMatchingKernel)
1134 exceptionMessage=""):
1135 """Raise NoWorkFound if template coverage < requiredTemplateFraction
1139 templateExposure : `lsst.afw.image.ExposureF`
1140 The template exposure to check
1141 logger : `logging.Logger`
1142 Logger for printing output.
1143 requiredTemplateFraction : `float`, optional
1144 Fraction of pixels of the science image required to have coverage
1146 exceptionMessage : `str`, optional
1147 Message to include in the exception raised if the template coverage
1152 templateCoverageFraction: `float`
1153 Fraction of pixels in the template with data.
1157 lsst.pipe.base.NoWorkFound
1158 Raised if fraction of good pixels, defined as not having NO_DATA
1159 set, is less than the requiredTemplateFraction
1163 pixNoData = np.count_nonzero(templateExposure.mask.array
1164 & templateExposure.mask.getPlaneBitMask(
'NO_DATA'))
1165 pixGood = templateExposure.getBBox().getArea() - pixNoData
1166 templateCoverageFraction = pixGood/templateExposure.getBBox().getArea()
1167 logger.info(
"template has %d good pixels (%.1f%%)", pixGood, 100*templateCoverageFraction)
1169 if templateCoverageFraction < requiredTemplateFraction:
1170 message = (
"Insufficient Template Coverage. (%.1f%% < %.1f%%)" % (
1171 100*templateCoverageFraction,
1172 100*requiredTemplateFraction))
1173 raise lsst.pipe.base.NoWorkFound(message +
" " + exceptionMessage)
1174 return templateCoverageFraction
1178 """Subtract template from science, propagating relevant metadata.
1182 science : `lsst.afw.Exposure`
1183 The input science image.
1184 template : `lsst.afw.Exposure`
1185 The template to subtract from the science image.
1186 backgroundModel : `lsst.afw.MaskedImage`, optional
1187 Differential background model
1191 difference : `lsst.afw.Exposure`
1192 The subtracted image.
1194 difference = science.clone()
1195 if backgroundModel
is not None:
1196 difference.maskedImage -= backgroundModel
1197 difference.maskedImage -= template.maskedImage
1202 """Determine that the PSF of ``exp1`` is not wider than that of ``exp2``.
1206 exp1 : `~lsst.afw.image.Exposure`
1207 Exposure with the reference point spread function (PSF) to evaluate.
1208 exp2 : `~lsst.afw.image.Exposure`
1209 Exposure with a candidate point spread function (PSF) to evaluate.
1210 fwhmExposureBuffer : `float`
1211 Fractional buffer margin to be left out of all sides of the image
1212 during the construction of the grid to compute mean PSF FWHM in an
1213 exposure, if the PSF is not available at its average position.
1214 fwhmExposureGrid : `int`
1215 Grid size to compute the mean FWHM in an exposure, if the PSF is not
1216 available at its average position.
1220 True if ``exp1`` has a PSF that is not wider than that of ``exp2`` in
1224 shape1 = getPsfFwhm(exp1.psf, average=
False)
1225 shape2 = getPsfFwhm(exp2.psf, average=
False)
1227 shape1 = evaluateMeanPsfFwhm(exp1,
1228 fwhmExposureBuffer=fwhmExposureBuffer,
1229 fwhmExposureGrid=fwhmExposureGrid
1231 shape2 = evaluateMeanPsfFwhm(exp2,
1232 fwhmExposureBuffer=fwhmExposureBuffer,
1233 fwhmExposureGrid=fwhmExposureGrid
1235 return shape1 <= shape2
1238 xTest = shape1[0] <= shape2[0]
1239 yTest = shape1[1] <= shape2[1]
1240 return xTest | yTest
1244 """Replace masked image pixels with interpolated values.
1248 maskedImage : `lsst.afw.image.MaskedImage`
1249 Image on which to perform interpolation.
1250 badMaskPlanes : `list` of `str`
1251 List of mask planes to interpolate over.
1252 fallbackValue : `float`, optional
1253 Value to set when interpolation fails.
1258 The number of masked pixels that were replaced.
1260 imgBadMaskPlanes = [
1261 maskPlane
for maskPlane
in badMaskPlanes
if maskPlane
in maskedImage.mask.getMaskPlaneDict()
1264 image = maskedImage.image.array
1265 badPixels = (maskedImage.mask.array & maskedImage.mask.getPlaneBitMask(imgBadMaskPlanes)) > 0
1266 image[badPixels] = np.nan
1267 if fallbackValue
is None:
1268 fallbackValue = np.nanmedian(image)
1271 image[badPixels] = fallbackValue
1272 return np.sum(badPixels)
Parameters to control convolution.
runPreconvolve(self, template, science, matchedScience, selectSources, preConvKernel)
run(self, template, science, sources, visitSummary=None)
_clearMask(self, mask, clearMaskPlanes=None)
_prepareInputs(self, template, science, visitSummary=None)
_calculateMagLim(self, exposure, nsigma=5.0, fallbackPsfSize=None)
runConvolveTemplate(self, template, science, selectSources)
_applyExternalCalibrations(self, exposure, visitSummary)
updateMasks(self, template, science)
_sourceSelector(self, sources, mask)
_convolveExposure(self, exposure, kernel, convolutionControl, bbox=None, psf=None, photoCalib=None, interpolateBadMaskPlanes=False)
runConvolveScience(self, template, science, selectSources)
run(self, template, science, sources, finalizedPsfApCorrCatalog=None, visitSummary=None)
finalize(self, template, science, difference, kernel, templateMatched=True, preConvMode=False, preConvKernel=None, spatiallyVarying=False)
_validateExposures(template, science)
_checkMask(mask, sources, excludeMaskPlanes)
_renameMaskPlanes(mask, maskPlane, newMaskPlane)
Provides consistent interface for LSST exceptions.
Reports invalid arguments.
Reports when the result of an operation cannot be represented by the destination type.
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)