23 """Support utilities for Measuring sources"""
49 objId = int((oid & 0xffff) - 1)
52 return dict(objId=objId)
56 def showSourceSet(sSet, xy0=(0, 0), frame=0, ctype=ds9.GREEN, symb=
"+", size=2):
57 """Draw the (XAstrom, YAstrom) positions of a set of Sources. Image has the given XY0"""
61 xc, yc = s.getXAstrom() - xy0[0], s.getYAstrom() - xy0[1]
64 ds9.dot(str(
splitId(s.getId(),
True)[
"objId"]), xc, yc, frame=frame, ctype=ctype, size=size)
66 ds9.dot(symb, xc, yc, frame=frame, ctype=ctype, size=size)
72 def showPsfSpatialCells(exposure, psfCellSet, nMaxPerCell=-1, showChi2=False, showMoments=False,
73 symb=
None, ctype=
None, ctypeUnused=
None, ctypeBad=
None, size=2, frame=
None):
74 """Show the SpatialCells. If symb is something that ds9.dot understands (e.g. "o"), the
75 top nMaxPerCell candidates will be indicated with that symbol, using ctype and size"""
78 origin = [-exposure.getMaskedImage().getX0(), -exposure.getMaskedImage().getY0()]
79 for cell
in psfCellSet.getCellList():
80 displayUtils.drawBBox(cell.getBBox(), origin=origin, frame=frame)
86 goodies = ctypeBad
is None
87 for cand
in cell.begin(goodies):
91 cand = algorithmsLib.cast_PsfCandidateF(cand)
93 xc, yc = cand.getXCenter() + origin[0], cand.getYCenter() + origin[1]
99 color = ctypeBad
if cand.isBad()
else ctype
107 ds9.dot(symb, xc, yc, frame=frame, ctype=ct, size=size)
109 source = cand.getSource()
112 rchi2 = cand.getChi2()
115 ds9.dot(
"%d %.1f" % (
splitId(source.getId(),
True)[
"objId"], rchi2),
116 xc - size, yc - size - 4, frame=frame, ctype=color, size=2)
119 ds9.dot(
"%.2f %.2f %.2f" % (source.getIxx(), source.getIxy(), source.getIyy()),
120 xc-size, yc + size + 4, frame=frame, ctype=color, size=size)
122 def showPsfCandidates(exposure, psfCellSet, psf=None, frame=None, normalize=True, showBadCandidates=True,
123 fitBasisComponents=
False, variance=
None, chi=
None):
124 """Display the PSF candidates.
126 If psf is provided include PSF model and residuals; if normalize is true normalize the PSFs
129 If chi is True, generate a plot of residuals/sqrt(variance), i.e. chi
131 If fitBasisComponents is true, also find the best linear combination of the PSF's components
135 if variance
is not None:
140 mos = displayUtils.Mosaic()
142 candidateCenters = []
143 candidateCentersBad = []
146 for cell
in psfCellSet.getCellList():
147 for cand
in cell.begin(
False):
148 cand = algorithmsLib.cast_PsfCandidateF(cand)
150 rchi2 = cand.getChi2()
154 if not showBadCandidates
and cand.isBad():
158 im_resid = displayUtils.Mosaic(gutter=0, background=-5, mode=
"x")
161 im = cand.getMaskedImage()
162 xc, yc = cand.getXCenter(), cand.getYCenter()
164 margin = 0
if True else 5
165 w, h = im.getDimensions()
169 bim = im.Factory(w + 2*margin, h + 2*margin)
174 var = bim.getVariance(); var.set(stdev**2); del var
176 sbim = im.Factory(bim, bbox)
180 xc += margin; yc += margin
182 im = im.Factory(im,
True)
183 im.setXY0(cand.getMaskedImage().getXY0())
188 im_resid.append(im.Factory(im,
True))
192 psfIm = mi.getImage()
193 config = measBase.SingleFrameMeasurementTask.ConfigClass()
194 config.slots.centroid =
"base_SdssCentroid"
196 schema = afwTable.SourceTable.makeMinimalSchema()
197 measureSources = measBase.SingleFrameMeasurementTask(schema, config=config)
201 miBig = mi.Factory(im.getWidth() + 2*extra, im.getHeight() + 2*extra)
202 miBig[extra:-extra, extra:-extra] = mi
203 miBig.setXY0(mi.getX0() - extra, mi.getY0() - extra)
204 mi = miBig; del miBig
212 footprintSet.makeSources(catalog)
214 if len(catalog) == 0:
215 raise RuntimeError(
"Failed to detect any objects")
217 measureSources.run(catalog, exp)
218 if len(catalog) == 1:
221 for i, s
in enumerate(catalog):
222 d = numpy.hypot(xc - s.getX(), yc - s.getY())
223 if i == 0
or d < dmin:
225 xc, yc = source.getCentroid()
229 chi2 = algorithmsLib.subtractPsf(psf, im, xc, yc)
236 resid = resid.getImage()
237 var = im.getVariance()
238 var = var.Factory(var,
True)
239 numpy.sqrt(var.getArray(), var.getArray())
242 im_resid.append(resid)
245 if fitBasisComponents:
246 im = cand.getMaskedImage()
248 im = im.Factory(im,
True)
249 im.setXY0(cand.getMaskedImage().getXY0())
252 noSpatialKernel = afwMath.cast_LinearCombinationKernel(psf.getKernel())
254 noSpatialKernel =
None
258 fit = algorithmsLib.fitKernelParamsToImage(noSpatialKernel, im, candCenter)
263 outImage = afwImage.ImageD(outputKernel.getDimensions())
264 outputKernel.computeImage(outImage,
False)
266 im -= outImage.convertF()
270 bim = im.Factory(w + 2*margin, h + 2*margin)
274 sbim = im.Factory(bim, bbox)
280 resid = resid.getImage()
283 im_resid.append(resid)
285 im = im_resid.makeMosaic()
287 im = cand.getMaskedImage()
292 objId =
splitId(cand.getSource().getId(),
True)[
"objId"]
294 lab =
"%d chi^2 %.1f" % (objId, rchi2)
295 ctype = ds9.RED
if cand.isBad()
else ds9.GREEN
297 lab =
"%d flux %8.3g" % (objId, cand.getSource().getPsfFlux())
300 mos.append(im, lab, ctype)
302 if False and numpy.isnan(rchi2):
303 ds9.mtv(cand.getMaskedImage().getImage(), title=
"candidate", frame=1)
304 print "amp", cand.getAmplitude()
306 im = cand.getMaskedImage()
307 center = (candidateIndex, xc - im.getX0(), yc - im.getY0())
310 candidateCentersBad.append(center)
312 candidateCenters.append(center)
315 title =
"chi(Psf fit)"
317 title =
"Stars & residuals"
318 mosaicImage = mos.makeMosaic(frame=frame, title=title)
320 with ds9.Buffering():
321 for centers, color
in ((candidateCenters, ds9.GREEN), (candidateCentersBad, ds9.RED)):
323 bbox = mos.getBBox(cen[0])
324 ds9.dot(
"+", cen[1] + bbox.getMinX(), cen[2] + bbox.getMinY(), frame=frame, ctype=color)
329 import matplotlib.pyplot
as plt
330 import matplotlib.colors
334 def makeSubplots(fig, nx=2, ny=2, Nx=1, Ny=1, plottingArea=(0.1, 0.1, 0.85, 0.80),
335 pxgutter=0.05, pygutter=0.05, xgutter=0.04, ygutter=0.04,
336 headroom=0.0, panelBorderWeight=0, panelColor=
'black'):
337 """Return a generator of a set of subplots, a set of Nx*Ny panels of nx*ny plots. Each panel is fully
338 filled by row (starting in the bottom left) before the next panel is started. If panelBorderWidth is
339 greater than zero a border is drawn around each panel, adjusted to enclose the axis labels.
342 subplots = makeSubplots(fig, 2, 2, Nx=1, Ny=1, panelColor='k')
343 ax = subplots.next(); ax.text(0.3, 0.5, '[0, 0] (0,0)')
344 ax = subplots.next(); ax.text(0.3, 0.5, '[0, 0] (1,0)')
345 ax = subplots.next(); ax.text(0.3, 0.5, '[0, 0] (0,1)')
346 ax = subplots.next(); ax.text(0.3, 0.5, '[0, 0] (1,1)')
349 @param fig The matplotlib figure to draw
350 @param nx The number of plots in each row of each panel
351 @param ny The number of plots in each column of each panel
352 @param Nx The number of panels in each row of the figure
353 @param Ny The number of panels in each column of the figure
354 @param plottingArea (x0, y0, x1, y1) for the part of the figure containing all the panels
355 @param pxgutter Spacing between columns of panels in units of (x1 - x0)
356 @param pygutter Spacing between rows of panels in units of (y1 - y0)
357 @param xgutter Spacing between columns of plots within a panel in units of (x1 - x0)
358 @param ygutter Spacing between rows of plots within a panel in units of (y1 - y0)
359 @param headroom Extra spacing above each plot for e.g. a title
360 @param panelBorderWeight Width of border drawn around panels
361 @param panelColor Colour of border around panels
367 except AttributeError:
368 fig.__show = fig.show
374 fig.show = types.MethodType(myShow, fig, fig.__class__)
382 Callback to draw the panel borders when the plots are drawn to the canvas
384 if panelBorderWeight <= 0:
387 for p
in axes.keys():
390 bboxes.append(ax.bbox.union([label.get_window_extent()
for label
in
391 ax.get_xticklabels() + ax.get_yticklabels()]))
398 bbox = ax.bbox.union(bboxes)
400 xy0, xy1 = ax.transData.inverted().transform(bbox)
401 x0, y0 = xy0; x1, y1 = xy1
402 w, h = x1 - x0, y1 - y0
404 x0 -= 0.02*w; w += 0.04*w
405 y0 -= 0.02*h; h += 0.04*h
409 rec = ax.add_patch(plt.Rectangle((x0, y0), w, h, fill=
False,
410 lw=panelBorderWeight, edgecolor=panelColor))
411 rec.set_clip_on(
False)
415 fig.canvas.mpl_connect(
'draw_event', on_draw)
419 x0, y0 = plottingArea[0:2]
420 W, H = plottingArea[2:4]
421 w = (W - (Nx - 1)*pxgutter - (nx*Nx - 1)*xgutter)/float(nx*Nx)
422 h = (H - (Ny - 1)*pygutter - (ny*Ny - 1)*ygutter)/float(ny*Ny)
426 for panel
in range(Nx*Ny):
430 for window
in range(nx*ny):
431 x = nx*px + window%nx
432 y = ny*py + window//nx
433 ax = fig.add_axes((x0 + xgutter + pxgutter + x*w + (px - 1)*pxgutter + (x - 1)*xgutter,
434 y0 + ygutter + pygutter + y*h + (py - 1)*pygutter + (y - 1)*ygutter,
435 w, h), frame_on=
True, axis_bgcolor=
'w')
436 axes[panel].append(ax)
440 matchKernelAmplitudes=
False, keepPlots=
True):
441 """Plot the PSF spatial model."""
444 print >> sys.stderr,
"Unable to import matplotlib"
447 noSpatialKernel = afwMath.cast_LinearCombinationKernel(psf.getKernel())
454 for cell
in psfCellSet.getCellList():
455 for cand
in cell.begin(
False):
456 cand = algorithmsLib.cast_PsfCandidateF(cand)
457 if not showBadCandidates
and cand.isBad():
461 im = cand.getMaskedImage()
465 fit = algorithmsLib.fitKernelParamsToImage(noSpatialKernel, im, candCenter)
469 for p, k
in zip(params, kernels):
470 amp += p * afwMath.cast_FixedKernel(k).getSum()
472 targetFits = badFits
if cand.isBad()
else candFits
473 targetPos = badPos
if cand.isBad()
else candPos
474 targetAmps = badAmps
if cand.isBad()
else candAmps
476 targetFits.append([x / amp
for x
in params])
477 targetPos.append(candCenter)
478 targetAmps.append(amp)
480 numCandidates = len(candFits)
481 numBasisFuncs = noSpatialKernel.getNBasisKernels()
483 xGood = numpy.array([pos.getX()
for pos
in candPos]) - exposure.getX0()
484 yGood = numpy.array([pos.getY()
for pos
in candPos]) - exposure.getY0()
485 zGood = numpy.array(candFits)
486 ampGood = numpy.array(candAmps)
488 xBad = numpy.array([pos.getX()
for pos
in badPos]) - exposure.getX0()
489 yBad = numpy.array([pos.getY()
for pos
in badPos]) - exposure.getY0()
490 zBad = numpy.array(badFits)
491 ampBad = numpy.array(badAmps)
494 xRange = numpy.linspace(0, exposure.getWidth(), num=numSample)
495 yRange = numpy.linspace(0, exposure.getHeight(), num=numSample)
497 kernel = psf.getKernel()
498 nKernelComponents = kernel.getNKernelParameters()
502 nPanelX = int(math.sqrt(nKernelComponents))
503 nPanelY = nKernelComponents//nPanelX
504 while nPanelY*nPanelX < nKernelComponents:
510 fig.canvas._tkcanvas._root().lift()
516 subplots =
makeSubplots(fig, 2, 2, Nx=nPanelX, Ny=nPanelY, xgutter=0.06, ygutter=0.06, pygutter=0.04)
518 for k
in range(nKernelComponents):
519 func = kernel.getSpatialFunction(k)
520 dfGood = zGood[:,k] - numpy.array([func(pos.getX(), pos.getY())
for pos
in candPos])
524 dfBad = zBad[:,k] - numpy.array([func(pos.getX(), pos.getY())
for pos
in badPos])
525 yMin = min([yMin, dfBad.min()])
526 yMax = max([yMax, dfBad.max()])
527 yMin -= 0.05 * (yMax - yMin)
528 yMax += 0.05 * (yMax - yMin)
533 fRange = numpy.ndarray((len(xRange), len(yRange)))
534 for j, yVal
in enumerate(yRange):
535 for i, xVal
in enumerate(xRange):
536 fRange[j][i] = func(xVal, yVal)
542 ax.set_autoscale_on(
False)
543 ax.set_xbound(lower=0, upper=exposure.getHeight())
544 ax.set_ybound(lower=yMin, upper=yMax)
545 ax.plot(yGood, dfGood,
'b+')
547 ax.plot(yBad, dfBad,
'r+')
549 ax.set_title(
'Residuals(y)')
555 if matchKernelAmplitudes
and k == 0:
562 norm = matplotlib.colors.Normalize(vmin=vmin, vmax=vmax)
563 im = ax.imshow(fRange, aspect=
'auto', origin=
"lower", norm=norm,
564 extent=[0, exposure.getWidth()-1, 0, exposure.getHeight()-1])
565 ax.set_title(
'Spatial poly')
566 plt.colorbar(im, orientation=
'horizontal', ticks=[vmin, vmax])
571 ax.set_autoscale_on(
False)
572 ax.set_xbound(lower=0, upper=exposure.getWidth())
573 ax.set_ybound(lower=yMin, upper=yMax)
574 ax.plot(xGood, dfGood,
'b+')
576 ax.plot(xBad, dfBad,
'r+')
578 ax.set_title(
'K%d Residuals(x)' % k)
585 ax.scatter(xGood, yGood, c=dfGood, marker=
'o')
586 ax.scatter(xBad, yBad, c=dfBad, marker=
'x')
587 ax.set_xbound(lower=0, upper=exposure.getWidth())
588 ax.set_ybound(lower=0, upper=exposure.getHeight())
589 ax.set_title(
'Spatial residuals')
590 plt.colorbar(im, orientation=
'horizontal')
592 calib = exposure.getCalib()
593 if calib.getFluxMag0()[0] <= 0:
594 calib = type(calib)()
595 calib.setFluxMag0(1.0)
598 ax.plot(calib.getMagnitude(candAmps), zGood[:,k],
'b+')
600 ax.plot(calib.getMagnitude(badAmps), zBad[:,k],
'r+')
602 ax.set_title(
'Flux variation')
607 if keepPlots
and not keptPlots:
610 print "%s: Please close plots when done." % __name__
615 print "Plots closed, exiting..."
617 atexit.register(show)
620 def showPsf(psf, eigenValues=None, XY=None, normalize=True, frame=None):
621 """Display a PSF's eigen images
623 If normalize is True, set the largest absolute value of each eigenimage to 1.0 (n.b. sum == 0.0 for i > 0)
629 coeffs = psf.getLocalKernel(
afwGeom.PointD(XY[0], XY[1])).getKernelParameters()
633 mos = displayUtils.Mosaic(gutter=2, background=-0.1)
634 for i, k
in enumerate(afwMath.cast_LinearCombinationKernel(psf.getKernel()).getKernelList()):
635 im = afwImage.ImageD(k.getDimensions())
636 k.computeImage(im,
False)
638 im /= numpy.max(numpy.abs(im.getArray()))
641 mos.append(im,
"%g" % (coeffs[i]/coeffs[0]))
645 mos.makeMosaic(frame=frame, title=
"Kernel Basis Functions")
650 showCenter=
True, showEllipticity=
False, showFwhm=
False,
651 stampSize=0, frame=
None, title=
None):
652 """Show a mosaic of Psf images. exposure may be an Exposure (optionally with PSF),
653 or a tuple (width, height)
655 If stampSize is > 0, the psf images will be trimmed to stampSize*stampSize
660 showEllipticity =
True
661 scale = 2*math.log(2)
663 mos = displayUtils.Mosaic()
666 width, height = exposure.getWidth(), exposure.getHeight()
667 x0, y0 = exposure.getXY0()
669 psf = exposure.getPsf()
670 except AttributeError:
672 width, height = exposure[0], exposure[1]
675 raise RuntimeError, (
"Unable to extract width/height from object of type %s" % type(exposure))
678 ny = int(nx*float(height)/width + 0.5)
682 centroidName =
"base_GaussianCentroid"
683 shapeName =
"base_SdssShape"
685 schema = afwTable.SourceTable.makeMinimalSchema()
686 schema.getAliasMap().set(
"slot_Centroid", centroidName)
687 schema.getAliasMap().set(
"slot_Centroid_flag", centroidName+
"_flag")
689 control = measBase.GaussianCentroidControl()
690 centroider = measBase.GaussianCentroidAlgorithm(control,centroidName,schema)
692 sdssShape = measBase.SdssShapeControl()
693 shaper = measBase.SdssShapeAlgorithm(sdssShape,shapeName,schema)
694 table = afwTable.SourceTable.make(schema)
696 table.defineCentroid(centroidName)
697 table.defineShape(shapeName)
702 if stampSize <= w
and stampSize <= h:
710 x = int(ix*(width-1)/(nx-1)) + x0
711 y = int(iy*(height-1)/(ny-1)) + y0
717 im = im.Factory(im, bbox)
718 lab =
"PSF(%d,%d)" % (x, y)
if False else ""
722 w, h = im.getWidth(), im.getHeight()
723 centerX = im.getX0() + w//2
724 centerY = im.getY0() + h//2
725 src = table.makeRecord()
726 src.set(centroidName+
"_x", centerX)
727 src.set(centroidName+
"_y", centerY)
729 src.setFootprint(foot)
731 centroider.measure(src, exp)
732 centers.append((src.getX() - im.getX0(), src.getY() - im.getY0()))
734 shaper.measure(src, exp)
735 shapes.append((src.getIxx(), src.getIxy(), src.getIyy()))
737 mos.makeMosaic(frame=frame, title=title
if title
else "Model Psf", mode=nx)
739 if centers
and frame
is not None:
741 with ds9.Buffering():
742 for cen, shape
in zip(centers, shapes):
743 bbox = mos.getBBox(i); i += 1
744 xc, yc = cen[0] + bbox.getMinX(), cen[1] + bbox.getMinY()
746 ds9.dot(
"+", xc, yc, ctype=ds9.BLUE, frame=frame)
749 ixx, ixy, iyy = shape
750 ixx *= scale; ixy *= scale; iyy *= scale
751 ds9.dot(
"@:%g,%g,%g" % (ixx, ixy, iyy), xc, yc, frame=frame, ctype=ds9.RED)
755 def showPsfResiduals(exposure, sourceSet, magType="psf", scale=10, frame=None, showAmps=False):
756 mimIn = exposure.getMaskedImage()
757 mimIn = mimIn.Factory(mimIn,
True)
759 psf = exposure.getPsf()
760 psfWidth, psfHeight = psf.getLocalKernel().getDimensions()
764 w, h = int(mimIn.getWidth()/scale), int(mimIn.getHeight()/scale)
766 im = mimIn.Factory(w + psfWidth, h + psfHeight)
770 x, y = s.getX(), s.getY()
772 sx, sy = int(x/scale + 0.5), int(y/scale + 0.5)
775 sim = smim.getImage()
780 elif magType ==
"model":
781 flux = s.getModelFlux()
782 elif magType ==
"psf":
783 flux = s.getPsfFlux()
785 raise RuntimeError(
"Unknown flux type %s" % magType)
787 algorithmsLib.subtractPsf(psf, mimIn, x, y, flux)
792 expIm = mimIn.getImage().Factory(mimIn.getImage(),
794 int(y) - psfHeight//2),
800 cenPos.append([x - expIm.getX0() + sx, y - expIm.getY0() + sy])
804 if frame
is not None:
805 ds9.mtv(im, frame=frame)
806 with ds9.Buffering():
808 ds9.dot(
"+", x, y, frame=frame)
814 xc = numpy.array([0, 1, 1, 0, 0])
815 yc = numpy.array([0, 0, 1, 1, 0])
818 for k
in range(len(xc)):
819 corners.append([psfWidth//2 + w/nx*(i + xc[k]), psfHeight//2 + h/ny*(j + yc[k])])
821 ds9.line(corners, frame=frame)
828 """Write the contents of a SpatialCellSet to a many-MEF fits file"""
831 for cell
in psfCellSet.getCellList():
832 for cand
in cell.begin(
False):
833 cand = algorithmsLib.cast_PsfCandidateF(cand)
840 md.set(
"CELL", cell.getLabel())
841 md.set(
"ID", cand.getId())
842 md.set(
"XCENTER", cand.getXCenter())
843 md.set(
"YCENTER", cand.getYCenter())
844 md.set(
"BAD", cand.isBad())
845 md.set(
"AMPL", cand.getAmplitude())
846 md.set(
"FLUX", cand.getSource().getPsfFlux())
847 md.set(
"CHI2", cand.getSource().getChi2())
849 im.writeFits(fileName, md, mode)
852 if frame
is not None:
853 ds9.mtv(im, frame=frame)
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())
std::pair< int, double > positionToIndex(double const pos, bool)
Convert image position to index (nearest integer and fractional parts)
A Threshold is used to pass a threshold value to detection algorithms.
An integer coordinate rectangle.
Custom catalog class for record/table subclasses that are guaranteed to have an ID, and should generally be sorted by that ID.
A kernel that is a linear combination of fixed basis kernels.
void randomGaussianImage(ImageT *image, Random &rand)
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.
Class for storing generic metadata.
Statistics makeStatistics(afwImage::Mask< afwImage::MaskPixel > const &msk, int const flags, StatisticsControl const &sctrl)
Specialization to handle Masks.
Exposure< ImagePixelT, MaskPixelT, VariancePixelT >::Ptr makeExposure(MaskedImage< ImagePixelT, MaskPixelT, VariancePixelT > &mimage, boost::shared_ptr< Wcs const > wcs=boost::shared_ptr< Wcs const >())
std::vector< boost::shared_ptr< Kernel > > KernelList