LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
LSSTDataManagementBasePackage
Classes | Functions
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]fpFootprint
[in]fmaskthe Mask plane for pixels in the Footprint
[in]pkthe Peak within the Footprint that we are going to fit a PSF model for
[in]pkFthe floating-point coordinates of this peak
[in,out]pkresa PerPeak object that will hold the results for this peak
[in]fbbthe bounding-box of this Footprint (a Box2I)
[in]peaksa list of all the Peak objects within this Footprint
[in]peaksFa python list of the floating-point (x,y) peak locations
[in,out]logpex Log object
[in]psfafw PSF object
[in]psffwhmFWHM of the PSF, in pixels
[in]imgthe Image in which this Footprint finds itself
[in]varimgVariance plane
[in]psfChisqCut*floats: cuts in chisq-per-pixel at which to accept the PSF model

Definition at line 581 of file baseline.py.

582  ):
583  """!
584  Fit a PSF + smooth background model (linear) to a small region around the peak.
585 
586  @param[in] fp Footprint
587  @param[in] fmask the Mask plane for pixels in the Footprint
588  @param[in] pk the Peak within the Footprint that we are going to fit a PSF model for
589  @param[in] pkF the floating-point coordinates of this peak
590  @param[in,out] pkres a PerPeak object that will hold the results for this peak
591  @param[in] fbb the bounding-box of this Footprint (a Box2I)
592  @param[in] peaks a list of all the Peak objects within this Footprint
593  @param[in] peaksF a python list of the floating-point (x,y) peak locations
594  @param[in,out] log pex Log object
595  @param[in] psf afw PSF object
596  @param[in] psffwhm FWHM of the PSF, in pixels
597  @param[in] img the Image in which this Footprint finds itself
598  @param[in] varimg Variance plane
599  @param[in] psfChisqCut* floats: cuts in chisq-per-pixel at which to accept the PSF model
600  """
601  import lsstDebug
602  # my __name__ is lsst.meas.deblender.baseline
603  debugPlots = lsstDebug.Info(__name__).plots
604  debugPsf = lsstDebug.Info(__name__).psf
605 
606  # The small region is a disk out to R0, plus a ramp with
607  # decreasing weight down to R1.
608  R0 = int(math.ceil(psffwhm*1.))
609  # ramp down to zero weight at this radius...
610  R1 = int(math.ceil(psffwhm*1.5))
611  cx, cy = pkF.getX(), pkF.getY()
612  psfimg = psf.computeImage(cx, cy)
613  # R2: distance to neighbouring peak in order to put it into the model
614  R2 = R1 + min(psfimg.getWidth(), psfimg.getHeight())/2.
615 
616  pbb = psfimg.getBBox()
617  pbb.clip(fbb)
618  px0, py0 = psfimg.getX0(), psfimg.getY0()
619 
620  # Make sure we haven't been given a substitute PSF that's nowhere near where we want, as may occur if
621  # "Cannot compute CoaddPsf at point (xx,yy); no input images at that point."
622  if not pbb.contains(afwGeom.Point2I(int(cx), int(cy))):
623  pkres.setOutOfBounds()
624  return
625 
626  # The bounding-box of the local region we are going to fit ("stamp")
627  xlo = int(math.floor(cx - R1))
628  ylo = int(math.floor(cy - R1))
629  xhi = int(math.ceil (cx + R1))
630  yhi = int(math.ceil (cy + R1))
631  stampbb = afwGeom.Box2I(afwGeom.Point2I(xlo, ylo), afwGeom.Point2I(xhi, yhi))
632  stampbb.clip(fbb)
633  xlo, xhi = stampbb.getMinX(), stampbb.getMaxX()
634  ylo, yhi = stampbb.getMinY(), stampbb.getMaxY()
635  if xlo > xhi or ylo > yhi:
636  log.logdebug('Skipping this peak: out of bounds')
637  pkres.setOutOfBounds()
638  return
639 
640  # drop tiny footprints too?
641  if tinyFootprintSize>0 and (((xhi - xlo) < tinyFootprintSize) or
642  ((yhi - ylo) < tinyFootprintSize)):
643  log.logdebug('Skipping this peak: tiny footprint / close to edge')
644  pkres.setTinyFootprint()
645  return
646 
647  # find other peaks within range...
648  otherpeaks = []
649  for pk2, pkF2 in zip(peaks, peaksF):
650  if pk2 == pk:
651  continue
652  if pkF.distanceSquared(pkF2) > R2**2:
653  continue
654  opsfimg = psf.computeImage(pkF2.getX(), pkF2.getY())
655  if not opsfimg.getBBox(afwImage.LOCAL).overlaps(stampbb):
656  continue
657  otherpeaks.append(opsfimg)
658  log.logdebug('%i other peaks within range' % len(otherpeaks))
659 
660  # Now we are going to do a least-squares fit for the flux in this
661  # PSF, plus a decenter term, a linear sky, and fluxes of nearby
662  # sources (assumed point sources). Build up the matrix...
663  # Number of terms -- PSF flux, constant sky, X, Y, + other PSF fluxes
664  NT1 = 4 + len(otherpeaks)
665  # + PSF dx, dy
666  NT2 = NT1 + 2
667  # Number of pixels -- at most
668  NP = (1 + yhi - ylo)*(1 + xhi - xlo)
669  # indices of columns in the "A" matrix.
670  I_psf = 0
671  I_sky = 1
672  I_sky_ramp_x = 2
673  I_sky_ramp_y = 3
674  # offset of other psf fluxes:
675  I_opsf = 4
676  I_dx = NT1 + 0
677  I_dy = NT1 + 1
678 
679  # Build the matrix "A", rhs "b" and weight "w".
680  ix0, iy0 = img.getX0(), img.getY0()
681  fx0, fy0 = fbb.getMinX(), fbb.getMinY()
682  fslice = (slice(ylo-fy0, yhi-fy0+1), slice(xlo-fx0, xhi-fx0+1))
683  islice = (slice(ylo-iy0, yhi-iy0+1), slice(xlo-ix0, xhi-ix0+1))
684  fmask_sub = fmask .getArray()[fslice]
685  var_sub = varimg.getArray()[islice]
686  img_sub = img .getArray()[islice]
687 
688  # Clip the PSF image to match its bbox
689  psfarr = psfimg.getArray()[pbb.getMinY()-py0: 1+pbb.getMaxY()-py0,
690  pbb.getMinX()-px0: 1+pbb.getMaxX()-px0]
691  px0, px1 = pbb.getMinX(), pbb.getMaxX()
692  py0, py1 = pbb.getMinY(), pbb.getMaxY()
693 
694  # Compute the "valid" pixels within our region-of-interest
695  valid = (fmask_sub > 0)
696  xx, yy = np.arange(xlo, xhi+1), np.arange(ylo, yhi+1)
697  RR = ((xx - cx)**2)[np.newaxis, :] + ((yy - cy)**2)[:, np.newaxis]
698  valid *= (RR <= R1**2)
699  valid *= (var_sub > 0)
700  NP = valid.sum()
701 
702  if NP == 0:
703  log.warn('Skipping peak at (%.1f, %.1f): no unmasked pixels nearby' % (cx, cy))
704  pkres.setNoValidPixels()
705  return
706 
707  # pixel coords of valid pixels
708  XX, YY = np.meshgrid(xx, yy)
709  ipixes = np.vstack((XX[valid] - xlo, YY[valid] - ylo)).T
710 
711  inpsfx = (xx >= px0)*(xx <= px1)
712  inpsfy = (yy >= py0)*(yy <= py1)
713  inpsf = np.outer(inpsfy, inpsfx)
714  indx = np.outer(inpsfy, (xx > px0)*(xx < px1))
715  indy = np.outer((yy > py0)*(yy < py1), inpsfx)
716 
717  del inpsfx
718  del inpsfy
719 
720  def _overlap(xlo, xhi, xmin, xmax):
721  assert((xlo <= xmax) and (xhi >= xmin) and
722  (xlo <= xhi) and (xmin <= xmax))
723  xloclamp = max(xlo, xmin)
724  Xlo = xloclamp - xlo
725  xhiclamp = min(xhi, xmax)
726  Xhi = Xlo + (xhiclamp - xloclamp)
727  assert(xloclamp >= 0)
728  assert(Xlo >= 0)
729  return (xloclamp, xhiclamp+1, Xlo, Xhi+1)
730 
731  A = np.zeros((NP, NT2))
732  # Constant term
733  A[:, I_sky] = 1.
734  # Sky slope terms: dx, dy
735  A[:, I_sky_ramp_x] = ipixes[:, 0] + (xlo-cx)
736  A[:, I_sky_ramp_y] = ipixes[:, 1] + (ylo-cy)
737 
738  # whew, grab the valid overlapping PSF pixels
739  px0, px1 = pbb.getMinX(), pbb.getMaxX()
740  py0, py1 = pbb.getMinY(), pbb.getMaxY()
741  sx1, sx2, sx3, sx4 = _overlap(xlo, xhi, px0, px1)
742  sy1, sy2, sy3, sy4 = _overlap(ylo, yhi, py0, py1)
743  dpx0, dpy0 = px0 - xlo, py0 - ylo
744  psf_y_slice = slice(sy3 - dpy0, sy4 - dpy0)
745  psf_x_slice = slice(sx3 - dpx0, sx4 - dpx0)
746  psfsub = psfarr[psf_y_slice, psf_x_slice]
747  vsub = valid[sy1-ylo: sy2-ylo, sx1-xlo: sx2-xlo]
748  A[inpsf[valid], I_psf] = psfsub[vsub]
749 
750  # PSF dx -- by taking the half-difference of shifted-by-one and
751  # shifted-by-minus-one.
752  oldsx = (sx1, sx2, sx3, sx4)
753  sx1, sx2, sx3, sx4 = _overlap(xlo, xhi, px0+1, px1-1)
754  psfsub = (psfarr[psf_y_slice, sx3 - dpx0 + 1: sx4 - dpx0 + 1] -
755  psfarr[psf_y_slice, sx3 - dpx0 - 1: sx4 - dpx0 - 1])/2.
756  vsub = valid[sy1-ylo: sy2-ylo, sx1-xlo: sx2-xlo]
757  A[indx[valid], I_dx] = psfsub[vsub]
758  # revert x indices...
759  (sx1, sx2, sx3, sx4) = oldsx
760 
761  # PSF dy
762  sy1, sy2, sy3, sy4 = _overlap(ylo, yhi, py0+1, py1-1)
763  psfsub = (psfarr[sy3 - dpy0 + 1: sy4 - dpy0 + 1, psf_x_slice] -
764  psfarr[sy3 - dpy0 - 1: sy4 - dpy0 - 1, psf_x_slice])/2.
765  vsub = valid[sy1-ylo: sy2-ylo, sx1-xlo: sx2-xlo]
766  A[indy[valid], I_dy] = psfsub[vsub]
767 
768  # other PSFs...
769  for j, opsf in enumerate(otherpeaks):
770  obb = opsf.getBBox()
771  ino = np.outer((yy >= obb.getMinY())*(yy <= obb.getMaxY()),
772  (xx >= obb.getMinX())*(xx <= obb.getMaxX()))
773  dpx0, dpy0 = obb.getMinX() - xlo, obb.getMinY() - ylo
774  sx1, sx2, sx3, sx4 = _overlap(xlo, xhi, obb.getMinX(), obb.getMaxX())
775  sy1, sy2, sy3, sy4 = _overlap(ylo, yhi, obb.getMinY(), obb.getMaxY())
776  opsfarr = opsf.getArray()
777  psfsub = opsfarr[sy3 - dpy0 : sy4 - dpy0, sx3 - dpx0: sx4 - dpx0]
778  vsub = valid[sy1-ylo: sy2-ylo, sx1-xlo: sx2-xlo]
779  A[ino[valid], I_opsf + j] = psfsub[vsub]
780 
781  b = img_sub[valid]
782 
783  # Weights -- from ramp and image variance map.
784  # Ramp weights -- from 1 at R0 down to 0 at R1.
785  rw = np.ones_like(RR)
786  ii = (RR > R0**2)
787  rr = np.sqrt(RR[ii])
788  rw[ii] = np.maximum(0, 1. - ((rr - R0)/(R1 - R0)))
789  w = np.sqrt(rw[valid]/var_sub[valid])
790  # save the effective number of pixels
791  sumr = np.sum(rw[valid])
792  log.log(-9, 'sumr = %g' % sumr)
793 
794  del ii
795 
796  Aw = A*w[:, np.newaxis]
797  bw = b*w
798 
799  if debugPlots:
800  import pylab as plt
801  plt.clf()
802  N = NT2 + 2
803  R, C = 2, (N+1)/2
804  for i in range(NT2):
805  im1 = np.zeros((1+yhi-ylo, 1+xhi-xlo))
806  im1[ipixes[:, 1], ipixes[:, 0]] = A[:, i]
807  plt.subplot(R, C, i+1)
808  plt.imshow(im1, interpolation='nearest', origin='lower')
809  plt.subplot(R, C, NT2+1)
810  im1 = np.zeros((1+yhi-ylo, 1+xhi-xlo))
811  im1[ipixes[:, 1], ipixes[:, 0]] = b
812  plt.imshow(im1, interpolation='nearest', origin='lower')
813  plt.subplot(R, C, NT2+2)
814  im1 = np.zeros((1+yhi-ylo, 1+xhi-xlo))
815  im1[ipixes[:, 1], ipixes[:, 0]] = w
816  plt.imshow(im1, interpolation='nearest', origin='lower')
817  plt.savefig('A.png')
818 
819  # We do fits with and without the decenter (dx,dy) terms.
820  # Since the dx,dy terms are at the end of the matrix,
821  # we can do that just by trimming off those elements.
822  #
823  # The SVD can fail if there are NaNs in the matrices; this should
824  # really be handled upstream
825  try:
826  # NT1 is number of terms without dx,dy;
827  # X1 is the result without decenter
828  X1, r1, rank1, s1 = np.linalg.lstsq(Aw[:, :NT1], bw)
829  # X2 is with decenter
830  X2, r2, rank2, s2 = np.linalg.lstsq(Aw, bw)
831  except np.linalg.LinAlgError, e:
832  log.log(log.WARN, "Failed to fit PSF to child: %s" % e)
833  pkres.setPsfFitFailed()
834  return
835 
836  log.log(-9, 'r1 r2 %s %s' % (r1, r2))
837 
838  # r is weighted chi-squared = sum over pixels: ramp * (model -
839  # data)**2/sigma**2
840  if len(r1) > 0:
841  chisq1 = r1[0]
842  else:
843  chisq1 = 1e30
844  if len(r2) > 0:
845  chisq2 = r2[0]
846  else:
847  chisq2 = 1e30
848  dof1 = sumr - len(X1)
849  dof2 = sumr - len(X2)
850  log.log(-9, 'dof1, dof2 %g %g' % (dof1, dof2))
851 
852  # This can happen if we're very close to the edge (?)
853  if dof1 <= 0 or dof2 <= 0:
854  log.logdebug('Skipping this peak: bad DOF %g, %g' % (dof1, dof2))
855  pkres.setBadPsfDof()
856  return
857 
858  q1 = chisq1/dof1
859  q2 = chisq2/dof2
860  log.logdebug('PSF fits: chisq/dof = %g, %g' % (q1, q2))
861  ispsf1 = (q1 < psfChisqCut1)
862  ispsf2 = (q2 < psfChisqCut2)
863 
864  pkres.psfFit1 = (chisq1, dof1)
865  pkres.psfFit2 = (chisq2, dof2)
866 
867  # check that the fit PSF spatial derivative terms aren't too big
868  if ispsf2:
869  fdx, fdy = X2[I_dx], X2[I_dy]
870  f0 = X2[I_psf]
871  # as a fraction of the PSF flux
872  dx = fdx/f0
873  dy = fdy/f0
874  ispsf2 = ispsf2 and (abs(dx) < 1. and abs(dy) < 1.)
875  log.logdebug('isPSF2 -- checking derivatives: dx,dy = %g, %g -> %s' %
876  (dx, dy, str(ispsf2)))
877  if not ispsf2:
878  pkres.psfFitBigDecenter = True
879 
880  # Looks like a shifted PSF: try actually shifting the PSF by that amount
881  # and re-evaluate the fit.
882  if ispsf2:
883  psfimg2 = psf.computeImage(cx + dx, cy + dy)
884  # clip
885  pbb2 = psfimg2.getBBox()
886  pbb2.clip(fbb)
887 
888  # Make sure we haven't been given a substitute PSF that's nowhere near where we want, as may occur if
889  # "Cannot compute CoaddPsf at point (xx,yy); no input images at that point."
890  if not pbb2.contains(afwGeom.Point2I(int(cx + dx), int(cy + dy))):
891  ispsf2 = False
892  else:
893  # clip image to bbox
894  px0, py0 = psfimg2.getX0(), psfimg2.getY0()
895  psfarr = psfimg2.getArray()[pbb2.getMinY()-py0:1+pbb2.getMaxY()-py0,
896  pbb2.getMinX()-px0:1+pbb2.getMaxX()-px0]
897  px0, py0 = pbb2.getMinX(), pbb2.getMinY()
898  px1, py1 = pbb2.getMaxX(), pbb2.getMaxY()
899 
900  # yuck! Update the PSF terms in the least-squares fit matrix.
901  Ab = A[:, :NT1]
902 
903  sx1, sx2, sx3, sx4 = _overlap(xlo, xhi, px0, px1)
904  sy1, sy2, sy3, sy4 = _overlap(ylo, yhi, py0, py1)
905  dpx0, dpy0 = px0 - xlo, py0 - ylo
906  psfsub = psfarr[sy3-dpy0:sy4-dpy0, sx3-dpx0:sx4-dpx0]
907  vsub = valid[sy1-ylo:sy2-ylo, sx1-xlo:sx2-xlo]
908  xx, yy = np.arange(xlo, xhi+1), np.arange(ylo, yhi+1)
909  inpsf = np.outer((yy >= py0)*(yy <= py1), (xx >= px0)*(xx <= px1))
910  Ab[inpsf[valid], I_psf] = psfsub[vsub]
911 
912  Aw = Ab*w[:, np.newaxis]
913  # re-solve...
914  Xb, rb, rankb, sb = np.linalg.lstsq(Aw, bw)
915  if len(rb) > 0:
916  chisqb = rb[0]
917  else:
918  chisqb = 1e30
919  dofb = sumr - len(Xb)
920  qb = chisqb/dofb
921  ispsf2 = (qb < psfChisqCut2b)
922  q2 = qb
923  X2 = Xb
924  log.logdebug('shifted PSF: new chisq/dof = %g; good? %s' %
925  (qb, ispsf2))
926  pkres.psfFit3 = (chisqb, dofb)
927 
928  # Which one do we keep?
929  if (((ispsf1 and ispsf2) and (q2 < q1)) or
930  (ispsf2 and not ispsf1)):
931  Xpsf = X2
932  chisq = chisq2
933  dof = dof2
934  log.log(-9, 'dof %g' % dof)
935  log.logdebug('Keeping shifted-PSF model')
936  cx += dx
937  cy += dy
938  pkres.psfFitWithDecenter = True
939  else:
940  # (arbitrarily set to X1 when neither fits well)
941  Xpsf = X1
942  chisq = chisq1
943  dof = dof1
944  log.log(-9, 'dof %g' % dof)
945  log.logdebug('Keeping unshifted PSF model')
946 
947  ispsf = (ispsf1 or ispsf2)
948 
949  # Save the PSF models in images for posterity.
950  if debugPsf:
951  SW, SH = 1+xhi-xlo, 1+yhi-ylo
952  psfmod = afwImage.ImageF(SW, SH)
953  psfmod.setXY0(xlo, ylo)
954  psfderivmodm = afwImage.MaskedImageF(SW, SH)
955  psfderivmod = psfderivmodm.getImage()
956  psfderivmod.setXY0(xlo, ylo)
957  model = afwImage.ImageF(SW, SH)
958  model.setXY0(xlo, ylo)
959  for i in range(len(Xpsf)):
960  for (x, y),v in zip(ipixes, A[:, i]*Xpsf[i]):
961  ix, iy = int(x), int(y)
962  model.set(ix, iy, model.get(ix, iy) + float(v))
963  if i in [I_psf, I_dx, I_dy]:
964  psfderivmod.set(ix, iy, psfderivmod.get(ix, iy) + float(v))
965  for ii in range(NP):
966  x, y = ipixes[ii, :]
967  psfmod.set(int(x), int(y), float(A[ii, I_psf]*Xpsf[I_psf]))
968  modelfp = afwDet.Footprint(fp.getPeaks().getSchema())
969  for (x, y) in ipixes:
970  modelfp.addSpan(int(y+ylo), int(x+xlo), int(x+xlo))
971  modelfp.normalize()
972 
973  pkres.psfFitDebugPsf0Img = psfimg
974  pkres.psfFitDebugPsfImg = psfmod
975  pkres.psfFitDebugPsfDerivImg = psfderivmod
976  pkres.psfFitDebugPsfModel = model
977  pkres.psfFitDebugStamp = img.Factory(img, stampbb, True)
978  pkres.psfFitDebugValidPix = valid # numpy array
979  pkres.psfFitDebugVar = varimg.Factory(varimg, stampbb, True)
980  ww = np.zeros(valid.shape, np.float)
981  ww[valid] = w
982  pkres.psfFitDebugWeight = ww # numpy
983  pkres.psfFitDebugRampWeight = rw
984 
985 
986  # Save things we learned about this peak for posterity...
987  pkres.psfFitR0 = R0
988  pkres.psfFitR1 = R1
989  pkres.psfFitStampExtent = (xlo, xhi, ylo, yhi)
990  pkres.psfFitCenter = (cx, cy)
991  log.log(-9, 'saving chisq,dof %g %g' % (chisq, dof))
992  pkres.psfFitBest = (chisq, dof)
993  pkres.psfFitParams = Xpsf
994  pkres.psfFitFlux = Xpsf[I_psf]
995  pkres.psfFitNOthers = len(otherpeaks)
996 
997  if ispsf:
998  pkres.setDeblendedAsPsf()
999 
1000  # replace the template image by the PSF + derivatives
1001  # image.
1002  log.logdebug('Deblending as PSF; setting template to PSF model')
1003 
1004  # Instantiate the PSF model and clip it to the footprint
1005  psfimg = psf.computeImage(cx, cy)
1006  # Scale by fit flux.
1007  psfimg *= Xpsf[I_psf]
1008  psfimg = psfimg.convertF()
1009 
1010  # Clip the Footprint to the PSF model image bbox.
1011  fpcopy = afwDet.Footprint(fp)
1012  psfbb = psfimg.getBBox()
1013  fpcopy.clipTo(psfbb)
1014  bb = fpcopy.getBBox()
1015 
1016  # Copy the part of the PSF model within the clipped footprint.
1017  psfmod = afwImage.ImageF(bb)
1018  afwDet.copyWithinFootprintImage(fpcopy, psfimg, psfmod)
1019  # Save it as our template.
1020  fpcopy.clipToNonzero(psfmod)
1021  fpcopy.normalize()
1022  pkres.setTemplate(psfmod, fpcopy)
1023 
1024  # DEBUG
1025  pkres.setPsfTemplate(psfmod, fpcopy)
1026 
An integer coordinate rectangle.
Definition: Box.h:53
A set of pixels in an Image.
Definition: Footprint.h:62
def lsst.meas.deblender.baseline._fitPsfs (   fp,
  peaks,
  fpres,
  log,
  psf,
  psffwhm,
  img,
  varimg,
  psfChisqCut1,
  psfChisqCut2,
  psfChisqCut2b,
  kwargs 
)
private

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.

Parameters
[in]fpFootprint containing the Peaks to model
[in]peaksa list of all the Peak objects within this Footprint
[in,out]fpresa PerFootprint result object where results will be placed
[in,out]logpex Log object
[in]psfafw PSF object
[in]psffwhmFWHM of the PSF, in pixels
[in]imgthe Image in which this Footprint finds itself
[in]varimgVariance plane
[in]psfChisqCut*floats: cuts in chisq-per-pixel at which to accept the PSF model

Results go into the fpres.peaks objects.

Definition at line 539 of file baseline.py.

540  **kwargs):
541  """!
542  Fit a PSF + smooth background models to a small region around each
543  peak (plus other peaks that are nearby) in the given list of
544  peaks.
545 
546  @param[in] fp Footprint containing the Peaks to model
547  @param[in] peaks a list of all the Peak objects within this Footprint
548  @param[in,out] fpres a PerFootprint result object where results will be placed
549  @param[in,out] log pex Log object
550  @param[in] psf afw PSF object
551  @param[in] psffwhm FWHM of the PSF, in pixels
552  @param[in] img the Image in which this Footprint finds itself
553  @param[in] varimg Variance plane
554  @param[in] psfChisqCut* floats: cuts in chisq-per-pixel at which to accept
555  the PSF model
556 
557  Results go into the fpres.peaks objects.
558  """
559  fbb = fp.getBBox()
560  cpsf = CachingPsf(psf)
561 
562  # create mask image for pixels within the footprint
563  fmask = afwImage.MaskU(fbb)
564  fmask.setXY0(fbb.getMinX(), fbb.getMinY())
565  afwDet.setMaskFromFootprint(fmask, fp, 1)
566 
567  # pk.getF() -- retrieving the floating-point location of the peak
568  # -- actually shows up in the profile if we do it in the loop, so
569  # grab them all here.
570  peakF = [pk.getF() for pk in peaks]
571 
572  for pki,(pk,pkres,pkF) in enumerate(zip(peaks, fpres.peaks, peakF)):
573  log.logdebug('Peak %i' % pki)
574  _fitPsf(fp, fmask, pk, pkF, pkres, fbb, peaks, peakF, log, cpsf, psffwhm,
575  img, varimg, psfChisqCut1, psfChisqCut2, psfChisqCut2b,
576  **kwargs)
577 
def _fitPsf
Fit a PSF + smooth background model (linear) to a small region around the peak.
Definition: baseline.py:581
In the PSF fitting code, we request PSF models for all peaks near the one being fit.
Definition: baseline.py:516
MaskT setMaskFromFootprint(lsst::afw::image::Mask< MaskT > *mask, Footprint const &footprint, MaskT const bitmask)
OR bitmask into all the Mask&#39;s pixels that are in the Footprint.
def lsst.meas.deblender.baseline._handle_flux_at_edge (   log,
  psffwhm,
  t1,
  tfoot,
  fp,
  maskedImage,
  x0,
  x1,
  y0,
  y1,
  psf,
  pk,
  sigma1,
  patchEdges 
)
private

Definition at line 1028 of file baseline.py.

1029  x0, x1, y0, y1, psf, pk, sigma1, patchEdges):
1030  # Import C++ routines
1031  import lsst.meas.deblender as deb
1032  butils = deb.BaselineUtilsF
1033 
1034  log.logdebug('Found significant flux at template edge.')
1035  # Compute the max of:
1036  # -symmetric-template-clipped image * PSF
1037  # -footprint-clipped image
1038  # Ie, extend the template by the PSF and "fill in" the footprint.
1039  # Then find the symmetric template of that image.
1040 
1041  # The size we'll grow by
1042  S = psffwhm*1.5
1043  # make it an odd integer
1044  S = (int(S + 0.5)/2)*2 + 1
1045 
1046  tbb = tfoot.getBBox()
1047  tbb.grow(S)
1048 
1049  # (footprint+margin)-clipped image;
1050  # we need the pixels OUTSIDE the footprint to be 0.
1051  fpcopy = afwDet.Footprint(fp)
1052  fpcopy = afwDet.growFootprint(fpcopy, S)
1053  fpcopy.clipTo(tbb)
1054  fpcopy.normalize()
1055  padim = maskedImage.Factory(tbb)
1056  afwDet.copyWithinFootprintMaskedImage(fpcopy, maskedImage, padim)
1057 
1058  # find pixels on the edge of the template
1059  edgepix = butils.getSignificantEdgePixels(t1, tfoot, -1e6)
1060 
1061  # instantiate PSF image
1062  xc = int((x0 + x1)/2)
1063  yc = int((y0 + y1)/2)
1064  psfim = psf.computeImage(afwGeom.Point2D(xc, yc))
1065  pbb = psfim.getBBox()
1066  # shift PSF image to by centered on zero
1067  lx, ly = pbb.getMinX(), pbb.getMinY()
1068  psfim.setXY0(lx - xc, ly - yc)
1069  pbb = psfim.getBBox()
1070  # clip PSF to S, if necessary
1071  Sbox = afwGeom.Box2I(afwGeom.Point2I(-S, -S), afwGeom.Extent2I(2*S+1, 2*S+1))
1072  if not Sbox.contains(pbb):
1073  # clip PSF image
1074  psfim = psfim.Factory(psfim, Sbox, afwImage.PARENT, True)
1075  pbb = psfim.getBBox()
1076  px0 = pbb.getMinX()
1077  px1 = pbb.getMaxX()
1078  py0 = pbb.getMinY()
1079  py1 = pbb.getMaxY()
1080 
1081  # Compute the ramped-down edge pixels
1082  ramped = t1.Factory(tbb)
1083  Tout = ramped.getArray()
1084  Tin = t1.getArray()
1085  tx0, ty0 = t1.getX0(), t1.getY0()
1086  ox0, oy0 = ramped.getX0(), ramped.getY0()
1087  P = psfim.getArray()
1088  P /= P.max()
1089  # For each edge pixel, Tout = max(Tout, edgepix * PSF)
1090  for span in edgepix.getSpans():
1091  y = span.getY()
1092  for x in range(span.getX0(), span.getX1()+1):
1093  slc = (slice(y+py0 - oy0, y+py1+1 - oy0),
1094  slice(x+px0 - ox0, x+px1+1 - ox0))
1095  Tout[slc] = np.maximum(Tout[slc], Tin[y-ty0, x-tx0]*P)
1096 
1097  # Fill in the "padim" (which has the right variance and
1098  # mask planes) with the ramped pixels, outside the footprint
1099  I = (padim.getImage().getArray() == 0)
1100  padim.getImage().getArray()[I] = ramped.getArray()[I]
1101 
1102  rtn, patched = butils.buildSymmetricTemplate(padim, fpcopy, pk, sigma1, True, patchEdges)
1103  # silly SWIG can't unpack pairs as tuples
1104  t2, tfoot2 = rtn[0], rtn[1]
1105  del rtn
1106 
1107  # This template footprint may extend outside the parent
1108  # footprint -- or the image. Clip it.
1109  # NOTE that this may make it asymmetric, unlike normal templates.
1110  imbb = maskedImage.getBBox()
1111  tfoot2.clipTo(imbb)
1112  tbb = tfoot2.getBBox()
1113  # clip template image to bbox
1114  t2 = t2.Factory(t2, tbb, afwImage.PARENT, True)
1115 
1116  return t2, tfoot2, patched
boost::shared_ptr< Footprint > growFootprint(Footprint const &foot, int nGrow, bool left, bool right, bool up, bool down)
Grow a Footprint in at least one of the cardinal directions, returning a new Footprint.
An integer coordinate rectangle.
Definition: Box.h:53
A set of pixels in an Image.
Definition: Footprint.h:62
def lsst.meas.deblender.baseline.deblend (   footprint,
  maskedImage,
  psf,
  psffwhm,
  psfChisqCut1 = 1.5,
  psfChisqCut2 = 1.5,
  psfChisqCut2b = 1.5,
  fitPsfs = True,
  medianSmoothTemplate = True,
  medianFilterHalfsize = 2,
  monotonicTemplate = True,
  weightTemplates = False,
  log = None,
  verbose = False,
  sigma1 = None,
  maxNumberOfPeaks = 0,
  findStrayFlux = True,
  assignStrayFlux = True,
  strayFluxToPointSources = 'necessary',
  strayFluxAssignment = 'r-to-peak',
  rampFluxAtEdge = False,
  patchEdges = False,
  tinyFootprintSize = 2,
  getTemplateSum = False,
  clipStrayFluxFraction = 0.001,
  clipFootprintToNonzero = True 
)

Deblend a single footprint in a maskedImage.

Each Peak in the Footprint will produce a deblended child.

psfChisqCut1, psfChisqCut2: used when deciding whether a given peak looks like a point source. We build a small local model of background + peak + neighboring peaks. These are the chi-squared-per-degree-of-freedom cuts (1=without and 2=with terms allowing the peak position to move); larger values will result in more peaks being classified as PSFs.

psfChisqCut2b: if the with-decenter fit is accepted, we then apply the computed dx,dy to the source position and re-fit without decentering. The model is accepted if its chi-squared-per-DOF is less than this value.

maxNumberOfPeaks: if positive, only deblend the brightest maxNumberOfPeaks peaks in the parent

Definition at line 233 of file baseline.py.

234  ):
235  """!
236  Deblend a single ``footprint`` in a ``maskedImage``.
237 
238  Each ``Peak`` in the ``Footprint`` will produce a deblended child.
239 
240  psfChisqCut1, psfChisqCut2: used when deciding whether a given
241  peak looks like a point source. We build a small local model of
242  background + peak + neighboring peaks. These are the
243  chi-squared-per-degree-of-freedom cuts (1=without and 2=with terms
244  allowing the peak position to move); larger values will result in
245  more peaks being classified as PSFs.
246 
247  psfChisqCut2b: if the with-decenter fit is accepted, we then
248  apply the computed dx,dy to the source position and re-fit without
249  decentering. The model is accepted if its chi-squared-per-DOF is
250  less than this value.
251 
252  maxNumberOfPeaks: if positive, only deblend the brightest
253  maxNumberOfPeaks peaks in the parent
254  """
255  # Import C++ routines
256  import lsst.meas.deblender as deb
257  butils = deb.BaselineUtilsF
258 
259  validStrayPtSrc = ['never', 'necessary', 'always']
260  validStrayAssign = ['r-to-peak', 'r-to-footprint', 'nearest-footprint', 'trim']
261  if not strayFluxToPointSources in validStrayPtSrc:
262  raise ValueError((('strayFluxToPointSources: value \"%s\" not in the '
263  + 'set of allowed values: ')
264  % strayFluxToPointSources) + str(validStrayPtSrc))
265  if not strayFluxAssignment in validStrayAssign:
266  raise ValueError((('strayFluxAssignment: value \"%s\" not in the '
267  + 'set of allowed values: ')
268  % strayFluxAssignment) + str(validStrayAssign))
269 
270  if log is None:
271  import lsst.pex.logging as pexLogging
272 
273  component = 'meas_deblender.baseline'
274  log = pexLogging.Log(pexLogging.Log.getDefaultLog(), component)
275 
276  if verbose: # pexLogging defines DEBUG in terms of "importance" == -verbosity
277  pexLogging.Trace.setVerbosity(component, -pexLogging.Log.DEBUG)
278 
279  img = maskedImage.getImage()
280  varimg = maskedImage.getVariance()
281  mask = maskedImage.getMask()
282 
283  fp = footprint
284  fp.normalize()
285  peaks = fp.getPeaks()
286  if maxNumberOfPeaks > 0:
287  peaks = peaks[:maxNumberOfPeaks]
288 
289  # Pull out the image bounds of the parent Footprint
290  imbb = img.getBBox()
291  bb = fp.getBBox()
292  if not imbb.contains(bb):
293  raise ValueError(('Footprint bounding-box %s extends outside image '
294  + 'bounding-box %s') % (str(bb), str(imbb)))
295  W, H = bb.getWidth(), bb.getHeight()
296  x0, y0 = bb.getMinX(), bb.getMinY()
297  x1, y1 = bb.getMaxX(), bb.getMaxY()
298 
299  # 'sigma1' is an estimate of the average noise level in this image.
300  if sigma1 is None:
301  # FIXME -- just within the bbox?
302  stats = afwMath.makeStatistics(varimg, mask, afwMath.MEDIAN)
303  sigma1 = math.sqrt(stats.getValue(afwMath.MEDIAN))
304  log.logdebug('Estimated sigma1 = %f' % sigma1)
305 
306  # get object that will hold our results
307  res = PerFootprint(fp, peaks=peaks)
308 
309  if fitPsfs:
310  # Find peaks that are well-fit by a PSF + background model.
311  _fitPsfs(fp, peaks, res, log, psf, psffwhm, img, varimg,
312  psfChisqCut1, psfChisqCut2, psfChisqCut2b,
313  tinyFootprintSize=tinyFootprintSize)
314 
315  # Create templates...
316  log.logdebug(('Creating templates for footprint at x0,y0,W,H = ' +
317  '(%i, %i, %i, %i)') % (x0, y0, W, H))
318  for peaki, pkres in enumerate(res.peaks):
319  log.logdebug('Deblending peak %i of %i' % (peaki, len(res.peaks)))
320  if pkres.skip or pkres.deblendedAsPsf:
321  continue
322  pk = pkres.peak
323  cx, cy = pk.getIx(), pk.getIy()
324  if not imbb.contains(afwGeom.Point2I(cx,cy)):
325  log.logdebug('Peak center is not inside image; skipping %i' % pkres.pki)
326  pkres.setOutOfBounds()
327  continue
328  log.logdebug('computing template for peak %i at (%i, %i)' % (pkres.pki, cx, cy))
329  S, patched = butils.buildSymmetricTemplate(maskedImage, fp, pk, sigma1, True, patchEdges)
330  # SWIG doesn't know how to unpack a std::pair into a 2-tuple...
331  # (in this case, a nested pair)
332  t1, tfoot = S[0], S[1]
333  del S
334 
335  if t1 is None:
336  log.logdebug(('Peak %i at (%i, %i): failed to build symmetric ' +
337  'template') % (pkres.pki, cx, cy))
338  pkres.setFailedSymmetricTemplate()
339  continue
340 
341  if patched:
342  pkres.setPatched()
343 
344  # possibly save the original symmetric template
345  pkres.setOrigTemplate(t1, tfoot)
346 
347  if rampFluxAtEdge:
348  log.logdebug('Checking for significant flux at edge: sigma1=%g' % sigma1)
349  if (rampFluxAtEdge and
350  butils.hasSignificantFluxAtEdge(t1, tfoot, 3*sigma1)):
351  log.logdebug("Template %i has significant flux at edge: ramping" % pkres.pki)
352  try:
353  (t2, tfoot2, patched) = _handle_flux_at_edge(log, psffwhm, t1, tfoot, fp,
354  maskedImage, x0, x1, y0, y1,
355  psf, pk, sigma1, patchEdges)
356  except lsst.pex.exceptions.Exception as exc:
357  if (isinstance(exc, lsst.pex.exceptions.InvalidParameterError) and "CoaddPsf" in str(exc)):
358  pkres.setOutOfBounds()
359  continue
360  raise
361  pkres.setRampedTemplate(t2, tfoot2)
362  if patched:
363  pkres.setPatched()
364  t1 = t2
365  tfoot = tfoot2
366 
367  if medianSmoothTemplate:
368  filtsize = medianFilterHalfsize*2 + 1
369  if t1.getWidth() >= filtsize and t1.getHeight() >= filtsize:
370  log.logdebug('Median filtering template %i' % pkres.pki)
371  # We want the output to go in "t1", so copy it into
372  # "inimg" for input
373  inimg = t1.Factory(t1, True)
374  butils.medianFilter(inimg, t1, medianFilterHalfsize)
375  # possible save this median-filtered template
376  pkres.setMedianFilteredTemplate(t1, tfoot)
377  else:
378  log.logdebug(('Not median-filtering template %i: size '
379  + '%i x %i smaller than required %i x %i') %
380  (pkres.pki, t1.getWidth(), t1.getHeight(), filtsize, filtsize))
381 
382  if monotonicTemplate:
383  log.logdebug('Making template %i monotonic' % pkres.pki)
384  butils.makeMonotonic(t1, pk)
385 
386  if clipFootprintToNonzero:
387  tfoot.clipToNonzero(t1)
388  tfoot.normalize()
389  if not tfoot.getBBox().isEmpty() and tfoot.getBBox() != t1.getBBox(afwImage.PARENT):
390  t1 = t1.Factory(t1, tfoot.getBBox(), afwImage.PARENT, True)
391 
392  pkres.setTemplate(t1, tfoot)
393 
394  # Prepare inputs to "apportionFlux" call.
395  # template maskedImages
396  tmimgs = []
397  # template footprints
398  tfoots = []
399  # deblended as psf
400  dpsf = []
401  # peak x,y
402  pkx = []
403  pky = []
404  # indices of valid templates
405  ibi = []
406 
407  for peaki, pkres in enumerate(res.peaks):
408  if pkres.skip:
409  continue
410  tmimgs.append(pkres.templateImage)
411  tfoots.append(pkres.templateFootprint)
412  # for stray flux...
413  dpsf.append(pkres.deblendedAsPsf)
414  pk = pkres.peak
415  pkx.append(pk.getIx())
416  pky.append(pk.getIy())
417  ibi.append(pkres.pki)
418 
419  if weightTemplates:
420  # Reweight the templates by doing a least-squares fit to the image
421  A = np.zeros((W*H, len(tmimgs)))
422  b = img.getArray()[y0:y1+1, x0:x1+1].ravel()
423 
424  for mim, i in zip(tmimgs, ibi):
425  ibb = mim.getBBox()
426  ix0, iy0 = ibb.getMinX(), ibb.getMinY()
427  pkres = res.peaks[i]
428  foot = pkres.templateFootprint
429  ima = mim.getImage().getArray()
430  for sp in foot.getSpans():
431  sy, sx0, sx1 = sp.getY(), sp.getX0(), sp.getX1()
432  imrow = ima[sy-iy0, sx0-ix0 : 1+sx1-ix0]
433  r0 = (sy-y0)*W
434  A[r0 + sx0-x0: r0+1+sx1-x0, i] = imrow
435 
436  X1, r1, rank1, s1 = np.linalg.lstsq(A, b)
437  del A
438  del b
439 
440  for mim, i, w in zip(tmimgs, ibi, X1):
441  im = mim.getImage()
442  im *= w
443  res.peaks[i].setTemplateWeight(w)
444 
445  # FIXME -- Remove templates that are too similar (via dot-product test)?
446 
447  # Now apportion flux according to the templates
448  log.logdebug('Apportioning flux among %i templates' % len(tmimgs))
449  sumimg = afwImage.ImageF(bb)
450  #.getDimensions())
451  #sumimg.setXY0(bb.getMinX(), bb.getMinY())
452 
453  strayflux = afwDet.HeavyFootprintPtrListF()
454 
455  strayopts = 0
456  if strayFluxAssignment == 'trim':
457  findStrayFlux = False
458  strayopts |= butils.STRAYFLUX_TRIM
459  if findStrayFlux:
460  strayopts |= butils.ASSIGN_STRAYFLUX
461  if strayFluxToPointSources == 'necessary':
462  strayopts |= butils.STRAYFLUX_TO_POINT_SOURCES_WHEN_NECESSARY
463  elif strayFluxToPointSources == 'always':
464  strayopts |= butils.STRAYFLUX_TO_POINT_SOURCES_ALWAYS
465 
466  if strayFluxAssignment == 'r-to-peak':
467  # this is the default
468  pass
469  elif strayFluxAssignment == 'r-to-footprint':
470  strayopts |= butils.STRAYFLUX_R_TO_FOOTPRINT
471  elif strayFluxAssignment == 'nearest-footprint':
472  strayopts |= butils.STRAYFLUX_NEAREST_FOOTPRINT
473 
474  portions = butils.apportionFlux(maskedImage, fp, tmimgs, tfoots, sumimg, dpsf,
475  pkx, pky, strayflux, strayopts, clipStrayFluxFraction)
476 
477  # Shrink parent to union of children
478  if strayFluxAssignment == 'trim':
479  fp.include(tfoots, True)
480 
481  if getTemplateSum:
482  res.setTemplateSum(sumimg)
483 
484  # Save the apportioned fluxes
485  ii = 0
486  for j, (pk, pkres) in enumerate(zip(peaks, res.peaks)):
487  if pkres.skip:
488  continue
489  pkres.setFluxPortion(portions[ii])
490 
491  if findStrayFlux:
492  # NOTE that due to a swig bug (https://github.com/swig/swig/issues/59)
493  # we CANNOT iterate over "strayflux", but must index into it.
494  stray = strayflux[ii]
495  else:
496  stray = None
497  ii += 1
498 
499  pkres.setStrayFlux(stray)
500 
501  # Set child footprints to contain the right number of peaks.
502  for j, (pk, pkres) in enumerate(zip(peaks, res.peaks)):
503  if pkres.skip:
504  continue
505 
506  for foot, add in [(pkres.templateFootprint, True), (pkres.origFootprint, True),
507  (pkres.strayFlux, False)]:
508  if foot is None:
509  continue
510  pks = foot.getPeaks()
511  pks.clear()
512  if add:
513  pks.append(pk)
514 
515  return res
a place to record messages and descriptions of the state of processing.
Definition: Log.h:154
def _fitPsfs
Fit a PSF + smooth background models to a small region around each peak (plus other peaks that are ne...
Definition: baseline.py:539
Statistics makeStatistics(afwImage::Mask< afwImage::MaskPixel > const &msk, int const flags, StatisticsControl const &sctrl)
Specialization to handle Masks.
Definition: Statistics.cc:1082