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
psfexPsfDeterminer.py
Go to the documentation of this file.
1# This file is part of meas_extensions_psfex.
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
22__all__ = ("PsfexPsfDeterminerConfig", "PsfexPsfDeterminerTask")
23
24import os
25import numpy as np
26
27import lsst.daf.base as dafBase
28import lsst.pex.config as pexConfig
29import lsst.pex.exceptions as pexExcept
30import lsst.geom as geom
31import lsst.afw.geom.ellipses as afwEll
32import lsst.afw.display as afwDisplay
33import lsst.afw.image as afwImage
34import lsst.afw.math as afwMath
35import lsst.meas.algorithms as measAlg
36import lsst.meas.algorithms.utils as maUtils
37import lsst.meas.extensions.psfex as psfex
38
39
40class PsfexPsfDeterminerConfig(measAlg.BasePsfDeterminerConfig):
41 spatialOrder = pexConfig.Field[int](
42 doc="specify spatial order for PSF kernel creation",
43 default=2,
44 check=lambda x: x >= 1,
45 )
46 sizeCellX = pexConfig.Field[int](
47 doc="size of cell used to determine PSF (pixels, column direction)",
48 default=256,
49 # minValue = 10,
50 check=lambda x: x >= 10,
51 )
52 sizeCellY = pexConfig.Field[int](
53 doc="size of cell used to determine PSF (pixels, row direction)",
54 default=sizeCellX.default,
55 # minValue = 10,
56 check=lambda x: x >= 10,
57 )
58 samplingSize = pexConfig.Field[float](
59 doc="Resolution of the internal PSF model relative to the pixel size; "
60 "e.g. 0.5 is equal to 2x oversampling",
61 default=0.5,
62 )
63 badMaskBits = pexConfig.ListField[str](
64 doc="List of mask bits which cause a source to be rejected as bad "
65 "N.b. INTRP is used specially in PsfCandidateSet; it means \"Contaminated by neighbour\"",
66 default=["INTRP", "SAT"],
67 )
68 psfexBasis = pexConfig.ChoiceField[str](
69 doc="BASIS value given to psfex. PIXEL_AUTO will use the requested samplingSize only if "
70 "the FWHM < 3 pixels. Otherwise, it will use samplingSize=1. PIXEL will always use the "
71 "requested samplingSize",
72 allowed={
73 "PIXEL": "Always use requested samplingSize",
74 "PIXEL_AUTO": "Only use requested samplingSize when FWHM < 3",
75 },
76 default='PIXEL_AUTO',
77 optional=False,
78 )
79 tolerance = pexConfig.Field[float](
80 doc="tolerance of spatial fitting",
81 default=1e-2,
82 )
83 lam = pexConfig.Field[float](
84 doc="floor for variance is lam*data",
85 default=0.05,
86 )
87 reducedChi2ForPsfCandidates = pexConfig.Field[float](
88 doc="for psf candidate evaluation",
89 default=2.0,
90 )
91 spatialReject = pexConfig.Field[float](
92 doc="Rejection threshold (stdev) for candidates based on spatial fit",
93 default=3.0,
94 )
95 recentroid = pexConfig.Field[bool](
96 doc="Should PSFEX be permitted to recentroid PSF candidates?",
97 default=False,
98 )
99 photometricFluxField = pexConfig.Field[str](
100 doc="Flux field to use for photometric normalization. This overrides the "
101 "``PHOTFLUX_KEY`` field for psfex. The associated flux error is "
102 "derived by appending ``Err`` to this field.",
103 default="base_CircularApertureFlux_9_0_instFlux",
104 )
105
106 def setDefaults(self):
107 super().setDefaults()
108 self.stampSize = 41
109
110
111class PsfexPsfDeterminerTask(measAlg.BasePsfDeterminerTask):
112 ConfigClass = PsfexPsfDeterminerConfig
113
114 def determinePsf(self, exposure, psfCandidateList, metadata=None, flagKey=None):
115 """Determine a PSFEX PSF model for an exposure given a list of PSF
116 candidates.
117
118 Parameters
119 ----------
120 exposure: `lsst.afw.image.Exposure`
121 Exposure containing the PSF candidates.
122 psfCandidateList: iterable of `lsst.meas.algorithms.PsfCandidate`
123 Sequence of PSF candidates typically obtained by detecting sources
124 and then running them through a star selector.
125 metadata: metadata, optional
126 A home for interesting tidbits of information.
127 flagKey: `lsst.afw.table.Key`, optional
128 Schema key used to mark sources actually used in PSF determination.
129
130 Returns
131 -------
132 psf: `lsst.meas.extensions.psfex.PsfexPsf`
133 The determined PSF.
134 """
135 psfCandidateList = self.downsampleCandidates(psfCandidateList)
136
137 import lsstDebug
138 display = lsstDebug.Info(__name__).display
139 displayExposure = display and \
140 lsstDebug.Info(__name__).displayExposure # display the Exposure + spatialCells
141 displayPsfComponents = display and \
142 lsstDebug.Info(__name__).displayPsfComponents # show the basis functions
143 showBadCandidates = display and \
144 lsstDebug.Info(__name__).showBadCandidates # Include bad candidates (meaningless, methinks)
145 displayResiduals = display and \
146 lsstDebug.Info(__name__).displayResiduals # show residuals
147 displayPsfMosaic = display and \
148 lsstDebug.Info(__name__).displayPsfMosaic # show mosaic of reconstructed PSF(x,y)
149 normalizeResiduals = lsstDebug.Info(__name__).normalizeResiduals
150 afwDisplay.setDefaultMaskTransparency(75)
151 # Normalise residuals by object amplitude
152
153 mi = exposure.getMaskedImage()
154
155 nCand = len(psfCandidateList)
156 if nCand == 0:
157 raise RuntimeError("No PSF candidates supplied.")
158 #
159 # How big should our PSF models be?
160 #
161 if display: # only needed for debug plots
162 # construct and populate a spatial cell set
163 bbox = mi.getBBox(afwImage.PARENT)
164 psfCellSet = afwMath.SpatialCellSet(bbox, self.config.sizeCellX, self.config.sizeCellY)
165 else:
166 psfCellSet = None
167
168 sizes = np.empty(nCand)
169 for i, psfCandidate in enumerate(psfCandidateList):
170 try:
171 if psfCellSet:
172 psfCellSet.insertCandidate(psfCandidate)
173 except Exception as e:
174 self.log.error("Skipping PSF candidate %d of %d: %s", i, len(psfCandidateList), e)
175 continue
176
177 source = psfCandidate.getSource()
178 quad = afwEll.Quadrupole(source.getIxx(), source.getIyy(), source.getIxy())
179 rmsSize = quad.getTraceRadius()
180 sizes[i] = rmsSize
181
182 pixKernelSize = self.config.stampSize
183 actualKernelSize = int(2*np.floor(0.5*pixKernelSize/self.config.samplingSize) + 1)
184
185 if display:
186 rms = np.median(sizes)
187 self.log.debug("Median PSF RMS size=%.2f pixels (\"FWHM\"=%.2f)",
188 rms, 2*np.sqrt(2*np.log(2))*rms)
189
190 self.log.trace("Psfex Kernel size=%.2f, Image Kernel Size=%.2f", actualKernelSize, pixKernelSize)
191
192 # -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- BEGIN PSFEX
193 #
194 # Insert the good candidates into the set
195 #
196 defaultsFile = os.path.join(os.environ["MEAS_EXTENSIONS_PSFEX_DIR"], "config", "default-lsst.psfex")
197 args_md = dafBase.PropertySet()
198 args_md.set("BASIS_TYPE", str(self.config.psfexBasis))
199 args_md.set("PSFVAR_DEGREES", str(self.config.spatialOrder))
200 args_md.set("PSF_SIZE", str(actualKernelSize))
201 args_md.set("PSF_SAMPLING", str(self.config.samplingSize))
202 args_md.set("PHOTFLUX_KEY", str(self.config.photometricFluxField))
203 args_md.set("PHOTFLUXERR_KEY", str(self.config.photometricFluxField) + "Err")
204 prefs = psfex.Prefs(defaultsFile, args_md)
205 prefs.setCommandLine([])
206 prefs.addCatalog("psfexPsfDeterminer")
207
208 prefs.use()
209 principalComponentExclusionFlag = bool(bool(psfex.Context.REMOVEHIDDEN)
210 if False else psfex.Context.KEEPHIDDEN)
211 context = psfex.Context(prefs.getContextName(), prefs.getContextGroup(),
212 prefs.getGroupDeg(), principalComponentExclusionFlag)
213 psfSet = psfex.Set(context)
214 psfSet.setVigSize(pixKernelSize, pixKernelSize)
215 psfSet.setFwhm(2*np.sqrt(2*np.log(2))*np.median(sizes))
216 psfSet.setRecentroid(self.config.recentroid)
217
218 catindex, ext = 0, 0
219 backnoise2 = afwMath.makeStatistics(mi.getImage(), afwMath.VARIANCECLIP).getValue()
220 ccd = exposure.getDetector()
221 if ccd:
222 gain = np.mean(np.array([a.getGain() for a in ccd]))
223 else:
224 gain = 1.0
225 self.log.warning("Setting gain to %g", gain)
226
227 contextvalp = []
228 for i, key in enumerate(context.getName()):
229 if key[0] == ':':
230 try:
231 contextvalp.append(exposure.getMetadata().getScalar(key[1:]))
232 except KeyError as e:
233 raise RuntimeError("%s parameter not found in the header of %s" %
234 (key[1:], prefs.getContextName())) from e
235 else:
236 try:
237 contextvalp.append(np.array([psfCandidateList[_].getSource().get(key)
238 for _ in range(nCand)]))
239 except KeyError as e:
240 raise RuntimeError("%s parameter not found" % (key,)) from e
241 psfSet.setContextname(i, key)
242
243 if display:
244 frame = 0
245 if displayExposure:
246 disp = afwDisplay.Display(frame=frame)
247 disp.mtv(exposure, title="psf determination")
248
249 badBits = mi.getMask().getPlaneBitMask(self.config.badMaskBits)
250 fluxName = prefs.getPhotfluxRkey()
251 fluxFlagName = "base_" + fluxName + "_flag"
252
253 xpos, ypos = [], []
254 for i, psfCandidate in enumerate(psfCandidateList):
255 source = psfCandidate.getSource()
256
257 # skip sources with bad centroids
258 xc, yc = source.getX(), source.getY()
259 if not np.isfinite(xc) or not np.isfinite(yc):
260 continue
261 # skip flagged sources
262 if fluxFlagName in source.schema and source.get(fluxFlagName):
263 continue
264 # skip nonfinite and negative sources
265 flux = source.get(fluxName)
266 if flux < 0 or not np.isfinite(flux):
267 continue
268
269 try:
270 pstamp = psfCandidate.getMaskedImage(pixKernelSize, pixKernelSize).clone()
272 self.log.warning("Could not get stamp image for psfCandidate: %s with kernel size: %s",
273 psfCandidate, pixKernelSize)
274 continue
275
276 # From this point, we're configuring the "sample" (PSFEx's version
277 # of a PSF candidate).
278 # Having created the sample, we must proceed to configure it, and
279 # then fini (finalize), or it will be malformed.
280 try:
281 sample = psfSet.newSample()
282 sample.setCatindex(catindex)
283 sample.setExtindex(ext)
284 sample.setObjindex(i)
285
286 imArray = pstamp.getImage().getArray()
287 imArray[np.where(np.bitwise_and(pstamp.getMask().getArray(), badBits))] = \
288 -2*psfex.BIG
289 sample.setVig(imArray)
290
291 sample.setNorm(flux)
292 sample.setBacknoise2(backnoise2)
293 sample.setGain(gain)
294 sample.setX(xc)
295 sample.setY(yc)
296 sample.setFluxrad(sizes[i])
297
298 for j in range(psfSet.getNcontext()):
299 sample.setContext(j, float(contextvalp[j][i]))
300 except Exception as e:
301 self.log.error("Exception when processing sample at (%f,%f): %s", xc, yc, e)
302 continue
303 else:
304 psfSet.finiSample(sample)
305
306 xpos.append(xc) # for QA
307 ypos.append(yc)
308
309 if displayExposure:
310 with disp.Buffering():
311 disp.dot("o", xc, yc, ctype=afwDisplay.CYAN, size=4)
312
313 if psfSet.getNsample() == 0:
314 raise RuntimeError("No good PSF candidates to pass to PSFEx")
315
316 # ---- Update min and max and then the scaling
317 for i in range(psfSet.getNcontext()):
318 cmin = contextvalp[i].min()
319 cmax = contextvalp[i].max()
320 psfSet.setContextScale(i, cmax - cmin)
321 psfSet.setContextOffset(i, (cmin + cmax)/2.0)
322
323 # Don't waste memory!
324 psfSet.trimMemory()
325
326 # =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- END PSFEX
327 #
328 # Do a PSFEX decomposition of those PSF candidates
329 #
330 fields = []
331 field = psfex.Field("Unknown")
332 field.addExt(exposure.getWcs(), exposure.getWidth(), exposure.getHeight(), psfSet.getNsample())
333 field.finalize()
334
335 fields.append(field)
336
337 sets = []
338 sets.append(psfSet)
339
340 psfex.makeit(fields, sets)
341 psfs = field.getPsfs()
342
343 # Flag which objects were actually used in psfex by
344 good_indices = []
345 for i in range(sets[0].getNsample()):
346 index = sets[0].getSample(i).getObjindex()
347 if index > -1:
348 good_indices.append(index)
349
350 if flagKey is not None:
351 for i, psfCandidate in enumerate(psfCandidateList):
352 source = psfCandidate.getSource()
353 if i in good_indices:
354 source.set(flagKey, True)
355
356 xpos = np.array(xpos)
357 ypos = np.array(ypos)
358 numGoodStars = len(good_indices)
359 avgX, avgY = np.mean(xpos), np.mean(ypos)
360
361 psf = psfex.PsfexPsf(psfs[0], geom.Point2D(avgX, avgY))
362
363 # If there are too few stars, the PSFEx psf model will reduce the order
364 # to 0, which the Science Pipelines code cannot handle (see
365 # https://github.com/lsst/meas_extensions_psfex/blob/f0d5218b5446faf5e39edc30e31d2e6f673ef294/src/PsfexPsf.cc#L118
366 # ). The easiest way to test for this condition is trying to compute
367 # the PSF kernel and checking for an InvalidParameterError.
368 try:
369 _ = psf.getKernel(psf.getAveragePosition())
371 raise RuntimeError("Failed to determine psfex psf: too few good stars.")
372
373 #
374 # Display code for debugging
375 #
376 if display:
377 assert psfCellSet is not None
378
379 if displayExposure:
380 maUtils.showPsfSpatialCells(exposure, psfCellSet, showChi2=True,
381 symb="o", ctype=afwDisplay.YELLOW, ctypeBad=afwDisplay.RED,
382 size=8, display=disp)
383 if displayResiduals:
384 disp4 = afwDisplay.Display(frame=4)
385 maUtils.showPsfCandidates(exposure, psfCellSet, psf=psf, display=disp4,
386 normalize=normalizeResiduals,
387 showBadCandidates=showBadCandidates)
388 if displayPsfComponents:
389 disp6 = afwDisplay.Display(frame=6)
390 maUtils.showPsf(psf, display=disp6)
391 if displayPsfMosaic:
392 disp7 = afwDisplay.Display(frame=7)
393 maUtils.showPsfMosaic(exposure, psf, display=disp7, showFwhm=True)
394 disp.scale('linear', 0, 1)
395 #
396 # Generate some QA information
397 #
398 # Count PSF stars
399 #
400 if metadata is not None:
401 metadata["spatialFitChi2"] = np.nan
402 metadata["numAvailStars"] = nCand
403 metadata["numGoodStars"] = numGoodStars
404 metadata["avgX"] = avgX
405 metadata["avgY"] = avgY
406
407 return psf, psfCellSet
408
409
410measAlg.psfDeterminerRegistry.register("psfex", PsfexPsfDeterminerTask)
int min
int max
A collection of SpatialCells covering an entire image.
Class for storing generic metadata.
Definition PropertySet.h:66
Represent a PSF as a linear combination of PSFEX (== Karhunen-Loeve) basis functions.
Definition PsfexPsf.h:40
determinePsf(self, exposure, psfCandidateList, metadata=None, flagKey=None)
Reports invalid arguments.
Definition Runtime.h:66
Reports attempts to exceed implementation-defined length limits for some classes.
Definition Runtime.h:76
Statistics makeStatistics(lsst::afw::image::Image< Pixel > const &img, lsst::afw::image::Mask< image::MaskPixel > const &msk, int const flags, StatisticsControl const &sctrl=StatisticsControl())
Handle a watered-down front-end to the constructor (no variance)
Definition Statistics.h:361