LSST Applications g04e9c324dd+8c5ae1fdc5,g134cb467dc+1b3060144d,g18429d2f64+f642bf4753,g199a45376c+0ba108daf9,g1fd858c14a+2dcf163641,g262e1987ae+7b8c96d2ca,g29ae962dfc+3bd6ecb08a,g2cef7863aa+aef1011c0b,g35bb328faa+8c5ae1fdc5,g3fd5ace14f+53e1a9e7c5,g4595892280+fef73a337f,g47891489e3+2efcf17695,g4d44eb3520+642b70b07e,g53246c7159+8c5ae1fdc5,g67b6fd64d1+2efcf17695,g67fd3c3899+b70e05ef52,g74acd417e5+317eb4c7d4,g786e29fd12+668abc6043,g87389fa792+8856018cbb,g89139ef638+2efcf17695,g8d7436a09f+3be3c13596,g8ea07a8fe4+9f5ccc88ac,g90f42f885a+a4e7b16d9b,g97be763408+ad77d7208f,g9dd6db0277+b70e05ef52,ga681d05dcb+a3f46e7fff,gabf8522325+735880ea63,gac2eed3f23+2efcf17695,gb89ab40317+2efcf17695,gbf99507273+8c5ae1fdc5,gd8ff7fe66e+b70e05ef52,gdab6d2f7ff+317eb4c7d4,gdc713202bf+b70e05ef52,gdfd2d52018+b10e285e0f,ge365c994fd+310e8507c4,ge410e46f29+2efcf17695,geaed405ab2+562b3308c0,gffca2db377+8c5ae1fdc5,w.2025.35
LSST Data Management Base Package
Loading...
Searching...
No Matches
modelPsfMatch.py
Go to the documentation of this file.
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/>.
21
22import numpy as np
23
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
37
38__all__ = (
39 "ModelPsfMatchTask",
40 "ModelPsfMatchConfig",
41 "WarpedPsfTransformTooBigError",
42 "PsfComputeShapeError",
43)
44
45sigma2fwhm = 2.*np.sqrt(2.*np.log(2.))
46
47
49 nextInt = int(np.ceil(x))
50 return nextInt + 1 if nextInt%2 == 0 else nextInt
51
52
53class WarpedPsfTransformTooBigError(pipeBase.AlgorithmError):
54 """Raised when the transform of a WarpedPsf is too large to compute
55 the FWHM of the PSF at a given position.
56 """
57 def metadata(self) -> dict:
58 return {}
59
60
61class PsfComputeShapeError(pipeBase.AlgorithmError):
62 def __init__(self, position):
63 message = f"Unable to compute the FWHM of the science Psf at {position}"
64 super().__init__(message)
65 self.position = position
66
67 def metadata(self) -> dict:
68 return {
69 "position": self.position,
70 }
71
72
73class ModelPsfMatchConfig(pexConfig.Config):
74 """Configuration for model-to-model Psf matching"""
75
76 kernel = pexConfig.ConfigChoiceField(
77 doc="kernel type",
78 typemap=dict(
79 AL=PsfMatchConfigAL,
80 ),
81 default="AL",
82 )
83 doAutoPadPsf = pexConfig.Field(
84 dtype=bool,
85 doc=("If too small, automatically pad the science Psf? "
86 "Pad to smallest dimensions appropriate for the matching kernel dimensions, "
87 "as specified by autoPadPsfTo. If false, pad by the padPsfBy config."),
88 default=True,
89 )
90 autoPadPsfTo = pexConfig.RangeField(
91 dtype=float,
92 doc=("Minimum Science Psf dimensions as a fraction of matching kernel dimensions. "
93 "If the dimensions of the Psf to be matched are less than the "
94 "matching kernel dimensions * autoPadPsfTo, pad Science Psf to this size. "
95 "Ignored if doAutoPadPsf=False."),
96 default=1.4,
97 min=1.0,
98 max=2.0
99 )
100 padPsfBy = pexConfig.Field(
101 dtype=int,
102 doc="Pixels (even) to pad Science Psf by before matching. Ignored if doAutoPadPsf=True",
103 default=0,
104 )
105
106 def setDefaults(self):
107 # No sigma clipping
108 self.kernel.active.singleKernelClipping = False
109 self.kernel.active.kernelSumClipping = False
110 self.kernel.active.spatialKernelClipping = False
111 self.kernel.active.checkConditionNumber = False
112
113 # Variance is ill defined
114 self.kernel.active.constantVarianceWeighting = True
115
116 # Do not change specified kernel size
117 self.kernel.active.scaleByFwhm = False
118
119
121 """Match two model Psfs, and application of the Psf-matching kernel
122 to an input Exposure.
123 """
124 ConfigClass = ModelPsfMatchConfig
125
126 def __init__(self, *args, **kwargs):
127 PsfMatchTask.__init__(self, *args, **kwargs)
128 self.kConfig = self.config.kernel.active
129
130 @timeMethod
131 def run(self, exposure, referencePsfModel, kernelSum=1.0):
132 """Psf-match an exposure to a model Psf.
133
134 Parameters
135 ----------
136 exposure : `lsst.afw.image.Exposure`
137 Exposure to Psf-match to the reference Psf model;
138 it must return a valid PSF model via exposure.getPsf()
139 referencePsfModel : `lsst.afw.detection.Psf`
140 The Psf model to match to
141 kernelSum : `float`, optional
142 A multipicative factor to apply to the kernel sum (default=1.0)
143
144 Returns
145 -------
146 result : `struct`
147 - ``psfMatchedExposure`` : the Psf-matched Exposure.
148 This has the same parent bbox, Wcs, PhotoCalib and
149 Filter as the input Exposure but no Psf.
150 In theory the Psf should equal referencePsfModel but
151 the match is likely not exact.
152 - ``psfMatchingKernel`` : the spatially varying Psf-matching kernel
153 - ``kernelCellSet`` : SpatialCellSet used to solve for the Psf-matching kernel
154 - ``referencePsfModel`` : Validated and/or modified reference model used
155
156 Raises
157 ------
158 RuntimeError
159 if the Exposure does not contain a Psf model
160 """
161 if not exposure.hasPsf():
162 raise RuntimeError("exposure does not contain a Psf model")
163
164 maskedImage = exposure.getMaskedImage()
165
166 self.log.info("compute Psf-matching kernel")
167 result = self._buildCellSet(exposure, referencePsfModel)
168 kernelCellSet = result.kernelCellSet
169 referencePsfModel = result.referencePsfModel
170 # TODO: This should be evaluated at (or close to) the center of the
171 # exposure's bounding box in DM-32756.
172 sciAvgPos = exposure.getPsf().getAveragePosition()
173 modelAvgPos = referencePsfModel.getAveragePosition()
174 try:
175 fwhmScience = exposure.getPsf().computeShape(sciAvgPos).getDeterminantRadius()*sigma2fwhm
176 except pexExceptions.RangeError:
178 f"Unable to compute the FWHM of the science Psf at {sciAvgPos}"
179 "due to an unexpectedly large transform."
180 )
181 except pexExceptions.Exception:
182 raise PsfComputeShapeError(sciAvgPos)
183 fwhmModel = referencePsfModel.computeShape(modelAvgPos).getDeterminantRadius()*sigma2fwhm
184
185 basisList = makeKernelBasisList(self.kConfig, fwhmScience, fwhmModel, metadata=self.metadata)
186 spatialSolution, psfMatchingKernel, backgroundModel = self._solve(kernelCellSet, basisList)
187
188 if psfMatchingKernel.isSpatiallyVarying():
189 sParameters = np.array(psfMatchingKernel.getSpatialParameters())
190 sParameters[0][0] = kernelSum
191 psfMatchingKernel.setSpatialParameters(sParameters)
192 else:
193 kParameters = np.array(psfMatchingKernel.getKernelParameters())
194 kParameters[0] = kernelSum
195 psfMatchingKernel.setKernelParameters(kParameters)
196
197 self.log.info("Psf-match science exposure to reference")
198 psfMatchedExposure = afwImage.ExposureF(exposure.getBBox(), exposure.getWcs())
199 psfMatchedExposure.info.id = exposure.info.id
200 psfMatchedExposure.setFilter(exposure.getFilter())
201 psfMatchedExposure.setPhotoCalib(exposure.getPhotoCalib())
202 psfMatchedExposure.getInfo().setVisitInfo(exposure.getInfo().getVisitInfo())
203 psfMatchedExposure.setPsf(referencePsfModel)
204 psfMatchedMaskedImage = psfMatchedExposure.getMaskedImage()
205
206 # Normalize the psf-matching kernel while convolving since its magnitude is meaningless
207 # when PSF-matching one model to another.
208 convolutionControl = afwMath.ConvolutionControl()
209 convolutionControl.setDoNormalize(True)
210 afwMath.convolve(psfMatchedMaskedImage, maskedImage, psfMatchingKernel, convolutionControl)
211
212 self.log.info("done")
213 return pipeBase.Struct(psfMatchedExposure=psfMatchedExposure,
214 psfMatchingKernel=psfMatchingKernel,
215 kernelCellSet=kernelCellSet,
216 metadata=self.metadata,
217 )
218
219 def _diagnostic(self, kernelCellSet, spatialSolution, spatialKernel, spatialBg):
220 """Print diagnostic information on spatial kernel and background fit
221
222 The debugging diagnostics are not really useful here, since the images we are matching have
223 no variance. Thus override the _diagnostic method to generate no logging information"""
224 return
225
226 def _buildCellSet(self, exposure, referencePsfModel):
227 """Build a SpatialCellSet for use with the solve method
228
229 Parameters
230 ----------
231 exposure : `lsst.afw.image.Exposure`
232 The science exposure that will be convolved; must contain a Psf
233 referencePsfModel : `lsst.afw.detection.Psf`
234 Psf model to match to
235
236 Returns
237 -------
238 result : `struct`
239 - ``kernelCellSet`` : a SpatialCellSet to be used by self._solve
240 - ``referencePsfModel`` : Validated and/or modified
241 reference model used to populate the SpatialCellSet
242
243 Notes
244 -----
245 If the reference Psf model and science Psf model have different dimensions,
246 adjust the referencePsfModel (the model to which the exposure PSF will be matched)
247 to match that of the science Psf. If the science Psf dimensions vary across the image,
248 as is common with a WarpedPsf, either pad or clip (depending on config.padPsf)
249 the dimensions to be constant.
250 """
251 sizeCellX = self.kConfig.sizeCellX
252 sizeCellY = self.kConfig.sizeCellY
253
254 scienceBBox = exposure.getBBox()
255 # Extend for proper spatial matching kernel all the way to edge, especially for narrow strips
256 scienceBBox.grow(geom.Extent2I(sizeCellX, sizeCellY))
257
258 sciencePsfModel = exposure.getPsf()
259
260 dimenR = referencePsfModel.getLocalKernel(scienceBBox.getCenter()).getDimensions()
261
262 regionSizeX, regionSizeY = scienceBBox.getDimensions()
263 scienceX0, scienceY0 = scienceBBox.getMin()
264
265 kernelCellSet = afwMath.SpatialCellSet(geom.Box2I(scienceBBox), sizeCellX, sizeCellY)
266
267 nCellX = regionSizeX//sizeCellX
268 nCellY = regionSizeY//sizeCellY
269
270 if nCellX == 0 or nCellY == 0:
271 raise ValueError("Exposure dimensions=%s and sizeCell=(%s, %s). Insufficient area to match" %
272 (scienceBBox.getDimensions(), sizeCellX, sizeCellY))
273
274 # Survey the PSF dimensions of the Spatial Cell Set
275 # to identify the minimum enclosed or maximum bounding square BBox.
276 widthList = []
277 heightList = []
278 for row in range(nCellY):
279 posY = sizeCellY*row + sizeCellY//2 + scienceY0
280 for col in range(nCellX):
281 posX = sizeCellX*col + sizeCellX//2 + scienceX0
282 widthS, heightS = sciencePsfModel.computeBBox(geom.Point2D(posX, posY)).getDimensions()
283 widthList.append(widthS)
284 heightList.append(heightS)
285
286 psfSize = max(max(heightList), max(widthList))
287
288 if self.config.doAutoPadPsf:
289 minPsfSize = nextOddInteger(self.kConfig.kernelSize*self.config.autoPadPsfTo)
290 paddingPix = max(0, minPsfSize - psfSize)
291 else:
292 if self.config.padPsfBy % 2 != 0:
293 raise ValueError("Config padPsfBy (%i pixels) must be even number." %
294 self.config.padPsfBy)
295 paddingPix = self.config.padPsfBy
296
297 if paddingPix > 0:
298 self.log.debug("Padding Science PSF from (%d, %d) to (%d, %d) pixels",
299 psfSize, psfSize, paddingPix + psfSize, paddingPix + psfSize)
300 psfSize += paddingPix
301
302 # Check that PSF is larger than the matching kernel
303 maxKernelSize = psfSize - 1
304 if maxKernelSize % 2 == 0:
305 maxKernelSize -= 1
306 if self.kConfig.kernelSize > maxKernelSize:
307 message = """
308 Kernel size (%d) too big to match Psfs of size %d.
309 Please reconfigure by setting one of the following:
310 1) kernel size to <= %d
311 2) doAutoPadPsf=True
312 3) padPsfBy to >= %s
313 """ % (self.kConfig.kernelSize, psfSize,
314 maxKernelSize, self.kConfig.kernelSize - maxKernelSize)
315 raise ValueError(message)
316
317 dimenS = geom.Extent2I(psfSize, psfSize)
318
319 if (dimenR != dimenS):
320 try:
321 referencePsfModel = referencePsfModel.resized(psfSize, psfSize)
322 self.log.info("Adjusted dimensions of reference PSF model from %s to %s", dimenR, dimenS)
323 except Exception as e:
324 self.log.warning("Zero padding or clipping the reference PSF model of type %s and dimensions"
325 " %s to the science Psf dimensions %s because: %s",
326 referencePsfModel.__class__.__name__, dimenR, dimenS, e)
327 dimenR = dimenS
328
329 ps = pexConfig.makePropertySet(self.kConfig)
330 for row in range(nCellY):
331 # place at center of cell
332 posY = sizeCellY*row + sizeCellY//2 + scienceY0
333
334 for col in range(nCellX):
335 # place at center of cell
336 posX = sizeCellX*col + sizeCellX//2 + scienceX0
337
338 getTraceLogger(self.log, 4).debug("Creating Psf candidate at %.1f %.1f", posX, posY)
339
340 # reference kernel image, at location of science subimage
341 referenceMI = self._makePsfMaskedImage(referencePsfModel, posX, posY, dimensions=dimenR)
342
343 # kernel image we are going to convolve
344 scienceMI = self._makePsfMaskedImage(sciencePsfModel, posX, posY, dimensions=dimenR)
345
346 # The image to convolve is the science image, to the reference Psf.
347 kc = diffimLib.makeKernelCandidate(posX, posY, scienceMI, referenceMI, ps)
348 kernelCellSet.insertCandidate(kc)
349
350 import lsstDebug
351 display = lsstDebug.Info(__name__).display
352 displaySpatialCells = lsstDebug.Info(__name__).displaySpatialCells
353 maskTransparency = lsstDebug.Info(__name__).maskTransparency
354 if not maskTransparency:
355 maskTransparency = 0
356 if display:
357 afwDisplay.setDefaultMaskTransparency(maskTransparency)
358 if display and displaySpatialCells:
359 dituils.showKernelSpatialCells(exposure.getMaskedImage(), kernelCellSet,
360 symb="o", ctype=afwDisplay.CYAN, ctypeUnused=afwDisplay.YELLOW,
361 ctypeBad=afwDisplay.RED, size=4, frame=lsstDebug.frame,
362 title="Image to be convolved")
363 lsstDebug.frame += 1
364 return pipeBase.Struct(kernelCellSet=kernelCellSet,
365 referencePsfModel=referencePsfModel,
366 )
367
368 def _makePsfMaskedImage(self, psfModel, posX, posY, dimensions=None):
369 """Return a MaskedImage of the a PSF Model of specified dimensions
370 """
371 rawKernel = psfModel.computeKernelImage(geom.Point2D(posX, posY)).convertF()
372 if dimensions is None:
373 dimensions = rawKernel.getDimensions()
374 if rawKernel.getDimensions() == dimensions:
375 kernelIm = rawKernel
376 else:
377 # make image of proper size
378 kernelIm = afwImage.ImageF(dimensions)
379 bboxToPlace = geom.Box2I(geom.Point2I((dimensions.getX() - rawKernel.getWidth())//2,
380 (dimensions.getY() - rawKernel.getHeight())//2),
381 rawKernel.getDimensions())
382 kernelIm.assign(rawKernel, bboxToPlace)
383
384 kernelMask = afwImage.Mask(dimensions, 0x0)
385 kernelVar = afwImage.ImageF(dimensions, 1.0)
386 return afwImage.MaskedImageF(kernelIm, kernelMask, kernelVar)
Represent a 2-dimensional array of bitmask pixels.
Definition Mask.h:82
Parameters to control convolution.
A collection of SpatialCells covering an entire image.
An integer coordinate rectangle.
Definition Box.h:55
run(self, exposure, referencePsfModel, kernelSum=1.0)
_buildCellSet(self, exposure, referencePsfModel)
_diagnostic(self, kernelCellSet, spatialSolution, spatialKernel, spatialBg)
_makePsfMaskedImage(self, psfModel, posX, posY, dimensions=None)
_solve(self, kernelCellSet, basisList)
Definition psfMatch.py:797
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.