22from astropy
import units
as u
29from lsst.ip.diffim.utils
import evaluateMeanPsfFwhm, getPsfFwhm
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",
106 dimensions=_dimensions,
107 defaultTemplates=_defaultTemplates):
108 scoreExposure = connectionTypes.Output(
109 doc=
"The maximum likelihood image, used for the detection of diaSources.",
110 dimensions=(
"instrument",
"visit",
"detector"),
111 storageClass=
"ExposureF",
112 name=
"{fakesType}{coaddName}Diff_scoreExp",
114 psfMatchingKernel = connectionTypes.Output(
115 doc=
"Kernel used to PSF match the science and template images.",
116 dimensions=(
"instrument",
"visit",
"detector"),
117 storageClass=
"MatchingKernel",
118 name=
"{fakesType}{coaddName}Diff_psfScoreMatchKernel",
128 target=MakeKernelTask,
129 doc=
"Task to construct a matching kernel for convolution.",
134 doc=
"Perform diffim decorrelation to undo pixel correlation due to A&L "
135 "kernel convolution? If True, also update the diffim PSF."
138 target=DecorrelateALKernelTask,
139 doc=
"Task to decorrelate the image difference.",
144 doc=
"Raise NoWorkFound and do not attempt image subtraction if template covers less than this "
145 " fraction of pixels. Setting to 0 will always attempt image subtraction."
150 doc=
"Raise NoWorkFound if PSF-matching fails and template covers less than this fraction of pixels."
151 " If the fraction of pixels covered by the template is less than this value (and greater than"
152 " requiredTemplateFraction) this task is attempted but failure is anticipated and tolerated."
157 doc=
"Scale variance of the image difference?"
160 target=ScaleVarianceTask,
161 doc=
"Subtask to rescale the variance of the template to the statistically expected level."
164 doc=
"Subtract the background fit when solving the kernel?",
170 "Replace science Exposure's calibration objects with those"
171 " in visitSummary. Ignored if `doApplyFinalizedPsf is True."
177 target=ScienceSourceSelectorTask,
178 doc=
"Task to select sources to be used for PSF matching.",
183 doc=
"Minimum signal to noise ratio of detected sources "
184 "to use for calculating the PSF matching kernel."
189 doc=
"Maximum signal to noise ratio of detected sources "
190 "to use for calculating the PSF matching kernel."
195 doc=
"Maximum number of sources to use for calculating the PSF matching kernel."
196 "Set to -1 to disable."
201 doc=
"Minimum number of sources needed for calculating the PSF matching kernel."
205 doc=
"Flags that, if set, the associated source should not "
206 "be used to determine the PSF matching kernel.",
207 default=(
"sky_source",
"slot_Centroid_flag",
208 "slot_ApFlux_flag",
"slot_PsfFlux_flag",
209 "base_PixelFlags_flag_interpolated",
210 "base_PixelFlags_flag_saturated",
211 "base_PixelFlags_flag_bad",
213 deprecated=
"This field is no longer used and will be removed after v28."
214 "Set the equivalent field for the sourceSelector subtask instead.",
218 default=(
"NO_DATA",
"BAD",
"SAT",
"EDGE",
"FAKE"),
219 doc=
"Mask planes to exclude when selecting sources for PSF matching."
223 default=(
"NO_DATA",
"BAD",
"SAT",
"EDGE"),
224 doc=
"Mask planes to interpolate over."
228 default=(
"NO_DATA",
"BAD",),
229 doc=
"Mask planes from the template to propagate to the image difference."
233 default=(
"SAT",
"INJECTED",
"INJECTED_CORE",),
234 doc=
"Mask planes from the template to propagate to the image difference"
235 "with '_TEMPLATE' appended to the name."
240 doc=
"Re-run source detection for kernel candidates if an error is"
241 " encountered while calculating the matching kernel."
247 self.
makeKernel.kernel.active.spatialKernelOrder = 1
248 self.
makeKernel.kernel.active.spatialBgOrder = 2
259 pipelineConnections=AlardLuptonSubtractConnections):
262 default=
"convolveTemplate",
263 allowed={
"auto":
"Choose which image to convolve at runtime.",
264 "convolveScience":
"Only convolve the science image.",
265 "convolveTemplate":
"Only convolve the template image."},
266 doc=
"Choose which image to convolve at runtime, or require that a specific image is convolved."
271 """Compute the image difference of a science and template image using
272 the Alard & Lupton (1998) algorithm.
274 ConfigClass = AlardLuptonSubtractConfig
275 _DefaultName =
"alardLuptonSubtract"
279 self.makeSubtask(
"decorrelate")
280 self.makeSubtask(
"makeKernel")
281 self.makeSubtask(
"sourceSelector")
282 if self.config.doScaleVariance:
283 self.makeSubtask(
"scaleVariance")
292 """Replace calibrations (psf, and ApCorrMap) on this exposure with
297 exposure : `lsst.afw.image.exposure.Exposure`
298 Input exposure to adjust calibrations.
299 visitSummary : `lsst.afw.table.ExposureCatalog`
300 Exposure catalog with external calibrations to be applied. Catalog
301 uses the detector id for the catalog id, sorted on id for fast
306 exposure : `lsst.afw.image.exposure.Exposure`
307 Exposure with adjusted calibrations.
309 detectorId = exposure.info.getDetector().getId()
311 row = visitSummary.find(detectorId)
313 self.
log.warning(
"Detector id %s not found in external calibrations catalog; "
314 "Using original calibrations.", detectorId)
317 apCorrMap = row.getApCorrMap()
319 self.
log.warning(
"Detector id %s has None for psf in "
320 "external calibrations catalog; Using original psf and aperture correction.",
322 elif apCorrMap
is None:
323 self.
log.warning(
"Detector id %s has None for apCorrMap in "
324 "external calibrations catalog; Using original psf and aperture correction.",
328 exposure.info.setApCorrMap(apCorrMap)
333 def run(self, template, science, sources, visitSummary=None):
334 """PSF match, subtract, and decorrelate two images.
338 template : `lsst.afw.image.ExposureF`
339 Template exposure, warped to match the science exposure.
340 science : `lsst.afw.image.ExposureF`
341 Science exposure to subtract from the template.
342 sources : `lsst.afw.table.SourceCatalog`
343 Identified sources on the science exposure. This catalog is used to
344 select sources in order to perform the AL PSF matching on stamp
346 visitSummary : `lsst.afw.table.ExposureCatalog`, optional
347 Exposure catalog with external calibrations to be applied. Catalog
348 uses the detector id for the catalog id, sorted on id for fast
353 results : `lsst.pipe.base.Struct`
354 ``difference`` : `lsst.afw.image.ExposureF`
355 Result of subtracting template and science.
356 ``matchedTemplate`` : `lsst.afw.image.ExposureF`
357 Warped and PSF-matched template exposure.
358 ``backgroundModel`` : `lsst.afw.math.Function2D`
359 Background model that was fit while solving for the
361 ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
362 Kernel used to PSF-match the convolved image.
367 If an unsupported convolution mode is supplied.
369 If there are too few sources to calculate the PSF matching kernel.
370 lsst.pipe.base.NoWorkFound
371 Raised if fraction of good pixels, defined as not having NO_DATA
372 set, is less then the configured requiredTemplateFraction
378 fwhmExposureBuffer = self.config.makeKernel.fwhmExposureBuffer
379 fwhmExposureGrid = self.config.makeKernel.fwhmExposureGrid
393 templatePsfSize = getPsfFwhm(template.psf)
394 sciencePsfSize = getPsfFwhm(science.psf)
396 self.
log.info(
"Unable to evaluate PSF at the average position. "
397 "Evaluting PSF on a grid of points."
399 templatePsfSize = evaluateMeanPsfFwhm(template,
400 fwhmExposureBuffer=fwhmExposureBuffer,
401 fwhmExposureGrid=fwhmExposureGrid
403 sciencePsfSize = evaluateMeanPsfFwhm(science,
404 fwhmExposureBuffer=fwhmExposureBuffer,
405 fwhmExposureGrid=fwhmExposureGrid
407 self.
log.info(
"Science PSF FWHM: %f pixels", sciencePsfSize)
408 self.
log.info(
"Template PSF FWHM: %f pixels", templatePsfSize)
409 self.metadata.add(
"sciencePsfSize", sciencePsfSize)
410 self.metadata.add(
"templatePsfSize", templatePsfSize)
413 maglim_science = self.
_calculateMagLim(science, fallbackPsfSize=sciencePsfSize)
414 fluxlim_science = (maglim_science*u.ABmag).to_value(u.nJy)
415 maglim_template = self.
_calculateMagLim(template, fallbackPsfSize=templatePsfSize)
416 if np.isnan(maglim_template):
417 self.
log.info(
"Cannot evaluate template limiting mag; adopting science limiting mag for diffim")
418 maglim_diffim = maglim_science
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.add(
"scienceLimitingMagnitude", maglim_science)
423 self.metadata.add(
"templateLimitingMagnitude", maglim_template)
424 self.metadata.add(
"diffimLimitingMagnitude", maglim_diffim)
426 if self.config.mode ==
"auto":
429 fwhmExposureBuffer=fwhmExposureBuffer,
430 fwhmExposureGrid=fwhmExposureGrid)
432 if sciencePsfSize < templatePsfSize:
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.add(
"convolvedExposure",
"Template")
456 self.metadata.add(
"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}")
469 return subtractResults
333 def run(self, template, science, sources, visitSummary=None):
…
472 """Convolve the template image with a PSF-matching kernel and subtract
473 from the science image.
477 template : `lsst.afw.image.ExposureF`
478 Template exposure, warped to match the science exposure.
479 science : `lsst.afw.image.ExposureF`
480 Science exposure to subtract from the template.
481 selectSources : `lsst.afw.table.SourceCatalog`
482 Identified sources on the science exposure. This catalog is used to
483 select sources in order to perform the AL PSF matching on stamp
488 results : `lsst.pipe.base.Struct`
490 ``difference`` : `lsst.afw.image.ExposureF`
491 Result of subtracting template and science.
492 ``matchedTemplate`` : `lsst.afw.image.ExposureF`
493 Warped and PSF-matched template exposure.
494 ``backgroundModel`` : `lsst.afw.math.Function2D`
495 Background model that was fit while solving for the PSF-matching kernel
496 ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
497 Kernel used to PSF-match the template to the science image.
500 kernelSources = self.makeKernel.selectKernelSources(template, science,
501 candidateList=selectSources,
503 kernelResult = self.makeKernel.run(template, science, kernelSources,
505 except Exception
as e:
506 if self.config.allowKernelSourceDetection:
507 self.
log.warning(
"Error encountered trying to construct the matching kernel"
508 f
" Running source detection and retrying. {e}")
509 kernelSources = self.makeKernel.selectKernelSources(template, science,
512 kernelResult = self.makeKernel.run(template, science, kernelSources,
517 matchedTemplate = self.
_convolveExposure(template, kernelResult.psfMatchingKernel,
519 bbox=science.getBBox(),
521 photoCalib=science.photoCalib)
524 backgroundModel=(kernelResult.backgroundModel
525 if self.config.doSubtractBackground
else None))
526 correctedExposure = self.
finalize(template, science, difference,
527 kernelResult.psfMatchingKernel,
528 templateMatched=
True)
530 return lsst.pipe.base.Struct(difference=correctedExposure,
531 matchedTemplate=matchedTemplate,
532 matchedScience=science,
533 backgroundModel=kernelResult.backgroundModel,
534 psfMatchingKernel=kernelResult.psfMatchingKernel)
537 """Convolve the science image with a PSF-matching kernel and subtract
542 template : `lsst.afw.image.ExposureF`
543 Template exposure, warped to match the science exposure.
544 science : `lsst.afw.image.ExposureF`
545 Science exposure to subtract from the template.
546 selectSources : `lsst.afw.table.SourceCatalog`
547 Identified sources on the science exposure. This catalog is used to
548 select sources in order to perform the AL PSF matching on stamp
553 results : `lsst.pipe.base.Struct`
555 ``difference`` : `lsst.afw.image.ExposureF`
556 Result of subtracting template and science.
557 ``matchedTemplate`` : `lsst.afw.image.ExposureF`
558 Warped template exposure. Note that in this case, the template
559 is not PSF-matched to the science image.
560 ``backgroundModel`` : `lsst.afw.math.Function2D`
561 Background model that was fit while solving for the PSF-matching kernel
562 ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
563 Kernel used to PSF-match the science image to the template.
565 bbox = science.getBBox()
566 kernelSources = self.makeKernel.selectKernelSources(science, template,
567 candidateList=selectSources,
569 kernelResult = self.makeKernel.run(science, template, kernelSources,
571 modelParams = kernelResult.backgroundModel.getParameters()
573 kernelResult.backgroundModel.setParameters([-p
for p
in modelParams])
575 kernelImage = lsst.afw.image.ImageD(kernelResult.psfMatchingKernel.getDimensions())
576 norm = kernelResult.psfMatchingKernel.computeImage(kernelImage, doNormalize=
False)
583 matchedScience.maskedImage /= norm
584 matchedTemplate = template.clone()[bbox]
585 matchedTemplate.maskedImage /= norm
586 matchedTemplate.setPhotoCalib(science.photoCalib)
589 backgroundModel=(kernelResult.backgroundModel
590 if self.config.doSubtractBackground
else None))
592 correctedExposure = self.
finalize(template, science, difference,
593 kernelResult.psfMatchingKernel,
594 templateMatched=
False)
596 return lsst.pipe.base.Struct(difference=correctedExposure,
597 matchedTemplate=matchedTemplate,
598 matchedScience=matchedScience,
599 backgroundModel=kernelResult.backgroundModel,
600 psfMatchingKernel=kernelResult.psfMatchingKernel,)
602 def finalize(self, template, science, difference, kernel,
603 templateMatched=True,
606 spatiallyVarying=False):
607 """Decorrelate the difference image to undo the noise correlations
608 caused by convolution.
612 template : `lsst.afw.image.ExposureF`
613 Template exposure, warped to match the science exposure.
614 science : `lsst.afw.image.ExposureF`
615 Science exposure to subtract from the template.
616 difference : `lsst.afw.image.ExposureF`
617 Result of subtracting template and science.
618 kernel : `lsst.afw.math.Kernel`
619 An (optionally spatially-varying) PSF matching kernel
620 templateMatched : `bool`, optional
621 Was the template PSF-matched to the science image?
622 preConvMode : `bool`, optional
623 Was the science image preconvolved with its own PSF
624 before PSF matching the template?
625 preConvKernel : `lsst.afw.detection.Psf`, optional
626 If not `None`, then the science image was pre-convolved with
627 (the reflection of) this kernel. Must be normalized to sum to 1.
628 spatiallyVarying : `bool`, optional
629 Compute the decorrelation kernel spatially varying across the image?
633 correctedExposure : `lsst.afw.image.ExposureF`
634 The decorrelated image difference.
636 if self.config.doDecorrelation:
637 self.
log.info(
"Decorrelating image difference.")
641 correctedExposure = self.decorrelate.run(science, template[science.getBBox()], difference, kernel,
642 templateMatched=templateMatched,
643 preConvMode=preConvMode,
644 preConvKernel=preConvKernel,
645 spatiallyVarying=spatiallyVarying).correctedExposure
647 self.
log.info(
"NOT decorrelating image difference.")
648 correctedExposure = difference
649 return correctedExposure
602 def finalize(self, template, science, difference, kernel,
…
652 """Calculate an exposure's limiting magnitude.
654 This method uses the photometric zeropoint together with the
655 PSF size from the average position of the exposure.
659 exposure : `lsst.afw.image.Exposure`
660 The target exposure to calculate the limiting magnitude for.
661 nsigma : `float`, optional
662 The detection threshold in sigma.
663 fallbackPsfSize : `float`, optional
664 PSF FWHM to use in the event the exposure PSF cannot be retrieved.
668 maglim : `astropy.units.Quantity`
669 The limiting magnitude of the exposure, or np.nan.
672 psf = exposure.getPsf()
673 psf_shape = psf.computeShape(psf.getAveragePosition())
675 if fallbackPsfSize
is not None:
676 self.
log.info(
"Unable to evaluate PSF, using fallback FWHM %f", fallbackPsfSize)
677 psf_area = np.pi*(fallbackPsfSize/2)**2
678 zeropoint = exposure.getPhotoCalib().instFluxToMagnitude(1)
679 maglim = zeropoint - 2.5*np.log10(nsigma*np.sqrt(psf_area))
681 self.
log.info(
"Unable to evaluate PSF, setting maglim to nan")
685 psf_area = np.pi*np.sqrt(psf_shape.getIxx()*psf_shape.getIyy())
686 zeropoint = exposure.getPhotoCalib().instFluxToMagnitude(1)
687 maglim = zeropoint - 2.5*np.log10(nsigma*np.sqrt(psf_area))
693 """Check that the WCS of the two Exposures match, the template bbox
694 contains the science bbox, and that the bands match.
698 template : `lsst.afw.image.ExposureF`
699 Template exposure, warped to match the science exposure.
700 science : `lsst.afw.image.ExposureF`
701 Science exposure to subtract from the template.
706 Raised if the WCS of the template is not equal to the science WCS,
707 if the science image is not fully contained in the template
708 bounding box, or if the bands do not match.
710 assert template.wcs == science.wcs, \
711 "Template and science exposure WCS are not identical."
712 templateBBox = template.getBBox()
713 scienceBBox = science.getBBox()
714 assert science.filter.bandLabel == template.filter.bandLabel, \
715 "Science and template exposures have different bands: %s, %s" % \
716 (science.filter, template.filter)
718 assert templateBBox.contains(scienceBBox), \
719 "Template bbox does not contain all of the science image."
725 interpolateBadMaskPlanes=False,
727 """Convolve an exposure with the given kernel.
731 exposure : `lsst.afw.Exposure`
732 exposure to convolve.
733 kernel : `lsst.afw.math.LinearCombinationKernel`
734 PSF matching kernel computed in the ``makeKernel`` subtask.
735 convolutionControl : `lsst.afw.math.ConvolutionControl`
736 Configuration for convolve algorithm.
737 bbox : `lsst.geom.Box2I`, optional
738 Bounding box to trim the convolved exposure to.
739 psf : `lsst.afw.detection.Psf`, optional
740 Point spread function (PSF) to set for the convolved exposure.
741 photoCalib : `lsst.afw.image.PhotoCalib`, optional
742 Photometric calibration of the convolved exposure.
746 convolvedExp : `lsst.afw.Exposure`
749 convolvedExposure = exposure.clone()
751 convolvedExposure.setPsf(psf)
752 if photoCalib
is not None:
753 convolvedExposure.setPhotoCalib(photoCalib)
754 if interpolateBadMaskPlanes
and self.config.badMaskPlanes
is not None:
756 self.config.badMaskPlanes)
757 self.metadata.add(
"nInterpolated", nInterp)
758 convolvedImage = lsst.afw.image.MaskedImageF(convolvedExposure.getBBox())
760 convolvedExposure.setMaskedImage(convolvedImage)
762 return convolvedExposure
764 return convolvedExposure[bbox]
767 """Select sources from a catalog that meet the selection criteria.
771 sources : `lsst.afw.table.SourceCatalog`
772 Input source catalog to select sources from.
773 mask : `lsst.afw.image.Mask`
774 The image mask plane to use to reject sources
775 based on their location on the ccd.
779 selectSources : `lsst.afw.table.SourceCatalog`
780 The input source catalog, with flagged and low signal-to-noise
786 If there are too few sources to compute the PSF matching kernel
787 remaining after source selection.
790 selected = self.sourceSelector.selectSources(sources).selected
791 nInitialSelected = np.count_nonzero(selected)
792 selected *= self.
_checkMask(mask, sources, self.config.excludeMaskPlanes)
793 nSelected = np.count_nonzero(selected)
794 self.
log.info(
"Rejecting %i candidate sources: an excluded template mask plane is set.",
795 nInitialSelected - nSelected)
796 selectSources = sources[selected].copy(deep=
True)
799 if (len(selectSources) > self.config.maxKernelSources) & (self.config.maxKernelSources > 0):
800 signalToNoise = selectSources.getPsfInstFlux()/selectSources.getPsfInstFluxErr()
801 indices = np.argsort(signalToNoise)
802 indices = indices[-self.config.maxKernelSources:]
803 selected = np.zeros(len(selectSources), dtype=bool)
804 selected[indices] =
True
805 selectSources = selectSources[selected].copy(deep=
True)
807 self.
log.info(
"%i/%i=%.1f%% of sources selected for PSF matching from the input catalog",
808 len(selectSources), len(sources), 100*len(selectSources)/len(sources))
809 if len(selectSources) < self.config.minKernelSources:
810 self.
log.error(
"Too few sources to calculate the PSF matching kernel: "
811 "%i selected but %i needed for the calculation.",
812 len(selectSources), self.config.minKernelSources)
813 raise RuntimeError(
"Cannot compute PSF matching kernel: too few sources selected.")
814 self.metadata.add(
"nPsfSources", len(selectSources))
819 def _checkMask(mask, sources, excludeMaskPlanes, checkAdjacent=True):
820 """Exclude sources that are located on or adjacent to masked pixels.
824 mask : `lsst.afw.image.Mask`
825 The image mask plane to use to reject sources
826 based on the location of their centroid on the ccd.
827 sources : `lsst.afw.table.SourceCatalog`
828 The source catalog to evaluate.
829 excludeMaskPlanes : `list` of `str`
830 List of the names of the mask planes to exclude.
834 flags : `numpy.ndarray` of `bool`
835 Array indicating whether each source in the catalog should be
836 kept (True) or rejected (False) based on the value of the
837 mask plane at its location.
839 setExcludeMaskPlanes = [
840 maskPlane
for maskPlane
in excludeMaskPlanes
if maskPlane
in mask.getMaskPlaneDict()
843 excludePixelMask = mask.getPlaneBitMask(setExcludeMaskPlanes)
845 xv = (np.rint(sources.getX() - mask.getX0())).astype(int)
846 yv = (np.rint(sources.getY() - mask.getY0())).astype(int)
848 flags = np.ones(len(sources), dtype=bool)
850 pixRange = (0, -1, 1)
855 mv = mask.array[yv + j, xv + i]
856 flags *= np.bitwise_and(mv, excludePixelMask) == 0
819 def _checkMask(mask, sources, excludeMaskPlanes, checkAdjacent=True):
…
860 """Perform preparatory calculations common to all Alard&Lupton Tasks.
864 template : `lsst.afw.image.ExposureF`
865 Template exposure, warped to match the science exposure. The
866 variance plane of the template image is modified in place.
867 science : `lsst.afw.image.ExposureF`
868 Science exposure to subtract from the template. The variance plane
869 of the science image is modified in place.
870 visitSummary : `lsst.afw.table.ExposureCatalog`, optional
871 Exposure catalog with external calibrations to be applied. Catalog
872 uses the detector id for the catalog id, sorted on id for fast
876 if visitSummary
is not None:
879 template[science.getBBox()], self.
log,
880 requiredTemplateFraction=self.config.requiredTemplateFraction,
881 exceptionMessage=
"Not attempting subtraction. To force subtraction,"
882 " set config requiredTemplateFraction=0"
884 self.metadata.add(
"templateCoveragePercent", 100*templateCoverageFraction)
886 if self.config.doScaleVariance:
890 templateVarFactor = self.scaleVariance.run(template.maskedImage)
891 sciVarFactor = self.scaleVariance.run(science.maskedImage)
892 self.
log.info(
"Template variance scaling factor: %.2f", templateVarFactor)
893 self.metadata.add(
"scaleTemplateVarianceFactor", templateVarFactor)
894 self.
log.info(
"Science variance scaling factor: %.2f", sciVarFactor)
895 self.metadata.add(
"scaleScienceVarianceFactor", sciVarFactor)
902 """Update the science and template mask planes before differencing.
906 template : `lsst.afw.image.Exposure`
907 Template exposure, warped to match the science exposure.
908 The template mask planes will be erased, except for a few specified
910 science : `lsst.afw.image.Exposure`
911 Science exposure to subtract from the template.
912 The DETECTED and DETECTED_NEGATIVE mask planes of the science image
915 self.
_clearMask(science.mask, clearMaskPlanes=[
"DETECTED",
"DETECTED_NEGATIVE"])
922 clearMaskPlanes = [mp
for mp
in template.mask.getMaskPlaneDict().keys()
923 if mp
not in self.config.preserveTemplateMask]
924 renameMaskPlanes = [mp
for mp
in self.config.renameTemplateMask
925 if mp
in template.mask.getMaskPlaneDict().keys()]
930 if "FAKE" in science.mask.getMaskPlaneDict().keys():
931 self.
log.info(
"Adding injected mask plane to science image")
933 if "FAKE" in template.mask.getMaskPlaneDict().keys():
934 self.
log.info(
"Adding injected mask plane to template image")
936 if "INJECTED" in renameMaskPlanes:
937 renameMaskPlanes.remove(
"INJECTED")
938 if "INJECTED_TEMPLATE" in clearMaskPlanes:
939 clearMaskPlanes.remove(
"INJECTED_TEMPLATE")
941 for maskPlane
in renameMaskPlanes:
943 self.
_clearMask(template.mask, clearMaskPlanes=clearMaskPlanes)
947 """Rename a mask plane by adding the new name and copying the data.
951 mask : `lsst.afw.image.Mask`
952 The mask image to update in place.
954 The name of the existing mask plane to copy.
956 The new name of the mask plane that will be added.
957 If the mask plane already exists, it will be updated in place.
959 mask.addMaskPlane(newMaskPlane)
960 originBitMask = mask.getPlaneBitMask(maskPlane)
961 destinationBitMask = mask.getPlaneBitMask(newMaskPlane)
962 mask.array |= ((mask.array & originBitMask) > 0)*destinationBitMask
965 """Clear the mask plane of an exposure.
969 mask : `lsst.afw.image.Mask`
970 The mask plane to erase, which will be modified in place.
971 clearMaskPlanes : `list` of `str`, optional
972 Erase the specified mask planes.
973 If not supplied, the entire mask will be erased.
975 if clearMaskPlanes
is None:
976 clearMaskPlanes = list(mask.getMaskPlaneDict().keys())
978 bitMaskToClear = mask.getPlaneBitMask(clearMaskPlanes)
979 mask &= ~bitMaskToClear
983 SubtractScoreOutputConnections):
988 pipelineConnections=AlardLuptonPreconvolveSubtractConnections):
993 """Subtract a template from a science image, convolving the science image
994 before computing the kernel, and also convolving the template before
997 ConfigClass = AlardLuptonPreconvolveSubtractConfig
998 _DefaultName =
"alardLuptonPreconvolveSubtract"
1000 def run(self, template, science, sources, visitSummary=None):
1001 """Preconvolve the science image with its own PSF,
1002 convolve the template image with a PSF-matching kernel and subtract
1003 from the preconvolved science image.
1007 template : `lsst.afw.image.ExposureF`
1008 The template image, which has previously been warped to the science
1009 image. The template bbox will be padded by a few pixels compared to
1011 science : `lsst.afw.image.ExposureF`
1012 The science exposure.
1013 sources : `lsst.afw.table.SourceCatalog`
1014 Identified sources on the science exposure. This catalog is used to
1015 select sources in order to perform the AL PSF matching on stamp
1017 visitSummary : `lsst.afw.table.ExposureCatalog`, optional
1018 Exposure catalog with complete external calibrations. Catalog uses
1019 the detector id for the catalog id, sorted on id for fast lookup.
1023 results : `lsst.pipe.base.Struct`
1024 ``scoreExposure`` : `lsst.afw.image.ExposureF`
1025 Result of subtracting the convolved template and science
1026 images. Attached PSF is that of the original science image.
1027 ``matchedTemplate`` : `lsst.afw.image.ExposureF`
1028 Warped and PSF-matched template exposure. Attached PSF is that
1029 of the original science image.
1030 ``matchedScience`` : `lsst.afw.image.ExposureF`
1031 The science exposure after convolving with its own PSF.
1032 Attached PSF is that of the original science image.
1033 ``backgroundModel`` : `lsst.afw.math.Function2D`
1034 Background model that was fit while solving for the
1036 ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
1037 Final kernel used to PSF-match the template to the science
1040 self.
_prepareInputs(template, science, visitSummary=visitSummary)
1043 scienceKernel = science.psf.getKernel()
1045 interpolateBadMaskPlanes=
True)
1046 self.metadata.add(
"convolvedExposure",
"Preconvolution")
1049 subtractResults = self.
runPreconvolve(template, science, matchedScience,
1050 selectSources, scienceKernel)
1053 self.
loglog.warning(
"Failed to match template. Checking coverage")
1056 self.config.minTemplateFractionForExpectedSuccess,
1057 exceptionMessage=
"Template coverage lower than expected to succeed."
1058 f
" Failure is tolerable: {e}")
1062 return subtractResults
1000 def run(self, template, science, sources, visitSummary=None):
…
1064 def runPreconvolve(self, template, science, matchedScience, selectSources, preConvKernel):
1065 """Convolve the science image with its own PSF, then convolve the
1066 template with a matching kernel and subtract to form the Score
1071 template : `lsst.afw.image.ExposureF`
1072 Template exposure, warped to match the science exposure.
1073 science : `lsst.afw.image.ExposureF`
1074 Science exposure to subtract from the template.
1075 matchedScience : `lsst.afw.image.ExposureF`
1076 The science exposure, convolved with the reflection of its own PSF.
1077 selectSources : `lsst.afw.table.SourceCatalog`
1078 Identified sources on the science exposure. This catalog is used to
1079 select sources in order to perform the AL PSF matching on stamp
1081 preConvKernel : `lsst.afw.math.Kernel`
1082 The reflection of the kernel that was used to preconvolve the
1083 `science` exposure. Must be normalized to sum to 1.
1087 results : `lsst.pipe.base.Struct`
1089 ``scoreExposure`` : `lsst.afw.image.ExposureF`
1090 Result of subtracting the convolved template and science
1091 images. Attached PSF is that of the original science image.
1092 ``matchedTemplate`` : `lsst.afw.image.ExposureF`
1093 Warped and PSF-matched template exposure. Attached PSF is that
1094 of the original science image.
1095 ``matchedScience`` : `lsst.afw.image.ExposureF`
1096 The science exposure after convolving with its own PSF.
1097 Attached PSF is that of the original science image.
1098 ``backgroundModel`` : `lsst.afw.math.Function2D`
1099 Background model that was fit while solving for the
1101 ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
1102 Final kernel used to PSF-match the template to the science
1105 bbox = science.getBBox()
1106 innerBBox = preConvKernel.shrinkBBox(bbox)
1108 kernelSources = self.makeKernel.selectKernelSources(template[innerBBox], matchedScience[innerBBox],
1109 candidateList=selectSources,
1111 kernelResult = self.makeKernel.run(template[innerBBox], matchedScience[innerBBox], kernelSources,
1114 matchedTemplate = self.
_convolveExposure(template, kernelResult.psfMatchingKernel,
1118 interpolateBadMaskPlanes=
True,
1119 photoCalib=science.photoCalib)
1121 backgroundModel=(kernelResult.backgroundModel
1122 if self.config.doSubtractBackground
else None))
1123 correctedScore = self.
finalize(template[bbox], science, score,
1124 kernelResult.psfMatchingKernel,
1125 templateMatched=
True, preConvMode=
True,
1126 preConvKernel=preConvKernel)
1128 return lsst.pipe.base.Struct(scoreExposure=correctedScore,
1129 matchedTemplate=matchedTemplate,
1130 matchedScience=matchedScience,
1131 backgroundModel=kernelResult.backgroundModel,
1132 psfMatchingKernel=kernelResult.psfMatchingKernel)
1064 def runPreconvolve(self, template, science, matchedScience, selectSources, preConvKernel):
…
1136 exceptionMessage=""):
1137 """Raise NoWorkFound if template coverage < requiredTemplateFraction
1141 templateExposure : `lsst.afw.image.ExposureF`
1142 The template exposure to check
1143 logger : `logging.Logger`
1144 Logger for printing output.
1145 requiredTemplateFraction : `float`, optional
1146 Fraction of pixels of the science image required to have coverage
1148 exceptionMessage : `str`, optional
1149 Message to include in the exception raised if the template coverage
1154 templateCoverageFraction: `float`
1155 Fraction of pixels in the template with data.
1159 lsst.pipe.base.NoWorkFound
1160 Raised if fraction of good pixels, defined as not having NO_DATA
1161 set, is less than the requiredTemplateFraction
1165 pixNoData = np.count_nonzero(templateExposure.mask.array
1166 & templateExposure.mask.getPlaneBitMask(
'NO_DATA'))
1167 pixGood = templateExposure.getBBox().getArea() - pixNoData
1168 templateCoverageFraction = pixGood/templateExposure.getBBox().getArea()
1169 logger.info(
"template has %d good pixels (%.1f%%)", pixGood, 100*templateCoverageFraction)
1171 if templateCoverageFraction < requiredTemplateFraction:
1172 message = (
"Insufficient Template Coverage. (%.1f%% < %.1f%%)" % (
1173 100*templateCoverageFraction,
1174 100*requiredTemplateFraction))
1175 raise lsst.pipe.base.NoWorkFound(message +
" " + exceptionMessage)
1176 return templateCoverageFraction
1180 """Subtract template from science, propagating relevant metadata.
1184 science : `lsst.afw.Exposure`
1185 The input science image.
1186 template : `lsst.afw.Exposure`
1187 The template to subtract from the science image.
1188 backgroundModel : `lsst.afw.MaskedImage`, optional
1189 Differential background model
1193 difference : `lsst.afw.Exposure`
1194 The subtracted image.
1196 difference = science.clone()
1197 if backgroundModel
is not None:
1198 difference.maskedImage -= backgroundModel
1199 difference.maskedImage -= template.maskedImage
1204 """Determine that the PSF of ``exp1`` is not wider than that of ``exp2``.
1208 exp1 : `~lsst.afw.image.Exposure`
1209 Exposure with the reference point spread function (PSF) to evaluate.
1210 exp2 : `~lsst.afw.image.Exposure`
1211 Exposure with a candidate point spread function (PSF) to evaluate.
1212 fwhmExposureBuffer : `float`
1213 Fractional buffer margin to be left out of all sides of the image
1214 during the construction of the grid to compute mean PSF FWHM in an
1215 exposure, if the PSF is not available at its average position.
1216 fwhmExposureGrid : `int`
1217 Grid size to compute the mean FWHM in an exposure, if the PSF is not
1218 available at its average position.
1222 True if ``exp1`` has a PSF that is not wider than that of ``exp2`` in
1226 shape1 = getPsfFwhm(exp1.psf, average=
False)
1227 shape2 = getPsfFwhm(exp2.psf, average=
False)
1229 shape1 = evaluateMeanPsfFwhm(exp1,
1230 fwhmExposureBuffer=fwhmExposureBuffer,
1231 fwhmExposureGrid=fwhmExposureGrid
1233 shape2 = evaluateMeanPsfFwhm(exp2,
1234 fwhmExposureBuffer=fwhmExposureBuffer,
1235 fwhmExposureGrid=fwhmExposureGrid
1237 return shape1 <= shape2
1240 xTest = shape1[0] <= shape2[0]
1241 yTest = shape1[1] <= shape2[1]
1242 return xTest | yTest
1246 """Replace masked image pixels with interpolated values.
1250 maskedImage : `lsst.afw.image.MaskedImage`
1251 Image on which to perform interpolation.
1252 badMaskPlanes : `list` of `str`
1253 List of mask planes to interpolate over.
1254 fallbackValue : `float`, optional
1255 Value to set when interpolation fails.
1260 The number of masked pixels that were replaced.
1262 imgBadMaskPlanes = [
1263 maskPlane
for maskPlane
in badMaskPlanes
if maskPlane
in maskedImage.mask.getMaskPlaneDict()
1266 image = maskedImage.image.array
1267 badPixels = (maskedImage.mask.array & maskedImage.mask.getPlaneBitMask(imgBadMaskPlanes)) > 0
1268 image[badPixels] = np.nan
1269 if fallbackValue
is None:
1270 fallbackValue = np.nanmedian(image)
1273 image[badPixels] = fallbackValue
1274 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)
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.
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)