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
139 displayExposure = display and \
141 displayPsfComponents = display and \
143 showBadCandidates = display and \
145 displayResiduals = display and \
147 displayPsfMosaic = display and \
150 afwDisplay.setDefaultMaskTransparency(75)
151
152
153 mi = exposure.getMaskedImage()
154
155 nCand = len(psfCandidateList)
156 if nCand == 0:
157 raise RuntimeError("No PSF candidates supplied.")
158
159
160
161 if display:
162
163 bbox = mi.getBBox(afwImage.PARENT)
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
193
194
195
196 defaultsFile = os.path.join(os.environ["MEAS_EXTENSIONS_PSFEX_DIR"], "config", "default-lsst.psfex")
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
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
258 xc, yc = source.getX(), source.getY()
259 if not np.isfinite(xc) or not np.isfinite(yc):
260 continue
261
262 if fluxFlagName in source.schema and source.get(fluxFlagName):
263 continue
264
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
277
278
279
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)
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
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
324 psfSet.trimMemory()
325
326
327
328
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
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
362
363
364
365
366
367
368 try:
369 _ = psf.getKernel(psf.getAveragePosition())
371 raise RuntimeError("Failed to determine psfex psf: too few good stars.")
372
373
374
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
397
398
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
A collection of SpatialCells covering an entire image.
Class for storing generic metadata.
Reports invalid arguments.
Reports attempts to exceed implementation-defined length limits for some classes.
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)