LSST Applications g0265f82a02+c6dfa2ddaf,g1162b98a3f+b2075782a9,g2079a07aa2+1b2e822518,g2bbee38e9b+c6dfa2ddaf,g337abbeb29+c6dfa2ddaf,g3ddfee87b4+a60788ef87,g50ff169b8f+2eb0e556e8,g52b1c1532d+90ebb246c7,g555ede804d+a60788ef87,g591dd9f2cf+ba8caea58f,g5ec818987f+864ee9cddb,g858d7b2824+9ee1ab4172,g876c692160+a40945ebb7,g8a8a8dda67+90ebb246c7,g8cdfe0ae6a+4fd9e222a8,g99cad8db69+5e309b7bc6,g9ddcbc5298+a1346535a5,ga1e77700b3+df8f93165b,ga8c6da7877+aa12a14d27,gae46bcf261+c6dfa2ddaf,gb0e22166c9+8634eb87fb,gb3f2274832+d0da15e3be,gba4ed39666+1ac82b564f,gbb8dafda3b+5dfd9c994b,gbeb006f7da+97157f9740,gc28159a63d+c6dfa2ddaf,gc86a011abf+9ee1ab4172,gcf0d15dbbd+a60788ef87,gdaeeff99f8+1cafcb7cd4,gdc0c513512+9ee1ab4172,ge79ae78c31+c6dfa2ddaf,geb67518f79+ba1859f325,geb961e4c1e+f9439d1e6f,gee10cc3b42+90ebb246c7,gf1cff7945b+9ee1ab4172,w.2024.12
LSST Data Management Base Package
Loading...
Searching...
No Matches
Functions | Variables
lsst.meas.algorithms.utils Namespace Reference

Functions

 splitId (oid, asDict=True)
 
 showSourceSet (sSet, xy0=(0, 0), display=None, ctype=afwDisplay.GREEN, symb="+", size=2)
 
 showPsfSpatialCells (exposure, psfCellSet, nMaxPerCell=-1, showChi2=False, showMoments=False, symb=None, ctype=None, ctypeUnused=None, ctypeBad=None, size=2, display=None)
 
 showPsfCandidates (exposure, psfCellSet, psf=None, display=None, normalize=True, showBadCandidates=True, fitBasisComponents=False, variance=None, chi=None)
 
 makeSubplots (fig, nx=2, ny=2, Nx=1, Ny=1, plottingArea=(0.1, 0.1, 0.85, 0.80), pxgutter=0.05, pygutter=0.05, xgutter=0.04, ygutter=0.04, headroom=0.0, panelBorderWeight=0, panelColor='black')
 
 plotPsfSpatialModel (exposure, psf, psfCellSet, showBadCandidates=True, numSample=128, matchKernelAmplitudes=False, keepPlots=True)
 
 showPsf (psf, eigenValues=None, XY=None, normalize=True, display=None)
 
 showPsfMosaic (exposure, psf=None, nx=7, ny=None, showCenter=True, showEllipticity=False, showFwhm=False, stampSize=0, display=None, title=None)
 
 showPsfResiduals (exposure, sourceSet, magType="psf", scale=10, display=None)
 
 saveSpatialCellSet (psfCellSet, fileName="foo.fits", display=None)
 

Variables

bool keptPlots = False
 
 _LOG = logging.getLogger(__name__)
 

Detailed Description

Support utilities for Measuring sources

Function Documentation

◆ makeSubplots()

lsst.meas.algorithms.utils.makeSubplots ( fig,
nx = 2,
ny = 2,
Nx = 1,
Ny = 1,
plottingArea = (0.1, 0.1, 0.85, 0.80),
pxgutter = 0.05,
pygutter = 0.05,
xgutter = 0.04,
ygutter = 0.04,
headroom = 0.0,
panelBorderWeight = 0,
panelColor = 'black' )
Return a generator of a set of subplots, a set of Nx*Ny panels of nx*ny plots.  Each panel is fully
filled by row (starting in the bottom left) before the next panel is started.  If panelBorderWidth is
greater than zero a border is drawn around each panel, adjusted to enclose the axis labels.

E.g.
subplots = makeSubplots(fig, 2, 2, Nx=1, Ny=1, panelColor='k')
ax = subplots.next(); ax.text(0.3, 0.5, '[0, 0] (0,0)')
ax = subplots.next(); ax.text(0.3, 0.5, '[0, 0] (1,0)')
ax = subplots.next(); ax.text(0.3, 0.5, '[0, 0] (0,1)')
ax = subplots.next(); ax.text(0.3, 0.5, '[0, 0] (1,1)')
fig.show()

Parameters
----------
fig : `matplotlib.pyplot.figure`
    The matplotlib figure to draw
nx : `int`
    The number of plots in each row of each panel
ny : `int`
    The number of plots in each column of each panel
Nx : `int`
    The number of panels in each row of the figure
Ny : `int`
    The number of panels in each column of the figure
plottingArea : `tuple`
    (x0, y0, x1, y1) for the part of the figure containing all the panels
pxgutter : `float`
    Spacing between columns of panels in units of (x1 - x0)
pygutter : `float`
    Spacing between rows of panels in units of (y1 - y0)
xgutter : `float`
    Spacing between columns of plots within a panel in units of (x1 - x0)
ygutter : `float`
    Spacing between rows of plots within a panel in units of (y1 - y0)
headroom : `float`
    Extra spacing above each plot for e.g. a title
panelBorderWeight : `int`
    Width of border drawn around panels
panelColor : `str`
    Colour of border around panels

Definition at line 342 of file utils.py.

344 headroom=0.0, panelBorderWeight=0, panelColor='black'):
345 """Return a generator of a set of subplots, a set of Nx*Ny panels of nx*ny plots. Each panel is fully
346 filled by row (starting in the bottom left) before the next panel is started. If panelBorderWidth is
347 greater than zero a border is drawn around each panel, adjusted to enclose the axis labels.
348
349 E.g.
350 subplots = makeSubplots(fig, 2, 2, Nx=1, Ny=1, panelColor='k')
351 ax = subplots.next(); ax.text(0.3, 0.5, '[0, 0] (0,0)')
352 ax = subplots.next(); ax.text(0.3, 0.5, '[0, 0] (1,0)')
353 ax = subplots.next(); ax.text(0.3, 0.5, '[0, 0] (0,1)')
354 ax = subplots.next(); ax.text(0.3, 0.5, '[0, 0] (1,1)')
355 fig.show()
356
357 Parameters
358 ----------
359 fig : `matplotlib.pyplot.figure`
360 The matplotlib figure to draw
361 nx : `int`
362 The number of plots in each row of each panel
363 ny : `int`
364 The number of plots in each column of each panel
365 Nx : `int`
366 The number of panels in each row of the figure
367 Ny : `int`
368 The number of panels in each column of the figure
369 plottingArea : `tuple`
370 (x0, y0, x1, y1) for the part of the figure containing all the panels
371 pxgutter : `float`
372 Spacing between columns of panels in units of (x1 - x0)
373 pygutter : `float`
374 Spacing between rows of panels in units of (y1 - y0)
375 xgutter : `float`
376 Spacing between columns of plots within a panel in units of (x1 - x0)
377 ygutter : `float`
378 Spacing between rows of plots within a panel in units of (y1 - y0)
379 headroom : `float`
380 Extra spacing above each plot for e.g. a title
381 panelBorderWeight : `int`
382 Width of border drawn around panels
383 panelColor : `str`
384 Colour of border around panels
385 """
386
387 log = _LOG.getChild("makeSubplots")
388 try:
389 import matplotlib.pyplot as plt
390 except ImportError as e:
391 log.warning("Unable to import matplotlib: %s", e)
392 return
393
394 # Make show() call canvas.draw() too so that we know how large the axis labels are. Sigh
395 try:
396 fig.__show
397 except AttributeError:
398 fig.__show = fig.show
399
400 def myShow(fig):
401 fig.__show()
402 fig.canvas.draw()
403
404 import types
405 fig.show = types.MethodType(myShow, fig)
406 #
407 # We can't get the axis sizes until after draw()'s been called, so use a callback Sigh^2
408 #
409 axes = {} # all axes in all the panels we're drawing: axes[panel][0] etc.
410 #
411
412 def on_draw(event):
413 """
414 Callback to draw the panel borders when the plots are drawn to the canvas
415 """
416 if panelBorderWeight <= 0:
417 return False
418
419 for p in axes.keys():
420 bboxes = []
421 for ax in axes[p]:
422 bboxes.append(ax.bbox.union([label.get_window_extent() for label in
423 ax.get_xticklabels() + ax.get_yticklabels()]))
424
425 ax = axes[p][0]
426
427 # this is the bbox that bounds all the bboxes, again in relative
428 # figure coords
429
430 bbox = ax.bbox.union(bboxes)
431
432 xy0, xy1 = ax.transData.inverted().transform(bbox)
433 x0, y0 = xy0
434 x1, y1 = xy1
435 w, h = x1 - x0, y1 - y0
436 # allow a little space around BBox
437 x0 -= 0.02*w
438 w += 0.04*w
439 y0 -= 0.02*h
440 h += 0.04*h
441 h += h*headroom
442 # draw BBox
443 ax.patches = [] # remove old ones
444 rec = ax.add_patch(plt.Rectangle((x0, y0), w, h, fill=False,
445 lw=panelBorderWeight, edgecolor=panelColor))
446 rec.set_clip_on(False)
447
448 return False
449
450 fig.canvas.mpl_connect('draw_event', on_draw)
451 #
452 # Choose the plotting areas for each subplot
453 #
454 x0, y0 = plottingArea[0:2]
455 W, H = plottingArea[2:4]
456 w = (W - (Nx - 1)*pxgutter - (nx*Nx - 1)*xgutter)/float(nx*Nx)
457 h = (H - (Ny - 1)*pygutter - (ny*Ny - 1)*ygutter)/float(ny*Ny)
458 #
459 # OK! Time to create the subplots
460 #
461 for panel in range(Nx*Ny):
462 axes[panel] = []
463 px = panel%Nx
464 py = panel//Nx
465 for window in range(nx*ny):
466 x = nx*px + window%nx
467 y = ny*py + window//nx
468 ax = fig.add_axes((x0 + xgutter + pxgutter + x*w + (px - 1)*pxgutter + (x - 1)*xgutter,
469 y0 + ygutter + pygutter + y*h + (py - 1)*pygutter + (y - 1)*ygutter,
470 w, h), frame_on=True, facecolor='w')
471 axes[panel].append(ax)
472 yield ax
473
474

◆ plotPsfSpatialModel()

lsst.meas.algorithms.utils.plotPsfSpatialModel ( exposure,
psf,
psfCellSet,
showBadCandidates = True,
numSample = 128,
matchKernelAmplitudes = False,
keepPlots = True )
Plot the PSF spatial model.

Definition at line 475 of file utils.py.

476 matchKernelAmplitudes=False, keepPlots=True):
477 """Plot the PSF spatial model."""
478
479 log = _LOG.getChild("plotPsfSpatialModel")
480 try:
481 import matplotlib.pyplot as plt
482 import matplotlib as mpl
483 except ImportError as e:
484 log.warning("Unable to import matplotlib: %s", e)
485 return
486
487 noSpatialKernel = psf.getKernel()
488 candPos = list()
489 candFits = list()
490 badPos = list()
491 badFits = list()
492 candAmps = list()
493 badAmps = list()
494 for cell in psfCellSet.getCellList():
495 for cand in cell.begin(False):
496 if not showBadCandidates and cand.isBad():
497 continue
498 candCenter = lsst.geom.PointD(cand.getXCenter(), cand.getYCenter())
499 try:
500 im = cand.getMaskedImage()
501 except Exception:
502 continue
503
504 fit = fitKernelParamsToImage(noSpatialKernel, im, candCenter)
505 params = fit[0]
506 kernels = fit[1]
507 amp = 0.0
508 for p, k in zip(params, kernels):
509 amp += p * k.getSum()
510
511 targetFits = badFits if cand.isBad() else candFits
512 targetPos = badPos if cand.isBad() else candPos
513 targetAmps = badAmps if cand.isBad() else candAmps
514
515 targetFits.append([x / amp for x in params])
516 targetPos.append(candCenter)
517 targetAmps.append(amp)
518
519 xGood = numpy.array([pos.getX() for pos in candPos]) - exposure.getX0()
520 yGood = numpy.array([pos.getY() for pos in candPos]) - exposure.getY0()
521 zGood = numpy.array(candFits)
522
523 xBad = numpy.array([pos.getX() for pos in badPos]) - exposure.getX0()
524 yBad = numpy.array([pos.getY() for pos in badPos]) - exposure.getY0()
525 zBad = numpy.array(badFits)
526 numBad = len(badPos)
527
528 xRange = numpy.linspace(0, exposure.getWidth(), num=numSample)
529 yRange = numpy.linspace(0, exposure.getHeight(), num=numSample)
530
531 kernel = psf.getKernel()
532 nKernelComponents = kernel.getNKernelParameters()
533 #
534 # Figure out how many panels we'll need
535 #
536 nPanelX = int(math.sqrt(nKernelComponents))
537 nPanelY = nKernelComponents//nPanelX
538 while nPanelY*nPanelX < nKernelComponents:
539 nPanelX += 1
540
541 fig = plt.figure(1)
542 fig.clf()
543 try:
544 fig.canvas._tkcanvas._root().lift() # == Tk's raise, but raise is a python reserved word
545 except Exception: # protect against API changes
546 pass
547 #
548 # Generator for axes arranged in panels
549 #
550 mpl.rcParams["figure.titlesize"] = "x-small"
551 subplots = makeSubplots(fig, 2, 2, Nx=nPanelX, Ny=nPanelY, xgutter=0.06, ygutter=0.06, pygutter=0.04)
552
553 for k in range(nKernelComponents):
554 func = kernel.getSpatialFunction(k)
555 dfGood = zGood[:, k] - numpy.array([func(pos.getX(), pos.getY()) for pos in candPos])
556 yMin = dfGood.min()
557 yMax = dfGood.max()
558 if numBad > 0:
559 dfBad = zBad[:, k] - numpy.array([func(pos.getX(), pos.getY()) for pos in badPos])
560 yMin = min([yMin, dfBad.min()])
561 yMax = max([yMax, dfBad.max()])
562 yMin -= 0.05 * (yMax - yMin)
563 yMax += 0.05 * (yMax - yMin)
564
565 yMin = -0.01
566 yMax = 0.01
567
568 fRange = numpy.ndarray((len(xRange), len(yRange)))
569 for j, yVal in enumerate(yRange):
570 for i, xVal in enumerate(xRange):
571 fRange[j][i] = func(xVal, yVal)
572
573 ax = next(subplots)
574
575 ax.set_autoscale_on(False)
576 ax.set_xbound(lower=0, upper=exposure.getHeight())
577 ax.set_ybound(lower=yMin, upper=yMax)
578 ax.plot(yGood, dfGood, 'b+')
579 if numBad > 0:
580 ax.plot(yBad, dfBad, 'r+')
581 ax.axhline(0.0)
582 ax.set_title('Residuals(y)')
583
584 ax = next(subplots)
585
586 if matchKernelAmplitudes and k == 0:
587 vmin = 0.0
588 vmax = 1.1
589 else:
590 vmin = fRange.min()
591 vmax = fRange.max()
592
593 norm = mpl.colors.Normalize(vmin=vmin, vmax=vmax)
594 im = ax.imshow(fRange, aspect='auto', origin="lower", norm=norm,
595 extent=[0, exposure.getWidth()-1, 0, exposure.getHeight()-1])
596 ax.set_title('Spatial poly')
597 plt.colorbar(im, orientation='horizontal', ticks=[vmin, vmax])
598
599 ax = next(subplots)
600 ax.set_autoscale_on(False)
601 ax.set_xbound(lower=0, upper=exposure.getWidth())
602 ax.set_ybound(lower=yMin, upper=yMax)
603 ax.plot(xGood, dfGood, 'b+')
604 if numBad > 0:
605 ax.plot(xBad, dfBad, 'r+')
606 ax.axhline(0.0)
607 ax.set_title('K%d Residuals(x)' % k)
608
609 ax = next(subplots)
610
611 photoCalib = exposure.getPhotoCalib()
612 # If there is no calibration factor, use 1.0.
613 if photoCalib.getCalibrationMean() <= 0:
614 photoCalib = afwImage.PhotoCalib(1.0)
615
616 ampMag = [photoCalib.instFluxToMagnitude(candAmp) for candAmp in candAmps]
617 ax.plot(ampMag, zGood[:, k], 'b+')
618 if numBad > 0:
619 badAmpMag = [photoCalib.instFluxToMagnitude(badAmp) for badAmp in badAmps]
620 ax.plot(badAmpMag, zBad[:, k], 'r+')
621
622 ax.set_title('Flux variation')
623
624 fig.show()
625
626 global keptPlots
627 if keepPlots and not keptPlots:
628 # Keep plots open when done
629 def show():
630 print("%s: Please close plots when done." % __name__)
631 try:
632 plt.show()
633 except Exception:
634 pass
635 print("Plots closed, exiting...")
636 import atexit
637 atexit.register(show)
638 keptPlots = True
639
640
int min
int max
The photometric calibration of an exposure.
Definition PhotoCalib.h:114

◆ saveSpatialCellSet()

lsst.meas.algorithms.utils.saveSpatialCellSet ( psfCellSet,
fileName = "foo.fits",
display = None )
Write the contents of a SpatialCellSet to a many-MEF fits file

Definition at line 846 of file utils.py.

846def saveSpatialCellSet(psfCellSet, fileName="foo.fits", display=None):
847 """Write the contents of a SpatialCellSet to a many-MEF fits file"""
848
849 mode = "w"
850 for cell in psfCellSet.getCellList():
851 for cand in cell.begin(False): # include bad candidates
852 dx = afwImage.positionToIndex(cand.getXCenter(), True)[1]
853 dy = afwImage.positionToIndex(cand.getYCenter(), True)[1]
854 im = afwMath.offsetImage(cand.getMaskedImage(), -dx, -dy, "lanczos5")
855
857 md.set("CELL", cell.getLabel())
858 md.set("ID", cand.getId())
859 md.set("XCENTER", cand.getXCenter())
860 md.set("YCENTER", cand.getYCenter())
861 md.set("BAD", cand.isBad())
862 md.set("AMPL", cand.getAmplitude())
863 md.set("FLUX", cand.getSource().getPsfInstFlux())
864 md.set("CHI2", cand.getSource().getChi2())
865
866 im.writeFits(fileName, md, mode)
867 mode = "a"
868
869 if display:
870 display.mtv(im, title="saveSpatialCellSet: image")
Class for storing generic metadata.
Definition PropertySet.h:66
int positionToIndex(double pos)
Convert image position to nearest integer index.
Definition ImageUtils.h:69
std::shared_ptr< ImageT > offsetImage(ImageT const &image, float dx, float dy, std::string const &algorithmName="lanczos5", unsigned int buffer=0)
Return an image offset by (dx, dy) using the specified algorithm.

◆ showPsf()

lsst.meas.algorithms.utils.showPsf ( psf,
eigenValues = None,
XY = None,
normalize = True,
display = None )
Display a PSF's eigen images

If normalize is True, set the largest absolute value of each eigenimage to 1.0 (n.b. sum == 0.0 for i > 0)

Definition at line 641 of file utils.py.

641def showPsf(psf, eigenValues=None, XY=None, normalize=True, display=None):
642 """Display a PSF's eigen images
643
644 If normalize is True, set the largest absolute value of each eigenimage to 1.0 (n.b. sum == 0.0 for i > 0)
645 """
646
647 if eigenValues:
648 coeffs = eigenValues
649 elif XY is not None:
650 coeffs = psf.getLocalKernel(lsst.geom.PointD(XY[0], XY[1])).getKernelParameters()
651 else:
652 coeffs = None
653
654 mos = displayUtils.Mosaic(gutter=2, background=-0.1)
655 for i, k in enumerate(psf.getKernel().getKernelList()):
656 im = afwImage.ImageD(k.getDimensions())
657 k.computeImage(im, False)
658 if normalize:
659 im /= numpy.max(numpy.abs(im.getArray()))
660
661 if coeffs:
662 mos.append(im, "%g" % (coeffs[i]/coeffs[0]))
663 else:
664 mos.append(im)
665
666 if not display:
667 display = afwDisplay.Display()
668 mos.makeMosaic(display=display, title="Kernel Basis Functions")
669
670 return mos
671
672

◆ showPsfCandidates()

lsst.meas.algorithms.utils.showPsfCandidates ( exposure,
psfCellSet,
psf = None,
display = None,
normalize = True,
showBadCandidates = True,
fitBasisComponents = False,
variance = None,
chi = None )
Display the PSF candidates.

If psf is provided include PSF model and residuals;  if normalize is true normalize the PSFs
(and residuals)

If chi is True, generate a plot of residuals/sqrt(variance), i.e. chi

If fitBasisComponents is true, also find the best linear combination of the PSF's components
(if they exist)

Definition at line 137 of file utils.py.

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: # old name for chi
154 chi = variance
155 #
156 # Show us the ccandidates
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): # include bad candidates
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() # copy of this object's image
178 xc, yc = cand.getXCenter(), cand.getYCenter()
179
180 margin = 0 if True else 5
181 w, h = im.getDimensions()
182 bbox = lsst.geom.BoxI(lsst.geom.PointI(margin, margin), im.getDimensions())
183
184 if margin > 0:
185 bim = im.Factory(w + 2*margin, h + 2*margin)
186
187 stdev = numpy.sqrt(afwMath.makeStatistics(im.getVariance(), afwMath.MEAN).getValue())
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: # tweak up centroids
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)
212 catalog = afwTable.SourceCatalog(schema)
213
214 extra = 10 # enough margin to run the sdss centroider
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
221 exp = afwImage.makeExposure(mi)
222 exp.setPsf(psf)
223
224 footprintSet = afwDet.FootprintSet(mi,
225 afwDet.Threshold(0.5*numpy.max(psfIm.getArray())),
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: # more than one source; find the once closest to (xc, yc)
236 dmin = None # an invalid value to catch logic errors
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 # residuals using spatial model
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()) # inplace sqrt
255 resid /= var
256
257 im_resid.append(resid)
258
259 # Fit the PSF components directly to the data (i.e. ignoring the spatial model)
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:
272 candCenter = lsst.geom.PointD(cand.getXCenter(), cand.getYCenter())
273 fit = fitKernelParamsToImage(noSpatialKernel, im, candCenter)
274 params = fit[0]
275 kernels = afwMath.KernelList(fit[1])
276 outputKernel = afwMath.LinearCombinationKernel(kernels, params)
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:
303 im /= afwMath.makeStatistics(im, afwMath.MAX).getValue()
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 set of Footprints, associated with a MaskedImage.
A Threshold is used to pass a threshold value to detection algorithms.
Definition Threshold.h:43
A kernel that is a linear combination of fixed basis kernels.
Definition Kernel.h:704
A class that can be used to generate sequences of random numbers according to a number of different a...
Definition Random.h:57
An integer coordinate rectangle.
Definition Box.h:55
daf::base::PropertySet * set
Definition fits.cc:931
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.
Definition Exposure.h:484
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)
Definition Statistics.h:361

◆ showPsfMosaic()

lsst.meas.algorithms.utils.showPsfMosaic ( exposure,
psf = None,
nx = 7,
ny = None,
showCenter = True,
showEllipticity = False,
showFwhm = False,
stampSize = 0,
display = None,
title = None )
Show a mosaic of Psf images.  exposure may be an Exposure (optionally with PSF),
or a tuple (width, height)

If stampSize is > 0, the psf images will be trimmed to stampSize*stampSize

Definition at line 673 of file utils.py.

674 showFwhm=False, stampSize=0, display=None, title=None):
675 """Show a mosaic of Psf images. exposure may be an Exposure (optionally with PSF),
676 or a tuple (width, height)
677
678 If stampSize is > 0, the psf images will be trimmed to stampSize*stampSize
679 """
680
681 scale = 1.0
682 if showFwhm:
683 showEllipticity = True
684 scale = 2*math.log(2) # convert sigma^2 to HWHM^2 for a Gaussian
685
686 mos = displayUtils.Mosaic()
687
688 try: # maybe it's a real Exposure
689 width, height = exposure.getWidth(), exposure.getHeight()
690 x0, y0 = exposure.getXY0()
691 if not psf:
692 psf = exposure.getPsf()
693 except AttributeError:
694 try: # OK, maybe a list [width, height]
695 width, height = exposure[0], exposure[1]
696 x0, y0 = 0, 0
697 except TypeError: # I guess not
698 raise RuntimeError("Unable to extract width/height from object of type %s" % type(exposure))
699
700 if not ny:
701 ny = int(nx*float(height)/width + 0.5)
702 if not ny:
703 ny = 1
704
705 centroidName = "SdssCentroid"
706 shapeName = "base_SdssShape"
707
708 schema = afwTable.SourceTable.makeMinimalSchema()
709 schema.getAliasMap().set("slot_Centroid", centroidName)
710 schema.getAliasMap().set("slot_Centroid_flag", centroidName+"_flag")
711
712 control = measBase.SdssCentroidControl()
713 centroider = measBase.SdssCentroidAlgorithm(control, centroidName, schema)
714
715 sdssShape = measBase.SdssShapeControl()
716 shaper = measBase.SdssShapeAlgorithm(sdssShape, shapeName, schema)
717 table = afwTable.SourceTable.make(schema)
718
719 table.defineCentroid(centroidName)
720 table.defineShape(shapeName)
721
722 bbox = None
723 if stampSize > 0:
724 w, h = psf.computeImage(lsst.geom.PointD(0, 0)).getDimensions()
725 if stampSize <= w and stampSize <= h:
726 bbox = lsst.geom.BoxI(lsst.geom.PointI((w - stampSize)//2, (h - stampSize)//2),
727 lsst.geom.ExtentI(stampSize, stampSize))
728
729 centers = []
730 shapes = []
731 for iy in range(ny):
732 for ix in range(nx):
733 x = int(ix*(width-1)/(nx-1)) + x0
734 y = int(iy*(height-1)/(ny-1)) + y0
735
736 im = psf.computeImage(lsst.geom.PointD(x, y)).convertF()
737 imPeak = psf.computePeak(lsst.geom.PointD(x, y))
738 im /= imPeak
739 if bbox:
740 im = im.Factory(im, bbox)
741 lab = "PSF(%d,%d)" % (x, y) if False else ""
742 mos.append(im, lab)
743
745 exp.setPsf(psf)
746 w, h = im.getWidth(), im.getHeight()
747 centerX = im.getX0() + w//2
748 centerY = im.getY0() + h//2
749 src = table.makeRecord()
750 spans = afwGeom.SpanSet(exp.getBBox())
751 foot = afwDet.Footprint(spans)
752 foot.addPeak(centerX, centerY, 1)
753 src.setFootprint(foot)
754
755 try:
756 centroider.measure(src, exp)
757 centers.append((src.getX() - im.getX0(), src.getY() - im.getY0()))
758
759 shaper.measure(src, exp)
760 shapes.append((src.getIxx(), src.getIxy(), src.getIyy()))
761 except Exception:
762 pass
763
764 if not display:
765 display = afwDisplay.Display()
766 mos.makeMosaic(display=display, title=title if title else "Model Psf", mode=nx)
767
768 if centers and display:
769 with display.Buffering():
770 for i, (cen, shape) in enumerate(zip(centers, shapes)):
771 bbox = mos.getBBox(i)
772 xc, yc = cen[0] + bbox.getMinX(), cen[1] + bbox.getMinY()
773 if showCenter:
774 display.dot("+", xc, yc, ctype=afwDisplay.BLUE)
775
776 if showEllipticity:
777 ixx, ixy, iyy = shape
778 ixx *= scale
779 ixy *= scale
780 iyy *= scale
781 display.dot("@:%g,%g,%g" % (ixx, ixy, iyy), xc, yc, ctype=afwDisplay.RED)
782
783 return mos
784
785
Class to describe the properties of a detected object from an image.
Definition Footprint.h:63
A compact representation of a collection of pixels.
Definition SpanSet.h:78
MaskedImage< ImagePixelT, MaskPixelT, VariancePixelT > * makeMaskedImage(typename std::shared_ptr< Image< ImagePixelT > > image, typename std::shared_ptr< Mask< MaskPixelT > > mask=Mask< MaskPixelT >(), typename std::shared_ptr< Image< VariancePixelT > > variance=Image< VariancePixelT >())
A function to return a MaskedImage of the correct type (cf.

◆ showPsfResiduals()

lsst.meas.algorithms.utils.showPsfResiduals ( exposure,
sourceSet,
magType = "psf",
scale = 10,
display = None )

Definition at line 786 of file utils.py.

786def showPsfResiduals(exposure, sourceSet, magType="psf", scale=10, display=None):
787 mimIn = exposure.getMaskedImage()
788 mimIn = mimIn.Factory(mimIn, True) # make a copy to subtract from
789
790 psf = exposure.getPsf()
791 psfWidth, psfHeight = psf.getLocalKernel(psf.getAveragePosition()).getDimensions()
792 #
793 # Make the image that we'll paste our residuals into. N.b. they can overlap the edges
794 #
795 w, h = int(mimIn.getWidth()/scale), int(mimIn.getHeight()/scale)
796
797 im = mimIn.Factory(w + psfWidth, h + psfHeight)
798
799 cenPos = []
800 for s in sourceSet:
801 x, y = s.getX(), s.getY()
802
803 sx, sy = int(x/scale + 0.5), int(y/scale + 0.5)
804
805 smim = im.Factory(im, lsst.geom.BoxI(lsst.geom.PointI(sx, sy),
806 lsst.geom.ExtentI(psfWidth, psfHeight)))
807 sim = smim.getImage()
808
809 try:
810 if magType == "ap":
811 flux = s.getApInstFlux()
812 elif magType == "model":
813 flux = s.getModelInstFlux()
814 elif magType == "psf":
815 flux = s.getPsfInstFlux()
816 else:
817 raise RuntimeError("Unknown flux type %s" % magType)
818
819 subtractPsf(psf, mimIn, x, y, flux)
820 except Exception as e:
821 print(e)
822
823 try:
824 expIm = mimIn.getImage().Factory(mimIn.getImage(),
825 lsst.geom.BoxI(lsst.geom.PointI(int(x) - psfWidth//2,
826 int(y) - psfHeight//2),
827 lsst.geom.ExtentI(psfWidth, psfHeight)),
828 )
829 except pexExcept.Exception:
830 continue
831
832 cenPos.append([x - expIm.getX0() + sx, y - expIm.getY0() + sy])
833
834 sim += expIm
835
836 if display:
837 display = afwDisplay.Display()
838 display.mtv(im, title="showPsfResiduals: image")
839 with display.Buffering():
840 for x, y in cenPos:
841 display.dot("+", x, y)
842
843 return im
844
845
Provides consistent interface for LSST exceptions.
Definition Exception.h:107

◆ showPsfSpatialCells()

lsst.meas.algorithms.utils.showPsfSpatialCells ( exposure,
psfCellSet,
nMaxPerCell = -1,
showChi2 = False,
showMoments = False,
symb = None,
ctype = None,
ctypeUnused = None,
ctypeBad = None,
size = 2,
display = None )
Show the SpatialCells.

If symb is something that afwDisplay.Display.dot() understands (e.g. "o"),
the top nMaxPerCell candidates will be indicated with that symbol, using
ctype and size.

Definition at line 81 of file utils.py.

82 symb=None, ctype=None, ctypeUnused=None, ctypeBad=None, size=2, display=None):
83 """Show the SpatialCells.
84
85 If symb is something that afwDisplay.Display.dot() understands (e.g. "o"),
86 the top nMaxPerCell candidates will be indicated with that symbol, using
87 ctype and size.
88 """
89
90 if not display:
91 display = afwDisplay.Display()
92 with display.Buffering():
93 origin = [-exposure.getMaskedImage().getX0(), -exposure.getMaskedImage().getY0()]
94 for cell in psfCellSet.getCellList():
95 displayUtils.drawBBox(cell.getBBox(), origin=origin, display=display)
96
97 if nMaxPerCell < 0:
98 nMaxPerCell = 0
99
100 i = 0
101 goodies = ctypeBad is None
102 for cand in cell.begin(goodies):
103 if nMaxPerCell > 0:
104 i += 1
105
106 xc, yc = cand.getXCenter() + origin[0], cand.getYCenter() + origin[1]
107
108 if i > nMaxPerCell:
109 if not ctypeUnused:
110 continue
111
112 color = ctypeBad if cand.isBad() else ctype
113
114 if symb:
115 if i > nMaxPerCell:
116 ct = ctypeUnused
117 else:
118 ct = ctype
119
120 display.dot(symb, xc, yc, ctype=ct, size=size)
121
122 source = cand.getSource()
123
124 if showChi2:
125 rchi2 = cand.getChi2()
126 if rchi2 > 1e100:
127 rchi2 = numpy.nan
128 display.dot("%d %.1f" % (splitId(source.getId(), True)["objId"], rchi2),
129 xc - size, yc - size - 4, ctype=color, size=2)
130
131 if showMoments:
132 display.dot("%.2f %.2f %.2f" % (source.getIxx(), source.getIxy(), source.getIyy()),
133 xc-size, yc + size + 4, ctype=color, size=size)
134 return display
135
136

◆ showSourceSet()

lsst.meas.algorithms.utils.showSourceSet ( sSet,
xy0 = (0, 0),
display = None,
ctype = afwDisplay.GREEN,
symb = "+",
size = 2 )
Draw the (XAstrom, YAstrom) positions of a set of Sources.  Image has the given XY0

Definition at line 62 of file utils.py.

62def showSourceSet(sSet, xy0=(0, 0), display=None, ctype=afwDisplay.GREEN, symb="+", size=2):
63 """Draw the (XAstrom, YAstrom) positions of a set of Sources. Image has the given XY0"""
64
65 if not display:
66 display = afwDisplay.Display()
67 with display.Buffering():
68 for s in sSet:
69 xc, yc = s.getXAstrom() - xy0[0], s.getYAstrom() - xy0[1]
70
71 if symb == "id":
72 display.dot(str(splitId(s.getId(), True)["objId"]), xc, yc, ctype=ctype, size=size)
73 else:
74 display.dot(symb, xc, yc, ctype=ctype, size=size)
75
76#
77# PSF display utilities
78#
79
80

◆ splitId()

lsst.meas.algorithms.utils.splitId ( oid,
asDict = True )

Definition at line 52 of file utils.py.

52def splitId(oid, asDict=True):
53
54 objId = int((oid & 0xffff) - 1) # Should be the value set by apps code
55
56 if asDict:
57 return dict(objId=objId)
58 else:
59 return [objId]
60
61

Variable Documentation

◆ _LOG

lsst.meas.algorithms.utils._LOG = logging.getLogger(__name__)
protected

Definition at line 49 of file utils.py.

◆ keptPlots

bool lsst.meas.algorithms.utils.keptPlots = False

Definition at line 45 of file utils.py.