lsst.meas.deblender.baseline Namespace Reference

Classes
class	PerFootprint
class	PerPeak
	Result of deblending a single Peak within a parent Footprint. More...
class	CachingPsf
	In the PSF fitting code, we request PSF models for all peaks near the one being fit. More...

Functions
def	deblend
	Deblend a single footprint in a maskedImage. More...
def	_fitPsfs
	Fit a PSF + smooth background models to a small region around each peak (plus other peaks that are nearby) in the given list of peaks. More...
def	_fitPsf
	Fit a PSF + smooth background model (linear) to a small region around the peak. More...
def	_handle_flux_at_edge

Function Documentation

def lsst.meas.deblender.baseline._fitPsf (   fp,   fmask,   pk,   pkF,   pkres,   fbb,   peaks,   peaksF,   log,   psf,   psffwhm,   img,   varimg,   psfChisqCut1,   psfChisqCut2,   psfChisqCut2b,   tinyFootprintSize = 2  )

private

Fit a PSF + smooth background model (linear) to a small region around the peak.

Parameters

[in]	fp	Footprint
[in]	fmask	the Mask plane for pixels in the Footprint
[in]	pk	the Peak within the Footprint that we are going to fit a PSF model for
[in]	pkF	the floating-point coordinates of this peak
[in,out]	pkres	a PerPeak object that will hold the results for this peak
[in]	fbb	the bounding-box of this Footprint (a Box2I)
[in]	peaks	a list of all the Peak objects within this Footprint
[in]	peaksF	a python list of the floating-point (x,y) peak locations
[in,out]	log	pex Log object
[in]	psf	afw PSF object
[in]	psffwhm	FWHM of the PSF, in pixels
[in]	img	the Image in which this Footprint finds itself
[in]	varimg	Variance plane
[in]	psfChisqCut*	floats: cuts in chisq-per-pixel at which to accept the PSF model

Definition at line 581 of file baseline.py.

            ):
    """!
    Fit a PSF + smooth background model (linear) to a small region around the peak.

    @param[in]     fp            Footprint
    @param[in]     fmask         the Mask plane for pixels in the Footprint
    @param[in]     pk            the Peak within the Footprint that we are going to fit a PSF model for
    @param[in]     pkF           the floating-point coordinates of this peak
    @param[in,out] pkres         a PerPeak object that will hold the results for this peak
    @param[in]     fbb           the bounding-box of this Footprint (a Box2I)
    @param[in]     peaks         a list of all the Peak objects within this Footprint
    @param[in]     peaksF        a python list of the floating-point (x,y) peak locations
    @param[in,out] log           pex Log object
    @param[in]     psf           afw PSF object
    @param[in]     psffwhm       FWHM of the PSF, in pixels
    @param[in]     img           the Image in which this Footprint finds itself
    @param[in]     varimg        Variance plane
    @param[in]     psfChisqCut*  floats: cuts in chisq-per-pixel at which to accept the PSF model
    """
    import lsstDebug
    # my __name__ is lsst.meas.deblender.baseline
    debugPlots = lsstDebug.Info(__name__).plots
    debugPsf = lsstDebug.Info(__name__).psf

    # The small region is a disk out to R0, plus a ramp with
    # decreasing weight down to R1.
    R0 = int(math.ceil(psffwhm*1.))
    # ramp down to zero weight at this radius...
    R1 = int(math.ceil(psffwhm*1.5))
    cx, cy = pkF.getX(), pkF.getY()
    psfimg = psf.computeImage(cx, cy)
    # R2: distance to neighbouring peak in order to put it into the model
    R2 = R1 + min(psfimg.getWidth(), psfimg.getHeight())/2.

    pbb = psfimg.getBBox()
    pbb.clip(fbb)
    px0, py0 = psfimg.getX0(), psfimg.getY0()

    # Make sure we haven't been given a substitute PSF that's nowhere near where we want, as may occur if
    # "Cannot compute CoaddPsf at point (xx,yy); no input images at that point."
    if not pbb.contains(afwGeom.Point2I(int(cx), int(cy))):
        pkres.setOutOfBounds()
        return

    # The bounding-box of the local region we are going to fit ("stamp")
    xlo = int(math.floor(cx - R1))
    ylo = int(math.floor(cy - R1))
    xhi = int(math.ceil (cx + R1))
    yhi = int(math.ceil (cy + R1))
    stampbb = afwGeom.Box2I(afwGeom.Point2I(xlo, ylo), afwGeom.Point2I(xhi, yhi))
    stampbb.clip(fbb)
    xlo, xhi = stampbb.getMinX(), stampbb.getMaxX()
    ylo, yhi = stampbb.getMinY(), stampbb.getMaxY()
    if xlo > xhi or ylo > yhi:
        log.logdebug('Skipping this peak: out of bounds')
        pkres.setOutOfBounds()
        return

    # drop tiny footprints too?
    if tinyFootprintSize>0 and (((xhi - xlo) < tinyFootprintSize) or
                                ((yhi - ylo) < tinyFootprintSize)):
        log.logdebug('Skipping this peak: tiny footprint / close to edge')
        pkres.setTinyFootprint()
        return

    # find other peaks within range...
    otherpeaks = []
    for pk2, pkF2 in zip(peaks, peaksF):
        if pk2 == pk:
            continue
        if pkF.distanceSquared(pkF2) > R2**2:
            continue
        opsfimg = psf.computeImage(pkF2.getX(), pkF2.getY())
        if not opsfimg.getBBox(afwImage.LOCAL).overlaps(stampbb):
            continue
        otherpeaks.append(opsfimg)
        log.logdebug('%i other peaks within range' % len(otherpeaks))

    # Now we are going to do a least-squares fit for the flux in this
    # PSF, plus a decenter term, a linear sky, and fluxes of nearby
    # sources (assumed point sources).  Build up the matrix...
    # Number of terms -- PSF flux, constant sky, X, Y, + other PSF fluxes
    NT1 = 4 + len(otherpeaks)
    # + PSF dx, dy
    NT2 = NT1 + 2
    # Number of pixels -- at most
    NP = (1 + yhi - ylo)*(1 + xhi - xlo)
    # indices of columns in the "A" matrix.
    I_psf = 0
    I_sky = 1
    I_sky_ramp_x = 2
    I_sky_ramp_y = 3
    # offset of other psf fluxes:
    I_opsf = 4
    I_dx = NT1 + 0
    I_dy = NT1 + 1

    # Build the matrix "A", rhs "b" and weight "w".
    ix0, iy0 = img.getX0(), img.getY0()
    fx0, fy0 = fbb.getMinX(), fbb.getMinY()
    fslice = (slice(ylo-fy0, yhi-fy0+1), slice(xlo-fx0, xhi-fx0+1))
    islice = (slice(ylo-iy0, yhi-iy0+1), slice(xlo-ix0, xhi-ix0+1))
    fmask_sub = fmask .getArray()[fslice]
    var_sub   = varimg.getArray()[islice]
    img_sub   = img   .getArray()[islice]

    # Clip the PSF image to match its bbox
    psfarr = psfimg.getArray()[pbb.getMinY()-py0: 1+pbb.getMaxY()-py0,
                               pbb.getMinX()-px0: 1+pbb.getMaxX()-px0]
    px0, px1 = pbb.getMinX(), pbb.getMaxX()
    py0, py1 = pbb.getMinY(), pbb.getMaxY()

    # Compute the "valid" pixels within our region-of-interest
    valid = (fmask_sub > 0)
    xx, yy = np.arange(xlo, xhi+1), np.arange(ylo, yhi+1)
    RR = ((xx - cx)**2)[np.newaxis, :] + ((yy - cy)**2)[:, np.newaxis]
    valid *= (RR <= R1**2)
    valid *= (var_sub > 0)
    NP = valid.sum()

    if NP == 0:
        log.warn('Skipping peak at (%.1f, %.1f): no unmasked pixels nearby' % (cx, cy))
        pkres.setNoValidPixels()
        return

    # pixel coords of valid pixels
    XX, YY = np.meshgrid(xx, yy)
    ipixes = np.vstack((XX[valid] - xlo, YY[valid] - ylo)).T

    inpsfx = (xx >= px0)*(xx <= px1)
    inpsfy = (yy >= py0)*(yy <= py1)
    inpsf = np.outer(inpsfy, inpsfx)
    indx = np.outer(inpsfy, (xx > px0)*(xx < px1))
    indy = np.outer((yy > py0)*(yy < py1), inpsfx)

    del inpsfx
    del inpsfy

    def _overlap(xlo, xhi, xmin, xmax):
        assert((xlo <= xmax) and (xhi >= xmin) and
               (xlo <= xhi)  and (xmin <= xmax))
        xloclamp = max(xlo, xmin)
        Xlo = xloclamp - xlo
        xhiclamp = min(xhi, xmax)
        Xhi = Xlo + (xhiclamp - xloclamp)
        assert(xloclamp >= 0)
        assert(Xlo >= 0)
        return (xloclamp, xhiclamp+1, Xlo, Xhi+1)

    A = np.zeros((NP, NT2))
    # Constant term
    A[:, I_sky] = 1.
    # Sky slope terms: dx, dy
    A[:, I_sky_ramp_x] = ipixes[:, 0] + (xlo-cx)
    A[:, I_sky_ramp_y] = ipixes[:, 1] + (ylo-cy)

    # whew, grab the valid overlapping PSF pixels
    px0, px1 = pbb.getMinX(), pbb.getMaxX()
    py0, py1 = pbb.getMinY(), pbb.getMaxY()
    sx1, sx2, sx3, sx4 = _overlap(xlo, xhi, px0, px1)
    sy1, sy2, sy3, sy4 = _overlap(ylo, yhi, py0, py1)
    dpx0, dpy0 = px0 - xlo, py0 - ylo
    psf_y_slice = slice(sy3 - dpy0, sy4 - dpy0)
    psf_x_slice = slice(sx3 - dpx0, sx4 - dpx0)
    psfsub = psfarr[psf_y_slice, psf_x_slice]
    vsub = valid[sy1-ylo: sy2-ylo, sx1-xlo: sx2-xlo]
    A[inpsf[valid], I_psf] = psfsub[vsub]

    # PSF dx -- by taking the half-difference of shifted-by-one and
    # shifted-by-minus-one.
    oldsx = (sx1, sx2, sx3, sx4)
    sx1, sx2, sx3, sx4 = _overlap(xlo, xhi, px0+1, px1-1)
    psfsub = (psfarr[psf_y_slice, sx3 - dpx0 + 1: sx4 - dpx0 + 1] -
              psfarr[psf_y_slice, sx3 - dpx0 - 1: sx4 - dpx0 - 1])/2.
    vsub = valid[sy1-ylo: sy2-ylo, sx1-xlo: sx2-xlo]
    A[indx[valid], I_dx] = psfsub[vsub]
    # revert x indices...
    (sx1, sx2, sx3, sx4) = oldsx

    # PSF dy
    sy1, sy2, sy3, sy4 = _overlap(ylo, yhi, py0+1, py1-1)
    psfsub = (psfarr[sy3 - dpy0 + 1: sy4 - dpy0 + 1, psf_x_slice] -
              psfarr[sy3 - dpy0 - 1: sy4 - dpy0 - 1, psf_x_slice])/2.
    vsub = valid[sy1-ylo: sy2-ylo, sx1-xlo: sx2-xlo]
    A[indy[valid], I_dy] = psfsub[vsub]

    # other PSFs...
    for j, opsf in enumerate(otherpeaks):
        obb = opsf.getBBox()
        ino = np.outer((yy >= obb.getMinY())*(yy <= obb.getMaxY()),
                       (xx >= obb.getMinX())*(xx <= obb.getMaxX()))
        dpx0, dpy0 = obb.getMinX() - xlo, obb.getMinY() - ylo
        sx1, sx2, sx3, sx4 = _overlap(xlo, xhi, obb.getMinX(), obb.getMaxX())
        sy1, sy2, sy3, sy4 = _overlap(ylo, yhi, obb.getMinY(), obb.getMaxY())
        opsfarr = opsf.getArray()
        psfsub = opsfarr[sy3 - dpy0 : sy4 - dpy0, sx3 - dpx0: sx4 - dpx0]
        vsub = valid[sy1-ylo: sy2-ylo, sx1-xlo: sx2-xlo]
        A[ino[valid], I_opsf + j] = psfsub[vsub]

    b = img_sub[valid]

    # Weights -- from ramp and image variance map.
    # Ramp weights -- from 1 at R0 down to 0 at R1.
    rw = np.ones_like(RR)
    ii = (RR > R0**2)
    rr = np.sqrt(RR[ii])
    rw[ii] = np.maximum(0, 1. - ((rr - R0)/(R1 - R0)))
    w = np.sqrt(rw[valid]/var_sub[valid])
    # save the effective number of pixels
    sumr = np.sum(rw[valid])
    log.log(-9, 'sumr = %g' % sumr)

    del ii

    Aw  = A*w[:, np.newaxis]
    bw  = b*w

    if debugPlots:
        import pylab as plt
        plt.clf()
        N = NT2 + 2
        R, C = 2, (N+1)/2
        for i in range(NT2):
            im1 = np.zeros((1+yhi-ylo, 1+xhi-xlo))
            im1[ipixes[:, 1], ipixes[:, 0]] = A[:, i]
            plt.subplot(R, C, i+1)
            plt.imshow(im1, interpolation='nearest', origin='lower')
        plt.subplot(R, C, NT2+1)
        im1 = np.zeros((1+yhi-ylo, 1+xhi-xlo))
        im1[ipixes[:, 1], ipixes[:, 0]] = b
        plt.imshow(im1, interpolation='nearest', origin='lower')
        plt.subplot(R, C, NT2+2)
        im1 = np.zeros((1+yhi-ylo, 1+xhi-xlo))
        im1[ipixes[:, 1], ipixes[:, 0]] = w
        plt.imshow(im1, interpolation='nearest', origin='lower')
        plt.savefig('A.png')

    # We do fits with and without the decenter (dx,dy) terms.
    # Since the dx,dy terms are at the end of the matrix,
    # we can do that just by trimming off those elements.
    #
    # The SVD can fail if there are NaNs in the matrices; this should
    # really be handled upstream
    try:
        # NT1 is number of terms without dx,dy;
        # X1 is the result without decenter
        X1, r1, rank1, s1 = np.linalg.lstsq(Aw[:, :NT1], bw)
        # X2 is with decenter
        X2, r2, rank2, s2 = np.linalg.lstsq(Aw, bw)
    except np.linalg.LinAlgError, e:
        log.log(log.WARN, "Failed to fit PSF to child: %s" % e)
        pkres.setPsfFitFailed()
        return

    log.log(-9, 'r1 r2 %s %s' % (r1, r2))

    # r is weighted chi-squared = sum over pixels: ramp * (model -
    # data)**2/sigma**2
    if len(r1) > 0:
        chisq1 = r1[0]
    else:
        chisq1 = 1e30
    if len(r2) > 0:
        chisq2 = r2[0]
    else:
        chisq2 = 1e30
    dof1 = sumr - len(X1)
    dof2 = sumr - len(X2)
    log.log(-9, 'dof1, dof2 %g %g' % (dof1, dof2))

    # This can happen if we're very close to the edge (?)
    if dof1 <= 0 or dof2 <= 0:
        log.logdebug('Skipping this peak: bad DOF %g, %g' % (dof1, dof2))
        pkres.setBadPsfDof()
        return

    q1 = chisq1/dof1
    q2 = chisq2/dof2
    log.logdebug('PSF fits: chisq/dof = %g, %g' % (q1, q2))
    ispsf1 = (q1 < psfChisqCut1)
    ispsf2 = (q2 < psfChisqCut2)

    pkres.psfFit1 = (chisq1, dof1)
    pkres.psfFit2 = (chisq2, dof2)

    # check that the fit PSF spatial derivative terms aren't too big
    if ispsf2:
        fdx, fdy = X2[I_dx], X2[I_dy]
        f0 = X2[I_psf]
        # as a fraction of the PSF flux
        dx = fdx/f0
        dy = fdy/f0
        ispsf2 = ispsf2 and (abs(dx) < 1. and abs(dy) < 1.)
        log.logdebug('isPSF2 -- checking derivatives: dx,dy = %g, %g -> %s' %
                     (dx, dy, str(ispsf2)))
        if not ispsf2:
            pkres.psfFitBigDecenter = True

    # Looks like a shifted PSF: try actually shifting the PSF by that amount
    # and re-evaluate the fit.
    if ispsf2:
        psfimg2 = psf.computeImage(cx + dx, cy + dy)
        # clip
        pbb2 = psfimg2.getBBox()
        pbb2.clip(fbb)

        # Make sure we haven't been given a substitute PSF that's nowhere near where we want, as may occur if
        # "Cannot compute CoaddPsf at point (xx,yy); no input images at that point."
        if not pbb2.contains(afwGeom.Point2I(int(cx + dx), int(cy + dy))):
            ispsf2 = False
        else:
            # clip image to bbox
            px0, py0 = psfimg2.getX0(), psfimg2.getY0()
            psfarr = psfimg2.getArray()[pbb2.getMinY()-py0:1+pbb2.getMaxY()-py0,
                                        pbb2.getMinX()-px0:1+pbb2.getMaxX()-px0]
            px0, py0 = pbb2.getMinX(), pbb2.getMinY()
            px1, py1 = pbb2.getMaxX(), pbb2.getMaxY()

            # yuck!  Update the PSF terms in the least-squares fit matrix.
            Ab = A[:, :NT1]

            sx1, sx2, sx3, sx4 = _overlap(xlo, xhi, px0, px1)
            sy1, sy2, sy3, sy4 = _overlap(ylo, yhi, py0, py1)
            dpx0, dpy0 = px0 - xlo, py0 - ylo
            psfsub = psfarr[sy3-dpy0:sy4-dpy0, sx3-dpx0:sx4-dpx0]
            vsub = valid[sy1-ylo:sy2-ylo, sx1-xlo:sx2-xlo]
            xx, yy = np.arange(xlo, xhi+1), np.arange(ylo, yhi+1)
            inpsf = np.outer((yy >= py0)*(yy <= py1), (xx >= px0)*(xx <= px1))
            Ab[inpsf[valid], I_psf] = psfsub[vsub]

            Aw  = Ab*w[:, np.newaxis]
            # re-solve...
            Xb, rb, rankb, sb = np.linalg.lstsq(Aw, bw)
            if len(rb) > 0:
                chisqb = rb[0]
            else:
                chisqb = 1e30
            dofb = sumr - len(Xb)
            qb = chisqb/dofb
            ispsf2 = (qb < psfChisqCut2b)
            q2 = qb
            X2 = Xb
            log.logdebug('shifted PSF: new chisq/dof = %g; good? %s' %
                         (qb, ispsf2))
            pkres.psfFit3 = (chisqb, dofb)

    # Which one do we keep?
    if (((ispsf1 and ispsf2) and (q2 < q1)) or
        (ispsf2 and not ispsf1)):
        Xpsf = X2
        chisq = chisq2
        dof = dof2
        log.log(-9, 'dof %g' % dof)
        log.logdebug('Keeping shifted-PSF model')
        cx += dx
        cy += dy
        pkres.psfFitWithDecenter = True
    else:
        # (arbitrarily set to X1 when neither fits well)
        Xpsf = X1
        chisq = chisq1
        dof = dof1
        log.log(-9, 'dof %g' % dof)
        log.logdebug('Keeping unshifted PSF model')

    ispsf = (ispsf1 or ispsf2)

    # Save the PSF models in images for posterity.
    if debugPsf:
        SW, SH = 1+xhi-xlo, 1+yhi-ylo
        psfmod = afwImage.ImageF(SW, SH)
        psfmod.setXY0(xlo, ylo)
        psfderivmodm = afwImage.MaskedImageF(SW, SH)
        psfderivmod = psfderivmodm.getImage()
        psfderivmod.setXY0(xlo, ylo)
        model = afwImage.ImageF(SW, SH)
        model.setXY0(xlo, ylo)
        for i in range(len(Xpsf)):
            for (x, y),v in zip(ipixes, A[:, i]*Xpsf[i]):
                ix, iy = int(x), int(y)
                model.set(ix, iy, model.get(ix, iy) + float(v))
                if i in [I_psf, I_dx, I_dy]:
                    psfderivmod.set(ix, iy, psfderivmod.get(ix, iy) + float(v))
        for ii in range(NP):
            x, y = ipixes[ii, :]
            psfmod.set(int(x), int(y), float(A[ii, I_psf]*Xpsf[I_psf]))
        modelfp = afwDet.Footprint(fp.getPeaks().getSchema())
        for (x, y) in ipixes:
            modelfp.addSpan(int(y+ylo), int(x+xlo), int(x+xlo))
        modelfp.normalize()

        pkres.psfFitDebugPsf0Img = psfimg
        pkres.psfFitDebugPsfImg = psfmod
        pkres.psfFitDebugPsfDerivImg = psfderivmod
        pkres.psfFitDebugPsfModel = model
        pkres.psfFitDebugStamp = img.Factory(img, stampbb, True)
        pkres.psfFitDebugValidPix = valid  # numpy array
        pkres.psfFitDebugVar = varimg.Factory(varimg, stampbb, True)
        ww = np.zeros(valid.shape, np.float)
        ww[valid] = w
        pkres.psfFitDebugWeight = ww # numpy
        pkres.psfFitDebugRampWeight = rw


    # Save things we learned about this peak for posterity...
    pkres.psfFitR0 = R0
    pkres.psfFitR1 = R1
    pkres.psfFitStampExtent = (xlo, xhi, ylo, yhi)
    pkres.psfFitCenter = (cx, cy)
    log.log(-9, 'saving chisq,dof %g %g' % (chisq, dof))
    pkres.psfFitBest = (chisq, dof)
    pkres.psfFitParams = Xpsf
    pkres.psfFitFlux = Xpsf[I_psf]
    pkres.psfFitNOthers = len(otherpeaks)

    if ispsf:
        pkres.setDeblendedAsPsf()

        # replace the template image by the PSF + derivatives
        # image.
        log.logdebug('Deblending as PSF; setting template to PSF model')

        # Instantiate the PSF model and clip it to the footprint
        psfimg = psf.computeImage(cx, cy)
        # Scale by fit flux.
        psfimg *= Xpsf[I_psf]
        psfimg = psfimg.convertF()

        # Clip the Footprint to the PSF model image bbox.
        fpcopy = afwDet.Footprint(fp)
        psfbb = psfimg.getBBox()
        fpcopy.clipTo(psfbb)
        bb = fpcopy.getBBox()

        # Copy the part of the PSF model within the clipped footprint.
        psfmod = afwImage.ImageF(bb)
        afwDet.copyWithinFootprintImage(fpcopy, psfimg, psfmod)
        # Save it as our template.
        fpcopy.clipToNonzero(psfmod)
        fpcopy.normalize()
        pkres.setTemplate(psfmod, fpcopy)

        # DEBUG
        pkres.setPsfTemplate(psfmod, fpcopy)