LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
Functions | Variables
lsst.meas.algorithms.utils Namespace Reference

Functions

def splitId (oid, asDict=True)
 
def showSourceSet (sSet, xy0=(0, 0), display=None, ctype=afwDisplay.GREEN, symb="+", size=2)
 
def showPsfSpatialCells (exposure, psfCellSet, nMaxPerCell=-1, showChi2=False, showMoments=False, symb=None, ctype=None, ctypeUnused=None, ctypeBad=None, size=2, display=None)
 
def showPsfCandidates (exposure, psfCellSet, psf=None, display=None, normalize=True, showBadCandidates=True, fitBasisComponents=False, variance=None, chi=None)
 
def 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')
 
def plotPsfSpatialModel (exposure, psf, psfCellSet, showBadCandidates=True, numSample=128, matchKernelAmplitudes=False, keepPlots=True)
 
def showPsf (psf, eigenValues=None, XY=None, normalize=True, display=None)
 
def showPsfMosaic (exposure, psf=None, nx=7, ny=None, showCenter=True, showEllipticity=False, showFwhm=False, stampSize=0, display=None, title=None)
 
def showPsfResiduals (exposure, sourceSet, magType="psf", scale=10, display=None)
 
def saveSpatialCellSet (psfCellSet, fileName="foo.fits", display=None)
 

Variables

bool keptPlots = False
 

Function Documentation

◆ makeSubplots()

def 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 
std::shared_ptr< FrameSet > append(FrameSet const &first, FrameSet const &second)
Construct a FrameSet that performs two transformations in series.
Definition: functional.cc:33

◆ plotPsfSpatialModel()

def 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
daf::base::PropertyList * list
Definition: fits.cc:913
def show(frame=None)
Definition: ds9.py:88
def 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')
Definition: utils.py:344
std::pair< std::vector< double >, afw::math::KernelList > fitKernelParamsToImage(afw::math::LinearCombinationKernel const &kernel, Image const &image, geom::Point2D const &pos)
Fit a LinearCombinationKernel to an Image, allowing the coefficients of the components to vary.

◆ saveSpatialCellSet()

def 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.

846 def 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 
856  md = dafBase.PropertySet()
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.
Definition: offsetImage.cc:41
def saveSpatialCellSet(psfCellSet, fileName="foo.fits", display=None)
Definition: utils.py:846

◆ showPsf()

def 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.

641 def 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 
def showPsf(psf, eigenValues=None, XY=None, normalize=True, display=None)
Definition: utils.py:641

◆ showPsfCandidates()

def 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.
Definition: FootprintSet.h:53
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:912
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:462
void randomGaussianImage(ImageT *image, Random &rand)
Set image to random numbers with a gaussian N(0, 1) distribution.
Definition: RandomImage.cc:130
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:359
def splitId(oid, asDict=True)
Definition: utils.py:52
double subtractPsf(afw::detection::Psf const &psf, ImageT *data, double x, double y, double psfFlux=std::numeric_limits< double >::quiet_NaN())

◆ showPsfMosaic()

def 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 
table::Key< int > type
Definition: Detector.cc:163
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.
Definition: MaskedImage.h:1240

◆ showPsfResiduals()

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

Definition at line 786 of file utils.py.

786 def 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().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
def showPsfResiduals(exposure, sourceSet, magType="psf", scale=10, display=None)
Definition: utils.py:786

◆ showPsfSpatialCells()

def 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()

def 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.

62 def 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 
def showSourceSet(sSet, xy0=(0, 0), display=None, ctype=afwDisplay.GREEN, symb="+", size=2)
Definition: utils.py:62

◆ splitId()

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

Definition at line 52 of file utils.py.

52 def 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

◆ keptPlots

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

Definition at line 45 of file utils.py.