138 fitBasisComponents=False, variance=None, chi=None):
139 """Display the PSF candidates.
140
141 If psf is provided include PSF model and residuals; if normalize is true normalize the PSFs
142 (and residuals)
143
144 If chi is True, generate a plot of residuals/sqrt(variance), i.e. chi
145
146 If fitBasisComponents is true, also find the best linear combination of the PSF's components
147 (if they exist)
148 """
149 if not display:
150 display = afwDisplay.Display()
151
152 if chi is None:
153 if variance is not None:
154 chi = variance
155
156
157
158 mos = displayUtils.Mosaic()
159
160 candidateCenters = []
161 candidateCentersBad = []
162 candidateIndex = 0
163
164 for cell in psfCellSet.getCellList():
165 for cand in cell.begin(False):
166 rchi2 = cand.getChi2()
167 if rchi2 > 1e100:
168 rchi2 = numpy.nan
169
170 if not showBadCandidates and cand.isBad():
171 continue
172
173 if psf:
174 im_resid = displayUtils.Mosaic(gutter=0, background=-5, mode="x")
175
176 try:
177 im = cand.getMaskedImage()
178 xc, yc = cand.getXCenter(), cand.getYCenter()
179
180 margin = 0 if True else 5
181 w, h = im.getDimensions()
183
184 if margin > 0:
185 bim = im.Factory(w + 2*margin, h + 2*margin)
186
189 bim.getVariance().set(stdev**2)
190
191 bim.assign(im, bbox)
192 im = bim
193 xc += margin
194 yc += margin
195
196 im = im.Factory(im, True)
197 im.setXY0(cand.getMaskedImage().getXY0())
198 except Exception:
199 continue
200
201 if not variance:
202 im_resid.append(im.Factory(im, True))
203
204 if True:
205 mi = im
206 psfIm = mi.getImage()
207 config = measBase.SingleFrameMeasurementTask.ConfigClass()
208 config.slots.centroid = "base_SdssCentroid"
209
210 schema = afwTable.SourceTable.makeMinimalSchema()
211 measureSources = measBase.SingleFrameMeasurementTask(schema, config=config)
213
214 extra = 10
215 miBig = mi.Factory(im.getWidth() + 2*extra, im.getHeight() + 2*extra)
216 miBig[extra:-extra, extra:-extra, afwImage.LOCAL] = mi
217 miBig.setXY0(mi.getX0() - extra, mi.getY0() - extra)
218 mi = miBig
219 del miBig
220
222 exp.setPsf(psf)
223
226 "DETECTED")
227 footprintSet.makeSources(catalog)
228
229 if len(catalog) == 0:
230 raise RuntimeError("Failed to detect any objects")
231
232 measureSources.run(catalog, exp)
233 if len(catalog) == 1:
234 source = catalog[0]
235 else:
236 dmin = None
237 for i, s in enumerate(catalog):
238 d = numpy.hypot(xc - s.getX(), yc - s.getY())
239 if i == 0 or d < dmin:
240 source, dmin = s, d
241 xc, yc = source.getCentroid()
242
243
244 try:
245 subtractPsf(psf, im, xc, yc)
246 except Exception:
247 continue
248
249 resid = im
250 if variance:
251 resid = resid.getImage()
252 var = im.getVariance()
253 var = var.Factory(var, True)
254 numpy.sqrt(var.getArray(), var.getArray())
255 resid /= var
256
257 im_resid.append(resid)
258
259
260 if fitBasisComponents:
261 im = cand.getMaskedImage()
262
263 im = im.Factory(im, True)
264 im.setXY0(cand.getMaskedImage().getXY0())
265
266 try:
267 noSpatialKernel = psf.getKernel()
268 except Exception:
269 noSpatialKernel = None
270
271 if noSpatialKernel:
273 fit = fitKernelParamsToImage(noSpatialKernel, im, candCenter)
274 params = fit[0]
277
278 outImage = afwImage.ImageD(outputKernel.getDimensions())
279 outputKernel.computeImage(outImage, False)
280
281 im -= outImage.convertF()
282 resid = im
283
284 if margin > 0:
285 bim = im.Factory(w + 2*margin, h + 2*margin)
287 bim *= stdev
288
289 bim.assign(resid, bbox)
290 resid = bim
291
292 if variance:
293 resid = resid.getImage()
294 resid /= var
295
296 im_resid.append(resid)
297
298 im = im_resid.makeMosaic()
299 else:
300 im = cand.getMaskedImage()
301
302 if normalize:
304
305 objId = splitId(cand.getSource().getId(), True)["objId"]
306 if psf:
307 lab = "%d chi^2 %.1f" % (objId, rchi2)
308 ctype = afwDisplay.RED if cand.isBad() else afwDisplay.GREEN
309 else:
310 lab = "%d flux %8.3g" % (objId, cand.getSource().getPsfInstFlux())
311 ctype = afwDisplay.GREEN
312
313 mos.append(im, lab, ctype)
314
315 if False and numpy.isnan(rchi2):
316 display.mtv(cand.getMaskedImage().getImage(), title="showPsfCandidates: candidate")
317 print("amp", cand.getAmplitude())
318
319 im = cand.getMaskedImage()
320 center = (candidateIndex, xc - im.getX0(), yc - im.getY0())
321 candidateIndex += 1
322 if cand.isBad():
323 candidateCentersBad.append(center)
324 else:
325 candidateCenters.append(center)
326
327 if variance:
328 title = "chi(Psf fit)"
329 else:
330 title = "Stars & residuals"
331 mosaicImage = mos.makeMosaic(display=display, title=title)
332
333 with display.Buffering():
334 for centers, color in ((candidateCenters, afwDisplay.GREEN), (candidateCentersBad, afwDisplay.RED)):
335 for cen in centers:
336 bbox = mos.getBBox(cen[0])
337 display.dot("+", cen[1] + bbox.getMinX(), cen[2] + bbox.getMinY(), ctype=color)
338
339 return mosaicImage
340
341
A Threshold is used to pass a threshold value to detection algorithms.
A kernel that is a linear combination of fixed basis kernels.
A class that can be used to generate sequences of random numbers according to a number of different a...
An integer coordinate rectangle.
std::shared_ptr< Exposure< ImagePixelT, MaskPixelT, VariancePixelT > > makeExposure(MaskedImage< ImagePixelT, MaskPixelT, VariancePixelT > &mimage, std::shared_ptr< geom::SkyWcs const > wcs=std::shared_ptr< geom::SkyWcs const >())
A function to return an Exposure of the correct type (cf.
void randomGaussianImage(ImageT *image, Random &rand)
Set image to random numbers with a gaussian N(0, 1) distribution.
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)