22"""Support utilities for Measuring sources"""
25__all__ = [
"DipoleTestImage",
"evaluateMeanPsfFwhm",
"getPsfFwhm",
"getKernelCenterDisplacement"]
41from lsst.utils.logging
import getLogger
42from .dipoleFitTask
import DipoleFitAlgorithm
43from .
import diffimLib
45afwDisplay.setDefaultMaskTransparency(75)
48_LOG = getLogger(__name__)
51def showSourceSet(sSet, xy0=(0, 0), frame=0, ctype=afwDisplay.GREEN, symb=
"+", size=2):
52 """Draw the (XAstrom, YAstrom) positions of a set of Sources.
54 Image has the given XY0.
56 disp = afwDisplay.afwDisplay(frame=frame)
57 with disp.Buffering():
59 xc, yc = s.getXAstrom() - xy0[0], s.getYAstrom() - xy0[1]
62 disp.dot(str(s.getId()), xc, yc, ctype=ctype, size=size)
64 disp.dot(symb, xc, yc, ctype=ctype, size=size)
72 ctype=None, ctypeUnused=None, ctypeBad=None, size=3,
73 frame=None, title="Spatial Cells"):
74 """Show the SpatialCells.
76 If symb is something that display.dot understands (e.g. "o"), the top
77 nMaxPerCell candidates will be indicated with that symbol, using ctype
80 disp = afwDisplay.Display(frame=frame)
81 disp.mtv(maskedIm, title=title)
82 with disp.Buffering():
83 origin = [-maskedIm.getX0(), -maskedIm.getY0()]
84 for cell
in kernelCellSet.getCellList():
85 afwDisplay.utils.drawBBox(cell.getBBox(), origin=origin, display=disp)
87 goodies = ctypeBad
is None
88 for cand
in cell.begin(goodies):
89 xc, yc = cand.getXCenter() + origin[0], cand.getYCenter() + origin[1]
90 if cand.getStatus() == afwMath.SpatialCellCandidate.BAD:
92 elif cand.getStatus() == afwMath.SpatialCellCandidate.GOOD:
94 elif cand.getStatus() == afwMath.SpatialCellCandidate.UNKNOWN:
100 disp.dot(symb, xc, yc, ctype=color, size=size)
103 rchi2 = cand.getChi2()
106 disp.dot(
"%d %.1f" % (cand.getId(), rchi2),
107 xc - size, yc - size - 4, ctype=color, size=size)
111 """Display Dia Sources.
117 disp = afwDisplay.Display(frame=frame)
118 for plane
in (
"BAD",
"CR",
"EDGE",
"INTERPOlATED",
"INTRP",
"SAT",
"SATURATED"):
119 disp.setMaskPlaneColor(plane, color=
"ignore")
121 mos = afwDisplay.utils.Mosaic()
122 for i
in range(len(sources)):
124 badFlag = isFlagged[i]
125 dipoleFlag = isDipole[i]
126 bbox = source.getFootprint().getBBox()
127 stamp = exposure.Factory(exposure, bbox,
True)
128 im = afwDisplay.utils.Mosaic(gutter=1, background=0, mode=
"x")
129 im.append(stamp.getMaskedImage())
130 lab =
"%.1f,%.1f:" % (source.getX(), source.getY())
132 ctype = afwDisplay.RED
135 ctype = afwDisplay.YELLOW
137 if not badFlag
and not dipoleFlag:
138 ctype = afwDisplay.GREEN
140 mos.append(im.makeMosaic(), lab, ctype)
141 title =
"Dia Sources"
142 mosaicImage = mos.makeMosaic(display=disp, title=title)
147 resids=False, kernels=False):
148 """Display the Kernel candidates.
150 If kernel is provided include spatial model and residuals;
151 If chi is True, generate a plot of residuals/sqrt(variance), i.e. chi.
157 mos = afwDisplay.utils.Mosaic(gutter=5, background=0)
159 mos = afwDisplay.utils.Mosaic(gutter=5, background=-1)
161 candidateCenters = []
162 candidateCentersBad = []
164 for cell
in kernelCellSet.getCellList():
165 for cand
in cell.begin(
False):
168 resid = cand.getDifferenceImage(diffimLib.KernelCandidateF.ORIG)
172 rchi2 = cand.getChi2()
176 if not showBadCandidates
and cand.isBad():
179 im_resid = afwDisplay.utils.Mosaic(gutter=1, background=-0.5, mode=
"x")
182 im = cand.getScienceMaskedImage()
183 im = im.Factory(im,
True)
184 im.setXY0(cand.getScienceMaskedImage().getXY0())
187 if (
not resids
and not kernels):
188 im_resid.append(im.Factory(im,
True))
190 im = cand.getTemplateMaskedImage()
191 im = im.Factory(im,
True)
192 im.setXY0(cand.getTemplateMaskedImage().getXY0())
195 if (
not resids
and not kernels):
196 im_resid.append(im.Factory(im,
True))
201 var = var.Factory(var,
True)
202 np.sqrt(var.array, var.array)
205 bbox = kernel.shrinkBBox(resid.getBBox())
206 resid = resid.Factory(resid, bbox, deep=
True)
208 kim = cand.getKernelImage(diffimLib.KernelCandidateF.ORIG).convertF()
209 resid = kim.Factory(kim,
True)
210 im_resid.append(resid)
213 ski = afwImage.ImageD(kernel.getDimensions())
214 kernel.computeImage(ski,
False, int(cand.getXCenter()), int(cand.getYCenter()))
218 sbg = background(int(cand.getXCenter()), int(cand.getYCenter()))
219 sresid = cand.getDifferenceImage(sk, sbg)
224 bbox = kernel.shrinkBBox(resid.getBBox())
225 resid = resid.Factory(resid, bbox, deep=
True)
228 resid = kim.Factory(kim,
True)
229 im_resid.append(resid)
231 im = im_resid.makeMosaic()
233 lab =
"%d chi^2 %.1f" % (cand.getId(), rchi2)
234 ctype = afwDisplay.RED
if cand.isBad()
else afwDisplay.GREEN
236 mos.append(im, lab, ctype)
238 if False and np.isnan(rchi2):
239 disp = afwDisplay.Display(frame=1)
240 disp.mtv(cand.getScienceMaskedImage.image, title=
"candidate")
241 print(
"rating", cand.getCandidateRating())
243 im = cand.getScienceMaskedImage()
244 center = (candidateIndex, cand.getXCenter() - im.getX0(), cand.getYCenter() - im.getY0())
247 candidateCentersBad.append(center)
249 candidateCenters.append(center)
256 title =
"Candidates & residuals"
258 disp = afwDisplay.Display(frame=frame)
259 mosaicImage = mos.makeMosaic(display=disp, title=title)
265 """Display a Kernel's basis images.
267 mos = afwDisplay.utils.Mosaic()
269 for k
in kernel.getKernelList():
270 im = afwImage.ImageD(k.getDimensions())
271 k.computeImage(im,
False)
274 disp = afwDisplay.Display(frame=frame)
275 mos.makeMosaic(display=disp, title=
"Kernel Basis Images")
283 numSample=128, keepPlots=True, maxCoeff=10):
284 """Plot the Kernel spatial model.
287 import matplotlib.pyplot
as plt
288 import matplotlib.colors
289 except ImportError
as e:
290 print(
"Unable to import numpy and matplotlib: %s" % e)
293 x0 = kernelCellSet.getBBox().getBeginX()
294 y0 = kernelCellSet.getBBox().getBeginY()
302 for cell
in kernelCellSet.getCellList():
303 for cand
in cell.begin(
False):
304 if not showBadCandidates
and cand.isBad():
306 candCenter =
geom.PointD(cand.getXCenter(), cand.getYCenter())
308 im = cand.getTemplateMaskedImage()
312 targetFits = badFits
if cand.isBad()
else candFits
313 targetPos = badPos
if cand.isBad()
else candPos
314 targetAmps = badAmps
if cand.isBad()
else candAmps
317 kp0 = np.array(cand.getKernel(diffimLib.KernelCandidateF.ORIG).getKernelParameters())
318 amp = cand.getCandidateRating()
320 targetFits = badFits
if cand.isBad()
else candFits
321 targetPos = badPos
if cand.isBad()
else candPos
322 targetAmps = badAmps
if cand.isBad()
else candAmps
324 targetFits.append(kp0)
325 targetPos.append(candCenter)
326 targetAmps.append(amp)
328 xGood = np.array([pos.getX()
for pos
in candPos]) - x0
329 yGood = np.array([pos.getY()
for pos
in candPos]) - y0
330 zGood = np.array(candFits)
332 xBad = np.array([pos.getX()
for pos
in badPos]) - x0
333 yBad = np.array([pos.getY()
for pos
in badPos]) - y0
334 zBad = np.array(badFits)
337 xRange = np.linspace(0, kernelCellSet.getBBox().getWidth(), num=numSample)
338 yRange = np.linspace(0, kernelCellSet.getBBox().getHeight(), num=numSample)
341 maxCoeff =
min(maxCoeff, kernel.getNKernelParameters())
343 maxCoeff = kernel.getNKernelParameters()
345 for k
in range(maxCoeff):
346 func = kernel.getSpatialFunction(k)
347 dfGood = zGood[:, k] - np.array([func(pos.getX(), pos.getY())
for pos
in candPos])
351 dfBad = zBad[:, k] - np.array([func(pos.getX(), pos.getY())
for pos
in badPos])
353 yMin =
min([yMin, dfBad.min()])
354 yMax =
max([yMax, dfBad.max()])
355 yMin -= 0.05*(yMax - yMin)
356 yMax += 0.05*(yMax - yMin)
358 fRange = np.ndarray((len(xRange), len(yRange)))
359 for j, yVal
in enumerate(yRange):
360 for i, xVal
in enumerate(xRange):
361 fRange[j][i] = func(xVal, yVal)
367 fig.canvas._tkcanvas._root().lift()
371 fig.suptitle(
'Kernel component %d' % k)
374 ax = fig.add_axes((0.1, 0.05, 0.35, 0.35))
377 norm = matplotlib.colors.Normalize(vmin=vmin, vmax=vmax)
378 im = ax.imshow(fRange, aspect=
'auto', norm=norm,
379 extent=[0, kernelCellSet.getBBox().getWidth() - 1,
380 0, kernelCellSet.getBBox().getHeight() - 1])
381 ax.set_title(
'Spatial polynomial')
382 plt.colorbar(im, orientation=
'horizontal', ticks=[vmin, vmax])
385 ax = fig.add_axes((0.1, 0.55, 0.35, 0.35))
386 ax.plot(-2.5*np.log10(candAmps), zGood[:, k],
'b+')
388 ax.plot(-2.5*np.log10(badAmps), zBad[:, k],
'r+')
389 ax.set_title(
"Basis Coefficients")
390 ax.set_xlabel(
"Instr mag")
391 ax.set_ylabel(
"Coeff")
394 ax = fig.add_axes((0.55, 0.05, 0.35, 0.35))
395 ax.set_autoscale_on(
False)
396 ax.set_xbound(lower=0, upper=kernelCellSet.getBBox().getHeight())
397 ax.set_ybound(lower=yMin, upper=yMax)
398 ax.plot(yGood, dfGood,
'b+')
400 ax.plot(yBad, dfBad,
'r+')
402 ax.set_title(
'dCoeff (indiv-spatial) vs. y')
405 ax = fig.add_axes((0.55, 0.55, 0.35, 0.35))
406 ax.set_autoscale_on(
False)
407 ax.set_xbound(lower=0, upper=kernelCellSet.getBBox().getWidth())
408 ax.set_ybound(lower=yMin, upper=yMax)
409 ax.plot(xGood, dfGood,
'b+')
411 ax.plot(xBad, dfBad,
'r+')
413 ax.set_title(
'dCoeff (indiv-spatial) vs. x')
418 if keepPlots
and not keptPlots:
421 print(
"%s: Please close plots when done." % __name__)
426 print(
"Plots closed, exiting...")
428 atexit.register(show)
433 """Plot the individual kernel candidate and the spatial kernel solution coefficients.
438 spatialKernel : `lsst.afw.math.LinearCombinationKernel`
439 The spatial spatialKernel solution model which is a spatially varying linear combination
440 of the spatialKernel basis functions.
441 Typically returned by `lsst.ip.diffim.SpatialKernelSolution.getSolutionPair()`.
443 kernelCellSet : `lsst.afw.math.SpatialCellSet`
444 The spatial cells that was used for solution for the spatialKernel. They contain the
445 local solutions of the AL kernel for the selected sources.
447 showBadCandidates : `bool`, optional
448 If True, plot the coefficient values for kernel candidates where the solution was marked
449 bad by the numerical algorithm. Defaults to False.
451 keepPlots: `bool`, optional
452 If True, sets ``plt.show()`` to be called before the task terminates, so that the plots
453 can be explored interactively. Defaults to True.
457 This function produces 3 figures per image subtraction operation.
458 * A grid plot of the local solutions. Each grid cell corresponds to a proportional area in
459 the image. In each cell, local kernel solution coefficients are plotted of kernel candidates (color)
460 that fall into this area as a function of the kernel basis function number.
461 * A grid plot of the spatial solution. Each grid cell corresponds to a proportional area in
462 the image. In each cell, the spatial solution coefficients are evaluated for the center of the cell.
463 * Histogram of the local solution coefficients. Red line marks the spatial solution value at
466 This function is called if ``lsst.ip.diffim.psfMatch.plotKernelCoefficients==True`` in lsstDebug. This
467 function was implemented as part of DM-17825.
470 import matplotlib.pyplot
as plt
471 except ImportError
as e:
472 print(
"Unable to import matplotlib: %s" % e)
476 imgBBox = kernelCellSet.getBBox()
477 x0 = imgBBox.getBeginX()
478 y0 = imgBBox.getBeginY()
479 wImage = imgBBox.getWidth()
480 hImage = imgBBox.getHeight()
481 imgCenterX = imgBBox.getCenterX()
482 imgCenterY = imgBBox.getCenterY()
494 fig.suptitle(
"Kernel candidate parameters on an image grid")
495 arrAx = fig.subplots(nrows=nY, ncols=nX, sharex=
True, sharey=
True, gridspec_kw=dict(
499 arrAx = arrAx[::-1, :]
502 for cell
in kernelCellSet.getCellList():
505 iX = int((cellBBox.getCenterX() - x0)//wCell)
506 iY = int((cellBBox.getCenterY() - y0)//hCell)
508 for cand
in cell.begin(
False):
510 kernel = cand.getKernel(cand.ORIG)
514 if not showBadCandidates
and cand.isBad():
517 nKernelParams = kernel.getNKernelParameters()
518 kernelParams = np.array(kernel.getKernelParameters())
519 allParams.append(kernelParams)
525 arrAx[iY, iX].plot(np.arange(nKernelParams), kernelParams,
'.-',
526 color=color, drawstyle=
'steps-mid', linewidth=0.1)
527 for ax
in arrAx.ravel():
528 ax.grid(
True, axis=
'y')
533 spatialFuncs = spatialKernel.getSpatialFunctionList()
534 nKernelParams = spatialKernel.getNKernelParameters()
537 fig.suptitle(
"Hist. of parameters marked with spatial solution at img center")
538 arrAx = fig.subplots(nrows=int(nKernelParams//nX)+1, ncols=nX)
539 arrAx = arrAx[::-1, :]
540 allParams = np.array(allParams)
541 for k
in range(nKernelParams):
542 ax = arrAx.ravel()[k]
543 ax.hist(allParams[:, k], bins=20, edgecolor=
'black')
544 ax.set_xlabel(
'P{}'.format(k))
545 valueParam = spatialFuncs[k](imgCenterX, imgCenterY)
546 ax.axvline(x=valueParam, color=
'red')
547 ax.text(0.1, 0.9,
'{:.1f}'.format(valueParam),
548 transform=ax.transAxes, backgroundcolor=
'lightsteelblue')
561 fig.suptitle(
"Spatial solution of kernel parameters on an image grid")
562 arrAx = fig.subplots(nrows=nY, ncols=nX, sharex=
True, sharey=
True, gridspec_kw=dict(
564 arrAx = arrAx[::-1, :]
565 kernelParams = np.zeros(nKernelParams, dtype=float)
572 kernelParams = [f(x, y)
for f
in spatialFuncs]
573 arrAx[iY, iX].plot(np.arange(nKernelParams), kernelParams,
'.-', drawstyle=
'steps-mid')
574 arrAx[iY, iX].grid(
True, axis=
'y')
577 if keepPlots
and not keptPlots:
580 print(
"%s: Please close plots when done." % __name__)
585 print(
"Plots closed, exiting...")
587 atexit.register(show)
592 showCenter=True, showEllipticity=True):
593 """Show a mosaic of Kernel images.
595 mos = afwDisplay.utils.Mosaic()
597 x0 = bbox.getBeginX()
598 y0 = bbox.getBeginY()
599 width = bbox.getWidth()
600 height = bbox.getHeight()
603 ny = int(nx*float(height)/width + 0.5)
607 schema = afwTable.SourceTable.makeMinimalSchema()
608 centroidName =
"base_SdssCentroid"
609 shapeName =
"base_SdssShape"
610 control = measBase.SdssCentroidControl()
611 schema.getAliasMap().set(
"slot_Centroid", centroidName)
612 schema.getAliasMap().set(
"slot_Centroid_flag", centroidName +
"_flag")
613 centroider = measBase.SdssCentroidAlgorithm(control, centroidName, schema)
614 sdssShape = measBase.SdssShapeControl()
615 shaper = measBase.SdssShapeAlgorithm(sdssShape, shapeName, schema)
616 table = afwTable.SourceTable.make(schema)
617 table.defineCentroid(centroidName)
618 table.defineShape(shapeName)
624 x = int(ix*(width - 1)/(nx - 1)) + x0
625 y = int(iy*(height - 1)/(ny - 1)) + y0
627 im = afwImage.ImageD(kernel.getDimensions())
628 ksum = kernel.computeImage(im,
False, x, y)
629 lab =
"Kernel(%d,%d)=%.2f" % (x, y, ksum)
if False else ""
635 w, h = im.getWidth(), im.getHeight()
636 centerX = im.getX0() + w//2
637 centerY = im.getY0() + h//2
638 src = table.makeRecord()
641 foot.addPeak(centerX, centerY, 1)
642 src.setFootprint(foot)
645 centroider.measure(src, exp)
646 centers.append((src.getX(), src.getY()))
648 shaper.measure(src, exp)
649 shapes.append((src.getIxx(), src.getIxy(), src.getIyy()))
653 disp = afwDisplay.Display(frame=frame)
654 mos.makeMosaic(display=disp, title=title
if title
else "Model Kernel", mode=nx)
656 if centers
and frame
is not None:
657 disp = afwDisplay.Display(frame=frame)
659 with disp.Buffering():
660 for cen, shape
in zip(centers, shapes):
661 bbox = mos.getBBox(i)
663 xc, yc = cen[0] + bbox.getMinX(), cen[1] + bbox.getMinY()
665 disp.dot(
"+", xc, yc, ctype=afwDisplay.BLUE)
668 ixx, ixy, iyy = shape
669 disp.dot(
"@:%g,%g,%g" % (ixx, ixy, iyy), xc, yc, ctype=afwDisplay.RED)
675 """Plot whisker diagram of astromeric offsets between results.matches.
677 refCoordKey = results.matches[0].first.getTable().getCoordKey()
678 inCentroidKey = results.matches[0].second.getTable().getCentroidSlot().getMeasKey()
679 positions = [m.first.get(refCoordKey)
for m
in results.matches]
680 residuals = [m.first.get(refCoordKey).getOffsetFrom(
681 newWcs.pixelToSky(m.second.get(inCentroidKey)))
for
682 m
in results.matches]
683 import matplotlib.pyplot
as plt
685 sp = fig.add_subplot(1, 1, 0)
686 xpos = [x[0].asDegrees()
for x
in positions]
687 ypos = [x[1].asDegrees()
for x
in positions]
688 xpos.append(0.02*(
max(xpos) -
min(xpos)) +
min(xpos))
689 ypos.append(0.98*(
max(ypos) -
min(ypos)) +
min(ypos))
690 xidxs = np.isfinite(xpos)
691 yidxs = np.isfinite(ypos)
692 X = np.asarray(xpos)[xidxs]
693 Y = np.asarray(ypos)[yidxs]
694 distance = [x[1].asArcseconds()
for x
in residuals]
696 distance = np.asarray(distance)[xidxs]
699 bearing = [x[0].asRadians()
for x
in residuals]
701 bearing = np.asarray(bearing)[xidxs]
702 U = (distance*np.cos(bearing))
703 V = (distance*np.sin(bearing))
704 sp.quiver(X, Y, U, V)
705 sp.set_title(
"WCS Residual")
710 """Utility class for dipole measurement testing.
712 Generate an image with simulated dipoles and noise; store the original
713 "pre-subtraction" images and catalogs as well.
714 Used to generate test data for DMTN-007 (http://dmtn-007.lsst.io).
717 def __init__(self, w=101, h=101, xcenPos=[27.], ycenPos=[25.], xcenNeg=[23.], ycenNeg=[25.],
718 psfSigma=2., flux=[30000.], fluxNeg=None, noise=10., gradientParams=None):
735 """Generate an exposure and catalog with the given dipole source(s).
744 dipole = posImage.clone()
745 di = dipole.getMaskedImage()
746 di -= negImage.getMaskedImage()
749 = dipole, posImage, posCatalog, negImage, negCatalog
751 def _makeStarImage(self, xc=[15.3], yc=[18.6], flux=[2500], schema=None, randomSeed=None):
752 """Generate an exposure and catalog with the given stellar source(s).
756 dataset = TestDataset(bbox, psfSigma=self.
psfSigma, threshold=1.)
758 for i
in range(len(xc)):
759 dataset.addSource(instFlux=flux[i], centroid=
geom.Point2D(xc[i], yc[i]))
762 schema = TestDataset.makeMinimalSchema()
763 exposure, catalog = dataset.realize(noise=self.
noise, schema=schema, randomSeed=randomSeed)
766 y, x = np.mgrid[:self.
w, :self.
h]
768 gradient = gp[0] + gp[1]*x + gp[2]*y
770 gradient += gp[3]*x*y + gp[4]*x*x + gp[5]*y*y
771 imgArr = exposure.image.array
774 return exposure, catalog
778 fitResult = alg.fitDipole(source, **kwds)
782 """Utility function for detecting dipoles.
784 Detect pos/neg sources in the diffim, then merge them. A
785 bigger "grow" parameter leads to a larger footprint which
786 helps with dipole measurement for faint dipoles.
791 Whether to merge the positive and negagive detections into a single
793 diffim : `lsst.afw.image.exposure.exposure.ExposureF`
794 Difference image on which to perform detection.
795 detectSigma : `float`
796 Threshold for object detection.
798 Number of pixels to grow the footprints before merging.
800 Minimum bin size for the background (re)estimation (only applies if
801 the default leads to min(nBinX, nBinY) < fit order so the default
802 config parameter needs to be decreased, but not to a value smaller
803 than ``minBinSize``, in which case the fitting algorithm will take
804 over and decrease the fit order appropriately.)
808 sources : `lsst.afw.table.SourceCatalog`
809 If doMerge=True, the merged source catalog is returned OR
810 detectTask : `lsst.meas.algorithms.SourceDetectionTask`
811 schema : `lsst.afw.table.Schema`
812 If doMerge=False, the source detection task and its schema are
819 schema = afwTable.SourceTable.makeMinimalSchema()
822 detectConfig = measAlg.SourceDetectionConfig()
823 detectConfig.returnOriginalFootprints =
False
826 detectConfig.thresholdPolarity =
"both"
827 detectConfig.thresholdValue = detectSigma
829 detectConfig.reEstimateBackground =
True
830 detectConfig.thresholdType =
"pixel_stdev"
831 detectConfig.excludeMaskPlanes = [
"EDGE"]
833 while ((
min(diffim.getWidth(), diffim.getHeight()))//detectConfig.background.binSize
834 < detectConfig.background.approxOrderX
and detectConfig.background.binSize > minBinSize):
835 detectConfig.background.binSize =
max(minBinSize, detectConfig.background.binSize//2)
838 detectTask = measAlg.SourceDetectionTask(schema, config=detectConfig)
840 table = afwTable.SourceTable.make(schema)
841 catalog = detectTask.run(table, diffim)
845 fpSet = catalog.positive
846 fpSet.merge(catalog.negative, grow, grow,
False)
848 fpSet.makeSources(sources)
853 return detectTask, schema
857 """Calculate the PSF matching kernel peak offset from the nominal
862 kernel : `~lsst.afw.math.LinearCombinationKernel`
863 The PSF matching kernel to evaluate.
865 The x position on the detector to evaluate the kernel
867 The y position on the detector to evaluate the kernel
868 image : `~lsst.afw.image.ImageD`
869 The image to use as base for computing kernel pixel values
874 The sum of the kernel on the desired location
876 The displacement of the kernel averaged peak, with respect to the
877 center of the extraction of the kernel
879 The displacement of the kernel averaged peak, with respect to the
880 center of the extraction of the kernel
882 The position angle in detector coordinates of the displacement
884 The displacement module of the kernel centroid in pixel units
888 image = afwImage.ImageD(kernel.getDimensions())
891 hsize = kernel.getWidth()//2
892 kernel_sum = kernel.computeImage(image, doNormalize=
False, x=x, y=y)
900 vx = data.sum(axis=0)
902 dx = np.dot(vx, xx) - hsize
904 vy = data.sum(axis=1)
906 dy = np.dot(vy, yy) - hsize
909 pos_angle = np.arctan2(dy, dx)
910 length = np.sqrt(dx**2 + dy**2)
912 return kernel_sum, dx, dy, pos_angle, length
916 """Directly calculate the horizontal and vertical widths
917 of a PSF at half its maximum value.
921 psf : `~lsst.afw.detection.Psf`
922 Point spread function (PSF) to evaluate.
923 average : `bool`, optional
924 Set to return the average width over Y and X axes.
925 position : `~lsst.geom.Point2D`, optional
926 The position at which to evaluate the PSF. If `None`, then the
927 average position is used.
931 psfSize : `float` | `tuple` [`float`]
932 The FWHM of the PSF computed at its average position.
933 Returns the widths along the Y and X axes,
934 or the average of the two if `average` is set.
941 position = psf.getAveragePosition()
942 shape = psf.computeShape(position)
943 sigmaToFwhm = 2*np.log(2*np.sqrt(2))
946 return sigmaToFwhm*shape.getTraceRadius()
948 return [sigmaToFwhm*np.sqrt(shape.getIxx()), sigmaToFwhm*np.sqrt(shape.getIyy())]
952 fwhmExposureBuffer: float, fwhmExposureGrid: int) -> float:
953 """Get the mean PSF FWHM by evaluating it on a grid within an exposure.
957 exposure : `~lsst.afw.image.Exposure`
958 The exposure for which the mean FWHM of the PSF is to be computed.
959 The exposure must contain a `psf` attribute.
960 fwhmExposureBuffer : `float`
961 Fractional buffer margin to be left out of all sides of the image
962 during the construction of the grid to compute mean PSF FWHM in an
964 fwhmExposureGrid : `int`
965 Grid size to compute the mean FWHM in an exposure.
970 The mean PSF FWHM on the exposure.
975 Raised if the PSF cannot be computed at any of the grid points.
985 bbox = exposure.getBBox()
986 xmax, ymax = bbox.getMax()
987 xmin, ymin = bbox.getMin()
989 xbuffer = fwhmExposureBuffer*(xmax-xmin)
990 ybuffer = fwhmExposureBuffer*(ymax-ymin)
993 for (x, y)
in itertools.product(np.linspace(xmin+xbuffer, xmax-xbuffer, fwhmExposureGrid),
994 np.linspace(ymin+ybuffer, ymax-ybuffer, fwhmExposureGrid)
998 fwhm = getPsfFwhm(psf, average=
True, position=pos)
999 except (InvalidParameterError, RangeError):
1000 _LOG.debug(
"Unable to compute PSF FWHM at position (%f, %f).", x, y)
1006 raise ValueError(
"Unable to compute PSF FWHM at any position on the exposure.")
1008 return np.nanmean(width)
1012 psfExposureBuffer: float, psfExposureGrid: int) -> afwImage.ImageD:
1013 """Get the average PSF by evaluating it on a grid within an exposure.
1017 exposure : `~lsst.afw.image.Exposure`
1018 The exposure for which the average PSF is to be computed.
1019 The exposure must contain a `psf` attribute.
1020 psfExposureBuffer : `float`
1021 Fractional buffer margin to be left out of all sides of the image
1022 during the construction of the grid to compute average PSF in an
1024 psfExposureGrid : `int`
1025 Grid size to compute the average PSF in an exposure.
1029 psfImage : `~lsst.afw.image.Image`
1030 The average PSF across the exposure.
1035 Raised if the PSF cannot be computed at any of the grid points.
1039 `evaluateMeanPsfFwhm`
1044 bbox = exposure.getBBox()
1045 xmax, ymax = bbox.getMax()
1046 xmin, ymin = bbox.getMin()
1048 xbuffer = psfExposureBuffer*(xmax-xmin)
1049 ybuffer = psfExposureBuffer*(ymax-ymin)
1053 for (x, y)
in itertools.product(np.linspace(xmin+xbuffer, xmax-xbuffer, psfExposureGrid),
1054 np.linspace(ymin+ybuffer, ymax-ybuffer, psfExposureGrid)
1058 singleImage = psf.computeKernelImage(pos)
1059 except InvalidParameterError:
1060 _LOG.debug(
"Unable to compute PSF image at position (%f, %f).", x, y)
1063 if psfArray
is None:
1064 psfArray = singleImage.array
1066 psfArray += singleImage.array
1069 if psfArray
is None:
1070 raise ValueError(
"Unable to compute PSF image at any position on the exposure.")
1072 psfImage = afwImage.ImageD(psfArray/nImg)
1077 """Minimal source detection wrapper suitable for unit tests.
1081 exposure : `lsst.afw.image.Exposure`
1082 Exposure on which to run detection/measurement
1083 The exposure is modified in place to set the 'DETECTED' mask plane.
1084 addMaskPlanes : `list` of `str`, optional
1085 Additional mask planes to add to the maskedImage of the exposure.
1090 Source catalog containing candidates
1092 if addMaskPlanes
is None:
1095 addMaskPlanes = [
"STREAK",
"INJECTED",
"INJECTED_TEMPLATE"]
1097 schema = afwTable.SourceTable.makeMinimalSchema()
1098 selectDetection = measAlg.SourceDetectionTask(schema=schema)
1099 selectMeasurement = measBase.SingleFrameMeasurementTask(schema=schema)
1100 table = afwTable.SourceTable.make(schema)
1102 detRet = selectDetection.run(
1108 for mp
in addMaskPlanes:
1109 exposure.mask.addMaskPlane(mp)
1111 selectSources = detRet.sources
1112 selectMeasurement.run(measCat=selectSources, exposure=exposure)
1114 return selectSources
1118 """Make a fake, affine Wcs.
1122 cdMatrix = np.array([[5.19513851e-05, -2.81124812e-07],
1123 [-3.25186974e-07, -5.19112119e-05]])
1128 noiseSeed=6, fluxLevel=500., fluxRange=2.,
1129 kernelSize=32, templateBorderSize=0,
1136 doApplyCalibration=False,
1140 clearEdgeMask=False,
1143 physicalFilter="g NotACamera"
1145 """Make a reproduceable PSF-convolved exposure for testing.
1149 seed : `int`, optional
1150 Seed value to initialize the random number generator for sources.
1151 nSrc : `int`, optional
1152 Number of sources to simulate.
1153 psfSize : `float`, optional
1154 Width of the PSF of the simulated sources, in pixels.
1155 noiseLevel : `float`, optional
1156 Standard deviation of the noise to add to each pixel.
1157 noiseSeed : `int`, optional
1158 Seed value to initialize the random number generator for noise.
1159 fluxLevel : `float`, optional
1160 Reference flux of the simulated sources.
1161 fluxRange : `float`, optional
1162 Range in flux amplitude of the simulated sources.
1163 kernelSize : `int`, optional
1164 Size in pixels of the kernel for simulating sources.
1165 templateBorderSize : `int`, optional
1166 Size in pixels of the image border used to pad the image.
1167 background : `lsst.afw.math.Chebyshev1Function2D`, optional
1168 Optional background to add to the output image.
1169 xSize, ySize : `int`, optional
1170 Size in pixels of the simulated image.
1171 x0, y0 : `int`, optional
1172 Origin of the image.
1173 calibration : `float`, optional
1174 Conversion factor between instFlux and nJy.
1175 doApplyCalibration : `bool`, optional
1176 Apply the photometric calibration and return the image in nJy?
1177 xLoc, yLoc : `list` of `float`, optional
1178 User-specified coordinates of the simulated sources.
1179 If specified, must have length equal to ``nSrc``
1180 flux : `list` of `float`, optional
1181 User-specified fluxes of the simulated sources.
1182 If specified, must have length equal to ``nSrc``
1183 clearEdgeMask : `bool`, optional
1184 Clear the "EDGE" mask plane after source detection.
1185 addMaskPlanes : `list` of `str`, optional
1186 Mask plane names to add to the image.
1190 modelExposure : `lsst.afw.image.Exposure`
1191 The model image, with the mask and variance planes. The DETECTED
1192 plane is filled in for the injected source footprints.
1193 sourceCat : `lsst.afw.table.SourceCatalog`
1194 Catalog of sources inserted in the model image.
1199 If `xloc`, `yloc`, or `flux` are supplied with inconsistant lengths.
1203 bufferSize = kernelSize/2 + templateBorderSize + 1
1206 if templateBorderSize > 0:
1207 bbox.grow(templateBorderSize)
1209 rng = np.random.RandomState(seed)
1210 rngNoise = np.random.RandomState(noiseSeed)
1211 x0, y0 = bbox.getBegin()
1212 xSize, ySize = bbox.getDimensions()
1214 xLoc = rng.rand(nSrc)*(xSize - 2*bufferSize) + bufferSize + x0
1216 if len(xLoc) != nSrc:
1217 raise ValueError(
"xLoc must have length equal to nSrc. %f supplied vs %f", len(xLoc), nSrc)
1219 yLoc = rng.rand(nSrc)*(ySize - 2*bufferSize) + bufferSize + y0
1221 if len(yLoc) != nSrc:
1222 raise ValueError(
"yLoc must have length equal to nSrc. %f supplied vs %f", len(yLoc), nSrc)
1225 flux = (rng.rand(nSrc)*(fluxRange - 1.) + 1.)*fluxLevel
1227 if len(flux) != nSrc:
1228 raise ValueError(
"flux must have length equal to nSrc. %f supplied vs %f", len(flux), nSrc)
1229 sigmas = [psfSize
for src
in range(nSrc)]
1230 injectList = list(zip(xLoc, yLoc, flux, sigmas))
1233 modelExposure = plantSources(bbox, kernelSize, skyLevel, injectList, addPoissonNoise=
False)
1235 noise = rngNoise.randn(ySize, xSize)*noiseLevel
1236 noise -= np.mean(noise)
1237 modelExposure.variance.array = np.sqrt(np.abs(modelExposure.image.array)) + noiseLevel**2
1238 modelExposure.image.array += noise
1243 modelExposure.mask &= ~modelExposure.mask.getPlaneBitMask(
"EDGE")
1245 if background
is not None:
1246 modelExposure.image += background
1247 modelExposure.maskedImage /= calibration
1249 modelExposure.info.setId(seed)
1250 if doApplyCalibration:
1251 modelExposure.maskedImage = modelExposure.photoCalib.calibrateImage(modelExposure.maskedImage)
1255 return modelExposure, truth
1259 """Make a schema for the truth catalog produced by `makeTestImage`.
1263 keys : `dict` [`str`]
1264 Fields added to the catalog, to make it easier to set them.
1265 schema : `lsst.afw.table.Schema`
1266 Schema to use to make a "truth" SourceCatalog.
1267 Calib, Ap, and Psf flux slots all are set to ``truth_instFlux``.
1269 schema = afwTable.SourceTable.makeMinimalSchema()
1272 keys[
"instFlux"] = schema.addField(
"truth_instFlux", type=np.float64,
1273 doc=
"true instFlux", units=
"count")
1274 keys[
"instFluxErr"] = schema.addField(
"truth_instFluxErr", type=np.float64,
1275 doc=
"true instFluxErr", units=
"count")
1276 keys[
"centroid"] = afwTable.Point2DKey.addFields(schema,
"truth",
"true simulated centroid",
"pixel")
1277 schema.addField(
"truth_flag",
"Flag",
"truth flux failure flag.")
1279 schema.addField(
"sky_source",
"Flag",
"testing flag.")
1280 schema.addField(
"base_PixelFlags_flag_interpolated",
"Flag",
"testing flag.")
1281 schema.addField(
"base_PixelFlags_flag_saturated",
"Flag",
"testing flag.")
1282 schema.addField(
"base_PixelFlags_flag_bad",
"Flag",
"testing flag.")
1283 schema.getAliasMap().set(
"slot_Centroid",
"truth")
1284 schema.getAliasMap().set(
"slot_CalibFlux",
"truth")
1285 schema.getAliasMap().set(
"slot_ApFlux",
"truth")
1286 schema.getAliasMap().set(
"slot_PsfFlux",
"truth")
1291 """Add injected sources to the truth catalog.
1295 injectList : `list` [`float`]
1296 Sources that were injected; tuples of (x, y, flux, size).
1300 catalog : `lsst.afw.table.SourceCatalog`
1301 Catalog with centroids and instFlux/instFluxErr values filled in and
1302 appropriate slots set.
1306 catalog.reserve(len(injectList))
1307 for x, y, flux, size
in injectList:
1308 record = catalog.addNew()
1310 keys[
"instFlux"].set(record, flux)
1312 keys[
"instFluxErr"].set(record, 20)
1315 footprint = afwDetection.Footprint(afwGeom.SpanSet.fromShape(circle))
1316 footprint.addPeak(x, y, flux)
1317 record.setFootprint(footprint)
1323 """Create a statistics control for configuring calculations on images.
1327 badMaskPlanes : `list` of `str`, optional
1328 List of mask planes to exclude from calculations.
1332 statsControl : ` lsst.afw.math.StatisticsControl`
1333 Statistics control object for configuring calculations on images.
1335 if badMaskPlanes
is None:
1336 badMaskPlanes = (
"INTRP",
"EDGE",
"DETECTED",
"SAT",
"CR",
1337 "BAD",
"NO_DATA",
"DETECTED_NEGATIVE")
1339 statsControl.setNumSigmaClip(3.)
1340 statsControl.setNumIter(3)
1341 statsControl.setAndMask(afwImage.Mask.getPlaneBitMask(badMaskPlanes))
1346 """Calculate a robust mean of the variance plane of an exposure.
1350 image : `lsst.afw.image.Image`
1351 Image or variance plane of an exposure to evaluate.
1352 mask : `lsst.afw.image.Mask`
1353 Mask plane to use for excluding pixels.
1354 statsCtrl : `lsst.afw.math.StatisticsControl`
1355 Statistics control object for configuring the calculation.
1356 statistic : `lsst.afw.math.Property`, optional
1357 The type of statistic to compute. Typical values are
1358 ``afwMath.MEANCLIP`` or ``afwMath.STDEVCLIP``.
1363 The result of the statistic calculated from the unflagged pixels.
1366 return statObj.getValue(statistic)
1370 """Compute the noise equivalent area for an image psf
1374 psf : `lsst.afw.detection.Psf`
1380 psfImg = psf.computeImage(psf.getAveragePosition())
1381 nea = 1./np.sum(psfImg.array**2)
1386 """Calculate the mean of an array of angles.
1391 An array of angles, in radians
1398 complexArray = [complex(np.cos(angle), np.sin(angle))
for angle
in angles]
1399 return (
geom.Angle(np.angle(np.mean(complexArray))))
1403 """Evaluate the fraction of masked pixels in a mask plane.
1407 mask : `lsst.afw.image.Mask`
1408 The mask to evaluate the fraction on
1410 The particular mask plane to evaluate the fraction
1415 The calculated fraction of masked pixels
1417 nMaskSet = np.count_nonzero((mask.array & mask.getPlaneBitMask(maskPlane)))
1418 return nMaskSet/mask.array.size
1422 """A custom CoaddPSF that overrides the getAveragePosition method.
1424 It intentionally moves the position off-image to cause a test failure.
A compact representation of a collection of pixels.
An ellipse defined by an arbitrary BaseCore and a center point.
A group of labels for a filter in an exposure or coadd.
The photometric calibration of an exposure.
A kernel created from an Image.
Pass parameters to a Statistics object.
A class representing an angle.
A floating-point coordinate rectangle geometry.
An integer coordinate rectangle.
Point in an unspecified spherical coordinate system.
detectDipoleSources(self, doMerge=True, diffim=None, detectSigma=5.5, grow=3, minBinSize=32)
_makeStarImage(self, xc=[15.3], yc=[18.6], flux=[2500], schema=None, randomSeed=None)
__init__(self, w=101, h=101, xcenPos=[27.], ycenPos=[25.], xcenNeg=[23.], ycenNeg=[25.], psfSigma=2., flux=[30000.], fluxNeg=None, noise=10., gradientParams=None)
fitDipoleSource(self, source, **kwds)
std::shared_ptr< SkyWcs > makeSkyWcs(daf::base::PropertySet &metadata, bool strip=false)
Construct a SkyWcs from FITS keywords.
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.
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.
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)
afwImage.ImageD computeAveragePsf(afwImage.Exposure exposure, float psfExposureBuffer, int psfExposureGrid)
plotKernelSpatialModel(kernel, kernelCellSet, showBadCandidates=True, numSample=128, keepPlots=True, maxCoeff=10)
detectTestSources(exposure, addMaskPlanes=None)
getKernelCenterDisplacement(kernel, x, y, image=None)
makeStats(badMaskPlanes=None)
showKernelMosaic(bbox, kernel, nx=7, ny=None, frame=None, title=None, showCenter=True, showEllipticity=True)
showSourceSet(sSet, xy0=(0, 0), frame=0, ctype=afwDisplay.GREEN, symb="+", size=2)
showKernelSpatialCells(maskedIm, kernelCellSet, showChi2=False, symb="o", ctype=None, ctypeUnused=None, ctypeBad=None, size=3, frame=None, title="Spatial Cells")
plotWhisker(results, newWcs)
showKernelBasis(kernel, frame=None)
_fillTruthCatalog(injectList)
showDiaSources(sources, exposure, isFlagged, isDipole, frame=None)
computePSFNoiseEquivalentArea(psf)
showKernelCandidates(kernelCellSet, kernel, background, frame=None, showBadCandidates=True, resids=False, kernels=False)
evaluateMaskFraction(mask, maskPlane)
plotKernelCoefficients(spatialKernel, kernelCellSet, showBadCandidates=False, keepPlots=True)
makeTestImage(seed=5, nSrc=20, psfSize=2., noiseLevel=5., noiseSeed=6, fluxLevel=500., fluxRange=2., kernelSize=32, templateBorderSize=0, background=None, xSize=256, ySize=256, x0=12345, y0=67890, calibration=1., doApplyCalibration=False, xLoc=None, yLoc=None, flux=None, clearEdgeMask=False, addMaskPlanes=None, band="g", physicalFilter="g NotACamera")
getPsfFwhm(psf, average=True, position=None)
float evaluateMeanPsfFwhm(afwImage.Exposure exposure, float fwhmExposureBuffer, int fwhmExposureGrid)
computeRobustStatistics(image, mask, statsCtrl, statistic=afwMath.MEANCLIP)