LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
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.pipe.base as pipeBase
31from lsst.utils.logging import getTraceLogger
32from lsst.utils.timer import timeMethod
33from .makeKernelBasisList import makeKernelBasisList
34from .psfMatch import PsfMatchTask, PsfMatchConfigAL
35from . import utils as dituils
36
37__all__ = ("ModelPsfMatchTask", "ModelPsfMatchConfig")
38
39sigma2fwhm = 2.*np.sqrt(2.*np.log(2.))
40
41
43 nextInt = int(np.ceil(x))
44 return nextInt + 1 if nextInt%2 == 0 else nextInt
45
46
47class ModelPsfMatchConfig(pexConfig.Config):
48 """Configuration for model-to-model Psf matching"""
49
50 kernel = pexConfig.ConfigChoiceField(
51 doc="kernel type",
52 typemap=dict(
53 AL=PsfMatchConfigAL,
54 ),
55 default="AL",
56 )
57 doAutoPadPsf = pexConfig.Field(
58 dtype=bool,
59 doc=("If too small, automatically pad the science Psf? "
60 "Pad to smallest dimensions appropriate for the matching kernel dimensions, "
61 "as specified by autoPadPsfTo. If false, pad by the padPsfBy config."),
62 default=True,
63 )
64 autoPadPsfTo = pexConfig.RangeField(
65 dtype=float,
66 doc=("Minimum Science Psf dimensions as a fraction of matching kernel dimensions. "
67 "If the dimensions of the Psf to be matched are less than the "
68 "matching kernel dimensions * autoPadPsfTo, pad Science Psf to this size. "
69 "Ignored if doAutoPadPsf=False."),
70 default=1.4,
71 min=1.0,
72 max=2.0
73 )
74 padPsfBy = pexConfig.Field(
75 dtype=int,
76 doc="Pixels (even) to pad Science Psf by before matching. Ignored if doAutoPadPsf=True",
77 default=0,
78 )
79
80 def setDefaults(self):
81 # No sigma clipping
82 self.kernel.active.singleKernelClipping = False
83 self.kernel.active.kernelSumClipping = False
84 self.kernel.active.spatialKernelClipping = False
85 self.kernel.active.checkConditionNumber = False
86
87 # Variance is ill defined
88 self.kernel.active.constantVarianceWeighting = True
89
90 # Do not change specified kernel size
91 self.kernel.active.scaleByFwhm = False
92
93
95 """Matching of two model Psfs, and application of the Psf-matching kernel to an input Exposure
96
97 """
98 ConfigClass = ModelPsfMatchConfig
99
100 def __init__(self, *args, **kwargs):
101 """Create a ModelPsfMatchTask
102
103 Parameters
104 ----------
105 *args
106 arguments to be passed to lsst.ip.diffim.PsfMatchTask.__init__
107 **kwargs
108 keyword arguments to be passed to lsst.ip.diffim.PsfMatchTask.__init__
109
110 Notes
111 -----
112 Upon initialization, the kernel configuration is defined by self.config.kernel.active. This Task
113 does have a run() method, which is the default way to call the Task.
114 """
115 PsfMatchTask.__init__(self, *args, **kwargs)
116 self.kConfigkConfig = self.config.kernel.active
117
118 @timeMethod
119 def run(self, exposure, referencePsfModel, kernelSum=1.0):
120 """Psf-match an exposure to a model Psf
121
122 Parameters
123 ----------
124 exposure : `lsst.afw.image.Exposure`
125 Exposure to Psf-match to the reference Psf model;
126 it must return a valid PSF model via exposure.getPsf()
127 referencePsfModel : `lsst.afw.detection.Psf`
128 The Psf model to match to
129 kernelSum : `float`, optional
130 A multipicative factor to apply to the kernel sum (default=1.0)
131
132 Returns
133 -------
134 result : `struct`
135 - ``psfMatchedExposure`` : the Psf-matched Exposure.
136 This has the same parent bbox, Wcs, PhotoCalib and
137 Filter as the input Exposure but no Psf.
138 In theory the Psf should equal referencePsfModel but
139 the match is likely not exact.
140 - ``psfMatchingKernel`` : the spatially varying Psf-matching kernel
141 - ``kernelCellSet`` : SpatialCellSet used to solve for the Psf-matching kernel
142 - ``referencePsfModel`` : Validated and/or modified reference model used
143
144 Raises
145 ------
146 RuntimeError
147 if the Exposure does not contain a Psf model
148 """
149 if not exposure.hasPsf():
150 raise RuntimeError("exposure does not contain a Psf model")
151
152 maskedImage = exposure.getMaskedImage()
153
154 self.log.info("compute Psf-matching kernel")
155 result = self._buildCellSet_buildCellSet(exposure, referencePsfModel)
156 kernelCellSet = result.kernelCellSet
157 referencePsfModel = result.referencePsfModel
158 # TODO: This should be evaluated at (or close to) the center of the
159 # exposure's bounding box in DM-32756.
160 sciAvgPos = exposure.getPsf().getAveragePosition()
161 modelAvgPos = referencePsfModel.getAveragePosition()
162 fwhmScience = exposure.getPsf().computeShape(sciAvgPos).getDeterminantRadius()*sigma2fwhm
163 fwhmModel = referencePsfModel.computeShape(modelAvgPos).getDeterminantRadius()*sigma2fwhm
164
165 basisList = makeKernelBasisList(self.kConfigkConfig, fwhmScience, fwhmModel, metadata=self.metadata)
166 spatialSolution, psfMatchingKernel, backgroundModel = self._solve(kernelCellSet, basisList)
167
168 if psfMatchingKernel.isSpatiallyVarying():
169 sParameters = np.array(psfMatchingKernel.getSpatialParameters())
170 sParameters[0][0] = kernelSum
171 psfMatchingKernel.setSpatialParameters(sParameters)
172 else:
173 kParameters = np.array(psfMatchingKernel.getKernelParameters())
174 kParameters[0] = kernelSum
175 psfMatchingKernel.setKernelParameters(kParameters)
176
177 self.log.info("Psf-match science exposure to reference")
178 psfMatchedExposure = afwImage.ExposureF(exposure.getBBox(), exposure.getWcs())
179 psfMatchedExposure.info.id = exposure.info.id
180 psfMatchedExposure.setFilter(exposure.getFilter())
181 psfMatchedExposure.setPhotoCalib(exposure.getPhotoCalib())
182 psfMatchedExposure.getInfo().setVisitInfo(exposure.getInfo().getVisitInfo())
183 psfMatchedExposure.setPsf(referencePsfModel)
184 psfMatchedMaskedImage = psfMatchedExposure.getMaskedImage()
185
186 # Normalize the psf-matching kernel while convolving since its magnitude is meaningless
187 # when PSF-matching one model to another.
188 convolutionControl = afwMath.ConvolutionControl()
189 convolutionControl.setDoNormalize(True)
190 afwMath.convolve(psfMatchedMaskedImage, maskedImage, psfMatchingKernel, convolutionControl)
191
192 self.log.info("done")
193 return pipeBase.Struct(psfMatchedExposure=psfMatchedExposure,
194 psfMatchingKernel=psfMatchingKernel,
195 kernelCellSet=kernelCellSet,
196 metadata=self.metadata,
197 )
198
199 def _diagnostic(self, kernelCellSet, spatialSolution, spatialKernel, spatialBg):
200 """Print diagnostic information on spatial kernel and background fit
201
202 The debugging diagnostics are not really useful here, since the images we are matching have
203 no variance. Thus override the _diagnostic method to generate no logging information"""
204 return
205
206 def _buildCellSet(self, exposure, referencePsfModel):
207 """Build a SpatialCellSet for use with the solve method
208
209 Parameters
210 ----------
211 exposure : `lsst.afw.image.Exposure`
212 The science exposure that will be convolved; must contain a Psf
213 referencePsfModel : `lsst.afw.detection.Psf`
214 Psf model to match to
215
216 Returns
217 -------
218 result : `struct`
219 - ``kernelCellSet`` : a SpatialCellSet to be used by self._solve
220 - ``referencePsfModel`` : Validated and/or modified
221 reference model used to populate the SpatialCellSet
222
223 Notes
224 -----
225 If the reference Psf model and science Psf model have different dimensions,
226 adjust the referencePsfModel (the model to which the exposure PSF will be matched)
227 to match that of the science Psf. If the science Psf dimensions vary across the image,
228 as is common with a WarpedPsf, either pad or clip (depending on config.padPsf)
229 the dimensions to be constant.
230 """
231 sizeCellX = self.kConfig.sizeCellX
232 sizeCellY = self.kConfig.sizeCellY
233
234 scienceBBox = exposure.getBBox()
235 # Extend for proper spatial matching kernel all the way to edge, especially for narrow strips
236 scienceBBox.grow(geom.Extent2I(sizeCellX, sizeCellY))
237
238 sciencePsfModel = exposure.getPsf()
239
240 dimenR = referencePsfModel.getLocalKernel(scienceBBox.getCenter()).getDimensions()
241
242 regionSizeX, regionSizeY = scienceBBox.getDimensions()
243 scienceX0, scienceY0 = scienceBBox.getMin()
244
245 kernelCellSet = afwMath.SpatialCellSet(geom.Box2I(scienceBBox), sizeCellX, sizeCellY)
246
247 nCellX = regionSizeX//sizeCellX
248 nCellY = regionSizeY//sizeCellY
249
250 if nCellX == 0 or nCellY == 0:
251 raise ValueError("Exposure dimensions=%s and sizeCell=(%s, %s). Insufficient area to match" %
252 (scienceBBox.getDimensions(), sizeCellX, sizeCellY))
253
254 # Survey the PSF dimensions of the Spatial Cell Set
255 # to identify the minimum enclosed or maximum bounding square BBox.
256 widthList = []
257 heightList = []
258 for row in range(nCellY):
259 posY = sizeCellY*row + sizeCellY//2 + scienceY0
260 for col in range(nCellX):
261 posX = sizeCellX*col + sizeCellX//2 + scienceX0
262 widthS, heightS = sciencePsfModel.computeBBox(geom.Point2D(posX, posY)).getDimensions()
263 widthList.append(widthS)
264 heightList.append(heightS)
265
266 psfSize = max(max(heightList), max(widthList))
267
268 if self.config.doAutoPadPsf:
269 minPsfSize = nextOddInteger(self.kConfig.kernelSize*self.config.autoPadPsfTo)
270 paddingPix = max(0, minPsfSize - psfSize)
271 else:
272 if self.config.padPsfBy % 2 != 0:
273 raise ValueError("Config padPsfBy (%i pixels) must be even number." %
274 self.config.padPsfBy)
275 paddingPix = self.config.padPsfBy
276
277 if paddingPix > 0:
278 self.log.debug("Padding Science PSF from (%d, %d) to (%d, %d) pixels",
279 psfSize, psfSize, paddingPix + psfSize, paddingPix + psfSize)
280 psfSize += paddingPix
281
282 # Check that PSF is larger than the matching kernel
283 maxKernelSize = psfSize - 1
284 if maxKernelSize % 2 == 0:
285 maxKernelSize -= 1
286 if self.kConfig.kernelSize > maxKernelSize:
287 message = """
288 Kernel size (%d) too big to match Psfs of size %d.
289 Please reconfigure by setting one of the following:
290 1) kernel size to <= %d
291 2) doAutoPadPsf=True
292 3) padPsfBy to >= %s
293 """ % (self.kConfig.kernelSize, psfSize,
294 maxKernelSize, self.kConfig.kernelSize - maxKernelSize)
295 raise ValueError(message)
296
297 dimenS = geom.Extent2I(psfSize, psfSize)
298
299 if (dimenR != dimenS):
300 try:
301 referencePsfModel = referencePsfModel.resized(psfSize, psfSize)
302 self.log.info("Adjusted dimensions of reference PSF model from %s to %s", dimenR, dimenS)
303 except Exception as e:
304 self.log.warning("Zero padding or clipping the reference PSF model of type %s and dimensions"
305 " %s to the science Psf dimensions %s because: %s",
306 referencePsfModel.__class__.__name__, dimenR, dimenS, e)
307 dimenR = dimenS
308
309 ps = pexConfig.makePropertySet(self.kConfig)
310 for row in range(nCellY):
311 # place at center of cell
312 posY = sizeCellY*row + sizeCellY//2 + scienceY0
313
314 for col in range(nCellX):
315 # place at center of cell
316 posX = sizeCellX*col + sizeCellX//2 + scienceX0
317
318 getTraceLogger(self.log, 4).debug("Creating Psf candidate at %.1f %.1f", posX, posY)
319
320 # reference kernel image, at location of science subimage
321 referenceMI = self._makePsfMaskedImage(referencePsfModel, posX, posY, dimensions=dimenR)
322
323 # kernel image we are going to convolve
324 scienceMI = self._makePsfMaskedImage(sciencePsfModel, posX, posY, dimensions=dimenR)
325
326 # The image to convolve is the science image, to the reference Psf.
327 kc = diffimLib.makeKernelCandidate(posX, posY, scienceMI, referenceMI, ps)
328 kernelCellSet.insertCandidate(kc)
329
330 import lsstDebug
331 display = lsstDebug.Info(__name__).display
332 displaySpatialCells = lsstDebug.Info(__name__).displaySpatialCells
333 maskTransparency = lsstDebug.Info(__name__).maskTransparency
334 if not maskTransparency:
335 maskTransparency = 0
336 if display:
337 afwDisplay.setDefaultMaskTransparency(maskTransparency)
338 if display and displaySpatialCells:
339 dituils.showKernelSpatialCells(exposure.getMaskedImage(), kernelCellSet,
340 symb="o", ctype=afwDisplay.CYAN, ctypeUnused=afwDisplay.YELLOW,
341 ctypeBad=afwDisplay.RED, size=4, frame=lsstDebug.frame,
342 title="Image to be convolved")
343 lsstDebug.frame += 1
344 return pipeBase.Struct(kernelCellSet=kernelCellSet,
345 referencePsfModel=referencePsfModel,
346 )
347
348 def _makePsfMaskedImage(self, psfModel, posX, posY, dimensions=None):
349 """Return a MaskedImage of the a PSF Model of specified dimensions
350 """
351 rawKernel = psfModel.computeKernelImage(geom.Point2D(posX, posY)).convertF()
352 if dimensions is None:
353 dimensions = rawKernel.getDimensions()
354 if rawKernel.getDimensions() == dimensions:
355 kernelIm = rawKernel
356 else:
357 # make image of proper size
358 kernelIm = afwImage.ImageF(dimensions)
359 bboxToPlace = geom.Box2I(geom.Point2I((dimensions.getX() - rawKernel.getWidth())//2,
360 (dimensions.getY() - rawKernel.getHeight())//2),
361 rawKernel.getDimensions())
362 kernelIm.assign(rawKernel, bboxToPlace)
363
364 kernelMask = afwImage.Mask(dimensions, 0x0)
365 kernelVar = afwImage.ImageF(dimensions, 1.0)
366 return afwImage.MaskedImageF(kernelIm, kernelMask, kernelVar)
int max
Represent a 2-dimensional array of bitmask pixels.
Definition Mask.h:77
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, returnOnExcept=False)
Definition psfMatch.py:894
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.