Coverage for tests / utils.py: 5%
681 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-05-07 08:35 +0000
1#
2# Developed for the LSST Data Management System.
3# This product includes software developed by the LSST Project
4# (https://www.lsst.org).
5# See the COPYRIGHT file at the top-level directory of this distribution
6# for details of code ownership.
7#
8# This program is free software: you can redistribute it and/or modify
9# it under the terms of the GNU General Public License as published by
10# the Free Software Foundation, either version 3 of the License, or
11# (at your option) any later version.
12#
13# This program is distributed in the hope that it will be useful,
14# but WITHOUT ANY WARRANTY; without even the implied warranty of
15# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16# GNU General Public License for more details.
17#
18# You should have received a copy of the GNU General Public License
19# along with this program. If not, see <https://www.gnu.org/licenses/>.
21"""Support utilities for Measuring sources"""
23# Export DipoleTestImage to expose fake image generating funcs
24__all__ = ["DipoleTestImage"]
27import numpy as np
28import lsst.geom as geom
29import lsst.afw.detection as afwDet
30import lsst.afw.display as afwDisplay
31import lsst.afw.detection as afwDetection
32import lsst.afw.geom as afwGeom
33import lsst.afw.image as afwImage
34import lsst.afw.math as afwMath
35import lsst.afw.table as afwTable
36from lsst.daf.butler import DataCoordinate, DimensionUniverse
37import lsst.meas.algorithms as measAlg
38import lsst.meas.base as measBase
39from lsst.meas.algorithms.testUtils import plantSources
40from lsst.ip.diffim.dipoleFitTask import DipoleFitAlgorithm
41from lsst.ip.diffim import diffimLib
43afwDisplay.setDefaultMaskTransparency(75)
44keptPlots = False # Have we arranged to keep spatial plots open?
47def 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.
51 """
52 disp = afwDisplay.afwDisplay(frame=frame)
53 with disp.Buffering():
54 for s in sSet:
55 xc, yc = s.getXAstrom() - xy0[0], s.getYAstrom() - xy0[1]
57 if symb == "id":
58 disp.dot(str(s.getId()), xc, yc, ctype=ctype, size=size)
59 else:
60 disp.dot(symb, xc, yc, ctype=ctype, size=size)
63# Kernel display utilities
64#
67def showKernelSpatialCells(maskedIm, kernelCellSet, showChi2=False, symb="o",
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
74 and size.
75 """
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:
87 color = ctypeBad
88 elif cand.getStatus() == afwMath.SpatialCellCandidate.GOOD:
89 color = ctype
90 elif cand.getStatus() == afwMath.SpatialCellCandidate.UNKNOWN:
91 color = ctypeUnused
92 else:
93 continue
95 if color:
96 disp.dot(symb, xc, yc, ctype=color, size=size)
98 if showChi2:
99 rchi2 = cand.getChi2()
100 if rchi2 > 1e100:
101 rchi2 = np.nan
102 disp.dot("%d %.1f" % (cand.getId(), rchi2),
103 xc - size, yc - size - 4, ctype=color, size=size)
106def showDiaSources(sources, exposure, isFlagged, isDipole, frame=None):
107 """Display Dia Sources.
108 """
109 #
110 # Show us the ccandidates
111 #
112 # Too many mask planes in diffims
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)):
119 source = sources[i]
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())
127 if badFlag:
128 ctype = afwDisplay.RED
129 lab += "BAD"
130 if dipoleFlag:
131 ctype = afwDisplay.YELLOW
132 lab += "DIPOLE"
133 if not badFlag and not dipoleFlag:
134 ctype = afwDisplay.GREEN
135 lab += "OK"
136 mos.append(im.makeMosaic(), lab, ctype)
137 title = "Dia Sources"
138 mosaicImage = mos.makeMosaic(display=disp, title=title)
139 return mosaicImage
142def showKernelCandidates(kernelCellSet, kernel, background, frame=None, showBadCandidates=True,
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.
148 """
149 #
150 # Show us the ccandidates
151 #
152 if kernels:
153 mos = afwDisplay.utils.Mosaic(gutter=5, background=0)
154 else:
155 mos = afwDisplay.utils.Mosaic(gutter=5, background=-1)
156 #
157 candidateCenters = []
158 candidateCentersBad = []
159 candidateIndex = 0
160 for cell in kernelCellSet.getCellList():
161 for cand in cell.begin(False): # include bad candidates
162 # Original difference image; if does not exist, skip candidate
163 try:
164 resid = cand.getDifferenceImage(diffimLib.KernelCandidateF.ORIG)
165 except Exception:
166 continue
168 rchi2 = cand.getChi2()
169 if rchi2 > 1e100:
170 rchi2 = np.nan
172 if not showBadCandidates and cand.isBad():
173 continue
175 im_resid = afwDisplay.utils.Mosaic(gutter=1, background=-0.5, mode="x")
177 try:
178 im = cand.getScienceMaskedImage()
179 im = im.Factory(im, True)
180 im.setXY0(cand.getScienceMaskedImage().getXY0())
181 except Exception:
182 continue
183 if (not resids and not kernels):
184 im_resid.append(im.Factory(im, True))
185 try:
186 im = cand.getTemplateMaskedImage()
187 im = im.Factory(im, True)
188 im.setXY0(cand.getTemplateMaskedImage().getXY0())
189 except Exception:
190 continue
191 if (not resids and not kernels):
192 im_resid.append(im.Factory(im, True))
194 # Difference image with original basis
195 if resids:
196 var = resid.variance
197 var = var.Factory(var, True)
198 np.sqrt(var.array, var.array) # inplace sqrt
199 resid = resid.image
200 resid /= var
201 bbox = kernel.shrinkBBox(resid.getBBox())
202 resid = resid.Factory(resid, bbox, deep=True)
203 elif kernels:
204 kim = cand.getKernelImage(diffimLib.KernelCandidateF.ORIG).convertF()
205 resid = kim.Factory(kim, True)
206 im_resid.append(resid)
208 # residuals using spatial model
209 ski = afwImage.ImageD(kernel.getDimensions())
210 kernel.computeImage(ski, False, int(cand.getXCenter()), int(cand.getYCenter()))
211 sk = afwMath.FixedKernel(ski)
212 sbg = 0.0
213 if background:
214 sbg = background(int(cand.getXCenter()), int(cand.getYCenter()))
215 sresid = cand.getDifferenceImage(sk, sbg)
216 resid = sresid
217 if resids:
218 resid = sresid.image
219 resid /= var
220 bbox = kernel.shrinkBBox(resid.getBBox())
221 resid = resid.Factory(resid, bbox, deep=True)
222 elif kernels:
223 kim = ski.convertF()
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.image, title="candidate")
237 print("rating", cand.getCandidateRating())
239 im = cand.getScienceMaskedImage()
240 center = (candidateIndex, cand.getXCenter() - im.getX0(), cand.getYCenter() - im.getY0())
241 candidateIndex += 1
242 if cand.isBad():
243 candidateCentersBad.append(center)
244 else:
245 candidateCenters.append(center)
247 if resids:
248 title = "chi Diffim"
249 elif kernels:
250 title = "Kernels"
251 else:
252 title = "Candidates & residuals"
254 disp = afwDisplay.Display(frame=frame)
255 mosaicImage = mos.makeMosaic(display=disp, title=title)
257 return mosaicImage
260def showKernelBasis(kernel, frame=None):
261 """Display a Kernel's basis images.
262 """
263 mos = afwDisplay.utils.Mosaic()
265 for k in kernel.getKernelList():
266 im = afwImage.ImageD(k.getDimensions())
267 k.computeImage(im, False)
268 mos.append(im)
270 disp = afwDisplay.Display(frame=frame)
271 mos.makeMosaic(display=disp, title="Kernel Basis Images")
273 return mos
275###############
278def plotKernelSpatialModel(kernel, kernelCellSet, showBadCandidates=True,
279 numSample=128, keepPlots=True, maxCoeff=10):
280 """Plot the Kernel spatial model.
281 """
282 try:
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)
287 return
289 x0 = kernelCellSet.getBBox().getBeginX()
290 y0 = kernelCellSet.getBBox().getBeginY()
292 candPos = list()
293 candFits = list()
294 badPos = list()
295 badFits = list()
296 candAmps = list()
297 badAmps = list()
298 for cell in kernelCellSet.getCellList():
299 for cand in cell.begin(False):
300 if not showBadCandidates and cand.isBad():
301 continue
302 candCenter = geom.PointD(cand.getXCenter(), cand.getYCenter())
303 try:
304 im = cand.getTemplateMaskedImage()
305 except Exception:
306 continue
308 targetFits = badFits if cand.isBad() else candFits
309 targetPos = badPos if cand.isBad() else candPos
310 targetAmps = badAmps if cand.isBad() else candAmps
312 # compare original and spatial kernel coefficients
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)
331 numBad = len(badPos)
333 xRange = np.linspace(0, kernelCellSet.getBBox().getWidth(), num=numSample)
334 yRange = np.linspace(0, kernelCellSet.getBBox().getHeight(), num=numSample)
336 if maxCoeff:
337 maxCoeff = min(maxCoeff, kernel.getNKernelParameters())
338 else:
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])
344 yMin = dfGood.min()
345 yMax = dfGood.max()
346 if numBad > 0:
347 dfBad = zBad[:, k] - np.array([func(pos.getX(), pos.getY()) for pos in badPos])
348 # Can really screw up the range...
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)
359 fig = plt.figure(k)
361 fig.clf()
362 try:
363 fig.canvas._tkcanvas._root().lift() # == Tk's raise, but raise is a python reserved word
364 except Exception: # protect against API changes
365 pass
367 fig.suptitle('Kernel component %d' % k)
369 # LL
370 ax = fig.add_axes((0.1, 0.05, 0.35, 0.35))
371 vmin = fRange.min() # - 0.05*np.fabs(fRange.min())
372 vmax = fRange.max() # + 0.05*np.fabs(fRange.max())
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])
380 # UL
381 ax = fig.add_axes((0.1, 0.55, 0.35, 0.35))
382 ax.plot(-2.5*np.log10(candAmps), zGood[:, k], 'b+')
383 if numBad > 0:
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")
389 # LR
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+')
395 if numBad > 0:
396 ax.plot(yBad, dfBad, 'r+')
397 ax.axhline(0.0)
398 ax.set_title('dCoeff (indiv-spatial) vs. y')
400 # UR
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+')
406 if numBad > 0:
407 ax.plot(xBad, dfBad, 'r+')
408 ax.axhline(0.0)
409 ax.set_title('dCoeff (indiv-spatial) vs. x')
411 fig.show()
413 global keptPlots
414 if keepPlots and not keptPlots:
415 # Keep plots open when done
416 def show():
417 print("%s: Please close plots when done." % __name__)
418 try:
419 plt.show()
420 except Exception:
421 pass
422 print("Plots closed, exiting...")
423 import atexit
424 atexit.register(show)
425 keptPlots = True
428def plotKernelCoefficients(spatialKernel, kernelCellSet, showBadCandidates=False, keepPlots=True):
429 """Plot the individual kernel candidate and the spatial kernel solution coefficients.
431 Parameters
432 ----------
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.
451 Notes
452 -----
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
460 center of the image.
462 This function is called if ``lsst.ip.diffim.psfMatch.plotKernelCoefficients==True`` in lsstDebug. This
463 function was implemented as part of DM-17825.
464 """
465 try:
466 import matplotlib.pyplot as plt
467 except ImportError as e:
468 print("Unable to import matplotlib: %s" % e)
469 return
471 # Image dimensions
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()
480 # Plot the local solutions
481 # ----
483 # Grid size
484 nX = 8
485 nY = 8
486 wCell = wImage / nX
487 hCell = hImage / nY
489 fig = plt.figure()
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(
492 wspace=0, hspace=0))
494 # Bottom left panel is for bottom left part of the image
495 arrAx = arrAx[::-1, :]
497 allParams = []
498 for cell in kernelCellSet.getCellList():
499 cellBBox = geom.Box2D(cell.getBBox())
500 # Determine which panel this spatial cell belongs to
501 iX = int((cellBBox.getCenterX() - x0)//wCell)
502 iY = int((cellBBox.getCenterY() - y0)//hCell)
504 for cand in cell.begin(False):
505 try:
506 kernel = cand.getKernel(cand.ORIG)
507 except Exception:
508 continue
510 if not showBadCandidates and cand.isBad():
511 continue
513 nKernelParams = kernel.getNKernelParameters()
514 kernelParams = np.array(kernel.getKernelParameters())
515 allParams.append(kernelParams)
517 if cand.isBad():
518 color = 'red'
519 else:
520 color = None
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')
526 # Plot histogram of the local parameters and the global solution at the image center
527 # ----
529 spatialFuncs = spatialKernel.getSpatialFunctionList()
530 nKernelParams = spatialKernel.getNKernelParameters()
531 nX = 8
532 fig = plt.figure()
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')
546 # Plot grid of the spatial solution
547 # ----
549 nX = 8
550 nY = 8
551 wCell = wImage / nX
552 hCell = hImage / nY
553 x0 += wCell / 2
554 y0 += hCell / 2
556 fig = plt.figure()
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(
559 wspace=0, hspace=0))
560 arrAx = arrAx[::-1, :]
561 kernelParams = np.zeros(nKernelParams, dtype=float)
563 for iX in range(nX):
564 for iY in range(nY):
565 x = x0 + iX * wCell
566 y = y0 + iY * hCell
567 # Evaluate the spatial solution functions for this x,y location
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')
572 global keptPlots
573 if keepPlots and not keptPlots:
574 # Keep plots open when done
575 def show():
576 print("%s: Please close plots when done." % __name__)
577 try:
578 plt.show()
579 except Exception:
580 pass
581 print("Plots closed, exiting...")
582 import atexit
583 atexit.register(show)
584 keptPlots = True
587def showKernelMosaic(bbox, kernel, nx=7, ny=None, frame=None, title=None,
588 showCenter=True, showEllipticity=True):
589 """Show a mosaic of Kernel images.
590 """
591 mos = afwDisplay.utils.Mosaic()
593 x0 = bbox.getBeginX()
594 y0 = bbox.getBeginY()
595 width = bbox.getWidth()
596 height = bbox.getHeight()
598 if not ny:
599 ny = int(nx*float(height)/width + 0.5)
600 if not ny:
601 ny = 1
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)
616 centers = []
617 shapes = []
618 for iy in range(ny):
619 for ix in range(nx):
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 ""
626 mos.append(im, lab)
628 # SdssCentroidAlgorithm.measure requires an exposure of floats
629 exp = afwImage.makeExposure(afwImage.makeMaskedImage(im.convertF()))
631 w, h = im.getWidth(), im.getHeight()
632 centerX = im.getX0() + w//2
633 centerY = im.getY0() + h//2
634 src = table.makeRecord()
635 spans = afwGeom.SpanSet(exp.getBBox())
636 foot = afwDet.Footprint(spans)
637 foot.addPeak(centerX, centerY, 1)
638 src.setFootprint(foot)
640 try: # The centroider requires a psf, so this will fail if none is attached to exp
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()))
646 except Exception:
647 pass
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)
654 i = 0
655 with disp.Buffering():
656 for cen, shape in zip(centers, shapes):
657 bbox = mos.getBBox(i)
658 i += 1
659 xc, yc = cen[0] + bbox.getMinX(), cen[1] + bbox.getMinY()
660 if showCenter:
661 disp.dot("+", xc, yc, ctype=afwDisplay.BLUE)
663 if showEllipticity:
664 ixx, ixy, iyy = shape
665 disp.dot("@:%g,%g,%g" % (ixx, ixy, iyy), xc, yc, ctype=afwDisplay.RED)
667 return mos
670def plotWhisker(results, newWcs):
671 """Plot whisker diagram of astromeric offsets between results.matches.
672 """
673 refCoordKey = results.matches[0].first.getTable().getCoordKey()
674 inCentroidKey = results.matches[0].second.getTable().getCentroidSlot().getMeasKey()
675 positions = [m.first.get(refCoordKey) for m in results.matches]
676 residuals = [m.first.get(refCoordKey).getOffsetFrom(
677 newWcs.pixelToSky(m.second.get(inCentroidKey))) for
678 m in results.matches]
679 import matplotlib.pyplot as plt
680 fig = plt.figure()
681 sp = fig.add_subplot(1, 1, 0)
682 xpos = [x[0].asDegrees() for x in positions]
683 ypos = [x[1].asDegrees() for x in positions]
684 xpos.append(0.02*(max(xpos) - min(xpos)) + min(xpos))
685 ypos.append(0.98*(max(ypos) - min(ypos)) + min(ypos))
686 xidxs = np.isfinite(xpos)
687 yidxs = np.isfinite(ypos)
688 X = np.asarray(xpos)[xidxs]
689 Y = np.asarray(ypos)[yidxs]
690 distance = [x[1].asArcseconds() for x in residuals]
691 distance.append(0.2)
692 distance = np.asarray(distance)[xidxs]
693 # NOTE: This assumes that the bearing is measured positive from +RA through North.
694 # From the documentation this is not clear.
695 bearing = [x[0].asRadians() for x in residuals]
696 bearing.append(0)
697 bearing = np.asarray(bearing)[xidxs]
698 U = (distance*np.cos(bearing))
699 V = (distance*np.sin(bearing))
700 sp.quiver(X, Y, U, V)
701 sp.set_title("WCS Residual")
702 plt.show()
705class DipoleTestImage:
706 """Utility class for dipole measurement testing.
708 Generate an image with simulated dipoles and noise; store the original
709 "pre-subtraction" images and catalogs as well.
710 Used to generate test data for DMTN-007 (http://dmtn-007.lsst.io).
711 """
713 def __init__(self, w=101, h=101, xcenPos=[27.], ycenPos=[25.], xcenNeg=[23.], ycenNeg=[25.],
714 psfSigma=2., flux=[30000.], fluxNeg=None, noise=10., gradientParams=None, edgeWidth=8):
715 self.w = w
716 self.h = h
717 self.xcenPos = xcenPos
718 self.ycenPos = ycenPos
719 self.xcenNeg = xcenNeg
720 self.ycenNeg = ycenNeg
721 self.psfSigma = psfSigma
722 self.flux = flux
723 self.fluxNeg = fluxNeg
724 if fluxNeg is None:
725 self.fluxNeg = self.flux
726 self.noise = noise
727 self.gradientParams = gradientParams
728 self.edgeWidth = edgeWidth
729 self._makeDipoleImage()
731 def _makeDipoleImage(self):
732 """Generate an exposure and catalog with the given dipole source(s).
733 """
734 # Must seed the pos/neg images with different values to ensure they get different noise realizations
735 posImage, posCatalog = self._makeStarImage(
736 xc=self.xcenPos, yc=self.ycenPos, flux=self.flux, randomSeed=111)
738 negImage, negCatalog = self._makeStarImage(
739 xc=self.xcenNeg, yc=self.ycenNeg, flux=self.fluxNeg, randomSeed=222)
741 dipole = posImage.clone()
742 di = dipole.getMaskedImage()
743 di -= negImage.getMaskedImage()
745 self.diffim, self.posImage, self.posCatalog, self.negImage, self.negCatalog \
746 = dipole, posImage, posCatalog, negImage, negCatalog
748 def _makeStarImage(self, xc=[15.3], yc=[18.6], flux=[2500], schema=None, randomSeed=None):
749 """Generate an exposure and catalog with the given stellar source(s).
750 """
751 from lsst.meas.base.tests import TestDataset
752 bbox = geom.Box2I(geom.Point2I(0, 0), geom.Point2I(self.w - 1, self.h - 1))
753 dataset = TestDataset(bbox, psfSigma=self.psfSigma, threshold=1.)
755 for i in range(len(xc)):
756 dataset.addSource(instFlux=flux[i], centroid=geom.Point2D(xc[i], yc[i]))
758 if schema is None:
759 schema = TestDataset.makeMinimalSchema()
760 exposure, catalog = dataset.realize(noise=self.noise, schema=schema, randomSeed=randomSeed)
762 # set EDGE by masking the whole exposure and un-masking an inner bbox
763 edgeMask = exposure.mask.getPlaneBitMask('EDGE')
764 exposure.mask.array |= edgeMask
765 inner_bbox = exposure.getBBox()
766 inner_bbox.grow(-self.edgeWidth)
767 exposure[inner_bbox].mask.array &= ~edgeMask
769 if self.gradientParams is not None:
770 y, x = np.mgrid[:self.w, :self.h]
771 gp = self.gradientParams
772 gradient = gp[0] + gp[1]*x + gp[2]*y
773 if len(self.gradientParams) > 3: # it includes a set of 2nd-order polynomial params
774 gradient += gp[3]*x*y + gp[4]*x*x + gp[5]*y*y
775 imgArr = exposure.image.array
776 imgArr += gradient
778 return exposure, catalog
780 def fitDipoleSource(self, source, **kwds):
781 alg = DipoleFitAlgorithm(self.diffim, self.posImage, self.negImage)
782 fitResult = alg.fitDipole(source, **kwds)
783 return fitResult
785 def detectDipoleSources(self, doMerge=True, diffim=None, detectSigma=5.5, grow=3, minBinSize=32):
786 """Utility function for detecting dipoles.
788 Detect pos/neg sources in the diffim, then merge them. A
789 bigger "grow" parameter leads to a larger footprint which
790 helps with dipole measurement for faint dipoles.
792 Parameters
793 ----------
794 doMerge : `bool`
795 Whether to merge the positive and negagive detections into a single
796 source table.
797 diffim : `lsst.afw.image.exposure.exposure.ExposureF`
798 Difference image on which to perform detection.
799 detectSigma : `float`
800 Threshold for object detection.
801 grow : `int`
802 Number of pixels to grow the footprints before merging.
803 minBinSize : `int`
804 Minimum bin size for the background (re)estimation (only applies if
805 the default leads to min(nBinX, nBinY) < fit order so the default
806 config parameter needs to be decreased, but not to a value smaller
807 than ``minBinSize``, in which case the fitting algorithm will take
808 over and decrease the fit order appropriately.)
810 Returns
811 -------
812 sources : `lsst.afw.table.SourceCatalog`
813 If doMerge=True, the merged source catalog is returned OR
814 detectTask : `lsst.meas.algorithms.SourceDetectionTask`
815 schema : `lsst.afw.table.Schema`
816 If doMerge=False, the source detection task and its schema are
817 returned.
818 """
819 if diffim is None:
820 diffim = self.diffim
822 # Start with a minimal schema - only the fields all SourceCatalogs need
823 schema = afwTable.SourceTable.makeMinimalSchema()
825 # Customize the detection task a bit (optional)
826 detectConfig = measAlg.SourceDetectionConfig()
827 detectConfig.returnOriginalFootprints = False # should be the default
829 # code from imageDifference.py:
830 detectConfig.thresholdPolarity = "both"
831 detectConfig.thresholdValue = detectSigma
832 # detectConfig.nSigmaToGrow = psfSigma
833 detectConfig.reEstimateBackground = True # if False, will fail often for faint sources on gradients?
834 detectConfig.thresholdType = "pixel_stdev"
835 detectConfig.excludeMaskPlanes = ["EDGE"]
836 # Test images are often quite small, so may need to adjust background binSize
837 while ((min(diffim.getWidth(), diffim.getHeight()))//detectConfig.background.binSize
838 < detectConfig.background.approxOrderX and detectConfig.background.binSize > minBinSize):
839 detectConfig.background.binSize = max(minBinSize, detectConfig.background.binSize//2)
841 # Create the detection task. We pass the schema so the task can declare a few flag fields
842 detectTask = measAlg.SourceDetectionTask(schema, config=detectConfig)
844 table = afwTable.SourceTable.make(schema)
845 catalog = detectTask.run(table, diffim)
847 # Now do the merge.
848 if doMerge:
849 fpSet = catalog.positive
850 fpSet.merge(catalog.negative, grow, grow, False)
851 sources = afwTable.SourceCatalog(table)
852 fpSet.makeSources(sources)
854 return sources
856 else:
857 return detectTask, schema
860def detectTestSources(exposure, addMaskPlanes=None):
861 """Minimal source detection wrapper suitable for unit tests.
863 Parameters
864 ----------
865 exposure : `lsst.afw.image.Exposure`
866 Exposure on which to run detection/measurement
867 The exposure is modified in place to set the 'DETECTED' mask plane.
868 addMaskPlanes : `list` of `str`, optional
869 Additional mask planes to add to the maskedImage of the exposure.
871 Returns
872 -------
873 selectSources
874 Source catalog containing candidates
875 """
876 if addMaskPlanes is None:
877 # And add empty INJECTED and INJECTED_TEMPLATE mask planes
878 addMaskPlanes = ["INJECTED", "INJECTED_TEMPLATE"]
880 schema = afwTable.SourceTable.makeMinimalSchema()
882 selectDetection = measAlg.SourceDetectionTask(schema=schema)
883 selectMeasurement = measBase.SingleFrameMeasurementTask(schema=schema)
884 table = afwTable.SourceTable.make(schema)
886 detRet = selectDetection.run(
887 table=table,
888 exposure=exposure,
889 sigma=None, # The appropriate sigma is calculated from the PSF
890 doSmooth=True
891 )
892 for mp in addMaskPlanes:
893 exposure.mask.addMaskPlane(mp)
895 selectSources = detRet.sources
896 selectMeasurement.run(measCat=selectSources, exposure=exposure)
898 return selectSources
901def makeFakeWcs():
902 """Make a fake, affine Wcs.
903 """
904 crpix = geom.Point2D(123.45, 678.9)
905 crval = geom.SpherePoint(0.1, 0.1, geom.degrees)
906 cdMatrix = np.array([[5.19513851e-05, -2.81124812e-07],
907 [-3.25186974e-07, -5.19112119e-05]])
908 return afwGeom.makeSkyWcs(crpix, crval, cdMatrix)
911def makeTestImage(seed=5, nSrc=20, psfSize=2., noiseLevel=5.,
912 noiseSeed=6, fluxLevel=500., fluxRange=2.,
913 kernelSize=31, templateBorderSize=0,
914 background=None,
915 xSize=256,
916 ySize=256,
917 x0=12345,
918 y0=67890,
919 calibration=1.,
920 doApplyCalibration=False,
921 xLoc=None,
922 yLoc=None,
923 flux=None,
924 clearEdgeMask=False,
925 addMaskPlanes=None,
926 band="g",
927 physicalFilter="g NotACamera"
928 ):
929 """Make a reproduceable PSF-convolved exposure for testing.
931 Parameters
932 ----------
933 seed : `int`, optional
934 Seed value to initialize the random number generator for sources.
935 nSrc : `int`, optional
936 Number of sources to simulate.
937 psfSize : `float`, optional
938 Width of the PSF of the simulated sources, in pixels.
939 noiseLevel : `float`, optional
940 Standard deviation of the noise to add to each pixel.
941 noiseSeed : `int`, optional
942 Seed value to initialize the random number generator for noise.
943 fluxLevel : `float`, optional
944 Reference flux of the simulated sources.
945 fluxRange : `float`, optional
946 Range in flux amplitude of the simulated sources.
947 kernelSize : `int`, optional
948 Size in pixels of the kernel for simulating sources. Must be odd.
949 templateBorderSize : `int`, optional
950 Size in pixels of the image border used to pad the image.
951 background : `lsst.afw.math.Chebyshev1Function2D`, optional
952 Optional background to add to the output image.
953 xSize, ySize : `int`, optional
954 Size in pixels of the simulated image.
955 x0, y0 : `int`, optional
956 Origin of the image.
957 calibration : `float`, optional
958 Conversion factor between instFlux and nJy.
959 doApplyCalibration : `bool`, optional
960 Apply the photometric calibration and return the image in nJy?
961 xLoc, yLoc : `list` of `float`, optional
962 User-specified coordinates of the simulated sources.
963 If specified, must have length equal to ``nSrc``
964 flux : `list` of `float`, optional
965 User-specified fluxes of the simulated sources.
966 If specified, must have length equal to ``nSrc``
967 clearEdgeMask : `bool`, optional
968 Clear the "EDGE" mask plane after source detection.
969 addMaskPlanes : `list` of `str`, optional
970 Mask plane names to add to the image.
972 Returns
973 -------
974 modelExposure : `lsst.afw.image.Exposure`
975 The model image, with the mask and variance planes. The DETECTED
976 plane is filled in for the injected source footprints.
977 sourceCat : `lsst.afw.table.SourceCatalog`
978 Catalog of sources inserted in the model image.
980 Raises
981 ------
982 ValueError
983 If `xloc`, `yloc`, or `flux` are supplied with inconsistant lengths,
984 or if ``kernelSize`` is even.
985 """
986 if kernelSize % 2 == 0:
987 raise ValueError(f"kernelSize must be odd, got {kernelSize}.")
988 # Distance from the inner edge of the bounding box to avoid placing test
989 # sources in the model images.
990 bufferSize = kernelSize//2 + templateBorderSize + 1
992 bbox = geom.Box2I(geom.Point2I(x0, y0), geom.Extent2I(xSize, ySize))
993 if templateBorderSize > 0:
994 bbox.grow(templateBorderSize)
996 rng = np.random.RandomState(seed)
997 rngNoise = np.random.RandomState(noiseSeed)
998 x0, y0 = bbox.getBegin()
999 xSize, ySize = bbox.getDimensions()
1000 if xLoc is None:
1001 xLoc = rng.rand(nSrc)*(xSize - 2*bufferSize) + bufferSize + x0
1002 else:
1003 if len(xLoc) != nSrc:
1004 raise ValueError("xLoc must have length equal to nSrc. %f supplied vs %f", len(xLoc), nSrc)
1005 if yLoc is None:
1006 yLoc = rng.rand(nSrc)*(ySize - 2*bufferSize) + bufferSize + y0
1007 else:
1008 if len(yLoc) != nSrc:
1009 raise ValueError("yLoc must have length equal to nSrc. %f supplied vs %f", len(yLoc), nSrc)
1011 if flux is None:
1012 flux = (rng.rand(nSrc)*(fluxRange - 1.) + 1.)*fluxLevel
1013 else:
1014 if len(flux) != nSrc:
1015 raise ValueError("flux must have length equal to nSrc. %f supplied vs %f", len(flux), nSrc)
1016 sigmas = [psfSize for src in range(nSrc)]
1017 injectList = list(zip(xLoc, yLoc, flux, sigmas))
1018 skyLevel = 0
1019 # Don't use the built in poisson noise: it modifies the global state of numpy random
1020 modelExposure = plantSources(bbox, kernelSize, skyLevel, injectList, addPoissonNoise=False)
1021 modelExposure.setWcs(makeFakeWcs())
1022 noise = rngNoise.randn(ySize, xSize)*noiseLevel
1023 noise -= np.mean(noise)
1024 modelExposure.variance.array = np.sqrt(np.abs(modelExposure.image.array)) + noiseLevel**2
1025 modelExposure.image.array += noise
1027 # Run source detection to set up the mask plane
1028 detectTestSources(modelExposure, addMaskPlanes=addMaskPlanes)
1029 if clearEdgeMask:
1030 modelExposure.mask &= ~modelExposure.mask.getPlaneBitMask("EDGE")
1031 modelExposure.setPhotoCalib(afwImage.PhotoCalib(calibration, 0., bbox))
1032 if background is not None:
1033 modelExposure.image += background
1034 modelExposure.maskedImage /= calibration
1035 modelExposure.setFilter(afwImage.FilterLabel(band, physicalFilter))
1036 modelExposure.info.setId(seed)
1037 if doApplyCalibration:
1038 modelExposure.maskedImage = modelExposure.photoCalib.calibrateImage(modelExposure.maskedImage)
1040 truth = _fillTruthCatalog(injectList)
1042 return modelExposure, truth
1045def _makeTruthSchema():
1046 """Make a schema for the truth catalog produced by `makeTestImage`.
1048 Returns
1049 -------
1050 keys : `dict` [`str`]
1051 Fields added to the catalog, to make it easier to set them.
1052 schema : `lsst.afw.table.Schema`
1053 Schema to use to make a "truth" SourceCatalog.
1054 Calib, Ap, and Psf flux slots all are set to ``truth_instFlux``.
1055 """
1056 schema = afwTable.SourceTable.makeMinimalSchema()
1057 keys = {}
1058 # Don't use a FluxResultKey so we can manage the flux and err separately.
1059 keys["instFlux"] = schema.addField("truth_instFlux", type=np.float64,
1060 doc="true instFlux", units="count")
1061 keys["instFluxErr"] = schema.addField("truth_instFluxErr", type=np.float64,
1062 doc="true instFluxErr", units="count")
1063 keys["centroid"] = afwTable.Point2DKey.addFields(schema, "truth", "true simulated centroid", "pixel")
1064 schema.addField("truth_flag", "Flag", "truth flux failure flag.")
1065 # Add the flag fields a source selector would need.
1066 schema.addField("sky_source", "Flag", "testing flag.")
1067 schema.addField("base_PixelFlags_flag_interpolated", "Flag", "testing flag.")
1068 schema.addField("base_PixelFlags_flag_saturated", "Flag", "testing flag.")
1069 schema.addField("base_PixelFlags_flag_bad", "Flag", "testing flag.")
1070 schema.addField("base_PixelFlags_flag_edge", "Flag", "testing flag.")
1071 schema.addField("base_PixelFlags_flag_nodata", "Flag", "testing flag.")
1072 schema.addField("base_PsfFlux_flag", "Flag", "testing flag.")
1073 schema.addField("base_ClassificationSizeExtendedness_value", "Flag", "testing flag.")
1074 schema.addField("deblend_nChild", "Flag", "testing flag.")
1075 schema.addField("detect_isPrimary", "Flag", "testing flag.")
1076 schema.addField("calib_psf_used", "Flag", "testing flag.")
1077 schema.getAliasMap().set("slot_Centroid", "truth")
1078 schema.getAliasMap().set("slot_CalibFlux", "truth")
1079 schema.getAliasMap().set("slot_ApFlux", "truth")
1080 schema.getAliasMap().set("slot_PsfFlux", "truth")
1081 schema.addField("glint_trail", "Flag", "testing flag.")
1082 return keys, schema
1085def _fillTruthCatalog(injectList):
1086 """Add injected sources to the truth catalog.
1088 Parameters
1089 ----------
1090 injectList : `list` [`float`]
1091 Sources that were injected; tuples of (x, y, flux, size).
1093 Returns
1094 -------
1095 catalog : `lsst.afw.table.SourceCatalog`
1096 Catalog with centroids and instFlux/instFluxErr values filled in and
1097 appropriate slots set.
1098 """
1099 keys, schema = _makeTruthSchema()
1100 catalog = afwTable.SourceCatalog(schema)
1101 catalog.reserve(len(injectList))
1102 for x, y, flux, size in injectList:
1103 record = catalog.addNew()
1104 keys["centroid"].set(record, geom.PointD(x, y))
1105 keys["instFlux"].set(record, flux)
1106 # Approximate injected errors
1107 keys["instFluxErr"].set(record, 20)
1108 # 5-sigma effective source width
1109 circle = afwGeom.Ellipse(afwGeom.ellipses.Axes(5*size, 5*size, 0), geom.Point2D(x, y))
1110 footprint = afwDetection.Footprint(afwGeom.SpanSet.fromShape(circle))
1111 footprint.addPeak(x, y, flux)
1112 record.setFootprint(footprint)
1114 # Set source records for isolated stars
1115 record["base_ClassificationSizeExtendedness_value"] = 0
1116 record["deblend_nChild"] = 0
1117 record["detect_isPrimary"] = True
1118 record["calib_psf_used"] = True
1120 return catalog
1123def makeStats(badMaskPlanes=None):
1124 """Create a statistics control for configuring calculations on images.
1126 Parameters
1127 ----------
1128 badMaskPlanes : `list` of `str`, optional
1129 List of mask planes to exclude from calculations.
1131 Returns
1132 -------
1133 statsControl : ` lsst.afw.math.StatisticsControl`
1134 Statistics control object for configuring calculations on images.
1135 """
1136 if badMaskPlanes is None:
1137 badMaskPlanes = ("INTRP", "EDGE", "DETECTED", "SAT", "CR",
1138 "BAD", "NO_DATA", "DETECTED_NEGATIVE")
1139 statsControl = afwMath.StatisticsControl()
1140 statsControl.setNumSigmaClip(3.)
1141 statsControl.setNumIter(3)
1142 statsControl.setAndMask(afwImage.Mask.getPlaneBitMask(badMaskPlanes))
1143 return statsControl
1146class CustomCoaddPsf(measAlg.CoaddPsf):
1147 """A custom CoaddPSF that overrides the getAveragePosition method.
1149 It intentionally moves the position off-image to cause a test failure.
1150 """
1151 def getAveragePosition(self):
1152 return geom.Point2D(-10000, -10000)
1155def generate_data_id(*,
1156 tract: int = 9813,
1157 patch: int = 42,
1158 cell_x: int = 4,
1159 cell_y: int = 2,
1160 band: str = "notR",
1161 ) -> DataCoordinate:
1162 """Generate a DataCoordinate instance to use as data_id.
1164 Modified from ``generate_data_id`` in ``lsst.cell_coadds.test_utils``
1166 Parameters
1167 ----------
1168 tract : `int`, optional
1169 Tract ID for the data_id
1170 patch : `int`, optional
1171 Patch ID for the data_id
1172 cell_x : `int`, optional
1173 X index of the cell this patch corresponds to
1174 cell_y : `int`, optional
1175 Y index of the cell this patch corresponds to
1176 band : `str`, optional
1177 Band for the data_id
1179 Returns
1180 -------
1181 data_id : `lsst.daf.butler.DataCoordinate`
1182 An expanded data_id instance.
1183 """
1184 universe = DimensionUniverse()
1186 instrument = universe["instrument"]
1187 instrument_record = instrument.RecordClass(
1188 name="DummyCam",
1189 class_name="lsst.obs.base.instrument_tests.DummyCam",
1190 )
1192 skymap = universe["skymap"]
1193 skymap_record = skymap.RecordClass(name="test_skymap")
1195 band_element = universe["band"]
1196 band_record = band_element.RecordClass(name=band)
1198 physical_filter = universe["physical_filter"]
1199 physical_filter_record = physical_filter.RecordClass(name=band, instrument="test", band=band)
1201 patch_element = universe["patch"]
1202 patch_record = patch_element.RecordClass(
1203 skymap="test_skymap", tract=tract, patch=patch, cell_x=cell_x, cell_y=cell_y
1204 )
1206 # A dictionary with all the relevant records.
1207 record = {
1208 "instrument": instrument_record,
1209 "patch": patch_record,
1210 "tract": 9813,
1211 "band": band_record.name,
1212 "skymap": skymap_record.name,
1213 "physical_filter": physical_filter_record,
1214 }
1216 # A dictionary with all the relevant recordIds.
1217 record_id = record.copy()
1218 for key in (
1219 "instrument",
1220 "physical_filter",
1221 ):
1222 record_id[key] = record_id[key].name
1224 data_id = DataCoordinate.standardize(record_id, universe=universe)
1225 return data_id.expanded(record)
1228def checkMask(mask, sources, excludeMaskPlanes, checkAdjacent=True):
1229 """Exclude sources that are located on or adjacent to masked pixels.
1231 Parameters
1232 ----------
1233 mask : `lsst.afw.image.Mask`
1234 The image mask plane to use to reject sources
1235 based on the location of their centroid on the ccd.
1236 sources : `lsst.afw.table.SourceCatalog`
1237 The source catalog to evaluate.
1238 excludeMaskPlanes : `list` of `str`
1239 List of the names of the mask planes to exclude.
1241 Returns
1242 -------
1243 flags : `numpy.ndarray` of `bool`
1244 Array indicating whether each source in the catalog should be
1245 kept (True) or rejected (False) based on the value of the
1246 mask plane at its location.
1247 """
1248 setExcludeMaskPlanes = [
1249 maskPlane for maskPlane in excludeMaskPlanes if maskPlane in mask.getMaskPlaneDict()
1250 ]
1252 excludePixelMask = mask.getPlaneBitMask(setExcludeMaskPlanes)
1254 xv = (np.rint(sources.getX() - mask.getX0())).astype(int)
1255 yv = (np.rint(sources.getY() - mask.getY0())).astype(int)
1257 flags = np.ones(len(sources), dtype=bool)
1258 if checkAdjacent:
1259 pixRange = (0, -1, 1)
1260 else:
1261 pixRange = (0,)
1262 for j in pixRange:
1263 for i in pixRange:
1264 mv = mask.array[yv + j, xv + i]
1265 flags *= np.bitwise_and(mv, excludePixelMask) == 0
1266 return flags