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