22 __all__ = [
"DeblenderPlugin",
"fitPsfs",
"buildSymmetricTemplates",
"rampFluxAtEdge",
23 "medianSmoothTemplates",
"makeTemplatesMonotonic",
"clipFootprintsToNonzero",
24 "weightTemplates",
"reconstructTemplates",
"apportionFlux"]
35 from .baselineUtils
import BaselineUtilsF
as bUtils
39 """Clips the given *Footprint* to the region in the *Image*
40 containing non-zero values.
42 The clipping drops spans that are
43 totally zero, and moves endpoints to non-zero; it does not
44 split spans that have internal zeros.
48 xImMax = x0 + image.getDimensions().getX()
49 yImMax = y0 + image.getDimensions().getY()
51 arr = image.getArray()
52 for span
in foot.spans:
54 if y < y0
or y > yImMax:
58 xMin = spanX0
if spanX0 >= x0
else x0
59 xMax = spanX1
if spanX1 <= xImMax
else xImMax
60 xarray = np.arange(xMin, xMax+1)[arr[y-y0, xMin-x0:xMax-x0+1] != 0]
65 foot.removeOrphanPeaks()
69 """Class to define plugins for the deblender.
71 The new deblender executes a series of plugins specified by the user.
72 Each plugin defines the function to be executed, the keyword arguments
73 required by the function, and whether or not certain portions of the
74 deblender might need to be rerun as a result of the function.
76 def __init__(self, func, onReset=None, maxIterations=50, **kwargs):
77 """Initialize a deblender plugin
82 Function to run when the plugin is executed.
83 The function should always take
84 `debResult`, a `DeblenderResult` that stores the deblender result,
85 and `log`, an `lsst.log`, as the first two arguments, as well as
86 any additional keyword arguments (that must
87 be specified in ``kwargs``). The function should also return
88 ``modified``, a `bool` that tells the deblender whether
89 or not any templates have been modified by the function.
90 If ``modified==True``, the deblender will go back to step
91 ``onReset``, unless the has already been run ``maxIterations``.
93 Index of the deblender plugin to return to if ``func`` modifies
94 any templates. The default is ``None``, which does not re-run
97 Maximum number of times the deblender will reset when
105 self.
kwargskwargs = kwargs
108 def run(self, debResult, log):
109 """Execute the current plugin
111 Once the plugin has finished, check to see if part of the deblender
112 must be executed again.
114 log.trace(
"Executing %s", self.
funcfunc.__name__)
115 reset = self.
funcfunc(debResult, log, **self.
kwargskwargs)
123 return (
"<Deblender Plugin: func={0}, kwargs={1}".
format(self.
funcfunc.__name__, self.
kwargskwargs))
129 def _setPeakError(debResult, log, pk, cx, cy, filters, msg, flag):
130 """Update the peak in each band with an error
132 This function logs an error that occurs during deblending and sets the
137 debResult: `lsst.meas.deblender.baseline.DeblenderResult`
138 Container for the final deblender results.
140 LSST logger for logging purposes.
142 Number of the peak that failed
144 x coordinate of the peak
146 y coordinate of the peak
148 List of filter names for the exposures
150 Message to display in log traceback
152 Name of the flag to set
158 log.trace(
"Peak %d at (%f,%f):%s", pk, cx, cy, msg)
159 for fidx, f
in enumerate(filters):
160 pkResult = debResult.deblendedParents[f].peaks[pk]
161 getattr(pkResult, flag)()
164 def fitPsfs(debResult, log, psfChisqCut1=1.5, psfChisqCut2=1.5, psfChisqCut2b=1.5, tinyFootprintSize=2):
165 """Fit a PSF + smooth background model (linear) to a small region
168 This function will iterate over all filters in deblender result but does
169 not compare results across filters.
170 DeblendedPeaks that pass the cuts have their templates modified to the
171 PSF + background model and their ``deblendedAsPsf`` property set
174 This will likely be replaced in the future with a function that compares
175 the psf chi-squared cuts so that peaks flagged as point sources will be
176 considered point sources in all bands.
180 debResult: `lsst.meas.deblender.baseline.DeblenderResult`
181 Container for the final deblender results.
183 LSST logger for logging purposes.
184 psfChisqCut*: `float`, optional
185 ``psfChisqCut1`` is the maximum chi-squared-per-degree-of-freedom
186 allowed for a peak to be considered a PSF match without recentering.
187 A fit is also made that includes terms to recenter the PSF.
188 ``psfChisqCut2`` is the same as ``psfChisqCut1`` except it
189 determines the restriction on the fit that includes
191 If the peak is a match for a re-centered PSF, the PSF is
192 repositioned at the new center and
193 the peak footprint is fit again, this time to the new PSF.
194 If the resulting chi-squared-per-degree-of-freedom is less than
195 ``psfChisqCut2b`` then it passes the re-centering algorithm.
196 If the peak passes both the re-centered and fixed position cuts,
197 the better of the two is accepted, but parameters for all three psf
198 fits are stored in the ``DebldendedPeak``.
199 The default for ``psfChisqCut1``, ``psfChisqCut2``, and
200 ``psfChisqCut2b`` is ``1.5``.
201 tinyFootprintSize: `float`, optional
202 The PSF model is shrunk to the size that contains the original
203 footprint. If the bbox of the clipped PSF model for a peak is
204 smaller than ``max(tinyFootprintSize,2)`` then ``tinyFootprint`` for
205 the peak is set to ``True`` and the peak is not fit. The default is 2.
210 If any templates have been assigned to PSF point sources then
211 ``modified`` is ``True``, otherwise it is ``False``.
213 from .baseline
import CachingPsf
216 for fidx
in debResult.filters:
217 dp = debResult.deblendedParents[fidx]
218 peaks = dp.fp.getPeaks()
223 fmask.setXY0(dp.bb.getMinX(), dp.bb.getMinY())
224 dp.fp.spans.setMask(fmask, 1)
229 peakF = [pk.getF()
for pk
in peaks]
231 for pki, (pk, pkres, pkF)
in enumerate(zip(peaks, dp.peaks, peakF)):
232 log.trace(
'Filter %s, Peak %i', fidx, pki)
233 ispsf = _fitPsf(dp.fp, fmask, pk, pkF, pkres, dp.bb, peaks, peakF, log, cpsf, dp.psffwhm,
234 dp.img, dp.varimg, psfChisqCut1, psfChisqCut2, psfChisqCut2b, tinyFootprintSize)
235 modified = modified
or ispsf
239 def _fitPsf(fp, fmask, pk, pkF, pkres, fbb, peaks, peaksF, log, psf, psffwhm,
240 img, varimg, psfChisqCut1, psfChisqCut2, psfChisqCut2b,
243 r"""Fit a PSF + smooth background model (linear) to a small region
246 See fitPsfs for a more thorough description, including all
247 parameters not described below.
251 fp: `afw.detection.Footprint`
252 Footprint containing the Peaks to model.
253 fmask: `afw.image.Mask`
254 The Mask plane for pixels in the Footprint
255 pk: `afw.detection.PeakRecord`
256 The peak within the Footprint that we are going to fit with PSF model
257 pkF: `afw.geom.Point2D`
258 Floating point coordinates of the peak.
259 pkres: `meas.deblender.DeblendedPeak`
260 Peak results object that will hold the results.
261 fbb: `afw.geom.Box2I`
262 Bounding box of ``fp``
263 peaks: `afw.detection.PeakCatalog`
264 Catalog of peaks contained in the parent footprint.
265 peaksF: list of `afw.geom.Point2D`
266 List of floating point coordinates of all of the peaks.
267 psf: list of `afw.detection.Psf`\ s
268 Psf of the ``maskedImage`` for each band.
269 psffwhm: list pf `float`\ s
270 FWHM of the ``maskedImage``\ 's ``psf`` in each band.
271 img: `afw.image.ImageF`
272 The image that contains the footprint.
273 varimg: `afw.image.ImageF`
274 The variance of the image that contains the footprint.
279 Whether or not the peak matches a PSF model.
289 R0 = int(np.ceil(psffwhm*1.))
291 R1 = int(np.ceil(psffwhm*1.5))
292 cx, cy = pkF.getX(), pkF.getY()
293 psfimg = psf.computeImage(cx, cy)
295 R2 = R1 +
min(psfimg.getWidth(), psfimg.getHeight())/2.
297 pbb = psfimg.getBBox()
299 px0, py0 = psfimg.getX0(), psfimg.getY0()
304 pkres.setOutOfBounds()
308 xlo = int(np.floor(cx - R1))
309 ylo = int(np.floor(cy - R1))
310 xhi = int(np.ceil(cx + R1))
311 yhi = int(np.ceil(cy + R1))
314 xlo, xhi = stampbb.getMinX(), stampbb.getMaxX()
315 ylo, yhi = stampbb.getMinY(), stampbb.getMaxY()
316 if xlo > xhi
or ylo > yhi:
317 log.trace(
'Skipping this peak: out of bounds')
318 pkres.setOutOfBounds()
322 if min(stampbb.getWidth(), stampbb.getHeight()) <=
max(tinyFootprintSize, 2):
325 log.trace(
'Skipping this peak: tiny footprint / close to edge')
326 pkres.setTinyFootprint()
331 for pk2, pkF2
in zip(peaks, peaksF):
334 if pkF.distanceSquared(pkF2) > R2**2:
336 opsfimg = psf.computeImage(pkF2.getX(), pkF2.getY())
337 if not opsfimg.getBBox().overlaps(stampbb):
339 otherpeaks.append(opsfimg)
340 log.trace(
'%i other peaks within range', len(otherpeaks))
346 NT1 = 4 + len(otherpeaks)
350 NP = (1 + yhi - ylo)*(1 + xhi - xlo)
362 ix0, iy0 = img.getX0(), img.getY0()
363 fx0, fy0 = fbb.getMinX(), fbb.getMinY()
364 fslice = (slice(ylo-fy0, yhi-fy0+1), slice(xlo-fx0, xhi-fx0+1))
365 islice = (slice(ylo-iy0, yhi-iy0+1), slice(xlo-ix0, xhi-ix0+1))
366 fmask_sub = fmask .getArray()[fslice]
367 var_sub = varimg.getArray()[islice]
368 img_sub = img.getArray()[islice]
371 psfarr = psfimg.getArray()[pbb.getMinY()-py0: 1+pbb.getMaxY()-py0,
372 pbb.getMinX()-px0: 1+pbb.getMaxX()-px0]
373 px0, px1 = pbb.getMinX(), pbb.getMaxX()
374 py0, py1 = pbb.getMinY(), pbb.getMaxY()
377 valid = (fmask_sub > 0)
378 xx, yy = np.arange(xlo, xhi+1), np.arange(ylo, yhi+1)
379 RR = ((xx - cx)**2)[np.newaxis, :] + ((yy - cy)**2)[:, np.newaxis]
380 valid *= (RR <= R1**2)
381 valid *= (var_sub > 0)
385 log.warning(
'Skipping peak at (%.1f, %.1f): no unmasked pixels nearby', cx, cy)
386 pkres.setNoValidPixels()
390 XX, YY = np.meshgrid(xx, yy)
391 ipixes = np.vstack((XX[valid] - xlo, YY[valid] - ylo)).T
393 inpsfx = (xx >= px0)*(xx <= px1)
394 inpsfy = (yy >= py0)*(yy <= py1)
395 inpsf = np.outer(inpsfy, inpsfx)
396 indx = np.outer(inpsfy, (xx > px0)*(xx < px1))
397 indy = np.outer((yy > py0)*(yy < py1), inpsfx)
402 def _overlap(xlo, xhi, xmin, xmax):
403 assert((xlo <= xmax)
and (xhi >= xmin)
404 and (xlo <= xhi)
and (xmin <= xmax))
405 xloclamp =
max(xlo, xmin)
407 xhiclamp =
min(xhi, xmax)
408 Xhi = Xlo + (xhiclamp - xloclamp)
409 assert(xloclamp >= 0)
411 return (xloclamp, xhiclamp+1, Xlo, Xhi+1)
413 A = np.zeros((NP, NT2))
417 A[:, I_sky_ramp_x] = ipixes[:, 0] + (xlo-cx)
418 A[:, I_sky_ramp_y] = ipixes[:, 1] + (ylo-cy)
421 px0, px1 = pbb.getMinX(), pbb.getMaxX()
422 py0, py1 = pbb.getMinY(), pbb.getMaxY()
423 sx1, sx2, sx3, sx4 = _overlap(xlo, xhi, px0, px1)
424 sy1, sy2, sy3, sy4 = _overlap(ylo, yhi, py0, py1)
425 dpx0, dpy0 = px0 - xlo, py0 - ylo
426 psf_y_slice = slice(sy3 - dpy0, sy4 - dpy0)
427 psf_x_slice = slice(sx3 - dpx0, sx4 - dpx0)
428 psfsub = psfarr[psf_y_slice, psf_x_slice]
429 vsub = valid[sy1-ylo: sy2-ylo, sx1-xlo: sx2-xlo]
430 A[inpsf[valid], I_psf] = psfsub[vsub]
434 oldsx = (sx1, sx2, sx3, sx4)
435 sx1, sx2, sx3, sx4 = _overlap(xlo, xhi, px0+1, px1-1)
436 psfsub = (psfarr[psf_y_slice, sx3 - dpx0 + 1: sx4 - dpx0 + 1]
437 - psfarr[psf_y_slice, sx3 - dpx0 - 1: sx4 - dpx0 - 1])/2.
438 vsub = valid[sy1-ylo: sy2-ylo, sx1-xlo: sx2-xlo]
439 A[indx[valid], I_dx] = psfsub[vsub]
441 (sx1, sx2, sx3, sx4) = oldsx
444 sy1, sy2, sy3, sy4 = _overlap(ylo, yhi, py0+1, py1-1)
445 psfsub = (psfarr[sy3 - dpy0 + 1: sy4 - dpy0 + 1, psf_x_slice]
446 - psfarr[sy3 - dpy0 - 1: sy4 - dpy0 - 1, psf_x_slice])/2.
447 vsub = valid[sy1-ylo: sy2-ylo, sx1-xlo: sx2-xlo]
448 A[indy[valid], I_dy] = psfsub[vsub]
451 for j, opsf
in enumerate(otherpeaks):
453 ino = np.outer((yy >= obb.getMinY())*(yy <= obb.getMaxY()),
454 (xx >= obb.getMinX())*(xx <= obb.getMaxX()))
455 dpx0, dpy0 = obb.getMinX() - xlo, obb.getMinY() - ylo
456 sx1, sx2, sx3, sx4 = _overlap(xlo, xhi, obb.getMinX(), obb.getMaxX())
457 sy1, sy2, sy3, sy4 = _overlap(ylo, yhi, obb.getMinY(), obb.getMaxY())
458 opsfarr = opsf.getArray()
459 psfsub = opsfarr[sy3 - dpy0: sy4 - dpy0, sx3 - dpx0: sx4 - dpx0]
460 vsub = valid[sy1-ylo: sy2-ylo, sx1-xlo: sx2-xlo]
461 A[ino[valid], I_opsf + j] = psfsub[vsub]
467 rw = np.ones_like(RR)
470 rw[ii] = np.maximum(0, 1. - ((rr - R0)/(R1 - R0)))
471 w = np.sqrt(rw[valid]/var_sub[valid])
473 sumr = np.sum(rw[valid])
474 log.debug(
'sumr = %g', sumr)
478 Aw = A*w[:, np.newaxis]
487 im1 = np.zeros((1+yhi-ylo, 1+xhi-xlo))
488 im1[ipixes[:, 1], ipixes[:, 0]] = A[:, i]
489 plt.subplot(R, C, i+1)
490 plt.imshow(im1, interpolation=
'nearest', origin=
'lower')
491 plt.subplot(R, C, NT2+1)
492 im1 = np.zeros((1+yhi-ylo, 1+xhi-xlo))
493 im1[ipixes[:, 1], ipixes[:, 0]] = b
494 plt.imshow(im1, interpolation=
'nearest', origin=
'lower')
495 plt.subplot(R, C, NT2+2)
496 im1 = np.zeros((1+yhi-ylo, 1+xhi-xlo))
497 im1[ipixes[:, 1], ipixes[:, 0]] = w
498 plt.imshow(im1, interpolation=
'nearest', origin=
'lower')
510 X1, r1, rank1, s1 = np.linalg.lstsq(Aw[:, :NT1], bw, rcond=-1)
512 X2, r2, rank2, s2 = np.linalg.lstsq(Aw, bw, rcond=-1)
513 except np.linalg.LinAlgError
as e:
514 log.warning(
"Failed to fit PSF to child: %s", e)
515 pkres.setPsfFitFailed()
518 log.debug(
'r1 r2 %s %s', r1, r2)
530 dof1 = sumr - len(X1)
531 dof2 = sumr - len(X2)
532 log.debug(
'dof1, dof2 %g %g', dof1, dof2)
535 if dof1 <= 0
or dof2 <= 0:
536 log.trace(
'Skipping this peak: bad DOF %g, %g', dof1, dof2)
542 log.trace(
'PSF fits: chisq/dof = %g, %g', q1, q2)
543 ispsf1 = (q1 < psfChisqCut1)
544 ispsf2 = (q2 < psfChisqCut2)
546 pkres.psfFit1 = (chisq1, dof1)
547 pkres.psfFit2 = (chisq2, dof2)
551 fdx, fdy = X2[I_dx], X2[I_dy]
556 ispsf2 = ispsf2
and (
abs(dx) < 1.
and abs(dy) < 1.)
557 log.trace(
'isPSF2 -- checking derivatives: dx,dy = %g, %g -> %s', dx, dy, str(ispsf2))
559 pkres.psfFitBigDecenter =
True
564 psfimg2 = psf.computeImage(cx + dx, cy + dy)
566 pbb2 = psfimg2.getBBox()
571 if not pbb2.contains(
geom.Point2I(int(cx + dx), int(cy + dy))):
575 px0, py0 = psfimg2.getX0(), psfimg2.getY0()
576 psfarr = psfimg2.getArray()[pbb2.getMinY()-py0:1+pbb2.getMaxY()-py0,
577 pbb2.getMinX()-px0:1+pbb2.getMaxX()-px0]
578 px0, py0 = pbb2.getMinX(), pbb2.getMinY()
579 px1, py1 = pbb2.getMaxX(), pbb2.getMaxY()
584 sx1, sx2, sx3, sx4 = _overlap(xlo, xhi, px0, px1)
585 sy1, sy2, sy3, sy4 = _overlap(ylo, yhi, py0, py1)
586 dpx0, dpy0 = px0 - xlo, py0 - ylo
587 psfsub = psfarr[sy3-dpy0:sy4-dpy0, sx3-dpx0:sx4-dpx0]
588 vsub = valid[sy1-ylo:sy2-ylo, sx1-xlo:sx2-xlo]
589 xx, yy = np.arange(xlo, xhi+1), np.arange(ylo, yhi+1)
590 inpsf = np.outer((yy >= py0)*(yy <= py1), (xx >= px0)*(xx <= px1))
591 Ab[inpsf[valid], I_psf] = psfsub[vsub]
593 Aw = Ab*w[:, np.newaxis]
595 Xb, rb, rankb, sb = np.linalg.lstsq(Aw, bw, rcond=-1)
600 dofb = sumr - len(Xb)
602 ispsf2 = (qb < psfChisqCut2b)
605 log.trace(
'shifted PSF: new chisq/dof = %g; good? %s', qb, ispsf2)
606 pkres.psfFit3 = (chisqb, dofb)
609 if (((ispsf1
and ispsf2)
and (q2 < q1))
610 or (ispsf2
and not ispsf1)):
614 log.debug(
'dof %g', dof)
615 log.trace(
'Keeping shifted-PSF model')
618 pkres.psfFitWithDecenter =
True
624 log.debug(
'dof %g', dof)
625 log.trace(
'Keeping unshifted PSF model')
627 ispsf = (ispsf1
or ispsf2)
631 SW, SH = 1+xhi-xlo, 1+yhi-ylo
632 psfmod = afwImage.ImageF(SW, SH)
633 psfmod.setXY0(xlo, ylo)
634 psfderivmodm = afwImage.MaskedImageF(SW, SH)
635 psfderivmod = psfderivmodm.getImage()
636 psfderivmod.setXY0(xlo, ylo)
637 model = afwImage.ImageF(SW, SH)
638 model.setXY0(xlo, ylo)
639 for i
in range(len(Xpsf)):
640 for (x, y), v
in zip(ipixes, A[:, i]*Xpsf[i]):
641 ix, iy = int(x), int(y)
642 model.set(ix, iy, model.get(ix, iy) + float(v))
643 if i
in [I_psf, I_dx, I_dy]:
644 psfderivmod.set(ix, iy, psfderivmod.get(ix, iy) + float(v))
647 psfmod.set(int(x), int(y), float(A[ii, I_psf]*Xpsf[I_psf]))
649 for (x, y)
in ipixes:
650 modelfp.addSpan(int(y+ylo), int(x+xlo), int(x+xlo))
653 pkres.psfFitDebugPsf0Img = psfimg
654 pkres.psfFitDebugPsfImg = psfmod
655 pkres.psfFitDebugPsfDerivImg = psfderivmod
656 pkres.psfFitDebugPsfModel = model
657 pkres.psfFitDebugStamp = img.Factory(img, stampbb,
True)
658 pkres.psfFitDebugValidPix = valid
659 pkres.psfFitDebugVar = varimg.Factory(varimg, stampbb,
True)
660 ww = np.zeros(valid.shape, np.float)
662 pkres.psfFitDebugWeight = ww
663 pkres.psfFitDebugRampWeight = rw
668 pkres.psfFitStampExtent = (xlo, xhi, ylo, yhi)
669 pkres.psfFitCenter = (cx, cy)
670 log.debug(
'saving chisq,dof %g %g', chisq, dof)
671 pkres.psfFitBest = (chisq, dof)
672 pkres.psfFitParams = Xpsf
673 pkres.psfFitFlux = Xpsf[I_psf]
674 pkres.psfFitNOthers = len(otherpeaks)
677 pkres.setDeblendedAsPsf()
681 log.trace(
'Deblending as PSF; setting template to PSF model')
684 psfimg = psf.computeImage(cx, cy)
686 psfimg *= Xpsf[I_psf]
687 psfimg = psfimg.convertF()
691 psfbb = psfimg.getBBox()
693 bb = fpcopy.getBBox()
696 psfmod = afwImage.ImageF(bb)
697 fpcopy.spans.copyImage(psfimg, psfmod)
700 pkres.setTemplate(psfmod, fpcopy)
703 pkres.setPsfTemplate(psfmod, fpcopy)
709 """Build a symmetric template for each peak in each filter
711 Given ``maskedImageF``, ``footprint``, and a ``DebldendedPeak``, creates
712 a symmetric template (``templateImage`` and ``templateFootprint``) around
713 the peak for all peaks not flagged as ``skip`` or ``deblendedAsPsf``.
717 debResult: `lsst.meas.deblender.baseline.DeblenderResult`
718 Container for the final deblender results.
720 LSST logger for logging purposes.
721 patchEdges: `bool`, optional
722 If True and if the parent Footprint touches pixels with the
723 ``EDGE`` bit set, then grow the parent Footprint to include
724 all symmetric templates.
729 If any peaks are not skipped or marked as point sources,
730 ``modified`` is ``True. Otherwise ``modified`` is ``False``.
734 for fidx
in debResult.filters:
735 dp = debResult.deblendedParents[fidx]
736 imbb = dp.img.getBBox()
737 log.trace(
'Creating templates for footprint at x0,y0,W,H = %i, %i, %i, %i)', dp.x0, dp.y0, dp.W, dp.H)
739 for peaki, pkres
in enumerate(dp.peaks):
740 log.trace(
'Deblending peak %i of %i', peaki, len(dp.peaks))
743 if pkres.skip
or pkres.deblendedAsPsf:
747 cx, cy = pk.getIx(), pk.getIy()
749 log.trace(
'Peak center is not inside image; skipping %i', pkres.pki)
750 pkres.setOutOfBounds()
752 log.trace(
'computing template for peak %i at (%i, %i)', pkres.pki, cx, cy)
753 timg, tfoot, patched = bUtils.buildSymmetricTemplate(dp.maskedImage, dp.fp, pk, dp.avgNoise,
756 log.trace(
'Peak %i at (%i, %i): failed to build symmetric template', pkres.pki, cx, cy)
757 pkres.setFailedSymmetricTemplate()
765 pkres.setOrigTemplate(timg, tfoot)
766 pkres.setTemplate(timg, tfoot)
771 r"""Adjust flux on the edges of the template footprints.
773 Using the PSF, a peak ``~afw.detection.Footprint`` with pixels on the edge
774 of ``footprint`` is grown by the ``psffwhm*1.5`` and filled in
775 with ramped pixels. The result is a new symmetric footprint
776 template for the peaks near the edge.
780 debResult: `lsst.meas.deblender.baseline.DeblenderResult`
781 Container for the final deblender results.
783 LSST logger for logging purposes.
784 patchEdges: `bool`, optional
785 If True and if the parent Footprint touches pixels with the
786 ``EDGE`` bit set, then grow the parent Footprint to include
787 all symmetric templates.
792 If any peaks have their templates modified to include flux at the
793 edges, ``modified`` is ``True``.
797 for fidx
in debResult.filters:
798 dp = debResult.deblendedParents[fidx]
799 log.trace(
'Checking for significant flux at edge: sigma1=%g', dp.avgNoise)
801 for peaki, pkres
in enumerate(dp.peaks):
802 if pkres.skip
or pkres.deblendedAsPsf:
804 timg, tfoot = pkres.templateImage, pkres.templateFootprint
805 if bUtils.hasSignificantFluxAtEdge(timg, tfoot, 3*dp.avgNoise):
806 log.trace(
"Template %i has significant flux at edge: ramping", pkres.pki)
808 (timg2, tfoot2, patched) = _handle_flux_at_edge(log, dp.psffwhm, timg, tfoot, dp.fp,
809 dp.maskedImage, dp.x0, dp.x1,
810 dp.y0, dp.y1, dp.psf, pkres.peak,
811 dp.avgNoise, patchEdges)
814 and "CoaddPsf" in str(exc)):
815 pkres.setOutOfBounds()
818 pkres.setRampedTemplate(timg2, tfoot2)
821 pkres.setTemplate(timg2, tfoot2)
826 def _handle_flux_at_edge(log, psffwhm, t1, tfoot, fp, maskedImage,
827 x0, x1, y0, y1, psf, pk, sigma1, patchEdges):
828 """Extend a template by the PSF to fill in the footprint.
830 Using the PSF, a footprint that touches the edge is passed to the
831 function and is grown by the psffwhm*1.5 and filled in with
837 LSST logger for logging purposes.
840 t1: `afw.image.ImageF`
841 The image template that contains the footprint to extend.
842 tfoot: `afw.detection.Footprint`
843 Symmetric Footprint to extend.
844 fp: `afw.detection.Footprint`
845 Parent Footprint that is being deblended.
846 maskedImage: `afw.image.MaskedImageF`
847 Full MaskedImage containing the parent footprint ``fp``.
849 Minimum x,y for the bounding box of the footprint ``fp``.
851 Maximum x,y for the bounding box of the footprint ``fp``.
852 psf: `afw.detection.Psf`
854 pk: `afw.detection.PeakRecord`
855 The peak within the Footprint whose footprint is being extended.
857 Estimated noise level in the image.
859 If ``patchEdges==True`` and if the footprint touches pixels with the
860 ``EDGE`` bit set, then for spans whose symmetric mirror are outside
861 the image, the symmetric footprint is grown to include them and their
862 pixel values are stored.
866 t2: `afw.image.ImageF`
867 Image of the extended footprint.
868 tfoot2: `afw.detection.Footprint`
871 If the footprint touches an edge pixel, ``patched`` will be set to
872 ``True``. Otherwise ``patched`` is ``False``.
874 log.trace(
'Found significant flux at template edge.')
884 S = int((S + 0.5)/2)*2 + 1
886 tbb = tfoot.getBBox()
893 fpcopy.setSpans(fpcopy.spans.clippedTo(tbb))
894 fpcopy.removeOrphanPeaks()
895 padim = maskedImage.Factory(tbb)
896 fpcopy.spans.clippedTo(maskedImage.getBBox()).copyMaskedImage(maskedImage, padim)
899 edgepix = bUtils.getSignificantEdgePixels(t1, tfoot, -1e6)
902 xc = int((x0 + x1)/2)
903 yc = int((y0 + y1)/2)
905 pbb = psfim.getBBox()
907 lx, ly = pbb.getMinX(), pbb.getMinY()
908 psfim.setXY0(lx - xc, ly - yc)
909 pbb = psfim.getBBox()
912 if not Sbox.contains(pbb):
914 psfim = psfim.Factory(psfim, Sbox, afwImage.PARENT,
True)
915 pbb = psfim.getBBox()
922 ramped = t1.Factory(tbb)
923 Tout = ramped.getArray()
925 tx0, ty0 = t1.getX0(), t1.getY0()
926 ox0, oy0 = ramped.getX0(), ramped.getY0()
930 for span
in edgepix.getSpans():
932 for x
in range(span.getX0(), span.getX1()+1):
933 slc = (slice(y+py0 - oy0, y+py1+1 - oy0),
934 slice(x+px0 - ox0, x+px1+1 - ox0))
935 Tout[slc] = np.maximum(Tout[slc], Tin[y-ty0, x-tx0]*P)
939 imZeros = (padim.getImage().getArray() == 0)
940 padim.getImage().getArray()[imZeros] = ramped.getArray()[imZeros]
942 t2, tfoot2, patched = bUtils.buildSymmetricTemplate(padim, fpcopy, pk, sigma1,
True, patchEdges)
947 imbb = maskedImage.getBBox()
949 tbb = tfoot2.getBBox()
951 t2 = t2.Factory(t2, tbb, afwImage.PARENT,
True)
953 return t2, tfoot2, patched
957 """Applying median smoothing filter to the template images for every
958 peak in every filter.
962 debResult: `lsst.meas.deblender.baseline.DeblenderResult`
963 Container for the final deblender results.
965 LSST logger for logging purposes.
966 medianFilterHalfSize: `int`, optional
967 Half the box size of the median filter, i.e. a
968 ``medianFilterHalfSize`` of 50 means that each output pixel will
969 be the median of the pixels in a 101 x 101-pixel box in the input
970 image. This parameter is only used when
971 ``medianSmoothTemplate==True``, otherwise it is ignored.
976 Whether or not any templates were modified.
977 This will be ``True`` as long as there is at least one source that
978 is not flagged as a PSF.
982 for fidx
in debResult.filters:
983 dp = debResult.deblendedParents[fidx]
984 for peaki, pkres
in enumerate(dp.peaks):
985 if pkres.skip
or pkres.deblendedAsPsf:
988 timg, tfoot = pkres.templateImage, pkres.templateFootprint
989 filtsize = medianFilterHalfsize*2 + 1
990 if timg.getWidth() >= filtsize
and timg.getHeight() >= filtsize:
991 log.trace(
'Median filtering template %i', pkres.pki)
994 inimg = timg.Factory(timg,
True)
995 bUtils.medianFilter(inimg, timg, medianFilterHalfsize)
997 pkres.setMedianFilteredTemplate(timg, tfoot)
999 log.trace(
'Not median-filtering template %i: size %i x %i smaller than required %i x %i',
1000 pkres.pki, timg.getWidth(), timg.getHeight(), filtsize, filtsize)
1001 pkres.setTemplate(timg, tfoot)
1006 """Make the templates monotonic.
1008 The pixels in the templates are modified such that pixels further
1009 from the peak will have values smaller than those closer to the peak.
1013 debResult: `lsst.meas.deblender.baseline.DeblenderResult`
1014 Container for the final deblender results.
1016 LSST logger for logging purposes.
1021 Whether or not any templates were modified.
1022 This will be ``True`` as long as there is at least one source that
1023 is not flagged as a PSF.
1027 for fidx
in debResult.filters:
1028 dp = debResult.deblendedParents[fidx]
1029 for peaki, pkres
in enumerate(dp.peaks):
1030 if pkres.skip
or pkres.deblendedAsPsf:
1033 timg, tfoot = pkres.templateImage, pkres.templateFootprint
1035 log.trace(
'Making template %i monotonic', pkres.pki)
1036 bUtils.makeMonotonic(timg, pk)
1037 pkres.setTemplate(timg, tfoot)
1042 r"""Clip non-zero spans in the template footprints for every peak in each filter.
1044 Peak ``Footprint``\ s are clipped to the region in the image containing
1045 non-zero values by dropping spans that are completely zero and moving
1046 endpoints to non-zero pixels (but does not split spans that have
1051 debResult: `lsst.meas.deblender.baseline.DeblenderResult`
1052 Container for the final deblender results.
1054 LSST logger for logging purposes.
1059 Whether or not any templates were modified.
1060 This will be ``True`` as long as there is at least one source that
1061 is not flagged as a PSF.
1064 for fidx
in debResult.filters:
1065 dp = debResult.deblendedParents[fidx]
1066 for peaki, pkres
in enumerate(dp.peaks):
1067 if pkres.skip
or pkres.deblendedAsPsf:
1069 timg, tfoot = pkres.templateImage, pkres.templateFootprint
1071 if not tfoot.getBBox().isEmpty()
and tfoot.getBBox() != timg.getBBox(afwImage.PARENT):
1072 timg = timg.Factory(timg, tfoot.getBBox(), afwImage.PARENT,
True)
1073 pkres.setTemplate(timg, tfoot)
1078 """Weight the templates to best fit the observed image in each filter
1080 This function re-weights the templates so that their linear combination
1081 best represents the observed image in that filter.
1082 In the future it may be useful to simultaneously weight all of the
1087 debResult: `lsst.meas.deblender.baseline.DeblenderResult`
1088 Container for the final deblender results.
1090 LSST logger for logging purposes.
1095 ``weightTemplates`` does not actually modify the ``Footprint``
1096 templates other than to add a weight to them, so ``modified``
1097 is always ``False``.
1100 log.trace(
'Weighting templates')
1101 for fidx
in debResult.filters:
1102 _weightTemplates(debResult.deblendedParents[fidx])
1106 def _weightTemplates(dp):
1107 """Weight the templates to best match the parent Footprint in a single
1110 This includes weighting both regular templates and point source templates
1114 dp: `DeblendedParent`
1115 The deblended parent to re-weight
1121 nchild = np.sum([pkres.skip
is False for pkres
in dp.peaks])
1122 A = np.zeros((dp.W*dp.H, nchild))
1123 parentImage = afwImage.ImageF(dp.bb)
1124 afwDet.copyWithinFootprintImage(dp.fp, dp.img, parentImage)
1125 b = parentImage.getArray().ravel()
1128 for pkres
in dp.peaks:
1131 childImage = afwImage.ImageF(dp.bb)
1132 afwDet.copyWithinFootprintImage(dp.fp, pkres.templateImage, childImage)
1133 A[:, index] = childImage.getArray().ravel()
1136 X1, r1, rank1, s1 = np.linalg.lstsq(A, b, rcond=-1)
1141 for pkres
in dp.peaks:
1144 pkres.templateImage *= X1[index]
1145 pkres.setTemplateWeight(X1[index])
1150 """Remove "degenerate templates"
1152 If galaxies have substructure, such as face-on spirals, the process of
1153 identifying peaks can "shred" the galaxy into many pieces. The templates
1154 of shredded galaxies are typically quite similar because they represent
1155 the same galaxy, so we try to identify these "degenerate" peaks
1156 by looking at the inner product (in pixel space) of pairs of templates.
1157 If they are nearly parallel, we only keep one of the peaks and reject
1158 the other. If only one of the peaks is a PSF template, the other template
1159 is used, otherwise the one with the maximum template value is kept.
1163 debResult: `lsst.meas.deblender.baseline.DeblenderResult`
1164 Container for the final deblender results.
1166 LSST logger for logging purposes.
1167 maxTempDotProd: `float`, optional
1168 All dot products between templates greater than ``maxTempDotProd``
1169 will result in one of the templates removed.
1174 If any degenerate templates are found, ``modified`` is ``True``.
1176 log.trace(
'Looking for degnerate templates')
1179 for fidx
in debResult.filters:
1180 dp = debResult.deblendedParents[fidx]
1181 nchild = np.sum([pkres.skip
is False for pkres
in dp.peaks])
1182 indexes = [pkres.pki
for pkres
in dp.peaks
if pkres.skip
is False]
1187 A = np.zeros((nchild, nchild))
1190 for pkres
in dp.peaks:
1194 afwImage.MaskedImageF(pkres.templateImage)))
1195 maxTemplate.append(np.max(pkres.templateImage.getArray()))
1197 for i
in range(nchild):
1198 for j
in range(i + 1):
1199 A[i, j] = heavies[i].
dot(heavies[j])
1202 for i
in range(nchild):
1204 norm = A[i, i]*A[j, j]
1208 A[i, j] /= np.sqrt(norm)
1213 for i
in range(nchild):
1216 if A[i, j] > currentMax:
1217 currentMax = A[i, j]
1218 if currentMax > maxTempDotProd:
1231 reject = indexes[rejectedIndex]
1232 if dp.peaks[keep].deblendedAsPsf
and dp.peaks[reject].deblendedAsPsf
is False:
1233 keep = indexes[rejectedIndex]
1235 elif dp.peaks[keep].deblendedAsPsf
is False and dp.peaks[reject].deblendedAsPsf:
1236 reject = indexes[rejectedIndex]
1239 if maxTemplate[rejectedIndex] > maxTemplate[i]:
1240 keep = indexes[rejectedIndex]
1242 log.trace(
'Removing object with index %d : %f. Degenerate with %d',
1243 reject, currentMax, keep)
1244 dp.peaks[reject].skip =
True
1245 dp.peaks[reject].degenerate =
True
1250 def apportionFlux(debResult, log, assignStrayFlux=True, strayFluxAssignment='r-to-peak',
1251 strayFluxToPointSources='necessary', clipStrayFluxFraction=0.001,
1252 getTemplateSum=False):
1253 """Apportion flux to all of the peak templates in each filter
1255 Divide the ``maskedImage`` flux amongst all of the templates based
1256 on the fraction of flux assigned to each ``template``.
1257 Leftover "stray flux" is assigned to peaks based on the other parameters.
1261 debResult: `lsst.meas.deblender.baseline.DeblenderResult`
1262 Container for the final deblender results.
1264 LSST logger for logging purposes.
1265 assignStrayFlux: `bool`, optional
1266 If True then flux in the parent footprint that is not covered by any
1267 of the template footprints is assigned to templates based on
1268 their 1/(1+r^2) distance.
1269 How the flux is apportioned is determined by ``strayFluxAssignment``.
1270 strayFluxAssignment: `string`, optional
1271 Determines how stray flux is apportioned.
1273 * ``trim``: Trim stray flux and do not include in any footprints
1274 * ``r-to-peak`` (default): Stray flux is assigned based on
1275 (1/(1+r^2) from the peaks
1276 * ``r-to-footprint``: Stray flux is distributed to the footprints
1277 based on 1/(1+r^2) of the minimum distance from the stray flux
1279 * ``nearest-footprint``: Stray flux is assigned to the footprint
1280 with lowest L-1 (Manhattan) distance to the stray flux
1282 strayFluxToPointSources: `string`, optional
1283 Determines how stray flux is apportioned to point sources
1285 * ``never``: never apportion stray flux to point sources
1286 * ``necessary`` (default): point sources are included only if there
1287 are no extended sources nearby
1288 * ``always``: point sources are always included in
1289 the 1/(1+r^2) splitting
1291 clipStrayFluxFraction: `float`, optional
1292 Minimum stray-flux portion.
1293 Any stray-flux portion less than ``clipStrayFluxFraction`` is
1295 getTemplateSum: `bool`, optional
1296 As part of the flux calculation, the sum of the templates is
1297 calculated. If ``getTemplateSum==True`` then the sum of the
1298 templates is stored in the result (a `DeblendedFootprint`).
1303 Apportion flux always modifies the templates, so ``modified`` is
1304 always ``True``. However, this should likely be the final step and
1305 it is unlikely that any deblender plugins will be re-run.
1307 validStrayPtSrc = [
'never',
'necessary',
'always']
1308 validStrayAssign = [
'r-to-peak',
'r-to-footprint',
'nearest-footprint',
'trim']
1309 if strayFluxToPointSources
not in validStrayPtSrc:
1310 raise ValueError(((
'strayFluxToPointSources: value \"%s\" not in the set of allowed values: ') %
1311 strayFluxToPointSources) + str(validStrayPtSrc))
1312 if strayFluxAssignment
not in validStrayAssign:
1313 raise ValueError(((
'strayFluxAssignment: value \"%s\" not in the set of allowed values: ') %
1314 strayFluxAssignment) + str(validStrayAssign))
1316 for fidx
in debResult.filters:
1317 dp = debResult.deblendedParents[fidx]
1330 bb = dp.fp.getBBox()
1332 for peaki, pkres
in enumerate(dp.peaks):
1335 tmimgs.append(pkres.templateImage)
1336 tfoots.append(pkres.templateFootprint)
1338 dpsf.append(pkres.deblendedAsPsf)
1340 pkx.append(pk.getIx())
1341 pky.append(pk.getIy())
1342 ibi.append(pkres.pki)
1345 log.trace(
'Apportioning flux among %i templates', len(tmimgs))
1346 sumimg = afwImage.ImageF(bb)
1351 if strayFluxAssignment ==
'trim':
1352 assignStrayFlux =
False
1353 strayopts |= bUtils.STRAYFLUX_TRIM
1355 strayopts |= bUtils.ASSIGN_STRAYFLUX
1356 if strayFluxToPointSources ==
'necessary':
1357 strayopts |= bUtils.STRAYFLUX_TO_POINT_SOURCES_WHEN_NECESSARY
1358 elif strayFluxToPointSources ==
'always':
1359 strayopts |= bUtils.STRAYFLUX_TO_POINT_SOURCES_ALWAYS
1361 if strayFluxAssignment ==
'r-to-peak':
1364 elif strayFluxAssignment ==
'r-to-footprint':
1365 strayopts |= bUtils.STRAYFLUX_R_TO_FOOTPRINT
1366 elif strayFluxAssignment ==
'nearest-footprint':
1367 strayopts |= bUtils.STRAYFLUX_NEAREST_FOOTPRINT
1369 portions, strayflux = bUtils.apportionFlux(dp.maskedImage, dp.fp, tmimgs, tfoots, sumimg, dpsf,
1370 pkx, pky, strayopts, clipStrayFluxFraction)
1373 if strayFluxAssignment ==
'trim':
1376 finalSpanSet = finalSpanSet.union(foot.spans)
1377 dp.fp.setSpans(finalSpanSet)
1381 debResult.setTemplateSums(sumimg, fidx)
1385 for j, (pk, pkres)
in enumerate(zip(dp.fp.getPeaks(), dp.peaks)):
1388 pkres.setFluxPortion(portions[ii])
1393 stray = strayflux[ii]
1398 pkres.setStrayFlux(stray)
1401 for j, (pk, pkres)
in enumerate(zip(dp.fp.getPeaks(), dp.peaks)):
1405 for foot, add
in [(pkres.templateFootprint,
True), (pkres.origFootprint,
True),
1406 (pkres.strayFlux,
False)]:
1409 pks = foot.getPeaks()
A range of pixels within one row of an Image.
A compact representation of a collection of pixels.
Represent a 2-dimensional array of bitmask pixels.
An integer coordinate rectangle.
def __init__(self, func, onReset=None, maxIterations=50, **kwargs)
def run(self, debResult, log)
Provides consistent interface for LSST exceptions.
Reports invalid arguments.
HeavyFootprint< ImagePixelT, MaskPixelT, VariancePixelT > makeHeavyFootprint(Footprint const &foot, lsst::afw::image::MaskedImage< ImagePixelT, MaskPixelT, VariancePixelT > const &img, HeavyFootprintCtrl const *ctrl=nullptr)
Create a HeavyFootprint with footprint defined by the given Footprint and pixel values from the given...
def dot(symb, c, r, frame=None, size=2, ctype=None, origin=afwImage.PARENT, *args, **kwargs)
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
def makeTemplatesMonotonic(debResult, log)
def fitPsfs(debResult, log, psfChisqCut1=1.5, psfChisqCut2=1.5, psfChisqCut2b=1.5, tinyFootprintSize=2)
def weightTemplates(debResult, log)
def buildSymmetricTemplates(debResult, log, patchEdges=False, setOrigTemplate=True)
def reconstructTemplates(debResult, log, maxTempDotProd=0.5)
def medianSmoothTemplates(debResult, log, medianFilterHalfsize=2)
def clipFootprintToNonzeroImpl(foot, image)
def apportionFlux(debResult, log, assignStrayFlux=True, strayFluxAssignment='r-to-peak', strayFluxToPointSources='necessary', clipStrayFluxFraction=0.001, getTemplateSum=False)
def clipFootprintsToNonzero(debResult, log)
def rampFluxAtEdge(debResult, log, patchEdges=False)
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Angle abs(Angle const &a)