31 from .baselineUtils
import BaselineUtilsF
as bUtils
36 Clips the given *Footprint* to the region in the *Image*
37 containing non-zero values. The clipping drops spans that are
38 totally zero, and moves endpoints to non-zero; it does not
39 split spans that have internal zeros.
43 xImMax = x0 + image.getDimensions().getX()
44 yImMax = y0 + image.getDimensions().getY()
46 arr = image.getArray()
47 for span
in foot.spans:
49 if y < y0
or y > yImMax:
53 xMin = spanX0
if spanX0 >= x0
else x0
54 xMax = spanX1
if spanX1 <= xImMax
else xImMax
55 xarray = np.arange(xMin, xMax+1)[arr[y-y0, xMin-x0:xMax-x0+1] != 0]
60 foot.removeOrphanPeaks()
64 """Class to define plugins for the deblender.
66 The new deblender executes a series of plugins specified by the user.
67 Each plugin defines the function to be executed, the keyword arguments required by the function,
68 and whether or not certain portions of the deblender might need to be rerun as a result of
71 def __init__(self, func, onReset=None, maxIterations=50, **kwargs):
72 """Initialize a deblender plugin
77 Function to run when the plugin is executed. The function should always take
78 `debResult`, a `DeblenderResult` that stores the deblender result, and
79 `log`, an `lsst.log`, as the first two arguments, as well as any additional
80 keyword arguments (that must be specified in ``kwargs``).
81 The function should also return ``modified``, a `bool` that tells the deblender whether
82 or not any templates have been modified by the function.
83 If ``modified==True``, the deblender will go back to step ``onReset``,
84 unless the has already been run ``maxIterations``.
86 Index of the deblender plugin to return to if ``func`` modifies any templates.
87 The default is ``None``, which does not re-run any plugins.
89 Maximum number of times the deblender will reset when the current plugin
99 def run(self, debResult, log):
100 """Execute the current plugin
102 Once the plugin has finished, check to see if part of the deblender must be executed again.
104 log.trace(
"Executing %s", self.
func.__name__)
105 reset = self.
func(debResult, log, **self.
kwargs)
113 return (
"<Deblender Plugin: func={0}, kwargs={1}".
format(self.
func.__name__, self.
kwargs))
119 def _setPeakError(debResult, log, pk, cx, cy, filters, msg, flag):
120 """Update the peak in each band with an error
122 This function logs an error that occurs during deblending and sets the
127 debResult: `lsst.meas.deblender.baseline.DeblenderResult`
128 Container for the final deblender results.
130 LSST logger for logging purposes.
132 Number of the peak that failed
134 x coordinate of the peak
136 y coordinate of the peak
138 List of filter names for the exposures
140 Message to display in log traceback
142 Name of the flag to set
148 log.trace(
"Peak {0} at ({1},{2}):{3}".
format(pk, cx, cy, msg))
149 for fidx, f
in enumerate(filters):
150 pkResult = debResult.deblendedParents[f].peaks[pk]
151 getattr(pkResult, flag)()
154 def fitPsfs(debResult, log, psfChisqCut1=1.5, psfChisqCut2=1.5, psfChisqCut2b=1.5, tinyFootprintSize=2):
155 """Fit a PSF + smooth background model (linear) to a small region around each peak
157 This function will iterate over all filters in deblender result but does not compare
158 results across filters.
159 DeblendedPeaks that pass the cuts have their templates modified to the PSF + background model
160 and their ``deblendedAsPsf`` property set to ``True``.
162 This will likely be replaced in the future with a function that compares the psf chi-squared cuts
163 so that peaks flagged as point sources will be considered point sources in all bands.
167 debResult: `lsst.meas.deblender.baseline.DeblenderResult`
168 Container for the final deblender results.
170 LSST logger for logging purposes.
171 psfChisqCut*: `float`, optional
172 ``psfChisqCut1`` is the maximum chi-squared-per-degree-of-freedom allowed for a peak to
173 be considered a PSF match without recentering.
174 A fit is also made that includes terms to recenter the PSF.
175 ``psfChisqCut2`` is the same as ``psfChisqCut1`` except it determines the restriction on the
176 fit that includes recentering terms.
177 If the peak is a match for a re-centered PSF, the PSF is repositioned at the new center and
178 the peak footprint is fit again, this time to the new PSF.
179 If the resulting chi-squared-per-degree-of-freedom is less than ``psfChisqCut2b`` then it
180 passes the re-centering algorithm.
181 If the peak passes both the re-centered and fixed position cuts, the better of the two is accepted,
182 but parameters for all three psf fits are stored in the ``DebldendedPeak``.
183 The default for ``psfChisqCut1``, ``psfChisqCut2``, and ``psfChisqCut2b`` is ``1.5``.
184 tinyFootprintSize: `float`, optional
185 The PSF model is shrunk to the size that contains the original footprint.
186 If the bbox of the clipped PSF model for a peak is smaller than ``max(tinyFootprintSize,2)``
187 then ``tinyFootprint`` for the peak is set to ``True`` and the peak is not fit.
193 If any templates have been assigned to PSF point sources then ``modified`` is ``True``,
194 otherwise it is ``False``.
196 from .baseline
import CachingPsf
199 for fidx
in debResult.filters:
200 dp = debResult.deblendedParents[fidx]
201 peaks = dp.fp.getPeaks()
206 fmask.setXY0(dp.bb.getMinX(), dp.bb.getMinY())
207 dp.fp.spans.setMask(fmask, 1)
212 peakF = [pk.getF()
for pk
in peaks]
214 for pki, (pk, pkres, pkF)
in enumerate(zip(peaks, dp.peaks, peakF)):
215 log.trace(
'Filter %s, Peak %i', fidx, pki)
216 ispsf = _fitPsf(dp.fp, fmask, pk, pkF, pkres, dp.bb, peaks, peakF, log, cpsf, dp.psffwhm,
217 dp.img, dp.varimg, psfChisqCut1, psfChisqCut2, psfChisqCut2b, tinyFootprintSize)
218 modified = modified
or ispsf
222 def _fitPsf(fp, fmask, pk, pkF, pkres, fbb, peaks, peaksF, log, psf, psffwhm,
223 img, varimg, psfChisqCut1, psfChisqCut2, psfChisqCut2b,
226 """Fit a PSF + smooth background model (linear) to a small region around a peak.
228 See fitPsfs for a more thorough description, including all parameters not described below.
232 fp: `afw.detection.Footprint`
233 Footprint containing the Peaks to model.
234 fmask: `afw.image.Mask`
235 The Mask plane for pixels in the Footprint
236 pk: `afw.detection.PeakRecord`
237 The peak within the Footprint that we are going to fit with PSF model
238 pkF: `afw.geom.Point2D`
239 Floating point coordinates of the peak.
240 pkres: `meas.deblender.DeblendedPeak`
241 Peak results object that will hold the results.
242 fbb: `afw.geom.Box2I`
243 Bounding box of ``fp``
244 peaks: `afw.detection.PeakCatalog`
245 Catalog of peaks contained in the parent footprint.
246 peaksF: list of `afw.geom.Point2D`
247 List of floating point coordinates of all of the peaks.
248 psf: list of `afw.detection.Psf`s
249 Psf of the ``maskedImage`` for each band.
250 psffwhm: list pf `float`s
251 FWHM of the ``maskedImage``'s ``psf`` in each band.
252 img: `afw.image.ImageF`
253 The image that contains the footprint.
254 varimg: `afw.image.ImageF`
255 The variance of the image that contains the footprint.
260 Whether or not the peak matches a PSF model.
270 R0 = int(np.ceil(psffwhm*1.))
272 R1 = int(np.ceil(psffwhm*1.5))
273 cx, cy = pkF.getX(), pkF.getY()
274 psfimg = psf.computeImage(cx, cy)
276 R2 = R1 +
min(psfimg.getWidth(), psfimg.getHeight())/2.
278 pbb = psfimg.getBBox()
280 px0, py0 = psfimg.getX0(), psfimg.getY0()
285 pkres.setOutOfBounds()
289 xlo = int(np.floor(cx - R1))
290 ylo = int(np.floor(cy - R1))
291 xhi = int(np.ceil(cx + R1))
292 yhi = int(np.ceil(cy + R1))
295 xlo, xhi = stampbb.getMinX(), stampbb.getMaxX()
296 ylo, yhi = stampbb.getMinY(), stampbb.getMaxY()
297 if xlo > xhi
or ylo > yhi:
298 log.trace(
'Skipping this peak: out of bounds')
299 pkres.setOutOfBounds()
303 if min(stampbb.getWidth(), stampbb.getHeight()) <=
max(tinyFootprintSize, 2):
306 log.trace(
'Skipping this peak: tiny footprint / close to edge')
307 pkres.setTinyFootprint()
312 for pk2, pkF2
in zip(peaks, peaksF):
315 if pkF.distanceSquared(pkF2) > R2**2:
317 opsfimg = psf.computeImage(pkF2.getX(), pkF2.getY())
318 if not opsfimg.getBBox().overlaps(stampbb):
320 otherpeaks.append(opsfimg)
321 log.trace(
'%i other peaks within range', len(otherpeaks))
327 NT1 = 4 + len(otherpeaks)
331 NP = (1 + yhi - ylo)*(1 + xhi - xlo)
343 ix0, iy0 = img.getX0(), img.getY0()
344 fx0, fy0 = fbb.getMinX(), fbb.getMinY()
345 fslice = (slice(ylo-fy0, yhi-fy0+1), slice(xlo-fx0, xhi-fx0+1))
346 islice = (slice(ylo-iy0, yhi-iy0+1), slice(xlo-ix0, xhi-ix0+1))
347 fmask_sub = fmask .getArray()[fslice]
348 var_sub = varimg.getArray()[islice]
349 img_sub = img.getArray()[islice]
352 psfarr = psfimg.getArray()[pbb.getMinY()-py0: 1+pbb.getMaxY()-py0,
353 pbb.getMinX()-px0: 1+pbb.getMaxX()-px0]
354 px0, px1 = pbb.getMinX(), pbb.getMaxX()
355 py0, py1 = pbb.getMinY(), pbb.getMaxY()
358 valid = (fmask_sub > 0)
359 xx, yy = np.arange(xlo, xhi+1), np.arange(ylo, yhi+1)
360 RR = ((xx - cx)**2)[np.newaxis, :] + ((yy - cy)**2)[:, np.newaxis]
361 valid *= (RR <= R1**2)
362 valid *= (var_sub > 0)
366 log.warn(
'Skipping peak at (%.1f, %.1f): no unmasked pixels nearby', cx, cy)
367 pkres.setNoValidPixels()
371 XX, YY = np.meshgrid(xx, yy)
372 ipixes = np.vstack((XX[valid] - xlo, YY[valid] - ylo)).T
374 inpsfx = (xx >= px0)*(xx <= px1)
375 inpsfy = (yy >= py0)*(yy <= py1)
376 inpsf = np.outer(inpsfy, inpsfx)
377 indx = np.outer(inpsfy, (xx > px0)*(xx < px1))
378 indy = np.outer((yy > py0)*(yy < py1), inpsfx)
383 def _overlap(xlo, xhi, xmin, xmax):
384 assert((xlo <= xmax)
and (xhi >= xmin)
385 and (xlo <= xhi)
and (xmin <= xmax))
386 xloclamp =
max(xlo, xmin)
388 xhiclamp =
min(xhi, xmax)
389 Xhi = Xlo + (xhiclamp - xloclamp)
390 assert(xloclamp >= 0)
392 return (xloclamp, xhiclamp+1, Xlo, Xhi+1)
394 A = np.zeros((NP, NT2))
398 A[:, I_sky_ramp_x] = ipixes[:, 0] + (xlo-cx)
399 A[:, I_sky_ramp_y] = ipixes[:, 1] + (ylo-cy)
402 px0, px1 = pbb.getMinX(), pbb.getMaxX()
403 py0, py1 = pbb.getMinY(), pbb.getMaxY()
404 sx1, sx2, sx3, sx4 = _overlap(xlo, xhi, px0, px1)
405 sy1, sy2, sy3, sy4 = _overlap(ylo, yhi, py0, py1)
406 dpx0, dpy0 = px0 - xlo, py0 - ylo
407 psf_y_slice = slice(sy3 - dpy0, sy4 - dpy0)
408 psf_x_slice = slice(sx3 - dpx0, sx4 - dpx0)
409 psfsub = psfarr[psf_y_slice, psf_x_slice]
410 vsub = valid[sy1-ylo: sy2-ylo, sx1-xlo: sx2-xlo]
411 A[inpsf[valid], I_psf] = psfsub[vsub]
415 oldsx = (sx1, sx2, sx3, sx4)
416 sx1, sx2, sx3, sx4 = _overlap(xlo, xhi, px0+1, px1-1)
417 psfsub = (psfarr[psf_y_slice, sx3 - dpx0 + 1: sx4 - dpx0 + 1]
418 - psfarr[psf_y_slice, sx3 - dpx0 - 1: sx4 - dpx0 - 1])/2.
419 vsub = valid[sy1-ylo: sy2-ylo, sx1-xlo: sx2-xlo]
420 A[indx[valid], I_dx] = psfsub[vsub]
422 (sx1, sx2, sx3, sx4) = oldsx
425 sy1, sy2, sy3, sy4 = _overlap(ylo, yhi, py0+1, py1-1)
426 psfsub = (psfarr[sy3 - dpy0 + 1: sy4 - dpy0 + 1, psf_x_slice]
427 - psfarr[sy3 - dpy0 - 1: sy4 - dpy0 - 1, psf_x_slice])/2.
428 vsub = valid[sy1-ylo: sy2-ylo, sx1-xlo: sx2-xlo]
429 A[indy[valid], I_dy] = psfsub[vsub]
432 for j, opsf
in enumerate(otherpeaks):
434 ino = np.outer((yy >= obb.getMinY())*(yy <= obb.getMaxY()),
435 (xx >= obb.getMinX())*(xx <= obb.getMaxX()))
436 dpx0, dpy0 = obb.getMinX() - xlo, obb.getMinY() - ylo
437 sx1, sx2, sx3, sx4 = _overlap(xlo, xhi, obb.getMinX(), obb.getMaxX())
438 sy1, sy2, sy3, sy4 = _overlap(ylo, yhi, obb.getMinY(), obb.getMaxY())
439 opsfarr = opsf.getArray()
440 psfsub = opsfarr[sy3 - dpy0: sy4 - dpy0, sx3 - dpx0: sx4 - dpx0]
441 vsub = valid[sy1-ylo: sy2-ylo, sx1-xlo: sx2-xlo]
442 A[ino[valid], I_opsf + j] = psfsub[vsub]
448 rw = np.ones_like(RR)
451 rw[ii] = np.maximum(0, 1. - ((rr - R0)/(R1 - R0)))
452 w = np.sqrt(rw[valid]/var_sub[valid])
454 sumr = np.sum(rw[valid])
455 log.debug(
'sumr = %g', sumr)
459 Aw = A*w[:, np.newaxis]
468 im1 = np.zeros((1+yhi-ylo, 1+xhi-xlo))
469 im1[ipixes[:, 1], ipixes[:, 0]] = A[:, i]
470 plt.subplot(R, C, i+1)
471 plt.imshow(im1, interpolation=
'nearest', origin=
'lower')
472 plt.subplot(R, C, NT2+1)
473 im1 = np.zeros((1+yhi-ylo, 1+xhi-xlo))
474 im1[ipixes[:, 1], ipixes[:, 0]] = b
475 plt.imshow(im1, interpolation=
'nearest', origin=
'lower')
476 plt.subplot(R, C, NT2+2)
477 im1 = np.zeros((1+yhi-ylo, 1+xhi-xlo))
478 im1[ipixes[:, 1], ipixes[:, 0]] = w
479 plt.imshow(im1, interpolation=
'nearest', origin=
'lower')
491 X1, r1, rank1, s1 = np.linalg.lstsq(Aw[:, :NT1], bw, rcond=-1)
493 X2, r2, rank2, s2 = np.linalg.lstsq(Aw, bw, rcond=-1)
494 except np.linalg.LinAlgError
as e:
495 log.warn(
"Failed to fit PSF to child: %s", e)
496 pkres.setPsfFitFailed()
499 log.debug(
'r1 r2 %s %s', r1, r2)
511 dof1 = sumr - len(X1)
512 dof2 = sumr - len(X2)
513 log.debug(
'dof1, dof2 %g %g', dof1, dof2)
516 if dof1 <= 0
or dof2 <= 0:
517 log.trace(
'Skipping this peak: bad DOF %g, %g', dof1, dof2)
523 log.trace(
'PSF fits: chisq/dof = %g, %g', q1, q2)
524 ispsf1 = (q1 < psfChisqCut1)
525 ispsf2 = (q2 < psfChisqCut2)
527 pkres.psfFit1 = (chisq1, dof1)
528 pkres.psfFit2 = (chisq2, dof2)
532 fdx, fdy = X2[I_dx], X2[I_dy]
537 ispsf2 = ispsf2
and (
abs(dx) < 1.
and abs(dy) < 1.)
538 log.trace(
'isPSF2 -- checking derivatives: dx,dy = %g, %g -> %s', dx, dy,
str(ispsf2))
540 pkres.psfFitBigDecenter =
True
545 psfimg2 = psf.computeImage(cx + dx, cy + dy)
547 pbb2 = psfimg2.getBBox()
552 if not pbb2.contains(
geom.Point2I(int(cx + dx), int(cy + dy))):
556 px0, py0 = psfimg2.getX0(), psfimg2.getY0()
557 psfarr = psfimg2.getArray()[pbb2.getMinY()-py0:1+pbb2.getMaxY()-py0,
558 pbb2.getMinX()-px0:1+pbb2.getMaxX()-px0]
559 px0, py0 = pbb2.getMinX(), pbb2.getMinY()
560 px1, py1 = pbb2.getMaxX(), pbb2.getMaxY()
565 sx1, sx2, sx3, sx4 = _overlap(xlo, xhi, px0, px1)
566 sy1, sy2, sy3, sy4 = _overlap(ylo, yhi, py0, py1)
567 dpx0, dpy0 = px0 - xlo, py0 - ylo
568 psfsub = psfarr[sy3-dpy0:sy4-dpy0, sx3-dpx0:sx4-dpx0]
569 vsub = valid[sy1-ylo:sy2-ylo, sx1-xlo:sx2-xlo]
570 xx, yy = np.arange(xlo, xhi+1), np.arange(ylo, yhi+1)
571 inpsf = np.outer((yy >= py0)*(yy <= py1), (xx >= px0)*(xx <= px1))
572 Ab[inpsf[valid], I_psf] = psfsub[vsub]
574 Aw = Ab*w[:, np.newaxis]
576 Xb, rb, rankb, sb = np.linalg.lstsq(Aw, bw, rcond=-1)
581 dofb = sumr - len(Xb)
583 ispsf2 = (qb < psfChisqCut2b)
586 log.trace(
'shifted PSF: new chisq/dof = %g; good? %s', qb, ispsf2)
587 pkres.psfFit3 = (chisqb, dofb)
590 if (((ispsf1
and ispsf2)
and (q2 < q1))
591 or (ispsf2
and not ispsf1)):
595 log.debug(
'dof %g', dof)
596 log.trace(
'Keeping shifted-PSF model')
599 pkres.psfFitWithDecenter =
True
605 log.debug(
'dof %g', dof)
606 log.trace(
'Keeping unshifted PSF model')
608 ispsf = (ispsf1
or ispsf2)
612 SW, SH = 1+xhi-xlo, 1+yhi-ylo
613 psfmod = afwImage.ImageF(SW, SH)
614 psfmod.setXY0(xlo, ylo)
615 psfderivmodm = afwImage.MaskedImageF(SW, SH)
616 psfderivmod = psfderivmodm.getImage()
617 psfderivmod.setXY0(xlo, ylo)
618 model = afwImage.ImageF(SW, SH)
619 model.setXY0(xlo, ylo)
620 for i
in range(len(Xpsf)):
621 for (x, y), v
in zip(ipixes, A[:, i]*Xpsf[i]):
622 ix, iy = int(x), int(y)
623 model.set(ix, iy, model.get(ix, iy) + float(v))
624 if i
in [I_psf, I_dx, I_dy]:
625 psfderivmod.set(ix, iy, psfderivmod.get(ix, iy) + float(v))
628 psfmod.set(int(x), int(y), float(A[ii, I_psf]*Xpsf[I_psf]))
630 for (x, y)
in ipixes:
631 modelfp.addSpan(int(y+ylo), int(x+xlo), int(x+xlo))
634 pkres.psfFitDebugPsf0Img = psfimg
635 pkres.psfFitDebugPsfImg = psfmod
636 pkres.psfFitDebugPsfDerivImg = psfderivmod
637 pkres.psfFitDebugPsfModel = model
638 pkres.psfFitDebugStamp = img.Factory(img, stampbb,
True)
639 pkres.psfFitDebugValidPix = valid
640 pkres.psfFitDebugVar = varimg.Factory(varimg, stampbb,
True)
641 ww = np.zeros(valid.shape, np.float)
643 pkres.psfFitDebugWeight = ww
644 pkres.psfFitDebugRampWeight = rw
649 pkres.psfFitStampExtent = (xlo, xhi, ylo, yhi)
650 pkres.psfFitCenter = (cx, cy)
651 log.debug(
'saving chisq,dof %g %g', chisq, dof)
652 pkres.psfFitBest = (chisq, dof)
653 pkres.psfFitParams = Xpsf
654 pkres.psfFitFlux = Xpsf[I_psf]
655 pkres.psfFitNOthers = len(otherpeaks)
658 pkres.setDeblendedAsPsf()
662 log.trace(
'Deblending as PSF; setting template to PSF model')
665 psfimg = psf.computeImage(cx, cy)
667 psfimg *= Xpsf[I_psf]
668 psfimg = psfimg.convertF()
672 psfbb = psfimg.getBBox()
674 bb = fpcopy.getBBox()
677 psfmod = afwImage.ImageF(bb)
678 fpcopy.spans.copyImage(psfimg, psfmod)
681 pkres.setTemplate(psfmod, fpcopy)
684 pkres.setPsfTemplate(psfmod, fpcopy)
690 """Build a symmetric template for each peak in each filter
692 Given ``maskedImageF``, ``footprint``, and a ``DebldendedPeak``, creates a symmetric template
693 (``templateImage`` and ``templateFootprint``) around the peak for all peaks not flagged as
694 ``skip`` or ``deblendedAsPsf``.
698 debResult: `lsst.meas.deblender.baseline.DeblenderResult`
699 Container for the final deblender results.
701 LSST logger for logging purposes.
702 patchEdges: `bool`, optional
703 If True and if the parent Footprint touches pixels with the ``EDGE`` bit set,
704 then grow the parent Footprint to include all symmetric templates.
709 If any peaks are not skipped or marked as point sources, ``modified`` is ``True.
710 Otherwise ``modified`` is ``False``.
714 for fidx
in debResult.filters:
715 dp = debResult.deblendedParents[fidx]
716 imbb = dp.img.getBBox()
717 log.trace(
'Creating templates for footprint at x0,y0,W,H = %i, %i, %i, %i)', dp.x0, dp.y0, dp.W, dp.H)
719 for peaki, pkres
in enumerate(dp.peaks):
720 log.trace(
'Deblending peak %i of %i', peaki, len(dp.peaks))
723 if pkres.skip
or pkres.deblendedAsPsf:
727 cx, cy = pk.getIx(), pk.getIy()
729 log.trace(
'Peak center is not inside image; skipping %i', pkres.pki)
730 pkres.setOutOfBounds()
732 log.trace(
'computing template for peak %i at (%i, %i)', pkres.pki, cx, cy)
733 timg, tfoot, patched = bUtils.buildSymmetricTemplate(dp.maskedImage, dp.fp, pk, dp.avgNoise,
736 log.trace(
'Peak %i at (%i, %i): failed to build symmetric template', pkres.pki, cx, cy)
737 pkres.setFailedSymmetricTemplate()
745 pkres.setOrigTemplate(timg, tfoot)
746 pkres.setTemplate(timg, tfoot)
751 """Adjust flux on the edges of the template footprints.
753 Using the PSF, a peak ``Footprint`` with pixels on the edge of ``footprint``
754 is grown by the ``psffwhm``*1.5 and filled in with ramped pixels.
755 The result is a new symmetric footprint template for the peaks near the edge.
759 debResult: `lsst.meas.deblender.baseline.DeblenderResult`
760 Container for the final deblender results.
762 LSST logger for logging purposes.
763 patchEdges: `bool`, optional
764 If True and if the parent Footprint touches pixels with the ``EDGE`` bit set,
765 then grow the parent Footprint to include all symmetric templates.
770 If any peaks have their templates modified to include flux at the edges,
771 ``modified`` is ``True``.
775 for fidx
in debResult.filters:
776 dp = debResult.deblendedParents[fidx]
777 log.trace(
'Checking for significant flux at edge: sigma1=%g', dp.avgNoise)
779 for peaki, pkres
in enumerate(dp.peaks):
780 if pkres.skip
or pkres.deblendedAsPsf:
782 timg, tfoot = pkres.templateImage, pkres.templateFootprint
783 if bUtils.hasSignificantFluxAtEdge(timg, tfoot, 3*dp.avgNoise):
784 log.trace(
"Template %i has significant flux at edge: ramping", pkres.pki)
786 (timg2, tfoot2, patched) = _handle_flux_at_edge(log, dp.psffwhm, timg, tfoot, dp.fp,
787 dp.maskedImage, dp.x0, dp.x1,
788 dp.y0, dp.y1, dp.psf, pkres.peak,
789 dp.avgNoise, patchEdges)
792 and "CoaddPsf" in str(exc)):
793 pkres.setOutOfBounds()
796 pkres.setRampedTemplate(timg2, tfoot2)
799 pkres.setTemplate(timg2, tfoot2)
804 def _handle_flux_at_edge(log, psffwhm, t1, tfoot, fp, maskedImage,
805 x0, x1, y0, y1, psf, pk, sigma1, patchEdges):
806 """Extend a template by the PSF to fill in the footprint.
808 Using the PSF, a footprint that touches the edge is passed to the function
809 and is grown by the psffwhm*1.5 and filled in with ramped pixels.
814 LSST logger for logging purposes.
817 t1: `afw.image.ImageF`
818 The image template that contains the footprint to extend.
819 tfoot: `afw.detection.Footprint`
820 Symmetric Footprint to extend.
821 fp: `afw.detection.Footprint`
822 Parent Footprint that is being deblended.
823 maskedImage: `afw.image.MaskedImageF`
824 Full MaskedImage containing the parent footprint ``fp``.
826 Minimum x,y for the bounding box of the footprint ``fp``.
828 Maximum x,y for the bounding box of the footprint ``fp``.
829 psf: `afw.detection.Psf`
831 pk: `afw.detection.PeakRecord`
832 The peak within the Footprint whose footprint is being extended.
834 Estimated noise level in the image.
836 If ``patchEdges==True`` and if the footprint touches pixels with the
837 ``EDGE`` bit set, then for spans whose symmetric mirror are outside the
838 image, the symmetric footprint is grown to include them and their
839 pixel values are stored.
843 t2: `afw.image.ImageF`
844 Image of the extended footprint.
845 tfoot2: `afw.detection.Footprint`
848 If the footprint touches an edge pixel, ``patched`` will be set to ``True``.
849 Otherwise ``patched`` is ``False``.
851 log.trace(
'Found significant flux at template edge.')
861 S = int((S + 0.5)/2)*2 + 1
863 tbb = tfoot.getBBox()
870 fpcopy.setSpans(fpcopy.spans.clippedTo(tbb))
871 fpcopy.removeOrphanPeaks()
872 padim = maskedImage.Factory(tbb)
873 fpcopy.spans.clippedTo(maskedImage.getBBox()).copyMaskedImage(maskedImage, padim)
876 edgepix = bUtils.getSignificantEdgePixels(t1, tfoot, -1e6)
879 xc = int((x0 + x1)/2)
880 yc = int((y0 + y1)/2)
882 pbb = psfim.getBBox()
884 lx, ly = pbb.getMinX(), pbb.getMinY()
885 psfim.setXY0(lx - xc, ly - yc)
886 pbb = psfim.getBBox()
889 if not Sbox.contains(pbb):
891 psfim = psfim.Factory(psfim, Sbox, afwImage.PARENT,
True)
892 pbb = psfim.getBBox()
899 ramped = t1.Factory(tbb)
900 Tout = ramped.getArray()
902 tx0, ty0 = t1.getX0(), t1.getY0()
903 ox0, oy0 = ramped.getX0(), ramped.getY0()
907 for span
in edgepix.getSpans():
909 for x
in range(span.getX0(), span.getX1()+1):
910 slc = (slice(y+py0 - oy0, y+py1+1 - oy0),
911 slice(x+px0 - ox0, x+px1+1 - ox0))
912 Tout[slc] = np.maximum(Tout[slc], Tin[y-ty0, x-tx0]*P)
916 imZeros = (padim.getImage().getArray() == 0)
917 padim.getImage().getArray()[imZeros] = ramped.getArray()[imZeros]
919 t2, tfoot2, patched = bUtils.buildSymmetricTemplate(padim, fpcopy, pk, sigma1,
True, patchEdges)
924 imbb = maskedImage.getBBox()
926 tbb = tfoot2.getBBox()
928 t2 = t2.Factory(t2, tbb, afwImage.PARENT,
True)
930 return t2, tfoot2, patched
934 """Applying median smoothing filter to the template images for every peak in every filter.
938 debResult: `lsst.meas.deblender.baseline.DeblenderResult`
939 Container for the final deblender results.
941 LSST logger for logging purposes.
942 medianFilterHalfSize: `int`, optional
943 Half the box size of the median filter, i.e. a ``medianFilterHalfSize`` of 50 means that
944 each output pixel will be the median of the pixels in a 101 x 101-pixel box in the input image.
945 This parameter is only used when ``medianSmoothTemplate==True``, otherwise it is ignored.
950 Whether or not any templates were modified.
951 This will be ``True`` as long as there is at least one source that is not flagged as a PSF.
955 for fidx
in debResult.filters:
956 dp = debResult.deblendedParents[fidx]
957 for peaki, pkres
in enumerate(dp.peaks):
958 if pkres.skip
or pkres.deblendedAsPsf:
961 timg, tfoot = pkres.templateImage, pkres.templateFootprint
962 filtsize = medianFilterHalfsize*2 + 1
963 if timg.getWidth() >= filtsize
and timg.getHeight() >= filtsize:
964 log.trace(
'Median filtering template %i', pkres.pki)
967 inimg = timg.Factory(timg,
True)
968 bUtils.medianFilter(inimg, timg, medianFilterHalfsize)
970 pkres.setMedianFilteredTemplate(timg, tfoot)
972 log.trace(
'Not median-filtering template %i: size %i x %i smaller than required %i x %i',
973 pkres.pki, timg.getWidth(), timg.getHeight(), filtsize, filtsize)
974 pkres.setTemplate(timg, tfoot)
979 """Make the templates monotonic.
981 The pixels in the templates are modified such that pixels further from the peak will
982 have values smaller than those closer to the peak.
986 debResult: `lsst.meas.deblender.baseline.DeblenderResult`
987 Container for the final deblender results.
989 LSST logger for logging purposes.
994 Whether or not any templates were modified.
995 This will be ``True`` as long as there is at least one source that is not flagged as a PSF.
999 for fidx
in debResult.filters:
1000 dp = debResult.deblendedParents[fidx]
1001 for peaki, pkres
in enumerate(dp.peaks):
1002 if pkres.skip
or pkres.deblendedAsPsf:
1005 timg, tfoot = pkres.templateImage, pkres.templateFootprint
1007 log.trace(
'Making template %i monotonic', pkres.pki)
1008 bUtils.makeMonotonic(timg, pk)
1009 pkres.setTemplate(timg, tfoot)
1014 """Clip non-zero spans in the template footprints for every peak in each filter.
1016 Peak ``Footprint``s are clipped to the region in the image containing non-zero values
1017 by dropping spans that are completely zero and moving endpoints to non-zero pixels
1018 (but does not split spans that have internal zeros).
1022 debResult: `lsst.meas.deblender.baseline.DeblenderResult`
1023 Container for the final deblender results.
1025 LSST logger for logging purposes.
1030 Whether or not any templates were modified.
1031 This will be ``True`` as long as there is at least one source that is not flagged as a PSF.
1034 for fidx
in debResult.filters:
1035 dp = debResult.deblendedParents[fidx]
1036 for peaki, pkres
in enumerate(dp.peaks):
1037 if pkres.skip
or pkres.deblendedAsPsf:
1039 timg, tfoot = pkres.templateImage, pkres.templateFootprint
1041 if not tfoot.getBBox().isEmpty()
and tfoot.getBBox() != timg.getBBox(afwImage.PARENT):
1042 timg = timg.Factory(timg, tfoot.getBBox(), afwImage.PARENT,
True)
1043 pkres.setTemplate(timg, tfoot)
1048 """Weight the templates to best fit the observed image in each filter
1050 This function re-weights the templates so that their linear combination best represents
1051 the observed image in that filter.
1052 In the future it may be useful to simultaneously weight all of the filters together.
1056 debResult: `lsst.meas.deblender.baseline.DeblenderResult`
1057 Container for the final deblender results.
1059 LSST logger for logging purposes.
1064 ``weightTemplates`` does not actually modify the ``Footprint`` templates other than
1065 to add a weight to them, so ``modified`` is always ``False``.
1068 log.trace(
'Weighting templates')
1069 for fidx
in debResult.filters:
1070 _weightTemplates(debResult.deblendedParents[fidx])
1074 def _weightTemplates(dp):
1075 """Weight the templates to best match the parent Footprint in a single filter
1077 This includes weighting both regular templates and point source templates
1081 dp: `DeblendedParent`
1082 The deblended parent to re-weight
1088 nchild = np.sum([pkres.skip
is False for pkres
in dp.peaks])
1089 A = np.zeros((dp.W*dp.H, nchild))
1090 parentImage = afwImage.ImageF(dp.bb)
1091 afwDet.copyWithinFootprintImage(dp.fp, dp.img, parentImage)
1092 b = parentImage.getArray().ravel()
1095 for pkres
in dp.peaks:
1098 childImage = afwImage.ImageF(dp.bb)
1099 afwDet.copyWithinFootprintImage(dp.fp, pkres.templateImage, childImage)
1100 A[:, index] = childImage.getArray().ravel()
1103 X1, r1, rank1, s1 = np.linalg.lstsq(A, b, rcond=-1)
1108 for pkres
in dp.peaks:
1111 pkres.templateImage *= X1[index]
1112 pkres.setTemplateWeight(X1[index])
1117 """Remove "degenerate templates"
1119 If galaxies have substructure, such as face-on spirals, the process of identifying peaks can
1120 "shred" the galaxy into many pieces. The templates of shredded galaxies are typically quite
1121 similar because they represent the same galaxy, so we try to identify these "degenerate" peaks
1122 by looking at the inner product (in pixel space) of pairs of templates.
1123 If they are nearly parallel, we only keep one of the peaks and reject the other.
1124 If only one of the peaks is a PSF template, the other template is used,
1125 otherwise the one with the maximum template value is kept.
1129 debResult: `lsst.meas.deblender.baseline.DeblenderResult`
1130 Container for the final deblender results.
1132 LSST logger for logging purposes.
1133 maxTempDotProd: `float`, optional
1134 All dot products between templates greater than ``maxTempDotProd`` will result in one
1135 of the templates removed.
1140 If any degenerate templates are found, ``modified`` is ``True``.
1142 log.trace(
'Looking for degnerate templates')
1145 for fidx
in debResult.filters:
1146 dp = debResult.deblendedParents[fidx]
1147 nchild = np.sum([pkres.skip
is False for pkres
in dp.peaks])
1148 indexes = [pkres.pki
for pkres
in dp.peaks
if pkres.skip
is False]
1153 A = np.zeros((nchild, nchild))
1156 for pkres
in dp.peaks:
1160 afwImage.MaskedImageF(pkres.templateImage)))
1161 maxTemplate.append(np.max(pkres.templateImage.getArray()))
1163 for i
in range(nchild):
1164 for j
in range(i + 1):
1165 A[i, j] = heavies[i].
dot(heavies[j])
1168 for i
in range(nchild):
1170 norm = A[i, i]*A[j, j]
1174 A[i, j] /= np.sqrt(norm)
1179 for i
in range(nchild):
1182 if A[i, j] > currentMax:
1183 currentMax = A[i, j]
1184 if currentMax > maxTempDotProd:
1197 reject = indexes[rejectedIndex]
1198 if dp.peaks[keep].deblendedAsPsf
and dp.peaks[reject].deblendedAsPsf
is False:
1199 keep = indexes[rejectedIndex]
1201 elif dp.peaks[keep].deblendedAsPsf
is False and dp.peaks[reject].deblendedAsPsf:
1202 reject = indexes[rejectedIndex]
1205 if maxTemplate[rejectedIndex] > maxTemplate[i]:
1206 keep = indexes[rejectedIndex]
1208 log.trace(
'Removing object with index %d : %f. Degenerate with %d' % (reject, currentMax,
1210 dp.peaks[reject].skip =
True
1211 dp.peaks[reject].degenerate =
True
1216 def apportionFlux(debResult, log, assignStrayFlux=True, strayFluxAssignment='r-to-peak',
1217 strayFluxToPointSources='necessary', clipStrayFluxFraction=0.001,
1218 getTemplateSum=False):
1219 """Apportion flux to all of the peak templates in each filter
1221 Divide the ``maskedImage`` flux amongst all of the templates based on the fraction of
1222 flux assigned to each ``template``.
1223 Leftover "stray flux" is assigned to peaks based on the other parameters.
1227 debResult: `lsst.meas.deblender.baseline.DeblenderResult`
1228 Container for the final deblender results.
1230 LSST logger for logging purposes.
1231 assignStrayFlux: `bool`, optional
1232 If True then flux in the parent footprint that is not covered by any of the
1233 template footprints is assigned to templates based on their 1/(1+r^2) distance.
1234 How the flux is apportioned is determined by ``strayFluxAssignment``.
1235 strayFluxAssignment: `string`, optional
1236 Determines how stray flux is apportioned.
1237 * ``trim``: Trim stray flux and do not include in any footprints
1238 * ``r-to-peak`` (default): Stray flux is assigned based on (1/(1+r^2) from the peaks
1239 * ``r-to-footprint``: Stray flux is distributed to the footprints based on 1/(1+r^2) of the
1240 minimum distance from the stray flux to footprint
1241 * ``nearest-footprint``: Stray flux is assigned to the footprint with lowest L-1 (Manhattan)
1242 distance to the stray flux
1243 strayFluxToPointSources: `string`, optional
1244 Determines how stray flux is apportioned to point sources
1245 * ``never``: never apportion stray flux to point sources
1246 * ``necessary`` (default): point sources are included only if there are no extended sources nearby
1247 * ``always``: point sources are always included in the 1/(1+r^2) splitting
1248 clipStrayFluxFraction: `float`, optional
1249 Minimum stray-flux portion.
1250 Any stray-flux portion less than ``clipStrayFluxFraction`` is clipped to zero.
1251 getTemplateSum: `bool`, optional
1252 As part of the flux calculation, the sum of the templates is calculated.
1253 If ``getTemplateSum==True`` then the sum of the templates is stored in the result
1254 (a `DeblendedFootprint`).
1259 Apportion flux always modifies the templates, so ``modified`` is always ``True``.
1260 However, this should likely be the final step and it is unlikely that
1261 any deblender plugins will be re-run.
1263 validStrayPtSrc = [
'never',
'necessary',
'always']
1264 validStrayAssign = [
'r-to-peak',
'r-to-footprint',
'nearest-footprint',
'trim']
1265 if strayFluxToPointSources
not in validStrayPtSrc:
1266 raise ValueError(((
'strayFluxToPointSources: value \"%s\" not in the set of allowed values: ') %
1267 strayFluxToPointSources) +
str(validStrayPtSrc))
1268 if strayFluxAssignment
not in validStrayAssign:
1269 raise ValueError(((
'strayFluxAssignment: value \"%s\" not in the set of allowed values: ') %
1270 strayFluxAssignment) +
str(validStrayAssign))
1272 for fidx
in debResult.filters:
1273 dp = debResult.deblendedParents[fidx]
1286 bb = dp.fp.getBBox()
1288 for peaki, pkres
in enumerate(dp.peaks):
1291 tmimgs.append(pkres.templateImage)
1292 tfoots.append(pkres.templateFootprint)
1294 dpsf.append(pkres.deblendedAsPsf)
1296 pkx.append(pk.getIx())
1297 pky.append(pk.getIy())
1298 ibi.append(pkres.pki)
1301 log.trace(
'Apportioning flux among %i templates', len(tmimgs))
1302 sumimg = afwImage.ImageF(bb)
1307 if strayFluxAssignment ==
'trim':
1308 assignStrayFlux =
False
1309 strayopts |= bUtils.STRAYFLUX_TRIM
1311 strayopts |= bUtils.ASSIGN_STRAYFLUX
1312 if strayFluxToPointSources ==
'necessary':
1313 strayopts |= bUtils.STRAYFLUX_TO_POINT_SOURCES_WHEN_NECESSARY
1314 elif strayFluxToPointSources ==
'always':
1315 strayopts |= bUtils.STRAYFLUX_TO_POINT_SOURCES_ALWAYS
1317 if strayFluxAssignment ==
'r-to-peak':
1320 elif strayFluxAssignment ==
'r-to-footprint':
1321 strayopts |= bUtils.STRAYFLUX_R_TO_FOOTPRINT
1322 elif strayFluxAssignment ==
'nearest-footprint':
1323 strayopts |= bUtils.STRAYFLUX_NEAREST_FOOTPRINT
1325 portions, strayflux = bUtils.apportionFlux(dp.maskedImage, dp.fp, tmimgs, tfoots, sumimg, dpsf,
1326 pkx, pky, strayopts, clipStrayFluxFraction)
1329 if strayFluxAssignment ==
'trim':
1332 finalSpanSet = finalSpanSet.union(foot.spans)
1333 dp.fp.setSpans(finalSpanSet)
1337 debResult.setTemplateSums(sumimg, fidx)
1341 for j, (pk, pkres)
in enumerate(zip(dp.fp.getPeaks(), dp.peaks)):
1344 pkres.setFluxPortion(portions[ii])
1349 stray = strayflux[ii]
1354 pkres.setStrayFlux(stray)
1357 for j, (pk, pkres)
in enumerate(zip(dp.fp.getPeaks(), dp.peaks)):
1361 for foot, add
in [(pkres.templateFootprint,
True), (pkres.origFootprint,
True),
1362 (pkres.strayFlux,
False)]:
1365 pks = foot.getPeaks()