Coverage for python / lsst / ip / diffim / modelPsfMatch.py: 17%
187 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-24 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/>.
22import numpy as np
24from . import diffimLib
25import lsst.afw.display as afwDisplay
26import lsst.afw.image as afwImage
27import lsst.afw.math as afwMath
28import lsst.geom as geom
29import lsst.pex.config as pexConfig
30import lsst.pex.exceptions as pexExceptions
31import lsst.pipe.base as pipeBase
32from lsst.utils.logging import getTraceLogger
33from lsst.utils.timer import timeMethod
34from .makeKernelBasisList import makeKernelBasisList
35from .psfMatch import PsfMatchTask, PsfMatchConfigAL
36from . import utils as dituils
38__all__ = (
39 "ModelPsfMatchTask",
40 "ModelPsfMatchConfig",
41 "WarpedPsfTransformTooBigError",
42 "PsfComputeShapeError",
43)
45sigma2fwhm = 2.*np.sqrt(2.*np.log(2.))
48def _safeEstimateRoughShape(psf, position, jitter=2.0, maxTries=5, log=None):
50 """
51 Try to compute the shape of `psf` at `position`.
53 If it fails, perturb the position until success or maxTries.
55 Parameters
56 ----------
57 psf : `lsst.afw.detection.Psf`
58 The PSF object.
59 position : `lsst.geom.Point2D`
60 The nominal position at which to evaluate the PSF shape.
61 jitter : `float`
62 Size of random offset in pixels to try when retrying.
63 maxTries : `int`
64 Maximum number of attempts (including the original).
65 log : `lsst.utils.logging.LsstLogAdapter` optional
66 logger to use for logging retries.
68 Returns
69 -------
70 shape : lsst.afw.geom.ellipses.Quadrupole
71 The shape of the PSF at the (possibly perturbed) position.
72 """
73 try:
74 return psf.computeShape(position)
75 except pexExceptions.Exception as err:
76 if log:
77 log.info(f"safeEstimateRoughShape: initial evaluation of PSF failed with message: {str(err)}."
78 " Retrying nearby")
79 firstErr = err
81 # random offsets must be deterministic
82 offsets = np.random.default_rng(seed=40).uniform(-jitter, jitter, size=(maxTries - 1, 2))
83 candidates = [position + geom.Extent2D(dx, dy) for dx, dy in offsets]
84 for i, candidate in enumerate(candidates):
85 try:
86 shape = psf.computeShape(candidate)
87 if log:
88 log.info(f"safeEstimateRoughShape succeeded on try {i + 2} at position {candidate}")
89 return shape
90 except pexExceptions.Exception:
91 continue
93 # If nothing worked, raise AlgorithmError from original position
94 raise PsfComputeShapeError(position) from firstErr
97def nextOddInteger(x):
98 nextInt = int(np.ceil(x))
99 return nextInt + 1 if nextInt%2 == 0 else nextInt
102class WarpedPsfTransformTooBigError(pipeBase.AlgorithmError):
103 """Raised when the transform of a WarpedPsf is too large to compute
104 the FWHM of the PSF at a given position.
105 """
106 @property
107 def metadata(self) -> dict:
108 return {}
111class PsfComputeShapeError(pipeBase.AlgorithmError):
112 def __init__(self, position):
113 message = f"Unable to evaluate Psf at ({float(position.getX()):.4f}, {float(position.getY()):.4f})"
114 super().__init__(message)
115 self.position = position
117 @property
118 def metadata(self) -> dict:
119 return {
120 # must have values `str`, `int`, `float`, `bool`or nested-dict
121 "position": {"x": float(self.position.getX()),
122 "y": float(self.position.getY())}
123 }
126class ModelPsfMatchConfig(pexConfig.Config):
127 """Configuration for model-to-model Psf matching"""
129 kernel = pexConfig.ConfigChoiceField(
130 doc="kernel type",
131 typemap=dict(
132 AL=PsfMatchConfigAL,
133 ),
134 default="AL",
135 )
136 doAutoPadPsf = pexConfig.Field(
137 dtype=bool,
138 doc=("If too small, automatically pad the science Psf? "
139 "Pad to smallest dimensions appropriate for the matching kernel dimensions, "
140 "as specified by autoPadPsfTo. If false, pad by the padPsfBy config."),
141 default=True,
142 )
143 autoPadPsfTo = pexConfig.RangeField(
144 dtype=float,
145 doc=("Minimum Science Psf dimensions as a fraction of matching kernel dimensions. "
146 "If the dimensions of the Psf to be matched are less than the "
147 "matching kernel dimensions * autoPadPsfTo, pad Science Psf to this size. "
148 "Ignored if doAutoPadPsf=False."),
149 default=1.4,
150 min=1.0,
151 max=2.0
152 )
153 padPsfBy = pexConfig.Field(
154 dtype=int,
155 doc="Pixels (even) to pad Science Psf by before matching. Ignored if doAutoPadPsf=True",
156 default=0,
157 )
159 def setDefaults(self):
160 # No sigma clipping
161 self.kernel.active.singleKernelClipping = False
162 self.kernel.active.kernelSumClipping = False
163 self.kernel.active.spatialKernelClipping = False
164 self.kernel.active.checkConditionNumber = False
166 # Variance is ill defined
167 self.kernel.active.constantVarianceWeighting = True
169 # Do not change specified kernel size
170 self.kernel.active.scaleByFwhm = False
173class ModelPsfMatchTask(PsfMatchTask):
174 """Match two model Psfs, and application of the Psf-matching kernel
175 to an input Exposure.
176 """
177 ConfigClass = ModelPsfMatchConfig
179 def __init__(self, *args, **kwargs):
180 PsfMatchTask.__init__(self, *args, **kwargs)
181 self.kConfig = self.config.kernel.active
183 @timeMethod
184 def run(self, exposure, referencePsfModel, kernelSum=1.0):
185 """Psf-match an exposure to a model Psf.
187 Parameters
188 ----------
189 exposure : `lsst.afw.image.Exposure`
190 Exposure to Psf-match to the reference Psf model;
191 it must return a valid PSF model via exposure.getPsf()
192 referencePsfModel : `lsst.afw.detection.Psf`
193 The Psf model to match to
194 kernelSum : `float`, optional
195 A multipicative factor to apply to the kernel sum (default=1.0)
197 Returns
198 -------
199 result : `struct`
200 - ``psfMatchedExposure`` : the Psf-matched Exposure.
201 This has the same parent bbox, Wcs, PhotoCalib and
202 Filter as the input Exposure but no Psf.
203 In theory the Psf should equal referencePsfModel but
204 the match is likely not exact.
205 - ``psfMatchingKernel`` : the spatially varying Psf-matching kernel
206 - ``kernelCellSet`` : SpatialCellSet used to solve for the Psf-matching kernel
207 - ``referencePsfModel`` : Validated and/or modified reference model used
209 Raises
210 ------
211 RuntimeError
212 if the Exposure does not contain a Psf model
213 """
214 if not exposure.hasPsf():
215 raise RuntimeError("exposure does not contain a Psf model")
217 maskedImage = exposure.getMaskedImage()
219 self.log.info("compute Psf-matching kernel")
220 result = self._buildCellSet(exposure, referencePsfModel)
221 kernelCellSet = result.kernelCellSet
222 referencePsfModel = result.referencePsfModel
223 # TODO: This should be evaluated at (or close to) the center of the
224 # exposure's bounding box in DM-32756.
225 sciAvgPos = exposure.getPsf().getAveragePosition()
226 modelAvgPos = referencePsfModel.getAveragePosition()
228 # TODO If DM-52537 eliminates the need for _safeEstimateRoughShape,
229 # the remove it, and replace next line with a direct call to
230 # computeShape in a try/except block.
231 fwhmScience = _safeEstimateRoughShape(exposure.getPsf(), sciAvgPos,
232 log=self.log).getDeterminantRadius()*sigma2fwhm
233 fwhmModel = referencePsfModel.computeShape(modelAvgPos).getDeterminantRadius()*sigma2fwhm
235 basisList = makeKernelBasisList(self.kConfig, fwhmScience, fwhmModel, metadata=self.metadata)
236 spatialSolution, psfMatchingKernel, backgroundModel = self._solve(kernelCellSet, basisList)
238 if psfMatchingKernel.isSpatiallyVarying():
239 sParameters = np.array(psfMatchingKernel.getSpatialParameters())
240 sParameters[0][0] = kernelSum
241 psfMatchingKernel.setSpatialParameters(sParameters)
242 else:
243 kParameters = np.array(psfMatchingKernel.getKernelParameters())
244 kParameters[0] = kernelSum
245 psfMatchingKernel.setKernelParameters(kParameters)
247 self.log.info("Psf-match science exposure to reference")
248 psfMatchedExposure = afwImage.ExposureF(exposure.getBBox(), exposure.getWcs())
249 psfMatchedExposure.info.id = exposure.info.id
250 psfMatchedExposure.setFilter(exposure.getFilter())
251 psfMatchedExposure.setPhotoCalib(exposure.getPhotoCalib())
252 psfMatchedExposure.getInfo().setVisitInfo(exposure.getInfo().getVisitInfo())
253 psfMatchedExposure.setPsf(referencePsfModel)
254 psfMatchedMaskedImage = psfMatchedExposure.getMaskedImage()
256 # Normalize the psf-matching kernel while convolving since its magnitude is meaningless
257 # when PSF-matching one model to another.
258 convolutionControl = afwMath.ConvolutionControl()
259 convolutionControl.setDoNormalize(True)
260 afwMath.convolve(psfMatchedMaskedImage, maskedImage, psfMatchingKernel, convolutionControl)
262 self.log.info("done")
263 return pipeBase.Struct(psfMatchedExposure=psfMatchedExposure,
264 psfMatchingKernel=psfMatchingKernel,
265 kernelCellSet=kernelCellSet,
266 metadata=self.metadata,
267 )
269 def _diagnostic(self, kernelCellSet, spatialSolution, spatialKernel, spatialBg):
270 """Print diagnostic information on spatial kernel and background fit
272 The debugging diagnostics are not really useful here, since the images we are matching have
273 no variance. Thus override the _diagnostic method to generate no logging information"""
274 return
276 def _buildCellSet(self, exposure, referencePsfModel):
277 """Build a SpatialCellSet for use with the solve method
279 Parameters
280 ----------
281 exposure : `lsst.afw.image.Exposure`
282 The science exposure that will be convolved; must contain a Psf
283 referencePsfModel : `lsst.afw.detection.Psf`
284 Psf model to match to
286 Returns
287 -------
288 result : `struct`
289 - ``kernelCellSet`` : a SpatialCellSet to be used by self._solve
290 - ``referencePsfModel`` : Validated and/or modified
291 reference model used to populate the SpatialCellSet
293 Notes
294 -----
295 If the reference Psf model and science Psf model have different dimensions,
296 adjust the referencePsfModel (the model to which the exposure PSF will be matched)
297 to match that of the science Psf. If the science Psf dimensions vary across the image,
298 as is common with a WarpedPsf, either pad or clip (depending on config.padPsf)
299 the dimensions to be constant.
300 """
301 sizeCellX = self.kConfig.sizeCellX
302 sizeCellY = self.kConfig.sizeCellY
304 scienceBBox = exposure.getBBox()
305 # Extend for proper spatial matching kernel all the way to edge, especially for narrow strips
306 scienceBBox.grow(geom.Extent2I(sizeCellX, sizeCellY))
308 sciencePsfModel = exposure.getPsf()
310 dimenR = referencePsfModel.getLocalKernel(scienceBBox.getCenter()).getDimensions()
312 regionSizeX, regionSizeY = scienceBBox.getDimensions()
313 scienceX0, scienceY0 = scienceBBox.getMin()
315 kernelCellSet = afwMath.SpatialCellSet(geom.Box2I(scienceBBox), sizeCellX, sizeCellY)
317 nCellX = regionSizeX//sizeCellX
318 nCellY = regionSizeY//sizeCellY
320 if nCellX == 0 or nCellY == 0:
321 raise ValueError("Exposure dimensions=%s and sizeCell=(%s, %s). Insufficient area to match" %
322 (scienceBBox.getDimensions(), sizeCellX, sizeCellY))
324 # Survey the PSF dimensions of the Spatial Cell Set
325 # to identify the minimum enclosed or maximum bounding square BBox.
326 widthList = []
327 heightList = []
328 for row in range(nCellY):
329 posY = sizeCellY*row + sizeCellY//2 + scienceY0
330 for col in range(nCellX):
331 posX = sizeCellX*col + sizeCellX//2 + scienceX0
332 try:
333 widthS, heightS = sciencePsfModel.computeBBox(geom.Point2D(posX, posY)).getDimensions()
334 except pexExceptions.InvalidParameterError as err:
335 raise PsfComputeShapeError(geom.Point2D(posX, posY)) from err
336 widthList.append(widthS)
337 heightList.append(heightS)
339 psfSize = max(max(heightList), max(widthList))
341 if self.config.doAutoPadPsf:
342 minPsfSize = nextOddInteger(self.kConfig.kernelSize*self.config.autoPadPsfTo)
343 paddingPix = max(0, minPsfSize - psfSize)
344 else:
345 if self.config.padPsfBy % 2 != 0:
346 raise ValueError("Config padPsfBy (%i pixels) must be even number." %
347 self.config.padPsfBy)
348 paddingPix = self.config.padPsfBy
350 if paddingPix > 0:
351 self.log.debug("Padding Science PSF from (%d, %d) to (%d, %d) pixels",
352 psfSize, psfSize, paddingPix + psfSize, paddingPix + psfSize)
353 psfSize += paddingPix
355 # Check that PSF is larger than the matching kernel
356 maxKernelSize = psfSize - 1
357 if maxKernelSize % 2 == 0:
358 maxKernelSize -= 1
359 if self.kConfig.kernelSize > maxKernelSize:
360 message = """
361 Kernel size (%d) too big to match Psfs of size %d.
362 Please reconfigure by setting one of the following:
363 1) kernel size to <= %d
364 2) doAutoPadPsf=True
365 3) padPsfBy to >= %s
366 """ % (self.kConfig.kernelSize, psfSize,
367 maxKernelSize, self.kConfig.kernelSize - maxKernelSize)
368 raise ValueError(message)
370 dimenS = geom.Extent2I(psfSize, psfSize)
372 if (dimenR != dimenS):
373 try:
374 referencePsfModel = referencePsfModel.resized(psfSize, psfSize)
375 self.log.info("Adjusted dimensions of reference PSF model from %s to %s", dimenR, dimenS)
376 except Exception as e:
377 self.log.warning("Zero padding or clipping the reference PSF model of type %s and dimensions"
378 " %s to the science Psf dimensions %s because: %s",
379 referencePsfModel.__class__.__name__, dimenR, dimenS, e)
380 dimenR = dimenS
382 ps = pexConfig.makePropertySet(self.kConfig)
383 for row in range(nCellY):
384 # place at center of cell
385 posY = sizeCellY*row + sizeCellY//2 + scienceY0
387 for col in range(nCellX):
388 # place at center of cell
389 posX = sizeCellX*col + sizeCellX//2 + scienceX0
391 getTraceLogger(self.log, 4).debug("Creating Psf candidate at %.1f %.1f", posX, posY)
393 # reference kernel image, at location of science subimage
394 referenceMI = self._makePsfMaskedImage(referencePsfModel, posX, posY, dimensions=dimenR)
396 # kernel image we are going to convolve
397 scienceMI = self._makePsfMaskedImage(sciencePsfModel, posX, posY, dimensions=dimenR)
399 # The image to convolve is the science image, to the reference Psf.
400 kc = diffimLib.makeKernelCandidate(posX, posY, scienceMI, referenceMI, ps)
401 kernelCellSet.insertCandidate(kc)
403 import lsstDebug
404 display = lsstDebug.Info(__name__).display
405 displaySpatialCells = lsstDebug.Info(__name__).displaySpatialCells
406 maskTransparency = lsstDebug.Info(__name__).maskTransparency
407 if not maskTransparency:
408 maskTransparency = 0
409 if display:
410 afwDisplay.setDefaultMaskTransparency(maskTransparency)
411 if display and displaySpatialCells:
412 dituils.showKernelSpatialCells(exposure.getMaskedImage(), kernelCellSet,
413 symb="o", ctype=afwDisplay.CYAN, ctypeUnused=afwDisplay.YELLOW,
414 ctypeBad=afwDisplay.RED, size=4, frame=lsstDebug.frame,
415 title="Image to be convolved")
416 lsstDebug.frame += 1
417 return pipeBase.Struct(kernelCellSet=kernelCellSet,
418 referencePsfModel=referencePsfModel,
419 )
421 def _makePsfMaskedImage(self, psfModel, posX, posY, dimensions=None):
422 """Return a MaskedImage of the a PSF Model of specified dimensions
423 """
424 rawKernel = psfModel.computeKernelImage(geom.Point2D(posX, posY)).convertF()
425 if dimensions is None:
426 dimensions = rawKernel.getDimensions()
427 if rawKernel.getDimensions() == dimensions:
428 kernelIm = rawKernel
429 else:
430 # make image of proper size
431 kernelIm = afwImage.ImageF(dimensions)
432 bboxToPlace = geom.Box2I(geom.Point2I((dimensions.getX() - rawKernel.getWidth())//2,
433 (dimensions.getY() - rawKernel.getHeight())//2),
434 rawKernel.getDimensions())
435 kernelIm.assign(rawKernel, bboxToPlace)
437 kernelMask = afwImage.Mask(dimensions, 0x0)
438 kernelVar = afwImage.ImageF(dimensions, 1.0)
439 return afwImage.MaskedImageF(kernelIm, kernelMask, kernelVar)