Coverage for python / lsst / ip / diffim / makeKernel.py: 22%
141 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-05-07 08:35 +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?
125 templateFwhmPix, scienceFwhmPix : `float` or `None`, optional
126 FWHM of the template or science PSF, in pixels.
127 Will be recalculated if not specified.
129 Returns
130 -------
131 results : `lsst.pipe.base.Struct`
133 ``psfMatchingKernel`` : `lsst.afw.math.LinearCombinationKernel`
134 Spatially varying Psf-matching kernel.
135 ``backgroundModel`` : `lsst.afw.math.Function2D`
136 Spatially varying background-matching function.
137 """
138 kernelCellSet = self._buildCellSet(template.maskedImage, science.maskedImage, kernelSources)
139 templateFwhmPix, scienceFwhmPix = self._checkPsfWidths(template, science,
140 templateFwhmPix, scienceFwhmPix,
141 preconvolved=preconvolved)
142 basisList = self.makeKernelBasisList(templateFwhmPix, scienceFwhmPix,
143 metadata=self.metadata)
144 spatialSolution, psfMatchingKernel, backgroundModel = self._solve(kernelCellSet, basisList)
145 return lsst.pipe.base.Struct(
146 psfMatchingKernel=psfMatchingKernel,
147 backgroundModel=backgroundModel,
148 )
150 def selectKernelSources(self, template, science, candidateList=None, preconvolved=False,
151 templateFwhmPix=None, scienceFwhmPix=None):
152 """Select sources from a list of candidates, and extract footprints.
154 Parameters
155 ----------
156 template : `lsst.afw.image.Exposure`
157 Exposure that will be convolved.
158 science : `lsst.afw.image.Exposure`
159 The exposure that will be matched.
160 candidateList : `lsst.afw.table.SourceCatalog`
161 Sources to check as possible kernel candidates.
162 preconvolved : `bool`, optional
163 Was the science image convolved with its own PSF?
164 templateFwhmPix, scienceFwhmPix : `float` or `None`, optional
165 FWHM of the template or science PSF, in pixels.
166 Will be recalculated if not specified.
168 Returns
169 -------
170 kernelSources : `lsst.afw.table.SourceCatalog`
171 Kernel candidates with appropriate sized footprints.
172 """
173 templateFwhmPix, scienceFwhmPix = self._checkPsfWidths(template, science,
174 templateFwhmPix, scienceFwhmPix,
175 preconvolved=preconvolved)
176 kernelSize = self.makeKernelBasisList(templateFwhmPix, scienceFwhmPix)[0].getWidth()
177 kernelSources = self.makeCandidateList(template, science, kernelSize,
178 candidateList=candidateList,
179 preconvolved=preconvolved)
180 return kernelSources
182 def getSelectSources(self, exposure, sigma=None, doSmooth=True, idFactory=None):
183 """Get sources to use for Psf-matching.
185 This method runs detection and measurement on an exposure.
186 The returned set of sources will be used as candidates for
187 Psf-matching.
189 Parameters
190 ----------
191 exposure : `lsst.afw.image.Exposure`
192 Exposure on which to run detection/measurement
193 sigma : `float`, optional
194 PSF sigma, in pixels, used for smoothing the image for detection.
195 If `None`, the PSF width will be used.
196 doSmooth : `bool`
197 Whether or not to smooth the Exposure with Psf before detection
198 idFactory : `lsst.afw.table.IdFactory`
199 Factory for the generation of Source ids
201 Returns
202 -------
203 selectSources :
204 source catalog containing candidates for the Psf-matching
205 """
206 if idFactory:
207 table = lsst.afw.table.SourceTable.make(self.selectSchema, idFactory)
208 else:
209 table = lsst.afw.table.SourceTable.make(self.selectSchema)
210 mi = exposure.getMaskedImage()
212 imArr = mi.image.array
213 maskArr = mi.mask.array
214 miArr = np.ma.masked_array(imArr, mask=maskArr)
215 try:
216 fitBg = self.background.fitBackground(mi)
217 bkgd = fitBg.getImageF(self.background.config.algorithm,
218 self.background.config.undersampleStyle)
219 except Exception:
220 self.log.warning("Failed to get background model. Falling back to median background estimation")
221 bkgd = np.ma.median(miArr)
223 # Take off background for detection
224 mi -= bkgd
225 try:
226 table.setMetadata(self.selectAlgMetadata)
227 detRet = self.selectDetection.run(
228 table=table,
229 exposure=exposure,
230 sigma=sigma,
231 doSmooth=doSmooth
232 )
233 selectSources = detRet.sources
234 self.selectMeasurement.run(measCat=selectSources, exposure=exposure)
235 finally:
236 # Put back on the background in case it is needed down stream
237 mi += bkgd
238 del bkgd
240 self.log.info("Selected %d sources via detection measurement.", len(selectSources))
241 return selectSources
243 def makeCandidateList(self, convolved, reference, kernelSize,
244 candidateList, preconvolved=False, sigma=None):
245 """Make a list of acceptable KernelCandidates.
247 Generate a list of candidate sources for Psf-matching, remove sources
248 with bad pixel masks set or that extend off the image.
250 Parameters
251 ----------
252 convolved : `lsst.afw.image.Exposure`
253 Exposure that will be convolved. This is typically the template
254 image, and may have a large bbox than the reference exposure.
255 reference : `lsst.afw.image.Exposure`
256 Exposure that will be matched-to. This is typically the science
257 image.
258 kernelSize : `float`
259 Dimensions of the Psf-matching Kernel, used to set detection
260 footprints.
261 candidateList : `lsst.afw.table.SourceCatalog`
262 List of Sources to examine for kernel candidacy.
263 preconvolved : `bool`, optional
264 Was the science exposure already convolved with its PSF?
266 Returns
267 -------
268 candidates : `lsst.afw.table.SourceCatalog`
269 Candidates with footprints extended to a ``kernelSize`` box.
271 Raises
272 ------
273 RuntimeError
274 If ``candidateList`` is empty after sub-selection.
275 """
276 if candidateList is None:
277 candidateList = self.getSelectSources(reference, doSmooth=not preconvolved, sigma=sigma)
278 if len(candidateList) < 1:
279 raise RuntimeError("No kernel candidates after detection and measurement.")
281 bitmask = reference.mask.getPlaneBitMask(self.config.badMaskPlanes)
282 good = np.ones(len(candidateList), dtype=bool)
283 # Make all candidates have the same size footprint, based on kernelSize.
284 for i, candidate in enumerate(candidateList):
285 # Only use the brightest peak; the input should be pre-deblended!
286 peak = candidate.getFootprint().getPeaks()[0]
287 size = 2*kernelSize + 1 # ensure the resulting box is odd
288 bbox = lsst.geom.Box2I.makeCenteredBox(candidate.getCentroid(),
289 lsst.geom.Extent2I(size, size))
290 boxFootprint = lsst.afw.detection.Footprint(lsst.afw.geom.SpanSet(bbox))
291 boxFootprint.addPeak(peak.getFx(), peak.getFy(), peak.getPeakValue())
292 candidate.setFootprint(boxFootprint)
294 # Reject footprints not contained in either image.
295 if not reference.getBBox().contains(bbox) or not convolved.getBBox().contains(bbox):
296 good[i] = False
297 continue
298 # Reject footprints with any bad mask bits set.
299 if (reference.subset(bbox).mask.array & bitmask).any():
300 good[i] = False
301 continue
303 # Reject footprints with any bad mask bits set.
304 if (convolved.subset(bbox).mask.array & bitmask).any():
305 good[i] = False
306 continue
307 candidates = candidateList[good].copy(deep=True)
309 self.log.info("Selected %d / %d sources as kernel candidates.", good.sum(), len(candidateList))
311 if len(candidates) < 1:
312 raise RuntimeError("No good kernel candidates available.")
314 return candidates
316 def makeKernelBasisList(self, targetFwhmPix=None, referenceFwhmPix=None,
317 basisDegGauss=None, basisSigmaGauss=None, metadata=None):
318 """Wrapper to set log messages for
319 `lsst.ip.diffim.makeKernelBasisList`.
321 Parameters
322 ----------
323 targetFwhmPix : `float`, optional
324 Passed on to `lsst.ip.diffim.generateAlardLuptonBasisList`.
325 Not used for delta function basis sets.
326 referenceFwhmPix : `float`, optional
327 Passed on to `lsst.ip.diffim.generateAlardLuptonBasisList`.
328 Not used for delta function basis sets.
329 basisDegGauss : `list` of `int`, optional
330 Passed on to `lsst.ip.diffim.generateAlardLuptonBasisList`.
331 Not used for delta function basis sets.
332 basisSigmaGauss : `list` of `int`, optional
333 Passed on to `lsst.ip.diffim.generateAlardLuptonBasisList`.
334 Not used for delta function basis sets.
335 metadata : `lsst.daf.base.PropertySet`, optional
336 Passed on to `lsst.ip.diffim.generateAlardLuptonBasisList`.
337 Not used for delta function basis sets.
339 Returns
340 -------
341 basisList: `list` of `lsst.afw.math.kernel.FixedKernel`
342 List of basis kernels.
343 """
344 basisList = makeKernelBasisList(self.kConfig,
345 targetFwhmPix=targetFwhmPix,
346 referenceFwhmPix=referenceFwhmPix,
347 basisDegGauss=basisDegGauss,
348 basisSigmaGauss=basisSigmaGauss,
349 metadata=metadata)
350 if targetFwhmPix == referenceFwhmPix:
351 self.log.info("Target and reference psf fwhms are equal, falling back to config values")
352 elif referenceFwhmPix > targetFwhmPix:
353 self.log.info("Reference psf fwhm is the greater, normal convolution mode")
354 else:
355 self.log.info("Target psf fwhm is the greater, deconvolution mode")
357 return basisList
359 def _buildCellSet(self, convolved, reference, candidateList):
360 """Build a SpatialCellSet for use with the solve method.
362 Parameters
363 ----------
364 convolved : `lsst.afw.image.MaskedImage`
365 MaskedImage to PSF-matched to reference.
366 reference : `lsst.afw.image.MaskedImage`
367 Reference MaskedImage.
368 candidateList : `lsst.afw.table.SourceCatalog`
369 Kernel candidate sources with footprints.
371 Returns
372 -------
373 kernelCellSet : `lsst.afw.math.SpatialCellSet`
374 A SpatialCellSet for use with self._solve.
375 """
376 sizeCellX, sizeCellY = self._adaptCellSize(candidateList)
378 imageBBox = convolved.getBBox()
379 imageBBox.clip(reference.getBBox())
380 # Object to store the KernelCandidates for spatial modeling
381 kernelCellSet = lsst.afw.math.SpatialCellSet(imageBBox, sizeCellX, sizeCellY)
383 candidateConfig = lsst.pex.config.makePropertySet(self.kConfig)
384 # Place candidates within the spatial grid
385 for candidate in candidateList:
386 bbox = candidate.getFootprint().getBBox()
387 templateCutout = lsst.afw.image.MaskedImageF(convolved, bbox)
388 scienceCutout = lsst.afw.image.MaskedImageF(reference, bbox)
390 kernelCandidate = diffimLib.makeKernelCandidate(candidate,
391 templateCutout,
392 scienceCutout,
393 candidateConfig)
395 self.log.debug("Candidate %d at %.2f, %.2f rating=%f",
396 kernelCandidate.getId(),
397 kernelCandidate.getXCenter(),
398 kernelCandidate.getYCenter(),
399 kernelCandidate.getCandidateRating())
400 kernelCellSet.insertCandidate(kernelCandidate)
402 return kernelCellSet
404 def _adaptCellSize(self, candidateList):
405 """NOT IMPLEMENTED YET.
407 Parameters
408 ----------
409 candidateList : `list`
410 A list of footprints/maskedImages for kernel candidates;
412 Returns
413 -------
414 sizeCellX, sizeCellY : `int`
415 New dimensions to use for the kernel.
416 """
417 return self.kConfig.sizeCellX, self.kConfig.sizeCellY
419 def _checkPsfWidths(self, template, science, templateFwhmPix, scienceFwhmPix, preconvolved=False):
420 """Check the science and template FWHM.
422 Parameters
423 ----------
424 template : `lsst.afw.image.Exposure`
425 Exposure that will be convolved.
426 science : `lsst.afw.image.Exposure`
427 The exposure that will be matched.
428 templateFwhmPix, scienceFwhmPix : `float` or `None`
429 FWHM of the template or science PSF, in pixels.
430 Will be recalculated if not specified.
431 preconvolved : `bool`, optional
432 Was the science image convolved with its own PSF?
434 Returns
435 -------
436 templateFwhmPix, scienceFwhmPix : `float`
437 The desired PSF FWHM
438 """
439 if (scienceFwhmPix is None) or (templateFwhmPix is None):
440 # Calling getPsfFwhm on template.psf fails on some rare occasions when
441 # the template has no input exposures at the average position of the
442 # stars. So we try getPsfFwhm first on template, and if that fails we
443 # evaluate the PSF on a grid specified by fwhmExposure* fields.
444 # To keep consistent definitions for PSF size on the template and
445 # science images, we use the same method for both.
446 try:
447 templateFwhmPix = getPsfFwhm(template.psf)
448 scienceFwhmPix = getPsfFwhm(science.psf)
449 except (InvalidParameterError, RangeError):
450 self.log.debug("Unable to evaluate PSF at the average position. "
451 "Evaluting PSF on a grid of points."
452 )
453 templateFwhmPix = evaluateMeanPsfFwhm(template,
454 fwhmExposureBuffer=self.config.fwhmExposureBuffer,
455 fwhmExposureGrid=self.config.fwhmExposureGrid
456 )
457 scienceFwhmPix = evaluateMeanPsfFwhm(science,
458 fwhmExposureBuffer=self.config.fwhmExposureBuffer,
459 fwhmExposureGrid=self.config.fwhmExposureGrid
460 )
461 # Apply the Gaussian approximation FWHM boost only when
462 # the size was auto-computed. An explicit ``scienceFwhmPix`` is
463 # assumed already to reflect the preconvolved effective width.
464 if preconvolved:
465 scienceFwhmPix *= np.sqrt(2)
466 return templateFwhmPix, scienceFwhmPix