LSSTApplications  10.0+286,10.0+36,10.0+46,10.0-2-g4f67435,10.1+152,10.1+37,11.0,11.0+1,11.0-1-g47edd16,11.0-1-g60db491,11.0-1-g7418c06,11.0-2-g04d2804,11.0-2-g68503cd,11.0-2-g818369d,11.0-2-gb8b8ce7
LSSTDataManagementBasePackage
Functions | Variables
lsst.meas.algorithms.utils Namespace Reference

Functions

def splitId
 
def showSourceSet
 
def showPsfSpatialCells
 
def showPsfCandidates
 
def makeSubplots
 
def plotPsfSpatialModel
 
def showPsf
 
def showPsfMosaic
 
def showPsfResiduals
 
def saveSpatialCellSet
 

Variables

 keptPlots = False
 
 plt = None
 

Function Documentation

def lsst.meas.algorithms.utils.makeSubplots (   fig,
  nx = 2,
  ny = 2,
  Nx = 1,
  Ny = 1,
  plottingArea = (0.1, 0.1,
  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()

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

Definition at line 336 of file utils.py.

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

441  matchKernelAmplitudes=False, keepPlots=True):
442  """Plot the PSF spatial model."""
443 
444  if not plt:
445  print >> sys.stderr, "Unable to import matplotlib"
446  return
447 
448  noSpatialKernel = afwMath.cast_LinearCombinationKernel(psf.getKernel())
449  candPos = list()
450  candFits = list()
451  badPos = list()
452  badFits = list()
453  candAmps = list()
454  badAmps = list()
455  for cell in psfCellSet.getCellList():
456  for cand in cell.begin(False):
457  cand = algorithmsLib.cast_PsfCandidateF(cand)
458  if not showBadCandidates and cand.isBad():
459  continue
460  candCenter = afwGeom.PointD(cand.getXCenter(), cand.getYCenter())
461  try:
462  im = cand.getMaskedImage()
463  except Exception:
464  continue
465 
466  fit = algorithmsLib.fitKernelParamsToImage(noSpatialKernel, im, candCenter)
467  params = fit[0]
468  kernels = fit[1]
469  amp = 0.0
470  for p, k in zip(params, kernels):
471  amp += p * afwMath.cast_FixedKernel(k).getSum()
472 
473  targetFits = badFits if cand.isBad() else candFits
474  targetPos = badPos if cand.isBad() else candPos
475  targetAmps = badAmps if cand.isBad() else candAmps
476 
477  targetFits.append([x / amp for x in params])
478  targetPos.append(candCenter)
479  targetAmps.append(amp)
480 
481  numCandidates = len(candFits)
482  numBasisFuncs = noSpatialKernel.getNBasisKernels()
483 
484  xGood = numpy.array([pos.getX() for pos in candPos]) - exposure.getX0()
485  yGood = numpy.array([pos.getY() for pos in candPos]) - exposure.getY0()
486  zGood = numpy.array(candFits)
487  ampGood = numpy.array(candAmps)
488 
489  xBad = numpy.array([pos.getX() for pos in badPos]) - exposure.getX0()
490  yBad = numpy.array([pos.getY() for pos in badPos]) - exposure.getY0()
491  zBad = numpy.array(badFits)
492  ampBad = numpy.array(badAmps)
493  numBad = len(badPos)
494 
495  xRange = numpy.linspace(0, exposure.getWidth(), num=numSample)
496  yRange = numpy.linspace(0, exposure.getHeight(), num=numSample)
497 
498  kernel = psf.getKernel()
499  nKernelComponents = kernel.getNKernelParameters()
500  #
501  # Figure out how many panels we'll need
502  #
503  nPanelX = int(math.sqrt(nKernelComponents))
504  nPanelY = nKernelComponents//nPanelX
505  while nPanelY*nPanelX < nKernelComponents:
506  nPanelX += 1
507 
508  fig = plt.figure(1)
509  fig.clf()
510  try:
511  fig.canvas._tkcanvas._root().lift() # == Tk's raise, but raise is a python reserved word
512  except: # protect against API changes
513  pass
514  #
515  # Generator for axes arranged in panels
516  #
517  subplots = makeSubplots(fig, 2, 2, Nx=nPanelX, Ny=nPanelY, xgutter=0.06, ygutter=0.06, pygutter=0.04)
518 
519  for k in range(nKernelComponents):
520  func = kernel.getSpatialFunction(k)
521  dfGood = zGood[:,k] - numpy.array([func(pos.getX(), pos.getY()) for pos in candPos])
522  yMin = dfGood.min()
523  yMax = dfGood.max()
524  if numBad > 0:
525  dfBad = zBad[:,k] - numpy.array([func(pos.getX(), pos.getY()) for pos in badPos])
526  yMin = min([yMin, dfBad.min()])
527  yMax = max([yMax, dfBad.max()])
528  yMin -= 0.05 * (yMax - yMin)
529  yMax += 0.05 * (yMax - yMin)
530 
531  yMin = -0.01
532  yMax = 0.01
533 
534  fRange = numpy.ndarray((len(xRange), len(yRange)))
535  for j, yVal in enumerate(yRange):
536  for i, xVal in enumerate(xRange):
537  fRange[j][i] = func(xVal, yVal)
538 
539  #-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
540 
541  ax = subplots.next()
542 
543  ax.set_autoscale_on(False)
544  ax.set_xbound(lower=0, upper=exposure.getHeight())
545  ax.set_ybound(lower=yMin, upper=yMax)
546  ax.plot(yGood, dfGood, 'b+')
547  if numBad > 0:
548  ax.plot(yBad, dfBad, 'r+')
549  ax.axhline(0.0)
550  ax.set_title('Residuals(y)')
551 
552  #-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
553 
554  ax = subplots.next()
555 
556  if matchKernelAmplitudes and k == 0:
557  vmin = 0.0
558  vmax = 1.1
559  else:
560  vmin = fRange.min()
561  vmax = fRange.max()
562 
563  norm = matplotlib.colors.Normalize(vmin=vmin, vmax=vmax)
564  im = ax.imshow(fRange, aspect='auto', origin="lower", norm=norm,
565  extent=[0, exposure.getWidth()-1, 0, exposure.getHeight()-1])
566  ax.set_title('Spatial poly')
567  plt.colorbar(im, orientation='horizontal', ticks=[vmin, vmax])
568 
569  #-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
570 
571  ax = subplots.next()
572  ax.set_autoscale_on(False)
573  ax.set_xbound(lower=0, upper=exposure.getWidth())
574  ax.set_ybound(lower=yMin, upper=yMax)
575  ax.plot(xGood, dfGood, 'b+')
576  if numBad > 0:
577  ax.plot(xBad, dfBad, 'r+')
578  ax.axhline(0.0)
579  ax.set_title('K%d Residuals(x)' % k)
580 
581  #-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
582 
583  ax = subplots.next()
584 
585  if False:
586  ax.scatter(xGood, yGood, c=dfGood, marker='o')
587  ax.scatter(xBad, yBad, c=dfBad, marker='x')
588  ax.set_xbound(lower=0, upper=exposure.getWidth())
589  ax.set_ybound(lower=0, upper=exposure.getHeight())
590  ax.set_title('Spatial residuals')
591  plt.colorbar(im, orientation='horizontal')
592  else:
593  calib = exposure.getCalib()
594  if calib.getFluxMag0()[0] <= 0:
595  calib = type(calib)()
596  calib.setFluxMag0(1.0)
597 
598  with CalibNoThrow():
599  ax.plot(calib.getMagnitude(candAmps), zGood[:,k], 'b+')
600  if numBad > 0:
601  ax.plot(calib.getMagnitude(badAmps), zBad[:,k], 'r+')
602 
603  ax.set_title('Flux variation')
604 
605  fig.show()
606 
607  global keptPlots
608  if keepPlots and not keptPlots:
609  # Keep plots open when done
610  def show():
611  print "%s: Please close plots when done." % __name__
612  try:
613  plt.show()
614  except:
615  pass
616  print "Plots closed, exiting..."
617  import atexit
618  atexit.register(show)
619  keptPlots = True
def lsst.meas.algorithms.utils.saveSpatialCellSet (   psfCellSet,
  fileName = "foo.fits",
  frame = None 
)
Write the contents of a SpatialCellSet to a many-MEF fits file

Definition at line 827 of file utils.py.

828 def saveSpatialCellSet(psfCellSet, fileName="foo.fits", frame=None):
829  """Write the contents of a SpatialCellSet to a many-MEF fits file"""
830 
831  mode = "w"
832  for cell in psfCellSet.getCellList():
833  for cand in cell.begin(False): # include bad candidates
834  cand = algorithmsLib.cast_PsfCandidateF(cand)
835 
836  dx = afwImage.positionToIndex(cand.getXCenter(), True)[1]
837  dy = afwImage.positionToIndex(cand.getYCenter(), True)[1]
838  im = afwMath.offsetImage(cand.getMaskedImage(), -dx, -dy, "lanczos5")
839 
840  md = dafBase.PropertySet()
841  md.set("CELL", cell.getLabel())
842  md.set("ID", cand.getId())
843  md.set("XCENTER", cand.getXCenter())
844  md.set("YCENTER", cand.getYCenter())
845  md.set("BAD", cand.isBad())
846  md.set("AMPL", cand.getAmplitude())
847  md.set("FLUX", cand.getSource().getPsfFlux())
848  md.set("CHI2", cand.getSource().getChi2())
849 
850  im.writeFits(fileName, md, mode)
851  mode = "a"
852 
853  if frame is not None:
854  ds9.mtv(im, frame=frame)
std::pair< int, double > positionToIndex(double const pos, bool)
Convert image position to index (nearest integer and fractional parts)
Definition: ImageUtils.h:100
ImageT::Ptr 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:55
Class for storing generic metadata.
Definition: PropertySet.h:82
def lsst.meas.algorithms.utils.showPsf (   psf,
  eigenValues = None,
  XY = None,
  normalize = True,
  frame = 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 620 of file utils.py.

621 def showPsf(psf, eigenValues=None, XY=None, normalize=True, frame=None):
622  """Display a PSF's eigen images
623 
624  If normalize is True, set the largest absolute value of each eigenimage to 1.0 (n.b. sum == 0.0 for i > 0)
625  """
626 
627  if eigenValues:
628  coeffs = eigenValues
629  elif XY is not None:
630  coeffs = psf.getLocalKernel(afwGeom.PointD(XY[0], XY[1])).getKernelParameters()
631  else:
632  coeffs = None
633 
634  mos = displayUtils.Mosaic(gutter=2, background=-0.1)
635  for i, k in enumerate(afwMath.cast_LinearCombinationKernel(psf.getKernel()).getKernelList()):
636  im = afwImage.ImageD(k.getDimensions())
637  k.computeImage(im, False)
638  if normalize:
639  im /= numpy.max(numpy.abs(im.getArray()))
640 
641  if coeffs:
642  mos.append(im, "%g" % (coeffs[i]/coeffs[0]))
643  else:
644  mos.append(im)
645 
646  mos.makeMosaic(frame=frame, title="Kernel Basis Functions")
647 
648  return mos
def lsst.meas.algorithms.utils.showPsfCandidates (   exposure,
  psfCellSet,
  psf = None,
  frame = 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 123 of file utils.py.

124  fitBasisComponents=False, variance=None, chi=None):
125  """Display the PSF candidates.
126 
127  If psf is provided include PSF model and residuals; if normalize is true normalize the PSFs
128  (and residuals)
129 
130  If chi is True, generate a plot of residuals/sqrt(variance), i.e. chi
131 
132  If fitBasisComponents is true, also find the best linear combination of the PSF's components
133  (if they exist)
134  """
135  if chi is None:
136  if variance is not None: # old name for chi
137  chi = variance
138  #
139  # Show us the ccandidates
140  #
141  mos = displayUtils.Mosaic()
142  #
143  candidateCenters = []
144  candidateCentersBad = []
145  candidateIndex = 0
146 
147  for cell in psfCellSet.getCellList():
148  for cand in cell.begin(False): # include bad candidates
149  cand = algorithmsLib.cast_PsfCandidateF(cand)
150 
151  rchi2 = cand.getChi2()
152  if rchi2 > 1e100:
153  rchi2 = numpy.nan
154 
155  if not showBadCandidates and cand.isBad():
156  continue
157 
158  if psf:
159  im_resid = displayUtils.Mosaic(gutter=0, background=-5, mode="x")
160 
161  try:
162  im = cand.getMaskedImage() # copy of this object's image
163  xc, yc = cand.getXCenter(), cand.getYCenter()
164 
165  margin = 0 if True else 5
166  w, h = im.getDimensions()
167  bbox = afwGeom.BoxI(afwGeom.PointI(margin, margin), im.getDimensions())
168 
169  if margin > 0:
170  bim = im.Factory(w + 2*margin, h + 2*margin)
171 
172  stdev = numpy.sqrt(afwMath.makeStatistics(im.getVariance(), afwMath.MEAN).getValue())
174  bim *= stdev
175  var = bim.getVariance(); var.set(stdev**2); del var
176 
177  sbim = im.Factory(bim, bbox)
178  sbim <<= im
179  del sbim
180  im = bim
181  xc += margin; yc += margin
182 
183  im = im.Factory(im, True)
184  im.setXY0(cand.getMaskedImage().getXY0())
185  except:
186  continue
187 
188  if not variance:
189  im_resid.append(im.Factory(im, True))
190 
191  if True: # tweak up centroids
192  mi = im
193  psfIm = mi.getImage()
194  config = measBase.SingleFrameMeasurementTask.ConfigClass()
195  config.slots.centroid = "base_SdssCentroid"
196 
197  schema = afwTable.SourceTable.makeMinimalSchema()
198  measureSources = measBase.SingleFrameMeasurementTask(schema, config=config)
199  catalog = afwTable.SourceCatalog(schema)
200 
201  extra = 10 # enough margin to run the sdss centroider
202  miBig = mi.Factory(im.getWidth() + 2*extra, im.getHeight() + 2*extra)
203  miBig[extra:-extra, extra:-extra] = mi
204  miBig.setXY0(mi.getX0() - extra, mi.getY0() - extra)
205  mi = miBig; del miBig
206 
207  exp = afwImage.makeExposure(mi)
208  exp.setPsf(psf)
209 
210  footprintSet = afwDet.FootprintSet(mi,
211  afwDet.Threshold(0.5*numpy.max(psfIm.getArray())),
212  "DETECTED")
213  footprintSet.makeSources(catalog)
214 
215  if len(catalog) == 0:
216  raise RuntimeError("Failed to detect any objects")
217 
218  measureSources.run(catalog, exp)
219  if len(catalog) == 1:
220  source = catalog[0]
221  else: # more than one source; find the once closest to (xc, yc)
222  for i, s in enumerate(catalog):
223  d = numpy.hypot(xc - s.getX(), yc - s.getY())
224  if i == 0 or d < dmin:
225  source, dmin = s, d
226  xc, yc = source.getCentroid()
227 
228  # residuals using spatial model
229  try:
230  chi2 = algorithmsLib.subtractPsf(psf, im, xc, yc)
231  except:
232  chi2 = numpy.nan
233  continue
234 
235  resid = im
236  if variance:
237  resid = resid.getImage()
238  var = im.getVariance()
239  var = var.Factory(var, True)
240  numpy.sqrt(var.getArray(), var.getArray()) # inplace sqrt
241  resid /= var
242 
243  im_resid.append(resid)
244 
245  # Fit the PSF components directly to the data (i.e. ignoring the spatial model)
246  if fitBasisComponents:
247  im = cand.getMaskedImage()
248 
249  im = im.Factory(im, True)
250  im.setXY0(cand.getMaskedImage().getXY0())
251 
252  try:
253  noSpatialKernel = afwMath.cast_LinearCombinationKernel(psf.getKernel())
254  except:
255  noSpatialKernel = None
256 
257  if noSpatialKernel:
258  candCenter = afwGeom.PointD(cand.getXCenter(), cand.getYCenter())
259  fit = algorithmsLib.fitKernelParamsToImage(noSpatialKernel, im, candCenter)
260  params = fit[0]
261  kernels = afwMath.KernelList(fit[1])
262  outputKernel = afwMath.LinearCombinationKernel(kernels, params)
263 
264  outImage = afwImage.ImageD(outputKernel.getDimensions())
265  outputKernel.computeImage(outImage, False)
266 
267  im -= outImage.convertF()
268  resid = im
269 
270  if margin > 0:
271  bim = im.Factory(w + 2*margin, h + 2*margin)
273  bim *= stdev
274 
275  sbim = im.Factory(bim, bbox)
276  sbim <<= resid
277  del sbim
278  resid = bim
279 
280  if variance:
281  resid = resid.getImage()
282  resid /= var
283 
284  im_resid.append(resid)
285 
286  im = im_resid.makeMosaic()
287  else:
288  im = cand.getMaskedImage()
289 
290  if normalize:
291  im /= afwMath.makeStatistics(im, afwMath.MAX).getValue()
292 
293  objId = splitId(cand.getSource().getId(), True)["objId"]
294  if psf:
295  lab = "%d chi^2 %.1f" % (objId, rchi2)
296  ctype = ds9.RED if cand.isBad() else ds9.GREEN
297  else:
298  lab = "%d flux %8.3g" % (objId, cand.getSource().getPsfFlux())
299  ctype = ds9.GREEN
300 
301  mos.append(im, lab, ctype)
302 
303  if False and numpy.isnan(rchi2):
304  ds9.mtv(cand.getMaskedImage().getImage(), title="candidate", frame=1)
305  print "amp", cand.getAmplitude()
306 
307  im = cand.getMaskedImage()
308  center = (candidateIndex, xc - im.getX0(), yc - im.getY0())
309  candidateIndex += 1
310  if cand.isBad():
311  candidateCentersBad.append(center)
312  else:
313  candidateCenters.append(center)
314 
315  if variance:
316  title = "chi(Psf fit)"
317  else:
318  title = "Stars & residuals"
319  mosaicImage = mos.makeMosaic(frame=frame, title=title)
320 
321  with ds9.Buffering():
322  for centers, color in ((candidateCenters, ds9.GREEN), (candidateCentersBad, ds9.RED)):
323  for cen in centers:
324  bbox = mos.getBBox(cen[0])
325  ds9.dot("+", cen[1] + bbox.getMinX(), cen[2] + bbox.getMinY(), frame=frame, ctype=color)
326 
327  return mosaicImage
328 
329 try:
330  import matplotlib.pyplot as plt
import matplotlib.colors
A Threshold is used to pass a threshold value to detection algorithms.
Definition: Threshold.h:44
An integer coordinate rectangle.
Definition: Box.h:53
Custom catalog class for record/table subclasses that are guaranteed to have an ID, and should generally be sorted by that ID.
Definition: fwd.h:55
A kernel that is a linear combination of fixed basis kernels.
Definition: Kernel.h:814
void randomGaussianImage(ImageT *image, Random &rand)
Definition: RandomImage.cc:155
A set of Footprints, associated with a MaskedImage.
Definition: FootprintSet.h:53
Statistics makeStatistics(afwImage::Mask< afwImage::MaskPixel > const &msk, int const flags, StatisticsControl const &sctrl)
Specialization to handle Masks.
Definition: Statistics.cc:1082
Exposure< ImagePixelT, MaskPixelT, VariancePixelT >::Ptr makeExposure(MaskedImage< ImagePixelT, MaskPixelT, VariancePixelT > &mimage, boost::shared_ptr< Wcs const > wcs=boost::shared_ptr< Wcs const >())
Definition: Exposure.h:308
std::vector< boost::shared_ptr< Kernel > > KernelList
Definition: Kernel.h:542
def lsst.meas.algorithms.utils.showPsfMosaic (   exposure,
  psf = None,
  nx = 7,
  ny = None,
  showCenter = True,
  showEllipticity = False,
  showFwhm = False,
  stampSize = 0,
  frame = 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 651 of file utils.py.

652  stampSize=0, frame=None, title=None):
653  """Show a mosaic of Psf images. exposure may be an Exposure (optionally with PSF),
654  or a tuple (width, height)
655 
656  If stampSize is > 0, the psf images will be trimmed to stampSize*stampSize
657  """
658 
659  scale = 1.0
660  if showFwhm:
661  showEllipticity = True
662  scale = 2*math.log(2) # convert sigma^2 to HWHM^2 for a Gaussian
663 
664  mos = displayUtils.Mosaic()
665 
666  try: # maybe it's a real Exposure
667  width, height = exposure.getWidth(), exposure.getHeight()
668  x0, y0 = exposure.getXY0()
669  if not psf:
670  psf = exposure.getPsf()
671  except AttributeError:
672  try: # OK, maybe a list [width, height]
673  width, height = exposure[0], exposure[1]
674  x0, y0 = 0, 0
675  except TypeError: # I guess not
676  raise RuntimeError, ("Unable to extract width/height from object of type %s" % type(exposure))
677 
678  if not ny:
679  ny = int(nx*float(height)/width + 0.5)
680  if not ny:
681  ny = 1
682 
683  centroidName = "base_GaussianCentroid"
684  shapeName = "base_SdssShape"
685 
686  schema = afwTable.SourceTable.makeMinimalSchema()
687  schema.getAliasMap().set("slot_Centroid", centroidName)
688  schema.getAliasMap().set("slot_Centroid_flag", centroidName+"_flag")
689 
690  control = measBase.GaussianCentroidControl()
691  centroider = measBase.GaussianCentroidAlgorithm(control,centroidName,schema)
692 
693  sdssShape = measBase.SdssShapeControl()
694  shaper = measBase.SdssShapeAlgorithm(sdssShape,shapeName,schema)
695  table = afwTable.SourceTable.make(schema)
696 
697  table.defineCentroid(centroidName)
698  table.defineShape(shapeName)
699 
700  bbox = None
701  if stampSize > 0:
702  w, h = psf.computeImage(afwGeom.PointD(0, 0)).getDimensions()
703  if stampSize <= w and stampSize <= h:
704  bbox = afwGeom.BoxI(afwGeom.PointI((w - stampSize)//2, (h - stampSize)//2),
705  afwGeom.ExtentI(stampSize, stampSize))
706 
707  centers = []
708  shapes = []
709  for iy in range(ny):
710  for ix in range(nx):
711  x = int(ix*(width-1)/(nx-1)) + x0
712  y = int(iy*(height-1)/(ny-1)) + y0
713 
714  im = psf.computeImage(afwGeom.PointD(x, y)).convertF()
715  imPeak = psf.computePeak(afwGeom.PointD(x, y))
716  im /= imPeak
717  if bbox:
718  im = im.Factory(im, bbox)
719  lab = "PSF(%d,%d)" % (x, y) if False else ""
720  mos.append(im, lab)
721 
723  w, h = im.getWidth(), im.getHeight()
724  centerX = im.getX0() + w//2
725  centerY = im.getY0() + h//2
726  src = table.makeRecord()
727  src.set(centroidName+"_x", centerX)
728  src.set(centroidName+"_y", centerY)
729  foot = afwDet.Footprint(exp.getBBox())
730  src.setFootprint(foot)
731 
732  centroider.measure(src, exp)
733  centers.append((src.getX() - im.getX0(), src.getY() - im.getY0()))
734 
735  shaper.measure(src, exp)
736  shapes.append((src.getIxx(), src.getIxy(), src.getIyy()))
737 
738  mos.makeMosaic(frame=frame, title=title if title else "Model Psf", mode=nx)
739 
740  if centers and frame is not None:
741  i = 0
742  with ds9.Buffering():
743  for cen, shape in zip(centers, shapes):
744  bbox = mos.getBBox(i); i += 1
745  xc, yc = cen[0] + bbox.getMinX(), cen[1] + bbox.getMinY()
746  if showCenter:
747  ds9.dot("+", xc, yc, ctype=ds9.BLUE, frame=frame)
748 
749  if showEllipticity:
750  ixx, ixy, iyy = shape
751  ixx *= scale; ixy *= scale; iyy *= scale
752  ds9.dot("@:%g,%g,%g" % (ixx, ixy, iyy), xc, yc, frame=frame, ctype=ds9.RED)
753 
754  return mos
MaskedImage< ImagePixelT, MaskPixelT, VariancePixelT > * makeMaskedImage(typename Image< ImagePixelT >::Ptr image, typename Mask< MaskPixelT >::Ptr mask=typename Mask< MaskPixelT >::Ptr(), typename Image< VariancePixelT >::Ptr variance=typename Image< VariancePixelT >::Ptr())
Definition: MaskedImage.h:1067
An integer coordinate rectangle.
Definition: Box.h:53
A set of pixels in an Image.
Definition: Footprint.h:62
Exposure< ImagePixelT, MaskPixelT, VariancePixelT >::Ptr makeExposure(MaskedImage< ImagePixelT, MaskPixelT, VariancePixelT > &mimage, boost::shared_ptr< Wcs const > wcs=boost::shared_ptr< Wcs const >())
Definition: Exposure.h:308
def lsst.meas.algorithms.utils.showPsfResiduals (   exposure,
  sourceSet,
  magType = "psf",
  scale = 10,
  frame = None,
  showAmps = False 
)

Definition at line 755 of file utils.py.

756 def showPsfResiduals(exposure, sourceSet, magType="psf", scale=10, frame=None, showAmps=False):
757  mimIn = exposure.getMaskedImage()
758  mimIn = mimIn.Factory(mimIn, True) # make a copy to subtract from
759 
760  psf = exposure.getPsf()
761  psfWidth, psfHeight = psf.getLocalKernel().getDimensions()
762  #
763  # Make the image that we'll paste our residuals into. N.b. they can overlap the edges
764  #
765  w, h = int(mimIn.getWidth()/scale), int(mimIn.getHeight()/scale)
766 
767  im = mimIn.Factory(w + psfWidth, h + psfHeight)
768 
769  cenPos = []
770  for s in sourceSet:
771  x, y = s.getX(), s.getY()
772 
773  sx, sy = int(x/scale + 0.5), int(y/scale + 0.5)
774 
775  smim = im.Factory(im, afwGeom.BoxI(afwGeom.PointI(sx, sy), afwGeom.ExtentI(psfWidth, psfHeight)))
776  sim = smim.getImage()
777 
778  try:
779  if magType == "ap":
780  flux = s.getApFlux()
781  elif magType == "model":
782  flux = s.getModelFlux()
783  elif magType == "psf":
784  flux = s.getPsfFlux()
785  else:
786  raise RuntimeError("Unknown flux type %s" % magType)
787 
788  algorithmsLib.subtractPsf(psf, mimIn, x, y, flux)
789  except Exception, e:
790  print e
791 
792  try:
793  expIm = mimIn.getImage().Factory(mimIn.getImage(),
794  afwGeom.BoxI(afwGeom.PointI(int(x) - psfWidth//2,
795  int(y) - psfHeight//2),
796  afwGeom.ExtentI(psfWidth, psfHeight)),
797  )
798  except pexExcept.Exception:
799  continue
800 
801  cenPos.append([x - expIm.getX0() + sx, y - expIm.getY0() + sy])
802 
803  sim += expIm
804 
805  if frame is not None:
806  ds9.mtv(im, frame=frame)
807  with ds9.Buffering():
808  for x, y in cenPos:
809  ds9.dot("+", x, y, frame=frame)
810 
811  if showAmps:
812  nx, ny = namp
813  for i in range(nx):
814  for j in range(ny):
815  xc = numpy.array([0, 1, 1, 0, 0])
816  yc = numpy.array([0, 0, 1, 1, 0])
817 
818  corners = []
819  for k in range(len(xc)):
820  corners.append([psfWidth//2 + w/nx*(i + xc[k]), psfHeight//2 + h/ny*(j + yc[k])])
821 
822  ds9.line(corners, frame=frame)
823 
824  return im
825 
826 #-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
An integer coordinate rectangle.
Definition: Box.h:53
def lsst.meas.algorithms.utils.showPsfSpatialCells (   exposure,
  psfCellSet,
  nMaxPerCell = -1,
  showChi2 = False,
  showMoments = False,
  symb = None,
  ctype = None,
  ctypeUnused = None,
  ctypeBad = None,
  size = 2,
  frame = None 
)
Show the SpatialCells.  If symb is something that ds9.dot understands (e.g. "o"), the
top nMaxPerCell candidates will be indicated with that symbol, using ctype and size

Definition at line 73 of file utils.py.

73 
74  symb=None, ctype=None, ctypeUnused=None, ctypeBad=None, size=2, frame=None):
75  """Show the SpatialCells. If symb is something that ds9.dot understands (e.g. "o"), the
76  top nMaxPerCell candidates will be indicated with that symbol, using ctype and size"""
77 
78  with ds9.Buffering():
79  origin = [-exposure.getMaskedImage().getX0(), -exposure.getMaskedImage().getY0()]
80  for cell in psfCellSet.getCellList():
81  displayUtils.drawBBox(cell.getBBox(), origin=origin, frame=frame)
82 
83  if nMaxPerCell < 0:
84  nMaxPerCell = 0
85 
86  i = 0
87  goodies = ctypeBad is None
88  for cand in cell.begin(goodies):
89  if nMaxPerCell > 0:
90  i += 1
91 
92  cand = algorithmsLib.cast_PsfCandidateF(cand)
93 
94  xc, yc = cand.getXCenter() + origin[0], cand.getYCenter() + origin[1]
95 
96  if i > nMaxPerCell:
97  if not ctypeUnused:
98  continue
99 
100  color = ctypeBad if cand.isBad() else ctype
101 
102  if symb:
103  if i > nMaxPerCell:
104  ct = ctypeUnused
105  else:
106  ct = ctype
107 
108  ds9.dot(symb, xc, yc, frame=frame, ctype=ct, size=size)
109 
110  source = cand.getSource()
111 
112  if showChi2:
113  rchi2 = cand.getChi2()
114  if rchi2 > 1e100:
115  rchi2 = numpy.nan
116  ds9.dot("%d %.1f" % (splitId(source.getId(), True)["objId"], rchi2),
117  xc - size, yc - size - 4, frame=frame, ctype=color, size=2)
118 
119  if showMoments:
120  ds9.dot("%.2f %.2f %.2f" % (source.getIxx(), source.getIxy(), source.getIyy()),
121  xc-size, yc + size + 4, frame=frame, ctype=color, size=size)
def lsst.meas.algorithms.utils.showSourceSet (   sSet,
  xy0 = (0, 0,
  frame = 0,
  ctype = ds9.GREEN,
  symb = "+",
  size = 2 
)
Draw the (XAstrom, YAstrom) positions of a set of Sources.  Image has the given XY0

Definition at line 56 of file utils.py.

56 
57 def showSourceSet(sSet, xy0=(0, 0), frame=0, ctype=ds9.GREEN, symb="+", size=2):
58  """Draw the (XAstrom, YAstrom) positions of a set of Sources. Image has the given XY0"""
59 
60  with ds9.Buffering():
61  for s in sSet:
62  xc, yc = s.getXAstrom() - xy0[0], s.getYAstrom() - xy0[1]
63 
64  if symb == "id":
65  ds9.dot(str(splitId(s.getId(), True)["objId"]), xc, yc, frame=frame, ctype=ctype, size=size)
66  else:
67  ds9.dot(symb, xc, yc, frame=frame, ctype=ctype, size=size)
68 
69 #-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
70 #
71 # PSF display utilities
#
def lsst.meas.algorithms.utils.splitId (   oid,
  asDict = True 
)

Definition at line 47 of file utils.py.

47 
48 def splitId(oid, asDict=True):
49 
50  objId = int((oid & 0xffff) - 1) # Should be the value set by apps code
51 
52  if asDict:
53  return dict(objId=objId)
54  else:
55  return [objId]

Variable Documentation

lsst.meas.algorithms.utils.keptPlots = False

Definition at line 42 of file utils.py.

lsst.meas.algorithms.utils.plt = None

Definition at line 332 of file utils.py.