22from astropy
import units
as u
35from .
import MakeKernelTask, DecorrelateALKernelTask
36from lsst.utils.timer
import timeMethod
38__all__ = [
"AlardLuptonSubtractConfig",
"AlardLuptonSubtractTask",
39 "AlardLuptonPreconvolveSubtractConfig",
"AlardLuptonPreconvolveSubtractTask"]
41_dimensions = (
"instrument",
"visit",
"detector")
42_defaultTemplates = {
"coaddName":
"deep",
"fakesType":
""}
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"
54 science = connectionTypes.Input(
55 doc=
"Input science exposure to subtract from.",
56 dimensions=(
"instrument",
"visit",
"detector"),
57 storageClass=
"ExposureF",
58 name=
"{fakesType}calexp"
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",
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",
78 if not config.doApplyExternalCalibrations:
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",
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",
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",
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"
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",
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",
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"
140 target=MakeKernelTask,
141 doc=
"Task to construct a matching kernel for convolution.",
146 doc=
"Perform diffim decorrelation to undo pixel correlation due to A&L "
147 "kernel convolution? If True, also update the diffim PSF."
150 target=DecorrelateALKernelTask,
151 doc=
"Task to decorrelate the image difference.",
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."
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."
169 doc=
"Scale variance of the image difference?"
172 target=ScaleVarianceTask,
173 doc=
"Subtask to rescale the variance of the template to the statistically expected level."
176 doc=
"Subtract the background fit when solving the kernel?",
182 "Replace science Exposure's calibration objects with those"
183 " in visitSummary. Ignored if `doApplyFinalizedPsf is True."
189 target=ScienceSourceSelectorTask,
190 doc=
"Task to select sources to be used for PSF matching.",
195 doc=
"Minimum signal to noise ratio of detected sources "
196 "to use for calculating the PSF matching kernel."
201 doc=
"Maximum signal to noise ratio of detected sources "
202 "to use for calculating the PSF matching kernel."
207 doc=
"Maximum number of sources to use for calculating the PSF matching kernel."
208 "Set to -1 to disable."
213 doc=
"Minimum number of sources needed for calculating the PSF matching kernel."
217 default=(
"NO_DATA",
"BAD",
"SAT",
"EDGE",
"FAKE"),
218 doc=
"Mask planes to exclude when selecting sources for PSF matching."
222 default=(
"NO_DATA",
"BAD",
"SAT",
"EDGE"),
223 doc=
"Mask planes to interpolate over."
227 default=(
"NO_DATA",
"BAD",),
228 doc=
"Mask planes from the template to propagate to the image difference."
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."
239 doc=
"Re-run source detection for kernel candidates if an error is"
240 " encountered while calculating the matching kernel."
247 self.
makeKernel.kernel.active.fitForBackground =
True
248 self.
makeKernel.kernel.active.spatialKernelOrder = 1
249 self.
makeKernel.kernel.active.spatialBgOrder = 2
260 pipelineConnections=AlardLuptonSubtractConnections):
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."
272 """Compute the image difference of a science and template image using
273 the Alard & Lupton (1998) algorithm.
275 ConfigClass = AlardLuptonSubtractConfig
276 _DefaultName =
"alardLuptonSubtract"
280 self.makeSubtask(
"decorrelate")
281 self.makeSubtask(
"makeKernel")
282 self.makeSubtask(
"sourceSelector")
283 if self.config.doScaleVariance:
284 self.makeSubtask(
"scaleVariance")
293 """Replace calibrations (psf, and ApCorrMap) on this exposure with
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
307 exposure : `lsst.afw.image.exposure.Exposure`
308 Exposure with adjusted calibrations.
310 detectorId = exposure.info.getDetector().getId()
312 row = visitSummary.find(detectorId)
314 self.
log.warning(
"Detector id %s not found in external calibrations catalog; "
315 "Using original calibrations.", detectorId)
318 apCorrMap = row.getApCorrMap()
320 self.
log.warning(
"Detector id %s has None for psf in "
321 "external calibrations catalog; Using original psf and aperture correction.",
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.",
329 exposure.info.setApCorrMap(apCorrMap)
334 def run(self, template, science, sources, visitSummary=None):
335 """PSF match, subtract, and decorrelate two images.
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
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
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
362 ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
363 Kernel used to PSF-match the convolved image.
368 If an unsupported convolution mode is supplied.
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
378 fwhmExposureBuffer = self.config.makeKernel.fwhmExposureBuffer
379 fwhmExposureGrid = self.config.makeKernel.fwhmExposureGrid
396 self.
log.info(
"Unable to evaluate PSF at the average position. "
397 "Evaluting PSF on a grid of points."
400 fwhmExposureBuffer=fwhmExposureBuffer,
401 fwhmExposureGrid=fwhmExposureGrid
404 fwhmExposureBuffer=fwhmExposureBuffer,
405 fwhmExposureGrid=fwhmExposureGrid
414 fluxlim_science = (maglim_science*u.ABmag).to_value(u.nJy)
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
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
426 if self.config.mode ==
"auto":
429 fwhmExposureBuffer=fwhmExposureBuffer,
430 fwhmExposureGrid=fwhmExposureGrid)
433 self.
log.info(
"Average template PSF size is greater, "
434 "but science PSF greater in one dimension: convolving template image.")
436 self.
log.info(
"Science PSF size is greater: convolving template image.")
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
446 raise RuntimeError(
"Cannot handle AlardLuptonSubtract mode: %s", self.config.mode)
449 sourceMask = science.mask.clone()
450 sourceMask.array |= template[science.getBBox()].mask.array
453 self.metadata[
"convolvedExposure"] =
"Template"
456 self.metadata[
"convolvedExposure"] =
"Science"
460 self.
log.warning(
"Failed to match template. Checking coverage")
463 self.config.minTemplateFractionForExpectedSuccess,
464 exceptionMessage=
"Template coverage lower than expected to succeed."
465 f
" Failure is tolerable: {e}")
471 return subtractResults
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
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.
491 The task metadata does not include docstrings, so descriptions of the
492 computed metrics are given here:
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.
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.
518 science_footprints = np.zeros(n)
519 difference_footprints = np.zeros(n)
521 for i, record
in enumerate(sources):
522 footprint = record.getFootprint()
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
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())
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"])
544 """Convolve the template image with a PSF-matching kernel and subtract
545 from the science image.
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
560 results : `lsst.pipe.base.Struct`
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.
572 kernelSources = self.makeKernel.selectKernelSources(template, science,
573 candidateList=selectSources,
575 kernelResult = self.makeKernel.run(template, science, kernelSources,
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(
583 sigmaToFwhm = 2*np.log(2*np.sqrt(2))
584 candidateList = self.makeKernel.makeCandidateList(template, science, kernelSize,
587 kernelSources = self.makeKernel.selectKernelSources(template, science,
588 candidateList=candidateList,
590 kernelResult = self.makeKernel.run(template, science, kernelSources,
595 matchedTemplate = self.
_convolveExposure(template, kernelResult.psfMatchingKernel,
597 bbox=science.getBBox(),
599 photoCalib=science.photoCalib)
602 backgroundModel=(kernelResult.backgroundModel
603 if self.config.doSubtractBackground
else None))
604 correctedExposure = self.
finalize(template, science, difference,
605 kernelResult.psfMatchingKernel,
606 templateMatched=
True)
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)
616 """Convolve the science image with a PSF-matching kernel and subtract
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
632 results : `lsst.pipe.base.Struct`
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.
644 bbox = science.getBBox()
645 kernelSources = self.makeKernel.selectKernelSources(science, template,
646 candidateList=selectSources,
648 kernelResult = self.makeKernel.run(science, template, kernelSources,
650 modelParams = kernelResult.backgroundModel.getParameters()
652 kernelResult.backgroundModel.setParameters([-p
for p
in modelParams])
654 kernelImage = lsst.afw.image.ImageD(kernelResult.psfMatchingKernel.getDimensions())
655 norm = kernelResult.psfMatchingKernel.computeImage(kernelImage, doNormalize=
False)
662 matchedScience.maskedImage /= norm
663 matchedTemplate = template.clone()[bbox]
664 matchedTemplate.maskedImage /= norm
665 matchedTemplate.setPhotoCalib(science.photoCalib)
668 backgroundModel=(kernelResult.backgroundModel
669 if self.config.doSubtractBackground
else None))
671 correctedExposure = self.
finalize(template, science, difference,
672 kernelResult.psfMatchingKernel,
673 templateMatched=
False)
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)
682 def finalize(self, template, science, difference, kernel,
683 templateMatched=True,
686 spatiallyVarying=False):
687 """Decorrelate the difference image to undo the noise correlations
688 caused by convolution.
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?
713 correctedExposure : `lsst.afw.image.ExposureF`
714 The decorrelated image difference.
716 if self.config.doDecorrelation:
717 self.
log.info(
"Decorrelating image difference.")
721 correctedExposure = self.decorrelate.run(science, template[science.getBBox()], difference, kernel,
722 templateMatched=templateMatched,
723 preConvMode=preConvMode,
724 preConvKernel=preConvKernel,
725 spatiallyVarying=spatiallyVarying).correctedExposure
727 self.
log.info(
"NOT decorrelating image difference.")
728 correctedExposure = difference
729 return correctedExposure
732 """Calculate an exposure's limiting magnitude.
734 This method uses the photometric zeropoint together with the
735 PSF size from the average position of the exposure.
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.
748 maglim : `astropy.units.Quantity`
749 The limiting magnitude of the exposure, or np.nan.
752 psf = exposure.getPsf()
753 psf_shape = psf.computeShape(psf.getAveragePosition())
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))
761 self.
log.info(
"Unable to evaluate PSF, setting maglim to nan")
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))
773 """Check that the WCS of the two Exposures match, the template bbox
774 contains the science bbox, and that the bands match.
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.
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.
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)
798 assert templateBBox.contains(scienceBBox), \
799 "Template bbox does not contain all of the science image."
805 interpolateBadMaskPlanes=False,
807 """Convolve an exposure with the given kernel.
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.
826 convolvedExp : `lsst.afw.Exposure`
829 convolvedExposure = exposure.clone()
831 convolvedExposure.setPsf(psf)
832 if photoCalib
is not None:
833 convolvedExposure.setPhotoCalib(photoCalib)
834 if interpolateBadMaskPlanes
and self.config.badMaskPlanes
is not None:
836 self.config.badMaskPlanes)
837 self.metadata[
"nInterpolated"] = nInterp
838 convolvedImage = lsst.afw.image.MaskedImageF(convolvedExposure.getBBox())
840 convolvedExposure.setMaskedImage(convolvedImage)
842 return convolvedExposure
844 return convolvedExposure[bbox]
847 """Select sources from a catalog that meet the selection criteria.
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.
859 selectSources : `lsst.afw.table.SourceCatalog`
860 The input source catalog, with flagged and low signal-to-noise
866 If there are too few sources to compute the PSF matching kernel
867 remaining after source selection.
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)
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)
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)
900 def _checkMask(mask, sources, excludeMaskPlanes, checkAdjacent=True):
901 """Exclude sources that are located on or adjacent to masked pixels.
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.
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.
920 setExcludeMaskPlanes = [
921 maskPlane
for maskPlane
in excludeMaskPlanes
if maskPlane
in mask.getMaskPlaneDict()
924 excludePixelMask = mask.getPlaneBitMask(setExcludeMaskPlanes)
926 xv = (np.rint(sources.getX() - mask.getX0())).astype(int)
927 yv = (np.rint(sources.getY() - mask.getY0())).astype(int)
929 flags = np.ones(len(sources), dtype=bool)
931 pixRange = (0, -1, 1)
936 mv = mask.array[yv + j, xv + i]
937 flags *= np.bitwise_and(mv, excludePixelMask) == 0
941 """Perform preparatory calculations common to all Alard&Lupton Tasks.
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
957 if visitSummary
is not None:
960 template[science.getBBox()], self.
log,
961 requiredTemplateFraction=self.config.requiredTemplateFraction,
962 exceptionMessage=
"Not attempting subtraction. To force subtraction,"
963 " set config requiredTemplateFraction=0"
965 self.metadata[
"templateCoveragePercent"] = 100*templateCoverageFraction
967 if self.config.doScaleVariance:
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
983 """Update the science and template mask planes before differencing.
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
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
996 self.
_clearMask(science.mask, clearMaskPlanes=[
"DETECTED",
"DETECTED_NEGATIVE"])
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()]
1011 if "FAKE" in science.mask.getMaskPlaneDict().keys():
1012 self.
log.info(
"Adding injected mask plane to science image")
1014 if "FAKE" in template.mask.getMaskPlaneDict().keys():
1015 self.
log.info(
"Adding injected mask plane to template image")
1017 if "INJECTED" in renameMaskPlanes:
1018 renameMaskPlanes.remove(
"INJECTED")
1019 if "INJECTED_TEMPLATE" in clearMaskPlanes:
1020 clearMaskPlanes.remove(
"INJECTED_TEMPLATE")
1022 for maskPlane
in renameMaskPlanes:
1024 self.
_clearMask(template.mask, clearMaskPlanes=clearMaskPlanes)
1028 """Rename a mask plane by adding the new name and copying the data.
1032 mask : `lsst.afw.image.Mask`
1033 The mask image to update in place.
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.
1040 mask.addMaskPlane(newMaskPlane)
1041 originBitMask = mask.getPlaneBitMask(maskPlane)
1042 destinationBitMask = mask.getPlaneBitMask(newMaskPlane)
1043 mask.array |= ((mask.array & originBitMask) > 0)*destinationBitMask
1046 """Clear the mask plane of an exposure.
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.
1056 if clearMaskPlanes
is None:
1057 clearMaskPlanes = list(mask.getMaskPlaneDict().keys())
1059 bitMaskToClear = mask.getPlaneBitMask(clearMaskPlanes)
1060 mask &= ~bitMaskToClear
1064 SubtractScoreOutputConnections):
1069 pipelineConnections=AlardLuptonPreconvolveSubtractConnections):
1074 """Subtract a template from a science image, convolving the science image
1075 before computing the kernel, and also convolving the template before
1078 ConfigClass = AlardLuptonPreconvolveSubtractConfig
1079 _DefaultName =
"alardLuptonPreconvolveSubtract"
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.
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
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
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.
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
1117 ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
1118 Final kernel used to PSF-match the template to the science
1121 self.
_prepareInputs(template, science, visitSummary=visitSummary)
1124 scienceKernel = science.psf.getKernel()
1126 interpolateBadMaskPlanes=
True)
1127 self.metadata[
"convolvedExposure"] =
"Preconvolution"
1130 subtractResults = self.
runPreconvolve(template, science, matchedScience,
1131 selectSources, scienceKernel)
1134 self.
log.warning(
"Failed to match template. Checking coverage")
1137 self.config.minTemplateFractionForExpectedSuccess,
1138 exceptionMessage=
"Template coverage lower than expected to succeed."
1139 f
" Failure is tolerable: {e}")
1143 return subtractResults
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
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
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.
1168 results : `lsst.pipe.base.Struct`
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
1182 ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
1183 Final kernel used to PSF-match the template to the science
1186 bbox = science.getBBox()
1187 innerBBox = preConvKernel.shrinkBBox(bbox)
1189 kernelSources = self.makeKernel.selectKernelSources(template[innerBBox], matchedScience[innerBBox],
1190 candidateList=selectSources,
1192 kernelResult = self.makeKernel.run(template[innerBBox], matchedScience[innerBBox], kernelSources,
1195 matchedTemplate = self.
_convolveExposure(template, kernelResult.psfMatchingKernel,
1199 interpolateBadMaskPlanes=
True,
1200 photoCalib=science.photoCalib)
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)
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)
1218 exceptionMessage=""):
1219 """Raise NoWorkFound if template coverage < requiredTemplateFraction
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
1230 exceptionMessage : `str`, optional
1231 Message to include in the exception raised if the template coverage
1236 templateCoverageFraction: `float`
1237 Fraction of pixels in the template with data.
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
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)
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
1262 """Subtract template from science, propagating relevant metadata.
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
1275 difference : `lsst.afw.Exposure`
1276 The subtracted image.
1278 difference = science.clone()
1279 if backgroundModel
is not None:
1280 difference.maskedImage -= backgroundModel
1281 difference.maskedImage -= template.maskedImage
1286 """Determine that the PSF of ``exp1`` is not wider than that of ``exp2``.
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.
1304 True if ``exp1`` has a PSF that is not wider than that of ``exp2`` in
1308 shape1 = getPsfFwhm(exp1.psf, average=
False)
1309 shape2 = getPsfFwhm(exp2.psf, average=
False)
1311 shape1 = evaluateMeanPsfFwhm(exp1,
1312 fwhmExposureBuffer=fwhmExposureBuffer,
1313 fwhmExposureGrid=fwhmExposureGrid
1315 shape2 = evaluateMeanPsfFwhm(exp2,
1316 fwhmExposureBuffer=fwhmExposureBuffer,
1317 fwhmExposureGrid=fwhmExposureGrid
1319 return shape1 <= shape2
1322 xTest = shape1[0] <= shape2[0]
1323 yTest = shape1[1] <= shape2[1]
1324 return xTest | yTest
1328 """Replace masked image pixels with interpolated values.
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.
1342 The number of masked pixels that were replaced.
1344 imgBadMaskPlanes = [
1345 maskPlane
for maskPlane
in badMaskPlanes
if maskPlane
in maskedImage.mask.getMaskPlaneDict()
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)
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)
_clearMask(self, mask, clearMaskPlanes=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)
_applyExternalCalibrations(self, exposure, visitSummary)
_checkMask(mask, sources, excludeMaskPlanes, checkAdjacent=True)
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)
finalize(self, template, science, difference, kernel, templateMatched=True, preConvMode=False, preConvKernel=None, spatiallyVarying=False)
_validateExposures(template, science)
_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.
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)