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
baseline.py
Go to the documentation of this file.
1 #!/usr/bin/env python
2 #
3 # LSST Data Management System
4 #
5 # Copyright 2008-2015 AURA/LSST.
6 #
7 # This product includes software developed by the
8 # LSST Project (http://www.lsst.org/).
9 #
10 # This program is free software: you can redistribute it and/or modify
11 # it under the terms of the GNU General Public License as published by
12 # the Free Software Foundation, either version 3 of the License, or
13 # (at your option) any later version.
14 #
15 # This program is distributed in the hope that it will be useful,
16 # but WITHOUT ANY WARRANTY; without even the implied warranty of
17 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 # GNU General Public License for more details.
19 #
20 # You should have received a copy of the LSST License Statement and
21 # the GNU General Public License along with this program. If not,
22 # see <https://www.lsstcorp.org/LegalNotices/>.
23 #
24 import math
25 import numpy as np
26 
28 import lsst.afw.image as afwImage
29 import lsst.afw.detection as afwDet
30 import lsst.afw.geom as afwGeom
31 import lsst.afw.math as afwMath
32 
33 # Result objects; will probably go away as we push more results
34 # into earlier-created Source objects, but useful for now for
35 # hauling debugging results around.
36 
37 class PerFootprint(object):
38  """Result of deblending a single parent Footprint."""
39 
40  def __init__(self, fp, peaks=None):
41  """!
42  Creates a result object for the given parent Footprint 'fp'.
43 
44  @param[in] fp Footprint
45  @param[in,out] peaks If specified, the list of peaks to use; default is
46  fp.getPeaks().
47  """
48  # PerPeak objects, one per Peak
49  self.peaks = []
50  if peaks is None:
51  peaks = fp.getPeaks()
52  for pki, pk in enumerate(peaks):
53  pkres = PerPeak()
54  pkres.peak = pk
55  pkres.pki = pki
56  self.peaks.append(pkres)
57 
58  self.templateSum = None
59 
60  def setTemplateSum(self, tsum):
61  self.templateSum = tsum
62 
63 class PerPeak(object):
64  """!
65  Result of deblending a single Peak within a parent Footprint.
66 
67  There is one of these objects for each Peak in the Footprint.
68  """
69  def __init__(self):
70  # Peak object
71  self.peak = None
72  # int, peak index number
73  self.pki = None
74  # union of all the ways of failing...
75  self.skip = False
76 
77  self.outOfBounds = False
78  self.tinyFootprint = False
79  self.noValidPixels = False
80  self.deblendedAsPsf = False
81 
82  # Field set during _fitPsf:
83  self.psfFitFailed = False
84  self.psfFitBadDof = False
85  # (chisq, dof) for PSF fit without decenter
86  self.psfFit1 = None
87  # (chisq, dof) for PSF fit with decenter
88  self.psfFit2 = None
89  # (chisq, dof) for PSF fit after applying decenter
90  self.psfFit3 = None
91  # decentered PSF fit wanted to move the center too much
92  self.psfFitBigDecenter = False
93  # was the fit with decenter better?
94  self.psfFitWithDecenter = False
95  #
96  self.psfFitR0 = None
97  self.psfFitR1 = None
98  self.psfFitStampExtent = None
99  self.psfFitCenter = None
100  self.psfFitBest = None
101  self.psfFitParams = None
102  self.psfFitFlux = None
103  self.psfFitNOthers = None
104 
105  # Things only set in _fitPsf when debugging is turned on:
106  self.psfFitDebugPsf0Img = None
107  self.psfFitDebugPsfImg = None
110 
112 
113  # The actual template Image and Footprint
114  self.templateImage = None
115  self.templateFootprint = None
116 
117  # The flux assigned to this template -- a MaskedImage
118  self.fluxPortion = None
119 
120  # The stray flux assigned to this template (may be None), a HeavyFootprint
121  self.strayFlux = None
122 
123  self.hasRampedTemplate = False
124 
125  self.patched = False
126 
127  # debug -- a copy of the original symmetric template
128  self.origTemplate = None
129  self.origFootprint = None
130  # MaskedImage
131  self.rampedTemplate = None
132  # MaskedImage
134 
135  # when least-squares fitting templates, the template weight.
136  self.templateWeight = 1.0
137 
138  def __str__(self):
139  return (('Per-peak deblend result: outOfBounds: %s, deblendedAsPsf: %s') %
140  (self.outOfBounds, self.deblendedAsPsf))
141 
142  @property
143  def psfFitChisq(self):
144  chisq, dof = self.psfFitBest
145  return chisq
146  @property
147  def psfFitDof(self):
148  chisq, dof = self.psfFitBest
149  return dof
150 
151  def getFluxPortion(self, strayFlux=True):
152  """!
153  Return a HeavyFootprint containing the flux apportioned to this peak.
154 
155  @param[in] strayFlux include stray flux also?
156  """
157  if self.templateFootprint is None or self.fluxPortion is None:
158  return None
160  self.fluxPortion)
161  if strayFlux:
162  if self.strayFlux is not None:
163  heavy.normalize()
164  self.strayFlux.normalize()
165  heavy = afwDet.mergeHeavyFootprintsF(heavy, self.strayFlux)
166 
167  return heavy
168 
169  def setStrayFlux(self, stray):
170  self.strayFlux = stray
171 
172  def setFluxPortion(self, mimg):
173  self.fluxPortion = mimg
174 
175  def setTemplateWeight(self, w):
176  self.templateWeight = w
177 
178  def setPatched(self):
179  self.patched = True
180 
181  # DEBUG
182  def setOrigTemplate(self, t, tfoot):
183  self.origTemplate = t.Factory(t, True)
184  self.origFootprint = tfoot
185  def setRampedTemplate(self, t, tfoot):
186  self.hasRampedTemplate = True
187  self.rampedTemplate = t.Factory(t, True)
188  def setMedianFilteredTemplate(self, t, tfoot):
189  self.medianFilteredTemplate = t.Factory(t, True)
190  def setPsfTemplate(self, tim, tfoot):
192  self.psfTemplate = tim.Factory(tim, True)
193 
194  def setOutOfBounds(self):
195  self.outOfBounds = True
196  self.skip = True
197 
198  def setTinyFootprint(self):
199  self.tinyFootprint = True
200  self.skip = True
201 
202  def setNoValidPixels(self):
203  self.noValidPixels = True
204  self.skip = True
205 
206  def setPsfFitFailed(self):
207  self.psfFitFailed = True
208 
209  def setBadPsfDof(self):
210  self.psfFitBadDof = True
211 
212  def setDeblendedAsPsf(self):
213  self.deblendedAsPsf = True
214 
216  self.failedSymmetricTemplate = True
217  self.skip = True
218 
219  def setTemplate(self, image, footprint):
220  self.templateImage = image
221  self.templateFootprint = footprint
222 
223 def deblend(footprint, maskedImage, psf, psffwhm,
224  psfChisqCut1=1.5, psfChisqCut2=1.5, psfChisqCut2b=1.5, fitPsfs=True,
225  medianSmoothTemplate=True, medianFilterHalfsize=2,
226  monotonicTemplate=True, weightTemplates=False,
227  log=None, verbose=False, sigma1=None, maxNumberOfPeaks=0,
228  findStrayFlux=True, assignStrayFlux=True,
229  strayFluxToPointSources='necessary', strayFluxAssignment='r-to-peak',
230  rampFluxAtEdge=False, patchEdges=False, tinyFootprintSize=2,
231  getTemplateSum=False,
232  clipStrayFluxFraction=0.001, clipFootprintToNonzero=True,
233  ):
234  """!
235  Deblend a single ``footprint`` in a ``maskedImage``.
236 
237  Each ``Peak`` in the ``Footprint`` will produce a deblended child.
238 
239  psfChisqCut1, psfChisqCut2: used when deciding whether a given
240  peak looks like a point source. We build a small local model of
241  background + peak + neighboring peaks. These are the
242  chi-squared-per-degree-of-freedom cuts (1=without and 2=with terms
243  allowing the peak position to move); larger values will result in
244  more peaks being classified as PSFs.
245 
246  psfChisqCut2b: if the with-decenter fit is accepted, we then
247  apply the computed dx,dy to the source position and re-fit without
248  decentering. The model is accepted if its chi-squared-per-DOF is
249  less than this value.
250 
251  maxNumberOfPeaks: if positive, only deblend the brightest
252  maxNumberOfPeaks peaks in the parent
253  """
254  # Import C++ routines
255  import lsst.meas.deblender as deb
256  butils = deb.BaselineUtilsF
257 
258  validStrayPtSrc = ['never', 'necessary', 'always']
259  validStrayAssign = ['r-to-peak', 'r-to-footprint', 'nearest-footprint', 'trim']
260  if not strayFluxToPointSources in validStrayPtSrc:
261  raise ValueError((('strayFluxToPointSources: value \"%s\" not in the '
262  + 'set of allowed values: ')
263  % strayFluxToPointSources) + str(validStrayPtSrc))
264  if not strayFluxAssignment in validStrayAssign:
265  raise ValueError((('strayFluxAssignment: value \"%s\" not in the '
266  + 'set of allowed values: ')
267  % strayFluxAssignment) + str(validStrayAssign))
268 
269  if log is None:
270  import lsst.pex.logging as pexLogging
271 
272  component = 'meas_deblender.baseline'
273  log = pexLogging.Log(pexLogging.Log.getDefaultLog(), component)
274 
275  if verbose: # pexLogging defines DEBUG in terms of "importance" == -verbosity
276  pexLogging.Trace.setVerbosity(component, -pexLogging.Log.DEBUG)
277 
278  img = maskedImage.getImage()
279  varimg = maskedImage.getVariance()
280  mask = maskedImage.getMask()
281 
282  fp = footprint
283  fp.normalize()
284  peaks = fp.getPeaks()
285  if maxNumberOfPeaks > 0:
286  peaks = peaks[:maxNumberOfPeaks]
287 
288  # Pull out the image bounds of the parent Footprint
289  imbb = img.getBBox()
290  bb = fp.getBBox()
291  if not imbb.contains(bb):
292  raise ValueError(('Footprint bounding-box %s extends outside image '
293  + 'bounding-box %s') % (str(bb), str(imbb)))
294  W, H = bb.getWidth(), bb.getHeight()
295  x0, y0 = bb.getMinX(), bb.getMinY()
296  x1, y1 = bb.getMaxX(), bb.getMaxY()
297 
298  # 'sigma1' is an estimate of the average noise level in this image.
299  if sigma1 is None:
300  # FIXME -- just within the bbox?
301  stats = afwMath.makeStatistics(varimg, mask, afwMath.MEDIAN)
302  sigma1 = math.sqrt(stats.getValue(afwMath.MEDIAN))
303  log.logdebug('Estimated sigma1 = %f' % sigma1)
304 
305  # get object that will hold our results
306  res = PerFootprint(fp, peaks=peaks)
307 
308  if fitPsfs:
309  # Find peaks that are well-fit by a PSF + background model.
310  _fitPsfs(fp, peaks, res, log, psf, psffwhm, img, varimg,
311  psfChisqCut1, psfChisqCut2, psfChisqCut2b,
312  tinyFootprintSize=tinyFootprintSize)
313 
314  # Create templates...
315  log.logdebug(('Creating templates for footprint at x0,y0,W,H = ' +
316  '(%i, %i, %i, %i)') % (x0, y0, W, H))
317  for peaki, pkres in enumerate(res.peaks):
318  log.logdebug('Deblending peak %i of %i' % (peaki, len(res.peaks)))
319  if pkres.skip or pkres.deblendedAsPsf:
320  continue
321  pk = pkres.peak
322  cx, cy = pk.getIx(), pk.getIy()
323  if not imbb.contains(afwGeom.Point2I(cx,cy)):
324  log.logdebug('Peak center is not inside image; skipping %i' % pkres.pki)
325  pkres.setOutOfBounds()
326  continue
327  log.logdebug('computing template for peak %i at (%i, %i)' % (pkres.pki, cx, cy))
328  S, patched = butils.buildSymmetricTemplate(maskedImage, fp, pk, sigma1, True, patchEdges)
329  # SWIG doesn't know how to unpack a std::pair into a 2-tuple...
330  # (in this case, a nested pair)
331  t1, tfoot = S[0], S[1]
332  del S
333 
334  if t1 is None:
335  log.logdebug(('Peak %i at (%i, %i): failed to build symmetric ' +
336  'template') % (pkres.pki, cx, cy))
337  pkres.setFailedSymmetricTemplate()
338  continue
339 
340  if patched:
341  pkres.setPatched()
342 
343  # possibly save the original symmetric template
344  pkres.setOrigTemplate(t1, tfoot)
345 
346  if rampFluxAtEdge:
347  log.logdebug('Checking for significant flux at edge: sigma1=%g' % sigma1)
348  if (rampFluxAtEdge and
349  butils.hasSignificantFluxAtEdge(t1, tfoot, 3*sigma1)):
350  log.logdebug("Template %i has significant flux at edge: ramping" % pkres.pki)
351  try:
352  (t2, tfoot2, patched) = _handle_flux_at_edge(log, psffwhm, t1, tfoot, fp,
353  maskedImage, x0, x1, y0, y1,
354  psf, pk, sigma1, patchEdges)
355  except lsst.pex.exceptions.Exception as exc:
356  if (isinstance(exc, lsst.pex.exceptions.InvalidParameterError) and "CoaddPsf" in str(exc)):
357  pkres.setOutOfBounds()
358  continue
359  raise
360  pkres.setRampedTemplate(t2, tfoot2)
361  if patched:
362  pkres.setPatched()
363  t1 = t2
364  tfoot = tfoot2
365 
366  if medianSmoothTemplate:
367  filtsize = medianFilterHalfsize*2 + 1
368  if t1.getWidth() >= filtsize and t1.getHeight() >= filtsize:
369  log.logdebug('Median filtering template %i' % pkres.pki)
370  # We want the output to go in "t1", so copy it into
371  # "inimg" for input
372  inimg = t1.Factory(t1, True)
373  butils.medianFilter(inimg, t1, medianFilterHalfsize)
374  # possible save this median-filtered template
375  pkres.setMedianFilteredTemplate(t1, tfoot)
376  else:
377  log.logdebug(('Not median-filtering template %i: size '
378  + '%i x %i smaller than required %i x %i') %
379  (pkres.pki, t1.getWidth(), t1.getHeight(), filtsize, filtsize))
380 
381  if monotonicTemplate:
382  log.logdebug('Making template %i monotonic' % pkres.pki)
383  butils.makeMonotonic(t1, pk)
384 
385  if clipFootprintToNonzero:
386  tfoot.clipToNonzero(t1)
387  tfoot.normalize()
388  if not tfoot.getBBox().isEmpty() and tfoot.getBBox() != t1.getBBox(afwImage.PARENT):
389  t1 = t1.Factory(t1, tfoot.getBBox(), afwImage.PARENT, True)
390 
391  pkres.setTemplate(t1, tfoot)
392 
393  # Prepare inputs to "apportionFlux" call.
394  # template maskedImages
395  tmimgs = []
396  # template footprints
397  tfoots = []
398  # deblended as psf
399  dpsf = []
400  # peak x,y
401  pkx = []
402  pky = []
403  # indices of valid templates
404  ibi = []
405 
406  for peaki, pkres in enumerate(res.peaks):
407  if pkres.skip:
408  continue
409  tmimgs.append(pkres.templateImage)
410  tfoots.append(pkres.templateFootprint)
411  # for stray flux...
412  dpsf.append(pkres.deblendedAsPsf)
413  pk = pkres.peak
414  pkx.append(pk.getIx())
415  pky.append(pk.getIy())
416  ibi.append(pkres.pki)
417 
418  if weightTemplates:
419  # Reweight the templates by doing a least-squares fit to the image
420  A = np.zeros((W*H, len(tmimgs)))
421  b = img.getArray()[y0:y1+1, x0:x1+1].ravel()
422 
423  for mim, i in zip(tmimgs, ibi):
424  ibb = mim.getBBox()
425  ix0, iy0 = ibb.getMinX(), ibb.getMinY()
426  pkres = res.peaks[i]
427  foot = pkres.templateFootprint
428  ima = mim.getImage().getArray()
429  for sp in foot.getSpans():
430  sy, sx0, sx1 = sp.getY(), sp.getX0(), sp.getX1()
431  imrow = ima[sy-iy0, sx0-ix0 : 1+sx1-ix0]
432  r0 = (sy-y0)*W
433  A[r0 + sx0-x0: r0+1+sx1-x0, i] = imrow
434 
435  X1, r1, rank1, s1 = np.linalg.lstsq(A, b)
436  del A
437  del b
438 
439  for mim, i, w in zip(tmimgs, ibi, X1):
440  im = mim.getImage()
441  im *= w
442  res.peaks[i].setTemplateWeight(w)
443 
444  # FIXME -- Remove templates that are too similar (via dot-product test)?
445 
446  # Now apportion flux according to the templates
447  log.logdebug('Apportioning flux among %i templates' % len(tmimgs))
448  sumimg = afwImage.ImageF(bb)
449  #.getDimensions())
450  #sumimg.setXY0(bb.getMinX(), bb.getMinY())
451 
452  strayflux = afwDet.HeavyFootprintPtrListF()
453 
454  strayopts = 0
455  if strayFluxAssignment == 'trim':
456  findStrayFlux = False
457  strayopts |= butils.STRAYFLUX_TRIM
458  if findStrayFlux:
459  strayopts |= butils.ASSIGN_STRAYFLUX
460  if strayFluxToPointSources == 'necessary':
461  strayopts |= butils.STRAYFLUX_TO_POINT_SOURCES_WHEN_NECESSARY
462  elif strayFluxToPointSources == 'always':
463  strayopts |= butils.STRAYFLUX_TO_POINT_SOURCES_ALWAYS
464 
465  if strayFluxAssignment == 'r-to-peak':
466  # this is the default
467  pass
468  elif strayFluxAssignment == 'r-to-footprint':
469  strayopts |= butils.STRAYFLUX_R_TO_FOOTPRINT
470  elif strayFluxAssignment == 'nearest-footprint':
471  strayopts |= butils.STRAYFLUX_NEAREST_FOOTPRINT
472 
473  portions = butils.apportionFlux(maskedImage, fp, tmimgs, tfoots, sumimg, dpsf,
474  pkx, pky, strayflux, strayopts, clipStrayFluxFraction)
475 
476  # Shrink parent to union of children
477  if strayFluxAssignment == 'trim':
478  fp.include(tfoots, True)
479 
480  if getTemplateSum:
481  res.setTemplateSum(sumimg)
482 
483  # Save the apportioned fluxes
484  ii = 0
485  for j, (pk, pkres) in enumerate(zip(peaks, res.peaks)):
486  if pkres.skip:
487  continue
488  pkres.setFluxPortion(portions[ii])
489 
490  if findStrayFlux:
491  # NOTE that due to a swig bug (https://github.com/swig/swig/issues/59)
492  # we CANNOT iterate over "strayflux", but must index into it.
493  stray = strayflux[ii]
494  else:
495  stray = None
496  ii += 1
497 
498  pkres.setStrayFlux(stray)
499 
500  # Set child footprints to contain the right number of peaks.
501  for j, (pk, pkres) in enumerate(zip(peaks, res.peaks)):
502  if pkres.skip:
503  continue
504 
505  for foot, add in [(pkres.templateFootprint, True), (pkres.origFootprint, True),
506  (pkres.strayFlux, False)]:
507  if foot is None:
508  continue
509  pks = foot.getPeaks()
510  pks.clear()
511  if add:
512  pks.append(pk)
513 
514  return res
515 
516 class CachingPsf(object):
517  """!
518  In the PSF fitting code, we request PSF models for all peaks near
519  the one being fit. This was turning out to be quite expensive in
520  some cases. Here, we cache the PSF models to bring the cost down
521  closer to O(N) rather than O(N^2).
522  """
523  def __init__(self, psf):
524  self.cache = {}
525  self.psf = psf
526  def computeImage(self, cx, cy):
527  im = self.cache.get((cx, cy), None)
528  if im is not None:
529  return im
530  try:
531  im = self.psf.computeImage(afwGeom.Point2D(cx, cy))
533  im = self.psf.computeImage()
534  self.cache[(cx, cy)] = im
535  return im
536 
537 def _fitPsfs(fp, peaks, fpres, log, psf, psffwhm, img, varimg,
538  psfChisqCut1, psfChisqCut2, psfChisqCut2b,
539  **kwargs):
540  """!
541  Fit a PSF + smooth background models to a small region around each
542  peak (plus other peaks that are nearby) in the given list of
543  peaks.
544 
545  @param[in] fp Footprint containing the Peaks to model
546  @param[in] peaks a list of all the Peak objects within this Footprint
547  @param[in,out] fpres a PerFootprint result object where results will be placed
548  @param[in,out] log pex Log object
549  @param[in] psf afw PSF object
550  @param[in] psffwhm FWHM of the PSF, in pixels
551  @param[in] img the Image in which this Footprint finds itself
552  @param[in] varimg Variance plane
553  @param[in] psfChisqCut* floats: cuts in chisq-per-pixel at which to accept
554  the PSF model
555 
556  Results go into the fpres.peaks objects.
557  """
558  fbb = fp.getBBox()
559  cpsf = CachingPsf(psf)
560 
561  # create mask image for pixels within the footprint
562  fmask = afwImage.MaskU(fbb)
563  fmask.setXY0(fbb.getMinX(), fbb.getMinY())
564  afwDet.setMaskFromFootprint(fmask, fp, 1)
565 
566  # pk.getF() -- retrieving the floating-point location of the peak
567  # -- actually shows up in the profile if we do it in the loop, so
568  # grab them all here.
569  peakF = [pk.getF() for pk in peaks]
570 
571  for pki,(pk,pkres,pkF) in enumerate(zip(peaks, fpres.peaks, peakF)):
572  log.logdebug('Peak %i' % pki)
573  _fitPsf(fp, fmask, pk, pkF, pkres, fbb, peaks, peakF, log, cpsf, psffwhm,
574  img, varimg, psfChisqCut1, psfChisqCut2, psfChisqCut2b,
575  **kwargs)
576 
577 
578 def _fitPsf(fp, fmask, pk, pkF, pkres, fbb, peaks, peaksF, log, psf, psffwhm,
579  img, varimg, psfChisqCut1, psfChisqCut2, psfChisqCut2b,
580  tinyFootprintSize=2,
581  ):
582  """!
583  Fit a PSF + smooth background model (linear) to a small region around the peak.
584 
585  @param[in] fp Footprint
586  @param[in] fmask the Mask plane for pixels in the Footprint
587  @param[in] pk the Peak within the Footprint that we are going to fit a PSF model for
588  @param[in] pkF the floating-point coordinates of this peak
589  @param[in,out] pkres a PerPeak object that will hold the results for this peak
590  @param[in] fbb the bounding-box of this Footprint (a Box2I)
591  @param[in] peaks a list of all the Peak objects within this Footprint
592  @param[in] peaksF a python list of the floating-point (x,y) peak locations
593  @param[in,out] log pex Log object
594  @param[in] psf afw PSF object
595  @param[in] psffwhm FWHM of the PSF, in pixels
596  @param[in] img the Image in which this Footprint finds itself
597  @param[in] varimg Variance plane
598  @param[in] psfChisqCut* floats: cuts in chisq-per-pixel at which to accept the PSF model
599  """
600  import lsstDebug
601  # my __name__ is lsst.meas.deblender.baseline
602  debugPlots = lsstDebug.Info(__name__).plots
603  debugPsf = lsstDebug.Info(__name__).psf
604 
605  # The small region is a disk out to R0, plus a ramp with
606  # decreasing weight down to R1.
607  R0 = int(math.ceil(psffwhm*1.))
608  # ramp down to zero weight at this radius...
609  R1 = int(math.ceil(psffwhm*1.5))
610  cx, cy = pkF.getX(), pkF.getY()
611  psfimg = psf.computeImage(cx, cy)
612  # R2: distance to neighbouring peak in order to put it into the model
613  R2 = R1 + min(psfimg.getWidth(), psfimg.getHeight())/2.
614 
615  pbb = psfimg.getBBox()
616  pbb.clip(fbb)
617  px0, py0 = psfimg.getX0(), psfimg.getY0()
618 
619  # Make sure we haven't been given a substitute PSF that's nowhere near where we want, as may occur if
620  # "Cannot compute CoaddPsf at point (xx,yy); no input images at that point."
621  if not pbb.contains(afwGeom.Point2I(int(cx), int(cy))):
622  pkres.setOutOfBounds()
623  return
624 
625  # The bounding-box of the local region we are going to fit ("stamp")
626  xlo = int(math.floor(cx - R1))
627  ylo = int(math.floor(cy - R1))
628  xhi = int(math.ceil (cx + R1))
629  yhi = int(math.ceil (cy + R1))
630  stampbb = afwGeom.Box2I(afwGeom.Point2I(xlo, ylo), afwGeom.Point2I(xhi, yhi))
631  stampbb.clip(fbb)
632  xlo, xhi = stampbb.getMinX(), stampbb.getMaxX()
633  ylo, yhi = stampbb.getMinY(), stampbb.getMaxY()
634  if xlo > xhi or ylo > yhi:
635  log.logdebug('Skipping this peak: out of bounds')
636  pkres.setOutOfBounds()
637  return
638 
639  # drop tiny footprints too?
640  if tinyFootprintSize>0 and (((xhi - xlo) < tinyFootprintSize) or
641  ((yhi - ylo) < tinyFootprintSize)):
642  log.logdebug('Skipping this peak: tiny footprint / close to edge')
643  pkres.setTinyFootprint()
644  return
645 
646  # find other peaks within range...
647  otherpeaks = []
648  for pk2, pkF2 in zip(peaks, peaksF):
649  if pk2 == pk:
650  continue
651  if pkF.distanceSquared(pkF2) > R2**2:
652  continue
653  opsfimg = psf.computeImage(pkF2.getX(), pkF2.getY())
654  if not opsfimg.getBBox(afwImage.LOCAL).overlaps(stampbb):
655  continue
656  otherpeaks.append(opsfimg)
657  log.logdebug('%i other peaks within range' % len(otherpeaks))
658 
659  # Now we are going to do a least-squares fit for the flux in this
660  # PSF, plus a decenter term, a linear sky, and fluxes of nearby
661  # sources (assumed point sources). Build up the matrix...
662  # Number of terms -- PSF flux, constant sky, X, Y, + other PSF fluxes
663  NT1 = 4 + len(otherpeaks)
664  # + PSF dx, dy
665  NT2 = NT1 + 2
666  # Number of pixels -- at most
667  NP = (1 + yhi - ylo)*(1 + xhi - xlo)
668  # indices of columns in the "A" matrix.
669  I_psf = 0
670  I_sky = 1
671  I_sky_ramp_x = 2
672  I_sky_ramp_y = 3
673  # offset of other psf fluxes:
674  I_opsf = 4
675  I_dx = NT1 + 0
676  I_dy = NT1 + 1
677 
678  # Build the matrix "A", rhs "b" and weight "w".
679  ix0, iy0 = img.getX0(), img.getY0()
680  fx0, fy0 = fbb.getMinX(), fbb.getMinY()
681  fslice = (slice(ylo-fy0, yhi-fy0+1), slice(xlo-fx0, xhi-fx0+1))
682  islice = (slice(ylo-iy0, yhi-iy0+1), slice(xlo-ix0, xhi-ix0+1))
683  fmask_sub = fmask .getArray()[fslice]
684  var_sub = varimg.getArray()[islice]
685  img_sub = img .getArray()[islice]
686 
687  # Clip the PSF image to match its bbox
688  psfarr = psfimg.getArray()[pbb.getMinY()-py0: 1+pbb.getMaxY()-py0,
689  pbb.getMinX()-px0: 1+pbb.getMaxX()-px0]
690  px0, px1 = pbb.getMinX(), pbb.getMaxX()
691  py0, py1 = pbb.getMinY(), pbb.getMaxY()
692 
693  # Compute the "valid" pixels within our region-of-interest
694  valid = (fmask_sub > 0)
695  xx, yy = np.arange(xlo, xhi+1), np.arange(ylo, yhi+1)
696  RR = ((xx - cx)**2)[np.newaxis, :] + ((yy - cy)**2)[:, np.newaxis]
697  valid *= (RR <= R1**2)
698  valid *= (var_sub > 0)
699  NP = valid.sum()
700 
701  if NP == 0:
702  log.warn('Skipping peak at (%.1f, %.1f): no unmasked pixels nearby' % (cx, cy))
703  pkres.setNoValidPixels()
704  return
705 
706  # pixel coords of valid pixels
707  XX, YY = np.meshgrid(xx, yy)
708  ipixes = np.vstack((XX[valid] - xlo, YY[valid] - ylo)).T
709 
710  inpsfx = (xx >= px0)*(xx <= px1)
711  inpsfy = (yy >= py0)*(yy <= py1)
712  inpsf = np.outer(inpsfy, inpsfx)
713  indx = np.outer(inpsfy, (xx > px0)*(xx < px1))
714  indy = np.outer((yy > py0)*(yy < py1), inpsfx)
715 
716  del inpsfx
717  del inpsfy
718 
719  def _overlap(xlo, xhi, xmin, xmax):
720  assert((xlo <= xmax) and (xhi >= xmin) and
721  (xlo <= xhi) and (xmin <= xmax))
722  xloclamp = max(xlo, xmin)
723  Xlo = xloclamp - xlo
724  xhiclamp = min(xhi, xmax)
725  Xhi = Xlo + (xhiclamp - xloclamp)
726  assert(xloclamp >= 0)
727  assert(Xlo >= 0)
728  return (xloclamp, xhiclamp+1, Xlo, Xhi+1)
729 
730  A = np.zeros((NP, NT2))
731  # Constant term
732  A[:, I_sky] = 1.
733  # Sky slope terms: dx, dy
734  A[:, I_sky_ramp_x] = ipixes[:, 0] + (xlo-cx)
735  A[:, I_sky_ramp_y] = ipixes[:, 1] + (ylo-cy)
736 
737  # whew, grab the valid overlapping PSF pixels
738  px0, px1 = pbb.getMinX(), pbb.getMaxX()
739  py0, py1 = pbb.getMinY(), pbb.getMaxY()
740  sx1, sx2, sx3, sx4 = _overlap(xlo, xhi, px0, px1)
741  sy1, sy2, sy3, sy4 = _overlap(ylo, yhi, py0, py1)
742  dpx0, dpy0 = px0 - xlo, py0 - ylo
743  psf_y_slice = slice(sy3 - dpy0, sy4 - dpy0)
744  psf_x_slice = slice(sx3 - dpx0, sx4 - dpx0)
745  psfsub = psfarr[psf_y_slice, psf_x_slice]
746  vsub = valid[sy1-ylo: sy2-ylo, sx1-xlo: sx2-xlo]
747  A[inpsf[valid], I_psf] = psfsub[vsub]
748 
749  # PSF dx -- by taking the half-difference of shifted-by-one and
750  # shifted-by-minus-one.
751  oldsx = (sx1, sx2, sx3, sx4)
752  sx1, sx2, sx3, sx4 = _overlap(xlo, xhi, px0+1, px1-1)
753  psfsub = (psfarr[psf_y_slice, sx3 - dpx0 + 1: sx4 - dpx0 + 1] -
754  psfarr[psf_y_slice, sx3 - dpx0 - 1: sx4 - dpx0 - 1])/2.
755  vsub = valid[sy1-ylo: sy2-ylo, sx1-xlo: sx2-xlo]
756  A[indx[valid], I_dx] = psfsub[vsub]
757  # revert x indices...
758  (sx1, sx2, sx3, sx4) = oldsx
759 
760  # PSF dy
761  sy1, sy2, sy3, sy4 = _overlap(ylo, yhi, py0+1, py1-1)
762  psfsub = (psfarr[sy3 - dpy0 + 1: sy4 - dpy0 + 1, psf_x_slice] -
763  psfarr[sy3 - dpy0 - 1: sy4 - dpy0 - 1, psf_x_slice])/2.
764  vsub = valid[sy1-ylo: sy2-ylo, sx1-xlo: sx2-xlo]
765  A[indy[valid], I_dy] = psfsub[vsub]
766 
767  # other PSFs...
768  for j, opsf in enumerate(otherpeaks):
769  obb = opsf.getBBox()
770  ino = np.outer((yy >= obb.getMinY())*(yy <= obb.getMaxY()),
771  (xx >= obb.getMinX())*(xx <= obb.getMaxX()))
772  dpx0, dpy0 = obb.getMinX() - xlo, obb.getMinY() - ylo
773  sx1, sx2, sx3, sx4 = _overlap(xlo, xhi, obb.getMinX(), obb.getMaxX())
774  sy1, sy2, sy3, sy4 = _overlap(ylo, yhi, obb.getMinY(), obb.getMaxY())
775  opsfarr = opsf.getArray()
776  psfsub = opsfarr[sy3 - dpy0 : sy4 - dpy0, sx3 - dpx0: sx4 - dpx0]
777  vsub = valid[sy1-ylo: sy2-ylo, sx1-xlo: sx2-xlo]
778  A[ino[valid], I_opsf + j] = psfsub[vsub]
779 
780  b = img_sub[valid]
781 
782  # Weights -- from ramp and image variance map.
783  # Ramp weights -- from 1 at R0 down to 0 at R1.
784  rw = np.ones_like(RR)
785  ii = (RR > R0**2)
786  rr = np.sqrt(RR[ii])
787  rw[ii] = np.maximum(0, 1. - ((rr - R0)/(R1 - R0)))
788  w = np.sqrt(rw[valid]/var_sub[valid])
789  # save the effective number of pixels
790  sumr = np.sum(rw[valid])
791  log.log(-9, 'sumr = %g' % sumr)
792 
793  del ii
794 
795  Aw = A*w[:, np.newaxis]
796  bw = b*w
797 
798  if debugPlots:
799  import pylab as plt
800  plt.clf()
801  N = NT2 + 2
802  R, C = 2, (N+1)/2
803  for i in range(NT2):
804  im1 = np.zeros((1+yhi-ylo, 1+xhi-xlo))
805  im1[ipixes[:, 1], ipixes[:, 0]] = A[:, i]
806  plt.subplot(R, C, i+1)
807  plt.imshow(im1, interpolation='nearest', origin='lower')
808  plt.subplot(R, C, NT2+1)
809  im1 = np.zeros((1+yhi-ylo, 1+xhi-xlo))
810  im1[ipixes[:, 1], ipixes[:, 0]] = b
811  plt.imshow(im1, interpolation='nearest', origin='lower')
812  plt.subplot(R, C, NT2+2)
813  im1 = np.zeros((1+yhi-ylo, 1+xhi-xlo))
814  im1[ipixes[:, 1], ipixes[:, 0]] = w
815  plt.imshow(im1, interpolation='nearest', origin='lower')
816  plt.savefig('A.png')
817 
818  # We do fits with and without the decenter (dx,dy) terms.
819  # Since the dx,dy terms are at the end of the matrix,
820  # we can do that just by trimming off those elements.
821  #
822  # The SVD can fail if there are NaNs in the matrices; this should
823  # really be handled upstream
824  try:
825  # NT1 is number of terms without dx,dy;
826  # X1 is the result without decenter
827  X1, r1, rank1, s1 = np.linalg.lstsq(Aw[:, :NT1], bw)
828  # X2 is with decenter
829  X2, r2, rank2, s2 = np.linalg.lstsq(Aw, bw)
830  except np.linalg.LinAlgError, e:
831  log.log(log.WARN, "Failed to fit PSF to child: %s" % e)
832  pkres.setPsfFitFailed()
833  return
834 
835  log.log(-9, 'r1 r2 %s %s' % (r1, r2))
836 
837  # r is weighted chi-squared = sum over pixels: ramp * (model -
838  # data)**2/sigma**2
839  if len(r1) > 0:
840  chisq1 = r1[0]
841  else:
842  chisq1 = 1e30
843  if len(r2) > 0:
844  chisq2 = r2[0]
845  else:
846  chisq2 = 1e30
847  dof1 = sumr - len(X1)
848  dof2 = sumr - len(X2)
849  log.log(-9, 'dof1, dof2 %g %g' % (dof1, dof2))
850 
851  # This can happen if we're very close to the edge (?)
852  if dof1 <= 0 or dof2 <= 0:
853  log.logdebug('Skipping this peak: bad DOF %g, %g' % (dof1, dof2))
854  pkres.setBadPsfDof()
855  return
856 
857  q1 = chisq1/dof1
858  q2 = chisq2/dof2
859  log.logdebug('PSF fits: chisq/dof = %g, %g' % (q1, q2))
860  ispsf1 = (q1 < psfChisqCut1)
861  ispsf2 = (q2 < psfChisqCut2)
862 
863  pkres.psfFit1 = (chisq1, dof1)
864  pkres.psfFit2 = (chisq2, dof2)
865 
866  # check that the fit PSF spatial derivative terms aren't too big
867  if ispsf2:
868  fdx, fdy = X2[I_dx], X2[I_dy]
869  f0 = X2[I_psf]
870  # as a fraction of the PSF flux
871  dx = fdx/f0
872  dy = fdy/f0
873  ispsf2 = ispsf2 and (abs(dx) < 1. and abs(dy) < 1.)
874  log.logdebug('isPSF2 -- checking derivatives: dx,dy = %g, %g -> %s' %
875  (dx, dy, str(ispsf2)))
876  if not ispsf2:
877  pkres.psfFitBigDecenter = True
878 
879  # Looks like a shifted PSF: try actually shifting the PSF by that amount
880  # and re-evaluate the fit.
881  if ispsf2:
882  psfimg2 = psf.computeImage(cx + dx, cy + dy)
883  # clip
884  pbb2 = psfimg2.getBBox()
885  pbb2.clip(fbb)
886 
887  # Make sure we haven't been given a substitute PSF that's nowhere near where we want, as may occur if
888  # "Cannot compute CoaddPsf at point (xx,yy); no input images at that point."
889  if not pbb2.contains(afwGeom.Point2I(int(cx + dx), int(cy + dy))):
890  ispsf2 = False
891  else:
892  # clip image to bbox
893  px0, py0 = psfimg2.getX0(), psfimg2.getY0()
894  psfarr = psfimg2.getArray()[pbb2.getMinY()-py0:1+pbb2.getMaxY()-py0,
895  pbb2.getMinX()-px0:1+pbb2.getMaxX()-px0]
896  px0, py0 = pbb2.getMinX(), pbb2.getMinY()
897  px1, py1 = pbb2.getMaxX(), pbb2.getMaxY()
898 
899  # yuck! Update the PSF terms in the least-squares fit matrix.
900  Ab = A[:, :NT1]
901 
902  sx1, sx2, sx3, sx4 = _overlap(xlo, xhi, px0, px1)
903  sy1, sy2, sy3, sy4 = _overlap(ylo, yhi, py0, py1)
904  dpx0, dpy0 = px0 - xlo, py0 - ylo
905  psfsub = psfarr[sy3-dpy0:sy4-dpy0, sx3-dpx0:sx4-dpx0]
906  vsub = valid[sy1-ylo:sy2-ylo, sx1-xlo:sx2-xlo]
907  xx, yy = np.arange(xlo, xhi+1), np.arange(ylo, yhi+1)
908  inpsf = np.outer((yy >= py0)*(yy <= py1), (xx >= px0)*(xx <= px1))
909  Ab[inpsf[valid], I_psf] = psfsub[vsub]
910 
911  Aw = Ab*w[:, np.newaxis]
912  # re-solve...
913  Xb, rb, rankb, sb = np.linalg.lstsq(Aw, bw)
914  if len(rb) > 0:
915  chisqb = rb[0]
916  else:
917  chisqb = 1e30
918  dofb = sumr - len(Xb)
919  qb = chisqb/dofb
920  ispsf2 = (qb < psfChisqCut2b)
921  q2 = qb
922  X2 = Xb
923  log.logdebug('shifted PSF: new chisq/dof = %g; good? %s' %
924  (qb, ispsf2))
925  pkres.psfFit3 = (chisqb, dofb)
926 
927  # Which one do we keep?
928  if (((ispsf1 and ispsf2) and (q2 < q1)) or
929  (ispsf2 and not ispsf1)):
930  Xpsf = X2
931  chisq = chisq2
932  dof = dof2
933  log.log(-9, 'dof %g' % dof)
934  log.logdebug('Keeping shifted-PSF model')
935  cx += dx
936  cy += dy
937  pkres.psfFitWithDecenter = True
938  else:
939  # (arbitrarily set to X1 when neither fits well)
940  Xpsf = X1
941  chisq = chisq1
942  dof = dof1
943  log.log(-9, 'dof %g' % dof)
944  log.logdebug('Keeping unshifted PSF model')
945 
946  ispsf = (ispsf1 or ispsf2)
947 
948  # Save the PSF models in images for posterity.
949  if debugPsf:
950  SW, SH = 1+xhi-xlo, 1+yhi-ylo
951  psfmod = afwImage.ImageF(SW, SH)
952  psfmod.setXY0(xlo, ylo)
953  psfderivmodm = afwImage.MaskedImageF(SW, SH)
954  psfderivmod = psfderivmodm.getImage()
955  psfderivmod.setXY0(xlo, ylo)
956  model = afwImage.ImageF(SW, SH)
957  model.setXY0(xlo, ylo)
958  for i in range(len(Xpsf)):
959  for (x, y),v in zip(ipixes, A[:, i]*Xpsf[i]):
960  ix, iy = int(x), int(y)
961  model.set(ix, iy, model.get(ix, iy) + float(v))
962  if i in [I_psf, I_dx, I_dy]:
963  psfderivmod.set(ix, iy, psfderivmod.get(ix, iy) + float(v))
964  for ii in range(NP):
965  x, y = ipixes[ii, :]
966  psfmod.set(int(x), int(y), float(A[ii, I_psf]*Xpsf[I_psf]))
967  modelfp = afwDet.Footprint(fp.getPeaks().getSchema())
968  for (x, y) in ipixes:
969  modelfp.addSpan(int(y+ylo), int(x+xlo), int(x+xlo))
970  modelfp.normalize()
971 
972  pkres.psfFitDebugPsf0Img = psfimg
973  pkres.psfFitDebugPsfImg = psfmod
974  pkres.psfFitDebugPsfDerivImg = psfderivmod
975  pkres.psfFitDebugPsfModel = model
976  pkres.psfFitDebugStamp = img.Factory(img, stampbb, True)
977  pkres.psfFitDebugValidPix = valid # numpy array
978  pkres.psfFitDebugVar = varimg.Factory(varimg, stampbb, True)
979  ww = np.zeros(valid.shape, np.float)
980  ww[valid] = w
981  pkres.psfFitDebugWeight = ww # numpy
982  pkres.psfFitDebugRampWeight = rw
983 
984 
985  # Save things we learned about this peak for posterity...
986  pkres.psfFitR0 = R0
987  pkres.psfFitR1 = R1
988  pkres.psfFitStampExtent = (xlo, xhi, ylo, yhi)
989  pkres.psfFitCenter = (cx, cy)
990  log.log(-9, 'saving chisq,dof %g %g' % (chisq, dof))
991  pkres.psfFitBest = (chisq, dof)
992  pkres.psfFitParams = Xpsf
993  pkres.psfFitFlux = Xpsf[I_psf]
994  pkres.psfFitNOthers = len(otherpeaks)
995 
996  if ispsf:
997  pkres.setDeblendedAsPsf()
998 
999  # replace the template image by the PSF + derivatives
1000  # image.
1001  log.logdebug('Deblending as PSF; setting template to PSF model')
1002 
1003  # Instantiate the PSF model and clip it to the footprint
1004  psfimg = psf.computeImage(cx, cy)
1005  # Scale by fit flux.
1006  psfimg *= Xpsf[I_psf]
1007  psfimg = psfimg.convertF()
1008 
1009  # Clip the Footprint to the PSF model image bbox.
1010  fpcopy = afwDet.Footprint(fp)
1011  psfbb = psfimg.getBBox()
1012  fpcopy.clipTo(psfbb)
1013  bb = fpcopy.getBBox()
1014 
1015  # Copy the part of the PSF model within the clipped footprint.
1016  psfmod = afwImage.ImageF(bb)
1017  afwDet.copyWithinFootprintImage(fpcopy, psfimg, psfmod)
1018  # Save it as our template.
1019  fpcopy.clipToNonzero(psfmod)
1020  fpcopy.normalize()
1021  pkres.setTemplate(psfmod, fpcopy)
1022 
1023  # DEBUG
1024  pkres.setPsfTemplate(psfmod, fpcopy)
1025 
1026 
1027 def _handle_flux_at_edge(log, psffwhm, t1, tfoot, fp, maskedImage,
1028  x0, x1, y0, y1, psf, pk, sigma1, patchEdges):
1029  # Import C++ routines
1030  import lsst.meas.deblender as deb
1031  butils = deb.BaselineUtilsF
1032 
1033  log.logdebug('Found significant flux at template edge.')
1034  # Compute the max of:
1035  # -symmetric-template-clipped image * PSF
1036  # -footprint-clipped image
1037  # Ie, extend the template by the PSF and "fill in" the footprint.
1038  # Then find the symmetric template of that image.
1039 
1040  # The size we'll grow by
1041  S = psffwhm*1.5
1042  # make it an odd integer
1043  S = (int(S + 0.5)/2)*2 + 1
1044 
1045  tbb = tfoot.getBBox()
1046  tbb.grow(S)
1047 
1048  # (footprint+margin)-clipped image;
1049  # we need the pixels OUTSIDE the footprint to be 0.
1050  fpcopy = afwDet.Footprint(fp)
1051  fpcopy = afwDet.growFootprint(fpcopy, S)
1052  fpcopy.clipTo(tbb)
1053  fpcopy.normalize()
1054  padim = maskedImage.Factory(tbb)
1055  afwDet.copyWithinFootprintMaskedImage(fpcopy, maskedImage, padim)
1056 
1057  # find pixels on the edge of the template
1058  edgepix = butils.getSignificantEdgePixels(t1, tfoot, -1e6)
1059 
1060  # instantiate PSF image
1061  xc = int((x0 + x1)/2)
1062  yc = int((y0 + y1)/2)
1063  psfim = psf.computeImage(afwGeom.Point2D(xc, yc))
1064  pbb = psfim.getBBox()
1065  # shift PSF image to by centered on zero
1066  lx, ly = pbb.getMinX(), pbb.getMinY()
1067  psfim.setXY0(lx - xc, ly - yc)
1068  pbb = psfim.getBBox()
1069  # clip PSF to S, if necessary
1070  Sbox = afwGeom.Box2I(afwGeom.Point2I(-S, -S), afwGeom.Extent2I(2*S+1, 2*S+1))
1071  if not Sbox.contains(pbb):
1072  # clip PSF image
1073  psfim = psfim.Factory(psfim, Sbox, afwImage.PARENT, True)
1074  pbb = psfim.getBBox()
1075  px0 = pbb.getMinX()
1076  px1 = pbb.getMaxX()
1077  py0 = pbb.getMinY()
1078  py1 = pbb.getMaxY()
1079 
1080  # Compute the ramped-down edge pixels
1081  ramped = t1.Factory(tbb)
1082  Tout = ramped.getArray()
1083  Tin = t1.getArray()
1084  tx0, ty0 = t1.getX0(), t1.getY0()
1085  ox0, oy0 = ramped.getX0(), ramped.getY0()
1086  P = psfim.getArray()
1087  P /= P.max()
1088  # For each edge pixel, Tout = max(Tout, edgepix * PSF)
1089  for span in edgepix.getSpans():
1090  y = span.getY()
1091  for x in range(span.getX0(), span.getX1()+1):
1092  slc = (slice(y+py0 - oy0, y+py1+1 - oy0),
1093  slice(x+px0 - ox0, x+px1+1 - ox0))
1094  Tout[slc] = np.maximum(Tout[slc], Tin[y-ty0, x-tx0]*P)
1095 
1096  # Fill in the "padim" (which has the right variance and
1097  # mask planes) with the ramped pixels, outside the footprint
1098  I = (padim.getImage().getArray() == 0)
1099  padim.getImage().getArray()[I] = ramped.getArray()[I]
1100 
1101  rtn, patched = butils.buildSymmetricTemplate(padim, fpcopy, pk, sigma1, True, patchEdges)
1102  # silly SWIG can't unpack pairs as tuples
1103  t2, tfoot2 = rtn[0], rtn[1]
1104  del rtn
1105 
1106  # This template footprint may extend outside the parent
1107  # footprint -- or the image. Clip it.
1108  # NOTE that this may make it asymmetric, unlike normal templates.
1109  imbb = maskedImage.getBBox()
1110  tfoot2.clipTo(imbb)
1111  tbb = tfoot2.getBBox()
1112  # clip template image to bbox
1113  t2 = t2.Factory(t2, tbb, afwImage.PARENT, True)
1114 
1115  return t2, tfoot2, patched
def _fitPsf
Fit a PSF + smooth background model (linear) to a small region around the peak.
Definition: baseline.py:581
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.
In the PSF fitting code, we request PSF models for all peaks near the one being fit.
Definition: baseline.py:516
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
def __init__
Creates a result object for the given parent Footprint &#39;fp&#39;.
Definition: baseline.py:40
An integer coordinate rectangle.
Definition: Box.h:53
def getFluxPortion
Return a HeavyFootprint containing the flux apportioned to this peak.
Definition: baseline.py:151
def deblend
Deblend a single footprint in a maskedImage.
Definition: baseline.py:233
A set of pixels in an Image.
Definition: Footprint.h:62
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.
Statistics makeStatistics(afwImage::Mask< afwImage::MaskPixel > const &msk, int const flags, StatisticsControl const &sctrl)
Specialization to handle Masks.
Definition: Statistics.cc:1082
Result of deblending a single Peak within a parent Footprint.
Definition: baseline.py:63
HeavyFootprint< ImagePixelT, MaskPixelT, VariancePixelT > makeHeavyFootprint(Footprint const &foot, lsst::afw::image::MaskedImage< ImagePixelT, MaskPixelT, VariancePixelT > const &img, HeavyFootprintCtrl const *ctrl=NULL)