39 spatialOrder = pexConfig.Field(
40 doc=
"specify spatial order for PSF kernel creation",
43 check=
lambda x: x >= 1,
45 sizeCellX = pexConfig.Field(
46 doc=
"size of cell used to determine PSF (pixels, column direction)",
50 check=
lambda x: x >= 10,
52 sizeCellY = pexConfig.Field(
53 doc=
"size of cell used to determine PSF (pixels, row direction)",
55 default=sizeCellX.default,
57 check=
lambda x: x >= 10,
59 samplingSize = pexConfig.Field(
60 doc=
"Resolution of the internal PSF model relative to the pixel size; "
61 "e.g. 0.5 is equal to 2x oversampling",
65 badMaskBits = pexConfig.ListField(
66 doc=
"List of mask bits which cause a source to be rejected as bad "
67 "N.b. INTRP is used specially in PsfCandidateSet; it means \"Contaminated by neighbour\"",
69 default=[
"INTRP",
"SAT"],
71 psfexBasis = pexConfig.ChoiceField(
72 doc=
"BASIS value given to psfex. PIXEL_AUTO will use the requested samplingSize only if "
73 "the FWHM < 3 pixels. Otherwise, it will use samplingSize=1. PIXEL will always use the "
74 "requested samplingSize",
77 "PIXEL":
"Always use requested samplingSize",
78 "PIXEL_AUTO":
"Only use requested samplingSize when FWHM < 3",
83 tolerance = pexConfig.Field(
84 doc=
"tolerance of spatial fitting",
88 lam = pexConfig.Field(
89 doc=
"floor for variance is lam*data",
93 reducedChi2ForPsfCandidates = pexConfig.Field(
94 doc=
"for psf candidate evaluation",
98 spatialReject = pexConfig.Field(
99 doc=
"Rejection threshold (stdev) for candidates based on spatial fit",
103 recentroid = pexConfig.Field(
104 doc=
"Should PSFEX be permitted to recentroid PSF candidates?",
108 kernelSize = pexConfig.Field(
109 doc=(
"Size of the postage stamp around each star that is extracted for fitting."
110 "Note: this reflects the oversampling setting of the psf, set by `samplingSize`;"
111 "e.g. `samplingSize=0.5` would require this value to be 2x what you expect."),
118 ConfigClass = PsfexPsfDeterminerConfig
120 def determinePsf(self, exposure, psfCandidateList, metadata=None, flagKey=None):
121 """Determine a PSFEX PSF model for an exposure given a list of PSF
126 exposure: `lsst.afw.image.Exposure`
127 Exposure containing the PSF candidates.
128 psfCandidateList: iterable of `lsst.meas.algorithms.PsfCandidate`
129 Sequence of PSF candidates typically obtained by detecting sources
130 and then running them through a star selector.
131 metadata: metadata, optional
132 A home for interesting tidbits of information.
133 flagKey: `lsst.afw.table.Key`, optional
134 Schema key used to mark sources actually used in PSF determination.
138 psf: `lsst.meas.extensions.psfex.PsfexPsf`
144 displayExposure = display
and \
146 displayPsfComponents = display
and \
148 showBadCandidates = display
and \
150 displayResiduals = display
and \
152 displayPsfMosaic = display
and \
155 afwDisplay.setDefaultMaskTransparency(75)
158 mi = exposure.getMaskedImage()
160 nCand = len(psfCandidateList)
162 raise RuntimeError(
"No PSF candidates supplied.")
168 bbox = mi.getBBox(afwImage.PARENT)
173 sizes = np.empty(nCand)
174 for i, psfCandidate
in enumerate(psfCandidateList):
177 psfCellSet.insertCandidate(psfCandidate)
178 except Exception
as e:
179 self.log.
error(
"Skipping PSF candidate %d of %d: %s", i, len(psfCandidateList), e)
182 source = psfCandidate.getSource()
183 quad = afwEll.Quadrupole(source.getIxx(), source.getIyy(), source.getIxy())
184 rmsSize = quad.getTraceRadius()
187 if self.config.kernelSize >= 15:
188 self.log.
warn(
"NOT scaling kernelSize by stellar quadrupole moment, but using absolute value")
189 actualKernelSize = self.config.kernelSize
191 actualKernelSize = 2 * int(self.config.kernelSize * np.sqrt(np.median(sizes)) + 0.5) + 1
192 if actualKernelSize < self.config.kernelSizeMin:
193 actualKernelSize = self.config.kernelSizeMin
194 if actualKernelSize > self.config.kernelSizeMax:
195 actualKernelSize = self.config.kernelSizeMax
197 rms = np.median(sizes)
198 msg =
"Median PSF RMS size=%.2f pixels (\"FWHM\"=%.2f)" % (rms, 2*np.sqrt(2*np.log(2))*rms)
203 pixKernelSize = actualKernelSize
204 if self.config.samplingSize > 0:
205 pixKernelSize = int(actualKernelSize*self.config.samplingSize)
206 if pixKernelSize % 2 == 0:
208 self.log.
trace(
"Psfex Kernel size=%.2f, Image Kernel Size=%.2f", actualKernelSize, pixKernelSize)
209 psfCandidateList[0].setHeight(pixKernelSize)
210 psfCandidateList[0].setWidth(pixKernelSize)
216 defaultsFile = os.path.join(os.environ[
"MEAS_EXTENSIONS_PSFEX_DIR"],
"config",
"default-lsst.psfex")
218 args_md.set(
"BASIS_TYPE", str(self.config.psfexBasis))
219 args_md.set(
"PSFVAR_DEGREES", str(self.config.spatialOrder))
220 args_md.set(
"PSF_SIZE", str(actualKernelSize))
221 args_md.set(
"PSF_SAMPLING", str(self.config.samplingSize))
222 prefs = psfex.Prefs(defaultsFile, args_md)
223 prefs.setCommandLine([])
224 prefs.addCatalog(
"psfexPsfDeterminer")
227 principalComponentExclusionFlag = bool(bool(psfex.Context.REMOVEHIDDEN)
228 if False else psfex.Context.KEEPHIDDEN)
229 context = psfex.Context(prefs.getContextName(), prefs.getContextGroup(),
230 prefs.getGroupDeg(), principalComponentExclusionFlag)
231 psfSet = psfex.Set(context)
232 psfSet.setVigSize(pixKernelSize, pixKernelSize)
233 psfSet.setFwhm(2*np.sqrt(2*np.log(2))*np.median(sizes))
234 psfSet.setRecentroid(self.config.recentroid)
238 ccd = exposure.getDetector()
240 gain = np.mean(np.array([a.getGain()
for a
in ccd]))
243 self.log.
warn(
"Setting gain to %g" % (gain,))
246 for i, key
in enumerate(context.getName()):
249 contextvalp.append(exposure.getMetadata().getScalar(key[1:]))
250 except KeyError
as e:
251 raise RuntimeError(
"%s parameter not found in the header of %s" %
252 (key[1:], prefs.getContextName()))
from e
255 contextvalp.append(np.array([psfCandidateList[_].getSource().get(key)
256 for _
in range(nCand)]))
257 except KeyError
as e:
258 raise RuntimeError(
"%s parameter not found" % (key,))
from e
259 psfSet.setContextname(i, key)
264 disp = afwDisplay.Display(frame=frame)
265 disp.mtv(exposure, title=
"psf determination")
267 badBits = mi.getMask().getPlaneBitMask(self.config.badMaskBits)
268 fluxName = prefs.getPhotfluxRkey()
269 fluxFlagName =
"base_" + fluxName +
"_flag"
272 for i, psfCandidate
in enumerate(psfCandidateList):
273 source = psfCandidate.getSource()
276 xc, yc = source.getX(), source.getY()
277 if not np.isfinite(xc)
or not np.isfinite(yc):
280 if fluxFlagName
in source.schema
and source.get(fluxFlagName):
283 flux = source.get(fluxName)
284 if flux < 0
or not np.isfinite(flux):
288 pstamp = psfCandidate.getMaskedImage().
clone()
299 sample = psfSet.newSample()
300 sample.setCatindex(catindex)
301 sample.setExtindex(ext)
302 sample.setObjindex(i)
304 imArray = pstamp.getImage().getArray()
305 imArray[np.where(np.bitwise_and(pstamp.getMask().getArray(), badBits))] = \
307 sample.setVig(imArray)
310 sample.setBacknoise2(backnoise2)
314 sample.setFluxrad(sizes[i])
316 for j
in range(psfSet.getNcontext()):
317 sample.setContext(j, float(contextvalp[j][i]))
318 except Exception
as e:
319 self.log.
error(
"Exception when processing sample at (%f,%f): %s", xc, yc, e)
322 psfSet.finiSample(sample)
328 with disp.Buffering():
329 disp.dot(
"o", xc, yc, ctype=afwDisplay.CYAN, size=4)
331 if psfSet.getNsample() == 0:
332 raise RuntimeError(
"No good PSF candidates to pass to PSFEx")
335 for i
in range(psfSet.getNcontext()):
336 cmin = contextvalp[i].
min()
337 cmax = contextvalp[i].
max()
338 psfSet.setContextScale(i, cmax - cmin)
339 psfSet.setContextOffset(i, (cmin + cmax)/2.0)
349 field = psfex.Field(
"Unknown")
350 field.addExt(exposure.getWcs(), exposure.getWidth(), exposure.getHeight(), psfSet.getNsample())
358 psfex.makeit(fields, sets)
359 psfs = field.getPsfs()
363 for i
in range(sets[0].getNsample()):
364 index = sets[0].getSample(i).getObjindex()
366 good_indices.append(index)
368 if flagKey
is not None:
369 for i, psfCandidate
in enumerate(psfCandidateList):
370 source = psfCandidate.getSource()
371 if i
in good_indices:
372 source.set(flagKey,
True)
374 xpos = np.array(xpos)
375 ypos = np.array(ypos)
376 numGoodStars = len(good_indices)
377 avgX, avgY = np.mean(xpos), np.mean(ypos)
385 assert psfCellSet
is not None
388 maUtils.showPsfSpatialCells(exposure, psfCellSet, showChi2=
True,
389 symb=
"o", ctype=afwDisplay.YELLOW, ctypeBad=afwDisplay.RED,
390 size=8, display=disp)
392 disp4 = afwDisplay.Display(frame=4)
393 maUtils.showPsfCandidates(exposure, psfCellSet, psf=psf, display=disp4,
394 normalize=normalizeResiduals,
395 showBadCandidates=showBadCandidates)
396 if displayPsfComponents:
397 disp6 = afwDisplay.Display(frame=6)
398 maUtils.showPsf(psf, display=disp6)
400 disp7 = afwDisplay.Display(frame=7)
401 maUtils.showPsfMosaic(exposure, psf, display=disp7, showFwhm=
True)
402 disp.scale(
'linear', 0, 1)
408 if metadata
is not None:
409 metadata.set(
"spatialFitChi2", np.nan)
410 metadata.set(
"numAvailStars", nCand)
411 metadata.set(
"numGoodStars", numGoodStars)
412 metadata.set(
"avgX", avgX)
413 metadata.set(
"avgY", avgY)
415 return psf, psfCellSet
418 measAlg.psfDeterminerRegistry.register(
"psfex", PsfexPsfDeterminerTask)
A collection of SpatialCells covering an entire image.
Class for storing generic metadata.
Represent a PSF as a linear combination of PSFEX (== Karhunen-Loeve) basis functions.
def determinePsf(self, exposure, psfCandidateList, metadata=None, flagKey=None)
Reports attempts to exceed implementation-defined length limits for some classes.
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
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)
Fit spatial kernel using approximate fluxes for candidates, and solving a linear system of equations.