Coverage for python / lsst / ip / diffim / makeKernel.py: 20%
147 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-25 08:36 +0000
1# This file is part of ip_diffim.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
22__all__ = ["MakeKernelConfig", "MakeKernelTask"]
24import numpy as np
26import lsst.afw.detection
27import lsst.afw.image
28import lsst.afw.math
29import lsst.afw.table
30import lsst.daf.base
31import lsst.geom
32from lsst.meas.algorithms import SourceDetectionTask, SubtractBackgroundTask
33from lsst.meas.base import SingleFrameMeasurementTask
34from lsst.pex.exceptions import InvalidParameterError, RangeError
35import lsst.pex.config
36import lsst.pipe.base
38from .makeKernelBasisList import makeKernelBasisList
39from .psfMatch import PsfMatchConfig, PsfMatchTask, PsfMatchConfigAL, PsfMatchConfigDF
41from . import diffimLib
42from .utils import evaluateMeanPsfFwhm, getPsfFwhm
45class MakeKernelConfig(PsfMatchConfig):
46 kernel = lsst.pex.config.ConfigChoiceField(
47 doc="kernel type",
48 typemap=dict(
49 AL=PsfMatchConfigAL,
50 DF=PsfMatchConfigDF
51 ),
52 default="AL",
53 )
54 selectDetection = lsst.pex.config.ConfigurableField(
55 target=SourceDetectionTask,
56 doc="Initial detections used to feed stars to kernel fitting",
57 )
58 selectMeasurement = lsst.pex.config.ConfigurableField(
59 target=SingleFrameMeasurementTask,
60 doc="Initial measurements used to feed stars to kernel fitting",
61 )
62 fwhmExposureGrid = lsst.pex.config.Field(
63 doc="Grid size to compute the average PSF FWHM in an exposure",
64 dtype=int,
65 default=10,
66 )
67 fwhmExposureBuffer = lsst.pex.config.Field(
68 doc="Fractional buffer margin to be left out of all sides of the image during construction"
69 "of grid to compute average PSF FWHM in an exposure",
70 dtype=float,
71 default=0.05,
72 )
74 def setDefaults(self):
75 # High sigma detections only
76 self.selectDetection.reEstimateBackground = False
77 self.selectDetection.thresholdValue = 10.0
78 self.selectDetection.excludeMaskPlanes = ["EDGE"]
80 # Minimal set of measurments for star selection
81 self.selectMeasurement.algorithms.names.clear()
82 self.selectMeasurement.algorithms.names = ('base_SdssCentroid', 'base_PsfFlux', 'base_PixelFlags',
83 'base_SdssShape', 'base_GaussianFlux', 'base_SkyCoord',
84 'base_ClassificationSizeExtendedness')
85 self.selectMeasurement.slots.modelFlux = None
86 self.selectMeasurement.slots.apFlux = None
87 self.selectMeasurement.slots.calibFlux = None
90class MakeKernelTask(PsfMatchTask):
91 """Construct a kernel for PSF matching two exposures.
92 """
94 ConfigClass = MakeKernelConfig
95 _DefaultName = "makeALKernel"
97 def __init__(self, *args, **kwargs):
98 super().__init__(*args, **kwargs)
99 # the background subtraction task uses a config from an unusual location,
100 # so cannot easily be constructed with makeSubtask
101 self.background = SubtractBackgroundTask(config=self.kConfig.afwBackgroundConfig, name="background",
102 parentTask=self)
104 self.selectSchema = lsst.afw.table.SourceTable.makeMinimalSchema()
105 self.selectAlgMetadata = lsst.daf.base.PropertyList()
106 self.makeSubtask("selectDetection", schema=self.selectSchema)
107 self.makeSubtask("selectMeasurement", schema=self.selectSchema, algMetadata=self.selectAlgMetadata)
109 def run(self, template, science, kernelSources, preconvolved=False,
110 templateFwhmPix=None, scienceFwhmPix=None):
111 """Solve for the kernel and background model that best match two
112 Exposures evaluated at the given source locations.
114 Parameters
115 ----------
116 template : `lsst.afw.image.Exposure`
117 Exposure that will be convolved.
118 science : `lsst.afw.image.Exposure`
119 The exposure that will be matched.
120 kernelSources : `lsst.afw.table.SourceCatalog`
121 Kernel candidate sources with appropriately sized footprints.
122 Typically the output of `MakeKernelTask.selectKernelSources`.
123 preconvolved : `bool`, optional
124 Was the science image convolved with its own PSF?
126 Returns
127 -------
128 results : `lsst.pipe.base.Struct`
130 ``psfMatchingKernel`` : `lsst.afw.math.LinearCombinationKernel`
131 Spatially varying Psf-matching kernel.
132 ``backgroundModel`` : `lsst.afw.math.Function2D`
133 Spatially varying background-matching function.
134 """
135 kernelCellSet = self._buildCellSet(template.maskedImage, science.maskedImage, kernelSources)
136 if (scienceFwhmPix is None) or (templateFwhmPix is None):
137 # Calling getPsfFwhm on template.psf fails on some rare occasions when
138 # the template has no input exposures at the average position of the
139 # stars. So we try getPsfFwhm first on template, and if that fails we
140 # evaluate the PSF on a grid specified by fwhmExposure* fields.
141 # To keep consistent definitions for PSF size on the template and
142 # science images, we use the same method for both.
143 try:
144 templateFwhmPix = getPsfFwhm(template.psf)
145 scienceFwhmPix = getPsfFwhm(science.psf)
146 except (InvalidParameterError, RangeError):
147 self.log.debug("Unable to evaluate PSF at the average position. "
148 "Evaluting PSF on a grid of points."
149 )
150 templateFwhmPix = evaluateMeanPsfFwhm(template,
151 fwhmExposureBuffer=self.config.fwhmExposureBuffer,
152 fwhmExposureGrid=self.config.fwhmExposureGrid
153 )
154 scienceFwhmPix = evaluateMeanPsfFwhm(science,
155 fwhmExposureBuffer=self.config.fwhmExposureBuffer,
156 fwhmExposureGrid=self.config.fwhmExposureGrid
157 )
159 if preconvolved:
160 scienceFwhmPix *= np.sqrt(2)
161 basisList = self.makeKernelBasisList(templateFwhmPix, scienceFwhmPix,
162 metadata=self.metadata)
163 spatialSolution, psfMatchingKernel, backgroundModel = self._solve(kernelCellSet, basisList)
164 return lsst.pipe.base.Struct(
165 psfMatchingKernel=psfMatchingKernel,
166 backgroundModel=backgroundModel,
167 )
169 def selectKernelSources(self, template, science, candidateList=None, preconvolved=False,
170 templateFwhmPix=None, scienceFwhmPix=None):
171 """Select sources from a list of candidates, and extract footprints.
173 Parameters
174 ----------
175 template : `lsst.afw.image.Exposure`
176 Exposure that will be convolved.
177 science : `lsst.afw.image.Exposure`
178 The exposure that will be matched.
179 candidateList : `lsst.afw.table.SourceCatalog`
180 Sources to check as possible kernel candidates.
181 preconvolved : `bool`, optional
182 Was the science image convolved with its own PSF?
184 Returns
185 -------
186 kernelSources : `lsst.afw.table.SourceCatalog`
187 Kernel candidates with appropriate sized footprints.
188 """
189 if (scienceFwhmPix is None) or (templateFwhmPix is None):
190 # Calling getPsfFwhm on template.psf fails on some rare occasions when
191 # the template has no input exposures at the average position of the
192 # stars. So we try getPsfFwhm first on template, and if that fails we
193 # evaluate the PSF on a grid specified by fwhmExposure* fields.
194 # To keep consistent definitions for PSF size on the template and
195 # science images, we use the same method for both.
196 try:
197 templateFwhmPix = getPsfFwhm(template.psf)
198 scienceFwhmPix = getPsfFwhm(science.psf)
199 except (InvalidParameterError, RangeError):
200 self.log.debug("Unable to evaluate PSF at the average position. "
201 "Evaluting PSF on a grid of points."
202 )
203 templateFwhmPix = evaluateMeanPsfFwhm(template,
204 fwhmExposureBuffer=self.config.fwhmExposureBuffer,
205 fwhmExposureGrid=self.config.fwhmExposureGrid
206 )
207 scienceFwhmPix = evaluateMeanPsfFwhm(science,
208 fwhmExposureBuffer=self.config.fwhmExposureBuffer,
209 fwhmExposureGrid=self.config.fwhmExposureGrid
210 )
211 if preconvolved:
212 scienceFwhmPix *= np.sqrt(2)
213 kernelSize = self.makeKernelBasisList(templateFwhmPix, scienceFwhmPix)[0].getWidth()
214 kernelSources = self.makeCandidateList(template, science, kernelSize,
215 candidateList=candidateList,
216 preconvolved=preconvolved)
217 return kernelSources
219 def getSelectSources(self, exposure, sigma=None, doSmooth=True, idFactory=None):
220 """Get sources to use for Psf-matching.
222 This method runs detection and measurement on an exposure.
223 The returned set of sources will be used as candidates for
224 Psf-matching.
226 Parameters
227 ----------
228 exposure : `lsst.afw.image.Exposure`
229 Exposure on which to run detection/measurement
230 sigma : `float`, optional
231 PSF sigma, in pixels, used for smoothing the image for detection.
232 If `None`, the PSF width will be used.
233 doSmooth : `bool`
234 Whether or not to smooth the Exposure with Psf before detection
235 idFactory : `lsst.afw.table.IdFactory`
236 Factory for the generation of Source ids
238 Returns
239 -------
240 selectSources :
241 source catalog containing candidates for the Psf-matching
242 """
243 if idFactory:
244 table = lsst.afw.table.SourceTable.make(self.selectSchema, idFactory)
245 else:
246 table = lsst.afw.table.SourceTable.make(self.selectSchema)
247 mi = exposure.getMaskedImage()
249 imArr = mi.image.array
250 maskArr = mi.mask.array
251 miArr = np.ma.masked_array(imArr, mask=maskArr)
252 try:
253 fitBg = self.background.fitBackground(mi)
254 bkgd = fitBg.getImageF(self.background.config.algorithm,
255 self.background.config.undersampleStyle)
256 except Exception:
257 self.log.warning("Failed to get background model. Falling back to median background estimation")
258 bkgd = np.ma.median(miArr)
260 # Take off background for detection
261 mi -= bkgd
262 try:
263 table.setMetadata(self.selectAlgMetadata)
264 detRet = self.selectDetection.run(
265 table=table,
266 exposure=exposure,
267 sigma=sigma,
268 doSmooth=doSmooth
269 )
270 selectSources = detRet.sources
271 self.selectMeasurement.run(measCat=selectSources, exposure=exposure)
272 finally:
273 # Put back on the background in case it is needed down stream
274 mi += bkgd
275 del bkgd
277 self.log.info("Selected %d sources via detection measurement.", len(selectSources))
278 return selectSources
280 def makeCandidateList(self, convolved, reference, kernelSize,
281 candidateList, preconvolved=False, sigma=None):
282 """Make a list of acceptable KernelCandidates.
284 Generate a list of candidate sources for Psf-matching, remove sources
285 with bad pixel masks set or that extend off the image.
287 Parameters
288 ----------
289 convolved : `lsst.afw.image.Exposure`
290 Exposure that will be convolved. This is typically the template
291 image, and may have a large bbox than the reference exposure.
292 reference : `lsst.afw.image.Exposure`
293 Exposure that will be matched-to. This is typically the science
294 image.
295 kernelSize : `float`
296 Dimensions of the Psf-matching Kernel, used to set detection
297 footprints.
298 candidateList : `lsst.afw.table.SourceCatalog`
299 List of Sources to examine for kernel candidacy.
300 preconvolved : `bool`, optional
301 Was the science exposure already convolved with its PSF?
303 Returns
304 -------
305 candidates : `lsst.afw.table.SourceCatalog`
306 Candidates with footprints extended to a ``kernelSize`` box.
308 Raises
309 ------
310 RuntimeError
311 If ``candidateList`` is empty after sub-selection.
312 """
313 if candidateList is None:
314 candidateList = self.getSelectSources(reference, doSmooth=not preconvolved, sigma=sigma)
315 if len(candidateList) < 1:
316 raise RuntimeError("No kernel candidates after detection and measurement.")
318 bitmask = reference.mask.getPlaneBitMask(self.config.badMaskPlanes)
319 good = np.ones(len(candidateList), dtype=bool)
320 # Make all candidates have the same size footprint, based on kernelSize.
321 for i, candidate in enumerate(candidateList):
322 # Only use the brightest peak; the input should be pre-deblended!
323 peak = candidate.getFootprint().getPeaks()[0]
324 size = 2*kernelSize + 1 # ensure the resulting box is odd
325 bbox = lsst.geom.Box2I.makeCenteredBox(candidate.getCentroid(),
326 lsst.geom.Extent2I(size, size))
327 boxFootprint = lsst.afw.detection.Footprint(lsst.afw.geom.SpanSet(bbox))
328 boxFootprint.addPeak(peak.getFx(), peak.getFy(), peak.getPeakValue())
329 candidate.setFootprint(boxFootprint)
331 # Reject footprints not contained in either image.
332 if not reference.getBBox().contains(bbox) or not convolved.getBBox().contains(bbox):
333 good[i] = False
334 continue
335 # Reject footprints with any bad mask bits set.
336 if (reference.subset(bbox).mask.array & bitmask).any():
337 good[i] = False
338 continue
340 # Reject footprints with any bad mask bits set.
341 if (convolved.subset(bbox).mask.array & bitmask).any():
342 good[i] = False
343 continue
344 candidates = candidateList[good].copy(deep=True)
346 self.log.info("Selected %d / %d sources as kernel candidates.", good.sum(), len(candidateList))
348 if len(candidates) < 1:
349 raise RuntimeError("No good kernel candidates available.")
351 return candidates
353 def makeKernelBasisList(self, targetFwhmPix=None, referenceFwhmPix=None,
354 basisDegGauss=None, basisSigmaGauss=None, metadata=None):
355 """Wrapper to set log messages for
356 `lsst.ip.diffim.makeKernelBasisList`.
358 Parameters
359 ----------
360 targetFwhmPix : `float`, optional
361 Passed on to `lsst.ip.diffim.generateAlardLuptonBasisList`.
362 Not used for delta function basis sets.
363 referenceFwhmPix : `float`, optional
364 Passed on to `lsst.ip.diffim.generateAlardLuptonBasisList`.
365 Not used for delta function basis sets.
366 basisDegGauss : `list` of `int`, optional
367 Passed on to `lsst.ip.diffim.generateAlardLuptonBasisList`.
368 Not used for delta function basis sets.
369 basisSigmaGauss : `list` of `int`, optional
370 Passed on to `lsst.ip.diffim.generateAlardLuptonBasisList`.
371 Not used for delta function basis sets.
372 metadata : `lsst.daf.base.PropertySet`, optional
373 Passed on to `lsst.ip.diffim.generateAlardLuptonBasisList`.
374 Not used for delta function basis sets.
376 Returns
377 -------
378 basisList: `list` of `lsst.afw.math.kernel.FixedKernel`
379 List of basis kernels.
380 """
381 basisList = makeKernelBasisList(self.kConfig,
382 targetFwhmPix=targetFwhmPix,
383 referenceFwhmPix=referenceFwhmPix,
384 basisDegGauss=basisDegGauss,
385 basisSigmaGauss=basisSigmaGauss,
386 metadata=metadata)
387 if targetFwhmPix == referenceFwhmPix:
388 self.log.info("Target and reference psf fwhms are equal, falling back to config values")
389 elif referenceFwhmPix > targetFwhmPix:
390 self.log.info("Reference psf fwhm is the greater, normal convolution mode")
391 else:
392 self.log.info("Target psf fwhm is the greater, deconvolution mode")
394 return basisList
396 def _buildCellSet(self, convolved, reference, candidateList):
397 """Build a SpatialCellSet for use with the solve method.
399 Parameters
400 ----------
401 convolved : `lsst.afw.image.MaskedImage`
402 MaskedImage to PSF-matched to reference.
403 reference : `lsst.afw.image.MaskedImage`
404 Reference MaskedImage.
405 candidateList : `lsst.afw.table.SourceCatalog`
406 Kernel candidate sources with footprints.
408 Returns
409 -------
410 kernelCellSet : `lsst.afw.math.SpatialCellSet`
411 A SpatialCellSet for use with self._solve.
412 """
413 sizeCellX, sizeCellY = self._adaptCellSize(candidateList)
415 imageBBox = convolved.getBBox()
416 imageBBox.clip(reference.getBBox())
417 # Object to store the KernelCandidates for spatial modeling
418 kernelCellSet = lsst.afw.math.SpatialCellSet(imageBBox, sizeCellX, sizeCellY)
420 candidateConfig = lsst.pex.config.makePropertySet(self.kConfig)
421 # Place candidates within the spatial grid
422 for candidate in candidateList:
423 bbox = candidate.getFootprint().getBBox()
424 templateCutout = lsst.afw.image.MaskedImageF(convolved, bbox)
425 scienceCutout = lsst.afw.image.MaskedImageF(reference, bbox)
427 kernelCandidate = diffimLib.makeKernelCandidate(candidate,
428 templateCutout,
429 scienceCutout,
430 candidateConfig)
432 self.log.debug("Candidate %d at %.2f, %.2f rating=%f",
433 kernelCandidate.getId(),
434 kernelCandidate.getXCenter(),
435 kernelCandidate.getYCenter(),
436 kernelCandidate.getCandidateRating())
437 kernelCellSet.insertCandidate(kernelCandidate)
439 return kernelCellSet
441 def _adaptCellSize(self, candidateList):
442 """NOT IMPLEMENTED YET.
444 Parameters
445 ----------
446 candidateList : `list`
447 A list of footprints/maskedImages for kernel candidates;
449 Returns
450 -------
451 sizeCellX, sizeCellY : `int`
452 New dimensions to use for the kernel.
453 """
454 return self.kConfig.sizeCellX, self.kConfig.sizeCellY