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."
246 self.
makeKernel.kernel.active.spatialKernelOrder = 1
247 self.
makeKernel.kernel.active.spatialBgOrder = 2
258 pipelineConnections=AlardLuptonSubtractConnections):
261 default=
"convolveTemplate",
262 allowed={
"auto":
"Choose which image to convolve at runtime.",
263 "convolveScience":
"Only convolve the science image.",
264 "convolveTemplate":
"Only convolve the template image."},
265 doc=
"Choose which image to convolve at runtime, or require that a specific image is convolved."
270 """Compute the image difference of a science and template image using
271 the Alard & Lupton (1998) algorithm.
273 ConfigClass = AlardLuptonSubtractConfig
274 _DefaultName =
"alardLuptonSubtract"
278 self.makeSubtask(
"decorrelate")
279 self.makeSubtask(
"makeKernel")
280 self.makeSubtask(
"sourceSelector")
281 if self.config.doScaleVariance:
282 self.makeSubtask(
"scaleVariance")
291 """Replace calibrations (psf, and ApCorrMap) on this exposure with
296 exposure : `lsst.afw.image.exposure.Exposure`
297 Input exposure to adjust calibrations.
298 visitSummary : `lsst.afw.table.ExposureCatalog`
299 Exposure catalog with external calibrations to be applied. Catalog
300 uses the detector id for the catalog id, sorted on id for fast
305 exposure : `lsst.afw.image.exposure.Exposure`
306 Exposure with adjusted calibrations.
308 detectorId = exposure.info.getDetector().getId()
310 row = visitSummary.find(detectorId)
312 self.
log.warning(
"Detector id %s not found in external calibrations catalog; "
313 "Using original calibrations.", detectorId)
316 apCorrMap = row.getApCorrMap()
318 self.
log.warning(
"Detector id %s has None for psf in "
319 "external calibrations catalog; Using original psf and aperture correction.",
321 elif apCorrMap
is None:
322 self.
log.warning(
"Detector id %s has None for apCorrMap in "
323 "external calibrations catalog; Using original psf and aperture correction.",
327 exposure.info.setApCorrMap(apCorrMap)
332 def run(self, template, science, sources, visitSummary=None):
333 """PSF match, subtract, and decorrelate two images.
337 template : `lsst.afw.image.ExposureF`
338 Template exposure, warped to match the science exposure.
339 science : `lsst.afw.image.ExposureF`
340 Science exposure to subtract from the template.
341 sources : `lsst.afw.table.SourceCatalog`
342 Identified sources on the science exposure. This catalog is used to
343 select sources in order to perform the AL PSF matching on stamp
345 visitSummary : `lsst.afw.table.ExposureCatalog`, optional
346 Exposure catalog with external calibrations to be applied. Catalog
347 uses the detector id for the catalog id, sorted on id for fast
352 results : `lsst.pipe.base.Struct`
353 ``difference`` : `lsst.afw.image.ExposureF`
354 Result of subtracting template and science.
355 ``matchedTemplate`` : `lsst.afw.image.ExposureF`
356 Warped and PSF-matched template exposure.
357 ``backgroundModel`` : `lsst.afw.math.Function2D`
358 Background model that was fit while solving for the
360 ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
361 Kernel used to PSF-match the convolved image.
366 If an unsupported convolution mode is supplied.
368 If there are too few sources to calculate the PSF matching kernel.
369 lsst.pipe.base.NoWorkFound
370 Raised if fraction of good pixels, defined as not having NO_DATA
371 set, is less then the configured requiredTemplateFraction
376 fwhmExposureBuffer = self.config.makeKernel.fwhmExposureBuffer
377 fwhmExposureGrid = self.config.makeKernel.fwhmExposureGrid
394 self.
log.info(
"Unable to evaluate PSF at the average position. "
395 "Evaluting PSF on a grid of points."
398 fwhmExposureBuffer=fwhmExposureBuffer,
399 fwhmExposureGrid=fwhmExposureGrid
402 fwhmExposureBuffer=fwhmExposureBuffer,
403 fwhmExposureGrid=fwhmExposureGrid
412 fluxlim_science = (maglim_science*u.ABmag).to_value(u.nJy)
414 if np.isnan(maglim_template):
415 self.
log.info(
"Cannot evaluate template limiting mag; adopting science limiting mag for diffim")
416 maglim_diffim = maglim_science
418 fluxlim_template = (maglim_template*u.ABmag).to_value(u.nJy)
419 maglim_diffim = (np.sqrt(fluxlim_science**2 + fluxlim_template**2)*u.nJy).to(u.ABmag).value
420 self.metadata[
"scienceLimitingMagnitude"] = maglim_science
421 self.metadata[
"templateLimitingMagnitude"] = maglim_template
422 self.metadata[
"diffimLimitingMagnitude"] = maglim_diffim
424 if self.config.mode ==
"auto":
427 fwhmExposureBuffer=fwhmExposureBuffer,
428 fwhmExposureGrid=fwhmExposureGrid)
431 self.
log.info(
"Average template PSF size is greater, "
432 "but science PSF greater in one dimension: convolving template image.")
434 self.
log.info(
"Science PSF size is greater: convolving template image.")
436 self.
log.info(
"Template PSF size is greater: convolving science image.")
437 elif self.config.mode ==
"convolveTemplate":
438 self.
log.info(
"`convolveTemplate` is set: convolving template image.")
439 convolveTemplate =
True
440 elif self.config.mode ==
"convolveScience":
441 self.
log.info(
"`convolveScience` is set: convolving science image.")
442 convolveTemplate =
False
444 raise RuntimeError(
"Cannot handle AlardLuptonSubtract mode: %s", self.config.mode)
447 sourceMask = science.mask.clone()
448 sourceMask.array |= template[science.getBBox()].mask.array
451 self.metadata[
"convolvedExposure"] =
"Template"
454 self.metadata[
"convolvedExposure"] =
"Science"
458 self.
log.warning(
"Failed to match template. Checking coverage")
461 self.config.minTemplateFractionForExpectedSuccess,
462 exceptionMessage=
"Template coverage lower than expected to succeed."
463 f
" Failure is tolerable: {e}")
469 return subtractResults
332 def run(self, template, science, sources, visitSummary=None):
…
472 r"""Compute quality metrics (saved to the task metadata) on the
473 difference image, at the locations of detected stars on the science
474 image. This restricts the metric to locations that should be
479 science : `lsst.afw.image.ExposureF`
480 Science exposure that was subtracted.
481 difference : `lsst.afw.image.ExposureF`
482 Result of subtracting template and science.
483 stars : `lsst.afw.table.SourceCatalog`
484 Good calibration sources detected on science image; these
485 footprints are what the metrics are computed on.
489 The task metadata does not include docstrings, so descriptions of the
490 computed metrics are given here:
492 differenceFootprintRatioMean
493 Mean of the ratio of the absolute value of the difference image
494 (with the mean absolute value of the sky regions on the difference
495 image removed) to the science image, computed in the footprints
496 of stars detected on the science image (the sums below are of the
497 pixels in each star or sky footprint):
498 :math:`\mathrm{mean}_{footprints}((\sum |difference| -
499 \mathrm{mean}(\sum |difference_{sky}|)) / \sum science)`
500 differenceFootprintRatioStdev
501 Standard Deviation across footprints of the above ratio.
502 differenceFootprintSkyRatioMean
503 Mean of the ratio of the absolute value of sky source regions on
504 the difference image to the science image (the sum below is of the
505 pixels in each sky source footprint):
506 :math:`\mathrm{mean}_{footprints}(\sum |difference_{sky}| / \sum science_{sky})`
507 differenceFootprintSkyRatioStdev
508 Standard Deivation across footprints of the above sky ratio.
510 def footprint_mean(sources, sky=0):
511 """Compute ratio of the absolute value of the diffim to the science
512 image, within each source footprint, subtracting the sky from the
513 diffim values if provided.
516 science_footprints = np.zeros(n)
517 difference_footprints = np.zeros(n)
519 for i, record
in enumerate(sources):
520 footprint = record.getFootprint()
523 science_footprints[i] = abs(heavy.getImageArray()).mean()
524 difference_footprints[i] = abs(heavy_diff.getImageArray()).mean()
525 ratio[i] = abs((difference_footprints[i] - sky)/science_footprints[i])
526 return science_footprints, difference_footprints, ratio
528 sky = stars[
"sky_source"]
529 sky_science, sky_difference, sky_ratio = footprint_mean(stars[sky])
530 science_footprints, difference_footprints, ratio = footprint_mean(stars[~sky], sky_difference.mean())
532 self.metadata[
"differenceFootprintRatioMean"] = ratio.mean()
533 self.metadata[
"differenceFootprintRatioStdev"] = ratio.std()
534 self.metadata[
"differenceFootprintSkyRatioMean"] = sky_ratio.mean()
535 self.metadata[
"differenceFootprintSkyRatioStdev"] = sky_ratio.std()
536 self.
log.info(
"Mean, stdev of ratio of difference to science "
537 "pixels in star footprints: %5.4f, %5.4f",
538 self.metadata[
"differenceFootprintRatioMean"],
539 self.metadata[
"differenceFootprintRatioStdev"])
542 """Convolve the template image with a PSF-matching kernel and subtract
543 from the science image.
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 and PSF-matched template exposure.
564 ``backgroundModel`` : `lsst.afw.math.Function2D`
565 Background model that was fit while solving for the PSF-matching kernel
566 ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
567 Kernel used to PSF-match the template to the science image.
570 kernelSources = self.makeKernel.selectKernelSources(template, science,
571 candidateList=selectSources,
573 kernelResult = self.makeKernel.run(template, science, kernelSources,
575 except Exception
as e:
576 if self.config.allowKernelSourceDetection:
577 self.
log.warning(
"Error encountered trying to construct the matching kernel"
578 f
" Running source detection and retrying. {e}")
579 kernelSize = self.makeKernel.makeKernelBasisList(
581 sigmaToFwhm = 2*np.log(2*np.sqrt(2))
582 candidateList = self.makeKernel.makeCandidateList(template, science, kernelSize,
585 kernelSources = self.makeKernel.selectKernelSources(template, science,
586 candidateList=candidateList,
588 kernelResult = self.makeKernel.run(template, science, kernelSources,
593 matchedTemplate = self.
_convolveExposure(template, kernelResult.psfMatchingKernel,
595 bbox=science.getBBox(),
597 photoCalib=science.photoCalib)
600 backgroundModel=(kernelResult.backgroundModel
601 if self.config.doSubtractBackground
else None))
602 correctedExposure = self.
finalize(template, science, difference,
603 kernelResult.psfMatchingKernel,
604 templateMatched=
True)
606 return lsst.pipe.base.Struct(difference=correctedExposure,
607 matchedTemplate=matchedTemplate,
608 matchedScience=science,
609 backgroundModel=kernelResult.backgroundModel,
610 psfMatchingKernel=kernelResult.psfMatchingKernel,
611 kernelSources=kernelSources)
614 """Convolve the science image with a PSF-matching kernel and subtract
619 template : `lsst.afw.image.ExposureF`
620 Template exposure, warped to match the science exposure.
621 science : `lsst.afw.image.ExposureF`
622 Science exposure to subtract from the template.
623 selectSources : `lsst.afw.table.SourceCatalog`
624 Identified sources on the science exposure. This catalog is used to
625 select sources in order to perform the AL PSF matching on stamp
630 results : `lsst.pipe.base.Struct`
632 ``difference`` : `lsst.afw.image.ExposureF`
633 Result of subtracting template and science.
634 ``matchedTemplate`` : `lsst.afw.image.ExposureF`
635 Warped template exposure. Note that in this case, the template
636 is not PSF-matched to the science image.
637 ``backgroundModel`` : `lsst.afw.math.Function2D`
638 Background model that was fit while solving for the PSF-matching kernel
639 ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
640 Kernel used to PSF-match the science image to the template.
642 bbox = science.getBBox()
643 kernelSources = self.makeKernel.selectKernelSources(science, template,
644 candidateList=selectSources,
646 kernelResult = self.makeKernel.run(science, template, kernelSources,
648 modelParams = kernelResult.backgroundModel.getParameters()
650 kernelResult.backgroundModel.setParameters([-p
for p
in modelParams])
652 kernelImage = lsst.afw.image.ImageD(kernelResult.psfMatchingKernel.getDimensions())
653 norm = kernelResult.psfMatchingKernel.computeImage(kernelImage, doNormalize=
False)
660 matchedScience.maskedImage /= norm
661 matchedTemplate = template.clone()[bbox]
662 matchedTemplate.maskedImage /= norm
663 matchedTemplate.setPhotoCalib(science.photoCalib)
666 backgroundModel=(kernelResult.backgroundModel
667 if self.config.doSubtractBackground
else None))
669 correctedExposure = self.
finalize(template, science, difference,
670 kernelResult.psfMatchingKernel,
671 templateMatched=
False)
673 return lsst.pipe.base.Struct(difference=correctedExposure,
674 matchedTemplate=matchedTemplate,
675 matchedScience=matchedScience,
676 backgroundModel=kernelResult.backgroundModel,
677 psfMatchingKernel=kernelResult.psfMatchingKernel,
678 kernelSources=kernelSources)
680 def finalize(self, template, science, difference, kernel,
681 templateMatched=True,
684 spatiallyVarying=False):
685 """Decorrelate the difference image to undo the noise correlations
686 caused by convolution.
690 template : `lsst.afw.image.ExposureF`
691 Template exposure, warped to match the science exposure.
692 science : `lsst.afw.image.ExposureF`
693 Science exposure to subtract from the template.
694 difference : `lsst.afw.image.ExposureF`
695 Result of subtracting template and science.
696 kernel : `lsst.afw.math.Kernel`
697 An (optionally spatially-varying) PSF matching kernel
698 templateMatched : `bool`, optional
699 Was the template PSF-matched to the science image?
700 preConvMode : `bool`, optional
701 Was the science image preconvolved with its own PSF
702 before PSF matching the template?
703 preConvKernel : `lsst.afw.detection.Psf`, optional
704 If not `None`, then the science image was pre-convolved with
705 (the reflection of) this kernel. Must be normalized to sum to 1.
706 spatiallyVarying : `bool`, optional
707 Compute the decorrelation kernel spatially varying across the image?
711 correctedExposure : `lsst.afw.image.ExposureF`
712 The decorrelated image difference.
714 if self.config.doDecorrelation:
715 self.
log.info(
"Decorrelating image difference.")
719 correctedExposure = self.decorrelate.run(science, template[science.getBBox()], difference, kernel,
720 templateMatched=templateMatched,
721 preConvMode=preConvMode,
722 preConvKernel=preConvKernel,
723 spatiallyVarying=spatiallyVarying).correctedExposure
725 self.
log.info(
"NOT decorrelating image difference.")
726 correctedExposure = difference
727 return correctedExposure
680 def finalize(self, template, science, difference, kernel,
…
730 """Calculate an exposure's limiting magnitude.
732 This method uses the photometric zeropoint together with the
733 PSF size from the average position of the exposure.
737 exposure : `lsst.afw.image.Exposure`
738 The target exposure to calculate the limiting magnitude for.
739 nsigma : `float`, optional
740 The detection threshold in sigma.
741 fallbackPsfSize : `float`, optional
742 PSF FWHM to use in the event the exposure PSF cannot be retrieved.
746 maglim : `astropy.units.Quantity`
747 The limiting magnitude of the exposure, or np.nan.
750 psf = exposure.getPsf()
751 psf_shape = psf.computeShape(psf.getAveragePosition())
753 if fallbackPsfSize
is not None:
754 self.
log.info(
"Unable to evaluate PSF, using fallback FWHM %f", fallbackPsfSize)
755 psf_area = np.pi*(fallbackPsfSize/2)**2
756 zeropoint = exposure.getPhotoCalib().instFluxToMagnitude(1)
757 maglim = zeropoint - 2.5*np.log10(nsigma*np.sqrt(psf_area))
759 self.
log.info(
"Unable to evaluate PSF, setting maglim to nan")
763 psf_area = np.pi*np.sqrt(psf_shape.getIxx()*psf_shape.getIyy())
764 zeropoint = exposure.getPhotoCalib().instFluxToMagnitude(1)
765 maglim = zeropoint - 2.5*np.log10(nsigma*np.sqrt(psf_area))
771 """Check that the WCS of the two Exposures match, the template bbox
772 contains the science bbox, and that the bands match.
776 template : `lsst.afw.image.ExposureF`
777 Template exposure, warped to match the science exposure.
778 science : `lsst.afw.image.ExposureF`
779 Science exposure to subtract from the template.
784 Raised if the WCS of the template is not equal to the science WCS,
785 if the science image is not fully contained in the template
786 bounding box, or if the bands do not match.
788 assert template.wcs == science.wcs, \
789 "Template and science exposure WCS are not identical."
790 templateBBox = template.getBBox()
791 scienceBBox = science.getBBox()
792 assert science.filter.bandLabel == template.filter.bandLabel, \
793 "Science and template exposures have different bands: %s, %s" % \
794 (science.filter, template.filter)
796 assert templateBBox.contains(scienceBBox), \
797 "Template bbox does not contain all of the science image."
803 interpolateBadMaskPlanes=False,
805 """Convolve an exposure with the given kernel.
809 exposure : `lsst.afw.Exposure`
810 exposure to convolve.
811 kernel : `lsst.afw.math.LinearCombinationKernel`
812 PSF matching kernel computed in the ``makeKernel`` subtask.
813 convolutionControl : `lsst.afw.math.ConvolutionControl`
814 Configuration for convolve algorithm.
815 bbox : `lsst.geom.Box2I`, optional
816 Bounding box to trim the convolved exposure to.
817 psf : `lsst.afw.detection.Psf`, optional
818 Point spread function (PSF) to set for the convolved exposure.
819 photoCalib : `lsst.afw.image.PhotoCalib`, optional
820 Photometric calibration of the convolved exposure.
824 convolvedExp : `lsst.afw.Exposure`
827 convolvedExposure = exposure.clone()
829 convolvedExposure.setPsf(psf)
830 if photoCalib
is not None:
831 convolvedExposure.setPhotoCalib(photoCalib)
832 if interpolateBadMaskPlanes
and self.config.badMaskPlanes
is not None:
834 self.config.badMaskPlanes)
835 self.metadata[
"nInterpolated"] = nInterp
836 convolvedImage = lsst.afw.image.MaskedImageF(convolvedExposure.getBBox())
838 convolvedExposure.setMaskedImage(convolvedImage)
840 return convolvedExposure
842 return convolvedExposure[bbox]
845 """Select sources from a catalog that meet the selection criteria.
849 sources : `lsst.afw.table.SourceCatalog`
850 Input source catalog to select sources from.
851 mask : `lsst.afw.image.Mask`
852 The image mask plane to use to reject sources
853 based on their location on the ccd.
857 selectSources : `lsst.afw.table.SourceCatalog`
858 The input source catalog, with flagged and low signal-to-noise
864 If there are too few sources to compute the PSF matching kernel
865 remaining after source selection.
868 selected = self.sourceSelector.selectSources(sources).selected
869 nInitialSelected = np.count_nonzero(selected)
870 selected *= self.
_checkMask(mask, sources, self.config.excludeMaskPlanes)
871 nSelected = np.count_nonzero(selected)
872 self.
log.info(
"Rejecting %i candidate sources: an excluded template mask plane is set.",
873 nInitialSelected - nSelected)
874 selectSources = sources[selected].copy(deep=
True)
877 if (len(selectSources) > self.config.maxKernelSources) & (self.config.maxKernelSources > 0):
878 signalToNoise = selectSources.getPsfInstFlux()/selectSources.getPsfInstFluxErr()
879 indices = np.argsort(signalToNoise)
880 indices = indices[-self.config.maxKernelSources:]
881 selected = np.zeros(len(selectSources), dtype=bool)
882 selected[indices] =
True
883 selectSources = selectSources[selected].copy(deep=
True)
885 self.
log.info(
"%i/%i=%.1f%% of sources selected for PSF matching from the input catalog",
886 len(selectSources), len(sources), 100*len(selectSources)/len(sources))
887 if len(selectSources) < self.config.minKernelSources:
888 self.
log.error(
"Too few sources to calculate the PSF matching kernel: "
889 "%i selected but %i needed for the calculation.",
890 len(selectSources), self.config.minKernelSources)
891 if not self.config.allowKernelSourceDetection:
892 raise RuntimeError(
"Cannot compute PSF matching kernel: too few sources selected.")
893 self.metadata[
"nPsfSources"] = len(selectSources)
898 def _checkMask(mask, sources, excludeMaskPlanes, checkAdjacent=True):
899 """Exclude sources that are located on or adjacent to masked pixels.
903 mask : `lsst.afw.image.Mask`
904 The image mask plane to use to reject sources
905 based on the location of their centroid on the ccd.
906 sources : `lsst.afw.table.SourceCatalog`
907 The source catalog to evaluate.
908 excludeMaskPlanes : `list` of `str`
909 List of the names of the mask planes to exclude.
913 flags : `numpy.ndarray` of `bool`
914 Array indicating whether each source in the catalog should be
915 kept (True) or rejected (False) based on the value of the
916 mask plane at its location.
918 setExcludeMaskPlanes = [
919 maskPlane
for maskPlane
in excludeMaskPlanes
if maskPlane
in mask.getMaskPlaneDict()
922 excludePixelMask = mask.getPlaneBitMask(setExcludeMaskPlanes)
924 xv = (np.rint(sources.getX() - mask.getX0())).astype(int)
925 yv = (np.rint(sources.getY() - mask.getY0())).astype(int)
927 flags = np.ones(len(sources), dtype=bool)
929 pixRange = (0, -1, 1)
934 mv = mask.array[yv + j, xv + i]
935 flags *= np.bitwise_and(mv, excludePixelMask) == 0
898 def _checkMask(mask, sources, excludeMaskPlanes, checkAdjacent=True):
…
939 """Perform preparatory calculations common to all Alard&Lupton Tasks.
943 template : `lsst.afw.image.ExposureF`
944 Template exposure, warped to match the science exposure. The
945 variance plane of the template image is modified in place.
946 science : `lsst.afw.image.ExposureF`
947 Science exposure to subtract from the template. The variance plane
948 of the science image is modified in place.
949 visitSummary : `lsst.afw.table.ExposureCatalog`, optional
950 Exposure catalog with external calibrations to be applied. Catalog
951 uses the detector id for the catalog id, sorted on id for fast
955 if visitSummary
is not None:
958 template[science.getBBox()], self.
log,
959 requiredTemplateFraction=self.config.requiredTemplateFraction,
960 exceptionMessage=
"Not attempting subtraction. To force subtraction,"
961 " set config requiredTemplateFraction=0"
963 self.metadata[
"templateCoveragePercent"] = 100*templateCoverageFraction
965 if self.config.doScaleVariance:
969 templateVarFactor = self.scaleVariance.run(template.maskedImage)
970 sciVarFactor = self.scaleVariance.run(science.maskedImage)
971 self.
log.info(
"Template variance scaling factor: %.2f", templateVarFactor)
972 self.metadata[
"scaleTemplateVarianceFactor"] = templateVarFactor
973 self.
log.info(
"Science variance scaling factor: %.2f", sciVarFactor)
974 self.metadata[
"scaleScienceVarianceFactor"] = sciVarFactor
981 """Update the science and template mask planes before differencing.
985 template : `lsst.afw.image.Exposure`
986 Template exposure, warped to match the science exposure.
987 The template mask planes will be erased, except for a few specified
989 science : `lsst.afw.image.Exposure`
990 Science exposure to subtract from the template.
991 The DETECTED and DETECTED_NEGATIVE mask planes of the science image
994 self.
_clearMask(science.mask, clearMaskPlanes=[
"DETECTED",
"DETECTED_NEGATIVE"])
1001 clearMaskPlanes = [mp
for mp
in template.mask.getMaskPlaneDict().keys()
1002 if mp
not in self.config.preserveTemplateMask]
1003 renameMaskPlanes = [mp
for mp
in self.config.renameTemplateMask
1004 if mp
in template.mask.getMaskPlaneDict().keys()]
1009 if "FAKE" in science.mask.getMaskPlaneDict().keys():
1010 self.
log.info(
"Adding injected mask plane to science image")
1012 if "FAKE" in template.mask.getMaskPlaneDict().keys():
1013 self.
log.info(
"Adding injected mask plane to template image")
1015 if "INJECTED" in renameMaskPlanes:
1016 renameMaskPlanes.remove(
"INJECTED")
1017 if "INJECTED_TEMPLATE" in clearMaskPlanes:
1018 clearMaskPlanes.remove(
"INJECTED_TEMPLATE")
1020 for maskPlane
in renameMaskPlanes:
1022 self.
_clearMask(template.mask, clearMaskPlanes=clearMaskPlanes)
1026 """Rename a mask plane by adding the new name and copying the data.
1030 mask : `lsst.afw.image.Mask`
1031 The mask image to update in place.
1033 The name of the existing mask plane to copy.
1034 newMaskPlane : `str`
1035 The new name of the mask plane that will be added.
1036 If the mask plane already exists, it will be updated in place.
1038 mask.addMaskPlane(newMaskPlane)
1039 originBitMask = mask.getPlaneBitMask(maskPlane)
1040 destinationBitMask = mask.getPlaneBitMask(newMaskPlane)
1041 mask.array |= ((mask.array & originBitMask) > 0)*destinationBitMask
1044 """Clear the mask plane of an exposure.
1048 mask : `lsst.afw.image.Mask`
1049 The mask plane to erase, which will be modified in place.
1050 clearMaskPlanes : `list` of `str`, optional
1051 Erase the specified mask planes.
1052 If not supplied, the entire mask will be erased.
1054 if clearMaskPlanes
is None:
1055 clearMaskPlanes = list(mask.getMaskPlaneDict().keys())
1057 bitMaskToClear = mask.getPlaneBitMask(clearMaskPlanes)
1058 mask &= ~bitMaskToClear
1062 SubtractScoreOutputConnections):
1067 pipelineConnections=AlardLuptonPreconvolveSubtractConnections):
1072 """Subtract a template from a science image, convolving the science image
1073 before computing the kernel, and also convolving the template before
1076 ConfigClass = AlardLuptonPreconvolveSubtractConfig
1077 _DefaultName =
"alardLuptonPreconvolveSubtract"
1079 def run(self, template, science, sources, visitSummary=None):
1080 """Preconvolve the science image with its own PSF,
1081 convolve the template image with a PSF-matching kernel and subtract
1082 from the preconvolved science image.
1086 template : `lsst.afw.image.ExposureF`
1087 The template image, which has previously been warped to the science
1088 image. The template bbox will be padded by a few pixels compared to
1090 science : `lsst.afw.image.ExposureF`
1091 The science exposure.
1092 sources : `lsst.afw.table.SourceCatalog`
1093 Identified sources on the science exposure. This catalog is used to
1094 select sources in order to perform the AL PSF matching on stamp
1096 visitSummary : `lsst.afw.table.ExposureCatalog`, optional
1097 Exposure catalog with complete external calibrations. Catalog uses
1098 the detector id for the catalog id, sorted on id for fast lookup.
1102 results : `lsst.pipe.base.Struct`
1103 ``scoreExposure`` : `lsst.afw.image.ExposureF`
1104 Result of subtracting the convolved template and science
1105 images. Attached PSF is that of the original science image.
1106 ``matchedTemplate`` : `lsst.afw.image.ExposureF`
1107 Warped and PSF-matched template exposure. Attached PSF is that
1108 of the original science image.
1109 ``matchedScience`` : `lsst.afw.image.ExposureF`
1110 The science exposure after convolving with its own PSF.
1111 Attached PSF is that of the original science image.
1112 ``backgroundModel`` : `lsst.afw.math.Function2D`
1113 Background model that was fit while solving for the
1115 ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
1116 Final kernel used to PSF-match the template to the science
1119 self.
_prepareInputs(template, science, visitSummary=visitSummary)
1122 scienceKernel = science.psf.getKernel()
1124 interpolateBadMaskPlanes=
True)
1125 self.metadata[
"convolvedExposure"] =
"Preconvolution"
1128 subtractResults = self.
runPreconvolve(template, science, matchedScience,
1129 selectSources, scienceKernel)
1132 self.
log.warning(
"Failed to match template. Checking coverage")
1135 self.config.minTemplateFractionForExpectedSuccess,
1136 exceptionMessage=
"Template coverage lower than expected to succeed."
1137 f
" Failure is tolerable: {e}")
1141 return subtractResults
1079 def run(self, template, science, sources, visitSummary=None):
…
1143 def runPreconvolve(self, template, science, matchedScience, selectSources, preConvKernel):
1144 """Convolve the science image with its own PSF, then convolve the
1145 template with a matching kernel and subtract to form the Score
1150 template : `lsst.afw.image.ExposureF`
1151 Template exposure, warped to match the science exposure.
1152 science : `lsst.afw.image.ExposureF`
1153 Science exposure to subtract from the template.
1154 matchedScience : `lsst.afw.image.ExposureF`
1155 The science exposure, convolved with the reflection of its own PSF.
1156 selectSources : `lsst.afw.table.SourceCatalog`
1157 Identified sources on the science exposure. This catalog is used to
1158 select sources in order to perform the AL PSF matching on stamp
1160 preConvKernel : `lsst.afw.math.Kernel`
1161 The reflection of the kernel that was used to preconvolve the
1162 `science` exposure. Must be normalized to sum to 1.
1166 results : `lsst.pipe.base.Struct`
1168 ``scoreExposure`` : `lsst.afw.image.ExposureF`
1169 Result of subtracting the convolved template and science
1170 images. Attached PSF is that of the original science image.
1171 ``matchedTemplate`` : `lsst.afw.image.ExposureF`
1172 Warped and PSF-matched template exposure. Attached PSF is that
1173 of the original science image.
1174 ``matchedScience`` : `lsst.afw.image.ExposureF`
1175 The science exposure after convolving with its own PSF.
1176 Attached PSF is that of the original science image.
1177 ``backgroundModel`` : `lsst.afw.math.Function2D`
1178 Background model that was fit while solving for the
1180 ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
1181 Final kernel used to PSF-match the template to the science
1184 bbox = science.getBBox()
1185 innerBBox = preConvKernel.shrinkBBox(bbox)
1187 kernelSources = self.makeKernel.selectKernelSources(template[innerBBox], matchedScience[innerBBox],
1188 candidateList=selectSources,
1190 kernelResult = self.makeKernel.run(template[innerBBox], matchedScience[innerBBox], kernelSources,
1193 matchedTemplate = self.
_convolveExposure(template, kernelResult.psfMatchingKernel,
1197 interpolateBadMaskPlanes=
True,
1198 photoCalib=science.photoCalib)
1200 backgroundModel=(kernelResult.backgroundModel
1201 if self.config.doSubtractBackground
else None))
1202 correctedScore = self.
finalize(template[bbox], science, score,
1203 kernelResult.psfMatchingKernel,
1204 templateMatched=
True, preConvMode=
True,
1205 preConvKernel=preConvKernel)
1207 return lsst.pipe.base.Struct(scoreExposure=correctedScore,
1208 matchedTemplate=matchedTemplate,
1209 matchedScience=matchedScience,
1210 backgroundModel=kernelResult.backgroundModel,
1211 psfMatchingKernel=kernelResult.psfMatchingKernel,
1212 kernelSources=kernelSources)
1143 def runPreconvolve(self, template, science, matchedScience, selectSources, preConvKernel):
…
1216 exceptionMessage=""):
1217 """Raise NoWorkFound if template coverage < requiredTemplateFraction
1221 templateExposure : `lsst.afw.image.ExposureF`
1222 The template exposure to check
1223 logger : `logging.Logger`
1224 Logger for printing output.
1225 requiredTemplateFraction : `float`, optional
1226 Fraction of pixels of the science image required to have coverage
1228 exceptionMessage : `str`, optional
1229 Message to include in the exception raised if the template coverage
1234 templateCoverageFraction: `float`
1235 Fraction of pixels in the template with data.
1239 lsst.pipe.base.NoWorkFound
1240 Raised if fraction of good pixels, defined as not having NO_DATA
1241 set, is less than the requiredTemplateFraction
1245 pixNoData = np.count_nonzero(templateExposure.mask.array
1246 & templateExposure.mask.getPlaneBitMask(
'NO_DATA'))
1247 pixGood = templateExposure.getBBox().getArea() - pixNoData
1248 templateCoverageFraction = pixGood/templateExposure.getBBox().getArea()
1249 logger.info(
"template has %d good pixels (%.1f%%)", pixGood, 100*templateCoverageFraction)
1251 if templateCoverageFraction < requiredTemplateFraction:
1252 message = (
"Insufficient Template Coverage. (%.1f%% < %.1f%%)" % (
1253 100*templateCoverageFraction,
1254 100*requiredTemplateFraction))
1255 raise lsst.pipe.base.NoWorkFound(message +
" " + exceptionMessage)
1256 return templateCoverageFraction
1260 """Subtract template from science, propagating relevant metadata.
1264 science : `lsst.afw.Exposure`
1265 The input science image.
1266 template : `lsst.afw.Exposure`
1267 The template to subtract from the science image.
1268 backgroundModel : `lsst.afw.MaskedImage`, optional
1269 Differential background model
1273 difference : `lsst.afw.Exposure`
1274 The subtracted image.
1276 difference = science.clone()
1277 if backgroundModel
is not None:
1278 difference.maskedImage -= backgroundModel
1279 difference.maskedImage -= template.maskedImage
1284 """Determine that the PSF of ``exp1`` is not wider than that of ``exp2``.
1288 exp1 : `~lsst.afw.image.Exposure`
1289 Exposure with the reference point spread function (PSF) to evaluate.
1290 exp2 : `~lsst.afw.image.Exposure`
1291 Exposure with a candidate point spread function (PSF) to evaluate.
1292 fwhmExposureBuffer : `float`
1293 Fractional buffer margin to be left out of all sides of the image
1294 during the construction of the grid to compute mean PSF FWHM in an
1295 exposure, if the PSF is not available at its average position.
1296 fwhmExposureGrid : `int`
1297 Grid size to compute the mean FWHM in an exposure, if the PSF is not
1298 available at its average position.
1302 True if ``exp1`` has a PSF that is not wider than that of ``exp2`` in
1306 shape1 = getPsfFwhm(exp1.psf, average=
False)
1307 shape2 = getPsfFwhm(exp2.psf, average=
False)
1309 shape1 = evaluateMeanPsfFwhm(exp1,
1310 fwhmExposureBuffer=fwhmExposureBuffer,
1311 fwhmExposureGrid=fwhmExposureGrid
1313 shape2 = evaluateMeanPsfFwhm(exp2,
1314 fwhmExposureBuffer=fwhmExposureBuffer,
1315 fwhmExposureGrid=fwhmExposureGrid
1317 return shape1 <= shape2
1320 xTest = shape1[0] <= shape2[0]
1321 yTest = shape1[1] <= shape2[1]
1322 return xTest | yTest
1326 """Replace masked image pixels with interpolated values.
1330 maskedImage : `lsst.afw.image.MaskedImage`
1331 Image on which to perform interpolation.
1332 badMaskPlanes : `list` of `str`
1333 List of mask planes to interpolate over.
1334 fallbackValue : `float`, optional
1335 Value to set when interpolation fails.
1340 The number of masked pixels that were replaced.
1342 imgBadMaskPlanes = [
1343 maskPlane
for maskPlane
in badMaskPlanes
if maskPlane
in maskedImage.mask.getMaskPlaneDict()
1346 image = maskedImage.image.array
1347 badPixels = (maskedImage.mask.array & maskedImage.mask.getPlaneBitMask(imgBadMaskPlanes)) > 0
1348 image[badPixels] = np.nan
1349 if fallbackValue
is None:
1350 fallbackValue = np.nanmedian(image)
1353 image[badPixels] = fallbackValue
1354 return np.sum(badPixels)
Parameters to control convolution.
runPreconvolve(self, template, science, matchedScience, selectSources, preConvKernel)
run(self, template, science, sources, visitSummary=None)
_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)