22 """Support utilities for Measuring sources"""
25 __all__ = [
"DipoleTestImage"]
39 from .dipoleFitTask
import DipoleFitAlgorithm
40 from .
import diffimLib
41 from .
import diffimTools
43 afwDisplay.setDefaultMaskTransparency(75)
47 def showSourceSet(sSet, xy0=(0, 0), frame=0, ctype=afwDisplay.GREEN, symb=
"+", size=2):
48 """Draw the (XAstrom, YAstrom) positions of a set of Sources.
50 Image has the given XY0.
52 disp = afwDisplay.afwDisplay(frame=frame)
53 with disp.Buffering():
55 xc, yc = s.getXAstrom() - xy0[0], s.getYAstrom() - xy0[1]
58 disp.dot(str(s.getId()), xc, yc, ctype=ctype, size=size)
60 disp.dot(symb, xc, yc, ctype=ctype, size=size)
68 ctype=None, ctypeUnused=None, ctypeBad=None, size=3,
69 frame=None, title="Spatial Cells"):
70 """Show the SpatialCells.
72 If symb is something that display.dot understands (e.g. "o"), the top
73 nMaxPerCell candidates will be indicated with that symbol, using ctype
76 disp = afwDisplay.Display(frame=frame)
77 disp.mtv(maskedIm, title=title)
78 with disp.Buffering():
79 origin = [-maskedIm.getX0(), -maskedIm.getY0()]
80 for cell
in kernelCellSet.getCellList():
81 afwDisplay.utils.drawBBox(cell.getBBox(), origin=origin, display=disp)
83 goodies = ctypeBad
is None
84 for cand
in cell.begin(goodies):
85 xc, yc = cand.getXCenter() + origin[0], cand.getYCenter() + origin[1]
86 if cand.getStatus() == afwMath.SpatialCellCandidate.BAD:
88 elif cand.getStatus() == afwMath.SpatialCellCandidate.GOOD:
90 elif cand.getStatus() == afwMath.SpatialCellCandidate.UNKNOWN:
96 disp.dot(symb, xc, yc, ctype=color, size=size)
99 rchi2 = cand.getChi2()
102 disp.dot(
"%d %.1f" % (cand.getId(), rchi2),
103 xc - size, yc - size - 4, ctype=color, size=size)
107 """Display Dia Sources.
113 disp = afwDisplay.Display(frame=frame)
114 for plane
in (
"BAD",
"CR",
"EDGE",
"INTERPOlATED",
"INTRP",
"SAT",
"SATURATED"):
115 disp.setMaskPlaneColor(plane, color=
"ignore")
117 mos = afwDisplay.utils.Mosaic()
118 for i
in range(len(sources)):
120 badFlag = isFlagged[i]
121 dipoleFlag = isDipole[i]
122 bbox = source.getFootprint().getBBox()
123 stamp = exposure.Factory(exposure, bbox,
True)
124 im = afwDisplay.utils.Mosaic(gutter=1, background=0, mode=
"x")
125 im.append(stamp.getMaskedImage())
126 lab =
"%.1f,%.1f:" % (source.getX(), source.getY())
128 ctype = afwDisplay.RED
131 ctype = afwDisplay.YELLOW
133 if not badFlag
and not dipoleFlag:
134 ctype = afwDisplay.GREEN
136 mos.append(im.makeMosaic(), lab, ctype)
137 title =
"Dia Sources"
138 mosaicImage = mos.makeMosaic(display=disp, title=title)
143 resids=False, kernels=False):
144 """Display the Kernel candidates.
146 If kernel is provided include spatial model and residuals;
147 If chi is True, generate a plot of residuals/sqrt(variance), i.e. chi.
153 mos = afwDisplay.utils.Mosaic(gutter=5, background=0)
155 mos = afwDisplay.utils.Mosaic(gutter=5, background=-1)
157 candidateCenters = []
158 candidateCentersBad = []
160 for cell
in kernelCellSet.getCellList():
161 for cand
in cell.begin(
False):
164 resid = cand.getDifferenceImage(diffimLib.KernelCandidateF.ORIG)
168 rchi2 = cand.getChi2()
172 if not showBadCandidates
and cand.isBad():
175 im_resid = afwDisplay.utils.Mosaic(gutter=1, background=-0.5, mode=
"x")
178 im = cand.getScienceMaskedImage()
179 im = im.Factory(im,
True)
180 im.setXY0(cand.getScienceMaskedImage().getXY0())
183 if (
not resids
and not kernels):
184 im_resid.append(im.Factory(im,
True))
186 im = cand.getTemplateMaskedImage()
187 im = im.Factory(im,
True)
188 im.setXY0(cand.getTemplateMaskedImage().getXY0())
191 if (
not resids
and not kernels):
192 im_resid.append(im.Factory(im,
True))
196 var = resid.getVariance()
197 var = var.Factory(var,
True)
198 np.sqrt(var.getArray(), var.getArray())
199 resid = resid.getImage()
201 bbox = kernel.shrinkBBox(resid.getBBox())
202 resid = resid.Factory(resid, bbox, deep=
True)
204 kim = cand.getKernelImage(diffimLib.KernelCandidateF.ORIG).convertF()
205 resid = kim.Factory(kim,
True)
206 im_resid.append(resid)
209 ski = afwImage.ImageD(kernel.getDimensions())
210 kernel.computeImage(ski,
False, int(cand.getXCenter()), int(cand.getYCenter()))
214 sbg =
background(int(cand.getXCenter()), int(cand.getYCenter()))
215 sresid = cand.getDifferenceImage(sk, sbg)
218 resid = sresid.getImage()
220 bbox = kernel.shrinkBBox(resid.getBBox())
221 resid = resid.Factory(resid, bbox, deep=
True)
224 resid = kim.Factory(kim,
True)
225 im_resid.append(resid)
227 im = im_resid.makeMosaic()
229 lab =
"%d chi^2 %.1f" % (cand.getId(), rchi2)
230 ctype = afwDisplay.RED
if cand.isBad()
else afwDisplay.GREEN
232 mos.append(im, lab, ctype)
234 if False and np.isnan(rchi2):
235 disp = afwDisplay.Display(frame=1)
236 disp.mtv(cand.getScienceMaskedImage.getImage(), title=
"candidate")
237 print(
"rating", cand.getCandidateRating())
239 im = cand.getScienceMaskedImage()
240 center = (candidateIndex, cand.getXCenter() - im.getX0(), cand.getYCenter() - im.getY0())
243 candidateCentersBad.append(center)
245 candidateCenters.append(center)
252 title =
"Candidates & residuals"
254 disp = afwDisplay.Display(frame=frame)
255 mosaicImage = mos.makeMosaic(display=disp, title=title)
261 """Display a Kernel's basis images.
263 mos = afwDisplay.utils.Mosaic()
265 for k
in kernel.getKernelList():
266 im = afwImage.ImageD(k.getDimensions())
267 k.computeImage(im,
False)
270 disp = afwDisplay.Display(frame=frame)
271 mos.makeMosaic(display=disp, title=
"Kernel Basis Images")
279 numSample=128, keepPlots=True, maxCoeff=10):
280 """Plot the Kernel spatial model.
283 import matplotlib.pyplot
as plt
284 import matplotlib.colors
285 except ImportError
as e:
286 print(
"Unable to import numpy and matplotlib: %s" % e)
289 x0 = kernelCellSet.getBBox().getBeginX()
290 y0 = kernelCellSet.getBBox().getBeginY()
298 for cell
in kernelCellSet.getCellList():
299 for cand
in cell.begin(
False):
300 if not showBadCandidates
and cand.isBad():
302 candCenter =
geom.PointD(cand.getXCenter(), cand.getYCenter())
304 im = cand.getTemplateMaskedImage()
308 targetFits = badFits
if cand.isBad()
else candFits
309 targetPos = badPos
if cand.isBad()
else candPos
310 targetAmps = badAmps
if cand.isBad()
else candAmps
313 kp0 = np.array(cand.getKernel(diffimLib.KernelCandidateF.ORIG).getKernelParameters())
314 amp = cand.getCandidateRating()
316 targetFits = badFits
if cand.isBad()
else candFits
317 targetPos = badPos
if cand.isBad()
else candPos
318 targetAmps = badAmps
if cand.isBad()
else candAmps
320 targetFits.append(kp0)
321 targetPos.append(candCenter)
322 targetAmps.append(amp)
324 xGood = np.array([pos.getX()
for pos
in candPos]) - x0
325 yGood = np.array([pos.getY()
for pos
in candPos]) - y0
326 zGood = np.array(candFits)
328 xBad = np.array([pos.getX()
for pos
in badPos]) - x0
329 yBad = np.array([pos.getY()
for pos
in badPos]) - y0
330 zBad = np.array(badFits)
333 xRange = np.linspace(0, kernelCellSet.getBBox().getWidth(), num=numSample)
334 yRange = np.linspace(0, kernelCellSet.getBBox().getHeight(), num=numSample)
337 maxCoeff =
min(maxCoeff, kernel.getNKernelParameters())
339 maxCoeff = kernel.getNKernelParameters()
341 for k
in range(maxCoeff):
342 func = kernel.getSpatialFunction(k)
343 dfGood = zGood[:, k] - np.array([func(pos.getX(), pos.getY())
for pos
in candPos])
347 dfBad = zBad[:, k] - np.array([func(pos.getX(), pos.getY())
for pos
in badPos])
349 yMin =
min([yMin, dfBad.min()])
350 yMax =
max([yMax, dfBad.max()])
351 yMin -= 0.05*(yMax - yMin)
352 yMax += 0.05*(yMax - yMin)
354 fRange = np.ndarray((len(xRange), len(yRange)))
355 for j, yVal
in enumerate(yRange):
356 for i, xVal
in enumerate(xRange):
357 fRange[j][i] = func(xVal, yVal)
363 fig.canvas._tkcanvas._root().lift()
367 fig.suptitle(
'Kernel component %d' % k)
370 ax = fig.add_axes((0.1, 0.05, 0.35, 0.35))
373 norm = matplotlib.colors.Normalize(vmin=vmin, vmax=vmax)
374 im = ax.imshow(fRange, aspect=
'auto', norm=norm,
375 extent=[0, kernelCellSet.getBBox().getWidth() - 1,
376 0, kernelCellSet.getBBox().getHeight() - 1])
377 ax.set_title(
'Spatial polynomial')
378 plt.colorbar(im, orientation=
'horizontal', ticks=[vmin, vmax])
381 ax = fig.add_axes((0.1, 0.55, 0.35, 0.35))
382 ax.plot(-2.5*np.log10(candAmps), zGood[:, k],
'b+')
384 ax.plot(-2.5*np.log10(badAmps), zBad[:, k],
'r+')
385 ax.set_title(
"Basis Coefficients")
386 ax.set_xlabel(
"Instr mag")
387 ax.set_ylabel(
"Coeff")
390 ax = fig.add_axes((0.55, 0.05, 0.35, 0.35))
391 ax.set_autoscale_on(
False)
392 ax.set_xbound(lower=0, upper=kernelCellSet.getBBox().getHeight())
393 ax.set_ybound(lower=yMin, upper=yMax)
394 ax.plot(yGood, dfGood,
'b+')
396 ax.plot(yBad, dfBad,
'r+')
398 ax.set_title(
'dCoeff (indiv-spatial) vs. y')
401 ax = fig.add_axes((0.55, 0.55, 0.35, 0.35))
402 ax.set_autoscale_on(
False)
403 ax.set_xbound(lower=0, upper=kernelCellSet.getBBox().getWidth())
404 ax.set_ybound(lower=yMin, upper=yMax)
405 ax.plot(xGood, dfGood,
'b+')
407 ax.plot(xBad, dfBad,
'r+')
409 ax.set_title(
'dCoeff (indiv-spatial) vs. x')
414 if keepPlots
and not keptPlots:
417 print(
"%s: Please close plots when done." % __name__)
422 print(
"Plots closed, exiting...")
424 atexit.register(show)
429 """Plot the individual kernel candidate and the spatial kernel solution coefficients.
434 spatialKernel : `lsst.afw.math.LinearCombinationKernel`
435 The spatial spatialKernel solution model which is a spatially varying linear combination
436 of the spatialKernel basis functions.
437 Typically returned by `lsst.ip.diffim.SpatialKernelSolution.getSolutionPair()`.
439 kernelCellSet : `lsst.afw.math.SpatialCellSet`
440 The spatial cells that was used for solution for the spatialKernel. They contain the
441 local solutions of the AL kernel for the selected sources.
443 showBadCandidates : `bool`, optional
444 If True, plot the coefficient values for kernel candidates where the solution was marked
445 bad by the numerical algorithm. Defaults to False.
447 keepPlots: `bool`, optional
448 If True, sets ``plt.show()`` to be called before the task terminates, so that the plots
449 can be explored interactively. Defaults to True.
453 This function produces 3 figures per image subtraction operation.
454 * A grid plot of the local solutions. Each grid cell corresponds to a proportional area in
455 the image. In each cell, local kernel solution coefficients are plotted of kernel candidates (color)
456 that fall into this area as a function of the kernel basis function number.
457 * A grid plot of the spatial solution. Each grid cell corresponds to a proportional area in
458 the image. In each cell, the spatial solution coefficients are evaluated for the center of the cell.
459 * Histogram of the local solution coefficients. Red line marks the spatial solution value at
462 This function is called if ``lsst.ip.diffim.psfMatch.plotKernelCoefficients==True`` in lsstDebug. This
463 function was implemented as part of DM-17825.
466 import matplotlib.pyplot
as plt
467 except ImportError
as e:
468 print(
"Unable to import matplotlib: %s" % e)
472 imgBBox = kernelCellSet.getBBox()
473 x0 = imgBBox.getBeginX()
474 y0 = imgBBox.getBeginY()
475 wImage = imgBBox.getWidth()
476 hImage = imgBBox.getHeight()
477 imgCenterX = imgBBox.getCenterX()
478 imgCenterY = imgBBox.getCenterY()
490 fig.suptitle(
"Kernel candidate parameters on an image grid")
491 arrAx = fig.subplots(nrows=nY, ncols=nX, sharex=
True, sharey=
True, gridspec_kw=dict(
495 arrAx = arrAx[::-1, :]
498 for cell
in kernelCellSet.getCellList():
501 iX = int((cellBBox.getCenterX() - x0)//wCell)
502 iY = int((cellBBox.getCenterY() - y0)//hCell)
504 for cand
in cell.begin(
False):
506 kernel = cand.getKernel(cand.ORIG)
510 if not showBadCandidates
and cand.isBad():
513 nKernelParams = kernel.getNKernelParameters()
514 kernelParams = np.array(kernel.getKernelParameters())
515 allParams.append(kernelParams)
521 arrAx[iY, iX].
plot(np.arange(nKernelParams), kernelParams,
'.-',
522 color=color, drawstyle=
'steps-mid', linewidth=0.1)
523 for ax
in arrAx.ravel():
524 ax.grid(
True, axis=
'y')
529 spatialFuncs = spatialKernel.getSpatialFunctionList()
530 nKernelParams = spatialKernel.getNKernelParameters()
533 fig.suptitle(
"Hist. of parameters marked with spatial solution at img center")
534 arrAx = fig.subplots(nrows=int(nKernelParams//nX)+1, ncols=nX)
535 arrAx = arrAx[::-1, :]
536 allParams = np.array(allParams)
537 for k
in range(nKernelParams):
538 ax = arrAx.ravel()[k]
539 ax.hist(allParams[:, k], bins=20, edgecolor=
'black')
540 ax.set_xlabel(
'P{}'.
format(k))
541 valueParam = spatialFuncs[k](imgCenterX, imgCenterY)
542 ax.axvline(x=valueParam, color=
'red')
543 ax.text(0.1, 0.9,
'{:.1f}'.
format(valueParam),
544 transform=ax.transAxes, backgroundcolor=
'lightsteelblue')
557 fig.suptitle(
"Spatial solution of kernel parameters on an image grid")
558 arrAx = fig.subplots(nrows=nY, ncols=nX, sharex=
True, sharey=
True, gridspec_kw=dict(
560 arrAx = arrAx[::-1, :]
561 kernelParams = np.zeros(nKernelParams, dtype=float)
568 kernelParams = [f(x, y)
for f
in spatialFuncs]
569 arrAx[iY, iX].
plot(np.arange(nKernelParams), kernelParams,
'.-', drawstyle=
'steps-mid')
570 arrAx[iY, iX].grid(
True, axis=
'y')
573 if keepPlots
and not keptPlots:
576 print(
"%s: Please close plots when done." % __name__)
581 print(
"Plots closed, exiting...")
583 atexit.register(show)
588 showCenter=True, showEllipticity=True):
589 """Show a mosaic of Kernel images.
591 mos = afwDisplay.utils.Mosaic()
593 x0 = bbox.getBeginX()
594 y0 = bbox.getBeginY()
595 width = bbox.getWidth()
596 height = bbox.getHeight()
599 ny = int(nx*float(height)/width + 0.5)
603 schema = afwTable.SourceTable.makeMinimalSchema()
604 centroidName =
"base_SdssCentroid"
605 shapeName =
"base_SdssShape"
606 control = measBase.SdssCentroidControl()
607 schema.getAliasMap().
set(
"slot_Centroid", centroidName)
608 schema.getAliasMap().
set(
"slot_Centroid_flag", centroidName +
"_flag")
609 centroider = measBase.SdssCentroidAlgorithm(control, centroidName, schema)
610 sdssShape = measBase.SdssShapeControl()
611 shaper = measBase.SdssShapeAlgorithm(sdssShape, shapeName, schema)
612 table = afwTable.SourceTable.make(schema)
613 table.defineCentroid(centroidName)
614 table.defineShape(shapeName)
620 x = int(ix*(width - 1)/(nx - 1)) + x0
621 y = int(iy*(height - 1)/(ny - 1)) + y0
623 im = afwImage.ImageD(kernel.getDimensions())
624 ksum = kernel.computeImage(im,
False, x, y)
625 lab =
"Kernel(%d,%d)=%.2f" % (x, y, ksum)
if False else ""
631 w, h = im.getWidth(), im.getHeight()
632 centerX = im.getX0() + w//2
633 centerY = im.getY0() + h//2
634 src = table.makeRecord()
637 foot.addPeak(centerX, centerY, 1)
638 src.setFootprint(foot)
641 centroider.measure(src, exp)
642 centers.append((src.getX(), src.getY()))
644 shaper.measure(src, exp)
645 shapes.append((src.getIxx(), src.getIxy(), src.getIyy()))
649 disp = afwDisplay.Display(frame=frame)
650 mos.makeMosaic(display=disp, title=title
if title
else "Model Kernel", mode=nx)
652 if centers
and frame
is not None:
653 disp = afwDisplay.Display(frame=frame)
655 with disp.Buffering():
656 for cen, shape
in zip(centers, shapes):
657 bbox = mos.getBBox(i)
659 xc, yc = cen[0] + bbox.getMinX(), cen[1] + bbox.getMinY()
661 disp.dot(
"+", xc, yc, ctype=afwDisplay.BLUE)
664 ixx, ixy, iyy = shape
665 disp.dot(
"@:%g,%g,%g" % (ixx, ixy, iyy), xc, yc, ctype=afwDisplay.RED)
671 kernel, background, testSources, config,
672 origVariance=False, nptsFull=1e6, keepPlots=True, titleFs=14):
673 """Plot diffim residuals for LOCAL and SPATIAL models.
679 for cell
in kernelCellSet.getCellList():
680 for cand
in cell.begin(
True):
682 if not (cand.getStatus() == afwMath.SpatialCellCandidate.GOOD):
685 diffim = cand.getDifferenceImage(diffimLib.KernelCandidateF.ORIG)
686 orig = cand.getScienceMaskedImage()
688 ski = afwImage.ImageD(kernel.getDimensions())
689 kernel.computeImage(ski,
False, int(cand.getXCenter()), int(cand.getYCenter()))
691 sbg =
background(int(cand.getXCenter()), int(cand.getYCenter()))
692 sdiffim = cand.getDifferenceImage(sk, sbg)
695 bbox = kernel.shrinkBBox(diffim.getBBox())
696 tdiffim = diffim.Factory(diffim, bbox)
697 torig = orig.Factory(orig, bbox)
698 tsdiffim = sdiffim.Factory(sdiffim, bbox)
701 candidateResids.append(np.ravel(tdiffim.getImage().getArray()
702 / np.sqrt(torig.getVariance().getArray())))
703 spatialResids.append(np.ravel(tsdiffim.getImage().getArray()
704 / np.sqrt(torig.getVariance().getArray())))
706 candidateResids.append(np.ravel(tdiffim.getImage().getArray()
707 / np.sqrt(tdiffim.getVariance().getArray())))
708 spatialResids.append(np.ravel(tsdiffim.getImage().getArray()
709 / np.sqrt(tsdiffim.getVariance().getArray())))
711 fullIm = diffExposure.getMaskedImage().getImage().getArray()
712 fullMask = diffExposure.getMaskedImage().getMask().getArray()
714 fullVar = exposure.getMaskedImage().getVariance().getArray()
716 fullVar = diffExposure.getMaskedImage().getVariance().getArray()
719 bitmaskBad |= afwImage.Mask.getPlaneBitMask(
'NO_DATA')
720 bitmaskBad |= afwImage.Mask.getPlaneBitMask(
'SAT')
721 idx = np.where((fullMask & bitmaskBad) == 0)
722 stride = int(len(idx[0])//nptsFull)
723 sidx = idx[0][::stride], idx[1][::stride]
724 allResids = fullIm[sidx]/np.sqrt(fullVar[sidx])
726 testFootprints = diffimTools.sourceToFootprintList(testSources, warpedTemplateExposure,
727 exposure, config, Log.getDefaultLogger())
728 for fp
in testFootprints:
729 subexp = diffExposure.Factory(diffExposure, fp[
"footprint"].getBBox())
730 subim = subexp.getMaskedImage().getImage()
732 subvar = afwImage.ExposureF(exposure, fp[
"footprint"].getBBox()).getMaskedImage().getVariance()
734 subvar = subexp.getMaskedImage().getVariance()
735 nonfitResids.append(np.ravel(subim.getArray()/np.sqrt(subvar.getArray())))
737 candidateResids = np.ravel(np.array(candidateResids))
738 spatialResids = np.ravel(np.array(spatialResids))
739 nonfitResids = np.ravel(np.array(nonfitResids))
743 from matplotlib.font_manager
import FontProperties
744 except ImportError
as e:
745 print(
"Unable to import pylab: %s" % e)
751 fig.canvas._tkcanvas._root().lift()
755 fig.suptitle(
"Diffim residuals: Normalized by sqrt(input variance)", fontsize=titleFs)
757 fig.suptitle(
"Diffim residuals: Normalized by sqrt(diffim variance)", fontsize=titleFs)
759 sp1 = pylab.subplot(221)
760 sp2 = pylab.subplot(222, sharex=sp1, sharey=sp1)
761 sp3 = pylab.subplot(223, sharex=sp1, sharey=sp1)
762 sp4 = pylab.subplot(224, sharex=sp1, sharey=sp1)
763 xs = np.arange(-5, 5.05, 0.1)
764 ys = 1./np.sqrt(2*np.pi)*np.exp(-0.5*xs**2)
766 sp1.hist(candidateResids, bins=xs, normed=
True, alpha=0.5, label=
"N(%.2f, %.2f)"
767 % (np.mean(candidateResids), np.var(candidateResids)))
768 sp1.plot(xs, ys,
"r-", lw=2, label=
"N(0,1)")
769 sp1.set_title(
"Candidates: basis fit", fontsize=titleFs - 2)
770 sp1.legend(loc=1, fancybox=
True, shadow=
True, prop=FontProperties(size=titleFs - 6))
772 sp2.hist(spatialResids, bins=xs, normed=
True, alpha=0.5, label=
"N(%.2f, %.2f)"
773 % (np.mean(spatialResids), np.var(spatialResids)))
774 sp2.plot(xs, ys,
"r-", lw=2, label=
"N(0,1)")
775 sp2.set_title(
"Candidates: spatial fit", fontsize=titleFs - 2)
776 sp2.legend(loc=1, fancybox=
True, shadow=
True, prop=FontProperties(size=titleFs - 6))
778 sp3.hist(nonfitResids, bins=xs, normed=
True, alpha=0.5, label=
"N(%.2f, %.2f)"
779 % (np.mean(nonfitResids), np.var(nonfitResids)))
780 sp3.plot(xs, ys,
"r-", lw=2, label=
"N(0,1)")
781 sp3.set_title(
"Control sample: spatial fit", fontsize=titleFs - 2)
782 sp3.legend(loc=1, fancybox=
True, shadow=
True, prop=FontProperties(size=titleFs - 6))
784 sp4.hist(allResids, bins=xs, normed=
True, alpha=0.5, label=
"N(%.2f, %.2f)"
785 % (np.mean(allResids), np.var(allResids)))
786 sp4.plot(xs, ys,
"r-", lw=2, label=
"N(0,1)")
787 sp4.set_title(
"Full image (subsampled)", fontsize=titleFs - 2)
788 sp4.legend(loc=1, fancybox=
True, shadow=
True, prop=FontProperties(size=titleFs - 6))
790 pylab.setp(sp1.get_xticklabels() + sp1.get_yticklabels(), fontsize=titleFs - 4)
791 pylab.setp(sp2.get_xticklabels() + sp2.get_yticklabels(), fontsize=titleFs - 4)
792 pylab.setp(sp3.get_xticklabels() + sp3.get_yticklabels(), fontsize=titleFs - 4)
793 pylab.setp(sp4.get_xticklabels() + sp4.get_yticklabels(), fontsize=titleFs - 4)
800 if keepPlots
and not keptPlots:
803 print(
"%s: Please close plots when done." % __name__)
808 print(
"Plots closed, exiting...")
810 atexit.register(show)
815 """Calculate first moment of a (kernel) image.
819 xarr = np.asarray([[el
for el
in range(x)]
for el2
in range(y)])
820 yarr = np.asarray([[el2
for el
in range(x)]
for el2
in range(y)])
823 centx = narr.sum()/sarrSum
825 centy = narr.sum()/sarrSum
830 """Calculate second moment of a (kernel) image.
835 xarr = np.asarray([[el
for el
in range(x)]
for el2
in range(y)])
836 yarr = np.asarray([[el2
for el
in range(x)]
for el2
in range(y)])
837 narr = sarr*np.power((xarr - centx), 2.)
839 xstd = np.sqrt(narr.sum()/sarrSum)
840 narr = sarr*np.power((yarr - centy), 2.)
841 ystd = np.sqrt(narr.sum()/sarrSum)
846 """Print differences in sky coordinates.
848 The difference is that between the source Position and its Centroid mapped
852 sCentroid = s.getCentroid()
853 sPosition = s.getCoord().getPosition(geom.degrees)
854 dra = 3600*(sPosition.getX() - wcs.pixelToSky(sCentroid).getPosition(geom.degrees).getX())/0.2
855 ddec = 3600*(sPosition.getY() - wcs.pixelToSky(sCentroid).getPosition(geom.degrees).getY())/0.2
856 if np.isfinite(dra)
and np.isfinite(ddec):
861 """Create regions file for display from input source list.
863 fh = open(outfilename,
"w")
864 fh.write(
"global color=red font=\"helvetica 10 normal\" "
865 "select=1 highlite=1 edit=1 move=1 delete=1 include=1 fixed=0 source\nfk5\n")
868 (ra, dec) = wcs.pixelToSky(s.getCentroid()).getPosition(geom.degrees)
870 (ra, dec) = s.getCoord().getPosition(geom.degrees)
871 if np.isfinite(ra)
and np.isfinite(dec):
872 fh.write(
"circle(%f,%f,2\")\n"%(ra, dec))
878 """Draw the (RA, Dec) positions of a set of Sources. Image has the XY0.
880 disp = afwDisplay.Display(frame=frame)
881 with disp.Buffering():
883 (xc, yc) = wcs.skyToPixel(s.getCoord().getRa(), s.getCoord().getDec())
886 disp.dot(symb, xc, yc, ctype=ctype, size=size)
890 """Plot whisker diagram of astromeric offsets between results.matches.
892 refCoordKey = results.matches[0].first.getTable().getCoordKey()
893 inCentroidKey = results.matches[0].second.getTable().getCentroidKey()
894 positions = [m.first.get(refCoordKey)
for m
in results.matches]
895 residuals = [m.first.get(refCoordKey).getOffsetFrom(
896 newWcs.pixelToSky(m.second.get(inCentroidKey)))
for
897 m
in results.matches]
898 import matplotlib.pyplot
as plt
900 sp = fig.add_subplot(1, 1, 0)
901 xpos = [x[0].asDegrees()
for x
in positions]
902 ypos = [x[1].asDegrees()
for x
in positions]
903 xpos.append(0.02*(
max(xpos) -
min(xpos)) +
min(xpos))
904 ypos.append(0.98*(
max(ypos) -
min(ypos)) +
min(ypos))
905 xidxs = np.isfinite(xpos)
906 yidxs = np.isfinite(ypos)
907 X = np.asarray(xpos)[xidxs]
908 Y = np.asarray(ypos)[yidxs]
909 distance = [x[1].asArcseconds()
for x
in residuals]
911 distance = np.asarray(distance)[xidxs]
914 bearing = [x[0].asRadians()
for x
in residuals]
916 bearing = np.asarray(bearing)[xidxs]
917 U = (distance*np.cos(bearing))
918 V = (distance*np.sin(bearing))
919 sp.quiver(X, Y, U, V)
920 sp.set_title(
"WCS Residual")
925 """Utility class for dipole measurement testing.
927 Generate an image with simulated dipoles and noise; store the original
928 "pre-subtraction" images and catalogs as well.
929 Used to generate test data for DMTN-007 (http://dmtn-007.lsst.io).
932 def __init__(self, w=101, h=101, xcenPos=[27.], ycenPos=[25.], xcenNeg=[23.], ycenNeg=[25.],
933 psfSigma=2., flux=[30000.], fluxNeg=None, noise=10., gradientParams=None):
949 def _makeDipoleImage(self):
950 """Generate an exposure and catalog with the given dipole source(s).
959 dipole = posImage.clone()
960 di = dipole.getMaskedImage()
961 di -= negImage.getMaskedImage()
965 posDetectedBits = posImage.getMaskedImage().getMask().getArray() == dm.getPlaneBitMask(
"DETECTED")
966 negDetectedBits = negImage.getMaskedImage().getMask().getArray() == dm.getPlaneBitMask(
"DETECTED")
967 pos_det = dm.addMaskPlane(
"DETECTED_POS")
968 neg_det = dm.addMaskPlane(
"DETECTED_NEG")
971 dma[:, :] = posDetectedBits*pos_det + negDetectedBits*neg_det
972 self.diffim, self.posImage, self.posCatalog, self.negImage, self.negCatalog \
973 = dipole, posImage, posCatalog, negImage, negCatalog
975 def _makeStarImage(self, xc=[15.3], yc=[18.6], flux=[2500], schema=None, randomSeed=None):
976 """Generate an exposure and catalog with the given stellar source(s).
980 dataset = TestDataset(bbox, psfSigma=self.
psfSigma, threshold=1.)
982 for i
in range(len(xc)):
983 dataset.addSource(instFlux=flux[i], centroid=
geom.Point2D(xc[i], yc[i]))
986 schema = TestDataset.makeMinimalSchema()
987 exposure, catalog = dataset.realize(noise=self.
noise, schema=schema, randomSeed=randomSeed)
990 y, x = np.mgrid[:self.
w, :self.
h]
992 gradient = gp[0] + gp[1]*x + gp[2]*y
994 gradient += gp[3]*x*y + gp[4]*x*x + gp[5]*y*y
995 imgArr = exposure.getMaskedImage().getArrays()[0]
998 return exposure, catalog
1002 fitResult = alg.fitDipole(source, **kwds)
1006 """Utility function for detecting dipoles.
1008 Detect pos/neg sources in the diffim, then merge them. A
1009 bigger "grow" parameter leads to a larger footprint which
1010 helps with dipole measurement for faint dipoles.
1015 Whether to merge the positive and negagive detections into a single
1017 diffim : `lsst.afw.image.exposure.exposure.ExposureF`
1018 Difference image on which to perform detection.
1019 detectSigma : `float`
1020 Threshold for object detection.
1022 Number of pixels to grow the footprints before merging.
1024 Minimum bin size for the background (re)estimation (only applies if
1025 the default leads to min(nBinX, nBinY) < fit order so the default
1026 config parameter needs to be decreased, but not to a value smaller
1027 than ``minBinSize``, in which case the fitting algorithm will take
1028 over and decrease the fit order appropriately.)
1032 sources : `lsst.afw.table.SourceCatalog`
1033 If doMerge=True, the merged source catalog is returned OR
1034 detectTask : `lsst.meas.algorithms.SourceDetectionTask`
1035 schema : `lsst.afw.table.Schema`
1036 If doMerge=False, the source detection task and its schema are
1040 diffim = self.diffim
1043 schema = afwTable.SourceTable.makeMinimalSchema()
1046 detectConfig = measAlg.SourceDetectionConfig()
1047 detectConfig.returnOriginalFootprints =
False
1049 psfSigma = diffim.getPsf().computeShape().getDeterminantRadius()
1052 detectConfig.thresholdPolarity =
"both"
1053 detectConfig.thresholdValue = detectSigma
1055 detectConfig.reEstimateBackground =
True
1056 detectConfig.thresholdType =
"pixel_stdev"
1058 while ((
min(diffim.getWidth(), diffim.getHeight()))//detectConfig.background.binSize
1059 < detectConfig.background.approxOrderX
and detectConfig.background.binSize > minBinSize):
1060 detectConfig.background.binSize =
max(minBinSize, detectConfig.background.binSize//2)
1063 detectTask = measAlg.SourceDetectionTask(schema, config=detectConfig)
1065 table = afwTable.SourceTable.make(schema)
1066 catalog = detectTask.run(table, diffim, sigma=psfSigma)
1070 fpSet = catalog.fpSets.positive
1071 fpSet.merge(catalog.fpSets.negative, grow, grow,
False)
1073 fpSet.makeSources(sources)
1078 return detectTask, schema