22"""Support utilities for Measuring sources"""
24__all__ = [
"splitId",
"showSourceSet",
"showPsfSpatialCells",
25 "showPsfCandidates",
"makeSubplots",
"plotPsfSpatialModel",
26 "showPsf",
"showPsfMosaic",
"showPsfResiduals",
"saveSpatialCellSet"]
43from .
import subtractPsf, fitKernelParamsToImage
47afwDisplay.setDefaultMaskTransparency(75)
49_LOG = logging.getLogger(__name__)
54 objId = int((oid & 0xffff) - 1)
57 return dict(objId=objId)
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"""
66 display = afwDisplay.Display()
67 with display.Buffering():
69 xc, yc = s.getXAstrom() - xy0[0], s.getYAstrom() - xy0[1]
72 display.dot(str(
splitId(s.getId(),
True)[
"objId"]), xc, yc, ctype=ctype, size=size)
74 display.dot(symb, xc, yc, ctype=ctype, size=size)
82 symb=None, ctype=None, ctypeUnused=None, ctypeBad=None, size=2, display=None):
83 """Show the SpatialCells.
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
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)
101 goodies = ctypeBad
is None
102 for cand
in cell.begin(goodies):
106 xc, yc = cand.getXCenter() + origin[0], cand.getYCenter() + origin[1]
112 color = ctypeBad
if cand.isBad()
else ctype
120 display.dot(symb, xc, yc, ctype=ct, size=size)
122 source = cand.getSource()
125 rchi2 = cand.getChi2()
128 display.dot(
"%d %.1f" % (
splitId(source.getId(),
True)[
"objId"], rchi2),
129 xc - size, yc - size - 4, ctype=color, size=2)
132 display.dot(
"%.2f %.2f %.2f" % (source.getIxx(), source.getIxy(), source.getIyy()),
133 xc-size, yc + size + 4, ctype=color, size=size)
137def showPsfCandidates(exposure, psfCellSet, psf=None, display=None, normalize=True, showBadCandidates=True,
138 fitBasisComponents=False, variance=None, chi=None):
139 """Display the PSF candidates.
141 If psf is provided include PSF model and residuals; if normalize is true normalize the PSFs
144 If chi is True, generate a plot of residuals/sqrt(variance), i.e. chi
146 If fitBasisComponents is true, also find the best linear combination of the PSF's components
150 display = afwDisplay.Display()
153 if variance
is not None:
158 mos = displayUtils.Mosaic()
160 candidateCenters = []
161 candidateCentersBad = []
164 for cell
in psfCellSet.getCellList():
165 for cand
in cell.begin(
False):
166 rchi2 = cand.getChi2()
170 if not showBadCandidates
and cand.isBad():
174 im_resid = displayUtils.Mosaic(gutter=0, background=-5, mode=
"x")
177 im = cand.getMaskedImage()
178 xc, yc = cand.getXCenter(), cand.getYCenter()
180 margin = 0
if True else 5
181 w, h = im.getDimensions()
185 bim = im.Factory(w + 2*margin, h + 2*margin)
189 bim.getVariance().
set(stdev**2)
196 im = im.Factory(im,
True)
197 im.setXY0(cand.getMaskedImage().getXY0())
202 im_resid.append(im.Factory(im,
True))
206 psfIm = mi.getImage()
207 config = measBase.SingleFrameMeasurementTask.ConfigClass()
208 config.slots.centroid =
"base_SdssCentroid"
210 schema = afwTable.SourceTable.makeMinimalSchema()
211 measureSources = measBase.SingleFrameMeasurementTask(schema, config=config)
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)
227 footprintSet.makeSources(catalog)
229 if len(catalog) == 0:
230 raise RuntimeError(
"Failed to detect any objects")
232 measureSources.run(catalog, exp)
233 if len(catalog) == 1:
237 for i, s
in enumerate(catalog):
238 d = numpy.hypot(xc - s.getX(), yc - s.getY())
239 if i == 0
or d < dmin:
241 xc, yc = source.getCentroid()
245 subtractPsf(psf, im, xc, yc)
251 resid = resid.getImage()
252 var = im.getVariance()
253 var = var.Factory(var,
True)
254 numpy.sqrt(var.getArray(), var.getArray())
257 im_resid.append(resid)
260 if fitBasisComponents:
261 im = cand.getMaskedImage()
263 im = im.Factory(im,
True)
264 im.setXY0(cand.getMaskedImage().getXY0())
267 noSpatialKernel = psf.getKernel()
269 noSpatialKernel =
None
273 fit = fitKernelParamsToImage(noSpatialKernel, im, candCenter)
278 outImage = afwImage.ImageD(outputKernel.getDimensions())
279 outputKernel.computeImage(outImage,
False)
281 im -= outImage.convertF()
285 bim = im.Factory(w + 2*margin, h + 2*margin)
289 bim.assign(resid, bbox)
293 resid = resid.getImage()
296 im_resid.append(resid)
298 im = im_resid.makeMosaic()
300 im = cand.getMaskedImage()
305 objId =
splitId(cand.getSource().getId(),
True)[
"objId"]
307 lab =
"%d chi^2 %.1f" % (objId, rchi2)
308 ctype = afwDisplay.RED
if cand.isBad()
else afwDisplay.GREEN
310 lab =
"%d flux %8.3g" % (objId, cand.getSource().getPsfInstFlux())
311 ctype = afwDisplay.GREEN
313 mos.append(im, lab, ctype)
315 if False and numpy.isnan(rchi2):
316 display.mtv(cand.getMaskedImage().getImage(), title=
"showPsfCandidates: candidate")
317 print(
"amp", cand.getAmplitude())
319 im = cand.getMaskedImage()
320 center = (candidateIndex, xc - im.getX0(), yc - im.getY0())
323 candidateCentersBad.append(center)
325 candidateCenters.append(center)
328 title =
"chi(Psf fit)"
330 title =
"Stars & residuals"
331 mosaicImage = mos.makeMosaic(display=display, title=title)
333 with display.Buffering():
334 for centers, color
in ((candidateCenters, afwDisplay.GREEN), (candidateCentersBad, afwDisplay.RED)):
336 bbox = mos.getBBox(cen[0])
337 display.dot(
"+", cen[1] + bbox.getMinX(), cen[2] + bbox.getMinY(), ctype=color)
342def makeSubplots(fig, nx=2, ny=2, Nx=1, Ny=1, plottingArea=(0.1, 0.1, 0.85, 0.80),
343 pxgutter=0.05, pygutter=0.05, xgutter=0.04, ygutter=0.04,
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.
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)')
359 fig : `matplotlib.pyplot.figure`
360 The matplotlib figure to draw
362 The number of plots in each row of each panel
364 The number of plots in each column of each panel
366 The number of panels in each row of the figure
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
372 Spacing between columns of panels in units of (x1 - x0)
374 Spacing between rows of panels in units of (y1 - y0)
376 Spacing between columns of plots within a panel in units of (x1 - x0)
378 Spacing between rows of plots within a panel in units of (y1 - y0)
380 Extra spacing above each plot for e.g. a title
381 panelBorderWeight : `int`
382 Width of border drawn around panels
384 Colour of border around panels
387 log = _LOG.getChild(
"makeSubplots")
389 import matplotlib.pyplot
as plt
390 except ImportError
as e:
391 log.warning(
"Unable to import matplotlib: %s", e)
397 except AttributeError:
398 fig.__show = fig.show
405 fig.show = types.MethodType(myShow, fig)
414 Callback to draw the panel borders when the plots are drawn to the canvas
416 if panelBorderWeight <= 0:
419 for p
in axes.keys():
422 bboxes.append(ax.bbox.union([label.get_window_extent()
for label
in
423 ax.get_xticklabels() + ax.get_yticklabels()]))
430 bbox = ax.bbox.union(bboxes)
432 xy0, xy1 = ax.transData.inverted().transform(bbox)
435 w, h = x1 - x0, y1 - y0
444 rec = ax.add_patch(plt.Rectangle((x0, y0), w, h, fill=
False,
445 lw=panelBorderWeight, edgecolor=panelColor))
446 rec.set_clip_on(
False)
450 fig.canvas.mpl_connect(
'draw_event', on_draw)
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)
461 for panel
in range(Nx*Ny):
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)
476 matchKernelAmplitudes=False, keepPlots=True):
477 """Plot the PSF spatial model."""
479 log = _LOG.getChild(
"plotPsfSpatialModel")
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)
487 noSpatialKernel = psf.getKernel()
494 for cell
in psfCellSet.getCellList():
495 for cand
in cell.begin(
False):
496 if not showBadCandidates
and cand.isBad():
500 im = cand.getMaskedImage()
504 fit = fitKernelParamsToImage(noSpatialKernel, im, candCenter)
508 for p, k
in zip(params, kernels):
509 amp += p * k.getSum()
511 targetFits = badFits
if cand.isBad()
else candFits
512 targetPos = badPos
if cand.isBad()
else candPos
513 targetAmps = badAmps
if cand.isBad()
else candAmps
515 targetFits.append([x / amp
for x
in params])
516 targetPos.append(candCenter)
517 targetAmps.append(amp)
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)
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)
528 xRange = numpy.linspace(0, exposure.getWidth(), num=numSample)
529 yRange = numpy.linspace(0, exposure.getHeight(), num=numSample)
531 kernel = psf.getKernel()
532 nKernelComponents = kernel.getNKernelParameters()
536 nPanelX = int(math.sqrt(nKernelComponents))
537 nPanelY = nKernelComponents//nPanelX
538 while nPanelY*nPanelX < nKernelComponents:
544 fig.canvas._tkcanvas._root().lift()
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)
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])
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)
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)
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+')
580 ax.plot(yBad, dfBad,
'r+')
582 ax.set_title(
'Residuals(y)')
586 if matchKernelAmplitudes
and k == 0:
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])
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+')
605 ax.plot(xBad, dfBad,
'r+')
607 ax.set_title(
'K%d Residuals(x)' % k)
611 photoCalib = exposure.getPhotoCalib()
613 if photoCalib.getCalibrationMean() <= 0:
616 ampMag = [photoCalib.instFluxToMagnitude(candAmp)
for candAmp
in candAmps]
617 ax.plot(ampMag, zGood[:, k],
'b+')
619 badAmpMag = [photoCalib.instFluxToMagnitude(badAmp)
for badAmp
in badAmps]
620 ax.plot(badAmpMag, zBad[:, k],
'r+')
622 ax.set_title(
'Flux variation')
627 if keepPlots
and not keptPlots:
630 print(
"%s: Please close plots when done." % __name__)
635 print(
"Plots closed, exiting...")
637 atexit.register(show)
641def showPsf(psf, eigenValues=None, XY=None, normalize=True, display=None):
642 """Display a PSF's eigen images
644 If normalize is True, set the largest absolute value of each eigenimage to 1.0 (n.b. sum == 0.0 for i > 0)
650 coeffs = psf.getLocalKernel(
lsst.geom.PointD(XY[0], XY[1])).getKernelParameters()
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)
659 im /= numpy.max(numpy.abs(im.getArray()))
662 mos.append(im,
"%g" % (coeffs[i]/coeffs[0]))
667 display = afwDisplay.Display()
668 mos.makeMosaic(display=display, title=
"Kernel Basis Functions")
673def showPsfMosaic(exposure, psf=None, nx=7, ny=None, showCenter=True, showEllipticity=False,
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)
678 If stampSize is > 0, the psf images will be trimmed to stampSize*stampSize
683 showEllipticity =
True
684 scale = 2*math.log(2)
686 mos = displayUtils.Mosaic()
689 width, height = exposure.getWidth(), exposure.getHeight()
690 x0, y0 = exposure.getXY0()
692 psf = exposure.getPsf()
693 except AttributeError:
695 width, height = exposure[0], exposure[1]
698 raise RuntimeError(
"Unable to extract width/height from object of type %s" % type(exposure))
701 ny = int(nx*float(height)/width + 0.5)
705 centroidName =
"SdssCentroid"
706 shapeName =
"base_SdssShape"
708 schema = afwTable.SourceTable.makeMinimalSchema()
709 schema.getAliasMap().
set(
"slot_Centroid", centroidName)
710 schema.getAliasMap().
set(
"slot_Centroid_flag", centroidName+
"_flag")
712 control = measBase.SdssCentroidControl()
713 centroider = measBase.SdssCentroidAlgorithm(control, centroidName, schema)
715 sdssShape = measBase.SdssShapeControl()
716 shaper = measBase.SdssShapeAlgorithm(sdssShape, shapeName, schema)
717 table = afwTable.SourceTable.make(schema)
719 table.defineCentroid(centroidName)
720 table.defineShape(shapeName)
725 if stampSize <= w
and stampSize <= h:
733 x = int(ix*(width-1)/(nx-1)) + x0
734 y = int(iy*(height-1)/(ny-1)) + y0
740 im = im.Factory(im, bbox)
741 lab =
"PSF(%d,%d)" % (x, y)
if False else ""
746 w, h = im.getWidth(), im.getHeight()
747 centerX = im.getX0() + w//2
748 centerY = im.getY0() + h//2
749 src = table.makeRecord()
752 foot.addPeak(centerX, centerY, 1)
753 src.setFootprint(foot)
756 centroider.measure(src, exp)
757 centers.append((src.getX() - im.getX0(), src.getY() - im.getY0()))
759 shaper.measure(src, exp)
760 shapes.append((src.getIxx(), src.getIxy(), src.getIyy()))
765 display = afwDisplay.Display()
766 mos.makeMosaic(display=display, title=title
if title
else "Model Psf", mode=nx)
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()
774 display.dot(
"+", xc, yc, ctype=afwDisplay.BLUE)
777 ixx, ixy, iyy = shape
781 display.dot(
"@:%g,%g,%g" % (ixx, ixy, iyy), xc, yc, ctype=afwDisplay.RED)
787 mimIn = exposure.getMaskedImage()
788 mimIn = mimIn.Factory(mimIn,
True)
790 psf = exposure.getPsf()
791 psfWidth, psfHeight = psf.getLocalKernel(psf.getAveragePosition()).getDimensions()
795 w, h = int(mimIn.getWidth()/scale), int(mimIn.getHeight()/scale)
797 im = mimIn.Factory(w + psfWidth, h + psfHeight)
801 x, y = s.getX(), s.getY()
803 sx, sy = int(x/scale + 0.5), int(y/scale + 0.5)
807 sim = smim.getImage()
811 flux = s.getApInstFlux()
812 elif magType ==
"model":
813 flux = s.getModelInstFlux()
814 elif magType ==
"psf":
815 flux = s.getPsfInstFlux()
817 raise RuntimeError(
"Unknown flux type %s" % magType)
819 subtractPsf(psf, mimIn, x, y, flux)
820 except Exception
as e:
824 expIm = mimIn.getImage().Factory(mimIn.getImage(),
826 int(y) - psfHeight//2),
832 cenPos.append([x - expIm.getX0() + sx, y - expIm.getY0() + sy])
837 display = afwDisplay.Display()
838 display.mtv(im, title=
"showPsfResiduals: image")
839 with display.Buffering():
841 display.dot(
"+", x, y)
847 """Write the contents of a SpatialCellSet to a many-MEF fits file"""
850 for cell
in psfCellSet.getCellList():
851 for cand
in cell.begin(
False):
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())
866 im.writeFits(fileName, md, mode)
870 display.mtv(im, title=
"saveSpatialCellSet: image")
A Threshold is used to pass a threshold value to detection algorithms.
A compact representation of a collection of pixels.
The photometric calibration of an exposure.
A kernel that is a linear combination of fixed basis kernels.
A class that can be used to generate sequences of random numbers according to a number of different a...
Class for storing generic metadata.
An integer coordinate rectangle.
Provides consistent interface for LSST exceptions.
daf::base::PropertySet * set
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.
int positionToIndex(double pos)
Convert image position to nearest integer index.
std::shared_ptr< Exposure< ImagePixelT, MaskPixelT, VariancePixelT > > makeExposure(MaskedImage< ImagePixelT, MaskPixelT, VariancePixelT > &mimage, std::shared_ptr< geom::SkyWcs const > wcs=std::shared_ptr< geom::SkyWcs const >())
A function to return an Exposure of the correct type (cf.
void randomGaussianImage(ImageT *image, Random &rand)
Set image to random numbers with a gaussian N(0, 1) distribution.
Statistics makeStatistics(lsst::afw::image::Image< Pixel > const &img, lsst::afw::image::Mask< image::MaskPixel > const &msk, int const flags, StatisticsControl const &sctrl=StatisticsControl())
Handle a watered-down front-end to the constructor (no variance)
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.
showPsfCandidates(exposure, psfCellSet, psf=None, display=None, normalize=True, showBadCandidates=True, fitBasisComponents=False, variance=None, chi=None)
showPsfMosaic(exposure, psf=None, nx=7, ny=None, showCenter=True, showEllipticity=False, showFwhm=False, stampSize=0, display=None, title=None)
showPsfSpatialCells(exposure, psfCellSet, nMaxPerCell=-1, showChi2=False, showMoments=False, symb=None, ctype=None, ctypeUnused=None, ctypeBad=None, size=2, display=None)
plotPsfSpatialModel(exposure, psf, psfCellSet, showBadCandidates=True, numSample=128, matchKernelAmplitudes=False, keepPlots=True)
showSourceSet(sSet, xy0=(0, 0), display=None, ctype=afwDisplay.GREEN, symb="+", size=2)
saveSpatialCellSet(psfCellSet, fileName="foo.fits", display=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')
showPsfResiduals(exposure, sourceSet, magType="psf", scale=10, display=None)
showPsf(psf, eigenValues=None, XY=None, normalize=True, display=None)
splitId(oid, asDict=True)